From ce6d551b895807c931340f35582200dbd5ca9372 Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Mon, 25 Nov 2024 14:54:57 +0100 Subject: [PATCH] Drop ConditionalDistribution --- doc/examples/plot_conditionaldistribution.py | 82 ----- doc/user_manual/user_manual.rst | 1 - otbenchmark/ConditionalDistribution.rst | 115 ------- otbenchmark/_ConditionalDistribution.py | 342 ------------------- otbenchmark/_CrossCutDistribution.py | 6 +- otbenchmark/__init__.py | 2 - requirements.txt | 4 +- tests/test_ConditionalDistribution.py | 165 --------- tests/test_CrossCutDistribution.py | 2 +- tests/test_CrossCutFunction.py | 2 +- 10 files changed, 7 insertions(+), 714 deletions(-) delete mode 100644 doc/examples/plot_conditionaldistribution.py delete mode 100644 otbenchmark/ConditionalDistribution.rst delete mode 100644 otbenchmark/_ConditionalDistribution.py delete mode 100644 tests/test_ConditionalDistribution.py diff --git a/doc/examples/plot_conditionaldistribution.py b/doc/examples/plot_conditionaldistribution.py deleted file mode 100644 index c111b2f..0000000 --- a/doc/examples/plot_conditionaldistribution.py +++ /dev/null @@ -1,82 +0,0 @@ -""" -Conditional distributions -========================= -""" - -# %% -# Conditioning is a way to reduce the dimensionnality of a multivariate distribution. - -# %% -import openturns as ot -import otbenchmark as otb -import openturns.viewer as otv - -# %% -# Conditional distribution of a three dimensional gaussian distribution -# --------------------------------------------------------------------- - -# %% -# The random variable is (X0, X1, X2). -distribution = ot.Normal(3) - -# %% -# We condition with respect to X1=mu1, i.e. we consider (X0, X1, X2) | X1=2. -conditionalIndices = [1] -conditionalReferencePoint = [2.0] -conditionalDistribution = ot.Distribution( - otb.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) -) - - -# %% -_ = otv.View(conditionalDistribution.drawPDF()) - -# %% -# Conditional distribution of a three dimensional mixture -# ------------------------------------------------------- - -# %% -# Create a Funky distribution -corr = ot.CorrelationMatrix(3) -corr[0, 1] = 0.2 -copula = ot.NormalCopula(corr) -x1 = ot.Normal(-1.0, 1.0) -x2 = ot.Normal(2.0, 1.0) -x3 = ot.Normal(1.0, 1.0) -x_funk = ot.ComposedDistribution([x1, x2, x3], copula) - - -# %% -# Create a Punk distribution -x1 = ot.Normal(1.0, 1.0) -x2 = ot.Normal(-2, 1.0) -x3 = ot.Normal(2.0, 1.0) -x_punk = ot.ComposedDistribution([x1, x2, x3], copula) - - -# %% -distribution = ot.Mixture([x_funk, x_punk], [0.5, 1.0]) - - -# %% -referencePoint = distribution.getMean() -referencePoint - - -# %% -conditionalIndices = [1] -conditionalReferencePoint = [-0.5] -conditionalDistribution = ot.Distribution( - otb.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) -) - - -# %% -_ = otv.View(conditionalDistribution.drawPDF()) - -# %% -otv.View.ShowAll() diff --git a/doc/user_manual/user_manual.rst b/doc/user_manual/user_manual.rst index 842653c..b695c80 100644 --- a/doc/user_manual/user_manual.rst +++ b/doc/user_manual/user_manual.rst @@ -52,7 +52,6 @@ Reliability methods LHS ReliabilityBenchmarkMetaAlgorithm ReliabilityBenchmarkResult - ConditionalDistribution CrossCutFunction CrossCutDistribution DrawEvent diff --git a/otbenchmark/ConditionalDistribution.rst b/otbenchmark/ConditionalDistribution.rst deleted file mode 100644 index 2f92088..0000000 --- a/otbenchmark/ConditionalDistribution.rst +++ /dev/null @@ -1,115 +0,0 @@ -Marginal and conditional distributions -====================================== - -Conditional distribution ------------------------- - -For a given pair of continuous random variables :math:`X` and :math:`Y` -with probability densities :math:`f_X` and :math:`f_Y`, the conditional -density of :math:`Y|X` is: - -.. math:: - - - f_{Y|X}(y|X = x) = \frac{f_{X,Y}(x,y)}{f_X(x)}, - -where :math:`f_X(x)` is the marginal distribution of :math:`X`: - -.. math:: - - - f_X(x) = \int_{\mathbb{R}^{|Y|}} f_{X,Y}(x,y) dy. - -In the previous equation, the integer :math:`|Y|` is the dimension of -the vector :math:`Y`. - -For a multivariate random variable :math:`X=(X_1, X_2, ..., X_d)` and we -consider conditional distribution which arise when we condition with -respect with several marginals of :math:`X`. More precisely, we assume -that a given set of conditioning indices - -.. math:: - - - C=(i_1, i_2, ..., i_c) - -are given, where :math:`0\leq c\leq d` is an integeger, -:math:`i_1, i_2, ..., i_c \in \{1, 2, ..., d\}` and are so that -:math:`i_1 < i_2 < ... >> import openturns as ot - >>> import otbenchmark as otb - >>> distribution = ot.Normal(3) - # We condition X = (X0, X1, X2) given X1=2.0, X2=3.0 - >>> conditionalIndices = [1, 2] - >>> conditionalRefencePoint = [2.0, 3.0] - >>> conditionalDistribution = ot.Distribution(otb.ConditionalDistribution( - ... distribution, conditionalIndices, conditionalRefencePoint)) - # PDF - >>> computed = conditionalDistribution.computePDF([1.0]) - """ - if len(conditionalIndices) != len(conditionalRefencePoint): - raise ValueError( - "The dimension of the indices is %d" - "but the dimension of the reference point is %s" - % (len(conditionalIndices), len(conditionalRefencePoint)) - ) - self.conditionalIndices = conditionalIndices - self.conditionalRefencePoint = conditionalRefencePoint - self.distribution = distribution - self.extendedDimension = distribution.getDimension() - self.conditionalDimension = distribution.getDimension() - len( - conditionalIndices - ) - self.distributionRange = ConditionalDistribution.ComputeRange( - distribution, conditionalIndices, conditionalRefencePoint - ) - self.verbose = verbose - self.useIntegration = useIntegration - self.sampleSizeForCDFBySampling = sampleSizeForCDFBySampling - # Compute integration factor - rule = ot.GaussKronrodRule(ot.GaussKronrodRule.G11K23) - univariateIntegration = ot.GaussKronrod(maximumSubIntervals, maximumError, rule) - self.quadrature = ot.IteratedQuadrature(univariateIntegration) - if len(self.conditionalIndices) == 0: - self.integrationFactor = 1.0 - else: - conditionedDistribution = self.distribution.getMarginal( - self.conditionalIndices - ) - self.integrationFactor = conditionedDistribution.computePDF( - self.conditionalRefencePoint - ) - # Compute unconditioned indices - dimension = self.distribution.getDimension() - self.unconditionedIndices = [] - for i in range(dimension): - if i not in self.conditionalIndices: - self.unconditionedIndices.append(i) - # Create the mask for extending inputs - self.mask = np.zeros(self.extendedDimension, dtype=bool) - self.mask[self.conditionalIndices] = True - super(ConditionalDistribution, self).__init__(self.conditionalDimension) - - def getRange(self): - """ - Returns the range. - - Returns - ------- - conditionalRange : ot.Interval - The range of the conditional distribution. - """ - return self.distributionRange - - def computeCDF(self, Y): - """ - Returns the CDF. - - Parameters - ---------- - Y : ot.Point(dimension_Y) - The input point. - - Returns - ------- - p : float - The CDF of the conditional distribution. - """ - if self.verbose: - print("computeCDF at", Y) - - if self.useIntegration: - p = self._compute_CDF_by_integration(Y) - else: - p = self._compute_CDF_by_sampling(Y) - return p - - def _compute_CDF_by_integration(self, Y): - """ - Compute the value of the CDF at point Y. - - Parameters - ---------- - Y : ot.Point(dimension_Y) - The point where the CDF must be evaluated. - - Returns - ------- - p : float, in [0, 1] - The value of the CDF. - - """ - - def DistributionPDFIntegrand(Y): - X = self.computeExtendedInput(Y) - pdf = self.distribution.computePDF(X) - return [pdf] - - distributionPDFIntegrandPy = ot.PythonFunction( - self.conditionalDimension, 1, DistributionPDFIntegrand - ) - # Bounds = (lowerBound, X) - lowerBound = self.distributionRange.getLowerBound() - integrationRange = ot.Interval(lowerBound, Y) - value = self.quadrature.integrate(distributionPDFIntegrandPy, integrationRange) - p = value[0] / self.integrationFactor - p = min(1.0, max(0.0, p)) - return p - - def _compute_CDF_by_sampling(self, Y): - """ - Compute the value of the CDF at point Y by sampling. - - Parameters - ---------- - Y : ot.Point(dimension_Y) - The point where to evaluate the CDF. - - Returns - ------- - p : float, in [0, 1] - The value of the CDF. - - """ - - # Compute F_Y(X) where Y is the conditioned random vector. - unconditionedDistribution = self.distribution.getMarginal( - self.unconditionedIndices - ) - - # Create the domain from [-INF] up to Y. - unconditioned_dimension = len(Y) - lower_bound = [-ot.SpecFunc.MaxScalar] * unconditioned_dimension - upper_bound = Y - domain = ot.Interval(lower_bound, upper_bound) - # Create a domain event - unconditioned_random_vector = ot.RandomVector(unconditionedDistribution) - event = ot.DomainEvent(unconditioned_random_vector, domain) - event_sample = event.getSample(self.sampleSizeForCDFBySampling) - p = event_sample.computeMean()[0] - return p - - def computePDF(self, Y): - """ - Returns the PDF. - - Returns - ------- - y : float - The PDF of the conditional distribution. - """ - if self.verbose: - print("computePDF at", Y) - X = self.computeExtendedInput(Y) - pdfConditioned = self.distribution.computePDF(X) - pdf = pdfConditioned / self.integrationFactor - return pdf - - def getDescription(self): - """ - Returns the description. - - Returns - ------- - description : ot.Description - The description of the conditional distribution. - """ - # Does not work, because of https://github.com/openturns/openturns/issues/1230 - description = [] - extendedDescription = self.distribution.getDescription() - conditionalIndex = 0 - for i in range(self.distribution.getDimension()): - if i not in self.conditionalIndices: - description.append(extendedDescription[conditionalIndex]) - conditionalIndex += 1 - return ot.Description(description) - - def computeExtendedInput(self, Y): - """ - Compute the extended input from the conditioned input. - - The conditioned input is Y. - The function returns the vector - - X = (X[0], X[1], ..., X[dimension - 1]) - - so that - - X[i] = X[i] if i is conditionned - conditionalRefencePoint[i] otherwise, - - for i = 0, 1, ..., dimension - 1. - - The following algorithm implements this: - - .. code-block:: python - - X = [] - referenceIndex = 0 - conditionalIndex = 0 - for i in range(self.extendedDimension): - if i in self.conditionalIndices: - X.append(self.conditionalRefencePoint[referenceIndex]) - referenceIndex += 1 - else: - X.append(X[conditionalIndex]) - conditionalIndex += 1 - - The actual code is, however, vectorized using numpy. - - Parameters - ---------- - Y : ot.Point - The conditioned input point. - extendedDimension : int - The extended dimension. - conditionalIndices : list - The indices of the marginal which are conditioned - conditionalRefencePoint : list - The values of the marginal which are conditioned - - Returns - ------- - X : ot.Point - The extended input point, as an input of the full distribution. - """ - if self.verbose: - Y = ot.Point(Y) - print("computeExtendedInput, Y =", Y) - X = np.zeros(self.extendedDimension) - X[self.conditionalIndices] = self.conditionalRefencePoint - X[~self.mask] = Y - X = ot.Point(X) - if self.verbose: - print("X=", X) - return X diff --git a/otbenchmark/_CrossCutDistribution.py b/otbenchmark/_CrossCutDistribution.py index d141df7..7ffbdd1 100644 --- a/otbenchmark/_CrossCutDistribution.py +++ b/otbenchmark/_CrossCutDistribution.py @@ -5,9 +5,9 @@ """ import openturns as ot +import openturns.experimental as otexp import pylab as pl import openturns.viewer as otv -import otbenchmark as otb class CrossCutDistribution: @@ -60,7 +60,7 @@ def drawConditionalPDF(self, referencePoint): crossCutIndices.append(k) crossCutReferencePoint.append(referencePoint[k]) conditionalDistribution = ot.Distribution( - otb.ConditionalDistribution( + otexp.PointConditionalDistribution( self.distribution, crossCutIndices, crossCutReferencePoint ) ) @@ -84,7 +84,7 @@ def drawConditionalPDF(self, referencePoint): crossCutIndices.append(k) crossCutReferencePoint.append(referencePoint[k]) conditionalDistribution = ot.Distribution( - otb.ConditionalDistribution( + otexp.PointConditionalDistribution( self.distribution, crossCutIndices, crossCutReferencePoint ) ) diff --git a/otbenchmark/__init__.py b/otbenchmark/__init__.py index fd0f2dd..3e23fb2 100644 --- a/otbenchmark/__init__.py +++ b/otbenchmark/__init__.py @@ -45,7 +45,6 @@ from ._GaussianSumSensitivity import GaussianSumSensitivity from ._GaussianProductSensitivity import GaussianProductSensitivity from ._GSobolSensitivity import GSobolSensitivity -from ._ConditionalDistribution import ConditionalDistribution from ._CrossCutFunction import CrossCutFunction from ._CrossCutDistribution import CrossCutDistribution from ._MorrisSensitivity import MorrisSensitivity @@ -104,7 +103,6 @@ "GaussianSumSensitivity", "GaussianProductSensitivity", "GSobolSensitivity", - "ConditionalDistribution", "CrossCutFunction", "CrossCutDistribution", "ProbabilitySimulationAlgorithmFactory", diff --git a/requirements.txt b/requirements.txt index 758555e..bbe503f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ -numpy<2 +numpy matplotlib -openturns>=1.23 +openturns>=1.24 black pre-commit flake8 diff --git a/tests/test_ConditionalDistribution.py b/tests/test_ConditionalDistribution.py deleted file mode 100644 index 01b567c..0000000 --- a/tests/test_ConditionalDistribution.py +++ /dev/null @@ -1,165 +0,0 @@ -# Copyright 2020 EDF. -""" -Test for ConditionalDistribution class. -""" -import openturns as ot -import otbenchmark -import unittest -import numpy as np -import openturns.viewer as otv - - -class CheckConditionalDistribution(unittest.TestCase): - def test_ComputeExtendedInput(self): - # The random variable is (X0, X1, X2) - distribution = ot.Normal(3) - # We condition with respect to X1=mu1, i.e. - # we consider (X0, X1, X2) | X1=2 - X = [1.1, 3.3] - conditionalIndices = [1] - conditionalReferencePoint = [2.2] - conditionalDistribution = otbenchmark.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) - extendedX = conditionalDistribution.computeExtendedInput(X) - np.testing.assert_equal(extendedX, ot.Point([1.1, 2.2, 3.3])) - - def test_ConditionalDistribution1(self): - # The random variable is (X0, X1, X2) - distribution = ot.Normal(3) - # We condition with respect to X1=mu1, i.e. - # we consider (X0, X1, X2) | X1=2 - conditionalIndices = [1] - conditionalReferencePoint = [2.0] - conditionalDistribution = ot.Distribution( - otbenchmark.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) - ) - # PDF - computed = conditionalDistribution.computePDF([1.0, 1.0]) - print("computed PDF=", computed) - conditionedPDF = distribution.computePDF( - [1.0, 1.0, conditionalReferencePoint[0]] - ) - n1 = ot.Normal(1) - conditioningPDF = n1.computePDF([conditionalReferencePoint[0]]) - exact = conditionedPDF / conditioningPDF - np.testing.assert_allclose(computed, exact) - # CDF - computed = conditionalDistribution.computeCDF([1.0, 1.0]) - print("computed CDF=", computed) - exact = 0.7078609817371252 - np.testing.assert_allclose(computed, exact) - # CDF when X0 = INF, X2 = INF - computed = conditionalDistribution.computeCDF([7.0, 7.0]) - print("computed CDF=", computed) - exact = 1.0 - np.testing.assert_allclose(computed, exact) - - def test_DrawPDF(self): - # The random variable is (X0, X1, X2) - distribution = ot.Normal(3) - # We condition with respect to X1=mu1, i.e. - # we consider (X0, X1, X2) | X1=2 - conditionalIndices = [1] - conditionalReferencePoint = [2.0] - conditionalDistribution = ot.Distribution( - otbenchmark.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) - ) - # Avoid failing on CircleCi - # _tkinter.TclError: no display name and no $DISPLAY environment variable - try: - graph = conditionalDistribution.drawPDF() - _ = otv.View(graph) - except Exception as e: - print(e) - - def test_ConditionalDistribution2(self): - # The random variable is (X0, X1, X2) - distribution = ot.Normal(3) - # We condition X = (X0, X1, X2) given X1=2.0, X2=3.0 - conditionalIndices = [1, 2] - conditionalReferencePoint = [2.0, 3.0] - for useIntegration in [True, False]: - conditionalDistribution = ot.Distribution( - otbenchmark.ConditionalDistribution( - distribution, - conditionalIndices, - conditionalReferencePoint, - useIntegration=useIntegration, - ) - ) - # Configure tolerance depending on algorithm - if useIntegration: - rtol = 1.0e-7 - else: - rtol = 1.0e-2 - # PDF - computed = conditionalDistribution.computePDF([1.0]) - print("computed PDF=", computed) - conditionedPDF = distribution.computePDF( - [1.0, conditionalReferencePoint[0], conditionalReferencePoint[1]] - ) - n2 = ot.Normal(2) - conditioningPDF = n2.computePDF( - [conditionalReferencePoint[0], conditionalReferencePoint[1]] - ) - exact = conditionedPDF / conditioningPDF - np.testing.assert_allclose(computed, exact) - # CDF - computed = conditionalDistribution.computeCDF([1.0]) - print("computed CDF=", computed) - exact = 0.8413447460685339 - np.testing.assert_allclose(computed, exact, rtol=rtol) - # CDF when X0 = INF, X2 = INF - computed = conditionalDistribution.computeCDF([7.0]) - print("computed CDF=", computed) - exact = 1.0 - np.testing.assert_allclose(computed, exact) - - def test_ConditionalDistribution3(self): - # Do not condition anything. - distribution = ot.Normal(3) - conditionalIndices = [] - conditionalReferencePoint = [] - conditionalDistribution = ot.Distribution( - otbenchmark.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) - ) - # PDF - computed = conditionalDistribution.computePDF([1.0, 1.0, 1.0]) - print("computed PDF=", computed) - exact = distribution.computePDF([1.0, 1.0, 1.0]) - np.testing.assert_allclose(computed, exact) - # CDF - computed = conditionalDistribution.computeCDF([1.0, 1.0, 1.0]) - print("computed CDF=", computed) - exact = distribution.computeCDF([1.0, 1.0, 1.0]) - np.testing.assert_allclose(computed, exact) - # CDF when X0 = INF, X2 = INF - computed = conditionalDistribution.computeCDF([7.0, 7.0, 7.0]) - print("computed CDF=", computed) - exact = 1.0 - np.testing.assert_allclose(computed, exact) - - def test_ConditionalDistribution4(self): - # Create a dimension 5 distribution - distribution = ot.Normal(5) - conditionalIndices = [1, 2] - conditionalReferencePoint = [2.0, 3.0] - conditionalDistribution = ot.Distribution( - otbenchmark.ConditionalDistribution( - distribution, conditionalIndices, conditionalReferencePoint - ) - ) - # PDF - computed = conditionalDistribution.computePDF([1.0] * 3) - print("computed PDF=", computed) - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_CrossCutDistribution.py b/tests/test_CrossCutDistribution.py index b2bef89..d0edb4e 100644 --- a/tests/test_CrossCutDistribution.py +++ b/tests/test_CrossCutDistribution.py @@ -1,6 +1,6 @@ # Copyright 2020 EDF. """ -Test for ConditionalDistribution class. +Test for CrossCutDistribution class. """ import otbenchmark import unittest diff --git a/tests/test_CrossCutFunction.py b/tests/test_CrossCutFunction.py index 8fa93c2..0e3604f 100644 --- a/tests/test_CrossCutFunction.py +++ b/tests/test_CrossCutFunction.py @@ -1,6 +1,6 @@ # Copyright 2020 EDF. """ -Test for ConditionalDistribution class. +Test for CrossCutDistribution class. """ import otbenchmark import unittest