Answer the referee report and clarify the paper: notation, an example, a conclusion and the organism's records

paper/README.md

@@ -13,7 +13,13 @@ runs/                                 the record written by each command
 
 Each command is `python src/<file>.py <command>` and writes its record to `runs/`;
 `paper.py rocq` needs Rocq 9.1 with the opam packages `coq-coquelicot` and
-`coq-interval`, and ocamlfind, on the PATH.
+`coq-interval`, and ocamlfind, on the PATH; the development was checked on OCaml 4.14.2 with
+these versions:
+
+```bash
+opam install rocq-core.9.1.1 rocq-stdlib.9.2.0 coq-coquelicot.3.4.5 coq-interval.4.11.4 \
+  rocq-mathcomp-ssreflect.2.4.0 coq-flocq.4.2.2 rocq-bignums.9.0.0+rocq9.1 ocamlfind.1.9.8
+```
 
 | Command | Record | Establishes |
 |---|---|---|
@@ -41,9 +47,9 @@ Each command is `python src/<file>.py <command>` and writes its record to `runs/
 | `twohead.py organism` | `paper_twohead.json` | the organism on the ring machine M_{3,160,3200}, six generations with the copies and the period (`organism run --generations 6 --S 160 --L 3200`); `organism rocq` writes `rocq/Organism.v` |
 | `twohead.py net` | `paper_twohead_net.json` | the netlist N_{K,S,L} against the machine on random states of three members (`net check`), one generation of the organism on it in C beside the machine (`net lockstep`), the netlists extracted from Rocq (`net extract`), and the rows of its levelization with the noise a generation of the stored host survives (`net lev`) |
 | `twohead.py stored` | `paper_stored_host.json` | the host compiled into records of the ring machine: one pass against the host's step and its device on random states (`stored check`), the startup, step, device and hand-over records on decision diagrams (`stored symbolic`), and six generations in C in a ring of three banks (`stored run 6`); `--width A` checks H_A in place of H_8, and `stored widths` tabulates A = 8 to 12 |
-| `noise.py host` | `paper_noise_*.json` | the margin, the profile bound in sum and product form, read noise, the deviation curve, the exact Poisson intervals of the measured rates (`host intervals`), weight error with leakage, the active-input counts of the orbit certificate |
+| `noise.py host` | `paper_noise_*.json` | the margin, the profile bound in sum and product form, read noise, the deviation curve, the exact Poisson intervals of the measured rates (`host intervals`), weight error with leakage, the active-input counts of the orbit certificate, and the steps in which some row fails against the product form b' with the fraction the state masks (`host masked`) |
 | `noise.py amplified` | `paper_redundant_*.json` | the amplified host N_3: one step on replicated states, the margin 3/2, exact self-reproduction in C, the run under noise, the deviation curve and its bound by layer, the failure rate of the state read by majority and the restoration bound, computed by `src/restore_bound.c` |
-| `independent.py` | `independent_*.json` | code written from the paper's definitions alone: H_A rebuilt from Construction 5.8 for A = 3 to 12, Lev(N_host) from sigma and the measurements of Section 8.1 and Appendix C with the noise simulator `src/nsim.c` |
+| `independent.py` | `independent_*.json` | code written from the paper's definitions alone: H_A rebuilt from Construction 5.8 for A = 3 to 12, Lev(N_host) from sigma and the measurements of Section 8.1 and Appendix E with the noise simulator `src/nsim.c` |
 
 Rebuild the typeset paper with any TeX installation:
 

paper/runs/paper_noise_masked.json

@@ -0,0 +1,126 @@
+{
+ "M": 41967,
+ "grid": [
+  0.08,
+  0.09,
+  0.1,
+  0.11,
+  0.12,
+  0.15,
+  0.2
+ ],
+ "rows": [
+  {
+   "sigma": 0.08,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 31,
+   "deviations": 30,
+   "masked": 1,
+   "masked_fraction": 0.03225806451612903,
+   "row_failure_rate": 7.568359375e-06,
+   "deviation_rate": 7.32421875e-06,
+   "product_form_mean": 6.339578158376758e-06,
+   "row_failures_expected": 25.9669121367112,
+   "row_failures_sd": 5.09575779301011,
+   "seconds": 208.27974448201712
+  },
+  {
+   "sigma": 0.09,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 1739,
+   "deviations": 1639,
+   "masked": 100,
+   "masked_fraction": 0.05750431282346176,
+   "row_failure_rate": 0.000424560546875,
+   "deviation_rate": 0.000400146484375,
+   "product_form_mean": 0.00042732896942724436,
+   "row_failures_expected": 1750.3394587739929,
+   "row_failures_sd": 41.82811662962366,
+   "seconds": 212.16002563602524
+  },
+  {
+   "sigma": 0.1,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 35940,
+   "deviations": 33671,
+   "masked": 2269,
+   "masked_fraction": 0.06313299944351697,
+   "row_failure_rate": 0.0087744140625,
+   "deviation_rate": 0.008220458984375,
+   "product_form_mean": 0.008815790265558641,
+   "row_failures_expected": 36109.476927728196,
+   "row_failures_sd": 189.185309214523,
+   "seconds": 300.7043810179457
+  },
+  {
+   "sigma": 0.11,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 333003,
+   "deviations": 313064,
+   "masked": 19939,
+   "masked_fraction": 0.05987633745041336,
+   "row_failure_rate": 0.081299560546875,
+   "deviation_rate": 0.076431640625,
+   "product_form_mean": 0.08118096356967355,
+   "row_failures_expected": 332517.22678138287,
+   "row_failures_sd": 552.7370818611178,
+   "seconds": 312.9539852250018
+  },
+  {
+   "sigma": 0.12,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 1556080,
+   "deviations": 1479361,
+   "masked": 76719,
+   "masked_fraction": 0.04930273507788803,
+   "row_failure_rate": 0.37990234375,
+   "deviation_rate": 0.361172119140625,
+   "product_form_mean": 0.3795909822472992,
+   "row_failures_expected": 1554804.6632849376,
+   "row_failures_sd": 982.1122694402926,
+   "seconds": 279.61808501696214
+  },
+  {
+   "sigma": 0.15,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 4095990,
+   "deviations": 4095979,
+   "masked": 11,
+   "masked_fraction": 2.6855534315269323e-06,
+   "row_failure_rate": 0.99999755859375,
+   "deviation_rate": 0.999994873046875,
+   "product_form_mean": 0.9999982229194877,
+   "row_failures_expected": 4095992.7210782217,
+   "row_failures_sd": 2.6979452350816384,
+   "seconds": 226.470565449039
+  },
+  {
+   "sigma": 0.2,
+   "copies": 1024,
+   "steps": 4000,
+   "copy_steps": 4096000,
+   "row_failure_steps": 4096000,
+   "deviations": 4096000,
+   "masked": 0,
+   "masked_fraction": 0.0,
+   "row_failure_rate": 1.0,
+   "deviation_rate": 1.0,
+   "product_form_mean": 1.0,
+   "row_failures_expected": 4096000.0,
+   "row_failures_sd": 0.0,
+   "seconds": 209.7733603650122
+  }
+ ]
+}

paper/universal_construction_ternary.pdf

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c812d6dcb630f1f6a521890a5bf09da1e46b29974b356a3d4f8040f335426c4d
-size 362405
+oid sha256:0889b461dfbb6eef038b940bbaea4141a3cd589b30c095aeed305175cb481164
+size 385800

src/noise.py

@@ -529,6 +529,73 @@ def run_certify(args) -> int:
     return 0
 
 
+def run_masked(args) -> int:
+    """The steps in which some row fails, against the product form b'(t), and the fraction of them
+    that leave the state on the orbit; every noisy copy starts each step from the exact state."""
+    import torch
+    from host import M_P, LevEvaluator, describe, read_host
+    sigma = read_host()
+    L = LevEvaluator(sigma, device=args.device, dense=True)
+    M = sum(int(W.shape[0]) for W in L.W)
+    r_H = describe(sigma)
+    B, T = args.batch, args.steps
+    grid = [float(x) for x in args.grid.split(",")]
+    print(f"[masked] M = {M:,}; {B} noisy copies restarted from the exact state at each of "
+          f"{T:,} steps", flush=True)
+    rows = []
+    for s in grid:
+        gen = torch.Generator(device=args.device).manual_seed(11 + int(1000 * s))
+        v = L._vec(0, M_P).unsqueeze(0).to(args.device)
+        D = Devices(L, [r_H], args.device)
+        failed = deviated = 0
+        expected = variance = 0.0
+        t0 = time.perf_counter()
+        steps = T
+        for n in range(T):
+            y = v.repeat(B + 1, 1)
+            bad = torch.zeros(B, dtype=torch.bool, device=args.device)
+            keep = torch.zeros((), dtype=torch.float64, device=args.device)
+            for W, b in zip(L.W, L.B):
+                pre = y @ W.T + b
+                # copy 0 is noise free and reads the orbit's values, so its pre-activations are the
+                # orbit's, and a row fails with probability Q(|z + 1/2| / s)
+                q = 0.5 * torch.erfc((pre[0].double() + 0.5).abs() / (s * math.sqrt(2.0)))
+                keep += torch.log1p(-q).sum()
+                noise = torch.randn(pre.shape, generator=gen, device=pre.device) * s
+                noise[0] = 0.0
+                y = ((pre + noise) >= -0.5).float()
+                bad |= (y[1:] != y[0]).any(dim=1)
+            off = (y[1:] != y[0]).any(dim=1)
+            failed += int(bad.sum())
+            deviated += int(off.sum())
+            bp = float(-torch.expm1(keep))          # b'(t): some row of the step fails
+            expected += B * bp
+            variance += B * bp * (1.0 - bp)
+            v = y[0:1]
+            if D.step(v)[0]:
+                steps = n + 1
+                break
+        dt = time.perf_counter() - t0
+        n_cs = B * steps
+        rows.append({"sigma": s, "copies": B, "steps": steps, "copy_steps": n_cs,
+                     "row_failure_steps": failed, "deviations": deviated,
+                     "masked": failed - deviated,
+                     "masked_fraction": (failed - deviated) / failed if failed else None,
+                     "row_failure_rate": failed / n_cs, "deviation_rate": deviated / n_cs,
+                     "product_form_mean": expected / n_cs,
+                     "row_failures_expected": expected,
+                     "row_failures_sd": math.sqrt(variance),
+                     "seconds": dt})
+        print(f"    s={s:.2f}: {failed:,} steps with a failing row against {expected:,.1f} "
+              f"(sd {math.sqrt(variance):.1f}) from b'; {deviated:,} left the orbit, "
+              f"{failed - deviated:,} masked"
+              + (f" ({(failed - deviated) / failed:.2%})" if failed else "")
+              + f" ({dt / 60:.1f} min)", flush=True)
+    json.dump({"M": M, "grid": grid, "rows": rows},
+              open(runs_path("paper_noise_masked.json"), "w"), indent=1)
+    return 0
+
+
 def run_intervals(args) -> int:
     """Exact Poisson intervals for the measured rates of the host and of the threefold host."""
     out = {"confidence": 0.95, "host": [], "amplified_r3": []}
@@ -549,7 +616,8 @@ def run_intervals(args) -> int:
 
 def main_host() -> int:
     ap = argparse.ArgumentParser()
-    ap.add_argument("mode", choices=["omega", "curve", "leak", "margin", "profile", "certify", "intervals"])
+    ap.add_argument("mode", choices=["omega", "curve", "leak", "margin", "profile", "certify", "intervals",
+                                     "masked"])
     ap.add_argument("--short", action="store_true", help="certify: the leak orbit only")
     ap.add_argument("--states", type=int, default=800)
     ap.add_argument("--device", default="cuda")
@@ -559,7 +627,7 @@ def main_host() -> int:
     args = ap.parse_args()
     return {"omega": run_omega, "curve": run_curve, "leak": run_leak,
             "margin": run_margin, "profile": run_profile,
-            "certify": run_certify, "intervals": run_intervals}[args.mode](args)
+            "certify": run_certify, "intervals": run_intervals, "masked": run_masked}[args.mode](args)
 
 SRC = os.path.join(REPO, "src")