<?php
/**
 * Classement des modèles : verify.php — chaque run est comparé à ce qui a été observé, puis les runs
 * d'un même modèle sont moyennés ; le vainqueur est le modèle dont l'écart moyen est le plus faible.
 *
 * Rafales (38 stations suivies) : pour chaque heure H de l'épisode déjà observée, rafale maxi prévue sur
 * l'heure (H-1, H] contre la plus forte rafale relevée sur la même heure (paquets 6 minutes, 8 sur 10
 * au moins). Deux écarts par run : sur la rafale maxi de chaque station (sur les heures que le run
 * couvre), et heure par heure (qui juge aussi le moment). Un run n'est jugé que sur ce qui suit son départ.
 * Centre de la dépression : position et pression prévues contre la trajectoire analysée (analyses des
 * modèles, à partir de 3 h après le départ du run).
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';

function dist_km(float $la, float $lo, float $la0, float $lo0): float
{
    return 111.19 * hypot($la - $la0, ($lo - $lo0) * cos(deg2rad($la0)));
}

/** Rafale observée maxi par station et par heure (fin d'heure, epoch) : [sid][H] = km/h. */
function observed_hourly(array $cache): array
{
    $acc = [];
    foreach ($cache as $k => [, $rows]) {
        $k = (string) $k;
        $t = gmmktime((int) substr($k, 8, 2), (int) substr($k, 10, 2), 0, (int) substr($k, 4, 2), (int) substr($k, 6, 2), (int) substr($k, 0, 4));
        $h = (int) (ceil($t / 3600) * 3600);          /* 6 min jusqu'à H inclus → heure H */
        foreach ($rows as $sid => $r) {
            if ($r[2] === null) {
                continue;
            }
            $sid = sprintf('%08d', (int) $sid);
            $a = $acc[$sid][$h] ?? [0, -1.0];
            $acc[$sid][$h] = [$a[0] + 1, max($a[1], (float) $r[2])];
        }
    }
    $out = [];
    foreach ($acc as $sid => $hours) {
        foreach ($hours as $h => [$n, $mx]) {
            if ($n >= 8) {
                $out[$sid][$h] = $mx;
            }
        }
    }
    return $out;
}

function mean(array $v): ?float { return $v ? array_sum($v) / count($v) : null; }

/**
 * Écarts d'une prévision de rafales sur les heures (max(h0, départ + 1 h) .. h1] : $fc(id, heure) donne la
 * rafale prévue ou null. Retourne null si moins de 10 stations sont comparables.
 */
function gust_score(callable $fc, int $rt, array $obs, array $ids, int $h0, int $h1): ?array
{
    $errH = $errMax = [];
    foreach ($ids as $sid) {
        if (!isset($obs[$sid])) {
            continue;
        }
        $fmax = $omax = null;
        for ($h = max($h0, $rt + 3600); $h <= $h1; $h += 3600) {
            $o = $obs[$sid][$h] ?? null;
            if ($o !== null) {
                $omax = max($omax ?? $o, $o);
            }
            $f = $fc($sid, $h);
            if ($f === null) {
                continue;
            }
            $fmax = max($fmax ?? $f, $f);
            if ($o !== null) {
                $errH[] = $f - $o;
            }
        }
        if ($fmax !== null && $omax !== null) {
            $errMax[] = $fmax - $omax;
        }
    }
    if (count($errMax) < 10) {
        return null;
    }
    return ['n_st' => count($errMax), 'n_h' => count($errH),
            'mae_max' => round(mean(array_map('abs', $errMax)), 1), 'bias_max' => round(mean($errMax), 1),
            'mae_h' => $errH ? round(mean(array_map('abs', $errH)), 1) : null, 'bias_h' => $errH ? round(mean($errH), 1) : null];
}

/** Série horaire d'un run : [id][heure] = rafale prévue. */
function run_series(string $m, string $r): array
{
    $d = json_in(data_dir() . "/derived/$m/$r.json");
    $out = [];
    foreach ($d['stations'] ?? [] as $sid => $v) {
        foreach ($d['gust_times'] as $i => $t) {
            if (($v['series'][$i] ?? null) !== null) {
                $out[sprintf('%08d', (int) $sid)][utc($t)] = $v['series'][$i];
            }
        }
    }
    return $out;
}

/** Moyenne de plusieurs séries de runs : d'abord par modèle, puis entre modèles. */
function blend_series(array $byModel): array
{
    $perModel = [];
    foreach ($byModel as $m => $runs) {
        $acc = [];
        foreach ($runs as $r) {
            foreach (run_series($m, $r) as $sid => $hours) {
                foreach ($hours as $h => $v) {
                    $acc[$sid][$h][] = $v;
                }
            }
        }
        foreach ($acc as $sid => $hours) {
            foreach ($hours as $h => $vs) {
                $perModel[$sid][$h][] = array_sum($vs) / count($vs);
            }
        }
    }
    return array_map(fn($hours) => array_map(fn($vs) => array_sum($vs) / count($vs), $hours), $perModel);
}

