/**
 * earth-now.co
 * Author:  Simon-Pierre Boucher
 * Contact: contact@spboucher.ai
 * File:    packages/counter/src/pchip.ts
 * Purpose: Monotone cubic (PCHIP, Fritsch–Carlson) evaluation — never overshoots, shared by client, server and badge
 */

/**
 * Compute PCHIP endpoint/interior slopes for strictly increasing xs.
 * Fritsch–Carlson: interior slopes are a weighted harmonic mean of adjacent
 * secants, zero at local extrema — this is what guarantees no overshoot
 * (a natural cubic spline could send a counter backwards; PCHIP cannot).
 */
export function pchipSlopes(xs: readonly number[], ys: readonly number[]): number[] {
  const n = xs.length;
  if (n !== ys.length) throw new Error("pchip: xs and ys length mismatch");
  if (n === 0) throw new Error("pchip: need at least one knot");
  if (n === 1) return [0];

  const h: number[] = new Array(n - 1);
  const delta: number[] = new Array(n - 1);
  for (let i = 0; i < n - 1; i++) {
    const dx = xs[i + 1]! - xs[i]!;
    if (dx <= 0) throw new Error("pchip: knots must be strictly increasing in time");
    h[i] = dx;
    delta[i] = (ys[i + 1]! - ys[i]!) / dx;
  }
  if (n === 2) return [delta[0]!, delta[0]!];

  const m: number[] = new Array(n).fill(0);
  for (let i = 1; i < n - 1; i++) {
    const d0 = delta[i - 1]!;
    const d1 = delta[i]!;
    if (d0 * d1 <= 0) {
      m[i] = 0;
    } else {
      const w1 = 2 * h[i]! + h[i - 1]!;
      const w2 = h[i]! + 2 * h[i - 1]!;
      m[i] = (w1 + w2) / (w1 / d0 + w2 / d1);
    }
  }

  m[0] = endpointSlope(h[0]!, h[1]!, delta[0]!, delta[1]!);
  m[n - 1] = endpointSlope(h[n - 2]!, h[n - 3]!, delta[n - 2]!, delta[n - 3]!);
  return m;
}

/** One-sided three-point endpoint slope with the standard monotonicity clamps. */
function endpointSlope(h0: number, h1: number, d0: number, d1: number): number {
  let m = ((2 * h0 + h1) * d0 - h0 * d1) / (h0 + h1);
  if (m * d0 <= 0) m = 0;
  else if (d0 * d1 < 0 && Math.abs(m) > 3 * Math.abs(d0)) m = 3 * d0;
  return m;
}

/**
 * Evaluate the PCHIP interpolant at x. Outside the knot range the value
 * extrapolates linearly with the endpoint derivative (bounded, no polynomial blowup).
 */
export function pchipEvaluate(
  xs: readonly number[],
  ys: readonly number[],
  slopes: readonly number[],
  x: number,
): number {
  const n = xs.length;
  if (n === 1) return ys[0]!;
  if (x <= xs[0]!) return ys[0]! + slopes[0]! * (x - xs[0]!);
  if (x >= xs[n - 1]!) return ys[n - 1]! + slopes[n - 1]! * (x - xs[n - 1]!);

  // Binary search for the segment containing x.
  let lo = 0;
  let hi = n - 1;
  while (hi - lo > 1) {
    const mid = (lo + hi) >> 1;
    if (xs[mid]! <= x) lo = mid;
    else hi = mid;
  }

  const h = xs[lo + 1]! - xs[lo]!;
  const t = (x - xs[lo]!) / h;
  const t2 = t * t;
  const t3 = t2 * t;
  const h00 = 2 * t3 - 3 * t2 + 1;
  const h10 = t3 - 2 * t2 + t;
  const h01 = -2 * t3 + 3 * t2;
  const h11 = t3 - t2;
  return (
    h00 * ys[lo]! + h10 * h * slopes[lo]! + h01 * ys[lo + 1]! + h11 * h * slopes[lo + 1]!
  );
}

/** First derivative of the PCHIP interpolant at x (endpoint slope outside the range). */
export function pchipDerivative(
  xs: readonly number[],
  ys: readonly number[],
  slopes: readonly number[],
  x: number,
): number {
  const n = xs.length;
  if (n === 1) return 0;
  if (x <= xs[0]!) return slopes[0]!;
  if (x >= xs[n - 1]!) return slopes[n - 1]!;

  let lo = 0;
  let hi = n - 1;
  while (hi - lo > 1) {
    const mid = (lo + hi) >> 1;
    if (xs[mid]! <= x) lo = mid;
    else hi = mid;
  }

  const h = xs[lo + 1]! - xs[lo]!;
  const t = (x - xs[lo]!) / h;
  const t2 = t * t;
  const dh00 = (6 * t2 - 6 * t) / h;
  const dh10 = 3 * t2 - 4 * t + 1;
  const dh01 = (-6 * t2 + 6 * t) / h;
  const dh11 = 3 * t2 - 2 * t;
  return dh00 * ys[lo]! + dh10 * slopes[lo]! + dh01 * ys[lo + 1]! + dh11 * slopes[lo + 1]!;
}