Fixes assignment with compressed graph and multiple classes (#223)

aequilibrae/paths/all_or_nothing.py

@@ -78,7 +78,6 @@ def execute(self):
         self.results.compact_link_loads = np.sum(self.aux_res.temp_link_loads, axis=2)
         assign_link_loads(self.results.link_loads, self.results.compact_link_loads,
                           self.results.crosswalk, self.results.cores)
-        self.results.total_flows()
         if pyqt:
             self.assignment.emit(["finished_threaded_procedure", None])
 

aequilibrae/paths/linear_approximation.py

@@ -114,7 +114,7 @@ def __init__(self, assig_spec, algorithm) -> None:
         for c in self.traffic_classes:
             r = AssignmentResults()
             r.prepare(c.graph, c.matrix)
-            self.step_direction[c.mode] = r
+            self.step_direction[c.__id__] = r
 
         if self.algorithm in ["cfw", "bfw"]:
 
@@ -123,19 +123,19 @@ def __init__(self, assig_spec, algorithm) -> None:
                 r.prepare(c.graph, c.matrix)
                 r.compact_link_loads = np.zeros([])
                 r.compact_total_link_loads = np.zeros([])
-                self.previous_step_direction[c.mode] = r
+                self.previous_step_direction[c.__id__] = r
 
                 r = AssignmentResults()
                 r.prepare(c.graph, c.matrix)
                 r.compact_link_loads = np.zeros([])
                 r.compact_total_link_loads = np.zeros([])
-                self.step_direction[c.mode] = r
+                self.step_direction[c.__id__] = r
 
                 r = AssignmentResults()
                 r.prepare(c.graph, c.matrix)
                 r.compact_link_loads = np.zeros([])
                 r.compact_total_link_loads = np.zeros([])
-                self.pre_previous_step_direction[c.mode] = r
+                self.pre_previous_step_direction[c.__id__] = r
 
     def calculate_conjugate_stepsize(self):
         self.vdf.apply_derivative(
@@ -148,7 +148,7 @@ def calculate_conjugate_stepsize(self):
         numerator = 0.0
         denominator = 0.0
         for c in self.traffic_classes:
-            stp_dir = self.step_direction[c.mode]
+            stp_dir = self.step_direction[c.__id__]
             prev_dir_minus_current_sol = np.sum(stp_dir.link_loads[:, :] - c.results.link_loads[:, :], axis=1)
             aon_minus_current_sol = np.sum(c._aon_results.link_loads[:, :] - c.results.link_loads[:, :], axis=1)
             aon_minus_prev_dir = np.sum(c._aon_results.link_loads[:, :] - stp_dir.link_loads[:, :], axis=1)
@@ -180,19 +180,18 @@ def calculate_biconjugate_direction(self):
         nu_denom = 0.0
         for c in self.traffic_classes:
             x_ = np.sum(
-                (
-                    self.step_direction[c.mode].link_loads[:, :] * self.stepsize
-                    + self.previous_step_direction[c.mode].link_loads[:, :] * (1.0 - self.stepsize)
-                    - c.results.link_loads[:, :]
-                ),
+                (self.step_direction[c.__id__].link_loads[:, :] * self.stepsize
+                 + self.previous_step_direction[c.__id__].link_loads[:, :] * (1.0 - self.stepsize)
+                 - c.results.link_loads[:, :]
+                 ),
                 axis=1,
             )
 
             y_ = np.sum(c._aon_results.link_loads[:, :] - c.results.link_loads[:, :], axis=1)
-            z_ = np.sum(self.step_direction[c.mode].link_loads[:, :] - c.results.link_loads[:, :], axis=1)
+            z_ = np.sum(self.step_direction[c.__id__].link_loads[:, :] - c.results.link_loads[:, :], axis=1)
             mu_numerator += x_ * y_
             mu_denominator += x_ * np.sum(
-                self.previous_step_direction[c.mode].link_loads - self.step_direction[c.mode].link_loads[:, :], axis=1
+                self.previous_step_direction[c.__id__].link_loads - self.step_direction[c.__id__].link_loads[:, :], axis=1
             )
             nu_nom += z_ * y_
             nu_denom += z_ * z_
@@ -224,32 +223,34 @@ def __calculate_step_direction(self):
         # 2nd iteration is a fw step. if the previous step replaced the aggregated
         # solution so far, we need to start anew.
         if (
-            (self.iter == 2)
-            or (self.stepsize == 1.0)
-            or (self.do_fw_step)
-            or (self.algorithm == "frank-wolfe")
-            or (self.algorithm == "msa")
+                (self.iter == 2)
+                or (self.stepsize == 1.0)
+                or (self.do_fw_step)
+                or (self.algorithm == "frank-wolfe")
+                or (self.algorithm == "msa")
         ):
             # logger.info("FW step")
             self.do_fw_step = False
             self.do_conjugate_step = True
             self.conjugate_stepsize = 0.0
             for c in self.traffic_classes:
                 aon_res = c._aon_results
-                stp_dir_res = self.step_direction[c.mode]
+                stp_dir_res = self.step_direction[c.__id__]
                 copy_two_dimensions(stp_dir_res.link_loads, aon_res.link_loads, self.cores)
+                stp_dir_res.total_flows()
                 if c.results.num_skims > 0:
                     copy_three_dimensions(stp_dir_res.skims.matrix_view, aon_res.skims.matrix_view, self.cores)
-                sd_flows.append(aon_res.total_link_loads * c.pce)
+                sd_flows.append(aon_res.total_link_loads)
 
         # 3rd iteration is cfw. also, if we had to reset direction search we need a cfw step before bfw
         elif (self.iter == 3) or (self.do_conjugate_step) or (self.algorithm == "cfw"):
             self.do_conjugate_step = False
             self.calculate_conjugate_stepsize()
             for c in self.traffic_classes:
-                stp_dr = self.step_direction[c.mode]
-                pre_previous = self.pre_previous_step_direction[c.mode]
+                stp_dr = self.step_direction[c.__id__]
+                pre_previous = self.pre_previous_step_direction[c.__id__]
                 copy_two_dimensions(pre_previous.link_loads, stp_dr.link_loads, self.cores)
+                pre_previous.total_flows()
                 if c.results.num_skims > 0:
                     copy_three_dimensions(pre_previous.skims.matrix_view, stp_dr.skims.matrix_view, self.cores)
 
@@ -265,18 +266,19 @@ def __calculate_step_direction(self):
                         self.conjugate_stepsize,
                         self.cores,
                     )
-
-                sd_flows.append(np.sum(stp_dr.link_loads, axis=1) * c.pce)
+                stp_dr.total_flows()
+                sd_flows.append(stp_dr.total_link_loads)
         # biconjugate
         else:
             self.calculate_biconjugate_direction()
             # deep copy because we overwrite step_direction but need it on next iteration
             for c in self.traffic_classes:
-                ppst = self.pre_previous_step_direction[c.mode]  # type: AssignmentResults
-                prev_stp_dir = self.previous_step_direction[c.mode]  # type: AssignmentResults
-                stp_dir = self.step_direction[c.mode]  # type: AssignmentResults
+                ppst = self.pre_previous_step_direction[c.__id__]  # type: AssignmentResults
+                prev_stp_dir = self.previous_step_direction[c.__id__]  # type: AssignmentResults
+                stp_dir = self.step_direction[c.__id__]  # type: AssignmentResults
 
                 copy_two_dimensions(ppst.link_loads, stp_dir.link_loads, self.cores)
+                ppst.total_flows()
                 if c.results.num_skims > 0:
                     copy_three_dimensions(ppst.skims.matrix_view, stp_dir.skims.matrix_view, self.cores)
 
@@ -288,7 +290,7 @@ def __calculate_step_direction(self):
                     self.betas,
                     self.cores,
                 )
