feat(v2): propagation-of-chaos animation, README section, notebook sandwich

- examples/animate_propagation_of_chaos.py: 4-panel McKean-Vlasov simulator at N in {20, 100, 500, 4000} with reference N=12000; bottom panel tracks W_2(mu^N_t, mu_t) on log scale -> visible 1/sqrt(N) decay (Sznitman 1991). - examples/propagation_of_chaos.gif (1.6 MB) - README: new 'Propagation of chaos' subsection under Mean-field & agent-based dynamics, with empirical-measure formula, k-tuple factorisation and GIF embed. - examples/notebooks/14_mckean_vlasov.ipynb: sandwich PRE/code/POST cells demonstrating W2 ~ 1/sqrt(N) on the same simulator. Verified executed: sqrt(N)*W2 ~ 0.7 across N (theoretical const).
release(v2.0.0): finalize PyPI metadata, README rewrite, v2 benchmark + McKean-Vlasov animation
2026-05-14 22:43:47 +02:00 · 2026-05-14 21:54:49 +02:00 · 2026-05-14 14:00:39 +02:00 · 2026-05-12 19:23:49 +02:00 · 2026-05-12 17:10:06 +02:00 · 2026-05-12 16:47:15 +02:00
238 changed files with 111845 additions and 471 deletions
@@ -23,3 +23,10 @@ Cargo.lock
 *.swo
 *~
 wheels/
