diff --git a/.claude/skills/comnumpy-tutorial/SKILL.md b/.claude/skills/comnumpy-tutorial/SKILL.md
index 0fe6378..1d3c763 100644
--- a/.claude/skills/comnumpy-tutorial/SKILL.md
+++ b/.claude/skills/comnumpy-tutorial/SKILL.md
@@ -54,7 +54,20 @@ 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: `monte_carlo(chain, param, values, metrics, stimulus, seed=)`.
+- Monte-Carlo: `monte_carlo(chain, param, values, metrics, stimulus, seed=)`
+  when one chain and standard metrics cover the study. When the study
+  is *around* a chain -- several detectors on one frame, several
+  receivers on one propagation -- write the loop out, on **one storage
+  convention**: declare the methods and the metrics first; pre-allocate
+  one `np.zeros(n_points)` per (metric, method), indexed by **name** on
+  both levels, never by position -- a 2-D `np.zeros((n_points,
+  n_methods))` puts a countable offset between parallel columns, which
+  is the classic silent bug; draw one child seed per point from a
+  master seed (`np.random.SeedSequence(seed).spawn`, D6/D35); let the
+  loop fill the arrays; display last, from the same dictionaries the
+  loop filled -- each inner dict is exactly the `curves` that
+  `print_data` and `plot_data` render. The loop is the pedagogy: it
+  stays visible.
 - **A swept result is one dictionary, shown two ways.**
   `data = {"x": snr_dB, "curves": {"ZF": ..., "ML": ...}}` — which is
   already the shape `monte_carlo` returns — then `print_data(data,
@@ -66,12 +79,17 @@ then learns the wrong lesson -- they came to see how the library is used.
 - Activate the style sheet once, at the top of the script, right after
   the imports: `style.use()`. The colours and the figure size are
   rcParams and a figure already created keeps the old ones.
-- Figures: **draw with matplotlib, decorate with `style.apply(ax, kind)`**.
-  A scatter of real against imaginary, a `semilogy` of a rate against an
-  SNR: the reader knows those calls, and hiding them behind a wrapper
-  teaches nothing. `style.apply(ax, "iq" | "error_rate" | "time" |
-  "spectrum")` fills the labels that are still empty, turns on the grid
-  and adds the legend; it never touches the data or the scales.
+- Figures: **the function when one exists, `style.apply(ax, kind)` when
+  none does.** `plot_iq` for a constellation, `plot_spectrum`,
+  `plot_error_rate` for measured curves paired with their closed forms:
+  a page that shows the same kind of figure three or seven times gains
+  nothing from spelling out `ax.plot(np.real(x), np.imag(x), ".")` each
+  time, and loses the equal aspect ratio the day someone forgets it.
+  When the figure is a one-off with no function behind it -- a runtime,
+  an effective SNR -- draw it with matplotlib and hand the axis to
+  `style.apply(ax, "iq" | "error_rate" | "time" | "spectrum")`, which
+  fills the labels still empty, turns on the grid and adds the legend;
+  it never touches the data or the scales.
   A kind names the *quantity*. A runtime in ms or an effective SNR in dB
   is none of them: give it its own labels and `ax.grid(True)`, and do not
   borrow `"error_rate"` because the axis happens to be logarithmic.
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9045779..84f544f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -45,6 +45,66 @@ 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`) |
 
+### Changed — the simulation sections read identically across the pages
+
+Every sweep in the tutorials now shows the same three markers, in the
+same order:
+
+```python
+# --- metrics, pre-allocated ---
+# --- simulation loop ---
+# --- results: tables and figures ---
+```
+
+The pass also removed what review found on the way: a list
+comprehension in `one_shot_ofdm.py` that the tutorial rules ban and
+that hid a `monte_carlo` call inside an expression; the last
+`append`-and-convert accumulations (`one_shot_ofdm`, `one_shot_mimo`,
+`one_shot_alamouti`, the by-hand loop of `monte_carlo_awgn`), replaced
+by arrays pre-allocated to zeros and filled by index; a hand-rolled
+aligned `print` loop in `one_shot_mimo` that `print_data` replaces; a
+dead five-line dict-to-array conversion in `one_shot_NLI`; and a
+`bits per symbol` column that is now `np.log2(orders).astype(int)`
+rather than a loop -- NumPy over iteration wherever the loop taught
+nothing.
+
+### Changed — one storage convention for every simulation loop
+
+Every study in the examples was the same sentence -- run the same
+simulation for several values of one parameter, keep what each run
+measured -- and every script spelled its storage out differently:
+`np.zeros((n_points, n_methods))` indexed by position, an accumulator
+per metric divided by `N_test` at the end, and, in the two MIMO
+detector sweeps, **no seed at all** -- the curves in the documentation
+were reproducible by nobody, their author included. One of those
+positional tables carried a `+1` column offset, the textbook silent
+bug.
+
+The four sweeps -- both MIMO detector comparisons, the
+chromatic-dispersion compensator study and the launch-power sweep of
+the DBP tutorial -- now follow one convention, written into the
+tutorial skill:
+
+* the methods and the metrics are **declared first**, as the ordered
+  dictionaries they are;
+* storage is **pre-allocated to zeros, one array per (metric, method),
+  indexed by name** on both levels -- a column never has to be counted
+  to be found, and a misplaced `+1` between parallel tables is not
+  expressible;
+* the simulation loop draws **one child seed per point** from a master
+  seed (`np.random.SeedSequence(seed).spawn`, decisions D6/D35) and
+  fills the arrays;
+* the display comes last, from the same dictionaries the loop filled --
+  each inner dictionary is exactly the `curves` that `print_data` and
+  `plot_data` render.
+
+An `Experiment` object that ran the loop behind a callback was written,
+measured and removed the same day, before any release: it made the
+scripts shorter nowhere, and it hid exactly the loop the tutorials are
+supposed to teach. What the tutorials needed was a storage convention,
+not an engine. The reproducibility it briefly carried stays: the seeds
+above are now in the scripts themselves.
+
 ### Changed — the figures finally use the style sheet that ships with the package
 
 `comnumpy.mplstyle` has shipped since D27b -- Okabe-Ito, colourblind-safe,
diff --git a/docs/getting_started/first_simulation.rst b/docs/getting_started/first_simulation.rst
index f889e51..3e63e2b 100644
--- a/docs/getting_started/first_simulation.rst
+++ b/docs/getting_started/first_simulation.rst
@@ -170,7 +170,7 @@ It also draws itself, next to what the channel did to it:
 
 .. literalinclude:: ../../examples/simple/one_shot_awgn.py
    :language: python
-   :lines: 89-101
+   :lines: 89-97
 
 .. image:: img/first_simulation_fig1.png
    :width: 100%
diff --git a/docs/tutorials/alamouti.rst b/docs/tutorials/alamouti.rst
index 891f4e7..db028b9 100644
--- a/docs/tutorials/alamouti.rst
+++ b/docs/tutorials/alamouti.rst
@@ -199,7 +199,7 @@ One-shot simulation
 
 .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py
    :language: python
-   :lines: 94-112
+   :lines: 94-107
 
 .. code::
 
@@ -241,7 +241,7 @@ inverts:
 
 .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py
    :language: python
-   :lines: 117-167
+   :lines: 112-169
 
 The three curves are three readings of one closed form,
 :func:`~comnumpy.core.metrics.compute_ser_rayleigh_psk`: :math:`L` branches
@@ -260,7 +260,7 @@ the simulation is what confronts them:
 
 .. literalinclude:: ../../examples/mimo/one_shot_alamouti.py
    :language: python
-   :lines: 169-208
+   :lines: 171-210
 
 .. code::
 
diff --git a/docs/tutorials/awgn.rst b/docs/tutorials/awgn.rst
index 434a5f0..a9bad99 100644
--- a/docs/tutorials/awgn.rst
+++ b/docs/tutorials/awgn.rst
@@ -110,7 +110,7 @@ outline marks a tapped block:
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 141-147
+   :lines: 144-150
 
 Monte Carlo Simulation
 """"""""""""""""""""""
@@ -122,7 +122,7 @@ against the transmitted symbols.
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 44-51
+   :lines: 44-53
 
 Three chain services appear there, and they are the ones every study is made
 of.
@@ -141,7 +141,7 @@ at every point:
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 53-64
+   :lines: 55-66
 
 .. code::
 
@@ -209,7 +209,7 @@ two expressions above transcribe line for line:
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 67-80
+   :lines: 69-82
 
 From the constellation
 """"""""""""""""""""""
@@ -222,7 +222,7 @@ by :math:`k` above happens inside.
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 82-87
+   :lines: 84-89
 
 .. code::
 
@@ -244,7 +244,7 @@ share a colour, so a pair reads as one statement.
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 89-95
+   :lines: 91-97
 
 .. image:: img/monte_carlo_awgn.png
    :width: 100%
@@ -272,7 +272,7 @@ not guess (decision D41), so it happens at the call site:
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 97-127
+   :lines: 99-128
 
 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
@@ -300,7 +300,7 @@ disagreeing with theory.
 
 .. literalinclude:: ../../examples/simple/monte_carlo_awgn.py
    :language: python
-   :lines: 129-139
+   :lines: 130-142
 
 .. image:: img/monte_carlo_awgn_orders.png
    :width: 100%
