From ed9068d9d8401752216b486c0505965b7fc5719b Mon Sep 17 00:00:00 2001 From: Cap Date: Thu, 13 Nov 2025 23:18:26 +0100 Subject: [PATCH 01/16] ongoing tests --- examples/quickstart-fr.py | 50 ++++++ mobility/choice_models/population_trips.py | 4 +- mobility/choice_models/results.py | 169 ++++++++++++++++++++- 3 files changed, 218 insertions(+), 5 deletions(-) create mode 100644 examples/quickstart-fr.py diff --git a/examples/quickstart-fr.py b/examples/quickstart-fr.py new file mode 100644 index 00000000..955e663b --- /dev/null +++ b/examples/quickstart-fr.py @@ -0,0 +1,50 @@ +import os +import dotenv +import logging + +import mobility + +dotenv.load_dotenv() + + +mobility.set_params( + # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"], + # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"] + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/test-09", +) + +# Using Foix (a small town) and a limited radius for quick results +transport_zones = mobility.TransportZones("fr-09122", radius = 10) + +# Using EMP, the latest national mobility survey for France +emp = mobility.EMPMobilitySurvey() + +# Creating a synthetic population of 1000 for the area +pop = mobility.Population(transport_zones, sample_size = 1000) + +# Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) +pop_trips = mobility.PopulationTrips( + pop, + [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], + [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], + [emp] + ) + +pop_trips_base = mobility.PopulationTrips( + pop, + [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], + [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], + [emp] + ) + +# You can get the weekday trips to inspect them +weekday_flows = pop_trips.get()["weekday_flows"].collect() + +# You can can also get global metrics that will be compared to the theoetical values for this population +tz = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base) +global_metrics = pop_trips.evaluate("global_metrics") + +# You can also plot the flows, with labels for the cities that are bigger than their neighbours +labels = pop_trips.get_prominent_cities() +pop_trips.plot_od_flows(labels=labels) diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index 761c7bb5..88c40544 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -100,7 +100,7 @@ def __init__( costs_aggregator = TravelCostsAggregator(modes) inputs = { - "version": 1, + "version": 2, "population": population, "costs_aggregator": costs_aggregator, "motives": motives, @@ -354,7 +354,7 @@ def run_model(self, is_weekday): costs = costs_aggregator.get_costs_by_od_and_mode( - ["distance", "time", "ghg_emissions"], + ["distance", "time", "ghg_emissions", "cost"], #To be tested congestion=True ) diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py index 41111b97..ee0f94d6 100644 --- a/mobility/choice_models/results.py +++ b/mobility/choice_models/results.py @@ -1,4 +1,5 @@ import json +import logging import polars as pl import numpy as np import plotly.express as px @@ -25,7 +26,7 @@ def __init__( weekend_chains, surveys, modes - ): + ): self.transport_zones = transport_zones self.demand_groups = demand_groups @@ -52,6 +53,7 @@ def __init__( "trip_count_by_demand_group": self.trip_count_by_demand_group, "distance_per_person": self.distance_per_person, "time_per_person": self.time_per_person, + "cost_per_person": self.cost_per_person, "immobility": self.immobility, "car_traffic": self.car_traffic, "travel_costs": self.travel_costs, @@ -744,8 +746,167 @@ def time_per_person( return time + + def cost_per_person( + self, + weekday: bool = True, + plot: bool = False, + mask_outliers: bool = False, + compare_with = None, + plot_delta = False + ): + """ + TODO + + Aggregate total travel time and time per person by demand group. + + Parameters + ---------- + weekday : bool, default True + Use weekday (True) or weekend (False) data. + plot : bool, default False + When True, shows a choropleth of average time per person by home zone. + + Returns + ------- + pl.DataFrame + Grouped by ['home_zone_id', 'csp', 'n_cars'] with: + - cost: sum(cost * n_persons) * 60.0 (minutes) + - n_persons: group size + - time_per_person: time / n_persons + """ + + states_steps = self.weekday_states_steps if weekday else self.weekend_states_steps + print("states_steps:") + print(type(states_steps)) + print(states_steps.columns) + print(states_steps) + costs = self.weekday_costs if weekday else self.weekend_costs + print("costs:") + print(costs.columns) + + transport_zones_df = ( + pl.DataFrame( + self.transport_zones.get().drop("geometry", axis=1) + ) + .filter(pl.col("is_inner_zone")) + .lazy() + ) + print("tzdf:") + print(transport_zones_df.columns) + + cost = ( + states_steps + .filter(pl.col("motive_seq_id") != 0) + .rename({"home_zone_id": "transport_zone_id"}) + .join( + transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), + on=["transport_zone_id"] + ) + .join(costs, on=["from", "to", "mode"]) + .group_by(["transport_zone_id", "csp", "n_cars"]) + .agg( + cost=(pl.col("cost")*pl.col("n_persons")).sum() + ) + .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) + .with_columns( + cost_per_person=pl.col("cost")/pl.col("n_persons") + ) + .collect(engine="streaming") + ) + + if compare_with is not None: + weekday_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekday_flows"]) + weekend_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekend_flows"]) + weekday_costs_comp = pl.scan_parquet(compare_with.cache_path["weekday_costs"]) + weekend_costs_comp = pl.scan_parquet(compare_with.cache_path["weekend_costs"]) + states_steps_comp = weekday_states_steps_comp if weekday else weekend_states_steps_comp + costs_comp = weekday_costs_comp if weekday else weekend_costs_comp + cost_comp = ( + states_steps_comp + .filter(pl.col("motive_seq_id") != 0) + .rename({"home_zone_id": "transport_zone_id"}) + .join( + transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), + on=["transport_zone_id"] + ) + .join(costs_comp, on=["from", "to", "mode"]) + .group_by(["transport_zone_id", "csp", "n_cars"]) + .agg( + cost_comp=(pl.col("cost")*pl.col("n_persons")).sum() + ) + .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) + .with_columns( + cost_per_person_comp=pl.col("cost_comp")/pl.col("n_persons") + ) + .collect(engine="streaming") + ) + cost = ( + cost + .join(cost_comp.select(["transport_zone_id", "csp", "n_cars", "cost_comp", "cost_per_person_comp"]), + on=["transport_zone_id", "csp", "n_cars"]) + .with_columns(delta=pl.col("cost_per_person")-pl.col("cost_per_person_comp")) + ) + + + if plot or plot_delta: + + tz = self.transport_zones.get() + tz = tz.to_crs(4326) + tz = tz.merge(transport_zones_df.collect().to_pandas(), on="transport_zone_id") + + tz = tz.merge( + ( + cost + .group_by(["transport_zone_id"]) + .agg( + cost_per_person=pl.col("cost").sum()/pl.col("n_persons").sum() + ) + .to_pandas() + ), + on="transport_zone_id", + how="left" + ) + + if plot_delta: + tz = tz.merge( + ( + cost + .group_by(["transport_zone_id"]) + .agg( + cost_per_person_comp=pl.col("cost_comp").sum()/pl.col("n_persons").sum() + # Assumption : always same number of persons in TZ + ) + .to_pandas() + ), + on="transport_zone_id", + how="left" + ) + tz["delta"] = tz["cost_per_person"] - tz["cost_per_person_comp"] + + + if plot: + tz["cost_per_person"] = tz["cost_per_person"].fillna(0.0) + + if mask_outliers: + tz["cost_per_person"] = self.mask_outliers(tz["cost_per_person"]) + + self.plot_map(tz, "cost_per_person") + + + if plot_delta: + tz["delta"] = tz["delta"].fillna(0.0) + + if mask_outliers: + tz["delta"] = self.mask_outliers(tz["delta"]) + + self.plot_map(tz, "delta") + + + return cost - def plot_map(self, tz, value: str = None, motive: str = None): + + def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = "browser"): """ Render a Plotly choropleth for a transport-zone metric. @@ -762,6 +923,8 @@ def plot_map(self, tz, value: str = None, motive: str = None): None Displays an interactive map in the browser. """ + logging.getLogger("kaleido").setLevel(logging.WARNING) + #plot_method="png" fig = px.choropleth( tz.drop(columns="geometry"), @@ -777,7 +940,7 @@ def plot_map(self, tz, value: str = None, motive: str = None): ) fig.update_geos(fitbounds="geojson", visible=False) fig.update_layout(margin=dict(l=0,r=0,t=0,b=0)) - fig.show("browser") + fig.show(plot_method) def mask_outliers(self, series): From 2467fa0845751ec6c866c86c4861715acacb6105 Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 17 Nov 2025 10:06:47 +0100 Subject: [PATCH 02/16] ongoing tests --- examples/quickstart-fr.py | 16 +- examples/quickstart-fr.py.py | 40 ---- mobility/choice_models/results.py | 358 ++++++++++++------------------ pyproject.toml | 3 +- 4 files changed, 153 insertions(+), 264 deletions(-) delete mode 100644 examples/quickstart-fr.py.py diff --git a/examples/quickstart-fr.py b/examples/quickstart-fr.py index 955e663b..ae9f0601 100644 --- a/examples/quickstart-fr.py +++ b/examples/quickstart-fr.py @@ -4,6 +4,8 @@ import mobility +import polars as pl + dotenv.load_dotenv() @@ -28,23 +30,31 @@ pop, [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], - [emp] + [emp], + parameters=mobility.PopulationTripsParameters(n_iterations=10) ) pop_trips_base = mobility.PopulationTrips( pop, [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], - [emp] + [emp], + parameters=mobility.PopulationTripsParameters(n_iterations=10, seed=49283092) ) +tz = transport_zones.get() + # You can get the weekday trips to inspect them weekday_flows = pop_trips.get()["weekday_flows"].collect() # You can can also get global metrics that will be compared to the theoetical values for this population -tz = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base) +cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base) +dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base) +time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base) +ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base) global_metrics = pop_trips.evaluate("global_metrics") + # You can also plot the flows, with labels for the cities that are bigger than their neighbours labels = pop_trips.get_prominent_cities() pop_trips.plot_od_flows(labels=labels) diff --git a/examples/quickstart-fr.py.py b/examples/quickstart-fr.py.py deleted file mode 100644 index ebc32e82..00000000 --- a/examples/quickstart-fr.py.py +++ /dev/null @@ -1,40 +0,0 @@ -import os -import dotenv - -import mobility - -dotenv.load_dotenv() - -mobility.set_params( - # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"], - # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"] - package_data_folder_path="D:/mobility-data", - project_data_folder_path="D:/test-09", -) - -# Using Foix (a small town) and a limited radius for quick results -transport_zones = mobility.TransportZones("fr-09122", radius = 10) - -# Using EMP, the latest national mobility survey for France -emp = mobility.EMPMobilitySurvey() - -# Creating a synthetic population of 1000 for the area -pop = mobility.Population(transport_zones, sample_size = 1000) - -# Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) -pop_trips = mobility.PopulationTrips( - pop, - [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], - [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], - [emp] - ) - -# You can get the weekday trips to inspect them -weekday_flows = pop_trips.get()["weekday_flows"].collect() - -# You can can also get global metrics that will be compared to the theoetical values for this population -global_metrics = pop_trips.evaluate("global_metrics") - -# You can also plot the flows, with labels for the cities that are bigger than their neighbours -labels = pop_trips.get_prominent_cities() -pop_trips.plot_od_flows(labels=labels) diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py index ee0f94d6..b8d4919b 100644 --- a/mobility/choice_models/results.py +++ b/mobility/choice_models/results.py @@ -52,6 +52,7 @@ def __init__( "sink_occupation": self.sink_occupation, "trip_count_by_demand_group": self.trip_count_by_demand_group, "distance_per_person": self.distance_per_person, + "ghg_per_person": self.ghg_per_person, "time_per_person": self.time_per_person, "cost_per_person": self.cost_per_person, "immobility": self.immobility, @@ -574,181 +575,10 @@ def trip_count_by_demand_group( self.plot_map(tz, "n_trips_per_person") return trip_count - - - def distance_per_person( - self, - weekday: bool = True, - plot: bool = False, - mask_outliers: bool = False - ): - """ - Aggregate total travel distance and distance per person by demand group. - - Parameters - ---------- - weekday : bool, default True - Use weekday (True) or weekend (False) data. - plot : bool, default False - When True, shows a choropleth of average distance per person by home zone. - - Returns - ------- - pl.DataFrame - Grouped by ['home_zone_id', 'csp', 'n_cars'] with: - - distance: sum(distance * n_persons) - - n_persons: group size - - distance_per_person: distance / n_persons - """ - - states_steps = self.weekday_states_steps if weekday else self.weekend_states_steps - costs = self.weekday_costs if weekday else self.weekend_costs - - transport_zones_df = ( - pl.DataFrame( - self.transport_zones.get().drop("geometry", axis=1) - ) - .filter(pl.col("is_inner_zone")) - .lazy() - ) - - distance = ( - states_steps - .filter(pl.col("motive_seq_id") != 0) - .rename({"home_zone_id": "transport_zone_id"}) - .join( - transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), - on=["transport_zone_id"] - ) - .join(costs, on=["from", "to", "mode"]) - .group_by(["transport_zone_id", "csp", "n_cars"]) - .agg( - distance=(pl.col("distance")*pl.col("n_persons")).sum() - ) - .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) - .with_columns( - distance_per_person=pl.col("distance")/pl.col("n_persons") - ) - .collect(engine="streaming") - ) - - if plot: - - tz = self.transport_zones.get() - tz = tz.to_crs(4326) - tz = tz.merge(transport_zones_df.collect().to_pandas(), on="transport_zone_id") - - tz = tz.merge( - ( - distance - .group_by(["transport_zone_id"]) - .agg( - distance_per_person=pl.col("distance").sum()/pl.col("n_persons").sum() - ) - .to_pandas() - ), - on="transport_zone_id", - how="left" - ) - - tz["distance_per_person"] = tz["distance_per_person"].fillna(0.0) - - if mask_outliers: - tz["distance_per_person"] = self.mask_outliers(tz["distance_per_person"]) - self.plot_map(tz, "distance_per_person") - - return distance - - - def time_per_person( - self, - weekday: bool = True, - plot: bool = False, - mask_outliers: bool = False - ): - """ - Aggregate total travel time and time per person by demand group. - - Parameters - ---------- - weekday : bool, default True - Use weekday (True) or weekend (False) data. - plot : bool, default False - When True, shows a choropleth of average time per person by home zone. - - Returns - ------- - pl.DataFrame - Grouped by ['home_zone_id', 'csp', 'n_cars'] with: - - time: sum(time * n_persons) * 60.0 (minutes) - - n_persons: group size - - time_per_person: time / n_persons - """ - - states_steps = self.weekday_states_steps if weekday else self.weekend_states_steps - costs = self.weekday_costs if weekday else self.weekend_costs - - transport_zones_df = ( - pl.DataFrame( - self.transport_zones.get().drop("geometry", axis=1) - ) - .filter(pl.col("is_inner_zone")) - .lazy() - ) - - time = ( - states_steps - .filter(pl.col("motive_seq_id") != 0) - .rename({"home_zone_id": "transport_zone_id"}) - .join( - transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), - on=["transport_zone_id"] - ) - .join(costs, on=["from", "to", "mode"]) - .group_by(["transport_zone_id", "csp", "n_cars"]) - .agg( - time=(pl.col("time")*pl.col("n_persons")).sum()*60.0 - ) - .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) - .with_columns( - time_per_person=pl.col("time")/pl.col("n_persons") - ) - .collect(engine="streaming") - ) - - if plot: - - tz = self.transport_zones.get() - tz = tz.to_crs(4326) - tz = tz.merge(transport_zones_df.collect().to_pandas(), on="transport_zone_id") - - tz = tz.merge( - ( - time - .group_by(["transport_zone_id"]) - .agg( - time_per_person=pl.col("time").sum()/pl.col("n_persons").sum() - ) - .to_pandas() - ), - on="transport_zone_id", - how="left" - ) - - tz["time_per_person"] = tz["time_per_person"].fillna(0.0) - - if mask_outliers: - tz["time_per_person"] = self.mask_outliers(tz["time_per_person"]) - - self.plot_map(tz, "time_per_person") - - - return time - - - def cost_per_person( + def metric_per_person( self, + metric: str, weekday: bool = True, plot: bool = False, mask_outliers: bool = False, @@ -777,13 +607,20 @@ def cost_per_person( """ states_steps = self.weekday_states_steps if weekday else self.weekend_states_steps - print("states_steps:") - print(type(states_steps)) - print(states_steps.columns) - print(states_steps) costs = self.weekday_costs if weekday else self.weekend_costs - print("costs:") - print(costs.columns) + + if compare_with is not None: + try: + compare_with.get() + weekday_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekday_flows"]) + weekend_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekend_flows"]) + weekday_costs_comp = pl.scan_parquet(compare_with.cache_path["weekday_costs"]) + weekend_costs_comp = pl.scan_parquet(compare_with.cache_path["weekend_costs"]) + states_steps_comp = weekday_states_steps_comp if weekday else weekend_states_steps_comp + costs_comp = weekday_costs_comp if weekday else weekend_costs_comp + except: + Exception("The PopulationsTrips to compare with did not work. Try to run them alone?") + transport_zones_df = ( pl.DataFrame( @@ -792,60 +629,58 @@ def cost_per_person( .filter(pl.col("is_inner_zone")) .lazy() ) - print("tzdf:") - print(transport_zones_df.columns) + + metric_per_person = metric + "_per_person" - cost = ( + metric_per_groups_and_transport_zones = ( states_steps .filter(pl.col("motive_seq_id") != 0) .rename({"home_zone_id": "transport_zone_id"}) - .join( - transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), - on=["transport_zone_id"] - ) .join(costs, on=["from", "to", "mode"]) .group_by(["transport_zone_id", "csp", "n_cars"]) .agg( - cost=(pl.col("cost")*pl.col("n_persons")).sum() + metric=(pl.col(metric)*pl.col("n_persons")).sum() ) - .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) + .join( + transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), + on=["transport_zone_id"] + ) .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) .with_columns( - cost_per_person=pl.col("cost")/pl.col("n_persons") + metric_per_person=pl.col("metric")/pl.col("n_persons") ) + .rename({"metric": metric, "metric_per_person": metric_per_person}) .collect(engine="streaming") ) - + if compare_with is not None: - weekday_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekday_flows"]) - weekend_states_steps_comp = pl.scan_parquet(compare_with.cache_path["weekend_flows"]) - weekday_costs_comp = pl.scan_parquet(compare_with.cache_path["weekday_costs"]) - weekend_costs_comp = pl.scan_parquet(compare_with.cache_path["weekend_costs"]) - states_steps_comp = weekday_states_steps_comp if weekday else weekend_states_steps_comp - costs_comp = weekday_costs_comp if weekday else weekend_costs_comp - cost_comp = ( + metric_comp = metric + "_comp" + metric_per_person_comp = metric + "_per_person_comp" + metric_per_groups_and_transport_zones_comp = ( states_steps_comp .filter(pl.col("motive_seq_id") != 0) .rename({"home_zone_id": "transport_zone_id"}) - .join( - transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), - on=["transport_zone_id"] - ) .join(costs_comp, on=["from", "to", "mode"]) .group_by(["transport_zone_id", "csp", "n_cars"]) .agg( - cost_comp=(pl.col("cost")*pl.col("n_persons")).sum() + metric=(pl.col(metric)*pl.col("n_persons")).sum() + ) + .join( + transport_zones_df.select(["transport_zone_id", "local_admin_unit_id"]), + on=["transport_zone_id"] ) .join(self.demand_groups.rename({"home_zone_id": "transport_zone_id"}), on=["transport_zone_id", "csp", "n_cars"]) .with_columns( - cost_per_person_comp=pl.col("cost_comp")/pl.col("n_persons") + metric_per_person=pl.col("metric")/pl.col("n_persons") ) + .rename({"metric": metric_comp, "metric_per_person": metric_per_person_comp}) .collect(engine="streaming") ) - cost = ( - cost - .join(cost_comp.select(["transport_zone_id", "csp", "n_cars", "cost_comp", "cost_per_person_comp"]), - on=["transport_zone_id", "csp", "n_cars"]) - .with_columns(delta=pl.col("cost_per_person")-pl.col("cost_per_person_comp")) + metric_per_groups_and_transport_zones = ( + metric_per_groups_and_transport_zones + .join(metric_per_groups_and_transport_zones_comp.select(["transport_zone_id", "csp", "n_cars", "n_persons", "local_admin_unit_id", + metric_comp, metric_per_person_comp]), + on=["transport_zone_id", "csp", "n_cars"]) # Same TZ ids? + .with_columns(delta=pl.col(metric_per_person)-pl.col(metric_per_person_comp)) ) @@ -857,11 +692,12 @@ def cost_per_person( tz = tz.merge( ( - cost + metric_per_groups_and_transport_zones .group_by(["transport_zone_id"]) .agg( - cost_per_person=pl.col("cost").sum()/pl.col("n_persons").sum() + metric_per_person=pl.col(metric).sum()/pl.col("n_persons").sum() ) + .rename({"metric_per_person": metric_per_person}) .to_pandas() ), on="transport_zone_id", @@ -871,27 +707,28 @@ def cost_per_person( if plot_delta: tz = tz.merge( ( - cost + metric_per_groups_and_transport_zones .group_by(["transport_zone_id"]) .agg( - cost_per_person_comp=pl.col("cost_comp").sum()/pl.col("n_persons").sum() + metric_per_person_comp=pl.col(metric_comp).sum()/pl.col("n_persons").sum() # Assumption : always same number of persons in TZ ) + .rename({"metric_per_person_comp": metric_per_person_comp}) .to_pandas() ), on="transport_zone_id", how="left" ) - tz["delta"] = tz["cost_per_person"] - tz["cost_per_person_comp"] + tz["delta"] = tz[metric_per_person] - tz[metric_per_person_comp] if plot: - tz["cost_per_person"] = tz["cost_per_person"].fillna(0.0) + tz[metric_per_person] = tz[metric_per_person].fillna(0.0) if mask_outliers: - tz["cost_per_person"] = self.mask_outliers(tz["cost_per_person"]) + tz[metric_per_person] = self.mask_outliers(tz[metric_per_person]) - self.plot_map(tz, "cost_per_person") + self.plot_map(tz, metric_per_person) if plot_delta: @@ -900,13 +737,93 @@ def cost_per_person( if mask_outliers: tz["delta"] = self.mask_outliers(tz["delta"]) - self.plot_map(tz, "delta") + self.plot_map(tz, "delta", color_continuous_scale="RdBu_r", color_continuous_midpoint=0) - return cost + return metric_per_groups_and_transport_zones + + + def distance_per_person(self, *args, **kwargs): + """ + Aggregate total travel distance and distance per person by demand group. + + Parameters + ---------- + weekday : bool, default True + Use weekday (True) or weekend (False) data. + plot : bool, default False + When True, shows a choropleth of average distance per person by home zone. + + Returns + ------- + pl.DataFrame + Grouped by ['home_zone_id', 'csp', 'n_cars'] with: + - distance: sum(distance * n_persons) + - n_persons: group size + - distance_per_person: distance / n_persons + """ + return self.metric_per_person("distance", *args, **kwargs) + + def ghg_per_person(self, *args, **kwargs): + return self.metric_per_person("ghg_emissions_per_trip", *args, **kwargs) + + def time_per_person(self, *args, **kwargs): + """ + Aggregate total travel time and time per person by demand group. + + Parameters + ---------- + weekday : bool, default True + Use weekday (True) or weekend (False) data. + plot : bool, default False + When True, shows a choropleth of average time per person by home zone. + + Returns + ------- + pl.DataFrame + Grouped by ['home_zone_id', 'csp', 'n_cars'] with: + - time: sum(time * n_persons) * 60.0 (minutes) + - n_persons: group size + - time_per_person: time / n_persons + """ + return self.metric_per_person("time", *args, **kwargs) + + + + def cost_per_person(self, *args, **kwargs): + """ self, + weekday: bool = True, + plot: bool = False, + mask_outliers: bool = False, + compare_with = None, + plot_delta = False + ): + + TODO + + Aggregate total travel time and time per person by demand group. + + Parameters + ---------- + weekday : bool, default True + Use weekday (True) or weekend (False) data. + plot : bool, default False + When True, shows a choropleth of average time per person by home zone. + + Returns + ------- + pl.DataFrame + Grouped by ['home_zone_id', 'csp', 'n_cars'] with: + - cost: sum(cost * n_persons) * 60.0 (minutes) + - n_persons: group size + - time_per_person: time / n_persons + """ + return self.metric_per_person("cost", *args, **kwargs) + - def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = "browser"): + def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = "browser", + color_continuous_scale="Viridis", color_continuous_midpoint=None): """ Render a Plotly choropleth for a transport-zone metric. @@ -933,7 +850,8 @@ def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = featureidkey="properties.transport_zone_id", color=value, hover_data=["transport_zone_id", value], - color_continuous_scale="Viridis", + color_continuous_scale=color_continuous_scale, + color_continuous_midpoint= color_continuous_midpoint, projection="mercator", title=motive, subtitle=motive diff --git a/pyproject.toml b/pyproject.toml index 66e06b17..dd0399eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,8 @@ dependencies = [ "networkx", "plotly", "scikit-learn", - "gtfs_kit" + "gtfs_kit", + "kaleido" ] requires-python = ">=3.11" From 76e88e59c744129519f98a09bea61734bbb3e47a Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 12:12:22 +0100 Subject: [PATCH 03/16] Parameter class New classes Parameter and ParameterSet to better manager parameters of the model --- mobility/__init__.py | 2 + .../population_trips_parameters.py | 110 +++++++++++++++-- mobility/model_parameters.py | 111 ++++++++++++++++++ 3 files changed, 210 insertions(+), 13 deletions(-) create mode 100644 mobility/model_parameters.py diff --git a/mobility/__init__.py b/mobility/__init__.py index 8a446432..65efae4c 100644 --- a/mobility/__init__.py +++ b/mobility/__init__.py @@ -1,5 +1,7 @@ from .set_params import set_params +from .model_parameters import Parameter, ParameterSet + from .study_area import StudyArea from .transport_zones import TransportZones diff --git a/mobility/choice_models/population_trips_parameters.py b/mobility/choice_models/population_trips_parameters.py index 85e8425d..0f5bdfe1 100644 --- a/mobility/choice_models/population_trips_parameters.py +++ b/mobility/choice_models/population_trips_parameters.py @@ -1,24 +1,108 @@ -from dataclasses import dataclass +from dataclasses import dataclass, fields +from typing import List, Union -@dataclass(frozen=True) -class PopulationTripsParameters: - n_iterations: int = 1 +from mobility import Parameter, ParameterSet + + +Number = Union[int, float] + + +@dataclass +class PopulationTripsParameters(ParameterSet): + """To add a new parameter: + - add name, type and default value here, + - add the description in the __post_init__ method + - if the parameter is linked to another one, maybe add a constraint + in _validate_param_interdependency() + """ + n_iterations: int = 4 alpha: float = 0.01 k_mode_sequences: int = 6 dest_prob_cutoff: float = 0.99 n_iter_per_cost_update: int = 3 cost_uncertainty_sd: float = 1.0 seed: int = 0 + seeds: List[int] = None mode_sequence_search_parallel: bool = True min_activity_time_constant: float = 1.0 simulate_weekend: bool = False + name: str = "" + + + def __post_init__(self): + self.parameters = { + "param_n_iterations": Parameter(name="Number of iterations", name_fr="Nombre d'itérations", + value=self.n_iterations, + description=("""Number of simulation iterations used to compute the population trips. + Increase this to get more diverse programmes and to allow congestion feedbacks to propagate"""), + parameter_type=int, min_value=1, max_value=42, interval=1, + source_default="experience", parameter_role="simulation_parameter"), + + "param_alpha" : Parameter(name="Alpha", name_fr="Alpha", description="Alpha", + value=self.alpha, + parameter_type=float, min_value=0.0, max_value=1.0, interval=0.005, + source_default="experience", parameter_role="simulation_parameter"), + "param_k_mode_sequences" : Parameter(name="Number k of sampled mode sequences", + value=self.k_mode_sequences, + name_fr="Nombre k de séquences échantillonnées", + description="Increase this to get more mode sequence options for each destination sequence, decrease to speed up the simulation", + parameter_type=int, min_value=1, interval=1, + source_default="experience", parameter_role="simulation_parameter"), + "param_dest_prob_cutoff" : Parameter(name="TBD", name_fr="TBD", + value=self.dest_prob_cutoff, + description="Increase this to get more destinations possibilities, decrease to speed up the simulation", + parameter_type=float, min_value=0.10, max_value=1.00, interval=0.05, + source_default="experience", parameter_role="simulation_parameter"), + "param_n_iter_per_cost_update" : Parameter(name="Number n of iteration between two cost updates", + value=self.n_iter_per_cost_update, + name_fr="Nombre n d'itérations entre deux mises à jours des coûts", + description="""Every n iterations, the model will update the costs: + set to zero to ignore congestion in the simulation, + set to 1 to update congestion at each iteration, + set to a higher level to speed up the simulation""", + parameter_type=int, min_value=0, max_value=99, interval=1, + source_default="experience", parameter_role="simulation_parameter"), + "param_cost_uncertainty_sd" : Parameter(name="TBD", name_fr="TBD", + value=self.cost_uncertainty_sd, + description="Increase this to make the location of opportunities more uncertain around the destination points", + parameter_type=float, min_value=0.1, interval=0.1, + source_default="experience", parameter_role="simulation_parameter"), + "param_seed" : Parameter(name="Seed", name_fr="Graine (seed)", + value=self.seed, + description="""Unique seed used for these trips. + Change this value to get new programmes and results for a given set of inputs, + keep this to a set value to get reproductible results.""", + parameter_type=int, min_value=0, + parameter_role="randomisation"), + "param_seeds" : Parameter(name="Set of seeds", name_fr="Ensemble de graines (seeds)", + value=self.seeds, + description="Set of seeds, used to improve resultats by minimising noise", + parameter_type=list, + parameter_role="randomisation"), + "param_mode_sequence_search_parallel" : Parameter(name="TBD", name_fr="TBD", + value=self.mode_sequence_search_parallel, + description="Set to False to debug or for small simulations, otherwise set to True to speed up the simulation", + parameter_type=bool, + source_default="tests", parameter_role="simulation_time_parameter"), + "param_min_activity_time_constant" : Parameter(name="TBD", name_fr="TBD", description="TBD", + value=self.min_activity_time_constant, + parameter_type=float, interval=0.1, + source_default="experience", parameter_role="simulation_parameter"), + "param_simulate_weekend" : Parameter(name="Simulate weekend", name_fr="Simulation du week-end", + value=self.simulate_weekend, + description="Simulate or not trips for the week-ends", + parameter_type=bool, + source_default="speed_optimisation", parameter_role="simulation_parameter"), + "param_name": Parameter(name="Name", name_fr="Nom", description="Name of the run with those parameters", + value=self.name, parameter_type=str)} + + def _validate_param_interdependency(self): + if self.seeds is not None: + for seed in self.seeds: + if (not isinstance(seed, int) or seed < 0): + raise TypeError(f"Seed must be a positive int and not {seed}") + if self.n_iter_per_cost_update > self.n_iterations: + raise ValueError(f"n_iter_per_cost_update ({self.n_iter_per_cost_update}) should not be higher than the number of iterations ({self.n_iterations})") + - def validate(self) -> None: - assert self.n_iterations >= 1 - assert 0.0 < self.dest_prob_cutoff <= 1.0 - assert self.alpha >= 0.0 - assert self.k_mode_sequences >= 1 - assert self.n_iter_per_cost_update >= 0 - assert self.cost_uncertainty_sd > 0.0 - assert self.seed >= 0 - assert self.min_activity_time_constant >= 0 + diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py new file mode 100644 index 00000000..d2ab7099 --- /dev/null +++ b/mobility/model_parameters.py @@ -0,0 +1,111 @@ +from dataclasses import dataclass, field, fields +from typing import List, Union + +Number = Union[int, float] + +@dataclass(frozen=True) +class Parameter: + name: str + name_fr: str + value: Number | bool + description: str + parameter_type: type + #default_value: Number | bool + possible_values: List[float] | List[int] | tuple = None + min_value: Number = None + max_value: Number = None + unit: str = None + interval: Number = None + source_default: str = "" + parameter_role: str = "" + + + + # def get(self): + # """Return parameter value.""" + # val = self.default_value + # self._validate(val) + # return val + + def set(self, new_value): + self.value = new_value + + def _validate(self): + if self.value is None: #todo: improve this! + return None + if self.parameter_type is not None: + if not isinstance(self.value, self.parameter_type): + t = type(self.value) + raise TypeError(f"Parameter '{self.name}' must be {self.parameter_type} (currently {t}).") + + if self.min_value is not None and self.value < self.min_value: + raise ValueError( + f"Parameter '{self.name}' below minimum {self.min_value}" + ) + + if self.max_value is not None and self.value > self.max_value: + raise ValueError( + f"Parameter '{self.name}' above maximum {self.max_value}" + ) + + def __repr__(self): + print("pop") + unit_str = f" [{self.unit}]" if self.unit else "" + return f"" + + + def get_values_for_sensitivity_analyses(self, i_max=10): + value = self.value + print(self) + print(value) + values = [value] + if self.interval is None: + raise ValueError("To run a sensitivity analysis, interval must be specified in the parameter configuration") + i = 0 + if self.max_value is not None: + print("hop") + while i < i_max and value < self.max_value: + value += self.interval + values.append(round(value,3)) + i += 1 + else: + while i < i_max: + value += self.interval + values.append(round(value,3)) + i += 1 + value = self.value + i = 0 + + if self.min_value is not None: + while i < i_max and value > self.min_value: + value -= self.interval + print(i, value) + values.append(round(value,3)) + i += 1 + else: + while i < i_max: + value -= self.interval + values.append(round(value,3)) + i += 1 + + values.sort() + print(values) + return values + +@dataclass +class ParameterSet: + parameters : dict = field(init=False) + + def _validate(self): + # for param in fields(self): + # print(param) + # getattr(self, param.name)._validate() + for param in fields(self)[1:]: + print(self.parameters) + param_name = "param_" + param.name + print(param_name) + self.parameters[param_name]._validate() + self._validate_param_interdependency() + + def _validate_param_interdependency(self): + pass \ No newline at end of file From a3d515f6a77965379d6bf002bfe2fe4150f27026 Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 14:31:10 +0100 Subject: [PATCH 04/16] correct method name --- mobility/model_parameters.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py index d2ab7099..9b9f6b6d 100644 --- a/mobility/model_parameters.py +++ b/mobility/model_parameters.py @@ -30,7 +30,7 @@ class Parameter: def set(self, new_value): self.value = new_value - def _validate(self): + def validate(self): if self.value is None: #todo: improve this! return None if self.parameter_type is not None: @@ -96,15 +96,12 @@ def get_values_for_sensitivity_analyses(self, i_max=10): class ParameterSet: parameters : dict = field(init=False) - def _validate(self): - # for param in fields(self): - # print(param) - # getattr(self, param.name)._validate() + def validate(self): for param in fields(self)[1:]: print(self.parameters) param_name = "param_" + param.name print(param_name) - self.parameters[param_name]._validate() + self.parameters[param_name].validate() self._validate_param_interdependency() def _validate_param_interdependency(self): From 3c2eea812ee24995c0a89f410b696db910472d0c Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 15:31:20 +0100 Subject: [PATCH 05/16] other improvements and experiments --- examples/quickstart-fr.py | 48 +++++++---- mobility/choice_models/population_trips.py | 32 ++++---- mobility/choice_models/results.py | 33 ++++++-- mobility/experiments/gtfs_experiments.py | 15 ++++ mobility/experiments/multiple_seeds.py | 96 ++++++++++++++++++++++ 5 files changed, 190 insertions(+), 34 deletions(-) create mode 100644 mobility/experiments/gtfs_experiments.py create mode 100644 mobility/experiments/multiple_seeds.py diff --git a/examples/quickstart-fr.py b/examples/quickstart-fr.py index ae9f0601..5d0d84ba 100644 --- a/examples/quickstart-fr.py +++ b/examples/quickstart-fr.py @@ -14,6 +14,7 @@ # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"] package_data_folder_path="D:/mobility-data", project_data_folder_path="D:/test-09", + debug=False ) # Using Foix (a small town) and a limited radius for quick results @@ -25,36 +26,53 @@ # Creating a synthetic population of 1000 for the area pop = mobility.Population(transport_zones, sample_size = 1000) +modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)] +surveys = [emp] +motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)] + + # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) pop_trips = mobility.PopulationTrips( pop, - [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], - [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], - [emp], - parameters=mobility.PopulationTripsParameters(n_iterations=10) + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=42) ) pop_trips_base = mobility.PopulationTrips( pop, - [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)], - [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)], - [emp], - parameters=mobility.PopulationTripsParameters(n_iterations=10, seed=49283092) + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=4, seed=49283092, k_mode_sequences=42) ) tz = transport_zones.get() + # You can get the weekday trips to inspect them weekday_flows = pop_trips.get()["weekday_flows"].collect() +weekday_flows_base = pop_trips_base.get()["weekday_flows"].collect() + +# You can also plot the flows, with labels for the cities that are bigger than their neighbours +labels=pop_trips.get_prominent_cities() # You can can also get global metrics that will be compared to the theoetical values for this population -cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base) -dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base) -time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base) -ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base) +cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=True, compare_with=pop_trips_base, labels=labels, plot=True) +dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True) +time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True) +ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base, plot=True) global_metrics = pop_trips.evaluate("global_metrics") -# You can also plot the flows, with labels for the cities that are bigger than their neighbours -labels = pop_trips.get_prominent_cities() -pop_trips.plot_od_flows(labels=labels) +pop_trips.plot_modal_share(mode="public_transport", labels=labels) +pop_trips.plot_modal_share(mode="bicycle") +pop_trips.plot_modal_share(mode="walk") +cms = pop_trips.plot_modal_share(mode="car") + +# OD flows between transport zones +pop_trips.plot_od_flows(mode="car", level_of_detail=1, labels=labels) +pop_trips.plot_od_flows(mode="walk", level_of_detail=1) +pop_trips.plot_od_flows(mode="bicycle", level_of_detail=1) +pop_trips.plot_od_flows(mode="public_transport") \ No newline at end of file diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index 88c40544..aafcd674 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -128,7 +128,9 @@ def __init__( } super().__init__(inputs, cache_path) - + + def has_multiple_seeds(self): + return (len(self.parameters.seeds) > 1) def resolve_parameters( self, @@ -459,7 +461,7 @@ def evaluate(self, metric, **kwargs): return evaluation - def plot_modal_share(self, zone="origin", mode="car", period="weekdays", + def plot_modal_share(self, zone="origin", mode="car", period="weekdays", plot=True, labels=None, labels_size=[10, 6, 4], labels_color="black"): """ Plot modal share for the given mode in the origin or destination zones during weekdays or weekends. @@ -480,7 +482,6 @@ def plot_modal_share(self, zone="origin", mode="car", period="weekdays", Mode share for the given mode in each transport zone. """ - logging.info(f"🗺️ Plotting {mode} modal share for {zone} zones during {period}") if period == "weekdays": population_df = self.get()["weekday_flows"].collect().to_pandas() @@ -508,18 +509,21 @@ def plot_modal_share(self, zone="origin", mode="car", period="weekdays", mode_name = mode.capitalize() mode_share = mode_share[mode_share["mode"] == mode] - transport_zones_df = self.population.transport_zones.get() - - gc = gpd.GeoDataFrame( - transport_zones_df.merge(mode_share, how="left", right_on=left_column, left_on="transport_zone_id", suffixes=('', '_z'))).fillna(0) - gcp = gc.plot("modal_share", legend=True) - gcp.set_axis_off() - plt.title(f"{mode_name} share per {zone} transport zone ({period})") + if plot: + logging.info(f"🗺️ Plotting {mode} modal share for {zone} zones during {period}") - if isinstance(labels, gpd.GeoDataFrame): - self.__show_labels__(labels, labels_size, labels_color) - - plt.show() + transport_zones_df = self.population.transport_zones.get() + + gc = gpd.GeoDataFrame( + transport_zones_df.merge(mode_share, how="left", right_on=left_column, left_on="transport_zone_id", suffixes=('', '_z'))).fillna(0) + gcp = gc.plot("modal_share", legend=True) + gcp.set_axis_off() + plt.title(f"{mode_name} share per {zone} transport zone ({period})") + + if isinstance(labels, gpd.GeoDataFrame): + self.__show_labels__(labels, labels_size, labels_color) + + plt.show() return mode_share diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py index b8d4919b..e7b676f9 100644 --- a/mobility/choice_models/results.py +++ b/mobility/choice_models/results.py @@ -3,6 +3,7 @@ import polars as pl import numpy as np import plotly.express as px +import plotly.graph_objects as go from typing import Literal from mobility.choice_models.evaluation.travel_costs_evaluation import TravelCostsEvaluation @@ -53,6 +54,7 @@ def __init__( "trip_count_by_demand_group": self.trip_count_by_demand_group, "distance_per_person": self.distance_per_person, "ghg_per_person": self.ghg_per_person, + "ghg_emissions_per_trip_per_person": self.ghg_emissions_per_trip, #Test "time_per_person": self.time_per_person, "cost_per_person": self.cost_per_person, "immobility": self.immobility, @@ -583,7 +585,8 @@ def metric_per_person( plot: bool = False, mask_outliers: bool = False, compare_with = None, - plot_delta = False + plot_delta = False, + labels = None ): """ TODO @@ -728,7 +731,7 @@ def metric_per_person( if mask_outliers: tz[metric_per_person] = self.mask_outliers(tz[metric_per_person]) - self.plot_map(tz, metric_per_person) + self.plot_map(tz, metric_per_person, labels=labels) if plot_delta: @@ -737,7 +740,7 @@ def metric_per_person( if mask_outliers: tz["delta"] = self.mask_outliers(tz["delta"]) - self.plot_map(tz, "delta", color_continuous_scale="RdBu_r", color_continuous_midpoint=0) + self.plot_map(tz, "delta", color_continuous_scale="RdBu_r", color_continuous_midpoint=0, labels=labels) return metric_per_groups_and_transport_zones @@ -767,6 +770,10 @@ def distance_per_person(self, *args, **kwargs): def ghg_per_person(self, *args, **kwargs): return self.metric_per_person("ghg_emissions_per_trip", *args, **kwargs) + def ghg_emissions_per_trip(self, *args, **kwargs): + """Test for compat""" + return self.metric_per_person("ghg_emissions_per_trip", *args, **kwargs) + def time_per_person(self, *args, **kwargs): """ Aggregate total travel time and time per person by demand group. @@ -823,7 +830,8 @@ def cost_per_person(self, *args, **kwargs): def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = "browser", - color_continuous_scale="Viridis", color_continuous_midpoint=None): + color_continuous_scale="Viridis", color_continuous_midpoint=None, + labels=None): """ Render a Plotly choropleth for a transport-zone metric. @@ -843,6 +851,8 @@ def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = logging.getLogger("kaleido").setLevel(logging.WARNING) #plot_method="png" + + fig = px.choropleth( tz.drop(columns="geometry"), geojson=json.loads(tz.to_json()), @@ -856,6 +866,19 @@ def plot_map(self, tz, value: str = None, motive: str = None, plot_method: str = title=motive, subtitle=motive ) + if labels is not None: + print(labels[["geometry"]]) + labels = labels.to_crs(4326) + print(labels[["geometry"]]) + for index, row in labels.iterrows(): + if row["prominence"] == 1: + fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=24, color="black"),)) + elif row["prominence"] <3: + fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=18, color="black"),)) + else: + fig.add_trace(go.Scattergeo(lon=[row["geometry"].centroid.x], lat=[row["geometry"].centroid.y], text=[row["local_admin_unit_name"]], mode="text", textposition="top center", textfont=dict(size=12, color="black"),)) + #plt.annotate(row["local_admin_unit_name"], (row["x"], row["y"]), size=size[0], ha="center", va="center", color=color) + fig.update_geos(fitbounds="geojson", visible=False) fig.update_layout(margin=dict(l=0,r=0,t=0,b=0)) fig.show(plot_method) @@ -897,5 +920,5 @@ def public_transport_network(self, *args, **kwargs): return PublicTransportNetworkEvaluation(self).get(*args, **kwargs) - + diff --git a/mobility/experiments/gtfs_experiments.py b/mobility/experiments/gtfs_experiments.py new file mode 100644 index 00000000..84cf3d00 --- /dev/null +++ b/mobility/experiments/gtfs_experiments.py @@ -0,0 +1,15 @@ +import gtfs_kit + +k = gtfs_kit.read_feed("https://app.mecatran.com/utw/ws/gtfsfeed/static/lio?apiKey=2b160d626f783808095373766f18714901325e45&type=gtfs_lio", "m") + +trips = k.get_trips() +dates = k.get_dates() +routes = k.get_routes() +shapes = k.get_shapes(as_gdf=True) +shapes.plot() + +t = gtfs_kit.routes.build_route_timetable(k, '452', ["20251020"]) + +ts = k.compute_trip_stats(["452"]) +rs = gtfs_kit.routes.compute_route_stats(feed=k, dates=["20251020"], trip_stats_subset=ts) + diff --git a/mobility/experiments/multiple_seeds.py b/mobility/experiments/multiple_seeds.py new file mode 100644 index 00000000..e89c9ff0 --- /dev/null +++ b/mobility/experiments/multiple_seeds.py @@ -0,0 +1,96 @@ +import os +import dotenv +import mobility +import polars as pl + +dotenv.load_dotenv() + + +mobility.set_params( + # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"], + # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"] + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/test-09", + debug=True +) + +transport_zones = mobility.TransportZones("fr-09122", radius = 40, level_of_detail=1) +tz = transport_zones.get() +emp = mobility.EMPMobilitySurvey() +pop = mobility.Population(transport_zones, sample_size = 10000) +modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones)] +surveys = [emp] +motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)] + +list_seeds = [0, 42, 64, 1000, 8888, 9999, 67654, 3333333333, 56563563208387, 7345678900000743] +print(len(list_seeds)) +list_pop = [] + +pop_trips_base = mobility.PopulationTrips( + pop, + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=6, seed=0, k_mode_sequences=6) + ) +labels=pop_trips_base.get_prominent_cities() + + +for seed in list_seeds: + print("\n SEED:", seed, "\n") + pop_trips = mobility.PopulationTrips( + pop, + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=6, seed=seed, k_mode_sequences=6) + ) + cost_per_person = pop_trips.evaluate("cost_per_person", plot_delta=False, compare_with=pop_trips_base, labels=labels) + #dist_per_person = pop_trips.evaluate("distance_per_person", plot_delta=True, compare_with=pop_trips_base) + #time_per_person = pop_trips.evaluate("time_per_person", plot_delta=True, compare_with=pop_trips_base) + #ghg_per_person = pop_trips.evaluate("ghg_per_person", plot_delta=True, compare_with=pop_trips_base) + list_pop.append(pop_trips) + + +cost_all = list_pop[0].evaluate("cost_per_person").select("transport_zone_id", "cost", "n_persons").group_by("transport_zone_id").agg(pl.col("cost").sum(), pl.col("n_persons").sum()) +i = 1 +for pop in list_pop[1:10]: + cost_add = pop.evaluate("cost_per_person").select("transport_zone_id", "cost", "n_persons").group_by("transport_zone_id").agg(pl.col("cost").sum(), pl.col("n_persons").sum()) + suffix = "_" + str(i) + cost_all = cost_all.join(cost_add, on="transport_zone_id", suffix=suffix) + i += 1 + +cost_all = cost_all.sort(by="cost") + +pandas_cost=cost_all.select(pl.col("^cost.*$")).to_pandas() +pandas_cost.plot() + + +distance_all = list_pop[0].evaluate("distance_per_person").select("transport_zone_id", "distance", "n_persons").group_by("transport_zone_id").agg(pl.col("distance").sum(), pl.col("n_persons").sum()) +i = 1 +for pop in list_pop[1:10]: + distance_add = pop.evaluate("distance_per_person").select("transport_zone_id", "distance", "n_persons").group_by("transport_zone_id").agg(pl.col("distance").sum(), pl.col("n_persons").sum()) + suffix = "_" + str(i) + distance_all = distance_all.join(distance_add, on="transport_zone_id", suffix=suffix) + i += 1 + +distance_all = distance_all.sort(by="distance") + +pandas_distance=distance_all.select(pl.col("^distance.*$")).to_pandas() +pandas_distance.plot() + +time_all = list_pop[0].evaluate("time_per_person").select("transport_zone_id", "time", "n_persons").group_by("transport_zone_id").agg(pl.col("time").sum(), pl.col("n_persons").sum()) +i = 1 +for pop in list_pop[1:10]: + time_add = pop.evaluate("time_per_person").select("transport_zone_id", "time", "n_persons").group_by("transport_zone_id").agg(pl.col("time").sum(), pl.col("n_persons").sum()) + suffix = "_" + str(i) + time_all = time_all.join(time_add, on="transport_zone_id", suffix=suffix) + i += 1 + +time_all = time_all.sort(by="time") + +pandas_time=time_all.select(pl.col("^time.*$")).to_pandas() +pandas_time.plot() + +for pop in list_pop: + pop.plot_od_flows() \ No newline at end of file From 36845b457d1c78ac47bdf26401d5f82007c7ba8f Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 16:06:42 +0100 Subject: [PATCH 06/16] exclude parameters config from comparison --- mobility/model_parameters.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py index 9b9f6b6d..7bc9373d 100644 --- a/mobility/model_parameters.py +++ b/mobility/model_parameters.py @@ -49,7 +49,6 @@ def validate(self): ) def __repr__(self): - print("pop") unit_str = f" [{self.unit}]" if self.unit else "" return f"" @@ -94,7 +93,7 @@ def get_values_for_sensitivity_analyses(self, i_max=10): @dataclass class ParameterSet: - parameters : dict = field(init=False) + parameters : dict = field(init=False, compare=False) def validate(self): for param in fields(self)[1:]: From 72dd752a17390e4ba7755a3eae923e0946f50911 Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 16:30:29 +0100 Subject: [PATCH 07/16] test --- mobility/choice_models/population_trips.py | 2 +- mobility/choice_models/results.py | 2 +- mobility/model_parameters.py | 8 ++++++-- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index aafcd674..8d024559 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -106,7 +106,7 @@ def __init__( "motives": motives, "modes": modes, "surveys": surveys, - "parameters": parameters + "parameters": parameters.to_hashable_dict() } project_folder = pathlib.Path(os.environ["MOBILITY_PROJECT_DATA_FOLDER"]) diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py index e7b676f9..6673b523 100644 --- a/mobility/choice_models/results.py +++ b/mobility/choice_models/results.py @@ -54,7 +54,7 @@ def __init__( "trip_count_by_demand_group": self.trip_count_by_demand_group, "distance_per_person": self.distance_per_person, "ghg_per_person": self.ghg_per_person, - "ghg_emissions_per_trip_per_person": self.ghg_emissions_per_trip, #Test + #"ghg_emissions_per_trip_per_person": self.ghg_emissions_per_trip, #Test "time_per_person": self.time_per_person, "cost_per_person": self.cost_per_person, "immobility": self.immobility, diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py index 7bc9373d..2ae4ea29 100644 --- a/mobility/model_parameters.py +++ b/mobility/model_parameters.py @@ -1,4 +1,4 @@ -from dataclasses import dataclass, field, fields +from dataclasses import dataclass, field, fields, asdict from typing import List, Union Number = Union[int, float] @@ -104,4 +104,8 @@ def validate(self): self._validate_param_interdependency() def _validate_param_interdependency(self): - pass \ No newline at end of file + pass + + def to_hashable_dict(self) -> dict: + """Return JSON-serializable values only.""" + return asdict(self) \ No newline at end of file From f9e961b15ca59352aeb85f57ad2c638a25061667 Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 19:02:58 +0100 Subject: [PATCH 08/16] other test --- mobility/choice_models/population_trips.py | 2 +- .../population_trips_parameters.py | 20 ++++++++++++++++++ mobility/model_parameters.py | 21 ++++++++++++++----- 3 files changed, 37 insertions(+), 6 deletions(-) diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index 8d024559..933c6f3e 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -106,7 +106,7 @@ def __init__( "motives": motives, "modes": modes, "surveys": surveys, - "parameters": parameters.to_hashable_dict() + "parameters": parameters.to_dict() } project_folder = pathlib.Path(os.environ["MOBILITY_PROJECT_DATA_FOLDER"]) diff --git a/mobility/choice_models/population_trips_parameters.py b/mobility/choice_models/population_trips_parameters.py index 0f5bdfe1..d6467310 100644 --- a/mobility/choice_models/population_trips_parameters.py +++ b/mobility/choice_models/population_trips_parameters.py @@ -96,6 +96,26 @@ def __post_init__(self): "param_name": Parameter(name="Name", name_fr="Nom", description="Name of the run with those parameters", value=self.name, parameter_type=str)} + def to_dict(self): + # Return a structured description of all parameters + return { + "n_iterations": self.n_iterations, + "alpha": self.alpha, + "k_mode_sequences": self.k_mode_sequences, + "dest_prob_cutoff": self.dest_prob_cutoff, + "n_iter_per_cost_update": self.n_iter_per_cost_update, + "cost_uncertainty_sd": self.cost_uncertainty_sd, + "seed": self.seed, + "seeds": self.seeds, + "mode_sequence_search_parallel": self.mode_sequence_search_parallel, + "min_activity_time_constant": self.min_activity_time_constant, + "simulate_weekend": self.simulate_weekend, + "name": self.name, + # Add parameters in structured form + "parameters": {key: param.to_dict() for key, param in self.parameters.items()} + } + + def _validate_param_interdependency(self): if self.seeds is not None: for seed in self.seeds: diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py index 2ae4ea29..8a8d9d03 100644 --- a/mobility/model_parameters.py +++ b/mobility/model_parameters.py @@ -19,7 +19,22 @@ class Parameter: source_default: str = "" parameter_role: str = "" - + def to_dict(self): + # Convert the parameter to a dictionary with serializable values + return { + "name": self.name, + "name_fr": self.name_fr, + "value": self.value, + "description": self.description, + "parameter_type": str(self.parameter_type), # Convert type to string + "possible_values": self.possible_values, + "min_value": self.min_value, + "max_value": self.max_value, + "unit": self.unit, + "interval": self.interval, + "source_default": self.source_default, + "parameter_role": self.parameter_role, + } # def get(self): # """Return parameter value.""" @@ -105,7 +120,3 @@ def validate(self): def _validate_param_interdependency(self): pass - - def to_hashable_dict(self) -> dict: - """Return JSON-serializable values only.""" - return asdict(self) \ No newline at end of file From 1ae982fa17eb54f1760ba0034dcfd871e40d7a01 Mon Sep 17 00:00:00 2001 From: Cap Date: Mon, 22 Dec 2025 19:30:18 +0100 Subject: [PATCH 09/16] access parameters via a dict --- mobility/choice_models/destination_sequence_sampler.py | 6 +++--- mobility/choice_models/population_trips.py | 6 +++--- mobility/choice_models/state_updater.py | 4 ++-- mobility/choice_models/top_k_mode_sequence_search.py | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/mobility/choice_models/destination_sequence_sampler.py b/mobility/choice_models/destination_sequence_sampler.py index 1f206814..5abc5b70 100644 --- a/mobility/choice_models/destination_sequence_sampler.py +++ b/mobility/choice_models/destination_sequence_sampler.py @@ -57,13 +57,13 @@ def run( transport_zones, remaining_sinks, costs, - parameters.cost_uncertainty_sd + parameters["cost_uncertainty_sd"] ) dest_prob = self.get_destination_probability( utilities, motives, - parameters.dest_prob_cutoff + parameters["dest_prob_cutoff"] ) chains = ( @@ -74,7 +74,7 @@ def run( ) chains = self.spatialize_anchor_motives(chains, dest_prob, seed) - chains = self.spatialize_other_motives(chains, dest_prob, costs, parameters.alpha, seed) + chains = self.spatialize_other_motives(chains, dest_prob, costs, parameters["alpha"], seed) dest_sequences = ( chains diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index 933c6f3e..aef716d2 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -221,7 +221,7 @@ def create_and_get_asset(self): demand_groups.write_parquet(self.cache_path["demand_groups"]) - if self.parameters.simulate_weekend: + if self.parameters["simulate_weekend"]: weekend_flows, weekend_sinks, demand_groups, weekend_costs, weekend_chains = self.run_model(is_weekday=False) @@ -292,7 +292,7 @@ def run_model(self, is_weekday): remaining_sinks = sinks.clone() - for iteration in range(1, parameters.n_iterations+1): + for iteration in range(1, parameters["n_iterations"]+1): logging.info(f"Iteration n°{iteration}") @@ -343,7 +343,7 @@ def run_model(self, is_weekday): costs = self.state_updater.get_new_costs( costs, iteration, - parameters.n_iter_per_cost_update, + parameters["n_iter_per_cost_update"], current_states_steps, costs_aggregator ) diff --git a/mobility/choice_models/state_updater.py b/mobility/choice_models/state_updater.py index 46f0bff3..b7e7faeb 100644 --- a/mobility/choice_models/state_updater.py +++ b/mobility/choice_models/state_updater.py @@ -62,7 +62,7 @@ def get_new_states( motive_dur, iteration, motives, - parameters.min_activity_time_constant, + parameters["min_activity_time_constant"], tmp_folders ) @@ -73,7 +73,7 @@ def get_new_states( home_night_dur, home_motive.value_of_time_stay_home, stay_home_state, - parameters.min_activity_time_constant + parameters["min_activity_time_constant"] ) transition_prob = self.get_transition_probabilities(current_states, possible_states_utility) diff --git a/mobility/choice_models/top_k_mode_sequence_search.py b/mobility/choice_models/top_k_mode_sequence_search.py index d4ee1556..c412d82e 100644 --- a/mobility/choice_models/top_k_mode_sequence_search.py +++ b/mobility/choice_models/top_k_mode_sequence_search.py @@ -103,7 +103,7 @@ def run(self, iteration, costs_aggregator, tmp_folders, parameters): leg_modes[(from_, to_)].append(mode) - if parameters.mode_sequence_search_parallel is False: + if parameters["mode_sequence_search_parallel"] is False: logging.info("Finding probable mode sequences for the spatialized trip chains...") @@ -142,7 +142,7 @@ def run(self, iteration, costs_aggregator, tmp_folders, parameters): "python", "-u", str(resources.files('mobility') / "transport_modes" / "compute_subtour_mode_probabilities.py"), - "--k_sequences", str(parameters.k_mode_sequences), + "--k_sequences", str(parameters["k_mode_sequences"]), "--location_chains_path", str(location_chains_path), "--costs_path", str(costs_path), "--leg_modes_path", str(leg_modes_path), From b24567735120f6ef35d94f14024ba97efb7e53bb Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:15:10 +0100 Subject: [PATCH 10/16] move old classes to legacy folder --- mobility/__init__.py | 2 +- mobility/{choice_models => legacy}/destination_choice_model.py | 0 .../leisure_destination_choice_model.py | 0 .../shopping_destination_choice_model.py | 0 .../studies_destination_choice_model.py | 0 5 files changed, 1 insertion(+), 1 deletion(-) rename mobility/{choice_models => legacy}/destination_choice_model.py (100%) rename mobility/{choice_models => legacy}/leisure_destination_choice_model.py (100%) rename mobility/{choice_models => legacy}/shopping_destination_choice_model.py (100%) rename mobility/{choice_models => legacy}/studies_destination_choice_model.py (100%) diff --git a/mobility/__init__.py b/mobility/__init__.py index 65efae4c..12126389 100644 --- a/mobility/__init__.py +++ b/mobility/__init__.py @@ -33,7 +33,7 @@ from .cost_of_time_parameters import CostOfTimeParameters -from .choice_models.transport_mode_choice_model import TransportModeChoiceModel +#from .choice_models.transport_mode_choice_model import TransportModeChoiceModel from .parsers import LocalAdminUnits diff --git a/mobility/choice_models/destination_choice_model.py b/mobility/legacy/destination_choice_model.py similarity index 100% rename from mobility/choice_models/destination_choice_model.py rename to mobility/legacy/destination_choice_model.py diff --git a/mobility/choice_models/leisure_destination_choice_model.py b/mobility/legacy/leisure_destination_choice_model.py similarity index 100% rename from mobility/choice_models/leisure_destination_choice_model.py rename to mobility/legacy/leisure_destination_choice_model.py diff --git a/mobility/choice_models/shopping_destination_choice_model.py b/mobility/legacy/shopping_destination_choice_model.py similarity index 100% rename from mobility/choice_models/shopping_destination_choice_model.py rename to mobility/legacy/shopping_destination_choice_model.py diff --git a/mobility/choice_models/studies_destination_choice_model.py b/mobility/legacy/studies_destination_choice_model.py similarity index 100% rename from mobility/choice_models/studies_destination_choice_model.py rename to mobility/legacy/studies_destination_choice_model.py From 5619e90bb733695869c5ad15d96c5f18a136aa4d Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:17:26 +0100 Subject: [PATCH 11/16] WorkHomeFlows_fr class to compare with INSEE flows --- mobility/parsers/__init__.py | 5 +- mobility/parsers/work_home_flows.py | 105 +++++++++++++++++++++------- 2 files changed, 83 insertions(+), 27 deletions(-) diff --git a/mobility/parsers/__init__.py b/mobility/parsers/__init__.py index d94ae8d7..7645625b 100644 --- a/mobility/parsers/__init__.py +++ b/mobility/parsers/__init__.py @@ -1,10 +1,11 @@ from .patch_openpyxl import patch_openpyxl patch_openpyxl() -from .work_home_flows import download_work_home_flows +#from .work_home_flows import download_work_home_flows from .city_legal_population import CityLegalPopulation from .census_localized_individuals import CensusLocalizedIndividuals from .mobility_survey import MobilitySurvey from .local_admin_units import LocalAdminUnits from .jobs_active_population_distribution import JobsActivePopulationDistribution -from .jobs_active_population_flows import JobsActivePopulationFlows \ No newline at end of file +from .jobs_active_population_flows import JobsActivePopulationFlows +from .work_home_flows import WorkHomeFlows_fr \ No newline at end of file diff --git a/mobility/parsers/work_home_flows.py b/mobility/parsers/work_home_flows.py index affc6c4d..3338b9d4 100644 --- a/mobility/parsers/work_home_flows.py +++ b/mobility/parsers/work_home_flows.py @@ -1,38 +1,93 @@ +import logging +import numpy as np import pandas as pd import os from pathlib import Path +import pathlib + import requests import zipfile +from mobility.file_asset import FileAsset +from mobility.parsers.download_file import download_file -def download_work_home_flows(proxies={}): - """ - Downloads (if needed) the INSEE census mobility flows data - (hosted on data.gouv.fr by the mobility project at - https://www.data.gouv.fr/fr/datasets/r/f3f22487-22d0-45f4-b250-af36fc56ccd0) - """ - data_folder_path = ( - Path(os.path.dirname(__file__)).parents[0] / "data/insee/work_home_flows" - ) +#script +import mobility - if data_folder_path.exists() is False: - os.makedirs(data_folder_path) - # Download the raw survey data from insee.fr if needed - path = data_folder_path / "rp2019-mobpro-csv.zip" - if path.exists() is False: - # Download the zip file - r = requests.get( - url="https://www.data.gouv.fr/fr/datasets/r/f3f22487-22d0-45f4-b250-af36fc56ccd0", - proxies=proxies, - ) - with open(path, "wb") as file: - file.write(r.content) +class WorkHomeFlows_fr(FileAsset): - # Unzip the content - with zipfile.ZipFile(path, "r") as zip_ref: - zip_ref.extractall(data_folder_path) + + def __init__(self, year="2021"): + + inputs = {"year": year} + + file_name = f"insee_mobpro_{year}.parquet" + cache_path = pathlib.Path(os.environ["MOBILITY_PACKAGE_DATA_FOLDER"]) / "insee" / "flows" / file_name + + super().__init__(inputs, cache_path) + + def get_cached_asset(self) -> pd.DataFrame: - return str(data_folder_path / "FD_MOBPRO_2019.csv") + logging.info("French home-work flows already prepared. Reusing the file : " + str(self.cache_path)) + flows = pd.read_parquet(self.cache_path) + + return flows + + + def create_and_get_asset(self) -> pd.DataFrame: + """ + Parse and format french home-work flows for the given year. + + Returns: + A pandas.DataFrame giving the french home-work flows + """ + + urls ={ + "2021" : "https://www.insee.fr/fr/statistiques/fichier/8201899/base-flux-mobilite-domicile-lieu-travail-2021-csv.zip", + "2020" : "https://www.insee.fr/fr/statistiques/fichier/7630376/base-flux-mobilite-domicile-lieu-travail-2020-csv.zip", + "2019" : "https://www.insee.fr/fr/statistiques/fichier/6454112/base-csv-flux-mobilite-domicile-lieu-travail-2019.zip", + "2018" : "https://www.insee.fr/fr/statistiques/fichier/5393835/base-csv-flux-mobilite-domicile-lieu-travail-2018.zip" + } + + year = self.year + + folder = pathlib.Path(os.environ["MOBILITY_PACKAGE_DATA_FOLDER"]) / "insee" / "flows" + if folder.exists() is False: + os.mkdir(folder) + path = folder / f"insee_mobpro_{year}.zip" + download_file(urls[year], path) + + with zipfile.ZipFile(path, "r") as zip_ref: + zip_ref.extractall(folder) + + match self.year: + case "2021": + file_name= "base-flux-mobilite-domicile-lieu-travail-2021.csv" + col_name = "NBFLUX_C21_ACTOCC15P" + case "2020": + file_name= "base-flux-mobilite-domicile-lieu-travail-2020.csv" + col_name = "NBFLUX_C20_ACTOCC15P" + case "2019": + file_name= "base-flux-mobilite-domicile-lieu-travail-2019.csv" + col_name = "NBFLUX_C19_ACTOCC15P" + case "2018": + file_name= "base-flux-mobilite-domicile-lieu-travail-2018.csv" + col_name = "NBFLUX_C18_ACTOCC15P" + + flows = pd.read_csv( + folder / file_name, + sep=";", + usecols=["CODGEO", "DCLT", col_name], + dtype={"CODGEO": str, "DCLT": str, col_name: np.float32} + ) + flows.columns = ["local_admin_unit_id_from", "local_admin_unit_id_to","insee_flows"] + + flows["local_admin_unit_id_from"] = "fr-" + flows["local_admin_unit_id_from"] + flows["local_admin_unit_id_to"] = "fr-" + flows["local_admin_unit_id_to"] + + flows.to_parquet(self.cache_path) + + return flows From 2aa54eeae4e66db09f7d10ce820f5fc0227f1735 Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:25:18 +0100 Subject: [PATCH 12/16] new default parameters --- mobility/model_parameters.py | 7 ------- mobility/motives/other.py | 2 +- mobility/motives/work.py | 2 +- mobility/parsers/work_home_flows.py | 13 +++++++++++++ mobility/transport_modes/bicycle/bicycle_mode.py | 4 ++-- mobility/transport_modes/car/car_mode.py | 4 ++-- mobility/transport_modes/modal_transfer.py | 6 +++--- .../public_transport/gtfs/gtfs_data.py | 9 ++++++++- .../public_transport/public_transport_mode.py | 14 ++++++++++---- mobility/transport_modes/walk/walk_mode.py | 4 ++-- 10 files changed, 42 insertions(+), 23 deletions(-) diff --git a/mobility/model_parameters.py b/mobility/model_parameters.py index 8a8d9d03..3afba3e5 100644 --- a/mobility/model_parameters.py +++ b/mobility/model_parameters.py @@ -70,14 +70,11 @@ def __repr__(self): def get_values_for_sensitivity_analyses(self, i_max=10): value = self.value - print(self) - print(value) values = [value] if self.interval is None: raise ValueError("To run a sensitivity analysis, interval must be specified in the parameter configuration") i = 0 if self.max_value is not None: - print("hop") while i < i_max and value < self.max_value: value += self.interval values.append(round(value,3)) @@ -93,7 +90,6 @@ def get_values_for_sensitivity_analyses(self, i_max=10): if self.min_value is not None: while i < i_max and value > self.min_value: value -= self.interval - print(i, value) values.append(round(value,3)) i += 1 else: @@ -103,7 +99,6 @@ def get_values_for_sensitivity_analyses(self, i_max=10): i += 1 values.sort() - print(values) return values @dataclass @@ -112,9 +107,7 @@ class ParameterSet: def validate(self): for param in fields(self)[1:]: - print(self.parameters) param_name = "param_" + param.name - print(param_name) self.parameters[param_name].validate() self._validate_param_interdependency() diff --git a/mobility/motives/other.py b/mobility/motives/other.py index 7e461d15..cc05c349 100644 --- a/mobility/motives/other.py +++ b/mobility/motives/other.py @@ -11,7 +11,7 @@ def __init__( value_of_time: float = 10.0, saturation_fun_ref_level: float = 1.5, saturation_fun_beta: float = 4.0, - radiation_lambda: float = 0.99986, + radiation_lambda: float = 0.9999999, opportunities: pd.DataFrame = None, population: Population = None ): diff --git a/mobility/motives/work.py b/mobility/motives/work.py index e8834c8d..4b400b77 100644 --- a/mobility/motives/work.py +++ b/mobility/motives/work.py @@ -14,7 +14,7 @@ def __init__( saturation_fun_ref_level: float = 1.5, saturation_fun_beta: float = 4.0, survey_ids: List[str] = ["9.91"], - radiation_lambda: float = 0.99986, + radiation_lambda: float = 0.9999999, opportunities: pd.DataFrame = None, utilities: pd.DataFrame = None, country_utilities: Dict = None diff --git a/mobility/parsers/work_home_flows.py b/mobility/parsers/work_home_flows.py index 3338b9d4..4200329f 100644 --- a/mobility/parsers/work_home_flows.py +++ b/mobility/parsers/work_home_flows.py @@ -91,3 +91,16 @@ def create_and_get_asset(self) -> pd.DataFrame: flows.to_parquet(self.cache_path) return flows + +if __name__ == "__main__": + + + mobility.set_params( + # package_data_folder_path=os.environ["MOBILITY_PACKAGE_DATA_FOLDER"], + # project_data_folder_path=os.environ["MOBILITY_PROJECT_DATA_FOLDER"] + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/test-09", + debug=True + ) + f = WorkHomeFlows_fr(year="2018") + f1 = f.get() \ No newline at end of file diff --git a/mobility/transport_modes/bicycle/bicycle_mode.py b/mobility/transport_modes/bicycle/bicycle_mode.py index 888fed22..b8ba332f 100644 --- a/mobility/transport_modes/bicycle/bicycle_mode.py +++ b/mobility/transport_modes/bicycle/bicycle_mode.py @@ -35,8 +35,8 @@ def __init__( if generalized_cost_parameters is None: generalized_cost_parameters = GeneralizedCostParameters( - cost_constant=0.0, - cost_of_distance=0.0, + cost_constant=25.0, + cost_of_distance=0.20, cost_of_time=CostOfTimeParameters() ) diff --git a/mobility/transport_modes/car/car_mode.py b/mobility/transport_modes/car/car_mode.py index 9f51e8ed..3f868066 100644 --- a/mobility/transport_modes/car/car_mode.py +++ b/mobility/transport_modes/car/car_mode.py @@ -45,8 +45,8 @@ def __init__( if generalized_cost_parameters is None: generalized_cost_parameters = GeneralizedCostParameters( - cost_constant=0.0, - cost_of_distance=0.1, + cost_constant=2.0, + cost_of_distance=0.05, cost_of_time=CostOfTimeParameters() ) diff --git a/mobility/transport_modes/modal_transfer.py b/mobility/transport_modes/modal_transfer.py index 9568a95b..a13cfd58 100644 --- a/mobility/transport_modes/modal_transfer.py +++ b/mobility/transport_modes/modal_transfer.py @@ -16,11 +16,11 @@ class IntermodalTransfer(): # Max travel time from or to the connection nodes, estimated from the crow # fly distance around connection nodes and an average speed - max_travel_time: float # hours - average_speed: float # km/h + max_travel_time: float = 20.0 / 60.0 # hours + average_speed: float = 4.0 # km/h # Average transfer time between the two connected modes - transfer_time: float # minutes + transfer_time: float = 1.0 # minutes # Optional shortcuts to make some connections faster than average shortcuts_transfer_time: float = None # minutes diff --git a/mobility/transport_modes/public_transport/gtfs/gtfs_data.py b/mobility/transport_modes/public_transport/gtfs/gtfs_data.py index 6b30190b..6e17d0f9 100644 --- a/mobility/transport_modes/public_transport/gtfs/gtfs_data.py +++ b/mobility/transport_modes/public_transport/gtfs/gtfs_data.py @@ -63,13 +63,20 @@ def is_gtfs_file_ok(self, path): if "e8f2aceaaaa2493f6041dc7f0251f325-5d7ae44c16ad373ca1afbc4590f53256_gtfs-2015-chamonix-mobilit" in path.name: logging.info("Manual exception, GTFS not used from path", path) return False + if "datasud" in path.name: + logging.info("Manual exception, GTFS not used from path", path) + return False + if "4a590cb87669fe2bc39328652ef1d2e9_gtfs_generic_eu" in path.name: + logging.info(f"Manual exception, Flixbus GTFS not used from path {path.name}") + return False + try: with zipfile.ZipFile(path, 'r') as zip_ref: zip_contents = zip_ref.namelist() has_an_agency = "agency.txt" in zip_contents if has_an_agency: - logging.info("Downloaded file is a proper GTFS zip which contains an agency file.") + logging.debug("Downloaded file is a proper GTFS zip which contains an agency file.") else: logging.info("Downloaded file is a proper GTFS zip but does not contain an agency file, it will not be used by Mobility.") return True diff --git a/mobility/transport_modes/public_transport/public_transport_mode.py b/mobility/transport_modes/public_transport/public_transport_mode.py index 2837f046..b695b770 100644 --- a/mobility/transport_modes/public_transport/public_transport_mode.py +++ b/mobility/transport_modes/public_transport/public_transport_mode.py @@ -4,6 +4,7 @@ from mobility.transport_zones import TransportZones from mobility.transport_modes.transport_mode import TransportMode +from mobility.transport_modes.walk import WalkMode from mobility.transport_modes.public_transport.public_transport_routing_parameters import PublicTransportRoutingParameters from mobility.transport_modes.public_transport.public_transport_travel_costs import PublicTransportTravelCosts from mobility.transport_modes.public_transport.public_transport_generalized_cost import PublicTransportGeneralizedCost @@ -32,8 +33,8 @@ def __init__( transport_zones: TransportZones, first_leg_mode: TransportMode = None, last_leg_mode: TransportMode = None, - first_intermodal_transfer: IntermodalTransfer = None, - last_intermodal_transfer: IntermodalTransfer = None, + first_intermodal_transfer: IntermodalTransfer = IntermodalTransfer(), + last_intermodal_transfer: IntermodalTransfer = IntermodalTransfer(), routing_parameters: PublicTransportRoutingParameters = PublicTransportRoutingParameters(), generalized_cost_parameters: GeneralizedCostParameters = None, survey_ids: List[str] = [ @@ -44,6 +45,11 @@ def __init__( ghg_intensity: float = 0.05 ): + if first_leg_mode is None: + first_leg_mode = WalkMode(transport_zones) + if last_leg_mode is None: + last_leg_mode = WalkMode(transport_zones) + travel_costs = PublicTransportTravelCosts( transport_zones, routing_parameters, @@ -57,8 +63,8 @@ def __init__( if generalized_cost_parameters is None: generalized_cost_parameters = GeneralizedCostParameters( - cost_constant=0.0, - cost_of_distance=0.1, + cost_constant=1.0, + cost_of_distance=0.0, cost_of_time=CostOfTimeParameters() ) diff --git a/mobility/transport_modes/walk/walk_mode.py b/mobility/transport_modes/walk/walk_mode.py index 7a936603..a954c6f2 100644 --- a/mobility/transport_modes/walk/walk_mode.py +++ b/mobility/transport_modes/walk/walk_mode.py @@ -34,8 +34,8 @@ def __init__( if generalized_cost_parameters is None: generalized_cost_parameters = GeneralizedCostParameters( - cost_constant=0.0, - cost_of_distance=0.0, + cost_constant=0.5, + cost_of_distance=0.1, cost_of_time=CostOfTimeParameters() ) From ae60101644e86b0d9a460fcd99bebd3dbe9538e6 Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:29:02 +0100 Subject: [PATCH 13/16] minor improvements on PopulationTrips * bug quickfix for small radiuses in get_prominent_cities * better legend in plot_od_flows (still have to improve it) * ability to save OD flows as an image (in local folder for now, to improve) --- mobility/choice_models/population_trips.py | 25 +++++++++++++--------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/mobility/choice_models/population_trips.py b/mobility/choice_models/population_trips.py index aef716d2..f04fb7db 100644 --- a/mobility/choice_models/population_trips.py +++ b/mobility/choice_models/population_trips.py @@ -529,7 +529,7 @@ def plot_modal_share(self, zone="origin", mode="car", period="weekdays", plot=Tr def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_detail=1, n_largest=2000, color="blue", transparency=0.2, zones_color="xkcd:light grey", - labels=None, labels_size=[10, 6, 4], labels_color="black"): + labels=None, labels_size=[10, 6, 4], labels_color="black", save=False): """ Plot flows between the different zones for the given mode, motive, period and level of detail. @@ -609,8 +609,8 @@ def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_de # Put a legend for width on bottom right, title on the top x_min = float(biggest_flows[["x"]].min().iloc[0]) y_min = float(biggest_flows[["y"]].min().iloc[0]) - plt.plot([x_min, x_min+4000], [y_min, y_min], linewidth=2, color=color) - plt.text(x_min+6000, y_min-1000, "1 000", color=color) + plt.plot([x_min-10000, x_min-8000], [y_min-2000, y_min-2000], linewidth=2, color=color) + plt.text(x_min-6000, y_min-3000, "1 000", color=color) plt.title(f"{mode_name} flows between transport zones on {period}") # Draw all origin-destinations @@ -620,6 +620,10 @@ def plot_od_flows(self, mode="all", motive="all", period="weekdays", level_of_de if isinstance(labels, gpd.GeoDataFrame): self.__show_labels__(labels, labels_size, labels_color) + + if save: + text = "plot-" + str(self.parameters["seed"]) + ".png" + plt.savefig(text) plt.show() @@ -681,13 +685,14 @@ def get_prominent_cities(self, n_cities=20, n_levels=3, distance_km=2): # If an admin unit is too close to a bigger admin unit, make it less prominent for i in range(n_cities//2): - coords = flows_per_commune.loc[i, "geometry"] - geoflows["dists"] = geoflows["geometry"].distance(coords) - # Use distance_km here - geoflows.loc[ - ((geoflows["dists"] < distance_km*1000) & (geoflows.index > i)), "prominence" - ] = geoflows["prominence"] + 2 - geoflows = geoflows.sort_values(by="prominence").reset_index(drop=True) + if i < len(flows_per_commune): + coords = flows_per_commune.loc[i, "geometry"] + geoflows["dists"] = geoflows["geometry"].distance(coords) + # Use distance_km here + geoflows.loc[ + ((geoflows["dists"] < distance_km*1000) & (geoflows.index > i)), "prominence" + ] = geoflows["prominence"] + 2 + geoflows = geoflows.sort_values(by="prominence").reset_index(drop=True) # Keep only the most prominent admin units and add their centroids geoflows = geoflows[geoflows["prominence"] <= n_levels] From 8197f95a2c2adb6d2b284b1e9ccabb8d56de9c0a Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:30:25 +0100 Subject: [PATCH 14/16] MultiSeedPopulationTrips (in progress) MultiSeedPopulationTrips with some test code --- .../choice_models/population_trips_multi.py | 235 ++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 mobility/choice_models/population_trips_multi.py diff --git a/mobility/choice_models/population_trips_multi.py b/mobility/choice_models/population_trips_multi.py new file mode 100644 index 00000000..c7f1a08b --- /dev/null +++ b/mobility/choice_models/population_trips_multi.py @@ -0,0 +1,235 @@ +import os +import pathlib +import logging +import shutil +import random +import warnings + +import geopandas as gpd +import matplotlib.pyplot as plt +import polars as pl + +import plotly.express as px +import plotly.colors as pc + + +from typing import List + +from mobility.population import Population +from mobility.choice_models.population_trips import PopulationTrips +from mobility.choice_models.population_trips_parameters import PopulationTripsParameters +from mobility.choice_models.results import Results +from mobility.motives import Motive +from mobility.transport_modes.transport_mode import TransportMode +from mobility.parsers.mobility_survey import MobilitySurvey + +#New library +from PIL import Image + +# Script ones +import dotenv +import mobility + +class MultiSeedPopulationTrips: + """Class created to manage multi-seed PopulationTrips. Parameters are the same than for PopulationTrips. """ + + def __init__( + + self, + population: Population, + modes: List[TransportMode] = None, + motives: List[Motive] = None, + surveys: List[MobilitySurvey] = None, + parameters: PopulationTripsParameters = None, + seeds: List[int] = [0] + + ): + + list_population_trips = [] + + for seed in seeds: + + print("init seed", seed) + parameters_seed = parameters + parameters_seed.seed = seed + pop_seed = PopulationTrips( + population = population, + modes = modes, + motives = motives, + surveys=surveys, + parameters = parameters_seed + ) + list_population_trips.append(pop_seed) + + self.trips = list_population_trips + + def get(self): + return self.trips + + def evaluate(self, metric, plot_diff=False, plot=True, compare_with=None, variable=None, **kwargs): + if metric in ["distance_per_person", "time_per_person", "cost_per_person"]: #todo: , "ghg_per_person" + metric_per_person = metric + metric = metric.replace("_per_person", "") + + metric_all = self.trips[0].evaluate(metric_per_person).select("transport_zone_id", metric, "n_persons").group_by("transport_zone_id").agg(pl.col(metric).sum(), pl.col("n_persons").sum()) + i = 1 + for pop in self.trips[1:]: + metric_add = pop.evaluate(metric_per_person).select("transport_zone_id", metric, "n_persons").group_by("transport_zone_id").agg(pl.col(metric).sum(), pl.col("n_persons").sum()) + suffix = "_" + str(i) + metric_all = metric_all.join(metric_add, on="transport_zone_id", suffix=suffix) + i += 1 + metric_all = metric_all.sort(by=metric) + regex = "^" + metric + ".*$" + metric_all = metric_all.with_columns(metric_mean=pl.mean_horizontal(regex), + n_persons_mean=pl.mean_horizontal("^n_persons.*$")) + metric_all = metric_all.with_columns(metric_per_person=pl.col("metric_mean")/pl.col("n_persons_mean")) + + if plot_diff: + logging.info(f"Plotting {metric}, each curve being a different seed. Curves should be close.") + pandas_metric=metric_all.select(pl.col(regex)).to_pandas() + pandas_metric.plot() + + if plot: + pass #todo : add the mean and median map for the metric + + if compare_with is not None: + metric_all = metric_all.join(compare_with.select("transport_zone_id", "metric_per_person"), on="transport_zone_id", suffix="_base") + metric_all = metric_all.with_columns(delta=pl.col("metric_per_person")-pl.col("metric_per_person_base")) + + return metric_all + + elif metric == "metrics_by_variable": + print(variable) + metric_all = self.trips[0].evaluate("metrics_by_variable", variable=variable, plot=False) + i = 1 + for pop in self.trips[1:]: + metric_add = pop.evaluate("metrics_by_variable", variable=variable, plot=False).select("variable", variable, "value", "delta", "delta_relative") + suffix = "_" + str(i) + metric_all = metric_all.join(metric_add, on=("variable", variable), suffix=suffix) + i += 1 + + if plot: + comparison_plot_df = ( + metric_all + .rename({"value":"value_0"}) + .select(["variable", variable, pl.col("^value.*$")]) + .melt(["variable", variable], variable_name="value_type") + .sort(variable) + ) + + value_types = comparison_plot_df["value_type"].unique() + value_types = [v for v in value_types if v != "value_ref"] + ["value_ref"] + + + blues = pc.sequential.Blues_r + blue_values = [v for v in value_types if v != "value_ref"] + + color_map = { + vt: blues[i % len(blues)] + for i, vt in enumerate(blue_values) + } + color_map["value_ref"] = "red" + + fig = px.bar( + comparison_plot_df, + x=variable, + y="value", + color="value_type", + facet_col="variable", + barmode="group", + facet_col_spacing=0.05, + color_discrete_map=color_map, + category_orders={"value_type": value_types} + ) + fig.update_yaxes(matches=None, showticklabels=True) + fig.show("browser") + + else: + return NotImplemented + + + def plot_modal_share(self, plot_mode="all", **kwargs): + if plot_mode == "all": + logging.info("Plotting all modal shares for the given seeds") + for t in self.trips: + t.plot_modal_share(**kwargs) + else: + return NotImplemented + + def plot_od_flows(self, plot_mode="all", **kwargs): + """ + Plots flows for the given seeds. Can show them all or create a GIF. + + Parameters + ---------- + plot_mode : TYPE, optional + Plot mode : "all" shows OD flows for all seeds. "gif" creates a gif. The default is "all". + **kwargs : TYPE + Same parameters than the function plot_od_flows from PopulationTrips + """ + if plot_mode == "all": + logging.info("Plotting all OD flows for the given seeds") + for t in self.trips: + t.plot_od_flows(**kwargs) + elif plot_mode == "gif": + frames = [] + for t in self.trips: + t.plot_od_flows(save=True, **kwargs) + text = "plot-" + str(t.parameters["seed"]) + ".png" + frames.append(Image.open(text)) + + frames[0].save( + "output.gif", + format="GIF", + append_images=frames[1:], + save_all=True, + duration=400, # milliseconds per frame + loop=0 # 0 = infinite loop + ) + + def get_prominent_cities(self, **kwargs): + """Returning the cities of the first simulation as it is consistent over simulations""" + return self.trips[0].get_prominent_cities(**kwargs) + +if __name__ == "__main__": + dotenv.load_dotenv() + os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2025/12/01" + mobility.set_params( + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/test-09", + debug=False + ) + + transport_zones = mobility.TransportZones("fr-29019", radius = 40, level_of_detail=1) + emp = mobility.EMPMobilitySurvey() + pop = mobility.Population(transport_zones, sample_size = 1000) + modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), mobility.BicycleMode(transport_zones), mobility.PublicTransportMode(transport_zones)] + surveys = [emp] + motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)] + + # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) + pop_trips = MultiSeedPopulationTrips( + pop, + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=6), + seeds = [0, 42, 69, 99, 123] + ) + + # for t in pop_trips.get(): + # g = t.get() + + #labels = pop_trips.get_prominent_cities() + #pop_trips.plot_modal_share(mode="walk") + #pop_trips.plot_od_flows(mode="gif", labels=labels) + d = pop_trips.evaluate("distance_per_person") + t = pop_trips.evaluate("time_per_person") + c = pop_trips.evaluate("cost_per_person") + g = pop_trips.evaluate("ghg_per_person") + + v = pop_trips.evaluate("metrics_by_variable", variable="mode", plot=True) + v2 = pop_trips.evaluate("metrics_by_variable", variable="motive", plot=True) + v3 = pop_trips.evaluate("metrics_by_variable", variable="time_bin", plot=True) + v4 = pop_trips.evaluate("metrics_by_variable", variable="distance_bin", plot=True) + From c9c48f501a0063de0da3635b61f52975cfe977f1 Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:31:15 +0100 Subject: [PATCH 15/16] SSI in results and use of that for radius sensitivity analysis --- mobility/choice_models/results.py | 81 ++++++++++++++++++++++++++++- mobility/experiments/sensitivity.py | 78 +++++++++++++++++++++++++++ 2 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 mobility/experiments/sensitivity.py diff --git a/mobility/choice_models/results.py b/mobility/choice_models/results.py index 6673b523..122021b5 100644 --- a/mobility/choice_models/results.py +++ b/mobility/choice_models/results.py @@ -11,6 +11,8 @@ from mobility.choice_models.evaluation.routing_evaluation import RoutingEvaluation from mobility.choice_models.evaluation.public_transport_network_evaluation import PublicTransportNetworkEvaluation +from mobility.parsers.work_home_flows import WorkHomeFlows_fr + class Results: def __init__( @@ -61,7 +63,8 @@ def __init__( "car_traffic": self.car_traffic, "travel_costs": self.travel_costs, "routing": self.routing, - "public_transport_network": self.public_transport_network + "public_transport_network": self.public_transport_network, + "ssi": self.ssi } @@ -919,6 +922,82 @@ def routing(self, *args, **kwargs): def public_transport_network(self, *args, **kwargs): return PublicTransportNetworkEvaluation(self).get(*args, **kwargs) + def ssi(self, threshold=20, *args, **kwargs): + transport_zones_df = pl.from_pandas(self.transport_zones.get()[["local_admin_unit_id", "transport_zone_id"]]) + od_flows = self.weekday_states_steps.collect().join( + transport_zones_df, left_on="from", right_on="transport_zone_id") + od_flows = od_flows.join( + transport_zones_df, left_on="to", right_on="transport_zone_id", suffix=('_to')) + + od_flows_work = ( + od_flows + .filter(pl.col("motive") == "work") + .group_by("local_admin_unit_id", "local_admin_unit_id_to") + .agg(pl.col("n_persons").sum()) + .rename({"local_admin_unit_id" : "local_admin_unit_id_from"}) + ) + + whf = pl.from_pandas(WorkHomeFlows_fr().get()) + od_flows_work = od_flows_work.join(whf, on=["local_admin_unit_id_from","local_admin_unit_id_to"]) + + od_flows_work_above_threshold = ( + od_flows_work.filter((pl.col("insee_flows") > threshold) | (pl.col("n_persons") > threshold)) + ) + + # Classic SSI formula + ssi = ( + od_flows_work_above_threshold + .select( + ( + 2 * pl.min_horizontal("n_persons", "insee_flows").sum() + / (pl.col("n_persons").sum() + pl.col("insee_flows").sum()) + ).alias(f"ssi-{threshold}") + ) + ) + + + # Modified SSI formula from Liu & Yan 2020 + # https://doi.org/10.1038/s41598-020-61613-y + # modified with threshold + # removes i -> i + ssi_liu_yan_2020 = ( + od_flows_work_above_threshold + # exclude i = j and flows under threshold + .filter(pl.col("local_admin_unit_id_from") != pl.col("local_admin_unit_id_to")) + # compute local Sørensen term + .with_columns( + ( + 2 * pl.min_horizontal("n_persons", "insee_flows") + / (pl.col("n_persons") + pl.col("insee_flows")) + ).alias("local_ssi") + ) + # mean + .select( + (pl.col("local_ssi").mean()).alias("ssi_liu_yan_2020") + ) + ) + + # Other SSI from Yan et al. 2014 + # https://doi.org/10.1098/rsif.2014.0834 + # modified with threshold + ssi_yan_2014 = ( + od_flows_work_above_threshold + # compute local Sørensen term + .with_columns( + ( + 2 * pl.min_horizontal("n_persons", "insee_flows") + / (pl.col("n_persons") + pl.col("insee_flows")) + ).alias("local_ssi") + ) + # mean + .select( + (pl.col("local_ssi").mean()).alias("ssi_yan_2014") + ) + ) + + ssis = pl.concat([ssi, ssi_liu_yan_2020, ssi_yan_2014], how="horizontal") + return ssis + diff --git a/mobility/experiments/sensitivity.py b/mobility/experiments/sensitivity.py new file mode 100644 index 00000000..349fa2c0 --- /dev/null +++ b/mobility/experiments/sensitivity.py @@ -0,0 +1,78 @@ +import os +import dotenv +import logging + +import mobility + +import polars as pl + +dotenv.load_dotenv() + +os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2025/12/01" + +mobility.set_params( + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/test-09", + debug=False +) + +global_metrics = pl.DataFrame() +ssis = pl.DataFrame() +ssis200 = pl.DataFrame() + +for radius in range(32, 51, 2): + if radius not in [52]: #blockage for this one + print(f"\nRADIUS: {radius}") + transport_zones = mobility.TransportZones("fr-87085", radius = radius, level_of_detail=1) + emp = mobility.EMPMobilitySurvey() + pop = mobility.Population(transport_zones, sample_size = 1000) + modes = [mobility.CarMode(transport_zones), mobility.WalkMode(transport_zones), + mobility.BicycleMode(transport_zones), mobility.PublicTransportMode(transport_zones)] + surveys = [emp] + motives = [mobility.HomeMotive(), mobility.WorkMotive(), mobility.OtherMotive(population=pop)] + + + # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) + pop_trips = mobility.PopulationTrips( + pop, + modes, + motives, + surveys, + parameters=mobility.PopulationTripsParameters(n_iterations=4, k_mode_sequences=42) + ) + + labels=pop_trips.get_prominent_cities() + + if global_metrics.is_empty(): + global_metrics = pop_trips.evaluate("global_metrics") + else: + suffix = str(radius) + global_metrics = global_metrics.join(pop_trips.evaluate("global_metrics"), on =["country", "variable"], suffix=suffix) + + if radius % 40 == 0: + metrics_by_mode = pop_trips.evaluate("metrics_by_variable", variable="mode", plot=True) + metrics_by_motive = pop_trips.evaluate("metrics_by_variable", variable="motive", plot=True) + + # # OD flows between transport zones + pop_trips.plot_od_flows(mode="car", level_of_detail=1, labels=labels) + pop_trips.plot_od_flows(mode="walk", level_of_detail=1, labels=labels) + pop_trips.plot_od_flows(mode="bicycle", level_of_detail=1, labels=labels) + pop_trips.plot_od_flows(mode="public_transport", labels=labels) + + if ssis.is_empty(): + rad = pl.DataFrame({"radius": radius}) + ssis = pl.concat([rad, pop_trips.evaluate("ssi")], how="horizontal") + ssis200 = pl.concat([rad, pop_trips.evaluate("ssi", threshold=200)], how="horizontal") + else: + rad = pl.DataFrame({"radius": radius}) + ssis = pl.concat([ssis, pl.concat([rad, pop_trips.evaluate("ssi")], how="horizontal")]) + ssis200 = pl.concat([ssis200, pl.concat([rad, pop_trips.evaluate("ssi", threshold=200)], how="horizontal")]) + + print(global_metrics) + print(ssis) + print(ssis200) +ssis.to_pandas().plot(x="radius") +ssis200.to_pandas().plot(x="radius") + + + \ No newline at end of file From 4abff2eaa215239f3971feba2ea9df86b33831c5 Mon Sep 17 00:00:00 2001 From: Cap Date: Wed, 31 Dec 2025 17:33:40 +0100 Subject: [PATCH 16/16] last use of parameters.something ? --- mobility/choice_models/top_k_mode_sequence_search.py | 2 +- mobility/experiments/sensitivity.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/mobility/choice_models/top_k_mode_sequence_search.py b/mobility/choice_models/top_k_mode_sequence_search.py index c412d82e..e4fe85c7 100644 --- a/mobility/choice_models/top_k_mode_sequence_search.py +++ b/mobility/choice_models/top_k_mode_sequence_search.py @@ -108,7 +108,7 @@ def run(self, iteration, costs_aggregator, tmp_folders, parameters): logging.info("Finding probable mode sequences for the spatialized trip chains...") compute_subtour_mode_probabilities_serial( - parameters.k_mode_sequences, + parameters["k_mode_sequences"], unique_location_chains, costs, leg_modes, diff --git a/mobility/experiments/sensitivity.py b/mobility/experiments/sensitivity.py index 349fa2c0..b0f0d316 100644 --- a/mobility/experiments/sensitivity.py +++ b/mobility/experiments/sensitivity.py @@ -20,8 +20,8 @@ ssis = pl.DataFrame() ssis200 = pl.DataFrame() -for radius in range(32, 51, 2): - if radius not in [52]: #blockage for this one +for radius in range(36, 51, 2): + if radius not in [36]: #blockage for this one print(f"\nRADIUS: {radius}") transport_zones = mobility.TransportZones("fr-87085", radius = radius, level_of_detail=1) emp = mobility.EMPMobilitySurvey()