easyvvuq.sampling.stochastic_collocation

  1from .base import BaseSamplingElement, Vary
  2import chaospy as cp
  3import numpy as np
  4import pickle
  5from itertools import product, chain
  6import logging
  7
  8__author__ = "Wouter Edeling"
  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 SCSampler(BaseSamplingElement, sampler_name="sc_sampler"):
 33    """
 34    Stochastic Collocation sampler
 35    """
 36
 37    def __init__(self,
 38                 vary=None,
 39                 polynomial_order=4,
 40                 quadrature_rule="G",
 41                 count=0,
 42                 growth=False,
 43                 sparse=False,
 44                 midpoint_level1=True,
 45                 dimension_adaptive=False):
 46        """
 47        Create the sampler for the Stochastic Collocation method.
 48
 49        Parameters
 50        ----------
 51        vary: dict or None
 52            keys = parameters to be sampled, values = distributions.
 53
 54        polynomial_order : int, optional
 55            The polynomial order, default is 4.
 56
 57        quadrature_rule : char, optional
 58            The quadrature method, default is Gaussian "G".
 59
 60        growth: bool, optional
 61             Sets the growth rule to exponential for Clenshaw Curtis quadrature,
 62             which makes it nested, and therefore more efficient for sparse grids.
 63             Default is False.
 64
 65        sparse : bool, optional
 66            If True use sparse grid instead of normal tensor product grid,
 67            default is False.
 68
 69        midpoint_level1 : bool, optional
 70            Used only for sparse=True (or for cp.DiscreteUniform distribution).
 71            Determines how many points the 1st level of a sparse grid will have.
 72            If True, order 0 quadrature will be generated,
 73            default is True.
 74
 75        dimension_adaptive : bool, optional
 76            Determines wether to use an insotropic sparse grid, or to
 77            adapt the levels in the sparse grid based on a hierachical
 78            error measure, default is False.
 79        """
 80
 81        self.vary = Vary(vary)
 82
 83        # List of the probability distributions of uncertain parameters
 84        self.params_distribution = list(self.vary.get_values())
 85        # N = number of uncertain parameters
 86        self.N = len(self.params_distribution)
 87
 88        logging.debug("param dist {}".format(self.params_distribution))
 89
 90
 91        self.joint_dist = cp.J(*self.params_distribution)
 92
 93        # The quadrature information: order, rule and sparsity
 94        if isinstance(polynomial_order, int):
 95            logging.debug('Received integer polynomial order, assuming isotropic grid')
 96            self.polynomial_order = [polynomial_order for i in range(self.N)]
 97        else:
 98            self.polynomial_order = polynomial_order
 99
