opt

Linear, mixed-integer, and convex quadratic programs built from cells in a sheet and solved through the HiGHS-backed _opt extension. See the Optimization guide for the user-facing commands.

opt

Sheet-level optimization.

Builds a linear, mixed-integer, or convex quadratic program from cells in a Grid and solves it via the HiGHS-backed _opt extension. The user-facing model is sheet-resident:

• One objective cell containing a linear formula (e.g. =3*A1+5*A2).
• A list of decision variable cells. They must hold numeric values (or be empty); formula cells are refused so the optimizer doesn't silently overwrite live computations.
• A list of constraint cells, each containing a comparison formula (e.g. =A1+A2<=10). Their current evaluated values (True/False) indicate live feasibility; the optimizer reads the underlying AST.

Linearity is enforced by walking gridcalc's formula AST. Cell references that resolve to decision variables become coefficients. A formula cell that reads a decision variable, directly or through other formulas, is inlined: its AST is walked in place of its stale value. Everything else is folded into the constant term using the cell's currently evaluated value, so non-decision cells act as parameters: edit them and re-run. A text parameter is refused.

Supported AST shapes

Number, CellRef, BinOp(+,-,*,/), UnaryOp(+,-), Percent, Call("SUM", RangeRef|expr), parenthesized expressions.

Anything else (Bool, String, ErrorLit, other Call/PyCall, RangeRef outside SUM, Name) raises NotLinear with a message naming the offending node.

OptModel dataclass

OptModel(
    sense: str,
    objective: str,
    vars: str,
    constraints: str,
    bounds: str = "",
    integers: str = "",
    binaries: str = "",
)

A persisted LP model definition stored in the workbook.

The fields hold the string specs the user typed for each component ("A4:A5", "D4:D6", "A1=-inf:10"), not pre-parsed cell coordinates. This preserves the user's range/list intent verbatim through save/load round-trips, mirrors how named ranges are stored, and defers cell-ref resolution (and any errors it would produce) to :opt run time.

SolveResult dataclass

SolveResult(
    status: int,
    status_name: str,
    objective: float,
    values: dict[CellKey, float],
    applied: bool,
    sensitivity: Sensitivity | None = None,
    conflict: list[CellKey] | None = None,
    unbounded: list[CellKey] | None = None,
    quadratic: bool = False,
)

Sensitivity dataclass

Sensitivity(
    variables: list[VarSensitivity],
    constraints: list[ConstraintSensitivity],
)

VarSensitivity dataclass

VarSensitivity(
    cell: CellKey,
    value: float,
    reduced_cost: float,
    obj_coef: float,
    obj_from: float,
    obj_till: float,
)

Sensitivity of the optimum to one decision variable.

reduced_cost is the amount the objective would change per unit if the variable were forced away from its bound; it is zero for any variable already in the basis. obj_from / obj_till bracket the range over which this variable's objective coefficient can move without changing the optimal basis (the values themselves would change, the choice of which variables are non-zero would not).

ConstraintSensitivity dataclass

ConstraintSensitivity(
    cell: CellKey,
    shadow_price: float,
    rhs: float,
    activity: float,
    slack: float,
    binding: bool,
    rhs_from: float,
    rhs_till: float,
)

Sensitivity of the optimum to one constraint.

shadow_price is the marginal change in the objective per unit relaxation of the right-hand side -- the value of one more unit of this resource, and the number a user actually wants when deciding what to buy more of. It is valid only within rhs_from..rhs_till; past those limits the basis changes and the price no longer applies.

SweepPoint dataclass

SweepPoint(
    rhs: float,
    status_name: str,
    objective: float,
    shadow_price: float | None,
    delta: float | None,
    breakpoint: bool,
)

One re-solve of the model at a substituted right-hand side.

LinearForm dataclass

LinearForm(
    coeffs: dict[CellKey, float] = dict(),
    constant: float = 0.0,
)

A sum of (coefficient * decision_var) terms plus a constant.

coeffs is sparse: missing keys are zero. Two LinearForms can be added, subtracted, scaled, and negated to compose larger expressions.

QuadForm dataclass

QuadForm(
    quad: dict[tuple[CellKey, CellKey], float] = dict(),
    linear: dict[CellKey, float] = dict(),
    constant: float = 0.0,
)

A degree-<=2 polynomial over decision variables.

quad is keyed by an ordered pair of cells: (k, k) is a squared term, (j, k) with j != k is a cross term. Cross terms are tracked rather than rejected on sight so the error message can name the pair that made the objective non-separable.

mul

mul(other: QuadForm) -> QuadForm

Multiply, refusing anything that would exceed degree 2.

Source code in src/gridcalc/opt.py

