improve diverse plan generator and models

input/v6.3/params/bg_diverse_9_fixed.yaml

@@ -0,0 +1,42 @@
+#  Results for model plan-choices-diverse_9_fixed_ascs
+#  Nbr of parameters:		4
+#  Sample size:			8259
+#  Excluded data:			0
+#  Null log likelihood:		-16584.31
+#  Final log likelihood:		-18747.14
+#  Likelihood ratio test (null):		-4325.659
+#  Rho square (null):			-0.13
+#  Rho bar square (null):			-0.131
+#  Akaike Information Criterion:	37502.28
+#  Bayesian Information Criterion:	37530.35
+scoring:
+  scoringParameters:
+    - modeParams:
+        - mode: walk
+          constant: 0
+        - mode: car
+          constant: -0.5341414592094356
+          dailyMonetaryConstant: -14.30
+          dailyUtilityConstant: 6.078136
+        - mode: pt
+          constant: 0.3971116
+        - mode: bike
+          constant: -1.3538876325
+        - mode: ride
+          constant: -1.23976957093642
+advancedScoring:
+  scoringParameters:
+    - subpopulation: person
+      modeParams:
+        - mode: car
+          deltaDailyConstant: 13.368402
+          varDailyConstant: truncatedNormal
+        - mode: bike
+          deltaConstant: 2.040308
+          varConstant: normal
+        - mode: pt
+          deltaConstant: 3.609201
+          varConstant: normal
+        - mode: ride
+          deltaConstant: 0.968879
+          varConstant: normal

input/v6.3/params/mxl_diverse_9.yaml

@@ -0,0 +1,34 @@
+#  Log-Likelihood= -14744.630
+#    AIC= 29505.259
+#    BIC= 29561.412
+scoring:
+  scoringParameters:
+    - modeParams:
+        - mode: walk
+          constant: 0
+        - mode: car
+          constant: 0
+          dailyMonetaryConstant: -14.30
+          dailyUtilityConstant: 8.9634862
+        - mode: pt
+          constant: 0.2118402
+        - mode: bike
+          constant: -1.4056730
+        - mode: ride
+          constant: -2.0662459
+advancedScoring:
+  scoringParameters:
+    - subpopulation: person
+      modeParams:
+        - mode: car
+          deltaDailyConstant: 5.2291379
+          varDailyConstant: truncatedNormal
+        - mode: bike
+          deltaConstant: 0.6887108
+          varConstant: normal
+        - mode: pt
+          deltaConstant: 1.7932632
+          varConstant: normal
+        - mode: ride
+          deltaConstant: 1.2994295
+          varConstant: normal

src/main/java/org/matsim/prepare/choices/ComputePlanChoices.java

