feat(v1.1.x): PyO3 bindings + executed companion notebooks for 5 new groups

Adds Python bindings (behind feature='python-bindings') for graph,
risk_measures, topology, volterra, signatures.

Companion notebooks under examples/notebooks/:
- 05_graph.ipynb           (Laplacians + spectral clustering)
- 06_risk_measures.ipynb   (VaR / CVaR + simplex projection)
- 07_topology.ipynb        (Vietoris-Rips + persistent homology)
- 08_volterra.ipynb        (fractional ODE, Markovian lift, Volterra,
                            Fourier inversion)
- 09_signatures.ipynb      (path / log / random / kernel signatures)

All notebooks executed end-to-end against analytic ground truth
(closed-form solutions, Mittag-Leffler, exp(-t), unit-circle homology,
identical-path signature kernel).

Built and validated via: maturin develop --release --features python-bindings.

Workflow generated by 5 parallel optimizRs subagents (.github/agents/).

src/graph/mod.rs

@@ -12,3 +12,6 @@ pub mod spectral_clustering;
 
 pub use laplacian::{combinatorial_laplacian, normalised_laplacian, random_walk_laplacian, LaplacianKind};
 pub use spectral_clustering::{spectral_cluster, SpectralClusterResult};
+
+#[cfg(feature = "python-bindings")]
+pub mod python_bindings;

src/graph/python_bindings.rs

@@ -0,0 +1,115 @@
+//! Python bindings for the `graph` module.
+//!
+//! Exposes graph Laplacian operators and the Ng--Jordan--Weiss
+//! spectral clustering routine.
+
+use ndarray::Array2;
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+use super::laplacian::{combinatorial_laplacian, normalised_laplacian, random_walk_laplacian};
+use super::spectral_clustering::spectral_cluster;
+
+fn vec_of_vec_to_array2(w: &[Vec<f64>]) -> PyResult<Array2<f64>> {
+    let n = w.len();
+    if n == 0 {
+        return Err(PyValueError::new_err(
+            "weight matrix must be non-empty",
+        ));
+    }
+    let m = w[0].len();
+    if m != n {
+        return Err(PyValueError::new_err(
+            "weight matrix must be square",
+        ));
+    }
+    let mut a = Array2::<f64>::zeros((n, n));
+    for (i, row) in w.iter().enumerate() {
+        if row.len() != n {
+            return Err(PyValueError::new_err(
+                "weight matrix must be square",
+            ));
+        }
+        for (j, &v) in row.iter().enumerate() {
+            a[[i, j]] = v;
+        }
+    }
+    Ok(a)
+}
+
+fn array2_to_vec_of_vec(a: &Array2<f64>) -> Vec<Vec<f64>> {
+    let n = a.nrows();
+    let m = a.ncols();
+    let mut out = Vec::with_capacity(n);
+    for i in 0..n {
+        let mut row = Vec::with_capacity(m);
+        for j in 0..m {
+            row.push(a[[i, j]]);
+        }
+        out.push(row);
+    }
+    out
+}
+
+/// Combinatorial Laplacian `L = D - W`.
+#[pyfunction]
+#[pyo3(signature = (w))]
+fn combinatorial_laplacian_py(w: Vec<Vec<f64>>) -> PyResult<Vec<Vec<f64>>> {
+    let arr = vec_of_vec_to_array2(&w)?;
+    let l = combinatorial_laplacian(arr.view())
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    Ok(array2_to_vec_of_vec(&l))
+}
+
+/// Symmetric normalised Laplacian `L_sym = I - D^{-1/2} W D^{-1/2}`.
+#[pyfunction]
+#[pyo3(signature = (w))]
+fn normalised_laplacian_py(w: Vec<Vec<f64>>) -> PyResult<Vec<Vec<f64>>> {
+    let arr = vec_of_vec_to_array2(&w)?;
+    let l = normalised_laplacian(arr.view())
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    Ok(array2_to_vec_of_vec(&l))
+}
+
+/// Random-walk Laplacian `L_rw = I - D^{-1} W`.
+#[pyfunction]
+#[pyo3(signature = (w))]
+fn random_walk_laplacian_py(w: Vec<Vec<f64>>) -> PyResult<Vec<Vec<f64>>> {
+    let arr = vec_of_vec_to_array2(&w)?;
+    let l = random_walk_laplacian(arr.view())
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    Ok(array2_to_vec_of_vec(&l))
+}
+
+/// Spectral clustering (Ng--Jordan--Weiss) on a non-negative symmetric
+/// similarity matrix.
+///
+/// Returns a dict with keys `labels`, `eigenvalues`, `fiedler_value`.
+#[pyfunction]
+#[pyo3(signature = (w, k, n_kmeans_iter=100, seed=0))]
+fn spectral_cluster_py(
+    py: Python<'_>,
+    w: Vec<Vec<f64>>,
+    k: usize,
+    n_kmeans_iter: usize,
+    seed: u64,
+) -> PyResult<PyObject> {
+    let arr = vec_of_vec_to_array2(&w)?;
+    let result = spectral_cluster(arr.view(), k, n_kmeans_iter, seed)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("labels", result.labels.clone())?;
+    dict.set_item("eigenvalues", result.eigenvalues.clone())?;
+    dict.set_item("fiedler_value", result.fiedler_value)?;
+    Ok(dict.into())
+}
+
+/// Register all graph functions with the Python module.
+pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(combinatorial_laplacian_py, m)?)?;
+    m.add_function(wrap_pyfunction!(normalised_laplacian_py, m)?)?;
+    m.add_function(wrap_pyfunction!(random_walk_laplacian_py, m)?)?;
+    m.add_function(wrap_pyfunction!(spectral_cluster_py, m)?)?;
+    Ok(())
+}

