Repository: PyVRP/PyVRP (python)

Issue:
Output solver progress
#### Is your feature request related to a problem? Please describe
There is currently no real-time feedback about the solver progress. This leaves the user in the dark about how the solver is doing; the solver might have converged already, and the user won't know until they investigate the statistics. 

It's common for many solvers like Gurobi and CPLEX to output the solver progress to the console like this:

<details>
<summary>Example</summary>

```
 ! --------------------------------------------------- CP Optimizer 22.1.1.0 --
 ! Minimization problem - 15586 variables, 31679 constraints
 ! Presolve      : 7195 extractables eliminated
 ! TimeLimit            = 100
 ! LogVerbosity         = Terse
 ! SearchType           = IterativeDiving
 ! Initial process time : 0.18s (0.18s extraction + 0.00s propagation)
 !  . Log search space  : 38467.7 (before), 38467.7 (after)
 !  . Memory usage      : 62.9 MB (before), 62.9 MB (after)
 ! Using parallel search with 8 workers.
 ! ----------------------------------------------------------------------------
 !          Best Branches  Non-fixed    W       Branch decision
                        0      15586                 -
 + New bound is 9702
 *         83128    20007  23.24s       1      (gap is 88.33%)
           83128    36858        921    6            -
 + New bound is 39749 (gap is 52.18%)
 *         83058    43066  31.68s       6      (gap is 52.14%)
 *         74197    36598  32.25s       3      (gap is 46.43%)
 ```


</details>

We already have access to data like incumbents, solver time and iterations since we collect statistics by default.

#### Describe the solution you'd like
I would like to have `Model.solve` interface output the solver progress. The desired behavior is to set some argument like `log=True` to output the solver progress. I even think logging should be the default.

The simplest version should output the following each time a new best solution is found:
- Incumbent (best feasible solution)
- Time (elapsed solving time)
- Iterations (elapsed number of iterations)

We can also have more verbose variants, which output more statistics or output logs on different events (e.g., population restart).

#### Describe alternatives you have considered
N/A

Pull request description:
Print solver progress
This PR adds support for printing solver progress. This is enabled by default when solving through the Model interface, but not for GA itself. As such, this PR:

- [x] Closes #455.
- [x] Adds and passes relevant tests.
- [x] Has well-formatted code and documentation.


Patch:
diff --git a/docs/source/setup/benchmarks.rst b/docs/source/setup/benchmarks.rst
index e32693f..9d5123f 100644
--- a/docs/source/setup/benchmarks.rst
+++ b/docs/source/setup/benchmarks.rst
@@ -86,6 +86,12 @@ The reported values represent the average gaps to the best-known solutions (at t
      - VRPTW
      - PC-VRPTW
      - MDVRPTW
+   * - 11 February 2024
+     - `Google OR-Tools <https://developers.google.com/optimization/routing>`_ [#note3]_
+     - 5.23%
+     - 10.86%
+     - 13.24%
+     -
    * - 28 February 2023
      - `HGS-CVRP <https://github.com/vidalt/HGS-CVRP>`_
      - 0.11%
@@ -117,3 +123,8 @@ The reported values represent the average gaps to the best-known solutions (at t
    Literature result from Table 8 of `Vidal et al. (2013) <https://doi.org/10.1016/j.cor.2012.07.018>`_.
    This is an average gap over five seeds, rather than the usual ten.
    Note that this paper introduced the MDVRPTW benchmark instances.
+
+.. [#note3]
+   Results obtained using Google OR-Tools `v9.8.3296 <https://pypi.org/project/ortools/9.8.3296/>`_.
+   The reported average gap is over only one seed, because there is no interface to set the seed.
+   The code used to benchmark Google OR-tools can be found in `this issue <https://github.com/PyVRP/PyVRP/issues/453>`_.
diff --git a/examples/basic_vrps.ipynb b/examples/basic_vrps.ipynb
index b0d6105..63b6676 100644
--- a/examples/basic_vrps.ipynb
+++ b/examples/basic_vrps.ipynb
@@ -127,7 +127,7 @@
    "outputs": [],
    "source": [
     "model = Model.from_data(INSTANCE)\n",
-    "result = model.solve(stop=MaxIterations(2000), seed=42)\n",
+    "result = model.solve(stop=MaxIterations(2000), seed=42, display=False)\n",
     "print(result)"
    ]
   },
@@ -261,7 +261,7 @@
    "outputs": [],
    "source": [
     "model = Model.from_data(INSTANCE)\n",
-    "result = model.solve(stop=MaxIterations(1000), seed=42)\n",
+    "result = model.solve(stop=MaxIterations(1000), seed=42, display=False)\n",
     "print(result)"
    ]
   },
@@ -424,7 +424,7 @@
    "outputs": [],
    "source": [
     "model = Model.from_data(INSTANCE)\n",
-    "result = model.solve(stop=MaxRuntime(30), seed=42)"
+    "result = model.solve(stop=MaxRuntime(30), seed=42, display=False)"
    ]
   },
   {
diff --git a/examples/quick_tutorial.ipynb b/examples/quick_tutorial.ipynb
index d81e320..82e678d 100644
--- a/examples/quick_tutorial.ipynb
+++ b/examples/quick_tutorial.ipynb
@@ -114,7 +114,7 @@
    "metadata": {},
    "source": [
     "The instance looks good, so we are ready to solve it.\n",
-    "Let's do so with a second of runtime."
+    "Let's do so with a second of runtime, and display the search progress using the `display` argument on `Model.solve`."
    ]
   },
   {
@@ -125,7 +125,26 @@
    "source": [
     "from pyvrp.stop import MaxRuntime\n",
     "\n",
-    "res = m.solve(stop=MaxRuntime(1))  # one second\n",
+    "res = m.solve(stop=MaxRuntime(1), display=True)  # one second"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By passing the `display` argument, PyVRP displays statistics about the solver progress and the instance being solved.\n",
+    "In particular, it outputs the sizes of the feasible and infeasible solution pools, their average objective values, and the objective of the best solutions in either pool.\n",
+    "A heuristic improvement is indicated by a `H` at the start of a line.\n",
+    "\n",
+    "Let's print the solution we have found to see the routes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "print(res)"
    ]
   },
@@ -265,7 +284,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "res = m.solve(stop=MaxRuntime(1))  # one second\n",
+    "res = m.solve(stop=MaxRuntime(1), display=False)  # one second\n",
     "print(res)"
    ]
   },