diff --git a/docs/tutorials/img/monte_carlo_mimo_fig1.png b/docs/tutorials/img/monte_carlo_mimo_fig1.png
index c6a8f76..a717f4b 100644
Binary files a/docs/tutorials/img/monte_carlo_mimo_fig1.png and b/docs/tutorials/img/monte_carlo_mimo_fig1.png differ
diff --git a/docs/tutorials/img/monte_carlo_mimo_fig2.png b/docs/tutorials/img/monte_carlo_mimo_fig2.png
index bc87503..a297d5c 100644
Binary files a/docs/tutorials/img/monte_carlo_mimo_fig2.png and b/docs/tutorials/img/monte_carlo_mimo_fig2.png differ
diff --git a/docs/tutorials/img/monte_carlo_mimo_fig4.png b/docs/tutorials/img/monte_carlo_mimo_fig4.png
index 6edaa9c..9f41e21 100644
Binary files a/docs/tutorials/img/monte_carlo_mimo_fig4.png and b/docs/tutorials/img/monte_carlo_mimo_fig4.png differ
diff --git a/docs/tutorials/img/nli_simulation_fig1.png b/docs/tutorials/img/nli_simulation_fig1.png
index e304c18..02d2c10 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 d105909..c9a1fa1 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_alamouti_fig2.png b/docs/tutorials/img/one_shot_alamouti_fig2.png
index 98d4073..15a0ded 100644
Binary files a/docs/tutorials/img/one_shot_alamouti_fig2.png and b/docs/tutorials/img/one_shot_alamouti_fig2.png differ
diff --git a/docs/tutorials/img/one_shot_nli_fig2.png b/docs/tutorials/img/one_shot_nli_fig2.png
index 55c7398..f832be6 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_fig4.png b/docs/tutorials/img/one_shot_nli_fig4.png
index af85e0e..c7556e7 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/img/one_shot_ofdm_fig2.png b/docs/tutorials/img/one_shot_ofdm_fig2.png
index 3bbb73f..3d08807 100644
Binary files a/docs/tutorials/img/one_shot_ofdm_fig2.png and b/docs/tutorials/img/one_shot_ofdm_fig2.png differ
diff --git a/docs/tutorials/img/one_shot_ofdm_fig3.png b/docs/tutorials/img/one_shot_ofdm_fig3.png
index de092c3..eb15bd0 100644
Binary files a/docs/tutorials/img/one_shot_ofdm_fig3.png and b/docs/tutorials/img/one_shot_ofdm_fig3.png differ
diff --git a/docs/tutorials/img/one_shot_ofdm_fig5.png b/docs/tutorials/img/one_shot_ofdm_fig5.png
index f94d76c..7d44490 100644
Binary files a/docs/tutorials/img/one_shot_ofdm_fig5.png and b/docs/tutorials/img/one_shot_ofdm_fig5.png differ
diff --git a/docs/tutorials/mermaid/gn_model.mmd b/docs/tutorials/mermaid/gn_model.mmd
index bb0e88d..e5f15b2 100644
--- a/docs/tutorials/mermaid/gn_model.mmd
+++ b/docs/tutorials/mermaid/gn_model.mmd
@@ -1,18 +1,18 @@
 flowchart LR
     source["SymbolGenerator"]
     tx["SymbolMapper"]
+    launch["Amplifier"]
     upsampler["Upsampler"]
     srrcfilter["SRRCFilter"]
-    launch["Amplifier"]
     fibre["FiberLink"]
     dbp["DBP"]
     srrcfilter_2["SRRCFilter"]
     downsampler["Downsampler"]
     source --> tx
-    tx --> upsampler
+    tx --> launch
+    launch --> upsampler
     upsampler --> srrcfilter
-    srrcfilter --> launch
-    launch --> fibre
+    srrcfilter --> fibre
     fibre --> dbp
     dbp --> srrcfilter_2
     srrcfilter_2 --> downsampler
diff --git a/docs/tutorials/mimo.rst b/docs/tutorials/mimo.rst
index a5b158c..06e6cbe 100644
--- a/docs/tutorials/mimo.rst
+++ b/docs/tutorials/mimo.rst
@@ -69,7 +69,7 @@ block:
 
 .. literalinclude:: ../../examples/mimo/one_shot_mimo.py
    :language: python
-   :lines: 183-186
+   :lines: 182-185
 
 Each chain is given the same seed before running, so the five numbers below
 differ by the detector alone -- same symbols, same noise, same channel:
@@ -89,7 +89,7 @@ Let us look at what each receive antenna sees, read from the ``"noise"`` tap:
 
 .. literalinclude:: ../../examples/mimo/one_shot_mimo.py
    :language: python
-   :lines: 65-73
+   :lines: 65-72
 
 .. image:: img/monte_carlo_mimo_fig1.png
    :width: 100%
@@ -103,7 +103,7 @@ Applying the pseudo-inverse of the channel matrix separates the streams:
 
 .. literalinclude:: ../../examples/mimo/one_shot_mimo.py
    :language: python
-   :lines: 75-87
+   :lines: 74-81
 
 .. image:: img/monte_carlo_mimo_fig2.png
    :width: 100%
@@ -233,15 +233,20 @@ the matrix the signal goes through **and** the one the detector inverts:
 
 .. literalinclude:: ../../examples/mimo/one_shot_mimo.py
    :language: python
-   :lines: 89-133
+   :lines: 83-129
 
 .. code::
 
-   ZF    0.3313 0.2106 0.1065 0.0428 0.0148 0.0047 0.0020
-   MMSE  0.2871 0.1762 0.0876 0.0345 0.0116 0.0036 0.0014
-   OSIC  0.2826 0.1610 0.0668 0.0192 0.0037 0.0007 0.0003
-   ML    0.2807 0.1560 0.0620 0.0159 0.0027 0.0006 0.0003
-   SD    0.2807 0.1560 0.0620 0.0159 0.0027 0.0006 0.0003
+   SER
+   snr_dB       ZF     MMSE       OSIC         ML         SD
+   ---------------------------------------------------------
+        0  0.33129  0.28710  2.826e-01  2.807e-01  2.807e-01
+        3  0.21064  0.17623  1.610e-01  1.560e-01  1.560e-01
+        6  0.10650  0.08759  6.679e-02  6.196e-02  6.196e-02
+        9  0.04285  0.03451  1.920e-02  1.588e-02  1.588e-02
+       12  0.01485  0.01159  3.738e-03  2.737e-03  2.737e-03
+       15  0.00474  0.00363  7.125e-04  5.500e-04  5.500e-04
+       18  0.00203  0.00142  3.375e-04  2.625e-04  2.625e-04
 
 .. image:: img/monte_carlo_mimo_fig3.png
    :width: 100%
@@ -276,16 +281,16 @@ where that product no longer fits. We therefore time both detectors on
 
 .. literalinclude:: ../../examples/mimo/one_shot_mimo.py
    :language: python
-   :lines: 135-181
+   :lines: 131-180
 
 .. code::
 
-     ML  0 dB    142.3 ms    65536 nodes   SER 0.7900
-     ML  9 dB    129.7 ms    65536 nodes   SER 0.5044
-     ML 18 dB    127.2 ms    65536 nodes   SER 0.0112
-     SD  0 dB    236.3 ms     91.8 nodes   SER 0.7900
-     SD  9 dB     49.5 ms     19.2 nodes   SER 0.5044
-     SD 18 dB     16.7 ms      5.8 nodes   SER 0.0112
+     ML  0 dB    160.3 ms    65536 nodes   SER 0.7900
+     ML  9 dB    120.5 ms    65536 nodes   SER 0.5044
+     ML 18 dB    109.9 ms    65536 nodes   SER 0.0112
+     SD  0 dB    203.0 ms     91.8 nodes   SER 0.7900
+     SD  9 dB     42.0 ms     19.2 nodes   SER 0.5044
+     SD 18 dB     14.3 ms      5.8 nodes   SER 0.0112
 
 .. image:: img/monte_carlo_mimo_fig4.png
    :width: 100%
diff --git a/docs/tutorials/ofdm.rst b/docs/tutorials/ofdm.rst
index 272a0d0..a499bed 100644
--- a/docs/tutorials/ofdm.rst
+++ b/docs/tutorials/ofdm.rst
@@ -136,7 +136,7 @@ Implementation
 
 .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py
    :language: python
-   :lines: 55-86
+   :lines: 55-77
 
 Results
 """""""
@@ -187,7 +187,7 @@ before the channel and a receiver block after it:
 
 .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py
    :language: python
-   :lines: 87-115
+   :lines: 78-101
 
 ``OFDMTransmitter`` and ``OFDMReceiver`` are themselves chains:
 serial-to-parallel conversion, subcarrier allocation, IFFT and cyclic prefix
@@ -218,21 +218,21 @@ range of SNR values:
 
 .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py
    :language: python
-   :lines: 116-138
+   :lines: 102-132
 
 .. code::
 
    SER
    SNR [dB]  single carrier    OFDM
    --------------------------------
-          6          0.6273  0.5234
-          8          0.5270  0.3953
-         10          0.4297  0.2973
-         12          0.2707  0.1980
-         14          0.1512  0.1426
-         16          0.0656  0.0672
-         18          0.0152  0.0383
-         20          0.0020  0.0195
+          6         0.62734  0.5234
+          8         0.52695  0.3953
+         10         0.42969  0.2973
+         12         0.27070  0.1980
+         14         0.15117  0.1426
+         16         0.06562  0.0672
+         18         0.01523  0.0383
+         20         0.00195  0.0195
 
 .. image:: img/one_shot_ofdm_fig4.png
    :width: 100%
@@ -265,15 +265,17 @@ block length grows:
 
 .. literalinclude:: ../../examples/ofdm/one_shot_ofdm.py
    :language: python
-   :lines: 139-174
+   :lines: 133-169
 
 .. code::
 
-        N   single carrier      OFDM     ratio
-      128           8.5 ms    1.01 ms        8
-      256          28.4 ms    0.82 ms       35
-      512         139.1 ms    0.84 ms      166
-     1024         901.0 ms    0.88 ms     1018
+   receiver runtime [ms], and their ratio
+   block length N  single carrier   OFDM   ratio
+   ---------------------------------------------
+              128            8.74  0.613    14.2
+              256           32.64  0.644    50.7
+              512          209.35  0.999   209.5
+             1024         1685.07  1.123  1500.4
 
 .. image:: img/one_shot_ofdm_fig5.png
    :width: 100%
diff --git a/docs/tutorials/optical_fiber_nonlinearity.rst b/docs/tutorials/optical_fiber_nonlinearity.rst
index 6d02690..e460cab 100644
--- a/docs/tutorials/optical_fiber_nonlinearity.rst
+++ b/docs/tutorials/optical_fiber_nonlinearity.rst
@@ -135,9 +135,9 @@ be read against.
 
    spans   measured   ASE only   the fibre      SER     phase     time
        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
+       5    25.93 dB   30.65 dB     4.72 dB   0.0000    -7.0 deg    0.3 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
+      15    18.56 dB   25.94 dB     7.38 dB   0.0023   -22.9 deg    1.1 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
 
@@ -190,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: 158-165
+   :lines: 158-160
 
 .. image:: img/one_shot_nli_fig1.png
    :width: 100%
@@ -204,7 +204,7 @@ three distances:
 
 .. literalinclude:: ../../examples/optical/one_shot_NLI.py
    :language: python
-   :lines: 167-197
+   :lines: 162-187
 
 .. image:: img/one_shot_nli_fig2.png
    :width: 100%
@@ -236,7 +236,7 @@ through, which says where it all went.
 
 .. literalinclude:: ../../examples/optical/one_shot_NLI.py
    :language: python
-   :lines: 199-204
+   :lines: 189-194
 
 .. code::
 
@@ -289,7 +289,7 @@ Results
 
 .. literalinclude:: ../../examples/optical/one_shot_NLI.py
    :language: python
-   :lines: 206-232
+   :lines: 196-218
 
 .. code::
 
@@ -363,21 +363,30 @@ instead of two, which is what ``polarizations=1`` is for:
 
 That is the prediction the sweep below has to land on.
 
+The sweep follows the storage convention every study of this library
+uses. The receivers and the metrics are declared first; each (metric,
+receiver) pair gets a pre-allocated array of zeros, indexed by **name**
+on both levels -- a column never has to be counted to be found, and a
+misplaced ``+1`` between parallel tables is not expressible. The
+simulation loop fills those arrays, one child seed per launch power so
+the master seed alone reproduces the study, and everything printed or
+drawn afterwards reads from the same dictionaries the loop filled.
+
 .. literalinclude:: ../../examples/optical/NLI_simulation.py
    :language: python
