Genetic speed optimization

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
 

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")
+                ])

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,28 +628,36 @@ 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=[
                     RandomPlateauMutation(config=config, constraints_list=constraints_list),
                     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.')

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.")
 

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
 

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.: