#include "field.h"

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void conf_read(const char *path, Conf *c)
{
    size_t len;
    char *buf = (char *)read_file(path, 1 << 20, &len);
    memset(c, 0, sizeof *c);
    int lineno = 0;
    char *cur_ = buf;
    for (char *line = next_line(&cur_); line; line = next_line(&cur_)) {
        lineno++;
        char key[16], id[16];
        long long t;
        double a, b;
        if (line[0] == '#' || line[0] == 0)
            continue;
        if (sscanf(line, "%15s", key) != 1)
            die("conf.tsv ligne %d illisible", lineno);
        if (!strcmp(key, "ep0") || !strcmp(key, "ep1") || !strcmp(key, "tr0") || !strcmp(key, "tr1")) {
            if (sscanf(line, "%*s %lld", &t) != 1)
                die("conf.tsv ligne %d : date attendue", lineno);
            int64_t *dst = key[0] == 'e' ? (key[2] == '0' ? &c->ep0 : &c->ep1) : (key[2] == '0' ? &c->tr0 : &c->tr1);
            *dst = t;
        } else if (!strcmp(key, "ref")) {
            if (c->nref >= MAX_REF || sscanf(line, "%*s %lld %lf %lf", &t, &a, &b) != 3)
                die("conf.tsv ligne %d : ref invalide", lineno);
            c->ref_t[c->nref] = t;
            c->ref_lat[c->nref] = a;
            c->ref_lon[c->nref] = b;
            c->nref++;
        } else if (!strcmp(key, "station")) {
            if (c->nst >= MAX_STATIONS || sscanf(line, "%*s %15s %lf %lf", id, &a, &b) != 3)
                die("conf.tsv ligne %d : station invalide", lineno);
            memcpy(c->st_id[c->nst], id, sizeof id);
            c->st_lat[c->nst] = a;
            c->st_lon[c->nst] = b;
            c->nst++;
        } else {
            die("conf.tsv ligne %d : directive inconnue %s", lineno, key);
        }
    }
    free(buf);
    if (!c->ep0 || !c->ep1 || !c->tr0 || !c->tr1 || c->nref == 0)
        die("conf.tsv incomplet");
}

void fill_nearest(const float *in, float *out, int nx, int ny)
{
    size_t n = (size_t)nx * (size_t)ny;
    int *queue = xmalloc(n, sizeof(int));
    size_t head = 0, tail = 0;
    for (size_t k = 0; k < n; k++) {
        out[k] = in[k];
        if (!isnan(in[k]))
            queue[tail++] = (int)k;
    }
    if (tail == 0) {
        free(queue);
        return;
    }
    while (head < tail) {
        int k = queue[head++], x = k % nx, y = k / nx;
        for (int dy = -1; dy <= 1; dy++)
            for (int dx = -1; dx <= 1; dx++) {
                int xx = x + dx, yy = y + dy;
                if ((dx || dy) && xx >= 0 && yy >= 0 && xx < nx && yy < ny) {
                    size_t kk = (size_t)yy * (size_t)nx + (size_t)xx;
                    if (isnan(out[kk])) {
                        out[kk] = out[k];
                        queue[tail++] = (int)kk;
                    }
                }
            }
    }
    free(queue);
}

static int kernel(double sigma, double **w)
{
    int r = (int)(4.0 * sigma + 0.5);
    *w = xmalloc((size_t)(2 * r + 1), sizeof(double));
    double s = 0;
    for (int i = -r; i <= r; i++) {
        (*w)[i + r] = exp(-0.5 * i * i / (sigma * sigma));
        s += (*w)[i + r];
    }
    for (int i = 0; i <= 2 * r; i++)
        (*w)[i] /= s;
    return r;
}

static int clampi(int v, int lo, int hi) { return v < lo ? lo : v > hi ? hi : v; }