-   :lines: 112-227
+   :lines: 112-241
 
 .. code::
 
    effective SNR [dB]
    launch power [dBm]       -6.0  -4.5  -3.0  -1.5   0.0   1.5   3.0   4.5
    -----------------------------------------------------------------------
-   amplifier noise only     15.9  17.4  18.8  20.3  21.9  23.3  24.8  26.4
-   dispersion compensation  15.8  17.1  18.2  18.8  18.4  16.9  14.6  11.7
-   DBP, 1 step/span         15.8  17.2  18.4  19.3  19.3  18.1  16.1  13.3
-   DBP, 2 steps/span        15.9  17.3  18.7  19.9  20.6  20.5  19.2  16.9
-   DBP, 4 steps/span        15.9  17.4  18.8  20.3  21.7  23.0  23.8  23.8
-   DBP, 50 steps/span       15.9  17.4  18.8  20.3  21.8  23.2  24.6  25.9
+   amplifier noise only     15.9  17.4  18.9  20.4  21.9  23.4  24.9  26.4
+   dispersion compensation  15.8  17.1  18.3  18.8  18.3  16.9  14.4  11.8
+   DBP, 1 step/span         15.8  17.2  18.5  19.3  19.2  18.2  16.0  13.4
+   DBP, 2 steps/span        15.8  17.3  18.7  19.9  20.7  20.6  19.1  16.9
+   DBP, 4 steps/span        15.9  17.3  18.9  20.3  21.7  23.0  23.8  23.8
+   DBP, 50 steps/span       15.9  17.4  18.9  20.4  21.8  23.3  24.6  25.9
 
 .. image:: img/nli_simulation_fig1.png
    :width: 100%
@@ -395,26 +404,28 @@ and it climbs as :math:`P/P_{\mathrm{ASE}}` forever -- 15.9 dB at
 defect in it.
 
 The maximum of each of the others is the operating point of that receiver,
-and back-propagation **moves it to the right**:
+and back-propagation **moves it to the right** -- hesitantly at one step per
+span, whose curve is nearly flat between :math:`-1.5` and :math:`0` dBm, then
+decisively:
 
 .. code::
 
    receiver                  best SNR   at power    total time
-   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
+   amplifier noise only      26.37 dB    4.5 dBm       0.4 s
+   dispersion compensation   18.85 dB   -1.5 dBm       0.3 s
+   DBP, 1 step/span          19.29 dB   -1.5 dBm       0.7 s
+   DBP, 2 steps/span         20.65 dB    0.0 dBm       1.3 s
+   DBP, 4 steps/span         23.82 dB    3.0 dBm       2.4 s
+   DBP, 50 steps/span        25.91 dB    4.5 dBm      28.4 s
 
 Start with the second row, because it is the one the closed form claims to
 predict:
 
 .. code::
 
-   GN model 18.77 dB at -1.35 dBm, dispersion compensation 18.78 dB at -1.5 dBm
+   GN model 18.77 dB at -1.35 dBm, dispersion compensation 18.85 dB at -1.5 dBm
 
-**0.01 dB.** The Gaussian-noise model of the previous tutorial ran in
+**0.08 dB.** The Gaussian-noise model of the previous tutorial ran in
 microseconds and knew nothing of this chain -- not the pulse shaping, not the
 split-step integration, not the matched filter. It described the receiver
 that leaves the nonlinearity in place, which is exactly the receiver its
@@ -431,7 +442,7 @@ Read the last two columns of those rows together. Going from one step per
 span to four buys **4.5 dB**; going from four to fifty buys **2.1 dB more**
 and costs twelve times the computation. The returns collapse because the
 step-size error falls with the number of steps while the ASE noise does not:
-at fifty steps the receiver is 0.41 dB from the noise-only bound, and no
+at fifty steps the receiver is 0.46 dB from the noise-only bound, and no
 number of steps will close that last gap.
 
 Against dispersion compensation alone, four steps per span are worth 5.0 dB
@@ -455,8 +466,8 @@ say so in errors rather than in a caveat:
 
    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
+   dispersion compensation        6 over 49152 symbols
+   DBP, 1 step/span               4 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
diff --git a/examples/mimo/monte_carlo_simulation_1.py b/examples/mimo/monte_carlo_simulation_1.py
index a80c382..c318607 100644
--- a/examples/mimo/monte_carlo_simulation_1.py
+++ b/examples/mimo/monte_carlo_simulation_1.py
@@ -1,6 +1,17 @@
+"""Block error rate of three MIMO detectors over Rayleigh fading.
+
+Reproduces figure 2 of [1]: a 2x2 QPSK link, maximum-likelihood
+detection against zero-forcing and ordered successive interference
+cancellation, one new channel draw per frame.
+
+[1] X. Li, H. C. Huang, A. Lozano and G. J. Foschini,
+"Reduced-complexity detection algorithms for systems using
+multi-element arrays", Globecom '00, San Francisco, 2000.
+"""
 import numpy as np
 import matplotlib.pyplot as plt
 
+from comnumpy import plot_data, print_data, style
 from comnumpy.core import Sequential
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper
@@ -8,83 +19,73 @@
 from comnumpy.core.utils import Constellation
 from comnumpy.mimo.channels import AWGN, FlatMIMOChannel
 from comnumpy.mimo.utils import rayleigh_channel
-from comnumpy.mimo.detectors import MaximumLikelihoodDetector, LinearDetector, OrderedSuccessiveInterferenceCancellationDetector
-from comnumpy import style
+from comnumpy.mimo.detectors import (
+    LinearDetector, MaximumLikelihoodDetector,
+    OrderedSuccessiveInterferenceCancellationDetector)
 
 style.use()
 
-# This script reproduces the figure 2 of the article [1]
-#
-# X. Li, H. C. Huang, A. Lozano and G. J. Foschini, "Reduced-complexity detection algorithms for systems using multi-element arrays,"
-# Globecom '00 - IEEE. Global Telecommunications Conference. Conference Record (Cat. No.00CH37137), San Francisco, CA, USA, 2000, pp. 1072-1076 vol.2,
-
-
-# Parameters
-N_test = 500  # increase the value for smoothing the ber
-N = 400
+# parameters
+N = 400          # symbols per frame
+N_test = 500     # channel draws per point; more smooths the curve
 N_r, N_t = 2, 2
 constellation = Constellation("PSK", 4)
-M = constellation.order
-
-H = rayleigh_channel(N_r=N_r, N_t=N_t)
-
-# construct chain
-chain = Sequential([SymbolGenerator(constellation.order, name="data_tx"),
-                    SymbolMapper(constellation),
-                    FlatMIMOChannel(H, name="channel"),
-                    AWGN(sigma2=0, name="noise")
-                    ], taps=["data_tx"])
-
-# prepare MC trial
-detector_names = ["ML", "ZF", "OSIC"]
-
-# compute snr list from snr per bit in dB
 snr_dB_list = np.arange(0, 45, 5)
-
-# perform simulation
-bler_data = np.zeros((len(snr_dB_list), len(detector_names)))
-sig_power = N_t
-
-for index_snr, snr_dB in enumerate(snr_dB_list):
-    sigma2 = sig_power*(10**(-snr_dB/10))
+seed = 42        # one master seed reproduces the whole figure
+detector_names = ("ML", "ZF", "OSIC")
+
+chain = Sequential([
+    SymbolGenerator(constellation.order, name="data_tx"),
+    SymbolMapper(constellation),
+    FlatMIMOChannel(rayleigh_channel(N_r=N_r, N_t=N_t), name="channel"),
+    AWGN(sigma2=1.0, name="noise"),
+    ], taps=["data_tx"])
+
+# --- metrics, pre-allocated ------------------------------------------
+# One array per detector, full of zeros, indexed by name -- a column
+# never has to be counted to be found, and the dictionary is what the
+# table and the figure render.
+bler = {}
+for name in detector_names:
+    bler[name] = np.zeros(len(snr_dB_list))
+
+# --- simulation loop -------------------------------------------------
+# One child seed per SNR point (D6/D35): the whole figure is reproduced
+# by the master seed alone.
+point_seeds = np.random.SeedSequence(seed).spawn(len(snr_dB_list))
+for index, snr_dB in enumerate(snr_dB_list):
+    rng = np.random.default_rng(int(point_seeds[index].generate_state(1)[0]))
+    sigma2 = N_t * 10 ** (-snr_dB / 10)
     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)
+        # new channel realization, decided by every detector: the
+        # comparison is over identical frames
+        H = rayleigh_channel(N_r=N_r, N_t=N_t, rng=rng)
+        chain.seed(int(rng.integers(2 ** 31)))
         chain.set_params(channel__H=H)
-
-        # generate data
         Y = chain((N_t, N))
         S_ref = chain.tap("data_tx")
 
-        # test detector
-        for index, detector_name in enumerate(detector_names):
-
-            match detector_name:
-                case "ML":
-                    detector = MaximumLikelihoodDetector(alphabet=constellation, H=H)
-                case "ZF":
-                    detector = LinearDetector(alphabet=constellation, H=H, method="zf")
-                case "OSIC":
-                    detector = OrderedSuccessiveInterferenceCancellationDetector(constellation, osic_type="sinr", H=H, sigma2=sigma2, name="OSIC")
-
-            # perform detection
-            S_est = detector(Y)
-            # evaluate metrics
-            bler_data[index_snr, index] += (compute_ber(S_ref, S_est, width=int(np.log2(M))) > 0)
-
-    bler_data[index_snr, :] /= N_test
-
-# plot figures
-for index, detector_name in enumerate(detector_names):
-    plt.semilogy(snr_dB_list, bler_data[:, index], label=detector_name)
-plt.ylabel("BLER")
-plt.xlabel("SNR (dB)")
-plt.xlim([0, 40])
-plt.ylim([10**-3, 1])
-plt.legend()
-plt.grid(True)
-plt.title("Performance Comparison of ZF and ML detector, 2*2 QPSK")
+        detectors = {
+            "ML": MaximumLikelihoodDetector(alphabet=constellation, H=H),
+            "ZF": LinearDetector(alphabet=constellation, H=H, method="zf"),
+            "OSIC": OrderedSuccessiveInterferenceCancellationDetector(
+                constellation, osic_type="sinr", H=H, sigma2=sigma2),
+        }
+        for name, detector in detectors.items():
+            ber = compute_ber(S_ref, detector(Y),
+                              width=constellation.bits_per_symbol)
+            # a frame is in error when at least one of its bits is
+            bler[name][index] += float(ber > 0) / N_test
+
+# --- results: table and figure ---------------------------------------
+data = {"x": snr_dB_list, "curves": bler}
+print_data(data, xlabel="snr_dB", ylabel="BLER")
+
+ax = plot_data(data, xlabel="SNR [dB]", ylabel="BLER", yscale="log",
+               marker="o", fillstyle="none")
+ax.set_xlim(0, 40)
+ax.set_ylim(1e-3, 1)
+ax.set_title("2x2 QPSK over Rayleigh fading, one draw per frame")
 plt.show()
