Repository: PyVRP/PyVRP
Pull request #240 (2023-06-12 14:13:11)

Only store non-empty routes in Individual
This PR:

- [x] Removes empty routes from Individual.
- [x] Adds and passes relevant tests.
- [x] Has well-formatted code and documentation.

### Benchmark

On VRPTW.

This PR (e25d6c35501095d17a90621dfa91be64ed4a7ad4; job ID 2530076; 10x seeds/1h with `configs/vrptw.toml`):
```
 Avg. objective: 333241.4 (min: 333117, max: 333411)
Avg. iterations: 115161.3 (min: 103903, max: 128743)
```


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) Report
Patch coverage: **`100.00`**% and project coverage change: **`-0.13`** :warning:
> Comparison is base [(`2feb931`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/2feb931848bf3720f9befce6f1a4c7a8c86d9589?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) 93.19% compared to head [(`5f198d7`)](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) 93.07%.

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


```diff
@@            Coverage Diff             @@
##             main     #240      +/-   ##
==========================================
- Coverage   93.19%   93.07%   -0.13%     
==========================================
  Files          61       61              
  Lines        2279     2280       +1     
  Branches      196      197       +1     
==========================================
- Hits         2124     2122       -2     
- Misses         64       66       +2     
- Partials       91       92       +1     
```


| [Impacted Files](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) | Coverage Δ | |
|---|---|---|
| [pyvrp/cpp/Individual.cpp](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP#diff-cHl2cnAvY3BwL0luZGl2aWR1YWwuY3Bw) | `91.39% <100.00%> (-1.86%)` | :arrow_down: |
| [pyvrp/cpp/educate/LocalSearch.cpp](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP#diff-cHl2cnAvY3BwL2VkdWNhdGUvTG9jYWxTZWFyY2guY3Bw) | `91.58% <100.00%> (ø)` | |
| [pyvrp/plotting/plot\_solution.py](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP#diff-cHl2cnAvcGxvdHRpbmcvcGxvdF9zb2x1dGlvbi5weQ==) | `100.00% <100.00%> (ø)` | |

... and [1 file with indirect coverage changes](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP)


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PyVRP/PyVRP/pull/240?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).

wouterkool: > This does not seem done to me yet. What about:
> 
> * SREX: sort asc angle no need to check empty
> * greedyRepair: numRoutes is just routes.size() (?)
> * Individual::evaluate: break unnecessary
> * Docstrings on e.g. Individual::numRoutes, ::getRoutes, and ::routes_, and possibly other places.
> 
> (there might be other places, but a search for e.g. num_routes/numRoutes and similar should turn that up)

You were right this was done a bit quick, thanks for the suggestions. I checked for num_routes/numRoutes and empty so I think all should be covered now. Note in the crossover greedy repair we may still have empty routes before repairing.
N-Wouda: Since this "routes have length `data.numVehicles()` thing has been an assumption for a very long time, I'm running a benchmark to be sure we got every odd corner case. Results should be in in 3h. If everything looks good I think we can merge it, no problemo.


Patch:
diff --git a/docs/source/api/pyvrp.rst b/docs/source/api/pyvrp.rst
index bc77bfa..9161b1d 100644
--- a/docs/source/api/pyvrp.rst
+++ b/docs/source/api/pyvrp.rst
@@ -15,6 +15,11 @@ This object can be used to obtain the best observed solution, and detailed runti
 
    Have a look at the examples to see how these classes relate!
 
+.. automodule:: pyvrp.Model
+
+   .. autoclass:: Model
+      :members:
+
 .. automodule:: pyvrp._CostEvaluator
    
    .. autoapiclass:: CostEvaluator
@@ -22,10 +27,10 @@ This object can be used to obtain the best observed solution, and detailed runti
 
 .. automodule:: pyvrp.GeneticAlgorithm
 
-   .. autoapiclass:: GeneticAlgorithmParams
+   .. autoclass:: GeneticAlgorithmParams
       :members:
 
-   .. autoapiclass:: GeneticAlgorithm
+   .. autoclass:: GeneticAlgorithm
       :members:
 
 .. automodule:: pyvrp._Individual
@@ -38,10 +43,10 @@ This object can be used to obtain the best observed solution, and detailed runti
 
 .. automodule:: pyvrp.PenaltyManager
 
-   .. autoapiclass:: PenaltyParams
+   .. autoclass:: PenaltyParams
       :members:
 
-   .. autoapiclass:: PenaltyManager
+   .. autoclass:: PenaltyManager
       :members: 
 
 .. automodule:: pyvrp.Population
diff --git a/examples/vrptw.ipynb b/examples/vrptw.ipynb
index 88cb6b6..7c5cf28 100644
--- a/examples/vrptw.ipynb
+++ b/examples/vrptw.ipynb
@@ -318,7 +318,6 @@
    ],
    "source": [
     "solution = result.best\n",
-    "routes = [route for route in solution.get_routes() if route]\n",
     "data = [\n",
     "    {\n",
     "        \"num_stops\": len(route),\n",
@@ -328,7 +327,7 @@
     "        \"time_warp\": route.time_warp(),\n",
     "        \"demand\": route.demand(),\n",
     "    }\n",
-    "    for route in routes\n",
+    "    for route in solution.get_routes()\n",
     "]\n",
     "\n",
     "tabulate([datum.values() for datum in data], data[0].keys(), tablefmt=\"html\")"
@@ -349,7 +348,8 @@
    "outputs": [],
    "source": [
     "cost_evaluator = CostEvaluator()\n",
-    "assert sum(r.distance() for r in routes) == cost_evaluator.cost(solution)"
+    "total_distance = sum(r.distance() for r in solution.get_routes())\n",
+    "assert total_distance == cost_evaluator.cost(solution)"
    ]
   },
   {
diff --git a/pyvrp/Model.py b/pyvrp/Model.py
new file mode 100644
index 0000000..e22c924
--- /dev/null
+++ b/pyvrp/Model.py
@@ -0,0 +1,283 @@
+from typing import List, Union
+from warnings import warn
+
+import numpy as np
+
+from pyvrp.GeneticAlgorithm import GeneticAlgorithm, GeneticAlgorithmParams
+from pyvrp.PenaltyManager import PenaltyManager
+from pyvrp.Population import Population, PopulationParams
+from pyvrp.Result import Result
+from pyvrp._Individual import Individual
+from pyvrp._ProblemData import Client, ProblemData
+from pyvrp._XorShift128 import XorShift128
+from pyvrp.constants import MAX_USER_VALUE, MAX_VALUE
+from pyvrp.crossover import selective_route_exchange as srex
+from pyvrp.diversity import broken_pairs_distance as bpd
+from pyvrp.exceptions import ScalingWarning
+from pyvrp.stop import StoppingCriterion
+
+Depot = Client
+
+
+class Edge:
+
+    __slots__ = ["frm", "to", "distance", "duration"]
+
+    def __init__(self, frm: Client, to: Client, distance: int, duration: int):
+        self.frm = frm
+        self.to = to
+        self.distance = distance
+        self.duration = duration
+
+
+class VehicleType:
+
+    __slots__ = ["number", "capacity"]
+
+    def __init__(self, number: int, capacity: int):
+        self.number = number
+        self.capacity = capacity
+
+
+class Model:
+    """
+    A simple interface for modelling vehicle routing problems with PyVRP.
+    """
+
+    def __init__(self):
+        self._clients: List[Client] = []
+        self._depots: List[Depot] = []
+        self._edges: List[Edge] = []
+        self._vehicle_types: List[VehicleType] = []
+
+    @property
+    def locations(self) -> List[Client]:
+        """
+        Returns all locations (depots and clients) in the current model. The
+        routes returned by :meth:`~solve` can be used to index these locations.
+        """
+        return self._depots + self._clients
+
+    @classmethod
+    def from_data(cls, data: ProblemData) -> "Model":
+        """
+        Constructs a model instance from the given data.
+
+        Parameters
+        ----------
+        data
+            Problem data to feed into the model.
+
+        Returns
+        -------
+        Model
+            A model instance representing the given data.
+        """
+        clients = [data.client(idx) for idx in range(data.num_clients + 1)]
+        edges = [
+            Edge(
+                clients[frm],
+                clients[to],
+                data.dist(frm, to),
+                data.duration(frm, to),
+            )
+            for frm in range(data.num_clients + 1)
+            for to in range(data.num_clients + 1)
+        ]
+
+        self = Model()
+        self._clients = clients[1:]
+        self._depots = clients[:1]
+        self._edges = edges
+        vehicle_types = [VehicleType(data.num_vehicles, data.vehicle_capacity)]
+        self._vehicle_types = vehicle_types
+
+        return self
+
+    def add_client(
+        self,
+        x: int,
+        y: int,
+        demand: int = 0,
+        service_duration: int = 0,
+        tw_early: int = 0,
+        tw_late: int = 0,
+        prize: int = 0,
+        required: bool = True,
+    ) -> Client:
+        """
+        Adds a client with the given attributes to the model. Returns the
+        created :class:`~pyvrp._ProblemData.Client` instance.
+        """
+        client = Client(
+            x,
+            y,
+            demand,
+            service_duration,
+            tw_early,
+            tw_late,
+            prize,
+            required,
+        )
+
+        self._clients.append(client)
+        return client
+
+    def add_depot(
+        self, x: int, y: int, tw_early: int = 0, tw_late: int = 0
+    ) -> Depot:
+        """
+        Adds a depot with the given attributes to the model. Returns the
+        created :class:`~pyvrp._ProblemData.Client` instance.
+
+        .. warning::
+
+           PyVRP does not yet support multi-depot VRPs. For now, only one depot
+           can be added to the model.
+        """
+        if len(self._depots) >= 1:
+            msg = "PyVRP does not yet support multi-depot VRPs."
+            raise ValueError(msg)
+
+        depot = Depot(x, y, tw_early=tw_early, tw_late=tw_late)
+        self._depots.append(depot)
+        return depot
+
+    def add_edge(
+        self,
+        frm: Union[Client, Depot],
+        to: Union[Client, Depot],
+        distance: int,
+        duration: int = 0,
+    ) -> Edge:
+        """
+        Adds an edge :math:`(i, j)` between ``frm`` (:math:`i`) and ``to``
+        (:math:`j`). The edge can be given distance and duration attributes.
+        Distance is required, but the default duration is zero. Returns the
+        created edge.
+
+        Raises
+        ------
+        ValueError
+            When either distance or duration is a negative value.
+        """
+        if distance < 0 or duration < 0:
+            raise ValueError("Cannot have negative edge distance or duration.")
+
+        edge = Edge(frm, to, distance, duration)
+        self._edges.append(edge)
+        return edge
+
+    def add_vehicle_type(self, number: int, capacity: int) -> VehicleType:
+        """
+        Adds a vehicle type with the given number of vehicles of given capacity
+        to the model. Returns the created vehicle type.
+
+        .. warning::
+
+           PyVRP does not yet support heterogeneous fleet VRPs. For now, only
+           one vehicle type can be added to the model.
+
+        Raises
+        ------
+        ValueError
+            When either the vehicle number or capacity is not a positive value.
+        """
+        if len(self._vehicle_types) >= 1:
+            msg = "PyVRP does not yet support heterogeneous fleet VRPs."
+            raise ValueError(msg)
+
+        if number <= 0:
+            raise ValueError("Must have positive number of vehicles.")
+
+        if capacity < 0:
+            raise ValueError("Cannot have negative vehicle capacity.")
+
+        vehicle_type = VehicleType(number, capacity)
+        self._vehicle_types.append(vehicle_type)
+        return vehicle_type
+
+    def data(self) -> ProblemData:
+        """
+        Creates and returns a :class:`~pyvrp._ProblemData.ProblemData` instance
+        from this model's attributes.
+        """
+        locs = self.locations
+        loc2idx = {id(loc): idx for idx, loc in enumerate(locs)}
+
+        num_vehicles = self._vehicle_types[0].number
+        vehicle_capacity = self._vehicle_types[0].capacity
+
+        max_data_value = max(max(e.distance, e.duration) for e in self._edges)
+        if max_data_value > MAX_USER_VALUE:
+            msg = """
+            The maximum distance or duration value is very large. This might
+            impact numerical stability. Consider rescaling your input data.
+            """
+            warn(msg, ScalingWarning)
+
+        # Default value is a sufficiently large value to make sure any edges
+        # 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)
+
+        for edge in self._edges:
+            frm = loc2idx[id(edge.frm)]
+            to = loc2idx[id(edge.to)]
+            distances[frm, to] = edge.distance
+            durations[frm, to] = edge.duration
+
+        return ProblemData(
+            locs,
+            num_vehicles,
+            vehicle_capacity,
+            distances,
+            durations,
+        )
+
+    def solve(self, stop: StoppingCriterion, seed: int = 0) -> Result:
+        """
+        Solve this model.
+
+        Parameters
+        ----------
+        stop
+            Stopping criterion to use.
+        seed, optional
+            Seed value to use for the PRNG, by default 0.
+
+        Returns
+        -------
+        Result
+            The solution result object, containing the best found solution.
+        """
+        # These cause a circular import, so the imports needed to be postponed
+        # to here (where they are actually used).
+        from pyvrp.educate import (
+            NODE_OPERATORS,
+            ROUTE_OPERATORS,
+            LocalSearch,
+            compute_neighbours,
+        )
+
+        data = self.data()
+        rng = XorShift128(seed=seed)
+        ls = LocalSearch(data, rng, compute_neighbours(data))
+
+        for op in NODE_OPERATORS:
+            ls.add_node_operator(op(data))
+
+        for op in ROUTE_OPERATORS:
+            ls.add_route_operator(op(data))
+
+        pm = PenaltyManager()
+        pop_params = PopulationParams()
+        pop = Population(bpd, pop_params)
+        init = [
+            Individual.make_random(data, rng)
+            for _ in range(pop_params.min_pop_size)
+        ]
+
+        gen_params = GeneticAlgorithmParams(collect_statistics=True)
+        algo = GeneticAlgorithm(data, pm, rng, pop, ls, srex, init, gen_params)
+        return algo.run(stop)
diff --git a/pyvrp/_Individual.pyi b/pyvrp/_Individual.pyi
index d6dd116..103f223 100644
--- a/pyvrp/_Individual.pyi
+++ b/pyvrp/_Individual.pyi
@@ -117,11 +117,9 @@ class Individual:
 
         .. note::
 
-           This list is of length
-           :py:attr:`~pyvrp._ProblemData.ProblemData.num_vehicles`, but there
-           could be a number of empty routes. These empty routes are all in the
-           higher indices (guarantee). Use :meth:`~num_routes` to determine
-           which of the lower indices contain non-empty routes.
+           The length of this list is equal to the number of non-empty routes,
+           which is at most equal to
+           :py:attr:`~pyvrp._ProblemData.ProblemData.num_vehicles`.
 
         Returns
         -------
@@ -208,12 +206,12 @@ class Individual:
         """
     def num_routes(self) -> int:
         """
-        Number of non-empty routes in this solution.
+        Number of (non-empty) routes in this solution.
 
         Returns
         -------
         int
-            Number of non-empty routes.
+            Number of (non-empty) routes.
         """
     def num_clients(self) -> int:
         """
diff --git a/pyvrp/__init__.py b/pyvrp/__init__.py
index fbab870..42b5771 100644
--- a/pyvrp/__init__.py
+++ b/pyvrp/__init__.py
@@ -1,4 +1,5 @@
 from .GeneticAlgorithm import GeneticAlgorithm, GeneticAlgorithmParams
+from .Model import Model
 from .PenaltyManager import PenaltyManager, PenaltyParams
 from .Population import Population, PopulationParams
 from .Result import Result
diff --git a/pyvrp/constants.py b/pyvrp/constants.py
new file mode 100644
index 0000000..c498a28
--- /dev/null
+++ b/pyvrp/constants.py
@@ -0,0 +1,2 @@
+MAX_USER_VALUE = 2 << 25
+MAX_VALUE = 2 << 28
diff --git a/pyvrp/cpp/Individual.cpp b/pyvrp/cpp/Individual.cpp
index 36f7149..aa895a1 100644
--- a/pyvrp/cpp/Individual.cpp
+++ b/pyvrp/cpp/Individual.cpp
@@ -17,11 +17,7 @@ void Individual::evaluate(ProblemData const &data)
 
     for (auto const &route : routes_)
     {
-        if (route.empty())  // First empty route. All subsequent routes are
-            break;          // empty as well.
-
         // Whole solution statistics.
-        numRoutes_++;
         numClients_ += route.size();
         prizes_ += route.prizes();
         distance_ += route.distance();
@@ -32,7 +28,7 @@ void Individual::evaluate(ProblemData const &data)
     uncollectedPrizes_ = allPrizes - prizes_;
 }
 
-size_t Individual::numRoutes() const { return numRoutes_; }
+size_t Individual::numRoutes() const { return routes_.size(); }
 
 size_t Individual::numClients() const { return numClients_; }
 
@@ -79,13 +75,13 @@ bool Individual::operator==(Individual const &other) const
     return distance_ == other.distance_
         && excessLoad_ == other.excessLoad_
         && timeWarp_ == other.timeWarp_
-        && numRoutes_ == other.numRoutes_
+        && routes_.size() == other.routes_.size()
         && neighbours == other.neighbours;
     // clang-format on
 }
 
 Individual::Individual(ProblemData const &data, XorShift128 &rng)
-    : routes_(data.numVehicles()), neighbours(data.numClients() + 1, {0, 0})
+    : neighbours(data.numClients() + 1, {0, 0})
 {
     // Shuffle clients (to create random routes)
     auto clients = std::vector<int>(data.numClients());
@@ -98,13 +94,15 @@ Individual::Individual(ProblemData const &data, XorShift128 &rng)
     auto const numClients = data.numClients();
     auto const perVehicle = std::max<size_t>(numClients / numVehicles, 1);
     auto const perRoute = perVehicle + (numClients % numVehicles != 0);
+    auto const numRoutes = (numClients + perRoute - 1) / perRoute;
 
-    std::vector<std::vector<Client>> routes(data.numVehicles());
+    std::vector<std::vector<Client>> routes(numRoutes);
     for (size_t idx = 0; idx != numClients; ++idx)
         routes[idx / perRoute].push_back(clients[idx]);
 
-    for (size_t idx = 0; idx != routes.size(); ++idx)
-        routes_[idx] = Route(data, routes[idx]);
+    routes_.reserve(numRoutes);
+    for (size_t idx = 0; idx != numRoutes; ++idx)
+        routes_.emplace_back(data, routes[idx]);
 
     makeNeighbours();
     evaluate(data);
@@ -112,7 +110,7 @@ Individual::Individual(ProblemData const &data, XorShift128 &rng)
 
 Individual::Individual(ProblemData const &data,
                        std::vector<std::vector<Client>> const &routes)
-    : routes_(data.numVehicles()), neighbours(data.numClients() + 1, {0, 0})
+    : neighbours(data.numClients() + 1, {0, 0})
 {
     if (routes.size() > data.numVehicles())
     {
@@ -142,14 +140,11 @@ Individual::Individual(ProblemData const &data,
         }
     }
 
+    // Only store non-empty routes
+    routes_.reserve(routes.size());
     for (size_t idx = 0; idx != routes.size(); ++idx)
-        routes_[idx] = Route(data, routes[idx]);
-
-    // a precedes b only when a is not empty and b is. Combined with a stable
-    // sort, this ensures we keep the original sorting as much as possible, but
-    // also make sure all empty routes are at the end of routes_.
-    auto comp = [](auto &a, auto &b) { return !a.empty() && b.empty(); };
-    std::stable_sort(routes_.begin(), routes_.end(), comp);
+        if (!routes[idx].empty())
+            routes_.emplace_back(data, routes[idx]);
 
     makeNeighbours();
     evaluate(data);
@@ -266,7 +261,7 @@ std::ostream &operator<<(std::ostream &out, Individual const &indiv)
 {
     auto const &routes = indiv.getRoutes();
 
-    for (size_t idx = 0; idx != indiv.numRoutes(); ++idx)
+    for (size_t idx = 0; idx != routes.size(); ++idx)
         out << "Route #" << idx + 1 << ": " << routes[idx] << '\n';
 
     out << "Distance: " << indiv.distance() << '\n';
diff --git a/pyvrp/cpp/Individual.h b/pyvrp/cpp/Individual.h
index 7ab86d5..d016f3a 100644
--- a/pyvrp/cpp/Individual.h
+++ b/pyvrp/cpp/Individual.h
@@ -68,7 +68,6 @@ public:
 private:
     using Routes = std::vector<Route>;
 
-    size_t numRoutes_ = 0;        // Number of routes
     size_t numClients_ = 0;       // Number of clients in the solution
     Distance distance_ = 0;       // Total distance
     Load excessLoad_ = 0;         // Total excess load over all routes
@@ -76,7 +75,7 @@ private:
     Cost uncollectedPrizes_ = 0;  // Total uncollected prize value
     Duration timeWarp_ = 0;       // Total time warp over all routes
 
-    Routes routes_;  // Routes - only the first numRoutes_ are non-empty
+    Routes routes_;  // Routes - only includes non-empty routes
     std::vector<std::pair<Client, Client>> neighbours;  // pairs of [pred, succ]
 
     // Determines the [pred, succ] pairs for each client.
@@ -87,9 +86,8 @@ private:
 
 public:
     /**
-     * Returns the number of non-empty routes in this individual's solution.
-     * Such non-empty routes are guaranteed to be in the lower indices of the
-     * routes returned by ``getRoutes``.
+     * Returns the number of (non-empty) routes in this individual's solution.
+     * Equal to the length of the vector of routes returned by ``getRoutes``.
      */
     [[nodiscard]] size_t numRoutes() const;
 
@@ -188,7 +186,7 @@ template <> struct hash<Individual>
     size_t operator()(Individual const &individual) const
     {
         size_t res = 17;
-        res = res * 31 + std::hash<size_t>()(individual.numRoutes_);
+        res = res * 31 + std::hash<size_t>()(individual.routes_.size());
         res = res * 31 + std::hash<Distance>()(individual.distance_);
         res = res * 31 + std::hash<Load>()(individual.excessLoad_);
         res = res * 31 + std::hash<Duration>()(individual.timeWarp_);
diff --git a/pyvrp/cpp/crossover/crossover.cpp b/pyvrp/cpp/crossover/crossover.cpp
index c398fea..af806f2 100644
--- a/pyvrp/cpp/crossover/crossover.cpp
+++ b/pyvrp/cpp/crossover/crossover.cpp
@@ -63,14 +63,11 @@ void crossover::greedyRepair(Routes &routes,
                              ProblemData const &data,
                              CostEvaluator const &costEvaluator)
 {
-    // Determine the index just past the last non-empty route.
-    auto const pred = [](Route const &route) { return !route.empty(); };
-    auto const last = std::find_if(routes.rbegin(), routes.rend(), pred);
-    auto const numRoutes = std::distance(last, routes.rend());
+    auto const numRoutes = routes.size();
 
     // Determine centroids of each route.
     std::vector<std::pair<double, double>> centroids(numRoutes, {0, 0});
-    for (auto rIdx = 0; rIdx != numRoutes; ++rIdx)
+    for (size_t rIdx = 0; rIdx != numRoutes; ++rIdx)
         for (Client client : routes[rIdx])
         {
             auto const size = static_cast<double>(routes[rIdx].size());
@@ -89,7 +86,7 @@ void crossover::greedyRepair(Routes &routes,
         // Determine non-empty route with centroid nearest to this client.
         auto bestDistance = std::numeric_limits<double>::max();
         auto bestRouteIdx = 0;
-        for (auto rIdx = 0; rIdx != numRoutes; ++rIdx)
+        for (size_t rIdx = 0; rIdx != numRoutes; ++rIdx)
         {
             if (routes[rIdx].empty())
                 continue;
diff --git a/pyvrp/cpp/crossover/crossover.h b/pyvrp/cpp/crossover/crossover.h
index 1dd5a4b..6aa7ce4 100644
--- a/pyvrp/cpp/crossover/crossover.h
+++ b/pyvrp/cpp/crossover/crossover.h
@@ -12,8 +12,8 @@
 namespace crossover
 {
 /**
- * Greedily inserts each unplanned client into the non-empty route that's
- * nearest to the client.
+ * Greedily inserts each unplanned client into the route that's nearest to the
+ * client.
  */
 void greedyRepair(std::vector<std::vector<int>> &routes,
                   std::vector<int> const &unplanned,
diff --git a/pyvrp/cpp/crossover/selective_route_exchange.cpp b/pyvrp/cpp/crossover/selective_route_exchange.cpp
index 063c226..b0f580f 100644
--- a/pyvrp/cpp/crossover/selective_route_exchange.cpp
+++ b/pyvrp/cpp/crossover/selective_route_exchange.cpp
@@ -15,8 +15,6 @@ namespace
 // Angle of the given route w.r.t. the centroid of all client locations.
 double routeAngle(ProblemData const &data, Route const &route)
 {
-    assert(!route.empty());
-
     // This computes a pseudo-angle that sorts roughly equivalently to the atan2
     // angle, but is much faster to compute. See the following post for details:
     // https://stackoverflow.com/a/16561333/4316405.
@@ -30,9 +28,7 @@ double routeAngle(ProblemData const &data, Route const &route)
 Routes sortByAscAngle(ProblemData const &data, Routes routes)
 {
     auto cmp = [&data](Route a, Route b) {
-        return (a.empty() || b.empty())
-                   ? !a.empty() && b.empty()
-                   : routeAngle(data, a) < routeAngle(data, b);
+        return routeAngle(data, a) < routeAngle(data, b);
     };
 
     std::sort(routes.begin(), routes.end(), cmp);
diff --git a/pyvrp/cpp/educate/LocalSearch.cpp b/pyvrp/cpp/educate/LocalSearch.cpp
index c6d80ab..dfc45b8 100644
--- a/pyvrp/cpp/educate/LocalSearch.cpp
+++ b/pyvrp/cpp/educate/LocalSearch.cpp
@@ -16,10 +16,6 @@ Individual LocalSearch::search(Individual &individual,
 {
     loadIndividual(individual);
 
-    // Shuffling the order beforehand adds diversity to the search
-    std::shuffle(orderNodes.begin(), orderNodes.end(), rng);
-    std::shuffle(nodeOps.begin(), nodeOps.end(), rng);
-
     if (nodeOps.empty())
         throw std::runtime_error("No known node operators.");
 
@@ -97,10 +93,6 @@ Individual LocalSearch::intensify(Individual &individual,
 
     auto const overlapTolerance = overlapToleranceDegrees * 65536;
 
-    // Shuffling the order beforehand adds diversity to the search
-    std::shuffle(orderRoutes.begin(), orderRoutes.end(), rng);
-    std::shuffle(routeOps.begin(), routeOps.end(), rng);
-
     if (routeOps.empty())
         throw std::runtime_error("No known route operators.");
 
@@ -146,6 +138,15 @@ Individual LocalSearch::intensify(Individual &individual,
     return exportIndividual();
 }
 
+void LocalSearch::shuffle(XorShift128 &rng)
+{
+    std::shuffle(orderNodes.begin(), orderNodes.end(), rng);
+    std::shuffle(nodeOps.begin(), nodeOps.end(), rng);
+
+    std::shuffle(orderRoutes.begin(), orderRoutes.end(), rng);
+    std::shuffle(routeOps.begin(), routeOps.end(), rng);
+}
+
 bool LocalSearch::applyNodeOps(Node *U,
                                Node *V,
                                CostEvaluator const &costEvaluator)
@@ -303,7 +304,7 @@ void LocalSearch::loadIndividual(Individual const &individual)
 
         Route *route = &routes[r];
 
-        if (!routesIndiv[r].empty())
+        if (r < routesIndiv.size())
         {
             Node *client = &clients[routesIndiv[r][0]];
             client->route = route;
@@ -383,13 +384,13 @@ void LocalSearch::setNeighbours(Neighbours neighbours)
     this->neighbours = neighbours;
 }
 
-LocalSearch::Neighbours LocalSearch::getNeighbours() { return neighbours; }
+LocalSearch::Neighbours const &LocalSearch::getNeighbours() const
+{
+    return neighbours;
+}
 
-LocalSearch::LocalSearch(ProblemData &data,
-                         XorShift128 &rng,
-                         Neighbours neighbours)
+LocalSearch::LocalSearch(ProblemData const &data, Neighbours neighbours)
     : data(data),
-      rng(rng),
       neighbours(data.numClients() + 1),
       orderNodes(data.numClients()),
       orderRoutes(data.numVehicles()),
diff --git a/pyvrp/cpp/educate/LocalSearch.h b/pyvrp/cpp/educate/LocalSearch.h
index d65c6d9..6ae7599 100644
--- a/pyvrp/cpp/educate/LocalSearch.h
+++ b/pyvrp/cpp/educate/LocalSearch.h
@@ -19,15 +19,14 @@ class LocalSearch
     using RouteOp = LocalSearchOperator<Route>;
     using Neighbours = std::vector<std::vector<int>>;
 
-    ProblemData &data;
-    XorShift128 &rng;
+    ProblemData const &data;
 
-    // Neighborhood restrictions: For each client, list of nearby clients (size
-    // nbClients + 1, but nothing stored for the depot!)
+    // Neighborhood restrictions: list of nearby clients for each client (size
+    // numClients + 1, but nothing stored for the depot!)
     Neighbours neighbours;
 
-    std::vector<int> orderNodes;   // random node order used in RI operators
-    std::vector<int> orderRoutes;  // random route order used in SWAP* operators
+    std::vector<int> orderNodes;   // node order used by LocalSearch::search
+    std::vector<int> orderRoutes;  // route order used by LocalSearch::intensify
 
     std::vector<int> lastModified;  // tracks when routes were last modified
 
@@ -85,7 +84,7 @@ public:
     /**
      * @return The neighbourhood structure currently in use.
      */
-    Neighbours getNeighbours();
+    Neighbours const &getNeighbours() const;
 
     /**
      * Performs regular (node-based) local search around the given individual,
@@ -95,15 +94,20 @@ public:
                       CostEvaluator const &costEvaluator);
 
     /**
-     * Performs a more intensive local search around the given individual,
-     * using route-based operators and subpath enumeration. Returns a new,
-     * hopefully improved individual.
+     * Performs a more intensive route-based local search around the given
+     * individual, and returns a new, hopefully improved individual.
      */
     Individual intensify(Individual &individual,
                          CostEvaluator const &costEvaluator,
                          int overlapToleranceDegrees = 0);
 
-    LocalSearch(ProblemData &data, XorShift128 &rng, Neighbours neighbours);
+    /**
+     * Shuffles the order in which the node and route pairs are evaluated, and
+     * the order in which node and route operators are applied.
+     */
+    void shuffle(XorShift128 &rng);
+
+    LocalSearch(ProblemData const &data, Neighbours neighbours);
 };
 
 #endif  // PYVRP_LOCALSEARCH_H
diff --git a/pyvrp/cpp/educate/LocalSearch_bindings.cpp b/pyvrp/cpp/educate/LocalSearch_bindings.cpp
index 0c69c3f..56aa782 100644
--- a/pyvrp/cpp/educate/LocalSearch_bindings.cpp
+++ b/pyvrp/cpp/educate/LocalSearch_bindings.cpp
@@ -9,14 +9,10 @@ namespace py = pybind11;
 PYBIND11_MODULE(_LocalSearch, m)
 {
     py::class_<LocalSearch>(m, "LocalSearch")
-        .def(py::init<ProblemData &,
-                      XorShift128 &,
-                      std::vector<std::vector<int>>>(),
+        .def(py::init<ProblemData const &, std::vector<std::vector<int>>>(),
              py::arg("data"),
-             py::arg("rng"),
              py::arg("neighbours"),
-             py::keep_alive<1, 2>(),  // keep data and rng alive at least until
-             py::keep_alive<1, 3>())  // local search is freed
+             py::keep_alive<1, 2>())  // keep data alive until LS is freed
         .def("add_node_operator",
              &LocalSearch::addNodeOperator,
              py::arg("op"),
@@ -39,5 +35,6 @@ PYBIND11_MODULE(_LocalSearch, m)
              &LocalSearch::intensify,
              py::arg("individual"),
              py::arg("cost_evaluator"),
-             py::arg("overlap_tolerance_degrees") = 0);
+             py::arg("overlap_tolerance_degrees") = 0)
+        .def("shuffle", &LocalSearch::shuffle, py::arg("rng"));
 }
diff --git a/pyvrp/educate/LocalSearch.py b/pyvrp/educate/LocalSearch.py
index a4631df..f7328ae 100644
--- a/pyvrp/educate/LocalSearch.py
+++ b/pyvrp/educate/LocalSearch.py
@@ -27,12 +27,10 @@ class LocalSearch:
     """
 
     def __init__(
-        self,
-        data: ProblemData,
-        rng: XorShift128,
-        neighbours: Neighbours,
+        self, data: ProblemData, rng: XorShift128, neighbours: Neighbours
     ):
-        self._ls = _LocalSearch(data, rng, neighbours)
+        self._ls = _LocalSearch(data, neighbours)
+        self._rng = rng
 
     def add_node_operator(self, op):
         """
@@ -170,6 +168,7 @@ class LocalSearch:
             The improved individual. This is not the same object as the
             individual that was passed in.
         """
+        self._ls.shuffle(self._rng)
         return self._ls.intensify(
             individual, cost_evaluator, overlap_tolerance_degrees
         )
@@ -200,4 +199,5 @@ class LocalSearch:
             The improved individual. This is not the same object as the
             individual that was passed in.
         """
+        self._ls.shuffle(self._rng)
         return self._ls.search(individual, cost_evaluator)
diff --git a/pyvrp/educate/_LocalSearch.pyi b/pyvrp/educate/_LocalSearch.pyi
index c30c231..ff7018f 100644
--- a/pyvrp/educate/_LocalSearch.pyi
+++ b/pyvrp/educate/_LocalSearch.pyi
@@ -11,13 +11,13 @@ class LocalSearch:
     def __init__(
         self,
         data: ProblemData,
-        rng: XorShift128,
         neighbours: Neighbours,
     ) -> None: ...
     def add_node_operator(self, op) -> None: ...
     def add_route_operator(self, op) -> None: ...
     def set_neighbours(self, neighbours: Neighbours) -> None: ...
     def get_neighbours(self) -> Neighbours: ...
+    def shuffle(self, rng: XorShift128) -> None: ...
     def intensify(
         self,
         individual: Individual,
diff --git a/pyvrp/exceptions.py b/pyvrp/exceptions.py
index 1c23845..1ff5b42 100644
--- a/pyvrp/exceptions.py
+++ b/pyvrp/exceptions.py
@@ -4,8 +4,6 @@ class StatisticsNotCollectedError(Exception):
     that can only be used when statistics are available.
     """
 
-    pass
-
 
 class EmptySolutionWarning(UserWarning):
     """
@@ -13,4 +11,9 @@ class EmptySolutionWarning(UserWarning):
     forbidden, per se, but very odd.
     """
 
-    pass
+
+class ScalingWarning(UserWarning):
+    """
+    Raised when the distance or duration values in the problem are very large,
+    which could cause the algorithm to start using forbidden edges as well.
+    """
diff --git a/pyvrp/plotting/plot_solution.py b/pyvrp/plotting/plot_solution.py
index 4c392c0..9300a7f 100644
--- a/pyvrp/plotting/plot_solution.py
+++ b/pyvrp/plotting/plot_solution.py
@@ -39,19 +39,18 @@ def plot_solution(
     ax.scatter(x_coords[0], y_coords[0], label="Depot", **kwargs)
 
     for idx, route in enumerate(solution.get_routes(), 1):
-        if route:
-            x = x_coords[route]
-            y = y_coords[route]
-
-            # Coordinates of customers served by this route.
-            if len(route) == 1 or plot_customers:
-                ax.scatter(x, y, label=f"Route {idx}", zorder=3, s=75)
-            ax.plot(x, y)
-
-            # Edges from and to the depot, very thinly dashed.
-            kwargs = dict(ls=(0, (5, 15)), linewidth=0.25, color="grey")
-            ax.plot([x_coords[0], x[0]], [y_coords[0], y[0]], **kwargs)
-            ax.plot([x[-1], x_coords[0]], [y[-1], y_coords[0]], **kwargs)
+        x = x_coords[route]
+        y = y_coords[route]
+
+        # Coordinates of customers served by this route.
+        if len(route) == 1 or plot_customers:
+            ax.scatter(x, y, label=f"Route {idx}", zorder=3, s=75)
+        ax.plot(x, y)
+
+        # Edges from and to the depot, very thinly dashed.
+        kwargs = dict(ls=(0, (5, 15)), linewidth=0.25, color="grey")
+        ax.plot([x_coords[0], x[0]], [y_coords[0], y[0]], **kwargs)
+        ax.plot([x[-1], x_coords[0]], [y[-1], y_coords[0]], **kwargs)
 
     ax.grid(color="grey", linestyle="solid", linewidth=0.2)
 
diff --git a/pyvrp/read.py b/pyvrp/read.py
index a182943..b255cd7 100644
--- a/pyvrp/read.py
+++ b/pyvrp/read.py
@@ -2,10 +2,14 @@ import functools
 import pathlib
 from numbers import Number
 from typing import Callable, Dict, List, Union
+from warnings import warn
 
 import numpy as np
 import vrplib
 
+from pyvrp.constants import MAX_USER_VALUE
+from pyvrp.exceptions import ScalingWarning
+
 from ._ProblemData import Client, ProblemData
 
 _Routes = List[List[int]]
@@ -165,6 +169,13 @@ def read(
         for idx in range(dimension)
     ]
 
+    if max(distances.max(), durations.max()) > MAX_USER_VALUE:
+        msg = """
+        The maximum distance or duration value is very large. This might
+        impact numerical stability. Consider rescaling your input data.
+        """
+        warn(msg, ScalingWarning)
+
     return ProblemData(
         clients,
         num_vehicles,


Test patch:
diff --git a/pyvrp/educate/tests/test_LocalSearch.py b/pyvrp/educate/tests/test_LocalSearch.py
index 2b678ee..9813e99 100644
--- a/pyvrp/educate/tests/test_LocalSearch.py
+++ b/pyvrp/educate/tests/test_LocalSearch.py
@@ -10,6 +10,7 @@ from pyvrp.educate import (
     Neighbours,
     compute_neighbours,
 )
+from pyvrp.educate._LocalSearch import LocalSearch as cpp_LocalSearch
 from pyvrp.tests.helpers import read
 
 
@@ -94,9 +95,7 @@ def test_local_search_set_get_neighbours(
     symmetric_neighbours: bool,
 ):
     data = read("data/RC208.txt", "solomon", round_func="trunc")
-
-    seed = 42
-    rng = XorShift128(seed=seed)
+    rng = XorShift128(seed=42)
 
     params = NeighbourhoodParams(nb_granular=1)
     prev_neighbours = compute_neighbours(data, params)
@@ -188,3 +187,27 @@ def test_prize_collecting():
     # ...but an optimised prize-collecting solution is likely not complete.
     assert_(improved.num_clients() < individual.num_clients())
     assert_(improved_cost < individual_cost)
+
+
+def test_cpp_shuffle_results_in_different_solution():
+    data = read("data/RC208.txt", "solomon", round_func="trunc")
+    rng = XorShift128(seed=42)
+
+    ls = cpp_LocalSearch(data, compute_neighbours(data))
+    ls.add_node_operator(Exchange10(data))
+    ls.add_node_operator(Exchange11(data))
+
+    cost_evaluator = CostEvaluator(1, 1)
+    individual = Individual.make_random(data, rng)
+
+    # LocalSearch::search is deterministic, so two calls with the same base
+    # individual should result in the same improved solution.
+    improved1 = ls.search(individual, cost_evaluator)
+    improved2 = ls.search(individual, cost_evaluator)
+    assert_(improved1 == improved2)
+
+    # But the shuffle method changes the order in which moves are evaluated,
+    # which should result in a very different search trajectory.
+    ls.shuffle(rng)
+    improved3 = ls.search(individual, cost_evaluator)
+    assert_(improved3 != improved1)
diff --git a/pyvrp/tests/data/ReallyLargeDistance.txt b/pyvrp/tests/data/ReallyLargeDistance.txt
new file mode 100644
index 0000000..f8880a3
--- /dev/null
+++ b/pyvrp/tests/data/ReallyLargeDistance.txt
@@ -0,0 +1,41 @@
+NAME : TEST
+COMMENT : Small problem instance with four customers
+DIMENSION : 5
+EDGE_WEIGHT_TYPE : EXPLICIT
+EDGE_WEIGHT_FORMAT : FULL_MATRIX
+CAPACITY : 10
+VEHICLES : 3
+NODE_COORD_SECTION
+1    2334   726
+2    226    1297
+3    590    530
+4    435    718
+5    1191   639
+EDGE_WEIGHT_SECTION
+0     1544  1944  1931  1000000000
+1726     0  1992  1427  1593
+1965  1975     0  621   1090
+2063  1433   647  0     818
+1475  1594  1090  828   0
+DEMAND_SECTION
+1    0
+2    5
+3    5
+4    3
+5    5
+TIME_WINDOW_SECTION
+1    0      45000
+2    15600  22500
+3    12000  19500
+4    8400   15300
+5    12000  19500
+SERVICE_TIME_SECTION
+1    0
+2    360
+3    360
+4    420
+5    360
+DEPOT_SECTION
+1
+-1
+EOF
diff --git a/pyvrp/tests/test_Individual.py b/pyvrp/tests/test_Individual.py
index 8279b9b..f094761 100644
--- a/pyvrp/tests/test_Individual.py
+++ b/pyvrp/tests/test_Individual.py
@@ -8,22 +8,19 @@ from pyvrp import Client, Individual, ProblemData, Route, XorShift128
 from pyvrp.tests.helpers import read
 
 
-def test_route_constructor_sorts_by_empty():
+def test_route_constructor_filters_empty():
     data = read("data/OkSmall.txt")
 
     indiv = Individual(data, [[3, 4], [], [1, 2]])
     routes = indiv.get_routes()
 
-    # num_routes() should show two non-empty routes. However, we passed in
-    # three routes, so len(routes) should not have changed.
+    # num_routes() and len(routes) should show two non-empty routes.
     assert_equal(indiv.num_routes(), 2)
-    assert_equal(len(routes), 3)
+    assert_equal(len(routes), 2)
 
-    # We expect Individual to sort the routes such that all non-empty routes
-    # are in the lower indices.
+    # The only two non-empty routes should now each have two clients.
     assert_equal(len(routes[0]), 2)
     assert_equal(len(routes[1]), 2)
-    assert_equal(len(routes[2]), 0)
 
 
 def test_random_constructor_cycles_over_routes():
@@ -37,9 +34,9 @@ def test_random_constructor_cycles_over_routes():
     routes = indiv.get_routes()
 
     assert_equal(indiv.num_routes(), 2)
-    assert_equal(len(routes), 3)
+    assert_equal(len(routes), 2)
 
-    for idx, size in enumerate([2, 2, 0]):
+    for idx, size in enumerate([2, 2]):
         assert_equal(len(routes[idx]), size)
 
 
@@ -48,10 +45,9 @@ def test_route_constructor_raises_too_many_vehicles():
 
     assert_equal(data.num_vehicles, 3)
 
-    # Only two routes should not raise. But we should always get num_vehicles
-    # routes back.
+    # Only two routes should not raise.
     individual = Individual(data, [[1, 2], [4, 3]])
-    assert_equal(len(individual.get_routes()), data.num_vehicles)
+    assert_equal(len(individual.get_routes()), 2)
 
     # Empty third route should not raise.
     Individual(data, [[1, 2], [4, 3], []])
@@ -314,16 +310,15 @@ def test_str_contains_essential_information():
         str_representation = str(individual).splitlines()
 
         routes = individual.get_routes()
-        num_routes = individual.num_routes()
 
-        # There should be no more than num_routes lines (each detailing a
+        # There should be no more than len(routes) lines (each detailing a
         # single route), and two final lines containing distance and prizes.
-        assert_equal(len(str_representation), num_routes + 2)
+        assert_equal(len(str_representation), len(routes) + 2)
 
-        # The first num_routes lines should each contain a route, where each
+        # The first len(routes) lines should each contain a route, where each
         # route should contain every client that is in the route as returned
         # by get_routes().
-        for route, str_route in zip(routes[:num_routes], str_representation):
+        for route, str_route in zip(routes, str_representation):
             for client in route:
                 assert_(str(client) in str_route)
 
diff --git a/pyvrp/tests/test_Model.py b/pyvrp/tests/test_Model.py
new file mode 100644
index 0000000..2de7419
--- /dev/null
+++ b/pyvrp/tests/test_Model.py
@@ -0,0 +1,242 @@
+from numpy.testing import assert_, assert_equal, assert_raises, assert_warns
+from pytest import mark
+
+from pyvrp import Model
+from pyvrp.constants import MAX_USER_VALUE, MAX_VALUE
+from pyvrp.exceptions import ScalingWarning
+from pyvrp.stop import MaxIterations
+from pyvrp.tests.helpers import read
+
+
+def test_model_data():
+    model = Model()
+
+    # Let's add some data: a single client, and edges from/to the depot.
+    depot = model.add_depot(0, 0)
+    client = model.add_client(0, 1, demand=1)
+    model.add_edge(depot, client, 1, 1)
+    model.add_edge(client, depot, 1, 1)
+    model.add_vehicle_type(1, 1)
+
+    # The model should now have one client.
+    data = model.data()
+    assert_equal(data.num_clients, 1)
+
+
+def test_add_depot_raises_more_than_one_depot():
+    model = Model()
+    model.add_depot(0, 0)  # first depot should be OK
+
+    with assert_raises(ValueError):
+        model.add_depot(0, 1)  # second (and more) should not be
+
+
+def test_add_edge_raises_negative_distance_or_duration():
+    model = Model()
+    depot = model.add_depot(0, 0)
+    client = model.add_client(0, 1)
+
+    model.add_edge(depot, client, distance=0, duration=0)  # zero should be OK
+
+    with assert_raises(ValueError):  # negative distance should not be OK
+        model.add_edge(client, depot, distance=-1, duration=0)
+
+    with assert_raises(ValueError):  # negative duration should also not be OK
+        model.add_edge(client, depot, distance=0, duration=-1)
+
+
+def test_add_vehicle_type_raises_more_than_one_type():
+    model = Model()
+    model.add_vehicle_type(1, 10)  # first type should be OK
+
+    with assert_raises(ValueError):
+        model.add_vehicle_type(2, 10)  # second (and more) should not be
+
+
+@mark.parametrize(
+    "number, capacity",
+    [
+        (0, 1),  # zero vehicles is not OK (but zero capacity is)
+        (-1, 1),  # negative vehicles is not OK
+        (1, -1),  # negative capacity is not OK
+    ],
+)
+def test_add_vehicle_type_raises_negative_number_or_capacity(number, capacity):
+    model = Model()
+
+    with assert_raises(ValueError):
+        model.add_vehicle_type(number, capacity)
+
+
+def test_add_client_attributes():
+    model = Model()
+    client = model.add_client(
+        x=1,
+        y=2,
+        demand=3,
+        service_duration=4,
+        tw_early=5,
+        tw_late=6,
+        prize=7,
+        required=False,
+    )
+
+    assert_equal(client.x, 1)
+    assert_equal(client.y, 2)
+    assert_equal(client.demand, 3)
+    assert_equal(client.service_duration, 4)
+    assert_equal(client.tw_early, 5)
+    assert_equal(client.tw_late, 6)
+    assert_equal(client.prize, 7)
+    assert_(not client.required)
+
+
+def test_add_depot_attributes():
+    model = Model()
+    depot = model.add_depot(x=1, y=0, tw_early=3, tw_late=5)
+
+    assert_equal(depot.x, 1)
+    assert_equal(depot.y, 0)
+    assert_equal(depot.tw_early, 3)
+    assert_equal(depot.tw_late, 5)
+
+
+def test_add_edge():
+    model = Model()
+    depot = model.add_depot(0, 0)
+    client = model.add_client(0, 1)
+    edge = model.add_edge(depot, client, distance=15, duration=49)
+
+    assert_(edge.frm is depot)
+    assert_(edge.to is client)
+    assert_equal(edge.distance, 15)
+    assert_equal(edge.duration, 49)
+
+
+def test_add_vehicle_type():
+    model = Model()
+    vehicle_type = model.add_vehicle_type(number=10, capacity=998)
+
+    assert_equal(vehicle_type.number, 10)
+    assert_equal(vehicle_type.capacity, 998)
+
+
+def test_from_data():
+    read_data = read("data/E-n22-k4.txt", round_func="dimacs")
+    model = Model.from_data(read_data)
+    model_data = model.data()
+
+    # We can first check if the overall problem dimension numbers agree.
+    assert_equal(model_data.num_clients, read_data.num_clients)
+    assert_equal(model_data.num_vehicles, read_data.num_vehicles)
+    assert_equal(model_data.vehicle_capacity, read_data.vehicle_capacity)
+
+    # It's a bit cumbersome to compare the whole matrices, so we use a few
+    # sample traces from the distance and duration matrices instead.
+    assert_equal(model_data.dist(3, 4), read_data.dist(3, 4))
+    assert_equal(model_data.duration(2, 1), read_data.duration(2, 1))
+
+
+def test_from_data_and_solve():
+    # Solve the small E-n22-k4 instance using the from_data constructor.
+    data = read("data/E-n22-k4.txt", round_func="dimacs")
+    model = Model.from_data(data)
+    res = model.solve(stop=MaxIterations(100), seed=0)
+    assert_equal(res.cost(), 3_743)
+    assert_(res.is_feasible())
+
+    data = read("data/OkSmall.txt")
+    model = Model.from_data(data)
+    res = model.solve(stop=MaxIterations(100), seed=0)
+    assert_equal(res.cost(), 9_155)
+    assert_(res.is_feasible())
+
+
+def test_model_and_solve():
+    # Solve the small OkSmall instance using the from_data constructor.
+    data = read("data/OkSmall.txt")
+    model = Model.from_data(data)
+    res = model.solve(stop=MaxIterations(100), seed=0)
+    assert_equal(res.cost(), 9_155)
+    assert_(res.is_feasible())
+
+    # Now do the same thing, but model the instance using the modelling API.
+    # This should of course result in the same solution.
+    model = Model()
+    model.add_vehicle_type(number=3, capacity=10)
+    depot = model.add_depot(x=2334, y=726, tw_early=0, tw_late=45000)
+    clients = [
+        model.add_client(226, 1297, 5, 360, 15600, 22500),
+        model.add_client(590, 530, 5, 360, 12000, 19500),
+        model.add_client(435, 718, 3, 420, 8400, 15300),
+        model.add_client(1191, 639, 5, 360, 12000, 19500),
+    ]
+
+    edge_weights = [
+        [0, 1544, 1944, 1931, 1476],
+        [1726, 0, 1992, 1427, 1593],
+        [1965, 1975, 0, 621, 1090],
+        [2063, 1433, 647, 0, 818],
+        [1475, 1594, 1090, 828, 0],
+    ]
+
+    for idx, client in enumerate(clients, 1):
+        from_depot = edge_weights[0][idx]
+        to_depot = edge_weights[idx][0]
+
+        model.add_edge(depot, client, from_depot, from_depot)
+        model.add_edge(client, depot, to_depot, to_depot)
+
+        for other_idx, other in enumerate(clients, 1):
+            from_client = edge_weights[idx][other_idx]
+            to_client = edge_weights[other_idx][idx]
+
+            model.add_edge(client, other, from_client, from_client)
+            model.add_edge(other, client, to_client, to_client)
+
+    res = model.solve(stop=MaxIterations(100), seed=0)
+    assert_equal(res.cost(), 9_155)
+    assert_(res.is_feasible())
+
+
+def test_partial_distance_duration_matrix():
+    model = Model()
+    depot = model.add_depot(0, 0)
+    clients = [model.add_client(0, 1), model.add_client(1, 0)]
+
+    # Only a subset of all possible edges exist. The model should be able to
+    # handle that.
+    model.add_edge(depot, clients[0], distance=1)
+    model.add_edge(clients[0], clients[1], distance=2)
+    model.add_edge(clients[1], depot, distance=1)
+
+    model.add_vehicle_type(number=1, capacity=0)
+
+    # These edges were not set, so their distance values should default to the
+    # maximum value we use for such edges.
+    data = model.data()
+    assert_equal(data.dist(0, 2), MAX_VALUE)
+    assert_equal(data.dist(1, 0), MAX_VALUE)
+
+    res = model.solve(stop=MaxIterations(100), seed=4)
+    assert_equal(res.best.num_routes(), 1)
+    assert_equal(res.cost(), 4)  # depot -> client 1 -> client 2 -> depot
+    assert_(res.is_feasible())
+
+
+def test_data_warns_about_scaling_issues(recwarn):
+    model = Model()
+    model.add_vehicle_type(number=1, capacity=0)
+    depot = model.add_depot(0, 0)
+    client = model.add_client(1, 1)
+
+    # Normal distance sizes; should all be OK.
+    for distance in [1, 10, 100, 1_000, 10_000, 100_000, MAX_USER_VALUE]:
+        model.add_edge(client, depot, distance=distance)
+        model.data()
+        assert_equal(len(recwarn), 0)
+
+    # But a value exceeding the maximum user value is not OK. This should warn.
+    model.add_edge(depot, client, distance=MAX_USER_VALUE + 1)
+    with assert_warns(ScalingWarning):
+        model.data()
diff --git a/pyvrp/tests/test_Result.py b/pyvrp/tests/test_Result.py
index f468f0d..e024cba 100644
--- a/pyvrp/tests/test_Result.py
+++ b/pyvrp/tests/test_Result.py
@@ -81,7 +81,6 @@ def test_str_contains_essential_information(routes):
     else:
         assert_("INFEASIBLE" in str_representation)
 
-    # And make sure that all non-empty routes are printed as well.
+    # And make sure that all (non-empty) routes are printed as well.
     for route in individual.get_routes():
-        if route:
-            assert_(str(route) in str_representation)
+        assert_(str(route) in str_representation)
diff --git a/pyvrp/tests/test_read.py b/pyvrp/tests/test_read.py
index 2d555e6..367affe 100644
--- a/pyvrp/tests/test_read.py
+++ b/pyvrp/tests/test_read.py
@@ -1,6 +1,7 @@
-from numpy.testing import assert_equal, assert_raises
+from numpy.testing import assert_equal, assert_raises, assert_warns
 from pytest import mark
 
+from pyvrp.exceptions import ScalingWarning
 from pyvrp.tests.helpers import read
 
 
@@ -166,4 +167,13 @@ def test_reading_RC208_instance():  # Solomon style instance
         assert_equal(data.client(client).service_duration, 100)
 
 
+def test_warns_about_scaling_issues(recwarn):
+    # But 100 million is too large. That's really close to INT_MAX, and might
+    # result in issues. So now a warning should be issued.
+    with assert_warns(ScalingWarning):
+        # The arc from the depot to client 4 is really large (one billion), so
+        # that should trigger a warning.
+        read("data/ReallyLargeDistance.txt")
+
+
 # TODO test round funcs