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 openturns#2067.
regislebrun · Apr 3, 2024 · 031f919 · 031f919
1 parent 5dbdd61
commit 031f919
Show file tree

Hide file tree

Showing 4 changed files with 127 additions and 5 deletions.
diff --git a/lib/etc/openturns.conf.in b/lib/etc/openturns.conf.in
@@ -542,6 +542,12 @@
   <LogNormalFactory-MaximumIteration  value_int="50" />
   <LogNormalFactory-EstimationMethod  value_int="0"  />
 
+  <!-- MarginalDistribution parameters -->
+  <MarginalDistribution-UsePDF              value_bool="true"         />
+  <MarginalDistribution-MaximumError        value_float="1.0e-7"      />
+  <MarginalDistribution-MaximumSubIntervals value_int="128"           />
+  <MarginalDistribution-Rule                value_str="G15K31"         />
+
   <!-- Meixner parameters -->
   <MeixnerDistribution-MaximumAbsoluteError      value_float="1.0e-12"         />
   <MeixnerDistribution-MaximumConstraintError    value_float="1.0e-12"         />

diff --git 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
@@ -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
@@ -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 */