18 rework metrics

bo_experiments.py

@@ -0,0 +1,63 @@
+import numpy as np
+import pickle
+import matplotlib.pyplot as plt
+import pandas as pd
+
+from GPy.kern import Matern52, RBF, RatQuad
+
+from emulation.emulator import Emulator
+from emulation.metrics import CompletedJourneysMetric, WaitTimeMetric
+from emulation.utils import run_simulation
+from plot import plot_metric_results
+from simulation_builder.scenarios import single_intersec_lop_2, single_intersec_bal_2, double_intersec_lop_2, \
+    double_intersec_bal_2, cambridge_scenario
+
+
+if __name__ == "__main__":
+    np.set_printoptions(formatter={'float': lambda x: "{0:0.3f}".format(x)})
+
+    interval = (10, 30)
+    N = 300
+    num_init_points = 50
+
+    # metric = WaitTimeMetric
+    # metric_name = "WT"
+    metric = CompletedJourneysMetric
+    metric_name = "CJ"
+
+    plt.style.use('ggplot')
+    plt.rc('font', family='serif')
+    file_id = f'BO_cambridge_lop_{metric_name}_{num_init_points}'
+
+    print(file_id)
+    g, strategy = cambridge_scenario()
+
+    kernel_func = RatQuad
+    kernel_kwargs = {'variance': 1, 'lengthscale': 2}
+
+    e = Emulator(g, strategy)
+
+    results, bo_model = e.bayes_opt(
+        kernel_func,
+        kernel_kwargs,
+        metric,
+        interval=interval,
+        max_iterations=N,
+        progress_N=300,
+        num_init_points=num_init_points)
+
+    with open(f'bo_models/{file_id}.obj', 'wb') as f:
+        pickle.dump(bo_model, f)
+
+    results.to_csv(f'csv_files/{file_id}.csv')
+    data = results[results.columns[-1]]
+
+    plt.plot(np.arange(len(data) - num_init_points) + num_init_points, data[num_init_points:], zorder=1)
+    # plt.scatter(num_init_points, data[num_init_points], s=100, zorder=2)
+    print(f"{metric_name} Results: ")
+    print(data[num_init_points:])
+    # plt.legend(fontsize="large")
+    plt.xlabel("Iterations")
+    plt.ylabel("Completed Journeys" if metric_name == "CJ" else "Average Wait Time")
+    plt.savefig(f"plots/kernel_tuning/cambridge_{metric_name}_", bbox_inches='tight', pad_inches=0.2)
+    plt.close()

emulation/emulator.py

@@ -1,22 +1,84 @@
 from typing import Tuple, Optional
 
 import numpy as np
+import pandas
 import scipy
-from emukit.core import ContinuousParameter
-from emukit.examples.gp_bayesian_optimization.single_objective_bayesian_optimization import GPBayesianOptimization
+
+from GPy.kern import Matern52, RBF
+from GPy.models import GPRegression
+
+from emukit.core import ContinuousParameter, ParameterSpace
+from emukit.core.loop.stopping_conditions import StoppingCondition
+from emukit.core.loop.user_function import UserFunctionWrapper, UserFunctionResult
+from emukit.core.loop.loop_state import LoopState
+from emukit.core.initial_designs import RandomDesign
+from emukit.bayesian_optimization.loops import BayesianOptimizationLoop
+from emukit.bayesian_optimization.acquisitions import ExpectedImprovement
+from emukit.sensitivity.monte_carlo import MonteCarloSensitivity
+from emukit.model_wrappers.gpy_model_wrappers import GPyModelWrapper
 
 from emulation.simulator import Simulator
 from emulation.utils import results_to_df
 
-from simulation_builder.flows import FlowStrategy
+from simulation_builder.flows import FlowStrategy, UniformFlowStrategy
 from simulation_builder.graph import Graph
 
 
+class ProgressStoppingCondition(StoppingCondition):
+    """
+    Stops after N iterations without improvement
+    """
+    def __init__(self, N: int, max_iterations: int) -> None:
+
+        self.N = N  
+        self.max_iterations = max_iterations
+
+        self.best = None
+        self.count = 0
+
+    def should_stop(self, loop_state: LoopState) -> bool:
+
+        # first iteration find best of start points
+        if self.best is None:
+            self.best = np.min(loop_state.Y)
+            print(f'start points: {loop_state.Y.flatten()}')
+            return False
+
+        # get current output
+        current_y = loop_state.Y[-1][0]
+
+        current_x = loop_state.X[-1]
+
+        # new best
+        if current_y < self.best:
+            print(f'iteration {loop_state.iteration}: {current_x} {current_y} - new best!')
+            self.count = 0
+            self.best = current_y
+
+        # not new best
+        else:
+            print(f'iteration {loop_state.iteration}: {current_x} {current_y}')
+            self.count += 1
+
+        # if reached max_iterations return True regardless
+        if loop_state.iteration > self.max_iterations:
+            print('exceeded max iterations, stopping')
+            return True
+
+        # exceeded N
+        elif self.count > self.N:
+            print(f'stopping due to {self.N} iterations without progress')
+            return True
+
+        else:
+            return False
+
+
 class Emulator:
     def __init__(self, graph: Graph, flow_strategy: FlowStrategy, simulation_iterations: int = 1000,
                  fixed_time_period: Optional[float] = None):
         self._g = graph
