Skip to main content

oxicuda_cs/
lib.rs

1//! `oxicuda-cs` — Compressed Sensing, Sparse Recovery, and Low-Rank Matrix Completion for OxiCUDA.
2//!
3//! # Architecture
4//!
5//! ```text
6//! oxicuda-cs
7//! ├── greedy/             — OMP, StOMP, ROMP, CoSaMP, Subspace Pursuit
8//! ├── thresholding/       — IHT, NIHT, HTP, Accelerated IHT
9//! ├── amp/                — AMP, VAMP, Empirical-Bayes AMP
10//! ├── basis_pursuit/      — Basis Pursuit (ADMM), BPDN, Dantzig Selector
11//! ├── lasso/              — Coord descent (Friedman et al.), LARS, FISTA-LASSO,
12//! │                         group/fused LASSO, Elastic Net
13//! ├── tv/                 — 1D/2D Chambolle Total Variation denoising
14//! ├── matrix_completion/  — SVT, Nuclear-norm minimisation, ADMM matrix completion
15//! ├── robust_pca/         — Principal Component Pursuit (PCP), GoDec
16//! ├── sparse_pca/         — Witten-Tibshirani-Hastie penalised matrix decomposition
17//! ├── sbl/                — Sparse Bayesian Learning, Fast Marginal Likelihood
18//! ├── dictionary/         — K-SVD, MOD, Online dictionary learning
19//! ├── measurement/        — Gaussian, Bernoulli, Partial Fourier matrices, RIP estimator
20//! ├── linalg/             — Private: Jacobi SVD, Householder QR, Cholesky, LSQR, normal equations
21//! └── metrics/            — sparsity, recovery error, support recovery rate, MSE, PSNR, SNR
22//! ```
23//!
24//! All algorithms are implemented in pure Rust with no external linear-algebra dependencies.
25//! Random sampling uses the workspace `LcgRng` (MMIX LCG with bit-32 boolean trick).
26
27#![forbid(unsafe_code)]
28#![allow(clippy::needless_borrows_for_generic_args)]
29#![allow(clippy::useless_vec)]
30#![allow(clippy::needless_range_loop)]
31#![allow(clippy::manual_div_ceil)]
32#![allow(clippy::manual_range_contains)]
33
34pub mod amp;
35pub mod basis_pursuit;
36pub mod dictionary;
37pub mod error;
38pub mod greedy;
39pub mod handle;
40pub mod lasso;
41pub mod linalg;
42pub mod matrix_completion;
43pub mod measurement;
44pub mod metrics;
45pub mod ptx_advanced;
46pub mod ptx_kernels;
47pub mod robust_pca;
48pub mod sbl;
49pub mod sparse_pca;
50pub mod thresholding;
51pub mod tv;
52
53pub use error::{CsError, CsResult};
54pub use handle::{CsHandle, LcgRng, SmVersion};
55
56// Wave AAA+58 re-exports.
57pub use greedy::{Lista, ListaConfig};
58pub use lasso::{Slope, SlopeConfig, sorted_l1_prox};
59pub use robust_pca::{RpcaGd, RpcaGdConfig};
60
61// Wave AAA+68 re-exports.
62pub use dictionary::{CoupledDictionary, CoupledDlConfig, couple_code, coupled_dl};
63pub use measurement::{CodedDiffraction, MaskKind, WirtingerConfig, phase_aligned_error};
64pub use sbl::{Rvm, RvmConfig, RvmFit, RvmKernel, rvm_fit_design};
65
66// Sequential compressed sensing via SMC (Ji, Xue & Carin 2008).
67pub use sbl::{SmcCs, SmcCsConfig};
68
69// Architecture-specialised PTX kernel variants + per-SM tile configuration.
70pub use ptx_advanced::{
71    TileConfig, correlate_fp8_ptx, correlate_tma_ptx, iht_step_cp_async_ptx,
72    svt_threshold_warpshuffle_ptx,
73};
74
75#[cfg(test)]
76mod e2e_tests;
77
78#[cfg(test)]
79mod wave_aaa58_e2e {
80    use super::*;
81
82    #[test]
83    fn e2e_lista_sparse_coding() {
84        let n = 8_usize;
85        let m = 12_usize;
86        let mut dict = vec![0.0_f64; m * n];
87        for i in 0..m {
88            for j in 0..n {
89                dict[i * n + j] = if (i + j) % 3 == 0 { 0.5 } else { -0.1 };
90            }
91        }
92        let cfg = ListaConfig {
93            n_measurements: n,
94            n_atoms: m,
95            threshold: 0.1,
96            n_layers: 10,
97        };
98        let lista = Lista::from_dict(&dict, cfg).expect("ok");
99        let y = vec![0.5_f64; n];
100        let code = lista.encode(&y).expect("ok");
101        assert_eq!(code.len(), m);
102        assert!(code.iter().all(|v| v.is_finite()), "all finite");
103        let nnz = code.iter().filter(|&&v| v.abs() > 1e-6).count();
104        assert!(nnz <= m, "nnz {nnz} <= {m}");
105    }
106
107    #[test]
108    fn e2e_rpca_gd_low_rank_recovery() {
109        let m = 10_usize;
110        let n = 8_usize;
111        let u_true: Vec<f64> = (0..m).map(|i| (i + 1) as f64 * 0.1).collect();
112        let v_true: Vec<f64> = (0..n).map(|j| (j + 1) as f64 * 0.1).collect();
113        let mut mat = vec![0.0_f64; m * n];
114        for i in 0..m {
115            for j in 0..n {
116                mat[i * n + j] = u_true[i] * v_true[j];
117            }
118        }
119        mat[0] += 1.0;
120
121        let cfg = RpcaGdConfig {
122            rank: 1,
123            sparsity_fraction: 0.05,
124            max_iter: 200,
125            lr: 0.01,
126            tol: 1e-4,
127        };
128        let mut rng = LcgRng::new(58);
129        let mut rpca = RpcaGd::new(m, n, cfg, &mut rng).expect("ok");
130        rpca.fit(&mat).expect("ok");
131        let l_hat = rpca.low_rank();
132        assert_eq!(l_hat.len(), m * n);
133        assert!(l_hat.iter().all(|v| v.is_finite()), "L finite");
134        let resid = rpca.residual_norm(&mat);
135        assert!(resid.is_finite(), "residual finite");
136    }
137}