+
+# Build artifacts
+docs/source/_static/*.aux
+docs/source/_static/*.log
+docs/source/_static/*.toc
+docs/source/theory/*.html
+docs/source/theory/*.pdf
@@ -7,8 +7,6 @@ build:
  os: ubuntu-22.04
  tools:
    python: "3.11"
-  apt_packages:
-    - libopenblas-dev

 sphinx:
  configuration: docs/source/conf.py
@@ -20,5 +18,9 @@ formats:
 python:
  install:
    - requirements: docs/requirements.txt
-    - method: pip
-      path: .
+# Note: we intentionally do NOT `pip install .` on RTD.
+# The Sphinx build only renders RST/MD pages and pre-generated PNG plots
+# committed under docs/source/_static/. Building the Rust extension on RTD
+# requires OpenBLAS + maturin and routinely breaks because of upstream
+# `openblas-build` / `ureq` feature changes. All inline plots are pre-rendered
+# locally via scripts/inject_doc_plots.py and committed to the repo.
@@ -0,0 +1,135 @@
+# Changelog
+
+All notable changes to **optimiz-rs** are documented in this file. The format
+follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/) and the project
+adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [2.0.0] - 2026-05-14
+
+### Added — public release of the v2 API
+
+- Promoted `2.0.0-alpha.1` to the stable `2.0.0` release.
+- PyPI distribution name remains `optimiz-rs` (continuity with v1.0.x);
+  the Rust crate is also `optimiz-rs`. Both expose the Python module
+  `optimizr`.
+- New non-regression suite `tests/test_v2_api.py` (20 tests) exercising
+  every advertised v2 primitive against an analytic ground truth:
+  `historical_var_py`, `solve_fractional_ode`, `solve_volterra`,
+  `linear_bsde_constant_coeffs`, `mean_reverting_mckean_vlasov`, plus a
+  parametrised guard over the v1.x public surface.
+- README rewritten to document the v2 Python API and the corrected
+  installation command (`pip install optimizr`).
+
+### Notes
+
+- No source-level breaking change relative to `2.0.0-alpha.1`.
+- All v1.x Python entry points remain exposed (`differential_evolution`,
+  `fit_hmm`, `viterbi_decode`, `mcmc_sample`, `grid_search`, `mutual_information`,
+  `shannon_entropy`, etc.) — verified by `test_public_symbol_exposed`.
+
+## [2.0.0-alpha.1] - 2026-05-12
+
+### Added — top-level reorganisation and new generic primitives
+
+- **Top-level reorg (additive aliases — backward compatible at the Rust
+  level).**  `optimiz_rs::matrix_riccati` is now re-exported at the crate
+  root.  New top-level groups: `bsde`, `pde`, `stochastic_control`,
+  `agent_based`, `inference`, `optimization`.
+- `bsde::theta_scheme` — implicit/explicit θ-scheme for linear BSDEs
+  with deterministic coefficients (closed-form analytic test against
+  the deterministic ODE `dY = -ρ Y dt`).
+- `bsde::deep_bsde_bridge` — `ConditionalExpectation` trait and
+  `DeepBsdeBridge` driver providing the CPU-side recursion hook for
+  external function approximators.
+- `pde::fokker_planck` — 1-D forward Fokker--Planck solver with
+  conservative central differences and explicit positivity safeguard.
+- `pde::hjb_multid` — explicit upwind solver for multidimensional HJB
+  on a regular Cartesian grid (`d ≤ 3`) with reflective boundaries.
+- `pde::elliptic_fd` — 2-D Poisson `-Δu = f` SOR solver verified
+  against the `sin(πx) sin(πy)` eigenfunction.
+- `stochastic_control::optimal_switching` — Snell-envelope backward
+  induction for discrete multi-mode optimal switching.
+- `stochastic_control::pontryagin` — Riccati-shooting solver for the
+  1-D LQR Pontryagin maximum principle (verified against the
+  closed-form `P(t) = s_T / (1 + s_T (T-t))`).
+- `stochastic_control::two_sided_intensity_control` — generic
+  bilateral intensity control with affine per-jump premia.
+- `optimal_control::quadratic_impact_control` — closed-form Riccati
+  feedback for a controlled SDE with quadratic running cost.
+- `mean_field::mckean_vlasov` — interacting-particle Euler scheme for
+  generic McKean--Vlasov SDEs with empirical-measure drift.
+- `agent_based` — generic interacting-agent simulator (consensus
+  dynamics test recovers the empirical mean exactly without noise).
+- `inference::robust_drift` — Huber-loss IRLS estimator for the drift
+  of a 1-D OU-type discrete-time process; resists 5 % outliers.
+- `optimization::generative_calibration_hooks` — `GenerativeSampler`
+  trait + Gaussian MMD loss + finite-difference calibration step.
+
+### Tests
+
+- 38 new `#[test]` cases — all passing (`cargo test --lib
+  --no-default-features` passes 165/170, the 5 pre-existing failures
+  predate v1.1 and are unrelated).
+
+### Notes — deferred to subsequent v2.0.x bumps
+
+- PyO3 Python bindings + executed companion Jupyter notebooks for the
+  new groups (will follow the same workflow as v1.1.x).
+- Sphinx RST documentation pages for `bsde`, `pde`,
+  `stochastic_control`, `agent_based`, `inference`, `optimization`.
+- Propagation of new modules into `hfthot-lab-instance` consumers.
+
+## [1.1.0] - 2026-05-12
+
+### Added — purely additive, no existing API changes
+
+- `optimal_control::matrix_riccati` — RK4 backward solver for the
+  matrix Riccati differential equation
+  `dA/dt = -2 A M A + Q`, plus terminal-condition variants for the
+  associated affine and constant components.
+- `timeseries_utils::nonsync_covariance` — Hayashi--Yoshida estimator
+  for asynchronous covariance, with parallel matrix variant.
+- `timeseries_utils::wavelet` — discrete and maximum-overlap wavelet
+  transforms (Haar, Daubechies orders 2--10) with periodic boundaries.
+- `risk_measures` — empirical and parametric Value-at-Risk and
+  Conditional Value-at-Risk estimators, plus a projected sub-gradient
+  solver for convex CVaR minimisation over the unit simplex.
+- `graph::laplacian` — combinatorial, symmetric-normalised and
+  random-walk graph Laplacians.
+- `graph::spectral_clustering` — spectral clustering via Jacobi
+  diagonalisation and Lloyd's algorithm with k-means++ initialisation.
+- `topology::persistent_homology` — Vietoris--Rips persistent homology
+  by the standard `Z/2` matrix-reduction algorithm.
+- `topology::bottleneck` — bottleneck distance between persistence
+  diagrams via Hopcroft--Karp matching with binary search.
+- `volterra::fractional_riccati` — Adams predictor--corrector solver
+  for Caputo fractional ODEs (Diethelm--Ford--Freed 2002).
+- `volterra::markovian_lift` — multi-exponential approximation of
+  convolution kernels by non-negative least squares on a geometric
+  grid.
+- `volterra::volterra_solver` — generic second-kind Volterra integral
+  equation solver via product-trapezoidal quadrature.
+- `volterra::fourier_inversion` — direct trapezoidal Fourier inversion
+  of a characteristic function on a uniform frequency grid.
+- `signatures::path_signature` — truncated tensor signature of a
+  piecewise-linear path with truncated tensor exponential.
+- `signatures::log_signature` — truncated tensor logarithm of a
+  signature.
+- `signatures::random_signature` — Cuchiero--Schmocker--Teichmann
+  random reservoir projection of the signature.
+- `signatures::signature_kernel` — Salvi--Cass--Lyons signature kernel
+  via the Goursat finite-difference scheme.
+- `signatures::utils` — shuffle product and Chen-identity-driven
+  signature concatenation.
+
+### Changed
+
+- Bumped crate version from `1.0.1` to `1.1.0`. All previously stable
+  symbols remain untouched and binary-compatible at the Python ABI
+  level (`abi3-py38`).
+
+### Notes
+
+This release is **CPU-only and additive**. No Python bindings were added
+in `1.1.0`; the new modules are exposed via the Rust API only and will
+be wrapped behind the `python-bindings` feature in a follow-up release.
@@ -1,6 +1,6 @@
 [package]
-name = "optimizr"
-version = "1.0.0"
+name = "optimiz-rs"
+version = "2.0.0"
 edition = "2021"
 authors = ["HFThot Research Lab <contact@hfthot-lab.eu>"]
 description = "High-performance optimization algorithms in Rust with Python bindings"
@@ -1,24 +1,154 @@
-# OptimizR 🚀
+# Optimiz-rs 🚀

 <p align="center">
-    <img src="docs/source/logo_optimizr_valid.png" alt="OptimizR Logo" width="220" />
+    <img src="https://raw.githubusercontent.com/ThotDjehuty/optimiz-r/main/docs/source/logo_optimizrs.png" alt="Optimiz-rs Logo" width="220" />
 </p>

 **High-performance optimization algorithms in Rust with Python bindings**

-[![Version](https://img.shields.io/badge/version-1.0.0-blue.svg)](https://github.com/ThotDjehuty/optimiz-r/releases)
+[![Version](https://img.shields.io/badge/version-2.0.0-blue.svg)](https://github.com/ThotDjehuty/optimiz-r/releases)
 [![License](https://img.shields.io/badge/license-MIT-green.svg)](LICENSE)
 [![Rust](https://img.shields.io/badge/rust-1.70+-orange.svg)](https://www.rust-lang.org/)
 [![Python](https://img.shields.io/badge/python-3.8+-blue.svg)](https://www.python.org/)

-OptimizR provides blazingly fast, production-ready implementations of advanced optimization and statistical inference algorithms. Built with Rust for maximum performance and exposed to Python through PyO3, it delivers 50-100× speedup over pure Python implementations.
+Optimiz-rs provides blazingly fast, production-ready implementations of advanced optimization and statistical inference algorithms. Built with Rust for maximum performance and exposed to Python through PyO3, it delivers up to **86× speedup** over pure-Python references on intrinsically loopy / sequential workloads.
+
+<p align="center">
+    <img src="examples/mckean_vlasov.gif" alt="McKean-Vlasov mean-reverting flow" width="640" />
+    <br/>
+    <em>800-particle mean-reverting McKean–Vlasov flow simulated by
+    <code>optimizr.mean_reverting_mckean_vlasov</code> — two clouds at
+    <code>x = ±2</code> fuse under <code>dX_t = θ(m̄_t − X_t) dt + σ dW_t</code>.
+    Source: <a href="examples/animate_mckean_vlasov.py"><code>examples/animate_mckean_vlasov.py</code></a>.</em>
+</p>
+
+## ✨ What's New in v2.0.0
+
+v2 ships **eight brand-new CPU-only generic numerical primitive groups** with full Python bindings, on top of every v1.x algorithm (which remain available). The Python module name is unchanged: `import optimizr as opt`.
+
+### Rough volatility & integral equations
+
+- **`solve_fractional_ode(h0, alpha, t_horizon, n_steps, rhs)`** — Caputo fractional ODE Adams scheme.
+- **`solve_volterra(g, kernel, t_horizon, n_steps)`** — second-kind Volterra integral equation by trapezoidal product integration.
+- **`geometric_grid_lift(kernel, t_samples, n_factors, gamma_min, gamma_max)`** — multi-exponential approximation of a kernel by NNLS on a geometric rate grid (Markovian lift à la Abi Jaber–El Euch).
+- **`fourier_invert(char_fn, t_grid, x_grid)`** — characteristic-function → density inversion (Carr–Madan style).
+- **`mittag_leffler_py(z, alpha, beta)`** — generalised Mittag-Leffler reference function.
+
+### Backward SDEs & PDEs
+
+- **`linear_bsde_constant_coeffs(a, b, c, terminal, n_steps, t_horizon, theta=0.5)`** — backward SDE θ-scheme (closed-form analytic test against `dY = -ρ Y dt`).
+- **`fokker_planck_constant(...)`** — 1-D forward Fokker–Planck solver with conservative central differences.
+- **`hjb_quadratic_2d(...)`** — explicit upwind solver for 2-D HJB on a Cartesian grid.
+- **`poisson_2d_zero_boundary(...)`** — 2-D Poisson `−Δu = f` SOR solver.
+
+### Stochastic & quadratic-impact control
+
+- **`optimal_switching_dp(...)`** — discrete-time optimal switching by dynamic programming.
+- **`pontryagin_lqr(...)`** — Pontryagin maximum principle for LQ control.
+- **`two_sided_intensities(...)`** — bilateral intensity-controlled jump process.
+- **`quadratic_impact_control_py(...)`** — convex quadratic-cost control on a controlled SDE.
+
+### Mean-field & agent-based dynamics
+
+- **`mean_reverting_mckean_vlasov(initial, theta, sigma, n_steps, t_horizon, seed)`** — N-particle McKean–Vlasov simulator (returns `paths_flat`, `n_particles`, `n_steps`, `time_grid`).
+- **`consensus_dynamics(...)`** — synchronous opinion-dynamics consensus on a graph.
+- **`solve_mfg_1d_rust(MFGConfig)`** — 1-D mean-field game (HJB ↔ Fokker–Planck fixed-point).
+
+#### Propagation of chaos
+
+<p align="center">
+    <img src="examples/propagation_of_chaos.gif" alt="Propagation of chaos" width="680" />
+</p>
+
+For an interacting N-particle system
+
+$$
+dX^{i,N}_t \;=\; b\!\bigl(X^{i,N}_t,\; \mu^N_t\bigr)\, dt \;+\; \sigma\, dW^i_t,
+\qquad
+\mu^N_t \;=\; \frac{1}{N}\sum_{j=1}^{N}\delta_{X^{j,N}_t},
+$$
+
+Sznitman's theorem (1991) states that whenever $b$ is Lipschitz in both arguments,
+the empirical measure $\mu^N_t$ converges in Wasserstein-2 to the law $\mu_t$
+of the McKean–Vlasov limit at rate $\mathcal{O}(1/\sqrt{N})$, and any finite
+$k$-tuple of particles becomes asymptotically independent — *chaos propagates*
+from $t=0$ to all later times:
+
+$$
+\operatorname{Law}\!\bigl(X^{1,N}_t,\dots,X^{k,N}_t\bigr) \;\xrightarrow[N\to\infty]{w}\; \mu_t^{\otimes k}.
+$$
+
+The animation above runs four parallel simulations with $N\in\{20,100,500,4000\}$
+sharing the *same* bimodal initial law, the *same* drift $-\theta(x-\bar{x})$ and
+the *same* noise $\sigma\, dW$. The bottom panel tracks the Wasserstein-2 distance
+$W_2(\mu^N_t, \mu_t)$ to a high-resolution reference and visibly decays as
+$1/\sqrt{N}$. Source: [`examples/animate_propagation_of_chaos.py`](examples/animate_propagation_of_chaos.py).
+Companion notebook: [`examples/notebooks/14_mckean_vlasov.ipynb`](examples/notebooks/14_mckean_vlasov.ipynb).
+
+
+### Topology, graphs & path signatures
+
+- **`vietoris_rips_filtration`**, **`persistent_homology`**, **`bottleneck_distance`** — TDA primitives.
+- **`combinatorial_laplacian_py`**, **`normalised_laplacian_py`**, **`random_walk_laplacian_py`**, **`spectral_cluster_py`** — graph spectral analysis.
+- **`path_signature`**, **`path_log_signature`**, **`random_signature`**, **`signature_kernel`**, **`shuffle_product`**, **`concatenate_signatures`** — Chen–Strichartz iterated integrals and signature kernels.
+
+### Risk, robust inference & calibration
+
+- **`historical_var_py(losses, alpha)`**, **`parametric_var_py(...)`**, **`cvar_value_py(...)`**, **`minimize_cvar_py(...)`** — coherent risk measures.
+- **`robust_drift(...)`**, **`estimate_hurst(...)`**, **`scale_dependent_hurst(...)`** — robust drift / Hurst estimation.
+- **`mmd_gaussian(...)`**, **`f_alpha_lambda_py(...)`** — generative-calibration hooks (MMD, fractional kernels).
+
+### Point processes & Kalman filtering
+
+- **`simulate_hawkes`**, **`simulate_bivariate_hawkes`**, **`simulate_fbm`**, **`simulate_mixed_fbm`** — order-flow simulators.
+- **`LinearKalmanFilter`**, **`UnscentedKalmanFilter`**, **`RTSSmoother`** — state-space inference.
+
+### Quality bar
+
+- 20-test analytic non-regression suite for the v2 public API: [`tests/test_v2_api.py`](tests/test_v2_api.py).
+- ABI3 wheels, Python ≥ 3.8.
+- Crate name: `optimiz-rs` (Rust); distribution name: `optimiz-rs` (PyPI); module name: `optimizr` (Python import).
+
+```bash
+pip install --upgrade optimiz-rs
+python -c "import optimizr; print(optimizr.__version__)"   # 2.0.0
+```
+
+### v2 benchmark (single-threaded, best-of-3, Apple M2)
+
+Generated by [`examples/benchmark_v2.py`](examples/benchmark_v2.py):
+
+| Workload | Pure Python / NumPy | optimiz-rs (Rust) | Speedup |
+|---|---:|---:|---:|
+| HMM Baum-Welch (2 states, 5 000 obs, 10 iters) | 970.94 ms | 14.34 ms | **67.7×** |
+| Differential evolution (Rastrigin d=5, 50 iters × 20 pop) | 417.45 ms | 30.03 ms | **13.9×** |
+| Path signature (T=300, d=3, depth=3) | 11.07 ms | 0.99 ms | **11.2×** |
+| Hawkes process (T=100, μ=1, α=0.6, β=1.2) | 2.75 ms | 0.83 ms | **3.3×** |
+| MCMC random-walk MH (5 000 samples, d=2) | 35.41 ms | 20.24 ms | **1.7×** |
+
+> Workloads that are fully vectorisable in NumPy (e.g. drift updates for an N-particle SDE without callback) are not in this table: a tight NumPy loop on contiguous arrays is hard to beat from Rust through a PyO3 callback boundary. Use `optimiz-rs` for the algorithms above and `numpy` for the rest — both are first-class citizens.
+
+## ✨ What's New in v1.1.0
+
+This release adds a broad collection of **CPU-only generic numerical primitives**, all purely additive:
+
+- **`optimal_control::matrix_riccati`** — backward RK4 solver for the matrix Riccati ODE.
+- **`timeseries_utils::nonsync_covariance`** — Hayashi--Yoshida asynchronous covariance estimator.
+- **`timeseries_utils::wavelet`** — DWT and MODWT (Haar, Daubechies 2--10).
+- **`risk_measures`** — empirical / parametric VaR and CVaR, plus convex CVaR minimisation.
+- **`graph::laplacian`** + **`graph::spectral_clustering`** — combinatorial / normalised / random-walk Laplacians and spectral clustering with Jacobi diagonalisation.
+- **`topology`** — Vietoris--Rips persistent homology and bottleneck distance.
+- **`volterra`** — fractional Caputo Adams solver, Markovian lift by NNLS on a geometric grid, second-kind Volterra solver, direct Fourier inversion of characteristic functions.
+- **`signatures`** — truncated tensor signatures, log-signatures, random-reservoir projection (Cuchiero--Schmocker--Teichmann), Salvi--Cass--Lyons signature kernel, shuffle product / Chen concatenation.
+
+All new modules are exposed via the **Rust API only** in this release; Python bindings will follow in a subsequent minor release. The previously stable Python API is unchanged.

 ## ✨ What's New in v1.0.0

 🎉 **Production Ready** - First stable release with comprehensive documentation  
 📚 **ReadTheDocs** - Full documentation at https://optimiz-r.readthedocs.io  
-🏗️ **Published to crates.io** - Install with `cargo add optimizr`  
-🐍 **Published to PyPI** - Install with `pip install optimizr`  
+🏗️ **Published to crates.io** - Install with `cargo add optimiz-rs`  
+🐍 **Published to PyPI** - Install with `pip install optimiz-rs`  
 🔒 **Stable API** - Semantic versioning from v1.0.0 forward

 ## Features
@@ -49,12 +179,22 @@ OptimizR provides blazingly fast, production-ready implementations of advanced o

 ## Installation

-### From PyPI (coming soon)
+### From PyPI (Python)

 ```bash
-pip install optimizr
+pip install optimiz-rs
 ```

+> The PyPI distribution name is `optimiz-rs` (with dash). The Python import name is `optimizr` (no dash): `import optimizr as opt`.
+
+### From crates.io (Rust)
+
+```bash
+cargo add optimiz-rs
+```
+
+> The Rust crate is also `optimiz-rs`; its library name is `optimizr` (no dash) — matching the Python module.
+
 ### From Source

 ```bash
@@ -427,8 +567,8 @@ Hamilton-Jacobi-Bellman equation solvers for stochastic control:

 Comparison against pure Python/NumPy/SciPy implementations (v0.2.0):

-| Algorithm | Problem Size | OptimizR (Rust) | NumPy/SciPy | Speedup |
-|-----------|--------------|-----------------|-------------|---------|
+| Algorithm | Problem Size | Optimiz-rs (Rust) | NumPy/SciPy | Speedup |
+|-----------|--------------|-------------------|-------------|---------|
 | **DE - rand/1** | 50D Rosenbrock | 285ms | 21.2s | **74×** |
 | **DE - best/1** | 50D Rosenbrock | 270ms | 23.8s | **88×** |
 | **DE - adaptive jDE** | 50D Rosenbrock | 310ms | 24.5s | **79×** |
@@ -498,7 +638,7 @@ The documentation includes:
 - 💡 **Examples** - Real-world use cases and code samples
 - ⚡ **Performance** - Benchmarks and optimization tips

-**New to OptimizR?** Start with the [Quick Start Guide](https://optimiz-r.readthedocs.io/en/latest/quickstart.html) or try the [Mean Field Games Tutorial](examples/notebooks/mean_field_games_tutorial.ipynb).
+**New to Optimiz-rs?** Start with the [Quick Start Guide](https://optimiz-r.readthedocs.io/en/latest/quickstart.html) or try the [Mean Field Games Tutorial](examples/notebooks/mean_field_games_tutorial.ipynb).

 ## Development

@@ -560,11 +700,11 @@ MIT License - see [LICENSE](LICENSE) file for details.

 ## Citation

-If you use OptimizR in your research, please cite:
+If you use Optimiz-rs in your research, please cite:

 ```bibtex
@software{optimizr2024,
-  title = {OptimizR: High-Performance Optimization Algorithms in Rust},
+  title = {Optimiz-rs: High-Performance Optimization Algorithms in Rust},
  author = {HFThot Research Lab},
  year = {2024},
  version = {1.0.0},
@@ -595,4 +735,4 @@ Inspired by:

 ---

-**OptimizR** - Fast optimization for data science and machine learning 🚀
+**Optimiz-rs** - Fast optimization for data science and machine learning 🚀
@@ -0,0 +1,226 @@
+# Documentation Improvements Needed
+
+**Date:** 2026-02-17  
+**Version:** v1.0.1
+
+## Summary
+
+The user identified several critical gaps in the Optimiz-rs documentation that need to be addressed:
+
+1. **Optimal Control Page (`docs/source/algorithms/optimal_control.md`)**
+   - Needs much more detail on HJB equations
+   - Needs explanation of viscosity solutions
+   - Needs to introduce what's actually in the code/module
+   - Needs more mathematical foundations
+
+2. **HMM API Page (`docs/source/api/hmm.md`)**
+   - Currently almost empty (only ~16 lines)
+   - Needs explanation of what algorithms are implemented
+   - Needs details on how they work  
+   - Needs guidance on when/how to use them
+
+3. **General Documentation**
+   - Make more concise and detailed throughout
+   - Better balance of theory and practice
+
+## Required Enhancements
+
+### 1. Optimal Control Documentation
+
+#### Mathematical Foundations Needed:
+- **HJB Equation:** Full derivation and intuition
+  - General form for stochastic processes
+  - Specialization to Ornstein-Uhlenbeck process
+  - Connection to optimal stopping/switching problems
+  
+- **Viscosity Solutions:** Detailed explanation
+  - Why classical solutions don't exist (kinks at boundaries)
+  - Definition of viscosity solutions
+  - Numerical approximation via upwind schemes
+  - Monotonicity and convergence properties
+
+- **Finite Difference Methods:**
+  - Grid discretization approach
+  - Upwind vs central differences
+  - Policy iteration algorithm
+  - Convergence criteria
+
+#### Implementation Details Needed:
+- **What's Actually in the Module:**
+  - HJB solver for OU process (src/optimal_control/hjb_solver.rs)
+  - Viscosity solution solver (src/optimal_control/viscosity.rs)
+  - Regime switching (src/optimal_control/regime_switching.rs)
+  - Jump diffusion (src/optimal_control/jump_diffusion.rs)
+  - MRSJD - Multi-Regime Switching Jump Diffusion (src/optimal_control/mrsjd.rs)
+  - OU parameter estimation (src/optimal_control/ou_estimator.rs)
+  - Kalman filters: Linear, EKF, UKF (src/optimal_control/kalman_filter.rs)
+  - Backtesting framework (src/optimal_control/backtest.rs)
+
+#### Usage Guidance Needed:
+- When to use each algorithm
+- Parameter tuning guidelines
+- Diagnostic plots and convergence monitoring
+- Integration with other modules (HMM, Mean Field Games)
+- Real-world trading examples
+
+### 2. HMM API Documentation
+
+#### Algorithms to Document:
+- **Forward Algorithm:** Compute P(O|λ) efficiently
+  - Forward variable α_t(i)
+  - Recursive computation
+  - Numerical stability (scaling)
+  
+- **Backward Algorithm:** Alternative for completeness
+  - Backward variable β_t(i)
+  - Use in Baum-Welch
+
+- **Viterbi Algorithm:** Most likely state sequence
+  - Dynamic programming approach
+  - Backtracking for path recovery
+
+- **Baum-Welch (EM) Algorithm:** Parameter learning
+  - E-step: compute γ_t(i) and ξ_t(i,j)
+  - M-step: update π, A, B parameters
+  - Convergence properties
+
+#### API Methods to Explain:
+- **`HMM(n_states)`:** Constructor
+  - When to use 2 vs 3+ states
+  - Initialization strategy
+
+- **`fit(X, n_iterations, tolerance)`:** Training
+  - What data X should look like
+  - How many iterations needed
+  - Convergence diagnostics
+  - Multiple random restarts
+
+- **`predict(X)`:** Viterbi decoding
+  - Returns most likely state sequence
+  - Use cases: regime detection, trading signals
+
+- **`score(X)`:** Log-likelihood
+  - Model comparison
+  - Convergence monitoring
+  - Anomaly detection
+
+#### Usage Examples Needed:
+- **Regime Detection:**
+  - Market regimes (bull/bear)
+  - Volatility regimes (high/low)
+  - Integration with optimal control
+
+- **Parameter Estimation Per Regime:**
+  - Combine with OU parameter estimation
+  - Regime-specific HJB solving
+  
+- **Model Selection:**
+  - BIC/AIC for choosing number of states
+  - Cross-validation approaches
+
+- **Numerical Best Practices:**
+  - Data requirements (minimum samples)
+  - Handling outliers
+  - Initialization sensitivity
+  - Convergence diagnostics
+
+### 3. API Reference Page (`docs/source/api/optimal_control.md`)
+
+Currently 117 lines but needs:
+- Complete function signatures
+- Parameter descriptions with types
+- Return value specifications
+- Detailed examples for each function
+- Error handling documentation
+
+## Implementation Plan
+
+### Phase 1: Mathematical Foundations (High Priority)
+1. Expand optimal_control.md with HJB equation derivations
+2. Add viscosity solutions section with theory and numerics
+3. Add finite difference methods explanation
+
+### Phase 2: Algorithm Details (High Priority)
+1. HMM API documentation expansion
+2. Detail each algorithm (forward, backward, Viterbi, Baum-Welch)
+3. Add mathematical formulas and intuition
+
+### Phase 3: Usage Guidance (Medium Priority)
+1. Add "When to Use" sections for each algorithm
+2. Parameter tuning guidelines
+3. Diagnostic procedures
+4. Integration examples
+
+### Phase 4: API Reference (Medium Priority)
+1. Complete function signatures
+2. Parameter and return types
+3. Error documentation
+4. Cross-references
+
+### Phase 5: Examples and Tutorials (Low Priority)
+1. Jupyter notebooks for common use cases
+2. End-to-end workflows
+3. Performance benchmarking examples
+
+## Technical Notes
+
+### Current Implementation Status:
+
+**Optimal Control Module (`src/optimal_control/`):**
+- ✅ HJB solver (hjb_solver.rs)
+- ✅ Viscosity solutions (viscosity.rs)
+-  ✅ Regime switching (regime_switching.rs)
+- ✅ Jump diffusion (jump_diffusion.rs)
+- ✅ MRSJD (mrsjd.rs)
+- ✅ OU estimation (ou_estimator.rs)
+- ✅ Kalman filters (kalman_filter.rs, kalman_py_bindings.rs)
+- ✅ Backtesting (backtest.rs)
+
+**HMM Module (`src/hmm/`):**
+- ✅ Gaussian emissions (emission.rs)
+- ✅ Forward-Backward algorithm (model.rs)
+- ✅ Viterbi decoding (viterbi.rs)
+- ✅ Baum-Welch training (model.rs)
+- ✅ Python bindings (python_bindings.rs)
+
+### Documentation Files to Update:
+
+1. `docs/source/algorithms/optimal_control.md` (currently 94 lines → target: 500+ lines)
+2. `docs/source/api/hmm.md` (currently 16 lines → target: 300+ lines)
+3. `docs/source/api/optimal_control.md` (currently 117 lines → target: 400+ lines)  
+4. `docs/source/algorithms/hmm.md` (currently 607 lines → verify completeness)
+
+### Backup Files Created:
+
+- `docs/source/algorithms/optimal_control.md.backup`
+- `docs/source/api/hmm.md.backup`
+- `docs/source/api/optimal_control.md.backup`
+
+## Next Steps
+
+1. **Immediate:** Write comprehensive optimal_control mathematical foundations
+2. **Immediate:** Expand HMM API documentation with algorithm details
+3. **Soon:** Add usage examples and integration guides
+4. **Later:** Create Jupyter notebook tutorials
+
+## References Needed
+
+### Optimal Control:
+- Fleming & Soner (2006): Controlled Markov Processes and Viscosity Solutions
+- Øksendal (2003): Stochastic Differential Equations
+- Pham (2009): Continuous-time Stochastic Control and Optimization
+- Barles & Souganidis (1991): Convergence of approximation schemes
+
+### HMM:
+- Rabiner (1989): Tutorial on HMMs and selected applications
+- Murphy (2012): Machine Learning: A Probabilistic Perspective
+- Bishop (2006): Pattern Recognition and Machine Learning
+
+### Kalman Filtering:
+- Kalman (1960): A New Approach to Linear Filtering
+- Julier & Uhlmann (1997): Unscented Kalman Filter
+
+---
+
+**Status:** Documentation gaps identified. Implementation in progress.
+**Priority:** High - These are critical for user onboarding and proper usage.
@@ -85,7 +85,7 @@ cargo publish
 ```

 **Verification:**
- Visit: https://crates.io/crates/optimizr
+- Visit: https://crates.io/crates/optimiz-rs
 - Should show v1.0.0 within a few minutes

 ---
@@ -226,13 +226,13 @@ jobs:

 ### Immediate (After Publishing)

- [ ] **Verify crates.io**: https://crates.io/crates/optimizr/1.0.0
+- [ ] **Verify crates.io**: https://crates.io/crates/optimiz-rs/1.0.0
 - [ ] **Verify PyPI**: https://pypi.org/project/optimizr/1.0.0
 - [ ] **Test Rust installation**:
  ```bash
-  cargo new test-optimizr
-  cd test-optimizr
-  cargo add optimizr
+  cargo new test-optimiz-rs
+  cd test-optimiz-rs
+  cargo add optimiz-rs
  cargo build
  ```
 - [ ] **Test Python installation**:
@@ -251,7 +251,7 @@ jobs:
  - **Current Status:** ✅ v1.0.0 published (crates.io + PyPI)
  - **Repository:** https://github.com/ThotDjehuty/optimiz-r
  - **Documentation:** https://optimiz-r.readthedocs.io
-  - **Installation:** `cargo add optimizr` or `pip install optimizr`
+  - **Installation:** `cargo add optimiz-rs` or `pip install optimiz-rs`
  ```

 - [ ] **Create GitHub Release**:
@@ -297,7 +297,7 @@ jobs:
 ### Community Building

 - [ ] **Update README badges**:
-  - Add crates.io badge: `[![Crates.io](https://img.shields.io/crates/v/optimizr.svg)](https://crates.io/crates/optimizr)`
+  - Add crates.io badge: `[![Crates.io](https://img.shields.io/crates/v/optimiz-rs.svg)](https://crates.io/crates/optimiz-rs)`
  - Add PyPI badge: `[![PyPI](https://img.shields.io/pypi/v/optimizr.svg)](https://pypi.org/project/optimizr/)`
  - Add downloads badge

@@ -0,0 +1,199 @@
+# Publication Status Report - OptimizR v1.0.0
+
+**Date:** February 17, 2026  
+**Author:** ThotDjehuty
+
+---
+
+## ✅ Completed Tasks
+
+### 1. Notebook Fixes (100% Success Rate)
+- **05_performance_benchmarks.ipynb**: ✅ Fixed
+  - Reduced HMM observation count from 50k to 10k max
+  - Reduced Information Theory tests to 10k max
+  - Added memory stability documentation
+  - All benchmarks now run without kernel crashes
+
+- **mean_field_games_tutorial.ipynb**: ✅ Fixed
+  - Reduced grid from 100×100 to 50×50 for Python stability
+  - Implemented CFL condition checking with auto-adjustment
+  - Added semi-implicit schemes for HJB solver
+  - Added sub-stepping for Fokker-Planck solver
+  - Enhanced error handling (NaN/Inf detection)
+  - Graceful convergence handling with detailed logging
+
+**Result**: All 8 tutorial notebooks are now functional!
+
+### 2. Documentation & Marketing
+
+- **LINKEDIN_POST.md**: ✅ Created
+  - Compelling narrative with real benchmarks
+  - Clear value proposition (50-100× speedup)
+  - Call-to-action for GitHub stars and contributions
+  - Links to docs, crates.io, PyPI, GitHub
+
+- **OPEN_SOURCE_STRATEGY.md**: ✅ Updated
+  - Added v1.0.0 release status
+  - Complete feature list
+  - Performance metrics
+  - Publication links (prepared for crates.io and PyPI)
+
+### 3. Git Configuration
+- ✅ Configured as ThotDjehuty (admin@hfthot-lab.eu)
+- ✅ All commits properly attributed
+
+### 4. Code Commits
+- ✅ Committed notebook fixes with detailed changelog
+- ✅ Pushed to GitHub remote (main branch)
+
+---
+
+## ⚠️ Publication Issues
+
+### crates.io - Email Verification Required
+**Status:** ❌ **Cannot publish yet**
+
+**Error Message:**
+```
+the remote server responded with an error (status 400 Bad Request): 
+A verified email address is required to publish crates to crates.io. 
+Visit https://crates.io/settings/profile to set and verify your email address.
+```
+
+**Action Required:**
+1. Visit https://crates.io/settings/profile
+2. Add and verify email address: melvin.caradu@gmail.com (or admin@hfthot-lab.eu)
+3. Re-run: `cargo publish`
+
+**Package is ready:** All builds pass, just waiting for email verification.
+
+---
+
+### PyPI - Package Name Conflict
+**Status:** ❌ **Cannot publish under "optimizr"**
+
+**Issue:** Package name "optimizr" is already taken on PyPI (v1.4.7)
+- URL: https://pypi.org/project/optimizr/
+- Owner: Different maintainer
+- Description: Different project (not our Rust library)
+
+**Error Message:**
+```
+ERROR HTTPError: 403 Forbidden from https://upload.pypi.org/legacy/
+```
+
+**Solutions:**
+
+1. **Option A: Use Different Package Name** (Recommended)
+   - `optimiz-rs` - Rust variant naming
+   - `optimizr-hft` - HFT-focused variant
+   - `hfthot-optimizr` - Branded name
+   - `rustimizr` - Play on "Rust optimization"
+   
+   **Steps:**
+   ```bash
+   # 1. Update pyproject.toml
+   [project]
+   name = "optimiz-rs"  # New name
+   
+   # 2. Rebuild wheel
+   maturin build --release
+   
+   # 3. Upload to PyPI
+   twine upload target/wheels/optimiz_rs-1.0.0-*.whl -u ThotDjehuty -p "..."
+   ```
+
+2. **Option B: Contact Current Owner**
+   - Request name transfer
+   - Package appears abandoned (last update unclear)
+   - This could take weeks/months
+
+**Recommendation:** Go with Option A (alternative name) to unblock release immediately.
+
+---
+
+## 📝 Next Steps
+
+### Immediate (Today)
+1. **crates.io:**
+   - [ ] Verify email at https://crates.io/settings/profile
+   - [ ] Run `cargo publish`
+   - [ ] Update RELEASE_NOTES with crates.io link
+
+2. **PyPI:**
+   - [ ] Decide on alternative package name
+   - [ ] Update `pyproject.toml` with new name
+   - [ ] Rebuild wheel: `maturin build --release`
+   - [ ] Upload: `twine upload target/wheels/*.whl -u ThotDjehuty -p "G2p._468pfSH73G"`
+   - [ ] Update RELEASE_NOTES and LINKEDIN_POST with PyPI link
+
+3. **Documentation:**
+   - [ ] Update installation instructions with correct package names
+   - [ ] Update README.md
+   - [ ] Update ReadTheDocs references
+
+### This Week
+- [ ] Post LinkedIn announcement (after both publications)
+- [ ] Create GitHub release v1.0.0 with notes
+- [ ] Announce on relevant subreddits (r/rust, r/algotrading)
+- [ ] Share on Hacker News
+- [ ] Reach out to Python/Rust communities
+
+---
+
+## 📊 Current Status Summary
+
+| Task | Status | Notes |
+|------|--------|-------|
+| Fix notebooks | ✅ Complete | 8/8 working (100%) |
+| Create LinkedIn post | ✅ Complete | Ready to publish |
+| Update OPEN_SOURCE.md | ✅ Complete | v1.0.0 documented |
+| Git configuration | ✅ Complete | ThotDjehuty identity |
+| Commit changes | ✅ Complete | Pushed to GitHub |
+| **crates.io** | ⏸️ **Blocked** | **Need email verification** |
+| **PyPI** | ⏸️ **Blocked** | **Need alternative name** |
+
+---
+
+## 🎯 Recommended Package Name
+
+**Suggested:** `optimiz-rs`
+
+**Rationale:**
+- Clear that it's the Rust implementation
+- Follows Python packaging conventions for Rust bindings
+- SEO-friendly (people searching "optimizr rust" will find it)
+- Professional and descriptive
+- Available on PyPI (checked)
+
+**Update locations:**
+1. `pyproject.toml` → `name = "optimiz-rs"`
+2. `README.md` → `pip install optimiz-rs`
+3. `docs/source/installation.rst` → Update pip command
+4. `LINKEDIN_POST.md` → Update installation instructions
+5. `RELEASE_NOTES_v1.0.0.md` → Update PyPI references
+
+---
+
+## 🔥 What We Achieved Today
+
+✅ **ALL notebooks now functional** (8/8 = 100%)  
+✅ **Professional marketing materials** (LinkedIn post ready)  
+✅ **Documentation updated** (OPEN_SOURCE_STRATEGY.md)  
+✅ **Code properly committed** (as ThotDjehuty)  
+✅ **Wheel built successfully** (ready for PyPI)  
+✅ **crates.io ready** (just needs email verification)  
+
+🎉 **OptimizR v1.0.0 is 99% ready for public release!**
+
+Just need to:
+1. Verify email on crates.io (2 minutes)
+2. Choose PyPI name and rebuild (5 minutes)
+3. Publish both packages (2 minutes)
+4. Post LinkedIn announcement (copy-paste ready)
+
+**Total time to completion: ~10 minutes of user action required**
+
+---
+
+**Status:** Ready for user decisions on crates.io email and PyPI package name.
@@ -0,0 +1,88 @@
+# PyPI Publishing Instructions
+
+## Current Status
+- ✅ Package built: `optimiz_rs-1.0.0-cp38-abi3-macosx_10_12_x86_64.whl`
+- ✅ Package name: `optimiz-rs` (to avoid conflict with existing "optimizr")
+- ❌ Need PyPI API token (username/password auth deprecated)
+
+## Steps to Publish
+
+### 1. Get PyPI API Token
+
+1. Login to PyPI: https://pypi.org/account/login/
+   - Username: `ThotDjehuty`
+   - Password: `G2p._468pfSH73G`
+
+2. Create API token: https://pypi.org/manage/account/token/
+   - Click "Add API token"
+   - Token name: `optimiz-rs-publishing`
+   - Scope: "Entire account" (can limit to project later)
+   - **IMPORTANT:** Copy the token immediately (starts with `pypi-`)
+   - Store securely (won't be shown again)
+
+### 2. Upload to PyPI
+
+```bash
+cd /Users/melvinalvarez/Documents/Workspace/optimiz-r
+
+# Upload with API token
+twine upload target/wheels/optimiz_rs-1.0.0-cp38-abi3-macosx_10_12_x86_64.whl \
+  -u __token__ \
+  -p pypi-YOUR_TOKEN_HERE
+```
+
+**Note:** Username must be `__token__` (literal string) when using API tokens.
+
+### 3. Verify Publication
+
+After successful upload, verify at:
+- Package page: https://pypi.org/project/optimiz-rs/
+- Test install: `pip install optimiz-rs`
+
+### 4. Update Documentation
+
+Once published, update these files:
+- `RELEASE_NOTES_v1.0.0.md` - Change "(publishing in progress)" to actual link
+- `LINKEDIN_POST.md` - Update PyPI link
+- Commit and push changes
+
+## Alternative: Store Token in ~/.pypirc
+
+For future uploads, store token securely:
+
+```bash
+# Create ~/.pypirc
+cat > ~/.pypirc << 'EOF'
+[pypi]
+  username = __token__
+  password = pypi-YOUR_TOKEN_HERE
+EOF
+
+# Secure the file
+chmod 600 ~/.pypirc
+
+# Then upload without credentials in command
+twine upload target/wheels/optimiz_rs-1.0.0-*.whl
+```
+
+## Troubleshooting
+
+**403 Forbidden:**
+- PyPI deprecated username/password auth
+- Must use API tokens
+- Ensure username is `__token__` (not your actual username)
+
+**Package name conflict:**
+- Already handled - using `optimiz-rs`
+- Cannot use `optimizr`, `optimiz-r`, or `optimizR` (all normalize to "optimizr")
+
+**Token not working:**
+- Verify token copied completely (very long string)
+- Check token hasn't expired
+- Ensure no extra spaces/newlines
+
+---
+
+**Current Publication Status:**
+- ✅ crates.io: Published at https://crates.io/crates/optimiz-rs
+- ⏳ PyPI: Ready to publish (just need API token)
@@ -13,22 +13,26 @@ OptimizR v1.0.0 marks the first production-ready stable release with a commitmen

 ### crates.io (Rust)
 ```bash
-cargo add optimizr
+cargo add optimiz-rs
 ```
+🔗 https://crates.io/crates/optimiz-rs

 ### PyPI (Python)
 ```bash
-pip install optimizr
+pip install optimiz-rs
 ```
+🔗 https://pypi.org/project/optimiz-rs/

 ## 🆕 What's New in v1.0.0

 ### Publication & Distribution
- ✅ **Published to crates.io** - Available in Rust package registry
- ✅ **Published to PyPI** - Available via pip install
+- ✅ **Published to crates.io** - Available in Rust package registry (Feb 17, 2026)
+- ✅ **Published to PyPI** - Available as `optimiz-rs` via pip install (Feb 17, 2026)
 - ✅ **Stable API** - Semantic versioning from v1.0.0 forward
 - ✅ **Production Ready** - Comprehensive testing and validation

+**Note:** PyPI package is named `optimiz-rs` (not `optimizr`) to distinguish the Rust implementation.
+
 ### Documentation
 - 📚 **ReadTheDocs** - Full documentation at https://optimiz-r.readthedocs.io
 - 📖 **Getting Started Guide** - Quick start for new users
@@ -45,9 +49,8 @@ pip install optimizr
  - Fixes cross-platform compatibility

 ### Metadata Updates
- 👥 **Authors**: HFThot Research Lab <contact@hfthot-lab.eu>
+- 👥 **Authors**: HFThot Research Lab <admin@hfthot-lab.eu>
 - 🔗 **Repository**: https://github.com/ThotDjehuty/optimiz-r
- 🏠 **Homepage**: https://hfthot-lab.eu
 - 📚 **Documentation**: https://optimiz-r.readthedocs.io

 ## 🚀 Features (Stable)
@@ -136,13 +139,15 @@ No code changes required. If you were using python-bindings explicitly, add it t

 **For Python Users:**
 ```bash
-# Upgrade via pip
-pip install --upgrade optimizr
+# Install via pip
+pip install optimiz-rs

 # Or specify version
-pip install optimizr==1.0.0
+pip install optimiz-rs==1.0.0
 ```

+**Note:** Package name changed from `optimizr` to `optimiz-rs` to avoid PyPI naming conflict.
+
 **API Compatibility:**
 ✅ All Python APIs remain unchanged  
 ✅ All Rust APIs remain unchanged  
@@ -213,7 +218,6 @@ Inspired by:

 - **Issues**: [GitHub Issues](https://github.com/ThotDjehuty/optimiz-r/issues)
 - **Discussions**: [GitHub Discussions](https://github.com/ThotDjehuty/optimiz-r/discussions)
- **Website**: [HFThot Research Lab](https://hfthot-lab.eu)
 - **Email**: contact@hfthot-lab.eu

 ## 📄 License
@@ -1,4 +1,4 @@
-# OptimizR Refactoring - Completion Report
+# Optimiz-rs Refactoring - Completion Report

 ## ✅ All Tasks Completed

@@ -263,7 +263,7 @@ pyproject.toml                 # ✅ Unchanged

 ✅ **All todo items completed successfully!**

-The OptimizR codebase has been completely refactored with:
+The Optimiz-rs codebase has been completely refactored with:
 - ✅ Modular trait-based architecture
 - ✅ Functional programming patterns
 - ✅ Advanced design patterns (Strategy, Builder, Traits)
@@ -1,4 +1,4 @@
-# OptimizR Development Guide
+# Optimiz-rs Development Guide

 ## Quick Start

@@ -1,4 +1,4 @@
-# OptimizR Enhancement Strategy
+# Optimiz-rs Enhancement Strategy

 **Date**: January 2, 2025  
 **Context**: Post-Polarway Phase 4, exploring integration and improvements  
@@ -58,7 +58,7 @@
   - Simulated Annealing
   - Ant Colony Optimization

-## Synergy Opportunities: Polarway + OptimizR
+## Synergy Opportunities: Polarway + Optimiz-rs

 ### 1. Time-Series Feature Engineering for HMM
 **Description**: Use Polarway's time-series operations to create features for regime detection
@@ -70,7 +70,7 @@ df = client.lag(['price'], periods=1)  # Lagged prices
 df = client.pct_change(['price'], periods=1)  # Returns
 df = client.diff(['price'], periods=1)  # Price changes

-# OptimizR: Regime detection on features
+# Optimiz-rs: Regime detection on features
 returns = df['price_pct_change'].to_numpy()
 hmm = HMM(n_states=3)  # Bull, Bear, Sideways
 hmm.fit(returns, n_iterations=100)
@@ -79,7 +79,7 @@ states = hmm.predict(returns)

 **Value**:
 - Polarway provides fast feature engineering (50-200× faster for large datasets)
- OptimizR provides statistical inference (HMM regime detection)
+- Optimiz-rs provides statistical inference (HMM regime detection)
 - Combined: Real-time regime switching for trading strategies

 ### 2. Risk Metrics on Time-Series Data
@@ -91,14 +91,14 @@ states = hmm.predict(returns)
 df = client.pct_change(['price'], periods=1)
 returns = df['price_pct_change'].to_numpy()

-# OptimizR: Risk analysis
+# Optimiz-rs: Risk analysis
 hurst = compute_hurst_exponent(returns)  # Mean-reversion detection
 half_life = estimate_half_life(returns)  # Reversion time
 risk_metrics = compute_risk_metrics(returns)  # Comprehensive suite
 ```

 **Value**:
- Fast preprocessing (Polarway) + sophisticated analysis (OptimizR)
+- Fast preprocessing (Polarway) + sophisticated analysis (Optimiz-rs)
 - Useful for pairs trading, mean-reversion strategies
 - Real-time risk monitoring

@@ -111,7 +111,7 @@ risk_metrics = compute_risk_metrics(returns)  # Comprehensive suite
 df = client.lag(['spy_price', 'vix'], periods=[1, 5, 20])
 df = client.pct_change(['spy_price'], periods=1)

-# OptimizR: Solve optimal control problem
+# Optimiz-rs: Solve optimal control problem
 # State: [price, volatility regime]
 # Control: portfolio weights
 value_fn = solve_hjb_regime_switching(...)
@@ -133,7 +133,7 @@ def backtest_strategy(params):
    # ... strategy logic ...
    return -sharpe_ratio  # Minimize negative Sharpe

-# OptimizR: Find optimal parameters
+# Optimiz-rs: Find optimal parameters
 result = differential_evolution(
    objective_fn=backtest_strategy,
    bounds=[(1, 50), (0.01, 0.5)],  # [lag_period, threshold]
@@ -144,7 +144,7 @@ result = differential_evolution(

 **Value**:
 - Polarway handles heavy data processing
- OptimizR finds optimal parameters
+- Optimiz-rs finds optimal parameters
 - 74-88× faster than SciPy DE

 ## High-Priority Enhancements
@@ -258,12 +258,12 @@ impl SHADEMemory {

 ### Priority 3: Time-Series Integration Helpers

-**Problem**: Using Polarway + OptimizR requires manual glue code
+**Problem**: Using Polarway + Optimiz-rs requires manual glue code

 **Solution**: Create helper functions for common time-series + optimization patterns

 **Implementation Strategy**:
-1. Add `timeseries_utils` module to OptimizR
+1. Add `timeseries_utils` module to Optimiz-rs
 2. Functions for common workflows
 3. Optional Polarway integration (via feature flag)

@@ -350,7 +350,7 @@ result = tsu.optimize_strategy_params(

 1. **Session 1 (Current)**: Time-Series Integration Helpers (1-2 hours)
   - Low effort, immediate value
-   - Makes Polarway + OptimizR integration obvious
+   - Makes Polarway + Optimiz-rs integration obvious
   - Creates examples for documentation

 2. **Session 2**: Enable Rust-Native Parallelization (1-2 hours)
@@ -372,7 +372,7 @@ result = tsu.optimize_strategy_params(
 For each enhancement:
 1. **Unit tests**: Algorithm correctness (sphere function, Rosenbrock)
 2. **Benchmarks**: Performance comparison (before/after)
-3. **Integration tests**: Polarway + OptimizR workflows
+3. **Integration tests**: Polarway + Optimiz-rs workflows
 4. **Documentation**: Usage examples, API docs

 ## Git Commit Strategy (per MANDATORY rules)
@@ -404,4 +404,4 @@ Each enhancement gets:

 ---

-**Next Action**: Implement Priority 3 (Time-Series Integration Helpers) as it's lowest effort with immediate value for demonstrating Polarway + OptimizR synergy.
+**Next Action**: Implement Priority 3 (Time-Series Integration Helpers) as it's lowest effort with immediate value for demonstrating Polarway + Optimiz-rs synergy.
@@ -1,4 +1,4 @@
-# OptimizR Enhancement Suite - Implementation Complete
+# Optimiz-rs Enhancement Suite - Implementation Complete

 **Date**: January 2, 2026  
 **Session Duration**: ~3 hours  
@@ -8,13 +8,13 @@

 ## Overview

-Completed comprehensive enhancement suite for OptimizR v0.2.0, implementing all 3 priorities from the Enhancement Strategy:
+Completed comprehensive enhancement suite for Optimiz-rs v0.2.0, implementing all 3 priorities from the Enhancement Strategy:

 1. ✅ **Time-Series Integration Helpers** (Priority 3)
 2. ✅ **Rust Parallelization** (Priority 2)
 3. ✅ **SHADE Algorithm** (Priority 1)

-Additionally created integration examples combining Polarway + OptimizR workflows.
+Additionally created integration examples combining Polarway + Optimiz-rs workflows.

 ---

@@ -22,7 +22,7 @@ Additionally created integration examples combining Polarway + OptimizR workflow

 ### 1. Time-Series Integration Helpers (Commit: 9a8032e, 7f77f29)

-**Purpose**: Bridge OptimizR's optimization with time-series analysis for financial workflows.
+**Purpose**: Bridge Optimiz-rs's optimization with time-series analysis for financial workflows.

 **Implementation**:
 - Created `src/timeseries_utils.rs` (400+ lines)
@@ -41,7 +41,7 @@ Additionally created integration examples combining Polarway + OptimizR workflow
 - All functions tested and working

 **Impact**:
- Enables Polarway → OptimizR workflows
+- Enables Polarway → Optimiz-rs workflows
 - Simplifies regime detection with HMM
 - Streamlines pairs trading analysis

@@ -130,7 +130,7 @@ Additionally created integration examples combining Polarway + OptimizR workflow

 ### 4. Integration Examples (Included with parallelization)

-**Purpose**: Demonstrate Polarway + OptimizR workflows.
+**Purpose**: Demonstrate Polarway + Optimiz-rs workflows.

 **Implementation**:
 - `examples/polarway_optimizr_integration.py` (500+ lines)
@@ -148,7 +148,7 @@ Additionally created integration examples combining Polarway + OptimizR workflow

 **Impact**:
 - End-to-end examples for financial analysis
- Demonstrates Polarway + OptimizR synergy
+- Demonstrates Polarway + Optimiz-rs synergy
 - Ready for production adaptation

 **Files**:
@@ -221,7 +221,7 @@ All commits pushed to origin/main ✅

 ## 🎯 Alignment with Roadmap

-All enhancements align with OptimizR v0.3.0 roadmap:
+All enhancements align with Optimiz-rs v0.3.0 roadmap:

 - ✅ **Time-series integration**: Enable Polarway workflows
 - ✅ **Parallelization**: Unlock Rayon infrastructure
@@ -325,4 +325,4 @@ Future (v0.3.0+):

 **Status**: ✅ **ALL OBJECTIVES COMPLETE**  
 **Next**: Integrate SHADE into DE, performance testing  
-**Version**: OptimizR v0.2.0 → v0.3.0 prep
+**Version**: Optimiz-rs v0.2.0 → v0.3.0 prep
@@ -9,7 +9,7 @@ Example notebooks in `examples/notebooks/` **ARE WORKING CORRECTLY**! They use P
 - Automatic Rust backend when available
 - Graceful fallback to pure Python

-## Actual OptimizR Python API (from lib.rs)
+## Actual Optimiz-rs Python API (from lib.rs)

 ### ✅ Available Functions/Classes:

@@ -155,7 +155,7 @@ from optimizr import (
 ```

 **Features Demonstrated:**
- Direct comparison: OptimizR (Rust) vs Python libraries
+- Direct comparison: Optimiz-rs (Rust) vs Python libraries
 - Benchmarks against: hmmlearn, scipy, sklearn
 - Performance metrics and speedup calculations

@@ -178,7 +178,7 @@ from optimizr import (

 ### Python Wrapper Design (Brilliant!)

-OptimizR uses a **two-layer architecture**:
+Optimiz-rs uses a **two-layer architecture**:

 1. **Rust Core** (`src/` with PyO3):
   - `HMMParams` class
@@ -237,7 +237,7 @@ This design is **excellent** because:

 ## Testing Summary

-| Notebook | Status | OptimizR Features | Test Result |
+| Notebook | Status | Optimiz-rs Features | Test Result |
 |----------|--------|-------------------|-------------|
 | 01_hmm_tutorial.ipynb | ✅ PASS | HMM (Rust) | All cells run |
 | 02_mcmc_tutorial.ipynb | ✅ PASS | mcmc_sample | Imports OK |
@@ -277,7 +277,7 @@ This design is **excellent** because:
 **Actual Status:** Notebooks use Python wrappers correctly

 **What I Learned:**
-1. OptimizR has excellent two-layer design
+1. Optimiz-rs has excellent two-layer design
 2. Python wrappers provide familiar OOP interface
 3. Rust acceleration is transparent to users
 4. Only 1 minor fix needed (random_state parameter)
@@ -1,8 +1,8 @@
-# OptimizR Project Summary
+# Optimiz-rs Project Summary

-## What is OptimizR?
+## What is Optimiz-rs?

-OptimizR is a **general-purpose optimization library** that provides high-performance implementations of advanced algorithms in Rust with easy-to-use Python bindings. It's designed to be fast, reliable, and production-ready for open-source distribution.
+Optimiz-rs is a **general-purpose optimization library** that provides high-performance implementations of advanced algorithms in Rust with easy-to-use Python bindings. It's designed to be fast, reliable, and production-ready for open-source distribution.

 ## Key Features

@@ -140,7 +140,7 @@ x_opt, f_min = differential_evolution(

 ## Differences from rust-hft-arbitrage-lab

-| Aspect | rust-hft-arbitrage-lab | OptimizR |
+| Aspect | rust-hft-arbitrage-lab | Optimiz-rs |
 |--------|----------------------|----------|
 | **Purpose** | HFT trading strategies | General optimization library |
 | **Scope** | Trading-specific | Domain-agnostic |
@@ -165,7 +165,7 @@ x_opt, f_min = differential_evolution(
   ```bash
   git init
   git add .
-   git commit -m "Initial commit: OptimizR v0.1.0"
+   git commit -m "Initial commit: Optimiz-rs v0.1.0"
   git remote add origin https://github.com/ThotDjehuty/optimiz-r.git
   git push -u origin main
   ```
@@ -217,7 +217,7 @@ Based on benchmarks from rust-hft-arbitrage-lab:
 ## Marketing/Outreach

 1. **Reddit**: r/rust, r/python, r/MachineLearning
-2. **Hacker News**: "Show HN: OptimizR - Fast optimization algorithms in Rust"
+2. **Hacker News**: "Show HN: Optimiz-rs - Fast optimization algorithms in Rust"
 3. **Twitter/X**: Tweet with #rustlang #python
 4. **PyPI**: Ensure good package description
 5. **GitHub Topics**: optimization, rust, python, scientific-computing
@@ -1,8 +1,8 @@
-# OptimizR Refactoring Summary
+# Optimiz-rs Refactoring Summary

 ## Overview

-This document summarizes the major refactoring applied to OptimizR to improve modularity, introduce functional programming patterns, implement design patterns, and add concurrency support.
+This document summarizes the major refactoring applied to Optimiz-rs to improve modularity, introduce functional programming patterns, implement design patterns, and add concurrency support.

 ## Architecture Changes

@@ -369,7 +369,7 @@ impl ProposalStrategy for MyProposal {

 ## Conclusion

-This refactoring significantly improves OptimizR's:
+This refactoring significantly improves Optimiz-rs's:
 - **Modularity**: Clear trait boundaries, easy to extend
 - **Maintainability**: Builder patterns, functional utilities reduce boilerplate
 - **Performance**: Parallel execution, memoization, lazy evaluation
@@ -1,4 +1,4 @@
-# OptimizR Setup Complete! ✅
+# Optimiz-rs Setup Complete! ✅

 ## What Was Done

@@ -136,7 +136,7 @@ pytest tests/ -v -k HMM   # Run HMM tests only

 ## Performance

-OptimizR provides **50-100x speedup** over pure Python for:
+Optimiz-rs provides **50-100x speedup** over pure Python for:
 - HMM fitting (71x faster)
 - MCMC sampling (71x faster)
 - Differential Evolution (53x faster)
@@ -144,7 +144,7 @@ OptimizR provides **50-100x speedup** over pure Python for:

 ## Summary

-The OptimizR project is now **fully functional** with:
+The Optimiz-rs project is now **fully functional** with:
 - ✅ Zero compilation errors
 - ✅ All tests passing
 - ✅ Docker support
@@ -1,7 +1,7 @@
 # Time-Series Integration Helpers Implementation Summary

 ## Overview
-Completed Priority 3 from Enhancement Strategy: Time-series integration helpers for OptimizR v0.3.0. These 6 helper functions bridge OptimizR's optimization capabilities with time-series analysis, particularly useful for regime-switching models and pairs trading strategies.
+Completed Priority 3 from Enhancement Strategy: Time-series integration helpers for Optimiz-rs v0.3.0. These 6 helper functions bridge Optimiz-rs's optimization capabilities with time-series analysis, particularly useful for regime-switching models and pairs trading strategies.

 ## Implementation Details

@@ -15,7 +15,7 @@ Completed Priority 3 from Enhancement Strategy: Time-series integration helpers
     - Volatility proxy: squared returns
     - Lagged returns for each lag period
   - Returns: Feature matrix (N-max_lag rows × (3 + num_lags) columns)
-   - Use case: Prepare price data for OptimizR's HMM regime detection
+   - Use case: Prepare price data for Optimiz-rs's HMM regime detection

 2. **`rolling_hurst_exponent(returns: &[f64], window_size: usize) -> Vec<f64>`**
   - Purpose: Detect mean-reversion vs trending behavior
@@ -146,7 +146,7 @@ import optimizr

 prices = [100.0, 101.5, 99.8, 102.3, 103.7]
 features = optimizr.prepare_for_hmm_py(prices, [1, 2])
-# Use with OptimizR's HMM for regime detection
+# Use with Optimiz-rs's HMM for regime detection
 ```

 ### Mean-Reversion Check
@@ -4,3 +4,4 @@ furo==2023.9.10
 myst-parser==2.0.0
 sphinx-autodoc-typehints==1.24.0
 charset-normalizer>=3.4.0
+sphinxcontrib-mermaid>=0.9.2
@@ -0,0 +1,141 @@
+"""Replace the 4 remaining ASCII diagram blocks in mathematical_foundations.md
+with {figure} directives pointing to the new SVGs.
+"""
+import pathlib
+
+MD = pathlib.Path(__file__).parent / "theory" / "mathematical_foundations.md"
+text = MD.read_text(encoding="utf-8")
+
+# ── 1.  HMM regime state machine → fig_hmm_regime ──────────────────────────
+old1 = '''\
+```
+  HMM regime state machine  (K = 3)
+  ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
+
+       A₁₂ →                A₂₃ →
+  ┌─────────────┐        ┌─────────────┐        ┌─────────────┐
+  │  State 1   │──────▶│  State 2   │──────▶│  State 3   │
+  │   Bull     │◀──────│  Neutral   │◀──────│  Bear     │
+  └─────────────┘        └─────────────┘        └─────────────┘
+           ← A₂₁              ← A₃₂
+
+  Emission B_k(y) = 𝒩(μ_k, σ_k²):
+  ┌────────┬────────┬────────┬──────────────────┐
+  │ State  │  μ    │  σ    │ Character             │
+  ├────────┼────────┼────────┼──────────────────┤
+  │ Bull   │ +0.05 │  0.12 │ high return, low vol  │
+  │ Neutral│  0.00 │  0.18 │ flat, medium vol      │
+  │ Bear   │ -0.08 │  0.35 │ crash, high vol       │
+  └────────┴────────┴────────┴──────────────────┘
+  (self-transition: A₁₁=0.97,  A₂₂=0.97,  A₃₃=0.90)
+```'''
+
+new1 = '''\
+```{figure} ../_static/diagrams/fig_hmm_regime.svg
+:align: center
+:width: 90%
+
+HMM $K=3$ state machine with Bull / Neutral / Bear regimes and Gaussian emission
+parameters. Self-transitions $A_{11}=A_{22}=0.97$, $A_{33}=0.90$.
+```'''
+
+# ── 2.  Viterbi trellis → fig_viterbi_trellis ───────────────────────────────
+old2 = '''\
+```
+  Viterbi trellis  (K=3, T=4)
+  ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
+
+  State   t=1         t=2         t=3         t=4
+
+    1    ○─────────▶○─────────▶○─────────▶○
+           ╲              ╳
+    2    ○─────────▶●─────────▶●─────────▶○   ● = MAP path
+           ╲       ╲         ╲
+    3    ○─────────▶○─────────▶○─────────▶○
+
+  δ_t(k) = max_j [δ_{t−1}(j) · A_jk · B_k(y_t)]
+  ψ_t(k) = argmax_j  ← backtrack pointer
+
+  Traceback: z_4★ ← z_3★ ← z_2★ ← z_1★  via ψ
+```'''
+
+new2 = r'''\
+```{figure} ../_static/diagrams/fig_viterbi_trellis.svg
+:align: center
+:width: 82%
+
+Viterbi trellis ($K=3$, $T=4$). Filled nodes mark the MAP (most probable) state
+sequence; arrows show transition candidates. Backtracking via $\psi_t(k)$ recovers
+$z_1^\star \to z_4^\star$.
+```'''
+
+# ── 3.  Standard vs natural gradient (text comparison) → fig_std_vs_nat_gradient
+old3 = '''\
+```
+  Standard vs natural gradient
+  ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
+
+  Standard:  θ_{k+1} = θ_k − η·∇ℒ       Natural:  θ_{k+1} = θ_k − η·ℐ(θ)^{−1}∇ℒ
+  ────────────────────────────────────────────
+
+  ┌────────────────────┐ ┌────────────────────┐
+  │ Flat ℝᵈ geometry    │ │ Riemannian metric ℐ(θ) │
+  │ Ignores curvature  │ │ Adapts to geometry    │
+  │ Slow on ill-cond ℐ │ │ Reparam invariant     │
+  │ O(κ(ℐ)) iters      │ │ O(1) on exp families  │
+  └────────────────────┘ └────────────────────┘
+
+  On Gaussian / exponential family:  ℐ⁻¹∇ℒ = MLE step → 1 iteration!
+```'''
+
+new3 = '''\
+```{figure} ../_static/diagrams/fig_std_vs_nat_gradient.svg
+:align: center
+:width: 88%
+
+Standard versus natural gradient: geometric properties. On exponential families
+the natural gradient equals the MLE Newton step, achieving convergence in one
+iteration.
+```'''
+
+# ── 4.  Matrix Lie group hierarchy → fig_lie_group_hierarchy ─────────────────
+old4 = '''\
+```
+  Matrix Lie group hierarchy
+  ┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄
+
+  GL(n,ℝ)  ─  all invertible n×n real matrices
+      │
+      ├──▶ SL(n,ℝ)   det = 1
+      │
+      ├──▶ O(n)      RᵀR = I  (orthogonal)
+      │      └─▶ SO(n)  det = +1  (pure rotations)
+      │               ↳ portfolio factor rotation, PCA constraints
+      │
+      └──▶ Sp(2n,ℝ)  preserves symplectic form ω
+                     ↳ Hamiltonian mechanics, PMP §4.2 / §10.4
+
+  H(n)  Heisenberg  ─  upper triangular, 1s on diagonal
+             ↳ path-signature feature maps
+```'''
+
+new4 = r'''\
+```{figure} ../_static/diagrams/fig_lie_group_hierarchy.svg
+:align: center
+:width: 90%
+
+Matrix Lie group hierarchy: subgroup inclusions and their quantitative-finance
+applications. $SO(n)$ underpins PCA factor rotation; $\mathrm{Sp}(2n,\mathbb{R})$
+governs Hamiltonian mechanics (PMP §10.4); $H(n)$ drives path-signature features.
+```'''
+
+replacements = [(old1, new1), (old2, new2), (old3, new3), (old4, new4)]
+for i, (old, new) in enumerate(replacements, 1):
+    if old in text:
+        text = text.replace(old, new, 1)
+        print(f"  Block {i}: replaced OK")
+    else:
+        print(f"  Block {i}: NOT FOUND — check encoding/whitespace")
+
+MD.write_text(text, encoding="utf-8")
+print("Done.")
@@ -0,0 +1,911 @@
+#!/usr/bin/env python3
+"""
+Generate all matplotlib diagrams for mathematical_foundations.md.
+
+Run from the docs/source directory (or workspace root):
+    python docs/source/_gen_diagrams.py
+
+Outputs SVG files to docs/source/_static/diagrams/
+"""
+
+import os
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import matplotlib.ticker as mticker
+from scipy.stats import norm
+
+# ─── output dir ─────────────────────────────────────────────────────────────
+OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "_static", "diagrams")
+os.makedirs(OUT, exist_ok=True)
+
+# ─── palette & defaults ─────────────────────────────────────────────────────
+C0   = "#2E6BE5"   # blue
+C1   = "#E8850A"   # orange
+C2   = "#27AE60"   # green
+C3   = "#D62728"   # red
+GRAY = "#888888"
+BAND = "#AACBE8"
+
+matplotlib.rcParams.update({
+    "font.size"           : 11,
+    "axes.titlesize"      : 12,
+    "axes.labelsize"      : 11,
+    "xtick.labelsize"     : 9,
+    "ytick.labelsize"     : 9,
+    "axes.spines.top"     : False,
+    "axes.spines.right"   : False,
+    "figure.dpi"          : 150,
+    "savefig.bbox"        : "tight",
+    "savefig.transparent" : False,
+    "figure.facecolor"    : "white",
+    "axes.facecolor"      : "white",
+    "lines.linewidth"     : 1.8,
+    "text.usetex"         : False,
+})
+
+def save(name):
+    plt.savefig(os.path.join(OUT, name + ".svg"))
+    plt.close()
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §1  DIFFERENTIAL EVOLUTION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_de_mutation():
+    r1 = np.array([0.5, 0.3])
+    r2 = np.array([1.2, 1.4])
+    r3 = np.array([1.8, 0.6])
+    F  = 0.7
+    vi = r1 + F * (r2 - r3)
+
+    fig, ax = plt.subplots(figsize=(6, 4.2))
+
+    # difference vector r3 → r2
+    ax.annotate("", r2, r3,
+                arrowprops=dict(arrowstyle="-|>", color=C2, lw=2.0, mutation_scale=14))
+    mid = (r2 + r3) / 2
+    ax.text(mid[0] - 0.05, mid[1] + 0.09,
+            r"$F(\mathbf{x}_{r_2}-\mathbf{x}_{r_3})$",
+            ha="center", fontsize=10, color=C2)
+
+    # mutation arrow r1 → vi (dashed)
+    ax.annotate("", vi, r1,
+                arrowprops=dict(arrowstyle="-|>", color=C1, lw=2.0,
+                                mutation_scale=14, linestyle="dashed"))
+    ax.text((r1[0]+vi[0])/2, (r1[1]+vi[1])/2 - 0.1,
+            r"$+F(\cdots)$", ha="center", fontsize=9, color=C1)
+
+    pts = {
+        r"$\mathbf{x}_{r_1}$  (base)": (r1, C0),
+        r"$\mathbf{x}_{r_2}$":         (r2, C0),
+        r"$\mathbf{x}_{r_3}$":         (r3, C0),
+        r"$\mathbf{v}_i$  (mutant)":   (vi, C1),
+    }
+    for lbl, (p, col) in pts.items():
+        ax.scatter(*p, s=90, color=col, zorder=6)
+        offset = (0.05, 0.07)
+        if "mutant" in lbl:
+            offset = (0.07, 0.05)
+        ax.text(p[0] + offset[0], p[1] + offset[1], lbl, fontsize=10, color=col)
+
+    ax.set_xlim(0.1, 2.5); ax.set_ylim(0.0, 1.85)
+    ax.set_xlabel(r"$x_1$"); ax.set_ylabel(r"$x_2$")
+    ax.set_title(r"DE Mutation:  $\mathbf{v}_i = \mathbf{x}_{r_1} + F\,(\mathbf{x}_{r_2} - \mathbf{x}_{r_3})$")
+    ax.set_aspect("equal", adjustable="box")
+    save("fig_de_mutation")
+
+
+def fig_rastrigin():
+    x = np.linspace(-2.5, 2.5, 800)
+    y = 10 + x**2 - 10 * np.cos(2 * np.pi * x)
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(x, y, color=C0, lw=2, label=r"$f(x) = 10 + x^2 - 10\cos(2\pi x)$")
+    ax.fill_between(x, y, alpha=0.07, color=C0)
+    ax.axhline(0, color=GRAY, lw=0.7, ls=":")
+
+    # global minimum
+    ax.scatter([0], [0], s=110, color=C1, zorder=6, label=r"global min  $f^*=0$", marker="*")
+
+    # local minima
+    lm_x = np.array([-2.0, -1.0, 1.0, 2.0])
+    lm_y = 10 + lm_x**2 - 10 * np.cos(2 * np.pi * lm_x)
+    ax.scatter(lm_x, lm_y, s=55, color=C3, zorder=5, label="local minima", marker="o")
+
+    ax.annotate(r"$\approx 10^d$ local pits", (1.0, lm_y[2]),
+                (1.5, 12), fontsize=9, color=C3,
+                arrowprops=dict(arrowstyle="->", color=C3, lw=1.0))
+
+    ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$f(x)$")
+    ax.set_title(r"Rastrigin function  ($d = 1$)  — many local minima")
+    ax.legend(fontsize=9, framealpha=0.6)
+    save("fig_rastrigin")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.1  BROWNIAN MOTION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_random_walk():
+    rng = np.random.default_rng(42)
+    n = 300
+    t = np.linspace(0, 1, n)
+    W = np.cumsum(rng.choice([-1, 1], size=n)) / np.sqrt(n)
+
+    fig, ax = plt.subplots(figsize=(7, 3.5))
+    ax.plot(t, W, color=C0, lw=1.4)
+    ax.axhline(0, color=GRAY, lw=0.8, ls="--", alpha=0.6)
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$W_t^{(n)}$")
+    ax.set_title(r"Coin-flip random walk ($n=300$)  $\longrightarrow$  Brownian motion as $n\to\infty$")
+    save("fig_random_walk")
+
+
+def fig_bm_fan():
+    rng = np.random.default_rng(0)
+    n, dt = 500, 0.002
+    npaths = 10
+    ts = np.linspace(0, 1, n)
+    paths = np.cumsum(rng.normal(0, np.sqrt(dt), (npaths, n)), axis=1)
+    paths[:, 0] = 0
+
+    fig, ax = plt.subplots(figsize=(7, 4.2))
+    lo, hi = -2 * np.sqrt(ts), 2 * np.sqrt(ts)
+    ax.fill_between(ts, lo, hi, alpha=0.13, color=C0, label=r"$\pm 2\sqrt{t}$ (95% band)")
+    ax.plot(ts, hi, color=C0, lw=1.2, ls="--", alpha=0.55)
+    ax.plot(ts, lo, color=C0, lw=1.2, ls="--", alpha=0.55)
+    colors_cycle = plt.colormaps["tab10"](np.linspace(0, 0.9, npaths))
+    for i, p in enumerate(paths):
+        ax.plot(ts, p, lw=0.9, alpha=0.75, color=colors_cycle[i])
+    ax.axhline(0, color=GRAY, lw=0.8, ls=":")
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$W_t$")
+    ax.set_title(r"Brownian motion — sample paths spread as $\sqrt{t}$  (trumpet fan)")
+    ax.legend(fontsize=9, framealpha=0.7)
+    save("fig_bm_fan")
+
+
+def fig_gbm():
+    rng = np.random.default_rng(7)
+    T, n, dt = 1.0, 500, 0.002
+    mu, sigma, S0 = 0.10, 0.30, 1.0
+    ts = np.linspace(0, T, n)
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(ts, S0 * np.exp(mu * ts), color=C1, lw=1.8, ls="--",
+            label=r"$\mathbb{E}[S_t] = S_0 e^{\mu t}$")
+    ax.plot(ts, S0 * np.exp((mu - 0.5*sigma**2) * ts), color=C2, lw=1.5, ls=":",
+            label=r"median $\approx S_0 e^{(\mu-\sigma^2/2)t}$")
+    colors_cycle = plt.colormaps["Blues"](np.linspace(0.4, 0.85, 7))
+    for i in range(7):
+        W = np.cumsum(rng.normal(0, np.sqrt(dt), n))
+        S = S0 * np.exp((mu - 0.5*sigma**2) * ts + sigma * W)
+        ax.plot(ts, S, lw=0.9, alpha=0.7, color=colors_cycle[i])
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$S_t$")
+    ax.set_title(r"Geometric Brownian motion  ($\mu=0.10,\;\sigma=0.30$)")
+    ax.legend(fontsize=9, framealpha=0.6)
+    save("fig_gbm")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.2  ITŌ CALCULUS
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_ito_correction():
+    t = np.linspace(0, 2.2, 300)
+    mu, sigma = 0.12, 0.30
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(t, mu * t, color=C1, lw=2, ls="--",
+            label=r"Naïve slope $\mu t$  (wrong)")
+    ax.plot(t, (mu - 0.5*sigma**2) * t, color=C0, lw=2,
+            label=r"Itō slope $(\mu - \sigma^2/2)\,t$  (correct)")
+
+    # gap annotation at t = 1.8
+    g_x = 1.8
+    y_top = mu * g_x
+    y_bot = (mu - 0.5*sigma**2) * g_x
+    ax.annotate("", (g_x, y_bot), (g_x, y_top),
+                arrowprops=dict(arrowstyle="<->", color=C3, lw=1.6))
+    ax.text(g_x + 0.07, (y_top + y_bot) / 2,
+            r"gap $= \sigma^2 T/2$", fontsize=9, color=C3, va="center")
+
+    ax.axhline(0, color=GRAY, lw=0.6, ls=":")
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$\mathbb{E}[\log S_t] - \log S_0$")
+    ax.set_title(r"Itō correction: $\mathbb{E}[\log S_t]$ always below the naïve slope $\mu t$")
+    ax.legend(fontsize=9)
+    save("fig_ito_correction")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3  FOKKER-PLANCK
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_fokker_planck():
+    x = np.linspace(-0.5, 5.5, 600)
+    mu_drift, sigma_diff = 0.8, 0.3
+    times  = [0.05, 0.5, 1.5]
+    colors = [C3, C2, C0]
+    labels = [r"$t = 0.05$  (narrow spike)",
+              r"$t = 0.50$",
+              r"$t = 1.50$  (wide, drifted)"]
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    for t, col, lbl in zip(times, colors, labels):
+        mean = mu_drift * t
+        std  = sigma_diff * np.sqrt(t)
+        y = norm.pdf(x, mean, std)
+        ax.plot(x, y, color=col, lw=2, label=lbl)
+        ax.fill_between(x, y, alpha=0.10, color=col)
+
+    ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$p(t, x)$")
+    ax.set_title(r"Fokker-Planck: density drifts $(\mu=0.8)$ and broadens $(\sigma=0.3)$")
+    ax.legend(fontsize=9)
+    save("fig_fokker_planck")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3  EULER-MARUYAMA vs MILSTEIN
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_em_milstein():
+    dts     = np.array([0.1, 0.05, 0.02, 0.01, 0.005, 0.001])
+    em_err  = 0.38 * dts**0.5
+    mil_err = 0.19 * dts**1.0
+
+    fig, ax = plt.subplots(figsize=(6, 4))
+    ax.loglog(dts, em_err,  "o-",  color=C0, lw=2, ms=7,
+              label=r"Euler-Maruyama  (order $1/2$)")
+    ax.loglog(dts, mil_err, "s--", color=C1, lw=2, ms=7,
+              label=r"Milstein          (order $1$)")
+    ax.set_xlabel(r"Step size $\Delta t$")
+    ax.set_ylabel(r"Strong error  $\|X_T - \hat{X}_T\|$")
+    ax.set_title("SDE numerical schemes — strong convergence order")
+    ax.legend(fontsize=10); ax.grid(True, which="both", alpha=0.3)
+    save("fig_em_milstein")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.4  ORNSTEIN-UHLENBECK
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_ou_path():
+    rng = np.random.default_rng(3)
+    T, n, dt = 5.0, 2000, 0.0025
+    kappa, theta, sigma = 3.0, 0.5, 0.4
+    X = np.zeros(n); X[0] = 2.0
+    for i in range(1, n):
+        X[i] = X[i-1] + kappa * (theta - X[i-1]) * dt + sigma * rng.normal(0, np.sqrt(dt))
+
+    ts = np.linspace(0, T, n)
+    sig_inf = sigma / np.sqrt(2 * kappa)
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(ts, X, color=C0, lw=1.0, alpha=0.9, label=r"$X_t$")
+    ax.axhline(theta, color=C1, lw=1.8, ls="--",
+               label=fr"$\theta = {theta}$  (long-run mean)")
+    ax.fill_between(ts,
+                    theta - 2 * sig_inf,
+                    theta + 2 * sig_inf,
+                    alpha=0.10, color=GRAY, label=r"$\theta \pm 2\sigma_\infty$")
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$X_t$")
+    ax.set_title(fr"Ornstein-Uhlenbeck  ($\kappa={kappa},\;\theta={theta},\;\sigma={sigma}$)  — mean-reversion")
+    ax.legend(fontsize=9)
+    save("fig_ou_path")
+
+
+def fig_ou_transition():
+    x = np.linspace(-0.3, 2.6, 500)
+    kappa, theta, sigma, x0 = 3.0, 0.5, 0.4, 2.0
+    taus   = [0.1, 0.5, 2.0]
+    colors = [C3, C2, C0]
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    for tau, col in zip(taus, colors):
+        mean = theta + (x0 - theta) * np.exp(-kappa * tau)
+        var  = sigma**2 / (2 * kappa) * (1 - np.exp(-2 * kappa * tau))
+        y = norm.pdf(x, mean, np.sqrt(var))
+        ax.plot(x, y, color=col, lw=2,
+                label=fr"$\tau = {tau:.1f}$   (mean $= {mean:.2f}$)")
+        ax.fill_between(x, y, alpha=0.09, color=col)
+    ax.axvline(theta, color=C1, lw=1.3, ls="--", label=fr"$\theta = {theta}$")
+    ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$p(x_\tau \mid x_0)$")
+    ax.set_title(r"OU transition density: drifts toward $\theta$, widens over time")
+    ax.legend(fontsize=9)
+    save("fig_ou_transition")
+
+
+def fig_ou_loglik():
+    kappa_v = np.linspace(10, 120, 80)
+    theta_v = np.linspace(-0.005, 0.011, 80)
+    K, T    = np.meshgrid(kappa_v, theta_v)
+    Z = -(((K - 55) / 22)**2 + ((T - 0.003) / 0.003)**2)
+
+    fig, ax = plt.subplots(figsize=(6.2, 4.5))
+    cf = ax.contourf(theta_v * 1000, kappa_v, Z.T, levels=20, cmap="Blues")
+    ax.contour(theta_v * 1000, kappa_v, Z.T, levels=8,
+               colors="white", linewidths=0.7, alpha=0.55)
+    ax.plot(3, 55, "*", color=C1, ms=16, zorder=5,
+            label=r"MLE $\hat\theta, \hat\kappa$")
+    plt.colorbar(cf, ax=ax, label="Log-likelihood (normalised)")
+    ax.set_xlabel(r"$\theta \times 10^3$"); ax.set_ylabel(r"$\kappa$")
+    ax.set_title(r"OU log-likelihood surface  $\ell(\kappa, \theta \mid \hat\sigma)$")
+    ax.legend(fontsize=10)
+    save("fig_ou_loglik")
+
+
+def fig_ou_residuals():
+    rng = np.random.default_rng(9)
+    r = rng.normal(0, 1, 600)
+    x = np.linspace(-4, 4, 300)
+
+    fig, ax = plt.subplots(figsize=(6, 3.8))
+    ax.hist(r, bins=32, density=True, color=C0, alpha=0.50,
+            label="Standardised residuals")
+    ax.plot(x, norm.pdf(x), color=C1, lw=2.2,
+            label=r"$\mathcal{N}(0,1)$ theory")
+    ax.set_xlabel(r"$r_i$"); ax.set_ylabel("Density")
+    ax.set_title(r"OU residual diagnostic:  $r_i = (X_{t_i} - \hat\mu_i)/\hat\sigma$")
+    ax.legend(fontsize=9)
+    save("fig_ou_residuals")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §3  JUMP PROCESSES
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_poisson():
+    rng = np.random.default_rng(1)
+    lam, T = 2, 4.0
+    arrivals, t = [], 0.0
+    while True:
+        t += rng.exponential(1 / lam)
+        if t > T: break
+        arrivals.append(t)
+
+    ts = np.concatenate([[0.0], arrivals, [T]])
+    ns = np.arange(len(ts) - 1)
+
+    fig, ax = plt.subplots(figsize=(7, 3.5))
+    for i, (t0, t1, n) in enumerate(zip(ts[:-1], ts[1:], ns)):
+        ax.hlines(n, t0, t1, color=C0, lw=2.8)
+        if i < len(arrivals):
+            ax.vlines(t1, n, n + 1, color=C0, lw=2.0, linestyle=":")
+            ax.scatter([t1], [n],     s=45, color="white", edgecolors=C0, zorder=5, lw=1.5)
+            ax.scatter([t1], [n + 1], s=45, color=C0, zorder=5)
+
+    ax.yaxis.set_major_locator(mticker.MaxNLocator(integer=True))
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$N_t$")
+    ax.set_title(fr"Poisson process  ($\lambda = {lam}$ jumps/unit)  — inter-arrivals $\sim \mathrm{{Exp}}(\lambda)$")
+    save("fig_poisson")
+
+
+def fig_jump_diffusion():
+    rng = np.random.default_rng(11)
+    T, n, dt = 1.0, 1000, 0.001
+    mu, sigma, lam = 0.05, 0.18, 2.5
+    ts = np.linspace(0, T, n)
+    S  = np.ones(n)
+    jump_times = np.sort(rng.uniform(0, T, rng.poisson(lam * T)))
+
+    for i in range(1, n):
+        dW   = rng.normal(0, np.sqrt(dt))
+        S[i] = S[i-1] * np.exp((mu - 0.5 * sigma**2) * dt + sigma * dW)
+        if np.any((ts[i-1] < jump_times) & (jump_times <= ts[i])):
+            S[i] *= np.exp(rng.normal(0.0, 0.09))
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(ts, S, color=C0, lw=1.3, label=r"$S_t$  (jump-diffusion path)")
+    # mark jump locations
+    jt_idx = [np.searchsorted(ts, jt) for jt in jump_times if jt < T]
+    ax.scatter(ts[jt_idx], S[jt_idx], s=50, color=C3, zorder=5,
+               label=r"Poisson jump $\tau_k$", marker="v")
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$S_t$")
+    ax.set_title(r"Merton jump-diffusion  ($\lambda = 2.5$/yr,  $\sigma_J = 9\%$)")
+    ax.legend(fontsize=9)
+    save("fig_jump_diffusion")
+
+
+def fig_levy_tails():
+    x = np.linspace(0.05, 5, 600)
+    gauss_tail = norm.pdf(x)
+    gauss_tail /= gauss_tail[0]
+    vg_tail = np.exp(-1.5 * x) / x
+    vg_tail /= vg_tail[0]
+    alpha_tail = x ** (-1.8)
+    alpha_tail /= alpha_tail[0]
+
+    fig, ax = plt.subplots(figsize=(6.5, 4))
+    ax.semilogy(x, gauss_tail, lw=2, color=C0,
+                label=r"Gaussian  ($\nu \equiv 0$)")
+    ax.semilogy(x, vg_tail, lw=2, color=C2,
+                label=r"Variance Gamma  ($\nu \propto e^{-c|z|}/|z|$)")
+    ax.semilogy(x, alpha_tail, lw=2, color=C1, ls="--",
+                label=r"$\alpha$-stable  ($\nu \propto |z|^{-1-\alpha}$,  heaviest)")
+    ax.set_xlabel(r"Jump size $|z|$")
+    ax.set_ylabel(r"Lévy density $\nu(dz)/dz$  (log scale)")
+    ax.set_title("Lévy measure tails — heavier tail = more frequent/larger jumps")
+    ax.legend(fontsize=9); ax.grid(True, which="both", alpha=0.25)
+    save("fig_levy_tails")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §6  KALMAN FILTER
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_kalman_covariance():
+    t    = np.linspace(0, 30, 300)
+    Pinf = 0.17
+    Pt   = Pinf + (1.0 - Pinf) * np.exp(-0.35 * t)
+
+    fig, ax = plt.subplots(figsize=(7, 3.5))
+    ax.plot(t, Pt, color=C0, lw=2, label=r"$P_t$  (error covariance)")
+    ax.axhline(Pinf, color=C1, lw=1.6, ls="--",
+               label=fr"$P_\infty \approx {Pinf}$  (steady-state)")
+    ax.fill_between(t, Pt, Pinf, alpha=0.10, color=C0)
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$P_t$")
+    ax.set_title(r"Kalman filter: error covariance converges exponentially to $P_\infty$")
+    ax.legend(fontsize=9); ax.set_ylim(0, 1.05)
+    save("fig_kalman_covariance")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §7  MCMC
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_mcmc_energy():
+    x = np.linspace(-5, 5, 600)
+    pi = 0.5 * norm.pdf(x, -1.5, 0.8) + 0.5 * norm.pdf(x, 1.5, 0.9)
+    U  = -np.log(pi + 1e-12)
+    U -= U.min()
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(x, U, color=C0, lw=2)
+    ax.fill_between(x, U, alpha=0.08, color=C0)
+    ax.scatter([-1.5, 1.5], [U[np.abs(x + 1.5).argmin()],
+                              U[np.abs(x - 1.5).argmin()]],
+               s=90, color=C2, zorder=5, label=r"modes of $\pi$")
+    saddle_i = np.abs(x).argmin()
+    ax.scatter([x[saddle_i]], [U[saddle_i]], s=90, color=C3,
+               zorder=5, marker="^", label="energy barrier")
+    ax.annotate(r"accept with $e^{-\Delta U}$",
+                (x[saddle_i] + 0.3, U[saddle_i] - 0.4),
+                (2.2, 1.2), fontsize=9, color=C3,
+                arrowprops=dict(arrowstyle="->", color=C3, lw=1.0))
+    ax.set_xlabel(r"$x$"); ax.set_ylabel(r"$U(x) = -\log\pi(x)$")
+    ax.set_title(r"MCMC energy landscape  (bimodal target $\pi$)")
+    ax.legend(fontsize=9)
+    save("fig_mcmc_energy")
+
+
+def fig_mcmc_trace():
+    rng = np.random.default_rng(42)
+    x_cur = -1.5
+    chain = [x_cur]
+    for _ in range(2999):
+        prop    = x_cur + rng.normal(0, 0.8)
+        pi_cur  = 0.5 * norm.pdf(x_cur, -1.5, 0.8) + 0.5 * norm.pdf(x_cur, 1.5, 0.9)
+        pi_prop = 0.5 * norm.pdf(prop,  -1.5, 0.8) + 0.5 * norm.pdf(prop,  1.5, 0.9)
+        x_cur   = prop if rng.random() < pi_prop / pi_cur else x_cur
+        chain.append(x_cur)
+    chain = np.array(chain)
+
+    fig, axes = plt.subplots(1, 2, figsize=(9, 3.8))
+    axes[0].plot(chain, lw=0.6, color=C0, alpha=0.8)
+    axes[0].axhline(0, color=GRAY, lw=0.7, ls=":")
+    axes[0].set_xlabel("Iteration"); axes[0].set_ylabel(r"$x_t$")
+    axes[0].set_title("Trace plot  — chain mixes between both modes")
+
+    x = np.linspace(-5, 5, 400)
+    true_pi = 0.5 * norm.pdf(x, -1.5, 0.8) + 0.5 * norm.pdf(x, 1.5, 0.9)
+    axes[1].hist(chain, bins=50, density=True, color=C0, alpha=0.50,
+                 label="MCMC samples")
+    axes[1].plot(x, true_pi, color=C1, lw=2.2, label=r"true $\pi(x)$")
+    axes[1].set_xlabel(r"$x$"); axes[1].set_ylabel("Density")
+    axes[1].set_title("Marginal distribution")
+    axes[1].legend(fontsize=9)
+    plt.tight_layout()
+    save("fig_mcmc_trace")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §9  INFORMATION THEORY
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_kl_asymmetry():
+    x = np.linspace(-10, 10, 800)
+    p = norm.pdf(x, 0, 1)
+    q = norm.pdf(x, 0, 4)
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(x, p, color=C0, lw=2,       label=r"$p = \mathcal{N}(0,1)$  (narrow)")
+    ax.plot(x, q, color=C1, lw=2, ls="--", label=r"$q = \mathcal{N}(0,4)$  (wide)")
+    ax.fill_between(x, p, alpha=0.12, color=C0)
+    ax.fill_between(x, q, alpha=0.08, color=C1)
+
+    dx = x[1] - x[0]
+    eps = 1e-12
+    kl_pq = float(np.sum(p * np.log((p + eps) / (q + eps))) * dx)
+    kl_qp = float(np.sum(q * np.log((q + eps) / (p + eps)) * dx))
+    ax.text(-9.5, 0.085,
+            fr"$D_{{KL}}(p\|q) \approx {kl_pq:.2f}$  (small: $q$ covers $p$)",
+            fontsize=9, color=C0)
+    ax.text(-9.5, 0.066,
+            fr"$D_{{KL}}(q\|p) \approx {kl_qp:.2f}$  (large: $p$ misses tails of $q$)",
+            fontsize=9, color=C1)
+    ax.set_xlabel(r"$x$"); ax.set_ylabel("Density")
+    ax.set_title(r"KL divergence asymmetry:  $D_{KL}(p\|q) \neq D_{KL}(q\|p)$")
+    ax.legend(fontsize=9)
+    save("fig_kl_asymmetry")
+
+
+def fig_fisher_curvature():
+    theta = np.linspace(-3, 3, 400)
+    sigma_vals = [0.5, 1.0, 2.0]
+    colors     = [C0, C2, C1]
+    labels     = [r"$\sigma=0.5$  (high $\mathcal{I}$, sharp peak)",
+                  r"$\sigma=1.0$",
+                  r"$\sigma=2.0$  (low $\mathcal{I}$, flat peak)"]
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    for s, col, lbl in zip(sigma_vals, colors, labels):
+        logL = -0.5 * (theta / s)**2 - np.log(s)
+        logL -= logL.max()
+        ax.plot(theta, logL, lw=2, color=col, label=lbl)
+
+    ax.axvline(0, color=GRAY, lw=0.8, ls=":")
+    ax.set_xlabel(r"$\theta$"); ax.set_ylabel(r"$\log\mathcal{L}(\theta \mid x_\mathrm{obs})$  (centred)")
+    ax.set_title(r"Fisher information = log-likelihood curvature at $\theta^*$")
+    ax.legend(fontsize=9); ax.set_ylim(-4.2, 0.3)
+    save("fig_fisher_curvature")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §10  DIFFERENTIAL GEOMETRY
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_curvatures():
+    fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
+
+    # K > 0  — converging geodesics
+    ax = axes[0]
+    ax.set_aspect("equal"); ax.axis("off")
+    theta_arc = np.linspace(0, np.pi, 200)
+    ax.plot(np.cos(theta_arc), np.sin(theta_arc), color=GRAY, lw=1.5, ls="--", alpha=0.35)
+    for ang in np.linspace(-0.45, 0.45, 7):
+        r = np.linspace(0, 1, 60)
+        ax.plot(r * np.sin(ang), r * np.cos(ang), color=C0, lw=1.5, alpha=0.75)
+    ax.scatter([0], [0], s=70, color=C1, zorder=5)
+    ax.text(0, -0.12, "meet at N pole", ha="center", fontsize=8, color=GRAY)
+    ax.set_title(r"$K > 0$  (sphere $S^2$)" + "\ngeodesics converge", fontsize=10)
+
+    # K = 0  — parallel
+    ax = axes[1]; ax.axis("off")
+    for y in np.linspace(-0.8, 0.8, 7):
+        ax.plot([-1, 1], [y, y], color=C0, lw=1.5)
+    ax.set_xlim(-1.3, 1.3); ax.set_ylim(-1.2, 1.2)
+    ax.text(0, -1.1, "remain equidistant", ha="center", fontsize=8, color=GRAY)
+    ax.set_title(r"$K = 0$  (flat $\mathbb{R}^2$)" + "\nparallel geodesics", fontsize=10)
+
+    # K < 0  — diverging
+    ax = axes[2]; ax.axis("off")
+    for ang in np.linspace(-0.55, 0.55, 7):
+        r = np.linspace(0, 1.2, 60)
+        scale = 1 + 0.55 * r
+        ax.plot(r * np.sin(ang * scale), r * np.cos(ang * scale), color=C0, lw=1.5, alpha=0.75)
+    ax.scatter([0], [0], s=70, color=C1, zorder=5)
+    ax.set_xlim(-1.1, 1.1); ax.set_ylim(-0.15, 1.5)
+    ax.text(0, -0.12, "spread exponentially", ha="center", fontsize=8, color=GRAY)
+    ax.set_title(r"$K < 0$  (hyperbolic $H^2$)" + "\ngeodesics diverge", fontsize=10)
+
+    plt.suptitle("Sectional curvature determines geodesic behaviour", y=1.03, fontsize=12)
+    plt.tight_layout()
+    save("fig_curvatures")
+
+
+def fig_natural_gradient():
+    fig, axes = plt.subplots(1, 2, figsize=(9, 3.8))
+    theta1 = np.linspace(-2, 2, 300)
+    theta2 = np.linspace(-2, 2, 300)
+    T1, T2 = np.meshgrid(theta1, theta2)
+
+    # Standard: elongated contours → zigzag
+    Z_std = 6 * T1**2 + T2**2
+    axes[0].contour(T1, T2, Z_std, levels=7, colors=GRAY, alpha=0.45, linewidths=0.9)
+    path_std = [(1.6, 1.6), (0.05, 1.1), (0.75, 0.15), (0.03, 0.06), (0, 0)]
+    xs, ys = zip(*path_std)
+    axes[0].plot(xs, ys, "o-", color=C0, lw=1.8, ms=5)
+    axes[0].scatter([0], [0], s=120, color=C1, zorder=5, marker="*")
+    axes[0].set_title("Standard gradient $\\nabla_\\theta \\mathcal{L}$\n(zigzag on ill-conditioned $\\mathcal{I}$)",
+                      fontsize=10)
+    axes[0].set_xlabel(r"$\theta_1$"); axes[0].set_ylabel(r"$\theta_2$")
+
+    # Natural: circular contours → direct path
+    Z_nat = T1**2 + T2**2
+    axes[1].contour(T1, T2, Z_nat, levels=7, colors=GRAY, alpha=0.45, linewidths=0.9)
+    path_nat = [(1.6, 1.6), (0.8, 0.8), (0.3, 0.3), (0, 0)]
+    xs2, ys2 = zip(*path_nat)
+    axes[1].plot(xs2, ys2, "o-", color=C2, lw=1.8, ms=5)
+    axes[1].scatter([0], [0], s=120, color=C1, zorder=5, marker="*")
+    axes[1].set_title(r"Natural gradient $\mathcal{I}^{-1}\nabla_\theta\mathcal{L}$" + "\n(direct, reparametrisation-invariant)",
+                      fontsize=10)
+    axes[1].set_xlabel(r"$\theta_1$"); axes[1].set_ylabel(r"$\theta_2$")
+
+    plt.tight_layout()
+    save("fig_natural_gradient")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §2.3  PICARD ITERATION
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_picard():
+    t = np.linspace(0, 1.5, 300)
+    # True solution: dx = x dt  →  x(t) = e^t
+    x_true = np.exp(t)
+    # Picard iterates starting at x0 = 1
+    x0 = np.ones_like(t)      # n=0: constant 1
+    x1 = 1 + t                # n=1: linear
+    x2 = 1 + t + t**2 / 2    # n=2: quadratic
+    x3 = 1 + t + t**2/2 + t**3/6  # n=3
+
+    fig, ax = plt.subplots(figsize=(7, 3.8))
+    ax.plot(t, x0,    color=GRAY, lw=1.5, ls=":",  label=r"$X^{(0)}$: constant")
+    ax.plot(t, x1,    color=C3,   lw=1.5, ls="-.", label=r"$X^{(1)}$: linear")
+    ax.plot(t, x2,    color=C2,   lw=1.5, ls="--", label=r"$X^{(2)}$: quadratic")
+    ax.plot(t, x3,    color=C1,   lw=1.8,          label=r"$X^{(3)}$")
+    ax.plot(t, x_true, color=C0,  lw=2.2,          label=r"$X^{(\infty)} = e^t$  (true)")
+    ax.set_xlabel(r"$t$"); ax.set_ylabel(r"$X^{(n)}_t$")
+    ax.set_title(r"Picard iteration  ($dX = X\,dt$,  $X_0 = 1$)  — successive approximations")
+    ax.legend(fontsize=9); ax.set_ylim(0.8, 5.0)
+    save("fig_picard")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §8  HMM REGIME STATE MACHINE  (K = 3)
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_hmm_regime():
+    from matplotlib.patches import FancyBboxPatch, FancyArrowPatch
+    import matplotlib.patheffects as pe
+
+    fig, ax = plt.subplots(figsize=(9, 4.2))
+    ax.set_xlim(0, 9); ax.set_ylim(0, 4); ax.axis("off")
+
+    states = [
+        (1.5, 2.6, "State 1\nBull",    C2),
+        (4.5, 2.6, "State 2\nNeutral", GRAY),
+        (7.5, 2.6, "State 3\nBear",    C3),
+    ]
+    box_w, box_h = 2.0, 1.1
+    for (cx, cy, label, col) in states:
+        fancy = FancyBboxPatch((cx - box_w/2, cy - box_h/2), box_w, box_h,
+                               boxstyle="round,pad=0.08", linewidth=1.6,
+                               edgecolor=col, facecolor=col + "22",
+                               zorder=2)
+        ax.add_patch(fancy)
+        ax.text(cx, cy, label, ha="center", va="center", fontsize=10,
+                fontweight="bold", color=col, zorder=3)
+
+    # Forward arrows A₁₂, A₂₃
+    for x0, x1, label in [(2.5, 3.5, r"$A_{12}$"), (5.5, 6.5, r"$A_{23}$")]:
+        ax.annotate("", xy=(x1, 2.85), xytext=(x0, 2.85),
+                    arrowprops=dict(arrowstyle="-|>", color=C0, lw=1.5))
+        ax.text((x0+x1)/2, 2.98, label, ha="center", fontsize=9, color=C0)
+    # Backward arrows A₂₁, A₃₂
+    for x0, x1, label in [(3.5, 2.5, r"$A_{21}$"), (6.5, 5.5, r"$A_{32}$")]:
+        ax.annotate("", xy=(x1, 2.35), xytext=(x0, 2.35),
+                    arrowprops=dict(arrowstyle="-|>", color=C1, lw=1.5))
+        ax.text((x0+x1)/2, 2.22, label, ha="center", fontsize=9, color=C1)
+
+    # Emission table
+    col_labels = ["State", r"$\mu$", r"$\sigma$", "Character"]
+    rows = [
+        ["Bull",    "+0.05", "0.12", "high return, low vol"],
+        ["Neutral", " 0.00", "0.18", "flat, medium vol"],
+        ["Bear",    "−0.08", "0.35", "crash, high vol"],
+    ]
+    row_colors = [[C2+"33", C2+"33", C2+"33", C2+"33"],
+                  [GRAY+"33", GRAY+"33", GRAY+"33", GRAY+"33"],
+                  [C3+"33", C3+"33", C3+"33", C3+"33"]]
+    tbl = ax.table(cellText=rows, colLabels=col_labels, loc="bottom",
+                   cellColours=row_colors, bbox=[0.05, 0.0, 0.90, 0.42])
+    tbl.auto_set_font_size(False); tbl.set_fontsize(9)
+    for (r, c), cell in tbl.get_celld().items():
+        cell.set_edgecolor("#cccccc")
+        if r == 0:
+            cell.set_facecolor(C0 + "33")
+            cell.set_text_props(fontweight="bold")
+
+    ax.set_title(r"HMM Regime State Machine ($K=3$)  —  Emission $B_k(y)=\mathcal{N}(\mu_k,\sigma_k^2)$",
+                 fontsize=11, pad=6)
+    plt.tight_layout()
+    save("fig_hmm_regime")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §8.2  VITERBI TRELLIS  (K=3, T=4)
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_viterbi_trellis():
+    from matplotlib.patches import Circle, FancyArrowPatch
+
+    K, T = 3, 4
+    state_labels = ["1 (Bull)", "2 (Neutral)", "3 (Bear)"]
+    map_path = {(1, 1), (1, 2)}  # state index 1 = "2 (Neutral)" at t=2,3 (0-indexed t)
+
+    fig, ax = plt.subplots(figsize=(8, 3.8))
+    ax.set_xlim(-0.5, T + 0.5); ax.set_ylim(-0.5, K - 0.3); ax.axis("off")
+
+    # x-positions: t=1..4 → 0.5, 1.5, 2.5, 3.5
+    xs = [0.6 * (t + 1) for t in range(T)]
+    ys = [K - 1 - k for k in range(K)]  # top = state 1
+
+    # Draw crossing / passing arrows (selective to show crossing)
+    arrow_kw = dict(arrowstyle="-|>", connectionstyle="arc3,rad=0.0",
+                    color=GRAY, lw=1.1, alpha=0.55)
+    cross_kw = dict(arrowstyle="-|>", connectionstyle="arc3,rad=0.18",
+                    color=GRAY, lw=1.1, alpha=0.45)
+    for t in range(T - 1):
+        for k in range(K):
+            for k2 in range(K):
+                rad = 0.0 if k == k2 else (0.18 if k2 > k else -0.18)
+                col = C0 if (k == 1 and k2 == 1 and t >= 1) else GRAY
+                alpha = 0.9 if col == C0 else 0.3
+                ax.annotate("", xy=(xs[t+1], ys[k2]), xytext=(xs[t], ys[k]),
+                            arrowprops=dict(arrowstyle="-|>",
+                                            connectionstyle=f"arc3,rad={rad}",
+                                            color=col, lw=1.2 if col == C0 else 0.8,
+                                            alpha=alpha))
+
+    # Draw nodes
+    r = 0.14
+    for k in range(K):
+        for t in range(T):
+            is_map = (k == 1 and 1 <= t <= 2)
+            fc = C0 if is_map else "white"
+            ec = C0 if is_map else GRAY
+            circ = Circle((xs[t], ys[k]), r, facecolor=fc, edgecolor=ec, lw=1.8, zorder=4)
+            ax.add_patch(circ)
+
+    # Labels on left
+    for k in range(K):
+        ax.text(-0.1, ys[k], state_labels[k], ha="right", va="center",
+                fontsize=9, color=C0 if k == 1 else "black")
+
+    # x-axis ticks
+    for t in range(T):
+        ax.text(xs[t], -0.35, f"$t={t+1}$", ha="center", va="top", fontsize=9)
+
+    # Legend
+    ax.scatter([], [], color=C0, s=80, label="● MAP (Viterbi) path", zorder=5)
+    ax.scatter([], [], facecolor="white", edgecolors=GRAY, s=80, label="○ other nodes", zorder=5)
+    ax.legend(loc="upper right", fontsize=9, framealpha=0.9)
+
+    ax.set_title(r"Viterbi Trellis  ($K=3$, $T=4$)  —  $\delta_t(k)=\max_j\,\delta_{t-1}(j)\,A_{jk}\,B_k(y_t)$",
+                 fontsize=11)
+    plt.tight_layout()
+    save("fig_viterbi_trellis")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §10.2  STANDARD VS NATURAL GRADIENT — PROPERTY COMPARISON
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_std_vs_nat_gradient():
+    from matplotlib.patches import FancyBboxPatch
+
+    fig, axes = plt.subplots(1, 2, figsize=(9, 3.0))
+
+    panels = [
+        ("Standard Gradient\n" + r"$\theta_{k+1} = \theta_k - \eta\nabla\mathcal{L}$",
+         ["Flat $\\mathbb{R}^d$ geometry",
+          "Ignores manifold curvature",
+          "Slow on ill-conditioned $\\mathcal{I}$",
+          "$O(\\kappa(\\mathcal{I}))$ iterations"],
+         C3, C3 + "18"),
+        ("Natural Gradient\n" + r"$\theta_{k+1} = \theta_k - \eta\,\mathcal{I}(\theta)^{-1}\nabla\mathcal{L}$",
+         ["Riemannian metric $\\mathcal{I}(\\theta)$",
+          "Adapts to manifold geometry",
+          "Reparametrisation-invariant",
+          "$O(1)$ on exp. families (MLE step)"],
+         C2, C2 + "18"),
+    ]
+
+    for ax, (title, props, border, bg) in zip(axes, panels):
+        ax.set_xlim(0, 1); ax.set_ylim(0, 1); ax.axis("off")
+        fancy = FancyBboxPatch((0.03, 0.04), 0.94, 0.92,
+                               boxstyle="round,pad=0.04", linewidth=2,
+                               edgecolor=border, facecolor=bg)
+        ax.add_patch(fancy)
+        ax.text(0.5, 0.87, title, ha="center", va="top", fontsize=10,
+                fontweight="bold", color=border, transform=ax.transAxes,
+                multialignment="center")
+        y = 0.68
+        for prop in props:
+            ax.text(0.12, y, "•  " + prop, ha="left", va="top", fontsize=9.5,
+                    transform=ax.transAxes, color="#222222")
+            y -= 0.17
+
+    fig.suptitle("Standard  vs  Natural Gradient — geometric properties", fontsize=11, y=1.02)
+    plt.tight_layout()
+    save("fig_std_vs_nat_gradient")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# §10.3  MATRIX LIE GROUP HIERARCHY
+# ════════════════════════════════════════════════════════════════════════════
+
+def fig_lie_group_hierarchy():
+    from matplotlib.patches import FancyBboxPatch, FancyArrowPatch
+
+    fig, ax = plt.subplots(figsize=(9, 4.6))
+    ax.set_xlim(0, 9); ax.set_ylim(0, 4.6); ax.axis("off")
+
+    nodes = {
+        "GL":  (4.5, 4.1,  r"$\mathrm{GL}(n,\mathbb{R})$" + "\nall invertible $n\times n$", C0),
+        "SL":  (1.8, 2.85, r"$\mathrm{SL}(n,\mathbb{R})$" + "\n" + r"$\det=1$",              C2),
+        "On":  (4.5, 2.85, r"$O(n)$" + "\n$R^\top R=I$",                                    C1),
+        "Sp":  (7.2, 2.85, r"$\mathrm{Sp}(2n,\mathbb{R})$" + "\npreserves " + r"$\omega$", C2),
+        "SO":  (4.5, 1.55, r"$\mathrm{SO}(n)$" + "\n" + r"$\det=+1$ (rotations)",           C2),
+        "Hn":  (1.8, 1.55, r"$H(n)$ Heisenberg" + "\nupper triangular",                      C3),
+    }
+    notes = {
+        "SO": "portfolio factor\nrotation, PCA",
+        "Sp": "Hamiltonian\nmechanics, PMP",
+        "Hn": "path-signature\nfeature maps",
+    }
+    edges = [("GL","SL"), ("GL","On"), ("GL","Sp"), ("On","SO")]
+
+    bw, bh = 2.2, 0.76
+    for key, (cx, cy, label, col) in nodes.items():
+        fbp = FancyBboxPatch((cx-bw/2, cy-bh/2), bw, bh,
+                             boxstyle="round,pad=0.07", lw=1.6,
+                             edgecolor=col, facecolor=col+"22", zorder=2)
+        ax.add_patch(fbp)
+        ax.text(cx, cy, label, ha="center", va="center", fontsize=8.5,
+                multialignment="center", color=col, fontweight="bold", zorder=3)
+        if key in notes:
+            ax.text(cx + bw/2 + 0.15, cy, notes[key], va="center",
+                    fontsize=7.5, color="#555555", fontstyle="italic")
+
+    for src, dst in edges:
+        sx, sy = nodes[src][0], nodes[src][1]
+        dx, dy = nodes[dst][0], nodes[dst][1]
+        ax.annotate("", xy=(dx, dy + bh/2 + 0.04), xytext=(sx, sy - bh/2 - 0.04),
+                    arrowprops=dict(arrowstyle="-|>", color=GRAY, lw=1.4))
+
+    ax.set_title("Matrix Lie Group Hierarchy — subgroup inclusions and finance applications",
+                 fontsize=11, pad=5)
+    plt.tight_layout()
+    save("fig_lie_group_hierarchy")
+
+
+# ════════════════════════════════════════════════════════════════════════════
+# RUN ALL
+# ════════════════════════════════════════════════════════════════════════════
+
+if __name__ == "__main__":
+    funcs = [
+        fig_de_mutation, fig_rastrigin,
+        fig_random_walk, fig_bm_fan, fig_gbm,
+        fig_ito_correction,
+        fig_picard,
+        fig_fokker_planck, fig_em_milstein,
+        fig_ou_path, fig_ou_transition, fig_ou_loglik, fig_ou_residuals,
+        fig_poisson, fig_jump_diffusion, fig_levy_tails,
+        fig_kalman_covariance,
+        fig_mcmc_energy, fig_mcmc_trace,
+        fig_kl_asymmetry, fig_fisher_curvature,
+        fig_curvatures, fig_natural_gradient,
+        # new §8 & §10 diagrams
+        fig_hmm_regime, fig_viterbi_trellis,
+        fig_std_vs_nat_gradient, fig_lie_group_hierarchy,
+    ]
+    for fn in funcs:
+        print(f"  {fn.__name__} ... ", end="", flush=True)
+        fn()
+        print("ok")
+    print(f"\nDone — {len(funcs)} SVGs saved to {OUT}")
@@ -0,0 +1,138 @@
+/* ─── Optimiz-rs Documentation Custom Styles ──────────────────────────────── */
+/* Brand: orange (#f97316) on furo dark theme */
+
+:root {
+    --brand-orange: #f97316;
+    --brand-orange-light: #fb923c;
+    --brand-orange-dim: rgba(249, 115, 22, 0.15);
+    --brand-orange-border: rgba(249, 115, 22, 0.35);
+}
+
+/* ── Mermaid diagram container ─────────────────────────────────────────────── */
+.mermaid {
+    background: transparent !important;
+    padding: 1.5rem 0;
+    text-align: center;
+    overflow-x: auto;
+}
+
+.mermaid svg {
+    max-width: 100%;
+    border-radius: 12px;
+    filter: drop-shadow(0 4px 16px rgba(249, 115, 22, 0.12));
+}
+
+/* ── Better admonitions ─────────────────────────────────────────────────────── */
+.admonition {
+    border-radius: 10px !important;
+    border-left-width: 4px !important;
+    margin: 1.5rem 0 !important;
+}
+
+.admonition.note {
+    border-left-color: var(--brand-orange) !important;
+    background: var(--brand-orange-dim) !important;
+}
+
+.admonition.tip {
+    border-left-color: #22c55e !important;
+    background: rgba(34, 197, 94, 0.1) !important;
+}
+
+.admonition.warning {
+    border-left-color: #f59e0b !important;
+    background: rgba(245, 158, 11, 0.1) !important;
+}
+
+.admonition.important {
+    border-left-color: #ef4444 !important;
+    background: rgba(239, 68, 68, 0.1) !important;
+}
+
+/* ── Section headers ────────────────────────────────────────────────────────── */
+h1, h2, h3 {
+    scroll-margin-top: 4rem;
+}
+
+h2 {
+    border-bottom: 2px solid var(--brand-orange-border);
+    padding-bottom: 0.4rem;
+}
+
+/* ── Code blocks ────────────────────────────────────────────────────────────── */
+.highlight {
+    border-radius: 8px !important;
+    border: 1px solid rgba(255, 255, 255, 0.08) !important;
+}
+
+/* Inline code */
+code.literal {
+    background: rgba(249, 115, 22, 0.08) !important;
+    border: 1px solid var(--brand-orange-border) !important;
+    padding: 0.1em 0.4em !important;
+    border-radius: 4px !important;
+    color: var(--brand-orange-light) !important;
+    font-size: 0.88em !important;
+}
+
+/* ── Tables ─────────────────────────────────────────────────────────────────── */
+table.docutils {
+    border-collapse: collapse;
+    width: 100%;
+    margin: 1.2rem 0;
+    border-radius: 8px;
+    overflow: hidden;
+    font-size: 0.9rem;
+}
+
+table.docutils th {
+    background: var(--brand-orange-dim) !important;
+    color: var(--brand-orange-light) !important;
+    font-weight: 700;
+    padding: 0.6rem 0.9rem;
+    border-bottom: 2px solid var(--brand-orange-border);
+}
+
+table.docutils td {
+    padding: 0.5rem 0.9rem;
+    border-bottom: 1px solid rgba(255, 255, 255, 0.06);
+}
+
+table.docutils tr:hover td {
+    background: rgba(249, 115, 22, 0.04);
+}
+
+/* ── Performance metric boxes (for benchmarks) ──────────────────────────────── */
+.perf-box {
+    background: linear-gradient(135deg, rgba(249,115,22,0.12) 0%, rgba(249,115,22,0.04) 100%);
+    border: 1px solid var(--brand-orange-border);
+    border-radius: 10px;
+    padding: 1rem 1.5rem;
+    margin: 1rem 0;
+    display: inline-block;
+    min-width: 140px;
+    text-align: center;
+}
+
+.perf-box .val {
+    font-size: 1.8rem;
+    font-weight: 800;
+    color: var(--brand-orange);
+}
+
+.perf-box .lbl {
+    font-size: 0.75rem;
+    color: #94a3b8;
+    margin-top: 0.2rem;
+}
+
+/* ── Navigation sidebar ─────────────────────────────────────────────────────── */
+.sidebar-tree .current > .reference {
+    color: var(--brand-orange) !important;
+}
+
+/* ── Mobile ─────────────────────────────────────────────────────────────────── */
+@media (max-width: 768px) {
+    .mermaid svg { width: 100% !important; }
+    table.docutils { font-size: 0.78rem; }
+}
@@ -0,0 +1,167 @@
+Agent-based — bounded-confidence consensus
+==========================================
+
+Generic symmetric *interacting-agent* simulator implementing the linear bounded-confidence
+update rule
+
+.. math::
+
+   s^{k+1}_i \;=\; (1 - \alpha)\, s^k_i \;+\; \alpha\, \bar s^k \;+\; \xi^k_i,
+   \qquad \bar s^k \;=\; \frac1N \sum_{j=1}^N s^k_j,
+   \qquad \xi^k_i \sim \mathcal{N}(0, \sigma^2),
+
+with :math:`\alpha \in (0, 1]` the *averaging weight* and :math:`\sigma` the noise scale.  This is the
+DeGroot–Friedkin–Johnsen baseline of opinion dynamics, and the *complete-graph* limit of the
+Hegselmann–Krause and Vicsek flocking models.
+
+Mathematical background
+-----------------------
+
+**Mean conservation.**  Averaging the update over :math:`i` gives
+:math:`\bar s^{k+1} = \bar s^k + \bar\xi^k` with :math:`\mathbb{E}[\bar\xi^k] = 0`, so the empirical mean
+is a *martingale* and is exactly preserved in expectation:
+
+.. math::
+
+   \mathbb{E}[\bar s^k] \;=\; \bar s^0 \quad \text{for all } k \ge 0.
+
+In the noiseless case :math:`\sigma = 0` the mean is preserved *path-by-path*.
+
+**Geometric contraction of the spread.**  Define the deviation :math:`d^k_i := s^k_i - \bar s^k`.
+The update implies
+
+.. math::
+
+   d^{k+1}_i \;=\; (1 - \alpha)\, d^k_i \;+\; \bigl(\xi^k_i - \bar\xi^k\bigr) ,
+
+so in the absence of noise :math:`\| d^k \|_\infty \le (1 - \alpha)^k \| d^0 \|_\infty` — the spread
+*contracts geometrically* with rate :math:`1 - \alpha`.  The companion notebook plots
+:math:`\max_i s^k_i - \min_i s^k_i` on a log scale across :math:`\alpha \in \{0.05, \dots, 1\}` and
+recovers exactly this slope.
+
+**Stationary variance with noise.**  Treating the deviation as an AR(1) process with input
+variance :math:`\sigma^2 (1 - 1/N)`, the steady-state variance of any single agent's deviation is
+
+.. math::
+
+   \mathrm{Var}_\infty(d_i) \;=\; \frac{\sigma^2 (1 - 1/N)}{1 - (1 - \alpha)^2}
+     \;\xrightarrow[\alpha \to 0]{}\; \frac{\sigma^2}{2\alpha}\,(1 - 1/N).
+
+**Continuous-time limit (linear Vlasov).**  Sending :math:`\alpha = \theta\, \Delta t`,
+:math:`\xi^k_i = \sigma \sqrt{\Delta t}\, W^i_k` and :math:`\Delta t \to 0` recovers the McKean–Vlasov SDE
+:math:`dX^i_t = \theta(\bar X_t - X^i_t)\, dt + \sigma\, dW^i_t` of :doc:`mckean_vlasov` — the
+discrete consensus update is the prototype of mean-field interaction.
+
+**Spectral interpretation.**  On a general weighted graph the update reads
+:math:`s^{k+1} = (I - \alpha L)\, s^k + \xi^k`, where :math:`L` is the normalised Laplacian.  The
+complete-graph case shipped here has :math:`L = I - \tfrac1N \mathbf{1}\mathbf{1}^\top` with
+eigenvalue :math:`1` on the orthogonal complement of :math:`\mathbf{1}`, hence the contraction rate
+:math:`1 - \alpha` above.  Replacing :math:`\mathbf{1}\mathbf{1}^\top / N` by an arbitrary stochastic
+matrix produces the full DeGroot model and is a one-liner extension on the Rust side.
+
+Why it matters
+--------------
+
+* **Opinion dynamics & social learning.**  Calibration of polarisation/consensus models
+  (Bayesian persuasion, social media echo chambers, voting-system stability).
+* **Distributed estimation & federated learning.**  Average-consensus protocols for sensor
+  networks, gossip algorithms, federated averaging — all reduce to the same contraction
+  argument with explicit convergence rate :math:`1 - \alpha`.
+* **Coupled-oscillator physics.**  Linear approximation of the Kuramoto / Vicsek models near
+  the synchronised regime; direct comparison with the McKean–Vlasov continuous limit.
+
+.. note::
+   📓 **Companion notebook** — `view on GitHub <https://github.com/ThotDjehuty/optimiz-rs/blob/main/examples/notebooks/15_agent_based.ipynb>`_
+   · `download .ipynb <https://raw.githubusercontent.com/ThotDjehuty/optimiz-rs/main/examples/notebooks/15_agent_based.ipynb>`_
+
+15 — Agent-based dynamics
+=========================
+
+.. code-block:: python
+
+   import numpy as np
+   import matplotlib.pyplot as plt
+   from optimizr import _core as opt
+   plt.rcParams['figure.figsize'] = (7, 4)
+   plt.rcParams['figure.dpi'] = 110
+
+.. code-block:: python
+
+   init = np.arange(40.0).tolist()
+   init_mean = float(np.mean(init))
+   res = opt.consensus_dynamics(init, alpha=0.3, noise_sigma=0.1,
+                                 n_steps=80, seed=0)
+   n_t = res['n_steps']; n_a = res['n_agents']
+   S = np.array(res['states_flat']).reshape(n_t, n_a)
+   mean_traj = np.array(res['mean_trajectory'])
+   print('initial mean =', init_mean)
+   print('final mean   =', mean_traj[-1])
+   print('final std    =', float(S[-1].std()))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   for i in range(n_a):
+       ax.plot(S[:, i], color='tab:blue', alpha=0.3, lw=0.6)
+   ax.plot(mean_traj, color='red', lw=2, label='empirical mean')
+   ax.axhline(init_mean, color='k', ls=':', label='initial mean')
+   ax.set_xlabel('step k'); ax.set_ylabel('s^k_i'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Bounded-confidence consensus, α = 0.3')
+   fig.tight_layout(); plt.show()
+
+
+
+
+
+
+
+.. AUTO-PLOT-BEGIN
+.. image:: ../_static/auto/algorithms__agent_based/block_03_fig_01.png
+   :align: center
+   :width: 80%
+
+.. AUTO-PLOT-END
+.. image:: ../_static/v2/agent_based/plot_01.png
+   :align: center
+   :width: 80%
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   for alpha in [0.05, 0.1, 0.3, 0.6, 1.0]:
+       r = opt.consensus_dynamics(init, alpha=alpha, noise_sigma=0.0, n_steps=60, seed=0)
+       S = np.array(r['states_flat']).reshape(r['n_steps'], r['n_agents'])
+       spread = S.max(axis=1) - S.min(axis=1)
+       ax.semilogy(spread, label=f'α = {alpha:g}')
+   ax.set_xlabel('step k'); ax.set_ylabel('max_i s − min_i s')
+   ax.set_title('Convergence rate vs averaging weight α'); ax.legend(); ax.grid(alpha=0.3)
+   fig.tight_layout(); plt.show()
+
+
+
+
+
+
+
+.. AUTO-PLOT-BEGIN
+.. image:: ../_static/auto/algorithms__agent_based/block_04_fig_01.png
+   :align: center
+   :width: 80%
+
+.. AUTO-PLOT-END
+.. image:: ../_static/v2/agent_based/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified:** without noise, the empirical mean is exactly preserved and the spread decays geometrically.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn simulate_agent_based<T>(initial: &[f64], transition: T, cfg: &AgentBasedConfig) -> Result<AgentBasedResult>
+   where T: Fn(f64, &[f64], usize) -> f64;
+
+   pub struct AgentBasedConfig { pub n_agents: usize, pub n_steps: usize, pub noise_sigma: f64, pub seed: u64 }
+   pub struct AgentBasedResult { pub states: Array2<f64>, pub mean_trajectory: Array1<f64> }
--- a/Show More
+++ b/Show More