added axes scales

src/uncertainpy/_version.py

@@ -1 +1 @@
-__version__ = "1.2.3"
+__version__ = "1.2.3"

src/uncertainpy/core/uncertainty_calculations.py

@@ -14,6 +14,8 @@
 from ..utils.logger import get_logger
 
 
+
+
 class UncertaintyCalculations(ParameterBase):
     """
     Perform the calculations for the uncertainty quantification and
@@ -139,7 +141,33 @@ def parameters(self, new_parameters):
 
         self.runmodel.parameters = self.parameters
 
-
+    def calc_sens(self, dim, poly, uhat):
+        variance = np.sum(uhat[1:]**2)
+        main = np.zeros(dim)
+        total = np.zeros(dim)
+        #Compute sens
+        for idx, pol in enumerate(poly[1:]):
+            add_total = [False,]*dim
+            add_main = [True,]*dim
+            for term in pol.exponents:
+                for var in range(dim):
+                    if term[var] > 0:
+                        add_total[var] = True
+                        add_main[var] = add_main[var] and True
+                        add_main[0:var] = [False,]*var
+                        try:
+                            add_main[var+1::] = [False,]*(dim-var)
+                        except IndexError:
+                            pass
+            for var in range(dim):
+                if add_main[var]:
+                    main[var] += uhat[idx+1]**2
+                if add_total[var]:
+                    total[var] += uhat[idx+1]**2
+
+        main = main/variance
+        total = total/variance
+        return main, total
 
 
     def convert_uncertain_parameters(self, uncertain_parameters=None):
@@ -392,7 +420,54 @@ def create_masked_nodes_weights(self, data, feature, nodes, weights):
         return masked_evaluations, mask, masked_nodes, masked_weights
 
 
+    def get_UQ_samples(self, method, uncertain_parameters, polynomial_order, quadrature_order=None, seed=None):
+        if method not in ["spectral", "collocation", "MC"]:
+            raise ValueError("The UQ method needs to be provided (spectral, collocation or MC)")
+        uncertain_parameters = self.convert_uncertain_parameters(uncertain_parameters)
+        distribution = self.create_distribution(uncertain_parameters=uncertain_parameters)
+        if method == "spectral":
+            if quadrature_order is None:
+                raise ValueError("The quadrature_order must be an integer")
+            P = cp.orth_ttr(polynomial_order, distribution)
+            nodes, weights = cp.generate_quadrature(quadrature_order,
+                                                    distribution,
+                                                    rule="J",
+                                                    sparse=True)
+            return uncertain_parameters, distribution, P, nodes, weights
+        elif method == "collocation":
+            P = cp.orth_ttr(polynomial_order, distribution, normed=True)
+            nr_collocation_nodes = quadrature_order
+            if nr_collocation_nodes is None:
+                nr_collocation_nodes = 2 * len(P) + 2
+            nodes = distribution.sample(nr_collocation_nodes, "M")
+            return uncertain_parameters, distribution, P, nodes, None
+        elif method == "MC":
+            nr_samples = polynomial_order
+            if nr_samples is None:
+                raise ValueError("The nr_samples must be an integer")
+
+            if seed is not None:
+                np.random.seed(seed)
+
+            problem = {
+                "num_vars": len(uncertain_parameters),
+                "names": uncertain_parameters,
+                "bounds": [[0,1]]*len(uncertain_parameters)
+            }
+
+            # Create the Multivariate normal distribution
+            dist_R = []
+            for parameter in uncertain_parameters:
+                dist_R.append(cp.Uniform())
+
+            dist_R = cp.J(*dist_R)
+
+            nr_sobol_samples = int(np.round(nr_samples/2.))
+
+            nodes_R = saltelli.sample(problem, nr_sobol_samples, calc_second_order=False)
 
+            nodes = distribution.inv(dist_R.fwd(nodes_R.transpose()))
+            return uncertain_parameters, distribution, None, nodes, nr_sobol_samples
 
 
     def create_PCE_spectral(self,
@@ -478,24 +553,14 @@ def create_PCE_spectral(self,
         uncertainpy.Data
         uncertainpy.Parameters
         """
