From f8422826c72d25e185e7344c719dbb0ff49b0c44 Mon Sep 17 00:00:00 2001 From: rakow Date: Wed, 7 Aug 2024 23:05:59 +0200 Subject: [PATCH] assign activity chains --- Makefile | 23 +-- pom.xml | 10 +- .../prepare/population/CreateDailyPlans.java | 191 ++++++++++++++++-- .../population/CreateKansaiPopulation.java | 52 +++-- 4 files changed, 225 insertions(+), 51 deletions(-) diff --git a/Makefile b/Makefile index 62a5262..a2cc27a 100644 --- a/Makefile +++ b/Makefile @@ -92,8 +92,9 @@ input/$V/$N-$V-facilities.xml.gz: input/$V/$N-$V-network.xml.gz input/facilities input/$V/$N-static-$V-10pct.plans.xml.gz: input/facilities.gpkg $(sc) prepare kansai-population\ --input $(kyoto)/data/census_kansai_region.csv\ - --shp $(kyoto)/data/kansai-region-epsg4612.gpkg --shp-crs $(CRS)\ - --facilities $< --facilities-attr resident\ + --shp $(kyoto)/data/kansai-region.gpkg\ + --postal-shp $(kyoto)/data/postalcodes.gpkg\ + --facilities $< --facilities-attr all\ --output $@ @@ -102,27 +103,11 @@ input/$V/$N-$V-10pct.plans-initial.xml.gz: input/$V/$N-static-$V-10pct.plans.xml $(sc) prepare create-daily-plans --input $< --output $@\ --persons src/main/python/table-persons.csv\ --activities src/main/python/table-activities.csv\ - --shp $(kyoto)/data/kansai-region-epsg4612.gpkg --shp-crs $(CRS)\ + --shp $(kyoto)/data/postalcodes.gpkg\ --facilities $(word 2,$^)\ --network $(word 3,$^)\ -# Assigns locations to the activities -$p/berlin-initial-$V-10pct.plans.xml.gz: $p/$N-activities-$V-25pct.plans.xml.gz input/$V/$N-$V-facilities.xml.gz input/$V/$N-$V-network.xml.gz - $(sc) prepare init-location-choice\ - --input $<\ - --output $@\ - --facilities $(word 2,$^)\ - --network $(word 3,$^)\ - --shp $(germany)/vg5000/vg5000_ebenen_0101/VG5000_GEM.shp\ - --commuter $(germany)/regionalstatistik/commuter.csv\ - - # For debugging and visualization - $(sc) prepare downsample-population $@\ - --sample-size 0.25\ - --samples 0.1 0.03 0.01\ - - # Aggregated target for input plans to calibration prepare: input/$V/$N-$V-10pct.plans-initial.xml.gz input/$V/$N-$V-network-with-pt.xml.gz echo "Done" \ No newline at end of file diff --git a/pom.xml b/pom.xml index 8580e4d..0715b63 100644 --- a/pom.xml +++ b/pom.xml @@ -67,7 +67,8 @@ com.github.matsim-scenarios matsim-berlin - ce81ea226c + 6.3-SNAPSHOT + @@ -103,6 +104,13 @@ + + maven-enforcer-plugin + + true + + + org.apache.maven.plugins maven-compiler-plugin diff --git a/src/main/java/org/matsim/prepare/population/CreateDailyPlans.java b/src/main/java/org/matsim/prepare/population/CreateDailyPlans.java index 74d00d5..932418c 100644 --- a/src/main/java/org/matsim/prepare/population/CreateDailyPlans.java +++ b/src/main/java/org/matsim/prepare/population/CreateDailyPlans.java @@ -1,18 +1,28 @@ package org.matsim.prepare.population; +import me.tongfei.progressbar.ProgressBar; +import org.apache.commons.csv.CSVFormat; +import org.apache.commons.csv.CSVParser; +import org.apache.commons.csv.CSVRecord; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; +import org.matsim.api.core.v01.Coord; import org.matsim.api.core.v01.population.Person; import org.matsim.api.core.v01.population.Population; import org.matsim.application.MATSimAppCommand; import org.matsim.application.options.ShpOptions; +import org.matsim.core.population.PersonUtils; import org.matsim.core.population.PopulationUtils; import org.matsim.core.population.algorithms.ParallelPersonAlgorithmUtils; import org.matsim.core.population.algorithms.PersonAlgorithm; import picocli.CommandLine; +import java.math.BigInteger; +import java.nio.charset.StandardCharsets; import java.nio.file.Path; -import java.util.SplittableRandom; +import java.util.*; +import java.util.stream.IntStream; +import java.util.stream.Stream; @CommandLine.Command( name = "create-daily-plans", @@ -21,37 +31,40 @@ public class CreateDailyPlans implements MATSimAppCommand, PersonAlgorithm { private static final Logger log = LogManager.getLogger(CreateDailyPlans.class); - + private final Map persons = new HashMap<>(); + private final Map> groups = new HashMap<>(); @CommandLine.Option(names = "--input", description = "Path to input population.") private Path input; - @CommandLine.Option(names = "--output", description = "Path to output population", required = true) private Path output; - @CommandLine.Option(names = "--persons", description = "Path to person table", required = true) private Path personsPath; - @CommandLine.Option(names = "--activities", description = "Path to activity table", required = true) private Path activityPath; - @CommandLine.Option(names = "--facilities", description = "Path to facilities file", required = true) private Path facilityPath; - @CommandLine.Option(names = "--network", description = "Path to network file", required = true) private Path networkPath; - @CommandLine.Option(names = "--seed", description = "Seed used to sample locations", defaultValue = "1") private long seed; - @CommandLine.Mixin private ShpOptions shp; - - private ThreadLocal rnd; + private Population population; + private Map> activities; + private ProgressBar pb; public static void main(String[] args) { new CreateDailyPlans().execute(args); } + /** + * Initializes random number generator with person specific seed. + */ + private SplittableRandom initRandomNumberGenerator(Person person) { + BigInteger i = new BigInteger(person.getId().toString().getBytes()); + return new SplittableRandom(i.longValue() + seed); + } + @Override public Integer call() throws Exception { @@ -60,9 +73,25 @@ public Integer call() throws Exception { return 2; } - Population population = PopulationUtils.readPopulation(input.toString()); + // TODO: map zones to shape file + + activities = RunActivitySampling.readActivities(activityPath); + + // Remove activities with missing leg duration + activities.values().removeIf( + rows -> rows.stream().anyMatch(act -> act.get("leg_duration").isBlank() || act.get("duration").isBlank()) + ); + + log.info("Got {} persons after cleaning", activities.size()); + + try (CSVParser csv = CSVParser.parse(personsPath, StandardCharsets.UTF_8, + CSVFormat.DEFAULT.builder().setHeader().setSkipHeaderRecord(true).build())) { + readPersons(csv); + } - rnd = ThreadLocal.withInitial(() -> new SplittableRandom(seed)); + population = PopulationUtils.readPopulation(input.toString()); + + pb = new ProgressBar("Creating daily plans", population.getPersons().size()); ParallelPersonAlgorithmUtils.run(population, 8, this); @@ -74,7 +103,141 @@ public Integer call() throws Exception { @Override public void run(Person person) { - // TODO: map zones to shape file + SplittableRandom rnd = initRandomNumberGenerator(person); + + Coord homeCoord = Attributes.getHomeCoord(person); + + String personId = matchPerson(rnd, createKey(person, Zone.FULL)); + if (personId == null) { + personId = matchPerson(rnd, createKey(person, Zone.CITY)); + } + if (personId == null) { + personId = matchPerson(rnd, createKey(person, Zone.NONE)); + } + + if (personId == null) { + log.warn("No matching person found for {}", createKey(person, Zone.NONE)); + return; + } + + List activities = Objects.requireNonNull(this.activities.get(personId), "No activities found"); + + RunActivitySampling.createPlan(homeCoord, activities, rnd, population.getFactory()); + + // TODO: select locations + + // TODO: add reference modes and weights for fully matched persons + + pb.step(); + } + + /** + * Create subpopulations for sampling. + */ + private void readPersons(CSVParser csv) { + + int i = 0; + int skipped = 0; + for (CSVRecord r : csv) { + + String idx = r.get("p_id"); + String gender = r.get("gender"); + int age = Integer.parseInt(r.get("age")); + + if (!activities.containsKey(idx)) { + skipped++; + continue; + } + + Stream keys = createKey(gender, age, r.get("location"), r.get("zone")); + keys.forEach(key -> { + groups.computeIfAbsent(key, (k) -> new ArrayList<>()).add(idx); + + // Alternative keys with different zones + groups.computeIfAbsent(new Key(key.gender, key.age, r.get("location")), (k) -> new ArrayList<>()).add(idx); + groups.computeIfAbsent(new Key(key.gender, key.age, null), (k) -> new ArrayList<>()).add(idx); + + }); + persons.put(idx, r); + i++; + } + + log.info("Read {} persons from csv. Skipped {} invalid persons", i, skipped); } + + /** + * Match person attributes with person from survey data. + * + * @return daily activities + */ + private String matchPerson(SplittableRandom rnd, Key key) { + List subgroup = groups.get(key); + if (subgroup == null) { + return null; + } + + if (subgroup.size() < 5) + return null; + + // TODO: needs to be weighted matching + // weights can be preprocessed after reading in + + return subgroup.get(rnd.nextInt(subgroup.size())); + } + + private Stream createKey(String gender, int age, String location, String zone) { + + String homeZone; + if (zone.isBlank() || zone.equals(location)) + homeZone = location; + else + // The last two digits of the postal code are not known + homeZone = location + "_" + zone.substring(0, zone.length() - 2); + + if (age <= 10) { + return IntStream.rangeClosed(0, 10).mapToObj(i -> new Key(null, i, homeZone)); + } + if (age < 18) { + return IntStream.rangeClosed(11, 18).mapToObj(i -> new Key(gender, i, homeZone)); + } + + int min = Math.max(18, age - 6); + int max = Math.min(65, age + 6); + + // larger groups for older people + if (age > 65) { + min = Math.max(66, age - 10); + max = Math.min(105, age + 10); + } + + return IntStream.rangeClosed(min, max).mapToObj(i -> new Key(gender, i, homeZone)); + } + + private Key createKey(Person person, Zone zone) { + + Integer age = PersonUtils.getAge(person); + String gender = PersonUtils.getSex(person); + if (age <= 10) + gender = null; + + String homeZone = switch (zone) { + case FULL -> person.getAttributes().getAttribute("city") + "_" + person.getAttributes().getAttribute("postal"); + case CITY -> Objects.toString(person.getAttributes().getAttribute("city")); + case NONE -> null; + }; + + return new Key(gender, age, homeZone); + } + + private enum Zone { + FULL, CITY, NONE + } + + /** + * Key used to match persons. + */ + public record Key(String gender, int age, String homeZone) { + } + } diff --git a/src/main/java/org/matsim/prepare/population/CreateKansaiPopulation.java b/src/main/java/org/matsim/prepare/population/CreateKansaiPopulation.java index 7cc2a96..c355f91 100644 --- a/src/main/java/org/matsim/prepare/population/CreateKansaiPopulation.java +++ b/src/main/java/org/matsim/prepare/population/CreateKansaiPopulation.java @@ -20,13 +20,11 @@ import org.matsim.core.population.PopulationUtils; import org.matsim.core.scenario.ProjectionUtils; import org.matsim.core.utils.geometry.CoordinateTransformation; -import org.matsim.core.utils.geometry.transformations.GeotoolsTransformation; -import org.matsim.run.OpenBerlinScenario; +import org.matsim.run.OpenKyotoScenario; import picocli.CommandLine; import java.nio.file.Files; import java.nio.file.Path; -import java.text.NumberFormat; import java.text.ParseException; import java.util.*; @@ -36,10 +34,7 @@ ) public class CreateKansaiPopulation implements MATSimAppCommand { - private static final NumberFormat FMT = NumberFormat.getInstance(Locale.GERMAN); - private static final Logger log = LogManager.getLogger(CreateBerlinPopulation.class); - private final CoordinateTransformation ct = new GeotoolsTransformation("EPSG:32653", "EPSG:32653"); @CommandLine.Option(names = "--input", description = "Path to input csv data", required = true) private Path input; @@ -50,6 +45,9 @@ public class CreateKansaiPopulation implements MATSimAppCommand { @CommandLine.Mixin private ShpOptions shp = new ShpOptions(); + @CommandLine.Option(names = "--postal-shp", description = "Path to postal code shape file", required = true) + private String postalShp; + @CommandLine.Option(names = "--output", description = "Path to output population", required = true) private Path output; @@ -59,6 +57,9 @@ public class CreateKansaiPopulation implements MATSimAppCommand { private Map zones; private SplittableRandom rnd; private Population population; + private ShpOptions.Index postalIndex; + private CoordinateTransformation ct; + public static void main(String[] args) { new CreateKansaiPopulation().execute(args); @@ -77,7 +78,7 @@ private static double parseDouble(String s) { public Integer call() throws Exception { if (!shp.isDefined()) { - log.error("Shape file with LOR zones is required."); + log.error("Shape file with zones is required."); return 2; } @@ -86,6 +87,7 @@ public Integer call() throws Exception { rnd = new SplittableRandom(0); zones = new HashMap<>(); population = PopulationUtils.createPopulation(ConfigUtils.createConfig()); + ct = shp.createInverseTransformation(OpenKyotoScenario.CRS); // Collect all LORs for (SimpleFeature ft : fts) { @@ -94,6 +96,8 @@ public Integer call() throws Exception { log.info("Found {} zones", zones.size()); + postalIndex = new ShpOptions(postalShp, null, null).createIndex(OpenKyotoScenario.CRS, "_"); + CSVFormat.Builder format = CSVFormat.DEFAULT.builder().setDelimiter(',').setHeader().setSkipHeaderRecord(true); try (CSVParser reader = new CSVParser(Files.newBufferedReader(input), format.build())) { @@ -102,7 +106,7 @@ public Integer call() throws Exception { try { processZone(row); } catch (RuntimeException e) { - log.error("Error processing lor", e); + log.error("Error processing zone", e); log.error(row.toString()); } } @@ -112,7 +116,7 @@ public Integer call() throws Exception { PopulationUtils.sortPersons(population); - ProjectionUtils.putCRS(population, OpenBerlinScenario.CRS); + ProjectionUtils.putCRS(population, OpenKyotoScenario.CRS); PopulationUtils.writePopulation(population, output.toString()); return 0; @@ -173,11 +177,13 @@ private void processZone(CSVRecord row) throws ParseException { double men = parseDouble(row.get("male")); double womenQuota = women / (men + women); - double unemployed = parseDouble(row.get("unemployed")) / (parseDouble(row.get("total_em"))); - EnumeratedAttributeDistribution sex = new EnumeratedAttributeDistribution<>(Map.of("f", womenQuota, "m", 1 - womenQuota)); - EnumeratedAttributeDistribution employment = new EnumeratedAttributeDistribution<>(Map.of(true, 1 - unemployed, false, unemployed)); + // Employment is not used, because it can not be used from survey data + // double unemployed = parseDouble(row.get("unemployed")) / (parseDouble(row.get("total_em"))); +// EnumeratedAttributeDistribution employment = new EnumeratedAttributeDistribution<>(Map.of(true, 1 - unemployed, false, unemployed)); + EnumeratedAttributeDistribution ageGroup = buildAgeDist(row); + EnumeratedAttributeDistribution sex = new EnumeratedAttributeDistribution<>(Map.of("f", womenQuota, "m", 1 - womenQuota)); MultiPolygon geom = zones.get(zone); @@ -185,6 +191,8 @@ private void processZone(CSVRecord row) throws ParseException { double inhabitants = n * sample; + int tries = 0; + while (inhabitants > 0) { // If the number of inhabitants is less than 1, we use it as a probability @@ -194,7 +202,7 @@ private void processZone(CSVRecord row) throws ParseException { } else inhabitants--; - Person person = f.createPerson(CreateBerlinPopulation.generateId(population, "", rnd)); + Person person = f.createPerson(CreateBerlinPopulation.generateId(population, "p", rnd)); PersonUtils.setSex(person, sex.sample()); PopulationUtils.putSubpopulation(person, "person"); @@ -202,18 +210,28 @@ private void processZone(CSVRecord row) throws ParseException { int age = group.max == Integer.MAX_VALUE ? 100 : rnd.nextInt(group.min, group.max + 1); - PersonUtils.setAge(person, age); - if (age >= 15) - PersonUtils.setEmployed(person, employment.sample()); - Coord coord = ct.transform(CreateBerlinPopulation.sampleHomeCoordinate(geom, "EPSG:32653", facilities, rnd)); + Coord coord = ct.transform(CreateBerlinPopulation.sampleHomeCoordinate(geom, OpenKyotoScenario.CRS, facilities, rnd, 100)); person.getAttributes().putAttribute(Attributes.HOME_X, coord.getX()); person.getAttributes().putAttribute(Attributes.HOME_Y, coord.getY()); person.getAttributes().putAttribute("city", parseInt(row.get("city code"))); + SimpleFeature ft = postalIndex.queryFeature(coord); + + // Skip and generate another person if no postal code is found + if (ft == null && tries++ < 10) { + inhabitants++; + continue; + } else { + String postal = ft == null ? "NA" : Objects.toString(ft.getAttribute("zip_pre")) + ft.getAttribute("zip_mid"); + person.getAttributes().putAttribute("postal", postal); + } + tries = 0; + + Plan plan = f.createPlan(); plan.addActivity(f.createActivityFromCoord("home", coord));