holos-tda 0.9.0

Vietoris-Rips persistence and checked degree-Rips modules
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

use rustc_hash::FxHashMap;

use crate::field::Z2;
use crate::reduce::{Engine, Pivots};
use crate::simplex::Simplex;
use crate::{Diagram, DistanceMatrix, RipsParams};

/// Every in-complex edge, in column order: diameter descending, index
/// ascending. The reducer takes any such list, so a test can hand it one
/// without running the dim-0 pass first.
fn edge_columns(dist: &DistanceMatrix, engine: &Engine<'_, Z2, DistanceMatrix>) -> Vec<Simplex> {
    let mut columns = Vec::new();
    for i in 1..dist.len() {
        for j in 0..i {
            let diameter = dist.get(i, j);
            if engine.in_complex(diameter) {
                columns.push(Simplex {
                    diameter,
                    index: engine.bt.get(i, 2) + j as u64,
                });
            }
        }
    }
    columns.sort_unstable_by(|a, b| {
        b.diameter
            .total_cmp(&a.diameter)
            .then(a.index.cmp(&b.index))
    });
    columns
}

fn points(seed: u64, n: usize, coord_dim: usize) -> Vec<Vec<f64>> {
    let mut x = seed | 1;
    let mut next = || {
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        (x >> 11) as f64 / (1u64 << 53) as f64
    };
    (0..n)
        .map(|_| (0..coord_dim).map(|_| next()).collect())
        .collect()
}

fn grid(side: usize) -> Vec<Vec<f64>> {
    (0..side)
        .flat_map(|a| (0..side).map(move |b| vec![a as f64, b as f64]))
        .collect()
}

/// The pivot registry the workers converge to must not depend on how many
/// of them there are. The check is wider than the diagram: it compares
/// the pivot index, the coefficient, the owning column, and the diameter
/// bits, and it compares the first three against the serial reducer as
/// well.
fn assert_registry_is_worker_invariant(dist: &DistanceMatrix, label: &str) {
    let mut serial_params = RipsParams::new(1);
    serial_params.threads = 1;
    let serial = Engine::new(dist, &serial_params, Z2).unwrap();
    let columns = edge_columns(dist, &serial);
    assert!(columns.len() > 64, "{label}: too few columns to be a gate");
    let empty: Pivots = FxHashMap::default();
    let mut diagram = Diagram::default();
    let want = serial.reduce_dimension(&columns, 1, &empty, &mut diagram);

    let mut first: Option<Vec<(u64, u64, usize, u64)>> = None;
    for budget in [2usize, 3, 4, 8] {
        let mut params = RipsParams::new(1);
        params.threads = budget;
        let engine = Engine::new(dist, &params, Z2).unwrap();
        let got = engine.parallel_pivot_registry(&columns, 1, &empty, budget);

        let by_index: Pivots = got
            .iter()
            .map(|&(index, coeff, col, _)| (index, (coeff, col)))
            .collect();
        assert_eq!(
            by_index, want,
            "{label}: {budget} workers disagree with the serial pivot registry"
        );
        match &first {
            None => first = Some(got),
            Some(want) => assert_eq!(
                &got, want,
                "{label}: {budget} workers disagree with 2 workers, diameter bits included"
            ),
        }
    }
}

#[test]
fn pivot_registry_is_worker_invariant_on_a_cloud() {
    let dist = DistanceMatrix::from_points(&points(20260818, 100, 3)).unwrap();
    assert_registry_is_worker_invariant(&dist, "cloud(n=100,d=3)");
}

#[test]
fn pivot_registry_is_worker_invariant_on_ties() {
    // A lattice puts many simplices at one diameter, which is where the
    // workers reorder the most and displace each other the most.
    let dist = DistanceMatrix::from_points(&grid(9)).unwrap();
    assert_registry_is_worker_invariant(&dist, "grid(9x9)");
}