def mul(self, other: QuadForm) -> QuadForm:
    """Multiply, refusing anything that would exceed degree 2."""
    if not self.is_linear and not other.is_constant:
        raise NotQuadratic("objective is degree 3 or higher")
    if not other.is_linear and not self.is_constant:
        raise NotQuadratic("objective is degree 3 or higher")

    out = QuadForm({}, {}, self.constant * other.constant)
    for pair, qv in self.quad.items():
        out.quad[pair] = out.quad.get(pair, 0.0) + qv * other.constant
    for pair, qv in other.quad.items():
        out.quad[pair] = out.quad.get(pair, 0.0) + qv * self.constant
    for cell, lv in self.linear.items():
        out.linear[cell] = out.linear.get(cell, 0.0) + lv * other.constant
    for cell, lv in other.linear.items():
        out.linear[cell] = out.linear.get(cell, 0.0) + lv * self.constant
    # The outer product of the two linear parts is where squares and
    # cross terms come from.
    for a, va in self.linear.items():
        for b, vb in other.linear.items():
            key = (a, b) if a <= b else (b, a)
            out.quad[key] = out.quad.get(key, 0.0) + va * vb
    return out

squares

squares() -> dict[CellKey, float]

Diagonal terms only. Raises if the objective has cross terms.

Source code in src/gridcalc/opt.py

def squares(self) -> dict[CellKey, float]:
    """Diagonal terms only. Raises if the objective has cross terms."""
    for (a, b), v in self.quad.items():
        if a != b and v != 0.0:
            raise NotQuadratic(f"objective couples {_cellname(*a)} and {_cellname(*b)}")
    return {a: v for (a, b), v in self.quad.items() if a == b and v != 0.0}

hessian

hessian(order: list[CellKey]) -> list[list[float]]

Dense lower triangle of Q for the 0.5 * x' Q x objective form.

Returned ragged and row-major: row i holds columns 0..i. The solver takes cross terms directly, so this no longer has to reject them -- q[(a, b)] * x_a * x_b with a != b contributes q to both symmetric halves, which in the 0.5 x'Qx convention means Q[a][b] = Q[b][a] = q. A squared term q * x_a^2 needs Q[a][a] = 2q for the same reason.

Source code in src/gridcalc/opt.py

def hessian(self, order: list[CellKey]) -> list[list[float]]:
    """Dense lower triangle of Q for the ``0.5 * x' Q x`` objective form.

    Returned ragged and row-major: row ``i`` holds columns ``0..i``. The
    solver takes cross terms directly, so this no longer has to reject
    them -- ``q[(a, b)] * x_a * x_b`` with ``a != b`` contributes ``q`` to
    both symmetric halves, which in the ``0.5 x'Qx`` convention means
    ``Q[a][b] = Q[b][a] = q``. A squared term ``q * x_a^2`` needs
    ``Q[a][a] = 2q`` for the same reason.
    """
    index = {cell: i for i, cell in enumerate(order)}
    tri: list[list[float]] = [[0.0] * (i + 1) for i in range(len(order))]
    for (a, b), v in self.quad.items():
        if v == 0.0:
            continue
        if a not in index or b not in index:
            continue
        ia, ib = index[a], index[b]
        if ia == ib:
            tri[ia][ia] += 2.0 * v
        else:
            lo, hi = (ia, ib) if ia < ib else (ib, ia)
            tri[hi][lo] += v
    return tri

OptError

Bases: Exception

Caller-facing error: malformed model, bad cell selection, etc.

NotLinear

Bases: OptError

A formula cannot be expressed as a linear combination of decision vars.

NotQuadratic

Bases: NotLinear

An objective is nonlinear in a way this optimizer cannot express.

Deliberately a subclass of :class:NotLinear rather than a sibling. Both walkers are the safety boundary for the optimizer -- each accepts a closed whitelist of AST nodes and rejects everything else, so nothing that could be a sandbox concern (Name, PyCall, attribute-style Call) reaches an evaluation path. Callers written against that guarantee catch NotLinear; widening the objective walker to degree 2 must not quietly slip past those handlers.

solve

solve(
    grid: Grid,
    objective_cell: CellKey,
    decision_vars: list[CellKey],
    constraint_cells: list[CellKey],
    *,
    maximize: bool = True,
    bounds: dict[CellKey, tuple[float, float]]
    | None = None,
    integer_vars: set[CellKey] | None = None,
    binary_vars: set[CellKey] | None = None,
    apply: bool = True,
    sensitivity: bool = False,
    diagnose: bool = False,
    rhs_override: dict[CellKey, float] | None = None,
) -> SolveResult

Build an LP (or MIP) from the named cells, solve, and (by default) write back.

The objective cell must contain a formula. Decision-variable cells must NOT contain formulas (they get overwritten on success). Each constraint cell must contain a formula whose root is a comparison operator.

integer_vars and binary_vars are subsets of decision_vars; cells in either set are flagged as integer or binary respectively, which routes the solve through branch-and-bound. Binary cells have their bounds clamped to [0,1] regardless of bounds; a cell appearing in both sets raises OptError.

