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
diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py
index 7f7e317..51ab193 100644
--- a/ets_fiber_assigner/netflow.py
+++ b/ets_fiber_assigner/netflow.py
@@ -1,3 +1,5 @@
+import array
+import logging
import numpy as np
from collections import defaultdict
from astropy.table import Table
@@ -267,7 +269,8 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None,
cobraFeatureFlags=None,
brokenCobrasMargin=0.,
targetCostOffset=None,
- avoidFiducials=True):
+ avoidFiducials=True,
+ solver=None, solverOptions=None):
"""Build the ILP problem for a given observation task
Parameters
@@ -316,9 +319,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
@@ -399,6 +404,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.
brokenCobrasMargin: float
defines the radius around broken Cobras, in which potential
@@ -435,7 +447,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()
@@ -746,6 +769,368 @@ def buildProblem(bench, targets, tpos, classdict, tvisit, vis_cost=None,
return prob
+class _RowBuffer(object):
+ """Rows collected by HighsProblem.add_constraint(), waiting for one addRows.
+
+ Five flat typed arrays (C int / double) rather than a list of expressions:
+ 4-8 bytes per entry, extend() runs in C, and np.frombuffer turns them into
+ numpy arrays without a copy. The expressions themselves live in
+ HighsProblem._constraintdict only.
+ """
+
+ def __init__(self):
+ self.clear()
+
+ def clear(self):
+ self.idx = array.array("i") # column indices, all rows concatenated
+ self.val = array.array("d") # matching coefficients
+ self.len = array.array("i") # nonzeros per row
+ self.lo = array.array("d") # row lower bounds
+ self.hi = array.array("d") # row upper bounds
+
+ def __len__(self):
+ return len(self.len)
+
+ def add(self, idxs, vals, lo, hi):
+ self.idx.extend(idxs)
+ self.val.extend(vals)
+ self.len.append(len(idxs))
+ self.lo.append(lo)
+ self.hi.append(hi)
+
+ def to_csr(self, ncols):
+ """Return (lo, hi, nnz, starts, idx, val) as Highs.addRows wants them.
+
+ Within each row the entries are ordered by column index, which is
+ what highs_linear_expression.unique_elements() does before addConstr
+ hands a row to HiGHS. Rows arrive here free of repeated columns
+ (add_constraint merges those), so a stable sort on a (row, column)
+ key reproduces that order exactly, and the matrix HiGHS stores --
+ and the .lp/.mps it writes -- is identical to the one-row-at-a-time
+ path. One integer key per entry sorts about 20x faster than
+ np.lexsort on the same data. Temporary memory is about 40 bytes per
+ nonzero.
+ """
+ m = len(self)
+ lens = np.frombuffer(self.len, dtype=np.intc)
+ idx = np.frombuffer(self.idx, dtype=np.intc)
+ val = np.frombuffer(self.val, dtype=np.float64)
+ lo = np.frombuffer(self.lo, dtype=np.float64)
+ hi = np.frombuffer(self.hi, dtype=np.float64)
+ # CSR row starts; addRows takes num_rows of them, no trailing sentinel.
+ # frombuffer views are read-only, so build starts as a fresh array.
+ starts = np.zeros(m, dtype=np.int32)
+ np.cumsum(lens[:-1], dtype=np.int32, out=starts[1:])
+ row_id = np.repeat(np.arange(m, dtype=np.int64), lens)
+ order = np.argsort(row_id * (ncols + 1) + idx, kind="stable")
+ return (lo, hi, idx.size, starts,
+ idx[order].astype(np.int32, copy=False), val[order])
+
+
+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.
+
+ Columns, rows and bound changes are buffered on the Python side and
+ handed to HiGHS in bulk (addCols / addRows / changeColsBounds) the first
+ time something needs the complete model: solve(), update() and dump().
+ Column and row names are passed to HiGHS in dump() alone, because
+ HiGHS's own log refers to rows and columns by index whether or not they
+ are named -- to see names, read the file dump() writes.
+
+ varBounds() and changeVarBounds() are instance methods here (they read
+ self._bounds), while GurobiProblem and PulpProblem define them as
+ staticmethods. Calling them through an instance works for all three
+ backends; calling them through the class (HighsProblem.varBounds(var))
+ does not work for this one.
+
+ Options: output_flag is set to False first, then extraOptions is applied,
+ so a caller can turn HiGHS's log back on. highspy 1.15.1 defaults worth
+ knowing when comparing against Gurobi: mip_rel_gap 1e-4 (Gurobi's MIPGap
+ default is also 1e-4), mip_abs_gap 1e-6, threads 0 (automatic),
+ mip_detect_symmetry True, time_limit inf (when a limit stops the search,
+ _checkSolved accepts a feasible incumbent), presolve "choose",
+ random_seed 0.
+
+ Requires highspy >= 1.15 (module-level highs_var).
+ """
+
+ 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
+ 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)
+
+ # Columns: reserved in _newCol(), created in _flush().
+ self._ncols = 0
+ self._pending = [] # is_integer per column awaiting _flush()
+ self._col_names = [] # one per column, passed to HiGHS in dump()
+ self._bounds = {} # column index -> (lb, ub): the bounds HiGHS
+ # gets at _flush(), and what varBounds() reports
+ self._bounds_changed = set() # created columns whose _bounds moved
+ self._colvals = None # cached solution vector, see value()
+
+ # Rows: collected in add_constraint(), created in _flushRows().
+ self._rows = _RowBuffer()
+ self._row_names = [] # one per row, passed to HiGHS in dump()
+ self._nrow_flushes = 0 # addRows calls made, for the log
+
+ # 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(is_integer)
+ self._col_names.append(name)
+ # 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
+
+ def _flush(self):
+ """Create every reserved column and apply pending bound changes, in
+ bulk. 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 and not self._bounds_changed:
+ return
+ if self._pending:
+ pending, self._pending = self._pending, []
+ n = len(pending)
+ first = self._ncols - n
+ lb, ub = self._boundArrays(range(first, self._ncols))
+ 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, is_int in enumerate(pending) if is_int),
+ dtype=np.int32)
+ if int_idx.size:
+ self._prob.changeColsIntegrality(
+ int_idx.size, int_idx,
+ np.full(int_idx.size, self._highs.HighsVarType.kInteger))
+ if self._bounds_changed:
+ # changeVarBounds() on columns HiGHS already had; applied here,
+ # in one call, like Gurobi applies pending changes in update().
+ idx = np.fromiter(sorted(self._bounds_changed), dtype=np.int32)
+ lb, ub = self._boundArrays(idx)
+ self._prob.changeColsBounds(idx.size, idx, lb, ub)
+ self._bounds_changed.clear()
+ self._colvals = None
+
+ def _boundArrays(self, indices):
+ lb = np.fromiter((self._bounds[i][0] for i in indices), dtype=np.float64)
+ ub = np.fromiter((self._bounds[i][1] for i in indices), dtype=np.float64)
+ return lb, ub
+
+ def _flushRows(self):
+ """Create every collected row in one addRows call. Idempotent.
+
+ Highs.addConstr costs about 19 us per row (unique_elements 9.5,
+ addRow 4.8, passRowName 0.9, wrapper 3), linear in the row count but
+ still 10-17 s per million rows. Handing the same rows to addRows as
+ one CSR block takes about 2 us per row including the sort.
+ """
+ self._flush() # rows refer to columns, so those must exist first
+ m = len(self._rows)
+ if m == 0:
+ return
+ lo, hi, nnz, starts, idx, val = self._rows.to_csr(self._ncols)
+ status = self._prob.addRows(m, lo, hi, nnz, starts, idx, val)
+ if status != self._highs.HighsStatus.kOk:
+ raise RuntimeError("HiGHS addRows failed: " + str(status))
+ self._rows.clear() # _row_names stays: dump() still needs the names
+ self._nrow_flushes += 1
+ logging.getLogger(__name__).info(
+ "HighsProblem: row flush #%d, %d rows, %d nonzeros",
+ self._nrow_flushes, m, nnz)
+ self._colvals = None
+
+ def _flushAll(self):
+ """Bring the HiGHS model up to date with everything added so far."""
+ self._flush()
+ self._flushRows()
+
+ def _passNames(self):
+ """Hand HiGHS every column and row name. Only dump() needs them.
+
+ There is no bulk-naming API, so this is one call per name, about
+ 1 us each; writeModel itself is far slower than that, so simply
+ repeating it on every dump() is cheaper than tracking what HiGHS
+ has already been told. Call after _flushAll(): a name can only be
+ attached to a column or row that exists.
+ """
+ for i, name in enumerate(self._col_names):
+ self._prob.passColName(i, name)
+ for i, name in enumerate(self._row_names):
+ self._prob.passRowName(i, name)
+
+ 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
+ # 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
+
+ def add_constraint(self, name, constraint):
+ bounds = constraint.bounds
+ if bounds is None:
+ # Same condition, and the same moment, at which Highs.addConstr
+ # would have refused the expression.
+ raise ValueError(
+ "Constraint bounds must be set via comparison (>=, ==, <=)")
+ self._constraintdict[name] = constraint
+ idxs, vals = constraint.idxs, constraint.vals
+ if len(set(idxs)) != len(idxs):
+ # A column appearing more than once in one row (buildProblem does
+ # not produce these, but keep the general case right). Let
+ # highspy's own unique_elements() merge them so the coefficients
+ # come out bit-identical to what addConstr would have stored.
+ u_idx, u_val = constraint.unique_elements()
+ idxs, vals = u_idx.tolist(), u_val.tolist()
+ # constant is ignored on purpose: the comparison that set `bounds`
+ # already moved it to the right-hand side, exactly as addConstr does.
+ self._rows.add(idxs, vals, bounds[0], bounds[1])
+ self._row_names.append(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.
+
+ Raises RuntimeError when there is no solution to read: before
+ solve(), after a solve() that failed, for a variable added since the
+ last solve(), and once a change -- new columns or rows, or
+ changeVarBounds() -- has reached HiGHS through update(), dump() or
+ solve(). Until that flush the previous solution stays readable. This
+ is the rule Gurobi follows as well: while a modification is pending,
+ var.X can still be read; once update() has applied it, var.X raises.
+ HiGHS itself would hand back 0.0 or the previous solution in every
+ one of these cases, which a caller cannot tell apart from a genuine
+ assignment.
+ """
+ if self._colvals is None or var.index >= self._colvals.size:
+ raise RuntimeError("HighsProblem.value(): no valid solution for "
+ "this variable; call solve() (again) first")
+ 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.
+ """
+ self._colvals = np.asarray(self._prob.getSolution().col_value)
+
+ 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._colvals = None # a failed solve must not leave an old vector
+ self._flushAll()
+ self._prob.minimize(self.cost)
+ self._checkSolved()
+ self._cacheSolution()
+
+ def update(self):
+ self._flushAll()
+
+ def dump(self, filename):
+ self._flushAll()
+ self._passNames()
+ # The objective reaches HiGHS in solve() only, so a file written
+ # before solve() has an empty objective. GurobiProblem.dump() behaves
+ # the same way (its setObjective is in solve() as well).
+ self._prob.writeModel(filename)
+
+ def varBounds(self, var):
+ return self._bounds[var.index]
+
+ def changeVarBounds(self, var, lower=None, upper=None):
+ """Record new bounds; they reach HiGHS at the next flush.
+
+ A column that is itself still pending is simply created with the new
+ bounds. One that HiGHS already has is queued for changeColsBounds()
+ in _flush(). Either way nothing is sent now, the same way a Gurobi
+ bound change stays pending until update() or optimize().
+ """
+ lb, ub = self._bounds[var.index]
+ if lower is not None:
+ lb = lower
+ if upper is not None:
+ ub = upper
+ self._bounds[var.index] = (float(lb), float(ub))
+ if var.index < self._ncols - len(self._pending):
+ self._bounds_changed.add(var.index)
+
+
class Telescope(object):
"""An object describing a telescope configuration to be used for observing
a target field.
diff --git a/misc/bench_highs_problem.py b/misc/bench_highs_problem.py
new file mode 100644
index 0000000..7d116f5
--- /dev/null
+++ b/misc/bench_highs_problem.py
@@ -0,0 +1,152 @@
+"""Task 0 benchmark: drive HighsProblem directly with a synthetic flow-like model.
+
+Usage: python misc/bench_highs_problem.py N [--no-solve] [--time-limit S] [--chunk K] [--threads T]
+
+Columns (N total, approximately):
+ T targets, C = T/2 cobras. Each target has 1..5 arcs to cobras (avg 3)
+ plus one sink arc -> ~4T columns. A handful of (0, None) overflow columns
+ exercise the None-bound path.
+Rows (M ~ 0.6 N):
+ collision pairs x_a + x_b <= 1 (T rows, 2 nnz) added "lazy"
+ cobra capacity sum(arcs into c) <= 1 (C rows, ~6 nnz)
+ target conservation sum(arcs of t) + sink == 1 (T rows, ~4 nnz)
+Objective: prob.cost += var * coef for every column.
+"""
+import argparse
+import sys
+import time
+
+import numpy as np
+
+import os
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
+from ets_fiber_assigner.netflow import HighsProblem # noqa: E402
+
+
+def build_and_time(N, cls, seed=1, solve=True, time_limit=600.0, chunk=20000,
+ threads=None, keep_prob=False):
+ rng = np.random.default_rng(seed)
+ T = max(1, N // 4)
+ C = max(1, T // 2)
+ opts = {"time_limit": float(time_limit)}
+ if threads is not None:
+ opts["threads"] = int(threads)
+ prob = cls(extraOptions=opts)
+ res = {"N_target": N}
+
+ # ---- variables + objective accumulation ------------------------------
+ t0 = time.perf_counter()
+ arcs_by_cobra = [[] for _ in range(C)]
+ arcs_by_target = []
+ all_arcs = []
+ for t in range(T):
+ k = int(rng.integers(1, 6))
+ cobras = rng.choice(C, size=min(k, C), replace=False)
+ tarcs = []
+ for c in cobras:
+ f = prob.addVar(f"Tv_Cv_{t}_{c}", 0, 1)
+ prob.cost += f * float(rng.random())
+ arcs_by_cobra[c].append(f)
+ tarcs.append(f)
+ all_arcs.append(f)
+ s = prob.addVar(f"ST_sink_{t}", 0, 1)
+ prob.cost += s * 10.0
+ tarcs.append(s)
+ arcs_by_target.append(tarcs)
+ for j in range(8):
+ f = prob.addVar(f"STC_sink_{j}", 0, None)
+ prob.cost += f * 1.0
+ t1 = time.perf_counter()
+ res["t_addVar_loop"] = t1 - t0 # Python-side, no HiGHS call yet
+ prob.update() # forces _flush (addCols [+names])
+ t2 = time.perf_counter()
+ res["t_flush"] = t2 - t1
+ res["t_addVar_total"] = t2 - t0
+ res["ncols"] = prob._prob.getNumCol()
+
+ # ---- constraints -----------------------------------------------------
+ # Build all expressions first so expression-building cost is separated
+ # from the add_constraint cost.
+ t3 = time.perf_counter()
+ exprs = []
+ narcs = len(all_arcs)
+ for p in range(T):
+ i, j = rng.integers(0, narcs, size=2)
+ if i == j:
+ continue
+ exprs.append((f"Coll_{p}", prob.sum([all_arcs[i], all_arcs[j]]) <= 1))
+ ncoll = len(exprs)
+ for c in range(C):
+ if arcs_by_cobra[c]:
+ exprs.append((f"Cvlim_{c}", prob.sum(arcs_by_cobra[c]) <= 1))
+ for t in range(T):
+ exprs.append((f"TvIO_{t}", prob.sum(arcs_by_target[t]) == 1))
+ t4 = time.perf_counter()
+ res["t_expr_build"] = t4 - t3
+ res["nrows_planned"] = len(exprs)
+
+ chunk_times = []
+ tc = time.perf_counter()
+ for k, (name, e) in enumerate(exprs):
+ if k < ncoll:
+ prob.add_lazy_constraint(name, e)
+ else:
+ prob.add_constraint(name, e)
+ if (k + 1) % chunk == 0:
+ now = time.perf_counter()
+ chunk_times.append(now - tc)
+ tc = now
+ t5 = time.perf_counter()
+ res["t_add_constraint"] = t5 - t4
+ res["chunk_times"] = chunk_times
+ # HighsProblem may buffer rows until something needs the full model;
+ # update() forces that, so the row flush is timed on its own here.
+ prob.update()
+ t5b = time.perf_counter()
+ res["t_flush_rows"] = t5b - t5
+ res["t_constraints_total"] = t5b - t4
+ res["nrows"] = prob._prob.getNumRow()
+ res["nnz"] = prob._prob.getNumNz()
+ res["nrow_flushes"] = getattr(prob, "_nrow_flushes", None)
+
+ if solve:
+ t6 = time.perf_counter()
+ prob.solve()
+ t7 = time.perf_counter()
+ res["t_solve"] = t7 - t6
+ res["status"] = prob._prob.modelStatusToString(prob._prob.getModelStatus())
+ res["objective"] = prob._prob.getObjectiveValue()
+ t8 = time.perf_counter()
+ vals = [prob.value(v) for v in prob._vardict.values()]
+ t9 = time.perf_counter()
+ res["t_value_all"] = t9 - t8
+ res["sum_values"] = float(sum(vals))
+ if keep_prob:
+ res["prob"] = prob
+ return res
+
+
+def main():
+ ap = argparse.ArgumentParser()
+ ap.add_argument("N", type=float)
+ ap.add_argument("--no-solve", action="store_true")
+ ap.add_argument("--time-limit", type=float, default=600.0)
+ ap.add_argument("--chunk", type=int, default=20000)
+ ap.add_argument("--threads", type=int, default=None)
+ a = ap.parse_args()
+ cls = HighsProblem
+ res = build_and_time(int(a.N), cls, solve=not a.no_solve,
+ time_limit=a.time_limit, chunk=a.chunk,
+ threads=a.threads)
+ print(f"=== {cls.__name__} N={int(a.N):,} ===")
+ for k, v in res.items():
+ if k == "chunk_times":
+ print(f" {k:18s} " + " ".join(f"{x:.2f}" for x in v))
+ elif isinstance(v, float):
+ print(f" {k:18s} {v:.4f}")
+ else:
+ print(f" {k:18s} {v}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/misc/compare_highs_real_instance.py b/misc/compare_highs_real_instance.py
new file mode 100644
index 0000000..288928f
--- /dev/null
+++ b/misc/compare_highs_real_instance.py
@@ -0,0 +1,331 @@
+"""Compare the HiGHS and Gurobi backends of buildProblem() on a real instance.
+
+Builds one buildProblem() instance from a target list and a few pointings
+the way pfs_target_uploader's PPP does (same Bench, classdict, cobraMoveCost
+and solver options), once with solver="gurobi" and once with solver="highs",
+and compares the two models by content: column bounds and integrality, rows
+as (lower, upper, coefficients keyed by column name) after normalising the
+sign convention, and the objective as accumulated on prob.cost. With --solve
+it also compares objective value, status, solve time and the set of Tv_Cv_*
+arcs that ended up at 1.
+
+buildProblem() draws from numpy's global random state (RandomTargetSelector
+in _get_vis_and_elbow), so the order in which it creates variables changes
+from one call to the next unless the state is reset; np.random.seed(--seed)
+is called before each build. The comparison itself is by name and does not
+depend on that order.
+
+Usage:
+ python misc/compare_highs_real_instance.py INPUT_DIR [--nvisit 4] [--solve]
+ [--gap 0.0] [--time-limit 600] [--out DIR] [--max-targets N] [--seed 20]
+INPUT_DIR holds target_.ecsv (ob_code, ra, dec, exptime, priority) and
+ppc_.ecsv (ppc_ra, ppc_dec, ppc_pa). Needs gurobipy, highspy, cobraOps,
+cobraCharmer and pfs.instdata importable; the spt_target_uploader venv has
+them all. gurobipy's pip wheel comes with a size-limited licence (2000
+variables / constraints), enough for --max-targets 300 --nvisit 1; larger
+instances need a full Gurobi licence.
+
+With --gap 0 both solvers return a proven optimum and the assignments should
+coincide whenever the optimum is unique; with a positive gap, only the
+objective values are expected to agree, to within twice the gap.
+"""
+import argparse
+import glob
+import logging
+import os
+import sys
+import tempfile
+import time
+
+import numpy as np
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+REPO_ROOT = os.path.dirname(HERE)
+sys.path.insert(0, REPO_ROOT)
+
+import ets_fiber_assigner.netflow as nf # noqa: E402
+
+BIG = 1e30 # anything beyond this is treated as infinite (Gurobi uses 1e100)
+
+
+def make_bench():
+ from pfs.instdata import setup_envvar
+ from ics.cobraOps.Bench import Bench
+ from ics.cobraCharmer.cobraCoach.cobraCoach import CobraCoach
+ setup_envvar()
+ with tempfile.TemporaryDirectory() as d:
+ cc = CobraCoach(loadModel=True, trajectoryMode=True, rootDir=d)
+ return Bench(cobraCoach=cc, blackDotsMargin=1.65)
+
+
+def classdict_like_ppp():
+ # NetflowPreparation() in pfs_target_uploader/utils/ppp.py
+ return {f"sci_P{p}": {"nonObservationCost": 100 - 10 * p,
+ "partialObservationCost": 200, "calib": False}
+ for p in range(10)}
+
+
+def observation_time(ra, dec):
+ try:
+ from pfs_target_uploader.utils.ppp import set_observation_time
+ return set_observation_time(ra, dec=dec)
+ except Exception as e: # noqa: BLE001
+ print(f"set_observation_time unavailable ({e}); using a fixed time")
+ return "2026-10-01T10:00:00Z"
+
+
+def load_inputs(input_dir, nvisit, max_targets):
+ from astropy.table import Table
+ tfile = glob.glob(os.path.join(input_dir, "target_*.ecsv"))[0]
+ pfile = glob.glob(os.path.join(input_dir, "ppc_*.ecsv"))[0]
+ tab = Table.read(tfile)
+ if max_targets and len(tab) > max_targets:
+ tab = tab[:max_targets]
+ ppc = Table.read(pfile)
+ seen, tel = set(), []
+ for row in ppc:
+ key = (float(row["ppc_ra"]), float(row["ppc_dec"]), float(row["ppc_pa"]))
+ if key in seen:
+ continue
+ seen.add(key)
+ tel.append(key)
+ if len(tel) == nvisit:
+ break
+ tgt = [nf.ScienceTarget(r["ob_code"], r["ra"], r["dec"], r["exptime"],
+ r["priority"], "sci") for r in tab]
+ return tgt, tel
+
+
+# ---------------------------------------------------------------------------
+# canonical model: backend-independent description by column name
+# ---------------------------------------------------------------------------
+def _inf(x):
+ x = float(x)
+ if x >= BIG:
+ return float("inf")
+ if x <= -BIG:
+ return float("-inf")
+ return x
+
+
+def _canonical_row(lo, hi, terms):
+ """terms: dict name -> coef. Fix the sign so the alphabetically first
+ column has a positive coefficient; flip and swap the bounds if not."""
+ terms = dict(terms)
+ if terms:
+ first = min(terms)
+ if terms[first] < 0:
+ terms = {k: -v for k, v in terms.items()}
+ lo, hi = -hi, -lo
+ return (round(_inf(lo), 9), round(_inf(hi), 9),
+ tuple(sorted((k, round(v, 9)) for k, v in terms.items())))
+
+
+def canonical_highs(prob):
+ import highspy
+ prob.update()
+ prob._passNames()
+ lp = prob._prob.getLp()
+ cnames = list(lp.col_names_)
+ integ = list(lp.integrality_) if len(lp.integrality_) else [None] * lp.num_col_
+ cols = {n: (_inf(lo), _inf(hi), it == highspy.HighsVarType.kInteger)
+ for n, lo, hi, it in zip(cnames, lp.col_lower_, lp.col_upper_, integ)}
+ start = np.asarray(lp.a_matrix_.start_)
+ index = np.asarray(lp.a_matrix_.index_)
+ value = np.asarray(lp.a_matrix_.value_)
+ entries = {r: {} for r in range(lp.num_row_)}
+ if lp.a_matrix_.format_ == highspy.MatrixFormat.kRowwise:
+ for r in range(lp.num_row_):
+ for k in range(start[r], start[r + 1]):
+ entries[r][cnames[index[k]]] = float(value[k])
+ else:
+ for c in range(lp.num_col_):
+ for k in range(start[c], start[c + 1]):
+ entries[index[k]][cnames[c]] = float(value[k])
+ rows = sorted(_canonical_row(lo, hi, entries[r])
+ for r, (lo, hi) in enumerate(zip(lp.row_lower_, lp.row_upper_)))
+ obj = {}
+ cost = prob.cost
+ if hasattr(cost, "idxs"):
+ for i, v in zip(cost.idxs, cost.vals):
+ obj[cnames[i]] = obj.get(cnames[i], 0.0) + float(v)
+ const = cost.constant or 0.0
+ else:
+ obj[cnames[cost.index]] = 1.0
+ const = 0.0
+ obj = {k: round(v, 9) for k, v in obj.items()}
+ return cols, rows, obj, float(const)
+
+
+def canonical_gurobi(prob):
+ import gurobipy as gbp
+ m = prob._prob
+ m.update()
+ cols = {v.VarName: (_inf(v.LB), _inf(v.UB), v.VType in (gbp.GRB.BINARY, gbp.GRB.INTEGER))
+ for v in m.getVars()}
+ rows = []
+ for c in m.getConstrs():
+ row = m.getRow(c)
+ terms = {}
+ for i in range(row.size()):
+ n = row.getVar(i).VarName
+ terms[n] = terms.get(n, 0.0) + float(row.getCoeff(i))
+ rhs = float(c.RHS)
+ if c.Sense == gbp.GRB.LESS_EQUAL:
+ lo, hi = float("-inf"), rhs
+ elif c.Sense == gbp.GRB.GREATER_EQUAL:
+ lo, hi = rhs, float("inf")
+ else:
+ lo, hi = rhs, rhs
+ rows.append(_canonical_row(lo, hi, terms))
+ rows.sort()
+ obj = {}
+ cost = prob.cost
+ if isinstance(cost, gbp.Var):
+ obj[cost.VarName] = 1.0
+ const = 0.0
+ else:
+ for i in range(cost.size()):
+ n = cost.getVar(i).VarName
+ obj[n] = obj.get(n, 0.0) + float(cost.getCoeff(i))
+ const = float(cost.getConstant())
+ obj = {k: round(v, 9) for k, v in obj.items()}
+ return cols, rows, obj, const
+
+
+def compare_canonical(a, b):
+ """Return a list of human-readable differences (empty when identical)."""
+ diffs = []
+ for what, da, db in (("columns", a[0], b[0]), ("objective", a[2], b[2])):
+ if set(da) != set(db):
+ diffs.append(f"{what}: name sets differ ({len(set(da) ^ set(db))} names)")
+ continue
+ bad = [k for k in da if da[k] != db[k]]
+ if bad:
+ diffs.append(f"{what}: {len(bad)} entries differ, e.g. {bad[0]}: {da[bad[0]]} vs {db[bad[0]]}")
+ if a[1] != b[1]:
+ sa, sb = set(a[1]), set(b[1])
+ ex = next(iter(sa ^ sb), None)
+ diffs.append(f"rows: {len(a[1])} vs {len(b[1])} rows, {len(sa ^ sb)} differ, e.g. {ex}")
+ if a[3] != b[3]:
+ diffs.append(f"objective constant {a[3]} vs {b[3]}")
+ return diffs
+
+
+# ---------------------------------------------------------------------------
+def build(solver, bench, tgt, tpos, classdict, nvisit, gap, time_limit):
+ if solver == "gurobi":
+ opts = {"MIPGap": gap, "Seed": 0, "OutputFlag": 0, "TimeLimit": float(time_limit)}
+ else:
+ opts = {"mip_rel_gap": gap, "random_seed": 0, "output_flag": False,
+ "time_limit": float(time_limit)}
+ t0 = time.perf_counter()
+ prob = nf.buildProblem(
+ bench, tgt, tpos, classdict, 900.0, [0] * nvisit,
+ cobraMoveCost=lambda d: 0.1 * d,
+ collision_distance=2.0, elbow_collisions=True,
+ solver=solver, solverOptions=opts,
+ alreadyObserved={}, forbiddenPairs=[[] for _ in range(nvisit)],
+ avoidFiducials=False, brokenCobrasMargin=0.0)
+ return prob, time.perf_counter() - t0
+
+
+def sizes(prob):
+ if isinstance(prob, nf.GurobiProblem):
+ m = prob._prob
+ m.update()
+ return m.NumVars, m.NumConstrs, m.NumNZs
+ prob.update()
+ h = prob._prob
+ return h.getNumCol(), h.getNumRow(), h.getNumNz()
+
+
+def status_and_objective(prob):
+ if isinstance(prob, nf.GurobiProblem):
+ import gurobipy as gbp
+ names = {getattr(gbp.GRB, n): n for n in ("OPTIMAL", "TIME_LIMIT", "INFEASIBLE",
+ "UNBOUNDED", "INTERRUPTED", "SUBOPTIMAL")}
+ return names.get(prob._prob.Status, str(prob._prob.Status)), prob._prob.ObjVal
+ h = prob._prob
+ return h.modelStatusToString(h.getModelStatus()), h.getObjectiveValue()
+
+
+def main():
+ ap = argparse.ArgumentParser()
+ ap.add_argument("input_dir")
+ ap.add_argument("--nvisit", type=int, default=4)
+ ap.add_argument("--solve", action="store_true")
+ ap.add_argument("--gap", type=float, default=0.0,
+ help="MIP relative gap for both solvers (default 0: proven optimum)")
+ ap.add_argument("--time-limit", type=float, default=600.0)
+ ap.add_argument("--out", default=None)
+ ap.add_argument("--max-targets", type=int, default=0)
+ ap.add_argument("--seed", type=int, default=20)
+ a = ap.parse_args()
+
+ logging.basicConfig(level=logging.INFO, format="%(levelname)s %(name)s: %(message)s")
+ for noisy in ("cobraCoach", "butler", "root"):
+ logging.getLogger(noisy).setLevel(logging.WARNING)
+
+ out = a.out or tempfile.mkdtemp(prefix="highs_vs_gurobi_")
+ os.makedirs(out, exist_ok=True)
+
+ tgt, tel = load_inputs(a.input_dir, a.nvisit, a.max_targets)
+ otime = observation_time(tel[0][0], tel[0][1])
+ telescopes = [nf.Telescope(ra, dec, pa, otime) for ra, dec, pa in tel]
+ print(f"{len(tgt)} targets, {len(telescopes)} pointings, otime {otime}")
+ bench = make_bench()
+ tpos = [t.get_fp_positions(tgt) for t in telescopes]
+ classdict = classdict_like_ppp()
+
+ results = {}
+ for solver in ("gurobi", "highs"):
+ np.random.seed(a.seed) # buildProblem's RandomTargetSelector uses np.random
+ prob, t_build = build(solver, bench, tgt, tpos, classdict, len(telescopes),
+ a.gap, a.time_limit)
+ r = {"t_build": t_build}
+ r["ncols"], r["nrows"], r["nnz"] = sizes(prob)
+ r["canon"] = canonical_gurobi(prob) if solver == "gurobi" else canonical_highs(prob)
+ prob.dump(os.path.join(out, f"{solver}.lp"))
+ if a.solve:
+ t0 = time.perf_counter()
+ prob.solve()
+ r["t_solve"] = time.perf_counter() - t0
+ r["status"], r["objective"] = status_and_objective(prob)
+ r["assigned"] = {k for k, v in prob._vardict.items()
+ if k.startswith("Tv_Cv_") and prob.value(v) > 0.5}
+ results[solver] = r
+ print(f"[{solver:6s}] build {t_build:.2f}s cols {r['ncols']} rows {r['nrows']} nnz {r['nnz']}"
+ + (f" solve {r['t_solve']:.1f}s {r['status']} obj {r['objective']:.6f} "
+ f"assigned {len(r['assigned'])}" if a.solve else ""))
+ del prob
+
+ G, H = results["gurobi"], results["highs"]
+ ok = (G["ncols"], G["nrows"], G["nnz"]) == (H["ncols"], H["nrows"], H["nnz"])
+ print(f"sizes identical: {ok}")
+ diffs = compare_canonical(G["canon"], H["canon"])
+ print(f"model identical by content (columns, rows, objective): {not diffs}")
+ for d in diffs:
+ print(" ", d)
+ ok &= not diffs
+ if a.solve:
+ og, oh = G["objective"], H["objective"]
+ tol = (2.0 * a.gap + 1e-9) * max(1.0, abs(og))
+ same_obj = abs(og - oh) <= tol
+ same_asg = G["assigned"] == H["assigned"]
+ print(f"objective agree within tolerance {tol:.3g}: {same_obj} (|diff| {abs(og - oh):.3g}) "
+ f"status gurobi {G['status']} / highs {H['status']} "
+ f"assignment identical: {same_asg} (sym. diff {len(G['assigned'] ^ H['assigned'])})")
+ ok &= same_obj
+ if a.gap == 0.0 and not same_asg:
+ print(" note: gap 0 but assignments differ -> the optimum is degenerate "
+ "(check objective agreement above)")
+ print(f"build time gurobi {G['t_build']:.2f}s / highs {H['t_build']:.2f}s"
+ + (f"; solve gurobi {G['t_solve']:.1f}s / highs {H['t_solve']:.1f}s" if a.solve else "")
+ + f"; files in {out}")
+ print("RESULT:", "PASS" if ok else "FAIL")
+ return 0 if ok else 1
+
+
+if __name__ == "__main__":
+ sys.exit(main())
diff --git a/tests/test_highs_problem.py b/tests/test_highs_problem.py
new file mode 100644
index 0000000..b5cd185
--- /dev/null
+++ b/tests/test_highs_problem.py
@@ -0,0 +1,589 @@
+"""Tests for HighsProblem, the HiGHS backend in ets_fiber_assigner/netflow.py.
+
+Three kinds of checks:
+
+1. HighsProblem's own behaviour: when value() may be read, what dump()
+ writes before and after solve(), names, row flushes, error cases.
+2. Model identity: HighsProblem buffers rows and hands them to HiGHS in one
+ addRows call. The model that produces must be the one highspy's own
+ one-row-at-a-time API (Highs.addVariable / Highs.addConstr) builds from
+ the same expressions -- matrix arrays and the written .lp/.mps files are
+ compared byte for byte.
+3. Solutions: small models with a known optimum, and agreement with
+ GurobiProblem on the same models (skipped when gurobipy is not
+ importable; the pip wheel's size-limited licence covers these models).
+
+Run with: python -m pytest tests -v
+"""
+import filecmp
+import os
+import sys
+
+import numpy as np
+import pytest
+
+highspy = pytest.importorskip("highspy")
+
+REPO_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.insert(0, REPO_ROOT)
+
+from ets_fiber_assigner.netflow import HighsProblem, GurobiProblem # noqa: E402
+
+try:
+ import gurobipy # noqa: F401
+ HAVE_GUROBI = True
+except ImportError:
+ HAVE_GUROBI = False
+
+needs_gurobi = pytest.mark.skipif(not HAVE_GUROBI, reason="gurobipy not installed")
+
+
+# ---------------------------------------------------------------------------
+# Reference: the same LPProblem-style interface on top of highspy's own
+# convenience API, one column and one row at a time. This is the path
+# HighsProblem replaced with bulk addCols/addRows, so the two must build
+# identical models.
+# ---------------------------------------------------------------------------
+class ReferenceHighs(object):
+ def __init__(self):
+ self._prob = highspy.Highs()
+ self._prob.setOptionValue("output_flag", False)
+ self.cost = self._prob.addVariable(0.0, highspy.kHighsInf, name="cost")
+ self.sum = self._prob.qsum
+ self._vardict = {}
+ self._constraintdict = {}
+ self._bounds = {}
+
+ def addVar(self, name, lo, hi):
+ lo = -highspy.kHighsInf if lo is None else lo
+ hi = highspy.kHighsInf if hi is None else hi
+ var = self._prob.addIntegral(lo, hi, name=name)
+ self._vardict[name] = var
+ self._bounds[var.index] = (float(lo), float(hi))
+ return var
+
+ def add_constraint(self, name, constraint):
+ self._constraintdict[name] = constraint
+ self._prob.addConstr(constraint, name=name)
+
+ add_lazy_constraint = add_constraint
+
+ def update(self):
+ pass
+
+ def dump(self, filename):
+ self._prob.writeModel(filename)
+
+ def solve(self):
+ self._prob.minimize(self.cost)
+
+ def value(self, var):
+ return self._prob.val(var)
+
+ def varBounds(self, var):
+ return self._bounds[var.index]
+
+ def changeVarBounds(self, var, lower=None, upper=None):
+ lb, ub = self._bounds[var.index]
+ lb = lb if lower is None else lower
+ ub = ub if upper is None else upper
+ self._bounds[var.index] = (float(lb), float(ub))
+ self._prob.changeColBounds(var.index, lb, ub)
+
+
+def make_gurobi():
+ return GurobiProblem(extraOptions={"OutputFlag": 0, "MIPGap": 0.0})
+
+
+def make_highs():
+ return HighsProblem(extraOptions={"mip_rel_gap": 0.0})
+
+
+# ---------------------------------------------------------------------------
+# Build sequences, written against the LPProblem interface so the same code
+# drives HighsProblem, ReferenceHighs and GurobiProblem. Each returns the
+# handles it created.
+# ---------------------------------------------------------------------------
+def scen_mixed(prob):
+ x = prob.addVar("x", 0, 1)
+ y = prob.addVar("y", 0, 1)
+ z = prob.addVar("z", 0, None)
+ w = prob.addVar("w", None, None)
+ u = prob.addVar("u", 2, 5)
+ a = [prob.addVar(f"a_{i}", 0, 1) for i in range(4)]
+ b = [prob.addVar(f"b_{i}", 0, 1) for i in range(3)]
+ for i, v in enumerate([x, y, z, w, u] + a + b):
+ prob.cost += v * float(i + 1)
+ prob.cost += w * 3.0 # same column twice in the objective
+ prob.add_constraint("le", prob.sum([x, y]) <= 1)
+ prob.add_constraint("ge", prob.sum([x, y, z]) >= 1)
+ prob.add_constraint("eq", prob.sum(a) == 2)
+ prob.add_constraint("lhs_const", x + 2 <= 5)
+ prob.add_constraint("both_sides", a[0] >= 0.25 * prob.sum(b))
+ prob.add_constraint("both_sides2", prob.sum(a) <= prob.sum(b) + 1)
+ prob.add_constraint("dup", prob.sum([a[2], a[3], a[2]]) <= 1)
+ prob.add_constraint("cancel", x - x + y == 0)
+ prob.add_constraint("empty", prob.sum([]) <= 3)
+ prob.add_constraint("free_lo", w >= -7)
+ prob.add_constraint("free_hi", w - z <= 4)
+ prob.add_constraint("neg", prob.sum([a[1]] + [-v for v in b]) == 0)
+ prob.add_lazy_constraint("lazy", prob.sum([a[2], b[2]]) <= 1)
+ prob.add_constraint("scaled",
+ prob.sum([v * t for v, t in zip(b, [900.0, 450.0, 1800.0])])
+ >= 900.0 * y)
+ prob.add_constraint("eq_var", x == 1)
+ return dict(x=x, y=y, z=z, w=w, u=u, a=a, b=b)
+
+
+def scen_no_constraints(prob):
+ vs = [prob.addVar(f"v_{i}", 0, 1) for i in range(3)]
+ for i, v in enumerate(vs):
+ prob.cost += v * float(i + 1)
+ return dict(vs=vs)
+
+
+def scen_single_var(prob):
+ x = prob.addVar("x", 0, 1)
+ prob.cost += x * 2.0
+ prob.add_constraint("c", x >= 1)
+ return dict(x=x)
+
+
+def scen_objective_untouched(prob):
+ # cost stays the bare variable handed out in __init__
+ x = prob.addVar("x", 0, 1)
+ prob.add_constraint("c", x <= 1)
+ return dict(x=x)
+
+
+def scen_bounds_before_solve(prob):
+ h = scen_mixed(prob)
+ # changeVarBounds in the middle of the build: the reference applies it to
+ # HiGHS at once, HighsProblem creates the still-pending column with the
+ # new bounds. Both must end up with the same model.
+ prob.changeVarBounds(h["u"], lower=3)
+ prob.add_constraint("late", prob.sum([h["a"][3], h["b"][0]]) <= 1)
+ prob.changeVarBounds(h["b"][1], upper=0)
+ return h
+
+
+def scen_assignment_small(prob):
+ """Three targets, two cobras, hand-checkable optimum.
+
+ Arc costs: t0-c0 1, t0-c1 5, t1-c0 2, t2-c1 3; leaving a target
+ unobserved (its sink arc) costs 10. Best: t0 on c0, t2 on c1, t1 unobserved
+ -> 1 + 3 + 10 = 14, and no other assignment reaches 14.
+ """
+ arcs = {}
+ for (t, c, cost) in [(0, 0, 1.0), (0, 1, 5.0), (1, 0, 2.0), (2, 1, 3.0)]:
+ f = prob.addVar(f"Tv_Cv_{t}_{c}", 0, 1)
+ prob.cost += f * cost
+ arcs[(t, c)] = f
+ sinks = []
+ for t in range(3):
+ s = prob.addVar(f"ST_sink_{t}", 0, 1)
+ prob.cost += s * 10.0
+ sinks.append(s)
+ for t in range(3):
+ flows = [f for (tt, c), f in arcs.items() if tt == t] + [sinks[t]]
+ prob.add_constraint(f"TvIO_{t}", prob.sum(flows) == 1)
+ for c in range(2):
+ flows = [f for (t, cc), f in arcs.items() if cc == c]
+ prob.add_constraint(f"Cvlim_{c}", prob.sum(flows) <= 1)
+ return dict(arcs=arcs, sinks=sinks, optimum=14.0,
+ assigned={"Tv_Cv_0_0", "Tv_Cv_2_1"})
+
+
+def scen_assignment_random(prob, ntargets=60, ncobras=25, seed=5):
+ """Flow-shaped model like buildProblem's: targets, cobras, collision pairs.
+
+ Random continuous costs make the optimum unique in practice, so two exact
+ solvers must agree on the assignment, not only on the objective.
+ """
+ rng = np.random.default_rng(seed)
+ arcs, by_cobra, by_target, all_arcs = {}, [[] for _ in range(ncobras)], [], []
+ for t in range(ntargets):
+ cobras = rng.choice(ncobras, size=int(rng.integers(1, 4)), replace=False)
+ tarcs = []
+ for c in cobras:
+ f = prob.addVar(f"Tv_Cv_{t}_{c}", 0, 1)
+ prob.cost += f * float(rng.random())
+ arcs[(t, int(c))] = f
+ by_cobra[c].append(f)
+ tarcs.append(f)
+ all_arcs.append(f)
+ s = prob.addVar(f"ST_sink_{t}", 0, 1)
+ prob.cost += s * 10.0
+ tarcs.append(s)
+ by_target.append(tarcs)
+ for p in range(ntargets // 2):
+ i, j = rng.integers(0, len(all_arcs), size=2)
+ if i != j:
+ prob.add_lazy_constraint(f"Coll_{p}", prob.sum([all_arcs[i], all_arcs[j]]) <= 1)
+ for c in range(ncobras):
+ if by_cobra[c]:
+ prob.add_constraint(f"Cvlim_{c}", prob.sum(by_cobra[c]) <= 1)
+ for t in range(ntargets):
+ prob.add_constraint(f"TvIO_{t}", prob.sum(by_target[t]) == 1)
+ return dict(arcs=arcs)
+
+
+def scen_random_rows(prob, n=300, seed=7):
+ """Random rows of every shape, duplicates included; exercises the sort."""
+ rng = np.random.default_rng(seed)
+ vs = [prob.addVar(f"v_{i}", 0, 1) for i in range(n)]
+ for v in vs:
+ prob.cost += v * float(rng.random())
+ for k in range(n):
+ m = int(rng.integers(1, 6))
+ idx = rng.integers(0, n, size=m) # duplicates allowed
+ expr = prob.sum([vs[i] * float(rng.integers(-2, 3)) for i in idx])
+ kind = k % 3
+ if kind == 0:
+ prob.add_constraint(f"r_{k}", expr <= float(rng.integers(0, 3)))
+ elif kind == 1:
+ prob.add_constraint(f"r_{k}", expr >= float(-rng.integers(0, 3)))
+ else:
+ prob.add_lazy_constraint(f"r_{k}", expr <= prob.sum(vs[:2]) + 1)
+ return dict(vs=vs)
+
+
+SCENARIOS = {
+ "mixed": scen_mixed,
+ "no_constraints": scen_no_constraints,
+ "single_var": scen_single_var,
+ "objective_untouched": scen_objective_untouched,
+ "bounds_before_solve": scen_bounds_before_solve,
+ "assignment_small": scen_assignment_small,
+ "assignment_random": scen_assignment_random,
+ "random_rows": scen_random_rows,
+}
+
+
+# ---------------------------------------------------------------------------
+# helpers
+# ---------------------------------------------------------------------------
+def assert_same_highs_model(ref, new):
+ """Element-by-element comparison of the two HiGHS-side models."""
+ a, b = ref._prob, new._prob
+ assert (a.getNumCol(), a.getNumRow(), a.getNumNz()) == \
+ (b.getNumCol(), b.getNumRow(), b.getNumNz())
+ la, lb = a.getLp(), b.getLp()
+ for attr in ("col_cost_", "col_lower_", "col_upper_", "row_lower_", "row_upper_"):
+ assert np.array_equal(np.asarray(getattr(la, attr)),
+ np.asarray(getattr(lb, attr))), attr
+ assert la.a_matrix_.format_ == lb.a_matrix_.format_
+ for attr in ("start_", "index_", "value_"):
+ assert np.array_equal(np.asarray(getattr(la.a_matrix_, attr)),
+ np.asarray(getattr(lb.a_matrix_, attr))), attr
+ assert list(la.integrality_) == list(lb.integrality_)
+ assert list(la.col_names_) == list(lb.col_names_)
+ assert list(la.row_names_) == list(lb.row_names_)
+
+
+def assert_same_dump(ref, new, tmp_path, tag):
+ for ext in ("lp", "mps"):
+ fa, fb = tmp_path / f"{tag}_ref.{ext}", tmp_path / f"{tag}_new.{ext}"
+ ref.dump(str(fa))
+ new.dump(str(fb))
+ assert filecmp.cmp(fa, fb, shallow=False), f".{ext} differs: {fa} {fb}"
+
+
+def objective(prob):
+ if isinstance(prob, GurobiProblem):
+ return prob._prob.ObjVal
+ return prob._prob.getObjectiveValue()
+
+
+def assigned(prob):
+ return {k for k, v in prob._vardict.items()
+ if k.startswith("Tv_Cv_") and prob.value(v) > 0.5}
+
+
+# ---------------------------------------------------------------------------
+# 2. model identity against highspy's own addConstr path
+# ---------------------------------------------------------------------------
+@pytest.mark.parametrize("name", sorted(SCENARIOS))
+def test_model_matches_highspy_addconstr_path(name, tmp_path):
+ ref, new = ReferenceHighs(), HighsProblem()
+ SCENARIOS[name](ref)
+ SCENARIOS[name](new)
+ assert_same_dump(ref, new, tmp_path, name) # before solve: no objective
+ assert_same_highs_model(ref, new) # names reached HiGHS in dump()
+ ref.solve()
+ new.solve()
+ assert_same_dump(ref, new, tmp_path, name + "_solved")
+ assert (ref._prob.modelStatusToString(ref._prob.getModelStatus())
+ == new._prob.modelStatusToString(new._prob.getModelStatus()))
+ assert abs(objective(ref) - objective(new)) <= 1e-9 * max(1.0, abs(objective(new)))
+
+
+def test_dump_twice_and_incremental(tmp_path):
+ ref, new = ReferenceHighs(), HighsProblem()
+ hr, hn = scen_mixed(ref), scen_mixed(new)
+ assert_same_dump(ref, new, tmp_path, "d1")
+ assert_same_dump(ref, new, tmp_path, "d2")
+ assert filecmp.cmp(tmp_path / "d1_new.lp", tmp_path / "d2_new.lp", shallow=False)
+ # add more columns and rows after a dump, then dump again
+ for prob, h in ((ref, hr), (new, hn)):
+ q = prob.addVar("q", 0, 3)
+ prob.cost += q * 0.5
+ prob.add_constraint("after_dump", prob.sum([q, h["x"]]) >= 2)
+ assert_same_dump(ref, new, tmp_path, "d3")
+ new.update()
+ assert_same_highs_model(ref, new)
+ # and the incremental build equals a one-shot build of the same model
+ fresh = HighsProblem()
+ h = scen_mixed(fresh)
+ q = fresh.addVar("q", 0, 3)
+ fresh.cost += q * 0.5
+ fresh.add_constraint("after_dump", fresh.sum([q, h["x"]]) >= 2)
+ fresh.dump(str(tmp_path / "fresh.lp"))
+ assert filecmp.cmp(tmp_path / "d3_new.lp", tmp_path / "fresh.lp", shallow=False)
+
+
+def test_change_bounds_after_solve(tmp_path):
+ ref, new = ReferenceHighs(), HighsProblem()
+ hr, hn = scen_mixed(ref), scen_mixed(new)
+ ref.solve()
+ new.solve()
+ obj1 = objective(new)
+ for prob, h in ((ref, hr), (new, hn)):
+ prob.changeVarBounds(h["u"], lower=4)
+ prob.changeVarBounds(h["a"][3], upper=0)
+ prob.changeVarBounds(h["w"], lower=-2, upper=10)
+ prob.solve()
+ assert_same_dump(ref, new, tmp_path, "cb")
+ assert abs(objective(ref) - objective(new)) <= 1e-9 * max(1.0, abs(objective(new)))
+ assert objective(new) != obj1
+ for x, y in zip(ref._vardict.values(), new._vardict.values()):
+ assert ref.varBounds(x) == new.varBounds(y)
+
+
+# ---------------------------------------------------------------------------
+# 3. solutions: known optimum, and agreement with Gurobi
+# ---------------------------------------------------------------------------
+def test_known_optimum():
+ new = make_highs()
+ h = scen_assignment_small(new)
+ new.solve()
+ assert abs(objective(new) - h["optimum"]) < 1e-9
+ assert assigned(new) == h["assigned"]
+ assert new.value(h["sinks"][1]) == 1.0
+ assert new.value(h["sinks"][0]) == 0.0 and new.value(h["sinks"][2]) == 0.0
+
+
+SOLVABLE = ["mixed", "no_constraints", "single_var", "objective_untouched",
+ "bounds_before_solve", "assignment_small", "assignment_random"]
+
+
+@needs_gurobi
+@pytest.mark.parametrize("name", SOLVABLE)
+def test_agrees_with_gurobi(name):
+ g, h = make_gurobi(), make_highs()
+ SCENARIOS[name](g)
+ SCENARIOS[name](h)
+ g.solve()
+ h.solve()
+ og, oh = objective(g), objective(h)
+ assert abs(og - oh) <= 1e-9 * max(1.0, abs(og)), (og, oh)
+ assert list(g._vardict) == list(h._vardict)
+ for a, b in zip(g._vardict.values(), h._vardict.values()):
+ assert g.varBounds(a) == h.varBounds(b)
+ if name.startswith("assignment"):
+ # random continuous costs: the optimum is unique, so the solvers must
+ # pick the same arcs, not just reach the same value
+ assert assigned(g) == assigned(h)
+
+
+@needs_gurobi
+def test_agrees_with_gurobi_after_bound_changes():
+ g, h = make_gurobi(), make_highs()
+ scen_assignment_random(g)
+ scen_assignment_random(h)
+ g.solve()
+ h.solve()
+ # forbid two arcs of the first solution and re-solve
+ forbid = sorted(assigned(h))[:2]
+ for prob in (g, h):
+ for name in forbid:
+ prob.changeVarBounds(prob.varByName(name), upper=0)
+ prob.solve()
+ assert abs(objective(g) - objective(h)) <= 1e-9 * max(1.0, abs(objective(g)))
+ assert assigned(g) == assigned(h)
+ assert not (assigned(h) & set(forbid))
+
+
+# ---------------------------------------------------------------------------
+# 1. HighsProblem behaviour
+# ---------------------------------------------------------------------------
+def test_infeasible_raises():
+ new = HighsProblem()
+ x = new.addVar("x", 0, 1)
+ y = new.addVar("y", 0, 1)
+ new.cost += x * 1.0 + y * 1.0
+ new.add_constraint("c", new.sum([x, y]) >= 3)
+ with pytest.raises(RuntimeError):
+ new.solve()
+ with pytest.raises(RuntimeError): # no all-zero "solution" either
+ new.value(x)
+
+
+def test_unbounded_expression_rejected():
+ new = HighsProblem()
+ x = new.addVar("x", 0, 1)
+ y = new.addVar("y", 0, 1)
+ with pytest.raises(Exception):
+ new.add_constraint("bad", new.sum([x, y]))
+ # nothing half-added: the model still builds and solves
+ new.add_constraint("ok", new.sum([x, y]) <= 1)
+ new.solve()
+ assert new._prob.getNumRow() == 1
+ assert "bad" not in new._constraintdict
+
+
+def test_names_reach_highs_in_order(tmp_path):
+ new = HighsProblem()
+ scen_mixed(new)
+ new.solve() # so the objective is in the file as well
+ new.dump(str(tmp_path / "n.lp"))
+ lp = new._prob.getLp()
+ assert list(lp.col_names_) == ["cost"] + list(new._vardict)
+ assert list(lp.row_names_) == list(new._constraintdict)
+ text = (tmp_path / "n.lp").read_text()
+ for name in ("both_sides2", "eq_var", "a_3"):
+ assert name in text
+ obj_line = text.split("obj:", 1)[1].split("\n", 1)[0]
+ assert "x" in obj_line and "w" in obj_line, obj_line
+
+
+def test_dump_objective_only_after_solve(tmp_path):
+ """The objective reaches HiGHS in solve(), so a dump written before
+ solve() has an empty objective -- exactly what GurobiProblem.dump()
+ writes before solve() ("Minimize 0 cost").
+ """
+ new = HighsProblem()
+ x = new.addVar("x", 0, 1)
+ y = new.addVar("y", 0, 1)
+ new.cost += x * 2.0 + y * 3.0
+ new.cost += x * 0.5
+ new.add_constraint("c", new.sum([x, y]) >= 1)
+ new.dump(str(tmp_path / "pre.lp"))
+ pre = (tmp_path / "pre.lp").read_text()
+ obj_line = pre.split("obj:", 1)[1].split("\n", 1)[0]
+ assert obj_line.strip() == "", obj_line
+ new.solve()
+ new.dump(str(tmp_path / "post.lp"))
+ post = (tmp_path / "post.lp").read_text()
+ # "+1 cost": the objective is accumulated onto the free `cost` column, so
+ # that column itself carries coefficient 1 -- same as with Gurobi.
+ assert "obj: +1 cost +2.5 x +3 y" in post, post
+ h = highspy.Highs()
+ h.setOptionValue("output_flag", False)
+ h.readModel(str(tmp_path / "post.lp"))
+ h.run()
+ assert abs(h.getObjectiveValue() - objective(new)) < 1e-9
+
+
+def test_row_flush_count():
+ new = HighsProblem()
+ scen_mixed(new)
+ assert new._nrow_flushes == 0
+ assert new._prob.getNumRow() == 0 # rows still buffered
+ new.solve()
+ assert new._nrow_flushes == 1
+ assert new._prob.getNumRow() == len(new._constraintdict)
+ new.solve() # nothing pending: no new flush
+ new.update()
+ assert new._nrow_flushes == 1
+
+
+def test_lookup_by_name():
+ ref, new = ReferenceHighs(), HighsProblem()
+ scen_mixed(ref)
+ scen_mixed(new)
+ for name in ("x", "w", "a_2"):
+ assert new.varByName(name) is new._vardict[name]
+ assert new.varByName(name).index == ref._vardict[name].index
+ for name in ("dup", "cancel", "empty"):
+ ca, cb = ref._constraintdict[name], new.constraintByName(name)
+ assert (ca.idxs, ca.vals, ca.bounds) == (cb.idxs, cb.vals, cb.bounds)
+
+
+def test_value_without_solution_raises():
+ """value() refuses whenever the cached solution is invalid.
+
+ Invalidation follows Gurobi's timing: a change that has reached the
+ solver (changeVarBounds, or a flush through update()/dump()/solve())
+ discards the solution; a change that is still buffered leaves the old
+ solution readable, just as a pending Gurobi modification leaves var.X.
+ """
+ new = HighsProblem()
+ x = new.addVar("x", 0, 1)
+ y = new.addVar("y", 0, 1)
+ new.cost += x * 1.0 + y * 2.0
+ new.add_constraint("c", new.sum([x, y]) >= 1)
+ with pytest.raises(RuntimeError): # nothing flushed yet
+ new.value(x)
+ new.update()
+ with pytest.raises(RuntimeError): # flushed, not solved
+ new.value(x)
+ new.solve()
+ assert new.value(x) == 1.0 and new.value(y) == 0.0
+ new.changeVarBounds(x, upper=0) # pending: old solution readable
+ assert new.value(x) == 1.0
+ new.update() # applied: stale
+ with pytest.raises(RuntimeError):
+ new.value(x)
+ new.solve()
+ assert new.value(x) == 0.0 and new.value(y) == 1.0
+ new.add_constraint("late", new.sum([x, y]) <= 5)
+ assert new.value(y) == 1.0 # row still buffered: readable
+ new.update() # row reached HiGHS: stale
+ with pytest.raises(RuntimeError):
+ new.value(y)
+ new.solve()
+ assert new.value(y) == 1.0
+ z = new.addVar("z", 0, 1) # column still buffered
+ new.cost += z * 1.0
+ assert new.value(y) == 1.0 # old column: readable
+ with pytest.raises(RuntimeError): # new column: no value yet
+ new.value(z)
+ new.update() # column reached HiGHS: stale
+ with pytest.raises(RuntimeError):
+ new.value(y)
+ new.solve()
+ assert new.value(z) == 0.0
+
+
+@needs_gurobi
+def test_value_rule_matches_gurobi():
+ """Same sequence on both backends: value() is readable or raises together."""
+ def probe(prob, var):
+ try:
+ return ("value", float(prob.value(var)))
+ except Exception: # GurobiError / RuntimeError
+ return ("raises",)
+
+ steps = []
+ for prob in (make_gurobi(), make_highs()):
+ x = prob.addVar("x", 0, 1)
+ y = prob.addVar("y", 0, 1)
+ prob.cost += x * 1.0 + y * 2.0
+ prob.add_constraint("c", prob.sum([x, y]) >= 1)
+ prob.update()
+ seq = [probe(prob, x)[0]] # before solve
+ prob.solve()
+ seq.append(probe(prob, x)) # after solve
+ prob.changeVarBounds(x, upper=0)
+ seq.append(probe(prob, x)[0]) # pending change
+ prob.update()
+ seq.append(probe(prob, x)[0]) # applied change
+ prob.solve()
+ seq.append(probe(prob, x)) # re-solved
+ prob.add_constraint("late", prob.sum([x, y]) <= 5)
+ seq.append(probe(prob, y)[0]) # pending row
+ prob.update()
+ seq.append(probe(prob, y)[0]) # applied row
+ steps.append(seq)
+ assert steps[0] == steps[1], steps