feat(swopp3): add CMA-ES diagnostics hook and stabilize FMS travel-time routing

Author: daniprec

.github/copilot-instructions.md

@@ -40,6 +40,6 @@ Write and reason in English.
 - Make sure you have the necessary permissions to push to the repository. If you do not have permissions, stop and ask for them, guiding the user to the appropriate process to gain access.
 - You can add, commit and push changes to this repository. Never commit to 'main' or 'swopp' branches directly.
 - If you are on 'main' or 'swopp', create a new branch for your changes and open a pull request for review.
-- Create tests before implementing new features or fixing bugs. Tests should be in the `tests/` directory and follow existing patterns.
+- Create tests before implementing new features or fixing bugs in `routetools` (you can ignore other directories such as `scripts` or `docs`). Tests should be in the `tests/` directory and follow existing patterns.
 - Make sure to run all tests and hooks before pushing your changes. If you encounter any issues, please fix them before pushing.
 - Do small commits, preferably one per logical change. This makes it easier to review and understand the history of changes.

docs/swopp3_sweep_results.md

@@ -1,8 +1,8 @@
 # SWOPP3 Parameter Sweep Results
 
-**Date:** 2026-03-19  
-**Branch:** `feat/swopp3-final-results`  
-**Route:** Pacific noWPS (no waypoint system)  
+**Date:** 2026-03-19
+**Branch:** `feat/swopp3-final-results`
+**Route:** Pacific noWPS (no waypoint system)
 **Operational constraints:** TWS ≤ 20 m/s, Hs ≤ 7.0 m
 
 ---
@@ -19,10 +19,10 @@ where $x_i$ is TWS or Hs at segment midpoint $i$.
 
 Two parameters control the penalty magnitude:
 
-| Parameter | Symbol | Role |
-|-----------|--------|------|
-| `weather_penalty_weight` (wpw) | outer multiplier | **Global scale** — multiplies the entire penalty term before it is added to the energy cost. Controls how much the penalty matters *relative to fuel cost*. |
-| `sharpness` | inner multiplier | **Violation sensitivity** — multiplies each individual squared excess *before* summation. Controls how steeply the penalty ramps up *per segment* as conditions worsen beyond the threshold. |
+| Parameter                      | Symbol           | Role                                                                                                                                                                                         |
+| ------------------------------ | ---------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `weather_penalty_weight` (wpw) | outer multiplier | **Global scale** — multiplies the entire penalty term before it is added to the energy cost. Controls how much the penalty matters _relative to fuel cost_.                                  |
+| `sharpness`                    | inner multiplier | **Violation sensitivity** — multiplies each individual squared excess _before_ summation. Controls how steeply the penalty ramps up _per segment_ as conditions worsen beyond the threshold. |
 
 They are mathematically interchangeable in a single-field scenario
 ($\texttt{wpw} \cdot \texttt{sharpness}$ acts as one effective multiplier),
@@ -31,6 +31,7 @@ scales each field's excess identically (same ramp steepness for TWS and Hs),
 while `wpw` scales the combined total against the fuel cost.
 
 **In practice for this sweep:**
+
 - Higher **wpw** → optimizer sees weather avoidance as more important than fuel savings → routes detour more.
 - Higher **sharpness** → small exceedances above the limit are punished more aggressively → optimizer reacts earlier to marginal violations.
 
@@ -42,38 +43,39 @@ while `wpw` scales the combined total against the fuel cost.
 over-spend (+30 to +99%) seen on stormy Pacific crossings.
 
 **Setup:**
+
 - 10 departures with worst delta%: [24, 337, 338, 339, 340, 341, 348, 349, 362, 364]
 - Grid: `wpw ∈ {5, 10, 20, 40, 100}` × `sharpness ∈ {1, 2, 5}` × `σ₀ ∈ {0.1, 0.2}`
 - 30 configs × 10 departures = **300 runs** completed in 622 s (~2 s/run)
 
 ### Summary by Configuration
 
