Repository: PyVRP/PyVRP (python)

Issue:
Order neighbours by proximity?
In LocalSearch, we shuffle the order of the nodes, but not the order of the neighbours since arguably this does not change much if we already shuffle the node. However, since the neighbours are sorted by index this makes lower index nodes much more often checked.

I suggest to try the suggestion by Team Kleopatra (see #181) to sort the neighbours by proximity, which they claim yields faster convergence as the best moves are checked and thus applied first.

Pull request description:
Sort granular neighbourhood by proximity
This PR:

- [x] Closes #201, part of #181.
- [x] Adds and passes relevant tests.
- [x] Has well-formatted code and documentation.

### Benchmarks

Current main (45e30f9bbc0cf6adaa1e7d5c95b0a9b43157664a, job ID 959880)
```
 Avg. objective: 333444.8 (min: 333296, max: 333585)
Avg. iterations: 102187.1 (min: 75036, max: 111424)
```
This branch (cdaf1f70c35947d1392aae326ec3cccd2fe3480f, job ID 960271)
```
 Avg. objective: 333434.3 (min: 333352, max: 333642)
Avg. iterations: 102605.7 (min: 83489, max: 108594)
```
This branch (30f468584e7bec313576fd33df18efb035674fcb, job ID 1041857)
```
 Avg. objective: 333556.2 (min: 333387, max: 333768)
Avg. iterations:  88005.9 (min: 73832, max: 104752)
```

<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/_ProblemData.pyi b/pyvrp/_ProblemData.pyi
index 4a4692b..a1ede86 100644
--- a/pyvrp/_ProblemData.pyi
+++ b/pyvrp/_ProblemData.pyi
@@ -1,36 +1,48 @@
-from typing import List, Tuple
+from typing import List
 
 from pyvrp._Matrix import Matrix
 
 class Client:
     """
-    Simple data object storing all client data as properties.
+    Simple data object storing all client data as (read-only) properties.
 
-    Attributes
+    Parameters
     ----------
-    demand
-        The amount this client's demanding.
-    service_duration
-        This client's service duration, that is, the amount of time we need to
-        visit the client for.
-    tw_early
-        Earliest time at which we can visit this client.
-    tw_late
-        Latest time at which we can visit this client.
     x
         Horizontal coordinate of this client, that is, the 'x' part of the
         client's (x, y) location tuple.
     y
         Vertical coordinate of this client, that is, the 'y' part of the
         client's (x, y) location tuple.
+    demand
+        The amount this client's demanding. Default 0.
+    service_duration
+        This client's service duration, that is, the amount of time we need to
+        visit the client for. Service should start (but not necessarily end)
+        within the [:py:attr:`~tw_early`, :py:attr:`~tw_late`] interval.
+        Default 0.
+    tw_early
+        Earliest time at which we can visit this client. Default 0.
+    tw_late
+        Latest time at which we can visit this client. Default 0.
     """
 
+    x: int
+    y: int
     demand: int
     service_duration: int
     tw_early: int
     tw_late: int
-    x: int
-    y: int
+
+    def __init__(
+        self,
+        x: int,
+        y: int,
+        demand: int = 0,
+        service_duration: int = 0,
+        tw_early: int = 0,
+        tw_late: int = 0,
+    ) -> None: ...
 
 class ProblemData:
     """
@@ -39,45 +51,24 @@ class ProblemData:
 
     Parameters
     ----------
-    coords
-        Array of (x, y) coordinates. The first coordinate at index 0 is assumed
-        to be the depot.
-    demands
-        Array of client demands. The demand at index 0 is assumed to be the
-        depot's demand, and should be zero.
+    clients
+        List of clients. The first client (at index 0) is assumed to be the
+        depot. The time window for the depot is assumed to describe the overall
+        time horizon. The depot should have 0 demand and 0 service duration.
     nb_vehicles
         The number of vehicles in this problem instance.
     vehicle_cap
         Homogenous vehicle capacity for all vehicles in the problem instance.
-    time_windows
-        Array of (early, late) time windows. The time window at index 0 is
-        assumed to be the depot's time window, and describes the overall time
-        horizon.
-    service_durations
-        Array of service durations, that is, the length of time needed to
-        service a customer upon visiting. The service duration at index 0 is
-        assumed to be the depot's service time, and should be zero.
-    distance_matrix
-        A matrix that gives the distances between clients (and the depot at
-        index 0).
     duration_matrix
         A matrix that gives the travel times between clients (and the depot at
         index 0).
-
-    Notes
-    -----
-    All array data assume that the data at or involving index 0 relates to the
-    depot, and all other indices specify client information.
     """
 
     def __init__(
         self,
-        coords: List[Tuple[int, int]],
-        demands: List[int],
+        clients: List[Client],
         nb_vehicles: int,
         vehicle_cap: int,
-        time_windows: List[Tuple[int, int]],
-        service_durations: List[int],
         distance_matrix: List[List[int]],
         duration_matrix: List[List[int]],
     ): ...
diff --git a/pyvrp/__init__.py b/pyvrp/__init__.py
index 7af461f..033d52f 100644
--- a/pyvrp/__init__.py
+++ b/pyvrp/__init__.py
@@ -6,7 +6,7 @@ from .Statistics import Statistics
 from ._CostEvaluator import CostEvaluator
 from ._Individual import Individual
 from ._Matrix import Matrix
-from ._ProblemData import ProblemData
+from ._ProblemData import Client, ProblemData
 from ._XorShift128 import XorShift128
 from .read import read, read_solution
 from .show_versions import show_versions
diff --git a/pyvrp/cpp/ProblemData.cpp b/pyvrp/cpp/ProblemData.cpp
index 1b29787..d178966 100644
--- a/pyvrp/cpp/ProblemData.cpp
+++ b/pyvrp/cpp/ProblemData.cpp
@@ -6,6 +6,25 @@
 #include <string>
 #include <vector>
 
+ProblemData::Client::Client(
+    int x, int y, int demand, int serviceDuration, int twEarly, int twLate)
+    : x(x),
+      y(y),
+      demand(demand),
+      serviceDuration(serviceDuration),
+      twEarly(twEarly),
+      twLate(twLate)
+{
+    if (demand < 0)
+        throw std::invalid_argument("demand must be >= 0");
+
+    if (serviceDuration < 0)
+        throw std::invalid_argument("service_duration must be >= 0");
+
+    if (twEarly > twLate)
+        throw std::invalid_argument("tw_early must be <= tw_late");
+}
+
 ProblemData::Client const &ProblemData::depot() const { return client(0); }
 
 Matrix<int> const &ProblemData::distanceMatrix() const { return dist_; }
@@ -18,28 +37,16 @@ size_t ProblemData::numVehicles() const { return numVehicles_; }
 
 size_t ProblemData::vehicleCapacity() const { return vehicleCapacity_; }
 
-ProblemData::ProblemData(std::vector<std::pair<int, int>> const &coords,
-                         std::vector<int> const &demands,
+ProblemData::ProblemData(std::vector<Client> const &clients,
                          size_t numVehicles,
                          size_t vehicleCap,
-                         std::vector<std::pair<int, int>> const &timeWindows,
-                         std::vector<int> const &servDurs,
                          std::vector<std::vector<int>> const &distMat,
                          std::vector<std::vector<int>> const &durMat)
     : dist_(distMat),
       dur_(durMat),
-      clients_(coords.size()),
-      numClients_(static_cast<int>(coords.size()) - 1),
+      clients_(clients),
+      numClients_(std::max(clients.size(), static_cast<size_t>(1)) - 1),
       numVehicles_(numVehicles),
       vehicleCapacity_(vehicleCap)
 {
-    // TODO argument checks
-
-    for (size_t idx = 0; idx <= static_cast<size_t>(numClients_); ++idx)
-        clients_[idx] = {coords[idx].first,
-                         coords[idx].second,
-                         servDurs[idx],
-                         demands[idx],
-                         timeWindows[idx].first,
-                         timeWindows[idx].second};
 }
diff --git a/pyvrp/cpp/ProblemData.h b/pyvrp/cpp/ProblemData.h
index ebfd8b5..96d28ab 100644
--- a/pyvrp/cpp/ProblemData.h
+++ b/pyvrp/cpp/ProblemData.h
@@ -14,10 +14,17 @@ public:
     {
         int x;                // Coordinate X
         int y;                // Coordinate Y
-        int serviceDuration;  // Service duration
         int demand;           // Demand
+        int serviceDuration;  // Service duration
         int twEarly;          // Earliest arrival (when using time windows)
         int twLate;           // Latest arrival (when using time windows)
+
+        Client(int x,
+               int y,
+               int demand = 0,
+               int serviceDuration = 0,
+               int twEarly = 0,
+               int twLate = 0);
     };
 
 private:
@@ -85,25 +92,19 @@ public:
     [[nodiscard]] size_t vehicleCapacity() const;
 
     /**
-     * Constructs a ProblemData object with the given data. Assumes the data
-     * contains the depot, such that each vector is one longer than the number
-     * of clients.
+     * Constructs a ProblemData object with the given data. Assumes the list of
+     * clients contains the depot, such that each vector is one longer than the
+     * number of clients.
      *
-     * @param coords       Coordinates as pairs of [x, y].
-     * @param demands      Client demands.
+     * @param clients      List of clients (including depot at index 0).
      * @param numVehicles  Number of vehicles.
      * @param vehicleCap   Vehicle capacity.
-     * @param timeWindows  Time windows as pairs of [early, late].
-     * @param servDurs     Service durations.
      * @param distMat      Distance matrix.
      * @param durMat       Duration matrix.
      */
-    ProblemData(std::vector<std::pair<int, int>> const &coords,
-                std::vector<int> const &demands,
+    ProblemData(std::vector<Client> const &clients,
                 size_t numVehicles,
                 size_t vehicleCap,
-                std::vector<std::pair<int, int>> const &timeWindows,
-                std::vector<int> const &servDurs,
                 std::vector<std::vector<int>> const &distMat,
                 std::vector<std::vector<int>> const &durMat);
 };
diff --git a/pyvrp/cpp/ProblemData_bindings.cpp b/pyvrp/cpp/ProblemData_bindings.cpp
index 239019b..c2f77bd 100644
--- a/pyvrp/cpp/ProblemData_bindings.cpp
+++ b/pyvrp/cpp/ProblemData_bindings.cpp
@@ -8,6 +8,13 @@ namespace py = pybind11;
 PYBIND11_MODULE(_ProblemData, m)
 {
     py::class_<ProblemData::Client>(m, "Client")
+        .def(py::init<int, int, int, int, int, int>(),
+             py::arg("x"),
+             py::arg("y"),
+             py::arg("demand") = 0,
+             py::arg("service_duration") = 0,
+             py::arg("tw_early") = 0,
+             py::arg("tw_late") = 0)
         .def_readonly("x", &ProblemData::Client::x)
         .def_readonly("y", &ProblemData::Client::y)
         .def_readonly("service_duration", &ProblemData::Client::serviceDuration)
@@ -16,20 +23,14 @@ PYBIND11_MODULE(_ProblemData, m)
         .def_readonly("tw_late", &ProblemData::Client::twLate);
 
     py::class_<ProblemData>(m, "ProblemData")
-        .def(py::init<std::vector<std::pair<int, int>> const &,
-                      std::vector<int> const &,
+        .def(py::init<std::vector<ProblemData::Client> const &,
                       int,
                       int,
-                      std::vector<std::pair<int, int>> const &,
-                      std::vector<int> const &,
                       std::vector<std::vector<int>> const &,
                       std::vector<std::vector<int>> const &>(),
-             py::arg("coords"),
-             py::arg("demands"),
+             py::arg("clients"),
              py::arg("nb_vehicles"),
              py::arg("vehicle_cap"),
-             py::arg("time_windows"),
-             py::arg("service_durations"),
              py::arg("distance_matrix"),
              py::arg("duration_matrix"))
         .def_property_readonly("num_clients", &ProblemData::numClients)
diff --git a/pyvrp/educate/neighbourhood.py b/pyvrp/educate/neighbourhood.py
index 72036ca..4aaae34 100644
--- a/pyvrp/educate/neighbourhood.py
+++ b/pyvrp/educate/neighbourhood.py
@@ -18,17 +18,24 @@ class NeighbourhoodParams:
     Attributes
     ----------
     weight_wait_time
-        TODO
+        Weight given to the minimum wait time aspect of the proximity
+        calculation. A large wait time indicates the clients are far apart
+        in duration/time.
     weight_time_warp
-        TODO
+        Weight given to the minimum time warp aspect of the proximity
+        calculation. A large time warp indicates the clients are far apart in
+        duration/time.
     nb_granular
         Number of other clients that are in each client's granular
         neighbourhood. This parameter determines the size of the overall
         neighbourhood.
     symmetric_proximity
-        TODO
+        Whether to calculate a symmetric proximity matrix. This ensures edge
+        :math:`(i, j)` is given the same weight as :math:`(j, i)`.
     symmetric_neighbours
-        TODO
+        Whether to symmetrise the neighbourhood structure. This ensures that
+        when edge :math:`(i, j)` is in, then so is :math:`(j, i)`. Note that
+        this is *not* the same as ``symmetric_proximity``.
 
     Raises
     ------
@@ -66,39 +73,43 @@ def compute_neighbours(
         A list of list of integers representing the neighbours for each client.
         The first element represents the depot and is an empty list.
     """
+    proximity = _compute_proximity(
+        data,
+        params.weight_wait_time,
+        params.weight_time_warp,
+    )
 
-    proximity = _compute_proximity(data, params)
+    if params.symmetric_proximity:
+        proximity = np.minimum(proximity, proximity.T)
 
-    # Mask self from proximities by adding a large constant
-    # and add the index with a very small weight as tie-breaker (to mimic
-    # original c++ version)
     n = len(proximity)
-    k = min(params.nb_granular, n - 1)
-    idcs = np.arange(n)
-    jitter = 1e-6 * idcs  # add a bit of jitter to break ties
-    proximity = proximity + jitter[None, :]
-    np.fill_diagonal(proximity, np.inf)  # exclude self from neighbourhood
-    idx_topk = np.argpartition(proximity, k, axis=-1)[:, :k]
+    k = min(params.nb_granular, n - 2)  # excl. depot and self
+
+    np.fill_diagonal(proximity, np.inf)  # cannot be in own neighbourhood
+    proximity[0, :] = np.inf  # depot has no neighbours
+    proximity[:, 0] = np.inf  # clients do not neighbour depot
 
-    if params.symmetric_neighbours:
-        # Convert into adjacency matrix that can be symmetrized
-        adj = np.zeros_like(proximity, dtype=bool)
-        adj[idcs[:, None], idx_topk] = True
+    top_k = np.argsort(proximity, axis=1, kind="stable")[1:, :k]  # excl. depot
 
-        # Add correlated vertex if correlated in one of two directions
-        adj = adj | adj.transpose()
+    if not params.symmetric_neighbours:
+        return [[]] + top_k.tolist()
 
-        # Append empty neighbours for depot and add 1 to offset depot indexing
-        return [[]] + [(np.flatnonzero(row) + 1).tolist() for row in adj]
-    else:
-        return [[]] + (np.sort(idx_topk, -1) + 1).tolist()
+    # Construct a symmetric adjacency matrix and return the adjacent clients
+    # as the neighbourhood structure.
+    adj = np.zeros_like(proximity, dtype=bool)
+    rows = np.expand_dims(np.arange(1, n), axis=1)
+    adj[rows, top_k] = True
+    adj = adj | adj.transpose()
+
+    return [np.flatnonzero(row).tolist() for row in adj]
 
 
 def _compute_proximity(
-    data: ProblemData, params: NeighbourhoodParams = NeighbourhoodParams()
+    data: ProblemData, weight_wait_time: float, weight_time_warp: float
 ) -> np.ndarray[float]:
     """
-    Computes proximity for neighborhood.
+    Computes proximity for neighborhood. Proximity is based on Vidal et al.
+    (2013).
 
     Parameters
     ----------
@@ -112,10 +123,15 @@ def _compute_proximity(
     np.ndarray[float]
         A numpy array of size n x n where n = data.num_clients containing
         the proximities values between all clients (depot excluded).
-    """
 
+    References
+    ----------
+    .. [1] Vidal, T., Crainic, T. G., Gendreau, M., and Prins, C. (2013). A
+           hybrid genetic algorithm with adaptive diversity management for a
+           large class of vehicle routing problems with time-windows.
+           *Computers & Operations Research*, 40(1), 475 - 489.
+    """
     dim = data.num_clients + 1
-
     clients = [data.client(idx) for idx in range(dim)]
 
     earliest = np.array([cli.tw_early for cli in clients])
@@ -129,6 +145,7 @@ def _compute_proximity(
     min_wait_time = np.maximum(
         earliest[None, :] - durations - service[:, None] - latest[:, None], 0
     )
+
     min_time_warp = np.maximum(
         earliest[:, None] + service[:, None] + durations - latest[None, :],
         0,
@@ -139,12 +156,8 @@ def _compute_proximity(
         dtype=float,
     )
 
-    prox = (
+    return (
         distances
-        + params.weight_wait_time * min_wait_time
-        + params.weight_time_warp * min_time_warp
+        + weight_wait_time * min_wait_time
+        + weight_time_warp * min_time_warp
     )
-
-    if params.symmetric_proximity:
-        prox = np.minimum(prox, prox.T)
-    return prox[1:, 1:]  # Skip depot
diff --git a/pyvrp/read.py b/pyvrp/read.py
index fb6b5b1..5a8f535 100644
--- a/pyvrp/read.py
+++ b/pyvrp/read.py
@@ -5,7 +5,7 @@ from typing import Callable, Dict, List, Union
 import numpy as np
 import vrplib
 
-from ._ProblemData import ProblemData
+from ._ProblemData import Client, ProblemData
 
 _Routes = List[List[int]]
 _RoundingFunc = Callable[[np.ndarray], np.ndarray]
@@ -71,9 +71,10 @@ def read(
     ProblemData
         Data instance constructed from the read data.
     """
-    if isinstance(round_func, str) and round_func in ROUND_FUNCS:
-        round_func = ROUND_FUNCS[round_func]
-    elif not callable(round_func):
+    if (key := str(round_func)) in ROUND_FUNCS:
+        round_func = ROUND_FUNCS[key]
+
+    if not callable(round_func):
         raise ValueError(
             f"round_func = {round_func} is not understood. Can be a function,"
             f" or one of {ROUND_FUNCS.keys()}."
@@ -83,41 +84,44 @@ def read(
 
     # A priori checks
     if "dimension" in instance:
-        num_clients = instance["dimension"]
+        dimension: int = instance["dimension"]
     else:
         if "demand" not in instance:
             raise ValueError("File should either contain dimension or demands")
-        num_clients = len(instance["demand"])
+        dimension = len(instance["demand"])
+
+    depots: np.ndarray = instance.get("depot", np.array([0]))
+    num_vehicles: int = instance.get("vehicles", dimension - 1)
+    capacity: int = instance.get("capacity", _INT_MAX)
 
-    depots = instance.get("depot", np.array([0]))
-    num_vehicles = instance.get("vehicles", num_clients - 1)
-    capacity = instance.get("capacity", _INT_MAX)
-    edge_weight = round_func(instance["edge_weight"])
+    distances: np.ndarray = round_func(instance["edge_weight"])
 
     if "demand" in instance:
-        demands = instance["demand"]
+        demands: np.ndarray = instance["demand"]
     else:
-        demands = np.zeros(num_clients, dtype=int)
+        demands = np.zeros(dimension, dtype=int)
 
     if "node_coord" in instance:
-        coords = round_func(instance["node_coord"])
+        coords: np.ndarray = round_func(instance["node_coord"])
     else:
-        coords = np.zeros((num_clients, 2), dtype=int)
+        coords = np.zeros((dimension, 2), dtype=int)
 
     if "service_time" in instance:
-        service_times = round_func(instance["service_time"])
+        service_times: np.ndarray = round_func(instance["service_time"])
     else:
-        service_times = np.zeros(num_clients, dtype=int)
+        service_times = np.zeros(dimension, dtype=int)
 
     if "time_window" in instance:
-        time_windows = round_func(instance["time_window"])
+        # VRPLIB instances typically do not have a duration data field, so we
+        # assume duration == distance if the instance has time windows.
+        durations = distances
+        time_windows: np.ndarray = round_func(instance["time_window"])
     else:
-        # The default value for the latest time window based on the maximum
-        # route duration. This ensures that the time window constraints are
-        # always satisfied.
-        bound = num_clients * (edge_weight.max() + service_times.max())
-        bound = min(bound, _INT_MAX)
-        time_windows = np.repeat([[0, bound]], num_clients, axis=0)
+        # No time window data. So the time window component is not relevant,
+        # and we set all time-related fields to zero.
+        durations = np.zeros_like(distances)
+        service_times = np.zeros(dimension, dtype=int)
+        time_windows = np.zeros((dimension, 2), dtype=int)
 
     # Checks
     if len(depots) != 1 or depots[0] != 0:
@@ -138,15 +142,24 @@ def read(
     if (time_windows[:, 0] > time_windows[:, 1]).any():
         raise ValueError("Time window cannot start after end")
 
+    clients = [
+        Client(
+            coords[idx][0],  # x
+            coords[idx][1],  # y
+            demands[idx],
+            service_times[idx],
+            time_windows[idx][0],  # TW early
+            time_windows[idx][1],  # TW late
+        )
+        for idx in range(dimension)
+    ]
+
     return ProblemData(
-        coords,
-        demands,
+        clients,
         num_vehicles,
         capacity,
-        time_windows,
-        service_times,
-        edge_weight,  # distance and duration are assumed to be the same in
-        edge_weight,  # the VRPLIB instances we read here.
+        distances,
+        durations,
     )
 
 


Test patch:
diff --git a/pyvrp/educate/tests/test_Exchange.py b/pyvrp/educate/tests/test_Exchange.py
index 80fddc0..bd8690b 100644
--- a/pyvrp/educate/tests/test_Exchange.py
+++ b/pyvrp/educate/tests/test_Exchange.py
@@ -1,7 +1,7 @@
 from numpy.testing import assert_, assert_equal
 from pytest import mark
 
-from pyvrp import CostEvaluator, Individual, ProblemData, XorShift128
+from pyvrp import Client, CostEvaluator, Individual, ProblemData, XorShift128
 from pyvrp.educate import LocalSearch, NeighbourhoodParams, compute_neighbours
 from pyvrp.educate._Exchange import (
     Exchange10,
@@ -213,13 +213,16 @@ def test_relocate_only_happens_when_distance_and_duration_allow_it():
     Tests that (1, 0)-exchange checks the duration matrix for time-window
     feasibility before applying a move that improves the traveled distance.
     """
+    clients = [
+        Client(x=0, y=0, demand=0, service_duration=0, tw_early=0, tw_late=10),
+        Client(x=1, y=0, demand=0, service_duration=0, tw_early=0, tw_late=5),
+        Client(x=2, y=0, demand=0, service_duration=0, tw_early=0, tw_late=5),
+    ]
+
     data = ProblemData(
-        coords=[(0, 0), (1, 0), (2, 0)],
-        demands=[0, 0, 0],
+        clients=clients,
         nb_vehicles=1,
         vehicle_cap=0,
-        time_windows=[(0, 10), (0, 5), (0, 5)],
-        service_durations=[0, 0, 0],
         distance_matrix=[  # distance-wise, the best route is 0 -> 1 -> 2 -> 0.
             [0, 1, 5],
             [5, 0, 1],
diff --git a/pyvrp/educate/tests/test_MoveTwoClientsReversed.py b/pyvrp/educate/tests/test_MoveTwoClientsReversed.py
index 19c0256..7084399 100644
--- a/pyvrp/educate/tests/test_MoveTwoClientsReversed.py
+++ b/pyvrp/educate/tests/test_MoveTwoClientsReversed.py
@@ -34,13 +34,13 @@ def test_single_route_OkSmall():
     assert_(improved_cost < current_cost)
 
     # (2, 3) was inserted after 1 as 1 -> 3 -> 2 -> 4. Then (1, 3) got inserted
-    # after 2 as 2 -> 3 -> 1 -> 4.
-    assert_equal(improved_individual.get_routes(), [[2, 3, 1, 4], [], []])
+    # after 4 as 2 -> 4 -> 3 -> 1.
+    assert_equal(improved_individual.get_routes(), [[2, 4, 3, 1], [], []])
 
     # These two-route solutions can all be created by MoveTwoClientsReversed
     # from the returned solution. So they must have a cost that's at best equal
     # to the returned solution's cost.
-    for routes in ([[3, 2], [1, 4]], [[2, 3], [4, 1]], [[2, 4], [1, 3]]):
+    for routes in ([[3, 1], [4, 2]], [[2, 1], [3, 4]], [[2, 4], [1, 3]]):
         other = Individual(data, routes)
         improved_cost = cost_evaluator.penalised_cost(improved_individual)
         other_cost = cost_evaluator.penalised_cost(other)
diff --git a/pyvrp/educate/tests/test_neighbourhood.py b/pyvrp/educate/tests/test_neighbourhood.py
index 1c3000c..721c925 100644
--- a/pyvrp/educate/tests/test_neighbourhood.py
+++ b/pyvrp/educate/tests/test_neighbourhood.py
@@ -1,5 +1,6 @@
-from typing import List
+from typing import Set
 
+import numpy as np
 from numpy.testing import assert_, assert_equal, assert_raises
 from pytest import mark
 
@@ -75,16 +76,16 @@ def test_neighbourhood_params_does_not_raise_for_valid_arguments(
     [
         # fmt: off
         (20, 20, 10, True, False, 2,
-         [1, 3, 4, 5, 6, 7, 8, 45, 46, 100]),
+         {1, 3, 4, 5, 6, 7, 8, 45, 46, 100}),
         (20, 20, 10, True, True, 2,
-         [1, 3, 4, 5, 6, 7, 8, 45, 46, 60, 70, 79, 100]),
+         {1, 3, 4, 5, 6, 7, 8, 45, 46, 60, 70, 79, 100}),
         # From original c++ implementation
         (18, 20, 34, True, False, 1,
-         [2, 3, 4, 5, 6, 7, 8, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 53,
-          54, 55, 60, 61, 68, 69, 70, 73, 78, 79, 81, 88, 90, 98, 100]),
+         {2, 3, 4, 5, 6, 7, 8, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 53,
+          54, 55, 60, 61, 68, 69, 70, 73, 78, 79, 81, 88, 90, 98, 100}),
         (18, 20, 34, True, False, 99,
-         [9, 10, 11, 12, 13, 14, 15, 16, 20, 22, 24, 47, 52, 53, 55, 56, 57,
-          58, 59, 60, 64, 65, 66, 69, 74, 80, 82, 83, 86, 87, 88, 90, 91, 98]),
+         {9, 10, 11, 12, 13, 14, 15, 16, 20, 22, 24, 47, 52, 53, 55, 56, 57,
+          58, 59, 60, 64, 65, 66, 69, 74, 80, 82, 83, 86, 87, 88, 90, 91, 98}),
         # fmt: on
     ],
 )
@@ -95,7 +96,7 @@ def test_compute_neighbours(
     symmetric_proximity: bool,
     symmetric_neighbours: bool,
     idx_check: int,
-    expected_neighbours_check: List[int],
+    expected_neighbours_check: Set[int],
 ):
     data = read("data/RC208.txt", "solomon", "trunc1")
     params = NeighbourhoodParams(
@@ -110,9 +111,10 @@ def test_compute_neighbours(
     assert_equal(len(neighbours), data.num_clients + 1)
     assert_equal(len(neighbours[0]), 0)
 
-    assert_equal(neighbours[idx_check], expected_neighbours_check)
+    # We compare sets because the expected data (from the old C++
+    # implementation) sorts by client ID (ascending), not proximity.
+    assert_equal(set(neighbours[idx_check]), expected_neighbours_check)
 
-    # Check that we have at least or exactly nb_granular depending on symmetric
     for neighb in neighbours[1:]:
         if symmetric_neighbours:
             assert_(len(neighb) >= nb_granular)
@@ -120,6 +122,40 @@ def test_compute_neighbours(
             assert_equal(len(neighb), nb_granular)
 
 
+def test_neighbours_are_sorted_by_proximity():
+    # E-n22-k4 is a CVRP instance, so the time component is irrelevant.
+    data = read("data/E-n22-k4.txt", round_func="trunc1")
+    params = NeighbourhoodParams(0, 0, data.num_clients)
+    neighbours = compute_neighbours(data, params)
+    clients = list(range(1, data.num_clients + 1))
+
+    for client in clients:
+        # Proximity is completely based on distance. We break ties by index
+        # (using stable sort). The test below checks that this is the same as
+        # what comes out of the granular neighbourhood calculation.
+        valid = np.array([other for other in clients if other != client])
+        dists = [data.dist(client, other) for other in valid]
+        by_proximity = valid[np.argsort(dists, kind="stable")]
+        assert_equal(by_proximity, neighbours[client])
+
+
+def test_symmetric_neighbours():
+    data = read("data/RC208.txt", "solomon", "trunc1")
+
+    # Symmetric neighbourhood structure: if (i, j) is in, then so is (j, i).
+    params = NeighbourhoodParams(symmetric_neighbours=True)
+    sym_neighbours = [set(n) for n in compute_neighbours(data, params)]
+
+    for client in range(data.num_clients + 1):
+        for neighbour in sym_neighbours[client]:
+            assert_(client in sym_neighbours[neighbour])
+
+    # But when neighbours are not symmetrised, this is typically not the case.
+    params = NeighbourhoodParams(symmetric_neighbours=False)
+    asym_neighbours = [set(n) for n in compute_neighbours(data, params)]
+    assert_(sym_neighbours != asym_neighbours)
+
+
 def test_more_neighbours_than_instance_size():
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     params = NeighbourhoodParams(nb_granular=data.num_clients)
diff --git a/pyvrp/tests/test_Individual.py b/pyvrp/tests/test_Individual.py
index 4d85f94..b946455 100644
--- a/pyvrp/tests/test_Individual.py
+++ b/pyvrp/tests/test_Individual.py
@@ -4,7 +4,7 @@ import numpy as np
 from numpy.testing import assert_, assert_equal, assert_raises
 from pytest import mark
 
-from pyvrp import Individual, ProblemData, XorShift128
+from pyvrp import Client, Individual, ProblemData, XorShift128
 from pyvrp.tests.helpers import read
 
 
@@ -155,12 +155,13 @@ def test_time_warp_for_a_very_constrained_problem(dist_mat):
     ]
 
     data = ProblemData(
-        coords=[(0, 0), (1, 0), (2, 0)],
-        demands=[0, 0, 0],
+        clients=[
+            Client(x=0, y=0, tw_late=10),
+            Client(x=1, y=0, tw_late=5),
+            Client(x=2, y=0, tw_late=5),
+        ],
         nb_vehicles=2,
         vehicle_cap=0,
-        time_windows=[(0, 10), (0, 5), (0, 5)],
-        service_durations=[0, 0, 0],
         distance_matrix=dist_mat,
         duration_matrix=dur_mat,
     )
diff --git a/pyvrp/tests/test_ProblemData.py b/pyvrp/tests/test_ProblemData.py
index 71b928d..670d9e5 100644
--- a/pyvrp/tests/test_ProblemData.py
+++ b/pyvrp/tests/test_ProblemData.py
@@ -1,8 +1,55 @@
-import numpy as np
 from numpy.random import default_rng
-from numpy.testing import assert_, assert_equal
+from numpy.testing import assert_, assert_allclose, assert_raises
+from pytest import mark
 
-from pyvrp import ProblemData
+from pyvrp import Client, ProblemData
+
+
+@mark.parametrize(
+    "x,y,demand,service_duration,tw_early,tw_late",
+    [
+        (1, 1, 1, 1, 0, 1),  # normal
+        (1, 1, 1, 0, 0, 1),  # zero duration
+        (1, 1, 0, 1, 0, 1),  # zero demand
+        (1, 1, 1, 1, 0, 0),  # zero length time interval
+        (-1, -1, 1, 1, 0, 1),  # negative coordinates
+    ],
+)
+def test_client_constructor_initialises_data_fields_correctly(
+    x: int,
+    y: int,
+    demand: int,
+    service_duration: int,
+    tw_early: int,
+    tw_late: int,
+):
+    client = Client(x, y, demand, service_duration, tw_early, tw_late)
+    assert_allclose(client.x, x)
+    assert_allclose(client.y, y)
+    assert_allclose(client.demand, demand)
+    assert_allclose(client.service_duration, service_duration)
+    assert_allclose(client.tw_early, tw_early)
+    assert_allclose(client.tw_late, tw_late)
+
+
+@mark.parametrize(
+    "x,y,demand,service_duration,tw_early,tw_late",
+    [
+        (1, 1, 1, 1, 1, 0),  # late < early
+        (1, 1, 1, -1, 0, 1),  # negative duration
+        (1, 1, -1, 1, 0, 1),  # negative demand
+    ],
+)
+def test_raises_for_invalid_client_data(
+    x: int,
+    y: int,
+    demand: int,
+    service_duration: int,
+    tw_early: int,
+    tw_late: int,
+):
+    with assert_raises(ValueError):
+        Client(x, y, demand, service_duration, tw_early, tw_late)
 
 
 def test_depot_is_first_client():
@@ -11,13 +58,11 @@ def test_depot_is_first_client():
     ``client(0)``.
     """
     mat = [[0, 1], [1, 0]]
+
     data = ProblemData(
-        coords=[(0, 0), (0, 1)],
-        demands=[0, 0],
+        clients=[Client(x=0, y=0), Client(x=0, y=1)],
         nb_vehicles=1,
         vehicle_cap=1,
-        time_windows=[(0, 10), (0, 10)],
-        service_durations=[0, 0],
         distance_matrix=mat,
         duration_matrix=mat,
     )
@@ -31,20 +76,20 @@ def test_matrix_access():
     ``distance_matrix()`` methods correctly index the underlying duration
     and distance matrices.
     """
-    size = 6
     gen = default_rng(seed=42)
+    size = 6
+
     dist_mat = gen.integers(500, size=(size, size))
     dur_mat = gen.integers(500, size=(size, size))
-
-    assert_(not np.allclose(dist_mat, dur_mat))
+    clients = [
+        Client(x=0, y=0, demand=0, service_duration=0, tw_early=0, tw_late=10)
+        for _ in range(size)
+    ]
 
     data = ProblemData(
-        coords=[(0, 0)] * size,
-        demands=[0] * size,
+        clients=clients,
         nb_vehicles=1,
         vehicle_cap=1,
-        time_windows=[(0, 10)] * size,
-        service_durations=[0] * size,
         distance_matrix=dist_mat,  # type: ignore
         duration_matrix=dur_mat,  # type: ignore
     )
@@ -54,8 +99,8 @@ def test_matrix_access():
 
     for frm in range(size):
         for to in range(size):
-            assert_equal(dur_mat[frm, to], data.duration(frm, to))
-            assert_equal(dist_mat[frm, to], data.dist(frm, to))
+            assert_allclose(dur_mat[frm, to], data.duration(frm, to))
+            assert_allclose(dist_mat[frm, to], data.dist(frm, to))
 
-            assert_equal(dur_mat_data[frm, to], dur_mat[frm, to])
-            assert_equal(dist_mat_data[frm, to], dist_mat[frm, to])
+            assert_allclose(dur_mat_data[frm, to], dur_mat[frm, to])
+            assert_allclose(dist_mat_data[frm, to], dist_mat[frm, to])
diff --git a/pyvrp/tests/test_read.py b/pyvrp/tests/test_read.py
index 1ca75d7..7e7f86e 100644
--- a/pyvrp/tests/test_read.py
+++ b/pyvrp/tests/test_read.py
@@ -128,19 +128,11 @@ def test_reading_En22k4_instance():  # instance from CVRPLIB
     assert_equal(data.dist(0, 1), 493)
     assert_equal(data.dist(1, 0), 493)
 
-    # These fields are not present in the data file, and should thus retain
-    # their default values.
-    max_distance = max(
-        data.dist(i, j)
-        for i in range(data.num_clients + 1)
-        for j in range(data.num_clients + 1)
-    )
-    bound = (data.num_clients + 1) * max_distance
-
+    # This is a CVRP instance, so time window data should all be zeroed out.
     for client in range(data.num_clients + 1):  # incl. depot
         assert_equal(data.client(client).service_duration, 0)
         assert_equal(data.client(client).tw_early, 0)
-        assert_equal(data.client(client).tw_late, bound)
+        assert_equal(data.client(client).tw_late, 0)
 
 
 def test_reading_RC208_instance():  # Solomon style instance