From 22ca880fa9b69859d050d5b2386983499dfd1104 Mon Sep 17 00:00:00 2001 From: regislebrun Date: Tue, 30 Apr 2024 21:46:07 +0200 Subject: [PATCH] Added the CombinationsDistribution class This multivariate discrete distribution is the uniform distribution over the set of strictly increasing functions from {0,...,k-1} to {0,...,n-1} where k and n are natural numbers such that k<=n. It is also the uniform distribution over the subsets of size k of a set of size n. --- .../Uncertainty/Distribution/CMakeLists.txt | 2 + .../Distribution/CombinationsDistribution.cxx | 468 ++++++++++++++++++ .../openturns/CombinationsDistribution.hxx | 161 ++++++ .../openturns/CombinationsDistributions.hxx | 153 ++++++ .../Distribution/openturns/OTDistribution.hxx | 1 + lib/test/CMakeLists.txt | 1 + lib/test/t_CombinationsDistribution_std.cxx | 114 +++++ .../doc/pyplots/CombinationsDistribution.py | 24 + .../user_manual/probabilistic_modelling.rst | 5 + python/src/CMakeLists.txt | 1 + python/src/CombinationsDistribution.i | 11 + python/src/CombinationsDistribution_doc.i.in | 85 ++++ python/src/dist_bundle1_module.i | 1 + 13 files changed, 1027 insertions(+) create mode 100644 lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx create mode 100644 lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx create mode 100644 lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx create mode 100644 lib/test/t_CombinationsDistribution_std.cxx create mode 100644 python/doc/pyplots/CombinationsDistribution.py create mode 100644 python/src/CombinationsDistribution.i create mode 100644 python/src/CombinationsDistribution_doc.i.in diff --git a/lib/src/Uncertainty/Distribution/CMakeLists.txt b/lib/src/Uncertainty/Distribution/CMakeLists.txt index 4b50220edb..d105a448f3 100644 --- a/lib/src/Uncertainty/Distribution/CMakeLists.txt +++ b/lib/src/Uncertainty/Distribution/CMakeLists.txt @@ -29,6 +29,7 @@ ot_add_source_file (ClaytonCopulaFactory.cxx) ot_add_source_file (BlockIndependentCopula.cxx) ot_add_source_file (JointDistribution.cxx) ot_add_source_file (CompositeDistribution.cxx) +ot_add_source_file (CombinationsDistribution.cxx) ot_add_source_file (ConditionalDistribution.cxx) ot_add_source_file (CumulativeDistributionNetwork.cxx) ot_add_source_file (Dirac.cxx) @@ -209,6 +210,7 @@ ot_install_header_file (ClaytonCopulaFactory.hxx) ot_install_header_file (BlockIndependentCopula.hxx) ot_install_header_file (JointDistribution.hxx) ot_install_header_file (CompositeDistribution.hxx) +ot_install_header_file (CombinationsDistribution.hxx) ot_install_header_file (ConditionalDistribution.hxx) ot_install_header_file (CumulativeDistributionNetwork.hxx) ot_install_header_file (DiracFactory.hxx) diff --git a/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx b/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx new file mode 100644 index 0000000000..6bb34e18b4 --- /dev/null +++ b/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx @@ -0,0 +1,468 @@ +// -*- C++ -*- +/** + * @brief The CombinationsDistribution distribution + * + * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca + * + * This library is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library. If not, see . + * + */ +#include +#include "openturns/Collection.hxx" +#include "openturns/CombinationsDistribution.hxx" +#include "openturns/UserDefined.hxx" +#include "openturns/Combinations.hxx" +#include "openturns/RandomGenerator.hxx" +#include "openturns/SpecFunc.hxx" +#include "openturns/Exception.hxx" +#include "openturns/PersistentObjectFactory.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +typedef Collection UnsignedIntegerCollection; + +CLASSNAMEINIT(CombinationsDistribution) + +static const Factory Factory_CombinationsDistribution; + +/* Default constructor */ +CombinationsDistribution::CombinationsDistribution() + : DiscreteDistribution() + , k_(0) + , n_(0) +{ + setName("CombinationsDistribution"); + setKN(1, 1); +} + +/* Parameters constructor */ +CombinationsDistribution::CombinationsDistribution(const UnsignedInteger k, + const UnsignedInteger n) + : DiscreteDistribution() + , k_(0) + , n_(0) +{ + setName("CombinationsDistribution"); + setKN(k, n); +} + +/* Comparison operator */ +Bool CombinationsDistribution::operator ==(const CombinationsDistribution & other) const +{ + if (this == &other) return true; + return (k_ == other.k_) && (n_ == other.n_); +} + +Bool CombinationsDistribution::equals(const DistributionImplementation & other) const +{ + const CombinationsDistribution* p_other = dynamic_cast(&other); + return p_other && (*this == *p_other); +} + +/* String converter */ +String CombinationsDistribution::__repr__() const +{ + OSS oss; + oss << "class=" << CombinationsDistribution::GetClassName() + << " name=" << getName() + << " dimension=" << getDimension() + << " k=" << k_ + << " n=" << n_; + return oss; +} + +String CombinationsDistribution::__str__(const String & ) const +{ + OSS oss; + oss << getClassName() << "(k = " << k_ << ", n = " << n_ << ")"; + return oss; +} + +/* Virtual constructor */ +CombinationsDistribution * CombinationsDistribution::clone() const +{ + return new CombinationsDistribution(*this); +} + +/* Compute the numerical range of the distribution given the parameters values */ +void CombinationsDistribution::computeRange() +{ + Point lowerBound(k_); + Point upperBound(k_); + for (UnsignedInteger i = 0; i < k_; ++i) + { + lowerBound[i] = i; + upperBound[i] = n_ - k_ + i; + } + const Interval::BoolCollection finiteLowerBound(k_, true); + const Interval::BoolCollection finiteUpperBound(k_, true); + setRange(Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound)); +} + +/* Get one realization of the distribution + See https://cs.stackexchange.com/questions/104930/efficient-n-choose-k-random-sampling +*/ +Point CombinationsDistribution::getRealization() const +{ + Indices flags(n_); + Indices integralRealization; + Bool done = false; + // Build the complementary subset if its cardinal is smaller than k_ + const UnsignedInteger actualK = (k_ > n_ / 2 ? n_ - k_ : k_); + while (integralRealization.getSize() < actualK) + { + const UnsignedInteger i = RandomGenerator::IntegerGenerate(n_); + if (flags[i] == 0) + { + integralRealization.add(i); + flags[i] = 1; + } + } // while + // Di I built the complementary set? + if (actualK != k_) + integralRealization = integralRealization.complement(n_); + // The realization must be sorted in ascending order + std::sort(integralRealization.begin(), integralRealization.end()); + Point realization(k_); + for (UnsignedInteger i = 0; i < k_; ++i) + realization[i] = integralRealization[i]; + return realization; +} + +/* Get the PDF of the distribution */ +Scalar CombinationsDistribution::computeLogPDF(const Point & point) const +{ + if (point.getDimension() != k_) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << k_ << ", here dimension=" << point.getDimension(); + Indices x(k_); + for (UnsignedInteger i = 0; i < k_; ++i) + { + const Scalar k = point[i]; + if ((k < i - supportEpsilon_) || (k > n_ - k_ + i + supportEpsilon_)) return SpecFunc::LowestScalar; + const UnsignedInteger ik = static_cast< UnsignedInteger > (round(k)); + if (std::abs(k - ik) > supportEpsilon_) return SpecFunc::LowestScalar; + x[i] = ik; + } + if (!x.isStrictlyIncreasing()) return SpecFunc::LowestScalar; + return logPDFValue_; +} + +Scalar CombinationsDistribution::computePDF(const Point & point) const +{ + const Scalar logPDF = computeLogPDF(point); + if (logPDF == SpecFunc::LowestScalar) return 0.0; + return std::exp(logPDF); +} + +Scalar CombinationsDistribution::exploreTree(const Scalar j, + const UnsignedInteger lower, + const UnsignedInteger upper, + const UnsignedInteger count, + const Point & xReduced) const +{ + const Scalar normalizationProba = (k_ - j) / (n_ - j); + // Upper branch of the tree + const UnsignedInteger lower1 = lower; + const UnsignedInteger upper1 = upper; + // The value lower1 == dimension is a guard, telling us that the lower bound is 0 and not a component of xReduced + const Scalar a1 = (lower1 < xReduced.getDimension()) ? xReduced[lower1] : -1.0; + const Scalar b1 = xReduced[upper1]; + const UnsignedInteger count1 = count; + const Scalar f1 = (b1 - a1 - count1) / (count1 + 1); + // Lower branch of the tree + const UnsignedInteger lower2 = upper; + const UnsignedInteger upper2 = j; + const Scalar a2 = xReduced[lower2]; + const Scalar b2 = xReduced[upper2]; + const UnsignedInteger count2 = 1; + const Scalar f2 = b2 - a2; + if (j == k_ - 1) + return (f1 + f2) * normalizationProba; + Scalar value = 0.0; + // This test allows one to cut upper parts of the tree + if (f1 > 0.0) + value += f1 * exploreTree(j + 1, lower1, upper1, count1 + 1, xReduced); + // This test allows one to cut lower parts of the tree + if (f2 > 0.0) + value += f2 * exploreTree(j + 1, lower2, upper2, count2, xReduced); + return value * normalizationProba; +} + +/* Get the CDF of the distribution */ +Scalar CombinationsDistribution::computeCDF(const Point & point) const +{ + if (point.getDimension() != k_) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << k_ << ", here dimension=" << point.getDimension(); + + Point sortedPoint(k_); + Scalar oldX = n_ - 1.0; + if (sortedPoint[k_ - 1] < -supportEpsilon_) return 0.0; + for (SignedInteger i = k_ - 1; i >= 0; --i) + { + const Scalar x = std::min(oldX, floor(point[i] + supportEpsilon_)); + if (x < -supportEpsilon_) return 0.0; + sortedPoint[i] = x; + oldX = x - 1.0; + } + // Explore the tree describing the domain of summation + // k = 1 + // lower = dimension_ (guard value to tell that the current interval is [0, xReduced[0]]) + // upper = 0 + // count = 1 + return ((sortedPoint[0] + 1.0) * k_) / n_ * exploreTree(1, k_, 0, 1, sortedPoint); +} + +/* Compute the scalar quantile of the 1D CombinationsDistribution distribution */ +Scalar CombinationsDistribution::computeScalarQuantile(const Scalar prob, + const Bool tail) const +{ + const UnsignedInteger i = static_cast< UnsignedInteger >(ceil(prob * (n_ - 1.0))); + return (tail ? n_ - 1.0 - i : i); +} // computeScalarQuantile + +/* Compute the quantile of the CombinationsDistribution distribution */ +Point CombinationsDistribution::computeQuantile(const Scalar prob, + const Bool tail, + Scalar & marginalProb) const +{ + const Scalar p = (tail ? 1.0 - prob : prob); + if (p <= 0.0) return Point(k_, 0.0); + if (p >= 1.0) return Point(k_, n_); + UnsignedInteger iMin = 0; + UnsignedInteger iMax = n_; + while (iMax > iMin + 1) + { + const UnsignedInteger iMiddle = (iMax + iMin) / 2; + const Scalar cdfMiddle = computeCDF(Point(k_, iMiddle)); + if (cdfMiddle < p) + { + iMin = iMiddle; + } + else + { + iMax = iMiddle; + } + } // while + marginalProb = computeScalarQuantile(prob, tail); + return Point(k_, iMax); +} // computeQuantile + +/* Get the i-th marginal distribution */ +Distribution CombinationsDistribution::getMarginal(const UnsignedInteger i) const +{ + if (i >= k_) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]"; + Sample support(n_ - k_ + 1, 1); + // First, compute the probabilities on a log scale + Point probabilities(n_ - k_ + 1, -SpecFunc::LogBinomialCoefficient(k_, n_)); + for (UnsignedInteger x = i; x <= n_ - k_ + i; ++x) + { + support(x - i, 0) = x; + probabilities[x - i] += SpecFunc::LogBinomialCoefficient(x, i) + SpecFunc::LogBinomialCoefficient(n_ - 1 - x, k_ - 1 - i); + } + // Then, go back to the [0, 1] interval + for (UnsignedInteger j = 0; j <= n_ - k_; ++j) + probabilities[j] = SpecFunc::Clip01(std::exp(probabilities[j])); + UserDefined marginal(support, probabilities); + marginal.setDescription(Description(1, getDescription()[i])); + return marginal; +} + +/* Get the distribution of the marginal distribution corresponding to indices dimensions */ +Distribution CombinationsDistribution::getMarginal(const Indices & indices) const +{ + const UnsignedInteger dimension = getDimension(); + if (!indices.check(dimension)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and must be different"; + // Special case for dimension 1 + if (dimension == 1) return clone(); + // Special case for indices of length 1 + if (indices.getSize() == 1) return getMarginal(indices[0]); + // General case + const UnsignedInteger outputDimension = indices.getSize(); + const Sample support(getSupport().getMarginal(indices)); + Point probabilities(support.getSize()); + for (UnsignedInteger j = 0; j < support.getSize(); ++j) + { + const UnsignedInteger j0 = indices[0]; + const UnsignedInteger x0 = static_cast(support(j, 0)); + const UnsignedInteger je = indices[outputDimension - 1]; + const UnsignedInteger xe = static_cast(support(j, outputDimension - 1)); + probabilities[j] += SpecFunc::LogBinomialCoefficient(x0, j0) + SpecFunc::LogBinomialCoefficient(n_ - xe - 1, k_ - je - 1); + for (UnsignedInteger i = 1; i < outputDimension - 1; ++i) + probabilities[j] += SpecFunc::LogBinomialCoefficient(xe, je); + probabilities[j] = SpecFunc::Clip01(std::exp(probabilities[j])); + } + UserDefined marginal(support, probabilities); + marginal.setDescription(getDescription().select(indices)); + return marginal; +} + +/* Get the support of a discrete distribution that intersect a given interval */ +Sample CombinationsDistribution::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."; + // Convert int values into float + const IndicesCollection intResult(Combinations(k_, n_).generate()); + const UnsignedInteger size = intResult.getSize(); + if (size == 0) return Sample(); + const Interval inter(interval.intersect(range_)); + // Common case: get the full support + if (inter == range_) + { + Sample result(size, dimension_); + for (UnsignedInteger i = 0; i < size; ++i) + for (UnsignedInteger j = 0; j < dimension_; ++j) + result(i, j) = intResult(i, j); + return result; + } + Sample result(0, dimension_); + for (UnsignedInteger i = 0; i < size; ++i) + { + Point point(dimension_); + for (UnsignedInteger j = 0; j < dimension_; ++j) + point[j] = intResult(i, j); + if (inter.contains(point)) + result.add(point); + } + return result; +} + +/* Compute the mean of the distribution */ +void CombinationsDistribution::computeMean() const +{ + mean_ = Point(dimension_); + for (UnsignedInteger i = 0; i < dimension_; ++i) + mean_[i] = getMarginal(i).getMean()[0]; + isAlreadyComputedMean_ = true; +} + +/* Compute the covariance of the distribution */ +void CombinationsDistribution::computeCovariance() const +{ + covariance_ = CovarianceMatrix(dimension_); + for (UnsignedInteger j = 0; j < dimension_; ++j) + { + for (UnsignedInteger i = 0; i < j; ++i) + { + const CovarianceMatrix covIJ(getMarginal({i, j}).getCovariance()); + covariance_(i, j) = covIJ(0, 1); + covariance_(i, i) = covIJ(0, 0); + } + } + isAlreadyComputedCovariance_ = true; +} + +/* Parameters value and description accessor */ +CombinationsDistribution::PointWithDescriptionCollection CombinationsDistribution::getParametersCollection() const +{ + const UnsignedInteger dimension = getDimension(); + PointWithDescriptionCollection parameters((dimension == 1 ? 1 : dimension + 1)); + for (UnsignedInteger i = 0; i < dimension; ++i) + { + PointWithDescription point(1); + point[0] = n_; + Description description(1); + description[0] = "n"; + point.setDescription(description); + point.setName(getDescription()[i]); + parameters[i] = point; + } + if (dimension > 1) + { + PointWithDescription point(2); + Description description = {"k", "n"}; + point[0] = k_; + point[1] = n_; + point.setDescription(description); + point.setName("dependence"); + parameters[dimension] = point; + } + return parameters; +} + +/* K accessor */ +void CombinationsDistribution::setK(const UnsignedInteger k) +{ + if (k == 0) throw InvalidArgumentException(HERE) << "Error: k must be > 0."; + if (k > n_) throw InvalidArgumentException(HERE) << "Error: k must be less or equal to n, here k=" << k << " and n=" << n_; + if (k != k_) + { + k_ = k; + logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1); + setDimension(k); + isAlreadyComputedMean_ = false; + isAlreadyComputedCovariance_ = false; + isAlreadyCreatedGeneratingFunction_ = false; + } +} + +/* K accessor */ +UnsignedInteger CombinationsDistribution::getK() const +{ + return k_; +} + +/* N accessor */ +void CombinationsDistribution::setN(const UnsignedInteger n) +{ + if (n == 0) throw InvalidArgumentException(HERE) << "Error: n must be > 0."; + if (n < k_) throw InvalidArgumentException(HERE) << "Error: n must be greater or equal to k, here n=" << n << " and k=" << k_; + if (n != n_) + { + n_ = n; + logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1); + isAlreadyComputedMean_ = false; + isAlreadyComputedCovariance_ = false; + computeRange(); + } +} + +UnsignedInteger CombinationsDistribution::getN() const +{ + return n_; +} + +/* K/N accessor */ +void CombinationsDistribution::setKN(const UnsignedInteger k, + const UnsignedInteger n) +{ + if (k == 0) throw InvalidArgumentException(HERE) << "Error: k must be > 0."; + if (n == 0) throw InvalidArgumentException(HERE) << "Error: n must be > 0."; + if (k > n) throw InvalidArgumentException(HERE) << "Error: k must be less or equal to n, here k=" << k << " and n=" << n; + k_ = k; + setDimension(k); + n_ = n; + logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1); + isAlreadyComputedMean_ = false; + isAlreadyComputedCovariance_ = false; + computeRange(); +} + +/* Method save() stores the object through the StorageManager */ +void CombinationsDistribution::save(Advocate & adv) const +{ + DiscreteDistribution::save(adv); + adv.saveAttribute( "k_", k_ ); + adv.saveAttribute( "n_", n_ ); + adv.saveAttribute( "logPDFValue_", logPDFValue_ ); +} + +/* Method load() reloads the object from the StorageManager */ +void CombinationsDistribution::load(Advocate & adv) +{ + DiscreteDistribution::load(adv); + adv.loadAttribute( "k_", k_ ); + adv.loadAttribute( "n_", n_ ); + adv.loadAttribute( "logPDFValue_", logPDFValue_ ); + computeRange(); +} + +END_NAMESPACE_OPENTURNS diff --git a/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx new file mode 100644 index 0000000000..5fa093a55c --- /dev/null +++ b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx @@ -0,0 +1,161 @@ +// -*- C++ -*- +/** + * @brief The CombinationsDistribution distribution + * + * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca + * + * This library is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library. If not, see . + * + */ +#ifndef OPENTURNS_COMBINATIONSDISTRIBUTION_HXX +#define OPENTURNS_COMBINATIONSDISTRIBUTION_HXX + +#include "openturns/OTprivate.hxx" +#include "openturns/DiscreteDistribution.hxx" +#include "openturns/Distribution.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +/** + * @class CombinationsDistribution + * + * The CombinationsDistribution distribution. + */ +class OT_API CombinationsDistribution + : public DiscreteDistribution +{ + CLASSNAME +public: + + /** Default constructor */ + CombinationsDistribution(); + + + /** Parameters constructor */ + CombinationsDistribution(const UnsignedInteger k, + const UnsignedInteger n); + + + /** Comparison operator */ + using DiscreteDistribution::operator ==; + Bool operator ==(const CombinationsDistribution & other) const; +protected: + Bool equals(const DistributionImplementation & other) const override; +public: + + /** String converter */ + String __repr__() const override; + String __str__(const String & offset = "") const override; + + + + /* Interface inherited from Distribution */ + + /** Virtual constructor */ + CombinationsDistribution * clone() const override; + + /** Get one realization of the distribution */ + Point getRealization() const override; + + /** Get the PDF of the distribution */ + using DiscreteDistribution::computePDF; + Scalar computePDF(const Point & point) const override; + using DiscreteDistribution::computeLogPDF; + Scalar computeLogPDF(const Point & point) const override; + + /** Get the CDF of the distribution */ + using DiscreteDistribution::computeCDF; + Scalar computeCDF(const Point & point) const override; + + /** Get the quantile of the distribution */ + using DiscreteDistribution::computeQuantile; + Point computeQuantile(const Scalar prob, + const Bool tail, + Scalar & marginalProb) const override; + + /** Get the i-th marginal distribution */ + using DiscreteDistribution::getMarginal; + Distribution getMarginal(const UnsignedInteger i) const override; + + /** Get the distribution of the marginal distribution corresponding to indices dimensions */ + Distribution getMarginal(const Indices & indices) const override; + + /** Get the support of a discrete distribution that intersect a given interval */ + using DistributionImplementation::getSupport; + Sample getSupport(const Interval & interval) const override; + + /** Parameters value and description accessor */ + PointWithDescriptionCollection getParametersCollection() const override; + + /* Interface specific to CombinationsDistribution */ + + /** K accessor */ + void setK(const UnsignedInteger k); + UnsignedInteger getK() const; + + /** N accessor */ + void setN(const UnsignedInteger n); + UnsignedInteger getN() const; + +private: + /** K/N accessor */ + void setKN(const UnsignedInteger k, + const UnsignedInteger n); + +public: + /** Method save() stores the object through the StorageManager */ + void save(Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(Advocate & adv) override; + +protected: + + +private: + + /** Compute the numerical range of the distribution given the parameters values */ + void computeRange() override; + + /** Quantile computation for dimension=1 */ + Scalar computeScalarQuantile(const Scalar prob, const Bool tail = false) const override; + + /** Compute the mean of the distribution */ + void computeMean() const override; + + /** Compute the covariance of the distribution */ + void computeCovariance() const override; + + /** Helper to compute the CDF */ + Scalar exploreTree(const Scalar j, + const UnsignedInteger lower, + const UnsignedInteger upper, + const UnsignedInteger count, + const Point & xReduced) const; + + /** Size of the permutations */ + UnsignedInteger k_; + + /** Size of the base set */ + UnsignedInteger n_; + + /** Log PDF value */ + Scalar logPDFValue_; + +}; /* class CombinationsDistribution */ + + +END_NAMESPACE_OPENTURNS + +#endif /* OPENTURNS_KPERMUTATIONSDISTRIBUTION_HXX */ diff --git a/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx new file mode 100644 index 0000000000..2555eae45a --- /dev/null +++ b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx @@ -0,0 +1,153 @@ +// -*- C++ -*- +/** + * @brief The CombinationsDistribution distribution + * + * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca + * + * This library is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library. If not, see . + * + */ +#ifndef OPENTURNS_COMBINATIONSDISTRIBUTION_HXX +#define OPENTURNS_COMBINATIONSDISTRIBUTION_HXX + +#include "openturns/OTprivate.hxx" +#include "openturns/DiscreteDistribution.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +/** + * @class CombinationsDistribution + * + * The CombinationsDistribution distribution. + */ +class OT_API CombinationsDistribution + : public DiscreteDistribution +{ + CLASSNAME +public: + + /** Default constructor */ + CombinationsDistribution(); + + + /** Parameters constructor */ + CombinationsDistribution(const UnsignedInteger k, + const UnsignedInteger n); + + + /** Comparison operator */ + using DiscreteDistribution::operator ==; + Bool operator ==(const CombinationsDistribution & other) const; +protected: + Bool equals(const DistributionImplementation & other) const override; +public: + + /** String converter */ + String __repr__() const override; + String __str__(const String & offset = "") const override; + + + + /* Interface inherited from Distribution */ + + /** Virtual constructor */ + CombinationsDistribution * clone() const override; + + /** Get one realization of the distribution */ + Point getRealization() const override; + + /** Get the PDF of the distribution */ + using DiscreteDistribution::computePDF; + Scalar computePDF(const Point & point) const override; + using DiscreteDistribution::computeLogPDF; + Scalar computeLogPDF(const Point & point) const override; + + /** Get the CDF of the distribution */ + using DiscreteDistribution::computeCDF; + Scalar computeCDF(const Point & point) const override; + + /** Get the quantile of the distribution */ + using DiscreteDistribution::computeQuantile; + Point computeQuantile(const Scalar prob, + const Bool tail, + Scalar & marginalProb) const override; + + /** Get the i-th marginal distribution */ + using DiscreteDistribution::getMarginal; + Distribution getMarginal(const UnsignedInteger i) const override; + + /** Get the distribution of the marginal distribution corresponding to indices dimensions */ + Distribution getMarginal(const Indices & indices) const override; + + /** Get the support of a discrete distribution that intersect a given interval */ + using DistributionImplementation::getSupport; + Sample getSupport(const Interval & interval) const override; + + /** Parameters value and description accessor */ + PointWithDescriptionCollection getParametersCollection() const override; + + /* Interface specific to CombinationsDistribution */ + + /** K accessor */ + void setK(const UnsignedInteger k); + UnsignedInteger getK() const; + + /** N accessor */ + void setN(const UnsignedInteger n); + UnsignedInteger getN() const; + +private: + /** K/N accessor */ + void setKN(const UnsignedInteger k, + const UnsignedInteger n); + +public: + /** Method save() stores the object through the StorageManager */ + void save(Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(Advocate & adv) override; + +protected: + + +private: + + /** Compute the numerical range of the distribution given the parameters values */ + void computeRange() override; + + /** Quantile computation for dimension=1 */ + Scalar computeScalarQuantile(const Scalar prob, const Bool tail = false) const override; + + /** Compute the mean of the distribution */ + void computeMean() const override; + + /** Compute the covariance of the distribution */ + void computeCovariance() const override; + + /** Size of the permutations */ + UnsignedInteger k_; + + /** Size of the base set */ + UnsignedInteger n_; + + /** Log PDF value */ + Scalar logPDFValue_; + +}; /* class CombinationsDistribution */ + + +END_NAMESPACE_OPENTURNS + +#endif /* OPENTURNS_KPERMUTATIONSDISTRIBUTION_HXX */ diff --git a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx index eaf18b0ac5..0f2eefc2c9 100644 --- a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx @@ -45,6 +45,7 @@ #include "openturns/ChiFactory.hxx" #include "openturns/ClaytonCopula.hxx" #include "openturns/ClaytonCopulaFactory.hxx" +#include "openturns/CombinationsDistribution.hxx" #include "openturns/BlockIndependentCopula.hxx" #include "openturns/JointDistribution.hxx" #include "openturns/CompositeDistribution.hxx" diff --git a/lib/test/CMakeLists.txt b/lib/test/CMakeLists.txt index ccf74e675f..f8a5b3e7b8 100644 --- a/lib/test/CMakeLists.txt +++ b/lib/test/CMakeLists.txt @@ -356,6 +356,7 @@ ot_check_test (ChiSquareFactory_std) ot_check_test (BlockIndependentCopula_std) ot_check_test (JointDistribution_std) ot_check_test (JointDistribution_large) +ot_check_test (CombinationsDistribution_std IGNOREOUT) ot_check_test (CompositeDistribution_std) ot_check_test (ConditionalDistribution_std) ot_check_test (CumulativeDistributionNetwork_std) diff --git a/lib/test/t_CombinationsDistribution_std.cxx b/lib/test/t_CombinationsDistribution_std.cxx new file mode 100644 index 0000000000..71b6df5246 --- /dev/null +++ b/lib/test/t_CombinationsDistribution_std.cxx @@ -0,0 +1,114 @@ +// -*- C++ -*- +/** + * @brief The test file of class CombinationsDistribution for standard methods + * + * Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca + * + * This library is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library. If not, see . + * + */ +#include "openturns/OT.hxx" +#include "openturns/OTtestcode.hxx" + +using namespace OT; +using namespace OT::Test; + +class TestObject : public CombinationsDistribution +{ +public: + TestObject() : CombinationsDistribution(5, 12) {} + virtual ~TestObject() {} +}; + + +int main(int, char *[]) +{ + TESTPREAMBLE; + OStream fullprint(std::cout); + setRandomGenerator(); + try + { + // Test basic functionnalities + checkClassWithClassName(); + + // Instantiate one distribution object + CombinationsDistribution distribution(5, 12); + fullprint << "Distribution " << distribution << std::endl; + std::cout << "Distribution " << distribution << std::endl; + UserDefined ref(distribution.getSupport()); + // Is this distribution elliptical ? + fullprint << "Elliptical = " << (distribution.isElliptical() ? "true" : "false") << std::endl; + + // Is this distribution continuous ? + fullprint << "Continuous = " << (distribution.isContinuous() ? "true" : "false") << std::endl; + + // Test for realization of distribution + Point oneRealization = distribution.getRealization(); + fullprint << "oneRealization=" << oneRealization << std::endl; + + // Test for sampling + UnsignedInteger size = 10000; + Sample oneSample = distribution.getSample( size ); + fullprint << "oneSample first=" << oneSample[0] << " last=" << oneSample[size - 1] << std::endl; + fullprint << "mean=" << oneSample.computeMean() << std::endl; + fullprint << "covariance=" << oneSample.computeCovariance() << std::endl; + // Define a point + Point point({1, 3, 6, 8, 10}); + fullprint << "Point= " << point << std::endl; + + // Show PDF and CDF of point + Scalar LPDF = distribution.computeLogPDF( point ); + Scalar PDF = distribution.computePDF( point ); + fullprint << "pdf =" << PDF << std::endl; + fullprint << "pdf ref =" << ref.computePDF(point) << std::endl; + Scalar CDF = distribution.computeCDF( point ); + fullprint << "cdf =" << CDF << std::endl; + fullprint << "cdf ref =" << ref.computeCDF(point) << std::endl; + Scalar CCDF = distribution.computeComplementaryCDF( point ); + fullprint << "ccdf =" << CCDF << std::endl; + fullprint << "ccdf ref =" << ref.computeComplementaryCDF(point) << std::endl; + // Scalar Survival = distribution.computeSurvivalFunction( point ); + // fullprint << "survival=" << Survival << std::endl; + Point quantile = distribution.computeQuantile( 0.95 ); + fullprint << "quantile =" << quantile << std::endl; + fullprint << "quantile ref =" << ref.computeQuantile(0.95) << std::endl; + fullprint << "cdf(quantile)=" << distribution.computeCDF(quantile) << std::endl; + fullprint << "entropy =" << distribution.computeEntropy() << std::endl; + fullprint << "entropy ref =" << ref.computeEntropy() << std::endl; + fullprint << "entropy (MC)=" << -distribution.computeLogPDF(distribution.getSample(1000000)).computeMean()[0] << std::endl; + Point mean = distribution.getMean(); + fullprint << "mean =" << mean << std::endl; + fullprint << "mean ref =" << ref.getMean() << std::endl; + CovarianceMatrix covariance = distribution.getCovariance(); + fullprint << "covariance =" << covariance << std::endl; + fullprint << "covariance ref =" << ref.getCovariance() << std::endl; + CovarianceMatrix correlation = distribution.getCorrelation(); + fullprint << "correlation =" << correlation << std::endl; + fullprint << "correlation ref =" << ref.getCorrelation() << std::endl; + // CovarianceMatrix spearman = distribution.getSpearmanCorrelation(); + // fullprint << "spearman=" << spearman << std::endl; + // CovarianceMatrix kendall = distribution.getKendallTau(); + // fullprint << "kendall=" << kendall << std::endl; + CombinationsDistribution::PointWithDescriptionCollection parameters = distribution.getParametersCollection(); + fullprint << "parameters=" << parameters << std::endl; + } + catch (TestFailed & ex) + { + std::cerr << ex << std::endl; + return ExitCode::Error; + } + + + return ExitCode::Success; +} diff --git a/python/doc/pyplots/CombinationsDistribution.py b/python/doc/pyplots/CombinationsDistribution.py new file mode 100644 index 0000000000..1f731509c0 --- /dev/null +++ b/python/doc/pyplots/CombinationsDistribution.py @@ -0,0 +1,24 @@ +import openturns as ot +import openturns.experimental as otexp +from openturns.viewer import View + +grid = ot.GridLayout(1, 2) +pdf_2d = ot.Graph("CombinationsDistribution, PDF", "x1", "x2", True) +cdf_2d = ot.Graph("CombinationsDistribution, CDF", "x1", "x2", True) + +distribution_2d = otexp.CombinationsDistribution(2, 10) + +xMin = [-0.1] * 2 +xMax = [10.1] * 2 + +pdf_2d.add(distribution_2d.drawPDF(xMin, xMax)) +pdf_2d.setLegends([""]) +cdf_2d.add(distribution_2d.drawCDF(xMin, xMax, [101] * 2)) +cdf_2d.setLegends([""]) +grid.setGraph(0, 0, pdf_2d) +grid.setGraph(0, 1, cdf_2d) +grid.setTitle("CombinationsDistribution") +grid.setLegendPosition("upper right") +v = View(grid) +fig = v.getFigure() +fig.axes[1].legend(loc="best") diff --git a/python/doc/user_manual/probabilistic_modelling.rst b/python/doc/user_manual/probabilistic_modelling.rst index c1e15eeecc..27e7387602 100644 --- a/python/doc/user_manual/probabilistic_modelling.rst +++ b/python/doc/user_manual/probabilistic_modelling.rst @@ -111,6 +111,11 @@ Discrete parametric distributions Bernoulli Binomial + + :template: classWithPlot.rst_t + experimental.CombinationsDistribution + + :template: Distribution.rst_t Dirac Geometric Hypergeometric diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt index c7a87967f0..e1070831f8 100644 --- a/python/src/CMakeLists.txt +++ b/python/src/CMakeLists.txt @@ -632,6 +632,7 @@ ot_add_python_module(dist_bundle1 dist_bundle1_module.i ChiSquareFactory.i ChiSquareFactory_doc.i.in Chi.i Chi_doc.i.in ChiFactory.i ChiFactory_doc.i.in + CombinationsDistribution.i CombinationsDistribution_doc.i.in CompositeDistribution.i CompositeDistribution_doc.i.in Dirac.i Dirac_doc.i.in DiracFactory.i DiracFactory_doc.i.in diff --git a/python/src/CombinationsDistribution.i b/python/src/CombinationsDistribution.i new file mode 100644 index 0000000000..0d595ebec4 --- /dev/null +++ b/python/src/CombinationsDistribution.i @@ -0,0 +1,11 @@ +// SWIG file CombinationsDistribution.i + +%{ +#include "openturns/CombinationsDistribution.hxx" +%} + +%include CombinationsDistribution_doc.i + +%copyctor OT::CombinationsDistribution; + +%include openturns/CombinationsDistribution.hxx diff --git a/python/src/CombinationsDistribution_doc.i.in b/python/src/CombinationsDistribution_doc.i.in new file mode 100644 index 0000000000..84f07dc2f7 --- /dev/null +++ b/python/src/CombinationsDistribution_doc.i.in @@ -0,0 +1,85 @@ +%feature("docstring") OT::CombinationsDistribution +"Combinations distribution. + +Parameters +---------- +k : int, :math:`k > 0` +n : int, :math:`n > 0` + +See also +-------- +:class:`~openturns.Combinations` + +Notes +----- +:class:`~openturns.experimental.CombinationsDistribution` is the discrete uniform +distribution on the set of strictly increasing functions :math:`(i_0, \hdots, i_{k_1})` +from :math:`\{0, \hdots, k-1\}` into :math:`\{0, \hdots, n-1\}`. +Its probability distribution function is defined as: + +.. math:: + + \Prob{\vect{X} = (i_0, \hdots, i_{k-1})} = \frac{1}{\binom{n}{k}} + +Its first moments are: + +.. math:: + :nowrap: + + \begin{eqnarray*} + \Expect{\vect{X}} & = & \frac{n - 1}{2}\\ + \Cov{\vect{X}} & = & \left\{ + \begin{array}{ll} + \displaystyle \frac{n^2-1}{12} \\ + \displaystyle -\frac{1 + n}{12} + \end{array} + \right. + \end{eqnarray*} + +Examples +-------- +Create a distribution: + +>>> import openturns as ot +>>> import openturns.experimental as otexp +>>> distribution = otexp.CombinationsDistribution(6, 8) + +Draw a sample: + +>>> sample = distribution.getSample(10)" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::CombinationsDistribution::getN +"Accessor to the parameter :math:`n`. + +Returns +------- +n : int" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::CombinationsDistribution::getK +"Accessor to the parameter :math:`k`. + +Returns +------- +k : int" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::CombinationsDistribution::setN +"Accessor to the parameter :math:`n`. + +Parameters +---------- +n : int, :math:`n > 0`" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::CombinationsDistribution::setK +"Accessor to the parameter :math:`k`. + +Parameters +---------- +k : int, :math:`k > 0`" diff --git a/python/src/dist_bundle1_module.i b/python/src/dist_bundle1_module.i index 4297efa083..95bd5bae6a 100644 --- a/python/src/dist_bundle1_module.i +++ b/python/src/dist_bundle1_module.i @@ -41,6 +41,7 @@ %include ChiSquare.i %include ChiSquareFactory.i %include CompositeDistribution.i +%include CombinationsDistribution.i %include Dirac.i %include DiracFactory.i %include Dirichlet.i