Repository: PyVRP/PyVRP (python)

Issue:
Evaluation of distance/duration of empty route in local search operators
#### Describe the bug

In `Solution.cpp`, `Route` does not compute any statistics if `visits_.empty()`, hence the distance and duration will be 0. However, when moves that result in an empty route are evaluated in the local search, the distance (at least in `Exchange`) is computed as `distanceMatrix(depot, depot)`. If this value is not 0, it results in a discrepancy in the cost evaluation in local search vs the actual cost computed by the solution. This happens for example in `test_model_and_solve` in `test_Model.py`, when a relocate move to an empty route is evaluated.

#### Code sample

```bash
pytest pyvrp/tests/test_Model.py -k test_model_and_solve
```

When running this test, the following line will evaluate to a large value as the distance depot -> depot is `MAX_VALUE` as this is the default in `Model.py`. https://github.com/PyVRP/PyVRP/blob/94b3954aa2072a1ef259cafe0d3ad9eec0bdc333/pyvrp/cpp/search/Exchange.h#L93-L94 
This results in a large negative delta cost here which is incorrect:
https://github.com/PyVRP/PyVRP/blob/94b3954aa2072a1ef259cafe0d3ad9eec0bdc333/pyvrp/cpp/search/Exchange.h#L102

<details>

**Note**: Please be sure you are using the latest released version of `pyvrp`, or a recent build of `main`.
If your problem has been fixed in an unreleased version, you might be able to use `main` until a new release occurs.

**Note**: If you are using a released version, have you verified that the bug exists in the master branch of this repository?

</details>

If the issue has not been resolved, please file it in the issue tracker.
Make sure to include the version of Python you are using, and the version of `pyvrp` you run.
You can copy-paste the output of `import pyvrp; pyvrp.show_versions()` to achieve this.

#### Expected Output

I expect operators to evaluate moves and define the cost of an empty route consistently with the definition of cost in a `Solution`. Most reasonable, I'd say this means that we should enforce `dist(depot,depot) = 0`. I suggest to enforce this for all locations just in case.

Pull request description:
Enforce zero matrix diagonal
This PR:

- [x] Closes #428.
- [x] Adds and passes relevant tests.
- [x] Has well-formatted code and documentation.

(pushed in new PR because I do not have write access to the ORTEC fork)

<details>

**Notes**:

- It is essential that you add a test 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.

</details>


Patch:
diff --git a/pyvrp/Model.py b/pyvrp/Model.py
index aead388..6e829c8 100644
--- a/pyvrp/Model.py
+++ b/pyvrp/Model.py
@@ -171,11 +171,15 @@ class Model:
         Raises
         ------
         ValueError
-            When either distance or duration is a negative value.
+            When either distance or duration is a negative value, or when self
+            loops have nonzero distance or duration values.
         """
         if distance < 0 or duration < 0:
             raise ValueError("Cannot have negative edge distance or duration.")
 
+        if frm == to and (distance != 0 or duration != 0):
+            raise ValueError("A self loop must have 0 distance and duration.")
+
         if max(distance, duration) > MAX_USER_VALUE:
             msg = """
             The given distance or duration value is very large. This may impact
@@ -246,6 +250,8 @@ class Model:
         # not set below are never traversed.
         distances = np.full((len(locs), len(locs)), MAX_VALUE, dtype=int)
         durations = np.full((len(locs), len(locs)), MAX_VALUE, dtype=int)
+        np.fill_diagonal(distances, 0)
+        np.fill_diagonal(durations, 0)
 
         for edge in self._edges:
             frm = loc2idx[id(edge.frm)]
