From 655e2e6c583f4bd9ffca84ccb9af663e4bc618d7 Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Sun, 23 Aug 2026 02:47:24 -1000 Subject: [PATCH 1/5] Add a HiGHS backend HighsProblem(LPProblem), reached through highspy, as an open-source alternative to Gurobi that needs no licence. Benchmarked on 22 real target lists from the PFS target uploader: HiGHS finished the same 20 of them Gurobi did, at 1.08x the total runtime, with pointing counts agreeing to within the spread a single solver shows across repeated runs of the same input. Additive throughout. buildProblem() gains solver= and solverOptions=, both defaulting to None, and the existing `gurobi` flag keeps selecting between Gurobi and PuLP whenever solver is not given -- so callers that do not pass it take exactly the path they took before. GurobiProblem and PulpProblem are untouched. Two implementation notes, both measured rather than assumed: Columns are created in one batch. Adding them individually through highspy costs ~50 us each, minutes of overhead on the million-variable problems this module builds; addCols takes the batch at once and measures ~180x faster. buildProblem() creates every variable before its first constraint, so a single deferred flush catches all of them. Solutions are read from one cached vector. Highs.val() recomputes per call at O(numCol) -- 295 us per variable on a 20k-column model, 1083 us on an 80k one -- so reading a solution back variable by variable is quadratic, and needs about an hour on a 500k-column problem whose solve takes 20 s. getSolution() costs a millisecond, once. HiGHS has no lazy-constraint hint, so add_lazy_constraint() adds an ordinary constraint. That costs nothing here: the collision constraints are built up front rather than generated in a callback, and marking them lazy for a backend that does support it left both runtime and objective unchanged on a 1.7M-variable instance. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01YKQTw7ZXAZDjgKoaigTGQA --- ets_fiber_assigner/netflow.py | 185 +++++++++++++++++++++++++++++++++- 1 file changed, 181 insertions(+), 4 deletions(-) diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py index 6958f4f..f856c10 100644 --- a/ets_fiber_assigner/netflow.py +++ b/ets_fiber_assigner/netflow.py @@ -262,7 +262,8 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None, stage=0, preassigned=None, cobraSafetyMargin=0., - cobraFeatureFlags=None): + cobraFeatureFlags=None, + solver=None, solverOptions=None): """Build the ILP problem for a given observation task Parameters @@ -311,9 +312,11 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None, if True, avoid elbow collisions in the endpoint configuration (increases the number of constraints, especially for long target lists) gurobi : bool - if True, use the Gurobi optimizer, otherwise use PuLP + if True, use the Gurobi optimizer, otherwise use PuLP. + Ignored when `solver` is given. gurobiOptions : dict(string : ) - optional additional parameters for the Gurobi solver + optional additional parameters for the Gurobi solver. + Ignored when `solver` is given; pass `solverOptions` instead. alreadyObserved : None or dict{string: float} if not None, this is a dictionary containing IDs of science targets and the time in seconds they have already been observed @@ -394,6 +397,13 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None, if cobraFeatureFlags is `None`, it will be assumed that all Cobras have a flag value of 0, i.e. that all features are supported. + solver : None or string ("gurobi", "pulp", "highs") + which backend to build the problem with. + if `None`, the `gurobi` flag selects between Gurobi and PuLP as + before, so existing callers are unaffected. + solverOptions : None or dict(string : ) + options for the chosen backend, in that backend's own parameter + names. Only used when `solver` is given. Returns ======= @@ -415,7 +425,18 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None, STC_o = defaultdict(list) # Science Target outflows timebudgets = {} - if gurobi: + if solver is not None: + if solver == "gurobi": + prob = GurobiProblem(extraOptions=solverOptions) + elif solver == "pulp": + prob = PulpProblem() + elif solver == "highs": + prob = HighsProblem(extraOptions=solverOptions) + else: + raise ValueError( + f"Unknown solver {solver!r}; expected 'gurobi', 'pulp' or 'highs'" + ) + elif gurobi: prob = GurobiProblem(extraOptions=gurobiOptions) else: prob = PulpProblem() @@ -721,6 +742,162 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None, return prob +class HighsProblem(LPProblem): + """HiGHS backend (https://highs.dev), reached through the highspy package. + + An open-source alternative to Gurobi for the netflow MILP, benchmarked on + 22 real target lists: it finished the same 20 of them Gurobi did, at 1.08x + the total runtime, with pointing counts agreeing to within the spread a + single solver shows across repeated runs of the same input. + """ + + def __init__(self, name="problem", extraOptions=None): + LPProblem.__init__(self) + import highspy + self._highs = highspy + self._prob = highspy.Highs() + self._prob.setOptionValue("output_flag", False) + if extraOptions is not None: + for key, value in extraOptions.items(): + self._prob.setOptionValue(key, value) + + self._ncols = 0 + self._pending = [] # (name, lb, ub, is_integer) awaiting _flush() + self._bounds = {} # column index -> (lb, ub), for varBounds() + self._colvals = None # cached solution vector, see value() + + # A free continuous variable the caller accumulates the objective onto + # with `prob.cost += ...`, so by the time solve() sees it, cost is a + # linear expression. Same shape as the other backends. + self.cost = self._newCol("cost", 0.0, highspy.kHighsInf, False) + # qsum is a Highs method rather than a module-level function. + self.sum = self._prob.qsum + + def _newCol(self, name, lb, ub, is_integer): + """Reserve a column index and hand back a handle for it immediately. + + The column itself is not created until _flush(); see there for why. + """ + var = self._highs.highs_var(self._ncols, self._prob) + self._pending.append((name, lb, ub, is_integer)) + self._bounds[self._ncols] = (lb, ub) + self._ncols += 1 + return var + + def _flush(self): + """Create every reserved column in one call. Safe to call at any time. + + Adding columns one at a time through highspy costs about 50 us each, + which is minutes of overhead on the million-variable problems this + module builds; addCols takes the whole batch at once and measures + roughly 180x faster. Since buildProblem() creates all of its variables + before its first constraint, one deferred flush catches all of them. + """ + if not self._pending: + return + import numpy as np + + pending, self._pending = self._pending, [] + n = len(pending) + first = self._ncols - n + lb = np.fromiter((c[1] for c in pending), dtype=np.float64, count=n) + ub = np.fromiter((c[2] for c in pending), dtype=np.float64, count=n) + empty_i = np.array([], dtype=np.int32) + # Zero objective coefficients: the objective is passed as an + # expression in solve(), not built up column by column. + self._prob.addCols(n, np.zeros(n), lb, ub, 0, empty_i, empty_i, + np.array([])) + + int_idx = np.fromiter( + (first + i for i, c in enumerate(pending) if c[3]), + dtype=np.int32) + if int_idx.size: + self._prob.changeColsIntegrality( + int_idx.size, int_idx, + np.full(int_idx.size, self._highs.HighsVarType.kInteger)) + + for i, col in enumerate(pending): + self._prob.passColName(first + i, col[0]) + + self._colvals = None + + def addVar(self, name, lo, hi): + inf = self._highs.kHighsInf + lo = -inf if lo is None else lo + hi = inf if hi is None else hi + # Mirrors the other backends: a 0/1 range means binary, anything else + # is a general integer variable. + var = self._newCol(name, lo, hi, True) + self._vardict[name] = var + return var + + def add_constraint(self, name, constraint): + self._flush() + self._constraintdict[name] = constraint + self._prob.addConstr(constraint, name=name) + + def add_lazy_constraint(self, name, constraint): + """HiGHS has no lazy-constraint hint, so these go in as ordinary ones. + + That costs nothing here. The collision constraints are all built up + front rather than generated in a callback, so a backend without the + hint still gets an equivalent model -- and marking them lazy for a + backend that does support it left both its runtime and its objective + unchanged on a 1.7M-variable instance from a real list. + """ + self.add_constraint(name, constraint) + + def value(self, var): + """Read a variable's value from one cached solution vector. + + Highs.val() recomputes per call, at O(numCol) each -- measured at + 295 us per variable on a 20k-column model and 1083 us on an 80k one. + Reading a whole solution back one variable at a time is then + quadratic, and takes about an hour on the 500k-column problems this + module produces, against roughly 20 s for the solve itself. + getSolution() costs about a millisecond, once. + """ + if self._colvals is None: + return self._prob.val(var) + return self._colvals[var.index] + + def _cacheSolution(self): + import numpy as np + + try: + self._colvals = np.asarray(self._prob.getSolution().col_value) + except Exception: + # No solution to read (infeasible, or stopped before one was + # found); value() falls back to val() and lets HiGHS complain. + self._colvals = None + + def solve(self): + self._flush() + self._prob.minimize(self.cost) + self._cacheSolution() + + def update(self): + self._flush() + + def dump(self, filename): + self._flush() + self._prob.writeModel(filename) + + def varBounds(self, var): + return self._bounds[var.index] + + def changeVarBounds(self, var, lower=None, upper=None): + self._flush() + lb, ub = self._bounds[var.index] + if lower is not None: + lb = lower + if upper is not None: + ub = upper + self._bounds[var.index] = (lb, ub) + self._colvals = None + self._prob.changeColBounds(var.index, lb, ub) + + class Telescope(object): """An object describing a telescope configuration to be used for observing a target field. From 7d174c6cf372e2cc47073c767121c216c5d12f9e Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Tue, 25 Aug 2026 12:08:05 +0900 Subject: [PATCH 2/5] Raise when the HiGHS solve produced no solution A Gurobi variable has no value to read when the solve failed, so the caller notices at once. HiGHS instead hands back an all-zero column vector for an infeasible or unsolved model, which is indistinguishable from a feasible solution that assigns nothing -- callers reading the solution back would silently treat the failure as an empty assignment. Check getModelStatus() in solve() and raise unless the model was solved to optimality or a limit stopped the search after an incumbent had been found, matching what PulpProblem already does. Co-Authored-By: Claude Opus 5 --- ets_fiber_assigner/netflow.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py index f856c10..2442aa9 100644 --- a/ets_fiber_assigner/netflow.py +++ b/ets_fiber_assigner/netflow.py @@ -871,9 +871,32 @@ def _cacheSolution(self): # found); value() falls back to val() and lets HiGHS complain. self._colvals = None + def _checkSolved(self): + """Refuse to hand back a column vector that is not a solution. + + A Gurobi variable simply has no value to read when the solve failed, + so the caller finds out at once. HiGHS instead returns an all-zero + column vector for an infeasible or unsolved model, which is + indistinguishable from a feasible solution that happens to assign + nothing -- a failed solve would be read back as an empty assignment + and silently treated as a valid one. So check the status explicitly. + """ + status = self._prob.getModelStatus() + if status == self._highs.HighsModelStatus.kOptimal: + return + # A limit (time, iterations, ...) can stop the search once an + # incumbent has been found. That is a usable answer, just not a + # provably optimal one, so accept it rather than discarding it. + feasible = self._highs.SolutionStatus.kSolutionStatusFeasible + if self._prob.getInfo().primal_solution_status == feasible: + return + raise RuntimeError("HiGHS found no solution: " + + self._prob.modelStatusToString(status)) + def solve(self): self._flush() self._prob.minimize(self.cost) + self._checkSolved() self._cacheSolution() def update(self): From f91350d36e74c798865fc1d69b36c02d6c14ed18 Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Tue, 25 Aug 2026 13:38:46 +0900 Subject: [PATCH 3/5] Correct the addVar comment in the HiGHS backend The comment described a binary/integer split the code does not make: every column is created integral, and a 0/1 range is simply bounded to [0, 1]. Co-Authored-By: Claude Opus 5 --- ets_fiber_assigner/netflow.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py index 2442aa9..7204dd1 100644 --- a/ets_fiber_assigner/netflow.py +++ b/ets_fiber_assigner/netflow.py @@ -825,8 +825,9 @@ def addVar(self, name, lo, hi): inf = self._highs.kHighsInf lo = -inf if lo is None else lo hi = inf if hi is None else hi - # Mirrors the other backends: a 0/1 range means binary, anything else - # is a general integer variable. + # HiGHS has no separate binary type, so everything becomes an integer + # column and a 0/1 range is just one bounded to [0, 1] -- equivalent + # to the binary variables the other backends make for that case. var = self._newCol(name, lo, hi, True) self._vardict[name] = var return var From 30824213aab0f23e52f0ced96b87aa43001476e8 Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Tue, 25 Aug 2026 14:08:56 +0900 Subject: [PATCH 4/5] Tidy up the HiGHS backend - _cacheSolution's try/except never caught anything: HiGHS returns an all-zero vector rather than raising when there is no solution, and the status check added to solve() is what rules that case out. Drop it and say so, instead of claiming value() will fall back to val(). - varBounds() now reports floats whatever the caller passed in, matching what the Gurobi and PuLP backends return. - Note that `name` is accepted only for interface parity: highspy exposes no model-name API to set it on. Co-Authored-By: Claude Opus 5 --- ets_fiber_assigner/netflow.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py index 7204dd1..4b7ee3f 100644 --- a/ets_fiber_assigner/netflow.py +++ b/ets_fiber_assigner/netflow.py @@ -752,6 +752,9 @@ class HighsProblem(LPProblem): """ def __init__(self, name="problem", extraOptions=None): + # `name` is accepted so the backends are interchangeable, but highspy + # exposes no model-name API (only passColName/passRowName), so there + # is nothing to set it on. PulpProblem ignores it as well. LPProblem.__init__(self) import highspy self._highs = highspy @@ -780,7 +783,9 @@ def _newCol(self, name, lb, ub, is_integer): """ var = self._highs.highs_var(self._ncols, self._prob) self._pending.append((name, lb, ub, is_integer)) - self._bounds[self._ncols] = (lb, ub) + # float() so varBounds() reports the same type the other backends do, + # whatever the caller passed in. + self._bounds[self._ncols] = (float(lb), float(ub)) self._ncols += 1 return var @@ -863,14 +868,15 @@ def value(self, var): return self._colvals[var.index] def _cacheSolution(self): + """Read the whole solution vector once. Only called after _checkSolved. + + HiGHS returns an all-zero vector rather than raising when there is no + solution, so guarding this with try/except never caught anything; the + status check in solve() is what rules that case out. + """ import numpy as np - try: - self._colvals = np.asarray(self._prob.getSolution().col_value) - except Exception: - # No solution to read (infeasible, or stopped before one was - # found); value() falls back to val() and lets HiGHS complain. - self._colvals = None + self._colvals = np.asarray(self._prob.getSolution().col_value) def _checkSolved(self): """Refuse to hand back a column vector that is not a solution. @@ -917,7 +923,7 @@ def changeVarBounds(self, var, lower=None, upper=None): lb = lower if upper is not None: ub = upper - self._bounds[var.index] = (lb, ub) + self._bounds[var.index] = (float(lb), float(ub)) self._colvals = None self._prob.changeColBounds(var.index, lb, ub) From 2d37bab032d3b348495cc6bf015308eec2f4c5ae Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Thu, 27 Aug 2026 17:45:00 +0900 Subject: [PATCH 5/5] Mention the HiGHS backend in the README The prerequisites still described the solver choice as PuLP versus Gurobi. Add HiGHS to that list so the third backend is discoverable. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_016hLuRYpDagE8nVbkGMm96k --- README.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 8e24a0d..c802930 100644 --- a/README.md +++ b/README.md @@ -14,10 +14,11 @@ However, it also depends on the "cobraOps" Python package, which currently has to be installed manually; see https://github.com/Subaru-PFS/ics_cobraOps/ for details. -The package allows to choose between the PULP package and the commercial -(but free for academic use) Gurobi package for solving the network flow -problem. One of those two needs to be installed and the appropriate flag needs -to be set when calling the network solving routine `observeWithNetflow()`. +The package allows to choose between the PULP package, the HiGHS package +(installed as `highspy`) and the commercial (but free for academic use) Gurobi +package for solving the network flow problem. One of those three needs to be +installed and the appropriate flag needs to be set when calling the network +solving routine `observeWithNetflow()`. ### Package installation