e213c1ebeadnet: Python port of the bead-percussion solver (numba), Zarr v2 output with every bead recorded, web components (run, controls, model, animation, waveform, spectrogram, energy, info), serve CLI, agent guide, examples, tests against the TS solverJeremy Magland 1import json
2import os
4import numpy as np
5import pytest
7import beadnet as bn
8from beadnet.kernels import hinge_angle
10HERE = os.path.dirname(__file__)
13def string_model(N=15, a=0.01, m=1e-3, f0=10.0, k=1e4, **kw):
14 pos, E, _ = bn.chain(N + 2, (N + 1) * a)
15 M = bn.Model(pos=pos, mass=m, pinned=[0, N + 1], springs=bn.springs(pos, E, k, f0))
16 M.hits = [bn.Hit(beads=4, velocity=[0, 0, 1e-4])]
17 M.duration = 0.2
18 for key, val in kw.items():
19 setattr(M, key, val)
20 return M
23def string_freq(kk, N=15, a=0.01, m=1e-3, f0=10.0):
24 return 2 * np.sqrt(f0 / (m * a)) * np.sin(kk * np.pi / (2 * (N + 1))) / (2 * np.pi)
27def peak_frequency(x, fs, f_lo, f_hi):
28 """Direct DFT magnitude with a Hann window on a 0.5 Hz grid, refined by a parabola."""
29 i = np.arange(len(x))
30 w = 0.5 - 0.5 * np.cos(2 * np.pi * i / (len(x) - 1))
31 xw = w * x
33 def mag(f):
34 return abs(np.sum(xw * np.exp(-2j * np.pi * f * i / fs)))
36 grid = np.arange(f_lo, f_hi, 0.5)
37 mags = np.array([mag(f) for f in grid])
38 f = grid[np.argmax(mags)]
39 a, b, c = mag(f - 0.5), mag(f), mag(f + 0.5)
40 return f + 0.25 * (a - c) / (a - 2 * b + c)
43def test_string_modes_analytic():
44 md = bn.modes(string_model())
45 assert md.kind == "transverse"
46 for kk in range(1, 16):
47 assert md.freqs[kk - 1] == pytest.approx(string_freq(kk), rel=1e-9)
50def test_string_rings_at_fundamental(tmp_path):
51 M = string_model()
52 M.hits = [bn.Hit(beads=8, velocity=[0, 0, 1e-4])]
53 r = bn.run(M, str(tmp_path / "s"), progress=False)
54 assert r.meta["aborted"] is None
55 f = peak_frequency(r.signal(8, "displacement"), r.sample_rate, 0.8 * string_freq(1), 1.2 * string_freq(1))
56 assert abs(f - string_freq(1)) / string_freq(1) < 2e-3
59def test_equilibrium_and_energy_conservation(tmp_path):
60 M = string_model(friction=0.0, method="yoshida4")
61 M.hits = [bn.Hit(beads=4, velocity=[0.5, 0.3, 2.0])]
62 r = bn.run(M, str(tmp_path / "s"), progress=False)
63 assert r.meta["residual_force"] < 1e-9
64 assert r.meta["energy_drift"] < 1e-6
65 assert np.all(r.dissipated == 0)
68def test_hinge_gradient_finite_difference():
69 rng = np.random.default_rng(1)
70 x = rng.normal(size=12)
71 g = np.zeros(12)
72 th = hinge_angle(x, 0, 1, 2, 3, g, 0)
73 eps = 1e-7
74 for q in range(12):
75 xp = x.copy()
76 xp[q] += eps
77 xm = x.copy()
78 xm[q] -= eps
79 d = (hinge_angle(xp, 0, 1, 2, 3, np.zeros(12), 0) - hinge_angle(xm, 0, 1, 2, 3, np.zeros(12), 0)) / (2 * eps)
80 assert d == pytest.approx(g[q], rel=1e-5, abs=1e-7)
81 assert np.isfinite(th)
85def test_decay_law_exact_at_small_amplitude(curved):
86 pos, tri, rim, trip = bn.hexdisk(0.06, 0.012)
87 if curved:
88 r = np.linalg.norm(pos[:, :2], axis=1)
89 pos[:, 2] = 0.01 * (1 - (r / 0.06) ** 2)
90 M = bn.Model(pos=pos, mass=1e-3, pinned=bn.near(pos, [0, 0, pos[:, 2].max()], 0.013),
91 springs=bn.springs(pos, bn.edges(tri), 2e6, 0), hinges=bn.hinges(tri, 50.0))
92 else:
93 M = bn.Model(pos=pos, mass=1e-3, pinned=rim, springs=bn.springs(pos, bn.edges(tri), 2e5, 20.0),
94 bends=bn.bends(trip, 1e3))
95 bn.decay(M, [(200, 2.0), (3000, 0.5)])
96 md = bn.modes(M)
97 sigma = bn.modal_damping(M, md)
98 s0, s1 = M.decay_law
99 expected = s0 + s1 * (2 * np.pi * md.freqs) ** 2
100 assert np.allclose(sigma, expected, rtol=1e-8, atol=1e-8 * expected.max())
103def test_geometry_helpers():
104 pos, tri, rim, trip = bn.hexdisk(0.15, 0.01)
105 assert pos.shape == (823, 3)
106 E = bn.edges(tri)
107 S = bn.springs(pos, E, 1e4, 10.0)
108 M = bn.Model(pos=pos, mass=1e-3, pinned=rim, springs=S)
109 prep = bn.prepare(bn.check(M))
110 assert prep["residual_force"] < 1e-9
111 q = bn.hinge_quads(tri)
112 assert q.shape[1] == 4 and len(q) == len(E) - len(rim) # interior edges
113 p2, cords, anchors = bn.anchor(pos, [0, 5], 500.0, 0.02)
114 assert p2.shape[0] == pos.shape[0] + 8 and cords.shape == (8, 4)
117def test_store_roundtrip(tmp_path):
118 from beadnet.store import Store
119 st = Store(str(tmp_path / "z"), "w")
120 a = np.arange(7 * 10, dtype=np.float32).reshape(7, 10)
121 st.write("a", a, chunks=(3, 4))
122 assert np.array_equal(st.read("a"), a)
123 assert np.array_equal(st.read("a", (2, slice(3, 9))), a[2, 3:9])
124 assert np.array_equal(st.read("a", (slice(1, 6), slice(None))), a[1:6])
125 zarr = pytest.importorskip("zarr")
126 z = zarr.open(str(tmp_path / "z"), mode="r")
127 assert np.array_equal(z["a"][:], a)
130def test_run_store_and_model_roundtrip(tmp_path):
131 M = string_model(record=("acceleration", "velocity", "displacement"), modes=True, duration=0.05)
132 bn.decay(M, [(30, 1.0), (300, 0.3)])
133 r = bn.run(M, str(tmp_path / "s"), progress=False)
134 acc = np.stack([r.signal(b) for b in range(M.n)])
135 assert np.allclose(acc.sum(axis=0), r.signal("sum"), rtol=1e-5, atol=1e-6 * np.abs(acc).max())
136 assert np.allclose(r.signal([3, 4]), acc[3] + acc[4], rtol=1e-6)
137 assert r.frames().shape == (r.meta["n_frames"], M.n, 3)
138 assert r.mode_energy.shape == (r.meta["n_frames"], 15)
139 M2 = r.model
140 r2 = bn.run(M2, str(tmp_path / "s2"), progress=False)
141 assert np.array_equal(r2.signal(4), r.signal(4))
142 r.write_wav(str(tmp_path / "x.wav"), 4)
145def _from_fixture(f):
146 n = f["n"]
147 M = bn.Model(pos=np.array(f["pos0"]).reshape(n, 3), mass=np.array(f["mass"]),
148 pinned=np.flatnonzero(np.array(f["pinned"])), friction=np.array(f["friction"]))
149 M.springs = np.array(f["springs"]).T.reshape(-1, 5)
150 M.bends = np.array(f["bends"]).T.reshape(-1, 5)
151 M.hinges = np.array(f["hinges"]).T.reshape(-1, 6)
152 M.dashpots = np.array(f["dashpots"]).T.reshape(-1, 4)
153 M.hits = [bn.Hit(beads=h["beads"], t=h["t"], velocity=None if h["duration"] > 0 else h["velocity"],
154 force=h["force"] if h["duration"] > 0 else None, duration=h["duration"]) for h in f["hits"]]
155 M.duration = f["duration"]
156 M.sample_rate = f["sampleRate"]
157 M.dt = f["dtGiven"]
158 M.method = f["method"]
159 M.record = ("acceleration", "velocity", "displacement")
160 return M
163with open(os.path.join(HERE, "fixtures", "ts_reference.json")) as _f:
164 TS_REF = json.load(_f)
168def test_matches_typescript_solver(name, tmp_path):
169 """Same model, same step: the per-bead signals agree with the TypeScript solver's listener."""
170 f = TS_REF[name]
171 M = _from_fixture(f)
172 r = bn.run(M, str(tmp_path / name), progress=False)
173 assert r.meta["dt"] == pytest.approx(f["dt"], rel=1e-12)
174 for q in ("acceleration", "velocity", "displacement"):
175 ref = np.array(f["audio"][q])
176 got = r.signal(f["listen"], q)
177 err = np.max(np.abs(got - ref)) / np.max(np.abs(ref))
178 assert err < 1e-5, (q, err)
181def test_blocks_give_the_same_signals(tmp_path):
182 """Several blocks (the last one partial) record the same signals as one block."""
183 M = string_model(duration=0.05)
184 bn.decay(M, [(30, 1.0), (300, 0.3)])
185 r1 = bn.run(M, str(tmp_path / "one"), progress=False)
186 r2 = bn.run(M, str(tmp_path / "many"), progress=False, block=1000)
187 for q in ("acceleration", "displacement"):
188 assert np.array_equal(r1.signal(5, q), r2.signal(5, q))
189 assert np.array_equal(r1.signal("sum", q), r2.signal("sum", q))
190 assert np.array_equal(r1.frames(), r2.frames())