src/lib.rs

@@ -133,5 +133,12 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     // Portfolio Optimization functions (CARA, Mean-Variance, ERC)
     portfolio_optimization::python_bindings::register_python_functions(m)?;
 
+    // ===== v1.1.0 additive bindings =====
+    graph::python_bindings::register_python_functions(m)?;
+    risk_measures::python_bindings::register_python_functions(m)?;
+    topology::python_bindings::register_python_functions(m)?;
+    volterra::python_bindings::register_python_functions(m)?;
+    signatures::python_bindings::register_python_functions(m)?;
+
     Ok(())
 }

src/risk_measures/mod.rs

@@ -15,6 +15,9 @@ pub mod cvar;
 pub use cvar::{cvar_value, minimize_cvar, CVaRConfig, CVaRResult};
 pub use var::{historical_var, parametric_var};
 
+#[cfg(feature = "python-bindings")]
+pub mod python_bindings;
+
 #[inline]
 pub(crate) fn check_alpha(alpha: f64) -> Result<()> {
     if !(0.0 < alpha && alpha < 1.0) {

src/risk_measures/python_bindings.rs

@@ -0,0 +1,92 @@
+//! Python bindings for empirical risk measures.
+//!
+//! Exposes Value-at-Risk and Conditional Value-at-Risk estimators
+//! together with the simplex-constrained CVaR minimiser.
+
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+use ndarray::Array2;
+
+use super::cvar::{cvar_value, minimize_cvar, CVaRConfig};
+use super::var::{historical_var, parametric_var};
+
+/// Empirical Value-at-Risk at confidence level `alpha`.
+#[pyfunction]
+#[pyo3(signature = (losses, alpha=0.95))]
+fn historical_var_py(losses: Vec<f64>, alpha: f64) -> PyResult<f64> {
+    historical_var(&losses, alpha).map_err(|e| PyValueError::new_err(format!("{}", e)))
+}
+
+/// Closed-form Gaussian Value-at-Risk `mu + sigma * Phi^{-1}(alpha)`.
+#[pyfunction]
+#[pyo3(signature = (mu, sigma, alpha=0.95))]
+fn parametric_var_py(mu: f64, sigma: f64, alpha: f64) -> PyResult<f64> {
+    parametric_var(mu, sigma, alpha).map_err(|e| PyValueError::new_err(format!("{}", e)))
+}
+
+/// Empirical Conditional Value-at-Risk at confidence level `alpha`.
+#[pyfunction]
+#[pyo3(signature = (losses, alpha=0.95))]
+fn cvar_value_py(losses: Vec<f64>, alpha: f64) -> PyResult<f64> {
+    cvar_value(&losses, alpha).map_err(|e| PyValueError::new_err(format!("{}", e)))
+}
+
+/// Minimise empirical CVaR of `L(w) = -<r^{(s)}, w>` over the unit
+/// simplex. `samples` has shape `(S, d)` (S samples, d decision
+/// components).
+#[pyfunction]
+#[pyo3(signature = (samples, alpha=0.95, n_iter=5000, step_size=0.01, tol=1e-8))]
+fn minimize_cvar_py(
+    py: Python<'_>,
+    samples: Vec<Vec<f64>>,
+    alpha: f64,
+    n_iter: usize,
+    step_size: f64,
+    tol: f64,
+) -> PyResult<PyObject> {
+    let s = samples.len();
+    if s == 0 {
+        return Err(PyValueError::new_err("samples must be non-empty"));
+    }
+    let d = samples[0].len();
+    if d == 0 {
+        return Err(PyValueError::new_err("samples rows must be non-empty"));
+    }
+    let mut flat = Vec::with_capacity(s * d);
+    for row in &samples {
+        if row.len() != d {
+            return Err(PyValueError::new_err(
+                "all sample rows must share the same length",
+            ));
+        }
+        flat.extend_from_slice(row);
+    }
+    let arr = Array2::from_shape_vec((s, d), flat)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+
+    let cfg = CVaRConfig {
+        alpha,
+        n_iter,
+        step_size,
+        tol,
+    };
+    let result = minimize_cvar(arr.view(), &cfg)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("w", result.w.to_vec())?;
+    dict.set_item("zeta", result.zeta)?;
+    dict.set_item("cvar", result.cvar)?;
+    dict.set_item("iterations", result.iterations)?;
+    Ok(dict.into())
+}
+
+/// Register all risk-measure functions with the Python module.
+pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(historical_var_py, m)?)?;
+    m.add_function(wrap_pyfunction!(parametric_var_py, m)?)?;
+    m.add_function(wrap_pyfunction!(cvar_value_py, m)?)?;
+    m.add_function(wrap_pyfunction!(minimize_cvar_py, m)?)?;
+    Ok(())
+}