-| wpw | sharp | σ₀ | mean Δ% | #viol TWS | #viol Hs | #any viol | mean TWS | mean Hs |
-|----:|------:|---:|--------:|----------:|---------:|----------:|---------:|--------:|
-|   5 |     1 | 0.2 | **-5.0** |     5 |    10 |    10 |     19.5 |    7.92 |
-|   5 |     2 | 0.2 |   -3.0 |     5 |     9 |     9 |     19.4 |    7.72 |
-|  10 |     1 | 0.2 |   -3.0 |     3 |     9 |     9 |     19.4 |    7.72 |
-|   5 |     1 | 0.1 |   -2.6 |     5 |     9 |     9 |     19.9 |    8.11 |
-|  10 |     1 | 0.1 |    0.0 |     5 |     8 |     9 |     19.7 |    7.66 |
-|   5 |     2 | 0.1 |    0.3 |     4 |     8 |     8 |     19.7 |    7.52 |
-| ... | ... | ... | ... | ... | ... | ... | ... | ... |
-| 100 |     5 | 0.1 | **+24.6** |     1 |     7 |     7 |     19.0 |    7.22 |
-| 100 |     2 | 0.1 |  +25.4 |     2 |     7 |     7 |     18.7 |    7.24 |
+| wpw | sharp |  σ₀ |   mean Δ% | #viol TWS | #viol Hs | #any viol | mean TWS | mean Hs |
+| --: | ----: | --: | --------: | --------: | -------: | --------: | -------: | ------: |
+|   5 |     1 | 0.2 |  **-5.0** |         5 |       10 |        10 |     19.5 |    7.92 |
+|   5 |     2 | 0.2 |      -3.0 |         5 |        9 |         9 |     19.4 |    7.72 |
+|  10 |     1 | 0.2 |      -3.0 |         3 |        9 |         9 |     19.4 |    7.72 |
+|   5 |     1 | 0.1 |      -2.6 |         5 |        9 |         9 |     19.9 |    8.11 |
+|  10 |     1 | 0.1 |       0.0 |         5 |        8 |         9 |     19.7 |    7.66 |
+|   5 |     2 | 0.1 |       0.3 |         4 |        8 |         8 |     19.7 |    7.52 |
+| ... |   ... | ... |       ... |       ... |      ... |       ... |      ... |     ... |
+| 100 |     5 | 0.1 | **+24.6** |         1 |        7 |         7 |     19.0 |    7.22 |
+| 100 |     2 | 0.1 |     +25.4 |         2 |        7 |         7 |     18.7 |    7.24 |
 
 ### Best Config vs Baseline — Per Departure
 
 | Dep | Baseline (wpw=100,s=5,σ=0.1) Δ% | Baseline viol | Best (wpw=5,s=1,σ=0.2) Δ% | Best viol |
-|----:|---:|---|---:|---|
-|  24 |  +28.2 | Hs | -21.6 | Hs |
-| 337 |  +18.2 | TWS+Hs | -5.7 | TWS+Hs |
-| 338 |  +34.9 | Hs |  +4.6 | Hs |
-| 339 |  +35.2 | Hs | +12.0 | TWS+Hs |
-| 340 |  +98.7 | Hs | +10.6 | TWS+Hs |
-| 341 |  +66.2 | Hs | +16.9 | Hs |
-| 348 |  +18.1 | **none** |  +7.0 | Hs |
-| 349 |   +8.2 | Hs |  +5.5 | Hs |
-| 362 |  -30.1 | none | -38.7 | TWS+Hs |
-| 364 |  -31.5 | none | -40.2 | TWS+Hs |
+| --: | ------------------------------: | ------------- | ------------------------: | --------- |
+|  24 |                           +28.2 | Hs            |                     -21.6 | Hs        |
+| 337 |                           +18.2 | TWS+Hs        |                      -5.7 | TWS+Hs    |
+| 338 |                           +34.9 | Hs            |                      +4.6 | Hs        |
+| 339 |                           +35.2 | Hs            |                     +12.0 | TWS+Hs    |
+| 340 |                           +98.7 | Hs            |                     +10.6 | TWS+Hs    |
+| 341 |                           +66.2 | Hs            |                     +16.9 | Hs        |
+| 348 |                           +18.1 | **none**      |                      +7.0 | Hs        |
+| 349 |                            +8.2 | Hs            |                      +5.5 | Hs        |
+| 362 |                           -30.1 | none          |                     -38.7 | TWS+Hs    |
+| 364 |                           -31.5 | none          |                     -40.2 | TWS+Hs    |
 
 ### Key Observations
 
@@ -106,25 +108,26 @@ over-spend (+30 to +99%) seen on stormy Pacific crossings.
 departures currently converging exactly to the great circle (0% delta).
 
 **Setup:**
