easyvvuq.sampling.qmc

  1"""This sampler is meant to be used with the QMC Analysis module.
  2"""
  3
  4import chaospy as cp
  5from SALib.sample import sobol
  6#from SALib.sample import saltelli
  7from .base import BaseSamplingElement, Vary
  8import logging
  9
 10__author__ = "Jalal Lakhlili"
 11__copyright__ = """
 12
 13    Copyright 2018 Robin A. Richardson, David W. Wright
 14
 15    This file is part of EasyVVUQ
 16
 17    EasyVVUQ is free software: you can redistribute it and/or modify
 18    it under the terms of the Lesser GNU General Public License as published by
 19    the Free Software Foundation, either version 3 of the License, or
 20    (at your option) any later version.
 21
 22    EasyVVUQ is distributed in the hope that it will be useful,
 23    but WITHOUT ANY WARRANTY; without even the implied warranty of
 24    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 25    Lesser GNU General Public License for more details.
 26
 27    You should have received a copy of the Lesser GNU General Public License
 28    along with this program.  If not, see <https://www.gnu.org/licenses/>.
 29
 30"""
 31__license__ = "LGPL"
 32
 33
 34class QMCSampler(BaseSamplingElement, sampler_name="QMC_sampler"):
 35    def __init__(self, vary, n_mc_samples, count=0):
 36        """Create a Quasi Monte Carlo sampler.
 37
 38        Parameters
 39        ----------
 40
 41        vary: dict
 42            Expects a dictionary where the keys are variable names
 43            (inputs for your simulation that you want to vary during
 44            sampling) and values are ChaosPy distributions you want to
 45            sample from.
 46        n_mc_samples : int
 47            An estimate for how many samples the monte carlo run will need.
 48        count : int
 49            This is used to resume sampling. It will skip the first
 50            count samples if this parameter is not zero.
 51        """
 52        if not isinstance(vary, dict):
 53            msg = ("'vary' must be a dictionary of the names of the "
 54                   "parameters you want to vary, and their corresponding "
 55                   "distributions.")
 56            raise RuntimeError(msg)
 57        if len(vary) == 0:
 58            msg = "'vary' cannot be empty."
 59            raise RuntimeError(msg)
 60
 61        discrete_input = [isinstance(p, cp.DiscreteUniform) for p in vary.values()]
 62        assert (True in discrete_input) == False, \
 63            "QMCSampler cannot handle DiscreteUniform, use MCSampler instead"
 64
 65        self.vary = Vary(vary)
 66        self.n_mc_samples = n_mc_samples
 67
 68        # List of the probability distributions of uncertain parameters
 69        params_distribution = list(vary.values())
 70
 71        # Multivariate distribution
 72        self.distribution = cp.J(*params_distribution)
 73
 74        # Generate samples
 75        self.n_params = len(vary)
 76
 77        dist_U = []
 78        for i in range(self.n_params):
 79            dist_U.append(cp.Uniform())
 80        dist_U = cp.J(*dist_U)
 81
 82        problem = {
 83            "num_vars": self.n_params,
 84            "names": list(vary.keys()),
 85            "bounds": [[0, 1]] * self.n_params
 86        }
 87
 88        nodes = sobol.sample(problem, n_mc_samples, calc_second_order=False,scramble=True)
 89
 90        self._samples = self.distribution.inv(dist_U.fwd(nodes.transpose()))
 91
 92        self._n_samples = n_mc_samples * (self.n_params + 2)
 93
 94        # Fast forward to specified count, if possible
 95        self.count = 0
 96        if count >= self._n_samples:
 97            msg = (f"Attempt to start sampler fastforwarded to count {count}, "
 98                   f"but sampler only has {self._n_samples} samples, therefore"
 99                   f"this sampler will not provide any more samples.")
100            logging.debug(msg)
101        else:
102            for _ in range(count):
103                self.__next__()
104
105    def is_finite(self):
106        """Can this sampler produce only a finite number of samples."""
107        return True
108
109    @property
110    def n_samples(self):
111        """Returns the number of samples in this sampler.
112
113        Returns
114        -------
115        This computed with the formula (d + 2) * N, where d is the number
116        of uncertain parameters and N is the (estimated) number of samples
117        for the Monte Carlo method.
118        """
119        return self._n_samples
120
121    @property
122    def analysis_class(self):
123        """Return a corresponding analysis class.
124        """
125        from easyvvuq.analysis import QMCAnalysis
126        return QMCAnalysis
127
128    def __next__(self):
129        if self.count < self.n_samples:
130            run_dict = {}
131            i_par = 0
132            for param_name in self.vary.get_keys():
133                run_dict[param_name] = self._samples.T[self.count][i_par]
134                i_par += 1
135            self.count += 1
136            return run_dict
137        else:
138            raise StopIteration