diff --git a/WeatherRoutingTool/algorithms/genetic/crossover.py b/WeatherRoutingTool/algorithms/genetic/crossover.py index 27bd093..d46cf86 100644 --- a/WeatherRoutingTool/algorithms/genetic/crossover.py +++ b/WeatherRoutingTool/algorithms/genetic/crossover.py @@ -217,66 +217,6 @@ def crossover(self, p1, p2): return r1, r2 -# FIXME: Adapt to continuous routing or eliminate. -class PMX(OffspringRejectionCrossover): - """Partially Mapped Crossover.""" - - def crossover(self, p1, p2): - if p1.shape != p2.shape: - logging.info("PMX — Not of equal length") - return p1, p2 - - N = min(p1.shape[0], p2.shape[0]) - - # Convert to lists of tuples - parent1 = [tuple(row) for row in p1] - parent2 = [tuple(row) for row in p2] - - # Choose crossover points - cx1, cx2 = sorted(np.random.choice(range(N), 2, replace=False)) - - # Initialize offspring placeholders - child1 = [None] * N - child2 = [None] * N - - # Copy the segment - for i in range(cx1, cx2): - child1[i] = parent2[i] - child2[i] = parent1[i] - - # Build mapping for the swapped segments - mapping12 = {parent2[i]: parent1[i] for i in range(cx1, cx2)} - mapping21 = {parent1[i]: parent2[i] for i in range(cx1, cx2)} - - def resolve(gene, segment, mapping): - # Keep resolving until gene is not in the given segment - while gene in segment: - gene = mapping[gene] - return gene - - # Fill remaining positions - for i in range(N): - if not (cx1 <= i < cx2): - g1 = parent1[i] - g2 = parent2[i] - - # If g1 is already in the swapped segment of child1, resolve via mapping12 - if g1 in child1[cx1:cx2]: - g1 = resolve(g1, child1[cx1:cx2], mapping12) - child1[i] = g1 - - # Likewise for child2 - if g2 in child2[cx1:cx2]: - g2 = resolve(g2, child2[cx1:cx2], mapping21) - child2[i] = g2 - - # Convert back to numpy arrays - c1 = np.array(child1, dtype=p1.dtype) - c2 = np.array(child2, dtype=p1.dtype) - - return c1, c2 - - # # ---------- class RandomizedCrossoversOrchestrator(CrossoverBase): diff --git a/WeatherRoutingTool/algorithms/genetic/mutation.py b/WeatherRoutingTool/algorithms/genetic/mutation.py index 0245a0c..8c8e47c 100644 --- a/WeatherRoutingTool/algorithms/genetic/mutation.py +++ b/WeatherRoutingTool/algorithms/genetic/mutation.py @@ -115,7 +115,8 @@ def _do(self, problem, X, **kw): :param problem: Routing problem. :type: RoutingProblem :param X: Route matrix in the form of ``np.array([[route_0], [route_1], ...])`` with - ``route_i=np.array([[lat_0, lon_0], [lat_1,lon_1], ...])``. X.shape = (n_routes, 1, n_waypoints, 2). + ``route_i=np.array([[lat_0, lon_0, v_0], [lat_1,lon_1, v_1], ...])``. + X.shape = (n_routes, 1, n_waypoints, 3). Access i'th route as ``X[i,0]`` and the j'th coordinate pair off the i'th route as ``X[i,0][j, :]``. :type X: np.array :return: Mutated route matrix. Same structure as for ``X``. @@ -215,10 +216,10 @@ def __init__( def random_walk( self, - point: tuple[float, float], + point: tuple[float, float, float], dist: float = 1e4, bearing: float = 45.0, - ) -> tuple[float, float]: + ) -> tuple[float, float, float]: """Pick an N4 neighbour of a waypoint. :param point: (lat, lon) in degrees. @@ -231,11 +232,11 @@ def random_walk( :rtype: tuple[float, float] """ - lat0, lon0 = point + lat0, lon0, speed = point result = Geodesic.WGS84.Direct(lat0, lon0, bearing, dist) lat2 = result["lat2"] lon2 = result["lon2"] - return lat2, lon2 + return lat2, lon2, speed def mutate(self, problem, rt, **kw): """ @@ -258,15 +259,15 @@ def mutate(self, problem, rt, **kw): :param problem: routing problem :type: RoutingProblem :params rt: route to be mutated - :type rt: np.array([[lat_0, lon_0], [lat_1,lon_1], ...]), + :type rt: np.array([[lat_0, lon_0, v_0], [lat_1,lon_1, v_1], ...]), :return: mutated route - :rtype: np.array([[lat_0, lon_0], [lat_1,lon_1], ...]), + :rtype: np.array([[lat_0, lon_0, v_0], [lat_1,lon_1, v_1], ...]), """ debug = False # test whether input route rt has the correct shape assert len(rt.shape) == 2 - assert rt.shape[1] == 2 + assert rt.shape[1] == 3 route_length = rt.shape[0] plateau_length = 2 * self.plateau_slope + self.plateau_size - 2 rt_new = np.full(rt.shape, -99.) @@ -421,7 +422,7 @@ def bezier_curve(control_points, n_points=100): control_points = np.array(control_points) n = len(control_points) - 1 # degree t = np.linspace(0, 1, n_points) - curve = np.zeros((n_points, 2)) + curve = np.zeros((n_points, 3)) for i in range(n + 1): bernstein = math.comb(n, i) * (t ** i) * ((1 - t) ** (n - i)) @@ -432,7 +433,7 @@ def bezier_curve(control_points, n_points=100): def mutate(self, problem, rt, **kw): # test shape of input route assert len(rt.shape) == 2 - assert rt.shape[1] == 2 + assert rt.shape[1] == 3 route_length = rt.shape[0] # only mutate routes that are long enough @@ -491,26 +492,26 @@ def __init__( def random_walk( self, - point: tuple[float, float], + point: tuple[float, float, float], dist: float = 1e4, bearing: float = 45.0, - ) -> tuple[float, float]: + ) -> tuple[float, float, float]: """Pick an N4 neighbour of a waypoint. :param point: (lat, lon) in degrees. - :type point: tuple[float, float] + :type point: tuple[float, float, float] :param dist: distance in meters :type dist: float :param bearing: Azimuth in degrees (clockwise from North) :type bearing: float :return: (lat, lon) in degrees. - :rtype: tuple[float, float] + :rtype: tuple[float, float, float] """ - lat0, lon0 = point + lat0, lon0, speed = point result = Geodesic.WGS84.Direct(lat0, lon0, bearing, dist) lat2 = result["lat2"] lon2 = result["lon2"] - return lat2, lon2 + return lat2, lon2, speed def mutate(self, problem, rt, **kw): for _ in range(self.n_updates): diff --git a/WeatherRoutingTool/algorithms/genetic/patcher.py b/WeatherRoutingTool/algorithms/genetic/patcher.py index c196a01..9fa2e6a 100644 --- a/WeatherRoutingTool/algorithms/genetic/patcher.py +++ b/WeatherRoutingTool/algorithms/genetic/patcher.py @@ -62,6 +62,8 @@ def __call__(cls, *args, **kwargs): class GreatCircleRoutePatcher(PatcherBase): """Produce a set of waypoints along the Great Circle Route between src and dst. + The same speed as the speed at `src` is added to every waypoint. + :param dist: Dist between each waypoint in the Great Circle Route :type dist: float """ @@ -73,20 +75,26 @@ def __init__(self, dist: float = 10_000.0): # variables self.dist = dist - def patch(self, src: tuple, dst: tuple, departure_time: datetime = None, npoints=None, ) -> np.ndarray: + def patch(self, + src: tuple[float, float, float], + dst: tuple[float, float, float], + departure_time: datetime = None, + npoints=None, + ) -> np.ndarray: """Generate equi-distant waypoints across the Great Circle Route from src to dst - :param src: Source waypoint as (lat, lon) pair - :type src: tuple[float, float] - :param dst: Destination waypoint as (lat, lon) pair - :type dst: tuple[float, float] - :return: List of waypoints along the great circle (lat, lon) - :rtype: np.array[tuple[float, float]] + :param src: Source waypoint as (lat, lon, v) triple + :type src: tuple[float, float, float] + :param dst: Destination waypoint as (lat, lon, v) triple + :type dst: tuple[float, float, float] + :return: List of waypoints along the great circle (lat, lon, v) + :rtype: np.array[tuple[float, float, float]] """ geod: Geodesic = Geodesic.WGS84 - line = geod.InverseLine(*src, *dst) + line = geod.InverseLine(*src[:-1], *dst[:-1]) + speed = src[2] if not npoints == None: self.dist = line.s13 / npoints @@ -97,7 +105,7 @@ def patch(self, src: tuple, dst: tuple, departure_time: datetime = None, npoints for i in range(npoints + 1): s = min(self.dist * i, line.s13) g = line.Position(s, Geodesic.STANDARD | Geodesic.LONG_UNROLL) - route.append((g['lat2'], g['lon2'])) + route.append((g['lat2'], g['lon2'], speed)) return np.array([src, *route[1:-1], dst]) @@ -256,14 +264,20 @@ def _setup_components(self) -> tuple[WeatherCond, Boat, WaterDepth, ConstraintsL return wt, boat, water_depth, constraints_list - def patch(self, src, dst, departure_time: datetime = None): + def patch(self, + src: tuple[float, float, float], + dst: tuple[float, float, float], + departure_time: datetime = None + ): """ Produce a set of waypoints between src and dst using the IsoFuel algorithm. - :param src: Source waypoint as (lat, lon) pair - :type src: tuple[float, float] - :param dst: Destination waypoint as (lat, lon) pair - :type dst: tuple[float, float] + The same speed as the speed at `src` is added to every waypoint. + + :param src: Source waypoint as (lat, lon, speed) triple + :type src: tuple[float, float, float] + :param dst: Destination waypoint as (lat, lon, speed) triple + :type dst: tuple[float, float, float] :param departure_time: departure time from src :type departure_time: datetime :return: List of waypoints or list of multiple routes connecting src and dst @@ -272,7 +286,7 @@ def patch(self, src, dst, departure_time: datetime = None): self.patch_count += 1 cfg = self.config.model_copy(update={ - "DEFAULT_ROUTE": [*src, *dst], + "DEFAULT_ROUTE": [*src[:-1], *dst[:-1]], "DEPARTURE_TIME": departure_time }) @@ -303,9 +317,11 @@ def patch(self, src, dst, departure_time: datetime = None): logger.debug('Falling back to gcr patching!') return self.patchfn_gcr.patch(src, dst, departure_time) + speed = np.full(min_fuel_route.lons_per_step.shape, src[2]) + # single route if self.n_routes == "single": - return np.stack([min_fuel_route.lats_per_step, min_fuel_route.lons_per_step], axis=1) + return np.stack([min_fuel_route.lats_per_step, min_fuel_route.lons_per_step, speed], axis=1) # list of routes if not alg.route_list: @@ -314,7 +330,7 @@ def patch(self, src, dst, departure_time: datetime = None): routes = [] for rt in alg.route_list: - routes.append(np.stack([rt.lats_per_step, rt.lons_per_step], axis=1)) + routes.append(np.stack([rt.lats_per_step, rt.lons_per_step, speed], axis=1)) return routes diff --git a/WeatherRoutingTool/algorithms/genetic/population.py b/WeatherRoutingTool/algorithms/genetic/population.py index df646d0..1796a48 100644 --- a/WeatherRoutingTool/algorithms/genetic/population.py +++ b/WeatherRoutingTool/algorithms/genetic/population.py @@ -3,6 +3,7 @@ import os.path from math import ceil +import astropy.units as u import numpy as np from geographiclib.geodesic import Geodesic from pymoo.core.sampling import Sampling @@ -25,7 +26,7 @@ class Population(Sampling): """Base Class for generating the initial population.""" - def __init__(self, default_route: list, constraints_list: list, pop_size: int): + def __init__(self, config: Config, default_route: list, constraints_list: ConstraintsList, pop_size: int): super().__init__() self.constraints_list = constraints_list @@ -35,12 +36,20 @@ def __init__(self, default_route: list, constraints_list: list, pop_size: int): self.src: tuple[float, float] = tuple(default_route[:-2]) self.dst: tuple[float, float] = tuple(default_route[-2:]) + self.departure_time = config.DEPARTURE_TIME + self.arrival_time = config.ARRIVAL_TIME + self.boat_speed = config.BOAT_SPEED * u.meter / u.second + + self.boat_speed_from_arrival_time = False + if self.boat_speed.value == -99.: + self.boat_speed_from_arrival_time = True + def _do(self, problem, n_samples, **kw): X = self.generate(problem, n_samples, **kw) for rt, in X: - assert tuple(rt[0]) == self.src, "Source waypoint not matching" - assert tuple(rt[-1]) == self.dst, "Destination waypoint not matching" + assert tuple(rt[0, :-1]) == self.src, "Source waypoint not matching" + assert tuple(rt[-1, :-1]) == self.dst, "Destination waypoint not matching" self.X = X return self.X @@ -71,6 +80,19 @@ def check_validity(self, routes): f"{self.n_constrained_routes} / {self.pop_size} constrained — " "More than 50% of the initial routes are constrained") + def recalculate_speed_for_route(self, rt): + """ + Recalculate speed at the waypoints if no `BOAT_SPEED` but an `ARRIVAL_TIME` is specified. + """ + bs = utils.get_speed_from_arrival_time( + lons=rt[:, 1], + lats=rt[:, 0], + departure_time=self.departure_time, + arrival_time=self.arrival_time, + ) + rt[:, 2] = np.full(rt[:, 1].shape, bs) + return rt + class GridBasedPopulation(GridMixin, Population): """Make initial population for genetic algorithm based on a grid and associated cost values @@ -81,8 +103,14 @@ class GridBasedPopulation(GridMixin, Population): - call `print(GridBasedPopulation.mro())` to see the method resolution order """ - def __init__(self, default_route, grid, constraints_list, pop_size): - super().__init__(default_route=default_route, grid=grid, constraints_list=constraints_list, pop_size=pop_size) + def __init__(self, config: Config, default_route, grid, constraints_list, pop_size): + super().__init__( + config=config, + default_route=default_route, + grid=grid, + constraints_list=constraints_list, + pop_size=pop_size + ) # update nan_mask with constraints_list # ---------- @@ -110,6 +138,10 @@ def generate(self, problem, n_samples, **kw): _, _, start_indices = self.coords_to_index([self.src]) _, _, end_indices = self.coords_to_index([self.dst]) + boat_speed = self.boat_speed + if self.boat_speed_from_arrival_time: + boat_speed = 6 * u.meter / u.second # dummy boat speed + for i in range(n_samples): shuffled_cost = self.get_shuffled_cost() @@ -122,10 +154,18 @@ def generate(self, problem, n_samples, **kw): # logger.debug(f"GridBasedPopulation._do: type(route)={type(route)}, route={route}") _, _, route = self.index_to_coords(route) + route = np.array(route) + speed_arr = np.full((route.shape[0], 1), boat_speed.value) + route = np.hstack((route, speed_arr)) # match first and last points to src and dst + src_speed = np.array(self.src + (boat_speed.value,)) + dst_speed = np.array(self.dst + (boat_speed.value,)) + if self.boat_speed_from_arrival_time: + route = self.recalculate_speed_for_route(route) + X[i, 0] = np.array([ - self.src, *route[1:-1], self.dst]) + src_speed, *route[1:-1], dst_speed]) return X @@ -142,8 +182,13 @@ class FromGeojsonPopulation(Population): :type routes_dir: str """ - def __init__(self, routes_dir: str, default_route, constraints_list, pop_size): - super().__init__(default_route=default_route, constraints_list=constraints_list, pop_size=pop_size) + def __init__(self, config: Config, routes_dir: str, default_route, constraints_list, pop_size): + super().__init__( + config=config, + default_route=default_route, + constraints_list=constraints_list, + pop_size=pop_size + ) if not os.path.exists(routes_dir) or not os.path.isdir(routes_dir): raise FileNotFoundError("Routes directory not found") @@ -167,9 +212,6 @@ def generate(self, problem, n_samples, **kw): else: route = utils.route_from_geojson_file(path) - assert np.array_equal(route[0], self.src), "Route not starting at source" - assert np.array_equal(route[-1], self.dst), "Route not ending at destination" - X[i, 0] = np.array(route) return X @@ -184,24 +226,49 @@ class IsoFuelPopulation(Population): produced route is repeated until the required number of individuals are met """ - def __init__(self, config: Config, boat: Boat, default_route, constraints_list, pop_size): + def __init__(self, config: Config, default_route, constraints_list, pop_size): super().__init__( + config=config, default_route=default_route, constraints_list=constraints_list, - pop_size=pop_size, ) - - self.departure_time = config.DEPARTURE_TIME + pop_size=pop_size, + ) self.patcher = PatchFactory.get_patcher(config=config, patch_type="isofuel_multiple_routes", application="initial population") def generate(self, problem, n_samples, **kw): - routes = self.patcher.patch(self.src, self.dst, self.departure_time) + """Generate the initial population. + + Calls the :py:class:`IsofuelPatcher` to patch + routes from the start coordinates to the destination. In case an `ARRIVAL_TIME` is specified, a dummy boat speed + is passed that is later on recalculated by `recalculate_speed_for_route`. If the number `n_samples` of routes + can not be provided by the patcher, the last route that can be provided is copied until the requested number of + routes has been achieved. + + :params problem: routing problem + :type problem: Problem + :params n_samples: number of routes for the initial population + :type n_samples: int + :return: Route matrix in the form of ``np.array([[route_0], [route_1], ...])`` with + ``route_i=np.array([[lat_0, lon_0, v_0], [lat_1,lon_1, v_1], ...])``. + X.shape = (n_routes, 1, n_waypoints, 3). + Access i'th route as ``X[i,0]`` and the j'th coordinate pair off the i'th route as ``X[i,0][j, :]``. + :rtype: np.array + """ + + boat_speed = self.boat_speed + if self.boat_speed_from_arrival_time: + boat_speed = 6 * u.meter / u.second # add dummy speed, will be recalculated + routes = self.patcher.patch(self.src + (boat_speed.value,), self.dst + (boat_speed.value,), self.departure_time) X = np.full((n_samples, 1), None, dtype=object) for i, rt in enumerate(routes): - X[i, 0] = np.array([self.src, *rt[1:-1], self.dst]) + if self.boat_speed_from_arrival_time: + rt = self.recalculate_speed_for_route(rt) + + X[i, 0] = rt # fallback: fill all other individuals with the same population as the last one for j in range(i + 1, n_samples): @@ -212,7 +279,13 @@ def generate(self, problem, n_samples, **kw): class GcrSliderPopulation(Population): def __init__(self, config: Config, default_route, constraints_list, pop_size): - super().__init__(default_route=default_route, constraints_list=constraints_list, pop_size=pop_size) + super().__init__( + config=config, + default_route=default_route, + constraints_list=constraints_list, + pop_size=pop_size + ) + self.algo = GcrSliderAlgorithm(config) def generate(self, problem, n_samples, **kw): @@ -227,7 +300,11 @@ def generate(self, problem, n_samples, **kw): distance used to move the point incrementally. """ # FIXME: how to handle already existing waypoints specified for the genetic algorithm? - route = self.create_route() + boat_speed = self.boat_speed + if self.boat_speed_from_arrival_time: + boat_speed = 6 * u.meter / u.second # dummy boat speed + + route = self.create_route(speed=boat_speed.value) routes = [] if route is not None: routes.append(route) @@ -235,18 +312,18 @@ def generate(self, problem, n_samples, **kw): line = geod.InverseLine(self.algo.start[0], self.algo.start[1], self.algo.finish[0], self.algo.finish[1]) wpt_increment_max = 0.5 * line.s13 wpt_increment = 0.05 * line.s13 - wpt_increment_steps_max = ceil(wpt_increment_max/wpt_increment) + wpt_increment_steps_max = ceil(wpt_increment_max / wpt_increment) element = 1 clockwise = True wpt_increment_step = 1 while len(routes) < n_samples: dist_fraction = self.van_der_corput_sequence(element) - g = line.Position(dist_fraction*line.s13, Geodesic.STANDARD | Geodesic.LONG_UNROLL) + g = line.Position(dist_fraction * line.s13, Geodesic.STANDARD | Geodesic.LONG_UNROLL) dist_orthogonal = wpt_increment_step * wpt_increment lat, lon = self.algo.move_point_orthogonally(g, dist_orthogonal, clockwise=clockwise) if not self.algo.is_land(lat, lon): - route = self.create_route(lat, lon) + route = self.create_route(lat, lon, boat_speed.value) if route is not None: routes.append(route) logger.info(f"Found {len(routes)} of {n_samples} routes for initial population.") @@ -264,14 +341,16 @@ def generate(self, problem, n_samples, **kw): X = np.full((n_samples, 1), None, dtype=object) for i, rt in enumerate(routes): - X[i, 0] = np.array([self.src, *rt[1:-1], self.dst]) + if self.boat_speed_from_arrival_time: + rt = self.recalculate_speed_for_route(rt) + X[i, 0] = rt # fallback: fill all other individuals with the same population as the last one for j in range(i + 1, n_samples): X[j, 0] = np.copy(X[j - 1, 0]) return X - def create_route(self, lat: float = None, lon: float = None): + def create_route(self, lat: float = None, lon: float = None, speed: float = None): """ :param lat: latitude of the waypoint :type lat: float @@ -286,7 +365,8 @@ def create_route(self, lat: float = None, lon: float = None): # import uuid # filename = f"{str(uuid.uuid4())}.geojson" # route.write_to_geojson(filename) - route = [[route.lats_per_step[i], route.lons_per_step[i]] for i in range(0, len(route.lats_per_step))] + route = [[route.lats_per_step[i], route.lons_per_step[i], speed] for i in + range(0, len(route.lats_per_step))] route = np.array(route) except Exception: pass @@ -330,6 +410,7 @@ def get_population( match population_type: case "grid_based": return GridBasedPopulation( + config=config, grid=wave_height, default_route=config.DEFAULT_ROUTE, constraints_list=constraints_list, @@ -338,13 +419,13 @@ def get_population( case "isofuel": return IsoFuelPopulation( config=config, - boat=boat, default_route=config.DEFAULT_ROUTE, constraints_list=constraints_list, pop_size=population_size, ) case "from_geojson": return FromGeojsonPopulation( + config=config, routes_dir=config.GENETIC_POPULATION_PATH, default_route=config.DEFAULT_ROUTE, constraints_list=constraints_list, diff --git a/WeatherRoutingTool/algorithms/genetic/problem.py b/WeatherRoutingTool/algorithms/genetic/problem.py index ea1bb5b..73a6397 100644 --- a/WeatherRoutingTool/algorithms/genetic/problem.py +++ b/WeatherRoutingTool/algorithms/genetic/problem.py @@ -48,19 +48,12 @@ def get_power(self, route): bs = self.boat_speed if self.boat_speed_from_arrival_time: - dummy_speed = 6 * u.meter / u.second - route_dict = RouteParams.get_per_waypoint_coords( - route[:, 1], - route[:, 0], - self.departure_time, - dummy_speed, ) - - full_travel_distance = np.sum(route_dict['dist']) - print('self.arrival_time: ', self.arrival_time) - print('self.departure_time: ', self.departure_time) - - time_diff = self.arrival_time - self.departure_time - bs = full_travel_distance / (time_diff.total_seconds() * u.second) + bs = utils.get_speed_from_arrival_time( + lons=route[:, 1], + lats=route[:, 0], + departure_time=self.departure_time, + arrival_time=self.arrival_time, + ) route_dict = RouteParams.get_per_waypoint_coords( route[:, 1], diff --git a/WeatherRoutingTool/algorithms/genetic/repair.py b/WeatherRoutingTool/algorithms/genetic/repair.py index a5b1a6e..19e6458 100644 --- a/WeatherRoutingTool/algorithms/genetic/repair.py +++ b/WeatherRoutingTool/algorithms/genetic/repair.py @@ -138,7 +138,7 @@ def repair_single_route(self, rt, patchfn, constrained): """ # check for correct input shape of rt assert len(rt.shape) == 2 - assert rt.shape[1] == 2 + assert rt.shape[1] == 3 debug = False prev_seg_end = 0 diff --git a/WeatherRoutingTool/algorithms/genetic/utils.py b/WeatherRoutingTool/algorithms/genetic/utils.py index fab16cf..9bbb95f 100644 --- a/WeatherRoutingTool/algorithms/genetic/utils.py +++ b/WeatherRoutingTool/algorithms/genetic/utils.py @@ -1,14 +1,13 @@ +import json import logging +from typing import Optional +import astropy.units as u +import numpy as np +from geographiclib.geodesic import Geodesic from pymoo.core.duplicate import ElementwiseDuplicateElimination -from geographiclib.geodesic import Geodesic -import numpy as np - -from typing import Optional -import functools -import json -import math +from WeatherRoutingTool.routeparams import RouteParams logger = logging.getLogger("WRT.genetic") @@ -26,7 +25,7 @@ def gcr_distance(src, dst) -> float: geod = Geodesic.WGS84 - rs = geod.Inverse(*src, *dst) + rs = geod.Inverse(*src[:-1], *dst[:-1]) return rs["s12"] @@ -126,10 +125,13 @@ def route_from_geojson(dt: dict) -> list[tuple[float, float]]: :rtype: list[tuple[float, float]] """ - route = [ - ft["geometry"]["coordinates"][::-1] - for ft in dt["features"] + waypoints = [ + ft["geometry"]["coordinates"][::-1] for ft in dt["features"] + ] + speed_info = [ + [ft["properties"]["speed"]["value"]] for ft in dt["features"] ] + route = np.hstack((waypoints, speed_info)) return route @@ -149,6 +151,36 @@ def route_from_geojson_file(path: str) -> list[tuple[float, float]]: return route_from_geojson(dt) +def get_speed_from_arrival_time(lons, lats, departure_time, arrival_time): + """ + Calculate boat speed based on coordinates, departure and arrival time for a route array. + + :param lons: longitudes + :type lons: np.array + :param lats: latitudes + :type lats: np.array + :param departure_time: departure time + :type departure_time: np.array of datetime objects + :param arrival_time: arrival time + :type arrival_time: datetime object + :return: array of boat speeds + :rtype: np.array + + """ + dummy_speed = 6 * u.meter / u.second + route_dict = RouteParams.get_per_waypoint_coords( + lons, + lats, + departure_time, + dummy_speed, ) + + full_travel_distance = np.sum(route_dict['dist']) + + time_diff = arrival_time - departure_time + bs = full_travel_distance / (time_diff.total_seconds() * u.second) + return bs + + # ---------- class RouteDuplicateElimination(ElementwiseDuplicateElimination): """Custom duplicate elimination strategy for routing problem.""" diff --git a/WeatherRoutingTool/algorithms/routingalg.py b/WeatherRoutingTool/algorithms/routingalg.py index 653c15e..03db1db 100644 --- a/WeatherRoutingTool/algorithms/routingalg.py +++ b/WeatherRoutingTool/algorithms/routingalg.py @@ -53,6 +53,8 @@ def __init__(self, config): self.boat_speed = config.BOAT_SPEED * u.meter/u.second def get_boat_speed(self, dists=None): + if self.boat_speed == -99: + return None return self.boat_speed def init_fig(self, **kwargs): diff --git a/requirements.test.txt b/requirements.test.txt index 1b4a4f8..3fd0ece 100644 --- a/requirements.test.txt +++ b/requirements.test.txt @@ -1,4 +1,5 @@ -r requirements.txt flake8 psycopg2 -pytest \ No newline at end of file +pytest +pytest-plt diff --git a/tests/test_genetic.py b/tests/test_genetic.py index c210569..8a6ac01 100644 --- a/tests/test_genetic.py +++ b/tests/test_genetic.py @@ -5,34 +5,31 @@ import cartopy.crs as ccrs import numpy as np -import matplotlib.pyplot as plt +import matplotlib.pyplot as pyplot +import pytest from astropy import units as u +from matplotlib.collections import LineCollection +from matplotlib.colors import Normalize import tests.basic_test_func as basic_test_func import WeatherRoutingTool.utils.graphics as graphics +from WeatherRoutingTool.algorithms.genetic.crossover import SinglePointCrossover from WeatherRoutingTool.algorithms.genetic.patcher import PatcherBase, GreatCircleRoutePatcher, IsofuelPatcher, \ GreatCircleRoutePatcherSingleton, IsofuelPatcherSingleton, PatchFactory +from WeatherRoutingTool.algorithms.genetic.population import IsoFuelPopulation from WeatherRoutingTool.algorithms.genetic.mutation import RandomPlateauMutation, RouteBlendMutation from WeatherRoutingTool.config import Config from WeatherRoutingTool.algorithms.genetic.repair import ConstraintViolationRepair from WeatherRoutingTool.ship.ship_config import ShipConfig from WeatherRoutingTool.utils.maps import Map -# FIXME: the following test functions fail if LaTeX is not installed: -# - tests/test_genetic.py::test_random_plateau_mutation -# - tests/test_genetic.py::test_bezier_curve_mutation -# - tests/test_genetic.py::test_constraint_violation_repair -# In the GH Actions workflow, we install the packages texlive, texlive-latex-extra and cm-super to make sure the -# tests are passing. However, this leads to additional traffic when running the workflow. It would be better to -# exclude plotting in the tests or adapt it so that LaTeX doesn't need to be installed. - def test_isofuelpatcher_singleton(): dirname = os.path.dirname(__file__) configpath = os.path.join(dirname, 'config.isofuel_single_route.json') config = Config.assign_config(Path(configpath)) - src = [38.851, 4.066] - dst = [37.901, 8.348] + src = [38.851, 4.066, 7] + dst = [37.901, 8.348, 7] departure_time = datetime(2025, 4, 1, 12, 11) pt_one = IsofuelPatcherSingleton(config) @@ -47,8 +44,8 @@ def test_isofuelpatcher_no_singleton(): dirname = os.path.dirname(__file__) configpath = os.path.join(dirname, 'config.isofuel_single_route.json') config = Config.assign_config(Path(configpath)) - src = [38.851, 4.066] - dst = [37.901, 8.348] + src = [38.851, 4.066, 7] + dst = [37.901, 8.348, 7] departure_time = datetime(2025, 4, 1, 12, 11) pt_one = IsofuelPatcher(config) @@ -61,28 +58,28 @@ def test_isofuelpatcher_no_singleton(): def get_dummy_route_input(length='long'): route1 = np.array([ - [35.199, 15.490], - [34.804, 16.759], - [34.447, 18.381], - [34.142, 18.763], - [33.942, 21.080], - [33.542, 23.024], - [33.408, 24.389], - [33.166, 26.300], - [32.937, 27.859], - [32.737, 28.859], + [35.199, 15.490, 10], + [34.804, 16.759, 10], + [34.447, 18.381, 10], + [34.142, 18.763, 10], + [33.942, 21.080, 10], + [33.542, 23.024, 10], + [33.408, 24.389, 10], + [33.166, 26.300, 10], + [32.937, 27.859, 10], + [32.737, 28.859, 10], ]) route2 = np.array([ - [35.199, 16.490], - [34.804, 17.759], - [34.447, 19.381], - [34.142, 19.763], - [33.942, 22.080], - [33.542, 23.024], - [33.408, 24.389], - [33.166, 25.300], - [32.937, 26.859], - [32.737, 27.859], + [35.199, 16.490, 20], + [34.804, 17.759, 20], + [34.447, 19.381, 20], + [34.142, 19.763, 20], + [33.942, 22.080, 20], + [33.542, 23.024, 20], + [33.408, 24.389, 20], + [33.166, 25.300, 20], + [32.937, 26.859, 20], + [32.737, 27.859, 20], ]) if length == "short": route1 = np.delete(route1, -1, 0) @@ -102,12 +99,27 @@ def get_dummy_route_input(length='long'): ''' -def test_random_plateau_mutation(): +def get_route_lc(X): + lats = X[:, 0] + lons = X[:, 1] + speed = X[:, 2] + + points = np.array([lons, lats]).T.reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + + norm = Normalize(vmin=10, vmax=20) + lc = LineCollection(segments, cmap='viridis', norm=norm, transform=ccrs.Geodetic()) + lc.set_array(speed) + lc.set_linewidth(3) + return lc + + +@pytest.mark.manual +def test_random_plateau_mutation(plt): dirname = os.path.dirname(__file__) configpath = os.path.join(dirname, 'config.isofuel_single_route.json') config = Config.assign_config(Path(configpath)) default_map = Map(32., 15, 36, 29) - input_crs = ccrs.PlateCarree() constraint_list = basic_test_func.generate_dummy_constraint_list() np.random.seed(1) @@ -127,10 +139,20 @@ def test_random_plateau_mutation(): show_depth=False, show_gcr=False ) - ax.plot(old_route[0, 0][:, 1], old_route[0, 0][:, 0], color="firebrick", transform=input_crs) - ax.plot(new_route[0, 0][:, 1], new_route[0, 0][:, 0], color="blue", transform=input_crs) - ax.plot(old_route[1, 0][:, 1], old_route[1, 0][:, 0], color="firebrick", transform=input_crs) - ax.plot(new_route[1, 0][:, 1], new_route[1, 0][:, 0], color="blue", transform=input_crs) + old_route_one_lc = get_route_lc(old_route[0, 0]) + old_route_two_lc = get_route_lc(old_route[1, 0]) + new_route_one_lc = get_route_lc(new_route[0, 0]) + new_route_two_lc = get_route_lc(new_route[1, 0]) + ax.add_collection(old_route_one_lc) + ax.add_collection(old_route_two_lc) + ax.add_collection(new_route_one_lc) + ax.add_collection(new_route_two_lc) + + cbar = fig.colorbar(old_route_one_lc, ax=ax, orientation='vertical', pad=0.15, shrink=0.7) + cbar.set_label('Geschwindigkeit ($m/s$)') + + pyplot.tight_layout() + plt.saveas = "test_random_plateau_mutation.png" assert old_route.shape == new_route.shape for i_route in range(old_route.shape[0]): @@ -167,12 +189,12 @@ def test_random_plateau_mutation_refusal(): ''' -def test_bezier_curve_mutation(): +@pytest.mark.manual +def test_bezier_curve_mutation(plt): dirname = os.path.dirname(__file__) configpath = os.path.join(dirname, 'config.isofuel_single_route.json') config = Config.assign_config(Path(configpath)) default_map = Map(32., 15, 36, 29) - input_crs = ccrs.PlateCarree() constraint_list = basic_test_func.generate_dummy_constraint_list() np.random.seed(2) @@ -192,10 +214,20 @@ def test_bezier_curve_mutation(): show_gcr=False ) - ax.plot(old_route[0, 0][:, 1], old_route[0, 0][:, 0], color="firebrick", transform=input_crs) - ax.plot(new_route[0, 0][:, 1], new_route[0, 0][:, 0], color="blue", transform=input_crs) - ax.plot(old_route[1, 0][:, 1], old_route[1, 0][:, 0], color="firebrick", transform=input_crs) - ax.plot(new_route[1, 0][:, 1], new_route[1, 0][:, 0], color="blue", transform=input_crs) + old_route_one_lc = get_route_lc(old_route[0, 0]) + old_route_two_lc = get_route_lc(old_route[1, 0]) + new_route_one_lc = get_route_lc(new_route[0, 0]) + new_route_two_lc = get_route_lc(new_route[1, 0]) + ax.add_collection(old_route_one_lc) + ax.add_collection(old_route_two_lc) + ax.add_collection(new_route_one_lc) + ax.add_collection(new_route_two_lc) + + cbar = fig.colorbar(old_route_one_lc, ax=ax, orientation='vertical', pad=0.15, shrink=0.7) + cbar.set_label('Geschwindigkeit ($m/s$)') + + pyplot.tight_layout() + plt.saveas = "test_bezier_curve_mutation.png" assert old_route.shape == new_route.shape for i_route in range(old_route.shape[0]): @@ -246,12 +278,12 @@ def test_configuration_isofuel_patcher(): assert config_ship.BOAT_UNDER_KEEL_CLEARANCE * u.meter == pt.boat.under_keel_clearance -def test_constraint_violation_repair(): +@pytest.mark.manual +def test_constraint_violation_repair(plt): dirname = os.path.dirname(__file__) configpath = os.path.join(dirname, 'config.isofuel_single_route.json') config = Config.assign_config(Path(configpath)) default_map = Map(32., 15, 36, 29) - input_crs = ccrs.PlateCarree() constraint_list = basic_test_func.generate_dummy_constraint_list() np.random.seed(2) @@ -276,9 +308,84 @@ def test_constraint_violation_repair(): show_depth=False, show_gcr=False ) + old_route_lc = get_route_lc(old_route[0, 0]) + new_route_lc = get_route_lc(new_route) + ax.add_collection(old_route_lc) + ax.add_collection(new_route_lc) + + cbar = fig.colorbar(old_route_lc, ax=ax, orientation='vertical', pad=0.15, shrink=0.7) + cbar.set_label('Geschwindigkeit ($m/s$)') + + pyplot.tight_layout() + plt.saveas = "test_constraint_violation_repair.png" - ax.plot(new_route[:, 1], new_route[:, 0], color="blue", transform=input_crs, marker='o') - ax.plot(old_route[0, 0][:, 1], old_route[0, 0][:, 0], color="firebrick", transform=input_crs, marker='o') assert np.array_equal(new_route[0], old_route[0, 0][0]) assert np.array_equal(new_route[-2], old_route[0, 0][-2]) assert np.array_equal(new_route[-1], old_route[0, 0][-1]) + + +def test_recalculate_speed_for_route(): + dirname = os.path.dirname(__file__) + configpath = os.path.join(dirname, 'config.isofuel_single_route.json') + config = Config.assign_config(Path(configpath)) + config.ARRIVAL_TIME = datetime(2025, 4, 2, 11, 11) + config.DEPARTURE_TIME = datetime(2025, 4, 1, 11, 11) + constraint_list = basic_test_func.generate_dummy_constraint_list() + + pop = IsoFuelPopulation( + config=config, + default_route=[35.199, 15.490, 32.737, 28.859], + constraints_list=constraint_list, + pop_size=1 + ) + rt = get_dummy_route_input() + rt = rt[0, 0] + new_route = copy.deepcopy(rt) + new_route = pop.recalculate_speed_for_route(new_route) + + dist_to_dest = 1262000 * u.meter + time_difference = config.ARRIVAL_TIME - config.DEPARTURE_TIME + bs_approx = dist_to_dest / (time_difference.total_seconds() * u.second) + + assert np.all((new_route[:, 2] - bs_approx.value) < 0.3) + + +@pytest.mark.skip(reason="Test needs modified route array.") +def test_single_point_crossover(plt): + dirname = os.path.dirname(__file__) + configpath = os.path.join(dirname, 'config.isofuel_single_route.json') + config = Config.assign_config(Path(configpath)) + default_map = Map(32., 15, 36, 29) + input_crs = ccrs.PlateCarree() + constraint_list = basic_test_func.generate_dummy_constraint_list() + departure_time = datetime(2025, 4, 1, 11, 11) + + np.random.seed(2) + + X = get_dummy_route_input() + old_route = copy.deepcopy(X) + + sp = SinglePointCrossover( + config=config, + constraints_list=constraint_list, + departure_time=departure_time + ) + # r1, r2 = sp.crossover(X[0,0], X[1,0]) + X = sp._do(problem=None, X=X) + + # plot figure with original and mutated routes + fig, ax = graphics.generate_basemap( + map=default_map.get_var_tuple(), + depth=None, + start=(35.199, 15.490), + finish=(32.737, 28.859), + title='', + show_depth=False, + show_gcr=False + ) + + ax.plot(X[0, 0][:, 1], old_route[0, 0][:, 0], color="green", transform=input_crs, marker='o') + ax.plot(old_route[0, 0][:, 1], old_route[0, 0][:, 0], color="green", transform=input_crs, marker='o') + ax.plot(old_route[1, 0][:, 1], old_route[0, 0][:, 0], color="orange", transform=input_crs, marker='o') + + plt.saveas = "test_single_point_crossoverr.png" diff --git a/tests/test_routeparams.py b/tests/test_routeparams.py index cfb0691..44b613d 100644 --- a/tests/test_routeparams.py +++ b/tests/test_routeparams.py @@ -1,6 +1,7 @@ import os from datetime import datetime, timedelta +from geovectorslib import geod import numpy as np from astropy import units as u @@ -12,7 +13,7 @@ def test_get_acc_variables(): lats = np.array([40, 50, 60, 70]) lons = np.array([4, 5, 6, 7]) - fuel_rate = np.array([1.12, 1.13, 1.15]) * u.kg/u.s + fuel_rate = np.array([1.12, 1.13, 1.15]) * u.kg / u.s dist = np.array([100, 200, 150]) * u.meter start_time = np.array([datetime(2022, 12, 19), datetime(2022, 12, 19) + timedelta(hours=1), @@ -25,7 +26,7 @@ def test_get_acc_variables(): fuel_rate=fuel_rate, power=dummy * u.Watt, rpm=dummy * u.Hz, - speed=dummy * u.m/u.s, + speed=dummy * u.m / u.s, r_calm=dummy * u.N, r_wind=dummy * u.N, r_waves=dummy * u.N, @@ -34,11 +35,11 @@ def test_get_acc_variables(): wave_height=dummy * u.m, wave_direction=dummy * u.rad, wave_period=dummy * u.second, - u_currents=dummy * u.m/u.s, - v_currents=dummy * u.m/u.s, - u_wind_speed=dummy * u.m/u.s, - v_wind_speed=dummy * u.m/u.s, - pressure=dummy * u.kg/u.meter/u.second**2, + u_currents=dummy * u.m / u.s, + v_currents=dummy * u.m / u.s, + u_wind_speed=dummy * u.m / u.s, + v_wind_speed=dummy * u.m / u.s, + pressure=dummy * u.kg / u.meter / u.second ** 2, air_temperature=dummy * u.deg_C, salinity=dummy * u.dimensionless_unscaled, water_temperature=dummy * u.deg_C, @@ -74,3 +75,46 @@ def test_get_acc_variables(): assert test_fuel == rp_test.get_full_fuel() assert np.allclose(test_dist, rp_test.get_full_dist()) assert test_time == rp_test.get_full_travel_time() + + +''' + Test whether parameters in `RouteParams.get_per_waypoint_coords` are calculated correctly. +''' + + +def test_get_waypoint_coords(): + bs = 6 * u.meter / u.second + start_time = datetime.strptime("2023-07-20T10:00Z", '%Y-%m-%dT%H:%MZ') + route_lats = np.array([54.9, 54.7, 54.5, 54.2]) + route_lons = np.array([13.2, 13.4, 13.7, 13.9]) + + start_lats_test = np.array([54.9, 54.7, 54.5]) + start_lons_test = np.array([13.2, 13.4, 13.7]) + dists_test = np.full(route_lats.shape[0] - 1, -99.) + start_times_test = np.full(route_lats.shape[0] - 1, start_time) + start_times_test[0] = start_time + travel_times_test = np.full(route_lats.shape[0] - 1, timedelta(seconds=0.)) + + for ipoint in range(3): + start_lat = route_lats[ipoint] + start_lon = route_lons[ipoint] + end_lat = route_lats[ipoint + 1] + end_lon = route_lons[ipoint + 1] + + dists_test[ipoint] = geod.inverse( + [start_lat], + [start_lon], + [end_lat], + [end_lon] + )['s12'][0] + travel_times_test[ipoint] = (dists_test[ipoint] * u.meter / bs).value + if ipoint < 2: + start_times_test[ipoint + 1] = start_times_test[ipoint] + timedelta(seconds=travel_times_test[ipoint]) + + waypoint_dict = RouteParams.get_per_waypoint_coords(route_lons, route_lats, start_time, bs) + + assert np.all(start_lats_test == waypoint_dict['start_lats']) + assert np.all(start_lons_test == waypoint_dict['start_lons']) + assert np.all(dists_test == waypoint_dict['dist'].value) + assert np.all(start_times_test == waypoint_dict['start_times']) + assert np.all(travel_times_test == waypoint_dict['travel_times'].value)