diff --git a/WeatherRoutingTool/algorithms/genetic/__init__.py b/WeatherRoutingTool/algorithms/genetic/__init__.py index 3324a4c..c514fb0 100644 --- a/WeatherRoutingTool/algorithms/genetic/__init__.py +++ b/WeatherRoutingTool/algorithms/genetic/__init__.py @@ -6,11 +6,9 @@ import cartopy.crs as ccrs import matplotlib.pyplot as plt import numpy as np -import seaborn as sns from astropy import units as u from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.core.result import Result -from pymoo.optimize import minimize from pymoo.termination import get_termination from pymoo.util.running_metric import RunningMetric diff --git a/WeatherRoutingTool/algorithms/genetic/crossover.py b/WeatherRoutingTool/algorithms/genetic/crossover.py index d46cf86..c525bf0 100644 --- a/WeatherRoutingTool/algorithms/genetic/crossover.py +++ b/WeatherRoutingTool/algorithms/genetic/crossover.py @@ -1,20 +1,22 @@ -from pymoo.core.crossover import Crossover - -import numpy as np - -from datetime import datetime import logging import random +from copy import deepcopy +from datetime import datetime +from math import ceil + +import numpy as np +from geographiclib.geodesic import Geodesic +from pymoo.core.crossover import Crossover from WeatherRoutingTool.constraints.constraints import ConstraintsList from WeatherRoutingTool.algorithms.genetic import utils from WeatherRoutingTool.config import Config -from WeatherRoutingTool.algorithms.genetic import patcher - from WeatherRoutingTool.algorithms.genetic.patcher import PatchFactory logger = logging.getLogger("WRT.genetic.crossover") +geod = Geodesic.WGS84 + # base classes # ---------- @@ -125,8 +127,7 @@ def crossover( p1: np.ndarray, p2: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: - """Sub-class' implementation of the crossover function.""" - + """Subclass' implementation of the crossover function.""" return p1, p2 def route_constraint_violations(self, route: np.ndarray) -> np.ndarray: @@ -240,28 +241,81 @@ def print_crossover_statistics(self): opt.print_crossover_statistics() +class SpeedCrossover(OffspringRejectionCrossover): + """ + Ship speed crossover class. + Crossover candidates are identified by finding points between the parents with a distance below a specified + threshold. For a specified percentage of these candidates speed values are swapped. + """ + + def __init__(self, **kw): + # for now, we don't want to allow repairing routes for speed crossover + config = deepcopy(kw['config']) + config.GENETIC_REPAIR_TYPE = ["no_repair"] + kw['config'] = config + super().__init__(**kw) + self.threshold = 50000 # in m + self.percentage = 0.5 + + def crossover( + self, + p1: np.ndarray, + p2: np.ndarray + ) -> tuple[np.ndarray, np.ndarray]: + # Find points between parents with a distance below the specified threshold. + # There should always be one candidate (source). The destination has to be ignored. + crossover_candidates = [] + for m in range(0, len(p1)-1): + coord1 = p1[m, 0:2] + for n in range(0, len(p2)-1): + coord2 = p2[n, 0:2] + d = geod.Inverse(coord1[0], coord1[1], coord2[0], coord2[1])["s12"] + if d < self.threshold: + crossover_candidates.append((m, n)) + # Swap speed values for a subset of candidate points + indices = random.sample(range(0, len(crossover_candidates)), ceil(self.percentage*len(crossover_candidates))) + for idx in indices: + c = crossover_candidates[idx] + speed1 = p1[c[0], -1] + p1[c[0], -1] = p2[c[1], -1] + p2[c[1], -1] = speed1 + return p1, p2 + + # factory # ---------- class CrossoverFactory: @staticmethod def get_crossover(config: Config, constraints_list: ConstraintsList): - # inputs departure_time = config.DEPARTURE_TIME - return RandomizedCrossoversOrchestrator( - opts=[ - TwoPointCrossover( - config=config, - patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton", - departure_time=departure_time, - constraints_list=constraints_list, - prob=.5, - crossover_type="TP crossover", ), - SinglePointCrossover( - config=config, - patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton", - departure_time=departure_time, - constraints_list=constraints_list, - prob=.5, - crossover_type="SP crossover", ), - ], ) + # FIXME: add exception for bad combinations (better do this on the Config) + + if config.GENETIC_CROSSOVER_TYPE == "speed": + logger.debug('Setting crossover type of genetic algorithm to "speed".') + return SpeedCrossover( + config=config, + departure_time=departure_time, + constraints_list=constraints_list, + prob=.5, + crossover_type="Speed crossover") + + if config.GENETIC_CROSSOVER_TYPE == "random": + logger.debug('Setting crossover type of genetic algorithm to "random".') + return RandomizedCrossoversOrchestrator( + opts=[ + TwoPointCrossover( + config=config, + patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton", + departure_time=departure_time, + constraints_list=constraints_list, + prob=.5, + crossover_type="TP crossover"), + SinglePointCrossover( + config=config, + patch_type=config.GENETIC_CROSSOVER_PATCHER + "_singleton", + departure_time=departure_time, + constraints_list=constraints_list, + prob=.5, + crossover_type="SP crossover") + ]) diff --git a/WeatherRoutingTool/algorithms/genetic/mutation.py b/WeatherRoutingTool/algorithms/genetic/mutation.py index 8c8e47c..d33f3af 100644 --- a/WeatherRoutingTool/algorithms/genetic/mutation.py +++ b/WeatherRoutingTool/algorithms/genetic/mutation.py @@ -1,7 +1,8 @@ -import copy import logging import math import os +import random +from operator import add, sub import cartopy.crs as ccrs import matplotlib.pyplot as plt @@ -548,6 +549,76 @@ def print_mutation_statistics(self): opt.print_mutation_statistics() +class RandomPercentageChangeSpeedMutation(MutationConstraintRejection): + """ + Ship speed mutation class. + Speed values are mutated by randomly adding or subtracting a percentage. The percentage is randomly chosen + between 0 and a fixed maximum percentage (20 %). + """ + n_updates: int + config: Config + + def __init__(self, n_updates: int = 10, **kw): + super().__init__( + mutation_type="RandomPercentageChangeSpeedMutation", + **kw + ) + self.n_updates = n_updates + self.change_percent_max = 0.2 + + def mutate(self, problem, rt, **kw): + try: + indices = random.sample(range(0, rt.shape[0] - 1), self.n_updates) + except ValueError: + indices = range(0, rt.shape[0] - 1) + ops = (add, sub) + for i in indices: + op = random.choice(ops) + change_percent = random.uniform(0.0, self.change_percent_max) + new = op(rt[i][2], change_percent * rt[i][2]) + if new < 0: + new = 0 + elif new > self.config.BOAT_SPEED_MAX: + new = self.config.BOAT_SPEED_MAX + rt[i][2] = new + return rt + + +class GaussianSpeedMutation(MutationConstraintRejection): + """ + Ship speed mutation class. + Speed values are updated by drawing random samples from a Gaussian distribution. The mean value of the distribution + is half of the maximum boat speed. The standard deviation is 1/6 of the maximum boat speed. + """ + n_updates: int + config: Config + + def __init__(self, n_updates: int = 10, **kw): + super().__init__( + mutation_type="GaussianSpeedMutation", + **kw + ) + self.n_updates = n_updates + # FIXME: these numbers should be carefully evaluated + # ~99.7 % in interval (0, BOAT_SPEED_MAX) + self.mu = 0.5 * self.config.BOAT_SPEED_MAX + self.sigma = self.config.BOAT_SPEED_MAX / 6 + + def mutate(self, problem, rt, **kw): + try: + indices = random.sample(range(0, rt.shape[0] - 1), self.n_updates) + except ValueError: + indices = range(0, rt.shape[0] - 1) + for i in indices: + new = random.normalvariate(self.mu, self.sigma) + if new < 0: + new = 0 + elif new > self.config.BOAT_SPEED_MAX: + new = self.config.BOAT_SPEED_MAX + rt[i][2] = new + return rt + + # factory # ---------- class MutationFactory: @@ -557,11 +628,11 @@ def get_mutation( constraints_list: None ) -> Mutation: - if "no_mutation" in config.GENETIC_MUTATION_TYPE: + if config.GENETIC_MUTATION_TYPE == "no_mutation": logger.debug('Setting mutation type of genetic algorithm to "no_mutation".') return NoMutation() - if "random" in config.GENETIC_MUTATION_TYPE: + if config.GENETIC_MUTATION_TYPE == "random": logger.debug('Setting mutation type of genetic algorithm to "random".') return RandomMutationsOrchestrator( opts=[ @@ -569,16 +640,24 @@ def get_mutation( RouteBlendMutation(config=config, constraints_list=constraints_list) ], ) - if "rndm_walk" in config.GENETIC_MUTATION_TYPE: + if config.GENETIC_MUTATION_TYPE == "rndm_walk": logger.debug('Setting mutation type of genetic algorithm to "random_walk".') return RandomWalkMutation(config=config, constraints_list=constraints_list) - if "rndm_plateau" in config.GENETIC_MUTATION_TYPE: + if config.GENETIC_MUTATION_TYPE == "rndm_plateau": logger.debug('Setting mutation type of genetic algorithm to "random_plateau".') return RandomPlateauMutation(config=config, constraints_list=constraints_list) - if "route_blend" in config.GENETIC_MUTATION_TYPE: + if config.GENETIC_MUTATION_TYPE == "route_blend": logger.debug('Setting mutation type of genetic algorithm to "route_blend".') return RouteBlendMutation(config=config, constraints_list=constraints_list) + if config.GENETIC_MUTATION_TYPE == "percentage_change_speed": + logger.debug('Setting mutation type of genetic algorithm to "percentage_change_speed".') + return RandomPercentageChangeSpeedMutation(config=config, constraints_list=constraints_list) + + if config.GENETIC_MUTATION_TYPE == "gaussian_speed": + logger.debug('Setting mutation type of genetic algorithm to "gaussian_speed".') + return GaussianSpeedMutation(config=config, constraints_list=constraints_list) + raise NotImplementedError(f'The mutation type {config.GENETIC_MUTATION_TYPE} is not implemented.') diff --git a/WeatherRoutingTool/algorithms/genetic/patcher.py b/WeatherRoutingTool/algorithms/genetic/patcher.py index 9fa2e6a..97573bf 100644 --- a/WeatherRoutingTool/algorithms/genetic/patcher.py +++ b/WeatherRoutingTool/algorithms/genetic/patcher.py @@ -28,13 +28,15 @@ class PatcherBase: def __init__(self, *args, **kwargs): pass - def patch(self, src: tuple, dst: tuple): + def patch(self, src: tuple, dst: tuple, departure_time: datetime = None): """Obtain waypoints between `src` and `dst`. :param src: Source coords as (lat, lon) :type src: tuple[float, float] :param dst: Destination coords as (lat, lon) :type dst: tuple[float, float] + :param departure_time: Departure time + :type departure_time: datetime """ raise NotImplementedError("This patching method is not implemented.") diff --git a/WeatherRoutingTool/algorithms/genetic/population.py b/WeatherRoutingTool/algorithms/genetic/population.py index 1796a48..eba4a11 100644 --- a/WeatherRoutingTool/algorithms/genetic/population.py +++ b/WeatherRoutingTool/algorithms/genetic/population.py @@ -1,7 +1,7 @@ import logging import os -import os.path from math import ceil +from re import match import astropy.units as u import numpy as np @@ -198,21 +198,33 @@ def generate(self, problem, n_samples, **kw): logger.debug(f"Population from geojson routes: {self.routes_dir}") # routes are expected to be named in the following format: - # route_{1..N}.json + # route_{1..N}.json (geojson extension is also possible) # example: route_1.json, route_2.json, route_3.json, ... + # FIXME: add test in config.py and raise exception depending on configuration (not only speed optimization...) + X = np.full((n_samples, 1), None, dtype=object) - for i in range(n_samples): - path = os.path.join(self.routes_dir, f"route_{i + 1}.json") + files = [] + for file in os.listdir(self.routes_dir): + if match(r"route_[0-9]+\.(json|geojson)$", file.lower()): + files.append(file) + + if len(files) == 0: + raise ValueError(f"Couldn't find any route in {self.routes_dir} for the initial population.") + for i, file in enumerate(files): + path = os.path.join(self.routes_dir, file) if not os.path.exists(path): - raise ValueError("The number of available routes for the initial population does not match the " - "population size.") + raise ValueError(f"Couldn't read route {path} for the initial population.") else: route = utils.route_from_geojson_file(path) + X[i, 0] = np.array(route) - X[i, 0] = np.array(route) + added_routes = len(files) + while added_routes < n_samples: + X[added_routes, 0] = np.copy(X[0, 0]) + added_routes += 1 return X diff --git a/WeatherRoutingTool/config.py b/WeatherRoutingTool/config.py index 3493a22..a190c4f 100644 --- a/WeatherRoutingTool/config.py +++ b/WeatherRoutingTool/config.py @@ -60,6 +60,7 @@ class Config(BaseModel): BOAT_TYPE: Literal['CBT', 'SAL', 'speedy_isobased', 'direct_power_method'] = 'direct_power_method' BOAT_SPEED: float = -99. # boat speed [m/s] + BOAT_SPEED_MAX: float = 10 # maximum possible boat speed [m/s] CONSTRAINTS_LIST: List[Literal[ 'land_crossing_global_land_mask', 'land_crossing_polygons', 'seamarks', 'water_depth', 'on_map', 'via_waypoints', 'status_error' @@ -105,8 +106,9 @@ class Config(BaseModel): 'waypoints_infill', 'constraint_violation', 'no_repair' ]] = ["waypoints_infill", "constraint_violation"] GENETIC_MUTATION_TYPE: Literal[ - 'random', 'rndm_walk', 'rndm_plateau', 'route_blend', 'no_mutation' + 'random', 'rndm_walk', 'rndm_plateau', 'route_blend', 'percentage_change_speed', 'gaussian_speed', 'no_mutation' ] = 'random' + GENETIC_CROSSOVER_TYPE: Literal['random', 'speed'] = 'random' GENETIC_CROSSOVER_PATCHER: Literal['gcr', 'isofuel'] = 'isofuel' GENETIC_FIX_RANDOM_SEED: bool = False @@ -217,7 +219,8 @@ def parse_and_validate_datetime(cls, v): except ValueError: raise ValueError("'DEPARTURE_TIME' must be in format YYYY-MM-DDTHH:MMZ") - @field_validator('COURSES_FILE', 'ROUTE_PATH', 'DIJKSTRA_MASK_FILE', mode='after') + @field_validator('COURSES_FILE', 'ROUTE_PATH', 'DIJKSTRA_MASK_FILE', 'GENETIC_POPULATION_PATH', + mode='after') @classmethod def validate_path_exists(cls, v, info: ValidationInfo): if info.field_name == 'COURSES_FILE': @@ -230,6 +233,11 @@ def validate_path_exists(cls, v, info: ValidationInfo): return v else: path = Path(v) + elif info.field_name == 'GENETIC_POPULATION_PATH': + if info.data.get('GENETIC_POPULATION_TYPE') != 'from_geojson': + return v + else: + path = Path(v) else: path = Path(v) if not path.exists(): @@ -468,8 +476,8 @@ def check_boat_speed(cls, v): @model_validator(mode='after') def check_speed_determination(self) -> Self: - print('arrival time: ', self.ARRIVAL_TIME) - print('speed: ', self.BOAT_SPEED) + logger.info(f'arrival time: {self.ARRIVAL_TIME}') + logger.info(f'speed: {self.BOAT_SPEED}') if self.ARRIVAL_TIME == '9999-99-99T99:99Z' and self.BOAT_SPEED == -99.: raise ValueError('Please specify either the boat speed or the arrival time') if not self.ARRIVAL_TIME == '9999-99-99T99:99Z' and not self.BOAT_SPEED == -99.: diff --git a/WeatherRoutingTool/constraints/route_postprocessing.py b/WeatherRoutingTool/constraints/route_postprocessing.py index c764816..7689491 100644 --- a/WeatherRoutingTool/constraints/route_postprocessing.py +++ b/WeatherRoutingTool/constraints/route_postprocessing.py @@ -551,7 +551,6 @@ def terminate( courses = route_dict['courses'] dists = route_dict['dist'] start_times = route_dict['start_times'] - # FIXME: boat_speed is a list or array arrival_time = start_times[-1] + timedelta(seconds=dists[-1].value / boat_speed[-1].value) travel_times = np.append(travel_times, -99 * u.second)