import { readFileSync, readdirSync } from "fs";
import { executeCode } from "numbl";
import { describe, expect, it } from "vitest";
import { computeModes, modalDamping } from "../src/sim/modes";
import { numblEig, runScript } from "../src/sim/runScript";
import { growthRate } from "../src/sim/analysis";
import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
import { runTwin } from "../src/sim/twin";
import type { BeadSystem } from "../src/sim/types";

const helpers = readdirSync("matlab/helpers").map(name => ({
  name,
  source: readFileSync(`matlab/helpers/${name}`, "utf8"),
}));
const build = (src: string) => runScript(src, helpers, () => {});

// A string: N free beads between two pinned ends, tension f0, spacing a.
// Transverse modes: omega_k = 2 sqrt(f0/(m a)) sin(k pi / (2(N+1))).
const STRING = (extra: string) => `
N = 15; a = 0.01; m = 1e-3; f0 = 10; k = 1e4;
[pos, E] = bp_chain(N + 2, (N + 1) * a);
sys.pos = pos;
sys.mass = m;
sys.pinned = [1, N + 2];
sys.springs = bp_springs(pos, E, k, f0);
hit.beads = 5; hit.velocity = [0 0 1e-4];
sys.hits = hit;
sys.listen = struct('beads', 4);
sys.duration = 0.2;
${extra}
`;

function peakFrequency(x: Float32Array, fs: number, fLo: number, fHi: number): number {
  // direct DFT magnitude on a fine grid, refined by a parabola
  let best = 0;
  let bestF = 0;
  const mag = (f: number) => {
    let re = 0;
    let im = 0;
    for (let i = 0; i < x.length; i++) {
      const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (x.length - 1));
      re += w * x[i] * Math.cos((2 * Math.PI * f * i) / fs);
      im -= w * x[i] * Math.sin((2 * Math.PI * f * i) / fs);
    }
    return Math.hypot(re, im);
  };
  for (let f = fLo; f <= fHi; f += 0.5) {
    const m = mag(f);
    if (m > best) {
      best = m;
      bestF = f;
    }
  }
  const [a, b, c] = [mag(bestF - 0.5), best, mag(bestF + 0.5)];
  return bestF + (0.5 * 0.5 * (a - c)) / (a - 2 * b + c);
}

describe("string", () => {
  const { sys } = build(STRING(""));
  const f1 = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin(Math.PI / 32)) / (2 * Math.PI);

  it("has the analytic transverse modes", () => {
    const modes = computeModes(sys, numblEig);
    for (let k = 1; k <= 15; k++) {
      const fk = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin((k * Math.PI) / 32)) / (2 * Math.PI);
      expect(modes.freqs[k - 1]).toBeCloseTo(fk, 6);
    }
  });

  it("rings at the fundamental for a small strike", () => {
    const s: BeadSystem = { ...sys, hits: [{ ...sys.hits[0], beads: Int32Array.of(8) }], listenBeads: Int32Array.of(8) };
    const r = simulate(s);
    expect(r.aborted).toBeNull();
    const f = peakFrequency(r.audio, r.sampleRate, f1 * 0.8, f1 * 1.2);
    expect(Math.abs(f - f1) / f1).toBeLessThan(2e-3);
  });

  it("is at equilibrium at rest", () => {
    const r = simulate({ ...sys, duration: 0.001 });
    expect(r.residualForce).toBeLessThan(1e-9);
    expect(r.warnings).toEqual([]);
  });
});

