Add NSGA2 and NRMSE

.gitignore

@@ -4,6 +4,8 @@
 
 ## private stuff
 .results/
+test_data/result
+test_data/01_WorldPositions_selection.csv
 .tmp/
 dist/
 ## data

carmodel_calibration/control_program/calibration_analysis.py

@@ -93,6 +93,7 @@ def _get_weighted_errors(param_sets, sim_result, ident, results_chunk,
         weighted_error = results_chunk["weightedError"].values[jdx-1]
         if weighted_error == 0 or force_simulated:
             weighted_error = objective_function(gt, pred)
+        weighted_error = np.sum(weighted_error)
         weighted_errors.append(weighted_error)
     return weighted_errors
 
@@ -167,6 +168,7 @@ def _plot_single(identification, results_chunk, simulation_result,
             pred, gt = _get_results(
                 {1: simulation_result}, identification, jdx-1)
             weighted_error = results_chunk["weightedError"].values[jdx-1]
+        weighted_error = np.sum(weighted_error)
         time = pred["time"].values
         covered_distance = pred["coveredDistanceFollower"].values
         distance = pred["distance"].values

carmodel_calibration/optimization.py

@@ -6,15 +6,18 @@
 
 import numpy as np
 import pandas as pd
+from pymoo.core.problem import Problem
 
 # from carmodel_calibration.exceptions import OptimizationFailed
 from carmodel_calibration.fileaccess.parameter import ModelParameters
 from carmodel_calibration.sumo.sumo_project import SumoProject
 
 _LOGGER = logging.getLogger(__name__)
 
+
 def measure_of_performance_factory(objectives=["distance"],
-                                   weights=None, gof="rmse", **__):
+                                   weights=None, gof="rmse", reduce_mmop2smop=True,
+                                   **__):
     """
     this function will return the appropiate objective function handle
     depending on the desired MOP and the calibration object
@@ -58,6 +61,7 @@ def measure_of_performance_factory(objectives=["distance"],
     objectives = weights_renamed
     for weights_value, objective in zip(weights_values, objectives):
         weights[objective] = weights_value
+
     def rmse(ground_truth, prediction):
         """
         RMSE(v)
@@ -70,7 +74,31 @@ def rmse(ground_truth, prediction):
                 prediction[sop].values
                 - ground_truth[sop].values))
             sums[idx] = np.sqrt(sums[idx] / len(ground_truth)) * weights[sop]
-        return np.sum(sums)
+        if reduce_mmop2smop:
+            return np.sum(sums)
+        else:
+            return sums
+
+    def nrmse(ground_truth, prediction):
+        """
+        NRMSE(v)
+        see 'About calibration of car-following dynamics of automated and
+        human-driven vehicles: Methodology, guidelines and codes'
+        """
+        sums = np.zeros(len(objectives))
+        for idx, sop in enumerate(objectives):
+            rmse_error = np.sum(np.square(
+                prediction[sop].values
+                - ground_truth[sop].values))
+            rmse_error = np.sqrt(rmse_error / len(ground_truth))
+            gt_root = np.sqrt(np.sum(np.square(ground_truth[sop].values))
+                              / len(ground_truth))
+            sums[idx] = rmse_error / (gt_root) * weights[sop]
+        if reduce_mmop2smop:
+            return np.sum(sums)
+        else:
+            return sums
+
     def rmspe(ground_truth, prediction):
         """
         RMSPE(v)
@@ -84,7 +112,11 @@ def rmspe(ground_truth, prediction):
                  - ground_truth[sop].values)
                 / ground_truth[sop].values))
             sums[idx] = np.sqrt(sums[idx] / len(ground_truth)) * weights[sop]
-        return np.sum(sums)
+        if reduce_mmop2smop:
+            return np.sum(sums)
+        else:
+            return sums
+
     def theils_u(ground_truth, prediction):
         """
         theils U(v)
@@ -102,16 +134,23 @@ def theils_u(ground_truth, prediction):
             gt_root = np.sqrt(np.sum(np.square(ground_truth[sop].values))
                               / len(ground_truth))
             sums[idx] = rmse_error / (sim_root + gt_root) * weights[sop]
-        return np.sum(sums)
+        if reduce_mmop2smop:
+            return np.sum(sums)
+        else:
+            return sums
 
     def model_output(_, prediction):
         """only sum model outputs"""
         sums = np.zeros(len(objectives))
         for idx, sop in enumerate(objectives):
             sums[idx] = prediction[sop].values[-1] * weights[sop]
-        return np.sum(sums)
-    gof_handles = {"rmse": rmse, "rmsep": rmspe, "modelOutput": model_output,
-                   "theils_u": theils_u}
+        if reduce_mmop2smop:
+            return np.sum(sums)
+        else:
+            return sums
+
+    gof_handles = {"rmse": rmse, "nrmse": nrmse, "rmsep": rmspe,
+                   "modelOutput": model_output, "theils_u": theils_u}
 
     def get_weighted_error(ground_truth, prediction):
         """calculate weighted error on case1"""