sensitivity=True additionally returns shadow prices, reduced costs, and ranging information in SolveResult.sensitivity. It is silently ignored for MIPs (the field stays None): branch-and-bound duals describe one LP relaxation rather than the integer problem, so reporting them would be actively misleading.

rhs_override replaces the right-hand side of the named constraint cells for this solve only, without touching the sheet. The constraint's coefficients still come from its formula; only the constant moves. This is what makes what-if analysis possible without rewriting cells and recalculating -- see sweep.

diagnose=True explains a failed solve. On INFEASIBLE it populates SolveResult.conflict with a minimal set of contradictory constraint cells (one extra solve per constraint); on UNBOUNDED it populates SolveResult.unbounded with the decision cells that can grow without limit (two extra solves). Neither runs on a successful solve.

Source code in src/gridcalc/opt.py

def solve(
    grid: Grid,
    objective_cell: CellKey,
    decision_vars: list[CellKey],
    constraint_cells: list[CellKey],
    *,
    maximize: bool = True,
    bounds: dict[CellKey, tuple[float, float]] | None = None,
    integer_vars: set[CellKey] | None = None,
    binary_vars: set[CellKey] | None = None,
    apply: bool = True,
    sensitivity: bool = False,
    diagnose: bool = False,
    rhs_override: dict[CellKey, float] | None = None,
) -> SolveResult:
    """Build an LP (or MIP) from the named cells, solve, and (by default) write back.

    The objective cell must contain a formula. Decision-variable cells must
    NOT contain formulas (they get overwritten on success). Each constraint
    cell must contain a formula whose root is a comparison operator.

    ``integer_vars`` and ``binary_vars`` are subsets of ``decision_vars``;
    cells in either set are flagged as integer or binary respectively, which
    routes the solve through branch-and-bound. Binary cells have their bounds
    clamped to [0,1] regardless of ``bounds``; a cell appearing in both sets
    raises ``OptError``.

    ``sensitivity=True`` additionally returns shadow prices, reduced costs,
    and ranging information in ``SolveResult.sensitivity``. It is silently
    ignored for MIPs (the field stays ``None``): branch-and-bound duals
    describe one LP relaxation rather than the integer problem, so reporting
    them would be actively misleading.

    ``rhs_override`` replaces the right-hand side of the named constraint
    cells for this solve only, without touching the sheet. The constraint's
    coefficients still come from its formula; only the constant moves. This
    is what makes what-if analysis possible without rewriting cells and
    recalculating -- see ``sweep``.

    ``diagnose=True`` explains a failed solve. On INFEASIBLE it populates
    ``SolveResult.conflict`` with a minimal set of contradictory constraint
    cells (one extra solve per constraint); on UNBOUNDED it populates
    ``SolveResult.unbounded`` with the decision cells that can grow without
    limit (two extra solves). Neither runs on a successful solve.
    """
    if not decision_vars:
        raise OptError("at least one decision variable is required")
    if len(set(decision_vars)) != len(decision_vars):
        raise OptError("decision variables must be unique")

    var_set = set(decision_vars)
    var_index = {v: i for i, v in enumerate(decision_vars)}
    n = len(decision_vars)

    # Reject formula decision cells up-front so the operator never silently
    # destroys live computation. Override for advanced use cases isn't
    # supported yet (would need a flag and an undo guarantee).
    for c, r in decision_vars:
        cell = grid.cells[c][r]
        if cell.type == FORMULA:
            raise OptError(
                f"decision cell {_cellname(c, r)} contains a formula; "
                "decision variables must hold values (or be empty)"
            )
        if cell.type not in (EMPTY, NUM):
            raise OptError(f"decision cell {_cellname(c, r)} must be numeric or empty")

    # Objective.
    obj_c, obj_r = objective_cell
    obj_cell = grid.cells[obj_c][obj_r]
    obj_ast = _cell_ast(obj_cell) if obj_cell.type == FORMULA else None
    if obj_cell.type != FORMULA or obj_ast is None:
        raise OptError(f"objective cell {_cellname(obj_c, obj_r)} must contain a formula")
    obj_quad = extract_quadratic(obj_ast, var_set, grid)
    obj_is_quadratic = obj_quad.has_quadratic()
    obj_form = LinearForm(dict(obj_quad.linear), obj_quad.constant)
    c_vec = [obj_form.coeffs.get(v, 0.0) for v in decision_vars]
    if not all(math.isfinite(x) for x in c_vec) or not math.isfinite(obj_form.constant):
        raise OptError(
            f"objective cell {_cellname(*objective_cell)} has a non-finite value; "
            "check the cells it references"
        )
    # The objective constant is dropped here: the solver is given only the
    # linear part. We add it back to the reported objective below.

    # Constraints.
    integral = (integer_vars or set()) | (binary_vars or set())
    A: list[list[float]] = []
    sense: list[int] = []
    rhs: list[float] = []
    for c, r in constraint_cells:
        cell = grid.cells[c][r]
        cell_ast = _cell_ast(cell) if cell.type == FORMULA else None
        if cell.type != FORMULA or cell_ast is None:
            raise OptError(f"constraint cell {_cellname(c, r)} must contain a comparison formula")
        coeffs, op_code, rhs_val = extract_constraint(cell_ast, var_set, grid)
        if rhs_override and (c, r) in rhs_override:
            rhs_val = float(rhs_override[(c, r)])
        row = [coeffs.get(v, 0.0) for v in decision_vars]
        # A constraint whose right-hand side reads an error cell (`=SQRT(-1)`,
        # a division by zero) arrives here as NaN. Bounds are already checked
        # for this; the constraint rows were not, and HiGHS rejects the model
        # with a bare "Highs_passLp failed" that names nothing.
        if not math.isfinite(rhs_val) or not all(math.isfinite(x) for x in row):
            raise OptError(
                f"constraint cell {_cellname(c, r)} has a non-finite value; "
                "check the cells it references"
            )
        if isinstance(cell_ast, BinOp) and cell_ast.op in ("<", ">"):
            # An integer row has integer activity, so `< b` is `<= ceil(b) - 1`.
            # A continuous row has no strict LP form and no nearest bound.
            if not all(
                x.is_integer() and (x == 0.0 or decision_vars[j] in integral)
                for j, x in enumerate(row)
            ):
                raise OptError(
                    f"constraint cell {_cellname(c, r)} uses strict '{cell_ast.op}', which "
                    f"LP cannot express; use '{cell_ast.op}=' or make its variables integer "
                    "with integer coefficients"
                )
            rhs_val = math.ceil(rhs_val) - 1.0 if cell_ast.op == "<" else math.floor(rhs_val) + 1.0
        A.append(row)
        sense.append(op_code)
        rhs.append(rhs_val)

    if rhs_override:
        unknown = sorted(set(rhs_override) - set(constraint_cells))
        if unknown:
            raise OptError(
                f"rhs_override names {_cellname(*unknown[0])} which is not a constraint cell"
            )

    # Bounds: default to [0, +inf) for each decision var, matching the
    # "amounts" intuition (no negative production levels).
    inf = float("inf")
    lb = [0.0] * n
    ub = [inf] * n
    if bounds:
        for cell_key, (lo, hi) in bounds.items():
            i = var_index.get(cell_key)
            if i is None:
                raise OptError(
                    f"bounds reference {_cellname(*cell_key)} which is not a decision variable"
                )
            # Validate here rather than letting the C++ bridge reject it.
            # The bridge raises ValueError("lb[j] > ub[j]") -- a column index
            # the user never sees, in an exception type callers of this
            # module do not expect, from a layer they cannot catch
            # meaningfully. Both conditions are reachable from a typed
            # `bounds A1=20:10` or `A1=nan:5`.
            lo_f, hi_f = float(lo), float(hi)
            name = _cellname(*cell_key)
            if math.isnan(lo_f) or math.isnan(hi_f):
                raise OptError(f"bounds for {name} are not numeric")
            if lo_f > hi_f:
                raise OptError(
                    f"bounds for {name} are reversed: lower {lo_f:g} exceeds upper {hi_f:g}"
                )
            lb[i] = lo_f
            ub[i] = hi_f

    # Integer / binary flags. Both must be subsets of decision_vars, and
    # they must be disjoint -- the C++ bridge re-checks for overlap but we
    # surface a clearer message here with the offending cell names.
    int_set = integer_vars or set()
    bin_set = binary_vars or set()
    for cell_key in int_set | bin_set:
        if cell_key not in var_index:
            raise OptError(
                f"integer/binary flag references {_cellname(*cell_key)} "
                "which is not a decision variable"
            )
    overlap = int_set & bin_set
    if overlap:
        c0, r0 = next(iter(overlap))
        raise OptError(f"cell {_cellname(c0, r0)} cannot be both integer and binary")
    int_indices = sorted(var_index[k] for k in int_set)
    bin_indices = sorted(var_index[k] for k in bin_set)

    # A quadratic objective goes to the solver as a Hessian. Cross terms are
    # supported; convexity is the solver's to check, and it rejects an
    # indefinite Hessian rather than returning a plausible wrong answer.
    hessian: list[list[float]] = []
    if obj_is_quadratic:
        if int_indices or bin_indices:
            raise OptError(
                "integer or binary variables cannot be combined with a quadratic objective"
            )
        hessian = obj_quad.hessian(decision_vars)
        check_convexity(hessian, maximize=maximize)
        # A quadratic model's duals do not carry the shadow-price reading the
        # sensitivity report describes. Withhold them, as for MIPs.
        sensitivity = False
        diagnose = False

    # Solve.
    sol = _solve_lp(
        c_vec,
        A,
        sense,
        rhs,
        lb,
        ub,
        maximize=maximize,
        integer_vars=int_indices,
        binary_vars=bin_indices,
        sensitivity=sensitivity,
        hessian=hessian,
    )

    # Add back the constant term that we dropped from the objective vector
    # so the user sees the formula's actual value at the optimum.
    solved_ok = sol.status in (_ext.OPTIMAL, _ext.SUBOPTIMAL)
    objective_total = sol.objective + obj_form.constant if solved_ok else 0.0

    values: dict[CellKey, float] = {}
    if solved_ok:
        # `sol.x` is longer than `decision_vars` when a quadratic relaxation
        # appended auxiliary columns; the trailing entries are not cells.
        for v, x in zip(decision_vars, sol.x[: len(decision_vars)], strict=True):
            values[v] = float(x)
        if obj_is_quadratic:
            # Recompute from the formula rather than trusting the solver's
            # objective: it is cheap, and it keeps the reported number tied to
            # what the user's cell actually says.
            objective_total = evaluate_quadratic(obj_quad, values)

    applied = False
    if apply and values:
        for (c, r), x in values.items():
            # `_ensure_cell`, not `grid.cells[c][r]`: decision cells are
            # allowed to be empty, and empty coordinates hand back a shared
            # placeholder rather than a stored Cell. Writing through that
            # placeholder used to corrupt every empty cell in the process.
            cell = grid._ensure_cell(c, r)
            cell.type = NUM
            cell.val = x
            cell.text = number_text(x)
            cell.ast = None
            cell.ast_text = ""
            cell.err = None
            cell.err_msg = None
        grid.recalc()
        applied = True

    sens: Sensitivity | None = None
    if sensitivity and solved_ok and sol.sensitivity_valid:
        sens = _build_sensitivity(
            decision_vars, constraint_cells, c_vec, A, rhs, sense, sol, values
        )

    runaway: list[CellKey] | None = None
    if diagnose and sol.status == _ext.UNBOUNDED:
        runaway = [
            decision_vars[j]
            for j in _unbounded_variables(
                c_vec, A, sense, rhs, lb, ub, maximize, int_indices, bin_indices
            )
        ]

    conflict: list[CellKey] | None = None
    if diagnose and sol.status == _ext.INFEASIBLE:
        conflict = [
            constraint_cells[i]
            for i in _irreducible_conflict(
                c_vec, A, sense, rhs, lb, ub, maximize, int_indices, bin_indices
            )
        ]

    return SolveResult(
        status=sol.status,
        status_name=_STATUS_NAMES.get(sol.status, f"UNKNOWN({sol.status})"),
        objective=objective_total,
        values=values,
        applied=applied,
        sensitivity=sens,
        conflict=conflict,
        unbounded=runaway,
        quadratic=obj_is_quadratic,
    )

sweep

sweep(
    grid: Grid,
    objective_cell: CellKey,
    decision_vars: list[CellKey],
    constraint_cells: list[CellKey],
    *,
    constraint: CellKey,
    lo: float,
    hi: float,
    steps: int = 10,
    maximize: bool = True,
    bounds: dict[CellKey, tuple[float, float]]
    | None = None,
    integer_vars: set[CellKey] | None = None,
    binary_vars: set[CellKey] | None = None,
) -> list[SweepPoint]

Re-solve the model across a range of right-hand sides for one constraint.

A shadow price answers "what is the next unit worth". It is valid only inside its ranging interval, so it cannot answer "how much more should I buy" -- past the interval edge the basis changes and the marginal value drops. Sweeping re-solves at each point and reports where that happens.

The sheet is never modified: each point substitutes the right-hand side via solve(rhs_override=...) with apply=False. The constraint's coefficients still come from its formula; only the constant moves.

steps is the number of intervals, so the result has steps + 1 points spanning lo..hi inclusive. Points where the model becomes infeasible or unbounded are included with their status rather than dropped -- discovering that a resource level is unattainable is a real answer to the question being asked.

Source code in src/gridcalc/opt.py

def sweep(
    grid: Grid,
    objective_cell: CellKey,
    decision_vars: list[CellKey],
    constraint_cells: list[CellKey],
    *,
    constraint: CellKey,
    lo: float,
    hi: float,
    steps: int = 10,
    maximize: bool = True,
    bounds: dict[CellKey, tuple[float, float]] | None = None,
    integer_vars: set[CellKey] | None = None,
    binary_vars: set[CellKey] | None = None,
) -> list[SweepPoint]:
    """Re-solve the model across a range of right-hand sides for one constraint.

    A shadow price answers "what is the next unit worth". It is valid only
    inside its ranging interval, so it cannot answer "how much more should I
    buy" -- past the interval edge the basis changes and the marginal value
    drops. Sweeping re-solves at each point and reports where that happens.

    The sheet is never modified: each point substitutes the right-hand side
    via ``solve(rhs_override=...)`` with ``apply=False``. The constraint's
    coefficients still come from its formula; only the constant moves.

    ``steps`` is the number of intervals, so the result has ``steps + 1``
    points spanning ``lo``..``hi`` inclusive. Points where the model becomes
    infeasible or unbounded are included with their status rather than
    dropped -- discovering that a resource level is unattainable is a real
    answer to the question being asked.
    """
    if steps < 1:
        raise OptError("sweep needs at least 1 step")
    if hi < lo:
        raise OptError(f"sweep range is reversed: {lo:g} to {hi:g}")
    if constraint not in constraint_cells:
        raise OptError(f"{_cellname(*constraint)} is not one of the constraint cells")

    points: list[SweepPoint] = []
    prev_obj: float | None = None
    prev_price: float | None = None

    for k in range(steps + 1):
        value = lo if steps == 0 else lo + (hi - lo) * k / steps
        result = solve(
            grid,
            objective_cell,
            decision_vars,
            constraint_cells,
            maximize=maximize,
            bounds=bounds,
            integer_vars=integer_vars,
            binary_vars=binary_vars,
            apply=False,
            sensitivity=True,
            rhs_override={constraint: value},
        )
        solved = result.status_name in ("OPTIMAL", "SUBOPTIMAL")

        price: float | None = None
        if result.sensitivity is not None:
            for c in result.sensitivity.constraints:
                if c.cell == constraint:
                    price = c.shadow_price
                    break

        delta = result.objective - prev_obj if (solved and prev_obj is not None) else None
        # Only call it a breakpoint when both prices are known; a None on
        # either side means "not comparable", not "changed".
        changed = price is not None and prev_price is not None and abs(price - prev_price) > 1e-9

        points.append(
            SweepPoint(
                rhs=value,
                status_name=result.status_name,
                objective=result.objective if solved else float("nan"),
                shadow_price=price,
                delta=delta,
                breakpoint=changed,
            )
        )
        if solved:
            prev_obj = result.objective
        prev_price = price

    return points

