Update Simpson Desert tests

Use single destination instead of Mercator grid.
Expand and add Javadoc comments.
Add better measure of distribution goodness-of-fit.
Apply percentile check and goodness-of-fit in every test.
Increase default Monte Carlo draws to get smoother histograms.

src/main/java/com/conveyal/r5/analyst/TravelTimeComputer.java

@@ -77,19 +77,21 @@ public OneOriginResult computeTravelTimes() {
 
         // Find the set of destinations for a travel time calculation, not yet linked to the street network, and with
         // no associated opportunities. By finding the extents and destinations up front, we ensure the exact same
-        // destination pointset is used for all steps below.
+        // destination PointSet is used for all steps below.
         // This reuses the logic for finding the appropriate grid size and linking, which is now in the NetworkPreloader.
         // We could change the preloader to retain these values in a compound return type, to avoid repetition here.
         PointSet destinations;
-
         if (request instanceof  RegionalTask
             && !request.makeTauiSite
             && request.destinationPointSets[0] instanceof FreeFormPointSet
         ) {
-            // Freeform; destination pointset was set by handleOneRequest in the main AnalystWorker
+            // Freeform destinations. Destination PointSet was set by handleOneRequest in the main AnalystWorker.
+            destinations = request.destinationPointSets[0];
+        } else if (request.destinationPointSets != null) {
+            LOG.warn("ONLY VALID IN TESTING: Using PointSet object embedded in request where this is not standard.");
             destinations = request.destinationPointSets[0];
         } else {
-            // Gridded (non-freeform) destinations. The extents are found differently in regional and single requests.
+            // Gridded (non-freeform) destinations. This method finds them differently for regional and single requests.
             WebMercatorExtents destinationGridExtents = request.getWebMercatorExtents();
             // Make a WebMercatorGridPointSet with the right extents, referring to the network's base grid and linkage.
             destinations = AnalysisWorkerTask.gridPointSetCache.get(destinationGridExtents, network.fullExtentGridPointSet);

src/test/java/com/conveyal/r5/analyst/network/Distribution.java

@@ -1,9 +1,12 @@
 package com.conveyal.r5.analyst.network;
 
 import com.conveyal.r5.analyst.cluster.TravelTimeResult;
+import com.google.common.base.Preconditions;
 
 import java.util.Arrays;
 
+import static java.lang.Math.pow;
+import static java.lang.Math.sqrt;
 import static org.junit.jupiter.api.Assertions.assertTrue;
 
 /**
@@ -90,6 +93,10 @@ public void normalize () {
         }
     }
 
+    /**
+     * Print a text-based representation of the distribution to standard out.
+     * There is another method to show the distribution in a graphical plot window.
+     */
     public void illustrate () {
         final int width = 50;
         double max = Arrays.stream(masses).max().getAsDouble();
@@ -102,6 +109,12 @@ public void illustrate () {
         }
     }
 
+    /**
+     * Given a percentile such as 25 or 50, find the x bin at which that percentile is situated in this Distribution,
+     * i.e. the lowest (binned or discretized) x value for which the cumulative probability is at least percentile.
+     * In common usage: find the lowest whole-minute travel time for which the cumulative probability is greater than
+     * the supplied percentile.
+     */
     public int findPercentile (int percentile) {
         double sum = 0;
         double threshold = percentile / 100d;
@@ -123,6 +136,10 @@ public static void main (String[] args) {
         out.illustrate();
     }
 
+    /**
+     * @return the probability mass situated at a particular x value (the probability density for a particular minute
+     *         when these are used in the usual way as 1-minute bins).
+     */
     public double probabilityOf (int x) {
         if (x < skip) {
             return 0;
@@ -200,10 +217,18 @@ public static Distribution fromTravelTimeResult (TravelTimeResult travelTimeResu
     }
 
     /**
-     * Find the probability mass of the overlapping region of the two distributions. The amount of "misplaced"
-     * probability is one minus overlap. Overlap is slightly more straightforward to calculate directly than mismatch.
+     * Find the probability mass of the overlapping region of the two distributions. This can be used to determine
+     * whether two distributions, often a theoretical one and an observed one, are sufficiently similar to one another.
+     * Overlapping here means in both dimensions, travel time (horizontal) and probability density (vertical).
+     * Proceeding bin by bin through both distributions in parallel, the smaller of the two values for each bin is
+     * accumulated into the total. The amount of "misplaced" probability (located in the wrong bin in the observed
+     * distribution relative to the theoretical one) is one minus overlap. Overlap is slightly more straightforward
+     * to calculate directly than mismatch. This method is not sensitive to how evenly the error is distributed
+     * across the domain. We should prefer using a measure that emphasizes larger errors and compensates for the
+     * magnitude of the predicted values.
      */
     public double overlap (Distribution other) {
+        // TODO This min is not necessary. The overlap is by definition fully within the domain of either Distribution.
         int iMin = Math.min(this.skip, other.skip);
         int iMax = Math.max(this.fullWidth(), other.fullWidth());
         double sum = 0;
@@ -216,12 +241,45 @@ public double overlap (Distribution other) {
         return sum;
     }
 
+    /**
+     * An ad-hoc measure of goodness of fit vaguely related to Pearson's chi-squared or root-mean-square error.
+     * Traditional measures used in fitting probability distributions like Pearson's have properties that deal poorly
+     * with our need to tolerate small horizontal shifts
+     * in the results (due to the destination grid being not precisely aligned with our street corner grid).
+     * Another way to deal with this would be to ensure there is no horizontal shift, by measuring travel times at
+     * exactly the right places instead of on a grid.
+     */
+    public double weightedSquaredError (Distribution other) {
+        double sum = 0;
+        // This is kind of ugly because we're only examining bins with nonzero probability (to avoid div by zero).
+        // Observed data in a region with predicted zero probability should be an automatic fail for the model.
+        for (int i = this.skip; i < this.fullWidth(); i++) {
+            double pe = this.probabilityOf(i);
+            double po = other.probabilityOf(i);
+            Preconditions.checkState(pe >= 0); // Ensure non-negative for good measure.
+            if (pe == 0) {
+                System.out.println("Zero (expected probability; skipping.");
+                continue;
+            }
+            // Errors are measured relative to the expected values, and stronger deviations emphasized by squaring.
+            // Measuring relative to expected density compensates for the case where it is spread over a wider domain.
+            sum += pow(po - pe, 2) / pe;
+        }
+        System.out.println("Squared error: " + sum);
+        System.out.println("Root Squared error: " + sqrt(sum));
+        return sum;
+    }
+
     public void assertSimilar (Distribution observed) {
+        double squaredError = this.weightedSquaredError(observed);
+        showChartsIfEnabled(observed);
+        assertTrue(squaredError < 0.02, String.format("Error metric too high at at %3f", squaredError));
+    }
+
+    public void showChartsIfEnabled (Distribution observed) {
         if (SHOW_CHARTS) {
             DistributionChart.showChart(this, observed);
         }
-        double overlapPercent = this.overlap(observed) * 100;
-        assertTrue(overlapPercent >= 95, String.format("Overlap less than 95%% at %3f", overlapPercent));
     }
 
     // This is ugly, it should be done some other way e.g. firstNonzero
@@ -249,4 +307,19 @@ public void trim () {
         }
         masses = Arrays.copyOfRange(masses, firstNonzero, lastNonzero + 1);
     }
+
+    /**
+     * Here we are performing two related checks for a bit of redundancy and to check different parts of the system:
+     * checking percentiles drawn from the observed distribution, as well as the full histogram of the distribution.
+     * This double comparison could be done automatically with a method like Distribution.assertSimilar(TravelTimeResult).
+     * @param destination the flattened 1D index into the pointset, which will be zero for single freeform points.
+     */
+    public void multiAssertSimilar(TravelTimeResult travelTimes, int destination) {
+        // Check a goodness-of-fit metric on the observed distribution relative to this distribution.
+        Distribution observed = Distribution.fromTravelTimeResult(travelTimes, destination);
+        this.assertSimilar(observed);
+        // Check that percentiles extracted from observed are similar to those predicted by this distribution.
+        int[] travelTimePercentiles = travelTimes.getTarget(destination);
+        DistributionTester.assertExpectedDistribution(this, travelTimePercentiles);
+    }
 }

src/test/java/com/conveyal/r5/analyst/network/DistributionChart.java

@@ -61,10 +61,13 @@ public JFreeChart createChart (Distribution... distributions) {
         return chart;
     }
 
+    // Note that the points are placed at the minute boundary, though the numbers represent densities over one minute.
+    // They should probably be represented as filled bars across the minute or as points midway across the minute.
     private static TimeSeriesCollection createTimeSeriesDataset (Distribution... distributions) {
         TimeSeriesCollection dataset = new TimeSeriesCollection();
+        int d = 0;
         for (Distribution distribution : distributions) {
-            TimeSeries ts = new TimeSeries("X");
+            TimeSeries ts = new TimeSeries("Distribution " + (d++));
             for (int i = distribution.skip(); i < distribution.fullWidth(); i++) {
                 double p = distribution.probabilityOf(i);
                 ts.add(new Minute(i, 0, 1, 1, 2000), p);

src/test/java/com/conveyal/r5/analyst/network/DistributionTester.java

@@ -24,6 +24,11 @@ public static void assertUniformlyDistributed (int[] sortedPercentiles, int min,
         }
     }
 
+    /**
+     * Given an expected distribution of travel times at a destination and the standard five percentiles of travel time
+     * at that same destination as computed by our router, check that the computed values seem to be drawn from the
+     * theoretically correct distribution.
+     */
     public static void assertExpectedDistribution (Distribution expectedDistribution, int[] values) {
         for (int p = 0; p < PERCENTILES.length; p++) {
             int expected = expectedDistribution.findPercentile(PERCENTILES[p]);

src/test/java/com/conveyal/r5/analyst/network/GridRoute.java

@@ -21,7 +21,7 @@ public class GridRoute {
     public Orientation orientation;
     public boolean bidirectional;
 
-    /** Explicit departure times from first stop; if set, startHour and endHour will be ignored*/
+    /** Explicit departure times from first stop; if set, startHour and endHour will be ignored. */
     public int[] startTimes;
 
     /** Override default hop times. Map of (trip, stopAtStartOfHop) to factor by which default hop is multiplied. */

src/test/java/com/conveyal/r5/analyst/network/GridSinglePointTaskBuilder.java

@@ -29,6 +29,7 @@
  */
 public class GridSinglePointTaskBuilder {
 
+    public static final int DEFAULT_MONTE_CARLO_DRAWS = 4800; // 40 per minute over a two hour window.
     private final GridLayout gridLayout;
     private final AnalysisWorkerTask task;
 
@@ -50,7 +51,7 @@ public GridSinglePointTaskBuilder (GridLayout gridLayout) {
         // In single point tasks all 121 cutoffs are required (there is a check).
         task.cutoffsMinutes = IntStream.rangeClosed(0, 120).toArray();
         task.decayFunction = new StepDecayFunction();
-        task.monteCarloDraws = 1200; // Ten per minute over a two hour window.
+        task.monteCarloDraws = DEFAULT_MONTE_CARLO_DRAWS;
         // By default, traverse one block in a round predictable number of seconds.
         task.walkSpeed = gridLayout.streetGridSpacingMeters / gridLayout.walkBlockTraversalTimeSeconds;
         // Record more detailed information to allow comparison to theoretical travel time distributions.
@@ -128,11 +129,29 @@ public GridSinglePointTaskBuilder uniformOpportunityDensity (double density) {
         return this;
     }
 
+    /**
+     * When trying to verify more complex distributions, the Monte Carlo approach may introduce too much noise.
+     * Increasing the number of draws will yield a better approximation of the true travel time distribution
+     * (while making the tests run slower).
+     */
     public GridSinglePointTaskBuilder monteCarloDraws (int draws) {
         task.monteCarloDraws = draws;
         return this;
     }
 
+    /**
+     * This eliminates any difficulty estimating the final segment of egress, walking from the street to a gridded
+     * travel time sample point. Although egress time is something we'd like to test too, it is not part of the transit
+     * routing we're concentrating on here, and will vary as the Simpson Desert street grid does not align with our
+     * web Mercator grid pixels. Using a single measurement point also greatly reduces the amount of travel time
+     * histograms that must be computed and retained, improving the memory and run time cost of tests.
+     */
+    public GridSinglePointTaskBuilder singleFreeformDestination(int x, int y) {
+        FreeFormPointSet ps = new FreeFormPointSet(gridLayout.getIntersectionLatLon(x, y));
+        task.destinationPointSets = new PointSet[] { ps };
+        return this;
+    }
+
     public AnalysisWorkerTask build () {
         return task;
     }