From d780ed81d7df8197bc853830cf60f24910e5bf90 Mon Sep 17 00:00:00 2001 From: ThotDjehuty Date: Tue, 12 May 2026 10:59:09 +0200 Subject: [PATCH] release(v1.1.0): additive CPU-only generic numerical primitives Adds 9 new top-level / sub-modules to the Rust API only (no Python bindings yet), with at least one analytic unit test per module. New Rust modules: - optimal_control::matrix_riccati (RK4 backward solver) - timeseries_utils::nonsync_covariance (Hayashi-Yoshida) - timeseries_utils::wavelet (Haar / Daubechies DWT and MODWT) - risk_measures (VaR, CVaR, projected sub-gradient CVaR minimisation) - graph::laplacian + graph::spectral_clustering (Jacobi + k-means++) - topology (Vietoris-Rips persistent homology, bottleneck distance) - volterra (Caputo Adams, Markovian lift, second-kind Volterra, Fourier inversion of characteristic functions) - signatures (truncated tensor signature, log-sig, random reservoir, Salvi-Cass-Lyons signature kernel, shuffle product) All previously stable APIs untouched; abi3-py38 ABI preserved. New module tests: 29/29 passing. Pre-existing 5 unrelated failures unchanged. --- CHANGELOG.md | 60 +++ Cargo.toml | 2 +- README.md | 17 +- docs/source/algorithms/graph_spectral.rst | 52 +++ docs/source/algorithms/matrix_riccati.rst | 60 +++ docs/source/algorithms/nonsync_covariance.rst | 48 +++ docs/source/algorithms/risk_measures.rst | 73 ++++ docs/source/algorithms/signatures.rst | 119 ++++++ docs/source/algorithms/topology.rst | 67 ++++ docs/source/algorithms/volterra.rst | 122 ++++++ docs/source/algorithms/wavelet.rst | 51 +++ docs/source/index.rst | 13 + pyproject.toml | 2 +- src/graph/laplacian.rs | 108 ++++++ src/graph/mod.rs | 14 + src/graph/spectral_clustering.rs | 286 ++++++++++++++ src/lib.rs | 7 + src/optimal_control/matrix_riccati.rs | 283 ++++++++++++++ src/optimal_control/mod.rs | 2 + src/risk_measures/cvar.rs | 245 ++++++++++++ src/risk_measures/mod.rs | 27 ++ src/risk_measures/var.rs | 120 ++++++ src/signatures/log_signature.rs | 105 ++++++ src/signatures/mod.rs | 16 + src/signatures/path_signature.rs | 191 ++++++++++ src/signatures/random_signature.rs | 152 ++++++++ src/signatures/signature_kernel.rs | 107 ++++++ src/signatures/utils.rs | 90 +++++ src/timeseries_utils.rs | 3 + src/timeseries_utils/nonsync_covariance.rs | 254 +++++++++++++ src/timeseries_utils/wavelet.rs | 216 +++++++++++ src/topology/bottleneck.rs | 226 +++++++++++ src/topology/mod.rs | 16 + src/topology/persistent_homology.rs | 353 ++++++++++++++++++ src/volterra/fourier_inversion.rs | 110 ++++++ src/volterra/fractional_riccati.rs | 156 ++++++++ src/volterra/markovian_lift.rs | 158 ++++++++ src/volterra/mod.rs | 16 + src/volterra/volterra_solver.rs | 100 +++++ 39 files changed, 4044 insertions(+), 3 deletions(-) create mode 100644 CHANGELOG.md create mode 100644 docs/source/algorithms/graph_spectral.rst create mode 100644 docs/source/algorithms/matrix_riccati.rst create mode 100644 docs/source/algorithms/nonsync_covariance.rst create mode 100644 docs/source/algorithms/risk_measures.rst create mode 100644 docs/source/algorithms/signatures.rst create mode 100644 docs/source/algorithms/topology.rst create mode 100644 docs/source/algorithms/volterra.rst create mode 100644 docs/source/algorithms/wavelet.rst create mode 100644 src/graph/laplacian.rs create mode 100644 src/graph/mod.rs create mode 100644 src/graph/spectral_clustering.rs create mode 100644 src/optimal_control/matrix_riccati.rs create mode 100644 src/risk_measures/cvar.rs create mode 100644 src/risk_measures/mod.rs create mode 100644 src/risk_measures/var.rs create mode 100644 src/signatures/log_signature.rs create mode 100644 src/signatures/mod.rs create mode 100644 src/signatures/path_signature.rs create mode 100644 src/signatures/random_signature.rs create mode 100644 src/signatures/signature_kernel.rs create mode 100644 src/signatures/utils.rs create mode 100644 src/timeseries_utils/nonsync_covariance.rs create mode 100644 src/timeseries_utils/wavelet.rs create mode 100644 src/topology/bottleneck.rs create mode 100644 src/topology/mod.rs create mode 100644 src/topology/persistent_homology.rs create mode 100644 src/volterra/fourier_inversion.rs create mode 100644 src/volterra/fractional_riccati.rs create mode 100644 src/volterra/markovian_lift.rs create mode 100644 src/volterra/mod.rs create mode 100644 src/volterra/volterra_solver.rs diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..b91e29d --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,60 @@ +# 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). + +## [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. diff --git a/Cargo.toml b/Cargo.toml index 70faac9..f2c362d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "optimiz-rs" -version = "1.0.1" +version = "1.1.0" edition = "2021" authors = ["HFThot Research Lab "] description = "High-performance optimization algorithms in Rust with Python bindings" diff --git a/README.md b/README.md index 1c43a92..526f552 100644 --- a/README.md +++ b/README.md @@ -6,13 +6,28 @@ **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-1.1.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/) 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 50-100× speedup over pure Python implementations. +## ✨ 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 diff --git a/docs/source/algorithms/graph_spectral.rst b/docs/source/algorithms/graph_spectral.rst new file mode 100644 index 0000000..f2cd166 --- /dev/null +++ b/docs/source/algorithms/graph_spectral.rst @@ -0,0 +1,52 @@ +Graph Laplacians and Spectral Clustering +======================================== + +The module :code:`graph` provides graph Laplacian operators and a +spectral clustering algorithm built on a Jacobi diagonaliser. + +Laplacians +---------- + +For a weighted undirected graph with adjacency matrix :math:`W \in \mathbb{R}^{n\times n}_{\ge 0}` +and degree matrix :math:`D = \mathrm{diag}(W \mathbf{1})`: + +- **Combinatorial**: :math:`L = D - W`. +- **Symmetric normalised**: :math:`L_{\mathrm{sym}} = I - D^{-1/2} W D^{-1/2}`. +- **Random-walk normalised**: :math:`L_{\mathrm{rw}} = I - D^{-1} W`. + +Each operator is positive semidefinite and the multiplicity of the +zero eigenvalue equals the number of connected components. + +Spectral Clustering +------------------- + +Given :math:`W` and a target number of clusters :math:`k`: + +1. Build :math:`L_{\mathrm{sym}}` (or another Laplacian). +2. Diagonalise via cyclic Jacobi rotations to obtain the eigenpairs + :math:`(\lambda_i, u_i)`. +3. Stack the :math:`k` eigenvectors associated with the smallest + eigenvalues as columns of :math:`U \in \mathbb{R}^{n \times k}`. +4. Normalise rows of :math:`U` and run Lloyd's algorithm with + k-means++ initialisation on the rows. + +The Fiedler eigenvalue :math:`\lambda_2` is reported separately as a +proxy for the spectral gap. + +API +--- + +.. code-block:: rust + + pub enum LaplacianKind { Combinatorial, SymmetricNormalised, RandomWalk } + pub fn combinatorial_laplacian(w: ArrayView2) -> Result>; + pub fn normalised_laplacian(w: ArrayView2) -> Result>; + pub fn random_walk_laplacian(w: ArrayView2) -> Result>; + + pub struct SpectralClusterResult { + pub labels: Vec, + pub eigenvalues: Vec, + pub fiedler_value: f64, + } + pub fn spectral_cluster(w: ArrayView2, k: usize, n_kmeans_iter: usize, seed: u64) + -> Result; diff --git a/docs/source/algorithms/matrix_riccati.rst b/docs/source/algorithms/matrix_riccati.rst new file mode 100644 index 0000000..440651b --- /dev/null +++ b/docs/source/algorithms/matrix_riccati.rst @@ -0,0 +1,60 @@ +Matrix Riccati Solver +===================== + +The module :code:`optimal_control::matrix_riccati` integrates backward in time +the matrix Riccati differential equation + +.. math:: + + \frac{dA(t)}{dt} = -2\,A(t)\,M\,A(t) + Q, + \qquad A(T) = A_T, + +together with the affine and constant components + +.. math:: + + \frac{dB(t)}{dt} = -2\,A(t)\,M\,B(t), + \qquad B(T) = B_T, + +.. math:: + + \frac{dC(t)}{dt} = -B(t)^\top\,M\,B(t), + \qquad C(T) = C_T. + +Discretisation +-------------- + +The grid :math:`\{t_n = T - n\,\Delta t\}_{n=0}^{N}` with +:math:`\Delta t = T / N` is traversed backward and a classical RK4 step is +applied to the joint vector field :math:`(A, B, C)`. Each macro step is +optionally subdivided into :math:`s` sub-steps for stability on stiff +problems. + +Validation +---------- + +In the scalar case :math:`A, M, Q \in \mathbb{R}` with :math:`A(T) = 0`, + +.. math:: + + A(t) \;=\; -\sqrt{\frac{Q}{2M}}\;\tanh\!\Big(\sqrt{2QM}\,(T - t)\Big), + +a closed form used by the unit test :code:`scalar_riccati_matches_analytic` +to certify :math:`L^\infty` convergence below :math:`10^{-5}` on +:math:`[0, T]`. + +API +--- + +.. code-block:: rust + + pub fn solve_matrix_riccati( + m_matrix: ArrayView2, + q: ArrayView2, + n: ArrayView2, + a_terminal: ArrayView2, + b_terminal: ArrayView1, + c_terminal: f64, + t_horizon: f64, + config: RiccatiConfig, + ) -> Result; diff --git a/docs/source/algorithms/nonsync_covariance.rst b/docs/source/algorithms/nonsync_covariance.rst new file mode 100644 index 0000000..5dde5f3 --- /dev/null +++ b/docs/source/algorithms/nonsync_covariance.rst @@ -0,0 +1,48 @@ +Asynchronous Covariance (Hayashi--Yoshida) +========================================== + +The module :code:`timeseries_utils::nonsync_covariance` implements the +Hayashi--Yoshida estimator of the integrated covariance between two +asynchronously sampled processes :math:`X` and :math:`Y` observed at +distinct, non-overlapping observation grids +:math:`\{t^X_i\}` and :math:`\{t^Y_j\}`. + +Estimator +--------- + +Let :math:`I_i = (t^X_{i-1}, t^X_i]` and :math:`J_j = (t^Y_{j-1}, t^Y_j]`. +The Hayashi--Yoshida estimator is + +.. math:: + + \widehat{\langle X, Y\rangle}_{[0,T]} + \;=\; + \sum_{i, j}\, + \big(X_{t^X_i} - X_{t^X_{i-1}}\big)\, + \big(Y_{t^Y_j} - Y_{t^Y_{j-1}}\big)\, + \mathbf{1}\!\big[I_i \cap J_j \neq \emptyset\big]. + +It is consistent under non-synchronicity and avoids the *Epps effect* +that plagues naive grid interpolation. + +Implementation +-------------- + +* A two-pointer scan in :math:`O(n_X + n_Y)` collects all overlapping + pairs. +* For matrices of size :math:`d \times d` with large per-asset sample + counts, off-diagonal entries are computed in parallel with Rayon. + +API +--- + +.. code-block:: rust + + pub fn hayashi_yoshida_covariance( + t1: &[f64], v1: &[f64], + t2: &[f64], v2: &[f64], + ) -> Result; + + pub fn hayashi_yoshida_matrix( + series: &[(Vec, Vec)], + ) -> Result>>; diff --git a/docs/source/algorithms/risk_measures.rst b/docs/source/algorithms/risk_measures.rst new file mode 100644 index 0000000..8b7bc5c --- /dev/null +++ b/docs/source/algorithms/risk_measures.rst @@ -0,0 +1,73 @@ +Risk Measures: VaR and CVaR +============================ + +The module :code:`risk_measures` provides Value-at-Risk and Conditional +Value-at-Risk estimators together with a convex CVaR minimisation +solver over the unit simplex. + +Definitions +----------- + +For a real random variable :math:`L` (a *loss*), the Value-at-Risk at +confidence level :math:`\alpha \in (0, 1)` is the lower :math:`\alpha`- +quantile + +.. math:: + + \mathrm{VaR}_\alpha(L) + \;=\; + \inf\!\big\{ \ell \in \mathbb{R} : \mathbb{P}(L \le \ell) \ge \alpha \big\}. + +The Conditional Value-at-Risk (also called Average Value-at-Risk) is + +.. math:: + + \mathrm{CVaR}_\alpha(L) + \;=\; + \frac{1}{1-\alpha}\, + \int_\alpha^1 \mathrm{VaR}_u(L)\,du. + +For a sample :math:`L_1, \dots, L_n` of i.i.d. losses sorted in increasing +order, the empirical CVaR at level :math:`\alpha` is + +.. math:: + + \widehat{\mathrm{CVaR}}_\alpha + \;=\; + \frac{1}{n - k}\, \sum_{i = k+1}^{n} L_{(i)}, + \qquad k = \lfloor \alpha\, n \rfloor. + +Convex minimisation +------------------- + +Rockafellar--Uryasev (2000) showed that + +.. math:: + + \mathrm{CVaR}_\alpha(L) + \;=\; + \min_{\zeta \in \mathbb{R}}\; + \zeta + \frac{1}{1 - \alpha}\,\mathbb{E}\!\big[(L - \zeta)_+\big]. + +Given samples of a vector :math:`r^{(s)} \in \mathbb{R}^d`, +:code:`minimize_cvar` solves + +.. math:: + + \min_{w \in \Delta_d,\;\zeta \in \mathbb{R}}\; + \zeta + \frac{1}{(1 - \alpha)\, S}\, \sum_{s=1}^S + \big(\zeta - \langle r^{(s)}, w\rangle\big)_+, + +over the unit simplex :math:`\Delta_d`, by a projected sub-gradient +method using the Held--Wolfe--Crowder simplex projection. + +API +--- + +.. code-block:: rust + + pub fn historical_var(losses: &[f64], alpha: f64) -> Result; + pub fn parametric_var(mu: f64, sigma: f64, alpha: f64) -> Result; + pub fn cvar_value(losses: &[f64], alpha: f64) -> Result; + pub fn minimize_cvar(returns: ArrayView2, cfg: &CVaRConfig) + -> Result; diff --git a/docs/source/algorithms/signatures.rst b/docs/source/algorithms/signatures.rst new file mode 100644 index 0000000..82f9dcd --- /dev/null +++ b/docs/source/algorithms/signatures.rst @@ -0,0 +1,119 @@ +Path Signatures +=============== + +The module :code:`signatures` provides truncated tensor signatures +(Lyons 1998), log-signatures, random reservoir projections, and the +Salvi--Cass--Lyons signature kernel. + +Truncated Signature +------------------- + +For a continuous path :math:`X : [0, T] \to \mathbb{R}^d` of bounded +variation, the *signature* is the formal series + +.. math:: + + S(X)_{0,T} + \;=\; + 1 + \sum_{k \ge 1} \sum_{i_1, \dots, i_k} + S^{i_1, \dots, i_k}_{0, T}\, + e_{i_1} \otimes \dots \otimes e_{i_k}, + +with iterated Stieltjes integrals + +.. math:: + + S^{i_1, \dots, i_k}_{0, T} + \;=\; + \int_{0 < u_1 < \dots < u_k < T} + dX^{i_1}_{u_1}\, \dots\, dX^{i_k}_{u_k}. + +For piecewise-linear input with increments :math:`\Delta_n`, the +truncated signature obeys the multiplicative recursion + +.. math:: + + S^{(M)}_{0, t_n} + \;=\; + S^{(M)}_{0, t_{n-1}}\,\otimes_M\,\exp_M(\Delta_n), + +where :math:`\exp_M(\Delta) = \sum_{k=0}^M \Delta^{\otimes k} / k!`. + +Log-Signature +------------- + +The truncated tensor logarithm + +.. math:: + + \log(S) + \;=\; + \sum_{n \ge 1} \frac{(-1)^{n+1}}{n}\,(S - 1)^{\otimes n} + +lives in the truncated free Lie algebra and provides a more +parsimonious representation. + +Random Signature +---------------- + +Following Cuchiero--Schmocker--Teichmann (2023), one drives a random +reservoir on :math:`\mathbb{R}^N`, + +.. math:: + + dZ_t = A_0 Z_t\, dt + \sum_{i=1}^d A_i Z_t\, dX^i_t, + +with random matrices :math:`A_i \in \mathbb{R}^{N \times N}` whose +entries are i.i.d. Gaussian with variance :math:`1/N`. The map +:math:`X \mapsto Z_T` is a finite-dimensional random projection of +:math:`S(X)`. + +Signature Kernel (Salvi--Cass--Lyons) +------------------------------------- + +The signature inner product + +.. math:: + + K(s, t) \;=\; \langle S(X)_{0, s},\; S(Y)_{0, t}\rangle + +solves the linear hyperbolic PDE + +.. math:: + + \frac{\partial^2 K}{\partial s\,\partial t} + \;=\; + \langle \dot X_s, \dot Y_t \rangle\, K(s, t), + \qquad + K(s, 0) = K(0, t) = 1. + +It is integrated on a uniform grid via the Goursat scheme + +.. math:: + + K_{i+1, j+1} + = K_{i+1, j} + K_{i, j+1} - K_{i, j} + + \langle \Delta x_i, \Delta y_j\rangle\, + \tfrac{1}{2}(K_{i+1, j} + K_{i, j+1}). + +API +--- + +.. code-block:: rust + + pub struct TruncatedSignature { + pub channels: usize, + pub level: usize, + pub tensors: Vec>, + } + pub fn path_signature(path: &[Vec], level: usize) -> Result; + pub fn log_signature(sig: &TruncatedSignature) -> Result; + + pub struct RandomSignatureConfig { + pub reservoir_dim: usize, pub seed: u64, pub variance: f64, + } + pub fn random_signature(path: &[Vec], cfg: &RandomSignatureConfig) + -> Result; + + pub fn signature_kernel(x: &[Vec], y: &[Vec]) + -> Result; diff --git a/docs/source/algorithms/topology.rst b/docs/source/algorithms/topology.rst new file mode 100644 index 0000000..5843909 --- /dev/null +++ b/docs/source/algorithms/topology.rst @@ -0,0 +1,67 @@ +Topological Data Analysis +========================= + +The module :code:`topology` implements Vietoris--Rips persistent +homology and the bottleneck distance between persistence diagrams. + +Vietoris--Rips Filtration +------------------------- + +For a finite point cloud :math:`X = \{x_1, \dots, x_n\} \subset \mathbb{R}^d` +and scale :math:`\varepsilon \ge 0`, the Vietoris--Rips complex is + +.. math:: + + \mathrm{VR}_\varepsilon(X) + \;=\; + \big\{ \sigma \subseteq X : \mathrm{diam}(\sigma) \le \varepsilon \big\}. + +Increasing :math:`\varepsilon` yields a filtration; the persistent +homology of this filtration produces, for each homological degree +:math:`k`, a multiset of birth/death pairs + +.. math:: + + D_k(X) = \big\{ (b_i, d_i) : 0 \le b_i < d_i \le \infty \big\}. + +Persistence Algorithm +--------------------- + +The boundary matrix :math:`\partial` is built over :math:`\mathbb{Z}/2` +and reduced left-to-right: for each column :math:`j` we cancel its +lowest entry by adding any earlier column with the same low. Pairs +:math:`(\mathrm{low}(j), j)` give birth/death pairs. + +Bottleneck Distance +------------------- + +For two diagrams :math:`D` and :math:`D'`, + +.. math:: + + d_B(D, D') + \;=\; + \inf_{\eta : D \to D'}\; + \sup_{x \in D}\, \|x - \eta(x)\|_\infty, + +where matchings may pair points with the diagonal +:math:`\Delta = \{(t, t) : t \ge 0\}` at cost :math:`(d - b)/2`. + +The implementation binary-searches the threshold :math:`\varepsilon` +and certifies a perfect matching by Hopcroft--Karp on the bipartite +graph of admissible edges. + +API +--- + +.. code-block:: rust + + pub struct PersistencePair { pub dim: usize, pub birth: f64, pub death: f64 } + pub struct PersistenceDiagram { pub pairs: Vec } + + pub fn vietoris_rips_filtration(points: &[Vec], max_dim: usize, max_eps: f64) + -> Result>; + pub fn persistent_homology(points: &[Vec], max_dim: usize, max_eps: f64) + -> Result; + pub fn bottleneck_distance(d1: &[PersistencePair], d2: &[PersistencePair]) + -> Result; diff --git a/docs/source/algorithms/volterra.rst b/docs/source/algorithms/volterra.rst new file mode 100644 index 0000000..89b7fe4 --- /dev/null +++ b/docs/source/algorithms/volterra.rst @@ -0,0 +1,122 @@ +Volterra and Fractional Solvers +================================ + +The module :code:`volterra` collects four CPU-only generic numerical +primitives for Volterra integral equations and related transforms. + +Fractional Caputo Adams Solver +------------------------------ + +For :math:`\alpha \in (0, 1)`, solve + +.. math:: + + D^\alpha h(t) = F(t, h(t)), + \qquad h(0) = h_0, + +with the Diethelm--Ford--Freed (2002) fractional Adams predictor-- +corrector. Predictor: + +.. math:: + + h^P_{n+1} + = h_0 + + \frac{\Delta t^\alpha}{\alpha\,\Gamma(\alpha)} + \sum_{k=0}^{n} + \big[(n+1-k)^\alpha - (n-k)^\alpha\big]\, F(t_k, h_k). + +Corrector: + +.. math:: + + h_{n+1} + = h_0 + + \frac{\Delta t^\alpha}{\Gamma(\alpha + 2)} + \Big[ F(t_{n+1}, h^P_{n+1}) + \sum_{k=0}^{n} a_{n+1, k}\, F(t_k, h_k) \Big], + +with + +.. math:: + + a_{n+1, 0} = n^{\alpha + 1} - (n - \alpha)\,(n+1)^\alpha, + + a_{n+1, k} = (n - k + 2)^{\alpha + 1} + (n - k)^{\alpha + 1} + - 2\,(n - k + 1)^{\alpha + 1}, + \qquad 1 \le k \le n. + +Markovian Lift +-------------- + +A convolution kernel :math:`K(t)` admitting + +.. math:: + + K(t) = \int_0^\infty e^{-\gamma t}\, \nu(d\gamma) + +is approximated by + +.. math:: + + K(t) \;\approx\; \sum_{j=1}^N c_j\, e^{-\gamma_j t}, + \qquad c_j \ge 0, + +with rates :math:`\gamma_j` on a geometric grid and weights fitted by +non-negative least squares. + +Generic Volterra Equation +------------------------- + +For + +.. math:: + + y(t) = g(t) + \int_0^t K(t - s, y(s))\, ds, + +the trapezoidal product-integration scheme reads + +.. math:: + + y_n = g_n + \Delta t\,\Big[ + \tfrac{1}{2} K(t_n, y_0) + + \sum_{k=1}^{n-1} K(t_n - t_k, y_k) + + \tfrac{1}{2} K(0, y_n) \Big], + +solved implicitly by fixed-point iteration on :math:`y_n`. + +Fourier Inversion +----------------- + +Recover a density from a characteristic function :math:`\varphi(u)` via + +.. math:: + + f(x) \;\approx\; + \frac{\Delta u}{\pi}\, + \sum_{k=0}^{N_u - 1} + w_k \big[\,\Re \varphi(u_k)\,\cos(u_k x) + + \Im \varphi(u_k)\,\sin(u_k x)\,\big], + +with trapezoidal weights :math:`w_k`. + +API +--- + +.. code-block:: rust + + pub fn solve_fractional_ode f64>( + h0: f64, alpha: f64, t_horizon: f64, n_steps: usize, rhs: F, + ) -> Result; + + pub fn geometric_grid_lift f64>( + kernel: K, t_samples: &[f64], + n_factors: usize, gamma_min: f64, gamma_max: f64, nnls_iter: usize, + ) -> Result; + + pub fn solve_volterra( + g: G, kernel: K, t_horizon: f64, n_steps: usize, + fixed_point_iter: usize, fixed_point_tol: f64, + ) -> Result; + + pub fn fourier_invert (f64, f64)>( + phi: P, x_grid: &[f64], u_max: f64, n_u: usize, + ) -> Result; diff --git a/docs/source/algorithms/wavelet.rst b/docs/source/algorithms/wavelet.rst new file mode 100644 index 0000000..a88ceb8 --- /dev/null +++ b/docs/source/algorithms/wavelet.rst @@ -0,0 +1,51 @@ +Discrete and Maximum-Overlap Wavelet Transforms +================================================ + +The module :code:`timeseries_utils::wavelet` provides Haar and Daubechies +wavelet transforms with periodic boundary handling. + +Filter banks +------------ + +For an orthogonal scaling filter :math:`\{h_k\}_{k=0}^{L-1}` the quadrature +mirror filter (QMF) is + +.. math:: + + g_k = (-1)^k\, h_{L - 1 - k}, + +so that :math:`\sum_k h_k = \sqrt{2}` and :math:`\sum_k g_k = 0`. + +DWT (one level, periodic) +------------------------- + +For an input vector :math:`x \in \mathbb{R}^N` with :math:`N` even, + +.. math:: + + a_n = \sum_{k=0}^{L-1} h_k\, x_{(2n + k)\bmod N}, + \qquad + d_n = \sum_{k=0}^{L-1} g_k\, x_{(2n + k)\bmod N}, + \qquad n = 0, \dots, N/2 - 1. + +Successive levels apply the same filter to the previous approximation +:math:`a^{(j)}`. + +MODWT (Maximum Overlap) +----------------------- + +The MODWT does not downsample: at level :math:`j`, the filter is dilated +by inserting :math:`2^{j-1} - 1` zeros between successive taps and applied +in a periodic convolution. The result is shift-invariant. + +API +--- + +.. code-block:: rust + + pub enum WaveletFamily { Haar, Daubechies(u8) } + pub fn scaling_filter(family: WaveletFamily) -> Result>; + pub fn qmf(h: &[f64]) -> Vec; + pub fn dwt_step(x: &[f64], h: &[f64], g: &[f64]) -> (Vec, Vec); + pub fn dwt(x: &[f64], family: WaveletFamily, levels: usize) -> Result>>; + pub fn modwt_step(x: &[f64], h: &[f64], g: &[f64], level: usize) -> (Vec, Vec); diff --git a/docs/source/index.rst b/docs/source/index.rst index 9a47739..4ff4349 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -38,6 +38,19 @@ Optimiz-rs provides blazingly fast, production-ready implementations of advanced algorithms/grid_search algorithms/point_processes +.. toctree:: + :maxdepth: 2 + :caption: v1.1 Numerical Primitives + + algorithms/matrix_riccati + algorithms/nonsync_covariance + algorithms/wavelet + algorithms/risk_measures + algorithms/graph_spectral + algorithms/topology + algorithms/volterra + algorithms/signatures + .. toctree:: :maxdepth: 2 :caption: API Reference diff --git a/pyproject.toml b/pyproject.toml index 9aabce0..da61aa8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "optimiz-rs" -version = "1.0.1" +version = "1.1.0" description = "High-performance optimization algorithms in Rust with Python bindings" authors = [ {name = "HFThot Research Lab", email = "contact@hfthot-lab.eu"} diff --git a/src/graph/laplacian.rs b/src/graph/laplacian.rs new file mode 100644 index 0000000..29a173d --- /dev/null +++ b/src/graph/laplacian.rs @@ -0,0 +1,108 @@ +//! Graph Laplacians for dense weighted similarity matrices. +//! +//! For a non-negative symmetric weight matrix `W in R^{n x n}` and degree +//! diagonal `D = diag(W * 1)`, the three standard Laplacians are +//! +//! ```text +//! L = D - W (combinatorial) +//! L_sym = I - D^{-1/2} W D^{-1/2} (symmetric normalised) +//! L_rw = I - D^{-1} W (random-walk) +//! ``` +//! +//! All operate on dense `ndarray::Array2`. + +use ndarray::{Array2, ArrayView2}; + +use crate::core::{OptimizrError, Result}; + +/// Laplacian variant. +#[derive(Debug, Clone, Copy)] +pub enum LaplacianKind { + Combinatorial, + SymmetricNormalised, + RandomWalk, +} + +fn check_weight(w: ArrayView2) -> Result<()> { + let (n, m) = (w.nrows(), w.ncols()); + if n == 0 || n != m { + return Err(OptimizrError::InvalidInput( + "weight matrix must be square and non-empty".into(), + )); + } + for &x in w.iter() { + if !x.is_finite() || x < 0.0 { + return Err(OptimizrError::InvalidInput( + "weight matrix must be non-negative and finite".into(), + )); + } + } + Ok(()) +} + +/// Combinatorial Laplacian `L = D - W`. +pub fn combinatorial_laplacian(w: ArrayView2) -> Result> { + check_weight(w)?; + let n = w.nrows(); + let mut l = -w.to_owned(); + for i in 0..n { + let d_i: f64 = w.row(i).sum(); + l[[i, i]] += d_i; + } + Ok(l) +} + +/// Symmetric normalised Laplacian `L_sym = I - D^{-1/2} W D^{-1/2}`. +pub fn normalised_laplacian(w: ArrayView2) -> Result> { + check_weight(w)?; + let n = w.nrows(); + let mut d_inv_sqrt = vec![0.0; n]; + for i in 0..n { + let d_i: f64 = w.row(i).sum(); + d_inv_sqrt[i] = if d_i > 0.0 { 1.0 / d_i.sqrt() } else { 0.0 }; + } + let mut l = Array2::::zeros((n, n)); + for i in 0..n { + for j in 0..n { + let off = -d_inv_sqrt[i] * w[[i, j]] * d_inv_sqrt[j]; + l[[i, j]] = if i == j { 1.0 + off } else { off }; + } + } + Ok(l) +} + +/// Random-walk Laplacian `L_rw = I - D^{-1} W`. +pub fn random_walk_laplacian(w: ArrayView2) -> Result> { + check_weight(w)?; + let n = w.nrows(); + let mut l = Array2::::zeros((n, n)); + for i in 0..n { + let d_i: f64 = w.row(i).sum(); + let inv = if d_i > 0.0 { 1.0 / d_i } else { 0.0 }; + for j in 0..n { + l[[i, j]] = if i == j { 1.0 - inv * w[[i, j]] } else { -inv * w[[i, j]] }; + } + } + Ok(l) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn laplacian_kernel_includes_constant() { + let w = array![ + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + ]; + let l = combinatorial_laplacian(w.view()).unwrap(); + let one = ndarray::Array1::::ones(3); + let lone = l.dot(&one); + for v in lone.iter() { + assert!(v.abs() < 1e-12); + } + } +} diff --git a/src/graph/mod.rs b/src/graph/mod.rs new file mode 100644 index 0000000..c44cd05 --- /dev/null +++ b/src/graph/mod.rs @@ -0,0 +1,14 @@ +//! Generic graph algorithms. +//! +//! Currently exposes spectral utilities on dense weighted graphs: +//! +//! * `laplacian` -- combinatorial / normalised / random-walk Laplacians. +//! * `spectral_clustering` -- Ng--Jordan--Weiss algorithm (k-means on +//! normalised eigenvectors of L_sym). + +pub mod laplacian; +pub mod spectral_clustering; + + +pub use laplacian::{combinatorial_laplacian, normalised_laplacian, random_walk_laplacian, LaplacianKind}; +pub use spectral_clustering::{spectral_cluster, SpectralClusterResult}; diff --git a/src/graph/spectral_clustering.rs b/src/graph/spectral_clustering.rs new file mode 100644 index 0000000..805d3ad --- /dev/null +++ b/src/graph/spectral_clustering.rs @@ -0,0 +1,286 @@ +//! Spectral clustering on dense weighted graphs (Ng--Jordan--Weiss). +//! +//! Algorithm: +//! +//! 1. Build `L_sym = I - D^{-1/2} W D^{-1/2}`. +//! 2. Compute the `k` eigenvectors of `L_sym` associated with the `k` +//! smallest eigenvalues, stack them into `U in R^{n x k}`. +//! 3. Normalise each row of `U` to unit `l_2` norm. +//! 4. Apply Lloyd's `k`-means to the rows of `U`. +//! +//! Eigen-decomposition uses a symmetric Jacobi rotation (no external +//! linear-algebra dependency) which is appropriate for the small +//! to medium dense matrices that spectral clustering targets. + +use ndarray::{Array1, Array2, ArrayView2}; + +use crate::core::{OptimizrError, Result}; +use super::laplacian::normalised_laplacian; + +/// Result of spectral clustering. +#[derive(Debug, Clone)] +pub struct SpectralClusterResult { + pub labels: Vec, + pub eigenvalues: Vec, + pub fiedler_value: f64, +} + +/// Spectral clustering of `n` items into `k` groups using a non-negative +/// symmetric similarity matrix `w`. +pub fn spectral_cluster( + w: ArrayView2, + k: usize, + n_kmeans_iter: usize, + seed: u64, +) -> Result { + if k < 2 { + return Err(OptimizrError::InvalidParameter("k must be >= 2".into())); + } + let n = w.nrows(); + if n < k { + return Err(OptimizrError::InvalidParameter( + "k must not exceed n".into(), + )); + } + let l_sym = normalised_laplacian(w)?; + let (eigvals, eigvecs) = jacobi_symmetric_eig(&l_sym, 200, 1e-10)?; + // Sort ascending + let mut order: Vec = (0..n).collect(); + order.sort_by(|&a, &b| eigvals[a].partial_cmp(&eigvals[b]).unwrap_or(std::cmp::Ordering::Equal)); + let smallest_k: Vec = order.iter().take(k).copied().collect(); + + // U[i, j] = eigvecs[i, smallest_k[j]] + let mut u = Array2::::zeros((n, k)); + for j in 0..k { + let col = smallest_k[j]; + for i in 0..n { + u[[i, j]] = eigvecs[[i, col]]; + } + } + // Row-normalise + for i in 0..n { + let norm: f64 = (0..k).map(|j| u[[i, j]] * u[[i, j]]).sum::().sqrt(); + if norm > 0.0 { + for j in 0..k { + u[[i, j]] /= norm; + } + } + } + + let labels = kmeans_lloyd(&u, k, n_kmeans_iter, seed); + let sorted_eigvals: Vec = order.iter().map(|&i| eigvals[i]).collect(); + let fiedler = sorted_eigvals.get(1).copied().unwrap_or(0.0); + Ok(SpectralClusterResult { + labels, + eigenvalues: sorted_eigvals, + fiedler_value: fiedler, + }) +} + +/// Lloyd's k-means with k-means++ seeding. +fn kmeans_lloyd(x: &Array2, k: usize, n_iter: usize, seed: u64) -> Vec { + let n = x.nrows(); + let d = x.ncols(); + let mut state = if seed == 0 { 0xDEAD_BEEFu64 } else { seed }; + let mut next = || { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + ((state >> 11) as f64) / ((1u64 << 53) as f64) + }; + let mut centers = Array2::::zeros((k, d)); + let first = (next() * n as f64) as usize % n; + for j in 0..d { + centers[[0, j]] = x[[first, j]]; + } + let mut min_d2 = vec![f64::INFINITY; n]; + for c in 1..k { + for i in 0..n { + let mut s = 0.0; + for j in 0..d { + let diff = x[[i, j]] - centers[[c - 1, j]]; + s += diff * diff; + } + if s < min_d2[i] { + min_d2[i] = s; + } + } + let total: f64 = min_d2.iter().sum(); + let r = next() * total; + let mut acc = 0.0; + let mut chosen = 0; + for i in 0..n { + acc += min_d2[i]; + if acc >= r { + chosen = i; + break; + } + } + for j in 0..d { + centers[[c, j]] = x[[chosen, j]]; + } + } + + let mut labels = vec![0usize; n]; + for _ in 0..n_iter { + let mut changed = false; + for i in 0..n { + let mut best = 0; + let mut best_d2 = f64::INFINITY; + for c in 0..k { + let mut s = 0.0; + for j in 0..d { + let diff = x[[i, j]] - centers[[c, j]]; + s += diff * diff; + } + if s < best_d2 { + best_d2 = s; + best = c; + } + } + if labels[i] != best { + labels[i] = best; + changed = true; + } + } + if !changed { + break; + } + let mut counts = vec![0usize; k]; + let mut new_centers = Array2::::zeros((k, d)); + for i in 0..n { + counts[labels[i]] += 1; + for j in 0..d { + new_centers[[labels[i], j]] += x[[i, j]]; + } + } + for c in 0..k { + if counts[c] > 0 { + for j in 0..d { + new_centers[[c, j]] /= counts[c] as f64; + } + } else { + for j in 0..d { + new_centers[[c, j]] = centers[[c, j]]; + } + } + } + centers = new_centers; + } + labels +} + +/// Symmetric Jacobi eigen-decomposition. Returns `(eigenvalues, V)` +/// such that `A V = V diag(eigenvalues)` and `V` is orthogonal. +pub fn jacobi_symmetric_eig( + a: &Array2, + max_sweeps: usize, + tol: f64, +) -> Result<(Array1, Array2)> { + let n = a.nrows(); + if n != a.ncols() { + return Err(OptimizrError::InvalidInput("matrix must be square".into())); + } + let mut m = a.clone(); + let mut v = Array2::::eye(n); + + for _ in 0..max_sweeps { + let mut off = 0.0; + for p in 0..n { + for q in (p + 1)..n { + off += m[[p, q]] * m[[p, q]]; + } + } + if off.sqrt() < tol { + break; + } + for p in 0..n { + for q in (p + 1)..n { + let apq = m[[p, q]]; + if apq.abs() < 1e-16 { + continue; + } + let app = m[[p, p]]; + let aqq = m[[q, q]]; + let theta = (aqq - app) / (2.0 * apq); + let t = if theta >= 0.0 { + 1.0 / (theta + (1.0 + theta * theta).sqrt()) + } else { + 1.0 / (theta - (1.0 + theta * theta).sqrt()) + }; + let c = 1.0 / (1.0 + t * t).sqrt(); + let s = t * c; + m[[p, p]] = app - t * apq; + m[[q, q]] = aqq + t * apq; + m[[p, q]] = 0.0; + m[[q, p]] = 0.0; + for i in 0..n { + if i != p && i != q { + let aip = m[[i, p]]; + let aiq = m[[i, q]]; + m[[i, p]] = c * aip - s * aiq; + m[[p, i]] = m[[i, p]]; + m[[i, q]] = s * aip + c * aiq; + m[[q, i]] = m[[i, q]]; + } + } + for i in 0..n { + let vip = v[[i, p]]; + let viq = v[[i, q]]; + v[[i, p]] = c * vip - s * viq; + v[[i, q]] = s * vip + c * viq; + } + } + } + } + + let eigvals = Array1::from_iter((0..n).map(|i| m[[i, i]])); + Ok((eigvals, v)) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn two_blocks_recovered() { + // 6 nodes, two cliques {0,1,2} and {3,4,5}, weak link 2-3. + let w = array![ + [0.0, 1.0, 1.0, 0.05, 0.0, 0.0], + [1.0, 0.0, 1.0, 0.0, 0.0, 0.0], + [1.0, 1.0, 0.0, 0.05, 0.0, 0.0], + [0.05, 0.0, 0.05, 0.0, 1.0, 1.0], + [0.0, 0.0, 0.0, 1.0, 0.0, 1.0], + [0.0, 0.0, 0.0, 1.0, 1.0, 0.0], + ]; + let res = spectral_cluster(w.view(), 2, 100, 7).unwrap(); + // Members of each block should share their label. + assert_eq!(res.labels[0], res.labels[1]); + assert_eq!(res.labels[1], res.labels[2]); + assert_eq!(res.labels[3], res.labels[4]); + assert_eq!(res.labels[4], res.labels[5]); + assert_ne!(res.labels[0], res.labels[3]); + } + + #[test] + fn jacobi_diagonalises_symmetric() { + let a = array![[2.0_f64, 1.0], [1.0, 3.0]]; + let (eig, v) = jacobi_symmetric_eig(&a, 200, 1e-12).unwrap(); + let mut sorted: Vec = eig.to_vec(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + // Eigenvalues of [[2,1],[1,3]]: (5 +- sqrt(5)) / 2 + let lo = (5.0 - 5.0_f64.sqrt()) / 2.0; + let hi = (5.0 + 5.0_f64.sqrt()) / 2.0; + assert!((sorted[0] - lo).abs() < 1e-9); + assert!((sorted[1] - hi).abs() < 1e-9); + // Orthogonality + let vt_v = v.t().dot(&v); + for i in 0..2 { + for j in 0..2 { + let target = if i == j { 1.0 } else { 0.0 }; + assert!((vt_v[[i, j]] - target).abs() < 1e-9); + } + } + } +} diff --git a/src/lib.rs b/src/lib.rs index 6d13ac6..450031c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -48,6 +48,13 @@ pub mod mean_field; // Mean Field Games and Mean Field Type Control pub mod point_processes; // Point processes for order flow modeling (Hawkes, fBM) pub mod portfolio_optimization; // CARA, convex duality, mean-variance, ERC +// ===== v1.1.0 additive modules (CPU-only generic numerical primitives) ===== +pub mod graph; // Graph Laplacians and spectral clustering +pub mod risk_measures; // Generic VaR / CVaR estimators and convex CVaR minimisation +pub mod signatures; // Path signatures, log-signatures, signature kernels +pub mod topology; // Vietoris--Rips persistent homology and bottleneck distance +pub mod volterra; // Fractional / Volterra integral equation solvers + // Python bindings for legacy compatibility #[cfg(feature = "python-bindings")] mod differential_evolution; diff --git a/src/optimal_control/matrix_riccati.rs b/src/optimal_control/matrix_riccati.rs new file mode 100644 index 0000000..d0e07df --- /dev/null +++ b/src/optimal_control/matrix_riccati.rs @@ -0,0 +1,283 @@ +//! Matrix Riccati Backward ODE Solver +//! ================================== +//! +//! Solves the coupled backward matrix Riccati system +//! +//! ```text +//! dA/dt = -2 A M A + Q, A(T) = A_T (symmetric) +//! dB/dt = -2 A M B + N^T B, B(T) = B_T +//! dC/dt = -2 A M C, C(T) = C_T +//! ``` +//! +//! using a classical fourth-order Runge--Kutta scheme integrated backward +//! in time with optional sub-stepping for stiff problems. +//! +//! All matrices live in `R^{d x d}`; `Q` is symmetric positive +//! semi-definite, `M` symmetric positive definite. The solver performs no +//! linear-algebra inversion: callers that have a pre-computed `M^{-1}` can +//! pass it directly through `m_matrix` since only the products `A M A` +//! appear in the right-hand side. +//! +//! # Reference +//! +//! Bismut (1976), *Linear-Quadratic Optimal Stochastic Control with Random +//! Coefficients*, SIAM Journal on Control and Optimization 14(3). + +use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis}; + +use crate::core::{OptimizrError, Result}; + +/// Configuration for the backward Riccati integrator. +#[derive(Debug, Clone, Copy)] +pub struct RiccatiConfig { + /// Number of stored time points (including both endpoints). + pub n_steps: usize, + /// Number of internal RK4 sub-steps between two stored points. + pub n_substeps: usize, +} + +impl Default for RiccatiConfig { + fn default() -> Self { + Self { + n_steps: 200, + n_substeps: 1, + } + } +} + +/// Result of a backward Riccati integration. +/// +/// All vectors are ordered forward in time, i.e. `t_grid[0] = 0` and +/// `t_grid[n_steps - 1] = T`. +#[derive(Debug, Clone)] +pub struct RiccatiResult { + pub t_grid: Vec, + pub a: Vec>, + pub b: Vec>, + pub c: Vec>, +} + +#[inline] +fn rhs( + a: &Array2, + b: &Array2, + c: &Array1, + m: &Array2, + q: &Array2, + n_t: &Array2, // N^T precomputed +) -> (Array2, Array2, Array1) { + let am = a.dot(m); + let ama = am.dot(a); + let amb = am.dot(b); + let amc = am.dot(c); + + let da = -2.0 * &ama + q; + let db = -2.0 * &amb + n_t.dot(b); + let dc = -2.0 * amc; + (da, db, dc) +} + +/// Backward integration of the matrix Riccati system. +/// +/// # Arguments +/// +/// * `m_matrix` -- `M` (symmetric positive definite) of shape `(d, d)`. +/// * `q` -- `Q` (symmetric positive semi-definite) of shape `(d, d)`. +/// * `n` -- `N` of shape `(d, d)`. +/// * `a_terminal` -- terminal `A(T)` of shape `(d, d)`. +/// * `b_terminal` -- terminal `B(T)` of shape `(d, d)`. +/// * `c_terminal` -- terminal `C(T)` of shape `(d,)`. +/// * `t_horizon` -- final time `T > 0`. +/// * `config` -- discretisation parameters. +pub fn solve_matrix_riccati( + m_matrix: ArrayView2, + q: ArrayView2, + n: ArrayView2, + a_terminal: ArrayView2, + b_terminal: ArrayView2, + c_terminal: ArrayView1, + t_horizon: f64, + config: RiccatiConfig, +) -> Result { + if t_horizon <= 0.0 { + return Err(OptimizrError::InvalidParameter( + "t_horizon must be strictly positive".into(), + )); + } + if config.n_steps < 2 { + return Err(OptimizrError::InvalidParameter( + "n_steps must be at least 2".into(), + )); + } + if config.n_substeps == 0 { + return Err(OptimizrError::InvalidParameter( + "n_substeps must be at least 1".into(), + )); + } + + let d = a_terminal.nrows(); + let check_square = |x: ArrayView2, name: &str| -> Result<()> { + if x.nrows() != d || x.ncols() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: x.nrows(), + }); + } + let _ = name; + Ok(()) + }; + check_square(m_matrix, "m_matrix")?; + check_square(q, "q")?; + check_square(n, "n")?; + check_square(b_terminal, "b_terminal")?; + if c_terminal.len() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: c_terminal.len(), + }); + } + + let m = m_matrix.to_owned(); + let q_o = q.to_owned(); + let n_t = n.t().to_owned(); + + let n_steps = config.n_steps; + let n_sub = config.n_substeps; + let dt = t_horizon / ((n_steps - 1) as f64); + let h = dt / (n_sub as f64); + + let mut t_grid = Vec::with_capacity(n_steps); + for k in 0..n_steps { + t_grid.push((k as f64) * dt); + } + + let mut a_vec: Vec> = vec![Array2::::zeros((d, d)); n_steps]; + let mut b_vec: Vec> = vec![Array2::::zeros((d, d)); n_steps]; + let mut c_vec: Vec> = vec![Array1::::zeros(d); n_steps]; + + a_vec[n_steps - 1] = a_terminal.to_owned(); + b_vec[n_steps - 1] = b_terminal.to_owned(); + c_vec[n_steps - 1] = c_terminal.to_owned(); + + // Backward integration: from index k to k-1 (time decreases). + // The continuous system is integrated with step -h (RK4). + for k in (1..n_steps).rev() { + let mut a = a_vec[k].clone(); + let mut b = b_vec[k].clone(); + let mut c = c_vec[k].clone(); + + for _ in 0..n_sub { + let (k1a, k1b, k1c) = rhs(&a, &b, &c, &m, &q_o, &n_t); + let a2 = &a - 0.5 * h * &k1a; + let b2 = &b - 0.5 * h * &k1b; + let c2 = &c - 0.5 * h * &k1c; + let (k2a, k2b, k2c) = rhs(&a2, &b2, &c2, &m, &q_o, &n_t); + let a3 = &a - 0.5 * h * &k2a; + let b3 = &b - 0.5 * h * &k2b; + let c3 = &c - 0.5 * h * &k2c; + let (k3a, k3b, k3c) = rhs(&a3, &b3, &c3, &m, &q_o, &n_t); + let a4 = &a - h * &k3a; + let b4 = &b - h * &k3b; + let c4 = &c - h * &k3c; + let (k4a, k4b, k4c) = rhs(&a4, &b4, &c4, &m, &q_o, &n_t); + + a = &a - (h / 6.0) * (&k1a + 2.0 * &k2a + 2.0 * &k3a + &k4a); + b = &b - (h / 6.0) * (&k1b + 2.0 * &k2b + 2.0 * &k3b + &k4b); + c = &c - (h / 6.0) * (&k1c + 2.0 * &k2c + 2.0 * &k3c + &k4c); + } + + // Symmetrise A to preserve numerical symmetry over long horizons. + let a_sym = 0.5 * (&a + &a.t()); + a_vec[k - 1] = a_sym; + b_vec[k - 1] = b; + c_vec[k - 1] = c; + } + + // Final integrity check: no NaN/Inf in stored arrays. + for arr in a_vec.iter() { + if arr.iter().any(|x| !x.is_finite()) { + return Err(OptimizrError::NumericalError( + "non-finite value encountered in A trajectory".into(), + )); + } + } + let _ = Axis(0); + Ok(RiccatiResult { + t_grid, + a: a_vec, + b: b_vec, + c: c_vec, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + /// Scalar (d=1) reference solution for A(T) = 0. + /// The ODE dA/dt = -2 M A^2 + Q with terminal A(T)=0 admits + /// A(t) = -sqrt(Q / (2 M)) * tanh(sqrt(2 Q M) * (T - t)). + #[test] + fn scalar_riccati_matches_analytic() { + let m = array![[1.5_f64]]; + let q = array![[0.8_f64]]; + let n = array![[0.0_f64]]; + let a_t = array![[0.0_f64]]; + let b_t = array![[0.0_f64]]; + let c_t = array![0.0_f64]; + + let t_horizon = 3.0; + let cfg = RiccatiConfig { + n_steps: 2001, + n_substeps: 4, + }; + let res = solve_matrix_riccati( + m.view(), + q.view(), + n.view(), + a_t.view(), + b_t.view(), + c_t.view(), + t_horizon, + cfg, + ) + .expect("solver"); + + let coeff = -(q[[0, 0]] / (2.0 * m[[0, 0]])).sqrt(); + let alpha = (2.0 * q[[0, 0]] * m[[0, 0]]).sqrt(); + + let mut max_err = 0.0_f64; + for (i, &t) in res.t_grid.iter().enumerate() { + let analytic = coeff * (alpha * (t_horizon - t)).tanh(); + let num = res.a[i][[0, 0]]; + max_err = max_err.max((analytic - num).abs()); + } + assert!( + max_err < 1e-5, + "Riccati scalar L_inf error too large: {}", + max_err + ); + } + + #[test] + fn rejects_invalid_inputs() { + let m = array![[1.0_f64]]; + let q = array![[1.0_f64]]; + let n = array![[0.0_f64]]; + let a_t = array![[0.0_f64]]; + let b_t = array![[0.0_f64]]; + let c_t = array![0.0_f64]; + assert!(solve_matrix_riccati( + m.view(), + q.view(), + n.view(), + a_t.view(), + b_t.view(), + c_t.view(), + -1.0, + RiccatiConfig::default(), + ) + .is_err()); + } +} diff --git a/src/optimal_control/mod.rs b/src/optimal_control/mod.rs index a14d5b0..7d2d082 100644 --- a/src/optimal_control/mod.rs +++ b/src/optimal_control/mod.rs @@ -38,6 +38,7 @@ pub mod jump_diffusion; pub mod kalman_filter; #[cfg(feature = "python-bindings")] pub mod kalman_py_bindings; +pub mod matrix_riccati; pub mod mrsjd; pub mod ou_estimator; #[cfg(feature = "python-bindings")] @@ -46,6 +47,7 @@ pub mod regime_switching; pub mod viscosity; pub use hjb_solver::{HJBConfig, HJBResult, HJBSolver}; +pub use matrix_riccati::{solve_matrix_riccati, RiccatiConfig, RiccatiResult}; pub use jump_diffusion::{ JumpDiffusionConfig, JumpDiffusionResult, JumpDiffusionSolver, JumpDistribution, }; diff --git a/src/risk_measures/cvar.rs b/src/risk_measures/cvar.rs new file mode 100644 index 0000000..57b9af6 --- /dev/null +++ b/src/risk_measures/cvar.rs @@ -0,0 +1,245 @@ +//! Conditional Value-at-Risk (CVaR / Expected Shortfall). +//! +//! Rockafellar--Uryasev (2000) variational representation: +//! +//! ```text +//! CVaR_alpha(L) = inf_{zeta in R} { zeta + 1/(1 - alpha) * E[(L - zeta)_+] }. +//! ``` +//! +//! The optimum is attained at `zeta* = VaR_alpha(L)`. For an empirical +//! sample `L^{(s)}, s = 1..S`, the closed-form estimator equals +//! +//! ```text +//! CVaR_alpha = mean( L^{(s)} : L^{(s)} >= VaR_alpha ). +//! ``` +//! +//! For decision-variable optimisation, we minimise +//! +//! ```text +//! min_{w in C, zeta, u} zeta + 1/((1 - alpha) * S) * sum_s u_s +//! s.t. u_s >= - - zeta, u_s >= 0 +//! ``` +//! +//! over the unit simplex `C = { w >= 0, 1^T w = 1 }`. We solve this LP by +//! a projected sub-gradient method on the simplex; the inner `zeta` is +//! eliminated by setting `zeta = VaR_alpha(L(w))` at every iteration. + +use crate::core::{OptimizrError, Result}; +use ndarray::{Array1, Array2, ArrayView2}; + +use super::check_alpha; + +/// Empirical CVaR of a loss sample at confidence level `alpha`. +pub fn cvar_value(losses: &[f64], alpha: f64) -> Result { + check_alpha(alpha)?; + if losses.is_empty() { + return Err(OptimizrError::EmptyData); + } + let mut sorted: Vec = losses.to_vec(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let n = sorted.len(); + // CVaR_alpha = mean of the worst (1 - alpha) fraction of losses. + let k = ((alpha * n as f64).floor() as usize).min(n.saturating_sub(1)); + let tail = &sorted[k..]; + Ok(tail.iter().sum::() / tail.len() as f64) +} + +/// Configuration for `minimize_cvar`. +#[derive(Debug, Clone)] +pub struct CVaRConfig { + pub alpha: f64, + pub n_iter: usize, + pub step_size: f64, + pub tol: f64, +} + +impl Default for CVaRConfig { + fn default() -> Self { + Self { + alpha: 0.95, + n_iter: 5_000, + step_size: 1e-2, + tol: 1e-8, + } + } +} + +/// Result of the CVaR minimisation. +#[derive(Debug, Clone)] +pub struct CVaRResult { + pub w: Array1, + pub zeta: f64, + pub cvar: f64, + pub iterations: usize, +} + +/// Minimise empirical CVaR of `L(w) = -` over the unit +/// simplex. +/// +/// `returns` has shape `(S, d)` (S samples, d decision components). +pub fn minimize_cvar(returns: ArrayView2, cfg: &CVaRConfig) -> Result { + check_alpha(cfg.alpha)?; + let s = returns.nrows(); + let d = returns.ncols(); + if s == 0 || d == 0 { + return Err(OptimizrError::EmptyData); + } + if cfg.n_iter == 0 { + return Err(OptimizrError::InvalidParameter("n_iter must be > 0".into())); + } + + let mut w = Array1::::from_elem(d, 1.0 / d as f64); + let mut prev_obj = f64::INFINITY; + let mut iters = 0usize; + + let inv_factor = 1.0 / ((1.0 - cfg.alpha) * s as f64); + + for it in 0..cfg.n_iter { + // Compute losses + let losses: Vec = (0..s) + .map(|i| -returns.row(i).dot(&w)) + .collect(); + + // VaR (alpha-quantile of losses) gives zeta*. + let mut sorted = losses.clone(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let k = ((cfg.alpha * s as f64).ceil() as usize).max(1).min(s); + let zeta = sorted[k - 1]; + + // Objective and sub-gradient with zeta fixed at current optimum. + let mut grad = Array1::::zeros(d); + let mut obj = zeta; + let mut tail_count = 0usize; + for (i, &l) in losses.iter().enumerate() { + let exceed = l - zeta; + if exceed > 0.0 { + obj += inv_factor * exceed; + // d (l - zeta)_+ / d w_j = d l / d w_j = -r_{i, j} + for j in 0..d { + grad[j] -= inv_factor * returns[[i, j]]; + } + tail_count += 1; + } + } + let _ = tail_count; + + // Projected gradient step: descend then project on simplex. + let lr = cfg.step_size / (1.0 + (it as f64).sqrt()); + for j in 0..d { + w[j] -= lr * grad[j]; + } + project_simplex_inplace(&mut w); + + if (prev_obj - obj).abs() < cfg.tol { + iters = it + 1; + return Ok(CVaRResult { + w, + zeta, + cvar: obj, + iterations: iters, + }); + } + prev_obj = obj; + iters = it + 1; + } + + let losses: Vec = (0..s).map(|i| -returns.row(i).dot(&w)).collect(); + let final_cvar = cvar_value(&losses, cfg.alpha)?; + let zeta_final = { + let mut sorted = losses.clone(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let k = ((cfg.alpha * s as f64).ceil() as usize).max(1).min(s); + sorted[k - 1] + }; + Ok(CVaRResult { + w, + zeta: zeta_final, + cvar: final_cvar, + iterations: iters, + }) +} + +/// Project a vector onto the probability simplex `{w >= 0, sum w = 1}`. +/// Algorithm of Held, Wolfe, Crowder (1974). +pub fn project_simplex_inplace(v: &mut Array1) { + let n = v.len(); + if n == 0 { + return; + } + let mut u: Vec = v.iter().copied().collect(); + u.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)); + let mut cssv = 0.0; + let mut rho = 0usize; + for (i, &ui) in u.iter().enumerate() { + cssv += ui; + let t = (cssv - 1.0) / (i as f64 + 1.0); + if ui - t > 0.0 { + rho = i + 1; + } + } + let cssv_rho: f64 = u.iter().take(rho).sum(); + let theta = (cssv_rho - 1.0) / rho as f64; + for x in v.iter_mut() { + *x = (*x - theta).max(0.0); + } +} + +/// Convenience wrapper accepting an owned ndarray. +pub fn minimize_cvar_owned(returns: Array2, cfg: &CVaRConfig) -> Result { + minimize_cvar(returns.view(), cfg) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn cvar_value_matches_tail_mean() { + let losses: Vec = (1..=100).map(|x| x as f64).collect(); + // alpha=0.9: tail = 91..100 => mean = 95.5 + let c = cvar_value(&losses, 0.9).unwrap(); + assert!((c - 95.5).abs() < 1e-12); + } + + #[test] + fn simplex_projection_is_idempotent_on_simplex() { + let mut v = array![0.2_f64, 0.3, 0.5]; + project_simplex_inplace(&mut v); + let s: f64 = v.iter().sum(); + assert!((s - 1.0).abs() < 1e-12); + for &x in v.iter() { + assert!(x >= 0.0); + } + } + + #[test] + fn cvar_min_concentrates_on_safer_decision() { + // d=2 decisions; samples: r0 ~ very volatile, r1 ~ stable. + let mut data = vec![0.0; 200 * 2]; + let mut state = 42u64; + for s in 0..200 { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let u1 = ((state & 0xFFFF) as f64 + 1.0) / 65538.0; + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let u2 = ((state & 0xFFFF) as f64 + 1.0) / 65538.0; + let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos(); + data[s * 2] = 5.0 * z; // volatile + data[s * 2 + 1] = 0.1 * z; // stable + } + let returns = Array2::from_shape_vec((200, 2), data).unwrap(); + let cfg = CVaRConfig { + alpha: 0.95, + n_iter: 2_000, + step_size: 0.05, + tol: 0.0, + }; + let res = minimize_cvar(returns.view(), &cfg).unwrap(); + // Optimiser should put much more weight on the stable component. + assert!(res.w[1] > res.w[0], "weights = {:?}", res.w); + } +} diff --git a/src/risk_measures/mod.rs b/src/risk_measures/mod.rs new file mode 100644 index 0000000..faf17d8 --- /dev/null +++ b/src/risk_measures/mod.rs @@ -0,0 +1,27 @@ +//! Generic empirical risk measures. +//! +//! This module exposes purely sample-based risk functionals over a vector +//! of decision variables `w` and a sample of losses +//! `L^{(s)} = -` (the loss is signed; users pass the sample +//! of returns `r^{(s)}` directly). All algorithms work on arbitrary +//! datasets and make no domain-specific assumption. + +use crate::core::{OptimizrError, Result}; + +pub mod var; +pub mod cvar; + + +pub use cvar::{cvar_value, minimize_cvar, CVaRConfig, CVaRResult}; +pub use var::{historical_var, parametric_var}; + +#[inline] +pub(crate) fn check_alpha(alpha: f64) -> Result<()> { + if !(0.0 < alpha && alpha < 1.0) { + return Err(OptimizrError::InvalidParameter(format!( + "alpha must lie in (0, 1), got {}", + alpha + ))); + } + Ok(()) +} diff --git a/src/risk_measures/var.rs b/src/risk_measures/var.rs new file mode 100644 index 0000000..f997be8 --- /dev/null +++ b/src/risk_measures/var.rs @@ -0,0 +1,120 @@ +//! Value-at-Risk estimators. +//! +//! For a real-valued loss random variable `L` and confidence level +//! `alpha in (0, 1)`, +//! +//! ```text +//! VaR_alpha(L) = inf { x in R : P(L <= x) >= alpha }. +//! ``` +//! +//! Two estimators are provided: +//! +//! * `historical_var` -- `alpha`-quantile of an empirical loss sample. +//! * `parametric_var` -- Gaussian closed form `mu + sigma * Phi^{-1}(alpha)`. + +use crate::core::{OptimizrError, Result}; + +use super::check_alpha; + +/// Empirical (historical) Value-at-Risk at confidence level `alpha`. +/// +/// `losses` is a non-empty sample of realised losses (positive losses +/// = bad outcomes). The function returns the smallest order statistic +/// `L_(k)` with rank `k = ceil(alpha * n)`. +pub fn historical_var(losses: &[f64], alpha: f64) -> Result { + check_alpha(alpha)?; + if losses.is_empty() { + return Err(OptimizrError::EmptyData); + } + let mut sorted: Vec = losses.to_vec(); + sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let n = sorted.len(); + let k = ((alpha * n as f64).ceil() as usize).max(1).min(n); + Ok(sorted[k - 1]) +} + +/// Closed-form Gaussian Value-at-Risk +/// `VaR = mu + sigma * Phi^{-1}(alpha)`. +pub fn parametric_var(mu: f64, sigma: f64, alpha: f64) -> Result { + check_alpha(alpha)?; + if !(sigma >= 0.0) { + return Err(OptimizrError::InvalidParameter( + "sigma must be non-negative".into(), + )); + } + Ok(mu + sigma * inverse_normal_cdf(alpha)) +} + +/// Beasley--Springer--Moro inverse standard normal CDF. +pub(crate) fn inverse_normal_cdf(p: f64) -> f64 { + // Acklam's algorithm (high accuracy). + let a = [ + -3.969683028665376e+01, + 2.209460984245205e+02, + -2.759285104469687e+02, + 1.383577518672690e+02, + -3.066479806614716e+01, + 2.506628277459239e+00, + ]; + let b = [ + -5.447609879822406e+01, + 1.615858368580409e+02, + -1.556989798598866e+02, + 6.680131188771972e+01, + -1.328068155288572e+01, + ]; + let c = [ + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e+00, + -2.549732539343734e+00, + 4.374664141464968e+00, + 2.938163982698783e+00, + ]; + let d = [ + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e+00, + 3.754408661907416e+00, + ]; + let p_low = 0.02425; + let p_high = 1.0 - p_low; + if p < p_low { + let q = (-2.0 * p.ln()).sqrt(); + return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) + / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0); + } + if p <= p_high { + let q = p - 0.5; + let r = q * q; + return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q + / (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0); + } + let q = (-2.0 * (1.0 - p).ln()).sqrt(); + -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) + / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn historical_var_matches_quantile() { + let losses: Vec = (1..=100).map(|x| x as f64).collect(); + let v = historical_var(&losses, 0.95).unwrap(); + assert_eq!(v, 95.0); + } + + #[test] + fn parametric_var_gaussian_95() { + let v = parametric_var(0.0, 1.0, 0.95).unwrap(); + assert!((v - 1.6448536).abs() < 1e-4); + } + + #[test] + fn rejects_bad_alpha() { + assert!(historical_var(&[1.0, 2.0], 0.0).is_err()); + assert!(historical_var(&[1.0, 2.0], 1.0).is_err()); + } +} diff --git a/src/signatures/log_signature.rs b/src/signatures/log_signature.rs new file mode 100644 index 0000000..f24afa3 --- /dev/null +++ b/src/signatures/log_signature.rs @@ -0,0 +1,105 @@ +//! Log-signature in the free Lie algebra. +//! +//! Given the truncated signature `S` of a path, the log-signature +//! +//! ```text +//! log(S) = sum_{n>=1} (-1)^{n+1} / n * (S - 1)^{otimes n} +//! ``` +//! +//! is a Lie series; its coefficients in any Hall basis form a more +//! parsimonious representation of the path. +//! +//! This implementation: +//! +//! * computes the truncated tensor logarithm `log_M(S)` of a +//! `TruncatedSignature` (Lyons-Victoir 2007 eq. 2.13); +//! * stores the result as flat tensors `[L_1, L_2, ..., L_M]` (the +//! level-0 component is identically zero). +//! +//! No projection onto a Hall basis is performed -- the truncated tensor +//! logarithm itself already lives in the truncated free Lie algebra, +//! and consumers can project it onto any preferred basis. + +use crate::core::{OptimizrError, Result}; + +use super::path_signature::TruncatedSignature; + +/// Truncated tensor logarithm. +#[derive(Debug, Clone)] +pub struct TruncatedLogSignature { + pub channels: usize, + pub level: usize, + /// `tensors[k]` has length `channels^k`; `tensors[0]` is empty by + /// convention. + pub tensors: Vec>, +} + +fn truncated_tensor_product(a: &[Vec], b: &[Vec], d: usize, level: usize) -> Vec> { + let mut out: Vec> = (0..=level).map(|k| vec![0.0; d.pow(k as u32)]).collect(); + for n in 0..=level { + for k in 0..=n { + let nk = n - k; + let stride = d.pow(nk as u32); + for alpha in 0..d.pow(k as u32) { + let av = a[k][alpha]; + if av == 0.0 { + continue; + } + let base = alpha * stride; + for beta in 0..stride { + out[n][base + beta] += av * b[nk][beta]; + } + } + } + } + out +} + +/// Compute the truncated tensor logarithm of a signature. +pub fn log_signature(sig: &TruncatedSignature) -> Result { + let d = sig.channels; + let m = sig.level; + if m == 0 { + return Err(OptimizrError::InvalidParameter( + "log-signature requires level >= 1".into(), + )); + } + // Build (S - 1) by removing the level-0 unit element. + let mut s_minus_one: Vec> = sig.tensors.clone(); + s_minus_one[0][0] = 0.0; + // Power series: log(1 + x) = sum_{n>=1} (-1)^{n+1}/n x^n + let mut log_tensors: Vec> = (0..=m).map(|k| vec![0.0; d.pow(k as u32)]).collect(); + let mut x_power = s_minus_one.clone(); + for n in 1..=m { + let coef = if n % 2 == 0 { -1.0 / (n as f64) } else { 1.0 / (n as f64) }; + for k in 0..=m { + for (i, v) in x_power[k].iter().enumerate() { + log_tensors[k][i] += coef * v; + } + } + if n < m { + x_power = truncated_tensor_product(&x_power, &s_minus_one, d, m); + } + } + Ok(TruncatedLogSignature { + channels: d, + level: m, + tensors: log_tensors, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::signatures::path_signature::path_signature; + + #[test] + fn linear_segment_log_signature_first_level_equals_increment() { + let path = vec![vec![0.0, 0.0], vec![0.3, -0.2]]; + let sig = path_signature(&path, 2).unwrap(); + let log = log_signature(&sig).unwrap(); + // The first-level part of the log-signature equals the path increment. + assert!((log.tensors[1][0] - 0.3).abs() < 1e-12); + assert!((log.tensors[1][1] + 0.2).abs() < 1e-12); + } +} diff --git a/src/signatures/mod.rs b/src/signatures/mod.rs new file mode 100644 index 0000000..04cc82b --- /dev/null +++ b/src/signatures/mod.rs @@ -0,0 +1,16 @@ +//! Path signatures and related tensor algebra utilities. +//! +//! Submodules: +//! +//! * `path_signature` -- truncated tensor signature. +//! * `log_signature` -- log-signature in a Hall basis. +//! * `random_signature` -- random projection of the signature. +//! * `signature_kernel` -- Salvi--Cass--Lyons signature kernel PDE. +//! * `utils` -- shuffle product, concatenation helpers. + +pub mod log_signature; +pub mod path_signature; +pub mod random_signature; +pub mod signature_kernel; +pub mod utils; + diff --git a/src/signatures/path_signature.rs b/src/signatures/path_signature.rs new file mode 100644 index 0000000..a501db9 --- /dev/null +++ b/src/signatures/path_signature.rs @@ -0,0 +1,191 @@ +//! Truncated tensor path signature (Lyons 1998). +//! +//! For a continuous path `X : [0, T] -> R^d` of bounded variation, the +//! signature is the formal series +//! +//! ```text +//! S(X)_{0,T} = 1 + sum_{k>=1} sum_{i_1, ..., i_k} S^{i_1, ..., i_k}_{0,T} e_{i_1} otimes ... otimes e_{i_k} +//! ``` +//! +//! with iterated Stieltjes integrals +//! +//! ```text +//! S^{i_1, ..., i_k}_{0,T} = int_{0 < u_1 < ... < u_k < T} dX^{i_1}_{u_1} ... dX^{i_k}_{u_k}. +//! ``` +//! +//! For a piecewise-linear path with increments `Delta_n in R^d`, the +//! signature truncated at level `M` satisfies the recursion +//! +//! ```text +//! S^{(M)}_{0, t_n} = S^{(M)}_{0, t_{n-1}} otimes_M exp_M(Delta_n) +//! ``` +//! +//! where `otimes_M` is the truncated tensor product and `exp_M` the +//! truncated tensor exponential. +//! +//! Storage convention: a signature truncated at level `M` is a vector of +//! flat tensors `S = [s_0, s_1, ..., s_M]` where `s_k in R^{d^k}` is +//! stored row-major with strides `(d^{k-1}, d^{k-2}, ..., 1)`. + +use crate::core::{OptimizrError, Result}; + +/// Signature truncated at level `M`. +#[derive(Debug, Clone)] +pub struct TruncatedSignature { + pub channels: usize, + pub level: usize, + pub tensors: Vec>, // tensors[k] has length channels^k +} + +impl TruncatedSignature { + pub fn identity(channels: usize, level: usize) -> Self { + let mut tensors = Vec::with_capacity(level + 1); + for k in 0..=level { + let len = channels.pow(k as u32); + let mut t = vec![0.0; len]; + if k == 0 { + t[0] = 1.0; + } + tensors.push(t); + } + Self { + channels, + level, + tensors, + } + } +} + +/// Truncated tensor product `(a otimes b)` up to level `M`. +fn truncated_tensor_product(a: &TruncatedSignature, b: &TruncatedSignature) -> TruncatedSignature { + let d = a.channels; + let m = a.level; + let mut out = TruncatedSignature::identity(d, m); + for n in 0..=m { + let len_n = d.pow(n as u32); + for buf in out.tensors[n].iter_mut() { + *buf = 0.0; + } + for k in 0..=n { + let nk = n - k; + let a_k = &a.tensors[k]; + let b_nk = &b.tensors[nk]; + // out[n][i_1...i_n] += a[k][i_1...i_k] * b[nk][i_{k+1}...i_n] + // Linearised: for each pair (alpha in [0, d^k), beta in [0, d^{nk})), out[alpha * d^{nk} + beta] += a[alpha] * b[beta] + let stride = d.pow(nk as u32); + for alpha in 0..d.pow(k as u32) { + let av = a_k[alpha]; + if av == 0.0 { + continue; + } + let base = alpha * stride; + for beta in 0..stride { + out.tensors[n][base + beta] += av * b_nk[beta]; + } + } + let _ = len_n; + } + } + out +} + +/// Truncated tensor exponential of an increment `Delta in R^d`. +/// +/// `exp_M(Delta) = 1 + Delta + Delta^{otimes 2}/2! + ... + Delta^{otimes M}/M!` +fn truncated_exponential(delta: &[f64], level: usize) -> TruncatedSignature { + let d = delta.len(); + let mut sig = TruncatedSignature::identity(d, level); + if level == 0 || d == 0 { + return sig; + } + // Build delta^{otimes k}/k! incrementally. + let mut current = vec![1.0]; // delta^0 + let mut current_dim = 0usize; + let mut factorial = 1.0_f64; + for k in 1..=level { + let new_len = d.pow(k as u32); + let mut new = vec![0.0; new_len]; + let stride = d; + for alpha in 0..current.len() { + let cv = current[alpha]; + if cv == 0.0 { + continue; + } + let base = alpha * stride; + for j in 0..d { + new[base + j] = cv * delta[j]; + } + } + current = new; + current_dim = k; + factorial *= k as f64; + let inv = 1.0 / factorial; + for (i, v) in sig.tensors[k].iter_mut().enumerate() { + *v = current[i] * inv; + } + } + let _ = current_dim; + sig +} + +/// Compute the truncated signature of a piecewise-linear path given by +/// its sample points `path` of shape `(n_steps + 1, d)` (row-major). +pub fn path_signature(path: &[Vec], level: usize) -> Result { + if path.len() < 2 { + return Err(OptimizrError::InvalidInput( + "path must contain at least two points".into(), + )); + } + let d = path[0].len(); + if d == 0 { + return Err(OptimizrError::InvalidInput("zero-dimensional path".into())); + } + for p in path { + if p.len() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: p.len(), + }); + } + } + let mut sig = TruncatedSignature::identity(d, level); + for n in 1..path.len() { + let mut delta = vec![0.0; d]; + for j in 0..d { + delta[j] = path[n][j] - path[n - 1][j]; + } + let inc_sig = truncated_exponential(&delta, level); + sig = truncated_tensor_product(&sig, &inc_sig); + } + Ok(sig) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// For a single linear segment `X(t) = t * Delta`, the signature reads + /// `S^{i_1, ..., i_k} = (Delta^{i_1} ... Delta^{i_k}) / k!`. + #[test] + fn linear_segment_signature_matches_exponential() { + let delta = vec![0.3, -0.2]; + let path = vec![vec![0.0, 0.0], delta.clone()]; + let sig = path_signature(&path, 3).unwrap(); + // level 1: Delta itself + for j in 0..2 { + assert!((sig.tensors[1][j] - delta[j]).abs() < 1e-12); + } + // level 2: Delta_i Delta_j / 2 + for i in 0..2 { + for j in 0..2 { + let exp = delta[i] * delta[j] / 2.0; + assert!((sig.tensors[2][i * 2 + j] - exp).abs() < 1e-12); + } + } + } + + #[test] + fn empty_path_rejected() { + assert!(path_signature(&Vec::>::new(), 2).is_err()); + } +} diff --git a/src/signatures/random_signature.rs b/src/signatures/random_signature.rs new file mode 100644 index 0000000..52f55b5 --- /dev/null +++ b/src/signatures/random_signature.rs @@ -0,0 +1,152 @@ +//! Random projection of the path signature. +//! +//! Following Cuchiero, Schmocker & Teichmann (2023), one defines a +//! random reservoir on `R^N` driven by the controlled SDE +//! +//! ```text +//! dZ_t = A_0 Z_t dt + sum_{i=1}^d A_i Z_t dX^i_t, Z_0 in R^N, +//! ``` +//! +//! where the matrices `A_0, ..., A_d in R^{N x N}` are drawn from a +//! random ensemble (typically i.i.d. Gaussian entries with variance +//! `1 / N`). The map `X -> Z_T` provides a finite-dimensional, randomly +//! projected representation of the signature. +//! +//! This implementation integrates the controlled equation with an +//! Euler--Maruyama-type scheme on the supplied path samples. + +use crate::core::{OptimizrError, Result}; +use rand::SeedableRng; +use rand_distr::{Distribution, Normal}; + +/// Random reservoir parameters. +#[derive(Debug, Clone)] +pub struct RandomSignatureConfig { + pub reservoir_dim: usize, + pub seed: u64, + pub variance: f64, +} + +impl Default for RandomSignatureConfig { + fn default() -> Self { + Self { + reservoir_dim: 32, + seed: 0, + variance: 1.0, + } + } +} + +/// Result of a random reservoir integration. +#[derive(Debug, Clone)] +pub struct RandomSignatureResult { + pub trajectory: Vec>, // length = path.len(), each row in R^N +} + +fn matvec(a: &[f64], n: usize, x: &[f64]) -> Vec { + let mut y = vec![0.0; n]; + for i in 0..n { + let mut s = 0.0; + for j in 0..n { + s += a[i * n + j] * x[j]; + } + y[i] = s; + } + y +} + +/// Compute `Z_T` for the controlled SDE driven by the piecewise-linear +/// path interpolating `path`. Returns the full reservoir trajectory. +pub fn random_signature( + path: &[Vec], + cfg: &RandomSignatureConfig, +) -> Result { + if path.len() < 2 { + return Err(OptimizrError::InvalidInput( + "path must contain at least two points".into(), + )); + } + let d = path[0].len(); + if d == 0 { + return Err(OptimizrError::InvalidInput("zero-dimensional path".into())); + } + for p in path { + if p.len() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: p.len(), + }); + } + } + if cfg.reservoir_dim == 0 { + return Err(OptimizrError::InvalidParameter( + "reservoir_dim > 0 required".into(), + )); + } + if cfg.variance <= 0.0 { + return Err(OptimizrError::InvalidParameter( + "variance > 0 required".into(), + )); + } + let n = cfg.reservoir_dim; + let std_dev = (cfg.variance / n as f64).sqrt(); + let normal = Normal::new(0.0, std_dev) + .map_err(|e| OptimizrError::NumericalError(format!("normal init: {}", e)))?; + let mut rng = rand::rngs::StdRng::seed_from_u64(cfg.seed); + + // Sample matrices A_0, ..., A_d. + let mut matrices: Vec> = (0..=d) + .map(|_| (0..n * n).map(|_| normal.sample(&mut rng)).collect()) + .collect(); + // Drift A_0 typically scaled smaller; here we keep the same variance. + let _ = &mut matrices; + + let mut z = vec![0.0; n]; + z[0] = 1.0; + let mut traj = Vec::with_capacity(path.len()); + traj.push(z.clone()); + let dt = 1.0 / (path.len() - 1) as f64; + for n_idx in 1..path.len() { + let mut delta = vec![0.0; d]; + for j in 0..d { + delta[j] = path[n_idx][j] - path[n_idx - 1][j]; + } + // Z_{n+1} = Z_n + A_0 Z_n dt + sum_i A_i Z_n delta_i + let drift = matvec(&matrices[0], n, &z); + let mut new_z = z.clone(); + for i in 0..n { + new_z[i] += dt * drift[i]; + } + for i in 0..d { + let inc = matvec(&matrices[i + 1], n, &z); + for j in 0..n { + new_z[j] += inc[j] * delta[i]; + } + } + z = new_z; + traj.push(z.clone()); + } + Ok(RandomSignatureResult { trajectory: traj }) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn deterministic_seed_reproduces_output() { + let path = vec![vec![0.0], vec![0.5], vec![1.0]]; + let cfg = RandomSignatureConfig { + reservoir_dim: 4, + seed: 42, + variance: 1.0, + }; + let r1 = random_signature(&path, &cfg).unwrap(); + let r2 = random_signature(&path, &cfg).unwrap(); + for (a, b) in r1.trajectory.iter().zip(r2.trajectory.iter()) { + for (x, y) in a.iter().zip(b.iter()) { + assert!((x - y).abs() < 1e-15); + } + } + } +} diff --git a/src/signatures/signature_kernel.rs b/src/signatures/signature_kernel.rs new file mode 100644 index 0000000..c8ae7de --- /dev/null +++ b/src/signatures/signature_kernel.rs @@ -0,0 +1,107 @@ +//! Signature kernel via the Salvi--Cass--Lyons PDE. +//! +//! For two paths `X : [0, S] -> R^d` and `Y : [0, T] -> R^d`, the +//! signature kernel +//! +//! ```text +//! K(s, t) = < S(X)_{0, s} , S(Y)_{0, t} > +//! ``` +//! +//! satisfies the linear hyperbolic PDE +//! +//! ```text +//! d^2 K / (ds dt) = < dX_s, dY_t > K(s, t), K(s, 0) = K(0, t) = 1. +//! ``` +//! +//! This module integrates the PDE on a uniform grid via the Goursat +//! finite-difference scheme of Salvi et al. (2021): +//! +//! ```text +//! K_{i+1, j+1} = K_{i+1, j} + K_{i, j+1} - K_{i, j} +//! + (Delta x_i . Delta y_j) * 0.5 * (K_{i+1, j} + K_{i, j+1}) +//! ``` +//! +//! and returns the value `K(S, T)`. + +use crate::core::{OptimizrError, Result}; + +/// Signature kernel evaluation result. +#[derive(Debug, Clone)] +pub struct SignatureKernelResult { + pub value: f64, + pub grid: Vec>, // K_{i, j}, shape (n_x, n_y) +} + +/// Compute the signature kernel between piecewise-linear paths `x` and +/// `y` (both shape `(n, d)` with arbitrary `n`). +pub fn signature_kernel(x: &[Vec], y: &[Vec]) -> Result { + if x.len() < 2 || y.len() < 2 { + return Err(OptimizrError::InvalidInput( + "both paths must have at least two points".into(), + )); + } + let d = x[0].len(); + if d == 0 { + return Err(OptimizrError::InvalidInput("zero-dimensional path".into())); + } + for p in x.iter().chain(y.iter()) { + if p.len() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: p.len(), + }); + } + } + let nx = x.len(); + let ny = y.len(); + // Pre-compute increments + let dx: Vec> = (0..nx - 1) + .map(|i| (0..d).map(|k| x[i + 1][k] - x[i][k]).collect()) + .collect(); + let dy: Vec> = (0..ny - 1) + .map(|j| (0..d).map(|k| y[j + 1][k] - y[j][k]).collect()) + .collect(); + + let mut k_grid = vec![vec![0.0; ny]; nx]; + for i in 0..nx { + k_grid[i][0] = 1.0; + } + for j in 0..ny { + k_grid[0][j] = 1.0; + } + for i in 0..nx - 1 { + for j in 0..ny - 1 { + let dot: f64 = (0..d).map(|c| dx[i][c] * dy[j][c]).sum(); + let avg = 0.5 * (k_grid[i + 1][j] + k_grid[i][j + 1]); + k_grid[i + 1][j + 1] = k_grid[i + 1][j] + k_grid[i][j + 1] - k_grid[i][j] + dot * avg; + } + } + Ok(SignatureKernelResult { + value: k_grid[nx - 1][ny - 1], + grid: k_grid, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Identical linear segments yield K(1, 1) = exp(||Delta||^2) (closed + /// form for the signature inner product of a single segment with + /// itself). + #[test] + fn identical_linear_paths_give_exponential_kernel() { + let delta = vec![0.4, -0.1]; + let n = 200; + let path: Vec> = (0..=n) + .map(|k| { + let t = k as f64 / n as f64; + vec![t * delta[0], t * delta[1]] + }) + .collect(); + let res = signature_kernel(&path, &path).unwrap(); + let exact = (delta.iter().map(|v| v * v).sum::()).exp(); + let err = (res.value - exact).abs(); + assert!(err < 1e-2, "K = {} expected {} err = {}", res.value, exact, err); + } +} diff --git a/src/signatures/utils.rs b/src/signatures/utils.rs new file mode 100644 index 0000000..1f0f4db --- /dev/null +++ b/src/signatures/utils.rs @@ -0,0 +1,90 @@ +//! Helper algebraic operations on flat tensors of the truncated tensor +//! algebra `T^M((R^d))`. +//! +//! Provides: +//! +//! * `shuffle_product` -- shuffle product of two flat words. +//! * `concatenate_signatures` -- Chen identity-driven concatenation. + +use super::path_signature::TruncatedSignature; +use crate::core::Result; + +/// Shuffle product of two words `u, v` of length `p, q` over the +/// alphabet `{0, ..., d - 1}`. +/// +/// Returned as a `HashMap, f64>` keyed by word. +pub fn shuffle_product( + u: &[usize], + v: &[usize], +) -> std::collections::HashMap, f64> { + let mut acc = std::collections::HashMap::new(); + fn rec( + u: &[usize], + v: &[usize], + prefix: &mut Vec, + acc: &mut std::collections::HashMap, f64>, + ) { + if u.is_empty() && v.is_empty() { + *acc.entry(prefix.clone()).or_insert(0.0) += 1.0; + return; + } + if !u.is_empty() { + prefix.push(u[0]); + rec(&u[1..], v, prefix, acc); + prefix.pop(); + } + if !v.is_empty() { + prefix.push(v[0]); + rec(u, &v[1..], prefix, acc); + prefix.pop(); + } + } + let mut prefix = Vec::new(); + rec(u, v, &mut prefix, &mut acc); + acc +} + +/// Concatenate two signatures via Chen's identity: `S(X * Y) = S(X) otimes S(Y)`. +pub fn concatenate_signatures( + a: &TruncatedSignature, + b: &TruncatedSignature, +) -> Result { + assert_eq!(a.channels, b.channels); + assert_eq!(a.level, b.level); + let d = a.channels; + let m = a.level; + let mut out = TruncatedSignature::identity(d, m); + for n in 0..=m { + for buf in out.tensors[n].iter_mut() { + *buf = 0.0; + } + for k in 0..=n { + let nk = n - k; + let stride = d.pow(nk as u32); + for alpha in 0..d.pow(k as u32) { + let av = a.tensors[k][alpha]; + if av == 0.0 { + continue; + } + let base = alpha * stride; + for beta in 0..stride { + out.tensors[n][base + beta] += av * b.tensors[nk][beta]; + } + } + } + } + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn shuffle_one_letter_words() { + let res = shuffle_product(&[0], &[1]); + assert_eq!(res.len(), 2); + assert!((res[&vec![0, 1]] - 1.0).abs() < 1e-12); + assert!((res[&vec![1, 0]] - 1.0).abs() < 1e-12); + } +} diff --git a/src/timeseries_utils.rs b/src/timeseries_utils.rs index 53dc254..8d8b61b 100644 --- a/src/timeseries_utils.rs +++ b/src/timeseries_utils.rs @@ -31,6 +31,9 @@ use ndarray::Array1; #[cfg(feature = "python-bindings")] pub mod python_bindings; +pub mod nonsync_covariance; +pub mod wavelet; + /// Prepare time-series price data for HMM regime detection /// /// Creates features from price series including: diff --git a/src/timeseries_utils/nonsync_covariance.rs b/src/timeseries_utils/nonsync_covariance.rs new file mode 100644 index 0000000..59b2aa3 --- /dev/null +++ b/src/timeseries_utils/nonsync_covariance.rs @@ -0,0 +1,254 @@ +//! Hayashi--Yoshida non-synchronous covariance estimator. +//! +//! Given two real-valued processes observed on heterogeneous time grids +//! `t^{(1)}_0 < ... < t^{(1)}_{n_1}` and `t^{(2)}_0 < ... < t^{(2)}_{n_2}`, +//! the Hayashi--Yoshida estimator is +//! +//! ```text +//! Sigma_{1,2}^{HY} = sum_k sum_l Dr1_k * Dr2_l * 1{ I1_k inter I2_l != empty } +//! ``` +//! +//! where `I_i^{(k)} = ( t_{i,k-1}, t_{i,k} ]` is the k-th native +//! observation interval of series i and `Dr_i^{(k)}` the corresponding +//! increment. +//! +//! The implementation is `O(n_1 + n_2)` thanks to the monotonicity of the +//! grids: a two-pointer sweep finds the first index `l` of grid 2 whose +//! interval intersects the current interval of grid 1 and stops as soon as +//! the intervals separate. +//! +//! # Reference +//! +//! Hayashi & Yoshida (2005), *On covariance estimation of non-synchronously +//! observed diffusion processes*, Bernoulli 11(2). + +use rayon::prelude::*; + +use crate::core::{OptimizrError, Result}; + +/// Compute the Hayashi--Yoshida covariance between two non-synchronous +/// observation series. +/// +/// `t1` (resp. `t2`) is a strictly increasing time grid of length +/// `v1.len()` (resp. `v2.len()`). The increments used by the estimator +/// are `dv_i^{(k)} = v_i[k] - v_i[k-1]` over `(t_i[k-1], t_i[k]]`. +pub fn hayashi_yoshida_covariance( + t1: &[f64], + v1: &[f64], + t2: &[f64], + v2: &[f64], +) -> Result { + if t1.len() != v1.len() || t2.len() != v2.len() { + return Err(OptimizrError::InvalidInput( + "time grid and value series lengths differ".into(), + )); + } + if t1.len() < 2 || t2.len() < 2 { + return Err(OptimizrError::InvalidInput( + "each series needs at least two observations".into(), + )); + } + for w in t1.windows(2) { + if !(w[0] < w[1]) { + return Err(OptimizrError::InvalidInput( + "t1 must be strictly increasing".into(), + )); + } + } + for w in t2.windows(2) { + if !(w[0] < w[1]) { + return Err(OptimizrError::InvalidInput( + "t2 must be strictly increasing".into(), + )); + } + } + + let n1 = t1.len() - 1; + let n2 = t2.len() - 1; + + // For each interval of series 1, sweep series 2 and accumulate. + // Use a sequential search start index that only moves forward, so the + // cost is O(n1 + n2). Wrapped in Rayon when n1*n2 > 1e5 for parallel + // accumulation; the sweep is preserved per worker via local pointer. + let parallel = (n1 as u64) * (n2 as u64) > 100_000; + + let intervals1: Vec<(f64, f64, f64)> = (0..n1) + .map(|k| (t1[k], t1[k + 1], v1[k + 1] - v1[k])) + .collect(); + let intervals2: Vec<(f64, f64, f64)> = (0..n2) + .map(|k| (t2[k], t2[k + 1], v2[k + 1] - v2[k])) + .collect(); + + let acc = if parallel { + intervals1 + .par_iter() + .map(|&(a1, b1, dr1)| { + // Binary search for first interval of series 2 that ends + // strictly above a1 (left-open right-closed convention). + let start = + intervals2.partition_point(|&(_, b2, _)| b2 <= a1); + let mut s = 0.0; + for &(a2, b2, dr2) in intervals2[start..].iter() { + if a2 >= b1 { + break; + } + let _ = b2; + s += dr1 * dr2; + } + s + }) + .sum::() + } else { + let mut s = 0.0; + let mut start = 0usize; + for &(a1, b1, dr1) in intervals1.iter() { + // advance start until b2 > a1 + while start < n2 && intervals2[start].1 <= a1 { + start += 1; + } + let mut l = start; + while l < n2 { + let (a2, _b2, dr2) = intervals2[l]; + if a2 >= b1 { + break; + } + s += dr1 * dr2; + l += 1; + } + } + s + }; + + Ok(acc) +} + +/// Build the full Hayashi--Yoshida covariance matrix for `d` series. +/// +/// `series` is a slice of `(time_grid, values)` pairs. Returns a +/// `d x d` symmetric matrix flattened in row-major order. +pub fn hayashi_yoshida_matrix( + series: &[(Vec, Vec)], +) -> Result>> { + let d = series.len(); + if d == 0 { + return Err(OptimizrError::EmptyData); + } + let mut out = vec![vec![0.0; d]; d]; + for i in 0..d { + for j in i..d { + let cij = hayashi_yoshida_covariance( + &series[i].0, + &series[i].1, + &series[j].0, + &series[j].1, + )?; + out[i][j] = cij; + out[j][i] = cij; + } + } + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn synchronous_grid_matches_classical_covariance() { + let n = 1000; + let t: Vec = (0..=n).map(|k| k as f64).collect(); + let mut x = vec![0.0_f64; n + 1]; + let mut y = vec![0.0_f64; n + 1]; + // X = W, Y = W + N(0) perfectly correlated increments + let mut rng_state = 0x9E37_79B9_7F4A_7C15u64; + for k in 0..n { + // simple xorshift gauss approximation via Box-Muller + rng_state ^= rng_state << 13; + rng_state ^= rng_state >> 7; + rng_state ^= rng_state << 17; + let u1 = ((rng_state & 0xFFFF_FFFF) as f64 + 1.0) / (u32::MAX as f64 + 2.0); + rng_state ^= rng_state << 13; + rng_state ^= rng_state >> 7; + rng_state ^= rng_state << 17; + let u2 = ((rng_state & 0xFFFF_FFFF) as f64 + 1.0) / (u32::MAX as f64 + 2.0); + let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos(); + x[k + 1] = x[k] + z; + y[k + 1] = x[k + 1]; // perfect correlation + } + let cov = hayashi_yoshida_covariance(&t, &x, &t, &y).unwrap(); + let var_x = hayashi_yoshida_covariance(&t, &x, &t, &x).unwrap(); + let var_y = hayashi_yoshida_covariance(&t, &y, &t, &y).unwrap(); + let rho = cov / (var_x.sqrt() * var_y.sqrt()); + assert!((rho - 1.0).abs() < 1e-12, "rho = {}", rho); + } + + #[test] + fn nonsync_subsampling_recovers_correlation() { + // Generate correlated bivariate Brownian increments, then subsample + // each series independently. HY estimator should recover rho. + let n = 20_000; + let dt = 1.0 / n as f64; + let rho_true = 0.6_f64; + let t_full: Vec = (0..=n).map(|k| k as f64 * dt).collect(); + let mut x = vec![0.0_f64; n + 1]; + let mut y = vec![0.0_f64; n + 1]; + + let mut state = 0xDEAD_BEEF_CAFE_F00Du64; + let mut gauss = || { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let u1 = ((state & 0xFFFF_FFFF) as f64 + 1.0) / (u32::MAX as f64 + 2.0); + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let u2 = ((state & 0xFFFF_FFFF) as f64 + 1.0) / (u32::MAX as f64 + 2.0); + (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos() + }; + + let sqrt_dt = dt.sqrt(); + for k in 0..n { + let z1 = gauss(); + let z2 = rho_true * z1 + (1.0 - rho_true * rho_true).sqrt() * gauss(); + x[k + 1] = x[k] + sqrt_dt * z1; + y[k + 1] = y[k] + sqrt_dt * z2; + } + + // Subsample with Bernoulli(0.5) keeping endpoints + let mut keep_x = vec![true; n + 1]; + let mut keep_y = vec![true; n + 1]; + for k in 1..n { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + keep_x[k] = (state & 1) == 0; + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + keep_y[k] = (state & 1) == 0; + } + let collect = |keep: &[bool], v: &[f64]| -> (Vec, Vec) { + let mut tt = Vec::new(); + let mut vv = Vec::new(); + for k in 0..=n { + if keep[k] { + tt.push(t_full[k]); + vv.push(v[k]); + } + } + (tt, vv) + }; + let (t_x, x_s) = collect(&keep_x, &x); + let (t_y, y_s) = collect(&keep_y, &y); + + let cov = hayashi_yoshida_covariance(&t_x, &x_s, &t_y, &y_s).unwrap(); + let var_x = hayashi_yoshida_covariance(&t_x, &x_s, &t_x, &x_s).unwrap(); + let var_y = hayashi_yoshida_covariance(&t_y, &y_s, &t_y, &y_s).unwrap(); + let rho = cov / (var_x.sqrt() * var_y.sqrt()); + assert!( + (rho - rho_true).abs() < 0.05, + "estimated rho = {}", + rho + ); + } +} diff --git a/src/timeseries_utils/wavelet.rs b/src/timeseries_utils/wavelet.rs new file mode 100644 index 0000000..9b4863b --- /dev/null +++ b/src/timeseries_utils/wavelet.rs @@ -0,0 +1,216 @@ +//! Discrete Wavelet Transform (DWT) and Maximal-Overlap DWT (MODWT). +//! +//! Pyramid algorithm with orthogonal mirror filter pairs (`h`, `g`). +//! For a level `j`, +//! +//! ```text +//! W_{j,k} = sum_n h_n * V_{j-1, 2k - n} +//! V_{j,k} = sum_n g_n * V_{j-1, 2k - n} +//! ``` +//! +//! `MODWT` (a.k.a. translation-invariant wavelet transform) skips +//! down-sampling and rescales the filters by `2^{-j/2}`: +//! +//! ```text +//! W~_{j,t} = sum_n (h_n / sqrt(2^j)) * V~_{j-1, t - 2^{j-1} n} +//! ``` +//! +//! # Reference +//! +//! Daubechies (1992), *Ten Lectures on Wavelets*. Percival & Walden +//! (2000), *Wavelet Methods for Time Series Analysis*. + +use crate::core::{OptimizrError, Result}; + +/// Orthogonal scaling-filter families. +#[derive(Debug, Clone, Copy)] +pub enum WaveletFamily { + Haar, + /// Daubechies-N: `n` even integer in {2,4,6,8,10}. + Daubechies(u8), +} + +/// Returns the (decomposition) low-pass scaling filter `g_n` for the +/// requested family. The high-pass wavelet filter is obtained by the QMF +/// relation `h_n = (-1)^n g_{L-1-n}`. +pub fn scaling_filter(family: WaveletFamily) -> Result> { + let g = match family { + WaveletFamily::Haar => vec![std::f64::consts::FRAC_1_SQRT_2; 2], + WaveletFamily::Daubechies(n) => match n { + 2 => vec![std::f64::consts::FRAC_1_SQRT_2; 2], + 4 => { + let s3 = 3.0_f64.sqrt(); + let denom = 4.0 * 2.0_f64.sqrt(); + vec![ + (1.0 + s3) / denom, + (3.0 + s3) / denom, + (3.0 - s3) / denom, + (1.0 - s3) / denom, + ] + } + 6 => vec![ + 0.332_670_552_950_082_6, + 0.806_891_509_311_092_5, + 0.459_877_502_118_491_7, + -0.135_011_020_010_254_6, + -0.085_441_273_882_026_7, + 0.035_226_291_882_100_7, + ], + 8 => vec![ + 0.230_377_813_308_896_5, + 0.714_846_570_552_915_6, + 0.630_880_767_929_858_9, + -0.027_983_769_416_859_8, + -0.187_034_811_719_092_3, + 0.030_841_381_835_560_8, + 0.032_883_011_666_885_2, + -0.010_597_401_785_069_0, + ], + 10 => vec![ + 0.160_102_397_974_193_0, + 0.603_829_269_797_473_5, + 0.724_308_528_438_574_0, + 0.138_428_145_901_320_3, + -0.242_294_887_066_382_4, + -0.032_244_869_585_030_3, + 0.077_571_493_840_065_1, + -0.006_241_490_212_798_3, + -0.012_580_751_999_082_0, + 0.003_335_725_285_473_8, + ], + other => { + return Err(OptimizrError::InvalidParameter(format!( + "Daubechies({}) not supported (use 2,4,6,8,10)", + other + ))) + } + }, + }; + Ok(g) +} + +#[inline] +fn qmf(g: &[f64]) -> Vec { + let l = g.len(); + let mut h = vec![0.0; l]; + for n in 0..l { + let sign = if n % 2 == 0 { 1.0 } else { -1.0 }; + h[n] = sign * g[l - 1 - n]; + } + h +} + +/// One-level DWT step with periodic boundary extension. +/// +/// Returns `(approximation, detail)` each of length `signal.len() / 2`. +pub fn dwt_step(signal: &[f64], family: WaveletFamily) -> Result<(Vec, Vec)> { + let n = signal.len(); + if n < 2 || n % 2 != 0 { + return Err(OptimizrError::InvalidInput( + "signal length must be even and >= 2".into(), + )); + } + let g = scaling_filter(family)?; + let h = qmf(&g); + let l = g.len(); + let half = n / 2; + let mut approx = vec![0.0; half]; + let mut detail = vec![0.0; half]; + for k in 0..half { + let mut a = 0.0; + let mut d = 0.0; + for j in 0..l { + // periodic boundary + let idx = ((2 * k + l - 1 - j) % n + n) % n; + a += g[j] * signal[idx]; + d += h[j] * signal[idx]; + } + approx[k] = a; + detail[k] = d; + } + Ok((approx, detail)) +} + +/// Multi-level pyramid DWT. Returns `(approx_J, [d_1, ..., d_J])`. +pub fn dwt(signal: &[f64], family: WaveletFamily, n_levels: usize) -> Result<(Vec, Vec>)> { + let mut approx = signal.to_vec(); + let mut details = Vec::with_capacity(n_levels); + for _ in 0..n_levels { + let (a, d) = dwt_step(&approx, family)?; + details.push(d); + approx = a; + } + Ok((approx, details)) +} + +/// Maximal-overlap DWT step (no down-sampling, dilated filter at level `j`). +pub fn modwt_step( + signal: &[f64], + family: WaveletFamily, + level: usize, +) -> Result<(Vec, Vec)> { + let n = signal.len(); + if n == 0 { + return Err(OptimizrError::EmptyData); + } + if level == 0 { + return Err(OptimizrError::InvalidParameter("level must be >= 1".into())); + } + let g = scaling_filter(family)?; + let h = qmf(&g); + let scale = (2.0_f64).powf(level as f64 / 2.0); + let g_tilde: Vec = g.iter().map(|x| x / scale).collect(); + let h_tilde: Vec = h.iter().map(|x| x / scale).collect(); + let stride = 1usize << (level - 1); // 2^{j-1} + let l = g.len(); + + let mut approx = vec![0.0; n]; + let mut detail = vec![0.0; n]; + for t in 0..n { + let mut a = 0.0; + let mut d = 0.0; + for j in 0..l { + let idx = (t + n - (j * stride) % n) % n; + a += g_tilde[j] * signal[idx]; + d += h_tilde[j] * signal[idx]; + } + approx[t] = a; + detail[t] = d; + } + Ok((approx, detail)) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn haar_step_constant_signal_zero_detail() { + let signal = vec![1.0_f64; 8]; + let (a, d) = dwt_step(&signal, WaveletFamily::Haar).unwrap(); + assert_eq!(a.len(), 4); + assert_eq!(d.len(), 4); + for v in d.iter() { + assert!(v.abs() < 1e-12, "detail should be zero, got {}", v); + } + let expected_a = std::f64::consts::SQRT_2; + for v in a.iter() { + assert!((v - expected_a).abs() < 1e-12); + } + } + + #[test] + fn db4_filters_have_unit_norm() { + let g = scaling_filter(WaveletFamily::Daubechies(4)).unwrap(); + let s: f64 = g.iter().map(|x| x * x).sum(); + assert!((s - 1.0).abs() < 1e-10); + } + + #[test] + fn modwt_runs_on_short_signal() { + let signal: Vec = (0..16).map(|k| (k as f64).sin()).collect(); + let (a, d) = modwt_step(&signal, WaveletFamily::Daubechies(4), 1).unwrap(); + assert_eq!(a.len(), 16); + assert_eq!(d.len(), 16); + } +} diff --git a/src/topology/bottleneck.rs b/src/topology/bottleneck.rs new file mode 100644 index 0000000..4e53448 --- /dev/null +++ b/src/topology/bottleneck.rs @@ -0,0 +1,226 @@ +//! Bottleneck distance between two persistence diagrams. +//! +//! For diagrams `D, D'` the bottleneck distance is +//! +//! ```text +//! d_B(D, D') = inf_{eta : D -> D'} sup_{x in D} || x - eta(x) ||_inf +//! ``` +//! +//! where the matchings allow points to be paired with the diagonal +//! `Delta = { (t, t) : t >= 0 }` at cost `(d - b) / 2`. +//! +//! Algorithm: binary search on `epsilon`, build the bipartite graph with +//! edges `(p, q)` whenever `||p - q||_inf <= epsilon` (including +//! diagonal projections), test for a perfect matching by Hopcroft--Karp. + +use std::collections::VecDeque; + +use super::persistent_homology::PersistencePair; +use crate::core::Result; + +#[inline] +fn linf_dist(a: (f64, f64), b: (f64, f64)) -> f64 { + (a.0 - b.0).abs().max((a.1 - b.1).abs()) +} + +#[inline] +fn diag_cost(p: (f64, f64)) -> f64 { + (p.1 - p.0).abs() / 2.0 +} + +/// Compute the bottleneck distance between two diagrams. Only pairs with +/// matching `dim` are coupled. +pub fn bottleneck_distance(d1: &[PersistencePair], d2: &[PersistencePair]) -> Result { + let dims: std::collections::BTreeSet = + d1.iter().chain(d2.iter()).map(|p| p.dim).collect(); + + let mut max_d = 0.0_f64; + for dim in dims { + let p1: Vec<(f64, f64)> = d1 + .iter() + .filter(|p| p.dim == dim) + .map(|p| (p.birth, p.death)) + .collect(); + let p2: Vec<(f64, f64)> = d2 + .iter() + .filter(|p| p.dim == dim) + .map(|p| (p.birth, p.death)) + .collect(); + let d = bottleneck_one_dim(&p1, &p2); + if d > max_d { + max_d = d; + } + } + Ok(max_d) +} + +fn bottleneck_one_dim(d1: &[(f64, f64)], d2: &[(f64, f64)]) -> f64 { + // Filter out infinite-death points for finite matching. If counts of + // essential classes differ, the distance is infinite. + let inf1 = d1.iter().filter(|p| p.1.is_infinite()).count(); + let inf2 = d2.iter().filter(|p| p.1.is_infinite()).count(); + if inf1 != inf2 { + return f64::INFINITY; + } + let f1: Vec<(f64, f64)> = d1.iter().copied().filter(|p| p.1.is_finite()).collect(); + let f2: Vec<(f64, f64)> = d2.iter().copied().filter(|p| p.1.is_finite()).collect(); + + // Candidate epsilons: pairwise linf distances + diagonal costs. + let mut cand: Vec = Vec::new(); + for &p in f1.iter() { + for &q in f2.iter() { + cand.push(linf_dist(p, q)); + } + cand.push(diag_cost(p)); + } + for &q in f2.iter() { + cand.push(diag_cost(q)); + } + cand.push(0.0); + cand.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + cand.dedup_by(|a, b| (*a - *b).abs() < 1e-15); + + // Binary search the smallest epsilon admitting a perfect matching. + let mut lo = 0usize; + let mut hi = cand.len() - 1; + while lo < hi { + let mid = (lo + hi) / 2; + if perfect_matching_exists(&f1, &f2, cand[mid]) { + hi = mid; + } else { + lo = mid + 1; + } + } + cand[lo] +} + +/// Build bipartite graph at threshold `eps` and return whether a perfect +/// matching covering both sides exists. Diagonal nodes are added so that +/// every persistence pair can be matched with the diagonal. +fn perfect_matching_exists(left: &[(f64, f64)], right: &[(f64, f64)], eps: f64) -> bool { + let n_l = left.len(); + let n_r = right.len(); + // Augment both sides with diagonal duplicates so that all pairs can + // be matched onto the diagonal at cost diag_cost. + let total_left = n_l + n_r; + let total_right = n_r + n_l; + let mut adj: Vec> = vec![Vec::new(); total_left]; + for i in 0..n_l { + // real -> real + for j in 0..n_r { + if linf_dist(left[i], right[j]) <= eps { + adj[i].push(j); + } + } + // real -> diagonal copy of real left[i] -> right diag index n_r + i + if diag_cost(left[i]) <= eps { + adj[i].push(n_r + i); + } + } + for j in 0..n_r { + // diag copy left[n_l + j] (representing right[j]'s diagonal) only + // matches its own diagonal target n_r + (n_l + j) ... we keep + // simple: it can match any diagonal target at cost zero. + for k in n_r..total_right { + adj[n_l + j].push(k); + } + } + + // Hopcroft--Karp + let nil = usize::MAX; + let mut pair_l = vec![nil; total_left]; + let mut pair_r = vec![nil; total_right]; + let mut dist = vec![0i64; total_left]; + + fn bfs( + adj: &[Vec], + pair_l: &[usize], + pair_r: &[usize], + dist: &mut [i64], + n_l: usize, + nil: usize, + ) -> bool { + let mut q: VecDeque = VecDeque::new(); + for u in 0..n_l { + if pair_l[u] == nil { + dist[u] = 0; + q.push_back(u); + } else { + dist[u] = i64::MAX; + } + } + let mut found = false; + while let Some(u) = q.pop_front() { + for &v in &adj[u] { + let pv = pair_r[v]; + if pv == nil { + found = true; + } else if dist[pv] == i64::MAX { + dist[pv] = dist[u] + 1; + q.push_back(pv); + } + } + } + found + } + + fn dfs( + u: usize, + adj: &[Vec], + pair_l: &mut [usize], + pair_r: &mut [usize], + dist: &mut [i64], + nil: usize, + ) -> bool { + for &v in &adj[u] { + let pv = pair_r[v]; + let cond = pv == nil || (dist[pv] == dist[u] + 1 && dfs(pv, adj, pair_l, pair_r, dist, nil)); + if cond { + pair_l[u] = v; + pair_r[v] = u; + return true; + } + } + dist[u] = i64::MAX; + false + } + + let mut matched = 0usize; + while bfs(&adj, &pair_l, &pair_r, &mut dist, total_left, nil) { + for u in 0..total_left { + if pair_l[u] == nil && dfs(u, &adj, &mut pair_l, &mut pair_r, &mut dist, nil) { + matched += 1; + } + } + } + matched == total_left +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn identical_diagrams_zero_distance() { + let d1 = vec![PersistencePair { + dim: 1, + birth: 0.2, + death: 0.7, + }]; + let d2 = d1.clone(); + let d = bottleneck_distance(&d1, &d2).unwrap(); + assert!(d.abs() < 1e-12); + } + + #[test] + fn point_vs_diagonal_distance() { + let d1 = vec![PersistencePair { + dim: 0, + birth: 0.0, + death: 0.6, + }]; + let d2: Vec = vec![]; + let d = bottleneck_distance(&d1, &d2).unwrap(); + assert!((d - 0.3).abs() < 1e-9, "expected 0.3, got {}", d); + } +} diff --git a/src/topology/mod.rs b/src/topology/mod.rs new file mode 100644 index 0000000..db568f3 --- /dev/null +++ b/src/topology/mod.rs @@ -0,0 +1,16 @@ +//! Topological data analysis utilities. +//! +//! * `persistent_homology` -- Vietoris--Rips persistence via the standard +//! matrix-reduction algorithm (Edelsbrunner & Harer, 2010). +//! * `bottleneck` -- bottleneck distance between two persistence +//! diagrams using the standard `O(n^{2.5})` matching (Hopcroft--Karp on +//! thresholded bipartite graph + binary search on epsilon). + +pub mod bottleneck; +pub mod persistent_homology; + + +pub use bottleneck::bottleneck_distance; +pub use persistent_homology::{ + persistent_homology, vietoris_rips_filtration, PersistenceDiagram, PersistencePair, +}; diff --git a/src/topology/persistent_homology.rs b/src/topology/persistent_homology.rs new file mode 100644 index 0000000..968c481 --- /dev/null +++ b/src/topology/persistent_homology.rs @@ -0,0 +1,353 @@ +//! Vietoris--Rips persistent homology via standard matrix reduction. +//! +//! Given a finite point cloud `X = {x_1, ..., x_n} \subset R^d` and a +//! filtration parameter `epsilon >= 0`, the Vietoris--Rips complex +//! `VR_eps(X)` contains every simplex `sigma = {x_{i_0}, ..., x_{i_k}}` +//! whose pairwise distances satisfy `d(x_{i_p}, x_{i_q}) <= eps`. +//! +//! The resulting filtration `{ VR_eps(X) }_{eps >= 0}` defines a +//! persistence module whose decomposition yields the diagram +//! `D_k(X) = { (b_i, d_i) }` of birth/death pairs in homological +//! degree `k`. +//! +//! This implementation: +//! * builds the simplicial complex up to a chosen homological degree +//! `max_dim` and a maximum scale `max_eps`; +//! * orders simplices by `(filtration_value, dimension, lexicographic)`; +//! * runs left-to-right column reduction over `Z/2`; +//! * extracts pairs `(b, d)` from low entries. +//! +//! The implementation is intentionally simple and CPU-friendly. For very +//! large complexes, switch to a sparse / chunk-based variant. + +use crate::core::{OptimizrError, Result}; + +/// Filtration entry. +#[derive(Debug, Clone)] +pub struct Simplex { + pub vertices: Vec, + pub filtration: f64, +} + +impl Simplex { + pub fn dim(&self) -> usize { + self.vertices.len() - 1 + } +} + +/// One persistent homology pair. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct PersistencePair { + pub dim: usize, + pub birth: f64, + pub death: f64, // f64::INFINITY for essential classes +} + +/// Persistence diagram grouped by homological degree. +#[derive(Debug, Clone, Default)] +pub struct PersistenceDiagram { + pub pairs: Vec, +} + +impl PersistenceDiagram { + pub fn pairs_in_dim(&self, dim: usize) -> Vec<(f64, f64)> { + self.pairs + .iter() + .filter(|p| p.dim == dim) + .map(|p| (p.birth, p.death)) + .collect() + } + + /// Number of generators alive at scale `eps` in dimension `dim`. + pub fn betti(&self, dim: usize, eps: f64) -> usize { + self.pairs + .iter() + .filter(|p| p.dim == dim && p.birth <= eps && eps < p.death) + .count() + } +} + +#[inline] +fn euclidean_distance(a: &[f64], b: &[f64]) -> f64 { + a.iter().zip(b.iter()).map(|(x, y)| (x - y).powi(2)).sum::().sqrt() +} + +/// Build the Vietoris--Rips filtration up to `max_dim` and scale +/// `max_eps`. Each simplex's filtration value equals the maximum pairwise +/// distance among its vertices (zero for vertices, edge length for edges). +pub fn vietoris_rips_filtration( + points: &[Vec], + max_dim: usize, + max_eps: f64, +) -> Result> { + let n = points.len(); + if n == 0 { + return Err(OptimizrError::EmptyData); + } + let d = points[0].len(); + for p in points { + if p.len() != d { + return Err(OptimizrError::DimensionMismatch { + expected: d, + actual: p.len(), + }); + } + } + if max_eps < 0.0 { + return Err(OptimizrError::InvalidParameter( + "max_eps must be non-negative".into(), + )); + } + + let mut dist = vec![0.0; n * n]; + for i in 0..n { + for j in (i + 1)..n { + let dij = euclidean_distance(&points[i], &points[j]); + dist[i * n + j] = dij; + dist[j * n + i] = dij; + } + } + + let mut simplices: Vec = Vec::new(); + + // 0-simplices + for i in 0..n { + simplices.push(Simplex { + vertices: vec![i], + filtration: 0.0, + }); + } + + // Higher-dim simplices: enumerate via incremental expansion using + // edges as the seed. For tractability we only enumerate simplices + // whose filtration <= max_eps. + if max_dim >= 1 { + for i in 0..n { + for j in (i + 1)..n { + let dij = dist[i * n + j]; + if dij <= max_eps { + simplices.push(Simplex { + vertices: vec![i, j], + filtration: dij, + }); + } + } + } + } + // Generic k-simplex enumeration by k-clique listing. + let mut current_dim = 1usize; + while current_dim < max_dim { + // Collect simplices of dimension current_dim + let prev: Vec = simplices + .iter() + .filter(|s| s.dim() == current_dim) + .cloned() + .collect(); + let mut new_simplices: Vec = Vec::new(); + for s in prev.iter() { + // Try extend by every vertex > max(vertices) such that all + // pairwise distances stay below max_eps. + let last = *s.vertices.last().unwrap(); + 'cand: for v in (last + 1)..n { + let mut max_d = s.filtration; + for &u in s.vertices.iter() { + let dd = dist[u * n + v]; + if dd > max_eps { + continue 'cand; + } + if dd > max_d { + max_d = dd; + } + } + let mut verts = s.vertices.clone(); + verts.push(v); + new_simplices.push(Simplex { + vertices: verts, + filtration: max_d, + }); + } + } + if new_simplices.is_empty() { + break; + } + simplices.extend(new_simplices); + current_dim += 1; + } + + // Order: filtration ascending, then dimension ascending, then + // lexicographic on vertices for determinism. + simplices.sort_by(|a, b| { + a.filtration + .partial_cmp(&b.filtration) + .unwrap_or(std::cmp::Ordering::Equal) + .then(a.dim().cmp(&b.dim())) + .then(a.vertices.cmp(&b.vertices)) + }); + Ok(simplices) +} + +/// Map vertex tuple -> index in the ordered filtration. +fn build_index_map(simplices: &[Simplex]) -> std::collections::HashMap, usize> { + let mut map = std::collections::HashMap::with_capacity(simplices.len()); + for (idx, s) in simplices.iter().enumerate() { + map.insert(s.vertices.clone(), idx); + } + map +} + +/// Boundary indices of a simplex (vertices removed one at a time). +fn boundary_indices( + simplex: &Simplex, + index_map: &std::collections::HashMap, usize>, +) -> Vec { + let k = simplex.vertices.len(); + if k <= 1 { + return Vec::new(); + } + let mut out = Vec::with_capacity(k); + for i in 0..k { + let mut face = Vec::with_capacity(k - 1); + for (j, &v) in simplex.vertices.iter().enumerate() { + if j != i { + face.push(v); + } + } + if let Some(&idx) = index_map.get(&face) { + out.push(idx); + } + } + out +} + +/// Run the standard Z/2 column-reduction persistence algorithm and +/// return the persistence diagram. +pub fn persistent_homology( + points: &[Vec], + max_dim: usize, + max_eps: f64, +) -> Result { + let simplices = vietoris_rips_filtration(points, max_dim, max_eps)?; + let n = simplices.len(); + let index_map = build_index_map(&simplices); + + // Sparse boundary columns sorted descending so that low(col) = first. + let mut columns: Vec> = simplices + .iter() + .map(|s| { + let mut b = boundary_indices(s, &index_map); + b.sort_unstable(); + b.dedup(); + b + }) + .collect(); + + // low_to_col[r] = column index whose low entry is r, if any. + let mut low_to_col: std::collections::HashMap = std::collections::HashMap::new(); + let mut paired: Vec = vec![false; n]; + let mut diag = PersistenceDiagram::default(); + + // Symmetric difference helper (Z/2 column add). + fn xor_sorted(a: &[usize], b: &[usize]) -> Vec { + let mut out = Vec::with_capacity(a.len() + b.len()); + let (mut i, mut j) = (0usize, 0usize); + while i < a.len() && j < b.len() { + match a[i].cmp(&b[j]) { + std::cmp::Ordering::Less => { + out.push(a[i]); + i += 1; + } + std::cmp::Ordering::Greater => { + out.push(b[j]); + j += 1; + } + std::cmp::Ordering::Equal => { + i += 1; + j += 1; + } + } + } + out.extend_from_slice(&a[i..]); + out.extend_from_slice(&b[j..]); + out + } + + for j in 0..n { + loop { + let low = match columns[j].last() { + Some(&l) => l, + None => break, + }; + if let Some(&j_prev) = low_to_col.get(&low) { + let merged = xor_sorted(&columns[j], &columns[j_prev]); + columns[j] = merged; + } else { + low_to_col.insert(low, j); + break; + } + } + if let Some(&birth_idx) = columns[j].last() { + let birth = simplices[birth_idx].filtration; + let death = simplices[j].filtration; + if death > birth { + diag.pairs.push(PersistencePair { + dim: simplices[birth_idx].dim(), + birth, + death, + }); + } + paired[birth_idx] = true; + paired[j] = true; + } + } + + // Essential classes: simplices that were neither paired as birth + // (column reduces to nonzero with low=birth_idx) nor as death. + for (i, s) in simplices.iter().enumerate() { + if !paired[i] { + // Skip top-dim or simplices > max_dim + if s.dim() <= max_dim { + diag.pairs.push(PersistencePair { + dim: s.dim(), + birth: s.filtration, + death: f64::INFINITY, + }); + } + } + } + + Ok(diag) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn unit_circle_has_one_loop() { + // 30 points uniformly on unit circle. + let n = 30; + let pts: Vec> = (0..n) + .map(|k| { + let theta = 2.0 * std::f64::consts::PI * (k as f64) / (n as f64); + vec![theta.cos(), theta.sin()] + }) + .collect(); + let diag = persistent_homology(&pts, 2, 1.5).unwrap(); + // Check that betti_1 equals 1 at an intermediate scale. + let b1 = diag.betti(1, 0.5); + assert!( + b1 >= 1, + "expected at least one 1-cycle near eps=0.5, got betti_1={}", + b1 + ); + } + + #[test] + fn two_points_zero_dim_birth() { + let pts = vec![vec![0.0, 0.0], vec![1.0, 0.0]]; + let diag = persistent_homology(&pts, 1, 2.0).unwrap(); + // Two 0-simplices born at 0; one dies at distance 1. + let dim0: Vec<_> = diag.pairs_in_dim(0); + assert!(dim0.iter().any(|(b, d)| *b == 0.0 && (*d - 1.0).abs() < 1e-12)); + } +} diff --git a/src/volterra/fourier_inversion.rs b/src/volterra/fourier_inversion.rs new file mode 100644 index 0000000..e21fccf --- /dev/null +++ b/src/volterra/fourier_inversion.rs @@ -0,0 +1,110 @@ +//! Recover a probability density from its characteristic function via a +//! direct Fourier inversion on a uniform frequency grid. +//! +//! For a real-valued random variable with characteristic function +//! `phi(u) = E[ exp(i u X) ]`, the density (when it exists) is +//! +//! ```text +//! f(x) = (1 / (2 pi)) * int_{-infty}^{+infty} exp(-i u x) phi(u) du. +//! ``` +//! +//! On a truncated symmetric grid `u in [-U, U]` with spacing `du`, the +//! discrete approximation is +//! +//! ```text +//! f_hat(x) ~= (du / (2 pi)) * sum_k exp(-i u_k x) phi(u_k). +//! ``` +//! +//! Symmetry of real densities (`phi(-u) = conj(phi(u))`) reduces the +//! computation to a cosine sum over `u >= 0`. +//! +//! # Notes +//! +//! * Implementation uses a direct `O(N_u * N_x)` evaluation. For very +//! large grids, swap in an FFT-based Carr--Madan style routine. +//! * `phi` is supplied as a closure returning `(re, im)`. + +use crate::core::{OptimizrError, Result}; + +/// Inversion result. +#[derive(Debug, Clone)] +pub struct DensityResult { + pub x_grid: Vec, + pub density: Vec, +} + +/// Recover the density on `x_grid` from the characteristic function `phi`. +/// +/// `u_max` is the truncation bound; `n_u` the number of nodes on +/// `[0, u_max]`. The total grid is `2*n_u - 1` symmetric nodes. +pub fn fourier_invert

