diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index b0804b50d..1b7c80608 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 feasible_xhat + run: | + coverage run $COV_ARGS -m pytest mpisppy/tests/test_feasible_xhat.py -v + - name: Test docs run: | cd ./doc/ @@ -1063,6 +1067,7 @@ jobs: mpisppy/tests/test_solver_spec.py \ mpisppy/tests/test_extensions.py \ mpisppy/tests/test_jensens.py \ + mpisppy/tests/test_feasible_xhat.py \ mpisppy/tests/test_proper_bundler.py \ mpisppy/tests/test_incumbent_writing.py diff --git a/doc/src/feasible_xhat.rst b/doc/src/feasible_xhat.rst new file mode 100644 index 000000000..940180db2 --- /dev/null +++ b/doc/src/feasible_xhat.rst @@ -0,0 +1,326 @@ +.. _feasible_xhat: + +Per-scenario-feasible candidate xhat +==================================== + +.. warning:: + + This is an advanced topic. Most ``mpi-sppy`` users do **not** need + ``feasible_xhat_creator`` to get started. However, if you know + of some way to get a good solution that is admissible and implementable, + this function could speed up the search for a good upper bound. + You can use any method to get this solution, but the documentation + here is biased toward a situation where you know how + to modify average scenario data so it will work as a way to get a good + xhat solution. There are some helper functions to facilitate that. + +The contract +------------ + +A scenario module that participates exposes: + +.. code-block:: python + + def feasible_xhat_creator(*, solver_name, + solver_options=None, + **scenario_creator_kwargs): + """Return {nodename: np.ndarray} -- a candidate first-stage + point that is feasible to fix in every real scenario's + per-scenario subproblem. Two-stage: only "ROOT" is populated.""" + +Discovery is via ``getattr(module, "feasible_xhat_creator", None)``, +parallel to ``average_scenario_creator`` (see :ref:`jensens`). The +returned dict is in the cache form consumed by +``Xhat_Eval._fix_nonants``; the ``np.ndarray`` is in +``_mpisppy_node_list[0].nonant_vardata_list`` order. + +Two-stage only, for now. Multi-stage extensions would need a candidate +per non-leaf node, plus a story for inter-stage feasibility coupling +that does not exist in the two-stage case. + +File layout +----------- + +By convention the model author's ``feasible_xhat_creator`` (and any +``average_scenario_creator``) lives in +``examples//_auxiliary.py``, **not** in the main +example file. The introductory example (``farmer.py``, ``netdes.py``, +...) stays focused on what a first-time reader needs: +``scenario_creator``, ``scenario_names_creator``, ``kw_creator``, +``inparser_adder``, ``sample_tree_scen_creator``, +``scenario_denouement``. Auxiliary functions used only by advanced +machinery go in the ``_auxiliary`` sibling. + +.. admonition:: Discovery + :class: note + + When one of the xhat-spoke flags below is set, + ``cfg_vanilla._find_feasible_xhat_creator`` first tries + ``getattr(scenario_module, "feasible_xhat_creator", None)`` on the + user's main scenario module; if that is ``None`` it imports + ``_auxiliary`` and looks there. Discovery is gated on + the flag, so the auxiliary import does not happen unless requested. + Downstream consumers (e.g. findW) that bypass the cylinder system + import the auxiliary module directly and call + ``feasible_xhat_creator`` themselves. + +Why this convention exists +-------------------------- + +The Jensen's xhat path (see :ref:`jensens`) already exposes the +average-scenario solution as a candidate first-stage. It tolerates +infeasibility: if the candidate xhat value, when fixed for one +or more scenarios is infeasible, then +``_jensens_evaluate_xhat`` returns ``None`` and the inner-bound spoke +silently moves on. That is the right behavior for an inner-bound +spoke whose only job is to opportunistically improve a bound. + +``feasible_xhat_creator`` strengthens the contract. Where Jensen's +xhat is opportunistic -- silently skipping any scenario in which the +candidate is not feasible to fix -- ``feasible_xhat_creator`` tries +to guarantee by construction that the candidate **is** feasible to +fix in every real scenario. The inner-bound spoke that consumes it +therefore produces a usable objective on each candidate +rather than potentially never updating the inner bound. + +That is a strictly stronger contract than "Jensen's xhat plus luck," +and it is the contract that ``feasible_xhat_creator`` aims to +provide. + +In-cylinder use: ``---try-feasible-xhat-first`` flags +-------------------------------------------------------------- + +The four xhat spokes that ship with ``mpi-sppy`` accept a +``feasible_xhat_creator`` candidate via a per-spoke flag, parallel to +``--*-try-jensens-first``: + +* ``--xhatshuffle-try-feasible-xhat-first`` +* ``--xhatxbar-try-feasible-xhat-first`` +* ``--xhatlooper-try-feasible-xhat-first`` +* ``--xhatspecific-try-feasible-xhat-first`` + +When set, the spoke calls the module's ``feasible_xhat_creator`` once +before entering its main loop, fixes the candidate as the first-stage +nonants, evaluates the expected objective across all real scenarios, +and -- if the evaluation is feasible -- sends that as its first inner +bound. Implementation lives in +``_JensensMixin._try_feasible_xhat`` in +``mpisppy/cylinders/_jensens_mixin.py``; the spoke ``main()`` methods +call it once after ``_try_average_scenario_xhat``. + +.. admonition:: Mutually exclusive with ``--*-try-jensens-first`` + :class: warning + + ``---try-jensens-first`` and + ``---try-feasible-xhat-first`` are mutually exclusive on + the same spoke. ``cfg_vanilla._maybe_attach_feasible_xhat`` raises + at spoke-setup time if both are enabled, with a message naming the + conflicting CLI options. + + The two pre-loop candidates serve overlapping purposes: Jensen's + often gives a tighter incumbent bound when its candidate happens to + be feasible everywhere, while ``feasible_xhat_creator`` is + guaranteed feasible by contract but can be a looser incumbent. Per + spoke, pick whichever fits the model's structure -- not both. + + Across spokes, mixing is fine: one xhat spoke can be configured + with ``--xhatshuffle-try-jensens-first`` while another runs with + ``--xhatxbar-try-feasible-xhat-first``. + +Average-data-based Methods +-------------------------- + +The two lp-based helpers in ``mpisppy.utils.xhat_helpers`` +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Some ``feasible_xhat_creator`` implementations are short. ``mpi-sppy`` +ships two reusable, lp-based engines that try to do the heavy lifting; the +implementations call one of them and apply a model-specific repair. + +``average_xhat_nonants(average_scenario_creator, *, solver_name, ...)`` + Builds the model returned by ``average_scenario_creator``, + optionally LP-relaxes it, solves it, and returns the ROOT first-stage + values as ``np.ndarray``. One deterministic solve over the average + data. + +``lp_xbar_nonants(scenario_creator, scenario_names, *, solver_name, ...)`` + For each real scenario, builds the model, applies + ``core.relax_integer_vars``, solves, and returns the + probability-weighted average of ROOT first-stage values across + scenarios. ``K`` LP solves, where ``K`` is the number of scenarios. + +These two are **not interchangeable** for models with binary first-stage: +averaging data and averaging solutions do not commute when the optimal +first-stage is not a continuous function of the data. The averaged-data +problem can omit first-stage activity that some real scenario individually +needs; the per-scenario LP-xbar instead carries any activity that any +scenario's LP wanted positive into the average, where a feasibility- +preserving rounding rule can promote it. For models with continuous +first-stage, the distinction collapses. + +Choosing between the two engines is the caller's responsibility. The +caller knows whether averaging data preserves enough information to +cover per-scenario feasibility; the framework cannot detect that from +the model. + +The rounding rule is also yours +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The output of either engine is a real-valued vector that has to be +turned into a feasible candidate. Whether the right rule is +``np.ceil``, ``np.floor``, ``np.round``, identity, or a per-component +try-and-check is a model-specific decision that depends on +**monotonicity of recourse feasibility in each first-stage variable**: + +* If raising :math:`x_e` from 0 to 1 only loosens recourse + constraints (as for "open the arc" binaries in netdes), ``np.ceil`` + is feasibility-preserving. +* If the variable indexes a covering decision (open the facility) and + more open never tightens recourse, ``np.round`` typically suffices. +* If recourse feasibility is non-monotone in the variable, neither + rule is safe and the implementation must do something model-specific + (a proof-of-feasibility per-component repair, an aggregation across + scenarios, etc.). + +``mpi-sppy`` does **not** ship an automatic rounder. Even within a +single model, different first-stage variables can need different +rules; per-component try-and-check degenerates into solving an +MIP-feasibility problem in itself. The repair belongs in the +``feasible_xhat_creator``, where the model author has the domain +knowledge. + +Worked example: farmer (continuous first-stage) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Farmer's first-stage variable ``DEVOTED_ACRES`` is bounded +``NonNegativeReals``, and farmer has relatively complete recourse via +the buy/sell variables (``QuantityPurchased``, +``QuantitySubQuotaSold``, ``QuantitySuperQuotaSold``), so any feasible +acreage allocation -- including the average-scenario optimum -- is +feasible to fix in every real scenario. No rounding is needed. + +``examples/farmer/farmer_auxiliary.py``: + +.. code-block:: python + + from mpisppy.utils.xhat_helpers import average_xhat_nonants + from farmer import average_scenario_creator + + + def feasible_xhat_creator(*, solver_name, solver_options=None, + **scenario_creator_kwargs): + arr = average_xhat_nonants( + average_scenario_creator, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": arr} + +This is the simplest case the convention has to handle, and it +illustrates an important point about the convention: callers always +go through ``feasible_xhat_creator`` rather than calling +``average_xhat_nonants`` directly. If a downstream model swap replaces +farmer with a binary-first-stage model, only the auxiliary file has +to change; the call site at the consumer (e.g., findW) is unchanged. + +Worked example: netdes (binary, arc-open monotonicity) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Netdes ``model.x[e]`` is ``Binary`` for each candidate arc. The +recourse constraint is :math:`y_e \le u_e \, x_e`; raising :math:`x_e` +from 0 to 1 only loosens this bound, and the flow-balance constraints +do not involve ``x``. So opening more arcs cannot make any per- +scenario subproblem less feasible -- ``np.ceil`` is feasibility- +preserving for the arc-open variables. + +The right *engine* for netdes is **not** ``average_xhat_nonants``. +The averaged-data problem can leave some :math:`x_e` at 0 because the +average demand pattern does not need that arc; a real scenario with +peakier demand may need it. The averaged-solution path is +``lp_xbar_nonants``: any arc that any scenario's LP wanted positive +contributes positively to the average, and ``np.ceil`` then promotes +it to 1. + +``examples/netdes/netdes_auxiliary.py``: + +.. code-block:: python + + import numpy as np + from mpisppy.utils.xhat_helpers import lp_xbar_nonants + from netdes import scenario_creator, scenario_names_creator + + + def feasible_xhat_creator(*, solver_name, solver_options=None, + num_scens=None, **scenario_creator_kwargs): + if num_scens is None: + from parse import parse + num_scens = parse(scenario_creator_kwargs["path"], + scenario_ix=None)["K"] + snames = scenario_names_creator(num_scens) + arr = lp_xbar_nonants( + scenario_creator, snames, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": np.ceil(arr - 1e-9)} + +The :math:`-10^{-9}` margin keeps integer-valued LP solutions from +being inadvertently bumped up by floating-point dust. + +Worked example: sslp (binary, set-covering monotonicity) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Sslp ``model.FacilityOpen[j]`` is ``Binary``. Opening more facilities +never tightens ``DemandConstraint`` (more capacity available) or +``ClientConstraint`` (the LHS does not involve ``FacilityOpen``). The +shipped model also carries a high-``Penalty`` ``Dummy`` slack, so any +fixed candidate is technically feasible; the rounded LP-xbar is +still a meaningful low-slack candidate for the inner-bound spoke +that consumes it. + +Sslp does not currently ship an ``average_scenario_creator``, so the +auxiliary skips the ``average_xhat_nonants`` engine entirely and goes +straight to ``lp_xbar_nonants``. The feasibility-preserving rule +chosen here is ``np.round``. + +``examples/sslp/sslp_auxiliary.py``: + +.. code-block:: python + + import numpy as np + from mpisppy.utils.xhat_helpers import lp_xbar_nonants + from sslp import scenario_creator, scenario_names_creator + + + def feasible_xhat_creator(*, solver_name, solver_options=None, + num_scens, **scenario_creator_kwargs): + snames = scenario_names_creator(num_scens) + arr = lp_xbar_nonants( + scenario_creator, snames, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": np.round(arr)} + +See also +-------- + +* :ref:`jensens` -- Jensen's bound and the + ``--*-try-jensens-first`` flags. Shares the + ``average_scenario_creator`` convention but uses it for a different + contract (silently-skip-on-infeasibility candidate xhat, plus an + outer-bound path that ``feasible_xhat_creator`` does not address). +* :ref:`scenario_creator` -- the core scenario-module conventions + (``scenario_creator``, ``scenario_names_creator``, ...) that are + prerequisites for everything in this document. + +Heuristics fixing methods +------------------------- + +For some problems, you might have heuristic ways to fix many of the +nonanticipative variables. Once they are fixed, the resulting problem +might solve fairly quickly, which can be the basis for a +``feasible_xhat_creator`` function. diff --git a/doc/src/index.rst b/doc/src/index.rst index 33c47fea3..e4e616473 100644 --- a/doc/src/index.rst +++ b/doc/src/index.rst @@ -54,6 +54,7 @@ MPI is used. properbundles.rst pickling.rst jensens.rst + feasible_xhat.rst xhat_from_file.rst smps.rst agnostic.rst diff --git a/doc/src/jensens.rst b/doc/src/jensens.rst index 7cf43d148..db88865d1 100644 --- a/doc/src/jensens.rst +++ b/doc/src/jensens.rst @@ -27,6 +27,20 @@ Jensen's Bound as potential starting bound evaluated across the real scenarios, so Jensen's convexity assumption never enters the validity argument. +.. note:: + + The Jensen's xhat path **silently skips** a candidate that is + infeasible to pin in some real scenario; that is the right behavior + for an opportunistic inner-bound spoke. Downstream consumers that + pin a candidate first-stage and solve a per-scenario subproblem + with no fallback (e.g. pin-dual rho-setting algorithms) need a + strictly stronger contract -- a candidate guaranteed to be feasible + to pin in *every* real scenario. See :ref:`feasible_xhat` for the + ``feasible_xhat_creator`` convention and the + ``---try-feasible-xhat-first`` flags. Per xhat spoke, + ``--*-try-jensens-first`` and ``--*-try-feasible-xhat-first`` are + mutually exclusive. + What the options do ------------------- @@ -177,6 +191,17 @@ would produce a *different* model whose solution would be solving a different problem — at best an LP-relaxation heuristic that happens to share the same shape, not a Jensen's bound. +Even when ``average_scenario_creator`` is not directly usable, the +underscore helpers above (``_scenario_data``, ``_average_scenario_data``, +``_build_model``) may still be worth shipping. A +``feasible_xhat_creator`` (see :ref:`feasible_xhat`) is free to +build an averaged-data model internally — even one with relaxed +Var domains — to derive a starting first-stage point, because its +output is then projected and verified feasible in every real +scenario rather than handed off as a Jensen's bound. The +data-only-averaging principle does not bind ``feasible_xhat_creator``; +it only binds ``average_scenario_creator``. + When the model isn't amenable: sslp ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/examples/farmer/farmer_auxiliary.py b/examples/farmer/farmer_auxiliary.py new file mode 100644 index 000000000..b7dc91291 --- /dev/null +++ b/examples/farmer/farmer_auxiliary.py @@ -0,0 +1,60 @@ +############################################################################### +# 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. +############################################################################### +"""Auxiliary functions for the farmer example. + +Kept out of ``farmer.py`` so first-time readers of the introductory +example are not confronted with helpers they do not need. Functions +here are consumed by downstream tools (e.g. findW's pin-dual algorithm) +that need a first-stage point feasible for every real scenario's +per-scenario subproblem. +""" + +from mpisppy.utils.xhat_helpers import average_xhat_nonants +from farmer import average_scenario_creator + + +def feasible_xhat_creator( + *, + solver_name, + solver_options=None, + **scenario_creator_kwargs, +): + """Return a candidate first-stage ``DEVOTED_ACRES`` for farmer that + is feasible to fix in every real scenario's per-scenario + subproblem. + + Returns ``{"ROOT": np.ndarray}`` in + ``_mpisppy_node_list[0].nonant_vardata_list`` order. + + Strategy: solve the average scenario; return its first-stage values + unrounded. + + Why no rounding is needed: ``DEVOTED_ACRES`` is bounded + ``NonNegativeReals``, so any value the solver returns is a valid + integer-free pin. Farmer also has relatively complete recourse via + the buy/sell second-stage variables (``QuantityPurchased``, + ``QuantitySubQuotaSold``, ``QuantitySuperQuotaSold``), so any + feasible first-stage acreage allocation -- including the + average-scenario optimum -- is feasible to pin in every real + scenario. + + The continuous-first-stage case is the simplest one for the + ``feasible_xhat_creator`` convention. Callers always go through this + function rather than calling ``average_xhat_nonants`` directly so + that switching the underlying model to one with integer first-stage + (where rounding *is* needed) does not require changes at the call + site. + """ + arr = average_xhat_nonants( + average_scenario_creator, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": arr} diff --git a/examples/netdes/netdes_auxiliary.py b/examples/netdes/netdes_auxiliary.py new file mode 100644 index 000000000..1f0903099 --- /dev/null +++ b/examples/netdes/netdes_auxiliary.py @@ -0,0 +1,74 @@ +############################################################################### +# 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. +############################################################################### +"""Auxiliary functions for the netdes example. + +Kept out of ``netdes.py`` so first-time readers of the introductory +example are not confronted with helpers they do not need. Functions +here are consumed by downstream tools (e.g. findW's pin-dual algorithm) +that need a first-stage point feasible for every real scenario's +per-scenario subproblem. +""" + +import numpy as np + +from mpisppy.utils.xhat_helpers import lp_xbar_nonants +from netdes import scenario_creator, scenario_names_creator + + +def feasible_xhat_creator( + *, + solver_name, + solver_options=None, + num_scens=None, + **scenario_creator_kwargs, +): + """Return a candidate first-stage x for netdes that is feasible to + fix in every real scenario's per-scenario subproblem. + + Returns ``{"ROOT": np.ndarray}`` in + ``_mpisppy_node_list[0].nonant_vardata_list`` order, ready to pass + to ``Xhat_Eval._fix_nonants``. + + Strategy: solve the LP relaxation of every real scenario, take the + probability-weighted average of arc-open values (the LP-xbar), and + apply ``np.ceil``. The recourse constraint is + ``y[e] - u[e]*x[e] <= 0``; raising ``x[e]`` to 1 only loosens this + bound, so opening more arcs never tightens any per-scenario + subproblem. Ceiling the LP-xbar opens every arc that any scenario's + LP wanted positive, producing a strict cover that is feasible + when pinned in every real scenario. + + The small ``- 1e-9`` margin keeps integer-valued LP solutions from + being inadvertently bumped up by floating-point dust. + + Args: + solver_name: pyomo solver name. + solver_options: solver options dict. + num_scens: number of real scenarios to aggregate over. If None, + it is read from the instance file via ``parse``. + **scenario_creator_kwargs: forwarded to + ``netdes.scenario_creator`` (notably ``path``). + """ + if num_scens is None: + from parse import parse + path = scenario_creator_kwargs.get("path") + if path is None: + raise ValueError( + "feasible_xhat_creator needs num_scens or path kwarg" + ) + num_scens = parse(path, scenario_ix=None)["K"] + snames = scenario_names_creator(num_scens) + arr = lp_xbar_nonants( + scenario_creator, + snames, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": np.ceil(arr - 1e-9)} diff --git a/examples/run_all.py b/examples/run_all.py index c7d5835ef..555dac1d2 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -257,6 +257,22 @@ def do_one_mmw(dirname, modname, runefstring, npyfile, mmwargstring): "--max-solver-threads=2 " "--solver-name={}".format(solver_name)) + # Smoke-test netdes's feasible_xhat_creator end-to-end via the + # xhat-feasible-first flag. Discovery: cfg_vanilla auto-imports + # netdes_auxiliary because no feasible_xhat_creator is on netdes + # itself. Per spoke, --*-try-feasible-xhat-first is mutually + # exclusive with --*-try-jensens-first, so this run uses Jensen's + # only on the lagrangian outer-bound side. + do_one("netdes", "../../mpisppy/generic_cylinders.py", 3, + "--module-name netdes --max-iterations=3 " + "--instance-name=network-10-20-L-01 --netdes-data-path ./data " + "--rel-gap=0.0 --default-rho=10000 --presolve " + "--lagrangian --xhatshuffle " + "--lagrangian-try-jensens-first " + "--xhatshuffle-try-feasible-xhat-first " + "--max-solver-threads=2 " + "--solver-name={}".format(solver_name)) + # -------- Second part: sizes, sslp, hydro, aircond, MMW -------- if run_second_part: # sizes is slow for xpress so try linearizing the proximal term. diff --git a/examples/sslp/sslp_auxiliary.py b/examples/sslp/sslp_auxiliary.py new file mode 100644 index 000000000..b71a3b54c --- /dev/null +++ b/examples/sslp/sslp_auxiliary.py @@ -0,0 +1,63 @@ +############################################################################### +# 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. +############################################################################### +"""Auxiliary functions for the sslp example. + +Kept out of ``sslp.py`` so first-time readers of the introductory +example are not confronted with helpers they do not need. Functions +here are consumed by downstream tools (e.g. findW's pin-dual algorithm) +that need a first-stage point feasible for every real scenario's +per-scenario subproblem. +""" + +import numpy as np + +from mpisppy.utils.xhat_helpers import lp_xbar_nonants +from sslp import scenario_creator, scenario_names_creator + + +def feasible_xhat_creator( + *, + solver_name, + solver_options=None, + num_scens, + **scenario_creator_kwargs, +): + """Return a candidate first-stage ``FacilityOpen`` for sslp. + + Returns ``{"ROOT": np.ndarray}`` in + ``_mpisppy_node_list[0].nonant_vardata_list`` order. + + Strategy: solve the LP relaxation of every real scenario, take the + probability-weighted average of ``FacilityOpen`` (sslp ships + ``_mpisppy_probability = "uniform"``), then ``np.round`` to + integer. + + Why ``np.round`` is feasibility-preserving for sslp: more open + facilities never tightens ``DemandConstraint`` (more capacity) or + ``ClientConstraint`` (does not involve ``FacilityOpen``). The + shipped model also carries a high-``Penalty`` ``Dummy`` slack, so + any pin is technically feasible; the rounded LP-xbar is a + meaningful low-slack candidate for the pin-dual machinery. + + Args: + solver_name: pyomo solver name. + solver_options: solver options dict. + num_scens: number of scenarios to aggregate over. + **scenario_creator_kwargs: forwarded to + ``sslp.scenario_creator`` (e.g. ``data_dir``, ``surrogate``). + """ + snames = scenario_names_creator(num_scens) + arr = lp_xbar_nonants( + scenario_creator, + snames, + solver_name=solver_name, + scenario_creator_kwargs=scenario_creator_kwargs, + solver_options=solver_options, + ) + return {"ROOT": np.round(arr)} diff --git a/mpisppy/cylinders/_jensens_mixin.py b/mpisppy/cylinders/_jensens_mixin.py index dab924f92..ed8b1c202 100644 --- a/mpisppy/cylinders/_jensens_mixin.py +++ b/mpisppy/cylinders/_jensens_mixin.py @@ -214,3 +214,48 @@ def _try_average_scenario_xhat(self): Eobj = self._jensens_evaluate_xhat(cache) if Eobj is not None: self.update_if_improving(Eobj) + + def _feasible_xhat_enabled(self): + return "feasible_xhat" in self.opt.options + + # NOTE: this method is not Jensen's-specific despite living on + # _JensensMixin. The containing class will be renamed in a + # subsequent PR to reflect that it covers the generic pre-loop + # xhat candidate machinery, not just Jensen's xhat. + def _try_feasible_xhat(self): + """One-shot helper for xhat spokes that opt in via + ``--*-try-feasible-xhat-first``. No-op when the flag is off. + + Mutually exclusive with ``--*-try-jensens-first`` (cfg_vanilla + raises if both are set on the same spoke). + + Unlike the Jensen's path, no average-scenario solve is needed: + the user's ``feasible_xhat_creator`` returns the candidate + directly, in the same dict-by-nodename cache form the rest of + the xhat path consumes. We pin it, evaluate it across all real + scenarios, and -- the contract of ``feasible_xhat_creator`` + guarantees feasibility -- expect ``no_incumbent_prob() == 0``. + We still go through ``_jensens_evaluate_xhat`` so an + unexpected per-scenario infeasibility surfaces as a silent skip + rather than a crash. + """ + if not self._feasible_xhat_enabled(): + return + f = self.opt.options["feasible_xhat"] + creator = f["feasible_xhat_creator"] + kwargs = f.get("scenario_creator_kwargs") or {} + cache = creator( + solver_name=self.opt.options["solver_name"], + solver_options=self.opt.options.get("iterk_solver_options") or None, + **kwargs, + ) + if not isinstance(cache, dict) or "ROOT" not in cache: + raise RuntimeError( + "feasible_xhat_creator must return a dict with at least " + "a 'ROOT' key; got " + f"{type(cache).__name__} with keys " + f"{list(cache.keys()) if isinstance(cache, dict) else 'n/a'}." + ) + Eobj = self._jensens_evaluate_xhat(cache) + if Eobj is not None: + self.update_if_improving(Eobj) diff --git a/mpisppy/cylinders/xhatlooper_bounder.py b/mpisppy/cylinders/xhatlooper_bounder.py index 6ada3ebe8..7d4ab4e44 100644 --- a/mpisppy/cylinders/xhatlooper_bounder.py +++ b/mpisppy/cylinders/xhatlooper_bounder.py @@ -32,8 +32,10 @@ def main(self): xhatter = self.xhat_prep() - # No-op unless --xhatlooper-try-jensens-first is set. + # No-ops unless --xhatlooper-try-jensens-first / + # --xhatlooper-try-feasible-xhat-first are set (mutually exclusive). self._try_average_scenario_xhat() + self._try_feasible_xhat() scen_limit = self.opt.options['xhat_looper_options']['scen_limit'] diff --git a/mpisppy/cylinders/xhatshufflelooper_bounder.py b/mpisppy/cylinders/xhatshufflelooper_bounder.py index 0059e49e8..f2e9c32b5 100644 --- a/mpisppy/cylinders/xhatshufflelooper_bounder.py +++ b/mpisppy/cylinders/xhatshufflelooper_bounder.py @@ -68,9 +68,11 @@ def main(self): self.xhat_prep() - # No-op unless --xhatshuffle-try-jensens-first is set. Tolerates - # integer recourse and per-scenario infeasibility (silent skip). + # No-ops unless --xhatshuffle-try-jensens-first / + # --xhatshuffle-try-feasible-xhat-first are set (mutually exclusive). + # Both tolerate per-scenario infeasibility via silent skip. self._try_average_scenario_xhat() + self._try_feasible_xhat() if "reverse" in self.opt.options["xhat_looper_options"]: self.reverse = self.opt.options["xhat_looper_options"]["reverse"] diff --git a/mpisppy/cylinders/xhatspecific_bounder.py b/mpisppy/cylinders/xhatspecific_bounder.py index 0d7518a48..a45f4b82a 100644 --- a/mpisppy/cylinders/xhatspecific_bounder.py +++ b/mpisppy/cylinders/xhatspecific_bounder.py @@ -43,8 +43,10 @@ def main(self): xhatter = self.xhat_prep() - # No-op unless --xhatspecific-try-jensens-first is set. + # No-ops unless --xhatspecific-try-jensens-first / + # --xhatspecific-try-feasible-xhat-first are set (mutually exclusive). self._try_average_scenario_xhat() + self._try_feasible_xhat() ib_iter = 1 # ib is for inner bound while (not self.got_kill_signal()): diff --git a/mpisppy/cylinders/xhatxbar_bounder.py b/mpisppy/cylinders/xhatxbar_bounder.py index c5f4ae202..6d0db7a89 100644 --- a/mpisppy/cylinders/xhatxbar_bounder.py +++ b/mpisppy/cylinders/xhatxbar_bounder.py @@ -62,8 +62,10 @@ def main(self): xhatter = self.xhat_prep() - # No-op unless --xhatxbar-try-jensens-first is set. + # No-ops unless --xhatxbar-try-jensens-first / + # --xhatxbar-try-feasible-xhat-first are set (mutually exclusive). self._try_average_scenario_xhat() + self._try_feasible_xhat() ib_iter = 1 # ib is for inner bound while (not self.got_kill_signal()): diff --git a/mpisppy/generic/decomp.py b/mpisppy/generic/decomp.py index d0822440e..418d950a5 100644 --- a/mpisppy/generic/decomp.py +++ b/mpisppy/generic/decomp.py @@ -63,10 +63,14 @@ def do_decomp(module, cfg, scenario_creator, scenario_creator_kwargs, average_scenario_creator = getattr(module, "average_scenario_creator", None) + from mpisppy.utils import cfg_vanilla as vanilla + feasible_xhat_creator = vanilla._find_feasible_xhat_creator(module, cfg) + list_of_spoke_dict = build_spoke_list(cfg, beans, scenario_creator_kwargs, rho_setter, all_nodenames, variable_probability=variable_probability, - average_scenario_creator=average_scenario_creator) + average_scenario_creator=average_scenario_creator, + feasible_xhat_creator=feasible_xhat_creator) # if the user dares, let them mess with the hubdict prior to solve if hasattr(module, 'hub_and_spoke_dict_callback'): diff --git a/mpisppy/generic/spokes.py b/mpisppy/generic/spokes.py index 72e4a5e54..123eee05b 100644 --- a/mpisppy/generic/spokes.py +++ b/mpisppy/generic/spokes.py @@ -16,7 +16,8 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, rho_setter, all_nodenames, variable_probability=None, - average_scenario_creator=None): + average_scenario_creator=None, + feasible_xhat_creator=None): """Build and return the list of spoke dicts for WheelSpinner. Args: @@ -28,6 +29,9 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, variable_probability: variable probability list or None (used by ADMM) average_scenario_creator: module's average_scenario_creator (or None), needed by any spoke whose --*-try-jensens-first flag is set. + feasible_xhat_creator: module's feasible_xhat_creator (or None), + needed by any xhat spoke whose --*-try-feasible-xhat-first + flag is set. Returns: list: list of spoke dicts @@ -108,7 +112,8 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, xhatshuffle_spoke = vanilla.xhatshuffle_spoke(*beans, scenario_creator_kwargs=scenario_creator_kwargs, all_nodenames=all_nodenames, - average_scenario_creator=average_scenario_creator) + average_scenario_creator=average_scenario_creator, + feasible_xhat_creator=feasible_xhat_creator) # special code for multi-stage (e.g., hydro) if cfg.get("stage2_ef_solver_name") is not None: xhatshuffle_spoke["opt_kwargs"]["options"]["stage2_ef_solver_name"] = cfg["stage2_ef_solver_name"] @@ -119,7 +124,8 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, scenario_creator_kwargs=scenario_creator_kwargs, variable_probability=variable_probability, all_nodenames=all_nodenames, - average_scenario_creator=average_scenario_creator) + average_scenario_creator=average_scenario_creator, + feasible_xhat_creator=feasible_xhat_creator) # reduced cost fixer options setup if cfg.reduced_costs: diff --git a/mpisppy/tests/test_feasible_xhat.py b/mpisppy/tests/test_feasible_xhat.py new file mode 100644 index 000000000..9fac1789d --- /dev/null +++ b/mpisppy/tests/test_feasible_xhat.py @@ -0,0 +1,436 @@ +############################################################################### +# 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. +############################################################################### +"""Tests for the ``feasible_xhat_creator`` convention and its +``average_xhat_nonants`` helper. + +Two prototypes are exercised: netdes (LP-relax + ceil, via the helper) +and sslp (LP-relax + average + round, rolls own). The convention's +contract is that the returned candidate must be feasible to fix in +every real scenario's per-scenario subproblem; the netdes test +verifies this directly. +""" + +import os +import sys +import unittest + +import numpy as np +import pyomo.environ as pyo +from pyomo.opt import TerminationCondition + +import mpisppy.utils.sputils as sputils +from mpisppy.utils.xhat_helpers import average_xhat_nonants +from mpisppy.tests.utils import get_solver + +_EXAMPLES_DIR = os.path.join( + os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))), + "examples", +) +sys.path.insert(0, os.path.join(_EXAMPLES_DIR, "farmer")) +sys.path.insert(0, os.path.join(_EXAMPLES_DIR, "netdes")) +sys.path.insert(0, os.path.join(_EXAMPLES_DIR, "sslp")) + +import farmer # noqa: E402 +import farmer_auxiliary # noqa: E402 +import netdes # noqa: E402 +import netdes_auxiliary # noqa: E402 +import sslp_auxiliary # noqa: E402 + +solver_available, solver_name, _, _ = get_solver() + +_NETDES_DATA = os.path.join( + _EXAMPLES_DIR, "netdes", "data", "network-10-20-L-01.dat" +) +_SSLP_DATA_DIR = os.path.join( + _EXAMPLES_DIR, "sslp", "data", "sslp_15_45_5", "scenariodata" +) + + +class TestAverageXhatNonantsContract(unittest.TestCase): + """Solver-free contract checks for the helper.""" + + def test_rejects_no_node_list(self): + def bad_creator(name): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.obj = pyo.Objective(expr=m.x) + return m + + with self.assertRaises(RuntimeError) as ctx: + average_xhat_nonants(bad_creator, solver_name="cplex") + self.assertIn("_mpisppy_node_list", str(ctx.exception)) + + def test_rejects_multi_stage(self): + def two_node_creator(name): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.fsc = pyo.Expression(expr=m.x) + m.obj = pyo.Objective(expr=m.x + m.y) + sputils.attach_root_node(m, m.fsc, [m.x]) + # tack on a fake second node + from mpisppy.scenario_tree import ScenarioNode + m._mpisppy_node_list.append( + ScenarioNode("STAGE2", 1.0, 2, m.fsc, [m.y], m, + parent_name="ROOT") + ) + return m + + with self.assertRaises(RuntimeError) as ctx: + average_xhat_nonants(two_node_creator, solver_name="cplex") + self.assertIn("two-stage only", str(ctx.exception)) + + +@unittest.skipIf(not solver_available, "no solver available") +class TestAverageXhatNonantsOnFarmer(unittest.TestCase): + """End-to-end on farmer (whose first stage is continuous).""" + + def test_returns_array_of_correct_length(self): + arr = average_xhat_nonants( + farmer.average_scenario_creator, + solver_name=solver_name, + scenario_creator_kwargs={"num_scens": 6}, + ) + self.assertIsInstance(arr, np.ndarray) + self.assertEqual(arr.shape, (3,)) # DEVOTED_ACRES over 3 crops + + def test_relax_integrality_is_no_op_on_continuous_first_stage(self): + a = average_xhat_nonants( + farmer.average_scenario_creator, + solver_name=solver_name, + scenario_creator_kwargs={"num_scens": 6}, + ) + b = average_xhat_nonants( + farmer.average_scenario_creator, + solver_name=solver_name, + scenario_creator_kwargs={"num_scens": 6}, + relax_integrality=True, + ) + np.testing.assert_allclose(a, b, atol=1e-6) + + +@unittest.skipIf(not solver_available, "no solver available") +class TestFarmerFeasibleXhatCreator(unittest.TestCase): + """Continuous first-stage: convention is satisfied by the + average-scenario optimum unchanged.""" + + def test_returns_root_array(self): + cache = farmer_auxiliary.feasible_xhat_creator( + solver_name=solver_name, num_scens=6, + ) + self.assertIn("ROOT", cache) + arr = cache["ROOT"] + self.assertEqual(arr.shape, (3,)) # DEVOTED_ACRES over 3 crops + self.assertTrue(np.all(np.isfinite(arr))) + self.assertTrue(np.all(arr >= -1e-9)) + + def test_kwargs_thread_through(self): + a = farmer_auxiliary.feasible_xhat_creator( + solver_name=solver_name, num_scens=6, seedoffset=0, + ) + b = farmer_auxiliary.feasible_xhat_creator( + solver_name=solver_name, num_scens=6, seedoffset=7, + ) + self.assertFalse(np.allclose(a["ROOT"], b["ROOT"]), + "seedoffset did not change the candidate") + + +@unittest.skipIf(not solver_available, "no solver available") +class TestNetdesFeasibleXhatCreator(unittest.TestCase): + """End-to-end check: netdes_auxiliary.feasible_xhat_creator returns + a candidate that is integer-valued and feasible to fix in every + real scenario.""" + + def setUp(self): + from parse import parse + full = parse(_NETDES_DATA, scenario_ix=None) + self.K = full["K"] + self.cache = netdes_auxiliary.feasible_xhat_creator( + solver_name=solver_name, path=_NETDES_DATA, + ) + + def test_returns_root_dict(self): + self.assertIn("ROOT", self.cache) + self.assertEqual(set(self.cache.keys()), {"ROOT"}) + + def test_values_are_integer(self): + arr = self.cache["ROOT"] + self.assertTrue(np.allclose(arr, np.round(arr), atol=1e-9), + f"ceil output is not integer-valued: {arr}") + + def test_fixes_are_feasible_in_every_real_scenario(self): + arr = self.cache["ROOT"] + for k in range(self.K): + sname = f"Scenario{k}" + scen = netdes.scenario_creator(sname, path=_NETDES_DATA) + root = scen._mpisppy_node_list[0] + for v, val in zip(root.nonant_vardata_list, arr): + v.fix(val) + solver = pyo.SolverFactory(solver_name) + if sputils.is_persistent(solver): + solver.set_instance(scen) + results = solver.solve(tee=False) + else: + results = solver.solve(scen, tee=False) + tc = results.solver.termination_condition + self.assertIn( + tc, (TerminationCondition.optimal, TerminationCondition.feasible), + f"netdes fix infeasible on {sname}: tc={tc}, xhat={arr}", + ) + + +@unittest.skipIf(not solver_available, "no solver available") +@unittest.skipIf(not os.path.isdir(_SSLP_DATA_DIR), "sslp data not found") +class TestSslpFeasibleXhatCreator(unittest.TestCase): + + def test_returns_integer_root_dict(self): + cache = sslp_auxiliary.feasible_xhat_creator( + solver_name=solver_name, + num_scens=5, + data_dir=_SSLP_DATA_DIR, + ) + self.assertIn("ROOT", cache) + arr = cache["ROOT"] + self.assertEqual(arr.shape, (15,)) # 15 servers in sslp_15_45_5 + self.assertTrue(np.allclose(arr, np.round(arr), atol=1e-9)) + self.assertTrue(np.all((arr >= 0) & (arr <= 1)), + f"FacilityOpen out of [0,1]: {arr}") + + +class TestMaybeAttachFeasibleXhat(unittest.TestCase): + """cfg_vanilla._maybe_attach_feasible_xhat behavior.""" + + def _cfg_with_flag(self, prefix, flag_on=False, jensens_on=False): + import mpisppy.utils.config as cfg_mod + cfg = cfg_mod.Config() + if prefix == "xhatshuffle": + cfg.xhatshuffle_args() + elif prefix == "xhatxbar": + cfg.xhatxbar_args() + elif prefix == "xhatlooper": + cfg.xhatlooper_args() + elif prefix == "xhatspecific": + cfg.xhatspecific_args() + if flag_on: + cfg[f"{prefix}_try_feasible_xhat_first"] = True + if jensens_on: + cfg[f"{prefix}_try_jensens_first"] = True + return cfg + + def test_noop_when_flag_off(self): + from mpisppy.utils.cfg_vanilla import _maybe_attach_feasible_xhat + cfg = self._cfg_with_flag("xhatshuffle", flag_on=False) + spoke_dict = {"opt_kwargs": {"options": {}}} + _maybe_attach_feasible_xhat(spoke_dict, cfg, "xhatshuffle", + feasible_xhat_creator=None, + scenario_creator_kwargs={}) + self.assertNotIn("feasible_xhat", spoke_dict["opt_kwargs"]["options"]) + + def test_raises_when_creator_is_none(self): + from mpisppy.utils.cfg_vanilla import _maybe_attach_feasible_xhat + cfg = self._cfg_with_flag("xhatshuffle", flag_on=True) + spoke_dict = {"opt_kwargs": {"options": {}}} + with self.assertRaises(RuntimeError) as ctx: + _maybe_attach_feasible_xhat(spoke_dict, cfg, "xhatshuffle", + feasible_xhat_creator=None, + scenario_creator_kwargs={}) + self.assertIn("feasible_xhat_creator", str(ctx.exception)) + + def test_raises_when_both_first_flags_set(self): + from mpisppy.utils.cfg_vanilla import _maybe_attach_feasible_xhat + cfg = self._cfg_with_flag("xhatshuffle", flag_on=True, jensens_on=True) + spoke_dict = {"opt_kwargs": {"options": {}}} + with self.assertRaises(RuntimeError) as ctx: + _maybe_attach_feasible_xhat(spoke_dict, cfg, "xhatshuffle", + feasible_xhat_creator=lambda **k: None, + scenario_creator_kwargs={}) + msg = str(ctx.exception) + self.assertIn("mutually exclusive", msg) + self.assertIn("try-jensens-first", msg) + self.assertIn("try-feasible-xhat-first", msg) + + def test_installs_dict_when_flag_on(self): + from mpisppy.utils.cfg_vanilla import _maybe_attach_feasible_xhat + cfg = self._cfg_with_flag("xhatshuffle", flag_on=True) + spoke_dict = {"opt_kwargs": {"options": {}}} + kwargs = {"path": "/x"} + _maybe_attach_feasible_xhat( + spoke_dict, cfg, "xhatshuffle", + feasible_xhat_creator=netdes_auxiliary.feasible_xhat_creator, + scenario_creator_kwargs=kwargs, + ) + f = spoke_dict["opt_kwargs"]["options"]["feasible_xhat"] + self.assertIs(f["feasible_xhat_creator"], + netdes_auxiliary.feasible_xhat_creator) + self.assertIs(f["scenario_creator_kwargs"], kwargs) + + +class TestFindFeasibleXhatCreator(unittest.TestCase): + """Discovery helper: tries main module, falls back to _auxiliary.""" + + def _cfg(self, **flags): + import mpisppy.utils.config as cfg_mod + cfg = cfg_mod.Config() + cfg.xhatshuffle_args() + cfg.xhatxbar_args() + cfg.xhatlooper_args() + cfg.xhatspecific_args() + for name, val in flags.items(): + cfg[name] = val + return cfg + + def test_returns_none_when_no_flag_set(self): + from mpisppy.utils.cfg_vanilla import _find_feasible_xhat_creator + cfg = self._cfg() + # netdes is the test module; even though it has an _auxiliary, + # no flag is set so we never reach for it + result = _find_feasible_xhat_creator(netdes, cfg) + self.assertIsNone(result) + + def test_falls_back_to_auxiliary(self): + from mpisppy.utils.cfg_vanilla import _find_feasible_xhat_creator + cfg = self._cfg(xhatshuffle_try_feasible_xhat_first=True) + # netdes does NOT export feasible_xhat_creator on the main + # module; only netdes_auxiliary does. + self.assertIsNone(getattr(netdes, "feasible_xhat_creator", None)) + result = _find_feasible_xhat_creator(netdes, cfg) + self.assertIs(result, netdes_auxiliary.feasible_xhat_creator) + + def test_finds_on_main_module_when_present(self): + # Synthesize a fake module that exposes feasible_xhat_creator + # directly, to confirm the main-module lookup wins. + import types + from mpisppy.utils.cfg_vanilla import _find_feasible_xhat_creator + marker = lambda **k: {"ROOT": np.array([0.0])} # noqa: E731 + fake = types.ModuleType("fake_module_with_feas") + fake.feasible_xhat_creator = marker + cfg = self._cfg(xhatxbar_try_feasible_xhat_first=True) + result = _find_feasible_xhat_creator(fake, cfg) + self.assertIs(result, marker) + + def test_raises_when_neither_has_it(self): + # farmer's main module has no feasible_xhat_creator either, + # but farmer_auxiliary does -- so swap to a module without an + # auxiliary file at all. + import types + from mpisppy.utils.cfg_vanilla import _find_feasible_xhat_creator + bare = types.ModuleType("bare_module_for_feas_test") + cfg = self._cfg(xhatlooper_try_feasible_xhat_first=True) + with self.assertRaises(RuntimeError) as ctx: + _find_feasible_xhat_creator(bare, cfg) + self.assertIn("feasible_xhat_creator", str(ctx.exception)) + + +class _StubFeasOpt: + """Stand-in for self.opt that captures _fix_nonants / solve_loop / + Eobjective calls so we can exercise _try_feasible_xhat without MPI.""" + + def __init__(self, *, infeas_prob, eobj, options=None): + self._infeas_prob = infeas_prob + self._eobj = eobj + self.options = options if options is not None else {} + self.fix_calls = [] + self.solve_loop_calls = [] + + def _fix_nonants(self, cache): + self.fix_calls.append(cache) + + def solve_loop(self, **kwargs): + self.solve_loop_calls.append(kwargs) + + def no_incumbent_prob(self): + return self._infeas_prob + + def Eobjective(self, verbose=False): + return self._eobj + + +class _StubFeasSpoke: + """Mixin host wired to capture update_if_improving calls.""" + + def __init__(self, opt): + from mpisppy.cylinders._jensens_mixin import _JensensMixin + self.opt = opt + self.bound_updates = [] + # bind the mixin methods + self._feasible_xhat_enabled = lambda: _JensensMixin._feasible_xhat_enabled(self) + self._jensens_evaluate_xhat = lambda c: _JensensMixin._jensens_evaluate_xhat(self, c) + self._try_feasible_xhat = lambda: _JensensMixin._try_feasible_xhat(self) + + def update_if_improving(self, Eobj): + self.bound_updates.append(Eobj) + + +class TestTryFeasibleXhat(unittest.TestCase): + """Unit tests for _JensensMixin._try_feasible_xhat.""" + + def test_noop_when_flag_off(self): + opt = _StubFeasOpt(infeas_prob=0.0, eobj=1.0) # no "feasible_xhat" key + sp = _StubFeasSpoke(opt) + sp._try_feasible_xhat() + self.assertEqual(sp.bound_updates, []) + self.assertEqual(opt.fix_calls, []) + + def test_updates_bound_when_feasible(self): + creator = lambda **kw: {"ROOT": np.array([1.0, 2.0])} # noqa: E731 + opt = _StubFeasOpt( + infeas_prob=0.0, eobj=42.0, + options={ + "feasible_xhat": { + "feasible_xhat_creator": creator, + "scenario_creator_kwargs": {}, + }, + "solver_name": "fake", + "iterk_solver_options": None, + }, + ) + sp = _StubFeasSpoke(opt) + sp._try_feasible_xhat() + self.assertEqual(sp.bound_updates, [42.0]) + self.assertEqual(len(opt.fix_calls), 1) + np.testing.assert_array_equal(opt.fix_calls[0]["ROOT"], + np.array([1.0, 2.0])) + + def test_skips_silently_when_infeasible(self): + creator = lambda **kw: {"ROOT": np.array([0.0])} # noqa: E731 + opt = _StubFeasOpt( + infeas_prob=0.5, eobj=999.0, # ignored on infeasible path + options={ + "feasible_xhat": { + "feasible_xhat_creator": creator, + "scenario_creator_kwargs": {}, + }, + "solver_name": "fake", + }, + ) + sp = _StubFeasSpoke(opt) + sp._try_feasible_xhat() + self.assertEqual(sp.bound_updates, []) + + def test_raises_on_bad_return_shape(self): + # creator returns a bare array, not a dict + creator = lambda **kw: np.array([1.0]) # noqa: E731 + opt = _StubFeasOpt( + infeas_prob=0.0, eobj=1.0, + options={ + "feasible_xhat": { + "feasible_xhat_creator": creator, + "scenario_creator_kwargs": {}, + }, + "solver_name": "fake", + }, + ) + sp = _StubFeasSpoke(opt) + with self.assertRaises(RuntimeError) as ctx: + sp._try_feasible_xhat() + self.assertIn("ROOT", str(ctx.exception)) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 4b53e056d..a83ed7b93 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -52,6 +52,82 @@ def _maybe_attach_jensens(spoke_dict, cfg, spoke_prefix, "scenario_creator_kwargs": scenario_creator_kwargs, } +def _find_feasible_xhat_creator(module, cfg): + """Look up ``feasible_xhat_creator`` only when at least one + ``---try-feasible-xhat-first`` flag is set. Tries the main + scenario module first; falls back to importing + ``_auxiliary`` and looking there. + + Returns ``None`` when no flag is set; raises if any flag is set + but the function cannot be located on the module or in + ``_auxiliary``. + """ + xhat_prefixes = ("xhatshuffle", "xhatxbar", "xhatlooper", "xhatspecific") + flags_set = [ + p for p in xhat_prefixes + if cfg.get(f"{p}_try_feasible_xhat_first", False) + ] + if not flags_set: + return None + fn = getattr(module, "feasible_xhat_creator", None) + if fn is not None: + return fn + import importlib + aux_name = f"{module.__name__}_auxiliary" + flag_cli = f"--{flags_set[0].replace('_', '-')}-try-feasible-xhat-first" + try: + aux = importlib.import_module(aux_name) + except ImportError: + raise RuntimeError( + f"{flag_cli} was set, but feasible_xhat_creator was not " + f"found on {module.__name__} and {aux_name} could not be " + f"imported. Define feasible_xhat_creator on the module or " + f"create {aux_name} with one, or turn the flag off." + ) + fn = getattr(aux, "feasible_xhat_creator", None) + if fn is None: + raise RuntimeError( + f"{flag_cli} was set, but neither {module.__name__} nor " + f"{aux_name} defines feasible_xhat_creator." + ) + return fn + + +def _maybe_attach_feasible_xhat(spoke_dict, cfg, spoke_prefix, + feasible_xhat_creator, scenario_creator_kwargs): + """Attach feasible_xhat options to spoke_dict when the corresponding + ---try-feasible-xhat-first flag is set. + + The spoke reads these at runtime via self.opt.options["feasible_xhat"] + and drives _JensensMixin._try_feasible_xhat. See doc/src/feasible_xhat.rst. + + Raises if the user enabled both ---try-jensens-first and + ---try-feasible-xhat-first on the same xhat spoke; the + two pre-loop xhat candidates are mutually exclusive per spoke. + """ + flag = f"{spoke_prefix}_try_feasible_xhat_first" + if not cfg.get(flag, False): + return + jensens_flag = f"{spoke_prefix}_try_jensens_first" + if cfg.get(jensens_flag, False): + raise RuntimeError( + f"--{spoke_prefix.replace('_', '-')}-try-jensens-first and " + f"--{spoke_prefix.replace('_', '-')}-try-feasible-xhat-first " + "were both set; they are mutually exclusive per xhat spoke." + ) + if feasible_xhat_creator is None: + raise RuntimeError( + f"--{spoke_prefix.replace('_', '-')}-try-feasible-xhat-first " + "was set, but the scenario module does not define " + "feasible_xhat_creator (looked up on the module and in " + "_auxiliary). Either implement feasible_xhat_creator, " + "or turn the flag off." + ) + spoke_dict["opt_kwargs"]["options"]["feasible_xhat"] = { + "feasible_xhat_creator": feasible_xhat_creator, + "scenario_creator_kwargs": scenario_creator_kwargs, + } + def shared_options(cfg, is_hub=False): shoptions = { "solver_name": cfg.solver_name, @@ -1135,6 +1211,7 @@ def xhatlooper_spoke( ph_extensions=None, extension_kwargs=None, average_scenario_creator=None, + feasible_xhat_creator=None, ): from mpisppy.cylinders.xhatlooper_bounder import XhatLooperInnerBound @@ -1157,6 +1234,8 @@ def xhatlooper_spoke( } _maybe_attach_jensens(xhatlooper_dict, cfg, "xhatlooper", average_scenario_creator, scenario_creator_kwargs) + _maybe_attach_feasible_xhat(xhatlooper_dict, cfg, "xhatlooper", + feasible_xhat_creator, scenario_creator_kwargs) return xhatlooper_dict @@ -1172,6 +1251,7 @@ def xhatxbar_spoke( extension_kwargs=None, all_nodenames=None, average_scenario_creator=None, + feasible_xhat_creator=None, ): from mpisppy.cylinders.xhatxbar_bounder import XhatXbarInnerBound xhatxbar_dict = _Xhat_Eval_spoke_foundation( @@ -1191,10 +1271,12 @@ def xhatxbar_spoke( "dump_prefix": "delme", "csvname": "looper.csv", } - + xhatxbar_dict["opt_kwargs"]["variable_probability"] = variable_probability _maybe_attach_jensens(xhatxbar_dict, cfg, "xhatxbar", average_scenario_creator, scenario_creator_kwargs) + _maybe_attach_feasible_xhat(xhatxbar_dict, cfg, "xhatxbar", + feasible_xhat_creator, scenario_creator_kwargs) return xhatxbar_dict @@ -1209,11 +1291,12 @@ def xhatshuffle_spoke( ph_extensions=None, extension_kwargs=None, average_scenario_creator=None, + feasible_xhat_creator=None, ): from mpisppy.cylinders.xhatshufflelooper_bounder import XhatShuffleInnerBound xhatshuffle_dict = _Xhat_Eval_spoke_foundation( - XhatShuffleInnerBound, + XhatShuffleInnerBound, cfg, scenario_creator, scenario_denouement, @@ -1234,6 +1317,8 @@ def xhatshuffle_spoke( xhatshuffle_dict["opt_kwargs"]["options"]["xhatshuffle_iter_step"] = cfg.xhatshuffle_iter_step _maybe_attach_jensens(xhatshuffle_dict, cfg, "xhatshuffle", average_scenario_creator, scenario_creator_kwargs) + _maybe_attach_feasible_xhat(xhatshuffle_dict, cfg, "xhatshuffle", + feasible_xhat_creator, scenario_creator_kwargs) return xhatshuffle_dict @@ -1249,6 +1334,7 @@ def xhatspecific_spoke( ph_extensions=None, extension_kwargs=None, average_scenario_creator=None, + feasible_xhat_creator=None, ): from mpisppy.cylinders.xhatspecific_bounder import XhatSpecificInnerBound @@ -1270,6 +1356,8 @@ def xhatspecific_spoke( } _maybe_attach_jensens(xhatspecific_dict, cfg, "xhatspecific", average_scenario_creator, scenario_creator_kwargs) + _maybe_attach_feasible_xhat(xhatspecific_dict, cfg, "xhatspecific", + feasible_xhat_creator, scenario_creator_kwargs) return xhatspecific_dict def xhatlshaped_spoke( diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 0c040fb68..3bf4b4b52 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -901,6 +901,16 @@ def xhatlooper_args(self): domain=bool, default=False) + self.add_to_config('xhatlooper_try_feasible_xhat_first', + description="before entering the xhatlooper main loop, take " + "the candidate first-stage from the scenario " + "module's feasible_xhat_creator (looked up on the " + "module or in _auxiliary) and try it as " + "an xhat. Two-stage only. Mutually exclusive " + "with --xhatlooper-try-jensens-first.", + domain=bool, + default=False) + def xhatshuffle_args(self): self.add_to_config('xhatshuffle', @@ -927,6 +937,16 @@ def xhatshuffle_args(self): domain=bool, default=False) + self.add_to_config('xhatshuffle_try_feasible_xhat_first', + description="before entering the xhatshuffle main loop, take " + "the candidate first-stage from the scenario " + "module's feasible_xhat_creator (looked up on the " + "module or in _auxiliary) and try it as " + "an xhat. Two-stage only. Mutually exclusive " + "with --xhatshuffle-try-jensens-first.", + domain=bool, + default=False) + self.add_stage2_ef_solver_name_arg() @@ -979,6 +999,16 @@ def xhatspecific_args(self): domain=bool, default=False) + self.add_to_config('xhatspecific_try_feasible_xhat_first', + description="before entering the xhatspecific main loop, take " + "the candidate first-stage from the scenario " + "module's feasible_xhat_creator (looked up on the " + "module or in _auxiliary) and try it as " + "an xhat. Two-stage only. Mutually exclusive " + "with --xhatspecific-try-jensens-first.", + domain=bool, + default=False) + def xhatxbar_args(self): @@ -996,6 +1026,16 @@ def xhatxbar_args(self): domain=bool, default=False) + self.add_to_config('xhatxbar_try_feasible_xhat_first', + description="before entering the xhatxbar main loop, take " + "the candidate first-stage from the scenario " + "module's feasible_xhat_creator (looked up on the " + "module or in _auxiliary) and try it as " + "an xhat. Two-stage only. Mutually exclusive " + "with --xhatxbar-try-jensens-first.", + domain=bool, + default=False) + def xhatlshaped_args(self): # we will not try to get the specification from the command line diff --git a/mpisppy/utils/xhat_helpers.py b/mpisppy/utils/xhat_helpers.py new file mode 100644 index 000000000..949920d8b --- /dev/null +++ b/mpisppy/utils/xhat_helpers.py @@ -0,0 +1,140 @@ +############################################################################### +# 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. +############################################################################### +"""Helpers for producing candidate xhat values from a deterministic proxy +of a stochastic program. + +Used by ``feasible_xhat_creator`` implementations in +``examples//_auxiliary.py`` and by downstream consumers +(e.g. findW's pin-dual algorithm) that need a first-stage point feasible +to fix in every real scenario's per-scenario subproblem. + +The Jensen's xhat path inside mpi-sppy tolerates a per-scenario-infeasible +candidate by silently skipping it; downstream consumers that pin a +candidate and solve a per-scenario MIP cannot, so they need a stronger +contract. ``feasible_xhat_creator`` is the entry point for that +contract; ``average_xhat_nonants`` is the common-case engine that the +implementations call when they have an ``average_scenario_creator``. + +The rounding/repair rule that turns the average-scenario solution into +a feasible candidate is model-specific (depends on monotonicity of +recourse feasibility in each first-stage variable) and lives in each +module's ``feasible_xhat_creator``, not here. +""" + +import numpy as np +import pyomo.environ as pyo +from pyomo.opt import TerminationCondition + +import mpisppy.utils.sputils as sputils + + +def _solve_and_extract_root(model, solver_name, solver_options): + solver = pyo.SolverFactory(solver_name) + for k, v in (solver_options or {}).items(): + solver.options[k] = v + if sputils.is_persistent(solver): + solver.set_instance(model) + results = solver.solve(tee=False) + else: + results = solver.solve(model, tee=False) + tc = results.solver.termination_condition + if tc not in (TerminationCondition.optimal, TerminationCondition.feasible): + raise RuntimeError( + f"solve in xhat_helpers failed with termination_condition={tc}" + ) + root = model._mpisppy_node_list[0] + return np.array( + [pyo.value(v) for v in root.nonant_vardata_list], dtype="d" + ) + + +def _check_two_stage(model, who): + if not hasattr(model, "_mpisppy_node_list"): + raise RuntimeError( + f"{who}: scenario must have _mpisppy_node_list attached " + "(via sputils.attach_root_node)." + ) + if len(model._mpisppy_node_list) != 1: + raise RuntimeError( + f"{who} is two-stage only; got " + f"{len(model._mpisppy_node_list)} tree nodes." + ) + + +def average_xhat_nonants( + average_scenario_creator, + *, + solver_name, + scenario_creator_kwargs=None, + solver_options=None, + relax_integrality=False, + scenario_name="AverageScenario", +): + """Build, optionally LP-relax, and solve the average scenario; return + its ROOT first-stage values as a 1-D ``np.ndarray`` in + ``_mpisppy_node_list[0].nonant_vardata_list`` order. + + Two-stage only. Returns the bare array; callers that need the + ``Xhat_Eval._fix_nonants`` cache form wrap as ``{"ROOT": arr}``. + + Note: solving the (LP-relaxed) average scenario is **not** the same + as the per-scenario LP-xbar. The average-scenario solution can + underestimate which first-stage variables need to be active when the + candidate must be feasible across every real scenario, because data + averaging can mask scenario-specific demand patterns. For a + feasibility-oriented candidate, prefer ``lp_xbar_nonants``. + """ + kwargs = scenario_creator_kwargs or {} + avg = average_scenario_creator(scenario_name, **kwargs) + _check_two_stage(avg, "average_xhat_nonants") + if relax_integrality: + pyo.TransformationFactory("core.relax_integer_vars").apply_to(avg) + return _solve_and_extract_root(avg, solver_name, solver_options) + + +def lp_xbar_nonants( + scenario_creator, + scenario_names, + *, + solver_name, + scenario_creator_kwargs=None, + solver_options=None, +): + """For each scenario, solve its LP relaxation; return the + probability-weighted average of ROOT nonants as a 1-D ``np.ndarray``. + + Two-stage only. Each scenario's ``_mpisppy_probability`` is used as + its weight; the literal string ``"uniform"`` is interpreted as + ``1/len(scenario_names)``. Weights are renormalized at the end so a + sub-list of scenarios produces a coherent average. + + This is the per-scenario LP-xbar that downstream consumers + (e.g. findW's pin-dual tests) round to obtain a candidate + first-stage feasible in every real scenario's per-scenario + subproblem. For models where opening more first-stage activity + only loosens recourse, ceiling the LP-xbar opens every variable + that any scenario's LP wanted positive, which makes the candidate + a strict cover. + """ + kwargs = scenario_creator_kwargs or {} + arr_sum = None + weight_sum = 0.0 + n = len(scenario_names) + for sname in scenario_names: + m = scenario_creator(sname, **kwargs) + _check_two_stage(m, "lp_xbar_nonants") + prob = getattr(m, "_mpisppy_probability", "uniform") + if prob == "uniform": + prob = 1.0 / n + pyo.TransformationFactory("core.relax_integer_vars").apply_to(m) + vals = _solve_and_extract_root(m, solver_name, solver_options) + contribution = prob * vals + arr_sum = contribution if arr_sum is None else arr_sum + contribution + weight_sum += prob + return arr_sum / weight_sum diff --git a/run_coverage.bash b/run_coverage.bash index 01d63dfc6..03bf82c4a 100755 --- a/run_coverage.bash +++ b/run_coverage.bash @@ -106,6 +106,9 @@ run_phase "test_generic_cylinders (serial)" \ run_phase "test_jensens (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_jensens.py -v +run_phase "test_feasible_xhat (serial)" \ + coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_feasible_xhat.py -v + run_phase "test_xhat_from_file (serial)" \ coverage run --rcfile=.coveragerc -m pytest mpisppy/tests/test_xhat_from_file.py -v