extract_linear

extract_linear(
    node: Node, decision_vars: set[CellKey], grid: Grid
) -> LinearForm

Reduce node to a LinearForm over decision_vars.

Cells in decision_vars contribute coefficients; all other cells are looked up in grid and folded into the constant term.

Source code in src/gridcalc/opt.py

def extract_linear(node: Node, decision_vars: set[CellKey], grid: Grid) -> LinearForm:
    """Reduce ``node`` to a LinearForm over ``decision_vars``.

    Cells in ``decision_vars`` contribute coefficients; all other cells are
    looked up in ``grid`` and folded into the constant term.
    """
    if isinstance(node, Number):
        return LinearForm({}, float(node.value))

    if isinstance(node, CellRef):
        _check_sheet(node.sheet, _active_sheet_name(grid))
        key: CellKey = (node.col, node.row)
        if key in decision_vars:
            return LinearForm({key: 1.0}, 0.0)
        helper = _helper_ast(grid, key, decision_vars)
        if helper is not None:
            return _inline(extract_linear, helper, key, decision_vars, grid)
        return LinearForm({}, _cell_value(grid, node.col, node.row))

    if isinstance(node, UnaryOp):
        inner = extract_linear(node.operand, decision_vars, grid)
        if node.op == "+":
            return inner
        if node.op == "-":
            return inner.neg()
        raise NotLinear(f"unsupported unary operator '{node.op}'")

    if isinstance(node, Percent):
        return extract_linear(node.operand, decision_vars, grid).scale(0.01)

    if isinstance(node, BinOp):
        if node.op == "+":
            return extract_linear(node.left, decision_vars, grid).add(
                extract_linear(node.right, decision_vars, grid)
            )
        if node.op == "-":
            return extract_linear(node.left, decision_vars, grid).sub(
                extract_linear(node.right, decision_vars, grid)
            )
        if node.op == "*":
            lhs = extract_linear(node.left, decision_vars, grid)
            rhs = extract_linear(node.right, decision_vars, grid)
            if lhs.is_constant:
                return rhs.scale(lhs.constant)
            if rhs.is_constant:
                return lhs.scale(rhs.constant)
            raise NotLinear("product of two decision-variable expressions is nonlinear")
        if node.op == "/":
            lhs = extract_linear(node.left, decision_vars, grid)
            rhs = extract_linear(node.right, decision_vars, grid)
            if not rhs.is_constant:
                raise NotLinear("division by a decision-variable expression is nonlinear")
            if rhs.constant == 0.0:
                raise NotLinear("division by zero in linear expression")
            return lhs.scale(1.0 / rhs.constant)
        # ^, &, comparisons, etc. -- not allowed inside an expression body
        raise NotLinear(f"unsupported operator '{node.op}' in linear expression")

    if isinstance(node, Call):
        if node.name.upper() == "SUM":
            total = LinearForm()
            for arg in node.args:
                total = total.add(_sum_arg(arg, decision_vars, grid))
            return total
        raise NotLinear(f"function '{node.name}' is not allowed in a linear expression")

    if isinstance(node, (Bool, String, ErrorLit, RangeRef, Name, PyCall)):
        raise NotLinear(f"{type(node).__name__} is not allowed in a linear expression")

    raise NotLinear(f"unhandled AST node: {type(node).__name__}")

