1 package org.opentrafficsim.draw.egtf;
2
3 import java.util.LinkedHashMap;
4 import java.util.LinkedHashSet;
5 import java.util.Map;
6 import java.util.NavigableMap;
7 import java.util.Objects;
8 import java.util.Optional;
9 import java.util.Set;
10 import java.util.SortedMap;
11 import java.util.TreeMap;
12 import java.util.stream.IntStream;
13
14 /**
15 * Extended Generalized Treiber-Helbing Filter (van Lint and Hoogendoorn, 2009). This is an extension of the Adaptive Smoothing
16 * Method (Treiber and Helbing, 2002). A fast filter for equidistant grids (Schreiter et al., 2010) is available. This fast
17 * implementation also supports multiple data sources.
18 * <p>
19 * To allow flexible usage the EGTF works with {@code DataSource}, {@code DataStream} and {@code Quantity}.
20 * <p>
21 * A {@code DataSource}, such as "loop detectors", "floating-car data" or "camera" is mostly an identifier, but can be requested
22 * to provide several data streams.
23 * <p>
24 * A {@code DataStream} is one {@code DataSource} supplying one {@code Quantity}. For instance "loop detectors" supplying
25 * "flow". In a {@code DataStream}, supplied by the {@code DataSource}, standard deviation of measurements in congestion and
26 * free flow are defined. These determine the reliability of the {@code Quantity} data from the given {@code DataSource}, and
27 * thus ultimately the weight of the data in the estimation of the quantity.
28 * <p>
29 * A {@code Quantity}, such as "flow" or "density" defines what is measured and what is requested as output. The output can be
30 * in typed format using a {@code Converter}. Default quantities are available under {@code SPEED_SI}, {@code FLOW_SI} and
31 * {@code DENSITY_SI}, all under {@code Quantity}.
32 * <p>
33 * Data can be added using several methods for point data, vector data (multiple independent location-time combinations) and
34 * grid data (data in a grid pattern). Data is added for a particular {@code DataStream}.
35 * <p>
36 * For simple use-cases where a single data source is used, data can be added directly with a {@code Quantity}, in which case a
37 * default {@code DataSource}, and default {@code DataStream} for each {@code Quantity} is internally used. All data should
38 * either be added using {@code Quantity}'s, or it should all be added using {@code DataSource}'s. Otherwise relative data
39 * reliability is undefined.
40 * <p>
41 * Output can be requested from the EGTF using a {@code Kernel}, a spatiotemporal pattern determining measurement weights. The
42 * {@code Kernel} defines an optional maximum spatial and temporal range for measurements to consider, and uses a {@code Shape}
43 * to determine the weight for a given distance and time from the estimated point. By default this is an exponential function. A
44 * Gaussian kernel is also available, while any other shape could be also be implemented.
45 * <p>
46 * Parameters from the EGTF are found in the following places:
47 * <ul>
48 * <li>{@code EGTF}: <i>cCong</i>, <i>cFree</i>, <i>deltaV</i> and <i>vc</i>, defining the overall traffic flow properties.</li>
49 * <li>{@code Kernel}: <i>tMax</i> and <i>xMax</i>, defining the maximum range to consider.</li>
50 * <li>{@code KernelShape}: <i>sigma</i> and <i>tau</i>, determining the decay of weights for further measurements in space and
51 * time. (Specifically {@code GaussKernelShape})</li>
52 * <li>{@code DataStream}: <i>thetaCong</i> and <i>thetaFree</i>, defining the reliability by the standard deviation of measured
53 * data in free flow and congestion from a particular data stream.</li>
54 * </ul>
55 * References:
56 * <ul>
57 * <li>van Lint, J. W. C. and Hoogendoorn, S. P. (2009). A robust and efficient method for fusing heterogeneous data from
58 * traffic sensors on freeways. Computer Aided Civil and Infrastructure Engineering, accepted for publication.</li>
59 * <li>Schreiter, T., van Lint, J. W. C., Treiber, M. and Hoogendoorn, S. P. (2010). Two fast implementations of the Adaptive
60 * Smoothing Method used in highway traffic state estimation. 13th International IEEE Conference on Intelligent Transportation
61 * Systems, 19-22 Sept. 2010, Funchal, Portugal.</li>
62 * <li>Treiber, M. and Helbing, D. (2002). Reconstructing the spatio-temporal traffic dynamics from stationary detector data.
63 * Cooper@tive Tr@nsport@tion Dyn@mics, 1:3.1–3.24.</li>
64 * </ul>
65 * <p>
66 * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
67 * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
68 * </p>
69 * @author <a href="https://github.com/wjschakel">Wouter Schakel</a>
70 */
71 public class Egtf
72 {
73
74 /** Default sigma value. */
75 private static final double DEFAULT_SIGMA = 300.0;
76
77 /** Default tau value. */
78 private static final double DEFAULT_TAU = 30.0;
79
80 /** Filter kernel. */
81 private Kernel kernel;
82
83 /** Shock wave speed in congestion. */
84 private final double cCong;
85
86 /** Shock wave speed in free flow. */
87 private final double cFree;
88
89 /** Speed range between congestion and free flow. */
90 private final double deltaV;
91
92 /** Flip-over speed below which we have congestion. */
93 private final double vc;
94
95 /** Data sources by label so we can return the same instances upon repeated request. */
96 private final Map<String, DataSource> dataSources = new LinkedHashMap<>();
97
98 /** Default data source for cases where a single data source is used. */
99 private DataSource defaultDataSource = null;
100
101 /** Default data streams for cases where a single data source is used. */
102 private Map<Quantity<?, ?>, DataStream<?>> defaultDataStreams = null;
103
104 /** True if data is currently being added using a quantity, in which case a check should not occur. */
105 private boolean addingByQuantity;
106
107 /** All point data sorted by space and time, and per data stream. */
108 private NavigableMap<Double, NavigableMap<Double, Map<DataStream<?>, Double>>> data = new TreeMap<>();
109
110 /** Whether the calculation was interrupted. */
111 private boolean interrupted = false;
112
113 /** Listeners. */
114 private Set<EgtfListener> listeners = new LinkedHashSet<>();
115
116 /**
117 * Constructor using cCong = -18km/h, cFree = 80km/h, deltaV = 10km/h and vc = 80km/h. A default kernel is set.
118 */
119 public Egtf()
120 {
121 this(-18.0, 80.0, 10.0, 80.0);
122 }
123
124 /**
125 * Constructor defining global settings. A default kernel is set.
126 * @param cCong shock wave speed in congestion [km/h]
127 * @param cFree shock wave speed in free flow [km/h]
128 * @param deltaV speed range between congestion and free flow [km/h]
129 * @param vc flip-over speed below which we have congestion [km/h]
130 */
131 public Egtf(final double cCong, final double cFree, final double deltaV, final double vc)
132 {
133 this.cCong = cCong / 3.6;
134 this.cFree = cFree / 3.6;
135 this.deltaV = deltaV / 3.6;
136 this.vc = vc / 3.6;
137 setKernel();
138 }
139
140 /**
141 * Convenience constructor that also sets a specified kernel.
142 * @param cCong shock wave speed in congestion [km/h]
143 * @param cFree shock wave speed in free flow [km/h]
144 * @param deltaV speed range between congestion and free flow [km/h]
145 * @param vc flip-over speed below which we have congestion [km/h]
146 * @param sigma spatial kernel size in [m]
147 * @param tau temporal kernel size in [s]
148 * @param xMax maximum spatial range in [m]
149 * @param tMax maximum temporal range in [s]
150 */
151 @SuppressWarnings("parameternumber")
152 public Egtf(final double cCong, final double cFree, final double deltaV, final double vc, final double sigma,
153 final double tau, final double xMax, final double tMax)
154 {
155 this(cCong, cFree, deltaV, vc);
156 setKernelSI(sigma, tau, xMax, tMax);
157 }
158
159 // ********************
160 // *** DATA METHODS ***
161 // ********************
162
163 /**
164 * Return a data source, which is created if necessary.
165 * @param name unique name for the data source
166 * @return data source
167 * @throws IllegalStateException when data has been added without a data source
168 */
169 public DataSource getDataSource(final String name)
170 {
171 if (this.defaultDataSource != null)
172 {
173 throw new IllegalStateException(
174 "Obtaining a (new) data source after data has been added without a data source is not allowed.");
175 }
176 return this.dataSources.computeIfAbsent(name, (key) -> new DataSource(key));
177 }
178
179 /**
180 * Removes all data from before the given time. This is useful in live usages of this class, where older data is no longer
181 * required.
182 * @param time time before which all data can be removed
183 */
184 public synchronized void clearDataBefore(final double time)
185 {
186 for (SortedMap<Double, Map<DataStream<?>, Double>> map : this.data.values())
187 {
188 map.subMap(Double.NEGATIVE_INFINITY, time).clear();
189 }
190 }
191
192 /**
193 * Adds point data.
194 * @param quantity quantity of the data
195 * @param location location in [m]
196 * @param time time in [s]
197 * @param value data value
198 * @throws IllegalStateException if data was added with a data stream previously
199 */
200 public synchronized void addPointDataSI(final Quantity<?, ?> quantity, final double location, final double time,
201 final double value)
202 {
203 this.addingByQuantity = true;
204 addPointDataSI(getDefaultDataStream(quantity), location, time, value);
205 this.addingByQuantity = false;
206 }
207
208 /**
209 * Adds point data.
210 * @param dataStream data stream of the data
211 * @param location location in [m]
212 * @param time time in [s]
213 * @param value data value
214 * @throws IllegalStateException if data was added with a quantity previously
215 */
216 public synchronized void addPointDataSI(final DataStream<?> dataStream, final double location, final double time,
217 final double value)
218 {
219 checkNoQuantityData();
220 Objects.requireNonNull(dataStream, "Datastream may not be null.");
221 if (!Double.isNaN(value))
222 {
223 getSpacioTemporalData(getSpatialData(location), time).put(dataStream, value);
224 }
225 }
226
227 /**
228 * Adds vector data.
229 * @param quantity quantity of the data
230 * @param location locations in [m]
231 * @param time times in [s]
232 * @param values data values in SI unit
233 * @throws IllegalStateException if data was added with a data stream previously
234 */
235 public synchronized void addVectorDataSI(final Quantity<?, ?> quantity, final double[] location, final double[] time,
236 final double[] values)
237 {
238 this.addingByQuantity = true;
239 addVectorDataSI(getDefaultDataStream(quantity), location, time, values);
240 this.addingByQuantity = false;
241 }
242
243 /**
244 * Adds vector data.
245 * @param dataStream data stream of the data
246 * @param location locations in [m]
247 * @param time times in [s]
248 * @param values data values in SI unit
249 * @throws IllegalStateException if data was added with a quantity previously
250 */
251 public synchronized void addVectorDataSI(final DataStream<?> dataStream, final double[] location, final double[] time,
252 final double[] values)
253 {
254 checkNoQuantityData();
255 Objects.requireNonNull(dataStream, "Datastream may not be null.");
256 Objects.requireNonNull(location, "Location may not be null.");
257 Objects.requireNonNull(time, "Time may not be null.");
258 Objects.requireNonNull(values, "Values may not be null.");
259 if (location.length != time.length || time.length != values.length)
260 {
261 throw new IllegalArgumentException(String.format("Unequal lengths: location %d, time %d, data %d.", location.length,
262 time.length, values.length));
263 }
264 for (int i = 0; i < values.length; i++)
265 {
266 if (!Double.isNaN(values[i]))
267 {
268 getSpacioTemporalData(getSpatialData(location[i]), time[i]).put(dataStream, values[i]);
269 }
270 }
271 }
272
273 /**
274 * Adds grid data.
275 * @param quantity quantity of the data
276 * @param location locations in [m]
277 * @param time times in [s]
278 * @param values data values in SI unit
279 * @throws IllegalStateException if data was added with a data stream previously
280 */
281 public synchronized void addGridDataSI(final Quantity<?, ?> quantity, final double[] location, final double[] time,
282 final double[][] values)
283 {
284 this.addingByQuantity = true;
285 addGridDataSI(getDefaultDataStream(quantity), location, time, values);
286 this.addingByQuantity = false;
287 }
288
289 /**
290 * Adds grid data.
291 * @param dataStream data stream of the data
292 * @param location locations in [m]
293 * @param time times in [s]
294 * @param values data values in SI unit
295 * @throws IllegalStateException if data was added with a quantity previously
296 */
297 public synchronized void addGridDataSI(final DataStream<?> dataStream, final double[] location, final double[] time,
298 final double[][] values)
299 {
300 checkNoQuantityData();
301 Objects.requireNonNull(dataStream, "Datastream may not be null.");
302 Objects.requireNonNull(location, "Location may not be null.");
303 Objects.requireNonNull(time, "Time may not be null.");
304 Objects.requireNonNull(values, "Values may not be null.");
305 if (values.length != location.length)
306 {
307 throw new IllegalArgumentException(
308 String.format("%d locations while length of data is %d", location.length, values.length));
309 }
310 for (int i = 0; i < location.length; i++)
311 {
312 if (values[i].length != time.length)
313 {
314 throw new IllegalArgumentException(
315 String.format("%d times while length of data is %d", time.length, values[i].length));
316 }
317 Map<Double, Map<DataStream<?>, Double>> spatialData = getSpatialData(location[i]);
318 for (int j = 0; j < time.length; j++)
319 {
320 if (!Double.isNaN(values[i][j]))
321 {
322 getSpacioTemporalData(spatialData, time[j]).put(dataStream, values[i][j]);
323 }
324 }
325 }
326 }
327
328 /**
329 * Check that no data was added using a quantity.
330 * @throws IllegalStateException if data was added with a quantity previously
331 */
332 private void checkNoQuantityData()
333 {
334 if (!this.addingByQuantity && this.defaultDataSource != null)
335 {
336 throw new IllegalStateException(
337 "Adding data with a data stream is not allowed after data has been added with a quantity.");
338 }
339 }
340
341 /**
342 * Returns a default data stream and checks that no data with a data stream was added.
343 * @param quantity quantity
344 * @return default data stream
345 * @throws IllegalStateException if data was added with a data stream previously
346 */
347 private DataStream<?> getDefaultDataStream(final Quantity<?, ?> quantity)
348 {
349 Objects.requireNonNull(quantity, "Quantity may not be null.");
350 if (!this.dataSources.isEmpty())
351 {
352 throw new IllegalStateException(
353 "Adding data with a quantity is not allowed after data has been added with a data stream.");
354 }
355 if (this.defaultDataSource == null)
356 {
357 this.defaultDataSource = new DataSource("default");
358 this.defaultDataStreams = new LinkedHashMap<>();
359 }
360 return this.defaultDataStreams.computeIfAbsent(quantity,
361 (key) -> this.defaultDataSource.addStreamSI(quantity, 1.0, 1.0));
362 }
363
364 /**
365 * Returns data from a specific location as a subset from all data. An empty map is returned if no such data exists.
366 * @param location location in [m]
367 * @return data from a specific location
368 */
369 private SortedMap<Double, Map<DataStream<?>, Double>> getSpatialData(final double location)
370 {
371 return this.data.computeIfAbsent(location, (key) -> new TreeMap<>());
372 }
373
374 /**
375 * Returns data from a specific time as a subset of data from a specific location. An empty map is returned if no such data
376 * exists.
377 * @param spatialData spatially selected data
378 * @param time time in [s]
379 * @return data from a specific time, from data from a specific location
380 */
381 private Map<DataStream<?>, Double> getSpacioTemporalData(final Map<Double, Map<DataStream<?>, Double>> spatialData,
382 final double time)
383 {
384 return spatialData.computeIfAbsent(time, (key) -> new LinkedHashMap<>());
385 }
386
387 // **********************
388 // *** KERNEL METHODS ***
389 // **********************
390
391 /**
392 * Sets a default exponential kernel with infinite range, sigma = 300m, and tau = 30s.
393 */
394 public void setKernel()
395 {
396 setKernelSI(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, new ExpKernelShape(DEFAULT_SIGMA, DEFAULT_TAU));
397 }
398
399 /**
400 * Sets an exponential kernel with infinite range.
401 * @param sigma spatial kernel size
402 * @param tau temporal kernel size
403 */
404 public void setKernelSI(final double sigma, final double tau)
405 {
406 setKernelSI(sigma, tau, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
407 }
408
409 /**
410 * Sets an exponential kernel with limited range.
411 * @param sigma spatial kernel size in [m]
412 * @param tau temporal kernel size in [s]
413 * @param xMax maximum spatial range in [m]
414 * @param tMax maximum temporal range in [s]
415 */
416 public void setKernelSI(final double sigma, final double tau, final double xMax, final double tMax)
417 {
418 setKernelSI(xMax, tMax, new ExpKernelShape(sigma, tau));
419 }
420
421 /**
422 * Sets a Gaussian kernel with infinite range.
423 * @param sigma spatial kernel size
424 * @param tau temporal kernel size
425 */
426 public void setGaussKernelSI(final double sigma, final double tau)
427 {
428 setGaussKernelSI(sigma, tau, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
429 }
430
431 /**
432 * Sets a Gaussian kernel with limited range.
433 * @param sigma spatial kernel size in [m]
434 * @param tau temporal kernel size in [s]
435 * @param xMax maximum spatial range in [m]
436 * @param tMax maximum temporal range in [s]
437 */
438 public void setGaussKernelSI(final double sigma, final double tau, final double xMax, final double tMax)
439 {
440 setKernelSI(xMax, tMax, new GaussKernelShape(sigma, tau));
441 }
442
443 /**
444 * Sets a kernel with limited range and provided shape. The shape allows using non-exponential kernels.
445 * @param xMax maximum spatial range
446 * @param tMax maximum temporal range
447 * @param shape shape of the kernel
448 */
449 public synchronized void setKernelSI(final double xMax, final double tMax, final KernelShape shape)
450 {
451 this.kernel = new Kernel(xMax, tMax, shape);
452 }
453
454 /**
455 * Returns the wave speed in congestion.
456 * @return wave speed in congestion
457 */
458 final double getWaveSpeedCongestion()
459 {
460 return this.cCong;
461 }
462
463 /**
464 * Returns the wave speed in free flow.
465 * @return wave speed in free flow
466 */
467 final double getWaveSpeedFreeFlow()
468 {
469 return this.cFree;
470 }
471
472 // **********************
473 // *** FILTER METHODS ***
474 // **********************
475
476 /**
477 * Executes filtering in parallel. The returned listener can be used to report progress and wait until the filtering is
478 * done. Finally, the filtering results can then be obtained from the listener.
479 * @param location location of output grid in [m]
480 * @param time time of output grid in [s]
481 * @param quantities quantities to calculate filtered data of
482 * @return listener to notify keep track of the progress
483 */
484 public EgtfParallelListener filterParallelSI(final double[] location, final double[] time,
485 final Quantity<?, ?>... quantities)
486 {
487 Objects.requireNonNull(location, "Location may not be null.");
488 Objects.requireNonNull(time, "Time may not be null.");
489 EgtfParallelListener listener = new EgtfParallelListener();
490 addListener(listener);
491 new Thread(new Runnable()
492 {
493 @Override
494 public void run()
495 {
496 Optional<Filter> filter = filterSI(location, time, quantities);
497 if (filter.isPresent())
498 {
499 listener.setFilter(filter.get());
500 }
501 removeListener(listener);
502 }
503 }, "Egtf calculation thread").start();
504 return listener;
505 }
506
507 /**
508 * Executes fast filtering in parallel. The returned listener can be used to report progress and wait until the filtering is
509 * done. Finally, the filtering results can then be obtained from the listener.
510 * @param xMin minimum location value of output grid [m]
511 * @param xStep location step of output grid [m]
512 * @param xMax maximum location value of output grid [m]
513 * @param tMin minimum time value of output grid [s]
514 * @param tStep time step of output grid [s]
515 * @param tMax maximum time value of output grid [s]
516 * @param quantities quantities to calculate filtered data of
517 * @return listener to notify keep track of the progress
518 */
519 public EgtfParallelListener filterParallelFastSI(final double xMin, final double xStep, final double xMax,
520 final double tMin, final double tStep, final double tMax, final Quantity<?, ?>... quantities)
521 {
522 EgtfParallelListener listener = new EgtfParallelListener();
523 addListener(listener);
524 new Thread(new Runnable()
525 {
526 @Override
527 public void run()
528 {
529 Optional<Filter> filter = filterFastSI(xMin, xStep, xMax, tMin, tStep, tMax, quantities);
530 if (filter.isPresent())
531 {
532 listener.setFilter(filter.get());
533 }
534 removeListener(listener);
535 }
536 }, "Egtf calculation thread").start();
537 return listener;
538 }
539
540 /**
541 * Returns filtered data. This is the standard EGTF implementation.
542 * @param location location of output grid in [m]
543 * @param time time of output grid in [s]
544 * @param quantities quantities to calculate filtered data of
545 * @return filtered data, empty when interrupted
546 */
547 @SuppressWarnings("methodlength")
548 public Optional<Filter> filterSI(final double[] location, final double[] time, final Quantity<?, ?>... quantities)
549 {
550 Objects.requireNonNull(location, "Location may not be null.");
551 Objects.requireNonNull(time, "Time may not be null.");
552
553 // initialize data
554 Map<Quantity<?, ?>, double[][]> map = new LinkedHashMap<>();
555 for (Quantity<?, ?> quantity : quantities)
556 {
557 map.put(quantity, new double[location.length][time.length]);
558 }
559
560 // loop grid locations
561 for (int i = 0; i < location.length; i++)
562 {
563 double xGrid = location[i];
564
565 // filter applicable data for location
566 Map<Double, NavigableMap<Double, Map<DataStream<?>, Double>>> spatialData =
567 this.data.subMap(this.kernel.fromLocation(xGrid), true, this.kernel.toLocation(xGrid), true);
568
569 // loop grid times
570 for (int j = 0; j < time.length; j++)
571 {
572 double tGrid = time[j];
573
574 // notify
575 if (notifyListeners((i + (double) j / time.length) / location.length))
576 {
577 return Optional.empty();
578 }
579
580 // initialize data per stream
581 // quantity z assuming congestion and free flow
582 Map<DataStream<?>, DualWeightedMean> zCongFree = new LinkedHashMap<>();
583
584 // filter and loop applicable data for time
585 for (Map.Entry<Double, NavigableMap<Double, Map<DataStream<?>, Double>>> xEntry : spatialData.entrySet())
586 {
587 double dx = xEntry.getKey() - xGrid;
588 Map<Double, Map<DataStream<?>, Double>> temporalData =
589 xEntry.getValue().subMap(this.kernel.fromTime(tGrid), true, this.kernel.toTime(tGrid), true);
590
591 for (Map.Entry<Double, Map<DataStream<?>, Double>> tEntry : temporalData.entrySet())
592 {
593 double dt = tEntry.getKey() - tGrid;
594 Map<DataStream<?>, Double> pData = tEntry.getValue();
595
596 double phiCong = this.kernel.weight(this.cCong, dx, dt);
597 double phiFree = this.kernel.weight(this.cFree, dx, dt);
598
599 // loop streams data at point
600 for (Map.Entry<DataStream<?>, Double> vEntry : pData.entrySet())
601 {
602 DataStream<?> stream = vEntry.getKey();
603 if (map.containsKey(stream.getQuantity()) || stream.getQuantity().isSpeed())
604 {
605 double v = vEntry.getValue();
606 DualWeightedMean zCongFreeOfStream =
607 zCongFree.computeIfAbsent(stream, (key) -> new DualWeightedMean());
608 zCongFreeOfStream.addCong(v, phiCong);
609 zCongFreeOfStream.addFree(v, phiFree);
610 }
611 }
612 }
613 }
614
615 // figure out the congestion level estimated for each data source
616 Map<DataSource, Double> w = new LinkedHashMap<>();
617 for (Map.Entry<DataStream<?>, DualWeightedMean> streamEntry : zCongFree.entrySet())
618 {
619 DataStream<?> dataStream = streamEntry.getKey();
620 if (dataStream.getQuantity().isSpeed()) // only one speed quantity allowed per data source
621 {
622 DualWeightedMean zCongFreeOfStream = streamEntry.getValue();
623 double u = Math.min(zCongFreeOfStream.getCong(), zCongFreeOfStream.getFree());
624 w.put(dataStream.getDataSource(), // 1 speed quantity per source allowed
625 .5 * (1.0 + Math.tanh((Egtf.this.vc - u) / Egtf.this.deltaV)));
626 continue;
627 }
628 }
629
630 // sum available data sources per quantity
631 Double wMean = null;
632 for (Map.Entry<Quantity<?, ?>, double[][]> qEntry : map.entrySet())
633 {
634 Quantity<?, ?> quantity = qEntry.getKey();
635 WeightedMean z = new WeightedMean();
636 for (Map.Entry<DataStream<?>, DualWeightedMean> zEntry : zCongFree.entrySet())
637 {
638 DataStream<?> dataStream = zEntry.getKey();
639 if (dataStream.getQuantity().equals(quantity))
640 {
641 // obtain congestion level
642 double wCong;
643 if (!w.containsKey(dataStream.getDataSource()))
644 {
645 // this data source has no speed data, but congestion level can be estimated from other sources
646 if (wMean == null)
647 {
648 // let's see if speed was estimated already
649 for (Quantity<?, ?> prevQuant : quantities)
650 {
651 if (prevQuant.equals(quantity))
652 {
653 // it was not, get mean of other data source
654 wMean = 0.0;
655 for (double ww : w.values())
656 {
657 wMean += ww / w.size();
658 }
659 break;
660 }
661 else if (prevQuant.isSpeed())
662 {
663 wMean = .5 * (1.0
664 + Math.tanh((Egtf.this.vc - map.get(prevQuant)[i][j]) / Egtf.this.deltaV));
665 break;
666 }
667 }
668 }
669 wCong = wMean;
670 }
671 else
672 {
673 wCong = w.get(dataStream.getDataSource());
674 }
675 // calculate estimated value z of this data source (no duplicate quantities per source allowed)
676 double wfree = 1.0 - wCong;
677 DualWeightedMean zCongFreej = zEntry.getValue();
678 double zStream = wCong * zCongFreej.getCong() + wfree * zCongFreej.getFree();
679 double weight;
680 if (w.size() > 1)
681 {
682 // data source more important if more and nearer measurements
683 double beta = wCong * zCongFreej.getDenominatorCong() + wfree * zCongFreej.getDenominatorFree();
684 // more important if more reliable (smaller standard deviation) at congestion level
685 double alpha = wCong / dataStream.getThetaCong() + wfree / dataStream.getThetaFree();
686 weight = alpha * beta;
687 }
688 else
689 {
690 weight = 1.0;
691 }
692 z.add(zStream, weight);
693 }
694 }
695 qEntry.getValue()[i][j] = z.get();
696 }
697 }
698 }
699 notifyListeners(1.0);
700
701 return Optional.of(new FilterDouble(location, time, map));
702 }
703
704 /**
705 * Returns filtered data that is processed using fast fourier transformation. This is much faster than the standard filter,
706 * at the cost that all input data is discretized to the output grid. The gain in computation speed is however such that
707 * finer output grids can be used to alleviate this. For discretization the output grid needs to be equidistant. It is
708 * recommended to set a Kernel with maximum bounds before using this method.
709 * <p>
710 * More than being a fast implementation of the Adaptive Smoothing Method, this implementation includes all data source like
711 * the Extended Generalized Treiber-Helbing Filter.
712 * @param xMin minimum location value of output grid [m]
713 * @param xStep location step of output grid [m]
714 * @param xMax maximum location value of output grid [m]
715 * @param tMin minimum time value of output grid [s]
716 * @param tStep time step of output grid [s]
717 * @param tMax maximum time value of output grid [s]
718 * @param quantities quantities to calculate filtered data of
719 * @return filtered data, empty when interrupted
720 */
721 @SuppressWarnings("methodlength")
722 public Optional<Filter> filterFastSI(final double xMin, final double xStep, final double xMax, final double tMin,
723 final double tStep, final double tMax, final Quantity<?, ?>... quantities)
724 {
725 if (xMin > xMax || xStep <= 0.0 || tMin > tMax || tStep <= 0.0)
726 {
727 throw new IllegalArgumentException(
728 "Ill-defined grid. Make sure that xMax >= xMin, dx > 0, tMax >= tMin and dt > 0");
729 }
730 if (notifyListeners(0.0))
731 {
732 return Optional.empty();
733 }
734
735 // initialize data
736 int n = 1 + (int) ((xMax - xMin) / xStep);
737 double[] location = new double[n];
738 IntStream.range(0, n).forEach(i -> location[i] = xMin + i * xStep);
739 n = 1 + (int) ((tMax - tMin) / tStep);
740 double[] time = new double[n];
741 IntStream.range(0, n).forEach(j -> time[j] = tMin + j * tStep);
742 Map<Quantity<?, ?>, double[][]> map = new LinkedHashMap<>();
743 Map<Quantity<?, ?>, double[][]> weights = new LinkedHashMap<>();
744 for (Quantity<?, ?> quantity : quantities)
745 {
746 map.put(quantity, new double[location.length][time.length]);
747 weights.put(quantity, new double[location.length][time.length]);
748 }
749
750 // discretize Kernel
751 double xFrom = this.kernel.fromLocation(0.0);
752 xFrom = Double.isInfinite(xFrom) ? 2.0 * (xMin - xMax) : xFrom;
753 double xTo = this.kernel.toLocation(0.0);
754 xTo = Double.isInfinite(xTo) ? 2.0 * (xMax - xMin) : xTo;
755 double[] dx = equidistant(xFrom, xStep, xTo);
756 double tFrom = this.kernel.fromTime(0.0);
757 tFrom = Double.isInfinite(tFrom) ? 2.0 * (tMin - tMax) : tFrom;
758 double tTo = this.kernel.toTime(0.0);
759 tTo = Double.isInfinite(tTo) ? 2.0 * (tMax - tMin) : tTo;
760 double[] dt = equidistant(tFrom, tStep, tTo);
761 double[][] phiCong = new double[dx.length][dt.length];
762 double[][] phiFree = new double[dx.length][dt.length];
763 for (int i = 0; i < dx.length; i++)
764 {
765 for (int j = 0; j < dt.length; j++)
766 {
767 phiCong[i][j] = this.kernel.weight(this.cCong, dx[i], dt[j]);
768 phiFree[i][j] = this.kernel.weight(this.cFree, dx[i], dt[j]);
769 }
770 }
771
772 // discretize data
773 Map<DataStream<?>, double[][]> dataSum = new LinkedHashMap<>();
774 Map<DataStream<?>, double[][]> dataCount = new LinkedHashMap<>(); // integer counts, must be double[][] for convolution
775 // loop grid locations
776 for (int i = 0; i < location.length; i++)
777 {
778 // filter applicable data for location
779 Map<Double, NavigableMap<Double, Map<DataStream<?>, Double>>> spatialData =
780 this.data.subMap(location[i] - 0.5 * xStep, true, location[i] + 0.5 * xStep, true);
781 // loop grid times
782 for (int j = 0; j < time.length; j++)
783 {
784 // filter and loop applicable data for time
785 for (NavigableMap<Double, Map<DataStream<?>, Double>> locationData : spatialData.values())
786 {
787 NavigableMap<Double, Map<DataStream<?>, Double>> temporalData =
788 locationData.subMap(time[j] - 0.5 * tStep, true, time[j] + 0.5 * tStep, true);
789 for (Map<DataStream<?>, Double> timeData : temporalData.values())
790 {
791 for (Map.Entry<DataStream<?>, Double> timeEntry : timeData.entrySet())
792 {
793 if (map.containsKey(timeEntry.getKey().getQuantity()) || timeEntry.getKey().getQuantity().isSpeed())
794 {
795 dataSum.computeIfAbsent(timeEntry.getKey(),
796 (key) -> new double[location.length][time.length])[i][j] += timeEntry.getValue();
797 dataCount.computeIfAbsent(timeEntry.getKey(),
798 (key) -> new double[location.length][time.length])[i][j]++;
799 }
800 }
801 }
802 }
803 }
804 }
805
806 // figure out the congestion level estimated for each data source
807 double steps = quantities.length + 1; // speed (for congestion level) and then all in quantities
808 double step = 0;
809 // store maps to prevent us from calculating the convolution for speed again later
810 Map<DataSource, double[][]> w = new LinkedHashMap<>();
811 Map<DataSource, double[][]> zCongSpeed = new LinkedHashMap<>();
812 Map<DataSource, double[][]> zFreeSpeed = new LinkedHashMap<>();
813 Map<DataSource, double[][]> nCongSpeed = new LinkedHashMap<>();
814 Map<DataSource, double[][]> nFreeSpeed = new LinkedHashMap<>();
815 for (Map.Entry<DataStream<?>, double[][]> zEntry : dataSum.entrySet())
816 {
817 DataStream<?> dataStream = zEntry.getKey();
818 if (dataStream.getQuantity().isSpeed()) // only one speed quantity allowed per data source
819 {
820 // notify
821 double[][] vCong = Convolution.convolution(phiCong, zEntry.getValue());
822 if (notifyListeners((step + 0.25) / steps))
823 {
824 return null;
825 }
826 double[][] vFree = Convolution.convolution(phiFree, zEntry.getValue());
827 if (notifyListeners((step + 0.5) / steps))
828 {
829 return null;
830 }
831 double[][] count = dataCount.get(dataStream);
832 double[][] nCong = Convolution.convolution(phiCong, count);
833 if (notifyListeners((step + 0.75) / steps))
834 {
835 return null;
836 }
837 double[][] nFree = Convolution.convolution(phiFree, count);
838 double[][] wSource = new double[vCong.length][vCong[0].length];
839 for (int i = 0; i < vCong.length; i++)
840 {
841 for (int j = 0; j < vCong[0].length; j++)
842 {
843 double u = Math.min(vCong[i][j] / nCong[i][j], vFree[i][j] / nFree[i][j]);
844 wSource[i][j] = .5 * (1.0 + Math.tanh((Egtf.this.vc - u) / Egtf.this.deltaV));
845 }
846 }
847 w.put(dataStream.getDataSource(), wSource);
848 zCongSpeed.put(dataStream.getDataSource(), vCong);
849 zFreeSpeed.put(dataStream.getDataSource(), vFree);
850 nCongSpeed.put(dataStream.getDataSource(), nCong);
851 nFreeSpeed.put(dataStream.getDataSource(), nFree);
852 }
853 }
854 step++;
855 if (notifyListeners(step / steps))
856 {
857 return null;
858 }
859
860 // sum available data sources per quantity
861 double[][] wMean = null;
862 for (Quantity<?, ?> quantity : quantities)
863 {
864 // gather place for this quantity
865 double[][] qData = map.get(quantity);
866 double[][] qWeights = weights.get(quantity);
867 // loop streams that provide this quantity
868 Set<Map.Entry<DataStream<?>, double[][]>> zEntries = new LinkedHashSet<>();
869 for (Map.Entry<DataStream<?>, double[][]> zEntry : dataSum.entrySet())
870 {
871 if (zEntry.getKey().getQuantity().equals(quantity))
872 {
873 zEntries.add(zEntry);
874 }
875 }
876 double streamCounter = 0;
877 for (Map.Entry<DataStream<?>, double[][]> zEntry : zEntries)
878 {
879 DataStream<?> dataStream = zEntry.getKey();
880
881 // obtain congestion level
882 double[][] wj;
883 if (!w.containsKey(dataStream.getDataSource()))
884 {
885 // this data source has no speed data, but congestion level can be estimated from other sources
886 if (wMean == null)
887 {
888 // let's see if speed was estimated already
889 for (Quantity<?, ?> prevQuant : quantities)
890 {
891 if (prevQuant.equals(quantity))
892 {
893 // it was not, get mean of other data source
894 wMean = new double[location.length][time.length];
895 for (double[][] ww : w.values())
896 {
897 for (int i = 0; i < location.length; i++)
898 {
899 for (int j = 0; j < time.length; j++)
900 {
901 wMean[i][j] += ww[i][j] / w.size();
902 }
903 }
904 }
905 break;
906 }
907 else if (prevQuant.isSpeed())
908 {
909 wMean = new double[location.length][time.length];
910 double[][] v = map.get(prevQuant);
911 for (int i = 0; i < location.length; i++)
912 {
913 for (int j = 0; j < time.length; j++)
914 {
915 wMean[i][j] = .5 * (1.0 + Math.tanh((Egtf.this.vc - v[i][j]) / Egtf.this.deltaV));
916 }
917 }
918 break;
919 }
920 }
921 }
922 wj = wMean;
923 }
924 else
925 {
926 wj = w.get(dataStream.getDataSource());
927 }
928
929 // convolutions of filters with discretized data and data counts
930 double[][] zCong;
931 double[][] zFree;
932 double[][] nCong;
933 double[][] nFree;
934 if (dataStream.getQuantity().isSpeed())
935 {
936 zCong = zCongSpeed.get(dataStream.getDataSource());
937 zFree = zFreeSpeed.get(dataStream.getDataSource());
938 nCong = nCongSpeed.get(dataStream.getDataSource());
939 nFree = nFreeSpeed.get(dataStream.getDataSource());
940 }
941 else
942 {
943 zCong = Convolution.convolution(phiCong, zEntry.getValue());
944 if (notifyListeners((step + (streamCounter + 0.25) / zEntries.size()) / steps))
945 {
946 return null;
947 }
948 zFree = Convolution.convolution(phiFree, zEntry.getValue());
949 if (notifyListeners((step + (streamCounter + 0.5) / zEntries.size()) / steps))
950 {
951 return null;
952 }
953 double[][] count = dataCount.get(dataStream);
954 nCong = Convolution.convolution(phiCong, count);
955 if (notifyListeners((step + (streamCounter + 0.75) / zEntries.size()) / steps))
956 {
957 return null;
958 }
959 nFree = Convolution.convolution(phiFree, count);
960 }
961
962 // loop grid to add to each weighted sum (weighted per data source)
963 for (int i = 0; i < location.length; i++)
964 {
965 for (int j = 0; j < time.length; j++)
966 {
967 double wCong = wj[i][j];
968 double wFree = 1.0 - wCong;
969 double value = wCong * zCong[i][j] / nCong[i][j] + wFree * zFree[i][j] / nFree[i][j];
970 // the fast filter supplies convoluted data counts, i.e. amount of data and filter proximity; this
971 // is exactly what the EGTF method needs to weigh data sources
972 double beta = wCong * nCong[i][j] + wFree * nFree[i][j];
973 double alpha = wCong / dataStream.getThetaCong() + wFree / dataStream.getThetaFree();
974 double weight = beta * alpha;
975 qData[i][j] += (value * weight);
976 qWeights[i][j] += weight;
977 }
978 }
979 streamCounter++;
980 if (notifyListeners((step + streamCounter / zEntries.size()) / steps))
981 {
982 return null;
983 }
984 }
985 for (int i = 0; i < location.length; i++)
986 {
987 for (int j = 0; j < time.length; j++)
988 {
989 qData[i][j] /= qWeights[i][j];
990 }
991 }
992 step++;
993 }
994
995 return Optional.of(new FilterDouble(location, time, map));
996 }
997
998 /**
999 * Returns an equidistant vector that includes 0.
1000 * @param from lowest value to include
1001 * @param step step
1002 * @param to highest value to include
1003 * @return equidistant vector that includes 0
1004 */
1005 private double[] equidistant(final double from, final double step, final double to)
1006 {
1007 int n1 = (int) (-from / step);
1008 int n2 = (int) (to / step);
1009 int n = n1 + n2 + 1;
1010 double[] array = new double[n];
1011 for (int i = 0; i < n; i++)
1012 {
1013 array[i] = i < n1 ? step * (-n1 + i) : step * (i - n1);
1014 }
1015 return array;
1016 }
1017
1018 // *********************
1019 // *** EVENT METHODS ***
1020 // *********************
1021
1022 /**
1023 * Interrupt the calculation.
1024 */
1025 public final void interrupt()
1026 {
1027 this.interrupted = true;
1028 }
1029
1030 /**
1031 * Add listener.
1032 * @param listener listener
1033 */
1034 public final void addListener(final EgtfListener listener)
1035 {
1036 this.listeners.add(listener);
1037 }
1038
1039 /**
1040 * Remove listener.
1041 * @param listener listener
1042 */
1043 public final void removeListener(final EgtfListener listener)
1044 {
1045 this.listeners.remove(listener);
1046 }
1047
1048 /**
1049 * Notify all listeners.
1050 * @param progress progress, a value in the range [0 ... 1]
1051 * @return whether the filter is interrupted
1052 */
1053 private boolean notifyListeners(final double progress)
1054 {
1055 if (!this.listeners.isEmpty())
1056 {
1057 EgtfEvent event = new EgtfEvent(this, progress);
1058 for (EgtfListener listener : this.listeners)
1059 {
1060 listener.notifyProgress(event);
1061 }
1062 }
1063 return this.interrupted;
1064 }
1065
1066 // **********************
1067 // *** HELPER CLASSES ***
1068 // **********************
1069
1070 /**
1071 * Small class to build up a weighted mean under the congestion and free flow assumption.
1072 */
1073 private final class DualWeightedMean
1074 {
1075 /** Cumulative congestion numerator of weighted mean fraction, i.e. weighted sum. */
1076 private double numeratorCong;
1077
1078 /** Cumulative free flow numerator of weighted mean fraction, i.e. weighted sum. */
1079 private double numeratorFree;
1080
1081 /** Cumulative congestion denominator of weighted mean fraction, i.e. sum of weights. */
1082 private double denominatorCong;
1083
1084 /** Cumulative free flow denominator of weighted mean fraction, i.e. sum of weights. */
1085 private double denominatorFree;
1086
1087 /**
1088 * Adds a congestion value with weight.
1089 * @param value value
1090 * @param weight weight
1091 */
1092 public void addCong(final double value, final double weight)
1093 {
1094 this.numeratorCong += value * weight;
1095 this.denominatorCong += weight;
1096 }
1097
1098 /**
1099 * Adds a free flow value with weight.
1100 * @param value value
1101 * @param weight weight
1102 */
1103 public void addFree(final double value, final double weight)
1104 {
1105 this.numeratorFree += value * weight;
1106 this.denominatorFree += weight;
1107 }
1108
1109 /**
1110 * Returns the weighted congestion mean of available data.
1111 * @return weighted mean of available data
1112 */
1113 public double getCong()
1114 {
1115 return this.numeratorCong / this.denominatorCong;
1116 }
1117
1118 /**
1119 * Returns the weighted free flow mean of available data.
1120 * @return weighted free flow mean of available data
1121 */
1122 public double getFree()
1123 {
1124 return this.numeratorFree / this.denominatorFree;
1125 }
1126
1127 /**
1128 * Returns the sum of congestion weights.
1129 * @return the sum of congestion weights
1130 */
1131 public double getDenominatorCong()
1132 {
1133 return this.denominatorCong;
1134 }
1135
1136 /**
1137 * Returns the sum of free flow weights.
1138 * @return the sum of free flow weights
1139 */
1140 public double getDenominatorFree()
1141 {
1142 return this.denominatorFree;
1143 }
1144
1145 @Override
1146 public String toString()
1147 {
1148 return "DualWeightedMean [numeratorCong=" + this.numeratorCong + ", numeratorFree=" + this.numeratorFree
1149 + ", denominatorCong=" + this.denominatorCong + ", denominatorFree=" + this.denominatorFree + "]";
1150 }
1151
1152 }
1153
1154 /**
1155 * Small class to build up a weighted mean.
1156 */
1157 private final class WeightedMean
1158 {
1159 /** Cumulative numerator of weighted mean fraction, i.e. weighted sum. */
1160 private double numerator;
1161
1162 /** Cumulative denominator of weighted mean fraction, i.e. sum of weights. */
1163 private double denominator;
1164
1165 /**
1166 * Adds a value with weight.
1167 * @param value value
1168 * @param weight weight
1169 */
1170 public void add(final double value, final double weight)
1171 {
1172 this.numerator += value * weight;
1173 this.denominator += weight;
1174 }
1175
1176 /**
1177 * Returns the weighted mean of available data.
1178 * @return weighted mean of available data
1179 */
1180 public double get()
1181 {
1182 return this.numerator / this.denominator;
1183 }
1184
1185 @Override
1186 public String toString()
1187 {
1188 return "WeightedMean [numerator=" + this.numerator + ", denominator=" + this.denominator + "]";
1189 }
1190
1191 }
1192
1193 @Override
1194 public String toString()
1195 {
1196 return "EGTF [kernel=" + this.kernel + ", cCong=" + this.cCong + ", cFree=" + this.cFree + ", deltaV=" + this.deltaV
1197 + ", vc=" + this.vc + ", dataSources=" + this.dataSources + ", data=" + this.data + ", interrupted="
1198 + this.interrupted + ", listeners=" + this.listeners + "]";
1199 }
1200
1201 }