From e594ea09ec6d92e09e7c93e9df9f2a677a01e1d3 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 27 Sep 2026 20:08:00 +0000 Subject: [PATCH 1/2] feat: a spec reports how its optimum moves with a parameter, which for a fixed decision is a Benders cut Spec.sensitivity(*names) adds a reported expression _sensitivity for each named parameter: the objective's own rate in it, plus dual(c) times the rate of rhs - lhs for each row that reads it. named.md now defines dual(c) as the rate at which the optimal objective rises with the right side of c. That is what specsolve returns. It gives an equality a sign, and it makes the rate the same expression under either sense. The old wording made a <= dual non-negative under minimize, which specsolve does not return. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01CvcU9JRnsinc9e9syvVcHS --- docs/reference/language/named.md | 15 +- src/mathspec/sensitivity.py | 236 +++++++++++++++++++++++++++++++ src/mathspec/spec.py | 35 +++++ tests/test_public_surface.py | 5 +- tests/test_sensitivity.py | 117 +++++++++++++++ 5 files changed, 403 insertions(+), 5 deletions(-) create mode 100644 src/mathspec/sensitivity.py create mode 100644 tests/test_sensitivity.py diff --git a/docs/reference/language/named.md b/docs/reference/language/named.md index 3df1763d..1e44fced 100644 --- a/docs/reference/language/named.md +++ b/docs/reference/language/named.md @@ -151,12 +151,21 @@ Constraint 'd': a dual exists only after a solve; the math cannot read one — keep the entry that carries it out of constraints, the objective, bounds and where. ``` -`dual(c)` is the rate at which the optimal objective improves as `c` is relaxed -in the direction its comparator points, under the spec's own `minimize` or -`maximize`. +`dual(c)` is the rate at which the optimal objective rises as the right side +of `c` rises. For `lhs <= rhs`, read `lhs <= rhs + b`: the dual is the rate in +`b` at `b = 0`. The same holds for `>=` and `==`, and under `minimize` and +`maximize` alike. So under `minimize`, the dual of a binding `<=` row is at +most zero, and the dual of a binding `>=` row is at least zero. A row that `c`'s `where:` deletes has no dual. +[`spec.sensitivity('q')`](../api.md#mathspec.Spec.sensitivity) writes this sum +as a reported expression, `q_sensitivity`: the rate at which the optimal +objective rises per unit of the parameter `q`. It is the objective's own rate +in `q`, plus `dual(c)` times the rate of `rhs - lhs` in `q` for each row `c` +that reads `q`. Of a variable that [`spec.fix`](../api.md#mathspec.Spec.fix) +made a parameter, it is the slope of a Benders cut. + ## `macros` A macro is a template that takes arguments and is substituted into an expression diff --git a/src/mathspec/sensitivity.py b/src/mathspec/sensitivity.py new file mode 100644 index 00000000..9d58716f --- /dev/null +++ b/src/mathspec/sensitivity.py @@ -0,0 +1,236 @@ +# SPDX-FileCopyrightText: mathspec Contributors +# +# SPDX-License-Identifier: MIT + +"""How the optimal objective moves with a parameter, written as a reported expression over the duals. + +At an optimum the objective's rate of change in a parameter ``q`` is what ``q`` +moves directly, plus, for every row that reads ``q``, the row's dual times how +far ``q`` raises the row's right side against its left. ``dual(c)`` is the rate +the optimal objective rises with the right side of ``c`` +(``docs/reference/language/named.md``), so a row ``lhs <= rhs`` contributes +``dual(c) * d(rhs - lhs)/dq`` whatever its comparator and whatever the sense. + +How far ``q`` moves a row is read off the row's tree from the root down, the +seed multiplied by each factor it passes and summed back over every dimension +the parameter does not carry. The walk spells what it multiplies by, so it +stops at an operator it cannot write back: a shift, a relation, a window, a +case. +""" + +from __future__ import annotations + +from dataclasses import dataclass, replace +from typing import TYPE_CHECKING + +from mathspec.dimensions import dims_of +from mathspec.errors import LanguageError, SchemaError, did_you_mean +from mathspec.program import ( + Add, + Constant, + Divide, + Dual, + Multiply, + Named, + Negate, + Parameter, + Power, + Sum, + Variable, + parameters_of, +) +from mathspec.sos import section +from mathspec.validation import to_spec + +if TYPE_CHECKING: + from mathspec.program import Expression + from mathspec.spec import Spec + + +def sensitivity(spec: Spec, names: tuple[str, ...]) -> Spec: + """*spec* with ``_sensitivity`` reported for each of *names*; see [`Spec.sensitivity`][mathspec.Spec.sensitivity].""" + program = spec.program + unknown = [ + f"sensitivity: '{name}' is not a parameter of this spec. " + + did_you_mean(name, program.parameters, label='Parameters') + for name in dict.fromkeys(names) + if name not in program.parameters + ] + if unknown: + raise SchemaError('\n'.join(unknown)) + raw = spec.to_dict() + for name in dict.fromkeys(names): + entry = f'{name}_sensitivity' + if entry in spec.expressions: + raise LanguageError(f"sensitivity: '{entry}' is already an expression of this spec.") + section(raw, 'expressions')[entry] = { + 'expression': _derivative(spec, name), + 'description': f"how fast the optimal objective rises per unit of '{name}', read from the duals", + } + return to_spec(raw) + + +def _derivative(spec: Spec, name: str) -> str: + """The rate of the optimal objective in parameter *name*, as an expression over the duals. + + Raises: + LanguageError: No objective; a bound or an operator the walk cannot + write back reads *name*; or nothing reads it. + """ + program = spec.program + if program.objective is None: + raise LanguageError(f"sensitivity of '{name}': the spec declares no objective to differentiate.") + if bounded := [v for v, d in program.variables.items() if name in parameters_of(*_bounds(d.lower, d.upper))]: + raise LanguageError( + f"sensitivity of '{name}': it bounds variable '{bounded[0]}', and a bound's price is a reduced cost, " + f'which no expression reads. Write the bound as a constraint.' + ) + terms = [ + term + for term, _ in _pulled(program.objective.expression, _Seed(1, None, frozenset()), name, spec, 'the objective') + ] + for label, row in program.constraints.items(): + if name not in parameters_of(row.lhs, row.rhs): + continue + seed = _Seed(1, f'dual({label})', frozenset(row.dims)) + found = [ + *_pulled(row.rhs, seed, name, spec, f"constraint '{label}'"), + *_pulled(row.lhs, seed.negated(), name, spec, f"constraint '{label}'"), + ] + if row.where is not None and any(pointwise for _, pointwise in found): + raise LanguageError( + f"sensitivity of '{name}': constraint '{label}' has a where:, and where it deletes a row the " + f'dual is absent and would take the whole rate with it. Read it where no mask deletes a row.' + ) + terms += [term for term, _ in found] + if not terms: + raise LanguageError(f"sensitivity of '{name}': neither the objective nor a constraint reads it.") + return _joined(terms) + + +def _joined(terms: list[str]) -> str: + """The terms as one sum, a negative term written as a subtraction.""" + text = terms[0] + for term in terms[1:]: + text += f' - {term[1:]}' if term.startswith('-') else f' + {term}' + return text + + +@dataclass(frozen=True) +class _Seed: + """The rate of a row in the node the walk stands at: a sign, the product so far, and the dims it carries. + + ``text`` is ``None`` while the product is still ``1``. + """ + + sign: int + text: str | None + carried: frozenset[str] + + def negated(self) -> _Seed: + return replace(self, sign=-self.sign) + + def times(self, factor: str, dims: frozenset[str]) -> _Seed: + return _Seed(self.sign, factor if self.text is None else f'{self.text} * {factor}', self.carried | dims) + + def over(self, divisor: str, dims: frozenset[str]) -> _Seed: + return _Seed(self.sign, f'{self.text or 1} / {divisor}', self.carried | dims) + + def spelled(self, summed: list[str]) -> str: + text = self.text or '1' + for dim in summed: + text = f'sum({text}, over={dim})' + return f'-{text}' if self.sign < 0 else text + + +def _bounds(*bounds: Expression | None) -> tuple[Expression, ...]: + return tuple(bound for bound in bounds if bound is not None) + + +def _pulled(node: Expression, seed: _Seed, name: str, spec: Spec, label: str) -> list[tuple[str, bool]]: + """What each read of *name* under *node* contributes, summed back to *name*'s dims, and whether it was pointwise.""" + if name not in parameters_of(node): + return [] + if isinstance(node, Parameter): + extra = sorted(seed.carried - frozenset(spec.program.parameters[name].dims)) + return [(seed.spelled(extra), not extra)] + if isinstance(node, Negate): + return _pulled(node.operand, seed.negated(), name, spec, label) + if isinstance(node, Add): + return [*_pulled(node.left, seed, name, spec, label), *_pulled(node.right, seed, name, spec, label)] + if isinstance(node, Multiply): + return [ + *_pulled(node.left, _times(seed, node.right, spec, label), name, spec, label), + *_pulled(node.right, _times(seed, node.left, spec, label), name, spec, label), + ] + if isinstance(node, Divide): + if name in parameters_of(node.divisor): + raise LanguageError( + f"sensitivity of '{name}': {label} divides by it, which this derivation does not carry." + ) + divided = seed.over(_factor(node.divisor, label), dims_of(node.divisor, spec, label)) + return _pulled(node.numerator, divided, name, spec, label) + if isinstance(node, Sum): + if clash := seed.carried & frozenset(node.over): + raise LanguageError( + f"sensitivity of '{name}': {label} sums over {min(clash)} inside a factor that carries it, " + f'and the one name cannot be both.' + ) + return _pulled(node.operand, seed, name, spec, label) + if isinstance(node, Named): + return _pulled(node.body, seed, name, spec, label) + raise LanguageError( + f"sensitivity of '{name}': {label} reads it through {type(node).__name__}, which this derivation does " + f'not carry yet: only arithmetic and sum(over=) lie between a row and the parameter.' + ) + + +def _times(seed: _Seed, factor: Expression, spec: Spec, label: str) -> _Seed: + """The seed multiplied by the other side of a product, a negated or negative factor folded into the sign.""" + sign = 1 + while isinstance(factor, Negate): + sign, factor = -sign, factor.operand + if isinstance(factor, Constant) and factor.value < 0: + sign, factor = -sign, Constant(-factor.value) + scaled = seed if sign == 1 else seed.negated() + if isinstance(factor, Constant) and factor.value == 1: + return scaled + return scaled.times(_factor(factor, label), dims_of(factor, spec, label)) + + +def _factor(node: Expression, label: str) -> str: + """*node* spelled as one factor of a product: bare where it is a name or a number, bracketed otherwise.""" + text = _spelled(node, label) + return text if isinstance(node, Parameter | Variable | Named | Constant | Dual | Sum) else f'({text})' + + +def _spelled(node: Expression, label: str) -> str: + """*node* written back as an expression string, every operation bracketed. + + Raises: + LanguageError: An operator other than arithmetic and ``sum(over=)``. + """ + if isinstance(node, Constant): + return repr(node.value) if node.value >= 0 else f'({node.value!r})' + if isinstance(node, Parameter | Variable | Named): + return node.name + if isinstance(node, Dual): + return f'dual({node.constraint})' + if isinstance(node, Negate): + return f'-({_spelled(node.operand, label)})' + if isinstance(node, Add | Multiply | Divide | Power): + left, right = ( + (node.left, node.right) + if isinstance(node, Add | Multiply) + else ((node.numerator, node.divisor) if isinstance(node, Divide) else (node.base, node.exponent)) + ) + symbol = {Add: '+', Multiply: '*', Divide: '/', Power: '**'}[type(node)] + return f'({_spelled(left, label)}) {symbol} ({_spelled(right, label)})' + if isinstance(node, Sum): + text = _spelled(node.operand, label) + for dim in node.over: + text = f'sum({text}, over={dim})' + return text + raise LanguageError( + f'sensitivity: {label} multiplies by a {type(node).__name__}, which this derivation cannot write back yet.' + ) diff --git a/src/mathspec/spec.py b/src/mathspec/spec.py index 1397c99c..fe7f30f6 100644 --- a/src/mathspec/spec.py +++ b/src/mathspec/spec.py @@ -897,6 +897,41 @@ def fix(self, *names: str) -> Spec: return fix(self, names) + def sensitivity(self, *names: str) -> Spec: + """This spec reporting how fast the optimal objective moves with each named parameter. + + Each name gets a reported expression, ``_sensitivity``, over the + parameter's dims: the objective's own rate in it, plus each row that + reads it, its dual times how far the parameter raises the row's right + side against its left. ``dual(c)`` is the rate the optimal objective + rises with the right side of ``c``, so the entry is the rise of the + optimal objective per unit of the parameter, under either sense and + for every comparator. Of a variable [`fix`][] made a parameter, it is + the slope of a Benders cut; of a load, the price of serving one more + unit. + + It holds where the solve has duals and they are unique: a linear + program, or a mixed-integer one's relaxation, away from degeneracy. A + ``where:`` or an assumption reading the parameter is a step the rate + does not see. + + Args: + names: The parameters to differentiate by. + + Returns: + A new spec with one reported expression per name. + + Raises: + SchemaError: A name that is not a parameter. + LanguageError: No objective; a bound reads the parameter, whose + price is a reduced cost; a row with a ``where:`` reads it + pointwise; or it is read through an operator other than + arithmetic and ``sum(over=)``. + """ + from mathspec.sensitivity import sensitivity + + return sensitivity(self, names) + @model_validator(mode='after') def _names_are_names(self) -> Spec: """Every declaration is keyed by something an expression could write. diff --git a/tests/test_public_surface.py b/tests/test_public_surface.py index 95aff005..c4703e9f 100644 --- a/tests/test_public_surface.py +++ b/tests/test_public_surface.py @@ -39,9 +39,10 @@ #: What `Spec` promises beyond the sections a file declares: the two ways back #: out, the verb that writes a formulation out, the verb that turns a decision -#: into a supplied number, and the program the file means. +#: into a supplied number, the verb that reports how the optimum moves with a +#: number, and the program the file means. #: A `model_`-prefixed name is pydantic's, not a contract this project keeps. -SPEC_SURFACE = frozenset({'to_dict', 'to_yaml', 'expand', 'fix', 'program'}) +SPEC_SURFACE = frozenset({'to_dict', 'to_yaml', 'expand', 'fix', 'sensitivity', 'program'}) #: The modules whose `__all__` a consumer imports from. MODULES = [ diff --git a/tests/test_sensitivity.py b/tests/test_sensitivity.py new file mode 100644 index 00000000..ec0fb82b --- /dev/null +++ b/tests/test_sensitivity.py @@ -0,0 +1,117 @@ +# SPDX-FileCopyrightText: mathspec Contributors +# +# SPDX-License-Identifier: MIT + +"""What `Spec.sensitivity` writes, and where it refuses to write anything. + +The rate is only as right as its sign, and the sign is what `dual(c)` means: +the rise of the optimal objective with the right side of `c`. The numbers +behind that reading were checked against finite differences on HiGHS; the +rows here pin the expression the derivation writes, which no solver is needed +to read. +""" + +from __future__ import annotations + +import pytest + +from mathspec import LanguageError, SchemaError +from tests.fixtures import override, schema_of +from tests.test_fixing import MONOLITH + +SUBPROBLEM = schema_of(MONOLITH).fix('cap').to_dict() + + +def _written(model, *names, **patch) -> dict[str, str]: + spec = schema_of(model, **patch).sensitivity(*names) + return {name: spec.expressions[f'{name}_sensitivity'].expression for name in names} + + +def test_the_cut_of_a_benders_subproblem(): + """The objective reads `cap` directly, and `capacity` reads it on its right side.""" + assert _written(SUBPROBLEM, 'cap', 'load') == { + 'cap': 'invest + sum(dual(capacity) * avail, over=snapshot)', + 'load': 'dual(balance)', + } + + +def test_the_rate_is_over_the_parameters_own_dims(): + program = schema_of(SUBPROBLEM).sensitivity('cap', 'load').program + assert program.expressions['cap_sensitivity'].dims == ('generator',), 'summed back over snapshot' + assert program.expressions['load_sensitivity'].dims == ('snapshot',) + assert not program.expressions['cap_sensitivity'].in_math, 'a reported expression: no solver sees it' + + +@pytest.mark.parametrize( + ('capacity', 'written'), + [ + pytest.param('p <= cap * avail', 'invest + sum(dual(capacity) * avail, over=snapshot)', id='on-the-right'), + pytest.param('cap * avail >= p', 'invest - sum(dual(capacity) * avail, over=snapshot)', id='on-the-left'), + pytest.param('p <= 2 * cap * avail', 'invest + sum(dual(capacity) * avail * 2.0, over=snapshot)', id='scaled'), + pytest.param('p <= -(-cap) * avail', 'invest + sum(dual(capacity) * avail, over=snapshot)', id='negated-twice'), + pytest.param('p == cap * avail', 'invest + sum(dual(capacity) * avail, over=snapshot)', id='an-equality'), + pytest.param( + 'p - cap * avail / 2 <= 0', 'invest + sum(dual(capacity) / 2.0 * avail, over=snapshot)', id='divided' + ), + ], +) +def test_the_rows_sides_decide_the_sign_and_its_comparator_does_not(capacity, written): + assert _written(SUBPROBLEM, 'cap', **{'constraints.capacity.expression': capacity})['cap'] == written + + +def test_maximize_takes_the_same_expression(): + """The dual is the rise of the objective with the right side under either sense, so nothing flips.""" + raw = override(SUBPROBLEM, objective={'sense': 'maximize', 'expression': '-sum(cap * invest) - sum(p * cost)'}) + assert _written(raw, 'cap')['cap'] == '-invest + sum(dual(capacity) * avail, over=snapshot)' + + +def test_a_masked_row_summed_back_is_one_summand_fewer(): + written = _written(SUBPROBLEM, 'cap', **{'constraints.capacity.where': 'avail > 0'}) + assert written['cap'] == 'invest + sum(dual(capacity) * avail, over=snapshot)' + + +@pytest.mark.parametrize( + ('names', 'patch', 'says'), + [ + pytest.param(('capp',), {}, "'capp' is not a parameter of this spec. Did you mean 'cap'", id='a-near-miss'), + pytest.param(('p',), {}, "'p' is not a parameter", id='a-variable'), + ], +) +def test_a_name_that_is_not_a_parameter_is_refused(names, patch, says): + with pytest.raises(SchemaError, match=says): + _written(SUBPROBLEM, *names, **patch) + + +@pytest.mark.parametrize( + ('name', 'patch', 'says'), + [ + pytest.param( + 'cap', + {'constraints.capacity.expression': 'p <= shift(cap, along=generator, offset=1, edge=0) * avail'}, + 'through Translate', + id='through-a-shift', + ), + pytest.param('cap', {'constraints.capacity.expression': 'p / cap <= avail'}, 'divides by it', id='a-divisor'), + pytest.param( + 'cap', + {'constraints.capacity.expression': 'p <= avail * sum(cap, over=generator)'}, + 'sums over generator inside a factor that carries it', + id='a-sum-over-a-dim-the-factor-outside-carries', + ), + pytest.param( + 'load', {'constraints.balance.where': 'load > 0'}, 'has a where:', id='a-masked-row-read-pointwise' + ), + pytest.param('avail', {'variables.p.bounds': {'lower': 0, 'upper': 'avail'}}, 'bounds variable', id='a-bound'), + pytest.param( + 'invest', {'objective.expression': 'sum(p * cost)'}, 'neither the objective nor a constraint', id='unread' + ), + ], +) +def test_what_the_derivation_cannot_carry_is_refused(name, patch, says): + with pytest.raises(LanguageError, match=says): + _written(SUBPROBLEM, name, **patch) + + +def test_an_entry_of_that_name_is_refused(): + with pytest.raises(LanguageError, match="'cap_sensitivity' is already an expression"): + _written(SUBPROBLEM, 'cap', expressions={'cap_sensitivity': 'sum(avail)'}) From 07dfd257f8fb6c6604657579ae8843c1ae139e5e Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 27 Sep 2026 20:08:30 +0000 Subject: [PATCH 2/2] docs: the changelog names #749 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01CvcU9JRnsinc9e9syvVcHS --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 985c66f3..953e8549 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ it releases that version ([RELEASING.md](https://github.com/energy-models/mathsp ## Upcoming version +- feat: a spec reports how its optimum moves with a parameter, which for a fixed decision is a Benders cut ([#749](https://github.com/energy-models/mathspec/pull/749)) - feat: a decided variable becomes a supplied number, so a subproblem is a call on the whole model ([#748](https://github.com/energy-models/mathspec/pull/748)) - feat: advice names the assumption that proves a quadratic objective or row convex ([#747](https://github.com/energy-models/mathspec/pull/747)) - feat: a program names its problem class and proves from the file whether its quadratic part is convex ([#746](https://github.com/energy-models/mathspec/pull/746))