@@ -398,7 +417,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "res = m.solve(stop=MaxRuntime(1))  # one second\n",
+    "res = m.solve(stop=MaxRuntime(1), display=False)  # one second\n",
     "print(res)"
    ]
   },
diff --git a/examples/using_pyvrp_components.ipynb b/examples/using_pyvrp_components.ipynb
index 9d73622..72bdca9 100644
--- a/examples/using_pyvrp_components.ipynb
+++ b/examples/using_pyvrp_components.ipynb
@@ -32,7 +32,7 @@
    "metadata": {},
    "source": [
     "We will implement a `solve()` function that will take a `stop` stopping criterion, and a `seed` for the random number generator.\n",
-    "This definition closely matches that of `Model.solve`.\n",
+    "This definition is very close to that of `Model.solve`.\n",
     "The signature is\n",
     "```python\n",
     "def solve(stop: StoppingCriterion, seed: int) -> Solution: ...\n",
diff --git a/pyvrp/GeneticAlgorithm.py b/pyvrp/GeneticAlgorithm.py
index bdeb506..7de2cab 100644
--- a/pyvrp/GeneticAlgorithm.py
+++ b/pyvrp/GeneticAlgorithm.py
@@ -4,6 +4,7 @@ import time
 from dataclasses import dataclass
 from typing import TYPE_CHECKING, Callable, Collection
 
+from pyvrp.ProgressPrinter import ProgressPrinter
 from pyvrp.Result import Result
 from pyvrp.Statistics import Statistics
 
@@ -128,7 +129,7 @@ class GeneticAlgorithm:
     def _cost_evaluator(self) -> CostEvaluator:
         return self._pm.get_cost_evaluator()
 
-    def run(self, stop: StoppingCriterion):
+    def run(self, stop: StoppingCriterion, display: bool = False):
         """
         Runs the genetic algorithm with the provided stopping criterion.
 
@@ -137,12 +138,18 @@ class GeneticAlgorithm:
         stop
             Stopping criterion to use. The algorithm runs until the first time
             the stopping criterion returns ``True``.
+        display
+            Whether to display information about the solver progress. Default
+            ``False``.
 
         Returns
         -------
         Result
             A Result object, containing statistics and the best found solution.
         """
+        print_progress = ProgressPrinter(should_print=display)
+        print_progress.start(self._data)
+
         start = time.perf_counter()
         stats = Statistics()
         iters = 0
@@ -155,6 +162,8 @@ class GeneticAlgorithm:
             iters += 1
 
             if iters_no_improvement == self._params.nb_iter_no_improvement:
+                print_progress.restart()
+
                 iters_no_improvement = 1
                 self._pop.clear()
 
@@ -177,9 +186,14 @@ class GeneticAlgorithm:
                 iters_no_improvement += 1
 
             stats.collect_from(self._pop, self._cost_evaluator)
+            print_progress.iteration(stats)
 
         end = time.perf_counter() - start
-        return Result(self._best, stats, iters, end)
+        res = Result(self._best, stats, iters, end)
+
+        print_progress.end(res)
+
+        return res
 
     def _improve_offspring(self, sol: Solution):
         def is_new_best(sol):
diff --git a/pyvrp/Model.py b/pyvrp/Model.py
index 6cca1af..d056205 100644
--- a/pyvrp/Model.py
+++ b/pyvrp/Model.py
@@ -267,7 +267,12 @@ class Model:
             durations,
         )
 
