// Is the large-pluck motion of the 5-bead string chaotic, or a solver artifact?
//
// Runs the post's large pluck (n=5, tension 5.001 N, bead 2, 10 mm, 2 s) with the app's own
// model and CPU solver (double precision), and measures
//   1. step size: runs at cfl 0.1, 0.05, ..., halving; each against the next finer one.
//      For chaos with a convergent solver, the runs agree up to a time that grows by about
//      ln(4)/lambda per halving (velocity Verlet is second order, so the error drops 4x).
//   2. perturbations: pluck amplitude times (1 + eps), at two step sizes. The separation
//      should grow exponentially at the same rate lambda (the largest Lyapunov exponent)
//      at both step sizes, and its divergence time should grow like ln(1/eps)/lambda.
//   3. energy balance error of every run (velocity Verlet nearly conserves energy).
//   4. statistics: band levels of the sound (octave bands) at every step size.
//   5. the gentle pluck (0.1 mm) against the linear modes.
//   6. aliasing: the solver point-samples the pickup at fs (48 kHz) without a low-pass
//      filter, so sound above fs/2 folds back into the spectrogram. Runs at 16 fs and
//      reports the fraction of the pickup's power above 24 kHz.
//
// Positions are saved at every audio sample (frameRate = fs), since the runs separate
// within milliseconds; the replay rate does not affect the dynamics.
//
//   npx tsx analysis/chaos-check.ts [out.json]      (about a minute per fine run)

import { writeFileSync } from 'node:fs';
import { CONFIGS } from '../src/instruments/configs';
import { parseSettings, simOptions } from '../src/ui/settings';
import { simulate } from '../src/sim/solver';
import { spectrum } from '../src/ui/fft';
import type { SimResult } from '../src/sim/types';

const BASE = 'n=5&tension=5.001&excBead=2';
const AMP = 0.01;

interface Run {
  label: string;
  cfl: number;
  amp: number;
  r: SimResult;
}

function run(label: string, amp: number, cfl: number, over: { fs?: number; frameRate?: number } = { frameRate: 48000 }): Run {
  const cfg = CONFIGS['bead-string'];
  const s = parseSettings(cfg, `${BASE}&pluckAmp=${amp}&cfl=${cfl}`);
  // set directly: a link only accepts the app's sample rates
  Object.assign(s.params, over);
  const built = cfg.setup(s.params).built!;
  const r = simulate(built.sys, built.exc, simOptions(s.params, built));
  const emax = Math.max(...r.energy);
  let bal = 0;
  for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
  console.error(`${label}: dt ${r.dt.toExponential(3)} s (${r.stepsPerSample}/sample), ${r.wallSeconds.toFixed(1)} s wall, energy error ${((bal / emax) * 100).toExponential(2)}% of peak`);
  return { label, cfl, amp, r };
}

// RMS over beads of the displacement difference at each frame (m)
function separation(a: SimResult, b: SimResult): Float64Array {
  const nf = a.frameTimes.length;
  const w = a.frames.length / nf;
  const d = new Float64Array(nf);
  for (let f = 0; f < nf; f++) {
    let s = 0;
    for (let i = 0; i < w; i++) {
      const x = a.frames[f * w + i] - b.frames[f * w + i];
      s += x * x;
    }
    d[f] = Math.sqrt(s / (w / 3));
  }
  return d;
}

// first time the separation exceeds thr (s), or null
function crossing(t: Float64Array, d: Float64Array, thr: number): number | null {
  for (let i = 0; i < d.length; i++) if (d[i] > thr) return t[i];
  return null;
}

// growth rate of the running maximum of ln d between lo and hi (1/s): least squares
function growthRate(t: Float64Array, d: Float64Array, lo: number, hi: number): { lambda: number; t0: number; t1: number } | null {
  const xs: number[] = [];
  const ys: number[] = [];
  let m = 0;
  for (let i = 0; i < d.length; i++) {
    m = Math.max(m, d[i]);
    if (m >= lo && m <= hi) {
      xs.push(t[i]);
      ys.push(Math.log(m));
    }
  }
  if (xs.length < 10) return null;
  const n = xs.length;
  const mx = xs.reduce((s, x) => s + x, 0) / n;
  const my = ys.reduce((s, y) => s + y, 0) / n;
  let sxy = 0;
  let sxx = 0;
  for (let i = 0; i < n; i++) {
    sxy += (xs[i] - mx) * (ys[i] - my);
    sxx += (xs[i] - mx) ** 2;
  }
  return { lambda: sxy / sxx, t0: xs[0], t1: xs[n - 1] };
}

// sound level in octave bands (dB re the loudest band), from 1 s on (the settled motion)
function bands(r: SimResult): { f: number; db: number }[] {
  const start = Math.round(1 * r.fs);
  const s = spectrum(r.audio.subarray(start), r.fs);
  const out: { f: number; db: number }[] = [];
  for (let fc = 125; fc <= 16000; fc *= 2) {
    let p = 0;
    for (let k = 0; k < s.freq.length; k++) if (s.freq[k] >= fc / Math.SQRT2 && s.freq[k] < fc * Math.SQRT2) p += 10 ** (s.db[k] / 10);
    out.push({ f: fc, db: 10 * Math.log10(p) });
  }
  const top = Math.max(...out.map((b) => b.db));
  return out.map((b) => ({ f: b.f, db: b.db - top }));
}

