import test from "node:test";
import assert from "node:assert/strict";
import { buildSetup } from "../src/model/setup.ts";
import { CpuSimulator } from "../src/cpu/simulator.ts";
import { fixture } from "./helpers.ts";

test("CPU time loop reproduces the Python run (dent-free, float64)", async () => {
  const setup = await buildSetup(fixture.config);
  assert.equal(setup.substeps, fixture.timestep.substeps);
  assert.ok(Math.abs(setup.dt / fixture.timestep.dt - 1) < 1e-12);
  const sim = new CpuSimulator(setup);
  const front: number[] = [], back: number[] = [];
  const py = fixture.run, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
  const nSamples = Math.round(last[0] * setup.fs);   // Python ran whole blocks, slightly past the duration
  const chunk = 147;
  for (let done = 0; done < nSamples; done += chunk) {
    const r = await sim.runChunk(Math.min(chunk, nSamples - done));
    front.push(...r.pressureFront); back.push(...r.pressureBack);
  }
  const d = await sim.diagnostics();
  assert.ok(Math.abs(d.time - last[0]) < 1e-9);
  const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
  assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-3, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
  assert.ok(Math.abs(d.hits[0].contactTime - hit.contact_time) <= 2 * setup.dt + 1e-12);
  assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-3, `rebound ${d.hits[0].velocity} vs ${hit.rebound_velocity}`);
  // the run is strongly nonlinear (max|w|/h ~ 10), so small scheme differences are amplified: 1 percent
  assert.ok(rel(d.kinetic, last[1]) < 1e-2 && rel(d.membrane, last[2]) < 1e-2 && rel(d.bending, last[3]) < 1e-2,
            `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
  assert.ok(rel(d.exciters, last[4]) < 1e-2 && rel(d.maxWoverH, last[5]) < 1e-2);
  // the lagged projection leaves a rigid velocity of one step's increment: tiny against the vibration
  assert.ok(sim.rigidContent() < 1e-3 * Math.sqrt(2 * d.kinetic), `rigid content ${sim.rigidContent()}`);
  // radiation: a different decimation and delay scheme, so compare loosely (rms and correlation)
  const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
  const n = Math.min(pu.length, ref.length) - 4;
  const rms = (a: number[]) => Math.sqrt(a.slice(0, n).reduce((s, x) => s + x * x, 0) / n);
  let best = 0;
  for (let lag = -3; lag <= 3; lag++) {
    let dot = 0;
    for (let i = 4; i < n; i++) dot += pu[i] * ref[i + lag];
    best = Math.max(best, dot / (rms(pu) * rms(ref) * (n - 4)));
  }
  assert.ok(Math.abs(rms(pu) / rms(ref) - 1) < 0.1, `rms ${rms(pu)} vs ${rms(ref)}`);
  assert.ok(best > 0.95, `correlation ${best}`);
});

test("linear run matches Python tightly (no chaotic amplification)", async () => {
  const cfg = structuredClone(fixture.config);
  cfg.numerics.nonlinear = false;
  const setup = await buildSetup(cfg);
  const sim = new CpuSimulator(setup);
  const py = fixture.run_linear, hit = py.hits[0], last = py.energy_log[py.energy_log.length - 1];
  const nSamples = Math.round(last[0] * setup.fs);
  const front: number[] = [], back: number[] = [];
  for (let done = 0; done < nSamples; done += 200) {
    const r = await sim.runChunk(Math.min(200, nSamples - done));
    front.push(...r.pressureFront); back.push(...r.pressureBack);
  }
  const d = await sim.diagnostics();
  const rel = (a: number, b: number) => Math.abs(a - b) / Math.abs(b);
  assert.ok(rel(d.hits[0].peakForce, hit.peak_force) < 1e-5, `peak ${d.hits[0].peakForce} vs ${hit.peak_force}`);
  assert.ok(rel(d.hits[0].velocity, hit.rebound_velocity) < 1e-5);
  assert.ok(rel(d.kinetic, last[1]) < 1e-4 && rel(d.membrane, last[2]) < 1e-4 && rel(d.bending, last[3]) < 1e-4,
            `energies ${[d.kinetic, d.membrane, d.bending]} vs ${last.slice(1, 4)}`);
  assert.ok(rel(d.maxWoverH, last[5]) < 1e-4);
  const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
  const n = Math.min(pu.length, ref.length) - 4;
  const rms = (a: number[]) => Math.sqrt(a.slice(0, n).reduce((s, x) => s + x * x, 0) / n);
  assert.ok(Math.abs(rms(pu) / rms(ref) - 1) < 0.1, `rms ${rms(pu)} vs ${rms(ref)}`);
});