-
Notifications
You must be signed in to change notification settings - Fork 26
/
mc_sampler.py
129 lines (107 loc) · 4.48 KB
/
mc_sampler.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
from . import RandomSampler
from .base import Vary
import logging
import numpy as np
import chaospy as cp
__author__ = "Wouter Edeling"
"""
Copyright 2018 Robin A. Richardson, David W. Wright
This file is part of EasyVVUQ
EasyVVUQ is free software: you can redistribute it and/or modify
it under the terms of the Lesser GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
EasyVVUQ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Lesser GNU General Public License for more details.
You should have received a copy of the Lesser GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
__license__ = "LGPL"
class MCSampler(RandomSampler, sampler_name='mc_sampler'):
"""
This is a Monte Carlo sampler, used to compute the Sobol indices, mean
and variance of the different QoIs.
"""
def __init__(self, vary, n_mc_samples, **kwargs):
"""
Parameters
----------
+ vary : a dictionary of chaospy input distributions
+ n_mc_samples : the number of MC samples. The total number of MC samples
required to compute the Sobol indices using a Saltelli sampling plan
will be n_mc_samples * (n_params + 1), where n_params is the number of
uncertain parameters in vary.
Returns
-------
None.
"""
super().__init__(vary=vary, count=0, max_num=n_mc_samples, **kwargs)
# the number of uncertain inputs
self.n_params = len(vary)
# the number of MC samples, for each of the n_params + 2 input matrices
self.n_mc_samples = n_mc_samples
# joint distribution
self.joint = cp.J(*list(vary.values()))
# create the Saltelli sampling plan
self.saltelli(n_mc_samples)
def __next__(self):
if self.is_finite():
if self.count >= self.max_num:
raise StopIteration
run_dict = {}
for idx, param_name in enumerate(self.vary.get_keys()):
run_dict[param_name] = self.xi_mc[self.count][idx]
self.count += 1
return run_dict
def saltelli(self, n_mc):
"""
Generates a Saltelli sampling plan of n_mc*(n_params + 2) input samples
needed to compute the Sobol indices. Stored in xi_mc.
Method: A. Saltelli, Making best use of model evaluations to compute
sensitivity indices, Computer Physics Communications, 2002.
Parameters
----------
n_mc : the number of Monte Carlo samples per input matrix. The total
number of samples is n_mc*(n_params + 2)
Returns
-------
None.
"""
logging.debug('Drawing input samples for Sobol index computation.')
# the number of MC samples required to compute the Sobol indices
self.max_num = n_mc * (self.n_params + 2)
logging.debug('Generating {} input samples spread over {} sample matrices.'.format(
self.max_num, self.n_params + 2))
# Matrix M1, the sample matrix
M_1 = self.joint.sample(n_mc).T
# Matrix M2, the resample matrix (see reference above)
M_2 = self.joint.sample(n_mc).T
# array which contains all samples
self.xi_mc = np.zeros([self.max_num, self.n_params])
# The order in which the inputs samples must be stored is
# [M2_1 N1_1, ..., Nd_1, M1_1, M2_2, N1_2, ...Nd_2, M1_2, M1_3 etc]
# number of different sampling matrices
step = self.n_params + 2
# store M2 first, with entries separated by step places
if M_2.ndim == 1:
M_2 = M_2.reshape([-1, 1])
self.xi_mc[0:self.max_num:step] = M_2
# store M1 entries last
if M_1.ndim == 1:
M_1 = M_1.reshape([-1, 1])
self.xi_mc[(step - 1):self.max_num:step] = M_1
# store N_i entries between M2 and M1
for i in range(self.n_params):
N_i = np.array(M_2)
# N_i = M2 with i-th colum from M1
N_i[:, i] = M_1[:, i]
self.xi_mc[(i + 1):self.max_num:step] = N_i
logging.debug('Done.')
@property
def analysis_class(self):
"""Return a corresponding analysis class.
"""
from easyvvuq.analysis import QMCAnalysis
return QMCAnalysis