easyvvuq.sampling.mc_sampler

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