legume-numeric 0.8.11

Numeric and ML foundation for the legume ecosystem (matrix, Leiden, candle, MCMC)
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

use super::*;
use nalgebra::DMatrix;
use rand::rngs::StdRng;
use rand::{RngExt, SeedableRng};

fn random_rows(n: usize, s: usize, seed: u64) -> DMatrix<f32> {
    let mut rng = StdRng::seed_from_u64(seed);
    DMatrix::from_fn(n, s, |_, _| rng.random_range(-1.0f32..1.0))
}

/// The shared oracle over every row: nearest first, ties by index, self
/// excluded.
fn brute(x: &DMatrix<f32>, k: usize) -> (Vec<Vec<usize>>, Vec<Vec<f32>>) {
    let points: Vec<Vec<f32>> = (0..x.nrows())
        .map(|i| x.row(i).iter().copied().collect())
        .collect();
    (0..x.nrows())
        .map(|q| crate::matrix::knn::tests::brute_others(&points, q, k))
        .unzip()
}

////////////////////////////////////////
// Exactness across block boundaries //
////////////////////////////////////////

#[test]
fn matches_brute_force_across_block_boundaries() {
    let x = random_rows(97, 6, 1);
    let k = 5;
    let (nb, ds) = knn_rows_l2_blocked(&x, k, 8);
    let (tnb, tds) = brute(&x, k);
    assert_eq!(nb.len(), 97);
    for i in 0..97 {
        assert_eq!(nb[i], tnb[i], "row {i} neighbours");
        for (a, b) in ds[i].iter().zip(&tds[i]) {
            assert!((a - b).abs() < 1e-4, "row {i}: {a} vs {b}");
        }
        for w in ds[i].windows(2) {
            assert!(w[0] <= w[1], "row {i} not nearest-first");
        }
    }
}

#[test]
fn default_block_matches_small_block() {
    let x = random_rows(150, 7, 2);
    assert_eq!(knn_rows_l2(&x, 4), knn_rows_l2_blocked(&x, 4, 5));
}

////////////////////////////
// Ties and degeneracies //
////////////////////////////

#[test]
fn ties_break_by_index_and_are_deterministic() {
    // Four copies of each of three distinct points: every copy's nearest
    // three are the other copies, at distance zero, in index order.
    let base = random_rows(3, 5, 3);
    let x = DMatrix::from_fn(12, 5, |i, c| base[(i % 3, c)]);
    let (nb, ds) = knn_rows_l2_blocked(&x, 3, 4);
    for i in 0..12 {
        let expect: Vec<usize> = (0..12).filter(|&j| j != i && j % 3 == i % 3).collect();
        assert_eq!(nb[i], expect, "row {i}");
        assert!(ds[i].iter().all(|&d| d == 0.0), "row {i}: {:?}", ds[i]);
    }
    assert_eq!(knn_rows_l2_blocked(&x, 3, 4), knn_rows_l2_blocked(&x, 3, 4));
}

#[test]
fn k_is_clamped_to_the_other_rows() {
    let x = random_rows(4, 3, 4);
    let (nb, ds) = knn_rows_l2(&x, 10);
    assert!(nb.iter().all(|v| v.len() == 3));
    assert!(ds.iter().all(|v| v.len() == 3));
    let one = random_rows(1, 3, 5);
    let (nb, ds) = knn_rows_l2(&one, 3);
    assert_eq!(nb, vec![Vec::<usize>::new()]);
    assert_eq!(ds, vec![Vec::<f32>::new()]);
    let none = DMatrix::<f32>::zeros(0, 3);
    assert!(knn_rows_l2(&none, 3).0.is_empty());
}

#[test]
fn non_finite_rows_sort_last() {
    let mut x = random_rows(20, 4, 6);
    x[(7, 2)] = f32::NAN;
    let (nb, ds) = knn_rows_l2_blocked(&x, 5, 6);
    for i in 0..20 {
        if i == 7 {
            assert!(ds[i].iter().all(|d| !d.is_finite()), "NaN row's distances");
        } else {
            assert!(!nb[i].contains(&7), "row {i} lists the NaN row");
            assert!(ds[i].iter().all(|d| d.is_finite()));
        }
    }
}

/////////////////////////////////////////////
// Precision: a large shared offset must   //
// not corrupt the neighbourhoods.         //
/////////////////////////////////////////////

#[test]
fn large_shared_offset_keeps_exact_neighbours() {
    // Gram-matrix distances cancel catastrophically when the norms dwarf the
    // separations; the returned neighbours and distances must still be the
    // direct-difference truth.
    let small = random_rows(64, 8, 7);
    let x = small.map(|v| 1.0e4 + 1.0e-2 * v);
    let (nb, ds) = knn_rows_l2_blocked(&x, 6, 16);
    let (tnb, tds) = brute(&x, 6);
    for i in 0..64 {
        assert_eq!(nb[i], tnb[i], "row {i} neighbours");
        for (a, b) in ds[i].iter().zip(&tds[i]) {
            assert!(((a - b) / b).abs() < 1e-3, "row {i}: {a} vs {b}");
        }
    }
}

////////////////////////////
// Thread-count invariance //
////////////////////////////

#[test]
fn independent_of_thread_count() {
    let x = random_rows(300, 6, 8);
    let reference = knn_rows_l2_blocked(&x, 8, 32);
    for threads in [1usize, 3] {
        let pool = rayon::ThreadPoolBuilder::new()
            .num_threads(threads)
            .build()
            .unwrap();
        let got = pool.install(|| knn_rows_l2_blocked(&x, 8, 32));
        assert_eq!(got, reference, "{threads} threads");
    }
}