diff --git a/lib/src/InverseFORM.cxx b/lib/src/InverseFORM.cxx index 6b877b9..0d033e3 100755 --- a/lib/src/InverseFORM.cxx +++ b/lib/src/InverseFORM.cxx @@ -228,43 +228,44 @@ Function InverseFORM::getG(const Scalar p) newFunction.setParameter(params); RandomVector antecedent(event_.getImplementation()->getAntecedent().getImplementation()->clone()); const Distribution distribution(antecedent.getDistribution()); - if(distribution.getImplementation()->getClassName() == "ComposedDistribution") +#if OPENTURNS_VERSION >= 102300 + const JointDistribution * p_joint = dynamic_cast(distribution.getImplementation().get()); +#else + const ComposedDistribution * p_joint = dynamic_cast(distribution.getImplementation().get()); +#endif + if (p_joint) { - const ComposedDistribution * p_composed = dynamic_cast(distribution.getImplementation().get()); - if (p_composed) + ComposedDistribution::DistributionCollection distributionCollection(p_joint->getDistributionCollection()); + for (UnsignedInteger i = 0; i < distributionCollection.getSize(); ++ i) { - ComposedDistribution::DistributionCollection distributionCollection(p_composed->getDistributionCollection()); - for (UnsignedInteger i = 0; i < distributionCollection.getSize(); ++ i) + if (distributionCollection[i].getImplementation()->getClassName() == "ConditionalDistribution") { - if (distributionCollection[i].getImplementation()->getClassName() == "ConditionalDistribution") + DistributionImplementation::PointWithDescriptionCollection parametersCollection(distributionCollection[i].getParametersCollection()); + for (UnsignedInteger j = 0; j < parametersCollection.getSize(); ++ j) { - 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 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 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_joint const CompositeRandomVector composite(newFunction, antecedent); const ThresholdEvent event(composite, op, threshold); diff --git a/lib/test/CMakeLists.txt b/lib/test/CMakeLists.txt index 0a1b447..efe0cce 100644 --- a/lib/test/CMakeLists.txt +++ b/lib/test/CMakeLists.txt @@ -48,6 +48,7 @@ ot_check_test (SequentialMonteCarloRobustAlgorithm_std) endif() ot_check_test (SubsetInverseSampling_std) ot_check_test (InverseFORM_std IGNOREOUT) +ot_check_test (InverseFORM_sphere IGNOREOUT) add_custom_target ( cppcheck COMMAND ${CMAKE_CTEST_COMMAND} -R "^cppcheck_" DEPENDS ${CHECK_TO_BE_RUN} diff --git a/lib/test/t_InverseFORM_sphere.cxx b/lib/test/t_InverseFORM_sphere.cxx new file mode 100644 index 0000000..3fea2d2 --- /dev/null +++ b/lib/test/t_InverseFORM_sphere.cxx @@ -0,0 +1,93 @@ +#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; + Log::Show(Log::INFO); + OStream fullprint(std::cout); + + try + { + // sphere under pressure + const UnsignedInteger dim = 3; + const SymbolicFunction function(Description({"R", "e", "mulog_e", "p"}), Description({"700.0-p*R/(2.*e)"})); + + 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 event(output, Less(), 0.); + + fullprint << "event=" << event << std::endl; + + // create the algorithm + InverseFORM algo(event, "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); + event = 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 form(cobyla, event, median); + form.run(); + const FORMResult resultFORM(form.getResult()); + fullprint << "result FORM="<