add feasible_xhat_creator convention for per-scenario-feasible candidates

.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/
@@ -1115,6 +1119,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 \
             mpisppy/tests/test_grad_rho_bundles.py \

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/<model>/<model>_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
+   ``<module_name>_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: ``--<xhatter>-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
+
+   ``--<xhatter>-try-jensens-first`` and
+   ``--<xhatter>-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.

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

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
+   ``--<xhatter>-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
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~