1 // -*- C++ -*-
2 /**
3 * @brief Result implements the results obtained from the First Order Reliability Method
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/FORM.hxx"
22 #include "openturns/Distribution.hxx"
23 #include "openturns/PersistentCollection.hxx"
24 #include "openturns/Point.hxx"
25 #include "openturns/PersistentObjectFactory.hxx"
26
27 BEGIN_NAMESPACE_OPENTURNS
28
29
30
31 CLASSNAMEINIT(FORMResult)
32
33 static const Factory<FORMResult> Factory_FORMResult;
34
35 typedef PersistentCollection<PointWithDescription> PersistentSensitivity;
36
37 /*
38 * @brief Standard constructor: the class is defined by an optimisation algorithm, a failure event and a physical starting point
39 */
FORMResult(const Point & standardSpaceDesignPoint,const RandomVector & limitStateVariable,const Bool isStandardPointOriginInFailureSpace)40 FORMResult::FORMResult(const Point & standardSpaceDesignPoint,
41 const RandomVector & limitStateVariable,
42 const Bool isStandardPointOriginInFailureSpace):
43 AnalyticalResult(standardSpaceDesignPoint, limitStateVariable, isStandardPointOriginInFailureSpace),
44 eventProbability_(0.),
45 generalisedReliabilityIndex_(0.),
46 eventProbabilitySensitivity_(Sensitivity(0)),
47 isAlreadyComputedEventProbabilitySensitivity_(false)
48 {
49 computeEventProbability();
50 computeGeneralisedReliabilityIndex();
51 }
52
53 /* Default constructor */
FORMResult()54 FORMResult::FORMResult():
55 AnalyticalResult(),
56 eventProbability_(0.),
57 generalisedReliabilityIndex_(0.),
58 eventProbabilitySensitivity_(Sensitivity(0)),
59 isAlreadyComputedEventProbabilitySensitivity_(false)
60 {
61 // Nothing to do
62 }
63
64 /* Virtual constructor */
clone() const65 FORMResult * FORMResult::clone() const
66 {
67 return new FORMResult(*this);
68 }
69
70 /* The function that actually evaluates the event probability with FORM approximation */
computeEventProbability() const71 void FORMResult::computeEventProbability() const
72 {
73 /* evaluate the event probability */
74 /* Be careful! computeCDF method takes an Point as an input argument */
75 /* in the standard space all marginals of the standard distribution are identical */
76 eventProbability_ = getLimitStateVariable().getImplementation()->getAntecedent().getDistribution().getStandardDistribution().getMarginal(0).computeCDF(Point(1, -getHasoferReliabilityIndex()));
77
78 if (getIsStandardPointOriginInFailureSpace())
79 {
80 // isStandardPointOriginInFailureSpace is true: unusual case
81 eventProbability_ = 1.0 - eventProbability_;
82 }
83 }
84
85 /* EventProbability accessor */
getEventProbability() const86 Scalar FORMResult::getEventProbability() const
87 {
88 return eventProbability_;
89 }
90
91 /* EventProbability accessor */
setEventProbability(const Scalar & eventProbability)92 void FORMResult::setEventProbability(const Scalar & eventProbability)
93 {
94 eventProbability_ = eventProbability;
95 }
96
97 /* The function that actually evaluates the generalised reliability index with FORM approximation */
computeGeneralisedReliabilityIndex() const98 void FORMResult::computeGeneralisedReliabilityIndex() const
99 {
100 /* GeneralisedReliabilityIndex is defined by : - Inverse standard marginal CDF (eventProbability). It
101 will thus be negative if the eventProbability is > 0.5. */
102 generalisedReliabilityIndex_ = -getLimitStateVariable().getImplementation()->getAntecedent().getDistribution().getStandardDistribution().getMarginal(0).computeQuantile(eventProbability_)[0];
103 }
104
105
106 /* GeneralisedReliabilityIndex accessor */
getGeneralisedReliabilityIndex() const107 Scalar FORMResult::getGeneralisedReliabilityIndex() const
108 {
109 return generalisedReliabilityIndex_;
110 }
111
112 /* GeneralisedReliabilityIndex accessor */
setGeneralisedReliabilityIndex(const Scalar & generalisedReliabilityIndex)113 void FORMResult::setGeneralisedReliabilityIndex(const Scalar & generalisedReliabilityIndex)
114 {
115 generalisedReliabilityIndex_ = generalisedReliabilityIndex;
116 }
117
118 /* The function that actually evaluates the event probability sensitivity with FORM approximation */
computeEventProbabilitySensitivity() const119 void FORMResult::computeEventProbabilitySensitivity() const
120 {
121 Point correctedReliabilityIndex(1, (getIsStandardPointOriginInFailureSpace() ? getHasoferReliabilityIndex() : -getHasoferReliabilityIndex()));
122 Distribution antecedent(getLimitStateVariable().getImplementation()->getAntecedent().getDistribution());
123 UnsignedInteger dimension = antecedent.getDimension();
124
125 /* Be carefull! computeCDF method takes an Point as an input argument */
126 /* in the standard space all marginals of the standard distribution are identical */
127 /* evaluate one marginal at the reliability index : the marginal is symmetric with respect to zero */
128 const Distribution standardMarginalDistribution(antecedent.getStandardDistribution().getMarginal(0));
129 Scalar correctedReliabilityIndexDensity = standardMarginalDistribution.computePDF(correctedReliabilityIndex);
130
131 if (! getIsStandardPointOriginInFailureSpace())
132 {
133 // isStandardPointOriginInFailureSpace is false : usual case
134 correctedReliabilityIndexDensity *= -1.0;
135 }
136 /* We initialize eventProbabilitySensitivity_ this way in order to inherit from the name and description of hasoferReliabilityIndexSensitivity */
137 eventProbabilitySensitivity_ = getHasoferReliabilityIndexSensitivity();
138 /* sensitivity with respect to the parameters influencing beta (beta is not sensitive to the parameters relative to the copula type) */
139 UnsignedInteger size = eventProbabilitySensitivity_.getSize();
140 for(UnsignedInteger i = 0; i < size; ++i) eventProbabilitySensitivity_[i] *= correctedReliabilityIndexDensity;
141 /* sensitivity with respect to the parameters of the density generator of the standard distribution */
142 PointWithDescriptionCollection standardMarginalParametersCollection(standardMarginalDistribution.getParametersCollection());
143 UnsignedInteger parametersDimension = standardMarginalParametersCollection[0].getDimension();
144 /* there 2 parameters generic to 1D elliptical distributions : mean and standard deviation */
145 UnsignedInteger genericParametersNumber = 2;
146 if (antecedent.getImplementation()->hasEllipticalCopula() && parametersDimension > genericParametersNumber)
147 {
148 const Point standardParametersGradient(standardMarginalDistribution.computeCDFGradient(correctedReliabilityIndex));
149 /* shift is the number of parameters of the correlation matrix (upper triangle) for elliptical copula*/
150 const UnsignedInteger shift = dimension * (dimension - 1) / 2;
151 for (UnsignedInteger index = genericParametersNumber; index < standardParametersGradient.getDimension(); ++index) eventProbabilitySensitivity_[dimension][index + shift - genericParametersNumber] = standardParametersGradient[index];
152 }
153 isAlreadyComputedEventProbabilitySensitivity_ = true;
154
155 }// end computeEventProbabilitySensitivity()
156
157
158 /* EventProbabilitySensitivity accessor */
getEventProbabilitySensitivity() const159 FORMResult::Sensitivity FORMResult::getEventProbabilitySensitivity() const
160 {
161 if (! isAlreadyComputedEventProbabilitySensitivity_) computeEventProbabilitySensitivity();
162 return eventProbabilitySensitivity_;
163 }
164
165
166 /* HasoferReliabilityIndexSensitivity Graph */
drawEventProbabilitySensitivity(Scalar width) const167 AnalyticalResult::GraphCollection FORMResult::drawEventProbabilitySensitivity(Scalar width) const
168 {
169 GraphCollection eventProbabilitySensitivityGraphCollection(0);
170 // To ensure that the eventProbabilitySensitivity_ are up to date
171 if (! isAlreadyComputedEventProbabilitySensitivity_) computeEventProbabilitySensitivity();
172 UnsignedInteger dimension = getStandardSpaceDesignPoint().getDimension();
173 UnsignedInteger size = eventProbabilitySensitivity_.getSize();
174
175 // The first graph shows the sensitivities with respect to the marginal parameters
176 Sensitivity marginalSensitivity(dimension);
177 for(UnsignedInteger i = 0; i < dimension; ++i) marginalSensitivity[i] = eventProbabilitySensitivity_[i];
178 Graph eventProbabilitySensitivityGraphMarginal(drawSensitivity(marginalSensitivity, width));
179 OSS oss1;
180 oss1 << "FORM - Event Probability Sensitivities - Marginal parameters - " << getLimitStateVariable().getName() ;
181 eventProbabilitySensitivityGraphMarginal.setTitle(oss1);
182 eventProbabilitySensitivityGraphCollection.add(eventProbabilitySensitivityGraphMarginal);
183
184 // The second graph shows the sensitivities with respect to the other parameters
185 if (size > dimension)
186 {
187 Sensitivity otherSensitivity(size - dimension);
188 for(UnsignedInteger i = dimension; i < size; ++i) otherSensitivity[i - dimension] = eventProbabilitySensitivity_[i];
189 Graph eventProbabilitySensitivityGraphOther(drawSensitivity(otherSensitivity, width));
190 OSS oss2;
191 oss2 << "FORM - Event Probability Sensitivities - Other parameters - " << getLimitStateVariable().getName() ;
192 eventProbabilitySensitivityGraphOther.setTitle(oss2);
193 eventProbabilitySensitivityGraphCollection.add(eventProbabilitySensitivityGraphOther);
194 }
195 return eventProbabilitySensitivityGraphCollection;
196 }
197
198
199 /* String converter */
__repr__() const200 String FORMResult::__repr__() const
201 {
202 OSS oss;
203 oss << "class=" << FORMResult::GetClassName()
204 << " " << AnalyticalResult::__repr__()
205 << " eventProbability=" << eventProbability_
206 << " generalisedReliabilityIndex=" << generalisedReliabilityIndex_
207 << " eventProbabilitySensitivity=" << eventProbabilitySensitivity_;
208 return oss;
209 }
210
211 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const212 void FORMResult::save(Advocate & adv) const
213 {
214 PersistentSensitivity sensitivity(eventProbabilitySensitivity_);
215 AnalyticalResult::save(adv);
216 adv.saveAttribute( "eventProbability_", eventProbability_ );
217 adv.saveAttribute( "generalisedReliabilityIndex_", generalisedReliabilityIndex_ );
218 adv.saveAttribute( "eventProbabilitySensitivity_", sensitivity );
219 adv.saveAttribute( "isAlreadyComputedEventProbabilitySensitivity_", isAlreadyComputedEventProbabilitySensitivity_ );
220 }
221
222 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)223 void FORMResult::load(Advocate & adv)
224 {
225 PersistentSensitivity sensitivity;
226 AnalyticalResult::load(adv);
227 adv.loadAttribute( "eventProbability_", eventProbability_ );
228 adv.loadAttribute( "generalisedReliabilityIndex_", generalisedReliabilityIndex_ );
229 adv.loadAttribute( "eventProbabilitySensitivity_", sensitivity );
230 adv.loadAttribute( "isAlreadyComputedEventProbabilitySensitivity_", isAlreadyComputedEventProbabilitySensitivity_ );
231 eventProbabilitySensitivity_ = sensitivity;
232 }
233
234 END_NAMESPACE_OPENTURNS
235
236