Remove hard-coded area-based bounds from Purpose (#459)

* Remove hard-coded area-based bounds from `Purpose`

Introduce more general sub-bounds that can define separate parameters
for arbitrary number of sub-regions instead of two

* Generalize transit trips per month accordingly

* Add helpful comment

Scripts/datatypes/purpose.py

@@ -1,6 +1,7 @@
 import numpy
 import pandas
 
+import parameters.zone as param
 from parameters.destination_choice import secondary_destination_threshold
 import models.logit as logit
 import models.generation as generation
@@ -35,25 +36,13 @@ def __init__(self, specification, zone_data, resultdata=None):
         self.dest = specification["dest"]
         self.area = specification["area"]
         self.sources = []
-        if self.area == "metropolitan":
-            l = 0
-            m = zone_data.first_surrounding_zone
-            u = zone_data.first_peripheral_zone
-        if self.area == "peripheral":
-            l = zone_data.first_peripheral_zone
-            m = None
-            u = zone_data.nr_zones
-        if self.area == "all":
-            l = 0
-            m = zone_data.first_surrounding_zone
-            u = zone_data.nr_zones
-        if self.area == "external":
-            l = zone_data.first_external_zone
-            m = None
-            u = None
-        self.bounds = slice(l, u)
-        self.lbounds = slice(l, m)
-        self.ubounds = slice(m, u)
+        zone_numbers = zone_data.zone_numbers
+        zone_intervals = param.purpose_areas[self.area]
+        self.bounds = slice(*zone_numbers.searchsorted(
+            [zone_intervals[0], zone_intervals[-1]]))
+        sub_intervals = zone_numbers[self.bounds].searchsorted(zone_intervals)
+        self.sub_bounds = [slice(sub_intervals[i-1], sub_intervals[i])
+            for i in range(1, len(sub_intervals))]
         self.zone_data = zone_data
         self.resultdata = resultdata
         self.model = None

Scripts/datatypes/tour.py

@@ -4,7 +4,7 @@
 import parameters.car as param
 import parameters.zone as zone_param
 from parameters.assignment import assignment_classes, vot_inv
-from parameters.impedance_transformation import divided_classes, trips_month
+from parameters.impedance_transformation import divided_classes, transit_trips_per_month
 
 
 class Tour:
@@ -200,14 +200,12 @@ def calc_cost(self, impedance):
 
     def _get_cost(self, impedance, mtx_type):
         """Get cost and time components from tour dest choice."""
-        if self.mode in divided_classes:
-            ass_class = "{}_{}".format(
-                self.mode, assignment_classes[self.purpose.name])
-        else:
-            ass_class = self.mode
+        demand_type = assignment_classes[self.purpose.name]
+        ass_class = ("{}_{}".format(self.mode, demand_type)
+            if self.mode in divided_classes else self.mode)
         cost = 0
         try:
-            if assignment_classes[self.purpose_name] == "work":
+            if demand_type == "work":
                 departure_imp = impedance["aht"][mtx_type][ass_class]
                 sec_dest_imp = impedance["iht"][mtx_type][ass_class]
                 return_imp = impedance["iht"][mtx_type][ass_class]
@@ -229,8 +227,10 @@ def _get_cost(self, impedance, mtx_type):
             pass
         # scale transit costs from month to day
         if self.mode == "transit" and mtx_type == "cost":
-            idx = int(self.position[0] > self.purpose.zone_data.first_surrounding_zone)
-            cost /= trips_month[ass_class][idx]
+            idx = numpy.searchsorted(
+                [bounds.stop for bounds in self.purpose.sub_bounds],
+                self.position[0], side="right")
+            cost /= transit_trips_per_month[self.purpose.area][demand_type][idx]
         return cost
 
     def __str__(self):

Scripts/models/logit.py

@@ -29,8 +29,7 @@ def __init__(self, zone_data, purpose, resultdata, is_agent_model):
         self.resultdata = resultdata
         self.purpose = purpose
         self.bounds = purpose.bounds
-        self.lbounds = purpose.lbounds
-        self.ubounds = purpose.ubounds
+        self.sub_bounds = purpose.sub_bounds
         self.zone_data = zone_data
         self.dest_exps = {}
         self.mode_exps = {}
@@ -112,13 +111,12 @@ def _add_constant(self, utility, b):
         """
         try: # If only one parameter
             utility += b
-        except ValueError: # Separate params for cap region and surrounding
-            if utility.ndim == 1: # 1-d array calculation
-                utility[self.lbounds] += b[0]
-                utility[self.ubounds] += b[1]
-            else: # 2-d matrix calculation
-                utility[self.lbounds, :] += b[0]
-                utility[self.ubounds, :] += b[1]
+        except ValueError: # Separate sub-region parameters
+            for i, bounds in enumerate(self.sub_bounds):
+                if utility.ndim == 1: # 1-d array calculation
+                    utility[bounds] += b[i]
+                else: # 2-d matrix calculation
+                    utility[bounds, :] += b[i]
 
     def _add_impedance(self, utility, impedance, b):
         """Adds simple linear impedances to utility.
@@ -139,9 +137,9 @@ def _add_impedance(self, utility, impedance, b):
         for i in b:
             try: # If only one parameter
                 utility += b[i] * impedance[i]
-            except ValueError: # Separate params for cap region and surrounding
-                utility[self.lbounds, :] += b[i][0] * impedance[i][self.lbounds, :]
-                utility[self.ubounds, :] += b[i][1] * impedance[i][self.ubounds, :]
+            except ValueError: # Separate sub-region parameters
+                for j, bounds in enumerate(self.sub_bounds):
+                    utility[bounds, :] += b[i][j] * impedance[i][bounds, :]
         return utility
 
     def _add_log_impedance(self, exps, impedance, b):
@@ -168,11 +166,10 @@ def _add_log_impedance(self, exps, impedance, b):
         for i in b:
             try: # If only one parameter
                 exps *= numpy.power(impedance[i] + 1, b[i])
-            except ValueError: # Separate params for cap region and surrounding
-                exps[self.lbounds, :] *= numpy.power(
-                    impedance[i][self.lbounds, :] + 1, b[i][0])
-                exps[self.ubounds, :] *= numpy.power(
-                    impedance[i][self.ubounds, :] + 1, b[i][1])
+            except ValueError: # Separate sub-region parameters
+                for j, bounds in enumerate(self.sub_bounds):
+                    exps[bounds, :] *= numpy.power(
+                        impedance[i][bounds, :] + 1, b[i][j])
         return exps
 
     def _add_zone_util(self, utility, b, generation=False):
@@ -196,15 +193,13 @@ def _add_zone_util(self, utility, b, generation=False):
         for i in b:
             try: # If only one parameter
                 utility += b[i] * zdata.get_data(i, self.bounds, generation)
-            except ValueError: # Separate params for cap region and surrounding
-                data_cap_region = zdata.get_data(i, self.lbounds, generation)
-                data_surrounding = zdata.get_data(i, self.ubounds, generation)
-                if utility.ndim == 1: # 1-d array calculation
-                    utility[self.lbounds] += b[i][0] * data_cap_region
-                    utility[self.ubounds] += b[i][1] * data_surrounding
-                else: # 2-d matrix calculation
-                    utility[self.lbounds, :] += b[i][0] * data_cap_region
-                    utility[self.ubounds, :] += b[i][1] * data_surrounding
+            except ValueError: # Separate sub-region parameters
+                for j, bounds in enumerate(self.sub_bounds):
+                    data = zdata.get_data(i, bounds, generation)
+                    if utility.ndim == 1: # 1-d array calculation
+                        utility[bounds] += b[i][j] * data
+                    else: # 2-d matrix calculation
+                        utility[bounds, :] += b[i][j] * data
         return utility
 
     def _add_sec_zone_util(self, utility, b, orig=None, dest=None):
@@ -344,8 +339,8 @@ def calc_individual_prob(self, mod_mode, dummy):
         try:
             self.mode_exps[mod_mode] *= numpy.exp(b)
         except ValueError:
-            self.mode_exps[mod_mode][self.lbounds] *= numpy.exp(b[0])
-            self.mode_exps[mod_mode][self.ubounds] *= numpy.exp(b[1])
+            for i, bounds in enumerate(self.sub_bounds):
+                self.mode_exps[mod_mode][bounds] *= numpy.exp(b[i])
         mode_expsum = numpy.zeros_like(self.mode_exps[mod_mode])
         for mode in self.mode_choice_param:
             mode_expsum += self.mode_exps[mode]
@@ -399,8 +394,9 @@ def calc_individual_mode_prob(self, is_car_user, zone):
             # Convert utility into euros
             money_utility = 1 / b
         except TypeError:
-            # Separate params for cap region and surrounding
-            money_utility = 1 / b[0] if self.lbounds.stop < zone else 1 / b[1]
+            # Separate sub-region parameters
+            money_utility = (1 / b[0] if self.sub_bounds[0].stop < zone
+                else 1 / b[1])
         money_utility /= self.mode_choice_param["car"]["log"]["logsum"]
         accessibility = -money_utility * logsum
         return probs, accessibility

Scripts/parameters/impedance_transformation.py

@@ -1,8 +1,18 @@
 ### IMPEDANCE TRANSFORMATION PARAMETERS ###
 
-trips_month = {
-    "transit_work": (60.0, 44.0),
-    "transit_leisure": (30.0, 30.0),
+transit_trips_per_month = {
+    "metropolitan": {
+        "work": (60.0, 44.0),
+        "leisure": (30.0, 30.0),
+    },
+    "peripheral": {
+        "work": (44.0,),
+        "leisure": (30.0,),
+    },
+    "all": {
+        "work": (60.0, 44.0),
+        "leisure": (30.0, 30.0),
+    },
 }
 
 impedance_share = {

Scripts/parameters/zone.py

@@ -93,6 +93,15 @@
         "area": "all",
     },
 )
+# Tour purpose zone intervals
+# Some demand models have separate sub-region parameters,
+# hence need sub-intervals defined.
+purpose_areas = {
+    "metropolitan": (0, 6000, 16000),
+    "peripheral": (16000, 31000),
+    "all": (0, 6000, 31000),
+    "external": (31031, 40000),
+}
 areas = {
     "helsinki_cbd": (0, 999),
     "helsinki_other": (1000, 1999),

Scripts/tests/unit/test_logit.py

@@ -45,8 +45,7 @@ class Purpose:
             },
         }
         pur.bounds = slice(0, 9)
