Read the HiGHS best bound from LpSolveStats

README.md

@@ -80,7 +80,7 @@ A schedule the solver stopped on at the limit is reported as `Feasible` rather t
 
 ### Solver
 
-`OPTIMIZER_SOLVER` selects the MILP solver: `cbc` (default), the CBC binary pulp resolves from `PATH` or `OPTIMIZER_CBC_PATH`, the 2.10.10 release build in the image; `highs`, HiGHS run in process through `highspy`; or `hybrid`, CBC for the joint probe and HiGHS for the cost, tie break and continuity stages. HiGHS is asked to prove to the same absolute gap as CBC, its own relative gap of 1e-4 is switched off, and its model is built in bulk (`optimizer/highs.py`, coin-or/pulp#957). On HiGHS the continuity stage is interrupted once it runs `CUTOFF_MARGIN` past its limit, keeping its incumbent; CBC is left to its own time limit.
+`OPTIMIZER_SOLVER` selects the MILP solver: `cbc` (default), the CBC binary pulp resolves from `PATH` or `OPTIMIZER_CBC_PATH`, the 2.10.10 release build in the image; `highs`, HiGHS run in process through `highspy`; or `hybrid`, CBC for the joint probe and HiGHS for the cost, tie break and continuity stages. HiGHS is asked to prove to the same absolute gap as CBC, its own relative gap of 1e-4 is switched off. On HiGHS the continuity stage is interrupted once it runs `CUTOFF_MARGIN` past its limit, keeping its incumbent; CBC is left to its own time limit.
 
 Measured over 170 captured requests, 150 of them slow ones: HiGHS alone is the fastest in sum but proves 8 of the 54 mid-sized requests, which CBC solves in 1 to 5 s, 2.5 to 5 times slower, overruns the probe and pays a 6 to 8 s split for each, which lifted production p99 from 6.4 to 8.5 s. `hybrid` matches CBC on those (p90 2.05 against 2.01 s) and keeps HiGHS's quality on the slow ones (11 CBC Not Solved become 1), at the price of the 19 slow requests HiGHS proves joint in about 1.3 s, which CBC cannot prove in the probe and the hybrid splits instead; its total time is CBC's.
 

pyproject.toml

@@ -15,7 +15,7 @@ dependencies = [
     "matplotlib",
     "numpy>=2.5.4",
     "pandas>=2.3.3",
-    "pulp>=4.0.0",
+    "pulp>=4.0.1",
     "pydantic>=2.14.1",
     "pydantic-settings>=2.15.0",
     "pyjwt>=2.15.1",

src/optimizer/optimizer.py

@@ -1,5 +1,4 @@
 import os
-import re
 import threading
 import time
 from contextlib import contextmanager
@@ -12,7 +11,6 @@
 import numpy as np
 import pulp
 
-from .highs import BulkHiGHS
 from .settings import OptimizerSettings
 
 
@@ -838,7 +836,7 @@ def _stage_solver(self) -> str:
     def _solver(self, tmpdir, **options):
         """The solver of the current stage with the shared settings. CBC writes its scratch files to tmpdir."""
         if self._stage_solver() == 'highs':
-            # the cost stage reads its gap from the stats, not a CBC log, see _stage_gap
+            # HiGHS reports its bound in the stats, no log to parse, see _stage_gap
             options.pop('logPath', None)
             # HiGHS also stops at a relative gap of 1e-4 by default, CBC only on the absolute one
             options.setdefault('gapRel', 0)
@@ -847,7 +845,7 @@ def _solver(self, tmpdir, **options):
                 options.update(callbackTuple=(self._interrupt, None),
                                callbacksToActivate=[types.kCallbackMipInterrupt, types.kCallbackSimplexInterrupt,
                                                     types.kCallbackIpmInterrupt])
-            return BulkHiGHS(msg=False, threads=self.settings.num_threads, **options)
+            return pulp.HiGHS(msg=False, threads=self.settings.num_threads, **options)
 
         # pulp resolves the binary unless one is configured: a cbc on PATH, the 2.10 release build
         # the image installs. The cbcbox wheel of the pulp cbc extra is deliberately not used, its
@@ -883,20 +881,20 @@ def _hard_limit(self, seconds: float):
         finally:
             self._cutoff = None
 
-    def _stage_gap(self, stats, log_path, scale) -> float | None:
-        """Distance between the cost stage's schedule and the solver's bound, in currency."""
-        if self._stage_solver() != 'highs':
-            return self._cbc_gap(log_path, scale)
-        # pulp 4 leaves best_bound unset for HiGHS, the solver model has it. A proven solve has
-        # nothing left open, the way CBC's log prints no bound for one
+    @staticmethod
+    def _stage_gap(stats, scale) -> float | None:
+        """Distance between the cost stage's schedule and the solver's bound, in currency.
+
+        HiGHS reports its bound natively, CBC's is parsed from the log pulp is handed as logPath.
+        A proven solve has nothing left open, the way CBC's log prints no bound for one.
+        """
         if not stats.has_solution:
             return None
         if stats.status in PROVEN:
             return 0.
-        info = self.problem.solverModel.getInfo()
-        if not (isfinite(info.objective_function_value) and isfinite(info.mip_dual_bound)):
+        if stats.best_bound is None:
             return None
-        return abs(info.mip_dual_bound - info.objective_function_value) / scale
+        return abs(stats.best_bound - stats.objective) / scale
 
     @contextmanager
     def _timed(self, stage):
@@ -913,26 +911,6 @@ def _timed(self, stage):
             self.stage_seconds[stage] = round(
                 self.stage_seconds.get(stage, 0.) + time.monotonic() - started, 4)
 
-    @staticmethod
-    def _cbc_gap(log_path, scale) -> float | None:
-        """Distance between CBC's schedule and its bound, from the log, in currency.
-
-        The solution file carries no bound, only the log does: a proven solve prints none and the
-        gap is zero, a stopped one prints the bound it reached.
-        """
-        # no log, no solver run behind the values, see the faked stages in the tests
-        if not os.path.exists(log_path):
-            return None
-        with open(log_path) as log:
-            text = log.read()
-        value = re.search(r'^Objective value:\s+(\S+)', text, re.MULTILINE)
-        if value is None:
-            return None
-        bound = re.search(r'^(?:Lower|Upper) bound:\s+(\S+)', text, re.MULTILINE)
-        if bound is None:
-            return 0.
-        return abs(float(bound.group(1)) - float(value.group(1))) / scale
-
     def _pin_integers(self):
         """Freeze every integer variable on the value it currently holds, undo data returned.
 
@@ -1271,10 +1249,10 @@ def _probe_then_split(self, tmpdir, deadline) -> None:
         reserve = (0. if self.settings.time_limit is None
                    else self.settings.time_limit * PREFERENCE_TIME_SHARE)
         remaining = None if deadline is None else max(min(COST_TIME_LIMIT, deadline - reserve - time.monotonic()), 0.1)
-        log_path = os.path.join(tmpdir, 'cost.log')
         with self._timed('cost'):
-            self.stats = self.problem.solve(self._solver(tmpdir, timeLimit=remaining, logPath=log_path, gapAbs=scaled_gap))
-        self.cost_stage_gap = self._stage_gap(self.stats, log_path, scale)
+            self.stats = self.problem.solve(self._solver(tmpdir, timeLimit=remaining, gapAbs=scaled_gap,
+                                                         logPath=os.path.join(tmpdir, 'cost.log')))
+        self.cost_stage_gap = self._stage_gap(self.stats, scale)
 
         if self.stats.has_solution and self._is_feasible():
             # the cost stage is allowed to stop on its gap, so proven is not required of it, only

tests/test_highs.py

@@ -34,15 +34,29 @@ def test_stored_optimum_without_gap(case):
     assert result['objective_value'] == pytest.approx(expected['objective_value'], abs=1e-4, rel=1e-5)
 
 
-def test_cost_stage_gap_comes_from_the_solver_model():
-    # pulp 4 leaves best_bound unset for HiGHS, so the gap is read from HiGHS itself
+def test_cost_stage_gap_comes_from_the_stats():
     optimizer = highs('026-attenuate-grid-peaks', probe_seconds=0)
     optimizer.solve()
 
     assert optimizer.solve_path == 'split'
     assert optimizer.cost_stage_gap == 0  # proven, nothing left open
 
 
+def test_stopped_cost_stage_reports_its_gap():
+    # the bound comes back in the problem's sense: a maximised cost objective with a bound above it
+    optimizer = highs('028-attenuate-grid-peaks-long-horizon')
+    optimizer.create_model()
+    scale = optimizer.objective_scale
+    optimizer.problem.setObjective(optimizer.cost_objective * scale)
+
+    with TemporaryDirectory() as tmpdir, optimizer._timed('cost'):
+        stats = optimizer.problem.solve(optimizer._solver(tmpdir, timeLimit=0.5))
+
+    assert stats.status == pulp.LpSolveStatus.TimeLimit
+    assert stats.best_bound > stats.objective
+    assert optimizer._stage_gap(stats, scale) == (stats.best_bound - stats.objective) / scale
+
+
 def test_overrunning_solve_is_interrupted_and_keeps_its_incumbent():
     # a long horizon HiGHS needs about 3 s for, cut off after a tenth of a second
     optimizer = highs('028-attenuate-grid-peaks-long-horizon')

tests/test_hybrid.py

@@ -7,7 +7,6 @@
 import pytest
 from test_objective_split import build
 
-from optimizer.highs import BulkHiGHS
 from optimizer.optimizer import Optimizer
 from optimizer.settings import OptimizerSettings
 
@@ -54,9 +53,9 @@ def test_split_runs_on_highs(solvers):
 
     assert optimizer.solve_path == 'split'
     stages = dict(solvers)
-    assert stages['cost'] is BulkHiGHS
-    assert stages['tie_break'] is BulkHiGHS
-    # the gap of a HiGHS cost stage comes from the solver model, not a CBC log
+    assert stages['cost'] is pulp.HiGHS
+    assert stages['tie_break'] is pulp.HiGHS
+    # the gap of a HiGHS cost stage comes from the stats, not a CBC log
     assert optimizer.cost_stage_gap == 0