Repository: PyVRP/PyVRP
Pull request #349 (2023-08-17 14:17:19)

Test descriptions
This PR adds a bundle of docstrings to our tests, removes duplicate tests, improves existing ones, and adds a few more functional tests in places that had open TODOs.

<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>


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PyVRP/PyVRP/pull/349?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.02%`** :tada:
> Comparison is base [(`7f75c2e`)](https://app.codecov.io/gh/PyVRP/PyVRP/commit/7f75c2eabdc13ccb7480d80e2a1755dc71aa011a?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) 95.26% compared to head [(`eda3def`)](https://app.codecov.io/gh/PyVRP/PyVRP/pull/349?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) 95.29%.

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


```diff
@@            Coverage Diff             @@
##             main     #349      +/-   ##
==========================================
+ Coverage   95.26%   95.29%   +0.02%     
==========================================
  Files          64       64              
  Lines        2493     2486       -7     
  Branches      217      217              
==========================================
- Hits         2375     2369       -6     
+ Misses         23       22       -1     
  Partials       95       95              
```


| [Files Changed](https://app.codecov.io/gh/PyVRP/PyVRP/pull/349?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP) | Coverage Δ | |
|---|---|---|
| [pyvrp/Model.py](https://app.codecov.io/gh/PyVRP/PyVRP/pull/349?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP#diff-cHl2cnAvTW9kZWwucHk=) | `100.00% <100.00%> (ø)` | |

... and [5 files with indirect coverage changes](https://app.codecov.io/gh/PyVRP/PyVRP/pull/349/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/349?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PyVRP).


Patch:
diff --git a/pyvrp/Model.py b/pyvrp/Model.py
index efb52c2..17af3e9 100644
--- a/pyvrp/Model.py
+++ b/pyvrp/Model.py
@@ -181,6 +181,13 @@ class Model:
         if distance < 0 or duration < 0:
             raise ValueError("Cannot have negative edge distance or duration.")
 
+        if max(distance, duration) > MAX_USER_VALUE:
+            msg = """
+            The given distance or duration value is very large. This may impact
+            numerical stability. Consider rescaling your input data.
+            """
+            warn(msg, ScalingWarning)
+
         edge = Edge(frm, to, distance, duration)
         self._edges.append(edge)
         return edge
@@ -216,16 +223,6 @@ class Model:
         locs = self.locations
         loc2idx = {id(loc): idx for idx, loc in enumerate(locs)}
 
-        if self._edges:
-            max_value = max(max(e.distance, e.duration) for e in self._edges)
-
-            if max_value > MAX_USER_VALUE:
-                msg = """
-                The maximum distance or duration value is very large. This may
-                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)


Test patch:
diff --git a/pyvrp/crossover/tests/test_selective_route_exchange.py b/pyvrp/crossover/tests/test_selective_route_exchange.py
index 5f7b9f6..6190270 100644
--- a/pyvrp/crossover/tests/test_selective_route_exchange.py
+++ b/pyvrp/crossover/tests/test_selective_route_exchange.py
@@ -31,6 +31,11 @@ def test_same_parents_same_offspring():
 
 
 def test_srex_empty_solution():
+    """
+    Tests that SREX returns the other parent when one of the solutions is
+    empty. This can occur during prize collecting, and in that case there are
+    no routes to exchange.
+    """
     data = read("data/p06-2-50.vrp", round_func="dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -59,6 +64,9 @@ def test_srex_empty_solution():
     ],
 )
 def test_raise_invalid_arguments(idx1, idx2, num_moved_routes):
+    """
+    Tests that SREX raises when some of the arguments are invalid.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
 
diff --git a/pyvrp/diversity/tests/test_broken_pairs_distance.py b/pyvrp/diversity/tests/test_broken_pairs_distance.py
index 4e06ea8..b6b47f1 100644
--- a/pyvrp/diversity/tests/test_broken_pairs_distance.py
+++ b/pyvrp/diversity/tests/test_broken_pairs_distance.py
@@ -6,6 +6,10 @@ from pyvrp.tests.helpers import read
 
 
 def test_broken_pairs_distance():
+    """
+    Test broken pairs distance (BPD) calculations on a number of solutions to
+    a small instance.
+    """
     data = read("data/OkSmall.txt")
 
     sol1 = Solution(data, [[1, 2, 3, 4]])
diff --git a/pyvrp/plotting/tests/test_plotting.py b/pyvrp/plotting/tests/test_plotting.py
index 2c3e2fc..c624cfe 100644
--- a/pyvrp/plotting/tests/test_plotting.py
+++ b/pyvrp/plotting/tests/test_plotting.py
@@ -18,6 +18,9 @@ IMG_KWARGS = dict(remove_text=True, tol=8, extensions=["png"], style="mpl20")
 
 @img_comp(["plot_solution", "plot_solution_with_customers"], **IMG_KWARGS)
 def test_plot_solution():
+    """
+    Compares plots generated by ``plot_solution``.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     bks = read_solution("data/RC208.sol")
 
@@ -29,6 +32,9 @@ def test_plot_solution():
 
 @img_comp(["plot_result"], **IMG_KWARGS)
 def test_plot_result():
+    """
+    Compares plots generated by ``plot_result``.
+    """
     num_iterations = 100
 
     data = read("data/RC208.txt", "solomon", round_func="trunc")
@@ -62,12 +68,18 @@ def test_plot_result():
 
 @img_comp(["plot_instance"], **IMG_KWARGS)
 def test_plot_instance():
+    """
+    Compares plots generated by ``plot_instance``.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     plotting.plot_instance(data)
 
 
 @img_comp(["plot_route_schedule"], **IMG_KWARGS)
 def test_plot_route_schedule():
+    """
+    Compares plots generated by ``plot_route_schedule``.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     bks = read_solution("data/RC208.sol")
     sol = Solution(data, bks)
diff --git a/pyvrp/search/tests/test_LocalSearch.py b/pyvrp/search/tests/test_LocalSearch.py
index 405415e..4dc71ad 100644
--- a/pyvrp/search/tests/test_LocalSearch.py
+++ b/pyvrp/search/tests/test_LocalSearch.py
@@ -21,20 +21,13 @@ from pyvrp.search._search import LocalSearch as cpp_LocalSearch
 from pyvrp.tests.helpers import make_heterogeneous, read
 
 
-def test_local_search_returns_same_solution_when_there_are_no_operators():
-    data = read("data/OkSmall.txt")
-    cost_evaluator = CostEvaluator(20, 6)
-    rng = RandomNumberGenerator(seed=42)
-
-    ls = LocalSearch(data, rng, compute_neighbours(data))
-    sol = Solution.make_random(data, rng)
-
-    # No operators have been added, so these calls should be no-ops.
-    assert_equal(ls.search(sol, cost_evaluator), sol)
-    assert_equal(ls.intensify(sol, cost_evaluator), sol)
-
-
 def test_local_search_returns_same_solution_with_empty_neighbourhood():
+    """
+    Tests that calling the local search when it only has node operators and
+    an empty neighbourhood is a no-op: since the node operators respect the
+    neighbourhood definition, they cannot do anything with an empty
+    neighbourhood.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -52,6 +45,10 @@ def test_local_search_returns_same_solution_with_empty_neighbourhood():
 
 @mark.parametrize("size", [1, 2, 3, 4, 6, 7])  # num_clients + 1 == 5
 def test_local_search_raises_when_neighbourhood_dimensions_do_not_match(size):
+    """
+    Tests that the local search raises when the neighbourhood size does not
+    correspond to the problem dimensions.
+    """
     data = read("data/OkSmall.txt")
     rng = RandomNumberGenerator(seed=42)
 
@@ -68,6 +65,10 @@ def test_local_search_raises_when_neighbourhood_dimensions_do_not_match(size):
 
 
 def test_local_search_raises_when_neighbourhood_contains_self_or_depot():
+    """
+    Tests that the local search raises when the granular neighbourhood contains
+    the depot (for any client) or the client is in its own neighbourhood.
+    """
     data = read("data/OkSmall.txt")
     rng = RandomNumberGenerator(seed=42)
 
@@ -100,6 +101,9 @@ def test_local_search_set_get_neighbours(
     symmetric_proximity: bool,
     symmetric_neighbours: bool,
 ):
+    """
+    Tests setting and getting neighbours on the local search instance.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=42)
 
@@ -196,6 +200,10 @@ def test_prize_collecting():
 
 
 def test_cpp_shuffle_results_in_different_solution():
+    """
+    Tests that calling shuffle changes the evaluation order, which can well
+    result in different solutions generated from the same initial solution.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=42)
 
@@ -220,11 +228,16 @@ def test_cpp_shuffle_results_in_different_solution():
 
 
 def test_route_vehicle_types_are_preserved_for_locally_optimal_solutions():
-    # This test tests that we will preserve vehicle types
+    """
+    Tests that a solution that is already locally optimal returns the same
+    solution, particularly w.r.t. the underlying vehicles. This exercises an
+    issue where loading the solution in the local search did not preserve the
+    vehicle types.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=42)
-
     neighbours = compute_neighbours(data)
+
     ls = cpp_LocalSearch(data, neighbours)
     ls.add_node_operator(Exchange10(data))
     ls.add_node_operator(Exchange11(data))
@@ -232,8 +245,6 @@ def test_route_vehicle_types_are_preserved_for_locally_optimal_solutions():
     cost_evaluator = CostEvaluator(1, 1)
     sol = Solution.make_random(data, rng)
 
-    # LocalSearch::search is deterministic, so two calls with the same base
-    # solution should result in the same improved solution.
     improved = ls.search(sol, cost_evaluator)
 
     # Now make the instance heterogeneous and update the local search
@@ -249,7 +260,7 @@ def test_route_vehicle_types_are_preserved_for_locally_optimal_solutions():
     improved = Solution(data, routes)
 
     # Doing the search should not find any further improvements thus not change
-    # the solution, especially not change the vehicle types
+    # the solution.
     further_improved = ls.search(improved, cost_evaluator)
     assert_equal(further_improved, improved)
 
@@ -279,6 +290,10 @@ def test_bugfix_vehicle_type_offsets():
 
 
 def test_intensify_overlap_tolerance():
+    """
+    Tests that the local search's intensifying route operators respect the
+    route overlap tolerance argument.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=42)
 
@@ -303,6 +318,10 @@ def test_intensify_overlap_tolerance():
 
 @mark.parametrize("tol", [-1.0, -0.01, 1.01, 10.9, 1000])
 def test_intensify_overlap_tolerance_raises_outside_unit_interval(tol):
+    """
+    Tests that calling ``intensify()`` raises when the overlap tolerance
+    argument is not in [0, 1].
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=42)
 
@@ -318,6 +337,11 @@ def test_intensify_overlap_tolerance_raises_outside_unit_interval(tol):
 
 
 def test_no_op_results_in_same_solution():
+    """
+    Tests that calling local search without first adding node or route
+    operators is a no-op, and returns the same solution as the one that was
+    given to it.
+    """
     data = read("data/OkSmall.txt")
     rng = RandomNumberGenerator(seed=42)
 
@@ -333,6 +357,10 @@ def test_no_op_results_in_same_solution():
 
 
 def test_intensify_can_improve_solution_further():
+    """
+    Tests that ``intensify()`` improves a solution further once ``search()`` is
+    stuck.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     rng = RandomNumberGenerator(seed=11)
 
@@ -365,6 +393,10 @@ def test_intensify_can_improve_solution_further():
 
 
 def test_local_search_raises_for_incomplete_solutions():
+    """
+    Tests that the local search object cannot (yet) improve solutions that are
+    incomplete. Passing an incomplete solution should raise.
+    """
     data = read("data/OkSmallPrizes.txt")
     rng = RandomNumberGenerator(seed=42)
 
diff --git a/pyvrp/search/tests/test_MoveTwoClientsReversed.py b/pyvrp/search/tests/test_MoveTwoClientsReversed.py
index 528f775..6b4f2d7 100644
--- a/pyvrp/search/tests/test_MoveTwoClientsReversed.py
+++ b/pyvrp/search/tests/test_MoveTwoClientsReversed.py
@@ -50,6 +50,9 @@ def test_single_route_OkSmall():
 
 @mark.parametrize("seed", [2643, 2742, 2941, 3457, 4299, 4497, 6178, 6434])
 def test_RC208_instance(seed: int):
+    """
+    Test a larger instance over several seeds.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=seed)
diff --git a/pyvrp/search/tests/test_Route.py b/pyvrp/search/tests/test_Route.py
index f1714ef..88bd087 100644
--- a/pyvrp/search/tests/test_Route.py
+++ b/pyvrp/search/tests/test_Route.py
@@ -10,6 +10,9 @@ from pyvrp.tests.helpers import make_heterogeneous, read
 
 @pytest.mark.parametrize("loc", [0, 1, 10])
 def test_node_init(loc: int):
+    """
+    Tests that after initialisation, a Node has index 0 and is not in a route.
+    """
     node = Node(loc=loc)
     assert_equal(node.client, loc)
     assert_equal(node.idx, 0)
@@ -18,6 +21,10 @@ def test_node_init(loc: int):
 
 @pytest.mark.parametrize(("idx", "vehicle_type"), [(0, 0), (1, 0), (1, 1)])
 def test_route_init(idx: int, vehicle_type: int):
+    """
+    Tests that after initialisation, a Route has the given index and vehicle
+    type.
+    """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(data, [VehicleType(1, 1), VehicleType(2, 2)])
 
@@ -28,11 +35,21 @@ def test_route_init(idx: int, vehicle_type: int):
 
 @pytest.mark.parametrize("loc", [0, 1, 10])
 def test_new_nodes_are_not_depots(loc: int):
+    """
+    Tests that new nodes cannot be depots: to be a depot in this context
+    requires route membership, and new nodes are not in routes yet.
+    """
     node = Node(loc=loc)
     assert_(not node.is_depot())
 
 
 def test_route_insert_and_remove_updates_node_idx_and_route_properties():
+    """
+    Tests that after a node is inserted into a route, its index and route
+    properties are updated to reflect its new position in the route. After
+    the node is removed from the route, these properties revert to their
+    defaults.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -53,6 +70,10 @@ def test_route_insert_and_remove_updates_node_idx_and_route_properties():
 
 
 def test_route_depots_are_depots():
+    """
+    Tests that the start and end depot nodes in a route known they are, in
+    fact, depots.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -65,6 +86,9 @@ def test_route_depots_are_depots():
 
 
 def test_route_append_increases_route_len():
+    """
+    Tests that appending nodes to a route increases the route's length.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
     assert_equal(len(route), 0)
@@ -81,18 +105,22 @@ def test_route_append_increases_route_len():
 
 
 def test_route_insert():
+    """
+    Tests that inserting and appending nodes works as expected: appending adds
+    to the end, inserting places at the given index.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
     assert_equal(len(route), 0)
 
-    # Insert a few customers so we have an actual route.
+    # Insert a few nodes so we have an actual route.
     route.append(Node(loc=1))
     route.append(Node(loc=2))
     assert_equal(len(route), 2)
     assert_equal(route[1].client, 1)
     assert_equal(route[2].client, 2)
 
-    # # Now insert a new customer at index 1.
+    # # Now insert a new nodes at index 1.
     route.insert(1, Node(loc=3))
     assert_equal(len(route), 3)
     assert_equal(route[1].client, 3)
@@ -101,6 +129,10 @@ def test_route_insert():
 
 
 def test_route_iter_returns_all_clients():
+    """
+    Tests that iterating over a route returns all clients in the route, but
+    not the depots.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -118,6 +150,9 @@ def test_route_iter_returns_all_clients():
 
 
 def test_route_add_and_delete_client_leaves_route_empty():
+    """
+    Tests that adding and removing a client leaves a route empty.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -162,6 +197,9 @@ def test_route_clear_empties_entire_route(num_nodes: int):
 
 
 def test_excess_load():
+    """
+    Tests that the route calculations excess load correctly.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -180,6 +218,10 @@ def test_excess_load():
 
 @pytest.mark.parametrize("client", [1, 2, 3, 4])
 def test_dist_and_load_for_single_client_routes(client: int):
+    """
+    Tests that the route calculates distance and load correctly for a
+    single-client route.
+    """
     data = read("data/OkSmall.txt")
 
     route = Route(data, idx=0, vehicle_type=0)
@@ -204,6 +246,9 @@ def test_dist_and_load_for_single_client_routes(client: int):
 
 
 def test_route_overlaps_with_self_no_matter_the_tolerance_value():
+    """
+    Tests that a route always overlaps with itself.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
     route.append(Node(loc=1))
@@ -217,6 +262,10 @@ def test_route_overlaps_with_self_no_matter_the_tolerance_value():
 
 
 def test_all_routes_overlap_with_maximum_tolerance_value():
+    """
+    Tests that any route overlaps with any other route with the maximum
+    tolerance value.
+    """
     data = read("data/OkSmall.txt")
 
     route1 = Route(data, idx=0, vehicle_type=0)
@@ -239,10 +288,11 @@ def test_all_routes_overlap_with_maximum_tolerance_value():
     assert_(route2.overlaps_with(route1, 1))
 
 
-# TODO test overlap with less extreme cases, including wrap around etc.
-
-
 def test_data_is_not_updated_until_update_call():
+    """
+    Tests that route statistics like distance, time window segments, and load
+    are not updated until ``update()`` has been called.
+    """
     data = read("data/OkSmall.txt")
     route = Route(data, idx=0, vehicle_type=0)
 
@@ -285,6 +335,10 @@ def test_str_contains_route(locs: List[int]):
 
 
 def test_route_tws_access():
+    """
+    Tests access to a client's or depot's time window segment, as tracked by
+    the route.
+    """
     data = read("data/OkSmall.txt")
 
     route = Route(data, idx=0, vehicle_type=0)
@@ -308,6 +362,10 @@ def test_route_tws_access():
 
 @pytest.mark.parametrize("loc", [1, 2, 3, 4])
 def test_tws_between_same_client_returns_node_tws(loc: int):
+    """
+    Tests that calling the ``tws_between()`` with the same start and end
+    arguments returns a node's time window segment data.
+    """
     data = read("data/OkSmall.txt")
     client = data.client(loc)
 
@@ -326,6 +384,10 @@ def test_tws_between_same_client_returns_node_tws(loc: int):
 
 
 def test_tws_between_same_as_tws_before_after_when_one_side_is_depot():
+    """
+    Tests that ``tws_between()`` returns the same value as ``tws_before()`` or
+    ``tws_after()`` when one side is the depot.
+    """
     data = read("data/OkSmall.txt")
 
     route = Route(data, idx=0, vehicle_type=0)
@@ -349,6 +411,10 @@ def test_tws_between_same_as_tws_before_after_when_one_side_is_depot():
 
 
 def test_tws_between_single_route_solution_has_time_warp_in_the_right_places():
+    """
+    Tests time window segment access on a single-route solution where we know
+    exactly where in the route time warp occurs.
+    """
     data = read("data/OkSmall.txt")
 
     route = Route(data, idx=0, vehicle_type=0)
diff --git a/pyvrp/search/tests/test_SwapStar.py b/pyvrp/search/tests/test_SwapStar.py
index e86a008..fccbefc 100644
--- a/pyvrp/search/tests/test_SwapStar.py
+++ b/pyvrp/search/tests/test_SwapStar.py
@@ -60,6 +60,9 @@ def test_swap_star_identifies_additional_moves_over_regular_swap():
 
 @mark.parametrize("seed", [2643, 2742, 2941, 3457, 4299, 4497, 6178, 6434])
 def test_swap_star_on_RC208_instance(seed: int):
+    """
+    Evaluate SWAP* on the RC208 instance, over a few seeds.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=seed)
diff --git a/pyvrp/search/tests/test_TwoOpt.py b/pyvrp/search/tests/test_TwoOpt.py
index 4d629e2..c56188a 100644
--- a/pyvrp/search/tests/test_TwoOpt.py
+++ b/pyvrp/search/tests/test_TwoOpt.py
@@ -21,6 +21,9 @@ from pyvrp.tests.helpers import make_heterogeneous, read
 
 
 def test_OkSmall_instance():
+    """
+    Test 2-OPT on the OkSmall instance where we know exactly what's going on.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -55,21 +58,19 @@ def test_OkSmall_instance():
         [VehicleType(8, 1), VehicleType(8, 1)],
     ],
 )
-def test_OkSmall_heterogeneous_capacity(vehicle_types: List[VehicleType]):
-    # This instance tests a two-opt move that is improving when disregarding
-    # the vehicle capacities. Depending on the (heterogeneous) capacities,
-    # the move may or may not be improving and should be applied or not
-    # The starting solution has routes [1, 3] and [2, 4] with demands 8, 10
-    # the 2-opt solution has routes [1, 4] and [2, 3] with demands 10, 8
-
+def test_OkSmall_multiple_vehicle_types(vehicle_types: List[VehicleType]):
+    """
+    This test evaluates a 2-OPT move that is improving for some of the vehicle
+    types, and not for others. Because the granular neighbourhood is very
+    constrained, only a single 2-OPT can be applied. When it's better to do so,
+    it should have been applied, and when it's not better, it should not.
+    """
     data = read("data/OkSmall.txt")
-    cost_evaluator = CostEvaluator(10000, 6)  # Large capacity penalty
-    rng = RandomNumberGenerator(seed=42)
-
-    # Now making it heterogenous, the same move should result in capacity
-    # penalties and thus not be applied
     data = make_heterogeneous(data, vehicle_types)
 
+    cost_evaluator = CostEvaluator(10_000, 6)  # large capacity penalty
+    rng = RandomNumberGenerator(seed=42)
+
     neighbours: List[List[int]] = [[], [2], [], [], []]  # only 1 -> 2
     ls = LocalSearch(data, rng, neighbours)
     ls.add_node_operator(TwoOpt(data))
@@ -81,7 +82,7 @@ def test_OkSmall_heterogeneous_capacity(vehicle_types: List[VehicleType]):
 
     assert_(not np.allclose(cost1, cost2))
 
-    # Using the local search, the result should not get worse
+    # Using the local search, the result should not get worse.
     improved_sol1 = ls.search(sol1, cost_evaluator)
     improved_sol2 = ls.search(sol2, cost_evaluator)
 
@@ -92,6 +93,9 @@ def test_OkSmall_heterogeneous_capacity(vehicle_types: List[VehicleType]):
 
 @mark.parametrize("seed", [2643, 2742, 2941, 3457, 4299, 4497, 6178, 6434])
 def test_RC208_instance(seed: int):
+    """
+    Test a larger instance over several seeds.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=seed)
diff --git a/pyvrp/search/tests/test_neighbourhood.py b/pyvrp/search/tests/test_neighbourhood.py
index 0dd4e07..bcac450 100644
--- a/pyvrp/search/tests/test_neighbourhood.py
+++ b/pyvrp/search/tests/test_neighbourhood.py
@@ -28,6 +28,9 @@ def test_neighbourhood_params_raises_for_invalid_arguments(
     symmetric_proximity: bool,
     symmetric_neighbours: bool,
 ):
+    """
+    Test that ``NeighbourhoodParams`` raises for invalid configurations.
+    """
     with assert_raises(ValueError):
         NeighbourhoodParams(
             weight_wait_time,
@@ -60,6 +63,9 @@ def test_neighbourhood_params_does_not_raise_for_valid_arguments(
     symmetric_proximity: bool,
     symmetric_neighbours: bool,
 ):
+    """
+    Tests that ``NeighbourhoodParams`` allows valid arguments and edge cases.
+    """
     NeighbourhoodParams(
         weight_wait_time,
         weight_time_warp,
@@ -85,7 +91,7 @@ def test_neighbourhood_params_does_not_raise_for_valid_arguments(
          {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}),
-        # From original c++ implementation
+        # 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}),
@@ -104,6 +110,9 @@ def test_compute_neighbours(
     idx_check: int,
     expected_neighbours_check: Set[int],
 ):
+    """
+    Tests ``compute_neighbours`` on several well-understood cases.
+    """
     data = read("data/RC208.txt", "solomon", "trunc1")
     params = NeighbourhoodParams(
         weight_wait_time,
@@ -129,6 +138,9 @@ def test_compute_neighbours(
 
 
 def test_neighbours_are_sorted_by_proximity():
+    """
+    Tests that the neighbourhood lists sort by their proximity: closest first.
+    """
     # 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)
@@ -146,6 +158,10 @@ def test_neighbours_are_sorted_by_proximity():
 
 
 def test_symmetric_neighbours():
+    """
+    Tests that when ``symmetric_neighbours`` is true, if an edge (i, j) is
+    in the granular neighbourhood, then so is (j, i).
+    """
     data = read("data/RC208.txt", "solomon", "trunc1")
 
     # Symmetric neighbourhood structure: if (i, j) is in, then so is (j, i).
@@ -163,6 +179,10 @@ def test_symmetric_neighbours():
 
 
 def test_more_neighbours_than_instance_size():
+    """
+    Tests that a value for ``nb_granular`` larger than the problem instance's
+    number of clients results in neighbourhoods of maximum size.
+    """
     data = read("data/RC208.txt", "solomon", round_func="trunc")
     params = NeighbourhoodParams(nb_granular=data.num_clients)
     neighbours = compute_neighbours(data, params)
@@ -172,6 +192,10 @@ def test_more_neighbours_than_instance_size():
 
 
 def test_proximity_with_prizes():
+    """
+    Tests that prizes factor into the neighbourhood structure, and offset
+    travel costs somewhat.
+    """
     data = read("data/p06-2-50.vrp", round_func="dimacs")
     params = NeighbourhoodParams(0, 0, nb_granular=10)
     neighbours = compute_neighbours(data, params)
diff --git a/pyvrp/stop/tests/test_MaxIterations.py b/pyvrp/stop/tests/test_MaxIterations.py
index ee33d83..eaed0ba 100644
--- a/pyvrp/stop/tests/test_MaxIterations.py
+++ b/pyvrp/stop/tests/test_MaxIterations.py
@@ -22,13 +22,21 @@ def test_does_not_raise(max_iterations: int):
 
 
 def test_before_max_iterations():
+    """
+    Tests that ``MaxIterations`` does not stop *before* the number of
+    iterations have passed.
+    """
     stop = MaxIterations(100)
 
-    for _ in range(100):
+    for _ in range(100):  # stop *after* 100 iterations, not before
         assert_(not stop(1))
 
 
 def test_after_max_iterations():
+    """
+    Tests that after the number of iterations have passed, ``MaxIterations``
+    stops.
+    """
     stop = MaxIterations(100)
 
     for _ in range(100):
diff --git a/pyvrp/stop/tests/test_MaxRuntime.py b/pyvrp/stop/tests/test_MaxRuntime.py
index c1aec4d..6f88045 100644
--- a/pyvrp/stop/tests/test_MaxRuntime.py
+++ b/pyvrp/stop/tests/test_MaxRuntime.py
@@ -24,6 +24,10 @@ def test_valid_parameters(max_runtime: float):
 
 @mark.parametrize("max_runtime", [0.01, 0.05, 0.10])
 def test_before_max_runtime(max_runtime):
+    """
+    Tests that ``MaxRuntime`` does not stop *before* the given runtime has
+    passed.
+    """
     stop = MaxRuntime(max_runtime)
 
     for _ in range(100):
@@ -32,6 +36,9 @@ def test_before_max_runtime(max_runtime):
 
 @mark.parametrize("max_runtime", [0.01, 0.05, 0.10])
 def test_after_max_runtime(max_runtime):
+    """
+    Tests that after the max runtime has passed, ``MaxRuntime`` stops.
+    """
     stop = MaxRuntime(max_runtime)
     assert_(not stop(1))  # trigger the first time measurement
 
diff --git a/pyvrp/stop/tests/test_TimedNoImprovement.py b/pyvrp/stop/tests/test_TimedNoImprovement.py
index 6d79b2e..f057c2c 100644
--- a/pyvrp/stop/tests/test_TimedNoImprovement.py
+++ b/pyvrp/stop/tests/test_TimedNoImprovement.py
@@ -25,12 +25,20 @@ def test_valid_parameters(max_iterations, max_runtime):
 
 
 def test_stops_if_zero_max_iterations():
+    """
+    Tests that ``TimedNoImprovement`` stops immediately if zero maximum
+    iterations are allowed.
+    """
     stop = TimedNoImprovement(0, 0.100)
     assert_(stop(1))
 
 
 @mark.parametrize("max_runtime", [0.01, 0.05, 0.10])
 def test_before_max_runtime(max_runtime):
+    """
+    Tests that ``TimedNoImprovement`` allows iterations if the maximum runtime
+    has not yet been violated.
+    """
     stop = TimedNoImprovement(101, max_runtime)
 
     for _ in range(100):
@@ -39,6 +47,10 @@ def test_before_max_runtime(max_runtime):
 
 @mark.parametrize("max_runtime", [0.01, 0.05, 0.10])
 def test_after_max_runtime(max_runtime):
+    """
+    Tests that ``TimedNoImprovement`` stops after maximum runtime, regardless
+    of maximum number of iterations.
+    """
     stop = TimedNoImprovement(101, max_runtime)
     assert_(not stop(1))  # trigger the first time measurement
 
diff --git a/pyvrp/tests/test_CostEvaluator.py b/pyvrp/tests/test_CostEvaluator.py
index 60a4422..65ffde9 100644
--- a/pyvrp/tests/test_CostEvaluator.py
+++ b/pyvrp/tests/test_CostEvaluator.py
@@ -7,6 +7,9 @@ from pyvrp.tests.helpers import read
 
 
 def test_load_penalty():
+    """
+    This test asserts that load penalty computations are correct.
+    """
     cost_evaluator = CostEvaluator(2, 1)
 
     assert_allclose(cost_evaluator.load_penalty(0, 1), 0)  # below capacity
@@ -29,6 +32,10 @@ def test_load_penalty():
 
 @mark.parametrize("capacity", [5, 15, 29, 51, 103])
 def test_load_penalty_always_zero_when_below_capacity(capacity: int):
+    """
+    This test asserts that load penalties are only applied to excess load, that
+    is, load in excess of the vehicle's capacity.
+    """
     load_penalty = 2
     cost_evaluator = CostEvaluator(load_penalty, 1)
 
@@ -47,6 +54,9 @@ def test_load_penalty_always_zero_when_below_capacity(capacity: int):
 
 
 def test_tw_penalty():
+    """
+    This test asserts that time window penalty computations are correct.
+    """
     cost_evaluator = CostEvaluator(1, 2)
 
     # Penalty per unit time warp is 2
@@ -63,6 +73,11 @@ def test_tw_penalty():
 
 
 def test_cost():
+    """
+    This test asserts that the cost is computed correctly for feasible
+    individuals, and is a large value (representing infinity) for infeasible
+    individuals.
+    """
     data = read("data/OkSmall.txt")
     default_cost_evaluator = CostEvaluator()
     cost_evaluator = CostEvaluator(20, 6)
@@ -88,6 +103,10 @@ def test_cost():
 
 
 def test_cost_with_prizes():
+    """
+    When solving a prize-collecting instance, the cost is equal to the distance
+    plus a prize term.
+    """
     data = read("data/p06-2-50.vrp", round_func="dimacs")
     cost_evaluator = CostEvaluator(1, 1)
 
@@ -106,6 +125,12 @@ def test_cost_with_prizes():
 
 
 def test_penalised_cost():
+    """
+    The penalised cost represents the smoothed objective, where constraint
+    violations are priced in using penalty terms. It can be computed for both
+    feasible and infeasible individuals. In case of the former, it is equal
+    to the actual cost: the penalty terms are all zero.
+    """
     data = read("data/OkSmall.txt")
     penalty_capacity = 20
     penalty_tw = 6
diff --git a/pyvrp/tests/test_DynamicBitset.py b/pyvrp/tests/test_DynamicBitset.py
index 9e2a382..28884cb 100644
--- a/pyvrp/tests/test_DynamicBitset.py
+++ b/pyvrp/tests/test_DynamicBitset.py
@@ -17,12 +17,18 @@ def test_size(num_bits, expected_size):
 
 
 def test_init_all_zero():
+    """
+    The bitset starts out empty, that is, is initialised to all zeros.
+    """
     bitset = DynamicBitset(64)
     assert_(all(not bitset[idx] for idx in range(len(bitset))))
     assert_equal(bitset.count(), 0)
 
 
 def test_eq():
+    """
+    Tests the equality operator.
+    """
     bitset1 = DynamicBitset(64)
     bitset2 = DynamicBitset(64)
     assert_(bitset1 == bitset2)  # both empty, equal size - should be the same
@@ -37,8 +43,11 @@ def test_eq():
 
 
 def test_get_set_item():
+    """
+    Tests that setting and retrieving an item from the bitset works correctly.
+    """
     bitset = DynamicBitset(128)
-    indices = [0, 1, 63, 64, 126, 127]
+    indices = [0, 1, 63, 64, 126, 127]  # at the bitset boundaries
 
     assert_equal(bitset.count(), 0)  # Bitset starts empty - no bits are set
 
@@ -58,6 +67,9 @@ def test_get_set_item():
 
 
 def test_bit_or():
+    """
+    Tests the union operator.
+    """
     bitset1 = DynamicBitset(128)
     bitset1[0] = True
     bitset1[64] = True
@@ -74,6 +86,9 @@ def test_bit_or():
 
 
 def test_bit_and():
+    """
+    Tests the intersection operator.
+    """
     bitset1 = DynamicBitset(128)
     bitset1[0] = True
     bitset1[64] = True
@@ -90,6 +105,9 @@ def test_bit_and():
 
 
 def test_bit_xor():
+    """
+    Tests the symmetric difference operator.
+    """
     bitset1 = DynamicBitset(128)
     bitset1[0] = True
     bitset1[64] = True
@@ -106,6 +124,9 @@ def test_bit_xor():
 
 
 def test_bit_not():
+    """
+    Tests the complement operator.
+    """
     bitset = DynamicBitset(128)
     assert_equal(bitset.count(), 0)
 
diff --git a/pyvrp/tests/test_GeneticAlgorithm.py b/pyvrp/tests/test_GeneticAlgorithm.py
index bf0d8e9..9dfbb98 100644
--- a/pyvrp/tests/test_GeneticAlgorithm.py
+++ b/pyvrp/tests/test_GeneticAlgorithm.py
@@ -1,3 +1,5 @@
+from types import MethodType
+
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 from pytest import mark
 
@@ -5,6 +7,7 @@ from pyvrp import (
     GeneticAlgorithm,
     GeneticAlgorithmParams,
     PenaltyManager,
+    PenaltyParams,
     Population,
     PopulationParams,
     RandomNumberGenerator,
@@ -137,6 +140,11 @@ def test_initial_solutions_added_when_restarting():
 
 
 def test_best_solution_improves_with_more_iterations():
+    """
+    Tests that additional iterations result in better solutions. This is a
+    smoke test that checks at least something's improving during the GA's
+    execution.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     rng = RandomNumberGenerator(seed=42)
     pm = PenaltyManager()
@@ -185,4 +193,39 @@ def test_best_initial_solution():
     assert_equal(result.best, bks)
 
 
-# TODO more functional tests
+def test_never_repairs_when_zero_repair_probability():
+    """
+    Tests that the genetic algorithm does not repair when the repair
+    probability parameter is set to zero.
+    """
+    data = read("data/RC208.txt", "solomon", "dimacs")
+    rng = RandomNumberGenerator(seed=42)
+    pm = PenaltyManager(PenaltyParams(repair_booster=10))
+    pop = Population(bpd)
+
+    ls = LocalSearch(data, rng, compute_neighbours(data))
+    ls.add_node_operator(Exchange10(data))
+
+    init = [Solution.make_random(data, rng) for _ in range(25)]
+
+    # Repair probability 100%, but we're still using the unpatched penalty
+    # manager. This should be OK.
+    ga_params = GeneticAlgorithmParams(repair_probability=1.0)
+    algo = GeneticAlgorithm(data, pm, rng, pop, ls, srex, init, ga_params)
+    algo.run(MaxIterations(50))
+
+    # Now we patch the penalty manager: when asked for a booster cost evaluator
+    # (as used during repair), this will now raise a runtime error. Since the
+    # repair probability is still 100%, this should certainly raise.
+    def raise_when_called(self):
+        raise RuntimeError
+
+    pm.get_booster_cost_evaluator = MethodType(raise_when_called, pm)
+    with assert_raises(RuntimeError):
+        algo.run(MaxIterations(50))
+
+    # But when we set the repair probability to 0%, the booster is no longer
+    # needed, and nothing should raise.
+    ga_params = GeneticAlgorithmParams(repair_probability=0.0)
+    algo = GeneticAlgorithm(data, pm, rng, pop, ls, srex, init, ga_params)
+    algo.run(MaxIterations(50))
diff --git a/pyvrp/tests/test_Model.py b/pyvrp/tests/test_Model.py
index e1c1128..4ef6f38 100644
--- a/pyvrp/tests/test_Model.py
+++ b/pyvrp/tests/test_Model.py
@@ -9,6 +9,10 @@ from pyvrp.tests.helpers import read
 
 
 def test_model_data():
+    """
+    Tests that calling the ``data()`` member on the model returns a data
+    instance representing the added data.
+    """
     model = Model()
 
     # Let's add some data: a single client, and edges from/to the depot.
@@ -18,12 +22,17 @@ def test_model_data():
     model.add_edge(client, depot, 1, 1)
     model.add_vehicle_type(capacity=1, num_available=1)
 
-    # The model should now have one client.
     data = model.data()
     assert_equal(data.num_clients, 1)
+    assert_equal(data.num_vehicle_types, 1)
+    assert_equal(data.num_vehicles, 1)
 
 
 def test_add_depot_raises_more_than_one_depot():
+    """
+    PyVRP does not currently support VRPs with multiple depots. Adding more
+    than one depot should raise.
+    """
     model = Model()
     model.add_depot(0, 0)  # first depot should be OK
 
@@ -32,6 +41,10 @@ def test_add_depot_raises_more_than_one_depot():
 
 
 def test_add_edge_raises_negative_distance_or_duration():
+    """
+    Negative distances or durations are not understood. Attempting to add
+    such edges should raise an error.
+    """
     model = Model()
     depot = model.add_depot(0, 0)
     client = model.add_client(0, 1)
@@ -54,12 +67,21 @@ def test_add_edge_raises_negative_distance_or_duration():
     ],
 )
 def test_add_vehicle_type_raises_negative_number_or_capacity(number, capacity):
+    """
+    A negative number of vehicles of a given type, or a negative capacity, is
+    not understood. Attempting to add such a vehicle type should raise an
+    error.
+    """
     model = Model()
     with assert_raises(ValueError):
         model.add_vehicle_type(capacity=capacity, num_available=number)
 
 
 def test_add_client_attributes():
+    """
+    Smoke test that checks, for a single client, that the model adds a client
+    whose attributes are the same as what was passed in.
+    """
     model = Model()
     client = model.add_client(
         x=1,
@@ -85,6 +107,10 @@ def test_add_client_attributes():
 
 
 def test_add_depot_attributes():
+    """
+    Smoke test that checks the depot attributes are the same as what was passed
+    in.
+    """
     model = Model()
     depot = model.add_depot(x=1, y=0, tw_early=3, tw_late=5)
 
@@ -95,6 +121,9 @@ def test_add_depot_attributes():
 
 
 def test_add_edge():
+    """
+    Smoke test that checks edge attributes are the same as what was passed in.
+    """
     model = Model()
     depot = model.add_depot(0, 0)
     client = model.add_client(0, 1)
@@ -107,6 +136,10 @@ def test_add_edge():
 
 
 def test_add_vehicle_type():
+    """
+    Smoke test that checks vehicle type attributes are the same as what was
+    passed in.
+    """
     model = Model()
     vehicle_type = model.add_vehicle_type(num_available=10, capacity=998)
 
@@ -115,6 +148,9 @@ def test_add_vehicle_type():
 
 
 def test_get_locations():
+    """
+    Checks that the ``locations`` property returns the depot and all clients.
+    """
     model = Model()
     client1 = model.add_client(0, 1)
     depot = model.add_depot(0, 0)
@@ -127,6 +163,9 @@ def test_get_locations():
 
 
 def test_get_vehicle_types():
+    """
+    Tests the ``vehicle_types`` property.
+    """
     model = Model()
     vehicle_type1 = model.add_vehicle_type(1, 2)
     vehicle_type2 = model.add_vehicle_type(1, 3)
@@ -137,6 +176,10 @@ def test_get_vehicle_types():
 
 
 def test_from_data():
+    """
+    Tests that initialising the model from a data instance results in a valid
+    model representation of that data instance.
+    """
     read_data = read("data/E-n22-k4.txt", round_func="dimacs")
     model = Model.from_data(read_data)
     model_data = model.data()
@@ -155,6 +198,10 @@ def test_from_data():
 
 
 def test_from_data_and_solve():
+    """
+    Tests that solving a model initialised from a data instance finds the
+    correct (known) solutions.
+    """
     # 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)
@@ -170,6 +217,10 @@ def test_from_data_and_solve():
 
 
 def test_model_and_solve():
+    """
+    Tests that solving a model initialised using the modelling interface
+    finds the correct (known) solutions.
+    """
     # Solve the small OkSmall instance using the from_data constructor.
     data = read("data/OkSmall.txt")
     model = Model.from_data(data)
@@ -217,6 +268,11 @@ def test_model_and_solve():
 
 
 def test_partial_distance_duration_matrix():
+    """
+    Tests that adding a full distance or duration matrix is not required. Any
+    "missing" arcs are given large default values, ensuring they are never
+    used.
+    """
     model = Model()
     depot = model.add_depot(0, 0)
     clients = [model.add_client(0, 1), model.add_client(1, 0)]
@@ -242,6 +298,12 @@ def test_partial_distance_duration_matrix():
 
 
 def test_data_warns_about_scaling_issues(recwarn):
+    """
+    Tests that the modelling interface warns when an edge is added whose
+    distance or duration values are too large: in that case, PyVRP might not be
+    able to detect missing edges. This situation is likely caused by scaling
+    issues, so a warning is appropriate.
+    """
     model = Model()
     model.add_vehicle_type(capacity=0, num_available=1)
     depot = model.add_depot(0, 0)
@@ -250,13 +312,15 @@ def test_data_warns_about_scaling_issues(recwarn):
     # 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)
+    # But a value exceeding the maximum user value is not OK. This should warn
+    # (both for distance and/or duration).
+    with assert_warns(ScalingWarning):
+        model.add_edge(depot, client, distance=MAX_USER_VALUE + 1)
+
     with assert_warns(ScalingWarning):
-        model.data()
+        model.add_edge(depot, client, distance=0, duration=MAX_USER_VALUE + 1)
 
 
 def test_model_solves_instance_with_zero_or_one_clients():
diff --git a/pyvrp/tests/test_PenaltyManager.py b/pyvrp/tests/test_PenaltyManager.py
index 9995fc5..6288f66 100644
--- a/pyvrp/tests/test_PenaltyManager.py
+++ b/pyvrp/tests/test_PenaltyManager.py
@@ -33,6 +33,9 @@ def test_constructor_throws_when_arguments_invalid(
     penalty_decrease: float,
     target_feasible: float,
 ):
+    """
+    Tests that invalid arguments are not accepted.
+    """
     with assert_raises(ValueError):
         PenaltyParams(
             init_load_penalty,
@@ -46,6 +49,10 @@ def test_constructor_throws_when_arguments_invalid(
 
 
 def test_repair_booster():
+    """
+    Tests that the booster evaluator returns a cost evaluator object that
+    penalises constraint violations much more severely.
+    """
     params = PenaltyParams(1, 1, 5, 1, 1, 1, 1)
     pm = PenaltyManager(params)
 
@@ -69,6 +76,10 @@ def test_repair_booster():
 
 
 def test_capacity_penalty_update_increase():
+    """
+    Tests that the load violation penalty is increased when too many load
+    infeasible solutions have been generated since the last update.
+    """
     num_registerations = 4
     params = PenaltyParams(1, 1, 1, num_registerations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params)
@@ -108,6 +119,10 @@ def test_capacity_penalty_update_increase():
 
 
 def test_capacity_penalty_update_decrease():
+    """
+    Tests that the load violation penalty is decreased when sufficiently many
+    load feasible solutions have been generated since the last update.
+    """
     num_registerations = 4
     params = PenaltyParams(4, 1, 1, num_registerations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params)
@@ -146,6 +161,10 @@ def test_capacity_penalty_update_decrease():
 
 
 def test_time_warp_penalty_update_increase():
+    """
+    Tests that the time warp violation penalty is increased when too many time
+    window infeasible solutions have been generated since the last update.
+    """
     num_registerations = 4
     params = PenaltyParams(1, 1, 1, num_registerations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params)
@@ -185,6 +204,11 @@ def test_time_warp_penalty_update_increase():
 
 
 def test_time_warp_penalty_update_decrease():
+    """
+    Tests that the time warp violation penalty is decreased when sufficently
+    many time window feasible solutions have been generated since the last
+    update.
+    """
     num_registerations = 4
     params = PenaltyParams(1, 4, 1, num_registerations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params)
@@ -223,6 +247,10 @@ def test_time_warp_penalty_update_decrease():
 
 
 def test_does_not_update_penalties_before_sufficient_registrations():
+    """
+    Tests that updates only happen every ``num_registrations`` times, not every
+    time a new value is registered.
+    """
     num_registerations = 4
     params = PenaltyParams(4, 4, 1, num_registerations, 1.1, 0.9, 0.5)
     pm = PenaltyManager(params)
diff --git a/pyvrp/tests/test_Population.py b/pyvrp/tests/test_Population.py
index baa8065..195920f 100644
--- a/pyvrp/tests/test_Population.py
+++ b/pyvrp/tests/test_Population.py
@@ -108,6 +108,11 @@ def test_params_constructor_does_not_raise_when_arguments_valid(
 
 
 def test_add_triggers_purge():
+    """
+    Tests that adding another solution to a population of maximum size triggers
+    survivor selection, that is, a purge that reduces the relevant population
+    back to minimum size.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -151,6 +156,10 @@ def test_add_triggers_purge():
 
 
 def test_select_returns_same_parents_if_no_other_option():
+    """
+    Tests that the ``select()`` method tries to return two different parents,
+    but will return the same solution twice if there is no other option.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=2_147_483_647)
@@ -187,10 +196,12 @@ def test_select_returns_same_parents_if_no_other_option():
     assert_(900 < different_parents < 1_000)
 
 
-# // TODO test more select() - diversity, feas/infeas pairs
-
-
 def test_population_is_empty_with_zero_min_pop_size_and_generation_size():
+    """
+    Tests that the population can never grow when it starts empty and the
+    generation size is set to 0: adding a new solution then immediately
+    triggers survivor selection.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=12)
@@ -210,6 +221,11 @@ def test_population_is_empty_with_zero_min_pop_size_and_generation_size():
 
 @mark.parametrize("nb_elite", [5, 25])
 def test_elite_solutions_are_not_purged(nb_elite: int):
+    """
+    Tests that elite solutions - those considered of such high quality that
+    they should be given special treatment - are not purged during survivor
+    selection.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     params = PopulationParams(nb_elite=nb_elite)
@@ -246,6 +262,10 @@ def test_elite_solutions_are_not_purged(nb_elite: int):
 
 @mark.parametrize("k", [2, 3])
 def test_tournament_ranks_by_fitness(k: int):
+    """
+    Tests that the tournament-based parent selection on average returns
+    solutions about as often as their relative fitness value would suggest.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -292,6 +312,9 @@ def test_tournament_ranks_by_fitness(k: int):
 
 @mark.parametrize("k", [-100, -1, 0])  # k must be strictly positive
 def test_tournament_raises_for_invalid_k(k: int):
+    """
+    Tests that k >= 0 is required for tournament-based selection.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -305,6 +328,10 @@ def test_tournament_raises_for_invalid_k(k: int):
 
 
 def test_purge_removes_duplicates():
+    """
+    Tests that purging/survivor selection first removes duplicate solutions,
+    before purging by (biased) fitness.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     params = PopulationParams(min_pop_size=20, generation_size=5)
@@ -344,6 +371,9 @@ def test_purge_removes_duplicates():
 
 
 def test_clear():
+    """
+    Tests that clearing the population reduces its size to zero.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -359,6 +389,11 @@ def test_clear():
 
 
 def test_add_emits_warning_when_solution_is_empty():
+    """
+    Tests that adding an empty solution to the population emits a warning. Such
+    solutions could be generated when solving a prize-collecting instance, and
+    typically indicate a scaling issue in the data.
+    """
     data = read("data/p06-2-50.vrp", round_func="dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     pop = Population(bpd)
diff --git a/pyvrp/tests/test_ProblemData.py b/pyvrp/tests/test_ProblemData.py
index 559404c..0510d13 100644
--- a/pyvrp/tests/test_ProblemData.py
+++ b/pyvrp/tests/test_ProblemData.py
@@ -38,6 +38,10 @@ def test_client_constructor_initialises_data_fields_correctly(
     release_time: int,
     prize: int,
 ):
+    """
+    Tests that the access properties return the data that was given to the
+    Client's constructor.
+    """
     client = Client(
         x, y, demand, service_duration, tw_early, tw_late, release_time, prize
     )
@@ -79,6 +83,9 @@ def test_raises_for_invalid_client_data(
     release_time: int,
     prize: int,
 ):
+    """
+    Tests that invalid configurations are not accepted.
+    """
     with assert_raises(ValueError):
         Client(x, y, demand, service, tw_early, tw_late, release_time, prize)
 
@@ -101,6 +108,9 @@ def test_raises_for_invalid_depot_data(
     release_time: int,
     prize: int,
 ):
+    """
+    Tests that an invalid depot configuration is not accepted.
+    """
     depot = Client(
         x, y, demand, service, tw_early, tw_late, release_time, prize
     )
@@ -163,6 +173,9 @@ def test_problem_data_raises_when_incorrect_matrix_dimensions(matrix):
 
 
 def test_centroid():
+    """
+    Tests the computation of the centroid of all clients in the data instance.
+    """
     data = read("data/OkSmall.txt")
 
     centroid = data.centroid()
@@ -205,6 +218,14 @@ def test_matrix_access():
 
 
 def test_matrices_are_not_writeable():
+    """
+    Tests that the data matrices provided by ``distance_matrix()`` and
+    ``duration_matrix()`` are not writeable. They can be read from, but
+    assigning new values should raise an error.
+
+    We require this because they're constant on the C++ side, and allowing
+    changes from Python causes undefined behaviour on the C++ side.
+    """
     data = ProblemData(
         clients=[Client(x=0, y=0)],
         vehicle_types=[VehicleType(1, 2)],
diff --git a/pyvrp/tests/test_RandomNumberGenerator.py b/pyvrp/tests/test_RandomNumberGenerator.py
index 659beed..aa244af 100644
--- a/pyvrp/tests/test_RandomNumberGenerator.py
+++ b/pyvrp/tests/test_RandomNumberGenerator.py
@@ -8,11 +8,19 @@ from pyvrp import RandomNumberGenerator
 
 
 def test_bounds():
+    """
+    Tests that the minimum and maximum integer value the RNG can produce
+    corresponds to the underlying type's maximum range.
+    """
     assert_equal(RandomNumberGenerator.min(), 0)
     assert_equal(RandomNumberGenerator.max(), np.iinfo(np.uint32).max)
 
 
 def test_call():
+    """
+    Calling the RNG as a function returns random integers, in a fixed sequence
+    depending on the seed.
+    """
     rng = RandomNumberGenerator(seed=42)
 
     assert_equal(rng(), 2386648076)
@@ -25,6 +33,10 @@ def test_call():
 
 
 def test_randint():
+    """
+    The ``randint(high)`` function returns a random integer between
+    ``[0, high]``. Internally, this relies on ``__call__`` (see previous test).
+    """
     rng = RandomNumberGenerator(seed=42)
 
     assert_equal(rng.randint(100), 2386648076 % 100)
@@ -33,6 +45,10 @@ def test_randint():
 
 @mark.parametrize("seed", [2, 10, 42])
 def test_rand(seed: int):
+    """
+    Tests that repeatedly calling ``rand()`` should result in a sample that is
+    approximately uniformly distributed.
+    """
     rng = RandomNumberGenerator(seed)
     sample = np.array([rng.rand() for _ in range(10_000)])
 
@@ -45,5 +61,9 @@ def test_rand(seed: int):
 
 @mark.parametrize("state", [[1, 2, 3, 4], [10, 14, 274, 83]])
 def test_rng_has_given_state(state: List[int]):
+    """
+    Tests that setting the RNG with a given state, and then requesting that
+    state, returns the same state.
+    """
     rng = RandomNumberGenerator(state=state)
     assert_equal(rng.state(), state)
diff --git a/pyvrp/tests/test_Result.py b/pyvrp/tests/test_Result.py
index 950cd85..65a9e74 100644
--- a/pyvrp/tests/test_Result.py
+++ b/pyvrp/tests/test_Result.py
@@ -16,6 +16,10 @@ from pyvrp.tests.helpers import read
     [([[1, 2], [3], [4]], 1, 1.5), ([[1, 2, 3, 4]], 100, 54.2)],
 )
 def test_fields_are_correctly_set(routes, num_iterations, runtime):
+    """
+    Tests that ``Result``'s data properties are correctly set after
+    initialisation completes.
+    """
     data = read("data/OkSmall.txt")
     sol = Solution(data, routes)
 
@@ -39,6 +43,9 @@ def test_fields_are_correctly_set(routes, num_iterations, runtime):
     ],
 )
 def test_init_raises_invalid_arguments(num_iterations, runtime):
+    """
+    Tests that invalid arguments are rejected.
+    """
     data = read("data/OkSmall.txt")
     sol = Solution(data, [[1, 2, 3, 4]])
 
@@ -48,6 +55,9 @@ def test_init_raises_invalid_arguments(num_iterations, runtime):
 
 @mark.parametrize("num_iterations", [0, 1, 10])
 def test_num_iterations(num_iterations: int):
+    """
+    Tests access to the ``num_iterations`` property.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -67,6 +77,10 @@ def test_num_iterations(num_iterations: int):
     [[[1, 2], [3], [4]], [[1, 2, 3, 4]]],
 )
 def test_str_contains_essential_information(routes):
+    """
+    Tests that printing (or, in general, calling ```str(result)``) returns a
+    bunch of useful information about the underlying solution.
+    """
     data = read("data/OkSmall.txt")
 
     sol = Solution(data, routes)
@@ -87,6 +101,10 @@ def test_str_contains_essential_information(routes):
 
 @mark.parametrize("num_iterations", [0, 1, 10])
 def test_result_can_be_pickled(num_iterations: int):
+    """
+    Tests that a ``Result`` object can be pickled: it can be serialised and
+    unserialised. This is useful for e.g. storing results to disk.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
diff --git a/pyvrp/tests/test_Solution.py b/pyvrp/tests/test_Solution.py
index 9976169..cc38d80 100644
--- a/pyvrp/tests/test_Solution.py
+++ b/pyvrp/tests/test_Solution.py
@@ -17,6 +17,9 @@ from pyvrp.tests.helpers import make_heterogeneous, read
 
 
 def test_route_constructor_raises_for_empty_routes():