extract_quadratic

extract_quadratic(
    node: Node, decision_vars: set[CellKey], grid: Grid
) -> QuadForm

Reduce node to a degree-<=2 form over decision_vars.

The quadratic counterpart of :func:extract_linear, and deliberately a separate walker: constraints must stay linear, so widening the shared one would have quietly admitted quadratic constraints the solver cannot take.

Source code in src/gridcalc/opt.py

def extract_quadratic(node: Node, decision_vars: set[CellKey], grid: Grid) -> QuadForm:
    """Reduce ``node`` to a degree-<=2 form over ``decision_vars``.

    The quadratic counterpart of :func:`extract_linear`, and deliberately a
    separate walker: constraints must stay linear, so widening the shared one
    would have quietly admitted quadratic constraints the solver cannot take.
    """
    if isinstance(node, Number):
        return QuadForm({}, {}, float(node.value))

    if isinstance(node, CellRef):
        _check_sheet(node.sheet, _active_sheet_name(grid))
        key: CellKey = (node.col, node.row)
        if key in decision_vars:
            return QuadForm({}, {key: 1.0}, 0.0)
        helper = _helper_ast(grid, key, decision_vars)
        if helper is not None:
            return _inline(extract_quadratic, helper, key, decision_vars, grid)
        return QuadForm({}, {}, _cell_value(grid, node.col, node.row))

    if isinstance(node, UnaryOp):
        inner = extract_quadratic(node.operand, decision_vars, grid)
        if node.op == "+":
            return inner
        if node.op == "-":
            return inner.neg()
        raise NotQuadratic(f"unsupported unary operator '{node.op}'")

    if isinstance(node, Percent):
        return extract_quadratic(node.operand, decision_vars, grid).scale(0.01)

    if isinstance(node, BinOp):
        if node.op == "+":
            return extract_quadratic(node.left, decision_vars, grid).add(
                extract_quadratic(node.right, decision_vars, grid)
            )
        if node.op == "-":
            return extract_quadratic(node.left, decision_vars, grid).sub(
                extract_quadratic(node.right, decision_vars, grid)
            )
        if node.op == "*":
            return extract_quadratic(node.left, decision_vars, grid).mul(
                extract_quadratic(node.right, decision_vars, grid)
            )
        if node.op == "^":
            base = extract_quadratic(node.left, decision_vars, grid)
            power = extract_quadratic(node.right, decision_vars, grid)
            if not power.is_constant:
                raise NotQuadratic("exponent must be a constant")
            p = power.constant
            if p == 0.0:
                return QuadForm({}, {}, 1.0)
            if p == 1.0:
                return base
            if p == 2.0:
                return base.mul(base)
            raise NotQuadratic(f"exponent {p:g} is not supported (only 0, 1, 2)")
        if node.op == "/":
            lhs = extract_quadratic(node.left, decision_vars, grid)
            rhs = extract_quadratic(node.right, decision_vars, grid)
            if not rhs.is_constant:
                raise NotQuadratic("division by a decision-variable expression")
            if rhs.constant == 0.0:
                raise NotQuadratic("division by zero")
            return lhs.scale(1.0 / rhs.constant)
        raise NotQuadratic(f"unsupported operator '{node.op}'")

    if isinstance(node, Call):
        if node.name.upper() == "SUM":
            total = QuadForm()
            for arg in node.args:
                lin = _sum_arg(arg, decision_vars, grid)
                total = total.add(QuadForm({}, dict(lin.coeffs), lin.constant))
            return total
        raise NotQuadratic(f"function '{node.name}' is not allowed")

    raise NotQuadratic(f"{type(node).__name__} is not allowed in an objective")