-        uncertain_parameters = self.convert_uncertain_parameters(uncertain_parameters)
-
-        distribution = self.create_distribution(uncertain_parameters=uncertain_parameters)
-
-        P = cp.orth_ttr(polynomial_order, distribution)
-
         if quadrature_order is None:
             quadrature_order = polynomial_order + 2
-
-
-        nodes, weights = cp.generate_quadrature(quadrature_order,
-                                                distribution,
-                                                rule="J",
-                                                sparse=True)
+        uncertain_parameters, distribution, P, nodes, weights = self.get_UQ_samples("spectral", uncertain_parameters, polynomial_order, quadrature_order=quadrature_order) 
 
         # Running the model
         data = self.runmodel.run(nodes, uncertain_parameters)
 
+        # TODO: import data here instead of running
 
         data.method = "polynomial chaos expansion with the pseudo-spectral method. polynomial_order={}, quadrature_order={}".format(polynomial_order, quadrature_order)
 
@@ -614,20 +679,13 @@ def create_PCE_collocation(self,
         uncertainpy.Parameters
         """
 
-        uncertain_parameters = self.convert_uncertain_parameters(uncertain_parameters)
-
-        distribution = self.create_distribution(uncertain_parameters=uncertain_parameters)
-
-        P = cp.orth_ttr(polynomial_order, distribution)
-        if nr_collocation_nodes is None:
-            nr_collocation_nodes = 2*len(P) + 2
-
-        nodes = distribution.sample(nr_collocation_nodes, "M")
-
+        uncertain_parameters, distribution, P, nodes, _ = self.get_UQ_samples("collocation", uncertain_parameters, polynomial_order, quadrature_order=nr_collocation_nodes) 
 
         # Running the model
         data = self.runmodel.run(nodes, uncertain_parameters)
 
+        # TODO: import data here instead of running
+
         data.method = "polynomial chaos expansion with point collocation. polynomial_order={}, nr_collocation_nodes={}".format(polynomial_order, nr_collocation_nodes)
 
         logger = get_logger(self)
@@ -643,8 +701,8 @@ def create_PCE_collocation(self,
             masked_evaluations, mask, masked_nodes = self.create_masked_nodes(data, feature, nodes)
 
             if (np.all(mask) or allow_incomplete) and sum(mask) > 0:
-                U_hat[feature] = cp.fit_regression(P, masked_nodes,
-                                                   masked_evaluations)
+                U_hat[feature], approx[feature] = cp.fit_regression(P, masked_nodes,
+                                                   masked_evaluations, retall=1)
             elif not allow_incomplete:
                 logger.warning("{}: not all parameter combinations give results.".format(feature) +
                                " No uncertainty quantification is performed since allow_incomplete=False")
@@ -656,7 +714,7 @@ def create_PCE_collocation(self,
             if not np.all(mask):
                 data.incomplete.append(feature)
 
-        return U_hat, distribution, data
+        return U_hat, approx, distribution, data, P
 
 
     def create_PCE_spectral_rosenblatt(self,
@@ -966,7 +1024,7 @@ def create_PCE_collocation_rosenblatt(self,
         return U_hat, dist_R, data
 
 
-    def analyse_PCE(self, U_hat, distribution, data, nr_samples=10**4):
+    def analyse_PCE(self, U_hat, approx, distribution, data, P, polynomial_order, nr_samples=10**4, save_samples=False):
         """
         Calculate the statistical metrics from the polynomial chaos
         approximation.
@@ -985,6 +1043,8 @@ def analyse_PCE(self, U_hat, distribution, data, nr_samples=10**4):
             Number of samples for the Monte Carlo sampling of the polynomial
             chaos approximation.
             Default is 10**4.
+        save_samples: bool, optional
+            Save samples in feature data to, for example, plot PDFs later.
 
         Returns
         -------
@@ -1024,26 +1084,33 @@ def analyse_PCE(self, U_hat, distribution, data, nr_samples=10**4):
             logger.info("Only 1 uncertain parameter. Sensitivities are not calculated")
 
         U_mc = {}
+        samples = distribution.sample(nr_samples, "M")
         for feature in tqdm(data,
                             desc="Calculating statistics from PCE",
                             total=len(data)):
             if feature in U_hat:
                 data[feature].mean = cp.E(U_hat[feature], distribution)
                 data[feature].variance = cp.Var(U_hat[feature], distribution)
 
-                samples = distribution.sample(nr_samples, "M")
+
 
                 if len(data.uncertain_parameters) > 1:
                     U_mc[feature] = U_hat[feature](*samples)
 
-                    data[feature].sobol_first = cp.Sens_m(U_hat[feature], distribution)
-                    data[feature].sobol_total = cp.Sens_t(U_hat[feature], distribution)
+                    #data[feature].sobol_first = cp.Sens_m(U_hat[feature], distribution)
+                    e = np.array(U_hat[feature].exponents,   dtype=int) 
+                    dim = e.shape[1]
+                    poly = P
+                    data[feature].sobol_first, data[feature].sobol_total = self.calc_sens(dim, poly, approx[feature])
+                    #data[feature].sobol_total = cp.Sens_t(U_hat[feature], distribution)
                     data = self.average_sensitivity(data, sensitivity="sobol_first")
                     data = self.average_sensitivity(data, sensitivity="sobol_total")
 
                 else:
                     U_mc[feature] = U_hat[feature](samples)
 
+                if save_samples:
+                    data[feature].samples = U_mc[feature]
                 data[feature].percentile_5 = np.percentile(U_mc[feature], 5, -1)
                 data[feature].percentile_95 = np.percentile(U_mc[feature], 95, -1)
 
@@ -1194,6 +1261,7 @@ def polynomial_chaos(self,
                          nr_pc_mc_samples=10**4,
                          allow_incomplete=True,
                          seed=None,
+                         save_samples=False,
                          **custom_kwargs):
         """
         Perform an uncertainty quantification and sensitivity analysis
@@ -1347,7 +1415,7 @@ def polynomial_chaos(self,
                                                            nr_collocation_nodes=nr_collocation_nodes,
                                                            allow_incomplete=allow_incomplete)
             else:
-                U_hat, distribution, data = \
+                U_hat, approx, distribution, data, P = \
                     self.create_PCE_collocation(uncertain_parameters=uncertain_parameters,
                                                 polynomial_order=polynomial_order,
                                                 nr_collocation_nodes=nr_collocation_nodes,
@@ -1380,18 +1448,20 @@ def polynomial_chaos(self,
         else:
             raise ValueError("No polynomial chaos method with name {}".format(method))
 
-        data = self.analyse_PCE(U_hat, distribution, data, nr_samples=nr_pc_mc_samples)
+        data = self.analyse_PCE(U_hat, approx, distribution, data, P, polynomial_order=polynomial_order,
+                                nr_samples=nr_pc_mc_samples,
+                                save_samples=save_samples)
 
         data.seed = seed
 
         return data
 
-
     def monte_carlo(self,
                     uncertain_parameters=None,
                     nr_samples=10**4,
                     seed=None,
-                    allow_incomplete=True):
+                    allow_incomplete=True,
+                    save_samples=False):
         """
         Perform an uncertainty quantification using the quasi-Monte Carlo method.
 
@@ -1411,6 +1481,8 @@ def monte_carlo(self,
             If the uncertainty quantification should be performed for features
             or models with incomplete evaluations.
             Default is True.
+        save_samples: bool, optional
+            Save samples in feature data to, for example, plot PDFs later.
 
         Returns
         -------
@@ -1471,34 +1543,7 @@ def monte_carlo(self,
         uncertainpy.Parameters
         """
 
-        if seed is not None:
-            np.random.seed(seed)
-
-        uncertain_parameters = self.convert_uncertain_parameters(uncertain_parameters)
-
-        distribution = self.create_distribution(uncertain_parameters=uncertain_parameters)
-
-        # nodes = distribution.sample(nr_samples, "M")
-
-        problem = {
-            "num_vars": len(uncertain_parameters),
-            "names": uncertain_parameters,
-            "bounds": [[0,1]]*len(uncertain_parameters)
-        }
-
-        # Create the Multivariate normal distribution
-        dist_R = []
-        for parameter in uncertain_parameters:
-            dist_R.append(cp.Uniform())
-
-        dist_R = cp.J(*dist_R)
-
-        nr_sobol_samples = int(np.round(nr_samples/2.))
-
-        nodes_R = saltelli.sample(problem, nr_sobol_samples, calc_second_order=False)
-
-        nodes = distribution.inv(dist_R.fwd(nodes_R.transpose()))
-
+        uncertain_parameters, distribution, _, nodes, nr_sobol_samples = self.get_UQ_samples("MC", uncertain_parameters, nr_samples)
 
         data = self.runmodel.run(nodes, uncertain_parameters)
 
@@ -1522,6 +1567,8 @@ def monte_carlo(self,
             logger = get_logger(self)
 
             if (np.all(mask) or allow_incomplete) and sum(mask) > 0:
+                if save_samples:
+                    data[feature].samples = masked_evaluations
                 data[feature].mean = np.mean(masked_evaluations, 0)
                 data[feature].variance = np.var(masked_evaluations, 0)
 
@@ -1694,4 +1741,4 @@ def average_sensitivity(self, data, sensitivity="sobol_first"):
                 data[feature][sensitivity + "_average"] = np.array(total_sense)
 
 
-        return data
+        return data

src/uncertainpy/data.py

@@ -2,7 +2,10 @@
 
 import six
 import os
-import collections
+try:
+    from collections import MutableMapping
+except ImportError:
+    from collections.abc import MutableMapping
 
 import numpy as np
 
@@ -11,7 +14,7 @@
 from ._version import __version__
 
 
-class DataFeature(collections.MutableMapping):
+class DataFeature(MutableMapping):
     """
     Store the results of each statistical metric calculated from the uncertainty
     quantification and sensitivity analysis for a single model/feature.
@@ -291,7 +294,7 @@ def ndim(self):
         return None
 
 
-class Data(collections.MutableMapping):
+class Data(MutableMapping):
     """
     Store the results of each statistical metric calculated from the uncertainty
     quantification and sensitivity analysis for each model/features.
@@ -696,7 +699,6 @@ def add_group(group, values, name="evaluation"):
 
         for feature in self.data:
             group = f.create_group(feature)
-
             for statistical_metric in self[feature]:
                 if statistical_metric in ["evaluations", "time"]:
                     if is_regular(self[feature][statistical_metric]):
@@ -810,6 +812,9 @@ def append_evaluations(evaluations, group):
         if "model ignore" in f.attrs:
             self.model_ignore = f.attrs["model ignore"]
 
+        if "samples" in f.attrs:
+            self.samples = f.attrs["samples"]
+
 
         for feature in f:
             self.add_features(str(feature))

src/uncertainpy/parameters.py

@@ -5,6 +5,10 @@
 import fileinput
 import sys
 import collections
+try:
+    from collections import MutableMapping
+except ImportError:
+    from collections.abc import MutableMapping
 
 import chaospy as cp
 
@@ -145,7 +149,7 @@ def __str__(self):
 
 
 
-class Parameters(collections.MutableMapping):
+class Parameters(MutableMapping):
     """
     A collection of parameters.