spb/earth-now Public License

1/**2 * earth-now.co3 * Author:  Simon-Pierre Boucher4 * Contact: contact@spboucher.ai5 * File:    packages/counter/src/pchip.ts6 * Purpose: Monotone cubic (PCHIP, Fritsch–Carlson) evaluation — never overshoots, shared by client, server and badge7 */89/**10 * Compute PCHIP endpoint/interior slopes for strictly increasing xs.11 * Fritsch–Carlson: interior slopes are a weighted harmonic mean of adjacent12 * secants, zero at local extrema — this is what guarantees no overshoot13 * (a natural cubic spline could send a counter backwards; PCHIP cannot).14 */15export function pchipSlopes(xs: readonly number[], ys: readonly number[]): number[] {16  const n = xs.length;17  if (n !== ys.length) throw new Error("pchip: xs and ys length mismatch");18  if (n === 0) throw new Error("pchip: need at least one knot");19  if (n === 1) return [0];2021  const h: number[] = new Array(n - 1);22  const delta: number[] = new Array(n - 1);23  for (let i = 0; i < n - 1; i++) {24    const dx = xs[i + 1]! - xs[i]!;25    if (dx <= 0) throw new Error("pchip: knots must be strictly increasing in time");26    h[i] = dx;27    delta[i] = (ys[i + 1]! - ys[i]!) / dx;28  }29  if (n === 2) return [delta[0]!, delta[0]!];3031  const m: number[] = new Array(n).fill(0);32  for (let i = 1; i < n - 1; i++) {33    const d0 = delta[i - 1]!;34    const d1 = delta[i]!;35    if (d0 * d1 <= 0) {36      m[i] = 0;37    } else {38      const w1 = 2 * h[i]! + h[i - 1]!;39      const w2 = h[i]! + 2 * h[i - 1]!;40      m[i] = (w1 + w2) / (w1 / d0 + w2 / d1);41    }42  }4344  m[0] = endpointSlope(h[0]!, h[1]!, delta[0]!, delta[1]!);45  m[n - 1] = endpointSlope(h[n - 2]!, h[n - 3]!, delta[n - 2]!, delta[n - 3]!);46  return m;47}4849/** One-sided three-point endpoint slope with the standard monotonicity clamps. */50function endpointSlope(h0: number, h1: number, d0: number, d1: number): number {51  let m = ((2 * h0 + h1) * d0 - h0 * d1) / (h0 + h1);52  if (m * d0 <= 0) m = 0;53  else if (d0 * d1 < 0 && Math.abs(m) > 3 * Math.abs(d0)) m = 3 * d0;54  return m;55}5657/**58 * Evaluate the PCHIP interpolant at x. Outside the knot range the value59 * extrapolates linearly with the endpoint derivative (bounded, no polynomial blowup).60 */61export function pchipEvaluate(62  xs: readonly number[],63  ys: readonly number[],64  slopes: readonly number[],65  x: number,66): number {67  const n = xs.length;68  if (n === 1) return ys[0]!;69  if (x <= xs[0]!) return ys[0]! + slopes[0]! * (x - xs[0]!);70  if (x >= xs[n - 1]!) return ys[n - 1]! + slopes[n - 1]! * (x - xs[n - 1]!);7172  // Binary search for the segment containing x.73  let lo = 0;74  let hi = n - 1;75  while (hi - lo > 1) {76    const mid = (lo + hi) >> 1;77    if (xs[mid]! <= x) lo = mid;78    else hi = mid;79  }8081  const h = xs[lo + 1]! - xs[lo]!;82  const t = (x - xs[lo]!) / h;83  const t2 = t * t;84  const t3 = t2 * t;85  const h00 = 2 * t3 - 3 * t2 + 1;86  const h10 = t3 - 2 * t2 + t;87  const h01 = -2 * t3 + 3 * t2;88  const h11 = t3 - t2;89  return (90    h00 * ys[lo]! + h10 * h * slopes[lo]! + h01 * ys[lo + 1]! + h11 * h * slopes[lo + 1]!91  );92}9394/** First derivative of the PCHIP interpolant at x (endpoint slope outside the range). */95export function pchipDerivative(96  xs: readonly number[],97  ys: readonly number[],98  slopes: readonly number[],99  x: number,100): number {101  const n = xs.length;102  if (n === 1) return 0;103  if (x <= xs[0]!) return slopes[0]!;104  if (x >= xs[n - 1]!) return slopes[n - 1]!;105106  let lo = 0;107  let hi = n - 1;108  while (hi - lo > 1) {109    const mid = (lo + hi) >> 1;110    if (xs[mid]! <= x) lo = mid;111    else hi = mid;112  }113114  const h = xs[lo + 1]! - xs[lo]!;115  const t = (x - xs[lo]!) / h;116  const t2 = t * t;117  const dh00 = (6 * t2 - 6 * t) / h;118  const dh10 = 3 * t2 - 4 * t + 1;119  const dh01 = (-6 * t2 + 6 * t) / h;120  const dh11 = 3 * t2 - 2 * t;121  return dh00 * ys[lo]! + dh10 * slopes[lo]! + dh01 * ys[lo + 1]! + dh11 * slopes[lo + 1]!;122}123