// UQO Éval — algèbre linéaire dense (Float64Array, ligne-major) pour l'estimation hédonique.
// Aucune dépendance : Cholesky, résolution, inverse, produits — assez pour k ≤ ~350.

export type Mat = { n: number; k: number; a: Float64Array }; // n lignes × k colonnes, a[i*k+j]

export const zeros = (n: number, k: number): Mat => ({ n, k, a: new Float64Array(n * k) });

/** G = X'X (symétrique k×k) et c = X'y. Coût n·k²/2. */
export function gram(X: Mat, y?: Float64Array, w?: Float64Array): { G: Float64Array; c: Float64Array } {
  const { n, k, a } = X;
  const G = new Float64Array(k * k);
  const c = new Float64Array(k);
  for (let i = 0; i < n; i++) {
    const off = i * k;
    const wi = w ? w[i] : 1;
    for (let j = 0; j < k; j++) {
      const xij = a[off + j] * wi;
      if (xij === 0) continue;
      const rowj = j * k;
      for (let l = j; l < k; l++) G[rowj + l] += xij * a[off + l];
      if (y) c[j] += xij * y[i];
    }
  }
  for (let j = 0; j < k; j++) for (let l = 0; l < j; l++) G[j * k + l] = G[l * k + j];
  return { G, c };
}

/** Cholesky in-place : A = L·L' (A symétrique définie positive, k×k). Retourne L (triangulaire inf.) ou null. */
export function cholesky(A: Float64Array, k: number, jitter = 0): Float64Array | null {
  const L = new Float64Array(k * k);
  for (let j = 0; j < k; j++) {
    let s = A[j * k + j] + jitter;
    for (let p = 0; p < j; p++) s -= L[j * k + p] * L[j * k + p];
    if (!(s > 1e-12)) return null;
    const ljj = Math.sqrt(s);
    L[j * k + j] = ljj;
    for (let i = j + 1; i < k; i++) {
      let t = A[i * k + j];
      for (let p = 0; p < j; p++) t -= L[i * k + p] * L[j * k + p];
      L[i * k + j] = t / ljj;
    }
  }
  return L;
}

/** Résout L L' x = b. */
export function cholSolve(L: Float64Array, k: number, b: Float64Array): Float64Array {
  const y = new Float64Array(k);
  for (let i = 0; i < k; i++) {
    let s = b[i];
    for (let p = 0; p < i; p++) s -= L[i * k + p] * y[p];
    y[i] = s / L[i * k + i];
  }
  const x = new Float64Array(k);
  for (let i = k - 1; i >= 0; i--) {
    let s = y[i];
    for (let p = i + 1; p < k; p++) s -= L[p * k + i] * x[p];
    x[i] = s / L[i * k + i];
  }
  return x;
}

/** Inverse d'une matrice SDP via Cholesky. */
export function cholInverse(L: Float64Array, k: number): Float64Array {
  const inv = new Float64Array(k * k);
  const e = new Float64Array(k);
  for (let j = 0; j < k; j++) {
    e.fill(0);
    e[j] = 1;
    const col = cholSolve(L, k, e);
    for (let i = 0; i < k; i++) inv[i * k + j] = col[i];
  }
  return inv;
}

/**
 * Résout (A + ridge·I) x = b avec repli : si Cholesky échoue (colinéarité parfaite),
 * on ajoute une régularisation infinitésimale croissante. Retourne aussi l'inverse.
 */
export function solveSPD(A: Float64Array, k: number, b: Float64Array, ridge = 0): { x: Float64Array; inv: Float64Array; jitter: number } {
  let jitter = ridge;
  for (let attempt = 0; attempt < 8; attempt++) {
    const L = cholesky(A, k, jitter);
    if (L) return { x: cholSolve(L, k, b), inv: cholInverse(L, k), jitter };
    // échelle du problème : trace moyenne
    let tr = 0;
    for (let j = 0; j < k; j++) tr += A[j * k + j];
    jitter = Math.max(jitter * 10, (tr / k) * 1e-10, 1e-10);
  }
  throw new Error("matrice X'X singulière — retirez une variable ou un effet fixe colinéaire");
}

export function matVec(X: Mat, beta: Float64Array): Float64Array {
  const { n, k, a } = X;
  const out = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    let s = 0;
    const off = i * k;
    for (let j = 0; j < k; j++) s += a[off + j] * beta[j];
    out[i] = s;
  }
  return out;
}

/** Sous-ensemble de lignes. */
export function takeRows(X: Mat, idx: Int32Array | number[]): Mat {
  const k = X.k;
  const out = zeros(idx.length, k);
  for (let r = 0; r < idx.length; r++) out.a.set(X.a.subarray(idx[r] * k, idx[r] * k + k), r * k);
  return out;
}

export function takeVec(v: Float64Array, idx: Int32Array | number[]): Float64Array {
  const out = new Float64Array(idx.length);
  for (let r = 0; r < idx.length; r++) out[r] = v[idx[r]];
  return out;
}

/** Colonnes sélectionnées. */
export function takeCols(X: Mat, cols: number[]): Mat {
  const out = zeros(X.n, cols.length);
  for (let i = 0; i < X.n; i++) for (let j = 0; j < cols.length; j++) out.a[i * cols.length + j] = X.a[i * X.k + cols[j]];
  return out;
}

/** Concatène deux matrices côte à côte. */
export function hcat(A: Mat, B: Mat): Mat {
  const k = A.k + B.k;
  const out = zeros(A.n, k);
  for (let i = 0; i < A.n; i++) {
    out.a.set(A.a.subarray(i * A.k, i * A.k + A.k), i * k);
    out.a.set(B.a.subarray(i * B.k, i * B.k + B.k), i * k + A.k);
  }
  return out;
}

export const mean = (v: ArrayLike<number>): number => {
  let s = 0;
  for (let i = 0; i < v.length; i++) s += v[i];
  return v.length ? s / v.length : NaN;
};

export function quantile(sorted: ArrayLike<number>, q: number): number {
  if (!sorted.length) return NaN;
  const pos = (sorted.length - 1) * q;
  const lo = Math.floor(pos);
  const hi = Math.ceil(pos);
  return sorted[lo] + (sorted[hi] - sorted[lo]) * (pos - lo);
}

export const sortedCopy = (v: ArrayLike<number>): Float64Array => Float64Array.from(v as ArrayLike<number>).sort();

export const median = (v: ArrayLike<number>): number => quantile(sortedCopy(v), 0.5);

export function sd(v: ArrayLike<number>): number {
  const m = mean(v);
  let s = 0;
  for (let i = 0; i < v.length; i++) s += (v[i] - m) ** 2;
  return Math.sqrt(s / Math.max(1, v.length - 1));
}