-        self._strategy = FlowStrategy() if flow_strategy is None else flow_strategy
+        self._strategy = UniformFlowStrategy() if flow_strategy is None else flow_strategy
         self._sim_iterations = simulation_iterations
         self._time_period = fixed_time_period
 
@@ -27,27 +89,106 @@ def __init__(self, graph: Graph, flow_strategy: FlowStrategy, simulation_iterati
         else:
             self._num_params = intersections * 3
 
-    def bayes_opt(self, metric, interval: Tuple[float, float], iterations: int):
-        np.random.seed(42)
+    def bayes_opt(
+            self, 
+            kernel_func,
+            kernel_kwargs,
+            metric, 
+            interval: Tuple[float, float],
+            max_iterations: int,
+            progress_N: int,
+            num_init_points: Optional[int] = 1):
+
+        print(f'\nbayesian optimisation on {metric().name}, interval {interval} with {num_init_points} init points')
 
         sim = Simulator(self._g, metric, self._strategy, self._time_period, self._sim_iterations)
-        x_init = np.random.uniform(*interval, size=(1, self._num_params))
-        y_init = sim.evaluate(x_init)
 
-        parameter_list = [ContinuousParameter(f"p{i}", *interval) for i in range(self._num_params)]
+        target_function = UserFunctionWrapper(sim.evaluate, extra_output_names=['raw metric'])
+
+        # parameter space
+        parameter_space = ParameterSpace([ContinuousParameter(f"x{i}", *interval) for i in range(self._num_params)])
 
-        bo_loop = GPBayesianOptimization(variables_list=parameter_list, X=x_init, Y=y_init, noiseless=True)
-        bo_loop.run_optimization(sim.evaluate, iterations)
+        # random sample init points
+        design = RandomDesign(parameter_space)
+        x_init = design.get_samples(num_init_points)
 
-        return results_to_df(bo_loop.model.X, bo_loop.model.Y, metric().name, self._time_period)
+        # evaluate at all init points, and prep for input to BOLoop
+        y_init = []
+        raw_metric = []
+        for i in range(x_init.shape[0]):
+            output = target_function(x_init[i:i+1])[0]
+            y_init.append(np.expand_dims(output.Y, axis=1))
+            raw_metric.append(output.extra_outputs['raw metric'])
+        raw_metric = np.stack(raw_metric)
+        y_init = np.concatenate(y_init, axis=0)
+
+        # init kernel
+        kernel = kernel_func(self._num_params, **kernel_kwargs)
+
+        # evaluate GP on initial points
+        gpmodel = GPRegression(x_init, y_init, kernel)
+        gpmodel.optimize()
+
+        # no noise in target_function
+        gpmodel.Gaussian_noise.constrain_fixed(0.001)
+
+        # wrap for emukit
+        model = GPyModelWrapper(gpmodel)
+
+        # create the BO loop
+        bo_loop = BayesianOptimizationLoop(
+                space=parameter_space,
+                model=model,
+                acquisition=ExpectedImprovement(model),
+                )
+
+        # put the inital raw metrics into the bo_loop results
+        for i, row in enumerate(bo_loop.loop_state.results):
+            row.extra_outputs['raw metric'] = raw_metric[i]
+
+        stopping_condition = ProgressStoppingCondition(N=progress_N, max_iterations=max_iterations)
+
+        # run optimisation
+        bo_loop.run_loop(target_function, stopping_condition)
+
+        # get x and raw metric values from loop state results
+        x = [step.X for step in bo_loop.loop_state.results]
+        raw_metric = [step.extra_outputs['raw metric'] for step in bo_loop.loop_state.results]
+
+        # convert into arrays
+        x = np.stack(x, axis=0)
+        raw_metric = np.concatenate(raw_metric)
+
+        return results_to_df(x, self._time_period, raw_metric, metric().name, num_init_points), bo_loop.model
+
+    def sensitivity(self, bo_model, interval: Tuple[float, float], num_mc: int = 50000):
+
+        parameter_list = [ContinuousParameter(f"x{i}", *interval) for i in range(self._num_params)]
+
+        senstivity = MonteCarloSensitivity(model=bo_model, input_domain=ParameterSpace(parameter_list))
+        main_effects, total_effects, _ = senstivity.compute_effects(num_monte_carlo_points=num_mc)
+
+        # converting from dict into arrays for results_df function
+        main_effects = np.fromiter(main_effects.values(), dtype=float)
+        main_effects = np.reshape(main_effects, (1, len(main_effects)))
+
+        total_effects = np.fromiter(total_effects.values(), dtype=float)
+        total_effects = np.reshape(total_effects, (1, len(total_effects)))
+
+        return results_to_df(main_effects, self._time_period), \
+            results_to_df(total_effects, self._time_period)
 
     def grid_search_opt(self, metric, interval: Tuple[float, float], steps_per_axis: int):
         """Evaluates target_function on all combinations of parameters taken from the same interval"""
-        np.random.seed(42)
+
+        print(f'\ngrid search on {metric().name}, interval {interval} with {steps_per_axis**self._num_params} grid points')
 
         sim = Simulator(self._g, metric, self._strategy, self._time_period, self._sim_iterations)
 
-        x_min, f_min, grid, results = scipy.optimize.brute(func=sim.evaluate,
+        # lambda function for selecting raw metric from outputs
+        target_function = lambda x: sim.evaluate(x)[1]
+
+        x_min, f_min, grid, results = scipy.optimize.brute(func=target_function,
                                                            ranges=(interval,) * self._num_params,
                                                            Ns=steps_per_axis,
                                                            full_output=True,
@@ -56,4 +197,4 @@ def grid_search_opt(self, metric, interval: Tuple[float, float], steps_per_axis:
         results = results.flatten()
         grid = np.moveaxis(grid, 0, self._num_params).reshape(-1, self._num_params)
 
-        return results_to_df(grid, results, metric().name, self._time_period)
+        return results_to_df(grid, self._time_period, results, metric().name)

emulation/metrics.py

@@ -0,0 +1,96 @@
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from numbers import Number
+from typing import List
+
+
+@dataclass
+class Report:
+    target_metric: Number  # 'normalised' for minimisation
+    raw_metric: Number
+
+    def __iter__(self):
+        return iter((self.target_metric, self.raw_metric))
+
+
+class Metric(ABC):
+    @abstractmethod
+    def update(self, eng):
+        pass
+
+    @abstractmethod
+    def report(self) -> Report:
+        pass
+
+    name: str
+
+
+class CompletedJourneysMetric(Metric):
+    def __init__(self):
+        self._total_vehicles = set() # unique vehicles in simulation
+        self.name = 'completed journeys'
+
+    def update(self, eng):
+
+        # get all vehicles currently in simulation
+        self._current_vehicles = set(eng.get_vehicles(include_waiting=True))
+
+        # add new vehicles to total count
+        self._total_vehicles |= self._current_vehicles
+
+    def report(self) -> Report:
+
+        # find the vehicles that had left simulator by last step
+        total_completed = len(self._total_vehicles - self._current_vehicles)
+
+        # optimisation target is 1 - (normalised total completed)
+        target_metric = 1 - total_completed / len(self._total_vehicles)
+
+        return Report(target_metric, total_completed)
+
+
+class WaitTimeMetric(Metric):
+    def __init__(self):
+        self._unique_vehicles = set()  # unique vehicles in simulation (should be constant for given graph + flow)
+        self._waiting_vehicle_steps = []  # a vehicle step is 1 vehicle waiting for 1 step
+        self.name = 'average steps waiting'
+
+    def update(self, eng):
+
+        # get all vehicles currently in simulation, including those waiting offscreen
+        current_vehicles = eng.get_vehicles(include_waiting=True)
+
+        # add new vehicles to set of unique 
+        self._unique_vehicles |= set(current_vehicles)
+
+        # find the waiting vehicles
+        waiting = 0
+        for v in current_vehicles:
+
+            # get vehicle info from engine
+            info = eng.get_vehicle_info(v)
+
+            # check if car on road
+            if info['running'] == '1':
+                # get speed of car
+                speed = float(info['speed'])
+            else:
+                # cars not on road have no key 'speed'
+                speed = 0.0
+
+            # check if speed is below threshold, add to count of waiting
+            if speed < 0.1:
+                waiting += 1
+
+        # add waiting vehicles to stepwise list
+        self._waiting_vehicle_steps.append(waiting)
+
+    def report(self) -> Report:
+
+        # computes the average proportion of time vehicles spent waiting in their journey
+        total_average = sum(self._waiting_vehicle_steps) / len(self._unique_vehicles)
+
+        # divide by number of steps (can't normalise as is not bounded)
+        target_metric = total_average / len(self._waiting_vehicle_steps)
+
+        return Report(target_metric, total_average)