Source code for gemseo.uncertainty.distributions.openturns.distribution

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# 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
#    OTHER AUTHORS   - MACROSCOPIC CHANGES

"""Class to create a probability distribution from the OpenTURNS library.

The :class:`.OTDistribution` class is a concrete class
inheriting from :class:`.Distribution` which is an abstract one.
OT stands for `OpenTURNS <http://www.openturns.org/>`_
which is the library it relies on.

The :class:`.OTDistribution` of a given uncertain variable is built
from mandatory arguments:

- a variable name,
- a distribution name recognized by OpenTURNS,
- a set of parameters provided as a tuple
  of positional arguments filled in the order
  specified by the OpenTURNS constructor of this distribution.

.. warning::

    The distribution parameters must be provided according to the signature
    of the openTURNS classes. `Access the openTURNS documentation
    <http://openturns.github.io/openturns/latest/user_manual/
    probabilistic_modelling.html>`_.

The constructor has also optional arguments:

- a variable dimension (default: 1),
- a standard representation of these parameters
  (default: use the parameters provided in the tuple),
- a transformation of the variable (default: no transformation),
- lower and upper bounds for truncation (default: no truncation),
- a threshold for the OpenTURNS truncation tool
  (`more details <http://openturns.github.io/openturns/latest/user_manual/
  _generated/openturns.TruncatedDistribution.html>`_).
"""

from __future__ import division, unicode_literals

import logging
from typing import Callable, Iterable, List, Optional

import openturns as ots
from numpy import array, inf, ndarray

from gemseo.uncertainty.distributions.distribution import (
    Distribution,
    ParametersType,
    StandardParametersType,
)
from gemseo.uncertainty.distributions.openturns.composed import OTComposedDistribution
from gemseo.utils.string_tools import MultiLineString

OT_WEBSITE = "http://openturns.github.io/openturns/latest/"
OT_WEBSITE += "user_manual/probabilistic_modelling.html"

LOGGER = logging.getLogger(__name__)


