psqr/
lib.rs

1#[derive(Default)]
2pub struct P2 {
3    quantile: f64,
4
5    heights: Vec<f64>,
6    pos: [i64; 5],
7    n_pos: [f64; 5],
8    dn: [f64; 5],
9
10    filled: bool,
11}
12
13impl P2 {
14    pub fn new(quantile: f64) -> Self {
15        if quantile < 0.0 || quantile > 1.0 {
16            panic!("quantile must be in [0, 1]");
17        }
18        let mut p2 = Self {
19            quantile,
20            n_pos: [
21                0.0,
22                2.0 * quantile,
23                4.0 * quantile,
24                2.0 + 2.0 * quantile,
25                4.0,
26            ],
27            dn: [0.0, quantile / 2.0, quantile, (1.0 + quantile) / 2.0, 1.0],
28            ..Default::default()
29        };
30
31        for i in 0..p2.pos.len() {
32            p2.pos[i] = i as i64;
33        }
34
35        p2
36    }
37
38    pub fn append(&mut self, data: f64) {
39        if self.heights.len() != 5 {
40            self.heights.push(data);
41            return;
42        }
43        if !self.filled {
44            self.filled = true;
45            self.heights.sort_by(|a, b| a.partial_cmp(b).unwrap());
46        }
47        self.append_data(data);
48    }
49
50    fn append_data(&mut self, data: f64) {
51        let l = self.heights.len() - 1;
52        let mut k: isize = -1;
53        if data < self.heights[0] {
54            k = 0;
55            self.heights[0] = data;
56        } else if self.heights[l] <= data {
57            k = l as isize - 1;
58            self.heights[l] = data;
59        } else {
60            for i in 1..=l {
61                if self.heights[i - 1] <= data && data < self.heights[i] {
62                    k = i as isize - 1;
63                    break;
64                }
65            }
66        }
67        for i in 0..self.pos.len() {
68            if i > k as usize {
69                self.pos[i] += 1;
70            }
71            self.n_pos[i] += self.dn[i];
72        }
73
74        self.adjust_heights();
75    }
76
77    fn adjust_heights(&mut self) {
78        for i in 1..self.heights.len() - 1 {
79            let n = self.pos[i];
80            let np1 = self.pos[i + 1];
81            let nm1 = self.pos[i - 1];
82
83            let d = self.n_pos[i] - n as f64;
84
85            if (d >= 1.0 && np1 - n > 1) || (d <= -1.0 && nm1 - n < -1) {
86                let d = if d >= 0.0 { 1.0 } else { -1.0 };
87
88                let h = self.heights[i];
89                let hp1 = self.heights[i + 1];
90                let hm1 = self.heights[i - 1];
91
92                let hi = parabolic(d, hp1, h, hm1, np1 as f64, n as f64, nm1 as f64);
93
94                if hm1 < hi && hi < hp1 {
95                    self.heights[i] = hi;
96                } else {
97                    let hd = self.heights[i + d as usize];
98                    let nd: i64 = self.pos[i + d as usize];
99                    self.heights[i] = h + d * (hd - h) / (nd - n) as f64;
100                }
101
102                self.pos[i] += d as i64;
103            }
104        }
105    }
106
107    pub fn value(&mut self) -> f64 {
108        if !self.filled {
109            let l = self.heights.len();
110            match l {
111                0 => return 0.0,
112                1 => return self.heights[0],
113                _ => self.heights.sort_by(|a, b| a.partial_cmp(b).unwrap()),
114            }
115            let rank = (self.quantile * l as f64) as usize;
116            return self.heights[rank];
117        }
118        self.heights[2]
119    }
120}
121
122fn parabolic(d: f64, qp1: f64, q: f64, qm1: f64, np1: f64, n: f64, nm1: f64) -> f64 {
123    let a = d / (np1 - nm1);
124    let b1 = (n - nm1 + d) * (qp1 - q) / (np1 - n);
125    let b2 = (np1 - n - d) * (q - qm1) / (n - nm1);
126    q + a * (b1 + b2)
127}
128
129#[cfg(test)]
130mod test {
131    use super::P2;
132
133    #[test]
134    fn test_p2() {
135        let mut p2 = P2::new(0.3);
136
137        for i in 1..=100 {
138            p2.append(i as f64);
139        }
140
141        let x = p2.value();
142        assert_eq!(x, 30.0);
143    }
144}
psqr/lib.rs

psqr/
lib.rs
