breakout 0.4.0

Breakout detection for Rust
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291

use crate::math::pow2;
use crate::multiset::Multiset;
use alloc::{vec, vec::Vec};
use core::cmp::Ordering;

#[derive(PartialEq)]
struct MaxItem(f64);

impl Eq for MaxItem {}

// opposite of edmx
impl PartialOrd for MaxItem {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for MaxItem {
    fn cmp(&self, other: &MaxItem) -> Ordering {
        other.0.partial_cmp(&self.0).unwrap()
    }
}

#[derive(PartialEq)]
struct MinItem(f64);

impl Eq for MinItem {}

// opposite of edmx
impl PartialOrd for MinItem {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for MinItem {
    fn cmp(&self, other: &MinItem) -> Ordering {
        self.0.partial_cmp(&other.0).unwrap()
    }
}

fn insert_element(m: &mut Multiset<MinItem>, m2: &mut Multiset<MaxItem>, x: f64) {
    if m.is_empty() || x < m.first().unwrap().0 {
        m2.insert(MaxItem(x));
    } else {
        m.insert(MinItem(x));
    }

    if m.len() > m2.len() + 1 {
        let i = m.first().unwrap().0;
        m2.insert(MaxItem(i));
        m.remove(MinItem(i));
    } else if m2.len() > m.len() + 1 {
        let i = m2.first().unwrap().0;
        m.insert(MinItem(i));
        m2.remove(MaxItem(i));
    }
}

fn remove_element(m: &mut Multiset<MinItem>, m2: &mut Multiset<MaxItem>, x: f64) {
    if x < m.first().unwrap().0 {
        m2.remove(MaxItem(x));
    } else {
        m.remove(MinItem(x));
    }

    if m.len() > m2.len() + 1 {
        let i = m.first().unwrap().0;
        m2.insert(MaxItem(i));
        m.remove(MinItem(i));
    } else if m2.len() > m.len() + 1 {
        let i = m2.first().unwrap().0;
        m.insert(MinItem(i));
        m2.remove(MaxItem(i));
    }
}

// given a pair of trees obtain the median
fn get_median(m: &Multiset<MinItem>, m2: &Multiset<MaxItem>) -> f64 {
    match m.len().cmp(&m2.len()) {
        Ordering::Greater => m.first().unwrap().0,
        Ordering::Less => m2.first().unwrap().0,
        Ordering::Equal => (m2.first().unwrap().0 + m.first().unwrap().0) / 2.0,
    }
}

fn linear(_x: f64) -> f64 {
    1.0
}

fn constant(_x: f64) -> f64 {
    0.0
}

fn quadratic(x: f64) -> f64 {
    2.0 * x + 1.0
}

pub fn edm_multi(z: &[f64], min_size: usize, beta: f64, degree: i32) -> Vec<usize> {
    // identify which type of penalization to use
    let g: fn(f64) -> f64 = match degree {
        1 => linear,
        2 => quadratic,
        _ => constant,
    };

    let n = z.len();
    let mut beta = beta;
    // assume that beta is a positive number
    if beta < 0.0 {
        beta = -beta;
    }
    let mut prev = vec![0; n + 1];
    let mut number = vec![0; n + 1];
    let mut f = vec![-3.0; n + 1];

    // trees used to store the "upper half" of the considered observations
    let mut right_min = Multiset::new();
    let mut left_min = Multiset::new();

    // trees used to store the "lower half" of the considered observations
    let mut right_max = Multiset::new();
    let mut left_max = Multiset::new();

    // iterate over possible locations for the last change
    for s in 2 * min_size..n + 1 {
        right_max.clear();
        right_min.clear();
        left_max.clear();
        left_min.clear();

        // initialize left and right trees to account for minimum segment size
        for i in prev[min_size - 1]..min_size - 1 {
            insert_element(&mut left_min, &mut left_max, z[i]);
        }
        for i in min_size - 1..s {
            insert_element(&mut right_min, &mut right_max, z[i]);
        }

        // iterate over possible locations for the penultimate change
        // modify limits to deal with min_size
        for t in min_size..s - min_size + 1 {
            // insert element into left tree
            insert_element(&mut left_min, &mut left_max, z[t - 1]);
            // remove element from right tree
            remove_element(&mut right_min, &mut right_max, z[t - 1]);

            // left tree now has { Z[prev[t-1]], ..., Z[t-1] }
            // right tree now has { Z[t], ..., Z[s-1] }

            // check to see if optimal position of previous change point has changed
            // if so update the left tree
            if prev[t] > prev[t - 1] {
                for i in prev[t - 1]..prev[t] {
                    remove_element(&mut left_min, &mut left_max, z[i]);
                }
            } else if prev[t] < prev[t - 1] {
                for i in prev[t]..prev[t - 1] {
                    insert_element(&mut left_min, &mut left_max, z[i]);
                }
            }

            // calculate statistic value
            let left_median = get_median(&left_min, &left_max);
            let right_median = get_median(&right_min, &right_max);
            let normalize = ((t - prev[t]) * (s - t)) as f64 / pow2((s - prev[t]) as f64);
            let tmp =
                f[t] + normalize * pow2(left_median - right_median) - beta * g(number[t] as f64);

            // check for improved optimal statistic value
            if tmp > f[s] {
                number[s] = number[t] + 1;
                f[s] = tmp;
                prev[s] = t;
            }
        }
    }

    // obtain list of optimal change point estimates
    let mut ret = Vec::new();
    let mut at = n;
    while at != 0 {
        // don't insert 0 as a change point estimate
        if prev[at] != 0 {
            ret.push(prev[at]);
        }
        at = prev[at];
    }
    ret.sort_unstable();

    ret
}

// Penalizes based on percent change in the statistic value.
// Linear penalty means that each new breakout must result in an at least X% increase
// Quadratic penalty means that each new breakout must result in at least an (X*k)% increase for k breakouts
pub fn edm_percent(z: &[f64], min_size: usize, percent: f64, degree: i32) -> Vec<usize> {
    // identify which type of penalization to use
    let g: fn(f64) -> f64 = match degree {
        1 => linear,
        2 => quadratic,
        _ => constant,
    };

    let n = z.len();
    let mut prev = vec![0; n + 1];
    let mut number = vec![0; n + 1];
    let mut f = vec![0.0; n + 1];

    // trees used to store the "upper half" of the considered observations
    let mut right_min = Multiset::new();
    let mut left_min = Multiset::new();

    // trees used to store the "lower half" of the considered observations
    let mut right_max = Multiset::new();
    let mut left_max = Multiset::new();

    // iterate over possible locations for the last change
    for s in 2 * min_size..n + 1 {
        right_max.clear();
        right_min.clear();
        left_max.clear();
        left_min.clear();

        // initialize left and right trees to account for minimum segment size
        for i in prev[min_size - 1]..min_size - 1 {
            insert_element(&mut left_min, &mut left_max, z[i]);
        }
        for i in min_size - 1..s {
            insert_element(&mut right_min, &mut right_max, z[i]);
        }

        // iterate over possible locations for the penultiamte change
        // modify limits to deal with min_size
        for t in min_size..s - min_size + 1 {
            // insert element into left tree
            insert_element(&mut left_min, &mut left_max, z[t - 1]);
            // remove element from right tree
            remove_element(&mut right_min, &mut right_max, z[t - 1]);

            // left tree now has { Z[prev[t-1]], ..., Z[t-1] }
            // right tree now has { Z[t], ..., Z[s-1] }

            // check to see if optimal position of previous change point has changed
            // if so update the left tree
            if prev[t] > prev[t - 1] {
                for i in prev[t - 1]..prev[t] {
                    remove_element(&mut left_min, &mut left_max, z[i]);
                }
            } else if prev[t] < prev[t - 1] {
                for i in prev[t]..prev[t - 1] {
                    insert_element(&mut left_min, &mut left_max, z[i]);
                }
            }

            // calculate statistic value
            let left_median = get_median(&left_min, &left_max);
            let right_median = get_median(&right_min, &right_max);
            let normalize = ((t - prev[t]) * (s - t)) as f64 / pow2((s - prev[t]) as f64);
            let tmp = f[t] + normalize * pow2(left_median - right_median);

            // find best location for change point. check % condition later
            if tmp > f[s] {
                number[s] = number[t] + 1;
                f[s] = tmp;
                prev[s] = t;
            }
        }

        // check to make sure we meet the percent change requirement
        if prev[s] != 0 && f[s] - f[prev[s]] < percent * g(number[prev[s]] as f64) * f[prev[s]] {
            number[s] = number[prev[s]];
            f[s] = f[prev[s]];
            prev[s] = prev[prev[s]];
        }
    }

    // obtain list of optimal change point estimates
    let mut ret = Vec::new();
    let mut at = n;
    while at != 0 {
        // don't insert 0 as a change point estimate
        if prev[at] != 0 {
            ret.push(prev[at]);
        }
        at = prev[at];
    }
    ret.sort_unstable();

    ret
}