nnls 0.5.0

Non-Negative Least Squares (nnls)
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

use approx::assert_relative_eq;
use ndarray::{arr1, arr2, Array1, Array2};

use nnls::nnls;

#[test]
fn test_nnls_1() {
    let a = arr2(&[[1.0, 0.0], [1.0, 0.0], [0.0, 1.0]]);
    let b = arr1(&[2.0, 1.0, 1.0]);
    let (x, r_norm) = nnls(a.view(), b.view());
    assert_relative_eq!(x, arr1(&[1.5, 1.0]), epsilon = 1e-9);
    assert_relative_eq!(r_norm, 0.7071067811865475, epsilon = 1e-9);
}

#[test]
fn test_nnls_2() {
    let a = arr2(&[[1.0, 0.5, 0.0], [1.0, 0.6, 0.0], [0.9, 0.5, 0.1], [0.0, 0.0, 1.0]]);
    let b = arr1(&[2.0, 1.0, 1.0, 0.5]);
    let (x, r_norm) = nnls(a.view(), b.view());
    assert_relative_eq!(x, arr1(&[1.37279152, 0.0, 0.47173145]), epsilon = 1e-5);
    assert_relative_eq!(r_norm, 0.7829905522718629, epsilon = 1e-5);
}

#[test]
fn test_trap() {
    let a = arr2(&[
        [9.6110e-8, 7.9576e-17, 2.7234e-44, 3.4490e-44],
        [7.9576e-17, 2.0449e-7, 7.4360e-44, 9.4173e-44],
        [2.7234e-44, 7.4360e-44, 3.1651e-7, 2.9014e-44],
        [3.4490e-44, 9.4173e-44, 2.9014e-44, 4.3443e-7],
    ]);
    let b = arr1(&[-2.4678e-17, -6.7118e-17, 3.6306e-44, 4.5923e-44]);
    let (x, r_norm) = nnls(a.view(), b.view());
    assert_relative_eq!(x, arr1(&[0.0, 0.0, 0.0, 0.0]), epsilon = 1e-5);
    assert_relative_eq!(r_norm, 0.0, epsilon = 1e-5);
}

#[test]
fn test_mixture_of_exponentials() {
    let mut t0 = Array1::zeros(19);
    t0[0] = 2.0;
    t0[1] = t0[0] * 2.0f64.sqrt();
    for i in 2..t0.len() {
        t0[i] = 2.0 * t0[i - 2];
    }

    let arg_limit = (2.0 * 2.2250738585072014e-308f64).ln();
    let a = Array2::from_shape_fn((100, 19), |(i, j)| {
        let i = (i + 1) as f64;
        let arg = -10.0 * i / t0[j];
        if arg > arg_limit {
            arg.exp()
        } else {
            0.0
        }
    });

    let b = Array1::from_shape_fn(100, |i| {
        let t = 10.0 * (i + 1) as f64;
        100.0 * ((-t / 5.0).exp() + (-t / 50.0).exp() + (-t / 500.0).exp()) - 0.5
    });

    let (x, r_norm) = nnls(a.view(), b.view());

    let mut x_gt = Array1::zeros(19);
    x_gt[2] = 146.18289932114504;
    x_gt[9] = 82.906587033240356;
    x_gt[10] = 18.699400056374593;
    x_gt[15] = 13.315689103364351;
    x_gt[16] = 87.197680367574506;
    assert_relative_eq!(x, x_gt, epsilon = 1e-5);
    assert_relative_eq!(r_norm, 0.88433322916388879, epsilon = 1e-5);
}