From 745f6abc14266c9bd5efa79c3e90116606f2fcd3 Mon Sep 17 00:00:00 2001 From: regislebrun Date: Wed, 3 Apr 2024 22:31:19 +0200 Subject: [PATCH] Improved MarginalDistribution Now it is possible to switch between a CDF-based approach to compute the marginal PDF or a PDF-based approach (by default). It results into a huge speed-up of the PDF computations, see https://openturns.discourse.group/t/skew-distributions-in-ot/348/2. It fixes #2067. --- lib/etc/openturns.conf.in | 6 ++ lib/src/Base/Common/ResourceMap.cxx | 6 ++ .../Distribution/MarginalDistribution.cxx | 83 ++++++++++++++++++- .../openturns/MarginalDistribution.hxx | 37 ++++++++- 4 files changed, 127 insertions(+), 5 deletions(-) diff --git a/lib/etc/openturns.conf.in b/lib/etc/openturns.conf.in index db4c2d6aef2..d38271d5570 100644 --- a/lib/etc/openturns.conf.in +++ b/lib/etc/openturns.conf.in @@ -542,6 +542,12 @@ + + + + + + diff --git a/lib/src/Base/Common/ResourceMap.cxx b/lib/src/Base/Common/ResourceMap.cxx index e3755efdb47..b492fb50378 100644 --- a/lib/src/Base/Common/ResourceMap.cxx +++ b/lib/src/Base/Common/ResourceMap.cxx @@ -1169,6 +1169,12 @@ void ResourceMap::loadDefaultConfiguration() addAsUnsignedInteger("LogNormalFactory-EstimationMethod", 0); addAsUnsignedInteger("LogNormalFactory-MaximumIteration", 50); + // MarginalDistribution parameters // + addAsBool("MarginalDistribution-UsePDF", true); + addAsScalar("MarginalDistribution-MaximumError", 1.0e-7); + addAsString("MarginalDistribution-Rule", "G15K31"); + addAsUnsignedInteger("MarginalDistribution-MaximumSubIntervals", 128); + // Meixner parameters // addAsScalar("MeixnerDistribution-MaximumAbsoluteError", 1.0e-12); addAsScalar("MeixnerDistribution-MaximumConstraintError", 1.0e-12); diff --git a/lib/src/Uncertainty/Distribution/MarginalDistribution.cxx b/lib/src/Uncertainty/Distribution/MarginalDistribution.cxx index 58bb6729c17..af7e4172f04 100644 --- a/lib/src/Uncertainty/Distribution/MarginalDistribution.cxx +++ b/lib/src/Uncertainty/Distribution/MarginalDistribution.cxx @@ -22,6 +22,7 @@ #include "openturns/MarginalDistribution.hxx" #include "openturns/Uniform.hxx" +#include "openturns/ParametricFunction.hxx" #include "openturns/PersistentObjectFactory.hxx" BEGIN_NAMESPACE_OPENTURNS @@ -141,10 +142,10 @@ void MarginalDistribution::setDistributionAndIndices(const Distribution & distri setDimension(dimension); // Compute the range // From the underlying distribution - Point distributionLowerBound(distribution.getRange().getLowerBound()); + const Point distributionLowerBound(distribution.getRange().getLowerBound()); lowerBound_ = distributionLowerBound; Interval::BoolCollection distributionFiniteLowerBound(distribution.getRange().getFiniteLowerBound()); - Point distributionUpperBound(distribution.getRange().getUpperBound()); + const Point distributionUpperBound(distribution.getRange().getUpperBound()); upperBound_ = distributionUpperBound; Interval::BoolCollection distributionFiniteUpperBound(distribution.getRange().getFiniteUpperBound()); // For the marginal distribution @@ -192,6 +193,36 @@ Sample MarginalDistribution::getSample(const UnsignedInteger size) const return distribution_.getSample(size).getMarginal(indices_); } +/* Get the PDF of the MarginalDistribution */ +Scalar MarginalDistribution::computePDF(const Point & point) const +{ + if (point.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: expected a point of dimension=" << getDimension() << ", got dimension=" << point.getDimension(); + // First, check if the marginal distribution is just a eordering of the underlying distribution + if (distribution_.getDimension() == getDimension()) + return distribution_.computePDF(expandPoint(point)); + if (usePDF_) + { + if (isContinuous()) + { + // Build the relevant parametric function to be integrated over the remaining parameters + const ParametricFunction kernel(PDFWrapper(distribution_.getImplementation()->clone()), indices_, point); + const Interval marginalInterval(distribution_.getRange().getMarginal(indices_.complement(distribution_.getDimension()))); + return integrationAlgorithm_.integrate(kernel, marginalInterval)[0]; + } + if (isDiscrete()) + { + const Point probabilities(distribution_.getProbabilities()); + // We modify the support in-place to speed-up the computation + Sample support(distribution_.getImplementation()->getSupport()); + for (UnsignedInteger i = 0; i < support.getSize(); ++i) + for (UnsignedInteger j = 0; j < indices_.getSize(); ++j) + support(i, indices_[j]) = point[j]; + return distribution_.computePDF(support).asPoint().dot(probabilities); + } + } // use PDF + return DistributionImplementation::computePDF(point); +} + /* Get the CDF of the MarginalDistribution */ Scalar MarginalDistribution::computeCDF(const Point & point) const { @@ -270,7 +301,7 @@ CorrelationMatrix MarginalDistribution::getSpearmanCorrelation() const return spearmanCorrelation; } -/* Get the Spearman correlation of the distribution */ +/* Get the Kendall tau of the distribution */ CorrelationMatrix MarginalDistribution::getKendallTau() const { const UnsignedInteger dimension = getDimension(); @@ -306,7 +337,9 @@ Distribution MarginalDistribution::getMarginal(const Indices & indices) const Indices marginalIndices(outputDimension); for (UnsignedInteger i = 0; i < outputDimension; ++i) marginalIndices[i] = indices_[indices[i]]; - return new MarginalDistribution(distribution_, marginalIndices); + MarginalDistribution marginal(distribution_, marginalIndices); + marginal.setIntegrationAlgorithm(getIntegrationAlgorithm()); + return marginal.clone(); } /* Get the isoprobabilistic transformation */ @@ -381,6 +414,48 @@ Bool MarginalDistribution::isIntegral() const return distribution_.isIntegral(); } +/* UsePDF accessor */ +void MarginalDistribution::setUsePDF(const Bool usePDF) +{ + usePDF_ = usePDF; +} + +Bool MarginalDistribution::getUsePDF() const +{ + return usePDF_; +} + +/* AlgoIntegration accessor */ +void MarginalDistribution::setIntegrationAlgorithm(const IntegrationAlgorithm & algo) +{ + integrationAlgorithm_ = algo; +} + +IntegrationAlgorithm MarginalDistribution::getIntegrationAlgorithm() const +{ + return integrationAlgorithm_; +} + +/* Get the support of a discrete distribution that intersect a given interval */ +Sample MarginalDistribution::getSupport(const Interval & interval) const +{ + if (interval.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: the given interval has a dimension that does not match the distribution dimension."; + const Interval extendedInterval(expandPoint(interval.getLowerBound(), false), expandPoint(interval.getUpperBound())); + return distribution_.getSupport(extendedInterval).getMarginal(indices_); +} + +/* Get the support on the whole range */ +Sample MarginalDistribution::getSupport() const +{ + return distribution_.getSupport().getMarginal(indices_); +} + +/* Get the discrete probability levels */ +Point MarginalDistribution::getProbabilities() const +{ + return distribution_.getProbabilities(); +} + /* Method to expand a given point in the marginal space to a point in the underlying distribution space */ Point MarginalDistribution::expandPoint(const Point & point, const Bool upper) const diff --git a/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx index 2b956035a27..c5d670b2417 100644 --- a/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx @@ -23,6 +23,9 @@ #include "openturns/Distribution.hxx" #include "openturns/DistributionImplementation.hxx" +#include "openturns/IntegrationAlgorithm.hxx" +#include "openturns/IteratedQuadrature.hxx" +#include "openturns/GaussKronrod.hxx" BEGIN_NAMESPACE_OPENTURNS @@ -67,7 +70,24 @@ public: void setIndices(const Indices & indices); Indices getIndices() const; -private: + /** UsePDF accessor */ + void setUsePDF(const Bool usePDF); + Bool getUsePDF() const; + + /** Integration algorithm accessor */ + void setIntegrationAlgorithm(const IntegrationAlgorithm & algo); + IntegrationAlgorithm getIntegrationAlgorithm() const; + + /** Get the support of a discrete distribution that intersect a given interval */ + using DistributionImplementation::getSupport; + Sample getSupport() const override; + Sample getSupport(const Interval & interval) const override; + + /** Get the discrete probability levels */ + using DistributionImplementation::getProbabilities; + Point getProbabilities() const override; + + private: /** Set the distribution and the indices in one shot */ void setDistributionAndIndices(const Distribution & distribution, const Indices & indices); @@ -82,6 +102,10 @@ public: Point getRealization() const override; Sample getSample(const UnsignedInteger size) const override; + /** Get the PDF of the MarginalDistribution */ + using DistributionImplementation::computePDF; + Scalar computePDF(const Point & point) const override; + /** Get the CDF of the MarginalDistribution */ using DistributionImplementation::computeCDF; Scalar computeCDF(const Point & point) const override; @@ -177,6 +201,17 @@ private: /** The upper bound of the underlying distribution */ Point upperBound_; + /** Flag to tell if the computations have to be done using the CDF or the PDF + of the underlying distribution */ + Bool usePDF_ = ResourceMap::GetAsBool("MarginalDistribution-UsePDF"); + + /** Integration algorithm used to compute the PDF. By default it is an IteratedQuadrature. */ + IntegrationAlgorithm integrationAlgorithm_ = IteratedQuadrature( + GaussKronrod(ResourceMap::GetAsUnsignedInteger("MarginalDistribution-MaximumSubIntervals"), + ResourceMap::GetAsScalar("MarginalDistribution-MaximumError"), + GaussKronrod::GetRuleFromName(ResourceMap::GetAsString("MarginalDistribution-Rule")) + )); + }; /* class MarginalDistribution */