slain-num 0.3.0

SLAIN's own numeric crate: dense matrices, least squares, polynomials, FFT, convolution, statistics and 1-D/N-D minimisation, with AVX2 kernels.
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

use crate::simd;

#[must_use]
pub fn filter_valid(x: &[f64], taps: &[f64]) -> Vec<f64> {
    if taps.is_empty() || x.len() < taps.len() {
        return Vec::new();
    }
    let n = x.len() - taps.len() + 1;
    let mut out = vec![0.0; n];
    for (k, &t) in taps.iter().enumerate() {
        simd::axpy(t, &x[k..k + n], &mut out);
    }
    out
}

#[must_use]
pub fn filter_same(x: &[f64], taps: &[f64]) -> Vec<f64> {
    assert!(taps.len() % 2 == 1, "centred filtering needs an odd tap count");
    let r = taps.len() / 2;
    if x.is_empty() {
        return Vec::new();
    }
    let mut padded = Vec::with_capacity(x.len() + 2 * r);
    padded.extend(std::iter::repeat_n(x[0], r));
    padded.extend_from_slice(x);
    padded.extend(std::iter::repeat_n(x[x.len() - 1], r));
    filter_valid(&padded, taps)
}

#[must_use]
pub fn filter_2d(img: &[f64], w: usize, h: usize, htaps: &[f64], vtaps: &[f64]) -> Vec<f64> {
    let mut tmp = vec![0.0; w * h];
    for r in 0..h {
        tmp[r * w..(r + 1) * w].copy_from_slice(&filter_same(&img[r * w..(r + 1) * w], htaps));
    }
    let rv = vtaps.len() / 2;
    let mut out = vec![0.0; w * h];
    for r in 0..h {
        let orow = &mut out[r * w..(r + 1) * w];
        for (k, &t) in vtaps.iter().enumerate() {
            let src = (r + k).saturating_sub(rv).min(h - 1);
            simd::axpy(t, &tmp[src * w..(src + 1) * w], orow);
        }
    }
    out
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn filters_match_naive() {
        let x: Vec<f64> = (0..37).map(|i| (i as f64 * 0.5).sin()).collect();
        let taps = [0.25, 0.5, 0.25];
        let got = filter_valid(&x, &taps);
        for (i, g) in got.iter().enumerate() {
            let want: f64 = taps.iter().enumerate().map(|(k, t)| t * x[i + k]).sum();
            assert!((g - want).abs() < 1e-12);
        }
        let same = filter_same(&x, &taps);
        assert_eq!(same.len(), x.len());
        assert!((same[0] - (0.75 * x[0] + 0.25 * x[1])).abs() < 1e-12);
        let img = vec![1.0; 6 * 5];
        let f = filter_2d(&img, 6, 5, &taps, &taps);
        assert!(f.iter().all(|v| (v - 1.0).abs() < 1e-12));
    }
}