diff --git a/CMakeLists.txt b/CMakeLists.txt index 893ba9b..e2eabc1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -138,7 +138,7 @@ endmacro ( ot_add_current_dir_to_include_dirs ) set ( CPACK_PACKAGE_NAME ${PACKAGE_NAME} ) set ( CPACK_PACKAGE_VERSION_MAJOR 0 ) -set ( CPACK_PACKAGE_VERSION_MINOR 12 ) +set ( CPACK_PACKAGE_VERSION_MINOR 13 ) set ( CPACK_PACKAGE_VERSION_PATCH ) set ( CPACK_SOURCE_GENERATOR "TGZ;TBZ2" ) set (CPACK_BINARY_STGZ OFF CACHE BOOL "STGZ") diff --git a/ChangeLog b/ChangeLog index e4bf413..59c3c90 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ += 0.13 release (in progress) + + * Add SubsetInverseSampling + * Add InverseFORM + = 0.12 release (2023-05-17) * Maintenance diff --git a/VERSION b/VERSION index c43e105..f304084 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.12 +0.13 diff --git a/distro/debian/changelog b/distro/debian/changelog index 4f9301d..547d452 100644 --- a/distro/debian/changelog +++ b/distro/debian/changelog @@ -1,4 +1,4 @@ -otrobopt (0.12-0.1) unstable; urgency=low +otrobopt (0.13-0.1) unstable; urgency=low * Non-maintainer upload. * Initial release. diff --git a/distro/rpm/otrobopt.spec b/distro/rpm/otrobopt.spec index 7232663..903724e 100644 --- a/distro/rpm/otrobopt.spec +++ b/distro/rpm/otrobopt.spec @@ -20,7 +20,7 @@ FFLAGS="${FFLAGS:-%optflags}" ; export FFLAGS ; \ -DBUILD_SHARED_LIBS:BOOL=ON Name: otrobopt -Version: 0.12 +Version: 0.13 Release: 0%{?dist} Summary: OpenTURNS module Group: System Environment/Libraries diff --git a/lib/include/otrobopt/OTRobOpt.hxx b/lib/include/otrobopt/OTRobOpt.hxx index faf47a3..8c89d71 100644 --- a/lib/include/otrobopt/OTRobOpt.hxx +++ b/lib/include/otrobopt/OTRobOpt.hxx @@ -32,5 +32,7 @@ #include "otrobopt/QuantileMeasure.hxx" #include "otrobopt/AggregatedMeasure.hxx" #include "otrobopt/MeasureFactory.hxx" +#include "otrobopt/SubsetInverseSampling.hxx" +#include "otrobopt/InverseFORM.hxx" #endif /* OTROBOPT_OTROBOPT_HXX */ diff --git a/lib/src/CMakeLists.txt b/lib/src/CMakeLists.txt index e17b232..53d5159 100644 --- a/lib/src/CMakeLists.txt +++ b/lib/src/CMakeLists.txt @@ -17,6 +17,10 @@ ot_add_source_file (MeasureFactory.cxx) ot_add_source_file (RobustOptimizationProblem.cxx) ot_add_source_file (RobustOptimizationAlgorithm.cxx) ot_add_source_file (SequentialMonteCarloRobustAlgorithm.cxx) +ot_add_source_file (SubsetInverseSamplingResult.cxx) +ot_add_source_file (SubsetInverseSampling.cxx) +ot_add_source_file (InverseFORMResult.cxx) +ot_add_source_file (InverseFORM.cxx) ot_install_header_file (MeasureEvaluation.hxx) ot_install_header_file (MeasureEvaluationImplementation.hxx) @@ -33,7 +37,10 @@ ot_install_header_file (MeasureFactory.hxx) ot_install_header_file (RobustOptimizationProblem.hxx) ot_install_header_file (RobustOptimizationAlgorithm.hxx) ot_install_header_file (SequentialMonteCarloRobustAlgorithm.hxx) - +ot_install_header_file (SubsetInverseSamplingResult.hxx) +ot_install_header_file (SubsetInverseSampling.hxx) +ot_install_header_file (InverseFORMResult.hxx) +ot_install_header_file (InverseFORM.hxx) include_directories ( ${INTERNAL_INCLUDE_DIRS} ) diff --git a/lib/src/InverseFORM.cxx b/lib/src/InverseFORM.cxx new file mode 100755 index 0000000..67c2fd2 --- /dev/null +++ b/lib/src/InverseFORM.cxx @@ -0,0 +1,339 @@ +// -*- C++ -*- +/** + * @brief InverseFORM + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 "otrobopt/InverseFORM.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace OT; + +namespace OTROBOPT { + + +CLASSNAMEINIT(InverseFORM); + +static Factory< InverseFORM > Factory_InverseFORM; + +/* Default constructor for the save/load mechanism */ +InverseFORM::InverseFORM() +: PersistentObject() +, targetBeta_(-Normal().computeQuantile(1e-2)[0]) +{ + // Nothing to do +} + +/* + * @brief Standard constructor: the class is defined by an optimisation algorithm, a failure event and a physical starting point + */ +InverseFORM::InverseFORM(const RandomVector & event, + const String & parameterName, + const Point & physicalStartingPoint) +: PersistentObject() +, event_(event) +, parameterName_(parameterName) +, physicalStartingPoint_(physicalStartingPoint) +, targetBeta_(-Normal().computeQuantile(1e-2)[0]) +{ + const Description parameterDescription(event_.getImplementation()->getFunction().getParameterDescription()); + parameterIndex_ = parameterDescription.find(parameterName); + if (parameterIndex_ >= parameterDescription.getSize()) + throw InvalidArgumentException(HERE) << "Cannot find parameter " << parameterName; +} + +/* Virtual constructor */ +InverseFORM * InverseFORM::clone() const +{ + return new InverseFORM(*this); +} + +/* Result accessor */ +InverseFORMResult InverseFORM::getResult() const +{ + return result_; +} + +/* Result accessor */ +void InverseFORM::setResult(const InverseFORMResult & formResult) +{ + result_ = formResult; +} + + +/* String converter */ +String InverseFORM::__repr__() const +{ + OSS oss; + oss << "class=" << InverseFORM::GetClassName() + << " result=" << result_; + return oss; +} + + +/* Merit function used to optimize the step */ +Scalar InverseFORM::meritFunction(const Point & u, const Scalar levelValue) const +{ + if (fixedStep_) + return -1.0; + + const Scalar unormSquare = u.normSquare(); + return 0.5*unormSquare + (targetBeta_*targetBeta_*unormSquare*1.0 / (variableConvergence_*variableConvergence_) + 1000.0)*fabs(levelValue); +} + + +/* Function that computes the design point by re-using the Analytical::run() and creates a InverseFORMResult */ +void InverseFORM::run() +{ + result_ = InverseFORMResult(); + + const Point u0 = event_.getImplementation()->getAntecedent().getDistribution().getIsoProbabilisticTransformation().operator()(physicalStartingPoint_); + + // sign of target beta + const Scalar signBetaT = targetBeta_ < 0. ? -1.0: 1.0; + + Point u = u0; + Point unew; + Point du; + + Scalar p = event_.getImplementation()->getFunction().getParameter()[parameterIndex_]; + Function G(getG(p)); + Scalar g = G(u0)[0]; + + Scalar g0 = g; + + Scalar gnew = -1.0; + Point dgdu; + Scalar pnew = -1.0; + Scalar dgdp = -1.0; + Scalar dp = -1.0; + const Scalar eps = ResourceMap::GetAsScalar("CenteredFiniteDifferenceGradient-DefaultEpsilon"); + Scalar m = 0.; + Scalar mnew = -1.0; + Scalar beta = signBetaT * u.norm(); + + const Scalar betaC = targetBeta_; + Bool convergence = false; + UnsignedInteger iteration = 0; + while ((!convergence) && (iteration < maximumIteration_) && (!stop_)) + { + // some initializations + event_.getImplementation()->getFunction().getImplementation()->setParameter(Point(1, p)); + m = meritFunction(u, g); + + // gradient / point + dgdu = G.gradient(u) * Point({1.0}); + Scalar dgdunorm = dgdu.norm(); + + // gradient / parameter + Scalar rel = fabs(eps * p); + Scalar delta = rel > SpecFunc::ScalarEpsilon ? rel : eps; + + G = getG(p + delta); + Scalar right = G(u)[0]; + G = getG(p - delta); + Scalar left = G(u)[0]; + + G = getG(p); + dgdp = (right - left) / (2.0 * delta); + + // direction + du = -(betaC / dgdunorm) * dgdu - u; + dp = (fabs(dgdp) < (SpecFunc::ScalarEpsilon * SpecFunc::ScalarEpsilon)) ? 0.0 : (u.dot(dgdu) - g + betaC * dgdunorm) / dgdp; + LOGINFO(OSS() << "InverseFORM::run i=" << iteration << " u="<(distribution.getImplementation().get()); + if (p_composed) + { + ComposedDistribution::DistributionCollection distributionCollection(p_composed->getDistributionCollection()); + for (UnsignedInteger i = 0; i < distributionCollection.getSize(); ++ i) + { + if (distributionCollection[i].getImplementation()->getClassName() == "ConditionalDistribution") + { + DistributionImplementation::PointWithDescriptionCollection parametersCollection(distributionCollection[i].getParametersCollection()); + for (UnsignedInteger j = 0; j < parametersCollection.getSize(); ++ j) + { + const String marginalName(parametersCollection[j].getName()); + const size_t pos = marginalName.find(parameterName_); + if (pos != String::npos) { + const String endName = marginalName.substr(pos, marginalName.size()); + if (endName == parameterName_) + { + parametersCollection[j] = params; + } + } + } + const ConditionalDistribution * p_conditional = dynamic_cast(distributionCollection[i].getImplementation().get()); + if (p_conditional) + { + Distribution conditioning(p_conditional->getConditioningDistribution()); + conditioning.setParametersCollection(parametersCollection); + ConditionalDistribution newConditional(p_conditional->getConditionedDistribution(), conditioning); + distributionCollection[i] = newConditional; + ComposedDistribution newDistribution(distributionCollection); + antecedent = RandomVector(newDistribution); + } // if p_conditional + } // if conditional + } // i + } // if p_composed + } // if composed + + const CompositeRandomVector composite(newFunction, antecedent); + const ThresholdEvent event(composite, op, threshold); + const StandardEvent standardEvent(event); + const Scalar gsign = op(1.0, 2.0) ? 1.0 : -1.0; + const LinearFunction handleOp(Point(1), Point(1, -1.0*gsign*threshold), gsign*IdentityMatrix(1)); + const ComposedFunction composed(handleOp, standardEvent.getImplementation()->getFunction()); + return composed; +} + + +void InverseFORM::setTargetBeta(Scalar targetBeta) { targetBeta_ = targetBeta; } +Scalar InverseFORM::getTargetBeta() const { return targetBeta_; } + +void InverseFORM::setFixedStep(Bool fixedStep) { fixedStep_ = fixedStep; } +Bool InverseFORM::getFixedStep() const { return fixedStep_; } + +void InverseFORM::setFixedStepValue(Scalar fixedStepValue) { fixedStepValue_ = fixedStepValue; } +Scalar InverseFORM::getFixedStepValue() const { return fixedStepValue_; } + +void InverseFORM::setVariableStepMaxIterations(UnsignedInteger variableStepMaxIterations) { variableStepMaxIterations_ = variableStepMaxIterations; } +UnsignedInteger InverseFORM::getVariableStepMaxIterations() const { return variableStepMaxIterations_; } + +void InverseFORM::setMaximumIteration(UnsignedInteger maximumIteration) { maximumIteration_ = maximumIteration; } +UnsignedInteger InverseFORM::getMaximumIteration() const { return maximumIteration_; } + +void InverseFORM::setVariableConvergence(Scalar variableConvergence) { variableConvergence_ = variableConvergence; } +Scalar InverseFORM::getVariableConvergence() const { return variableConvergence_; } + +void InverseFORM::setBetaConvergence(Scalar betaConvergence) { betaConvergence_ = betaConvergence; } +Scalar InverseFORM::getBetaConvergence() const { return betaConvergence_; } + +void InverseFORM::setLimitStateConvergence(Scalar limitStateConvergence) { limitStateConvergence_ = limitStateConvergence; } +Scalar InverseFORM::getLimitStateConvergence() const { return limitStateConvergence_; } + + +/* Method save() stores the object through the StorageManager */ +void InverseFORM::save(Advocate & adv) const +{ + PersistentObject::save(adv); + adv.saveAttribute("event_", event_); + adv.saveAttribute("parameterName_", parameterName_); + adv.saveAttribute("physicalStartingPoint_", physicalStartingPoint_); + adv.saveAttribute("fixedStep_", fixedStep_); + adv.saveAttribute("fixedStepValue_", fixedStepValue_); + adv.saveAttribute("variableStepMaxIterations_", variableStepMaxIterations_); + adv.saveAttribute("maximumIteration_", maximumIteration_); + adv.saveAttribute("variableConvergence_", variableConvergence_); + adv.saveAttribute("betaConvergence_", betaConvergence_); + adv.saveAttribute("limitStateConvergence_", limitStateConvergence_); + adv.saveAttribute("result_", result_); + +} + +/* Method load() reloads the object from the StorageManager */ +void InverseFORM::load(Advocate & adv) +{ + PersistentObject::load(adv); + adv.loadAttribute("event_", event_); + adv.loadAttribute("parameterName_", parameterName_); + adv.loadAttribute("physicalStartingPoint_", physicalStartingPoint_); + adv.loadAttribute("fixedStep_", fixedStep_); + adv.loadAttribute("fixedStepValue_", fixedStepValue_); + adv.loadAttribute("variableStepMaxIterations_", variableStepMaxIterations_); + adv.loadAttribute("maximumIteration_", maximumIteration_); + adv.loadAttribute("variableConvergence_", variableConvergence_); + adv.loadAttribute("betaConvergence_", betaConvergence_); + adv.loadAttribute("limitStateConvergence_", limitStateConvergence_); + adv.loadAttribute("result_", result_); +} + +} diff --git a/lib/src/InverseFORMResult.cxx b/lib/src/InverseFORMResult.cxx new file mode 100755 index 0000000..87693b0 --- /dev/null +++ b/lib/src/InverseFORMResult.cxx @@ -0,0 +1,115 @@ +// -*- C++ -*- +/** + * @brief InverseFORMResult + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 "otrobopt/InverseFORM.hxx" +#include +#include +#include + +using namespace OT; + +namespace OTROBOPT { + +CLASSNAMEINIT(InverseFORMResult); + +static Factory Factory_InverseFORMResult; + +typedef PersistentCollection PersistentSensitivity; + +/* + * @brief Standard constructor: the class is defined by an optimisation algorithm, a failure event and a physical starting point + */ +InverseFORMResult::InverseFORMResult (const Point & standardSpaceDesignPoint, + const RandomVector & event, + const Bool isStandardPointOriginInFailureSpace, + const Point & parameter, + const Description & parameterDescription, + const OT::Point & convergenceCriteriaOptimResult) +: FORMResult(standardSpaceDesignPoint, event, isStandardPointOriginInFailureSpace) +, parameter_(parameter) +, parameterDescription_(parameterDescription) +, convergenceCriteria_(convergenceCriteriaOptimResult) +{ +} + +/* Default constructor */ +InverseFORMResult::InverseFORMResult() +: FORMResult() +{ + // Nothing to do +} + +/* Virtual constructor */ +InverseFORMResult * InverseFORMResult::clone() const +{ + return new InverseFORMResult(*this); +} + + +Description InverseFORMResult::getParameterDescription() const +{ + return parameterDescription_; +} + + +Point InverseFORMResult::getParameter() const +{ + return parameter_; +} + + +OT::Point InverseFORMResult::getConvergenceCriteria() const +{ + return convergenceCriteria_; +} + + +/* String converter */ +String InverseFORMResult::__repr__() const +{ + OSS oss; + oss << "class=" << InverseFORMResult::GetClassName() + << " " << FORMResult::__repr__() + << " parameter_=" << parameter_ + << " parameterDescription_=" << parameterDescription_ + << " convergenceCriteria=" << convergenceCriteria_; + return oss; +} + +/* Method save() stores the object through the StorageManager */ +void InverseFORMResult::save(Advocate & adv) const +{ + FORMResult::save(adv); + adv.saveAttribute( "parameter_", parameter_ ); + adv.saveAttribute( "parameterDescription_", parameterDescription_ ); + adv.saveAttribute( "convergenceCriteria_", convergenceCriteria_); +} + +/* Method load() reloads the object from the StorageManager */ +void InverseFORMResult::load(Advocate & adv) +{ + FORMResult::load(adv); + adv.loadAttribute( "parameter_", parameter_ ); + adv.loadAttribute( "parameterDescription_", parameterDescription_ ); + adv.loadAttribute( "convergenceCriteria_", convergenceCriteria_); +} + +} diff --git a/lib/src/SubsetInverseSampling.cxx b/lib/src/SubsetInverseSampling.cxx new file mode 100644 index 0000000..ba05a09 --- /dev/null +++ b/lib/src/SubsetInverseSampling.cxx @@ -0,0 +1,687 @@ +// -*- C++ -*- +/** + * @brief SubsetInverseSampling + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 "otrobopt/SubsetInverseSampling.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace OT; + +namespace OTROBOPT +{ + +CLASSNAMEINIT(SubsetInverseSampling); + +static Factory Factory_SubsetInverseSampling; + +const UnsignedInteger SubsetInverseSampling::DefaultMaximumOuterSampling = 10000; +const Scalar SubsetInverseSampling::DefaultConditionalProbability = 0.1; +const Scalar SubsetInverseSampling::DefaultProposalRange = 2.0; +const Scalar SubsetInverseSampling::DefaultBetaMin = 2.0; + + +/* Default constructor */ +SubsetInverseSampling::SubsetInverseSampling() +: EventSimulation() +, proposalRange_(0.) +, conditionalProbability_(0.) +, iSubset_(false) +, betaMin_(0.) +, keepEventSample_(false) +, targetProbability_(0.) +, numberOfSteps_(0) +{ +} + + +/* Constructor with parameters */ +SubsetInverseSampling::SubsetInverseSampling(const RandomVector & event, + const Scalar targetProbability, + const Scalar proposalRange, + const Scalar conditionalProbability) +: EventSimulation(event) +, proposalRange_(proposalRange) +, conditionalProbability_(conditionalProbability) +, iSubset_(false) +, betaMin_(DefaultBetaMin) +, keepEventSample_(false) +, targetProbability_(targetProbability) +, numberOfSteps_(0) +{ + setMaximumOuterSampling( DefaultMaximumOuterSampling );// overide simulation default outersampling + UnsignedInteger outputDimension = event.getFunction().getOutputDimension(); + if ( outputDimension > 1 ) + throw InvalidArgumentException(HERE) << "Output dimension for SubsetInverseSampling cannot be greater than 1, here output dimension=" << outputDimension; +} + + +/* Virtual constructor */ +SubsetInverseSampling * SubsetInverseSampling::clone() const +{ + return new SubsetInverseSampling(*this); +} + + +/* Performs the actual computation. */ +void SubsetInverseSampling::run() +{ + // First, initialize some parameters + convergenceStrategy_.clear(); + numberOfSteps_ = 0; + thresholdPerStep_.clear(); + thresholdCoefficientOfVariationPerStep_.clear(); + gammaPerStep_.clear(); + coefficientOfVariationPerStep_.clear(); + probabilityEstimatePerStep_.clear(); + eventInputSample_.clear(); + eventOutputSample_.clear(); + allPointSample_.clear(); + allLevelSample_.clear(); + + dimension_ = getEvent().getAntecedent().getDimension(); + + if ( getMaximumCoefficientOfVariation() != ResourceMap::GetAsScalar( "SimulationAlgorithm-DefaultMaximumCoefficientOfVariation" ) ) + Log::Warn(OSS() << "The maximum coefficient of variation was set. It won't be used as termination criteria."); + + if ( conditionalProbability_ * getMaximumOuterSampling() * getBlockSize() < 1 ) + throw InvalidArgumentException(HERE) << "The number of samples (maximumOuterSampling x blockSize : " << getMaximumOuterSampling() << " x " << getBlockSize() << ") should be >= " << ceil( 1.0 / conditionalProbability_ ); + + if ( getMaximumOuterSampling() * getBlockSize() <= 100 ) + Log::Warn(OSS() << "The number of samples per step is very low : " << getMaximumOuterSampling()*getBlockSize() << "."); + + // perform isoprobabilistic transformation (the study is done in the standard space): + standardEvent_ = StandardEvent(getEvent()); + + Scalar currentCoVsquare = 0.0; + Scalar varianceEstimate = 0.0; + Scalar coefficientOfVariationSquare = 0.0; + Scalar tempCoefficientOfVariation = 0.0; + Scalar thresholdVariance = 0.0; + Scalar thresholdCoefficientOfVariationSquare = 0.0; + Scalar gamma = 0.0; + + // allocate input/output samples + const UnsignedInteger maximumOuterSampling = getMaximumOuterSampling(); + const UnsignedInteger blockSize = getBlockSize(); + const UnsignedInteger N = maximumOuterSampling * blockSize; + currentPointSample_ = Sample( N, dimension_ ); + currentLevelSample_ = Sample( N, getEvent().getFunction().getOutputDimension() ); + + // Step 1: sampling + for ( UnsignedInteger i = 0; i < maximumOuterSampling; ++ i ) + { + Sample inputSample; + if (!iSubset_) { + // crude MC + inputSample = standardEvent_.getAntecedent().getDistribution().getSample( blockSize ); + } + else { + // conditional sampling + TruncatedDistribution truncatedChiSquare(ChiSquare(dimension_), betaMin_ * betaMin_, TruncatedDistribution::LOWER); + Normal normal(dimension_); + inputSample = Sample(0, dimension_); + for ( UnsignedInteger j = 0; j < blockSize; ++ j) + { + Point direction = normal.getRealization(); + Scalar norm = direction.norm(); + Scalar radius = sqrt(truncatedChiSquare.getRealization()[0]); + if (fabs(norm) > 1e-12) + { + radius *= 1.0 / norm; + } + inputSample.add(direction * radius); + } + } + Sample blockSample( standardEvent_.getFunction()( inputSample ) ); + for ( UnsignedInteger j = 0 ; j < blockSize; ++ j ) + { + currentPointSample_[ i*blockSize+j ] = inputSample[j]; + currentLevelSample_[ i*blockSize+j ] = blockSample[j]; + } + } + ++ numberOfSteps_; + + // Stop if the wanted probability is greater than the target probability per step + Bool stop = targetProbability_ >= conditionalProbability_; + + if (stop) + { + setConditionalProbability(targetProbability_); + } + + // computation of the first intermediate threshold with the sample create with a normal distribution */ + Scalar currentThreshold = computeThreshold(); + + // compute monte carlo probability estimate + Scalar probabilityEstimate = computeProbability(1.0, currentThreshold); + + if (iSubset_) + { + Scalar correction = 1.0 - ChiSquare(standardEvent_.getImplementation()->getAntecedent().getDistribution().getDimension()).computeCDF(betaMin_ * betaMin_); + probabilityEstimate *= correction; + } + + // compute the coefficient of variation + if (probabilityEstimate > 0.0) + { + // ... compute coefficient of variation + coefficientOfVariationSquare = (1.0 - probabilityEstimate) / (probabilityEstimate * currentLevelSample_.getSize() * 1.0); + // ... compute variance estimate + varianceEstimate = coefficientOfVariationSquare * probabilityEstimate * probabilityEstimate; + + // compute coefficient of variation of the threshold + Distribution ks = KernelSmoothing().build(currentLevelSample_); + thresholdVariance = probabilityEstimate * (1. - probabilityEstimate) / currentLevelSample_.getSize() / pow(ks.computePDF(currentThreshold), 2); + thresholdCoefficientOfVariationSquare = thresholdVariance / (currentThreshold * currentThreshold); + } + + allPointSample_.add(currentPointSample_); + allLevelSample_.add(currentLevelSample_); + thresholdPerStep_.add( currentThreshold ); + gammaPerStep_.add(0.); + probabilityEstimatePerStep_.add(probabilityEstimate); + coefficientOfVariationPerStep_.add(sqrt(coefficientOfVariationSquare)); + thresholdCoefficientOfVariationPerStep_.add(sqrt(thresholdCoefficientOfVariationSquare)); + + // as long as the conditional failure domain do not overlap the global one + while ( !stop ) + { + // prepare new seeds + initializeSeed( currentThreshold ); + + // create new points using MCMC + generatePoints( currentThreshold ); + + // save the level sample + allPointSample_.add(currentPointSample_); + allLevelSample_.add(currentLevelSample_); + + // compute new threshold + currentThreshold = computeThreshold(); + + // compute probability estimate on the current sample and group seeds at the beginning of the work sample + Scalar currentProbabilityEstimate = computeProbability( probabilityEstimate, currentThreshold ); + + // update probability estimate + probabilityEstimate *= currentProbabilityEstimate; + + // update coefficient of variation + gamma = computeVarianceGamma( currentProbabilityEstimate, currentThreshold ); + currentCoVsquare = (1.0 - currentProbabilityEstimate) / (currentProbabilityEstimate * currentLevelSample_.getSize() * 1.0); + tempCoefficientOfVariation = sqrt(coefficientOfVariationSquare + (1.0 + gamma) * currentCoVsquare); + + // update stopping criterion : target proba > inf bound of confidence interval of pf + stop = targetProbability_ >= probabilityEstimate * (1 - tempCoefficientOfVariation * DistFunc::qNormal(0.95)) ; + + if (stop) + { + // change the target probability of the final step + setConditionalProbability(targetProbability_ / probabilityEstimatePerStep_[numberOfSteps_-1]); + // compute the final threshold + currentThreshold = computeThreshold(); + // compute the current probability estimate + Scalar currentProbabilityEstimate = computeProbability(probabilityEstimatePerStep_[numberOfSteps_-1], currentThreshold); + // update probability estimate + probabilityEstimate = probabilityEstimatePerStep_[numberOfSteps_-1] * currentProbabilityEstimate; + + // update a second time the conditional probability to be as close as possible of the target probability + // change the target probability of the final step + setConditionalProbability(targetProbability_ * targetProbability_ / (probabilityEstimatePerStep_[numberOfSteps_-1] * probabilityEstimate)); + // compute the final threshold + currentThreshold = computeThreshold(); + // compute the current probability estimate + currentProbabilityEstimate = computeProbability( probabilityEstimatePerStep_[numberOfSteps_-1], currentThreshold ); + // update probability estimate + probabilityEstimate = probabilityEstimatePerStep_[numberOfSteps_-1] * currentProbabilityEstimate; + } + + // update coefficient of variation + gamma = computeVarianceGamma(currentProbabilityEstimate, currentThreshold); + currentCoVsquare = (1.0 - currentProbabilityEstimate) / (currentProbabilityEstimate * currentLevelSample_.getSize() * 1.0); + coefficientOfVariationSquare += (1.0 + gamma) * currentCoVsquare; + // update threshold coefficient of variation + Distribution ks = KernelSmoothing().build(currentLevelSample_); + thresholdVariance = currentProbabilityEstimate * (1. - currentProbabilityEstimate) / currentLevelSample_.getSize() / pow(ks.computePDF(currentThreshold), 2); + thresholdCoefficientOfVariationSquare += (1.0 + gamma) * thresholdVariance / (currentThreshold * currentThreshold); + + thresholdPerStep_.add( currentThreshold ); + gammaPerStep_.add(gamma); + probabilityEstimatePerStep_.add(probabilityEstimate); + coefficientOfVariationPerStep_.add(sqrt(coefficientOfVariationSquare)); + thresholdCoefficientOfVariationPerStep_.add(sqrt(thresholdCoefficientOfVariationSquare)); + + // stop if the number of subset steps is too high, else results are not numerically defined anymore + if ( fabs( pow( probabilityEstimate, 2.) ) < SpecFunc::MinScalar ) + throw NotDefinedException(HERE) << "Probability estimate too small: " << probabilityEstimate; + + // compute variance estimate + varianceEstimate = coefficientOfVariationSquare * pow( probabilityEstimate, 2. ); + thresholdVariance = thresholdCoefficientOfVariationSquare * pow( currentThreshold, 2. ); + + ++ numberOfSteps_; + + if (stopCallback_.first && stopCallback_.first(stopCallback_.second)) + throw InternalException(HERE) << "User stopped simulation"; + } + + // define the threshold distribution + thresholdDistribution_ = Normal(currentThreshold, sqrt(thresholdVariance)); + + //update the event with the final threshold + RandomVector modified_event = ThresholdEvent(CompositeRandomVector(getEvent().getFunction(), getEvent().getAntecedent()), getEvent().getOperator(), currentThreshold); + + setResult( SubsetInverseSamplingResult(modified_event, probabilityEstimate, varianceEstimate, numberOfSteps_ * getMaximumOuterSampling(), getBlockSize(), sqrt( coefficientOfVariationSquare ), currentThreshold) ); + + // keep the event sample if requested + if (keepEventSample_) + { + eventInputSample_ = Sample(0, dimension_); + eventOutputSample_ = Sample (0, getEvent().getFunction().getOutputDimension()); + for (UnsignedInteger i = 0; i < currentPointSample_.getSize(); ++ i) + { + if (getEvent().getOperator()(currentLevelSample_(i, 0), currentThreshold)) + { + eventInputSample_.add(getEvent().getAntecedent().getDistribution().getInverseIsoProbabilisticTransformation()(currentPointSample_[i])); + eventOutputSample_.add(currentLevelSample_[i]); + } + } + } + + // free work samples + currentLevelSample_.clear(); + currentPointSample_.clear(); +} + + +/* Compute the block sample */ +Sample SubsetInverseSampling::computeBlockSample() +{ + return Sample(); +} + + +/* Compute the new threshold corresponding to the target failure probability */ +Scalar SubsetInverseSampling::computeThreshold() +{ + // compute the quantile according to the event operator + const Scalar ratio = getEvent().getOperator()(1.0, 2.0) ? conditionalProbability_ : 1.0 - conditionalProbability_; + + const Scalar currentThreshold = currentLevelSample_.computeQuantile( ratio )[0]; + + return currentThreshold; +} + + +Scalar SubsetInverseSampling::computeProbability(Scalar probabilityEstimateFactor, Scalar threshold) +{ + const UnsignedInteger maximumOuterSampling = getMaximumOuterSampling(); + const UnsignedInteger blockSize = getBlockSize(); + Scalar probabilityEstimate = 0.0; + Scalar varianceEstimate = 0.0; + + for (UnsignedInteger i = 0; i < maximumOuterSampling; ++ i) + { + const Scalar size = i + 1.0; + Scalar meanBlock = 0.0; + Scalar varianceBlock = 0.0; + for (UnsignedInteger j = 0 ; j < blockSize; ++ j) + { + if (getEvent().getOperator()( currentLevelSample_(i*blockSize+j, 0), threshold)) + { + // update local mean and variance + meanBlock += 1.0 / blockSize; + varianceBlock += 1.0 * 1.0 / blockSize; + } + } + varianceBlock -= pow( meanBlock, 2.0 ); + + // update global mean and variance + varianceEstimate = (varianceBlock + (size - 1.0) * varianceEstimate) / size + (1.0 - 1.0 / size) * (probabilityEstimate - meanBlock) * (probabilityEstimate - meanBlock) / size; + probabilityEstimate = std::min(1.0, (meanBlock + (size - 1.0) * probabilityEstimate) / size); + + // store convergence at each block + Point convergencePoint(2); + convergencePoint[0] = probabilityEstimate * probabilityEstimateFactor; + convergencePoint[1] = varianceEstimate * probabilityEstimateFactor * probabilityEstimateFactor / size; + convergenceStrategy_.store(convergencePoint); + } + + // cannot determine next subset domain if no variance + const Scalar epsilon = ResourceMap::GetAsScalar( "SpecFunc-Precision" ); + if ( fabs( varianceEstimate ) < epsilon ) + throw NotDefinedException(HERE) << "Null output variance"; + + return probabilityEstimate; +} + + +/* Sort new seeds */ +void SubsetInverseSampling::initializeSeed(Scalar threshold) +{ + UnsignedInteger seedIndex = 0; + const UnsignedInteger maximumOuterSampling = getMaximumOuterSampling(); + const UnsignedInteger blockSize = getBlockSize(); + for ( UnsignedInteger i = 0; i < maximumOuterSampling; ++ i ) + { + for ( UnsignedInteger j = 0 ; j < blockSize; ++ j ) + { + if ( getEvent().getOperator()( currentLevelSample_[ i*blockSize+j ][0], threshold ) ) + { + // initialize seeds : they're grouped at the beginning of the sample + currentPointSample_[ seedIndex ] = currentPointSample_[ i*blockSize+j ]; + currentLevelSample_[ seedIndex ] = currentLevelSample_[ i*blockSize+j ]; + ++ seedIndex; + } + } + } +} + + +/* Compute the correlation on markov chains at the current state of the algorithm */ +Scalar SubsetInverseSampling::computeVarianceGamma(Scalar currentFailureProbability, Scalar threshold) +{ + const UnsignedInteger N = currentPointSample_.getSize(); + const UnsignedInteger Nc = std::max(1, conditionalProbability_ * N); + Matrix IndicatorMatrice( Nc, N / Nc ); + Point correlationSequence( N / Nc - 1 ); + Scalar currentFailureProbability2 = pow( currentFailureProbability, 2. ); + for ( UnsignedInteger i = 0; i < N / Nc; ++ i ) + { + for ( UnsignedInteger j = 0; j < Nc; ++ j ) + { + IndicatorMatrice(j, i) = getEvent().getOperator()(currentLevelSample_[ i*Nc+j ][0], threshold); + } + } + for ( UnsignedInteger k = 0; k < N / Nc - 1; ++ k ) + { + for ( UnsignedInteger j = 0; j < Nc; ++ j ) + { + for ( UnsignedInteger l = 0; l < N / Nc - k - 1; ++ l ) + { + correlationSequence[k] += 1.0 * IndicatorMatrice(j, l) * IndicatorMatrice(j, l + (k + 1)); + } + } + correlationSequence[k] /= 1.0 * N - 1.0 * (k + 1) * Nc; + correlationSequence[k] -= currentFailureProbability2; + } + const Scalar R0 = currentFailureProbability * ( 1.0 - currentFailureProbability ); + Point rho = ((1.0 / R0) * correlationSequence); + Scalar gamma = 0.0; + for ( UnsignedInteger k = 0; k < N / Nc - 1; ++ k ) + { + gamma += 2.0 * (1.0 - (k + 1) * 1.0 * Nc / N) * rho[k]; + } + return gamma; +} + + +/* Iterate one step of the algorithm */ +void SubsetInverseSampling::generatePoints(Scalar threshold) +{ + UnsignedInteger maximumOuterSampling = getMaximumOuterSampling(); + UnsignedInteger blockSize = getBlockSize(); + Distribution randomWalk(ComposedDistribution(ComposedDistribution::DistributionCollection(dimension_, Uniform(-0.5*proposalRange_, 0.5*proposalRange_)))); + UnsignedInteger N = currentPointSample_.getSize(); // total sample size + UnsignedInteger Nc = conditionalProbability_ * N; //number of seeds (also = maximumOuterSampling*blockSize) + + for ( UnsignedInteger i = 0; i < maximumOuterSampling; ++ i ) + { + Sample inputSample( blockSize, dimension_ ); + for ( UnsignedInteger j = 0; j < blockSize; ++ j ) + { + // assign the new point to the seed, seed points being regrouped at the beginning of the sample + if ( i*blockSize+j >= Nc ) + { + currentPointSample_[ i*blockSize+j ] = currentPointSample_[ i*blockSize+j-Nc ]; + currentLevelSample_[ i*blockSize+j ] = currentLevelSample_[ i*blockSize+j-Nc ]; + } + + // generate a new point + Point oldPoint( currentPointSample_[ i*blockSize+j ] ); + Point newPoint( oldPoint + randomWalk.getRealization() ); + + // 1. accept / reject new components + Point uniform( RandomGenerator::Generate(dimension_) ); + for (UnsignedInteger k = 0; k < dimension_; ++ k) + { + // compute ratio + const Scalar ratio = std::min(1.0, exp(0.5 * (oldPoint[k] * oldPoint[k] - newPoint[k] * newPoint[k]))); + + // accept new point with probability ratio + if (ratio < uniform[k]) + { + newPoint[k] = oldPoint[k]; + } + } + + inputSample[j] = newPoint; + } + + Sample blockSample( standardEvent_.getFunction()( inputSample ) ); + + for ( UnsignedInteger j = 0; j < getBlockSize(); ++ j ) + { + // 2. accept the new point if in the failure domain + if ( getEvent().getOperator()( blockSample[j][0], threshold ) ) + { + currentPointSample_[ i*blockSize+j ] = inputSample[j]; + currentLevelSample_[ i*blockSize+j ] = blockSample[j]; + } + } + } +} + + +/* Confidence Length of the threshold */ +Scalar SubsetInverseSampling::getThresholdConfidenceLength(const Scalar level) const +{ + Scalar thresholdInf = thresholdDistribution_.computeQuantile((1 - level)/2)[0]; + Scalar thresholdSup = thresholdDistribution_.computeQuantile(1.-(1 - level)/2)[0]; + return Scalar (std::max(thresholdSup, thresholdInf) - std::min(thresholdSup, thresholdInf)); +} + + +/* Markov parameter accessor */ +void SubsetInverseSampling::setProposalRange(Scalar proposalRange) +{ + proposalRange_ = proposalRange; +} + + +Scalar SubsetInverseSampling::getProposalRange() const +{ + return proposalRange_; +} + + +/* Ratio accessor */ +void SubsetInverseSampling::setConditionalProbability(Scalar conditionalProbability) +{ + if ( (conditionalProbability <= 0.) || (conditionalProbability >= 1.) ) throw InvalidArgumentException(HERE) << "In setConditionalProbability::Probability should be in (0, 1)"; + conditionalProbability_ = conditionalProbability; +} + +Scalar SubsetInverseSampling::getConditionalProbability() const +{ + return conditionalProbability_; +} + +/* target probability accessor */ +void SubsetInverseSampling::setTargetProbability(Scalar targetProbability) +{ + if ( (targetProbability <= 0.) || (targetProbability >= 1.) ) throw InvalidArgumentException(HERE) << "In setTargetProbability::Probability should be in (0, 1)"; + targetProbability_ = targetProbability; +} + +Scalar SubsetInverseSampling::getTargetProbability() const +{ + return targetProbability_; +} + +UnsignedInteger SubsetInverseSampling::getNumberOfSteps() +{ + return numberOfSteps_; +} + + +OT::Point SubsetInverseSampling::getGammaPerStep() const +{ + return gammaPerStep_; +} + + +OT::Point SubsetInverseSampling::getCoefficientOfVariationPerStep() const +{ + return coefficientOfVariationPerStep_; +} + + +OT::Point SubsetInverseSampling::getProbabilityEstimatePerStep() const +{ + return probabilityEstimatePerStep_; +} + + +OT::Point SubsetInverseSampling::getThresholdPerStep() const +{ + return thresholdPerStep_; +} + +OT::Point SubsetInverseSampling::getThresholdCoefficientOfVariationPerStep() const +{ + return thresholdCoefficientOfVariationPerStep_; +} + +void SubsetInverseSampling::setKeepEventSample(bool keepEventSample) +{ + keepEventSample_ = keepEventSample; +} + + +Sample SubsetInverseSampling::getEventInputSample() const +{ + return eventInputSample_; +} + + +Sample SubsetInverseSampling::getEventOutputSample() const +{ + return eventOutputSample_; +} + +SubsetInverseSampling::SampleCollection SubsetInverseSampling::getInputSample() const +{ + return allPointSample_; +} + +SubsetInverseSampling::SampleCollection SubsetInverseSampling::getOutputSample() const +{ + return allLevelSample_; +} + +void SubsetInverseSampling::setISubset(OT::Bool iSubset) +{ + iSubset_ = iSubset; +} + +void SubsetInverseSampling::setBetaMin(Scalar betaMin) +{ + if (betaMin <= 0.) throw InvalidArgumentException(HERE) << "Beta min should be positive"; + betaMin_ = betaMin; +} + +void SubsetInverseSampling::setResult(const SubsetInverseSamplingResult & result) +{ + result_ = result; +} + +SubsetInverseSamplingResult SubsetInverseSampling::getResult() const +{ + return result_; +} + +/* String converter */ +String SubsetInverseSampling::__repr__() const +{ + OSS oss; + oss << "class=" << getClassName() + << " derived from " << EventSimulation::__repr__() + << " event=" << event_ + << " targetProbability=" << targetProbability_ + << " proposalRange=" << proposalRange_ + << " conditionalProbability=" << conditionalProbability_ + << " keepEventSample_=" << keepEventSample_; + return oss; +} + + +/* Method save() stores the object through the StorageManager */ +void SubsetInverseSampling::save(Advocate & adv) const +{ + EventSimulation::save(adv); + adv.saveAttribute("targetProbability", targetProbability_); + adv.saveAttribute("proposalRange_", proposalRange_); + adv.saveAttribute("conditionalProbability_", conditionalProbability_); + adv.saveAttribute("iSubset_", iSubset_); + adv.saveAttribute("betaMin_", betaMin_); + adv.saveAttribute("keepEventSample_", keepEventSample_); + + adv.saveAttribute("numberOfSteps_", numberOfSteps_); + adv.saveAttribute("thresholdPerStep_", thresholdPerStep_); + adv.saveAttribute("gammaPerStep_", gammaPerStep_); + adv.saveAttribute("coefficientOfVariationPerStep_", coefficientOfVariationPerStep_); + adv.saveAttribute("probabilityEstimatePerStep_", probabilityEstimatePerStep_); +} + + +/* Method load() reloads the object from the StorageManager */ +void SubsetInverseSampling::load(Advocate & adv) +{ + EventSimulation::load(adv); + adv.loadAttribute("targetProbability", targetProbability_); + adv.loadAttribute("proposalRange_", proposalRange_); + adv.loadAttribute("conditionalProbability_", conditionalProbability_); + adv.loadAttribute("keepEventSample_", keepEventSample_); + adv.loadAttribute("iSubset_", iSubset_); + adv.loadAttribute("betaMin_", betaMin_); + + adv.loadAttribute("numberOfSteps_", numberOfSteps_); + adv.loadAttribute("thresholdPerStep_", thresholdPerStep_); + adv.loadAttribute("gammaPerStep_", gammaPerStep_); + adv.loadAttribute("coefficientOfVariationPerStep_", coefficientOfVariationPerStep_); + adv.loadAttribute("probabilityEstimatePerStep_", probabilityEstimatePerStep_); +} + + + + +} /* namespace OTROBOPT */ diff --git a/lib/src/SubsetInverseSamplingResult.cxx b/lib/src/SubsetInverseSamplingResult.cxx new file mode 100644 index 0000000..edc0641 --- /dev/null +++ b/lib/src/SubsetInverseSamplingResult.cxx @@ -0,0 +1,99 @@ +// -*- C++ -*- +/** + * @brief SubsetSamplingImplementation + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 "otrobopt/SubsetInverseSamplingResult.hxx" + +using namespace OT; + +namespace OTROBOPT +{ + +/* + * @class SubsetInverseSamplingResult + */ + +CLASSNAMEINIT(SubsetInverseSamplingResult); + +static Factory Factory_SubsetInverseSamplingResult; + + +/* Default constructor */ +SubsetInverseSamplingResult::SubsetInverseSamplingResult() +: ProbabilitySimulationResult() +{ +} + + +/* Constructor with parameters */ +SubsetInverseSamplingResult::SubsetInverseSamplingResult(const RandomVector & event, + const Scalar probabilityEstimate, + const Scalar varianceEstimate, + const UnsignedInteger outerSampling, + const UnsignedInteger blockSize, + const Scalar coefficientOfVariation, + const Scalar threshold) +: ProbabilitySimulationResult(event, probabilityEstimate, varianceEstimate, outerSampling, blockSize), + coefficientOfVariation_(coefficientOfVariation), + threshold_(threshold) +{ +} + + +/* Virtual constructor */ +SubsetInverseSamplingResult * SubsetInverseSamplingResult::clone() const +{ + return new SubsetInverseSamplingResult(*this); +} + + +/* Coefficient of variation estimate accessor */ +Scalar SubsetInverseSamplingResult::getCoefficientOfVariation() const +{ + return coefficientOfVariation_; +} + +/* String converter */ +String SubsetInverseSamplingResult::__repr__() const +{ + OSS oss; + oss << ProbabilitySimulationResult::__repr__(); + oss << " threshold=" << threshold_; + return oss; +} + + +/* Method save() stores the object through the StorageManager */ +void SubsetInverseSamplingResult::save(Advocate & adv) const +{ + ProbabilitySimulationResult::save(adv); + adv.saveAttribute("coefficientOfVariation_", coefficientOfVariation_); +} + + +/* Method load() reloads the object from the StorageManager */ +void SubsetInverseSamplingResult::load(Advocate & adv) +{ + ProbabilitySimulationResult::load(adv); + adv.loadAttribute("coefficientOfVariation_", coefficientOfVariation_); +} + + +} /* namespace OTROBOPT */ diff --git a/lib/src/otrobopt/InverseFORM.hxx b/lib/src/otrobopt/InverseFORM.hxx new file mode 100755 index 0000000..16a6b20 --- /dev/null +++ b/lib/src/otrobopt/InverseFORM.hxx @@ -0,0 +1,122 @@ +// -*- C++ -*- +/** + * @brief InverseFORM + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 OTROBOPT_INVERSEFORM_HXX +#define OTROBOPT_INVERSEFORM_HXX + +#include "otrobopt/InverseFORMResult.hxx" + +#include +#include + +namespace OTROBOPT { + +/** + * @class InverseFORM + * InverseFORM implements the First Order Reliability Method + * and compute the results of such kind of analyses + */ + +class OTROBOPT_API InverseFORM + : public OT::PersistentObject +{ + + CLASSNAME; +public: + + /** Default constructor */ + InverseFORM(); + + /** Constructor with parameters */ + InverseFORM(const OT::RandomVector & event, + const OT::String & parameterName, + const OT::Point & physicalStartingPoint); + + /** Virtual constructor */ + virtual InverseFORM * clone() const override; + + void setTargetBeta(OT::Scalar targetBeta); + OT::Scalar getTargetBeta() const; + + void setFixedStep(OT::Bool fixedStep); + OT::Bool getFixedStep() const; + + void setFixedStepValue(OT::Scalar fixedStepValue); + OT::Scalar getFixedStepValue() const; + + void setVariableStepMaxIterations(OT::UnsignedInteger variableStepMaxIterations); + OT::UnsignedInteger getVariableStepMaxIterations() const; + + void setMaximumIteration(OT::UnsignedInteger maximumIteration); + OT::UnsignedInteger getMaximumIteration() const; + + void setVariableConvergence(OT::Scalar variableConvergence); + OT::Scalar getVariableConvergence() const; + + void setBetaConvergence(OT::Scalar betaConvergence); + OT::Scalar getBetaConvergence() const; + + void setLimitStateConvergence(OT::Scalar limitStateConvergence); + OT::Scalar getLimitStateConvergence() const; + + /** Result accessor */ + void setResult(const InverseFORMResult & result); + InverseFORMResult getResult() const; + + /** String converter */ + OT::String __repr__() const override; + + /** Function that computes the design point by re-using the Analytical::run() and creates a InverseFORM::Result */ + void run(); + + /** Method save() stores the object through the StorageManager */ + void save(OT::Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(OT::Advocate & adv) override; + +private: + OT::Scalar meritFunction(const OT::Point & u, const OT::Scalar levelValue) const; + OT::Function getG(const OT::Scalar p); + + + OT::RandomVector event_; + OT::String parameterName_; + OT::UnsignedInteger parameterIndex_ = 0; + OT::Point physicalStartingPoint_; + + InverseFORMResult result_; + OT::Scalar targetBeta_ = 0.0; + OT::Bool fixedStep_ = true; + OT::Scalar fixedStepValue_ = 1.0; + OT::UnsignedInteger variableStepMaxIterations_ = 5; + OT::UnsignedInteger maximumIteration_ = 100; + OT::Scalar variableConvergence_ = 1e-3; + OT::Scalar betaConvergence_ = 1e-2; + OT::Scalar limitStateConvergence_ = 1e-3; + OT::Bool stop_ = false; +}; + +} + +#endif + + + diff --git a/lib/src/otrobopt/InverseFORMResult.hxx b/lib/src/otrobopt/InverseFORMResult.hxx new file mode 100755 index 0000000..364c186 --- /dev/null +++ b/lib/src/otrobopt/InverseFORMResult.hxx @@ -0,0 +1,79 @@ +// -*- C++ -*- +/** + * @brief InverseFORMResult + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 OTROBOPT_INVERSEFORMRESULT_HXX +#define OTROBOPT_INVERSEFORMRESULT_HXX + +#include +#include "otrobopt/OTRobOptprivate.hxx" + +namespace OTROBOPT { + + +/** + * @class InverseFORMResult + */ +class OTROBOPT_API InverseFORMResult: + public OT::FORMResult +{ + CLASSNAME; +public: + + /** Standard constructor */ + InverseFORMResult(const OT::Point & standardSpaceDesignPoint, + const OT::RandomVector & event, + const OT::Bool isStandardPointOriginInFailureSpace, + const OT::Point & parameter, + const OT::Description & parameterDescription, + const OT::Point & convergenceCriteria); + + /* Default constructor */ + InverseFORMResult(); + + /** Virtual constructor */ + virtual InverseFORMResult * clone() const override; + + /** String converter */ + OT::String __repr__() const override; + + /** Method save() stores the object through the StorageManager */ + void save(OT::Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(OT::Advocate & adv) override; + + // parameter accessor + OT::Point getParameter() const; + OT::Description getParameterDescription() const; + + OT::Point getConvergenceCriteria() const; + +protected: + OT::Point parameter_; + OT::Description parameterDescription_; + OT::Point convergenceCriteria_; +}; + +} + +#endif + + + diff --git a/lib/src/otrobopt/SubsetInverseSampling.hxx b/lib/src/otrobopt/SubsetInverseSampling.hxx new file mode 100644 index 0000000..734035f --- /dev/null +++ b/lib/src/otrobopt/SubsetInverseSampling.hxx @@ -0,0 +1,169 @@ +// -*- C++ -*- +/** + * @brief SubsetInverseSampling + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 OTROBOPT_SUBSETINVERSESAMPLING_HXX +#define OTROBOPT_SUBSETINVERSESAMPLING_HXX + +#include +#include +#include "otrobopt/SubsetInverseSamplingResult.hxx" + +namespace OTROBOPT +{ + + +class OTROBOPT_API SubsetInverseSampling +: public OT::EventSimulation +{ +CLASSNAME +public: + + typedef OT::Collection SampleCollection; + typedef OT::PersistentCollection SamplePersistentCollection; + + /** Default Parameters */ + static const OT::UnsignedInteger DefaultMaximumOuterSampling; + static const OT::Scalar DefaultProposalRange; + static const OT::Scalar DefaultConditionalProbability; + static const OT::Scalar DefaultBetaMin; + + /** Default Constructor */ + SubsetInverseSampling(); + + /** Constructor with parameters */ + SubsetInverseSampling(const OT::RandomVector & event, + const OT::Scalar targetProbability, + const OT::Scalar proposalRange = DefaultProposalRange, + const OT::Scalar conditionalProbability = DefaultConditionalProbability); + + /** Virtual constructor */ + SubsetInverseSampling * clone() const override; + + /** Result accessor */ + SubsetInverseSamplingResult getResult() const; + + /** The range of the uniform proposal pdf */ + void setProposalRange(OT::Scalar proposalRange); + OT::Scalar getProposalRange() const; + + /** Ratio parameter */ + void setConditionalProbability(OT::Scalar conditionalProbability); + OT::Scalar getConditionalProbability() const; + + /** final target probability */ + void setTargetProbability(OT::Scalar targetProbability); + OT::Scalar getTargetProbability() const; + + /** Accessor to the achieved number of steps */ + OT::UnsignedInteger getNumberOfSteps(); + + OT::Scalar getThresholdConfidenceLength(const OT::Scalar level = OT::ResourceMap::GetAsScalar("ProbabilitySimulationResult-DefaultConfidenceLevel")) const; + + /** Stepwise result accessors */ + OT::Point getThresholdPerStep() const; + OT::Point getGammaPerStep() const; + OT::Point getCoefficientOfVariationPerStep() const; + OT::Point getProbabilityEstimatePerStep() const; + OT::Point getThresholdCoefficientOfVariationPerStep() const; + + /** Keep event sample */ + void setKeepEventSample(bool keepEventSample); + + /** Event input/output sample accessor */ + OT::Sample getEventInputSample() const; + OT::Sample getEventOutputSample() const; + + /** All level sample accessor*/ + SampleCollection getOutputSample() const; + SampleCollection getInputSample() const; + + /** i-subset */ + void setISubset(OT::Bool iSubset); + void setBetaMin(OT::Scalar betaMin); + + /** Performs the actual computation. */ + void run() override; + + /** String converter */ + OT::String __repr__() const override; + + /** Method save() stores the object through the StorageManager */ + void save(OT::Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(OT::Advocate & adv) override; + +private: + /** Result accessor */ + void setResult(const SubsetInverseSamplingResult & result); + + /** Compute the block sample */ + OT::Sample computeBlockSample() override; + + /** Compute the new threshold corresponding to the conditional failure probability */ + OT::Scalar computeThreshold(); + + /** compute probability estimate on the current sample */ + OT::Scalar computeProbability(OT::Scalar probabilityEstimate, OT::Scalar threshold); + + /** Sort new seeds */ + void initializeSeed(OT::Scalar threshold); + + /** Compute the correlation on markov chains at the current state of the algorithm */ + OT::Scalar computeVarianceGamma(OT::Scalar currentFailureProbability, OT::Scalar threshold); + + /** Generate new points in the conditional failure domain */ + void generatePoints(OT::Scalar threshold); + + // The result + SubsetInverseSamplingResult result_; + + // some parameters + OT::Scalar proposalRange_;// width of the proposal pdf + OT::Scalar conditionalProbability_;// conditional probability at each subset + OT::Bool iSubset_;// conditional pre-sampling + OT::Scalar betaMin_;// pre-sampling hypersphere exclusion radius + OT::Bool keepEventSample_;// do we keep the event sample ? + OT::Scalar targetProbability_;// final target probability + + // some results + OT::UnsignedInteger numberOfSteps_;// number of subset steps + OT::Point thresholdPerStep_;// intermediate thresholds + OT::Point gammaPerStep_;// intermediate gammas + OT::Point coefficientOfVariationPerStep_;// intermediate COVS + OT::Point probabilityEstimatePerStep_;// intermediate PFs + OT::Point thresholdCoefficientOfVariationPerStep_;// intermediate threshold COVs + OT::Sample eventInputSample_;// event input sample + OT::Sample eventOutputSample_;// event output sample + SamplePersistentCollection allLevelSample_; // all event output sample + SamplePersistentCollection allPointSample_; // all event output sample + OT::Distribution thresholdDistribution_;//distribution of the final threshold value + + // attributes used for conveniency, not to be saved/loaded + OT::StandardEvent standardEvent_;// the algorithm happens in U + OT::UnsignedInteger dimension_;// input dimension + OT::Sample currentPointSample_;// X + OT::Sample currentLevelSample_;//f(X) + +}; /* class SubsetInverseSampling */ + +} /* namespace OTROBOPT */ + +#endif /* OTROBOPT_SUBSETINVERSESAMPLING_HXX */ diff --git a/lib/src/otrobopt/SubsetInverseSamplingResult.hxx b/lib/src/otrobopt/SubsetInverseSamplingResult.hxx new file mode 100644 index 0000000..c54de7c --- /dev/null +++ b/lib/src/otrobopt/SubsetInverseSamplingResult.hxx @@ -0,0 +1,73 @@ +// -*- C++ -*- +/** + * @brief SubsetInverseSamplingResult + * + * Copyright 2005-2023 Airbus-EDF-IMACS-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 OTROBOPT_SUBSETINVERSESAMPLINGRESULT_HXX +#define OTROBOPT_SUBSETINVERSESAMPLINGRESULT_HXX + +#include +#include "otrobopt/OTRobOptprivate.hxx" + +namespace OTROBOPT +{ + +class OTROBOPT_API SubsetInverseSamplingResult +: public OT::ProbabilitySimulationResult +{ + CLASSNAME; +public: + + /** Constructor with parameters */ + SubsetInverseSamplingResult(); + + /** Constructor with parameters */ + SubsetInverseSamplingResult(const OT::RandomVector & event, + const OT::Scalar probabilityEstimate, + const OT::Scalar varianceEstimate, + const OT::UnsignedInteger outerSampling, + const OT::UnsignedInteger blockSize, + const OT::Scalar coefficientOfVariation = 0.0, + const OT::Scalar threshold = 0.0); + + /** Virtual constructor */ + SubsetInverseSamplingResult * clone() const override; + + /** Coefficient of variation estimate accessor */ + OT::Scalar getCoefficientOfVariation() const override; + + /** String converter */ + OT::String __repr__() const override; + + /** Method save() stores the object through the StorageManager */ + void save(OT::Advocate & adv) const override; + + /** Method load() reloads the object from the StorageManager */ + void load(OT::Advocate & adv) override; + +protected: + OT::Scalar coefficientOfVariation_; + OT::Scalar threshold_; + +private: + +}; /* class SubsetInverseSamplingResult */ + +} /* namespace OTROBOPT */ + +#endif /* OTROBOPT_SUBSETINVERSESAMPLINGRESULT_HXX */ diff --git a/lib/test/CMakeLists.txt b/lib/test/CMakeLists.txt index 139c0db..b529eb7 100644 --- a/lib/test/CMakeLists.txt +++ b/lib/test/CMakeLists.txt @@ -73,6 +73,8 @@ ot_check_test (MeasureEvaluation_std) if (NOT WIN32) ot_check_test (SequentialMonteCarloRobustAlgorithm_std) endif() +ot_check_test (SubsetInverseSampling_std) +ot_check_test (InverseFORM_std IGNOREOUT) add_custom_target ( cppcheck COMMAND ${CMAKE_CTEST_COMMAND} -R "^cppcheck_" DEPENDS ${CHECK_TO_BE_RUN} diff --git a/lib/test/t_InverseFORM_std.cxx b/lib/test/t_InverseFORM_std.cxx new file mode 100644 index 0000000..2268e74 --- /dev/null +++ b/lib/test/t_InverseFORM_std.cxx @@ -0,0 +1,156 @@ +#include "openturns/OT.hxx" +#include "openturns/OTtestcode.hxx" +#include "otrobopt/OTRobOpt.hxx" + +using namespace OT; +using namespace OT::Test; +using namespace OTROBOPT; + + +int main() +{ + TESTPREAMBLE; + + try { + Log::Show(Log::INFO); + + OStream fullprint(std::cout); + { + // beam model + UnsignedInteger dim = 4; + + SymbolicFunction function(Description({"E", "F", "L", "b", "h"}), Description({"F*L^3/(48.*E*b*h^3/12.)"})); + ParametricFunction parametric(function, Indices({2}), Point({5.0})); + + ComposedDistribution::DistributionCollection coll; + coll.add(LogNormalMuSigmaOverMu(30000., 0.12, 0.).getDistribution());//E + coll.add(LogNormalMuSigmaOverMu(0.1, 0.2, 0.).getDistribution());//F + coll.add(LogNormalMuSigmaOverMu(0.2, 0.05, 0.).getDistribution());//b + coll.add(LogNormalMuSigmaOverMu(0.4, 0.05, 0.).getDistribution());//h + + const ComposedDistribution myDistribution(coll); + + + Point median(dim); + for(UnsignedInteger i = 0; i < dim; ++ i) + median[i] = myDistribution.getMarginal(i).computeQuantile(0.5)[0]; + + /* We create a 'usual' RandomVector from the Distribution */ + const RandomVector vect(myDistribution); + + /* We create a composite random vector */ + CompositeRandomVector output(parametric, vect); + + /* We create an Event from this RandomVector */ + ThresholdEvent myEvent(output, Greater(), 0.02); + + fullprint << "event=" << myEvent << std::endl; + + // create the algorithm + InverseFORM algo(myEvent, "L", median); + + // print result + algo.run(); + const InverseFORMResult result(algo.getResult()); + fullprint << "Inverse FORM result=" << result << std::endl; + assert_almost_equal(result.getParameter(), Point({5.44343})); + + // FORM must yield the same probability on the limit state with parameter set to the optimum + parametric.setParameter(result.getParameter()); + output = CompositeRandomVector(parametric, vect); + myEvent = ThresholdEvent(output, Greater(), 0.02); + + Cobyla cobyla; + cobyla.setMaximumIterationNumber(100); + cobyla.setMaximumAbsoluteError(1.0e-10); + cobyla.setMaximumRelativeError(1.0e-10); + cobyla.setMaximumResidualError(1.0e-10); + cobyla.setMaximumConstraintError(1.0e-10); + + FORM myFORM(cobyla, myEvent, median); + myFORM.run(); + const FORMResult resultFORM(myFORM.getResult()); + fullprint << "FORM result="<< resultFORM << std::endl; + assert_almost_equal(resultFORM.getHasoferReliabilityIndex(), algo.getTargetBeta()); + } + + { + // sphere under pressure + UnsignedInteger dim = 3; + + SymbolicFunction function(Description({"R", "e", "mulog_e", "p"}), Description({"700.0-p*R/(2.*e)"})); + Indices set = {2}; + + const Scalar L0 = -4.715; + ParametricFunction parametric(function, Indices({2}), Point({L0})); + + Dirac mulog_eDist(L0); + mulog_eDist.setDescription(Description(1, "mulog_e")); + ComposedDistribution::DistributionCollection eColl; eColl.add(mulog_eDist); eColl.add(Dirac(0.1)); eColl.add(Dirac(0.)); + ComposedDistribution eParams(eColl); + + ComposedDistribution::DistributionCollection coll; + coll.add(Beta(0.117284, 0.117284, 2.9, 3.1));//R + ConditionalDistribution eDist(LogNormal(L0, 0.1, 0.), eParams); + + coll.add(eDist);//e + coll.add(WeibullMin(3.16471, 9.21097, 0.0));//p + ComposedDistribution myDistribution(coll); + + Point median(dim); + for(UnsignedInteger i = 0; i < dim; ++ i) + median[i] = myDistribution.getMarginal(i).computeQuantile(0.5)[0]; + + /* We create a 'usual' RandomVector from the Distribution */ + RandomVector vect(myDistribution); + + /* We create a composite random vector */ + CompositeRandomVector output(parametric, vect); + + /* We create an Event from this RandomVector */ + ThresholdEvent myEvent(output, Less(), 0.); + + fullprint << "event=" << myEvent << std::endl; + + // create the algorithm + InverseFORM algo(myEvent, "mulog_e", median); + + // print result + algo.run(); + const InverseFORMResult result(algo.getResult()); + fullprint << "Inverse FORM result=" << result << std::endl; + assert_almost_equal(result.getParameter(), Point({-4.69056}), 0.0, 0.03); + + // FORM must yield the same probability on the limit state with parameter set to the optimum + eColl[0] = Dirac(result.getParameter()[0]); + eParams = ComposedDistribution(eColl); + eDist = ConditionalDistribution(LogNormal(result.getParameter()[0], 0.1, 0.0), eParams); + coll[1] = eDist; + myDistribution = ComposedDistribution(coll); + vect = RandomVector(myDistribution); + parametric.setParameter(result.getParameter()); + output = CompositeRandomVector(parametric, vect); + myEvent = ThresholdEvent(output, Less(), 0.); + + Cobyla cobyla; + cobyla.setMaximumIterationNumber(100); + cobyla.setMaximumAbsoluteError(1.0e-10); + cobyla.setMaximumRelativeError(1.0e-10); + cobyla.setMaximumResidualError(1.0e-10); + cobyla.setMaximumConstraintError(1.0e-10); + + FORM myFORM(cobyla, myEvent, median); + myFORM.run(); + const FORMResult resultFORM(myFORM.getResult()); + fullprint << "result FORM="< +#include +#include "otrobopt/SubsetInverseSampling.hxx" + +using namespace OT; +using namespace OTROBOPT; + +int main() +{ + Distribution myDistribution = Normal(); + RandomVector vect(myDistribution); + + SymbolicFunction limitState("X","2*X"); + CompositeRandomVector output(limitState, vect); + + ThresholdEvent myEvent(output, Less(), 0.0); + + Scalar finalProbability(0.0001); + + SubsetInverseSampling ssi(myEvent, finalProbability); + std::cout << ssi << std::endl; + ssi.run(); + std::cout << ssi.getResult() << std::endl; + return 0; +} diff --git a/lib/test/t_SubsetInverseSampling_std.expout b/lib/test/t_SubsetInverseSampling_std.expout new file mode 100644 index 0000000..8e68f7d --- /dev/null +++ b/lib/test/t_SubsetInverseSampling_std.expout @@ -0,0 +1,2 @@ +class=SubsetInverseSampling derived from class=EventSimulation event=class=RandomVector implementation=class=ThresholdEventImplementation antecedent=class=CompositeRandomVector function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[X,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] antecedent=class=UsualRandomVector distribution=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] operator=class=Less name=Unnamed threshold=0 maximumOuterSampling=10000 maximumCoefficientOfVariation=0.1 maximumStandardDeviation=0 blockSize=1 event=class=RandomVector implementation=class=ThresholdEventImplementation antecedent=class=CompositeRandomVector function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[X,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[X] outputVariablesNames=[y0] formulas=[2*X] antecedent=class=UsualRandomVector distribution=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] operator=class=Less name=Unnamed threshold=0 targetProbability=0.0001 proposalRange=2 conditionalProbability=0.1 keepEventSample_=false +probabilityEstimate=9.999990e-05 varianceEstimate=9.138324e-11 standard deviation=9.56e-06 coefficient of variation=9.56e-02 confidenceLength(0.95)=3.75e-05 outerSampling=40000 blockSize=1 threshold=-7.56294 diff --git a/python/doc/examples/plot_example2.py b/python/doc/examples/plot_example2.py index c4796da..e658f36 100755 --- a/python/doc/examples/plot_example2.py +++ b/python/doc/examples/plot_example2.py @@ -5,7 +5,7 @@ # %% # Problem statement: -# +# # .. math:: # # \begin{aligned} @@ -15,7 +15,7 @@ # & & 2x_1 + 4x_0 - 120 \leq 0 \\ # & & & \Theta \thicksim \cN(1, 3) # \end{aligned} -# +# # Solution: :math:`\hat{x} = [15, 30]` # %% @@ -49,8 +49,8 @@ algo.setInitialSearch(100) algo.run() result = algo.getResult() -#openturns.testing.assert_almost_equal( - #result.getOptimalPoint(), [46.2957, 46.634], 1e-3) +# openturns.testing.assert_almost_equal( +# result.getOptimalPoint(), [46.2957, 46.634], 1e-3) print('x*=', result.getOptimalPoint(), 'J(x*)=', result.getOptimalValue(), 'iteration=', result.getIterationNumber()) diff --git a/python/doc/user_manual/user_manual.rst b/python/doc/user_manual/user_manual.rst index 2314a3a..84cfdb5 100644 --- a/python/doc/user_manual/user_manual.rst +++ b/python/doc/user_manual/user_manual.rst @@ -117,3 +117,19 @@ Solve a robust optimization problem RobustOptimizationAlgorithm SequentialMonteCarloRobustAlgorithm + +.. FIXME: sphinx.errors.SphinxWarning: .../otrobopt.py:docstring of openturns.analytical.AnalyticalResult.getHasoferReliabilityIndexSensitivity:4:undefined label: sensitivity_form +.. _sensitivity_form: +.. _importance_form: + +Solve an inverse reliability problem +------------------------------------ + +.. autosummary:: + :toctree: _generated/ + :template: class.rst_t + + SubsetInverseSampling + SubsetInverseSamplingResult + InverseFORM + InverseFORMResult diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt index 2f0af07..e6cd867 100644 --- a/python/src/CMakeLists.txt +++ b/python/src/CMakeLists.txt @@ -93,6 +93,10 @@ ot_add_python_module( ${PACKAGE_NAME} ${PACKAGE_NAME}_module.i RobustOptimizationProblem.i RobustOptimizationProblem_doc.i.in RobustOptimizationAlgorithm.i RobustOptimizationAlgorithm_doc.i.in SequentialMonteCarloRobustAlgorithm.i SequentialMonteCarloRobustAlgorithm_doc.i.in + SubsetInverseSampling.i SubsetInverseSampling_doc.i.in + SubsetInverseSamplingResult.i SubsetInverseSamplingResult_doc.i.in + InverseFORMResult.i InverseFORMResult_doc.i.in + InverseFORM.i InverseFORM_doc.i.in ) diff --git a/python/src/InverseFORM.i b/python/src/InverseFORM.i new file mode 100644 index 0000000..1008529 --- /dev/null +++ b/python/src/InverseFORM.i @@ -0,0 +1,12 @@ +// SWIG file InverseFORM.i + +%{ +#include "otrobopt/InverseFORM.hxx" +%} + +%include InverseFORM_doc.i + +%include otrobopt/InverseFORM.hxx + +namespace OTROBOPT { %extend InverseFORM { InverseFORM(const InverseFORM & other) { return new OTROBOPT::InverseFORM(other); } } } + diff --git a/python/src/InverseFORMResult.i b/python/src/InverseFORMResult.i new file mode 100644 index 0000000..16527da --- /dev/null +++ b/python/src/InverseFORMResult.i @@ -0,0 +1,11 @@ +// SWIG file InverseFORMResult.i + +%{ +#include "otrobopt/InverseFORMResult.hxx" +%} + +%include InverseFORMResult_doc.i + +%include otrobopt/InverseFORMResult.hxx + +namespace OTROBOPT { %extend InverseFORMResult { InverseFORMResult(const InverseFORMResult & other) { return new OTROBOPT::InverseFORMResult(other); } } } diff --git a/python/src/InverseFORMResult_doc.i.in b/python/src/InverseFORMResult_doc.i.in new file mode 100644 index 0000000..b2a9488 --- /dev/null +++ b/python/src/InverseFORMResult_doc.i.in @@ -0,0 +1,32 @@ +%feature("docstring") OTROBOPT::InverseFORMResult +"Result of Inverse FORM." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORMResult::getParameter +"Parameter accessor. + +Returns +------- +parameter : :py:class:`openturns.Point` + Parameter value." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORMResult::getParameterDescription +"Parameter description accessor. + +Returns +------- +description : :py:class:`openturns.Description` + Parameter description." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORMResult::getConvergenceCriteria +"Convergence criteria accessor. + +Returns +------- +description : :py:class:`openturns.Point` + Parameter description." diff --git a/python/src/InverseFORM_doc.i.in b/python/src/InverseFORM_doc.i.in new file mode 100644 index 0000000..4ee4c03 --- /dev/null +++ b/python/src/InverseFORM_doc.i.in @@ -0,0 +1,112 @@ +%feature("docstring") OTROBOPT::InverseFORM +"Inverse FORM. + +Parameters +---------- +event : :class:`~openturns.RandomVector` + Event we are computing the probability of. +parameterName : string + The optimized parameter +physicalStartingPoint : sequence of float + Optimization starting point +" + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getResult +"Result accessor. + +Returns +------- +result : :class:`~otrobopt.InverseFORMResult` + Result." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getVariableConvergence +"Variable convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getVariableStepMaxIterations +"Line search maximum iteration accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::run +"Run algorithm." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setBetaConvergence +"Beta convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getBetaConvergence +"Beta convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setFixedStep +"Fixed-step mode accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getFixedStep +"Fixed-step mode accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setFixedStepValue +"Fixed-step value accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getFixedStepValue +"Fixed-step value accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setLimitStateConvergence +"Limit-state convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getLimitStateConvergence +"Limit-state convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setMaximumIteration +"Maximum iteration accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getMaximumIteration +"Maximum iteration accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setResult +"Result accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setTargetBeta +"Target beta accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::getTargetBeta +"Target beta accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setVariableConvergence +"Variable convergence criterion accessor." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::InverseFORM::setVariableStepMaxIterations +"Line search maximum iteration accessor." diff --git a/python/src/SubsetInverseSampling.i b/python/src/SubsetInverseSampling.i new file mode 100644 index 0000000..0e5146f --- /dev/null +++ b/python/src/SubsetInverseSampling.i @@ -0,0 +1,10 @@ +// SWIG file SubsetInverseSampling.i + +%{ +#include "otrobopt/SubsetInverseSampling.hxx" +%} + +%include SubsetInverseSampling_doc.i + +%include otrobopt/SubsetInverseSampling.hxx +namespace OTROBOPT { %extend SubsetInverseSampling { SubsetInverseSampling(const SubsetInverseSampling & other) { return new OTROBOPT::SubsetInverseSampling(other); } } } diff --git a/python/src/SubsetInverseSamplingResult.i b/python/src/SubsetInverseSamplingResult.i new file mode 100644 index 0000000..5f9b88d --- /dev/null +++ b/python/src/SubsetInverseSamplingResult.i @@ -0,0 +1,8 @@ +// SWIG file SubsetInverseSamplingResult.i + +%{ +#include "otrobopt/SubsetInverseSamplingResult.hxx" +%} + +%include otrobopt/SubsetInverseSamplingResult.hxx +namespace OTROBOPT { %extend SubsetInverseSamplingResult { SubsetInverseSamplingResult(const SubsetInverseSamplingResult & other) { return new OTROBOPT::SubsetInverseSamplingResult(other); } } } diff --git a/python/src/SubsetInverseSamplingResult_doc.i.in b/python/src/SubsetInverseSamplingResult_doc.i.in new file mode 100644 index 0000000..298b196 --- /dev/null +++ b/python/src/SubsetInverseSamplingResult_doc.i.in @@ -0,0 +1,44 @@ +//%define OTSUBSETINVERSE_SubsetSampling_doc +//"SubsetSampling class. +// +//Examples +//-------- +//>>> import otsubsetinverse +//>>> a = otsubsetinverse.SubsetSampling() +// +//Notes +//----- +//Template module class." +//%enddef +//%feature("docstring") OTSUBSETINVERSE::SubsetSamplingResult +//OTSUBSETINVERSE_SubsetSampling_doc +// +//// --------------------------------------------------------------------- +// +//%define OTSUBSETINVERSE_SubsetSampling_square_doc +//"Square of a point. +// +//Parameters +//---------- +//point : sequence of float +// Input point +// +//Returns +//------- +//square : :py:class:`openturns.NumericalPoint` +// Square of point +// +//Examples +//-------- +//>>> import openturns as ot +//>>> import otsubsetinverse +//>>> a = otsubsetinverse.SubsetSampling() +//>>> point = ot.NumericalPoint([1.0, 2.0, 3.0]) +//>>> a.square(point) +//class=NumericalPoint name=Unnamed dimension=3 values=[1,4,9]" +//%enddef +//%feature("docstring") OTSUBSETINVERSE::SubsetSamplingResult::square +//OTSUBSETINVERSE_SubsetSampling_square_doc +// +// +// \ No newline at end of file diff --git a/python/src/SubsetInverseSampling_doc.i.in b/python/src/SubsetInverseSampling_doc.i.in new file mode 100644 index 0000000..5c9b6de --- /dev/null +++ b/python/src/SubsetInverseSampling_doc.i.in @@ -0,0 +1,322 @@ +%feature("docstring") OTROBOPT::SubsetInverseSampling +"Subset inverse simulation. + +Parameters +---------- +event : :class:`~openturns.Event` + Event we are computing the probability of. The threshold of the event is + not used. +targetProbability : float + The wanted final probability. +proposalRange : float, optional + Proposal range length +conditionalProbability : float, optional + Value of :math:`P(F_i|F_{i-1})` between successive steps + +Notes +----- +The goal is to estimate the threshold of the following target probability : + +.. math:: + + P_f = \int_{\mathcal D_f} f_{\uX}(\ux)\di{\ux}\\ + = \int_{\mathbb R^{n_X}} \mathbf{1}_{\{g(\ux,\underline{d}) \:\leq 0\: \}}f_{\uX}(\ux)\di{\ux}\\ + = \Prob {\{g(\uX,\underline{d}) \leq q\}} + + +The idea of the subset simulation method is to replace simulating a +rare failure event in the original probability space by a sequence of +simulations of more frequent conditional events :math:`F_i` : + +.. math:: + + F_1 \supset F_2 \supset \dots \supset F_m = F + + +The original probability estimate rewrites : + +.. math:: + + P_f = P(F_m) = P(\bigcap \limits_{i=1}^m F_i) = P(F_1) \prod_{i=2}^m P(F_i|F_{i-1}) + + +And each conditional subset failure region is chosen by setting the threshold +:math:`g_i` so that :math:`P(F_i|F_{i-1})` leads to a conditional failure +probability of order :math:`0.1` : + +.. math:: + + F_i =\Prob {\{g(\uX,\underline{d}) \leq g_i\}} + + +The conditional samples are generated by the means of Markov Chains, +using the Metropolis Hastings algorithm. + +:math:`N` being the number of simulations per subset, and :math:`p_{0i}` the +conditional probability of each subset event, and :math:`\gamma_i` the +autocorrelation between Markov chain samples. + +.. math:: + + \delta^2 = \sum_{i=1}^m \delta^2_i = \sum_{i=1}^m (1+\gamma_i) \frac{1-p_{0i}}{p_{0i}N} + + +The first event :math:`F_1` not being conditional, :math:`\delta^2_1` +expresses as the classic Monte Carlo c.o.v. + +Examples +-------- +>>> import openturns as ot +>>> import otrobopt + +>>> ot.RandomGenerator.SetSeed(0) + +Create a performance function with an associated distribution. + +>>> limitState = ot.SymbolicFunction(['u1', 'u2'], ['u1-u2']) +>>> dim = limitState.getInputDimension() +>>> mean = ot.Point([7., 2.]) +>>> sigma = ot.Point(dim, 1.0) +>>> R = ot.IdentityMatrix(dim) +>>> myDistribution = ot.Normal(mean, sigma, R) +>>> vect = ot.RandomVector(myDistribution) +>>> output = ot.CompositeRandomVector(limitState, vect) + +Create an event with a fictional threshold value which will not be used. + +>>> myEvent = ot.ThresholdEvent(output, ot.Less(), 0.) + +Define the target probability for which the threshold will be computed. + +>>> targetProbability = 0.0002 +>>> mySS = otrobopt.SubsetInverseSampling(myEvent, targetProbability) +>>> mySS.setMaximumOuterSampling(10000) +>>> mySS.run() + +Get some results. + +>>> resultSS = mySS.getResult() +>>> print('Pf = {}'.format(resultSS.getProbabilityEstimate())) +Pf = 0.0002... +>>> print('Threshold = {}'.format(mySS.getThresholdPerStep()[-1])) +Threshold = 0.02027... +>>> print('Threshold confidence length = {}'.format( +... mySS.getThresholdConfidenceLength(0.90))) +Threshold confidence length = 0.0245... + +See also +-------- +openturns.Simulation" + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getCoefficientOfVariationPerStep +"Coefficient of variation per step accessor. + +Returns +------- +coef : :py:class:`openturns.Point` + Coefficient of variation at each subset step." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getThresholdCoefficientOfVariationPerStep +"Threshold coefficient of variation per step accessor. + +Returns +------- +coef : :py:class:`openturns.Point` + Coefficient of variation at each subset step." +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setConditionalProbability +"Conditional probability accessor. + +Value of :math:`P(F_i|F_{i-1})` between successive steps. + +Parameters +---------- +prob : float + Conditional probability value." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getConditionalProbability +"Conditional probability accessor. + +Value of :math:`P(F_i|F_{i-1})` between successive steps. + +Returns +------- +prob : float + Conditional probability value." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setTargetProbability +"Final target probability accessor. + +Value of :math:`P(F_m)`. + +Parameters +---------- +prob : float + Final target probability value." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getThresholdConfidenceLength +"Accessor to the confidence length of the threshold. + +Parameters +---------- +level : float, :math:`level \in ]0, 1[` + Confidence level. By default, it is :math:`0.95`. + +Returns +------- +confidenceLength : float + Length of the confidence interval at the confidence level *level*." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getTargetProbability +"Final target probability accessor. + +Value of :math:`P(F_m)`. + +Returns +------- +prob : float + Final target probability value." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setKeepEventSample +"Sample storage accessor. + +Parameters +---------- +prob : bool + Whether to keep the event samples." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getEventInputSample +"Input sample accessor. + +Returns +------- +inputSample : :py:class:`openturns.Sample` + Input sample." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getEventOutputSample +"Output sample accessor. + +Returns +------- +outputSample : :py:class:`openturns.Sample` + Ouput sample." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getInputSample +"All input sample accessor. + +Returns +------- +inputSample : :py:class:`openturns.Sample` + Input sample." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getOutputSample +"All output sample accessor. + +Returns +------- +outputSample : :py:class:`openturns.Sample` + Output sample." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getThresholdPerStep +"Threshold accessor. + +Returns +------- +threshold : :py:class:`openturns.Point` + Threshold values." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getGammaPerStep +"Autocorrelation accessor. + +Returns +------- +prob : :py:class:`openturns.Point` + Autocorrelation values." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getProbabilityEstimatePerStep +"Probability estimate accessor. + +Returns +------- +prob : :py:class:`openturns.Point` + Probability estimate values." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getNumberOfSteps +"Subset steps number accesor. + +Returns +------- +n : int + Number of subset steps." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setProposalRange +"Proposal range length accessor. + +Parameters +---------- +range : float + Range length." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::getProposalRange +"Proposal range length accessor. + +Returns +------- +range : float + Range length." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setBetaMin +"Radius of the hypershere accessor. + +Parameters +---------- +beta : float + Radius value of the exclusion hypershere when the conditional simulation is activated." + +// --------------------------------------------------------------------------- + +%feature("docstring") OTROBOPT::SubsetInverseSampling::setISubset +"Conditonal simulation activation accessor. + +Parameters +---------- +isubset : bool + Activate or not the conditional simulation for the first step of the + simulation." diff --git a/python/src/__init__.py b/python/src/__init__.py index 15d4bc1..1b5070f 100644 --- a/python/src/__init__.py +++ b/python/src/__init__.py @@ -15,4 +15,4 @@ from .otrobopt import * -__version__ = '0.12' +__version__ = '0.13' diff --git a/python/src/otrobopt_module.i b/python/src/otrobopt_module.i index 3be8d0d..fe9b24f 100644 --- a/python/src/otrobopt_module.i +++ b/python/src/otrobopt_module.i @@ -33,4 +33,7 @@ %include RobustOptimizationProblem.i %include RobustOptimizationAlgorithm.i %include SequentialMonteCarloRobustAlgorithm.i - +%include SubsetInverseSamplingResult.i +%include SubsetInverseSampling.i +%include InverseFORMResult.i +%include InverseFORM.i diff --git a/python/test/CMakeLists.txt b/python/test/CMakeLists.txt index 4693afb..5311a53 100644 --- a/python/test/CMakeLists.txt +++ b/python/test/CMakeLists.txt @@ -66,6 +66,9 @@ ot_pyinstallcheck_test (MeasureEvaluation_std) ot_pyinstallcheck_test (MeasureFactory_std) ot_pyinstallcheck_test (MeasureFunction_std) ot_pyinstallcheck_test (SequentialMonteCarloRobustAlgorithm_std) +ot_pyinstallcheck_test (SubsetInverseSampling_R-S) +ot_pyinstallcheck_test (SubsetInverseSampling_Waarts_system_series) +ot_pyinstallcheck_test (InverseFORM_std IGNOREOUT) ot_pyinstallcheck_test (docstring) if (MATPLOTLIB_FOUND) diff --git a/python/test/sequential_mc.py b/python/test/sequential_mc.py index 255dc95..e01aaaa 100644 --- a/python/test/sequential_mc.py +++ b/python/test/sequential_mc.py @@ -117,6 +117,7 @@ def solve(self): 0, self.J_.getInputDimension() - self.distributionXi_.getDimension()) currentSampleXi = ot.Sample(0, distributionXi.getDimension()) currentN = self.N0_ + currentPoint = None for iteration in range(self.robustIteration_): if self.verbose_: print("start iteration", iteration) diff --git a/python/test/t_InverseFORM_std.py b/python/test/t_InverseFORM_std.py new file mode 100755 index 0000000..6d38abe --- /dev/null +++ b/python/test/t_InverseFORM_std.py @@ -0,0 +1,33 @@ +#! /usr/bin/env python + +import openturns as ot +import openturns.testing as ott +import otrobopt + +ot.Log.Show(ot.Log.INFO) + +# beam model, what else +full = ot.SymbolicFunction(['E', 'F', 'L', 'b', 'h'], ['F*L^3/(48.*E*b*h^3/12.)']) +g = ot.ParametricFunction(full, [2], [5.0]) # freeze L=5.0 + +coll = [] +coll.append(ot.LogNormalMuSigmaOverMu(30000., 0.12, 0.).getDistribution()) # E +coll.append(ot.LogNormalMuSigmaOverMu(0.1, 0.2, 0.).getDistribution()) # F +coll.append(ot.LogNormalMuSigmaOverMu(0.2, 0.05, 0.).getDistribution()) # b +coll.append(ot.LogNormalMuSigmaOverMu(0.4, 0.05, 0.).getDistribution()) # h + +distribution = ot.ComposedDistribution(coll) + +x0 = [coll[i].computeQuantile(0.5)[0] for i in range(len(coll))] + +vect = ot.RandomVector(distribution) + +output = ot.CompositeRandomVector(g, vect) + +event = ot.ThresholdEvent(output, ot.Greater(), 0.02) + +algo = otrobopt.InverseFORM(event, 'L', x0) +algo.run() +result = algo.getResult() +print(result) +ott.assert_almost_equal(result.getParameter(), [5.44343]) diff --git a/python/test/t_MeasureEvaluation_std.py b/python/test/t_MeasureEvaluation_std.py index 712e579..d44c033 100755 --- a/python/test/t_MeasureEvaluation_std.py +++ b/python/test/t_MeasureEvaluation_std.py @@ -56,7 +56,7 @@ otrobopt.IndividualChanceMeasure( f, thetaDist, ot.GreaterOrEqual(), [0.5]), otrobopt.MeanStandardDeviationTradeoffMeasure(f, thetaDist, [0.8]), - #otrobopt.QuantileMeasure(f, thetaDist, 0.5) + # otrobopt.QuantileMeasure(f, thetaDist, 0.5) ] aggregated = otrobopt.AggregatedMeasure(measures) measures.append(aggregated) diff --git a/python/test/t_SubsetInverseSampling_R-S.expout b/python/test/t_SubsetInverseSampling_R-S.expout new file mode 100644 index 0000000..d23f17c --- /dev/null +++ b/python/test/t_SubsetInverseSampling_R-S.expout @@ -0,0 +1,25 @@ + +******************************************************************* +***************************** MONTE CARLO ************************* +******************************************************************* +Pf estimation = 2.00000e-04 +Pf Variance estimation = 1.99960e-10 +CoV = 0.07070 +90% Confidence Interval = 4.65188e-05 +CI at 90% =[ 1.76741e-04 ; 2.23259e-04 ] +Threshold = 0.00000e+00 +Limit state calls = 1000000 +******************************************************************* + +******************************************************************* +************************** SUBSET SAMPLING ************************ +******************************************************************* +Pf estimation = 2.00000e-04 +Pf Variance estimation = 3.55352e-10 +CoV = 0.09425 +90% Confidence Interval = 6.20136e-05 +CI at 90% =[ 1.68993e-04 ; 2.31007e-04 ] +Threshold = -4.97416e-03 +Limit state calls = 40000 +******************************************************************* + diff --git a/python/test/t_SubsetInverseSampling_R-S.py b/python/test/t_SubsetInverseSampling_R-S.py new file mode 100755 index 0000000..c1b2ada --- /dev/null +++ b/python/test/t_SubsetInverseSampling_R-S.py @@ -0,0 +1,131 @@ +#! /usr/bin/env python + +# This script shows the use of the module SubsetInverse on the model G=R-S, +# with R~N, S~N + +import openturns as ot +import otrobopt + +ot.TESTPREAMBLE() + + +def cleanScalar(inScalar): + if abs(inScalar) < 1.0e-10: + inScalar = 0.0 + return inScalar + + +########################################################################### +# Physical model +########################################################################### + +limitState = ot.SymbolicFunction(["u1", "u2"], ["u1-u2"]) +dim = limitState.getInputDimension() + +########################################################################### +# Probabilistic model +########################################################################### + +mean = [7.0, 2.0] +sigma = [1.0] * 2 + +R = ot.IdentityMatrix(dim) +myDistribution = ot.Normal(mean, sigma, R) + +start = myDistribution.getMean() + +########################################################################### +# Limit state +########################################################################### + +vect = ot.RandomVector(myDistribution) + +output = ot.CompositeRandomVector(limitState, vect) + +threshold = 0.0 +myEvent = ot.ThresholdEvent(output, ot.Less(), threshold) + +########################################################################### +# Computation Monte Carlo +########################################################################### +bs = 1 + +# Monte Carlo +experiment = ot.MonteCarloExperiment() +myMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment) +myMC.setMaximumOuterSampling(int(1e6) // bs) +myMC.setBlockSize(bs) +myMC.setMaximumCoefficientOfVariation(-1) +myMC.run() + +########################################################################### +# Results Monte Carlo +########################################################################### + +ResultMC = myMC.getResult() +PFMC = ResultMC.getProbabilityEstimate() +CVMC = ResultMC.getCoefficientOfVariation() +Variance_PF_MC = ResultMC.getVarianceEstimate() +length90MC = ResultMC.getConfidenceLength(0.90) +N_MC = ResultMC.getOuterSampling() * ResultMC.getBlockSize() + +########################################################################### +# Computation SubsetSampling +########################################################################### + +finalTargetProbability = 0.0002 +mySS = otrobopt.SubsetInverseSampling(myEvent, finalTargetProbability) +mySS.setMaximumOuterSampling(10000 // bs) +mySS.setBlockSize(bs) +mySS.run() + +########################################################################### +# Results SubsetSampling +########################################################################### + +ResultSS = mySS.getResult() +PFSS = ResultSS.getProbabilityEstimate() +CVSS = ResultSS.getCoefficientOfVariation() +Variance_PF_SS = ResultSS.getVarianceEstimate() +length90SS = ResultSS.getConfidenceLength(0.90) +N_SS = ResultSS.getOuterSampling() * ResultSS.getBlockSize() +thresholdSS = mySS.getThresholdPerStep()[-1] +thLengthSS = mySS.getThresholdConfidenceLength(0.90) + +########################################################################### + +print("") +print("*******************************************************************") +print("***************************** MONTE CARLO *************************") +print("*******************************************************************") +print("Pf estimation = %.5e" % PFMC) +print("Pf Variance estimation = %.5e" % Variance_PF_MC) +print("CoV = %.5f" % CVMC) +print("90% Confidence Interval =", "%.5e" % length90MC) +print( + "CI at 90% =[", + "%.5e" % (PFMC - 0.5 * length90MC), + "; %.5e" % (PFMC + 0.5 * length90MC), + "]", +) +print("Threshold = %.5e" % threshold) +print("Limit state calls =", N_MC) +print("*******************************************************************") +print("") +print("*******************************************************************") +print("************************** SUBSET SAMPLING ************************") +print("*******************************************************************") +print("Pf estimation = %.5e" % PFSS) +print("Pf Variance estimation = %.5e" % Variance_PF_SS) +print("CoV = %.5f" % CVSS) +print("90% Confidence Interval =", "%.5e" % length90SS) +print( + "CI at 90% =[", + "%.5e" % (PFSS - 0.5 * length90SS), + "; %.5e" % (PFSS + 0.5 * length90SS), + "]", +) +print("Threshold = %.5e" % thresholdSS) +print("Limit state calls =", N_SS) +print("*******************************************************************") +print("") diff --git a/python/test/t_SubsetInverseSampling_Waarts_system_series.expout b/python/test/t_SubsetInverseSampling_Waarts_system_series.expout new file mode 100644 index 0000000..a2c8b01 --- /dev/null +++ b/python/test/t_SubsetInverseSampling_Waarts_system_series.expout @@ -0,0 +1,25 @@ + +********************************************************************* +********************** MONTE CARLO ********************************** +********************************************************************* +Pf estimation = 2.21600e-03 +Pf Variance estimation = 2.21109e-09 +CoV = 0.02122 +90% Confidence Interval = 1.54689e-04 +CI at 90% =[ 2.13866e-03 ; 2.29334e-03 ] +Threshold = 0.00000e+00 +Limit state calls = 1000000 +********************************************************************* + +********************************************************************* +******************** SUBSET SAMPLING ******************************** +********************************************************************* +Pf estimation = 2.21500e-03 +Pf Variance estimation = 2.73671e-08 +CoV = 0.07469 +90% Confidence Interval = 5.44217e-04 +CI at 90% =[ 1.94289e-03 ; 2.48711e-03 ] +Threshold = 1.56174e-02 +Limit state calls = 30000 +********************************************************************* + diff --git a/python/test/t_SubsetInverseSampling_Waarts_system_series.py b/python/test/t_SubsetInverseSampling_Waarts_system_series.py new file mode 100755 index 0000000..267514e --- /dev/null +++ b/python/test/t_SubsetInverseSampling_Waarts_system_series.py @@ -0,0 +1,133 @@ +#! /usr/bin/env python + +# This script shows the use of the module SubsetInverse on the Waarts System +# series benchmark +# Taken from P. H. Waarts, Structural reliability using finite element analysis, 2000. + +import openturns as ot +import otrobopt + +ot.TESTPREAMBLE() + + +def cleanScalar(inScalar): + if abs(inScalar) < 1.0e-10: + inScalar = 0.0 + return inScalar + + +########################################################################### +# Physical model +########################################################################### + +formulas = ["min(0.1 * (u1 - u2)^2.0 - (u1 + u2) / sqrt(2.0) + 3.0, 0.1 * (u1 - u2)^2.0 \ + + (u1 + u2) / sqrt(2.0) + 3.0, u1 - u2 + 3.5 * sqrt(2.0), -u1 + u2 + 3.5 * sqrt(2.0))"] +limitState = ot.SymbolicFunction(["u1", "u2"], formulas) +dim = limitState.getInputDimension() + +########################################################################### +# Probabilistic model +########################################################################### + +mean = [0.0] * 2 +sigma = [1.0] * 2 +R = ot.IdentityMatrix(dim) +myDistribution = ot.Normal(mean, sigma, R) + +start = myDistribution.getMean() + +########################################################################### +# Limit state +########################################################################### + +vect = ot.RandomVector(myDistribution) + +output = ot.CompositeRandomVector(limitState, vect) + +threshold = 0.0 +myEvent = ot.ThresholdEvent(output, ot.Less(), threshold) + +########################################################################### +# Computation Monte Carlo +########################################################################### +bs = 1 + +# Monte Carlo +experiment = ot.MonteCarloExperiment() +myMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment) +myMC.setMaximumOuterSampling(int(1e6) // bs) +myMC.setBlockSize(bs) +myMC.setMaximumCoefficientOfVariation(-1.0) +myMC.run() + +########################################################################### +# Results Monte Carlo +########################################################################### + +ResultMC = myMC.getResult() +PFMC = ResultMC.getProbabilityEstimate() +CVMC = ResultMC.getCoefficientOfVariation() +Variance_PF_MC = ResultMC.getVarianceEstimate() +length90MC = ResultMC.getConfidenceLength(0.90) +N_MC = ResultMC.getOuterSampling() * ResultMC.getBlockSize() + +########################################################################### +# Computation SubsetSampling +########################################################################### + +finalTargetProbability = PFMC +mySS = otrobopt.SubsetInverseSampling(myEvent, finalTargetProbability) +mySS.setMaximumOuterSampling(10000 // bs) +mySS.setBlockSize(bs) +mySS.run() + +########################################################################### +# Results SubsetSampling +########################################################################### + +ResultSS = mySS.getResult() +PFSS = ResultSS.getProbabilityEstimate() +CVSS = ResultSS.getCoefficientOfVariation() +Variance_PF_SS = ResultSS.getVarianceEstimate() +length90SS = ResultSS.getConfidenceLength(0.90) +N_SS = ResultSS.getOuterSampling() * ResultSS.getBlockSize() +thresholdSS = mySS.getThresholdPerStep()[-1] +thLengthSS = mySS.getThresholdConfidenceLength(0.90) + +########################################################################### + +print("") +print("*********************************************************************") +print("********************** MONTE CARLO **********************************") +print("*********************************************************************") +print("Pf estimation = %.5e" % PFMC) +print("Pf Variance estimation = %.5e" % Variance_PF_MC) +print("CoV = %.5f" % CVMC) +print("90% Confidence Interval =", "%.5e" % length90MC) +print( + "CI at 90% =[", + "%.5e" % (PFMC - 0.5 * length90MC), + "; %.5e" % (PFMC + 0.5 * length90MC), + "]", +) +print("Threshold = %.5e" % threshold) +print("Limit state calls =", N_MC) +print("*********************************************************************") +print("") +print("*********************************************************************") +print("******************** SUBSET SAMPLING ********************************") +print("*********************************************************************") +print("Pf estimation = %.5e" % PFSS) +print("Pf Variance estimation = %.5e" % Variance_PF_SS) +print("CoV = %.5f" % CVSS) +print("90% Confidence Interval =", "%.5e" % length90SS) +print( + "CI at 90% =[", + "%.5e" % (PFSS - 0.5 * length90SS), + "; %.5e" % (PFSS + 0.5 * length90SS), + "]", +) +print("Threshold = %.5e" % thresholdSS) +print("Limit state calls =", N_SS) +print("*********************************************************************") +print("")