+
 - 20 sampled zero-delta departures: [44, 45, 46, 111, 113, 116, 118, 125,
   143, 159, 160, 169, 174, 177, 185, 186, 232, 233, 239, 274]
 - Grid: `σ₀ ∈ {0.1, 0.3, 0.5}` × `popsize ∈ {200, 400}` × `maxfevals ∈ {25k, 50k}` × `K ∈ {10, 15}`
 - 24 configs × 20 departures = **480 runs** completed in 1294 s (~2.7 s/run)
 
 ### Summary by Configuration
 
-| σ₀ | pop | maxfevals | K | mean Δ% | median Δ% | min Δ% | #any viol |
-|---:|----:|----------:|--:|--------:|----------:|-------:|----------:|
-| 0.3 | 200 | 50000 | 15 | **-0.7** | 0.0 | -7.6 | 0 |
-| 0.5 | 200 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
-| 0.3 | 400 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
-| 0.5 | 400 | 50000 | 15 | -0.7 | 0.0 | -7.4 | 0 |
-| 0.3 | 200 | 25000 | 15 | -0.7 | 0.0 | -7.5 | 0 |
-| 0.5 | 200 | 25000 | 15 | -0.7 | 0.0 | -7.3 | 0 |
-| 0.1 | 200 | 25000 | 15 | -0.4 | 0.0 | -7.6 | 0 |
-| ... | ... | ... | ... | ... | ... | ... | ... |
-| 0.1 | 200 | 25000 | 10 | **0.0** | 0.0 | 0.0 | 0 |
-| 0.1 | 400 | 25000 | 10 | 0.0 | 0.0 | 0.0 | 0 |
+|  σ₀ | pop | maxfevals |   K |  mean Δ% | median Δ% | min Δ% | #any viol |
+| --: | --: | --------: | --: | -------: | --------: | -----: | --------: |
+| 0.3 | 200 |     50000 |  15 | **-0.7** |       0.0 |   -7.6 |         0 |
+| 0.5 | 200 |     50000 |  15 |     -0.7 |       0.0 |   -7.4 |         0 |
+| 0.3 | 400 |     50000 |  15 |     -0.7 |       0.0 |   -7.4 |         0 |
+| 0.5 | 400 |     50000 |  15 |     -0.7 |       0.0 |   -7.4 |         0 |
+| 0.3 | 200 |     25000 |  15 |     -0.7 |       0.0 |   -7.5 |         0 |
+| 0.5 | 200 |     25000 |  15 |     -0.7 |       0.0 |   -7.3 |         0 |
+| 0.1 | 200 |     25000 |  15 |     -0.4 |       0.0 |   -7.6 |         0 |
+| ... | ... |       ... | ... |      ... |       ... |    ... |       ... |
+| 0.1 | 200 |     25000 |  10 |  **0.0** |       0.0 |    0.0 |         0 |
+| 0.1 | 400 |     25000 |  10 |      0.0 |       0.0 |    0.0 |         0 |
 
 ### Key Observations
 
@@ -152,14 +155,14 @@ departures currently converging exactly to the great circle (0% delta).
 
 ### Parameter Selection
 
-| Parameter | Recommended | Rationale |
-|-----------|-------------|-----------|
-| **K** | **15** | Consistently improves zero-delta departures; no downside |
-| **σ₀** | **0.3** | Best tradeoff: helps Stage A detours and Stage B exploration |
-| **popsize** | **200** | 400 shows no benefit; saves compute |
-| **maxfevals** | **25000** | 50k shows no benefit on these departures |
-| **wpw** | **10** | Compromise: mean Δ% ≈ 0% (vs +24.6% baseline), TWS violations drop from 5 to 3 |
-| **sharpness** | **1** | Lower sharpness reduces over-reaction to marginal exceedances |
+| Parameter     | Recommended | Rationale                                                                      |
+| ------------- | ----------- | ------------------------------------------------------------------------------ |
+| **K**         | **15**      | Consistently improves zero-delta departures; no downside                       |
+| **σ₀**        | **0.3**     | Best tradeoff: helps Stage A detours and Stage B exploration                   |
+| **popsize**   | **200**     | 400 shows no benefit; saves compute                                            |
+| **maxfevals** | **25000**   | 50k shows no benefit on these departures                                       |
+| **wpw**       | **10**      | Compromise: mean Δ% ≈ 0% (vs +24.6% baseline), TWS violations drop from 5 to 3 |
+| **sharpness** | **1**       | Lower sharpness reduces over-reaction to marginal exceedances                  |
 
 ### Tradeoff: Energy vs Compliance
 

routetools/cmaes.py