diff --git a/examples/mimo/monte_carlo_simulation_2.py b/examples/mimo/monte_carlo_simulation_2.py
index a16a6a9..115cadd 100644
--- a/examples/mimo/monte_carlo_simulation_2.py
+++ b/examples/mimo/monte_carlo_simulation_2.py
@@ -1,6 +1,16 @@
+"""BER of three OSIC orderings over Rayleigh fading, 4x4 16-QAM.
+
+Reproduces figure 11.3 of [1]: ordered successive interference
+cancellation, the ordering criterion being the column norm, the
+post-detection SNR or the post-detection SINR.
+
+[1] Cho, Yong Soo, et al. MIMO-OFDM wireless communications with
+MATLAB. John Wiley & Sons, 2010.
+"""
 import numpy as np
 import matplotlib.pyplot as plt
 
+from comnumpy import plot_data, print_data, style
 from comnumpy.core import Sequential
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper
@@ -8,83 +18,79 @@
 from comnumpy.core.utils import Constellation
 from comnumpy.mimo.channels import AWGN, FlatMIMOChannel
 from comnumpy.mimo.utils import rayleigh_channel
-from comnumpy.mimo.detectors import OrderedSuccessiveInterferenceCancellationDetector
-from comnumpy import style
+from comnumpy.mimo.detectors import (
+    OrderedSuccessiveInterferenceCancellationDetector)
 
 style.use()
 
-# This script shows the performance of several OSIC MIMO detectors. It reproduces the figure 11.3 of the book [1]
-# [1] Cho, Yong Soo, et al. MIMO-OFDM wireless communications with MATLAB. John Wiley & Sons, 2010.
-
-# Parameters
-N_test = 5000  # increase the value for smoothing the ber
-N = 100
+# parameters
+N = 100          # symbols per frame
+N_test = 5000    # channel draws per point; more smooths the curve
 N_r, N_t = 4, 4
 constellation = Constellation("QAM", 16)
-M = constellation.order
-
-H = rayleigh_channel(N_r=N_r, N_t=N_t)
-
-# construct chain
-chain = Sequential([SymbolGenerator(constellation.order, name="data_tx"),
-                    SymbolMapper(constellation),
-                    FlatMIMOChannel(H, name="channel"),
-                    AWGN(sigma2=0, name="noise"),
-                    ], taps=["data_tx"])
-
-# prepare MC trial
-detector_names = ["colnorm", "snr", "sinr"]
-
-# compute snr list from snr per bit in dB
 snr_dB_list = np.arange(0, 45, 5)
-
-# 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))
+seed = 42        # one master seed reproduces the whole figure
+orderings = ("colnorm", "snr", "sinr")
+
+chain = Sequential([
+    SymbolGenerator(constellation.order, name="data_tx"),
+    SymbolMapper(constellation),
+    FlatMIMOChannel(rayleigh_channel(N_r=N_r, N_t=N_t), name="channel"),
+    AWGN(sigma2=1.0, name="noise"),
+    ], taps=["data_tx"])
+
+# --- metrics, pre-allocated ------------------------------------------
+# One array per ordering, indexed by name. The error counts get one
+# too -- they say where the estimate runs out of samples, so they are
+# part of the record even though they stay out of the figure.
+ber = {}
+errors = {}
+for ordering in orderings:
+    ber[ordering] = np.zeros(len(snr_dB_list))
+    errors[ordering] = np.zeros(len(snr_dB_list))
+
+# --- simulation loop -------------------------------------------------
+# One child seed per SNR point (D6/D35): the whole figure is reproduced
+# by the master seed alone.
+point_seeds = np.random.SeedSequence(seed).spawn(len(snr_dB_list))
+for index, snr_dB in enumerate(snr_dB_list):
+    rng = np.random.default_rng(int(point_seeds[index].generate_state(1)[0]))
+    sigma2 = N_t * 10 ** (-snr_dB / 10)
     chain.set_params(noise__sigma2=sigma2)
 
-    for _ in range(N_test):
+    # one ErrorCounter per ordering: each BER is a ratio of error and
+    # bit totals over the whole point, not a mean of per-frame rates
+    counters = {}
+    for ordering in orderings:
+        counters[ordering] = ErrorCounter(
+            width=constellation.bits_per_symbol)
 
-        # new channel realization
-        H = rayleigh_channel(N_r=N_r, N_t=N_t)
+    for _ in range(N_test):
+        H = rayleigh_channel(N_r=N_r, N_t=N_t, rng=rng)
+        chain.seed(int(rng.integers(2 ** 31)))
         chain.set_params(channel__H=H)
-
-        # generate data
         Y = chain((N_t, N))
         S_ref = chain.tap("data_tx")
 
-        # test detector
-        for index, detector_name in enumerate(detector_names):
-            # create detector
-            detector = OrderedSuccessiveInterferenceCancellationDetector(alphabet=constellation, osic_type=detector_name, H=H, sigma2=sigma2)
-            # perform detection
-            S_est = detector(Y)
-            # evaluate metrics
-            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
-
-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):
-    plt.semilogy(snr_dB_list, ber_data[:, index], label=detector_name)
-plt.ylabel("BER")
-plt.xlabel("SNR (dB)")
-plt.xlim([0, 40])
-plt.ylim([10**-4, 1])
-plt.legend()
-plt.grid(True)
-plt.title("Performance of OSIC methods (NT=4, Nr=4, 16QAM)")
+        for ordering, counter in counters.items():
+            detector = OrderedSuccessiveInterferenceCancellationDetector(
+                alphabet=constellation, osic_type=ordering, H=H,
+                sigma2=sigma2)
+            counter.update(S_ref, detector(Y))
+
+    for ordering, counter in counters.items():
+        ber[ordering][index] = counter.rate
+        errors[ordering][index] = counter.n_errors
+
+# --- results: tables and figure --------------------------------------
+print_data({"x": snr_dB_list, "curves": ber}, xlabel="snr_dB", ylabel="BER")
+print()
+print_data({"x": snr_dB_list, "curves": errors},
+           xlabel="snr_dB", ylabel="errors behind each point")
+
+ax = plot_data({"x": snr_dB_list, "curves": ber}, xlabel="SNR [dB]",
+               ylabel="BER", yscale="log", marker="o", fillstyle="none")
+ax.set_xlim(0, 40)
+ax.set_ylim(1e-4, 1)
+ax.set_title("OSIC orderings, 4x4 16-QAM over Rayleigh fading")
 plt.show()
diff --git a/examples/mimo/one_shot_CMA.py b/examples/mimo/one_shot_CMA.py
index e83b6cf..d1f9858 100644
--- a/examples/mimo/one_shot_CMA.py
+++ b/examples/mimo/one_shot_CMA.py
@@ -3,7 +3,7 @@
 from dataclasses import dataclass
 
 from comnumpy import style
-from comnumpy.core import Sequential
+from comnumpy.core import Sequential, plot_iq
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper
 from comnumpy.core.filters import SRRCFilter
@@ -90,11 +90,7 @@ def process_phase_correction(self, Y_sub):
 Y = chain((N_t, N))
 
 data = chain.tap("filter")  # full CMA output (before tail removal)
-symbols = chain.tap("tail")
-_, ax = plt.subplots()
-ax.plot(np.real(symbols), np.imag(symbols), ".")
-ax.set_title("after CMA convergence")
-style.apply(ax, "iq")
+plot_iq(chain.tap("tail"), title="after CMA convergence")
 
 # compute CMA loss
 kernel_size = 100
diff --git a/examples/mimo/one_shot_alamouti.py b/examples/mimo/one_shot_alamouti.py
index 8745794..43af67f 100644
--- a/examples/mimo/one_shot_alamouti.py
+++ b/examples/mimo/one_shot_alamouti.py
@@ -13,7 +13,7 @@
 from comnumpy.core.metrics import compute_ser, compute_ser_rayleigh_psk
 from comnumpy.core.processors import Amplifier
 from comnumpy.core.utils import Constellation
-from comnumpy.core.visualizers import plot_error_rate
+from comnumpy.core.visualizers import plot_error_rate, plot_iq
 from comnumpy.mimo.channels import AWGN, FlatMIMOChannel
 from comnumpy.mimo.coding import SpaceTimeDecoder, SpaceTimeEncoder, get_code
 from comnumpy.mimo.detectors import LinearDetector
@@ -99,15 +99,10 @@ def get_chain(kind, H, sigma2=sigma2):
 received = alamouti.tap("noise")[0]
 combined = alamouti.tap("detector") / power
 fig2, (ax_left, ax_right) = plt.subplots(nrows=1, ncols=2, figsize=(9, 4.2))
-ax_left.plot(np.real(received), np.imag(received), ".")
+plot_iq(received, marker=".", ax=ax_left)
 ax_left.set_title("tap('noise'): the single receive antenna")
-style.apply(ax_left, "iq")
-ax_right.plot(np.real(combined), np.imag(combined), ".", label="combined")
-ax_right.plot(np.real(constellation.alphabet),
-              np.imag(constellation.alphabet), "kx", markersize=9,
-              label="transmitted alphabet")
+plot_iq(combined, reference=constellation, ax=ax_right)
 ax_right.set_title("tap('detector'): after Alamouti combining")
-style.apply(ax_right, "iq")
 plt.tight_layout()
 plt.savefig(f"{img_dir}/one_shot_alamouti_fig2.png")
 
@@ -144,12 +139,19 @@ def average_ser(chain, n_rx, n_tx, snr_dB, stimulus, n_channels, seed=0):
     "Alamouti, 2 Tx, 1 Rx": (alamouti, 1, 2, n_symbols),
     "MRC, 1 Tx, 2 Rx": (mrc, 2, 1, (1, n_symbols)),
 }
+# --- metrics, pre-allocated ------------------------------------------
 curves = {}
+for name in links:
+    curves[name] = np.zeros(len(snr_dB_list))
+
+# --- simulation loop -------------------------------------------------
 for name, (chain, n_rx, n_tx, stimulus) in links.items():
-    values = []
-    for snr_dB, count in zip(snr_dB_list, draws, strict=True):
-        values.append(average_ser(chain, n_rx, n_tx, snr_dB, stimulus, count))
-    curves[name] = values
+    for index, (snr_dB, count) in enumerate(zip(snr_dB_list, draws,
+                                                strict=True)):
+        curves[name][index] = average_ser(chain, n_rx, n_tx, snr_dB,
+                                          stimulus, count)
+
+# --- results: figure against the closed forms ------------------------
 
 fine = np.linspace(snr_dB_list[0], snr_dB_list[-1], 100)
 per_bit = 10 ** (fine / 10) / np.log2(M)
