diff --git a/doc/designs/bootsp_merge_design.md b/doc/designs/bootsp_merge_design.md index ae9ee8c65..72cd9fa21 100644 --- a/doc/designs/bootsp_merge_design.md +++ b/doc/designs/bootsp_merge_design.md @@ -8,7 +8,12 @@ extended 2026-07-03 to state the end goal PR-2 (statdist + smoothed) and PR-3 (the `generic_cylinders` integration, the end goal — `--boot-*` group, `do_boot`, the positional layer with a strictly disjoint M/N split, the `K = 1` batch executor) are implemented on the stacked -branches `bootsp-pr-b` / `bootsp-pr-c`. +branches `bootsp-pr-b` / `bootsp-pr-c`. PR-4 (the `K > 1` batch executor: a +`BatchExecutor` that groups the ranks and a wheel per group from +`--boot-batch-config-file`, retiring the interim `--boot-solver-*`) is +implemented on `bootsp-pr-d`; both endpoints are validated end to end — the +`G = 1` checkpoint (one group, a wheel per batch in sequence) and `G > 1` +(np = 4, K = 2). **Author:** dlw (captured with Claude Code assistance) **Last updated:** 2026-07-03 @@ -279,11 +284,13 @@ Behavior-preserving unless noted. pool stream, and the sequential seed offsets used by the coverage simulation reuse batch streams across replications (replication k on rank r has the same stream as replication k+1 on rank r-1). The port - seeds every stream with a `(seed_offset, word)` pair - (`boot_sp._pool_rng` / `_batch_rng`), which numpy's `SeedSequence` - hashes into independent streams; no two streams coincide within a - run or across seed offsets, and `_extended_resample`'s special - `+ my_rank + 1` offset is retired. + seeds every stream with a `(seed_offset, word)` pair, which numpy's + `SeedSequence` hashes into independent streams; no two streams + coincide within a run or across seed offsets, and + `_extended_resample`'s special `+ my_rank + 1` offset is retired. + The pool/center stream is `boot_sp._pool_rng` (word 0) and the batch + streams come from `BatchExecutor.group_seed` (word `group_index + 1`, + which at `K = 1` is the rank, so the per-rank streams are unchanged). 12. **Real json booleans.** boot-sp's `cfg_from_json` accepts only the strings `"True"`/`"False"` for bool options and crashes with an `AttributeError` on a real json `true`/`false`; a json missing @@ -613,10 +620,11 @@ given. The **batch config** governs solve (2). It is *singular across batches* — you resample the data, not re-tune per batch, so it is one config, never a per-batch or `boot_*`-prefixed copy of the whole option surface — but it is **distinct -from the xhat-solve config**: a batch has `N` (or the subsample size) scenarios, -usually far more than the `M` candidate records, so its rho, iteration count, -and spoke mix are a different problem and must be set independently rather than -inherited from the xhat solve. Because it is a *full* `generic_cylinders` +from the xhat-solve config**: a batch is a resample of the data, with its own +scenario count (`N` for the classical and extended methods, the subsample size +for subsampling and bagging) set independently of the `M` candidate records, so +its rho, iteration count, and spoke mix are a different problem and must be set +independently rather than inherited from the xhat solve. Because it is a *full* `generic_cylinders` configuration (solver, rho, which spokes, convergence, relative gap), it is supplied as a **file** — `--boot-batch-config-file` (§9.5), parsed by the same `Config` machinery — not as a growing set of `boot_*`-prefixed CLI flags or an @@ -643,8 +651,22 @@ Note the two endpoints are the same mechanism: `K = 1` is `G = R` (one rank per group, batches spread one-per-rank), and "serial cylinders per batch" is `K = R` → `G = 1` (one group of all `R` ranks, batches in sequence). So PR-4 is a single executor with `--boot-ranks-per-batch` as its only new rank knob, and the -`G = 1` case is a development checkpoint, not a separate PR. The prerequisite -(#782) has landed, so PR-4 is unblocked. +`G = 1` case is a development checkpoint, not a separate PR. + +*Implemented (PR-4).* A `BatchExecutor` +(`mpisppy/confidence_intervals/bootsp/batch_executor.py`) owns the rank +arithmetic and the collectives: it splits `R` ranks into `G = R // K` groups, +exposes each group's communicator (for the wheel / xhat-evaluation) and a +leaders-only communicator (for the cross-group `Gatherv`), and reproduces the +`K = 1` per-rank behavior bit-for-bit as its degenerate case (so the standalone +drivers are unchanged). The empirical estimators route their parallelism through +it. For `K > 1` the per-group wheel is built from `--boot-batch-config-file` +(`mpisppy/generic/boot_batch.py`), and its outer (decomposition) bound is read +back as `L_b`. Validated end to end on the `schultz_data` MIP: the `G = 1` +checkpoint and `G > 1` (np = 4, K = 2) both reproduce the exact value-at-`xhat` +of the `K = 1` EF path (the xhat-evaluation is K-invariant) while their outer +bound sits at or below the EF optimum, so the reported gap is conservative, +exactly as §9.4.1 predicts. ### 9.4.1 The per-batch value: inner bound minus outer bound @@ -683,8 +705,10 @@ the estimators report a gap that the drivers floor at 0 (`ci_gap[0] = max(0, …)`), which would turn a maximization interval into `[0, 0]`, and the coverage harness's one-sided check assumes the same orientation. Per the repo-wide rule that maximization either works or raises, this raises: -`boot_sp._require_minimization` is called from `solve_routine`, which every -batch goes through while every batch is an extensive form. +`boot_sp._require_minimization` is called from `solve_routine` (every extensive +form) and once from `do_boot` on a probe scenario, which is what covers +`K > 1`, where a batch is solved by a wheel and no extensive form is ever built +for `solve_routine` to inspect. Supporting it properly means choosing how the gap is *reported* — mpi-sppy's MMW estimator takes the magnitude, which would keep the gap non-negative in both senses and leave the floors and coverage checks untouched — and threading the diff --git a/doc/src/boot_sp.rst b/doc/src/boot_sp.rst index a6a82025b..1c11e71ab 100644 --- a/doc/src/boot_sp.rst +++ b/doc/src/boot_sp.rst @@ -6,10 +6,11 @@ Bootstrap Confidence Intervals The ``mpisppy.confidence_intervals.bootsp`` subpackage provides bootstrap and bagging confidence intervals for the optimality gap (and for the optimal value and the value at a candidate solution) of *data-based*, two-stage -stochastic programs. Unlike the other confidence-interval methods in -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]_. +stochastic programs. The estimators +work directly from data. The methods and software are described in +[ChenWoodruff2023]_ and [ChenWoodruff2024]_. Although both bootstrap and +bagging are supported, we often just refer to the methods collectively +as "bootstrap" for ease of exposition. The package has two families of estimators. The *empirical* methods (classical, extended, subsampling, and bagging) resample the observed data @@ -20,55 +21,24 @@ resample from the fitted distribution; they need `scipy imports lazily. If scipy is not installed, the empirical methods still work and a smoothed method fails with an informative import error. -Modes ------ - -There are two modes, each runnable with ``python -m``: - -*User mode* (``user_boot``) computes a confidence interval for one problem -instance. A long list of arguments is supplied on the command line, so users -usually put the command in a shell script: - -.. code-block:: bash - - $ python -m mpisppy.confidence_intervals.bootsp.user_boot module arguments - -Here ``module`` is the name of an importable Python module (without ``.py``) -that supplies the scenario creator and helper functions, and ``arguments`` is -the list of double-dash options described below. - -*Simulation mode* (``simulate_boot``) estimates the coverage rate of a method -over many replications; it is aimed at researchers. All options come from a -json file: - -.. code-block:: bash - - $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot instance.json - -The model module ----------------- - -The named module must provide the usual mpi-sppy scenario-creation contract -plus a few helpers used by the bootstrap code: +The rest of this page is organized around how a bootstrap run is put together. +The **Background** section describes the pieces every run shares: the estimator *methods*, +the *contract* a model module must satisfy, and how the *optimality gap* is +defined. There are then two ways to compute an interval — the **standalone +drivers** (``user_boot`` for a single interval, ``simulate_boot`` for coverage +studies), and the everyday **generic_cylinders** driver, which offers the +data-based bootstrap as a first-class option alongside its other +confidence-interval methods. The **smoothed methods** and their ``statdist`` +dependency come last; they are currently available only through the standalone +drivers. + +Background +---------- -* ``scenario_creator(scenario_name, ...)`` — build a Pyomo model for one - (data) scenario, annotated as usual for mpi-sppy; -* ``scenario_names_creator(num_scens, start=None)`` — the list of scenario - names; -* ``kw_creator(cfg)`` — keyword arguments for the scenario creator; -* ``inparser_adder(cfg)`` — add any model-specific options; -* ``xhat_generator(scenario_names, solver_name=None, ...)`` — solve for a - candidate solution ``xhat`` when none is supplied. The bootstrap code looks - for this fixed name first and falls back to the legacy - ``xhat_generator_``. 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. +These pieces are shared by every bootstrap run, whichever driver you use. Methods -------- +~~~~~~~ The ``--boot-method`` (json ``boot_method``) option selects the estimator: @@ -106,11 +76,88 @@ 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. +The model module +~~~~~~~~~~~~~~~~~~ + +The named module must provide the usual mpi-sppy scenario-creation contract +plus a few helpers used by the bootstrap code: + +* ``scenario_creator(scenario_name, ...)`` — build a Pyomo model for one + (data) scenario, annotated as usual for mpi-sppy; +* ``scenario_names_creator(num_scens, start=None)`` — the full list of scenario + names; +* ``kw_creator(cfg)`` — keyword arguments for the scenario creator; +* ``inparser_adder(cfg)`` — add any model-specific options; +* ``xhat_generator(scenario_names, solver_name=None, ...)`` (optional and + seldom provided) — solve for a candidate solution ``xhat`` when none is + supplied. The bootstrap code looks for this fixed name first and falls back + to the legacy ``xhat_generator_``. If a precomputed ``xhat`` + file is given (``--xhat-fname``), the generator is not called. It is also + ignored by ``generic_cylinders``. +* ``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. + +The optimality gap +~~~~~~~~~~~~~~~~~~~ + +The optimality gap of each batch is its value at ``xhat`` minus the batch's +optimal. For the optimal the estimators use the solver's **best bound** (an +outer bound), not the incumbent objective: a mixed-integer batch is solved only +to a MIP gap, so its incumbent would understate the gap. (The solver's +incumbent is used only where no bound is reported.) + +What that buys is a *point* estimate that never understates: the reported gap, +for each batch and for the pool as a whole, is at least the gap an exact solve +would report. The interval **endpoints** carry no such guarantee. Each endpoint +combines the point estimate with a width built from the *spread* of the batch +gaps, and a bound slack perturbs that spread in either direction; the pivotal +methods (``Classical_quantile``, ``Subsampling``, ``Extended``) reflect the +bootstrap quantiles about the point estimate, so a batch slack can move both of +their endpoints *down*. So do not read the reported interval as conservative. + +All of the methods converge to the exact-solve interval as the batch solves +tighten, so when the interval itself matters, control the batch optimality gap +with ``--rel-gap`` in the batch config file (see below) rather than relying on +the outer bound to err in a safe direction. + +Standalone drivers +------------------ + +The legacy standalone drivers compute a bootstrap CI directly, outside the ``generic_cylinders`` +driver. Both standalone drivers are run with ``python -m``. + +Modes +~~~~~ + +There are two modes: + +*User mode* (``user_boot``) computes a confidence interval for one problem +instance. A long list of arguments is supplied on the command line, so users +usually put the command in a shell script: + +.. code-block:: bash + + $ python -m mpisppy.confidence_intervals.bootsp.user_boot module arguments + +Here ``module`` is the name of an importable Python module (without ``.py``) +that supplies the scenario creator and helper functions, and ``arguments`` is +the list of double-dash options described below. + +*Simulation mode* (``simulate_boot``) estimates the coverage rate of a method +over many replications; it is aimed at researchers. All options come from a +json file: + +.. code-block:: bash + + $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot instance.json + Arguments ---------- +~~~~~~~~~ Simulation and user modes use almost the same options; simulation mode reads -them from json (some with underscores), while user mode takes them on the +them from json (some with underscores), while user mode takes them from the command line (with dashes). The main options are: * ``max_count`` / ``--max-count`` — total sample size (integer). @@ -144,36 +191,19 @@ the json, for the empirical methods): There may also be model-specific options added by ``inparser_adder``. Batch parallelism ------------------ +~~~~~~~~~~~~~~~~~~ The bootstrap batches are split across MPI ranks and reassembled on rank 0 with ``Gatherv``, so a run can be accelerated with, e.g., ``mpiexec -np 2 python -m mpi4py -m mpisppy.confidence_intervals.bootsp.user_boot ...``. The estimate on rank 0 depends on the number of ranks because each rank seeds -its own bootstrap stream. - -The optimality gap of each batch is its value at ``xhat`` minus the batch's -optimal. For the optimal the estimators use the solver's **best bound** (an -outer bound), not the incumbent objective: a mixed-integer batch is solved only -to a MIP gap, so its incumbent would understate the gap. (The solver's -incumbent is used only where no bound is reported.) - -What that buys is a *point* estimate that never understates: the reported gap, -for each batch and for the pool as a whole, is at least the gap an exact solve -would report. The interval **endpoints** carry no such guarantee. Each endpoint -combines the point estimate with a width built from the *spread* of the batch -gaps, and a bound slack perturbs that spread in either direction; the pivotal -methods (``Classical_quantile``, ``Subsampling``, ``Extended``) reflect the -bootstrap quantiles about the point estimate, so a batch slack can move both of -their endpoints *down*. So do not read the reported interval as conservative. - -All of the methods converge to the exact-solve interval as the batch solves -tighten, so when the interval itself matters, tighten the batch solves (e.g. a -smaller ``mipgap`` via ``--boot-solver-options``) rather than relying on the -outer bound to err in a safe direction. +its own bootstrap stream. In the standalone drivers each rank solves its own +batches as extensive forms; the ``generic_cylinders`` integration generalizes +this to groups of ``K`` ranks per batch (see `MPI ranks: groups of K`_ below), +of which this is the ``K = 1`` case. boot_general_prep ------------------ +~~~~~~~~~~~~~~~~~~ ``boot_general_prep`` writes the two npy files (a candidate ``xhat`` and a presumed optimal value) that a simulation can reuse: @@ -182,32 +212,15 @@ presumed optimal value) that a simulation can reuse: $ python -m mpisppy.confidence_intervals.bootsp.boot_general_prep instance.json -Example -------- +Examples +~~~~~~~~ -The ``examples/bootsp/schultz`` directory has a small two-stage example whose -data is a deterministic function of the scenario number, so its results are -reproducible across solvers. From that directory: - -.. code-block:: bash +From a dataset file +^^^^^^^^^^^^^^^^^^^^ - $ python -m mpisppy.confidence_intervals.bootsp.user_boot unique_schultz \ - --max-count 50 --candidate-sample-size 1 --sample-size 30 \ - --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 100 \ - --solver-name gurobi_direct --boot-method Classical_quantile - - $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot unique_schultz.json - -See ``examples/bootsp/schultz/schultz.bash`` for a serial run, a parallel run, -and a coverage simulation. - -Working from a dataset file ---------------------------- - -The ``schultz`` example above generates its data arithmetically from the -scenario number. The companion example ``examples/bootsp/schultz_data`` shows -the more typical *data-based* setup: the same model, but each scenario reads -one observation (one row) from a committed dataset, ``schultz_data.csv``: +The typical *data-based* workflow reads each scenario's observation from a +dataset. ``examples/bootsp/schultz_data`` is a small two-stage model wired that +way: each scenario reads one row from a committed dataset, ``schultz_data.csv``: .. code-block:: text @@ -239,34 +252,65 @@ 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. +On-the-fly data +^^^^^^^^^^^^^^^^ + +An unusual example that is useful for testing does not have to read a data file, +``examples/bootsp/schultz`` is the same two-stage model as ``schultz_data`` +above, except each scenario computes its own observation from the scenario +number instead of reading a row — so it is self-contained and needs no dataset. +From that directory: + +.. code-block:: bash + + $ python -m mpisppy.confidence_intervals.bootsp.user_boot unique_schultz \ + --max-count 50 --candidate-sample-size 1 --sample-size 30 \ + --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 100 \ + --solver-name gurobi_direct --boot-method Classical_quantile + + $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot unique_schultz.json + +See ``examples/bootsp/schultz/schultz.bash`` for a serial run, a parallel run, +and a coverage simulation. + In generic_cylinders -------------------- -The everyday driver ``generic_cylinders`` can report a data-based bootstrap -confidence interval as a first-class option, the data-based analog of its -``--mmw-*`` (MMW) group. Given a dataset, it finds a candidate solution -``xhat`` with whatever the command line configured (an extensive form for small -instances, or the PH cylinder system for large ones) and then reports a -bootstrap/bagging CI on the optimality gap of ``xhat``, computed from a part of -the data that is held **strictly disjoint** from the records that produced -``xhat``. That disjointness is a correctness requirement: a gap CI is only -meaningful when estimated on data that did not choose ``xhat``. - -The workflow is a hold-out split over the *positions* in the dataset. The +The standard driver ``generic_cylinders`` (see :ref:`generic_cylinders`) can report a +data-based bootstrap confidence interval as a first-class option, the +data-based analog of its ``--mmw-*`` (MMW) group. Given a dataset, it finds a +candidate solution ``xhat`` with whatever the command line configured (an +extensive form for small instances, or the PH cylinder system for large ones) +and then reports a bootstrap/bagging CI on the optimality gap of ``xhat``, +computed from a part of the data that is held **strictly disjoint** from the +records that produced ``xhat``. That disjointness is a correctness requirement: +a gap CI is only meaningful when estimated on data that did not choose ``xhat``. + +The workflow involves a hold-out split over the *positions* in the dataset. The driver treats ``scenario_names_creator(None)`` as the whole dataset (one scenario name per record), reserves the first ``--boot-candidate-sample-size`` (``M``) records as the candidate block that ``xhat`` came from, and resamples the next ``--boot-sample-size`` (``N``) records — disjoint from the candidate -block — for the CI. Because a model stays name-based while bootstrap resampling +block — for the CI. Since they cannot overlap, ``M + N`` cannot exceed the +dataset size; you normally make them add up to it, using every record — the +first ``M`` to find ``xhat`` and the remaining ``N`` for the CI. Because a model +stays name-based while bootstrap resampling is positional, the driver owns the position/name reconciliation: a model only has to follow the usual mpi-sppy naming and map its own names to its own data. -The dataset is interpreted by the **model**, not the framework: the model owns +The dataset is interpreted by the **model** python module: the model owns loading and any data-source option (e.g. ``--data-file``), and reports the dataset by returning every implied scenario name from ``scenario_names_creator(None)``. There is no dataset-size or data-source ``--boot-*`` flag. +Because it works from a fixed dataset, a bootstrap run is mutually exclusive +with the distribution-sampling CI methods (MMW and sequential sampling). At +present, it is two-stage only. + +Options +~~~~~~~ + The ``--boot-*`` options are: .. list-table:: @@ -284,7 +328,8 @@ The ``--boot-*`` options are: * - ``--boot-sample-size`` - ``N``: records resampled for the CI (the disjoint pool) * - ``--boot-subsample-size`` - - subsample size (subsampling and bagging methods) + - the number of records in each subsample (or bag) drawn from the + ``N``-record pool (subsampling and bagging methods) * - ``--boot-nB`` - number of bootstrap/bagging batches * - ``--boot-alpha`` @@ -292,36 +337,102 @@ The ``--boot-*`` options are: * - ``--boot-seed-offset`` - RNG offset for replication * - ``--boot-xhat-input-file-name`` - - optional precomputed ``xhat`` (the no-wheel path) - * - ``--boot-solver-name`` - - solver for the batch solves (falls back to the generic ``--solver-name``) - * - ``--boot-solver-options`` - - options string for the batch solver, e.g. ``mipgap=0.01`` + - optional precomputed ``xhat`` supplied from a file (not found by a solve) + * - ``--boot-batch-config-file`` + - **required**: a file of ``generic_cylinders`` flags configuring how each + resampled batch is solved (see below) * - ``--boot-ranks-per-batch`` - - ``K``: ranks cooperating on one batch solve; only ``K = 1`` (a per-rank - extensive form) is supported so far + - ``K``: MPI ranks cooperating on one batch solve. ``K = 1`` solves each + batch as a per-rank extensive form; ``K > 1`` solves it with cylinders + on a group of ``K`` ranks. ``K`` must divide the rank count and, for + ``K > 1``, be a multiple of the number of cylinders in the batch config + +The batch config file +~~~~~~~~~~~~~~~~~~~~~~~ + +A batch solve is *different* from the ``xhat`` solve — a batch is a resample of +the data, with its own scenario count (``N`` for the classical and extended +methods, the subsample size for subsampling and bagging) set independently of +the ``M`` candidate records — so its solver, rho, spokes, convergence and +relative gap are configured separately, by role. The ``xhat`` solve uses the +ordinary ``generic_cylinders`` command line; the batch solves are configured +entirely by ``--boot-batch-config-file``, which is literally a file of +``generic_cylinders`` flags (``#`` starts a comment). The two cases below are +alternatives — one file per run, not both together. A ``K = 1`` file (a direct +extensive form) need only name a solver:: + + --solver-name gurobi + +while a ``K > 1`` file is the group's full cylinder configuration — a hub and +one or more spokes, e.g. a PH hub with a Lagrangian bounding spoke:: + + --solver-name gurobi --lagrangian --default-rho 1.0 --max-iterations 50 --max-solver-threads 2 + +The framework injects only the batch scenario set (its count and the positional +sample→record mapping); the file must not set the scenario-formation options +(those are the ``--boot-*`` flags above). A ``K > 1`` batch must yield an +**outer** (decomposition) bound on the batch optimal — via a Lagrangian or +subgradient spoke, or the subgradient hub — since that outer bound is the +batch's optimal ``L_b``. + +MPI ranks: groups of K +~~~~~~~~~~~~~~~~~~~~~~~~ + +A batch is solved by a group of ``K`` ranks (``--boot-ranks-per-batch``). The +bootstrap partitions all ``R`` ranks into ``G = R // K`` groups that run +concurrently, each solving its share of the batches in sequence, with the +results gathered to rank 0. The two endpoints are the same mechanism: ``K = 1`` +is ``G = R`` (each rank its own group, a direct extensive form per batch), and +``K = R`` is ``G = 1`` (one group of all ranks solving the batches in sequence, +each solved with cylinders). The xhat-evaluation solves are spread across the +group's ranks the same way. + +``K`` obeys two divisibility rules: it must divide the rank count ``R`` (so the +groups partition the ranks with none left idle), and for ``K > 1`` it must be a +multiple of the number of cylinders in the batch config (so a group's ``K`` +ranks split evenly among the batch's hub and spokes). A two-cylinder batch — a +PH hub and one bounding spoke, as above — therefore needs an even ``K``. + +Because the ``G`` groups run at the same time, many batch solves are in flight +at once. Cap the threads each solver may take with ``--max-solver-threads`` in +the batch config (``2`` above) so the concurrent solves do not oversubscribe the +cores. -Because it works from a fixed dataset, a bootstrap run is mutually exclusive -with the distribution-sampling CI methods (MMW and sequential sampling) and is -two-stage only. ``examples/bootsp/schultz_data/schultz_data_boot.bash`` runs -the whole workflow end to end: +Example +~~~~~~~ + +``examples/bootsp/schultz_data/schultz_data_boot.bash`` runs the whole workflow +end to end: .. code-block:: bash + $ echo "--solver-name gurobi_direct" > batch_config.txt $ mpiexec -np 3 python -m mpisppy.generic_cylinders \ --module-name schultz_data --num-scens 5 \ --max-iterations 20 --default-rho 1.0 --solver-name gurobi_direct \ - --xhatshuffle --lagrangian \ + --xhatshuffle --lagrangian --max-solver-threads 2 \ --boot-method Classical_quantile \ --boot-candidate-sample-size 5 --boot-sample-size 100 \ --boot-subsample-size 20 --boot-nB 20 --boot-alpha 0.1 \ - --boot-seed-offset 100 + --boot-seed-offset 100 \ + --boot-batch-config-file batch_config.txt Here the main run finds ``xhat`` from the first 5 dataset records and the -bootstrap resamples the next 100 (disjoint) records for the gap CI. +bootstrap resamples the next 100 (disjoint) records for the gap CI. The batch +config file names the solver for the ``K = 1`` batch extensive forms; to solve +each batch with cylinders instead, raise ``--boot-ranks-per-batch`` and give the +group's cylinder configuration in the batch config file. +``examples/bootsp/schultz_data/schultz_data_boot_cylinders.bash`` is a runnable +``K > 1`` demo: 6 ranks find ``xhat`` together, then re-form into two groups of +three that solve the batches concurrently, each batch by a single-cylinder +subgradient solve (configured in ``schultz_wheel_batch.txt``). + +Smoothed methods +---------------- -Smoothed methods and statdist ------------------------------ +At present the smoothed methods are available only through the standalone +drivers (``user_boot`` and ``simulate_boot``); ``generic_cylinders`` offers the +empirical methods only. The smoothed methods (the ``Smoothed_*`` tokens) fit a univariate distribution to the sampled data and then resample from the *fitted* distribution rather diff --git a/examples/bootsp/schultz_data/schultz_data_boot.bash b/examples/bootsp/schultz_data/schultz_data_boot.bash index ce0ca819b..dd63fac27 100644 --- a/examples/bootsp/schultz_data/schultz_data_boot.bash +++ b/examples/bootsp/schultz_data/schultz_data_boot.bash @@ -8,16 +8,26 @@ # N = 100 records (kept strictly disjoint from the candidate records) for the CI. # # Pass a solver name as the first argument (default: gurobi_direct). +# +# How each resampled batch is solved is configured by a separate file of +# generic_cylinders flags, --boot-batch-config-file, because a batch (N +# scenarios) is a different problem from the xhat solve (M candidate records). +# Here K = 1 (--boot-ranks-per-batch defaults to 1), so each batch is a direct +# extensive form and the batch config need only name the solver. SOLVER=${1:-gurobi_direct} +BATCH_CONFIG=$(mktemp) +echo "--solver-name ${SOLVER}" > "${BATCH_CONFIG}" + BOOT="--boot-method Classical_quantile \ --boot-candidate-sample-size 5 \ --boot-sample-size 100 \ --boot-subsample-size 20 \ --boot-nB 20 \ --boot-alpha 0.1 \ - --boot-seed-offset 100" + --boot-seed-offset 100 \ + --boot-batch-config-file ${BATCH_CONFIG}" echo "Find xhat with a PH hub + xhatshuffle/lagrangian spokes, then a bootstrap" echo "CI on its optimality gap from a disjoint part of the dataset." @@ -29,4 +39,7 @@ mpiexec -np 3 python -m mpisppy.generic_cylinders \ --default-rho 1.0 \ --solver-name ${SOLVER} \ --xhatshuffle --lagrangian \ + --max-solver-threads 2 \ ${BOOT} + +rm -f "${BATCH_CONFIG}" diff --git a/examples/bootsp/schultz_data/schultz_data_boot_cylinders.bash b/examples/bootsp/schultz_data/schultz_data_boot_cylinders.bash new file mode 100755 index 000000000..31a4018ca --- /dev/null +++ b/examples/bootsp/schultz_data/schultz_data_boot_cylinders.bash @@ -0,0 +1,47 @@ +#!/bin/bash +# Data-based bootstrap CI where each resampled batch is solved WITH CYLINDERS +# (--boot-ranks-per-batch K > 1). This is the K > 1 generalization of +# schultz_data_boot.bash (which solves each batch as a direct extensive form). +# +# The 6 ranks first find xhat together (PH hub + xhatshuffle + lagrangian over +# the M = 5 candidate records), then re-form into G = 6 // K groups of K ranks; +# each group solves its share of the bootstrap batches with a wheel (configured +# by schultz_wheel_batch.txt) and the per-group results are gathered to rank 0. +# +# Here K = 3, so G = 2 groups of 3 ranks run concurrently -- one batch each at a +# time -- until all the batches are done. +# +# Pass a solver name as the first argument (default: gurobi_direct). + +SOLVER=${1:-gurobi_direct} +HERE=$(dirname "$0") + +# The batch config file names gurobi_direct; substitute the chosen solver into a +# temporary copy so the demo runs with any solver. +BATCH_CONFIG=$(mktemp) +sed "s/gurobi_direct/${SOLVER}/" "${HERE}/schultz_wheel_batch.txt" > "${BATCH_CONFIG}" + +BOOT="--boot-method Classical_quantile \ + --boot-candidate-sample-size 5 \ + --boot-sample-size 20 \ + --boot-subsample-size 8 \ + --boot-nB 12 \ + --boot-alpha 0.1 \ + --boot-seed-offset 100 \ + --boot-ranks-per-batch 3 \ + --boot-batch-config-file ${BATCH_CONFIG}" + +echo "Find xhat with 6 ranks, then solve each bootstrap batch with a 3-rank" +echo "cylinder wheel (2 concurrent groups) for its optimality-gap CI." +echo +mpiexec -np 6 python -m mpisppy.generic_cylinders \ + --module-name schultz_data \ + --num-scens 5 \ + --max-iterations 20 \ + --default-rho 1.0 \ + --solver-name ${SOLVER} \ + --xhatshuffle --lagrangian \ + --max-solver-threads 2 \ + ${BOOT} + +rm -f "${BATCH_CONFIG}" diff --git a/examples/bootsp/schultz_data/schultz_wheel_batch.txt b/examples/bootsp/schultz_data/schultz_wheel_batch.txt new file mode 100644 index 000000000..14519c3bb --- /dev/null +++ b/examples/bootsp/schultz_data/schultz_wheel_batch.txt @@ -0,0 +1,19 @@ +# Batch config for solving each resampled bootstrap batch WITH CYLINDERS (K > 1). +# +# This is a file of generic_cylinders flags (# starts a comment). It configures +# how ONE resampled batch is solved. A batch has N scenarios -- a different, and +# usually larger, problem than the xhat solve over the M candidate records -- so +# its solver, rho, spokes and convergence are set here, independently of the +# main command line. The framework injects only the batch's scenario set. +# +# A K > 1 batch must yield an OUTER (decomposition) bound on the batch optimal +# (that bound is the batch's optimal L_b). The subgradient hub is a single +# cylinder that does exactly that, so it works for any group size K. A PH hub +# plus a Lagrangian or subgradient spoke would work too, when K is a multiple of +# the cylinder count. +--solver-name gurobi_direct +--subgradient-hub +--max-iterations 20 +--default-rho 1.0 +# cap threads per solver so the concurrent batch groups do not oversubscribe the cores +--max-solver-threads 2 diff --git a/mpisppy/confidence_intervals/bootsp/batch_executor.py b/mpisppy/confidence_intervals/bootsp/batch_executor.py new file mode 100644 index 000000000..226a6e008 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/batch_executor.py @@ -0,0 +1,160 @@ +############################################################################### +# 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. +############################################################################### +"""Rank grouping for the bootstrap batch solves (design section 9.4). + +A bootstrap run has two things that want the MPI ranks and pull in opposite +directions: a *solve* wants ranks arranged as a hub plus spokes (a wheel), +while the bootstrap is embarrassingly parallel across its ``nB`` batches and +wants ranks as independent batch workers. The reconciliation is a single knob +``K`` (``--boot-ranks-per-batch``): a batch is solved by a *group* of ``K`` +ranks, and the ``R`` ranks are partitioned into ``G = R // K`` groups that run +concurrently, each solving its share of the batches in sequence, with the +results gathered to global rank 0. + +The two endpoints are the same mechanism: + +* ``K = 1`` -> ``G = R``: every rank is its own group and a batch is a direct + extensive form. This reproduces the original standalone behavior exactly, so + ``user_boot`` / ``simulate_boot`` are unaffected. +* ``K = R`` -> ``G = 1``: one group of all ranks solves the batches in sequence, + each by a full wheel on the whole communicator. This is the development + checkpoint for the ``K > 1`` cylinders path. + +This module owns only the rank arithmetic and the collectives. The two solves a +batch needs -- the optimal (outer-bound) solve and the xhat evaluation -- are +supplied by the caller: the xhat evaluation is just an ``Xhat_Eval`` on the +group communicator (scenarios spread across the group's ranks), and the optimal +solve is a direct EF for ``K = 1`` or a wheel (via ``batch_optimal_solver``) for +``K > 1``. +""" + +import numpy as np + +import mpisppy.MPI as MPI + + +def slice_lens_for(nB, nslices): + """Split ``nB`` items into ``nslices`` contiguous shares as evenly as possible. + + Returns a list of length ``nslices`` summing to ``nB`` (the historical + ``boot_sp.slice_lens`` split, generalized to an arbitrary slice count so it + can apportion batches over groups instead of over individual ranks). + """ + avg = nB / nslices + lens = [int((i + 1) * avg) - int(i * avg) for i in range(nslices)] + assert sum(lens) == nB + return lens + + +class BatchExecutor: + """Partition ``comm`` into ``G = size // K`` groups of ``K`` ranks. + + Attributes (all valid on every rank unless noted): + K (int): ranks per group (the group size). + n_groups (int): G, the number of groups. + group_index (int): this rank's group, 0..G-1. + group_rank (int): this rank's position within its group, 0..K-1. + is_group_leader (bool): group_rank == 0 (the rank that reports the + group's batch results into the cross-group gather). + is_root (bool): global rank 0 -- where the estimator does its analysis. + groupcomm (MPI comm): the K ranks of this group; the communicator a + batch solve / xhat evaluation runs on. + leadercomm (MPI comm): the G group leaders (MPI.COMM_NULL on non-leaders + under real MPI); the communicator the batch results are gathered on. + batch_optimal_solver (callable or None): for K > 1, a callable + ``(scenario_names, sample_mapping, groupcomm) -> outer_bound`` that + runs a wheel on the group; None selects the K = 1 direct-EF path. + """ + + def __init__(self, K=1, comm=None, batch_optimal_solver=None): + if comm is None: + comm = MPI.COMM_WORLD + self.comm = comm + R = comm.Get_size() + rank = comm.Get_rank() + + if K is None or K < 1: + raise ValueError(f"boot_ranks_per_batch (K) must be a positive integer, got {K}") + if K > R: + raise ValueError( + f"boot_ranks_per_batch (K={K}) exceeds the number of MPI ranks " + f"({R}); a group cannot have more ranks than the world.") + # Leftover-rank handling (the R mod K ranks that would sit out) is a + # scheduled refinement; for now the partition must be exact so no rank + # is silently idle. The two checkpoints (K=1 and K=R) both divide. + if R % K != 0: + raise ValueError( + f"boot_ranks_per_batch (K={K}) must divide the number of MPI " + f"ranks ({R}) evenly; got remainder {R % K}. Choose K in " + f"{[k for k in range(1, R + 1) if R % k == 0]}.") + + self.K = K + self.n_groups = R // K + self.group_index = rank // K + self.group_rank = rank % K + self.is_group_leader = (self.group_rank == 0) + self.is_root = (rank == 0) + self.batch_optimal_solver = batch_optimal_solver + + # The group communicator: K ranks cooperate on one batch here. At K=1 + # this is a private single-rank comm (the old rankcomm); at K=R it is + # the whole world (the G=1 checkpoint). + self.groupcomm = comm.Split(color=self.group_index, key=self.group_rank) + + # The leader communicator: one rank per group, used only to gather the + # per-group batch results. Non-leaders are excluded (COMM_NULL under + # real MPI) and never touch it. At K=1 every rank is a leader, so this + # is the whole world -- exactly the old COMM_WORLD gather. + leader_color = 0 if self.is_group_leader else MPI.UNDEFINED + self.leadercomm = comm.Split(color=leader_color, key=self.group_index) + + @property + def uses_cylinders(self): + """True when a batch is solved by a wheel (K > 1) rather than an EF.""" + return self.K > 1 + + def batch_share(self, nB): + """This group's number of batches (all ranks in a group agree).""" + return slice_lens_for(nB, self.n_groups)[self.group_index] + + def group_seed(self, seed_offset): + """RNG seed (a SeedSequence pair) for this group's batch stream. + + All ranks in a group share it so they resample the *same* batches (and + then cooperate on solving them); distinct groups get distinct seeds so + they cover different batches. The seed is the ``[seed_offset, word]`` + pair numpy's SeedSequence hashes, with ``word = group_index + 1`` so the + batch stream never coincides with the pool/center stream (word 0; see + ``boot_sp._pool_rng``) or across seed_offsets -- summing the two instead + can collide. At K=1 the group index is the rank, so this is + ``[seed_offset, rank + 1]``, exactly the original per-rank batch stream + (the former ``boot_sp._batch_rng``). + """ + return [seed_offset, self.group_index + 1] + + def gather(self, local, nB, item_len=1): + """Gather each group's local results into the full array. + + ``local`` is this group's ``batch_share(nB) * item_len`` results + (identical on all ranks of the group; only the leader contributes). + ``item_len`` is the number of floats each batch contributes -- 1 for a + scalar per batch, or e.g. ``sample_size`` for the bagging count rows. + Returns the assembled length-``nB * item_len`` array on global rank 0 + and None elsewhere. + """ + local = np.ascontiguousarray(local, dtype=np.float64) + if not self.is_group_leader: + return None + lenlist = [ell * item_len for ell in slice_lens_for(nB, self.n_groups)] + if self.is_root: + full = np.empty(nB * item_len, dtype=np.float64) + else: + full = None + self.leadercomm.Gatherv(sendbuf=local, recvbuf=(full, lenlist), root=0) + return full diff --git a/mpisppy/confidence_intervals/bootsp/boot_sp.py b/mpisppy/confidence_intervals/bootsp/boot_sp.py index 740380ae8..2db4f0fa1 100644 --- a/mpisppy/confidence_intervals/bootsp/boot_sp.py +++ b/mpisppy/confidence_intervals/bootsp/boot_sp.py @@ -21,8 +21,14 @@ import mpisppy.utils.sputils as sputils import mpisppy.utils.xhat_eval as xhat_eval import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils +from mpisppy.confidence_intervals.bootsp.batch_executor import BatchExecutor # The communicators live in boot_utils so there is a single source of truth. +# These module globals are the K = 1 view (each rank its own batch worker) and +# are still used by the smoothed methods and the simulation prep, which are +# K = 1 only. The empirical estimators below instead route their parallelism +# through a BatchExecutor so they generalize to K > 1 (a wheel per batch); at +# K = 1 the executor reproduces exactly what these globals describe. comm = boot_utils.comm n_proc = boot_utils.n_proc my_rank = boot_utils.my_rank @@ -44,8 +50,10 @@ def _require_minimization(is_minimizing, what): """Refuse a maximization model instead of reporting a wrong interval. Per the repo-wide rule that maximization either works or raises, this is - the raise. It is checked in solve_routine, which every extensive form goes - through. + the raise. It is checked in solve_routine, which every extensive-form path + goes through, and once in generic_cylinders' do_boot, which also covers + K > 1 -- that path solves each batch with a wheel and so never builds an + extensive form for solve_routine to inspect. """ if not is_minimizing: raise ValueError(_MAXIMIZATION_MSG.format(what=what)) @@ -76,6 +84,11 @@ def _best_bound(ef, results): objective value. That fallback makes the reported optimality gap read optimistically (design 9.4.1: the incumbent is an inner bound), which is the wrong direction, so warn rather than let it pass silently. + + The K > 1 sibling, boot_batch._outer_bound_over_group, raises on the same + condition instead of falling back. The asymmetry is deliberate: an + extensive-form solve has an incumbent to fall back to, while a wheel run for + its outer bound has nothing to offer in its place. """ try: prob = results.problem[0] @@ -116,6 +129,43 @@ def _ef_optimal_value(ef): return pyo.value(ef.EF_Obj) if val is None else val +def _sample_names_and_mapping(cfg, scenarios, duplication): + """Build the (distinct) sample scenario names and their position mapping. + + A resample can select the same record more than once, but an mpi-sppy + extensive form needs *distinct* scenario names, so with ``duplication`` we + mint fresh ``SampleScenario{i}`` names and map each back to its record's + canonical name (via the positional resolver). Without duplication the + scenarios are distinct positions and can be named directly. Shared by the + optimal solve and the xhat evaluation so the two always agree. + """ + name_of_pos = _name_of_position_fn(cfg) + if duplication: + names = ['SampleScenario' + str(i) for i in range(len(scenarios))] + mapping = {names[i]: name_of_pos(scenarios[i]) for i in range(len(scenarios))} + else: + names = [name_of_pos(s) for s in scenarios] + mapping = None + return names, mapping + + +def _batch_optimal_value(cfg, module, scenarios, executor, duplication): + """The batch optimal (outer bound) ``L_b`` for one resampled batch. + + For K = 1 this is a direct extensive form solved to its best (outer) bound + (``_ef_optimal_value``). For K > 1 the batch is solved by a wheel on the + group's communicator via ``executor.batch_optimal_solver``, which returns + the wheel's decomposition (outer) bound. Both play the same role in the + optimality-gap estimators (design section 9.4.1), so the estimator code is + identical across K. + """ + if executor.uses_cylinders: + names, mapping = _sample_names_and_mapping(cfg, scenarios, duplication) + return executor.batch_optimal_solver(names, mapping, executor.groupcomm) + ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=duplication) + return _ef_optimal_value(ef) + + def _scenario_creator_w_mapping(scenario_name, module=None, mapping=None, **kwargs): """ A wrapper to allow for bootstrap samples to map to actual samples Args: @@ -185,11 +235,6 @@ def _pool_rng(cfg): return default_rng([cfg.seed_offset, 0]) -def _batch_rng(cfg): - """ The per-rank random stream for batch resampling (see _pool_rng). """ - return default_rng([cfg.seed_offset, my_rank + 1]) - - def process_optimal(cfg, module): """ For simulations we need a known or assumed z* Args: @@ -237,13 +282,8 @@ def solve_routine(cfg, module, scenarios, num_threads=None, duplication=False): scenario_creator_kwargs = module.kw_creator(cfg) # we get a new one every time... scenario_creator_kwargs['module'] = module # we are going to call a wrapper - name_of_pos = _name_of_position_fn(cfg) - if duplication: - scenario_names = ['SampleScenario' + str(i) for i in range(len(scenarios))] - scenario_creator_kwargs['mapping'] = {'SampleScenario' + str(i): name_of_pos(scenarios[i]) for i in range(len(scenarios))} - else: - scenario_names = [name_of_pos(s) for s in scenarios] - scenario_creator_kwargs['mapping'] = None + scenario_names, scenario_creator_kwargs['mapping'] = \ + _sample_names_and_mapping(cfg, scenarios, duplication) ef = sputils.create_EF( scenario_names, @@ -274,7 +314,7 @@ def solve_routine(cfg, module, scenarios, num_threads=None, duplication=False): return ef -def evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping): +def evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping, mpicomm=None): """ evaluate a given xhat over given scenario names Args: @@ -285,10 +325,16 @@ def evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping): scenario_names (list of str): the scenario number will be peeled off the ends sample_mapping (dict): If not None, maps the scenario_name argument to a scenario sent to the module scenario creator + mpicomm (MPI comm): the communicator to spread the evaluation scenarios + over; defaults to the private single-rank rankcomm (the K = 1 view). + The empirical estimators pass the batch's group communicator so an + xhat evaluation is spread across the group's ranks (design 9.4). Returns: zhat (float): the computed expected value """ + if mpicomm is None: + mpicomm = rankcomm # optional batch-solver options (see solve_routine); None for the standalone # drivers, so the xhat-evaluation solves are unchanged there solver_options = cfg.get("solver_options", None) @@ -308,7 +354,7 @@ def evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping): ev = xhat_eval.Xhat_Eval(xhat_eval_options, scenario_names, scenario_creator, - mpicomm=rankcomm, + mpicomm=mpicomm, scenario_creator_kwargs=scenario_creator_kwargs ) @@ -317,7 +363,7 @@ def evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping): return zhat -def evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True): +def evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True, mpicomm=None): """ evaluate xhat using a list of (sampled) scenario numbers Args: @@ -336,18 +382,12 @@ def evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True): # If need mapping, create a set of scenario names and a mapping function that maps the scenario names to the original ones # Return the function value evaluated for a given xhat - name_of_pos = _name_of_position_fn(cfg) - if duplication: - scenario_names = ['SampleScenario' + str(i) for i in range(len(scenarios))] - sample_mapping = {'SampleScenario' + str(i): name_of_pos(scenarios[i]) for i in range(len(scenarios))} - else: - scenario_names = [name_of_pos(s) for s in scenarios] - sample_mapping = None + scenario_names, sample_mapping = _sample_names_and_mapping(cfg, scenarios, duplication) - return evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping) + return evaluate_routine(cfg, module, xhat, scenario_names, sample_mapping, mpicomm=mpicomm) -def _bootstrap_resample(cfg, module, scenario_pool, xhat, serial=False): +def _bootstrap_resample(cfg, module, scenario_pool, xhat, executor): """ Get gaps and optimal values for classic bootstrap. Args: cfg (Config): parameters @@ -355,33 +395,31 @@ def _bootstrap_resample(cfg, module, scenario_pool, xhat, serial=False): scenario_pool (iterable; e.g., list): scenario numbers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) - serial (bool): indicates that only one MPI rank should be used + executor (BatchExecutor): the rank grouping (design 9.4) Returns: - numpy arrays (vector) with gaps and optimal values that are *local* if serial is False + numpy arrays (vector) with this group's gaps, optimal values, and uppers """ - # loop over batches + # loop over this group's share of the batches; all ranks of a group share + # the seed (so they resample the same batch) and cooperate on the solve. - rng = _batch_rng(cfg) - if serial: - local_nB = cfg.nB - else: - local_nB = slice_lens(cfg.nB)[my_rank] + rng = default_rng(executor.group_seed(cfg.seed_offset)) + local_nB = executor.batch_share(cfg.nB) local_boot_gaps = np.empty(local_nB, dtype=np.float64) local_boot_optimals = np.empty(local_nB, dtype=np.float64) local_boot_uppers = np.empty(local_nB, dtype=np.float64) for iter in range(local_nB): scenarios = rng.choice(scenario_pool, size=cfg.sample_size, replace=True) - boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True) - boot_ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=True) - local_boot_optimals[iter] = _ef_optimal_value(boot_ef) + boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True, + mpicomm=executor.groupcomm) + local_boot_optimals[iter] = _batch_optimal_value(cfg, module, scenarios, executor, duplication=True) local_boot_uppers[iter] = boot_ev local_boot_gaps[iter] = local_boot_uppers[iter] - local_boot_optimals[iter] return local_boot_gaps, local_boot_optimals, local_boot_uppers -def classical_bootstrap(cfg, module, xhat, quantile=True): +def classical_bootstrap(cfg, module, xhat, quantile=True, executor=None): """ perform a classic bootstrap estimation of confidence intervals Args: @@ -390,49 +428,40 @@ def classical_bootstrap(cfg, module, xhat, quantile=True): xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) quantile (bool): use the quantile method (else the gaussian method) + executor (BatchExecutor or None): the rank grouping; None => K = 1 Returns: tuple with confidence interval if on MPI rank 0 """ + if executor is None: + executor = BatchExecutor(1) rng = _pool_rng(cfg) scenario_pool = rng.choice(eligible_scenarios(cfg), size=cfg.sample_size, replace=False) - dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False) - dag_ef = solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=False) - - dag_optimal = _ef_optimal_value(dag_ef) + dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False, + mpicomm=executor.groupcomm) + dag_optimal = _batch_optimal_value(cfg, module, scenario_pool, executor, duplication=False) dag_gap = dag_upper - dag_optimal # this is gamma(D) in the note # tron is a "secret" way to turn on internal trace information if cfg.get("tron", False): print(f"rank {my_rank} at dag barrier", flush=True) - comm.Barrier() + executor.comm.Barrier() # bootstrap from pool - local_boot_gaps, local_boot_optimals, local_boot_uppers = _bootstrap_resample(cfg, module, scenario_pool, xhat, serial=False) + local_boot_gaps, local_boot_optimals, local_boot_uppers = _bootstrap_resample(cfg, module, scenario_pool, xhat, executor) - comm.Barrier() + executor.comm.Barrier() - # do analysis only on rank 0 - if my_rank == 0: - boot_gaps = np.empty(cfg.nB, dtype=np.float64) - boot_optimals = np.empty(cfg.nB, dtype=np.float64) - boot_uppers = np.empty(cfg.nB, dtype=np.float64) - else: - boot_gaps = None - boot_optimals = None - boot_uppers = None - - # but everyone needs to send to the gather - lenlist = slice_lens(cfg.nB) - comm.Gatherv(sendbuf=local_boot_gaps, recvbuf=(boot_gaps, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_optimals, recvbuf=(boot_optimals, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_uppers, recvbuf=(boot_uppers, lenlist), root=0) - if cfg.get("tron", False) and my_rank == 0: + # gather every group's batches to global rank 0 for analysis + boot_gaps = executor.gather(local_boot_gaps, cfg.nB) + boot_optimals = executor.gather(local_boot_optimals, cfg.nB) + boot_uppers = executor.gather(local_boot_uppers, cfg.nB) + if cfg.get("tron", False) and executor.is_root: print("*** rank 0 ends gather", flush=True) - if my_rank == 0: + if executor.is_root: if quantile: alpha = cfg.alpha / 2 ci_optimal = np.quantile(2 * dag_optimal - boot_optimals, [alpha, 1 - alpha]) @@ -453,7 +482,7 @@ def classical_bootstrap(cfg, module, xhat, quantile=True): return None, None, None, None, None, None -def _sub_resample(cfg, module, scenario_pool, xhat, serial=False): +def _sub_resample(cfg, module, scenario_pool, xhat, executor): """ Get gaps and optimal values for subsampling method. Args: cfg (Config): parameters @@ -461,37 +490,31 @@ def _sub_resample(cfg, module, scenario_pool, xhat, serial=False): scenario_pool (iterable; e.g., list): scenario numbers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) - serial (bool): indicates that only one MPI rank should be used + executor (BatchExecutor): the rank grouping (design 9.4) Returns: - numpy arrays (vector) with gaps and optimal values that are *local* if serial is False + numpy arrays (vector) with this group's gaps, optimal values, and uppers """ - # loop over batches - # only difference between this and bootstrap_sampling is the size of the scenarios: one is subsample_size, the other is sample_size + # loop over this group's share of the batches + # only difference from _bootstrap_resample is the batch size: subsample_size here, sample_size there - rng = _batch_rng(cfg) - if serial: - local_nB = cfg.nB - else: - local_nB = slice_lens(cfg.nB)[my_rank] + rng = default_rng(executor.group_seed(cfg.seed_offset)) + local_nB = executor.batch_share(cfg.nB) local_boot_gaps = np.empty(local_nB, dtype=np.float64) local_boot_optimals = np.empty(local_nB, dtype=np.float64) local_boot_uppers = np.empty(local_nB, dtype=np.float64) for iter in range(local_nB): scenarios = rng.choice(scenario_pool, size=cfg.subsample_size, replace=False) - boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True) - boot_ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=True) - if cfg.get("tron", False) and my_rank == 0: - print(f"_sub_resample using EF_obj: {pyo.value(boot_ef.EF_Obj)}") - print(f" using evaluation: {boot_ev}") - local_boot_optimals[iter] = _ef_optimal_value(boot_ef) + boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True, + mpicomm=executor.groupcomm) + local_boot_optimals[iter] = _batch_optimal_value(cfg, module, scenarios, executor, duplication=True) local_boot_uppers[iter] = boot_ev local_boot_gaps[iter] = local_boot_uppers[iter] - local_boot_optimals[iter] return local_boot_gaps, local_boot_optimals, local_boot_uppers -def subsampling(cfg, module, xhat): +def subsampling(cfg, module, xhat, executor=None): """ perform a subsampling estimation of confidence intervals Args: @@ -499,42 +522,34 @@ def subsampling(cfg, module, xhat): module (Python module): contains the scenario creator function and helpers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) + executor (BatchExecutor or None): the rank grouping; None => K = 1 Returns: tuple with confidence interval if on MPI rank 0 """ + if executor is None: + executor = BatchExecutor(1) rng = _pool_rng(cfg) scenario_pool = rng.choice(eligible_scenarios(cfg), size=cfg.sample_size, replace=False) - dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False) - dag_ef = solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=False) - dag_optimal = _ef_optimal_value(dag_ef) + dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False, + mpicomm=executor.groupcomm) + dag_optimal = _batch_optimal_value(cfg, module, scenario_pool, executor, duplication=False) dag_gap = dag_upper - dag_optimal # this is gamma(D) in the note - comm.Barrier() + executor.comm.Barrier() # subsampling from pool - local_boot_gaps, local_boot_optimals, local_boot_uppers = _sub_resample(cfg, module, scenario_pool, xhat, serial=False) - comm.Barrier() - - # do analysis only on rank 0 - if my_rank == 0: - boot_gaps = np.empty(cfg.nB, dtype=np.float64) - boot_optimals = np.empty(cfg.nB, dtype=np.float64) - boot_uppers = np.empty(cfg.nB, dtype=np.float64) - else: - boot_gaps = None - boot_optimals = None - boot_uppers = None + local_boot_gaps, local_boot_optimals, local_boot_uppers = _sub_resample(cfg, module, scenario_pool, xhat, executor) + executor.comm.Barrier() - # but everyone needs to send to the gather - lenlist = slice_lens(cfg.nB) - comm.Gatherv(sendbuf=local_boot_gaps, recvbuf=(boot_gaps, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_optimals, recvbuf=(boot_optimals, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_uppers, recvbuf=(boot_uppers, lenlist), root=0) + # gather every group's batches to global rank 0 for analysis + boot_gaps = executor.gather(local_boot_gaps, cfg.nB) + boot_optimals = executor.gather(local_boot_optimals, cfg.nB) + boot_uppers = executor.gather(local_boot_uppers, cfg.nB) - if my_rank == 0: + if executor.is_root: alpha = cfg.alpha / 2 err_optimal = np.sqrt(cfg.subsample_size / cfg.sample_size) * np.quantile(boot_optimals - dag_optimal, [1 - alpha, alpha]) ci_optimal = dag_optimal - err_optimal @@ -550,25 +565,22 @@ def subsampling(cfg, module, xhat): return None, None, None, None, None, None -def _extended_resample(cfg, module, xhat, serial=False): +def _extended_resample(cfg, module, xhat, executor): """ Get gaps and optimal values differences for extended bootstrap. Args: cfg (Config): parameters module (Python module): contains the scenario creator function and helpers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) - serial (bool): indicates that only one MPI rank should be used + executor (BatchExecutor): the rank grouping (design 9.4) Returns: - numpy arrays (vector) with gaps and optimal values differences that are *local* if serial is False + numpy arrays (vector) with this group's gap / optimal / upper differences """ - # loop over batches + # loop over this group's share of the batches - rng = _batch_rng(cfg) - if serial: - local_nB = cfg.nB - else: - local_nB = slice_lens(cfg.nB)[my_rank] + rng = default_rng(executor.group_seed(cfg.seed_offset)) + local_nB = executor.batch_share(cfg.nB) local_boot_optimals_diff = np.empty(local_nB, dtype=np.float64) local_boot_uppers_diff = np.empty(local_nB, dtype=np.float64) @@ -577,14 +589,16 @@ def _extended_resample(cfg, module, xhat, serial=False): eligible = eligible_scenarios(cfg) for iter in range(local_nB): scenario_pool = rng.choice(eligible, size=cfg.sample_size, replace=True) - dag_optimal_ef = solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=True) - dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=True) + dag_optimal = _batch_optimal_value(cfg, module, scenario_pool, executor, duplication=True) + dag_upper = evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=True, + mpicomm=executor.groupcomm) scenarios = rng.choice(scenario_pool, size=cfg.sample_size, replace=True) - boot_optimal_ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=True) - boot_upper = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True) + boot_optimal = _batch_optimal_value(cfg, module, scenarios, executor, duplication=True) + boot_upper = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True, + mpicomm=executor.groupcomm) - local_boot_optimals_diff[iter] = _ef_optimal_value(boot_optimal_ef) - _ef_optimal_value(dag_optimal_ef) + local_boot_optimals_diff[iter] = boot_optimal - dag_optimal local_boot_uppers_diff[iter] = boot_upper - dag_upper local_boot_gaps_diff[iter] = local_boot_uppers_diff[iter] - local_boot_optimals_diff[iter] @@ -592,7 +606,7 @@ def _extended_resample(cfg, module, xhat, serial=False): return local_boot_gaps_diff, local_boot_optimals_diff, local_boot_uppers_diff -def extended_bootstrap(cfg, module, xhat): +def extended_bootstrap(cfg, module, xhat, executor=None): """ perform an extended bootstrap estimation of confidence intervals Args: @@ -600,64 +614,62 @@ def extended_bootstrap(cfg, module, xhat): module (Python module): contains the scenario creator function and helpers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) + executor (BatchExecutor or None): the rank grouping; None => K = 1 Returns: tuple with confidence interval if on MPI rank 0 """ + if executor is None: + executor = BatchExecutor(1) rng = _pool_rng(cfg) # extended bootstrap - local_boot_gaps_diff, local_boot_optimals_diff, local_boot_uppers_diff = _extended_resample(cfg, module, xhat, serial=False) - comm.Barrier() + local_boot_gaps_diff, local_boot_optimals_diff, local_boot_uppers_diff = _extended_resample(cfg, module, xhat, executor) + executor.comm.Barrier() - # do analysis only on rank 0 - if my_rank == 0: - boot_gaps_diff = np.empty(cfg.nB, dtype=np.float64) - boot_optimals_diff = np.empty(cfg.nB, dtype=np.float64) - boot_uppers_diff = np.empty(cfg.nB, dtype=np.float64) - else: - boot_gaps_diff = None - boot_optimals_diff = None - boot_uppers_diff = None - - # but everyone needs to send to the gather - lenlist = slice_lens(cfg.nB) - comm.Gatherv(sendbuf=local_boot_gaps_diff, recvbuf=(boot_gaps_diff, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_optimals_diff, recvbuf=(boot_optimals_diff, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_uppers_diff, recvbuf=(boot_uppers_diff, lenlist), root=0) + # gather every group's batch differences to global rank 0 for analysis + boot_gaps_diff = executor.gather(local_boot_gaps_diff, cfg.nB) + boot_optimals_diff = executor.gather(local_boot_optimals_diff, cfg.nB) + boot_uppers_diff = executor.gather(local_boot_uppers_diff, cfg.nB) - if my_rank == 0: + # The center is a fresh set of solves; unlike the per-batch differences it + # is computed once, on group 0 (all of whose ranks must cooperate on the + # solves when K > 1). At K = 1 group 0 is exactly global rank 0, so this is + # the original "on rank 0" behavior. The final CI needs the gathered + # differences, which land only on global rank 0 (is_root, itself in group 0). + if executor.group_index == 0: # get center eligible = eligible_scenarios(cfg) scenarios = rng.choice(eligible, size=cfg.sample_size, replace=True) - dag_optimal_ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=True) - dag_optimal = _ef_optimal_value(dag_optimal_ef) - dag_upper = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True) + dag_optimal = _batch_optimal_value(cfg, module, scenarios, executor, duplication=True) + dag_upper = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=True, + mpicomm=executor.groupcomm) scenarios_ = rng.choice(eligible, size=cfg.sample_size, replace=True) scenarios_combined = np.concatenate([scenarios, scenarios_]) - dag_optimal_ef_combined = solve_routine(cfg, module, scenarios_combined, num_threads=2, duplication=True) - dag_optimal_combined = _ef_optimal_value(dag_optimal_ef_combined) - dag_upper_combined = evaluate_scenarios(cfg, module, scenarios_combined, xhat, duplication=True) + dag_optimal_combined = _batch_optimal_value(cfg, module, scenarios_combined, executor, duplication=True) + dag_upper_combined = evaluate_scenarios(cfg, module, scenarios_combined, xhat, duplication=True, + mpicomm=executor.groupcomm) center_optimal = 2 * dag_optimal_combined - dag_optimal center_upper = 2 * dag_upper_combined - dag_upper center_gap = center_upper - center_optimal - alpha = cfg.alpha / 2 - ci_optimal = center_optimal - np.quantile(boot_optimals_diff, [1 - alpha, alpha]) - ci_upper = center_upper - np.quantile(boot_uppers_diff, [1 - alpha, alpha]) - ci_gap = center_gap - np.quantile(boot_gaps_diff, [1 - alpha, alpha]) + if executor.is_root: + alpha = cfg.alpha / 2 + ci_optimal = center_optimal - np.quantile(boot_optimals_diff, [1 - alpha, alpha]) + ci_upper = center_upper - np.quantile(boot_uppers_diff, [1 - alpha, alpha]) + ci_gap = center_gap - np.quantile(boot_gaps_diff, [1 - alpha, alpha]) - return ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap - else: - return None, None, None, None, None, None + return ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap + + return None, None, None, None, None, None -def _bagging_resample(cfg, module, scenario_pool, xhat, serial=False, replacement=True): +def _bagging_resample(cfg, module, scenario_pool, xhat, executor, replacement=True): """ Get gaps and optimal values differences for bagging. Args: cfg (Config): parameters @@ -665,18 +677,16 @@ def _bagging_resample(cfg, module, scenario_pool, xhat, serial=False, replacemen scenario_pool (iterable; e.g., list): scenario numbers xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) - serial (bool): indicates that only one MPI rank should be used + executor (BatchExecutor): the rank grouping (design 9.4) + replacement (bool): sample the subsample with replacement Returns: - numpy arrays (vector) with gaps, optimal values, and boot counts that are *local* if serial is False + numpy arrays (vector) with this group's gaps, optimal values, uppers, and boot counts """ - # loop over batches + # loop over this group's share of the batches - rng = _batch_rng(cfg) - if serial: - local_nB = cfg.nB - else: - local_nB = slice_lens(cfg.nB)[my_rank] + rng = default_rng(executor.group_seed(cfg.seed_offset)) + local_nB = executor.batch_share(cfg.nB) local_boot_gaps = np.empty(local_nB, dtype=np.float64) local_boot_optimals = np.empty(local_nB, dtype=np.float64) local_boot_uppers = np.empty(local_nB, dtype=np.float64) @@ -684,10 +694,9 @@ def _bagging_resample(cfg, module, scenario_pool, xhat, serial=False, replacemen for iter in range(local_nB): scenarios_index = rng.choice(len(scenario_pool), size=cfg.subsample_size, replace=replacement) scenarios = [scenario_pool[index] for index in scenarios_index] - boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=replacement) - boot_ef = solve_routine(cfg, module, scenarios, num_threads=2, duplication=replacement) - - local_boot_optimals[iter] = _ef_optimal_value(boot_ef) + boot_ev = evaluate_scenarios(cfg, module, scenarios, xhat, duplication=replacement, + mpicomm=executor.groupcomm) + local_boot_optimals[iter] = _batch_optimal_value(cfg, module, scenarios, executor, duplication=replacement) local_boot_uppers[iter] = boot_ev local_boot_gaps[iter] = local_boot_uppers[iter] - local_boot_optimals[iter] @@ -699,7 +708,7 @@ def _bagging_resample(cfg, module, scenario_pool, xhat, serial=False, replacemen return local_boot_gaps, local_boot_optimals, local_boot_uppers, local_boot_counts -def bagging_bootstrap(cfg, module, xhat, replacement=True): +def bagging_bootstrap(cfg, module, xhat, replacement=True, executor=None): """ perform a bagging-based estimation of confidence intervals Args: @@ -708,40 +717,28 @@ def bagging_bootstrap(cfg, module, xhat, replacement=True): xhat (dict): keys are scenario tree node names (e.g. ROOT) and values are mpi-sppy nonant vectors (i.e. the specification of a candidate solution) replacement (bool): sample the subsample with replacement + executor (BatchExecutor or None): the rank grouping; None => K = 1 Returns: tuple with confidence interval if on MPI rank 0 """ - + if executor is None: + executor = BatchExecutor(1) rng = _pool_rng(cfg) scenario_pool = rng.choice(eligible_scenarios(cfg), size=cfg.sample_size, replace=False) # bootstrap from pool - local_boot_gaps, local_boot_optimals, local_boot_uppers, local_boot_counts = _bagging_resample(cfg, module, scenario_pool, xhat, serial=False, replacement=replacement) - comm.Barrier() - - # do analysis only on rank 0 - if my_rank == 0: - boot_gaps = np.empty(cfg.nB, dtype=np.float64) - boot_optimals = np.empty(cfg.nB, dtype=np.float64) - boot_uppers = np.empty(cfg.nB, dtype=np.float64) - boot_counts = np.empty(cfg.nB * cfg.sample_size, dtype=np.float64) - else: - boot_gaps = None - boot_optimals = None - boot_uppers = None - boot_counts = None - - # but everyone needs to send to the gather - lenlist = slice_lens(cfg.nB) - comm.Gatherv(sendbuf=local_boot_gaps, recvbuf=(boot_gaps, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_optimals, recvbuf=(boot_optimals, lenlist), root=0) - comm.Gatherv(sendbuf=local_boot_uppers, recvbuf=(boot_uppers, lenlist), root=0) + local_boot_gaps, local_boot_optimals, local_boot_uppers, local_boot_counts = _bagging_resample(cfg, module, scenario_pool, xhat, executor, replacement=replacement) + executor.comm.Barrier() - receive_len = [x * cfg.sample_size for x in lenlist] - comm.Gatherv(sendbuf=local_boot_counts, recvbuf=(boot_counts, receive_len), root=0) + # gather every group's batches (counts are sample_size floats per batch) to + # global rank 0 for analysis + boot_gaps = executor.gather(local_boot_gaps, cfg.nB) + boot_optimals = executor.gather(local_boot_optimals, cfg.nB) + boot_uppers = executor.gather(local_boot_uppers, cfg.nB) + boot_counts = executor.gather(local_boot_counts, cfg.nB, item_len=cfg.sample_size) - if my_rank == 0: + if executor.is_root: center_gap = np.mean(boot_gaps) center_optimal = np.mean(boot_optimals) center_upper = np.mean(boot_uppers) @@ -776,13 +773,19 @@ def bagging_bootstrap(cfg, module, xhat, replacement=True): return None, None, None, None, None, None -def compute_ci(cfg, module, xhat): +def compute_ci(cfg, module, xhat, executor=None): """ Dispatch to the requested bootstrap method and return its result. Args: cfg (Config): parameters (cfg.boot_method selects the method) module (Python module): contains the scenario creator function and helpers xhat (dict): a candidate solution in mpi-sppy nonant format + executor (BatchExecutor or None): the rank grouping for the batch solves + (design 9.4). None selects the default K = 1 grouping (each rank its + own batch worker, a direct EF per batch) used by the standalone + user_boot / simulate_boot drivers; generic_cylinders' do_boot passes + an executor with K = --boot-ranks-per-batch and a wheel-based batch + solver. Returns: (ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap); @@ -801,18 +804,20 @@ def compute_ci(cfg, module, xhat): f"boot_method={method} is a smoothed method; it is dispatched by " "smoothed_boot_sp.compute_smoothed_ci, not boot_sp.compute_ci " "(the drivers route smoothed methods automatically).") + if executor is None: + executor = BatchExecutor(1) if method == "Extended": - return extended_bootstrap(cfg, module, xhat) + return extended_bootstrap(cfg, module, xhat, executor=executor) elif method == "Bagging_with_replacement": - return bagging_bootstrap(cfg, module, xhat, replacement=True) + return bagging_bootstrap(cfg, module, xhat, replacement=True, executor=executor) elif method == "Bagging_without_replacement": - return bagging_bootstrap(cfg, module, xhat, replacement=False) + return bagging_bootstrap(cfg, module, xhat, replacement=False, executor=executor) elif method == "Classical_quantile": - return classical_bootstrap(cfg, module, xhat, quantile=True) + return classical_bootstrap(cfg, module, xhat, quantile=True, executor=executor) elif method == "Classical_gaussian": - return classical_bootstrap(cfg, module, xhat, quantile=False) + return classical_bootstrap(cfg, module, xhat, quantile=False, executor=executor) elif method == "Subsampling": - return subsampling(cfg, module, xhat) + return subsampling(cfg, module, xhat, executor=executor) else: raise ValueError(f"boot_method={method} is not supported.") diff --git a/mpisppy/generic/boot.py b/mpisppy/generic/boot.py index 25dc13d1e..d12317b18 100644 --- a/mpisppy/generic/boot.py +++ b/mpisppy/generic/boot.py @@ -61,6 +61,17 @@ def boot_requested(cfg): "boot_candidate_sample_size (M) and boot_xhat_input_file_name are " "mutually exclusive: give a positive M to find xhat from the " "candidate records, or an xhat file to read it (with M omitted or 0).") + + # The batch config file is required whenever a bootstrap run is requested, + # for K = 1 as well as K > 1 (design 9.5): one uniform mechanism, no + # simple-case shortcut. It configures how each resampled batch is solved. + if cfg.get("boot_batch_config_file") is None: + raise ValueError( + "a bootstrap run requires --boot-batch-config-file: a file of " + "generic_cylinders flags configuring how each resampled batch is " + "solved (its solver, and for K>1 its rho/spokes/convergence). For " + "K=1 it need only name a solver, e.g. a one-line file " + "'--solver-name gurobi'.") return True @@ -82,16 +93,8 @@ def _check_compatibility(cfg): raise ValueError( f"The bootstrap CI cannot be combined with --{opt.replace('_', '-')}.") - # K > 1 (a wheel per group of ranks on a sub-communicator) is a scheduled - # follow-on; the first integration ships K = 1 (a per-rank extensive form). - K = cfg.get("boot_ranks_per_batch", 1) - if K is not None and K != 1: - raise ValueError( - f"boot_ranks_per_batch={K}: only K=1 (a per-rank extensive form) " - "is supported so far.") - -def _estimator_cfg(module_basename, module, cfg, N, M, pool_names): +def _estimator_cfg(module_basename, module, cfg, batch_cfg, N, M, pool_names): """Build the Config the bootsp estimator expects from the boot_* options. The estimator (boot_sp) reads its own historical option names (max_count, @@ -100,6 +103,11 @@ def _estimator_cfg(module_basename, module, cfg, N, M, pool_names): addresses records by position 0..N-1; we set max_count = N so its resampling pool *is* the disjoint N-record block, and install a resolver mapping each position to its canonical scenario name. + + The batch solver name and options come from the parsed ``batch_cfg`` (the + --boot-batch-config-file), not the xhat-solve command line: they govern the + K=1 direct-EF solve and the xhat-evaluation solves, so both agree with the + K>1 wheel, which reads batch_cfg directly (design 9.4 / 9.5). """ import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils @@ -134,16 +142,14 @@ def _estimator_cfg(module_basename, module, cfg, N, M, pool_names): boot_cfg.alpha = cfg.get("boot_alpha") boot_cfg.seed_offset = cfg.get("boot_seed_offset") - # the batch ("boot") solver role, falling back to the generic solver_name - # (and its options), so the batch EF solves are independent of any xhat-EF - # solver - _, boot_solver_name, boot_solver_options = solver_specification(cfg, ["boot", ""]) - boot_cfg.solver_name = boot_solver_name + # the batch solver name and options, from the batch config file + _, batch_solver_name, batch_solver_options = solver_specification(batch_cfg, "") + boot_cfg.solver_name = batch_solver_name boot_cfg.add_to_config( "solver_options", description="options dict for the bootstrap batch solver", domain=None, default=None, argparse=False) - boot_cfg.solver_options = boot_solver_options + boot_cfg.solver_options = batch_solver_options # the positional resolver: estimator position p -> canonical pool name boot_cfg.add_to_config( @@ -245,9 +251,36 @@ def do_boot(module_fname, cfg, wheel=None): except OSError: pass - boot_cfg = _estimator_cfg(module_basename, module, cfg, N, M, pool_names) + # The batch config file (required) governs how each resampled batch is + # solved. For K=1 it feeds the estimator's direct-EF and xhat-eval solves; + # for K>1 it also builds the per-group wheel (design 9.4 / 9.5). + from mpisppy.generic import boot_batch + from mpisppy.confidence_intervals.bootsp.batch_executor import BatchExecutor - result = boot_sp.compute_ci(boot_cfg, module, xhat) + batch_cfg = boot_batch.parse_batch_config_file( + cfg.boot_batch_config_file, module) + + # Build the executor first so it validates K against the rank count before + # we build the (K>1) wheel solver, which validates the batch config's rho. + K = cfg.get("boot_ranks_per_batch", 1) + executor = BatchExecutor(K, comm=global_comm) + if executor.uses_cylinders: + executor.batch_optimal_solver = boot_batch.make_batch_optimal_solver( + batch_cfg, module) + # else K=1: each batch is a direct extensive form (no wheel solver) + + boot_cfg = _estimator_cfg(module_basename, module, cfg, batch_cfg, N, M, pool_names) + + # Refuse a maximization model before any batch is solved. boot_sp's + # solve_routine makes the same check, but only K=1 builds an extensive form + # for it to look at: a K>1 batch is solved by a wheel, so probe one scenario + # here instead. It is built exactly as the batches build theirs. + probe = module.scenario_creator(pool_names[0], **module.kw_creator(boot_cfg)) + boot_sp._require_minimization( + sputils.find_active_objective(probe).is_minimizing(), + f"module {module_basename}") + + result = boot_sp.compute_ci(boot_cfg, module, xhat, executor=executor) if global_rank == 0: ci_optimal, ci_upper, ci_gap, c_optimal, c_upper, c_gap = result diff --git a/mpisppy/generic/boot_batch.py b/mpisppy/generic/boot_batch.py new file mode 100644 index 000000000..811581a96 --- /dev/null +++ b/mpisppy/generic/boot_batch.py @@ -0,0 +1,143 @@ +############################################################################### +# 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. +############################################################################### +"""Solve one bootstrap batch with cylinders (design section 9.4, K > 1). + +For ``--boot-ranks-per-batch`` (K) greater than 1 a resampled batch is not a +direct extensive form but a full hub/spoke wheel run on the batch's group of K +ranks. How that wheel is configured -- solver, rho, which spokes, convergence, +the relative gap -- is a *different* problem from the xhat solve (a batch is a +resample of the data, with its own scenario count set independently of the M +candidate records), so it is supplied +separately as ``--boot-batch-config-file``: literally a file of +``generic_cylinders`` flags. This module parses that file into a Config and +turns it into a callable that the estimator invokes per batch, returning the +wheel's outer (decomposition) bound as the batch optimal ``L_b`` (design 9.4.1). + +The batch config must produce an *outer* bound (a Lagrangian/subgradient spoke, +or the subgradient hub); a bare PH hub with no bounding spoke leaves the outer +bound undefined and is rejected at solve time with a clear message. +""" + +import shlex + +from mpisppy.spin_the_wheel import WheelSpinner +from mpisppy.generic.parsing import register_generic_args +from mpisppy.generic.hub import build_hub_dict +from mpisppy.generic.spokes import build_spoke_list +from mpisppy.generic.extensions import configure_extensions +from mpisppy.generic import decomp +import mpisppy.utils.cfg_vanilla as vanilla +import mpisppy.utils.config as config +from mpisppy.confidence_intervals.bootsp.boot_sp import _scenario_creator_w_mapping + + +def _batch_scenario_denouement(rank, name, scenario): + """No-op denouement for the batch wheels (nothing to report per scenario).""" + pass + + +def parse_batch_config_file(path, module): + """Parse a batch config file (generic_cylinders flags) into a Config. + + The file is read as if its contents were the generic_cylinders command line + (``--solver-name gurobi --lagrangian --default-rho 1.0`` ...); ``#`` starts a + comment. It is parsed by the same Config machinery as the main run, so it is + exactly a batch generic_cylinders configuration. The framework -- not the + file -- supplies the batch's scenario set (its count and the positional + sample->record mapping), so the file must not set the scenario-formation + options; those are the ``--boot-*`` flags on the main command line. + + Args: + path (str): path to the batch config file. + module: the model module (its inparser_adder registers model options). + + Returns: + Config: the parsed batch configuration. + """ + cfg = config.Config() + register_generic_args(cfg, module) + parser = cfg.create_parser("boot batch config") + with open(path, "r") as f: + tokens = shlex.split(f.read(), comments=True) + args = parser.parse_args(tokens) + cfg.import_argparse(args) + cfg.checker() # same inconsistency checks as the main run + return cfg + + +def _outer_bound_over_group(wheel, groupcomm): + """The wheel's outer bound, made available on every rank of the group. + + ``WheelSpinner`` populates ``BestOuterBound`` only on the hub ranks + (strata_rank 0); the rest carry None. The group leader (group_rank 0), which + reports the batch result into the cross-group gather, is not necessarily a + hub rank, so share the bound across the whole group and let every rank return + the same value. + """ + candidates = [b for b in groupcomm.allgather(wheel.BestOuterBound) if b is not None] + if not candidates: + raise RuntimeError( + "boot batch solve produced no outer bound. The --boot-batch-config-file " + "must configure a wheel that yields an outer (decomposition) bound on " + "the batch optimal -- e.g. add a Lagrangian or subgradient spoke, or " + "use the subgradient hub. A bare PH hub does not provide one.") + bound = float(candidates[0]) + if bound in (float("inf"), float("-inf")): + raise RuntimeError( + "boot batch solve reported an infinite outer bound (no bounding " + "progress). Check the --boot-batch-config-file solver/spoke settings.") + return bound + + +def make_batch_optimal_solver(batch_cfg, module): + """Return a callable that solves one batch with cylinders for its outer bound. + + The returned callable has the signature expected by the estimator's + ``BatchExecutor.batch_optimal_solver``: + + solver(scenario_names, sample_mapping, groupcomm) -> outer_bound (float) + + It builds a fresh wheel over the batch's sample scenarios (mapped back to + their resampled records) using ``batch_cfg`` and runs it on ``groupcomm`` + (the K ranks assigned to this batch), returning the wheel's outer bound. + """ + rho_setter = decomp._get_rho_setter(module, batch_cfg) + ph_converger = decomp._get_converger(batch_cfg) + average_scenario_creator = getattr(module, "average_scenario_creator", None) + feasible_xhat_creator = vanilla._find_feasible_xhat_creator(module, batch_cfg) + + def solver(scenario_names, sample_mapping, groupcomm): + # a fresh kwargs dict each call: the mapping (and thus the batch) changes + scenario_creator_kwargs = module.kw_creator(batch_cfg) + scenario_creator_kwargs["module"] = module + scenario_creator_kwargs["mapping"] = sample_mapping + + # a batch has this many scenarios; set num_scens so a model that reads it + # for uniform probabilities matches the direct-EF (K=1) semantics. + batch_cfg.num_scens = len(scenario_names) + + beans = (batch_cfg, _scenario_creator_w_mapping, + _batch_scenario_denouement, scenario_names) + + hub_dict = build_hub_dict(batch_cfg, beans, scenario_creator_kwargs, + rho_setter, None, ph_converger) + configure_extensions(hub_dict, module, batch_cfg) + if batch_cfg.reduced_costs: + vanilla.add_reduced_costs_fixer(hub_dict, batch_cfg) + + list_of_spoke_dict = build_spoke_list( + batch_cfg, beans, scenario_creator_kwargs, rho_setter, None, + average_scenario_creator=average_scenario_creator, + feasible_xhat_creator=feasible_xhat_creator) + + wheel = WheelSpinner(hub_dict, list_of_spoke_dict) + wheel.run(comm_world=groupcomm) + return _outer_bound_over_group(wheel, groupcomm) + + return solver diff --git a/mpisppy/generic/parsing.py b/mpisppy/generic/parsing.py index a5361d0a8..bcff001b0 100644 --- a/mpisppy/generic/parsing.py +++ b/mpisppy/generic/parsing.py @@ -152,9 +152,20 @@ def add_decomp_args(cfg): default=None) -def parse_args(m): - """Parse CLI args given the model module m. Returns a Config object.""" - cfg = config.Config() +def register_generic_args(cfg, m): + """Register the full generic_cylinders option set on cfg (no parsing). + + Factored out of parse_args so a second consumer -- the bootstrap batch + config file (--boot-batch-config-file), which is literally a batch + generic_cylinders configuration -- can build an identical Config and then + read its options from a file rather than the command line. parse_args calls + this and then parses argv; the batch reader calls this and then imports the + file. Keeping a single registration point means the two never drift. + + Args: + cfg (Config): a fresh Config to populate. + m (module or class): the model, for its inparser_adder. + """ cfg.proper_bundle_config() cfg.pickle_scenarios_config() cfg.pre_pickle_args() @@ -193,6 +204,12 @@ def parse_args(m): from mpisppy.generic.admm import admm_args admm_args(cfg) + +def parse_args(m): + """Parse CLI args given the model module m. Returns a Config object.""" + cfg = config.Config() + register_generic_args(cfg, m) + cfg.parse_command_line(f"mpi-sppy for {cfg.module_name}") cfg.checker() # looks for inconsistencies diff --git a/mpisppy/tests/test_boot_generic.py b/mpisppy/tests/test_boot_generic.py index 9b0a6e663..8fb450461 100644 --- a/mpisppy/tests/test_boot_generic.py +++ b/mpisppy/tests/test_boot_generic.py @@ -51,6 +51,20 @@ MODULE_NAME = "schultz_data" +# A bootstrap run now requires --boot-batch-config-file (a file of +# generic_cylinders flags for the batch solves). For the K=1 tests it need only +# name the solver. One file per rank avoids a write race under mpiexec. +_batch_cfg_fd, BATCH_CFG_PATH = tempfile.mkstemp(prefix=f"bootbatch{my_rank}", suffix=".txt") +with os.fdopen(_batch_cfg_fd, "w") as _f: + _f.write(f"--solver-name {solver_name or 'gurobi'}\n") + +# For the K>1 (wheel-per-batch) test: a subgradient hub is a single cylinder +# (works at any rank count) that yields an outer (dual) bound on each batch. +_wheel_cfg_fd, WHEEL_BATCH_CFG_PATH = tempfile.mkstemp(prefix=f"bootwheel{my_rank}", suffix=".txt") +with os.fdopen(_wheel_cfg_fd, "w") as _f: + _f.write(f"--solver-name {solver_name or 'gurobi'}\n" + "--subgradient-hub\n--max-iterations 5\n--default-rho 1.0\n") + # A fixed, feasible candidate solution so the CI depends only on the # (deterministic) bootstrap draws over the dataset rows, not on which optimum a # given solver returns. @@ -94,10 +108,23 @@ def _make_cfg(method="Classical_quantile"): cfg.boot_nB = 20 cfg.boot_alpha = 0.1 cfg.boot_seed_offset = 100 + cfg.boot_batch_config_file = BATCH_CFG_PATH cfg.data_file = "schultz_data.csv" return cfg +def _make_small_cfg(K, batch_cfg_path): + """A cheap cfg for the wheel-vs-EF comparison: a small pool and few batches + so the per-batch wheel solves stay fast.""" + cfg = _make_cfg("Classical_quantile") + cfg.boot_sample_size = 8 + cfg.boot_subsample_size = 4 + cfg.boot_nB = 4 + cfg.boot_batch_config_file = batch_cfg_path + cfg.boot_ranks_per_batch = K + return cfg + + #***************************************************************************** class Test_boot_generic(unittest.TestCase): """Test do_boot and boot_requested through the generic_cylinders surface. @@ -145,22 +172,71 @@ def test_do_boot_xhat_from_solve(self): res = do_boot(MODULE_NAME, cfg, wheel=_FakeWheel()) self._assert_gap(res, locked_ci_gap_xhat_from_solve, locked_center_gap_xhat_from_solve) - def test_boot_solver_role_resolution(self): - # the batch ("boot") solver role resolves its own name and options, and - # falls back to the generic solver_name when --boot-solver-name is unset + def test_batch_solver_from_config_file(self): + # the estimator's solver name and options come from the batch config + # file, not the xhat-solve command line (PR-4 retired --boot-solver-*). + from mpisppy.generic import boot_batch import schultz_data pool = schultz_data.scenario_names_creator(3) + bf = tempfile.mkstemp(suffix=".txt")[1] + with open(bf, "w") as f: + f.write("--solver-name myboot_solver --solver-options mipgap=0.01\n") + try: + batch_cfg = boot_batch.parse_batch_config_file(bf, schultz_data) + cfg = _make_cfg("Classical_quantile") + boot_cfg = _estimator_cfg("schultz_data", schultz_data, cfg, batch_cfg, 3, 0, pool) + self.assertEqual(boot_cfg.solver_name, "myboot_solver") + self.assertIn("mipgap", boot_cfg.solver_options) + finally: + if os.path.exists(bf): + os.remove(bf) + + def test_boot_requested_requires_batch_config(self): + # a bootstrap run must name a batch config file cfg = _make_cfg("Classical_quantile") - cfg.boot_solver_name = "myboot_solver" - cfg.boot_solver_options = "mipgap=0.01" - boot_cfg = _estimator_cfg("schultz_data", schultz_data, cfg, 3, 0, pool) - self.assertEqual(boot_cfg.solver_name, "myboot_solver") - self.assertIn("mipgap", boot_cfg.solver_options) + cfg.boot_batch_config_file = None + with self.assertRaises(ValueError): + boot_requested(cfg) - cfg2 = _make_cfg("Classical_quantile") # no boot_solver_name - boot_cfg2 = _estimator_cfg("schultz_data", schultz_data, cfg2, 3, 0, pool) - self.assertEqual(boot_cfg2.solver_name, solver_name) + @unittest.skipIf(not solver_available, "no solver is available") + @unittest.skipIf(n_proc < 2, "the K>1 wheel path needs at least 2 ranks") + def test_do_boot_g1_wheel_matches_ef(self): + # The G=1 checkpoint (design 9.4): K = n_proc, so one group of all ranks + # solves each batch by a wheel in sequence. Compare that wheel path to + # the K=1 direct-EF path over the *same* pool (both deterministic, file + # xhat). The value at xhat (center_upper) is solver-exact and does not + # depend on K, so it must match exactly; the wheel's optimal is an outer + # (dual) bound, so it must sit at or below the EF optimum, making the + # reported gap conservative (>= the EF gap). + xf = tempfile.mkstemp(prefix=f"xhatw{my_rank}", suffix=".npy")[1] + ciutils.write_xhat(XHAT, path=xf) + try: + cfg_ef = _make_small_cfg(K=1, batch_cfg_path=BATCH_CFG_PATH) + cfg_ef.boot_xhat_input_file_name = xf + res_ef = do_boot(MODULE_NAME, cfg_ef) + + cfg_wheel = _make_small_cfg(K=n_proc, batch_cfg_path=WHEEL_BATCH_CFG_PATH) + cfg_wheel.boot_xhat_input_file_name = xf + res_wheel = do_boot(MODULE_NAME, cfg_wheel) + + if my_rank == 0: + co_ef, cu_ef, cg_ef = res_ef[3], res_ef[4], res_ef[5] + co_w, cu_w, cg_w = res_wheel[3], res_wheel[4], res_wheel[5] + ci_gap_w = list(res_wheel[2]) + tol = 1e-4 * (1 + abs(cu_ef)) + # value at xhat is K-invariant -> exact match + self.assertAlmostEqual(cu_w, cu_ef, delta=tol) + # wheel optimal is an outer bound -> at or below the EF optimum + self.assertLessEqual(co_w, co_ef + tol) + # so the reported gap over-states (is conservative) + self.assertGreaterEqual(cg_w, cg_ef - tol) + # and the CI is ordered and finite + self.assertLessEqual(ci_gap_w[0], ci_gap_w[1]) + self.assertTrue(all(abs(v) < float("inf") for v in ci_gap_w)) + finally: + if os.path.exists(xf): + os.remove(xf) def test_boot_requested_none(self): cfg = _make_cfg() @@ -196,14 +272,23 @@ def test_do_boot_disjoint_overflow_raises(self): with self.assertRaises(ValueError): do_boot(MODULE_NAME, cfg, wheel=_FakeWheel()) - def test_do_boot_ranks_per_batch_gt_1_raises(self): + def test_do_boot_bad_ranks_per_batch_raises(self): + # K must divide the number of MPI ranks; serially (1 rank) any K > 1 is + # invalid (K > R), so the executor rejects it. cfg = _make_cfg("Classical_quantile") cfg.boot_candidate_sample_size = 5 cfg.num_scens = 5 - cfg.boot_ranks_per_batch = 2 + cfg.boot_ranks_per_batch = n_proc + 1 with self.assertRaises(ValueError): do_boot(MODULE_NAME, cfg, wheel=_FakeWheel()) +def tearDownModule(): + # remove the per-rank batch config files created at import time + for path in (BATCH_CFG_PATH, WHEEL_BATCH_CFG_PATH): + if os.path.exists(path): + os.remove(path) + + if __name__ == '__main__': unittest.main() diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index d71188fb4..7d2214af8 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -1851,24 +1851,25 @@ def boot_args(self): default=None, ) self.add_to_config( - "boot_solver_name", - description="Solver for the bootstrap batch solves; falls back to the" - " generic solver_name when unset (default None)", - domain=str, - default=None, - ) - self.add_to_config( - "boot_solver_options", - description="Options string for the bootstrap batch solver, e.g." - " 'mipgap=0.01' (default None)", + "boot_batch_config_file", + description="Required for a bootstrap run: a file of generic_cylinders" + " flags (e.g. '--solver-name gurobi --lagrangian --default-rho 1.0')" + " that configures how each resampled batch is solved. It is a separate" + " problem from the xhat solve (a batch is a resample of the data, sized" + " independently of the M candidate records), so its solver, rho, spokes," + " convergence and relative gap are set here, not inherited. For K=1 it" + " need only name a" + " solver (a direct EF); for K>1 it is the group's full cylinder" + " configuration (default None)", domain=str, default=None, ) self.add_to_config( "boot_ranks_per_batch", - description="K: number of MPI ranks that cooperate on one batch solve." - " Only K=1 (a per-rank extensive form) is supported so far" - " (default 1)", + description="K: number of MPI ranks that cooperate on one batch solve" + " (design 9.4). K=1 solves each batch as a per-rank extensive form;" + " K>1 runs a wheel per group of K ranks using boot_batch_config_file." + " K must divide the number of MPI ranks (default 1)", domain=int, default=1, )