100        self.quad_rule = quadrature_rule
101        self.sparse = sparse
102        self.midpoint_level1 = midpoint_level1
103        self.dimension_adaptive = dimension_adaptive
104        self.nadaptations = 0
105        self.growth = growth
106        self.check_max_quad_level()
107
108        # determine if a nested sparse grid is used
109        if self.sparse is True and self.growth is True and \
110           (self.quad_rule == "C" or self.quad_rule == "clenshaw_curtis"):
111            self.nested = True
112        elif self.sparse is True and self.growth is False and self.quad_rule == "gauss_patterson":
113            self.nested = True
114        elif self.sparse is True and self.growth is True and self.quad_rule == "newton_cotes":
115            self.nested = True
116        else:
117            self.nested = False
118
119        # L = level of (sparse) grid
120        self.L = np.max(self.polynomial_order)
121
122        # compute the 1D collocation points (and quad weights)
123        self.compute_1D_points_weights(self.L, self.N)
124
125        # compute N-dimensional collocation points
126        if not self.sparse:
127
128            # generate collocation as a standard tensor product
129            l_norm = np.array([self.polynomial_order])
130            self.xi_d = self.generate_grid(l_norm)
131
132        else:
133            self.l_norm = self.compute_sparse_multi_idx(self.L, self.N)
134            # create sparse grid of dimension N and level q using the 1d
135            # rules in self.xi_1d
136            self.xi_d = self.generate_grid(self.l_norm)
137
138        self._n_samples = self.xi_d.shape[0]
139
140        self.count = 0
141
142    @property
143    def analysis_class(self):
144        """Return a corresponding analysis class.
145        """
146        from easyvvuq.analysis import SCAnalysis
147        return SCAnalysis
148
149    def compute_1D_points_weights(self, L, N):
150        """
151        Computes 1D collocation points and quad weights,
152        and stores this in self.xi_1d, self.wi_1d.
153
154        Parameters
155        ----------
156        L : (int) the max polynomial order of the (sparse) grid
157        N : (int) the number of uncertain parameters
158
159        Returns
160        -------
161        None.
162
163        """
164        # for every dimension (parameter), create a hierachy of 1D
165        # quadrature rules of increasing order
166        self.xi_1d = [{} for n in range(N)]
167        self.wi_1d = [{} for n in range(N)]
168
169        if self.sparse:
170
171            # if level one of the sparse grid is a midpoint rule, generate
172            # the quadrature with order 0 (1 quad point). Else set order at
173            # level 1 to 1
174            if self.midpoint_level1:
175                j = 0
176            else:
177                j = 1
178
179            for n in range(N):
180                # check if input is discrete uniform, in which case the
181                # rule and growth flag must be modified
182                if isinstance(self.params_distribution[n], cp.DiscreteUniform):
183                    rule = "discrete"
184                else:
185                    rule = self.quad_rule
186
187                for i in range(L):
188                    xi_i, wi_i = cp.generate_quadrature(i + j,
189                                                        self.params_distribution[n],
190                                                        rule=rule,
191                                                        growth=self.growth)
192                    self.xi_1d[n][i + 1] = xi_i[0]
193                    self.wi_1d[n][i + 1] = wi_i
194        else:
195            for n in range(N):
196                # check if input is discrete uniform, in which case the
197                # rule flag must be modified
198                if isinstance(self.params_distribution[n], cp.DiscreteUniform):
199                    rule = "discrete"
200                else:
201                    rule = self.quad_rule
202
203                xi_i, wi_i = cp.generate_quadrature(self.polynomial_order[n],
204                                                    self.params_distribution[n],
205                                                    rule=rule,
206                                                    growth=self.growth) 
207
208                self.xi_1d[n][self.polynomial_order[n]] = xi_i[0]
209                self.wi_1d[n][self.polynomial_order[n]] = wi_i
210
211    def check_max_quad_level(self):
212        """
213
214        If a discrete variable is specified, there is the possibility of
215        non unique collocation points if the quadrature order is high enough.
216        This subroutine prevents that.
217
218        NOTE: Only detects cp.DiscreteUniform thus far
219
220        The max quad orders are stores in self.max_quad_order
221
222        Returns
223        -------
224        None
225
226        """
227        # assume no maximum by default
228        self.max_level = np.ones(self.N) * 1000
229        for n in range(self.N):
230
231            # if a discrete uniform is specified check max order
232            if isinstance(self.params_distribution[n], cp.DiscreteUniform):
233
234                #TODO: it is assumed that self.sparse=True, but this assumption
235                #does not have to hold here!!!
236
237                # if level one of the sparse grid is a midpoint rule, generate
238                # the quadrature with order 0 (1 quad point). Else set order at
239                # level 1 to 1
240                if self.midpoint_level1:
241                    j = 0
242                else:
243                    j = 1
244
245                number_of_points = 0
246                for order in range(1000):
247                    xi_i, wi_i = cp.generate_quadrature(order + j,
248                                                        self.params_distribution[n],
249                                                        growth=self.growth)
250                    # if the quadrature points no longer grow with the quad order,
251                    # then the max order has been reached
252                    if xi_i.size == number_of_points:
253                        break
254                    number_of_points = xi_i.size
255
256                logging.debug("Input %d is discrete, setting max quadrature order to %d"
257                              % (n, order - 1))
258                # level 1 = order 0 etc
259                self.max_level[n] = order
260
261    def next_level_sparse_grid(self):
262        """
263        Adds the points of the next level for isotropic hierarchical sparse grids.
264
265        Returns
266        -------
267        None.
268
269        """
270
271        if self.nested is False:
272            print('Only works for nested sparse grids')
273            return
274
275        self.look_ahead(self.l_norm)
276        self.l_norm = np.append(self.l_norm, self.admissible_idx, axis=0)
277
278    def look_ahead(self, current_multi_idx):
279        """
280        The look-ahead step in (dimension-adaptive) sparse grid sampling. Allows
281        for anisotropic sampling plans.
282
283        Computes the admissible forward neighbors with respect to the current level
284        multi-indices. The admissible level indices l are added to self.admissible_idx.
285        The code will be evaluated next iteration at the new collocation points
286        corresponding to the levels in admissble_idx.
287
288        Source: Gerstner, Griebel, "Numerical integration using sparse grids"
289
290        Parameters
291        ----------
292        current_multi_idx : array of the levels in the current iteration
293        of the sparse grid.
294
295        Returns
296        -------
297        None.
298
299        """
300
301        # compute all forward neighbors for every l in current_multi_idx
302        forward_neighbor = []
303        e_n = np.eye(self.N, dtype=int)
304        for l in current_multi_idx:
305            for n in range(self.N):
306                # the forward neighbor is l plus a unit vector
307                forward_neighbor.append(l + e_n[n])
308
309        # remove duplicates
310        forward_neighbor = np.unique(np.array(forward_neighbor), axis=0)
311        # remove those which are already in the grid
312        forward_neighbor = setdiff2d(forward_neighbor, current_multi_idx)
313        # make sure the final candidates are admissible (all backward neighbors
314        # must be in the current multi indices)
315        logging.debug('Computing admissible levels...')
316        admissible_idx = []
317        for l in forward_neighbor:
318            admissible = True
319            for n in range(self.N):
320                backward_neighbor = l - e_n[n]
321                # find backward_neighbor in current_multi_idx
322                idx = np.where((backward_neighbor == current_multi_idx).all(axis=1))[0]
323                # if backward neighbor is not in the current index set
324                # and it is 'on the interior' (contains no 0): not admissible
325                if idx.size == 0 and backward_neighbor[n] != 0:
326                    admissible = False
327                    break
328            # if all backward neighbors are in the current index set: l is admissible
329            if admissible:
330                admissible_idx.append(l)
331        logging.debug('done')
332
333        self.admissible_idx = np.array(admissible_idx)
334        # make sure that all entries of each index are <= the max quadrature order
335        # The max quad order can be low for discrete input variables
336        idx = np.where((self.admissible_idx <= self.max_level).all(axis=1))[0]
337        self.admissible_idx = self.admissible_idx[idx]
338        logging.debug('Admissible multi-indices:\n%s', self.admissible_idx)
339
340        # determine the maximum level L of the new index set L = |l| - N + 1
341        # self.L = np.max(np.sum(self.admissible_idx, axis=1) - self.N + 1)
342        self.L = np.max(self.admissible_idx)
343        # recompute the 1D weights and collocation points
344        self.compute_1D_points_weights(self.L, self.N)
345        # compute collocation grid based on the admissible level indices
346        admissible_grid = self.generate_grid(self.admissible_idx)
347        # remove collocation points which have already been computed
348        if not hasattr(self, 'xi_d'):
349            self.xi_d = self.generate_grid(self.admissible_idx)
350            self._n_samples = self.xi_d.shape[0]
351        new_points = setdiff2d(admissible_grid, self.xi_d)
352
353        logging.debug('%d new points added' % new_points.shape[0])
354
355        # keep track of the number of points added per iteration
356        if not hasattr(self, 'n_new_points'):
357            self.n_new_points = []
358        self.n_new_points.append(new_points.shape[0])
359
360        # update the number of samples
361        self._n_samples += new_points.shape[0]
362
363        # update the N-dimensional sparse grid if unsampled points are added
364        if new_points.shape[0] > 0:
365            self.xi_d = np.concatenate((self.xi_d, new_points))
366
367        # count the number of times the dimensions were adapted
368        self.nadaptations += 1
369
370    def is_finite(self):
371        return True
372
373    @property
374    def n_samples(self):
375        """
376        Number of samples (Ns) of SC method.
377        - When using tensor quadrature: Ns = (p + 1)**d
378        - When using sparid: Ns = bigO((p + 1)*log(p + 1)**(d-1))
379        Where: p is the polynomial degree and d is the number of
380        uncertain parameters.
381
382        Ref: Eck et al. 'A guide to uncertainty quantification and
383        sensitivity analysis for cardiovascular applications' [2016].
384        """
385        return self._n_samples
386
387    # SC collocations points are not random, generate_runs simply returns
388    # one collocation point from the tensor product after the other
389    def __next__(self):
390        if self.count < self._n_samples:
391            run_dict = {}
392            i_par = 0
393            for param_name in self.vary.get_keys():
394                # the current input parameter
395                current_param = self.xi_d[self.count][i_par]
396                # all parameters self.xi_d will store floats. If current param is
397                # DiscreteUniform, convert e.g. 2.0 to 2 before running
398                # the simulation.
399                if isinstance(self.params_distribution[i_par], cp.DiscreteUniform):
400                    current_param = int(current_param)
401                run_dict[param_name] = current_param
402                i_par += 1
403            self.count += 1
404            return run_dict
405        else:
406            raise StopIteration
407
408    def save_state(self, filename):
409        logging.debug("Saving sampler state to %s" % filename)
410        file = open(filename, 'wb')
411        pickle.dump(self.__dict__, file)
412        file.close()
413
414    def load_state(self, filename):
415        logging.debug("Loading sampler state from %s" % filename)
416        file = open(filename, 'rb')
417        self.__dict__ = pickle.load(file)
418        file.close()
419
420    """
421    =========================
422    (SPARSE) GRID SUBROUTINES
423    =========================
424    """
425
426    def generate_grid(self, l_norm):
427
428        dimensions = range(self.N)
429        H_L_N = []
430        # loop over all multi indices i
431        for l in l_norm:
432            # compute the tensor product of nodes indexed by i
433            X_l = [self.xi_1d[n][l[n]] for n in dimensions]
434            H_L_N.append(list(product(*X_l)))
435        # flatten the list of lists
436        H_L_N = np.array(list(chain(*H_L_N)))
437        # return unique nodes
438        return np.unique(H_L_N, axis=0)
439
440    # L : (int) max polynomial order
441    # N : (int) the number of uncertain parameters
442    def compute_sparse_multi_idx(self, L, N):
443        """
444        computes all N dimensional multi-indices l = (l1,...,lN) such that
445        |l| <= L + N - 1, i.e. a simplex set:
446        3    *
447        2    *    *          (L=3 and N=2)
448        1    *    *    *
449              1    2    3
450        Here |l| is the internal sum of i (l1+...+lN)
451        """
452        # old implementation: does not scale well
453        # P = np.array(list(product(range(1, L + 1), repeat=N)))
454        # multi_idx = P[np.where(np.sum(P, axis=1) <= L + N - 1)[0]]
455
456        # use the look_ahead subroutine to build an isotropic sparse grid (faster)
457        multi_idx = np.array([np.ones(self.N, dtype='int')])
458        for l in range(self.L - 1):
459            self.look_ahead(multi_idx)
460            # accept all admissible indices to build an isotropic grid
461            multi_idx = np.append(multi_idx, self.admissible_idx, axis=0)
462
463        return multi_idx
464
465
466def setdiff2d(X, Y):
467    """
468    Computes the difference of two 2D arrays X and Y
469
470    Parameters
471    ----------
472    X : 2D numpy array
473    Y : 2D numpy array
474
475    Returns
476    -------
477    The difference X \\ Y as a 2D array
478
479    """
480    diff = set(map(tuple, X)) - set(map(tuple, Y))
481    return np.array(list(diff))