// UQO Éval — modélisation hédonique : qualité d'ajustement, tests de spécification, VIF.
import { gram, mean, median, solveSPD, sortedCopy, type Mat } from "./linalg";
import { chi2Sf, normalQuantile } from "./stats";

export interface FitMetrics {
  r2: number;
  adjR2: number;
  rmse: number; // en $
  mae: number; // en $
  mdape: number; // fraction
  within10: number;
  within20: number;
  aic: number;
  bic: number;
  sigma: number; // écart-type des résidus (échelle de la variable dépendante)
  smear: number; // facteur de Duan (log) ou 1
  n: number;
  k: number;
}

/** Prédictions en dollars (retransformation de Duan pour le log). */
export function toDollars(fitted: Float64Array, logModel: boolean, smear: number): Float64Array {
  const out = new Float64Array(fitted.length);
  for (let i = 0; i < fitted.length; i++) out[i] = logModel ? Math.exp(fitted[i]) * smear : fitted[i];
  return out;
}

export function smearing(resid: Float64Array, logModel: boolean): number {
  if (!logModel) return 1;
  let s = 0;
  for (let i = 0; i < resid.length; i++) s += Math.exp(resid[i]);
  return s / resid.length;
}

export function fitMetrics(y: Float64Array, fitted: Float64Array, resid: Float64Array, prices: Float64Array, logModel: boolean, k: number, smearOverride?: number): FitMetrics {
  const n = y.length;
  const smear = smearOverride ?? smearing(resid, logModel);
  const pred = toDollars(fitted, logModel, smear);
  const ym = mean(y);
  let sst = 0;
  let sse = 0;
  for (let i = 0; i < n; i++) {
    sst += (y[i] - ym) ** 2;
    sse += resid[i] ** 2;
  }
  const r2 = sst > 0 ? 1 - sse / sst : NaN;
  let se$ = 0;
  let ae$ = 0;
  const ape: number[] = [];
  let w10 = 0;
  let w20 = 0;
  for (let i = 0; i < n; i++) {
    const d = pred[i] - prices[i];
    se$ += d * d;
    ae$ += Math.abs(d);
    const a = Math.abs(d) / prices[i];
    ape.push(a);
    if (a <= 0.1) w10++;
    if (a <= 0.2) w20++;
  }
  return {
    r2,
    adjR2: 1 - ((1 - r2) * (n - 1)) / Math.max(1, n - k),
    rmse: Math.sqrt(se$ / n),
    mae: ae$ / n,
    mdape: median(ape),
    within10: w10 / n,
    within20: w20 / n,
    aic: n * Math.log(sse / n) + 2 * k,
    bic: n * Math.log(sse / n) + k * Math.log(n),
    sigma: Math.sqrt(sse / Math.max(1, n - k)),
    smear,
    n,
    k,
  };
}

/** Breusch-Pagan (version Koenker, robuste à la non-normalité) : n·R² de e² sur X. */
export function breuschPagan(X: Mat, resid: Float64Array): { stat: number; df: number; p: number } {
  const n = X.n;
  const e2 = new Float64Array(n);
  for (let i = 0; i < n; i++) e2[i] = resid[i] * resid[i];
  const { G, c } = gram(X, e2);
  let tr = 0;
  for (let j = 0; j < X.k; j++) tr += G[j * X.k + j];
  const { x: b } = solveSPD(G, X.k, c, (tr / X.k) * 1e-10);
  const m = mean(e2);
  let sst = 0;
  let sse = 0;
  for (let i = 0; i < n; i++) {
    let f = 0;
    for (let j = 0; j < X.k; j++) f += X.a[i * X.k + j] * b[j];
    sst += (e2[i] - m) ** 2;
    sse += (e2[i] - f) ** 2;
  }
  const r2 = sst > 0 ? 1 - sse / sst : 0;
  const stat = n * r2;
  const df = X.k - 1;
  return { stat, df, p: chi2Sf(stat, df) };
}

export function jarqueBera(resid: Float64Array): { stat: number; p: number; skew: number; kurt: number } {
  const n = resid.length;
  const m = mean(resid);
  let m2 = 0, m3 = 0, m4 = 0;
  for (let i = 0; i < n; i++) {
    const d = resid[i] - m;
    m2 += d * d;
    m3 += d * d * d;
    m4 += d * d * d * d;
  }
  m2 /= n; m3 /= n; m4 /= n;
  const skew = m3 / Math.pow(m2, 1.5);
  const kurt = m4 / (m2 * m2);
  const stat = (n / 6) * (skew * skew + ((kurt - 3) ** 2) / 4);
  return { stat, p: chi2Sf(stat, 2), skew, kurt };
}

/** VIF des colonnes demandées : VIF_j = [(X'X)^-1]_jj · Σ(x_ij − x̄_j)². */
export function vif(X: Mat, XtXinv: Float64Array, cols: number[]): number[] {
  const { n, k, a } = X;
  return cols.map((j) => {
    let s = 0;
    for (let i = 0; i < n; i++) s += a[i * k + j];
    const m = s / n;
    let ss = 0;
    for (let i = 0; i < n; i++) ss += (a[i * k + j] - m) ** 2;
    return XtXinv[j * k + j] * ss;
  });
}

/** Points du diagramme quantile-quantile (résidus standardisés). */
export function qqPoints(resid: Float64Array, sigma: number, maxPts = 300): { q: number; e: number }[] {
  const s = sortedCopy(resid);
  const n = s.length;
  const step = Math.max(1, Math.floor(n / maxPts));
  const out: { q: number; e: number }[] = [];
  for (let i = 0; i < n; i += step) {
    const p = (i + 0.5) / n;
    out.push({ q: normalQuantile(p), e: s[i] / sigma });
  }
  return out;
}

export function histogram(v: Float64Array, bins = 30): { lo: number; hi: number; n: number }[] {
  const s = sortedCopy(v);
  const lo = s[Math.floor(s.length * 0.005)];
  const hi = s[Math.min(s.length - 1, Math.ceil(s.length * 0.995))];
  const w = (hi - lo) / bins || 1;
  const out = Array.from({ length: bins }, (_, i) => ({ lo: lo + i * w, hi: lo + (i + 1) * w, n: 0 }));
  for (let i = 0; i < v.length; i++) {
    const b = Math.min(bins - 1, Math.max(0, Math.floor((v[i] - lo) / w)));
    out[b].n++;
  }
  return out;
}