diff --git a/lib/src/Uncertainty/Distribution/CMakeLists.txt b/lib/src/Uncertainty/Distribution/CMakeLists.txt index d105a448f3..3ad97100da 100644 --- a/lib/src/Uncertainty/Distribution/CMakeLists.txt +++ b/lib/src/Uncertainty/Distribution/CMakeLists.txt @@ -107,6 +107,7 @@ ot_add_source_file (LogUniform.cxx) ot_add_source_file (LogUniformFactory.cxx) ot_add_source_file (MarshallOlkinCopula.cxx) ot_add_source_file (MarginalDistribution.cxx) +ot_add_source_file (MarginalUniformOrderStatistics.cxx) ot_add_source_file (MaximumDistribution.cxx) ot_add_source_file (MaximumEntropyOrderStatisticsDistribution.cxx) ot_add_source_file (MaximumEntropyOrderStatisticsCopula.cxx) @@ -283,6 +284,7 @@ ot_install_header_file (LogUniform.hxx) ot_install_header_file (LogUniformFactory.hxx) ot_install_header_file (MarshallOlkinCopula.hxx) ot_install_header_file (MarginalDistribution.hxx) +ot_install_header_file (MarginalUniformOrderStatistics.hxx) ot_install_header_file (MaximumDistribution.hxx) ot_install_header_file (MaximumEntropyOrderStatisticsDistribution.hxx) ot_install_header_file (MaximumEntropyOrderStatisticsCopula.hxx) diff --git a/lib/src/Uncertainty/Distribution/MarginalUniformOrderStatistics.cxx b/lib/src/Uncertainty/Distribution/MarginalUniformOrderStatistics.cxx new file mode 100644 index 0000000000..f21a615380 --- /dev/null +++ b/lib/src/Uncertainty/Distribution/MarginalUniformOrderStatistics.cxx @@ -0,0 +1,169 @@ +// -*- C++ -*- +/** + * @brief The MarginalUniformOrderStatistics 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 "openturns/MarginalUniformOrderStatistics.hxx" +#include "openturns/PersistentObjectFactory.hxx" +#include "openturns/SimplicialCubature.hxx" +#include "openturns/CubaIntegration.hxx" +#include "openturns/GaussLegendre.hxx" +#include "openturns/IntervalMesher.hxx" +#include "openturns/UniformOrderStatistics.hxx" +#include "openturns/DistFunc.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +CLASSNAMEINIT(MarginalUniformOrderStatistics) + +static const Factory Factory_MarginalUniformOrderStatistics; + + +/* Default constructor */ +MarginalUniformOrderStatistics::MarginalUniformOrderStatistics() + : MarginalDistribution(UniformOrderStatistics(1), {0}) +{ + setName("MarginalUniformOrderStatistics"); +} + +/* Parameters constructor */ +MarginalUniformOrderStatistics::MarginalUniformOrderStatistics(const UnsignedInteger n, + const Indices & indices) + : MarginalDistribution(UniformOrderStatistics(n), indices) +{ + setName("MarginalUniformOrderStatistics"); + // compute the log-PDF normalization factor + // Waiting for a more accurate computation... + const UnsignedInteger m = indices.getSize(); + if (m == 0) throw InvalidArgumentException(HERE) << "Error: cannot build a MarginalUniformOrderStatistics based on an empty Indice"; + logNormalizationFactor_ = SpecFunc::LogGamma(n + 1) - SpecFunc::LogGamma(indices[0] + 1); + for (UnsignedInteger i = 1; i < m; ++i) + logNormalizationFactor_ -= SpecFunc::LogGamma(indices[i] - indices[i - 1]); + logNormalizationFactor_ -= SpecFunc::LogGamma(n - indices[indices.getSize() - 1]); + simplex_ = IntervalMesher(Indices(m, 1)).build(Interval(m)); +} + +/* Comparison operator */ +Bool MarginalUniformOrderStatistics::operator ==(const MarginalUniformOrderStatistics & other) const +{ + if (this == &other) return true; + return (getN() == other.getN()) && (indices_ == other.indices_); +} + +Bool MarginalUniformOrderStatistics::equals(const DistributionImplementation & other) const +{ + const MarginalUniformOrderStatistics* p_other = dynamic_cast(&other); + if (p_other) return (*this == *p_other); + return MarginalDistribution::equals(other); +} + +/* Virtual constructor */ +MarginalUniformOrderStatistics * MarginalUniformOrderStatistics::clone() const +{ + return new MarginalUniformOrderStatistics(*this); +} + +/* Get the PDF of the distribution */ +Scalar MarginalUniformOrderStatistics::computePDF(const Point & point) const +{ + const Scalar logPDF = computeLogPDF(point); + if (logPDF == SpecFunc::LowestScalar) return 0.0; + return std::exp(logPDF); +} + +Scalar MarginalUniformOrderStatistics::computeLogPDF(const Point & point) const +{ + if (point.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << getDimension() << ", here dimension=" << point.getDimension(); + if (!point.isIncreasing()) return SpecFunc::LowestScalar; + // Now, we know that the components of point are in increasing order + if ((point[0] <= 0.0) || (point[dimension_ - 1] >= 1.0)) return SpecFunc::LowestScalar; + // First term, only u_{i0} + Scalar logPDF = logNormalizationFactor_ + indices_[0] * std::log(point[0]); + // Central terms, functions of (u_{i_j}-u_{i_{j-1}}) + for (UnsignedInteger j = 1; j < dimension_; ++j) + logPDF += (indices_[j] - indices_[j - 1] - 1.0) * std::log(point[j] - point[j - 1]); + // Last term, function of u_{i_{n-1}} + logPDF += (getN() - indices_[dimension_ - 1] - 1.0) * std::log1p(-point[dimension_ - 1]); + return logPDF; +} + +/* Get the CDF of the distribution */ +Scalar MarginalUniformOrderStatistics::computeCDF(const Point & point) const +{ + if (point.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << getDimension() << ", here dimension=" << point.getDimension(); + // If the marginal distribution is a permutation of the underlying distribution + if (getN() == dimension_) return MarginalDistribution::computeCDF(point); + // 1D case is trivial + if (dimension_ == 1) + return DistFunc::pBeta(indices_[0] + 1.0, getN() - indices_[0], point[0]); + // Large N case leads to a stack overflow + if (getN() > ResourceMap::GetAsUnsignedInteger("MarginalUniformOrderStatistics-LargeCaseCDF")) + { + // std::cerr << "In the large N case, N=" << getN() << std::endl; + //SimplicialCubature algo; + //algo.setMaximumCallsNumber(1000000); + //const Mesh domain(simplex_.getVertices() * point, simplex_.getSimplices()); + //std::cerr << "domain=" << domain << std::endl; + //return SpecFunc::Clip01(algo.integrate(Function(PDFWrapper(this)), domain)[0]); + CubaIntegration algo("cuhre"); + return SpecFunc::Clip01(algo.integrate(Function(PDFWrapper(this)), Interval(Point(dimension_), point))[0]); + } + // std::cerr << "In the small N case, N=" << getN() << std::endl; + return MarginalDistribution::computeCDF(point); +} + +/* Get the distribution of the marginal distribution corresponding to indices dimensions */ +Distribution MarginalUniformOrderStatistics::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"; + if (dimension == 1) return clone(); + // Build the indices associated to the marginal of the marginal + const UnsignedInteger outputDimension = indices.getSize(); + Indices marginalIndices(outputDimension); + for (UnsignedInteger i = 0; i < outputDimension; ++i) + marginalIndices[i] = indices_[indices[i]]; + MarginalUniformOrderStatistics marginal(distribution_.getDimension(), marginalIndices); + return marginal.clone(); +} // getMarginal(Indices) + +/* N accessor */ +UnsignedInteger MarginalUniformOrderStatistics::getN() const +{ + return distribution_.getDimension(); +} + +/* Method save() stores the object through the StorageManager */ +void MarginalUniformOrderStatistics::save(Advocate & adv) const +{ + MarginalDistribution::save(adv); + adv.saveAttribute("logNormalizationFactor_", logNormalizationFactor_); + adv.saveAttribute("simplex_", simplex_); +} + +/* Method load() reloads the object from the StorageManager */ +void MarginalUniformOrderStatistics::load(Advocate & adv) +{ + MarginalDistribution::load(adv); + adv.loadAttribute("logNormalizationFactor_", logNormalizationFactor_); + adv.loadAttribute("simplex_", simplex_); + computeRange(); +} + + +END_NAMESPACE_OPENTURNS diff --git a/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx index 33d50e4a79..4d95f814f3 100644 --- a/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx @@ -175,7 +175,7 @@ public: void load(Advocate & adv) override; -private: +protected: /** Compute the mean of the distribution */ void computeMean() const override; diff --git a/lib/src/Uncertainty/Distribution/openturns/MarginalUniformOrderStatistics.hxx b/lib/src/Uncertainty/Distribution/openturns/MarginalUniformOrderStatistics.hxx new file mode 100644 index 0000000000..20b5fc552a --- /dev/null +++ b/lib/src/Uncertainty/Distribution/openturns/MarginalUniformOrderStatistics.hxx @@ -0,0 +1,99 @@ +// -*- C++ -*- +/** + * @brief The MarginalUniformOrderStatistics 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_MARGINALUNIFORMORDERSTATISTICS_HXX +#define OPENTURNS_MARGINALUNIFORMORDERSTATISTICS_HXX + +#include "openturns/OTprivate.hxx" +#include "openturns/MarginalDistribution.hxx" +#include "openturns/Mesh.hxx" +#include "openturns/Distribution.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +/** + * @class MarginalUniformOrderStatistics + * + * The MarginalUniformOrderStatistics distribution. + */ +class OT_API MarginalUniformOrderStatistics + : public MarginalDistribution +{ + CLASSNAME + +public: + + /** Default constructor */ + MarginalUniformOrderStatistics(); + + /** Parameters constructor */ + MarginalUniformOrderStatistics(const UnsignedInteger n, + const Indices & indices); + + /** N accessor */ + UnsignedInteger getN() const; + + /** Comparison operator */ + using MarginalDistribution::operator ==; + Bool operator ==(const MarginalUniformOrderStatistics & other) const; +protected: + Bool equals(const DistributionImplementation & other) const override; +public: + + /* Interface inherited from Distribution */ + + /** Virtual constructor */ + MarginalUniformOrderStatistics * clone() const override; + + /** Compute the PDF of the distribution */ + using MarginalDistribution::computePDF; + Scalar computePDF(const Point & point) const override; + + using MarginalDistribution::computeLogPDF; + Scalar computeLogPDF(const Point & point) const override; + + /** Compute the PDF of the distribution */ + using MarginalDistribution::computeCDF; + Scalar computeCDF(const Point & point) const override; + + /** Get the distribution of the marginal distribution corresponding to indices dimensions */ + using MarginalDistribution::getMarginal; + Distribution getMarginal(const Indices & indices) const override; + + /** 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; + +private: + + /** Log-normalization factor for the PDF */ + Scalar logNormalizationFactor_ = 0.0; + + /** Actual domain of the distribution */ + Mesh simplex_; + +}; /* class MarginalUniformOrderStatistics */ + + +END_NAMESPACE_OPENTURNS + +#endif /* OPENTURNS_MARGINALUNIFORMORDERSTATISTICS_HXX */ diff --git a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx index 0f2eefc2c9..09867c0b89 100644 --- a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx @@ -116,6 +116,7 @@ #include "openturns/OrdinalSumCopula.hxx" #include "openturns/MarshallOlkinCopula.hxx" #include "openturns/MarginalDistribution.hxx" +#include "openturns/MarginalUniformOrderStatistics.hxx" #include "openturns/MaximumDistribution.hxx" #include "openturns/MaximumEntropyOrderStatisticsDistribution.hxx" #include "openturns/MaximumEntropyOrderStatisticsCopula.hxx" diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt index e1070831f8..5f01611cf7 100644 --- a/python/src/CMakeLists.txt +++ b/python/src/CMakeLists.txt @@ -632,7 +632,6 @@ 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 @@ -747,7 +746,7 @@ ot_add_python_module(dist_bundle3 dist_bundle3_module.i Uniform.i Uniform_doc.i.in UniformFactory.i UniformFactory_doc.i.in UniformMuSigma.i UniformMuSigma_doc.i.in - UniformOverMesh.i UniformOverMesh_doc.i.in + UniformOverMesh.i UniformOverMesh_doc.i.in UserDefined.i UserDefined_doc.i.in UserDefinedFactory.i UserDefinedFactory_doc.i.in VonMises.i VonMises_doc.i.in @@ -1034,6 +1033,7 @@ ot_add_python_module(experimental experimental_module.i TruncatedOverMesh.i TruncatedOverMesh_doc.i.in PosteriorDistribution.i PosteriorDistribution_doc.i.in UniformOrderStatistics.i UniformOrderStatistics_doc.i.in + MarginalUniformOrderStatistics.i MarginalUniformOrderStatistics_doc.i.in GeneralizedExtremeValueValidation.i GeneralizedExtremeValueValidation_doc.i.in GeneralizedParetoValidation.i GeneralizedParetoValidation_doc.i.in BoundaryMesher.i BoundaryMesher_doc.i.in @@ -1041,6 +1041,7 @@ ot_add_python_module(experimental experimental_module.i RankSobolSensitivityAlgorithm.i RankSobolSensitivityAlgorithm_doc.i.in CubaIntegration.i CubaIntegration_doc.i.in ExperimentIntegration.i ExperimentIntegration_doc.i.in + CombinationsDistribution.i CombinationsDistribution_doc.i.in ) set (OPENTURNS_PYTHON_MODULES ${OPENTURNS_PYTHON_MODULES} PARENT_SCOPE) # for the docstring test diff --git a/python/src/MarginalUniformOrderStatistics.i b/python/src/MarginalUniformOrderStatistics.i new file mode 100644 index 0000000000..ebbbeb2626 --- /dev/null +++ b/python/src/MarginalUniformOrderStatistics.i @@ -0,0 +1,11 @@ +// SWIG file MarginalUniformOrderStatistics.i + +%{ +#include "openturns/MarginalUniformOrderStatistics.hxx" +%} + +%include MarginalUniformOrderStatistics_doc.i + +%copyctor OT::MarginalUniformOrderStatistics; + +%include openturns/MarginalUniformOrderStatistics.hxx diff --git a/python/src/MarginalUniformOrderStatistics_doc.i.in b/python/src/MarginalUniformOrderStatistics_doc.i.in new file mode 100644 index 0000000000..f60aeac318 --- /dev/null +++ b/python/src/MarginalUniformOrderStatistics_doc.i.in @@ -0,0 +1,100 @@ +%feature("docstring") OT::MarginalUniformOrderStatistics +"Marginal uniformOrderStatistics. + +Allows one to extract marginals of a uniformOrderStatistics when it does not natively +support it. + +Parameters +---------- +uniformOrderStatistics : :class:`~openturns.UniformOrderStatistics` + The underlying uniformOrderStatistics. +indices : sequence of int + Marginal indices. +" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::setUniformOrderStatistics +"Accessor to the uniformOrderStatistics. + +Parameters +---------- +uniformOrderStatistics : :class:`~openturns.UniformOrderStatistics` + The underlying uniformOrderStatistics." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::getUniformOrderStatistics +"Accessor to the uniformOrderStatistics. + +Returns +------- +uniformOrderStatistics : :class:`~openturns.UniformOrderStatistics` + The underlying uniformOrderStatistics." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::setIndices +"Accessor to the marginal indices. + +Parameters +---------- +indices : sequence of int + Marginal indices." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::getIndices +"Accessor to the marginal indices. + +Returns +------- +indices : :class:`~openturns.Indices` + Marginal indices." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::setIntegrationAlgorithm +"Accessor to the integration algorithm used to compute the PDF. + +Parameters +---------- +algo : :class:`~openturns.IntegrationAlgorithm` + The integration algorithm used to marginalize the unwanted + components of the PDF." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::getIntegrationAlgorithm +"Accessor to the integration algorithm used to compute the PDF. + +Returns +------- +algo : :class:`~openturns.IntegrationAlgorithm` + The integration algorithm used to marginalize the unwanted + components of the PDF." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::setUsePDF +"Accessor to the flag telling how the PDF is computed. + +Parameters +---------- +flag : bool + Flag telling if the marginal PDF is computed using an integration + of the PDF of the underlying uniformOrderStatistics or if it is computed + using finite differences of the underlying CDF." + +// --------------------------------------------------------------------- + +%feature("docstring") OT::MarginalUniformOrderStatistics::getUsePDF +"Accessor to the flag telling how the PDF is computed. + +Returns +------- +flag : bool + Flag telling if the marginal PDF is computed using an integration + of the PDF of the underlying uniformOrderStatistics or if it is computed + using finite differences of the underlying CDF." + diff --git a/python/src/dist_bundle1_module.i b/python/src/dist_bundle1_module.i index 95bd5bae6a..4297efa083 100644 --- a/python/src/dist_bundle1_module.i +++ b/python/src/dist_bundle1_module.i @@ -41,7 +41,6 @@ %include ChiSquare.i %include ChiSquareFactory.i %include CompositeDistribution.i -%include CombinationsDistribution.i %include Dirac.i %include DiracFactory.i %include Dirichlet.i diff --git a/python/src/experimental_module.i b/python/src/experimental_module.i index 835e9bf7fd..04908c56f9 100644 --- a/python/src/experimental_module.i +++ b/python/src/experimental_module.i @@ -55,31 +55,33 @@ %include simulation_module.i /* Base/Algo */ -%include SimplicialCubature.i %include CubaIntegration.i +%include SimplicialCubature.i /* Uncertainty/Algorithm/Metamodel */ -%include UserDefinedMetropolisHastings.i %include FieldFunctionalChaosResult.i %include FieldToPointFunctionalChaosAlgorithm.i %include FieldFunctionalChaosSobolIndices.i %include PointToFieldFunctionalChaosAlgorithm.i +%include UserDefinedMetropolisHastings.i /* Uncertainty/Algorithm/EventSimulation */ %include CrossEntropyResult.i %include CrossEntropyImportanceSampling.i -%include StandardSpaceCrossEntropyImportanceSampling.i %include PhysicalSpaceCrossEntropyImportanceSampling.i +%include StandardSpaceCrossEntropyImportanceSampling.i /*Uncertainty/Algorithm/Sensitivity */ %include RankSobolSensitivityAlgorithm.i /* Uncertainty/Distribution */ +%include CombinationsDistribution.i +%include GeneralizedExtremeValueValidation.i +%include GeneralizedParetoValidation.i +%include MarginalUniformOrderStatistics.i +%include PosteriorDistribution.i %include SmoothedUniformFactory.i %include StudentCopula.i %include StudentCopulaFactory.i %include TruncatedOverMesh.i -%include PosteriorDistribution.i %include UniformOrderStatistics.i -%include GeneralizedExtremeValueValidation.i -%include GeneralizedParetoValidation.i diff --git a/python/test/t_CombinationsDistribution_std.py b/python/test/t_CombinationsDistribution_std.py new file mode 100755 index 0000000000..0b8178e026 --- /dev/null +++ b/python/test/t_CombinationsDistribution_std.py @@ -0,0 +1,49 @@ +#! /usr/bin/env python + +import openturns as ot + +ot.TESTPREAMBLE() + +# Instantiate one distribution object +distribution = ot.KPermutationsDistribution(5, 12) +print("Distribution ", distribution) + +# Is this distribution elliptical ? +print("Elliptical = ", distribution.isElliptical()) + +# Is this distribution continuous ? +print("Continuous = ", distribution.isContinuous()) + +# Test for realization of distribution +oneRealization = distribution.getRealization() +print("oneRealization=", oneRealization) + +# Test for sampling +size = 10000 +oneSample = distribution.getSample(size) +print("oneSample first=", oneSample[0], " last=", oneSample[size - 1]) +print("mean=", oneSample.computeMean()) +print("covariance=", oneSample.computeCovariance()) +# Define a point +point = ot.Point(distribution.getDimension(), 4.0) +print("Point= ", point) + +# Show PDF and CDF of point +LPDF = distribution.computeLogPDF(point) +assert (LPDF == ot.SpecFunc.LowestScalar) +PDF = distribution.computePDF(point) +print("pdf =", PDF) +CDF = distribution.computeCDF(point) +print("cdf=%.6f" % CDF) +CCDF = distribution.computeComplementaryCDF(point) +print("ccdf=%.6f" % CCDF) +quantile = distribution.computeQuantile(0.95) +print("quantile=", quantile) +print("cdf(quantile)=", distribution.computeCDF(quantile)) +print("entropy=%.6f" % distribution.computeEntropy()) +mean = distribution.getMean() +print("mean=", mean) +covariance = distribution.getCovariance() +print("covariance=", covariance) +parameters = distribution.getParametersCollection() +print("parameters=", parameters)