extract_constraint

extract_constraint(
    node: Node, decision_vars: set[CellKey], grid: Grid
) -> tuple[dict[CellKey, float], int, float]

Reduce a comparison-rooted formula to (coeffs, sense, rhs) form.

Both sides are walked as linear forms; variables move to the left and constants to the right, so the LP sees a single row a^T x OP b.

Source code in src/gridcalc/opt.py

def extract_constraint(
    node: Node,
    decision_vars: set[CellKey],
    grid: Grid,
) -> tuple[dict[CellKey, float], int, float]:
    """Reduce a comparison-rooted formula to (coeffs, sense, rhs) form.

    Both sides are walked as linear forms; variables move to the left and
    constants to the right, so the LP sees a single row ``a^T x OP b``.
    """
    if not isinstance(node, BinOp) or node.op not in _SENSE:
        if isinstance(node, BinOp) and node.op == "<>":
            raise OptError("'<>' is not a valid LP constraint operator")
        raise OptError("constraint formula must be a comparison (<=, >=, =, <, >)")
    lhs = extract_linear(node.left, decision_vars, grid)
    rhs = extract_linear(node.right, decision_vars, grid)
    diff = lhs.sub(rhs)  # coeffs * x + (lhs.const - rhs.const) OP 0
    rhs_value = -diff.constant  # move constant to RHS
    return diff.coeffs, _SENSE[node.op], rhs_value

