/**
 * earth-now.co
 * Author:  Simon-Pierre Boucher
 * Contact: contact@spboucher.ai
 * File:    packages/models/src/linalg.ts
 * Purpose: Minimal dense linear algebra — least-squares via normal equations with partial-pivot Gaussian elimination
 */

/** Solve A x = b (A: n×n, row-major) with partial pivoting. Throws on singular systems. */
export function solveLinearSystem(A: number[][], b: number[]): number[] {
  const n = b.length;
  // Augmented copy so callers' arrays are never mutated.
  const M = A.map((row, i) => [...row, b[i]!]);

  for (let col = 0; col < n; col++) {
    let pivot = col;
    for (let r = col + 1; r < n; r++) {
      if (Math.abs(M[r]![col]!) > Math.abs(M[pivot]![col]!)) pivot = r;
    }
    if (Math.abs(M[pivot]![col]!) < 1e-12) throw new Error("solveLinearSystem: singular matrix");
    if (pivot !== col) {
      const tmp = M[col]!;
      M[col] = M[pivot]!;
      M[pivot] = tmp;
    }
    for (let r = col + 1; r < n; r++) {
      const f = M[r]![col]! / M[col]![col]!;
      for (let c = col; c <= n; c++) M[r]![c]! -= f * M[col]![c]!;
    }
  }

  const x = new Array<number>(n).fill(0);
  for (let r = n - 1; r >= 0; r--) {
    let acc = M[r]![n]!;
    for (let c = r + 1; c < n; c++) acc -= M[r]![c]! * x[c]!;
    x[r] = acc / M[r]![r]!;
  }
  return x;
}

/**
 * Ordinary least squares: minimize ||X β − y||² via normal equations XᵀX β = Xᵀy.
 * X is m×k row-major (m observations, k basis functions).
 */
export function leastSquares(X: number[][], y: number[]): number[] {
  const m = X.length;
  if (m === 0 || m !== y.length) throw new Error("leastSquares: dimension mismatch");
  const k = X[0]!.length;
  const XtX: number[][] = Array.from({ length: k }, () => new Array<number>(k).fill(0));
  const Xty = new Array<number>(k).fill(0);
  for (let i = 0; i < m; i++) {
    const row = X[i]!;
    for (let a = 0; a < k; a++) {
      Xty[a]! += row[a]! * y[i]!;
      for (let b = a; b < k; b++) XtX[a]![b]! += row[a]! * row[b]!;
    }
  }
  for (let a = 0; a < k; a++) for (let b = 0; b < a; b++) XtX[a]![b] = XtX[b]![a]!;
  return solveLinearSystem(XtX, Xty);
}