@@ -124,6 +163,7 @@ def get_weighted_error(ground_truth, prediction):
 
     return get_weighted_error
 
+
 def factory_wrapper(factory_kwargs):
     """invokes factory from results data"""
     if isinstance(factory_kwargs, pd.DataFrame):
@@ -142,31 +182,33 @@ def factory_wrapper(factory_kwargs):
     else:
         raise TypeError(
             "`factory_kwargs` must either be of type dict or pandas DattaFrame"
-            )
+        )
+
 
 def _get_results(simulation_results, identification, lane):
     """returns the simulation result for specific identification"""
     edge_name = f"B{lane}A{lane}"
     condition = (
-        (simulation_results[1]["prediction"]["edgeName"]==edge_name)
+        (simulation_results[1]["prediction"]["edgeName"] == edge_name)
     )
     prediction_result = simulation_results[1]["prediction"][condition]
     leader = identification[0]
     recording_id = identification[2]
     follower_condition = (
-        simulation_results[1]["ground_truth"]["follower"]==identification[1])
+        simulation_results[1]["ground_truth"]["follower"] == identification[1])
     condition = (
-        (simulation_results[1]["ground_truth"]["leader"]==leader)
+        (simulation_results[1]["ground_truth"]["leader"] == leader)
         & (follower_condition)
         & (simulation_results[1]["ground_truth"]
-            ["recordingId"]==recording_id)
+            ["recordingId"] == recording_id)
         & (simulation_results[1]["ground_truth"]
-            ["counter"]==lane)
+            ["counter"] == lane)
     )
     ground_truth_result = (
         simulation_results[1]["ground_truth"][condition])
     return prediction_result, ground_truth_result
 
+
 def _run_sumo(identification, sumo_interface, param_sets, project_path, model, timestep):
     cfmodels = []
     route_count = SumoProject.get_number_routes(
@@ -177,7 +219,7 @@ def _run_sumo(identification, sumo_interface, param_sets, project_path, model, t
         sumo_interface.start_simulation_module()
     for idx, param_set in enumerate(param_sets):
         cfmodel = ModelParameters.create_parameter_set(f"set{idx}.json", model,
-                                                        **param_set)
+                                                       **param_set)
         cfmodels.append(cfmodel)
     SumoProject.write_followers_leader(
         project_path / "calibration_routes.rou.xml",
@@ -186,13 +228,15 @@ def _run_sumo(identification, sumo_interface, param_sets, project_path, model, t
         identification=identification)
     return simulation_results
 
+
 def _calculate_performance(idx, simulation_results, identification,
-                          objective_function):
+                           objective_function):
     prediction, ground_truth = _get_results(
-            simulation_results, identification, idx)
+        simulation_results, identification, idx)
     objective_function = factory_wrapper(objective_function)
     return objective_function(ground_truth, prediction)
 
+
 def _vectorized_target(params, *data):
     # params.shape = (Number params, number of solutions)
     params = np.atleast_2d(params).T
@@ -206,7 +250,7 @@ def _vectorized_target(params, *data):
     # keys of the paramters that are optimized
     param_names = data["param_names"]
     param_sets = []
-    if params.shape[1]==1:
+    if params.shape[1] == 1:
         params = params.T
     for solution in params:
         params_dict = data["default_parameters"].copy()
@@ -215,22 +259,36 @@ def _vectorized_target(params, *data):
         param_sets.append(params_dict)
     simulation_results = _run_sumo(
         identification, sumo_interface, param_sets, project_path, model, timestep)
-    with Pool(os.cpu_count()//2) as pool:
+    with Pool(os.cpu_count() // 2) as pool:
         results = []
         for idx in range(len(params)):
             results.append((idx, simulation_results, identification,
                             objective_function))
         performance = list(pool.starmap(_calculate_performance, results))
-    #performance = []
-    #for idx in range(len(params)):
-    #    performance.append(_calculate_performance(idx, simulation_results, identification,
-    #                    objective_function))
+    performance = []
+    for idx in range(len(params)):
+        performance.append(_calculate_performance(idx, simulation_results, identification,
+                                                  objective_function))
     return performance
 
 
 def target_factory(data, invert_error=False):
     """factory for creating a target callback"""
-    def _vectorized_wrapper(params, solution_indexes):
+    def _vectorized_wrapper(ga_instance, params, solution_indexes):
+        solutions = _vectorized_target(params.T, data)
+        # OutComment this when using adaptive mutation
+        # if solution_indexes is None:
+        #    return None
+        if invert_error:
+            return [1 / solution for solution in solutions]
+        else:
+            return solutions
+    return _vectorized_wrapper
+
+
+def target_factory_nsga2(data, invert_error=False):
+    """factory for creating a target callback"""
+    def _vectorized_wrapper(ga_instance, params, solution_indexes):
         solutions = _vectorized_target(params.T, data)
         if solution_indexes is None:
             return None
@@ -239,3 +297,14 @@ def _vectorized_wrapper(params, solution_indexes):
         else:
             return solutions
     return _vectorized_wrapper
+
+
+class target_factory_mo_de(Problem):
+
+    def __init__(self, data, **kwargs):
+        super().__init__(**kwargs)
+        self.data = data
+
+    def _evaluate(self, X, out, *args, **kwargs):
+        # store the function values and return them.
+        out["F"] = np.array(_vectorized_target(X.T, self.data))

carmodel_calibration/sumo/sumo_project.py

@@ -15,9 +15,12 @@
 configure_logging()
 _LOGGER = logging.getLogger(__name__)
 
+
 def _get_edge_number(edge):
     return int(float(re.findall(r"([\+\-]*\d*\.*\d+)",
-                               str(edge.attrib["id"]))[0]))
+                                str(edge.attrib["id"]))[0]))
+
+
 class SumoProject:
     """class used to build sumo project for calibration purposes"""
 
