package org.opentrafficsim.animation.graphs;
   
   import java.util.ArrayList;
   import java.util.Arrays;
   import java.util.LinkedHashMap;
   import java.util.LinkedHashSet;
   import java.util.List;
   import java.util.Map;
   import java.util.Optional;
  import java.util.Set;
  
  import org.djunits.unit.SpeedUnit;
  import org.djunits.value.vdouble.scalar.Duration;
  import org.djunits.value.vdouble.scalar.Frequency;
  import org.djunits.value.vdouble.scalar.Length;
  import org.djunits.value.vdouble.scalar.LinearDensity;
  import org.djunits.value.vdouble.scalar.Speed;
  import org.djutils.event.EventType;
  import org.djutils.exceptions.Throw;
  import org.djutils.math.means.ArithmeticMean;
  import org.djutils.metadata.MetaData;
  import org.djutils.metadata.ObjectDescriptor;
  import org.opentrafficsim.animation.egtf.Converter;
  import org.opentrafficsim.animation.egtf.DataSource;
  import org.opentrafficsim.animation.egtf.DataStream;
  import org.opentrafficsim.animation.egtf.Egtf;
  import org.opentrafficsim.animation.egtf.Filter;
  import org.opentrafficsim.animation.egtf.Quantity;
  import org.opentrafficsim.animation.egtf.typed.TypedQuantity;
  import org.opentrafficsim.animation.graphs.AbstractContourPlot.ContourPaintState;
  import org.opentrafficsim.animation.graphs.GraphPath.Section;
  import org.opentrafficsim.base.logger.Logger;
  import org.opentrafficsim.kpi.interfaces.LaneData;
  import org.opentrafficsim.kpi.sampling.SamplerData;
  import org.opentrafficsim.kpi.sampling.Trajectory;
  import org.opentrafficsim.kpi.sampling.Trajectory.SpaceTimeView;
  import org.opentrafficsim.kpi.sampling.TrajectoryGroup;
  
  /**
   * Class that contains data for contour plots. One data source can be shared between contour plots, in which case the
   * granularity, path, sampler, update interval, and whether the data is smoothed (EGTF) are equal between the plots.
   * <p>
   * By default the source contains traveled time and traveled distance per cell.
   * <p>
   * Copyright (c) 2013-2026 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
   * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
   * </p>
   * @author Alexander Verbraeck
   * @author Peter Knoppers
   * @author Wouter Schakel
   */
  public class ContourDataSource extends PlotDelegate<ContourPaintState, AbstractContourPlot<?>>
  {
  
      // *******************
      // *** EVENT TYPES ***
      // *******************
  
      /** Granularity changed. */
      public static final EventType GRANULARITY = new EventType("GRANULARITY",
              new MetaData("Granularity", "Granularity changed.", new ObjectDescriptor("Axis", "Axis", Dimension.class),
                      new ObjectDescriptor("Granularity", "Granularity", Double.class)));
  
      /** Interpolation changed. */
      public static final EventType INTERPOLATE = new EventType("INTERPOLATE", new MetaData("Interpolate", "Interpolate changed.",
              new ObjectDescriptor("Interpolate", "Interpolate", Boolean.class)));
  
      /** Smooth changed. */
      public static final EventType SMOOTH = new EventType("SMOOTH",
              new MetaData("Smooth", "Smooth changed.", new ObjectDescriptor("Smooth", "Smooth", Boolean.class)));
  
      // *************************
      // *** GLOBAL PROPERTIES ***
      // *************************
  
      /** Space granularities. */
      protected static final PlotSetting<Length> SPACE_GRANULARITIES =
              PlotSetting.of(new double[] {10.0, 20.0, 50.0, 100.0, 200.0, 500.0, 1000.0}, Length::ofSI, 3);
  
      /** Time granularities. */
      protected static final PlotSetting<Duration> TIME_GRANULARITIES =
              PlotSetting.of(new double[] {1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 60.0, 120.0, 300.0, 600.0}, Duration::ofSI, 3);
  
      /** Initial lower bound for the time scale. */
      protected static final Duration DEFAULT_LOWER_TIME_BOUND = Duration.ZERO;
  
      /**
       * Total kernel size relative to sigma and tau. This factor is determined through -log(1 - p) with p ~= 99%. This means that
       * the cumulative exponential distribution has 99% at 5 times sigma or tau. Note that due to a coordinate change in the
       * Adaptive Smoothing Method, the actual cumulative distribution is slightly different. Hence, this is just a heuristic.
       */
      private static final int KERNEL_FACTOR = 5;
  
      /** Maximum free flow propagation speed. */
      private static final Speed MAX_C_FREE = new Speed(80.0, SpeedUnit.KM_PER_HOUR);
  
      /** Factor on speed limit to determine vc, the flip over speed between congestion and free flow. */
      private static final double VC_FACRTOR = 0.8;
  
     /** Congestion propagation speed. */
     private static final Speed C_CONG = new Speed(-18.0, SpeedUnit.KM_PER_HOUR);
 
     /** Delta v, speed transition region around threshold. */
     private static final Speed DELTA_V = new Speed(10.0, SpeedUnit.KM_PER_HOUR);
 
     // *****************************
     // *** CONTEXTUAL PROPERTIES ***
     // *****************************
 
     /** Sampler data. */
     private final SamplerData<?> samplerData;
 
     /** Path. */
     private final GraphPath<? extends LaneData<?>> path;
 
     /** Space axis. */
     private final Axis spaceAxis;
 
     /** Time axis. */
     private final Axis timeAxis;
 
     /** Data types. */
     private final Set<ContourDataType<?>> dataTypes = new LinkedHashSet<>();
 
     // *****************
     // *** PLOT DATA ***
     // *****************
 
     /** Total distance traveled per cell. */
     private float[][] distance;
 
     /** Total time traveled per cell. */
     private float[][] time;
 
     /** Data of other types. */
     private final Map<ContourAdditionalDataType<?, ?>, float[][]> additionalData = new LinkedHashMap<>();
 
     // ****************************
     // *** SMOOTHING PROPERTIES ***
     // ****************************
 
     /** Free flow propagation speed. */
     private Speed cFree;
 
     /** Flip-over speed between congestion and free flow. */
     private Speed vc;
 
     /** Smoothing filter. */
     private Egtf egtf;
 
     /** Data stream for speed. */
     private DataStream<Speed> speedStream;
 
     /** Data stream for travel time. */
     private DataStream<Duration> travelTimeStream;
 
     /** Data stream for travel distance. */
     private DataStream<Length> travelDistanceStream;
 
     /** Quantity for travel time. */
     private final Quantity<Duration, double[][]> travelTimeQuantity = new Quantity<>("travel time", Converter.SI);
 
     /** Quantity for travel distance. */
     private final Quantity<Length, double[][]> travelDistanceQuantity = new Quantity<>("travel distance", Converter.SI);
 
     /** Data streams for any additional data. */
     private Map<ContourAdditionalDataType<?, ?>, DataStream<?>> additionalStreams = new LinkedHashMap<>();
 
     // *****************************
     // *** CONTINUITY PROPERTIES ***
     // *****************************
 
     /** Time up to which to determine data. This is a multiple of the update interval, which is now, or recent on a redo. */
     private double toTime = 0.0;
 
     /** Number of items that are ready. To return NaN values if not ready, and for operations between consecutive updates. */
     private int readyItems = -1;
 
     /** Whether to smooth data. */
     private boolean smooth = false;
 
     // ********************
     // *** CONSTRUCTORS ***
     // ********************
 
     /**
      * Constructor using default granularities.
      * @param samplerData sampler data
      * @param path path
      * @param plotScheduler plot scheduler
      */
     public ContourDataSource(final SamplerData<?> samplerData, final GraphPath<? extends LaneData<?>> path,
             final PlotScheduler plotScheduler)
     {
         this(samplerData, Duration.ofSI(1.0), path, plotScheduler, SPACE_GRANULARITIES, TIME_GRANULARITIES,
                 DEFAULT_LOWER_TIME_BOUND, AbstractPlot.DEFAULT_INITIAL_UPPER_TIME_BOUND);
     }
 
     /**
      * Constructor for non-default input.
      * @param samplerData sampler data
      * @param delay delay so critical future events have occurred, e.g. GTU's next move's to extend trajectories
      * @param path path
      * @param plotScheduler plot scheduler
      * @param spaceGranularities granularity options for space dimension
      * @param timeGranularities granularity options for time dimension
      * @param start start time
      * @param initialEnd initial end time of plots, will be expanded if simulation time exceeds it
      */
     @SuppressWarnings("parameternumber")
     public ContourDataSource(final SamplerData<?> samplerData, final Duration delay,
             final GraphPath<? extends LaneData<?>> path, final PlotScheduler plotScheduler,
             final PlotSetting<Length> spaceGranularities, final PlotSetting<Duration> timeGranularities, final Duration start,
             final Duration initialEnd)
     {
         super(Duration.ofSI(timeGranularities.defaultValueIndex()), delay, plotScheduler);
         this.samplerData = samplerData;
         this.path = path;
         this.spaceAxis = new Axis(0.0, path.getTotalLength().si, spaceGranularities.getDefaultValue().si,
                 spaceGranularities.values().stream().mapToDouble((len) -> len.si).toArray());
         this.timeAxis = new Axis(start.si, initialEnd.si, timeGranularities.getDefaultValue().si,
                 timeGranularities.values().stream().mapToDouble((len) -> len.si).toArray());
 
         // get length-weighted mean speed limit from path to determine cFree and Vc for smoothing
         this.cFree = Speed.min(path.getSpeedLimit(), MAX_C_FREE);
         this.vc = Speed.min(path.getSpeedLimit().times(VC_FACRTOR), MAX_C_FREE);
     }
 
     @Override
     public String toString()
     {
         return "ContourDataSource []";
     }
 
     // ************************************
     // *** PLOT INTERFACING AND GETTERS ***
     // ************************************
 
     /**
      * Returns the path for an {@link AbstractContourPlot} using this {@link ContourDataSource}.
      * @return the path
      */
     GraphPath<? extends LaneData<?>> getPath()
     {
         return this.path;
     }
 
     @Override
     public void addPlot(final AbstractContourPlot<?> contourPlot)
     {
         ContourDataType<?> contourDataType = contourPlot.getContourDataType();
         if (contourDataType instanceof ContourAdditionalDataType<?, ?> type)
         {
             this.additionalData.put(type, null);
         }
         this.dataTypes.add(contourDataType);
         super.addPlot(contourPlot);
     }
 
     /**
      * Returns the available granularities that a linked plot may use.
      * @param dimension space or time
      * @return available granularities that a linked plot may use
      */
     public double[] getGranularities(final Dimension dimension)
     {
         return dimension.getAxis(this).granularities;
     }
 
     /**
      * Returns the selected granularity that a linked plot should use.
      * @param dimension space or time
      * @return granularity that a linked plot should use
      */
     public double getGranularity(final Dimension dimension)
     {
         return dimension.getAxis(this).granularity;
     }
 
     /**
      * Sets the granularity of the plot. This will invalidate the plot triggering a redraw.
      * @param dimension space or time
      * @param granularity granularity in space or time (SI unit)
      */
     public void setGranularity(final Dimension dimension, final double granularity)
     {
         synchronized (this)
         {
             dimension.getAxis(this).setGranularity(granularity);
             if (dimension.equals(Dimension.TIME))
             {
                 getPlots().forEach((p) -> p.offerUpdateInterval(Duration.ofSI(granularity)));
             }
             invalidateTimeSpan();
         }
         fireEvent(GRANULARITY, new Object[] {dimension, granularity});
     }
 
     /**
      * Sets bi-linear interpolation enabled or disabled. This will invalidate the plot triggering a redraw.
      * @param interpolate whether to enable interpolation
      */
     public void setInterpolate(final boolean interpolate)
     {
         boolean did = false;
         synchronized (this)
         {
             if (this.timeAxis.interpolate != interpolate)
             {
                 did = true;
                 this.timeAxis.setInterpolate(interpolate);
                 this.spaceAxis.setInterpolate(interpolate);
             }
             invalidateTimeSpan();
         }
         if (did)
         {
             fireEvent(INTERPOLATE, interpolate);
         }
     }
 
     /**
      * Sets the adaptive smoothing enabled or disabled. This will invalidate the plot triggering a redraw.
      * @param smooth whether to smooth the plot
      */
     public void setSmooth(final boolean smooth)
     {
         boolean did = false;
         synchronized (this)
         {
             if (this.smooth != smooth)
             {
                 did = true;
                 this.smooth = smooth;
             }
             invalidateTimeSpan();
         }
         if (did)
         {
             fireEvent(SMOOTH, smooth);
         }
     }
 
     // ************************
     // *** UPDATING METHODS ***
     // ************************
 
     /**
      * Heart of the data pool. This method is invoked regularly by the worker thread of a plot, either for a scheduled update or
      * due to user input. No two invocations can happen at the same time.
      * <p>
      * This method regularly checks conditions that indicate the update should be interrupted as for example a setting has
      * changed and repainting is required. Whenever a new invalidation causes {@link #invalidateTimeSpan} to be invoked, this
      * method can stop as the full data needs to be recalculated. This can be set by any change of e.g. granularity or
      * smoothing, during the update.
      * @param t time up to which to show data
      */
     @Override
     protected void calculatePaintStateUnsafe(final Duration t)
     {
         Throw.when(getPlots().isEmpty(), IllegalStateException.class, "ContourDataSource is used, but not by a contour plot!");
 
         // Get consistent update context
         UpdateContext uc = snapshotAndPrepare(t);
 
         // Reset data arrays and clear filter upon a redo
         resetDataOnRedo(uc);
 
         // Setup or clear filter data
         configureFilter(uc);
 
         // Ensure capacity
         ensureTimeCapacity(uc.toTimeIndex());
 
         // Process the data
         for (int j = uc.fromTimeIndex(); j <= uc.toTimeIndex(); j++)
         {
             if (!processTimeSlice(j, uc))
             {
                 return; // early stop requested
             }
         }
 
         // Smoothing filter
         applySmoothingIfNeeded(uc);
     }
 
     /**
      * Takes a snapshot of the update context with which the worker thread can run. Then other threads can set properties while
      * the worker is running with a consistent context.
      * @param now current time
      * @return update context
      */
     private UpdateContext snapshotAndPrepare(final Duration now)
     {
         /**
          * This method is executed once at a time by a plot worker thread. Many properties, such as the data, are maintained by
          * this method. Other properties, which other methods can change, are read first in a synchronized block, while those
          * methods are also synchronized.
          */
         boolean redo;
         double timeGranularity;
         double spaceGranularity;
         boolean smooth0;
         boolean interpolate0;
         double timeKernelSize;
         double spaceKernelSize;
         double[] spaceTicks;
         double[] timeTicks;
         int fromSpaceIndex = 0;
         int fromTimeIndex = 0;
         int toTimeIndex;
         double tFromEgtf = 0;
         int skipTime = 0;
         int tSliceFromEgtf = 0;
         double snappedToTime;
         synchronized (this)
         {
             timeGranularity = this.timeAxis.granularity;
             spaceGranularity = this.spaceAxis.granularity;
 
             redo = getAndResetInvalidTimeSpan();
             // snap to granularity
             this.toTime = (double) (timeGranularity * ((int) (now.si / timeGranularity)));
             snappedToTime = this.toTime;
             if (snappedToTime > this.timeAxis.maxValue)
             {
                 this.timeAxis.setMaxValue(snappedToTime);
             }
 
             // save local copies so commands given during this execution can change it for the next execution
             smooth0 = this.smooth && snappedToTime > timeGranularity;
             interpolate0 = this.timeAxis.interpolate;
             // kernel size based on granularity
             timeKernelSize = timeGranularity * 2 * KERNEL_FACTOR;
             spaceKernelSize = spaceGranularity * 2 * KERNEL_FACTOR;
             spaceTicks = Arrays.copyOf(this.spaceAxis.getTicks(), this.spaceAxis.getTicks().length);
             timeTicks = Arrays.copyOf(this.timeAxis.getTicks(), this.timeAxis.getTicks().length);
             if (!redo)
             {
                 // remember where we started, readyItems will be updated but we need to know where we started during the update
                 fromSpaceIndex = (this.readyItems + 1) % this.spaceAxis.getSliceCount();
                 fromTimeIndex = (this.readyItems + 1) / this.spaceAxis.getSliceCount();
             }
             toTimeIndex = ((int) (snappedToTime / timeGranularity)) - (interpolate0 ? 0 : 1);
             if (smooth0)
             {
                 // time of current slice - kernel size, get slice of that time, get time (middle) of that slice
                 tFromEgtf = this.timeAxis.getSliceValue(redo ? 0
                         : this.timeAxis.getValueSlice(this.timeAxis.getSliceValue(fromTimeIndex) - 2 * timeKernelSize));
                 tSliceFromEgtf = this.timeAxis.getValueSlice(tFromEgtf);
 
                 /*
                  * The above time is based on twice the kernel size because the fast implementation only accounts for data on
                  * (and within range of) the output grid. To make sure all data within the kernel size (now-kernel : now) is
                  * correct (given the data available up to now), we need all data in twice that size (now-2*kernel : now). Only
                  * the second half of that (now-kernel : now) should be written in the output data. The value of skipTime makes
                  * overwriteSmoothed() skip the first half (now-2*kernel : now-kernel).
                  */
                 double tFromEgtf2 = this.timeAxis.getSliceValue(
                         redo ? 0 : this.timeAxis.getValueSlice(this.timeAxis.getSliceValue(fromTimeIndex) - timeKernelSize));
                 int nFromEgtf2 = this.timeAxis.getValueSlice(tFromEgtf2);
                 skipTime = nFromEgtf2 - tSliceFromEgtf;
             }
 
             if (redo)
             {
                 this.readyItems = -1;
             }
         }
         return new UpdateContext(redo, timeGranularity, spaceGranularity, smooth0, interpolate0, timeKernelSize,
                 spaceKernelSize, spaceTicks, timeTicks, fromSpaceIndex, fromTimeIndex, toTimeIndex, tFromEgtf, tSliceFromEgtf,
                 skipTime, snappedToTime);
     }
 
     /**
      * Resets data and filter upon a redo.
      * @param uc update context
      */
     private void resetDataOnRedo(final UpdateContext uc)
     {
         if (!uc.redo())
         {
             return;
         }
         int nSpace = uc.spaceTicks().length - 1;
         int nTime = uc.timeTicks().length - 1;
         this.distance = new float[nSpace][nTime];
         this.time = new float[nSpace][nTime];
         for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
         {
             this.additionalData.put(type, new float[nSpace][nTime]);
         }
         this.egtf = null;
     }
 
     /**
      * Configure filter; setting up a filter or clearing any existent filter objects.
      * @param uc update context
      */
     private void configureFilter(final UpdateContext uc)
     {
         if (uc.smooth() && this.egtf == null)
         {
             setupFilter(uc.timeGranularity(), uc.timeKernelSize(), uc.spaceKernelSize());
         }
         else if (!uc.smooth())
         {
             // discard smoothing state
             this.egtf = null;
             this.speedStream = null;
             this.travelTimeStream = null;
             this.travelDistanceStream = null;
             this.additionalStreams.clear();
         }
     }
 
     /**
      * Setup the filter.
      * @param timeGranularity time granularity
      * @param timeKernelSize time kernel size
      * @param spaceKernelSize space kernel size
      */
     private void setupFilter(final double timeGranularity, final double timeKernelSize, final double spaceKernelSize)
     {
         // create the filter
         this.egtf = new Egtf(C_CONG.si, this.cFree.si, DELTA_V.si, this.vc.si);
 
         // create data source and its data streams for speed, distance traveled, time traveled, and additional
         DataSource generic = this.egtf.getDataSource("generic");
         generic.addStream(TypedQuantity.SPEED, Speed.ofSI(1.0), Speed.ofSI(1.0));
         generic.addStreamSI(this.travelTimeQuantity, 1.0, 1.0);
         generic.addStreamSI(this.travelDistanceQuantity, 1.0, 1.0);
         this.speedStream = generic.getStream(TypedQuantity.SPEED);
         this.travelTimeStream = generic.getStream(this.travelTimeQuantity);
         this.travelDistanceStream = generic.getStream(this.travelDistanceQuantity);
         for (ContourAdditionalDataType<?, ?> contourDataType : this.additionalData.keySet())
         {
             this.additionalStreams.put(contourDataType, generic.addStreamSI(contourDataType.getQuantity(), 1.0, 1.0));
         }
 
         // for maximum space and time range, increase sigma and tau by KERNEL_FACTOR, beyond which both kernels diminish
         this.egtf.setKernelSI(spaceKernelSize / KERNEL_FACTOR, timeKernelSize / KERNEL_FACTOR, spaceKernelSize,
                 timeGranularity);
 
         // add listener to provide a filter status update and to possibly stop the filter when the plot is invalidated
         this.egtf.addListener((event) ->
         {
             // check stop (explicit use of property, not locally stored value)
             if (isInvalidTimeSpan())
             {
                 // plots need to be redone
                 Logger.ots().debug("Interrupting EGTF");
                 event.interrupt(); // stop the EGTF
             }
         });
     }
 
     /**
      * Ensure capacity.
      * @param toTimeIndex to time index
      */
     private void ensureTimeCapacity(final int toTimeIndex)
     {
         for (int i = 0; i < this.distance.length; i++)
         {
             this.distance[i] = GraphUtil.ensureCapacity(this.distance[i], toTimeIndex + 1);
             this.time[i] = GraphUtil.ensureCapacity(this.time[i], toTimeIndex + 1);
             for (float[][] add : this.additionalData.values())
             {
                 add[i] = GraphUtil.ensureCapacity(add[i], toTimeIndex + 1);
             }
         }
     }
 
     /**
      * Processing a time slice.
      * @param j time slice index
      * @param uc update context
      * @return false if the processing should be aborted
      */
     private boolean processTimeSlice(final int j, final UpdateContext uc)
     {
         final Duration tFrom = Duration.ofSI(uc.timeTicks()[j]);
         final Duration tTo = Duration.ofSI(uc.timeTicks()[j + 1]);
 
         int fromSpaceIndex = uc.fromSpaceIndex(); // local copy; set to 0 after first cell
 
         for (int i = fromSpaceIndex; i < uc.spaceTicks().length - 1; i++)
         {
             if (handleInterpolationEdges(i, j, uc.interpolate()))
             {
                 this.readyItems++;
                 if (isInvalidTimeSpan())
                 {
                     return false;
                 }
                 continue;
             }
 
             // in next time slice, all of space needs to be processed
             fromSpaceIndex = 0;
 
             // define cell
             Length xFrom = Length.ofSI(uc.spaceTicks()[i]);
             Length xTo = Length.ofSI(Math.min(uc.spaceTicks()[i + 1], this.path.getTotalLength().si));
             CellWindow window = new CellWindow(i, j, xFrom, xTo, tFrom, tTo);
 
             // compute cell totals
             CellTotals totals = aggregateCell(window);
 
             // write cell data
             this.distance[i][j] = (float) totals.distance();
             this.time[i][j] = (float) totals.time();
             for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
             {
                 this.additionalData.get(type)[i][j] = finalizeAdditional(totals.additional(), type);
             }
 
             feedFilterIfNeeded(window, totals, uc);
 
             if (isInvalidTimeSpan())
             {
                 return false; // early stop
             }
             this.readyItems++;
         }
 
         // offer time slice result
         offerPaintState(uc);
 
         return true;
     }
 
     /**
      * Handle edge cases for interpolation.
      * @param i space slice index
      * @param j time slice index
      * @param interpolate whether the data will be interpolated
      * @return whether the cell was processed as an edge case for interpolation
      */
     private boolean handleInterpolationEdges(final int i, final int j, final boolean interpolate)
     {
         if ((j == 0 || i == 0) && interpolate)
         {
             this.distance[i][j] = Float.NaN;
             this.time[i][j] = Float.NaN;
             for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
             {
                 this.additionalData.get(type)[i][j] = Float.NaN;
             }
             return true;
         }
         return false;
     }
 
     /**
      * Aggregate data in a single cell.
      * @param cell cell window
      * @return cell totals
      */
     private CellTotals aggregateCell(final CellWindow cell)
     {
         double totalDistance = 0.0;
         double totalTime = 0.0;
 
         Map<ContourAdditionalDataType<?, ?>, Object> additionalIntermediate = new LinkedHashMap<>();
         for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
         {
             additionalIntermediate.put(type, type.identity());
         }
 
         int nSeries = this.path.getNumberOfSeries();
         for (int series = 0; series < nSeries; series++)
         {
             // gather groups for series
             List<TrajectoryGroup<?>> groups = groupsForSeries(series);
 
             // filter groups for cell
             List<TrajectoryGroup<?>> included = new ArrayList<>();
             List<Length> xStart = new ArrayList<>();
             List<Length> xEnd = new ArrayList<>();
             includedGroupsForCell(groups, cell, included, xStart, xEnd);
 
             // accumulate data
             DistTime distTime = accumulateDistanceAndTime(cell, included, xStart, xEnd);
             totalDistance += distTime.distance();
             totalTime += distTime.time();
             for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
             {
                 addAdditional(additionalIntermediate, type, included, xStart, xEnd, cell.tFrom(), cell.tTo());
             }
         }
 
         // normalize to full cell on single lane so EGTF compares apples to apples
         double length = cell.xTo().si - cell.xFrom().si;
         double norm = this.spaceAxis.granularity / length / nSeries;
         totalDistance *= norm;
         totalTime *= norm;
 
         return new CellTotals(totalDistance, totalTime, additionalIntermediate);
     }
 
     /**
      * Returns trajectory groups for the series.
      * @param series series index
      * @return trajectory groups for the series
      */
     private List<TrajectoryGroup<?>> groupsForSeries(final int series)
     {
         List<TrajectoryGroup<?>> groups = new ArrayList<>();
         for (Section<? extends LaneData<?>> section : getPath().getSections())
         {
             TrajectoryGroup<?> group = this.samplerData.getTrajectoryGroup(section.getSource(series)).orElse(null);
             if (group == null)
             {
                 Logger.ots().error("trajectoryGroup {} is null", series);
             }
             groups.add(group);
         }
         return groups;
     }
 
     /**
      * Filter the trajectory groups regarding the cell. The results are added to the last three input parameters.
      * @param trajectories trajectory groups
      * @param cell cell window
      * @param included list for included trajectories to be stored in
      * @param xStart list of start coordinates for included trajectories
      * @param xEnd list of end coordinates for included trajectories
      */
     private void includedGroupsForCell(final List<TrajectoryGroup<?>> trajectories, final CellWindow cell,
             final List<TrajectoryGroup<?>> included, final List<Length> xStart, final List<Length> xEnd)
     {
         for (int k = 0; k < trajectories.size(); k++)
         {
             TrajectoryGroup<?> tg = trajectories.get(k);
             LaneData<?> lane = tg.getLane();
             Length startDistance = this.path.getStartDistance(this.path.get(k));
             double secStart = startDistance.si;
             double secEnd = secStart + this.path.get(k).length().si;
 
             if (secEnd > cell.xFrom().si && secStart < cell.xTo().si)
             {
                 included.add(tg);
                 double scale = this.path.get(k).length().si / lane.getLength().si;
                 xStart.add(Length.max(cell.xFrom().minus(startDistance).divide(scale), Length.ZERO));
                 xEnd.add(Length.min(cell.xTo().minus(startDistance).divide(scale), tg.getLane().getLength()));
             }
         }
     }
 
     /**
      * Accumulate distance and time of included trajectories.
      * @param cell cell window
      * @param included included trajectories
      * @param xStart list of start coordinates for included trajectories
      * @param xEnd list of end coordinates for included trajectories
      * @return accumulated distance and time
      */
     private DistTime accumulateDistanceAndTime(final CellWindow cell, final List<TrajectoryGroup<?>> included,
             final List<Length> xStart, final List<Length> xEnd)
     {
         double dist = 0.0;
         double tim = 0.0;
         for (int k = 0; k < included.size(); k++)
         {
             TrajectoryGroup<?> tg = included.get(k);
             for (Trajectory<?> tr : tg.getTrajectories())
             {
                 if (!GraphUtil.considerTrajectory(tr, cell.tFrom(), cell.tTo()))
                 {
                     continue;
                 }
                 try
                 {
                     SpaceTimeView v = tr.getSpaceTimeView(xStart.get(k), xEnd.get(k), cell.tFrom(), cell.tTo());
                     dist += v.distance().si;
                     tim += v.time().si;
                 }
                 catch (IllegalArgumentException ex)
                 {
                     Logger.ots().debug(ex, "Unable to generate space-time view x={}..{}, t={}..{}", xStart.get(k), xEnd.get(k),
                             cell.tFrom(), cell.tTo());
                 }
             }
         }
         return new DistTime(dist, tim);
     }
 
     /**
      * Add additional data to stored intermediate result.
      * @param additionalIntermediate intermediate storage map
      * @param contourDataType additional data type
      * @param included trajectories
      * @param xStart start distance per trajectory group
      * @param xEnd end distance per trajectory group
      * @param tFrom start time
      * @param tTo end time
      * @param <I> intermediate data type
      */
     @SuppressWarnings("unchecked")
     private <I> void addAdditional(final Map<ContourAdditionalDataType<?, ?>, Object> additionalIntermediate,
             final ContourAdditionalDataType<?, ?> contourDataType, final List<TrajectoryGroup<?>> included,
             final List<Length> xStart, final List<Length> xEnd, final Duration tFrom, final Duration tTo)
     {
         additionalIntermediate.put(contourDataType, ((ContourAdditionalDataType<?, I>) contourDataType)
                 .processSeries((I) additionalIntermediate.get(contourDataType), included, xStart, xEnd, tFrom, tTo));
     }
 
     /**
      * Stores a finalized result for additional data.
      * @param additionalIntermediate intermediate storage map
      * @param contourDataType additional data type
      * @return finalized results for a cell
      * @param <I> intermediate data type
      */
     @SuppressWarnings("unchecked")
     private <I> float finalizeAdditional(final Map<ContourAdditionalDataType<?, ?>, Object> additionalIntermediate,
             final ContourAdditionalDataType<?, ?> contourDataType)
     {
         return ((ContourAdditionalDataType<?, I>) contourDataType).finalize((I) additionalIntermediate.get(contourDataType))
                 .floatValue();
     }
 
     /**
      * Feed data into the filter if we are smoothing.
      * @param cell cell window
      * @param totals totals in cell
      * @param uc update context
      */
     private void feedFilterIfNeeded(final CellWindow cell, final CellTotals totals, final UpdateContext uc)
     {
         if (!uc.smooth())
         {
             return;
         }
         double xDat = (cell.xFrom().si + cell.xTo().si) / 2.0;
         double tDat = (cell.tFrom().si + cell.tTo().si) / 2.0;
 
         if (this.path.isCircular())
         {
             double pathLength = this.path.getTotalLength().si;
             if (xDat < uc.spaceKernelSize())
             {
                 setDataInEgtf(pathLength + xDat, tDat, totals.distance(), totals.time(), cell.i(), cell.j());
             }
             if (xDat > pathLength - uc.spaceKernelSize())
             {
                 setDataInEgtf(xDat - pathLength, tDat, totals.distance(), totals.time(), cell.i(), cell.j());
             }
         }
         setDataInEgtf(xDat, tDat, totals.distance(), totals.time(), cell.i(), cell.j());
     }
 
     /**
      * Sets data in the EGTF for filtering.
      * @param xDat position of data
      * @param tDat time of data
      * @param totalDistance total distance traveled
      * @param totalTime total time traveled
      * @param i space index in data grid
      * @param j time index in data grid
      */
     private void setDataInEgtf(final double xDat, final double tDat, final double totalDistance, final double totalTime,
             final int i, final int j)
     {
         // speed data is implicit as totalDistance/totalTime, but the EGTF needs it explicitly
         this.egtf.addPointDataSI(this.speedStream, xDat, tDat, totalDistance / totalTime);
         this.egtf.addPointDataSI(this.travelDistanceStream, xDat, tDat, totalDistance);
         this.egtf.addPointDataSI(this.travelTimeStream, xDat, tDat, totalTime);
         for (ContourAdditionalDataType<?, ?> contourDataType : this.additionalStreams.keySet())
         {
             this.egtf.addPointDataSI(this.additionalStreams.get(contourDataType), xDat, tDat,
                     this.additionalData.get(contourDataType)[i][j]);
         }
     }
 
     /**
      * Apply smoothing filter.
      * @param uc update context
      */
     private void applySmoothingIfNeeded(final UpdateContext uc)
     {
         if (!uc.smooth())
         {
             return;
         }
 
         // gather quantities
         Set<Quantity<?, ?>> quantities = new LinkedHashSet<>();
         quantities.add(this.travelDistanceQuantity);
         quantities.add(this.travelTimeQuantity);
         this.additionalData.keySet().forEach((type) -> quantities.add(type.getQuantity()));
 
         // size of space to skip as this space was only used to provide data around edges
         int skipSpace = this.path.isCircular() ? (int) Math.ceil(uc.spaceKernelSize() / uc.spaceGranularity()) : 0;
 
         // do the filtering
         double tTo = uc.snappedToTime();
         if (tTo <= uc.tSliceFromEgtf())
         {
             return;
         }
         Optional<Filter> filter = this.egtf.filterFastSI(uc.spaceTicks()[0] + (0.5 - skipSpace) * uc.spaceGranularity(),
                 uc.spaceGranularity(), uc.spaceTicks()[0] + (-1.5 + uc.spaceTicks().length + skipSpace) * uc.spaceGranularity(),
                 uc.tSliceFromEgtf(), uc.timeGranularity(), tTo, quantities.toArray(new Quantity<?, ?>[0]));
         if (filter.isEmpty())
         {
             return;
         }
 
         // overwrite data with smoothed data
         overwriteSmoothed(this.distance, uc.nFromEgtf(), filter.get().getSI(this.travelDistanceQuantity), uc.skipTime(),
                 skipSpace);
         overwriteSmoothed(this.time, uc.nFromEgtf(), filter.get().getSI(this.travelTimeQuantity), uc.skipTime(), skipSpace);
         for (ContourAdditionalDataType<?, ?> type : this.additionalData.keySet())
         {
             overwriteSmoothed(this.additionalData.get(type), uc.nFromEgtf(), filter.get().getSI(type.getQuantity()),
                     uc.skipTime(), skipSpace);
         }
 
         // notify filter result
         offerPaintState(uc);
     }
 
     /**
      * Offers the current state of data for painting.
      * @param uc update context
      */
     private void offerPaintState(final UpdateContext uc)
     {
         final int nTimeSlices = uc.timeTicks().length - 1;
         final int nSpaceSlices = uc.spaceTicks().length - 1;
         final int n = nTimeSlices * nSpaceSlices;
         final Map<ContourDataType<?>, float[]> dataMap = new LinkedHashMap<>();
         final double area = uc.timeGranularity() * uc.spaceGranularity();
         final int limit = Math.min(this.readyItems + 1, n);
 
         for (ContourDataType<?> dataType : this.dataTypes)
         {
             float[] data = new float[n];
             int i = 0;
             if (dataType instanceof ContourEdieDataType edieType)
             {
                 OUTER: for (int timeSlice = 0; timeSlice < nTimeSlices; timeSlice++)
                 {
                     for (int spaceSlice = 0; spaceSlice < nSpaceSlices; spaceSlice++, i++)
                     {
                         if (i == limit)
                         {
                             break OUTER;
                         }
                         else
                         {
                             data[i] = (float) edieType.calculate(this.time[spaceSlice][timeSlice],
                                     this.distance[spaceSlice][timeSlice], area);
                         }
                     }
                 }
             }
             else
             {
                 ContourAdditionalDataType<?, ?> additionalType = (ContourAdditionalDataType<?, ?>) dataType;
                 float[][] matrix = this.additionalData.get(dataType);
                 float scale = (float) (additionalType.normalize() ? area : 1.0);
                 OUTER: for (int timeSlice = 0; timeSlice < nTimeSlices; timeSlice++)
                 {
                     for (int spaceSlice = 0; spaceSlice < nSpaceSlices; spaceSlice++, i++)
                     {
                         if (i == limit)
                         {
                             break OUTER;
                         }
                         else
                         {
                             data[i] = matrix[i % nSpaceSlices][i / nSpaceSlices] / scale;
                         }
                     }
                 }
             }
             if (i < n - 1)
             {
                 Arrays.fill(data, i, n, Float.NaN);
             }
             dataMap.put(dataType, data);
         }
 
         for (AbstractContourPlot<?> plot : getPlots())
         {
             plot.offerPaintState(new ContourPaintState(dataMap.get(plot.getContourDataType()), uc.spaceGranularity(),
                     uc.timeGranularity(), nSpaceSlices, uc.interpolate(), Duration.ofSI(uc.snappedToTime())));
         }
     }
 
     /**
      * Helper method to fill smoothed data in to raw data.
      * @param raw the raw non-smoothed data
      * @param rawCol column from which onward to fill smoothed data in to the raw data which is used for plotting
      * @param smoothed smoothed data returned by {@code EGTF}
      * @param skipTime slices to skip as this was only part of the smoothed data to include the kernel size
      * @param skipSpace slices to ignore at start and end because of circular graph path (i.e. this was only included for data)
      */
     private void overwriteSmoothed(final float[][] raw, final int rawCol, final double[][] smoothed, final int skipTime,
             final int skipSpace)
     {
         for (int i = 0; i < raw.length; i++)
         {
             int ii = i + skipSpace;
             // can't use System.arraycopy due to float vs double
             for (int j = skipTime; j < smoothed[ii].length; j++)
             {
                 raw[i][j + rawCol] = (float) smoothed[ii][j];
             }
         }
     }
 
     // **********************
     // *** HELPER CLASSES ***
     // **********************
 
     /**
      * Update context.
      * @param redo redo whole time window
      * @param timeGranularity time granularity
      * @param spaceGranularity space granularity
      * @param smooth smooth data
      * @param interpolate interpolate data
      * @param timeKernelSize time kernel size
      * @param spaceKernelSize space kernel size
      * @param spaceTicks space ticks
      * @param timeTicks time ticks
      * @param fromSpaceIndex from space index
      * @param fromTimeIndex from time index
      * @param toTimeIndex tom time index
      * @param tSliceFromEgtf start slice for filter; only meaningful if {@code smooth==true}
      * @param nFromEgtf n start filter; only meaningful if {@code smooth==true}
      * @param skipTime skip time for filter; only meaningful if {@code smooth==true}
      * @param snappedToTime to time that adheres time granularity
      */
     private record UpdateContext(boolean redo, double timeGranularity, double spaceGranularity, boolean smooth,
             boolean interpolate, double timeKernelSize, double spaceKernelSize, double[] spaceTicks, double[] timeTicks,
             int fromSpaceIndex, int fromTimeIndex, int toTimeIndex, double tSliceFromEgtf, int nFromEgtf, int skipTime,
             double snappedToTime)
     {
     }
 
     /**
      * Defines a cell for processing.
      * @param i space slice index
      * @param j time slice index
      * @param xFrom space start coordinate
      * @param xTo space end coordinate
      * @param tFrom time start coordinate
      * @param tTo time end coordinate
      */
     private record CellWindow(int i, int j, Length xFrom, Length xTo, Duration tFrom, Duration tTo)
     {
     }
 
     /**
      * Totals computed in a cell.
      * @param distance distance in cell
      * @param time time in cell
      * @param additional additional data in cell
      */
     private record CellTotals(double distance, double time, Map<ContourAdditionalDataType<?, ?>, Object> additional)
     {
     }
 
     /**
      * Intermediate data storage to accumulate values in cell.
      * @param distance distance in cell from one lane
      * @param time time in cell from one lane
      */
     private record DistTime(double distance, double time)
     {
     }
 
     /**
      * Enum to refer to either the distance or time axis.
      */
     public enum Dimension
     {
         /** Distance axis. */
         DISTANCE
         {
             @Override
             protected Axis getAxis(final ContourDataSource dataPool)
             {
                 return dataPool.spaceAxis;
             }
         },
 
         /** Time axis. */
         TIME
         {
             @Override
             protected Axis getAxis(final ContourDataSource dataPool)
             {
                 return dataPool.timeAxis;
             }
         };
 
         /**
          * Returns the {@code Axis} object.
          * @param dataPool data pool
          * @return axis
          */
         protected abstract Axis getAxis(ContourDataSource dataPool);
     }
 
     /**
      * Class to store and determine axis information such as granularity, ticks, and range.
      */
     static class Axis
     {
         /** Minimum value. */
         private final double minValue;
 
         /** Maximum value. */
         private double maxValue;
 
         /** Selected granularity. */
         private double granularity;
 
         /** Possible granularities. */
         private final double[] granularities;
 
         /** Whether the data pool is set to interpolate. */
         private boolean interpolate = true;
 
         /** Tick values. */
         private double[] ticks;
 
         /**
          * Constructor.
          * @param minValue minimum value
          * @param maxValue maximum value
          * @param granularity initial granularity
          * @param granularities possible granularities
          */
         Axis(final double minValue, final double maxValue, final double granularity, final double[] granularities)
         {
             this.minValue = minValue;
             this.maxValue = maxValue;
             this.granularity = granularity;
             this.granularities = granularities;
         }
 
         /**
          * Sets the maximum value.
          * @param maxValue maximum value
          */
         void setMaxValue(final double maxValue)
         {
             if (this.maxValue != maxValue)
             {
                 this.maxValue = maxValue;
                 this.ticks = null;
             }
         }
 
         /**
          * Sets the granularity.
          * @param granularity granularity
          */
         void setGranularity(final double granularity)
         {
             if (this.granularity != granularity)
             {
                 this.granularity = granularity;
                 this.ticks = null;
             }
         }
 
         /**
          * Returns the ticks, which are calculated if needed.
          * @return ticks
          */
         double[] getTicks()
         {
             if (this.ticks == null)
             {
                 int n = getSliceCount() + 1;
                 this.ticks = new double[n];
                 int di = this.interpolate ? 1 : 0;
                 for (int i = 0; i < n; i++)
                 {
                     if (i == n - 1)
                     {
                         this.ticks[i] = Math.min((i - di) * this.granularity, this.maxValue);
                     }
                     else
                     {
                         this.ticks[i] = (i - di) * this.granularity;
                     }
                 }
             }
             return this.ticks;
         }
 
         /**
          * Calculates the number of slices.
          * @return number of slices
          */
         int getSliceCount()
         {
             return (int) Math.ceil((this.maxValue - this.minValue) / this.granularity) + (this.interpolate ? 1 : 0);
         }
 
         /**
          * Calculates the center value of a slices.
          * @param slice slice number
          * @return center value of the slice
          */
         double getSliceValue(final int slice)
         {
             return this.minValue + (0.5 + slice - (this.interpolate ? 1 : 0)) * this.granularity;
         }
 
         /**
          * Looks up the slice number of the value.
          * @param value value
          * @return slice number
          */
         int getValueSlice(final double value)
         {
             getTicks();
             if (value > this.ticks[this.ticks.length - 1])
             {
                 return this.ticks.length - 1;
             }
             int i = 0;
             while (i < this.ticks.length - 1 && this.ticks[i + 1] < value + 1e-9)
             {
                 i++;
             }
             return i;
         }
 
         /**
          * Sets interpolation, important is it required the data to have an additional row or column.
          * @param interpolate interpolation
          */
         void setInterpolate(final boolean interpolate)
         {
             if (this.interpolate != interpolate)
             {
                 this.interpolate = interpolate;
                 this.ticks = null;
             }
         }
 
         /**
          * Retrieve the interpolate flag.
          * @return true if interpolation is on; false if interpolation is off
          */
         public boolean isInterpolate()
         {
             return this.interpolate;
         }
 
         @Override
         public String toString()
         {
             return "Axis [minValue=" + this.minValue + ", maxValue=" + this.maxValue + ", granularity=" + this.granularity
                     + ", granularities=" + Arrays.toString(this.granularities) + ", interpolate=" + this.interpolate
                     + ", ticks=" + Arrays.toString(this.ticks) + "]";
         }
 
     }
 
     /**
      * Contour data type.
      * @param <Z> value type
      */
     public sealed interface ContourDataType<Z extends Number> permits ContourEdieDataType, ContourAdditionalDataType
     {
 
         /**
          * Adds weighted values to the mean, where the weight of each value in {@code values} is equal to the respective delta
          * in {@code weightDimension}. The last value in {@code values} is ignored. Argument {@code weightDimension} is
          * typically space or time to produce a space-mean or a time-mean.
          * @param values values
          * @param weightDimension weight dimension
          * @param mean mean
          */
         static void weighted(final float[] values, final float[] weightDimension, final ArithmeticMean<Double, Double> mean)
         {
             for (int i = 0; i < values.length - 1; i++)
             {
                 mean.add((double) values[i], (double) (weightDimension[i + 1] - weightDimension[i]));
             }
         }
 
         /**
          * Adds weighted values to the mean, where the weight of each value in {@code values} is equal to the respective delta
          * in {@code weightDimension}. NaN values and the last value in {@code values} are ignored. Argument
          * {@code weightDimension} is typically space or time to produce a space-mean or a time-mean.
          * @param values values
          * @param weightDimension weight dimension
          * @param mean mean
          */
         static void weightedNaN(final float[] values, final float[] weightDimension, final ArithmeticMean<Double, Double> mean)
         {
             for (int i = 0; i < values.length - 1; i++)
             {
                 if (!Float.isNaN(values[i]))
                 {
                     mean.add((double) values[i], (double) (weightDimension[i + 1] - weightDimension[i]));
                 }
             }
         }
 
     }
 
     /**
      * Edie's contour data source types. These are the standard flow, density and speed, calculated based on total distance,
      * total time, and area "time x space" of a cell.
      * @param <Z> value type
      */
     public non-sealed interface ContourEdieDataType<Z extends Number> extends ContourDataType<Z>
     {
 
         /** Contour data type for flow. */
         ContourEdieDataType<Frequency> FLOW = new ContourEdieDataType<>()
         {
             @Override
             public double calculate(final double totalTime, final double totalDistance, final double area)
             {
                 return totalDistance / area;
             }
         };
 
         /** Contour data type for density. */
         ContourEdieDataType<LinearDensity> DENSITY = new ContourEdieDataType<>()
         {
             @Override
             public double calculate(final double totalTime, final double totalDistance, final double area)
             {
                 return totalTime / area;
             }
         };
 
         /** Contour data type for speed. */
         ContourEdieDataType<Speed> SPEED = new ContourEdieDataType<>()
         {
             @Override
             public double calculate(final double totalTime, final double totalDistance, final double area)
             {
                 return totalDistance / totalTime;
             }
         };
 
         /**
          * Calculate the value.
          * @param totalTime total trajectory time in area
          * @param totalDistance total trajectory distance in area
          * @param area area "time x space" of space-time cell
          * @return calculated value
          */
         double calculate(double totalTime, double totalDistance, double area);
 
     }
 
     /**
      * Interface for data types of which a contour plot can be made. Using this class, the data pool can determine and store
      * cell values for a variable set of additional data types (besides total distance, total time and speed).
      * @param <Z> value type
      * @param <I> intermediate data type
      */
     public non-sealed interface ContourAdditionalDataType<Z extends Number, I> extends ContourDataType<Z>
     {
 
         /**
          * Returns the initial value for intermediate result.
          * @return I, initial intermediate value
          */
         I identity();
 
         /**
          * Calculate value from provided trajectories that apply to a single grid cell on a single series (lane).
          * @param intermediate intermediate value of previous series, starts as the identity
          * @param trajectories trajectories, all groups overlap the requested space-time
          * @param xFrom start location of cell on the section
          * @param xTo end location of cell on the section
          * @param tFrom start time of cell
          * @param tTo end time of cell
          * @return intermediate value
          */
         I processSeries(I intermediate, List<TrajectoryGroup<?>> trajectories, List<Length> xFrom, List<Length> xTo,
                 Duration tFrom, Duration tTo);
 
         /**
          * Returns the final value of the intermediate result after all lanes.
          * @param intermediate intermediate result after all lanes
          * @return final value
          */
         Z finalize(I intermediate);
 
         /**
          * Returns the quantity that is being plotted on the z-axis for the EGTF filter.
          * @return quantity that is being plotted on the z-axis for the EGTF filter
          */
         Quantity<Z, ?> getQuantity();
 
         /**
          * Returns whether the data type needs normalization by the area "space x time".
          * @return whether the data type needs normalization by the area "space x time"
          */
         boolean normalize();
 
     }
 
 }