Add make docs commands for html build

docs/Makefile

@@ -0,0 +1,19 @@
+# Minimal makefile for Sphinx documentation
+SHELL := /bin/sh
+
+SPHINXOPTS ?=
+SPHINXBUILD ?= sphinx-build
+SOURCEDIR = source
+BUILDDIR = build
+
+.PHONY: help clean html
+
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS)
+
+clean:
+	rm -rf "$(BUILDDIR)"
+
+html:
+	$(SPHINXBUILD) -b html "$(SOURCEDIR)" "$(BUILDDIR)/html" $(SPHINXOPTS)
+	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."

docs/source/algorithms/differential_evolution.md

@@ -76,22 +76,38 @@ print(f"Best solution: {best_x}")
 
 ## Advanced Features
 
-### Adaptive jDE
-
-Enable self-adaptive F and CR parameters:
+### Adaptive control (jDE, SHADE-ready)
 
 ```python
 de = DifferentialEvolution(
     bounds=[(-5, 5)] * 20,
-    adaptive=True,  # Enable jDE
-    tau_F=0.1,      # F adaptation rate
-    tau_CR=0.1      # CR adaptation rate
+    adaptive=True,   # jDE by default
+    tau_F=0.1,
+    tau_CR=0.1,
 )
 ```
 
-### Constraint Handling
+- jDE is enabled when `adaptive=True` (self-adapts F, CR).
+- SHADE and L-SHADE are implemented in Rust (`shade.rs`) and ready to be wired into the Python API in an upcoming release; see `SHADE_IMPLEMENTATION.md` for details.
 
-For constrained optimization:
+### Parallel evaluation (Rust-only objectives)
+
+For pure Rust benchmarks or when you avoid Python callbacks, you can turn on data-parallel evaluation (Rayon-based) via the Rust entry point:
+
+```python
+from optimizr import parallel_differential_evolution_rust
+
+best = parallel_differential_evolution_rust(
+    objective_name="rastrigin",  # sphere, rosenbrock, ackley, griewank
+    bounds=[(-5, 5)] * 20,
+    maxiter=500,
+    parallel=True,
+)
+```
+
+This yields 10–100× speedups on multi-core for built-in objectives (no GIL contention).
+
+### Constraint handling
 
 ```python
 def constraints(x):
@@ -122,6 +138,10 @@ de = DifferentialEvolution(
    - Smooth, unimodal → `best/1/bin`
    - Multimodal, deceptive → `current-to-best/1/bin`
 
+### Pipeline integrations
+- **Time-series workflows**: couple DE with `timeseries_utils` (rolling Hurst/half-life) to optimize strategy thresholds.
+- **Grid search fallback**: for separable problems, try `grid_search` first; switch to DE when interactions matter.
+
 ## Benchmarks
 
 Performance on standard test functions (10D, 500 iterations):

docs/source/algorithms/hmm.md

@@ -20,6 +20,17 @@ states = model.predict(returns)
 print(np.unique(states, return_counts=True))
 ```
 
+### With time-series helpers
+
+```python
+from optimizr import prepare_for_hmm_py
+
+features = prepare_for_hmm_py(prices, lag_periods=[1, 5, 20])
+hmm = HMM(n_states=3).fit(features, n_iterations=120)
+```
+
+Use rolling Hurst/half-life from `timeseries_utils` as additional features for richer regime classification.
+
 ## Notes
 - Uses Rust backend when available; falls back to Python.
 - `fit` runs Baum-Welch; `predict` runs Viterbi.

docs/source/algorithms/mean_field_games.md

@@ -1,34 +1,50 @@
 # Mean Field Games
 
-Solve 1D Mean Field Games and Mean Field-Type Control problems using the Rust backend.
+Solve 1D Mean Field Games and mean-field–type control problems with the Rust backend. The solver couples a backward HJB equation with a forward Fokker–Planck equation using a fixed-point loop and implicit diffusion for stability.
 