function main(): int
{
    $pub = pub_dir();
    $cachePath = data_dir() . '/obs_cache.json';
    if (!is_file($cachePath) || !is_file("$pub/data.json")) {
        return 0;
    }
    $obs = observed_hourly(json_in($cachePath));
    $data = json_in("$pub/data.json");
    $analysis = array_values(array_filter($data['analysis'] ?? [], fn($a) => empty($a['f'])));
    $anaAt = array_column($analysis, null, 't');
    $lastObs = 0;
    foreach ($obs as $hours) {
        $lastObs = max($lastObs, max(array_keys($hours)));
    }
    $h0 = ep0() + 3600;
    $h1 = min(ep1(), $lastObs);
    $ids = array_column(STATIONS, 0);
    $models = [];
    foreach (array_keys(MODELS) as $m) {
        $runs = [];
        foreach (glob(data_dir() . "/derived/$m/*.json") ?: [] as $p) {
            $r = basename($p, '.json');
            if (!valid_run($r)) {
                continue;
            }
            $rt = run_time($r);
            $d = json_in($p);
            $row = ['run' => $r];
            /* ---- rafales */
            if (MODELS[$m]['gust'] && isset($d['stations']) && $h1 >= $h0) {
                $ser = run_series($m, $r);
                $sc = gust_score(fn($sid, $h) => $ser[$sid][$h] ?? null, $rt, $obs, $ids, $h0, $h1);
                if ($sc) {
                    $row += $sc;
                }
            }
            /* ---- centre de la dépression */
            $dk = $dp = [];
            foreach ($d['track'] as $pt) {
                $t = utc($pt['t']);
                $a = $anaAt[$pt['t']] ?? null;
                if ($a && $t >= $rt + 3 * 3600) {
                    $dk[] = dist_km($pt['lat'], $pt['lon'], $a['lat'], $a['lon']);
                    $dp[] = $pt['p'] - $a['p'];
                }
            }
            if ($dk) {
                $row += ['n_c' => count($dk), 'km' => round(mean($dk)), 'dp' => round(mean($dp), 1), 'adp' => round(mean(array_map('abs', $dp)), 1)];
            }
            if (count($row) > 1) {
                $runs[] = $row;
            }
        }
        usort($runs, fn($a, $b) => strcmp($a['run'], $b['run']));
        $g = array_values(array_filter($runs, fn($x) => isset($x['mae_max'])));
        $c = array_values(array_filter($runs, fn($x) => isset($x['km'])));
        $models[$m] = ['runs' => $runs, 'n_runs' => count($g),
            'mae_max' => $g ? round(mean(array_column($g, 'mae_max')), 1) : null,
            'bias_max' => $g ? round(mean(array_column($g, 'bias_max')), 1) : null,
            'mae_h' => $g ? round(mean(array_filter(array_column($g, 'mae_h'), fn($v) => $v !== null)) ?? 0, 1) : null,
            'n_track' => count($c), 'km' => $c ? round(mean(array_column($c, 'km'))) : null,
            'adp' => $c ? round(mean(array_column($c, 'adp')), 1) : null, 'dp' => $c ? round(mean(array_column($c, 'dp')), 1) : null];
    }
    /* deux concurrents de plus : la moyenne des modèles (derniers runs partis avant l'épisode, comme la carte)
       et la moyenne de tous les runs partis avant l'épisode, chaque modèle moyenné sur ses runs d'abord */
    if ($h1 >= $h0) {
        $fc = array_filter($data['selection']['forecast'] ?? [], fn($m) => MODELS[$m]['gust'] ?? false, ARRAY_FILTER_USE_KEY);
        foreach (['mean' => array_map(fn($r) => [$r], $fc), 'allruns' => $data['allruns'] ?? []] as $key => $byModel) {
            if (!$byModel) {
                continue;
            }
            $ser = blend_series($byModel);
            $sc = gust_score(fn($sid, $h) => $ser[$sid][$h] ?? null, ep0(), $obs, $ids, $h0, $h1);
            if ($sc) {
                $models[$key] = $sc + ['runs' => [], 'n_runs' => array_sum(array_map('count', $byModel)), 'n_track' => 0,
                                       'km' => null, 'adp' => null, 'dp' => null, 'pseudo' => true];
            }
        }
    }
    $ranked = array_keys(array_filter($models, fn($v) => $v['mae_max'] !== null));
    usort($ranked, fn($a, $b) => $models[$a]['mae_max'] <=> $models[$b]['mae_max']);
    $rankedC = array_keys(array_filter($models, fn($v) => $v['km'] !== null));
    usort($rankedC, fn($a, $b) => ($models[$a]['km'] + 20 * $models[$a]['adp']) <=> ($models[$b]['km'] + 20 * $models[$b]['adp']));
    json_out("$pub/verif.json", ['updated' => gmdate('Y-m-d\TH:i:s\Z'),
        'window' => $h1 >= $h0 ? [iso($h0 - 3600), iso($h1)] : null, 'hours' => $h1 >= $h0 ? intdiv($h1 - $h0, 3600) + 1 : 0,
        'analyses' => count($analysis), 'models' => $models, 'ranking' => $ranked, 'ranking_center' => $rankedC,
        'winner' => $ranked[0] ?? null]);
    return 0;
}

exit(main());