+    """
+    Tests that constructing a ``Solution`` with empty routes fails.
+    """
     data = read("data/OkSmall.txt")
 
     with assert_raises(RuntimeError):
@@ -25,28 +28,32 @@ def test_route_constructor_raises_for_empty_routes():
         Solution(data, [[3, 4], [], [1, 2]])
 
 
-def test_route_constructor_heterogeneous():
+def test_route_constructor_with_different_vehicle_types():
+    """
+    Tests that Solution's route constructor respects vehicle types.
+    """
     data = read("data/OkSmall.txt")
-    # Test heterogeneous case
     data = make_heterogeneous(data, [VehicleType(10, 1), VehicleType(20, 2)])
-    sol = Solution(data, [Route(data, [3, 4], 0), Route(data, [1, 2], 1)])
 
-    # num_routes() should show two non-empty routes.
-    assert_equal(sol.num_routes(), 2)
+    sol = Solution(data, [Route(data, [3, 4], 0), Route(data, [1, 2], 1)])
 
-    # We expect Solution to remove empty routes and return routes with the
-    # correct vehicle types.
+    # We expect Solution to return routes with the correct vehicle types.
     routes = sol.get_routes()
     assert_equal(len(routes), 2)
+
     assert_equal(routes[0].visits(), [3, 4])
     assert_equal(routes[0].vehicle_type(), 0)
     assert_equal(routes[0], Route(data, [3, 4], 0))
