spb/earth-now Public License

1/**2 * earth-now.co3 * Author:  Simon-Pierre Boucher4 * Contact: contact@spboucher.ai5 * File:    packages/models/src/fourier.ts6 * Purpose: Least-squares fit of trend + Fourier harmonics in the VALUE domain, converted to rate harmonics for CounterModel7 */89import {10  type Harmonic,11  DAY_SECONDS,12  SEASONAL_EPOCH_MS,13  WEEK_SECONDS,14  YEAR_SECONDS,15  parseIsoUtc,16} from "@earth-now/counter";17import { leastSquares } from "./linalg.js";1819export interface FourierFitInput {20  /** Observations: ISO time → value. Needs >= 2·harmonics + 2 points. */21  observations: Array<{ time: string; value: number }>;22  /** Harmonic spec to fit, e.g. [{period:"year", order:1}, {period:"year", order:2}]. */23  harmonics: Array<{ period: "year" | "week" | "day"; order: number }>;24}2526export interface FourierFitResult {27  /** Value at the seasonal epoch (intercept). */28  intercept: number;29  /** Linear trend of the value, per second. Becomes rateFn.base. */30  trendPerSecond: number;31  /** Rate-domain harmonics ready for a `seasonal` RateFunction. */32  rateHarmonics: Harmonic[];33  /** Evaluate the fitted VALUE curve at t (ms). */34  valueAt: (tMs: number) => number;35  /** Root-mean-square residual of the fit. */36  rmse: number;37}3839function omegaOf(period: "year" | "week" | "day", order: number): number {40  const seconds =41    period === "year" ? YEAR_SECONDS : period === "week" ? WEEK_SECONDS : DAY_SECONDS;42  return (2 * Math.PI * order) / seconds;43}4445/**46 * Fit value(τ) = a + b·τ + Σ_h [ c_h·cos(ω_h τ) + s_h·sin(ω_h τ) ], τ in seconds47 * since the seasonal epoch. The derivative of each fitted value harmonic is an48 * exact rate harmonic, so the returned Harmonics reproduce the fitted curve when49 * integrated by packages/counter — the client animates the very curve we fitted.50 */51export function fitFourier(input: FourierFitInput): FourierFitResult {52  const { observations, harmonics } = input;53  const k = 2 + 2 * harmonics.length;54  if (observations.length < k) {55    throw new Error(56      `fitFourier: need at least ${k} observations for ${harmonics.length} harmonics, got ${observations.length}`,57    );58  }5960  const taus = observations.map((o) => (parseIsoUtc(o.time) - SEASONAL_EPOCH_MS) / 1000);61  // Center τ for conditioning; un-center the intercept afterwards.62  const tauMean = taus.reduce((s, t) => s + t, 0) / taus.length;63  const omegas = harmonics.map((h) => omegaOf(h.period, h.order));6465  const X = taus.map((tau) => {66    const row = [1, tau - tauMean];67    for (const omega of omegas) row.push(Math.cos(omega * tau), Math.sin(omega * tau));68    return row;69  });70  const y = observations.map((o) => o.value);71  const beta = leastSquares(X, y);7273  const b = beta[1]!;74  const a = beta[0]! - b * tauMean;75  const coefs = harmonics.map((h, i) => ({76    ...h,77    c: beta[2 + 2 * i]!,78    s: beta[3 + 2 * i]!,79    omega: omegas[i]!,80  }));8182  const valueAt = (tMs: number): number => {83    const tau = (tMs - SEASONAL_EPOCH_MS) / 1000;84    let v = a + b * tau;85    for (const { c, s, omega } of coefs) v += c * Math.cos(omega * tau) + s * Math.sin(omega * tau);86    return v;87  };8889  let sse = 0;90  for (let i = 0; i < taus.length; i++) {91    const r = valueAt(parseIsoUtc(observations[i]!.time)) - y[i]!;92    sse += r * r;93  }9495  // d/dτ [c·cos(ωτ) + s·sin(ωτ)] = (s·ω)·cos(ωτ) − (c·ω)·sin(ωτ)96  //                              = A·ω·cos(ωτ + φ) with A = √(c²+s²), φ = atan2(c, s)... derived below.97  const rateHarmonics: Harmonic[] = coefs.map(({ period, order, c, s, omega }) => {98    const amplitude = Math.hypot(c, s) * omega;99    // Write derivative as R·cos(ωτ + φ): R·cosφ = s·ω, R·sinφ = c·ω ⇒ φ = atan2(c, s).100    const phase = Math.atan2(c, s);101    return { period, order, amplitude, phase };102  });103104  return {105    intercept: a,106    trendPerSecond: b,107    rateHarmonics,108    valueAt,109    rmse: Math.sqrt(sse / taus.length),110  };111}112