Source code for gemseo.mlearning.regression.pce

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry,
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

# Contributors:
#    INITIAL AUTHORS - initial API and implementation and/or initial
#                         documentation
#        :author: Matthias De Lozzo
Polynomial chaos expansion

The polynomial chaos expansion (PCE) discipline expresses the model output
as a weighted sum of polynomial functions which are orthonormal
in the stochastic input space spanned by the random input variables:

.. math::

    Y = w_0 + w_1\phi_1(X) + w_2\phi_2(X) + ... + w_K\phi_K(X)

where :math:`\phi_i(x)=\phi_{\tau_1(i),1}(x_1)\times\ldots\times

Enumerating strategy

The choice of the function :math:`\tau=(\tau_1,\ldots,\tau_d)` is an
enumerating strategy and :math:`\tau_j(i)` is the polynomial degree of


PCE are stochastic models whose inputs are random variables
and are often used to deal with uncertainty quantification problems.

If :math:`X_j` is a Gaussian random variable,
:math:`(\psi_{ij})_{i\geq 0}` is the Legendre basis.
If :math:`X_j` is an uniform random variable,
:math:`(\psi_{ij})_{i\geq 0}` is the Hermite basis.

When the problem is deterministic, we can still use PCE
under the assumptions that
the random variables are independent uniform random variables.
Then, the orthonormal basis function is the Hermite basis.


The degree :math:`P` of a PCE is defined
in such a way that :math:`\text{degree}(\psi_i)=\sum_{j=1}^d\tau_j(i)\leq P`.


The coefficients :math:`(w_1, w_2, ..., w_K)` and the intercept
:math:`w_0` are estimated either by least square regression,
sparse least square regression or quadrature.

The PCE model relies on the FunctionalChaosAlgorithm class
of the `openturns library <
from __future__ import absolute_import, division, unicode_literals

