fast-nnt 0.4.0

Rust implementation of the NeighborNet phylogenetic network algorithm.
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

use crate::splits::asplit::ASplit;
use fixedbitset::FixedBitSet;
use ndarray::Array2;

/// Least-squares fit of the splits to the distances, as a percentage:
/// 100 * (1 - sum((s_ij - d_ij)^2) / sum(d_ij^2)) over i<j.
pub fn compute_least_squares_fit(distances: &Array2<f64>, splits: &[ASplit]) -> f32 {
    let n = distances.nrows();
    assert_eq!(n, distances.ncols(), "distances must be square");

    let tri_len = n * (n - 1) / 2;
    let mut split_dist = vec![0.0_f64; tri_len];
    for s in splits {
        let w = s.weight;
        let a: &FixedBitSet = s.get_a();
        let b: &FixedBitSet = s.get_b();

        for i1 in a.ones() {
            if i1 == 0 || i1 > n {
                continue;
            }
            let ii = i1 - 1;
            for j1 in b.ones() {
                if j1 == 0 || j1 > n {
                    continue;
                }
                let jj = j1 - 1;
                let (i, j) = if ii < jj { (ii, jj) } else { (jj, ii) };
                let idx = pair_index(i, j, n);
                split_dist[idx] += w;
            }
        }
    }

    let mut sum_diff_sq = 0.0_f64;
    let mut sum_d_sq = 0.0_f64;
    for i in 0..n - 1 {
        let start = row_start(i, n);
        for j in (i + 1)..n {
            let idx = start + (j - i - 1);
            let sij = split_dist[idx];
            let dij = distances[[i, j]];
            let diff = sij - dij;
            sum_diff_sq += diff * diff;
            sum_d_sq += dij * dij;
        }
    }

    let fit = if sum_d_sq > 0.0 {
        100.0 * (1.0 - (sum_diff_sq / sum_d_sq))
    } else {
        0.0
    };

    fit as f32
}

#[inline]
fn row_start(i: usize, n: usize) -> usize {
    i * (2 * n - i - 1) / 2
}

#[inline]
fn pair_index(i: usize, j: usize, n: usize) -> usize {
    debug_assert!(i < j && j < n);
    row_start(i, n) + (j - i - 1)
}