Handle dual access grids

- Switch to use a thresholds array instead of a single threshold in `AnalysisWorkerTask`s.
- Check that only a single grid is being created (they currently would have the same filename)
- Create a `MultiGridResultWriter` for dual access grids. Store dual access results in `accessibilityValues` in `RegionalWorkResult`.

src/main/java/com/conveyal/analysis/controllers/RegionalAnalysisController.java

@@ -39,7 +39,6 @@
 import java.io.FileOutputStream;
 import java.io.IOException;
 import java.io.InputStream;
-import java.io.OutputStream;
 import java.net.URI;
 import java.nio.file.FileSystem;
 import java.nio.file.FileSystems;
@@ -599,6 +598,19 @@ private RegionalAnalysis createRegionalAnalysis (Request req, Response res) thro
                 "recordTimes can only be used with a single destination pointset, which must be freeform (non-grid)."
             );
         }
+        if (task.includeTemporalDensity) {
+            checkArgument(
+                    task.dualAccessibilityThresholds != null &&
+                            task.dualAccessibilityThresholds.length > 0,
+                    "dualAccessibilityThresholds not specified when includeTemporalDensity is enabled."
+            );
+            if (task.originPointSet == null) {
+                checkArgument(
+                        !task.recordAccessibility,
+                        "Accessibility and dual accessibility grids cannot be created simultaneously."
+                );
+            }
+        }
 
         // TODO remove duplicate fields from RegionalAnalysis that are already in the nested task.
         // The RegionalAnalysis should just be a minimal wrapper around the template task, adding the origin point set.

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

