MATSim Application: Improve network params estimators

contribs/application/src/main/java/org/matsim/application/analysis/traffic/CountComparisonAnalysis.java

@@ -76,6 +76,16 @@ private static int[] sum(int[] a, int[] b) {
 		return counts;
 	}
 
+	/**
+	 * Calculate the geh value for simulated and reference count
+	 */
+	private static double geh(double simulated, double observed) {
+		final double diff = simulated - observed;
+		final double sum = simulated + observed;
+
+		return Math.sqrt(2 * diff * diff / sum);
+	}
+
 	@Override
 	public Integer call() throws Exception {
 
@@ -114,14 +124,16 @@ private Table writeOutput(Counts<Link> counts, Network network, VolumesAnalyzer
 			StringColumn.create("road_type"),
 			IntColumn.create("hour"),
 			DoubleColumn.create("observed_traffic_volume"),
-			DoubleColumn.create("simulated_traffic_volume")
+			DoubleColumn.create("simulated_traffic_volume"),
+			DoubleColumn.create("geh")
 		);
 
 		Table dailyTrafficVolume = Table.create(StringColumn.create("link_id"),
 			StringColumn.create("name"),
 			StringColumn.create("road_type"),
 			DoubleColumn.create("observed_traffic_volume"),
-			DoubleColumn.create("simulated_traffic_volume")
+			DoubleColumn.create("simulated_traffic_volume"),
+			DoubleColumn.create("geh")
 		);
 
 		for (Map.Entry<Id<Link>, MeasurementLocation<Link>> entry : counts.getMeasureLocations().entrySet()) {
@@ -170,6 +182,7 @@ private Table writeOutput(Counts<Link> counts, Network network, VolumesAnalyzer
 					row.setInt("hour", hour);
 					row.setDouble("observed_traffic_volume", observedTrafficVolumeAtHour);
 					row.setDouble("simulated_traffic_volume", simulatedTrafficVolumeAtHour);
+					row.setDouble("geh", geh(simulatedTrafficVolumeAtHour, observedTrafficVolumeAtHour));
 				}
 			} else {
 				// Get the daily values
@@ -183,6 +196,7 @@ private Table writeOutput(Counts<Link> counts, Network network, VolumesAnalyzer
 			row.setString("road_type", type);
 			row.setDouble("observed_traffic_volume", observedTrafficVolumeByDay);
 			row.setDouble("simulated_traffic_volume", simulatedTrafficVolumeByDay);
+			row.setDouble("geh", geh(simulatedTrafficVolumeByDay, observedTrafficVolumeByDay));
 		}
 
 		DoubleColumn relError = dailyTrafficVolume.doubleColumn("simulated_traffic_volume")

contribs/application/src/main/java/org/matsim/application/prepare/Predictor.java

@@ -11,6 +11,7 @@ public interface Predictor {
 
 	/**
 	 * Predict value from given features.
+	 * @return predicted value, maybe NaN if no prediction is possible.
 	 */
 	double predict(Object2DoubleMap<String> features, Object2ObjectMap<String, String> categories);
 

contribs/application/src/main/java/org/matsim/application/prepare/network/SampleNetwork.java

@@ -144,6 +144,10 @@ public Integer call() throws Exception {
 
 				List<? extends Node> list = e.getValue();
 
+				// for now only consider traffic lights
+				if (!e.getKey().equals("traffic_light"))
+					continue;
+
 				log.info("Sampling {} out of {} intersections for type {}", sample, list.size(), e.getKey());
 
 				for (int i = 0; i < sample && !list.isEmpty(); i++) {
@@ -193,7 +197,7 @@ public Integer call() throws Exception {
 
 		RandomizedTravelTime tt = new RandomizedTravelTime(rnd);
 
-		LeastCostPathCalculator router = createRandomizedRouter(network, tt);
+		LeastCostPathCalculator router = createRandomizedRouter(cityNetwork, tt);
 
 		sampleCityRoutes(cityNetwork, router, tt, rnd);
 
@@ -223,9 +227,16 @@ private void sampleCityRoutes(Network network, LeastCostPathCalculator router, R
 
 				Link to = NetworkUtils.getNearestLink(network, dest);
 
+				// Links could be on the very edge so that nodes are outside the network
+				if (to == null || !network.getNodes().containsKey(link.getFromNode().getId()) ||
+					!network.getNodes().containsKey(to.getToNode().getId())) {
+					i--;
+					continue;
+				}
+
 				LeastCostPathCalculator.Path path = router.calcLeastCostPath(link.getFromNode(), to.getToNode(), 0, null, null);
 
-				if (path.nodes.size() < 2) {
+				if (path == null || path.nodes.size() < 2) {
 					i--;
 					continue;
 				}

contribs/application/src/main/java/org/matsim/application/prepare/network/params/ApplyNetworkParams.java

@@ -4,7 +4,6 @@
 import com.fasterxml.jackson.annotation.JsonInclude;
 import com.fasterxml.jackson.annotation.PropertyAccessor;
 import com.fasterxml.jackson.databind.ObjectMapper;
-import it.unimi.dsi.fastutil.objects.Object2DoubleMap;
 import org.apache.logging.log4j.LogManager;
 import org.apache.logging.log4j.Logger;
 import org.matsim.api.core.v01.Id;
@@ -14,12 +13,16 @@
 import org.matsim.application.MATSimAppCommand;
 import org.matsim.application.options.InputOptions;
 import org.matsim.application.options.OutputOptions;
+import org.matsim.application.prepare.Predictor;
 import org.matsim.application.prepare.network.params.NetworkParamsOpt.Feature;
 import org.matsim.core.network.NetworkUtils;
 import org.matsim.core.utils.io.IOUtils;
 import picocli.CommandLine;
 
 import java.io.BufferedReader;
+import java.math.BigDecimal;
+import java.math.RoundingMode;
+import java.util.List;
 import java.util.Map;
 import java.util.Set;
 
@@ -40,6 +43,7 @@ public class ApplyNetworkParams implements MATSimAppCommand {
 	@CommandLine.Mixin
 	private final OutputOptions output = OutputOptions.ofCommand(ApplyNetworkParams.class);
 
+
 	@CommandLine.Parameters(arity = "1..*", description = "Type of parameters to apply. Available: ${COMPLETION-CANDIDATES}")
 	private Set<NetworkAttribute> params;
 
@@ -52,6 +56,19 @@ public class ApplyNetworkParams implements MATSimAppCommand {
 	@CommandLine.Option(names = "--factor-bounds", split = ",", description = "Speed factor limits (lower,upper bound)", defaultValue = NetworkParamsOpt.DEFAULT_FACTOR_BOUNDS)
 	private double[] speedFactorBounds;
 
+	@CommandLine.Option(names = "--capacity-bounds", split = ",", defaultValue = "0.4,0.6,0.8",
+		description = "Minimum relative capacity against theoretical max (traffic light,right before left, priority)")
+	private List<Double> capacityBounds;
+
+	@CommandLine.Option(names = "--road-types", split = ",", description = "Road types to apply changes to")
+	private Set<String> roadTypes;
+
+	@CommandLine.Option(names = "--junction-types", split = ",", description = "Junction types to apply changes to")
+	private Set<String> junctionTypes;
+
+	@CommandLine.Option(names = "--decrease-only", description = "Only set values if the are lower than the current value", defaultValue = "false")
+	private boolean decrease;
+
 	private NetworkModel model;
 	private NetworkParams paramsOpt;
 
@@ -62,14 +79,14 @@ public static void main(String[] args) {
 	}
 
 	/**
-	 * Theoretical capacity.
+	 * Theoretical capacity assuming fixed car length and headway. This should usually not be exceeded.
 	 */
-	private static double capacityEstimate(double v) {
+	public static double capacityEstimate(double v) {
 
 		// headway
 		double tT = 1.2;
 
-		// car length
+		// car length + buffer
 		double lL = 7.0;
 
 		double Qc = v / (v * tT + lL);
@@ -94,13 +111,19 @@ public Integer call() throws Exception {
 			}
 		}
 
-		Map<Id<Link>, Feature> features = NetworkParamsOpt.readFeatures(input.getPath("features.csv"), network.getLinks().size());
+		Map<Id<Link>, Feature> features = NetworkParamsOpt.readFeatures(input.getPath("features.csv"), network.getLinks());
 
 		for (Link link : network.getLinks().values()) {
 			Feature ft = features.get(link.getId());
 
+			if (roadTypes != null && !roadTypes.isEmpty() && !roadTypes.contains(ft.highwayType()))
+				continue;
+
+			if (junctionTypes != null && !junctionTypes.isEmpty() && !junctionTypes.contains(ft.junctionType()))
+				continue;
+
 			try {
-				applyChanges(link, ft.junctionType(), ft.features());
+				applyChanges(link, ft);
 			} catch (IllegalArgumentException e) {
 				warn++;
 				log.warn("Error processing link {}", link.getId(), e);
@@ -117,67 +140,126 @@ public Integer call() throws Exception {
 	/**
 	 * Apply speed and capacity models and apply changes.
 	 */
-	private void applyChanges(Link link, String junctionType, Object2DoubleMap<String> features) {
-
-		String type = NetworkUtils.getHighwayType(link);
+	private void applyChanges(Link link, Feature ft) {
 
 		boolean modified = false;
 
 		if (params.contains(NetworkAttribute.capacity)) {
+			modified = applyCapacity(link, ft);
+		}
+
+		if (params.contains(NetworkAttribute.freespeed)) {
+			modified |= applyFreeSpeed(link, ft);
+		}
 
-			FeatureRegressor capacity = model.capacity(junctionType);
+		if (modified)
+			warn++;
+	}
 
-			double perLane = capacity.predict(features);
+	private boolean applyCapacity(Link link, Feature ft) {
 
-			double cap = capacityEstimate(features.getDouble("speed"));
+		Predictor capacity = model.capacity(ft.junctionType(), ft.highwayType());
+		// No operation performed if not supported
+		if (capacity == null) {
+			return false;
+		}
 
-			// Minimum thresholds
-			double threshold = switch (junctionType) {
-				// traffic light can reduce capacity at least to 50% (with equal green split)
-				case "traffic_light" -> 0.4;
-				case "right_before_left" -> 0.6;
-				// Motorways are kept at their max theoretical capacity
-				case "priority" -> type.startsWith("motorway") ? 1 : 0.8;
-				default -> 0;
-			};
+		double perLane = capacity.predict(ft.features(), ft.categories());
+		if (Double.isNaN(perLane)) {
+			return true;
+		}
 
-			if (perLane < cap * threshold) {
-				log.warn("Increasing capacity per lane on {} ({}, {}) from {} to {}", link.getId(), type, junctionType, perLane, cap * threshold);
-				perLane = cap * threshold;
-				modified = true;
-			}
+		double cap = capacityEstimate(ft.features().getDouble("speed"));
 
-			link.setCapacity(link.getNumberOfLanes() * perLane);
+		if (capacityBounds.isEmpty())
+			capacityBounds.add(0.0);
+
+		// Fill up to 3 elements if not provided
+		if (capacityBounds.size() < 3) {
+			capacityBounds.add(capacityBounds.get(0));
+			capacityBounds.add(capacityBounds.get(1));
 		}
 
+		// Minimum thresholds
+		double threshold = switch (ft.junctionType()) {
+			case "traffic_light" -> capacityBounds.get(0);
+			case "right_before_left" -> capacityBounds.get(1);
+			case "priority" -> capacityBounds.get(2);
+			default -> 0;
+		};
 
-		if (params.contains(NetworkAttribute.freespeed)) {
+		boolean modified = false;
 
-			double speedFactor = 1.0;
-			FeatureRegressor speedModel = model.speedFactor(junctionType);
+		if (perLane < cap * threshold) {
+			log.warn("Increasing capacity per lane on {} ({}, {}) from {} to {}",
+				link.getId(), ft.highwayType(), ft.junctionType(), perLane, cap * threshold);
+			perLane = cap * threshold;
+			modified = true;
+		}
 
-			speedFactor =  paramsOpt != null ?
-				speedModel.predict(features, paramsOpt.getParams(junctionType)) :
-				speedModel.predict(features);
+		if (perLane > cap) {
+			log.warn("Reducing capacity per lane on {} ({}, {}) from {} to {}",
+				link.getId(), ft.highwayType(), ft.junctionType(), perLane, cap);
+			perLane = cap;
+			modified = true;
+		}
 
-			if (speedFactor > speedFactorBounds[1]) {
-				log.warn("Reducing speed factor on {} from {} to {}", link.getId(), speedFactor, speedFactorBounds[1]);
-				speedFactor = speedFactorBounds[1];
-				modified = true;
-			}
+		if (ft.features().getOrDefault("num_lanes", link.getNumberOfLanes()) != link.getNumberOfLanes())
+			log.warn("Number of lanes for link {} does not match the feature file", link.getId());
 
-			// Threshold for very low speed factors
-			if (speedFactor < speedFactorBounds[0]) {
-				log.warn("Increasing speed factor on {} from {} to {}", link, speedFactor, speedFactorBounds[0]);
-				speedFactor = speedFactorBounds[0];
-				modified = true;
-			}
+		int totalCap = BigDecimal.valueOf(link.getNumberOfLanes() * perLane).setScale(0, RoundingMode.HALF_UP).intValue();
+
+		if (decrease && totalCap > link.getCapacity())
+			return false;
+
+		link.setCapacity(totalCap);
+
+		return modified;
+	}
+
+	private boolean applyFreeSpeed(Link link, Feature ft) {
 
-			link.setFreespeed((double) link.getAttributes().getAttribute("allowed_speed") * speedFactor);
-			link.getAttributes().putAttribute("speed_factor", speedFactor);
+		Predictor speedModel = model.speedFactor(ft.junctionType(), ft.highwayType());
+
+		// No operation performed if not supported
+		if (speedModel == null) {
+			return false;
 		}
 
-		if (modified)
-			warn++;
+		double speedFactor = paramsOpt != null ?
+			speedModel.predict(ft.features(), ft.categories(), paramsOpt.getParams(ft.junctionType())) :
+			speedModel.predict(ft.features(), ft.categories());
+
+		if (Double.isNaN(speedFactor)) {
+			return false;
+		}
+
+		boolean modified = false;
+
+		if (speedFactor > speedFactorBounds[1]) {
+			log.warn("Reducing speed factor on {} from {} to {}", link.getId(), speedFactor, speedFactorBounds[1]);
+			speedFactor = speedFactorBounds[1];
+			modified = true;
+		}
+
+		// Threshold for very low speed factors
+		if (speedFactor < speedFactorBounds[0]) {
+			log.warn("Increasing speed factor on {} from {} to {}", link, speedFactor, speedFactorBounds[0]);
+			speedFactor = speedFactorBounds[0];
+			modified = true;
+		}
+
+		double freeSpeed = (double) link.getAttributes().getAttribute("allowed_speed") * speedFactor;
+
+		freeSpeed = BigDecimal.valueOf(freeSpeed).setScale(3, RoundingMode.HALF_EVEN).doubleValue();
+
+		if (decrease && freeSpeed > link.getFreespeed())
+			return false;
+
+		link.setFreespeed(freeSpeed);
+		link.getAttributes().putAttribute("speed_factor", freeSpeed);
+
+		return modified;
 	}
+
 }