describe("energy", () => {
  const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
  const drift = (method: "verlet" | "yoshida4", dt: number) =>
    simulate({ ...sys, method, dt, duration: 0.05 }).energyDrift;

  it("is conserved to second order with Verlet", () => {
    const d1 = drift("verlet", 1 / 48000 / 8);
    const d2 = drift("verlet", 1 / 48000 / 16);
    expect(d1).toBeLessThan(1e-3);
    expect(d1 / d2).toBeGreaterThan(3);
    expect(d1 / d2).toBeLessThan(5);
  });

  it("is conserved to fourth order with Yoshida", () => {
    const d1 = drift("yoshida4", 1 / 48000 / 8);
    const d2 = drift("yoshida4", 1 / 48000 / 16);
    expect(d1 / d2).toBeGreaterThan(12);
  });

  it("balances with the dissipated energy when damped", () => {
    const { sys: damped } = build(
      STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 2e-3; sys.dashpots = [E, 1e-3*ones(size(E,1),1)];")
    );
    // The damping step's energy loss is counted exactly, so E + dissipated
    // drifts only by the conservative integrator's O(h^2) error.
    const r1 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 8 });
    const r2 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 16 });
    expect(r1.dissipated[r1.dissipated.length - 1]).toBeGreaterThan(0.3 * r1.energy[0]);
    expect(r1.energyDrift).toBeLessThan(1e-4);
    expect(r1.energyDrift / r2.energyDrift).toBeGreaterThan(3);
  });
});

describe("damping", () => {
  it("matches friction when a dashpot ties a bead to a pinned bead", () => {
    // one free bead between two pinned beads: friction c0 on it is the same
    // force as two dashpots of c0/2 to its pinned neighbours
    const base = STRING("").replace("N = 15", "N = 1");
    const { sys: a } = build(base + "sys.friction = [0; 0.02; 0]; sys.hits.beads = 2; sys.listen.beads = 2;");
    const { sys: b } = build(base + "sys.dashpots = [1 2 0.01; 2 3 0.01]; sys.hits.beads = 2; sys.listen.beads = 2;");
    const ra = simulate({ ...a, duration: 0.05 });
    const rb = simulate({ ...b, duration: 0.05 });
    let err = 0;
    let peak = 0;
    for (let i = 0; i < ra.audio.length; i++) {
      err = Math.max(err, Math.abs(ra.audio[i] - rb.audio[i]));
      peak = Math.max(peak, Math.abs(ra.audio[i]));
    }
    expect(err / peak).toBeLessThan(1e-6);
  });

  it("decays a mode at the predicted rate", () => {
    // one free bead, friction c0: amplitude decays as exp(-c0/(2m) t)
    const base = STRING("").replace("N = 15", "N = 1");
    const { sys } = build(base + "sys.friction = 0.02; sys.hits.beads = 2; sys.listen.beads = 2;");
    const r = simulate({ ...sys, duration: 0.2 });
    const rms = (from: number, to: number) => {
      let s = 0;
      for (let i = from; i < to; i++) s += r.audio[i] ** 2;
      return Math.sqrt(s / (to - from));
    };
    const fs = r.sampleRate;
    const ratio = rms(0.15 * fs, 0.2 * fs) / rms(0.05 * fs, 0.1 * fs);
    const expected = Math.exp((-0.02 / (2 * 1e-3)) * 0.1);
    expect(Math.abs(ratio - expected) / expected).toBeLessThan(0.02);
  });
});

describe("integrator", () => {
  it("reports Verlet's known frequency error", () => {
    // Verlet: omega_num h = 2 asin(omega h / 2)
    const w = 2 * Math.PI * 5000;
    const h = 1e-5;
    const exact = 1200 * Math.log2((2 * Math.asin((w * h) / 2)) / (w * h));
    expect(frequencyErrorCents("verlet", w, h)).toBeCloseTo(exact, 9);
  });

  it("has the expected stability limits", () => {
    expect(stabilityLimit("verlet")).toBe(2);
    expect(stabilityLimit("yoshida4")).toBeGreaterThan(1.5);
    expect(stabilityLimit("yoshida4")).toBeLessThan(1.6);
  });
});

describe("script errors", () => {
  it("names the field and bead", () => {
    expect(() => build(STRING("sys.pinned = [1 99];"))).toThrow(/sys.pinned: bead index 99/);
    expect(() => build(STRING("sys.springs = [1 2 3];"))).toThrow(/4 or 5 columns/);
  });
});

