From 1af4768fc2c91affcddb5100fbfed229be650e90 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 12:25:25 -0700 Subject: [PATCH 01/17] boot-sp PR-2: trimmed statdist univariate distribution library Add mpisppy/confidence_intervals/bootsp/statdist/, a trimmed port of the statdist library holding only the univariate distributions the smoothed bootstrap methods need (uniform, normal, student-t, Gaussian kernel, epi-spline, empirical, discrete) plus their support modules. The multivariate machinery (copula.py, vine.py, bicop.py and the multivariate distribution classes) is dropped; that also removes the `from scipy.stats import mvn` import (gone in scipy 1.14) and the gosm hook. scipy is imported lazily via pyomo.common.dependencies so the empirical bootstrap path stays scipy-free, and matplotlib is imported inside the plotting methods so it stays optional. distribution_factory now scans the registry once. .ruff.toml relaxes the legacy-style rules this ported library predates, scoped to the statdist subpackage. Co-Authored-By: Claude Opus 4.8 --- .ruff.toml | 15 + .../bootsp/statdist/README.md | 34 + .../bootsp/statdist/__init__.py | 13 + .../bootsp/statdist/base_distribution.py | 1099 +++++++++++++++++ .../bootsp/statdist/distribution_factory.py | 105 ++ .../bootsp/statdist/distributions.py | 854 +++++++++++++ .../bootsp/statdist/sampler.py | 34 + .../bootsp/statdist/splines.py | 543 ++++++++ .../bootsp/statdist/utilities.py | 194 +++ 9 files changed, 2891 insertions(+) create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/README.md create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/__init__.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/distribution_factory.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/distributions.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/sampler.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/splines.py create mode 100644 mpisppy/confidence_intervals/bootsp/statdist/utilities.py 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/mpisppy/confidence_intervals/bootsp/statdist/README.md b/mpisppy/confidence_intervals/bootsp/statdist/README.md new file mode 100644 index 000000000..4f5758798 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/README.md @@ -0,0 +1,34 @@ +# 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`, and the +multivariate distribution classes in `distributions.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..ca3dc001d --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py @@ -0,0 +1,1099 @@ +############################################################################### +# 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 +from functools import wraps +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 + self.instantiated = value is 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 + + 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) + + +class MultivariateDistribution(BaseDistribution): + """ + This class is an abstract base class for all multivariate distributions + TODO: This docstring should be improved greatly!! + """ + __metaclass__ = ABCMeta + + def __init__(self, dimension, dimkeys=None, parameters=None, lower=None, + upper=None, bounds=None): + """ + Args: + dimension (int): The dimension of the distribution + dimkeys (List): A list of the names of the dimensions, by default, + these will just be the indices. If passed in, this will enable + you to refer to values by the dimension name in certain + functions + parameters (list[Parameter]): A list of parameters for the + distribution + lower (list[float]): A list of the lower bounds of the support + of the distribution + upper (list[float]): A list of the upper bounds of the support + of the distribution + bounds (list|dict): A colection of bounds on the support for + each dimension. We assume the support is on a rectangular + region. If it is passed as a dictionary it should map + dimension names to ordered pairs of lower and upper bounds. A + None indicates that there is no lower or upper bound for a + given dimension. + """ + BaseDistribution.__init__(self, dimension, parameters) + self.ndim = dimension + if dimkeys is None: + # We default to using the integers if no dimkeys are passed in. + self.dimkeys = list(range(self.ndim)) + else: + self.dimkeys = dimkeys + + if lower is None: + self.lower = [None for _ in range(dimension)] + else: + self.lower = lower + if upper is None: + self.upper = [None for _ in range(dimension)] + else: + self.upper = upper + + if bounds is None: + self.bounds = [(-np.inf, np.inf)] * dimension + elif isinstance(bounds, list): + self.bounds = bounds + elif isinstance(bounds, dict): + self.bounds = [bounds[dim] for dim in self.dimkeys] + + def pdf(self, *xs): + raise NotImplementedError + + def log_likelihood(self, data): + """ + This method will return the log likelihood of the observed data + given the fitted model. + + This method just naively computes the pdf and applies the logarithm. + It would be more efficient in subclasses to find an expression for + the log-likelihood. + + The argument data can either be a list of vectors for each dimension + of the data or it can be a dictionary mapping dimension names to the + corresponding vector of data. + + Args: + data (list[list[float]] | dict[list[float]]): The observed values + Returns: + float: The computed log-likelihood + """ + if isinstance(data, dict): + vects = [data[dimkey] for dimkey in dimkeys] + else: + vects = data + + return sum(np.log(self.pdf(*xs)) for xs in zip(*vects)) + + def plot(self, func, lower=None, upper=None): + """ + Args: + func (str): The function to plot, either 'pdf' or 'cdf' + lower (list[float]): A list of the lower bounds for the plot, + will default to the lower bounds of the support if None + upper (list[float]): A list of the upper bounds of the plot + will default to the upper bounds of the support if None + """ + if self.dimension != 2: + raise ValueError("This plot method is only functional for 2-d " + "distributions.") + + if lower is None: + lower = self.lower + if lower[0] is None: + lower[0] = -5 + if lower[1] is None: + lower[1] = 5 + if upper is None: + upper = self.upper + if upper[0] is None: + upper[0] = -5 + if upper[1] is None: + upper[1] = 5 + + import matplotlib.pyplot as plt + fig = plt.figure() + ax = fig.gca(projection='3d') + + X = np.arange(lower[0], upper[0], 0.1) + Y = np.arange(lower[1], upper[1], 0.1) + X, Y = np.meshgrid(X, Y) + + Z = np.zeros_like(X) + for i, row in enumerate(X): + for j, x in enumerate(row): + y = Y[i,j] + if func == 'pdf': + z = self.pdf(x, y) + elif func == 'cdf': + z = self.cdf(x, y) + + Z[i,j] = z + + ax.plot_surface(X, Y, Z) + ax.set_xlim(lower[0], upper[0]) + ax.set_ylim(lower[1], upper[1]) + return ax + + def contour_plot(self, func, lower=None, upper=None): + """ + Args: + func (str): The function to plot, either 'pdf' or 'cdf' + lower (list[float]): A list of the lower bounds for the plot, + will default to the lower bounds of the support if None + upper (list[float]): A list of the upper bounds of the plot + will default to the upper bounds of the support if None + """ + if self.dimension != 2: + raise ValueError("This plot method is only functional for 2-d " + "distributions.") + + if lower is None: + lower = self.lower + if upper is None: + upper = self.upper + + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + + X = np.arange(lower[0], upper[0], 0.1) + Y = np.arange(lower[1], upper[1], 0.1) + X, Y = np.meshgrid(X, Y) + + Z = np.zeros_like(X) + for i, row in enumerate(X): + for j, x in enumerate(row): + y = Y[i,j] + if func == 'pdf': + z = self.pdf(x, y) + elif func == 'cdf': + z = self.cdf(x, y) + + Z[i,j] = z + + ax.contour(X, Y, Z) + ax.set_xlim(lower[0], upper[0]) + ax.set_ylim(lower[1], upper[1]) + return ax + + @memoize_method + def rect_prob(self, lowerdict, upperdict): + tempdict = dict.fromkeys(self.dimkeys) + def f(n): + + # recursive function that will calculate the cdf + # It has a structure of binary tree + if n == 0: + return self.cdf(tempdict) + else: + tempdict[self.dimkeys[n - 1]] = upperdict[self.dimkeys[n - 1]] + leftresult = f(n - 1) + tempdict[self.dimkeys[n - 1]] = lowerdict[self.dimkeys[n - 1]] + rightresult = f(n - 1) + return leftresult - rightresult + + return f(self.dimension) + + def marginal(self, ys, bounds = None, error_tolerance=None): + """ + This function will evaluate the marginal distribution of the joint + cdf which is composed of the dimensions passed in through ys. + It will evaluate it at the point ys. + + Args: + ys (dict): A dictionary mapping dimension names to their + corresponding values + error_tolerance (int): Value to increase the error tolerance + of the integration process by powers of 10 + Returns: + float: The value of the marginal + """ + + other_dims = [(i, dim) for i, dim in enumerate(self.dimkeys) + if dim not in ys] + + point_dict = ys.copy() + + def pdf_x(*xs): + for (_, dim), x in zip(other_dims, xs): + point_dict[dim] = x + #print(self.pdf(point_dict)) + return self.pdf(point_dict) + if bounds == None: + bounds = [self.bounds[i] for i, _ in other_dims] + + if error_tolerance: + tol = error_tolerance + else: + tol = 0 + + return scipy.integrate.nquad(pdf_x, bounds, opts={'epsabs': (1.49e-08 * (10**tol)), 'epsrel': (1.49e-08 * (10**tol))} )[0] + + def conditional_pdf(self, xs, cond_xs, marginal_cdf=None): + """ + This will evaluate the conditional pdf at the point xs given + that the dimensions in cond_names + + Args: + xs (dict): A dictionary mapping dimension names to values + cond_xs (dict): A dictionary mapping the dimension names + of the conditioned variables to their values + Returns: + float: The value fo the conditional pdf + """ + if marginal_cdf == None: + marg = self.marginal(cond_xs) + else: + marg = marginal_cdf + + point_dict = {} + for dim, x in xs.items(): + point_dict[dim] = x + for dim, x in cond_xs.items(): + point_dict[dim] = x + return self.pdf(point_dict) / marg + + def conditional_cdf(self, xs, cond_xs, marginal_cdf = None): + """ + This will evaluate the conditional cdf at the point xs given + that the dimensions in cond_names are set to the values in cond_xs. + + Args: + xs (dict): A dictionary mapping dimension names to values + cond_xs (dict): A dictionary mapping the dimension names + of the conditioned variables to their values + Returns: + float: The value fo the conditional cdf + """ + bounds = [] + + dimkeys = list(xs.keys()) + + for dim in dimkeys: + dim_index = self.dimkeys.index(dim) + lower_bound = self.bounds[dim_index][0] + bounds.append([lower_bound, xs[dim]]) + + point_dict = cond_xs.copy() + def f(*xs): + for dim, x in zip(dimkeys, xs): + point_dict[dim] = x + return self.pdf(point_dict) + + if marginal_cdf == None: + marg = self.marginal(cond_xs) + else: + marg = marginal_cdf + + try: + return scipy.integrate.nquad(f, bounds)[0] / marg + except: + return 0 + + def conditional_cdf_inverse(self, cond_xs, cdf_value, dim, marginal, + capacity = 4000, n = 100, + method = 'combination', xtol = 0.001): + """ + This function computes the inverse of a conditional cdf value + conditioned on a given point. Therefore 3 different methods are + provided: + - linear interpolation: The conditional cdf is evaluated at several + points in a given interval. Since two points which conditional cdfs + wrap the cdf_value, a linear interpolation between these two + points is used to compute the inverse of the cdf_value. + - bisection: The bisection method from scipy is used to solve the + equation 0 = conditional_cdf - cdf_value. + - combination of both: First a bisection method is used to find the + two wrapping points like in the linear interpolation method. After + that a linear interpolation is used to compute the inverse. + + Args: + cond_xs (dict): A dictionary mapping dimension the dimension + names of the conditioned variables to their values. + cdf_value: The value you want to compute the inverse for. + dim (int or str): The name of the dimension you want to get the + inverse for (e.g. F(X|Y=500) = 0.2: You want to compute the + value of X under the condition that Y=500, so that F equals + 0.2. In that case "dim" equals X.). + marginal (distribution like): The marginal of dimension dim. + capacity (float): The capacity for that day. + n (int): The number of intersection of the interval + [-capacity, capacity], which specify the points which are + evaluated for the linear interpolation method. The number + specifies also a break criteria for the bisection part in the + combined method. + method (str): The method you want to use. "default" refers to the + linear interpolation method, "bisect" to the bisection method + and "combination" to the combined method. + xtol (float): The tolerance for the bisection method. + (break criteria) + + Returns: + The inverse value of the passed in cdf_value conditioned on the + point cond_xs. + """ + marginal_cdf = self.marginal(cond_xs) # This value is needed a lot. So + # it is computed here once. + + if method == 'default': + """ + For the default or linear interpolation method first a list of + points are created. These points are evaluated one after the other + with the conditional_cdf function. After that it is checked, if the + given cdf_value is wrapped by two consecutive points' conditional + cdf. If thats the case, these two points and there conditional cdfs + are used to compute a linear interpolation between them. This + linear interpolation then is used to compute the inverse of the + given cdf_value. Because the conditional_cdf lives in the copula + space (which is [0,1]^n), the points have to be converted to [0,1]. + For the purpose of getting power values as a return, the computed + inverse values have to be transformed back in the end. + """ + points = np.linspace(-capacity, capacity, n) + x = [] + for point in points: + x.append(marginal.cdf(point)) + y = [] + j = 0 + for i in x: + xs = {dim: i} + yi = self.conditional_cdf(xs, cond_xs) + y.append(yi) + if (j==0) and (yi > cdf_value): + inverse = marginal.cdf(-capacity) + break + elif (j != 0): + if y[j-1] <= cdf_value <= y[j]: #linear interpolation + lin = scipy.interpolate.interp1d([y[j-1], y[j]], [x[j-1], x[j]]) + inverse = lin(cdf_value) #computing the inverse + break + j += 1 + else: + inverse = marginal.cdf(capacity) + return marginal.cdf_inverse(inverse) + elif method == 'bisect': + """ + In this method a help function is defined. After that the root + of this function is computed using the bisection mehtod from + scipy. For more information see the scipy documentation. + The transformation of the values is done for the same reason like + above. + """ + def help(d): + dict = {dim: d} + return self.conditional_cdf(dict, cond_xs, + marginal_cdf=marginal_cdf) \ + - cdf_value + if help(marginal.cdf(-capacity)) > 0: + inverse = -capacity + elif help(marginal.cdf(capacity)) < 0: + inverse = capacity + else: + inverse = marginal.cdf_inverse(scipy.optimize.bisect(help, + marginal.cdf(-capacity), + marginal.cdf(capacity), + xtol=xtol)) + + return inverse + + + elif method == 'combination': + """ + In this method not every single point is evaluated. There is some- + thing like a bisection method used to find faster the wrapping + points. After they are found, the linear interpolation is used + to compute the inverse value. + The transformation of the values is done for the same reason like + above. + """ + if capacity is None: + capacity = 0 + l = -capacity + u = capacity + l_cdf = self.conditional_cdf({dim: marginal.cdf(l)}, cond_xs, + marginal_cdf=marginal_cdf) + u_cdf = self.conditional_cdf({dim: marginal.cdf(u)}, cond_xs, + marginal_cdf=marginal_cdf) + if l_cdf >= cdf_value: + print('cdf', cdf_value) + print('lower', l_cdf) + return l + elif u_cdf <= cdf_value: + print('cdf', cdf_value) + print('upper', u_cdf) + return u + tol = (capacity * 2) / n + k = 0 + while ((u - l) > tol) and (k < n): + m = (u + l) / 2 + m_cdf = self.conditional_cdf({dim: marginal.cdf(m)}, cond_xs, + marginal_cdf=marginal_cdf) + if m_cdf < cdf_value: + l = m + l_cdf = m_cdf + elif m_cdf > cdf_value: + u = m + u_cdf = m_cdf + else: + return m + k = k + 1 + lin = scipy.interpolate.interp1d([l_cdf, u_cdf], [marginal.cdf(l), marginal.cdf(u)]) + return marginal.cdf_inverse(lin(cdf_value)) + + +def fit_wrapper(method): + """ + This is a function decorator which will wrap the fit method for + multivariate distributions. It will allow for data to be passed using + a dictionary mapping names to lists of data. + + Internally this transforms the data into a lists of lists and then fits + the distribution to that data. Then it assigns to the dimkeys attribute + the list of names. + + Args: + method: The class method fit of a multivariate distribution + Returns: + method: The modified method to handle dictionaries of input data + """ + @wraps(method) + def fit(cls, data, dimkeys=None, **kwargs): + """ + This function converts the dictionary into a list, passes it to the + fit method and then assigns to the distribution the dimkeys attribute. + """ + vectors = [] + if isinstance(data, dict): + for key in dimkeys: + vectors.append(data[key]) + else: + vectors = data + + distribution = method(cls, vectors, dimkeys, **kwargs) + return distribution + + return fit + + +def accepts_dict(method): + """ + This function decorator will allow any of the methods which accept separate + values for each dimension to also accept a dictionary which has keys + mapping to each dimension. + + For example, the pdf for any distribution generally has the prototype + def pdf(self, *x): + This decorator will unpack the dictionary into its respective dimensions + and pass it to the function. + + The function that this decorator is applied to must have a prototype like + def f(self, *x) + + This will enable you to call a function in the following three ways. + + Suppose distr is a Distribution with distr.dimkeys = ['foo', 'bar', 'baz'] + If pdf is decorated with accepts_dict, we can call it like so + 1) distr.pdf(1, 2, 3) + 2) distr.pdf(foo=1, bar=2, baz=3) + 3) distr.pdf({'foo': 1, 'bar': 2, 'baz': 3}) + + Args: + method: The method accepting the different values for each dimensions + Returns: + method: The modified method to handle dictionaries of input data + """ + @wraps(method) + def f(self, *xs, **kwargs): + if xs: + # If xs is passed in, we check if the user passed it as each + # dimension separately or as a dictionary + if isinstance(xs[0], dict): + # If the first element of xs, is a dict, assume only element. + value_dict = xs[0] + values = [value_dict[key] for key in self.dimkeys] + else: + # Otherwise, it is a list of the values at each dimension + values = xs + else: + # Otherwise, we expect the values to be passed as keyword args. + values = [kwargs[key] for key in self.dimkeys] + return method(self, *values) + + return f + + +def returns_dict(method): + """ + This function decorator will allow descendants of MultivariateDistribution + which have methods which return values for each dimension to instead + return a dictionary of values mapping dimension name to value. + + This adds an as_dict argument which if set to True, will pack the return + value into a dictionary assuming the order is in that of the dimkeys + attribute of the distribution. + + The as_dict argument must be passed by keyword. + + Args: + method: The method which returns a list of values for each dimension + Returns: + method: The modified method to return a dictionary if specified to + """ + + @wraps(method) + def f(self, *pargs, as_dict=False, **kwargs): + values = method(self, *pargs, **kwargs) + if as_dict: + output = {key: value for key, value in zip(self.dimkeys, values)} + return output + else: + return values + + return f + + +def params_as_args(arg_names): + def decorator(method): + """ + + """ + @wraps(method) + def f(cls, x, params=None): + if params is None: + params = {} + for name in arg_names: + value = getattr(cls, name).value + if value is None: + message = """The {} parameter is unset. To use this method + it must be either called from an instance of + the distribution class or it must be called + directly from the class with a dictionary of + the parameters passed + with the params keyword.""".format(name) + raise ValueError(message) + params[name] = value + return method(cls, x, params) + return f + return decorator + +def params_as_args2(arg_names): + def decorator(method): + """ + + """ + @wraps(method) + def f(cls, x, y, params=None): + if params is None: + params = {} + for name in arg_names: + value = getattr(cls, name).value + if value is None: + message = """The {} parameter is unset. To use this method + it must be either called from an instance of + the distribution class or it must be called + directly from the class with a dictionary of + the parameters passed + with the params keyword.""".format(name) + raise ValueError(message) + params[name] = value + return method(cls, x, y, params) + return f + return decorator 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..1c5350504 --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -0,0 +1,854 @@ +############################################################################### +# 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 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) + + 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: + 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!') + # 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] + j = n - 1 + while self.input_data[j] == self.input_data[j - 1]: + j -= 1 + if j - 1 == -len(self.input_data): + print("Warning: all input values are the same (", + self.input_data[j], ")") + return self.input_data[j] + # 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 = self.mean*self.mean - Esqsum + if 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 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..7e64f6a4f --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/statdist/utilities.py @@ -0,0 +1,194 @@ +############################################################################### +# 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. + + In general, any objects passed to a memoized method should be hashable, + however this will convert any lists or dictionaries passed in to hashable + tuples to store their values in the cache. + + 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 = {} + + key = (self.func, pargs[1:], frozenset(kwargs)) + + try: + value = cache[key] + except KeyError: + value = cache[key] = 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_ From 683154d704ef5af5f896ee0903523c7eb13351ba Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 12:25:37 -0700 Subject: [PATCH 02/17] boot-sp PR-2: smoothed bootstrap/bagging engine and driver dispatch Add smoothed_boot_sp.py (smoothed bootstrap and bagging that fit a statdist univariate distribution to the sampled data and resample from it) with compute_smoothed_ci as the single smoothed-dispatch point, the counterpart to boot_sp.compute_ci for the empirical methods. Wire the drivers to route smoothed vs empirical by boot_method: user_boot returns the empirical 6-tuple or the smoothed (ci_gap, center_gap) pair; simulate_boot gains smoothed_main_routine. This fixes the design-doc section 4.3 latent bug: boot-sp's smoothed simulation called the undefined fit_resample_utils.compute_xhat; it now calls boot_utils.compute_xhat. The PR-1 "smoothed not yet merged" error is removed. Non-root ranks return a matching-arity (None, None) so callers can unpack safely under mpiexec. The only scipy use (norm.ppf) is replaced by statistics.NormalDist().inv_cdf, matching boot_sp.py. Co-Authored-By: Claude Opus 4.8 --- .../confidence_intervals/bootsp/__init__.py | 4 +- .../confidence_intervals/bootsp/boot_sp.py | 12 +- .../confidence_intervals/bootsp/boot_utils.py | 19 +- .../bootsp/simulate_boot.py | 71 +++- .../bootsp/smoothed_boot_sp.py | 308 ++++++++++++++++++ .../confidence_intervals/bootsp/user_boot.py | 59 +++- 6 files changed, 432 insertions(+), 41 deletions(-) create mode 100644 mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py 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..2b28a1a3b 100644 --- a/mpisppy/confidence_intervals/bootsp/boot_sp.py +++ b/mpisppy/confidence_intervals/bootsp/boot_sp.py @@ -684,14 +684,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..4ab51c4c1 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,77 @@ 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 + seed_list = [i * cfg.nB * 100 + 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..f90ded0ad --- /dev/null +++ b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py @@ -0,0 +1,308 @@ +############################################################################### +# 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, mpicomm): + """ Estimate the CI center (the optimality gap) from the fitted distribution. """ + 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). """ + 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) + + boot_cfg = cfg() # for ephemeral changes to deal with seed_offset + boot_cfg.use_fitted = True + + for iter in range(local_nB): + # seed_offset makes unique samples + if serial: + seed_offset = iter + else: + seed_offset = sum(boot_sp.slice_lens(boot_cfg.nB)[:my_rank]) + iter + boot_cfg.seed_offset = cfg.seed_offset + seed_offset + + scenario_pool = list(range(boot_cfg.seed_offset, + boot_cfg.seed_offset + cfg.subsample_size)) + + 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): + """ use the original data to estimate the center, then perform a smoothed estimation of width of confidence intervals + 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) + + # estimation of CI center + dag_gap = center_smoothed(cfg, module, xhat, mpicomm=comm) + comm.Barrier() + cfg.use_fitted = True + + # conduct an m out of n bootstrap, with B = cfg.nB + 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 + + assert cfg.smoothed_B_I is not None, "B_I required for smoothed bagging" + + 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) + + s1 = np.var(avg_gaps) + s2 = np.var(all_gaps) + 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/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__': From cf574255f85b95cdb2b07d5c38ceda65adc96aa4 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 12:25:49 -0700 Subject: [PATCH 03/17] boot-sp PR-2: farmer, cvar, and multi_knapsack examples Add the three examples that build their scenario data with statdist and so could not ship in the statdist-free PR-1: farmer (crop yields perturbed by a univariate distribution), cvar (Lam & Qian, standard normal data), and multi_knapsack (Vaagen & Wallace, deterministic data from a json file). Each is self-contained: a fixed-name xhat_generator builds the EF from the module's own scenario_creator (no external amalgamator coupling), with an empirical json/bash and a smoothed json. multi_knapsack's data_sampler reads its deterministic data from the file when the smoothed driver has not stashed it, so the empirical path works too, and resolves the json relative to the module directory. Co-Authored-By: Claude Opus 4.8 --- examples/bootsp/cvar/cvar.bash | 19 + examples/bootsp/cvar/cvar.json | 16 + examples/bootsp/cvar/cvar.py | 174 ++++++++ examples/bootsp/cvar/smoothed_cvar.json | 18 + examples/bootsp/farmer/farmer.bash | 27 ++ examples/bootsp/farmer/farmer.json | 19 + examples/bootsp/farmer/farmer.py | 399 ++++++++++++++++++ examples/bootsp/farmer/smoothed_farmer.json | 21 + .../bootsp/multi_knapsack/multi_knapsack.bash | 22 + .../bootsp/multi_knapsack/multi_knapsack.json | 17 + .../bootsp/multi_knapsack/multi_knapsack.py | 234 ++++++++++ .../multi_knapsack/multi_knapsack_data.json | 85 ++++ .../smoothed_multi_knapsack.json | 19 + 13 files changed, 1070 insertions(+) create mode 100644 examples/bootsp/cvar/cvar.bash create mode 100644 examples/bootsp/cvar/cvar.json create mode 100644 examples/bootsp/cvar/cvar.py create mode 100644 examples/bootsp/cvar/smoothed_cvar.json create mode 100644 examples/bootsp/farmer/farmer.bash create mode 100644 examples/bootsp/farmer/farmer.json create mode 100644 examples/bootsp/farmer/farmer.py create mode 100644 examples/bootsp/farmer/smoothed_farmer.json create mode 100644 examples/bootsp/multi_knapsack/multi_knapsack.bash create mode 100644 examples/bootsp/multi_knapsack/multi_knapsack.json create mode 100644 examples/bootsp/multi_knapsack/multi_knapsack.py create mode 100644 examples/bootsp/multi_knapsack/multi_knapsack_data.json create mode 100644 examples/bootsp/multi_knapsack/smoothed_multi_knapsack.json 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" +} From f61dd6b59ac78c9437d65e5b0600135a4da4d998 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 12:26:10 -0700 Subject: [PATCH 04/17] boot-sp PR-2: smoothed tests, CI/coverage wiring, run_all examples Add test_boot_sp_smoothed.py: direct tests of the statdist univariate distributions, the empirical farmer/cvar tests that PR-1 could not host (they need statdist), and the smoothed methods (kernel/bagging serial and under mpiexec -np 2; epi-spline tests skip without a nonlinear solver). Wire it into the confidence-intervals CI job and run_coverage.bash, serial and np=2, in the same commit. Update the two PR-1 test files whose "smoothed not yet merged" assertions no longer hold. Register a farmer (empirical) and a cvar (smoothed) run in part 1 of run_all.py. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/test_pr_and_main.yml | 2 + examples/run_all.py | 10 + mpisppy/tests/test_boot_sp.py | 13 +- mpisppy/tests/test_boot_sp_simulate.py | 7 - mpisppy/tests/test_boot_sp_smoothed.py | 335 +++++++++++++++++++++++++ run_coverage.bash | 6 + 6 files changed, 357 insertions(+), 16 deletions(-) create mode 100644 mpisppy/tests/test_boot_sp_smoothed.py diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index dc3fb82b0..04717a28a 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -863,12 +863,14 @@ jobs: cd mpisppy/tests coverage run $COV_ARGS test_boot_sp.py coverage run $COV_ARGS test_boot_sp_simulate.py + coverage run $COV_ARGS test_boot_sp_smoothed.py - name: run bootstrap CI tests (mpiexec -np 2) timeout-minutes: 10 run: | cd mpisppy/tests mpiexec -np 2 coverage run $COV_ARGS -m mpi4py test_boot_sp_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/examples/run_all.py b/examples/run_all.py index ecfb49bfa..e90198e48 100644 --- a/examples/run_all.py +++ b/examples/run_all.py @@ -154,6 +154,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/tests/test_boot_sp.py b/mpisppy/tests/test_boot_sp.py index 23a945435..1a2ffceed 100644 --- a/mpisppy/tests/test_boot_sp.py +++ b/mpisppy/tests/test_boot_sp.py @@ -205,9 +205,11 @@ 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_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()) @@ -257,13 +259,6 @@ def test_user_boot_main_routine(self): self._assert_close_list(res[0], locked_ci_optimal["Classical_quantile"]) self.assertGreaterEqual(res[2][0], 0.0) # ci_gap[0] clamped to >= 0 - @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..da25fc5c1 --- /dev/null +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -0,0 +1,335 @@ +############################################################################### +# 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 warnings +import unittest + +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 +from mpisppy.confidence_intervals.bootsp.statdist.distribution_factory import ( + distribution_factory, +) + +sputils.disable_tictoc_output() + +# statdist integrates array-valued pdfs through scipy.integrate.quad, which +# emits a NumPy>=1.25 "array to scalar" DeprecationWarning many thousands of +# times; silence just that one so the (large) CI logs stay readable. The +# numerics are unchanged. +warnings.filterwarnings( + "ignore", + message="Conversion of an array with ndim > 0 to a scalar is deprecated", + category=DeprecationWarning, +) + +solver_available, solver_name, persistent_available, persistent_solver_name = get_solver() +ipopt_available = pyo.SolverFactory("ipopt").available(exception_flag=False) + +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_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 + import importlib + import mpisppy.confidence_intervals.bootsp.statdist.distributions as dmod + importlib.reload # (noop reference; module already imported) + self.assertTrue(hasattr(dmod, "UnivariateGaussianKernelDistribution")) + + 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_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_kernel_fit_inverse(self): + # the kernel-density fit backs Smoothed_boot_kernel and Smoothed_bagging + import numpy as np + 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_empirical_fit_inverse(self): + import numpy as np + 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)) + + @unittest.skipIf(not ipopt_available, "ipopt (nonlinear solver) not available") + def test_epispline_fit_inverse(self): + import numpy as np + 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 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"]: + res = boot_sp.compute_ci(_make_farmer_cfg(method), module, xhat) + self.assertEqual(len(res), 6) + for ci in res[:3]: + self.assertLessEqual(ci[0], ci[1], msg=f"{method}: {ci}") + + @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) + for ci in res[:3]: + self.assertLessEqual(ci[0], ci[1]) + + @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) + r1 = boot_sp.compute_ci(_make_cvar_cfg("Classical_gaussian"), module, xhat) + r2 = boot_sp.compute_ci(_make_cvar_cfg("Classical_gaussian"), module, xhat) + 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 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 199858c32..fffae4e7e 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 @@ -246,6 +249,9 @@ run_phase "test_with_cylinders (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 From 980b261c599fef4f0681070f6c50cccecfebf648 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 12:26:20 -0700 Subject: [PATCH 05/17] boot-sp PR-2: document the smoothed methods and statdist Replace the "smoothed methods merged separately" note with a description of the empirical vs smoothed families, add the five Smoothed_* tokens to the methods table, document data_sampler and the smoothed-only options, and add a section on the smoothed methods, the bundled statdist library, and the farmer/cvar/multi_knapsack examples. Co-Authored-By: Claude Opus 4.8 --- doc/src/boot_sp.rst | 81 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 75 insertions(+), 6 deletions(-) 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 ---------- From e9c3d56bd0b7fb2311f8259cc5b7dd857d2f8f61 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 13:17:05 -0700 Subject: [PATCH 06/17] boot-sp PR-2: guard np=2 rank assertions in the empirical example tests The Test_empirical_examples tests inspected compute_ci's result on every rank, but only rank 0 gets the values (non-root ranks get a tuple of Nones). Subscripting None raised on the non-root rank and aborted the test mid-loop, so that rank skipped the next compute_ci while rank 0 entered its collective and deadlocked. Guard the value checks on rank 0 (and assert the None-tuple off-rank) so both ranks always make the same collective calls. Found by running the suite under mpiexec -np 2. Co-Authored-By: Claude Opus 4.8 --- mpisppy/tests/test_boot_sp_smoothed.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index da25fc5c1..145de3e9e 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -192,10 +192,14 @@ def test_farmer_empirical_wellformed(self): 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) - for ci in res[:3]: - self.assertLessEqual(ci[0], ci[1], msg=f"{method}: {ci}") + 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): @@ -205,18 +209,23 @@ def test_cvar_empirical_wellformed(self): self.assertIn("ROOT", xhat) res = boot_sp.compute_ci(_make_cvar_cfg("Classical_quantile"), module, xhat) self.assertEqual(len(res), 6) - for ci in res[:3]: - self.assertLessEqual(ci[0], ci[1]) + 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) - for a, b in zip(list(r1[0]), list(r2[0])): - self.assertEqual(round_pos_sig(a, 6), round_pos_sig(b, 6)) + 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)) #***************************************************************************** From 40e6fb59ef27c8a46e8c949c667e19b91073aaa8 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Fri, 3 Jul 2026 13:40:54 -0700 Subject: [PATCH 07/17] boot-sp PR-2: return a scalar from the kernel pdf (fix scipy quad warning) The Gaussian-kernel distribution's pdf returned gaussian_kde.evaluate(x), a size-1 numpy array. The base-class cdf feeds pdf to scipy.integrate.quad, which converts the integrand to a scalar and, on NumPy>=1.25, emitted an "array to scalar" DeprecationWarning tens of thousands of times during a smoothed run. Return float(...[0]) so the integrand is a scalar; the numerics are unchanged. This removes the need for the warnings filter that was in test_boot_sp_smoothed.py. Co-Authored-By: Claude Opus 4.8 --- .../bootsp/statdist/distributions.py | 5 ++++- mpisppy/tests/test_boot_sp_smoothed.py | 11 ----------- 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py index 1c5350504..aa96000df 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -311,7 +311,10 @@ def pdf(self, x): (float) The value of the probability density function of this distribution on x. """ - return self.kernel.evaluate(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): """ diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index 145de3e9e..4867bede0 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -22,7 +22,6 @@ import os import sys import math -import warnings import unittest import pyomo.environ as pyo @@ -40,16 +39,6 @@ sputils.disable_tictoc_output() -# statdist integrates array-valued pdfs through scipy.integrate.quad, which -# emits a NumPy>=1.25 "array to scalar" DeprecationWarning many thousands of -# times; silence just that one so the (large) CI logs stay readable. The -# numerics are unchanged. -warnings.filterwarnings( - "ignore", - message="Conversion of an array with ndim > 0 to a scalar is deprecated", - category=DeprecationWarning, -) - solver_available, solver_name, persistent_available, persistent_solver_name = get_solver() ipopt_available = pyo.SolverFactory("ipopt").available(exception_flag=False) From ad06446a9aa3df4393ecb561b9aa7bcd474a810f Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 10:49:21 -0700 Subject: [PATCH 08/17] boot-sp PR-2: draw the smoothed bootstrap batches and center correctly Two defects in the smoothed bootstrap, both measured against Chen & Woodruff (2024), "Distributions and Bootstrap for Data-based Stochastic Programming". They are fixed together because the second is what forces the index-space reservation the first introduces. 1. 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 rather than independent resamples. Algorithm 3 draws a fresh set of N points from the fitted distribution for each of the B batches. Stride by the batch size, as smoothed_bagging already did, so the blocks are pairwise disjoint. The estimated spread was badly understated: on the cvar example (N=20, nB=10) the interval widens about fivefold. The dead boot_cfg copy goes with it -- it was mutated but never passed to the solves, which is what hid the missing stride. 2. 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. Set use_fitted before estimating the center (smoothed_bagging already did), and reserve the center's block ahead of the batch blocks, since both now sample the same distribution and must not reuse each other's draws. simulate_boot's coverage replications are spaced to match. The ported spacing was a hard-coded nB * 100, unrelated to how many record numbers a replication actually consumes; it is now exactly that footprint. The empirical harness needs no spacing at all, because it gets independent numpy streams from default_rng([seed_offset, word]); the smoothed path addresses draws by record number, since that is what the model seeds each draw with, so its replications are separated by giving each a disjoint block. The new test captures the actual scenario pools and asserts the batch blocks are pairwise disjoint and disjoint from the center's, and that use_fitted holds when the center is computed; it fails on the old code for both reasons. No locked values change: the smoothed tests assert only structure. Co-Authored-By: Claude Opus 5 --- doc/designs/bootsp_merge_design.md | 27 ++++++++++ .../bootsp/simulate_boot.py | 13 ++++- .../bootsp/smoothed_boot_sp.py | 45 +++++++++++------ mpisppy/tests/test_boot_sp_smoothed.py | 49 +++++++++++++++++++ 4 files changed, 117 insertions(+), 17 deletions(-) diff --git a/doc/designs/bootsp_merge_design.md b/doc/designs/bootsp_merge_design.md index 0c7cfa0a1..6432052e4 100644 --- a/doc/designs/bootsp_merge_design.md +++ b/doc/designs/bootsp_merge_design.md @@ -286,6 +286,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/mpisppy/confidence_intervals/bootsp/simulate_boot.py b/mpisppy/confidence_intervals/bootsp/simulate_boot.py index 4ab51c4c1..b73533209 100644 --- a/mpisppy/confidence_intervals/bootsp/simulate_boot.py +++ b/mpisppy/confidence_intervals/bootsp/simulate_boot.py @@ -107,7 +107,18 @@ def smoothed_main_routine(cfg, module): ci_len = [] run_time = [] seed_offset = cfg.seed_offset # store the original offset - seed_list = [i * cfg.nB * 100 + seed_offset for i in range(cfg.coverage_replications)] + # 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 diff --git a/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py index f90ded0ad..7e9ef266a 100644 --- a/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py +++ b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py @@ -71,7 +71,18 @@ def center_smoothed(cfg, module, xhat, mpicomm): 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). """ + 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: @@ -79,19 +90,14 @@ def smoothed_resample_helper(cfg, module, xhat, serial=False): local_boot_gaps = np.empty(local_nB, dtype=np.float64) - boot_cfg = cfg() # for ephemeral changes to deal with seed_offset - boot_cfg.use_fitted = True + 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): - # seed_offset makes unique samples - if serial: - seed_offset = iter - else: - seed_offset = sum(boot_sp.slice_lens(boot_cfg.nB)[:my_rank]) + iter - boot_cfg.seed_offset = cfg.seed_offset + seed_offset - - scenario_pool = list(range(boot_cfg.seed_offset, - boot_cfg.seed_offset + cfg.subsample_size)) + 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) @@ -102,7 +108,8 @@ def smoothed_resample_helper(cfg, module, xhat, serial=False): def smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline', quantile=False, serial=False): - """ use the original data to estimate the center, then perform a smoothed estimation of width of confidence intervals + """ 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 @@ -121,13 +128,19 @@ def smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline', qua 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 - # estimation of CI center + # 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, mpicomm=comm) comm.Barrier() - cfg.use_fitted = True - # conduct an m out of n bootstrap, with B = cfg.nB + # 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() diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index 4867bede0..e02f3e8f8 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -284,6 +284,55 @@ def test_simulate_smoothed_coverage(self): 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 + + def spy_eval(cfg_, module_, scenarios, xhat_, duplication=True, mpicomm=None): + pools.append(list(scenarios)) + return real_eval(cfg_, module_, scenarios, xhat_, + duplication=duplication, mpicomm=mpicomm) + + def spy_center(cfg_, module_, xhat_, mpicomm): + fitted_at_center.append(cfg_.use_fitted) + return real_center(cfg_, module_, xhat_, mpicomm) + + 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): From 7a82a8ef17bc9c19020ce017309d3d240d17d765 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 11:50:56 -0700 Subject: [PATCH 09/17] boot-sp PR-2: sample variances in smoothed bagging, and a dead parameter np.var defaults to ddof=0, but the algorithm specifies sample variances and the empirical estimators already use ddof=1 for their Gaussian half-width. ddof=0 understates s1 by (B_I-1)/B_I -- a third at the B_I=3 used in run_all.py -- so the bagging interval came out too narrow. Since s1 is the variance among the B_I per-seed-point averages, B_I < 2 leaves it undefined and was silently producing a nan; that is now a clear error. center_smoothed took an mpicomm it never used; dropped, along with the argument at its one call site. run_all.py's do_one_boot still described its examples as statdist-free with the others "merged separately"; it registers farmer and cvar now. Co-Authored-By: Claude Opus 5 --- examples/run_all.py | 5 ++-- .../bootsp/smoothed_boot_sp.py | 27 ++++++++++++++----- mpisppy/tests/test_boot_sp_smoothed.py | 11 +++++--- 3 files changed, 31 insertions(+), 12 deletions(-) diff --git a/examples/run_all.py b/examples/run_all.py index e90198e48..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}") diff --git a/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py index 7e9ef266a..43fad1daa 100644 --- a/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py +++ b/mpisppy/confidence_intervals/bootsp/smoothed_boot_sp.py @@ -51,8 +51,13 @@ def fit_distribution(sample_data, distr_type='univariate-epispline'): return fitted_distr -def center_smoothed(cfg, module, xhat, mpicomm): - """ Estimate the CI center (the optimality gap) from the fitted distribution. """ +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, @@ -136,7 +141,7 @@ def smoothed_bootstrap(cfg, module, xhat, distr_type='univariate-epispline', qua # 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, mpicomm=comm) + dag_gap = center_smoothed(cfg, module, xhat) comm.Barrier() # each batch is a fresh set of cfg.sample_size draws from the same fitted @@ -206,7 +211,14 @@ def smoothed_bagging(cfg, module, xhat, distr_type='univariate-kernel', serial=F all_gaps = None avg_gaps = None - assert cfg.smoothed_B_I is not None, "B_I required for smoothed bagging" + # 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): @@ -234,8 +246,11 @@ def smoothed_bagging(cfg, module, xhat, distr_type='univariate-kernel', serial=F dag_gap = np.mean(avg_gaps) - s1 = np.var(avg_gaps) - s2 = np.var(all_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 diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index e02f3e8f8..ba96e1b2e 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -303,14 +303,17 @@ def test_smoothed_bootstrap_draws_are_disjoint_and_fitted(self): real_eval = boot_sp.evaluate_scenarios real_center = smoothed_boot_sp.center_smoothed - def spy_eval(cfg_, module_, scenarios, xhat_, duplication=True, mpicomm=None): + # 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, mpicomm=mpicomm) + duplication=duplication) - def spy_center(cfg_, module_, xhat_, mpicomm): + def spy_center(cfg_, module_, xhat_): fitted_at_center.append(cfg_.use_fitted) - return real_center(cfg_, module_, xhat_, mpicomm) + return real_center(cfg_, module_, xhat_) boot_sp.evaluate_scenarios = spy_eval smoothed_boot_sp.center_smoothed = spy_center From fcea64ae1059baa34655f8639ce5e43ea5fada73 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 11:53:57 -0700 Subject: [PATCH 10/17] boot-sp PR-2: refuse a maximization model instead of reporting [0, 0] The estimators form every gap as (value at xhat) minus the optimal. For a maximization the value at xhat is a lower bound on the optimal, so that quantity is non-positive -- and user_boot floors the reported interval with ci_gap[0] = max(0, ci_gap[0]). A maximization run would therefore not merely be wrong, it would report [0, 0]. Per the repo-wide rule that maximization either works or raises, this is the raise, checked in solve_routine, which every extensive form goes through. The message says what supporting maximization would take rather than just refusing. The test builds one trivial model twice, once each sense, and checks that the maximization raises while the minimization gets past the guard and fails later on a deliberately bogus solver name, so the guard is not what stopped it. Co-Authored-By: Claude Opus 5 --- .../confidence_intervals/bootsp/boot_sp.py | 23 ++++++++++++ mpisppy/tests/test_boot_sp.py | 36 +++++++++++++++++++ 2 files changed, 59 insertions(+) diff --git a/mpisppy/confidence_intervals/bootsp/boot_sp.py b/mpisppy/confidence_intervals/bootsp/boot_sp.py index 2b28a1a3b..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 diff --git a/mpisppy/tests/test_boot_sp.py b/mpisppy/tests/test_boot_sp.py index 1a2ffceed..f648611ee 100644 --- a/mpisppy/tests/test_boot_sp.py +++ b/mpisppy/tests/test_boot_sp.py @@ -205,6 +205,42 @@ def test_compute_xhat_requires_generator(self): self.assertIn("xhat_generator", msg) self.assertIn("xhat_generator_no_generator_module", msg) + 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 From fe0dff6e3c70136daec99b06a77fa1781b41086c Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 15:40:15 -0700 Subject: [PATCH 11/17] boot-sp PR-2: fix eight statdist defects the new unit tests expose Writing direct tests for the trimmed statdist univariate distributions turned up eight things that were wrong in the ported code: - Parameter.instantiated was inverted: it was True exactly when the parameter had no value. set_value now keeps it up to date too. - memoize_method left the keyword *values* out of the cache key, so cdf(x, epsabs=1e-4) and cdf(x, epsabs=1e-9) shared one answer. - memoize_method also died with "unhashable type: 'list'" when a method was handed a list or a dict, hiding the method's own behavior (region_expectation and region_probability validate their argument and say what is wrong). Such a call cannot be cached, so it is now passed through; the docstring, which claimed the arguments were converted, is corrected to match. - UnivariateDiscrete.var had the sign backwards (mean**2 - E[X**2]), so the variance of every non-degenerate discrete distribution came out negative. - UnivariateDiscrete accepted breakpoints whose probabilities summed to less than one: the check was one-sided. - UnivariateDiscrete.cdf_inverse's error message named cdf instead. - UnivariateGaussianKernelDistribution.fit dropped the bw_method it was given and always passed None, so the caller could not widen the bandwidth. - The kernel distribution's draw method passed [n, seed] to gaussian_kde.resample as the sample size, which raises; the size and the seed are two arguments. It is also named generates_X now, like its siblings. - UnivariateEmpiricalDistribution.cdf and cdf_inverse raised IndexError on a constant sample (a one-record sample included) while looking for a second point to extrapolate along. Every quantile of a constant sample is that value, and the cdf steps from 0 to 1 there. Co-Authored-By: Claude Opus 5 --- .../bootsp/statdist/base_distribution.py | 5 ++- .../bootsp/statdist/distributions.py | 36 +++++++++++++------ .../bootsp/statdist/utilities.py | 14 +++++--- 3 files changed, 40 insertions(+), 15 deletions(-) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py index ca3dc001d..1e7fd375d 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py @@ -50,7 +50,9 @@ def __init__(self, name, value=None, bounds=(None, None), kind=float): """ self.name = name self.value = value - self.instantiated = value is None + # 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 @@ -61,6 +63,7 @@ def set_value(self, value): 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) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py index aa96000df..5fbdb313b 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -300,7 +300,7 @@ def fit(cls, data, bw_method=None): Returns: UnivariateEmpiricalDistribution: The fitted distribution """ - return UnivariateGaussianKernelDistribution(data, bw_method=None) + return UnivariateGaussianKernelDistribution(data, bw_method=bw_method) def pdf(self, x): """ @@ -328,8 +328,10 @@ def _cdf(self, x): #use the cdf_inverse function in base_distribution - def generate_X(self, n=1, seed=0): - return self.kernel.resample([n, seed]) + 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) @@ -585,6 +587,11 @@ def cdf(self, x, lower_bound=None, upper_bound=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) @@ -605,6 +612,10 @@ def cdf(self, x, lower_bound=None, upper_bound=None): 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 @@ -654,6 +665,11 @@ def cdf_inverse(self, x, lower_bound=None, upper_bound=None): 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( @@ -682,13 +698,11 @@ def cdf_inverse(self, x, lower_bound=None, upper_bound=None): # (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 - if j - 1 == -len(self.input_data): - print("Warning: all input values are the same (", - self.input_data[j], ")") - return self.input_data[j] # 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] @@ -787,8 +801,10 @@ def __init__(self, breakpoints): if val < lastval: raise RuntimeError("DiscreteDistribution dict must be ordered by val:"+str(val)+" < "+str(lastval)) lastval = val - self.var = self.mean*self.mean - Esqsum - if sumprob - 1 > tol: # could use gosm_options.cdf_tolerance + 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__() @@ -825,7 +841,7 @@ 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 called for a discrete distribution.") + raise RuntimeError("cdf_inverse called for a discrete distribution.") def sample_one(self): """ diff --git a/mpisppy/confidence_intervals/bootsp/statdist/utilities.py b/mpisppy/confidence_intervals/bootsp/statdist/utilities.py index 7e64f6a4f..a7c6d30ac 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/utilities.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/utilities.py @@ -94,9 +94,9 @@ class memoize_method: 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. - In general, any objects passed to a memoized method should be hashable, - however this will convert any lists or dictionaries passed in to hashable - tuples to store their values in the cache. + 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 @@ -150,12 +150,18 @@ def __call__(self, *pargs, **kwargs): else: cache = obj.__cache = {} - key = (self.func, pargs[1:], frozenset(kwargs)) + # 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 From f94657fd70a91724635db747b5441f78a80ff5ef Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 15:40:26 -0700 Subject: [PATCH 12/17] boot-sp PR-2: unit test the statdist univariate distributions statdist arrived with the smoothed methods as its only exercise, which left it at 27% coverage. These are direct tests, all solver-free, of what the library promises: the plotting-position contract of the empirical cdf and its inverse, the closed forms for the uniform, normal and Student's t, the kernel density's padded domain and scalar pdf, the discrete distribution's moments and step function, and the registry's naming rules. The generic machinery in base_distribution.py (numeric cdf by quadrature, cdf inversion by refinement, expectations, sampling) is tested against a local distribution with a closed form for all of it -- density 2x on [0, 1] -- so the assertions are exact rather than self-referential. utilities.py gets tests for both memoizers and the argv context manager. Coverage of the reachable statdist code: distributions.py 37% -> 87%, utilities.py 48% -> 98%, base_distribution.py 24% -> 42% (only three univariate statements are left uncovered there: two abstract method bodies and the interactive plt.show(); the rest of the miss is the MultivariateDistribution class, which the trim left behind and which nothing can reach). Everything lands in the existing test file, so no harness wiring changes. Two tests are skipped in the usual environments: the plot smoke test needs matplotlib, and the check that importing statdist does not import scipy needs to spawn a plain python subprocess, which cannot be done from inside an mpiexec launch. Co-Authored-By: Claude Opus 5 --- mpisppy/tests/test_boot_sp_smoothed.py | 573 ++++++++++++++++++++++++- 1 file changed, 565 insertions(+), 8 deletions(-) diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index ba96e1b2e..6e0caa497 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -22,8 +22,13 @@ import os import sys import math +import importlib.util +import subprocess +import tempfile import unittest +from collections import OrderedDict +import numpy as np import pyomo.environ as pyo import mpisppy.utils.sputils as sputils from mpisppy.tests.utils import get_solver, round_pos_sig @@ -33,6 +38,12 @@ 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, ) @@ -41,6 +52,8 @@ 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 @@ -124,13 +137,29 @@ def test_factory_drops_multivariate(self): 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 - import importlib - import mpisppy.confidence_intervals.bootsp.statdist.distributions as dmod - importlib.reload # (noop reference; module already imported) - self.assertTrue(hasattr(dmod, "UnivariateGaussianKernelDistribution")) + # 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) @@ -138,14 +167,98 @@ def test_uniform_inverse(self): 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_uses_the_documented_df(self): + data = list(np.random.RandomState(4).normal(0.0, 2.0, size=500)) + var = float(np.var(data)) + self.assertGreater(var, 1.0) + st = distribution_factory("univariate-student").fit(data) + self.assertAlmostEqual(st.mean, float(np.mean(data))) + self.assertAlmostEqual(st.var, var) + self.assertAlmostEqual(st.df, 2*var/(var - 1)) + + def test_student_fit_falls_back_to_one_df(self): + # 2v/(v-1) is not a usable number of degrees of freedom once v <= 1 + data = list(np.random.RandomState(5).normal(0.0, 0.1, size=200)) + self.assertLessEqual(float(np.var(data)), 1.0) + st = distribution_factory("univariate-student").fit(data) + self.assertEqual(st.df, 1) + + 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 location and scale as a normal, but with the fatter tails + norm = distribution_factory("univariate-normal")(var=4.0, mean=1.0) + 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 - import numpy as np data = list(np.random.RandomState(0).normal(0, 1, size=200)) kde = distribution_factory("univariate-kernel").fit(data) lo = kde.cdf_inverse(0.25) @@ -153,20 +266,464 @@ def test_kernel_fit_inverse(self): 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): - import numpy as np 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): - import numpy as np 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). From e2e9ba9666bc0c78473e16a6f1cd8743d34930c3 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 15:56:41 -0700 Subject: [PATCH 13/17] boot-sp PR-2: finish the univariate trim of statdist The port dropped copula.py, vine.py, bicop.py and the multivariate distribution classes in distributions.py, but MultivariateDistribution itself stayed behind in base_distribution.py along with the five decorators that exist only to serve it (fit_wrapper, accepts_dict, returns_dict, params_as_args, params_as_args2). With no multivariate subclasses left, none of it is reachable: nothing in mpisppy, examples or doc names any of those symbols. So this removes 596 lines that statdist/README.md already said were not here, and the README now says so about base_distribution.py too. The univariate half of the file goes from 42% to 98% covered (what is left is the two abstract method bodies and the interactive plt.show()). Co-Authored-By: Claude Opus 5 --- .../bootsp/statdist/README.md | 6 +- .../bootsp/statdist/base_distribution.py | 597 ------------------ 2 files changed, 4 insertions(+), 599 deletions(-) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/README.md b/mpisppy/confidence_intervals/bootsp/statdist/README.md index 4f5758798..2c2146554 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/README.md +++ b/mpisppy/confidence_intervals/bootsp/statdist/README.md @@ -17,8 +17,10 @@ Only the **univariate** distributions and their support modules: ## What was dropped -The multivariate machinery — `copula.py`, `vine.py`, `bicop.py`, and the -multivariate distribution classes in `distributions.py` — is **not** included. +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 diff --git a/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py index 1e7fd375d..9f9a6708d 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/base_distribution.py @@ -10,7 +10,6 @@ This abstract base class is the parent class of all distribution classes. """ from abc import ABCMeta, abstractmethod -from functools import wraps import os import numpy as np @@ -504,599 +503,3 @@ def log_likelihood(self, data): """ return sum(np.log(self.pdf(x)) for x in data) - -class MultivariateDistribution(BaseDistribution): - """ - This class is an abstract base class for all multivariate distributions - TODO: This docstring should be improved greatly!! - """ - __metaclass__ = ABCMeta - - def __init__(self, dimension, dimkeys=None, parameters=None, lower=None, - upper=None, bounds=None): - """ - Args: - dimension (int): The dimension of the distribution - dimkeys (List): A list of the names of the dimensions, by default, - these will just be the indices. If passed in, this will enable - you to refer to values by the dimension name in certain - functions - parameters (list[Parameter]): A list of parameters for the - distribution - lower (list[float]): A list of the lower bounds of the support - of the distribution - upper (list[float]): A list of the upper bounds of the support - of the distribution - bounds (list|dict): A colection of bounds on the support for - each dimension. We assume the support is on a rectangular - region. If it is passed as a dictionary it should map - dimension names to ordered pairs of lower and upper bounds. A - None indicates that there is no lower or upper bound for a - given dimension. - """ - BaseDistribution.__init__(self, dimension, parameters) - self.ndim = dimension - if dimkeys is None: - # We default to using the integers if no dimkeys are passed in. - self.dimkeys = list(range(self.ndim)) - else: - self.dimkeys = dimkeys - - if lower is None: - self.lower = [None for _ in range(dimension)] - else: - self.lower = lower - if upper is None: - self.upper = [None for _ in range(dimension)] - else: - self.upper = upper - - if bounds is None: - self.bounds = [(-np.inf, np.inf)] * dimension - elif isinstance(bounds, list): - self.bounds = bounds - elif isinstance(bounds, dict): - self.bounds = [bounds[dim] for dim in self.dimkeys] - - def pdf(self, *xs): - raise NotImplementedError - - def log_likelihood(self, data): - """ - This method will return the log likelihood of the observed data - given the fitted model. - - This method just naively computes the pdf and applies the logarithm. - It would be more efficient in subclasses to find an expression for - the log-likelihood. - - The argument data can either be a list of vectors for each dimension - of the data or it can be a dictionary mapping dimension names to the - corresponding vector of data. - - Args: - data (list[list[float]] | dict[list[float]]): The observed values - Returns: - float: The computed log-likelihood - """ - if isinstance(data, dict): - vects = [data[dimkey] for dimkey in dimkeys] - else: - vects = data - - return sum(np.log(self.pdf(*xs)) for xs in zip(*vects)) - - def plot(self, func, lower=None, upper=None): - """ - Args: - func (str): The function to plot, either 'pdf' or 'cdf' - lower (list[float]): A list of the lower bounds for the plot, - will default to the lower bounds of the support if None - upper (list[float]): A list of the upper bounds of the plot - will default to the upper bounds of the support if None - """ - if self.dimension != 2: - raise ValueError("This plot method is only functional for 2-d " - "distributions.") - - if lower is None: - lower = self.lower - if lower[0] is None: - lower[0] = -5 - if lower[1] is None: - lower[1] = 5 - if upper is None: - upper = self.upper - if upper[0] is None: - upper[0] = -5 - if upper[1] is None: - upper[1] = 5 - - import matplotlib.pyplot as plt - fig = plt.figure() - ax = fig.gca(projection='3d') - - X = np.arange(lower[0], upper[0], 0.1) - Y = np.arange(lower[1], upper[1], 0.1) - X, Y = np.meshgrid(X, Y) - - Z = np.zeros_like(X) - for i, row in enumerate(X): - for j, x in enumerate(row): - y = Y[i,j] - if func == 'pdf': - z = self.pdf(x, y) - elif func == 'cdf': - z = self.cdf(x, y) - - Z[i,j] = z - - ax.plot_surface(X, Y, Z) - ax.set_xlim(lower[0], upper[0]) - ax.set_ylim(lower[1], upper[1]) - return ax - - def contour_plot(self, func, lower=None, upper=None): - """ - Args: - func (str): The function to plot, either 'pdf' or 'cdf' - lower (list[float]): A list of the lower bounds for the plot, - will default to the lower bounds of the support if None - upper (list[float]): A list of the upper bounds of the plot - will default to the upper bounds of the support if None - """ - if self.dimension != 2: - raise ValueError("This plot method is only functional for 2-d " - "distributions.") - - if lower is None: - lower = self.lower - if upper is None: - upper = self.upper - - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - - X = np.arange(lower[0], upper[0], 0.1) - Y = np.arange(lower[1], upper[1], 0.1) - X, Y = np.meshgrid(X, Y) - - Z = np.zeros_like(X) - for i, row in enumerate(X): - for j, x in enumerate(row): - y = Y[i,j] - if func == 'pdf': - z = self.pdf(x, y) - elif func == 'cdf': - z = self.cdf(x, y) - - Z[i,j] = z - - ax.contour(X, Y, Z) - ax.set_xlim(lower[0], upper[0]) - ax.set_ylim(lower[1], upper[1]) - return ax - - @memoize_method - def rect_prob(self, lowerdict, upperdict): - tempdict = dict.fromkeys(self.dimkeys) - def f(n): - - # recursive function that will calculate the cdf - # It has a structure of binary tree - if n == 0: - return self.cdf(tempdict) - else: - tempdict[self.dimkeys[n - 1]] = upperdict[self.dimkeys[n - 1]] - leftresult = f(n - 1) - tempdict[self.dimkeys[n - 1]] = lowerdict[self.dimkeys[n - 1]] - rightresult = f(n - 1) - return leftresult - rightresult - - return f(self.dimension) - - def marginal(self, ys, bounds = None, error_tolerance=None): - """ - This function will evaluate the marginal distribution of the joint - cdf which is composed of the dimensions passed in through ys. - It will evaluate it at the point ys. - - Args: - ys (dict): A dictionary mapping dimension names to their - corresponding values - error_tolerance (int): Value to increase the error tolerance - of the integration process by powers of 10 - Returns: - float: The value of the marginal - """ - - other_dims = [(i, dim) for i, dim in enumerate(self.dimkeys) - if dim not in ys] - - point_dict = ys.copy() - - def pdf_x(*xs): - for (_, dim), x in zip(other_dims, xs): - point_dict[dim] = x - #print(self.pdf(point_dict)) - return self.pdf(point_dict) - if bounds == None: - bounds = [self.bounds[i] for i, _ in other_dims] - - if error_tolerance: - tol = error_tolerance - else: - tol = 0 - - return scipy.integrate.nquad(pdf_x, bounds, opts={'epsabs': (1.49e-08 * (10**tol)), 'epsrel': (1.49e-08 * (10**tol))} )[0] - - def conditional_pdf(self, xs, cond_xs, marginal_cdf=None): - """ - This will evaluate the conditional pdf at the point xs given - that the dimensions in cond_names - - Args: - xs (dict): A dictionary mapping dimension names to values - cond_xs (dict): A dictionary mapping the dimension names - of the conditioned variables to their values - Returns: - float: The value fo the conditional pdf - """ - if marginal_cdf == None: - marg = self.marginal(cond_xs) - else: - marg = marginal_cdf - - point_dict = {} - for dim, x in xs.items(): - point_dict[dim] = x - for dim, x in cond_xs.items(): - point_dict[dim] = x - return self.pdf(point_dict) / marg - - def conditional_cdf(self, xs, cond_xs, marginal_cdf = None): - """ - This will evaluate the conditional cdf at the point xs given - that the dimensions in cond_names are set to the values in cond_xs. - - Args: - xs (dict): A dictionary mapping dimension names to values - cond_xs (dict): A dictionary mapping the dimension names - of the conditioned variables to their values - Returns: - float: The value fo the conditional cdf - """ - bounds = [] - - dimkeys = list(xs.keys()) - - for dim in dimkeys: - dim_index = self.dimkeys.index(dim) - lower_bound = self.bounds[dim_index][0] - bounds.append([lower_bound, xs[dim]]) - - point_dict = cond_xs.copy() - def f(*xs): - for dim, x in zip(dimkeys, xs): - point_dict[dim] = x - return self.pdf(point_dict) - - if marginal_cdf == None: - marg = self.marginal(cond_xs) - else: - marg = marginal_cdf - - try: - return scipy.integrate.nquad(f, bounds)[0] / marg - except: - return 0 - - def conditional_cdf_inverse(self, cond_xs, cdf_value, dim, marginal, - capacity = 4000, n = 100, - method = 'combination', xtol = 0.001): - """ - This function computes the inverse of a conditional cdf value - conditioned on a given point. Therefore 3 different methods are - provided: - - linear interpolation: The conditional cdf is evaluated at several - points in a given interval. Since two points which conditional cdfs - wrap the cdf_value, a linear interpolation between these two - points is used to compute the inverse of the cdf_value. - - bisection: The bisection method from scipy is used to solve the - equation 0 = conditional_cdf - cdf_value. - - combination of both: First a bisection method is used to find the - two wrapping points like in the linear interpolation method. After - that a linear interpolation is used to compute the inverse. - - Args: - cond_xs (dict): A dictionary mapping dimension the dimension - names of the conditioned variables to their values. - cdf_value: The value you want to compute the inverse for. - dim (int or str): The name of the dimension you want to get the - inverse for (e.g. F(X|Y=500) = 0.2: You want to compute the - value of X under the condition that Y=500, so that F equals - 0.2. In that case "dim" equals X.). - marginal (distribution like): The marginal of dimension dim. - capacity (float): The capacity for that day. - n (int): The number of intersection of the interval - [-capacity, capacity], which specify the points which are - evaluated for the linear interpolation method. The number - specifies also a break criteria for the bisection part in the - combined method. - method (str): The method you want to use. "default" refers to the - linear interpolation method, "bisect" to the bisection method - and "combination" to the combined method. - xtol (float): The tolerance for the bisection method. - (break criteria) - - Returns: - The inverse value of the passed in cdf_value conditioned on the - point cond_xs. - """ - marginal_cdf = self.marginal(cond_xs) # This value is needed a lot. So - # it is computed here once. - - if method == 'default': - """ - For the default or linear interpolation method first a list of - points are created. These points are evaluated one after the other - with the conditional_cdf function. After that it is checked, if the - given cdf_value is wrapped by two consecutive points' conditional - cdf. If thats the case, these two points and there conditional cdfs - are used to compute a linear interpolation between them. This - linear interpolation then is used to compute the inverse of the - given cdf_value. Because the conditional_cdf lives in the copula - space (which is [0,1]^n), the points have to be converted to [0,1]. - For the purpose of getting power values as a return, the computed - inverse values have to be transformed back in the end. - """ - points = np.linspace(-capacity, capacity, n) - x = [] - for point in points: - x.append(marginal.cdf(point)) - y = [] - j = 0 - for i in x: - xs = {dim: i} - yi = self.conditional_cdf(xs, cond_xs) - y.append(yi) - if (j==0) and (yi > cdf_value): - inverse = marginal.cdf(-capacity) - break - elif (j != 0): - if y[j-1] <= cdf_value <= y[j]: #linear interpolation - lin = scipy.interpolate.interp1d([y[j-1], y[j]], [x[j-1], x[j]]) - inverse = lin(cdf_value) #computing the inverse - break - j += 1 - else: - inverse = marginal.cdf(capacity) - return marginal.cdf_inverse(inverse) - elif method == 'bisect': - """ - In this method a help function is defined. After that the root - of this function is computed using the bisection mehtod from - scipy. For more information see the scipy documentation. - The transformation of the values is done for the same reason like - above. - """ - def help(d): - dict = {dim: d} - return self.conditional_cdf(dict, cond_xs, - marginal_cdf=marginal_cdf) \ - - cdf_value - if help(marginal.cdf(-capacity)) > 0: - inverse = -capacity - elif help(marginal.cdf(capacity)) < 0: - inverse = capacity - else: - inverse = marginal.cdf_inverse(scipy.optimize.bisect(help, - marginal.cdf(-capacity), - marginal.cdf(capacity), - xtol=xtol)) - - return inverse - - - elif method == 'combination': - """ - In this method not every single point is evaluated. There is some- - thing like a bisection method used to find faster the wrapping - points. After they are found, the linear interpolation is used - to compute the inverse value. - The transformation of the values is done for the same reason like - above. - """ - if capacity is None: - capacity = 0 - l = -capacity - u = capacity - l_cdf = self.conditional_cdf({dim: marginal.cdf(l)}, cond_xs, - marginal_cdf=marginal_cdf) - u_cdf = self.conditional_cdf({dim: marginal.cdf(u)}, cond_xs, - marginal_cdf=marginal_cdf) - if l_cdf >= cdf_value: - print('cdf', cdf_value) - print('lower', l_cdf) - return l - elif u_cdf <= cdf_value: - print('cdf', cdf_value) - print('upper', u_cdf) - return u - tol = (capacity * 2) / n - k = 0 - while ((u - l) > tol) and (k < n): - m = (u + l) / 2 - m_cdf = self.conditional_cdf({dim: marginal.cdf(m)}, cond_xs, - marginal_cdf=marginal_cdf) - if m_cdf < cdf_value: - l = m - l_cdf = m_cdf - elif m_cdf > cdf_value: - u = m - u_cdf = m_cdf - else: - return m - k = k + 1 - lin = scipy.interpolate.interp1d([l_cdf, u_cdf], [marginal.cdf(l), marginal.cdf(u)]) - return marginal.cdf_inverse(lin(cdf_value)) - - -def fit_wrapper(method): - """ - This is a function decorator which will wrap the fit method for - multivariate distributions. It will allow for data to be passed using - a dictionary mapping names to lists of data. - - Internally this transforms the data into a lists of lists and then fits - the distribution to that data. Then it assigns to the dimkeys attribute - the list of names. - - Args: - method: The class method fit of a multivariate distribution - Returns: - method: The modified method to handle dictionaries of input data - """ - @wraps(method) - def fit(cls, data, dimkeys=None, **kwargs): - """ - This function converts the dictionary into a list, passes it to the - fit method and then assigns to the distribution the dimkeys attribute. - """ - vectors = [] - if isinstance(data, dict): - for key in dimkeys: - vectors.append(data[key]) - else: - vectors = data - - distribution = method(cls, vectors, dimkeys, **kwargs) - return distribution - - return fit - - -def accepts_dict(method): - """ - This function decorator will allow any of the methods which accept separate - values for each dimension to also accept a dictionary which has keys - mapping to each dimension. - - For example, the pdf for any distribution generally has the prototype - def pdf(self, *x): - This decorator will unpack the dictionary into its respective dimensions - and pass it to the function. - - The function that this decorator is applied to must have a prototype like - def f(self, *x) - - This will enable you to call a function in the following three ways. - - Suppose distr is a Distribution with distr.dimkeys = ['foo', 'bar', 'baz'] - If pdf is decorated with accepts_dict, we can call it like so - 1) distr.pdf(1, 2, 3) - 2) distr.pdf(foo=1, bar=2, baz=3) - 3) distr.pdf({'foo': 1, 'bar': 2, 'baz': 3}) - - Args: - method: The method accepting the different values for each dimensions - Returns: - method: The modified method to handle dictionaries of input data - """ - @wraps(method) - def f(self, *xs, **kwargs): - if xs: - # If xs is passed in, we check if the user passed it as each - # dimension separately or as a dictionary - if isinstance(xs[0], dict): - # If the first element of xs, is a dict, assume only element. - value_dict = xs[0] - values = [value_dict[key] for key in self.dimkeys] - else: - # Otherwise, it is a list of the values at each dimension - values = xs - else: - # Otherwise, we expect the values to be passed as keyword args. - values = [kwargs[key] for key in self.dimkeys] - return method(self, *values) - - return f - - -def returns_dict(method): - """ - This function decorator will allow descendants of MultivariateDistribution - which have methods which return values for each dimension to instead - return a dictionary of values mapping dimension name to value. - - This adds an as_dict argument which if set to True, will pack the return - value into a dictionary assuming the order is in that of the dimkeys - attribute of the distribution. - - The as_dict argument must be passed by keyword. - - Args: - method: The method which returns a list of values for each dimension - Returns: - method: The modified method to return a dictionary if specified to - """ - - @wraps(method) - def f(self, *pargs, as_dict=False, **kwargs): - values = method(self, *pargs, **kwargs) - if as_dict: - output = {key: value for key, value in zip(self.dimkeys, values)} - return output - else: - return values - - return f - - -def params_as_args(arg_names): - def decorator(method): - """ - - """ - @wraps(method) - def f(cls, x, params=None): - if params is None: - params = {} - for name in arg_names: - value = getattr(cls, name).value - if value is None: - message = """The {} parameter is unset. To use this method - it must be either called from an instance of - the distribution class or it must be called - directly from the class with a dictionary of - the parameters passed - with the params keyword.""".format(name) - raise ValueError(message) - params[name] = value - return method(cls, x, params) - return f - return decorator - -def params_as_args2(arg_names): - def decorator(method): - """ - - """ - @wraps(method) - def f(cls, x, y, params=None): - if params is None: - params = {} - for name in arg_names: - value = getattr(cls, name).value - if value is None: - message = """The {} parameter is unset. To use this method - it must be either called from an instance of - the distribution class or it must be called - directly from the class with a dictionary of - the parameters passed - with the params keyword.""".format(name) - raise ValueError(message) - params[name] = value - return method(cls, x, y, params) - return f - return decorator From c88e5b88732295bd0a1db77e8919ab67994dd071 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 16:01:31 -0700 Subject: [PATCH 14/17] boot-sp PR-2: install ipopt in the confidence intervals CI job The epi-spline distributions fit with a nonlinear solver, so with no ipopt on the runner the two epi-spline tests skip and splines.py never runs at all -- Smoothed_boot_epi ships with no verification behind it. This gets ipopt the way Pyomo's own CI does: apt-get the BLAS/LAPACK/gfortran libraries the binary links against, then unpack the IDAES idaes-ext solvers release (30MB) onto the PATH. `ipopt -v` at the end of the step means a bad download fails the job instead of quietly leaving the tests skipped. Verified locally by installing the same tarball: both epi-spline tests pass (serially and under mpiexec -np 2), splines.py goes from 6% to 39% covered and distributions.py from 87% to 97%, since the epi-spline class body now runs. Co-Authored-By: Claude Opus 5 --- .github/workflows/test_pr_and_main.yml | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index 04717a28a..8c544d0f8 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 . From e9ea6b661b84ea0214c92255b216a685c5317e47 Mon Sep 17 00:00:00 2001 From: Dave Woodruff Date: Sun, 26 Jul 2026 16:28:48 -0700 Subject: [PATCH 15/17] boot-sp PR-2: give the fitted student's t the variance it was asked for UnivariateStudentDistribution named its third argument the variance, and the fit docstring promised "the mean and variance of the distribution as the mean and variance of the data", but the constructor handed sqrt(var) to scipy as the *scale*. A t of scale s has variance s**2 * df/(df-2), so the distribution's variance came out as var * df/(df-2) instead of var -- and under the fit's own df rule of 2v/(v-1) that factor is exactly v, so fitting to data with variance 3.88 produced a t with variance 15.08. The parameter was a scale wearing the name of a variance. The constructor now solves that relation for the scale, so the variance argument is the variance for any df, and df <= 2 is refused: such a t has no finite variance at all, so it cannot be given one. That also settles what the df rule was for. 2v/(v-1) is the df at which a t of *scale one* has variance v, so with the scale derived the fitted scale comes out as one and the rule is self-consistent. What it cannot do is fit a sample variance of one or less -- a t of scale one always has variance above one -- which is what the old "Impossible to define a student distribution" print was about before it went on to build a df=1 Cauchy, a distribution with no variance whatsoever. That case is now a ValueError telling the caller to choose a df themselves or fit something else. Nothing in the repository fits this distribution: the smoothed methods only ever ask for univariate-epispline or univariate-kernel. Co-Authored-By: Claude Opus 5 --- .../bootsp/statdist/distributions.py | 41 +++++++++++++------ mpisppy/tests/test_boot_sp_smoothed.py | 36 ++++++++++++---- 2 files changed, 57 insertions(+), 20 deletions(-) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py index 5fbdb313b..8a1ff6adc 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -183,24 +183,34 @@ def __init__(self, df, mean, var): degrees of freedom, the mean, and the variance of the distribution. Args: - df (int): The number of degrees of freedom + df (int): The number of degrees of freedom, which must be more + than 2: a t distribution with df <= 2 has no finite variance, + so it cannot be given one here mean (float): The mean parameter var (float): The variance parameter """ + if df <= 2: + raise ValueError("A student's t distribution with df={} has no " + "finite variance, so it cannot be given the " + "variance {}.".format(df, var)) self.df = df self.mean = mean self.var = var - # We make the lower bound a very small number to exclude the - # possibility of 0 for the degrees of freedom. + # The degrees of freedom have to exceed 2 for the variance to exist. params = [Parameter('mean', mean), Parameter('variance', var, (0, None)), - Parameter('df', df, bounds=(epsilon, None))] + Parameter('df', df, bounds=(2, None))] UnivariateDistribution.__init__(self, params) + # scipy parameterizes the t by a scale, not by its variance: a t with + # scale s has variance s**2 * df/(df-2), so solve that for s. Passing + # sqrt(var) as the scale instead would give a distribution whose + # variance is var * df/(df-2), not var. + scale = np.sqrt(self.var * (self.df - 2) / self.df) self.distribution = scipy.stats.t(df=self.df, loc=self.mean, - scale=np.sqrt(self.var)) + scale=scale) @classmethod def fit(cls, data): @@ -209,19 +219,26 @@ def fit(cls, data): This will estimate the mean and variance of the distribution as the mean and variance of the data. + The degrees of freedom are taken to be 2v/(v-1) for a sample variance + v, which is the value at which a t of scale one has variance v; the + constructor then derives that scale of one. A sample variance of one or + less has no such value (a t of scale one always has variance above + one), and it is the caller who has to decide what to do about that, so + it is an error rather than a silent substitution. + Args: data (List[float]): The list of values to fit the data to Returns: UnivariateStudentDistribution: The fitted student's t distribution """ - var = np.var(data) + var = np.var(data) if var <= 1: - print('input_data gives Var < 1: ' - 'Impossible to define a student distribution') - print('Degree of freedom is by default set to 1') - df=1 - else: - df = 2*var/(var-1) + raise ValueError( + "The data have variance {}, and a student's t distribution " + "cannot be fit to a variance of one or less by this rule. " + "Either construct one directly with a chosen df > 2, or fit a " + "different distribution.".format(var)) + df = 2*var/(var-1) mean = np.mean(data) return UnivariateStudentDistribution(df, mean, var) diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index 6e0caa497..350c3fe21 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -226,29 +226,49 @@ def test_normal_generates_X(self): self.assertEqual(len(draws), 1000) self.assertLess(abs(float(np.mean(draws))), 0.25) - def test_student_fit_uses_the_documented_df(self): + def test_student_fit_matches_the_data_moments(self): + # the fit promises the distribution's mean and variance are the data's; + # passing sqrt(var) to scipy as the *scale* would instead give a + # variance of var*df/(df-2), i.e. var**2 under the df rule below data = list(np.random.RandomState(4).normal(0.0, 2.0, size=500)) var = float(np.var(data)) self.assertGreater(var, 1.0) st = distribution_factory("univariate-student").fit(data) - self.assertAlmostEqual(st.mean, float(np.mean(data))) - self.assertAlmostEqual(st.var, var) + self.assertAlmostEqual(st.distribution.mean(), float(np.mean(data))) + self.assertAlmostEqual(st.distribution.var(), var) self.assertAlmostEqual(st.df, 2*var/(var - 1)) + # that df is exactly the one at which a t of scale one has variance var + self.assertAlmostEqual(st.df / (st.df - 2), var) + + 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_falls_back_to_one_df(self): - # 2v/(v-1) is not a usable number of degrees of freedom once v <= 1 + def test_student_fit_refuses_low_variance_data(self): + # 2v/(v-1) is not a usable number of degrees of freedom once v <= 1, + # and what to do instead is the caller's decision data = list(np.random.RandomState(5).normal(0.0, 0.1, size=200)) self.assertLessEqual(float(np.var(data)), 1.0) - st = distribution_factory("univariate-student").fit(data) - self.assertEqual(st.df, 1) + with self.assertRaises(ValueError): + distribution_factory("univariate-student").fit(data) 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 location and scale as a normal, but with the fatter tails + # 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)) From ebd77be5c5b7483eddc89bf0c1fcfba6393e23c8 Mon Sep 17 00:00:00 2001 From: David L Woodruff Date: Mon, 27 Jul 2026 15:47:32 -0700 Subject: [PATCH 16/17] statdist: fit the student's t df by kurtosis, not 2v/(v-1) (#819 part 2) Now that the constructor derives scipy's scale from (var, df), any df > 2 reproduces the requested variance, so the old df = 2v/(v-1) is just a free tail-heaviness knob -- and a poor one: it is not scale-free (df changes when the data are rescaled), its coupling runs backwards (larger variance -> df toward 2 -> heavier tails), and it cannot fit a sample variance <= 1. Replace it with method of moments on the excess kurtosis: a t with df > 4 has excess kurtosis 6/(df-4), so df = 4 + 6/excess_kurtosis. This is scale-free and data-driven. A sample with excess kurtosis at or below zero (tails no heavier than the normal) has no finite-variance t that matches it, so df falls back to a large "effectively normal" value instead of raising. The v <= 1 ValueError is gone -- any variance v > 0 fits. Update Test_statdist: drop the 2v/(v-1) df assertion and the low-variance refusal test; add tests that df is recovered from kurtosis (data from a t with df=5 recovers df near 5), that light-tailed data falls back to the normal, and that low-variance data now fits with its moments honored. Nothing in the repository fits this distribution (compute_smoothed_ci only requests epispline/kernel), so the blast radius is confined to these tests. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../bootsp/statdist/distributions.py | 46 ++++++++++------- mpisppy/tests/test_boot_sp_smoothed.py | 50 ++++++++++++++----- 2 files changed, 67 insertions(+), 29 deletions(-) diff --git a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py index 8a1ff6adc..cdc0ff7c0 100644 --- a/mpisppy/confidence_intervals/bootsp/statdist/distributions.py +++ b/mpisppy/confidence_intervals/bootsp/statdist/distributions.py @@ -212,34 +212,46 @@ def __init__(self, df, mean, var): self.distribution = scipy.stats.t(df=self.df, loc=self.mean, scale=scale) + # df selection for fit(): a student's t with df > 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): """ - This will fit a student's distribution to the passed in data. - This will estimate the mean and variance of the distribution as - the mean and variance of the 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). - The degrees of freedom are taken to be 2v/(v-1) for a sample variance - v, which is the value at which a t of scale one has variance v; the - constructor then derives that scale of one. A sample variance of one or - less has no such value (a t of scale one always has variance above - one), and it is the caller who has to decide what to do about that, so - it is an error rather than a silent substitution. + 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 """ - var = np.var(data) - if var <= 1: - raise ValueError( - "The data have variance {}, and a student's t distribution " - "cannot be fit to a variance of one or less by this rule. " - "Either construct one directly with a chosen df > 2, or fit a " - "different distribution.".format(var)) - df = 2*var/(var-1) 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) diff --git a/mpisppy/tests/test_boot_sp_smoothed.py b/mpisppy/tests/test_boot_sp_smoothed.py index 350c3fe21..90862af19 100644 --- a/mpisppy/tests/test_boot_sp_smoothed.py +++ b/mpisppy/tests/test_boot_sp_smoothed.py @@ -29,6 +29,7 @@ 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 @@ -228,17 +229,38 @@ def test_normal_generates_X(self): def test_student_fit_matches_the_data_moments(self): # the fit promises the distribution's mean and variance are the data's; - # passing sqrt(var) to scipy as the *scale* would instead give a - # variance of var*df/(df-2), i.e. var**2 under the df rule below + # 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)) - self.assertGreater(var, 1.0) st = distribution_factory("univariate-student").fit(data) self.assertAlmostEqual(st.distribution.mean(), float(np.mean(data))) self.assertAlmostEqual(st.distribution.var(), var) - self.assertAlmostEqual(st.df, 2*var/(var - 1)) - # that df is exactly the one at which a t of scale one has variance var - self.assertAlmostEqual(st.df / (st.df - 2), 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): @@ -253,13 +275,17 @@ def test_student_needs_a_df_that_has_a_variance(self): distribution_factory("univariate-student")( df=df, mean=0.0, var=1.0) - def test_student_fit_refuses_low_variance_data(self): - # 2v/(v-1) is not a usable number of degrees of freedom once v <= 1, - # and what to do instead is the caller's decision + 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)) - self.assertLessEqual(float(np.var(data)), 1.0) - with self.assertRaises(ValueError): - distribution_factory("univariate-student").fit(data) + 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) From 0f271ec848aeaf15aa88ca7589c5e9ad4732083a Mon Sep 17 00:00:00 2001 From: David L Woodruff Date: Tue, 28 Jul 2026 09:34:06 -0700 Subject: [PATCH 17/17] Design doc: refresh the boot-sp merge status (PR-1 and its test fix merged) PR-1 (#783) merged 2026-07-24 and the test_boot_sp.py np=2 fix (#820) merged 2026-07-28, so the Status block no longer describes PR-1 as an open draft. PR-2 is this branch (#818). Co-Authored-By: Claude Opus 5 (1M context) --- doc/designs/bootsp_merge_design.md | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/designs/bootsp_merge_design.md b/doc/designs/bootsp_merge_design.md index 6432052e4..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