Repository: PyVRP/PyVRP
Pull request #711 (2024-12-04 20:10:30)

Make `CostEvaluator` penalties doubles
Discussed with @leonlan today. This might help the adjusting penalty mechanism function a bit more effectively.

### Benchmarks on VRPTW

This PR (c27abce16605b43d3df6d7332cf6f543d8969cbc; job ID 8936788; 10x seed / 1h):
```
               MEAN     MIN     MAX  COV
      group
obj   C1    41394.6 41331.6 41474.4  0.1
      C2    16196.0 16188.1 16207.0  0.0
      R1    47322.2 47175.1 47457.1  0.2
      R2    27789.0 27715.7 27876.2  0.2
      RC1   44302.1 44179.5 44433.8  0.2
      RC2   23277.3 23211.1 23367.4  0.2
iters C1    41153.1 38454.0 44936.0  5.3
      C2    68980.4 66037.7 73756.2  3.5
      R1    30716.4 28265.4 33251.8  5.5
      R2    46666.1 41527.8 52342.2  7.6
      RC1   31535.2 28756.7 34437.4  5.8
      RC2   48554.0 42973.0 54788.7  8.0
```

Current main (70855082b0b9e5eff849df506c9fc9e55b71e8b3; job ID 8935548; 10x seed / 1h):
```
               MEAN     MIN     MAX  COV
      group
obj   C1    41392.2 41333.2 41464.0  0.1
      C2    16203.4 16194.3 16214.9  0.0
      R1    47352.5 47210.1 47512.2  0.2
      R2    27834.9 27729.2 27936.5  0.3
      RC1   44333.5 44212.1 44473.0  0.2
      RC2   23311.6 23228.3 23396.8  0.2
iters C1    54374.0 50054.1 59133.9  5.1
      C2    81281.6 75760.5 88082.0  5.1
      R1    28758.4 26047.9 31823.8  6.1
      R2    51628.1 46260.3 58864.5  7.7
      RC1   29438.7 26725.9 32977.2  6.7
      RC2   54580.8 47919.0 63167.9  8.5
```

<details>

**Notes**:

Please read our [contributing guidelines](https://pyvrp.org/dev/contributing.html) first.
In particular:

- You must add tests when making code changes.
  This keeps the code coverage level up, and helps ensure the changes work as intended.
- When fixing a bug, you must add a test that would produce the bug in the master branch, and then show that it is fixed with the new code. 
- New code additions must be well formatted. Changes should pass the pre-commit workflow, which you can set up locally using [pre-commit](https://pre-commit.com/#intro). 
- Docstring additions must render correctly, including escapes and LaTeX.
- Finally, it is essential that all contributions in this PR are license-compatible with PyVRP's MIT license.
  Please check that this PR can be included into PyVRP under the MIT license.

</details>


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PyVRP/PyVRP/pull/711?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.98%. Comparing base [(`7085508`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/70855082b0b9e5eff849df506c9fc9e55b71e8b3?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) to head [(`277bc22`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/277bc22ea7828c23cc96eeef912722b18b4dac1f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #711      +/-   ##
==========================================
- Coverage   96.68%   94.98%   -1.70%     
==========================================
  Files          76       74       -2     
  Lines        3737     3433     -304     
  Branches      993      886     -107     
==========================================
- Hits         3613     3261     -352     
+ Misses         14       13       -1     
- Partials      110      159      +49     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PyVRP/PyVRP/pull/711?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).

codspeed-hq[bot]: ## [CodSpeed Performance Report](https://codspeed.io/PyVRP/PyVRP/branches/make-penalties-floats)

### Merging #711 will **not alter performance**

<sub>Comparing <code>make-penalties-floats</code> (277bc22) with <code>main</code> (7085508)</sub>



### Summary

`✅ 42` untouched benchmarks  




N-Wouda: I think this warrants a benchmark. I'll set that up later today.

N-Wouda: @leonlan looks pretty good right? 😄 
leonlan: > @leonlan looks pretty good right? 😄

Not bad 😃 

Patch:
diff --git a/pyvrp/PenaltyManager.py b/pyvrp/PenaltyManager.py
index 437e03b..a6123b0 100644
--- a/pyvrp/PenaltyManager.py
+++ b/pyvrp/PenaltyManager.py
@@ -116,14 +116,14 @@ class PenaltyManager:
         MAX_PENALTY]``.
     """
 
-    MIN_PENALTY = 1
-    MAX_PENALTY = 100_000
-    FEAS_TOL = 0.05
+    MIN_PENALTY: float = 0.1
+    MAX_PENALTY: float = 100_000.0
+    FEAS_TOL: float = 0.05
 
     def __init__(
         self,
         params: PenaltyParams = PenaltyParams(),
-        initial_penalties: tuple[int, int, int] = (20, 6, 6),
+        initial_penalties: tuple[float, float, float] = (20, 6, 6),
     ):
         self._params = params
         self._penalties = np.clip(
@@ -178,12 +178,12 @@ class PenaltyManager:
 
         # Initial penalty parameters are meant to weigh an average increase
         # in the relevant value by the same amount as the average edge cost.
-        init_load = round(avg_cost / max(avg_load, 1))
-        init_tw = round(avg_cost / max(avg_duration, 1))
-        init_dist = round(avg_cost / max(avg_distance, 1))
+        init_load = avg_cost / max(avg_load, 1)
+        init_tw = avg_cost / max(avg_duration, 1)
+        init_dist = avg_cost / max(avg_distance, 1)
         return cls(params, (init_load, init_tw, init_dist))
 
-    def _compute(self, penalty: int, feas_percentage: float) -> int:
+    def _compute(self, penalty: float, feas_percentage: float) -> float:
         # Computes and returns the new penalty value, given the current value
         # and the percentage of feasible solutions since the last update.
         diff = self._params.target_feasible - feas_percentage
@@ -191,13 +191,12 @@ class PenaltyManager:
         if abs(diff) < self.FEAS_TOL:
             return penalty
 
-        # +/- 1 to ensure we do not get stuck at the same integer values.
         if diff > 0:
-            new_penalty = self._params.penalty_increase * penalty + 1
+            new_penalty = self._params.penalty_increase * penalty
         else:
-            new_penalty = self._params.penalty_decrease * penalty - 1
+            new_penalty = self._params.penalty_decrease * penalty
 
-        clipped = int(np.clip(new_penalty, self.MIN_PENALTY, self.MAX_PENALTY))
+        clipped = np.clip(new_penalty, self.MIN_PENALTY, self.MAX_PENALTY)
 
         if clipped == self.MAX_PENALTY:
             msg = """
@@ -211,7 +210,7 @@ class PenaltyManager:
 
         return clipped
 
-    def _register(self, feas_list: list[bool], penalty: int, is_feas: bool):
+    def _register(self, feas_list: list[bool], penalty: float, is_feas: bool):
         feas_list.append(is_feas)
 
         if len(feas_list) != self._params.solutions_between_updates:
diff --git a/pyvrp/_pyvrp.pyi b/pyvrp/_pyvrp.pyi
index 107c1fa..b14b599 100644
--- a/pyvrp/_pyvrp.pyi
+++ b/pyvrp/_pyvrp.pyi
@@ -5,9 +5,9 @@ import numpy as np
 class CostEvaluator:
     def __init__(
         self,
-        load_penalty: int,
-        tw_penalty: int,
-        dist_penalty: int,
+        load_penalty: float,
+        tw_penalty: float,
+        dist_penalty: float,
     ) -> None: ...
     def load_penalty(self, load: int, capacity: int) -> int: ...
     def tw_penalty(self, time_warp: int) -> int: ...
diff --git a/pyvrp/cpp/CostEvaluator.cpp b/pyvrp/cpp/CostEvaluator.cpp
index 4312172..8e8a1d8 100644
--- a/pyvrp/cpp/CostEvaluator.cpp
+++ b/pyvrp/cpp/CostEvaluator.cpp
@@ -2,7 +2,9 @@
 
 using pyvrp::CostEvaluator;
 
-CostEvaluator::CostEvaluator(Cost loadPenalty, Cost twPenalty, Cost distPenalty)
+CostEvaluator::CostEvaluator(double loadPenalty,
+                             double twPenalty,
+                             double distPenalty)
     : loadPenalty_(loadPenalty),
       twPenalty_(twPenalty),
       distPenalty_(distPenalty)
diff --git a/pyvrp/cpp/CostEvaluator.h b/pyvrp/cpp/CostEvaluator.h
index 99c6828..3631187 100644
--- a/pyvrp/cpp/CostEvaluator.h
+++ b/pyvrp/cpp/CostEvaluator.h
@@ -42,7 +42,7 @@ concept DeltaCostEvaluatable = requires(T arg, size_t dimension) {
 };
 
 /**
- * CostEvaluator(load_penalty: int, tw_penalty: int, dist_penalty: int)
+ * CostEvaluator(load_penalty: float, tw_penalty: float, dist_penalty: float)
  *
  * Creates a CostEvaluator instance.
  *
@@ -61,12 +61,12 @@ concept DeltaCostEvaluatable = requires(T arg, size_t dimension) {
  */
 class CostEvaluator
 {
-    Cost loadPenalty_;
-    Cost twPenalty_;
-    Cost distPenalty_;
+    double loadPenalty_;
+    double twPenalty_;
+    double distPenalty_;
 
 public:
-    CostEvaluator(Cost loadPenalty, Cost twPenalty, Cost distPenalty);
+    CostEvaluator(double loadPenalty, double twPenalty, double distPenalty);
 
     /**
      * Computes the total excess load penalty for the given load and vehicle
@@ -170,7 +170,7 @@ public:
 Cost CostEvaluator::loadPenalty(Load load, Load capacity) const
 {
     auto const excessLoad = std::max<Load>(load - capacity, 0);
-    return static_cast<Cost>(excessLoad) * loadPenalty_;
+    return static_cast<Cost>(excessLoad.get() * loadPenalty_);
 }
 
 Cost CostEvaluator::twPenalty([[maybe_unused]] Duration timeWarp) const
@@ -178,14 +178,14 @@ Cost CostEvaluator::twPenalty([[maybe_unused]] Duration timeWarp) const
 #ifdef PYVRP_NO_TIME_WINDOWS
     return 0;
 #else
-    return static_cast<Cost>(timeWarp) * twPenalty_;
+    return static_cast<Cost>(timeWarp.get() * twPenalty_);
 #endif
 }
 
 Cost CostEvaluator::distPenalty(Distance distance, Distance maxDistance) const
 {
     auto const excessDistance = std::max<Distance>(distance - maxDistance, 0);
-    return static_cast<Cost>(excessDistance) * distPenalty_;
+    return static_cast<Cost>(excessDistance.get() * distPenalty_);
 }
 
 template <CostEvaluatable T>
diff --git a/pyvrp/cpp/bindings.cpp b/pyvrp/cpp/bindings.cpp
index d4ed8d9..ca1cf79 100644
--- a/pyvrp/cpp/bindings.cpp
+++ b/pyvrp/cpp/bindings.cpp
@@ -728,7 +728,7 @@ PYBIND11_MODULE(_pyvrp, m)
              });
 
     py::class_<CostEvaluator>(m, "CostEvaluator", DOC(pyvrp, CostEvaluator))
-        .def(py::init<pyvrp::Cost, pyvrp::Cost, pyvrp::Cost>(),
+        .def(py::init<double, double, double>(),
              py::arg("load_penalty"),
              py::arg("tw_penalty"),
              py::arg("dist_penalty"))


Test patch:
diff --git a/tests/test_PenaltyManager.py b/tests/test_PenaltyManager.py
index c0c6b53..d139a53 100644
--- a/tests/test_PenaltyManager.py
+++ b/tests/test_PenaltyManager.py
@@ -1,6 +1,12 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_equal, assert_raises, assert_warns
+from numpy.testing import (
+    assert_,
+    assert_allclose,
+    assert_equal,
+    assert_raises,
+    assert_warns,
+)
 
 from pyvrp import PenaltyManager, PenaltyParams, Solution, VehicleType
 from pyvrp.exceptions import PenaltyBoundWarning
@@ -95,22 +101,21 @@ def test_load_penalty_update_increase(ok_small):
         pm.register(sol)
     assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
 
-    # Below targetFeasible, so should increase the loadPenalty by +1 (normally
-    # to 1.1 due to penaltyIncrease, but we should not end up at the same int).
+    # Below targetFeasible, so should increase the loadPenalty to 1.1 due to
+    # penaltyIncrease, and then int((2 - 1) * 1.1) = 1.
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 2)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
 
     # Now we start from a much bigger initial loadPenalty. Here we want the
-    # penalty to increase by 10% due to penaltyIncrease = 1.1, and +1 due to
-    # double -> int.
+    # penalty to increase by 10% due to penaltyIncrease = 1.1.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params, initial_penalties=(100, 1, 1))
 
     assert_equal(pm.cost_evaluator().load_penalty(2, 1), 100)
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 111)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 110)
 
 
 def test_load_penalty_update_decrease(ok_small):
@@ -135,10 +140,10 @@ def test_load_penalty_update_decrease(ok_small):
     assert_equal(pm.cost_evaluator().load_penalty(2, 1), 4)
 
     # Above targetFeasible, so should decrease the loadPenalty to 90%, and -1
-    # from the bounds check. So 0.9 * 4 = 3.6, 3.6 - 1 = 2.6, (int) 2.6 = 2
+    # from the bounds check. So 0.9 * 4 = 3.6, and int(3.6) = 3.
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 2)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 3)
 
     # Now we start from a much bigger initial loadPenalty. Here we want the
     # penalty to decrease by 10% due to penaltyDecrease = 0.9, and -1 due to
@@ -149,16 +154,16 @@ def test_load_penalty_update_decrease(ok_small):
     assert_equal(pm.cost_evaluator().load_penalty(2, 1), 100)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 89)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 90)
 
-    # Test that the penalty cannot decrease beyond 1, its minimum value.
+    # Test that the penalty cannot decrease beyond MIN_PENALTY.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 1, 1))
+    pm = PenaltyManager(params, initial_penalties=(0.1, 1, 1))
 
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(11, 1), 1)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 1)
+    assert_equal(pm.cost_evaluator().load_penalty(11, 1), 1)
 
 
 def test_time_warp_penalty_update_increase(ok_small):
@@ -182,11 +187,11 @@ def test_time_warp_penalty_update_increase(ok_small):
         pm.register(sol)
     assert_equal(pm.cost_evaluator().tw_penalty(1), 1)
 
-    # Below targetFeasible, so should increase the tw penalty by +1 (normally
-    # to 1.1 due to penaltyIncrease, but we should not end up at the same int).
+    # Below targetFeasible, so should increase the tw penalty to 1.1 due to
+    # penaltyIncrease, and int(1 * 1.1) = 1.
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 2)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), 1)
 
     # Now we start from a much bigger initial tw penalty. Here we want the
     # penalty to increase by 10% due to penaltyIncrease = 1.1, and +1 due
@@ -197,7 +202,7 @@ def test_time_warp_penalty_update_increase(ok_small):
     assert_equal(pm.cost_evaluator().tw_penalty(1), 100)
     for sol in [infeas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 111)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), 110)
 
 
 def test_time_warp_penalty_update_decrease(ok_small):
@@ -222,11 +227,11 @@ def test_time_warp_penalty_update_decrease(ok_small):
         pm.register(sol)
     assert_equal(pm.cost_evaluator().tw_penalty(1), 4)
 
-    # Above targetFeasible, so should decrease the twCapacity to 90%, and -1
-    # from the bounds check. So 0.9 * 4 = 3.6, 3.6 - 1 = 2.6, (int) 2.6 = 2.
+    # Above targetFeasible, so should decrease the twCapacity to 90%. So
+    # 0.9 * 4 = 3.6, and int(3.6) = 3.
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 2)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), 3)
 
     # Now we start from a much bigger initial twCapacity. Here we want the
     # penalty to decrease by 10% due to penaltyDecrease = 0.9, and -1 due
@@ -237,16 +242,16 @@ def test_time_warp_penalty_update_decrease(ok_small):
     assert_equal(pm.cost_evaluator().tw_penalty(1), 100)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 89)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), 90)
 
-    # Test that the penalty cannot decrease beyond 1, its minimum value.
+    # Test that the penalty cannot decrease beyond MIN_PENALTY.
     params = PenaltyParams(1, num_registrations, 1.1, 0.9, 0.5)
-    pm = PenaltyManager(params, initial_penalties=(1, 1, 1))
+    pm = PenaltyManager(params, initial_penalties=(1, 0.1, 1))
 
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 1)
+    assert_equal(pm.cost_evaluator().tw_penalty(10), 1)
     for sol in [feas] * num_registrations:
         pm.register(sol)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 1)
+    assert_equal(pm.cost_evaluator().tw_penalty(10), 1)
 
 
 def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
@@ -283,9 +288,9 @@ def test_does_not_update_penalties_before_sufficient_registrations(ok_small):
     # Register a fourth time. Now the penalties should change. Since there are
     # more feasible registrations than desired, the penalties should decrease.
     pm.register(feas)
-    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 2)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), 2)
-    assert_equal(pm.cost_evaluator().dist_penalty(1, 0), 2)
+    assert_equal(pm.cost_evaluator().load_penalty(2, 1), 3)
+    assert_equal(pm.cost_evaluator().tw_penalty(1), 3)
+    assert_equal(pm.cost_evaluator().dist_penalty(1, 0), 3)
 
 
 @pytest.mark.filterwarnings("ignore::pyvrp.exceptions.PenaltyBoundWarning")
@@ -321,7 +326,9 @@ def test_max_min_penalty(ok_small):
     # should drop to zero due to the penalty_decrease parameter. But penalty
     # parameters cannot drop below MIN_PENALTY.
     pm.register(feas)
-    assert_equal(pm.cost_evaluator().tw_penalty(1), PenaltyManager.MIN_PENALTY)
+    assert_allclose(
+        pm.cost_evaluator().tw_penalty(10), 10 * PenaltyManager.MIN_PENALTY
+    )
 
 
 def test_warns_max_penalty_value(ok_small):
@@ -350,7 +357,7 @@ def test_init_from_load_penalty_value(ok_small):
 
     avg_cost = ok_small.distance_matrix(0).mean()
     avg_load = np.mean([c.delivery for c in ok_small.clients()])
-    assert_equal(cost_eval.load_penalty(1, 0), round(avg_cost / avg_load))
+    assert_equal(cost_eval.load_penalty(1, 0), int(avg_cost / avg_load))
 
 
 def test_init_from_tw_penalty_value(ok_small):
@@ -399,6 +406,9 @@ def test_init_clips_penalties():
     pm = PenaltyManager(initial_penalties=penalties)
 
     cost_eval = pm.cost_evaluator()
-    assert_equal(cost_eval.load_penalty(1, 0), 1)  # set to MIN
+    assert_allclose(  # MIN
+        cost_eval.load_penalty(10, 0),
+        10 * PenaltyManager.MIN_PENALTY,
+    )
     assert_equal(cost_eval.tw_penalty(1), PenaltyManager.MAX_PENALTY)  # MAX
     assert_equal(cost_eval.dist_penalty(1, 0), 2)  # already OK, so unchanged