diff --git a/examples/mimo/one_shot_mimo.py b/examples/mimo/one_shot_mimo.py
index 884832b..daf0dec 100644
--- a/examples/mimo/one_shot_mimo.py
+++ b/examples/mimo/one_shot_mimo.py
@@ -3,13 +3,13 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from comnumpy import monte_carlo, style
+from comnumpy import monte_carlo, print_data, style
 from comnumpy.core import Sequential
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper
 from comnumpy.core.metrics import compute_ser
 from comnumpy.core.utils import Constellation
-from comnumpy.core.visualizers import plot_error_rate
+from comnumpy.core.visualizers import plot_error_rate, plot_iq
 from comnumpy.mimo.channels import AWGN, FlatMIMOChannel
 from comnumpy.mimo.detectors import (
     LinearDetector, MaximumLikelihoodDetector,
@@ -65,25 +65,19 @@ def get_chain(detector):
 Y = chains["ZF"].tap("noise")
 fig1, axes1 = plt.subplots(nrows=1, ncols=N_r, figsize=(4 * N_r, 4))
 for index in range(N_r):
-    axes1[index].plot(np.real(Y[index, :]), np.imag(Y[index, :]), ".")
+    plot_iq(Y[index, :], ax=axes1[index])
     axes1[index].set_title(f"Received signal (antenna {index + 1})")
     axes1[index].set_xlim([-2, 2])
     axes1[index].set_ylim([-2, 2])
-    style.apply(axes1[index], "iq")
 plt.savefig(f"{img_dir}/monte_carlo_mimo_fig1.png")
 
 Z = detectors["ZF"].linear_estimator(Y)
 fig2, axes2 = plt.subplots(nrows=1, ncols=N_t, figsize=(4 * N_t, 4))
 for index in range(N_t):
-    axes2[index].plot(np.real(Z[index, :]), np.imag(Z[index, :]), ".",
-                      label="estimated")
-    axes2[index].plot(np.real(constellation.alphabet),
-                      np.imag(constellation.alphabet), "kx",
-                      markersize=9, label="transmitted alphabet")
+    plot_iq(Z[index, :], reference=constellation, ax=axes2[index])
     axes2[index].set_title(f"Estimated signal (stream {index + 1})")
     axes2[index].set_xlim([-2, 2])
     axes2[index].set_ylim([-2, 2])
-    style.apply(axes2[index], "iq")
 plt.savefig(f"{img_dir}/monte_carlo_mimo_fig2.png")
 
 snr_dB_list = np.arange(0, 20, 3)
@@ -112,17 +106,19 @@ def average_ser(name, chain, snr_dB, seed=0):
     return float(np.mean(results["ser"]))
 
 
+# --- metrics, pre-allocated ------------------------------------------
 curves = {}
+for name in chains:
+    curves[name] = np.zeros(len(snr_dB_list))
+
+# --- simulation loop -------------------------------------------------
 for name, chain in chains.items():
-    values = []
-    for snr_dB in snr_dB_list:
-        values.append(average_ser(name, chain, snr_dB))
-    curves[name] = values
-for name, values in curves.items():
-    line = f"{name:5s} "
-    for value in values:
-        line += f"{value:.4f} "
-    print(line)
+    for index, snr_dB in enumerate(snr_dB_list):
+        curves[name][index] = average_ser(name, chain, snr_dB)
+
+# --- results: table and figure ---------------------------------------
+print_data({"x": snr_dB_list, "curves": curves},
+           xlabel="snr_dB", ylabel="SER")
 
 ax = plot_error_rate(snr_dB_list, curves, ylabel="SER",
                      title=f"{N_r}x{N_t} MIMO, {constellation.order}-"
@@ -143,9 +139,12 @@ def average_ser(name, chain, snr_dB, seed=0):
     "ML": MaximumLikelihoodDetector(big_constellation, H=big_H, name="detector"),
     "SD": big_decoder,
 }
+# --- metrics, pre-allocated ------------------------------------------
 elapsed = {}
 for name in big_detectors:
-    elapsed[name] = []
+    elapsed[name] = np.zeros(len(snr_dB_list))
+
+# --- simulation loop -------------------------------------------------
 for name, detector in big_detectors.items():
     big_chain = Sequential([
         SymbolGenerator(big_constellation.order, name="tx"),
@@ -154,16 +153,16 @@ def average_ser(name, chain, snr_dB, seed=0):
         AWGN(sigma2=1.0, name="noise"),
         detector,
     ], taps=["tx"], name=f"4x4 MIMO, {name}")
-    for snr_dB in snr_dB_list:
+    for index, snr_dB in enumerate(snr_dB_list):
         big_chain.seed(4)
         big_chain.set_params(noise__sigma2=4 * 10 ** (-snr_dB / 10))
         start = time.perf_counter()
         detected = big_chain((4, 400))
-        elapsed[name].append((time.perf_counter() - start) * 1e3)
+        elapsed[name][index] = (time.perf_counter() - start) * 1e3
         errors = compute_ser(big_chain.tap("tx"), detected)
         nodes = (f"{big_decoder.nodes_:7.1f} nodes" if name == "SD"
                  else f"{16 ** 4:7d} nodes")
-        print(f"  {name:2s} {snr_dB:2d} dB {elapsed[name][-1]:8.1f} ms  "
+        print(f"  {name:2s} {snr_dB:2d} dB {elapsed[name][index]:8.1f} ms  "
               f"{nodes}   SER {errors:.4f}")
 
 # A runtime is not an error rate, so it does not go through
diff --git a/examples/ofdm/one_shot_ofdm.py b/examples/ofdm/one_shot_ofdm.py
index 364bca9..8ee4c4c 100644
--- a/examples/ofdm/one_shot_ofdm.py
+++ b/examples/ofdm/one_shot_ofdm.py
@@ -16,7 +16,7 @@
 from comnumpy.core.mappers import SymbolDemapper, SymbolMapper
 from comnumpy.core.metrics import compute_ser
 from comnumpy.core.utils import Constellation
-from comnumpy.core.visualizers import plot_error_rate
+from comnumpy.core.visualizers import plot_error_rate, plot_iq
 from comnumpy.ofdm.chains import OFDMReceiver, OFDMTransmitter
 
 style.use()
@@ -68,19 +68,10 @@
 sc_ser = compute_ser(sc_chain.tap("data_tx"), detected)
 print(f"single carrier: SER {sc_ser:.4f}, {sc_chain.elapsed_ * 1e3:.0f} ms")
 
-# An IQ plane is a scatter of the real part against the imaginary one;
-# style.apply gives it the axis labels, the grid and the equal aspect
-# ratio without which a constellation is not the constellation.
-alphabet = constellation.alphabet
 fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4.2))
 for ax, tap, name in [(axes[0], "data_rx", "received"),
                       (axes[1], "data_rx_eq", "after ZF equalization")]:
-    symbols = sc_chain.tap(tap)
-    ax.plot(np.real(symbols), np.imag(symbols), ".", label="symbols")
-    ax.plot(np.real(alphabet), np.imag(alphabet), "kx", markersize=9,
-            label="alphabet")
-    ax.set_title(name)
-    style.apply(ax, "iq")
+    plot_iq(sc_chain.tap(tap), reference=constellation, title=name, ax=ax)
 plt.tight_layout()
 plt.savefig(f"{img_dir}/one_shot_ofdm_fig2.png")
 
@@ -104,13 +95,8 @@
 print(f"OFDM          : SER {ofdm_ser:.4f}, {ofdm_chain.elapsed_ * 1e3:.2f} ms "
       f"({speedup:.0f} times faster)")
 
-equalized = ofdm_chain.tap("data_rx_eq")
-fig3, ax3 = plt.subplots()
-ax3.plot(np.real(equalized), np.imag(equalized), ".", label="symbols")
-ax3.plot(np.real(alphabet), np.imag(alphabet), "kx", markersize=9,
-         label="alphabet")
-ax3.set_title("OFDM, after one-tap equalization")
-style.apply(ax3, "iq")
+plot_iq(ofdm_chain.tap("data_rx_eq"), reference=constellation,
+        title="OFDM, after one-tap equalization")
 plt.savefig(f"{img_dir}/one_shot_ofdm_fig3.png")
 
 # --- error rate ------------------------------------------------------
@@ -118,16 +104,24 @@
 # independent noise seeds because 1280 symbols only resolve down to
 # 8e-4; the channel is the same throughout, by construction.
 snr_list = np.arange(6, 22, 2)
+n_trials = 2
+
+# --- metrics, pre-allocated ------------------------------------------
 measured = {}
+for name in ("single carrier", "OFDM"):
+    measured[name] = np.zeros(len(snr_list))
+
+# --- simulation loop -------------------------------------------------
 for name, chain in (("single carrier", sc_chain), ("OFDM", ofdm_chain)):
-    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)]
+    runs = np.zeros((n_trials, len(snr_list)))
+    for trial in range(n_trials):
+        runs[trial] = monte_carlo(chain, "data_rx.snr_dB", snr_list,
+                                  {"ser": compute_ser}, N,
+                                  reference="data_tx", seed=trial + 1)["ser"]
     measured[name] = np.mean(runs, axis=0)
 
-# A swept result is an abscissa and one series per curve, which is what
-# monte_carlo already returns. Written down once, it is printed and
-# plotted from the same object rather than restated for each.
+# --- results: table and figure ---------------------------------------
+
 ser_data = {"x": snr_list, "curves": measured}
 print()
 print_data(ser_data, xlabel="SNR [dB]", ylabel="SER")
@@ -139,35 +133,36 @@
 # --- what it costs ---------------------------------------------------
 # A chain records the wall time of its last pass in `elapsed_`, so the
 # run is also the measurement.
-lengths = [128, 256, 512, 1024]
-runtime = {"single carrier": [], "OFDM": []}
-for length in lengths:
+lengths = np.array([128, 256, 512, 1024])
+
+# --- metrics, pre-allocated ------------------------------------------
+runtime = {}
+for name in ("single carrier", "OFDM"):
+    runtime[name] = np.zeros(len(lengths))
+
+# --- simulation loop -------------------------------------------------
+for index, length in enumerate(lengths):
     for name, chain in (("single carrier", sc_chain), ("OFDM", ofdm_chain)):
         chain.seed(1)
-        chain(length)
-        runtime[name].append(1e3 * chain.elapsed_)
-
-runtime_data = {
-    "x": lengths,
-    "curves": {
-        "single carrier": np.array(runtime["single carrier"]),
-        "OFDM": np.array(runtime["OFDM"]),
-        "ratio": (np.array(runtime["single carrier"])
-                  / np.array(runtime["OFDM"])),
-    },
-}
+        chain(int(length))
+        runtime[name][index] = 1e3 * chain.elapsed_
+
+# --- results: table and figure ---------------------------------------
+
 print()