-
+                stp_dir.total_flows()
                 if c.results.num_skims > 0:
                     triple_linear_combination_skims(
                         stp_dir.skims.matrix_view,
@@ -299,9 +301,10 @@ def __calculate_step_direction(self):
                         self.cores,
                     )
 
-                sd_flows.append(np.sum(stp_dir.link_loads, axis=1) * c.pce)
+                sd_flows.append(np.sum(stp_dir.link_loads, axis=1))
 
                 copy_two_dimensions(prev_stp_dir.link_loads, ppst.link_loads, self.cores)
+                prev_stp_dir.total_flows()
                 if c.results.num_skims > 0:
                     copy_three_dimensions(prev_stp_dir.skims.matrix_view, ppst.skims.matrix_view, self.cores)
 
@@ -323,21 +326,26 @@ def execute(self):
                 self.equilibration.emit(["iterations", self.iter])
 
             aon_flows = []
-            for c in self.traffic_classes:
+            for c in self.traffic_classes:  # type: TrafficClass
+                aggregate_link_costs(self.congested_time, c.graph.compact_cost, c.results.crosswalk)
                 aon = allOrNothing(c.matrix, c.graph, c._aon_results)
                 if pyqt:
                     aon.assignment.connect(self.signal_handler)
                 aon.execute()