parse_cells

parse_cells(spec: str) -> list[tuple[int, int]]

Expand a cell-list spec like A1:B3 or A1,A2,B5 into (col,row)s.

Returns the cells in row-major order within each range and in spec order across comma-separated parts. Duplicate-detection is the caller's job. Raises ValueError on a malformed ref.

Source code in src/gridcalc/opt.py

def parse_cells(spec: str) -> list[tuple[int, int]]:
    """Expand a cell-list spec like ``A1:B3`` or ``A1,A2,B5`` into (col,row)s.

    Returns the cells in row-major order within each range and in spec order
    across comma-separated parts. Duplicate-detection is the caller's job.
    Raises ``ValueError`` on a malformed ref.
    """
    out: list[tuple[int, int]] = []
    for part in spec.split(","):
        part = part.strip()
        if not part:
            continue
        if ":" in part:
            a_str, b_str = part.split(":", 1)
            a = ref(a_str.strip())
            b = ref(b_str.strip())
            if not a or not b:
                raise ValueError(f"bad cell range: {part}")
            _, c1, r1 = a
            _, c2, r2 = b
            c1, c2 = sorted((c1, c2))
            r1, r2 = sorted((r1, r2))
            for c in range(c1, c2 + 1):
                for r in range(r1, r2 + 1):
                    out.append((c, r))
        else:
            m = ref(part)
            if not m:
                raise ValueError(f"bad cell ref: {part}")
            _, c, r = m
            out.append((c, r))
    return out

parse_bounds

parse_bounds(
    spec: str,
) -> dict[tuple[int, int], tuple[float, float]]

Parse A1=lo:hi,B2=lo:hi into a bounds dict. Raises ValueError.

Source code in src/gridcalc/opt.py

def parse_bounds(spec: str) -> dict[tuple[int, int], tuple[float, float]]:
    """Parse ``A1=lo:hi,B2=lo:hi`` into a bounds dict. Raises ``ValueError``."""
    out: dict[tuple[int, int], tuple[float, float]] = {}
    for part in spec.split(","):
        part = part.strip()
        if not part:
            continue
        if "=" not in part:
            raise ValueError(f"bounds entry missing '=': {part}")
        cellref_str, range_str = part.split("=", 1)
        m = ref(cellref_str.strip())
        if not m:
            raise ValueError(f"bad cell ref in bounds: {cellref_str}")
        _, c, r = m
        if ":" not in range_str:
            raise ValueError(f"bounds range needs 'lo:hi': {range_str}")
        lo_s, hi_s = range_str.split(":", 1)
        try:
            out[(c, r)] = (
                _parse_bound_value(lo_s, positive=False),
                _parse_bound_value(hi_s, positive=True),
            )
        except ValueError as e:
            raise ValueError(f"bad bounds for {cellref_str.strip()}: {e}") from e
    return out

cells_to_spec

cells_to_spec(cells: list[tuple[int, int]]) -> str

Render a cell list as a comma-separated spec string -- the inverse of :func:parse_cells.

Inferred models are stored as specs, exactly like typed ones, so they round-trip through the workbook JSON and can be re-run after reopening. A contiguous single-column or single-row run collapses to range syntax so the saved model stays readable.

Source code in src/gridcalc/opt.py

def cells_to_spec(cells: list[tuple[int, int]]) -> str:
    """Render a cell list as a comma-separated spec string -- the inverse of
    :func:`parse_cells`.

    Inferred models are stored as specs, exactly like typed ones, so they
    round-trip through the workbook JSON and can be re-run after reopening. A
    contiguous single-column or single-row run collapses to range syntax so the
    saved model stays readable.
    """
    if not cells:
        return ""
    cols = {c for c, _ in cells}
    rows = {r for _, r in cells}
    if len(cols) == 1:
        rs = sorted(rows)
        if rs == list(range(rs[0], rs[-1] + 1)) and len(rs) > 1:
            c = next(iter(cols))
            return f"{_cellname(c, rs[0])}:{_cellname(c, rs[-1])}"
    if len(rows) == 1:
        cs = sorted(cols)
        if cs == list(range(cs[0], cs[-1] + 1)) and len(cs) > 1:
            r = next(iter(rows))
            return f"{_cellname(cs[0], r)}:{_cellname(cs[-1], r)}"
    return ",".join(_cellname(c, r) for c, r in cells)