@@ -164,17 +164,16 @@ public class AnalysisRequest {
 
     /**
      * 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.
+     * to the broker. Requires a destination pointset to be specified. If an origin pointset is specified the results
+     * will be collated to CSV.
      */
     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
+     * 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;
+    public int[] dualAccessibilityThresholds;
 
     /**
      * Freeform (untyped) flags for enabling experimental, undocumented, or arcane behavior in backend or workers.
@@ -289,7 +288,7 @@ public void populateTask (AnalysisWorkerTask task, UserPermissions userPermissio
         }
 
         task.includeTemporalDensity = includeTemporalDensity;
-        task.dualAccessibilityThreshold = dualAccessibilityThreshold;
+        task.dualAccessibilityThresholds = dualAccessibilityThresholds;
         task.flags = flags;
         task.csvResultOptions = csvResultOptions;
     }

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

@@ -139,15 +139,21 @@ public MultiOriginAssembler (RegionalAnalysis regionalAnalysis, Job job, FileSto
 
             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.");
+                    // Gridded origins.
+                    resultWriters.add(new MultiGridResultWriter(
+                            regionalAnalysis,
+                            job.templateTask,
+                            job.templateTask.dualAccessibilityThresholds.length,
+                            fileStorage
+                    ));
                 } else {
                     // Freeform origins.
                     // Output includes temporal density of opportunities and optionally dual accessibility.
-                    resultWriters.add(new TemporalDensityCsvResultWriter(job.templateTask, fileStorage));
+                    resultWriters.add(new TemporalDensityCsvResultWriter(
+                            job.templateTask,
+                            job.templateTask.dualAccessibilityThresholds[0],
+                            fileStorage
+                    ));
                 }
             }
 

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

@@ -15,11 +15,11 @@
  */
 public class TemporalDensityCsvResultWriter extends CsvResultWriter {
 
-    private final int dualThreshold;
+    private final int dualAccessibilityThreshold;
 
-    public TemporalDensityCsvResultWriter(RegionalTask task, FileStorage fileStorage) throws IOException {
+    public TemporalDensityCsvResultWriter(RegionalTask task, int dualAccessibilityThreshold, FileStorage fileStorage) throws IOException {
         super(task, CsvResultType.TDENSITY, fileStorage);
-        dualThreshold = task.dualAccessibilityThreshold;
+        this.dualAccessibilityThreshold = dualAccessibilityThreshold;
     }
 
     @Override
@@ -68,7 +68,7 @@ public Iterable<String[]> rowValues (RegionalWorkResult workResult) {
                 int m = 0;
                 double sum = 0;
                 // Find smallest integer M such that we have already reached D destinations after M minutes of travel.
-                while (sum < dualThreshold && m < 120) {
+                while (sum < dualAccessibilityThreshold && m < 120) {
                     sum += densitiesPerMinute[m];
                     m += 1;
                 }

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

@@ -26,8 +26,7 @@ public class TemporalDensityResult {
     // Internal state fields
 
     private final PointSet[] destinationPointSets;
-    private final int nPercentiles;
-    private final int opportunityThreshold;
+    private final int[] dualAccessibilityThresholds;
 
     // Externally visible fields for accumulating results
 
@@ -44,16 +43,15 @@ public TemporalDensityResult(AnalysisWorkerTask task) {
             "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];
+        this.dualAccessibilityThresholds = task.dualAccessibilityThresholds;
+        this.opportunitiesPerMinute = new double[destinationPointSets.length][task.percentiles.length][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++) {
+            for (int p = 0; p < opportunitiesPerMinute.length; p++) {
                 if (travelTimePercentilesSeconds[p] == UNREACHED) {
                     break; // If any percentile is unreached, all higher ones are also unreached.
                 }
@@ -66,31 +64,30 @@ public void recordOneTarget (int target, int[] travelTimePercentilesSeconds) {
     }
 
     /**
-     * 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.
+     * Writes dual accessibility values (in minutes) to our standard access grid format. The value returned (for
+     * an origin) is the number of minutes required to reach a threshold number of opportunities (specified by
+     * task.dualAccessibilityThresholds) in the specified destination layer at a given percentile of
+     * travel time. If the threshold cannot be reached in less than 120 minutes, returns 0.
      */
-    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;
+    public int[][][] calculateDualAccessibilityGrid() {
+        int nPointSets = opportunitiesPerMinute.length;
+        int nPercentiles = opportunitiesPerMinute.length > 0 ? opportunitiesPerMinute[0].length : 0;
+        int nThresholds = dualAccessibilityThresholds.length;
+        int[][][] dualAccessibilityGrid = new int[nPointSets][nPercentiles][nThresholds];
+        for (int i = 0; i < nPointSets; i++) {
+            for (int j = 0; j < nPercentiles; j++) {
+                for (int k = 0; k < nThresholds; k++) {
+                    int threshold = dualAccessibilityThresholds[k];
+                    int minute = 0;
+                    double sum = 0;
+                    while (sum < threshold && minute < 120) {
+                        sum += opportunitiesPerMinute[i][j][minute];
+                        minute += 1;
                     }
+                    dualAccessibilityGrid[i][j][k] = minute;
                 }
             }
         }
-        return result;
-    }
-
-    public int[][] minutesToReachOpportunities() {
-        return minutesToReachOpportunities(opportunityThreshold);
+        return dualAccessibilityGrid;
     }
-
 }

src/main/java/com/conveyal/r5/analyst/cluster/AnalysisWorkerTask.java

@@ -114,6 +114,7 @@ public abstract class AnalysisWorkerTask extends ProfileRequest {
      * opportunities from this origin (known technically as "dual accessibility").
      */
     public int dualAccessibilityThreshold = 0;
+    public int[] dualAccessibilityThresholds;
 
     /** Whether to build a histogram of travel times to each destination, generally used in testing and debugging. */
     public boolean recordTravelTimeHistograms = false;

src/main/java/com/conveyal/r5/analyst/cluster/RegionalWorkResult.java

@@ -65,7 +65,12 @@ public RegionalWorkResult(OneOriginResult result, RegionalTask task) {
         this.travelTimeValues = result.travelTimes == null ? null : result.travelTimes.values;
         this.accessibilityValues = result.accessibility == null ? null : result.accessibility.getIntValues();
         this.pathResult = result.paths == null ? null : result.paths.summarizeIterations(PathResult.Stat.MINIMUM);
-        this.opportunitiesPerMinute = result.density == null ? null : result.density.opportunitiesPerMinute;
+        if (result.density != null) {
+            this.opportunitiesPerMinute = result.density.opportunitiesPerMinute;
+            if (task.originPointSet == null) {
+                this.accessibilityValues = result.density.calculateDualAccessibilityGrid();
+            }
+        }
         // TODO checkTravelTimeInvariants, checkAccessibilityInvariants to verify that values are monotonically increasing
     }