Skip to content

Commit

Permalink
InverseFORM: Fix code with 1.23
Browse files Browse the repository at this point in the history
  • Loading branch information
@jschueller
jschueller committed May 2, 2024
1 parent 1b17747 commit 38db2f5
Show file tree
Hide file tree
Showing 7 changed files with 132 additions and 116 deletions.
63 changes: 32 additions & 31 deletions lib/src/InverseFORM.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<JointDistribution *>(distribution.getImplementation().get());
#else
const ComposedDistribution * p_joint = dynamic_cast<ComposedDistribution *>(distribution.getImplementation().get());
#endif
if (p_joint)
{
const ComposedDistribution * p_composed = dynamic_cast<ComposedDistribution *>(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<ConditionalDistribution *>(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<ConditionalDistribution *>(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);
Expand Down
1 change: 1 addition & 0 deletions lib/test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
93 changes: 93 additions & 0 deletions lib/test/t_InverseFORM_sphere.cxx
View file
Original file line number Diff line number Diff line change
@@ -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="<<resultFORM << std::endl;
assert_almost_equal(resultFORM.getHasoferReliabilityIndex(), algo.getTargetBeta(), 0.0, 0.1);
}
catch (const TestFailed & ex)
{
std::cerr << ex << std::endl;
return ExitCode::Error;
}
return ExitCode::Success;
}

81 changes: 4 additions & 77 deletions lib/test/t_InverseFORM_std.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,16 +10,15 @@ using namespace OTROBOPT;
int main()
{
TESTPREAMBLE;

try {
Log::Show(Log::INFO);

OStream fullprint(std::cout);

try
{
// beam model
UnsignedInteger dim = 4;
const UnsignedInteger dim = 4;

SymbolicFunction function(Description({"E", "F", "L", "b", "h"}), Description({"F*L^3/(48.*E*b*h^3/12.)"}));
const 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;
Expand All @@ -30,7 +29,6 @@ int main()

const ComposedDistribution myDistribution(coll);


Point median(dim);
for(UnsignedInteger i = 0; i < dim; ++ i)
median[i] = myDistribution.getMarginal(i).computeQuantile(0.5)[0];
Expand Down Expand Up @@ -73,77 +71,6 @@ int main()
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)"}));

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="<<resultFORM << std::endl;
assert_almost_equal(resultFORM.getHasoferReliabilityIndex(), algo.getTargetBeta(), 0.0, 0.1);
}
}
catch (TestFailed & ex)
{
std::cerr << ex << std::endl;
Expand Down
4 changes: 1 addition & 3 deletions python/src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,5 @@ endif ()

# wheel files
configure_file (METADATA.in METADATA @ONLY)
configure_file (WHEEL.in WHEEL @ONLY)
install (DIRECTORY DESTINATION ${OTROBOPT_PYTHON_MODULE_PATH}/${PACKAGE_NAME}-${PACKAGE_VERSION}.dist-info)
install (FILES ${CMAKE_CURRENT_BINARY_DIR}/METADATA ${CMAKE_CURRENT_BINARY_DIR}/WHEEL
install (FILES ${CMAKE_CURRENT_BINARY_DIR}/METADATA
DESTINATION ${OTROBOPT_PYTHON_MODULE_PATH}/${PACKAGE_NAME}-${PACKAGE_VERSION}.dist-info)
2 changes: 1 addition & 1 deletion python/src/InverseFORM_doc.i.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ Parameters
----------
event : :class:`~openturns.RandomVector`
Event we are computing the probability of.
parameterName : string
parameterName : str
The optimized parameter
physicalStartingPoint : sequence of float
Optimization starting point
Expand Down
4 changes: 0 additions & 4 deletions python/src/WHEEL.in
View file

This file was deleted.

0 comments on commit 38db2f5

Please sign in to comment.