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";

/** Relative rms difference over the common length. */
function relRms(a: ArrayLike<number>, b: ArrayLike<number>): number {
  let num = 0, den = 0;
  for (let i = 0; i < Math.min(a.length, b.length); i++) { num += (a[i] - b[i]) ** 2; den += b[i] ** 2; }
  return Math.sqrt(num / den);
}

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, { recordAcceleration: true });
  const front: number[] = [], back: number[] = [], acc: Float32Array[] = [], N = setup.mesh.n;
  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);
    for (let i = 0; i < r.pressureFront.length; i++) acc.push(r.acceleration!.slice(i * N, (i + 1) * N));
  }
  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 and the acceleration record: the same decimation kernel and normal component as Python;
  // the delay is applied at sub-step resolution there and by a 16-tap interpolator here
  const pu = front.map((f, i) => f - back[i]), ref: number[] = py.pressure_unbaffled;
  assert.ok(relRms(pu, ref) < 1e-3, `unbaffled rel rms error ${relRms(pu, ref)}`);
  assert.ok(relRms(front, py.pressure_rayleigh) < 1e-3, `baffled rel rms error ${relRms(front, py.pressure_rayleigh)}`);
  py.acceleration_nodes.forEach((nd: number, k: number) => {
    const mine = acc.map((row) => row[nd]);
    assert.ok(relRms(mine, py.acceleration[k]) < 1e-3, `acceleration of node ${nd}: rel rms error ${relRms(mine, py.acceleration[k])}`);
  });
});

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]);
  assert.ok(relRms(pu, py.pressure_unbaffled) < 1e-3, `unbaffled rel rms error ${relRms(pu, py.pressure_unbaffled)}`);
});