const downsample = (x: Float64Array, k: number) => Array.from(x.filter((_, i) => i % k === 0));

const out: Record<string, unknown> = { settings: `${BASE}&pluckAmp=${AMP}` };

// 1. step size
const CFLS = [0.1, 0.05, 0.025, 0.0125, 0.00625, 0.003125];
const byCfl = CFLS.map((c) => run(`cfl ${c}`, AMP, c));
const t = byCfl[0].r.frameTimes;
const DS = 24; // keep every 24th frame (0.5 ms) in the output
out.time = downsample(t, DS);
out.stepSize = byCfl.slice(0, -1).map((a, i) => {
  const b = byCfl[i + 1];
  const d = separation(a.r, b.r);
  return {
    cfl: a.cfl,
    against: b.cfl,
    dt: a.r.dt,
    divergeAt1um: crossing(t, d, 1e-6),
    divergeAt1mm: crossing(t, d, 1e-3),
    d: downsample(d, DS),
  };
});

// 2. perturbations at two step sizes
const EPS = [1e-5, 1e-7, 1e-9, 1e-11];
const pert: Record<string, unknown>[] = [];
for (const [ci, cfl] of [[0, 0.1], [2, 0.025]] as const) {
  const base = byCfl[ci];
  for (const eps of EPS) {
    const p = run(`cfl ${cfl}, eps ${eps}`, AMP * (1 + eps), cfl);
    const d = separation(base.r, p.r);
    pert.push({
      cfl,
      eps,
      divergeAt1mm: crossing(t, d, 1e-3),
      growth: growthRate(t, d, 1e-8, 1e-4),
      bands: bands(p.r),
      d: downsample(d, DS),
    });
  }
}
out.perturbation = pert;

// 3. energy balance and 4. sound statistics per step size
out.runs = byCfl.map(({ cfl, r }) => {
  const emax = Math.max(...r.energy);
  let bal = 0;
  for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
  return { cfl, dt: r.dt, stepsPerSample: r.stepsPerSample, wallSeconds: r.wallSeconds, energyError: bal / emax, maxStretch: r.maxStretch, bands: bands(r) };
});

// 5. gentle pluck vs the linear modes: f_k = (1/pi) sqrt(T / (m a)) sin(k pi / 12)
const gentle = run('gentle 0.1 mm', 0.0001, 0.1);
{
  const s = spectrum(gentle.r.audio, gentle.r.fs);
  const T = 5.001;
  const m = 5.08e-5;
  const a = 0.0127;
  const linear = [1, 2, 3, 4, 5].map((k) => ({ k, hz: (1 / Math.PI) * Math.sqrt(T / (m * a)) * Math.sin((k * Math.PI) / 12) }));
  const top = Math.max(...s.db);
  out.gentle = linear.map(({ k, hz }) => {
    // strongest bin within 3% of the predicted frequency
    let best = -1;
    for (let i = 0; i < s.freq.length; i++) if (Math.abs(s.freq[i] - hz) < 0.03 * hz && (best < 0 || s.db[i] > s.db[best])) best = i;
    return { k, predictedHz: hz, peakHz: s.freq[best], levelDb: s.db[best] - top };
  });
}

// 6. aliasing
out.aliasing = [0.001, AMP].map((amp) => {
  const { r } = run(`${amp * 1000} mm at 16 fs`, amp, 0.1, { fs: 768000, frameRate: 2000 });
  const s = spectrum(r.audio, r.fs);
  let total = 0;
  let above = 0;
  let above20 = 0;
  for (let k = 0; k < s.freq.length; k++) {
    const p = 10 ** (s.db[k] / 10);
    total += p;
    if (s.freq[k] > 24000) above += p;
    if (s.freq[k] > 20000) above20 += p;
  }
  // octave bands up to 384 kHz, dB re the loudest
  const b: { f: number; db: number }[] = [];
  for (let fc = 125; fc <= 256000; fc *= 2) {
    let p = 0;
    for (let k = 0; k < s.freq.length; k++) if (s.freq[k] >= fc / Math.SQRT2 && s.freq[k] < fc * Math.SQRT2) p += 10 ** (s.db[k] / 10);
    b.push({ f: fc, db: 10 * Math.log10(p) });
  }
  const top = Math.max(...b.map((x) => x.db));
  return { amp, fractionAbove24k: above / total, fractionAbove20k: above20 / total, bands: b.map((x) => ({ f: x.f, db: x.db - top })) };
});

const file = process.argv[2] ?? 'chaos-check.json';
writeFileSync(file, JSON.stringify(out));
console.error(`wrote ${file}`);