+
     assert_equal(routes[1].visits(), [1, 2])
     assert_equal(routes[1].vehicle_type(), 1)
     assert_equal(routes[1], Route(data, [1, 2], 1))
 
 
 def test_route_eq():
+    """
+    Tests ``Route``'s equality operator.
+    """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(data, [VehicleType(10, 1), VehicleType(20, 2)])
 
@@ -68,6 +75,10 @@ def test_route_eq():
 
 
 def test_random_constructor_cycles_over_routes():
+    """
+    Tests that a randomly constructed solution fills all available vehicles
+    in turn.
+    """
     # This instance has four clients and three vehicles. Since 1 client per
     # vehicle would not work (insufficient vehicles), each route is given two
     # clients (and the last route should be empty).
@@ -85,6 +96,10 @@ def test_random_constructor_cycles_over_routes():
 
 
 def test_route_constructor_raises_too_many_vehicles():
+    """
+    Tests that constructing a solution with more routes than available in the
+    data instance raises.
+    """
     data = read("data/OkSmall.txt")
 
     assert_equal(data.num_vehicles, 3)
@@ -100,7 +115,7 @@ def test_route_constructor_raises_too_many_vehicles():
     with assert_raises(RuntimeError):
         Solution(data, [[1], [2], [3], [4]])
 