diff --git a/pyvrp/cpp/ProblemData.cpp b/pyvrp/cpp/ProblemData.cpp
index 0dca380..ec727c8 100644
--- a/pyvrp/cpp/ProblemData.cpp
+++ b/pyvrp/cpp/ProblemData.cpp
@@ -233,6 +233,15 @@ ProblemData::ProblemData(std::vector<Client> const &clients,
         throw std::invalid_argument("Duration matrix shape does not match the "
                                     "problem size.");
 
+    for (size_t idx = 0; idx != numLocations(); ++idx)
+    {
+        if (dist_(idx, idx) != 0)
+            throw std::invalid_argument("Distance matrix diagonal must be 0.");
+
+        if (dur_(idx, idx) != 0)
+            throw std::invalid_argument("Duration matrix diagonal must be 0.");
+    }
+
     for (auto const &depot : depots_)
     {
         if (depot.demand != 0)


Test patch:
diff --git a/pyvrp/search/tests/test_LocalSearch.py b/pyvrp/search/tests/test_LocalSearch.py
index 32d32be..d31e2fc 100644
--- a/pyvrp/search/tests/test_LocalSearch.py
+++ b/pyvrp/search/tests/test_LocalSearch.py
@@ -430,7 +430,7 @@ def test_local_search_does_not_remove_required_clients():
         ],
         depots=[Client(x=0, y=0)],
         vehicle_types=[VehicleType(1, capacity=50)],
-        distance_matrix=np.full((3, 3), fill_value=10),
+        distance_matrix=np.where(np.eye(3), 0, 10),
         duration_matrix=np.zeros((3, 3), dtype=int),
     )
 
diff --git a/pyvrp/tests/test_Model.py b/pyvrp/tests/test_Model.py
index e9932a5..ce36230 100644
--- a/pyvrp/tests/test_Model.py
+++ b/pyvrp/tests/test_Model.py
@@ -50,6 +50,31 @@ def test_add_edge_raises_negative_distance_or_duration():
         model.add_edge(client, depot, distance=0, duration=-1)
 
 
+def test_add_edge_raises_self_connection_nonzero_distance_or_duration():
+    """
+    Edges connecting a node to itself must have 0 distance and duration.
+    """
+    model = Model()
+    depot = model.add_depot(0, 0)
+    client = model.add_client(0, 1)
+
+    model.add_edge(depot, client, distance=1, duration=1)  # should be OK
+    model.add_edge(depot, depot, distance=0, duration=0)  # should be OK
+    model.add_edge(client, client, distance=0, duration=0)  # should be OK
+
+    with assert_raises(ValueError):  # self loop with nonzero distance not OK
+        model.add_edge(depot, depot, distance=1, duration=0)
+
+    with assert_raises(ValueError):  # self loop with nonzero duration not OK
+        model.add_edge(depot, depot, distance=0, duration=1)
+
+    with assert_raises(ValueError):  # self loop with nonzero distance not OK
+        model.add_edge(client, client, distance=1, duration=0)
+
+    with assert_raises(ValueError):  # self loop with nonzero duration not OK
+        model.add_edge(client, client, distance=0, duration=1)
+
+
 def test_add_client_attributes():
     """
     Smoke test that checks, for a single client, that the model adds a client
@@ -382,7 +407,8 @@ def test_model_solves_small_instance_with_fixed_costs():
 
     for frm in m.locations:
         for to in m.locations:
-            m.add_edge(frm, to, distance=0, duration=5)
+            if frm != to:
+                m.add_edge(frm, to, distance=0, duration=5)
 
     res = m.solve(stop=MaxIterations(100))
     assert_(res.is_feasible())
@@ -424,7 +450,8 @@ def test_model_solves_small_instance_with_shift_durations():
 
     for frm in m.locations:
         for to in m.locations:
-            m.add_edge(frm, to, distance=0, duration=5)
+            if frm != to:
+                m.add_edge(frm, to, distance=0, duration=5)
 
     res = m.solve(stop=MaxIterations(100))
     assert_(res.is_feasible())
diff --git a/pyvrp/tests/test_ProblemData.py b/pyvrp/tests/test_ProblemData.py
index 2751aee..a536a80 100644
--- a/pyvrp/tests/test_ProblemData.py
+++ b/pyvrp/tests/test_ProblemData.py
@@ -191,6 +191,27 @@ def test_problem_data_raises_when_incorrect_matrix_dimensions(matrix):
         ProblemData(clients, depots, vehicle_types, other_matrix, matrix)
 
 
+@pytest.mark.parametrize(
+    ("dist_mat", "dur_mat"),
+    [
+        (np.eye(2, dtype=int), np.zeros((2, 2), dtype=int)),  # distance diag
+        (np.zeros((2, 2), dtype=int), np.eye(2, dtype=int)),  # duration diag
+        (np.eye(2, dtype=int), np.eye(2, dtype=int)),  # both diags nonzero
+    ],
+)
+def test_problem_data_raises_matrix_diagonal_nonzero(dist_mat, dur_mat):
+    """
+    Tests that the ``ProblemData`` constructor raises a ``ValueError`` when
+    the distance or duration matrix has a non-zero value on the diagonal.
+    """
+    clients = [Client(x=0, y=0)]
+    depots = [Client(x=0, y=0)]
+    vehicle_types = [VehicleType(2, capacity=1)]
+
+    with assert_raises(ValueError):
+        ProblemData(clients, depots, vehicle_types, dist_mat, dur_mat)
+
+
 def test_problem_data_replace_no_changes():
     """
     Tests that when using ``ProblemData.replace()`` without any arguments
