Bootstrap/bagging confidence intervals (boot-sp merge PR-2: statdist + smoothed methods)

.github/workflows/test_pr_and_main.yml

@@ -829,6 +829,25 @@ jobs:
           conda install mpi4py "numpy" setuptools
           pip install pyomo pandas xpress cplex scipy sympy dill packaging coverage
 
+      # the bootstrap epi-spline distributions fit with a nonlinear solver, so
+      # without ipopt those tests skip; this is how Pyomo's own CI gets ipopt
+      - name: Install Ipopt
+        run: |
+          # the ipopt binary links against these; Pyomo's CI installs the same
+          sudo apt-get update -q
+          sudo apt-get install -y libopenblas-dev gfortran liblapack-dev
+          IPOPT_DIR=$HOME/ipopt
+          mkdir -p "$IPOPT_DIR"
+          echo "$IPOPT_DIR" >> $GITHUB_PATH
+          echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$IPOPT_DIR" >> $GITHUB_ENV
+          URL=https://github.com/IDAES/idaes-ext
+          VER=$(curl -sL -H 'Accept: application/json' $URL/releases/latest \
+                | sed -e 's/.*"tag_name":"\([^"]*\)".*/\1/')
+          if test -z "$VER"; then echo "FAILED identifying the ipopt release"; exit 1; fi
+          curl -fL $URL/releases/download/$VER/idaes-solvers-ubuntu2204-x86_64.tar.gz \
+            | tar -xz -C "$IPOPT_DIR"
+          "$IPOPT_DIR"/ipopt -v
+
       - name: setup the program
         run: |
           pip install -e .
@@ -863,13 +882,15 @@ jobs:
           cd mpisppy/tests
           coverage run $COV_ARGS test_boot_sp.py
           coverage run $COV_ARGS test_boot_sp_simulate.py
+          coverage run $COV_ARGS test_boot_sp_smoothed.py
 
       - name: run bootstrap CI tests (mpiexec -np 2)
         timeout-minutes: 10
         run: |
           cd mpisppy/tests
           mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp.py
           mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp_simulate.py
+          mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp_smoothed.py
 
       - name: Upload coverage data
         if: always()

.ruff.toml

@@ -6,3 +6,18 @@ extend-exclude = [
     "./examples/hydro/hydro.py",
     "./examples/sizes/models/ExpressionModel.py",
 ]
+
+[lint.per-file-ignores]
+# The bootsp/statdist/ subpackage is a faithful port of the (legacy) statdist
+# distribution library; relax the style rules it predates rather than rewrite
+# its numerics. splines.py builds a Pyomo model with `from pyomo.environ import *`
+# (like the excluded Pyomo model files above), hence F403/F405.
+"mpisppy/confidence_intervals/bootsp/statdist/*" = [
+    "E711",  # comparison to None
+    "E722",  # bare except
+    "E731",  # lambda assignment
+    "E741",  # ambiguous variable name (math notation, e.g. l)
+    "F403",  # star import (pyomo.environ)
+    "F405",  # name may be from star import
+    "F821",  # undefined name (in retained multivariate base class)
+]

doc/designs/bootsp_merge_design.md

@@ -1,12 +1,16 @@
 # Bootstrap/bagging for data-based stochastic programming in mpi-sppy — design
 
