easyvvuq.analysis.fd_analysis

  1"""Analysis element for polynomial chaos expansion (PCE). We use ChaosPy
  2under the hood for this functionality.
  3"""
  4import logging
  5import chaospy as cp
  6import numpy as np
  7import numpoly
  8import warnings
  9from easyvvuq import OutputType
 10from .base import BaseAnalysisElement
 11from .results import AnalysisResults
 12from .qmc_analysis import QMCAnalysisResults
 13
 14__author__ = 'Jalal Lakhlili'
 15__license__ = "LGPL"
 16
 17logger = logging.getLogger(__name__)
 18
 19
 20class PCEAnalysisResults(QMCAnalysisResults):
 21    """Analysis results for the FDAnalysis class.
 22    """
 23
 24    def _get_derivatives_first(self, qoi, input_):
 25        """Returns the first order derivative-based index for a given qoi wrt input variable.
 26
 27        Parameters
 28        ----------
 29        qoi : str
 30           Quantity of interest
 31        input_ : str
 32           Input variable
 33
 34        Returns
 35        -------
 36        float
 37            First order derivative-based index.
 38        """
 39
 40        raw_dict = AnalysisResults._keys_to_tuples(self.raw_data['derivatives_first'])
 41        return raw_dict[AnalysisResults._to_tuple(qoi)][input_]
 42
 43    def _get_sobols_first(self, qoi, input_):
 44        """Returns the first order sobol index for a given qoi wrt input variable.
 45
 46        Parameters
 47        ----------
 48        qoi : str
 49           Quantity of interest
 50        input_ : str
 51           Input variable
 52
 53        Returns
 54        -------
 55        float
 56            First order sobol index.
 57        """
 58        raw_dict = AnalysisResults._keys_to_tuples(self.raw_data['sobols_first'])
 59        return raw_dict[AnalysisResults._to_tuple(qoi)][input_]
 60
 61    def _get_sobols_second(self, qoi, input_):
 62        """Returns the second order sobol index for a given qoi wrt input variable.
 63
 64        Parameters
 65        ----------
 66        qoi : str
 67           Quantity of interest
 68        input_ : str
 69           Input variable
 70
 71        Returns
 72        -------
 73        float
 74            Second order sobol index.
 75        """
 76        raw_dict = AnalysisResults._keys_to_tuples(self.raw_data['sobols_second'])
 77        return dict([(in_, raw_dict[AnalysisResults._to_tuple(qoi)][input_][i])
 78                     for i, in_ in enumerate(self.inputs) if in_ != input_])
 79
 80    def _get_sobols_total(self, qoi, input_):
 81        """Returns the total order sobol index for a given qoi wrt input variable.
 82
 83        Parameters
 84        ----------
 85        qoi : str
 86           Quantity of interest
 87        input_ : str
 88           Input variable
 89
 90        Returns
 91        -------
 92        float
 93            Total order sobol index.
 94        """
 95        raw_dict = AnalysisResults._keys_to_tuples(self.raw_data['sobols_total'])
 96        return raw_dict[AnalysisResults._to_tuple(qoi)][input_]
 97
 98    def supported_stats(self):
 99        """Types of statistics supported by the describe method.
100
101        Returns
102        -------
103        list of str
104        """
105        return ['min', 'max', '10%', '90%', '1%', '99%', 'median',
106                'mean', 'var', 'std']
107
108    def _describe(self, qoi, statistic):
109        """Returns descriptive statistics, similar to pandas describe.
110
111        Parameters
112        ----------
113        qoi : str
114            Name of quantity of interest.
115        statistic : str
116            One of 'min', 'max', '10%', '90%', 'median', 'mean', 'var', 'std'
117
118        Returns
119        -------
120        float
121            Value of the requested statistic.
122        """
123        if statistic not in self.supported_stats():
124            raise NotImplementedError
125        if statistic == 'min':
126            return np.array([v.lower[0] for _, v in enumerate(
127                self.raw_data['output_distributions'][qoi])])
128        elif statistic == 'max':
129            return np.array([v.upper[0] for _, v in enumerate(
130                self.raw_data['output_distributions'][qoi])])
131        elif statistic == '1%':
132            return self.raw_data['percentiles'][qoi]['p01']
133        elif statistic == '10%':
134            return self.raw_data['percentiles'][qoi]['p10']
135        elif statistic == '90%':
136            return self.raw_data['percentiles'][qoi]['p90']
137        elif statistic == '99%':
138            return self.raw_data['percentiles'][qoi]['p99']
139        elif statistic == 'median':
140            return self.raw_data['percentiles'][qoi]['p50']
141        else:
142            try:
143                return self.raw_data['statistical_moments'][qoi][statistic]
144            except KeyError:
145                raise NotImplementedError
146
147    def surrogate(self):
148        """Return a PCE surrogate model.
149
150        Returns
151        -------
152        A function that takes a dictionary of parameter - value pairs and returns
153        a dictionary with the results (same output as decoder).
154        """
155        def surrogate_fn(inputs):
156            def swap(x):
157                if len(x) > 1:
158                    return list(x)
159                else:
160                    return x[0]
161            values = np.array([inputs[key] for key in self.inputs])
162            results = dict([(qoi, swap((self.raw_data['fit'][qoi](*values)).T)) for qoi in self.qois])
163            return results
164        return surrogate_fn
165
166    def get_distribution(self, qoi):
167        """Returns a distribution for the given qoi.
168
169        Parameters
170        ----------
171        qoi: str
172            QoI name
173
174        Returns
175        -------
176        A ChaosPy PDF
177        """
178        if qoi not in self.qois:
179            raise RuntimeError('no such quantity of interest - {}'.format(qoi))
180        return self.raw_data['output_distributions'][qoi]
181
182
183class FDAnalysis(BaseAnalysisElement):
184
185    def __init__(self, sampler=None, qoi_cols=None):
186        """Analysis element for polynomial chaos expansion (PCE).
187
188        Parameters
189        ----------
190        sampler : PCESampler
191            Sampler used to initiate the PCE analysis.
192        qoi_cols : list or None
193            Column names for quantities of interest (for which analysis is
194            performed).
195        """
196
197        if sampler is None:
198            msg = 'FD analysis requires a paired sampler to be passed'
199            raise RuntimeError(msg)
200
201        # Flag specifing if we should scale the runs with the nominal base run
202        self.relative_analysis = sampler.relative_analysis
203
204        if qoi_cols is None:
205            raise RuntimeError("Analysis element requires a list of "
206                               "quantities of interest (qoi)")
207
208        self.qoi_cols = qoi_cols
209        self.output_type = OutputType.SUMMARY
210        self.sampler = sampler
211
212    def element_name(self):
213        """Name for this element for logging purposes.
214
215        Returns
216        -------
217        str
218            "FD_Analysis"
219        """
220        return "FD_Analysis"
221
222    def element_version(self):
223        """Version of this element for logging purposes.
224
225        Returns
226        -------
227        str
228            Element version.
229        """
230        return "0.6"
231
232    def analyse(self, data_frame=None):
233        """Perform PCE analysis on input `data_frame`.
234
235        Parameters
236        ----------
237        data_frame : pandas DataFrame
238            Input data for analysis.
239
240        Returns
241        -------
242        PCEAnalysisResults
243            Use it to get the sobol indices and other information.
244        """
245
246        if data_frame is None:
247            raise RuntimeError("Analysis element needs a data frame to "
248                               "analyse")
249        elif data_frame.empty:
250            raise RuntimeError(
251                "No data in data frame passed to analyse element")
252
253        qoi_cols = self.qoi_cols
254        T = len(data_frame[qoi_cols[0]].values[-1])
255
256        results = {'statistical_moments': {k: {'mean':np.zeros(T),
257                                               'var':np.zeros(T),
258                                               'std':np.zeros(T)} for k in qoi_cols},
259                   'percentiles': {k: {'p01': np.zeros(T),
260                                       'p10': np.zeros(T),
261                                       'p50': np.zeros(T),
262                                       'p90': np.zeros(T),
263                                       'p99': np.zeros(T)} for k in qoi_cols},
264                   'sobols_first': {k: {p: np.zeros(T) for p in self.sampler.vary.vary_dict} for k in qoi_cols},
265                   'sobols_second': {k: {p: np.zeros(T) for p in self.sampler.vary.vary_dict} for k in qoi_cols},
266                   'sobols_total': {k: {p: np.zeros(T) for p in self.sampler.vary.vary_dict} for k in qoi_cols},
267                   'correlation_matrices': {k: {} for k in qoi_cols},
268                   'output_distributions': {k: {} for k in qoi_cols},
269                   'fit': {k: cp.polynomial(np.zeros(T)) for k in qoi_cols},
270                   'Fourier_coefficients': {k: {p: np.zeros(T) for p in self.sampler.vary.vary_dict} for k in qoi_cols},
271                   'derivatives_first': {k: {p: np.zeros(T) for p in self.sampler.vary.vary_dict} for k in qoi_cols},
272                   }
273
274        # Get sampler informations
275        nodes = self.sampler._nodes
276        perturbations = self.sampler._perturbations
277        if self.sampler._is_dependent:
278            nodes_dep = self.sampler._nodes_dep
279            #perturbations_dep = self.sampler._perturbations_dep
280            
281
282        for k in qoi_cols:
283
284            base = data_frame[k].values[0]
285            if self.relative_analysis:
286                if np.all(np.array(base) == 0):
287                    warnings.warn(f"Removing QoI {k} from the analysis, contains all zeros", RuntimeWarning)
288                    continue
289                if np.any(np.array(base) == 0):
290                    warnings.warn(f"Removing QoI {k} from the analysis, contains some zeros", RuntimeWarning)
291                    continue
292
293            results['statistical_moments'][k] = {'mean': np.mean(data_frame[k].values, axis=0),
294                                                     'var': np.var(data_frame[k].values, axis=0),
295                                                     'std': np.std(data_frame[k].values, axis=0)}
296
297            # Get the QoI value for the base value of the parameters
298            y_base = data_frame[k].values[0]
299
300            # Compute FD approximation
301            offset = 1
302            for pi, p in enumerate(self.sampler.vary.vary_dict):
303                
304                # assumes ordering of the nodes [0, ..., +delta, -delta, ...]
305                y_pos = data_frame[k].values[offset]
306                y_neg = data_frame[k].values[offset+1]
307                
308                if self.relative_analysis:
309                    d_pos = perturbations[pi][offset]
310                    d_neg = perturbations[pi][offset+1]
311                    #d_pos = nodes[pi][offset]/nodes[pi][0] - 1
312                    #d_neg = nodes[pi][offset+1]/nodes[pi][0] - 1
313
314                    results["derivatives_first"][k][p] = 0.5*(y_pos/y_base-1)/(d_pos) + 0.5*(y_neg/y_base - 1)/(d_neg)
315
316                    # scale the derivatives to the absolute values
317                    x_base = nodes[pi][0] # base value of the parameter
318                    scaling_factor = y_base/x_base
319                    results["derivatives_first"][k][p] *= scaling_factor
320                else:
321                    d_pos = nodes[pi][offset] - nodes[pi][0]
322                    d_neg = nodes[pi][offset+1] - nodes[pi][0]
323
324                    # norm([dg, 0, 0]) = delta_g
325                    results["derivatives_first"][k][p] = 0.5*(y_pos - y_base)/(d_pos) + 0.5*(y_neg - y_base)/(d_neg)
326
327                offset = offset + 2
328
329        return PCEAnalysisResults(raw_data=results, samples=data_frame,
330                                  qois=self.qoi_cols, inputs=list(self.sampler.vary.get_keys()))