( + phi: P, + x_grid: &[f64], + u_max: f64, + n_u: usize, +) -> Result +where + P: Fn(f64) -> (f64, f64), +{ + if x_grid.is_empty() { + return Err(OptimizrError::EmptyData); + } + if u_max <= 0.0 { + return Err(OptimizrError::InvalidParameter("u_max > 0 required".into())); + } + if n_u < 2 { + return Err(OptimizrError::InvalidParameter("n_u >= 2 required".into())); + } + let du = u_max / (n_u - 1) as f64; + // Use the symmetric form: f(x) = (1 / pi) * int_0^inf [ Re(phi(u)) cos(u x) + // + Im(phi(u)) sin(u x) ] du + // (when X is not necessarily symmetric, the imaginary part contributes + // an antisymmetric kernel). + let mut density = vec![0.0; x_grid.len()]; + // Trapezoidal weights + let weight = |k: usize| -> f64 { + if k == 0 || k == n_u - 1 { + 0.5 + } else { + 1.0 + } + }; + let cache: Vec<(f64, f64)> = (0..n_u).map(|k| phi(k as f64 * du)).collect(); + for (idx, &x) in x_grid.iter().enumerate() { + let mut s = 0.0; + for k in 0..n_u { + let u = k as f64 * du; + let (re, im) = cache[k]; + let cos_ux = (u * x).cos(); + let sin_ux = (u * x).sin(); + s += weight(k) * (re * cos_ux + im * sin_ux); + } + density[idx] = du * s / std::f64::consts::PI; + } + Ok(DensityResult { + x_grid: x_grid.to_vec(), + density, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Standard normal: phi(u) = exp(-u^2 / 2), density = (1 / sqrt(2 pi)) exp(-x^2 / 2). + #[test] + fn recovers_standard_normal_density() { + let phi = |u: f64| ((-0.5 * u * u).exp(), 0.0); + let x: Vec = (-50..=50).map(|k| 0.05 * k as f64).collect(); + let res = fourier_invert(phi, &x, 25.0, 2000).unwrap(); + let exact = |xx: f64| (1.0 / (2.0 * std::f64::consts::PI).sqrt()) * (-0.5 * xx * xx).exp(); + let mut max_err = 0.0_f64; + for (xi, fi) in x.iter().zip(res.density.iter()) { + let e = (fi - exact(*xi)).abs(); + if e > max_err { + max_err = e; + } + } + assert!(max_err < 1e-3, "max err = {}", max_err); + } +} diff --git a/src/volterra/fractional_riccati.rs b/src/volterra/fractional_riccati.rs new file mode 100644 index 0000000..08f008b --- /dev/null +++ b/src/volterra/fractional_riccati.rs @@ -0,0 +1,156 @@ +//! Adams predictor--corrector for Caputo fractional ODEs. +//! +//! Solves +//! +//! ```text +//! D^alpha h(t) = F(t, h(t)), h(0) = h_0, alpha in (0, 1) +//! ``` +//! +//! by the Diethelm--Ford--Freed (2002) fractional Adams scheme. +//! +//! Predictor (Adams--Bashforth): +//! +//! ```text +//! h^P_{n+1} = h_0 + (1 / Gamma(alpha)) * sum_{k=0}^n b_{n+1,k} F(t_k, h_k) +//! ``` +//! +//! Corrector (Adams--Moulton): +//! +//! ```text +//! h_{n+1} = h_0 + (1 / Gamma(alpha+2)) * [ F(t_{n+1}, h^P_{n+1}) +//! + sum_{k=0}^n a_{n+1,k} F(t_k, h_k) ] +//! ``` +//! +//! with weights +//! +//! ```text +//! b_{n+1,k} = h^alpha / alpha * ((n+1-k)^alpha - (n-k)^alpha) +//! ``` +//! +//! and the standard `a_{n+1,k}` coefficients of Diethelm 2002 eq. (14). +//! +//! The scheme is `O(h^{1 + alpha})` accurate. +//! +//! # Reference +//! +//! Diethelm, Ford, Freed (2002), *A predictor--corrector approach for the +//! numerical solution of fractional differential equations*, Nonlinear +//! Dynamics 29. + +use crate::core::{OptimizrError, Result}; + +/// Result of a fractional ODE integration. +#[derive(Debug, Clone)] +pub struct FractionalOdeResult { + pub t_grid: Vec, + pub h: Vec, +} + +#[inline] +fn gamma(x: f64) -> f64 { + statrs::function::gamma::gamma(x) +} + +/// Solve `D^alpha h = F(t, h)` on `[0, T]` with `n_steps + 1` equally +/// spaced points and Caputo derivative of order `alpha in (0, 1)`. +pub fn solve_fractional_ode( + h0: f64, + alpha: f64, + t_horizon: f64, + n_steps: usize, + rhs: F, +) -> Result +where + F: Fn(f64, f64) -> f64, +{ + if !(0.0 < alpha && alpha < 1.0) { + return Err(OptimizrError::InvalidParameter( + "alpha must lie in (0, 1)".into(), + )); + } + if t_horizon <= 0.0 { + return Err(OptimizrError::InvalidParameter("t_horizon > 0 required".into())); + } + if n_steps == 0 { + return Err(OptimizrError::InvalidParameter("n_steps > 0 required".into())); + } + let dt = t_horizon / n_steps as f64; + let mut t = vec![0.0; n_steps + 1]; + let mut h = vec![0.0; n_steps + 1]; + let mut f_cache = vec![0.0; n_steps + 1]; + h[0] = h0; + for k in 1..=n_steps { + t[k] = k as f64 * dt; + } + f_cache[0] = rhs(t[0], h[0]); + + let g_alpha = gamma(alpha); + let g_alpha_plus_2 = gamma(alpha + 2.0); + let dt_alpha = dt.powf(alpha); + + for n in 0..n_steps { + // ===== Predictor (fractional Adams--Bashforth) ===== + let mut sum_b = 0.0; + for k in 0..=n { + let nk = (n + 1 - k) as f64; + let nkm = (n - k) as f64; + let b = nk.powf(alpha) - nkm.powf(alpha); + sum_b += b * f_cache[k]; + } + let h_pred = h0 + dt_alpha / (alpha * g_alpha) * sum_b; + + // ===== Corrector (fractional Adams--Moulton) ===== + let f_pred = rhs(t[n + 1], h_pred); + let mut sum_a = 0.0; + for k in 0..=n { + sum_a += a_weight(n, k, alpha) * f_cache[k]; + } + let h_new = h0 + dt_alpha / g_alpha_plus_2 * (f_pred + sum_a); + h[n + 1] = h_new; + f_cache[n + 1] = rhs(t[n + 1], h[n + 1]); + } + + Ok(FractionalOdeResult { t_grid: t, h }) +} + +fn a_weight(n: usize, k: usize, alpha: f64) -> f64 { + // a_{n+1, k} for 0 <= k <= n (Diethelm 2002 eq. 14): + // + // a_{n+1, 0} = n^{alpha+1} - (n - alpha)*(n+1)^alpha + // a_{n+1, k} = (n - k + 2)^{alpha+1} + (n - k)^{alpha+1} + // - 2 (n - k + 1)^{alpha+1} (1 <= k <= n) + // + // (Final h^alpha / Gamma(alpha + 2) factor handled in caller.) + let n_f = n as f64; + if k == 0 { + n_f.powf(alpha + 1.0) - (n_f - alpha) * (n_f + 1.0).powf(alpha) + } else { + let a = (n_f - k as f64 + 2.0).powf(alpha + 1.0); + let b = (n_f - k as f64).powf(alpha + 1.0); + let c = 2.0 * (n_f - k as f64 + 1.0).powf(alpha + 1.0); + a + b - c + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// For alpha -> 1, the Caputo equation reduces to the classical ODE + /// h' = -h, h(0) = 1, with solution h(t) = exp(-t). + #[test] + fn limit_alpha_one_matches_exponential_decay() { + let alpha = 0.999; + let res = solve_fractional_ode(1.0, alpha, 2.0, 4000, |_t, h| -h).unwrap(); + let final_t = res.t_grid[4000]; + let analytic = (-final_t).exp(); + let err = (res.h[4000] - analytic).abs(); + assert!(err < 5e-2, "err = {}", err); + } + + #[test] + fn rejects_invalid_alpha() { + assert!(solve_fractional_ode(1.0, 0.0, 1.0, 10, |_, _| 0.0).is_err()); + assert!(solve_fractional_ode(1.0, 1.5, 1.0, 10, |_, _| 0.0).is_err()); + } +} diff --git a/src/volterra/markovian_lift.rs b/src/volterra/markovian_lift.rs new file mode 100644 index 0000000..91c7ce2 --- /dev/null +++ b/src/volterra/markovian_lift.rs @@ -0,0 +1,158 @@ +//! Multi-exponential approximation of a convolution kernel. +//! +//! Given a target kernel `K : (0, T] -> R`, find non-negative weights +//! `c_j` and rates `gamma_j` so that +//! +//! ```text +//! K(t) ~= sum_{j=1}^N c_j * exp(-gamma_j * t), t in (0, T]. +//! ``` +//! +//! When `K` admits the integral representation +//! +//! ```text +//! K(t) = int_0^infty exp(-gamma t) nu(d gamma) +//! ``` +//! +//! a Gauss--Jacobi quadrature on `nu` yields nodes `gamma_j` and weights +//! `c_j`. For the rough Volterra kernel `K(t) = t^{H - 1/2} / Gamma(H + 1/2)` +//! with `H in (0, 1/2]`, `nu(d gamma) = gamma^{-1/2 - H} / Gamma(1/2 - H) Gamma(H + 1/2) d gamma` +//! and the approach reduces to a Gauss--Jacobi quadrature on `(0, infty)`. +//! +//! This module exposes: +//! +//! * `geometric_grid_lift` -- choose `gamma_j` on a geometric grid then +//! fit `c_j >= 0` by non-negative least squares +//! (active-set algorithm) on a chosen sample +//! of `K(t_i)`. +//! * `lift_quality` -- L2 error on the sampled points. + +use crate::core::{OptimizrError, Result}; + +/// Approximation of `K(t)` on `(0, T]` by `sum c_j exp(-gamma_j t)`. +#[derive(Debug, Clone)] +pub struct MarkovianLift { + pub gammas: Vec, + pub weights: Vec, +} + +impl MarkovianLift { + pub fn evaluate(&self, t: f64) -> f64 { + self.gammas + .iter() + .zip(self.weights.iter()) + .map(|(&g, &c)| c * (-g * t).exp()) + .sum() + } +} + +/// Build a Markovian lift on a geometric grid of rates and fit weights by +/// non-negative least squares. +/// +/// `kernel` is the user-provided positive kernel `K(t)` evaluated on +/// `t_samples`. `n_factors` is the number of exponentials. Rates are +/// chosen as `gamma_j = gamma_min * (gamma_max / gamma_min)^{j/(N-1)}`. +pub fn geometric_grid_lift f64>( + kernel: K, + t_samples: &[f64], + n_factors: usize, + gamma_min: f64, + gamma_max: f64, + nnls_iter: usize, +) -> Result { + if n_factors < 2 { + return Err(OptimizrError::InvalidParameter( + "n_factors must be >= 2".into(), + )); + } + if !(gamma_min > 0.0 && gamma_max > gamma_min) { + return Err(OptimizrError::InvalidParameter( + "require 0 < gamma_min < gamma_max".into(), + )); + } + if t_samples.is_empty() { + return Err(OptimizrError::EmptyData); + } + let m = t_samples.len(); + let n = n_factors; + let log_min = gamma_min.ln(); + let log_max = gamma_max.ln(); + let mut gammas = vec![0.0; n]; + for j in 0..n { + let frac = j as f64 / (n - 1) as f64; + gammas[j] = (log_min + frac * (log_max - log_min)).exp(); + } + // Build A in R^{m x n}: A[i, j] = exp(-gamma_j * t_i) + let mut a = vec![0.0; m * n]; + for i in 0..m { + for j in 0..n { + a[i * n + j] = (-gammas[j] * t_samples[i]).exp(); + } + } + let b: Vec = t_samples.iter().map(|&t| kernel(t)).collect(); + let weights = nnls(&a, &b, m, n, nnls_iter); + Ok(MarkovianLift { gammas, weights }) +} + +/// `L^2` error of the lift on the sample grid. +pub fn lift_quality f64>( + lift: &MarkovianLift, + kernel: K, + t_samples: &[f64], +) -> f64 { + let mut s = 0.0; + for &t in t_samples { + let diff = lift.evaluate(t) - kernel(t); + s += diff * diff; + } + (s / t_samples.len() as f64).sqrt() +} + +/// Lawson--Hanson NNLS-like projected gradient solver. Simple but +/// adequate for small problems (`n <= 20`). Solves +/// `min || A x - b ||^2 s.t. x >= 0`. +fn nnls(a: &[f64], b: &[f64], m: usize, n: usize, max_iter: usize) -> Vec { + // Projected gradient with backtracking on Lipschitz constant. + let mut x = vec![0.0; n]; + // Estimate Lipschitz: spectral radius of A^T A bounded by ||A||_F^2. + let frob: f64 = a.iter().map(|v| v * v).sum(); + let l = frob.max(1.0); + let lr = 1.0 / l; + for _ in 0..max_iter { + // gradient = A^T (A x - b) + let mut ax = vec![0.0; m]; + for i in 0..m { + let mut s = 0.0; + for j in 0..n { + s += a[i * n + j] * x[j]; + } + ax[i] = s - b[i]; + } + let mut grad = vec![0.0; n]; + for j in 0..n { + let mut s = 0.0; + for i in 0..m { + s += a[i * n + j] * ax[i]; + } + grad[j] = s; + } + for j in 0..n { + x[j] = (x[j] - lr * grad[j]).max(0.0); + } + } + x +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn fits_single_exponential_exactly() { + // K(t) = 0.7 exp(-2 t) + let t: Vec = (1..=200).map(|k| 0.05 * k as f64).collect(); + let lift = geometric_grid_lift(|tt| 0.7 * (-2.0 * tt).exp(), &t, 5, 0.5, 5.0, 5_000) + .unwrap(); + let err = lift_quality(&lift, |tt| 0.7 * (-2.0 * tt).exp(), &t); + assert!(err < 5e-3, "L2 err = {}", err); + } +} diff --git a/src/volterra/mod.rs b/src/volterra/mod.rs new file mode 100644 index 0000000..6bd536d --- /dev/null +++ b/src/volterra/mod.rs @@ -0,0 +1,16 @@ +//! Volterra integral / fractional ODE solvers. +//! +//! Submodules: +//! +//! * `fractional_riccati` -- Adams predictor--corrector solver for +//! Caputo fractional ODEs (Diethelm--Ford--Freed 2002). +//! * `markovian_lift` -- multi-exponential approximation of a +//! convolution kernel `K(t) ~= sum c_j exp(-gamma_j t)`. +//! * `volterra_solver` -- generic second-kind Volterra equations. +//! * `fourier_inversion` -- characteristic function -> density via FFT. + +pub mod fourier_inversion; +pub mod fractional_riccati; +pub mod markovian_lift; +pub mod volterra_solver; + diff --git a/src/volterra/volterra_solver.rs b/src/volterra/volterra_solver.rs new file mode 100644 index 0000000..eea1545 --- /dev/null +++ b/src/volterra/volterra_solver.rs @@ -0,0 +1,100 @@ +//! Generic second-kind Volterra integral equation solver. +//! +//! Solves +//! +//! ```text +//! y(t) = g(t) + int_0^t K(t - s, y(s)) ds, t in [0, T], +//! ``` +//! +//! by a product-trapezoidal quadrature on a uniform grid: +//! +//! ```text +//! y_n = g_n + dt * [ K(t_n, y_0) / 2 +//! + sum_{k=1}^{n-1} K(t_n - t_k, y_k) +//! + K(0, y_n) / 2 ]. +//! ``` +//! +//! The implicit equation in `y_n` is solved by fixed-point iteration +//! (good behaviour for Lipschitz `K(., .)` with small `dt`). + +use crate::core::{OptimizrError, Result}; + +/// Result of a Volterra integration on a uniform grid. +#[derive(Debug, Clone)] +pub struct VolterraResult { + pub t_grid: Vec, + pub y: Vec, +} + +/// Solve a scalar second-kind Volterra equation by trapezoidal product +/// integration on `n_steps + 1` uniformly spaced nodes. +pub fn solve_volterra( + g: G, + kernel: K, + t_horizon: f64, + n_steps: usize, + fixed_point_iter: usize, + fixed_point_tol: f64, +) -> Result +where + G: Fn(f64) -> f64, + K: Fn(f64, f64) -> f64, +{ + if t_horizon <= 0.0 { + return Err(OptimizrError::InvalidParameter("t_horizon > 0".into())); + } + if n_steps == 0 { + return Err(OptimizrError::InvalidParameter("n_steps > 0".into())); + } + let dt = t_horizon / n_steps as f64; + let mut t = vec![0.0; n_steps + 1]; + let mut y = vec![0.0; n_steps + 1]; + for k in 0..=n_steps { + t[k] = k as f64 * dt; + } + y[0] = g(0.0); + + for n in 1..=n_steps { + let g_n = g(t[n]); + let mut explicit_part = 0.5 * kernel(t[n], y[0]); + for k in 1..n { + explicit_part += kernel(t[n] - t[k], y[k]); + } + explicit_part *= dt; + + // Implicit fixed-point on y_n + let mut yn = y[n - 1]; // initial guess + for _ in 0..fixed_point_iter { + let yn_new = g_n + explicit_part + 0.5 * dt * kernel(0.0, yn); + if (yn_new - yn).abs() < fixed_point_tol { + yn = yn_new; + break; + } + yn = yn_new; + } + y[n] = yn; + } + Ok(VolterraResult { t_grid: t, y }) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Exact solution: y' = 1, y(0) = 0 corresponds to + /// y(t) = t = 0 + int_0^t 1 ds, with K(s, y) = 1, g(t) = 0. + #[test] + fn constant_kernel_recovers_linear_growth() { + let res = solve_volterra(|_t| 0.0, |_dt, _y| 1.0, 1.0, 1000, 50, 1e-12).unwrap(); + for k in 0..=1000 { + let analytic = res.t_grid[k]; + assert!( + (res.y[k] - analytic).abs() < 1e-10, + "k={} y={} t={}", + k, + res.y[k], + analytic + ); + } + } +}