void gauss_nearest(const float *in, float *out, int nx, int ny, double sigma)
{
    double *w;
    int r = kernel(sigma, &w);
    size_t n = (size_t)nx * (size_t)ny;
    double *tmp = xmalloc(n, sizeof(double));
    /* noyau séparable : passe en latitude puis en longitude */
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double s = 0;
            for (int k = -r; k <= r; k++)
                s += w[k + r] * in[(size_t)clampi(y + k, 0, ny - 1) * (size_t)nx + (size_t)x];
            tmp[(size_t)y * (size_t)nx + (size_t)x] = s;
        }
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double s = 0;
            for (int k = -r; k <= r; k++)
                s += w[k + r] * tmp[(size_t)y * (size_t)nx + (size_t)clampi(x + k, 0, nx - 1)];
            out[(size_t)y * (size_t)nx + (size_t)x] = (float)s;
        }
    free(tmp);
    free(w);
}

void gauss_masked(const float *in, float *out, int nx, int ny, double sigma)
{
    size_t n = (size_t)nx * (size_t)ny;
    float *v = xmalloc(n, sizeof(float)), *m = xmalloc(n, sizeof(float));
    float *sv = xmalloc(n, sizeof(float)), *sm = xmalloc(n, sizeof(float));
    for (size_t k = 0; k < n; k++) {
        int ok = !isnan(in[k]);
        v[k] = ok ? in[k] : 0.0f;
        m[k] = ok ? 1.0f : 0.0f;
    }
    gauss_nearest(v, sv, nx, ny, sigma);
    gauss_nearest(m, sm, nx, ny, sigma);
    for (size_t k = 0; k < n; k++)
        out[k] = isnan(in[k]) ? NAN : sv[k] / (sm[k] > 1e-6f ? sm[k] : 1e-6f);
    free(v);
    free(m);
    free(sv);
    free(sm);
}

void city_block(const uint8_t *mask, int *dist, int nx, int ny)
{
    const int inf = 1 << 28;
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            size_t k = (size_t)y * (size_t)nx + (size_t)x;
            if (!mask[k]) {
                dist[k] = 0;
                continue;
            }
            int up = y > 0 ? dist[k - (size_t)nx] : 0, left = x > 0 ? dist[k - 1] : 0;
            int d = (up < left ? up : left) + 1;
            dist[k] = d < inf ? d : inf;
        }
    for (int y = ny - 1; y >= 0; y--)
        for (int x = nx - 1; x >= 0; x--) {
            size_t k = (size_t)y * (size_t)nx + (size_t)x;
            if (!dist[k])
                continue;
            int dn = y < ny - 1 ? dist[k + (size_t)nx] : 0, rt = x < nx - 1 ? dist[k + 1] : 0;
            int d = (dn < rt ? dn : rt) + 1;
            if (d < dist[k])
                dist[k] = d;
        }
}

double bilinear(const Grid *g, const float *f, double lat, double lon)
{
    double fy = (lat - g->lat0) / g->dlat, fx = (lon - g->lon0) / g->dlon;
    if (fx < 0 || fy < 0 || fx > g->nx - 1 || fy > g->ny - 1)
        return NAN;
    int x0 = (int)fx, y0 = (int)fy;
    if (x0 >= g->nx - 1)
        x0 = g->nx - 2;
    if (y0 >= g->ny - 1)
        y0 = g->ny - 2;
    if (x0 < 0 || y0 < 0)
        return NAN;
    double ax = fx - x0, ay = fy - y0;
    size_t k = (size_t)y0 * (size_t)g->nx + (size_t)x0;
    double a = f[k], b = f[k + 1], c = f[k + (size_t)g->nx], d = f[k + (size_t)g->nx + 1];
    return (1 - ay) * ((1 - ax) * a + ax * b) + ay * ((1 - ax) * c + ax * d);
}

double dist_km(double lat, double lon, double la0, double lo0)
{
    return KM_PER_DEG * hypot(lat - la0, (lon - lo0) * cos(la0 * DEG));
}