From 655e2e6c583f4bd9ffca84ccb9af663e4bc618d7 Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Sun, 23 Aug 2026 02:47:24 -1000 Subject: [PATCH 1/6] 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/6] 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/6] 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/6] 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/6] 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 From 80434108c79f0d54ae4041ad2eeff3b42850dbcc Mon Sep 17 00:00:00 2001 From: Masato Onodera Date: Tue, 8 Sep 2026 17:25:29 -1000 Subject: [PATCH 6/6] Batch row additions in the HiGHS backend and match Gurobi's read rules HighsProblem.add_constraint used to hand every row to Highs.addConstr, about 19 us per row (highspy's unique_elements 9.5, addRow 4.8, passRowName 0.9, wrapper 3), i.e. 10-17 s per million rows on the real instances. Rows are now collected in flat typed arrays (_RowBuffer) and handed to HiGHS in one addRows call at the next flush, ordered within each row like unique_elements does, so the model HiGHS stores and the .lp/.mps it writes are byte-identical to the one-row-at-a-time path. About 2 us per row; a 625k-row synthetic model builds its rows in 1.2 s instead of 12.5 s, a 147k-row real instance in 18.8 s instead of 21.6 s. Column and row names are passed to HiGHS in dump() only (HiGHS's own log refers to rows and columns by index either way). Two behaviours now follow GurobiProblem: value() raises when there is no valid solution -- before solve(), after a failed solve(), or once a change (new columns or rows, changeVarBounds) has been applied by update(), dump() or solve(); a still-pending change leaves the previous solution readable, as a pending Gurobi modification leaves var.X. changeVarBounds itself is pending until the next flush, like a Gurobi bound change until update(). dump() before solve() writes an empty objective, as GurobiProblem.dump() does. Docstring lists the HiGHS options relevant here with their defaults, notes that varBounds/changeVarBounds are instance methods on this backend, and that highspy >= 1.15 is required. tests/test_highs_problem.py (pytest) checks the behaviour, identity with highspy's own addVariable/addConstr path, known optima, and agreement with GurobiProblem (skipped without gurobipy). misc/bench_highs_problem.py is the synthetic benchmark; misc/compare_highs_real_instance.py compares solver="gurobi" and solver="highs" through buildProblem on a real instance. Co-Authored-By: Claude Fable 5.1 --- ets_fiber_assigner/netflow.py | 258 ++++++++++-- misc/bench_highs_problem.py | 152 +++++++ misc/compare_highs_real_instance.py | 331 ++++++++++++++++ tests/test_highs_problem.py | 589 ++++++++++++++++++++++++++++ 4 files changed, 1290 insertions(+), 40 deletions(-) create mode 100644 misc/bench_highs_problem.py create mode 100644 misc/compare_highs_real_instance.py create mode 100644 tests/test_highs_problem.py diff --git a/ets_fiber_assigner/netflow.py b/ets_fiber_assigner/netflow.py index eb23c9f..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 @@ -767,6 +769,64 @@ 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. @@ -774,6 +834,29 @@ class HighsProblem(LPProblem): 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): @@ -789,11 +872,20 @@ def __init__(self, name="problem", extraOptions=None): for key, value in extraOptions.items(): self._prob.setOptionValue(key, value) + # Columns: reserved in _newCol(), created in _flush(). self._ncols = 0 - self._pending = [] # (name, lb, ub, is_integer) awaiting _flush() - self._bounds = {} # column index -> (lb, ub), for varBounds() + 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. @@ -807,7 +899,8 @@ def _newCol(self, name, lb, ub, is_integer): 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._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)) @@ -815,7 +908,8 @@ def _newCol(self, name, lb, ub, is_integer): return var def _flush(self): - """Create every reserved column in one call. Safe to call at any time. + """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 @@ -823,34 +917,81 @@ def _flush(self): 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: + if not self._pending and not self._bounds_changed: 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]) + 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 @@ -863,9 +1004,25 @@ def addVar(self, name, lo, hi): return var def add_constraint(self, name, constraint): - self._flush() + 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 - self._prob.addConstr(constraint, name=name) + 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. @@ -887,9 +1044,21 @@ def value(self, var): 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: - return self._prob.val(var) + 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): @@ -899,8 +1068,6 @@ def _cacheSolution(self): 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 - self._colvals = np.asarray(self._prob.getSolution().col_value) def _checkSolved(self): @@ -926,31 +1093,42 @@ def _checkSolved(self): + self._prob.modelStatusToString(status)) def solve(self): - self._flush() + 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._flush() + self._flushAll() def dump(self, filename): - self._flush() + 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): - self._flush() + """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)) - self._colvals = None - self._prob.changeColBounds(var.index, lb, ub) + if var.index < self._ncols - len(self._pending): + self._bounds_changed.add(var.index) class Telescope(object): 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