-        pur.lbounds = slice(0, 7)
-        pur.ubounds = slice(7, 9)
+        pur.sub_bounds = [slice(0, 7), slice(7, 9)]
         pur.zone_numbers = METROPOLITAN_ZONES
         for i in ("hw", "hc", "hu", "hs", "ho"):
             pur.name = i

Scripts/transform/impedance_transformer.py

@@ -1,7 +1,7 @@
 from collections import defaultdict
 import numpy
 
-from parameters.impedance_transformation import impedance_share, divided_classes, trips_month
+import parameters.impedance_transformation as param
 from parameters.assignment import assignment_classes
 
 
@@ -32,9 +32,10 @@ def transform(self, purpose, impedance):
         cols = (purpose.bounds if purpose.name == "hoo"
             else slice(0, purpose.zone_data.nr_zones))
         day_imp = {}
+        impedance_share = param.impedance_share
         for mode in impedance_share[purpose.name]:
             day_imp[mode] = defaultdict(float)
-            if mode in divided_classes:
+            if mode in param.divided_classes:
                 ass_class = "{}_{}".format(
                     mode, assignment_classes[purpose.name])
             else:
@@ -47,9 +48,11 @@ def transform(self, purpose, impedance):
                         day_imp[mode][mtx_type] += share[0] * imp[rows, cols]
                         day_imp[mode][mtx_type] += share[1] * imp[cols, rows].T
         # transit cost to eur per day
-        transit_class = "{}_{}".format("transit", assignment_classes[purpose.name])
+        trips_month = (param.transit_trips_per_month
+            [purpose.area][assignment_classes[purpose.name]])
         trips_per_month = numpy.full_like(
-            day_imp["transit"]["cost"], trips_month[transit_class][0])
-        trips_per_month[purpose.ubounds, :] = trips_month[transit_class][1]
+            day_imp["transit"]["cost"], trips_month[0])
+        for i in range(1, len(purpose.sub_bounds)):
+            trips_per_month[purpose.sub_bounds[i], :] = trips_month[i]
         day_imp["transit"]["cost"] /= trips_per_month
         return day_imp