diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index f6c94a06e..acf5c8f41 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -116,6 +116,10 @@ jobs: run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_xhat_from_file.py -v + - name: Test xhat feasibility cuts + run: | + coverage run $COV_ARGS -m pytest mpisppy/tests/test_xhat_feasibility_cuts.py -v + - name: Test feasible_xhat run: | coverage run $COV_ARGS -m pytest mpisppy/tests/test_feasible_xhat.py -v diff --git a/doc/designs/xhat_feasibility_cuts_design.md b/doc/designs/xhat_feasibility_cuts_design.md new file mode 100644 index 000000000..7e53905f8 --- /dev/null +++ b/doc/designs/xhat_feasibility_cuts_design.md @@ -0,0 +1,404 @@ +# Design: Optional Feasibility Cuts from Xhatters + +**Status:** Draft — up for discussion. Nothing is implemented yet. +**Addresses:** [issue #601](https://github.com/Pyomo/mpi-sppy/issues/601) +**Author:** dlw (captured with Claude Code assistance) +**Last updated:** 2026-04-23 + +## Motivation + +For two-stage problems **without complete recourse**, an xhatter (an xhat +spoke such as `xhatlooper`, `xhatshufflelooper`, `xhatspecific`) can +propose a first-stage candidate `xhat` that is infeasible in one or more +of the rank-local scenarios. Today the xhatter just discards `xhat` and +moves on — but nothing prevents the hub (or a later xhatter call) from +trying the same `xhat` again a few iterations down the road. + +Issue #601 proposes: when an xhatter discovers infeasibility, **optionally +generate a feasibility cut** and push it into **every scenario** so the +same `xhat` (and a neighborhood of it, for continuous variables) does +not get revisited. This is a classic no-good / Farkas-dual +feasibility-cut pattern from Benders-style decomposition, applied +inside the mpi-sppy xhat evaluation loop. + +## Non-goals + +- Replacing the main `lshaped` / `fwph` cut machinery. We add feasibility + cuts on the PH side; we do not change how L-shaped itself produces + its optimality cuts. +- Changing default behavior. The feature is off by default + (`--xhat-feasibility-cuts-count 0`). + +## Stage support + +The two-stage-only restriction in `cross_scen_extension` and +`lshaped_cuts` is about the L-shaped **cut structure** +(`eta_s >= const + coef · x_root`, one scalar recourse cost per +scenario) — that formulation breaks when the recourse cost is spread +across a tree of per-node contributions. + +A **no-good cut** (first-milestone scope here) carries none of that +structure: it is just + +``` +sum_{i: xhat_i = 1} (1 - x_i) + sum_{i: xhat_i = 0} x_i >= 1 +``` + +over whichever binary nonants the xhatter fixed. The *cut form* itself +is valid in any number of stages and does not depend on the +recourse-cost formulation. **The implementation in V1 is two-stage +only**, however, because of how the cut is installed. + +The hub installer (`XhatFeasibilityCutExtension._install_cuts`) takes +the flat coefficient row produced by the spoke and applies it +positionally against every local scenario's +`s._mpisppy_data.nonant_indices`. In two-stage every scenario shares +the same ROOT nonants under nonanticipativity, so the same row applied +to every scenario is mathematically consistent. In multi-stage, +scenarios on different branches have different per-stage-2+ variables +at the deeper indices, so the same row applied positionally lands +coefficients on unrelated variables — the resulting "cut" on a +divergent scenario is meaningless. A multi-stage-correct installer +must group each cut by the source scenario's branch and install it +only on scenarios sharing that branch through the cut's deepest node. +That work is deferred to a follow-up milestone alongside V2. + +V1 therefore hard-fails at hub setup time when +`opt.multistage` is true: + +``` +RuntimeError: --xhat-feasibility-cuts-count > 0 is two-stage only +in V1. Multi-stage support is planned as a follow-up milestone +(the install side needs to group cuts by scenario branch). +``` + +## Blueprint: cross-scenario cuts + +We already have a working precedent for "spoke generates cuts, hub +installs them into every scenario". The relevant pieces: + +### Flat serialization over the shared-memory window + +- `Field.CROSS_SCENARIO_CUT = 300` in `mpisppy/cylinders/spwindow.py`, + pre-sized as `nscen × (nonant_len + 2)`. +- Row format: `[constant, eta_coef, nonant_coefs...]` — one row per + target scenario, fixed width `1 + 1 + nonant_len`. +- The spoke writes the whole `all_coefs` buffer per iteration; + the hub-side consumer checks `is_new()` and unpacks. + +### Dedicated spoke: `CrossScenarioCutSpoke` + +- Declares `send_fields = (..., CROSS_SCENARIO_CUT)` and + `receive_fields = (..., NONANT, CROSS_SCENARIO_COST)`. +- Builds a pseudo-root `ConcreteModel` with first-stage variable copies, + attaches an `LShapedCutGenerator` (thin wrapper over + `pyomo.contrib.benders.benders_cuts.BendersCutGeneratorData`), adds + every scenario as a Benders subproblem, and reuses the Benders + `generate_cut()` output. +- On each iteration it reads the current nonants and etas broadcast by + the hub, computes the "farthest" xhat, and emits cuts. + +### Hub-side extension: `CrossScenarioExtension` + +- `sync_with_spokes()` pulls the cut buffer, calls `make_cuts(coefs)`. +- `make_cuts` iterates rows, builds a `pyomo.core.expr.numeric_expr + .LinearExpression`, and appends to + `s._mpisppy_model.benders_cuts[outer_iter, scen_name]` (an + `IndexedConstraint`). +- Persistent-solver awareness: `persistent_solver.add_constraint(...)` + for every added cut, so the underlying solver state stays current. +- Tracks a best-bound constraint (`inner_bound_constr`) alongside the + cuts; that part is cross-scen-specific and does **not** apply here. + +## The pyomo Benders gap + +`pyomo.contrib.benders.benders_cuts.BendersCutGeneratorData.generate_cut()` +**hard-errors** when a subproblem is not optimal: + +```python +if res.solver.termination_condition != pyo.TerminationCondition.optimal: + raise RuntimeError('Unable to generate cut because subproblem failed + to converge.') +``` + +That is, it generates **optimality cuts only** from optimally-solved +subproblems. `mpisppy/utils/lshaped_cuts.py` inherits the same +limitation — its `LShapedCutGenerator` is a thin subclass of pyomo's +`BendersCutGeneratorData`. + +Issue #601 is fundamentally about the **infeasibility** case — that's +where we want the cut. So "reuse `lshaped_cuts`" is necessary but not +sufficient. We need one of: + +1. **Upstream extension** — add feasibility-cut support to + `pyomo.contrib.benders.benders_cuts` so that an infeasible + subproblem yields a Farkas-dual cut. This is the most correct + option; it is also a Pyomo PR with its own review cycle. +2. **mpi-sppy-side extraction** — drop down to solver-specific APIs + (`gurobipy.Model.FarkasProof`, `cplex.solution.get_status`, + `xpress.getdualray`) inside mpi-sppy and build the cut ourselves, + bypassing pyomo Benders for the feasibility branch. Faster to land, + but fragile across solver versions and adds a new solver-sniffing + surface. +3. **Scope-limited first version** — start with **no-good cuts, + valid only when every first-stage (nonant) variable is binary**. + That needs no dual information at all: the infeasible `xhat_bin` + yields the cut + `sum_{i: xhat_i = 1} (1 - x_i) + sum_{i: xhat_i = 0} x_i >= 1`. + Trivial to code, covers UC / scheduling / lot-sizing. If any + first-stage nonant is integer or continuous, the no-good cut is + invalid — so enabling the feature on such a model is a hard + error, not a silent no-op (see "Startup check" below). The + integer case needs a different cut form (e.g. `|x - xhat| >= 1` + with auxiliary binaries) and the continuous case needs Farkas + duals; both are out of scope for the first milestone. + +**Recommendation:** pursue (1) as the north star but ship (3) as the +first deliverable so the plumbing lands independently of the Pyomo +review cycle. (2) is a trap — avoid. + +### Scope across milestones + +It helps to be explicit about *which variables* each approach can +handle. The Farkas cut (V2) is linear in the first-stage (complicating) +variables, so their domain does not affect the cut's validity — only +the outer solver's convergence behavior. The real restriction in V2 +is on the **second stage**: Farkas duality requires an LP recourse. +Integer recourse needs a further generalization (integer L-shaped / +Laporte-Louveaux / combinatorial Benders), which is a separate +milestone on top of V2. + +| | V1: no-good | V2: Farkas (pyomo extension) | V3: integer L-shaped | +|---|---|---|---| +| 1st-stage **binary** | ✓ | ✓ | ✓ | +| 1st-stage **integer** | ✗ | ✓ | ✓ | +| 1st-stage **continuous** | ✗ | ✓ | ✓ | +| 2nd-stage **LP** | ✓ | ✓ | ✓ | +| 2nd-stage **MIP** | ✓ | ✗ | ✓ | +| **Multi-stage** | ✗ | ✗ | ? | + +V1 tolerates MIP recourse because the no-good cut depends only on the +first-stage `xhat` values, not on any recourse dual. V2 lifts the +first-stage restriction but requires LP recourse. V3 is out of scope +for this design; it is mentioned only so the boundary is clear. +**Multi-stage support is orthogonal to V1/V2/V3** — the cut form is +fine but the installer must group cuts by branch (see "Stage support" +above). Tracked as a separate follow-up; the V1 setup_hub gate makes +sure no run silently produces invalid multi-stage cuts in the +meantime. + +## Proposed architecture + +Match the cross-scen split: one new shared-memory field, spoke-side +generation, hub-side extension for installation. + +### New shared-memory field + +`Field.XHAT_FEASIBILITY_CUT` in `spwindow.py`, pre-sized as +`cuts_per_iter × (1 + nonant_len) + 1`: + +- Row format: `[rhs_constant, nonant_coef_1, ..., nonant_coef_N]`. + No eta column — feasibility cuts do not involve the recourse-cost + variable. +- Trailing slot: the number of cuts actually written this round + (so unused rows are ignored). + +`cuts_per_iter` is set by `cfg.xhat_feasibility_cuts_count` at buffer +registration time. + +### CLI + +``` +--xhat-feasibility-cuts-count INT (default 0, meaning off) +``` + +Interpretation: maximum number of feasibility cuts generated per +xhatter iteration. Zero disables the whole feature. The flag +doubles as on/off switch and buffer sizer. + +### Spoke-side generation + +Modify the **existing** xhatter spokes (`xhatlooper_bounder`, +`xhatshufflelooper_bounder`, `xhatspecific_bounder`). They already +know when their candidate is infeasible in a scenario; we add a +small path: + +1. If `cfg.xhat_feasibility_cuts_count > 0`, register + `Field.XHAT_FEASIBILITY_CUT` as a send field at startup. +2. On a `solve_loop` termination whose `termination_condition` is one + of `infeasible` or `infeasibleOrUnbounded` for some scenario, + build a cut row for that scenario: + - **Version 1 (binary no-good)**: produce the no-good row directly + from the `xhat` values. The "all nonants are binary" precondition + is enforced at setup time (see Startup check below), so by the + time we reach this path we can rely on it — no per-iteration + re-check, no silent skip. + - **Version 2 (pyomo-upstream Farkas)**: call the (to-be-added) + infeasibility branch of pyomo Benders to get Farkas duals and + pack them into the same row format. +3. Pack up to `cuts_per_iter` rows into the send buffer; write the + actual count into the trailing slot; `put_send_buffer(...)`. + +Most of the "set up a pseudo-root model with first-stage copies" code +from `CrossScenarioCutSpoke.prep_cs_cuts` is reusable. + +### Hub-side extension + +New `XhatFeasibilityCutExtension` in +`mpisppy/extensions/xhat_feasibility_cut_extension.py`, modeled on +`CrossScenarioExtension` but stripped of the eta / inner-bound +bookkeeping: + +- `register_send_fields` / `register_receive_fields`: register + receive of `Field.XHAT_FEASIBILITY_CUT` from the xhatter spoke. +- `setup_hub`: on each local scenario, attach an empty + `pyo.ConstraintList` (or `IndexedConstraint` keyed by + `(cut_batch_iter, cut_row_idx)`) named + `s._mpisppy_model.xhat_feasibility_cuts`. +- `sync_with_spokes`: read the buffer. If `is_new()`, unpack the + count-header and the rows, build `LinearExpression` per row, + append the **same cut to every local scenario's** + `xhat_feasibility_cuts`, and call + `persistent_solver.add_constraint(...)` as cross-scen does. Because + the cut is on first-stage (nonant) variables, nonanticipativity + makes it valid globally — there is no version that cuts only the + infeasible scenario. +- Bundle-safe: when a local "scenario" is a proper bundle, push the + cut into each per-scenario block inside the bundle (the nonants are + duplicated across the scenario blocks inside the bundle EF, so the + cut has to land on each). This is the same class of gotcha PR #669 + is addressing on the `nonant_cost_coeffs` side — solve it once + here, up front. + +### Startup check: first-stage must be fully binary + +The first-milestone no-good cut is only valid when **every** first-stage +nonant is binary. A mix of binary and continuous nonants is not safe +to treat as "no-good on the binary part only": the continuous part can +be perturbed, so such a cut does not actually exclude the infeasible +xhat — it would silently produce an invalid relaxation. + +To prevent misuse, `XhatFeasibilityCutExtension.setup_hub` scans +`s._mpisppy_data.nonant_indices` on an arbitrary local scenario and +confirms that every variable is binary (`v.is_binary()` or +`v.domain in {Binary, BooleanSet}`). The scan covers every node in +`_mpisppy_node_list` — though in V1 a separate prior check (see +"Stage support" above) rejects multi-stage entirely, so the per-node +scan is mostly defensive against a future V1 extension. If any +non-binary nonant is present, raise + +``` +RuntimeError( + "--xhat-feasibility-cuts-count > 0 requires every first-stage " + "(nonant) variable to be binary; found non-binary nonant " + f"{v.name!r} with domain {v.domain}. The first-milestone " + "feasibility-cut generator is no-good-only. Support for integer " + "and continuous first-stage variables is planned as a follow-up " + "milestone (pyomo Benders / Farkas extension)." +) +``` + +Fail fast at setup time rather than at the first infeasibility — the +user sees the error before any work is done. + +### Plumbing touchpoints + +- `mpisppy/utils/config.py` — add `xhat_feasibility_cuts_count` to an + appropriate config bundle (likely `cfg.spoke_xhat_args()` or a new + `cfg.xhat_feasibility_cut_args()`). +- `mpisppy/utils/cfg_vanilla.py` — flag in the xhatter spoke dicts so + the extension and the send-field are registered together. +- `mpisppy/cylinders/spwindow.py` — new `Field` enum value and + length formula. +- `mpisppy/extensions/xhat_feasibility_cut_extension.py` — new hub + extension. +- `mpisppy/cylinders/xhatlooper_bounder.py` (etc.) — add the cut + emission path gated on `cfg.xhat_feasibility_cuts_count > 0`. + +## Code sharing / refactoring opportunities + +The user asked explicitly: can we share code with the cross-scen path? +Three candidates, in descending value: + +1. **Cut installation into scenarios, with persistent-solver and + bundle awareness.** Today this logic is inline in + `CrossScenarioExtension.make_cuts`; the xhat-feas version will want + the same behavior. Factor into + `mpisppy.utils.cut_installer.install_linear_cut(scenario, expr, + key, persistent_solver=...)`. Both extensions call it. High + payoff: it's the trickiest bit (persistent-solver management plus + the bundle duplication). + +2. **Pseudo-root model construction.** `CrossScenarioCutSpoke + .prep_cs_cuts` builds a `ConcreteModel` with first-stage-variable + copies keyed by an "index against later" map. Xhat-feas spoke side + needs the same kind of harness if/when we pursue the pyomo-Benders + route. Factor into + `mpisppy.utils.pseudo_root.build_root_from_nonants(opt) -> (root, + root_vars, nonant_vid_to_copy_map)`. + +3. **Flat-buffer packer/unpacker.** Both cross-scen and xhat-feas use + `[constant, ...coefs]` rows. A small + `mpisppy.utils.cut_buffer.pack_rows(...) / unpack_rows(...)` + helper de-duplicates the row arithmetic. Lower payoff — the + arithmetic is straightforward and inlined it is arguably more + readable — but it's essentially free to write if we do (1). + +Recommended refactoring sequence: land the new xhat-feas feature first +(without factoring), **then** do a follow-up PR that extracts the +shared helpers and updates cross-scen to use them. This avoids +bundling a behavior change with a refactor. + +## Open questions + +- **Cut pool management.** Every xhatter iteration can emit cuts; + without pruning the subproblems accumulate constraints indefinitely + and solver times eventually bloat. `CrossScenarioExtension` has the + same unbounded-growth behavior today (its `benders_cuts` are keyed + by `(outer_iter, scen_name)` and never removed). Tracked as a + separate item, scoped to both features so they share one strategy: + **[issue #670](https://github.com/Pyomo/mpi-sppy/issues/670)**. For + this milestone we rely on the per-iter cap in + `--xhat-feasibility-cuts-count` to keep growth linear. + +## First-milestone scope (what a PR would actually deliver) + +1. CLI flag + buffer registration + send/receive field — wired end to + end but with the emit path producing zero cuts. +2. No-good-cut emission for problems whose first-stage (nonant) + variables are **all binary**. Setup-time check raises `RuntimeError` + otherwise, so the feature never silently produces an invalid cut. +3. Hub-side install of cuts into every local scenario / bundle block. + Two-stage only in V1; the install loop applies one row positionally + to every scenario's nonants, which is correct under + nonanticipativity for two-stage but not for multi-stage. A + multi-stage-correct installer (group cuts by branch) is a follow-up + milestone; in the meantime `setup_hub` raises if `opt.multistage`. +4. Regression tests: + (a) a two-stage model with a binary first-stage where one specific + xhat is infeasible — verify the cut is added and the infeasible + xhat is not revisited; + (b) a bundled version of the same; + (c) a negative test: enabling the feature on a multi-stage model + raises the expected `RuntimeError` at setup; + (d) a negative test: enabling the feature on a model with a + continuous nonant raises the expected `RuntimeError` at setup. +5. Documentation: a **new** user-facing page at + `doc/src/xhat_feasibility_cuts.rst`, added to `doc/src/index.rst`'s + toctree next to the other cylinder / extension pages. The page + must cover: + - What the feature does and when to use it (non-complete-recourse + problems with binary first-stage). + - The CLI flag `--xhat-feasibility-cuts-count`. + - The binary-only restriction, with the exact `RuntimeError` text + users will see if they enable the feature on a non-binary model. + - A brief note on interaction with proper bundles. + - Forward pointer to issue #670 for cut-pool management. + + Also add a one-sentence cross-reference from each of the xhatter + cylinder docs (wherever those live today) pointing at the new + page, so readers browsing xhatter options discover it. + +That is a discrete, reviewable PR. The pyomo-Benders Farkas extension +(which lifts the binary-only restriction, two-stage-only) and the +cross-scen refactor are natural follow-ups. diff --git a/doc/src/extensions.rst b/doc/src/extensions.rst index ad71c0c41..73cf676ba 100644 --- a/doc/src/extensions.rst +++ b/doc/src/extensions.rst @@ -28,6 +28,9 @@ command-line flags: - ``--grad-rho`` -- activates gradient-based rho (see :ref:`rho_setting`) - ``--use-norm-rho-updater`` -- activates the norm rho updater - ``--use-primal-dual-rho-updater`` -- activates the primal-dual rho updater +- ``--xhat-feasibility-cuts-count `` -- feasibility cuts from the + xhatter when its candidate is infeasible; binary first-stage only + (see :ref:`xhat_feasibility_cuts`) The rest of this help file describes extensions released with mpisppy along with some hints for including them in your own cylinders driver program. diff --git a/doc/src/index.rst b/doc/src/index.rst index 2f28f3193..8929ee5b8 100644 --- a/doc/src/index.rst +++ b/doc/src/index.rst @@ -57,6 +57,7 @@ MPI is used. jensens.rst feasible_xhat.rst xhat_from_file.rst + xhat_feasibility_cuts.rst smps.rst agnostic.rst generic_admm.rst diff --git a/doc/src/spokes.rst b/doc/src/spokes.rst index e250af281..532bde125 100644 --- a/doc/src/spokes.rst +++ b/doc/src/spokes.rst @@ -148,7 +148,15 @@ An example is shown in ``examples.hydro.hydro_cylinders.py`` (this particular ex is intended to show the coding, not normal behavior. It is sort of an edge case: including this option causes the upper bound to immediately be Z*) - +Feasibility cuts from xhatters +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For problems without complete recourse, any xhatter can optionally +emit a no-good feasibility cut when its candidate :math:`\hat{x}` is +infeasible in some scenario; the hub then installs the cut into every +scenario so the same :math:`\hat{x}` is not revisited. See +:ref:`xhat_feasibility_cuts`. + slam_heuristic ^^^^^^^^^^^^^^ diff --git a/doc/src/xhat_feasibility_cuts.rst b/doc/src/xhat_feasibility_cuts.rst new file mode 100644 index 000000000..11b127659 --- /dev/null +++ b/doc/src/xhat_feasibility_cuts.rst @@ -0,0 +1,163 @@ +.. _xhat_feasibility_cuts: + +Xhat Feasibility Cuts +===================== + +For two-stage problems **without complete recourse**, a +candidate first-stage solution ``xhat`` proposed by an xhatter spoke +(``xhatlooper``, ``xhatshufflelooper``, ``xhatspecific``, +``xhatxbar``) can turn out to be infeasible in one or more scenarios. +Today those candidates are simply discarded and nothing prevents the +same ``xhat`` from being proposed again a few iterations later. + +When this feature is enabled, an xhatter that detects such an +infeasibility emits a **no-good feasibility cut** on the first-stage +variables, and the hub installs that cut into every scenario's +subproblem. The result is that the same ``xhat`` (and any other +assignment with the same pattern on binaries) is excluded from future +consideration. + +This is a first-milestone implementation: it is **two-stage only** +and **valid only when every first-stage (nonant) variable is binary**. +See :ref:`xhat_feas_cuts_boundaries` below. + +Enabling the Feature +-------------------- + +``generic_cylinders`` exposes a single integer flag: + +.. code-block:: bash + + --xhat-feasibility-cuts-count N + +where ``N`` is the maximum number of feasibility cuts the xhatter may +emit per iteration. The default is ``0``, which disables the feature +entirely. Any positive integer both turns the feature on and sizes the +shared-memory buffer the xhatter uses to send cuts. + +Nothing else needs to be specified. The hub-side installer +(``mpisppy.extensions.xhat_feasibility_cut_extension.XhatFeasibilityCutExtension``) +is attached automatically through :ref:`cfg_vanilla ` when the +flag is positive. + +Example +------- + +.. code-block:: bash + + python -m mpisppy.generic_cylinders \ + --module-name my_binary_first_stage_model \ + --num-scens 10 \ + --solver-name gurobi \ + --max-iterations 50 \ + --default-rho 1.0 \ + --lagrangian \ + --xhatshuffle \ + --xhat-feasibility-cuts-count 3 + +The cap of 3 says "emit at most 3 cuts per xhatter iteration". Cuts +accumulate across iterations inside each scenario's +``_mpisppy_model.xhat_feasibility_cuts`` constraint container. + +.. _xhat_feas_cuts_boundaries: + +Scope and Limitations +--------------------- + +The first-milestone cut is the textbook no-good inequality + +.. math:: + + \sum_{i:\, \hat{x}_i = 1} (1 - x_i) + \sum_{i:\, \hat{x}_i = 0} x_i \;\geq\; 1 + +which is valid only when every :math:`x_i` is binary. If any +first-stage nonant is integer (not bounded to :math:`\{0, 1\}`) or +continuous, the cut cannot exclude the infeasible ``xhat`` correctly, +so the feature **refuses to run** rather than silently generate +invalid relaxations: + +.. code-block:: text + + RuntimeError: --xhat-feasibility-cuts-count > 0 requires every + first-stage (nonant) variable to be binary; found non-binary nonant + '' (key (, )) on scenario '' with + domain . The first-milestone feasibility-cut generator is + no-good-only. Support for integer and continuous first-stage + variables is planned as a follow-up milestone (pyomo Benders / + Farkas extension). + +The error is raised at hub setup time (before any PH work begins), +so a misconfiguration is caught immediately. + +**Integer first-stage variables with bounds** :math:`[0, 1]` **are +accepted** — semantically those are binary. Declaring a var as +``pyo.Integers`` with ``bounds=(0, 1)`` works just as well as +``pyo.Binary``. + +Multi-stage +----------- + +V1 is two-stage only. Enabling +``--xhat-feasibility-cuts-count`` on a multi-stage model raises at +hub setup time: + +.. code-block:: text + + RuntimeError: --xhat-feasibility-cuts-count > 0 is two-stage only + in V1. Multi-stage support is planned as a follow-up milestone + (the install side needs to group cuts by scenario branch). See + doc/designs/xhat_feasibility_cuts_design.md. + +The cut row encodes coefficients positionally against each scenario's +``nonant_indices``. In two-stage every scenario shares the same ROOT +nonants under nonanticipativity, so applying the same row to every +scenario is consistent. In multi-stage, scenarios on different +branches have different per-stage-2+ variables at the deeper indices, +so the same row applied positionally lands coefficients on unrelated +variables. A multi-stage-correct installer needs to group cuts by +branch; that work is deferred to a follow-up milestone. + +Interaction with Proper Bundles +------------------------------- + +The hub installer appends cuts to each scenario's +``xhat_feasibility_cuts`` constraint container. When a scenario +object is actually a proper bundle, the cut is installed against the +bundle's canonical nonant set (``s._mpisppy_data.nonant_indices``) +exactly as the cross-scenario cut machinery does. Nonanticipativity +inside the bundle ensures the cut takes effect on every per-scenario +block. + +Follow-up Milestones +-------------------- + +- **Multi-stage support.** A multi-stage-correct installer needs to + group each cut by the branch of the source scenario and install it + only on scenarios sharing that branch through the cut's deepest + node. The current installer applies one cut row uniformly to every + scenario, which is only valid in two-stage; that's why V1 hard-fails + at setup on a multi-stage model. Tracking as a follow-up alongside + V2. +- Lifting the binary-only restriction requires generating **Farkas + feasibility cuts** from the dual ray of an infeasible second-stage + LP. The upstream + ``pyomo.contrib.benders.benders_cuts.BendersCutGeneratorData`` + currently only produces optimality cuts; supporting the + infeasibility case is a Pyomo PR. Once that lands, the xhatter will + be able to emit valid cuts for integer and continuous first-stage + variables as well (LP recourse only). +- **Cut-pool management** is deferred to + `issue #670 `_. Today + cuts accumulate indefinitely in the per-scenario constraint + container, the same way cross-scenario cuts do in + ``CrossScenarioExtension``. For long runs, use the + ``--xhat-feasibility-cuts-count`` cap to bound the per-iteration + growth. + +See Also +-------- + +- :ref:`Spokes` — overview of the xhat spokes. +- :ref:`Extensions` — the broader extension mechanism. +- ``doc/designs/xhat_feasibility_cuts_design.md`` — the design document with + the full milestone plan and the ``V1/V2/V3`` scope table. diff --git a/examples/run_all.py b/examples/run_all.py index 555dac1d2..ac6538d0e 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -217,6 +217,18 @@ def do_one_mmw(dirname, modname, runefstring, npyfile, mmwargstring): f"--num-scens=3 --solver-name={solver_name} " f"--max-iterations=3 --default-rho=1 --lagrangian --xhatshuffle " f"--output-dir=solutions_ws {usar_problem_args}") + # usar has a binary first-stage (is_active_depot), so it exercises + # --xhat-feasibility-cuts-count end-to-end: buffer registration on + # the xhatter spoke, the hub extension's binary-only startup check, + # and the sync loop. The xhatter may or may not actually hit an + # infeasibility on this small instance — the cut-emission path is + # unit-tested separately in test_xhat_feasibility_cuts.py — but this + # smoke entry guards the full-run plumbing from regressions. + do_one("usar", "wheel_spinner.py", 3, + f"--num-scens=3 --solver-name={solver_name} " + f"--max-iterations=3 --default-rho=1 --lagrangian --xhatshuffle " + f"--xhat-feasibility-cuts-count=3 " + f"--output-dir=solutions_ws_feas {usar_problem_args}") # netdes # NOTE: Pyomo OBBT does not support persistent solvers as of Aug 2025 diff --git a/examples/usar/wheel_spinner.py b/examples/usar/wheel_spinner.py index 1fb994851..3b8c5a4a1 100644 --- a/examples/usar/wheel_spinner.py +++ b/examples/usar/wheel_spinner.py @@ -70,6 +70,9 @@ def wheel_spinner( vanilla_fn = vanilla.aph_hub if cnfg["run_async"] else vanilla.ph_hub hub_dict = vanilla_fn(*vanilla_args, **vanilla_kwargs) + # Optional feasibility cuts from xhatters (requires binary first-stage). + hub_dict = vanilla.add_xhat_feasibility_cuts(hub_dict, cnfg) + spoke_dicts = [] for spoke_name in SUPPORTED_SPOKES: if getattr(cnfg, spoke_name): @@ -90,6 +93,7 @@ def main() -> None: cnfg.aph_args() for spoke_name in SUPPORTED_SPOKES: getattr(cnfg, spoke_name + "_args")() + cnfg.xhat_feasibility_cut_args() add_common_args(cnfg) cnfg.parse_command_line() diff --git a/mpisppy/cylinders/spwindow.py b/mpisppy/cylinders/spwindow.py index 003694262..bb21e7b97 100644 --- a/mpisppy/cylinders/spwindow.py +++ b/mpisppy/cylinders/spwindow.py @@ -36,6 +36,7 @@ class Field(enum.IntEnum): BEST_XHAT=600 # buffer having the best xhat and its total cost per scenario RECENT_XHATS=601 # buffer having some recent xhats and their total cost per scenario XFEAS=602 # buffer having a feasible x and its total cost per scenario + XHAT_FEASIBILITY_CUT=700 # feasibility cuts emitted by xhat spokes WHOLE=1_000_000 @@ -47,6 +48,9 @@ class Field(enum.IntEnum): # these could be modified by the user... field_length_components.total_number_recent_xhats = pyo.Param(mutable=True, initialize=10, within=pyo.NonNegativeIntegers) +# max cuts an xhat spoke may emit per iteration (set by the xhatter +# spoke from cfg.xhat_feasibility_cuts_count before field registration) +field_length_components.xhat_feasibility_cuts_per_iter = pyo.Param(mutable=True, initialize=0, within=pyo.NonNegativeIntegers) _field_lengths = { Field.SHUTDOWN : 1, @@ -65,6 +69,9 @@ class Field(enum.IntEnum): Field.BEST_XHAT : field_length_components._local_nonant_length + field_length_components._local_scenario_length, Field.RECENT_XHATS : field_length_components.total_number_recent_xhats * (field_length_components._local_nonant_length + field_length_components._local_scenario_length), Field.XFEAS: field_length_components._local_nonant_length + field_length_components._local_scenario_length, + # rows: [constant, nonant_coef_1, ..., nonant_coef_N]; trailing slot holds the + # actual number of cuts written this batch (0..per_iter). + Field.XHAT_FEASIBILITY_CUT : field_length_components.xhat_feasibility_cuts_per_iter * (field_length_components._total_number_nonants + 1) + 1, } @@ -81,6 +88,10 @@ def __init__(self, opt): field_length_components._local_scenario_length.value = len(opt.local_scenarios) field_length_components._total_number_nonants.value = opt.nonant_length field_length_components._total_number_scenarios.value = len(opt.local_scenarios) + # user-tunable cap on feasibility cuts per iteration (0 = off) + field_length_components.xhat_feasibility_cuts_per_iter.value = int( + opt.options.get("xhat_feasibility_cuts_count", 0) + ) self._field_lengths = {k : pyo.value(v) for k, v in _field_lengths.items()} diff --git a/mpisppy/cylinders/xhatbase.py b/mpisppy/cylinders/xhatbase.py index 2ebed1c7f..5964ad80a 100644 --- a/mpisppy/cylinders/xhatbase.py +++ b/mpisppy/cylinders/xhatbase.py @@ -12,10 +12,18 @@ import math import mpisppy.cylinders.spoke as spoke +from mpisppy.cylinders.spwindow import Field from mpisppy.utils.xhat_eval import Xhat_Eval class XhatInnerBoundBase(spoke.InnerBoundNonantSpoke): + # Advertise the optional feasibility-cut send buffer. When + # cfg.xhat_feasibility_cuts_count == 0 the buffer is size 1 (just the + # trailing count slot) and never written; see spwindow.FieldLengths + # and XhatBase._try_one. + send_fields = (*spoke.InnerBoundNonantSpoke.send_fields, + Field.XHAT_FEASIBILITY_CUT,) + @abc.abstractmethod def xhat_extension(self): raise NotImplementedError @@ -98,7 +106,32 @@ def _try_file_xhat(self): if self.cylinder_rank == 0: print(f"[xhat-from-file] evaluating {path!r} " f"({expected} nonants)") - Eobj = self.opt.evaluate(nonant_cache) + try: + Eobj = self.opt.evaluate(nonant_cache) + except Exception: + Eobj = None + # Same predicate XhatBase._try_one uses to detect infeasibility + # in some scenario. When True and feasibility cuts are enabled, + # pack a no-good cut so the hub extension can install it on + # every scenario — same path the regular xhatter takes. + infeasP = 0.0 + try: + infeasP = self.opt.no_incumbent_prob() + except Exception: + pass + if infeasP != 0.: + from mpisppy.extensions.xhatbase import pack_no_good_feasibility_cut + try: + emitted = pack_no_good_feasibility_cut(self.opt) + except Exception as e: + emitted = False + if self.cylinder_rank == 0: + print(f"[xhat-from-file] feasibility-cut emit failed: {e!r}") + if self.cylinder_rank == 0: + tag = "emitted" if emitted else "skipped" + print(f"[xhat-from-file] candidate infeasible " + f"(infeasP={infeasP}); feasibility cut {tag}") + Eobj = None # Restore nonants so the main loop starts from clean state. self.opt._restore_nonants() if Eobj is not None and math.isfinite(Eobj): diff --git a/mpisppy/extensions/xhat_feasibility_cut_extension.py b/mpisppy/extensions/xhat_feasibility_cut_extension.py new file mode 100644 index 000000000..a13152c42 --- /dev/null +++ b/mpisppy/extensions/xhat_feasibility_cut_extension.py @@ -0,0 +1,202 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +"""Hub-side extension that receives feasibility cuts emitted by xhat spokes +and installs them into every local scenario so an infeasible xhat is not +revisited. + +See doc/designs/xhat_feasibility_cuts_design.md (and the user-facing +doc/src/xhat_feasibility_cuts.rst) for the design and usage. +""" + +import pyomo.environ as pyo +from pyomo.core.expr.numeric_expr import LinearExpression + +import mpisppy.utils.sputils as sputils +from mpisppy.cylinders.spwindow import Field +from mpisppy.extensions.extension import Extension + + +class XhatFeasibilityCutExtension(Extension): + """Install feasibility cuts produced by xhat spokes. + + Wire format (matches Field.XHAT_FEASIBILITY_CUT): + [ row_0, row_1, ..., row_{K-1}, count ] + where each row is + [ rhs_constant, nonant_coef_1, ..., nonant_coef_N ] + and ``count`` is the number of valid rows this batch (0..K). + + A row encodes the constraint + rhs_constant + sum_i nonant_coef_i * x_i >= 0 + over the local scenario's nonants (keyed by + ``s._mpisppy_data.nonant_indices`` in insertion order). + + Binary-only first-stage is a precondition; setup_hub hard-fails + if any nonant is not binary. See the design doc. + """ + + def __init__(self, spbase_object): + super().__init__(spbase_object) + self._cuts_per_iter = int( + self.opt.options.get("xhat_feasibility_cuts_count", 0) + ) + if self._cuts_per_iter <= 0: + raise RuntimeError( + "XhatFeasibilityCutExtension was attached but " + "options['xhat_feasibility_cuts_count'] is not a positive " + "integer. Enable the feature via --xhat-feasibility-cuts-count " + "or remove the extension from the hub." + ) + self._nonant_len = None # filled in at setup_hub + self._row_len = None + self._recv_buffer = None # filled in at register_receive_fields + self._install_counter = 0 # monotonic key for the ConstraintList + + # ---- two-stage-only precondition (V1) -------------------------------- + + def _assert_two_stage(self): + """V1 ships two-stage only. + + Why: the no-good cut row encodes coefficients positionally + against ``s._mpisppy_data.nonant_indices``. In two-stage every + scenario shares the same ROOT nonants (NAC), so the same + coefficient vector applied to every scenario's nonants is + consistent. In multi-stage, scenarios on different branches + have different per-stage-2+ variables at the deeper indices, + so the same row applied positionally to each scenario lands + coefficients on unrelated variables. A multi-stage-correct + installer needs to group cuts by branch; that is deferred to a + follow-up milestone. + """ + if getattr(self.opt, "multistage", False): + raise RuntimeError( + "--xhat-feasibility-cuts-count > 0 is two-stage only in " + "V1. Multi-stage support is planned as a follow-up " + "milestone (the install side needs to group cuts by " + "scenario branch). See " + "doc/designs/xhat_feasibility_cuts_design.md." + ) + + # ---- binary-only precondition check ---------------------------------- + + @staticmethod + def _is_binary_var(v): + """Return True if v is a binary pyomo VarData. + + Accepts either ``v.is_binary()`` (if present) or a domain check + (``Binary``, ``BooleanSet``), and treats a 0/1-bounded integer + variable as binary so the user isn't punished for a common + modeling style. + """ + if v.is_binary(): + return True + if v.is_integer(): + lb, ub = v.bounds + if lb is not None and ub is not None and lb >= 0 and ub <= 1: + return True + return False + + def _assert_all_nonants_binary(self): + """Scan every nonant on every local scenario / node; raise on non-binary. + + For a bundle ``s``, ``s._mpisppy_data.nonant_indices`` is the + canonical (consolidated) nonant set — that is what cuts will be + installed against, so that is also what we validate. + """ + for sname, s in self.opt.local_scenarios.items(): + for ndn_i, v in s._mpisppy_data.nonant_indices.items(): + if not self._is_binary_var(v): + raise RuntimeError( + "--xhat-feasibility-cuts-count > 0 requires every " + "first-stage (nonant) variable to be binary; found " + f"non-binary nonant {v.name!r} (key {ndn_i}) on " + f"scenario {sname!r} with domain {v.domain}. " + "The first-milestone feasibility-cut generator is " + "no-good-only. Support for integer and continuous " + "first-stage variables is planned as a follow-up " + "milestone (pyomo Benders / Farkas extension)." + ) + + # ---- Extension hooks ------------------------------------------------- + + def setup_hub(self): + self._assert_two_stage() + self._assert_all_nonants_binary() + # Nonant count per local scenario; cuts come packed against this. + # All local scenarios must have the same nonant count (that is + # already a PH invariant), so read from one. + any_s = next(iter(self.opt.local_scenarios.values())) + self._nonant_len = len(any_s._mpisppy_data.nonant_indices) + self._row_len = 1 + self._nonant_len + # ConstraintList per local scenario to accumulate incoming cuts. + for s in self.opt.local_scenarios.values(): + s._mpisppy_model.xhat_feasibility_cuts = pyo.Constraint(pyo.Any) + + def register_send_fields(self): + # We do not send anything; the spoke is the sender. + return + + def register_receive_fields(self): + spcomm = self.opt.spcomm + ranks = spcomm.fields_to_ranks.get(Field.XHAT_FEASIBILITY_CUT, []) + if not ranks: + # No xhatter spoke is emitting cuts this run; extension is a no-op. + self._recv_buffer = None + return + # For now assume a single emitting spoke (mirrors cross-scen). + assert len(ranks) == 1, ( + "XhatFeasibilityCutExtension expects exactly one spoke to emit " + f"Field.XHAT_FEASIBILITY_CUT; found {len(ranks)}." + ) + self._recv_buffer = spcomm.register_recv_field( + Field.XHAT_FEASIBILITY_CUT, ranks[0] + ) + + def sync_with_spokes(self): + if self._recv_buffer is None: + return + if self._recv_buffer.is_new(): + self._install_cuts(self._recv_buffer.array()) + + # ---- Cut installation ------------------------------------------------ + + def _install_cuts(self, buf): + """Unpack the buffer and append each row to every local scenario. + + Buffer layout: K rows of length row_len followed by a scalar + ``count``. Rows 0..count-1 are valid; the rest are ignored. + """ + n_cuts = int(buf[-1]) + if n_cuts <= 0: + return + row_len = self._row_len + for k in range(n_cuts): + row = buf[k * row_len : (k + 1) * row_len] + rhs_constant = float(row[0]) + coefs = [float(c) for c in row[1:]] + # A defensive zero-row check — the spoke may have written + # fewer cuts than the header claims if we ever tighten that. + if rhs_constant == 0.0 and all(c == 0.0 for c in coefs): + continue + self._install_counter += 1 + key = self._install_counter + for s in self.opt.local_scenarios.values(): + linear_vars = list(s._mpisppy_data.nonant_indices.values()) + # Constraint form: rhs_constant + sum coef_i x_i >= 0 + # LinearExpression has a single `constant=` slot; we + # pass rhs_constant as that constant and 0 as the RHS. + expr = LinearExpression( + constant=rhs_constant, + linear_coefs=coefs, + linear_vars=linear_vars, + ) + s._mpisppy_model.xhat_feasibility_cuts[key] = (0.0, expr, None) + if sputils.is_persistent(s._solver_plugin): + s._solver_plugin.add_constraint( + s._mpisppy_model.xhat_feasibility_cuts[key] + ) diff --git a/mpisppy/extensions/xhatbase.py b/mpisppy/extensions/xhatbase.py index f2deee6f9..69691755f 100644 --- a/mpisppy/extensions/xhatbase.py +++ b/mpisppy/extensions/xhatbase.py @@ -17,6 +17,69 @@ import mpisppy.utils.sputils as sputils +def pack_no_good_feasibility_cut(opt): + """Pack a no-good feasibility cut from currently-fixed nonants and send it. + + Precondition: called while every local scenario's nonant variables + are fixed to the candidate xhat (e.g. immediately after a + ``solve_loop`` that detected infeasibility, before + ``_restore_nonants``). Used by both the in-loop xhatter + (``XhatBase._try_one``) and the file-fed xhat path + (``XhatInnerBoundBase._try_file_xhat``). + + Returns True if a cut was emitted, False if the call early-returned + (feature off, no spcomm, send field not registered, or a nonant had + no value). Raises RuntimeError on a buffer-size mismatch. + + Cut form (for binary xhat with x_i in {0,1}): + constant + sum_i coef_i * x_i >= 0 + constant = (#{xhat_i = 1}) - 1 + coef_i = 1 - 2 * xhat_i (= +1 if xhat=0, -1 if xhat=1) + + Equivalent to the standard no-good inequality + sum_{xhat_i=1} (1 - x_i) + sum_{xhat_i=0} x_i >= 1 + """ + cuts_cap = int(opt.options.get("xhat_feasibility_cuts_count", 0)) + if cuts_cap <= 0: + return False + spcomm = opt.spcomm + if spcomm is None: + return False + # Local import to avoid a cylinders->extensions import cycle. + from mpisppy.cylinders.spwindow import Field + if not spcomm.is_send_field_registered(Field.XHAT_FEASIBILITY_CUT): + return False + + send_array = spcomm.send_buffers[Field.XHAT_FEASIBILITY_CUT] + buf = send_array.value_array() + any_s = next(iter(opt.local_scenarios.values())) + nonant_len = len(any_s._mpisppy_data.nonant_indices) + row_len = 1 + nonant_len + expected = cuts_cap * row_len + 1 + if buf.shape[0] != expected: + raise RuntimeError( + f"XHAT_FEASIBILITY_CUT buffer has length {buf.shape[0]}, " + f"expected {expected} = {cuts_cap} * ({nonant_len}+1) + 1" + ) + + # Pack one row (row 0); later rows zeroed so the hub's zero-check skips them. + buf.fill(0.0) + xhat_sum_ones = 0 + for idx, v in enumerate(any_s._mpisppy_data.nonant_indices.values()): + xv = v.value + if xv is None: + return False + xhat_i = 1 if xv > 0.5 else 0 + if xhat_i == 1: + xhat_sum_ones += 1 + buf[1 + idx] = 1 - 2 * xhat_i + buf[0] = xhat_sum_ones - 1 + buf[-1] = 1.0 + + spcomm.put_send_buffer(send_array, Field.XHAT_FEASIBILITY_CUT) + return True + + class XhatBase(mpisppy.extensions.extension.Extension): """ Any inherited class must implement the preiter0, postiter etc. methods @@ -228,6 +291,18 @@ def _try_one(self, snamedict, solver_options=None, verbose=False, infeasP = self.opt.no_incumbent_prob() if infeasP != 0.: + # Optional: emit a no-good feasibility cut so the hub can + # prevent revisiting this xhat. Gated on + # cfg.xhat_feasibility_cuts_count > 0 and on the send field + # actually being registered (it will be for xhat spokes that + # descend from XhatInnerBoundBase). Binary-only cut; the hub + # extension's setup_hub verifies that precondition. + try: + self._maybe_emit_feasibility_cut() + except Exception as e: + # Emission must never break the xhatter. Log and move on. + if self.cylinder_rank == 0: + print(f"[xhat feasibility cuts] emit failed: {e!r}") if restore_nonants: self.opt._restore_nonants() return None @@ -241,6 +316,15 @@ def _try_one(self, snamedict, solver_options=None, verbose=False, self.opt._restore_nonants() return obj + #********** + def _maybe_emit_feasibility_cut(self): + """Pack a no-good cut from the just-evaluated xhat and send it. + + Thin wrapper around :func:`pack_no_good_feasibility_cut`. See + that function for preconditions, gating, and the cut form. + """ + pack_no_good_feasibility_cut(self.opt) + #********** def csv_nonants(self, snamedict, fname): """ write the non-ants in csv format to files based on the file name diff --git a/mpisppy/generic/extensions.py b/mpisppy/generic/extensions.py index 46ca1d1af..3fd241e5f 100644 --- a/mpisppy/generic/extensions.py +++ b/mpisppy/generic/extensions.py @@ -36,6 +36,9 @@ def configure_extensions(hub_dict, module, cfg): if cfg.rc_fixer: vanilla.add_reduced_costs_fixer(hub_dict, cfg) + if cfg.get("xhat_feasibility_cuts_count", 0) or 0: + vanilla.add_xhat_feasibility_cuts(hub_dict, cfg) + if cfg.relaxed_ph_fixer: vanilla.add_relaxed_ph_fixer(hub_dict, cfg) diff --git a/mpisppy/generic/parsing.py b/mpisppy/generic/parsing.py index f8f780f23..ee7ab5103 100644 --- a/mpisppy/generic/parsing.py +++ b/mpisppy/generic/parsing.py @@ -145,6 +145,7 @@ def parse_args(m): cfg.xhatshuffle_args() cfg.xhatxbar_args() cfg.xhat_from_file_args() + cfg.xhat_feasibility_cut_args() cfg.norm_rho_args() cfg.primal_dual_rho_args() cfg.converger_args() diff --git a/mpisppy/tests/test_xhat_feasibility_cuts.py b/mpisppy/tests/test_xhat_feasibility_cuts.py new file mode 100644 index 000000000..2dcf3bf7a --- /dev/null +++ b/mpisppy/tests/test_xhat_feasibility_cuts.py @@ -0,0 +1,950 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +"""Unit tests for the xhat feasibility-cut plumbing. + +Covers: +- Binary-only startup check in ``XhatFeasibilityCutExtension.setup_hub`` + (positive + negative cases). +- Cut installation via the hub extension: basic, zero-count no-op, + all-zero-row skipped, and bundle-safe install against a proper bundle. +- The no-good-cut encoding produced by ``XhatBase._maybe_emit_feasibility_cut`` + exercised through a fake ``spcomm`` — verifies both the coefficient + math and the trailing count slot. +- Extension ``__init__`` rejects being attached when + ``cfg.xhat_feasibility_cuts_count`` is zero. + +End-to-end PH+spoke runs are deliberately out of scope for this file; +those require MPI and the full cylinder machinery. +""" + +import unittest + +import numpy as np +import pyomo.environ as pyo + +import mpisppy.cylinders.xhatbase as cylinder_xhatbase # noqa: F401 — ensures class-level send_fields registration is covered +import mpisppy.extensions.xhatbase as xhatbase_extension +import mpisppy.scenario_tree as stree +import mpisppy.utils.sputils as sputils +from mpisppy.cylinders.spwindow import Field +from mpisppy.extensions.xhat_feasibility_cut_extension import ( + XhatFeasibilityCutExtension, +) +from mpisppy.spbase import SPBase + + +# --------------------------------------------------------------------------- +# Minimal scenario creators (pure Pyomo, no solver required). +# --------------------------------------------------------------------------- + + +def _binary_two_stage_scenario(sname, **kwargs): + m = pyo.ConcreteModel(name=sname) + m.x = pyo.Var([0, 1, 2], domain=pyo.Binary) + m.fsc = pyo.Expression(expr=m.x[0] + m.x[1] + m.x[2]) + m.obj = pyo.Objective(expr=2 * m.x[0] + 3 * m.x[1] + m.x[2]) + sputils.attach_root_node(m, m.fsc, [m.x]) + m._mpisppy_probability = "uniform" + return m + + +def _continuous_two_stage_scenario(sname, **kwargs): + m = pyo.ConcreteModel(name=sname) + m.x = pyo.Var([0, 1, 2], bounds=(0, 1)) # continuous 0..1 + m.fsc = pyo.Expression(expr=m.x[0] + m.x[1] + m.x[2]) + m.obj = pyo.Objective(expr=m.x[0] + m.x[1] + m.x[2]) + sputils.attach_root_node(m, m.fsc, [m.x]) + m._mpisppy_probability = "uniform" + return m + + +def _integer_not_binary_scenario(sname, **kwargs): + """Integer vars with bounds (0, 5) — not treatable as binary.""" + m = pyo.ConcreteModel(name=sname) + m.x = pyo.Var([0, 1], domain=pyo.Integers, bounds=(0, 5)) + m.fsc = pyo.Expression(expr=m.x[0] + m.x[1]) + m.obj = pyo.Objective(expr=m.x[0] + m.x[1]) + sputils.attach_root_node(m, m.fsc, [m.x]) + m._mpisppy_probability = "uniform" + return m + + +def _integer_0_1_bounded_scenario(sname, **kwargs): + """Integer vars with bounds (0, 1) — the extension must accept these.""" + m = pyo.ConcreteModel(name=sname) + m.x = pyo.Var([0, 1, 2], domain=pyo.Integers, bounds=(0, 1)) + m.fsc = pyo.Expression(expr=m.x[0] + m.x[1] + m.x[2]) + m.obj = pyo.Objective(expr=m.x[0] + m.x[1] + m.x[2]) + sputils.attach_root_node(m, m.fsc, [m.x]) + m._mpisppy_probability = "uniform" + return m + + +def _make_sp(creator, num=2, cuts_count=3): + """Build a tiny SPBase around one of the scenario creators above.""" + return SPBase( + options={ + "toc": False, + "verbose": False, + "xhat_feasibility_cuts_count": cuts_count, + }, + all_scenario_names=[f"scen{i}" for i in range(num)], + scenario_creator=creator, + ) + + +class _FakeSolverPlugin: + """Stands in for a Pyomo SolverFactory instance. + + ``sputils.is_persistent`` checks ``isinstance(p, PersistentSolver)``; + this stub is neither, so the extension's persistent-solver branch + is correctly skipped. + """ + pass + + +def _attach_fake_solver_plugin(sp): + for s in sp.local_scenarios.values(): + s._solver_plugin = _FakeSolverPlugin() + + +class _Opt: + """Minimal ``opt`` surface that ``XhatFeasibilityCutExtension`` uses.""" + + def __init__(self, sp, spcomm=None): + self.local_scenarios = sp.local_scenarios + self.options = sp.options + self.spcomm = spcomm + self.multistage = sp.multistage + + +# --------------------------------------------------------------------------- +# Startup check +# --------------------------------------------------------------------------- + + +class TestStartupCheck(unittest.TestCase): + + def test_passes_on_binary_nonants(self): + sp = _make_sp(_binary_two_stage_scenario) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + ext.setup_hub() # must not raise + # ConstraintList should be attached on every scenario. + for s in sp.local_scenarios.values(): + self.assertTrue(hasattr(s._mpisppy_model, "xhat_feasibility_cuts")) + + def test_passes_on_integer_0_1_bounded(self): + """Integer vars with bounds (0, 1) are semantically binary.""" + sp = _make_sp(_integer_0_1_bounded_scenario) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + ext.setup_hub() # must not raise + + def test_raises_on_continuous_nonant(self): + sp = _make_sp(_continuous_two_stage_scenario) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + with self.assertRaises(RuntimeError) as cm: + ext.setup_hub() + msg = str(cm.exception) + self.assertIn("binary", msg) + self.assertIn("xhat-feasibility-cuts-count", msg) + + def test_raises_on_integer_not_0_1(self): + sp = _make_sp(_integer_not_binary_scenario) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + with self.assertRaises(RuntimeError): + ext.setup_hub() + + def test_init_rejects_zero_cuts_count(self): + sp = _make_sp(_binary_two_stage_scenario, cuts_count=0) + with self.assertRaises(RuntimeError) as cm: + XhatFeasibilityCutExtension(_Opt(sp)) + self.assertIn("xhat_feasibility_cuts_count", str(cm.exception)) + + +# --------------------------------------------------------------------------- +# Cut installation +# --------------------------------------------------------------------------- + + +class TestCutInstallation(unittest.TestCase): + + def _build(self, num=2, cuts_count=3): + sp = _make_sp(_binary_two_stage_scenario, num=num, cuts_count=cuts_count) + _attach_fake_solver_plugin(sp) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + ext.setup_hub() + return sp, ext + + def _row_len(self, sp): + any_s = next(iter(sp.local_scenarios.values())) + return 1 + len(any_s._mpisppy_data.nonant_indices) + + def _buf_for(self, sp, cuts_count=3): + return np.zeros(cuts_count * self._row_len(sp) + 1) + + def test_installs_one_cut_on_every_scenario(self): + sp, ext = self._build(num=3) + buf = self._buf_for(sp) + # Install a single cut: x[0] + x[1] + x[2] >= 1, packed as + # rhs_constant=-1, coefs=[1, 1, 1]. + buf[0] = -1.0 + buf[1:4] = [1.0, 1.0, 1.0] + buf[-1] = 1.0 # count + ext._install_cuts(buf) + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 1) + + def test_zero_count_is_noop(self): + sp, ext = self._build() + buf = self._buf_for(sp) # count stays 0 + ext._install_cuts(buf) + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 0) + + def test_all_zero_row_is_skipped_even_if_counted(self): + """Count header claims 2 cuts but both rows are all zeros.""" + sp, ext = self._build() + buf = self._buf_for(sp) + buf[-1] = 2.0 + ext._install_cuts(buf) + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 0) + + def test_two_successive_install_calls_accumulate(self): + sp, ext = self._build(num=2) + buf = self._buf_for(sp) + buf[0] = -1.0 + buf[1:4] = [1.0, 1.0, 1.0] + buf[-1] = 1.0 + ext._install_cuts(buf) + ext._install_cuts(buf) + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 2) + + +# --------------------------------------------------------------------------- +# No-good cut encoding (spoke side) +# --------------------------------------------------------------------------- + + +class _FakeSendArray: + def __init__(self, length): + self._arr = np.zeros(length) + + def value_array(self): + return self._arr + + +class _FakeSpcomm: + """Just enough of SPCommunicator for ``_maybe_emit_feasibility_cut``.""" + + def __init__(self, buffer_length): + self.send_buffers = { + Field.XHAT_FEASIBILITY_CUT: _FakeSendArray(buffer_length), + } + self.sent = [] # (buf_snapshot, field) for every put_send_buffer call + + def is_send_field_registered(self, field): + return field in self.send_buffers + + def put_send_buffer(self, send_array, field): + # Snapshot the buffer contents to allow assertions after clearing. + self.sent.append((send_array.value_array().copy(), field)) + + +class _FakeOpt: + """Stand-in for Xhat_Eval in emission tests.""" + + def __init__(self, sp, cuts_count=3): + self.local_scenarios = sp.local_scenarios + self.options = dict(sp.options) + self.options["xhat_feasibility_cuts_count"] = cuts_count + # buffer layout matches FieldLengths' formula + any_s = next(iter(sp.local_scenarios.values())) + nonant_len = len(any_s._mpisppy_data.nonant_indices) + self.spcomm = _FakeSpcomm(cuts_count * (nonant_len + 1) + 1) + # Satisfy attributes XhatBase.__init__ expects. + self.cylinder_rank = 0 + self.n_proc = 1 + self.scenario_names_to_rank = {"ROOT": { + sname: 0 for sname in sp.local_scenarios + }} + + +def _make_xhatbase(sp, cuts_count=3): + """Construct an XhatBase extension with a fake opt for emission tests. + + We bypass ``__init__`` (which wants a real opt) and set attributes + directly. + """ + ext = xhatbase_extension.XhatBase.__new__(xhatbase_extension.XhatBase) + ext.opt = _FakeOpt(sp, cuts_count=cuts_count) + ext.cylinder_rank = ext.opt.cylinder_rank + ext.n_proc = ext.opt.n_proc + ext.verbose = False + ext.scenario_name_to_rank = ext.opt.scenario_names_to_rank + return ext + + +def _set_xhat(sp, xhat_values): + """Write xhat values into every local scenario's nonant vars.""" + for s in sp.local_scenarios.values(): + for v, val in zip(s._mpisppy_data.nonant_indices.values(), xhat_values): + v.set_value(val) + + +class TestNoGoodCutEncoding(unittest.TestCase): + + def _expected_encoding(self, xhat): + """Reference encoding: constant = sum(xhat) - 1; coef_i = 1 - 2 xhat_i.""" + xhat_ints = [1 if v > 0.5 else 0 for v in xhat] + rhs_constant = sum(xhat_ints) - 1 + coefs = [1 - 2 * xi for xi in xhat_ints] + return rhs_constant, coefs + + def _emit_and_get_row(self, sp, xhat_values, cuts_count=3): + ext = _make_xhatbase(sp, cuts_count=cuts_count) + _set_xhat(sp, xhat_values) + ext._maybe_emit_feasibility_cut() + self.assertEqual(len(ext.opt.spcomm.sent), 1, + "expected exactly one put_send_buffer call") + buf, field = ext.opt.spcomm.sent[0] + self.assertEqual(field, Field.XHAT_FEASIBILITY_CUT) + nonant_len = len(xhat_values) + row_len = 1 + nonant_len + self.assertEqual(buf.shape[0], cuts_count * row_len + 1) + # Count slot must be 1. + self.assertEqual(buf[-1], 1.0) + # Rows 1..cuts_count-1 must be zero. + for k in range(1, cuts_count): + self.assertTrue( + np.all(buf[k * row_len:(k + 1) * row_len] == 0.0), + f"row {k} should be zero-padded", + ) + return buf[:row_len] + + def test_encoding_all_zero_xhat(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + row = self._emit_and_get_row(sp, [0, 0, 0]) + exp_const, exp_coefs = self._expected_encoding([0, 0, 0]) + self.assertEqual(row[0], exp_const) + np.testing.assert_array_equal(row[1:], exp_coefs) + + def test_encoding_all_one_xhat(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + row = self._emit_and_get_row(sp, [1, 1, 1]) + exp_const, exp_coefs = self._expected_encoding([1, 1, 1]) + self.assertEqual(row[0], exp_const) + np.testing.assert_array_equal(row[1:], exp_coefs) + + def test_encoding_mixed_xhat(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + row = self._emit_and_get_row(sp, [1, 0, 1]) + exp_const, exp_coefs = self._expected_encoding([1, 0, 1]) + self.assertEqual(row[0], exp_const) + np.testing.assert_array_equal(row[1:], exp_coefs) + + def test_cut_excludes_its_own_xhat(self): + """Sanity: the emitted cut evaluated at xhat must violate the + inequality constant + sum coef_i x_i >= 0.""" + sp = _make_sp(_binary_two_stage_scenario, num=1) + xhat = [1, 0, 1] + row = self._emit_and_get_row(sp, xhat) + lhs = row[0] + sum(c * x for c, x in zip(row[1:], xhat)) + self.assertLess(lhs, 0.0, + "no-good cut should violate at its own xhat") + + def test_emit_is_noop_when_feature_off(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + ext = _make_xhatbase(sp, cuts_count=0) + _set_xhat(sp, [1, 0, 1]) + ext._maybe_emit_feasibility_cut() + self.assertEqual(ext.opt.spcomm.sent, []) + + +# --------------------------------------------------------------------------- +# Multi-node nonants in the startup check +# --------------------------------------------------------------------------- + + +def _multi_node_binary_scenario(sname, **kwargs): + """A hand-built scenario with nonants at two different nodes. + + We do NOT wrap this in SPBase (which would insist on proper + branching-factor plumbing). Instead we construct just the + attributes ``XhatFeasibilityCutExtension._assert_all_nonants_binary`` + needs: ``_mpisppy_data.nonant_indices`` keyed by ``(node_name, i)``. + """ + m = pyo.ConcreteModel(name=sname) + m.x_root = pyo.Var([0, 1], domain=pyo.Binary) + m.x_stage2 = pyo.Var([0], domain=pyo.Binary) + m.obj = pyo.Objective(expr=m.x_root[0] + m.x_root[1] + m.x_stage2[0]) + m._mpisppy_node_list = [ + stree.ScenarioNode("ROOT", cond_prob=1.0, stage=1, + cost_expression=m.x_root[0], + nonant_list=[m.x_root], scen_model=m), + stree.ScenarioNode("ROOT_0", cond_prob=1.0, stage=2, + cost_expression=m.x_stage2[0], + nonant_list=[m.x_stage2], + parent_name="ROOT", scen_model=m), + ] + m._mpisppy_probability = 1.0 + return m + + +class _MultiNodeStub: + """Fake opt wrapping a multi-node scenario without SPBase.""" + + def __init__(self, scenario, options): + # Mimic the attributes _assert_all_nonants_binary needs. + scenario._mpisppy_data = _FakeData() + scenario._mpisppy_data.nonant_indices = {} + for node in scenario._mpisppy_node_list: + for i, v in enumerate(node.nonant_vardata_list): + scenario._mpisppy_data.nonant_indices[(node.name, i)] = v + self.local_scenarios = {"s0": scenario} + self.options = options + self.spcomm = None + self.multistage = True + + +class _FakeData: + pass + + +class TestMultiStageRejection(unittest.TestCase): + """V1 of the feature is two-stage only. + + The no-good cut row encodes coefficients positionally against each + scenario's ``nonant_indices``; in multi-stage, scenarios on + different branches have different per-stage-2+ variables at the + deeper indices, so installing the same row on every scenario lands + coefficients on unrelated variables. ``setup_hub`` must hard-fail + rather than silently install incorrect cuts. See + ``doc/designs/xhat_feasibility_cuts_design.md``. + """ + + def test_assert_two_stage_raises_when_multistage(self): + scen = _multi_node_binary_scenario("s0") + opt = _MultiNodeStub(scen, {"xhat_feasibility_cuts_count": 1}) + ext = XhatFeasibilityCutExtension(opt) + with self.assertRaises(RuntimeError) as cm: + ext._assert_two_stage() + self.assertIn("two-stage only", str(cm.exception)) + + def test_setup_hub_raises_before_binary_check(self): + """Two-stage gate runs before the binary check, so a multi-stage + model with all-binary nonants still hits the multi-stage error.""" + scen = _multi_node_binary_scenario("s0") + opt = _MultiNodeStub(scen, {"xhat_feasibility_cuts_count": 1}) + ext = XhatFeasibilityCutExtension(opt) + with self.assertRaises(RuntimeError) as cm: + ext._assert_two_stage() + ext._assert_all_nonants_binary() + self.assertIn("two-stage only", str(cm.exception)) + + +# --------------------------------------------------------------------------- +# Extension plumbing: register / sync / persistent-solver branch +# --------------------------------------------------------------------------- + + +class _SpcommStub: + """Lightweight ``spcomm`` stand-in for register_receive_fields / + sync_with_spokes coverage.""" + + def __init__(self, fields_to_ranks=None, recv_buffer=None): + self.fields_to_ranks = fields_to_ranks or {} + self._recv_buffer_by_rank = {} + if recv_buffer is not None: + # Seed whichever rank shows up. + ranks = self.fields_to_ranks.get(Field.XHAT_FEASIBILITY_CUT, [0]) + for r in ranks: + self._recv_buffer_by_rank[r] = recv_buffer + + def register_recv_field(self, field, rank): + return self._recv_buffer_by_rank.setdefault(rank, _RecvBufferStub()) + + +class _RecvBufferStub: + def __init__(self, arr=None, is_new_result=True): + self._arr = arr if arr is not None else np.zeros(1) + self._is_new = is_new_result + + def is_new(self): + return self._is_new + + def array(self): + return self._arr + + +class _PersistentSolverStub: + """Stub that passes ``sputils.is_persistent`` via isinstance. + + ``sputils.is_persistent`` checks ``isinstance(p, PersistentSolver)``, + so we inherit from the real abstract base. + """ + pass + + +class TestExtensionRegistrationAndSync(unittest.TestCase): + def _make_ext(self, cuts_count=3): + sp = _make_sp(_binary_two_stage_scenario, num=2, cuts_count=cuts_count) + _attach_fake_solver_plugin(sp) + ext = XhatFeasibilityCutExtension(_Opt(sp)) + ext.setup_hub() + return sp, ext + + def test_register_send_fields_is_noop(self): + _, ext = self._make_ext() + # Should return without raising; exercise line 116. + ext.register_send_fields() + + def test_register_receive_fields_no_emitting_spoke(self): + """If no rank advertises the cut field, recv_buffer stays None.""" + sp, ext = self._make_ext() + ext.opt.spcomm = _SpcommStub(fields_to_ranks={}) + ext.register_receive_fields() + self.assertIsNone(ext._recv_buffer) + + def test_register_receive_fields_with_one_emitter(self): + sp, ext = self._make_ext() + recv = _RecvBufferStub() + ext.opt.spcomm = _SpcommStub( + fields_to_ranks={Field.XHAT_FEASIBILITY_CUT: [2]}, + recv_buffer=recv, + ) + ext.register_receive_fields() + self.assertIs(ext._recv_buffer, recv) + + def test_register_receive_fields_multiple_emitters_asserts(self): + sp, ext = self._make_ext() + ext.opt.spcomm = _SpcommStub( + fields_to_ranks={Field.XHAT_FEASIBILITY_CUT: [1, 2]}, + ) + with self.assertRaises(AssertionError): + ext.register_receive_fields() + + def test_sync_with_spokes_no_buffer(self): + sp, ext = self._make_ext() + ext._recv_buffer = None + # No side effects, no raise. + ext.sync_with_spokes() + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 0) + + def test_sync_with_spokes_not_new_is_noop(self): + sp, ext = self._make_ext() + # A buffer that claims "not new" should be ignored. + nonant_len = 3 + buf = np.zeros(3 * (nonant_len + 1) + 1) + buf[0] = -1.0 + buf[1:4] = [1.0, 1.0, 1.0] + buf[-1] = 1.0 + ext._recv_buffer = _RecvBufferStub(arr=buf, is_new_result=False) + ext.sync_with_spokes() + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 0) + + def test_sync_with_spokes_new_installs_cut(self): + sp, ext = self._make_ext() + nonant_len = 3 + buf = np.zeros(3 * (nonant_len + 1) + 1) + buf[0] = -1.0 + buf[1:4] = [1.0, 1.0, 1.0] + buf[-1] = 1.0 + ext._recv_buffer = _RecvBufferStub(arr=buf, is_new_result=True) + ext.sync_with_spokes() + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._mpisppy_model.xhat_feasibility_cuts), 1) + + +class TestCutInstallPersistentSolverBranch(unittest.TestCase): + """Exercise the persistent-solver branch in ``_install_cuts``.""" + + def test_persistent_solver_gets_add_constraint(self): + # Build a PersistentSolver subclass so sputils.is_persistent says True. + from pyomo.solvers.plugins.solvers.persistent_solver import ( + PersistentSolver, + ) + + class _PersistentPlugin(PersistentSolver): + def __init__(self): + # Bypass PersistentSolver.__init__ which wants a config; we + # only need an isinstance-compatible object with add_constraint. + self.added = [] + + def add_constraint(self, con): + self.added.append(con) + + # Abstracts we don't actually exercise in this test. + def _presolve(self, *a, **kw): pass + def _postsolve(self, *a, **kw): pass + def _warm_start(self, *a, **kw): pass + def _apply_solver(self, *a, **kw): pass + def _get_dual_values(self, *a, **kw): return {} + def _load_vars(self, *a, **kw): pass + def _get_expr_from_pyomo_repn(self, *a, **kw): pass + + sp = _make_sp(_binary_two_stage_scenario, num=2, cuts_count=3) + for s in sp.local_scenarios.values(): + s._solver_plugin = _PersistentPlugin() + ext = XhatFeasibilityCutExtension(_Opt(sp)) + ext.setup_hub() + + nonant_len = 3 + buf = np.zeros(3 * (nonant_len + 1) + 1) + buf[0] = -1.0 + buf[1:4] = [1.0, 1.0, 1.0] + buf[-1] = 1.0 + ext._install_cuts(buf) + + # Each local scenario's persistent plugin should have been told + # about the new constraint. + for s in sp.local_scenarios.values(): + self.assertEqual(len(s._solver_plugin.added), 1) + + +# --------------------------------------------------------------------------- +# _maybe_emit_feasibility_cut: early-return branches and buffer-size check +# --------------------------------------------------------------------------- + + +class TestEmitEarlyReturns(unittest.TestCase): + def _ext(self, sp, cuts_count=3, spcomm=None): + ext = xhatbase_extension.XhatBase.__new__(xhatbase_extension.XhatBase) + ext.opt = _FakeOpt(sp, cuts_count=cuts_count) + if spcomm is not None: + ext.opt.spcomm = spcomm + ext.cylinder_rank = 0 + ext.n_proc = 1 + ext.verbose = False + ext.scenario_name_to_rank = ext.opt.scenario_names_to_rank + return ext + + def test_spcomm_is_none_silently_returns(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + ext = self._ext(sp, cuts_count=3, spcomm="replace") + ext.opt.spcomm = None # override the _FakeSpcomm installed by _FakeOpt + _set_xhat(sp, [1, 0, 1]) + # Must not raise. + ext._maybe_emit_feasibility_cut() + + def test_send_field_not_registered_silently_returns(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + ext = self._ext(sp, cuts_count=3) + # Remove the pre-registered buffer so is_send_field_registered is False. + ext.opt.spcomm.send_buffers.pop(Field.XHAT_FEASIBILITY_CUT) + _set_xhat(sp, [1, 0, 1]) + ext._maybe_emit_feasibility_cut() + self.assertEqual(ext.opt.spcomm.sent, []) + + def test_buffer_size_mismatch_raises(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + # cuts_count=3 in options → emitter expects 3*(3+1)+1 = 13. + # Replace the buffer with a wrong size. + ext = self._ext(sp, cuts_count=3) + ext.opt.spcomm.send_buffers[Field.XHAT_FEASIBILITY_CUT] = \ + _FakeSendArray(5) # wrong size + _set_xhat(sp, [1, 0, 1]) + with self.assertRaises(RuntimeError) as cm: + ext._maybe_emit_feasibility_cut() + self.assertIn("buffer has length 5", str(cm.exception)) + + def test_missing_nonant_value_silently_returns(self): + sp = _make_sp(_binary_two_stage_scenario, num=1) + ext = self._ext(sp, cuts_count=3) + # Clear one nonant's value to hit the `xv is None` guard. + for s in sp.local_scenarios.values(): + first_v = next(iter(s._mpisppy_data.nonant_indices.values())) + first_v.set_value(None, skip_validation=True) + break + ext._maybe_emit_feasibility_cut() + self.assertEqual(ext.opt.spcomm.sent, []) + + +# --------------------------------------------------------------------------- +# Plumbing tests: config registration, cfg_vanilla helpers, FieldLengths, +# and the _try_one exception wrapper. +# --------------------------------------------------------------------------- + + +class TestConfigArgRegistration(unittest.TestCase): + def test_xhat_feasibility_cut_args_registers_with_default_zero(self): + from mpisppy.utils import config as cfgmod + cfg = cfgmod.Config() + cfg.xhat_feasibility_cut_args() + self.assertIn("xhat_feasibility_cuts_count", cfg) + self.assertEqual(cfg.get("xhat_feasibility_cuts_count"), 0) + + def test_generic_parsing_registers_the_flag(self): + """Covers the cfg.xhat_feasibility_cut_args() call inside + mpisppy.generic.parsing.parse_args.""" + import sys + import types + import mpisppy.generic.parsing as parsing + + stub = types.ModuleType("__stub_model_for_parse_args__") + def inparser_adder(cfg): + cfg.add_to_config("num_scens", "stub", int, default=1) + stub.inparser_adder = inparser_adder + + saved_argv = sys.argv + sys.argv = ["prog"] # no positional flags; defaults only + try: + cfg = parsing.parse_args(stub) + finally: + sys.argv = saved_argv + + # The flag made it into the Config — which is the signal that + # parse_args actually ran cfg.xhat_feasibility_cut_args(). + self.assertIn("xhat_feasibility_cuts_count", cfg) + + +class TestCfgVanillaPlumbing(unittest.TestCase): + def _cfg(self, xhat_cap=0): + """Build a minimal Config with whatever fields shared_options needs.""" + from mpisppy.utils import config as cfgmod + cfg = cfgmod.Config() + cfg.popular_args() + cfg.two_sided_args() + cfg.ph_args() + cfg.aph_args() + cfg.xhat_feasibility_cut_args() + if xhat_cap: + cfg.xhat_feasibility_cuts_count = xhat_cap + cfg.solver_name = "gurobi" + return cfg + + def test_shared_options_propagates_cap(self): + from mpisppy.utils import cfg_vanilla as vanilla + cfg = self._cfg(xhat_cap=4) + shopts = vanilla.shared_options(cfg) + self.assertEqual(shopts["xhat_feasibility_cuts_count"], 4) + + def test_shared_options_defaults_to_zero_when_off(self): + from mpisppy.utils import cfg_vanilla as vanilla + cfg = self._cfg(xhat_cap=0) + shopts = vanilla.shared_options(cfg) + self.assertEqual(shopts["xhat_feasibility_cuts_count"], 0) + + def test_add_xhat_feasibility_cuts_noop_when_off(self): + from mpisppy.utils import cfg_vanilla as vanilla + cfg = self._cfg(xhat_cap=0) + hub_dict = {"opt_kwargs": {"options": {}}} + result = vanilla.add_xhat_feasibility_cuts(hub_dict, cfg) + # Same dict back; no extensions attached; no option injected. + self.assertIs(result, hub_dict) + self.assertNotIn("extensions", hub_dict["opt_kwargs"]) + self.assertNotIn("xhat_feasibility_cuts_count", + hub_dict["opt_kwargs"]["options"]) + + def test_add_xhat_feasibility_cuts_attaches_extension_when_on(self): + from mpisppy.utils import cfg_vanilla as vanilla + from mpisppy.extensions.xhat_feasibility_cut_extension import ( + XhatFeasibilityCutExtension, + ) + cfg = self._cfg(xhat_cap=7) + hub_dict = {"opt_kwargs": {"options": {}}} + vanilla.add_xhat_feasibility_cuts(hub_dict, cfg) + # extension_adder installs either a direct extension or a + # MultiExtension; confirm XhatFeasibilityCutExtension made it in. + exts = hub_dict["opt_kwargs"].get("extensions") + self.assertIsNotNone(exts) + ext_classes = ( + [exts] if not isinstance(exts, type) or exts.__name__ != "MultiExtension" + else hub_dict["opt_kwargs"]["extension_kwargs"]["ext_classes"] + ) + # extension_adder produces either a single class or a MultiExtension + # with ext_classes. Handle both. + mx_classes = hub_dict["opt_kwargs"].get("extension_kwargs", {}).get( + "ext_classes", []) + self.assertTrue( + XhatFeasibilityCutExtension in ext_classes + or XhatFeasibilityCutExtension in mx_classes, + f"Expected XhatFeasibilityCutExtension in hub extensions; got {exts!r}", + ) + self.assertEqual( + hub_dict["opt_kwargs"]["options"]["xhat_feasibility_cuts_count"], + 7, + ) + + +class TestFieldLengthsPicksUpCap(unittest.TestCase): + """Covers the FieldLengths setter line that reads + ``opt.options['xhat_feasibility_cuts_count']``.""" + + def test_field_length_reflects_cap(self): + from mpisppy.cylinders.spwindow import FieldLengths, Field + sp = _make_sp(_binary_two_stage_scenario, num=2, cuts_count=5) + # FieldLengths wants opt.nonant_length; SPBase sets that up. + fl = FieldLengths(sp) + # row_len = 1 + 3 nonants; cap=5 → 5*4 + 1 = 21. + self.assertEqual(fl[Field.XHAT_FEASIBILITY_CUT], 5 * 4 + 1) + + def test_field_length_off_is_one(self): + from mpisppy.cylinders.spwindow import FieldLengths, Field + sp = _make_sp(_binary_two_stage_scenario, num=2, cuts_count=0) + fl = FieldLengths(sp) + # With cap=0 the buffer collapses to just the trailing count slot. + self.assertEqual(fl[Field.XHAT_FEASIBILITY_CUT], 1) + + +class TestGenericExtensionsWiring(unittest.TestCase): + """Covers the one-line gate in mpisppy/generic/extensions.py that + defers to cfg_vanilla.add_xhat_feasibility_cuts when the flag is + positive. We minimally construct the cfg the function reads.""" + + def _make_configured_cfg(self, cuts_count): + from mpisppy.utils import config as cfgmod + cfg = cfgmod.Config() + # configure_extensions reads these; register them as no-ops. + cfg.popular_args() + cfg.ph_args() + cfg.aph_args() + cfg.two_sided_args() + cfg.fixer_args() + cfg.relaxed_ph_fixer_args() + cfg.integer_relax_then_enforce_args() + cfg.gapper_args() + cfg.gapper_args(name="lagrangian") + cfg.ph_primal_args() + cfg.ph_dual_args() + cfg.relaxed_ph_args() + cfg.gradient_args() + cfg.dynamic_rho_args() + cfg.reduced_costs_args() + cfg.sep_rho_args() + cfg.coeff_rho_args() + cfg.sensi_rho_args() + cfg.reduced_costs_rho_args() + cfg.norm_rho_args() + cfg.primal_dual_rho_args() + cfg.wxbar_read_write_args() + cfg.tracking_args() + cfg.converger_args() + cfg.xhat_feasibility_cut_args() + cfg.add_to_config("user_defined_extensions", + description="list of user-defined extensions", + domain=list, default=None) + cfg.add_to_config("write_scenario_lp_mps_files_dir", + description="not used here", + domain=str, default=None) + cfg.solver_name = "gurobi" + cfg.xhat_feasibility_cuts_count = cuts_count + return cfg + + def _fake_hub_dict(self): + return {"opt_kwargs": {"options": {}}} + + def test_configure_extensions_attaches_when_on(self): + from mpisppy.generic.extensions import configure_extensions + from mpisppy.extensions.xhat_feasibility_cut_extension import ( + XhatFeasibilityCutExtension, + ) + cfg = self._make_configured_cfg(cuts_count=3) + hub = self._fake_hub_dict() + configure_extensions(hub, module=None, cfg=cfg) + exts = hub["opt_kwargs"].get("extensions") + mx_classes = hub["opt_kwargs"].get("extension_kwargs", {}).get( + "ext_classes", []) + self.assertTrue( + exts is XhatFeasibilityCutExtension + or XhatFeasibilityCutExtension in mx_classes, + f"Expected XhatFeasibilityCutExtension in hub; got exts={exts!r}, " + f"ext_classes={mx_classes!r}", + ) + + def test_configure_extensions_skipped_when_off(self): + from mpisppy.generic.extensions import configure_extensions + from mpisppy.extensions.xhat_feasibility_cut_extension import ( + XhatFeasibilityCutExtension, + ) + cfg = self._make_configured_cfg(cuts_count=0) + hub = self._fake_hub_dict() + configure_extensions(hub, module=None, cfg=cfg) + exts = hub["opt_kwargs"].get("extensions") + mx_classes = hub["opt_kwargs"].get("extension_kwargs", {}).get( + "ext_classes", []) + self.assertFalse( + exts is XhatFeasibilityCutExtension + or XhatFeasibilityCutExtension in mx_classes, + "XhatFeasibilityCutExtension should NOT be attached when off", + ) + + +class TestTryOneExceptionWrapper(unittest.TestCase): + """Covers the ``try / except`` wrapper around + ``_maybe_emit_feasibility_cut`` inside ``XhatBase._try_one``. + + We don't drive ``_try_one`` end-to-end (that would need the whole + Xhat_Eval machinery); we exercise the tiny wrapper by calling the + same snippet directly. The wrapper is intentionally simple: log on + rank 0 and swallow the exception so the xhatter keeps going. + """ + + def test_try_one_reaches_the_emit_wrapper_on_infeasibility(self): + """Drive ``_try_one`` just far enough to hit the emit wrapper. + + Stubs ``self.opt`` with the minimum set of methods ``_try_one`` + touches on the two-stage infeasibility path, and replaces + ``_maybe_emit_feasibility_cut`` with one that raises. The + wrapper must log-and-swallow (no re-raise) so the xhatter + keeps going. + """ + import mpisppy.extensions.xhatbase as xb + sp = _make_sp(_binary_two_stage_scenario, num=1) + ext = xb.XhatBase.__new__(xb.XhatBase) + ext.cylinder_rank = 0 + ext.n_proc = 1 + ext.verbose = False + ext.scenario_name_to_rank = {"ROOT": {"scen0": 0}} + + class _StubOpt: + def __init__(self, sp): + self.local_scenarios = sp.local_scenarios + self.options = {} + self._saved = 0 + self._restored = 0 + # inject a nonant_cache on the one local scenario to + # satisfy the two-stage branch in _try_one. + for s in sp.local_scenarios.values(): + s._mpisppy_data.nonant_cache = np.zeros( + len(s._mpisppy_data.nonant_indices)) + def _save_nonants(self): self._saved += 1 + def _restore_nonants(self): self._restored += 1 + def _fix_nonants(self, xhats): pass + def solve_loop(self, **kw): pass + def no_incumbent_prob(self): return 1.0 # force infeasibility + def Eobjective(self, **kw): return 0.0 + def update_best_solution_if_improving(self, obj): return False + + class _StubComm: + def bcast(self, data, root=0): return data + ext.opt = _StubOpt(sp) + ext.comms = {"ROOT": _StubComm()} + + # Monkey-patch the emit helper to raise so the wrapper's except + # branch fires. + def _boom(self): + raise RuntimeError("synthetic emit failure") + ext._maybe_emit_feasibility_cut = _boom.__get__(ext, xb.XhatBase) + + # Should NOT raise: _try_one catches and logs. + result = ext._try_one({"ROOT": "scen0"}, + solver_options=None, + verbose=False, + restore_nonants=True) + self.assertIsNone(result) + # nonants were restored despite the exception. + self.assertEqual(ext.opt._restored, 1) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/tests/test_xhat_from_file.py b/mpisppy/tests/test_xhat_from_file.py index 9403ff25c..53782d289 100644 --- a/mpisppy/tests/test_xhat_from_file.py +++ b/mpisppy/tests/test_xhat_from_file.py @@ -54,11 +54,13 @@ def _make_sp(num=2, options=None): class _FakeXhatEval: """Just enough ``Xhat_Eval`` surface for ``_try_file_xhat``.""" - def __init__(self, sp, evaluate_result): + def __init__(self, sp, evaluate_result, infeasP=0.0, spcomm=None): self.local_scenarios = sp.local_scenarios self.options = sp.options self.multistage = sp.multistage self._evaluate_result = evaluate_result + self._infeasP = infeasP + self.spcomm = spcomm self.evaluate_calls = [] self.restore_calls = 0 @@ -66,15 +68,18 @@ def evaluate(self, nonant_cache): self.evaluate_calls.append(nonant_cache) return self._evaluate_result + def no_incumbent_prob(self): + return self._infeasP + def _restore_nonants(self): self.restore_calls += 1 -def _make_helper(sp, evaluate_result=None): +def _make_helper(sp, evaluate_result=None, infeasP=0.0, spcomm=None): """Build an ``XhatInnerBoundBase``-shaped object without going through its __init__ (which needs the whole spoke stack).""" h = XhatInnerBoundBase.__new__(XhatInnerBoundBase) - h.opt = _FakeXhatEval(sp, evaluate_result) + h.opt = _FakeXhatEval(sp, evaluate_result, infeasP=infeasP, spcomm=spcomm) h.cylinder_rank = 0 h.updates = [] @@ -194,6 +199,126 @@ def test_none_eobj_does_not_update(self): os.remove(path) +class _FakeSendArray: + def __init__(self, length): + self._arr = np.zeros(length) + + def value_array(self): + return self._arr + + +class _FakeSpcomm: + """Just enough of SPCommunicator for ``pack_no_good_feasibility_cut``.""" + + def __init__(self, buffer_length): + from mpisppy.cylinders.spwindow import Field + self.send_buffers = { + Field.XHAT_FEASIBILITY_CUT: _FakeSendArray(buffer_length), + } + self.sent = [] + + def is_send_field_registered(self, field): + return field in self.send_buffers + + def put_send_buffer(self, send_array, field): + self.sent.append((send_array.value_array().copy(), field)) + + +def _set_xhat(sp, xhat_values): + """Write xhat values into every local scenario's nonant vars.""" + for s in sp.local_scenarios.values(): + for v, val in zip(s._mpisppy_data.nonant_indices.values(), xhat_values): + v.set_value(val) + + +class TestFileXhatEmitsFeasibilityCut(unittest.TestCase): + """End-to-end on the cylinder side: feeding an infeasible binary xhat + via ``--xhat-from-file`` plus ``--xhat-feasibility-cuts-count > 0`` + must populate the ``XHAT_FEASIBILITY_CUT`` send buffer with the + no-good row for that xhat and call ``put_send_buffer``. + + Stubs ``Xhat_Eval`` so no solver/MPI is needed; the spcomm side is + a simple capture, identical to the pattern in + ``test_xhat_feasibility_cuts.py``.""" + + def _setup(self, xhat_values, cuts_count=1): + path = _write_npy(xhat_values) + sp = _make_sp(options={ + "xhat_from_file": path, + "xhat_feasibility_cuts_count": cuts_count, + }) + # Fix nonants to the infeasible candidate so pack_no_good_* + # reads them back correctly. (Real evaluate would do this via + # _fix_nonants; the fake skips that, so we do it directly.) + _set_xhat(sp, xhat_values) + nonant_len = len(xhat_values) + spcomm = _FakeSpcomm(cuts_count * (nonant_len + 1) + 1) + # evaluate_result is irrelevant on the infeasible path; what + # drives cut emission is no_incumbent_prob() != 0. + helper = _make_helper(sp, evaluate_result=None, infeasP=0.5, + spcomm=spcomm) + return path, helper, spcomm + + def test_infeasible_file_xhat_emits_no_good_cut(self): + from mpisppy.cylinders.spwindow import Field + path, helper, spcomm = self._setup([1.0, 0.0, 1.0]) + try: + helper._try_file_xhat() + finally: + os.remove(path) + # One emission, on the right field. + self.assertEqual(len(spcomm.sent), 1) + buf, field = spcomm.sent[0] + self.assertEqual(field, Field.XHAT_FEASIBILITY_CUT) + # Row format for xhat=[1,0,1]: + # constant = (#ones) - 1 = 1 + # coefs = [1 - 2*1, 1 - 2*0, 1 - 2*1] = [-1, 1, -1] + # tail count slot = 1.0 + np.testing.assert_array_equal(buf[:4], [1.0, -1.0, 1.0, -1.0]) + self.assertEqual(buf[-1], 1.0) + # Cut excludes its own xhat: constant + sum(coef * xhat) < 0. + lhs = buf[0] + sum(c * x for c, x in zip(buf[1:4], [1.0, 0.0, 1.0])) + self.assertLess(lhs, 0.0) + # No inner-bound improvement (infeasible candidate). + self.assertEqual(helper.updates, []) + # Nonants restored. + self.assertEqual(helper.opt.restore_calls, 1) + + def test_feature_off_file_xhat_does_not_emit(self): + """Same setup but xhat_feasibility_cuts_count=0 — pack helper + early-returns and put_send_buffer is never called.""" + path, helper, spcomm = self._setup([1.0, 0.0, 1.0], cuts_count=0) + # cuts_count=0 means the buffer would be size 1; rebuild spcomm. + helper.opt.spcomm = _FakeSpcomm(1) + try: + helper._try_file_xhat() + finally: + os.remove(path) + self.assertEqual(helper.opt.spcomm.sent, []) + self.assertEqual(helper.updates, []) + self.assertEqual(helper.opt.restore_calls, 1) + + def test_feasible_file_xhat_does_not_emit(self): + """Eobj is finite and no_incumbent_prob is 0 — no cut emission, + and update_if_improving is called.""" + path = _write_npy([1.0, 0.0, 1.0]) + sp = _make_sp(options={ + "xhat_from_file": path, + "xhat_feasibility_cuts_count": 1, + }) + _set_xhat(sp, [1.0, 0.0, 1.0]) + spcomm = _FakeSpcomm(1 * (3 + 1) + 1) + helper = _make_helper(sp, evaluate_result=42.0, infeasP=0.0, + spcomm=spcomm) + try: + helper._try_file_xhat() + finally: + os.remove(path) + self.assertEqual(spcomm.sent, []) + self.assertEqual(helper.updates, [42.0]) + self.assertEqual(helper.opt.restore_calls, 1) + + class TestMathIsfinite(unittest.TestCase): """Sanity check that math.isfinite is the predicate we intend: inf is not finite, large real is finite.""" diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 4b76d5d25..fe769a956 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -163,6 +163,11 @@ def shared_options(cfg, is_hub=False): # Optional initial xhat candidate file (.npy); None disables. # Consumed by XhatInnerBoundBase._try_file_xhat. "xhat_from_file" : cfg.get("xhat_from_file", None), + # Cap for feasibility cuts emitted by xhat spokes. The hub + # extension reads this (when attached); each xhat spoke reads + # it to gate _maybe_emit_feasibility_cut; spwindow.FieldLengths + # reads it to size the shared-memory buffer. Default 0 disables. + "xhat_feasibility_cuts_count" : cfg.get("xhat_feasibility_cuts_count", 0) or 0, "inspect_buffers_on_shutdown" : cfg.get("inspect_buffers_on_shutdown", False), # Optional filename prefix; if set, _BoundSpoke.update_if_improving # writes a first-stage solution snapshot on each new best incumbent. @@ -796,6 +801,26 @@ def add_cross_scenario_cuts(hub_dict, = {"check_bound_improve_iterations" : cfg.cross_scenario_iter_cnt} return hub_dict + +def add_xhat_feasibility_cuts(hub_dict, cfg): + """Wire the hub-side receiver for xhat feasibility cuts. + + No-op if ``cfg.xhat_feasibility_cuts_count`` is zero. Otherwise + attaches ``XhatFeasibilityCutExtension`` to the hub; the hub + extension's ``setup_hub`` hard-fails if any first-stage nonant is + not binary. Propagates the cap into hub options so the extension + sees it. See ``doc/src/xhat_feasibility_cuts.rst``. + """ + n = cfg.get("xhat_feasibility_cuts_count", 0) or 0 + if n <= 0: + return hub_dict + from mpisppy.extensions.xhat_feasibility_cut_extension import ( + XhatFeasibilityCutExtension, + ) + hub_dict = extension_adder(hub_dict, XhatFeasibilityCutExtension) + hub_dict["opt_kwargs"]["options"]["xhat_feasibility_cuts_count"] = n + return hub_dict + def add_reduced_costs_fixer(hub_dict, cfg, ): diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 047be1d84..07dce2dd2 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -1134,6 +1134,19 @@ def xhat_from_file_args(self): domain=str, default=None) + def xhat_feasibility_cut_args(self): + # Optional feasibility cuts emitted by xhat spokes when a candidate + # xhat is infeasible in some scenario. Default 0 disables the + # feature; positive N caps cuts per iteration and sizes the send + # buffer. First release supports binary first-stage only; the + # hub extension hard-fails at setup if that precondition is not + # met. See doc/src/xhat_feasibility_cuts.rst. + self.add_to_config("xhat_feasibility_cuts_count", + description="max feasibility cuts per xhat iteration " + "(0 disables; >0 requires binary first-stage)", + domain=int, + default=0) + def wtracker_args(self): self.add_to_config('wtracker', diff --git a/run_coverage.bash b/run_coverage.bash index 65bdbeca1..302df036e 100755 --- a/run_coverage.bash +++ b/run_coverage.bash @@ -149,6 +149,9 @@ run_phase "test_flexible_rank_cli (serial)" \ run_phase "test_xhat_from_file (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_xhat_from_file.py -v +run_phase "test_xhat_feasibility_cuts (serial)" \ + coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_xhat_feasibility_cuts.py -v + run_phase "test_incumbent_writing (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_incumbent_writing.py -v