From de573c5ecbcd66247741369a352eef88abaf6353 Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 13:41:50 +0100 Subject: [PATCH 01/10] Setup for primal integral --- bench/solve.py | 41 ++++++++++++++++++++++++++++++++++------- pyproject.toml | 6 +----- 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index 78a11f7..76ffb9c 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -79,6 +79,9 @@ class SolveResult(NamedTuple): gap The gap to the best-known solution if there is one, otherwise ``float('nan')``. + primal_integral + The primal integral of the solver run if a best-known solution is + provided and statistics are collected. Otherwise, ``float('nan')``. """ instance: str @@ -87,6 +90,7 @@ class SolveResult(NamedTuple): num_iterations: int runtime: float gap: float + primal_integral: float def _solve( @@ -189,12 +193,24 @@ def _solve( write_solution(sol_dir / (instance_name + ".sol"), data, result) gap = float("nan") + primal_integral = float("nan") + if bks_loc: sol = read_solution(bks_loc, data) cost_eval = CostEvaluator([0] * data.num_load_dimensions, 0, 0) bks = cost_eval.cost(sol) gap = 100 * (result.cost() - bks) / bks + # TODO COMPUTE PRIMAL INTEGRAL + stats = result.stats + if stats.is_collecting: + # [datum.best_cost for datum in stats.data] + # [datum.best_feas for datum in stats.data] + + # np.trapz over de stats, gedeeld door runtime * bks value, -1 + area = np.trapezoid(stats.best_costs, stats.runtimes) + primal_integral = area / (result.runtime * bks) - 1 + return SolveResult( instance_name, "Y" if result.is_feasible() else "N", @@ -202,6 +218,7 @@ def _solve( result.num_iterations, round(result.runtime, 3), round(gap, 2), + round(primal_integral, 3), ) @@ -246,15 +263,24 @@ def benchmark( ("iters", int), ("time", float), ("gap", float), + ("pi", float), ] data = np.asarray(res, dtype=dtypes) - headers = ["Instance", "OK", "Obj.", "Iters. (#)", "Time (s)", "Gap (%)"] + headers = [ + "Instance", + "OK", + "Obj.", + "Iters. (#)", + "Time (s)", + "Gap (%)", + "Primal Int.", + ] - exclude_gap = solutions is None - if exclude_gap: + exclude_headers = solutions is None + if exclude_headers: data = data[["inst", "ok", "obj", "iters", "time"]] - headers = headers[:-1] + headers = headers[:-2] print("\n", tabulate(headers, data), "\n", sep="") print(f" Avg. objective: {data['obj'].mean():.0f}") @@ -262,8 +288,9 @@ def benchmark( print(f" Avg. run-time: {data['time'].mean():.2f}s") print(f" Total not OK: {np.count_nonzero(data['ok'] == 'N')}") - if not exclude_gap: + if not exclude_headers: print(f" Avg. gap: {data['gap'].mean():.2f}%") + print(f" Avg. primal int.: {data['pi'].mean():.3f}") def setup_parser(subparser): @@ -280,8 +307,8 @@ def setup_parser(subparser): msg = """ Optional paths to best-known solutions in VRPLIB format, used to calculate - gaps. If provided, it must match the number of instances. Instances and - solutions are paired in the given order. + gaps and primal integrals. If provided, it must match the number of + instances. Instances and solutions are paired in the given order. """ parser.add_argument("--solutions", nargs="+", type=Path, help=msg) diff --git a/pyproject.toml b/pyproject.toml index 778d2cf..8445459 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,11 +13,7 @@ readme = "README.md" [tool.poetry.dependencies] python = ">=3.10,<4.0" -numpy = [ - # Numpy 1.26 is the first version of numpy that supports Python 3.12+. - { version = ">=1.15.2", python = "<3.12" }, - { version = ">=1.26.0", python = ">=3.12" } -] +numpy = ">=2.0.0" tqdm = "^4.64.1" tomli = "^2.0.1" From 3e88e703168e907e6c190c8fb73c8f9c999bf224 Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 16:25:21 +0100 Subject: [PATCH 02/10] Think this works? --- bench/solve.py | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index 76ffb9c..495998d 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -201,15 +201,29 @@ def _solve( bks = cost_eval.cost(sol) gap = 100 * (result.cost() - bks) / bks - # TODO COMPUTE PRIMAL INTEGRAL stats = result.stats - if stats.is_collecting: - # [datum.best_cost for datum in stats.data] - # [datum.best_feas for datum in stats.data] + primal_integral = 1 + if stats.is_collecting and result.is_feasible(): + first_feas = np.argmax([d.best_feas for d in stats.data]) - # np.trapz over de stats, gedeeld door runtime * bks value, -1 - area = np.trapezoid(stats.best_costs, stats.runtimes) - primal_integral = area / (result.runtime * bks) - 1 + # Start the primal integral from the first feasible solution. + best_costs = np.array([d.best_cost for d in stats.data]) + best_costs = best_costs[first_feas:] + + idcs = np.flatnonzero(np.diff(best_costs, prepend=0) != 0) + incumbents = best_costs[idcs] + + # Calculate gaps. + denom = np.maximum(np.abs(incumbents), abs(bks)) + gaps = abs(bks - incumbents) / denom + gaps = np.append(1, gaps) # gap is 100% at time 0 + + runtimes = np.cumsum(stats.runtimes)[first_feas:] + times = runtimes[idcs] + times = np.append(times, result.runtime) + times = np.diff(times, prepend=0) + + primal_integral = sum(times * gaps) / result.runtime * 100 return SolveResult( instance_name, @@ -218,7 +232,7 @@ def _solve( result.num_iterations, round(result.runtime, 3), round(gap, 2), - round(primal_integral, 3), + round(primal_integral, 2), ) @@ -274,7 +288,7 @@ def benchmark( "Iters. (#)", "Time (s)", "Gap (%)", - "Primal Int.", + "PI (%)", ] exclude_headers = solutions is None @@ -290,7 +304,7 @@ def benchmark( if not exclude_headers: print(f" Avg. gap: {data['gap'].mean():.2f}%") - print(f" Avg. primal int.: {data['pi'].mean():.3f}") + print(f" Avg. PI.: {data['pi'].mean():.2f}%") def setup_parser(subparser): From 3fba99dbe59ed6a3a0460eae1fae74199a47f707 Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:02:32 +0100 Subject: [PATCH 03/10] Simplified --- bench/solve.py | 38 +++++++++++++++++--------------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index 495998d..d1a458e 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -202,27 +202,23 @@ def _solve( gap = 100 * (result.cost() - bks) / bks stats = result.stats - primal_integral = 1 - if stats.is_collecting and result.is_feasible(): - first_feas = np.argmax([d.best_feas for d in stats.data]) - - # Start the primal integral from the first feasible solution. - best_costs = np.array([d.best_cost for d in stats.data]) - best_costs = best_costs[first_feas:] - - idcs = np.flatnonzero(np.diff(best_costs, prepend=0) != 0) - incumbents = best_costs[idcs] - - # Calculate gaps. - denom = np.maximum(np.abs(incumbents), abs(bks)) - gaps = abs(bks - incumbents) / denom - gaps = np.append(1, gaps) # gap is 100% at time 0 - - runtimes = np.cumsum(stats.runtimes)[first_feas:] - times = runtimes[idcs] - times = np.append(times, result.runtime) - times = np.diff(times, prepend=0) - + primal_integral = 100 + if stats.data and result.is_feasible(): + # Filter all cumulative runtimes and best costs from the first + # feasible solution. + feas = next(idx for idx, d in enumerate(stats.data) if d.best_feas) + runtimes = np.cumsum(stats.runtimes)[feas:] + best_costs = np.array([d.best_cost for d in stats.data[feas:]]) + + # Store the cost for each incumbent, along with when it was found. + _, idcs = np.unique(best_costs, return_index=True) + idcs.sort() # in decreasing order of costs + costs = best_costs[idcs] + found_at = runtimes[idcs] + + gaps = abs(bks - costs) / np.maximum(abs(costs), abs(bks)) + gaps = np.concatenate([[1], gaps]) + times = np.diff([0, *found_at, result.runtime]) primal_integral = sum(times * gaps) / result.runtime * 100 return SolveResult( From 1081a553ca5455b475b05cc339d2fb900533a75e Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:25:12 +0100 Subject: [PATCH 04/10] Simplify --- bench/solve.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index d1a458e..0e65468 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -82,6 +82,13 @@ class SolveResult(NamedTuple): primal_integral The primal integral of the solver run if a best-known solution is provided and statistics are collected. Otherwise, ``float('nan')``. + See [1]_ for details. + + References + ---------- + .. [1] Berthold, T. (2013). Measuring the impact of primal heuristics. + *Operations Research Letters*, 41(6): 611-614. + https://doi.org/10.1016/j.orl.2013.08.007. """ instance: str @@ -204,22 +211,16 @@ def _solve( stats = result.stats primal_integral = 100 if stats.data and result.is_feasible(): - # Filter all cumulative runtimes and best costs from the first - # feasible solution. + # Compute primal gaps starting from the first feasible solution. feas = next(idx for idx, d in enumerate(stats.data) if d.best_feas) - runtimes = np.cumsum(stats.runtimes)[feas:] - best_costs = np.array([d.best_cost for d in stats.data[feas:]]) - - # Store the cost for each incumbent, along with when it was found. - _, idcs = np.unique(best_costs, return_index=True) - idcs.sort() # in decreasing order of costs - costs = best_costs[idcs] - found_at = runtimes[idcs] - + costs = np.array([d.best_cost for d in stats.data[feas:]]) gaps = abs(bks - costs) / np.maximum(abs(costs), abs(bks)) - gaps = np.concatenate([[1], gaps]) - times = np.diff([0, *found_at, result.runtime]) - primal_integral = sum(times * gaps) / result.runtime * 100 + gaps = np.concatenate([[1], gaps]) # primal gap is 1 at time 0 + + # Do the same for time intervals between two solutions. + runtimes = np.cumsum(stats.runtimes)[feas:] + times = np.diff(prepend=0, a=runtimes, append=result.runtime) + primal_integral = sum(gaps * times) / result.runtime * 100 return SolveResult( instance_name, From 47907ba39b1b7cf4bbaff4138ac6bac4faa2b62e Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:30:22 +0100 Subject: [PATCH 05/10] Fix small alignment --- bench/solve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bench/solve.py b/bench/solve.py index 0e65468..57b405f 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -301,7 +301,7 @@ def benchmark( if not exclude_headers: print(f" Avg. gap: {data['gap'].mean():.2f}%") - print(f" Avg. PI.: {data['pi'].mean():.2f}%") + print(f" Avg. PI: {data['pi'].mean():.2f}%") def setup_parser(subparser): From a14e96a238d7a1cce4af4341e5909183a2a7462f Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:40:43 +0100 Subject: [PATCH 06/10] Rephrase --- bench/solve.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index 57b405f..ee12782 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -211,15 +211,16 @@ def _solve( stats = result.stats primal_integral = 100 if stats.data and result.is_feasible(): - # Compute primal gaps starting from the first feasible solution. + # Compute primal gap from first feasible solution onward. feas = next(idx for idx, d in enumerate(stats.data) if d.best_feas) costs = np.array([d.best_cost for d in stats.data[feas:]]) gaps = abs(bks - costs) / np.maximum(abs(costs), abs(bks)) - gaps = np.concatenate([[1], gaps]) # primal gap is 1 at time 0 + gaps = np.concatenate([[1], gaps]) # gap is 1 till first incumbent - # Do the same for time intervals between two solutions. - runtimes = np.cumsum(stats.runtimes)[feas:] - times = np.diff(prepend=0, a=runtimes, append=result.runtime) + # Time intervals between consecutive incumbent solutions, including + # border cases for the first and last found solutions. + found_at = np.cumsum(stats.runtimes)[feas:] + times = np.diff(prepend=0, a=found_at, append=result.runtime) primal_integral = sum(gaps * times) / result.runtime * 100 return SolveResult( From 4bc3ed36229b95ae44b167c8595c6a9e4d77e75a Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:47:47 +0100 Subject: [PATCH 07/10] Revert numpy changes --- pyproject.toml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8445459..05716af 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,8 +12,12 @@ readme = "README.md" [tool.poetry.dependencies] -python = ">=3.10,<4.0" -numpy = ">=2.0.0" +python = ">=3.11" +numpy = [ + # Numpy 1.26 is the first version of numpy that supports Python 3.12+. + { version = ">=1.15.2", python = "<3.12" }, + { version = ">=1.26.0", python = ">=3.12" } +] tqdm = "^4.64.1" tomli = "^2.0.1" From 47b2b41f30d7012e07a66a1b906d8567c1eef0e1 Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Mon, 29 Dec 2025 17:51:42 +0100 Subject: [PATCH 08/10] Rephrase --- bench/solve.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index ee12782..e6a1a4c 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -211,14 +211,14 @@ def _solve( stats = result.stats primal_integral = 100 if stats.data and result.is_feasible(): - # Compute primal gap from first feasible solution onward. + # Compute primal gap from the first feasible solution onward. feas = next(idx for idx, d in enumerate(stats.data) if d.best_feas) costs = np.array([d.best_cost for d in stats.data[feas:]]) gaps = abs(bks - costs) / np.maximum(abs(costs), abs(bks)) gaps = np.concatenate([[1], gaps]) # gap is 1 till first incumbent - # Time intervals between consecutive incumbent solutions, including - # border cases for the first and last found solutions. + # Time intervals between incumbent solutions (including start and + # end of the solver progress). found_at = np.cumsum(stats.runtimes)[feas:] times = np.diff(prepend=0, a=found_at, append=result.runtime) primal_integral = sum(gaps * times) / result.runtime * 100 From b5bba77999c18b9ab6a7584f998b0d2b9ad5b8ea Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Tue, 30 Dec 2025 10:52:24 +0100 Subject: [PATCH 09/10] Simplify PI --- bench/solve.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/bench/solve.py b/bench/solve.py index e6a1a4c..00379f2 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -60,6 +60,20 @@ def write_solution(where: Path, data, result): fh.write(f"Cost: {round(result.cost(), 2)}\n") +def pi(stats, bks_value: int) -> float: + if len(stats.data) == 0: + return 100 + + bks_value = min(bks_value, stats.data[-1].best_cost) + best_values = np.array([datum.best_cost for datum in stats], dtype=float) + gaps = (best_values - bks_value) / best_values + + is_feas = np.array([datum.best_feas for datum in stats], dtype=bool) + gaps[~is_feas] = 1 + + return 100 * np.sum(gaps * stats.runtimes) / sum(stats.runtimes) + + class SolveResult(NamedTuple): """ Named tuple to store the results of a single solver run. @@ -207,21 +221,7 @@ def _solve( cost_eval = CostEvaluator([0] * data.num_load_dimensions, 0, 0) bks = cost_eval.cost(sol) gap = 100 * (result.cost() - bks) / bks - - stats = result.stats - primal_integral = 100 - if stats.data and result.is_feasible(): - # Compute primal gap from the first feasible solution onward. - feas = next(idx for idx, d in enumerate(stats.data) if d.best_feas) - costs = np.array([d.best_cost for d in stats.data[feas:]]) - gaps = abs(bks - costs) / np.maximum(abs(costs), abs(bks)) - gaps = np.concatenate([[1], gaps]) # gap is 1 till first incumbent - - # Time intervals between incumbent solutions (including start and - # end of the solver progress). - found_at = np.cumsum(stats.runtimes)[feas:] - times = np.diff(prepend=0, a=found_at, append=result.runtime) - primal_integral = sum(gaps * times) / result.runtime * 100 + primal_integral = pi(result.stats, bks) return SolveResult( instance_name, From a32c09ebe9137b176db1c58376ddc42d6bfc60df Mon Sep 17 00:00:00 2001 From: Leon Lan Date: Tue, 30 Dec 2025 11:01:27 +0100 Subject: [PATCH 10/10] Add docstring --- bench/solve.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/bench/solve.py b/bench/solve.py index 00379f2..0d9e42c 100644 --- a/bench/solve.py +++ b/bench/solve.py @@ -61,6 +61,10 @@ def write_solution(where: Path, data, result): def pi(stats, bks_value: int) -> float: + """ + Computes the primal integral over the given statistics, using the provided + best-known solution value. + """ if len(stats.data) == 0: return 100