@@ -181,7 +181,7 @@ public Integer call() throws Exception {
 
 		String out = output.toString().replace(".csv", "-%s_%d.csv".formatted(planCandidates, topK));
 
-		if (timeUtil && planCandidates == PlanCandidates.bestK) {
+		if (timeUtil && (planCandidates == PlanCandidates.bestK || planCandidates == PlanCandidates.diverse)) {
 			out = out.replace(".csv", "-tt-only.csv");
 		}
 

src/main/java/org/matsim/prepare/choices/DiversePlanGenerator.java

@@ -16,6 +16,7 @@ public class DiversePlanGenerator implements CandidateGenerator {
 
 	private final int topK;
 	private final TopKChoicesGenerator gen;
+	private final SplittableRandom rnd = new SplittableRandom(0);
 
 	DiversePlanGenerator(int topK, TopKChoicesGenerator generator) {
 		this.topK = topK;
@@ -33,20 +34,42 @@ public List<PlanCandidate> generate(PlanModel planModel, @Nullable Set<String> c
 
 		List<PlanCandidate> candidates = new ArrayList<>();
 		boolean carUser = PersonUtils.canUseCar(planModel.getPerson());
-
 		Set<String> modes = new HashSet<>(consideredModes);
-		modes.remove(carUser ? "ride": "car");
-
 
 		for (String mode : modes) {
+			if (!carUser && mode.equals("car"))
+				continue;
+
 			List<PlanCandidate> tmp = gen.generate(planModel, Set.of(mode), mask);
 			if (!tmp.isEmpty())
 				candidates.add(tmp.get(0));
 		}
 
-		Collections.sort(candidates);
-		candidates.add(0, existing);
+		candidates.addFirst(existing);
+
+		// Add combination of modes as well
+		addToCandidates(candidates, gen.generate(planModel, modes, mask), carUser ? "car" : "ride", 1);
+		addToCandidates(candidates, gen.generate(planModel, consideredModes, mask), null, 1);
 
-		return candidates.stream().distinct().limit(topK).toList();
+		// Remove the primary mode to generate remaining alternatives
+		modes.remove(carUser ? "car" : "ride");
+		addToCandidates(candidates, gen.generate(planModel, modes, mask), null, 2);
+
+		return candidates.stream().distinct().limit(this.topK).toList();
 	}
+
+	private void addToCandidates(List<PlanCandidate> candidates, List<PlanCandidate> topK, String requireMode, int n) {
+
+		topK.removeIf(candidates::contains);
+
+		if (requireMode != null) {
+			topK.removeIf(c -> !c.containsMode(requireMode));
+		}
+
+		for (int i = 0; i < n; i++) {
+			if (topK.size() > 1)
+				candidates.add(topK.remove(rnd.nextInt(topK.size())));
+		}
+	}
+
 }

src/main/java/org/matsim/prepare/choices/RandomPlanGenerator.java

@@ -33,6 +33,8 @@ public List<PlanCandidate> generate(PlanModel planModel, @Nullable Set<String> c
 		PlanCandidate existing = gen.generatePredefined(planModel, chosen).get(0);
 
 		// This changes the internal state to randomize the estimates
+		// TODO random selection is biased because of mass conservation
+		// due to that, this class should not be used
 		for (Map.Entry<String, List<ModeEstimate>> entry : planModel.getEstimates().entrySet()) {
 			for (ModeEstimate est : entry.getValue()) {
 				double[] utils = est.getEstimates();

src/main/java/org/matsim/prepare/network/ModifyNetwork.java

@@ -0,0 +1,134 @@
+package org.matsim.prepare.network;
+
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+import org.geotools.api.feature.simple.SimpleFeature;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.LineString;
+import org.locationtech.jts.geom.Point;
+import org.locationtech.jts.index.strtree.STRtree;
+import org.matsim.api.core.v01.Id;
+import org.matsim.api.core.v01.network.Link;
+import org.matsim.api.core.v01.network.Network;
+import org.matsim.api.core.v01.network.Node;
+import org.matsim.application.MATSimAppCommand;
+import org.matsim.application.options.ShpOptions;
+import org.matsim.core.network.NetworkUtils;
+import org.matsim.core.utils.geometry.geotools.MGC;
+import picocli.CommandLine;
+
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.util.Comparator;
+import java.util.List;
+import java.util.function.ToDoubleFunction;
+
+@CommandLine.Command(name = "modify-network", description = "Remove or add network elements")
+public class ModifyNetwork implements MATSimAppCommand {
+
+	private static final Logger log = LogManager.getLogger(ModifyNetwork.class);
+
+	@CommandLine.Option(names = {"--network"}, description = "Path to the network file", required = true)
+	private String networkPath;
+
+	@CommandLine.Option(names = {"--remove-links"}, description = "Path to the CSV file with link ids to remove")
+	private Path removeCSV;
+
+	@CommandLine.Option(names = {"--matching-distance"}, description = "Distance in meters to match links from shapefile to network nodes", defaultValue = "10.0")
+	private double matchingDistance;
+
+	@CommandLine.Option(names = "--output", description = "Path to the output network file", required = true)
+	private Path output;
+
+	@CommandLine.Mixin
+	private ShpOptions shp = new ShpOptions();
+
+	public static void main(String[] args) {
+		new ModifyNetwork().execute(args);
+	}
+
+	@Override
+	public Integer call() throws Exception {
+
+		Network network = NetworkUtils.readNetwork(networkPath);
+
+		if (removeCSV != null) {
+			removeLinks(network);
+		}
+
+		if (shp.isDefined()) {
+			addLinksFromShape(network);
+		}
+
+		return 0;
+	}
+
+	private void removeLinks(Network network) throws IOException {
+		for (String line : Files.readAllLines(removeCSV)) {
+			Id<Link> linkId = Id.createLinkId(line);
+			network.removeLink(linkId);
+		}
+
+		// Remove links that are not connected
+		List<Id<Node>> toRemove = network.getNodes().values().stream().filter(n -> n.getInLinks().isEmpty() && n.getOutLinks().isEmpty())
+			.map(Node::getId)
+			.toList();
+
+		toRemove.forEach(network::removeNode);
+	}
+
+	/**
+	 * Read line geometry from a shapefile and add links to the network.
+	 */
+	private void addLinksFromShape(Network network) {
+
+		STRtree nodeIndex = new STRtree();
+		for (Node node : network.getNodes().values()) {
+			nodeIndex.insert(MGC.coord2Point(node.getCoord()).getEnvelopeInternal(), node);
+		}
+		nodeIndex.build();
+
+		List<SimpleFeature> features = shp.readFeatures();
+		for (SimpleFeature feature : features) {
+			// Find the nearest node
+			LineString geom = (LineString) feature.getDefaultGeometry();
+
+			Id<Link> linkId = Id.createLinkId(feature.getID());
+			Node fromNode = matchNode(nodeIndex, geom.getCoordinateN(0));
+			Node toNode = matchNode(nodeIndex, geom.getCoordinateN(geom.getNumPoints() - 1));
+
+			if (fromNode == toNode) {
+				throw new IllegalArgumentException("Link %s would have the same start and end node. Check the shapefile and improve matching to zones.".formatted(linkId));
+			}
+
+			Link link = network.getFactory().createLink(linkId, fromNode, toNode);
+
+			// TODO: Set link attributes
+
+			network.addLink(link);
+		}
+	}
+
+	@SuppressWarnings("unchecked")
+	private Node matchNode(STRtree index, Coordinate coord) {
+
+		Point point = MGC.coordinate2Point(coord);
+
+		ToDoubleFunction<Object> distance = n -> NetworkUtils.getEuclideanDistance(((Node) n).getCoord(), MGC.coordinate2Coord(coord));
+
+		List<Node> result = index.query(point.buffer(matchingDistance).getEnvelopeInternal())
+			.stream()
+			.map(Node.class::cast)
+			.filter(n -> distance.applyAsDouble(n) < matchingDistance)
+			.sorted(Comparator.comparingDouble(distance))
+			.toList();
+
+		if (result.isEmpty()) {
+			throw new IllegalArgumentException("No node found for coordinate %s".formatted(coord));
+		}
+
+		return result.getFirst();
+	}
+
+}