-    # Now test heterogeneous case
+    # Now test the case with multiple vehicle types.
     data = make_heterogeneous(data, [VehicleType(10, 2), VehicleType(20, 1)])
 
     # Only two routes (of type 0) should not raise.
@@ -128,6 +143,9 @@ def test_route_constructor_raises_too_many_vehicles():
 
 
 def test_route_constructor_raises_when_clients_are_visited_more_than_once():
+    """
+    Tests that visiting the same client more than once raises.
+    """
     data = read("data/OkSmall.txt")
     with assert_raises(RuntimeError):
         Solution(data, [[1, 2], [1, 3, 4]])  # client 1 is visited twice
@@ -137,6 +155,12 @@ def test_route_constructor_raises_when_clients_are_visited_more_than_once():
 
 
 def test_route_constructor_allows_incomplete_solutions():
+    """
+    Tests that not visiting a client at all is allowed, but turns the solution
+    infeasible (unless the client is not required). Allowing this is useful
+    for both prize collecting, and in LNS settings where an incomplete solution
+    is subsequently repaired.
+    """
     data = read("data/OkSmallPrizes.txt")
 
     # Client 1 is required but not visited.
@@ -156,11 +180,17 @@ def test_route_constructor_allows_incomplete_solutions():
 
 
 def test_get_neighbours():