-## Key Concepts
-- Coupled HJB–Fokker-Planck PDE system
-- Density evolution over time/space
-- Congestion and noise parameters control stability
+## What this module provides
+- **Rust solver with PyO3 bindings**: `solve_mfg_1d_rust` and `MFGConfig` exposed to Python.
+- **Stable numerics**: upwind transport + implicit diffusion, mass renormalization each iteration.
+- **Performance**: ~0.4 s for a 100×100 grid on laptop-class CPUs (measured in the tutorial notebook).
+- **Coverage**: Congestion term, relaxation `alpha`, configurable domain and viscosity `nu`.
 
 ## Usage
 
 ```python
+import numpy as np
 from optimizr import MFGConfig, solve_mfg_1d_rust
 
+x = np.linspace(0, 1, 100)
+m0 = np.exp(-50 * (x - 0.3) ** 2)
+m0 /= np.trapz(m0, x)
+
+u_terminal = 0.5 * (x - 0.7) ** 2
 config = MFGConfig(
-    nx=64,
-    nt=40,
-    x_min=-3.0,
-    x_max=3.0,
+    nx=100,
+    nt=100,
+    x_min=0.0,
+    x_max=1.0,
     T=1.0,
-    epsilon=0.1,
-    kappa=1.0,
+    nu=0.01,
+    max_iter=50,
+    tol=1e-5,
+    alpha=0.5,
 )
 
-solution = solve_mfg_1d_rust(config)
-print("Converged:", solution.converged)
-print("Value function grid:", solution.value_function.shape)
-print("Density grid:", solution.density.shape)
+u, m, iters = solve_mfg_1d_rust(m0, u_terminal, config, lambda_congestion=0.5)
+print(f"converged in {iters} iterations: u{u.shape}, m{m.shape}")
 ```
 
-## Tips
-- Increase `nx`/`nt` for smoother solutions; expect higher compute.
-- Reduce `epsilon` for sharper dynamics; increase if unstable.
-- Use `kappa` to control congestion cost.
+## What to monitor
+- Residuals: track `||m^{k+1}-m^k||_1` and `||u^{k+1}-u^k||_inf`; stop when both flatten.
+- Mass conservation: integrate `m` after each iteration; values close to 1.0 indicate stable transport.
+- CFL sanity: if oscillations appear, reduce `dt` (increase `nt`) or raise `nu` slightly.
+
+## Practical tips
+- Grid resolution: start with `nx=64, nt=40`; move to 100×100 for publication-quality plots.
+- Congestion: increase `lambda_congestion` to avoid density spikes; decrease for freer flow.
+- Relaxation: `alpha=0.5` is a stable default; lower if the fixed-point loop jitters.
+
+## Notebook and audit
+- Full walkthrough: `examples/notebooks/mean_field_games_tutorial.ipynb` (all cells validated).
+- Audit notes: the notebook renders convergence plots, 3D density/value surfaces, and time-slice snapshots; runs cleanly with the Rust backend (see `docs/MFG_TUTORIAL_COMPLETE.md`).

docs/source/algorithms/optimal_control.md

@@ -1,28 +1,94 @@
 # Optimal Control
 
-High-level bindings for HJB-style optimal control and Kalman filtering utilities.
+Hamilton–Jacobi–Bellman (HJB) solvers, regime-switching thresholds, OU parameter estimation, and Kalman filtering utilities backed by Rust.
 
-## Kalman Filter (sensor fusion)
+## HJB switching boundaries (OU process)
+
+```python
+from optimizr import solve_hjb_py, solve_hjb_full_py
+
+lower, upper, residual, iters = solve_hjb_py(
+    kappa=3.0,
+    theta=0.0,
+    sigma=0.2,
+    rho=0.04,
+    transaction_cost=0.001,
+    n_points=400,
+    max_iter=4000,
+    tolerance=1e-7,
+    n_std=5.0,
+)
+print(f"bounds=({lower:.3f}, {upper:.3f}), residual={residual:.2e}, iters={iters}")
+```
+
+- Model: $dX_t = \kappa(\theta - X_t)\,dt + \sigma\,dW_t$ with quadratic transaction costs.
+- Output: optimal buy/sell thresholds; `solve_hjb_full_py` also returns $V, V_x, V_{xx}$ for diagnostics.
+- Diagnostics: plot $V_x$ for smoothness near thresholds; monitor `residual` and increase `max_iter` if not converged.
+
+## Backtesting optimal switching
+
+```python
+from optimizr import backtest_optimal_switching_py
+
+metrics = backtest_optimal_switching_py(
+    spread=spread,
+    lower_bound=lower,
+    upper_bound=upper,
+    transaction_cost=0.001,
+)
+(
+    total_return,
+    sharpe,
+    max_dd,
+    n_trades,
+    win_rate,
+    pnl_path,
+) = metrics
+```
+
+Inspect `win_rate` vs `max_dd` to tune aggressiveness; combine with HMM regimes for state-aware controls.
+
+## OU parameter estimation
 
 ```python
 import numpy as np
-from optimizr import maths_toolkit
+from optimizr import estimate_ou_params_py
 
-# maths_toolkit is provided by the Rust extension (_core)
-F = np.eye(2)  # state transition
-H = np.eye(2)  # observation
-Q = 0.01 * np.eye(2)
-R = 0.1 * np.eye(2)
-
-if maths_toolkit is not None:
-    kf = maths_toolkit.init_kalman_filter(F.tolist(), H.tolist(), Q.tolist(), R.tolist())
-    state = maths_toolkit.kalman_predict(kf, [0.0, 0.0])
-    print(state)
-else:
-    print("Rust backend not available; install with `pip install .`.")
+spread = np.random.randn(10_000)
+kappa, theta, sigma, half_life = estimate_ou_params_py(spread, dt=1/252)
 ```
 
