From 2667bf3d7ed971ee27a7774d18f64fab6332cc70 Mon Sep 17 00:00:00 2001 From: regislebrun Date: Thu, 7 Nov 2024 13:08:00 +0100 Subject: [PATCH] Improved PointConditionalDistribution The following special cases have been implemented: BernsteinCopula, BlockIndependentCopula, BlockIndependentDistribution, Mixture, KernelMixture There is also an overloading of the getIsoProbabilisticTransformation()/getInverseIsoProbabilisticTransfrmation() methods to use the simplified version. The log-normalization factor is computed only if no simplified version has been found. It avoid a potentially costly call to computeLogPDF() using a marginal extracted the generic way (it involves numerical integration in possibly high dimension) --- .../PointConditionalDistribution.cxx | 239 +++++++++++++++--- .../PointConditionalDistribution.hxx | 9 + .../t_PointConditionalDistribution_std.py | 25 ++ 3 files changed, 244 insertions(+), 29 deletions(-) diff --git a/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx b/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx index 27cca80960..75abaf8691 100644 --- a/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx +++ b/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx @@ -28,6 +28,11 @@ #include "openturns/SimplicialCubature.hxx" #include "openturns/CubaIntegration.hxx" #include "openturns/Tuples.hxx" +#include "openturns/Beta.hxx" +#include "openturns/BlockIndependentCopula.hxx" +#include "openturns/BlockIndependentDistribution.hxx" +#include "openturns/Dirichlet.hxx" +#include "openturns/EmpiricalBernsteinCopula.hxx" #include "openturns/JointDistribution.hxx" #include "openturns/KernelMixture.hxx" #include "openturns/Mixture.hxx" @@ -38,6 +43,7 @@ #include "openturns/DomainEvent.hxx" #include "openturns/PlatformInfo.hxx" #include "openturns/GaussKronrod.hxx" +#include "openturns/SobolSequence.hxx" #include "openturns/DistFunc.hxx" #include "openturns/SobolSequence.hxx" @@ -144,6 +150,8 @@ PointConditionalDistribution * PointConditionalDistribution::clone() const return new PointConditionalDistribution(*this); } +namespace +{ class PointConditionalDistributionUBoundEvaluation : public EvaluationImplementation { public: @@ -230,21 +238,79 @@ class PointConditionalDistributionVBoundEvaluation : public EvaluationImplementa }; +} // anonymous namespace + +void PointConditionalDistribution::dispatchConditioning(const Collection & distributions, + Distribution & simplified) const +{ + Collection conditioningIndicesBlocks(distributions.getSize()); + Collection conditioningValuesBlocks(distributions.getSize()); + // Sort both the conditioning indices and values in increasing order + const UnsignedInteger conditioningSize = conditioningIndices_.getSize(); + Collection< std::pair > sortedConditioningPairs(conditioningSize); + for (UnsignedInteger i = 0; i < conditioningSize; ++i) + sortedConditioningPairs[i] = std::pair(conditioningIndices_[i], conditioningValues_[i]); + std::sort(sortedConditioningPairs.begin(), sortedConditioningPairs.end()); + UnsignedInteger beginBlock = 0; + UnsignedInteger endBlock = distributions[0].getDimension(); + UnsignedInteger blockCounter = 0; + for (UnsignedInteger conditioningCounter = 0; conditioningCounter < conditioningSize; ++conditioningCounter) + { + const UnsignedInteger conditioningIndex = sortedConditioningPairs[conditioningCounter].first; + // Find the block the conditioning index belongs to + while (conditioningIndex >= endBlock) + { + ++blockCounter; + // Should never go there + if (blockCounter == distributions.getSize()) + break; + beginBlock = endBlock; + endBlock += distributions[blockCounter].getDimension(); + } + conditioningIndicesBlocks[blockCounter].add(conditioningIndex - beginBlock); + conditioningValuesBlocks[blockCounter].add(sortedConditioningPairs[conditioningCounter].second); + } // for conditioningCounter + // Now, decide what to do for each block + Collection newBlocks; + for (UnsignedInteger i = 0; i < distributions.getSize(); ++i) + { + // If the block is fully conditioned, skip it + if (conditioningIndicesBlocks[i].getSize() == distributions[i].getDimension()) + continue; + const PointConditionalDistribution conditionalBlock(distributions[i], conditioningIndicesBlocks[i], conditioningValuesBlocks[i]); + Distribution simplifiedBlock; + if (conditionalBlock.useSimplifiedVersion_) + newBlocks.add(conditionalBlock.getSimplifiedVersion()); + else + newBlocks.add(conditionalBlock); + } // for i + // Here we return a BlockIndependentDistribution even if the initial collection of distributions was + // made of copulas only, as the conditioning breaks the copula property + simplified = (newBlocks.getSize() == 1 ? newBlocks[0] : BlockIndependentDistribution(newBlocks)); +} + void PointConditionalDistribution::update() { const UnsignedInteger fullDimension = distribution_.getDimension(); - if (conditioningIndices_.getSize()) - logNormalizationFactor_ = distribution_.getMarginal(conditioningIndices_).computeLogPDF(conditioningValues_); - if (!(logNormalizationFactor_ > std::log(getPDFEpsilon()))) - throw InvalidArgumentException(HERE) << "Conditioning vector log PDF value is too low (" << logNormalizationFactor_ << ")"; nonConditioningIndices_ = conditioningIndices_.complement(fullDimension); setDimension(nonConditioningIndices_.getSize()); + if (getDimension() == 0) + throw InvalidArgumentException(HERE) << "Cannot define a conditional distribution by fixing the value of all the components"; setDescription(distribution_.getDescription().select(nonConditioningIndices_)); // enable simplified path useSimplifiedVersion_ = hasSimplifiedVersion(simplifiedVersion_); LOGDEBUG(OSS() << "useSimplifiedVersion_=" << useSimplifiedVersion_); + // We can postpone the computation of the normalization factor here as we will not need it if there is a simplified version (and it can be costly due to the marginal extraction) + if (!useSimplifiedVersion_) + { + if (conditioningIndices_.getSize()) + logNormalizationFactor_ = distribution_.getMarginal(conditioningIndices_).computeLogPDF(conditioningValues_); + if (!(logNormalizationFactor_ > std::log(getPDFEpsilon()))) + throw InvalidArgumentException(HERE) << "Conditioning vector log PDF value is too low (" << logNormalizationFactor_ << ")"; + } + // cache marginal for reuse if (!useSimplifiedVersion_) marginalConditionedDistribution_ = distribution_.getMarginal(nonConditioningIndices_); @@ -293,7 +359,18 @@ void PointConditionalDistribution::update() problemU.setBounds(bounds); OptimizationAlgorithm algo(OptimizationAlgorithm::GetByName(ResourceMap::GetAsString("PointConditionalDistribution-OptimizationAlgorithm"))); algo.setProblem(problemU); - algo.setStartingPoint(start); + const SobolSequence generator(dimension); + Bool found = false; + Point candidate; + while (!found) + { + candidate = lb; + const Point u(generator.generate()); + for (UnsignedInteger i = 0; i < dimension; ++i) + candidate[i] += (ub[i] - lb[i]) * u[i]; + found = computePDF(candidate) > 0.0; + } // found + algo.setStartingPoint(candidate); algo.run(); supU_ = std::exp(algo.getResult().getOptimalValue()[0]); LOGDEBUG(OSS() << "supU_=" << supU_ << " u*=" << algo.getResult().getOptimalPoint()); @@ -312,40 +389,53 @@ void PointConditionalDistribution::update() { problemVI.setBounds(Interval(zero, ub)); algo.setProblem(problemVI); - Point startUb(start); - for (UnsignedInteger j = 0; j < dimension; ++ j) - startUb[j] = std::max(start[j], 0.0); - algo.setStartingPoint(startUb); + Bool found = false; + while (!found) + { + candidate = zero; + const Point u(generator.generate()); + for (UnsignedInteger i = 0; i < dimension; ++i) + candidate[i] = ub[i] * u[i]; + found = computePDF(candidate) > 0.0; + } // found + algo.setStartingPoint(candidate); algo.run(); supV_[i] = std::exp(algo.getResult().getOptimalValue()[0]); - LOGDEBUG(OSS() << "supV_[" << i << "]=" << supV_[i] << " v*=" << algo.getResult().getOptimalPoint()); - } + } // ub[i] > 0.0 if (lb[i] < 0.0) { problemVI.setBounds(Interval(lb, zero)); algo.setProblem(problemVI); - Point startLb(start); - for (UnsignedInteger j = 0; j < dimension; ++ j) - startLb[j] = std::min(startLb[j], 0.0); - algo.setStartingPoint(startLb); + Bool found = false; + while (!found) + { + candidate = zero; + const Point u(generator.generate()); + for (UnsignedInteger i = 0; i < dimension; ++i) + candidate[i] = lb[i] * u[i]; + found = computePDF(candidate) > 0.0; + } // found + algo.setStartingPoint(candidate); algo.run(); infV_[i] = -std::exp(algo.getResult().getOptimalValue()[0]); - LOGDEBUG(OSS() << "infV_[" << i << "]=" << infV_[i] << " v*=" << algo.getResult().getOptimalPoint()); - } + } // lb[i] < 0.0 } } - // cache reordered marginals - Indices indices(conditioningIndices_); - indices.add(nonConditioningIndices_); // initialized in update() - reorderedDistribution_ = distribution_.getMarginal(indices); - - // cache qI - Point x(conditioningValues_); - x.add(getRange().getLowerBound()); - conditioningCDF_ = reorderedDistribution_.computeSequentialConditionalCDF(x); - conditioningCDF_.resize(conditioningIndices_.getSize()); - LOGDEBUG(OSS() << "conditioningCDF_=" << conditioningCDF_); + if (!useSimplifiedVersion_) + { + // cache reordered marginals + Indices indices(conditioningIndices_); + indices.add(nonConditioningIndices_); // initialized in update() + reorderedDistribution_ = distribution_.getMarginal(indices); + + // cache qI + Point x(conditioningValues_); + x.add(getRange().getLowerBound()); + conditioningCDF_ = reorderedDistribution_.computeSequentialConditionalCDF(x); + conditioningCDF_.resize(conditioningIndices_.getSize()); + LOGDEBUG(OSS() << "conditioningCDF_=" << conditioningCDF_); + } } @@ -463,13 +553,86 @@ Bool PointConditionalDistribution::hasSimplifiedVersion(Distribution & simplifie return true; } + // EmpiricalBernsteinCopula + EmpiricalBernsteinCopula *p_empirical_bernstein_copula = dynamic_cast(distribution_.getImplementation().get()); + if (p_empirical_bernstein_copula) + { + const Sample copulaSample(p_empirical_bernstein_copula->getCopulaSample()); + const UnsignedInteger sampleSize = copulaSample.getSize(); + const UnsignedInteger binNumber = p_empirical_bernstein_copula->getBinNumber(); + Collection atoms(sampleSize); + Point weights(sampleSize, 0.0); + const UnsignedInteger dimension = getDimension(); + const UnsignedInteger conditioningDimension = conditioningIndices_.getSize(); + for (UnsignedInteger i = 0; i < sampleSize; ++i) + { + Collection atomComponents(dimension); + for (UnsignedInteger j = 0; j < dimension; ++j) + { + const UnsignedInteger newJ = nonConditioningIndices_[j]; + const Scalar r = std::ceil(binNumber * copulaSample(i, newJ)); + atomComponents[j] = Beta(r, binNumber - r + 1.0, 0.0, 1.0); + } // j + atoms[i] = JointDistribution(atomComponents); + for (UnsignedInteger j = 0; j < conditioningDimension; ++j) + { + const UnsignedInteger newJ = conditioningIndices_[j]; + const Scalar r = std::ceil(binNumber * copulaSample(i, newJ)); + const Scalar xJ = conditioningValues_[j]; + weights[i] += -SpecFunc::LogBeta(r, binNumber - r + 1.0) + (r - 1.0) * std::log(xJ) + (binNumber - r) * std::log1p(-xJ); + } // j + weights[i] = std::exp(weights[i]); + } // i + simplified = Mixture(atoms, weights); + return true; + } + + // Dirichlet + Dirichlet *p_dirichlet = dynamic_cast(distribution_.getImplementation().get()); + if (p_dirichlet) + { + const Point theta(p_dirichlet->getTheta()); + Scalar weight = 1.0; + for (UnsignedInteger i = 0; i < conditioningValues_.getSize(); ++i) + weight -= conditioningValues_[i]; + const UnsignedInteger dimension = getDimension(); + Point newTheta(dimension + 1, theta[p_dirichlet->getDimension()]); + for (UnsignedInteger i = 0; i < nonConditioningIndices_.getSize(); ++i) + { + const UnsignedInteger newI = nonConditioningIndices_[i]; + newTheta[i] =theta[i]; + } + // Unfortunately we don't know how to multiply a multivariate distribution by a scalar + // simplified = Distribution(Dirichlet(newTheta)) * weight; + // return true; + return false; + } + + // BlockIndependentDistribution + std::cerr << "Try to cast into a BlockIndependentDistribution" << std::endl; + BlockIndependentDistribution *p_block_independent_distribution = dynamic_cast(distribution_.getImplementation().get()); + if (p_block_independent_distribution) + { + std::cerr << "Start dispatching" << std::endl; + dispatchConditioning(p_block_independent_distribution->getDistributionCollection(), simplified); + return true; + } + + // BlockIndependentCopula + BlockIndependentCopula *p_block_independent_copula = dynamic_cast(distribution_.getImplementation().get()); + if (p_block_independent_copula) + { + dispatchConditioning(p_block_independent_copula->getCopulaCollection(), simplified); + return true; + } + // Joint JointDistribution *p_joint = dynamic_cast(distribution_.getImplementation().get()); if (p_joint) { const Collection marginals(p_joint->getDistributionCollection()); Point coreConditioniningValues(conditioningIndices_.getSize()); - for (UnsignedInteger i = 0; i < conditioningIndices_.getSize(); ++ i) + for (UnsignedInteger i = 0; i < conditioningIndices_.getSize(); ++i) coreConditioniningValues[i] = marginals[conditioningIndices_[i]].computeCDF(conditioningValues_[i]); const PointConditionalDistribution conditionalCore(p_joint->getCore(), conditioningIndices_, coreConditioniningValues); simplified = JointDistribution(marginals.select(nonConditioningIndices_), conditionalCore); @@ -901,6 +1064,24 @@ Point PointConditionalDistribution::computeSequentialConditionalQuantile(const P return result; } +/* Get the isoprobabilist transformation */ +PointConditionalDistribution::IsoProbabilisticTransformation PointConditionalDistribution::getIsoProbabilisticTransformation() const +{ + if (useSimplifiedVersion_) + return simplifiedVersion_.getIsoProbabilisticTransformation(); + + return DistributionImplementation::getIsoProbabilisticTransformation(); +} + +/* Get the inverse isoprobabilist transformation */ +PointConditionalDistribution::InverseIsoProbabilisticTransformation PointConditionalDistribution::getInverseIsoProbabilisticTransformation() const +{ + if (useSimplifiedVersion_) + return simplifiedVersion_.getInverseIsoProbabilisticTransformation(); + + return DistributionImplementation::getInverseIsoProbabilisticTransformation(); +} + Point PointConditionalDistribution::getConditioningValues() const { return conditioningValues_; diff --git a/lib/src/Uncertainty/Distribution/openturns/PointConditionalDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/PointConditionalDistribution.hxx index e46048ca57..4ae4373330 100644 --- a/lib/src/Uncertainty/Distribution/openturns/PointConditionalDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/PointConditionalDistribution.hxx @@ -101,6 +101,12 @@ public: Scalar computeConditionalQuantile(const Scalar q, const Point & y) const override; Point computeSequentialConditionalQuantile(const Point & q) const override; + /** Get the isoprobabilist transformation */ + IsoProbabilisticTransformation getIsoProbabilisticTransformation() const override; + + /** Get the inverse isoprobabilist transformation */ + InverseIsoProbabilisticTransformation getInverseIsoProbabilisticTransformation() const override; + /** Get the quantile of the distribution */ Point computeQuantile(const Scalar prob, const Bool tail = false) const override; Scalar computeScalarQuantile(const Scalar prob, const Bool tail = false) const override; @@ -185,6 +191,9 @@ private: /** Expand the given marginal point to the underlying distribution argument point */ Point expandPoint(const Point & point) const; + /* Get the simplified version */ + void dispatchConditioning(const Collection & distributions, Distribution & simplified) const; + /* Get the simplified version */ Bool hasSimplifiedVersion(Distribution & simplified) const; diff --git a/python/test/t_PointConditionalDistribution_std.py b/python/test/t_PointConditionalDistribution_std.py index 845fff1665..e6843f6812 100644 --- a/python/test/t_PointConditionalDistribution_std.py +++ b/python/test/t_PointConditionalDistribution_std.py @@ -112,3 +112,28 @@ joint = ot.JointDistribution([ot.Exponential()] * 3, core) distribution = otexp.PointConditionalDistribution(joint, [1], [2.0]) sample = distribution.getSample(10) + +# special case for EmpiricalBernsteinCopula +sample = normal.getSample(10) +bernstein = ot.EmpiricalBernsteinCopula(sample, 5) +distribution = otexp.PointConditionalDistribution(bernstein, [1], [0.5]) +simplified = distribution.getSimplifiedVersion() +print(simplified) +assert simplified.getName() == "Mixture", "wrong type" + +# special case for BlockIndependentDistribution +print("special case for BlockIndependentDistribution") +distributions = [normal, student, mixture] +print("Create BlockIndependentDistribution") +blockIndep = ot.BlockIndependentDistribution(distributions) +# test a conditioning which remove the first block, modify the second block and let the last block unchanged +distribution = otexp.PointConditionalDistribution(blockIndep, [1, 3, 0, 2], [0.5, 0.5, 0.5, 0.5]) +simplified = distribution.getSimplifiedVersion() +print(simplified) +assert simplified.getName() == "BlockIndependentDistribution", "wrong type" +# test a conditioning which remove all but the last block +# As the last block can be simplified, it is its actual simplification which is used +distribution = otexp.PointConditionalDistribution(blockIndep, [0, 1, 2, 3, 4, 5, 6], [0.5]*7) +simplified = distribution.getSimplifiedVersion() +print(simplified) +assert simplified.getName() == "Mixture", "wrong type"