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(all(test, feature = "gpu-tests"))]
79mod gpu_tests;
80
81#[cfg(test)]
82mod wave_aaa58_e2e {
83    use super::*;
84
85    #[test]
86    fn e2e_lista_sparse_coding() {
87        let n = 8_usize;
88        let m = 12_usize;
89        let mut dict = vec![0.0_f64; m * n];
90        for i in 0..m {
91            for j in 0..n {
92                dict[i * n + j] = if (i + j) % 3 == 0 { 0.5 } else { -0.1 };
93            }
94        }
95        let cfg = ListaConfig {
96            n_measurements: n,
97            n_atoms: m,
98            threshold: 0.1,
99            n_layers: 10,
100        };
101        let lista = Lista::from_dict(&dict, cfg).expect("ok");
102        let y = vec![0.5_f64; n];
103        let code = lista.encode(&y).expect("ok");
104        assert_eq!(code.len(), m);
105        assert!(code.iter().all(|v| v.is_finite()), "all finite");
106        let nnz = code.iter().filter(|&&v| v.abs() > 1e-6).count();
107        assert!(nnz <= m, "nnz {nnz} <= {m}");
108    }
109
110    #[test]
111    fn e2e_rpca_gd_low_rank_recovery() {
112        let m = 10_usize;
113        let n = 8_usize;
114        let u_true: Vec<f64> = (0..m).map(|i| (i + 1) as f64 * 0.1).collect();
115        let v_true: Vec<f64> = (0..n).map(|j| (j + 1) as f64 * 0.1).collect();
116        let mut mat = vec![0.0_f64; m * n];
117        for i in 0..m {
118            for j in 0..n {
119                mat[i * n + j] = u_true[i] * v_true[j];
120            }
121        }
122        mat[0] += 1.0;
123
124        let cfg = RpcaGdConfig {
125            rank: 1,
126            sparsity_fraction: 0.05,
127            max_iter: 200,
128            lr: 0.01,
129            tol: 1e-4,
130        };
131        let mut rng = LcgRng::new(58);
132        let mut rpca = RpcaGd::new(m, n, cfg, &mut rng).expect("ok");
133        rpca.fit(&mat).expect("ok");
134        let l_hat = rpca.low_rank();
135        assert_eq!(l_hat.len(), m * n);
136        assert!(l_hat.iter().all(|v| v.is_finite()), "L finite");
137        let resid = rpca.residual_norm(&mat);
138        assert!(resid.is_finite(), "residual finite");
139    }
140}