@@ -1,3 +1,26 @@
+"""CMA-ES route optimization utilities.
+
+This module parameterizes vessel routes as Bezier curves and optimizes their
+free control points with `pycma`. The main `optimize` entry point converts an
+initial straight line or seed curve into control-point space, evaluates each
+population as full routes, applies optional land and weather penalties, and
+returns the best route in waypoint form together with summary metadata.
+
+The default path minimizes `routetools.cost.cost_function`, but callers may
+inject a custom batched `cost_fn` to align CMA-ES with another objective, such
+as SWOPP3's RISE energy model. Land intersection penalties, distance-to-land
+repulsion, and weather penalties are layered outside the custom objective so
+the population loop can be reused across synthetic and real-weather scenarios.
+
+New in this workflow is the optional `snapshot_callback` hook exposed on both
+`_cma_evolution_strategy` and `optimize`. When provided, it receives one
+dictionary per CMA-ES iteration containing the whole candidate population,
+their costs, the generation-best route, and the global-best route so far. This
+is intended for controlled visualizations and diagnostics, including the
+single-departure SWOPP3 GIF script, without changing the optimizer's default
+behavior or return type.
+"""
+
 import time
 import warnings
 from collections.abc import Callable
@@ -19,6 +42,9 @@
 from routetools.weather import weather_penalty as _weather_penalty
 from routetools.weather import weather_penalty_smooth as _weather_penalty_smooth
 
+type CmaesSnapshot = dict[str, Any]
+type CmaesSnapshotCallback = Callable[[CmaesSnapshot], None]
+
 
 @jit  # type: ignore[misc]
 def batch_bezier(t: jnp.ndarray, control: jnp.ndarray) -> jnp.ndarray:
@@ -316,6 +342,7 @@ def _cma_evolution_strategy(
     land_margin: int = 0,
     verbose: bool = True,
     bounds: list[list[float]] | None = None,
+    snapshot_callback: CmaesSnapshotCallback | None = None,
     **kwargs: dict[str, Any],
 ) -> cma.CMAEvolutionStrategy:
     curve: jnp.ndarray
@@ -429,10 +456,36 @@ def _cma_evolution_strategy(
                     cost = cost.at[idx_worst[i]].set(top_costs[i])
                     curve = curve.at[idx_worst[i], ...].set(top_curves[i, ...])
 
+        generation_best_index = int(jnp.argmin(cost))
+        generation_best_curve = curve[generation_best_index]
+        generation_best_cost = float(cost[generation_best_index])
+
         es.tell(X, cost.tolist())  # update the optimizer
         if verbose:
             es.disp()
 
+        if snapshot_callback is not None:
+            best_curve = control_to_curve(
+                jnp.asarray(es.best.x),
+                src,
+                dst,
+                L=L,
+                num_pieces=num_pieces,
+                force_L_multiple_of_num_pieces=force_L_multiple_of_num_pieces,
+            )
+            snapshot_callback(
+                {
+                    "iteration": es.countiter,
+                    "population_curves": curve,
+                    "population_costs": cost,
+                    "generation_best_index": generation_best_index,
+                    "generation_best_curve": generation_best_curve,
+                    "generation_best_cost": generation_best_cost,
+                    "best_curve": best_curve,
+                    "best_cost": float(es.best.f),
+                }
+            )
+
         # Save the top solutions (use cost as fitness)
         if keep_top > 0:
             idx_sorted = jnp.argsort(cost)
@@ -487,6 +540,7 @@ def optimize(
     land_margin: int = 0,
     verbose: bool = True,
     bounds: list[list[float]] | None = None,
+    snapshot_callback: CmaesSnapshotCallback | None = None,
 ) -> tuple[jnp.ndarray, dict[str, Any]]:
     """
     Solve the vessel routing problem for a given vector field.
@@ -570,6 +624,10 @@ def optimize(
         Per-dimension ``[lower, upper]`` bounds for CMA-ES control-point
         parameters.  Each list has length ``2*(K-2)``.  ``None`` disables
         bounds (default).
+    snapshot_callback : callable, optional
+        Called after every CMA-ES iteration with the current population,
+        per-candidate costs, the generation-best route, and the global-best
+        route so far. Default is ``None``.
 
     Returns
     -------
@@ -662,6 +720,7 @@ def optimize(
         land_margin=land_margin,
         verbose=verbose,
         bounds=bounds,
+        snapshot_callback=snapshot_callback,
     )
     time_end = time.time()
     if verbose: