Revise travel day

config/base.toml

@@ -20,7 +20,7 @@ nts_regions = [
 ]
 # nts day of the week to use
 # 1: Monday, 2: Tuesday, 3: Wednesday, 4: Thursday, 5: Friday, 6: Saturday, 7: Sunday
-nts_day_of_week = 3
+nts_days_of_week = [3]
 # what crs do we want the output to be in? (just add the number, e.g. 3857)
 output_crs = 3857
 
@@ -37,6 +37,13 @@ optional_columns = [
 ]
 n_matches = 10 # What is the maximum number of NTS matches we want for each SPC household?
 
+[feasible_assignment]
+# detour factor when converting euclidian distance to actual travel distance
+detour_factor = 1.56
+# decay rate when converting euclidian to travel distance (0.0001 is a good value)
+# actual_distance = distance * (1 + ((detour_factor - 1) * np.exp(-decay_rate * distance)))
+decay_rate = 0.0001
+
 [work_assignment]
 commute_level = "OA"
 use_percentages = true # if true, optimization problem will try to minimize percentage difference at OD level (not absolute numbers). Recommended to set it to true
@@ -48,3 +55,15 @@ max_zones = 8          # maximum number of feasible zones to include in the opti
 [postprocessing]
 pam_jitter = 30
 pam_min_duration = 10
+# for get_pt_subscription: everyone above this age has a subscription (pensioners get free travel)
+# TODO: more sophisticated approach
+pt_subscription_age = 66
+# to define if a person is a student:
+# eveyone below this age is a student
+student_age_base = 16
+# everyone below this age that has at least one "education" activity is a student
+student_age_upper = 30
+# eveyone who uses one of the modes below is classified as a passenger (isPassenger = True)
+modes_passenger =  ['car_passenger', 'taxi']
+# yearly state pension: for getting hhlIncome of pensioners
+state_pension = 11502

config/base_all_msoa.toml

@@ -0,0 +1,78 @@
+[parameters]
+seed = 0
+ # this is used to query poi data from osm and to load in SPC data
+region = "leeds"
+# how many people from the SPC do we want to run the model for? Comment out if you want to run the analysis on the entire SPC populaiton
+number_of_households = 25000
+# "OA21CD": OA level, "MSOA11CD": MSOA level
+zone_id = "MSOA21CD"
+ # Only set to true if you have travel time matrix at the level specified in boundary_geography
+travel_times = false
+boundary_geography = "MSOA"
+# NTS years to use
+nts_years = [2019, 2021, 2022]
+# NTS regions to use
+nts_regions = [
+    'Yorkshire and the Humber',
+    'North West',
+    'North East',
+    'East Midlands',
+    'West Midlands',
+    'East of England',
+    'South East',
+    'South West']
+# nts day of the week to use
+# 1: Monday, 2: Tuesday, 3: Wednesday, 4: Thursday, 5: Friday, 6: Saturday, 7: Sunday
+nts_day_of_week = 3
+# what crs do we want the output to be in? (just add the number, e.g. 3857)
+output_crs = 3857
+
+[matching]
+# for optional and required columns, see the [iterative_match_categorical](https://github.com/Urban-Analytics-Technology-Platform/acbm/blob/ca181c54d7484ebe44706ff4b43c26286b22aceb/src/acbm/matching.py#L110) function
+# Do not add any column not listed below. You can only move a column from optional to require (or vise versa)
+required_columns = [
+    "number_adults",
+    "number_children",
+    "num_pension_age",
+]
+optional_columns = [
+    "number_cars",
+    "rural_urban_2_categories",
+    "employment_status",
+    "tenure_status",
+]
+# What is the maximum number of NTS matches we want for each SPC household?
+n_matches = 10
+
+[feasible_assignment]
+# detour factor when converting euclidian distance to actual travel distance
+detour_factor = 1.56
+# decay rate when converting euclidian to travel distance (0.0001 is a good value)
+# actual_distance = distance * (1 + ((detour_factor - 1) * np.exp(-decay_rate * distance)))
+decay_rate = 0.0001
+
+[work_assignment]
+commute_level = "MSOA"
+# if true, optimization problem will try to minimize percentage difference at OD level (not absolute numbers). Recommended to set it to true
+use_percentages = true
+# weights to add for each objective in the optimization problem
+weight_max_dev = 0.2
+weight_total_dev = 0.8
+# maximum number of feasible zones to include in the optimization problem (less zones makes problem smaller - so faster, but at the cost of a better solution)
+max_zones = 10
+
+[postprocessing]
+pam_jitter = 30
+pam_min_duration = 10
+# for get_pt_subscription: everyone above this age has a subscription (pensioners get free travel)
+# TODO: more sophisticated approach
+pt_subscription_age = 66
+# to define if a person is a student:
+# eveyone below this age is a student
+student_age_base = 16
+# everyone below this age that has at least one "education" activity is a student
+student_age_upper = 30
+# eveyone who uses one of the modes below is classified as a passenger (isPassenger = True)
+modes_passenger =  ['car_passenger', 'taxi']
+# yearly state pension: for getting hhlIncome of pensioners
+state_pension = 11502

scripts/2_match_households_and_individuals.py

@@ -7,7 +7,7 @@
 from acbm.assigning.utils import cols_for_assignment_all
 from acbm.cli import acbm_cli
 from acbm.config import load_and_setup_config
-from acbm.matching import MatcherExact, match_individuals
+from acbm.matching import MatcherExact, match_individuals, match_remaining_individuals
 from acbm.preprocessing import (
     count_per_group,
     nts_filter_by_region,
@@ -16,6 +16,10 @@
     transform_by_group,
     truncate_values,
 )
+from acbm.utils import (
+    households_with_common_travel_days,
+    households_with_travel_days_in_nts_weeks,
+)
 
 
 @acbm_cli
@@ -222,23 +226,48 @@ def get_interim_path(
 
     logger.info("Filtering NTS data by specified year(s)")
 
+    logger.info(f"Total NTS households: {nts_households.shape[0]:,.0f}")
     years = config.parameters.nts_years
 
     nts_individuals = nts_filter_by_year(nts_individuals, psu, years)
     nts_households = nts_filter_by_year(nts_households, psu, years)
     nts_trips = nts_filter_by_year(nts_trips, psu, years)
 
+    logger.info(
+        f"Total NTS households (after year filtering): {nts_households.shape[0]:,.0f}"
+    )
     # #### Filter by geography
-    #
 
     regions = config.parameters.nts_regions
 
     nts_individuals = nts_filter_by_region(nts_individuals, psu, regions)
     nts_households = nts_filter_by_region(nts_households, psu, regions)
     nts_trips = nts_filter_by_region(nts_trips, psu, regions)
 
-    # Create dictionaries of key value pairs
+    logger.info(
+        f"Total NTS households (after region filtering): {nts_households.shape[0]:,.0f}"
+    )
+
+    # Ensure that the households have at least one day in `nts_days_of_week` that
+    # all household members have trips for
+    if config.parameters.common_household_day:
+        hids = households_with_common_travel_days(
+            nts_trips, config.parameters.nts_days_of_week
+        )
+    else:
+        hids = households_with_travel_days_in_nts_weeks(
+            nts_trips, config.parameters.nts_days_of_week
+        )
+
+    # Subset individuals and households given filtering of trips
+    nts_trips = nts_trips[
+        nts_trips["HouseholdID"].isin(hids)
+        & nts_trips["TravDay"].isin(config.parameters.nts_days_of_week)
+    ]
+    nts_individuals = nts_individuals[nts_individuals["HouseholdID"].isin(hids)]
+    nts_households = nts_households[nts_households["HouseholdID"].isin(hids)]
 
+    # Create dictionaries of key value pairs
     """
     guide to the dictionaries:
 
@@ -924,6 +953,19 @@ def get_interim_path(
             show_progress=True,
         )
 
+        # match remaining individuals
+        remaining_ids = spc_edited.loc[
+            ~spc_edited.index.isin(matches_ind.keys()), "id"
+        ].to_list()
+        matches_remaining_ind = match_remaining_individuals(
+            df1=spc_edited,
+            df2=nts_individuals,
+            matching_columns=["age_group", "sex"],
+            remaining_ids=remaining_ids,
+            show_progress=True,
+        )
+        matches_ind.update(matches_remaining_ind)
+
         # save random sample
         with open(
             get_interim_path("matches_ind_level_categorical_random_sample.pkl"), "wb"

scripts/3.1_assign_primary_feasible_zones.py

@@ -7,6 +7,7 @@
 from acbm.assigning.utils import (
     activity_chains_for_assignment,
     get_activities_per_zone,
+    get_chosen_day,
     intrazone_time,
     replace_intrazonal_travel_time,
     zones_to_time_matrix,
@@ -28,11 +29,15 @@ def main(config_file):
     activity_chains = activity_chains_for_assignment(config)
     logger.info("Activity chains loaded")
 
-    # Filter to a specific day of the week
     logger.info("Filtering activity chains to a specific day of the week")
-    activity_chains = activity_chains[
-        activity_chains["TravDay"] == config.parameters.nts_day_of_week
-    ]
+
+    # Generate random sample of days by household
+    get_chosen_day(config).to_parquet(
+        config.output_path / "interim" / "assigning" / "chosen_trav_day.parquet"
+    )
+
+    # Filter to chosen day
+    activity_chains = activity_chains_for_assignment(config, subset_to_chosen_day=True)
 
     # --- Study area boundaries
 
@@ -74,7 +79,11 @@ def main(config_file):
     logger.info("Creating estimated travel times matrix")
     # Create a new travel time matrix based on distances between zones
     travel_time_estimates = zones_to_time_matrix(
-        zones=boundaries, id_col=config.zone_id, time_units="m"
+        zones=boundaries,
+        id_col=config.zone_id,
+        time_units="m",
+        detour_factor=config.feasible_assignment.detour_factor,
+        decay_rate=config.feasible_assignment.decay_rate,
     )
     logger.info("Travel time estimates created")
 
@@ -203,6 +212,8 @@ def main(config_file):
         time_tolerance=config.parameters.tolerance_edu
         if config.parameters.tolerance_edu is not None
         else 0.3,
+        detour_factor=config.feasible_assignment.detour_factor,
+        decay_rate=config.feasible_assignment.decay_rate,
     )
 
     logger.info("Saving feasible zones for education activities")
@@ -230,6 +241,8 @@ def main(config_file):
         time_tolerance=config.parameters.tolerance_work
         if config.parameters.tolerance_work is not None
         else 0.3,
+        detour_factor=config.feasible_assignment.detour_factor,
+        decay_rate=config.feasible_assignment.decay_rate,
     )
 
     logger.info("Saving feasible zones for work activities")

scripts/3.2.1_assign_primary_zone_edu.py

@@ -51,11 +51,8 @@ def main(config_file):
     logger.info("Loading activity chains")
 
     activity_chains = activity_chains_for_assignment(
-        config, columns=cols_for_assignment_edu()
+        config, columns=cols_for_assignment_edu(), subset_to_chosen_day=True
     )
-    activity_chains = activity_chains[
-        activity_chains["TravDay"] == config.parameters.nts_day_of_week
-    ]
 
     logger.info("Filtering activity chains for trip purpose: education")
     activity_chains_edu = activity_chains[activity_chains["dact"] == "education"]

scripts/3.2.2_assign_primary_zone_work.py

@@ -42,11 +42,9 @@ def main(config_file):
 
     # --- Activity chains
     logger.info("Loading activity chains")
-
-    activity_chains = activity_chains_for_assignment(config, cols_for_assignment_work())
-    activity_chains = activity_chains[
-        activity_chains["TravDay"] == config.parameters.nts_day_of_week
-    ]
+    activity_chains = activity_chains_for_assignment(
+        config, cols_for_assignment_work(), subset_to_chosen_day=True
+    )
 
     logger.info("Filtering activity chains for trip purpose: work")
     activity_chains_work = activity_chains[activity_chains["dact"] == "work"]
@@ -100,7 +98,7 @@ def main(config_file):
 
         logger.info("Step 4: Filtering rows and dropping unnecessary columns")
         travel_demand_clipped = travel_demand[
-            travel_demand["Place of work indicator (4 categories) code"].isin([1, 3])
+            travel_demand["Place of work indicator (4 categories) code"].isin([3])
         ]
         travel_demand_clipped = travel_demand_clipped.drop(
             columns=[
@@ -141,7 +139,7 @@ def main(config_file):
 
         logger.info("Step 2: Filtering rows and dropping unnecessary columns")
         travel_demand_clipped = travel_demand[
-            travel_demand["Place of work indicator (4 categories) code"].isin([1, 3])
+            travel_demand["Place of work indicator (4 categories) code"].isin([3])
         ]
         travel_demand_clipped = travel_demand_clipped.drop(
             columns=[
@@ -202,7 +200,9 @@ def main(config_file):
     #### ASSIGN TO ZONE FROM FEASIBLE ZONES ####
 
     zone_assignment = WorkZoneAssignment(
-        activities_to_assign=possible_zones_work, actual_flows=travel_demand_dict_nomode
+        activities_to_assign=possible_zones_work,
+        actual_flows=travel_demand_dict_nomode,
+        scaling=config.parameters.part_time_work_prob,
     )
 
     assignments_df = zone_assignment.select_work_zone_optimization(

scripts/3.2.3_assign_secondary_zone.py

@@ -38,10 +38,7 @@ def main(config_file):
     # --- Load in the data
     logger.info("Loading: activity chains")
 
-    activity_chains = activity_chains_for_assignment(config)
-    activity_chains = activity_chains[
-        activity_chains["TravDay"] == config.parameters.nts_day_of_week
-    ]
+    activity_chains = activity_chains_for_assignment(config, subset_to_chosen_day=True)
 
     # TODO: remove obsolete comment
     # --- Add OA21CD to the data

scripts/3.3_assign_facility_all.py

@@ -307,7 +307,10 @@ def main(config_file):
             x_col="TripDisIncSW",
             y_col="length",
             x_label="Reported Travel Distance (km)",
-            y_label="Actual Distance - Euclidian (km)",
+            y_label="Actual Distance - Estimated (km)",
+            detour_factor=config.feasible_assignment.detour_factor,
+            decay_rate=config.feasible_assignment.decay_rate,
+            x_axis_max=50,
             crs=f"EPSG:{config.output_crs}",
             title_prefix=f"Scatter plot of TripDisIncSW vs. Length for {activity_type}",
             save_dir=config.output_path / "plots/assigning/",
@@ -330,8 +333,11 @@ def main(config_file):
             activity_type_col="destination activity",
             x_col="TripTotalTime",
             y_col="length",
-            x_label="Reported Travel TIme (min)",
-            y_label="Actual Distance - Euclidian (km)",
+            x_label="Reported Travel Time (min)",
+            y_label="Actual Distance - Estimated (km)",
+            detour_factor=config.feasible_assignment.detour_factor,
+            decay_rate=config.feasible_assignment.decay_rate,
+            x_axis_max=180,
             crs=f"EPSG:{config.output_crs}",
             title_prefix="Scatter plot of TripTotalTime vs. Length",
             save_dir=config.output_path / "plots/assigning/",
@@ -355,6 +361,9 @@ def main(config_file):
             y_col="time",
             x_label="Reported Travel TIme (min)",
             y_label="Modelled time (min)",
+            detour_factor=config.feasible_assignment.detour_factor,
+            decay_rate=config.feasible_assignment.decay_rate,
+            x_axis_max=180,
             crs=f"EPSG:{config.output_crs}",
             title_prefix="Scatter plot of TripTotalTime vs. Modelled time",
             save_dir=config.output_path / "plots/assigning/",

scripts/4_validation.py

@@ -29,7 +29,7 @@ def main(config_file):
     # NTS data
     legs_nts = pd.read_parquet(config.output_path / "nts_trips.parquet")
 
-    legs_nts = legs_nts[legs_nts["TravDay"] == config.parameters.nts_day_of_week]
+    legs_nts = legs_nts[legs_nts["TravDay"].isin(config.parameters.nts_days_of_week)]
 
     # Model outputs
     legs_acbm = pd.read_csv(config.output_path / "legs.csv")
@@ -217,6 +217,8 @@ def main(config_file):
         start_wkt_col="start_location_geometry_wkt",
         end_wkt_col="end_location_geometry_wkt",
         crs_epsg=config.output_crs,
+        detour_factor=1.56,
+        decay_rate=0.0001,
     )
 
     # Plot: Aggregate