easyvvuq.sampling.fd
1#from hashlib import shake_128 2import logging 3import chaospy as cp 4import numpy as np 5import random 6from .base import BaseSamplingElement, Vary 7from .transformations import Transformations 8 9 10# DEBUG USI 11from os import stat, path 12from time import ctime 13import json 14 15__author__ = "Jalal Lakhlili" 16__copyright__ = """ 17 18 Copyright 2018 Robin A. Richardson, David W. Wright 19 20 This file is part of EasyVVUQ 21 22 EasyVVUQ is free software: you can redistribute it and/or modify 23 it under the terms of the Lesser GNU General Public License as published by 24 the Free Software Foundation, either version 3 of the License, or 25 (at your option) any later version. 26 27 EasyVVUQ is distributed in the hope that it will be useful, 28 but WITHOUT ANY WARRANTY; without even the implied warranty of 29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 30 Lesser GNU General Public License for more details. 31 32 You should have received a copy of the Lesser GNU General Public License 33 along with this program. If not, see <https://www.gnu.org/licenses/>. 34 35""" 36__license__ = "LGPL" 37 38 39class FDSampler(BaseSamplingElement, sampler_name="FD_sampler"): 40 def __init__(self, 41 vary=None, 42 distribution=None, 43 perturbation=0.05, 44 nominal_value=None, 45 count=0, 46 relative_analysis=False): 47 """ 48 Create the sampler for the Polynomial Chaos Expansion using 49 pseudo-spectral projection or regression (Point Collocation). 50 51 Parameters 52 ---------- 53 vary: dict or None 54 keys = parameters to be sampled, values = distributions. 55 56 distribution: cp.Distribution or matrix-like 57 Joint distribution specifying dependency between the parameters or 58 correlation matrix of the parameters. Depending on the type of the argument 59 either Rosenblatt or Cholesky transformation will be used to handle the 60 dependent parameters. 61 62 perturbation: float 63 Perturbation of the parameters used in the finite difference scheme 64 65 nominal_value : dict, optional 66 FD sampler perturbs the base value (NV) of the parameters by the value specified. 67 It should be a dict with the keys which are present in vary. 68 In case the nominal_value is None, the mean of the distribution is used (assuming cp.Normal). 69 70 count : int, optional 71 Specified counter for Fast forward, default is 0. 72 73 relative_analysis : bool, optional 74 Default False, the nodes are perturbed by the value specified in perturbation 75 such that n = NV + perturbation. If True, the nodes are perturbed by the relative 76 value specified in perturbation such that n = NB * (1 + perturbation). 77 78 """ 79 # Create and initialize the logger 80 self.logger = logging.getLogger(__name__) 81 self.logger.setLevel(logging.DEBUG) 82 83 # Logger is already configured, remove all handlers 84 if self.logger.hasHandlers(): 85 self.logger.handlers = [] 86 87 formatter = logging.Formatter('%(asctime)s:%(name)s:%(levelname)s:%(message)s') 88 89 file_handler = logging.FileHandler('FD.log') 90 file_handler.setLevel(logging.DEBUG) 91 file_handler.setFormatter(formatter) 92 93 stream_handler = logging.StreamHandler() 94 stream_handler.setFormatter(formatter) 95 96 self.logger.addHandler(file_handler) 97 self.logger.addHandler(stream_handler) 98 99 if vary is None: 100 msg = ("'vary' cannot be None. RandomSampler must be passed a " 101 "dict of the names of the parameters you want to vary, " 102 "and their corresponding distributions.") 103 self.logger.error(msg) 104 raise Exception(msg) 105 if not isinstance(vary, dict): 106 msg = ("'vary' must be a dictionary of the names of the " 107 "parameters you want to vary, and their corresponding " 108 "distributions.") 109 self.logger.error(msg) 110 raise Exception(msg) 111 if len(vary) == 0: 112 msg = "'vary' cannot be empty." 113 self.logger.error(msg) 114 raise Exception(msg) 115 116 self.vary = Vary(vary) 117 118 # List of the probability distributions of uncertain parameters 119 params_distribution = list(vary.values()) 120 params_num = len(params_distribution) 121 122 # Remember whether to add the extra run 123 self.logger.info(f"Performing relative analysis: {relative_analysis}") 124 self.relative_analysis = relative_analysis 125 self._perturbation = perturbation 126 127 # Perturbation of the parameters 128 if nominal_value is None: 129 self.logger.info(f"Performing perturbation of the nodes, base value = mean, with delta = {perturbation}") 130 # Assumes that v is cp.Normal() 131 assert(all([type(v) == type(cp.Normal()) for v in vary.values()])) 132 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 133 else: 134 if (len(nominal_value) != params_num): 135 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 136 self.logger.error(msg) 137 raise ValueError(msg) 138 self.logger.info(f"Performing perturbation of the nodes, base value = {nominal_value}, with delta = {perturbation}") 139 140 # Generate the perturbed values of the parameters for the FD 141 #FD = 0.5*(y_pos/y_base-1)/(delta) + 0.5*(y_neg/y_base - 1)/(-delta) 142 self.generate_nodes(nominal_value, vary, distribution) 143 144 # Fast forward to specified count, if possible 145 self.count = 0 146 if self.count >= self._n_samples: 147 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 148 f"but sampler only has {self.n_samples} samples, therefore" 149 f"this sampler will not provide any more samples.") 150 self.logger.warning(msg) 151 else: 152 for i in range(count): 153 self.__next__() 154 155 156 def permute_problem(self, perm, mean, sigma, cov): 157 # Apply selected permutation to the problem 158 nParams = len(mean) 159 # cyclic permutation i = perm_step, N = nParams 160 # perm = [i i+1 ... N-1 0 1 ... i-1] 161 assert(len(perm) == nParams) 162 # create permutation matrix 163 P = np.eye(nParams)[perm] 164 165 # Model properties that need to be permuted 166 mean_ = P@mean 167 sigma_ = P@sigma 168 cov_ = np.matmul(P, np.matmul(cov, P.transpose())) 169 170 return mean_, sigma_, cov_ 171 172 def generate_nodes(self, nominal_value, vary, distribution): 173 params_num = len(nominal_value) 174 self._n_samples = 2*params_num + 1 175 176 # Multivariate distribution, the behaviour changes based on the 177 # 'distribution' argument, which can be: 178 # None - use default joint 179 # cp.Distribution - use Rosenblatt if the distribution is dependent 180 # matrix-like - use Cholesky 181 self._is_dependent = False 182 self._transformation = None 183 self.distribution_dep = None 184 if 'distributions' in str(type(distribution)): 185 if not distribution.stochastic_dependent: 186 raise ValueError("User provided joint distribution needs to contain dependency between the parameters") 187 if not isinstance(distribution, cp.MvNormal): 188 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 189 if not len(distribution._parameters['mean']) == params_num: 190 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 191 self.logger.info("Using user provided joint distribution with Rosenblatt transformation") 192 193 self._is_dependent = True 194 self._transformation = "Rosenblatt" 195 self.distribution_dep = distribution 196 197 mu = distribution._parameters['mean'] 198 cov = distribution._covariance 199 200 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 201 if not len(distribution) == params_num: 202 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 203 for i in range(params_num): 204 if not distribution[i][i] == 1.0: 205 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 206 self.logger.info("Using user provided correlation matrix for Cholesky transformation") 207 208 self._is_dependent = True 209 self._transformation = "Cholesky" 210 self.distribution_dep = np.array(distribution) 211 elif distribution is None: 212 pass 213 else: 214 raise ValueError("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 215 216 # Set up independent perturbations 217 #G: [0, d, -d, 0, 0, 0, 0] 218 #C: [0, 0, 0, d, -d, 0, 0] 219 #E: [0, 0, 0, 0, 0, d, -d] 220 self._perturbations = np.zeros((params_num, self._n_samples)) 221 offset = 1 #the first sample is the nominal value at x0 222 for p in range(params_num): 223 self._perturbations[p][offset] = + self._perturbation 224 self._perturbations[p][offset+1] = - self._perturbation 225 offset = offset + 2 226 227 # Convert perturbation to absolute values 228 self._nodes = np.array([ nominal_value[p] * np.ones(self._n_samples) for p in vary.keys() ]) 229 offset = 1 #the first sample is the nominal value at x0 230 for p in range(params_num): 231 232 if self.relative_analysis: 233 self._nodes[p][offset] = (1 + self._perturbations[p][offset]) * self._nodes[p][offset] 234 self._nodes[p][offset+1] = (1 + self._perturbations[p][offset+1]) * self._nodes[p][offset+1] 235 else: 236 self._nodes[p][offset] = self._nodes[p][offset] + self._perturbations[p][offset] 237 self._nodes[p][offset+1] = self._nodes[p][offset+1] + self._perturbations[p][offset+1] 238 239 offset = offset + 2 240 241 self.logger.info(f"Generated {offset}/{self._n_samples} samples for the FD scheme") 242 243 # Create perturbed values with correlations 244 # dependent Nodes, where di is the induced movement of the parameter i caused by movement in d 245 #G: [0, -d, d, 0, 0, 0, 0] 246 #C: [0, -di, di, -d, d, 0, 0] 247 #E: [0, -di, di, -di, di, -d, d] 248 if self._is_dependent: 249 250 # Assume permutation [0,1,2] 251 perm = [(i) % params_num for i in range(params_num)] 252 vary_ = {x: vary[x] for i, x in enumerate([list(vary.keys())[i] for i in perm])} 253 254 if self._transformation == "Rosenblatt": 255 self.logger.info("Performing Rosenblatt transformation") 256 257 # Create the dependent distribution 258 mean_, _, cov_ = self.permute_problem(perm, mu, np.zeros(params_num), cov) 259 distribution_dep_ = cp.MvNormal(mean_, cov_) 260 261 # Create the independent distribution 262 params_distribution = [vary_dist for vary_dist in vary_.values()] 263 distribution_ = cp.J(*params_distribution) 264 265 # This assumes that the order of the parameters in distribution and distribution_dep is the same 266 # and the distribution type is cp.Normal 267 for id_v, v in enumerate(vary_): 268 assert(type(vary_[v]) == type(cp.Normal())) 269 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_dep_._parameters['mean'][id_v]) 270 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_[id_v].get_mom_parameters()['shift'][0]) 271 self.logger.debug(f"The independent distribution consists of: {distribution_}") 272 self.logger.debug(f"Using parameter permutation: {list(vary_.keys())}") 273 274 self._nodes_dep = Transformations.rosenblatt(self._nodes, distribution_, distribution_dep_) 275 #self._perturbations_dep = Transformations.rosenblatt(self._perturbations, distribution_, distribution_dep_) 276 elif self._transformation == "Cholesky": 277 self.logger.info("Performing Cholesky transformation") 278 279 _, _, distribution_ = self.permute_problem(perm, np.zeros(params_num), np.zeros(params_num), distribution) 280 self._nodes_dep = Transformations.cholesky(self._nodes, vary_, distribution_) 281 #self._perturbations_dep = Transformations.cholesky(self._perturbations, vary_, distribution_) 282 else: 283 self.logger.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 284 exit() 285 286 return 287 288 289 def is_finite(self): 290 return True 291 292 @property 293 def n_samples(self): 294 """ 295 Number of samples (Ns) of PCE method. 296 - When using pseudo-spectral projection method with tensored 297 quadrature: Ns = (p + 1)**d 298 - When using pseudo-spectral projection method with sparce grid 299 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 300 - When using regression method: Ns = 2*(p + d)!/p!*d! 301 Where: p is the polynomial degree and d is the number of 302 uncertain parameters. 303 304 Ref: Eck et al. 'A guide to uncertainty quantification and 305 sensitivity analysis for cardiovascular applications' [2016]. 306 """ 307 return self._n_samples 308 309 @property 310 def analysis_class(self): 311 """Return a corresponding analysis class. 312 """ 313 from easyvvuq.analysis import FDAnalysis 314 return FDAnalysis 315 316 def __next__(self): 317 if self.count < self._n_samples: #base Train samples used to evaluate the PCE 318 run_dict = {} 319 for i, param_name in enumerate(self.vary.vary_dict): 320 # These are nodes that need to be returned as samples o be used for the model execution, 321 # for the SA in EasyVVUQ we will use only the raw independent nodes 322 if self._is_dependent: 323 # Return transformed nodes reflecting the dependencies 324 run_dict[param_name] = self._nodes_dep[i][self.count] 325 else: 326 run_dict[param_name] = self._nodes[i][self.count] 327 self.count += 1 328 return run_dict 329 else: 330 raise StopIteration
40class FDSampler(BaseSamplingElement, sampler_name="FD_sampler"): 41 def __init__(self, 42 vary=None, 43 distribution=None, 44 perturbation=0.05, 45 nominal_value=None, 46 count=0, 47 relative_analysis=False): 48 """ 49 Create the sampler for the Polynomial Chaos Expansion using 50 pseudo-spectral projection or regression (Point Collocation). 51 52 Parameters 53 ---------- 54 vary: dict or None 55 keys = parameters to be sampled, values = distributions. 56 57 distribution: cp.Distribution or matrix-like 58 Joint distribution specifying dependency between the parameters or 59 correlation matrix of the parameters. Depending on the type of the argument 60 either Rosenblatt or Cholesky transformation will be used to handle the 61 dependent parameters. 62 63 perturbation: float 64 Perturbation of the parameters used in the finite difference scheme 65 66 nominal_value : dict, optional 67 FD sampler perturbs the base value (NV) of the parameters by the value specified. 68 It should be a dict with the keys which are present in vary. 69 In case the nominal_value is None, the mean of the distribution is used (assuming cp.Normal). 70 71 count : int, optional 72 Specified counter for Fast forward, default is 0. 73 74 relative_analysis : bool, optional 75 Default False, the nodes are perturbed by the value specified in perturbation 76 such that n = NV + perturbation. If True, the nodes are perturbed by the relative 77 value specified in perturbation such that n = NB * (1 + perturbation). 78 79 """ 80 # Create and initialize the logger 81 self.logger = logging.getLogger(__name__) 82 self.logger.setLevel(logging.DEBUG) 83 84 # Logger is already configured, remove all handlers 85 if self.logger.hasHandlers(): 86 self.logger.handlers = [] 87 88 formatter = logging.Formatter('%(asctime)s:%(name)s:%(levelname)s:%(message)s') 89 90 file_handler = logging.FileHandler('FD.log') 91 file_handler.setLevel(logging.DEBUG) 92 file_handler.setFormatter(formatter) 93 94 stream_handler = logging.StreamHandler() 95 stream_handler.setFormatter(formatter) 96 97 self.logger.addHandler(file_handler) 98 self.logger.addHandler(stream_handler) 99 100 if vary is None: 101 msg = ("'vary' cannot be None. RandomSampler must be passed a " 102 "dict of the names of the parameters you want to vary, " 103 "and their corresponding distributions.") 104 self.logger.error(msg) 105 raise Exception(msg) 106 if not isinstance(vary, dict): 107 msg = ("'vary' must be a dictionary of the names of the " 108 "parameters you want to vary, and their corresponding " 109 "distributions.") 110 self.logger.error(msg) 111 raise Exception(msg) 112 if len(vary) == 0: 113 msg = "'vary' cannot be empty." 114 self.logger.error(msg) 115 raise Exception(msg) 116 117 self.vary = Vary(vary) 118 119 # List of the probability distributions of uncertain parameters 120 params_distribution = list(vary.values()) 121 params_num = len(params_distribution) 122 123 # Remember whether to add the extra run 124 self.logger.info(f"Performing relative analysis: {relative_analysis}") 125 self.relative_analysis = relative_analysis 126 self._perturbation = perturbation 127 128 # Perturbation of the parameters 129 if nominal_value is None: 130 self.logger.info(f"Performing perturbation of the nodes, base value = mean, with delta = {perturbation}") 131 # Assumes that v is cp.Normal() 132 assert(all([type(v) == type(cp.Normal()) for v in vary.values()])) 133 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 134 else: 135 if (len(nominal_value) != params_num): 136 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 137 self.logger.error(msg) 138 raise ValueError(msg) 139 self.logger.info(f"Performing perturbation of the nodes, base value = {nominal_value}, with delta = {perturbation}") 140 141 # Generate the perturbed values of the parameters for the FD 142 #FD = 0.5*(y_pos/y_base-1)/(delta) + 0.5*(y_neg/y_base - 1)/(-delta) 143 self.generate_nodes(nominal_value, vary, distribution) 144 145 # Fast forward to specified count, if possible 146 self.count = 0 147 if self.count >= self._n_samples: 148 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 149 f"but sampler only has {self.n_samples} samples, therefore" 150 f"this sampler will not provide any more samples.") 151 self.logger.warning(msg) 152 else: 153 for i in range(count): 154 self.__next__() 155 156 157 def permute_problem(self, perm, mean, sigma, cov): 158 # Apply selected permutation to the problem 159 nParams = len(mean) 160 # cyclic permutation i = perm_step, N = nParams 161 # perm = [i i+1 ... N-1 0 1 ... i-1] 162 assert(len(perm) == nParams) 163 # create permutation matrix 164 P = np.eye(nParams)[perm] 165 166 # Model properties that need to be permuted 167 mean_ = P@mean 168 sigma_ = P@sigma 169 cov_ = np.matmul(P, np.matmul(cov, P.transpose())) 170 171 return mean_, sigma_, cov_ 172 173 def generate_nodes(self, nominal_value, vary, distribution): 174 params_num = len(nominal_value) 175 self._n_samples = 2*params_num + 1 176 177 # Multivariate distribution, the behaviour changes based on the 178 # 'distribution' argument, which can be: 179 # None - use default joint 180 # cp.Distribution - use Rosenblatt if the distribution is dependent 181 # matrix-like - use Cholesky 182 self._is_dependent = False 183 self._transformation = None 184 self.distribution_dep = None 185 if 'distributions' in str(type(distribution)): 186 if not distribution.stochastic_dependent: 187 raise ValueError("User provided joint distribution needs to contain dependency between the parameters") 188 if not isinstance(distribution, cp.MvNormal): 189 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 190 if not len(distribution._parameters['mean']) == params_num: 191 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 192 self.logger.info("Using user provided joint distribution with Rosenblatt transformation") 193 194 self._is_dependent = True 195 self._transformation = "Rosenblatt" 196 self.distribution_dep = distribution 197 198 mu = distribution._parameters['mean'] 199 cov = distribution._covariance 200 201 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 202 if not len(distribution) == params_num: 203 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 204 for i in range(params_num): 205 if not distribution[i][i] == 1.0: 206 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 207 self.logger.info("Using user provided correlation matrix for Cholesky transformation") 208 209 self._is_dependent = True 210 self._transformation = "Cholesky" 211 self.distribution_dep = np.array(distribution) 212 elif distribution is None: 213 pass 214 else: 215 raise ValueError("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 216 217 # Set up independent perturbations 218 #G: [0, d, -d, 0, 0, 0, 0] 219 #C: [0, 0, 0, d, -d, 0, 0] 220 #E: [0, 0, 0, 0, 0, d, -d] 221 self._perturbations = np.zeros((params_num, self._n_samples)) 222 offset = 1 #the first sample is the nominal value at x0 223 for p in range(params_num): 224 self._perturbations[p][offset] = + self._perturbation 225 self._perturbations[p][offset+1] = - self._perturbation 226 offset = offset + 2 227 228 # Convert perturbation to absolute values 229 self._nodes = np.array([ nominal_value[p] * np.ones(self._n_samples) for p in vary.keys() ]) 230 offset = 1 #the first sample is the nominal value at x0 231 for p in range(params_num): 232 233 if self.relative_analysis: 234 self._nodes[p][offset] = (1 + self._perturbations[p][offset]) * self._nodes[p][offset] 235 self._nodes[p][offset+1] = (1 + self._perturbations[p][offset+1]) * self._nodes[p][offset+1] 236 else: 237 self._nodes[p][offset] = self._nodes[p][offset] + self._perturbations[p][offset] 238 self._nodes[p][offset+1] = self._nodes[p][offset+1] + self._perturbations[p][offset+1] 239 240 offset = offset + 2 241 242 self.logger.info(f"Generated {offset}/{self._n_samples} samples for the FD scheme") 243 244 # Create perturbed values with correlations 245 # dependent Nodes, where di is the induced movement of the parameter i caused by movement in d 246 #G: [0, -d, d, 0, 0, 0, 0] 247 #C: [0, -di, di, -d, d, 0, 0] 248 #E: [0, -di, di, -di, di, -d, d] 249 if self._is_dependent: 250 251 # Assume permutation [0,1,2] 252 perm = [(i) % params_num for i in range(params_num)] 253 vary_ = {x: vary[x] for i, x in enumerate([list(vary.keys())[i] for i in perm])} 254 255 if self._transformation == "Rosenblatt": 256 self.logger.info("Performing Rosenblatt transformation") 257 258 # Create the dependent distribution 259 mean_, _, cov_ = self.permute_problem(perm, mu, np.zeros(params_num), cov) 260 distribution_dep_ = cp.MvNormal(mean_, cov_) 261 262 # Create the independent distribution 263 params_distribution = [vary_dist for vary_dist in vary_.values()] 264 distribution_ = cp.J(*params_distribution) 265 266 # This assumes that the order of the parameters in distribution and distribution_dep is the same 267 # and the distribution type is cp.Normal 268 for id_v, v in enumerate(vary_): 269 assert(type(vary_[v]) == type(cp.Normal())) 270 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_dep_._parameters['mean'][id_v]) 271 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_[id_v].get_mom_parameters()['shift'][0]) 272 self.logger.debug(f"The independent distribution consists of: {distribution_}") 273 self.logger.debug(f"Using parameter permutation: {list(vary_.keys())}") 274 275 self._nodes_dep = Transformations.rosenblatt(self._nodes, distribution_, distribution_dep_) 276 #self._perturbations_dep = Transformations.rosenblatt(self._perturbations, distribution_, distribution_dep_) 277 elif self._transformation == "Cholesky": 278 self.logger.info("Performing Cholesky transformation") 279 280 _, _, distribution_ = self.permute_problem(perm, np.zeros(params_num), np.zeros(params_num), distribution) 281 self._nodes_dep = Transformations.cholesky(self._nodes, vary_, distribution_) 282 #self._perturbations_dep = Transformations.cholesky(self._perturbations, vary_, distribution_) 283 else: 284 self.logger.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 285 exit() 286 287 return 288 289 290 def is_finite(self): 291 return True 292 293 @property 294 def n_samples(self): 295 """ 296 Number of samples (Ns) of PCE method. 297 - When using pseudo-spectral projection method with tensored 298 quadrature: Ns = (p + 1)**d 299 - When using pseudo-spectral projection method with sparce grid 300 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 301 - When using regression method: Ns = 2*(p + d)!/p!*d! 302 Where: p is the polynomial degree and d is the number of 303 uncertain parameters. 304 305 Ref: Eck et al. 'A guide to uncertainty quantification and 306 sensitivity analysis for cardiovascular applications' [2016]. 307 """ 308 return self._n_samples 309 310 @property 311 def analysis_class(self): 312 """Return a corresponding analysis class. 313 """ 314 from easyvvuq.analysis import FDAnalysis 315 return FDAnalysis 316 317 def __next__(self): 318 if self.count < self._n_samples: #base Train samples used to evaluate the PCE 319 run_dict = {} 320 for i, param_name in enumerate(self.vary.vary_dict): 321 # These are nodes that need to be returned as samples o be used for the model execution, 322 # for the SA in EasyVVUQ we will use only the raw independent nodes 323 if self._is_dependent: 324 # Return transformed nodes reflecting the dependencies 325 run_dict[param_name] = self._nodes_dep[i][self.count] 326 else: 327 run_dict[param_name] = self._nodes[i][self.count] 328 self.count += 1 329 return run_dict 330 else: 331 raise StopIteration
Baseclass for all EasyVVUQ sampling elements.
Attributes
- sampler_name (str): Name of the particular sampler.
FDSampler( vary=None, distribution=None, perturbation=0.05, nominal_value=None, count=0, relative_analysis=False)
41 def __init__(self, 42 vary=None, 43 distribution=None, 44 perturbation=0.05, 45 nominal_value=None, 46 count=0, 47 relative_analysis=False): 48 """ 49 Create the sampler for the Polynomial Chaos Expansion using 50 pseudo-spectral projection or regression (Point Collocation). 51 52 Parameters 53 ---------- 54 vary: dict or None 55 keys = parameters to be sampled, values = distributions. 56 57 distribution: cp.Distribution or matrix-like 58 Joint distribution specifying dependency between the parameters or 59 correlation matrix of the parameters. Depending on the type of the argument 60 either Rosenblatt or Cholesky transformation will be used to handle the 61 dependent parameters. 62 63 perturbation: float 64 Perturbation of the parameters used in the finite difference scheme 65 66 nominal_value : dict, optional 67 FD sampler perturbs the base value (NV) of the parameters by the value specified. 68 It should be a dict with the keys which are present in vary. 69 In case the nominal_value is None, the mean of the distribution is used (assuming cp.Normal). 70 71 count : int, optional 72 Specified counter for Fast forward, default is 0. 73 74 relative_analysis : bool, optional 75 Default False, the nodes are perturbed by the value specified in perturbation 76 such that n = NV + perturbation. If True, the nodes are perturbed by the relative 77 value specified in perturbation such that n = NB * (1 + perturbation). 78 79 """ 80 # Create and initialize the logger 81 self.logger = logging.getLogger(__name__) 82 self.logger.setLevel(logging.DEBUG) 83 84 # Logger is already configured, remove all handlers 85 if self.logger.hasHandlers(): 86 self.logger.handlers = [] 87 88 formatter = logging.Formatter('%(asctime)s:%(name)s:%(levelname)s:%(message)s') 89 90 file_handler = logging.FileHandler('FD.log') 91 file_handler.setLevel(logging.DEBUG) 92 file_handler.setFormatter(formatter) 93 94 stream_handler = logging.StreamHandler() 95 stream_handler.setFormatter(formatter) 96 97 self.logger.addHandler(file_handler) 98 self.logger.addHandler(stream_handler) 99 100 if vary is None: 101 msg = ("'vary' cannot be None. RandomSampler must be passed a " 102 "dict of the names of the parameters you want to vary, " 103 "and their corresponding distributions.") 104 self.logger.error(msg) 105 raise Exception(msg) 106 if not isinstance(vary, dict): 107 msg = ("'vary' must be a dictionary of the names of the " 108 "parameters you want to vary, and their corresponding " 109 "distributions.") 110 self.logger.error(msg) 111 raise Exception(msg) 112 if len(vary) == 0: 113 msg = "'vary' cannot be empty." 114 self.logger.error(msg) 115 raise Exception(msg) 116 117 self.vary = Vary(vary) 118 119 # List of the probability distributions of uncertain parameters 120 params_distribution = list(vary.values()) 121 params_num = len(params_distribution) 122 123 # Remember whether to add the extra run 124 self.logger.info(f"Performing relative analysis: {relative_analysis}") 125 self.relative_analysis = relative_analysis 126 self._perturbation = perturbation 127 128 # Perturbation of the parameters 129 if nominal_value is None: 130 self.logger.info(f"Performing perturbation of the nodes, base value = mean, with delta = {perturbation}") 131 # Assumes that v is cp.Normal() 132 assert(all([type(v) == type(cp.Normal()) for v in vary.values()])) 133 nominal_value = {k: v.get_mom_parameters()['shift'][0] for k,v in vary.items()} #Set nominal_value to the mean_of_the_parameters 134 else: 135 if (len(nominal_value) != params_num): 136 msg = ("'nominal_value' must be a 1D array of the same size as the number of parameters.") 137 self.logger.error(msg) 138 raise ValueError(msg) 139 self.logger.info(f"Performing perturbation of the nodes, base value = {nominal_value}, with delta = {perturbation}") 140 141 # Generate the perturbed values of the parameters for the FD 142 #FD = 0.5*(y_pos/y_base-1)/(delta) + 0.5*(y_neg/y_base - 1)/(-delta) 143 self.generate_nodes(nominal_value, vary, distribution) 144 145 # Fast forward to specified count, if possible 146 self.count = 0 147 if self.count >= self._n_samples: 148 msg = (f"Attempt to start sampler fastforwarded to count {self.count}, " 149 f"but sampler only has {self.n_samples} samples, therefore" 150 f"this sampler will not provide any more samples.") 151 self.logger.warning(msg) 152 else: 153 for i in range(count): 154 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.
- perturbation (float): Perturbation of the parameters used in the finite difference scheme
- nominal_value (dict, optional): FD sampler perturbs the base value (NV) of the parameters by the value specified. It should be a dict with the keys which are present in vary. In case the nominal_value is None, the mean of the distribution is used (assuming cp.Normal).
- count (int, optional): Specified counter for Fast forward, default is 0.
- relative_analysis (bool, optional): Default False, the nodes are perturbed by the value specified in perturbation such that n = NV + perturbation. If True, the nodes are perturbed by the relative value specified in perturbation such that n = NB * (1 + perturbation).
def
permute_problem(self, perm, mean, sigma, cov):
157 def permute_problem(self, perm, mean, sigma, cov): 158 # Apply selected permutation to the problem 159 nParams = len(mean) 160 # cyclic permutation i = perm_step, N = nParams 161 # perm = [i i+1 ... N-1 0 1 ... i-1] 162 assert(len(perm) == nParams) 163 # create permutation matrix 164 P = np.eye(nParams)[perm] 165 166 # Model properties that need to be permuted 167 mean_ = P@mean 168 sigma_ = P@sigma 169 cov_ = np.matmul(P, np.matmul(cov, P.transpose())) 170 171 return mean_, sigma_, cov_
def
generate_nodes(self, nominal_value, vary, distribution):
173 def generate_nodes(self, nominal_value, vary, distribution): 174 params_num = len(nominal_value) 175 self._n_samples = 2*params_num + 1 176 177 # Multivariate distribution, the behaviour changes based on the 178 # 'distribution' argument, which can be: 179 # None - use default joint 180 # cp.Distribution - use Rosenblatt if the distribution is dependent 181 # matrix-like - use Cholesky 182 self._is_dependent = False 183 self._transformation = None 184 self.distribution_dep = None 185 if 'distributions' in str(type(distribution)): 186 if not distribution.stochastic_dependent: 187 raise ValueError("User provided joint distribution needs to contain dependency between the parameters") 188 if not isinstance(distribution, cp.MvNormal): 189 raise ValueError("User provided joint distribution needs to be a cp.MvNormal") 190 if not len(distribution._parameters['mean']) == params_num: 191 raise ValueError("User provided joint distribution does not contain all the parameters listed in vary") 192 self.logger.info("Using user provided joint distribution with Rosenblatt transformation") 193 194 self._is_dependent = True 195 self._transformation = "Rosenblatt" 196 self.distribution_dep = distribution 197 198 mu = distribution._parameters['mean'] 199 cov = distribution._covariance 200 201 elif 'list' in str(type(distribution)) or 'ndarray' in str(type(distribution)): 202 if not len(distribution) == params_num: 203 raise ValueError("User provided correlation matrix does not contain all the parameters listed in vary") 204 for i in range(params_num): 205 if not distribution[i][i] == 1.0: 206 raise ValueError("User provided correlation matrix is not a correlation matrix (diagonal elements are not 1.0)") 207 self.logger.info("Using user provided correlation matrix for Cholesky transformation") 208 209 self._is_dependent = True 210 self._transformation = "Cholesky" 211 self.distribution_dep = np.array(distribution) 212 elif distribution is None: 213 pass 214 else: 215 raise ValueError("Unsupported type of the distribution argument. It should be either cp.Distribution or a matrix-like array") 216 217 # Set up independent perturbations 218 #G: [0, d, -d, 0, 0, 0, 0] 219 #C: [0, 0, 0, d, -d, 0, 0] 220 #E: [0, 0, 0, 0, 0, d, -d] 221 self._perturbations = np.zeros((params_num, self._n_samples)) 222 offset = 1 #the first sample is the nominal value at x0 223 for p in range(params_num): 224 self._perturbations[p][offset] = + self._perturbation 225 self._perturbations[p][offset+1] = - self._perturbation 226 offset = offset + 2 227 228 # Convert perturbation to absolute values 229 self._nodes = np.array([ nominal_value[p] * np.ones(self._n_samples) for p in vary.keys() ]) 230 offset = 1 #the first sample is the nominal value at x0 231 for p in range(params_num): 232 233 if self.relative_analysis: 234 self._nodes[p][offset] = (1 + self._perturbations[p][offset]) * self._nodes[p][offset] 235 self._nodes[p][offset+1] = (1 + self._perturbations[p][offset+1]) * self._nodes[p][offset+1] 236 else: 237 self._nodes[p][offset] = self._nodes[p][offset] + self._perturbations[p][offset] 238 self._nodes[p][offset+1] = self._nodes[p][offset+1] + self._perturbations[p][offset+1] 239 240 offset = offset + 2 241 242 self.logger.info(f"Generated {offset}/{self._n_samples} samples for the FD scheme") 243 244 # Create perturbed values with correlations 245 # dependent Nodes, where di is the induced movement of the parameter i caused by movement in d 246 #G: [0, -d, d, 0, 0, 0, 0] 247 #C: [0, -di, di, -d, d, 0, 0] 248 #E: [0, -di, di, -di, di, -d, d] 249 if self._is_dependent: 250 251 # Assume permutation [0,1,2] 252 perm = [(i) % params_num for i in range(params_num)] 253 vary_ = {x: vary[x] for i, x in enumerate([list(vary.keys())[i] for i in perm])} 254 255 if self._transformation == "Rosenblatt": 256 self.logger.info("Performing Rosenblatt transformation") 257 258 # Create the dependent distribution 259 mean_, _, cov_ = self.permute_problem(perm, mu, np.zeros(params_num), cov) 260 distribution_dep_ = cp.MvNormal(mean_, cov_) 261 262 # Create the independent distribution 263 params_distribution = [vary_dist for vary_dist in vary_.values()] 264 distribution_ = cp.J(*params_distribution) 265 266 # This assumes that the order of the parameters in distribution and distribution_dep is the same 267 # and the distribution type is cp.Normal 268 for id_v, v in enumerate(vary_): 269 assert(type(vary_[v]) == type(cp.Normal())) 270 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_dep_._parameters['mean'][id_v]) 271 assert(vary_[v].get_mom_parameters()['shift'][0] == distribution_[id_v].get_mom_parameters()['shift'][0]) 272 self.logger.debug(f"The independent distribution consists of: {distribution_}") 273 self.logger.debug(f"Using parameter permutation: {list(vary_.keys())}") 274 275 self._nodes_dep = Transformations.rosenblatt(self._nodes, distribution_, distribution_dep_) 276 #self._perturbations_dep = Transformations.rosenblatt(self._perturbations, distribution_, distribution_dep_) 277 elif self._transformation == "Cholesky": 278 self.logger.info("Performing Cholesky transformation") 279 280 _, _, distribution_ = self.permute_problem(perm, np.zeros(params_num), np.zeros(params_num), distribution) 281 self._nodes_dep = Transformations.cholesky(self._nodes, vary_, distribution_) 282 #self._perturbations_dep = Transformations.cholesky(self._perturbations, vary_, distribution_) 283 else: 284 self.logger.critical("Error: How did this happen? We are transforming the nodes but not with Rosenblatt nor Cholesky") 285 exit() 286 287 return
n_samples
293 @property 294 def n_samples(self): 295 """ 296 Number of samples (Ns) of PCE method. 297 - When using pseudo-spectral projection method with tensored 298 quadrature: Ns = (p + 1)**d 299 - When using pseudo-spectral projection method with sparce grid 300 quadratue: Ns = bigO((p + 1)*log(p + 1)**(d-1)) 301 - When using regression method: Ns = 2*(p + d)!/p!*d! 302 Where: p is the polynomial degree and d is the number of 303 uncertain parameters. 304 305 Ref: Eck et al. 'A guide to uncertainty quantification and 306 sensitivity analysis for cardiovascular applications' [2016]. 307 """ 308 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].