// UQO Éval — modélisation hédonique : voisinage spatial (k plus proches voisins), matrice W, Moran.
import { normalCdf } from "./stats";

/** Index spatial par grille pour des requêtes kNN rapides (coordonnées → km locaux). */
export class KnnIndex {
  private xs: Float64Array;
  private ys: Float64Array;
  private cell: number;
  private buckets = new Map<string, number[]>();
  readonly n: number;

  constructor(lat: ArrayLike<number>, lng: ArrayLike<number>) {
    this.n = lat.length;
    let mlat = 0;
    for (let i = 0; i < this.n; i++) mlat += lat[i];
    mlat /= Math.max(1, this.n);
    const cos = Math.cos((mlat * Math.PI) / 180);
    this.cosLat = cos;
    this.xs = new Float64Array(this.n);
    this.ys = new Float64Array(this.n);
    let minx = Infinity, maxx = -Infinity, miny = Infinity, maxy = -Infinity;
    for (let i = 0; i < this.n; i++) {
      this.xs[i] = lng[i] * 111 * cos;
      this.ys[i] = lat[i] * 111;
      minx = Math.min(minx, this.xs[i]); maxx = Math.max(maxx, this.xs[i]);
      miny = Math.min(miny, this.ys[i]); maxy = Math.max(maxy, this.ys[i]);
    }
    // ~25 points par cellule
    const area = Math.max(1e-6, (maxx - minx) * (maxy - miny));
    this.cell = Math.max(0.05, Math.sqrt((area / Math.max(1, this.n)) * 25));
    for (let i = 0; i < this.n; i++) {
      const key = this.key(this.xs[i], this.ys[i]);
      const b = this.buckets.get(key);
      if (b) b.push(i);
      else this.buckets.set(key, [i]);
    }
  }
  private key(x: number, y: number) {
    return `${Math.floor(x / this.cell)}|${Math.floor(y / this.cell)}`;
  }
  private cosLat = 1;
  /** k voisins les plus proches d'un point (km locaux), en excluant éventuellement un index. */
  queryXY(x: number, y: number, k: number, exclude = -1): { idx: Int32Array; dist: Float64Array } {
    const cx = Math.floor(x / this.cell);
    const cy = Math.floor(y / this.cell);
    const best: { i: number; d: number }[] = [];
    let ring = 0;
    const maxRing = 400;
    while (ring <= maxRing) {
      // les cellules de l'anneau `ring`
      for (let dx = -ring; dx <= ring; dx++) {
        for (let dy = -ring; dy <= ring; dy++) {
          if (Math.max(Math.abs(dx), Math.abs(dy)) !== ring) continue;
          const b = this.buckets.get(`${cx + dx}|${cy + dy}`);
          if (!b) continue;
          for (const i of b) {
            if (i === exclude) continue;
            const d = Math.hypot(this.xs[i] - x, this.ys[i] - y);
            if (best.length < k) {
              best.push({ i, d });
              if (best.length === k) best.sort((a, b2) => a.d - b2.d);
            } else if (d < best[k - 1].d) {
              best[k - 1] = { i, d };
              best.sort((a, b2) => a.d - b2.d);
            }
          }
        }
      }
      // arrêt : k trouvés et l'anneau suivant ne peut pas contenir plus proche
      if (best.length >= k && ring * this.cell >= best[k - 1].d) break;
      if (best.length >= k && ring > 60) break;
      ring++;
    }
    best.sort((a, b2) => a.d - b2.d);
    const m = Math.min(k, best.length);
    const idx = new Int32Array(m);
    const dist = new Float64Array(m);
    for (let j = 0; j < m; j++) {
      idx[j] = best[j].i;
      dist[j] = best[j].d;
    }
    return { idx, dist };
  }
  queryLatLng(lat: number, lng: number, k: number, exclude = -1) {
    return this.queryXY(lng * 111 * this.cosLat, lat * 111, k, exclude);
  }
  /** Voisinage de chaque point de l'index (sans lui-même). */
  selfNeighbours(k: number): Int32Array[] {
    const out: Int32Array[] = new Array(this.n);
    for (let i = 0; i < this.n; i++) out[i] = this.queryXY(this.xs[i], this.ys[i], k, i).idx;
    return out;
  }
}

/** W·v avec W = kNN ligne-standardisée (moyenne des voisins). */
export function lag(nb: Int32Array[], v: ArrayLike<number>): Float64Array {
  const out = new Float64Array(nb.length);
  for (let i = 0; i < nb.length; i++) {
    const ids = nb[i];
    if (!ids.length) continue;
    let s = 0;
    for (let j = 0; j < ids.length; j++) s += v[ids[j]];
    out[i] = s / ids.length;
  }
  return out;
}

export interface Moran {
  I: number;
  expected: number;
  z: number;
  p: number;
  k: number;
}

/** I de Moran des résidus sous W kNN ligne-standardisée (test de normalité asymptotique). */
export function moranI(nb: Int32Array[], e: Float64Array): Moran {
  const n = e.length;
  const k = nb[0]?.length ?? 0;
  if (n < 10 || !k) return { I: NaN, expected: NaN, z: NaN, p: NaN, k };
  let m = 0;
  for (let i = 0; i < n; i++) m += e[i];
  m /= n;
  const We = lag(nb, e);
  let num = 0;
  let den = 0;
  for (let i = 0; i < n; i++) {
    num += (e[i] - m) * (We[i] - m); // Σ_i (e_i−m)·Σ_j w_ij (e_j−m) (lag conserve la moyenne)
    den += (e[i] - m) ** 2;
  }
  const S0 = n; // lignes standardisées
  const I = (n / S0) * (num / den);
  // S1, S2
  const sets = nb.map((ids) => new Set(Array.from(ids)));
  const colSum = new Float64Array(n); // c_i = Σ_j w_ji
  for (let i = 0; i < n; i++) for (const j of nb[i]) colSum[j] += 1 / k;
  let S1 = 0;
  for (let i = 0; i < n; i++)
    for (const j of nb[i]) {
      const wji = sets[j].has(i) ? 1 / k : 0;
      S1 += (1 / k + wji) ** 2;
    }
  // chaque paire symétrique a été visitée deux fois (depuis i et depuis j) : S1 = Σ_visites − ½·Σ_sym
  let sym = 0;
  for (let i = 0; i < n; i++) for (const j of nb[i]) if (sets[j].has(i)) sym += (2 / k) ** 2;
  S1 = S1 - sym / 2;
  let S2 = 0;
  for (let i = 0; i < n; i++) S2 += (1 + colSum[i]) ** 2;
  const EI = -1 / (n - 1);
  const varI = (n * n * S1 - n * S2 + 3 * S0 * S0) / ((n * n - 1) * S0 * S0) - EI * EI;
  const z = varI > 0 ? (I - EI) / Math.sqrt(varI) : NaN;
  const p = Number.isFinite(z) ? 2 * (1 - normalCdf(Math.abs(z))) : NaN;
  return { I, expected: EI, z, p, k };
}