+    """
+    Tests that accessing the neighbour structur of (pred, succ) pairs for each
+    client in the solution works correctly.
+    """
     data = read("data/OkSmall.txt")
+    assert_equal(data.num_clients, 4)
 
-    sol = Solution(data, [[3], [1, 2]])
-    neighbours = sol.get_neighbours()
+    sol = Solution(data, [[3], [1, 2]])  # client 4 is not in the solution
+    assert_(not sol.is_complete())
 
+    neighbours = sol.get_neighbours()
     expected = [
         None,  # 0: is depot
         (0, 2),  # 1: between depot (0) to 2
@@ -169,13 +199,15 @@ def test_get_neighbours():
         None,  # 4: unassigned
     ]
 
-    assert_equal(data.num_clients, 4)
-
     for client in range(data.num_clients + 1):  # incl. depot
         assert_equal(neighbours[client], expected[client])
 
 
 def test_feasibility():
+    """
+    Tests that solutions are infeasible when they have load or time window
+    violations.
+    """
     data = read("data/OkSmall.txt")
 
     # This solution is infeasible due to both load and time window violations.
@@ -197,6 +229,10 @@ def test_feasibility():
 
 
 def test_feasibility_release_times():
+    """
+    Tests solutions can be infeasible due to release time violations, which
+    adds time warp.
+    """
     data = read("data/OkSmallReleaseTimes.txt")
 
     # Client 1 is released at 20'000, but client 2 time window ends at 19'500,