-    def solve(self, stop: StoppingCriterion, seed: int = 0) -> Result:
+    def solve(
+        self,
+        stop: StoppingCriterion,
+        seed: int = 0,
+        display: bool = True,
+    ) -> Result:
         """
         Solve this model.
 
@@ -276,7 +281,10 @@ class Model:
         stop
             Stopping criterion to use.
         seed
-            Seed value to use for the PRNG, by default 0.
+            Seed value to use for the random number stream. Default 0.
+        display
+            Whether to display information about the solver progress. Default
+            ``True``.
 
         Returns
         -------
@@ -315,4 +323,4 @@ class Model:
 
         gen_args = (data, pm, rng, pop, ls, crossover, init)
         algo = GeneticAlgorithm(*gen_args)  # type: ignore
-        return algo.run(stop)
+        return algo.run(stop, display)
diff --git a/pyvrp/ProgressPrinter.py b/pyvrp/ProgressPrinter.py
new file mode 100644
index 0000000..5ce0143
--- /dev/null
+++ b/pyvrp/ProgressPrinter.py
@@ -0,0 +1,125 @@
+from importlib.metadata import version
+
+from pyvrp._pyvrp import ProblemData
+
+from .Result import Result
+from .Statistics import Statistics
+
+# Templates for various different outputs.
+_ITERATION = (
+    "{special} {iters:>7} {elapsed:>6}s | "
+    "{feas_size:>3} {feas_avg:>8} {feas_best:>8} | "
+    "{infeas_size:>3} {infeas_avg:>8} {infeas_best:>8}"
+)
+
+_START = """PyVRP v{version}
+
+Solving an instance with:
+    {depot_text}
+    {client_text}
+    {vehicle_text} ({vehicle_type_text})
+
+                  |       Feasible        |      Infeasible
+    Iters    Time |   #      Avg     Best |   #      Avg     Best"""
+
+_END = """
+Search terminated in {runtime:.2f}s after {iters} iterations.
+Best-found solution has cost {best_cost}.
+"""
+
+_RESTART = "R                 |        restart        |        restart"
+
+
+class ProgressPrinter:
+    """
+    A helper class that prints relevant solver progress information to the
+    console, if desired.
+
+    Parameters
+    ----------
+    should_print
+        Whether to print information to the console. When ``False``, nothing is
+        printed.
+    """
+
+    def __init__(self, should_print: bool):
+        self._print = should_print
+        self._best_cost = float("inf")
+
+    def iteration(self, stats: Statistics):
+        """
+        Outputs relevant information every few hundred iterations. The output
+        contains information about the feasible and infeasible populations,
+        whether a new best solution has been found, and the search duration.
+        """
+        if not self._print or stats.num_iterations % 500 != 0:
+            return
+
+        feas = stats.feas_stats[-1]
+        infeas = stats.infeas_stats[-1]
+
+        msg = _ITERATION.format(
+            special="H" if feas.best_cost < self._best_cost else " ",
+            iters=stats.num_iterations,
+            elapsed=round(sum(stats.runtimes)),
+            feas_size=feas.size,
+            feas_avg=round(feas.avg_cost),
+            feas_best=round(feas.best_cost),
+            infeas_size=infeas.size,
+            infeas_avg=round(infeas.avg_cost),
+            infeas_best=round(infeas.best_cost),
+        )
+        print(msg)
+
+        if feas.best_cost < self._best_cost:
+            self._best_cost = feas.best_cost
+
+    def start(self, data: ProblemData):
+        """
+        Outputs information about PyVRP and the data instance that is being
+        solved.
+        """
+        if not self._print:
+            return
+
+        num_d = data.num_depots
+        depot_text = f"{num_d} depot{'s' if num_d > 1 else ''}"
+
+        num_c = data.num_clients
+        client_text = f"{num_c} client{'s' if num_c > 1 else ''}"
+
+        num_v = data.num_vehicles
+        vehicle_text = f"{num_v} vehicle{'s' if num_v > 1 else ''}"
+
+        num_vt = data.num_vehicle_types
+        vehicle_type_text = f"{num_vt} vehicle type{'s' if num_vt > 1 else ''}"
+
+        msg = _START.format(
+            version=version("pyvrp"),
+            depot_text=depot_text,
+            client_text=client_text,
+            vehicle_text=vehicle_text,
+            vehicle_type_text=vehicle_type_text,
+        )
+        print(msg)
+
+    def end(self, result: Result):
+        """
+        Outputs information about the search duration and the best-found
+        solution.
+        """
+        if self._print:
+            msg = _END.format(
+                iters=result.num_iterations,
+                runtime=result.runtime,
+                best_cost=round(result.cost(), 2),
+            )
+            print(msg)
+
+    def restart(self):
+        """
+        Indicates in the progress information that the algorithm has restarted
+        the search.
+        """
+        if self._print:
+            print(_RESTART)
diff --git a/pyvrp/cli.py b/pyvrp/cli.py
index ecf69e3..5ad5371 100644
--- a/pyvrp/cli.py
+++ b/pyvrp/cli.py
@@ -237,10 +237,10 @@ def benchmark(instances: list[Path], num_procs: int = 1, **kwargs):
     headers = ["Instance", "OK", "Obj.", "Iters. (#)", "Time (s)"]
 
     print("\n", tabulate(headers, data), "\n", sep="")
