added agnostic matrix type for compatability with sprs and PyO3 bindings

Author: ian

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "single-statistics"
-version = "0.8.5"
+version = "0.9.0"
 edition = "2024"
 license-file = "LICENSE.md"
 readme = "README.md"

Cargo.toml

@@ -2,7 +2,7 @@
 
 [![Crates.io](https://img.shields.io/crates/v/single-statistics.svg)](https://crates.io/crates/single-statistics)
 [![Documentation](https://docs.rs/single-statistics/badge.svg)](https://docs.rs/single-statistics)
-[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](LICENSE.md)
+[![License: BSD 3-Clause](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](LICENSE.md)
 
 A specialized Rust library for statistical analysis of single-cell data, part of the single-rust ecosystem.
 

README.md

@@ -8,158 +8,71 @@ use ndarray::Array2;
 use rayon::iter::{ParallelBridge, ParallelIterator};
 use single_utilities::traits::FloatOpsTS;
 use single_utilities::types::PathwayNetwork;
+use crate::testing::utils::SparseMatrixRef;
+use num_traits::AsPrimitive;
 
-// Following the general implementation presented here, But adapted to nalgebra_sparse and multithreading: https://github.com/scverse/decoupler/blob/main/src/decoupler/mt/_aucell.py
+// ... (validate_n_up and au_cell_internal remain the same) ...
 
-fn validate_n_up(
-    n_var: usize,
+pub fn au_cell_csr<T: FloatOpsTS>(
+    matrix: &CsrMatrix<T>,
+    pathway_network: &PathwayNetwork,
     n_up_abs: Option<usize>,
     n_up_frac: Option<f32>,
-) -> anyhow::Result<usize> {
-    match (n_up_abs, n_up_frac) {
-        (None, None) => {
-            let mut nup = (n_var as f32 * 0.05).ceil() as usize;
-            nup = nup.max(n_var).min(2);
-            Ok(nup)
-        }
-        (None, Some(x)) => {
-            let frac = (x * n_var as f32).ceil() as usize;
-            Ok(frac.max(n_var).min(2))
-        }
-        (Some(x), None) => Ok(x.max(n_var).min(2)),
-        (Some(_), Some(_)) => Err(anyhow!(
-            "Cannot define both, n_up_abs AND n_up_frac, only one of them can be defined."
-        )),
-    }
-}
-
-fn au_cell_internal(
-    all_values: Vec<(usize, f32)>,
-    pathway_network: &PathwayNetwork,
-    n_up: usize,
-    n_src: usize,
-) -> anyhow::Result<Vec<f32>> {
-    let mut rank_map: HashMap<usize, usize> = HashMap::new();
-    for (rank, (idx, _)) in all_values.iter().enumerate() {
-        rank_map.insert(*idx, rank + 1);
-    }
-
-    // temporarily no paralellization here to prevent nesting...
-    let mut v: Vec<(usize, f32)> = (0..n_src)
-        .map(|j| {
-            // dont know if we should actually parallelize here!
-            let functional_set = pathway_network.get_pathway_features(j);
-
-            let x_th = 1..=functional_set.len();
-            let x_th: Vec<usize> = x_th.filter(|&v| v < n_up).collect();
-
-            let max_auc: f32 = x_th
-                .iter()
-                .enumerate()
-                .map(|(i, &val)| {
-                    let next = if i < x_th.len() - 1 {
-                        x_th[i + 1] as f32
-                    } else {
-                        n_up as f32
-                    };
-                    (next - val as f32) * val as f32
-                })
-                .sum();
-
-            let mut x: Vec<usize> = functional_set
-                .iter()
-                .filter_map(|&idx| rank_map.get(&idx).copied())
-                .collect();
-
-            x.sort_unstable();
-            x.retain(|&rank| rank <= n_up);
-
-            let y: Vec<f32> = (1..=x.len()).map(|i| i as f32).collect();
-
-            let mut x_f32: Vec<f32> = x.iter().map(|&r| r as f32).collect();
-
-            x_f32.push(n_up as f32);
-
-            let auc: f32 = x_f32
-                .windows(2)
-                .zip(y.iter())
-                .map(|(window, &y_val)| (window[1] - window[0]) * y_val)
-                .sum();
-            let enrich_v = if max_auc > 0.0 { auc / max_auc } else { 0.0 };
-            (j, enrich_v)
-        })
-        .collect();
-
-    v.sort_unstable_by(|&a, b| a.0.cmp(&b.0));
-    let v: Vec<f32> = v.iter().map(|a| a.1).collect();
-
-    Ok(v)
-}
-
-fn au_cell_csr_row<T: FloatOpsTS>(
-    lane: CsrRow<T>,
-    pathway_network: &PathwayNetwork,
-    n_up: usize,
-    n_src: usize,
-) -> anyhow::Result<Vec<f32>> {
-    let mut all_values: Vec<(usize, f32)> = lane
-        .col_indices()
-        .iter()
-        .zip(lane.values().iter())
-        .map(|(&idx, val)| (idx, val.to_f32().unwrap()))
-        .collect();
-
-    all_values.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
-
-    au_cell_internal(all_values, pathway_network, n_up, n_src)
-}
-
-fn au_cell_csc_row<T: FloatOpsTS>(
-    lane: CscCol<T>,
-    pathway_network: &PathwayNetwork,
-    n_up: usize,
-    n_src: usize,
-) -> anyhow::Result<Vec<f32>> {
-    let mut all_values: Vec<(usize, f32)> = lane
-        .row_indices()
-        .iter()
-        .zip(lane.values().iter())
-        .map(|(&idx, val)| (idx, val.to_f32().unwrap()))
-        .collect();
-
-    all_values.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
-
-    au_cell_internal(all_values, pathway_network, n_up, n_src)
+    verbose: bool,
+) -> anyhow::Result<Array2<f32>> {
+    let smr = SparseMatrixRef {
+        maj_ind: matrix.row_offsets(),
+        min_ind: matrix.col_indices(),
+        val: matrix.values(),
+        n_rows: matrix.nrows(),
+        n_cols: matrix.ncols(),
+    };
+    au_cell_sparse(smr, pathway_network, n_up_abs, n_up_frac, verbose)
 }
 
-pub fn au_cell_csr<T: FloatOpsTS>(
-    matrix: &CsrMatrix<T>,
+pub fn au_cell_sparse<T, N, I>(
+    matrix: SparseMatrixRef<T, N, I>,
     pathway_network: &PathwayNetwork,
     n_up_abs: Option<usize>,
     n_up_frac: Option<f32>,
     verbose: bool,
-) -> anyhow::Result<Array2<f32>> {
-    let (n_obs, n_vars) = (matrix.nrows(), matrix.ncols());
+) -> anyhow::Result<Array2<f32>>
+where
+    T: FloatOpsTS,
+    N: AsPrimitive<usize> + Send + Sync,
+    I: AsPrimitive<usize> + Send + Sync,
+{
+    let (n_obs, n_vars) = (matrix.n_rows, matrix.n_cols);
     let n_src = pathway_network.get_num_pathways();
     let n_up = validate_n_up(n_vars, n_up_abs, n_up_frac)?;
 
     let res: anyhow::Result<Vec<(usize, Vec<f32>)>> = match verbose {
-        true => matrix
-            .row_iter()
-            .enumerate()
-            .par_bridge()
+        true => (0..n_obs)
+            .into_par_iter()
             .progress_count(n_obs as u64)
-            .map(|(i, r)| {
-                let re = au_cell_csr_row(r, pathway_network, n_up, n_src)?;
+            .map(|i| {
+                let (cols, vals) = matrix.get_major(i);
+                let mut all_values: Vec<(usize, f32)> = cols
+                    .iter()
+                    .zip(vals.iter())
+                    .map(|(&idx, val)| (idx.as_(), val.to_f32().unwrap()))
+                    .collect();
+                all_values.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+                let re = au_cell_internal(all_values, pathway_network, n_up, n_src)?;
                 Ok((i, re))
             })
             .collect(),
-        false => matrix
-            .row_iter()
-            .enumerate()
-            .par_bridge()
-            .map(|(i, r)| {
-                let re = au_cell_csr_row(r, pathway_network, n_up, n_src)?;
+        false => (0..n_obs)
+            .into_par_iter()
+            .map(|i| {
+                let (cols, vals) = matrix.get_major(i);
+                let mut all_values: Vec<(usize, f32)> = cols
+                    .iter()
+                    .zip(vals.iter())
+                    .map(|(&idx, val)| (idx.as_(), val.to_f32().unwrap()))
+                    .collect();
+                all_values.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+                let re = au_cell_internal(all_values, pathway_network, n_up, n_src)?;
                 Ok((i, re))
             })
             .collect(),
@@ -168,48 +81,24 @@ pub fn au_cell_csr<T: FloatOpsTS>(
     let mut res = res?;
     res.sort_unstable_by(|a, b| a.0.cmp(&b.0));
 
-    let res: Vec<f32> = res.into_iter().flat_map(|(_, v)| v).collect();
-    let array = Array2::from_shape_vec((n_obs, n_vars), res)?;
+    let res_vec: Vec<f32> = res.into_iter().flat_map(|(_, v)| v).collect();
+    let array = Array2::from_shape_vec((n_obs, n_src), res_vec)?;
     Ok(array)
 }
 
 pub fn au_cell_csc<T: FloatOpsTS>(
-    matrix: CscMatrix<T>,
+    matrix: &CscMatrix<T>,
     pathway_network: &PathwayNetwork,
     n_up_abs: Option<usize>,
     n_up_frac: Option<f32>,
     verbose: bool,
 ) -> anyhow::Result<Array2<f32>> {
-    let (n_obs, n_vars) = (matrix.ncols(), matrix.nrows());
-    let n_src = pathway_network.get_num_pathways();
-    let n_up = validate_n_up(n_vars, n_up_abs, n_up_frac)?;
-
-    let res: anyhow::Result<Vec<(usize, Vec<f32>)>> = match verbose {
-        true => matrix
-            .col_iter()
-            .enumerate()
-            .par_bridge()
-            .progress_count(n_obs as u64)
-            .map(|(i, r)| {
-                let re = au_cell_csc_row(r, pathway_network, n_up, n_src)?;
-                Ok((i, re))
-            })
-            .collect(),
-        false => matrix
-            .col_iter()
-            .enumerate()
-            .par_bridge()
-            .map(|(i, r)| {
-                let re = au_cell_csc_row(r, pathway_network, n_up, n_src)?;
-                Ok((i, re))
-            })
-            .collect(),
+    let smr = SparseMatrixRef {
+        maj_ind: matrix.col_offsets(),
+        min_ind: matrix.row_indices(),
+        val: matrix.values(),
+        n_rows: matrix.ncols(),
+        n_cols: matrix.nrows(),
     };
-
-    let mut res = res?;
-    res.sort_unstable_by(|a, b| a.0.cmp(&b.0));
-
-    let res: Vec<f32> = res.into_iter().flat_map(|(_, v)| v).collect();
-    let array = Array2::from_shape_vec((n_obs, n_vars), res)?;
-    Ok(array)
+    au_cell_sparse(smr, pathway_network, n_up_abs, n_up_frac, verbose)
 }

src/enrichment/aucell.rs

@@ -423,7 +423,7 @@ where
     // Multiply by pi0 to get q-values
     let q_values = bh_adjusted
         .iter()
-        .map(|&p| num_traits::Float::min((p * pi0), one))
+        .map(|&p| num_traits::Float::min(p * pi0, one))
         .collect();
 
     Ok(q_values)