massive upgrade to library, removed duplicate steps

Author: Ian

.github/workflows/publish.yml

@@ -6,7 +6,22 @@ on:
     tags:
       - "v*"
 jobs:
+  test:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+      - name: Install system dependencies
+        run: |
+          sudo apt-get update
+      - name: Run tests
+        run: cargo test --verbose
+      
   publish:
+    needs: test
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v2

.gitignore

@@ -1,3 +1,4 @@
 /target
 .idea
-.idea/
+.idea/
+.fleet/

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "single-statistics"
-version = "0.5.0"
+version = "0.6.0"
 edition = "2024"
 license-file = "LICENSE.md"
 readme = "README.md"
@@ -13,6 +13,7 @@ homepage = "https://singlerust.com"
 [dependencies]
 anyhow = "1.0.98"
 nalgebra-sparse = "0.10.0"
+ndarray = {version = "0.16.1", features = ["rayon"]}
 num-traits = "0.2.19"
 rayon = "1.10.0"
 single-utilities = "0.6.0"

Cargo.toml

@@ -0,0 +1,107 @@
+// 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
+
+use std::mem::offset_of;
+
+use nalgebra_sparse::CsrMatrix;
+use ndarray::Array2;
+use single_utilities::traits::{FloatOps, FloatOpsTS};
+
+fn validate_n_up<T>(n_genes: usize, n_up: Option<T>) -> anyhow::Result<usize> 
+where 
+T: FloatOps {
+    let n_up_value = match n_up {
+        Some(val) => {
+            let n_up_int = num_traits::Float::ceil(val).to_usize().unwrap();
+            n_up_int
+        },
+        None => {
+            let n_up_float = T::from(0.05).unwrap() * T::from(n_genes).unwrap();
+            let n_up_ceil = num_traits::Float::ceil(n_up_float);
+            let n_up_int = n_up_ceil.to_usize().unwrap();
+            n_up_int.clamp(2, n_genes)
+        },
+    };
+
+    if n_up_value <= 1 || n_up_value > n_genes {
+        return Err(anyhow::anyhow!(
+            "For n_genes={}, n_up={} must be between 1 and {}",
+            n_genes, n_up_value, n_genes
+        ));
+    }
+    Ok(n_up_value)
+}
+
+fn get_gene_set(
+    connectivity: &[usize],
+    starts: &[usize],
+    offsets: &[usize],
+    source_idx: usize
+) -> Vec<usize> {
+    let start = starts[source_idx];
+    let offset = offsets[source_idx];
+    connectivity[start..start + offset].to_vec()
+}
+
+fn rank_data<T>(values: &[T]) -> Vec<usize> 
+where
+    T: FloatOps 
+    {
+        let mut indexed_values: Vec<(usize, T)> = values
+        .iter()
+        .enumerate()
+        .map(|(i, &v)| (i,v))
+        .collect();
+
+        indexed_values.sort_by(|a, b| {
+            b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)
+        });
+
+        let mut ranks = vec![0; values.len()];
+        for (rank, (original_idx, _)) in indexed_values.iter().enumerate() {
+            ranks[*original_idx] = rank + 1;
+        }
+
+        ranks
+    }
+
+
+
+fn compute_chunk_size(n_cells: usize, n_genes: usize) -> usize {
+    let n_cores = rayon::current_num_threads();
+
+    let base_chunk_size = if n_genes > 20000 {
+        200
+    } else if n_genes > 10000 {
+        500
+    } else {
+        1000
+    };
+
+    let min_chunks = n_cores;
+    let max_chunk_size = (n_cells + min_chunks - 1) / min_chunks;
+    base_chunk_size.min(max_chunk_size).max(1)
+}
+
+pub(crate) fn aucell_compute<T>(
+    matrix: &CsrMatrix<T>,
+    connectivity: &[usize],
+    starts: &[usize],
+    offsets: &[usize],
+    n_up: Option<T>
+) -> anyhow::Result<Array2<T>>
+where 
+    T: FloatOpsTS {
+        let n_cells = matrix.nrows();
+        let n_genes = matrix.ncols();
+        let n_sources = starts.len();
+
+        if connectivity.is_empty() || starts.is_empty() || offsets.is_empty() {
+            return Err(anyhow::anyhow!("Connectivity arrays cannot be empty!!"))
+        }
+
+        if starts.len() != offsets.len() {
+            return Err(anyhow::anyhow!("Starts and offsets must have the same length!"))
+        }
+
+        todo!()
+    }

src/enrichment/aucell.rs

@@ -0,0 +1,3 @@
+mod gsea;
+mod aucell;
+mod ora;

src/enrichment/gsea.rs

@@ -1 +1,2 @@
-pub mod testing;
+pub mod testing;
+pub mod enrichment;

src/enrichment/mod.rs

@@ -42,7 +42,7 @@ where
     let n_t = T::from(n).unwrap();
     let adjusted = p_values
         .iter()
-        .map(|&p| num_traits::Float::min((p * n_t), T::one()))
+        .map(|&p| num_traits::Float::min(p * n_t, T::one()))
         .collect();
 
     Ok(adjusted)
@@ -149,7 +149,7 @@ where
 
         // Calculate adjustment and take minimum of current and previous
         let adjustment =
-            num_traits::Float::min((p_val * c_n * n_f64 / T::from(rank).unwrap()), T::one());
+            num_traits::Float::min(p_val * c_n * n_f64 / T::from(rank).unwrap(), T::one());
         current_min = num_traits::Float::min(adjustment, current_min);
         adjusted_p_values[orig_idx] = current_min;
     }
@@ -264,7 +264,7 @@ where
         let denominator_t = T::from(n - i).unwrap_or(one);
 
         let hochberg_value =
-            num_traits::Float::min((indexed_p_values[i].1 * n_t / denominator_t), one);
+            num_traits::Float::min(indexed_p_values[i].1 * n_t / denominator_t, one);
         adjusted_p_values[current_index] =
             num_traits::Float::min(hochberg_value, adjusted_p_values[prev_index]);
     }