Audit the whole library: 37 verified defects found and fixed

.claude/skills/comnumpy-tutorial/SKILL.md

@@ -95,6 +95,13 @@ then learns the wrong lesson -- they came to see how the library is used.
   loop filled -- each inner dict is exactly the `curves` that
   `print_data` and `plot_data` render. The loop is the pedagogy: it
   stays visible.
+- **One vocabulary across the pages.** The swept dict a page displays
+  is `curves` -- never `measured`, `collected` or `results` -- and its
+  closed-form companion is `theory`; a domain name (`snr_dB`,
+  `runtime`, `ccdf_curves`) is allowed only when it says strictly more.
+  The sweep axis carries its unit: `snr_dB_list`, `ebn0_dB_list`,
+  `dBm_list`. `reference` is reserved for the library's meaning -- the
+  known signal a data-aided block compares against.
 - **A swept result is one dictionary, shown two ways.**
   `data = {"x": snr_dB, "curves": {"ZF": ..., "ML": ...}}` — which is
   already the shape `monte_carlo` returns — then `print_data(data,

CHANGELOG.md

@@ -45,6 +45,84 @@ one release; there is no compatibility layer.
 | `core.metrics.calculate_acpr` | `compute_acpr` — it was the only `calculate_*` in the library, against 17 `compute_*` |
 | `core.metrics.compute_effective_SNR`, `ofdm.metrics.compute_PAPR` | `compute_effective_snr`, `compute_papr` — the two capitalized outliers among functions otherwise all lowercase (`compute_ser`, `compute_ber`, `compute_evm`, `compute_ccdf`, `compute_mi`) |
 
+### Fixed — a full-library audit: 37 findings, 37 fixed
+
+Five independent reviews swept the whole library -- core signal path,
+compensators/metrics/utils, optical, MIMO+FEC, and the chain machinery
+with OFDM -- hunting for wrong axes, dtype loss, RNG misuse, numerical
+blow-ups and validation gaps, each finding verified by running a minimal
+reproduction before being counted. The numerical cores came back clean
+against independent references (closed-form SER curves cross-checked by
+Monte-Carlo and numeric fading integration; Viterbi against exhaustive
+ML; LDPC min-sum against brute force; ML against the sphere decoder
+bit-for-bit; CD/DBP round trips to 1e-14; EDFA and Raman against their
+closed forms; the (133,171) distance spectrum against the published
+tables). The confirmed defects, all fixed and locked by
+`tests/test_audit_regressions.py`:
+
+- **Real mis-computations.** `Sequential.seed()` never reached blocks
+  inside a nested `Sequential` or a wrapper chain, so a seeded composite
+  chain was not reproducible; it now recurses. The AMP detector
+  evaluated its posterior-variance term `G` at variance `tau^2` instead
+  of `sigma^2 (1 + tau^2)` (Jeon et al., Algorithm 2) and its state went
+  NaN at high SNR -- silently hard-projected to symbol 0;
+  `soft_projector` also gained the log-sum-exp guard, which fixes the
+  same divergence in OAMP. `ErrorCounter` raveled before truncating, so
+  a batched reference against a delay-truncated detection counted
+  phantom errors on every row after the first. `CarrierAllocator`
+  allocated its buffer with the data dtype and silently zeroed complex
+  pilots on real data. `OFDMReceiver` read `N_fft` as `len(mask)`, which
+  is the *period* of a 2-D scattered mask. `compute_acpr` counted a bin
+  sitting exactly on a band edge in both bands. `SRRCFilter(rho=0)` -- 
+  the documented sinc case -- crashed on `1/(4 rho)`. A logarithmic
+  split-step schedule with `alpha_dB=0` (the default) was all-NaN; a
+  lossless fiber now gets the uniform spacing it is the limit of. The
+  FIR CD compensator applied a `sqrt(1/(4 pi K))` gain when its design
+  collapsed to one tap (short fiber, `z=0` crashed outright); below one
+  tap of memory it is now the identity it tends to. `Sequential.summary`
+  ran the chain without feeding `wiring` or recording taps.
+- **State that lied.** `FrameField` froze the *caller's* array in place;
+  it now freezes its own copy. `Deframer` recorded views of the input
+  buffer; the recorded fields now own their data. `tap()` served signals
+  from a previous run or a previous `taps=` configuration; the record is
+  cleared at each pass. `DataExtractor.forward` overwrote its declared
+  tuple selector with a slice. `Downsampler`/`Upsampler(use_filter=True)`
+  silently turned a real signal complex through `BWFilter`, which now
+  keeps a real signal real (its mask is symmetric).
+- **Contract mismatches.** A nested chain exported to JSON but failed to
+  import (`Sequential` was missing from the registry); int-keyed
+  callbacks were documented but never fired; the encoder/decoder pair
+  broke exactly at the empty message (`ViterbiDecoder` now accepts the
+  tail-only codeword); the distribution matcher pair died on a
+  zero-length stream; the Optical 90 hybrid and the MZM drew their local
+  oscillator with `len(x)` and crashed on any `(2, N)` input despite
+  advertising shape-agnostic behavior; `print_data` died on empty or
+  curveless tables instead of rendering the header.
+- **Validation gaps** (wrong or late errors): `compute_ser_awgn_qam`
+  accepted cross constellations (32, 128) and returned a wrong closed
+  form; negative `AWGN` variance produced NaN with only a warning;
+  `Serial2Parallel` validated `method=` only when the length did not
+  divide; `Complex2Real` leaked `UnboundLocalError`; `Upsampler(phase)`,
+  `SampleRemover` windows, `CyclicPrefixer(N_cp > N)`,
+  `Parallel2Serial` on 1-D input and `SRRCFilter(method='fft')` on a
+  signal shorter than the filter all died deep in numpy instead of
+  saying what was wrong; the generators rejected `np.int64` sizes and
+  accepted `True`; OAMP's `get_vt_2` ordered a complex trace.
+
+Also tightened along the way: `compute_ser` on an empty comparison now
+raises a clear `ValueError` instead of a bare `ZeroDivisionError`.
+
+### Changed — one vocabulary across the tutorial scripts
+
+The displayed sweep dict is `curves` everywhere, its closed-form
+companion `theory`, and the axis carries its unit (`snr_dB_list`);
+`measured`, `collected` and bare `reference` (outside the library's
+data-aided meaning) are gone, with domain names (`snr_dB`,
+`ccdf_curves`) kept only where they say strictly more. The canon is in
+the tutorial skill. The optical pages were already on the structural
+doctrine: `one_shot_pdm` inline, and the two DBP scripts keeping the
+one sanctioned chain split and their per-span progress table.
+
 ### Changed — every tutorial script held to the writing doctrine
 
 The conventions sweep closed over the whole tutorial set. The coding

examples/ofdm/monte_carlo_ofdm_papr.py

@@ -120,9 +120,9 @@
 ofdm_1024.seed(1)
 papr_1024 = compute_papr(ofdm_1024(n_symbols * N_sc), unit="dB", axis=-1)
 
-measured = {"$N_{sc}$ = 256": compute_ccdf(papr_256),
+ccdf_curves = {"$N_{sc}$ = 256": compute_ccdf(papr_256),
             "$N_{sc}$ = 1024": compute_ccdf(papr_1024)}
-reference = {
+theory = {
     "$N_{sc}$ = 256": compute_papr_ccdf_theo(threshold_dB, 256,
                                              oversampling=os, unit="dB"),
     "$N_{sc}$ = 1024": compute_papr_ccdf_theo(threshold_dB, 1024,
@@ -134,11 +134,11 @@
 levels = np.logspace(0, -4, 25)
 
 fig, ax = plt.subplots(figsize=(7, 4.5))
-for name, (sorted_dB, ccdf) in measured.items():
+for name, (sorted_dB, ccdf) in ccdf_curves.items():
     shown = np.clip((ccdf.size * (1 - levels)).astype(int), 0, ccdf.size - 1)
     line, = ax.plot(sorted_dB[shown], ccdf[shown], "o", fillstyle="none",
                     label=f"{name}, simulated")
-    ax.plot(threshold_dB, reference[name], "-", color=line.get_color(),
+    ax.plot(threshold_dB, theory[name], "-", color=line.get_color(),
             label=f"{name}, theory")
 ax.set_yscale("log")
 ax.set_ylim(1e-4, 1)

examples/ofdm/one_shot_ofdm.py

@@ -106,25 +106,25 @@
 # pools the n_trials frames of each sweep point -- exactly the mean of
 # the per-trial rates. The channel stays the same throughout, by
 # construction.
-snr_list = np.arange(6, 22, 2)
+snr_dB_list = np.arange(6, 22, 2)
 n_trials = 2
 
 # --- the sweep: no simulation loop -----------------------------------
-measured = {}
-measured["single carrier"] = monte_carlo(
-    sc_chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser},
+curves = {}
+curves["single carrier"] = monte_carlo(
+    sc_chain, "data_rx.snr_dB", snr_dB_list, {"ser": compute_ser},
     (n_trials, N), reference="data_tx", seed=1)["ser"]
-measured["OFDM"] = monte_carlo(
-    ofdm_chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser},
+curves["OFDM"] = monte_carlo(
+    ofdm_chain, "data_rx.snr_dB", snr_dB_list, {"ser": compute_ser},
     (n_trials, N), reference="data_tx", seed=1)["ser"]
 
 # --- results: table and figure ---------------------------------------
 
-ser_data = {"x": snr_list, "curves": measured}
+ser_data = {"x": snr_dB_list, "curves": curves}
 print()
 print_data(ser_data, xlabel="SNR [dB]", ylabel="SER")
 
-plot_error_rate(snr_list, measured, ylabel="SER",
+plot_error_rate(snr_dB_list, curves, ylabel="SER",
                 title="16-QAM over one EPA realization")
 plt.savefig(f"{img_dir}/one_shot_ofdm_fig4.png")
 

examples/optical/one_shot_NLI.py

@@ -169,14 +169,14 @@ def score(chain):
 
 # the dicts are keyed by span count in sweep order, so their values
 # are already the curves
-measured = np.array(list(snr_per_span.values()))
+snr_measured = np.array(list(snr_per_span.values()))
 ase_only = np.array(list(snr_ase_only.values()))
 
 # An effective SNR is not an error rate: linear ordinate, two curves,
 # nothing a helper would say better than the four lines that draw it.
 fig3, ax = plt.subplots()
 ax.plot(spans, ase_only, "o-", fillstyle="none", label="amplifier noise only")
-ax.plot(spans, measured, "s-", fillstyle="none",
+ax.plot(spans, snr_measured, "s-", fillstyle="none",
         label="dispersion compensation only")
 ax.set_xlabel("spans travelled")
 ax.set_ylabel("effective SNR [dB]")

examples/simple/monte_carlo_awgn.py

@@ -95,8 +95,8 @@ def qfunc(x):
 N_compare = 100000
 orders = np.array([4, 16, 64, 256])
 bits_per_symbol = np.log2(orders).astype(int)     # vectorized, not a loop
-measured_ber = {}
-theory_ber = {}
+curves = {}
+theory = {}
 for index, order in enumerate(orders):
     other = Constellation("QAM", int(order))
     bits = int(bits_per_symbol[index])
@@ -110,10 +110,10 @@ def qfunc(x):
                             {"ber": partial(compute_ber, width=bits)},
                             N_compare, reference="tx", seed=1)
     order_label = f"{order}-QAM"
-    measured_ber[order_label] = collected["ber"]
-    theory_ber[order_label] = other.metrics(ebn0_dB_list, per="bit")["ber"]
+    curves[order_label] = collected["ber"]
+    theory[order_label] = other.metrics(ebn0_dB_list, per="bit")["ber"]
 
-ax = plot_error_rate(ebn0_dB_list, measured_ber, theory=theory_ber,
+ax = plot_error_rate(ebn0_dB_list, curves, theory=theory,
                      xlabel="$E_b/N_0$ [dB]", ylabel="BER",
                      title="Square QAM at equal energy per bit")
 # The closed forms run down to 1e-12; the measurements cannot follow them

src/comnumpy/core/channels.py

@@ -99,6 +99,13 @@ def noise_rvs(self, x: np.ndarray) -> None:
             P_x = float(np.mean(np.abs(x) ** 2))
             sigma2n = P_x * 10 ** (-self.snr_dB / 10)
 
+        if sigma2n < 0:
+            # a swept sigma2 gone negative would otherwise turn into
+            # sqrt(<0) and an all-NaN pass with only a RuntimeWarning
+            raise ValueError(
+                f"AWGN: the noise variance is negative (sigma2={sigma2n}); "
+                f"a variance is non-negative -- check the sweep values.")
+
         shape = x.shape
         if np.iscomplexobj(x):
             scale = np.sqrt(sigma2n / 2)

src/comnumpy/core/compensators.py

@@ -24,8 +24,13 @@
 ]
 
 
-@dataclass(slots=True)
 class DataAidedMixin():
+    # a plain class, not a dataclass: it has no fields of its own, and
+    # dataclass(slots=True) would *recreate* the class, which multiple
+    # inheritance with the equally-recreated Processor cannot express.
+    # The empty slots keep the concrete slotted subclasses dict-free.
+    __slots__ = ()
+
     if TYPE_CHECKING:
         # declared by every concrete class; annotating it here would
         # make it a field of the mixin and reorder theirs

src/comnumpy/core/filters.py

@@ -2,6 +2,7 @@
 from dataclasses import dataclass, field
 from typing import Literal, Optional
 from comnumpy.core.generics import Processor
+from comnumpy.exceptions import ShapeError
 
 __all__ = ["SRRCFilter", "BWFilter"]
 
@@ -144,7 +145,9 @@ def h(self, t: Optional[np.ndarray] = None) -> np.ndarray:
             if t_temp == 0:
                 h_temp = (1 + rho*((4/np.pi)-1))
             else:
-                if np.abs(t_temp) == (1/(4*rho)):
+                # rho = 0 is the sinc pulse: the singularity at 1/(4*rho)
+                # moves to infinity and the general branch reduces to it
+                if rho != 0 and np.abs(t_temp) == (1/(4*rho)):
                     term1 = (1 + (2/np.pi))
                     term2 = (1 - (2/np.pi))
                     coef = rho / np.sqrt(2)
@@ -169,6 +172,11 @@ def H(self, NFFT: int) -> np.ndarray:
         from comnumpy._backend import fft  # local import (D36), cupy-compatible (D3)
         # see hager code on LDBP
         h = self.h()
+        if NFFT < len(h):
+            raise ShapeError(
+                f"SRRCFilter(method='fft') needs at least as many samples as "
+                f"filter taps ({len(h)} = 2*N_h*oversampling + 1), got {NFFT} "
+                f"-- use method='lfilter' or a longer signal")
         filter_delay = self.oversampling*self.N_h
         H_tmp = np.concatenate((h, np.zeros(NFFT-len(h))))
         H_tmp = np.roll(H_tmp, -filter_delay)
@@ -295,4 +303,8 @@ def forward(self, x: np.ndarray) -> np.ndarray:
         H = (abs(w) <= self.wn / 2).astype(float)
         fft_x = fft(x, NFFT)
         y = ifft(H*fft_x, NFFT)
+        if np.isrealobj(x):
+            # the mask is symmetric in +/-f, so a real signal stays real:
+            # keep its dtype instead of leaking the ifft's complex128
+            y = np.real(y)
         return y

src/comnumpy/core/frames.py

@@ -67,7 +67,9 @@ def __post_init__(self) -> None:
         # the enum, so `role.name` works and an invalid role fails here
         object.__setattr__(self, "role", FieldRole(self.role))
         if self.values is not None:
-            values = np.asarray(self.values)
+            # np.array, not asarray: freezing must apply to a copy this
+            # field owns, never to the caller's own buffer
+            values = np.array(self.values)
             values.setflags(write=False)
             object.__setattr__(self, "values", values)
             if self.length is not None and self.length != len(values):
@@ -351,8 +353,10 @@ def get_field(self, name: str) -> np.ndarray:
         return self.received_fields[name]
 
     def forward(self, X: np.ndarray) -> np.ndarray:
+        # copies, not views: the recorded fields are state that outlives
+        # the pass, and must not change if the caller reuses the buffer
         self.received_fields = {
-            f.name: X[..., self.frame.slice_of(f.name)]
+            f.name: X[..., self.frame.slice_of(f.name)].copy()
             for f in self.frame.fields}
         payload_name = self.frame.fields_by_role(FieldRole.PAYLOAD)[0].name
         return self.received_fields[payload_name]

src/comnumpy/core/generators.py

@@ -113,8 +113,10 @@ def __post_init__(self):
         self.rng = np.random.default_rng(self.seed)
 
     def forward(self, X: object) -> np.ndarray:
-        if isinstance(X, int):
-            size = (X,)
+        # np.integer too: sizes come out of shape arithmetic. bool is an
+        # int subclass and would slip through as 0 or 1.
+        if isinstance(X, (int, np.integer)) and not isinstance(X, bool):
+            size = (int(X),)
         elif isinstance(X, (tuple, list)):
             size = tuple(X)
         else:
@@ -180,8 +182,10 @@ def __post_init__(self):
         self.rng = np.random.default_rng(self.seed)
 
     def forward(self, X: object) -> np.ndarray:
-        if isinstance(X, int):
-            size = (X,)
+        # np.integer too: sizes come out of shape arithmetic. bool is an
+        # int subclass and would slip through as 0 or 1.
+        if isinstance(X, (int, np.integer)) and not isinstance(X, bool):
+            size = (int(X),)
         elif isinstance(X, (tuple, list)):
             size = tuple(X)
         else: