bixverse-rs 0.4.3

A Rust library with various utilities for computational biology and single cell.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

//! Pre-processing methods for NMF on sparse and dense data.

use faer::{Mat, MatRef};
use rayon::prelude::*;

use crate::core::math::matrix_helpers::col_sds;
use crate::core::math::sparse::sparse_col_sds_csc;
use crate::prelude::*;

//////////
// Enum //
//////////

use crate::prelude::BixverseFloat;

/// Options to do NmfPreprocessing
#[derive(Clone, Copy, Debug, Default)]
pub enum NmfPreprocessing {
    /// No pre-processing.
    #[default]
    None,
    /// Scale each column (features) by the standard deviation.
    SdScaling,
    /// Scale each column (features) by the squared standard deviation.
    SqrtSdScaling,
}

/// Parse the NMF pre-processing
///
/// ### Params
///
/// * `s` - String to parse. `"none"` -> leave data as is; `"sd" | "scaling"`
///   -> divide the respective features by their standard deviation; `"sqrt_sd"`
///   -> divide the respective features by the square root of their standard
///   deviation.
///
/// ### Returns
///
/// The option of [NmfPreprocessing]
pub fn parse_nmf_processing(s: &str) -> Option<NmfPreprocessing> {
    match s.to_lowercase().as_str() {
        "none" => Some(NmfPreprocessing::None),
        "sd" | "scaling" => Some(NmfPreprocessing::SdScaling),
        "sqrt_sd" => Some(NmfPreprocessing::SqrtSdScaling),
        _ => None,
    }
}

///////////
// Dense //
///////////

/// Process dense matrices of samples x features
///
/// ### Params
///
/// * `x` - Matrix to preprocess
/// * `process` - Which approach of NMF pre-processing to apply, see
///   [NmfPreprocessing]
///
/// ### Returns
///
/// The pre-processed matrix
pub fn nmf_process_dense<T: BixverseFloat>(x: MatRef<T>, process: &NmfPreprocessing) -> Mat<T> {
    let epsilon = T::from_f32(1e-6).unwrap();

    match process {
        NmfPreprocessing::None => x.to_owned(),
        NmfPreprocessing::SdScaling => {
            let mut sds = col_sds(x);
            // safe guard against tiny values...
            sds.par_iter_mut()
                .for_each(|x| *x = if *x <= epsilon { T::one() } else { *x });

            Mat::from_fn(x.nrows(), x.ncols(), |i, j| x[(i, j)] / sds[j])
        }
        NmfPreprocessing::SqrtSdScaling => {
            let mut sds = col_sds(x);
            sds.par_iter_mut()
                .for_each(|x| *x = if *x <= epsilon { T::one() } else { x.sqrt() });

            Mat::from_fn(x.nrows(), x.ncols(), |i, j| x[(i, j)] / sds[j])
        }
    }
}

////////////
// Sparse //
////////////

/// Process sparse matrices of samples x features (need to be in CSC format!)
///
/// ### Params
///
/// * `x` - Sparse matrix to preprocess, see [CompressedSparseData2].
/// * `process` - Which approach of NMF pre-processing to apply, see
///   [NmfPreprocessing].
/// * `use_second_layer` - Use the second data layer.
///
/// ### Returns
///
/// The pre-processed matrix
pub fn nmf_process_sparse<T>(
    x: &CompressedSparseData2<T>,
    process: &NmfPreprocessing,
    use_second_layer: bool,
) -> Result<CompressedSparseData2<T>, BixverseErrors>
where
    T: BixverseFloat + BixverseSimd,
{
    let epsilon = T::from_f32(1e-6).unwrap();

    match process {
        NmfPreprocessing::None => Ok(x.clone()),
        NmfPreprocessing::SdScaling | NmfPreprocessing::SqrtSdScaling => {
            let use_sqrt = matches!(process, NmfPreprocessing::SqrtSdScaling);
            let mut sds = sparse_col_sds_csc(x, use_second_layer)?;

            sds.par_iter_mut().for_each(|s| {
                *s = if *s <= epsilon {
                    T::one()
                } else if use_sqrt {
                    s.sqrt()
                } else {
                    *s
                };
            });

            let mut out = x.clone();
            let active_data: &mut Vec<T> = if use_second_layer {
                out.data_2
                    .as_mut()
                    .ok_or(BixverseErrors::Data2NotAvailable)?
            } else {
                &mut out.data
            };

            let ncols = x.shape().1;
            for j in 0..ncols {
                let start = x.indptr[j] as usize;
                let end = x.indptr[j + 1] as usize;
                let divisor = sds[j];
                for v in &mut active_data[start..end] {
                    *v /= divisor;
                }
            }

            Ok(out)
        }
    }
}