easyvvuq.sampling.pce
1import logging 2import chaospy as cp 3import numpy as np 4import random 5from .base import BaseSamplingElement, Vary 6from .transformations import Transformations 7 8__author__ = "Jalal Lakhlili" 9__copyright__ = """ 10 11 Copyright 2018 Robin A. Richardson, David W. Wright 12 13 This file is part of EasyVVUQ 14 15 EasyVVUQ is free software: you can redistribute it and/or modify 16 it under the terms of the Lesser GNU General Public License as published by 17 the Free Software Foundation, either version 3 of the License, or 18 (at your option) any later version. 19 20 EasyVVUQ is distributed in the hope that it will be useful, 21 but WITHOUT ANY WARRANTY; without even the implied warranty of 22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 23 Lesser GNU General Public License for more details. 24 25 You should have received a copy of the Lesser GNU General Public License 26 along with this program. If not, see <https://www.gnu.org/licenses/>. 27 28""" 29__license__ = "LGPL" 30 31 32class PCESampler(BaseSamplingElement, sampler_name="PCE_sampler"): 33 def __init__(self, 34 vary=None, 35 distribution=None, 36 count=0, 37 polynomial_order=4, 38 regression=False, 39 rule="G", 40 sparse=False, 41 growth=False, 42 relative_analysis=False, 43 nominal_value=None): 44 """ 45 Create the sampler for the Polynomial Chaos Expansion using 46 pseudo-spectral projection or regression (Point Collocation). 47 48 Parameters 49 ---------- 50 vary: dict or None 51 keys = parameters to be sampled, values = distributions. 52 53 distribution: cp.Distribution or matrix-like 54 Joint distribution specifying dependency between the parameters or 55 correlation matrix of the parameters. Depending on the type of the argument 56 either Rosenblatt or Cholesky transformation will be used to handle the 57 dependent parameters. 58 59 count : int, optional 60 Specified counter for Fast forward, default is 0. 61 62 polynomial_order : int, optional 63 The polynomial order, default is 4. 64 65 regression : bool, optional 66 If True, regression variante (point collecation) will be used, 67 otherwise projection variante (pseud-spectral) will be used. 68 Default value is False. 69 70 rule : char, optional 71 The quadrature method, in case of projection (default is Gaussian "G"). 72 The sequence sampler in case of regression (default is Hammersley "M") 73 74 sparse : bool, optional 75 If True, use Smolyak sparse grid instead of normal tensor product 76 grid. Default value is False. 77 78 growth (bool, None), optional 79 If True, quadrature point became nested. 80 81 relative_analysis (bool, None), optional 82 If True, we add one additional sample with all parameters having zero (nominal) value. 83 This is used in the relative analysis, where the model output is represented 84 relative to the nominal output, and similarly, the parameters represent the delta of 85 the parameter nominal value (i.e. zero represents parameter's nominal value, nominal + delta*nominal) 86 87 nominal_value : dict, optional 88 Evaluate derivative of the model at the nominal value of the parameters. 89 It should be a dict with the keys which are present in vary. 90 In case the base_value is None, the mean of the distribution is used (assuming cp.Normal). 91 """ 92 93 if vary is None: 94 msg = ("'vary' cannot be None. RandomSampler must be passed a " 95 "dict of the names of the parameters you want to vary, " 96 "and their corresponding distributions.") 97 logging.error(msg) 98 raise Exception(msg) 99 if not isinstance(vary, dict): 100 msg = ("'vary' must be a dictionary of the names of the " 101 "parameters you want to vary, and their corresponding " 102 "distributions.") 103 logging.error(msg) 104 raise Exception(msg) 105 if len(vary) == 0: 106 msg = "'vary' cannot be empty." 107 logging.error(msg) 108 raise Exception(msg) 109 110 # Remember whether to add the extra run using the base value of the parameters (0 corresponding to the mean) 111 logging.info(f"Performing relative analysis: {relative_analysis}") 112 self.relative_analysis = relative_analysis 113 114 self.vary = Vary(vary) 115 self.polynomial_order = polynomial_order 116 117 # List of the probability distributions of uncertain parameters 118 self.params_distribution = list(vary.values()) 119 params_num = len(self.params_distribution) 120 121 # Nominal value of the parameters 122 if nominal_value is None: 123 # Assumes that v is cp.Normal() 124 if (all([type(v) == type(cp.Normal()) for v in vary.values()])): 125 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 126 logging.info(f"Using parameter mean for the relative analysis {nominal_value}") 127 else: 128 if (len(nominal_value) != params_num): 129 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 130 logging.error(msg) 131 raise ValueError(msg) 132 logging.info(f"Using user-provided nominal value for relative analysis {nominal_value}") 133 self.nominal_value = nominal_value 134 135 # Multivariate distribution, the behaviour changes based on the 136 # 'distribution' argument, which can be: 137 # None - use default joint 138 # cp.Distribution - use Rosenblatt if the distribution is dependent 139 # matrix-lie - use Cholesky 140 self._is_dependent = False 141 self._transformation = None 142 self.distribution_dep = None 143 if distribution is None: 144 logging.info("Using default joint distribution") 145 self.distribution = cp.J(*self.params_distribution) 146 elif 'distributions' in str(type(distribution)): 147 if distribution.stochastic_dependent: 148 if not isinstance(distribution, cp.MvNormal): 149 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 150 if not len(distribution._parameters['mean']) == params_num: 151 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 152 logging.info("Using user provided joint distribution with Rosenblatt transformation") 153 154 self._is_dependent = True 155 self._transformation = "Rosenblatt" 156 self.distribution_dep = distribution 157 else: 158 logging.info("Using user provided joint distribution without any transformation") 159 self.distribution = distribution 160 assert(self._is_dependent == False) 161 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 162 if not len(distribution) == params_num: 163 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 164 for i in range(params_num): 165 if not distribution[i][i] == 1.0: 166 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 167 logging.info("Using user provided correlation matrix for Cholesky transformation") 168 169 self._is_dependent = True 170 self._transformation = "Cholesky" 171 self.distribution_dep = np.array(distribution) 172 else: 173 logging.error("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 174 exit() 175 176 # Build independent joint multivariate distribution considering each uncertain paramter 177 if not self.distribution_dep is None: 178 179 params_distribution = [vary_dist for vary_dist in vary.values()] 180 self.distribution = cp.J(*params_distribution) 181 182 # This assumes that the order of the parameters in distribution and distribution_dep is the same 183 # and the distribution type is cp.Normal 184 for id_v, v in enumerate(vary): 185 assert(type(vary[v]) == type(cp.Normal())) 186 if self._transformation == "Rosenblatt": 187 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution_dep._parameters['mean'][id_v]) 188 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution[id_v].get_mom_parameters()['shift'][0]) 189 logging.debug(f"The independent distribution consists of: {self.distribution}") 190 logging.debug(f"Using parameter permutation: {list(vary.keys())}") 191 192 # The orthogonal polynomials corresponding to the joint distribution 193 self.P = cp.expansion.stieltjes(polynomial_order, self.distribution, normed=True) 194 195 # The quadrature information 196 self.quad_sparse = sparse 197 self.rule = rule 198 199 # Clenshaw-Curtis should be nested if sparse (#139 chaospy issue) 200 self.quad_growth = growth 201 cc = ['c', 'C', 'clenshaw_curtis', 'Clenshaw_Curtis'] 202 if sparse and rule in cc: 203 self.quad_growth = True 204 205 # To determinate the PCE vraint to use 206 self.regression = regression 207 208 # indices of cp.DiscreteUniform parameters 209 idx_discrete = np.where([isinstance(p, cp.DiscreteUniform) for p in self.params_distribution])[0] 210 211 # Regression variante (Point collocation method) 212 if regression: 213 logging.info(f"Using point collocation method to create PCE") 214 # Change the default rule 215 if rule == "G": 216 self.rule = "M" 217 218 # Generates samples 219 self._n_samples = 2 * len(self.P) 220 logging.info(f"Generating {self._n_samples} samples using {self.rule} rule") 221 self._nodes = cp.generate_samples(order=self._n_samples, 222 domain=self.distribution, 223 rule=self.rule) 224 225 # Transform relative nodes to absolute nodes 226 if self.relative_analysis: 227 for pi,p in enumerate(vary.keys()): 228 self._nodes[pi] = (1.0 + self._nodes[pi]) * nominal_value[p] 229 230 self._weights = None 231 232 # Nodes transformation 233 if self._is_dependent: 234 if self._transformation == "Rosenblatt": 235 logging.info("Performing Rosenblatt transformation") 236 self._nodes_dep = Transformations.rosenblatt(self._nodes, self.distribution, self.distribution_dep, regression) 237 elif self._transformation == "Cholesky": 238 logging.info("Performing Cholesky transformation") 239 self._nodes_dep = Transformations.cholesky(self._nodes, self.vary, self.distribution_dep, regression) 240 else: 241 logging.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 242 exit() 243 244 # Projection variante (Pseudo-spectral method) 245 else: 246 if type(self.rule) is str and idx_discrete.size > 0: 247 tmp = [] 248 for i in range(params_num): 249 if i in idx_discrete: 250 tmp.append("discrete") 251 else: 252 tmp.append(rule) 253 self.rule = tmp 254 255 logging.info(f"Using pseudo-spectral method to create PCE") 256 # Nodes and weights for the integration 257 self._nodes, self._weights = cp.generate_quadrature(order=polynomial_order, 258 dist=self.distribution, 259 rule=self.rule, 260 sparse=sparse, 261 growth=self.quad_growth) 262 # Number of samples 263 self._n_samples = len(self._nodes[0]) 264 logging.info(f"Generated {self._n_samples} nodes/weights pairs using {self.rule} rule") 265 266 # Nodes transformation 267 if self._is_dependent: 268 # Scale the independent nodes 269 raise NotImplementedError(f'Transformation of the independent nodes not supported with {regression = }') 270 271 if self.relative_analysis: 272 raise NotImplementedError(f'Transformation of the relative nodes not supported with {regression = }') 273 274 # Fast forward to specified count, if possible 275 self.count = 0 276 if self.count >= self._n_samples: 277 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 278 f"but sampler only has {self.n_samples} samples, therefore" 279 f"this sampler will not provide any more samples.") 280 logging.warning(msg) 281 else: 282 for i in range(count): 283 self.__next__() 284 285 def is_finite(self): 286 return True 287 288 @property 289 def n_samples(self): 290 """ 291 Number of samples (Ns) of PCE method. 292 - When using pseudo-spectral projection method with tensored 293 quadrature: Ns = (p + 1)**d 294 - When using pseudo-spectral projection method with sparce grid 295 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 296 - When using regression method: Ns = 2*(p + d)!/p!*d! 297 Where: p is the polynomial degree and d is the number of 298 uncertain parameters. 299 300 Ref: Eck et al. 'A guide to uncertainty quantification and 301 sensitivity analysis for cardiovascular applications' [2016]. 302 """ 303 return self._n_samples 304 305 @property 306 def analysis_class(self): 307 """Return a corresponding analysis class. 308 """ 309 from easyvvuq.analysis import PCEAnalysis 310 return PCEAnalysis 311 312 def __next__(self): 313 if self.count < self._n_samples: #base Train samples used to evaluate the PCE 314 run_dict = {} 315 for i, param_name in enumerate(self.vary.vary_dict): 316 # These are nodes that need to be returned as samples o be used for the model execution, 317 # for the SA in EasyVVUQ we will use only the raw independent nodes 318 if self._is_dependent: 319 # Return transformed nodes reflecting the dependencies 320 run_dict[param_name] = self._nodes_dep[i][self.count] 321 else: 322 current_param = self._nodes[i][self.count] 323 # all parameters self.xi_d will store floats. If current param is 324 # DiscreteUniform, convert e.g. 2.0 to 2 before running 325 # the simulation. 326 if isinstance(self.params_distribution[i], cp.DiscreteUniform): 327 current_param = int(current_param) 328 run_dict[param_name] = current_param 329 self.count += 1 330 return run_dict 331 elif self.relative_analysis and self.count == self._n_samples: #extra sample for the nominal case 332 run_dict = self.nominal_value 333 self.count += 1 334 return run_dict 335 else: 336 raise StopIteration
33class PCESampler(BaseSamplingElement, sampler_name="PCE_sampler"): 34 def __init__(self, 35 vary=None, 36 distribution=None, 37 count=0, 38 polynomial_order=4, 39 regression=False, 40 rule="G", 41 sparse=False, 42 growth=False, 43 relative_analysis=False, 44 nominal_value=None): 45 """ 46 Create the sampler for the Polynomial Chaos Expansion using 47 pseudo-spectral projection or regression (Point Collocation). 48 49 Parameters 50 ---------- 51 vary: dict or None 52 keys = parameters to be sampled, values = distributions. 53 54 distribution: cp.Distribution or matrix-like 55 Joint distribution specifying dependency between the parameters or 56 correlation matrix of the parameters. Depending on the type of the argument 57 either Rosenblatt or Cholesky transformation will be used to handle the 58 dependent parameters. 59 60 count : int, optional 61 Specified counter for Fast forward, default is 0. 62 63 polynomial_order : int, optional 64 The polynomial order, default is 4. 65 66 regression : bool, optional 67 If True, regression variante (point collecation) will be used, 68 otherwise projection variante (pseud-spectral) will be used. 69 Default value is False. 70 71 rule : char, optional 72 The quadrature method, in case of projection (default is Gaussian "G"). 73 The sequence sampler in case of regression (default is Hammersley "M") 74 75 sparse : bool, optional 76 If True, use Smolyak sparse grid instead of normal tensor product 77 grid. Default value is False. 78 79 growth (bool, None), optional 80 If True, quadrature point became nested. 81 82 relative_analysis (bool, None), optional 83 If True, we add one additional sample with all parameters having zero (nominal) value. 84 This is used in the relative analysis, where the model output is represented 85 relative to the nominal output, and similarly, the parameters represent the delta of 86 the parameter nominal value (i.e. zero represents parameter's nominal value, nominal + delta*nominal) 87 88 nominal_value : dict, optional 89 Evaluate derivative of the model at the nominal value of the parameters. 90 It should be a dict with the keys which are present in vary. 91 In case the base_value is None, the mean of the distribution is used (assuming cp.Normal). 92 """ 93 94 if vary is None: 95 msg = ("'vary' cannot be None. RandomSampler must be passed a " 96 "dict of the names of the parameters you want to vary, " 97 "and their corresponding distributions.") 98 logging.error(msg) 99 raise Exception(msg) 100 if not isinstance(vary, dict): 101 msg = ("'vary' must be a dictionary of the names of the " 102 "parameters you want to vary, and their corresponding " 103 "distributions.") 104 logging.error(msg) 105 raise Exception(msg) 106 if len(vary) == 0: 107 msg = "'vary' cannot be empty." 108 logging.error(msg) 109 raise Exception(msg) 110 111 # Remember whether to add the extra run using the base value of the parameters (0 corresponding to the mean) 112 logging.info(f"Performing relative analysis: {relative_analysis}") 113 self.relative_analysis = relative_analysis 114 115 self.vary = Vary(vary) 116 self.polynomial_order = polynomial_order 117 118 # List of the probability distributions of uncertain parameters 119 self.params_distribution = list(vary.values()) 120 params_num = len(self.params_distribution) 121 122 # Nominal value of the parameters 123 if nominal_value is None: 124 # Assumes that v is cp.Normal() 125 if (all([type(v) == type(cp.Normal()) for v in vary.values()])): 126 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 127 logging.info(f"Using parameter mean for the relative analysis {nominal_value}") 128 else: 129 if (len(nominal_value) != params_num): 130 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 131 logging.error(msg) 132 raise ValueError(msg) 133 logging.info(f"Using user-provided nominal value for relative analysis {nominal_value}") 134 self.nominal_value = nominal_value 135 136 # Multivariate distribution, the behaviour changes based on the 137 # 'distribution' argument, which can be: 138 # None - use default joint 139 # cp.Distribution - use Rosenblatt if the distribution is dependent 140 # matrix-lie - use Cholesky 141 self._is_dependent = False 142 self._transformation = None 143 self.distribution_dep = None 144 if distribution is None: 145 logging.info("Using default joint distribution") 146 self.distribution = cp.J(*self.params_distribution) 147 elif 'distributions' in str(type(distribution)): 148 if distribution.stochastic_dependent: 149 if not isinstance(distribution, cp.MvNormal): 150 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 151 if not len(distribution._parameters['mean']) == params_num: 152 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 153 logging.info("Using user provided joint distribution with Rosenblatt transformation") 154 155 self._is_dependent = True 156 self._transformation = "Rosenblatt" 157 self.distribution_dep = distribution 158 else: 159 logging.info("Using user provided joint distribution without any transformation") 160 self.distribution = distribution 161 assert(self._is_dependent == False) 162 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 163 if not len(distribution) == params_num: 164 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 165 for i in range(params_num): 166 if not distribution[i][i] == 1.0: 167 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 168 logging.info("Using user provided correlation matrix for Cholesky transformation") 169 170 self._is_dependent = True 171 self._transformation = "Cholesky" 172 self.distribution_dep = np.array(distribution) 173 else: 174 logging.error("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 175 exit() 176 177 # Build independent joint multivariate distribution considering each uncertain paramter 178 if not self.distribution_dep is None: 179 180 params_distribution = [vary_dist for vary_dist in vary.values()] 181 self.distribution = cp.J(*params_distribution) 182 183 # This assumes that the order of the parameters in distribution and distribution_dep is the same 184 # and the distribution type is cp.Normal 185 for id_v, v in enumerate(vary): 186 assert(type(vary[v]) == type(cp.Normal())) 187 if self._transformation == "Rosenblatt": 188 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution_dep._parameters['mean'][id_v]) 189 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution[id_v].get_mom_parameters()['shift'][0]) 190 logging.debug(f"The independent distribution consists of: {self.distribution}") 191 logging.debug(f"Using parameter permutation: {list(vary.keys())}") 192 193 # The orthogonal polynomials corresponding to the joint distribution 194 self.P = cp.expansion.stieltjes(polynomial_order, self.distribution, normed=True) 195 196 # The quadrature information 197 self.quad_sparse = sparse 198 self.rule = rule 199 200 # Clenshaw-Curtis should be nested if sparse (#139 chaospy issue) 201 self.quad_growth = growth 202 cc = ['c', 'C', 'clenshaw_curtis', 'Clenshaw_Curtis'] 203 if sparse and rule in cc: 204 self.quad_growth = True 205 206 # To determinate the PCE vraint to use 207 self.regression = regression 208 209 # indices of cp.DiscreteUniform parameters 210 idx_discrete = np.where([isinstance(p, cp.DiscreteUniform) for p in self.params_distribution])[0] 211 212 # Regression variante (Point collocation method) 213 if regression: 214 logging.info(f"Using point collocation method to create PCE") 215 # Change the default rule 216 if rule == "G": 217 self.rule = "M" 218 219 # Generates samples 220 self._n_samples = 2 * len(self.P) 221 logging.info(f"Generating {self._n_samples} samples using {self.rule} rule") 222 self._nodes = cp.generate_samples(order=self._n_samples, 223 domain=self.distribution, 224 rule=self.rule) 225 226 # Transform relative nodes to absolute nodes 227 if self.relative_analysis: 228 for pi,p in enumerate(vary.keys()): 229 self._nodes[pi] = (1.0 + self._nodes[pi]) * nominal_value[p] 230 231 self._weights = None 232 233 # Nodes transformation 234 if self._is_dependent: 235 if self._transformation == "Rosenblatt": 236 logging.info("Performing Rosenblatt transformation") 237 self._nodes_dep = Transformations.rosenblatt(self._nodes, self.distribution, self.distribution_dep, regression) 238 elif self._transformation == "Cholesky": 239 logging.info("Performing Cholesky transformation") 240 self._nodes_dep = Transformations.cholesky(self._nodes, self.vary, self.distribution_dep, regression) 241 else: 242 logging.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 243 exit() 244 245 # Projection variante (Pseudo-spectral method) 246 else: 247 if type(self.rule) is str and idx_discrete.size > 0: 248 tmp = [] 249 for i in range(params_num): 250 if i in idx_discrete: 251 tmp.append("discrete") 252 else: 253 tmp.append(rule) 254 self.rule = tmp 255 256 logging.info(f"Using pseudo-spectral method to create PCE") 257 # Nodes and weights for the integration 258 self._nodes, self._weights = cp.generate_quadrature(order=polynomial_order, 259 dist=self.distribution, 260 rule=self.rule, 261 sparse=sparse, 262 growth=self.quad_growth) 263 # Number of samples 264 self._n_samples = len(self._nodes[0]) 265 logging.info(f"Generated {self._n_samples} nodes/weights pairs using {self.rule} rule") 266 267 # Nodes transformation 268 if self._is_dependent: 269 # Scale the independent nodes 270 raise NotImplementedError(f'Transformation of the independent nodes not supported with {regression = }') 271 272 if self.relative_analysis: 273 raise NotImplementedError(f'Transformation of the relative nodes not supported with {regression = }') 274 275 # Fast forward to specified count, if possible 276 self.count = 0 277 if self.count >= self._n_samples: 278 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 279 f"but sampler only has {self.n_samples} samples, therefore" 280 f"this sampler will not provide any more samples.") 281 logging.warning(msg) 282 else: 283 for i in range(count): 284 self.__next__() 285 286 def is_finite(self): 287 return True 288 289 @property 290 def n_samples(self): 291 """ 292 Number of samples (Ns) of PCE method. 293 - When using pseudo-spectral projection method with tensored 294 quadrature: Ns = (p + 1)**d 295 - When using pseudo-spectral projection method with sparce grid 296 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 297 - When using regression method: Ns = 2*(p + d)!/p!*d! 298 Where: p is the polynomial degree and d is the number of 299 uncertain parameters. 300 301 Ref: Eck et al. 'A guide to uncertainty quantification and 302 sensitivity analysis for cardiovascular applications' [2016]. 303 """ 304 return self._n_samples 305 306 @property 307 def analysis_class(self): 308 """Return a corresponding analysis class. 309 """ 310 from easyvvuq.analysis import PCEAnalysis 311 return PCEAnalysis 312 313 def __next__(self): 314 if self.count < self._n_samples: #base Train samples used to evaluate the PCE 315 run_dict = {} 316 for i, param_name in enumerate(self.vary.vary_dict): 317 # These are nodes that need to be returned as samples o be used for the model execution, 318 # for the SA in EasyVVUQ we will use only the raw independent nodes 319 if self._is_dependent: 320 # Return transformed nodes reflecting the dependencies 321 run_dict[param_name] = self._nodes_dep[i][self.count] 322 else: 323 current_param = self._nodes[i][self.count] 324 # all parameters self.xi_d will store floats. If current param is 325 # DiscreteUniform, convert e.g. 2.0 to 2 before running 326 # the simulation. 327 if isinstance(self.params_distribution[i], cp.DiscreteUniform): 328 current_param = int(current_param) 329 run_dict[param_name] = current_param 330 self.count += 1 331 return run_dict 332 elif self.relative_analysis and self.count == self._n_samples: #extra sample for the nominal case 333 run_dict = self.nominal_value 334 self.count += 1 335 return run_dict 336 else: 337 raise StopIteration
Baseclass for all EasyVVUQ sampling elements.
Attributes
- sampler_name (str): Name of the particular sampler.
PCESampler( vary=None, distribution=None, count=0, polynomial_order=4, regression=False, rule='G', sparse=False, growth=False, relative_analysis=False, nominal_value=None)
34 def __init__(self, 35 vary=None, 36 distribution=None, 37 count=0, 38 polynomial_order=4, 39 regression=False, 40 rule="G", 41 sparse=False, 42 growth=False, 43 relative_analysis=False, 44 nominal_value=None): 45 """ 46 Create the sampler for the Polynomial Chaos Expansion using 47 pseudo-spectral projection or regression (Point Collocation). 48 49 Parameters 50 ---------- 51 vary: dict or None 52 keys = parameters to be sampled, values = distributions. 53 54 distribution: cp.Distribution or matrix-like 55 Joint distribution specifying dependency between the parameters or 56 correlation matrix of the parameters. Depending on the type of the argument 57 either Rosenblatt or Cholesky transformation will be used to handle the 58 dependent parameters. 59 60 count : int, optional 61 Specified counter for Fast forward, default is 0. 62 63 polynomial_order : int, optional 64 The polynomial order, default is 4. 65 66 regression : bool, optional 67 If True, regression variante (point collecation) will be used, 68 otherwise projection variante (pseud-spectral) will be used. 69 Default value is False. 70 71 rule : char, optional 72 The quadrature method, in case of projection (default is Gaussian "G"). 73 The sequence sampler in case of regression (default is Hammersley "M") 74 75 sparse : bool, optional 76 If True, use Smolyak sparse grid instead of normal tensor product 77 grid. Default value is False. 78 79 growth (bool, None), optional 80 If True, quadrature point became nested. 81 82 relative_analysis (bool, None), optional 83 If True, we add one additional sample with all parameters having zero (nominal) value. 84 This is used in the relative analysis, where the model output is represented 85 relative to the nominal output, and similarly, the parameters represent the delta of 86 the parameter nominal value (i.e. zero represents parameter's nominal value, nominal + delta*nominal) 87 88 nominal_value : dict, optional 89 Evaluate derivative of the model at the nominal value of the parameters. 90 It should be a dict with the keys which are present in vary. 91 In case the base_value is None, the mean of the distribution is used (assuming cp.Normal). 92 """ 93 94 if vary is None: 95 msg = ("'vary' cannot be None. RandomSampler must be passed a " 96 "dict of the names of the parameters you want to vary, " 97 "and their corresponding distributions.") 98 logging.error(msg) 99 raise Exception(msg) 100 if not isinstance(vary, dict): 101 msg = ("'vary' must be a dictionary of the names of the " 102 "parameters you want to vary, and their corresponding " 103 "distributions.") 104 logging.error(msg) 105 raise Exception(msg) 106 if len(vary) == 0: 107 msg = "'vary' cannot be empty." 108 logging.error(msg) 109 raise Exception(msg) 110 111 # Remember whether to add the extra run using the base value of the parameters (0 corresponding to the mean) 112 logging.info(f"Performing relative analysis: {relative_analysis}") 113 self.relative_analysis = relative_analysis 114 115 self.vary = Vary(vary) 116 self.polynomial_order = polynomial_order 117 118 # List of the probability distributions of uncertain parameters 119 self.params_distribution = list(vary.values()) 120 params_num = len(self.params_distribution) 121 122 # Nominal value of the parameters 123 if nominal_value is None: 124 # Assumes that v is cp.Normal() 125 if (all([type(v) == type(cp.Normal()) for v in vary.values()])): 126 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 127 logging.info(f"Using parameter mean for the relative analysis {nominal_value}") 128 else: 129 if (len(nominal_value) != params_num): 130 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 131 logging.error(msg) 132 raise ValueError(msg) 133 logging.info(f"Using user-provided nominal value for relative analysis {nominal_value}") 134 self.nominal_value = nominal_value 135 136 # Multivariate distribution, the behaviour changes based on the 137 # 'distribution' argument, which can be: 138 # None - use default joint 139 # cp.Distribution - use Rosenblatt if the distribution is dependent 140 # matrix-lie - use Cholesky 141 self._is_dependent = False 142 self._transformation = None 143 self.distribution_dep = None 144 if distribution is None: 145 logging.info("Using default joint distribution") 146 self.distribution = cp.J(*self.params_distribution) 147 elif 'distributions' in str(type(distribution)): 148 if distribution.stochastic_dependent: 149 if not isinstance(distribution, cp.MvNormal): 150 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 151 if not len(distribution._parameters['mean']) == params_num: 152 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 153 logging.info("Using user provided joint distribution with Rosenblatt transformation") 154 155 self._is_dependent = True 156 self._transformation = "Rosenblatt" 157 self.distribution_dep = distribution 158 else: 159 logging.info("Using user provided joint distribution without any transformation") 160 self.distribution = distribution 161 assert(self._is_dependent == False) 162 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 163 if not len(distribution) == params_num: 164 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 165 for i in range(params_num): 166 if not distribution[i][i] == 1.0: 167 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 168 logging.info("Using user provided correlation matrix for Cholesky transformation") 169 170 self._is_dependent = True 171 self._transformation = "Cholesky" 172 self.distribution_dep = np.array(distribution) 173 else: 174 logging.error("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 175 exit() 176 177 # Build independent joint multivariate distribution considering each uncertain paramter 178 if not self.distribution_dep is None: 179 180 params_distribution = [vary_dist for vary_dist in vary.values()] 181 self.distribution = cp.J(*params_distribution) 182 183 # This assumes that the order of the parameters in distribution and distribution_dep is the same 184 # and the distribution type is cp.Normal 185 for id_v, v in enumerate(vary): 186 assert(type(vary[v]) == type(cp.Normal())) 187 if self._transformation == "Rosenblatt": 188 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution_dep._parameters['mean'][id_v]) 189 assert(vary[v].get_mom_parameters()['shift'][0] == self.distribution[id_v].get_mom_parameters()['shift'][0]) 190 logging.debug(f"The independent distribution consists of: {self.distribution}") 191 logging.debug(f"Using parameter permutation: {list(vary.keys())}") 192 193 # The orthogonal polynomials corresponding to the joint distribution 194 self.P = cp.expansion.stieltjes(polynomial_order, self.distribution, normed=True) 195 196 # The quadrature information 197 self.quad_sparse = sparse 198 self.rule = rule 199 200 # Clenshaw-Curtis should be nested if sparse (#139 chaospy issue) 201 self.quad_growth = growth 202 cc = ['c', 'C', 'clenshaw_curtis', 'Clenshaw_Curtis'] 203 if sparse and rule in cc: 204 self.quad_growth = True 205 206 # To determinate the PCE vraint to use 207 self.regression = regression 208 209 # indices of cp.DiscreteUniform parameters 210 idx_discrete = np.where([isinstance(p, cp.DiscreteUniform) for p in self.params_distribution])[0] 211 212 # Regression variante (Point collocation method) 213 if regression: 214 logging.info(f"Using point collocation method to create PCE") 215 # Change the default rule 216 if rule == "G": 217 self.rule = "M" 218 219 # Generates samples 220 self._n_samples = 2 * len(self.P) 221 logging.info(f"Generating {self._n_samples} samples using {self.rule} rule") 222 self._nodes = cp.generate_samples(order=self._n_samples, 223 domain=self.distribution, 224 rule=self.rule) 225 226 # Transform relative nodes to absolute nodes 227 if self.relative_analysis: 228 for pi,p in enumerate(vary.keys()): 229 self._nodes[pi] = (1.0 + self._nodes[pi]) * nominal_value[p] 230 231 self._weights = None 232 233 # Nodes transformation 234 if self._is_dependent: 235 if self._transformation == "Rosenblatt": 236 logging.info("Performing Rosenblatt transformation") 237 self._nodes_dep = Transformations.rosenblatt(self._nodes, self.distribution, self.distribution_dep, regression) 238 elif self._transformation == "Cholesky": 239 logging.info("Performing Cholesky transformation") 240 self._nodes_dep = Transformations.cholesky(self._nodes, self.vary, self.distribution_dep, regression) 241 else: 242 logging.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 243 exit() 244 245 # Projection variante (Pseudo-spectral method) 246 else: 247 if type(self.rule) is str and idx_discrete.size > 0: 248 tmp = [] 249 for i in range(params_num): 250 if i in idx_discrete: 251 tmp.append("discrete") 252 else: 253 tmp.append(rule) 254 self.rule = tmp 255 256 logging.info(f"Using pseudo-spectral method to create PCE") 257 # Nodes and weights for the integration 258 self._nodes, self._weights = cp.generate_quadrature(order=polynomial_order, 259 dist=self.distribution, 260 rule=self.rule, 261 sparse=sparse, 262 growth=self.quad_growth) 263 # Number of samples 264 self._n_samples = len(self._nodes[0]) 265 logging.info(f"Generated {self._n_samples} nodes/weights pairs using {self.rule} rule") 266 267 # Nodes transformation 268 if self._is_dependent: 269 # Scale the independent nodes 270 raise NotImplementedError(f'Transformation of the independent nodes not supported with {regression = }') 271 272 if self.relative_analysis: 273 raise NotImplementedError(f'Transformation of the relative nodes not supported with {regression = }') 274 275 # Fast forward to specified count, if possible 276 self.count = 0 277 if self.count >= self._n_samples: 278 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 279 f"but sampler only has {self.n_samples} samples, therefore" 280 f"this sampler will not provide any more samples.") 281 logging.warning(msg) 282 else: 283 for i in range(count): 284 self.__next__()
Create the sampler for the Polynomial Chaos Expansion using pseudo-spectral projection or regression (Point Collocation).
Parameters
- vary (dict or None): keys = parameters to be sampled, values = distributions.
- distribution (cp.Distribution or matrix-like): Joint distribution specifying dependency between the parameters or correlation matrix of the parameters. Depending on the type of the argument either Rosenblatt or Cholesky transformation will be used to handle the dependent parameters.
- count (int, optional): Specified counter for Fast forward, default is 0.
- polynomial_order (int, optional): The polynomial order, default is 4.
- regression (bool, optional): If True, regression variante (point collecation) will be used, otherwise projection variante (pseud-spectral) will be used. Default value is False.
- rule (char, optional): The quadrature method, in case of projection (default is Gaussian "G"). The sequence sampler in case of regression (default is Hammersley "M")
- sparse (bool, optional): If True, use Smolyak sparse grid instead of normal tensor product grid. Default value is False.
- growth (bool, None), optional: If True, quadrature point became nested.
- relative_analysis (bool, None), optional: If True, we add one additional sample with all parameters having zero (nominal) value. This is used in the relative analysis, where the model output is represented relative to the nominal output, and similarly, the parameters represent the delta of the parameter nominal value (i.e. zero represents parameter's nominal value, nominal + delta*nominal)
- nominal_value (dict, optional): Evaluate derivative of the model at the nominal value of the parameters. It should be a dict with the keys which are present in vary. In case the base_value is None, the mean of the distribution is used (assuming cp.Normal).
n_samples
289 @property 290 def n_samples(self): 291 """ 292 Number of samples (Ns) of PCE method. 293 - When using pseudo-spectral projection method with tensored 294 quadrature: Ns = (p + 1)**d 295 - When using pseudo-spectral projection method with sparce grid 296 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 297 - When using regression method: Ns = 2*(p + d)!/p!*d! 298 Where: p is the polynomial degree and d is the number of 299 uncertain parameters. 300 301 Ref: Eck et al. 'A guide to uncertainty quantification and 302 sensitivity analysis for cardiovascular applications' [2016]. 303 """ 304 return self._n_samples
Number of samples (Ns) of PCE method.
- When using pseudo-spectral projection method with tensored quadrature: Ns = (p + 1)**d
- When using pseudo-spectral projection method with sparce grid quadratue: Ns = bigO((p + 1)log(p + 1)*(d-1))
- When using regression method: Ns = 2(p + d)!/p!d! Where: p is the polynomial degree and d is the number of uncertain parameters.
Ref: Eck et al. 'A guide to uncertainty quantification and sensitivity analysis for cardiovascular applications' [2016].