release(v2.0.0-alpha.2): PyO3 bindings + executed companion notebooks + Sphinx RST with inline plots

PyO3 abi3 bindings for the 13 v2.0.0 functions across 8 module groups:
  bsde, pde, stochastic_control, optimal_control::quadratic_impact_control,
  mean_field::mckean_vlasov, agent_based, inference, optimization.

8 executed companion notebooks under examples/notebooks/10_bsde.ipynb …
17_generative_calibration.ipynb (cell outputs and matplotlib figures
preserved as proof-of-work; verified against analytic ground truths).

8 Sphinx RST pages under docs/source/algorithms/{bsde,pde,stochastic_control,
quadratic_impact_control,mckean_vlasov,agent_based,robust_drift,
generative_calibration_hooks}.rst with .. math:: derivations and inline
.. image:: directives placed immediately after each .. code-block:: python
so each plot appears directly under the code that produced it.

18 PNG plot assets under docs/source/_static/v2/<group>/.

index.rst extended with a new 'v2.0 Generic Stochastic Control & PDE'
toctree caption.

Forbidden-vocabulary audit on new src/, docs/source/algorithms/ and
binding files: zero matches.

All previously stable APIs untouched; v2.0.0 is additive at the binding
level — no v1.x function signature was changed.

docs/source/algorithms/agent_based.rst

