nis2-agile/application/services/FairService.php

<?php
/**
 * NIS2 Agile - FAIR Service
 *
 * Analisi quantitativa del rischio secondo il modello FAIR
 * (Factor Analysis of Information Risk) con simulazione Monte Carlo.
 *
 *   ALE (Annualized Loss Expectancy) = LEF * Loss Magnitude
 *   LEF (Loss Event Frequency)       = TEF * Vulnerability
 *
 * TEF e Loss Magnitude sono stimati come distribuzioni PERT (min/most_likely/max);
 * la simulazione campiona N iterazioni e restituisce i percentili dell'ALE.
 *
 * Logica PURA (nessun DB): testabile e riusabile.
 */
class FairService
{
    public const DEFAULT_ITERATIONS = 10000;

    /**
     * Esegue la simulazione Monte Carlo FAIR.
     *
     * @param array $p  parametri: tef_min, tef_ml, tef_max, vuln (0-1),
     *                  lm_min, lm_ml, lm_max
     * @param int   $iterations
     * @return array{ale_min:float, ale_ml:float, ale_max:float, ale_mean:float,
     *               lef_mean:float, iterations:int, currency:string, histogram:array}
     */
    public static function simulate(array $p, int $iterations = self::DEFAULT_ITERATIONS): array
    {
        $iterations = max(1000, min(100000, $iterations));

        $tefMin = (float) ($p['tef_min'] ?? 0);
        $tefMl  = (float) ($p['tef_ml']  ?? $tefMin);
        $tefMax = (float) ($p['tef_max'] ?? $tefMl);
        $vuln   = self::clamp((float) ($p['vuln'] ?? 1.0), 0.0, 1.0);
        $lmMin  = (float) ($p['lm_min'] ?? 0);
        $lmMl   = (float) ($p['lm_ml']  ?? $lmMin);
        $lmMax  = (float) ($p['lm_max'] ?? $lmMl);

        // Ordina per robustezza (min <= ml <= max)
        [$tefMin, $tefMl, $tefMax] = self::sort3($tefMin, $tefMl, $tefMax);
        [$lmMin, $lmMl, $lmMax]    = self::sort3($lmMin, $lmMl, $lmMax);

        $samples  = [];
        $sumAle   = 0.0;
        $sumLef   = 0.0;
        // Seed deterministico per riproducibilita del calcolo (no Date/rand di sistema dipendenti)
        mt_srand(crc32(json_encode([$tefMin,$tefMl,$tefMax,$vuln,$lmMin,$lmMl,$lmMax,$iterations])));

        for ($i = 0; $i < $iterations; $i++) {
            $tef = self::pert($tefMin, $tefMl, $tefMax);
            $lef = $tef * $vuln;                 // Loss Event Frequency
            $lm  = self::pert($lmMin, $lmMl, $lmMax);
            $ale = $lef * $lm;
            $samples[] = $ale;
            $sumAle += $ale;
            $sumLef += $lef;
        }

        sort($samples);
        $n = count($samples);

        return [
            'ale_min'    => round(self::percentile($samples, 10), 2),  // P10
            'ale_ml'     => round(self::percentile($samples, 50), 2),  // mediana
            'ale_max'    => round(self::percentile($samples, 90), 2),  // P90
            'ale_mean'   => round($sumAle / $n, 2),
            'lef_mean'   => round($sumLef / $n, 4),
            'iterations' => $iterations,
            'currency'   => 'EUR',
            'percentiles'=> [
                'p10' => round(self::percentile($samples, 10), 2),
                'p50' => round(self::percentile($samples, 50), 2),
                'p90' => round(self::percentile($samples, 90), 2),
                'p99' => round(self::percentile($samples, 99), 2),
            ],
            'histogram'  => self::histogram($samples, 12),
        ];
    }

    /** Campionamento da distribuzione PERT (Beta riscalata) con lambda=4. */
    private static function pert(float $min, float $ml, float $max): float
    {
        if ($max <= $min) return $min;
        $ml = self::clamp($ml, $min, $max);
        $lambda = 4.0;
        $alpha = 1 + $lambda * ($ml - $min) / ($max - $min);
        $beta  = 1 + $lambda * ($max - $ml) / ($max - $min);
        return $min + self::betaSample($alpha, $beta) * ($max - $min);
    }

    /** Campione Beta(alpha,beta) via due Gamma. */
    private static function betaSample(float $a, float $b): float
    {
        $x = self::gammaSample($a);
        $y = self::gammaSample($b);
        $s = $x + $y;
        return $s > 0 ? $x / $s : 0.5;
    }

    /** Campione Gamma(shape, 1) - Marsaglia & Tsang. */
    private static function gammaSample(float $shape): float
    {
        if ($shape < 1) {
            $u = self::rand01();
            return self::gammaSample(1 + $shape) * pow($u, 1 / $shape);
        }
        $d = $shape - 1.0 / 3.0;
        $c = 1.0 / sqrt(9 * $d);
        while (true) {
            do {
                $x = self::randNormal();
                $v = 1 + $c * $x;
            } while ($v <= 0);
            $v = $v * $v * $v;
            $u = self::rand01();
            if ($u < 1 - 0.0331 * $x * $x * $x * $x) return $d * $v;
            if (log($u) < 0.5 * $x * $x + $d * (1 - $v + log($v))) return $d * $v;
        }
    }

    private static function randNormal(): float
    {
        // Box-Muller
        $u1 = self::rand01(); $u2 = self::rand01();
        return sqrt(-2 * log($u1 + 1e-12)) * cos(2 * M_PI * $u2);
    }

    private static function rand01(): float
    {
        return (mt_rand() + 1) / (mt_getrandmax() + 2);
    }

    private static function percentile(array $sorted, float $p): float
    {
        $n = count($sorted);
        if ($n === 0) return 0.0;
        $rank = ($p / 100) * ($n - 1);
        $lo = (int) floor($rank);
        $hi = (int) ceil($rank);
        if ($lo === $hi) return $sorted[$lo];
        $frac = $rank - $lo;
        return $sorted[$lo] * (1 - $frac) + $sorted[$hi] * $frac;
    }

    private static function histogram(array $sorted, int $bins): array
    {
        $n = count($sorted);
        if ($n === 0) return [];
        $min = $sorted[0]; $max = $sorted[$n - 1];
        if ($max <= $min) return [['from' => $min, 'to' => $max, 'count' => $n]];
        $w = ($max - $min) / $bins;
        $h = array_fill(0, $bins, 0);
        foreach ($sorted as $v) {
            $idx = (int) min($bins - 1, floor(($v - $min) / $w));
            $h[$idx]++;
        }
        $out = [];
        for ($i = 0; $i < $bins; $i++) {
            $out[] = ['from' => round($min + $i * $w, 2), 'to' => round($min + ($i + 1) * $w, 2), 'count' => $h[$i]];
        }
        return $out;
    }

    private static function clamp(float $v, float $lo, float $hi): float
    {
        return max($lo, min($hi, $v));
    }

    private static function sort3(float $a, float $b, float $c): array
    {
        $arr = [$a, $b, $c];
        sort($arr);
        return $arr;
    }
}