-print_data(runtime_data, xlabel="block length N",
+print_data({"x": lengths,
+            "curves": {"single carrier": runtime["single carrier"],
+                       "OFDM": runtime["OFDM"],
+                       "ratio": (runtime["single carrier"]
+                                 / runtime["OFDM"])}},
+           xlabel="block length N",
            ylabel="receiver runtime [ms], and their ratio")
 
-# The same dictionary, drawn. The ratio is dimensionless and does not
+# The same arrays, drawn. The ratio is dimensionless and does not
 # belong on an axis in milliseconds, so the figure takes the two
-# runtimes; everything else comes from the object the table printed.
-timed = {"x": lengths, "curves": {}}
-for name, values in runtime_data["curves"].items():
-    if name != "ratio":
-        timed["curves"][name] = values
-ax5 = plot_data(timed, xlabel="block length $N$",
+# runtimes alone.
+ax5 = plot_data({"x": lengths, "curves": runtime},
+                xlabel="block length $N$",
                 ylabel="receiver runtime [ms]",
                 xscale="log", yscale="log", marker="o", fillstyle="none")
 ax5.set_title("what equalization costs")
diff --git a/examples/optical/CD_compensation_part1.py b/examples/optical/CD_compensation_part1.py
index da21e27..9b2b0b9 100644
--- a/examples/optical/CD_compensation_part1.py
+++ b/examples/optical/CD_compensation_part1.py
@@ -1,6 +1,14 @@
+"""BER of two chromatic-dispersion FIR compensators, three modulations.
+
+Reproduces figures 3 and 4 of Eghbali et al., "Optimal Least-Squares
+FIR Digital Filters for Compensation of Chromatic Dispersion in Digital
+Coherent Optical Receivers": the truncated closed form of Savory
+against the least-squares design, on system 1 of their Table 1.
+"""
 import numpy as np
 import matplotlib.pyplot as plt
 
+from comnumpy import print_data, style
 from comnumpy.core import Sequential
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper, SymbolDemapper
@@ -9,122 +17,103 @@
 from comnumpy.core.channels import AWGN
 from comnumpy.core.metrics import compute_ber
 from comnumpy.core.filters import SRRCFilter
-from comnumpy.optical.compensators import ChromaticDispersionFIRCompensator, ChromaticDispersionLSFIRCompensator
+from comnumpy.optical.compensators import (
+    ChromaticDispersionFIRCompensator, ChromaticDispersionLSFIRCompensator)
 from comnumpy.optical.channels import ChromaticDispersion
-from comnumpy import style
 
 style.use()
 
-# This script reproduces the figures 3 and 4 in the paper Optimal Least-Squares FIR Digital Filters
-# for Compensation of Chromatic Dispersion in Digital Coherent Optical Receivers
-# see: https://www.researchgate.net/profile/Amir-Eghbali-2/publication/260642085_Optimal_Least-Squares_FIR_Digital_Filters_for_Compensation_of_ChromaticDispersion_in_Digital_Coherent_Optical_Receivers/links/02e7e534adbd8b0d77000000/Optimal-Least-Squares-FIR-Digital-Filters-for-Compensation-of-Chromatic-Dispersion-in-Digital-Coherent-Optical-Receivers.pdf
-
-system = 1
-type = "QAM"
-N = 200000  # increase number for smoothing MC curves
-N_h = 1000  # number of delay samples for the SRRC filter
-oversampling = 2  # oversampling factor
-rolloff = 0.25  # rolloff factor
-k_vect = [4, 6, 8]  # bits per symbol; the chain is built with the first one
-
-# system selection (see Table 1)
-if system == 1:
-    z = 4000  # km
-    fs = 21.4*(10**9)
-
-if system == 2:
-    z = 1000  # km
-    fs = 80*(10**9)
+# parameters
+N = 200000             # symbols per point; more smooths the curves
+z = 4000               # km -- system 1 of Table 1
+fs = 21.4e9
+N_h = 1000             # SRRC filter delay, in symbols
+oversampling = 2
+rolloff = 0.25
+k_vect = [4, 6, 8]     # bits per symbol; the chain is reconfigured per k
+snr_bit_dB_list = np.arange(10, 31)
+seed = 1               # one master seed reproduces the whole figure
+curve_names = ("theory", "Savory", "LS")
 
 comp_ref = ChromaticDispersionFIRCompensator(z, fs=fs)
-K = comp_ref.K  # compute length of filter
-N_filter = len(comp_ref.h)  # compute length of filter
-print("K: {}".format(K))
-print("filter length: {}".format(N_filter))
-total_delay = int(2*N_h*oversampling + np.floor(N_filter/2))  # total delay of the filters (need to be compensated before sampling)
-
-plt.figure()
+N_filter = len(comp_ref.h)
+print(f"K: {comp_ref.K}")
+print(f"filter length: {N_filter}")
+# total delay of the filters, compensated before sampling
+total_delay = int(2 * N_h * oversampling + np.floor(N_filter / 2))
 
-# create your chain and compensator list
 chain = Sequential([
-            SymbolGenerator(M=4, name="generator"),
-            SymbolMapper(Constellation(type, 2 ** k_vect[0]), name="mapper"),
-            Upsampler(oversampling),
-            SRRCFilter(rolloff, oversampling, N_h=N_h),
-            ChromaticDispersion(z, fs=fs),
-            AWGN(sigma2=0, name="noise")
-            ], taps=["generator"])
-
-full_compensator1 = Sequential([
-            ChromaticDispersionFIRCompensator(z, fs=fs),
-            SRRCFilter(rolloff, oversampling, N_h=N_h),
-            Downsampler(oversampling, phase=total_delay),
-            SymbolDemapper(Constellation(type, 2 ** k_vect[0]), name="demapper")
-            ])
-
-full_compensator2 = Sequential([
-            ChromaticDispersionLSFIRCompensator(z, N_filter, fs=fs, w_vect=[-np.pi, np.pi]),
-            SRRCFilter(rolloff, oversampling, N_h=N_h),
-            Downsampler(oversampling, phase=total_delay),
-            SymbolDemapper(Constellation(type, 2 ** k_vect[0]), name="demapper")
-            ])
-
-compensator_list = [full_compensator1, full_compensator2]
-
-# Monte carlo simulation
-ber_list_names = ["theoretical", "Savory", "LS"]
-SNR_vect = range(10, 31, 1)
-
+    SymbolGenerator(M=4, name="generator"),
+    SymbolMapper(Constellation("QAM", 2 ** k_vect[0]), name="mapper"),
+    Upsampler(oversampling),
+    SRRCFilter(rolloff, oversampling, N_h=N_h),
+    ChromaticDispersion(z, fs=fs),
+    AWGN(sigma2=1.0, name="noise"),
+    ], taps=["generator"])
+
+compensators = {
+    "Savory": Sequential([
+        ChromaticDispersionFIRCompensator(z, fs=fs),
+        SRRCFilter(rolloff, oversampling, N_h=N_h),
+        Downsampler(oversampling, phase=total_delay),
+        SymbolDemapper(Constellation("QAM", 2 ** k_vect[0]), name="demapper"),
+        ]),
+    "LS": Sequential([
+        ChromaticDispersionLSFIRCompensator(z, N_filter,
+                                            fs=fs,
+                                            w_vect=[-np.pi, np.pi]),
+        SRRCFilter(rolloff, oversampling, N_h=N_h),
+        Downsampler(oversampling, phase=total_delay),
+        SymbolDemapper(Constellation("QAM", 2 ** k_vect[0]), name="demapper"),
+        ]),
+}
+
+
+fig, ax = plt.subplots()
+line_style = {"theory": "k-", "Savory": "--", "LS": ":"}
 for k in k_vect:
-    constellation = Constellation(type, 2 ** k)
-    print("Modulation: {}{}".format(type, constellation.order))
-    epsilon_b = constellation.energy / k
-
-    # update chain and compensator parameters
-    # 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
+    constellation = Constellation("QAM", 2 ** k)
+    # Reconfigured through set_params, never by assignment (D50):
+    # SymbolMapper coerces its alphabet in __post_init__, and a direct
+    # write skips that.
     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 = []
-    ber_exp = []
-    ber_list = np.zeros((len(SNR_vect), len(ber_list_names)))
-
-    for index_SNR, SNR_bitdB in enumerate(SNR_vect):
-        snr_per_bit = 10 ** (SNR_bitdB/10)
-
-        # compute theoretical ber
-        ber_list[index_SNR, 0] = constellation.metrics(SNR_bitdB)["ber"]
-
-        # perform MC simulations
-        N0 = epsilon_b/snr_per_bit  # snr_bit = epsilon_b/N0 = (epsilon_s/log2(order))/N0
-        chain.set_params(noise__sigma2=N0)
+    for compensator in compensators.values():
+        compensator.set_params(demapper__alphabet=constellation)
+
+    # --- metrics, pre-allocated --------------------------------------
+    # One array per curve, indexed by name.
+    ber = {}
+    for name in curve_names:
+        ber[name] = np.zeros(len(snr_bit_dB_list))
+
+    # --- simulation loop ---------------------------------------------
+    # One child seed per SNR point (D6/D35).
+    point_seeds = np.random.SeedSequence(seed).spawn(len(snr_bit_dB_list))
+    for index, snr_bit_dB in enumerate(snr_bit_dB_list):
+        ber["theory"][index] = constellation.metrics(snr_bit_dB)["ber"]
+
+        sigma2 = constellation.energy / k / 10 ** (snr_bit_dB / 10)
+        chain.set_params(noise__sigma2=sigma2)
+        chain.seed(int(point_seeds[index].generate_state(1)[0]))
         y = chain(N)
-
-        # evaluate metric
         s = chain.tap("generator")
 
-        for index_comp, full_compensator in enumerate(compensator_list):
-            s_est = full_compensator(y)
-            ber_exp_temp = compute_ber(s_est, s[:len(s_est)], width=k)
-            ber_list[index_SNR, index_comp+1] = ber_exp_temp
-
-    # plot curves
-    color_list = ["k-", "r-", "b-"]
-    for index in range(3):
-        name = "{} ({}{})".format(ber_list_names[index], type,
-                                  constellation.order)
-        plt.semilogy(SNR_vect, ber_list[:, index], color_list[index], label=name)
-
-
-plt.ylabel("BER")
-plt.xlabel("SNR per bit (dB)")
-plt.ylim([10**-6, 0.2])
-plt.xlim([SNR_vect[0], SNR_vect[-1]])
-plt.legend()
-plt.grid()
+        for name, compensator in compensators.items():
+            s_est = compensator(y)
+            ber[name][index] = compute_ber(s_est, s[:len(s_est)], width=k)
+
+    # --- results: table and figure ---------------------------------------
+    print_data({"x": snr_bit_dB_list, "curves": ber},
+               xlabel="snr_bit_dB", ylabel=f"BER, QAM{2 ** k}")
+    print()
+    for name in curve_names:
+        ax.semilogy(snr_bit_dB_list, ber[name], line_style[name],
+                    label=f"{name} (QAM{2 ** k})")
+
+ax.set_xlabel("SNR per bit (dB)")
+ax.set_ylabel("BER")
+ax.set_ylim(1e-6, 0.2)
+ax.set_xlim(10, 30)
+ax.legend()
 plt.show()