[docs]class OTDistribution(Distribution): """OpenTURNS probability distribution. Create a probability distribution for an uncertain variable from its dimension and distribution names and properties. Example: >>> from gemseo.uncertainty.distributions.openturns.distribution import ( ... OTDistribution ... ) >>> distribution = OTDistribution('x', 'Exponential', (3, 2)) >>> print(distribution) Exponential(3, 2) """ _COMPOSED_DISTRIBUTION = OTComposedDistribution def __init__( self, variable, # type: str interfaced_distribution, # type: str parameters, # type: ParametersType dimension=1, # type: int standard_parameters=None, # type: Optional[StandardParametersType] transformation=None, # type: Optional[str] lower_bound=None, # type: Optional[float] upper_bound=None, # type: Optional[float] threshold=0.5, # type: float ): # noqa: D205,D212,D415 # type: (...) -> None """ Args: variable: The name of the random variable. distribution: The name of the probability distribution, typically the name of a class wrapped from an external library, such as 'Normal' for OpenTURNS or 'norm' for SciPy. parameters: The parameters of the probability distribution. dimension: The dimension of the random variable. standard_parameters: The standard representation of the parameters of the probability distribution. transformation: A transformation applied to the random variable, e.g. 'sin(x)'. If None, no transformation. lower_bound: A lower bound to truncate the distribution. If None, no lower truncation. upper_bound: An upper bound to truncate the distribution. If None, no upper truncation. threshold: A threshold in [0,1]. """ super(OTDistribution, self).__init__( variable, interfaced_distribution, parameters, dimension, standard_parameters, ) self.marginals = self.__create_distributions( self.distribution_name, self.parameters, transformation, lower_bound, upper_bound, threshold, ) self.distribution = ots.ComposedDistribution(self.marginals) msg = MultiLineString() msg.indent() msg.add("Mathematical support: {}", self.support) msg.add("Numerical range: {}", self.range) msg.add("Transformation: {}", self.transformation) LOGGER.info("%s", msg)
[docs] def compute_samples( self, n_samples=1, # type: int ): # noqa: D102 # type: (...) -> ndarray sample = array(self.distribution.getSample(n_samples)) return sample
[docs] def compute_cdf( self, vector, # type: Iterable[float] ): # noqa: D102 # type: (...) -> ndarray return array( [ self.marginals[index].computeCDF(ots.Point([value])) for index, value in enumerate(vector) ] )
[docs] def compute_inverse_cdf( self, vector, # type: Iterable[float] ): # noqa: D102 # type: (...) -> ndarray return array( [ self.marginals[index].computeQuantile(value)[0] for index, value in enumerate(vector) ] )
def _pdf( self, index, # type: int ): # noqa: D102 # type: (...) -> Callable def pdf( point, # type: float ): # type: (...) -> float """Probability Density Function (PDF). Args: point: An evaluation point. Returns: The PDF value at the evaluation point. """ return self.marginals[index].computePDF(point) return pdf def _cdf( self, index, # type: int ): # noqa: D102 # type: (...) -> Callable def cdf( level, # type: float ): # type: (...) -> float """Cumulative Density Function (CDF). Args: level: A probability level. Returns: The CDF value for the probability level. """ return self.marginals[index].computeCDF(level) return cdf @property def mean(self): # noqa: D102 # type: (...) -> ndarray return array(self.distribution.getMean()) @property def standard_deviation(self): # noqa: D102 # type: (...) -> ndarray return array(self.distribution.getStandardDeviation()) def __create_distributions( self, distribution, # type: str parameters, # type:ParametersType transformation, # type: str lower_bound, # type:float upper_bound, # type:float threshold, # type:float ): # type: (...) -> List[ots.Distribution] """Instantiate an OpenTURNS distribution for each random variable component. Args: distribution: The name of the distribution. parameters: The parameters of the distribution. transformation: A transformation applied to the random variable, e.g. 'sin(x)'. If None, no transformation. lower_bound: A lower truncation bound. If None, no lower truncation. upper_bound: An upper truncation bound. If None, no upper truncation. threshold: A threshold value in [0,1]. Returns: The marginal OpenTURNS distributions. """ try: ot_dist = getattr(ots, distribution) except Exception: raise ValueError( "{} is an unknown OpenTURNS distribution.".format(distribution) ) try: ot_dist = [ot_dist(*parameters)] * self.dimension except Exception: args = ", ".join([str(val) for val in parameters]) raise ValueError( "Arguments are wrong in {}({}). " "More details on: {}.".format(distribution, args, OT_WEBSITE) ) self.__set_bounds(ot_dist) if transformation is not None: ot_dist = self.__transform_marginal_dist(ot_dist, transformation) self.__set_bounds(ot_dist) if lower_bound is not None or upper_bound is not None: ot_dist = self.__truncate_marginal_dist( ot_dist, lower_bound, upper_bound, threshold ) self.__set_bounds(ot_dist) return ot_dist def __transform_marginal_dist( self, marginals, # type: Iterable[ots.Distribution], transformation, # type: str ): # type: (...) -> List[ots.Distribution] """Apply the standard transformations on the marginals. Examples of transformations: -, +, *, **, sin, exp, log, ... Args: marginals: The marginal distributions. transformation: A transformation applied to the random variable, e.g. 'sin(x)'. If None, no transformation. Returns: The transformed marginal distributions. """ variable_name = self.variable_name transformation = transformation.replace(" ", "") symbolic_function = ots.SymbolicFunction([self.variable_name], [transformation]) marginals = [ ots.CompositeDistribution(symbolic_function, marginal) for marginal in marginals ] prev = self.transformation transformation = transformation.replace(variable_name, "({})".format(prev)) self.transformation = transformation return marginals def __truncate_marginal_dist( self, distributions, # type: Iterable[ots.Distribution] lower_bound, # type: float upper_bound, # type: float threshold=0.5, # type: float ): # type: (...) -> List[ots.Distribution] """Truncate the distribution of a random variable. Args: distributions: The distributions. lower_bound: A lower bound to truncate the distributions. If None, no lower truncation. upper_bound: An upper bound to truncate the distributions. If None, no upper truncation. threshold: A threshold in [0,1]. Returns: The transformed distributions. """ marginals = [ self.__truncate_distribution( dist, index, lower_bound, upper_bound, threshold ) for index, dist in enumerate(distributions) ] prev = self.transformation self.transformation = "Trunc({})".format(prev) return marginals def __truncate_distribution( self, distributions, # type: Iterable[ots.Distribution] index, # type: int lower_bound, # type: float upper_bound, # type: float threshold=0.5, # type: float ): # type: (...) -> List[ots.Distribution] """Truncate a distribution with lower bound, upper bound or both. Args: distributions: The distributions. index: A random variable component. lower_bound: A lower bound to truncate the distribution. If None, no lower truncation. upper_bound: An upper bound to truncate the distribution. If None, no upper truncation. threshold: A threshold in [0,1]. Returns: The truncated distributions. """ if lower_bound is None: LOGGER.info( "Truncate distribution of component %s above %s.", index, upper_bound ) upper = ots.TruncatedDistribution.UPPER current_u_b = self.math_upper_bound[index] if upper_bound > current_u_b: raise ValueError("u_b is greater than the current upper bound.") distributions = ots.TruncatedDistribution( distributions, upper_bound, upper, threshold ) elif upper_bound is None: LOGGER.info( "Truncate distribution of component %s" " below %s.", index, lower_bound ) lower = ots.TruncatedDistribution.LOWER current_l_b = self.math_lower_bound[index] if lower_bound < current_l_b: raise ValueError("l_b is lower " "than the current lower bound.") distributions = ots.TruncatedDistribution( distributions, lower_bound, lower, threshold ) else: LOGGER.info( "Truncate distribution of component %s" " below %s and above %s.", index, lower_bound, upper_bound, ) current_l_b = self.math_lower_bound[index] current_u_b = self.math_upper_bound[index] if lower_bound < current_l_b: raise ValueError("l_b is lower " "than the current lower bound.") if upper_bound > current_u_b: raise ValueError("u_b is greater " "than the current upper bound.") distributions = ots.TruncatedDistribution( distributions, lower_bound, upper_bound, threshold ) return distributions def __set_bounds( self, distributions # type: Iterable[ots.Distribution] ): # type: (...) -> None """Set the mathematical and numerical bounds (= support and range). Args: distributions: The distributions. """ self.math_lower_bound = [] self.math_upper_bound = [] self.num_lower_bound = [] self.num_upper_bound = [] for distribution in distributions: dist_range = distribution.getRange() l_b = dist_range.getLowerBound()[0] u_b = dist_range.getUpperBound()[0] self.num_lower_bound.append(l_b) self.num_upper_bound.append(u_b) if not dist_range.getFiniteLowerBound()[0]: l_b = -inf if not dist_range.getFiniteUpperBound()[0]: u_b = inf self.math_lower_bound.append(l_b) self.math_upper_bound.append(u_b) self.math_lower_bound = array(self.math_lower_bound) self.math_upper_bound = array(self.math_upper_bound) self.num_lower_bound = array(self.num_lower_bound) self.num_upper_bound = array(self.num_upper_bound)