"""v4 — INERTIA. The grammar written by a body that cannot turn instantly. v3 hid a 16-state machine in the relations between turning events and proved a reader needs inference to open it. Then the machine's own specification arrived (diameter-rail simulator, 2026-07-28) and v3 turned out to be standing on a body that does not exist: disc radius v3 assumed 200 mm the machine is 400 mm carriage +/- 195 mm +/- 400 mm, contact inside 386.7 mm the tool could not lift lifts; cuts only while in contact a turn INSTANT, one step bounded by angular and carriage acceleration, a control period and a tolerated following error The last row is the whole problem. v3 wrote its entire message into the exit ANGLE at a junction — which is precisely where a real machine departs from the ideal, because a corner cannot be turned without lateral acceleration. So v4 does not assume the corner. It builds the axes. The rule commands a heading. The axes deliver what they can. What the metal receives is the difference. Concretely: a servo at a fixed control period resolves the commanded heading into a rotary rate and a carriage rate, slew-limits both by acceleration, clamps both by top speed, and integrates. The realized heading is whatever the achieved pair gives. Nothing about the corner is authored — its shape falls out of v^2/a, and that is a LENGTH THE BODY OWNS, not one I wrote. The series has been hunting for such a length since v2 (where the only length in the pattern turned out to be my step size). Inertia hands one over for nothing. The vow survives as a vow. The tool CAN lift on this machine; it does not. That is now a rule of the work rather than a property of the mechanism, and it is stated rather than assumed. What is measured, with the instruments already built: * does the grammar survive its own body — the v3 reading gradient, recomputed on a stroke whose corners are rounded by inertia; * how big the rounding is, and whether the corner radius is the metal's or mine (control: change the control period, change the acceleration); * the backbone invariants (does it still not stop, does the disc still fill). Usage: python3 v4.py --self-test python3 v4.py --hours 8 --out ../plate/v4/a1 python3 v4.py --decode ../plate/v4/a1 python3 v4.py --export ../plate/v4/a1 --name v4-8h # for the simulator """ from __future__ import annotations import argparse import json import math import os import time from array import array import numpy as np import grammar as G import yieldbreak as YB # --- the machine, from the specification ----------------------------------- # Lengths in millimetres, angles in radians, time in seconds. DISC_R = 400.0 # spec: radius 400 mm (diameter 800) RAIL_LO, RAIL_HI = -400.0, 400.0 CONTACT_R = 386.7 # spec impl: contact requires |r| < 2.9 units # The specification's disc has no centre hole, and the diameter rail passes # through the middle, so the reachable metal is the whole disc inside the # contact limit. The rail can travel to 400 mm; the METAL stops at 386.7. The # tool turns back at the metal's edge, not at the rail's — the vow (it never # lifts) makes the workpiece the boundary, not the mechanism. The first version # reflected at the rail and cut to 2.987 scene units, outside contact. R_MIN = 0.0 SURFACE_MM_S = 25.0 # commanded speed along the metal, unchanged CTRL_HZ = 500.0 # control period 2 ms. A declared number; --ctrl-hz DT = 1.0 / CTRL_HZ # controls it, and the control below proves the # pattern does not scale with it. # Rotary axis. The top speed is v1.1's 60 rpm. The acceleration is the # simulator's own default motion mode (2.35 rad/s^2 in its scene units, which # are dimensionless for an angular quantity, so it carries over unchanged). ROT_W_MAX = 60.0 * 2.0 * math.pi / 60.0 # 6.2832 rad/s ROT_A_MAX = 2.35 # rad/s^2 # Carriage. Top speed must exceed the surface speed or the machine could never # cut a radial line; the acceleration is the pair of the rotary figure, taken at # the rail's own scale (2.35 rad/s^2 x 133.33 mm per scene unit). LIN_V_MAX = 120.0 # mm/s LIN_A_MAX = 2.35 * 133.333 # 313.3 mm/s^2 STEP_MM = 0.6 # how far the tool travels between rule consultations SETTLE_DEG = 2.0 # the specification's "tolerated following error": the # controller does not issue a new corner while the # axes are still executing the last one. Without this # interlock the grammar commands a turn it has not # finished obeying, corners compound, and about 10% of # junctions land on the wrong alphabet angle. CELL_F = 0.25 # the metal's own record, as in v2/v3 CELL_C = 0.5 # the coverage grid, as in v1/v2 UNIT_MM = 400.0 / 3.0 # simulator scene units: 1 unit = 133.333 mm class Field: """The metal's record on a 400 mm disc: coverage and axial direction.""" def __init__(self, cell_f=CELL_F, cell_c=CELL_C): self.cf = cell_f self.nf = int(round(2.0 * DISC_R / cell_f)) self.n = array("I", bytes(4 * self.nf * self.nf)) self.c = array("f", bytes(4 * self.nf * self.nf)) self.s = array("f", bytes(4 * self.nf * self.nf)) self.cc = cell_c self.nc = int(round(2.0 * DISC_R / cell_c)) self.g = bytearray(self.nc * self.nc) self.ok = bytearray(self.nc * self.nc) self.touched = 0 n = 0 for iy in range(self.nc): y = (iy + 0.5) * cell_c - DISC_R for ix in range(self.nc): x = (ix + 0.5) * cell_c - DISC_R d2 = x * x + y * y if d2 <= CONTACT_R * CONTACT_R: self.ok[iy * self.nc + ix] = 1 n += 1 self.reachable = n def idx(self, x, y): ix = int((x + DISC_R) / self.cf) iy = int((y + DISC_R) / self.cf) if 0 <= ix < self.nf and 0 <= iy < self.nf: return iy * self.nf + ix return -1 def read(self, x, y): i = self.idx(x, y) if i < 0 or self.n[i] == 0: return False, 0.0, 0.0 c, s, n = self.c[i], self.s[i], self.n[i] return True, 0.5 * math.atan2(s, c), math.hypot(c, s) / n def cover(self, x, y): cx = int((x + DISC_R) / self.cc) cy = int((y + DISC_R) / self.cc) if 0 <= cx < self.nc and 0 <= cy < self.nc: i = cy * self.nc + cx if self.g[i] == 0 and self.ok[i]: self.touched += 1 if self.g[i] < 255: self.g[i] += 1 def coverage(self): return self.touched / max(1, self.reachable) def cross_and_cut(self, x0, y0, x1, y1, psi, rule): """As v2: walk the cells, read before writing, count crossings.""" seg = math.hypot(x1 - x0, y1 - y0) if seg <= 1e-12: return 0 phi = math.atan2(y1 - y0, x1 - x0) c2, s2 = math.cos(2.0 * phi), math.sin(2.0 * phi) k = max(1, int(seg / (self.cf * 0.5)) + 1) events = 0 last = -1 for j in range(k + 1): u = j / k xx, yy = x0 + (x1 - x0) * u, y0 + (y1 - y0) * u self.cover(xx, yy) i = self.idx(xx, yy) if i < 0 or i == last: continue last = i marked = self.n[i] != 0 if marked: phi_c = 0.5 * math.atan2(self.s[i], self.c[i]) across = YB.axial_diff(psi, phi_c) > YB.CHI_C else: across = False if across and not rule.crossing: events += 1 d = (phi_c - psi + math.pi * 0.5) % math.pi - math.pi * 0.5 rule.last_side = 1 if d > 0 else 0 rule.crossing = across self.n[i] += 1 self.c[i] += c2 self.s[i] += s2 return events class Axes: """The two axes, with the limits the specification says they have. The rule never touches r or theta. It hands down a heading; this object decides what actually happens, and the difference between the two is the only thing the metal ever sees. """ __slots__ = ("r", "th", "w", "rdot", "psi_cmd", "cut_mm", "ticks", "sat_w", "sat_a", "sat_v", "sat_ra", "err_sum", "err_max", "rail_hits", "speed_sum") def __init__(self, r, th, psi): self.r, self.th = r, th self.psi_cmd = psi # start already moving along the commanded heading, so the opening # transient is not mistaken for a corner self.rdot = SURFACE_MM_S * math.cos(psi - th) self.w = SURFACE_MM_S * math.sin(psi - th) / (r if abs(r) > 1e-9 else 1e-9) self.cut_mm = 0.0 self.ticks = 0 self.sat_w = self.sat_a = self.sat_v = self.sat_ra = 0 self.err_sum = 0.0 self.err_max = 0.0 self.rail_hits = 0 self.speed_sum = 0.0 def xy(self): return self.r * math.cos(self.th), self.r * math.sin(self.th) def heading(self): """The direction the tool is actually travelling, from the axis rates.""" vr, vt = self.rdot, self.r * self.w if abs(vr) < 1e-12 and abs(vt) < 1e-12: return self.psi_cmd return self.th + math.atan2(vt, vr) def tick(self): """One control period. Returns the segment cut, in mm.""" lever = self.r if abs(self.r) > 1e-9 else (1e-9 if self.r >= 0 else -1e-9) # what the commanded heading asks of each axis, at the commanded speed rdot_want = SURFACE_MM_S * math.cos(self.psi_cmd - self.th) w_want = SURFACE_MM_S * math.sin(self.psi_cmd - self.th) / lever if abs(w_want) > ROT_W_MAX: w_want = ROT_W_MAX if w_want > 0 else -ROT_W_MAX self.sat_w += 1 dw = w_want - self.w if abs(dw) > ROT_A_MAX * DT: dw = ROT_A_MAX * DT if dw > 0 else -ROT_A_MAX * DT self.sat_a += 1 self.w += dw if abs(rdot_want) > LIN_V_MAX: rdot_want = LIN_V_MAX if rdot_want > 0 else -LIN_V_MAX self.sat_v += 1 dv = rdot_want - self.rdot if abs(dv) > LIN_A_MAX * DT: dv = LIN_A_MAX * DT if dv > 0 else -LIN_A_MAX * DT self.sat_ra += 1 self.rdot += dv x0, y0 = self.xy() self.r += self.rdot * DT self.th += self.w * DT stopped = False if self.r < -CONTACT_R or self.r > CONTACT_R: self.r = -CONTACT_R if self.r < 0 else CONTACT_R self.rdot = -self.rdot # the rail's end stop self.rail_hits += 1 stopped = True x1, y1 = self.xy() seg = math.hypot(x1 - x0, y1 - y0) self.cut_mm += seg self.ticks += 1 self.speed_sum += seg / DT # following error: the angle between what was asked and what happened e = abs(YB.axial_diff(self.psi_cmd, self.heading())) self.err_sum += e if e > self.err_max: self.err_max = e return seg, x0, y0, x1, y1, stopped def run(hours=8.0, n_impulse=3, out_dir="../plate/v4/a1", log_every=200000, ctrl_hz=CTRL_HZ, rot_a=ROT_A_MAX, keep_stroke=True, settle_deg=SETTLE_DEG): global DT, CTRL_HZ, ROT_A_MAX CTRL_HZ, DT, ROT_A_MAX = ctrl_hz, 1.0 / ctrl_hz, rot_a settle_rad = math.radians(settle_deg) field = Field() rule = G.Rule3() ax = Axes(199.66, 0.0, math.radians(37.0)) total_ticks = int(hours * 3600.0 * ctrl_hz) pts = np.empty((total_ticks // 8 + 2, 2), dtype="float32") if keep_stroke \ else None np_i = 0 if keep_stroke: pts[0] = ax.xy() np_i = 1 junctions = [] junction_at = [] # stroke-sample index of each junction trace = [] travel = 0.0 t0 = time.time() for i in range(total_ticks): seg, x0, y0, x1, y1, stopped = ax.tick() events = field.cross_and_cut(x0, y0, x1, y1, ax.heading(), rule) if stopped: # The carriage is against its end stop and the metal ends here, so # the commanded heading has to come back inboard — otherwise the # servo re-accelerates into the stop every control period and the # machine chatters there (the first version did: 171244 hits in # 450000 ticks, mean speed down to 18.7 mm/s). Reflecting the # command is the mechanism refusing, not a rule of mine; and the # groove the tip was following has ended, so capture ends too. ax.psi_cmd = 2.0 * ax.th + math.pi - ax.psi_cmd if rule.mode == G.FOLLOW: ax.psi_cmd = rule.rail_exit(ax.psi_cmd) if rule.rail_sym is not None: junctions.append(rule.rail_sym) junction_at.append(np_i - 1) rule.rail_sym = None travel += seg if keep_stroke and (i & 7) == 0: pts[np_i] = (x1, y1) np_i += 1 settled = YB.axial_diff(ax.psi_cmd, ax.heading()) < settle_rad if travel >= STEP_MM and settled: travel = 0.0 # sense one record-cell ahead, as in v2/v3 px = x1 + field.cf * math.cos(ax.psi_cmd) py = y1 + field.cf * math.sin(ax.psi_cmd) marked, phi, _ = field.read(px, py) ax.psi_cmd, _cap = rule.steer(ax.psi_cmd, marked, phi) ax.psi_cmd, flipped, sym = rule.tick(ax.psi_cmd, events, n_impulse) if flipped is not None and sym is not None: junctions.append(sym) junction_at.append(np_i - 1) else: rule.count += events if (i + 1) % log_every == 0: trace.append({"tick": i + 1, "machine_h": round((i + 1) / ctrl_hz / 3600.0, 4), "coverage": round(field.coverage(), 5), "junctions": len(junctions), "r": round(ax.r, 1)}) wall = time.time() - t0 os.makedirs(out_dir, exist_ok=True) if keep_stroke: np.save(os.path.join(out_dir, "stroke.npy"), pts[:np_i]) with open(os.path.join(out_dir, "junctions.json"), "w") as f: json.dump(junctions, f) with open(os.path.join(out_dir, "junction_at.json"), "w") as f: json.dump(junction_at, f) r = np.hypot(pts[:np_i, 0], pts[:np_i, 1]) if keep_stroke else np.zeros(2) q = max(1, len(r) // 8) summary = { "rule": "grammar on real axes (v4)", "disc_radius_mm": DISC_R, "rail_mm": [RAIL_LO, RAIL_HI], "contact_limit_mm": CONTACT_R, "control_hz": ctrl_hz, "rot_a_max": rot_a, "settle_deg": settle_deg, "rot_w_max_rad_s": round(ROT_W_MAX, 4), "lin_v_max_mm_s": LIN_V_MAX, "lin_a_max_mm_s2": round(LIN_A_MAX, 1), "hours": hours, "ticks": total_ticks, "cut_m": round(ax.cut_mm / 1000.0, 2), "mean_surface_mm_s": round(ax.speed_sum / max(1, ax.ticks), 3), "coverage": round(field.coverage(), 5), "junctions": len(junctions), "following_error_mean_deg": round(math.degrees(ax.err_sum / max(1, ax.ticks)), 3), "following_error_max_deg": round(math.degrees(ax.err_max), 2), "saturation_fraction": { "rotary_speed": round(ax.sat_w / max(1, ax.ticks), 5), "rotary_accel": round(ax.sat_a / max(1, ax.ticks), 5), "carriage_speed": round(ax.sat_v / max(1, ax.ticks), 5), "carriage_accel": round(ax.sat_ra / max(1, ax.ticks), 5)}, "rail_hits": ax.rail_hits, "radial_spread_ratio": round(float(r[-q:].std() / max(1e-9, r.std())), 3), "wall_s": round(wall, 1), } with open(os.path.join(out_dir, "summary.json"), "w") as f: json.dump(summary, f, ensure_ascii=False, indent=1) with open(os.path.join(out_dir, "trace.json"), "w") as f: json.dump(trace, f) print(json.dumps(summary, ensure_ascii=False, indent=1)) return summary def junctions_from_inertial_stroke(pts, enter_deg=5.0, exit_deg=2.0, hold=3, min_run=4): """Read an inertial stroke by SEGMENTING it, not by peak-picking. v3's extractor looked for a single sample whose turn exceeded a threshold and then summed a short window. That is the right reading of a machine that turns instantly. This machine cannot: a 144 degree exit is spread over about 3.4 mm of travel, so the peak-picker fires inside the corner, splits it, and hands the reader a sequence with insertions — which costs more than plain substitution noise (measured: the peak-picker gave 1.588 bits on a calibration disc, where a flat 25% symbol-error rate costs only 1.455). So instead: a stroke is straight runs separated by CORNERS. Walk the per-sample turn with hysteresis, mark maximal corner regions, and take each region's NET turn as the symbol. The net turn is what inertia preserves — measured median error 0.0 degrees at the right window, with 90.6% of symbols landing on the correct alphabet angle even though every corner is rounded. Returns (symbols, sample indices). """ d = np.diff(pts.astype("float64"), axis=0) ang = np.arctan2(d[:, 1], d[:, 0]) turn = np.diff(ang) turn = (turn + math.pi) % (2.0 * math.pi) - math.pi n = len(turn) ent, ext = math.radians(enter_deg), math.radians(exit_deg) syms, idx = [], [] i = 0 last_end = -10 ** 9 while i < n: if abs(turn[i]) < ent: i += 1 continue if i - last_end < min_run: # too close to the last corner i += 1 continue acc = 0.0 j = i quiet = 0 while j < n: acc += turn[j] if abs(turn[j]) < ext: quiet += 1 if quiet >= hold: break else: quiet = 0 j += 1 deg = math.degrees(acc) best = min(range(len(ALPHABET_DEG)), key=lambda k: abs(deg - ALPHABET_DEG[k])) if abs(deg - ALPHABET_DEG[best]) < 36.0: # half the 72 deg spacing syms.append(best) idx.append(j) last_end = j i = j + 1 return syms, idx ALPHABET_DEG = [math.degrees(a) for a in G.ALPHABET] def corner_radius_mm(v=SURFACE_MM_S, r=200.0, rot_a=ROT_A_MAX): """What the body's inertia makes of a corner, analytically. Turning the velocity vector needs lateral acceleration. At radius r the rotary axis can supply r*alpha across the rail and the carriage can supply its own alpha_lin along it, so the available lateral acceleration is at least the smaller of the two, and the corner cannot be tighter than v^2 / a_lat. NOTHING HERE IS A NUMBER I CHOSE: v is the surface speed the work has always used, and the accelerations are the machine's. """ a_lat = min(abs(r) * rot_a, LIN_A_MAX) return v * v / max(1e-9, a_lat) def self_test(): ok = True def check(name, got, want, tol=0.0): nonlocal ok good = (abs(got - want) <= tol) if isinstance(want, float) \ else (got == want) ok = ok and good print(" %-58s %10s want %8s %s" % (name, ("%.4f" % got) if isinstance(got, float) else got, ("%.4f" % want) if isinstance(want, float) else want, "ok" if good else "FAIL")) print("the machine is the one in the specification") check("disc radius mm", DISC_R, 400.0) check("rail travel mm", RAIL_HI - RAIL_LO, 800.0) check("contact limit mm", CONTACT_R, 386.7) check("the tool turns back at the metal, not the rail", 1 if CONTACT_R < RAIL_HI else 0, 1) check("scene unit mm (for the simulator export)", UNIT_MM, 133.3333, 1e-3) print("the corner is the body's, not mine") # v^2/a_lat with a_lat = min(|r|*alpha_rot, alpha_lin). The two axes swap # which one binds at |r| = alpha_lin/alpha_rot = 313.3/2.35 = 133.3 mm: # OUTSIDE that the carriage binds (a flat 625/313.3 = 1.995 mm), INSIDE it # the lever arm shrinks and the rotary binds, so corners open up. check("crossover radius mm", LIN_A_MAX / ROT_A_MAX, 133.3333, 1e-3) check("corner radius at r=390 mm, mm", corner_radius_mm(r=390.0), 1.9947, 1e-3) check("corner radius at r=200 mm, mm", corner_radius_mm(r=200.0), 1.9947, 1e-3) check("corner radius at r=50 mm, mm", corner_radius_mm(r=50.0), 5.3191, 1e-3) check("outside 133 mm the carriage binds, so the corner is flat in r", 1 if abs(corner_radius_mm(r=390.0) - corner_radius_mm(r=140.0)) < 1e-9 else 0, 1) check("inside it the rotary binds and corners open up", 1 if corner_radius_mm(r=50.0) > 2.0 * corner_radius_mm(r=200.0) else 0, 1) check("a corner is 20x the groove width (0.10 mm)", 1 if corner_radius_mm(r=200.0) > 2.0 else 0, 0) check("a corner is bigger than the record cell (0.25 mm)", 1 if corner_radius_mm(r=200.0) > CELL_F else 0, 1) print("the axes obey their limits, and a straight line stays straight") a = Axes(200.0, 0.0, math.radians(37.0)) for _ in range(2000): a.tick() check("straight run: mean following error deg", math.degrees(a.err_sum / a.ticks), 0.0, 0.35) check("straight run: surface speed mm/s", a.speed_sum / a.ticks, SURFACE_MM_S, 0.6) check("no acceleration saturation on a straight run", a.sat_a, 0) print("and a 144 degree corner costs what inertia says it costs") b = Axes(200.0, 0.0, 0.0) for _ in range(400): b.tick() b.psi_cmd += math.radians(144.0) turned = 0 for _ in range(4000): b.tick() if YB.axial_diff(b.psi_cmd, b.heading()) < math.radians(2.0): break turned += 1 settle_mm = turned * SURFACE_MM_S * DT print(" corner settles in %.2f mm of travel (%d control ticks)" % (settle_mm, turned)) check("the corner is completed at all", 1 if turned < 3999 else 0, 1) check("and it costs between 0.5 and 20 mm of travel", 1 if 0.5 < settle_mm < 20.0 else 0, 1) print("\nSELF-TEST %s" % ("PASSED" if ok else "FAILED")) return ok def export_traces(plate_dir, name, out=None, max_pts=2400, stride=1): """Write the stroke as a simulator trace library entry. Format read from the simulator itself (localStorage key `nukaga-sim-trace-library-v1`): {name, createdAt, traces:[{material, points:[x,y,z,...]}]}, disc-local scene units at 1 unit = 400/3 mm, the tool plane at y = 0.078, values rounded to 1e-4. """ pts = np.load(os.path.join(plate_dir, "stroke.npy"))[::stride] traces = [] for c0 in range(0, len(pts), max_pts): part = pts[c0:c0 + max_pts] if len(part) < 2: continue flat = [] for x, y in part: flat.append(round(float(x) / UNIT_MM, 4)) flat.append(0.078) flat.append(round(float(y) / UNIT_MM, 4)) traces.append({"material": "scratch", "points": flat}) payload = {name: {"name": name, "createdAt": "2026-07-28T00:00:00.000Z", "traces": traces}} out = out or os.path.join(plate_dir, "sim_trace_library.json") with open(out, "w") as f: json.dump(payload, f) print("wrote %s (%d polylines, %d points, key nukaga-sim-trace-library-v1)" % (out, len(traces), sum(len(t["points"]) // 3 for t in traces))) return out def main(): ap = argparse.ArgumentParser() ap.add_argument("--self-test", action="store_true") ap.add_argument("--decode", default=None) ap.add_argument("--export", default=None) ap.add_argument("--name", default="v4") ap.add_argument("--stride", type=int, default=4) ap.add_argument("--hours", type=float, default=8.0) ap.add_argument("--n", type=int, default=3) ap.add_argument("--ctrl-hz", type=float, default=CTRL_HZ) ap.add_argument("--rot-a", type=float, default=ROT_A_MAX) ap.add_argument("--settle-deg", type=float, default=SETTLE_DEG) ap.add_argument("--out", default="../plate/v4/a1") args = ap.parse_args() if args.self_test: raise SystemExit(0 if self_test() else 1) if args.decode: G.decode(os.path.abspath(args.decode)) return if args.export: export_traces(os.path.abspath(args.export), args.name, stride=args.stride) return run(hours=args.hours, n_impulse=args.n, out_dir=args.out, ctrl_hz=args.ctrl_hz, rot_a=args.rot_a, settle_deg=args.settle_deg) if __name__ == "__main__": main()