@@ -30,7 +33,7 @@ def create_sumo(project_path: Path, model: str, number_network_lanes: int, times
         """
         default_params = ModelParameters.get_defaults_dict()
         cfmodel = ModelParameters.create_parameter_set(".tmp/params", model,
-                                                        **default_params)
+                                                       **default_params)
         if not project_path.exists():
             project_path.mkdir(parents=True)
         network_path = project_path / "calibration_network.xml"
@@ -42,7 +45,8 @@ def create_sumo(project_path: Path, model: str, number_network_lanes: int, times
                                    network_path)
         SumoProject.create_routes_from_network(network_path, routes_path,
                                                number_network_lanes, timestep)
-        SumoProject.create_config(network_path, routes_path, config_path, timestep)
+        SumoProject.create_config(
+            network_path, routes_path, config_path, timestep)
         SumoProject.write_followers_leader(routes_path,
                                            [cfmodel] * number_network_lanes)
 
@@ -63,8 +67,13 @@ def create_network(x_number, y_number, file_path):
         try:
             outs, errs = proc.communicate(timeout=10)
             _LOGGER.debug("%s \noutput:%s",
-                         " ".join(cmd),
-                         str(outs))
+                          " ".join(cmd),
+                          str(outs))
+            if proc.returncode:
+                _LOGGER.error("Network creation failed with %s and error %s",
+                              outs,
+                              errs)
+                raise RuntimeError("Network creation failed")
         except TimeoutExpired as exc:
             proc.kill()
             outs, errs = proc.communicate()
@@ -91,8 +100,8 @@ def create_routes_from_network(network_path: Path, out_path: Path,
         for edge in net_root.findall('edge'):
             if re.match(pattern, str(edge.attrib['id'])):
                 num_edges += 1
-        edges =  net_root.findall('edge')
-        edges = sorted(edges, key= lambda x: _get_edge_number(edge))
+        edges = net_root.findall('edge')
+        edges = sorted(edges, key=lambda x: _get_edge_number(edge))
         for edge in edges:
             if re.match(pattern, str(edge.attrib['id'])):
                 source = str(edge.attrib['id'])
@@ -119,7 +128,7 @@ def create_routes_from_network(network_path: Path, out_path: Path,
                             # <route id="r_TOP" edges="B249A249 A249A248 A248A249\
                             # A249B249"/>
                             edges = [f"{source}", f"A{num_edges-1}B{num_edges-2}",
-                                    dest]
+                                     dest]
                         else:
                             # <route id="r_MIDDLE" edges="B145A145 A145B145"/>
                             edges = [source, dest]
@@ -146,9 +155,9 @@ def create_config(network_file: Path, routes_file: Path, out_path: Path, timeste
                 file_type = "routes"
                 file = str(routes_file)
             _LOGGER.error("%s-file not found at %s"
-                            "\ncreate network-file first",
-                            file_type,
-                            str(file))
+                          "\ncreate network-file first",
+                          file_type,
+                          str(file))
             raise FileNotFoundError
         cmd = [
             "sumo",
@@ -171,8 +180,8 @@ def create_config(network_file: Path, routes_file: Path, out_path: Path, timeste
         try:
             outs, errs = proc.communicate(timeout=10)
             _LOGGER.debug("%s \noutput:%s",
-                         " ".join(cmd),
-                         str(outs))
+                          " ".join(cmd),
+                          str(outs))
         except TimeoutExpired as exc:
             proc.kill()
             outs, errs = proc.communicate()
@@ -270,7 +279,7 @@ def get_network_info(network_file):
         for edge in network_root.findall("edge"):
             if re.match(pattern, str(edge.attrib['id'])):
                 edge_numbers = (re.findall(r"([\+\-]*\d*\.*\d+)",
-                               str(edge.attrib["id"])))
+                                           str(edge.attrib["id"])))
                 edge_numbers = [int(float(item)) for item in edge_numbers]
                 if edge_numbers[0] != edge_numbers[1]:
                     print(str(edge.attrib['id']))

requirements.txt

@@ -3,10 +3,12 @@ scipy
 numpy==1.20.0
 matplotlib
 lxml
-pygad==2.19.0
+pygad==3.3.1
 regex
 argparse
 entrypoints>=0.2.2
 SALib==1.4.5
 tqdm
-cloudpickle==2.2.1
+cloudpickle==2.2.1
+pymoo
+pymoode