904 lines
42 KiB
Java
904 lines
42 KiB
Java
package propertymap;
|
||
|
||
import com.conveyal.r5.OneOriginResult;
|
||
import com.conveyal.r5.analyst.FreeFormPointSet;
|
||
import com.conveyal.r5.analyst.LinkageCache;
|
||
import com.conveyal.r5.analyst.PointSet;
|
||
import com.conveyal.r5.analyst.StreetTimesAndModes;
|
||
import com.conveyal.r5.analyst.TravelTimeComputer;
|
||
import com.conveyal.r5.analyst.WebMercatorExtents;
|
||
import com.conveyal.r5.analyst.cluster.PathResult;
|
||
import com.conveyal.r5.analyst.cluster.RegionalTask;
|
||
import com.conveyal.r5.analyst.cluster.TravelTimeResult;
|
||
import com.conveyal.r5.api.util.LegMode;
|
||
import com.conveyal.r5.api.util.TransitModes;
|
||
import com.conveyal.r5.kryo.KryoNetworkSerializer;
|
||
import com.conveyal.r5.profile.StreetMode;
|
||
import com.conveyal.r5.streets.StreetLayer;
|
||
import com.conveyal.r5.streets.VertexStore;
|
||
import com.conveyal.gtfs.model.Stop;
|
||
import com.conveyal.r5.transit.TransitLayer;
|
||
import com.conveyal.r5.transit.TransportNetwork;
|
||
import com.conveyal.r5.transit.path.RouteSequence;
|
||
import com.google.common.collect.Multimap;
|
||
import org.locationtech.jts.geom.Coordinate;
|
||
import org.locationtech.jts.geom.Envelope;
|
||
import org.locationtech.jts.index.strtree.STRtree;
|
||
|
||
import java.io.File;
|
||
import java.time.LocalDate;
|
||
import java.util.ArrayList;
|
||
import java.util.Arrays;
|
||
import java.util.BitSet;
|
||
import java.util.EnumSet;
|
||
import java.util.List;
|
||
|
||
/** R5 routing: network loading, spatial filtering, travel time computation. */
|
||
public class Router {
|
||
|
||
// A stop counts as "should have linked but didn't" only when a walkable street
|
||
// sits within R5's own link radius of it (streetWithinLinkRadius). This is
|
||
// deliberately NOT a GB bounding box: the OSM extract is England-only while the
|
||
// BODS GTFS is GB-wide, so ~59k Welsh/Scottish stops legitimately have no street
|
||
// nearby and must not be counted as linking failures.
|
||
//
|
||
// MAX_UNLINKED_STOP_FRACTION bounds the fraction of stops that DO have a nearby
|
||
// street yet still fail to link — corrupt coordinates or a whole mode displaced
|
||
// (the Tottenham Court Road failure mode). It needs per-stop coordinates and the
|
||
// street index, so it only runs on a fresh build.
|
||
private static final double MAX_UNLINKED_STOP_FRACTION = 0.10;
|
||
// CATASTROPHIC_UNLINKED_FRACTION bounds the raw unlinked fraction over ALL stops
|
||
// and runs on every load (fresh or cached). The legitimate England-vs-GB residue
|
||
// is ~19%; only a grossly wrong or truncated extract (almost nothing links) blows
|
||
// past this — the one failure a coordinate-less cached load can still catch.
|
||
private static final double CATASTROPHIC_UNLINKED_FRACTION = 0.50;
|
||
|
||
private static final int ZOOM = 9; // R5 enforces range 9-12
|
||
private static final int MAX_GRID_CELLS = 4_900_000; // under R5's 5M limit
|
||
// 30-minute peak window: RAPTOR cost is linear in (toTime-fromTime)/60.
|
||
// best_minutes (5th percentile) is the best of these 30 minute-shifted departures.
|
||
private static final int DEPARTURE_FROM_TIME = 7 * 3600 + 45 * 60; // 07:45
|
||
private static final int DEPARTURE_TO_TIME = 8 * 3600 + 15 * 60; // 08:15
|
||
private static final int MAX_TRIP_DURATION_MINUTES = 90;
|
||
// Transit-only: cap walk access/egress at 20 min to shrink the egress
|
||
// street subgraph PerTargetPropagater walks per stop.
|
||
private static final int TRANSIT_MAX_WALK_TIME_MIN = 20;
|
||
|
||
// Hard R5 limit when path recording is enabled (PathResult internals).
|
||
// Larger is better here: each chunk forces R5 to rebuild the egress cost
|
||
// table (~334k stop linkages), so fewer chunks per origin = fewer rebuilds.
|
||
private static final int PATH_MAX_DESTINATIONS = 5000;
|
||
|
||
// Per-chunk destination cap for non-transit direct modes (car/bicycle/walking).
|
||
// London car origins filter to ~1M postcodes within 150km. Without a cap, each
|
||
// chunk's per-task LinkedPointSet + FreeFormPointSet allocate ~50-100 MB and
|
||
// R5's StreetRouter scratch state stacks across concurrent workers, OOMing the
|
||
// heap. 150k caps per-chunk transient memory at ~5-10 MB; chunk count for
|
||
// London goes from 1 to ~7, adding ~10-20% wall-clock per origin via repeated
|
||
// Dijkstra. Walking has so few dests this is a no-op.
|
||
private static final int DIRECT_MAX_DESTINATIONS = 150_000;
|
||
|
||
// Percentile indices in R5 result arrays (order must match task.percentiles in buildTask)
|
||
private static final int PERCENTILE_BEST = 0; // 5th percentile (transit only)
|
||
private static final int PERCENTILE_MEDIAN = 1; // 50th percentile (transit: index 1, others: index 0)
|
||
|
||
/** Result of computing travel times for a single origin with spatial pre-filtering. */
|
||
record FilteredResult(int[] originalIndices, short[] times, short[] bestTimes, String[] journeys) {}
|
||
|
||
/**
|
||
* Global transit tile: a destination subset bundled with the FreeFormPointSet
|
||
* R5 routes against. Reused across origins so R5's LinkageCache (and the
|
||
* expensive EgressCostTable) is built once per tile, not once per origin × chunk.
|
||
*/
|
||
record PostcodeTile(
|
||
FreeFormPointSet pointSet,
|
||
WebMercatorExtents extents,
|
||
int[] originalIndices,
|
||
double minLat, double maxLat, double minLon, double maxLon
|
||
) {}
|
||
|
||
/** True for any transit variant (transit, transit-no-bus, transit-no-change, …). */
|
||
static boolean isTransitMode(String mode) {
|
||
return mode.startsWith("transit");
|
||
}
|
||
|
||
/** Max plausible travel radius in km for {@link #MAX_TRIP_DURATION_MINUTES}-minute trips. */
|
||
static double maxRadiusKm(String mode) {
|
||
if (isTransitMode(mode)) return 150;
|
||
return switch (mode) {
|
||
case "car" -> 150;
|
||
case "bicycle" -> 60;
|
||
case "walking" -> 12;
|
||
default -> throw new IllegalArgumentException("Unknown mode: " + mode);
|
||
};
|
||
}
|
||
|
||
/**
|
||
* Transit variant configuration. {@code maxRides} is the number of transit legs:
|
||
* 1 = walk-transit-walk (no change), 2 = one change, 3 = two changes.
|
||
* {@code excludeBus} drops {@link TransitModes#BUS} from the allowed mode set.
|
||
*/
|
||
private record TransitConfig(int maxRides, boolean excludeBus) {}
|
||
|
||
private static TransitConfig transitConfigFor(String mode) {
|
||
return switch (mode) {
|
||
case "transit" -> new TransitConfig(3, false);
|
||
case "transit-no-bus" -> new TransitConfig(3, true);
|
||
case "transit-no-change" -> new TransitConfig(1, false);
|
||
case "transit-no-change-no-bus" -> new TransitConfig(1, true);
|
||
case "transit-one-change" -> new TransitConfig(2, false);
|
||
case "transit-one-change-no-bus" -> new TransitConfig(2, true);
|
||
default -> throw new IllegalArgumentException("Unknown transit mode: " + mode);
|
||
};
|
||
}
|
||
|
||
/**
|
||
* Load or build the transport network with Kryo caching.
|
||
* The returned network is read-only after buildDistanceTables, safe for concurrent use.
|
||
*
|
||
* The evictable LinkageCache is left small (32 entries) because non-transit modes
|
||
* create one huge per-origin LinkedPointSet each (~1M dests for car @ 150km radius).
|
||
* Caching 1024 such entries OOMs the heap. Transit tile linkages instead go into
|
||
* the unevictable {@code linkageMap} via {@link #preloadTransitTileLinkages} after
|
||
* tiles are built. That map has no count limit and is checked first on lookup.
|
||
*/
|
||
static TransportNetwork loadNetwork(String dataDir, String cacheDir) throws Exception {
|
||
// Must be set BEFORE the TransportNetwork is deserialized, since its LinkageCache
|
||
// is constructed (and sized) during that deserialization. 32 fits the working
|
||
// set of non-transit per-origin linkages without exhausting heap.
|
||
LinkageCache.LINKAGE_CACHE_SIZE = 32;
|
||
|
||
System.err.println("Loading transport network...");
|
||
File cacheFile = new File(cacheDir, "network.dat");
|
||
TransportNetwork network;
|
||
|
||
if (cacheFile.exists()) {
|
||
System.err.println(" Loading cached network from " + cacheFile);
|
||
network = KryoNetworkSerializer.read(cacheFile);
|
||
} else {
|
||
System.err.println(" Building network (first run, takes a few minutes)...");
|
||
network = TransportNetwork.fromDirectory(new File(dataDir));
|
||
new File(cacheDir).mkdirs();
|
||
KryoNetworkSerializer.write(network, cacheFile);
|
||
System.err.println(" Cached to " + cacheFile);
|
||
}
|
||
|
||
validateTransitNetwork(network, cacheFile, dataDir);
|
||
|
||
System.err.println(" Building distance tables...");
|
||
network.transitLayer.buildDistanceTables(null);
|
||
System.err.println(" Network ready");
|
||
return network;
|
||
}
|
||
|
||
private static void validateTransitNetwork(TransportNetwork network, File cacheFile, String dataDir) {
|
||
TransitLayer transitLayer = network.transitLayer;
|
||
int stops = transitLayer == null ? 0 : transitLayer.getStopCount();
|
||
int routes = transitLayer == null || transitLayer.routes == null ? 0 : transitLayer.routes.size();
|
||
int patterns = transitLayer == null || transitLayer.tripPatterns == null ? 0 : transitLayer.tripPatterns.size();
|
||
int services = transitLayer == null || transitLayer.services == null ? 0 : transitLayer.services.size();
|
||
|
||
if (stops == 0 || routes == 0 || patterns == 0) {
|
||
throw new IllegalStateException(String.format(
|
||
"R5 network has no usable transit data (stops=%d, routes=%d, patterns=%d). "
|
||
+ "The cache at %s was likely built without GTFS. Ensure %s contains GTFS .zip files, "
|
||
+ "then delete %s and rerun.",
|
||
stops, routes, patterns, cacheFile.getPath(), dataDir, cacheFile.getPath()));
|
||
}
|
||
|
||
System.err.printf(" Transit: %,d stops, %,d routes, %,d patterns, %,d services%n",
|
||
stops, routes, patterns, services);
|
||
|
||
validateStopLinkage(transitLayer, network.streetLayer);
|
||
}
|
||
|
||
/**
|
||
* Verify transit stops are connected to the walkable street network. A stop
|
||
* whose coordinate lands nowhere near a street gets streetVertexForStop == -1
|
||
* (R5's "unlinked" sentinel); it is then silently unroutable — no access,
|
||
* egress, or distance table — so trains pass through but nobody can board or
|
||
* alight (the Tottenham Court Road failure mode, upstream of the coordinate
|
||
* fix).
|
||
*
|
||
* The signal for a systemic break is NOT the raw unlinked count: the BODS feed
|
||
* is GB-wide while the OSM extract is England-only, so every Welsh/Scottish stop
|
||
* (~59k, ~19%) is legitimately unlinkable and must be ignored. Instead we count
|
||
* only stops that HAVE a walkable street within R5's link radius yet still fail
|
||
* to link — the coordinates-corrupt / mode-displaced failure mode — and bound
|
||
* that fraction (MAX_UNLINKED_STOP_FRACTION). A separate, coarser bound on the
|
||
* raw unlinked fraction (CATASTROPHIC_UNLINKED_FRACTION) still catches a grossly
|
||
* wrong or truncated extract where almost nothing links.
|
||
*
|
||
* streetVertexForStop is persisted in the network cache, so the raw unlinked
|
||
* count is available on both fresh and cached loads. Per-stop coordinates and the
|
||
* street spatial index are only present on a fresh fromDirectory build, so the
|
||
* precise nearby-street check runs only then; a cached load applies just the
|
||
* catastrophic bound (the network was already checked precisely when first built).
|
||
*/
|
||
private static void validateStopLinkage(TransitLayer transitLayer, StreetLayer streetLayer) {
|
||
int n = transitLayer.getStopCount();
|
||
// stopForIndex is transient: on a cached (Kryo) load it comes back empty
|
||
// (size 0), not null, so require it to be fully populated (fresh build)
|
||
// before indexing into it; otherwise fall back to the catastrophic bound.
|
||
boolean haveCoords = transitLayer.stopForIndex != null
|
||
&& transitLayer.stopForIndex.size() == n;
|
||
int unlinked = 0; // all unlinked stops, wherever they are
|
||
int linkableStops = 0; // stops with a street within link radius (should link)
|
||
int unlinkedNearStreets = 0; // ... of those, the ones that still failed to link
|
||
List<String> sample = new ArrayList<>();
|
||
for (int i = 0; i < n; i++) {
|
||
boolean linked = transitLayer.streetVertexForStop.get(i) != -1;
|
||
if (!linked) {
|
||
unlinked++;
|
||
}
|
||
if (!haveCoords) {
|
||
continue;
|
||
}
|
||
Stop stop = transitLayer.stopForIndex.get(i);
|
||
if (stop == null) {
|
||
continue;
|
||
}
|
||
// A linked stop necessarily had a street nearby; only probe the index for
|
||
// unlinked ones (also far cheaper than probing all ~316k stops).
|
||
boolean streetNear = linked
|
||
|| streetWithinLinkRadius(streetLayer, stop.stop_lat, stop.stop_lon);
|
||
if (streetNear) {
|
||
linkableStops++;
|
||
}
|
||
if (!linked && streetNear) {
|
||
unlinkedNearStreets++;
|
||
if (sample.size() < 25) {
|
||
String id = transitLayer.stopIdForIndex.get(i);
|
||
String name = transitLayer.stopNames != null && i < transitLayer.stopNames.size()
|
||
? transitLayer.stopNames.get(i) : "";
|
||
sample.add(String.format("%s '%s' (%.5f, %.5f)",
|
||
id, name, stop.stop_lat, stop.stop_lon));
|
||
}
|
||
}
|
||
}
|
||
|
||
double globalFraction = n == 0 ? 0.0 : (double) unlinked / n;
|
||
double anomalyFraction = (haveCoords && linkableStops > 0)
|
||
? (double) unlinkedNearStreets / linkableStops : 0.0;
|
||
|
||
if (haveCoords) {
|
||
System.err.printf(
|
||
" Stop linkage: %,d/%,d stops unlinked (%.2f%%); of the %,d with a street "
|
||
+ "within %.0f m, %,d failed to link (%.2f%%). Stops with no nearby street "
|
||
+ "(e.g. Welsh/Scottish stops absent from an England-only OSM extract) are "
|
||
+ "expected-unlinked and excluded.%n",
|
||
unlinked, n, globalFraction * 100.0, linkableStops,
|
||
StreetLayer.LINK_RADIUS_METERS, unlinkedNearStreets, anomalyFraction * 100.0);
|
||
for (String s : sample) {
|
||
System.err.println(" unlinked near streets: " + s);
|
||
}
|
||
} else {
|
||
System.err.printf(
|
||
" Stop linkage: %,d/%,d stops unlinked (%.2f%%) (coords unavailable on cached "
|
||
+ "load; the precise nearby-street check ran at build time)%n",
|
||
unlinked, n, globalFraction * 100.0);
|
||
}
|
||
|
||
// Catastrophic bound, applied on every load: a grossly wrong or truncated
|
||
// extract leaves the great majority of stops unlinked. The legitimate
|
||
// England-vs-GB residue (~19%) is well under this.
|
||
if (globalFraction > CATASTROPHIC_UNLINKED_FRACTION) {
|
||
throw new IllegalStateException(String.format(
|
||
"R5 left %.1f%% of transit stops unlinked (%,d of %,d), past the catastrophic "
|
||
+ "bound of %.0f%%. The OSM extract is likely wrong or truncated (almost "
|
||
+ "nothing links); unlinked stops are silently unroutable.",
|
||
globalFraction * 100.0, unlinked, n, CATASTROPHIC_UNLINKED_FRACTION * 100.0));
|
||
}
|
||
|
||
// Precise bound, fresh build only: stops that HAD a nearby street but still
|
||
// didn't link — corrupt coordinates or a displaced mode (the TCR failure mode).
|
||
if (haveCoords && anomalyFraction > MAX_UNLINKED_STOP_FRACTION) {
|
||
throw new IllegalStateException(String.format(
|
||
"R5 failed to link %.1f%% of the transit stops that have a street within %.0f m "
|
||
+ "(%,d of %,d). This is a systemic linking failure (corrupt stop coordinates, "
|
||
+ "or a whole mode displaced), not the expected residue of stops outside the "
|
||
+ "OSM extract; unlinked stops are silently unroutable.",
|
||
anomalyFraction * 100.0, StreetLayer.LINK_RADIUS_METERS,
|
||
unlinkedNearStreets, linkableStops));
|
||
}
|
||
}
|
||
|
||
/**
|
||
* True if a walkable street edge sits within R5's stop link radius of (lat, lon).
|
||
* Mirrors the spatial query R5's own linker makes, so "no edge here" means R5
|
||
* could not have linked the stop either (it lies outside the OSM extract), while
|
||
* "edge here but the stop is unlinked" is a genuine anomaly. The probe envelope is
|
||
* built in WGS84 degrees, then converted to the fixed-point space the edge index
|
||
* uses — findEdgesInEnvelope does no conversion of its own.
|
||
*/
|
||
private static boolean streetWithinLinkRadius(StreetLayer streetLayer, double lat, double lon) {
|
||
double degLat = StreetLayer.LINK_RADIUS_METERS / 111_320.0;
|
||
double cosLat = Math.max(Math.cos(Math.toRadians(lat)), 1e-6);
|
||
double degLon = StreetLayer.LINK_RADIUS_METERS / (111_320.0 * cosLat);
|
||
Envelope env = new Envelope(lon - degLon, lon + degLon, lat - degLat, lat + degLat);
|
||
return !streetLayer.findEdgesInEnvelope(VertexStore.envelopeToFixed(env)).isEmpty();
|
||
}
|
||
|
||
static void validateTransitServices(TransportNetwork network, LocalDate date) {
|
||
BitSet activeServices = network.transitLayer.getActiveServicesForDate(date);
|
||
if (activeServices.cardinality() == 0) {
|
||
throw new IllegalStateException("R5 network has transit data, but no active services on "
|
||
+ date + ". Rebuild property-data/transit from current feeds or choose a date covered by GTFS.");
|
||
}
|
||
System.err.printf(" Active transit services on %s: %,d%n", date, activeServices.cardinality());
|
||
}
|
||
|
||
/**
|
||
* Estimate routing workload for an origin: count of postcodes within mode radius.
|
||
* Cheap STRtree bbox query; used as the LPT sort key for scheduling.
|
||
*/
|
||
@SuppressWarnings("unchecked")
|
||
static int estimateWorkload(STRtree postcodeIndex, double originLat, double originLon, double maxRadiusKm) {
|
||
double degLat = maxRadiusKm / 111.0;
|
||
double degLon = maxRadiusKm / (111.0 * Math.cos(Math.toRadians(originLat)));
|
||
Envelope env = new Envelope(
|
||
originLon - degLon, originLon + degLon,
|
||
originLat - degLat, originLat + degLat);
|
||
return postcodeIndex.query(env).size();
|
||
}
|
||
|
||
/**
|
||
* Build an STRtree spatial index over postcode points. Forces an initial query
|
||
* to trigger lazy build, so the returned tree is safe for concurrent queries.
|
||
*/
|
||
static STRtree buildPostcodeIndex(double[] lats, double[] lons) {
|
||
STRtree tree = new STRtree();
|
||
for (int i = 0; i < lats.length; i++) {
|
||
tree.insert(new Envelope(lons[i], lons[i], lats[i], lats[i]), Integer.valueOf(i));
|
||
}
|
||
// Force build (otherwise the first concurrent query races on lazy init)
|
||
tree.query(new Envelope(0, 0, 0, 0));
|
||
return tree;
|
||
}
|
||
|
||
/**
|
||
* Build WALK linkages for every transit tile and store them as unevictable on
|
||
* the network's LinkageCache. Subsequent transit routing calls get cache hits
|
||
* regardless of how many per-origin car/bike/walk linkages cycle through the
|
||
* evictable LRU. Logs progress because this is multi-minute work.
|
||
*/
|
||
static void preloadTransitTileLinkages(TransportNetwork network, List<PostcodeTile> tiles) {
|
||
System.err.printf("Pre-building WALK linkages for %,d transit tiles (unevictable)...%n", tiles.size());
|
||
long t0 = System.currentTimeMillis();
|
||
int n = tiles.size();
|
||
for (int i = 0; i < n; i++) {
|
||
network.linkageCache.buildUnevictableLinkage(
|
||
tiles.get(i).pointSet(), network.streetLayer, StreetMode.WALK);
|
||
if ((i + 1) % 25 == 0 || i + 1 == n) {
|
||
double secs = (System.currentTimeMillis() - t0) / 1000.0;
|
||
System.err.printf(" %,d/%,d tile linkages built (%.1fs elapsed)%n", i + 1, n, secs);
|
||
}
|
||
}
|
||
}
|
||
|
||
/**
|
||
* Build global transit tiles ONCE from all postcodes. Each tile holds a
|
||
* FreeFormPointSet that is reused across every transit origin, so R5's
|
||
* LinkageCache hits and the EgressCostTable is built once per tile rather
|
||
* than once per (origin × chunk). This is the dominant transit speed-up.
|
||
*
|
||
* Tiles are sized to the R5 path-result hard limit (PATH_MAX_DESTINATIONS=5000)
|
||
* so the same tiles serve both path-recording and non-path transit requests.
|
||
*/
|
||
static List<PostcodeTile> buildGlobalTransitTiles(double[] lats, double[] lons) {
|
||
int n = lats.length;
|
||
int[] sorted = sortIndicesByLat(lats);
|
||
|
||
// Global lon span sets gridWidth: all tiles share the same horizontal extent
|
||
// bound, so each tile is a horizontal band of postcodes.
|
||
double minLon = Double.MAX_VALUE, maxLon = -Double.MAX_VALUE;
|
||
for (double lon : lons) {
|
||
minLon = Math.min(minLon, lon);
|
||
maxLon = Math.max(maxLon, lon);
|
||
}
|
||
int totalPixels = 256 << ZOOM;
|
||
int gridWidth = lonToPixel(maxLon, totalPixels) - lonToPixel(minLon, totalPixels) + 1;
|
||
int maxHeight = MAX_GRID_CELLS / gridWidth;
|
||
|
||
List<PostcodeTile> tiles = new ArrayList<>();
|
||
int start = 0;
|
||
while (start < n) {
|
||
int end = start + 1;
|
||
int topPixel = latToPixel(lats[sorted[start]], totalPixels);
|
||
|
||
while (end < n) {
|
||
if (end - start >= PATH_MAX_DESTINATIONS) break;
|
||
int bottomPixel = latToPixel(lats[sorted[end]], totalPixels);
|
||
if (Math.abs(bottomPixel - topPixel) + 1 > maxHeight) break;
|
||
end++;
|
||
}
|
||
|
||
tiles.add(buildTile(lats, lons, sorted, start, end));
|
||
start = end;
|
||
}
|
||
return tiles;
|
||
}
|
||
|
||
/**
|
||
* Filter destinations by distance, build chunks, compute travel times for one origin.
|
||
* Returns only the filtered subset indices and their travel times.
|
||
*
|
||
* Transit uses {@code globalTiles} (shared FreeFormPointSets → linkage cache hits);
|
||
* other modes use per-origin filter+chunk (no expensive linkage to amortize, and
|
||
* tiling would force routing to many irrelevant destinations).
|
||
*/
|
||
static FilteredResult computeForOrigin(
|
||
TransportNetwork network,
|
||
STRtree postcodeIndex,
|
||
List<PostcodeTile> globalTiles,
|
||
double[] allLats, double[] allLons,
|
||
double originLat, double originLon,
|
||
String mode, LocalDate date, boolean enablePaths) {
|
||
|
||
if (isTransitMode(mode)) {
|
||
return computeTransit(network, globalTiles, originLat, originLon, mode, date, enablePaths);
|
||
}
|
||
|
||
double maxRadius = maxRadiusKm(mode);
|
||
|
||
// 1. Filter destinations by bounding box (STRtree query)
|
||
int[] filtered = filterByDistance(postcodeIndex, originLat, originLon, maxRadius);
|
||
if (filtered.length == 0) {
|
||
return new FilteredResult(new int[0], new short[0], null, null);
|
||
}
|
||
|
||
// 2. Extract filtered coordinate arrays
|
||
double[] fLats = new double[filtered.length];
|
||
double[] fLons = new double[filtered.length];
|
||
for (int i = 0; i < filtered.length; i++) {
|
||
fLats[i] = allLats[filtered[i]];
|
||
fLons[i] = allLons[filtered[i]];
|
||
}
|
||
|
||
// 3. Build per-origin chunks. Cap at DIRECT_MAX_DESTINATIONS so car at 150km
|
||
// radius (~1M dests for London) gets split into ~7 manageable chunks
|
||
// instead of one giant LinkedPointSet allocation.
|
||
List<DestinationChunk> chunks = buildDestinationChunks(fLats, fLons, DIRECT_MAX_DESTINATIONS);
|
||
|
||
// 4. Compute travel times
|
||
short[][] allTimes = computeTravelTimesDirect(network, chunks, originLat, originLon, mode, fLats.length, date);
|
||
return new FilteredResult(filtered, allTimes[0], null, null);
|
||
}
|
||
|
||
/**
|
||
* Transit routing path: route from origin to every global tile whose bbox intersects
|
||
* the origin's max-radius bbox. Reuses tile FreeFormPointSets for R5 LinkageCache hits.
|
||
* {@code mode} selects the transit variant (rides cap, bus exclusion).
|
||
*/
|
||
private static FilteredResult computeTransit(
|
||
TransportNetwork network, List<PostcodeTile> globalTiles,
|
||
double originLat, double originLon, String mode, LocalDate date, boolean enablePaths) {
|
||
|
||
double maxRadius = maxRadiusKm(mode);
|
||
double degLat = maxRadius / 111.0;
|
||
double degLon = maxRadius / (111.0 * Math.cos(Math.toRadians(originLat)));
|
||
double oMinLat = originLat - degLat, oMaxLat = originLat + degLat;
|
||
double oMinLon = originLon - degLon, oMaxLon = originLon + degLon;
|
||
|
||
List<PostcodeTile> selected = new ArrayList<>();
|
||
int totalDests = 0;
|
||
for (PostcodeTile tile : globalTiles) {
|
||
if (tile.maxLat() < oMinLat || tile.minLat() > oMaxLat) continue;
|
||
if (tile.maxLon() < oMinLon || tile.minLon() > oMaxLon) continue;
|
||
selected.add(tile);
|
||
totalDests += tile.originalIndices().length;
|
||
}
|
||
|
||
if (selected.isEmpty()) {
|
||
return new FilteredResult(new int[0], new short[0], new short[0],
|
||
enablePaths ? new String[0] : null);
|
||
}
|
||
|
||
int[] outIndices = new int[totalDests];
|
||
short[] medianTimes = new short[totalDests];
|
||
short[] bestTimes = new short[totalDests];
|
||
Arrays.fill(medianTimes, (short) -1);
|
||
Arrays.fill(bestTimes, (short) -1);
|
||
String[] journeys = enablePaths ? new String[totalDests] : null;
|
||
|
||
int offset = 0;
|
||
for (PostcodeTile tile : selected) {
|
||
int tileLen = tile.originalIndices().length;
|
||
System.arraycopy(tile.originalIndices(), 0, outIndices, offset, tileLen);
|
||
|
||
RegionalTask task = buildTaskForTile(tile, originLat, originLon, mode, date, enablePaths);
|
||
TravelTimeComputer computer = new TravelTimeComputer(task, network);
|
||
OneOriginResult result = computer.computeTravelTimes();
|
||
|
||
TravelTimeResult tt = result.travelTimes;
|
||
if (tt == null) {
|
||
throw new RuntimeException("R5 returned null travelTimes for tile with " + tileLen + " destinations");
|
||
}
|
||
int[][] values = tt.getValues();
|
||
if (values.length < 2) {
|
||
throw new RuntimeException("R5 returned " + values.length + " percentiles, expected 2");
|
||
}
|
||
for (int i = 0; i < tileLen; i++) {
|
||
if (values[PERCENTILE_BEST][i] != Integer.MAX_VALUE) {
|
||
bestTimes[offset + i] = (short) values[PERCENTILE_BEST][i];
|
||
}
|
||
if (values[PERCENTILE_MEDIAN][i] != Integer.MAX_VALUE) {
|
||
medianTimes[offset + i] = (short) values[PERCENTILE_MEDIAN][i];
|
||
}
|
||
}
|
||
|
||
if (enablePaths && result.paths != null) {
|
||
extractPathsIntoOffset(result.paths, tileLen, offset, network.transitLayer, journeys);
|
||
}
|
||
|
||
offset += tileLen;
|
||
}
|
||
|
||
return new FilteredResult(outIndices, medianTimes, bestTimes, journeys);
|
||
}
|
||
|
||
/**
|
||
* Filter destination indices to those within a bounding box of maxRadiusKm from origin.
|
||
* Uses degree-based approximation. Slightly overestimates at corners, which is fine.
|
||
* Backed by STRtree: O(log n + k) per query instead of O(n) scan.
|
||
*/
|
||
@SuppressWarnings("unchecked")
|
||
private static int[] filterByDistance(
|
||
STRtree postcodeIndex,
|
||
double originLat, double originLon,
|
||
double maxRadiusKm) {
|
||
|
||
double degLat = maxRadiusKm / 111.0;
|
||
double degLon = maxRadiusKm / (111.0 * Math.cos(Math.toRadians(originLat)));
|
||
|
||
Envelope queryEnv = new Envelope(
|
||
originLon - degLon, originLon + degLon,
|
||
originLat - degLat, originLat + degLat);
|
||
|
||
List<Integer> hits = postcodeIndex.query(queryEnv);
|
||
int[] result = new int[hits.size()];
|
||
for (int i = 0; i < hits.size(); i++) result[i] = hits.get(i);
|
||
return result;
|
||
}
|
||
|
||
/**
|
||
* Split destinations into geographic chunks that each fit within R5's grid cell limit
|
||
* and optionally a maximum destination count (required for path recording).
|
||
* Sorts by latitude and splits into bands.
|
||
*/
|
||
private static List<DestinationChunk> buildDestinationChunks(
|
||
double[] lats, double[] lons, int maxDestsPerChunk) {
|
||
int n = lats.length;
|
||
|
||
// Sort indices by latitude for geographic chunking: primitive long sort to
|
||
// avoid Integer[] autoboxing per origin (millions of Integer allocs at scale).
|
||
// Pack: high 32 bits = lat as sortable int, low 32 bits = original index.
|
||
int[] sorted = sortIndicesByLat(lats);
|
||
|
||
// Determine grid width (longitude span is the same for all chunks)
|
||
double minLon = Double.MAX_VALUE, maxLon = -Double.MAX_VALUE;
|
||
for (double lon : lons) {
|
||
minLon = Math.min(minLon, lon);
|
||
maxLon = Math.max(maxLon, lon);
|
||
}
|
||
int totalPixels = 256 << ZOOM;
|
||
int gridWidth = lonToPixel(maxLon, totalPixels) - lonToPixel(minLon, totalPixels) + 1;
|
||
int maxHeight = MAX_GRID_CELLS / gridWidth;
|
||
|
||
// Greedily build chunks: extend each band until it would exceed maxHeight or maxDestsPerChunk
|
||
List<DestinationChunk> chunks = new ArrayList<>();
|
||
int start = 0;
|
||
while (start < n) {
|
||
int end = start + 1;
|
||
int topPixel = latToPixel(lats[sorted[start]], totalPixels);
|
||
|
||
while (end < n) {
|
||
if (end - start >= maxDestsPerChunk) break;
|
||
int bottomPixel = latToPixel(lats[sorted[end]], totalPixels);
|
||
if (Math.abs(bottomPixel - topPixel) + 1 > maxHeight) break;
|
||
end++;
|
||
}
|
||
|
||
chunks.add(buildChunk(lats, lons, sorted, start, end));
|
||
start = end;
|
||
}
|
||
|
||
return chunks;
|
||
}
|
||
|
||
/**
|
||
* Compute travel times for a non-transit mode (single percentile, no paths).
|
||
* Result is indexed into the per-origin filtered subset (via chunk.originalIndices).
|
||
*/
|
||
private static short[][] computeTravelTimesDirect(
|
||
TransportNetwork network, List<DestinationChunk> chunks,
|
||
double originLat, double originLon, String mode, int nDest, LocalDate date) {
|
||
|
||
short[][] allTimes = new short[1][nDest];
|
||
Arrays.fill(allTimes[0], (short) -1);
|
||
|
||
for (DestinationChunk chunk : chunks) {
|
||
RegionalTask task = buildTask(chunk, originLat, originLon, mode, date, false);
|
||
TravelTimeComputer computer = new TravelTimeComputer(task, network);
|
||
OneOriginResult result = computer.computeTravelTimes();
|
||
|
||
TravelTimeResult tt = result.travelTimes;
|
||
if (tt == null) {
|
||
throw new RuntimeException("R5 returned null travelTimes for chunk with "
|
||
+ chunk.originalIndices.length + " destinations");
|
||
}
|
||
int[][] values = tt.getValues();
|
||
if (values.length < 1) {
|
||
throw new RuntimeException("R5 returned 0 percentiles, expected 1");
|
||
}
|
||
for (int i = 0; i < chunk.originalIndices.length; i++) {
|
||
if (values[0][i] != Integer.MAX_VALUE) {
|
||
allTimes[0][chunk.originalIndices[i]] = (short) values[0][i];
|
||
}
|
||
}
|
||
}
|
||
return allTimes;
|
||
}
|
||
|
||
/**
|
||
* Extract the most common journey pattern for each destination in a tile,
|
||
* writing into a combined output array at the given offset.
|
||
*/
|
||
@SuppressWarnings("unchecked")
|
||
private static void extractPathsIntoOffset(
|
||
PathResult paths, int tileLen, int offset, TransitLayer transitLayer,
|
||
String[] journeysOut) {
|
||
Multimap<RouteSequence, PathResult.Iteration>[] allPaths = paths.iterationsForPathTemplates;
|
||
int n = Math.min(tileLen, allPaths.length);
|
||
for (int i = 0; i < n; i++) {
|
||
Multimap<RouteSequence, PathResult.Iteration> destPaths = allPaths[i];
|
||
if (destPaths == null || destPaths.isEmpty()) continue;
|
||
|
||
// Find the RouteSequence used by the most departure-time iterations
|
||
RouteSequence bestRoute = null;
|
||
int maxCount = 0;
|
||
for (RouteSequence rs : destPaths.keySet()) {
|
||
int count = destPaths.get(rs).size();
|
||
if (count > maxCount) {
|
||
maxCount = count;
|
||
bestRoute = rs;
|
||
}
|
||
}
|
||
if (bestRoute == null) continue;
|
||
|
||
journeysOut[offset + i] = buildJourneyJson(bestRoute, transitLayer);
|
||
}
|
||
}
|
||
|
||
/** Build a JSON array of journey legs from a RouteSequence. */
|
||
private static String buildJourneyJson(RouteSequence routeSequence, TransitLayer transitLayer) {
|
||
StringBuilder sb = new StringBuilder("[");
|
||
boolean first = true;
|
||
|
||
// Access leg (walk/bike to first transit stop)
|
||
StreetTimesAndModes.StreetTimeAndMode access = routeSequence.stopSequence.access;
|
||
if (access != null && access.time > 0) {
|
||
int mins = (access.time + 30) / 60; // round to nearest minute
|
||
sb.append("{\"mode\":\"").append(access.mode.name().toLowerCase())
|
||
.append("\",\"minutes\":").append(mins).append("}");
|
||
first = false;
|
||
}
|
||
|
||
// Transit legs
|
||
for (RouteSequence.TransitLeg leg : routeSequence.transitLegs(transitLayer)) {
|
||
if (!first) sb.append(",");
|
||
sb.append("{\"mode\":\"").append(escapeJson(leg.route))
|
||
.append("\",\"from\":\"").append(escapeJson(leg.board))
|
||
.append("\",\"to\":\"").append(escapeJson(leg.alight))
|
||
.append("\",\"minutes\":").append(Math.round(leg.inVehicleTime))
|
||
.append("}");
|
||
first = false;
|
||
}
|
||
|
||
// Egress leg (walk/bike from last transit stop)
|
||
StreetTimesAndModes.StreetTimeAndMode egress = routeSequence.stopSequence.egress;
|
||
if (egress != null && egress.time > 0) {
|
||
if (!first) sb.append(",");
|
||
int mins = (egress.time + 30) / 60;
|
||
sb.append("{\"mode\":\"").append(egress.mode.name().toLowerCase())
|
||
.append("\",\"minutes\":").append(mins).append("}");
|
||
}
|
||
|
||
sb.append("]");
|
||
return sb.toString();
|
||
}
|
||
|
||
private static String escapeJson(String s) {
|
||
if (s == null) return "";
|
||
return s.replace("\\", "\\\\").replace("\"", "\\\"");
|
||
}
|
||
|
||
// --- Private helpers ---
|
||
|
||
private record DestinationChunk(FreeFormPointSet pointSet, WebMercatorExtents extents, int[] originalIndices) {}
|
||
|
||
/**
|
||
* Sort destination indices by latitude using a primitive long-packed sort.
|
||
* Encodes lat as a fixed-point microdeg int (+offset to keep it non-negative
|
||
* for any plausible lat) so high 32 bits of the packed long give a monotonic
|
||
* sort key. Low 32 bits hold the original index, breaking ties deterministically.
|
||
*/
|
||
private static int[] sortIndicesByLat(double[] lats) {
|
||
int n = lats.length;
|
||
long[] packed = new long[n];
|
||
// Offset by 90° so any lat in [-90, 90] maps to a non-negative key
|
||
long offset = 900_000_000L;
|
||
for (int i = 0; i < n; i++) {
|
||
long latKey = (long) Math.round(lats[i] * 10_000_000L) + offset;
|
||
packed[i] = (latKey << 32) | (i & 0xFFFFFFFFL);
|
||
}
|
||
Arrays.sort(packed);
|
||
int[] sorted = new int[n];
|
||
for (int i = 0; i < n; i++) sorted[i] = (int) (packed[i] & 0xFFFFFFFFL);
|
||
return sorted;
|
||
}
|
||
|
||
private static DestinationChunk buildChunk(
|
||
double[] lats, double[] lons, int[] sorted, int start, int end) {
|
||
int size = end - start;
|
||
int[] originalIndices = new int[size];
|
||
Coordinate[] coords = new Coordinate[size];
|
||
double minLat = Double.MAX_VALUE, maxLat = -Double.MAX_VALUE;
|
||
double minLon = Double.MAX_VALUE, maxLon = -Double.MAX_VALUE;
|
||
|
||
for (int i = 0; i < size; i++) {
|
||
int idx = sorted[start + i];
|
||
originalIndices[i] = idx;
|
||
double lat = lats[idx], lon = lons[idx];
|
||
coords[i] = new Coordinate(lon, lat); // x=lon, y=lat
|
||
minLat = Math.min(minLat, lat);
|
||
maxLat = Math.max(maxLat, lat);
|
||
minLon = Math.min(minLon, lon);
|
||
maxLon = Math.max(maxLon, lon);
|
||
}
|
||
|
||
FreeFormPointSet pointSet = new FreeFormPointSet(coords);
|
||
int totalPixels = 256 << ZOOM;
|
||
int west = lonToPixel(minLon, totalPixels);
|
||
int north = latToPixel(maxLat, totalPixels);
|
||
int width = lonToPixel(maxLon, totalPixels) - west + 1;
|
||
int height = latToPixel(minLat, totalPixels) - north + 1;
|
||
WebMercatorExtents extents = new WebMercatorExtents(west, north, width, height, ZOOM);
|
||
|
||
return new DestinationChunk(pointSet, extents, originalIndices);
|
||
}
|
||
|
||
/** Like {@link #buildChunk} but produces a {@link PostcodeTile} with bbox + global indices. */
|
||
private static PostcodeTile buildTile(
|
||
double[] lats, double[] lons, int[] sorted, int start, int end) {
|
||
int size = end - start;
|
||
int[] originalIndices = new int[size];
|
||
Coordinate[] coords = new Coordinate[size];
|
||
double minLat = Double.MAX_VALUE, maxLat = -Double.MAX_VALUE;
|
||
double minLon = Double.MAX_VALUE, maxLon = -Double.MAX_VALUE;
|
||
|
||
for (int i = 0; i < size; i++) {
|
||
int idx = sorted[start + i];
|
||
originalIndices[i] = idx;
|
||
double lat = lats[idx], lon = lons[idx];
|
||
coords[i] = new Coordinate(lon, lat);
|
||
minLat = Math.min(minLat, lat);
|
||
maxLat = Math.max(maxLat, lat);
|
||
minLon = Math.min(minLon, lon);
|
||
maxLon = Math.max(maxLon, lon);
|
||
}
|
||
|
||
FreeFormPointSet pointSet = new FreeFormPointSet(coords);
|
||
int totalPixels = 256 << ZOOM;
|
||
int west = lonToPixel(minLon, totalPixels);
|
||
int north = latToPixel(maxLat, totalPixels);
|
||
int width = lonToPixel(maxLon, totalPixels) - west + 1;
|
||
int height = latToPixel(minLat, totalPixels) - north + 1;
|
||
WebMercatorExtents extents = new WebMercatorExtents(west, north, width, height, ZOOM);
|
||
|
||
return new PostcodeTile(pointSet, extents, originalIndices, minLat, maxLat, minLon, maxLon);
|
||
}
|
||
|
||
/** Build a transit RegionalTask that targets one global tile, configured by {@code mode}. */
|
||
private static RegionalTask buildTaskForTile(
|
||
PostcodeTile tile, double originLat, double originLon, String mode, LocalDate date, boolean recordPaths) {
|
||
RegionalTask task = new RegionalTask();
|
||
task.fromLat = originLat;
|
||
task.fromLon = originLon;
|
||
task.date = date;
|
||
task.percentiles = new int[]{5, 50};
|
||
task.recordTimes = true;
|
||
task.destinationPointSets = new PointSet[]{tile.pointSet()};
|
||
task.zoom = tile.extents().zoom;
|
||
task.west = tile.extents().west;
|
||
task.north = tile.extents().north;
|
||
task.width = tile.extents().width;
|
||
task.height = tile.extents().height;
|
||
task.fromTime = DEPARTURE_FROM_TIME;
|
||
task.toTime = DEPARTURE_TO_TIME;
|
||
task.maxTripDurationMinutes = MAX_TRIP_DURATION_MINUTES;
|
||
// TfL GTFS uses frequency-based service patterns. With the default
|
||
// monteCarloDraws=220 R5 runs 8 iters/min (~240 iters per 30-min window).
|
||
// Set to 0 to use HALF_HEADWAY mode → 1 iter/min, deterministic, 8x cheaper.
|
||
task.monteCarloDraws = 0;
|
||
|
||
if (recordPaths) {
|
||
task.includePathResults = true;
|
||
task.nPathsPerTarget = 1;
|
||
}
|
||
|
||
configureMode(task, mode);
|
||
return task;
|
||
}
|
||
|
||
private static RegionalTask buildTask(
|
||
DestinationChunk chunk, double originLat, double originLon, String mode, LocalDate date,
|
||
boolean recordPaths) {
|
||
RegionalTask task = new RegionalTask();
|
||
task.fromLat = originLat;
|
||
task.fromLon = originLon;
|
||
task.date = date;
|
||
task.percentiles = mode.equals("transit") ? new int[]{5, 50} : new int[]{50};
|
||
task.recordTimes = true;
|
||
task.destinationPointSets = new PointSet[]{chunk.pointSet};
|
||
task.zoom = chunk.extents.zoom;
|
||
task.west = chunk.extents.west;
|
||
task.north = chunk.extents.north;
|
||
task.width = chunk.extents.width;
|
||
task.height = chunk.extents.height;
|
||
task.fromTime = DEPARTURE_FROM_TIME;
|
||
task.toTime = DEPARTURE_TO_TIME;
|
||
task.maxTripDurationMinutes = MAX_TRIP_DURATION_MINUTES;
|
||
|
||
if (recordPaths) {
|
||
task.includePathResults = true;
|
||
// We only use the most common RouteSequence (see extractPaths), so 1 path
|
||
// per target is sufficient and cuts path memory by ~67% vs 3.
|
||
task.nPathsPerTarget = 1;
|
||
}
|
||
|
||
configureMode(task, mode);
|
||
return task;
|
||
}
|
||
|
||
private static void configureMode(RegionalTask task, String mode) {
|
||
if (isTransitMode(mode)) {
|
||
TransitConfig config = transitConfigFor(mode);
|
||
task.maxRides = config.maxRides();
|
||
task.maxWalkTime = TRANSIT_MAX_WALK_TIME_MIN;
|
||
// R5 requires directModes ⊆ accessModes. BICYCLE egress is too expensive
|
||
// (builds cost tables from 59k stops × N destinations), so keep WALK only
|
||
// for egress and match access/direct to avoid the R5 validation error.
|
||
task.accessModes = EnumSet.of(LegMode.WALK);
|
||
task.egressModes = EnumSet.of(LegMode.WALK);
|
||
task.directModes = EnumSet.of(LegMode.WALK);
|
||
EnumSet<TransitModes> transitModes = EnumSet.allOf(TransitModes.class);
|
||
if (config.excludeBus()) transitModes.remove(TransitModes.BUS);
|
||
task.transitModes = transitModes;
|
||
return;
|
||
}
|
||
switch (mode) {
|
||
case "car" -> setDirectMode(task, LegMode.CAR);
|
||
case "bicycle" -> setDirectMode(task, LegMode.BICYCLE);
|
||
case "walking" -> setDirectMode(task, LegMode.WALK);
|
||
default -> throw new IllegalArgumentException("Unknown mode: " + mode);
|
||
}
|
||
}
|
||
|
||
private static void setDirectMode(RegionalTask task, LegMode legMode) {
|
||
task.maxRides = 0;
|
||
task.accessModes = EnumSet.of(legMode);
|
||
task.egressModes = EnumSet.of(legMode);
|
||
task.directModes = EnumSet.of(legMode);
|
||
task.transitModes = EnumSet.noneOf(TransitModes.class);
|
||
}
|
||
|
||
private static int lonToPixel(double lon, int totalPixels) {
|
||
return (int) Math.floor(totalPixels * (lon + 180.0) / 360.0);
|
||
}
|
||
|
||
private static int latToPixel(double lat, int totalPixels) {
|
||
double latRad = Math.toRadians(lat);
|
||
return (int) Math.floor(totalPixels * (1.0 - Math.log(Math.tan(latRad) + 1.0 / Math.cos(latRad)) / Math.PI) / 2.0);
|
||
}
|
||
}
|