/* ensemble.c : « austen ensemble <departements.bin> <sortie> modèle=gmax.grd ... »
 *
 * Interpole chaque carte de rafales maxi sur une grille commune de 0,025° (celle d'AROME) ; un modèle
 * peut recevoir plusieurs runs (« modèle=f1,f2 »), moyennés d'abord entre eux. Puis
 * calcule point par point la moyenne, le minimum, le maximum et le nombre de modèles >= 100 km/h.
 * Écrit <sortie>.grd (couches : moyenne, mini, maxi, n100, puis un modèle par couche) et
 * <sortie>.json (noms des couches ; par département : maximum, médiane et 90e centile de la
 * moyenne, maximum de chaque modèle).
 *
 * departements.bin (écrit par la partie PHP) : "AUPO", u32 nb ; pour chaque département :
 * u8 longueur du code, code, u32 nb d'anneaux extérieurs ; par anneau : u32 nb de points, f64 lon, lat.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "field.h"
#include "json.h"
#include "util.h"

#define LON0 -7.5
#define LON1 10.0
#define LAT0 41.2
#define LAT1 52.4
#define STEP 0.025
#define MAX_MODELS 16

typedef struct {
    char code[8];
    int nring;
    int *npt;
    double **xy;      /* lon, lat alternés */
    double bb[4];     /* lon mini, lat mini, lon maxi, lat maxi */
} Dept;

static uint32_t rd32(const uint8_t *b, size_t len, size_t *p)
{
    if (*p + 4 > len)
        die("departements.bin tronqué");
    uint32_t v;
    memcpy(&v, b + *p, 4);
    *p += 4;
    return v;
}

static Dept *read_depts(const char *path, int *n)
{
    size_t len, p = 4;
    uint8_t *b = read_file(path, 64 << 20, &len);
    if (len < 8 || memcmp(b, "AUPO", 4) != 0)
        die("departements.bin : format inconnu");
    uint32_t nd = rd32(b, len, &p);
    if (nd > 1000)
        die("departements.bin : trop de départements");
    Dept *d = xcalloc(nd, sizeof *d);
    for (uint32_t k = 0; k < nd; k++) {
        if (p >= len || b[p] == 0 || b[p] >= sizeof d[k].code || p + 1 + b[p] > len)
            die("departements.bin : code invalide");
        memcpy(d[k].code, b + p + 1, b[p]);
        p += 1 + (size_t)b[p];
        uint32_t nr = rd32(b, len, &p);
        if (nr > 10000)
            die("departements.bin : trop d'anneaux");
        d[k].nring = (int)nr;
        d[k].npt = xmalloc(nr, sizeof(int));
        d[k].xy = xmalloc(nr, sizeof(double *));
        d[k].bb[0] = d[k].bb[1] = INFINITY;
        d[k].bb[2] = d[k].bb[3] = -INFINITY;
        for (uint32_t r = 0; r < nr; r++) {
            uint32_t np = rd32(b, len, &p);
            if (np > 10000000 || p + (size_t)np * 16 > len)
                die("departements.bin : anneau tronqué");
            d[k].npt[r] = (int)np;
            d[k].xy[r] = xmalloc((size_t)np * 2, sizeof(double));
            memcpy(d[k].xy[r], b + p, (size_t)np * 16);
            p += (size_t)np * 16;
            for (uint32_t q = 0; q < np; q++) {
                double x = d[k].xy[r][2 * q], y = d[k].xy[r][2 * q + 1];
                d[k].bb[0] = fmin(d[k].bb[0], x);
                d[k].bb[1] = fmin(d[k].bb[1], y);
                d[k].bb[2] = fmax(d[k].bb[2], x);
                d[k].bb[3] = fmax(d[k].bb[3], y);
            }
        }
    }
    free(b);
    *n = (int)nd;
    return d;
}

/* point dans un polygone : règle pair-impair (rayon horizontal) */
static int inside_ring(const double *xy, int n, double x, double y)
{
    int in = 0;
    for (int i = 0, j = n - 1; i < n; j = i++) {
        double xi = xy[2 * i], yi = xy[2 * i + 1], xj = xy[2 * j], yj = xy[2 * j + 1];
        if ((yi > y) != (yj > y) && x < (xj - xi) * (y - yi) / (yj - yi) + xi)
            in = !in;
    }
    return in;
}