@@ -246,7 +267,7 @@ def test_problem_data_replace_with_changes():
     new = original.replace(
         clients=[Client(x=1, y=1)],
         vehicle_types=[VehicleType(3, 4), VehicleType(5, 6)],
-        distance_matrix=np.ones((2, 2), dtype=int),
+        distance_matrix=np.where(np.eye(2), 0, 2),
     )
 
     assert_(new is not original)
@@ -288,16 +309,16 @@ def test_problem_data_replace_raises_mismatched_argument_shapes():
         data.replace(clients=[])  # matrices are 2x2
 
     with assert_raises(ValueError):
-        data.replace(distance_matrix=np.ones((3, 3), dtype=int))  # two clients
+        data.replace(distance_matrix=np.where(np.eye(3), 0, 1))  # two clients
 
     with assert_raises(ValueError):
-        data.replace(duration_matrix=np.ones((3, 3), dtype=int))  # two clients
+        data.replace(duration_matrix=np.where(np.eye(3), 0, 1))  # two clients
 
     with assert_raises(ValueError):
         data.replace(
             clients=[Client(x=1, y=1)],
-            distance_matrix=np.ones((3, 3), dtype=int),
-            duration_matrix=np.ones((3, 3), dtype=int),
+            distance_matrix=np.where(np.eye(3), 0, 1),
+            duration_matrix=np.where(np.eye(3), 0, 1),
         )
 
 
@@ -323,6 +344,9 @@ def test_matrix_access():
 
     dist_mat = gen.integers(500, size=(size, size))
     dur_mat = gen.integers(500, size=(size, size))
+    np.fill_diagonal(dist_mat, 0)
+    np.fill_diagonal(dur_mat, 0)
+
     clients = [
         Client(x=0, y=0, demand=0, service_duration=0, tw_early=0, tw_late=10)
         for _ in range(size)
diff --git a/pyvrp/tests/test_Solution.py b/pyvrp/tests/test_Solution.py
index 5b2d499..a2d27a1 100644
--- a/pyvrp/tests/test_Solution.py
+++ b/pyvrp/tests/test_Solution.py
@@ -570,9 +570,9 @@ def test_route_release_time():
 @mark.parametrize(
     "dist_mat",
     [
-        np.full((3, 3), fill_value=100, dtype=int),
-        np.full((3, 3), fill_value=1, dtype=int),
-        np.full((3, 3), fill_value=1000, dtype=int),
+        np.where(np.eye(3), 0, 100),
+        np.where(np.eye(3), 0, 1),
+        np.where(np.eye(3), 0, 1000),
     ],
 )
 def test_time_warp_for_a_very_constrained_problem(dist_mat):