-    print(f"      Avg. objective: {data['obj'].mean():.0f}")
-    print(f"     Avg. iterations: {data['iters'].mean():.0f}")
-    print(f"   Avg. run-time (s): {data['time'].mean():.2f}")
-    print(f"        Total not OK: {np.count_nonzero(data['ok'] == 'N')}")
+    print(f"     Avg. objective: {data['obj'].mean():.0f}")
+    print(f"    Avg. iterations: {data['iters'].mean():.0f}")
+    print(f"      Avg. run-time: {data['time'].mean():.2f}s")
+    print(f"       Total not OK: {np.count_nonzero(data['ok'] == 'N')}")
 
 
 def main():


Test patch:
diff --git a/tests/test_Model.py b/tests/test_Model.py
index 8dee084..0ec9bd0 100644
--- a/tests/test_Model.py
+++ b/tests/test_Model.py
@@ -306,6 +306,35 @@ def test_model_and_solve(ok_small):
     assert_equal(res.cost(), 9_155)
 
 
+def test_model_solve_display_argument(ok_small, capsys):
+    """
+    Tests that solving a model displays solver progress when the ``display``
+    argument is ``True``.
+    """
+    model = Model.from_data(ok_small)
+
+    # First solve with display turned off. We should not see any output in this
+    # case.
+    model.solve(stop=MaxIterations(1000), seed=0, display=False)
+    printed = capsys.readouterr().out
+    assert_equal(printed, "")
+
+    # Now solve with display turned on. We should see output now.
+    res = model.solve(stop=MaxIterations(1000), seed=0, display=True)
+    printed = capsys.readouterr().out
+
+    # Check that some of the header data is in the output.
+    assert_("PyVRP" in printed)
+    assert_("Time" in printed)
+    assert_("Iters" in printed)
+    assert_("Feasible" in printed)
+    assert_("Infeasible" in printed)
+
+    # Check that we include the cost and total runtime in the output somewhere.
+    assert_(str(round(res.cost())) in printed)
+    assert_(str(round(res.runtime)) in printed)
+
+
 def test_partial_distance_duration_matrix():
     """
     Tests that adding a full distance or duration matrix is not required. Any
diff --git a/tests/test_ProgressPrinter.py b/tests/test_ProgressPrinter.py
new file mode 100644
index 0000000..cb04464
--- /dev/null
+++ b/tests/test_ProgressPrinter.py
@@ -0,0 +1,149 @@
+from numpy.testing import assert_, assert_equal
+
+from pyvrp import (
+    CostEvaluator,
+    Population,
+    RandomNumberGenerator,
+    Result,
+    Solution,
+    Statistics,
+)
+from pyvrp.ProgressPrinter import ProgressPrinter
+from pyvrp.diversity import broken_pairs_distance as bpd
+
+
+def test_start(ok_small, capsys):
+    """
+    Tests that the progress printer outputs some statistics about the given
+    instance when calling start.
+    """
+    printer = ProgressPrinter(should_print=True)
+    printer.start(ok_small)
+
+    out = capsys.readouterr().out
+    assert_(f"{ok_small.num_depots} depot" in out)
+    assert_(f"{ok_small.num_clients} clients" in out)
+    assert_(f"{ok_small.num_vehicles} vehicles" in out)
+
+
+def test_start_multiple_depot_plural(ok_small_multi_depot, capsys):
+    """
+    Tests that "depots" is plural when there are multiple depots.
+    """
+    printer = ProgressPrinter(should_print=True)
+    printer.start(ok_small_multi_depot)
+
+    out = capsys.readouterr().out
+    assert_(f"{ok_small_multi_depot.num_depots} depots" in out)
+
+
+def test_end(ok_small, capsys):
+    """
+    Tests that the progress printer outputs some summary statistics about the
+    best-found solution and solver run when calling end.
+    """
+    best = Solution(ok_small, [[1, 2], [3, 4]])
+    res = Result(best, Statistics(), 25, 0.05)
+
+    printer = ProgressPrinter(should_print=True)
+    printer.end(res)
+
+    out = capsys.readouterr().out
+    assert_(str(round(res.cost())) in out)
+    assert_(str(res.num_iterations) in out)
+    assert_(str(round(res.runtime)) in out)
+
+
+def test_restart(capsys):
+    """
+    Tests that calling restart outputs a line indicating such.
+    """
+    printer = ProgressPrinter(should_print=True)
+    printer.restart()
+
+    out = capsys.readouterr().out
+    assert_("R" in out)
+    assert_("restart" in out)
+
+
+def test_iteration(ok_small, capsys):
+    """
+    Tests that calling iteration prints a line about the current feasible and
+    infeasible populations, but only every 500 iterations.
+    """
+    pop = Population(bpd)
+    rng = RandomNumberGenerator(seed=42)
+    cost_eval = CostEvaluator(1, 1)
+
+    for _ in range(10):
+        pop.add(Solution.make_random(ok_small, rng), cost_eval)
+
+    assert_(pop.num_feasible() > 0)
+    assert_(pop.num_infeasible() > 0)
+
+    stats = Statistics()
+    stats.collect_from(pop, cost_eval)
+
+    printer = ProgressPrinter(should_print=True)
+
+    # Output is printed only every once in a while. We do not print after a
+    # single iteration.
+    printer.iteration(stats)
+    out = capsys.readouterr().out
+    assert_(stats.num_iterations, 1)
+    assert_equal(out, "")
+
+    # But only every 500 iterations. Thus, we should output results now.
+    stats.num_iterations = 500
+    printer.iteration(stats)
+    out = capsys.readouterr().out
+
+    # Statistics about solver progress.
+    assert_(str(stats.num_iterations) in out)
+    assert_(str(round(sum(stats.runtimes))) in out)
+
+    # Statistics about feasible population.
+    feas = stats.feas_stats[-1]
+    assert_(str(feas.size) in out)
+    assert_(str(round(feas.avg_cost)) in out)
+    assert_(str(round(feas.best_cost)) in out)
+
+    # Statistics about infeasible population.
+    infeas = stats.infeas_stats[-1]
+    assert_(str(infeas.size) in out)
+    assert_(str(round(infeas.avg_cost)) in out)
+    assert_(str(round(infeas.best_cost)) in out)
+
+
+def test_should_print_false_no_output(ok_small, capsys):
+    """
+    Tests that disabling printing works.
+    """
+    pop = Population(bpd)
+    rng = RandomNumberGenerator(seed=42)
+    cost_eval = CostEvaluator(1, 1)
+
+    for _ in range(10):
+        pop.add(Solution.make_random(ok_small, rng), cost_eval)
+
+    best = Solution(ok_small, [[1, 2], [3, 4]])  # BKS
+    pop.add(best, cost_eval)
+
+    # Checks that there are both feasible and infeasible solutions in this
+    # population.
+    assert_(pop.num_feasible() > 0)
+    assert_(pop.num_infeasible() > 0)
+
+    stats = Statistics()
+    stats.collect_from(pop, cost_eval)
+
+    # Set up the progress printer, call all its methods, and then check that
+    # the output captured on stdout is empty.
+    printer = ProgressPrinter(should_print=False)
+    printer.start(ok_small)
+    printer.iteration(stats=stats)
+    printer.restart()
+    printer.end(result=Result(best, stats, 25, 0.05))
+
+    out = capsys.readouterr().out
+    assert_equal(out, "")