describe("chaos example", () => {
  const src = readFileSync("matlab/examples/chaos.m", "utf8");
  const rate = (strength: number) => {
    const { sys } = build(src.replace(/^s = 1;/m, `s = ${strength};`));
    const s1 = { ...sys, duration: 1, frameRate: 1000 };
    const saved: Float64Array[] = [];
    const r = simulate(s1, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
    const { sep } = runTwin(s1, saved, r.dt, 1e-9);
    return growthRate(r.frameTimes, sep, 1e-9);
  };

  it("is chaotic when struck hard", () => {
    const g = rate(1);
    expect(g).not.toBeNull();
    expect(g!).toBeGreaterThan(50); // e-folding time under 20 ms
  });

  it("is regular when struck softly", () => {
    expect(rate(0.025)).toBeNull();
  });
});

describe("bp_decay", () => {
  // String struck in the shape of mode p; its amplitude should decay at
  // sigma0 + sigma1 w_p^2, the law through T60 = 0.5 s at 100 Hz and 0.1 s at 1 kHz.
  const shaped = (p: number) =>
    build(
      STRING(`
sys.hits = [];
for i = 1:N
  sys.hits(i).t = 0;
  sys.hits(i).beads = i + 1;
  sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
end
sys = bp_decay(sys, [100 0.5; 1000 0.1]);
sys.listen.beads = 2:N+1;
sys.listen.dir = [0 0 1];
`).replace("sys.listen = struct('beads', 4);", "")
    );

  it("gives each mode the decay rate of the chosen law", () => {
    const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
    const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
    const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
    const sigma0 = s[0] - sigma1 * w2[0];
    const { sys } = shaped(1);
    const modes = computeModes(sys, numblEig);
    const md = modalDamping(sys, modes);
    for (let k = 0; k < modes.freqs.length; k++) {
      const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
      expect(Math.abs(md[k] - expected) / expected).toBeLessThan(1e-9);
    }
    for (const p of [1, 3]) {
      const { sys: sp } = shaped(p);
      const r = simulate({ ...sp, duration: 0.3, frameRate: 1000 });
      // energy decays at 2 sigma; compare energy at 0.05 s and 0.25 s (averaged over a period)
      const avg = (t: number) => {
        let e = 0;
        let c = 0;
        for (let i = 0; i < r.frameTimes.length; i++) if (Math.abs(r.frameTimes[i] - t) < 0.01) (e += r.energy[i]), c++;
        return e / c;
      };
      const rate = Math.log(avg(0.05) / avg(0.25)) / 0.2 / 2;
      const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[p - 1]) ** 2;
      expect(Math.abs(rate - expected) / expected).toBeLessThan(0.02);
    }
  });

  it("rejects a law the dashpots cannot make", () => {
    expect(() => build(STRING("sys = bp_decay(sys, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
  });
});

describe("default damping", () => {
  it("makes every mode fall 60 dB in 2.5 s when the script sets none", () => {
    const { sys, warnings } = build(STRING(""));
    expect(warnings.join(" ")).toMatch(/default/);
    const modes = computeModes(sys, numblEig);
    for (const s of modalDamping(sys, modes)) expect(Math.log(1000) / s).toBeCloseTo(2.5, 9);
  });

  it("is off when the script sets friction = 0", () => {
    const { sys } = build(STRING("sys.friction = 0;"));
    expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
  });
});

describe("bending", () => {
  // Stiff string: tension plus bending triplets. With theta = k pi/(N+1),
  // omega_k^2 = [g (2 - 2 cos theta) + kb (2 - 2 cos theta)^2] / m exactly.
  const STIFF = (extra: string) =>
    STRING(`
[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
${extra}`);

  it("gives the stiff string its exact modes", () => {
    const { sys } = build(STIFF(""));
    const modes = computeModes(sys, numblEig);
    const g = 10 / 0.01;
    for (let k = 1; k <= 15; k++) {
      const c = 2 - 2 * Math.cos((k * Math.PI) / 16);
      const fk = Math.sqrt((g * c + 50 * c * c) / 1e-3) / (2 * Math.PI);
      expect(modes.freqs[k - 1] / fk).toBeCloseTo(1, 10);
    }
  });

  it("conserves energy to second order", () => {
    const { sys } = build(STIFF("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
    const d1 = simulate({ ...sys, dt: 1 / 48000 / 8, duration: 0.05 }).energyDrift;
    const d2 = simulate({ ...sys, dt: 1 / 48000 / 16, duration: 0.05 }).energyDrift;
    expect(d1).toBeLessThan(1e-3);
    expect(d1 / d2).toBeGreaterThan(3);
    expect(d1 / d2).toBeLessThan(5);
  });

  it("follows the bp_decay law with bending damping", () => {
    const { sys } = build(STIFF("sys = bp_decay(sys, [100 0.5; 1000 0.1]);"));
    const modes = computeModes(sys, numblEig);
    const md = modalDamping(sys, modes);
    const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
    const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
    const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
    const sigma0 = s[0] - sigma1 * w2[0];
    for (let k = 0; k < modes.freqs.length; k++) {
      const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
      expect(md[k] / expected).toBeCloseTo(1, 9);
    }
    // and the solver's damping step agrees: energy balance still exact
    const r1 = simulate({ ...sys, hits: [{ ...sys.hits[0], velocity: [0, 0.3, 2] }], dt: 1 / 48000 / 8, duration: 0.05 });
    expect(r1.energyDrift).toBeLessThan(1e-3);
  });

  it("leaves a straight chain in equilibrium after a rigid rotation", () => {
    const { sys } = build(
      STIFF("R3 = [cos(0.7) -sin(0.7) 0; sin(0.7) cos(0.7) 0; 0 0 1] * [1 0 0; 0 cos(0.4) -sin(0.4); 0 sin(0.4) cos(0.4)]; sys.pos = pos * R3';")
    );
    const r = simulate({ ...sys, duration: 0.001 });
    expect(r.residualForce).toBeLessThan(1e-9);
  });
});

describe("free objects", () => {
  it("bp_anchor holds a bead with the chosen out-of-plane stiffness, at equilibrium", () => {
    const { sys } = build(`
[pos, springs, anchors] = bp_anchor([0 0 0], 1, 400, 0.02);
sys.pos = pos; sys.mass = 1e-3; sys.pinned = anchors; sys.springs = springs;
hit.beads = 1; hit.velocity = [0 0 1e-4]; sys.hits = hit; sys.listen = struct('beads', 1);
sys.friction = 0; sys.duration = 0.01;`);
    const modes = computeModes(sys, numblEig);
    expect(modes.freqs[0]).toBeCloseTo(Math.sqrt(400 / 1e-3) / (2 * Math.PI), 8);
    expect(simulate(sys).residualForce).toBeLessThan(1e-9);
  });

  it("bar example has the overtones of a free bar", () => {
    const { sys } = build(readFileSync("matlab/examples/bar.m", "utf8"));
    const f = computeModes(sys, numblEig).freqs;
    // two low modes are the bar bouncing on its supports
    expect(f[1]).toBeLessThan(100);
    const f1 = f[2];
    [2.756, 5.404, 8.933].forEach((ratio, q) => expect(Math.abs(f[3 + q] / f1 / ratio - 1)).toBeLessThan(0.01));
  });
});

describe("listener quantity", () => {
  // Stiff, heavily damped string (dashpots and bending damping), listened to
  // as displacement, velocity and acceleration at the same bead.
  const run = (q: string) => {
    const { sys } = build(
      STRING(`
[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
sys = bp_decay(sys, [100 0.05; 1000 0.02]);
sys.hits.velocity = [0 0 1e-4];
sys.listen.quantity = '${q}';`)
    );
    return simulate({ ...sys, duration: 0.05 });
  };

  it("defaults to acceleration", () => {
    expect(build(STRING("")).sys.listenQuantity).toBe("acceleration");
  });

  it("gives acceleration and velocity consistent with each other and the displacement", () => {
    const x = run("displacement").audio;
    const v = run("velocity").audio;
    const acc = run("acceleration").audio;
    const fs = 48000;
    const rel = (p: (i: number) => number, ref: Float32Array) => {
      let num = 0;
      let den = 0;
      for (let i = 1; i < ref.length - 1; i++) {
        num += (p(i) - ref[i]) ** 2;
        den += ref[i] ** 2;
      }
      return Math.sqrt(num / den);
    };
    // central differences of the recorded signals; their own error is second
    // order in 1/fs (about 5e-4 for v -> a here, and 1e-5 for x -> v, which is
    // limited by the single-precision audio)
    expect(rel(i => ((x[i + 1] - x[i - 1]) * fs) / 2, v)).toBeLessThan(1e-4);
    expect(rel(i => ((v[i + 1] - v[i - 1]) * fs) / 2, acc)).toBeLessThan(1e-3);
  });

  it("rejects an unknown quantity", () => {
    expect(() => build(STRING("sys.listen.quantity = 'jerk';"))).toThrow(/quantity/);
  });
});

describe("hinges", () => {
  // A small dome: hexagonal lattice lifted to a paraboloid, clamped at the top,
  // with stretching springs (no tension) and hinges across interior edges.
  const DOME = (extra: string, H = 0.01) => `
R = 0.06; a = 0.012; H = ${H}; m = 1e-3; kh = 0.5; k = 2e4;
[pos, tri] = bp_hexdisk(R, a);
pos(:, 3) = H * (1 - sum(pos(:, 1:2).^2, 2) / R^2);
E = bp_edges(tri);
sys.pos = pos; sys.mass = m;
sys.pinned = bp_near(pos, [0 0 H], 1.05 * a);
sys.springs = bp_springs(pos, E, k, 0);
hh = bp_hinges(tri);
sys.hinges = [hh, kh * ones(size(hh, 1), 1)];
hit.beads = bp_near(pos, [0.04 0 0], 0); hit.velocity = [0 0 1e-3]; sys.hits = hit;
sys.listen = struct('beads', bp_near(pos, [-0.03 0.02 0], 0));
sys.friction = 0; sys.duration = 0.05;
${extra}`;

  it("bp_hinges lists every interior edge once, with the two opposite beads", () => {
    const r = executeCode(
      `[pos, tri] = bp_hexdisk(0.06, 0.012); E = bp_edges(tri); H = bp_hinges(tri);
       nint = size(E, 1) - sum(accumarray(sort([tri(:,[1 2]); tri(:,[2 3]); tri(:,[3 1])], 2) * [1e4; 1], 1) == 1);`,
      {},
      helpers,
      "t.m"
    );
    const H = r.variableValues.H as unknown as { data: Float64Array; shape: number[] };
    const tri = r.variableValues.tri as unknown as { data: Float64Array; shape: number[] };
    const nh = H.shape[0];
    const nt = tri.shape[0];
    const tris = new Set<string>();
    for (let t = 0; t < nt; t++) tris.add([0, 1, 2].map(c => tri.data[t + c * nt]).sort((p, q) => p - q).join(","));
    for (let h = 0; h < nh; h++) {
      const [i, j, k, l] = [0, 1, 2, 3].map(c => H.data[h + c * nh]);
      expect(tris.has([i, j, k].sort((p, q) => p - q).join(","))).toBe(true);
      expect(tris.has([l, j, k].sort((p, q) => p - q).join(","))).toBe(true);
    }
    expect(nh).toBe(r.variableValues.nint);
  });

  it("leaves a dome, and a rotated dome, unstressed at rest", () => {
    for (const rot of ["", "Q = [cos(0.7) 0 sin(0.7); 0 1 0; -sin(0.7) 0 cos(0.7)]; sys.pos = sys.pos * Q';"]) {
      const { sys } = build(DOME(rot));
      const r = simulate({ ...sys, duration: 0.001 });
      expect(r.residualForce).toBeLessThan(1e-9);
      expect(r.warnings).toEqual([]);
    }
  });

  it("conserves energy to second order", () => {
    const { sys } = build(DOME("sys.hits.velocity = [1 0.5 3];"));
    const d1 = simulate({ ...sys, dt: 1 / 48000 / 8 }).energyDrift;
    const d2 = simulate({ ...sys, dt: 1 / 48000 / 16 }).energyDrift;
    expect(d1).toBeLessThan(1e-3);
    expect(d1 / d2).toBeGreaterThan(3);
    expect(d1 / d2).toBeLessThan(5);
  });

  it("gives a flat net the same transverse modes in the fast and the full analysis", () => {
    const { sys } = build(DOME("", 0));
    const fz = computeModes(sys, numblEig, "transverse").freqs;
    const full = Array.from(computeModes(sys, numblEig, "full").freqs);
    for (const f of fz) expect(Math.min(...full.map(g => Math.abs(g / f - 1)))).toBeLessThan(1e-8);
  });

  it("gives every 3D mode of a dome the bp_decay law", () => {
    const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
    const modes = computeModes(sys, numblEig, "full");
    const sigma = modalDamping(sys, modes);
    const s = [Math.log(1000) / 0.4, Math.log(1000) / 0.1];
    const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 3000) ** 2];
    const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
    const sigma0 = s[0] - sigma1 * w2[0];
    modes.freqs.forEach((f, k) => expect(sigma[k] / (sigma0 + sigma1 * (2 * Math.PI * f) ** 2)).toBeCloseTo(1, 8));
  });

  it.each([6, -3])("keeps a dome's energy in one 3D mode (%i), decaying at the predicted rate", kSel => {
    const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
    const modes = computeModes(sys, numblEig, "full");
    const sigma = modalDamping(sys, modes);
    const k = kSel >= 0 ? kSel : modes.freqs.length + kSel; // -3: one of the highest, stretching modes
    const nd = modes.dof.length;
    // initial velocity in the shape of mode k: v = eps * psi_k / sqrt(m)
    const v = new Float64Array(3 * sys.n);
    for (let r = 0; r < nd; r++) v[modes.dof[r]] = (1e-4 * modes.psi[k * nd + r]) / Math.sqrt(sys.mass[(modes.dof[r] / 3) | 0]);
    const hits = Array.from({ length: sys.n }, (_, b) => ({
      t: 0,
      beads: Int32Array.of(b),
      velocity: [v[3 * b], v[3 * b + 1], v[3 * b + 2]] as [number, number, number],
      force: [0, 0, 0] as [number, number, number],
      duration: 0,
    })).filter(h => !sys.pinned[h.beads[0]]);
    const r = simulate({ ...sys, hits, duration: 0.1, frameRate: 1000, dt: 1 / 48000 / 8 }, { modes });
    const nm = modes.freqs.length;
    const at = (t: number) => r.frameTimes.findIndex(x => x >= t);
    const share = (f: number) => r.modeEnergy![f * nm + k] / r.modeEnergy!.slice(f * nm, f * nm + nm).reduce((p, q) => p + q, 0);
    expect(share(at(0.05))).toBeGreaterThan(0.999);
    const e = (f: number) => r.modeEnergy![f * nm + k];
    const rate = Math.log(e(at(0.02)) / e(at(0.08))) / (r.frameTimes[at(0.08)] - r.frameTimes[at(0.02)]) / 2;
    expect(Math.abs(rate / sigma[k] - 1)).toBeLessThan(0.02);
  });
});