diff --git a/examples/optical/NLI_simulation.py b/examples/optical/NLI_simulation.py
index ce43b32..1da34c6 100644
--- a/examples/optical/NLI_simulation.py
+++ b/examples/optical/NLI_simulation.py
@@ -7,7 +7,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from comnumpy import print_data
+from comnumpy import plot_data, print_data
 from comnumpy.core import Sequential
 from comnumpy.core.compensators import DataAidedPhaseCompensator
 from comnumpy.core.filters import SRRCFilter
@@ -61,9 +61,9 @@ def get_channel():
     return Sequential([
         SymbolGenerator(constellation.order, name="data_tx"),
         SymbolMapper(constellation, name="signal_tx"),
+        Amplifier(1.0, name="launch"),
         Upsampler(oversampling_sim, scale=np.sqrt(oversampling_sim)),
         SRRCFilter(rolloff, oversampling_sim, method="fft"),
-        Amplifier(1.0, name="launch"),
         FiberLink(N_spans=N_span, L_span=L_span, StPS=StPS, NF_dB=NF_dB,
                   fs=fs, fiber=fiber, name="link"),
         Downsampler(oversampling_ratio, use_filter=True),
@@ -122,55 +122,73 @@ def get_receiver(steps, linear_only, *, gain, reference):
 }
 
 channel = get_channel()
+
+# --- metrics, pre-allocated ------------------------------------------
+# One array per (metric, receiver), indexed by name on both levels -- a
+# column never has to be counted to be found, and each inner dictionary
+# is exactly what print_data and plot_data render.
 snr = {}
-counters = {}
-elapsed = {}
+ser = {}
+errors = {}
+times = {}
 for name in receivers:
     snr[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
-
+    ser[name] = np.zeros(len(dBm_list))
+    errors[name] = np.zeros(len(dBm_list))
+    times[name] = np.zeros(len(dBm_list))
+
+# --- simulation loop -------------------------------------------------
+# One child seed per launch power (D6/D35): the whole study is
+# reproduced by the master seed alone.
+seed = 0
+point_seeds = np.random.SeedSequence(seed).spawn(len(dBm_list))
 for index, dBm in enumerate(dBm_list):
+    point_seed = int(point_seeds[index].generate_state(1)[0])
     amp = launch_amplitude(dbm_to_watt(dBm))
+
+    # one counter per receiver: the rate is a ratio of error and symbol
+    # totals over the trials, not a mean of per-trial rates, and the
+    # count that says whether a point means anything is kept
+    counters = {}
+    for name in receivers:
+        counters[name] = ErrorCounter()
+
     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
+        # the two propagations of a trial share their seed, so they
+        # carry the same symbols and the same amplifier noise and
+        # differ only by the fibre's Kerr term
         fields = {}
         for use_only_linear in (True, False):
-            channel.seed(index * N_trial + trial)
+            channel.seed(point_seed + trial)
             channel.set_params(launch__gain=amp,
                                link__use_only_linear=use_only_linear)
             fields[use_only_linear] = channel(N_s)
-        symbols, reference = channel.tap("signal_tx"), channel.tap("data_tx")
+        symbols = channel.tap("signal_tx")
+        reference = channel.tap("data_tx")
 
+        # each receiver decides the field that matches its claim: the
+        # bound receives the linear one
         for name, (steps, linear_only) in receivers.items():
             receiver = get_receiver(steps, linear_only, gain=1 / amp,
                                     reference=symbols)
-            bound = name == "amplifier noise only"
-            detected = receiver(fields[bound])
-            elapsed[name] += receiver.elapsed_
-            estimate = receiver.tap("phase")
-            snr[name][index] += compute_effective_snr(symbols, estimate) / N_trial
-            counters[name][index].update(reference, detected)
+            detected = receiver(fields[name == "amplifier noise only"])
+            snr[name][index] += compute_effective_snr(
+                symbols, receiver.tap("phase")) / N_trial
+            counters[name].update(reference, detected)
+            times[name][index] += receiver.elapsed_
 
-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):
+    for name, counter in counters.items():
         ser[name][index] = counter.rate
+        errors[name][index] = counter.n_errors
 
-# The sweep result, written down once. It is printed here and drawn
-# below from the same object; transposed because six receiver names as
-# column headers would make the table 160 characters wide.
+# --- results: tables and figures -------------------------------------
+# The SNR averages in linear and is read in dB.
+snr_dB = {}
+for name in receivers:
+    snr_dB[name] = 10 * np.log10(snr[name])
+
+# transposed: six receiver names as column headers would make the table
+# 160 characters wide
 snr_data = {"x": dBm_list, "curves": snr_dB}
 print_data(snr_data, xlabel="launch power [dBm]",
            ylabel="effective SNR [dB]", transpose=True)
@@ -178,8 +196,9 @@ def get_receiver(steps, linear_only, *, gain, reference):
 print("\nreceiver                  best SNR   at power    total time")
 for name, values in snr_dB.items():
     best = int(np.argmax(values))
+    total = float(np.sum(times[name]))
     print(f"{name:24s} {values[best]:6.2f} dB {dBm_list[best]:6.1f} dBm "
-          f"{elapsed[name]:9.1f} s")
+          f"{total:9.1f} s")
 
 reference = snr_dB["dispersion compensation"]
 best = int(np.argmax(reference))
@@ -191,31 +210,26 @@ def get_receiver(steps, linear_only, *, gain, reference):
 # 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.
+symbols_per_point = N_s * N_trial
 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")
+for name in receivers:
+    best = int(np.argmax(snr_dB[name]))
+    count = int(errors[name][best])
+    print(f"{name:24s} {count:7d} over {symbols_per_point} 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
 # receiver that removes part of it.
 fine_powers = dbm_to_watt(np.linspace(dBm_list[0], dBm_list[-1], 200))
-# An effective SNR in dB is not an error rate -- a linear ordinate, six
-# receivers and one closed form -- so it is drawn here rather than routed
-# through plot_error_rate, whose name would then be describing something
-# it is not.
-fig1, ax = plt.subplots()
-for name, values in snr_dB.items():
-    ax.plot(dBm_list, values, "o-", fillstyle="none", label=name)
+# The same object the table printed, drawn, with the closed form on top.
+ax = plot_data(snr_data, xlabel="launch power [dBm]",
+               ylabel="effective SNR [dB]", marker="o", fillstyle="none")
 ax.plot(watt_to_dbm(fine_powers),
         10 * np.log10(gn_model_snr(ase_W, eta, fine_powers)), "k:",
         label="GN model, single polarization")
-ax.set_xlabel("launch power [dBm]")
-ax.set_ylabel("effective SNR [dB]")
 ax.set_title(f"{N_span} x {L_span} km, {constellation.order}-"
              f"{constellation.family} at {R_s / 1e9:.0f} GBd")
-ax.grid(True)
 ax.legend()
 plt.tight_layout()
 plt.savefig(f"{img_dir}/nli_simulation_fig1.png")
diff --git a/examples/optical/gn_model.py b/examples/optical/gn_model.py
index 6170a1d..1ca5c27 100644
--- a/examples/optical/gn_model.py
+++ b/examples/optical/gn_model.py
@@ -125,9 +125,9 @@ def get_chain(order=16, power_W=1e-3):
                SymbolMapper(Constellation("QAM", order), name="tx")])
     return Sequential([
         *source,
+        Amplifier(launch_amplitude(power_W, polarizations=2), name="launch"),
         Upsampler(OS, scale=np.sqrt(OS)),
         SRRCFilter(ROLLOFF, OS, N_h=40, method="fft"),
-        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,
diff --git a/examples/optical/one_shot_NLI.py b/examples/optical/one_shot_NLI.py
index 1511b5b..122e386 100644
--- a/examples/optical/one_shot_NLI.py
+++ b/examples/optical/one_shot_NLI.py
@@ -6,7 +6,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from comnumpy.core import Sequential
+from comnumpy.core import Sequential, plot_iq
 from comnumpy.core.compensators import DataAidedPhaseCompensator
 from comnumpy.core.filters import SRRCFilter
 from comnumpy.core.generators import SymbolGenerator
@@ -62,9 +62,9 @@ def get_chain(n_spans, *, steps=1, linear_only=True):
     return Sequential([
         SymbolGenerator(constellation.order, name="data_tx"),
         SymbolMapper(constellation, name="signal_tx"),
+        Amplifier(amp),
         Upsampler(oversampling_sim, scale=np.sqrt(oversampling_sim)),
         SRRCFilter(rolloff, oversampling_sim, method="fft"),
-        Amplifier(amp),
         FiberLink(N_spans=n_spans, L_span=L_span, StPS=StPS, NF_dB=NF_dB,
                   fs=fs, fiber=fiber, name="link"),
         Downsampler(oversampling_ratio, use_filter=True, name="rx_field"),
@@ -155,32 +155,22 @@ def score(chain):
     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")
 
-# An IQ plane is two matplotlib calls; style.apply gives it the labels,
-# the grid and the equal aspect ratio a constellation has to be read on.
-alphabet = constellation.alphabet
-fig1, ax1 = plt.subplots()
-ax1.plot(np.real(received_field), np.imag(received_field), ".")
-ax1.set_title(f"received field ({dBm} dBm, {N_span} spans)")
-style.apply(ax1, "iq")
+plot_iq(received_field,
+        title=f"received field ({dBm} dBm, {N_span} spans)")
 plt.savefig(f"{img_dir}/one_shot_nli_fig1.png")
 
 fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 4))
 for ax, n_spans in zip(axes, (1, 10, 25), strict=True):
-    symbols = estimates[n_spans]
-    ax.plot(np.real(symbols), np.imag(symbols), ".", label="symbols")
-    ax.plot(np.real(alphabet), np.imag(alphabet), "kx", markersize=9,
-            label="alphabet")
+    plot_iq(estimates[n_spans], reference=constellation, ax=ax)
     ax.set_title(f"after {n_spans} span{'s' if n_spans > 1 else ''}, "
                  f"SNR {snr_per_span[n_spans]:.1f} dB")
-    style.apply(ax, "iq")
 plt.tight_layout()
 plt.savefig(f"{img_dir}/one_shot_nli_fig2.png")
 
-measured = np.zeros(len(spans))
-ase_only = np.zeros(len(spans))
-for index, n_spans in enumerate(spans):
-    measured[index] = snr_per_span[n_spans]
-    ase_only[index] = snr_ase_only[n_spans]
+# the dicts are keyed by span count in sweep order, so their values
+# are already the curves
+measured = np.array(list(snr_per_span.values()))
+ase_only = np.array(list(snr_ase_only.values()))
 
 # An effective SNR is not an error rate: linear ordinate, two curves,
 # nothing a helper would say better than the four lines that draw it.
@@ -223,10 +213,6 @@ def score(chain):
          f"dispersion compensation\nSNR {snr_per_span[N_span]:.2f} dB"),
         (axes[1], back_propagated.tap("phase"),
          f"digital back-propagation\nSNR {snr_dbp:.2f} dB")]:
-    ax.plot(np.real(symbols), np.imag(symbols), ".", label="symbols")
-    ax.plot(np.real(alphabet), np.imag(alphabet), "kx", markersize=9,
-            label="alphabet")
-    ax.set_title(name)
-    style.apply(ax, "iq")
+    plot_iq(symbols, reference=constellation, title=name, ax=ax)
 plt.tight_layout()
 plt.savefig(f"{img_dir}/one_shot_nli_fig4.png")
diff --git a/examples/simple/monte_carlo_awgn.py b/examples/simple/monte_carlo_awgn.py
index f218372..5660ff2 100644
--- a/examples/simple/monte_carlo_awgn.py
+++ b/examples/simple/monte_carlo_awgn.py
@@ -41,14 +41,16 @@ def get_chain(constellation):
 
 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:
-by_hand = []
-for snr_dB in snr_dB_list[::8]:
+# A Monte Carlo simulation is three parts, always in this order: the
+# metric pre-allocated, the loop that fills it, the display at the end.
+# Written out, at three SNR values:
+probed_snr_dB = snr_dB_list[::8]
+ser_by_hand = np.zeros(len(probed_snr_dB))
+for index, snr_dB in enumerate(probed_snr_dB):
     chain.seed(1)
     chain.set_params(awgn_channel__snr_dB=snr_dB)
     detected = chain(N)                  # run first: the tap is filled by it
-    by_hand.append(compute_ser(chain.tap("tx"), detected))
+    ser_by_hand[index] = compute_ser(chain.tap("tx"), detected)
 
 # monte_carlo (decision D35) is that loop, and nothing more: the same
 # reconfigure-reseed-run-collect, over every point, in one call.
@@ -58,8 +60,8 @@ def get_chain(constellation):
 
 # An abscissa and one series per curve is what a sweep produces, and
 # print_data is what shows it: the same dictionary would draw the figure.
-print_data({"x": snr_dB_list[::8],
-            "curves": {"loop": np.array(by_hand),
+print_data({"x": probed_snr_dB,
+            "curves": {"loop": ser_by_hand,
                        "monte_carlo": ser_array[::8]}},
            xlabel="SNR [dB]", ylabel="SER")
 
@@ -101,13 +103,14 @@ def qfunc(x):
 ebn0_dB_list = np.arange(0, 25, 2)
 N_compare = 100000
 fine_dB_list = np.arange(0, 25, 0.05)
+orders = np.array([4, 16, 64, 256])
+bits_per_symbol = np.log2(orders).astype(int)     # vectorized, not a loop
 measured_ber = {}
 theory_ber = {}
-required_rows = {"bits per symbol": [],
-                 "Eb/N0 at BER=1e-3 [dB]": []}
-for order in [4, 16, 64, 256]:
-    other = Constellation("QAM", order)
-    bits = int(np.log2(order))
+required_dB = np.zeros(len(orders))
+for index, order in enumerate(orders):
+    other = Constellation("QAM", int(order))
+    bits = int(bits_per_symbol[index])
     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,
@@ -120,14 +123,14 @@ def qfunc(x):
     # 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["bits per symbol"].append(bits)
-    required_rows["Eb/N0 at BER=1e-3 [dB]"].append(required_dB)
+    required_dB[index] = np.interp(
+        -3.0, np.log10(fine_ber[above_underflow])[::-1],
+        fine_dB_list[above_underflow][::-1])
 
 print()
-print_data({"x": [4, 16, 64, 256], "curves": required_rows},
+print_data({"x": orders,
+            "curves": {"bits per symbol": bits_per_symbol,
+                       "Eb/N0 at BER=1e-3 [dB]": required_dB}},
            xlabel="QAM order")
 
 ax = plot_error_rate(ebn0_dB_list, measured_ber, theory=theory_ber,
diff --git a/examples/simple/one_shot_awgn.py b/examples/simple/one_shot_awgn.py
index b38b8b8..1a1352c 100644
--- a/examples/simple/one_shot_awgn.py
+++ b/examples/simple/one_shot_awgn.py
@@ -8,7 +8,7 @@
 from scipy.stats import norm
 
 from comnumpy import style
-from comnumpy.core import Sequential
+from comnumpy.core import Sequential, plot_iq
 from comnumpy.core.channels import AWGN
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolDemapper, SymbolMapper
@@ -89,11 +89,7 @@
 # --- 4. the figures ---------------------------------------------------
 fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4.2))
 constellation.plot(ax=axes[0])
-received = chain.tap("awgn")
-axes[1].plot(np.real(received), np.imag(received), ".")
-axes[1].set_title(f"received at {snr_dB} dB")
-style.apply(axes[1], "iq")
-axes[0].grid(True)
+plot_iq(chain.tap("awgn"), title=f"received at {snr_dB} dB", ax=axes[1])
 plt.tight_layout()
 plt.savefig(f"{img_dir}/first_simulation_fig1.png")
 
diff --git a/examples/simple/one_shot_cfo_iq_compensation.py b/examples/simple/one_shot_cfo_iq_compensation.py
index 5837c97..eea7c85 100644
--- a/examples/simple/one_shot_cfo_iq_compensation.py
+++ b/examples/simple/one_shot_cfo_iq_compensation.py
@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 
 from comnumpy import style
-from comnumpy.core import Sequential
+from comnumpy.core import Sequential, plot_iq
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper, SymbolDemapper
 from comnumpy.core.impairments import IQImbalance, CFO
@@ -48,17 +48,11 @@
 # print metric and plot
 print(f"after: SER={ser_after}")
 
-# Four IQ planes, one per stage: a scatter of real against imaginary,
-# and style.apply for the labels, the grid and the equal aspect ratio.
 for tap, name in [("awgn", "received signal"),
                   ("gsop", "after GSOP"),
                   ("cfo_comp", "after GSOP+CFO comp"),
                   ("phase_comp", "after GSOP + CFO comp + phase correction")]:
-    symbols = chain.tap(tap)
-    _, ax = plt.subplots()
-    ax.plot(np.real(symbols), np.imag(symbols), ".")
-    ax.set_title(name)
-    style.apply(ax, "iq")
+    plot_iq(chain.tap(tap), title=name)
 
 # show evolution of the angular frequency estimate
 w0_history = chain["cfo_comp"].history
diff --git a/examples/simple/one_shot_phase_compensation.py b/examples/simple/one_shot_phase_compensation.py
index dbdfa03..f35aee8 100644
--- a/examples/simple/one_shot_phase_compensation.py
+++ b/examples/simple/one_shot_phase_compensation.py
@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 
 from comnumpy import style
-from comnumpy.core import Sequential
+from comnumpy.core import Sequential, plot_iq
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper, SymbolDemapper
 from comnumpy.core.channels import AWGN
@@ -49,9 +49,5 @@
 
 for tap, name in [("awgn", "received data"),
                   ("phase_compensation", "after phase correction")]:
-    symbols = chain.tap(tap)
-    _, ax = plt.subplots()
-    ax.plot(np.real(symbols), np.imag(symbols), ".")
-    ax.set_title(name)
-    style.apply(ax, "iq")
+    plot_iq(chain.tap(tap), title=name)
 plt.show()
diff --git a/examples/simple/one_shot_srrc_awgn.py b/examples/simple/one_shot_srrc_awgn.py
index 30a6b40..d90e1d6 100644
--- a/examples/simple/one_shot_srrc_awgn.py
+++ b/examples/simple/one_shot_srrc_awgn.py
@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 
 from comnumpy import style
-from comnumpy.core import Sequential, plot_spectrum
+from comnumpy.core import Sequential, plot_spectrum, plot_iq
 from comnumpy.core.generators import SymbolGenerator
 from comnumpy.core.mappers import SymbolMapper, SymbolDemapper
 from comnumpy.core.processors import Upsampler, Downsampler, DataExtractor
@@ -44,11 +44,7 @@
 
 # plot the extracted signals
 plot_spectrum(chain.tap("awgn_channel"), title="received signal")
-symbols = chain.tap("extractor")
-_, ax = plt.subplots()
-ax.plot(np.real(symbols), np.imag(symbols), ".")
-ax.set_title("after SRRC+downsampling+extractor")
-style.apply(ax, "iq")
+plot_iq(chain.tap("extractor"), title="after SRRC+downsampling+extractor")
 
 # plot error distribution
 N_min = np.min([len(data_tx), len(y)])
diff --git a/tests/test_examples_run.py b/tests/test_examples_run.py
index 8aa9840..51b3cd9 100644
--- a/tests/test_examples_run.py
+++ b/tests/test_examples_run.py
@@ -55,11 +55,13 @@
 # row. The old figure is not reproducible here and no change to the
 # script explains it, so rather than keep a script out of the smoke test
 # on a number nobody can obtain, it goes back in.
+#
+# `mimo/monte_carlo_simulation_1.py` left the same way on 2026-08-13:
+# listed at 32 s, measured twice at 4.4 s after its Experiment rewrite.
 SLOW = {
-    "mimo/monte_carlo_simulation_1.py":
-        "Monte-Carlo BER sweep, 4 detectors -- 32 s measured 2026-08-09",
     "mimo/monte_carlo_simulation_2.py":
-        "Monte-Carlo BER sweep, OSIC variants -- 143 s measured 2026-08-09",
+        "Monte-Carlo BER sweep, OSIC variants -- 80 s measured 2026-08-13 "
+        "after the Experiment rewrite",
     "mimo/run_all_scripts.py":
         "runner, not an example: re-executes the other mimo scripts -- "
         "~350 s by construction",
@@ -68,7 +70,8 @@
         "matrix, plus a runtime sweep -- 33 s wall measured 2026-08-12",
     "optical/CD_compensation_part1.py":
         "BER vs SNR for 3 modulations, 200 000 symbols through a 249-tap "
-        "compensator -- 227 s measured 2026-08-09",
+        "compensator -- 332 s measured 2026-08-13 after the Experiment "
+        "rewrite, alongside a docs build",
     "optical/NLI_simulation.py":
         "eight launch powers x four trials x two fibre configurations of "
         "split-step propagation, then six receivers on each -- 144 s "