diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 8bd672f7e..cf0e7374d 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -829,6 +829,25 @@ jobs: conda install mpi4py "numpy" setuptools pip install pyomo pandas xpress cplex scipy sympy dill packaging coverage + # the bootstrap epi-spline distributions fit with a nonlinear solver, so + # without ipopt those tests skip; this is how Pyomo's own CI gets ipopt + - name: Install Ipopt + run: | + # the ipopt binary links against these; Pyomo's CI installs the same + sudo apt-get update -q + sudo apt-get install -y libopenblas-dev gfortran liblapack-dev + IPOPT_DIR=$HOME/ipopt + mkdir -p "$IPOPT_DIR" + echo "$IPOPT_DIR" >> $GITHUB_PATH + echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$IPOPT_DIR" >> $GITHUB_ENV + URL=https://github.com/IDAES/idaes-ext + VER=$(curl -sL -H 'Accept: application/json' $URL/releases/latest \ + | sed -e 's/.*"tag_name":"\([^"]*\)".*/\1/') + if test -z "$VER"; then echo "FAILED identifying the ipopt release"; exit 1; fi + curl -fL $URL/releases/download/$VER/idaes-solvers-ubuntu2204-x86_64.tar.gz \ + | tar -xz -C "$IPOPT_DIR" + "$IPOPT_DIR"/ipopt -v + - name: setup the program run: | pip install -e . @@ -863,6 +882,7 @@ jobs: cd mpisppy/tests coverage run $COV_ARGS test_boot_sp.py coverage run $COV_ARGS test_boot_sp_simulate.py + coverage run $COV_ARGS test_boot_sp_smoothed.py - name: run bootstrap CI tests (mpiexec -np 2) timeout-minutes: 10 @@ -870,6 +890,7 @@ jobs: cd mpisppy/tests mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp.py mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp_simulate.py + mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp_smoothed.py - name: Upload coverage data if: always() diff --git a/.ruff.toml b/.ruff.toml index 06bbb4bac..9721fd735 100644 --- a/.ruff.toml +++ b/.ruff.toml @@ -6,3 +6,18 @@ extend-exclude = [ "./examples/hydro/hydro.py", "./examples/sizes/models/ExpressionModel.py", ] + +[lint.per-file-ignores] +# The bootsp/statdist/ subpackage is a faithful port of the (legacy) statdist +# distribution library; relax the style rules it predates rather than rewrite +# its numerics. splines.py builds a Pyomo model with `from pyomo.environ import *` +# (like the excluded Pyomo model files above), hence F403/F405. +"mpisppy/confidence_intervals/bootsp/statdist/*" = [ + "E711", # comparison to None + "E722", # bare except + "E731", # lambda assignment + "E741", # ambiguous variable name (math notation, e.g. l) + "F403", # star import (pyomo.environ) + "F405", # name may be from star import + "F821", # undefined name (in retained multivariate base class) +] diff --git a/doc/designs/bootsp_merge_design.md b/doc/designs/bootsp_merge_design.md index 0c7cfa0a1..00956c27d 100644 --- a/doc/designs/bootsp_merge_design.md +++ b/doc/designs/bootsp_merge_design.md @@ -1,12 +1,16 @@ # Bootstrap/bagging for data-based stochastic programming in mpi-sppy — design -**Status:** design captured and decisions ratified 2026-07-02; PR-1 -(empirical core + schultz, incl. a data-file example) implemented and open -upstream as draft [Pyomo/mpi-sppy#783](https://github.com/Pyomo/mpi-sppy/pull/783); -extended 2026-07-03 to state the end goal -(`generic_cylinders` integration) and a stacked, multi-PR roadmap (§6, §9). +**Status:** design captured and decisions ratified 2026-07-02; extended +2026-07-03 to state the end goal (`generic_cylinders` integration) and a +stacked, multi-PR roadmap (§6, §9). PR-1 (empirical core + schultz, incl. a +data-file example) merged 2026-07-24 as +[Pyomo/mpi-sppy#783](https://github.com/Pyomo/mpi-sppy/pull/783), followed by +a `test_boot_sp.py` `np=2` fix merged 2026-07-28 as +[#820](https://github.com/Pyomo/mpi-sppy/pull/820). PR-2 (statdist + smoothed +methods) is this branch, open upstream as +[#818](https://github.com/Pyomo/mpi-sppy/pull/818); PR-3 not yet started. **Author:** dlw (captured with Claude Code assistance) -**Last updated:** 2026-07-03 +**Last updated:** 2026-07-28 **Ultimate goal.** The end state this design builds toward is *bootstrap and bagging confidence intervals, computed from a given dataset, available @@ -286,6 +290,33 @@ Behavior-preserving unless noted. `boot_method` dies with a raw `TypeError` instead of a friendly message. The port accepts real json booleans (keeping the strings for boot-sp files) and reports a missing `boot_method` clearly. +13. **Independent smoothed-bootstrap batches (behavior change, from + broken to working).** `smoothed_resample_helper` advanced the record + index by one per batch while taking a block of `subsample_size` + consecutive records, so consecutive batches shared all but one of + their draws — a sliding window, not independent resamples. Chen & + Woodruff (2024, Algorithm 3) draws a fresh set of `N` points from the + fitted distribution for each of the `B` batches. The port strides by + the batch size, as `smoothed_bagging` already did, so the blocks are + pairwise disjoint. The estimated spread was badly understated before: + on the `cvar` example (`N = 20`, `nB = 10`) the interval widens about + fivefold. `simulate_boot` spaces its coverage replications to match: + the ported spacing was a hard-coded `nB * 100` unrelated to how many + record numbers a replication actually consumes, and is now exactly + that footprint. (The empirical harness needs no such spacing at all, + because item 11 gave it independent numpy streams; the smoothed path + addresses its draws by record number, since that is what the model + seeds each draw with, so its replications are separated by giving + each one a disjoint block of record numbers.) +14. **Smoothed center from the fitted distribution (behavior change, + from broken to working).** `smoothed_bootstrap` called + `center_smoothed` while `use_fitted` was still `False`, so the center + was the purely empirical gap. Algorithm 3 takes the center from + Algorithm 2 run on the *fitted* distribution, and that smoothed center + is the paper's leading conclusion, so `use_fitted` is now set before + the center is estimated (`smoothed_bagging` already did this). The + center block of the index space is reserved ahead of the batch blocks, + since both now sample the same fitted distribution. --- diff --git a/doc/src/boot_sp.rst b/doc/src/boot_sp.rst index 0f1230af8..6e3bc38a5 100644 --- a/doc/src/boot_sp.rst +++ b/doc/src/boot_sp.rst @@ -11,12 +11,14 @@ mpi-sppy, no distribution of the uncertain data is assumed: the estimators work directly from sampled data. The methods and software are described in [ChenWoodruff2023]_ and [ChenWoodruff2024]_. -.. note:: - - This is the empirical (numpy-only) part of the package: the classical, - extended, subsampling, and bagging methods. The *smoothed* methods, which - depend on a distribution-fitting library, are merged separately; asking for - a ``Smoothed_*`` method raises an informative error until then. +The package has two families of estimators. The *empirical* methods +(classical, extended, subsampling, and bagging) resample the observed data +directly and need only numpy. The *smoothed* methods fit a univariate +distribution to the sampled data (using the bundled ``statdist`` library) and +resample from the fitted distribution; they need `scipy +`_, which mpi-sppy treats as an optional dependency and +imports lazily. If scipy is not installed, the empirical methods still work +and a smoothed method fails with an informative import error. Modes ----- @@ -60,6 +62,10 @@ plus a few helpers used by the bootstrap code: for this fixed name first and falls back to the legacy ``xhat_generator_``. If a precomputed ``xhat`` file is given (``--xhat-fname``) the generator is not called. +* ``data_sampler(record_num, cfg)`` — return the data for one record (a scalar, + or a dict keyed by variable name for multivariate data). This is used by the + *smoothed* methods to build the sample that a distribution is fitted to; the + empirical methods do not need it. Methods ------- @@ -84,6 +90,21 @@ The ``--boot-method`` (json ``boot_method``) option selects the estimator: - Bagging with replacement [lam2018]_ * - ``Bagging_without_replacement`` - Bagging without replacement [lam2018]_ + * - ``Smoothed_boot_epi`` + - Smoothed bootstrap, epi-spline fit, Gaussian interval [ChenWoodruff2024]_ + * - ``Smoothed_boot_kernel`` + - Smoothed bootstrap, kernel-density fit, Gaussian interval [ChenWoodruff2024]_ + * - ``Smoothed_boot_epi_quantile`` + - Smoothed bootstrap, epi-spline fit, quantile interval [ChenWoodruff2024]_ + * - ``Smoothed_boot_kernel_quantile`` + - Smoothed bootstrap, kernel-density fit, quantile interval [ChenWoodruff2024]_ + * - ``Smoothed_bagging`` + - Smoothed bagging, kernel-density fit [ChenWoodruff2024]_ + +The ``Smoothed_*`` tokens are the smoothed methods; the others are empirical. +The epi-spline fit builds a small Pyomo nonlinear program, so those two methods +additionally need a nonlinear solver (e.g. ``ipopt``); the kernel methods do +not. Arguments --------- @@ -112,6 +133,14 @@ command line (with dashes). The main options are: * ``coverage_replications`` (simulation only) — number of coverage replications. * ``boot_method`` / ``--boot-method`` — one of the tokens above. +The smoothed methods use two additional options (ignored, and not required in +the json, for the empirical methods): + +* ``smoothed_center_sample_size`` / ``--smoothed-center-sample-size`` — number + of points drawn from the fitted distribution to estimate the gap center. +* ``smoothed_B_I`` / ``--smoothed-B-I`` — number of outer replications for + smoothed bagging. + There may also be model-specific options added by ``inparser_adder``. Batch parallelism @@ -190,6 +219,46 @@ from the confidence-interval sampling, so ``sample_size`` plus makes it reproducible); replace it with your own two-column dataset, or point ``--data-file`` at another file, to run the bootstrap on your own data. +Smoothed methods and statdist +----------------------------- + +The smoothed methods (the ``Smoothed_*`` tokens) fit a univariate distribution +to the sampled data and then resample from the *fitted* distribution rather +than from the data directly. The distribution fitting is provided by the +bundled ``statdist`` library +(``mpisppy.confidence_intervals.bootsp.statdist``), a trimmed port of the +univariate distributions from the statdist package; ``statdist`` uses scipy, +which is imported lazily so that the empirical methods remain scipy-free. + +To use a smoothed method the model module must supply ``data_sampler`` (see +above): the smoothed estimator calls it for each sampled record to assemble the +data that ``statdist`` fits. The kernel-density methods +(``Smoothed_boot_kernel``, ``Smoothed_boot_kernel_quantile``, +``Smoothed_bagging``) fit with a Gaussian kernel and need only scipy; the +epi-spline methods (``Smoothed_boot_epi``, ``Smoothed_boot_epi_quantile``) fit +by solving a small Pyomo nonlinear program and additionally need a nonlinear +solver such as ``ipopt``. + +Three examples that need statdist ship in ``examples/bootsp``: + +* ``farmer`` — the scalable farmer, with crop yields perturbed by a fitted + (or, empirically, a uniform) distribution; +* ``cvar`` — a CVaR example (Lam & Qian) with standard-normal data; +* ``multi_knapsack`` — a multi-product knapsack (Vaagen & Wallace) whose + deterministic data is read from a json file (``--deterministic-data-json``). + +Each has an empirical json/bash and a ``smoothed_*.json``; for instance, from +``examples/bootsp/cvar``: + +.. code-block:: bash + + $ python -m mpisppy.confidence_intervals.bootsp.user_boot cvar \ + --max-count 3000 --candidate-sample-size 10 --sample-size 75 \ + --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 0 \ + --solver-name cplex_direct --boot-method Bagging_with_replacement + + $ python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_cvar.json + References ---------- diff --git a/examples/bootsp/cvar/cvar.bash b/examples/bootsp/cvar/cvar.bash new file mode 100644 index 000000000..4141224ba --- /dev/null +++ b/examples/bootsp/cvar/cvar.bash @@ -0,0 +1,19 @@ +#!/bin/bash +# Run the CVaR bootstrap example (needs the statdist library). +# Pass a solver name as the first argument (default: cplex_direct). + +SOLVER=${1:-cplex_direct} +BOOT="python -m mpisppy.confidence_intervals.bootsp.user_boot" +COMMON="--max-count 3000 --candidate-sample-size 10 --sample-size 75 \ + --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 0 \ + --solver-name ${SOLVER}" + +echo "Serial, compute xhat within user_boot (empirical Bagging_with_replacement)" +echo +time ${BOOT} cvar ${COMMON} --boot-method Bagging_with_replacement +echo +echo "========================" +echo +echo "Smoothed coverage simulation from a json file (Smoothed_bagging)" +echo +time python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_cvar.json diff --git a/examples/bootsp/cvar/cvar.json b/examples/bootsp/cvar/cvar.json new file mode 100644 index 000000000..683450e94 --- /dev/null +++ b/examples/bootsp/cvar/cvar.json @@ -0,0 +1,16 @@ +{ + "module_name": "cvar", + "max_count": 3000, + "candidate_sample_size": 10, + "sample_size": 75, + "subsample_size": 10, + "nB": 20, + "alpha": 0.1, + "seed_offset": 0, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "boot_method": "Bagging_with_replacement", + "trace_fname": "None", + "coverage_replications": 5 +} diff --git a/examples/bootsp/cvar/cvar.py b/examples/bootsp/cvar/cvar.py new file mode 100644 index 000000000..6b7c501b6 --- /dev/null +++ b/examples/bootsp/cvar/cvar.py @@ -0,0 +1,174 @@ +############################################################################### +# 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. +############################################################################### +# A CVaR example (as in the Lam & Qian paper) for the bootstrap +# confidence-interval code. The scenario data are draws from a standard normal +# (empirical path) or from a statdist distribution fitted to the sample data +# (smoothed path), so importing this example needs the statdist library. + +import pyomo.environ as pyo +import mpisppy.scenario_tree as scenario_tree +import mpisppy.utils.sputils as sputils +import numpy as np +# importing Sampler pulls in the statdist package; the smoothed path fits a +# distribution (upstream) and samples it here, so cvar itself needs only Sampler +from mpisppy.confidence_intervals.bootsp.statdist.sampler import Sampler + +# Use this random stream: +sstream = np.random.RandomState(1) + + +def make_model(xi, num_scens, alpha=0.1): + + # Create the concrete model object + model = pyo.ConcreteModel("Lam_CVaR") + + model.nu = pyo.Var(within=pyo.NonNegativeReals) # second stage (xi - x)+ in L&Q + model.eta = pyo.Var(within=pyo.Reals) # first stage (x in Lam and Qian) + + model.Obj1 = pyo.Expression(expr=model.eta + (model.nu/alpha)) + + model.obj = pyo.Objective(expr=model.Obj1) + + def excess_rule(m): + return m.nu >= xi - m.eta + model.excess_constraint = pyo.Constraint(rule=excess_rule) + + # Create the list of nodes associated with the scenario (for two stage, + # there is only one node associated with the scenario--leaf nodes are + # ignored). + model._mpisppy_node_list = [ + scenario_tree.ScenarioNode( + name="ROOT", + cond_prob=1.0, + stage=1, + cost_expression=model.Obj1, + nonant_list=[model.eta], + scen_model=model, + ) + ] + + # Add the probability of the scenario + if num_scens is not None: + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" + return model + + +def data_sampler(record_num, cfg): + # return a single point from a sample + # Note: we are syncronizing using the seed + sstream.seed(record_num + cfg.seed_offset) + xi = sstream.normal(0, 1) + return xi + + +def scenario_creator(scenario_name, cfg): + """ Create a CVaR scenario. + + Args: + scenario_name (str): + Name of the scenario to construct. + cfg (Config): the control parameters + """ + # scenario_name has the form e.g. scen12, foobar7 + # The digits are scraped off the right of scenario_name using regex. + scennum = sputils.extract_num(scenario_name) + sstream.seed(scennum + cfg.seed_offset) # allows for resampling easily + + if getattr(cfg, "use_fitted", False): + # sampler works with a list + sampler = Sampler([cfg.fitted_distribution], sstream) + xi = sampler.sample_one()[0] + else: + xi = sstream.normal(0, 1) + num_scens = cfg.get('num_scens', None) + return make_model(xi, num_scens, alpha=0.1) + + +#========= +def scenario_names_creator(num_scens, start=None): + # (only for Amalgamator): return the full list of num_scens scenario names + # if start!=None, the list starts with the 'start' labeled scenario + if (start is None): + start = 0 + return [f"scen{i}" for i in range(start, start+num_scens)] + + +#========= +def inparser_adder(cfg): + # add options unique to the model + pass + + +#========= +def kw_creator(cfg): + # linked to the scenario_creator and inparser_adder + kwargs = {"cfg": cfg} + return kwargs + + +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + """ Create a scenario within a sample tree. Mainly for multi-stage and simple for two-stage. + (this function supports zhat and confidence interval code) + Args: + sname (string): scenario name to be created + stage (int >=1 ): for stages > 1, fix data based on sname in earlier stages + sample_branching_factors (list of ints): branching factors for the sample tree + seed (int): To allow random sampling (for some problems, it might be scenario offset) + given_scenario (Pyomo concrete model): if not None, use this to get data for ealier stages + scenario_creator_kwargs (dict): keyword args for the standard scenario creator funcion + Returns: + scenario (Pyomo concrete model): A scenario for sname with data in stages < stage determined + by the arguments + """ + # Since this is a two-stage problem, we don't have to do much. + sca = scenario_creator_kwargs.copy() + sca["seed_offset"] = seed + sca["num_scens"] = sample_branching_factors[0] # two-stage problem + return scenario_creator(sname, **sca) + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + pass + + +#============================ +def xhat_generator(scenario_names, solver_name=None, solver_options=None, cfg=None): + """ Solve the extensive form over the given scenarios and return xhat. + + This is the fixed-name generator the bootstrap code calls when no xhat file + is supplied (see boot_utils.compute_xhat). It builds the EF directly from + this module's scenario_creator so the example is self-contained. + + Args: + scenario_names (list of str): scenarios to build the EF from + solver_name (str): solver to use + solver_options (dict, optional): options passed to the solver + cfg (Config): control parameters + Returns: + xhat (dict): the first-stage nonants keyed by tree node (e.g. ROOT) + """ + ef = sputils.create_EF( + scenario_names, + scenario_creator, + scenario_creator_kwargs={"cfg": cfg}, + ) + solver = pyo.SolverFactory(solver_name) + if solver_options is not None: + for k, v in solver_options.items(): + solver.options[k] = v + if 'persistent' in solver_name: + solver.set_instance(ef, symbolic_solver_labels=True) + solver.solve(tee=False) + else: + solver.solve(ef, tee=False, symbolic_solver_labels=True) + return sputils.nonant_cache_from_ef(ef) diff --git a/examples/bootsp/cvar/smoothed_cvar.json b/examples/bootsp/cvar/smoothed_cvar.json new file mode 100644 index 000000000..50dc9b5b5 --- /dev/null +++ b/examples/bootsp/cvar/smoothed_cvar.json @@ -0,0 +1,18 @@ +{ + "module_name": "cvar", + "max_count": 3000, + "candidate_sample_size": 10, + "sample_size": 20, + "subsample_size": 3, + "smoothed_B_I": 5, + "smoothed_center_sample_size": 40, + "nB": 5, + "alpha": 0.05, + "seed_offset": 11, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "trace_fname": "None", + "boot_method": "Smoothed_bagging", + "coverage_replications": 5 +} diff --git a/examples/bootsp/farmer/farmer.bash b/examples/bootsp/farmer/farmer.bash new file mode 100644 index 000000000..1129d9ad8 --- /dev/null +++ b/examples/bootsp/farmer/farmer.bash @@ -0,0 +1,27 @@ +#!/bin/bash +# Run the farmer bootstrap example (needs the statdist library). +# Pass a solver name as the first argument (default: cplex_direct). + +SOLVER=${1:-cplex_direct} +BOOT="python -m mpisppy.confidence_intervals.bootsp.user_boot" +COMMON="--max-count 300 --candidate-sample-size 5 --sample-size 50 \ + --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 100 \ + --crops-multiplier 1 --yield-cv 0.1 --solver-name ${SOLVER}" + +echo "Serial, compute xhat within user_boot (empirical Bagging_with_replacement)" +echo +time ${BOOT} farmer ${COMMON} --boot-method Bagging_with_replacement +echo +echo "========================" +echo +echo "Parallel batches with mpiexec -np 2 (empirical Bagging_with_replacement)" +echo +time mpiexec -np 2 python -m mpi4py \ + -m mpisppy.confidence_intervals.bootsp.user_boot \ + farmer ${COMMON} --boot-method Bagging_with_replacement +echo +echo "========================" +echo +echo "Smoothed coverage simulation from a json file (Smoothed_bagging)" +echo +time python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_farmer.json diff --git a/examples/bootsp/farmer/farmer.json b/examples/bootsp/farmer/farmer.json new file mode 100644 index 000000000..2730ca164 --- /dev/null +++ b/examples/bootsp/farmer/farmer.json @@ -0,0 +1,19 @@ +{ + "module_name": "farmer", + "max_count": 300, + "candidate_sample_size": 5, + "sample_size": 50, + "subsample_size": 10, + "nB": 20, + "alpha": 0.1, + "seed_offset": 100, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "boot_method": "Bagging_with_replacement", + "trace_fname": "None", + "coverage_replications": 5, + "crops_multiplier": 1, + "farmer_with_integers": "False", + "yield_cv": "0.1" +} diff --git a/examples/bootsp/farmer/farmer.py b/examples/bootsp/farmer/farmer.py new file mode 100644 index 000000000..4f8a4842c --- /dev/null +++ b/examples/bootsp/farmer/farmer.py @@ -0,0 +1,399 @@ +############################################################################### +# 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. +############################################################################### +# A scalable "farmer" example for the bootstrap confidence-interval code. The +# crop yields fluctuate around a baseline according to a statdist univariate +# distribution (unif(0,1) by default, or a fitted distribution on the smoothed +# path), so importing this example needs the statdist library. + +import pyomo.environ as pyo +import numpy as np +import mpisppy.scenario_tree as scenario_tree +import mpisppy.utils.sputils as sputils +import mpisppy.confidence_intervals.bootsp.statdist as statdist +from mpisppy.confidence_intervals.bootsp.statdist.sampler import Sampler + +# Use this random stream: +farmerstream = np.random.RandomState() + + +def _get_distr_dict(cfg): + + def _get_b(c, cv): + # c is approximately the lower bound of crop yield, cv is approx coefficient of variation + # if no specified yield_cv, use the original scalable farmer unif(0,1) + if cv is None: + return 1 + else: + return c*cv/(1/np.sqrt(12) - cv/2) + + if not getattr(cfg, "use_fitted", False): + uunif = statdist.distribution_factory('univariate-unif') + distr_dict = {} + for i in range(cfg.crops_multiplier): + distr_dict[f"WHEAT{i}"] = uunif(0, _get_b(2.5, cfg.yield_cv)) + distr_dict[f"CORN{i}"] = uunif(0, _get_b(3, cfg.yield_cv)) + distr_dict[f"SUGAR_BEETS{i}"] = uunif(0, _get_b(20, cfg.yield_cv)) + else: + distr_dict = cfg.fitted_distribution + return distr_dict + + +def scenario_creator( + scenario_name, cfg, sense=pyo.minimize, seed_offset=None +): + """ Create a scenario for the (scalable) farmer example. + Args: + scenario_name (str): + Name of the scenario to construct. + cfg (Config): + control parameters + sense (int, optional): + Model sense (minimization or maximization). Must be either + pyo.minimize or pyo.maximize. Default is pyo.minimize. + seed_offset (int): used by confidence interval code + Note: + if cfg.yield_cv is None, give the behavior of the original scalable farmer + """ + # scenario_name has the form e.g. scen12, foobar7 + # The digits are scraped off the right of scenario_name using regex then + # converted mod 3 into one of the below avg./avg./above avg. scenarios + scennum = sputils.extract_num(scenario_name) + basenames = ['BelowAverageScenario', 'AverageScenario', 'AboveAverageScenario'] + basenum = scennum % 3 + groupnum = scennum // 3 + scenname = basenames[basenum]+str(groupnum) + + # The RNG is seeded with the scenario number so that it is + # reproducible when used with multiple threads. + # NOTE: if you want to do replicates, you will need to pass a seed + # as a kwarg to scenario_creator then use seed+scennum as the seed argument. + seed_offset = cfg.get("seed_offset", 0) if seed_offset is None else seed_offset + + farmerstream.seed(scennum+seed_offset) + + use_integer = cfg.get('use_integer', False) + crops_multiplier = cfg.get('crops_multiplier', 1) + num_scens = cfg.get('num_scens', None) + + # Check for minimization vs. maximization + if sense not in [pyo.minimize, pyo.maximize]: + raise ValueError("Model sense Not recognized") + + distr_dict = _get_distr_dict(cfg) + + # Create the concrete model object + model = pysp_instance_creation_callback( + scenname, + use_integer=use_integer, + sense=sense, + crops_multiplier=crops_multiplier, + distr_dict=distr_dict, + num_scens=num_scens + ) + + return model + + +def data_sampler(record_num, cfg): + # return the fluctuation data around the baseline from a sample + # Note: we are syncronizing using the seed + # yield as in "crop yield" + + distr_dict = _get_distr_dict(cfg) + farmerstream.seed(record_num+cfg.seed_offset) + groupnum = record_num // 3 + + sampler_dict = {} + for i in range(cfg.crops_multiplier): + if groupnum != 0: + sampler_dict[f"WHEAT{i}"] = Sampler([distr_dict[f"WHEAT{i}"]], farmerstream) + sampler_dict[f"CORN{i}"] = Sampler([distr_dict[f"CORN{i}"]], farmerstream) + sampler_dict[f"SUGAR_BEETS{i}"] = Sampler([distr_dict[f"SUGAR_BEETS{i}"]], farmerstream) + + data = {} + for i in range(cfg.crops_multiplier): + if groupnum != 0: + data[f"WHEAT{i}"] = sampler_dict[f"WHEAT{i}"].sample_one()[0] + data[f"CORN{i}"] = sampler_dict[f"CORN{i}"].sample_one()[0] + data[f"SUGAR_BEETS{i}"] = sampler_dict[f"SUGAR_BEETS{i}"].sample_one()[0] + else: + data[f"WHEAT{i}"] = 0 + data[f"CORN{i}"] = 0 + data[f"SUGAR_BEETS{i}"] = 0 + return data + + +def pysp_instance_creation_callback( + scenario_name, use_integer=False, sense=pyo.minimize, crops_multiplier=1, distr_dict=None, num_scens=None +): + # long function to create the entire model + # scenario_name is a string (e.g. AboveAverageScenario0) + # + # Returns a concrete model for the specified scenario + + # scenarios come in groups of three + scengroupnum = sputils.extract_num(scenario_name) + scenario_base_name = scenario_name.rstrip("0123456789") + + model = pyo.ConcreteModel() + + def crops_init(m): + retval = [] + for i in range(crops_multiplier): + retval.append("WHEAT"+str(i)) + retval.append("CORN"+str(i)) + retval.append("SUGAR_BEETS"+str(i)) + return retval + + model.CROPS = pyo.Set(initialize=crops_init) + + # + # Parameters + # + + model.TOTAL_ACREAGE = 500.0 * crops_multiplier + + def _scale_up_data(indict): + outdict = {} + for i in range(crops_multiplier): + for crop in ['WHEAT', 'CORN', 'SUGAR_BEETS']: + outdict[crop+str(i)] = indict[crop] + return outdict + + model.PriceQuota = _scale_up_data( + {'WHEAT': 100000.0, 'CORN': 100000.0, 'SUGAR_BEETS': 6000.0}) + + model.SubQuotaSellingPrice = _scale_up_data( + {'WHEAT': 170.0, 'CORN': 150.0, 'SUGAR_BEETS': 36.0}) + + model.SuperQuotaSellingPrice = _scale_up_data( + {'WHEAT': 0.0, 'CORN': 0.0, 'SUGAR_BEETS': 10.0}) + + model.CattleFeedRequirement = _scale_up_data( + {'WHEAT': 200.0, 'CORN': 240.0, 'SUGAR_BEETS': 0.0}) + + model.PurchasePrice = _scale_up_data( + {'WHEAT': 238.0, 'CORN': 210.0, 'SUGAR_BEETS': 100000.0}) + + model.PlantingCostPerAcre = _scale_up_data( + {'WHEAT': 150.0, 'CORN': 230.0, 'SUGAR_BEETS': 260.0}) + + # + # Stochastic Data + # + Yield = {} + Yield['BelowAverageScenario'] = \ + {'WHEAT': 2.0, 'CORN': 2.4, 'SUGAR_BEETS': 16.0} + Yield['AverageScenario'] = \ + {'WHEAT': 2.5, 'CORN': 3.0, 'SUGAR_BEETS': 20.0} + Yield['AboveAverageScenario'] = \ + {'WHEAT': 3.0, 'CORN': 3.6, 'SUGAR_BEETS': 24.0} + + def Yield_init(m, cropname): + # yield as in "crop yield" + sampler = Sampler([distr_dict[cropname]], farmerstream) + crop_base_name = cropname.rstrip("0123456789") + if scengroupnum != 0: + pertubation = sampler.sample_one()[0] + return Yield[scenario_base_name][crop_base_name] + pertubation + else: + return Yield[scenario_base_name][crop_base_name] + + model.Yield = pyo.Param(model.CROPS, + within=pyo.NonNegativeReals, + initialize=Yield_init, + mutable=True) + + # + # Variables + # + + if (use_integer): + model.DevotedAcreage = pyo.Var(model.CROPS, + within=pyo.NonNegativeIntegers, + bounds=(0.0, model.TOTAL_ACREAGE)) + else: + model.DevotedAcreage = pyo.Var(model.CROPS, + bounds=(0.0, model.TOTAL_ACREAGE)) + + model.QuantitySubQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantitySuperQuotaSold = pyo.Var(model.CROPS, bounds=(0.0, None)) + model.QuantityPurchased = pyo.Var(model.CROPS, bounds=(0.0, None)) + + # + # Constraints + # + + def ConstrainTotalAcreage_rule(model): + return pyo.sum_product(model.DevotedAcreage) <= model.TOTAL_ACREAGE + + model.ConstrainTotalAcreage = pyo.Constraint(rule=ConstrainTotalAcreage_rule) + + def EnforceCattleFeedRequirement_rule(model, i): + return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i] + + model.EnforceCattleFeedRequirement = pyo.Constraint(model.CROPS, rule=EnforceCattleFeedRequirement_rule) + + def LimitAmountSold_rule(model, i): + return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0 + + model.LimitAmountSold = pyo.Constraint(model.CROPS, rule=LimitAmountSold_rule) + + def EnforceQuotas_rule(model, i): + return (0.0, model.QuantitySubQuotaSold[i], model.PriceQuota[i]) + + model.EnforceQuotas = pyo.Constraint(model.CROPS, rule=EnforceQuotas_rule) + + # Stage-specific cost computations; + + def ComputeFirstStageCost_rule(model): + return pyo.sum_product(model.PlantingCostPerAcre, model.DevotedAcreage) + model.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule) + + def ComputeSecondStageCost_rule(model): + expr = pyo.sum_product(model.PurchasePrice, model.QuantityPurchased) + expr -= pyo.sum_product(model.SubQuotaSellingPrice, model.QuantitySubQuotaSold) + expr -= pyo.sum_product(model.SuperQuotaSellingPrice, model.QuantitySuperQuotaSold) + return expr + model.SecondStageCost = pyo.Expression(rule=ComputeSecondStageCost_rule) + + def total_cost_rule(model): + if (sense == pyo.minimize): + return model.FirstStageCost + model.SecondStageCost + return -model.FirstStageCost - model.SecondStageCost + model.Total_Cost_Objective = pyo.Objective(rule=total_cost_rule, + sense=sense) + + # Create the list of nodes associated with the scenario (for two stage, + # there is only one node associated with the scenario--leaf nodes are + # ignored). + model._mpisppy_node_list = [ + scenario_tree.ScenarioNode( + name="ROOT", + cond_prob=1.0, + stage=1, + cost_expression=model.FirstStageCost, + nonant_list=[model.DevotedAcreage], + scen_model=model, + ) + ] + + # Add the probability of the scenario + if num_scens is not None: + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" + + return model + + +# begin functions not needed by farmer_cylinders +# (but needed by special codes such as confidence intervals) +#========= +def scenario_names_creator(num_scens, start=None): + # (only for Amalgamator): return the full list of num_scens scenario names + # if start!=None, the list starts with the 'start' labeled scenario + if (start is None): + start = 0 + return [f"scen{i}" for i in range(start, start+num_scens)] + + +#========= +def inparser_adder(cfg): + # add options unique to farmer + #cfg.num_scens_required() Not on the command line for bootstrap. + cfg.add_to_config("crops_multiplier", + description="number of crops will be three times this (default 1)", + domain=int, + default=1) + + cfg.add_to_config("farmer_with_integers", + description="make the version that has integers (default False)", + domain=bool, + default=False) + cfg.add_to_config("yield_cv", + description="approximate farmer crop yield coefficient of variation (default None for unif(0,1) )", + domain=float, + default=None) + + +#========= +def kw_creator(cfg): + # (for Amalgamator): linked to the scenario_creator and inparser_adder + kwargs = {"cfg": cfg} + return kwargs + + +def sample_tree_scen_creator(sname, stage, sample_branching_factors, seed, + given_scenario=None, **scenario_creator_kwargs): + """ Create a scenario within a sample tree. Mainly for multi-stage and simple for two-stage. + (this function supports zhat and confidence interval code) + Args: + sname (string): scenario name to be created + stage (int >=1 ): for stages > 1, fix data based on sname in earlier stages + sample_branching_factors (list of ints): branching factors for the sample tree + seed (int): To allow random sampling (for some problems, it might be scenario offset) + given_scenario (Pyomo concrete model): if not None, use this to get data for ealier stages + scenario_creator_kwargs (dict): keyword args for the standard scenario creator funcion + Returns: + scenario (Pyomo concrete model): A scenario for sname with data in stages < stage determined + by the arguments + """ + # Since this is a two-stage problem, we don't have to do much. + sca = scenario_creator_kwargs.copy() + sca["seed_offset"] = seed + sca["num_scens"] = sample_branching_factors[0] # two-stage problem + return scenario_creator(sname, **sca) + + +# end functions not needed by farmer_cylinders + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + sname = scenario_name + s = scenario + if sname == 'scen0': + print("Arbitrary sanity checks:") + print("SUGAR_BEETS0 for scenario", sname, "is", + pyo.value(s.DevotedAcreage["SUGAR_BEETS0"])) + print("FirstStageCost for scenario", sname, "is", pyo.value(s.FirstStageCost)) + + +#============================ +def xhat_generator(scenario_names, solver_name=None, solver_options=None, cfg=None): + """ Solve the extensive form over the given scenarios and return xhat. + + This is the fixed-name generator the bootstrap code calls when no xhat file + is supplied (see boot_utils.compute_xhat). It builds the EF directly from + this module's scenario_creator so the example is self-contained. + + Args: + scenario_names (list of str): scenarios to build the EF from + solver_name (str): solver to use + solver_options (dict, optional): options passed to the solver + cfg (Config): control parameters (crops_multiplier, yield_cv, ...) + Returns: + xhat (dict): the first-stage nonants keyed by tree node (e.g. ROOT) + """ + ef = sputils.create_EF( + scenario_names, + scenario_creator, + scenario_creator_kwargs={"cfg": cfg}, + ) + solver = pyo.SolverFactory(solver_name) + if solver_options is not None: + for k, v in solver_options.items(): + solver.options[k] = v + if 'persistent' in solver_name: + solver.set_instance(ef, symbolic_solver_labels=True) + solver.solve(tee=False) + else: + solver.solve(ef, tee=False, symbolic_solver_labels=True) + return sputils.nonant_cache_from_ef(ef) diff --git a/examples/bootsp/farmer/smoothed_farmer.json b/examples/bootsp/farmer/smoothed_farmer.json new file mode 100644 index 000000000..82fd55190 --- /dev/null +++ b/examples/bootsp/farmer/smoothed_farmer.json @@ -0,0 +1,21 @@ +{ + "module_name": "farmer", + "max_count": 300, + "candidate_sample_size": 5, + "sample_size": 20, + "subsample_size": 5, + "smoothed_B_I": 5, + "smoothed_center_sample_size": 40, + "nB": 10, + "alpha": 0.05, + "seed_offset": 111, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "trace_fname": "None", + "boot_method": "Smoothed_bagging", + "coverage_replications": 5, + "crops_multiplier": 1, + "farmer_with_integers": "False", + "yield_cv": "0.1" +} diff --git a/examples/bootsp/multi_knapsack/multi_knapsack.bash b/examples/bootsp/multi_knapsack/multi_knapsack.bash new file mode 100644 index 000000000..4b36ab5e9 --- /dev/null +++ b/examples/bootsp/multi_knapsack/multi_knapsack.bash @@ -0,0 +1,22 @@ +#!/bin/bash +# Run the multi-knapsack bootstrap example (needs the statdist library). +# The sample sizes here are small; this is just a demonstration. +# NOTE: do not be alarmed by infeasibility messages during the confidence +# interval calculations. Pass a solver name as the first argument. + +SOLVER=${1:-cplex_direct} +BOOT="python -m mpisppy.confidence_intervals.bootsp.user_boot" +COMMON="--max-count 300 --candidate-sample-size 5 --sample-size 50 \ + --subsample-size 10 --nB 20 --alpha 0.1 --seed-offset 100 \ + --deterministic-data-json multi_knapsack_data.json \ + --solver-name ${SOLVER}" + +echo "Serial, compute xhat within user_boot (empirical Bagging_with_replacement)" +echo +time ${BOOT} multi_knapsack ${COMMON} --boot-method Bagging_with_replacement +echo +echo "========================" +echo +echo "Smoothed coverage simulation from a json file (Smoothed_bagging)" +echo +time python -m mpisppy.confidence_intervals.bootsp.simulate_boot smoothed_multi_knapsack.json diff --git a/examples/bootsp/multi_knapsack/multi_knapsack.json b/examples/bootsp/multi_knapsack/multi_knapsack.json new file mode 100644 index 000000000..3680837af --- /dev/null +++ b/examples/bootsp/multi_knapsack/multi_knapsack.json @@ -0,0 +1,17 @@ +{ + "module_name": "multi_knapsack", + "max_count": 300, + "candidate_sample_size": 5, + "sample_size": 50, + "subsample_size": 10, + "nB": 20, + "alpha": 0.1, + "seed_offset": 100, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "boot_method": "Bagging_with_replacement", + "trace_fname": "None", + "coverage_replications": 5, + "deterministic_data_json": "multi_knapsack_data.json" +} diff --git a/examples/bootsp/multi_knapsack/multi_knapsack.py b/examples/bootsp/multi_knapsack/multi_knapsack.py new file mode 100644 index 000000000..09223bf15 --- /dev/null +++ b/examples/bootsp/multi_knapsack/multi_knapsack.py @@ -0,0 +1,234 @@ +############################################################################### +# 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. +############################################################################### +# A multi-product knapsack example (Vaagen & Wallace, IJPE 2007; the model +# version is from chapter 6 of the King/Wallace book) for the bootstrap +# confidence-interval code. Deterministic data come from a json file named by +# --deterministic-data-json; the random demands are drawn from statdist +# univariate-normal distributions (empirical path) or from distributions fitted +# to the sample data (smoothed path), so importing this example needs statdist. + +import os +import json +import numpy as np +import pyomo.environ as pyo +import mpisppy.scenario_tree as scenario_tree # noqa: F401 (kept for parity/attach_root_node users) +import mpisppy.utils.sputils as sputils +import mpisppy.confidence_intervals.bootsp.statdist as statdist +from mpisppy.confidence_intervals.bootsp.statdist.sampler import Sampler + +# Use this random stream: +sstream = np.random.RandomState(1) + + +def _read_detdata(cfg): + # deterministic data; resolve the file relative to this module if it is not + # found relative to the current working directory + json_fname = cfg.deterministic_data_json + if not os.path.isabs(json_fname) and not os.path.exists(json_fname): + here = os.path.dirname(os.path.abspath(__file__)) + candidate = os.path.join(here, json_fname) + if os.path.exists(candidate): + json_fname = candidate + try: + with open(json_fname, "r") as read_file: + detdata = json.load(read_file) + except Exception: + print(f"Could not read the json file: {json_fname}") + raise + return detdata + + +def _detdata_for(cfg): + # the smoothed driver stashes the parsed data on cfg.detdata; otherwise read + # it from the file (this makes the empirical path work without that stash) + if "detdata" in cfg and cfg.detdata is not None: + return cfg.detdata + return _read_detdata(cfg) + + +def _get_distr_dict(cfg, detdata): + if not getattr(cfg, "use_fitted", False): + unorm = statdist.distribution_factory('univariate-normal') + varset = pyo.RangeSet(detdata["num_prods"]) + distr_dict = {} + for i in varset: + distr_dict[i] = { + "high": unorm(var=(detdata["stdev_d"]["high"])**2, mean=detdata["mean_d"]["high"]), + "low": unorm(var=(detdata["stdev_d"]["low"])**2, mean=detdata["mean_d"]["low"]) + } + else: + distr_dict = cfg.fitted_distribution + return distr_dict + + +def data_sampler(record_num, cfg): + detdata = _detdata_for(cfg) + + distr_dict = _get_distr_dict(cfg, detdata) + sstream.seed(record_num+cfg.seed_offset) + + # this part of the code is the same as in the scenario creator + data = {} + varset = pyo.RangeSet(detdata["num_prods"]) + if getattr(cfg, "use_fitted", False): + for i in varset: + sampler = Sampler([distr_dict[i]], sstream) + data[i] = max(0, int(sampler.sample_one()[0])) + else: + for i in varset: + state = 'high' if sstream.uniform() < 0.5 else 'low' + sampler = Sampler([distr_dict[i][state]], sstream) + data[i] = max(0, int(sampler.sample_one()[0])) + return data + + +def scenario_creator(scenario_name, cfg=None, seed_offset=None, num_scens=None): + """ Create a multi-knapsack scenario. + + Args: + scenario_name (str): + Name of the scenario to construct. + cfg (Config): the control parameters + seed_offset (int): used by confidence interval code + Returns: + model (ConcreteModel): the Pyomo model + """ + # scenario_name has the form e.g. scen12, foobar7 + # The digits are scraped off the right of scenario_name using regex. + scennum = sputils.extract_num(scenario_name) + + seed_offset = cfg.get("seed_offset", 0) if seed_offset is None else seed_offset + sstream.seed(scennum+seed_offset) # allows for resampling easily + num_scens = cfg.get('num_scens', None) + + # Create the concrete model object + model = pyo.ConcreteModel(f"multi-knapsack {scenario_name}") + + detdata = _detdata_for(cfg) + v = detdata["v"] + c = detdata["c"] + g = detdata["g"] + alpha = detdata["alpha"] # a dict of lists + + # use the same variable names as in chapter 6 of the King/Wallace book + # item numbers start at 1 + model.I = pyo.RangeSet(detdata["num_prods"]) + + model.x = pyo.Var(model.I, within=pyo.NonNegativeReals, initialize=0) + model.y = pyo.Var(model.I, within=pyo.NonNegativeReals, initialize=0) + model.z = pyo.Var(model.I, model.I, within=pyo.NonNegativeReals, initialize=0) + model.zt = pyo.Var(model.I, within=pyo.NonNegativeReals, initialize=0) + model.w = pyo.Var(model.I, within=pyo.NonNegativeReals, initialize=0) + + d = data_sampler(scennum, cfg) + + # note: the json indexes are strings + + def d_rule(m, i): + return m.y[i] + sum(m.z[j, i] for j in model.I if j != i) <= d[i] + model.d_constraint = pyo.Constraint(model.I, rule=d_rule) + + def z_rule(m, i, j): + # note that alpha is a dict of lists + if i == j: + return pyo.Constraint.Skip + else: + return m.z[i, j] <= alpha[str(i)][j-1] * (d[j]-m.y[j]) + model.z_constraint = pyo.Constraint(model.I, model.I, rule=z_rule) + + def zt_rule(m, i): + return m.zt[i] == sum(m.z[i, j] for j in model.I if j != i) + model.zt_constraint = pyo.Constraint(model.I, rule=zt_rule) + + def w_rule(m, i): + return m.w[i] == m.x[i] - (m.y[i]+m.zt[i]) + model.w_constraint = pyo.Constraint(model.I, rule=w_rule) + + m = model # typing aid + model.Obj1 = pyo.Expression(expr=-sum(v[str(i)]*(m.y[i]+m.zt[i]) + + g[str(i)]*m.w[i] + - c[str(i)]*m.x[i] for i in m.I)) + + model.obj = pyo.Objective(expr=model.Obj1, sense=pyo.minimize) + + # Create the list of nodes associated with the scenario (for two stage, + # there is only one node associated with the scenario--leaf nodes are + # ignored). + varlist = [model.x] + sputils.attach_root_node(model, model.Obj1, varlist) + + # Add the probability of the scenario + if num_scens is not None: + model._mpisppy_probability = 1/num_scens + else: + model._mpisppy_probability = "uniform" + return model + + +#========= +def scenario_names_creator(num_scens, start=None): + # (only for Amalgamator): return the full list of num_scens scenario names + # if start!=None, the list starts with the 'start' labeled scenario + if (start is None): + start = 0 + return [f"scen{i}" for i in range(start, start+num_scens)] + + +#========= +def inparser_adder(cfg): + # add options unique to the model + cfg.add_to_config("deterministic_data_json", + description="file name for json file with determinstic data", + domain=str, + default=None) + + +#========= +def kw_creator(cfg): + # linked to the scenario_creator and inparser_adder + kwargs = {"cfg": cfg} + return kwargs + + +#============================ +def scenario_denouement(rank, scenario_name, scenario): + pass + + +#============================ +def xhat_generator(scenario_names, solver_name=None, solver_options=None, cfg=None): + """ Solve the extensive form over the given scenarios and return xhat. + + This is the fixed-name generator the bootstrap code calls when no xhat file + is supplied (see boot_utils.compute_xhat). It builds the EF directly from + this module's scenario_creator so the example is self-contained. + + Args: + scenario_names (list of str): scenarios to build the EF from + solver_name (str): solver to use + solver_options (dict, optional): options passed to the solver + cfg (Config): control parameters (includes deterministic_data_json) + Returns: + xhat (dict): the first-stage nonants keyed by tree node (e.g. ROOT) + """ + ef = sputils.create_EF( + scenario_names, + scenario_creator, + scenario_creator_kwargs={"cfg": cfg}, + ) + solver = pyo.SolverFactory(solver_name) + if solver_options is not None: + for k, v in solver_options.items(): + solver.options[k] = v + if 'persistent' in solver_name: + solver.set_instance(ef, symbolic_solver_labels=True) + solver.solve(tee=False) + else: + solver.solve(ef, tee=False, symbolic_solver_labels=True) + return sputils.nonant_cache_from_ef(ef) diff --git a/examples/bootsp/multi_knapsack/multi_knapsack_data.json b/examples/bootsp/multi_knapsack/multi_knapsack_data.json new file mode 100644 index 000000000..88a98f9f5 --- /dev/null +++ b/examples/bootsp/multi_knapsack/multi_knapsack_data.json @@ -0,0 +1,85 @@ +{ + "c": { + "1": 6.270388712278181, + "3": 3.8127756332188465, + "2": 2.4265018086314196, + "5": 5.112747213686085, + "4": 2.5891675029296337, + "6": 4.037494846539122 + }, + "g": { + "1": 1.351074962440077, + "3": 0.672914529329352, + "2": 1.212727044704484, + "5": 0.8180395541897737, + "4": 0.4142668004687414, + "6": 0.647894619920663 + }, + "mean_d": { + "high": 1160, + "low": 116 + }, + "num_prods": 6, + "v": { + "1": 16.888437030500963, + "3": 8.4114316166169, + "2": 15.159088058806049, + "5": 10.22549442737217, + "4": 5.178335005859267, + "6": 8.098682749008287 + }, + "stdev_d": { + "high": 74, + "low": 96 + }, + "alpha": { + "1": [ + 0, + 0.1, + 0.1, + 0.1, + 0.1, + 0.1 + ], + "3": [ + 0.1, + 0.1, + 0, + 0.1, + 0.1, + 0.1 + ], + "2": [ + 0.1, + 0, + 0.1, + 0.1, + 0.1, + 0.1 + ], + "5": [ + 0.1, + 0.1, + 0.1, + 0.1, + 0, + 0.1 + ], + "4": [ + 0.1, + 0.1, + 0.1, + 0, + 0.1, + 0.1 + ], + "6": [ + 0.1, + 0.1, + 0.1, + 0.1, + 0.1, + 0 + ] + } +} \ No newline at end of file diff --git a/examples/bootsp/multi_knapsack/smoothed_multi_knapsack.json b/examples/bootsp/multi_knapsack/smoothed_multi_knapsack.json new file mode 100644 index 000000000..4b7be189c --- /dev/null +++ b/examples/bootsp/multi_knapsack/smoothed_multi_knapsack.json @@ -0,0 +1,19 @@ +{ + "module_name": "multi_knapsack", + "max_count": 300, + "candidate_sample_size": 5, + "sample_size": 20, + "subsample_size": 5, + "smoothed_B_I": 5, + "smoothed_center_sample_size": 40, + "nB": 10, + "alpha": 0.05, + "seed_offset": 111, + "optimal_fname": "None", + "xhat_fname": "None", + "solver_name": "cplex_direct", + "trace_fname": "None", + "boot_method": "Smoothed_bagging", + "coverage_replications": 5, + "deterministic_data_json": "multi_knapsack_data.json" +} diff --git a/examples/run_all.py b/examples/run_all.py index ecfb49bfa..502579a88 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -135,8 +135,9 @@ def do_one_mmw(dirname, modname, runefstring, npyfile, mmwargstring): os.chdir("..") # moved to CI directory def do_one_boot(dirname, module, boot_method, size_args, np=2): - # A small bootstrap confidence-interval run on a statdist-free example - # (the other bootstrap examples need statdist, which is merged separately). + # A small bootstrap confidence-interval run. schultz/schultz_data need only + # numpy; farmer/cvar/multi_knapsack pull in statdist (and so scipy), and a + # Smoothed_* method fits a distribution with it. # xhat is computed by the model's xhat_generator (no npy file needed). argstring = (f"{module} {size_args} --alpha 0.1 --seed-offset 100 " f"--solver-name {solver_name} --boot-method {boot_method}") @@ -154,6 +155,16 @@ def do_one_boot(dirname, module, boot_method, size_args, np=2): do_one_boot("schultz_data", "schultz_data", "Bagging_with_replacement", "--max-count 200 --candidate-sample-size 5 --sample-size 100 " "--subsample-size 20 --nB 20", np=2) + # farmer: an empirical run on a statdist-dependent example (crop yields are + # perturbed by a statdist univariate distribution) + do_one_boot("farmer", "farmer", "Bagging_with_replacement", + "--max-count 200 --candidate-sample-size 5 --sample-size 30 " + "--subsample-size 10 --nB 8 --crops-multiplier 1 --yield-cv 0.1", np=2) + # cvar: a smoothed run (statdist fits a kernel density to the sampled data) + do_one_boot("cvar", "cvar", "Smoothed_bagging", + "--max-count 200 --candidate-sample-size 5 --sample-size 20 " + "--subsample-size 5 --nB 8 --smoothed-B-I 3 " + "--smoothed-center-sample-size 20", np=2) do_one("farmer/CI", "farmer_ef.py", 1, "1 3 {}".format(solver_name)) # for farmer_cylinders, the first arg is num_scens and is required diff --git a/mpisppy/confidence_intervals/bootsp/__init__.py b/mpisppy/confidence_intervals/bootsp/__init__.py index 51220c4b5..95385ec80 100644 --- a/mpisppy/confidence_intervals/bootsp/__init__.py +++ b/mpisppy/confidence_intervals/bootsp/__init__.py @@ -7,5 +7,5 @@ # full copyright and license information. ############################################################################### # Bootstrap and bagging confidence intervals for data-based, two-stage -# stochastic programs (the empirical methods; smoothed methods and the -# statdist distribution library arrive in a follow-on merge). +# stochastic programs: the empirical methods (numpy only) and the smoothed +# methods (which fit a distribution with the bundled statdist library). diff --git a/mpisppy/confidence_intervals/bootsp/boot_sp.py b/mpisppy/confidence_intervals/bootsp/boot_sp.py index 699c86cd1..e262cca3b 100644 --- a/mpisppy/confidence_intervals/bootsp/boot_sp.py +++ b/mpisppy/confidence_intervals/bootsp/boot_sp.py @@ -27,6 +27,28 @@ rankcomm = boot_utils.rankcomm +_MAXIMIZATION_MSG = ( + "the bootstrap confidence intervals are minimization-only, but {what} has a " + "maximization objective. The estimators form every gap as (value at xhat) " + "minus the optimal, which is non-positive for a maximization, and the " + "drivers floor the reported interval's lower end at 0 -- so a maximization " + "run would not merely be wrong, it would report [0, 0]. Supporting " + "maximization means deciding how the gap is reported (mpi-sppy's MMW " + "estimator reports its magnitude) and threading the sense through every " + "estimator, the interval floors and the coverage checks.") + + +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. + """ + if not is_minimizing: + raise ValueError(_MAXIMIZATION_MSG.format(what=what)) + + 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: @@ -160,6 +182,7 @@ def solve_routine(cfg, module, scenarios, num_threads=None, duplication=False): scenario_creator, scenario_creator_kwargs=scenario_creator_kwargs, ) + _require_minimization(ef.EF_Obj.sense == pyo.minimize, "this model") solver = pyo.SolverFactory(cfg.solver_name) solver.options["threads"] = num_threads @@ -684,14 +707,18 @@ def compute_ci(cfg, module, xhat): the ci_* entries are None on MPI ranks other than 0. Note: - This is the single dispatch point shared by user_boot and - simulate_boot. A smoothed method raises a friendly "not yet merged" - error (the smoothed methods land in a follow-on merge). + This is the empirical dispatch point shared by user_boot and + simulate_boot. The smoothed methods have a different (gap-only) return + signature and are dispatched by smoothed_boot_sp.compute_smoothed_ci; + a smoothed method reaching here is an error. """ method = cfg.boot_method boot_utils.BootMethods.check_for_it(method) if boot_utils.is_smoothed(method): - boot_utils.smoothed_not_yet_merged(method) + raise ValueError( + 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 method == "Extended": return extended_bootstrap(cfg, module, xhat) elif method == "Bagging_with_replacement": diff --git a/mpisppy/confidence_intervals/bootsp/boot_utils.py b/mpisppy/confidence_intervals/bootsp/boot_utils.py index 3f460b80f..018e98846 100644 --- a/mpisppy/confidence_intervals/bootsp/boot_utils.py +++ b/mpisppy/confidence_intervals/bootsp/boot_utils.py @@ -58,24 +58,13 @@ def is_smoothed(boot_method): def empirical_members(): - """ The BootMethods tokens that are available now (the empirical ones). """ + """ The BootMethods tokens that use only the empirical (statdist-free) code. """ return [m for m in BootMethods.list_of_members() if not is_smoothed(m)] -def smoothed_not_yet_merged(boot_method): - """ Raise a friendly error for a smoothed method that is not merged yet. - - The smoothed bootstrap methods depend on the statdist distribution - library and are being merged separately. Until they land, the empirical - methods are available here and the full set lives in the boot-sp package. - """ - raise RuntimeError( - f"boot_method={boot_method} is a smoothed method, which is not yet " - "available in mpi-sppy (it arrives in a follow-on merge along with " - "the statdist distribution library). Use one of the empirical " - f"methods {empirical_members()} here, or the smoothed methods in the " - "separate boot-sp package (https://github.com/boot-sp/boot-sp)." - ) +def smoothed_members(): + """ The BootMethods tokens that use the statdist smoothed code. """ + return [m for m in BootMethods.list_of_members() if is_smoothed(m)] def module_name_to_module(module_name): diff --git a/mpisppy/confidence_intervals/bootsp/simulate_boot.py b/mpisppy/confidence_intervals/bootsp/simulate_boot.py index 75179866b..b73533209 100644 --- a/mpisppy/confidence_intervals/bootsp/simulate_boot.py +++ b/mpisppy/confidence_intervals/bootsp/simulate_boot.py @@ -12,9 +12,11 @@ # python -m mpisppy.confidence_intervals.bootsp.simulate_boot import sys +import time import mpisppy.confidence_intervals.ciutils as ciutils import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils import mpisppy.confidence_intervals.bootsp.boot_sp as boot_sp +import mpisppy.confidence_intervals.bootsp.smoothed_boot_sp as smoothed_boot_sp my_rank = boot_utils.my_rank @@ -73,14 +75,88 @@ def empirical_main_routine(cfg, module): return None, None +def smoothed_main_routine(cfg, module): + """ The smoothed-method coverage harness; called by main() and test drivers. + + Args: + cfg (Config): parameters + module (Python module): contains the scenario creator function and helpers + Returns: + (coverage_two_sided, coverage_one_sided, ci_lengths, run_times); + all None on MPI ranks other than 0. + + Note: + The smoothed estimators report only the optimality-gap interval, so the + coverage counts are against opt_gap (from process_optimal) rather than + against z* as in the empirical harness. + """ + if my_rank == 0: + # only opt_gap is used by the smoothed coverage counting + _, opt_gap = boot_sp.process_optimal(cfg, module) + else: + opt_gap = None + + if cfg["xhat_fname"] is not None and cfg["xhat_fname"] != "None": + xhat = ciutils.read_xhat(cfg["xhat_fname"]) + else: + # boot-sp called an undefined fit_resample_utils.compute_xhat here; the + # intended call is boot_utils.compute_xhat (design doc section 4.3). + xhat = boot_utils.compute_xhat(cfg, module) + + coverage_cnt_one_sided, coverage_cnt_two_sided = 0, 0 + ci_len = [] + run_time = [] + seed_offset = cfg.seed_offset # store the original offset + # Replications must not share draws. Unlike the empirical methods, which + # get independent streams straight from numpy (default_rng([seed_offset, + # word]) and so can step the offset by 1), the smoothed methods address + # their draws by record number -- the model seeds each draw with it -- so + # replications are separated by giving each one its own block of record + # numbers. The stride is therefore exactly what one replication consumes: + # the center block plus one block per batch (smoothed bootstrap), or B_I + # groups of nB bags (smoothed bagging). Take the larger of the two so the + # stride covers whichever method is running. + stride = max((cfg.smoothed_center_sample_size or 0) + cfg.nB * cfg.sample_size, + (cfg.smoothed_B_I or 1) * cfg.nB * cfg.subsample_size) + seed_list = [i * stride + seed_offset for i in range(cfg.coverage_replications)] + + for seed in seed_list: + cfg.seed_offset = seed + if my_rank == 0: + st_time = time.time() + ci_gap_two_sided, _ = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + if my_rank == 0: + en_time = time.time() + if cfg.trace_fname is not None: + with open(cfg.trace_fname, "a+") as f: + f.write(f"seed: {cfg.seed_offset}\n") + f.write(f"optimality gap: {opt_gap}\n") + f.write(f"ci for optimality gap: {ci_gap_two_sided}\n") + if (ci_gap_two_sided[0] <= opt_gap) and (opt_gap <= ci_gap_two_sided[1]): + coverage_cnt_two_sided += 1 + if (opt_gap <= ci_gap_two_sided[1]): + coverage_cnt_one_sided += 1 + ci_len.append(ci_gap_two_sided[1] - ci_gap_two_sided[0]) + run_time.append(en_time - st_time) + + if my_rank == 0: + assert cfg.coverage_replications != 0 + return (coverage_cnt_two_sided / cfg.coverage_replications, + coverage_cnt_one_sided / cfg.coverage_replications, + ci_len, run_time) + else: + return None, None, None, None + + def main(cfg, module): """ Dispatch to the appropriate coverage harness for cfg.boot_method. - A smoothed method raises a friendly "not yet merged" error; the empirical - methods run the empirical coverage harness. + The empirical methods run the empirical coverage harness (returns a + (rate, length) pair); the smoothed methods run the smoothed harness + (returns a (cov_two, cov_one, ci_lengths, run_times) tuple). """ if boot_utils.is_smoothed(cfg.boot_method): - boot_utils.smoothed_not_yet_merged(cfg.boot_method) + return smoothed_main_routine(cfg, module) return empirical_main_routine(cfg, module) diff --git a/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py new file mode 100644 index 000000000..43fad1daa --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py @@ -0,0 +1,336 @@ +############################################################################### +# 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. +############################################################################### +# Smoothed bootstrap/bagging for data-based, two-stage stochastic programs. +# These methods fit a (univariate) distribution to the sampled data using the +# statdist library and then resample from the fitted distribution. They are the +# counterpart to the empirical methods in boot_sp.py. + +import json + +import numpy as np +from numpy.random import default_rng +from statistics import NormalDist +import pyomo.environ as pyo + +from mpisppy import global_toc +import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils +import mpisppy.confidence_intervals.bootsp.boot_sp as boot_sp +import mpisppy.confidence_intervals.bootsp.statdist as statdist + +# The communicators live in boot_utils so there is a single source of truth. +comm = boot_utils.comm +n_proc = boot_utils.n_proc +my_rank = boot_utils.my_rank +rankcomm = boot_utils.rankcomm + + +def fit_distribution(sample_data, distr_type='univariate-epispline'): + """ Fit a (univariate) distribution to sample data. + + Args: + sample_data (list or list of dict): a list of scalars (one variable) or + a list of dicts (multivariate, keyed by variable name) + distr_type (str): a statdist univariate distribution token + Returns: + the fitted distribution (or a dict of them, keyed as the input dicts) + """ + distr_func = statdist.distribution_factory(distr_type) + if isinstance(sample_data[0], (float, int)): # 1-dim + fitted_distr = distr_func.fit(sample_data) + else: + fitted_distr = {} + for key in sample_data[0]: + data = [data_dict[key] for data_dict in sample_data] + fitted_distr[key] = distr_func.fit(data) + return fitted_distr + + +def center_smoothed(cfg, module, xhat): + """ Estimate the CI center (the optimality gap) from the fitted distribution. + + The smoothed methods are single-rank-per-solve, so the solves here go to the + module globals' view; there is no communicator to thread through (an earlier + signature took one and ignored it). + """ + assert cfg.smoothed_center_sample_size is not None, \ + "need a sample size for smoothed bootstrap center estimation" + scenario_pool = list(range(cfg.seed_offset, + cfg.seed_offset + cfg.smoothed_center_sample_size)) + + center_upper = boot_sp.evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False) + center_ef = boot_sp.solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=False) + center_optimal = pyo.value(center_ef.EF_Obj) + center_gap = center_upper - center_optimal + + if my_rank == 0: + return center_gap + else: + return None + + +def smoothed_resample_helper(cfg, module, xhat, serial=False): + """ Get local gaps for the smoothed bootstrap (the fitted-distribution + analog of boot_sp._bootstrap_resample). + + Every batch is an *independent* set of cfg.subsample_size draws from the + fitted distribution, so the batches take disjoint blocks of the draw index + space: batch b covers [start + b*m, start + (b+1)*m). A record number is + the draw's seed, so striding by anything less than m (the batch size) would + hand consecutive batches most of the same draws and collapse the spread the + interval is built from. The block the center estimate draws from, + [seed_offset, seed_offset + smoothed_center_sample_size) (see + center_smoothed), is reserved ahead of the batches because it samples the + same fitted distribution and must not reuse their draws. + """ + if serial: + local_nB = cfg.nB + else: + local_nB = boot_sp.slice_lens(cfg.nB)[my_rank] + + local_boot_gaps = np.empty(local_nB, dtype=np.float64) + + m = cfg.subsample_size + start = cfg.seed_offset + (cfg.smoothed_center_sample_size or 0) + # this rank's first batch in the global 0..nB-1 numbering + first_batch = 0 if serial else sum(boot_sp.slice_lens(cfg.nB)[:my_rank]) + + for iter in range(local_nB): + b = first_batch + iter + scenario_pool = list(range(start + b * m, start + (b + 1) * m)) + + local_boot_upper = boot_sp.evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False) + local_boot_ef = boot_sp.solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=False) + local_boot_optimal = pyo.value(local_boot_ef.EF_Obj) + local_boot_gaps[iter] = local_boot_upper - local_boot_optimal + + return local_boot_gaps + + +def smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline', quantile=False, serial=False): + """ fit a distribution to the sample, then draw both the center and the + batches from it to get a smoothed point estimate and interval width + 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) + distr_type (str): a statdist univariate distribution token to fit + quantile (bool): use the quantile method (else the gaussian method) + serial (bool): indicates that only one MPI rank should be used + Returns: + tuple (ci_gap_two_sided, center_gap) if on MPI rank 0, else None + + """ + rng = default_rng(cfg.seed_offset) + scenario_pool = rng.choice(cfg.max_count, size=cfg.sample_size, replace=False) + + cfg.use_fitted = False + sample_data = [module.data_sampler(scenario, cfg) for scenario in scenario_pool] + cfg.fitted_distribution = fit_distribution(sample_data, distr_type=distr_type) + # From here on both the center and the batches draw from the fitted + # distribution. Estimating the center from the raw sample instead would + # make it the purely empirical point estimate and give up the smoothing + # that is the whole point of these methods. + cfg.use_fitted = True + + # the center: one replication at a large resample size + # (smoothed_center_sample_size) drawn from the fitted distribution + dag_gap = center_smoothed(cfg, module, xhat) + comm.Barrier() + + # each batch is a fresh set of cfg.sample_size draws from the same fitted + # distribution, so the bootstrap batch size is the full sample size + cfg.subsample_size = cfg.sample_size + local_boot_gaps = smoothed_resample_helper(cfg, module, xhat, serial) + comm.Barrier() + + # do analysis only on rank 0 + if my_rank == 0: + boot_gap = np.empty(cfg.nB, dtype=np.float64) + else: + boot_gap = None + + # but everyone needs to send to the gather + lenlist = boot_sp.slice_lens(cfg.nB) + comm.Gatherv(sendbuf=local_boot_gaps, recvbuf=(boot_gap, lenlist), root=0) + + if my_rank == 0: + global_toc("Done smoothed bootstrap") + + if not quantile: + s_g = np.std(boot_gap, ddof=1) + ppf = NormalDist().inv_cdf(1 - cfg.alpha / 2) + error = s_g * ppf + ci_gap_two_sided = [dag_gap - error, dag_gap + error] + else: + alpha = cfg.alpha / 2 + eps = np.quantile(boot_gap - dag_gap, [alpha, 1 - alpha]) + ci_gap_two_sided = [dag_gap - eps[1], dag_gap - eps[0]] + print(f"{ci_gap_two_sided = }") + return ci_gap_two_sided, dag_gap + else: + # non-root ranks return a matching arity so callers can unpack safely + return None, None + + +def smoothed_bagging(cfg, module, xhat, distr_type='univariate-kernel', serial=False): + """ perform a bagging-based estimation of confidence intervals using a fitted distribution + 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) + distr_type (str): a statdist univariate distribution token to fit + serial (bool): indicates that only one MPI rank should be used + Returns: + tuple (ci_gap_two_sided, center_gap) if on MPI rank 0, else None + """ + rng = default_rng(cfg.seed_offset) + scenario_pool = rng.choice(cfg.max_count, size=cfg.sample_size, replace=False) + + cfg.use_fitted = False + sample_data = [module.data_sampler(scenario, cfg) for scenario in scenario_pool] + cfg.fitted_distribution = fit_distribution(sample_data, distr_type=distr_type) + cfg.use_fitted = True + + local_nB = boot_sp.slice_lens(cfg.nB)[my_rank] + local_gaps = np.empty(local_nB, dtype=np.float64) + + if my_rank == 0: + bagging_gap = np.empty(cfg.nB, dtype=np.float64) + all_gaps = [] + avg_gaps = [] + else: + bagging_gap = None + all_gaps = None + avg_gaps = None + + # B_I is the number of initial seed points; s1 below is the variance *among* + # their averages, so fewer than two of them leaves it undefined + if cfg.smoothed_B_I is None or cfg.smoothed_B_I < 2: + raise ValueError( + "smoothed_B_I (the number of initial seed points) must be at least " + f"2 for smoothed bagging; got {cfg.smoothed_B_I}. The variance of " + "the per-seed-point averages is what estimates the between-point " + "term of the interval width.") + + B_I = cfg.smoothed_B_I + for i in range(B_I): + seed_offset_base = cfg.seed_offset + cfg.nB * cfg.subsample_size * i + + for j in range(local_nB): + seed_offset = seed_offset_base + (sum(boot_sp.slice_lens(cfg.nB)[:my_rank]) + j) * cfg.subsample_size + scenario_pool = list(range(seed_offset, seed_offset + cfg.subsample_size)) + scenario_pool[0] = seed_offset_base + + local_upper = boot_sp.evaluate_scenarios(cfg, module, scenario_pool, xhat, duplication=False) + local_ef = boot_sp.solve_routine(cfg, module, scenario_pool, num_threads=2, duplication=False) + local_optimal = pyo.value(local_ef.EF_Obj) + local_gaps[j] = local_upper - local_optimal + comm.Barrier() + lenlist = boot_sp.slice_lens(cfg.nB) + comm.Gatherv(sendbuf=local_gaps, recvbuf=(bagging_gap, lenlist), root=0) + + if my_rank == 0: + all_gaps = all_gaps + bagging_gap.tolist() + avg_gaps.append(np.mean(bagging_gap)) + + if my_rank == 0: + global_toc("Done Smoothed Bagging") + + dag_gap = np.mean(avg_gaps) + + # sample variances (ddof=1), as the algorithm specifies and as the + # empirical estimators already use for their Gaussian half-width; + # ddof=0 understates s1 by (B_I-1)/B_I, which is a third at B_I=3 + s1 = np.var(avg_gaps, ddof=1) + s2 = np.var(all_gaps, ddof=1) + ppf = NormalDist().inv_cdf(1 - cfg.alpha / 2) + s_g_2 = (cfg.subsample_size**2) * s1 / cfg.sample_size + s2 / (B_I * cfg.nB) + error = np.sqrt(s_g_2) * ppf + ci_gap_two_sided = [dag_gap - error, dag_gap + error] + + print(f"{ci_gap_two_sided = }") + return ci_gap_two_sided, dag_gap + else: + # non-root ranks return a matching arity so callers can unpack safely + return None, None + + +def _ensure_smoothed_cfg(cfg): + """ Idempotently attach the run-time config entries the smoothed methods need. + + The smoothed estimators toggle ``use_fitted`` and stash a + ``fitted_distribution`` on the cfg; a module may also supply deterministic + data via a json file named by ``deterministic_data_json``. This may be + called repeatedly (e.g. once per replication in a coverage simulation), so + every add is guarded. + """ + if "use_fitted" not in cfg: + cfg.add_to_config(name="use_fitted", + description="a boolean to control use of fitted distribution", + domain=bool, + default=None, + argparse=False) + cfg.use_fitted = False + if "fitted_distribution" not in cfg: + cfg.add_to_config(name="fitted_distribution", + description="a fitted distribution from sample data", + domain=None, + default=None, + argparse=False) + if "deterministic_data_json" in cfg and "detdata" not in cfg: + json_fname = cfg.deterministic_data_json + try: + with open(json_fname, "r") as read_file: + detdata = json.load(read_file) + except Exception: + print(f"Could not read the json file: {json_fname}") + raise + cfg.add_to_config("detdata", + description="deterministic data from json file", + domain=dict, + default=detdata) + + +def compute_smoothed_ci(cfg, module, xhat): + """ Dispatch to the requested smoothed bootstrap/bagging method. + + 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 + Returns: + (ci_gap_two_sided, center_gap) on MPI rank 0, else None + + Note: + This is the single smoothed-dispatch point shared by user_boot and + simulate_boot (the counterpart to boot_sp.compute_ci for the empirical + methods). + """ + _ensure_smoothed_cfg(cfg) + method = cfg.boot_method + boot_utils.BootMethods.check_for_it(method) + if method == "Smoothed_boot_epi": + return smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline') + elif method == "Smoothed_boot_kernel": + return smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-kernel') + elif method == "Smoothed_boot_epi_quantile": + return smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline', quantile=True) + elif method == "Smoothed_boot_kernel_quantile": + return smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-kernel', quantile=True) + elif method == "Smoothed_bagging": + return smoothed_bagging(cfg, module, xhat, distr_type='univariate-kernel') + else: + raise ValueError(f"boot_method={method} is not a smoothed method.") + + +if __name__ == "__main__": + print("smoothed_boot_sp contains only functions and is not directly runnable.") + print("Try, e.g., user_boot.py") diff --git a/mpisppy/confidence_intervals/bootsp/statdist/README.md b/mpisppy/confidence_intervals/bootsp/statdist/README.md new file mode 100644 index 000000000..2c2146554 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/README.md @@ -0,0 +1,36 @@ +# statdist (trimmed) + +This is a trimmed port of the **statdist** statistical-distribution library, +bundled here for the smoothed bootstrap methods in +`mpisppy.confidence_intervals.bootsp`. + +## What is here + +Only the **univariate** distributions and their support modules: + +- `base_distribution.py` — the distribution base classes and helpers +- `distributions.py` — the univariate distribution classes (uniform, normal, + student-t, Gaussian kernel, epi-spline, empirical, discrete) +- `distribution_factory.py` — name → class registry (`distribution_factory`) +- `splines.py` — epi-spline fitting (builds a small Pyomo model) +- `utilities.py`, `sampler.py` — memoization/context helpers and the sampler + +## What was dropped + +The multivariate machinery — `copula.py`, `vine.py`, `bicop.py`, the +multivariate distribution classes in `distributions.py`, and the +`MultivariateDistribution` base class and its decorators in +`base_distribution.py` — is **not** included. +Dropping it also removes the `from scipy.stats import mvn` import (removed in +scipy 1.14) and the optional `gosm` hook, neither of which the smoothed +bootstrap methods use. scipy is imported lazily (via +`pyomo.common.dependencies`) so the empirical bootstrap path stays scipy-free. + +The full library, including the multivariate code, lives in the archived +boot-sp repository: https://github.com/boot-sp/boot-sp + +## Provenance + +statdist was developed under separate funding, always intended to be +open-source, and shares lineage with the GOSM/Prescient scenario-generation +tools. diff --git a/mpisppy/confidence_intervals/bootsp/statdist/__init__.py b/mpisppy/confidence_intervals/bootsp/statdist/__init__.py new file mode 100644 index 000000000..d176e5f0e --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/__init__.py @@ -0,0 +1,13 @@ +############################################################################### +# 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. +############################################################################### +# Trimmed statdist: univariate distributions only (see README.md). The +# distribution_factory re-export lets callers write statdist.distribution_factory(...). + +from mpisppy.confidence_intervals.bootsp.statdist.distribution_factory import distribution_factory # noqa: F401 +from mpisppy.confidence_intervals.bootsp.statdist.distributions import * # noqa: F401,F403 diff --git a/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py new file mode 100644 index 000000000..9f9a6708d --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py @@ -0,0 +1,505 @@ +############################################################################### +# 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. +############################################################################### +""" +This abstract base class is the parent class of all distribution classes. +""" +from abc import ABCMeta, abstractmethod +import os + +import numpy as np +# scipy is an optional dependency; import it lazily so the empirical +# bootstrap path stays scipy-free. matplotlib (also optional) is imported +# locally in the plotting methods below. +from pyomo.common.dependencies import scipy + +from mpisppy.confidence_intervals.bootsp.statdist.utilities import memoize_method + +class Parameter: + """ + This class will encode the information to fully specify a parameter for + a distribution. It will have a name, a value, bounds on what the value + can be, and the type of the value. + + Attributes: + name (str): The name of the parameter + value (float): The value of the parameter, if None, the parameter + is not instantiated + bounds (tuple): An ordered pair (a, b) specifying the lower and upper + bounds of the value inclusive, either may be None to specify a + lack of bound. + kind (type): The type of value the parameter has + """ + def __init__(self, name, value=None, bounds=(None, None), kind=float): + """ + Args: + name (str): The name of the parameter + value (float): The value of the parameter, if None, the parameter + is not instantiated + bounds (tuple): An ordered pair (a, b) specifying the lower and + upper bounds of the value inclusive. Either may be None to + specify a lack of bound. + kind (type): The type of value the parameter has + + """ + self.name = name + self.value = value + # a parameter is instantiated once it has a value (the docstring above: + # "if None, the parameter is not instantiated") + self.instantiated = value is not None + self.bounds = bounds + self.kind = kind + + def set_value(self, value): + """ + Sets the value attribute. + Args: + value: The value to set the parameter to + """ + self.value = value + self.instantiated = value is not None + + def __repr__(self): + return "Parameter({},{})".format(self.name, self.value) + + __str__ = __repr__ + + +class BaseDistribution(object): + __metaclass__ = ABCMeta + + # -------------------------------------------------------------------- + # Abstract methods (have to be implemented within the subclass) + # -------------------------------------------------------------------- + + @abstractmethod + def __init__(self, dimension=0, parameters=None): + """ + Initializes the distribution. + + Args: + dimension (int): the dimension of the distribution + parameters (list[Parameter]): A list of parameters for the + distribution + """ + self.name = self.__class__.__name__ + self.dimension = dimension + + self.parameters = parameters if parameters else [] + + @abstractmethod + def pdf(self, x): + """ + Evaluates the probability density function at a given point x. + + Args: + x (float): the point at which the pdf is to be evaluated + + Returns: + float: the value of the pdf + """ + pass + + @classmethod + @abstractmethod + def fit(cls, data): + """ + This function will fit the distribution of this class to the passed + in data. This will return an instance of the class. + + Args: + data (List[float]): The data the distribution is to be fit to + Returns: + baseDistribution: The fitted distribution + """ + pass + + @staticmethod + def seed_reset(seed=None): + """ + Resets the random seed for sampling. + If no argument is passed, the current time is used. + + Args: + seed: the random seed + """ + np.random.seed(seed) + + def __str__(self): + string = self.name + ': ' + for parameter in self.parameters: + string += '\n{}: {}'.format(parameter.name, parameter.value) + return string + + def __repr__(self): + return "Distribution({})".format(self.name) + + +class UnivariateDistribution(BaseDistribution): + """ + This is the base for all univariate distributions. It will have specialized + pdf and cdf methods which take a single argument. + """ + __metaclass__ = ABCMeta + + def __init__(self, parameters=None, lower=None, upper=None): + """ + Args: + parameters (list[Parameter]): A list of parameters for the + distribution + """ + if lower is None: + self.lower = -np.inf + else: + self.lower = lower + + if upper is None: + self.upper = np.inf + else: + self.upper = upper + + BaseDistribution.__init__(self, 1, parameters) + + def plot(self, plot_pdf=True, plot_cdf=True, output_file=None, title=None, + xlabel=None, ylabel=None, output_directory='.'): + """ + Plots the pdf/cdf within the interval [alpha, beta]. + If no output file is specified, the plots are shown at + runtime. + + Args: + plot_pdf (bool): True if the plot should include the pdf + plot_cdf (bool): True if the plot should include the cdf + output_file (str): name of an output file to save the plot + title (str): the title of the plot + xlabel (str): the name of the x-axis + ylabel (str): the name of the y-axis + output_directory (str): The name of the directory to save the + files, defaults to the current working directory + """ + if self.lower == -np.inf: + lower = -5 + else: + lower = self.lower + + if self.upper == np.inf: + upper = 5 + else: + upper = self.upper + + directory = output_directory + try: + os.makedirs(directory) + except FileExistsError: + pass + + x_range = np.linspace(lower, upper, 100) + import matplotlib.pyplot as plt + fig = plt.figure() + + # Plot the pdf if required. + if plot_pdf: + y_range = [] + for x in x_range: + y_range.append(self.pdf(x)) + plt.plot(x_range, y_range, label='PDF', color='blue') + + # Plot the cdf if required. + if plot_cdf: + y_range = [] + for x in x_range: + y_range.append(self.cdf(x)) + plt.plot(x_range, y_range, label='CDF', color='red') + + # Display a legend. + lgd = plt.legend(loc='lower center', bbox_to_anchor=(0.5, -0.25), + ncol=3, shadow=True) + + # Display a grid and the axes. + plt.grid(True, which='both') + plt.axhline(y=0, color='k') + plt.axvline(x=0, color='k') + + # Name the axes. + plt.xlabel(xlabel) + plt.ylabel(ylabel) + + plt.title(title, y=1.08) + + if output_file is None: + # Display the plot. + plt.show() + else: + # Save the plot. + plt.savefig(directory + os.sep + output_file, + bbox_extra_artists=(lgd,), bbox_inches='tight') + + plt.close(fig) + + @memoize_method + def cdf(self, x, epsabs=1e-4): + """ + Evaluates the cumulative distribution function at a given point x. + + Args: + x (float): the point at which the cdf is to be evaluated + epsabs (float): The accuracy to which the cdf is to be calculated + + Returns: + float: the value of the cdf + """ + if x <= self.alpha: + return 0 + elif x >= self.beta: + return 1 + else: + return scipy.integrate.quad(self.pdf, self.alpha, x, epsabs=epsabs)[0] + + @memoize_method + def cdf_inverse(self, x, cdf_inverse_tolerance=1e-4, + cdf_inverse_max_refinements=10, + cdf_tolerance=1e-4): + """ + Evaluates the inverse cumulative distribution function at a given + point x. + + TODO: Explain better how this is calculated + + Args: + x (float): the point at which the inverse cdf is to be evaluated + cdf_inverse_tolerance (float): The accuracy which the inverse + cdf is to be calculated to + cdf_inverse_max_refinements (int): The number of times the + the partition on the x-domain will be made finer + cdf_tolerance (float): The accuracy to which the cdf is calculated + to + Returns: + float: the value of the inverse cdf + """ + + # For ease in calculating the cdf, we define this temp function. + cdf = lambda x: self.cdf(x, epsabs=cdf_tolerance) + + # This method calculates the cdf of start and then increases + # (if the cdf value is less than or equal x) or decreases + # (if the cdf value is greater than x) start iteratively by one + # stepsize until x is passed. It returns the increased (or decreased) + # start value and its cdf value. + def approximate_inverse_value(start): + cdf_val = cdf(start) + if x >= cdf_val: + while x >= cdf_val: + start += stepsize + cdf_val = cdf(start) + else: + while x <= cdf_val: + start -= stepsize + cdf_val = cdf(start) + return cdf_val, start + + # Handle some special cases. + if x < 0 or x > 1: + return None + elif abs(x) <= cdf_inverse_tolerance: + return self.alpha + elif abs(x-1) <= cdf_inverse_tolerance: + return self.beta + else: + + # Initialize variables. + approx_x = 0 + result = None + number_of_refinement = 0 + + # The starting stepsize was chosen arbitrarily. + stepsize = (self.beta - self.alpha)/10 + + while abs(approx_x - x) > cdf_inverse_tolerance \ + and number_of_refinement <= cdf_inverse_max_refinements: + + # If this is the first iteration, start at one of the bounds + # of the domain. + if number_of_refinement == 0: + + # If x is greater than or equal 0.5, start the + # approximation at the upper bound of the domain. + if x >= 0.5: + approx_x, result = approximate_inverse_value(self.beta) + + # If x is less than 0.5, start the approximation at + # the lower bound of the domain. + else: + approx_x, result = approximate_inverse_value( + self.alpha) + else: + + # If this is not the first iteration, halve the stepsize + # and call the approximation method. + stepsize /= 2 + approx_x, result = approximate_inverse_value(result) + + number_of_refinement += 1 + + return result + + def mean(self): + """ + Computes the mean value (expectation) of the distribution. + + Returns: + float: the mean value + """ + + # Use region_expectation to compute the mean value. + return self.region_expectation((self.alpha, self.beta)) + + @memoize_method + def region_expectation(self, region): + """ + Computes the mean value (expectation) of a specified region. + + Args: + region: the region (tuple of dimension 2) of which the expectation + is to be computed + + Returns: + float: the expectation + """ + + # Check whether region is a tuple of dimension 2. + if isinstance(region, tuple) and len(region) == 2: + a, b = region + if a > b: + raise ValueError("Error: The upper bound of 'region' can't be " + "less than the lower bound.") + else: + raise TypeError("Error: Parameter 'region' must be a tuple of " + "dimension 2.") + + integral, _ = scipy.integrate.quad(lambda x: x * self.pdf(x), a, b) + + return integral + + @memoize_method + def region_probability(self, region): + """ + Computes the probability of a specified region. + + Args: + region: the region of which the probability is to be computed + + Returns: + float: the probability + """ + + # Compute the region's probability by integration, + + # Check whether region is a tuple of dimension 2. + if isinstance(region, tuple) and len(region) == 2: + a, b = region + integral, _ = scipy.integrate.quad(self.pdf, a, b) + else: + raise ValueError("Error: Parameter 'region' must be a tuple of" + " dimension 2.") + + return integral + + def conditional_expectation(self, interval, cdf_inverse_tolerance=1e-4, + cdf_inverse_max_refinements=10, + cdf_tolerance=1e-4): + """ + This computes the conditional expectation of the distribution + conditioned on being in the hyperrectangle passed in. + The hyperrectangle will actually for this be just an interval contained + in [0, 1] potentially with some cutouts. This will work for + 1-dimensional hyperrectangles, the multivariate distribution subclass + should implement a different version of this. + + If the region is (a, b), this will compute the expectation on + [cdf^-1(a), cdf^-1(b)] and divide it by (b-a). + + Args: + Interval (Interval): An interval on which + the conditional expectation is to be computed on + cdf_inverse_tolerance (float): The accuracy which the inverse + cdf is to be calculated to + cdf_inverse_max_refinements (int): The number of times the + the partition on the x-domain will be made finer + cdf_tolerance (float): The accuracy to which the cdf is calculated + to + """ + a, b = interval.a, interval.b + cdf_inverse = lambda x: self.cdf_inverse(x, cdf_inverse_tolerance, + cdf_inverse_max_refinements, + cdf_tolerance) + + lower, upper = cdf_inverse(a), cdf_inverse(b) + expectation = self.region_expectation((lower, upper)) + probability = b-a + + # A hyperrectangle may subtract some intervals from the larger interval + if hasattr(interval, 'cutouts'): + for cutout in interval.cutouts: + a, b = cutout.a, cutout.b + lower, upper = cdf_inverse(a), cdf_inverse(b) + expectation -= self.region_expectation((lower, upper)) + probability -= b-a + return expectation / probability + + def sample_one(self): + """ + Returns a single sample of the distribution + + Returns: + float: the sample + """ + + return self.cdf_inverse(np.random.uniform()) + + def sample_on_interval(self, a, b): + """ + This samples from the distribution conditioned on X being in [a, b]. + This does this by sampling uniformly on [F(a), F(b)] and then applying + the inverse transform to the result. + + Args: + a (float): The lower limit of the interval + b (float): The upper limit of the interval + Returns: + float: The sampled value in the interval + """ + return self.sample_between_quantiles(self.cdf(a), self.cdf(b)) + + def sample_between_quantiles(self, a, b): + """ + This samples from the distribution conditioned on the quantile of the + point being between a and b, i.e., it generates X given that + a < F(X) < b. It does this by sampling from a uniform distribution on + (a,b) and then applying the inverse transform to the point. + + Args: + a (float): The lower quantile, must be between 0 and 1. + b (float): The upper quantile, must be between 0 and 1. + Returns: + float: The sampled value + """ + y = np.random.uniform(a, b) + return self.cdf_inverse(y) + + def log_likelihood(self, data): + """ + This method will return the log likelihood of the observed data + given the fitted model. + + Args: + data (list[float]): A list of observed values + Returns: + float: The computed log-likelihood + """ + return sum(np.log(self.pdf(x)) for x in data) + diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distribution_factory.py b/mpisppy/confidence_intervals/bootsp/statdist/distribution_factory.py new file mode 100644 index 000000000..2377db479 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/distribution_factory.py @@ -0,0 +1,105 @@ +############################################################################### +# 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. +############################################################################### +""" +distribution_factory.py + +This module will export a distribution_factory function which should essentially +accept a name for a distribution and return the class associated with that distribution. +This will work by performing a scan through all the modules in this directory and finding +all the classes that are "registered" as distributions. + +Registering a entails using the class decorator register which is also +exported from this module. +""" + + +distribution_registry = {} + + +# Right now, this function simply scans through all modules in the current working directory +# It is perhaps more desirable (and safer) to have a list of modules to scan through +# This would be even easier to implement + +def import_all_classes(): + """ + Imports all classes in the current directory and stores + all registered distribution classes in the distribution_registry object + """ + if distribution_registry: + # already populated; the registry is stable, so scan only once + return + from . import base_distribution + from . import distributions + for mod in (base_distribution, distributions): + for var in mod.__dict__: + obj = getattr(mod, var) + if (hasattr(obj, "is_registered_distribution") and + getattr(obj, "is_registered_distribution")): + + distribution_registry[obj.registered_name] = obj + + +def register_distribution(name, ndim=None): + """ + Class decorator with arguments to register a class for use in distribution_factory + One must specify the name of the distribution to register under and the number of + dimensions the input argument must be (if there is a requirement). + + This will add the attributes is_registered_distribution, registered_name, and registered_ndim to the class + + Names will be case insensitive. + + Examples: + To register a function, simply use this as a decorator before any + class to be registered + + @register(name='clayton', ndim=2) + class ClaytonCopula(CopulaBase): + ... + + Args: + name (str): The name the distribution will be registered under + ndim (:obj: `int`, optional): The required dimensionality of input. + Defaults to None signaling no requirement + + Returns: + A class decorator to register a class as a distribution + """ + + def class_decorator(cls): + cls.is_registered_distribution = True + cls.registered_name = name.lower() + cls.registered_ndim = ndim + return cls + + return class_decorator + + +def distribution_factory(name): + """ + This function will accept a name and return the associated + distribution class raising an error if the name is unrecognized. + Names should be case insensitive + + Args: + name (str): The name of the distribution wanted + Returns: + The distribution class associated with name + """ + + import_all_classes() + + try: + distribution = distribution_registry[name.lower()] + except KeyError: + possible_names = '\n'.join(sorted(distribution_registry.keys(), key=str.lower)) + raise NameError("The specified distribution {} is unrecognized. " + "Possible distributions are:\n{}".format(name, possible_names)) + + return distribution diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py new file mode 100644 index 000000000..cdc0ff7c0 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -0,0 +1,902 @@ +############################################################################### +# 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. +############################################################################### +""" +distributions.py + +This module houses a host of distribution classes which all adhere to the +interface defined in base_distribution.py +""" + +import math +from collections import OrderedDict + +import numpy as np +# scipy is an optional dependency; import it lazily so the empirical +# bootstrap path stays scipy-free (mmw_ci.py uses the same pattern). +from pyomo.common.dependencies import scipy + +from mpisppy.confidence_intervals.bootsp.statdist.distribution_factory import register_distribution +from mpisppy.confidence_intervals.bootsp.statdist.base_distribution import Parameter +from mpisppy.confidence_intervals.bootsp.statdist.base_distribution import UnivariateDistribution +from mpisppy.confidence_intervals.bootsp.statdist.utilities import memoize_method +from mpisppy.confidence_intervals.bootsp.statdist import splines + +epsilon = 1e-12 + +@register_distribution(name="univariate-unif",ndim=1) +class UnivariateUniformDistribution(UnivariateDistribution): + """ + This class creates a univariate uniform distribution in the segment [a,b] + + Attributes: + a (float): The lower bound of the support of the distribution + b (float): The upper bound of the support of the distribution + """ + + def __init__(self, a, b): + """ + To construct a UnivariateUniformDistribution object, one must pass in + the lower and upper bounds for the support of the distribution. + These are passed in through a and b. + + Args: + a (float): The lower bound of the support of the distribution + b (float): The upper bound of the support of the distribution + """ + if a==b: + raise ValueError("The bounds should be different") + self.a=a + self.b=b + params = [Parameter('a', a), Parameter('b', b)] + UnivariateDistribution.__init__(self, params) + + @classmethod + def fit(cls, data): + """ + This method will fit a uniform distribution to the data. This will + set the lower bound of the distribution to the minimum of the data + and the upper bound to the maximum. + + Args: + data (List[float]): The list of values to fit the data to + Returns: + UnivariateUniformDistribution: The fitted uniform distribution + """ + return UnivariateUniformDistribution(min(data), max(data)) + + def pdf(self, x): + """ + Args: + x (float): The values where you want to compute the pdf + + Returns: + (float) The value of the probability density function of this + distribution on x. + """ + if xself.b: + return 0 + else: + return 1/(self.b-self.a) + + def cdf(self, x): + """ + Args: + x (float): The values where you want to compute the cdf + + Returns: + (float) The value of the cumulative density function + """ + if x 4 has excess kurtosis + # 6/(df-4), so method of moments gives df = 4 + 6/excess_kurtosis. Data + # with excess kurtosis at or below zero (tails no heavier than the normal) + # cannot be matched by any finite-variance t, so df falls back to this + # large value, which is numerically indistinguishable from the normal. + _FIT_DF_MAX = 1.0e6 + _FIT_KURT_TOL = 1.0e-12 + + @classmethod + def fit(cls, data): + """ + Fit a student's t distribution to the passed-in data. + + The mean and variance are taken to be the sample mean and variance. + The degrees of freedom are estimated by method of moments on the + excess kurtosis: a t with df > 4 has excess kurtosis 6/(df-4), so + df = 4 + 6/excess_kurtosis. The constructor then derives the scale + from (variance, df), so this rule is scale-free -- unlike the old + 2v/(v-1), it does not change when the data are rescaled, and any + variance v > 0 can be fit (the old v <= 1 restriction is gone). + + A sample whose excess kurtosis is zero or negative (tails no heavier + than the normal) has no finite-variance t that matches it, so df + falls back to a large value (effectively the normal). Because sample + excess kurtosis is always finite, this rule always yields df > 4 and + so cannot reach the heavy-tailed 2 < df <= 4 regime; construct such a + distribution directly if it is needed. + + Args: + data (List[float]): The list of values to fit the data to + Returns: + UnivariateStudentDistribution: The fitted student's t distribution + """ + mean = np.mean(data) + var = np.var(data) + excess_kurt = scipy.stats.kurtosis(data, fisher=True) + if not np.isfinite(excess_kurt) or excess_kurt <= cls._FIT_KURT_TOL: + df = cls._FIT_DF_MAX + else: + df = min(4.0 + 6.0 / excess_kurt, cls._FIT_DF_MAX) + + return UnivariateStudentDistribution(df, mean, var) + + def pdf(self, x): + """ + The probability distribution function of the student's t distribution. + + Args: + x (float): The values where you want to compute the pdf + + Returns: + float: The value of the probability density function of this + distribution on x. + """ + return self.distribution.pdf(x) + + def cdf(self, x): + """ + The cumulative distribution function for a student's t distribution. + + Args: + x (float): The values where you want to compute the cdf + + Returns: + float: The value of the cumulative density function of this + distribution on x. + """ + return self.distribution.cdf(x) + + def cdf_inverse(self, x): + """ + The inverse cumulative distribution function for a student's t + distribution. + + Args: + x (float): The values where you want to compute the inverse cdf + + Returns: + float: The value of the cumulative density function of this + distribution on x. + """ + return self.distribution.ppf(x) + + def generates_X(self, n=1): + return self.distribution.rvs(n) + + +@register_distribution(name="univariate-kernel", ndim=1) +class UnivariateGaussianKernelDistribution(UnivariateDistribution): + """ + This class will fit an Gaussian kernel density estimtion to a vector of data. + """ + def __init__(self, input_data, bw_method=None, dom_std=1): + """ + Initializes the distribution. + + Args: + input_data: list of data points + """ + + # Check the type of the input data and sort it. + self.kernel = scipy.stats.gaussian_kde(input_data, bw_method) + self.alpha = min(input_data) - dom_std*np.std(input_data) + self.beta = max(input_data) +dom_std*np.std(input_data) + UnivariateDistribution.__init__(self) + + + @classmethod + def fit(cls, data, bw_method=None): + """ + This function will fit an empirical distribution to the data. + + Args: + data (List[float]): The data to fit the distribution to + Returns: + UnivariateEmpiricalDistribution: The fitted distribution + """ + return UnivariateGaussianKernelDistribution(data, bw_method=bw_method) + + def pdf(self, x): + """ + Args: + x (float): The values where you want to compute the pdf + + Returns: + (float) The value of the probability density function of this + distribution on x. + """ + # gaussian_kde.evaluate returns an array; the base-class cdf feeds this + # to scipy.integrate.quad, which needs a scalar integrand (a size-1 + # array triggers a NumPy>=1.25 "array to scalar" DeprecationWarning). + return float(self.kernel.evaluate(x)[0]) + + def _cdf(self, x): + """ + Args: + x (float): The values where you want to compute the cdf + + Returns: + (float) The value of the cumulative density function + """ + return self.kernel.integrate_box_1d(self.lower, x) + + #use the cdf_inverse function in base_distribution + + def generates_X(self, n=1, seed=None): + # gaussian_kde.resample takes the sample size and the seed as two + # separate arguments + return self.kernel.resample(n, seed) + + + + +@register_distribution(name="univariate-epispline", ndim=1) +class UnivariateEpiSplineDistribution(UnivariateDistribution): + def __init__(self, input_data, error_distribution_domain='4,min,max', + specific_prob_constraint=None, seg_N=20, seg_kappa=100, + non_negativity_constraint_distributions=0, + probability_constraint_of_distributions=1, + nonlinear_solver='ipopt'): + """ + Initializes the distribution. + + Args: + input_data: list, dict or OrderedDict of data + dom: A number (int or float) specifying how many standard + deviations we want to consider as a domain of the distribution + or a string that defines the sign of the domain (pos for + positive and neg for negative). + specific_prob_constraint: either a tuple or a list of length 2 + with values for alpha and beta + seg_N (int): An integer specifying the number of knots + seg_kappa (float): A bound on the curvature of the spline + non_negativity_constraint_distributions: Set to 1 if u and w should + be nonnegative + probability_constraint_of_distributions: Set to 1 if integral of + the distribution should sum to 1 + nonlinear_solver (str): String specifying which solver to use + """ + + self.dom = error_distribution_domain + + self.specific_prob_constraint = specific_prob_constraint + + self.seg_N = seg_N + self.seg_kappa = seg_kappa + self.non_negativity_constraint_distributions = \ + non_negativity_constraint_distributions + self.probability_constraint_of_distributions = \ + probability_constraint_of_distributions + self.nonlinear_solver = nonlinear_solver + + # Fits the epi-spline distribution and computes the lower and upper + # bounds of the domain of the distribution + model, self.alpha, self.beta = splines.fit_distribution( + input_data, error_distribution_domain, specific_prob_constraint, + seg_N, seg_kappa, non_negativity_constraint_distributions, + probability_constraint_of_distributions, nonlinear_solver) + + # We only store the relevant parameters of the model and not the model + # itself, as the model takes a significant amount of memory. + self.tau = {i: model.tau[i] for i in model.tau} + self.a = {i: model.a[i].value for i in model.a} + self.delta = model.delta.value + self.w0 = model.w0.value + self.u0 = model.u0.value + + # Construction of Parameter objects + params = [] + params.extend([Parameter('tau_{}'.format(i), self.tau[i]) for i in self.tau]) + params.extend([Parameter('a_{}'.format(i), self.a[i]) for i in self.a]) + params.append(Parameter('delta', self.delta)) + params.append(Parameter('w0', self.w0)) + params.append(Parameter('u0', self.u0)) + + # This is an approximation of the integral of the normalized pdf + # over the entire domain + self.area = self.beta - self.alpha + + UnivariateDistribution.__init__(self, params, self.alpha, self.beta) + + @classmethod + def fit(cls, data, **distr_options): + """ + This function will fit a Univariate EpiSpline Distribution to the data + passed in. See the docstring for UnivariateEpiSplineDistribution + for a description of all the options that can be passed in. + + Args: + data (List[float]): A list of the points to fit the distribution to + Returns: + UnivariateEpiSplineDistribution: The fitted distribution + """ + return UnivariateEpiSplineDistribution(data, **distr_options) + + def _cdf_inverse(self, y): + """ + prescient_1.0-style inverse calculation + + Integrates the normalized pdf to get cdfs. + Then interpolates this cdf to get approximation + of cdf_inverse. The inverse is in the interval [0,1] + since it acts on the normalized pdf + Args: + y (float): number in [0,1] to calculate inverse of, + or list of values + Returns: + Either a list of computed inverses or a single inverse depending + on what has been passed in + """ + xs = np.linspace(0,1,100) + cdfs = [scipy.integrate.quad(self._normalized_pdf, 0, x)[0] for x in xs] + if isinstance(y, float): + return np.interp([y], cdfs, xs) + else: + return np.interp(y, cdfs, xs) + + def _normalized_pdf(self, x): + """ + Evaluates the pdf as if it were defined on [0,1] + Args: + x (float): a number in [0,1] + Returns: + float: the pdf of the unnormalized data + """ + k = int(math.ceil(self.seg_N * x)) + if k > 1: + tauk = self.tau[k - 1] + else: + tauk = 0 + k = 1 # avoids errors when i = 0 + summation = sum((x - self.tau[j] + 0.5 * self.delta) + * self.a[j] for j in range(1, k)) + + w = (self.w0 + self.u0 * x + + self.delta * summation + + 0.5 * self.a[k] * (x - tauk) ** 2) + + return math.exp(-w) + + @memoize_method + def _normalized_cdf(self, x): + return scipy.integrate.quad(self._normalized_pdf, 0, x)[0] + + def pdf(self, x): + """ + Evaluates the probability density function at a given point x. + + Args: + x (float): the point at which the pdf is to be evaluated + + Returns: + float: the value of the pdf + + Note: + The pdf values are calculated for the original scale of the data + (i.e. not normalized to [0,1]). + """ + # Set the pdf to 0 if the variable is out of bounds. + if x > self.beta or x < self.alpha: + return 0 + + # Noramlize x for using the normalized model. + norm_x = (x - self.alpha)/(self.beta - self.alpha) + + # Scale the return value to the original data. + return self._normalized_pdf(norm_x) / self.area + + +@register_distribution(name="univariate-empirical", ndim=1) +class UnivariateEmpiricalDistribution(UnivariateDistribution): + """ + This class will fit an empirical distribution to a vector of data. + """ + def __init__(self, input_data): + """ + Initializes the distribution. + + Args: + input_data: list of data points + """ + + # Check the type of the input data and sort it. + input_data = sorted(input_data) + + if len(input_data) == 0: + raise ValueError("You must provide at least one value to fit an " + "empirical distribution to the data.") + + self.alpha = input_data[0] + self.beta = input_data[len(input_data)-1] + self.input_data = input_data + UnivariateDistribution.__init__(self) + + @classmethod + def fit(cls, data): + """ + This function will fit an empirical distribution to the data. + + Args: + data (List[float]): The data to fit the distribution to + Returns: + UnivariateEmpiricalDistribution: The fitted distribution + """ + return UnivariateEmpiricalDistribution(data) + + def pdf(self, x): + """ + Evaluates the discrete probability of a given point x. + + Args: + x (float): the point at which the probability is to be evaluated + + Returns: + float: the probability + """ + + # Count all self.input_data that are equal to x. + number = sum(1 for y in self.input_data if x == y) + + return number/len(self.input_data) + + def cdf(self, x, lower_bound=None, upper_bound=None): + """ + This method calculates a empirical cdf, which is fitted to the data by + interpolation. If a lower bound is provided, any point smaller will + have cdf value 0. If an upper bound is provided, any point larger will + have cdf value 1. If either is not provided the value is estimated + using the line between the nearest two self.input_data. + + Args: + x (float): the point at which the cdf is to be evaluated + lower_bound (float): the lower bound + upper_bound (float): the upper bound + + Returns: + float: the value of the cdf + + Notes: + This method was copied from PINT's distributions class. + """ + + n = len(self.input_data) + lower_neighbor = None + lower_neighbor_index = None + upper_neighbor = None + upper_neighbor_index = None + for index in range(n): + if self.input_data[index] <= x: + lower_neighbor = self.input_data[index] + lower_neighbor_index = index + if self.input_data[index] > x: + upper_neighbor = self.input_data[index] + upper_neighbor_index = index + break + + if lower_neighbor == x: + cdf_x = (lower_neighbor_index + 1) / (n + 1) + + elif lower_neighbor is None: # x is smaller than all of the values + if lower_bound is None: + x1 = self.input_data[0] + index1 = self._count_less_than_or_equal(self.input_data, x1) + + if index1 == n: + # every record has the same value, so there is no second + # point to extrapolate along; x is below the lone mass point + return 0 + + x2 = self.input_data[index1] + index2 = self._count_less_than_or_equal(self.input_data, x2) + + y1 = index1 / (n + 1) + y2 = index2 / (n + 1) + interpolating_line = interpolate_line(x1, y1, x2, y2) + cdf_x = max(0, interpolating_line(x)) + else: + if lower_bound > x: + cdf_x = 0 + else: + x1 = lower_bound + x2 = upper_neighbor + y1 = 0 + y2 = 1 / (n + 1) + interpolating_line = interpolate_line(x1, y1, x2, y2) + cdf_x = interpolating_line(x) + + elif upper_neighbor is None: # x is greater than all of the values + if upper_bound is None: + if self.input_data[0] == self.input_data[n - 1]: + # every record has the same value, so there is no second + # point to extrapolate along; x is above the lone mass point + return 1 + j = n - 1 + while self.input_data[j] == self.input_data[n - 1]: + j -= 1 + x1 = self.input_data[j] + x2 = self.input_data[n - 1] + y1 = (j+1) / (n + 1) + y2 = n / (n + 1) + interpolating_line = interpolate_line(x1, y1, x2, y2) + cdf_x = min(1, interpolating_line(x)) + else: + if upper_bound < x: + cdf_x = 1 + else: + x1 = lower_neighbor + x2 = upper_bound + y1 = n / (n + 1) + y2 = 1 + interpolating_line = interpolate_line(x1, y1, x2, y2) + cdf_x = interpolating_line(x) + else: + x1 = lower_neighbor + x2 = upper_neighbor + y1 = (lower_neighbor_index + 1) / (n + 1) + y2 = (upper_neighbor_index + 1) / (n + 1) + interpolating_line = interpolate_line(x1, y1, x2, y2) + cdf_x = interpolating_line(x) + + return cdf_x + + def cdf_inverse(self, x, lower_bound=None, upper_bound=None): + """ + This method calculates a empirical inverse cdf, which is fitted to the + data by interpolation. + + Args: + x (float): the point at which the inverse cdf is to be evaluated + lower_bound (float): the lower bound + upper_bound (float): the upper bound + + Returns: + float: the value of the inverse cdf + + Notes: + This method was copied from PINT's distributions class. + """ + + n = len(self.input_data) + if x < 0 or x > 1: + raise ValueError('x must be between 0 and 1!') + if self.input_data[0] == self.input_data[n - 1]: + # every record has the same value (this includes a one-record + # sample): that value is every quantile, and neither extrapolation + # below nor above has a second point to find a slope from + return self.input_data[0] + # compute 'index' of this x + index = x * (n + 1) - 1 + first_index = self._count_less_than_or_equal( + self.input_data, self.input_data[0]) - 1 + + if index < first_index: + if lower_bound is None: + # take linear function through (0, self.input_data[0]) and + # (1, self.input_data[1]) + # input_data[0]) could occur several times, + # so find highest index j with input_data[j] = input_data[0] + first_index += 1 + second_index = self._count_less_than_or_equal( + self.input_data, self.input_data[first_index]) + interpolating_line = interpolate_line( + first_index / (n + 1), self.input_data[0], + second_index / (n + 1), self.input_data[first_index]) + + return interpolating_line(x) + else: + return lower_bound * (1 / (n + 1) - x) / (1 / (n + 1)) + \ + self.input_data[0] * x / (1 / (n + 1)) + elif index > n - 1: + if upper_bound is None: + # take linear function through (n-2, input_data[n-2]) and + # (n-1, input_data[n-1]) + # NOTE: input_data[n-1] could occur several times, + # so find lowest index j with input_data[j] = input_data[n-1] + # the all-equal case returned above, so this walk down the run + # of largest values always stops with j >= 1 + j = n - 1 + while self.input_data[j] == self.input_data[j - 1]: + j -= 1 + # g(x) = a*x + b + a = self.input_data[j] - self.input_data[j - 1] + b = self.input_data[j - 1] - (self.input_data[j] + - self.input_data[j-1]) * (j-1) + return a * index + b + else: + return self.input_data[n - 1] * \ + (1 - x) / (1 - n / (n + 1)) + \ + upper_bound * (x - n / (n + 1)) / (1 - n / (n + 1)) + else: + if math.floor(index) == index: + return self.input_data[math.floor(index)] + else: + interpolating_line = interpolate_line( + x1=math.floor(index), + y1=self.input_data[math.floor(index)], + x2=math.ceil(index), y2=self.input_data[math.ceil(index)]) + return interpolating_line(index) + + def _count_less_than_or_equal(self, xs, x): + """ + Counts the number of elements less than or equal to x in + a sorted list xs + + Args: + xs: A sorted list of elements + x: An element that you wish to find the number of elements less + than it + + Returns: + int: The number of elements in xs less than or equal to x + """ + count = 0 + for elem in xs: + if elem <= x: + count += 1 + else: + break + return count + + +def interpolate_line(x1, y1, x2, y2): + """ + This functions accepts two points (passed in as four arguments) + and returns the function of the line which passes through the points. + + Args: + x1 (float): x-value of point 1 + y1 (float): y-value of point 1 + x2 (float): x-value of point 2 + y2 (float): y-value of point 2 + + Returns: + callable: the function of the line + """ + + if x1 == x2: + raise ValueError("x1 and x2 must be different values") + + def f(x): + slope = (y2 - y1) / (x2 - x1) + return slope * (x - x1) + y1 + + return f + +#========= +@register_distribution(name="univariate-discrete", ndim=1) +class UnivariateDiscrete(UnivariateDistribution): + """ + This class creates a discrete univariate distribution. + The constructor takes an ordered dict of breakpoints. + """ + + def __init__(self, breakpoints): + """ + Univariate Discrete distribution constructor. + args: + breakpoints (OrderedDict): [value] := probability, + which need to be in increasing value and with prob that sums to 1. + Written for 3.x+ + """ + if not isinstance(breakpoints, OrderedDict): + raise RuntimeError("DiscreteDistribution expecting breakpoints to be a dict") + + self.breakpoints = breakpoints + # check the breakpoints + tol = 1e-6 + sumprob = 0 + self.mean = 0 + Esqsum = 0 + lastval, prob = list(self.breakpoints.items())[0] + for val, prob in self.breakpoints.items(): + sumprob += prob + self.mean += prob * val + Esqsum += prob * val * val + if val < lastval: + raise RuntimeError("DiscreteDistribution dict must be ordered by val:"+str(val)+" < "+str(lastval)) + lastval = val + self.var = Esqsum - self.mean*self.mean + # the probabilities have to sum to one from *either* side: a set of + # breakpoints summing to (say) 0.5 is not a distribution + if abs(sumprob - 1) > tol: # could use gosm_options.cdf_tolerance + raise ValueError("Discrete distribution with total prob=" + +str(sumprob)+" tolerance="+str(tol)) + super(UnivariateDiscrete, self).__init__() + + def pdf(self, x): + raise RuntimeError("pdf called for a discrete distribution.") + + def cdf(self, x): + """ + Cummulative Distribution Function: prob(X < x), which is weird + Args: + x (float): The value where you want to compute the cdf + + Returns: + (float) The value of the cumulative density function of this distribution on x. + """ + lastval, prob = list(self.breakpoints.items())[0] + if x < lastval: + return 0 + elif x == lastval: + return prob + sumprob = 0 + for val, prob in self.breakpoints.items(): + sumprob += prob + if x == val: + return sumprob + if x > lastval and x < val: + return sumprob - prob + lastval = val + return sumprob # should be one if we got this far + + def cdf_inverse(self, x): + """ + Evaluates the inverse of the cdf at probability value x, but + that does not really fly for discrete distrs... + """ + raise RuntimeError("cdf_inverse called for a discrete distribution.") + + def sample_one(self): + """ + Returns a single sample from the distribution + + Returns: + float or int: the sample + """ + p = np.random.uniform() + sumprob = 0 + for val, prob in self.breakpoints.items(): + sumprob += prob + if sumprob >= p: + return val + # if the probs dont' quite sum to one... + val, prob = list(self.breakpoints.items())[-1] + return val + + def rect_prob(self,down,up): + """ + + Args: + up (float): the upper values where you want to compute the probability + down (float): the upper values where you want to compute the probability + + Returns: the probability of being between up and down + + """ + return (self.cdf(up)-self.cdf(down)) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/sampler.py b/mpisppy/confidence_intervals/bootsp/statdist/sampler.py new file mode 100644 index 000000000..18bdb4159 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/sampler.py @@ -0,0 +1,34 @@ +############################################################################### +# 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. +############################################################################### +# pseudo-random numbers from distributions + + +class Sampler: + """" + This class enables generation of pseudo random numbers from distributions + args: + distributions (list of BaseDistribution): we sample from inverse of the cdf; len implies sample dimension + stream (np.random): should be seeded and reseeded by the caller + """ + def __init__(self, distributions, stream): + self.distributions = distributions + self.stream = stream + + def sample_one(self): + """ + Return a single sample from the distribution as a list + """ + # independent variables + retval = [] + + for distr in self.distributions: + unorm = self.stream.uniform(0,1) + # print(f"{unorm=}") + retval.append(distr.cdf_inverse(unorm)) + return retval \ No newline at end of file diff --git a/mpisppy/confidence_intervals/bootsp/statdist/splines.py b/mpisppy/confidence_intervals/bootsp/statdist/splines.py new file mode 100644 index 000000000..92803f7af --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/splines.py @@ -0,0 +1,543 @@ +############################################################################### +# 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. +############################################################################### +""" +splines.py + +This module should house all of the functions related +to fitting and evaluating splines. +""" + +import math +from collections import OrderedDict + +import numpy as np +from pyomo.environ import * + +class Spline: + """ + This fits a epi-spline to the data passed in the lists x and y. + This has functions for evaluating the spline and computing the derivative + of the spline, evaluate and derivative, respectively. + + Args: + x (List[float]): A list of numbers + y (List[float]): A list of numbers where y = f(x) + positiveness_constraint (bool): Set to True if spline values should be + positive + increasingness_constraint (bool): Set to True if spline should be + increasing + seg_N (int): The desired number of knots for the spline + seg_kappa (float): The bound on the curvature of the spline + L1Linf_solver (str): The solver for the L1 norm minimizer + L2Norm_solver (str): The solver for the L2 norm minimizer + """ + def __init__(self, x, y, positiveness_constraint=False, + epifit_error_norm='L2', + seg_N=20, seg_kappa=100, L1Linf_solver='gurobi', + increasingness_constraint=False, L2Norm_solver='gurobi'): + self.model = fit_epispline(x, y, positiveness_constraint, + epifit_error_norm, + seg_N, seg_kappa, L1Linf_solver, + increasingness_constraint, L2Norm_solver) + self.alpha = self.model.alpha.value + self.beta = self.model.beta.value + self.delta = self.model.delta.value + + def _interval_index(self, x): + """ + Compute the index of the interval x is in in the spline. + + Args: + x (float): The value x + """ + l = int(math.ceil(float(x-self.alpha)/self.delta)) + if l == 0: + l = 1 + + return l + + def evaluate(self, x): + """ + Evaluates the spline at a point x + + Args: + x (float): The point to evaluate the spline at + """ + if x < self.alpha or x > self.beta: + raise ValueError("This spline is only defined on [{}, {}]".format( + self.alpha, self.beta)) + + m = self.model + + s0 = value(m.s0) + v0 = value(m.v0) + delta = value(m.delta) + + # We find what interval x is in + l = self._interval_index(x) + + return (s0 + v0 * x + delta * sum( + (x - j * delta + 0.5 * delta) + * value(m.a[j]) for j in range(1, l)) + + 0.5 * value(m.a[l]) * (x - (l - 1) * delta) ** 2) + + __call__ = evaluate + + def derivative(self, x): + """ + Evaluates the derivative of the spline at a point x + + Args: + x (float): The point to evaluate the derivative at + """ + l = self._interval_index(x) + m = self.model + + v0 = value(m.v0) + delta = value(m.delta) + + + return (v0 + delta*sum(value(m.a[j]) for j in range(1,l)) + + value(m.a[l])*(x-(l-1)*delta)) + + +def fit_epispline(x, y, positiveness_constraint=False, error_norm='L2', + seg_N=20, seg_kappa=100, L1Linf_solver='gurobi', + increasingness_constraint=False, L2Norm_solver='gurobi'): + """ + This functions fits an epispline to the function based on passed in input. + This approximates the function f(x) = y where x and y are passed in lists + of data. + + Args: + x (List[float]): A list of numbers + y (List[float]): A list of numbers where y = f(x) + positiveness_constraint (bool): Set to True if spline values should be + positive + increasingness_constraint (bool): Set to True if spline should be + increasing + seg_N (int): The desired number of knots for the spline + seg_kappa (float): The bound on the curvature of the spline + L1Linf_solver (str): The solver for the L1 norm minimizer + L2Norm_solver (str): The solver for the L2 norm minimizer + """ + if len(x) != len(y): + raise RuntimeError('***ERROR: x and y must have the same length.') + + # We first create a new model + model = ConcreteModel() + + # Sets + model.I = Set(initialize=list(range(len(x)))) + model.intervals = RangeSet(int(seg_N)) + + # Parameters + model.N = Param(initialize=int(seg_N)) + model.kappa = Param(initialize=float(seg_kappa)) + + model.alpha = Param(initialize=min(x)) + model.beta = Param(initialize=max(x)) + + def x_init(m, i): + return x[i] + + model.x = Param(model.I, initialize=x_init) + + def fx_init(m, i): + return y[i] + + model.fx = Param(model.I, initialize=fx_init) + + model.delta = Param(initialize=float(model.beta - model.alpha) / model.N) + + def k_init(m, i): + aux = int(math.ceil(float(m.x[i] - m.alpha.value) / m.delta)) + if aux == 0: + aux = 1 + return aux + + model.k = Param(model.I, initialize=k_init) + + # Variables + model.e = Var(model.I, within=Reals) + model.s = Var(model.I, within=Reals) + + model.s0 = Var(within=Reals, initialize=0.0) + model.v0 = Var(within=Reals, initialize=0.0) + model.a = Var(model.intervals, bounds=(-model.kappa, model.kappa), initialize=0.0) + + # Constraints + def compute_spline(m, i): + return m.s[i] == m.s0 + m.v0 * m.x[i] + m.delta * sum( + (m.x[i] - j * m.delta + 0.5 * m.delta) * m.a[j] for j in range(1, m.k[i])) \ + + 0.5 * m.a[m.k[i]] * (m.x[i] - (m.k[i] - 1) * m.delta) ** 2 + + model.ComputeSpline = Constraint(model.I, rule=compute_spline) + + # Positiveness + if positiveness_constraint is True: + eps = 0.01 + w = [float(i) for i in np.arange(min(x), max(x), eps)] + model.J = Set(initialize=list(range(len(w)))) + + def positive_spline(m, i): + l = int(math.ceil(float(w[i] - m.alpha) / m.delta)) + if l == 0: + l = 1 + return m.s0 + m.v0 * w[i] + m.delta * sum( + (w[i] - j * m.delta + 0.5 * m.delta) * m.a[j] for j in + range(1, l)) + 0.5 * m.a[l] * (w[i] - (l - 1) * m.delta) ** 2 >= 0 + + model.PositiveSpline = Constraint(model.J, rule=positive_spline) + + # Increasingness + if increasingness_constraint is True: + + # First derivative + def increasing_spline(m, i): + + l = int(math.ceil(float(w[i] - m.alpha) / m.delta)) + if l == 0: + l = 1 + return m.v0 + m.delta * sum(m.a[j] for j in + range(1, l)) + m.a[l] * (w[i] - 2 * (l - 1) * m.delta) >= 0 + + model.IncreasingSpline = Constraint(model.J, rule=increasing_spline) + + if error_norm == "L1": + def ePositiveSide_rule(m, i): + return m.e[i] >= m.fx[i] - m.s[i] + + model.eDefPos = Constraint(model.I, rule=ePositiveSide_rule) + + def eNegativeSide_rule(m, i): + return m.e[i] >= - m.fx[i] + m.s[i] + + model.eDefNeg = Constraint(model.I, rule=eNegativeSide_rule) + elif error_norm == "L2": + def compute_error_rule(m, i): + return m.e[i] == m.fx[i] - m.s[i] + + model.ComputeError = Constraint(model.I, rule=compute_error_rule) + else: + raise RuntimeError("***ERROR: Unknown error norm=" + error_norm + " selected") + + # Objective function + if error_norm == 'L1': + def Obj_rule(m): + return summation(m.e) + + model.Obj = Objective(rule=Obj_rule) + elif error_norm == 'L2': + def Obj_rule(m): + return sum(m.e[i] ** 2 for i in m.I) + + model.Obj = Objective(rule=Obj_rule, sense=minimize) + else: + raise RuntimeError("***ERROR: Unknown error norm=" + error_norm + " selected") + + # Instance creation and optimization + model.preprocess() + if error_norm == "L1": + opt = SolverFactory(L1Linf_solver) + opt.options.mip_tolerances_absmipgap = 0 + opt.options.mip_tolerances_mipgap = 0 + opt.options.mip_tolerances_integrality = 1e-9 + elif error_norm == 'L2': + opt = SolverFactory(L2Norm_solver) + else: + raise RuntimeError("***ERROR: Unknown error norm=" + error_norm + " selected") + + opt.solve(model, tee=False) + return model + + +def error_domain(e, dom=None): + """ + This computes the parameters alpha and beta + from a list or dictionary, errors, and a + string dom which specifies how to compute alpha and beta + + alpha and beta will act as the bound on the domain of error distribution. + If no domain is specified then these will be the minimum and maximum + of the data + + dom should be a string of the following form ",,...," + where is replaced by one of the following: + 1. A number specifying how many standard deviations away from mean + you want alpha and beta to be set to + 2. pos which fixes alpha to 0 if alpha was prior set to a negative value + 3. neg which sets beta to 0 if beta was prior set to a positive value + 4. min which sets alpha to min + 5. max which sets beta to max + These fields are processed in order and set alpha and beta to subsequent + values. This will set alpha (beta) to be the smallest (largest) + value found while processing each field. + + Args: + e (List[float]): A list of error values + dom (str): The specified error string + + Returns (alpha, beta) + """ + data = e + mu = np.mean(data) + sigma = np.std(data, ddof=1) + _min = min(data) + _max = max(data) + + pos_error = ('***Error: You set the domain to be positive and there are ' + + 'some data with negative values') + neg_error = ('***Error: You set the domain to be negative and there are ' + + 'some data with positive values') + + if dom is None: + return _min, _max + elif isinstance(dom, (int, float)): + return mu - dom*sigma, mu + dom*sigma + elif isinstance(dom, str): + # We set alpha (beta) to max (min) and decrease (increase) as we + # process each field to ensure we get the smallest (largest) value + # from all the fields + alpha, beta = _max, _min + fields = dom.split(',') + for i, field in enumerate(fields): + if is_number(field): + a = mu-float(field)*sigma + if a < alpha: + alpha = a + a = mu+float(field)*sigma + if a > beta: + beta = a + else: + if field == 'pos' and _min < 0: + raise RuntimeError(pos_error) + elif field == 'neg' and _max > 0: + raise RuntimeError(neg_error) + elif field == 'pos' and alpha < 0: + alpha = 0 + elif field == 'neg' and beta > 0: + beta = 0 + elif field == 'min' and alpha > _min: + alpha = _min + elif field == 'max' and beta < _max: + beta = _max + + return alpha, beta + else: + raise RuntimeError("Unrecognized data type for domain") + + +def is_number(n): + """ + This function checks if n can be coerced to a floating point. + + Args: + n (str): Possibly a number string + """ + try: + float(n) + return True + except: + return False + + +def fit_distribution(x, dom=None, specific_prob_constraint=None, + seg_N=20, seg_kappa=100, + non_negativity_constraint_distributions=0, + probability_constraint_of_distributions=1, + nonlinear_solver=None): + """ + Fits a univariate epi-spline distribution to the given data. + The additional parameter dom defines special characteristics of the support + of the distribution. It can be pos (positive domain), neg (negative domain) + or it can be also a float that defines how many standard deviations from + the mean define the support. + + Args: + x: list, dict or OrderedDict of data + dom: A number (int or float) specifying how many standard deviations we + want to consider as a domain of the distribution or a string that + defines the sign of the domain (pos for positive and neg + for negative). + specific_prob_constraint: either a tuple or a list of length 2 + with values for alpha and beta + seg_N (int): An integer specifying the number of knots + seg_kappa (float): A bound on the curvature of the spline + non_negativity_constraint_distributions: Set to 1 if u and w should be + nonnegative + probability_constraint_of_distributions: Set to 1 if integral should + sum to 1 + nonlinear_solver (str): String specifying which solver to use + Returns: + (AbstractModel, float, float): tuple consisting of an instance of the + model, alpha and beta + + Note: + The data in the model is normalized to [0,1]. + """ + + N = int(seg_N) + kappa = float(seg_kappa) + + # ------------------------------------------------------- + # Model construction + # ------------------------------------------------------- + + model = AbstractModel() + + # ------------------------------------------------------- + # Parameters + # ------------------------------------------------------- + + if isinstance(x, OrderedDict) or isinstance(x, dict): + days = list(x.keys()) + elif isinstance(x, list): + days = list(range(len(x))) + elif isinstance(x, np.ndarray): + days = list(range(len(x))) + else: + raise RuntimeError('***ERROR: Unknown type of input data.') + + intervals = list(range(1, N + 1)) + delta = float(1) / float(N) + + model.N = Param(within=PositiveReals, initialize=N) + model.delta = Param(within=PositiveReals, initialize=delta) + + if specific_prob_constraint is None: + alpha, beta = error_domain(x, dom) + else: + if isinstance(specific_prob_constraint, tuple): + alpha, beta = specific_prob_constraint + alpha = float(alpha) + beta = float(beta) # avoid numpy64 + elif isinstance(specific_prob_constraint, list): + if len(specific_prob_constraint) == 2: + alpha = specific_prob_constraint[0] + beta = specific_prob_constraint[1] + else: + raise RuntimeError('***ERROR: The list specific_prob_constraint has to have a length of 2.') + + elif isinstance(specific_prob_constraint, str): + alpha, beta = error_domain(x, dom) + else: + raise RuntimeError('***ERROR: specific_prob_constraint has either to be a tuple or a list of length 2.') + + if alpha == beta: # this means there is only a CONSTANT bias + return model, alpha, beta + + # Here we normalize the data. Then, m.et is in [0,1]. + def et_init(modelo, j, k=None): + if k != None: + val = float(x[j, k] - alpha) / (beta - alpha) + if val < 0.0: + return 0.0 + if val > 1.0: + return 1.0 + return val + else: + val = float(x[j] - alpha) / (beta - alpha) + if val < 0.0: + return 0.0 + if val > 1.0: + return 1.0 + return val + + model.et = Param(days, initialize=et_init) + + def tau_init(modelo, i): + return i * delta + + model.tau = Param(intervals, initialize=tau_init) + + # -------------------------------------------------------- + # Variables + # -------------------------------------------------------- + + if non_negativity_constraint_distributions == 1: + model.w0 = Var(within=NonNegativeReals) + model.u0 = Var(within=NonNegativeReals) + else: + model.w0 = Var() + model.u0 = Var() + model.a = Var(intervals, bounds=(0, kappa)) + + # -------------------------------------------------------- + # Constraints + # -------------------------------------------------------- + + if probability_constraint_of_distributions == 1: + def prob_rule(modelo): # The sum of the probabilities over all the domain must be 1. + if specific_prob_constraint is not None: + s = 0.01 + samp = numpy.arange(0.0, 1.0 + s, s) + aux = 0 + for x in samp: + x = float(x) + k = int(math.ceil(N * x)) + if k > 1: + tauk = modelo.tau[k - 1] + else: + tauk = 0 + k = 1 # avoids erros when i = 0 + aux += s * exp(-(modelo.w0 + modelo.u0 * x \ + + delta * sum( + (x - modelo.tau[j] + 0.5 * delta) * modelo.a[j] for j in range(1, k)) \ + + 0.5 * modelo.a[k] * (x - tauk) ** 2)) + else: + aux = delta * exp(-modelo.w0) + for i in intervals: + aux += delta * exp(-(modelo.w0 + modelo.u0 * modelo.tau[i] \ + + delta * sum( + (modelo.tau[i] - modelo.tau[j] + 0.5 * delta) * modelo.a[j] for j in range(1, i)) \ + + 0.5 * modelo.a[i] * delta ** 2)) + return aux == 1 + + model.prob = Constraint(rule=prob_rule) + + # ------------------------------------------------------- + # Objective function + # ------------------------------------------------------- + + def fobj_rule(modelo): # appending _rule we don't need to define rule=rulename + aux = 0 + for d in days: + k = int(math.ceil(N * modelo.et[d])) + if k > 1: + tauk = modelo.tau[k - 1] + else: + tauk = 0 + k = 1 # avoids erros when i = 0 + aux += modelo.w0 + modelo.u0 * modelo.et[d] \ + + delta * sum((modelo.et[d] - modelo.tau[j] + 0.5 * delta) * modelo.a[j] for j in range(1, k)) \ + + 0.5 * modelo.a[k] * (modelo.et[d] - tauk) ** 2 + aux /= len(days) + + # We add the integral + if probability_constraint_of_distributions != 1: + aux += delta * exp(-modelo.w0) + for i in intervals: + aux += delta * exp(-(modelo.w0 + modelo.u0 * modelo.tau[i] \ + + delta * sum( + (modelo.tau[i] - modelo.tau[j] + 0.5 * delta) * modelo.a[j] for j in range(1, i)) \ + + 0.5 * modelo.a[i] * delta ** 2)) + return aux + + model.fobj = Objective(rule=fobj_rule, sense=minimize) + + # ------------------------------------------------------- + # Instance creation and optimization + # ------------------------------------------------------- + instance = model.create_instance() + opt = SolverFactory(nonlinear_solver) + opt.solve(instance, tee=False) + + return instance, alpha, beta + diff --git a/mpisppy/confidence_intervals/bootsp/statdist/utilities.py b/mpisppy/confidence_intervals/bootsp/statdist/utilities.py new file mode 100644 index 000000000..a7c6d30ac --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/utilities.py @@ -0,0 +1,200 @@ +############################################################################### +# 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. +############################################################################### +""" +utilities.py + +This module will contain any miscellaneous utilities for processing data +or enhancing functions or anything else. + +This currently exports tools for memoizing functions and a context manager +which enables the use of changing the program level arguments. +""" + +import sys +import inspect +from functools import partial, wraps +from contextlib import contextmanager + +def normalize_args(func, pargs, kwargs): + """ + This function puts the arguments into a dictionary mapping + keywords to arguments. To do this it must look up the function spec + for positional arguments. + """ + + # This should be a list of the names of the arguments + spec = inspect.getargs(func.__code__).args + + # Convert pargs to a list temporarily if need to change any mutable + # types to immutable types + pargs = list(pargs) + # We normalize any list or dictionary arguments to tuples + for i, parg in enumerate(pargs): + if isinstance(parg, list): + pargs[i] = tuple(parg) + elif isinstance(parg, dict): + pargs[i] = tuple(sorted(parg.items())) + + for key, value in kwargs.items(): + if isinstance(value, list): + kwargs[key] = tuple(value) + elif isinstance(value, dict): + kwargs[key] = tuple(sorted(value.items())) + + return dict(list(kwargs.items()) + list(zip(spec, pargs))) + +def memoize(func): + """ + This function implements memoization of a function by internally + storing a dictionary which stores argument-return value pairs. This + is to be used as a function decorator. + + Note that this only works with functions which has hashable types as + arguments. This function is designed in particular + to work with functions which have referential transparency and thus, the + calculation of a function with the same arguments should be the same every + time. + + This will convert any list or dictionary arguments to tuples so that + they can be stored in a dictionary + + Warning: If this function is used over a long period of time with a variety + of arguments, it can use up a large amount of memory. Do not use this + with class methods as the cache will exist beyond the life of the + instance. + + Args: + func: The function to be memoized + """ + + results = {} + + @wraps(func) + def f(*pargs, **kwargs): + args = normalize_args(func, pargs, kwargs) + arg_key = tuple(sorted(args.items())) + if arg_key not in results: + results[arg_key] = func(*pargs, **kwargs) + return results[arg_key] + + return f + + +class memoize_method: + """ + This class will be used as a class method decorator to internally cache + the results of a method in an instance-level dictionary. This differs + from the function decorator memoize in that it will store any results + with the instance meaning that once the instance goes out of scope, the + cache will be garbage collected and this will not lead to memory leaks. + + Any objects passed to a memoized method should be hashable; a call with an + unhashable argument (a list or a dictionary, say) cannot be cached, so it + is simply passed through to the method every time. + + This will internally store in any object which has a method decorated + with this class a dictionary with the name _memoize_method__cache which + maps functions and their arguments to the corresponding values. + + Example Usage: + class Obj: + @memoize_method + def super_expensive_function(self, arg): + ... + + obj = Obj() + obj.super_expensive_function(1) # This time, it will be computed + obj.super_expensive_function(1) # This time, it will be faster + + This will only compute the function on the first call. On any + subsequent call, it will look it up in the instance cache. + """ + def __init__(self, func): + self.func = func + + def __get__(self, instance, cls): + """ + This method will turn the decorator into a descriptor. This means + that trying to access the memoized method will not return the normal + method, but a slightly modified method. + + In this case, if an instance is calling the method, it will return + the partially applied __call__ method to the instance. If a class + is calling the method, it will just return the method. + """ + if instance is None: + # This means we are calling it from the class directly + # We need to pass in all the arguments including instance + # This is not memoized + return self.func + else: + # Calling from the instance, just need arguments, not instance + # This will call __call__ and replace the first element of pargs + # with instance. + return partial(self, instance) + + def __call__(self, *pargs, **kwargs): + # The first argument to any instance method is always the instance + obj = pargs[0] + + # Because the attribute is __cache, the real attribute name is mangled + # to have the callable name first (in this case memoize_method) + if hasattr(obj, '_memoize_method__cache'): + cache = obj.__cache + else: + cache = obj.__cache = {} + + # the keyword *values* have to be part of the key: cdf(x, epsabs=1e-4) + # and cdf(x, epsabs=1e-9) are different questions + key = (self.func, pargs[1:], frozenset(kwargs.items())) + + try: + value = cache[key] + except KeyError: + value = cache[key] = self.func(*pargs, **kwargs) + except TypeError: + # an unhashable argument: call through rather than hiding the + # method's own behavior behind "unhashable type: 'list'" + value = self.func(*pargs, **kwargs) + return value + + +@contextmanager +def set_arguments(args): + """ + This function will act as a context manager and will set the sys.argv + variable to the list of arguments passed in. This will enable calling + other scripts from within python as if they were called from the command + line. + + Example: + Say you had a script which simply printed out the system arguments + defined like such in the file print_args.py: + import sys + def main(): + print(sys.argv) + + Then in a separate file, you could call this function and set the + system arguments to whatever you want for the entirety of the with + block and the arguments would be restored at the end. + + In call_print_args.py called like python call_print_args.py 1 2, + import print_args + if __name__ == '__main__': + print(sys.argv) # ['call_print_args.py', '1', '2'] + with set_arguments(['arg1', 'arg2', 'arg3']): + print_args.main() # ['arg1', 'arg2', 'arg3'] + print(sys.argv) # ['call_print_args.py', '1', '2'] + Args: + args (List[str]): A list of strings which will become the arguments + """ + sys.argv_ = sys.argv + sys.argv = args + yield + sys.argv = sys.argv_ diff --git a/mpisppy/confidence_intervals/bootsp/user_boot.py b/mpisppy/confidence_intervals/bootsp/user_boot.py index efcad3b5c..b3becf153 100644 --- a/mpisppy/confidence_intervals/bootsp/user_boot.py +++ b/mpisppy/confidence_intervals/bootsp/user_boot.py @@ -14,10 +14,45 @@ import mpisppy.confidence_intervals.ciutils as ciutils import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils import mpisppy.confidence_intervals.bootsp.boot_sp as boot_sp +import mpisppy.confidence_intervals.bootsp.smoothed_boot_sp as smoothed_boot_sp my_rank = boot_utils.my_rank +def _empirical_report(cfg, module, xhat): + """ Run and print an empirical bootstrap CI; return the 6-tuple. """ + ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap = \ + boot_sp.compute_ci(cfg, module, xhat) + + if my_rank == 0: + # print result + print(f"point estimator for optimal function value: {center_optimal}") + print(f"point estimator for function value at xhat: {center_upper}") + print(f"point estimator for optimality gap: {center_gap}") + ci_gap[0] = max(0, ci_gap[0]) + print(f"ci for optimal function value: {ci_optimal}") + print(f"ci for function value at xhat: {ci_upper}") + print(f"ci for optimality gap: {ci_gap}") + + return ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap + + +def _smoothed_report(cfg, module, xhat): + """ Run and print a smoothed bootstrap/bagging CI; return (ci_gap, center_gap). + + The smoothed methods estimate only the optimality-gap interval, so the + return signature differs from the empirical 6-tuple. + """ + result = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + if my_rank == 0: + ci_gap_two_sided, center_gap = result + ci_gap_two_sided[0] = max(0, ci_gap_two_sided[0]) + print(f"point estimator for the optimality gap: {center_gap}") + print(f"two-sided CI for optimality gap: {ci_gap_two_sided}") + return ci_gap_two_sided, center_gap + return result + + def main_routine(cfg, module): """ The top level of user_boot; called by __main__ and by test drivers. @@ -25,31 +60,21 @@ def main_routine(cfg, module): cfg (Config): parameters module (Python module): contains the scenario creator function and helpers Returns: + For an empirical boot_method, the 6-tuple (ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap); - the ci_* entries are None on MPI ranks other than 0. + for a smoothed boot_method, the pair (ci_gap_two_sided, center_gap). + The ci_* entries are None on MPI ranks other than 0. Note: - Prints the confidence-interval results to the terminal on rank 0. A - smoothed boot_method raises a friendly "not yet merged" error. + Prints the confidence-interval results to the terminal on rank 0. """ if cfg["xhat_fname"] is not None and cfg["xhat_fname"] != "None": xhat = ciutils.read_xhat(cfg["xhat_fname"]) else: xhat = boot_utils.compute_xhat(cfg, module) - ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap = \ - boot_sp.compute_ci(cfg, module, xhat) - - if my_rank == 0: - # print result - print(f"point estimator for optimal function value: {center_optimal}") - print(f"point estimator for function value at xhat: {center_upper}") - print(f"point estimator for optimality gap: {center_gap}") - ci_gap[0] = max(0, ci_gap[0]) - print(f"ci for optimal function value: {ci_optimal}") - print(f"ci for function value at xhat: {ci_upper}") - print(f"ci for optimality gap: {ci_gap}") - - return ci_optimal, ci_upper, ci_gap, center_optimal, center_upper, center_gap + if boot_utils.is_smoothed(cfg.boot_method): + return _smoothed_report(cfg, module, xhat) + return _empirical_report(cfg, module, xhat) if __name__ == '__main__': diff --git a/mpisppy/tests/test_boot_sp.py b/mpisppy/tests/test_boot_sp.py index bd6c3e1bb..a831af1d1 100644 --- a/mpisppy/tests/test_boot_sp.py +++ b/mpisppy/tests/test_boot_sp.py @@ -226,9 +226,47 @@ def test_compute_xhat_requires_generator(self): self.assertIn("xhat_generator", msg) self.assertIn("xhat_generator_no_generator_module", msg) - def test_smoothed_not_yet_merged_boot_sp(self): + def test_solve_routine_rejects_maximization(self): + # The estimators form the gap as (value at xhat) - (optimal), which is + # non-positive for a maximization, and the drivers floor the reported + # interval at 0 -- so a maximization run would silently report [0, 0]. + # It must raise instead (maximization either works or errors). + import pyomo.environ as pyo + import mpisppy.scenario_tree as scenario_tree + + def _make(sense): + def scenario_creator(scenario_name, **kwargs): + m = pyo.ConcreteModel() + m.x = pyo.Var(within=pyo.NonNegativeReals, bounds=(0, 1)) + m.obj_expr = pyo.Expression(expr=m.x) + m.obj = pyo.Objective(expr=m.obj_expr, sense=sense) + m._mpisppy_probability = "uniform" + m._mpisppy_node_list = [scenario_tree.ScenarioNode( + name="ROOT", cond_prob=1.0, stage=1, + cost_expression=m.obj_expr, nonant_list=[m.x], scen_model=m)] + return m + fake = types.ModuleType(f"sense_module_{sense}") + fake.scenario_creator = scenario_creator + fake.kw_creator = lambda cfg: {} + return fake + + cfg = _make_cfg() + cfg.solver_name = "nosolver" # must raise before any solve is attempted + with self.assertRaises(ValueError) as ctx: + boot_sp.solve_routine(cfg, _make(pyo.maximize), range(2)) + self.assertIn("minimization-only", str(ctx.exception)) + + # the same model as a minimization gets past the guard (and then fails + # on the bogus solver name, proving the guard was not what stopped it) + with self.assertRaises(Exception) as ctx: + boot_sp.solve_routine(cfg, _make(pyo.minimize), range(2)) + self.assertNotIn("minimization-only", str(ctx.exception)) + + def test_compute_ci_rejects_smoothed(self): + # compute_ci is the empirical dispatch; a smoothed method is routed to + # smoothed_boot_sp.compute_smoothed_ci instead and must be rejected here cfg = _make_cfg("Smoothed_boot_kernel") - with self.assertRaises(RuntimeError) as ctx: + with self.assertRaises(ValueError) as ctx: boot_sp.compute_ci(cfg, None, {"ROOT": [0.0, 5.0]}) self.assertIn("smoothed", str(ctx.exception).lower()) @@ -292,13 +330,6 @@ def test_user_boot_main_routine(self): else: self.assertEqual(res, (None, None, None, None, None, None)) - @unittest.skipIf(not solver_available, "no solver is available") - def test_user_boot_smoothed_raises(self): - module = boot_utils.module_name_to_module(MODULE_NAME) - cfg = _make_cfg("Smoothed_bagging") - with self.assertRaises(RuntimeError): - user_boot.main_routine(cfg, module) - #***************************************************************************** class Test_boot_sp_data(unittest.TestCase): diff --git a/mpisppy/tests/test_boot_sp_simulate.py b/mpisppy/tests/test_boot_sp_simulate.py index 54f2ec332..62cce73d4 100644 --- a/mpisppy/tests/test_boot_sp_simulate.py +++ b/mpisppy/tests/test_boot_sp_simulate.py @@ -135,13 +135,6 @@ def test_bagging_gatherv(self): else: self.assertEqual(res, (None, None, None, None, None, None)) - def test_smoothed_not_yet_merged(self): - # no solver needed; the guard fires before any solve - module = boot_utils.module_name_to_module(MODULE_NAME) - cfg = _make_cfg("Smoothed_boot_kernel") - with self.assertRaises(RuntimeError): - simulate_boot.main(cfg, module) - if __name__ == '__main__': unittest.main() diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py new file mode 100644 index 000000000..90862af19 --- /dev/null +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -0,0 +1,988 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +# Tests for the smoothed bootstrap/bagging code (bootsp) and the statdist +# univariate distributions, plus the empirical farmer/cvar examples that need +# statdist (and so could not live in test_boot_sp.py). Run serially: +# +# python -m pytest mpisppy/tests/test_boot_sp_smoothed.py +# Parallel (exercises the smoothed Gatherv batch split across ranks): +# mpiexec -np 2 python -m mpi4py mpisppy/tests/test_boot_sp_smoothed.py +# +# The smoothed methods fit a distribution with statdist (scipy); the kernel and +# bagging methods need only an LP/MIP solver, while the epi-spline methods also +# need a nonlinear solver (ipopt), so those tests are skipped when ipopt is +# absent. + +import os +import sys +import math +import importlib.util +import subprocess +import tempfile +import unittest +from collections import OrderedDict + +import numpy as np +from pyomo.common.dependencies import scipy +import pyomo.environ as pyo +import mpisppy.utils.sputils as sputils +from mpisppy.tests.utils import get_solver, round_pos_sig + +import mpisppy.confidence_intervals.bootsp.boot_utils as boot_utils +import mpisppy.confidence_intervals.bootsp.boot_sp as boot_sp +import mpisppy.confidence_intervals.bootsp.smoothed_boot_sp as smoothed_boot_sp +import mpisppy.confidence_intervals.bootsp.user_boot as user_boot +import mpisppy.confidence_intervals.bootsp.simulate_boot as simulate_boot +import mpisppy.confidence_intervals.bootsp.statdist.distributions as statdist_distributions +import mpisppy.confidence_intervals.bootsp.statdist.utilities as statdist_utilities +from mpisppy.confidence_intervals.bootsp.statdist.base_distribution import ( + Parameter, + UnivariateDistribution, +) +from mpisppy.confidence_intervals.bootsp.statdist.distribution_factory import ( + distribution_factory, +) + +sputils.disable_tictoc_output() + +solver_available, solver_name, persistent_available, persistent_solver_name = get_solver() +ipopt_available = pyo.SolverFactory("ipopt").available(exception_flag=False) +# matplotlib is the optional [plot] extra; only the plotting test needs it +matplotlib_available = importlib.util.find_spec("matplotlib") is not None + +comm = boot_utils.comm +n_proc = boot_utils.n_proc +my_rank = boot_utils.my_rank + +module_dir = os.path.dirname(os.path.abspath(__file__)) +bootsp_examples = os.path.join(module_dir, "..", "..", "examples", "bootsp") +for _sub in ("farmer", "cvar", "multi_knapsack"): + _d = os.path.join(bootsp_examples, _sub) + if not os.path.exists(_d): + raise RuntimeError(f"Directory not found: {_d}") + if _d not in sys.path: + sys.path.insert(0, _d) + +MK_DATA = os.path.abspath( + os.path.join(bootsp_examples, "multi_knapsack", "multi_knapsack_data.json")) + +univariate_tokens = ["univariate-unif", "univariate-normal", "univariate-student", + "univariate-kernel", "univariate-epispline", + "univariate-empirical", "univariate-discrete"] + + +def _make_cvar_cfg(method="Smoothed_bagging", seed=42, reps=2): + cfg = boot_utils._process_module("cvar") + cfg.module_name = "cvar" + cfg.max_count = 200 + cfg.candidate_sample_size = 5 + cfg.sample_size = 20 + cfg.subsample_size = 5 + cfg.nB = 8 + cfg.alpha = 0.1 + cfg.seed_offset = seed + cfg.xhat_fname = "None" + cfg.optimal_fname = "None" + cfg.trace_fname = None + cfg.coverage_replications = reps + cfg.solver_name = solver_name + cfg.boot_method = method + cfg.smoothed_B_I = 3 + cfg.smoothed_center_sample_size = 20 + return cfg + + +def _make_farmer_cfg(method="Classical_quantile", seed=100): + cfg = boot_utils._process_module("farmer") + cfg.module_name = "farmer" + cfg.max_count = 200 + cfg.candidate_sample_size = 5 + cfg.sample_size = 30 + cfg.subsample_size = 10 + cfg.nB = 8 + cfg.alpha = 0.1 + cfg.seed_offset = seed + cfg.xhat_fname = "None" + cfg.optimal_fname = "None" + cfg.trace_fname = None + cfg.coverage_replications = 2 + cfg.solver_name = solver_name + cfg.boot_method = method + cfg.crops_multiplier = 1 + cfg.yield_cv = 0.1 + return cfg + + +#***************************************************************************** +class Test_statdist(unittest.TestCase): + """ Direct tests of the trimmed statdist univariate distributions. """ + + def test_factory_resolves_univariate(self): + for token in univariate_tokens: + cls = distribution_factory(token) + self.assertTrue(hasattr(cls, "fit") or callable(cls), msg=token) + + def test_factory_rejects_unknown(self): + with self.assertRaises(NameError): + distribution_factory("not-a-distribution") + + def test_factory_drops_multivariate(self): + # the multivariate/copula distributions were trimmed out of the port + for token in ["multivariate-normal", "gaussian-copula"]: + with self.assertRaises(NameError): + distribution_factory(token) + + def test_registry_metadata(self): + # every univariate distribution registers under its lower-cased name + # and declares one dimension; the lookup itself is case insensitive + for token in univariate_tokens: + cls = distribution_factory(token) + self.assertEqual(cls.registered_name, token) + self.assertEqual(cls.registered_ndim, 1) + self.assertIs(distribution_factory(token.upper()), cls) + + # the subprocess imports mpi-sppy, so it initializes MPI; under an mpiexec + # launch it would inherit this job's environment and join it + @unittest.skipIf(n_proc > 1, "spawns a plain (non-MPI) python subprocess") + def test_scipy_not_imported_at_module_import(self): + # statdist defers scipy so the empirical path stays scipy-free; the + # distributions module must not pull scipy in merely on import, so ask + # a fresh interpreter (this one has scipy loaded by the tests below) + code = ("import sys;" + "import mpisppy.confidence_intervals.bootsp.statdist.distributions;" + "print('scipy loaded:', 'scipy' in sys.modules)") + done = subprocess.run([sys.executable, "-c", code], + capture_output=True, text=True) + self.assertEqual(done.returncode, 0, msg=done.stderr) + self.assertIn("scipy loaded: False", done.stdout) + + def test_uniform_inverse(self): + uunif = distribution_factory("univariate-unif")(0, 1) + mid = uunif.cdf_inverse(0.5) + self.assertAlmostEqual(mid, 0.5, places=6) + self.assertLessEqual(uunif.cdf_inverse(0.25), uunif.cdf_inverse(0.75)) + + def test_uniform_rejects_degenerate_support(self): + with self.assertRaises(ValueError): + distribution_factory("univariate-unif")(1.0, 1.0) + + def test_uniform_density_and_cdf(self): + unif = distribution_factory("univariate-unif")(2.0, 6.0) + self.assertEqual(unif.pdf(1.0), 0) # outside the support + self.assertEqual(unif.pdf(7.0), 0) + self.assertAlmostEqual(unif.pdf(4.0), 0.25) + self.assertEqual(unif.cdf(1.0), 0) + self.assertEqual(unif.cdf(2.0), 0) + self.assertAlmostEqual(unif.cdf(3.0), 0.25) + self.assertEqual(unif.cdf(6.0), 1) + self.assertEqual(unif.cdf(9.0), 1) + for q in (0.1, 0.5, 0.9): + self.assertAlmostEqual(unif.cdf(unif.cdf_inverse(q)), q) + + def test_uniform_fit_spans_the_data(self): + data = [3.0, -1.0, 2.5, 7.25] + unif = distribution_factory("univariate-unif").fit(data) + self.assertEqual((unif.a, unif.b), (min(data), max(data))) + self.assertEqual([p.value for p in unif.parameters], + [min(data), max(data)]) + + def test_uniform_generates_X_in_the_support(self): + unif = distribution_factory("univariate-unif")(2.0, 6.0) + unif.seed_reset(13) + draws = unif.generates_X(50) + self.assertEqual(len(draws), 50) + self.assertTrue(all(2.0 <= d <= 6.0 for d in draws)) + + def test_normal_inverse(self): + unorm = distribution_factory("univariate-normal")(mean=3.0, var=4.0) + self.assertAlmostEqual(unorm.cdf_inverse(0.5), 3.0, places=4) + self.assertLess(unorm.cdf_inverse(0.25), unorm.cdf_inverse(0.75)) + + def test_normal_fit_recovers_the_moments(self): + data = list(np.random.RandomState(3).normal(5.0, 2.0, size=500)) + norm = distribution_factory("univariate-normal").fit(data) + self.assertAlmostEqual(norm.mean, float(np.mean(data))) + self.assertAlmostEqual(norm.var, float(np.var(data))) + + def test_normal_matches_the_closed_form(self): + norm = distribution_factory("univariate-normal")(var=4.0, mean=3.0) + self.assertAlmostEqual(norm.cdf(3.0), 0.5) + self.assertAlmostEqual(norm.pdf(3.0), 1.0/math.sqrt(2*math.pi*4.0)) + self.assertAlmostEqual(norm.pdf(1.0), norm.pdf(5.0)) # symmetry + # the mass within one standard deviation of the mean + self.assertAlmostEqual(norm.cdf(5.0) - norm.cdf(1.0), 0.6826894921, + places=6) + self.assertAlmostEqual(norm.cdf(norm.cdf_inverse(0.975)), 0.975) + + def test_normal_generates_X(self): + norm = distribution_factory("univariate-normal")(var=1.0, mean=0.0) + norm.seed_reset(42) + draws = norm.generates_X(1000) + self.assertEqual(len(draws), 1000) + self.assertLess(abs(float(np.mean(draws))), 0.25) + + def test_student_fit_matches_the_data_moments(self): + # the fit promises the distribution's mean and variance are the data's; + # the constructor derives scipy's scale from (var, df), so the reported + # variance is var -- not var*df/(df-2), which is what passing sqrt(var) + # as the scale would have produced + data = list(np.random.RandomState(4).normal(0.0, 2.0, size=500)) + var = float(np.var(data)) + st = distribution_factory("univariate-student").fit(data) + self.assertAlmostEqual(st.distribution.mean(), float(np.mean(data))) + self.assertAlmostEqual(st.distribution.var(), var) + + def test_student_fit_sets_df_from_kurtosis(self): + # df is method of moments on the excess kurtosis: a t with df > 4 has + # excess kurtosis 6/(df-4), so df = 4 + 6/excess_kurtosis. Data drawn + # from a t with df=5 (excess kurtosis 6) should recover a df near 5. + data = list(scipy.stats.t(df=5).rvs(size=3000, random_state=7)) + ek = float(scipy.stats.kurtosis(data, fisher=True)) + self.assertGreater(ek, 0.0) # heavier than normal + st = distribution_factory("univariate-student").fit(data) + self.assertAlmostEqual(st.df, 4.0 + 6.0/ek) # the rule, exactly + self.assertGreater(st.df, 4.0) + self.assertLess(st.df, 10.0) # sane recovery of df=5 + self.assertAlmostEqual(st.distribution.var(), float(np.var(data))) + + def test_student_fit_falls_back_to_normal_for_light_tails(self): + # data with tails no heavier than the normal (here uniform, whose + # excess kurtosis is negative) has no finite-variance t that matches + # it, so df falls back to the large "effectively normal" value + data = list(np.random.RandomState(5).uniform(0.0, 1.0, size=1000)) + self.assertLess(float(scipy.stats.kurtosis(data, fisher=True)), 0.0) + st = distribution_factory("univariate-student").fit(data) + self.assertEqual(st.df, statdist_distributions. + UnivariateStudentDistribution._FIT_DF_MAX) + self.assertAlmostEqual(st.distribution.var(), float(np.var(data))) + + def test_student_variance_is_honored_for_any_df(self): + for df in (2.5, 4.0, 30.0): + st = distribution_factory("univariate-student")( + df=df, mean=-2.0, var=9.0) + self.assertAlmostEqual(st.distribution.var(), 9.0, msg=f"df={df}") + self.assertAlmostEqual(st.distribution.mean(), -2.0) + + def test_student_needs_a_df_that_has_a_variance(self): + for df in (1, 2.0): + with self.assertRaises(ValueError): + distribution_factory("univariate-student")( + df=df, mean=0.0, var=1.0) + + def test_student_fit_accepts_low_variance_data(self): + # because the scale is derived from (var, df), a variance of one or + # less is no longer special: the old 2v/(v-1) rule could not fit it, + # but the kurtosis rule can, and the fitted moments still match + data = list(np.random.RandomState(5).normal(0.0, 0.1, size=200)) + var = float(np.var(data)) + self.assertLessEqual(var, 1.0) + st = distribution_factory("univariate-student").fit(data) + self.assertGreater(st.df, 2.0) + self.assertAlmostEqual(st.distribution.mean(), float(np.mean(data))) + self.assertAlmostEqual(st.distribution.var(), var) + + def test_student_is_symmetric_with_heavier_tails(self): + st = distribution_factory("univariate-student")(df=3.0, mean=1.0, var=4.0) + self.assertAlmostEqual(st.cdf(1.0), 0.5) + self.assertAlmostEqual(st.pdf(0.0), st.pdf(2.0)) # symmetry + self.assertAlmostEqual(st.cdf(st.cdf_inverse(0.9)), 0.9) + # same mean and variance as a normal, but with the fatter tails + norm = distribution_factory("univariate-normal")(var=4.0, mean=1.0) + self.assertAlmostEqual(st.distribution.var(), norm.var) + self.assertLess(st.cdf_inverse(0.01), norm.cdf_inverse(0.01)) + self.assertGreater(st.cdf_inverse(0.99), norm.cdf_inverse(0.99)) + + def test_student_generates_X(self): + st = distribution_factory("univariate-student")(df=5.0, mean=0.0, var=1.0) + st.seed_reset(7) + self.assertEqual(len(st.generates_X(500)), 500) + + def test_kernel_fit_inverse(self): + # the kernel-density fit backs Smoothed_boot_kernel and Smoothed_bagging + data = list(np.random.RandomState(0).normal(0, 1, size=200)) + kde = distribution_factory("univariate-kernel").fit(data) + lo = kde.cdf_inverse(0.25) + hi = kde.cdf_inverse(0.75) + self.assertTrue(math.isfinite(lo) and math.isfinite(hi)) + self.assertLess(lo, hi) + + def test_kernel_pdf_returns_a_python_float(self): + # gaussian_kde.evaluate answers with an array, but the base-class cdf + # hands the density to scipy.integrate.quad, which wants a scalar + data = list(np.random.RandomState(6).normal(0, 1, size=100)) + kde = distribution_factory("univariate-kernel").fit(data) + self.assertIsInstance(kde.pdf(0.0), float) + + def test_kernel_honors_bw_method(self): + data = list(np.random.RandomState(7).normal(0, 1, size=100)) + cls = distribution_factory("univariate-kernel") + default = cls.fit(data) + wide = cls.fit(data, bw_method=2.0) + self.assertAlmostEqual(wide.kernel.factor, 2.0) + self.assertNotAlmostEqual(default.kernel.factor, wide.kernel.factor) + + def test_kernel_pads_the_domain(self): + data = [1.0, 2.0, 3.0, 4.0] + kde = distribution_factory("univariate-kernel")(data, dom_std=2) + sd = float(np.std(data)) + self.assertAlmostEqual(kde.alpha, 1.0 - 2*sd) + self.assertAlmostEqual(kde.beta, 4.0 + 2*sd) + + def test_kernel_cdf_is_a_distribution(self): + # the kernel class leans on the base-class cdf, which integrates the + # density numerically between alpha and beta + data = list(np.random.RandomState(8).normal(0, 1, size=60)) + kde = distribution_factory("univariate-kernel").fit(data) + self.assertEqual(kde.cdf(kde.alpha), 0) + self.assertEqual(kde.cdf(kde.beta), 1) + grid = [float(x) for x in np.linspace(kde.alpha, kde.beta, 12)] + values = [kde.cdf(x) for x in grid] + for lo, hi in zip(values, values[1:]): + self.assertLessEqual(lo, hi + 1e-6) + self.assertTrue(all(kde.pdf(x) >= 0 for x in grid)) + # the padded domain holds essentially all of the mass + self.assertGreater(kde.cdf(grid[-2]), 0.9) + + def test_kernel_generates_X_is_seeded(self): + data = list(np.random.RandomState(9).normal(0, 1, size=40)) + kde = distribution_factory("univariate-kernel").fit(data) + drawn = kde.generates_X(5, seed=11) + self.assertEqual(np.shape(drawn), (1, 5)) + np.testing.assert_allclose(drawn, kde.generates_X(5, seed=11)) + + def test_empirical_fit_inverse(self): + data = list(np.random.RandomState(1).normal(0, 1, size=200)) + emp = distribution_factory("univariate-empirical").fit(data) + self.assertLessEqual(emp.cdf_inverse(0.25), emp.cdf_inverse(0.75)) + + def test_empirical_rejects_empty_data(self): + with self.assertRaises(ValueError): + distribution_factory("univariate-empirical").fit([]) + + def test_empirical_uses_plotting_positions(self): + # the contract of the interpolated empirical cdf: the i-th smallest of + # n records sits at quantile (i+1)/(n+1), and cdf_inverse inverts that + data = [4.0, 1.0, 3.0, 2.0, 5.0] + emp = distribution_factory("univariate-empirical").fit(data) + n = len(data) + for i, value in enumerate(sorted(data)): + self.assertAlmostEqual(emp.cdf(value), (i + 1)/(n + 1)) + self.assertAlmostEqual(emp.cdf_inverse((i + 1)/(n + 1)), value) + # and it interpolates linearly between two records + self.assertAlmostEqual(emp.cdf(2.5), 2.5/(n + 1)) + + def test_empirical_cdf_is_monotone_and_bounded(self): + data = list(np.random.RandomState(10).normal(0, 1, size=30)) + emp = distribution_factory("univariate-empirical").fit(data) + values = [emp.cdf(float(x)) + for x in np.linspace(min(data) - 1, max(data) + 1, 40)] + self.assertTrue(all(0 <= v <= 1 for v in values)) + for lo, hi in zip(values, values[1:]): + self.assertLessEqual(lo, hi) + + def test_empirical_pdf_is_the_relative_frequency(self): + emp = distribution_factory("univariate-empirical").fit( + [1.0, 2.0, 2.0, 3.0]) + self.assertAlmostEqual(emp.pdf(2.0), 0.5) + self.assertAlmostEqual(emp.pdf(1.0), 0.25) + self.assertEqual(emp.pdf(9.0), 0) + + def test_empirical_respects_explicit_bounds(self): + emp = distribution_factory("univariate-empirical").fit([1.0, 2.0, 3.0]) + self.assertEqual(emp.cdf(-1.0, lower_bound=0.0), 0) # past the bound + self.assertEqual(emp.cdf(9.0, upper_bound=4.0), 1) + # inside the bound the cdf interpolates towards it + self.assertGreater(emp.cdf(0.5, lower_bound=0.0), 0) + self.assertLess(emp.cdf(3.5, upper_bound=4.0), 1) + self.assertGreaterEqual(emp.cdf_inverse(0.01, lower_bound=0.0), 0.0) + self.assertLessEqual(emp.cdf_inverse(0.99, upper_bound=4.0), 4.0) + + def test_empirical_extrapolates_below_the_smallest_record(self): + # with no lower bound given, a quantile under 1/(n+1) follows the line + # through the first two plotting positions, (0.25, 1.0) and (0.5, 2.0) + emp = distribution_factory("univariate-empirical").fit([1.0, 2.0, 3.0]) + self.assertAlmostEqual(emp.cdf_inverse(0.1), 0.4) + self.assertLess(emp.cdf_inverse(0.2), 1.0) # below the smallest record + + def test_empirical_cdf_inverse_rejects_bad_quantiles(self): + emp = distribution_factory("univariate-empirical").fit([1.0, 2.0, 3.0]) + for bad in (-0.1, 1.1): + with self.assertRaises(ValueError): + emp.cdf_inverse(bad) + + def test_empirical_handles_a_degenerate_sample(self): + # a resample can easily come out constant (or hold a single record): + # every quantile is then that value, and neither tail has a second + # point to take a slope from + for data in ([5.0], [5.0, 5.0, 5.0]): + emp = distribution_factory("univariate-empirical").fit(data) + for q in (0.0, 0.1, 0.5, 0.9, 1.0): + self.assertEqual(emp.cdf_inverse(q), 5.0, msg=f"{data}: {q}") + self.assertEqual(emp.cdf(4.0), 0) + self.assertEqual(emp.cdf(6.0), 1) + + def test_empirical_extrapolates_past_a_repeated_extreme(self): + # the largest value is repeated, so the upper extrapolation has to walk + # down to the bottom of that run to find a slope + emp = distribution_factory("univariate-empirical").fit( + [1.0, 2.0, 3.0, 3.0]) + self.assertGreater(emp.cdf_inverse(0.99), 3.0) + self.assertEqual(emp.cdf(9.0), 1) # the extrapolated line is clamped + self.assertEqual(emp.cdf(-9.0), 0) + + def test_interpolate_line(self): + line = statdist_distributions.interpolate_line(0.0, 1.0, 2.0, 5.0) + self.assertAlmostEqual(line(1.0), 3.0) + with self.assertRaises(ValueError): + statdist_distributions.interpolate_line(1.0, 0.0, 1.0, 5.0) + + def _discrete(self, pairs): + return distribution_factory("univariate-discrete")(OrderedDict(pairs)) + + def test_discrete_moments(self): + # a fair two-point distribution on {0, 2}: mean 1, variance 1 + fair = self._discrete([(0.0, 0.5), (2.0, 0.5)]) + self.assertAlmostEqual(fair.mean, 1.0) + self.assertAlmostEqual(fair.var, 1.0) + # and a three-point one, against E[X^2] - E[X]^2 + d = self._discrete([(1.0, 0.2), (2.0, 0.3), (5.0, 0.5)]) + mean = 0.2*1 + 0.3*2 + 0.5*5 + self.assertAlmostEqual(d.mean, mean) + self.assertAlmostEqual(d.var, 0.2*1 + 0.3*4 + 0.5*25 - mean**2) + self.assertGreater(d.var, 0.0) + + def test_discrete_validates_its_breakpoints(self): + with self.assertRaises(RuntimeError): # not a dict at all + distribution_factory("univariate-discrete")([(0.0, 1.0)]) + with self.assertRaises(RuntimeError): # values out of order + self._discrete([(2.0, 0.5), (1.0, 0.5)]) + for bad in ([(0.0, 0.25), (1.0, 0.25)], [(0.0, 0.9), (1.0, 0.9)]): + with self.assertRaises(ValueError): # probabilities not one + self._discrete(bad) + + def test_discrete_cdf_is_a_step_function(self): + d = self._discrete([(1.0, 0.2), (2.0, 0.3), (5.0, 0.5)]) + self.assertEqual(d.cdf(0.0), 0) + self.assertAlmostEqual(d.cdf(1.0), 0.2) + self.assertAlmostEqual(d.cdf(1.5), 0.2) # flat between breakpoints + self.assertAlmostEqual(d.cdf(2.0), 0.5) + self.assertAlmostEqual(d.cdf(4.9), 0.5) + self.assertAlmostEqual(d.cdf(5.0), 1.0) + self.assertAlmostEqual(d.cdf(6.0), 1.0) + self.assertAlmostEqual(d.rect_prob(1.0, 5.0), 0.8) + + def test_discrete_has_no_density_or_inverse(self): + d = self._discrete([(1.0, 0.5), (2.0, 0.5)]) + with self.assertRaises(RuntimeError): + d.pdf(1.0) + with self.assertRaises(RuntimeError): + d.cdf_inverse(0.5) + + def test_discrete_sample_one_draws_from_the_breakpoints(self): + d = self._discrete([(1.0, 0.25), (2.0, 0.75)]) + d.seed_reset(4) + draws = [d.sample_one() for _ in range(400)] + self.assertEqual(set(draws), {1.0, 2.0}) + self.assertAlmostEqual(draws.count(2.0)/len(draws), 0.75, places=1) + + @unittest.skipIf(not ipopt_available, "ipopt (nonlinear solver) not available") + def test_epispline_fit_inverse(self): + data = list(np.random.RandomState(2).normal(0, 1, size=100)) + epi = distribution_factory("univariate-epispline").fit(data) + self.assertLessEqual(epi.cdf_inverse(0.25), epi.cdf_inverse(0.75)) + + +#***************************************************************************** +class _RampDistribution(UnivariateDistribution): + """ A closed-form distribution for testing the base-class machinery. + + The density is 2x on [0, 1], so cdf(x) = x**2, cdf_inverse(q) = sqrt(q), + and the mean is 2/3. + """ + + def __init__(self, declare_support=True): + self.alpha = 0.0 + self.beta = 1.0 + params = [Parameter("slope", 2.0)] + if declare_support: + UnivariateDistribution.__init__(self, params, self.alpha, self.beta) + else: + UnivariateDistribution.__init__(self, params) + + @classmethod + def fit(cls, data): + return cls() + + def pdf(self, x): + if x < self.alpha or x > self.beta: + return 0.0 + return 2.0 * x + + +class _Interval: + """ The minimal interval protocol conditional_expectation expects. """ + + def __init__(self, a, b, cutouts=None): + self.a = a + self.b = b + if cutouts is not None: + self.cutouts = cutouts + + +class Test_statdist_base(unittest.TestCase): + """ The generic univariate machinery in base_distribution.py: the numeric + cdf and its inversion, expectations, sampling and parameter bookkeeping. """ + + def setUp(self): + self.d = _RampDistribution() + + def test_support_defaults_to_unbounded(self): + self.assertEqual((self.d.lower, self.d.upper), (0.0, 1.0)) + undeclared = _RampDistribution(declare_support=False) + self.assertEqual((undeclared.lower, undeclared.upper), + (-np.inf, np.inf)) + self.assertEqual(undeclared.dimension, 1) + + def test_cdf_integrates_the_density(self): + for x in (0.1, 0.5, 0.9): + self.assertAlmostEqual(self.d.cdf(x), x**2, places=5) + self.assertEqual(self.d.cdf(self.d.alpha), 0) + self.assertEqual(self.d.cdf(-1.0), 0) + self.assertEqual(self.d.cdf(self.d.beta), 1) + self.assertEqual(self.d.cdf(2.0), 1) + + def test_cdf_inverse_inverts_the_cdf(self): + for q in (0.1, 0.25, 0.5, 0.81): + self.assertAlmostEqual(self.d.cdf_inverse(q), math.sqrt(q), + places=3) + # the ends of the support, and quantiles that are not quantiles + self.assertEqual(self.d.cdf_inverse(0.0), self.d.alpha) + self.assertEqual(self.d.cdf_inverse(1.0), self.d.beta) + self.assertIsNone(self.d.cdf_inverse(-0.1)) + self.assertIsNone(self.d.cdf_inverse(1.1)) + + def test_cdf_is_cached_per_tolerance(self): + # the cdf is memoized, and a different accuracy is a different question + self.assertAlmostEqual(self.d.cdf(0.5, epsabs=1e-2), 0.25, places=2) + self.assertAlmostEqual(self.d.cdf(0.5, epsabs=1e-12), 0.25, places=9) + + def test_mean_and_region_expectation(self): + self.assertAlmostEqual(self.d.mean(), 2/3, places=5) + self.assertAlmostEqual(self.d.region_expectation((0.0, 1.0)), 2/3, + places=5) + self.assertAlmostEqual(self.d.region_expectation((0.0, 0.5)), 1/12, + places=5) + self.assertAlmostEqual(self.d.region_probability((0.0, 0.5)), 0.25, + places=5) + self.assertAlmostEqual(self.d.region_probability((0.0, 1.0)), 1.0, + places=5) + + def test_region_arguments_are_validated(self): + with self.assertRaises(ValueError): + self.d.region_expectation((0.75, 0.25)) # upper below lower + # a region has to be a tuple, and the complaint about that has to + # survive the memoization wrapper (a list is not hashable) + for not_a_region in ([0.0, 1.0], "region"): + with self.assertRaises(TypeError): + self.d.region_expectation(not_a_region) + with self.assertRaises(ValueError): + self.d.region_probability(not_a_region) + + def test_conditional_expectation(self): + # conditioning on the whole support is just the mean + self.assertAlmostEqual( + self.d.conditional_expectation(_Interval(0.0, 1.0)), 2/3, places=3) + # cutting the lower half out conditions on the upper half, which pulls + # the expectation up; E[X | X > median] = (2/3)(1 - 0.5**1.5)/0.5 + upper_half = self.d.conditional_expectation( + _Interval(0.0, 1.0, cutouts=[_Interval(0.0, 0.5)])) + self.assertAlmostEqual(upper_half, (2/3)*(1 - 0.5**1.5)/0.5, places=3) + self.assertGreater(upper_half, 2/3) + + def test_log_likelihood(self): + data = [0.25, 0.5, 0.75] + self.assertAlmostEqual(self.d.log_likelihood(data), + sum(math.log(2*x) for x in data)) + + def test_sampling_stays_in_the_support(self): + # the inversion is numeric, so allow it a little slack at the ends + slack = 1e-3 + self.d.seed_reset(12) + for _ in range(20): + drawn = self.d.sample_one() + self.assertGreaterEqual(drawn, self.d.alpha - slack) + self.assertLessEqual(drawn, self.d.beta + slack) + for _ in range(10): + drawn = self.d.sample_on_interval(0.25, 0.75) + self.assertGreaterEqual(drawn, 0.25 - slack) + self.assertLessEqual(drawn, 0.75 + slack) + # a quantile range maps to the matching range of values + between = self.d.sample_between_quantiles(0.1, 0.2) + self.assertGreaterEqual(between, math.sqrt(0.1) - slack) + self.assertLessEqual(between, math.sqrt(0.2) + slack) + + def test_str_and_repr_name_the_parameters(self): + self.assertIn("slope", str(self.d)) + self.assertIn("2.0", str(self.d)) + self.assertEqual(repr(self.d), "Distribution(_RampDistribution)") + + def test_parameter_bookkeeping(self): + p = Parameter("mean", 3.0, bounds=(0, None)) + self.assertTrue(p.instantiated) # it has a value + self.assertEqual(p.bounds, (0, None)) + self.assertIs(p.kind, float) + self.assertEqual(repr(p), "Parameter(mean,3.0)") + self.assertEqual(str(p), repr(p)) + unset = Parameter("variance") + self.assertFalse(unset.instantiated) # and this one does not + unset.set_value(2.5) + self.assertEqual(unset.value, 2.5) + self.assertTrue(unset.instantiated) + + @unittest.skipIf(not matplotlib_available, "matplotlib is not installed") + def test_plot_writes_a_file(self): + import matplotlib + matplotlib.use("Agg") # no display in a test run + with tempfile.TemporaryDirectory() as tmpdir: + plot_dir = os.path.join(tmpdir, "plots") + self.d.plot(output_file="ramp.png", title="ramp", xlabel="x", + ylabel="density", output_directory=plot_dir) + self.assertTrue(os.path.exists(os.path.join(plot_dir, "ramp.png"))) + # an unbounded support falls back to a [-5, 5] window, and the + # directory this time already exists + _RampDistribution(declare_support=False).plot( + plot_cdf=False, output_file="unbounded.png", + output_directory=plot_dir) + self.assertTrue( + os.path.exists(os.path.join(plot_dir, "unbounded.png"))) + + +#***************************************************************************** +class Test_statdist_utilities(unittest.TestCase): + """ The memoization helpers and the argv context manager in + statdist/utilities.py. """ + + def test_memoize_caches_by_value(self): + calls = [] + + @statdist_utilities.memoize + def total(xs, offset=0): + calls.append(1) + return sum(xs) + offset + + self.assertEqual(total([1, 2, 3]), 6) + self.assertEqual(total([1, 2, 3]), 6) + self.assertEqual(len(calls), 1) # the second call was cached + # an unhashable list argument normalizes to the tuple's key + self.assertEqual(total((1, 2, 3)), 6) + self.assertEqual(len(calls), 1) + self.assertEqual(total([1, 2, 3], offset=10), 16) + self.assertEqual(len(calls), 2) + + def test_memoize_normalizes_dictionary_arguments(self): + calls = [] + + @statdist_utilities.memoize + def size(mapping): + calls.append(1) + return len(mapping) + + self.assertEqual(size({"a": 1, "b": 2}), 2) + self.assertEqual(size({"b": 2, "a": 1}), 2) # equal dict, new object + self.assertEqual(len(calls), 1) + + def test_normalize_args_maps_positionals_to_names(self): + def f(a, b, c=0): + return a + + args = statdist_utilities.normalize_args(f, (1, [2, 3]), {"c": {"k": 4}}) + self.assertEqual(args["a"], 1) + self.assertEqual(args["b"], (2, 3)) # list -> tuple + self.assertEqual(args["c"], (("k", 4),)) # dict -> sorted pairs + + def test_memoize_method_caches_per_instance(self): + class Counter: + def __init__(self): + self.calls = 0 + + @statdist_utilities.memoize_method + def squared(self, x): + self.calls += 1 + return x * x + + one, two = Counter(), Counter() + self.assertEqual(one.squared(3), 9) + self.assertEqual(one.squared(3), 9) + self.assertEqual(one.calls, 1) + self.assertEqual(two.squared(3), 9) # a cache of its own + self.assertEqual(two.calls, 1) + # reached through the class the method is the undecorated one + self.assertEqual(Counter.squared(one, 4), 16) + self.assertEqual(one.calls, 2) + + def test_memoize_method_keys_on_keyword_values(self): + class Rounder: + def __init__(self): + self.calls = 0 + + @statdist_utilities.memoize_method + def value(self, x, places=2): + self.calls += 1 + return round(x, places) + + r = Rounder() + self.assertEqual(r.value(1.23456, places=2), 1.23) + self.assertEqual(r.value(1.23456, places=4), 1.2346) + self.assertEqual(r.calls, 2) # not one answer for both + self.assertEqual(r.value(1.23456, places=4), 1.2346) + self.assertEqual(r.calls, 2) # and now it is cached + + def test_memoize_method_passes_unhashable_arguments_through(self): + class Sizer: + def __init__(self): + self.calls = 0 + + @statdist_utilities.memoize_method + def size(self, thing): + self.calls += 1 + if not isinstance(thing, tuple): + raise TypeError("tuples only") + return len(thing) + + s = Sizer() + self.assertEqual(s.size((1, 2, 3)), 3) + # an unhashable argument cannot be cached, but the method still runs + # and its own error is what comes back + with self.assertRaises(TypeError): + s.size([1, 2, 3]) + self.assertEqual(s.calls, 2) + + def test_set_arguments_restores_argv(self): + saved = list(sys.argv) + with statdist_utilities.set_arguments(["prog", "--flag"]): + self.assertEqual(sys.argv, ["prog", "--flag"]) + self.assertEqual(sys.argv, saved) + + +#***************************************************************************** +class Test_empirical_examples(unittest.TestCase): + """ Empirical methods on the statdist-dependent examples (farmer, cvar). + + These could not live in test_boot_sp.py because importing farmer/cvar pulls + in statdist; the methods themselves are the empirical ones. + """ + + @unittest.skipIf(not solver_available, "no solver is available") + def test_farmer_empirical_wellformed(self): + module = boot_utils.module_name_to_module("farmer") + xhat = boot_utils.compute_xhat(_make_farmer_cfg(), module) + self.assertIn("ROOT", xhat) + for method in ["Classical_quantile", "Bagging_with_replacement"]: + # every rank participates in the collective inside compute_ci + res = boot_sp.compute_ci(_make_farmer_cfg(method), module, xhat) + self.assertEqual(len(res), 6) + if my_rank == 0: + for ci in res[:3]: + self.assertLessEqual(ci[0], ci[1], msg=f"{method}: {ci}") + else: + self.assertEqual(res, (None, None, None, None, None, None)) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_empirical_wellformed(self): + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Classical_quantile") + xhat = boot_utils.compute_xhat(cfg, module) + self.assertIn("ROOT", xhat) + res = boot_sp.compute_ci(_make_cvar_cfg("Classical_quantile"), module, xhat) + self.assertEqual(len(res), 6) + if my_rank == 0: + for ci in res[:3]: + self.assertLessEqual(ci[0], ci[1]) + else: + self.assertEqual(res, (None, None, None, None, None, None)) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_empirical_deterministic(self): + # same cfg twice must give the same interval (seeded streams) + module = boot_utils.module_name_to_module("cvar") + xhat = boot_utils.compute_xhat(_make_cvar_cfg("Classical_gaussian"), module) + # both runs are collectives on every rank; only rank 0 gets real values + r1 = boot_sp.compute_ci(_make_cvar_cfg("Classical_gaussian"), module, xhat) + r2 = boot_sp.compute_ci(_make_cvar_cfg("Classical_gaussian"), module, xhat) + if my_rank == 0: + for a, b in zip(list(r1[0]), list(r2[0])): + self.assertEqual(round_pos_sig(a, 6), round_pos_sig(b, 6)) + + +#***************************************************************************** +class Test_smoothed(unittest.TestCase): + """ Smoothed methods (kernel/bagging need no nonlinear solver). """ + + def _check_gap_ci(self, result, method): + # rank-0 result is (ci_gap_two_sided, center_gap); non-root is (None, None) + if my_rank == 0: + ci_gap, center_gap = result + self.assertEqual(len(ci_gap), 2) + self.assertTrue(math.isfinite(center_gap), msg=method) + self.assertLessEqual(ci_gap[0], ci_gap[1], msg=f"{method}: {ci_gap}") + else: + self.assertEqual(result, (None, None)) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_smoothed_bagging(self): + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_bagging") + xhat = boot_utils.compute_xhat(cfg, module) + result = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + self._check_gap_ci(result, "Smoothed_bagging") + + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_smoothed_kernel(self): + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_boot_kernel") + xhat = boot_utils.compute_xhat(cfg, module) + result = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + self._check_gap_ci(result, "Smoothed_boot_kernel") + + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_smoothed_kernel_quantile(self): + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_boot_kernel_quantile") + xhat = boot_utils.compute_xhat(cfg, module) + result = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + self._check_gap_ci(result, "Smoothed_boot_kernel_quantile") + + @unittest.skipIf(not solver_available, "no solver is available") + def test_user_boot_smoothed(self): + # the end-user entry point routes smoothed methods and clamps ci_gap[0] + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_bagging") + result = user_boot.main_routine(cfg, module) + if my_rank == 0: + ci_gap, center_gap = result + self.assertGreaterEqual(ci_gap[0], 0.0) + self.assertLessEqual(ci_gap[0], ci_gap[1]) + else: + self.assertEqual(result, (None, None)) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_simulate_smoothed_coverage(self): + # the smoothed coverage harness (this exercises the section-4.3 + # compute_xhat fix: no xhat file, so it computes xhat internally) + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_bagging", reps=2) + result = simulate_boot.main(cfg, module) + if my_rank == 0: + cov_two, cov_one, ci_len, run_time = result + self.assertGreaterEqual(cov_two, 0.0) + self.assertLessEqual(cov_two, 1.0) + self.assertGreaterEqual(cov_one, cov_two) # one-sided covers at least as often + self.assertEqual(len(ci_len), cfg.coverage_replications) + else: + self.assertEqual(result, (None, None, None, None)) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_smoothed_bootstrap_draws_are_disjoint_and_fitted(self): + # Two properties of the smoothed bootstrap that are easy to lose: + # (1) every batch is an independent set of draws from the fitted + # distribution, so the per-batch record blocks are pairwise + # disjoint and disjoint from the center's block. Overlapping + # blocks reuse draws and collapse the estimated spread. + # (2) the center is drawn from the *fitted* distribution, not from + # the raw sample; drawing it raw makes it the purely empirical + # point estimate. + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_boot_kernel") + xhat = boot_utils.compute_xhat(cfg, module) + + pools = [] + fitted_at_center = [] + real_eval = boot_sp.evaluate_scenarios + real_center = smoothed_boot_sp.center_smoothed + + # the smoothed callers never pass a communicator, so the spy does not + # need one either -- and not taking one keeps this working whether or + # not evaluate_scenarios has grown an mpicomm argument + def spy_eval(cfg_, module_, scenarios, xhat_, duplication=True): + pools.append(list(scenarios)) + return real_eval(cfg_, module_, scenarios, xhat_, + duplication=duplication) + + def spy_center(cfg_, module_, xhat_): + fitted_at_center.append(cfg_.use_fitted) + return real_center(cfg_, module_, xhat_) + + boot_sp.evaluate_scenarios = spy_eval + smoothed_boot_sp.center_smoothed = spy_center + try: + smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + finally: + boot_sp.evaluate_scenarios = real_eval + smoothed_boot_sp.center_smoothed = real_center + + self.assertEqual(fitted_at_center, [True]) # (2) + + # pools[0] is the center; the rest are this rank's batches + center_pool, batch_pools = set(pools[0]), [set(p) for p in pools[1:]] + self.assertEqual(len(center_pool), cfg.smoothed_center_sample_size) + for i, bp in enumerate(batch_pools): # (1) + self.assertEqual(len(bp), cfg.sample_size) + self.assertEqual(bp & center_pool, set(), + msg=f"batch {i} reuses the center's draws") + for j, other in enumerate(batch_pools[i + 1:], start=i + 1): + self.assertEqual(bp & other, set(), + msg=f"batches {i} and {j} share draws") + + @unittest.skipIf(not ipopt_available, "ipopt (nonlinear solver) not available") + @unittest.skipIf(not solver_available, "no solver is available") + def test_cvar_smoothed_epi(self): + module = boot_utils.module_name_to_module("cvar") + cfg = _make_cvar_cfg("Smoothed_boot_epi") + xhat = boot_utils.compute_xhat(cfg, module) + result = smoothed_boot_sp.compute_smoothed_ci(cfg, module, xhat) + self._check_gap_ci(result, "Smoothed_boot_epi") + + +#***************************************************************************** +class Test_multi_knapsack(unittest.TestCase): + """ Smoke test the multi_knapsack example (deterministic-data-json path). """ + + def test_import_and_data(self): + module = boot_utils.module_name_to_module("multi_knapsack") + self.assertTrue(hasattr(module, "scenario_creator")) + self.assertTrue(hasattr(module, "data_sampler")) + self.assertTrue(hasattr(module, "xhat_generator")) + + @unittest.skipIf(not solver_available, "no solver is available") + def test_multi_knapsack_empirical(self): + module = boot_utils.module_name_to_module("multi_knapsack") + cfg = boot_utils._process_module("multi_knapsack") + cfg.module_name = "multi_knapsack" + cfg.max_count = 60 + cfg.candidate_sample_size = 3 + cfg.sample_size = 15 + cfg.subsample_size = 5 + cfg.nB = 6 + cfg.alpha = 0.1 + cfg.seed_offset = 100 + cfg.xhat_fname = "None" + cfg.optimal_fname = "None" + cfg.trace_fname = None + cfg.coverage_replications = 2 + cfg.solver_name = solver_name + cfg.boot_method = "Bagging_with_replacement" + cfg.deterministic_data_json = MK_DATA + xhat = boot_utils.compute_xhat(cfg, module) + self.assertIn("ROOT", xhat) + res = boot_sp.compute_ci(cfg, module, xhat) + self.assertEqual(len(res), 6) + + +if __name__ == '__main__': + unittest.main() diff --git a/run_coverage.bash b/run_coverage.bash index 781de1cac..f2bdd86c8 100755 --- a/run_coverage.bash +++ b/run_coverage.bash @@ -226,6 +226,9 @@ run_phase "test_boot_sp (serial)" \ run_phase "test_boot_sp_simulate (serial)" \ coverage run --rcfile=.coveragerc mpisppy/tests/test_boot_sp_simulate.py +run_phase "test_boot_sp_smoothed (serial)" \ + coverage run --rcfile=.coveragerc mpisppy/tests/test_boot_sp_smoothed.py + run_phase "test_gradient_rho (spawns mpiexec)" \ coverage run --rcfile=.coveragerc mpisppy/tests/test_gradient_rho.py @@ -249,6 +252,9 @@ run_phase "test_boot_sp (mpiexec -np 2)" \ run_phase "test_boot_sp_simulate (mpiexec -np 2)" \ mpiexec -np 2 coverage run --rcfile="$PROJ_DIR/.coveragerc" -m mpi4py mpisppy/tests/test_boot_sp_simulate.py +run_phase "test_boot_sp_smoothed (mpiexec -np 2)" \ + mpiexec -np 2 coverage run --rcfile="$PROJ_DIR/.coveragerc" -m mpi4py mpisppy/tests/test_boot_sp_smoothed.py + run_phase "test_cg_main (serial)" \ coverage run --rcfile=.coveragerc mpisppy/tests/test_cg_main.py