src/signatures/mod.rs

@@ -14,3 +14,6 @@ pub mod random_signature;
 pub mod signature_kernel;
 pub mod utils;
 
+#[cfg(feature = "python-bindings")]
+pub mod python_bindings;
+

src/signatures/python_bindings.rs

@@ -0,0 +1,168 @@
+//! Python bindings for the `signatures` module.
+//!
+//! Exposes truncated path signatures, log-signatures, random reservoir
+//! projections, the Salvi--Cass--Lyons signature kernel, and the
+//! shuffle product / Chen concatenation utilities.
+
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+use super::log_signature::log_signature as rs_log_signature;
+use super::path_signature::{path_signature as rs_path_signature, TruncatedSignature};
+use super::random_signature::{
+    random_signature as rs_random_signature, RandomSignatureConfig,
+};
+use super::signature_kernel::signature_kernel as rs_signature_kernel;
+use super::utils::{concatenate_signatures as rs_concatenate, shuffle_product as rs_shuffle};
+
+fn signature_to_dict(py: Python<'_>, sig: &TruncatedSignature) -> PyResult<PyObject> {
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("channels", sig.channels)?;
+    dict.set_item("level", sig.level)?;
+    dict.set_item("tensors", sig.tensors.clone())?;
+    Ok(dict.into())
+}
+
+fn dict_to_signature(channels: usize, level: usize, tensors: Vec<Vec<f64>>) -> PyResult<TruncatedSignature> {
+    if tensors.len() != level + 1 {
+        return Err(PyValueError::new_err(format!(
+            "tensors length {} does not match level + 1 = {}",
+            tensors.len(),
+            level + 1
+        )));
+    }
+    for (k, t) in tensors.iter().enumerate() {
+        let expected = channels.pow(k as u32);
+        if t.len() != expected {
+            return Err(PyValueError::new_err(format!(
+                "tensors[{}] has length {}, expected channels^k = {}",
+                k,
+                t.len(),
+                expected
+            )));
+        }
+    }
+    Ok(TruncatedSignature {
+        channels,
+        level,
+        tensors,
+    })
+}
+
+/// Truncated tensor signature of a piecewise-linear multivariate path.
+#[pyfunction]
+#[pyo3(signature = (path, level))]
+fn path_signature(py: Python<'_>, path: Vec<Vec<f64>>, level: usize) -> PyResult<PyObject> {
+    let sig = rs_path_signature(&path, level)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    signature_to_dict(py, &sig)
+}
+
+/// Truncated tensor log-signature of a piecewise-linear multivariate path.
+#[pyfunction]
+#[pyo3(signature = (path, level))]
+fn path_log_signature(py: Python<'_>, path: Vec<Vec<f64>>, level: usize) -> PyResult<PyObject> {
+    let sig = rs_path_signature(&path, level)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    let log = rs_log_signature(&sig).map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("channels", log.channels)?;
+    dict.set_item("level", log.level)?;
+    dict.set_item("tensors", log.tensors.clone())?;
+    Ok(dict.into())
+}
+
+/// Random reservoir projection of the path signature.
+#[pyfunction]
+#[pyo3(signature = (path, reservoir_dim=32, seed=0, variance=1.0))]
+fn random_signature(
+    py: Python<'_>,
+    path: Vec<Vec<f64>>,
+    reservoir_dim: usize,
+    seed: u64,
+    variance: f64,
+) -> PyResult<PyObject> {
+    let cfg = RandomSignatureConfig {
+        reservoir_dim,
+        seed,
+        variance,
+    };
+    let res = rs_random_signature(&path, &cfg)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("trajectory", res.trajectory.clone())?;
+    Ok(dict.into())
+}
+
+/// Salvi--Cass--Lyons signature kernel between two multivariate paths.
+#[pyfunction]
+#[pyo3(signature = (x, y))]
+fn signature_kernel(
+    py: Python<'_>,
+    x: Vec<Vec<f64>>,
+    y: Vec<Vec<f64>>,
+) -> PyResult<PyObject> {
+    let res = rs_signature_kernel(&x, &y)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    let dict = pyo3::types::PyDict::new_bound(py);
+    dict.set_item("value", res.value)?;
+    dict.set_item("grid", res.grid.clone())?;
+    Ok(dict.into())
+}
+
+/// Shuffle product of two words `u`, `v` over `{0, ..., d-1}`.
+///
+/// Returned as a list of `(word, multiplicity)` tuples to remain
+/// hashable-key agnostic on the Python side.
+#[pyfunction]
+#[pyo3(signature = (u, v))]
+fn shuffle_product(
+    py: Python<'_>,
+    u: Vec<usize>,
+    v: Vec<usize>,
+) -> PyResult<PyObject> {
+    let map = rs_shuffle(&u, &v);
+    let list = pyo3::types::PyList::empty_bound(py);
+    for (word, mult) in map.into_iter() {
+        let tup = pyo3::types::PyTuple::new_bound(py, &[word.into_py(py), mult.into_py(py)]);
+        list.append(tup)?;
+    }
+    Ok(list.into())
+}
+
+/// Concatenate two truncated signatures via Chen's identity.
+///
+/// Inputs are passed as `(channels, level, tensors)` triples matching
+/// the dict layout returned by :func:`path_signature`.
+#[pyfunction]
+#[pyo3(signature = (a_channels, a_level, a_tensors, b_channels, b_level, b_tensors))]
+fn concatenate_signatures(
+    py: Python<'_>,
+    a_channels: usize,
+    a_level: usize,
+    a_tensors: Vec<Vec<f64>>,
+    b_channels: usize,
+    b_level: usize,
+    b_tensors: Vec<Vec<f64>>,
+) -> PyResult<PyObject> {
+    if a_channels != b_channels || a_level != b_level {
+        return Err(PyValueError::new_err(
+            "signatures must share channels and level",
+        ));
+    }
+    let a = dict_to_signature(a_channels, a_level, a_tensors)?;
+    let b = dict_to_signature(b_channels, b_level, b_tensors)?;
+    let out = rs_concatenate(&a, &b).map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    signature_to_dict(py, &out)
+}
+
+/// Register all signatures functions with the Python module.
+pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(path_signature, m)?)?;
+    m.add_function(wrap_pyfunction!(path_log_signature, m)?)?;
+    m.add_function(wrap_pyfunction!(random_signature, m)?)?;
+    m.add_function(wrap_pyfunction!(signature_kernel, m)?)?;
+    m.add_function(wrap_pyfunction!(shuffle_product, m)?)?;
+    m.add_function(wrap_pyfunction!(concatenate_signatures, m)?)?;
+    Ok(())
+}

