altos: Perform time comparisons using 16-bit arithmetic to handle wrap
[fw/altos] / altoslib / AltosTimeSeries.java
1 /*
2  * Copyright © 2017 Keith Packard <keithp@keithp.com>
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful, but
10  * WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12  * General Public License for more details.
13  */
14
15 package org.altusmetrum.altoslib_12;
16
17 import java.util.*;
18
19 public class AltosTimeSeries implements Iterable<AltosTimeValue>, Comparable<AltosTimeSeries> {
20         public String                   label;
21         public AltosUnits               units;
22         ArrayList<AltosTimeValue>       values;
23
24         public int compareTo(AltosTimeSeries other) {
25                 return label.compareTo(other.label);
26         }
27
28         public void add(AltosTimeValue tv) {
29                 values.add(tv);
30         }
31
32         public void add(double time, double value) {
33                 add(new AltosTimeValue(time, value));
34         }
35
36         public AltosTimeValue get(int i) {
37                 return values.get(i);
38         }
39
40         private double lerp(AltosTimeValue v0, AltosTimeValue v1, double t) {
41                 /* degenerate case */
42                 if (v0.time == v1.time)
43                         return (v0.value + v1.value) / 2;
44
45                 return (v0.value * (v1.time - t) + v1.value * (t - v0.time)) / (v1.time - v0.time);
46         }
47
48         private int after_index(double time) {
49                 int     lo = 0;
50                 int     hi = values.size() - 1;
51
52                 while (lo <= hi) {
53                         int mid = (lo + hi) / 2;
54
55                         if (values.get(mid).time < time)
56                                 lo = mid + 1;
57                         else
58                                 hi = mid - 1;
59                 }
60                 return lo;
61         }
62
63         /* Compute a value for an arbitrary time */
64         public double value(double time) {
65                 int after = after_index(time);
66                 double ret;
67
68                 if (after == 0)
69                         ret = values.get(0).value;
70                 else if (after == values.size())
71                         ret = values.get(after - 1).value;
72                 else {
73                         AltosTimeValue b = values.get(after-1);
74                         AltosTimeValue a = values.get(after);
75                         ret = lerp(b, a, time);
76                 }
77                 return ret;
78         }
79
80         /* Find the value just before an arbitrary time */
81         public double value_before(double time) {
82                 int after = after_index(time);
83
84                 if (after == 0)
85                         return values.get(0).value;
86                 return values.get(after-1).value;
87         }
88
89         /* Find the value just after an arbitrary time */
90         public double value_after(double time) {
91                 int after = after_index(time);
92
93                 if (after == values.size())
94                         return values.get(after-1).value;
95                 return values.get(after).value;
96         }
97
98         public double time_of(double value) {
99                 double  last = AltosLib.MISSING;
100                 for (AltosTimeValue v : values) {
101                         if (v.value >= value)
102                                 return v.time;
103                         last = v.time;
104                 }
105                 return last;
106         }
107
108         public int size() {
109                 return values.size();
110         }
111
112         public Iterator<AltosTimeValue> iterator() {
113                 return values.iterator();
114         }
115
116         public AltosTimeValue max() {
117                 AltosTimeValue max = null;
118                 for (AltosTimeValue tv : values)
119                         if (max == null || tv.value > max.value)
120                                 max = tv;
121                 return max;
122         }
123
124         public AltosTimeValue max(double start_time, double end_time) {
125                 AltosTimeValue max = null;
126                 for (AltosTimeValue tv : values) {
127                         if (start_time <= tv.time && tv.time <= end_time)
128                                 if (max == null || tv.value > max.value)
129                                         max = tv;
130                 }
131                 return max;
132         }
133
134         public AltosTimeValue min() {
135                 AltosTimeValue min = null;
136                 for (AltosTimeValue tv : values) {
137                         if (min == null || tv.value < min.value)
138                                 min = tv;
139                 }
140                 return min;
141         }
142
143         public AltosTimeValue min(double start_time, double end_time) {
144                 AltosTimeValue min = null;
145                 for (AltosTimeValue tv : values) {
146                         if (start_time <= tv.time && tv.time <= end_time)
147                                 if (min == null || tv.value < min.value)
148                                         min = tv;
149                 }
150                 return min;
151         }
152
153         public AltosTimeValue first() {
154                 return values.get(0);
155         }
156
157         public AltosTimeValue last() {
158                 return values.get(values.size() - 1);
159         }
160
161         public double average() {
162                 double total_value = 0;
163                 double total_time = 0;
164                 AltosTimeValue prev = null;
165                 for (AltosTimeValue tv : values) {
166                         if (prev != null) {
167                                 total_value += (tv.value + prev.value) / 2 * (tv.time - prev.time);
168                                 total_time += (tv.time - prev.time);
169                         }
170                         prev = tv;
171                 }
172                 if (total_time == 0)
173                         return AltosLib.MISSING;
174                 return total_value / total_time;
175         }
176
177         public double average(double start_time, double end_time) {
178                 double total_value = 0;
179                 double total_time = 0;
180                 AltosTimeValue prev = null;
181                 for (AltosTimeValue tv : values) {
182                         if (start_time <= tv.time && tv.time <= end_time) {
183                                 if (prev != null) {
184                                         total_value += (tv.value + prev.value) / 2 * (tv.time - start_time);
185                                         total_time += (tv.time - start_time);
186                                 }
187                                 start_time = tv.time;
188                         }
189                         prev = tv;
190                 }
191                 if (total_time == 0)
192                         return AltosLib.MISSING;
193                 return total_value / total_time;
194         }
195
196         public AltosTimeSeries integrate(AltosTimeSeries integral) {
197                 double  value = 0.0;
198                 double  pvalue = 0.0;
199                 double  time = 0.0;
200                 boolean start = true;
201
202                 for (AltosTimeValue v : values) {
203                         if (start) {
204                                 value = 0.0;
205                                 start = false;
206                         } else {
207                                 value += (pvalue + v.value) / 2.0 * (v.time - time);
208                         }
209                         pvalue = v.value;
210                         time = v.time;
211                         integral.add(time, value);
212
213                 }
214                 return integral;
215         }
216
217         public AltosTimeSeries differentiate(AltosTimeSeries diff) {
218                 double value = 0.0;
219                 double time = 0.0;
220                 boolean start = true;
221
222                 for (AltosTimeValue v: values) {
223                         if (start) {
224                                 value = v.value;
225                                 time = v.time;
226                                 start = false;
227                         } else {
228                                 double  dx = v.time - time;
229                                 double  dy = v.value - value;
230
231                                 if (dx != 0)
232                                         diff.add(time, dy/dx);
233
234                                 time = v.time;
235                                 value = v.value;
236                         }
237                 }
238                 return diff;
239         }
240
241         private int find_left(int i, double dt) {
242                 int j;
243                 double t = values.get(i).time - dt;
244                 for (j = i; j >= 0; j--)        {
245                         if (values.get(j).time < t)
246                                 break;
247                 }
248                 return j + 1;
249
250         }
251
252         private int find_right(int i, double dt) {
253                 int j;
254                 double t = values.get(i).time + dt;
255                 for (j = i; j < values.size(); j++) {
256                         if (values.get(j).time > t)
257                                 break;
258                 }
259                 return j - 1;
260
261         }
262
263         private double filter_coeff(double dist, double width) {
264                 double ratio = dist / (width / 2);
265
266                 return Math.cos(ratio * Math.PI / 2);
267         }
268
269         public AltosTimeSeries filter(AltosTimeSeries f, double width) {
270                 double  half_width = width/2;
271                 for (int i = 0; i < values.size(); i++) {
272                         double  center_time = values.get(i).time;
273                         double  left_time = center_time - half_width;
274                         double  right_time = center_time + half_width;
275                         double  total_coeff = 0.0;
276                         double  total_value = 0.0;
277
278                         int     left = find_left(i, half_width);
279                         int     right = find_right(i, half_width);
280
281                         for (int j = left; j <= right; j++) {
282                                 double  j_time = values.get(j).time;
283
284                                 if (left_time <= j_time && j_time <= right_time) {
285                                         double  j_left = j == left ? left_time : values.get(j-1).time;
286                                         double  j_right = j == right ? right_time : values.get(j+1).time;
287                                         double  interval = (j_right - j_left) / 2.0;
288                                         double  coeff = filter_coeff(j_time - center_time, width) * interval;
289                                         double  value = values.get(j).value;
290                                         double  partial = value * coeff;
291
292                                         total_coeff += coeff;
293                                         total_value += partial;
294                                 }
295                         }
296                         if (total_coeff != 0.0)
297                                 f.add(center_time, total_value / total_coeff);
298                 }
299                 return f;
300         }
301
302         public AltosTimeSeries(String label, AltosUnits units) {
303                 this.label = label;
304                 this.units = units;
305                 this.values = new ArrayList<AltosTimeValue>();
306         }
307 }