-**Status:** design captured and decisions ratified 2026-07-02; PR-1
-(empirical core + schultz, incl. a data-file example) implemented and open
-upstream as draft [Pyomo/mpi-sppy#783](https://github.com/Pyomo/mpi-sppy/pull/783);
-extended 2026-07-03 to state the end goal
-(`generic_cylinders` integration) and a stacked, multi-PR roadmap (§6, §9).
+**Status:** design captured and decisions ratified 2026-07-02; extended
+2026-07-03 to state the end goal (`generic_cylinders` integration) and a
+stacked, multi-PR roadmap (§6, §9). PR-1 (empirical core + schultz, incl. a
+data-file example) merged 2026-07-24 as
+[Pyomo/mpi-sppy#783](https://github.com/Pyomo/mpi-sppy/pull/783), followed by
+a `test_boot_sp.py` `np=2` fix merged 2026-07-28 as
+[#820](https://github.com/Pyomo/mpi-sppy/pull/820). PR-2 (statdist + smoothed
+methods) is this branch, open upstream as
+[#818](https://github.com/Pyomo/mpi-sppy/pull/818); PR-3 not yet started.
 **Author:** dlw (captured with Claude Code assistance)
-**Last updated:** 2026-07-03
+**Last updated:** 2026-07-28
 
 **Ultimate goal.** The end state this design builds toward is *bootstrap and
 bagging confidence intervals, computed from a given dataset, available
@@ -286,6 +290,33 @@ Behavior-preserving unless noted.
     `boot_method` dies with a raw `TypeError` instead of a friendly
     message. The port accepts real json booleans (keeping the strings
     for boot-sp files) and reports a missing `boot_method` clearly.
+13. **Independent smoothed-bootstrap batches (behavior change, from
+    broken to working).** `smoothed_resample_helper` advanced the record
+    index by one per batch while taking a block of `subsample_size`
+    consecutive records, so consecutive batches shared all but one of
+    their draws — a sliding window, not independent resamples. Chen &
+    Woodruff (2024, Algorithm 3) draws a fresh set of `N` points from the
+    fitted distribution for each of the `B` batches. The port strides by
+    the batch size, as `smoothed_bagging` already did, so the blocks are
+    pairwise disjoint. The estimated spread was badly understated before:
+    on the `cvar` example (`N = 20`, `nB = 10`) the interval widens about
+    fivefold. `simulate_boot` spaces its coverage replications to match:
+    the ported spacing was a hard-coded `nB * 100` unrelated to how many
+    record numbers a replication actually consumes, and is now exactly
+    that footprint. (The empirical harness needs no such spacing at all,
+    because item 11 gave it independent numpy streams; the smoothed path
+    addresses its draws by record number, since that is what the model
+    seeds each draw with, so its replications are separated by giving
+    each one a disjoint block of record numbers.)
+14. **Smoothed center from the fitted distribution (behavior change,
+    from broken to working).** `smoothed_bootstrap` called
+    `center_smoothed` while `use_fitted` was still `False`, so the center
+    was the purely empirical gap. Algorithm 3 takes the center from
+    Algorithm 2 run on the *fitted* distribution, and that smoothed center
+    is the paper's leading conclusion, so `use_fitted` is now set before
+    the center is estimated (`smoothed_bagging` already did this). The
+    center block of the index space is reserved ahead of the batch blocks,
+    since both now sample the same fitted distribution.
 
 ---
 

doc/src/boot_sp.rst

@@ -11,12 +11,14 @@ mpi-sppy, no distribution of the uncertain data is assumed: the estimators
 work directly from sampled data. The methods and software are described in
 [ChenWoodruff2023]_ and [ChenWoodruff2024]_.
 
-.. note::
-
-   This is the empirical (numpy-only) part of the package: the classical,
-   extended, subsampling, and bagging methods. The *smoothed* methods, which
-   depend on a distribution-fitting library, are merged separately; asking for
-   a ``Smoothed_*`` method raises an informative error until then.
+The package has two families of estimators. The *empirical* methods
+(classical, extended, subsampling, and bagging) resample the observed data
+directly and need only numpy. The *smoothed* methods fit a univariate
+distribution to the sampled data (using the bundled ``statdist`` library) and
+resample from the fitted distribution; they need `scipy
+<https://scipy.org>`_, which mpi-sppy treats as an optional dependency and
+imports lazily. If scipy is not installed, the empirical methods still work
+and a smoothed method fails with an informative import error.
 
 Modes
 -----
@@ -60,6 +62,10 @@ plus a few helpers used by the bootstrap code:
   for this fixed name first and falls back to the legacy
   ``xhat_generator_<module_name>``. If a precomputed ``xhat`` file is given
   (``--xhat-fname``) the generator is not called.
+* ``data_sampler(record_num, cfg)`` — return the data for one record (a scalar,
+  or a dict keyed by variable name for multivariate data). This is used by the
+  *smoothed* methods to build the sample that a distribution is fitted to; the
+  empirical methods do not need it.
 
 Methods
 -------
@@ -84,6 +90,21 @@ The ``--boot-method`` (json ``boot_method``) option selects the estimator:
      - Bagging with replacement [lam2018]_
    * - ``Bagging_without_replacement``
      - Bagging without replacement [lam2018]_
+   * - ``Smoothed_boot_epi``
+     - Smoothed bootstrap, epi-spline fit, Gaussian interval [ChenWoodruff2024]_
+   * - ``Smoothed_boot_kernel``
+     - Smoothed bootstrap, kernel-density fit, Gaussian interval [ChenWoodruff2024]_
+   * - ``Smoothed_boot_epi_quantile``
+     - Smoothed bootstrap, epi-spline fit, quantile interval [ChenWoodruff2024]_
+   * - ``Smoothed_boot_kernel_quantile``
+     - Smoothed bootstrap, kernel-density fit, quantile interval [ChenWoodruff2024]_
+   * - ``Smoothed_bagging``
+     - Smoothed bagging, kernel-density fit [ChenWoodruff2024]_
+
+The ``Smoothed_*`` tokens are the smoothed methods; the others are empirical.
+The epi-spline fit builds a small Pyomo nonlinear program, so those two methods
+additionally need a nonlinear solver (e.g. ``ipopt``); the kernel methods do
+not.
 
 Arguments
 ---------
@@ -112,6 +133,14 @@ command line (with dashes). The main options are:
 * ``coverage_replications`` (simulation only) — number of coverage replications.
 * ``boot_method`` / ``--boot-method`` — one of the tokens above.
 
+The smoothed methods use two additional options (ignored, and not required in
+the json, for the empirical methods):
+
+* ``smoothed_center_sample_size`` / ``--smoothed-center-sample-size`` — number
+  of points drawn from the fitted distribution to estimate the gap center.
+* ``smoothed_B_I`` / ``--smoothed-B-I`` — number of outer replications for
+  smoothed bagging.
+
 There may also be model-specific options added by ``inparser_adder``.
 
 Batch parallelism
@@ -190,6 +219,46 @@ from the confidence-interval sampling, so ``sample_size`` plus
 makes it reproducible); replace it with your own two-column dataset, or point
 ``--data-file`` at another file, to run the bootstrap on your own data.
 
+Smoothed methods and statdist
+-----------------------------
+
+The smoothed methods (the ``Smoothed_*`` tokens) fit a univariate distribution
+to the sampled data and then resample from the *fitted* distribution rather
+than from the data directly. The distribution fitting is provided by the
+bundled ``statdist`` library
+(``mpisppy.confidence_intervals.bootsp.statdist``), a trimmed port of the
+univariate distributions from the statdist package; ``statdist`` uses scipy,
+which is imported lazily so that the empirical methods remain scipy-free.
+
+To use a smoothed method the model module must supply ``data_sampler`` (see
+above): the smoothed estimator calls it for each sampled record to assemble the
+data that ``statdist`` fits. The kernel-density methods
+(``Smoothed_boot_kernel``, ``Smoothed_boot_kernel_quantile``,
+``Smoothed_bagging``) fit with a Gaussian kernel and need only scipy; the
+epi-spline methods (``Smoothed_boot_epi``, ``Smoothed_boot_epi_quantile``) fit
+by solving a small Pyomo nonlinear program and additionally need a nonlinear
+solver such as ``ipopt``.
+
+Three examples that need statdist ship in ``examples/bootsp``:
+
+* ``farmer`` — the scalable farmer, with crop yields perturbed by a fitted
+  (or, empirically, a uniform) distribution;
+* ``cvar`` — a CVaR example (Lam & Qian) with standard-normal data;
+* ``multi_knapsack`` — a multi-product knapsack (Vaagen & Wallace) whose
+  deterministic data is read from a json file (``--deterministic-data-json``).
+
+Each has an empirical json/bash and a ``smoothed_*.json``; for instance, from
+``examples/bootsp/cvar``:
+
+.. code-block:: bash
+
+   $ python -m mpisppy.confidence_intervals.bootsp.user_boot cvar \
+       --max-count 3000 --candidate-sample-size 10 --sample-size 75 \
+       --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 0 \
+       --solver-name cplex_direct --boot-method Bagging_with_replacement
+
+   $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_cvar.json
+
 References
 ----------
 

examples/bootsp/cvar/cvar.bash

@@ -0,0 +1,19 @@
+#!/bin/bash
+# Run the CVaR bootstrap example (needs the statdist library).
+# Pass a solver name as the first argument (default: cplex_direct).
+
+SOLVER=${1:-cplex_direct}
+BOOT="python -m mpisppy.confidence_intervals.bootsp.user_boot"
+COMMON="--max-count 3000 --candidate-sample-size 10 --sample-size 75 \
+        --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 0 \
+        --solver-name ${SOLVER}"
+
+echo "Serial, compute xhat within user_boot (empirical Bagging_with_replacement)"
+echo
+time ${BOOT} cvar ${COMMON} --boot-method Bagging_with_replacement
+echo
+echo "========================"
+echo
+echo "Smoothed coverage simulation from a json file (Smoothed_bagging)"
+echo
+time python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_cvar.json

examples/bootsp/cvar/cvar.json

@@ -0,0 +1,16 @@
+{
+    "module_name": "cvar",
+    "max_count": 3000,
+    "candidate_sample_size": 10,
+    "sample_size": 75,
+    "subsample_size": 10,
+    "nB": 20,
+    "alpha": 0.1,
+    "seed_offset": 0,
+    "optimal_fname": "None",
+    "xhat_fname": "None",
+    "solver_name": "cplex_direct",
+    "boot_method": "Bagging_with_replacement",
+    "trace_fname": "None",
+    "coverage_replications": 5
+}