static int cmp_f(const void *a, const void *b)
{
    float x = *(const float *)a, y = *(const float *)b;
    return x < y ? -1 : x > y;
}

/* centile à interpolation linéaire, comme numpy.percentile (méthode par défaut) */
static double percentile(const float *sorted, size_t n, double q)
{
    double pos = q / 100.0 * (double)(n - 1);
    size_t lo = (size_t)floor(pos), hi = lo + 1 < n ? lo + 1 : lo;
    return sorted[lo] + (pos - (double)lo) * (sorted[hi] - sorted[lo]);
}

int cmd_ensemble(int argc, char **argv)
{
    if (argc < 3)
        die("usage : austen ensemble <departements.bin> <sortie> modèle=gmax.grd ...");
    int nm = argc - 2;
    if (nm > MAX_MODELS)
        die("trop de modèles");
    Grid cg = {(int)lround((LON1 - LON0) / STEP) + 1, (int)lround((LAT1 - LAT0) / STEP) + 1, LAT0, STEP, LON0, STEP};
    size_t np = grid_points(&cg);
    int nl = 4 + nm;
    float *layers = xmalloc((size_t)nl * np, sizeof(float));
    char names[MAX_MODELS][32];
    for (int m = 0; m < nm; m++) {
        const char *arg = argv[2 + m], *eq = strchr(arg, '=');
        if (!eq || eq == arg || (size_t)(eq - arg) >= sizeof names[m])
            die("argument modèle=fichier attendu : %s", arg);
        memcpy(names[m], arg, (size_t)(eq - arg));
        names[m][eq - arg] = 0;
        /* « modèle=f1,f2,... » : plusieurs runs du même modèle, moyennés point par point avant la
         * moyenne entre modèles, pour qu'un modèle à 8 runs par jour ne pèse pas plus qu'un autre */
        float *dst = layers + (size_t)(4 + m) * np;
        float *sum = xcalloc(np, sizeof(float));
        int *cnt = xcalloc(np, sizeof(int));
        char *list = xstrdup(eq + 1), *cur = list;
        int nfiles = 0;
        while (cur && *cur) {
            char *comma = strchr(cur, ',');
            if (comma)
                *comma = 0;
            Grd *d = grd_read(cur);
            if (d->nt_gust != 1)
                die("%s : une seule couche attendue", cur);
            for (int y = 0; y < cg.ny; y++)
                for (int x = 0; x < cg.nx; x++) {
                    double v = bilinear(&d->g, d->gust, cg.lat0 + y * cg.dlat, cg.lon0 + x * cg.dlon);
                    size_t q = (size_t)y * (size_t)cg.nx + (size_t)x;
                    if (!isnan(v)) {
                        sum[q] += (float)v;
                        cnt[q]++;
                    }
                }
            grd_free(d);
            nfiles++;
            cur = comma ? comma + 1 : NULL;
        }
        if (!nfiles)
            die("aucun fichier pour %s", names[m]);
        for (size_t q = 0; q < np; q++)
            dst[q] = cnt[q] ? sum[q] / (float)cnt[q] : NAN;
        free(list);
        free(sum);
        free(cnt);
    }
    float *mean = layers, *mn = layers + np, *mx = layers + 2 * np, *n100 = layers + 3 * np;
    for (size_t k = 0; k < np; k++) {
        double s = 0;
        int c = 0, over = 0;
        float lo = NAN, hi = NAN;
        for (int m = 0; m < nm; m++) {
            float v = layers[(size_t)(4 + m) * np + k];
            if (isnan(v))
                continue;
            s += v;
            c++;
            over += v >= 100.0f;
            lo = isnan(lo) || v < lo ? v : lo;
            hi = isnan(hi) || v > hi ? v : hi;
        }
        mean[k] = c ? (float)(s / c) : NAN;
        mn[k] = lo;
        mx[k] = hi;
        n100[k] = (float)over;
    }
    char path[4096];
    if (snprintf(path, sizeof path, "%s.grd", argv[1]) >= (int)sizeof path)
        die("chemin trop long");
    int64_t *lt = xmalloc((size_t)nl, sizeof(int64_t));
    for (int l = 0; l < nl; l++)
        lt[l] = l;
    Grd out = {cg, 0, nl, NULL, lt, NULL, layers};
    grd_write(path, &out);

    /* départements */
    int nd;
    Dept *dep = read_depts(argv[0], &nd);
    Json j;
    js_init(&j);
    js_obj(&j);
    js_key(&j, "layers");
    js_arr(&j);
    js_str(&j, "mean");
    js_str(&j, "min");
    js_str(&j, "max");
    js_str(&j, "n100");
    for (int m = 0; m < nm; m++)
        js_str(&j, names[m]);
    js_end_arr(&j);
    js_key(&j, "mean_max");
    double mm = -INFINITY;
    for (size_t k = 0; k < np; k++)
        if (!isnan(mean[k]) && mean[k] > mm)
            mm = mean[k];
    js_num(&j, mm, 1);
    js_key(&j, "depts");
    js_arr(&j);
    float *vals = xmalloc(np, sizeof(float));
    for (int k = 0; k < nd; k++) {
        size_t nv = 0;
        double per[MAX_MODELS];
        for (int m = 0; m < nm; m++)
            per[m] = NAN;
        int y0 = (int)floor((dep[k].bb[1] - cg.lat0) / cg.dlat), y1 = (int)ceil((dep[k].bb[3] - cg.lat0) / cg.dlat);
        int x0 = (int)floor((dep[k].bb[0] - cg.lon0) / cg.dlon), x1 = (int)ceil((dep[k].bb[2] - cg.lon0) / cg.dlon);
        for (int y = y0 < 0 ? 0 : y0; y <= y1 && y < cg.ny; y++)
            for (int x = x0 < 0 ? 0 : x0; x <= x1 && x < cg.nx; x++) {
                double lo = cg.lon0 + x * cg.dlon, la = cg.lat0 + y * cg.dlat;
                int in = 0;
                for (int r = 0; r < dep[k].nring && !in; r++)
                    in = inside_ring(dep[k].xy[r], dep[k].npt[r], lo, la);
                if (!in)
                    continue;
                size_t q = (size_t)y * (size_t)cg.nx + (size_t)x;
                if (!isnan(mean[q]))
                    vals[nv++] = mean[q];
                for (int m = 0; m < nm; m++) {
                    float v = layers[(size_t)(4 + m) * np + q];
                    if (!isnan(v) && (isnan(per[m]) || v > per[m]))
                        per[m] = v;
                }
            }
        if (nv < 5)
            continue;
        qsort(vals, nv, sizeof(float), cmp_f);
        js_obj(&j);
        js_key(&j, "code");
        js_str(&j, dep[k].code);
        js_key(&j, "max");
        js_num(&j, round(vals[nv - 1]), 0);
        js_key(&j, "med");
        js_num(&j, round(percentile(vals, nv, 50)), 0);
        js_key(&j, "p90");
        js_num(&j, round(percentile(vals, nv, 90)), 0);
        js_key(&j, "models");
        js_obj(&j);
        for (int m = 0; m < nm; m++)
            if (!isnan(per[m])) {
                js_key(&j, names[m]);
                js_num(&j, round(per[m]), 0);
            }
        js_end_obj(&j);
        js_end_obj(&j);
    }
    js_end_arr(&j);
    js_end_obj(&j);
    if (snprintf(path, sizeof path, "%s.json", argv[1]) >= (int)sizeof path)
        die("chemin trop long");
    js_save(&j, path);
    js_free(&j);
    free(vals);
    for (int k = 0; k < nd; k++) {
        for (int r = 0; r < dep[k].nring; r++)
            free(dep[k].xy[r]);
        free(dep[k].xy);
        free(dep[k].npt);
    }
    free(dep);
    free(lt);
    free(layers);
    return 0;
}