-## Notes
-- Rust backend (`optimizr._core`) must be present for control utilities.
-- The API is thin and intentionally low-level; matrices are passed as lists.
-- For 1D Mean Field Games, use the dedicated guide in `mean_field_games.md`.
+Method-of-moments / MLE fit returns $(\kappa, \theta, \sigma, \text{half-life})$. Use a few thousand samples for stability; winsorize heavy tails if needed.
+
+## Kalman filtering (linear, EKF, UKF)
+
+```python
+import numpy as np
+from optimizr import LinearKalmanFilter
+
+F = [[1.0, 1.0], [0.0, 1.0]]
+H = [[1.0, 0.0]]
+Q = [[1e-4, 0.0], [0.0, 1e-4]]
+R = [[1e-2]]
+
+kf = LinearKalmanFilter(
+    f_matrix=F,
+    h_matrix=H,
+    q_matrix=Q,
+    r_matrix=R,
+    initial_state=[0.0, 0.0],
+    initial_covariance=[[1.0, 0.0], [0.0, 1.0]],
+)
+
+kf.predict(control=[0.0, 0.0])
+kf.update(observation=[1.2])
+state = kf.get_state()
+```
+
+- Interfaces: `LinearKalmanFilter`, `UnscentedKalmanFilter`, and `KalmanState` for batch `filter` and smoothing.
+- Concept: prediction (dynamics prior) + correction (measurement residual); RTS smoother refines past states.
+
+## Practical notes
+- Rust backend (`optimizr._core`) must be present for control utilities; install from source if wheels are unavailable.
+- Grids: for HJB, `n_points≈400` is stable; widen `n_std` for volatile spreads.
+- Combine with Mean Field Games: see `mean_field_games.md` for population dynamics; use Kalman estimates as control inputs if needed.

docs/source/algorithms/risk_metrics.md

@@ -1,8 +1,8 @@
 # Risk Metrics
 
-Time-series utilities for risk analysis and mean-reversion signals.
+Time-series utilities for risk analysis, mean-reversion detection, and bootstrapped P&L distributions.
 
-## Quick Start
+## Quick start
 
 ```python
 import numpy as np
@@ -10,17 +10,27 @@ from optimizr import (
     hurst_exponent_py,
     estimate_half_life_py,
     bootstrap_returns_py,
+    compute_risk_metrics_py,
 )
 
 returns = np.random.randn(2000) * 0.01
 print("Hurst:", hurst_exponent_py(returns))
 print("Half-life:", estimate_half_life_py(returns))
 
+metrics = compute_risk_metrics_py(returns)
+print(metrics)  # mean, std, skew, kurtosis, sharpe
+
 bootstrapped = bootstrap_returns_py(returns, n_samples=1000)
-print("Bootstrap sample shape:", len(bootstrapped))
+print("Bootstrap samples:", len(bootstrapped))
 ```
 
-## Notes
-- Input arrays should be 1D NumPy arrays of returns.
-- Half-life is useful for calibrating mean-reversion strategies.
-- Bootstrap utilities help estimate drawdown and VaR distributions.
+## Rolling and integration helpers
+
+- Use `rolling_hurst_exponent_py` and `rolling_half_life_py` (from `timeseries_utils`) for sliding-window diagnostics on trading pairs.
+- Combine with HMM: feed rolling statistics as features for regime detection.
+- Pair with DE/Grid search: optimize strategy thresholds while computing half-life inside the objective.
+
+## Practical guidance
+- Input should be 1D NumPy arrays of returns; winsorize extreme tails before estimating Hurst/half-life for stability.
+- Half-life helps size holding periods for mean-reversion trades; revisit whenever volatility regime changes.
+- Bootstrap outputs can feed VaR/ES estimates; increase `n_samples` for tighter confidence bands.

docs/source/algorithms/sparse_optimization.md

@@ -1,6 +1,6 @@
 # Sparse Optimization
 
-Sparse PCA and Elastic Net utilities with Rust speed.
+Sparse PCA, Elastic Net, and Box–Tao decomposition with Rust speed.
 
 ## Sparse PCA
 
@@ -13,6 +13,9 @@ components = sparse_pca_py(X, n_components=5, l1_ratio=0.15)
 print(components.shape)  # (5, 20)
 ```
 
+- Output: component matrix `(n_components, n_features)`; rows are sparse loadings.
+- Tuning: increase `l1_ratio` for harder sparsity; decrease to retain variance.
+
 ## Elastic Net
 
 ```python
