1 //                                               -*- C++ -*-
2 /**
3  *  @brief The test file of FORM class
4  *
5  *  Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6  *
7  *  This library is free software: you can redistribute it and/or modify
8  *  it under the terms of the GNU Lesser General Public License as published by
9  *  the Free Software Foundation, either version 3 of the License, or
10  *  (at your option) any later version.
11  *
12  *  This library is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU Lesser General Public License for more details.
16  *
17  *  You should have received a copy of the GNU Lesser General Public License
18  *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
19  *
20  */
21 #include "openturns/OT.hxx"
22 #include "openturns/OTtestcode.hxx"
23 
24 using namespace OT;
25 using namespace OT::Test;
26 
printPoint(const Point & point,const UnsignedInteger digits)27 String printPoint(const Point & point, const UnsignedInteger digits)
28 {
29   OSS oss;
30   oss << "[";
31   Scalar eps = pow(0.1, 1.0 * digits);
32   for (UnsignedInteger i = 0; i < point.getDimension(); i++)
33   {
34     oss << std::fixed << std::setprecision(digits) << (i == 0 ? "" : ",") << Bulk<double>((std::abs(point[i]) < eps) ? std::abs(point[i]) : point[i]);
35   }
36   oss << "]";
37   return oss;
38 }
39 
main(int,char * [])40 int main(int, char *[])
41 {
42   TESTPREAMBLE;
43   OStream fullprint(std::cout);
44 
45   try
46   {
47 
48     /* We create a numerical math function */
49     Description input(4);
50     input[0] = "E";
51     input[1] = "F";
52     input[2] = "L";
53     input[3] = "I";
54     SymbolicFunction myFunction(input, Description(1, "-F*L^3/(3*E*I)"));
55 
56     UnsignedInteger dim = myFunction.getInputDimension();
57     /* We create a normal distribution point of dimension 1 */
58     Point mean(dim, 0.0);
59     mean[0] = 50.0; // E
60     mean[1] =  1.0; // F
61     mean[2] = 10.0; // L
62     mean[3] =  5.0; // I
63     Point sigma(dim, 1.0);
64     IdentityMatrix R(dim);
65     Normal myDistribution(mean, sigma, R);
66 
67     /* We create a 'usual' RandomVector from the Distribution */
68     RandomVector vect(myDistribution);
69 
70     /* We create a composite random vector */
71     CompositeRandomVector output(myFunction, vect);
72 
73     /* We create an Event from this RandomVector */
74     ThresholdEvent myEvent(output, Less(), -3.0);
75 
76     /* We create a NearestPoint algorithm */
77     // Test function operator ()
78     input[0] = "x1";
79     input[1] = "x2";
80     input[2] = "x3";
81     input[3] = "x4";
82     SymbolicFunction levelFunction(input, Description(1, "x1+2*x2-3*x3+4*x4"));
83     Point startingPoint(4, 1.0);
84     AbdoRackwitz myAlgorithm(NearestPointProblem(levelFunction, 3.0));
85     myAlgorithm.setStartingPoint(startingPoint);
86     myAlgorithm.setMaximumIterationNumber(100);
87     myAlgorithm.setMaximumAbsoluteError(1.0e-10);
88     myAlgorithm.setMaximumRelativeError(1.0e-10);
89     myAlgorithm.setMaximumResidualError(1.0e-10);
90     myAlgorithm.setMaximumConstraintError(1.0e-10);
91 
92     /* We create a FORM algorithm */
93     /* The first parameter is an OptimizationAlgorithm */
94     /* The second parameter is an event */
95     /* The third parameter is a starting point for the design point research */
96     FORM myAlgo(myAlgorithm, myEvent, mean);
97 
98     fullprint << "FORM=" << myAlgo << std::endl;
99 
100     /* Perform the simulation */
101     myAlgo.run();
102 
103     /* Stream out the result */
104     FORMResult result(myAlgo.getResult());
105     UnsignedInteger digits = 5;
106     fullprint << "event probability=" << result.getEventProbability() << std::endl;
107     fullprint << "generalized reliability index=" << std::setprecision(digits) << result.getGeneralisedReliabilityIndex() << std::endl;
108     fullprint << "standard space design point=" << printPoint(result.getStandardSpaceDesignPoint(), digits) << std::endl;
109     fullprint << "physical space design point=" << printPoint(result.getPhysicalSpaceDesignPoint(), digits) << std::endl;
110     fullprint << "is standard point origin in failure space? " << (result.getIsStandardPointOriginInFailureSpace() ? "true" : "false") << std::endl;
111     fullprint << "importance factors=" << printPoint(result.getImportanceFactors(), digits) << std::endl;
112     fullprint << "importance factors (classical)=" << printPoint(result.getImportanceFactors(AnalyticalResult::CLASSICAL), digits) << std::endl;
113     fullprint << "Hasofer reliability index=" << std::setprecision(digits) << result.getHasoferReliabilityIndex() << std::endl;
114   }
115   catch (TestFailed & ex)
116   {
117     std::cerr << ex << std::endl;
118     return ExitCode::Error;
119   }
120 
121 
122   return ExitCode::Success;
123 }
124