@@ -217,6 +253,10 @@ def test_feasibility_release_times():
 
 
 def test_distance_calculation():
+    """
+    Tests that route distance calculations are correct, and that the overall
+    Solution's distance is the sum of the route distances.
+    """
     data = read("data/OkSmall.txt")
 
     sol = Solution(data, [[1, 2], [3], [4]])
@@ -241,6 +281,9 @@ def test_distance_calculation():
 
 
 def test_excess_load_calculation():
+    """
+    Tests the Solution's excess load calculation on a single-route case.
+    """
     data = read("data/OkSmall.txt")
 
     sol = Solution(data, [[4, 3, 1, 2]])
@@ -252,7 +295,11 @@ def test_excess_load_calculation():
     assert_allclose(sol.excess_load(), 18 - data.vehicle_type(0).capacity)
 
 
-def test_heterogeneous_capacity_excess_load_calculation():
+def test_excess_load_calculation_with_multiple_vehicle_capacities():
+    """
+    Tests that vehicles of different capacities result in different (excess)
+    load calaculations.
+    """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(
         data, vehicle_types=[VehicleType(10, 2), VehicleType(20, 1)]
@@ -272,6 +319,9 @@ def test_heterogeneous_capacity_excess_load_calculation():
 
 
 def test_route_access_methods():
+    """
+    Tests that accessing route statistics returns the correct numbers.
+    """
     data = read("data/OkSmall.txt")
     sol = Solution(data, [[1, 3], [2, 4]])
     routes = sol.get_routes()
@@ -306,6 +356,9 @@ def test_route_access_methods():
 
 
 def test_route_time_warp_calculations():
+    """
+    Tests route time warp calculations.
+    """
     data = read("data/OkSmall.txt")
     sol = Solution(data, [[1, 3], [2, 4]])
     routes = sol.get_routes()
@@ -327,6 +380,9 @@ def test_route_time_warp_calculations():
 
 
 def test_route_wait_time_calculations():
+    """
+    Tests route wait time and slack calculations.
+    """
     data = read("data/OkSmallWaitTime.txt")
     sol = Solution(data, [[1, 3], [2, 4]])
     routes = sol.get_routes()
@@ -372,6 +428,10 @@ def test_route_wait_time_calculations():
 
 
 def test_route_start_and_end_time_calculations():
+    """
+    Tests route start time, slack, and end time calculations for cases with
+    and without time warp or wait duration.
+    """
     data = read("data/OkSmall.txt")
     sol = Solution(data, [[1, 3], [2, 4]])
     routes = sol.get_routes()
@@ -414,6 +474,10 @@ def test_route_start_and_end_time_calculations():
 
 
 def test_route_release_time():
+    """
+    Tests that routes return the correct release times, and, when feasible,
+    start after the release time.
+    """
     data = read("data/OkSmallReleaseTimes.txt")
     sol = Solution(data, [[1, 3], [2, 4]])
     routes = sol.get_routes()
@@ -502,10 +566,12 @@ def test_time_warp_return_to_depot():
     assert_allclose(sol.time_warp(), 1)
 
 
-# TODO test all time warp cases
-
-
-def test_num_routes_calculation():
+def tests_that_not_specifying_the_vehicle_type_assumes_a_default():
+    """
+    Not specifying the vehicle type when providing a list of visits uses the
+    first vehicle type to complete the routes. That could result in a solution
+    using too many vehicles of the first type.
+    """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(
         data, vehicle_types=[VehicleType(10, 2), VehicleType(20, 1)]
@@ -525,14 +591,14 @@ def test_num_routes_calculation():
         # 3 routes of type 0 whereas we only have 2 available.
         sol = Solution(data, [[1], [2], [3, 4]])
 
-    # It works if we specify the correct vehicle types
+    # It works if we specify the correct vehicle types.
     sol = Solution(
         data,
         [Route(data, [1], 0), Route(data, [2], 0), Route(data, [3, 4], 1)],
     )
     assert_equal(sol.num_routes(), 3)
 
-    # But not if we violate the qty available per vehicle type
+    # But not if we violate the number of available vehicles of that type.
     with assert_raises(RuntimeError):
         sol = Solution(
             data,
@@ -541,22 +607,29 @@ def test_num_routes_calculation():
 
 
 def test_copy():
+    """
+    Tests that copied solutions are equal to the original solution, but not
+    the exact same object.
+    """
     data = read("data/OkSmall.txt")
 
     sol = Solution(data, [[1, 2, 3, 4]])
     copy_sol = copy(sol)
     deepcopy_sol = deepcopy(sol)
 
-    # Copied solutions are equal to the original solution
+    # Copied solutions are equal to the original solution.
     assert_(sol == copy_sol)
     assert_(sol == deepcopy_sol)
 
-    # But they are not the same object
+    # But they are not the same object.
     assert_(sol is not copy_sol)
     assert_(sol is not deepcopy_sol)
 
 
 def test_eq():
+    """
+    Tests the solution's equality operator.
+    """
     data = read("data/OkSmall.txt")
 
     sol1 = Solution(data, [[1, 2, 3, 4]])
@@ -581,7 +654,7 @@ def test_eq():
     assert_(sol5 != "cd")
 
 
-def test_eq_heterogeneous_vehicle():
+def test_eq_with_multiple_vehicle_types():
     """
     Tests that two solutions are not considered equal if they have the same
     routes (orders of clients) but served by different vehicle types.
@@ -639,8 +712,8 @@ def test_eq_unassigned():
 
 def test_duplicate_vehicle_types():
     """
-    Tests that even though it is not useful it is allowed to have duplicate
-    vehicle types which will be considered different.
+    Tests that it is allowed to have duplicate vehicle types. These will be
+    considered completely different during optimisation.
     """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(data, [VehicleType(10, 1), VehicleType(10, 1)])
@@ -656,6 +729,9 @@ def test_duplicate_vehicle_types():
     [[VehicleType(10, 3)], [VehicleType(10, 2), VehicleType(20, 1)]],
 )
 def test_str_contains_routes(vehicle_types):
+    """
+    Tests that the Solution's string representation contains each route.
+    """
     data = read("data/OkSmall.txt")
     data = make_heterogeneous(data, vehicle_types)
 
@@ -678,19 +754,19 @@ def test_str_contains_routes(vehicle_types):
 
 
 def test_hash():
+    """
+    Tests that solutions that compare the same have the same hash.
+    """
     data = read("data/OkSmall.txt")
     rng = RandomNumberGenerator(seed=2)
 
     sol1 = Solution.make_random(data, rng)
     sol2 = Solution.make_random(data, rng)
 
-    hash1 = hash(sol1)
-    hash2 = hash(sol2)
-
     # Two random solutions. They're not the same, so the hashes should not be
-    # the same either.
+    # the same either (unless there's a collision, which is not the case here).
     assert_(sol1 != sol2)
-    assert_(hash1 != hash2)
+    assert_(hash(sol1) != hash(sol2))
 
     sol3 = deepcopy(sol2)  # is a direct copy
 
@@ -700,6 +776,10 @@ def test_hash():
 
 
 def test_route_centroid():
+    """
+    Tests that each route's center point is the center point of all clients
+    visited by that route.
+    """
     data = read("data/OkSmall.txt")
     x = np.array([data.client(client).x for client in range(5)])
     y = np.array([data.client(client).y for client in range(5)])
@@ -713,6 +793,9 @@ def test_route_centroid():
 
 
 def test_solution_can_be_pickled():
+    """
+    Tests that a solution can be serialised and unserialised.
+    """
     data = read("data/OkSmall.txt")
     rng = RandomNumberGenerator(seed=2)
 
@@ -724,6 +807,9 @@ def test_solution_can_be_pickled():
 
 
 def test_route_can_be_pickled():
+    """
+    Tests that individual routes can be serialised and unserialised.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     rng = RandomNumberGenerator(seed=2)
     sol = Solution.make_random(data, rng)
diff --git a/pyvrp/tests/test_Statistics.py b/pyvrp/tests/test_Statistics.py
index 8b8cb62..100593f 100644
--- a/pyvrp/tests/test_Statistics.py
+++ b/pyvrp/tests/test_Statistics.py
@@ -13,6 +13,10 @@ from pyvrp.tests.helpers import read
 
 
 def test_csv_serialises_correctly(tmp_path):
+    """
+    Tests that writing a CSV of a ``Statistics`` object and then reading that
+    CSV again returns the same object.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -39,6 +43,10 @@ def test_csv_serialises_correctly(tmp_path):
 
 @pytest.mark.parametrize("num_iterations", [0, 5])
 def test_stats_collects_a_data_point_per_iteration(num_iterations: int):
+    """
+    Tests that the statistics object collects on feasible and infeasible data
+    point every time ``collect_from`` is called.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -55,6 +63,9 @@ def test_stats_collects_a_data_point_per_iteration(num_iterations: int):
 
 @pytest.mark.parametrize("num_iterations", [0, 1, 10])
 def test_eq(num_iterations: int):
+    """
+    Tests the equality operator.
+    """
     data = read("data/OkSmall.txt")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=42)
@@ -81,6 +92,9 @@ def test_eq(num_iterations: int):
 
 
 def test_more_eq():
+    """
+    Tests the equality operator for the same population trajectory.
+    """
     cost_evaluator = CostEvaluator(20, 6)
     pop = Population(broken_pairs_distance)
     assert_equal(len(pop), 0)
diff --git a/pyvrp/tests/test_SubPopulation.py b/pyvrp/tests/test_SubPopulation.py
index cde3dcd..ce55e76 100644
--- a/pyvrp/tests/test_SubPopulation.py
+++ b/pyvrp/tests/test_SubPopulation.py
@@ -15,6 +15,10 @@ from pyvrp.tests.helpers import read
 
 @mark.parametrize("nb_close", [5, 10, 25])
 def test_avg_distance_closest_is_same_up_to_nb_close(nb_close: int):
+    """
+    Tests that the average distance of a solution to other solutions only looks
+    at the nearest ``nb_close`` solutions.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=5)
@@ -49,6 +53,10 @@ def test_avg_distance_closest_is_same_up_to_nb_close(nb_close: int):
 
 
 def test_avg_distance_closest_for_single_route_solutions():
+    """
+    Tests that the closest computations are exactly right for a simple, single
+    route solution where it's easy to reason about what's going on.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     params = PopulationParams(min_pop_size=0, nb_close=10)
@@ -82,6 +90,10 @@ def test_avg_distance_closest_for_single_route_solutions():
 
 
 def test_fitness_is_purely_based_on_cost_when_only_elites():
+    """
+    Tests than when all solutions are considered elite, the fitness values
+    are completely determines by the solutions' cost.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=51)
@@ -114,6 +126,10 @@ def test_fitness_is_purely_based_on_cost_when_only_elites():
 
 
 def test_fitness_is_average_of_cost_and_diversity_when_no_elites():
+    """
+    When there are no elite solutions, the fitness ranking averages the cost
+    and diversity rank.
+    """
     data = read("data/RC208.txt", "solomon", "dimacs")
     cost_evaluator = CostEvaluator(20, 6)
     rng = RandomNumberGenerator(seed=52)
diff --git a/pyvrp/tests/test_TimeWindowSegment.py b/pyvrp/tests/test_TimeWindowSegment.py
index 7913998..5ab2eab 100644
--- a/pyvrp/tests/test_TimeWindowSegment.py
+++ b/pyvrp/tests/test_TimeWindowSegment.py
@@ -1,17 +1,24 @@
 import numpy as np
-from numpy.testing import assert_equal
+from numpy.testing import assert_, assert_equal
 from pytest import mark
 
 from pyvrp._pyvrp import TimeWindowSegment
 
 
 @mark.parametrize("existing_time_warp", [2, 5, 10])
-def test_total_time_warp(existing_time_warp):
+def test_total_time_warp_when_there_is_existing_time_warp(existing_time_warp):
+    """
+    Tests that the ``total_time_warp()`` returns existing time warp when no
+    segments have been merged yet.
+    """
     tws1 = TimeWindowSegment(0, 0, 0, existing_time_warp, 0, 0, 0)
     assert_equal(tws1.total_time_warp(), existing_time_warp)
 
 
 def test_merge_two():
+    """
+    Tests merging two time window segments.
+    """
     tws1 = TimeWindowSegment(0, 0, 5, 0, 0, 5, 0)
     tws2 = TimeWindowSegment(1, 1, 0, 5, 3, 6, 0)
 
@@ -33,7 +40,10 @@ def test_merge_two():
     assert_equal(merged.total_time_warp(), 11)
 
 
-def test_merge_multiple():
+def test_merge_three():
+    """
+    Tests merging three time window segments.
+    """
     tws1 = TimeWindowSegment(0, 0, 5, 0, 0, 5, 0)
     tws2 = TimeWindowSegment(1, 1, 0, 0, 3, 6, 0)
     tws3 = TimeWindowSegment(2, 2, 0, 0, 2, 3, 2)
@@ -52,5 +62,45 @@ def test_merge_multiple():
     assert_equal(merged3.total_time_warp(), 12)
 
 
-# TODO test two previously merged segments, e.g. merge(merge(1, 2),
-#  merge(3, 4)), each with time warp
+def test_merging_two_previously_merged_tws():
+    """
+    This test evaluates what happens when we merge two previously merged
+    segments, when both have time warp.
+    """
+    time_warp = 1
+    tws1 = TimeWindowSegment(0, 0, 5, time_warp, 0, 5, 0)  # depot
+    tws2 = TimeWindowSegment(1, 1, 1, time_warp, 3, 6, 0)  # client #1
+
+    # Each of these segments has some initial time warp.
+    assert_equal(tws1.total_time_warp(), 1)
+    assert_equal(tws2.total_time_warp(), 1)
+
+    mat = np.asarray([[0, 4], [3, 0]])
+    merged12 = TimeWindowSegment.merge(mat, tws1, tws2)  # depot -> client
+    merged21 = TimeWindowSegment.merge(mat, tws2, tws1)  # client -> depot
+
+    # The order in which time window segments are merged matters, so although
+    # these two segments cover the same clients, the total time warp is not
+    # the same.
+    assert_(merged12.total_time_warp() != merged21.total_time_warp())
+
+    # Both segments start with 1 initial time warp. Going from depot -> client
+    # happens at t = 5 - 1, and takes 4 time units to complete. So we arrive at
+    # t = 8, which is three units after the time window closes. This adds 2
+    # time warp to the segment, which, together with the initial time warps
+    # makes for 2 + 1 + 1 = 4 total time warp.
+    assert_equal(merged12.total_time_warp(), 4)
+
+    # We can leave at the earliest at t = 3, the start of the client's time
+    # window. We then arrive at t = 6, which adds one unit of time warp.
+    # Combined with the initial time warp, this results in 1 + 1 + 1 = 3 units
+    # of total time warp.
+    assert_equal(merged21.total_time_warp(), 3)
+
+    # This merged TWS represents the route plan 0 -> 1 -> 1 -> 0. We leave 1
+    # at t = 10 - 4, and travel takes no time. We do get the 1 unit of
+    # existing time warp. We then travel to the depot, arriving at t = 9. This
+    # is 4 time units after the time window closes, which adds 4 time warp
+    # (plus the existing unit). So we get 4 + 1 + 1 + 4 = 10 time warp.
+    merged = TimeWindowSegment.merge(mat, merged12, merged21)
+    assert_equal(merged.total_time_warp(), 10)
diff --git a/pyvrp/tests/test_read.py b/pyvrp/tests/test_read.py
index 1343473..bab31d6 100644
--- a/pyvrp/tests/test_read.py
+++ b/pyvrp/tests/test_read.py
@@ -1,3 +1,5 @@
+from math import sqrt
+
 from numpy.testing import assert_equal, assert_raises, assert_warns
 from pytest import mark
 
@@ -20,11 +22,17 @@ from pyvrp.tests.helpers import read
     ],
 )
 def test_raises_invalid_file(where: str, exception: Exception):
+    """
+    Tests that ``read()`` raises when there are issues with the given file.
+    """
     with assert_raises(exception):
         read(where)
 
 
 def test_raises_unknown_round_func():
+    """
+    Tests that ``read()`` raises when the rounding function is not known.
+    """
     with assert_raises(TypeError):
         # Unknown round_func, so should raise.
         read("data/OkSmall.txt", round_func="asdbsadfas")
@@ -34,6 +42,9 @@ def test_raises_unknown_round_func():
 
 
 def test_reading_OkSmall_instance():
+    """
+    Tests that the parsed data from the "OkSmall" instance is correct.
+    """
     data = read("data/OkSmall.txt")
 
     # From the DIMENSION, VEHICLES, and CAPACITY fields in the file.
@@ -98,6 +109,9 @@ def test_reading_OkSmall_instance():
 
 
 def test_reading_En22k4_instance():  # instance from CVRPLIB
+    """
+    Tests that the small E-n22-k4 instance from CVRPLIB is correctly parsed.
+    """
     data = read("data/E-n22-k4.txt", round_func="trunc1")
 
     assert_equal(data.num_clients, 21)
@@ -132,6 +146,9 @@ def test_reading_En22k4_instance():  # instance from CVRPLIB
 
 
 def test_reading_RC208_instance():  # Solomon style instance
+    """
+    Tests that a Solomon-style VRPTW instance is correctly parsed.
+    """
     data = read(
         "data/RC208.txt", instance_format="solomon", round_func="trunc1"
     )
@@ -174,13 +191,51 @@ def test_reading_RC208_instance():  # Solomon style instance
         assert_equal(data.client(client).required, True)
 
 
-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.
+def test_warns_about_scaling_issues():
+    """
+    Tests that ``read()`` warns about scaling issues when a distance value is
+    very large.
+    """
     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
+def test_round_func_trunc1_and_dimacs_are_same():
+    """
+    Tests that the DIMACS convention is equivalent to truncating to the first
+    decimal.
+    """
+    trunc1 = read("data/RC208.txt", "solomon", "trunc1")
+    dimacs = read("data/RC208.txt", "solomon", "dimacs")
+
+    trunc1_dist = trunc1.distance_matrix()
+    dimacs_dist = dimacs.distance_matrix()
+    assert_equal(trunc1_dist, dimacs_dist)
+
+    trunc1_dur = trunc1.duration_matrix()
+    dimacs_dur = trunc1.duration_matrix()
+    assert_equal(trunc1_dur, dimacs_dur)
+
+
+def test_round_func_round_nearest():
+    """
+    Tests rounding to the nearest integer works well for the RC208 instance,
+    which has Euclidean distances computed from integer coordinates. Since the
+    instance is large, we'll test one particular distance.
+    """
+    data = read("data/RC208.txt", "solomon", "round")
+
+    # We're going to test dist(0, 1) and dist(1, 0), which should be the same
+    # since the distances are symmetric/Euclidean.
+    assert_equal(data.client(0).x, 40)
+    assert_equal(data.client(0).y, 50)
+
+    assert_equal(data.client(1).x, 25)
+    assert_equal(data.client(1).y, 85)
+
+    # Compute the distance, and assert that it is indeed correctly rounded.
+    dist = sqrt((40 - 25) ** 2 + (85 - 50) ** 2)
+    assert_equal(data.dist(0, 1), round(dist))
+    assert_equal(data.dist(1, 0), round(dist))