easyvvuq.sampling.transformations

  1import chaospy as cp
  2import numpy as np
  3
  4__author__ = "Juraj Kardos"
  5__copyright__ = """
  6
  7    Copyright 2022 Juraj Kardos
  8
  9    This file is part of EasyVVUQ
 10
 11    EasyVVUQ is free software: you can redistribute it and/or modify
 12    it under the terms of the Lesser GNU General Public License as published by
 13    the Free Software Foundation, either version 3 of the License, or
 14    (at your option) any later version.
 15
 16    EasyVVUQ is distributed in the hope that it will be useful,
 17    but WITHOUT ANY WARRANTY; without even the implied warranty of
 18    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 19    Lesser GNU General Public License for more details.
 20
 21    You should have received a copy of the Lesser GNU General Public License
 22    along with this program.  If not, see <https://www.gnu.org/licenses/>.
 23
 24"""
 25__license__ = "LGPL"
 26
 27
 28class Transformations:
 29    def __init__(self):
 30        pass
 31
 32    # Applies Rosenblatt transformation
 33    # to the independent nodes.
 34    # Returns: The transformed nodes or transformed (weights, nodes)
 35    # Args:
 36    #   @Nodes - Independent nodes to be transformed
 37    #   @Distribution - PDF of the independent nodes
 38    #   @Distribution_dep - PDF of the correlated nodes
 39    #   @regression - see PCESampler(regression) parameter
 40    @staticmethod
 41    def rosenblatt(nodes, distribution, distribution_dep, regression=True):
 42
 43        # Input nodes are expected to be in sape (ndim x nsamples),
 44        # but user might have have provided the transposed array
 45        do_transpose = False
 46        if not nodes.shape[0] == len(distribution):
 47            do_transpose = True
 48            nodes = nodes.T
 49            if not nodes.shape[0] == len(distribution):
 50                raise ValueError("Input nodes have wrong shape.")
 51        
 52        transformed_nodes = []
 53
 54        transformed_nodes = distribution_dep.inv(distribution.fwd(nodes))
 55
 56        if do_transpose:
 57            transformed_nodes = transformed_nodes.T
 58            nodes = nodes.T # need to transpose back (args. are passed by reference)
 59
 60        # Transform node weights in the pseudo-spectral method
 61        if not regression:
 62            # The transformed weights are not used
 63            transformed_weights = None
 64            # TODO: need to add weights argument
 65            # transformed_weights = weights * distribution_dep.pdf(transformed_nodes)/distribution.pdf(nodes)
 66            return (transformed_weights, transformed_nodes)
 67
 68        return transformed_nodes
 69
 70    # Applies Cholesky transformation
 71    # to the independent nodes.
 72    # Returns: The transformed nodes or transformed (weights, nodes)
 73    # Args:
 74    #   @Nodes - Independent nodes to be transformed
 75    #   @vary - Vary object containing the PDF of the parameters
 76    #   @correlation - correlation matrix
 77    #   @regression - see PCESampler(regression) parameter
 78    @staticmethod
 79    def cholesky(nodes, vary, correlation, regression=True):
 80
 81        if str(type(vary)) == "<class 'dict'>":
 82            items = vary.items() # simple dictionary
 83        else:
 84            items = vary.get_items() # Vary object
 85        nparams = len(items)
 86
 87        # Input nodes are expected to be in shape (ndim x nsamples),
 88        # but user might have have provided the transposed array
 89        do_transpose = False
 90        if not nodes.shape[0] == nparams:
 91            do_transpose = True
 92            nodes = nodes.T
 93            if not nodes.shape[0] == nparams:
 94                raise ValueError("Input nodes have wrong shape.")
 95
 96        # Shift and stretch the nodes to a unit normal distribution
 97        # Until now we have samples from a general non-unit normal distribution
 98        nodes_unit = np.zeros(nodes.shape)
 99        for i, (param, distribution) in enumerate(items):
100            if type(distribution).__name__ == "Uniform":
101                a = distribution._parameters['lower'] #lower
102                b = distribution._parameters['upper'] #upper
103                nodes_unit[i] = (nodes[i] - a) / (b-a)
104            elif type(distribution).__name__ == "Normal":
105                a = distribution._parameters['shift'] #mu
106                b = distribution._parameters['scale'] #sigma
107                nodes_unit[i] = (nodes[i] - a) / b
108
109        transformed_nodes = []
110        L = np.linalg.cholesky(correlation)
111        transformed_nodes = np.matmul(L, nodes_unit)
112
113        # Shift and stretch the transformed nodes to the target distr.
114        # Until now we had samples from unit uniform (or normal) distributions
115        for i, (param, distribution) in enumerate(items):
116            if type(distribution).__name__ == "Uniform":
117                a = distribution._parameters['lower'] #lower
118                b = distribution._parameters['upper'] #upper
119                transformed_nodes[i] = a + (b-a)*transformed_nodes[i]
120            elif type(distribution).__name__ == "Normal":
121                a = distribution._parameters['shift'] #mu
122                b = distribution._parameters['scale'] #sigma
123                transformed_nodes[i] = a + b*transformed_nodes[i]
124
125        if do_transpose:
126            transformed_nodes = transformed_nodes.T
127            nodes = nodes.T # need to transpose back (args. are passed by reference)
128
129        # Tested & implemented only with the point collocation!
130        # For spectral projection we need to work also with
131        # the node weights, which requires some additional care
132        # assert(regression)
133        if not regression:
134            transformed_weights = None
135            return (transformed_weights, transformed_nodes)
136
137        return transformed_nodes