src/topology/mod.rs

@@ -14,3 +14,6 @@ pub use bottleneck::bottleneck_distance;
 pub use persistent_homology::{
     persistent_homology, vietoris_rips_filtration, PersistenceDiagram, PersistencePair,
 };
+
+#[cfg(feature = "python-bindings")]
+pub mod python_bindings;

src/topology/python_bindings.rs

@@ -0,0 +1,105 @@
+//! Python bindings for topological data analysis.
+//!
+//! Exposes Vietoris--Rips filtration construction, persistent homology,
+//! and bottleneck distance between persistence diagrams.
+
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+
+use super::bottleneck::bottleneck_distance as rs_bottleneck_distance;
+use super::persistent_homology::{
+    persistent_homology as rs_persistent_homology,
+    vietoris_rips_filtration as rs_vietoris_rips_filtration, PersistencePair,
+};
+
+fn diagram_to_pylist(py: Python<'_>, pairs: &[PersistencePair]) -> PyResult<PyObject> {
+    let list = pyo3::types::PyList::empty_bound(py);
+    for p in pairs {
+        let d = pyo3::types::PyDict::new_bound(py);
+        d.set_item("dim", p.dim)?;
+        d.set_item("birth", p.birth)?;
+        d.set_item("death", p.death)?;
+        list.append(d)?;
+    }
+    Ok(list.into())
+}
+
+fn pylist_to_diagram(pairs: Vec<Bound<'_, pyo3::types::PyDict>>) -> PyResult<Vec<PersistencePair>> {
+    let mut out = Vec::with_capacity(pairs.len());
+    for d in pairs {
+        let dim: usize = d
+            .get_item("dim")?
+            .ok_or_else(|| PyValueError::new_err("missing key 'dim'"))?
+            .extract()?;
+        let birth: f64 = d
+            .get_item("birth")?
+            .ok_or_else(|| PyValueError::new_err("missing key 'birth'"))?
+            .extract()?;
+        let death: f64 = d
+            .get_item("death")?
+            .ok_or_else(|| PyValueError::new_err("missing key 'death'"))?
+            .extract()?;
+        out.push(PersistencePair { dim, birth, death });
+    }
+    Ok(out)
+}
+
+/// Build the Vietoris--Rips filtration up to ``max_dim`` and scale
+/// ``max_eps``. Returns a list of simplices as dicts
+/// ``{"vertices": [..], "filtration": float, "dim": int}``.
+#[pyfunction]
+#[pyo3(signature = (points, max_dim=1, max_eps=1.0))]
+fn vietoris_rips_filtration(
+    py: Python<'_>,
+    points: Vec<Vec<f64>>,
+    max_dim: usize,
+    max_eps: f64,
+) -> PyResult<PyObject> {
+    let simplices = rs_vietoris_rips_filtration(&points, max_dim, max_eps)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    let list = pyo3::types::PyList::empty_bound(py);
+    for s in simplices {
+        let d = pyo3::types::PyDict::new_bound(py);
+        d.set_item("vertices", s.vertices.clone())?;
+        d.set_item("filtration", s.filtration)?;
+        d.set_item("dim", s.dim())?;
+        list.append(d)?;
+    }
+    Ok(list.into())
+}
+
+/// Compute the Vietoris--Rips persistence diagram. Returns a list of
+/// dicts ``[{"dim": int, "birth": float, "death": float}, ...]``.
+#[pyfunction]
+#[pyo3(signature = (points, max_dim=1, max_eps=1.0))]
+fn persistent_homology(
+    py: Python<'_>,
+    points: Vec<Vec<f64>>,
+    max_dim: usize,
+    max_eps: f64,
+) -> PyResult<PyObject> {
+    let diag = rs_persistent_homology(&points, max_dim, max_eps)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    diagram_to_pylist(py, &diag.pairs)
+}
+
+/// Bottleneck distance between two persistence diagrams. Each diagram
+/// is a list of dicts ``{"dim": int, "birth": float, "death": float}``.
+#[pyfunction]
+#[pyo3(signature = (diagram_a, diagram_b))]
+fn bottleneck_distance(
+    diagram_a: Vec<Bound<'_, pyo3::types::PyDict>>,
+    diagram_b: Vec<Bound<'_, pyo3::types::PyDict>>,
+) -> PyResult<f64> {
+    let a = pylist_to_diagram(diagram_a)?;
+    let b = pylist_to_diagram(diagram_b)?;
+    rs_bottleneck_distance(&a, &b).map_err(|e| PyValueError::new_err(format!("{}", e)))
+}
+
+/// Register all topology functions with the Python module.
+pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(vietoris_rips_filtration, m)?)?;
+    m.add_function(wrap_pyfunction!(persistent_homology, m)?)?;
+    m.add_function(wrap_pyfunction!(bottleneck_distance, m)?)?;
+    Ok(())
+}