@@ -0,0 +1,74 @@
+Agent-based — bounded-confidence consensus
+==========================================
+
+Generic interacting-agent simulator (`consensus_dynamics`) — linear bounded-confidence rule $s_i^{k+1} = (1-α) s_i^k + α \bar s^k + ξ_i$.
+
+.. note:: Companion executed notebook: `15_agent_based.ipynb <../../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()
+
+.. 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()
+
+.. 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> }

docs/source/algorithms/bsde.rst

@@ -0,0 +1,99 @@
+BSDE — θ-scheme and deep-BSDE bridge
+====================================
+
+This notebook exercises `optimizr.linear_bsde_constant_coeffs`, the Crank–Nicolson θ-scheme for the BSDE
+`-dY = (a Y + b Z + c) dt - Z dW` with constant coefficients, and verifies the discrete trajectory against the analytic solution `Y_t = exp(-ρ (T - t))`.
+
+.. note:: Companion executed notebook: `10_bsde.ipynb <../../examples/notebooks/10_bsde.ipynb>`_
+
+10 — BSDE θ-scheme
+==================
+
+Generic CPU-only Crank–Nicolson scheme for linear backward stochastic differential equations.  Reference doc page: [bsde.rst](../../docs/source/algorithms/bsde.rst).
+
+.. 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
+
+Exponential ground-truth check
+------------------------------
+
+With $a(t) \equiv -\rho$, $b = c = 0$ and $Y_T = 1$ the analytic deterministic solution is $Y_t = e^{-\rho (T-t)}$.
+
+.. code-block:: python
+
+   rho = 0.3
+   T   = 1.0
+   res = opt.linear_bsde_constant_coeffs(
+       a_const=-rho, b_const=0.0, c_const=0.0,
+       terminal=1.0, n_steps=200, t_horizon=T, theta=0.5,
+   )
+   tg = np.array(res['time_grid'])
+   yg = np.array(res['y'])
+   analytic = np.exp(-rho * (T - tg))
+   print('Y0 =', yg[0], '   exp(-rho T) =', analytic[0])
+   print('max abs error =', float(np.max(np.abs(yg - analytic))))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.plot(tg, yg, label='θ-scheme', lw=2)
+   ax.plot(tg, analytic, '--', label='analytic exp(-ρ(T-t))')
+   ax.set_xlabel('t'); ax.set_ylabel('Y_t')
+   ax.set_title('Linear BSDE — Crank–Nicolson vs analytic')
+   ax.legend(); ax.grid(alpha=0.3)
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/bsde/plot_01.png
+   :align: center
+   :width: 80%
+
+Convergence rate study
+----------------------
+
+Crank–Nicolson is second-order in `Δt`.
+
+.. code-block:: python
+
+   errs = []
+   ns = [25, 50, 100, 200, 400, 800]
+   for n in ns:
+       r = opt.linear_bsde_constant_coeffs(-rho, 0.0, 0.0, 1.0, n, T, 0.5)
+       errs.append(abs(r['y'][0] - np.exp(-rho * T)))
+   print(list(zip(ns, errs)))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.loglog(ns, errs, 'o-')
+   ax.loglog(ns, [errs[0] * (ns[0] / n) ** 2 for n in ns],
+             ':', label='O(Δt²) reference')
+   ax.set_xlabel('n_steps'); ax.set_ylabel('|Y0 − analytic|')
+   ax.set_title('Crank–Nicolson convergence'); ax.grid(which='both', alpha=0.3); ax.legend()
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/bsde/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified against analytic ground truth:** `Y_t = exp(-ρ (T - t))` — relative error at `t = 0` below `1e-3` for `n_steps = 200`.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn solve_linear_bsde<A, B, C>(
+       a: A, b: B, c: C, terminal: f64, cfg: &ThetaSchemeConfig
+   ) -> Result<ThetaSchemeResult>
+   where A: Fn(f64) -> f64, B: Fn(f64) -> f64, C: Fn(f64) -> f64;
+
+   pub struct ThetaSchemeConfig { pub n_steps: usize, pub t_horizon: f64, pub theta: f64 }
+   pub struct ThetaSchemeResult { pub y: Array1<f64>, pub z: Array1<f64>, pub time_grid: Array1<f64> }
+
+   pub trait ConditionalExpectation { /* deep-BSDE bridge */ }
+   pub struct DeepBsdeBridge { /* ... */ }

docs/source/algorithms/generative_calibration_hooks.rst

@@ -0,0 +1,66 @@
+Generative calibration — Gaussian MMD loss
+==========================================
+
+Maximum-Mean-Discrepancy distance with Gaussian kernel (`mmd_gaussian`).  Self-distance is exactly zero; the metric grows monotonically with sample shift.
+
+.. note:: Companion executed notebook: `17_generative_calibration.ipynb <../../examples/notebooks/17_generative_calibration.ipynb>`_
+
+17 — MMD calibration loss
+=========================
+
+.. 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
+
+   x = np.linspace(0.0, 5.0, 80)
+   shifts = np.linspace(0.0, 6.0, 40)
+   d = [opt.mmd_gaussian(x.tolist(), (x + s).tolist(), 1.0) for s in shifts]
+   print('MMD self =', d[0])
+   print('MMD at shift 6.0 =', d[-1])
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.plot(shifts, d, lw=2)
+   ax.set_xlabel('translation Δ'); ax.set_ylabel('MMD(P, P + Δ)')
+   ax.set_title('Gaussian-kernel MMD vs translation (σ = 1)')
+   ax.grid(alpha=0.3); fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/generative_calibration_hooks/plot_01.png
+   :align: center
+   :width: 80%
+
+Bandwidth dependence
+--------------------
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   for sigma in [0.25, 0.5, 1.0, 2.0]:
+       d = [opt.mmd_gaussian(x.tolist(), (x + s).tolist(), sigma) for s in shifts]
+       ax.plot(shifts, d, label=f'σ = {sigma:g}')
+   ax.set_xlabel('translation Δ'); ax.set_ylabel('MMD'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('MMD as a function of kernel bandwidth')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/generative_calibration_hooks/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified:** `MMD(x, x) = 0`; metric is strictly monotonic in shift.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn mmd_distance(x: &[f64], y: &[f64], loss: &MmdLoss) -> Result<f64>;
+   pub fn calibration_step<S: GenerativeSampler>(sampler: &mut S, target: &[f64], loss: &MmdLoss, lr: f64) -> Result<f64>;
+   pub trait GenerativeSampler { fn sample(&self, n: usize, seed: u64) -> Vec<f64>; fn parameters(&self) -> Vec<f64>; fn perturb(&mut self, deltas: &[f64]); }
+   pub struct MmdLoss { pub sigma: f64 }

docs/source/algorithms/mckean_vlasov.rst

@@ -0,0 +1,72 @@
+McKean–Vlasov — propagation of chaos
+====================================
+
+Interacting-particle Euler scheme for $dX_t = θ(\bar X_t - X_t) dt + σ dW_t$ (`mean_reverting_mckean_vlasov`).  The empirical mean is preserved; the empirical variance approaches the diffusion-only equilibrium.
+
+.. note:: Companion executed notebook: `14_mckean_vlasov.ipynb <../../examples/notebooks/14_mckean_vlasov.ipynb>`_
+
+14 — McKean–Vlasov mean-reverting 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.linspace(-2.0, 2.0, 200).tolist()
+   init_mean = float(np.mean(init))
+   res = opt.mean_reverting_mckean_vlasov(
+       initial=init, theta=1.0, sigma=0.1,
+       n_steps=1000, t_horizon=1.0, seed=42,
+   )
+   n_t = res['n_steps']; n_p = res['n_particles']
+   X   = np.array(res['paths_flat']).reshape(n_t, n_p)
+   tg  = np.array(res['time_grid'])
+   print('initial mean =', init_mean)
+   print('final  mean  =', float(X[-1].mean()))
+   print('final  std   =', float(X[-1].std()))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.plot(tg, X[:, ::20], color='tab:blue', alpha=0.2, lw=0.6)
+   ax.plot(tg, X.mean(axis=1), color='red', lw=2, label='empirical mean')
+   ax.axhline(init_mean, color='k', ls=':', label='initial mean')
+   ax.set_xlabel('t'); ax.set_ylabel('X^i_t'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Mean-reverting McKean–Vlasov — 200 particles')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/mckean_vlasov/plot_01.png
+   :align: center
+   :width: 80%
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.hist(X[0],  bins=30, alpha=0.5, label='t = 0',  density=True)
+   ax.hist(X[-1], bins=30, alpha=0.5, label='t = T',  density=True)
+   ax.set_xlabel('x'); ax.set_ylabel('empirical density'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Marginal density at t = 0 and t = T')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/mckean_vlasov/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified:** empirical mean stays within `0.05` of the initial mean.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn simulate_mckean_vlasov<B>(initial: &[f64], drift: B, cfg: &McKeanVlasovConfig) -> Result<McKeanVlasovResult>
+   where B: Fn(f64, &[f64]) -> f64;
+
+   pub struct McKeanVlasovConfig { pub n_particles: usize, pub n_steps: usize, pub t_horizon: f64, pub sigma: f64, pub seed: u64 }
+   pub struct McKeanVlasovResult { pub paths: Array2<f64>, pub time_grid: Array1<f64> }

docs/source/algorithms/pde.rst

@@ -0,0 +1,125 @@
+PDE — Fokker–Planck, HJB, elliptic Poisson
+==========================================
+
+Three CPU-only finite-difference solvers: 1-D forward Fokker–Planck (`fokker_planck_constant`), 2-D explicit HJB (`hjb_quadratic_2d`) and 2-D Poisson SOR (`poisson_2d_zero_boundary`).  Each routine is verified against an analytic ground truth.
+
+.. note:: Companion executed notebook: `11_pde.ipynb <../../examples/notebooks/11_pde.ipynb>`_
+
+11 — PDE solvers
+================
+
+Fokker–Planck, HJB, Poisson.
+
+.. 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
+
+Pure-diffusion Fokker–Planck
+----------------------------
+
+$\partial_t m = \tfrac12 \partial_{xx} m$ with Gaussian initial density should remain centred and approximately Gaussian.
+
+.. code-block:: python
+
+   res = opt.fokker_planck_constant(
+       mu=0.0, sigma_sq=1.0, init_sigma=1.0,
+       x_min=-8.0, x_max=8.0, n_x=401,
+       t_horizon=0.5, n_t=8000,
+   )
+   x = np.array(res['x_grid'])
+   t = np.array(res['time_grid'])
+   nx = res['n_x']; nt = res['n_t']
+   M = np.array(res['density']).reshape(nt + 1, nx)
+   print('total mass at t=0:',  np.trapezoid(M[0], x))
+   print('total mass at t=T:',  np.trapezoid(M[-1], x))
+   print('mean   at t=T:',      np.trapezoid(x * M[-1], x))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   for k in [0, nt // 4, nt // 2, 3 * nt // 4, nt]:
+       ax.plot(x, M[k], label=f't = {t[k]:.2f}')
+   ax.set_xlim(-5, 5); ax.set_xlabel('x'); ax.set_ylabel('m(x, t)')
+   ax.set_title('Pure-diffusion Fokker–Planck'); ax.grid(alpha=0.3); ax.legend()
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/pde/plot_01.png
+   :align: center
+   :width: 80%
+
+2-D Poisson eigenfunction
+-------------------------
+
+$-\Delta u = 2\pi^2 \sin(\pi x)\sin(\pi y)$ on the unit square with zero Dirichlet boundary admits the exact solution $u(x,y) = \sin(\pi x)\sin(\pi y)$.
+
+.. code-block:: python
+
+   n = 65
+   xs = np.linspace(0, 1, n); ys = np.linspace(0, 1, n)
+   X, Y = np.meshgrid(xs, ys, indexing='ij')
+   F = 2 * np.pi ** 2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
+   res = opt.poisson_2d_zero_boundary(F.flatten().tolist(), n, n)
+   U = np.array(res['u']).reshape(n, n)
+   U_exact = np.sin(np.pi * X) * np.sin(np.pi * Y)
+   print('iterations =', res['iterations'])
+   print('residual   =', res['residual'])
+   print('max error  =', float(np.max(np.abs(U - U_exact))))
+
+.. code-block:: python
+
+   fig, axes = plt.subplots(1, 2, figsize=(11, 4))
+   im0 = axes[0].imshow(U.T, origin='lower', extent=(0, 1, 0, 1), cmap='viridis')
+   axes[0].set_title('SOR solution'); plt.colorbar(im0, ax=axes[0])
+   im1 = axes[1].imshow((U - U_exact).T, origin='lower', extent=(0, 1, 0, 1), cmap='RdBu_r')
+   axes[1].set_title('error vs analytic'); plt.colorbar(im1, ax=axes[1])
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/pde/plot_02.png
+   :align: center
+   :width: 80%
+
+2-D HJB with quadratic terminal
+-------------------------------
+
+Heat-only relaxation ($H = 0$, σ² > 0) preserves a constant value, while a quadratic terminal $g(x) = ½(x²+y²)$ smooths.
+
+.. code-block:: python
+
+   res = opt.hjb_quadratic_2d(n_per_dim=21, x_min=-1.0, x_max=1.0,
+                               n_t=200, t_horizon=0.2, sigma_sq=0.1)
+   ax_x = np.array(res['axis']); npd = res['n_per_dim']
+   V = np.array(res['value']).reshape(npd, npd)
+   print('V(0,0)   =', V[npd // 2, npd // 2])
+   print('V(±1,±1) =', V[0, 0], V[-1, -1])
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   im = ax.imshow(V.T, origin='lower', extent=(-1, 1, -1, 1), cmap='magma')
+   ax.set_title('HJB value V(0, x, y) — quadratic terminal')
+   plt.colorbar(im, ax=ax)
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/pde/plot_03.png
+   :align: center
+   :width: 80%
+
+**Verified:** Poisson max-error vs analytic eigenfunction below `5e-3`; Fokker–Planck mean stays at 0 within `0.05`.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn solve_fokker_planck_1d<F, G, H>(drift: F, diffusion_sq: G, initial_density: H, cfg: &FokkerPlanckConfig) -> Result<FokkerPlanckResult>
+   where F: Fn(f64) -> f64, G: Fn(f64) -> f64, H: Fn(f64) -> f64;
+
+   pub fn solve_hjb_multid<H, G>(hamiltonian: H, terminal: G, cfg: &HjbMultidConfig) -> Result<HjbMultidResult>
+   where H: Fn(&[f64], &[f64]) -> f64, G: Fn(&[f64]) -> f64;
+
+   pub fn solve_poisson_2d<F, G>(rhs: F, boundary: G, cfg: &EllipticFdConfig) -> Result<EllipticFdResult>
+   where F: Fn(f64, f64) -> f64, G: Fn(f64, f64) -> f64;

docs/source/algorithms/quadratic_impact_control.rst

@@ -0,0 +1,76 @@
+Quadratic-impact control — closed-form Riccati
+==============================================
+
+Closed-form Riccati feedback for a controlled 1-D SDE with quadratic running cost (`quadratic_impact_control_py`).
+
+.. note:: Companion executed notebook: `13_quadratic_impact.ipynb <../../examples/notebooks/13_quadratic_impact.ipynb>`_
+
+13 — Quadratic-impact controlled SDE
+====================================
+
+.. 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
+
+Riccati fixed-point check
+-------------------------
+
+$h'(t) = h(t)^2/γ - φ$ with $h(T) = A$.  When $γ = φ = A = 1$ the right-hand side is $h^2 - 1 = 0$ at $h = 1$, so `h ≡ 1`.
+
+.. code-block:: python
+
+   res = opt.quadratic_impact_control_py(
+       gamma=1.0, phi=1.0, a_terminal=1.0,
+       t_horizon=0.5, n_steps=500,
+   )
+   tg = np.array(res['time_grid'])
+   h  = np.array(res['h']); k = np.array(res['feedback_gain'])
+   print('h drift from 1:', float(np.max(np.abs(h - 1.0))))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.plot(tg, h, label='h(t)')
+   ax.plot(tg, k, '--', label='k(t) = h(t)/γ')
+   ax.axhline(1.0, color='k', alpha=0.3, ls=':', label='fixed point')
+   ax.set_xlabel('t'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Riccati fixed point  γ=φ=A=1')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/quadratic_impact_control/plot_01.png
+   :align: center
+   :width: 80%
+
+Sensitivity to the terminal weight
+----------------------------------
+
+Vary $A$, fix $γ = 1$, $φ = 0.25$, $T = 1$.
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   for A in [0.0, 0.25, 0.5, 1.0, 2.0, 5.0]:
+       r = opt.quadratic_impact_control_py(1.0, 0.25, A, 1.0, 1000)
+       ax.plot(r['time_grid'], r['h'], label=f'A = {A:g}')
+   ax.set_xlabel('t'); ax.set_ylabel('h(t)'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Riccati sensitivity to terminal weight')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/quadratic_impact_control/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified:** `h ≡ 1` with `max|h - 1| < 1e-9` at the fixed point.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn solve_quadratic_impact_control(cfg: &QuadraticImpactConfig) -> Result<QuadraticImpactResult>;
+   pub struct QuadraticImpactConfig { pub gamma: f64, pub phi: f64, pub a_terminal: f64, pub t_horizon: f64, pub n_steps: usize }
+   pub struct QuadraticImpactResult { pub time_grid: Array1<f64>, pub h: Array1<f64>, pub feedback_gain: Array1<f64> }

docs/source/algorithms/robust_drift.rst

@@ -0,0 +1,87 @@
+Inference — Huber-IRLS drift estimator
+======================================
+
+Robust drift estimator (`robust_drift`) for $x_{k+1} = x_k + (a + b x_k) Δt + σ ε_k$ via Huber IRLS — resists 5 % heavy-tailed innovations.
+
+.. note:: Companion executed notebook: `16_robust_drift.ipynb <../../examples/notebooks/16_robust_drift.ipynb>`_
+
+16 — Robust drift estimation
+============================
+
+.. 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
+
+Synthetic stationary process with 5 % outliers
+----------------------------------------------
+
+.. code-block:: python
+
+   rng = np.random.default_rng(7)
+   true_a, true_b = 1.0, -0.5
+   dt, n = 0.01, 5000
+   x = [0.0]
+   for k in range(n):
+       if k % 20 == 0:
+           eps = rng.uniform(-2.0, 2.0)
+       else:
+           eps = rng.uniform(-0.1, 0.1)
+       x.append(x[-1] + (true_a + true_b * x[-1]) * dt + eps * np.sqrt(dt))
+   x = np.array(x)
+   print('observation length =', len(x))
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.plot(x, lw=0.6)
+   ax.axhline(true_a / -true_b, color='red', ls='--', label='OU level a/(-b) = 2')
+   ax.set_xlabel('k'); ax.set_ylabel('x_k'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Synthetic series with heavy-tailed innovations')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/robust_drift/plot_01.png
+   :align: center
+   :width: 80%
+
+.. code-block:: python
+
+   res = opt.robust_drift(x.tolist(), dt=dt)
+   print(f'a (true 1.0)  ->  {res["a"]:.4f}')
+   print(f'b (true -0.5) ->  {res["b"]:.4f}')
+   print('IRLS iterations =', res['iterations'])
+
+.. code-block:: python
+
+   # Compare against a naïve OLS that is broken by outliers.
+   y = (x[1:] - x[:-1]) / dt
+   X = np.vstack([np.ones_like(x[:-1]), x[:-1]]).T
+   ols_ab, *_ = np.linalg.lstsq(X, y, rcond=None)
+   print('OLS a, b =', ols_ab)
+   fig, ax = plt.subplots()
+   labels = ['true', 'OLS', 'robust']
+   vals_a = [true_a, ols_ab[0], res['a']]
+   vals_b = [true_b, ols_ab[1], res['b']]
+   ax.bar(np.arange(3) - 0.2, vals_a, width=0.4, label='a')
+   ax.bar(np.arange(3) + 0.2, vals_b, width=0.4, label='b')
+   ax.set_xticks(range(3)); ax.set_xticklabels(labels)
+   ax.legend(); ax.grid(alpha=0.3); ax.set_title('Robust vs OLS drift estimate')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/robust_drift/plot_02.png
+   :align: center
+   :width: 80%
+
+**Verified:** Huber IRLS recovers `(a, b)` within `0.2` even with 5 % heavy outliers.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn estimate_robust_drift(observations: &[f64], cfg: &RobustDriftConfig) -> Result<RobustDriftResult>;
+   pub struct RobustDriftConfig { pub dt: f64, pub huber_delta: f64, pub max_iterations: usize, pub tolerance: f64 }
+   pub struct RobustDriftResult { pub a: f64, pub b: f64, pub iterations: usize }

docs/source/algorithms/stochastic_control.rst

@@ -0,0 +1,114 @@
+Stochastic control — switching, Pontryagin, two-sided intensities
+=================================================================
+
+Three primitives: discrete-time optimal switching (`optimal_switching_dp`), 1-D Pontryagin LQR shooting (`pontryagin_lqr`) and the bilateral intensity controller (`two_sided_intensities`).
+
+.. note:: Companion executed notebook: `12_stochastic_control.ipynb <../../examples/notebooks/12_stochastic_control.ipynb>`_
+
+12 — Stochastic control
+=======================
+
+.. 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
+
+Optimal switching (Snell envelope)
+----------------------------------
+
+Two modes; only mode 1 pays a unit reward.  Free switching should give `V_0(0) = N - 1` and `V_0(1) = N`.
+
+.. code-block:: python
+
+   n_steps, n_modes = 5, 2
+   stage = np.zeros((n_steps, n_modes)); stage[:, 1] = 1.0
+   cost = [0.0] * (n_modes * n_modes)
+   res = opt.optimal_switching_dp(stage.flatten().tolist(),
+                                   [0.0] * n_modes, cost,
+                                   n_modes, n_steps)
+   value  = np.array(res['value']).reshape(n_steps + 1, n_modes)
+   policy = np.array(res['policy']).reshape(n_steps + 1, n_modes)
+   print('V_0 =', value[0])
+   print('Optimal next mode at each (k, i):'); print(policy)
+
+.. code-block:: python
+
+   fig, ax = plt.subplots()
+   ax.step(range(n_steps + 1), value[:, 0], where='post', label='V_k(mode 0)')
+   ax.step(range(n_steps + 1), value[:, 1], where='post', label='V_k(mode 1)')
+   ax.set_xlabel('k'); ax.set_ylabel('value'); ax.legend(); ax.grid(alpha=0.3)
+   ax.set_title('Snell envelope — free switching')
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/stochastic_control/plot_01.png
+   :align: center
+   :width: 80%
+
+Pontryagin 1-D LQR
+------------------
+
+Closed-form Riccati for $a=q=0$, $b=r=s_T=1$, $T=1$ is $P(t) = 1/(1 + (T - t))$, hence $P(0) = 0.5$.
+
+.. code-block:: python
+
+   res = opt.pontryagin_lqr(a=0.0, b=1.0, q=0.0, r=1.0,
+                             s_terminal=1.0, x0=1.0,
+                             t_horizon=1.0, n_steps=2000)
+   tg = np.array(res['time_grid'])
+   P = np.array(res['riccati'])
+   x = np.array(res['state']); u = np.array(res['control'])
+   P_an = 1.0 / (1.0 + (1.0 - tg))
+   print('P(0) =', P[0], '   analytic =', P_an[0])
+   print('cost =', res['cost'])
+
+.. code-block:: python
+
+   fig, axes = plt.subplots(1, 3, figsize=(13, 4))
+   axes[0].plot(tg, P, label='numeric'); axes[0].plot(tg, P_an, '--', label='analytic')
+   axes[0].set_title('Riccati P(t)'); axes[0].set_xlabel('t'); axes[0].legend(); axes[0].grid(alpha=0.3)
+   axes[1].plot(tg, x); axes[1].set_title('state x(t)'); axes[1].set_xlabel('t'); axes[1].grid(alpha=0.3)
+   axes[2].plot(tg[:-1], u); axes[2].set_title('feedback u(t) = -(b/r) P(t) x(t)'); axes[2].set_xlabel('t'); axes[2].grid(alpha=0.3)
+   fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/stochastic_control/plot_02.png
+   :align: center
+   :width: 80%
+
+Two-sided intensity control
+---------------------------
+
+Affine premium $δ_±(λ) = α_± + κ_± λ$.  First-order condition: $\lambda^*_\pm = \max(0, (α_\pm - ΔV_\pm) / (2 κ_\pm))$.
+
+.. code-block:: python
+
+   deltas = np.linspace(-2.0, 2.0, 41)
+   lam_plus = []
+   for dv in deltas:
+       r = opt.two_sided_intensities(1.0, 1.0, 0.5, 0.5, dv, -dv)
+       lam_plus.append(r['lambda_plus'])
+   lam_plus = np.array(lam_plus)
+   fig, ax = plt.subplots()
+   ax.plot(deltas, lam_plus, lw=2)
+   ax.set_xlabel('ΔV_+'); ax.set_ylabel('λ*_+')
+   ax.set_title('Optimal upward intensity vs value-function gradient')
+   ax.grid(alpha=0.3); fig.tight_layout(); plt.show()
+
+.. image:: ../_static/v2/stochastic_control/plot_03.png
+   :align: center
+   :width: 80%
+
+**Verified:** switching `V_0` matches analytic recursion exactly; Pontryagin `P(0) = 0.4999` against analytic `0.5`.
+
+API
+---
+
+.. code-block:: rust
+
+   pub fn solve_optimal_switching<R, T>(stage_reward: R, terminal_payoff: T, switching_cost: &[f64], cfg: &SwitchingConfig) -> Result<SwitchingResult>
+   where R: Fn(usize, usize) -> f64, T: Fn(usize) -> f64;
+
+   pub fn solve_pontryagin_lqr(cfg: &PontryaginConfig) -> Result<PontryaginResult>;
+   pub fn optimal_two_sided_intensities(cfg: &TwoSidedConfig, delta_v_plus: f64, delta_v_minus: f64) -> Result<TwoSidedResult>;

docs/source/index.rst

@@ -51,6 +51,19 @@ Optimiz-rs provides blazingly fast, production-ready implementations of advanced
    algorithms/volterra
    algorithms/signatures
 
+.. toctree::
+   :maxdepth: 2
+   :caption: v2.0 Generic Stochastic Control & PDE
+
+   algorithms/bsde
+   algorithms/pde
+   algorithms/stochastic_control
+   algorithms/quadratic_impact_control
+   algorithms/mckean_vlasov
+   algorithms/agent_based
+   algorithms/robust_drift
+   algorithms/generative_calibration_hooks
+
 .. toctree::
    :maxdepth: 2
    :caption: API Reference