diff --git a/bo_experiments.py b/bo_experiments.py new file mode 100644 index 0000000..9462188 --- /dev/null +++ b/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() diff --git a/emulation/emulator.py b/emulation/emulator.py index b8d9afd..c66c3b8 100644 --- a/emulation/emulator.py +++ b/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) diff --git a/emulation/metrics.py b/emulation/metrics.py new file mode 100644 index 0000000..fb32ce1 --- /dev/null +++ b/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) diff --git a/emulation/simulator.py b/emulation/simulator.py index baa8749..e2f27da 100644 --- a/emulation/simulator.py +++ b/emulation/simulator.py @@ -6,7 +6,7 @@ import numpy as np import cityflow as cf -from simulation_builder.flows import FlowStrategy, graph_to_flow +from simulation_builder.flows import FlowStrategy, graph_to_flow, UniformFlowStrategy from simulation_builder.graph import Graph from simulation_builder.roadnets import graph_to_roadnet @@ -26,7 +26,7 @@ def __init__(self, g: Graph, metric, strategy=None, timing_period: Optional[int] self.g = g self.metric = metric - self.strategy = FlowStrategy() if strategy is None else strategy + self.strategy = UniformFlowStrategy() if strategy is None else strategy self.timing_period = timing_period self.steps = steps @@ -36,6 +36,7 @@ def __init__(self, g: Graph, metric, strategy=None, timing_period: Optional[int] with open("cityflow_config/flows/auto_flow.json", 'w') as f: f.write(json.dumps(flow, indent=4)) + def evaluate(self, x): """ Parameters @@ -72,51 +73,51 @@ def evaluate(self, x): eng.next_step() metric.update(eng) - aggregate, _ = metric.report() + target_metric, raw_metric = metric.report() - return np.array([[aggregate]]) + return np.array([[target_metric]]), np.array([[raw_metric]]) - def multithreaded_evaluate(self, x): - """ - Parameters - ---------- - x: A 1D numpy array of traffic light phase timings + # def multithreaded_evaluate(self, x): + # """ + # Parameters + # ---------- + # x: A 1D numpy array of traffic light phase timings - Returns - ------- - The resulting aggregate metric calculated after N simulation iterations, with traffic light timings x. - """ - x = np.array(np.array_split(x.flatten(), len(self.intersections))) - if self.timing_period is not None: - x3 = self.timing_period - x.sum(axis=1) - x = np.insert(x, x.shape[1] - 1, x3, axis=1) + # Returns + # ------- + # The resulting aggregate metric calculated after N simulation iterations, with traffic light timings x. + # """ + # x = np.array(np.array_split(x.flatten(), len(self.intersections))) + # if self.timing_period is not None: + # x3 = self.timing_period - x.sum(axis=1) + # x = np.insert(x, x.shape[1] - 1, x3, axis=1) - traffic_light_phases = {intersection: timing for (intersection, timing) in zip(self.intersections, x.tolist())} - roadnet = graph_to_roadnet(self.g, traffic_light_phases, intersection_width=50, lane_width=8) + # traffic_light_phases = {intersection: timing for (intersection, timing) in zip(self.intersections, x.tolist())} + # roadnet = graph_to_roadnet(self.g, traffic_light_phases, intersection_width=50, lane_width=8) - # Generate UUID tag to ensure each thread's config file is unique - name = str(uuid.uuid1())[:8] + # # Generate UUID tag to ensure each thread's config file is unique + # name = str(uuid.uuid1())[:8] - roadnet_file = f"cityflow_config/roadnets/auto_roadnet_{name}.json" - with open(roadnet_file, 'w') as f: - f.write(json.dumps(roadnet, indent=4)) - config_file = f"cityflow_config/config_{name}.json" - with open(f"cityflow_config/config.json", 'r') as f: - config = json.loads(f.read()) + # roadnet_file = f"cityflow_config/roadnets/auto_roadnet_{name}.json" + # with open(roadnet_file, 'w') as f: + # f.write(json.dumps(roadnet, indent=4)) + # config_file = f"cityflow_config/config_{name}.json" + # with open(f"cityflow_config/config.json", 'r') as f: + # config = json.loads(f.read()) - config["roadnetFile"] = f"roadnets/auto_roadnet_{name}.json" - with open(config_file, 'w') as f: - f.write(json.dumps(config, indent=4)) + # config["roadnetFile"] = f"roadnets/auto_roadnet_{name}.json" + # with open(config_file, 'w') as f: + # f.write(json.dumps(config, indent=4)) - eng = cf.Engine(config_file, thread_num=1) + # eng = cf.Engine(config_file, thread_num=1) - metric = self.metric() + # metric = self.metric() - for _ in range(self.steps): - eng.next_step() - metric.update(eng) + # for _ in range(self.steps): + # eng.next_step() + # metric.update(eng) - aggregate, _ = metric.report() - os.remove(roadnet_file) - os.remove(config_file) - return np.array([[aggregate]]) + # aggregate, _ = metric.report() + # os.remove(roadnet_file) + # os.remove(config_file) + # return np.array([[aggregate]]) diff --git a/emulation/utils.py b/emulation/utils.py index 442272d..edb6391 100644 --- a/emulation/utils.py +++ b/emulation/utils.py @@ -1,10 +1,19 @@ -from typing import Optional +import json +from typing import List, Optional, Dict import numpy as np import pandas +import cityflow as cf -def results_to_df(x, y, metric_name: str, time_period: Optional[float]): +from emulation.metrics import CompletedJourneysMetric +from simulation_builder.flows import FlowStrategy, graph_to_flow +from simulation_builder.graph import Graph +from simulation_builder.roadnets import graph_to_roadnet + + +def results_to_df(x, time_period: Optional[float], y: Optional[List[float]] = None, metric_name: Optional[str] = None, + num_init_points: Optional[int] = None): """ Parameters ---------- @@ -19,22 +28,44 @@ def results_to_df(x, y, metric_name: str, time_period: Optional[float]): """ d = {} + if num_init_points != None: + eval_type = ['init'] * num_init_points + ['BO'] * (x.shape[0] - num_init_points) + d['eval type'] = eval_type + for i in range(x.shape[1]): if time_period is not None: node, phase = i // 3, i % 3 d[f"x_{node}_{phase}"] = np.array(x[:, i]) if i % 3 == 2: - d[f"x_{node}_3"] = time_period - x[:, i - 2:i + 1].sum(axis=1) + d[f"x_{node}_3 (i)"] = time_period - x[:, i - 2:i + 1].sum(axis=1) else: node, phase = i // 4, i % 4 d[f"x_{node}_{phase}"] = np.array(x[:, i]) - d[metric_name] = np.array(y).flatten() + if y is not None: + d[metric_name] = np.array(y).flatten() df = pandas.DataFrame(data=d) return df +def run_simulation(g: Graph, strategy: FlowStrategy, n=1000, traffic_light_phases: Optional[Dict] = None): + roadnet = graph_to_roadnet(g, traffic_light_phases, intersection_width=50, lane_width=8) + flow = graph_to_flow(g, strategy) + with open("cityflow_config/roadnets/auto_roadnet.json", 'w') as f: + f.write(json.dumps(roadnet, indent=4)) + + with open("cityflow_config/flows/auto_flow.json", 'w') as f: + f.write(json.dumps(flow, indent=4)) + eng = cf.Engine("cityflow_config/config.json", thread_num=1) + + metric = CompletedJourneysMetric() + for _ in range(n): + eng.next_step() + metric.update(eng) + print(metric.report()) + + # Test functions def forrester(x): return (6 * x - 2) ** 2 * np.sin(12 * x - 4) diff --git a/evaluate_kernels.py b/evaluate_kernels.py new file mode 100644 index 0000000..f23ca74 --- /dev/null +++ b/evaluate_kernels.py @@ -0,0 +1,114 @@ +import pandas +import os +import re + +def mean_best(files, metric): + + data = {'RBF': [], 'RQ': [], 'M52': []} + + assert not len(files) % len(data) + + files_per_kernel = (len(files) / len(data)) + + for f in files: + df = pandas.read_csv('csv_files/'+f, index_col=0) + + if metric == 'WT': + best = df['average steps waiting'].min() + else: + best = df['completed journeys'].max() + + if 'RBF' in f: + data['RBF'] += [best] + elif 'RQ' in f: + data['RQ'] += [best] + elif 'M52' in f: + data['M52'] += [best] + + return(data) + + +def convergence(files, metric, eps, N): + + data = {'RBF': [], 'RQ': [], 'M52': []} + + assert not len(files) % len(data) + + files_per_kernel = (len(files) / len(data)) + + for f in files: + df = pandas.read_csv('csv_files/'+f, index_col=0) + + converged = False + prev = None + count = 0 + + for i, row in df.iterrows(): + + if metric == 'WT': + current = row['average steps waiting'] + else: + current = row['completed journeys'] + #print('curr', current) + #print('prev', prev) + + if prev is not None: + + diff = abs(current - prev) + + if diff < eps: + count += 1 + #print('diff', diff) + #print('count', count) + #breakpoint() + else: + count = 0 + + if count >= N: + if 'RBF' in f: + data['RBF'] += [i] + elif 'RQ' in f: + data['RQ'] += [i] + elif 'M52' in f: + data['M52'] += [i] + + converged = True + break + + prev = current + + if not converged: + if 'RBF' in f: + data['RBF'] += ['DNC'] + elif 'RQ' in f: + data['RQ'] += ['DNC'] + elif 'M52' in f: + data['M52'] += ['DNC'] + + return(data) + + +csv_files = sorted(os.listdir('csv_files/')) + +DL2_files = [f for f in csv_files if 'DL2' in f] +DB2_files = [f for f in csv_files if 'DB2' in f] + +DL2_CJ_files = [f for f in DL2_files if 'CJ' in f] +DB2_CJ_files = [f for f in DB2_files if 'CJ' in f] + +DL2_WT_files = [f for f in DL2_files if 'WT' in f] +DB2_WT_files = [f for f in DB2_files if 'WT' in f] + +print('DL2 CJ', convergence(DL2_CJ_files, 'CJ', 3, 3)) +print('DB2 CJ', convergence(DB2_CJ_files, 'CJ', 3, 3)) +print('DL2 WT', convergence(DL2_WT_files, 'WT', 0.2, 3)) +print('DB2 WT', convergence(DB2_WT_files, 'WT', 1, 3)) + +print() + +print('DL2 CJ', mean_best(DL2_CJ_files, 'CJ')) +print('DB2 CJ', mean_best(DB2_CJ_files, 'CJ')) +print('DL2 WT', mean_best(DL2_WT_files, 'WT')) +print('DB2 WT', mean_best(DB2_WT_files, 'WT')) +breakpoint() + diff --git a/main.py b/main.py index c585c7c..bf4c26c 100644 --- a/main.py +++ b/main.py @@ -1,32 +1,101 @@ +import matplotlib.pyplot as plt import numpy as np +import pickle + +import pandas as pd +from GPy.kern import Matern52, RBF, RatQuad from emulation.emulator import Emulator -from simulation_builder.flows import CustomEndpointFlowStrategy, FlowStrategy -from simulation_builder.graph import Graph, I_graph +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, cambridge_scenario, \ + double_intersec_bal_2 + + + +def plot_cum(file_name, new_name=""): + df = pd.read_csv(file_name) + + met_name = df.columns[-1] + res = df[met_name] + + file = f'plots/cambridge/cum_{new_name}.png' + + if False: + cuma = res.cummin() + else: + cuma = res.cummax() + + plt.figure() + plt.style.use('ggplot') + plt.rc('font', family='serif') + plt.plot(np.arange(len(res)), res, ) + plt.plot(np.arange(len(res)), cuma) + plt.legend(['Actual', 'Cumulative']) + plt.xlabel("Iterations") + plt.ylabel("Completed Journeys") + # plt.title(file_id) + plt.savefig(file, bbox_inches='tight', pad_inches=0.2) -from metrics.metrics import CompletedJourneysMetric if __name__ == "__main__": np.set_printoptions(formatter={'float': lambda x: "{0:0.3f}".format(x)}) - # g = Graph([(0, -400), (0, 0), (0, 400), (-400, 0), (400, 0)], - # [((0, -400), (0, 0)), ((0, 400), (0, 0)), ((-400, 0), (0, 0)), ((400, 0), (0, 0))]) - # - # - # - # strategy = CustomEndpointFlowStrategy(start_flows={(0, -400): 1, - # (0, 400): 240, - # (-400, 0): 240, - # (400, 0): 240}, - # end_flows={(0, -400): 240, - # (0, 400): 1, - # (-400, 0): 240, - # (400, 0): 240}) + interval = (0.1, 30) - g = I_graph() - strategy = FlowStrategy() + np.random.seed(42) - e = Emulator(g, strategy, fixed_time_period=60) + # g, strategy = double_intersec_bal_2() + g, strategy = cambridge_scenario() + run_simulation(g, strategy) + # plot_cum("csv_files/BO_cambridge_lop_CJ_50.csv", "CJ_lop") - # print(e.bayes_opt(CompletedJourneysMetric, interval=(0.1, 20), iterations=5)) - print(e.grid_search_opt(CompletedJourneysMetric, interval=(0.1, 20), steps_per_axis=2)) + # for scenario in ['DB2']: + # for kernel_name in ['RQ', 'RBF']: + # + # # define id + # scenario = 'SL2' + # num_init_points = 1 + # metric_name = 'WT' + # variance = 2 + # + # # id -> config + # if scenario == 'SL2': + # g, strategy = single_intersec_lop_2() + # if scenario == 'SB2': + # g, strategy = single_intersec_bal_2() + # + # if kernel_name == 'M52': + # kernel_func = Matern52 + # elif kernel_name == 'RBF': + # kernel_func = RBF + # else: + # kernel_func = RatQuad + # + # if metric_name == 'WT': + # metric = WaitTimeMetric + # else: + # metric = CompletedJourneysMetric + # + # kernel_kwargs = {'variance': variance} + # + # e = Emulator(g, strategy) + # + # results, bo_model = e.bayes_opt( + # kernel_func, + # kernel_kwargs, + # metric, + # interval=interval, + # max_iterations=250, + # progress_N=500, + # num_init_points=num_init_points) + # + # file_id = f'BO_{scenario}_{metric_name}_{kernel_name}_{variance}_{num_init_points}' + # + # results.to_csv(f'csv_files/{file_id}.csv') + # + # plot_metric_results(results, f'plots/{file_id}.png') + # + # with open(f'bo_models/{file_id}.obj', 'wb') as f: + # pickle.dump(bo_model, f) diff --git a/metrics/metrics.py b/metrics/metrics.py deleted file mode 100644 index 6870c5f..0000000 --- a/metrics/metrics.py +++ /dev/null @@ -1,65 +0,0 @@ -from abc import ABC, abstractmethod -from dataclasses import dataclass -from numbers import Number -from typing import List - - -@dataclass -class Report: - aggregate: Number - data: List[Number] - - def __iter__(self): - return iter((self.aggregate, self.data)) - - -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() - self._prev_step = set() - self._completed_journeys = [0] - self.name = 'completed journeys' - - def update(self, eng): - curr_step = set(eng.get_vehicles(include_waiting=True)) - self._total_vehicles |= curr_step - self._completed_journeys.append(len(self._prev_step - curr_step) + self._completed_journeys[-1]) - self._prev_step = curr_step - - def report(self) -> Report: - return Report(1 - self._completed_journeys[-1]/len(self._total_vehicles), self._completed_journeys) - - -class WaitTimeMetric(Metric): - """ - Reports the overall average waiting time, and the proportion of cars waiting at each time step. - """ - - def __init__(self): - self._waiting_vehicles = [] - self._total_vehicles = [] - self.name = 'wait time' - - def update(self, eng): - vehicles = eng.get_vehicles(include_waiting=True) - wait_time = sum([float(eng.get_vehicle_info(v)['speed']) < 0.1 for v in vehicles]) - self._total_vehicles.append(len(vehicles)) - self._waiting_vehicles.append(wait_time) - - def report(self) -> Report: - total_average = sum(self._waiting_vehicles) / sum(self._total_vehicles) - proportion_waiting = [wait / total if total > 0 else 0 for wait, total in - zip(self._waiting_vehicles, self._total_vehicles)] - return Report(total_average, proportion_waiting) diff --git a/metrics/useful_metrics.py b/metrics/useful_metrics.py deleted file mode 100644 index 58d2c3c..0000000 --- a/metrics/useful_metrics.py +++ /dev/null @@ -1,73 +0,0 @@ -import numpy as np - -set_ids = set() -last_step = None -counter = 0 - -def completed_journeys(eng, online = False): - """Returns online or total completed journeys""" - global last_step - # List of vehicle IDs currently on road - current_ids = eng.get_vehicles(include_waiting=False) - - # If want completed journeys at each time step: - if online: - # If function is called for the first time - if last_step is None: - last_step = current_ids - return - # Yields the elements in last_step that are NOT in current_ids - completed_ids = np.setdiff1d(last_step, current_ids) - completed_journeys = len(completed_ids) - - last_step = current_ids - return completed_journeys - - # total = set of car ids - number of cars on road at last step - set_ids.update(current_ids) - total_completed_journeys = len(set_ids) - len(eng.get_vehicles()) - - return total_completed_journeys - -def count_waiting_cars(eng, current_ids, counter): - for car in current_ids: - # CityFlow defines a waiting car - # as speed < 0.1 mph - if float(eng.get_vehicle_info(car)['speed']) < 0.1: - counter += 1 - - -def wait_time(eng, _ = False, steps = False, online = False): - """Returns online or total average wait""" - #List of IDs of cars currently on road - current_ids = eng.get_vehicles(include_waiting=False) - - # If want average wait across entire simulation - global counter - for car in current_ids: - # CityFlow defines a waiting car - # as speed < 0.1 mph - if float(eng.get_vehicle_info(car)['speed']) < 0.1: - counter += 1 - set_ids.update(current_ids) - # If last step of simulation, calculate total avg - if _ == steps - 1: - # avg = num of times cars were idle in sim / num cars in sim - total_avg = counter / len(set_ids) - return total_avg - - # If want average wait at each time step - if online: - counter = 0 - for car in current_ids: - # CityFlow defines a waiting car - # as speed < 0.1 mph - if float(eng.get_vehicle_info(car)['speed']) < 0.1: - counter += 1 - if counter == 0: - return - # avg = num cars that were idle / num cars on road - avg = counter / len(eng.get_vehicles()) - return avg - - diff --git a/plot.py b/plot.py new file mode 100644 index 0000000..be9d8f4 --- /dev/null +++ b/plot.py @@ -0,0 +1,48 @@ +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +def plot_metric_results(df, file_id, minimisation=True): + metric_name = df.columns[-1] + results = df[metric_name] + + file_name = f'plots/{file_id}.png' + + if minimisation: + cuma = results.cummin() + else: + cuma = results.cummax() + + plt.figure() + plt.style.use('ggplot') + plt.rc('font', family='serif') + plt.plot(np.arange(len(results)), results,) + # plt.plot(np.arange(len(results)), cuma) + plt.legend(['actual', 'cumulative']) + # plt.title(file_id) + plt.savefig(file_name) + + return plt + +def plot_sensitivity(main_effects, total_effects, file_name): + + params = main_effects.columns + main_effects = main_effects.values.tolist()[0] + total_effects = total_effects.values.tolist()[0] + + plt.figure() + plt.style.use('ggplot') + plt.rc('font', family='serif') + + X_axis = np.arange(len(params)) + + plt.bar(X_axis - 0.2, main_effects, 0.4, label='Main Effects') + plt.bar(X_axis + 0.2, [-x/2 if x < 0 else x for x in total_effects], 0.4, label='Total Effects') + + # plt.xticks(X_axis, params) + plt.ylabel('Sensitivity') + + plt.savefig(file_name, bbox_inches='tight', pad_inches=0.2) + + return plt diff --git a/plots/.gitkeep b/plots/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/sensitivity_experiments.py b/sensitivity_experiments.py new file mode 100644 index 0000000..66daded --- /dev/null +++ b/sensitivity_experiments.py @@ -0,0 +1,38 @@ +import numpy as np +import pickle + +from emulation.emulator import Emulator +from simulation_builder.scenarios import double_intersec_bal_2, double_intersec_lop_2, cambridge_scenario + +from plot import plot_sensitivity + +if __name__ == "__main__": + np.set_printoptions(formatter={'float': lambda x: "{0:0.3f}".format(x)}) + + np.random.seed(42) + + interval = (10, 30) + + for kernel_name in ['RQ']: + g, strategy = cambridge_scenario() + + e = Emulator(g, strategy) + + # define id + lengthscale = 2 + num_init_points = 50 + metric_name = 'CJ' + variance = 2 + + file_id = f'BO_cambridge_lop_CJ_50' + + with open(f'bo_models/{file_id}.obj', 'rb') as f: + bo_model = pickle.load(f) + + main_effects, total_effects = e.sensitivity(bo_model, interval=interval) + print(main_effects) + print(total_effects) + + plt = plot_sensitivity(main_effects, total_effects, f'sensitivity_plots/{file_id}.png') + + diff --git a/simulation_builder/flows.py b/simulation_builder/flows.py index 0df2fc2..a9554f3 100644 --- a/simulation_builder/flows.py +++ b/simulation_builder/flows.py @@ -1,4 +1,5 @@ import math +from abc import ABC, abstractmethod from logging import warning from numpy import random @@ -34,10 +35,21 @@ def json(self) -> Dict: } -class FlowStrategy: - # Default flow strategy creates one flow per route. +class FlowStrategy(ABC): + @abstractmethod def gen_flows(self, route: List[Tuple[int, int]]) -> List[Flow]: - return [Flow(route)] + pass + + +class UniformFlowStrategy(FlowStrategy): + """ + Creates exactly one flow per route, all initialised with the same interval + """ + def __init__(self, interval=5.0): + self._interval = interval + + def gen_flows(self, route: List[Tuple[int, int]]) -> List[Flow]: + return [Flow(route, interval=self._interval)] class RandomFlowStrategy(FlowStrategy): @@ -82,7 +94,36 @@ def gen_flows(self, route: List[Tuple[int, int]]) -> List[Flow]: return [Flow(route, interval=self._default)] -def graph_to_flow(g: Graph, strategy: FlowStrategy = FlowStrategy()) -> List[Dict]: +class ManualFlowStrategy(FlowStrategy): + """ + Takes a dictionary mapping routes defined by start/end pairs to flow intervals. Assumes flows are uniquely defined + by their start and end points. + """ + def __init__(self, flows: Dict[Tuple[Tuple, Tuple], float]): + self._flows = flows + + def gen_flows(self, route: List[Tuple[int, int]]) -> List[Flow]: + start, end = route[0], route[-1] + if (start, end) in self._flows: + return [Flow(route, interval=self._flows[(start, end)])] + return [] + + +class CompositeFlowStrategy(FlowStrategy): + """ + Takes a list of flow strategies and combines them into one. + """ + def __init__(self, strategies: List[FlowStrategy]): + self._strategies = strategies + + def gen_flows(self, route: List[Tuple[int, int]]) -> List[Flow]: + flows = [] + for strategy in self._strategies: + flows += strategy.gen_flows(route) + return flows + + +def graph_to_flow(g: Graph, strategy: FlowStrategy = UniformFlowStrategy()) -> List[Dict]: paths = all_pairs_shortest_paths(g) flows = [] for start in paths: diff --git a/simulation_builder/roadnets.py b/simulation_builder/roadnets.py index cfdff0b..10e13d7 100644 --- a/simulation_builder/roadnets.py +++ b/simulation_builder/roadnets.py @@ -157,8 +157,9 @@ def gen_intersections(g: Graph, traffic_light_phases: Optional[Dict], intersecti processed_light_phases = [{"time": 30, "availableRoadLinks": list(phase)} for phase in light_phases if len(phase) > 0] else: + timings = traffic_light_phases[u] if u in traffic_light_phases else [10, 10, 10, 10] processed_light_phases = [{"time": t, "availableRoadLinks": list(phase)} - for phase, t in zip(light_phases, traffic_light_phases[u]) + for phase, t in zip(light_phases, timings) if len(phase) > 0] intersection["trafficLight"] = { diff --git a/simulation_builder/scenarios.py b/simulation_builder/scenarios.py new file mode 100644 index 0000000..343526a --- /dev/null +++ b/simulation_builder/scenarios.py @@ -0,0 +1,110 @@ +from typing import Tuple + +from simulation_builder.graph import Graph +from simulation_builder.flows import CustomEndpointFlowStrategy, FlowStrategy, ManualFlowStrategy, UniformFlowStrategy, \ + CompositeFlowStrategy + + +def single_intersec_lop_2() -> Tuple[Graph, FlowStrategy]: + strategy = ManualFlowStrategy({ + ((0, -400), (0, 400)): 4, + ((-400, 0), (400, 0)): 100, + ((400, 0), (-400, 0)): 100 + }) + + return single_intersec_g(), strategy + + +def single_intersec_bal_2() -> Tuple[Graph, FlowStrategy]: + strategy = UniformFlowStrategy(interval=10) + return single_intersec_g(), strategy + + +def single_intersec_g() -> Graph: + vertices = [ + (0, -400), + (0, 0), + (0, 400), + (-400, 0), + (400, 0), + ] + + edges = [ + ((0, -400), (0, 0)), + ((0, 400), (0, 0)), + ((-400, 0), (0, 0)), + ((400, 0), (0, 0)), + ] + + return Graph(vertices=vertices, edges=edges) + + +def double_intersec_bal_2() -> Tuple[Graph, FlowStrategy]: + strategy = UniformFlowStrategy(interval=20) + return double_intersec_g(), strategy + + +def double_intersec_lop_2() -> Tuple[Graph, FlowStrategy]: + strategy = ManualFlowStrategy({ + ((-400, 0), (800, 0)): 4, + ((800, 0), (-400, 0)): 4, + ((0, -400), (0, 400)): 100, + ((400, 400), (400, -400)): 100, + }) + + return double_intersec_g(), strategy + + +def double_intersec_g() -> Graph: + vertices = [ + (-400, 0), + (0, 0), + (400, 0), + (800, 0), + (0, -400), + (400, -400), + (0, 400), + (400, 400), + ] + + edges = [ + ((-400, 0), (0, 0)), + ((0, 0), (400, 0)), + ((400, 0), (800, 0)), + ((0, 0), (0, 400)), + ((400, 0), (400, 400)), + ((0, 0), (0, -400)), + ((400, 0), (400, -400)) + ] + + return Graph(vertices=vertices, edges=edges) + + +def cambridge_scenario() -> Tuple[Graph, FlowStrategy]: + vertices = [(0, 0), (400, 0), (0, 200), (400, 200), (800, 200), (1000, 200), (0, 900), (400, 900), + (800, 900), (1000, 900), (400, 1300)] + + edges = [ + ((0, 200), (400, 200)), + ((400, 200), (800, 200)), + ((800, 200), (1000, 200)), + ((0, 900), (400, 900)), + ((400, 900), (800, 900)), + ((800, 900), (1000, 900)), + + ((0, 0), (0, 200)), + ((400, 0), (400, 200)), + ((0, 200), (0, 900)), + ((400, 200), (400, 900)), + ((800, 200), (800, 900)), + ((400, 900), (400, 1300)), + ] + + strategy = CompositeFlowStrategy([ + UniformFlowStrategy(interval=20), + # ManualFlowStrategy({ + # ((400, 0), (400, 1300)): 20 + # }) + ]) + + return Graph(vertices, edges), strategy diff --git a/simulation_builder/utils.py b/simulation_builder/utils.py deleted file mode 100644 index 8c3bb80..0000000 --- a/simulation_builder/utils.py +++ /dev/null @@ -1,122 +0,0 @@ -import json -import os - -from CityFlow.tools.generator.generate_json_from_grid import gridToRoadnet - - -def generate_roadnet(rowNum: int, colNum: int, rowDistance: int = 300, columnDistance: int = 300, - intersectionWidth: int = 100, numLeftLanes: int = 1, numStraightLanes: int = 1, - numRightLanes: int = 1, laneMaxSpeed: float = 16.67, vehLen: float = 5.0, vehWidth: float = 2.0, - vehMaxPosAcc: float = 2.0, vehMaxNegAcc: float = 4.5, vehUsualPosAcc: float = 2.0, - vehUsualNegAcc: float = 4.5, - vehMinGap: float = 2.5, vehMaxSpeed: float = 16.67, vehHeadwayTime: float = 1.5, - directory: str = "cityflow_config", roadnetFile: str = None, turn: bool = False, - tlPlan: bool = False, interval: float = 2.0, - flowFile: str = None): - """ - Generates grid-shaped roadnet and flow JSON files, and updates config.json file. - - (Code adapted from CityFlows/tools/generator directory) - """ - - # Helper function - def generate_route(rowNum, colNum, turn=False): - routes = [] - move = [(1, 0), (0, 1), (-1, 0), (0, -1)] - - def get_straight_route(start, direction, step): - x, y = start - route = [] - for _ in range(step): - route.append("road_%d_%d_%d" % (x, y, direction)) - x += move[direction][0] - y += move[direction][1] - return route - - for i in range(1, rowNum + 1): - routes.append(get_straight_route((0, i), 0, colNum + 1)) - routes.append(get_straight_route((colNum + 1, i), 2, colNum + 1)) - for i in range(1, colNum + 1): - routes.append(get_straight_route((i, 0), 1, rowNum + 1)) - routes.append(get_straight_route((i, rowNum + 1), 3, rowNum + 1)) - - if turn: - def get_turn_route(start, direction): - if direction[0] % 2 == 0: - step = min(rowNum * 2, colNum * 2 + 1) - else: - step = min(colNum * 2, rowNum * 2 + 1) - x, y = start - route = [] - cur = 0 - for _ in range(step): - route.append("road_%d_%d_%d" % (x, y, direction[cur])) - x += move[direction[cur]][0] - y += move[direction[cur]][1] - cur = 1 - cur - return route - - routes.append(get_turn_route((1, 0), (1, 0))) - routes.append(get_turn_route((0, 1), (0, 1))) - routes.append(get_turn_route((colNum + 1, rowNum), (2, 3))) - routes.append(get_turn_route((colNum, rowNum + 1), (3, 2))) - routes.append(get_turn_route((0, rowNum), (0, 3))) - routes.append(get_turn_route((1, rowNum + 1), (3, 0))) - routes.append(get_turn_route((colNum + 1, 1), (2, 1))) - routes.append(get_turn_route((colNum, 0), (1, 2))) - - return routes - - if roadnetFile is None: - roadnetFile = "roadnet_%d_%d%s.json" % (rowNum, colNum, "_turn" if turn else "") - if flowFile is None: - flowFile = "flow_%d_%d%s.json" % (rowNum, colNum, "_turn" if turn else "") - - grid = { - "rowNumber": rowNum, - "columnNumber": colNum, - "rowDistances": [rowDistance] * (colNum - 1), - "columnDistances": [columnDistance] * (rowNum - 1), - "outRowDistance": rowDistance, - "outColumnDistance": columnDistance, - "intersectionWidths": [[intersectionWidth] * colNum] * rowNum, - "numLeftLanes": numLeftLanes, - "numStraightLanes": numStraightLanes, - "numRightLanes": numRightLanes, - "laneMaxSpeed": laneMaxSpeed, - "tlPlan": tlPlan - } - - json.dump(gridToRoadnet(**grid), open(os.path.join(directory, "roadnets/", roadnetFile), "w"), indent=4) - - vehicle_template = { - "length": vehLen, - "width": vehWidth, - "maxPosAcc": vehMaxPosAcc, - "maxNegAcc": vehMaxNegAcc, - "usualPosAcc": vehUsualPosAcc, - "usualNegAcc": vehUsualNegAcc, - "minGap": vehMinGap, - "maxSpeed": vehMaxSpeed, - "headwayTime": vehHeadwayTime - } - routes = generate_route(rowNum, colNum, turn) - flow = [] - for route in routes: - flow.append({ - "vehicle": vehicle_template, - "route": route, - "interval": interval, - "startTime": 0, - "endTime": -1 - }) - json.dump(flow, open(os.path.join(directory, "flows/", flowFile), "w"), indent=4) - - with open("cityflow_config/config.json", "r") as f: - config_file = json.load(f) - - config_file["roadnetFile"] = f"roadnets/{roadnetFile}" - config_file["flowFile"] = f"flows/{flowFile}" - - with open("cityflow_config/config.json", "w") as f: - json.dump(config_file, f, indent=4)