Add speed information to all chromosomes of the genetic algorithm

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

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

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