-                aon_flows.append(c._aon_results.total_link_loads * c.pce)
+                c._aon_results.link_loads *= c.pce
+                c._aon_results.total_flows()
+                aon_flows.append(c._aon_results.total_link_loads)
+
             self.aon_total_flow = np.sum(aon_flows, axis=0)
+
             flows = []
             if self.iter == 1:
                 for c in self.traffic_classes:
                     copy_two_dimensions(c.results.link_loads, c._aon_results.link_loads, self.cores)
-                    copy_one_dimension(c.results.total_link_loads, c._aon_results.total_link_loads, self.cores)
+                    c.results.total_flows()
                     if c.results.num_skims > 0:
                         copy_three_dimensions(c.results.skims.matrix_view, c._aon_results.skims.matrix_view, self.cores)
-                    flows.append(c.results.total_link_loads * c.pce)
+                    flows.append(c.results.total_link_loads)
 
                 if self.algorithm == "all-or-nothing":
                     break
@@ -346,7 +354,7 @@ def execute(self):
                 self.__calculate_step_direction()
                 self.calculate_stepsize()
                 for c in self.traffic_classes:
-                    stp_dir = self.step_direction[c.mode]
+                    stp_dir = self.step_direction[c.__id__]
                     cls_res = c.results
                     linear_combination(
                         cls_res.link_loads, stp_dir.link_loads, cls_res.link_loads, self.stepsize, self.cores
@@ -360,14 +368,13 @@ def execute(self):
                             self.cores,
                         )
                     cls_res.total_flows()
-                    flows.append(cls_res.total_link_loads * c.pce)
+                    flows.append(cls_res.total_link_loads)
             self.fw_total_flow = np.sum(flows, axis=0)
 
             # Check convergence
             # This needs to be done with the current costs, and not the future ones
-            converged = False
-            if self.iter > 1:
-                converged = self.check_convergence()
+            converged = self.check_convergence() if self.iter > 1 else False
+
             self.convergence_report["iteration"].append(self.iter)
             self.convergence_report["rgap"].append(self.rgap)
             self.convergence_report["warnings"].append("; ".join(self.iteration_issue))
@@ -393,12 +400,18 @@ def execute(self):
                 *self.vdf_parameters,
                 self.cores,
             )
+
             for c in self.traffic_classes:
-                aggregate_link_costs(self.congested_time, c.graph.compact_cost, c.results.crosswalk)
-                if self.time_field in c.graph.skim_fields:
-                    idx = c.graph.skim_fields.index(self.time_field)
-                    c.graph.skims[:, idx] = self.congested_time[:]
                 c._aon_results.reset()
+                if self.time_field not in c.graph.skim_fields:
+                    continue
+                idx = c.graph.skim_fields.index(self.time_field)
+                c.graph.skims[:, idx] = self.congested_time[:]
+
+        for c in self.traffic_classes:
+            c.results.link_loads /= c.pce
+            c.results.total_flows()
+
         if (self.rgap > self.rgap_target) and (self.algorithm != "all-or-nothing"):
             logger.error(f"Desired RGap of {self.rgap_target} was NOT reached")
         logger.info(f"{self.algorithm} Assignment finished. {self.iter} iterations and {self.rgap} final gap")
@@ -426,11 +439,11 @@ def derivative_of_objective(stepsize):
                 min_res = root_scalar(derivative_of_objective, bracket=[0, 1])
                 self.stepsize = min_res.root
                 if not min_res.converged:
-                    logger.warning("Descent direction stepsize finder is not converged")
+                    logger.warning("Descent direction stepsize finder has not converged")
             else:
                 min_res = root_scalar(derivative_of_objective, 1 / self.iter)
                 if not min_res.success:
-                    logger.warning("Descent direction stepsize finder is not converged")
+                    logger.warning("Descent direction stepsize finder has not converged")
                 self.stepsize = min_res.x[0]
                 if self.stepsize <= 0.0 or self.stepsize >= 1.0:
                     raise ValueError("wrong root")
@@ -448,7 +461,7 @@ def derivative_of_objective(stepsize):
             if derivative_of_objective(0.0) < derivative_of_objective(1.0):
                 if self.algorithm == "frank-wolfe" or self.conjugate_failed:
                     msa_step = 1.0 / self.iter
-                    logger.warning(f"# Alert: Adding {msa_step} to stepsize to make it non-zero. {e.args}")
+                    logger.warning(f"# Alert: Adding {msa_step} as step size to make it non-zero. {e.args}")
                     self.stepsize = msa_step
                 else:
                     self.stepsize = 0.0

aequilibrae/paths/parallel_numpy.pyx