@@ -25,7 +28,19 @@ coeffs = elastic_net_py(X, y, l1_ratio=0.3, alpha=0.01)
 print(coeffs)
 ```
 
-## Notes
-- Inputs should be NumPy arrays; data is copied to Rust.
-- `l1_ratio` balances sparsity vs ridge penalty.
-- Standardize features before calling for stable solutions.
+- Handles collinearity better than pure Lasso; use for factor shrinkage.
+- Sweep `alpha` on a log scale (e.g., $10^{-3}$ to $10^{-1}$) and pick via validation.
+
+## Box–Tao decomposition
+
+```python
+from optimizr import box_tao_decomposition_py
+solution = box_tao_decomposition_py(X)
+```
+
+Useful for constrained sparse decomposition problems; the Rust backend keeps iterations fast.
+
+## Practical notes
+- Inputs must be NumPy arrays; standardize features for stable conditioning.
+- For high dimensional data, start with fewer components/features to avoid over-regularization.
+- Combine with risk metrics: use sparse loadings to build interpretable factors, then evaluate with `compute_risk_metrics_py`.

docs/source/api/differential_evolution.md

@@ -8,20 +8,38 @@ best_x, best_fx = differential_evolution(
     bounds,
     popsize=15,
     maxiter=1000,
-    f=None,
-    cr=None,
-    strategy="rand1",
+    f=None,                 # mutation factor (auto if None)
+    cr=None,                # crossover rate (auto if None)
+    strategy="rand1",      # rand1, best1, currenttobest1, rand2, best2
     seed=None,
     tol=1e-6,
     atol=1e-8,
-    track_history=False,
-    parallel=False,
-    adaptive=False,
+    track_history=False,    # keep per-iter best
+    parallel=False,         # Python callbacks stay sequential; see Rust path below
+    adaptive=False,         # jDE when True
     constraint_penalty=1000.0,
 )
 ```
 
 - `objective_fn`: callable `f(x: np.ndarray) -> float`
 - `bounds`: list of `(min, max)` tuples
-- Strategies: `rand1`, `best1`, `currenttobest1`, `rand2`, `best2`
 - Returns `(best_x: np.ndarray, best_fx: float)`
+
+## Parallel Rust entry point
+
+For built-in benchmark objectives (no Python callbacks), use the Rust-native path with Rayon:
+
+```python
+from optimizr import parallel_differential_evolution_rust
+
+result = parallel_differential_evolution_rust(
+    objective_name="rastrigin",  # sphere, rosenbrock, ackley, griewank
+    bounds=[(-5, 5)] * 20,
+    maxiter=500,
+    parallel=True,
+)
+```
+
+## Notes
+- Adaptive control uses jDE in the current Python API; SHADE/L-SHADE live in Rust and will surface in a future release.
+- Use `track_history=True` to export convergence curves for benchmarking.

docs/source/index.rst

@@ -72,6 +72,8 @@ Features
 - **Sparse Optimization**: Sparse PCA, Box-Tao decomposition, Elastic Net, ADMM
 - **Risk Metrics**: Hurst exponent, half-life estimation, time series analysis
 - **Information Theory**: Mutual information, Shannon entropy, feature selection
+- **Time-Series Helpers**: Rolling Hurst/half-life, feature prep for HMM, lagged feature builders
+- **Parallelization**: Rust-native population evaluation with Rayon for built-in objectives
 
 🚀 **Performance:**
 

docs/source/installation.md

@@ -83,6 +83,9 @@ Requires GCC or Clang:
 # Ubuntu/Debian
 sudo apt-get install build-essential
 
+# If you see OpenBLAS link errors during wheels/docs build
+sudo apt-get install libopenblas-dev
+
 # Fedora/RHEL
 sudo dnf install gcc gcc-c++
 ```
@@ -109,3 +112,7 @@ rustup update stable
 ```bash
 maturin develop --release -i python3.10  # Replace with your Python version
 ```
+
+**Issue**: BLAS/LAPACK linkage errors on Linux
+
+**Solution**: Install OpenBLAS headers (see Linux section above) and rebuild with `maturin develop --release`.

docs/source/quickstart.md

@@ -93,3 +93,9 @@ print(f"Converged: {solution.converged}")
 - See [Getting Started](getting-started.md) for environment setup and verification.
 - Browse [Examples](examples.md) for code snippets per optimizer.
 - Deep dive into algorithms in [Algorithms](algorithms/differential_evolution.md).
+
+## Notebook status and reproducibility
+
+- Audit (2025-01-04): 6/7 notebooks execute cleanly; `03_optimal_control_tutorial.ipynb` is theory-only by design.
+- Fully validated: `01_hmm_tutorial`, `02_mcmc_tutorial`, `04_real_world_applications`, `05_performance_benchmarks`, `mean_field_games_tutorial`.
+- Differential Evolution tutorial works with current API; enable `track_history=True` to capture convergence curves during runs.