diff --git a/.claude/skills/comnumpy-tutorial/SKILL.md b/.claude/skills/comnumpy-tutorial/SKILL.md index 152f6bc..b5e998a 100644 --- a/.claude/skills/comnumpy-tutorial/SKILL.md +++ b/.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, diff --git a/CHANGELOG.md b/CHANGELOG.md index 03bd35d..3f660e2 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/examples/ofdm/monte_carlo_ofdm_papr.py b/examples/ofdm/monte_carlo_ofdm_papr.py index 17b0845..6e77261 100644 --- a/examples/ofdm/monte_carlo_ofdm_papr.py +++ b/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) diff --git a/examples/ofdm/one_shot_ofdm.py b/examples/ofdm/one_shot_ofdm.py index 348cd40..f9af515 100644 --- a/examples/ofdm/one_shot_ofdm.py +++ b/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") diff --git a/examples/optical/one_shot_NLI.py b/examples/optical/one_shot_NLI.py index 122e386..ebc2137 100644 --- a/examples/optical/one_shot_NLI.py +++ b/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]") diff --git a/examples/simple/monte_carlo_awgn.py b/examples/simple/monte_carlo_awgn.py index 89e8e1f..e894d23 100644 --- a/examples/simple/monte_carlo_awgn.py +++ b/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 diff --git a/src/comnumpy/core/channels.py b/src/comnumpy/core/channels.py index 562d594..c2a6f24 100644 --- a/src/comnumpy/core/channels.py +++ b/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) diff --git a/src/comnumpy/core/compensators.py b/src/comnumpy/core/compensators.py index e987df3..550b103 100644 --- a/src/comnumpy/core/compensators.py +++ b/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 diff --git a/src/comnumpy/core/filters.py b/src/comnumpy/core/filters.py index 44f7d52..f240c7e 100644 --- a/src/comnumpy/core/filters.py +++ b/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 diff --git a/src/comnumpy/core/frames.py b/src/comnumpy/core/frames.py index 0a8e505..e855092 100644 --- a/src/comnumpy/core/frames.py +++ b/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] diff --git a/src/comnumpy/core/generators.py b/src/comnumpy/core/generators.py index edceae4..52180c3 100644 --- a/src/comnumpy/core/generators.py +++ b/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: diff --git a/src/comnumpy/core/generics.py b/src/comnumpy/core/generics.py index a6728c6..8557fda 100644 --- a/src/comnumpy/core/generics.py +++ b/src/comnumpy/core/generics.py @@ -335,8 +335,10 @@ def seed(self, seed: int) -> "Sequential": A ``numpy.random.SeedSequence`` spawned from ``seed`` gives each block that declares a ``seed`` field its own independent child seed, then re-runs its parametric initialization so the RNG is - rebuilt. Same chain + same seed = same signal, whatever the - number or order of stochastic blocks. + rebuilt. A nested :class:`Sequential` -- in the module list or + held by a wrapper block -- is seeded recursively from its own + child seed. Same chain + same seed = same signal, whatever the + number, order or nesting of stochastic blocks. Examples -------- @@ -349,18 +351,30 @@ def seed(self, seed: int) -> "Sequential": True """ seed_sequence = np.random.SeedSequence(seed) - stochastic = [ - module for module in self.module_list - if dataclasses.is_dataclass(module) - and "seed" in {f.name for f in dataclasses.fields(module)} - ] - for module, child in zip(stochastic, seed_sequence.spawn(len(stochastic)), + targets: List[object] = [] + for module in self.module_list: + if isinstance(module, Sequential): + targets.append(module) + elif dataclasses.is_dataclass(module): + if "seed" in {f.name for f in dataclasses.fields(module)}: + targets.append(module) + else: + # a wrapper block holding an inner chain (e.g. the + # OFDM chains) is seeded through that chain + for f in dataclasses.fields(module): + value = getattr(module, f.name, None) + if isinstance(value, Sequential): + targets.append(value) + for target, child in zip(targets, seed_sequence.spawn(len(targets)), strict=True): - # setattr, not module.seed: the field was discovered by name + child_seed = int(child.generate_state(1)[0]) + if isinstance(target, Sequential): + target.seed(child_seed) + continue + # setattr, not target.seed: the field was discovered by name # above and no base class declares it (noqa: B010 is the price) - setattr(module, "seed", # noqa: B010 - int(child.generate_state(1)[0])) - post_init = getattr(module, "__post_init__", None) + setattr(target, "seed", child_seed) # noqa: B010 + post_init = getattr(target, "__post_init__", None) if post_init is not None: post_init() return self @@ -464,14 +478,26 @@ def summary(self, X: Any, print_out: bool = True) -> List[_SummaryRow]: """ Run the chain on ``X`` and tabulate each block's output shape, dtype and execution time (decision D33b). Returns the rows. + + The pass is the ordinary one: taps are recorded and the data + edges declared in ``wiring`` are fed, so the chain summarized is + the chain that runs. """ + ids, recorded, feeds = self._resolve_edges() + if ids is None: + ids = self.block_ids() + self.tapped_.clear() rows: List[_SummaryRow] = [] Y = X for index, (block_id, module) in enumerate( - zip(self.block_ids(), self.module_list, strict=True)): + zip(ids, self.module_list, strict=True)): + for param, source in feeds.get(index, ()): + setattr(module, param, self.tapped_[source]) start_time = time.time() Y = module(Y) elapsed_ms = (time.time() - start_time) * 1e3 + if block_id in recorded: + self.tapped_[block_id] = Y shape = getattr(Y, "shape", None) dtype = getattr(Y, "dtype", type(Y).__name__) rows.append((index, type(module).__name__, block_id, @@ -599,6 +625,7 @@ def profile_execution_time(self, X: np.ndarray) -> Dict[str, float]: ids, recorded, feeds = self._resolve_edges() if ids is None: ids = self.block_ids() + self.tapped_.clear() Y = X time_elapsed: Dict[str, float] = {} @@ -660,6 +687,9 @@ def forward(self, X: np.ndarray) -> np.ndarray: Process the input data through all modules in the sequence. """ ids, recorded, feeds = self._resolve_edges() + # a tap is a record of *this* pass: stale entries from an earlier + # run or an earlier taps= configuration must not survive it + self.tapped_.clear() Y = X start = time.perf_counter() @@ -681,6 +711,9 @@ def forward(self, X: np.ndarray) -> np.ndarray: key = getattr(processor, 'name', None) if key is not None and key in callbacks: callbacks[key](Y) + elif index in callbacks: + # the documented alternative: callbacks keyed by position + callbacks[index](Y) self.elapsed_ = time.perf_counter() - start return Y diff --git a/src/comnumpy/core/metrics.py b/src/comnumpy/core/metrics.py index 9455a05..e9f7427 100644 --- a/src/comnumpy/core/metrics.py +++ b/src/comnumpy/core/metrics.py @@ -183,6 +183,12 @@ def compute_ser_awgn_qam(order: int, snr_per_bit: np.ndarray | float from scipy.stats import norm # local import (D36) gamma_b = snr_per_bit + root = np.sqrt(order) + if order < 4 or root != int(root): + raise ValueError( + f"compute_ser_awgn_qam: expected a square QAM order (4, 16, 64, " + f"256, ...), got order={order} -- the two-independent-PAM " + f"decomposition below does not hold for cross constellations.") # see book Proakis "Digital communication", p 280 M = order @@ -648,6 +654,10 @@ def compute_ser(X_target: np.ndarray, X_detected: np.ndarray, # of a batch after the first -- a silent offset, not an error N = min(x_target.shape[-1], x_detected.shape[-1]) difference = x_target[..., :N] - x_detected[..., :N] + if difference.size == 0: + raise ValueError( + "compute_ser: no symbols to compare (an input is empty, or " + "the two share no overlap along the last axis)") return np.count_nonzero(difference) / difference.size else: N = X_target.shape[axis] @@ -1123,10 +1133,15 @@ def compute_acpr(signal: np.ndarray, bandwidth: float, [-27.6193 -27.6193] """ - def band_power(lower_freq: float, upper_freq: float) -> float: + def band_power(lower_freq: float, upper_freq: float, *, + open_lower: bool = False, open_upper: bool = False) -> float: # Parseval: the power in a band is the energy of the bins it - # contains, divided by the number of samples - mask = (freq_axis >= lower_freq) & (freq_axis <= upper_freq) + # contains, divided by the number of samples. The edge a band + # shares with its neighbour is open on one side, so a bin sitting + # exactly on it is counted once, never in both bands. + above = freq_axis > lower_freq if open_lower else freq_axis >= lower_freq + below = freq_axis < upper_freq if open_upper else freq_axis <= upper_freq + mask = above & below return float(np.sum(np.abs(spectrum[mask]) ** 2) / signal.size) signal = np.asarray(signal) @@ -1144,10 +1159,13 @@ def band_power(lower_freq: float, upper_freq: float) -> float: adj_lower_left = main_lower - bandwidth adj_upper_left = main_lower - # Calculate power in main and adjacent channels + # Calculate power in main and adjacent channels; the main channel + # owns its edges, the adjacent bands start just past them main_power = band_power(main_lower, main_upper) - adj_power_right = band_power(adj_lower_right, adj_upper_right) - adj_power_left = band_power(adj_lower_left, adj_upper_left) + adj_power_right = band_power(adj_lower_right, adj_upper_right, + open_lower=True) + adj_power_left = band_power(adj_lower_left, adj_upper_left, + open_upper=True) # Calculate ACPR in dB acpr_right = 10 * np.log10(adj_power_right / main_power) @@ -1319,8 +1337,14 @@ def update(self, reference: np.ndarray, shorter of the two sets the length, which is what a chain with a filter delay produces. """ - target = np.asarray(reference).ravel() - estimate = np.asarray(detected).ravel() + target = np.atleast_1d(np.asarray(reference)) + estimate = np.atleast_1d(np.asarray(detected)) + # truncate along the *last* axis before raveling, as compute_ser + # does: pooling a batch first would misalign every row after the + # first against its per-row-truncated detection + N = min(target.shape[-1], estimate.shape[-1]) + target = target[..., :N].ravel() + estimate = estimate[..., :N].ravel() length = min(target.size, estimate.size) if length == 0: raise ValueError( diff --git a/src/comnumpy/core/processors.py b/src/comnumpy/core/processors.py index d6304e2..e2be2e8 100644 --- a/src/comnumpy/core/processors.py +++ b/src/comnumpy/core/processors.py @@ -115,6 +115,11 @@ class Upsampler(Processor): filter: Optional[BWFilter] = field(init=False, repr=False, default_factory=lambda: None) def __post_init__(self): + if not 0 <= self.phase < self.L: + raise ValueError( + f"Upsampler: expected 0 <= phase < L, got phase={self.phase} " + f"with L={self.L} -- the phase places the sample inside one " + f"group of L output slots.") self.filter = BWFilter(1/self.L) def forward(self, X: np.ndarray) -> np.ndarray: @@ -326,6 +331,12 @@ class Serial2Parallel(Processor): def __post_init__(self): if not (self.N_sub > 0): raise ValueError("N_sub must be a positive number.") + # validated here, not in the remainder branch of forward: a typo'd + # method must not pass silently just because the lengths divide + if self.method not in ("zero-padding", "truncate"): + raise ValueError( + f"{type(self).__name__}: unknown method {self.method!r}, " + f"expected 'zero-padding' or 'truncate'.") def set_N_sub(self, N_sub: int) -> None: self.N_sub = N_sub @@ -341,15 +352,8 @@ def forward(self, X: np.ndarray) -> np.ndarray: new_shape = X.shape[:-1] + (N_sub*M,) X_processed = np.zeros(new_shape, dtype=X.dtype) X_processed[..., :N] = X - elif self.method == "truncate": + else: # "truncate" -- validated in __post_init__ X_processed = X[..., :M * N_sub] - else: - raise ValueError( - f"{type(self).__name__}: unknown method " - f"{self.method!r}, expected 'zero-padding' or " - f"'truncate' -- the input length {N} is not a " - f"multiple of N_sub={N_sub}, so one of the two has " - f"to say what to do with the remainder.") else: X_processed = X @@ -411,6 +415,13 @@ class Parallel2Serial(Processor): """ name: str = "P2S" + def prepare(self, X: np.ndarray): + if X.ndim < 2: + raise ShapeError( + f"Parallel2Serial expects the Block layout (..., T, F) with " + f"at least two axes, got {X.shape} -- the signal is already " + f"serial, or Serial2Parallel is missing upstream.") + def forward(self, X: np.ndarray) -> np.ndarray: # Pure C-order reshape flattening the last two axes new_shape = X.shape[:-2] + (X.shape[-2] * X.shape[-1],) @@ -624,6 +635,10 @@ def forward(self, X: np.ndarray) -> np.ndarray: if self.validate_input and (not np.isclose(np.real(np.ravel(X)), 0, atol=10**-7).all()): raise ValueError("the input data is not imaginary since the real part is non zero ") Y = np.imag(X) + case _: + raise ValueError( + f"Complex2Real: unknown part {self.part!r}, expected " + f"'real' or 'imag'.") return Y @@ -848,6 +863,11 @@ def prepare(self, x: np.ndarray) -> None: f"{x.ndim}D {x.shape} -- flatten the stream with " f"Parallel2Serial, or apply the block to each stream " f"separately.") + if self.N_start + self.length > len(x): + raise ShapeError( + f"SampleRemover removes samples [{self.N_start}, " + f"{self.N_start + self.length}), which extends past the end " + f"of the {len(x)}-sample input.") def forward(self, x: np.ndarray) -> np.ndarray: y = np.zeros(len(x) - self.length, dtype=x.dtype) @@ -1050,11 +1070,13 @@ def forward(self, x: np.ndarray) -> np.ndarray: if self.selector is None: return x - # turn a plain tuple into the slice it describes - if isinstance(self.selector, tuple): - self.selector = slice(*self.selector) + # turn a plain tuple into the slice it describes -- locally, so + # the configured selector stays the tuple the caller declared + selector = self.selector + if isinstance(selector, tuple): + selector = slice(*selector) - return x[self.selector] + return x[selector] @dataclass(slots=True) diff --git a/src/comnumpy/core/shaping.py b/src/comnumpy/core/shaping.py index 506bcf5..d3eccad 100644 --- a/src/comnumpy/core/shaping.py +++ b/src/comnumpy/core/shaping.py @@ -1200,6 +1200,9 @@ def forward(self, X: np.ndarray) -> np.ndarray: f"not a multiple of it. Pad the stream, or choose a " f"shaper whose n_bits divides it.") flat = bits.reshape(-1, step) + if flat.shape[0] == 0: + # zero blocks is a legal multiple; np.stack refuses an empty list + return np.zeros(bits.shape[:-1] + (0,), dtype=np.int64) blocks = np.stack([self.shaper.encode(row) for row in flat]) return blocks.reshape(bits.shape[:-1] + (-1,)) @@ -1257,6 +1260,9 @@ def forward(self, X: np.ndarray) -> np.ndarray: f"and got {indices.shape[-1]} along the last axis, which " f"is not a multiple of it.") flat = indices.reshape(-1, step) + if flat.shape[0] == 0: + # zero blocks is a legal multiple; np.stack refuses an empty list + return np.zeros(indices.shape[:-1] + (0,), dtype=np.int64) blocks = np.stack([self.shaper.decode(row) for row in flat]) return blocks.reshape(indices.shape[:-1] + (-1,)) diff --git a/src/comnumpy/core/utils.py b/src/comnumpy/core/utils.py index 7d4d7b5..4ee7ceb 100644 --- a/src/comnumpy/core/utils.py +++ b/src/comnumpy/core/utils.py @@ -762,7 +762,12 @@ def soft_projector(z: np.ndarray, alphabet: np.ndarray, sigma2: float, alphabet = np.asarray(alphabet).reshape(1, -1) z = z.reshape(-1, 1) - term1 = np.exp(-(1 / np.real(sigma2)) * np.abs(alphabet - z) ** 2) + # log-sum-exp: subtracting the row-wise best exponent changes nothing + # analytically (it cancels between num and den) but keeps the weights + # finite when every |z - a_m|^2 / sigma2 is large -- the high-SNR + # regime the AMP detectors call this in, where the naive form is 0/0 + exponent = -(1 / np.real(sigma2)) * np.abs(alphabet - z) ** 2 + term1 = np.exp(exponent - np.max(exponent, axis=1, keepdims=True)) if kernel is None: kernel = alphabet diff --git a/src/comnumpy/data.py b/src/comnumpy/data.py index bfc0bc6..4313d32 100644 --- a/src/comnumpy/data.py +++ b/src/comnumpy/data.py @@ -164,7 +164,8 @@ def _transposed(x: np.ndarray, curves: Dict[str, np.ndarray], *, row is what is read across, and a row whose cells change format is unreadable. `x` keeps its own, since it is the header. """ - grid = np.concatenate([values for values in curves.values()]) + stacked = [values for values in curves.values()] + grid = np.concatenate(stacked) if stacked else np.array([]) template = _column_format(grid) heads = _cells(x) body = {} @@ -174,10 +175,12 @@ def _transposed(x: np.ndarray, curves: Dict[str, np.ndarray], *, row.append(_MISSING if math.isnan(value) else template.format(value)) body[name] = row - name_width = max(len(xlabel), *(len(name) for name in body)) - cell_width = max(len(head) for head in heads) + # list form, not star-args: a degenerate table (no curves, empty x) + # must still render its header instead of dying inside max() + name_width = max([len(xlabel)] + [len(name) for name in body]) + cell_width = max([0] + [len(head) for head in heads]) for row in body.values(): - cell_width = max(cell_width, *(len(cell) for cell in row)) + cell_width = max([cell_width] + [len(cell) for cell in row]) def line(first: str, cells: list[str]) -> str: return (first.ljust(name_width) + " " @@ -276,7 +279,7 @@ def format_data(data: Mapping[str, Any], *, headers = [xlabel] + list(curves) columns = [_cells(x)] + [_cells(values) for values in curves.values()] - widths = [max(len(header), *(len(cell) for cell in column)) + widths = [max([len(header)] + [len(cell) for cell in column]) for header, column in zip(headers, columns, strict=True)] def row(fields: list[str]) -> str: diff --git a/src/comnumpy/fec/convolutional.py b/src/comnumpy/fec/convolutional.py index fe2c9fa..5fc6c4a 100644 --- a/src/comnumpy/fec/convolutional.py +++ b/src/comnumpy/fec/convolutional.py @@ -221,7 +221,9 @@ def forward(self, r: np.ndarray) -> np.ndarray: f"got {r.shape} -- check the encoder rate.") T = r.shape[-1] // n n_tail = self.K - 1 if self.terminated else 0 - if T <= n_tail: + # strictly shorter: T == n_tail is the tail alone, the codeword the + # encoder produces for an empty message, and decodes to an empty one + if T < n_tail: raise ShapeError( f"coded length {r.shape[-1]} is shorter than the tail " f"({n_tail} steps of {n} bits) -- nothing to decode.") diff --git a/src/comnumpy/mimo/detectors.py b/src/comnumpy/mimo/detectors.py index 9760886..dce6ccc 100644 --- a/src/comnumpy/mimo/detectors.py +++ b/src/comnumpy/mimo/detectors.py @@ -844,7 +844,10 @@ def fit(self, y: np.ndarray) -> np.ndarray: z_t, self.alphabet, sigma2_t ) # F function in the original publication kernel = np.abs(self.alphabet.reshape(1, -1) - x_t.reshape(-1, 1)) ** 2 - G = soft_projector(z_t, self.alphabet, tau_2, kernel) + # same effective variance as F: G is the posterior variance of + # the *same* distribution (Jeon et al., Algorithm 2); tau_2 + # alone made the weights collapse and the state go NaN + G = soft_projector(z_t, self.alphabet, sigma2_t, kernel) tau_2_old = tau_2 tau_2 = (beta / sigma2) * np.mean(G) term1 = tau_2 / (1 + tau_2_old) @@ -991,7 +994,9 @@ def get_vt_2(self, error: np.ndarray, epsilon: float = 0.001) -> float: R = _required_noise_variance(self.sigma2, block) * np.eye(N_r) H_H = np.conjugate(np.transpose(H)) num = np.sum(np.abs(error) ** 2) - np.trace(R) - den = np.trace(np.matmul(H_H, H)) + # trace(H^H H) is real analytically but complex-typed: take the + # real part rather than lean on numpy's complex ordering + den = np.real(np.trace(np.matmul(H_H, H))) return float(max(num / den, epsilon)) def get_tau_2(self, B: np.ndarray, W: np.ndarray, vt_2: float) -> float: diff --git a/src/comnumpy/ofdm/chains.py b/src/comnumpy/ofdm/chains.py index f6c4c25..a17cb63 100644 --- a/src/comnumpy/ofdm/chains.py +++ b/src/comnumpy/ofdm/chains.py @@ -171,7 +171,9 @@ def __post_init__(self) -> None: N_carriers = self.carrier_type.N_fft else: self.carrier_type = np.asarray(self.carrier_type) - N_carriers = len(self.carrier_type) + # shape[-1], not len(): a 2-D scattered mask is (T_period, N_fft) + # and its first axis is the period, not the FFT size + N_carriers = self.carrier_type.shape[-1] self.h = np.asarray(self.h) self.chain = Sequential([ diff --git a/src/comnumpy/ofdm/processors.py b/src/comnumpy/ofdm/processors.py index af9331b..ccc0139 100644 --- a/src/comnumpy/ofdm/processors.py +++ b/src/comnumpy/ofdm/processors.py @@ -87,6 +87,10 @@ def __post_init__(self): def prepare(self, X: np.ndarray): """Initialize the mask based on the signal shape""" input_length = X.shape[-1] + if self.N_cp > input_length: + raise ShapeError( + f"CyclicPrefixer(N_cp={self.N_cp}) needs blocks of at least " + f"N_cp samples to copy from, got {input_length}") output_mask_length = input_length + self.N_cp input_copy_mask = np.zeros(input_length) @@ -524,7 +528,12 @@ def prepare(self, X: np.ndarray): def forward(self, X: np.ndarray) -> np.ndarray: new_shape = X.shape[:-1] + (self.N,) - Y = np.zeros(new_shape, dtype=X.dtype) + # the buffer must hold the pilots too: real data symbols with + # complex pilots would otherwise silently drop the imaginary part + pilots = np.asarray(self.pilots if self.pilots is not None else []) + dtype = np.result_type(X.dtype, pilots.dtype) if pilots.size \ + else X.dtype + Y = np.zeros(new_shape, dtype=dtype) if self.period == 1: row = self.mask[0] diff --git a/src/comnumpy/optical/compensators.py b/src/comnumpy/optical/compensators.py index 79b236b..00391a9 100644 --- a/src/comnumpy/optical/compensators.py +++ b/src/comnumpy/optical/compensators.py @@ -127,11 +127,17 @@ def __post_init__(self) -> None: beta2_ps2_per_km = compute_beta2(self.lamb, self.D, self.c) beta2 = ((10**-12)**2)*beta2_ps2_per_km # convert into s^2/km K = - beta2 * self.z * (self.fs**2) / 2 - N = int(2 * np.floor(2 * K * np.pi) + 1) - bound = int(np.floor(N / 2)) - n_vect = np.arange(-bound, bound + 1) - coef = np.sqrt(1j / (4 * K * np.pi)) - self.h = coef * np.exp(-1j * (n_vect**2) / (4 * K)) + bound = int(np.floor(2 * K * np.pi)) + if bound == 0: + # below one tap of memory the Savory design is outside its + # validity: its sqrt(1/(4*pi*K)) amplitude diverges as K -> 0 + # although the channel tends to identity. The dispersion is + # not resolvable at this fs, so the compensator is a pass-through. + self.h = np.ones(1, dtype=complex) + else: + n_vect = np.arange(-bound, bound + 1) + coef = np.sqrt(1j / (4 * K * np.pi)) + self.h = coef * np.exp(-1j * (n_vect**2) / (4 * K)) self.K = K def forward(self, x: np.ndarray) -> np.ndarray: diff --git a/src/comnumpy/optical/devices.py b/src/comnumpy/optical/devices.py index d5b48fe..3100409 100644 --- a/src/comnumpy/optical/devices.py +++ b/src/comnumpy/optical/devices.py @@ -163,9 +163,11 @@ def forward(self, x: np.ndarray) -> np.ndarray: y = x else: if self.laser_in: - E_laser = self.laser_in(len(x)) + # shape[-1], not len(x): the local oscillator runs along the + # sample axis and broadcasts over any leading batch axes + E_laser = self.laser_in(x.shape[-1]) else: - E_laser = np.ones(len(x)) + E_laser = np.ones(x.shape[-1]) y = self.sensitivity * x * np.conj(E_laser) y = np.real(y) + 1j*np.imag(y) @@ -393,9 +395,11 @@ def forward(self, x: np.ndarray) -> np.ndarray: m = max(np.max(np.abs(np.real(x))), np.max(np.abs(np.imag(x)))) # normalization coeff if self.laser_in is not None: - E_laser = self.laser_in(len(x)) + # shape[-1], not len(x): the carrier runs along the sample axis + # and broadcasts over any leading batch axes + E_laser = self.laser_in(x.shape[-1]) else: - E_laser = np.ones(len(x)) # ideal laser + E_laser = np.ones(x.shape[-1]) # ideal laser if self.is_ideal: y = -np.pi / 2 * self.Vpp / (2 * self.Vpi) * E_laser * x diff --git a/src/comnumpy/optical/utils.py b/src/comnumpy/optical/utils.py index 2d17548..9231976 100644 --- a/src/comnumpy/optical/utils.py +++ b/src/comnumpy/optical/utils.py @@ -531,6 +531,10 @@ def get_logarithmic_step_size(L_span: float, StPS: int, alpha_dB: float = 0, """ alpha = (np.log(10)/10) * (alpha_dB) alpha_adj = step_log_factor*alpha + if alpha_adj == 0: + # lossless limit of the Sinkin spacing: without attenuation the + # nonlinear phase accrues uniformly, and the formula below is 0/0 + return get_linear_step_size(L_span, StPS) delta = (1-np.exp(-alpha_adj*L_span))/StPS n_vect = 1 + np.arange(StPS) z = -(1/alpha_adj)*np.log((1-n_vect*delta)/(1-(n_vect-1)*delta)) diff --git a/src/comnumpy/serialization.py b/src/comnumpy/serialization.py index b08781b..aecf95b 100644 --- a/src/comnumpy/serialization.py +++ b/src/comnumpy/serialization.py @@ -84,6 +84,9 @@ def _registry() -> dict[str, type]: if (isinstance(attr, type) and issubclass(attr, Processor) and attr is not Processor): registry.setdefault(attr.__name__, attr) + # a nested Sequential is a block too: the encoder recurses into it as + # a dataclass, so the decoder must know how to build one back + registry.setdefault("Sequential", Sequential) registry.update(_EXTRA_BLOCKS) return registry diff --git a/tests/test_audit_regressions.py b/tests/test_audit_regressions.py new file mode 100644 index 0000000..6626d24 --- /dev/null +++ b/tests/test_audit_regressions.py @@ -0,0 +1,311 @@ +r"""Regressions from the v1 code audit. + +Every test here pins a defect found by the full-library audit: five +reviewers swept the code hunting for wrong axes, dtype loss, RNG misuse, +validation gaps and numerical blow-ups, and each confirmed finding was +fixed and locked by one test below. The names say what used to go wrong. +""" +import unittest +import warnings + +import numpy as np + +from comnumpy.core import Sequential +from comnumpy.core.channels import AWGN +from comnumpy.core.filters import BWFilter, SRRCFilter +from comnumpy.core.frames import Deframer, FieldRole, FrameField, FrameStructure +from comnumpy.core.generators import GaussianGenerator, SymbolGenerator +from comnumpy.core.metrics import (ErrorCounter, compute_acpr, compute_ser, + compute_ser_awgn_qam) +from comnumpy.core.processors import (Complex2Real, DataExtractor, Downsampler, + Parallel2Serial, SampleRemover, + Serial2Parallel, Upsampler) +from comnumpy.core.shaping import (DistributionDematcher, DistributionMatcher, + SphereShaper) +from comnumpy.core.utils import get_alphabet, soft_projector +from comnumpy.data import format_data +from comnumpy.exceptions import ShapeError +from comnumpy.serialization import from_json, to_json + + +class TestChainMachinery(unittest.TestCase): + + def test_seed_reaches_a_nested_sequential(self): + chain = Sequential([SymbolGenerator(4), + Sequential([AWGN(sigma2=0.1)], name="inner")]) + first = chain.seed(7)(8) + second = chain.seed(7)(8) + np.testing.assert_array_equal(first, second) + + def test_seed_reaches_a_wrapper_holding_a_chain(self): + from comnumpy.ofdm.chains import OFDMTransmitter + chain = Sequential([SymbolGenerator(4), OFDMTransmitter(4, 1)]) + first = chain.seed(3)(8) + second = chain.seed(3)(8) + np.testing.assert_array_equal(first, second) + + def test_seed_of_a_flat_chain_is_unchanged(self): + # the recursion must not move the child seeds of the existing scheme + chain = Sequential([SymbolGenerator(4), AWGN(sigma2=0.1)]) + chain.seed(7) + expected = [int(c.generate_state(1)[0]) + for c in np.random.SeedSequence(7).spawn(2)] + self.assertEqual([chain[0].seed, chain[1].seed], expected) + + def test_summary_feeds_wiring_and_records_taps(self): + from comnumpy.core.compensators import DataAidedPhaseCompensator + wired = Sequential([SymbolGenerator(4, name="generator"), + AWGN(sigma2=0.01), + DataAidedPhaseCompensator(name="comp")], + wiring={"comp.reference": "generator"}) + rows = wired.summary(64, print_out=False) + self.assertEqual(len(rows), 3) + + tapped = Sequential([SymbolGenerator(4, name="generator"), + AWGN(sigma2=0.1)], taps=["generator"]) + tapped.summary(8, print_out=False) + self.assertEqual(tapped.tap("generator").shape, (8,)) + + def test_an_int_keyed_callback_fires(self): + fired = [] + Sequential([SymbolGenerator(4, seed=1)], + callbacks={0: fired.append})(5) + self.assertEqual(len(fired), 1) + + def test_a_tap_does_not_survive_its_reconfiguration(self): + chain = Sequential([SymbolGenerator(4, name="generator"), + AWGN(sigma2=0.1)], taps=["generator", "awgn"]) + chain(8) + chain.taps = ["awgn"] + chain(8) + with self.assertRaises(KeyError): + chain.tap("generator") + + def test_a_nested_chain_round_trips_through_json(self): + chain = Sequential([SymbolGenerator(4, seed=1), + Sequential([AWGN(sigma2=0.1, seed=2)], + name="inner")]) + rebuilt = from_json(to_json(chain)) + np.testing.assert_array_equal(rebuilt(8), chain(8)) + + +class TestCoreBlocks(unittest.TestCase): + + def test_framefield_freezes_a_copy_not_the_callers_array(self): + values = np.array([1.0, -1.0]) + field = FrameField("S", FieldRole.SYNC, values) + values[0] = 5.0 # must stay writable + assert field.values is not None + self.assertEqual(field.values[0], 1.0) + + def test_deframer_fields_own_their_data(self): + frame = FrameStructure(( + FrameField("SYNC", FieldRole.SYNC, np.array([1.0, -1.0])), + FrameField("PAYLOAD", FieldRole.PAYLOAD, length=3)), + standard="example") + deframer = Deframer(frame) + buffer = np.array([[1.0, -1.0, 1.0, 2.0, 3.0]]) + deframer(buffer) + buffer[:] = 999.0 + np.testing.assert_array_equal(deframer.get_field("SYNC"), + [[1.0, -1.0]]) + + def test_complex2real_rejects_an_unknown_part(self): + with self.assertRaisesRegex(ValueError, "real.*imag"): + Complex2Real(part="magnitude")(np.array([1 + 2j])) # type: ignore[arg-type] + + def test_serial2parallel_validates_method_at_construction(self): + with self.assertRaisesRegex(ValueError, "zero-padding"): + Serial2Parallel(3, method="trunc") # type: ignore[arg-type] + + def test_parallel2serial_refuses_a_serial_signal(self): + with self.assertRaises(ShapeError): + Parallel2Serial()(np.arange(5)) + + def test_upsampler_validates_the_phase(self): + with self.assertRaisesRegex(ValueError, "phase"): + Upsampler(L=2, phase=2) + + def test_sampleremover_validates_the_window(self): + with self.assertRaises(ShapeError): + SampleRemover(N_start=4, length=3)(np.arange(5)) + + def test_dataextractor_keeps_its_declared_selector(self): + extractor = DataExtractor((2, 8)) + extractor(np.arange(10)) + self.assertEqual(extractor.selector, (2, 8)) + + def test_generators_accept_numpy_sizes_and_refuse_bool(self): + self.assertEqual(SymbolGenerator(4, seed=1)(np.int64(5)).shape, (5,)) + with self.assertRaises(ValueError): + GaussianGenerator()(True) + + def test_awgn_rejects_a_negative_variance(self): + with self.assertRaisesRegex(ValueError, "negative"): + AWGN(sigma2=-0.1, seed=0)(np.zeros(4, dtype=complex)) + + def test_the_matcher_pair_handles_an_empty_stream(self): + shaper = SphereShaper(np.array([1.0, 3.0]), length=8, max_energy=3) + empty_bits = np.array([], dtype=int) + self.assertEqual(DistributionMatcher(shaper)(empty_bits).shape, (0,)) + self.assertEqual(DistributionDematcher(shaper)(empty_bits).shape, (0,)) + + def test_the_antialias_filter_keeps_a_real_signal_real(self): + x = np.cos(2 * np.pi * 0.01 * np.arange(64)) + self.assertFalse(np.iscomplexobj(Downsampler(L=4, use_filter=True)(x))) + self.assertFalse(np.iscomplexobj(BWFilter(0.5)(x))) + + +class TestFiltersAndNumerics(unittest.TestCase): + + def test_srrc_at_zero_rolloff_is_the_sinc_pulse(self): + h = SRRCFilter(0.0, 4).h() + n = np.arange(-40, 41) + sinc = np.sinc(n / 4) + sinc = sinc / np.sqrt(np.sum(sinc ** 2)) + np.testing.assert_allclose(h, sinc, atol=1e-12) + + def test_srrc_fft_method_names_the_minimum_length(self): + with self.assertRaisesRegex(ShapeError, "method='lfilter'"): + SRRCFilter(0.25, 4, method="fft")(np.ones(64, dtype=complex)) + + def test_soft_projector_survives_the_high_snr_regime(self): + alphabet = get_alphabet("QAM", 4) + with warnings.catch_warnings(): + warnings.simplefilter("error") + estimate = soft_projector(np.array([5 + 5j]), alphabet, 0.01) + self.assertFalse(np.isnan(estimate).any()) + + +class TestMetrics(unittest.TestCase): + + def test_errorcounter_truncates_along_the_last_axis(self): + reference = np.tile(np.arange(10), (2, 1)) + detected = reference[:, :9] # one symbol of filter delay + counter = ErrorCounter().update(reference, detected) + self.assertEqual(counter.n_errors, 0) + + def test_ser_awgn_qam_rejects_cross_constellations(self): + with self.assertRaisesRegex(ValueError, "square"): + compute_ser_awgn_qam(32, 10.0) + + def test_compute_ser_refuses_an_empty_comparison(self): + with self.assertRaisesRegex(ValueError, "no symbols"): + compute_ser(np.array([]), np.array([])) + + def test_acpr_counts_an_edge_bin_once(self): + # a tone exactly on the shared band edge belongs to the main + # channel only; counted twice it read as acpr = 0 dB + tone = np.exp(2j * np.pi * 10 * np.arange(1000) / 100) + acpr_right, _ = compute_acpr(tone, 20.0, 100.0) + self.assertLess(acpr_right, -100.0) + + +class TestOFDM(unittest.TestCase): + + def test_receiver_reads_nfft_off_a_scattered_mask(self): + from comnumpy.ofdm.chains import OFDMReceiver, OFDMTransmitter + mask = np.array([[1, 2, 1, 1], [1, 1, 2, 1]]) + transmitter = OFDMTransmitter(3, 1, carrier_type=mask, + pilots=np.array([1.0 + 0j])) + receiver = OFDMReceiver(3, 1, carrier_type=mask) + x = np.arange(12) + 0j + np.testing.assert_allclose(receiver(transmitter(x)), x, atol=1e-12) + + def test_allocator_keeps_complex_pilots_on_real_data(self): + from comnumpy.ofdm.processors import CarrierAllocator + allocator = CarrierAllocator(carrier_type=np.array([1, 2, 1]), + pilots=np.array([1j])) + np.testing.assert_array_equal(allocator(np.array([1.0, 2.0])), + np.array([1.0, 1j, 2.0])) + + def test_prefixer_refuses_a_prefix_longer_than_the_block(self): + from comnumpy.ofdm.processors import CyclicPrefixer + with self.assertRaises(ShapeError): + CyclicPrefixer(N_cp=6)(np.arange(4)) + + +class TestDataDisplay(unittest.TestCase): + + def test_degenerate_tables_render_their_header(self): + empty = {"x": [], "curves": {"a": []}} + self.assertIn("a", format_data(empty)) + self.assertIn("a", format_data(empty, transpose=True)) + self.assertIn("x", format_data({"x": [1, 2], "curves": {}}, + transpose=True)) + + +class TestMIMOAndFEC(unittest.TestCase): + + def test_amp_no_longer_diverges_to_nan(self): + from comnumpy.mimo.detectors import ApproximateMessagePassingDetector + alphabet = get_alphabet("QAM", 4) + alphabet = alphabet / np.sqrt(np.mean(np.abs(alphabet) ** 2)) + for seed in range(10): + rng = np.random.default_rng(seed) + H = (rng.normal(size=(4, 4)) + + 1j * rng.normal(size=(4, 4))) / np.sqrt(8) + sent = alphabet[rng.integers(0, 4, size=(4, 10))] + Y = H @ sent + np.sqrt(0.25) * (rng.normal(size=(4, 10)) + + 1j * rng.normal(size=(4, 10))) + detector = ApproximateMessagePassingDetector( + alphabet, H=H, sigma2=0.5, N_it=20) + detector(Y) + for n in range(10): + soft = detector.fit(Y[:, n]) + self.assertFalse(np.isnan(soft).any(), + f"AMP state went NaN at seed {seed}") + + def test_oamp_mmse_raises_no_complex_warning(self): + from comnumpy.mimo.detectors import ( + OrthogonalApproximateMessagePassingDetector) + alphabet = get_alphabet("QAM", 4) + alphabet = alphabet / np.sqrt(np.mean(np.abs(alphabet) ** 2)) + rng = np.random.default_rng(1) + H = (rng.normal(size=(4, 4)) + 1j * rng.normal(size=(4, 4))) / 2 + Y = H @ alphabet[rng.integers(0, 4, size=(4, 5))] + detector = OrthogonalApproximateMessagePassingDetector( + alphabet, H=H, sigma2=0.5, type="MMSE", N_it=5) + with warnings.catch_warnings(): + warnings.simplefilter("error") + detector(Y) + + def test_viterbi_round_trips_an_empty_message(self): + from comnumpy.fec.convolutional import (ConvolutionalEncoder, + ViterbiDecoder) + encoder = ConvolutionalEncoder((0o5, 0o7)) + decoder = ViterbiDecoder((0o5, 0o7)) + decoded = decoder(encoder(np.array([], dtype=int))) + self.assertEqual(decoded.shape, (0,)) + + +class TestOptical(unittest.TestCase): + + def test_logarithmic_steps_of_a_lossless_fiber_are_uniform(self): + from comnumpy.optical.utils import get_logarithmic_step_size + steps = get_logarithmic_step_size(80.0, 5, alpha_dB=0.0) + np.testing.assert_allclose(steps, 16.0) + + def test_the_fir_cd_compensator_degrades_to_identity(self): + from comnumpy.optical.compensators import ( + ChromaticDispersionFIRCompensator) + for z_km in (0.0, 1.0, 10.0): + compensator = ChromaticDispersionFIRCompensator(z_km, fs=20e9) + assert compensator.h is not None + self.assertEqual(len(compensator.h), 1) + self.assertAlmostEqual(float(np.sum(np.abs(compensator.h) ** 2)), + 1.0) + + def test_the_lo_broadcasts_over_a_polarization_pair(self): + from comnumpy.optical.devices import (Laser, MachZehnderModulator, + Optical90HybridCircuit) + pair = np.ones((2, 8), dtype=complex) + self.assertEqual(Optical90HybridCircuit(is_ideal=False)(pair).shape, + (2, 8)) + modulator = MachZehnderModulator(is_ideal=True, + laser_in=Laser(0, theta0=0.0)) + self.assertEqual(modulator(pair).shape, (2, 8)) + + +if __name__ == "__main__": + unittest.main()