@@ -194,8 +194,8 @@ cpdef void copy_two_dimensions_cython(double[:, :] target,
     cdef long long l = target.shape[0]
     cdef long long k = target.shape[1]
 
-    for i in prange(l, nogil=True, num_threads=cores):
-        for j in range(k):
+    for j in range(k):
+        for i in prange(l, nogil=True, num_threads=cores):
             target[i, j] = source[i, j]
 
 

aequilibrae/paths/traffic_assignment.py

@@ -182,9 +182,9 @@ def set_classes(self, classes: List[TrafficClass]) -> None:
             classes (:obj:`List[TrafficClass]`:) List of Traffic classes for assignment
         """
 
-        ids = set([x._id for x in classes])
+        ids = set([x.__id__ for x in classes])
         if len(ids) < len(classes):
-            raise Exception("Classes need to be unique. Your list of classes has repeated items")
+            raise Exception("Classes need to be unique. Your list of classes has repeated items/IDs")
         self.classes = classes  # type: List[TrafficClass]
 
     def add_class(self, traffic_class: TrafficClass) -> None:
@@ -195,7 +195,7 @@ def add_class(self, traffic_class: TrafficClass) -> None:
             traffic_class (:obj:`TrafficClass`:) Traffic class
         """
 
-        ids = [x._id for x in self.classes if x._id == traffic_class._id]
+        ids = [x.__id__ for x in self.classes if x.__id__ == traffic_class.__id__]
         if len(ids) > 0:
             raise Exception("Traffic class already in the assignment")
 

aequilibrae/paths/traffic_class.py

@@ -16,6 +16,7 @@ class TrafficClass():
         tc = TrafficClass(graph, demand_matrix)
         tc.set_pce(1.3)
     """
+
     def __init__(self, graph: Graph, matrix: AequilibraeMatrix) -> None:
         """
         Instantiates the class
@@ -41,7 +42,7 @@ def __init__(self, graph: Graph, matrix: AequilibraeMatrix) -> None:
         self.results.reset()
         self._aon_results = AssignmentResults()
         self._aon_results.prepare(self.graph, self.matrix)
-        self._id = uuid4().hex
+        self.__id__ = uuid4().hex
 
     def set_pce(self, pce: Union[float, int]) -> None:
         """Sets Passenger Car equivalent
@@ -52,3 +53,9 @@ def set_pce(self, pce: Union[float, int]) -> None:
         if not isinstance(pce, (float, int)):
             raise ValueError('PCE needs to be either integer or float ')
         self.pce = pce
+
+    def __setattr__(self, key, value):
+
+        if key not in ['graph', 'matrix', 'pce', 'mode', 'class_flow', 'results', '_aon_results', '__id__']:
+            raise KeyError('Traffic Class does not have that element')
+        self.__dict__[key] = value

tests/aequilibrae/paths/test_traffic_assignment.py

@@ -167,14 +167,14 @@ def test_execute_and_save_results(self):
         correl = np.corrcoef(self.assigclass.results.total_link_loads, results.volume.values)[0, 1]
         self.assertLess(0.8, correl)
 
-        self.assignment.max_iter = 30
+        self.assignment.max_iter = 50
         self.assignment.set_algorithm("msa")
         self.assignment.execute()
         msa25 = self.assignment.assignment.rgap
 
         self.assigclass.results.total_flows()
         correl = np.corrcoef(self.assigclass.results.total_link_loads, results.volume)[0, 1]
-        self.assertLess(0.95, correl)
+        self.assertLess(0.98, correl)
 
         self.assignment.set_algorithm("frank-wolfe")
         self.assignment.execute()
@@ -183,15 +183,15 @@ def test_execute_and_save_results(self):
 
         self.assigclass.results.total_flows()
         correl = np.corrcoef(self.assigclass.results.total_link_loads, results.volume)[0, 1]
-        self.assertLess(0.97, correl)
+        self.assertLess(0.99, correl)
 
         self.assignment.set_algorithm("cfw")
         self.assignment.execute()
         cfw25 = self.assignment.assignment.rgap
 
         self.assigclass.results.total_flows()
         correl = np.corrcoef(self.assigclass.results.total_link_loads, results.volume)[0, 1]
-        self.assertLess(0.98, correl)
+        self.assertLess(0.999, correl)
 
         # For the last algorithm, we set skimming
         self.car_graph.set_skimming(["free_flow_time", "distance"])
@@ -204,7 +204,7 @@ def test_execute_and_save_results(self):
 
         self.assigclass.results.total_flows()
         correl = np.corrcoef(self.assigclass.results.total_link_loads, results.volume)[0, 1]
-        self.assertLess(0.99, correl)
+        self.assertLess(0.9998, correl)
 
         self.assertLess(msa25, msa10)
         self.assertLess(fw25, msa25)