from future import standard_library
from numpy import all as np_all
from numpy import array, concatenate, isin, zeros, zeros_like
from openturns import (

from gemseo.mlearning.regression.regression import MLRegressionAlgo


from gemseo import LOGGER

[docs]class PCERegression(MLRegressionAlgo): """ Polynomial chaos expansion. """ LIBRARY = "openturns" ABBR = "PCE" LS_STRATEGY = "LS" QUAD_STRATEGY = "Quad" SPARSE_STRATEGY = "SparseLS" AVAILABLE_STRATEGIES = [LS_STRATEGY, QUAD_STRATEGY, SPARSE_STRATEGY] def __init__( self, data, probability_space, discipline=None, transformer=None, input_names=None, output_names=None, strategy=LS_STRATEGY, degree=2, n_quad=None, stieltjes=True, sparse_param=None, ): """Constructor. :param data: learning dataset :type data: Dataset :param probability_space: probability space. :type probability_space: ParameterSpace :param discipline: discipline to evaluate if strategy='Quad' and data is empty. :type discipline: MDODiscipline :param transformer: transformation strategy for data groups. If None, do not transform data. Default: None. :type transformer: dict(str) :param input_names: names of the input variables. :type input_names: list(str) :param output_names: names of the output variables. :type output_names: list(str) :param strategy: strategy to compute the parameters of the PCE, either 'LS', 'Quad' or 'SparseLS'. Default: 'LS'. :type strategy: str :param degree: polynomial degree of the PCE :type degree: int :param n_quad: number of quadrature points :type n_quad: int :param stieltjes: stieltjes :type stieltjes: bool :param sparse_param: Parameters for the Sparse Cleaning Truncation Strategy and/or hyperbolic truncation of the initial basis: - **max_considered_terms** (int) -- Maximum Considered Terms, - **most_significant** (int), Most Significant number to retain, - **significance_factor** (float), Significance Factor, - **hyper_factor** (float), factor for hyperbolic truncation strategy. """ prob_space = probability_space super(PCERegression, self).__init__( data, transformer=transformer, input_names=input_names, output_names=output_names, probability_space=prob_space, strategy=strategy, degree=degree, n_quad=n_quad, stieltjes=stieltjes, sparse_param=sparse_param, ) self.prob_space = prob_space self.discipline = discipline if transformer is not None and transformer != {}: raise ValueError("PCERegression does not support transformers.") try: if data: u_names = set(prob_space.variables_names) assert set(self.input_names) == u_names else: self.input_names = prob_space.variables_names except Exception: raise ValueError( "Data inputs names are %s " "while probability distributions are defined " "%s." % (self.input_names, list(prob_space.marginals.keys())) ) self.distributions = prob_space.marginals self.sparse_param = sparse_param or {} self.input_dim = sum([dist.dimension for _, dist in self.distributions.items()]) self.strategy = strategy if self.strategy not in self.AVAILABLE_STRATEGIES: strategies = " ".join(self.AVAILABLE_STRATEGIES) raise ValueError( "{} is a wrong strategy." " Available ones are: {}".format(self.strategy, strategies) ) self.stieltjes = stieltjes = degree self.n_quad = n_quad self.ot_distributions = [ self.distributions[name].distribution for name in self.input_names ] self.dist = ComposedDistribution(self.ot_distributions) hyper_factor = self.sparse_param.get("hyper_factor", 1.0) self.enumerate_function = EnumerateFunction(self.input_dim, hyper_factor) self.basis = self._get_basis() self.n_basis = self._get_basis_size() self.trunc_strategy = self._get_trunc_strategy() if self.strategy == self.QUAD_STRATEGY: ( self.sample, self.weights, self.proj_strategy, ) = self._get_quadrature_points() else: self.sample = None self.weights = None def _fit(self, input_data, output_data): """Fit the regression model. :param ndarray input_data: input data (2D) :param ndarray output_data: output data (2D) """ self.proj_strategy = self._get_proj_strategy(input_data, output_data) weights = self._get_weights(input_data) self.algo = self._build_pce(input_data, weights, output_data) def _predict(self, input_data): """Predict output. :param ndarray input_data: input data (2D). :return: output prediction (2D). :rtype: ndarray """ return array(self.algo.getMetaModel()(input_data))
[docs] def compute_sobol(self): """Compute first and total Sobol' indices :returns: first and total SObol' indices :rtype: list, list """ sensitivity_analysis = FunctionalChaosSobolIndices(self.algo) dimension = len(self.input_names) first_order = [sensitivity_analysis.getSobolIndex(i) for i in range(dimension)] total_order = [ sensitivity_analysis.getSobolTotalIndex(i) for i in range(dimension) ] return first_order, total_order
def _build_pce(self, x_learn, weights, y_learn): """Build PCE""" pce_algo = FunctionalChaosAlgorithm( x_learn, weights, y_learn, self.dist, self.trunc_strategy ) pce_algo.setProjectionStrategy(self.proj_strategy) return pce_algo.getResult() def _get_basis(self): """Get basis function for PCE construction""" if self.stieltjes: # Tend to result in performance issue basis = OrthogonalProductPolynomialFactory( [ StandardDistributionPolynomialFactory( AdaptiveStieltjesAlgorithm(marginal) ) for marginal in self.ot_distributions ], self.enumerate_function, ) else: basis = OrthogonalProductPolynomialFactory( [ StandardDistributionPolynomialFactory(margin) for margin in self.ot_distributions ], self.enumerate_function, ) return basis def _get_basis_size(self): """Get basis size for PCE construction""" return self.enumerate_function.getStrataCumulatedCardinal( def _get_quadrature_points(self): """Get quadrature points for PCE construction""" measure = self.basis.getMeasure() if self.n_quad is not None: degree_by_dim = int(self.n_quad ** (1.0 / self.input_dim)) else: degree_by_dim = + 1 degrees = [degree_by_dim] * self.input_dim proj_strategy = IntegrationStrategy(GaussProductExperiment(measure, degrees)) sample, weights = proj_strategy.getExperiment().generateWithWeights() if not self.stieltjes: transformation = Function( MarginalTransformationEvaluation( [measure.getMarginal(i) for i in range(self.input_dim)], self.ot_distributions, False, ) ) sample = transformation(sample) sample = array(sample) if not self.learning_set: inputs_names = self.prob_space.variables_names in_grp = self.learning_set.INPUT_GROUP self.learning_set.set_from_array( sample, self.prob_space.variables_names, self.prob_space.variables_sizes, {name: in_grp for name in inputs_names}, ) if self.discipline is not None: n_samples = len(self.learning_set) outputs = {name: [] for name in self.discipline.get_output_data_names()} for data in self.learning_set: output_data = self.discipline.execute(data) for name in self.discipline.get_output_data_names(): outputs[name] += list(output_data[name]) for name in self.discipline.get_output_data_names(): outputs[name] = array(outputs[name]).reshape((n_samples, -1)) outputs_names = list(self.discipline.get_output_data_names()) outputs = [outputs[name] for name in outputs_names] data = concatenate(outputs, axis=1) out_grp = self.learning_set.OUTPUT_GROUP sizes = { name: len(self.discipline.local_data[name]) for name in outputs_names } self.learning_set.add_group( out_grp, data, outputs_names, sizes, cache_as_input=False ) self.output_names = outputs_names return sample, weights, proj_strategy def _get_trunc_strategy(self): """Get truncation strategy for PCE construction""" if self.strategy in [self.LS_STRATEGY, self.QUAD_STRATEGY]: trunc_strategy = FixedStrategy(self.basis, self.n_basis) elif self.strategy == self.SPARSE_STRATEGY: sparse_param = self.sparse_param sparse_param = {} if sparse_param is None else sparse_param max_considered_terms = sparse_param.get("max_considered_terms", 120) most_significant = sparse_param.get("most_significant", 30) significance_factor = sparse_param.get("significance_factor", 1e-3) trunc_strategy = CleaningStrategy( OrthogonalBasis(self.basis), max_considered_terms, most_significant, significance_factor, True, ) return trunc_strategy def _get_proj_strategy(self, x_learn, y_learn): """Get projection strategy for PCE construction :param x_learn: input data :type x_learn: array :param y_learn: output data :type y_learn: array """ if self.strategy == self.QUAD_STRATEGY: proj_strategy = self.proj_strategy elif self.strategy == self.LS_STRATEGY: proj_strategy = LeastSquaresStrategy(x_learn, y_learn) else: app = LeastSquaresMetaModelSelectionFactory(LARS(), CorrectedLeaveOneOut()) proj_strategy = LeastSquaresStrategy(x_learn, y_learn, app) return proj_strategy def _get_ls_weights(self): """Get LS weights for PCE construction""" _, weights = self.proj_strategy.getExperiment().generateWithWeights() return weights def _get_quad_weights(self, x_learn): """Get quadrature weights for PCE construction""" sample = zeros_like(self.sample) common_len = len(x_learn) sample[:common_len] = x_learn sample[0] = self.sample[0] sample_arg = np_all(isin(sample, self.sample), axis=1) new_weights = array(self.weights)[sample_arg] return new_weights def _get_weights(self, x_learn): """Get weights for PCE construction)""" if self.strategy == self.QUAD_STRATEGY: weights = self._get_quad_weights(x_learn) else: weights = self._get_ls_weights() return weights def _predict_jacobian(self, input_data): """Predict Jacobian of the regression model for the given input data. :param ndarray input_data: input_data (2D). :return: Jacobian matrices (3D, one for each sample). :rtype: ndarray """ input_shape, output_shape = self._get_raw_shapes() gradient = self.algo.getMetaModel().gradient jac = zeros((input_data.shape[0], int(output_shape), int(input_shape))) for index, data in enumerate(input_data): jac[index] = array(gradient(Point(data))).T return jac