1 // -*- C++ -*-
2 /**
3 * @brief The test file of FunctionalChaosAlgorithm 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
sobol(const Indices & indices,const Point & a)27 Scalar sobol(const Indices & indices,
28 const Point & a)
29 {
30 Scalar value = 1.0;
31 for (UnsignedInteger i = 0; i < indices.getSize(); ++i)
32 {
33 value *= 1.0 / (3.0 * pow(1.0 + a[indices[i]], 2.0));
34 }
35 return value;
36 }
37
main(int,char * [])38 int main(int, char *[])
39 {
40 TESTPREAMBLE;
41 OStream fullprint(std::cout);
42 setRandomGenerator();
43
44 try
45 {
46
47 // Problem parameters
48 UnsignedInteger dimension = 5;
49
50 // Reference analytical values
51 Scalar meanTh = 1.0;
52 Scalar covTh = 1.0;
53 Point a(dimension);
54 // Create the gSobol function
55 Description inputVariables(dimension);
56 Description formula(1);
57 formula[0] = "1.0";
58 for (UnsignedInteger i = 0; i < dimension; ++i)
59 {
60 a[i] = 0.5 * i;
61 covTh *= 1.0 + 1.0 / (3.0 * pow(1.0 + a[i], 2.0));
62 inputVariables[i] = (OSS() << "xi" << i);
63 formula[0] = (OSS() << formula[0] << " * ((abs(4.0 * xi" << i << " - 2.0) + " << a[i] << ") / (1.0 + " << a[i] << "))");
64 }
65 --covTh;
66
67 SymbolicFunction model(inputVariables, formula);
68
69 // Create the input distribution
70 Collection<Distribution> marginals(dimension);
71 for (UnsignedInteger i = 0; i < dimension; ++i)
72 {
73 marginals[i] = Uniform(0.0, 1.0);
74 }
75 ComposedDistribution distribution(marginals);
76
77 // Create the orthogonal basis
78 Collection<OrthogonalUniVariatePolynomialFamily> polynomialCollection(dimension);
79 for (UnsignedInteger i = 0; i < dimension; ++i)
80 {
81 polynomialCollection[i] = LegendreFactory();
82 }
83 LinearEnumerateFunction enumerateFunction(dimension);
84 OrthogonalProductPolynomialFactory productBasis(polynomialCollection, enumerateFunction);
85
86 // Create the adaptive strategy
87 // We can choose amongst several strategies
88 // First, the most efficient (but more complex!) strategy
89 Collection<AdaptiveStrategy> listAdaptiveStrategy(0);
90 UnsignedInteger degree = 4;
91 UnsignedInteger indexMax = enumerateFunction.getStrataCumulatedCardinal(degree);
92 UnsignedInteger basisDimension = enumerateFunction.getStrataCumulatedCardinal(degree / 2);
93 Scalar threshold = 1.0e-6;
94 listAdaptiveStrategy.add(CleaningStrategy(productBasis, indexMax, basisDimension, threshold, false));
95 // Second, the most used (and most basic!) strategy
96 listAdaptiveStrategy.add(FixedStrategy(productBasis, enumerateFunction.getStrataCumulatedCardinal(degree)));
97 // Third, a slight enhancement with respect to the basic strategy
98 listAdaptiveStrategy.add(SequentialStrategy(productBasis, enumerateFunction.getStrataCumulatedCardinal(degree / 2), false));
99
100 for(UnsignedInteger adaptiveStrategyIndex = 0; adaptiveStrategyIndex < listAdaptiveStrategy.getSize(); ++adaptiveStrategyIndex)
101 {
102 AdaptiveStrategy adaptiveStrategy(listAdaptiveStrategy[adaptiveStrategyIndex]);
103 // Create the projection strategy
104 UnsignedInteger samplingSize = 250;
105 Collection<ProjectionStrategy> listProjectionStrategy(0);
106 // The least squares strategy
107 // Monte Carlo sampling
108 listProjectionStrategy.add(LeastSquaresStrategy(MonteCarloExperiment(samplingSize)));
109 // LHS sampling
110 listProjectionStrategy.add(LeastSquaresStrategy(LHSExperiment(samplingSize)));
111 // Low Discrepancy sequence
112 listProjectionStrategy.add(LeastSquaresStrategy(LowDiscrepancyExperiment(LowDiscrepancySequence(SobolSequence()), samplingSize)));
113 // The integration strategy
114 // Monte Carlo sampling
115 listProjectionStrategy.add(IntegrationStrategy(MonteCarloExperiment(samplingSize)));
116 // LHS sampling
117 listProjectionStrategy.add(IntegrationStrategy(LHSExperiment(samplingSize)));
118 // Low Discrepancy sequence
119 listProjectionStrategy.add(IntegrationStrategy(LowDiscrepancyExperiment(LowDiscrepancySequence(SobolSequence()), samplingSize)));
120 for(UnsignedInteger projectionStrategyIndex = 0; projectionStrategyIndex < listProjectionStrategy.getSize(); ++projectionStrategyIndex)
121 {
122 ProjectionStrategy projectionStrategy(listProjectionStrategy[projectionStrategyIndex]);
123 // Create the polynomial chaos algorithm
124 Scalar maximumResidual = 1.0e-10;
125 FunctionalChaosAlgorithm algo(model, distribution, adaptiveStrategy, projectionStrategy);
126 algo.setMaximumResidual(maximumResidual);
127 // Reinitialize the RandomGenerator to see the effect of the sampling method only
128 RandomGenerator::SetSeed(0);
129
130 algo.run();
131
132 // Examine the results
133 FunctionalChaosResult result(algo.getResult());
134 fullprint << "//////////////////////////////////////////////////////////////////////" << std::endl;
135 fullprint << algo.getAdaptiveStrategy() << std::endl;
136 fullprint << algo.getProjectionStrategy() << std::endl;
137 Point residuals(result.getResiduals());
138 fullprint << "residuals=" << std::fixed << std::setprecision(5) << residuals << std::endl;
139 Point relativeErrors(result.getRelativeErrors());
140 fullprint << "relative errors=" << std::fixed << std::setprecision(5) << relativeErrors << std::endl;
141
142 // Post-process the results
143 FunctionalChaosRandomVector vector(result);
144 Scalar mean = vector.getMean()[0];
145 fullprint << "mean=" << std::fixed << std::setprecision(5) << mean << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(mean - meanTh) << std::endl;
146 Scalar variance = vector.getCovariance()(0, 0);
147 fullprint << "variance=" << std::fixed << std::setprecision(5) << variance << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(variance - covTh) << std::endl;
148 FunctionalChaosSobolIndices sensitivity(result);
149 Indices indices(1);
150 for(UnsignedInteger i = 0; i < dimension; ++i)
151 {
152 indices[0] = i;
153 Scalar value = sensitivity.getSobolIndex(i);
154 fullprint << "Sobol index " << i << " = " << std::fixed << std::setprecision(5) << value << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(value - sobol(indices, a) / covTh) << std::endl;
155 }
156 indices = Indices(2);
157 UnsignedInteger k = 0;
158 for (UnsignedInteger i = 0; i < dimension; ++i)
159 {
160 indices[0] = i;
161 for (UnsignedInteger j = i + 1; j < dimension; ++j)
162 {
163 indices[1] = j;
164 Scalar value = sensitivity.getSobolIndex(indices);
165 fullprint << "Sobol index " << indices << " =" << std::fixed << std::setprecision(5) << value << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(value - sobol(indices, a) / covTh) << std::endl;
166 k = k + 1;
167 }
168 }
169 indices = Indices(3);
170 indices[0] = 0;
171 indices[1] = 1;
172 indices[2] = 2;
173 Scalar value = sensitivity.getSobolIndex(indices);
174 fullprint << "Sobol index " << indices << " =" << std::fixed << std::setprecision(5) << value << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(value - sobol(indices, a) / covTh) << std::endl;
175 // First order grouped indice
176 indices = Indices(2);
177 indices[0] = 0;
178 indices[1] = 1;
179 Scalar exactS = sobol(indices, a) / covTh;
180 value = sensitivity.getSobolGroupedIndex(indices);
181 fullprint << "Grouped First Sobol index " << indices << " =" << std::fixed << std::setprecision(5) << value << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(value - exactS) << std::endl;
182 // Total order grouped indice
183 indices = Indices(2);
184 indices[0] = 0;
185 indices[1] = 1;
186 value = sensitivity.getSobolGroupedTotalIndex(indices);
187 Indices complementaryindices = Indices(3);
188 complementaryindices[0] = 2;
189 complementaryindices[1] = 3;
190 complementaryindices[2] = 4;
191 exactS = 1 - sensitivity.getSobolGroupedIndex(complementaryindices);
192 fullprint << "Grouped Total Sobol index " << indices << " =" << std::fixed << std::setprecision(5) << value << " absolute error=" << std::scientific << std::setprecision(1) << std::abs(value - exactS) << std::endl;
193 }
194 }
195 }
196 catch (TestFailed & ex)
197 {
198 std::cerr << ex << std::endl;
199 return ExitCode::Error;
200 }
201
202
203 return ExitCode::Success;
204 }
205