Skip to content

Commit

Permalink
Improved PointConditionalDistribution
Browse files Browse the repository at this point in the history 
The following special cases have been implemented: BernsteinCopula, BlockIndependentCopula, BlockIndependentDistribution, Mixture, KernelMixture
  • Loading branch information
@regislebrun
regislebrun committed Nov 9, 2024
1 parent f1ecda3 commit 4a465b9
Show file tree
Hide file tree
Showing 2 changed files with 170 additions and 14 deletions.
181 changes: 167 additions & 14 deletions lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"

Expand Down Expand Up @@ -144,6 +150,8 @@ PointConditionalDistribution * PointConditionalDistribution::clone() const
return new PointConditionalDistribution(*this);
}

namespace
{
class PointConditionalDistributionUBoundEvaluation : public EvaluationImplementation
{
public:
Expand Down Expand Up @@ -230,6 +238,57 @@ class PointConditionalDistributionVBoundEvaluation : public EvaluationImplementa
};


} // anonymous namespace

void PointConditionalDistribution::dispatchConditioning(const Collection<Distribution> & distributions,
Distribution & simplified) const
{
Collection<Indices> conditioningIndicesBlocks(distributions.getSize());
Collection<Point> conditioningValuesBlocks(distributions.getSize());
// Sort both the conditioning indices and values in increasing order
const UnsignedInteger conditioningSize = conditioningIndices_.getSize();
Collection< std::pair<UnsignedInteger, Scalar> > sortedConditioningPairs(conditioningSize);
for (UnsignedInteger i = 0; i < conditioningSize; ++i)
sortedConditioningPairs[i] = std::pair<UnsignedInteger, Scalar>(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<Distribution> 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 simplified;
if (conditionalBlock.hasSimplifiedVersion(simplified))
newBlocks.add(simplified);
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();
Expand All @@ -239,6 +298,8 @@ void PointConditionalDistribution::update()
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
Expand Down Expand Up @@ -293,7 +354,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());
Expand All @@ -312,26 +384,36 @@ 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
}
}

Expand Down Expand Up @@ -463,13 +545,84 @@ Bool PointConditionalDistribution::hasSimplifiedVersion(Distribution & simplifie
return true;
}

// EmpiricalBernsteinCopula
EmpiricalBernsteinCopula *p_empirical_bernstein_copula = dynamic_cast<EmpiricalBernsteinCopula *>(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<Distribution> atoms(sampleSize);
Point weights(sampleSize, 0.0);
const UnsignedInteger dimension = getDimension();
const UnsignedInteger conditioningDimension = conditioningIndices_.getSize();
for (UnsignedInteger i = 0; i < sampleSize; ++i)
{
Collection<Distribution> 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<Dirichlet *>(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
BlockIndependentDistribution *p_block_independent_distribution = dynamic_cast<BlockIndependentDistribution *>(distribution_.getImplementation().get());
if (p_block_independent_distribution)
{
dispatchConditioning(p_block_independent_distribution->getDistributionCollection(), simplified);
return true;
}

// BlockIndependentCopula
BlockIndependentCopula *p_block_independent_copula = dynamic_cast<BlockIndependentCopula *>(distribution_.getImplementation().get());
if (p_block_independent_copula)
{
dispatchConditioning(p_block_independent_copula->getCopulaCollection(), simplified);
return true;
}

// Joint
JointDistribution *p_joint = dynamic_cast<JointDistribution *>(distribution_.getImplementation().get());
if (p_joint)
{
const Collection<Distribution> 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);
Expand Down
3 changes: 3 additions & 0 deletions lib/src/Uncertainty/Distribution/openturns/PointConditionalDistribution.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,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<Distribution> & distributions, Distribution & simplified) const;

/* Get the simplified version */
Bool hasSimplifiedVersion(Distribution & simplified) const;

Expand Down

0 comments on commit 4a465b9

Please sign in to comment.