From c9b937d4f0c0c0740f70c427a0bf22f06ff5a52c Mon Sep 17 00:00:00 2001 From: jasherma Date: Fri, 22 Nov 2024 20:30:01 -0500 Subject: [PATCH] Add chi-squared confidence level attr to `EllipsoidalSet` --- .../pyros/tests/test_uncertainty_sets.py | 82 +++++++++++++++++ pyomo/contrib/pyros/uncertainty_sets.py | 90 ++++++++++++++++--- 2 files changed, 159 insertions(+), 13 deletions(-) diff --git a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py index 13195b3270e..e28e557fabd 100644 --- a/pyomo/contrib/pyros/tests/test_uncertainty_sets.py +++ b/pyomo/contrib/pyros/tests/test_uncertainty_sets.py @@ -20,6 +20,7 @@ attempt_import, numpy as np, numpy_available, + scipy as sp, scipy_available, ) from pyomo.environ import SolverFactory @@ -1863,6 +1864,11 @@ def test_normal_construction_and_update(self): np.testing.assert_allclose( scale, eset.scale, err_msg="EllipsoidalSet scale not as expected" ) + np.testing.assert_allclose( + sp.stats.chi2.cdf(scale ** 0.5, df=2), + eset.chi_sq_conf_lvl, + err_msg="EllipsoidalSet chi-squared confidence level not as expected", + ) # check attributes update new_center = [-1, -3] @@ -1886,6 +1892,40 @@ def test_normal_construction_and_update(self): np.testing.assert_allclose( new_scale, eset.scale, err_msg="EllipsoidalSet scale update not as expected" ) + np.testing.assert_allclose( + sp.stats.chi2.cdf(new_scale ** 0.5, df=2), + eset.chi_sq_conf_lvl, + err_msg="EllipsoidalSet chi-square confidence level update not as expected" + ) + + def test_normal_construction_and_update_chi_sq_conf_lvl(self): + """ + Test EllipsoidalSet constructor and setter + work normally when arguments are appropriate. + """ + init_conf_lvl = 0.95 + eset = EllipsoidalSet( + center=[0, 0, 0], + shape_matrix=np.eye(3), + scale=None, + chi_sq_conf_lvl=init_conf_lvl, + ) + + self.assertEqual(eset.chi_sq_conf_lvl, init_conf_lvl) + np.testing.assert_allclose( + sp.stats.chi2.isf(q=1 - init_conf_lvl, df=eset.dim), + eset.scale, + err_msg="EllipsoidalSet scale not as expected", + ) + + new_conf_lvl = 0.99 + eset.chi_sq_conf_lvl = new_conf_lvl + self.assertEqual(eset.chi_sq_conf_lvl, new_conf_lvl) + np.testing.assert_allclose( + sp.stats.chi2.isf(q=1 - new_conf_lvl, df=eset.dim), + eset.scale, + err_msg="EllipsoidalSet scale not as expected", + ) def test_error_on_ellipsoidal_dim_change(self): """ @@ -1926,6 +1966,48 @@ def test_error_on_neg_scale(self): with self.assertRaisesRegex(ValueError, exc_str): eset.scale = neg_scale + def test_error_invalid_chi_sq_conf_lvl(self): + """ + Test error when attempting to initialize with chi-squared + confidence level outside range. + """ + center = [0, 0] + shape_matrix = [[1, 0], [0, 2]] + invalid_conf_lvl = 1.001 + + exc_str = r"Ensure the confidence level is a value in \[0, 1\)." + + # error on construction + with self.assertRaisesRegex(ValueError, exc_str): + EllipsoidalSet( + center=center, + shape_matrix=shape_matrix, + scale=None, + chi_sq_conf_lvl=invalid_conf_lvl, + ) + + # error on updating valid ellipsoid + eset = EllipsoidalSet(center, shape_matrix, scale=None, chi_sq_conf_lvl=0.95) + with self.assertRaisesRegex(ValueError, exc_str): + eset.chi_sq_conf_lvl = invalid_conf_lvl + + # negative confidence level + eset = EllipsoidalSet(center, shape_matrix, scale=None, chi_sq_conf_lvl=0.95) + with self.assertRaisesRegex(ValueError, exc_str): + eset.chi_sq_conf_lvl = -0.1 + + def test_error_scale_chi_sq_conf_lvl_construction(self): + """ + Test exception raised if neither or both of + `scale` and `chi_sq_conf_lvl` are None. + """ + exc_str = r"Exactly one of `scale` and `chi_sq_conf_lvl` should be None" + with self.assertRaisesRegex(ValueError, exc_str): + EllipsoidalSet([0], [[1]], scale=None, chi_sq_conf_lvl=None) + + with self.assertRaisesRegex(ValueError, exc_str): + EllipsoidalSet([0], [[1]], scale=1, chi_sq_conf_lvl=0.95) + def test_error_on_shape_matrix_with_wrong_size(self): """ Test error in event EllipsoidalSet shape matrix diff --git a/pyomo/contrib/pyros/uncertainty_sets.py b/pyomo/contrib/pyros/uncertainty_sets.py index ed03aa7553a..3da37648528 100644 --- a/pyomo/contrib/pyros/uncertainty_sets.py +++ b/pyomo/contrib/pyros/uncertainty_sets.py @@ -2374,31 +2374,37 @@ class EllipsoidalSet(UncertaintySet): center : (N,) array-like Center of the ellipsoid. shape_matrix : (N, N) array-like - A positive definite matrix characterizing the shape - and orientation of the ellipsoid. - scale : numeric type, optional + A symmetric positive definite matrix characterizing + the shape and orientation of the ellipsoid. + scale : numeric type or None, optional Square of the factor by which to scale the semi-axes of the ellipsoid (i.e. the eigenvectors of the shape matrix). The default is `1`. + chi_sq_conf_lvl : numeric type or None, optional + (Fractional) chi-squared confidence level of the multivariate + normal distribution with mean `center` and covariance + matrix `shape_matrix`. + Exactly one of `scale` and `chi_sq_conf_lvl` should be + None; otherwise, an exception is raised. Examples -------- - 3D origin-centered unit hypersphere: + A 3D origin-centered unit ball: >>> from pyomo.contrib.pyros import EllipsoidalSet >>> import numpy as np - >>> hypersphere = EllipsoidalSet( + >>> ball = EllipsoidalSet( ... center=[0, 0, 0], ... shape_matrix=np.eye(3), ... scale=1, ... ) - >>> hypersphere.center + >>> ball.center array([0, 0, 0]) - >>> hypersphere.shape_matrix + >>> ball.shape_matrix array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) - >>> hypersphere.scale + >>> ball.scale 1 A 2D ellipsoid with custom rotation and scaling: @@ -2416,13 +2422,42 @@ class EllipsoidalSet(UncertaintySet): >>> rotated_ellipsoid.scale 0.5 + A 4D 95% confidence ellipsoid: + + >>> conf_ellipsoid = EllipsoidalSet( + ... center=np.zeros(4), + ... shape_matrix=np.diag(range(1, 5)), + ... scale=None, + ... chi_sq_conf_lvl=0.95, + ... ) + >>> conf_ellipsoid.center + array([0, 0, 0, 0]) + >>> conf_ellipsoid.shape_matrix + array([[1, 0, 0, 0]], + [0, 2, 0, 0]], + [0, 0, 3, 0]], + [0, 0, 0. 4]]) + >>> conf_ellipsoid.scale + ...9.4877... + >>> conf_ellipsoid.chi_sq_conf_lvl + 0.95 + """ - def __init__(self, center, shape_matrix, scale=1): + def __init__(self, center, shape_matrix, scale=1, chi_sq_conf_lvl=None): """Initialize self (see class docstring).""" self.center = center self.shape_matrix = shape_matrix - self.scale = scale + + if scale is not None and chi_sq_conf_lvl is None: + self.scale = scale + elif scale is None and chi_sq_conf_lvl is not None: + self.chi_sq_conf_lvl = chi_sq_conf_lvl + else: + raise ValueError( + "Exactly one of `scale` and `chi_sq_conf_lvl` should be " + f"None (got {scale=}, {chi_sq_conf_lvl=})" + ) @property def type(self): @@ -2456,7 +2491,7 @@ def center(self, val): if val_arr.size != self.dim: raise ValueError( "Attempting to set attribute 'center' of " - f"AxisAlignedEllipsoidalSet of dimension {self.dim} " + f"{type(self).__name__} of dimension {self.dim} " f"to value of dimension {val_arr.size}" ) @@ -2535,7 +2570,7 @@ def shape_matrix(self, val): if hasattr(self, "_center"): if not all(size == self.dim for size in shape_mat_arr.shape): raise ValueError( - f"EllipsoidalSet attribute 'shape_matrix' " + f"{type(self).__name__} attribute 'shape_matrix' " f"must be a square matrix of size " f"{self.dim} to match set dimension " f"(provided matrix with shape {shape_mat_arr.shape})" @@ -2558,12 +2593,41 @@ def scale(self, val): validate_arg_type("scale", val, valid_num_types, "a valid numeric type", False) if val < 0: raise ValueError( - "EllipsoidalSet attribute " + f"{type(self).__name__} attribute " f"'scale' must be a non-negative real " f"(provided value {val})" ) self._scale = val + self._chi_sq_conf_lvl = sp.stats.chi2.cdf(x=val ** 0.5, df=self.dim) + + @property + def chi_sq_conf_lvl(self): + """ + numeric type : (Fractional) chi-squared confidence level of the + multivariate normal distribution with mean ``self.origin`` + and covariance ``self.shape_matrix`` for ellipsoidal region + with square magnification factor ``self.scale``. + """ + return self._chi_sq_conf_lvl + + @chi_sq_conf_lvl.setter + def chi_sq_conf_lvl(self, val): + validate_arg_type( + "chi_sq_conf_lvl", + val, + valid_num_types, + "a valid numeric type", + False, + ) + self._chi_sq_conf_lvl = val + self._scale = sp.stats.chi2.isf(q=1 - val, df=self.dim) + if np.isnan(self._scale) or np.isinf(self._scale): + raise ValueError( + f"Squared scaling factor calculation for confidence level {val} " + f"and set dimension {self.dim} returned {self._scale}. " + "Ensure the confidence level is a value in [0, 1)." + ) @property def dim(self):