src/volterra/mod.rs

@@ -14,3 +14,6 @@ pub mod fractional_riccati;
 pub mod markovian_lift;
 pub mod volterra_solver;
 
+#[cfg(feature = "python-bindings")]
+pub mod python_bindings;
+

src/volterra/python_bindings.rs

@@ -0,0 +1,224 @@
+//! Python bindings for the `volterra` module.
+//!
+//! Exposes the four core primitives:
+//!
+//! * `solve_fractional_ode`  -- Caputo fractional ODE Adams scheme.
+//! * `geometric_grid_lift`   -- multi-exponential approximation of a kernel.
+//! * `solve_volterra`        -- generic second-kind Volterra equation.
+//! * `fourier_invert`        -- characteristic function -> density.
+//!
+//! Python callables are accepted as `&Bound<'_, PyAny>` and called
+//! through `.call1(...)?.extract::<...>()?`. Errors of type
+//! [`crate::core::OptimizrError`] are mapped to [`PyValueError`].
+//!
+//! The underlying Rust solvers require an immutable [`Fn`] callback,
+//! so the first error raised by the Python callable is captured via a
+//! [`RefCell`] and re-raised after the solver returns.
+
+use std::cell::RefCell;
+
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+use pyo3::types::PyDict;
+
+use super::fourier_inversion::fourier_invert as rs_fourier_invert;
+use super::fractional_riccati::solve_fractional_ode as rs_solve_fractional_ode;
+use super::markovian_lift::geometric_grid_lift as rs_geometric_grid_lift;
+use super::volterra_solver::solve_volterra as rs_solve_volterra;
+
+/// Solve the Caputo fractional ODE `D^alpha h = rhs(t, h)` on `[0, t_horizon]`.
+///
+/// `rhs` is a Python callable `(t: float, h: float) -> float`.
+///
+/// Returns a dict `{"t_grid": [...], "h": [...]}`.
+#[pyfunction]
+#[pyo3(signature = (h0, alpha, t_horizon, n_steps, rhs))]
+fn solve_fractional_ode(
+    py: Python<'_>,
+    h0: f64,
+    alpha: f64,
+    t_horizon: f64,
+    n_steps: usize,
+    rhs: &Bound<'_, PyAny>,
+) -> PyResult<PyObject> {
+    let cb_err: RefCell<Option<PyErr>> = RefCell::new(None);
+    let rust_rhs = |t: f64, h: f64| -> f64 {
+        if cb_err.borrow().is_some() {
+            return 0.0;
+        }
+        match rhs.call1((t, h)).and_then(|v| v.extract::<f64>()) {
+            Ok(v) => v,
+            Err(e) => {
+                *cb_err.borrow_mut() = Some(e);
+                0.0
+            }
+        }
+    };
+    let result = rs_solve_fractional_ode(h0, alpha, t_horizon, n_steps, rust_rhs)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    if let Some(e) = cb_err.into_inner() {
+        return Err(e);
+    }
+    let dict = PyDict::new_bound(py);
+    dict.set_item("t_grid", result.t_grid)?;
+    dict.set_item("h", result.h)?;
+    Ok(dict.into())
+}
+
+/// Build a Markovian lift `K(t) ~= sum c_j exp(-gamma_j t)` on a
+/// geometric grid of rates with non-negative least-squares weights.
+///
+/// `kernel` is a Python callable `(t: float) -> float`.
+///
+/// Returns a dict `{"gammas": [...], "weights": [...]}`.
+#[pyfunction]
+#[pyo3(signature = (kernel, t_samples, n_factors, gamma_min, gamma_max, nnls_iter=5000))]
+fn geometric_grid_lift(
+    py: Python<'_>,
+    kernel: &Bound<'_, PyAny>,
+    t_samples: Vec<f64>,
+    n_factors: usize,
+    gamma_min: f64,
+    gamma_max: f64,
+    nnls_iter: usize,
+) -> PyResult<PyObject> {
+    let cb_err: RefCell<Option<PyErr>> = RefCell::new(None);
+    let rust_kernel = |t: f64| -> f64 {
+        if cb_err.borrow().is_some() {
+            return 0.0;
+        }
+        match kernel.call1((t,)).and_then(|v| v.extract::<f64>()) {
+            Ok(v) => v,
+            Err(e) => {
+                *cb_err.borrow_mut() = Some(e);
+                0.0
+            }
+        }
+    };
+    let lift = rs_geometric_grid_lift(
+        rust_kernel,
+        &t_samples,
+        n_factors,
+        gamma_min,
+        gamma_max,
+        nnls_iter,
+    )
+    .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    if let Some(e) = cb_err.into_inner() {
+        return Err(e);
+    }
+    let dict = PyDict::new_bound(py);
+    dict.set_item("gammas", lift.gammas)?;
+    dict.set_item("weights", lift.weights)?;
+    Ok(dict.into())
+}
+
+/// Solve the scalar second-kind Volterra equation
+/// `y(t) = g(t) + int_0^t K(t - s, y(s)) ds` by trapezoidal product
+/// integration on `n_steps + 1` equispaced nodes.
+///
+/// `g` is a Python callable `(t: float) -> float`.
+/// `kernel` is a Python callable `(dt: float, y: float) -> float`.
+///
+/// Returns a dict `{"t_grid": [...], "y": [...]}`.
+#[pyfunction]
+#[pyo3(signature = (g, kernel, t_horizon, n_steps, fixed_point_iter=50, fixed_point_tol=1e-12))]
+fn solve_volterra(
+    py: Python<'_>,
+    g: &Bound<'_, PyAny>,
+    kernel: &Bound<'_, PyAny>,
+    t_horizon: f64,
+    n_steps: usize,
+    fixed_point_iter: usize,
+    fixed_point_tol: f64,
+) -> PyResult<PyObject> {
+    let cb_err: RefCell<Option<PyErr>> = RefCell::new(None);
+    let rust_g = |t: f64| -> f64 {
+        if cb_err.borrow().is_some() {
+            return 0.0;
+        }
+        match g.call1((t,)).and_then(|v| v.extract::<f64>()) {
+            Ok(v) => v,
+            Err(e) => {
+                *cb_err.borrow_mut() = Some(e);
+                0.0
+            }
+        }
+    };
+    let rust_k = |dt: f64, y: f64| -> f64 {
+        if cb_err.borrow().is_some() {
+            return 0.0;
+        }
+        match kernel.call1((dt, y)).and_then(|v| v.extract::<f64>()) {
+            Ok(v) => v,
+            Err(e) => {
+                *cb_err.borrow_mut() = Some(e);
+                0.0
+            }
+        }
+    };
+    let result = rs_solve_volterra(
+        rust_g,
+        rust_k,
+        t_horizon,
+        n_steps,
+        fixed_point_iter,
+        fixed_point_tol,
+    )
+    .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    if let Some(e) = cb_err.into_inner() {
+        return Err(e);
+    }
+    let dict = PyDict::new_bound(py);
+    dict.set_item("t_grid", result.t_grid)?;
+    dict.set_item("y", result.y)?;
+    Ok(dict.into())
+}
+
+/// Recover a probability density on `x_grid` from a characteristic
+/// function `phi`.
+///
+/// `phi` is a Python callable `(u: float) -> (re: float, im: float)`.
+///
+/// Returns a dict `{"x_grid": [...], "density": [...]}`.
+#[pyfunction]
+#[pyo3(signature = (phi, x_grid, u_max, n_u))]
+fn fourier_invert(
+    py: Python<'_>,
+    phi: &Bound<'_, PyAny>,
+    x_grid: Vec<f64>,
+    u_max: f64,
+    n_u: usize,
+) -> PyResult<PyObject> {
+    let cb_err: RefCell<Option<PyErr>> = RefCell::new(None);
+    let rust_phi = |u: f64| -> (f64, f64) {
+        if cb_err.borrow().is_some() {
+            return (0.0, 0.0);
+        }
+        match phi.call1((u,)).and_then(|v| v.extract::<(f64, f64)>()) {
+            Ok(v) => v,
+            Err(e) => {
+                *cb_err.borrow_mut() = Some(e);
+                (0.0, 0.0)
+            }
+        }
+    };
+    let result = rs_fourier_invert(rust_phi, &x_grid, u_max, n_u)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    if let Some(e) = cb_err.into_inner() {
+        return Err(e);
+    }
+    let dict = PyDict::new_bound(py);
+    dict.set_item("x_grid", result.x_grid)?;
+    dict.set_item("density", result.density)?;
+    Ok(dict.into())
+}
+
+/// Register all Volterra-related functions with the Python module.
+pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(solve_fractional_ode, m)?)?;
+    m.add_function(wrap_pyfunction!(geometric_grid_lift, m)?)?;
+    m.add_function(wrap_pyfunction!(solve_volterra, m)?)?;
+    m.add_function(wrap_pyfunction!(fourier_invert, m)?)?;
+    Ok(())
+}