pelt 0.3.1

Changepoint detection with Pruned Exact Linear Time
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

//! L1 cost functions.

use std::ops::Range;

use ndarray::{ArrayView1, ArrayView2};

use crate::cost::tree::KthSmallestTree;

/// Precalculation output.
pub struct L1Cost1D {
    /// Tree for finding the mean (middle K-th smallest).
    kth_smallest_tree: KthSmallestTree,
}

impl L1Cost1D {
    /// Build the K-th smallest tree.
    #[inline]
    pub(crate) fn precalculate(signal: &ArrayView1<f64>) -> Self {
        // Build the tree from the signal
        let kth_smallest_tree = KthSmallestTree::build(signal);

        Self { kth_smallest_tree }
    }

    /// Calculate the loss.
    #[inline]
    pub(crate) fn loss(&self, total_loss: &mut f64, signal: &ArrayView1<f64>, range: Range<usize>) {
        // Calculate the median for the segment
        let median = self.median(range.clone());

        // Take the sub slice of the 2D object
        signal
            .slice(ndarray::s!(range))
            // Calculate the absolute difference for each point with the median
            .iter()
            // Sum the values
            .for_each(|signal| *total_loss += (*signal - median).abs());
    }

    /// Get the median of a range in the signal.
    #[inline]
    fn median(&self, range: Range<usize>) -> f64 {
        // Find the midpoint of the range
        let len = range.len();

        // Convert the range to an inclusive one
        let range_inclusive = range.start..=(range.end - 1);

        // Calculate median based one one or two variables if it's even
        if len.is_multiple_of(2) {
            // Get the middle value, offset by 1 for the two points
            let kth = len / 2;

            // Get the average of the two median
            let median1 = self.kth_smallest_tree.kth(range_inclusive.clone(), kth);
            let median2 = self.kth_smallest_tree.kth(range_inclusive, kth + 1);

            median1.midpoint(median2)
        } else {
            // Get the middle value
            let kth = len / 2 + 1;

            // Get the median
            self.kth_smallest_tree.kth(range_inclusive, kth)
        }
    }
}

/// Precalculation output.
pub struct L1Cost2D {
    /// Precalculated per column.
    columns: Vec<L1Cost1D>,
}

impl L1Cost2D {
    /// Precalculate the sum queries.
    #[inline]
    pub fn precalculate(signal: &ArrayView2<f64>) -> Self {
        let columns = signal
            .columns()
            .into_iter()
            .map(|column| L1Cost1D::precalculate(&column))
            .collect();

        Self { columns }
    }

    /// Calculate the loss.
    ///
    /// Calculated using Welford's algorithm.
    #[inline]
    pub(crate) fn loss(&self, total_loss: &mut f64, signal: &ArrayView2<f64>, range: Range<usize>) {
        // Calculate total loss
        self.columns
            .iter()
            .zip(signal.columns())
            .for_each(|(column, signal_column)| {
                column.loss(total_loss, &signal_column, range.clone())
            })
    }
}

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

    /// Check the L1 cost function.
    #[test]
    fn cost_1d() {
        let array_1d = ndarray::array![10.0, 30.0, 20.0];
        let cost = L1Cost1D::precalculate(&array_1d.view());
        let mut loss = 0.0;
        cost.loss(&mut loss, &array_1d.view(), 0..3);
        assert_eq!(loss, 20.0);
    }

    /// Check the L1 cost function.
    #[test]
    fn cost_2d() {
        let array_2d = ndarray::array![[10.0], [30.0], [20.0]];
        let cost = L1Cost2D::precalculate(&array_2d.view());
        let mut loss = 0.0;
        cost.loss(&mut loss, &array_2d.view(), 0..3);
        assert_eq!(loss, 20.0);
    }
}