Temporal Opportunity Density and Dual Accessibility

src/main/java/com/conveyal/analysis/models/AnalysisRequest.java

@@ -161,6 +161,21 @@ public class AnalysisRequest {
      */
     public ChaosParameters injectFault;
 
+    /**
+     * Whether to include the number of opportunities reached during each minute of travel in results sent back
+     * to the broker. Requires both an origin and destination pointset to be specified, and in the case of regional
+     * analyses the origins must be non-gridded, and results will be collated to CSV.
+     * It should be possible to enable regional results for gridded origins as well.
+     */
+    public boolean includeTemporalDensity = false;
+
+    /**
+     * If this is set to a value above zero, report the amount of time needed to reach the given number of
+     * opportunities from this origin (known technically as "dual accessibility").
+     */
+    public int dualAccessibilityThreshold = 0;
+
+
     /**
      * Create the R5 `Scenario` from this request.
      */
@@ -265,6 +280,9 @@ public void populateTask (AnalysisWorkerTask task, UserPermissions userPermissio
                 throw new IllegalArgumentException("Must be admin user to inject faults.");
             }
         }
+
+        task.includeTemporalDensity = includeTemporalDensity;
+        task.dualAccessibilityThreshold = dualAccessibilityThreshold;
     }
 
     private EnumSet<LegMode> getEnumSetFromString (String s) {

src/main/java/com/conveyal/analysis/results/CsvResultType.java

@@ -5,5 +5,5 @@
  * do serve to enumerate the acceptable parameters coming over the HTTP API.
  */
 public enum CsvResultType {
-    ACCESS, TIMES, PATHS
+    ACCESS, TIMES, PATHS, TDENSITY
 }

src/main/java/com/conveyal/analysis/results/CsvResultWriter.java

@@ -12,6 +12,8 @@
 import java.io.FileWriter;
 import java.io.IOException;
 
+import static com.google.common.base.Preconditions.checkArgument;
+import static com.google.common.base.Preconditions.checkNotNull;
 import static com.google.common.base.Preconditions.checkState;
 
 /**
@@ -55,6 +57,7 @@ public abstract class CsvResultWriter extends BaseResultWriter implements Region
      */
     CsvResultWriter (RegionalTask task, FileStorage fileStorage) throws IOException {
         super(fileStorage);
+        checkArgument(task.originPointSet != null, "CsvResultWriters require FreeFormPointSet origins.");
         super.prepare(task.jobId);
         this.fileName = task.jobId + "_" + resultType() +".csv";
         BufferedWriter bufferedWriter = new BufferedWriter(new FileWriter(bufferFile));

src/main/java/com/conveyal/analysis/results/MultiOriginAssembler.java

@@ -109,8 +109,10 @@ public MultiOriginAssembler (RegionalAnalysis regionalAnalysis, Job job, FileSto
 
             if (job.templateTask.recordAccessibility) {
                 if (job.templateTask.originPointSet != null) {
+                    // Freeform origins - create CSV regional analysis results
                     resultWriters.add(new AccessCsvResultWriter(job.templateTask, fileStorage));
                 } else {
+                    // Gridded origins - create gridded regional analysis results
                     resultWriters.add(new MultiGridResultWriter(regionalAnalysis, job.templateTask, fileStorage));
                 }
             }
@@ -123,6 +125,20 @@ public MultiOriginAssembler (RegionalAnalysis regionalAnalysis, Job job, FileSto
                 resultWriters.add(new PathCsvResultWriter(job.templateTask, fileStorage));
             }
 
+            if (job.templateTask.includeTemporalDensity) {
+                if (job.templateTask.originPointSet == null) {
+                    // Gridded origins. The full temporal density information is probably too voluminous to be useful.
+                    // We might want to record a grid of dual accessibility values, but this will require some serious
+                    // refactoring of the GridResultWriter.
+                    // if (job.templateTask.dualAccessibilityThreshold > 0) { ... }
+                    throw new IllegalArgumentException("Temporal density of opportunities cannot be recorded for gridded origin points.");
+                } else {
+                    // Freeform origins.
+                    // Output includes temporal density of opportunities and optionally dual accessibility.
+                    resultWriters.add(new TemporalDensityCsvResultWriter(job.templateTask, fileStorage));
+                }
+            }
+
             checkArgument(job.templateTask.makeTauiSite || notNullOrEmpty(resultWriters),
                 "A non-Taui regional analysis should always create at least one grid or CSV file.");
 

src/main/java/com/conveyal/analysis/results/TemporalDensityCsvResultWriter.java

@@ -0,0 +1,90 @@
+package com.conveyal.analysis.results;
+
+import com.conveyal.file.FileStorage;
+import com.conveyal.r5.analyst.cluster.RegionalTask;
+import com.conveyal.r5.analyst.cluster.RegionalWorkResult;
+
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ * This handles collating regional results into CSV files containing temporal opportunity density
+ * (number of opportunities reached in each one-minute interval, the derivative of step-function accessibility)
+ * as well as "dual" accessibility (the amount of time needed to reach n opportunities).
+ */
+public class TemporalDensityCsvResultWriter extends CsvResultWriter {
+
+    private final int dualThreshold;
+
+    public TemporalDensityCsvResultWriter(RegionalTask task, FileStorage fileStorage) throws IOException {
+        super(task, fileStorage);
+        dualThreshold = task.dualAccessibilityThreshold;
+    }
+
+    @Override
+    public CsvResultType resultType () {
+        return CsvResultType.TDENSITY;
+    }
+
+    @Override
+    public String[] columnHeaders () {
+        List<String> headers = new ArrayList<>();
+        // The ids of the freeform origin point and destination set
+        headers.add("originId");
+        headers.add("destId");
+        headers.add("percentile");
+        for (int m = 0; m < 120; m += 1) {
+            // The opportunity density over travel minute m
+            headers.add(Integer.toString(m));
+        }
+        // The number of minutes needed to reach d destination opportunities
+        headers.add("D" + dualThreshold);
+        return headers.toArray(new String[0]);
+    }
+
+    @Override
+    protected void checkDimension (RegionalWorkResult workResult) {
+        checkDimension(
+            workResult, "destination pointsets",
+            workResult.opportunitiesPerMinute.length, task.destinationPointSetKeys.length
+        );
+        for (double[][] percentilesForPointset : workResult.opportunitiesPerMinute) {
+            checkDimension(workResult, "percentiles", percentilesForPointset.length, task.percentiles.length);
+            for (double[] minutesForPercentile : percentilesForPointset) {
+                checkDimension(workResult, "minutes", minutesForPercentile.length, 120);
+            }
+        }
+    }
+
+    @Override
+    public Iterable<String[]> rowValues (RegionalWorkResult workResult) {
+        List<String[]> rows = new ArrayList<>();
+        String originId = task.originPointSet.getId(workResult.taskId);
+        for (int d = 0; d < task.destinationPointSetKeys.length; d++) {
+            double[][] percentilesForDestPointset = workResult.opportunitiesPerMinute[d];
+            for (int p = 0; p < task.percentiles.length; p++) {
+                List<String> row = new ArrayList<>(125);
+                row.add(originId);
+                row.add(task.destinationPointSetKeys[d]);
+                row.add(Integer.toString(p));
+                // One density value for each of 120 minutes
+                double[] densitiesPerMinute = percentilesForDestPointset[p];
+                for (int m = 0; m < 120; m++) {
+                    row.add(Double.toString(densitiesPerMinute[m]));
+                }
+                // One dual accessibility value
+                int m = 0;
+                double sum = 0;
+                while (sum < dualThreshold && m < 120) {
+                    sum += densitiesPerMinute[m];
+                    m += 1;
+                }
+                row.add(Integer.toString(m >= 120 ? -1 : m));
+                rows.add(row.toArray(new String[row.size()]));
+            }
+        }
+        return rows;
+    }
+
+}

src/main/java/com/conveyal/r5/OneOriginResult.java

@@ -1,6 +1,7 @@
 package com.conveyal.r5;
 
 import com.conveyal.r5.analyst.AccessibilityResult;
+import com.conveyal.r5.analyst.TemporalDensityResult;
 import com.conveyal.r5.analyst.cluster.PathResult;
 import com.conveyal.r5.analyst.cluster.TravelTimeResult;
 
@@ -20,10 +21,18 @@ public class OneOriginResult {
 
     public final PathResult paths;
 
-    public OneOriginResult(TravelTimeResult travelTimes, AccessibilityResult accessibility, PathResult paths) {
+    public final TemporalDensityResult density;
+
+    public OneOriginResult(
+            TravelTimeResult travelTimes,
+            AccessibilityResult accessibility,
+            PathResult paths,
+            TemporalDensityResult density
+    ) {
         this.travelTimes = travelTimes;
         this.accessibility = accessibility;
         this.paths = paths;
+        this.density = density;
     }
 
 }

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

@@ -0,0 +1,95 @@
+package com.conveyal.r5.analyst;
+
+import com.conveyal.r5.analyst.cluster.AnalysisWorkerTask;
+import com.google.common.base.Preconditions;
+
+import static com.conveyal.r5.common.Util.notNullOrEmpty;
+import static com.conveyal.r5.profile.FastRaptorWorker.UNREACHED;
+
+/**
+ * An instance of this is included in a OneOriginResult for reporting how many opportunities are encountered during each
+ * minute of travel. If we use more than one destination point set they are already constrained to all be aligned with
+ * the same number of destinations.
+ *
+ * The data retained here feed into three different kinds of results: "Dual" accessibility (the number of opportunities
+ * reached in a given number of minutes of travel time); temporal opportunity density (analogous to a probability density
+ * function, how many opportunities are encountered during each minute of travel, whose integral is the cumulative
+ * accessibility curve).
+ *
+ * Originally this class was tracking the identity of the N nearest points rather than just binning them by travel time.
+ * This is more efficient in cases where N is small, and allows retaining the one-second resolution. However currently
+ * there does not seem to be much demand among users for this level of detail, so it has been removed in the interest
+ * of simplicity and maintainability. See issue 884 for more comments on implementation trade-offs.
+ */
+public class TemporalDensityResult {
+
+    // Internal state fields
+
+    private final PointSet[] destinationPointSets;
+    private final int nPercentiles;
+    private final int opportunityThreshold;
+
+    // Externally visible fields for accumulating results
+
+    /**
+     * The temporal density of opportunities. For each destination set, for each percentile, for each minute of
+     * travel from 0 to 120, the number of opportunities reached in travel times from i (inclusive) to i+1 (exclusive).
+     */
+    public final double[][][] opportunitiesPerMinute;
+
+    public TemporalDensityResult(AnalysisWorkerTask task) {
+        Preconditions.checkArgument(
+            notNullOrEmpty(task.destinationPointSets),
+            "Dual accessibility requires at least one destination pointset."
+        );
+        this.destinationPointSets = task.destinationPointSets;
+        this.nPercentiles = task.percentiles.length;
+        this.opportunityThreshold = task.dualAccessibilityThreshold;
+        this.opportunitiesPerMinute = new double[destinationPointSets.length][nPercentiles][120];
+    }
+
+    public void recordOneTarget (int target, int[] travelTimePercentilesSeconds) {
+        // Increment histogram bin for the number of minutes of travel by the number of opportunities at the target.
+        for (int d = 0; d < destinationPointSets.length; d++) {
+            PointSet dps = destinationPointSets[d];
+            for (int p = 0; p < nPercentiles; p++) {
+                if (travelTimePercentilesSeconds[p] == UNREACHED) {
+                    break; // If any percentile is unreached, all higher ones are also unreached.
+                }
+                int m = travelTimePercentilesSeconds[p] / 60;
+                if (m <= 120) {
+                    opportunitiesPerMinute[d][p][m] += dps.getOpportunityCount(target);
+                }
+            }
+        }
+    }
+
+    /**
+     * Calculate "dual" accessibility from the accumulated temporal opportunity density array.
+     * @param n the threshold quantity of opportunities
+     * @return the minimum whole number of minutes necessary to reach n opportunities,
+     *         for each destination set and percentile of travel time.
+     */
+    public int[][] minutesToReachOpportunities(int n) {
+        int[][] result = new int[destinationPointSets.length][nPercentiles];
+        for (int d = 0; d < destinationPointSets.length; d++) {
+            for (int p = 0; p < nPercentiles; p++) {
+                result[d][p] = -1;
+                double count = 0;
+                for (int m = 0; m < 120; m++) {
+                    count += opportunitiesPerMinute[d][p][m];
+                    if (count >= n) {
+                        result[d][p] = m + 1;
+                        break;
+                    }
+                }
+            }
+        }
+        return result;
+    }
+
+    public int[][] minutesToReachOpportunities() {
+        return minutesToReachOpportunities(opportunityThreshold);
+    }
+
+}

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

@@ -46,6 +46,9 @@ public class TravelTimeReducer {
     /** Retains the paths to one or all destinations, for recording in CSV or reporting in the UI. */
     private PathResult pathResult = null;
 
+    /** Represents how many destinations are reached in each minute of travel from this origin. */
+    private TemporalDensityResult temporalDensityResult = null;
+
     /** If we are calculating accessibility, the PointSets containing opportunities. */
     private PointSet[] destinationPointSets;
 
@@ -121,6 +124,8 @@ public TravelTimeReducer (AnalysisWorkerTask task, TransportNetwork network) {
         this.destinationPointSets = task.destinationPointSets;
         if (task instanceof TravelTimeSurfaceTask) {
             calculateTravelTimes = true;
+            // In single-point analyses, destination pointsets may be missing if the user has not selected one in the
+            // UI, or if the user has selected the step decay function instead of one of the other decay functions.
             calculateAccessibility = notNullOrEmpty(task.destinationPointSets);
         } else {
             // Maybe we should define recordAccessibility and recordTimes on the common superclass AnalysisWorkerTask.
@@ -140,6 +145,9 @@ public TravelTimeReducer (AnalysisWorkerTask task, TransportNetwork network) {
         if (task.includePathResults) {
             pathResult = new PathResult(task, network.transitLayer);
         }
+        if (task.includeTemporalDensity) {
+            temporalDensityResult = new TemporalDensityResult(task);
+        }
 
         // Validate and copy the travel time cutoffs, converting them to seconds to avoid repeated multiplication
         // in tight loops. Also find the points where the decay function reaches zero for these cutoffs.
@@ -271,6 +279,9 @@ private void recordTravelTimePercentilesForTarget (int target, int[] travelTimeP
                 }
             }
         }
+        if (temporalDensityResult != null) {
+            temporalDensityResult.recordOneTarget(target, travelTimePercentilesSeconds);
+        }
     }
 
     /**
@@ -329,7 +340,7 @@ private int convertToMinutes (int timeSeconds) {
      * origin point is not connected to the street network.
      */
     public OneOriginResult finish () {
-        return new OneOriginResult(travelTimeResult, accessibilityResult, pathResult);
+        return new OneOriginResult(travelTimeResult, accessibilityResult, pathResult, temporalDensityResult);
     }
 
     /**