From c7c0ee36d064406d4aa2413aad8233d89329bd0b Mon Sep 17 00:00:00 2001 From: Cap Date: Fri, 5 Jun 2026 17:35:51 +0200 Subject: [PATCH] add SSI computation for France --- experiments/sensitivity-transitionprob.py | 463 +++++++++++++++++ experiments/sensitivity.py | 491 ++++++++++++++++++ mobility/surveys/france/__init__.py | 3 +- mobility/surveys/france/work_home_flows.py | 106 ++++ .../trips/group_day_trips/core/results.py | 86 ++- 5 files changed, 1146 insertions(+), 3 deletions(-) create mode 100644 experiments/sensitivity-transitionprob.py create mode 100644 experiments/sensitivity.py create mode 100644 mobility/surveys/france/work_home_flows.py diff --git a/experiments/sensitivity-transitionprob.py b/experiments/sensitivity-transitionprob.py new file mode 100644 index 00000000..cfd01a17 --- /dev/null +++ b/experiments/sensitivity-transitionprob.py @@ -0,0 +1,463 @@ +import os +import dotenv + +import mobility +from mobility import CarMode, WalkMode, BicycleMode, PublicTransportMode, GeneralizedCostParameters, CostOfTimeParameters, WorkParameters, OtherParameters + + +import polars as pl +import matplotlib.pyplot as plt + +from mobility.trips.group_day_trips import Parameters + + +dotenv.load_dotenv() + +os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2026/04/30" + +mobility.set_params( + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/sensitivity/2026", + debug=False, + # r_timeout_seconds=600, + # r_max_retries=2, + # r_retry_delay_seconds=10, + # r_heartbeat_interval_seconds=30, +) + +global_metrics = pl.DataFrame() +modal_shares = pl.DataFrame() +ssis = pl.DataFrame() +ssis200 = pl.DataFrame() + + +# VARIABLES +congestion_flows_scaling_factors = [0.16, 0.18, 0.20, 0.25] +radiations = [0.999, 0.9999, 0.99999, 0.999999, 0.9999999, 0.99999999] +car_constant_costs = [0.0, 2.0, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0] +car_distance_costs = [0.0, 0.05, 0.1, 0.15, 0.2, 0.3] +walk_constant_costs = [0.0, 0.5, 1.0, 2.0] +walk_distance_costs = [0.0, 0.01, 0.02, 0.05, 0.1, 0.2] +bicycle_constant_costs = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0] +bicycle_distance_costs = walk_distance_costs +pt_constant_costs = [0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0] +pt_distance_costs = [0.0, 0.02, 0.04, 0.06, 0.1, 0.2] +costs_of_time = [i * 2 for i in range (13)] # 0 to 24 +values_of_time = costs_of_time +values_of_time_home = [i * 2 for i in range(5)] # 0 to 8 +n_iterations = range(14, 15) + +# 1. Choose just below which variable to investigate +# 2. Change the name of the legend just below +# 3. Then change the line where it is used by using the parameter 'factor' (do not forget to remove it where it was used elsewhere) +factors = n_iterations +legend = "n_iterations (congestion at every iteration)" +# legend = "radiation_lambda for both work and other motives" +# legend = "radiation_lambda for work (other=0.9999)" +# legend = "radiation_lambda for other (work=0.99999)" + +radius_exceptions = [] + +#Rennes 20- +# radiuses = range(32, 60, 12) +# insee = "fr-35238" +# radius_exceptions = [] + +# Bordeaux 20- | blocage edge 68 +insee = "fr-33063" +city_name = "Bordeaux" + +#Toulouse 20- +# insee = "fr-31555" +# city_name = "Toulouse" + +#Bayonne 20- +# insee = "fr-64102" +# city_name = "Bayonne" + +radius = 20 +transitions = [0.05, 0.1, 0.25, 0.5, 0.8, 1.0] +# speed things up when n_iterations > 10: +transitions = [0.1, 0.5, 0.8, 1.0] +#tmp +# transitions = [0.8, 1.0] + + +for transition in transitions: + for i, factor in enumerate(factors): + if radius not in radius_exceptions: + # if i >= 1: + # factor_other = factors[i-1] + print(f"\n\n\nRADIUS: {radius}, TRANSITION PROBABILITY: {transition}, FACTOR: {factor}\n") + transport_zones = mobility.TransportZones(insee, radius = radius, level_of_detail=1) + + emp = mobility.EMPMobilitySurvey() + pop = mobility.Population(transport_zones, sample_size = 1000) + car_mode = CarMode(transport_zones, + congestion=True, + congestion_flows_scaling_factor=0.2, #0.2 by default + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=4.0, #4.0 by default + cost_of_distance=0.15, #0.15 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + walk_mode = WalkMode(transport_zones, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=0.0, #0.0 by default + cost_of_distance=0.05, #0.05 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + bicycle_mode = BicycleMode(transport_zones, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=5.0, #5.0 by default + cost_of_distance=0.04, #0.04 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + pt_mode = PublicTransportMode(transport_zones, + first_leg_mode= walk_mode, last_leg_mode= walk_mode, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=1.5, #1.5 by default + cost_of_distance=0.06, #0.06 by default + cost_of_time=CostOfTimeParameters(intercept=8.0))) #8.0 by default + modes = [car_mode, walk_mode, bicycle_mode, pt_mode] + surveys = [emp] + activities = [mobility.HomeActivity(value_of_time=3.0, #3.0 by default + value_of_time_stay_home=2.0), #2.0 by default + mobility.WorkActivity(parameters=WorkParameters(radiation_lambda=0.99999, #0.99999 by default + value_of_time=14.0)), #14.0 by default + mobility.OtherActivity(population=pop, parameters=OtherParameters(radiation_lambda=0.9999, #0.9999 by default + value_of_time=12.0))] #12.0 by default + + # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) + population_trips = mobility.PopulationGroupDayTrips( + pop, + modes, + activities, + surveys, + parameters=Parameters( + n_iterations=factor, # 4 by default + n_iter_per_cost_update=5, #3 by default + mode_sequence_search_parallel=False, + transition_revision_probability=transition, # 1 b y default + use_rust_mode_sequence_search=True, + persist_iteration_artifacts=False, + + ), + ) + + # get it back + results = population_trips.weekday_run.results() + labels=results.get_prominent_cities() + + # You can get weekday plan steps to inspect them + weekday_plan_steps = population_trips.get()["weekday_plan_steps"].collect() + + + + if global_metrics.is_empty(): + rad = pl.DataFrame({"transition": [transition, transition], + "factor": [factor, factor], + "type": ["value", "value_ref"]}) + global_metrics = pl.concat([rad, population_trips.weekday_run.evaluate("global_metrics").select(["value", "value_ref"]).transpose()], how="horizontal") + else: + rad = pl.DataFrame({"transition": [transition, transition], + "factor": [factor, factor], + "type": ["value", "value_ref"]}) + global_metrics = pl.concat([global_metrics, pl.concat([rad, population_trips.weekday_run.evaluate("global_metrics").select(["value", "value_ref"]).transpose()], how="horizontal")]) + #global_metrics = global_metrics.join(population_trips.weekday_run.evaluate("global_metrics"), on =["country", "variable"], suffix=suffix) + + if modal_shares.is_empty(): + rad = pl.DataFrame({"transition": [transition for i in range(16)], + "factor": [factor for i in range(16)]}) + mp = population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False) + modal_shares = pl.concat([rad, population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False).select(["variable", "mode", "value", "value_ref"])], how="horizontal") + else: + rad = pl.DataFrame({"transition": [transition for i in range(16)], + "factor": [factor for i in range(16)]}) + local_modal_shares = pl.concat([rad, population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False).select(["variable", "mode", "value", "value_ref"])], how="horizontal") + modal_shares = pl.concat([modal_shares, local_modal_shares]) + + + if radius % 50 == 0: + metrics_by_mode = population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=True) + metrics_by_motive = population_trips.weekday_run.evaluate("metrics_by_variable", variable="motive", plot=True) + + # # OD flows between transport zones and modal shares + if radius % 12 == 0: + results.plot_od_flows(mode="car", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="walk", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="bicycle", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="public_transport", labels=labels) + results.plot_modal_share(mode="public_transport", labels=labels) + results.plot_modal_share(mode="bicycle", labels=labels) + results.plot_modal_share(mode="walk", labels=labels) + cms = results.plot_modal_share(mode="car", labels=labels) + + if radius % 20 == 0: + # Congestion (visualization possible in QGIS) + population_trips.weekday_run.evaluate("car_traffic") + + + if ssis.is_empty(): + + rad = pl.DataFrame({"transition": transition, "factor": factor, "type": "value"}) + ssis = pl.concat([rad, population_trips.weekday_run.evaluate("ssi")], how="horizontal") + ssis200 = pl.concat([rad, population_trips.weekday_run.evaluate("ssi", threshold=200)], how="horizontal") + else: + rad = pl.DataFrame({"transition": transition, "factor": factor, "type": "value"}) + ssis = pl.concat([ssis, pl.concat([rad, population_trips.weekday_run.evaluate("ssi")], how="horizontal")]) + ssis200 = pl.concat([ssis200, pl.concat([rad, population_trips.weekday_run.evaluate("ssi", threshold=200)], how="horizontal")]) + + + +global_metrics = global_metrics.rename({"column_0": "n_trips", + "column_1": "time", + "column_2": "distance"}) + +print(global_metrics) +print(ssis) +print(ssis200) + + +def plot_df(df, metrics, city_name, legend, plot_ref_values=True): + """Plot a graph with the result of the sensitivity analysis""" + + fig, axes = plt.subplots( + nrows=len(metrics), + figsize=(10, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + for ax, metric in zip(axes, metrics): + + for i, transition in enumerate(transitions): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("transition") == transition) + & (pl.col("type") == "value") + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values[metric].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"transition_prob={transition}", + ) + + if plot_ref_values: + # reference values + values_ref = ( + df.filter( + (pl.col("transition") == transition) + & (pl.col("type") == "value_ref") + ) + .sort("factor") + ) + + x_ref = values_ref["factor"].to_numpy() + y_ref = values_ref[metric].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"r={transition} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(transitions), 6), + ) + + plt.suptitle(f"Sensitivity around {city_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + +def plot_mode(df, mode, metrics, city_name, legend, plot_ref_values=True, mode_name=None): + """Plot a graph with the result of the sensitivity analysis""" + + if mode_name is None: + mode_name = mode + + fig, axes = plt.subplots( + nrows=len(metrics), + figsize=(10, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + for ax, metric in zip(axes, metrics): + + for i, radius in enumerate(transitions): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("transition") == transition) + & (pl.col("mode") == mode) + & (pl.col("variable") == metric) + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values["value"].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"t_prob={transition}", + ) + + if plot_ref_values: + # reference values + + x_ref = values["factor"].to_numpy() + y_ref = values["value_ref"].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"t_prob={transition} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(transitions), 6), + ) + + plt.suptitle(f"{city_name}. {mode_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + +def plot_all_modes(df, modes, metrics, city_name, legend, plot_ref_values=True): + """Plot a graph with the result of the sensitivity analysis""" + + fig, axess = plt.subplots( + nrows=len(modes), + ncols=len(metrics), + figsize=(20, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + + for axes, mode in zip(axess, modes): + for ax, metric in zip(axes, metrics): + + for i, transition in enumerate(transitions): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("transition") == transition) + & (pl.col("mode") == mode) + & (pl.col("variable") == metric) + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values["value"].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"t_prob={transition}", + ) + + if plot_ref_values: + # reference values + + x_ref = values["factor"].to_numpy() + y_ref = values["value_ref"].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"t_prob={transition} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(transitions), 6), + ) + + plt.suptitle(f"{city_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + + +plot_df(global_metrics, ["n_trips", "time", "distance"], city_name, legend) +plot_df(ssis, ["ssi-20", "ssi_liu_yan_2020", "ssi_yan_2014"], city_name, legend, plot_ref_values=False) +# plot_mode(modal_shares, "car", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚗 Car") +# plot_mode(modal_shares, "walk", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚶 Walk") +# plot_mode(modal_shares, "bicycle", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚲 Bicyle") +# plot_mode(modal_shares, "walk/public_transport/walk", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚎 Public transport") +plot_all_modes(modal_shares, ["car", "walk", "bicycle", "walk/public_transport/walk"], ["n_trips", "time", "distance"], city_name, legend) +# plot_df(ssis200, ["ssi-200", "ssi_liu_yan_2020", "ssi_yan_2014"], city_name, legend, plot_ref_values=False) diff --git a/experiments/sensitivity.py b/experiments/sensitivity.py new file mode 100644 index 00000000..e9f53a2c --- /dev/null +++ b/experiments/sensitivity.py @@ -0,0 +1,491 @@ +import os +import dotenv + +import mobility +from mobility import CarMode, WalkMode, BicycleMode, PublicTransportMode, GeneralizedCostParameters, CostOfTimeParameters, WorkParameters, OtherParameters + + +import polars as pl +import pandas as pd +import matplotlib.pyplot as plt + +from mobility.trips.group_day_trips import Parameters + +"""A script to explore sensitivity of variables within Mobility""" + + +dotenv.load_dotenv() + +os.environ["MOBILITY_GTFS_DOWNLOAD_DATE"] = "2026/04/30" + +mobility.set_params( + package_data_folder_path="D:/mobility-data", + project_data_folder_path="D:/sensitivity/2026", + debug=False, + # r_timeout_seconds=600, + # r_max_retries=2, + # r_retry_delay_seconds=10, + # r_heartbeat_interval_seconds=30, +) + +global_metrics = pl.DataFrame() +modal_shares = pl.DataFrame() +ssis = pl.DataFrame() +ssis200 = pl.DataFrame() + + +# VARIABLES +congestion_flows_scaling_factors = [0.16, 0.18, 0.20, 0.25] +radiations = [0.999, 0.9999, 0.99999, 0.999999, 0.9999999, 0.99999999] +car_constant_costs = [0.0, 2.0, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0] +car_distance_costs = [0.0, 0.05, 0.1, 0.15, 0.2, 0.3] +walk_constant_costs = [0.0, 0.5, 1.0, 2.0] +walk_distance_costs = [0.0, 0.01, 0.02, 0.05, 0.1, 0.2] +bicycle_constant_costs = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0] +bicycle_distance_costs = walk_distance_costs +pt_constant_costs = [0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0] +pt_distance_costs = [0.0, 0.02, 0.04, 0.06, 0.1, 0.2] +costs_of_time = [i * 2 for i in range (13)] # 0 to 24 +values_of_time = costs_of_time +values_of_time_home = [i * 2 for i in range(5)] # 0 to 8 +n_iterations = range(1, 17) +k_activity_sequences = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 42] +k_mode_sequences = [1, 2, 4, 6, 8, 10, 12, 15, 20, 42] +alphas = [0.0, 0.01, 0.05, 0.1, 0.2, 0.5, 0.8] +dest_prob_cutoffs = [0.9, 0.95, 0.99, 0.995, 0.999] +transition_logit_scales = [0.08, 0.1, 0.25, 0.5, 0.75, 1.0] # +transition_utility_pruning_deltas = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] +transition_distance_frictions = [0.0, 0.25, 0.5, 0.75, 1.0] +cost_uncertainty_sds = [0.1, 0.5, 1.0, 2.0, 4.0] + +# 1. Choose just below which variable to investigate +# 2. Change the name of the legend just below +# 3. Then change the line where it is used by using the parameter 'factor' (do not forget to remove it where it was used elsewhere) +factors = cost_uncertainty_sds +legend = "cost_uncertainty_sd" +# legend = "radiation_lambda for both work and other motives" +# legend = "radiation_lambda for work (other=0.9999)" +# legend = "radiation_lambda for other (work=0.99999)" + + +#Default +#radiuses = range(20, 62, 4) +radius_exceptions = [] + +#Rennes 20- +# radiuses = range(32, 60, 12) +# insee = "fr-35238" +# radius_exceptions = [] + +# Bordeaux 20- | blocage edge 68 +radiuses = range(20, 57, 12) +insee = "fr-33063" +radius_exceptions = [] +city_name = "Bordeaux" + +#Toulouse 20- +radiuses = range(20, 44, 12) +insee = "fr-31555" +city_name = "Toulouse" + +#Bayonne 20- +radiuses = range(20, 32, 12) +insee = "fr-64102" +city_name = "Bayonne" + +#Strasbourg +# insee = "fr-67482" +# radiuses = range(20, 45, 12) + +#Lyon 20 (32 trop de mémoire) +# radiuses = range(20, 57, 12) +# insee = "fr-69381" + + +for radius in radiuses: + for i, factor in enumerate(factors): + if radius not in radius_exceptions: + # if i >= 1: + # factor_other = factors[i-1] + print(f"\nRADIUS: {radius}, FACTOR: {factor}\n\n") + transport_zones = mobility.TransportZones(insee, radius = radius, level_of_detail=1) + + emp = mobility.EMPMobilitySurvey() + pop = mobility.Population(transport_zones, sample_size = 1000) + car_mode = CarMode(transport_zones, + congestion=True, + congestion_flows_scaling_factor=0.2, #0.2 by default + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=4.0, #4.0 by default + cost_of_distance=0.15, #0.15 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + walk_mode = WalkMode(transport_zones, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=0.0, #0.0 by default + cost_of_distance=0.05, #0.05 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + bicycle_mode = BicycleMode(transport_zones, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=5.0, #5.0 by default + cost_of_distance=0.04, #0.04 by default + cost_of_time=CostOfTimeParameters(intercept=20.0))) #20.0 by default + pt_mode = PublicTransportMode(transport_zones, + first_leg_mode= walk_mode, last_leg_mode= walk_mode, + generalized_cost_parameters=GeneralizedCostParameters(cost_constant=1.5, #1.5 by default + cost_of_distance=0.06, #0.06 by default + cost_of_time=CostOfTimeParameters(intercept=8.0))) #8.0 by default + modes = [car_mode, walk_mode, bicycle_mode, pt_mode] + surveys = [emp] + activities = [mobility.HomeActivity(value_of_time=3.0, #3.0 by default + value_of_time_stay_home=2.0), #2.0 by default + mobility.WorkActivity(parameters=WorkParameters(radiation_lambda=0.99999, #0.99999 by default + value_of_time=14.0)), #14.0 by default + mobility.OtherActivity(population=pop, parameters=OtherParameters(radiation_lambda=0.9999, #0.9999 by default + value_of_time=12.0))] #12.0 by default + + # Simulating the trips for this population for three modes : car, walk and bicyle, and only home and work motives (OtherMotive is mandatory) + population_trips = mobility.PopulationGroupDayTrips( + pop, + modes, + activities, + surveys, + parameters=Parameters( + n_iterations=4, # 4 by default + n_iter_per_cost_update=3, #3 by default + mode_sequence_search_parallel=False, + transition_revision_probability=0.5, # 0.5, new default + use_rust_mode_sequence_search=True, + persist_iteration_artifacts=True, + k_activity_sequences=10, #10, new default + k_mode_sequences=10, #10, new default + alpha = 0.01, + dest_prob_cutoff=0.99, + transition_logit_scale=1.0, + transition_utility_pruning_delta=3.0, + transition_distance_friction=1.0, + cost_uncertainty_sd=factor + ), + ) + + # get it back + results = population_trips.weekday_run.results() + labels=results.get_prominent_cities() + + # You can get weekday plan steps to inspect them + weekday_plan_steps = population_trips.get()["weekday_plan_steps"].collect() + + if global_metrics.is_empty(): + rad = pl.DataFrame({"radius": [radius, radius], + "factor": [factor, factor], + "type": ["value", "value_ref"]}) + global_metrics = pl.concat([rad, population_trips.weekday_run.evaluate("global_metrics").select(["value", "value_ref"]).transpose()], how="horizontal") + else: + rad = pl.DataFrame({"radius": [radius, radius], + "factor": [factor, factor], + "type": ["value", "value_ref"]}) + global_metrics = pl.concat([global_metrics, pl.concat([rad, population_trips.weekday_run.evaluate("global_metrics").select(["value", "value_ref"]).transpose()], how="horizontal")]) + #global_metrics = global_metrics.join(population_trips.weekday_run.evaluate("global_metrics"), on =["country", "variable"], suffix=suffix) + + if modal_shares.is_empty(): + rad = pl.DataFrame({"radius": [radius for i in range(16)], + "factor": [factor for i in range(16)]}) + mp = population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False) + modal_shares = pl.concat([rad, population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False).select(["variable", "mode", "value", "value_ref"])], how="horizontal") + else: + rad = pl.DataFrame({"radius": [radius for i in range(16)], + "factor": [factor for i in range(16)]}) + local_modal_shares = pl.concat([rad, population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=False).select(["variable", "mode", "value", "value_ref"])], how="horizontal") + modal_shares = pl.concat([modal_shares, local_modal_shares]) + + + if radius % 50 == 0: + metrics_by_mode = population_trips.weekday_run.evaluate("metrics_by_variable", variable="mode", plot=True) + metrics_by_motive = population_trips.weekday_run.evaluate("metrics_by_variable", variable="motive", plot=True) + + # # OD flows between transport zones and modal shares + if radius % 12 == 0: + results.plot_od_flows(mode="car", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="walk", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="bicycle", level_of_detail=1, labels=labels) + results.plot_od_flows(mode="public_transport", labels=labels) + results.plot_modal_share(mode="public_transport", labels=labels) + results.plot_modal_share(mode="bicycle", labels=labels) + results.plot_modal_share(mode="walk", labels=labels) + cms = results.plot_modal_share(mode="car", labels=labels) + # Congestion (visualization possible in QGIS) + + + if factor == factors[-1]: + population_trips.weekday_run.evaluate("car_traffic") + + + if ssis.is_empty(): + + rad = pl.DataFrame({"radius": radius, "factor": factor, "type": "value"}) + ssis = pl.concat([rad, population_trips.weekday_run.evaluate("ssi")], how="horizontal") + ssis200 = pl.concat([rad, population_trips.weekday_run.evaluate("ssi", threshold=200)], how="horizontal") + else: + rad = pl.DataFrame({"radius": radius, "factor": factor, "type": "value"}) + ssis = pl.concat([ssis, pl.concat([rad, population_trips.weekday_run.evaluate("ssi")], how="horizontal")]) + ssis200 = pl.concat([ssis200, pl.concat([rad, population_trips.weekday_run.evaluate("ssi", threshold=200)], how="horizontal")]) + + + +global_metrics = global_metrics.rename({"column_0": "n_trips", + "column_1": "time", + "column_2": "distance"}) + +print(global_metrics) +print(ssis) +print(ssis200) + + +def plot_df(df, metrics, city_name, legend, plot_ref_values=True): + """Plot a graph with the result of the sensitivity analysis""" + + fig, axes = plt.subplots( + nrows=len(metrics), + figsize=(10, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + for ax, metric in zip(axes, metrics): + + for i, radius in enumerate(radiuses): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("radius") == radius) + & (pl.col("type") == "value") + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values[metric].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"r={radius}", + ) + + if plot_ref_values: + # reference values + values_ref = ( + df.filter( + (pl.col("radius") == radius) + & (pl.col("type") == "value_ref") + ) + .sort("factor") + ) + + x_ref = values_ref["factor"].to_numpy() + y_ref = values_ref[metric].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"r={radius} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(radiuses), 6), + ) + + plt.suptitle(f"Sensitivity around {city_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + +def plot_mode(df, mode, metrics, city_name, legend, plot_ref_values=True, mode_name=None): + """Plot a graph with the result of the sensitivity analysis""" + + if mode_name is None: + mode_name = mode + + fig, axes = plt.subplots( + nrows=len(metrics), + figsize=(10, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + for ax, metric in zip(axes, metrics): + + for i, radius in enumerate(radiuses): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("radius") == radius) + & (pl.col("mode") == mode) + & (pl.col("variable") == metric) + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values["value"].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"r={radius}", + ) + + if plot_ref_values: + # reference values + + x_ref = values["factor"].to_numpy() + y_ref = values["value_ref"].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"r={radius} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(radiuses), 6), + ) + + plt.suptitle(f"{city_name}. {mode_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + +def plot_all_modes(df, modes, metrics, city_name, legend, plot_ref_values=True): + """Plot a graph with the result of the sensitivity analysis""" + + fig, axess = plt.subplots( + nrows=len(modes), + ncols=len(metrics), + figsize=(20, 4 * len(metrics)), + sharex=True, + ) + + cmap = plt.get_cmap("tab10") + + print(axess) + + for axes, mode in zip(axess, modes): + print(axes) + for ax, metric in zip(axes, metrics): + print(ax) + + for i, radius in enumerate(radiuses): + + color = cmap(i) + + # values varying with factor + values = ( + df.filter( + (pl.col("radius") == radius) + & (pl.col("mode") == mode) + & (pl.col("variable") == metric) + ) + .sort("factor") + ) + + x = values["factor"].to_numpy() + y = values["value"].to_numpy() + + ax.plot( + x, + y, + marker="o", + color=color, + label=f"r={radius}", + ) + + if plot_ref_values: + # reference values + + x_ref = values["factor"].to_numpy() + y_ref = values["value_ref"].to_numpy() + + ax.plot( + x_ref, + y_ref, + marker="o", + color=color, + linestyle="--", + alpha=0.7, + label=f"r={radius} (ref)", + ) + + + ax.set_title(metric) + ax.set_ylabel("value") + ax.grid(alpha=0.3) + + axes[-1].set_xlabel(legend) + #axes[-1].set_xscale("logit") # only for radiation_lambda + + # single legend + handles, labels = axes[0].get_legend_handles_labels() + fig.legend( + handles, + labels, + loc="upper right", + ncol=min(len(radiuses), 6), + ) + + plt.suptitle(f"{city_name}", x=0.10, size='xx-large', ha='right', va='center') + # plt.tight_layout(rect=[0, 0, 1, 0.95]) + plt.tight_layout(rect=[0, 0, 1, 1]) + plt.show() + + +plot_df(global_metrics, ["n_trips", "time", "distance"], city_name, legend) +plot_df(ssis, ["ssi-20", "ssi_liu_yan_2020", "ssi_yan_2014"], city_name, legend, plot_ref_values=False) +# plot_mode(modal_shares, "car", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚗 Car") +# plot_mode(modal_shares, "walk", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚶 Walk") +# plot_mode(modal_shares, "bicycle", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚲 Bicyle") +# plot_mode(modal_shares, "walk/public_transport/walk", ["n_trips", "time", "distance"], city_name, legend, mode_name="🚎 Public transport") +plot_all_modes(modal_shares, ["car", "walk", "bicycle", "walk/public_transport/walk"], ["n_trips", "time", "distance"], city_name, legend) +# plot_df(ssis200, ["ssi-200", "ssi_liu_yan_2020", "ssi_yan_2014"], city_name, legend, plot_ref_values=False) diff --git a/mobility/surveys/france/__init__.py b/mobility/surveys/france/__init__.py index 91c4c486..9cf83697 100644 --- a/mobility/surveys/france/__init__.py +++ b/mobility/surveys/france/__init__.py @@ -1,2 +1,3 @@ from .emp import EMPMobilitySurvey -from .entd import ENTDMobilitySurvey \ No newline at end of file +from .entd import ENTDMobilitySurvey +from .work_home_flows import WorkHomeFlows_fr \ No newline at end of file diff --git a/mobility/surveys/france/work_home_flows.py b/mobility/surveys/france/work_home_flows.py new file mode 100644 index 00000000..50dcfb45 --- /dev/null +++ b/mobility/surveys/france/work_home_flows.py @@ -0,0 +1,106 @@ +import logging +import numpy as np +import pandas as pd +import os +from pathlib import Path +import pathlib + +import requests +import zipfile + +from mobility.runtime.assets.file_asset import FileAsset +from mobility.runtime.io.download_file import download_file + + +#script +import mobility + + + +class WorkHomeFlows_fr(FileAsset): + + + 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: + + logging.debug("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 + +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/trips/group_day_trips/core/results.py b/mobility/trips/group_day_trips/core/results.py index faba9545..92de2a76 100644 --- a/mobility/trips/group_day_trips/core/results.py +++ b/mobility/trips/group_day_trips/core/results.py @@ -17,6 +17,8 @@ from ..evaluation.travel_costs_evaluation import TravelCostsEvaluation from ..transitions.transition_metrics import state_waterfall as _state_waterfall +from ....surveys.france import WorkHomeFlows_fr + class RunResults: """Run-scoped analysis helper for one day-type output set.""" @@ -69,6 +71,7 @@ def __init__( "travel_costs": self.travel_costs, "routing": self.routing, "public_transport_network": self.public_transport_network, + "ssi": self.ssi, } @property @@ -787,8 +790,8 @@ def plot_od_flows( 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-6500, y_min-2500, "1 000", color=color) plt.title(f"{mode_name} flows between transport zones on {self.period}") for _, row in biggest_flows.iterrows(): @@ -835,6 +838,85 @@ def get_prominent_cities(self, n_cities=20, n_levels=3, distance_km=2): xy_coords = geoflows["geometry"].centroid.get_coordinates() return geoflows.merge(xy_coords, left_index=True, right_index=True) + + 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.plan_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("activity") == "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 + + @staticmethod def _show_labels(labels, size, color): """Annotate a matplotlib axes with place labels."""