diff --git a/.claude/skills/comnumpy-tutorial/SKILL.md b/.claude/skills/comnumpy-tutorial/SKILL.md index 1a26062..a6d25a8 100644 --- a/.claude/skills/comnumpy-tutorial/SKILL.md +++ b/.claude/skills/comnumpy-tutorial/SKILL.md @@ -54,7 +54,7 @@ then learns the wrong lesson -- they came to see how the library is used. `compute_ccdf`, `compute_metric_awgn_theo` and `compute_metric_rayleigh_theo` for the closed forms, `constellation_capacity` and `bicm_capacity` for the rates. -- Monte-Carlo: `sweep(chain, param, values, metrics, stimulus, seed=)`. +- Monte-Carlo: `monte_carlo(chain, param, values, metrics, stimulus, seed=)`. - Figures: `plot_iq`, `plot_error_rate`, `plot_spectrum`, `plot_time`, `plot_channel_response`, and the channels' own `info()` / `plot()`. - If the library makes something awkward, that is a finding about the @@ -63,9 +63,39 @@ then learns the wrong lesson -- they came to see how the library is used. ## Structuring the example script -- One `get_*` function per assembly the page talks about, so that what - varies is an argument. `get_full_chain(n_spans)` beats building a chain - inline six times. +### Naming the chain + +Three rules, and they are about where the reader looks and what the +signature promises. + +1. **A function only when the chain is built more than once** -- several + values of a parameter, several variants, a Monte-Carlo. A chain built + once and used once is written inline: wrapping four blocks in a + function called once, with no argument, is ceremony, and it puts an + indirection between the reader and their first chain. +2. **Always `get_()`**, one prefix, the noun saying what comes + back. Which noun follows from how the page is organized: + + | the page… | names | + |---|---| + | has one chain | `get_chain(...)` | + | cuts one chain in two | `get_channel()` and `get_receiver()` | + | compares two structurally different chains | `get_uncoded_chain()`, `get_coded_chain(soft)` | + | builds a sub-assembly, not a whole chain | the part: `get_transmitter()` | + + Two chains that differ by a block or two are **one** function with a + parameter, not two functions. Two chains that differ in structure -- + four blocks against six -- are two functions, because merging them + would hide behind a flag exactly what the page is comparing. +3. **What varies is a parameter, never a module global.** Constants of + the page (the fibre, the constellation, the roll-off) may stay + global; anything the page *varies* goes through the signature. A + function whose signature announces three arguments and silently reads + eleven more is a closure wearing a `def`, and the reader cannot tell + from the call site what changed. + +### The rest of the script + - **Prefer one chain, transmitter to decision.** A full chain is what the reader will write. Cutting the chain into pieces is justified only when a Monte-Carlo re-runs an expensive stage that could be run once -- and then diff --git a/.github/scripts/run_slow_examples.py b/.github/scripts/run_slow_examples.py new file mode 100644 index 0000000..8f8e2c7 --- /dev/null +++ b/.github/scripts/run_slow_examples.py @@ -0,0 +1,94 @@ +"""Run the examples `tests/test_examples_run.py` marks as SLOW. + +The smoke test skips them so that a pull request does not wait minutes +per script. Skipping is the right call there and the wrong one +everywhere: for a while nothing ran them at all, and a regression lived +in `CD_compensation_part1.py` for two commits -- it raised `TypeError: +'Constellation' object is not subscriptable` on its first pass while +every check stayed green. + +The list is not repeated here. It is imported from the test module, so +adding an entry to `SLOW` is enough to have it covered, and an entry +that is removed stops being run without anyone editing a workflow. + +Each script runs in its own process, from its own directory, with the +same sandbox the smoke test uses: the examples hardcode +`../../docs/tutorials/img/`, and a scheduled job has no business +rewriting figures that are committed. Every script runs even if an +earlier one failed -- one broken example must not hide the state of the +other eight. +""" +import os +import shutil +import subprocess +import sys +import tempfile +import time +from pathlib import Path + +REPO = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(REPO)) + +from tests.test_examples_run import SLOW # noqa: E402 -- after sys.path + +TIMEOUT_S = 1800 + + +def sandbox() -> Path: + """A throwaway copy of `examples/`, with the directories it writes to.""" + root = Path(tempfile.mkdtemp(prefix="comnumpy_slow_")) + shutil.copytree(REPO / "examples", root / "examples", + ignore=shutil.ignore_patterns("__pycache__")) + for pattern in ("docs/**/img", "docs/**/mermaid"): + for directory in REPO.glob(pattern): + (root / directory.relative_to(REPO)).mkdir(parents=True, + exist_ok=True) + (root / "validation" / "figures").mkdir(parents=True, exist_ok=True) + return root + + +def main() -> int: + only = os.environ.get("ONLY", "").strip() + names = sorted(name for name in SLOW if only in name) + if not names: + print(f"no SLOW example matches {only!r}; known: {sorted(SLOW)}") + return 1 + + root = sandbox() + env = os.environ.copy() + env["MPLBACKEND"] = "Agg" + env["PYTHONPATH"] = os.pathsep.join( + p for p in (str(REPO / "src"), env.get("PYTHONPATH", "")) if p) + + failures = [] + for name in names: + script = root / "examples" / name + print(f"::group::{name}", flush=True) + start = time.perf_counter() + try: + done = subprocess.run([sys.executable, script.name], + cwd=script.parent, env=env, + capture_output=True, text=True, + timeout=TIMEOUT_S) + code, output = done.returncode, done.stdout + done.stderr + except subprocess.TimeoutExpired: + code, output = 1, f"did not finish within {TIMEOUT_S} s" + elapsed = time.perf_counter() - start + print(output[-4000:]) + print("::endgroup::", flush=True) + status = "ok" if code == 0 else "FAILED" + print(f"{status:6s} {name:42s} {elapsed:7.1f} s", flush=True) + if code != 0: + failures.append((name, output[-2000:])) + + if failures: + print(f"\n{len(failures)} of {len(names)} slow examples failed:") + for name, output in failures: + print(f"\n--- {name} ---\n{output}") + return 1 + print(f"\nall {len(names)} slow examples ran") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/.github/workflows/slow-examples.yml b/.github/workflows/slow-examples.yml new file mode 100644 index 0000000..2bbf6b3 --- /dev/null +++ b/.github/workflows/slow-examples.yml @@ -0,0 +1,46 @@ +name: Slow examples + +# The nine scripts in `SLOW` are skipped by the smoke test, which is what +# `SLOW` is for: they cost minutes each and a pull request should not wait +# for them. But nothing else ran them either, so for a while nothing ran +# them at all -- and a regression lived there for two commits, invisible. +# `CD_compensation_part1.py` raised `TypeError: 'Constellation' object is +# not subscriptable` on its first pass and every check stayed green. +# +# So they run here instead: once a week, and on demand. Weekly is the +# right cadence for a guard against drift rather than against a specific +# change -- the whole set is about 15 minutes of runner time, and a +# failure names the script and the traceback. + +on: + schedule: + - cron: "17 4 * * 1" # Monday, 04:17 UTC + workflow_dispatch: + inputs: + example: + description: "Substring of one example to run; empty runs all of them" + required: false + default: "" + +jobs: + slow-examples: + runs-on: ubuntu-latest + timeout-minutes: 45 + + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Install package + run: | + python -m pip install --upgrade pip + pip install -e .[test] + + - name: Run the examples the smoke test skips + env: + MPLBACKEND: Agg + ONLY: ${{ inputs.example }} + run: python .github/scripts/run_slow_examples.py diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index 4b263b9..b3ade90 100644 --- a/ARCHITECTURE.md +++ b/ARCHITECTURE.md @@ -159,7 +159,7 @@ Le détail fichier par fichier est en annexe A. | # | Décision | Motif | Alternatives rejetées | Statut | |---|----------|-------|-----------------------|--------| -| D22 | **Cycle de vie inspiré de scikit-learn.** Signature unique `fit(X, y=None)` où `y` est la référence connue (préambule, pilotes) ; `y=None` = estimation aveugle. `fit` retourne `self`. `partial_fit(X)` pour les algorithmes adaptatifs (CMA, RDE, décision-dirigé). `forward()` reste l'interface unique du `Processor`. *Amendé v0.5* : **vocabulaire arrêté** — le signal connu contre lequel un estimateur se cale s'appelle `reference` partout (paramètre du bloc, `sweep(reference=…)`, `wiring`), et les blocs qui en demandent un sont `DataAided*`, par opposition à la famille `Blind*` déjà nommée ainsi. Les classes `TrainedBased*` sont renommées en conséquence | La convention sklearn est la plus répandue en Python scientifique, et le couple supervisé/non-supervisé se superpose exactement à *data-aided*/*blind*. L'audit montre que `fit()` existe déjà sur 8 classes, mais que le second argument y désigne tour à tour `x_target`, `x_preamble` et `w0` (une initialisation d'algorithme, pas une donnée). L'amendement v0.5 clôt le même problème côté *nom* : `target_data` était un emprunt à scikit-learn, où la cible est ce qu'on **prédit** — ici c'est le signal **connu**. `reference` est le terme du domaine (3GPP *reference signal*), et surtout celui que `sweep()` employait déjà pour l'exacte même chose. Note : les *métriques* gardent `X_target` / `X_detected`, qui est la paire standard d'une mesure d'erreur (`y_true`/`y_pred` chez sklearn) — rôle différent, mot différent, les deux justes | Renommer `forward` en `transform` (fragmenterait `Sequential`) ; garder l'état actuel | **Amendé v0.5** | +| D22 | **Cycle de vie inspiré de scikit-learn.** Signature unique `fit(X, y=None)` où `y` est la référence connue (préambule, pilotes) ; `y=None` = estimation aveugle. `fit` retourne `self`. `partial_fit(X)` pour les algorithmes adaptatifs (CMA, RDE, décision-dirigé). `forward()` reste l'interface unique du `Processor`. *Amendé v0.5* : **vocabulaire arrêté** — le signal connu contre lequel un estimateur se cale s'appelle `reference` partout (paramètre du bloc, `monte_carlo(reference=…)`, `wiring`), et les blocs qui en demandent un sont `DataAided*`, par opposition à la famille `Blind*` déjà nommée ainsi. Les classes `TrainedBased*` sont renommées en conséquence | La convention sklearn est la plus répandue en Python scientifique, et le couple supervisé/non-supervisé se superpose exactement à *data-aided*/*blind*. L'audit montre que `fit()` existe déjà sur 8 classes, mais que le second argument y désigne tour à tour `x_target`, `x_preamble` et `w0` (une initialisation d'algorithme, pas une donnée). L'amendement v0.5 clôt le même problème côté *nom* : `target_data` était un emprunt à scikit-learn, où la cible est ce qu'on **prédit** — ici c'est le signal **connu**. `reference` est le terme du domaine (3GPP *reference signal*), et surtout celui que `monte_carlo()` employait déjà pour l'exacte même chose. Note : les *métriques* gardent `X_target` / `X_detected`, qui est la paire standard d'une mesure d'erreur (`y_true`/`y_pred` chez sklearn) — rôle différent, mot différent, les deux justes | Renommer `forward` en `transform` (fragmenterait `Sequential`) ; garder l'état actuel | **Amendé v0.5** | | D23 | **Underscore final pour toute grandeur estimée à partir de données** : `theta_`, `h_`, `gain_`, `H_`, `n_iter_`. Les quantités dérivées des seuls paramètres n'en portent pas. `NotFittedError` si `forward` précède `fit` | P7 ; l'ambiguïté existe déjà et mord : le canal *estimé* `self.H` du compensateur MIMO est homonyme du canal *configuré* `H` de `FlatMIMOChannel`. Aujourd'hui un `forward` avant `fit` applique silencieusement `theta = 0` | Préfixe `est_` ; documentation seule | **Acté** | | D24 | **Frontière explicite avec scikit-learn** : on emprunte les noms de méthodes et le cycle de vie, **jamais la disposition des données**. Le layout `(n_samples, n_features)` est incompatible avec D2 et n'entre pas dans la lib. Pas de `score()` | sklearn met les échantillons en lignes (axe 0) ; D2 met le temps en dernier axe. Sans clause explicite, « faisons comme sklearn » ferait rentrer la mauvaise convention par la fenêtre. Et `score()` est un « plus grand est meilleur », alors que BER/SER/EVM/MSE sont des « plus petit est meilleur » | Adoption intégrale de l'API sklearn | **Acté** | @@ -416,12 +416,12 @@ Sequential( | # | Décision | Motif | Alternatives rejetées | Statut | |---|----------|-------|-----------------------|--------| | D34 | **Toute chaîne est reconfigurable après construction** : `chain.set_params(**{"awgn.value": 12})` modifie le paramètre désigné en notation pointée `id_bloc.champ` et redéclenche le précalcul paramétrique du bloc (`__post_init__` ou équivalent). Emprunt du `set_params` de scikit-learn, adressage par les `id` de D31. Corollaire rendu possible par la restriction sur `__init__` de la §4.3 : puisque le constructeur ne fait ni travail dépendant des données ni création de RNG, reconstruire l'état d'un bloc après modification d'un champ est sûr et bon marché | Tout balayage de paramètre (BER vs SNR, pénalité vs distance) exige de modifier un paramètre et relancer. Aujourd'hui la seule voie est de reconstruire la chaîne entière à la main, en dupliquant sa définition — la classe de bugs « configurations divergentes » de D18, appliquée au temps : la chaîne du point de mesure n s'écarte silencieusement de celle du point 0. `set_params` + `Sequential.seed` (D6) forment le couple minimal : reconfigurer, re-semer, relancer | Mutation directe des attributs (saute le précalcul → état incohérent, exactement le piège que D23 nomme) ; reconstruction depuis le JSON D31 avec surcharge (valide, mais plus coûteux à écrire pour l'usage interactif ; reste la voie recommandée pour les scénarios sérialisés) | **Acté** | -| D35 | **Le balayage est une fonction, pas un framework.** Forme cible : `sweep(chain, param, values, metrics, seed, n_frames)` retournant un dict de tableaux — le patron `GridSearchCV`/`cross_validate` de sklearn, cohérent avec l'emprunt D22–D24. **L'implémentation est différée** : elle n'entre dans `src/` qu'après que trois scripts `validation/` (D7) auront exhibé le même squelette recopié à l'identique, et sa forme finale sera extraite de ces scripts, pas conçue a priori. Le patron *Trainer* de PyTorch Lightning (inversion de contrôle par hooks) est **explicitement rejeté** | Un balayage Monte-Carlo est une boucle `for` sans état, sans optimiseur, embarrassingly parallel : l'inversion de contrôle rendrait la boucle invisible, ce qui contredit P2 frontalement — un étudiant doit pouvoir lire *où* le SNR est balayé. À l'inverse, tout coder à la main pour chaque figure recopie le triplet construire/balayer/collecter. La fonction sklearn-like est le point d'équilibre : la boucle reste visible dans son implémentation (quinze lignes lisibles), l'appelant n'écrit qu'une ligne. D31/D32/D6 fournissent déjà sérialisation et reproductibilité ; un « scénario » n'est que chaîne JSON + spec de balayage + graine + noms de métriques — un format de fichier avant d'être du code | *Trainer* à hooks (inversion de contrôle injustifiée, boucle invisible) ; classe `Simulation` avec cycle de vie (sur-abstraction pour une boucle sans état) ; implémentation immédiate (violerait le principe de complexité à la demande : la forme doit sortir d'un motif constaté dans `validation/`, pas d'une anticipation) | **Acté — implémentation différée, déclencheur : 3ᵉ script `validation/` au squelette identique** | +| D35 | **Le balayage est une fonction, pas un framework.** Forme cible : `monte_carlo(chain, param, values, metrics, seed, n_frames)` retournant un dict de tableaux — le patron `GridSearchCV`/`cross_validate` de sklearn, cohérent avec l'emprunt D22–D24. **L'implémentation est différée** : elle n'entre dans `src/` qu'après que trois scripts `validation/` (D7) auront exhibé le même squelette recopié à l'identique, et sa forme finale sera extraite de ces scripts, pas conçue a priori. Le patron *Trainer* de PyTorch Lightning (inversion de contrôle par hooks) est **explicitement rejeté** | Un balayage Monte-Carlo est une boucle `for` sans état, sans optimiseur, embarrassingly parallel : l'inversion de contrôle rendrait la boucle invisible, ce qui contredit P2 frontalement — un étudiant doit pouvoir lire *où* le SNR est balayé. À l'inverse, tout coder à la main pour chaque figure recopie le triplet construire/balayer/collecter. La fonction sklearn-like est le point d'équilibre : la boucle reste visible dans son implémentation (quinze lignes lisibles), l'appelant n'écrit qu'une ligne. D31/D32/D6 fournissent déjà sérialisation et reproductibilité ; un « scénario » n'est que chaîne JSON + spec de balayage + graine + noms de métriques — un format de fichier avant d'être du code | *Trainer* à hooks (inversion de contrôle injustifiée, boucle invisible) ; classe `Simulation` avec cycle de vie (sur-abstraction pour une boucle sans état) ; implémentation immédiate (violerait le principe de complexité à la demande : la forme doit sortir d'un motif constaté dans `validation/`, pas d'une anticipation) | **Acté — implémentation différée, déclencheur : 3ᵉ script `validation/` au squelette identique** | **Esquisse (D35), non normative tant que le déclencheur n'est pas atteint :** ```python -results = sweep( +results = monte_carlo( chain, param="awgn.snr_dB", # dotted path, resolved by D34 values=np.arange(0, 21, 2), @@ -452,7 +452,7 @@ vérifiable en CI, pas une intention. | D38 | **Hiérarchie d'exceptions à trois membres, pas plus** : `ComnumpyError(Exception)` racine ; `ShapeError(ComnumpyError, ValueError)` pour tout refus de `prepare()` ; `NotFittedError(ComnumpyError, RuntimeError)` pour D23. Le double héritage préserve les `except ValueError` existants. Le gabarit de message de la §8 devient **normatif** : *constat, attendu, action* — « got X, expected Y, do Z ». Jamais d'`assert` pour la validation d'entrée utilisateur (`python -O` les supprime) | P5 exige de refuser explicitement ; ceci outille le refus. `except ComnumpyError` sépare « j'ai mal configuré ma chaîne » de « bug numpy sous-jacent », distinction qu'un `ValueError` nu ne fait pas. Trois classes suffisent : chaque exception supplémentaire est une décision de plus imposée à l'utilisateur au moment du `except` | Une exception par module (taxonomie sans usage) ; codes d'erreur ; `ValueError` nus (état actuel — indistinguables de numpy) | **Acté** | | D39 | **Un seul outil de qualité : `ruff`** (lint **et** format, remplace flake8 + black + isort), configuré dans `pyproject.toml`, exécuté par pre-commit et en CI. Couverture mesurée avec **seuil cliquet** : le pourcentage constaté devient le minimum exigé, il ne peut que monter. Complète C1–C5 : matrice 3.11/3.12/3.13, doctests bloquants (D8), étape `pip install .` + import nu | Zéro débat de style pour un projet à mainteneur unique qui espère des contributions JOSS : le formateur tranche, la revue parle d'architecture. Le cliquet évite le double piège des seuils : l'objectif arbitraire (90 % ?) et les tests de remplissage écrits pour l'atteindre — on interdit seulement la régression, ce qui est le vrai risque constaté (`optical/` à 0 %) | black + flake8 + isort (trois outils, trois configs, conflits connus) ; seuil de couverture fixe ; pylint (lent, verbeux) | **Acté** | | D40 | **Python moderne comme filet, budget de simplicité comme contrat.** (a) Toutes les dataclasses passent en `slots=True` ; les paramètres optionnels sont `kw_only=True` ; les objets valeur (`CarrierAllocation`, `FrameStructure`, `FrameField`) restent `frozen`. (b) Dans toute signature publique, seul le premier argument (la donnée ou le paramètre principal) est positionnel ; le reste est keyword-only. (c) **Budget d'ergonomie normatif, vérifié par doctest** : l'exemple du README tient en ≤ 8 lignes imports compris ; aucune tâche du tutoriel n'importe un module préfixé `_` ; toute chaîne de la doc se construit sans variable intermédiaire autre que les blocs eux-mêmes | `slots` transforme la classe de défauts S1/S2 (typo d'attribut, `__post__init__` jamais appelé silencieusement) d'un bug silencieux en `AttributeError` immédiate — c'est le même mouvement que P5, appliqué à la lib elle-même. Le kw-only rend les appels auto-documentés (`AWGN(snr_dB=15)` se lit sans la doc, P7) et rend l'ajout de paramètres non cassant. Le budget (c) rend « simple à utiliser » **falsifiable** : si l'exemple ne tient plus en 8 lignes, c'est l'API qui a tort, pas l'exemple | `attrs` (dépendance, viole P1) ; `pydantic` (validation runtime coûteuse sur des tableaux, dépendance lourde) ; budget en prose non testé (vœu pieux — exactement ce que P3 interdit ailleurs) | **Acté** | -| D41 | **Paramétrages alternatifs : un kwarg nommé par paramétrage, mutuellement exclusifs — jamais de couple `value`/`unit`.** Un bloc admettant plusieurs façons équivalentes d'être configuré expose un keyword par voie (`AWGN(snr_dB=15)` **ou** `AWGN(sigma2=0.01)`), tous `None` par défaut, avec validation « exactement un » dans `__post_init__` et message énumérant les choix. Le nom du kwarg porte l'unité (`snr_dB`, `sigma2`, `length_km`) : il n'existe **aucun** paramètre `unit` dans la lib | Le couple `value`/`unit` actuel est de la *primitive obsession* : deux paramètres sans sens séparés, où la faute de frappe (`"snr_db"`) traverse jusqu'au runtime. Surtout, `snr_dB` et `sigma2` ne sont pas la même grandeur dans deux unités — l'un est relatif à la puissance du signal, l'autre absolu ; ils changent le *calcul* de `forward`, pas une conversion. Un kwarg par voie rend la sémantique visible dans l'appel (P7), la faute de frappe devient `TypeError` immédiate (D40a/b), le JSON D31 sérialise `{"snr_dB": 15}` sans champ discriminant, et `sweep(param="awgn.snr_dB")` (D35) devient auto-descriptif là où `awgn.value` était muet. Précédent : `timedelta(hours=…, minutes=…)` | `Literal["snr_dB", "sigma2"]` sur `unit` (corrige le typage, garde le couple, masque la différence relatif/absolu) ; `IntEnum NoiseUnit` (verbeux à l'appel, un import de plus, conversion JSON nécessaire — les enums D15/D21 se justifient pour des *données*, pas un sélecteur binaire) ; `pint`/objets quantité (viole P1) ; surcharge par type (indécidable : tout est `float`) | **Acté** | +| D41 | **Paramétrages alternatifs : un kwarg nommé par paramétrage, mutuellement exclusifs — jamais de couple `value`/`unit`.** Un bloc admettant plusieurs façons équivalentes d'être configuré expose un keyword par voie (`AWGN(snr_dB=15)` **ou** `AWGN(sigma2=0.01)`), tous `None` par défaut, avec validation « exactement un » dans `__post_init__` et message énumérant les choix. Le nom du kwarg porte l'unité (`snr_dB`, `sigma2`, `length_km`) : il n'existe **aucun** paramètre `unit` dans la lib | Le couple `value`/`unit` actuel est de la *primitive obsession* : deux paramètres sans sens séparés, où la faute de frappe (`"snr_db"`) traverse jusqu'au runtime. Surtout, `snr_dB` et `sigma2` ne sont pas la même grandeur dans deux unités — l'un est relatif à la puissance du signal, l'autre absolu ; ils changent le *calcul* de `forward`, pas une conversion. Un kwarg par voie rend la sémantique visible dans l'appel (P7), la faute de frappe devient `TypeError` immédiate (D40a/b), le JSON D31 sérialise `{"snr_dB": 15}` sans champ discriminant, et `monte_carlo(param="awgn.snr_dB")` (D35) devient auto-descriptif là où `awgn.value` était muet. Précédent : `timedelta(hours=…, minutes=…)` | `Literal["snr_dB", "sigma2"]` sur `unit` (corrige le typage, garde le couple, masque la différence relatif/absolu) ; `IntEnum NoiseUnit` (verbeux à l'appel, un import de plus, conversion JSON nécessaire — les enums D15/D21 se justifient pour des *données*, pas un sélecteur binaire) ; `pint`/objets quantité (viole P1) ; surcharge par type (indécidable : tout est `float`) | **Acté** | **Piège d'implémentation (D40a, à connaître avant le refactor).** Avec `slots=True`, tout attribut créé dans `__post_init__` doit être déclaré @@ -974,6 +974,23 @@ qui dit ce qu'un code récupérera. --- +### 4.19 Où vit l'état dérivé, et comment on change un paramètre (nouveau) + +| # | Décision | Motif | Alternatives rejetées | Statut | +|---|----------|-------|-----------------------|--------| +| D50 | **Un bloc précalcule à deux endroits, et un seul chemin change un paramètre.** (a) `__post_init__` valide, normalise et précalcule ce qui ne dépend **que des paramètres du bloc** — coefficients de filtre, alphabet coercé en tableau, masque d'allocation. (b) `prepare(x)` dérive ce qui dépend du **signal** : sa longueur, son rythme, sa forme, son nombre de trajets — et c'est là que se fait la validation d'axe des blocs de catégorie 3 (D2). (c) Un paramètre se change par `Sequential.set_params` (D34), **jamais par affectation directe** sur le bloc. (d) Corollaire normatif qui rend (c) vrai : tout ce que `__post_init__` calcule doit être recalculable par un second appel à `__post_init__`, sans effet de bord — c'est exactement ce que `set_params` fait sur chaque bloc touché. | **La majorité des blocs précalcule au constructeur** : sur les 95 sous-classes de `Processor`, 50 dérivent leur état dans `__post_init__` seul, 7 dans les deux, 10 dans `prepare()` seul — soit **57 sur 67** de ceux qui dérivent quelque chose. Une affectation directe laisse donc la majorité des blocs dans un état incohérent, **en silence** : mesuré sur `Downsampler(4, use_filter=True)`, écrire `.L = 2` laisse le filtre anti-repliement à une coupure de 1/4 alors que la décimation vaut 2, et `set_params(downsampler__L=2)` le reconstruit à 1/2. Ce n'est pas un risque théorique — le défaut corrigé dans les exemples optiques était **le même couple désaccordé**, une coupure et une décimation écrites à la main au lieu d'être dérivées l'une de l'autre, et il a posé un plancher d'ISI à 22,9 dB sous toutes les courbes d'un tutoriel sans rien lever. **La direction est asymétrique**, et c'est ce qui justifie deux endroits plutôt qu'un : mettre de l'état dépendant du signal dans `__post_init__` est **faux** (le constructeur ne connaît pas le signal) ; mettre de l'état purement paramétrique dans `prepare()` est **seulement du travail répété**. `DBP` est dans ce second cas — son `prepare()` ne lit pas le signal du tout — et `FiberLink` y est pour tout sauf `manakov_`. Négligeable devant le SSFM, donc laissé tel quel, mais cela veut dire que « dérive dans `prepare()` » ne prouve pas la dépendance au signal, et la règle doit être écrite plutôt que lue dans le code. | **Geler les blocs** pour que l'affectation lève : fermé par la bibliothèque elle-même — `wiring` (D42c) écrit `reference` sur un bloc **pendant** la passe, et tout estimateur écrit son `theta_` (D23). Un `Processor` doit rester mutable, ce qui distingue sa nature de celle des objets valeur (`FiberSpec`, `Constellation`, `CarrierAllocation`), gelés à juste titre. La discipline ne peut donc pas être portée par le type. **Un `__setattr__` qui reconstruit tout seul** : rendrait `.L = 2` correct, mais fait qu'une affectation déclenche du calcul en silence, exige un drapeau de ré-entrance (`__post_init__` affecte des champs) et un second pour « construction terminée » (la dataclass affecte champ par champ), et surtout **cache la règle au lieu de l'enseigner**. **Tout dériver dans `prepare()`** : supprime l'ambiguïté, mais reconstruirait les 16 001 coefficients d'un SRRC à chaque passe, ou imposerait un cache invalidé à la main — le même problème déplacé d'un cran. **Ne rien écrire** : c'est l'état d'avant, celui qui laisse un bloc neuf dériver sans que rien ne le dise. | **Acté** | + +**Ce que cette décision ne fait pas.** Elle énonce l'invariant, elle ne +le vérifie pas. Le cliquet correspondant — construire un bloc en A, +`set_params` vers B, et exiger la même sortie qu'un bloc construit +directement en B — reste à écrire, et c'est lui qui rendrait un +`__post_init__` incomplet impossible à livrer. Il demande une table de +cas par classe (on ne devine pas des paramètres valides génériquement), +avec les blocs non couverts listés en `skip` explicite pour que +l'omission se voie. + +--- + ### 4.15 Fibre et simulation, séparées (nouveau) | # | Décision | Motif | Alternatives rejetées | Statut | @@ -1138,7 +1155,7 @@ versions de Python annoncées. export JSON et test d'aller-retour. Peu coûteux et immédiatement utile : `summary()` sert le refactor D2, le JSON sert D6 et le papier *(2 après-midi)* 5ter. **D35** — *si et seulement si* le déclencheur est atteint (3ᵉ script - `validation/` au squelette identique) : extraction de `sweep()` depuis + `validation/` au squelette identique) : extraction de `monte_carlo()` depuis ces scripts *(1 après-midi)* ; sinon, reporté post-1.0 6. **D6, D10, D11, D23, D25, D27, D30, D37** — graine de chaîne, docstrings au gabarit §4.10 (validation `numpydoc` en CI dès le @@ -1245,7 +1262,7 @@ ou à corriger. longueur dépend de l'en-tête). À spécifier avec le point précédent. - *Format de scénario sérialisé* — extension naturelle de D31 + D35 : un fichier décrivant chaîne (JSON D31), balayage (`param`, `values`), - graine (D6) et métriques, exécutable par `sweep()` sans code. Utile + graine (D6) et métriques, exécutable par `monte_carlo()` sans code. Utile pour le BANC (rejouer une campagne de mesures) et comme argument de reproductibilité JOSS, mais c'est un format de fichier, pas une API : à spécifier après que D35 existe, jamais avant. Danger identifié à ne diff --git a/CHANGELOG.md b/CHANGELOG.md index ab0c3f3..e4a0e18 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ one release; there is no compatibility layer. | Before (0.91) | After (1.0.0) | |---|---| +| `sweep(chain, param, values, metrics, stimulus)` | `monte_carlo(...)`, same signature — the only bare verb of the root namespace, and a name that said nothing about the points being random draws. Every point reseeds and redraws, so `seed=` is what makes a curve reproducible; the new name says so | | `Serial2Parallel(N, order="F")` → shape `(N_sub, M)` | `Serial2Parallel(N)` → Block layout `(..., T, F)` = `(..., M, N_sub)`, pure C-order reshape; `order` removed | | `Parallel2Serial(order="F")` | `Parallel2Serial()` — C-order flatten of `(..., T, F)`; `order` removed | | OFDM blocks operate on axis 0 by default | OFDM blocks operate on the block content axis -1 (`FFTProcessor`/`IFFTProcessor` hardcode it; `axis` removed there) | @@ -34,7 +35,7 @@ one release; there is no compatibility layer. | `Scope(scope_type="iq", ...)`, `TimeScope`, `SpectrumScope`, `IQScope`, `KDEScope`, `WelchScope` blocks | removed — functions `plot_iq`, `plot_time`, `plot_spectrum`, `plot_kde`, `plot_welch` applied to a tapped signal | | `ofdm.visualizers.FFTMonitor` block | `ofdm.visualizers.plot_subcarrier_amplitude(X, ax=...)` | | `CarrierExtractor(..., pilot_recorder=rec)` | pilot content exposed as the estimated attribute `extractor.pilots_` (D23) | -| `TrainedBased*(target_data=recorder)` | `DataAided*(reference=…)` — the class family is renamed after the standard pair of the field (*data-aided* vs *blind*, which the `Blind*` classes already used), and the known signal an estimator compares against is a `reference`, the same word `sweep(reference=…)` uses. It takes a plain array; when the reference is produced by the chain itself, declare `wiring={"comp.reference": "source"}` | +| `TrainedBased*(target_data=recorder)` | `DataAided*(reference=…)` — the class family is renamed after the standard pair of the field (*data-aided* vs *blind*, which the `Blind*` classes already used), and the known signal an estimator compares against is a `reference`, the same word `monte_carlo(reference=…)` uses. It takes a plain array; when the reference is produced by the chain itself, declare `wiring={"comp.reference": "source"}` | | `DataAidedFIRCompensator(h, reference=…)` | `DataAidedFIRCompensator(reference=…)` — `h` was only an initial value that `fit` overwrote from scratch, i.e. a purely estimated quantity; it is now `h_` and not a constructor parameter | | `Normalizer(gain, method, …)` (inherited from `Amplifier`) | `Normalizer(method, …)` — `gain` is no longer constructible, so `Normalizer('max')` finally means what it reads; the measured gain is `gain_` (D23) | | `Amplifier(gain, axis=…)` | `Amplifier(gain)` — `axis` implemented no defensible model (it scaled only entries at index `axis` of the *last* axis, whatever axis was asked); use `WeightAmplifier` for a per-branch gain | @@ -44,6 +45,36 @@ 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`) | +### Documentation — the AWGN tutorial derives its reference curve + +The Monte-Carlo page compared a measurement against +`constellation.metrics(...)` and left it there, so the closed form +arrived as a number out of an object: the reader saw *that* theory and +simulation agree without being shown *what* the theory is. The page now +states the expression -- the two-PAM decomposition of a square QAM, the +`Q` argument, the reference to Proakis and Salehi section 4.3 -- +evaluates it in plain NumPy, and only then shows that the constellation +returns the same array. The two agree to 3.3e-16, and the reason to +prefer the method afterwards is stated: it is one line rather than five, +and it cannot describe a modulation other than the one the chain +transmits. + +It also no longer stops at one constellation. A closing section sweeps +4-, 16-, 64- and 256-QAM at equal energy per **bit** -- the comparison +`ebn0_to_snr_dB` exists for -- and prints what the density costs: 6.8, +10.5, 14.8 and 19.4 dB of Eb/N0 for a BER of 1e-3. The page says where +the Gray-mapping approximation `P_b ~ P_s / k` stops holding (the +markers run a factor 2.2 above the curve for 256-QAM at 0 dB, and settle +onto it below 1e-2), and where the estimator stops: 100 000 symbols +per point cannot resolve a 4-QAM BER under 5e-6, so the axis stops +there. + +Two defects surfaced while doing it. The page still named `sweep` in two +cross-references after the rename, and the figure it displayed was +written by no script -- `monte_carlo_awgn.py` ended on `plt.show()` and +saved nothing, so `img/monte_carlo_awgn.png` had been orphaned since the +last time someone regenerated it by hand. + ### Added — a chain says when a rate change carries its filter by hand `Sequential` now inspects its module list at construction and warns when @@ -247,6 +278,30 @@ and wiring was not fed, and profiling a wired chain raised. It now shares the edge plan with `forward` and keys its result by block id (`block_ids()`), so two blocks with the same name get one entry each. +### Documentation — one name for the function that builds a chain + +Seven example scripts wrapped their chain in a factory, under four +spellings for the same idea: `get_link`, `link`, `link_chain`, +`get_full_chain`, `uncoded_chain`. A reader who learnt one had to +re-learn it on the next page. + +They now follow one rule, written into the tutorial skill. A factory +exists only when the chain is built more than once — a chain built once +and used once stays inline, because wrapping four blocks in a function +called once with no argument is an indirection between the reader and +their first chain. When it exists it is `get_()`: `get_chain()` +for the page's one chain, `get_channel()` and `get_receiver()` for the +two halves when a page cuts one in two, `get_transmitter()` when what +comes back is a sub-assembly rather than a chain, and +`get_uncoded_chain()` / `get_coded_chain(soft)` when the page compares +two chains that differ in structure rather than in a parameter. + +The third rule is the one that motivated the other two: what the page +varies is a parameter of the function, never a module global read +behind the signature's back. `get_full_chain(n_spans, steps, +linear_only)` also read eleven module-level names, so the call site did +not say what changed between two calls. + ### Documentation — the tutorials on one plan The OFDM, MIMO, Alamouti and PAPR tutorials are rewritten on a single @@ -561,7 +616,7 @@ Substantive changes rather than reorganisation: the two overlaid on one figure (0.25 dB apart on the optimal launch power, 0.17 dB on the peak SNR). -- `sweep(n_jobs=...)`: the points of a sweep are independent by +- `monte_carlo(n_jobs=...)`: the points of a sweep are independent by construction -- each one reconfigures and reseeds the chain from scratch -- so they can run at once. The curve is **identical** value for value, not merely statistically equivalent, because every point @@ -1224,7 +1279,7 @@ Substantive changes rather than reorganisation: to `cupyx.scipy.fft`. CuPy is imported only when a CuPy array is seen and remains a non-dependency. -- `comnumpy.sweep(chain, param, values, metrics, stimulus, ...)` (D35): +- `comnumpy.monte_carlo(chain, param, values, metrics, stimulus, ...)` (D35): the parameter-sweep loop shared by the validation scripts, extracted after the third script needed it. Dotted `set_params` addressing, per-point child seeds, zip semantics for multi-parameter sweeps. @@ -1249,7 +1304,7 @@ Substantive changes rather than reorganisation: `comnumpy.core.metrics.signal_report(x)`, which returns the statistics as a dict for the caller to log, tabulate or assert on. A CI guard test keeps these names out of the public surface. -- `sweep(..., reference="tx")` now names a tapped block and declares the +- `monte_carlo(..., reference="tx")` now names a tapped block and declares the tap itself if needed; blocks no longer hold references to other blocks (`reference=` takes a plain array). - `Sequential(..., wiring={"comp.reference": "ref"})`: declares the diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index aac997b..292b628 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -39,6 +39,21 @@ Run the tests to make sure everything works: make test ``` +### Type checking against the right NumPy + +`pyright` is a blocking check, and its verdict depends on which NumPy is +installed: the 1.26 stubs require `ndarray` to carry type arguments, the +2.x stubs do not. CI lints against the newest NumPy, so an environment +pinned to the oldest supported one reports hundreds of errors that are +not real -- and, worse, can stay silent on one that is. + +Keep a second environment for the check: + +```bash +python -m venv .venv-lint && .venv-lint/bin/pip install -e . +pyright --pythonpath .venv-lint/bin/python +``` + ## Guidelines - **Code style**: Follow PEP 8 formatting. Run `make lint` to check. diff --git a/CONVENTIONS.md b/CONVENTIONS.md index 7a84996..27a1da7 100644 --- a/CONVENTIONS.md +++ b/CONVENTIONS.md @@ -115,11 +115,51 @@ D31: a second input to a block, declared by the chain. A general DAG is still out of scope. **Naming.** The known signal an estimator compares against is the -`reference` — the same word `sweep(reference=...)` uses for the same +`reference` — the same word `monte_carlo(reference=...)` uses for the same idea. Blocks that need one are `DataAided*`, by opposition to the `Blind*` family; that is the standard pair of the field (Proakis, §8: *data-aided* vs *blind* estimation). +## Where derived state lives, and how a parameter changes (D50) + +A block precomputes in two places, and the split is by *what the value +depends on*: + +* **`__post_init__`** — validate, normalize, and precompute whatever + depends **only on the block's own parameters**: filter taps, an + alphabet coerced to an array, an allocation mask. +* **`prepare(x)`** — derive whatever depends on the **signal**: its + length, its rate, its shape, how many paths it carries. Axis + validation for category-3 blocks happens here. + +The direction is asymmetric, which is why there are two places and not +one. Signal-dependent state in `__post_init__` is **wrong** — the +constructor has not seen the signal. Parameter-only state in `prepare()` +is merely **repeated work**; `FiberLink` and `DBP` do this and it is +negligible against the split-step cost, so "derived in `prepare()`" is +not evidence that a value depends on the signal. + +**A parameter changes through `Sequential.set_params` (D34), never by +assignment.** `set_params` re-runs `__post_init__` on every block it +touches; a direct assignment does not, and 57 of the 67 blocks that +derive anything derive it there: + +```python +chain.set_params(downsampler__L=2) # rebuilds the anti-alias filter +chain["downsampler"].L = 2 # does not -- the filter stays at 1/4 +``` + +Blocks cannot be frozen to make that assignment raise: `wiring` writes +`reference` on a block *during* the pass, and every estimator writes its +own `theta_`. Mutability is what separates a `Processor` from a value +object such as `FiberSpec` or `Constellation`. + +**Writing a new block.** Whatever `__post_init__` computes must be +recomputable by calling `__post_init__` again, with no side effect — that +second call is exactly what `set_params` makes. A block that appends to a +list, opens a file or consumes a generator there breaks `set_params` for +everything downstream of it. + ## Error messages Validation failures raise `comnumpy.ShapeError` and follow the template diff --git a/docs/documentation/toplevel/index.rst b/docs/documentation/toplevel/index.rst index 5539a54..0d62aca 100644 --- a/docs/documentation/toplevel/index.rst +++ b/docs/documentation/toplevel/index.rst @@ -11,5 +11,5 @@ import, parameter sweeps, and the plotting style. exceptions serialization - sweep + monte_carlo style diff --git a/docs/documentation/toplevel/sweep.rst b/docs/documentation/toplevel/monte_carlo.rst similarity index 56% rename from docs/documentation/toplevel/sweep.rst rename to docs/documentation/toplevel/monte_carlo.rst index eb89e49..8cd61a0 100644 --- a/docs/documentation/toplevel/sweep.rst +++ b/docs/documentation/toplevel/monte_carlo.rst @@ -1,5 +1,5 @@ Parameter sweeps ================ -.. automodule:: comnumpy.sweep +.. automodule:: comnumpy.monte_carlo :members: diff --git a/docs/getting_started/first_simulation.rst b/docs/getting_started/first_simulation.rst index aa47b73..eabea60 100644 --- a/docs/getting_started/first_simulation.rst +++ b/docs/getting_started/first_simulation.rst @@ -181,7 +181,7 @@ The chain answers the same way: ``seed`` for reproducibility, :meth:`~comnumpy.core.generics.Sequential.set_params` to reconfigure a named block after construction, ``elapsed_`` for the time the last pass took, and -:func:`~comnumpy.sweep.sweep` to run the whole thing over a range of one +:func:`~comnumpy.monte_carlo.sweep` to run the whole thing over a range of one parameter. None of those exist for a script of loose arrays -- not because they would be hard to write, but because there is nothing to attach them to. diff --git a/docs/getting_started/who_is_comnumpy_for.rst b/docs/getting_started/who_is_comnumpy_for.rst index a1570b7..04e417c 100644 --- a/docs/getting_started/who_is_comnumpy_for.rst +++ b/docs/getting_started/who_is_comnumpy_for.rst @@ -55,7 +55,7 @@ the ones that matter in daily research work. A chain can be:: chain.seed(42) # every stochastic block, reproducibly chain.set_params(noise__snr_dB=12) # reconfigured after construction - results = sweep(chain, {"noise.snr_dB": range(0, 20, 2)}) + results = monte_carlo(chain, {"noise.snr_dB": range(0, 20, 2)}) text = to_json(chain) # the experiment, as a file print(chain.to_mermaid()) # the experiment, as a picture diff --git a/docs/tutorials/alamouti.rst b/docs/tutorials/alamouti.rst index afde28b..9976fdd 100644 --- a/docs/tutorials/alamouti.rst +++ b/docs/tutorials/alamouti.rst @@ -234,7 +234,7 @@ Monte Carlo Evaluation ^^^^^^^^^^^^^^^^^^^^^^ Averaging over fading means running the chain once per channel realization, -which is a sweep whose parameter is the channel. :func:`~comnumpy.sweep.sweep` +which is a sweep whose parameter is the channel. :func:`~comnumpy.monte_carlo.sweep` takes several dotted parameter names at once and zips them, so one sweep point sets the channel the signal goes through **and** the channel the detector inverts: diff --git a/docs/tutorials/awgn.rst b/docs/tutorials/awgn.rst index 3883929..09c65fb 100644 --- a/docs/tutorials/awgn.rst +++ b/docs/tutorials/awgn.rst @@ -7,9 +7,9 @@ that repetition is what "Monte Carlo simulation" means, and there is nothing more to it than a loop. We will write that loop by hand first, because it is worth seeing exactly -what a sweep is made of. Then we will replace it with :func:`~comnumpy.sweep.sweep`, -which does the same four things in one call, and use *that* in every tutorial -that follows. +what a sweep is made of. Then we will replace it with +:func:`~comnumpy.monte_carlo.monte_carlo`, which does the same four things in +one call, and use *that* in every tutorial that follows. .. note:: @@ -22,10 +22,12 @@ that follows. - How to reconfigure a chain between runs with ``set_params``, and make each run reproducible with ``seed``. -- How to write a Monte Carlo loop, and then how to replace it with ``sweep``. -- How to compare a measured curve with the closed form it should follow. +- How to write a Monte Carlo loop, and then how to replace it with ``monte_carlo``. +- How to write the closed form a measured curve should follow, evaluate it in + plain NumPy, and then get the same numbers from the constellation itself. - How to draw the result with ``plot_error_rate``, the figure this library uses for every error rate. +- How to compare modulation orders fairly, at equal energy per bit. Introduction @@ -39,6 +41,7 @@ Make sure you have the following Python libraries installed: .. code:: numpy + scipy matplotlib comnumpy @@ -49,7 +52,7 @@ We start by importing the necessary libraries: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 1-10 + :lines: 1-14 Define Parameters @@ -62,7 +65,7 @@ closed form we compare against below. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 13-16 + :lines: 16-19 AWGN Communication Chain @@ -73,26 +76,27 @@ Define Chain We define the communication chain using the ``Sequential`` object. The chain includes symbol generation, mapping, transmission over an AWGN channel, -and symbol demapping. +and symbol demapping. It is written as a function because the last section of +this page builds one per constellation order. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 17-24 + :lines: 22-39 The processors are: -- ``SymbolGenerator`` +- ``SymbolGenerator`` Generates a stream of integer-valued symbols to transmit. It is named ``"tx"`` so the sweep can compare the chain output against the transmitted symbols (see ``reference="tx"`` below). -- ``SymbolMapper`` +- ``SymbolMapper`` Maps integers to QAM constellation points. -- ``AWGN`` +- ``AWGN`` Simulates the effect of noise for a given SNR value (here expressed in dB). -- ``SymbolDemapper`` +- ``SymbolDemapper`` Maps received noisy constellation points back to integers. The chain, as the chain itself describes it: @@ -106,7 +110,7 @@ outline marks a tapped block: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 63-69 + :lines: 137-143 Monte Carlo Simulation """""""""""""""""""""" @@ -118,7 +122,7 @@ against the transmitted symbols. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 26-34 + :lines: 41-48 Three chain services appear there, and they are the ones every study is made of. @@ -130,43 +134,99 @@ construction, with the dotted notation ``"awgn_channel.snr_dB"``; this is why blocks are named. And the tap returns what the transmitter produced, so the metric has something to compare against -- after the run, never before. -Now the same thing in one call. :func:`~comnumpy.sweep.sweep` takes the chain, -the dotted name of what varies, the values it takes, and the metrics to -collect -- and does the reseed, the reconfigure, the run and the measurement +Now the same thing in one call. :func:`~comnumpy.monte_carlo.monte_carlo` takes +the chain, the dotted name of what varies, the values it takes, and the metrics +to collect -- and does the reseed, the reconfigure, the run and the measurement at every point: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 36-48 + :lines: 50-62 .. code:: - loop : 7.406e-01 3.535e-01 7.237e-03 - sweep : 7.410e-01 3.533e-01 7.117e-03 + loop : 7.406e-01 3.535e-01 7.237e-03 + monte_carlo : 7.410e-01 3.533e-01 7.117e-03 The two are the same computation. They do not print the same digits because -``sweep`` gives each point its own child seed rather than reseeding every +``monte_carlo`` gives each point its own child seed rather than reseeding every point to the same value, so the noise realizations differ; the gap is the Monte Carlo error of a million symbols, not a difference of method. -**From here on, the other tutorials use** ``sweep`` **without rewriting the +**From here on, the other tutorials use** ``monte_carlo`` **without rewriting the loop.** When you see it sweeping a channel matrix rather than an SNR (in :doc:`mimo` and :doc:`alamouti`), it is still these four steps -- only the parameter that varies has changed. Theoretical SER +^^^^^^^^^^^^^^^ + +A measured curve is worth little on its own. What makes it a result is the +curve it is supposed to lie on. + +The closed form """"""""""""""" -For comparison, we ask the constellation for its own closed form. It -knows its family and its order, so the theory cannot describe a -modulation other than the one the chain transmits; ``per="symbol"`` says -that the swept SNR is the symbol SNR ``AWGN(snr_dB=)`` takes, and the -conversion to :math:`E_b/N_0` happens inside. +A square :math:`M`-QAM constellation is the Cartesian product of two +independent :math:`\sqrt{M}`-PAM constellations, one on each quadrature. With +:math:`k = \log_2 M` bits per symbol, the signal-to-noise ratio per bit +:math:`\gamma_b = E_b/N_0` and the Gaussian tail function + +.. math:: + + Q(x) = \frac{1}{\sqrt{2\pi}} \int_x^{\infty} e^{-t^2/2}\, dt + +the error probability of one PAM component is + +.. math:: + + P_{\sqrt{M}} = 2 \left(1 - \frac{1}{\sqrt{M}}\right) + Q\!\left(\sqrt{\frac{3 k \gamma_b}{M-1}}\right) + +A symbol is received correctly only when **both** components are, so the +symbol error probability is + +.. math:: + + P_s = 1 - \left(1 - P_{\sqrt{M}}\right)^2 + +This is the standard result; see Proakis and Salehi [Proakis2008]_, Section 4.3. +Note which SNR parameterizes it: the expression is written against +:math:`\gamma_b`, the ratio per *bit*, whereas ``AWGN(snr_dB=)`` takes a symbol +SNR :math:`\gamma_s = k\gamma_b`. Getting that wrong moves the curve by +:math:`10\log_{10} k` dB -- 6 dB for 16-QAM -- with nothing to signal it. + +In NumPy, written out +""""""""""""""""""""" + +There is no machinery in this. ``scipy.special.erfc`` gives :math:`Q`, and the +two expressions above transcribe line for line: .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 54-54 + :lines: 65-78 + +From the constellation +"""""""""""""""""""""" + +The same expression is already in the library, and the constellation knows +which one applies to it: it carries its family and its order, so the theory +cannot end up describing a modulation other than the one the chain transmits. +``per="symbol"`` says that the swept SNR is the symbol SNR, and the division +by :math:`k` above happens inside. + +.. literalinclude:: ../../examples/simple/monte_carlo_awgn.py + :language: python + :lines: 80-85 + +.. code:: + + largest gap between the two closed forms: 3.331e-16 + +Machine precision -- they are the same formula. From here on we call the +method, which is one line instead of five and cannot describe the wrong +modulation. Results and Visualization @@ -180,27 +240,109 @@ share a colour, so a pair reads as one statement. .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py :language: python - :lines: 56-62 + :lines: 87-93 .. image:: img/monte_carlo_awgn.png :width: 100% :align: center +The markers sit on the line over five decades. The last point, at 21 dB, is an +SER of :math:`10^{-6}` measured from :math:`10^6` symbols -- one error. It is +on the curve by luck as much as by physics, and it is the point at which this +sweep has run out of samples rather than out of accuracy. + + +Comparing modulation orders +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +One constellation against its own theory says the simulation is right. It does +not say what a constellation is *worth*. For that, several orders have to be +put on the same axes -- and the choice of that axis is the whole question. + +At equal :math:`E_s/N_0`, a dense constellation is unduly favoured: it spends +the same energy per symbol while carrying more bits. The fair comparison is at +equal energy per **bit**, :math:`E_b/N_0`, which is why +:func:`~comnumpy.core.utils.ebn0_to_snr_dB` exists. The conversion needs +:math:`k` -- chain-level knowledge the channel block does not have and should +not guess (decision D41), so it happens at the call site: + +.. literalinclude:: ../../examples/simple/monte_carlo_awgn.py + :language: python + :lines: 95-123 + +The metric changes with the axis. Against :math:`E_b/N_0` the natural quantity +is the bit error rate, and ``compute_ber`` needs the symbol width, which +``functools.partial`` supplies once per order. + +The reference curve is no longer exact here, and the figure shows it. The BER +of the closed form is the Gray-mapping approximation +:math:`P_b \simeq P_s / k`, which assumes a symbol error corrupts exactly one +bit out of :math:`k` -- true when the wrong decision is a nearest neighbour, +false when the noise is large enough to reach further. So the markers sit +**above** the line at the left of the figure, by 17 % for 16-QAM and by a +factor 2.2 for 256-QAM at 0 dB, and settle onto it as the error rate falls: +below :math:`10^{-2}` the two agree to within a few percent for every order. +This is an approximation behaving as documented, not a simulation +disagreeing with theory. + +.. code:: + + order bits/symbol Eb/N0 at BER=1e-3 + 4-QAM 2 6.8 dB + 16-QAM 4 10.5 dB + 64-QAM 6 14.8 dB + 256-QAM 8 19.4 dB + +.. literalinclude:: ../../examples/simple/monte_carlo_awgn.py + :language: python + :lines: 125-135 + +.. image:: img/monte_carlo_awgn_orders.png + :width: 100% + :align: center + +Going from 2 to 8 bits per symbol costs **12.6 dB** of energy per bit at +:math:`10^{-3}`, spent in three roughly equal steps -- 3.7, 4.3 then 4.6 dB +per doubling of :math:`k`, the price rising slowly with the order. That is +what the rest of +a communication system is built to negotiate: coding buys some of it back +(:doc:`coding`), shaping buys a little more (:doc:`shaping`), and a channel +that is not simply AWGN changes the terms entirely (:doc:`ofdm`). + +The measured curves stop where the estimator does. Each point here is 100 000 +symbols, so a 4-QAM run has 200 000 bits and cannot resolve a BER below about +:math:`5 \times 10^{-6}`; ``plot_error_rate`` drops the points where no error +was seen rather than drawing them at zero, and the axis stops at +:math:`10^{-6}` for the same reason. Reading a curve past its own floor is the +most common way to publish a wrong figure. + Conclusion ^^^^^^^^^^ -You have turned one run into a curve, and you have seen that ``sweep`` is not -a new concept but the loop you already wrote, packaged. +You have turned one run into a curve, checked it against a closed form you +wrote yourself, and used it to price four modulation orders against each +other. You have learned how to: - Reconfigure and reseed a chain between runs. - Collect a metric over a range of parameter values, by hand and with - ``sweep``. -- Read a measured curve against a closed form. + ``monte_carlo``. +- Write a theoretical error rate in NumPy, and get the same numbers from the + constellation. +- Read a measured curve against a closed form, and know where it stops being + readable. +- Compare modulations at equal energy per bit. From here, :doc:`ofdm` keeps the chain and changes the channel: instead of one noise term, the signal arrives several times, and the receiver has to undo that. Two ways of doing so, at very different prices. + +References +^^^^^^^^^^ + +.. [Proakis2008] J. G. Proakis and M. Salehi, *Digital Communications*, + 5th ed., McGraw-Hill, 2008. Section 4.3, symbol and bit error + probabilities for PAM, PSK and QAM over the AWGN channel. diff --git a/docs/tutorials/coding.rst b/docs/tutorials/coding.rst index 2c3d52d..b7128ff 100644 --- a/docs/tutorials/coding.rst +++ b/docs/tutorials/coding.rst @@ -17,7 +17,7 @@ tutorial is plotted against :math:`E_b/N_0` rather than SNR. .. note:: - **Before you start.** :doc:`awgn` introduced ``sweep`` and + **Before you start.** :doc:`awgn` introduced ``monte_carlo`` and ``plot_error_rate``, both used here without comment. Nothing else is assumed: the codes are built from scratch. @@ -318,7 +318,7 @@ independent, so they can be run at the same time: .. code:: python - results = sweep(chain, "noise.snr_dB", snr_dB(ebn0_dB, rate), + results = monte_carlo(chain, "noise.snr_dB", snr_dB(ebn0_dB, rate), {"ber": compute_ser}, n_bits, reference="tx", seed=4, n_jobs=4) diff --git a/docs/tutorials/gn_model.rst b/docs/tutorials/gn_model.rst index f26b14d..069dcce 100644 --- a/docs/tutorials/gn_model.rst +++ b/docs/tutorials/gn_model.rst @@ -107,27 +107,29 @@ of this page asks it several. the GN model's convention, and the first of three places a factor of two hides in this tutorial. -``analytic_ase`` is the second. It calls -:func:`~comnumpy.optical.utils.compute_erbium_doped_fiber_N_ase`, which is -*the same function* :class:`~comnumpy.optical.links.FiberLink` uses to -generate the noise it adds during propagation. That is deliberate: in Part 2 -this prediction is compared against a measurement through the whole chain, -and had the prediction used a formula retyped from a textbook, the comparison -would only be checking the transcription. Using the library's own function -makes the comparison test the *chain* instead -- the spectral density, the -bandwidth, the two polarizations, the matched filter. - -The factor two in ``analytic_ase`` is the polarization: the function returns -the density for one, and an amplifier emits into both, each with its own -independent noise. +``analytic_ase`` is the second, and it does not compute anything: it asks +the link. :meth:`~comnumpy.optical.links.FiberLink.budget` returns the noise +that same block will add during propagation, integrated over the bandwidth a +receiver keeps. That is deliberate: in Part 2 this prediction is compared +against a measurement through the whole chain, and had the prediction used a +formula retyped from a textbook, the comparison would only be checking the +transcription. Asking the link makes the comparison test the *chain* instead +-- the spectral density, the bandwidth, the two polarizations, the matched +filter. + +``polarizations=2`` is where the factor lives. A spectral density is per +polarization, and an amplifier emits into both, each with its own +independent noise; writing that product by hand at every call site is how +two pages end up quoting budgets 3 dB apart. .. warning:: **The 3 dB that becomes 9 dB.** ``power_W`` is the power of the *channel*, summed over both polarizations, because that is the convention the GN model uses. Each polarization therefore carries half of it, which - is what ``launch_gain`` is for in Part 2. Give each polarization the full - ``power_W`` instead and you launch 3 dB more than you think; since the + is what ``launch_amplitude(..., polarizations=2)`` is for in Part 2. Give + each polarization the full ``power_W`` instead and you launch 3 dB more + than you think; since the nonlinear interference goes as the cube of the power, your simulation then reports **9 dB** more of it than the model predicts, and the model looks broken when it is not. @@ -162,7 +164,7 @@ roll-off smearing the very guard band the figure exists to show. eta = 1223 /W^2 (GN model, one channel) P_ASE = -20.88 dBm (5 spans, NF = 6 dB, both polarizations) - optimum = +1.74 dBm SNR = 20.86 dB + optimum = +1.74 dBm SNR = 20.86 dB check : eta P^3 / (P_ASE/2) = 1.000000 Two things read off the figure. The **cut** is the filled channel: the one @@ -248,17 +250,18 @@ Everything above is a prediction. Now we propagate samples and find out. .. literalinclude:: ../../examples/optical/gn_model.py :language: python - :lines: 115-153 + :lines: 114-149 .. mermaid:: mermaid/gn_model.mmd -``launch_gain`` is the factor of two of the warning above, made explicit: -:math:`\sqrt{P/2}` per polarization for a channel power :math:`P`. +``launch_amplitude(power_W, polarizations=2)`` is the factor of two of the +warning above, named rather than written: :math:`\sqrt{P/2}` of field per +polarization for a channel power :math:`P`. ``name="launch"`` and ``name="fibre"`` are what let the rest of the page *reconfigure* the link instead of rebuilding it: ``chain.set_params`` reaches a block by its identifier and re-runs its precomputation (decision D34), and -:func:`~comnumpy.sweep` drives ``"launch.gain"`` over a list of values the +:func:`~comnumpy.monte_carlo` drives ``"launch.gain"`` over a list of values the same way. One chain, built once, answers every question below. ``taps=["tx"]`` marks the mapper output as readable from outside the chain -- @@ -293,7 +296,7 @@ form. At 0 dBm launch the measured SNR *is* :math:`-P_{\mathrm{ASE}}` in dBm. .. literalinclude:: ../../examples/optical/gn_model.py :language: python - :lines: 155-165 + :lines: 151-161 .. code:: @@ -320,7 +323,7 @@ Now the nonlinear term, at fourteen launch powers. .. literalinclude:: ../../examples/optical/gn_model.py :language: python - :lines: 166-193 + :lines: 162-189 .. code:: @@ -376,7 +379,7 @@ modulation formats are not, and the model is blind to the difference. .. literalinclude:: ../../examples/optical/gn_model.py :language: python - :lines: 195-204 + :lines: 191-200 .. code:: @@ -401,7 +404,7 @@ decibel for 16QAM, less as the format grows. .. literalinclude:: ../../examples/optical/gn_model.py :language: python - :lines: 206-210 + :lines: 202-206 Going further diff --git a/docs/tutorials/img/monte_carlo_awgn.png b/docs/tutorials/img/monte_carlo_awgn.png index 38ee84a..2a9c745 100644 Binary files a/docs/tutorials/img/monte_carlo_awgn.png and b/docs/tutorials/img/monte_carlo_awgn.png differ diff --git a/docs/tutorials/img/monte_carlo_awgn_orders.png b/docs/tutorials/img/monte_carlo_awgn_orders.png new file mode 100644 index 0000000..a9a4584 Binary files /dev/null and b/docs/tutorials/img/monte_carlo_awgn_orders.png differ diff --git a/docs/tutorials/img/nli_simulation_fig1.png b/docs/tutorials/img/nli_simulation_fig1.png index 2cf3699..a22c9ce 100644 Binary files a/docs/tutorials/img/nli_simulation_fig1.png and b/docs/tutorials/img/nli_simulation_fig1.png differ diff --git a/docs/tutorials/img/nli_simulation_fig2.png b/docs/tutorials/img/nli_simulation_fig2.png index 6c1f283..a9ea90f 100644 Binary files a/docs/tutorials/img/nli_simulation_fig2.png and b/docs/tutorials/img/nli_simulation_fig2.png differ diff --git a/docs/tutorials/img/one_shot_nli_fig1.png b/docs/tutorials/img/one_shot_nli_fig1.png index 93ae3ee..f2b8e58 100644 Binary files a/docs/tutorials/img/one_shot_nli_fig1.png and b/docs/tutorials/img/one_shot_nli_fig1.png differ diff --git a/docs/tutorials/img/one_shot_nli_fig2.png b/docs/tutorials/img/one_shot_nli_fig2.png index 408fdf0..95dd2cd 100644 Binary files a/docs/tutorials/img/one_shot_nli_fig2.png and b/docs/tutorials/img/one_shot_nli_fig2.png differ diff --git a/docs/tutorials/img/one_shot_nli_fig3.png b/docs/tutorials/img/one_shot_nli_fig3.png index e791909..a9595bd 100644 Binary files a/docs/tutorials/img/one_shot_nli_fig3.png and b/docs/tutorials/img/one_shot_nli_fig3.png differ diff --git a/docs/tutorials/img/one_shot_nli_fig4.png b/docs/tutorials/img/one_shot_nli_fig4.png index 1e9b103..28ee3dd 100644 Binary files a/docs/tutorials/img/one_shot_nli_fig4.png and b/docs/tutorials/img/one_shot_nli_fig4.png differ diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 89a9b6f..7115774 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -19,8 +19,9 @@ and nothing is explained twice. * - 1 - :doc:`awgn` - A Monte Carlo study: first as an ordinary ``for`` loop, then the same - thing with ``sweep``. Also ``seed``, ``set_params`` and - ``plot_error_rate``. + thing with ``monte_carlo``, read against a closed form written out in + NumPy, and ending on four QAM orders at equal energy per bit. Also + ``seed``, ``set_params`` and ``plot_error_rate``. * - 2 - :doc:`ofdm` - A frequency-selective channel, and the two ways to equalize it: diff --git a/docs/tutorials/mimo.rst b/docs/tutorials/mimo.rst index 2463011..9d19ed7 100644 --- a/docs/tutorials/mimo.rst +++ b/docs/tutorials/mimo.rst @@ -7,7 +7,7 @@ compare on it. .. note:: - **Before you start.** :doc:`awgn` introduced ``sweep``, which is used + **Before you start.** :doc:`awgn` introduced ``monte_carlo``, which is used here to average over channel realizations rather than over noise. The chain is the same object; only the number of antennas changes. @@ -227,7 +227,7 @@ Monte Carlo Evaluation A single channel realization proves nothing: over fading, the error rate is an *average*, dominated by the rare draws where the matrix is nearly singular. Averaging means running the chain once per realization, which is a -sweep whose parameter is the channel. :func:`~comnumpy.sweep.sweep` takes +sweep whose parameter is the channel. :func:`~comnumpy.monte_carlo.sweep` takes several dotted parameter names at once and zips them, so one sweep point sets the matrix the signal goes through **and** the one the detector inverts: diff --git a/docs/tutorials/ofdm.rst b/docs/tutorials/ofdm.rst index af19ece..5392e5a 100644 --- a/docs/tutorials/ofdm.rst +++ b/docs/tutorials/ofdm.rst @@ -14,7 +14,7 @@ channel (**OFDM**). .. note:: **Before you start.** :doc:`awgn` introduced the chain, ``seed``, - ``sweep`` and ``plot_error_rate``; they are used here without being + ``monte_carlo`` and ``plot_error_rate``; they are used here without being re-explained. What is new is the channel: it is no longer a single coefficient. diff --git a/docs/tutorials/optical_fiber_nonlinearity.rst b/docs/tutorials/optical_fiber_nonlinearity.rst index b9366b7..113a0f0 100644 --- a/docs/tutorials/optical_fiber_nonlinearity.rst +++ b/docs/tutorials/optical_fiber_nonlinearity.rst @@ -71,16 +71,36 @@ of variance :math:`N_{ase} = (G-1)h\nu\,n_{sp}`, with Implementation """""""""""""" -We reproduce the setup of Häger and Pfister: 25 spans of 80 km, 16-QAM at -10.7 GBd, root-raised-cosine pulses, and 200 split steps per span, which is -enough that the propagation itself is not in question -- 500 steps give the -same effective SNR to three decimals. +We follow the setup of Häger and Pfister: 25 spans of 80 km, 16-QAM at +10.7 GBd, root-raised-cosine pulses, and 50 split steps per span. Their own +reproduction, with a Gaussian stimulus and the paper's 500 steps, is in +``validation/optical_dbp_hager.py``; this page carries a constellation so +it can also show an error rate. + +Fifty steps is not a guess. The step-size error is worst where the +nonlinearity is strongest, so it is measured there -- at the highest launch +power of the sweep on the next page -- by refining until the answer stops +moving: + +.. code:: + + StPS = 200 SNR = 11.758 dB +0.000 dB 38.0 s + StPS = 100 SNR = 11.759 dB +0.001 dB 19.3 s + StPS = 50 SNR = 11.763 dB +0.005 dB 9.4 s + StPS = 25 SNR = 11.776 dB +0.018 dB 5.0 s + StPS = 12 SNR = 11.845 dB +0.087 dB 2.5 s + +Fifty steps sit **0.005 dB** from the converged answer and cost a quarter of +what two hundred cost. Every number this page reports is quoted to a +hundredth of a decibel at best, so the discretization is two orders of +magnitude below what is claimed from it -- which is the only sense in which +a step count is ever "enough". The whole system is one chain, and the number of spans is an argument: .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 1-89 + :lines: 1-88 Two things in that chain are worth naming. @@ -109,17 +129,17 @@ be read against. .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 92-128 + :lines: 91-127 .. code:: spans measured ASE only the fibre SER phase time - 1 36.09 dB 37.11 dB 1.02 dB 0.0000 -1.2 deg 1.1 s - 5 25.93 dB 30.65 dB 4.72 dB 0.0000 -7.0 deg 4.0 s - 10 20.99 dB 27.68 dB 6.69 dB 0.0007 -14.8 deg 7.4 s - 15 18.55 dB 25.94 dB 7.38 dB 0.0023 -22.9 deg 14.2 s - 20 16.73 dB 24.72 dB 7.99 dB 0.0091 -31.0 deg 14.9 s - 25 15.30 dB 23.70 dB 8.39 dB 0.0238 -39.2 deg 18.5 s + 1 36.09 dB 37.11 dB 1.02 dB 0.0000 -1.2 deg 0.1 s + 5 25.93 dB 30.65 dB 4.72 dB 0.0000 -7.0 deg 0.4 s + 10 20.99 dB 27.68 dB 6.69 dB 0.0007 -14.8 deg 0.7 s + 15 18.56 dB 25.94 dB 7.38 dB 0.0023 -22.9 deg 1.0 s + 20 16.73 dB 24.72 dB 7.99 dB 0.0091 -31.0 deg 1.4 s + 25 15.31 dB 23.70 dB 8.39 dB 0.0238 -39.2 deg 1.7 s The third column is the subject of this tutorial. It is the price of the Kerr effect, measured rather than argued: **1.02 dB after one span, 8.39 dB @@ -138,7 +158,7 @@ are cheap enough to leave in the script: .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 130-155 + :lines: 129-154 .. code:: @@ -170,7 +190,7 @@ two noises that close add. The check reports it rather than hiding it. .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 157-159 + :lines: 156-158 .. image:: img/one_shot_nli_fig1.png :width: 100% @@ -184,7 +204,7 @@ three distances: .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 161-181 + :lines: 160-180 .. image:: img/one_shot_nli_fig2.png :width: 100% @@ -210,31 +230,32 @@ span, 39 degrees over the link. Left in, that rotation would dominate the error rate, which is why the compensator sits in the chain and not in a comment. -The last column is the one to keep. It grows with the number of spans, and -``profile_execution_time`` says why -- it runs the chain and times each -block on the way through: +The last column is the one to keep, and it is not the receiver's: +``profile_execution_time`` runs the chain and times each block on the way +through, which says where it all went. .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 183-188 + :lines: 182-187 .. code:: block time data_tx 0.1 ms signal_tx 0.0 ms - upsampler 0.3 ms - srrcfilter 2.6 ms - signal_amplifier 0.2 ms - link 17961.4 ms - rx_field 1.7 ms - dbp 13.5 ms - srrcfilter_2 0.7 ms + upsampler 0.2 ms + srrcfilter 1.9 ms + signal_amplifier 0.1 ms + link 1721.7 ms + rx_field 0.9 ms + dbp 3.9 ms + srrcfilter_2 0.5 ms downsampler 0.0 ms signal_amplifier_2 0.0 ms - phase 0.2 ms + phase 0.1 ms + data_rx 0.4 ms -Twelve blocks, and one of them is **99.9 %** of the run. The split-step +Thirteen blocks, and one of them is **99.5 %** of the run. The split-step propagation is 25 spans of 200 steps, each an FFT pair and a pointwise rotation; everything else is a handful of milliseconds. Keep that ratio in mind -- it is what the Monte-Carlo section below has to work around. @@ -261,19 +282,19 @@ removed by the backward pass. Dispersion compensation alone is the same algorithm with the nonlinear term switched off -- one step per span instead of :math:`\mathrm{StPS}` -- so the -two strategies are one argument of ``get_full_chain`` apart. +two strategies are one argument of ``get_chain`` apart. Results """"""" .. literalinclude:: ../../examples/optical/one_shot_NLI.py :language: python - :lines: 190-211 + :lines: 189-210 .. code:: - dispersion compensation SNR=15.30 dB receiver 13.5 ms - digital back-propagation SNR=23.48 dB SER=0.0000 receiver 3037.2 ms residual phase=-0.1 deg + dispersion compensation SNR=15.31 dB receiver 3.9 ms + digital back-propagation SNR=23.48 dB SER=0.0000 receiver 1045.1 ms residual phase=-0.1 deg .. image:: img/one_shot_nli_fig4.png :width: 100% @@ -307,8 +328,8 @@ What it costs The last column of the table is the reason DBP is not simply switched on everywhere. Dispersion compensation is one FFT pair for the whole link; DBP at :math:`\mathrm{StPS}` steps per span is :math:`N_{sp} \times \mathrm{StPS}` -FFT pairs plus as many pointwise phase rotations. Here that is 13.5 ms -against 3037 ms -- **225 times** -- for 8.2 dB. +FFT pairs plus as many pointwise phase rotations. Here that is 3.9 ms +against 1045 ms -- **268 times** -- for 8.2 dB. That ratio is what the literature on low-complexity back-propagation exists to improve, and it is also why the useful question is not "DBP or not" but *how @@ -325,7 +346,7 @@ the amplifier noise alone would allow. .. literalinclude:: ../../examples/optical/NLI_simulation.py :language: python - :lines: 1-84 + :lines: 1-86 Before propagating anything, the closed form of :doc:`gn_model` says where the optimum should fall. It describes this link with two changes from the @@ -334,7 +355,7 @@ instead of two, which is what ``polarizations=1`` is for: .. literalinclude:: ../../examples/optical/NLI_simulation.py :language: python - :lines: 87-104 + :lines: 89-105 .. code:: @@ -344,7 +365,7 @@ That is the prediction the sweep below has to land on. .. literalinclude:: ../../examples/optical/NLI_simulation.py :language: python - :lines: 107-198 + :lines: 108-219 .. code:: @@ -377,12 +398,12 @@ and back-propagation **moves it to the right**: .. code:: receiver best SNR at power total time - amplifier noise only 26.35 dB 4.5 dBm 0.8 s - dispersion compensation 18.78 dB -1.5 dBm 0.7 s - DBP, 1 step/span 19.28 dB 0.0 dBm 1.7 s - DBP, 2 steps/span 20.64 dB 0.0 dBm 3.4 s - DBP, 4 steps/span 23.81 dB 3.0 dBm 6.5 s - DBP, 50 steps/span 25.94 dB 4.5 dBm 76.8 s + amplifier noise only 26.35 dB 4.5 dBm 0.4 s + dispersion compensation 18.78 dB -1.5 dBm 0.4 s + DBP, 1 step/span 19.28 dB 0.0 dBm 0.7 s + DBP, 2 steps/span 20.64 dB 0.0 dBm 1.4 s + DBP, 4 steps/span 23.82 dB 3.0 dBm 2.6 s + DBP, 50 steps/span 25.94 dB 4.5 dBm 30.7 s Start with the second row, because it is the one the closed form claims to predict: @@ -425,11 +446,25 @@ adds the shape the SNR curve hides: a bathtub, whose left wall is the amplifier noise and whose right wall is the nonlinearity. Back-propagation does not move the left wall -- it cannot -- and pushes the right one out. -The flat bottom at :math:`1.2 \times 10^{-4}` is the estimator, not the link: -8192 symbols per point cannot resolve fewer than one error, so the two -best receivers are floor-limited from 0 dBm on. That is why the comparison -above is made in effective SNR, which uses every symbol rather than only the -wrong ones. +The bottom of the figure is the estimator, not the link, and the counters +say so in errors rather than in a caveat: + +.. code:: + + receiver errors at its best power + amplifier noise only 0 over 49152 symbols + dispersion compensation 15 over 49152 symbols + DBP, 1 step/span 5 over 49152 symbols + DBP, 2 steps/span 0 over 49152 symbols + DBP, 4 steps/span 0 over 49152 symbols + DBP, 50 steps/span 0 over 49152 symbols + +Four receivers of six saw **no error at all** at their own best power. +Those points are not measurements of an error rate; they are the +statement that it is below one in 49152, i.e. :math:`2 \times 10^{-5}`. +That is why the comparison above is made in effective SNR, which uses +every symbol rather than only the wrong ones -- and why a curve that +stops has to say how many errors it stopped on. Conclusion diff --git a/examples/mimo/monte_carlo_simulation_1.py b/examples/mimo/monte_carlo_simulation_1.py index 52f187d..d342bf8 100644 --- a/examples/mimo/monte_carlo_simulation_1.py +++ b/examples/mimo/monte_carlo_simulation_1.py @@ -44,13 +44,13 @@ for index_snr, snr_dB in enumerate(snr_dB_list): sigma2 = sig_power*(10**(-snr_dB/10)) - chain["noise"].sigma2 = sigma2 + chain.set_params(noise__sigma2=sigma2) for _ in range(N_test): # new channel realization H = rayleigh_channel(N_r=N_r, N_t=N_t) - chain["channel"].H = H + chain.set_params(channel__H=H) # generate data Y = chain((N_t, N)) diff --git a/examples/mimo/monte_carlo_simulation_2.py b/examples/mimo/monte_carlo_simulation_2.py index fa6303f..1f19162 100644 --- a/examples/mimo/monte_carlo_simulation_2.py +++ b/examples/mimo/monte_carlo_simulation_2.py @@ -4,7 +4,7 @@ from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolMapper -from comnumpy.core.metrics import compute_ber +from comnumpy.core.metrics import ErrorCounter from comnumpy.core.utils import Constellation from comnumpy.mimo.channels import AWGN, FlatMIMOChannel from comnumpy.mimo.utils import rayleigh_channel @@ -35,18 +35,23 @@ # compute snr list from snr per bit in dB snr_dB_list = np.arange(0, 45, 5) -# perform simulation -ber_data = np.zeros((len(snr_dB_list), len(detector_names))) +# perform simulation. One counter per (SNR, detector): it accumulates +# errors and bits over the trials, so the rate is a ratio of totals and +# the count that makes it credible stays available. +counters = {} +for index_snr in range(len(snr_dB_list)): + for index in range(len(detector_names)): + counters[index_snr, index] = ErrorCounter(width=int(np.log2(M))) for index_snr, snr_dB in enumerate(snr_dB_list): sigma2 = N_t*(10**(-snr_dB/10)) - chain["noise"].sigma2 = sigma2 + chain.set_params(noise__sigma2=sigma2) for _ in range(N_test): # new channel realization H = rayleigh_channel(N_r=N_r, N_t=N_t) - chain["channel"].H = H + chain.set_params(channel__H=H) # generate data Y = chain((N_t, N)) @@ -59,9 +64,15 @@ # perform detection S_est = detector(Y) # evaluate metrics - ber_data[index_snr, index] += compute_ber(S_ref, S_est, width=int(np.log2(M))) + counters[index_snr, index].update(S_ref, S_est) + +ber_data = np.zeros((len(snr_dB_list), len(detector_names))) +for (index_snr, index), counter in counters.items(): + ber_data[index_snr, index] = counter.rate - ber_data[index_snr, :] /= N_test +worst = min(counters.values(), key=lambda counter: counter.n_errors) +print(f"fewest errors at any point: {worst.n_errors} over " + f"{worst.n_symbols * worst.width} bits") # plot figures for index, detector_name in enumerate(detector_names): diff --git a/examples/mimo/one_shot_alamouti.py b/examples/mimo/one_shot_alamouti.py index e06a945..2c51ab1 100644 --- a/examples/mimo/one_shot_alamouti.py +++ b/examples/mimo/one_shot_alamouti.py @@ -6,7 +6,7 @@ import matplotlib.pyplot as plt import numpy as np -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolDemapper, SymbolMapper @@ -55,7 +55,7 @@ # --- the three links -------------------------------------------------- -def get_link(kind, H, sigma2=sigma2): +def get_chain(kind, H, sigma2=sigma2): """Build the chain of one scheme, on a given channel realization. The three schemes differ by two blocks -- what is put on the @@ -84,9 +84,9 @@ def get_link(kind, H, sigma2=sigma2): raise ValueError(f"unknown link {kind!r}") -alamouti = get_link("alamouti", rayleigh_channel(1, 2, seed=42)) -siso = get_link("linear", rayleigh_channel(1, 1, seed=1)) -mrc = get_link("linear", rayleigh_channel(2, 1, seed=2)) +alamouti = get_chain("alamouti", rayleigh_channel(1, 2, seed=42)) +siso = get_chain("linear", rayleigh_channel(1, 1, seed=1)) +mrc = get_chain("linear", rayleigh_channel(2, 1, seed=2)) # --- one shot --------------------------------------------------------- alamouti.seed(7) # every stochastic block, reproducibly @@ -110,7 +110,7 @@ def get_link(kind, H, sigma2=sigma2): def average_ser(chain, n_rx, n_tx, snr_dB, stimulus, n_channels, seed=0): """Average one chain over independent quasi-static fading draws. - The chain is built once. ``sweep`` takes several dotted parameters + The chain is built once. ``monte_carlo`` takes several dotted parameters at once and zips them, so one sweep point sets the channel the signal goes through *and* the channel the detector inverts -- which is exactly what a fading realization is. @@ -121,9 +121,9 @@ def average_ser(chain, n_rx, n_tx, snr_dB, stimulus, n_channels, seed=0): H = rayleigh_channel(n_rx, n_tx, rng=rng) draws.append((H, H)) chain.set_params(noise__sigma2=10 ** (-snr_dB / 10)) - results = sweep(chain, ("channel.H", "detector.H"), draws, - {"ser": compute_ser}, stimulus=stimulus, - reference="tx", seed=seed) + results = monte_carlo(chain, ("channel.H", "detector.H"), draws, + {"ser": compute_ser}, stimulus=stimulus, + reference="tx", seed=seed) return float(np.mean(results["ser"])) diff --git a/examples/mimo/one_shot_mimo.py b/examples/mimo/one_shot_mimo.py index f610dd1..7628568 100644 --- a/examples/mimo/one_shot_mimo.py +++ b/examples/mimo/one_shot_mimo.py @@ -3,7 +3,7 @@ import numpy as np import matplotlib.pyplot as plt -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolMapper @@ -25,7 +25,7 @@ H = rayleigh_channel(N_r, N_t, seed=0) -def link(detector): +def get_chain(detector): """The same link, closed by one detector or another. Only the last block changes, which is the point: a detector is a @@ -53,7 +53,7 @@ def link(detector): } chains = {} for name, detector in detectors.items(): - chains[name] = link(detector) + chains[name] = get_chain(detector) for name, chain in chains.items(): chain.seed(0) @@ -100,9 +100,9 @@ def average_ser(name, chain, snr_dB, seed=0): if name in NOISE_AWARE: params["detector.sigma2"] = noise_variance chain.set_params(**params) - results = sweep(chain, ("channel.H", "detector.H"), draws, - {"ser": compute_ser}, stimulus=(N_t, n_symbols), - reference="tx", seed=seed) + results = monte_carlo(chain, ("channel.H", "detector.H"), draws, + {"ser": compute_ser}, stimulus=(N_t, n_symbols), + reference="tx", seed=seed) return float(np.mean(results["ser"])) diff --git a/examples/ofdm/one_shot_ofdm.py b/examples/ofdm/one_shot_ofdm.py index 375f033..bb94b21 100644 --- a/examples/ofdm/one_shot_ofdm.py +++ b/examples/ofdm/one_shot_ofdm.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import numpy as np -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.channels import AWGN, FIRChannel, TappedDelayLineChannel from comnumpy.core.compensators import LinearEqualizer @@ -105,8 +105,8 @@ snr_list = np.arange(6, 22, 2) measured = {} for name, chain in (("single carrier", sc_chain), ("OFDM", ofdm_chain)): - runs = [sweep(chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser}, N, - reference="data_tx", seed=trial)["ser"] + runs = [monte_carlo(chain, "data_rx.snr_dB", snr_list, {"ser": compute_ser}, N, + reference="data_tx", seed=trial)["ser"] for trial in range(1, 3)] measured[name] = np.mean(runs, axis=0) diff --git a/examples/optical/CD_compensation_part1.py b/examples/optical/CD_compensation_part1.py index ef4d42d..ec072ff 100644 --- a/examples/optical/CD_compensation_part1.py +++ b/examples/optical/CD_compensation_part1.py @@ -78,10 +78,13 @@ epsilon_b = constellation.energy / k # update chain and compensator parameters - chain["generator"].M = constellation.order - chain["mapper"].alphabet = constellation - full_compensator1["demapper"].alphabet = constellation - full_compensator2["demapper"].alphabet = constellation + # through set_params, never by assignment (D50): SymbolMapper coerces + # its alphabet in __post_init__, and a direct write skips that -- the + # block then holds a Constellation where its forward indexes an array + chain.set_params(generator__M=constellation.order, + mapper__alphabet=constellation) + full_compensator1.set_params(demapper__alphabet=constellation) + full_compensator2.set_params(demapper__alphabet=constellation) # perform simulation ber_theo = [] @@ -96,7 +99,7 @@ # perform MC simulations N0 = epsilon_b/snr_per_bit # snr_bit = epsilon_b/N0 = (epsilon_s/log2(order))/N0 - chain["noise"].sigma2 = N0 # change noise variance + chain.set_params(noise__sigma2=N0) y = chain(N) # evaluate metric diff --git a/examples/optical/NLI_simulation.py b/examples/optical/NLI_simulation.py index c7fde91..09cb81b 100644 --- a/examples/optical/NLI_simulation.py +++ b/examples/optical/NLI_simulation.py @@ -11,28 +11,29 @@ from comnumpy.core.compensators import DataAidedPhaseCompensator from comnumpy.core.filters import SRRCFilter from comnumpy.core.generators import SymbolGenerator -from comnumpy.core.mappers import SymbolMapper -from comnumpy.core.metrics import compute_effective_snr, compute_ser +from comnumpy.core.mappers import SymbolDemapper, SymbolMapper +from comnumpy.core.metrics import ErrorCounter, compute_effective_snr from comnumpy.core.processors import Amplifier, Downsampler, Upsampler -from comnumpy.core.utils import Constellation, hard_projector +from comnumpy.core.utils import Constellation from comnumpy.core.visualizers import plot_error_rate from comnumpy.optical.dbp import DBP from comnumpy.optical.fiber import FiberSpec from comnumpy.optical.gn_model import (gn_model_nli_power, gn_model_snr, optimal_launch_power) from comnumpy.optical.links import FiberLink -from comnumpy.optical.utils import (compute_erbium_doped_fiber_N_ase, - dbm_to_watt, watt_to_dbm) +from comnumpy.optical.utils import (dbm_to_watt, launch_amplitude, + watt_to_dbm) img_dir = "../../docs/tutorials/img/" constellation = Constellation("QAM", 16) -N_s = 2**11 # 2048 per trial, 4 trials: SER floor 1.2e-4 +N_s = 2**11 * 6 # 12288 symbols per trial; over 4 trials + # the SER floor is one error in 49152, 2.0e-5 oversampling_sim = 6 oversampling_dsp = 2 NF_dB = 5 rolloff = 0.1 -StPS = 200 # forward propagation, the reference +StPS = 50 # forward: converged, see the tutorial R_s = 32e9 L_span = 100 N_span = 10 @@ -51,7 +52,7 @@ # receiver would re-do 192 split-step propagations instead of 32, and # the propagation is where all the time goes. So the chain is cut where # the physics ends and the DSP begins. -def get_unprocessed_chain(): +def get_channel(): """Symbols to the field the receiver sees, launch power included.""" return Sequential([ SymbolGenerator(constellation.order, name="data_tx"), @@ -81,7 +82,8 @@ def get_receiver(steps, linear_only, *, gain, reference): Downsampler(oversampling_dsp), Amplifier(gain), DataAidedPhaseCompensator(reference, name="phase"), - ]) + SymbolDemapper(constellation, name="data_rx"), + ], taps=["phase"]) # --- what the closed form expects of this link ------------------------- @@ -91,8 +93,7 @@ def get_receiver(steps, linear_only, *, gain, reference): # and it is *single*-polarization, so the scalar NLSE applies and the # weights are the ones of `polarizations=1` -- 5.3 dB more interference # at equal power. The ASE is over one polarization for the same reason. -ase_W = N_span * compute_erbium_doped_fiber_N_ase( - fiber.alpha_dB, L_span, NF_dB, nu=fiber.carrier_frequency_Hz) * R_s +ase_W = get_channel()["link"].budget(R_s)["ase_power_W"] eta = gn_model_nli_power(fiber, span_length_km=L_span, n_spans=N_span, powers_W=np.array([1e-3]), frequencies_Hz=np.array([fiber.carrier_frequency_Hz]), @@ -116,17 +117,22 @@ def get_receiver(steps, linear_only, *, gain, reference): "DBP, 50 steps/span": (50, False), } -channel = get_unprocessed_chain() +channel = get_channel() snr = {} -ser = {} +counters = {} elapsed = {} for name in receivers: snr[name] = np.zeros(len(dBm_list)) - ser[name] = np.zeros(len(dBm_list)) + # one counter per launch power: it accumulates errors and symbols + # over the four trials rather than averaging four rates, so the + # count that says whether a point means anything stays readable + counters[name] = [] + for _ in dBm_list: + counters[name].append(ErrorCounter()) elapsed[name] = 0.0 for index, dBm in enumerate(dBm_list): - amp = np.sqrt(dbm_to_watt(dBm)) + amp = launch_amplitude(dbm_to_watt(dBm)) for trial in range(N_trial): # the two fields differ only by the fibre's nonlinearity, and # the same seed gives them the same symbols and the same noise @@ -135,23 +141,29 @@ def get_receiver(steps, linear_only, *, gain, reference): channel.seed(index * N_trial + trial) channel.set_params(launch__gain=amp, link__use_only_linear=use_only_linear) - fields[use_only_linear] = channel(N_s * oversampling_sim) + fields[use_only_linear] = channel(N_s) symbols, reference = channel.tap("signal_tx"), channel.tap("data_tx") for name, (steps, linear_only) in receivers.items(): receiver = get_receiver(steps, linear_only, gain=1 / amp, reference=symbols) bound = name == "amplifier noise only" - estimate = receiver(fields[bound]) + detected = receiver(fields[bound]) elapsed[name] += receiver.elapsed_ - detected, _ = hard_projector(estimate, constellation) + estimate = receiver.tap("phase") snr[name][index] += compute_effective_snr(symbols, estimate) / N_trial - ser[name][index] += compute_ser(reference, detected) / N_trial + counters[name][index].update(reference, detected) snr_dB = {} for name, values in snr.items(): snr_dB[name] = 10 * np.log10(values) +ser = {} +for name, points in counters.items(): + ser[name] = np.zeros(len(points)) + for index, counter in enumerate(points): + ser[name][index] = counter.rate + header = "launch power [dBm] " for value in dBm_list: header += f"{value:7.1f}" @@ -174,6 +186,15 @@ def get_receiver(steps, linear_only, *, gain, reference): f"{watt_to_dbm(best_power):+.2f} dBm, dispersion compensation " f"{reference[best]:.2f} dB at {dBm_list[best]:+.1f} dBm") +# What the error-rate figure can resolve, in counts rather than in a +# caveat. At its own best power each receiver is where its curve bottoms +# out, and a point that saw no error has not measured a rate: it has run +# out of symbols. The counters kept the numbers, so the page can say so. +print("\nreceiver errors at its best power") +for name, points in counters.items(): + counter = points[int(np.argmax(snr_dB[name]))] + print(f"{name:24s} {counter.n_errors:7d} over {counter.n_symbols} symbols") + # The prediction, on the same axes as the measurement. It describes the # receiver that only undoes the dispersion -- the GN model counts the # nonlinear interference as noise, so it has nothing to say about a diff --git a/examples/optical/gn_model.py b/examples/optical/gn_model.py index 8bcf97c..3e73c72 100644 --- a/examples/optical/gn_model.py +++ b/examples/optical/gn_model.py @@ -1,7 +1,7 @@ import numpy as np import matplotlib.pyplot as plt -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.compensators import DataAidedComplexGainCompensator from comnumpy.core.filters import SRRCFilter @@ -15,8 +15,8 @@ from comnumpy.optical.gn_model import (gn_model_nli_power, gn_model_snr, optimal_launch_power) from comnumpy.optical.links import FiberLink -from comnumpy.optical.utils import (compute_erbium_doped_fiber_N_ase, - dbm_to_watt, watt_to_dbm) +from comnumpy.optical.utils import (dbm_to_watt, launch_amplitude, + watt_to_dbm) from comnumpy.optical.wdm import WDMGrid img_dir = "../../docs/tutorials/img/" @@ -59,13 +59,13 @@ def comb_eta(n_channels): def analytic_ase(n_spans=N_SPANS, bandwidth_Hz=BAUD): """Amplifier noise in the channel, over both polarizations. - ``compute_erbium_doped_fiber_N_ase`` returns the one-sided spectral - density of a single amplifier **per polarization**; the factor two - is the second polarization, which carries its own independent noise - and is the term the channel-power convention makes easy to forget. + The link budgets its own noise. ``polarizations=2`` is the term the + channel-power convention makes easy to forget: the spectral density + of an amplifier is **per polarization**, and a coherent link carries + two independent copies of it. """ - density = compute_erbium_doped_fiber_N_ase(SMF.alpha_dB, SPAN_KM, NF_dB) - return 2 * n_spans * density * bandwidth_Hz + link = FiberLink(n_spans, L_span=SPAN_KM, fs=FS, fiber=SMF, NF_dB=NF_dB) + return link.budget(bandwidth_Hz, polarizations=2)["ase_power_W"] grid = comb(9) @@ -115,12 +115,7 @@ def analytic_ase(n_spans=N_SPANS, bandwidth_Hz=BAUD): # Part 2 -- the simulation, which is here to disagree # =========================================================================== -def launch_gain(power_W): - """Amplifier gain for a channel power split over two polarizations.""" - return np.sqrt(np.asarray(power_W) / 2) - - -def link_chain(order=16, power_W=1e-3): +def get_chain(order=16, power_W=1e-3): """The whole link, from symbols to equalized symbols, as one chain.""" source = ([GaussianGenerator(1.0, name="tx")] if order is None else [SymbolGenerator(order, name="source"), @@ -129,7 +124,7 @@ def link_chain(order=16, power_W=1e-3): *source, Upsampler(OS, scale=np.sqrt(OS)), SRRCFilter(ROLLOFF, OS, N_h=40, method="fft"), - Amplifier(launch_gain(power_W), name="launch"), + Amplifier(launch_amplitude(power_W, polarizations=2), name="launch"), FiberLink(N_SPANS, L_span=SPAN_KM, StPS=40, fs=FS, fiber=SMF, NF_dB=NF_dB, name="fibre"), DBP(N_SPANS, L_span=SPAN_KM, StPS=1, fs=FS, fiber=SMF, @@ -147,12 +142,13 @@ def effective_snr_dB(sent, received): def measure(chain, powers_W, seed=3): - return sweep(chain, "launch.gain", launch_gain(powers_W), - {"snr_dB": effective_snr_dB}, STIMULUS, - reference="tx", seed=seed)["snr_dB"] + return monte_carlo(chain, "launch.gain", + launch_amplitude(powers_W, polarizations=2), + {"snr_dB": effective_snr_dB}, STIMULUS, + reference="tx", seed=seed)["snr_dB"] -chain = link_chain() +chain = get_chain() chain.summary(STIMULUS) chain.set_params(fibre__use_only_linear=True) @@ -198,13 +194,13 @@ def measure(chain, powers_W, seed=3): print("stimulus nonlinear SNR above the model") for order, label in ((4, "QPSK"), (16, "16QAM"), (64, "64QAM"), (256, "256QAM"), (None, "Gaussian")): - noiseless = link_chain(order) + noiseless = get_chain(order) noiseless.set_params(fibre__noise_scaling=0.0) value = measure(noiseless, [1e-3])[0] print(f"{label:12s} {value:10.2f} dB {value - nli_only_dB:+10.2f} dB") mermaid_dir = "../../docs/tutorials/mermaid/" with open(f"{mermaid_dir}/gn_model.mmd", "w") as stream: - stream.write(link_chain().to_mermaid()) + stream.write(get_chain().to_mermaid()) plt.show() diff --git a/examples/optical/one_shot_NLI.py b/examples/optical/one_shot_NLI.py index bfdd197..cf44567 100644 --- a/examples/optical/one_shot_NLI.py +++ b/examples/optical/one_shot_NLI.py @@ -10,26 +10,25 @@ from comnumpy.core.compensators import DataAidedPhaseCompensator from comnumpy.core.filters import SRRCFilter from comnumpy.core.generators import SymbolGenerator -from comnumpy.core.mappers import SymbolMapper +from comnumpy.core.mappers import SymbolDemapper, SymbolMapper from comnumpy.core.metrics import compute_effective_snr, compute_ser from comnumpy.core.processors import Amplifier, Downsampler, Upsampler -from comnumpy.core.utils import Constellation, hard_projector +from comnumpy.core.utils import Constellation from comnumpy.core.visualizers import plot_error_rate, plot_iq from comnumpy.optical.dbp import DBP from comnumpy.optical.fiber import FiberSpec from comnumpy.optical.links import FiberLink -from comnumpy.optical.utils import (compute_erbium_doped_fiber_N_ase, - dbm_to_watt) +from comnumpy.optical.utils import dbm_to_watt, launch_amplitude img_dir = "../../docs/tutorials/img/" constellation = Constellation("QAM", 16) -N_s = 2**9 # symbols per run +N_s = 2**9 * 6 # symbols per run oversampling_sim = 6 # samples per symbol in the channel oversampling_dsp = 2 # samples per symbol in the receiver NF_dB = 5 # amplifier noise figure in dB rolloff = 0.1 -StPS = 200 # split steps per span, forward +StPS = 50 # split steps per span, forward StPS_DBP = 100 # split steps per span, backward R_s = 10.7e9 # baud rate L_span = 80 # span length in km @@ -37,14 +36,13 @@ dBm = -3 fiber = FiberSpec() # the standard fibre, named so the budget sees it -N = N_s * oversampling_sim fs = R_s * oversampling_sim oversampling_ratio = oversampling_sim // oversampling_dsp -amp = np.sqrt(dbm_to_watt(dBm)) +amp = launch_amplitude(dbm_to_watt(dBm)) # --- the chain -------------------------------------------------------- -def get_full_chain(n_spans, *, steps=1, linear_only=True): +def get_chain(n_spans, *, steps=1, linear_only=True): """The whole link, from the symbols to the decisions. Transmitter, ``n_spans`` of fibre and amplifier, and the receiver @@ -76,17 +74,16 @@ def get_full_chain(n_spans, *, steps=1, linear_only=True): Downsampler(oversampling_dsp), Amplifier(1 / amp), DataAidedPhaseCompensator(name="phase"), - ], taps=["data_tx", "signal_tx", "rx_field", "phase"], + SymbolDemapper(constellation, name="data_rx"), + ], taps=["data_tx", "signal_tx", "rx_field", "phase", "data_rx"], wiring={"phase.reference": "signal_tx"}) def score(chain): """Effective SNR in dB and symbol error rate of a finished run.""" - estimate = chain.tap("phase") - detected, _ = hard_projector(estimate, constellation) - return (10 * np.log10(compute_effective_snr(chain.tap("signal_tx"), - estimate)), - compute_ser(chain.tap("data_tx"), detected)) + return (compute_effective_snr(chain.tap("signal_tx"), + chain.tap("phase"), unit="dB"), + compute_ser(chain.tap("data_tx"), chain.tap("data_rx"))) # --- 1. what the link does to the signal ------------------------------ @@ -109,8 +106,9 @@ def score(chain): snr_ase_only = {} print("spans measured ASE only the fibre SER phase time") for n_spans in spans: - chain = get_full_chain(n_spans).seed(0) - estimates[n_spans] = chain(N) + chain = get_chain(n_spans).seed(0) + chain(N_s) + estimates[n_spans] = chain.tap("phase") snr, ser = score(chain) snr_per_span[n_spans] = snr elapsed = chain.elapsed_ @@ -120,7 +118,7 @@ def score(chain): received_field = chain.tap("rx_field") chain.seed(0).set_params(link__use_only_linear=True) - chain(N) + chain(N_s) snr_ase_only[n_spans] = score(chain)[0] print(f"{n_spans:5d} {snr:8.2f} dB {snr_ase_only[n_spans]:7.2f} dB " @@ -136,21 +134,22 @@ def score(chain): # pulse shaping, resampling, matched filtering. It has to sit far above # every number in the table, or the chain is measuring its own filters # rather than the fibre. -floor = get_full_chain(1).seed(0) +floor = get_chain(1).seed(0) floor.set_params(link__use_only_linear=True, link__noise_scaling=0.0) -floor(N) +floor(N_s) print(f"\ndistortion floor of the chain, no noise and no fibre: " f"{score(floor)[0]:.1f} dB") -# The second is the amplifier noise against its closed form. N spans pile -# up N times the noise of one, in the symbol-rate bandwidth the matched -# filter collects, and that prediction owes nothing to the simulation: -# it comes from the noise figure and the span loss alone. -n_ase = compute_erbium_doped_fiber_N_ase(fiber.alpha_dB, L_span, NF_dB, - nu=fiber.carrier_frequency_Hz) +# The second is the amplifier noise against its closed form. The link +# budgets its own noise: `budget` asks it how much ASE it accumulates in +# the bandwidth a matched filter keeps -- the symbol rate, not the +# simulated `fs` -- and that prediction owes the simulation nothing, it +# comes from the noise figure and the span loss alone. print("\nspans ASE only P / P_ASE gap") for n_spans in spans: - predicted = 10 * np.log10(dbm_to_watt(dBm) / (n_spans * n_ase * R_s)) + link = get_chain(n_spans)["link"] + ase_W = link.budget(R_s)["ase_power_W"] + predicted = 10 * np.log10(dbm_to_watt(dBm) / ase_W) print(f"{n_spans:5d} {snr_ase_only[n_spans]:8.2f} dB {predicted:8.2f} dB " f"{snr_ase_only[n_spans] - predicted:+8.2f} dB") @@ -182,7 +181,7 @@ def score(chain): # Where that time goes. `profile_execution_time` runs the chain and times # each block on the way through, so the same pass answers the question. -profile = get_full_chain(N_span).seed(0).profile_execution_time(N) +profile = get_chain(N_span).seed(0).profile_execution_time(N_s) print("\nblock time") for block_id, elapsed in profile.items(): print(f"{block_id:22s} {1e3 * elapsed:8.1f} ms") @@ -191,9 +190,9 @@ def score(chain): # The same chain with the nonlinear term switched back on in the # receiver. Everything else -- the seed, the link, the phase correction # -- is unchanged, so the comparison is over the same realization. -back_propagated = get_full_chain(N_span, steps=StPS_DBP, +back_propagated = get_chain(N_span, steps=StPS_DBP, linear_only=False).seed(0) -profile_dbp = back_propagated.profile_execution_time(N) +profile_dbp = back_propagated.profile_execution_time(N_s) snr_dbp, ser_dbp = score(back_propagated) print(f"\ndispersion compensation SNR={snr_per_span[N_span]:5.2f} dB " f"receiver {1e3 * profile['dbp']:7.1f} ms") diff --git a/examples/simple/channel_coding.py b/examples/simple/channel_coding.py index e63c0bb..b21e58a 100644 --- a/examples/simple/channel_coding.py +++ b/examples/simple/channel_coding.py @@ -3,7 +3,7 @@ import numpy as np import matplotlib.pyplot as plt -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.channels import AWGN from comnumpy.core.generators import SymbolGenerator @@ -52,7 +52,7 @@ def row(label, values): print(row("beta_d ", spectrum.beta_d)) -def uncoded_chain(): +def get_uncoded_chain(): return Sequential([ SymbolGenerator(2, name="tx"), SymbolMapper(BPSK), @@ -61,7 +61,7 @@ def uncoded_chain(): ], taps=["tx"], name="uncoded BPSK") -def coded_chain(soft): +def get_coded_chain(soft): return Sequential([ SymbolGenerator(2, name="tx"), ConvolutionalEncoder((0o133, 0o171)), @@ -73,12 +73,12 @@ def coded_chain(soft): curves = {} -for label, chain, code_rate in (("uncoded", uncoded_chain(), 1.0), - ("hard-decision Viterbi", coded_chain(False), rate), - ("soft-decision Viterbi", coded_chain(True), rate)): +for label, chain, code_rate in (("uncoded", get_uncoded_chain(), 1.0), + ("hard-decision Viterbi", get_coded_chain(False), rate), + ("soft-decision Viterbi", get_coded_chain(True), rate)): start = time.perf_counter() - results = sweep(chain, "noise.snr_dB", snr_dB(ebn0_dB, code_rate), - {"ber": compute_ser}, n_bits, reference="tx", seed=4) + results = monte_carlo(chain, "noise.snr_dB", snr_dB(ebn0_dB, code_rate), + {"ber": compute_ser}, n_bits, reference="tx", seed=4) curves[label] = results["ber"] line = f"{label:24s} " for value in results["ber"]: @@ -115,9 +115,9 @@ def coded_chain(soft): SymbolDemapper(BPSK, soft=True), LDPCDecoder(H, n_iter=n_iter), ], taps=["tx"], name=f"LDPC, {n_iter} iterations") - results = sweep(link, "noise.snr_dB", snr_dB(ebn0_dB, ldpc_encoder.rate), - {"ber": compute_ser}, (n_frames, ldpc_encoder.k), - reference="tx", seed=6) + results = monte_carlo(link, "noise.snr_dB", snr_dB(ebn0_dB, ldpc_encoder.rate), + {"ber": compute_ser}, (n_frames, ldpc_encoder.k), + reference="tx", seed=6) ldpc_curves[f"LDPC (2040, {ldpc_encoder.k}), {n_iter} iterations"] = results["ber"] line = f"LDPC {n_iter:2d} iterations " for value in results["ber"]: @@ -134,7 +134,7 @@ def coded_chain(soft): plt.savefig(f"{img_dir}/channel_coding_fig2.png") mermaid_dir = "../../docs/tutorials/mermaid/" -for diagram_name, diagram_chain in [("channel_coding", coded_chain(True))]: +for diagram_name, diagram_chain in [("channel_coding", get_coded_chain(True))]: with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: stream.write(diagram_chain.to_mermaid()) diff --git a/examples/simple/monte_carlo_awgn.py b/examples/simple/monte_carlo_awgn.py index be0b9f1..5586120 100644 --- a/examples/simple/monte_carlo_awgn.py +++ b/examples/simple/monte_carlo_awgn.py @@ -1,31 +1,45 @@ import numpy as np import matplotlib.pyplot as plt -from comnumpy import sweep +from functools import partial +from scipy.special import erfc +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.generators import SymbolGenerator from comnumpy.core.mappers import SymbolMapper, SymbolDemapper -from comnumpy.core.utils import Constellation +from comnumpy.core.utils import Constellation, ebn0_to_snr_dB from comnumpy.core.channels import AWGN -from comnumpy.core.metrics import compute_ser +from comnumpy.core.metrics import compute_ser, compute_ber from comnumpy.core.visualizers import plot_error_rate +img_dir = "../../docs/tutorials/img/" # parameters N = 1000000 constellation = Constellation("QAM", 16) snr_dB_list = np.arange(0, 22) -# create chain -chain = Sequential([ - SymbolGenerator(constellation.order, name="tx"), - SymbolMapper(constellation), - AWGN(snr_dB=0, name="awgn_channel"), - SymbolDemapper(constellation), - ]) + +def get_chain(constellation): + """Symbols in, decided symbols out, over an AWGN channel. + + The last section of this page runs one of these per constellation + order, so the constellation is an argument rather than a global. The + transmitter is named because a sweep has to compare against what it + produced, and the channel is named because a sweep has to reconfigure + it. + """ + return Sequential([ + SymbolGenerator(constellation.order, name="tx"), + SymbolMapper(constellation), + AWGN(snr_dB=0, name="awgn_channel"), + SymbolDemapper(constellation), + ], taps=["tx"]) + + +chain = get_chain(constellation) # A Monte Carlo simulation is a stateless loop: reconfigure the chain, # run it, measure, collect. Written out, at three SNR values: -chain.taps = ["tx"] by_hand = [] for snr_dB in snr_dB_list[::8]: chain.seed(1) @@ -33,25 +47,42 @@ detected = chain(N) # run first: the tap is filled by it by_hand.append(compute_ser(chain.tap("tx"), detected)) -# sweep (decision D35) is that loop, and nothing more: the same +# monte_carlo (decision D35) is that loop, and nothing more: the same # reconfigure-reseed-run-collect, over every point, in one call. -results = sweep(chain, "awgn_channel.snr_dB", snr_dB_list, - {"ser": compute_ser}, N, reference="tx", seed=1) +results = monte_carlo(chain, "awgn_channel.snr_dB", snr_dB_list, + {"ser": compute_ser}, N, reference="tx", seed=1) ser_array = results["ser"] -loop_line = "loop :" +loop_line = "loop :" for value in by_hand: loop_line += f" {value:.3e}" print(loop_line) -sweep_line = "sweep :" +sweep_line = "monte_carlo :" for value in ser_array[::8]: sweep_line += f" {value:.3e}" print(sweep_line) -# The closed form comes from the constellation itself: it knows its -# family and its order, so the theory cannot describe a modulation other -# than the one the chain transmits. AWGN(snr_dB=) is a symbol SNR, hence -# per="symbol". -ser_theo_array = constellation.metrics(snr_dB_list, per="symbol")["ser"] + +def qfunc(x): + """Gaussian tail function, the one the textbooks write as Q.""" + return 0.5 * erfc(x / np.sqrt(2)) + + +# The closed form of Proakis, written out. A square M-QAM is two +# independent sqrt(M)-PAM constellations, so a symbol is correct only +# when both components are. AWGN(snr_dB=) is a symbol SNR and the +# expression is parameterized by the SNR per bit, hence the division by k. +M = constellation.order +k = int(np.log2(M)) +gamma_b = 10 ** (snr_dB_list / 10) / k +P_pam = 2 * (1 - 1 / np.sqrt(M)) * qfunc(np.sqrt(3 * k * gamma_b / (M - 1))) +ser_theo_array = 1 - (1 - P_pam) ** 2 + +# The same expression is already in the library, and the constellation +# knows which one applies to it. per="symbol" says that the swept SNR is +# the symbol SNR, and the conversion to Eb/N0 above happens inside. +from_object = constellation.metrics(snr_dB_list, per="symbol")["ser"] +gap = np.max(np.abs(ser_theo_array - from_object)) +print(f"largest gap between the two closed forms: {gap:.3e}") # plot_error_rate draws the measurements as markers and the closed form # as a line of the same colour: the figure every sweep ends with @@ -59,6 +90,49 @@ plot_error_rate(snr_dB_list, {label: ser_array}, theory={label: ser_theo_array}, ylabel="SER", title=f"SER performance for {label}") +plt.savefig(f"{img_dir}/monte_carlo_awgn.png") + +# Comparing orders is only fair at equal energy per *bit*: at equal +# energy per symbol a dense constellation is favoured, because it spends +# the same energy on more bits. ebn0_to_snr_dB does that conversion, and +# it needs k, which is chain-level knowledge the channel does not have. +ebn0_dB_list = np.arange(0, 25, 2) +N_compare = 100000 +fine_dB_list = np.arange(0, 25, 0.05) +measured_ber = {} +theory_ber = {} +required_rows = [] +for order in [4, 16, 64, 256]: + other = Constellation("QAM", order) + bits = int(np.log2(order)) + snr_dB_values = ebn0_to_snr_dB(ebn0_dB_list, bits_per_symbol=bits) + collected = monte_carlo(get_chain(other), "awgn_channel.snr_dB", + snr_dB_values, + {"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"] + # What the density costs, read off the closed form rather than the + # measurement: the sweep grid is 2 dB coarse and noisy at 1e-3. + fine_ber = other.metrics(fine_dB_list, per="bit")["ber"] + above_underflow = fine_ber > 0 # the tail reaches 0 in float64 + required_dB = float( + np.interp(-3.0, np.log10(fine_ber[above_underflow])[::-1], + fine_dB_list[above_underflow][::-1])) + required_rows.append((order_label, bits, required_dB)) + +print(f"{'order':>8s} {'bits/symbol':>12s} {'Eb/N0 at BER=1e-3':>19s}") +for order_label, bits, ebn0_dB in required_rows: + print(f"{order_label:>8s} {bits:>12d} {ebn0_dB:>16.1f} dB") + +ax = plot_error_rate(ebn0_dB_list, measured_ber, theory=theory_ber, + 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 +# past a few 1e-6, so the axis stops where the simulation stops speaking. +ax.set_ylim(1e-6, 1) +plt.savefig(f"{img_dir}/monte_carlo_awgn_orders.png") # The chain diagrams this tutorial shows are exported from the chains # themselves (D33c), so the picture cannot drift from the code -- the @@ -67,5 +141,3 @@ for diagram_name, diagram_chain in [("awgn_chain", chain)]: with open(f"{mermaid_dir}/{diagram_name}.mmd", "w") as stream: stream.write(diagram_chain.to_mermaid()) - -plt.show() diff --git a/pyright_all_tmp.json b/pyright_all_tmp.json deleted file mode 100644 index 32b4a1b..0000000 --- a/pyright_all_tmp.json +++ /dev/null @@ -1,62 +0,0 @@ -{ - "include": [ - "src/comnumpy/_backend.py", - "src/comnumpy/core/capacity.py", - "src/comnumpy/core/channels.py", - "src/comnumpy/core/compensators.py", - "src/comnumpy/core/devices.py", - "src/comnumpy/core/fading.py", - "src/comnumpy/core/filters.py", - "src/comnumpy/core/frames.py", - "src/comnumpy/core/generators.py", - "src/comnumpy/core/generics.py", - "src/comnumpy/core/impairments.py", - "src/comnumpy/core/mappers.py", - "src/comnumpy/core/metrics.py", - "src/comnumpy/core/processors.py", - "src/comnumpy/core/sequences.py", - "src/comnumpy/core/utils.py", - "src/comnumpy/core/validators.py", - "src/comnumpy/core/visualizers.py", - "src/comnumpy/exceptions.py", - "src/comnumpy/fec/analysis.py", - "src/comnumpy/fec/convolutional.py", - "src/comnumpy/fec/ldpc.py", - "src/comnumpy/mimo/channels.py", - "src/comnumpy/mimo/compensators.py", - "src/comnumpy/mimo/detectors.py", - "src/comnumpy/mimo/utils.py", - "src/comnumpy/mimo/validators.py", - "src/comnumpy/ofdm/allocation.py", - "src/comnumpy/ofdm/chains.py", - "src/comnumpy/ofdm/compensators.py", - "src/comnumpy/ofdm/metrics.py", - "src/comnumpy/ofdm/predistorders.py", - "src/comnumpy/ofdm/processors.py", - "src/comnumpy/ofdm/utils.py", - "src/comnumpy/ofdm/visualizers.py", - "src/comnumpy/optical/channels.py", - "src/comnumpy/optical/compensators.py", - "src/comnumpy/optical/constants.py", - "src/comnumpy/optical/dbp.py", - "src/comnumpy/optical/devices.py", - "src/comnumpy/optical/fiber.py", - "src/comnumpy/optical/links.py", - "src/comnumpy/optical/raman.py", - "src/comnumpy/optical/utils.py", - "src/comnumpy/optical/wdm.py", - "src/comnumpy/serialization.py", - "src/comnumpy/style.py", - "src/comnumpy/sweep.py" - ], - "typeCheckingMode": "strict", - "pythonVersion": "3.11", - "reportUnknownArgumentType": false, - "reportUnknownMemberType": false, - "reportUnknownVariableType": false, - "reportUnknownParameterType": false, - "reportUnknownLambdaType": false, - "reportMissingTypeStubs": false, - "reportConstantRedefinition": false, - "reportUnusedFunction": false -} \ No newline at end of file diff --git a/pyrightconfig.json b/pyrightconfig.json index 3e330be..ca784b0 100644 --- a/pyrightconfig.json +++ b/pyrightconfig.json @@ -20,7 +20,7 @@ "src/comnumpy/fec/ldpc.py", "src/comnumpy/fec/analysis.py", "src/comnumpy/core/devices.py", - "src/comnumpy/sweep.py", + "src/comnumpy/monte_carlo.py", "src/comnumpy/core/filters.py", "src/comnumpy/ofdm/compensators.py", "src/comnumpy/ofdm/predistorders.py", diff --git a/src/comnumpy/__init__.py b/src/comnumpy/__init__.py index 2c4a51f..e6f8f64 100644 --- a/src/comnumpy/__init__.py +++ b/src/comnumpy/__init__.py @@ -9,7 +9,7 @@ from .exceptions import ComnumpyError, ShapeError, NotFittedError from .serialization import to_json, from_json -from .sweep import sweep +from .monte_carlo import monte_carlo from .core import ( Processor, Sequential, SymbolGenerator, SymbolMapper, SymbolDemapper, @@ -23,5 +23,5 @@ "SymbolGenerator", "SymbolMapper", "SymbolDemapper", "AWGN", "compute_ser", "Constellation", "get_alphabet", "ebn0_to_snr_dB", "esn0_to_snr_dB", - "to_json", "from_json", "sweep", + "to_json", "from_json", "monte_carlo", ] diff --git a/src/comnumpy/core/generics.py b/src/comnumpy/core/generics.py index 64b7eec..bcd08e1 100644 --- a/src/comnumpy/core/generics.py +++ b/src/comnumpy/core/generics.py @@ -332,7 +332,7 @@ def _split_parameter(self, key: str) -> tuple[str, str]: """Split ``"block.field"`` or ``"block__field"`` (decision D34). The dotted form is the addressing used everywhere a parameter is - held in a variable -- :func:`~comnumpy.sweep.sweep` builds those + held in a variable -- :func:`~comnumpy.monte_carlo.sweep` builds those strings. A dot is not a valid character in a keyword argument, though, so writing one by hand costs a ``**{...}`` wrapper. The double underscore is the same address spelled as an identifier, diff --git a/src/comnumpy/core/metrics.py b/src/comnumpy/core/metrics.py index 24d55ff..7447890 100644 --- a/src/comnumpy/core/metrics.py +++ b/src/comnumpy/core/metrics.py @@ -1,8 +1,10 @@ +from dataclasses import dataclass, field from typing import Optional, Sequence import numpy as np from comnumpy.core.utils import sym_2_bin # single definition (annex A.5) +from comnumpy.exceptions import NotFittedError __all__ = [ "compute_ser_awgn_psk", "compute_ser_awgn_qam", "compute_metric_awgn_theo", @@ -10,6 +12,7 @@ "compute_metric_rayleigh_theo", "compute_ser", "compute_ber", "compute_evm", "compute_effective_snr", "compute_power", "compute_ccdf", "compute_acpr", "signal_report", + "ErrorCounter", ] # Craig's finite-range form of the Q function turns every average over @@ -1204,3 +1207,145 @@ def signal_report(x: np.ndarray, compute_papr: bool = False, from comnumpy.ofdm.metrics import compute_papr as papr_of report["papr"] = float(papr_of(abs_x, unit=papr_unit)) return report + + +@dataclass(slots=True) +class ErrorCounter: + r"""Accumulate errors over the trials of one Monte-Carlo point. + + Signal Model + ------------ + A Monte-Carlo point runs the chain several times and reports one + error rate. The rate is a ratio of **totals**, not a mean of ratios: + + .. math:: + + \mathrm{rate} = \frac{\sum_t N_{e,t}}{\sum_t N_t} + + Averaging the per-trial rates instead gives the same answer only + when every trial has the same length, and it throws away the number + the point's credibility rests on -- how many errors were actually + seen. A rate of :math:`10^{-4}` from three errors and one from three + hundred are different claims, and only the counts tell them apart. + + That is what this object is for. It does not own the loop: a + Monte-Carlo over a fading channel, over several detectors or over a + split-step link has a different shape each time, and that shape is + what the reader of an example came to see. It owns the counting, + which is the same everywhere and the part that hides mistakes. + + Parameters + ---------- + width : int, optional + Bits per symbol. The default 1 counts **symbol** errors, one per + mismatched entry. Give the constellation's ``bits_per_symbol`` + to count **bit** errors instead, expanded by + :func:`~comnumpy.core.utils.sym_2_bin` exactly as + :func:`compute_ber` does. + + Attributes + ---------- + n_errors : int + Errors accumulated so far -- symbols or bits, per ``width``. + n_symbols : int + Symbols compared so far. + n_trials : int + Number of :meth:`update` calls. + + Raises + ------ + ValueError + If ``width`` is not positive. + + Examples + -------- + Three trials of a hundred symbols, one error in the second: + + >>> import numpy as np + >>> counter = ErrorCounter() + >>> sent = np.zeros(100, dtype=int) + >>> for trial in range(3): + ... received = sent.copy() + ... if trial == 1: + ... received[7] = 1 + ... _ = counter.update(sent, received) + >>> print(counter) + ErrorCounter(1 error / 300 symbols, rate 3.33e-03) + + The counts are what makes the rate readable, and they are what an + example should print beside it: + + >>> print(counter.n_errors, counter.n_symbols, counter.n_trials) + 1 300 3 + + ``width`` switches to bits, on the same convention as + :func:`compute_ber`: + + >>> bits = ErrorCounter(width=4) + >>> _ = bits.update(np.array([0, 0]), np.array([0, 15])) + >>> print(bits.n_errors, bits.n_symbols, f"{bits.rate:.3f}") + 4 2 0.500 + """ + width: int = 1 + n_errors: int = field(init=False, default=0) + n_symbols: int = field(init=False, default=0) + n_trials: int = field(init=False, default=0) + + def __post_init__(self) -> None: + if self.width < 1: + raise ValueError( + f"a symbol carries at least one bit, got width={self.width}. " + f"Leave it at 1 to count symbol errors, or pass the " + f"constellation's bits_per_symbol to count bit errors.") + + def update(self, reference: np.ndarray, + detected: np.ndarray) -> "ErrorCounter": + """Count one trial, and return self so calls can be chained. + + Both arrays carry symbol **indices**, as :func:`compute_ser` + requires: the output of a ``SymbolGenerator`` against that of a + ``SymbolDemapper``, not complex constellation points. The + 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() + length = min(target.size, estimate.size) + if length == 0: + raise ValueError( + "nothing to count: one of the two sequences is empty.") + target, estimate = target[:length], estimate[:length] + if self.width == 1: + errors = int(np.count_nonzero(target - estimate)) + else: + errors = int(np.count_nonzero( + sym_2_bin(target, self.width) - sym_2_bin(estimate, self.width))) + self.n_errors += errors + self.n_symbols += length + self.n_trials += 1 + return self + + @property + def rate(self) -> float: + """Errors over the units compared -- symbols, or bits if ``width``.""" + if self.n_symbols == 0: + raise NotFittedError( + "no trial has been counted yet, so there is no rate; call " + "update(reference, detected) at least once.") + return self.n_errors / (self.n_symbols * self.width) + + def reset(self) -> "ErrorCounter": + """Forget every trial, keeping ``width``. Returns self.""" + self.n_errors = 0 + self.n_symbols = 0 + self.n_trials = 0 + return self + + def __str__(self) -> str: + unit = "symbol" if self.width == 1 else "bit" + plural = "" if self.n_errors == 1 else "s" + if self.n_symbols == 0: + return "ErrorCounter(nothing counted yet)" + return (f"ErrorCounter({self.n_errors} error{plural} / " + f"{self.n_symbols * self.width} {unit}s, " + f"rate {self.rate:.2e})") diff --git a/src/comnumpy/sweep.py b/src/comnumpy/monte_carlo.py similarity index 88% rename from src/comnumpy/sweep.py rename to src/comnumpy/monte_carlo.py index 150aed0..46e9ef5 100644 --- a/src/comnumpy/sweep.py +++ b/src/comnumpy/monte_carlo.py @@ -14,6 +14,7 @@ """ from __future__ import annotations +import logging import multiprocessing import os import pickle @@ -25,7 +26,10 @@ from comnumpy.core.generics import Sequential -__all__ = ["sweep"] +__all__ = ["monte_carlo"] + +logger = logging.getLogger(__name__) + # Every worker runs one sweep point at a time, so a BLAS that spawns one # thread per core inside each of them oversubscribes the machine by @@ -75,15 +79,15 @@ def _evaluate(chain: Sequential, params: list[str], point: Any, for name, metric in metrics.items()} -def sweep(chain: Sequential, - param: str | Sequence[str], - values: Sequence[Any], - metrics: Mapping[str, Callable[..., Any]], - stimulus: Any, - *, - reference: Optional[str] = None, - seed: Optional[int] = None, - n_jobs: int = 1) -> dict[str, np.ndarray]: +def monte_carlo(chain: Sequential, + param: str | Sequence[str], + values: Sequence[Any], + metrics: Mapping[str, Callable[..., Any]], + stimulus: Any, + *, + reference: Optional[str] = None, + seed: Optional[int] = None, + n_jobs: int = 1) -> dict[str, np.ndarray]: """Run a chain over a range of parameter values and collect metrics. Parameters @@ -174,14 +178,25 @@ def sweep(chain: Sequential, >>> from comnumpy import AWGN, Sequential, SymbolGenerator >>> from comnumpy.core.metrics import compute_ser >>> chain = Sequential([SymbolGenerator(4), AWGN(snr_dB=0, name="noise")]) - >>> out = sweep(chain, "noise.snr_dB", [0, 10], {"power": lambda y: float(np.mean(np.abs(y)**2))}, - ... stimulus=1000, seed=1) + >>> out = monte_carlo(chain, "noise.snr_dB", [0, 10], + ... {"power": lambda y: float(np.mean(np.abs(y)**2))}, + ... stimulus=1000, seed=1) >>> print(out["power"].shape) (2,) """ params = [param] if isinstance(param, str) else list(param) if reference is not None and reference not in (chain.taps or []): chain.taps = (chain.taps or []) + [reference] + if seed is None: + # Not an error: a first look at a curve does not need to be + # reproducible. But a curve nobody can reproduce is not a result, + # and nothing else would ever say so -- the parallel path refuses + # outright, the serial one used to stay silent. + logger.warning( + "monte_carlo() ran without seed=, so every point drew freely and " + "this curve cannot be reproduced -- not by you, not by a reader " + "of the page it lands on. Pass seed= to give each point its own " + "child seed (decision D6/D35).") seeds = (np.random.SeedSequence(seed).spawn(len(values)) if seed is not None else [None] * len(values)) @@ -208,7 +223,7 @@ def _sweep_parallel(chain: Sequential, params: list[str], """The same points, on several processes. Same values, out of order.""" if any(child is None for child in children): raise ValueError( - f"sweep(n_jobs={n_jobs}) needs seed= as well: the points are " + f"monte_carlo(n_jobs={n_jobs}) needs seed= as well: the points are " f"then drawn from independent child seeds and the curve is the " f"same as the serial one. Without it, what a point draws would " f"depend on which worker ran it.") @@ -217,7 +232,7 @@ def _sweep_parallel(chain: Sequential, params: list[str], pickle.dumps(payload) except Exception as error: # noqa: BLE001 -- re-raised raise ValueError( - f"sweep(n_jobs={n_jobs}) runs the points in worker " + f"monte_carlo(n_jobs={n_jobs}) runs the points in worker " f"processes, and the {label} cannot be sent to one: " f"{error}. A lambda or a closure cannot be pickled -- " f"define the metric at module level, or keep n_jobs=1." @@ -250,7 +265,7 @@ def _sweep_parallel(chain: Sequential, params: list[str], # under a __main__ guard runs itself again inside every # worker. raise RuntimeError( - "sweep(n_jobs>1) starts worker processes, and each one " + "monte_carlo(n_jobs>1) starts worker processes, and each one " "re-imports the module it was started from. If you are " "calling this from a script, its body must sit under\n\n" " if __name__ == \"__main__\":\n\n" diff --git a/src/comnumpy/optical/dbp.py b/src/comnumpy/optical/dbp.py index 6b54a4f..e2d1511 100644 --- a/src/comnumpy/optical/dbp.py +++ b/src/comnumpy/optical/dbp.py @@ -4,7 +4,8 @@ from comnumpy.core import Processor from .fiber import FiberSpec from .utils import (get_linear_step_size, get_logarithmic_step_size, compute_erbium_doped_fiber_amplifier_gain, - apply_chromatic_dispersion, apply_kerr_nonlinearity, + step_transfers, apply_frequency_response, TransferKey, + apply_kerr_nonlinearity, is_polarization_pair, manakov_kerr) __all__ = ["DBP"] @@ -45,7 +46,13 @@ class DBP(Processor): cannot end up integrating different equations. Axes: *declared axis* -- a full-field signal ``(N,)`` or ``(..., P, N)`` - with ``P`` in ``{1, 2}`` polarizations. + with ``P`` in ``{1, 2}`` polarizations. The leading axes are a + **batch**: ``(B, P, N)`` propagates ``B`` realizations in one call, + which is 1.8x faster than the same ``B`` calls at 12288 samples and + 3.1x at 1024 -- the split-step loop, the per-span amplifier and the + parameter precomputation are paid once instead of ``B`` times. A + batch of single-polarization fields is ``(B, 1, N)``; ``(B, N)`` + reads ``B`` as polarizations and is refused. Parameters ---------- @@ -133,6 +140,9 @@ class DBP(Processor): beta2: Optional[float] = field(init=False, repr=False, default_factory=lambda: None) gain: Optional[float] = field(init=False, repr=False, default_factory=lambda: None) manakov_: bool = field(init=False, repr=False, default=False) + transfer_: Optional[np.ndarray] = field(init=False, repr=False, default=None) + transfer_index_: Optional[np.ndarray] = field(init=False, repr=False, default=None) + transfer_key_: Optional[TransferKey] = field(init=False, repr=False, default=None) def prepare(self, x: np.ndarray) -> None: self.manakov_ = is_polarization_pair(x, type(self).__name__) @@ -148,6 +158,7 @@ def prepare(self, x: np.ndarray) -> None: self.beta2 = self.fiber.beta2 self.gain = 1/edfa_gain self.step_size = step_size[::-1] # reverse order + self._build_transfer(x.shape[-1]) def _apply_kerr(self, y: np.ndarray, dz: float) -> np.ndarray: """Backward Kerr step: scalar NLSE, or Manakov for a pair.""" @@ -155,26 +166,47 @@ def _apply_kerr(self, y: np.ndarray, dz: float) -> np.ndarray: return apply_kerr_nonlinearity(y, dz, gamma, direction=-1, intensity=intensity) + def _build_transfer(self, n_samples: int) -> None: + """Cache one linear-step transfer function per distinct length.""" + assert self.step_size is not None and self.beta2 is not None + if self.use_only_linear: + lengths = np.array([self.L_span], dtype=float) + elif self.step_method == "symmetric": + lengths = np.asarray(self.step_size, dtype=float) / 2 + else: + lengths = np.asarray(self.step_size, dtype=float) + previous = None + if (self.transfer_ is not None and self.transfer_index_ is not None + and self.transfer_key_ is not None): + previous = (self.transfer_, self.transfer_index_, + self.transfer_key_) + self.transfer_, self.transfer_index_, self.transfer_key_ = step_transfers( + n_samples, lengths, beta2=self.beta2, fs=self.fs, + alpha_dB=self.fiber.alpha_dB, direction=-1, previous=previous) + def forward(self, x: np.ndarray) -> np.ndarray: y = x # set by prepare(), which the Processor base always runs first assert (self.gain is not None and self.beta2 is not None - and self.step_size is not None) + and self.step_size is not None + and self.transfer_ is not None + and self.transfer_index_ is not None) for _ in range(self.N_spans): y = self.gain * y # correct for edfa gain if self.use_only_linear: - y = apply_chromatic_dispersion(y, self.L_span, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=-1) + y = apply_frequency_response(y, self.transfer_[0]) else: for num_step in range(self.StPS): dz = self.step_size[num_step] + H = self.transfer_[self.transfer_index_[num_step]] if self.step_method == "symmetric": - y = apply_chromatic_dispersion(y, dz/2, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=-1) + y = apply_frequency_response(y, H) y = self._apply_kerr(y, dz) - y = apply_chromatic_dispersion(y, dz/2, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=-1) + y = apply_frequency_response(y, H) if self.step_method == "asymetric": - y = apply_chromatic_dispersion(y, dz, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=-1) + y = apply_frequency_response(y, H) y = self._apply_kerr(y, dz) return y diff --git a/src/comnumpy/optical/links.py b/src/comnumpy/optical/links.py index d4cecad..afe623e 100644 --- a/src/comnumpy/optical/links.py +++ b/src/comnumpy/optical/links.py @@ -2,16 +2,17 @@ import numpy as np from dataclasses import dataclass, field -from typing import Optional, Literal, Callable, Dict -from comnumpy._backend import fft, fftfreq, ifft # cupy-compatible (D3) +from typing import Any, Callable, Dict, Literal, Optional +from comnumpy._backend import fftfreq # cupy-compatible (D3) from comnumpy.core import Processor from .devices import ErbiumDopedFiberAmplifier from .fiber import FiberSpec from .raman import RamanSolution from .constants import PLANCK_CONSTANT from .utils import (get_linear_step_size, get_logarithmic_step_size, compute_erbium_doped_fiber_amplifier_gain, - compute_erbium_doped_fiber_N_ase, apply_chromatic_dispersion, apply_kerr_nonlinearity, - is_polarization_pair, manakov_kerr) + compute_erbium_doped_fiber_N_ase, apply_kerr_nonlinearity, + is_polarization_pair, manakov_kerr, watt_to_dbm, + step_transfers, apply_frequency_response, TransferKey) __all__ = ["FiberLink"] @@ -74,7 +75,13 @@ class FiberLink(Processor): polarization and keeps the scalar NLSE with :math:`\gamma`. Axes: *declared axis* -- a full-field signal ``(N,)`` or ``(..., P, N)`` - with ``P`` in ``{1, 2}`` polarizations. + with ``P`` in ``{1, 2}`` polarizations. The leading axes are a + **batch**: ``(B, P, N)`` propagates ``B`` realizations in one call, + which is 1.8x faster than the same ``B`` calls at 12288 samples and + 3.1x at 1024 -- the split-step loop, the per-span amplifier and the + parameter precomputation are paid once instead of ``B`` times. A + batch of single-polarization fields is ``(B, 1, N)``; ``(B, N)`` + reads ``B`` as polarizations and is refused. Parameters ---------- @@ -206,6 +213,9 @@ class FiberLink(Processor): raman_tilt_: Optional[np.ndarray] = field(init=False, repr=False, default=None) manakov_: bool = field(init=False, repr=False, default=False) raman_sigma2_: float = field(init=False, repr=False, default=0.0) + transfer_: Optional[np.ndarray] = field(init=False, repr=False, default=None) + transfer_index_: Optional[np.ndarray] = field(init=False, repr=False, default=None) + transfer_key_: Optional[TransferKey] = field(init=False, repr=False, default=None) rng_: Optional[np.random.Generator] = field(init=False, repr=False, default=None) def prepare(self, x: np.ndarray) -> None: @@ -253,10 +263,7 @@ def prepare(self, x: np.ndarray) -> None: # bandwidth, spread over the simulated one. Averaged over the # channels and added flat across the band: the gain is shaped # in frequency, this noise is not. - if self.raman.bandwidth_Hz > 0: - ase_W = float(np.mean(np.atleast_2d(self.raman.ase_W)[:, -1])) - self.raman_sigma2_ = (self.noise_scaling * ase_W - * self.fs / self.raman.bandwidth_Hz) + self.raman_sigma2_ = self.fs * self._ase_densities()[2] if residual_dB < 0: logger.warning( "Raman over-compensates the span: %.2f dB of on-off gain " @@ -267,11 +274,154 @@ def prepare(self, x: np.ndarray) -> None: # the EDFA makes up whatever the Raman gain did not, so a span # stays transparent whether or not it is Raman-pumped - equivalent_alpha_dB = residual_dB / self.L_span + equivalent_alpha_dB, edfa_density, _ = self._ase_densities() self.edfa_gain = compute_erbium_doped_fiber_amplifier_gain(equivalent_alpha_dB, self.L_span) - self.edfa_N_ase = self.noise_scaling * self.fs * compute_erbium_doped_fiber_N_ase( + # a density becomes a power in the simulated bandwidth; `budget` + # multiplies the same density by the bandwidth a receiver keeps + self.edfa_N_ase = self.fs * edfa_density + self._build_transfer(x.shape[-1]) + + def _build_transfer(self, n_samples: int) -> None: + """Cache one linear-step transfer function per distinct length.""" + assert self.step_size is not None and self.beta2 is not None + if self.use_only_linear: + lengths = np.array([self.L_span], dtype=float) + elif self.step_method == "symmetric": + lengths = np.asarray(self.step_size, dtype=float) / 2 + else: + lengths = np.asarray(self.step_size, dtype=float) + previous = None + if (self.transfer_ is not None and self.transfer_index_ is not None + and self.transfer_key_ is not None): + previous = (self.transfer_, self.transfer_index_, + self.transfer_key_) + self.transfer_, self.transfer_index_, self.transfer_key_ = step_transfers( + n_samples, lengths, beta2=self.beta2, fs=self.fs, + alpha_dB=self.fiber.alpha_dB, direction=1, previous=previous) + + def _ase_densities(self) -> tuple[float, float, float]: + """One span's noise, before any bandwidth is chosen. + + Returns the equivalent attenuation the EDFA has to make up, and + the one-sided ASE spectral densities of the EDFA and of the + distributed Raman amplifier, in W/Hz **per polarization**, both + already scaled by ``noise_scaling``. + + It exists so that :meth:`prepare` and :meth:`budget` cannot + disagree: the first turns these densities into a variance over + ``fs``, the second into a power over the bandwidth a receiver + collects, and the physics is written once. + """ + residual_dB = self.fiber.alpha_dB * self.L_span + raman_density = 0.0 + if self.raman is not None: + residual_dB -= float(np.mean(np.atleast_1d(self.raman.on_off_gain_dB))) + if self.raman.bandwidth_Hz > 0: + ase_W = float(np.mean(np.atleast_2d(self.raman.ase_W)[:, -1])) + raman_density = (self.noise_scaling * ase_W + / self.raman.bandwidth_Hz) + equivalent_alpha_dB = residual_dB / self.L_span + edfa_density = self.noise_scaling * compute_erbium_doped_fiber_N_ase( equivalent_alpha_dB, self.L_span, self.NF_dB, h=PLANCK_CONSTANT, nu=self.fiber.carrier_frequency_Hz) + return equivalent_alpha_dB, edfa_density, raman_density + + def budget(self, bandwidth_Hz: float, *, + polarizations: int = 1) -> Dict[str, Any]: + r"""The amplifier noise this link accumulates, in closed form. + + Signal Model + ------------ + Each of the :math:`N_s` spans adds one amplifier's worth of + spontaneous emission, and the noise a receiver sees is that + density integrated over the bandwidth it keeps -- the symbol + rate, for a matched filter: + + .. math:: + + P_{\mathrm{ASE}} = P \, N_s \, N_{\mathrm{ASE}} \, B, + \qquad + N_{\mathrm{ASE}} = (G - 1) h \nu \, n_{sp}, + \qquad + n_{sp} = \frac{\mathrm{NF}/2}{1 - 1/G} + + with :math:`P` the number of polarizations. That factor is the + whole reason this method exists. :math:`N_{\mathrm{ASE}}` is a + density **per polarization**, and a dual-polarization link + carries two independent copies of it; quoting a channel power + against a single-polarization budget is a 3 dB error that looks + exactly like a modelling disagreement. Writing the product by + hand at each call site is how the two conventions drift apart. + + The result is the link's *own* prediction: it owes nothing to a + simulation, so it is the reference a measured effective SNR is + checked against. ``noise_scaling`` is applied, so a link with + its noise switched off budgets zero. + + Parameters + ---------- + bandwidth_Hz : float + Bandwidth the receiver collects, in Hz. For a matched filter + this is the symbol rate, not the simulated ``fs``. + polarizations : int, optional, keyword-only + 1 (default) for a single-polarization link, 2 for the + dual-polarization link a coherent system uses -- the same + convention as + :func:`~comnumpy.optical.gn_model.gn_model_nli_power`. + + Returns + ------- + dict + ``ase_power_W`` and ``ase_power_dBm``, the accumulated noise; + ``ase_density_W_per_Hz``, one span's density per + polarization; and the inputs it was built from + (``n_spans``, ``NF_dB``, ``span_gain_dB``, ``bandwidth_Hz``, + ``polarizations``). + + Raises + ------ + ValueError + If the bandwidth is not positive, or if ``polarizations`` is + neither 1 nor 2: a fibre carries one or two. + + Examples + -------- + Ten 100 km spans with a 5 dB noise figure, one polarization, + read in a 32 GBd matched filter: + + >>> link = FiberLink(10, L_span=100.0, NF_dB=5.0, fs=192e9) + >>> budget = link.budget(32e9) + >>> print(f"{budget['ase_power_dBm']:.2f} dBm") + -21.88 dBm + + The second polarization carries its own noise, and nothing else + changes: + + >>> print(f"{link.budget(32e9, polarizations=2)['ase_power_dBm']:.2f} dBm") + -18.87 dBm + """ + if bandwidth_Hz <= 0: + raise ValueError( + f"a noise power needs a bandwidth to be quoted in, got " + f"bandwidth_Hz={bandwidth_Hz}. For a matched filter this " + f"is the symbol rate, not the simulated fs.") + if polarizations not in (1, 2): + raise ValueError( + f"a fibre carries one or two polarizations, got " + f"{polarizations}.") + equivalent_alpha_dB, edfa_density, raman_density = self._ase_densities() + density = edfa_density + raman_density + power = polarizations * self.N_spans * density * bandwidth_Hz + return { + "ase_power_W": power, + "ase_power_dBm": float(watt_to_dbm(power)) if power > 0 else -np.inf, + "ase_density_W_per_Hz": density, + "span_gain_dB": equivalent_alpha_dB * self.L_span, + "n_spans": self.N_spans, + "NF_dB": self.NF_dB, + "bandwidth_Hz": bandwidth_Hz, + "polarizations": polarizations, + } def _raman_step_gain(self) -> np.ndarray: """Amplitude gain of each SSFM step, from the Raman profile.""" @@ -285,7 +435,7 @@ def _raman_step_gain(self) -> np.ndarray: f"describe this fibre") # Sampled at the *half*-step boundaries, not the step boundaries. # The gain belongs to the linear operator, exactly like the loss - # that apply_chromatic_dispersion already applies over dz/2 in + # that the linear step already applies over dz/2 in # each half; applying a whole step's gain at one point breaks the # symmetry of the symmetric split-step and drops it from second # order to first. Measured on the SPM phase of a CW field: the @@ -347,18 +497,20 @@ def _apply_raman(self, y: np.ndarray, index: int) -> np.ndarray: return y if self.raman_tilt_ is None: return self.raman_step_gain_[index] * y - return ifft(self.raman_tilt_[index] * fft(y)) + return apply_frequency_response(y, self.raman_tilt_[index]) def forward(self, x: np.ndarray) -> np.ndarray: # perform SSFM y = x # set by prepare(), which the Processor base always runs first assert (self.beta2 is not None and self.step_size is not None - and self.edfa_gain is not None and self.edfa_N_ase is not None) + and self.edfa_gain is not None and self.edfa_N_ase is not None + and self.transfer_ is not None + and self.transfer_index_ is not None) for num_span in range(self.N_spans): # perform for each span if self.use_only_linear: - y = apply_chromatic_dispersion(y, self.L_span, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=1) + y = apply_frequency_response(y, self.transfer_[0]) if self.raman_step_gain_ is not None: # no step loop here, so the whole profile applies at # once -- the span must stay transparent in this mode @@ -366,23 +518,25 @@ def forward(self, x: np.ndarray) -> np.ndarray: if self.raman_tilt_ is None: y = float(np.prod(self.raman_step_gain_)) * y else: - y = ifft(np.prod(self.raman_tilt_, axis=0) * fft(y)) + y = apply_frequency_response( + y, np.prod(self.raman_tilt_, axis=0)) else: for num_step in range(self.StPS): dz = self.step_size[num_step] # the two half-step Raman gains bracket the Kerr term, # like the two half-step dispersion operators + H = self.transfer_[self.transfer_index_[num_step]] if self.step_method == "symmetric": y = self._apply_raman(y, 2*num_step) - y = apply_chromatic_dispersion(y, dz/2, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=1) + y = apply_frequency_response(y, H) y = self._apply_kerr(y, dz) y = self._apply_raman(y, 2*num_step+1) - y = apply_chromatic_dispersion(y, dz/2, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=1) + y = apply_frequency_response(y, H) if self.step_method == "asymetric": y = self._apply_raman(self._apply_raman(y, 2*num_step), 2*num_step+1) y = self._apply_kerr(y, dz) - y = apply_chromatic_dispersion(y, dz, self.beta2, alpha_dB=self.fiber.alpha_dB, fs=self.fs, direction=1) + y = apply_frequency_response(y, H) # the ASE the distributed amplifier generated over this span, # already amplified by the gain downstream of where it was born diff --git a/src/comnumpy/optical/utils.py b/src/comnumpy/optical/utils.py index d9466af..ea62349 100644 --- a/src/comnumpy/optical/utils.py +++ b/src/comnumpy/optical/utils.py @@ -1,17 +1,29 @@ +import logging + import numpy as np from comnumpy._backend import fft, ifft, fftfreq # cupy-compatible (D3) -from typing import Optional +from typing import Optional, Tuple from comnumpy.exceptions import ShapeError from .constants import PLANCK_CONSTANT, OPTICAL_CARRIER_FREQUENCY +logger = logging.getLogger(__name__) + +# What `step_transfers` is keyed on: the FFT length, the sampling +# frequency, the dispersion, the attenuation, the direction, and the +# distinct step lengths. Named because it travels: the blocks hold one +# to know whether their table is still the right one. +TransferKey = Tuple[int, float, float, Optional[float], int, Tuple[float, ...]] +TransferTable = Tuple[np.ndarray, np.ndarray, TransferKey] + __all__ = [ "compute_beta2", "apply_chromatic_dispersion", "apply_kerr_nonlinearity", "compute_erbium_doped_fiber_amplifier_gain", "compute_erbium_doped_fiber_N_ase", "get_linear_step_size", "get_logarithmic_step_size", "itu_grid_frequency", "is_polarization_pair", - "manakov_kerr", "dbm_to_watt", "watt_to_dbm", + "manakov_kerr", "dbm_to_watt", "watt_to_dbm", "launch_amplitude", + "linear_step_transfer", "apply_frequency_response", "step_transfers", ] @@ -120,20 +132,237 @@ def apply_chromatic_dispersion(x: np.ndarray, z: float, beta2: float, optical receivers." Journal of lightwave technology 32.8 (2014): 1449-1456. """ - if alpha_dB: - alpha = (np.log(10)/10) * alpha_dB # convert dB to linear factor - gain = np.exp(-(alpha/2) * z * direction) # see text before equation 6 in https://arxiv.org/pdf/2010.14258.pdf - else: - gain = 1 + H = linear_step_transfer(x.shape[-1], z, beta2, fs=fs, alpha_dB=alpha_dB, + direction=direction, like=x) + return apply_frequency_response(x, H) + + +def linear_step_transfer(n_samples: int, z: float, beta2: float, *, + fs: float = 1, alpha_dB: Optional[float] = None, + direction: int = 1, + like: Optional[np.ndarray] = None) -> np.ndarray: + r"""Transfer function of one linear step of the split-step method. + + Signal Model + ------------ + The linear operator of the NLSE is diagonal in frequency, and it is + the whole of it -- attenuation and dispersion together: + + .. math:: + + H[k] = e^{-\frac{\alpha}{2} z d}\, + e^{j \frac{\beta_2}{2} z \omega_k^2 d} + + with :math:`\alpha = \frac{\ln 10}{10}\alpha_{dB}`, + :math:`\omega_k = 2\pi f_s f[k]` over the FFT grid, and + :math:`d = \pm 1` the propagation direction. The two factors are + the same operator :math:`D`, one real and constant across the band, + the other complex and frequency-dependent; separating them would + make a caller carry a scalar alongside the array for nothing. + + **Why this is a function of its own.** It depends only on + :math:`(n, z, \beta_2, f_s, \alpha_{dB}, d)` -- never on the + signal -- while a split-step loop applies it thousands of times. + Building it costs 0.42 ms at 12288 samples against 0.77 ms for the + FFT pair it serves, of which the complex exponential alone is 88 %, + so recomputing it per step is most of the cost of propagating. + :class:`~comnumpy.optical.links.FiberLink` builds one per distinct + step length in ``prepare`` and reuses them. + + Axes: *element-wise* -- returns a 1-D array of ``n_samples`` bins, + which broadcasts against any field shaped ``(..., n_samples)``. + + Parameters + ---------- + n_samples : int + Length of the FFT, i.e. of the field's last axis. + z : float + Step length in km. + beta2 : float + Group-velocity dispersion in ps^2/km. + fs : float, optional, keyword-only + Sampling frequency in Hz. Default 1. + alpha_dB : float, optional, keyword-only + Attenuation in dB/km. Default None, a lossless step. + direction : int, optional, keyword-only + 1 forward, -1 backward (back-propagation). Default 1. + like : np.ndarray, optional, keyword-only + Array whose type the FFT grid should follow, so that a cupy + field gets a cupy transfer function (D3). + + Returns + ------- + np.ndarray + Transfer function on the FFT grid, in FFT order. + + Examples + -------- + A lossless step is unit modulus -- dispersion moves phase, not + energy: + >>> H = linear_step_transfer(8, 10.0, -21.7, fs=100e9) + >>> print(np.round(np.abs(H), 12)) + [1. 1. 1. 1. 1. 1. 1. 1.] + + Backward is the conjugate of forward, which is what makes + back-propagation exact in the absence of noise: + + >>> back = linear_step_transfer(8, 10.0, -21.7, fs=100e9, direction=-1) + >>> print(np.round(np.max(np.abs(back - np.conj(H))), 12)) + 0.0 + """ beta2_s2_per_km = ((10**-12)**2) * beta2 # convert into s^2/km - NFFT = x.shape[-1] # the field may carry a polarization axis - w = (2*np.pi*fs)*fftfreq(NFFT, d=1, like=x) + w = (2*np.pi*fs)*fftfreq(n_samples, d=1, like=like) H = np.exp(1j * (beta2_s2_per_km/2) * z * (w**2) * direction) # see equation 4 - fftx = fft(x) - ffty = H * fftx - y = gain * ifft(ffty) - return y + if alpha_dB: + alpha = (np.log(10)/10) * alpha_dB # convert dB to linear factor + # see text before equation 6 in https://arxiv.org/pdf/2010.14258.pdf + H = H * np.exp(-(alpha/2) * z * direction) + return H + + +def step_transfers(n_samples: int, lengths: np.ndarray, *, beta2: float, + fs: float, alpha_dB: Optional[float] = None, + direction: int = 1, + previous: Optional[TransferTable] = None) -> TransferTable: + r"""Transfer functions of a split-step schedule, one per distinct length. + + Signal Model + ------------ + A split step applies :func:`linear_step_transfer` once per step, and + the operator depends only on the step length. A **linear** schedule + makes every step the same length, so one transfer function serves + all of them whatever ``StPS`` is; a **logarithmic** schedule makes + them all different and needs ``StPS``. Both cases are found the same + way -- by asking which lengths are distinct -- so the propagation + loop indexes a table and never branches on the schedule. + + Rebuilding costs a complex exponential per distinct length, which is + 88 % of the cost of a transfer function and 60 % of the FFT pair it + serves, so the table is keyed on everything it depends on and + returned unchanged when the key matches. A Monte-Carlo runs the same + block at the same length over and over, and then it is built once + for the whole sweep instead of once per pass. + + Axes: *element-wise* -- the table is ``(n_distinct, n_samples)`` and + each row broadcasts against a field shaped ``(..., n_samples)``. + + Parameters + ---------- + n_samples : int + Length of the FFT, i.e. of the field's last axis. + lengths : np.ndarray + Length of every step, in km, in the order they are applied. + beta2 : float, keyword-only + Group-velocity dispersion in ps^2/km. + fs : float, keyword-only + Sampling frequency in Hz. + alpha_dB : float, optional, keyword-only + Attenuation in dB/km. Default None, a lossless step. + direction : int, optional, keyword-only + 1 forward, -1 for back-propagation. Default 1. + previous : tuple, optional, keyword-only + The ``(table, index, key)`` returned by the last call, returned + as is when nothing it depends on has changed. + + Returns + ------- + tuple + ``(table, index, key)``: the ``(n_distinct, n_samples)`` + transfer functions, the row each step reads, and the key the + three were built from. + + Examples + -------- + A linear schedule collapses to one row, however many steps: + + >>> table, index, _ = step_transfers(64, np.full(50, 1.6), beta2=-21.7, + ... fs=100e9) + >>> print(table.shape, index[:4], index[-1]) + (1, 64) [0 0 0 0] 0 + + A logarithmic one keeps them apart: + + >>> table, index, _ = step_transfers(64, np.array([0.5, 1.5, 4.0]), + ... beta2=-21.7, fs=100e9) + >>> print(table.shape, index) + (3, 64) [0 1 2] + + Asking again with the same schedule returns the same object, which + is what keeps a sweep from rebuilding it: + + >>> again = step_transfers(64, np.array([0.5, 1.5, 4.0]), beta2=-21.7, + ... fs=100e9, previous=(table, index, _)) + >>> print(again[0] is table) + True + """ + lengths = np.asarray(lengths, dtype=float) + distinct, index = np.unique(lengths, return_inverse=True) + key = (n_samples, fs, beta2, alpha_dB, direction, + tuple(distinct.tolist())) + if previous is not None and previous[2] == key: + return previous + # A logarithmic schedule stores one array per step. That is the price + # of steps that differ, and logarithmic steps exist so that the step + # count stays small: 0.8 MB at StPS=4, 590 MB at StPS=500 -- and at + # five hundred steps a linear schedule converges just as well and + # caches a single row. + megabytes = distinct.size * n_samples * 16 / 1e6 + if megabytes > 64: + logger.warning( + "caching %d distinct step lengths at %d samples needs %.0f MB; " + "a linear schedule needs one row whatever the step count.", + distinct.size, n_samples, megabytes) + rows = [] + for z in distinct: + rows.append(linear_step_transfer(n_samples, float(z), beta2, fs=fs, + alpha_dB=alpha_dB, + direction=direction)) + return np.stack(rows), index, key + + +def apply_frequency_response(x: np.ndarray, H: np.ndarray, *, + axis: int = -1) -> np.ndarray: + r"""Filter a signal by a transfer function, in the frequency domain. + + Signal Model + ------------ + .. math:: + + y = \mathcal{F}^{-1}\left\{ H \odot \mathcal{F}\{x\} \right\} + + The multiplication is **periodic**: no zero padding, so the result + is a circular convolution and ``y`` has the length of ``x``. That is + what a split-step wants, and what a channel mask or a Raman tilt + wants; a filter whose impulse response must not wrap around -- + pulse shaping -- pads instead, which is a different operation and + lives in :mod:`comnumpy.core.filters`. + + Axes: *declared axis* -- filters along ``axis`` and broadcasts over + the rest, so a ``(B, P, N)`` field is filtered by a length-``N`` + response without a loop. + + Parameters + ---------- + x : np.ndarray + Signal to filter. + H : np.ndarray + Transfer function on the FFT grid, broadcastable against ``x``. + axis : int, optional, keyword-only + Axis to filter along. Default -1. + + Returns + ------- + np.ndarray + Filtered signal, same shape as ``x``. + + Examples + -------- + >>> x = np.array([1.0, 0.0, 0.0, 0.0]) + >>> print(np.round(np.real(apply_frequency_response(x, np.ones(4))), 12)) + [1. 0. 0. 0.] + """ + return ifft(H * fft(x, axis=axis), axis=axis) def apply_kerr_nonlinearity(x: np.ndarray, z: float, gamma: float, @@ -367,6 +596,14 @@ def is_polarization_pair(x: np.ndarray, block: str) -> bool: refused, because a pointwise Kerr step applied row by row would describe parallel fibres, not one fibre carrying several signals. + **The leading axes are a batch.** Only ``-2`` and ``-1`` are read, + so ``(B, P, N)`` propagates ``B`` independent realizations in one + call and every block downstream broadcasts over them -- the same + rule the MIMO blocks follow on their antenna axis (D2). That is + also why the polarization axis must be written even when there is + one of them: ``(B, N)`` is indistinguishable from ``(P, N)``, so a + batch of single-polarization fields is ``(B, 1, N)``. + Parameters ---------- x : np.ndarray @@ -394,8 +631,12 @@ def is_polarization_pair(x: np.ndarray, block: str) -> bool: f"(..., P, N) with P in {{1, 2}}; got {x.shape}, i.e. " f"{polarizations} on the polarization axis. A pointwise Kerr step " f"row by row would describe {polarizations} separate fibres, with " - f"no XPM and no FWM between them. Multiplex WDM channels into one " - f"field first with comnumpy.optical.WDMMultiplexer (decision D44).") + f"no XPM and no FWM between them. Two shapes are confused with " + f"this one. For {polarizations} WDM channels, multiplex them into " + f"one field first with comnumpy.optical.WDMMultiplexer (D44). For " + f"{polarizations} independent realizations propagated at once, " + f"write the polarization axis: " + f"({polarizations}, 1, {x.shape[-1]}).") def manakov_kerr(x: np.ndarray, gamma: float, @@ -526,3 +767,71 @@ def watt_to_dbm(power_W: np.ndarray | float) -> np.ndarray | float: f"positive power in watts; got a minimum of {float(np.min(power))}.") result = 10 * np.log10(power / 1e-3) return result if np.ndim(power_W) else float(result) + + +def launch_amplitude(power_W: np.ndarray | float, *, + polarizations: int = 1) -> np.ndarray | float: + r"""Field amplitude that launches a given optical power. + + Signal Model + ------------ + An :class:`~comnumpy.core.processors.Amplifier` multiplies the + **field**, and a power is the squared modulus of it, so a launch + power is set through a square root: + + .. math:: + + a = \sqrt{\frac{P}{P_{\mathrm{pol}}}} + + The division is the point. A dual-polarization signal carries the + channel power split between its two polarizations, so each one is + launched at :math:`\sqrt{P/2}` -- and that factor, written by hand + at each call site, is the same one that makes an ASE budget or a + nonlinear coefficient disagree by 3 dB (see + :meth:`~comnumpy.optical.links.FiberLink.budget`). + + Axes: *element-wise* -- an array of powers gives an array of + amplitudes, which is what a launch-power sweep needs. + + Parameters + ---------- + power_W : float or np.ndarray + Channel power in watts, summed over the polarizations. + polarizations : int, optional, keyword-only + 1 (default) for a single-polarization field, 2 for a + polarization pair. + + Returns + ------- + float or np.ndarray + Amplitude gain to give an ``Amplifier``. + + Raises + ------ + ValueError + If a power is negative, or if ``polarizations`` is neither 1 + nor 2. + + Examples + -------- + >>> print(f"{launch_amplitude(dbm_to_watt(0.0)):.6f}") + 0.031623 + + Two polarizations share the channel power, so each carries half of + it -- 3 dB down, which is :math:`\sqrt{2}` in amplitude: + + >>> single = launch_amplitude(1e-3) + >>> pair = launch_amplitude(1e-3, polarizations=2) + >>> print(f"{single / pair:.6f}") + 1.414214 + """ + if polarizations not in (1, 2): + raise ValueError( + f"a fibre carries one or two polarizations, got {polarizations}.") + power = np.asarray(power_W, dtype=float) + if np.any(power < 0): + raise ValueError( + f"a launch power is not negative; got a minimum of " + f"{float(np.min(power))} W.") + result = np.sqrt(power / polarizations) + return result if np.ndim(power_W) else float(result) diff --git a/src/comnumpy/optical/wdm.py b/src/comnumpy/optical/wdm.py index 4518761..b3576f4 100644 --- a/src/comnumpy/optical/wdm.py +++ b/src/comnumpy/optical/wdm.py @@ -32,6 +32,8 @@ import numpy as np +from .utils import apply_frequency_response + from comnumpy.core.generics import Processor from comnumpy.exceptions import ShapeError @@ -652,4 +654,4 @@ def forward(self, x: np.ndarray) -> np.ndarray: mixed = x[..., None, :] * np.conj(phasors) else: mixed = x * np.conj(phasors) - return np.fft.ifft(np.fft.fft(mixed, axis=-1) * self.mask_, axis=-1) + return apply_frequency_response(mixed, self.mask_) diff --git a/tests/core/test_error_counter.py b/tests/core/test_error_counter.py new file mode 100644 index 0000000..026b7f5 --- /dev/null +++ b/tests/core/test_error_counter.py @@ -0,0 +1,106 @@ +"""Counting errors across the trials of a Monte-Carlo point. + +Four example scripts accumulate errors over trials, each with its own +loop, and three of them wrote ``+= compute_ser(...) / n_trials``. That +averages *rates*, which is only the same answer when every trial has the +same length, and it discards the count the point's credibility rests on. +The counter keeps the totals; the loop stays where it was. +""" +import unittest + +import numpy as np + +from comnumpy.core.metrics import ErrorCounter, compute_ber, compute_ser +from comnumpy.exceptions import NotFittedError + + +class TestErrorCounter(unittest.TestCase): + + def test_one_trial_agrees_with_compute_ser(self): + rng = np.random.default_rng(0) + sent = rng.integers(0, 16, 500) + received = sent.copy() + received[rng.choice(500, 20, replace=False)] ^= 1 + counter = ErrorCounter().update(sent, received) + self.assertAlmostEqual(counter.rate, compute_ser(sent, received)) + + def test_bits_agree_with_compute_ber(self): + rng = np.random.default_rng(1) + sent = rng.integers(0, 16, 500) + received = rng.integers(0, 16, 500) + counter = ErrorCounter(width=4).update(sent, received) + self.assertAlmostEqual(counter.rate, compute_ber(sent, received, 4)) + + def test_it_is_a_ratio_of_totals_not_a_mean_of_rates(self): + """The regression this exists for: unequal trials. + + One error in ten symbols, then none in a thousand. The honest + rate is 1/1010; averaging the two rates gives 1/20, fifty times + larger, and that is what a loop dividing by the trial count + computes. + """ + counter = ErrorCounter() + counter.update(np.zeros(10, dtype=int), + np.array([1] + [0] * 9)) + counter.update(np.zeros(1000, dtype=int), np.zeros(1000, dtype=int)) + self.assertAlmostEqual(counter.rate, 1 / 1010) + mean_of_rates = (0.1 + 0.0) / 2 + self.assertNotAlmostEqual(counter.rate, mean_of_rates) + + def test_equal_trials_reproduce_the_loop_being_replaced(self): + """Where the examples stand today, the answer must not move.""" + rng = np.random.default_rng(2) + counter = ErrorCounter() + by_hand = 0.0 + n_trials = 4 + for _ in range(n_trials): + sent = rng.integers(0, 4, 250) + received = rng.integers(0, 4, 250) + counter.update(sent, received) + by_hand += compute_ser(sent, received) / n_trials + self.assertAlmostEqual(counter.rate, by_hand) + + def test_counts_are_reported_beside_the_rate(self): + counter = ErrorCounter() + counter.update(np.zeros(100, dtype=int), np.zeros(100, dtype=int)) + counter.update(np.zeros(100, dtype=int), + np.array([1, 1] + [0] * 98)) + self.assertEqual((counter.n_errors, counter.n_symbols, + counter.n_trials), (2, 200, 2)) + + def test_a_shorter_estimate_sets_the_length(self): + """A chain with a filter delay returns fewer symbols than it got.""" + counter = ErrorCounter().update(np.zeros(100, dtype=int), + np.zeros(90, dtype=int)) + self.assertEqual(counter.n_symbols, 90) + + def test_update_chains_and_reset_forgets(self): + counter = ErrorCounter(width=2) + counter.update(np.zeros(8, dtype=int), np.zeros(8, dtype=int)) + self.assertIs(counter.reset(), counter) + self.assertEqual((counter.n_errors, counter.n_trials), (0, 0)) + self.assertEqual(counter.width, 2) # reset keeps the unit + + def test_a_rate_before_any_trial_is_refused(self): + """An empty average is a NaN waiting to reach a figure (D38).""" + with self.assertRaises(NotFittedError) as ctx: + _ = ErrorCounter().rate + self.assertIn("update", str(ctx.exception)) + + def test_str_says_counts_and_rate(self): + counter = ErrorCounter() + self.assertIn("nothing counted", str(counter)) + counter.update(np.zeros(1000, dtype=int), + np.array([1] + [0] * 999)) + self.assertEqual(str(counter), + "ErrorCounter(1 error / 1000 symbols, rate 1.00e-03)") + + def test_refusals(self): + with self.assertRaises(ValueError): + ErrorCounter(width=0) + with self.assertRaises(ValueError): + ErrorCounter().update(np.zeros(0, dtype=int), np.zeros(5, dtype=int)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/core/test_set_params.py b/tests/core/test_set_params.py index 0d16816..cbabadc 100644 --- a/tests/core/test_set_params.py +++ b/tests/core/test_set_params.py @@ -6,7 +6,7 @@ The double underscore is the same address spelled as an identifier, which is what scikit-learn uses, so it can be a plain keyword argument. -The dotted form stays: `sweep` builds its addresses as strings, and a +The dotted form stays: `monte_carlo` builds its addresses as strings, and a string is what you want when the address is computed rather than typed. """ import unittest diff --git a/tests/optical/test_batch_propagation.py b/tests/optical/test_batch_propagation.py new file mode 100644 index 0000000..27c4706 --- /dev/null +++ b/tests/optical/test_batch_propagation.py @@ -0,0 +1,128 @@ +"""Propagating several realizations in one call. + +`is_polarization_pair` reads axes -2 and -1 only, so the leading axes of +a field are free and `FiberLink` broadcasts over them. That has always +been true and nothing pinned it, which for a capability is the same as +not having it: the next refactor of `prepare()` could take it away and +every test would still pass. + +The claim is stronger than "it runs". A batch must give *the same +answer* as the loop it replaces, span for span and step for step, or it +is a different simulation wearing the same name. With the amplifier +noise off that is bit-for-bit; with it on it cannot be, because one draw +of `B x N` samples is not `B` draws of `N`, and the test says which of +the two it checks. +""" +import unittest + +import numpy as np + +from comnumpy.exceptions import ShapeError +from comnumpy.optical.dbp import DBP +from comnumpy.optical.links import FiberLink + +N = 512 + + +def field(shape, seed=0): + rng = np.random.default_rng(seed) + return 1e-3 * (rng.normal(size=shape) + 1j * rng.normal(size=shape)) + + +class TestBatchIsTheSameSimulation(unittest.TestCase): + + def link(self, **kwargs): + params = dict(L_span=80.0, StPS=8, fs=10e9, noise_scaling=0.0) + params.update(kwargs) + return FiberLink(3, **params) + + def test_a_batch_of_one_polarization_equals_the_loop(self): + """Bit for bit, which is what makes the batch a shortcut.""" + rows = [field(N, seed=index) for index in range(4)] + batch = np.stack(rows)[:, None, :] # (4, 1, N) + link = self.link() + together = link(batch.copy()) + self.assertEqual(together.shape, (4, 1, N)) + for index, row in enumerate(rows): + with self.subTest(item=index): + np.testing.assert_array_equal(together[index, 0], + link(row.copy())) + + def test_a_batch_of_pairs_equals_the_loop(self): + """The Manakov coupling stays inside an item, never across items.""" + pairs = [field((2, N), seed=index) for index in range(3)] + link = self.link() + together = link(np.stack(pairs).copy()) + self.assertEqual(together.shape, (3, 2, N)) + for index, pair in enumerate(pairs): + with self.subTest(item=index): + np.testing.assert_array_equal(together[index], + link(pair.copy())) + + def test_one_item_of_a_batch_equals_the_bare_signal(self): + """`(1, 1, N)` and `(N,)` are the same field written two ways.""" + signal = field(N) + link = self.link() + np.testing.assert_array_equal( + link(signal.reshape(1, 1, N).copy())[0, 0], link(signal.copy())) + + def test_the_kerr_term_does_not_leak_between_items(self): + """The nonlinearity is pointwise, so a loud item must not be felt. + + Two realizations, one of them a hundred times stronger. If the + intensity were summed over the batch instead of over the + polarizations, the quiet one would come back with the loud + one's phase. + """ + quiet, loud = field(N, seed=1), 100 * field(N, seed=2) + link = self.link() + batch = np.stack([quiet, loud])[:, None, :] + together = link(batch.copy()) + np.testing.assert_array_equal(together[0, 0], link(quiet.copy())) + + def test_dbp_batches_the_same_way(self): + """A receiver has to accept what the link produced.""" + rows = [field(N, seed=index) for index in range(3)] + dbp = DBP(3, L_span=80.0, StPS=8, fs=10e9) + together = dbp(np.stack(rows)[:, None, :].copy()) + for index, row in enumerate(rows): + with self.subTest(item=index): + np.testing.assert_array_equal(together[index, 0], + dbp(row.copy())) + + def test_the_round_trip_survives_a_batch(self): + """Propagate a batch, back-propagate it, recover every item.""" + rows = np.stack([field(N, seed=index) for index in range(3)])[:, None, :] + link = self.link(StPS=200) + back = DBP(3, L_span=80.0, StPS=200, fs=10e9) + recovered = back(link(rows.copy())) + np.testing.assert_allclose(recovered, rows, rtol=0, atol=1e-9) + + def test_noise_is_drawn_for_the_whole_batch(self): + """With ASE on the batch is not the loop, and that is not a bug. + + One draw of `B x N` Gaussian samples is not `B` draws of `N`, so + the realizations differ. What must hold is that each item still + carries the noise the link budgets for it. + """ + link = self.link(noise_scaling=1.0, NF_dB=5.0, seed=0) + rows = np.zeros((64, 1, 4096), dtype=complex) + noise = link(rows) + per_item = np.var(noise, axis=-1).ravel() + expected = link.budget(link.fs)["ase_power_W"] + self.assertAlmostEqual(float(np.mean(per_item)) / expected, 1.0, + delta=0.05) + # every item got noise, none was left silent + self.assertGreater(float(np.min(per_item)) / expected, 0.5) + + def test_a_two_dimensional_batch_is_refused_and_says_how(self): + """`(B, N)` is what a user reaches for first; the message helps.""" + with self.assertRaises(ShapeError) as ctx: + self.link()(np.zeros((32, N), dtype=complex)) + message = str(ctx.exception) + self.assertIn(f"(32, 1, {N})", message) # the shape to write + self.assertIn("WDMMultiplexer", message) # the other reading + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/optical/test_link_budget.py b/tests/optical/test_link_budget.py new file mode 100644 index 0000000..c703603 --- /dev/null +++ b/tests/optical/test_link_budget.py @@ -0,0 +1,146 @@ +"""The link budgets its own amplifier noise. + +Three tutorials used to write the same product by hand -- span count, +spectral density, bandwidth, and a factor for the polarizations -- with +three different spellings of the last term. That factor is the one whose +convention made two pages of this documentation appear to disagree by +several decibels, so the formula now lives next to the block that +generates the noise, and these tests hold it there. +""" +import unittest + +import numpy as np + +from comnumpy.optical.fiber import FiberSpec +from comnumpy.optical.links import FiberLink +from comnumpy.optical.utils import (compute_erbium_doped_fiber_N_ase, + dbm_to_watt, launch_amplitude, + watt_to_dbm) + + +class TestBudget(unittest.TestCase): + + def link(self, **kwargs): + params = dict(L_span=100.0, NF_dB=5.0, fs=192e9) + params.update(kwargs) + n_spans = params.pop("N_spans", 10) + return FiberLink(n_spans, **params) + + def test_matches_the_formula_it_replaces(self): + """The three examples wrote this product; it must not move.""" + link = self.link() + density = compute_erbium_doped_fiber_N_ase( + link.fiber.alpha_dB, link.L_span, link.NF_dB, + nu=link.fiber.carrier_frequency_Hz) + by_hand = link.N_spans * density * 32e9 + self.assertAlmostEqual(link.budget(32e9)["ase_power_W"], by_hand, + delta=1e-18) + + def test_the_second_polarization_doubles_it(self): + """The 3 dB that made two pages look like they disagreed.""" + link = self.link() + single = link.budget(32e9)["ase_power_dBm"] + pair = link.budget(32e9, polarizations=2)["ase_power_dBm"] + self.assertAlmostEqual(pair - single, 3.0103, places=3) + + def test_published_values_of_the_two_tutorials(self): + """Both pages quote a budget; both come from here now.""" + # the back-propagation page: 10 x 100 km, NF 5 dB, one polarization + self.assertAlmostEqual(self.link().budget(32e9)["ase_power_dBm"], + -21.88, places=2) + # the GN-model page: 5 x 100 km, NF 6 dB, both polarizations + smf = FiberSpec(0.2, gamma=1.3, cd_coefficient=17.0, + wavelength_nm=1550.0) + page = self.link(N_spans=5, NF_dB=6.0, fiber=smf) + self.assertAlmostEqual(page.budget(32e9, polarizations=2)["ase_power_dBm"], + -20.88, places=2) + + def test_it_scales_with_the_spans_and_the_bandwidth(self): + one = self.link(N_spans=1).budget(32e9)["ase_power_W"] + self.assertAlmostEqual(self.link(N_spans=10).budget(32e9)["ase_power_W"], + 10 * one, delta=1e-18) + self.assertAlmostEqual(self.link(N_spans=1).budget(64e9)["ase_power_W"], + 2 * one, delta=1e-18) + + def test_it_agrees_with_what_the_link_actually_adds(self): + """`prepare` and `budget` read one formula, so they cannot drift. + + The block stores its per-span variance over the *simulated* + bandwidth; the budget quotes the accumulated power over the + bandwidth a receiver keeps. Same density, two bandwidths. + """ + link = self.link() + link(np.zeros(64, dtype=complex)) # prepare() runs on the call + density = link.budget(32e9)["ase_density_W_per_Hz"] + self.assertAlmostEqual(link.edfa_N_ase, link.fs * density, delta=1e-18) + + def test_noise_scaling_is_honoured(self): + """A link with its noise off budgets nothing, and says so.""" + budget = self.link(noise_scaling=0.0).budget(32e9) + self.assertEqual(budget["ase_power_W"], 0.0) + self.assertEqual(budget["ase_power_dBm"], -np.inf) + + def test_the_link_really_adds_what_it_budgets(self): + """The claim the tutorials make: this predicts the simulation. + + A transparent span returns the power it was given, so whatever + power the output carries above the input is the noise the + amplifiers put there -- measured over the whole simulated band, + which is the bandwidth to ask the budget for. Comparing powers + rather than a sample-wise difference keeps the dispersion out of + it: it reshapes the waveform and conserves the power. + """ + from comnumpy.core import Sequential + from comnumpy.core.generators import GaussianGenerator + + link = FiberLink(5, L_span=80.0, NF_dB=5.0, fs=64e9, + use_only_linear=True, name="link") + chain = Sequential([GaussianGenerator(dbm_to_watt(0.0), name="tx"), + link], taps=["tx"]).seed(0) + y = chain(400000) + measured = float(np.var(y) - np.var(chain.tap("tx"))) + predicted = link.budget(link.fs)["ase_power_W"] + self.assertAlmostEqual(10 * np.log10(measured / predicted), 0.0, + delta=0.1) + + def test_refusals(self): + link = self.link() + for bandwidth in (0.0, -1.0): + with self.subTest(bandwidth=bandwidth): + with self.assertRaises(ValueError) as ctx: + link.budget(bandwidth) + self.assertIn("symbol rate", str(ctx.exception)) + with self.assertRaises(ValueError): + link.budget(32e9, polarizations=3) + + +class TestLaunchAmplitude(unittest.TestCase): + + def test_it_is_the_square_root_of_the_power(self): + self.assertAlmostEqual(launch_amplitude(1e-4) ** 2, 1e-4) + + def test_the_pair_splits_the_channel_power(self): + """Each polarization carries half, i.e. 3 dB less.""" + pair = launch_amplitude(1e-3, polarizations=2) + self.assertAlmostEqual(watt_to_dbm(pair ** 2), watt_to_dbm(5e-4)) + + def test_an_array_of_powers_gives_an_array(self): + """A launch-power sweep passes the whole axis at once.""" + powers = np.array([1e-4, 1e-3, 1e-2]) + amplitudes = launch_amplitude(powers) + self.assertEqual(amplitudes.shape, powers.shape) + np.testing.assert_allclose(amplitudes ** 2, powers) + + def test_a_scalar_stays_a_scalar(self): + """Scalar in, scalar out -- the library's rule everywhere else.""" + self.assertIsInstance(launch_amplitude(1e-3), float) + + def test_refusals(self): + with self.assertRaises(ValueError): + launch_amplitude(-1e-3) + with self.assertRaises(ValueError): + launch_amplitude(1e-3, polarizations=0) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/optical/test_step_transfer.py b/tests/optical/test_step_transfer.py new file mode 100644 index 0000000..fcf1266 --- /dev/null +++ b/tests/optical/test_step_transfer.py @@ -0,0 +1,155 @@ +"""The linear step, split into what it is and what it does. + +`apply_chromatic_dispersion` built its transfer function on every call +and threw it away: an `fftfreq`, a `w**2` and a length-N complex +exponential, a thousand times per propagation, for something that +depends only on `(n, dz, beta2, fs, alpha_dB, direction)`. Splitting it +into `linear_step_transfer` (what the step is) and +`apply_frequency_response` (what it does to a signal) lets the blocks +build the first once and call the second in the loop. + +The composed function keeps its name and its behaviour, so nothing +downstream had to move; these tests hold both halves and their +composition, and pin the property that makes the cache legitimate -- +that a linear schedule has exactly one distinct step length. +""" +import unittest + +import numpy as np + +from comnumpy.optical.utils import (apply_chromatic_dispersion, + apply_frequency_response, + linear_step_transfer, step_transfers) + +N = 256 +BETA2, FS, ALPHA = -21.7, 64e9, 0.2 + + +class TestLinearStepTransfer(unittest.TestCase): + + def test_lossless_dispersion_moves_phase_not_energy(self): + H = linear_step_transfer(N, 20.0, BETA2, fs=FS) + np.testing.assert_allclose(np.abs(H), 1.0, atol=1e-12) + + def test_attenuation_is_folded_in(self): + """The step is the whole linear operator, loss included.""" + H = linear_step_transfer(N, 20.0, BETA2, fs=FS, alpha_dB=ALPHA) + expected = 10 ** (-ALPHA * 20.0 / 20) # amplitude, not power + np.testing.assert_allclose(np.abs(H), expected, rtol=1e-12) + + def test_backward_is_the_conjugate_of_forward(self): + """What makes noiseless back-propagation exact.""" + forward = linear_step_transfer(N, 20.0, BETA2, fs=FS) + backward = linear_step_transfer(N, 20.0, BETA2, fs=FS, direction=-1) + np.testing.assert_allclose(backward, np.conj(forward), atol=1e-15) + + def test_steps_compose_by_adding_their_lengths(self): + """Two half-steps are one step -- the split-step's own premise.""" + half = linear_step_transfer(N, 10.0, BETA2, fs=FS, alpha_dB=ALPHA) + whole = linear_step_transfer(N, 20.0, BETA2, fs=FS, alpha_dB=ALPHA) + np.testing.assert_allclose(half * half, whole, rtol=1e-12) + + def test_zero_dispersion_and_zero_length_are_identity(self): + np.testing.assert_allclose(linear_step_transfer(N, 0.0, BETA2, fs=FS), + np.ones(N), atol=1e-15) + + +class TestApplyFrequencyResponse(unittest.TestCase): + + def test_a_flat_response_returns_the_signal(self): + x = np.random.default_rng(0).normal(size=N) + 0j + np.testing.assert_allclose(apply_frequency_response(x, np.ones(N)), x, + atol=1e-12) + + def test_it_broadcasts_over_the_leading_axes(self): + """A batch of polarization pairs is filtered without a loop.""" + rng = np.random.default_rng(1) + x = rng.normal(size=(4, 2, N)) + 1j * rng.normal(size=(4, 2, N)) + H = linear_step_transfer(N, 20.0, BETA2, fs=FS) + together = apply_frequency_response(x, H) + self.assertEqual(together.shape, x.shape) + for item in range(4): + for pol in range(2): + with self.subTest(item=item, pol=pol): + np.testing.assert_allclose( + together[item, pol], + apply_frequency_response(x[item, pol], H), atol=1e-15) + + +class TestTheCompositionIsTheOldFunction(unittest.TestCase): + + def test_apply_chromatic_dispersion_is_the_two_halves(self): + rng = np.random.default_rng(2) + x = rng.normal(size=N) + 1j * rng.normal(size=N) + for z, direction in ((2.0, 1), (80.0, -1)): + with self.subTest(z=z, direction=direction): + H = linear_step_transfer(N, z, BETA2, fs=FS, alpha_dB=ALPHA, + direction=direction) + np.testing.assert_allclose( + apply_chromatic_dispersion(x, z, BETA2, alpha_dB=ALPHA, + fs=FS, direction=direction), + apply_frequency_response(x, H), rtol=1e-12) + + +class TestStepTransfers(unittest.TestCase): + + def test_a_linear_schedule_needs_one_row(self): + """The property the cache rests on: equal steps, one operator.""" + table, index, _ = step_transfers(N, np.full(200, 0.4), beta2=BETA2, + fs=FS) + self.assertEqual(table.shape, (1, N)) + self.assertEqual(set(index.tolist()), {0}) + + def test_a_logarithmic_schedule_needs_one_row_each(self): + lengths = np.array([0.2, 0.6, 1.8, 5.4]) + table, index, _ = step_transfers(N, lengths, beta2=BETA2, fs=FS) + self.assertEqual(table.shape, (4, N)) + np.testing.assert_array_equal(index, [0, 1, 2, 3]) + for row, z in enumerate(lengths): + with self.subTest(z=z): + np.testing.assert_allclose( + table[row], linear_step_transfer(N, z, BETA2, fs=FS), + atol=1e-15) + + def test_the_rows_follow_the_order_the_steps_are_applied(self): + """Unsorted lengths must still map each step to its own row.""" + lengths = np.array([5.0, 1.0, 5.0, 2.0]) + table, index, _ = step_transfers(N, lengths, beta2=BETA2, fs=FS) + self.assertEqual(table.shape, (3, N)) # three distinct + for step, z in enumerate(lengths): + with self.subTest(step=step): + np.testing.assert_allclose( + table[index[step]], + linear_step_transfer(N, z, BETA2, fs=FS), atol=1e-15) + + def test_an_unchanged_schedule_is_not_rebuilt(self): + """What makes a Monte-Carlo pay for this once instead of per pass.""" + previous = step_transfers(N, np.full(8, 0.4), beta2=BETA2, fs=FS) + again = step_transfers(N, np.full(8, 0.4), beta2=BETA2, fs=FS, + previous=previous) + self.assertIs(again[0], previous[0]) + + def test_every_dependency_invalidates_the_cache(self): + lengths = np.full(8, 0.4) + previous = step_transfers(N, lengths, beta2=BETA2, fs=FS) + changes = { + "n_samples": dict(n_samples=2 * N), + "fs": dict(fs=2 * FS), + "beta2": dict(beta2=2 * BETA2), + "alpha_dB": dict(alpha_dB=ALPHA), + "direction": dict(direction=-1), + "lengths": dict(lengths=np.full(8, 0.8)), + } + for name, override in changes.items(): + with self.subTest(changed=name): + params = dict(n_samples=N, lengths=lengths, beta2=BETA2, + fs=FS) + params.update(override) + rebuilt = step_transfers(params.pop("n_samples"), + params.pop("lengths"), + previous=previous, **params) + self.assertIsNot(rebuilt[0], previous[0]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_docs_coverage.py b/tests/test_docs_coverage.py index 9819935..270d927 100644 --- a/tests/test_docs_coverage.py +++ b/tests/test_docs_coverage.py @@ -3,7 +3,7 @@ The documentation described 25 of the 46 modules in ``src/``. Everything added after the first milestone -- capacity, fading, WDM, Raman, the whole ``fec`` package, allocation, frames, sequences, serialization, -sweep, the exception hierarchy -- existed with full course-material +monte_carlo, the exception hierarchy -- existed with full course-material docstrings that no generated page ever showed. Sphinx does not notice: a module with no ``automodule`` directive is not an error, it is simply absent, which is the worst way for documentation to be wrong. diff --git a/tests/test_examples_run.py b/tests/test_examples_run.py index 61b2eb2..9527413 100644 --- a/tests/test_examples_run.py +++ b/tests/test_examples_run.py @@ -63,12 +63,12 @@ "compensator -- 227 s measured 2026-08-09", "optical/NLI_simulation.py": "eight launch powers x four trials x two fibre configurations of " - "split-step propagation, then six receivers on each -- ~1400 s " - "measured 2026-08-12", + "split-step propagation, then six receivers on each -- 144 s " + "measured 2026-08-13 at StPS=50, with the step transfer cached", "optical/one_shot_NLI.py": "the whole chain re-run at six span counts, each once nonlinear and " "once linear, then twice more for the profile and the receiver " - "comparison -- 104 s measured 2026-08-12", + "comparison -- 13 s measured 2026-08-13", "simple/one_shot_srrc_awgn.py": "16 001-tap SRRC filter at oversampling 8 -- 113 s wall / 316 s " "CPU measured 2026-08-09", diff --git a/tests/test_monte_carlo.py b/tests/test_monte_carlo.py new file mode 100644 index 0000000..208dffc --- /dev/null +++ b/tests/test_monte_carlo.py @@ -0,0 +1,52 @@ +"""monte_carlo(): reconfigure, reseed, run, collect (decision D35).""" +import unittest + +import numpy as np + +from comnumpy import AWGN, Sequential, SymbolDemapper, \ + SymbolGenerator, SymbolMapper, compute_ser, get_alphabet, monte_carlo + + +class TestMonteCarlo(unittest.TestCase): + + def build(self): + alphabet = get_alphabet("QAM", 4) + return Sequential([ + SymbolGenerator(4, name="tx"), + SymbolMapper(alphabet), + AWGN(snr_dB=0, name="noise"), + SymbolDemapper(alphabet), + ]) + + def test_results_aligned_and_monotonic(self): + values = np.arange(0, 14, 4) + out = monte_carlo(self.build(), "noise.snr_dB", values, + {"ser": compute_ser}, 20_000, reference="tx", seed=1) + self.assertEqual(out["ser"].shape, values.shape) + self.assertTrue(np.all(np.diff(out["ser"]) < 0)) # SER falls with SNR + + def test_reproducible_across_calls(self): + values = [0, 6] + a = monte_carlo(self.build(), "noise.snr_dB", values, + {"ser": compute_ser}, 5_000, reference="tx", seed=3) + b = monte_carlo(self.build(), "noise.snr_dB", values, + {"ser": compute_ser}, 5_000, reference="tx", seed=3) + np.testing.assert_array_equal(a["ser"], b["ser"]) + + def test_multi_param_zip(self): + chain = self.build() + out = monte_carlo(chain, ["noise.snr_dB", "noise.name"], + [(0, "noise"), (10, "noise")], + {"snr_set": lambda t, y: chain["noise"].snr_dB}, + 1_000, reference="tx", seed=2) + np.testing.assert_array_equal(out["snr_set"], [0, 10]) + + def test_metric_without_reference(self): + out = monte_carlo(self.build(), "noise.snr_dB", [0, 20], + {"power": lambda y: float(np.mean(np.abs(y) ** 2))}, + 2_000, seed=4) + self.assertEqual(len(out["power"]), 2) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_sweep_parallel.py b/tests/test_monte_carlo_parallel.py similarity index 65% rename from tests/test_sweep_parallel.py rename to tests/test_monte_carlo_parallel.py index d047f54..25badd6 100644 --- a/tests/test_sweep_parallel.py +++ b/tests/test_monte_carlo_parallel.py @@ -13,7 +13,7 @@ import numpy as np -from comnumpy import sweep +from comnumpy import monte_carlo from comnumpy.core import Sequential from comnumpy.core.channels import AWGN from comnumpy.core.generators import SymbolGenerator @@ -42,51 +42,51 @@ def mean_power(y): class TestParallelIsTheSameSweep(unittest.TestCase): def test_the_curve_is_identical_value_for_value(self): - serial = sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"ser": compute_ser}, 2000, reference="tx", seed=5) - parallel = sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"ser": compute_ser}, 2000, reference="tx", seed=5, - n_jobs=2) + serial = monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"ser": compute_ser}, 2000, reference="tx", seed=5) + parallel = monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"ser": compute_ser}, 2000, reference="tx", seed=5, + n_jobs=2) np.testing.assert_array_equal(parallel["ser"], serial["ser"]) def test_the_order_of_the_points_is_preserved(self): """Workers finish out of order; the arrays must not.""" - serial = sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"power": mean_power}, 500, seed=1) - parallel = sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"power": mean_power}, 500, seed=1, n_jobs=-1) + serial = monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"power": mean_power}, 500, seed=1) + parallel = monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"power": mean_power}, 500, seed=1, n_jobs=-1) np.testing.assert_array_equal(parallel["power"], serial["power"]) self.assertEqual(parallel["power"].shape, (len(SNR_DB),)) def test_more_workers_than_points_is_not_an_error(self): - out = sweep(build_chain(), "noise.snr_dB", [3.0, 6.0], - {"power": mean_power}, 200, seed=2, n_jobs=8) + out = monte_carlo(build_chain(), "noise.snr_dB", [3.0, 6.0], + {"power": mean_power}, 200, seed=2, n_jobs=8) self.assertEqual(out["power"].shape, (2,)) class TestWhatCannotCrossIntoAWorker(unittest.TestCase): - def test_a_sweep_without_a_seed_is_refused(self): + def test_a_run_without_a_seed_is_refused(self): """Without a seed, what a point draws would depend on the worker.""" with self.assertRaises(ValueError) as ctx: - sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"power": mean_power}, 100, n_jobs=2) + monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"power": mean_power}, 100, n_jobs=2) self.assertIn("seed=", str(ctx.exception)) def test_a_lambda_metric_is_named_rather_than_left_to_pickle(self): with self.assertRaises(ValueError) as ctx: - sweep(build_chain(), "noise.snr_dB", SNR_DB, - {"power": lambda y: float(np.mean(np.abs(y)))}, 100, - seed=3, n_jobs=2) + monte_carlo(build_chain(), "noise.snr_dB", SNR_DB, + {"power": lambda y: float(np.mean(np.abs(y)))}, 100, + seed=3, n_jobs=2) message = str(ctx.exception) self.assertIn("metrics", message) self.assertIn("lambda", message) def test_the_serial_path_still_accepts_a_lambda(self): """n_jobs=1 has no boundary to cross, and must stay permissive.""" - out = sweep(build_chain(), "noise.snr_dB", [5.0], - {"power": lambda y: float(np.mean(np.abs(y)))}, 100, - seed=3) + out = monte_carlo(build_chain(), "noise.snr_dB", [5.0], + {"power": lambda y: float(np.mean(np.abs(y)))}, 100, + seed=3) self.assertEqual(out["power"].shape, (1,)) diff --git a/tests/test_sweep.py b/tests/test_sweep.py deleted file mode 100644 index 530c169..0000000 --- a/tests/test_sweep.py +++ /dev/null @@ -1,52 +0,0 @@ -"""sweep(): reconfigure, reseed, run, collect (decision D35).""" -import unittest - -import numpy as np - -from comnumpy import AWGN, Sequential, SymbolDemapper, \ - SymbolGenerator, SymbolMapper, compute_ser, get_alphabet, sweep - - -class TestSweep(unittest.TestCase): - - def build(self): - alphabet = get_alphabet("QAM", 4) - return Sequential([ - SymbolGenerator(4, name="tx"), - SymbolMapper(alphabet), - AWGN(snr_dB=0, name="noise"), - SymbolDemapper(alphabet), - ]) - - def test_results_aligned_and_monotonic(self): - values = np.arange(0, 14, 4) - out = sweep(self.build(), "noise.snr_dB", values, - {"ser": compute_ser}, 20_000, reference="tx", seed=1) - self.assertEqual(out["ser"].shape, values.shape) - self.assertTrue(np.all(np.diff(out["ser"]) < 0)) # SER falls with SNR - - def test_reproducible_across_calls(self): - values = [0, 6] - a = sweep(self.build(), "noise.snr_dB", values, - {"ser": compute_ser}, 5_000, reference="tx", seed=3) - b = sweep(self.build(), "noise.snr_dB", values, - {"ser": compute_ser}, 5_000, reference="tx", seed=3) - np.testing.assert_array_equal(a["ser"], b["ser"]) - - def test_multi_param_zip(self): - chain = self.build() - out = sweep(chain, ["noise.snr_dB", "noise.name"], - [(0, "noise"), (10, "noise")], - {"snr_set": lambda t, y: chain["noise"].snr_dB}, - 1_000, reference="tx", seed=2) - np.testing.assert_array_equal(out["snr_set"], [0, 10]) - - def test_metric_without_reference(self): - out = sweep(self.build(), "noise.snr_dB", [0, 20], - {"power": lambda y: float(np.mean(np.abs(y) ** 2))}, - 2_000, seed=4) - self.assertEqual(len(out["power"]), 2) - - -if __name__ == "__main__": - unittest.main() diff --git a/validation/README.md b/validation/README.md index 3bbf1fb..e1fc842 100644 --- a/validation/README.md +++ b/validation/README.md @@ -28,6 +28,7 @@ python validation/core_awgn_ser.py | `optical_dbp_roundtrip.py` | Noiseless DBP with matched steps inverts the fiber to machine precision; CD-only equalization leaves a nonlinear residual growing with launch power | Ip & Kahn, JLT 26(20), 2008 | | `optical_raman.py` | Distributed Raman: on-off gain against the undepleted closed form `exp(g P_p L_eff)` for co, counter and bidirectional pumping; photon-number conservation `P_s/nu_s + P_p/nu_p` in the lossless limit, which holds *under depletion* where no closed form does; the counter/co gap opening monotonically with the pump; and the BVP path reproducing the IVP path | Agrawal, *Nonlinear Fiber Optics*, ch. 8; Islam, IEEE JSTQE 8(3), 2002 | | `optical_wdm_opticommpy.py` | A published 11 x 32 GBd PDM-16QAM WDM transmission over 700 km, reproduced against a **second implementation** rather than a closed form: the comb power matches the printed 8.41 dBm exactly, the ASE-limited SNR matches its own closed form to 0.02 dB, a fixed 0.5 km split step is shown to be 0.8 dB short of converged, the MI, GMI and normalized GMI land on the published 4.00 / 4.00 / 1.00, and the full link lands 1.0 dB above the published 20.63 dB -- the sign every impairment the notebook models and this one does not would give | OptiCommPy, `test_WDM_transmission.ipynb`; Essiambre et al., JLT 28(4), 2010; Marcuse et al., JLT 15(9), 1997 | +| `optical_dbp_hager.py` | The whole split-step + DBP stack against the **published effective-SNR curves** of Häger & Pfister (arXiv:2010.14258), on their own setup: Gaussian stimulus, 500 steps per span forward, logarithmic steps for the one- and two-step back-propagators. The amplifier-noise reference spans 16.0 to 29.8 dB over -6 to +8 dBm, against the 16.5-29.5 dB of the paper's figure for the same link; more steps buy a higher peak at a higher launch power, in order, and the best back-propagator comes within 3.5 dB of the noise-only bound while moving the optimum 8 dB to the right. `--full` runs both links at StPS = 500 (~1 h) | Häger & Pfister, IEEE JSAC 39(1), 2021; the reference implementation at github.com/chaeger/LDBP | | `optical_gn_model.py` | The Gaussian Noise model against **a published measurement and an independent split-step solver**: the closed form reproduces the 15-channel, 5 x 100 km `a_NL = -23.5` dB of Serena & Bononi to 0.20 dB; this library's own SSFM agrees to 0.01 dB on a five-channel WDM link and to better than 1 dB single-channel across 1-20 spans (the drift being the coherence exponent, measured at 0.21); the cubic law comes out at n = 3.006; the 27/8 Manakov-vs-scalar polarization factor of Table I is measured at 5.29 dB against 5.28; and the modulation-format gap the GN model is blind to is measured (16QAM 1.4 dB, QPSK 2.2 dB) against the published 15-channel 1.6 / 2.8 dB. `--full` reproduces the whole published link by SSFM (-23.15 dB, ~15 min) | Serena & Bononi, JLT 33(7), 2015; Poggiolini, JLT 30(24), 2012; GNPy `NliSolver` | | `optical_raman_gnpy.py` | Counter-pumped Raman against **a second implementation**, on the one thing the analytic checks cannot see: the gain *shape*. 96 channels over 4.8 THz amplified by two counter-propagating pumps, against GNPy's own expected output and its measured SSMF profile. The table normalizes to 0.4396 /W/km, the textbook SSMF figure; the profile agrees to **0.04 dB at worst** and the tilt to 6.53 dB against 6.59; the coupling coefficients match GNPy's exactly on the gain side, element by element. The effective-area correction is what buys it: tilt error 0.72 -> 0.06 dB, worst channel 0.92 -> 0.04 dB. GNPy's own solver was re-run at perturbative orders 1-4 and `numerical` to rule it out as an explanation; all five agree to the last digit **on this case** -- which is a property of counter-pumping, not of the solver, and does not carry over (see the next row) | GNPy `RamanSolver`, BSD-3-Clause; D'Amico et al., JLT 40, 3499-3511 (2022) §III.D | | `optical_raman_gnpy_regimes.py` | The same second implementation, over **four regimes and the spontaneous emission**, to answer what the residual is made of. Counter-pumped, co-propagating (the pump falls 25 dB while feeding the signals), both pumps at 600 mW, and 202 channels over 10 THz. Under GNPy's own conservation convention all four agree to **0.011 dB**, against a reference whose own 20 m / 5 m convergence is 0.030 dB. Under this library's convention the residual is 0.054 / 0.140 / 0.204 / 0.320 dB -- *ordered by depletion*, which is asserted, because a wrong coefficient would not sort itself that way and the phonon term would. Which convention is right is settled elsewhere, by the closed form: `tests/optical/test_raman.py` matches the photon one to 1e-11 and rejects the power-conserving one by 6.5%, the ratio of the two saturation limits. Spontaneous ASE agrees to **0.002 dB** on the source GNPy models; the +0.17 dB left over is the inter-channel spontaneous emission GNPy omits, and the decomposition shows it accounts for all of it. Generating the reference required pinning `method='numerical'` -- GNPy's default second-order expansion is 0.4 dB off co-propagating -- and matching the fact that a direct `RamanSolver` call applies the connector to the pumps but not the signals | GNPy `RamanSolver`, BSD-3-Clause; Agrawal, *Nonlinear Fiber Optics*, §8.1 (Manley-Rowe) | diff --git a/validation/core_awgn_ser.py b/validation/core_awgn_ser.py index 2ee8ace..d6e5a14 100644 --- a/validation/core_awgn_ser.py +++ b/validation/core_awgn_ser.py @@ -11,7 +11,7 @@ from comnumpy import (AWGN, Sequential, SymbolDemapper, SymbolGenerator, SymbolMapper, compute_ser, Constellation) -from comnumpy.sweep import sweep +from comnumpy.monte_carlo import monte_carlo FIG_DIR = pathlib.Path(__file__).parent / "figures" @@ -27,9 +27,9 @@ def simulate(order, seed): AWGN(snr_dB=0, name="noise"), SymbolDemapper(constellation), ]) - results = sweep(chain, "noise.snr_dB", SNR_DB_RANGE, - {"ser": compute_ser}, N_SYMBOLS, - reference="tx", seed=seed) + results = monte_carlo(chain, "noise.snr_dB", SNR_DB_RANGE, + {"ser": compute_ser}, N_SYMBOLS, + reference="tx", seed=seed) return results["ser"] diff --git a/validation/fec_coding_gain.py b/validation/fec_coding_gain.py index 3ccdcc7..3279ffe 100644 --- a/validation/fec_coding_gain.py +++ b/validation/fec_coding_gain.py @@ -7,13 +7,20 @@ the shared build/sweep/collect skeleton (decision D35). """ import pathlib +from functools import partial import numpy as np from comnumpy import AWGN, Sequential, SymbolGenerator +from comnumpy.core.metrics import compute_ber from comnumpy.core.mappers import SymbolDemapper, SymbolMapper from comnumpy.fec import ConvolutionalEncoder, ViterbiDecoder -from comnumpy.sweep import sweep +from comnumpy.monte_carlo import monte_carlo + +# one bit per symbol: the chain carries bits, not symbol indices. Bound +# through partial rather than a lambda, which a worker process could not +# unpickle if this run were parallelised. +BER_BITS = partial(compute_ber, width=1) FIG_DIR = pathlib.Path(__file__).parent / "figures" @@ -25,7 +32,7 @@ BPSK = np.array([1.0 + 0.0j, -1.0 + 0.0j]) # bit 0 -> +1, bit 1 -> -1 -def coded_chain(soft): +def get_coded_chain(soft): bpsk = BPSK demapper = (SymbolDemapper(bpsk, soft=True, name="demap") if soft else SymbolDemapper(bpsk, name="demap")) @@ -39,7 +46,7 @@ def coded_chain(soft): ]) -def uncoded_chain(): +def get_uncoded_chain(): bpsk = BPSK return Sequential([ SymbolGenerator(2, name="tx"), @@ -54,15 +61,15 @@ def main(): snr_coded = EBN0_DB_RANGE + 10 * np.log10(CODE_RATE) ber = {} - ber["uncoded"] = sweep(uncoded_chain(), "noise.snr_dB", EBN0_DB_RANGE, - {"ber": compute_ber_bits}, N_BITS, - reference="tx", seed=10)["ber"] - ber["hard"] = sweep(coded_chain(soft=False), "noise.snr_dB", snr_coded, - {"ber": compute_ber_bits}, N_BITS, - reference="tx", seed=20)["ber"] - ber["soft"] = sweep(coded_chain(soft=True), "noise.snr_dB", snr_coded, - {"ber": compute_ber_bits}, N_BITS, - reference="tx", seed=30)["ber"] + ber["uncoded"] = monte_carlo(get_uncoded_chain(), "noise.snr_dB", EBN0_DB_RANGE, + {"ber": BER_BITS}, N_BITS, + reference="tx", seed=10)["ber"] + ber["hard"] = monte_carlo(get_coded_chain(soft=False), "noise.snr_dB", snr_coded, + {"ber": BER_BITS}, N_BITS, + reference="tx", seed=20)["ber"] + ber["soft"] = monte_carlo(get_coded_chain(soft=True), "noise.snr_dB", snr_coded, + {"ber": BER_BITS}, N_BITS, + reference="tx", seed=30)["ber"] # Hard decisions operate 10*log10(1/R) = 3 dB below the uncoded channel # SNR, so the hard curve crosses the uncoded one around 3-4 dB (Proakis @@ -96,10 +103,5 @@ def main(): f"soft {ber['soft'][-1]:.2e}") -def compute_ber_bits(bits_tx, bits_rx): - """Bit error rate between two bit streams (decoded length = tx length).""" - return float(np.mean(bits_tx != bits_rx)) - - if __name__ == "__main__": main() diff --git a/validation/ofdm_awgn_ser.py b/validation/ofdm_awgn_ser.py index 573bada..1716891 100644 --- a/validation/ofdm_awgn_ser.py +++ b/validation/ofdm_awgn_ser.py @@ -20,7 +20,7 @@ from comnumpy.ofdm.processors import (CarrierAllocator, CarrierExtractor, CyclicPrefixer, CyclicPrefixRemover, FFTProcessor, IFFTProcessor) -from comnumpy.sweep import sweep +from comnumpy.monte_carlo import monte_carlo FIG_DIR = pathlib.Path(__file__).parent / "figures" @@ -52,9 +52,9 @@ def ofdm_chain(): def main(): - results = sweep(ofdm_chain(), "noise.snr_dB", SNR_DB_RANGE, - {"ser": compute_ser_truncated}, N_SYMBOLS, - reference="tx", seed=7) + results = monte_carlo(ofdm_chain(), "noise.snr_dB", SNR_DB_RANGE, + {"ser": compute_ser}, N_SYMBOLS, + reference="tx", seed=7) ser_sim = results["ser"] # The AWGN block measures the *time-domain* power. With unit-power @@ -90,11 +90,5 @@ def main(): f"over {mask.sum()} SNR points") -def compute_ser_truncated(tx, rx): - """SER over the transmitted length (S/P zero-padding adds trailing symbols).""" - n = len(tx) - return compute_ser(tx, rx[:n]) - - if __name__ == "__main__": main() diff --git a/validation/optical_dbp_hager.py b/validation/optical_dbp_hager.py new file mode 100644 index 0000000..8a808d2 --- /dev/null +++ b/validation/optical_dbp_hager.py @@ -0,0 +1,239 @@ +"""Digital back-propagation against Häger and Pfister's published curves. + +Decision D7 asks that a claim be confronted with something outside this +library. The claim here is the whole split-step + DBP stack, and the +confrontation is the effective-SNR curves of + + C. Häger and H. D. Pfister, "Physics-based deep learning for fiber-optic + communication systems", IEEE JSAC 39(1), 2021 (arXiv:2010.14258), + +whose reference implementation is +https://github.com/chaeger/LDBP/blob/master/ldbp/ldbp.py. + +**The setup is theirs, not a teaching variant of it.** A *Gaussian* +stimulus, because the paper's curves are Gaussian-input effective SNR +and a real constellation would sit above them by the modulation-format +gap the GN model is blind to (measured in ``optical_gn_model.py``); 500 +split steps per span forward, so the channel itself is not in question; +and the two links of the paper -- 25 x 80 km at 10.7 GBd, and 10 x 100 km +at 32 GBd. The low-step back-propagators use **logarithmic** steps, which +is what the paper uses and what makes one or two steps per span worth +anything: a logarithmic grid puts the short steps where the power, and +therefore the nonlinearity, is highest. + +The tutorial ``examples/optical/NLI_simulation.py`` shares none of that +on purpose -- it carries 16-QAM so it can show a symbol error rate, and +fewer steps so it runs in minutes. It teaches; this reproduces. + +**What is asserted, and what is only printed.** The paper's numbers are +read off figures, so this script asserts the structure those figures +show rather than values transcribed by eye: + +* the linear reference climbs monotonically -- it is P/P_ASE, and + nothing bends it; +* every other receiver peaks, because the nonlinearity eventually grows + faster than the signal; +* more steps buy a higher peak at a higher launch power, in order; +* the best back-propagator stays within a few decibels of the + reference -- 3.48 dB on the default run, the rest being the + amplifier noise, which back-propagation amplifies rather than + removes. + +The table itself is printed in full so the curves can be laid over the +paper's figures, which is the check a reader actually wants. + +**One departure from the original script, and it matters.** That script +band-limited with ``BWFilter(1 / oversampling_sim)`` before decimating by +``oversampling_sim // oversampling_dsp`` -- a brick wall at half the +width the decimation needed, cutting the roll-off shoulders off the +signal it was meant to pass. It put a deterministic distortion floor +under every curve, which is why the reference saturated instead of +climbing. ``Downsampler(use_filter=True)`` applies the mask the +decimation actually calls for. + +The default run is the paper's second link at StPS = 60 with two trials +per point, and takes about ten seconds -- 512 Gaussian symbols per trial +is all an effective-SNR estimate needs. ``--full`` runs both links at the +paper's StPS = 500 with ten trials, which is of the order of an hour. + +Measured on the default run, against the axes of the paper's own figure +for this link (its y range is 16.5 to 29.5 dB over -6 to +8 dBm): + + launch [dBm] -6.0 -4.0 -2.0 0.0 2.0 4.0 6.0 8.0 + linear link 16.04 17.71 19.81 21.98 23.75 25.99 27.72 29.75 + linear 15.97 17.32 18.22 17.87 15.99 10.61 8.62 4.09 + DBP 1 15.96 17.50 18.83 19.06 17.35 12.70 10.40 5.52 + DBP 2 (log) 16.04 17.64 19.62 21.55 22.11 20.12 17.36 11.08 + DBP 4 (log) 16.03 17.67 19.77 21.87 23.58 24.75 24.68 19.68 + DBP full 16.03 17.68 19.78 21.89 23.68 25.47 26.27 25.71 +""" +import sys +import time + +import numpy as np + +from comnumpy.core import Sequential +from comnumpy.core.compensators import DataAidedComplexGainCompensator +from comnumpy.core.filters import SRRCFilter +from comnumpy.core.generators import GaussianGenerator +from comnumpy.core.metrics import compute_effective_snr +from comnumpy.core.processors import Amplifier, Downsampler, Upsampler +from comnumpy.optical.dbp import DBP +from comnumpy.optical.utils import dbm_to_watt, launch_amplitude +from comnumpy.optical.links import FiberLink + +FULL = "--full" in sys.argv + +# The two links of the paper. `dBm` is the launch-power axis of its figures. +LINKS = { + "25 x 80 km, 10.7 GBd": dict(baud_Hz=10.7e9, span_km=80.0, n_spans=25, + dBm=np.arange(-10.0, 6.0, 2.0)), + "10 x 100 km, 32 GBd": dict(baud_Hz=32e9, span_km=100.0, n_spans=10, + dBm=np.arange(-6.0, 9.0, 2.0)), +} + +# The receivers of the paper, in the order its legend lists them. The +# first is the reference: the same link with its Kerr term switched off, +# so it measures the amplifier noise alone. +RECEIVERS = [ + ("linear link", 1, True, "linear"), + ("linear", 1, True, "linear"), + ("DBP 1", 1, False, "linear"), + ("DBP 2 (log)", 2, False, "logarithmic"), + ("DBP 4 (log)", 4, False, "logarithmic"), + ("DBP full", 500 if FULL else 60, False, "linear"), +] + +N_SYMBOLS = 2 ** 9 +OVERSAMPLING_SIM = 6 +OVERSAMPLING_DSP = 2 +NF_DB = 5.0 +ROLLOFF = 0.1 +STPS = 500 if FULL else 60 +N_TRIALS = 10 if FULL else 2 + + +def get_channel(link, *, linear_only): + """Gaussian symbols to the field the receiver sees, launch included.""" + fs = link["baud_Hz"] * OVERSAMPLING_SIM + ratio = OVERSAMPLING_SIM // OVERSAMPLING_DSP + return Sequential([ + GaussianGenerator(1.0, name="tx"), + Amplifier(1.0, name="launch"), + Upsampler(OVERSAMPLING_SIM, scale=np.sqrt(OVERSAMPLING_SIM)), + SRRCFilter(ROLLOFF, OVERSAMPLING_SIM, method="fft"), + FiberLink(N_spans=link["n_spans"], L_span=link["span_km"], StPS=STPS, + NF_dB=NF_DB, fs=fs, use_only_linear=linear_only, + name="link"), + Downsampler(ratio, use_filter=True), + ], taps=["tx"]) + + +def get_receiver(link, steps, linear_only, step_type, reference): + """Back-propagation, matched filter, and the complex gain removed. + + The paper corrects amplitude and phase with a single complex scalar + fitted on the transmitted signal; that is exactly + ``DataAidedComplexGainCompensator(shared=True)``, so the estimate + comes from the library rather than from four lines of algebra. + """ + fs = link["baud_Hz"] * OVERSAMPLING_SIM + ratio = OVERSAMPLING_SIM // OVERSAMPLING_DSP + return Sequential([ + DBP(link["n_spans"], L_span=link["span_km"], StPS=steps, + step_type=step_type, fs=fs / ratio, use_only_linear=linear_only, + name="dbp"), + SRRCFilter(ROLLOFF, OVERSAMPLING_DSP, method="fft", + scale=1 / np.sqrt(OVERSAMPLING_DSP)), + Downsampler(OVERSAMPLING_DSP), + DataAidedComplexGainCompensator(reference, shared=True, name="gain"), + ]) + + +def run(name, link): + """Effective SNR against launch power, one column per receiver.""" + powers_dBm = link["dBm"] + snr_dB = np.zeros((len(powers_dBm), len(RECEIVERS))) + linear = get_channel(link, linear_only=True) + nonlinear = get_channel(link, linear_only=False) + + for index, dBm in enumerate(powers_dBm): + amplitude = launch_amplitude(dbm_to_watt(float(dBm))) + for trial in range(N_TRIALS): + fields = {} + for is_linear, channel in ((True, linear), (False, nonlinear)): + channel.seed(index * N_TRIALS + trial) + channel.set_params(launch__gain=amplitude) + fields[is_linear] = channel(N_SYMBOLS) + sent = linear.tap("tx") + + for column, (_, steps, only_linear, step_type) in enumerate(RECEIVERS): + receiver = get_receiver(link, steps, only_linear, step_type, + sent) + estimate = receiver(fields[column == 0]) + snr_dB[index, column] += ( + compute_effective_snr(sent, estimate) / N_TRIALS) + + return powers_dBm, 10 * np.log10(snr_dB) + + +def report(name, powers_dBm, snr_dB): + header = "launch [dBm] " + for value in powers_dBm: + header += f"{value:7.1f}" + print(f"\n{name} ({STPS} steps/span forward, {N_TRIALS} trials)") + print(header) + for column, (label, _, _, _) in enumerate(RECEIVERS): + line = f"{label:13s} " + for value in snr_dB[:, column]: + line += f"{value:7.2f}" + print(line) + + +def check(name, powers_dBm, snr_dB): + """The structure the paper's figures show, asserted.""" + reference = snr_dB[:, 0] + increments = np.diff(reference) + assert np.all(increments > 0), ( + f"{name}: the amplifier-noise reference must climb with the launch " + f"power -- it is P/P_ASE and nothing bends it; got {reference}") + + peaks, best_dBm = [], [] + for column in range(1, len(RECEIVERS)): + curve = snr_dB[:, column] + index = int(np.argmax(curve)) + peaks.append(curve[index]) + best_dBm.append(powers_dBm[index]) + assert index < len(curve) - 1, ( + f"{name}: {RECEIVERS[column][0]} has not turned over inside the " + f"swept range; the nonlinearity should eventually win") + + for earlier in range(len(peaks) - 1): + assert peaks[earlier + 1] >= peaks[earlier] - 0.05, ( + f"{name}: more steps must not lower the peak -- " + f"{RECEIVERS[earlier + 1][0]} {peaks[earlier]:.2f} dB then " + f"{RECEIVERS[earlier + 2][0]} {peaks[earlier + 1]:.2f} dB") + assert best_dBm[-1] >= best_dBm[0], ( + f"{name}: back-propagation must move the optimum to the right, " + f"got {best_dBm[0]:+.1f} dBm then {best_dBm[-1]:+.1f} dBm") + + gap = float(np.max(reference) - peaks[-1]) + print(f" best back-propagation is {gap:.2f} dB from the noise-only " + f"reference, and moved the optimum by " + f"{best_dBm[-1] - best_dBm[0]:+.1f} dB") + assert gap < 6.0, f"{name}: {gap:.2f} dB left on the table" + + +def main(): + links = LINKS if FULL else {"10 x 100 km, 32 GBd": LINKS["10 x 100 km, 32 GBd"]} + start = time.perf_counter() + for name, link in links.items(): + powers_dBm, snr_dB = run(name, link) + report(name, powers_dBm, snr_dB) + check(name, powers_dBm, snr_dB) + print(f"\ntotal {time.perf_counter() - start:.0f} s" + f"{'' if FULL else ' (run with --full for both links at StPS=500)'}") + + +if __name__ == "__main__": + main()