1 // -*- C++ -*-
2 /**
3 * @brief The test file of class RandomMixture for standard methods
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
main(int,char * [])27 int main(int, char *[])
28 {
29 TESTPREAMBLE;
30 OStream fullprint(std::cout);
31 setRandomGenerator();
32
33 ResourceMap::SetAsUnsignedInteger("RandomMixture-DefaultMaxSize", 4000000);
34
35 try
36 {
37 // Create a collection of test-cases and the associated references
38 UnsignedInteger numberOfTests = 3;
39 Collection< Collection< Distribution > > testCases(numberOfTests);
40 Collection< Distribution > references(numberOfTests);
41 testCases[0] = Collection<Distribution>(2);
42 testCases[0][0] = Uniform(-1.0, 3.0);
43 testCases[0][1] = Uniform(-1.0, 3.0);
44 references[0] = Triangular(-2.0, 2.0, 6.0);
45 testCases[1] = Collection<Distribution>(3);
46 testCases[1][0] = Normal();
47 testCases[1][1] = Normal(1.0, 2.0);
48 testCases[1][2] = Normal(-2.0, 2.0);
49 references[1] = Normal(-1.0, 3.0);
50 testCases[2] = Collection<Distribution>(3);
51 testCases[2][0] = Exponential();
52 testCases[2][1] = Exponential();
53 testCases[2][2] = Exponential();
54 references[2] = Gamma(3.0, 1.0, 0.0);
55 fullprint << "testCases=" << testCases << std::endl;
56 fullprint << "references=" << references << std::endl;
57 for (UnsignedInteger testIndex = 0; testIndex < testCases.getSize(); ++testIndex)
58 {
59 // Instanciate one distribution object
60 RandomMixture distribution(testCases[testIndex]);
61 distribution.setBlockMin(5);
62 distribution.setBlockMax(20);
63 Distribution distributionReference(references[testIndex]);
64 fullprint << "Distribution " << distribution << std::endl;
65 std::cout << "Distribution " << distribution << std::endl;
66
67 // Is this distribution elliptical ?
68 fullprint << "Elliptical = " << (distribution.isElliptical() ? "true" : "false") << std::endl;
69
70 // Is this distribution continuous ?
71 fullprint << "Continuous = " << (distribution.isContinuous() ? "true" : "false") << std::endl;
72
73 // Test for realization of distribution
74 Point oneRealization = distribution.getRealization();
75 fullprint << "oneRealization=" << oneRealization << std::endl;
76
77 // Test for sampling
78 UnsignedInteger size = 10000;
79 Sample oneSample = distribution.getSample( size );
80 fullprint << "oneSample first=" << oneSample[0] << " last=" << oneSample[size - 1] << std::endl;
81 fullprint << "mean=" << oneSample.computeMean() << std::endl;
82 fullprint << "covariance=" << oneSample.computeCovariance() << std::endl;
83
84 size = 100;
85 for (UnsignedInteger i = 0; i < 2; ++i)
86 {
87 fullprint << "Kolmogorov test for the generator, sample size=" << size << " is " << (FittingTest::Kolmogorov(distribution.getSample(size), distribution).getBinaryQualityMeasure() ? "accepted" : "rejected") << std::endl;
88 size *= 10;
89 }
90
91 // Define a point
92 Point point(distribution.getDimension(), 0.5);
93 fullprint << "Point= " << point << std::endl;
94
95 // Show PDF and CDF of point
96 Scalar eps = 1e-5;
97 UnsignedInteger oldPrecision = PlatformInfo::GetNumericalPrecision();
98 PlatformInfo::SetNumericalPrecision(5);
99 Point DDF = distribution.computeDDF(point);
100 fullprint << "ddf =" << DDF << std::endl;
101 fullprint << "ddf (ref)=" << distributionReference.computeDDF(point) << std::endl;
102 PlatformInfo::SetNumericalPrecision(oldPrecision);
103 Scalar PDF = distribution.computePDF(point);
104 fullprint << "pdf =" << PDF << std::endl;
105 fullprint << "pdf (FD)=" << (distribution.computeCDF( point + Point(1, eps) ) - distribution.computeCDF( point + Point(1, -eps) )) / (2.0 * eps) << std::endl;
106 fullprint << "pdf (ref)=" << distributionReference.computePDF(point) << std::endl;
107 Scalar CDF = distribution.computeCDF( point );
108 fullprint << "cdf =" << CDF << std::endl;
109 fullprint << "cdf (ref)=" << distributionReference.computeCDF(point) << std::endl;
110 Complex CF = distribution.computeCharacteristicFunction( point[0] );
111 fullprint << "characteristic function=" << CF << std::endl;
112 Complex LCF = distribution.computeLogCharacteristicFunction( point[0] );
113 fullprint << "log characteristic function=" << LCF << std::endl;
114 Point quantile = distribution.computeQuantile( 0.95 );
115 fullprint << "quantile =" << quantile << std::endl;
116 fullprint << "quantile (ref)=" << distributionReference.computeQuantile(0.95) << std::endl;
117 fullprint << "cdf(quantile)=" << distribution.computeCDF(quantile) << std::endl;
118 Point quantileComp = distribution.computeQuantile( 0.95, true );
119 fullprint << "quantile comp.=" << quantileComp << std::endl;
120 fullprint << "cdfComp(quantileComp)=" << distribution.computeComplementaryCDF(quantileComp) << std::endl;
121 fullprint << "entropy =" << distribution.computeEntropy() << std::endl;
122 fullprint << "entropy (ref)=" << distributionReference.computeEntropy() << std::endl;
123 fullprint << "entropy (MC)=" << -distribution.computeLogPDF(distribution.getSample(1000000)).computeMean()[0] << std::endl;
124 Point mean = distribution.getMean();
125 fullprint << "mean =" << mean << std::endl;
126 fullprint << "mean (ref)=" << distributionReference.getMean() << std::endl;
127 Point standardDeviation = distribution.getStandardDeviation();
128 fullprint << "standard deviation =" << standardDeviation << std::endl;
129 fullprint << "standard deviation (ref)=" << distributionReference.getStandardDeviation() << std::endl;
130 Point skewness = distribution.getSkewness();
131 fullprint << "skewness =" << skewness << std::endl;
132 fullprint << "skewness (ref)=" << distributionReference.getSkewness() << std::endl;
133 Point kurtosis = distribution.getKurtosis();
134 fullprint << "kurtosis =" << kurtosis << std::endl;
135 fullprint << "kurtosis (ref)=" << distributionReference.getKurtosis() << std::endl;
136 CovarianceMatrix covariance = distribution.getCovariance();
137 fullprint << "covariance =" << covariance << std::endl;
138 fullprint << "covariance (ref)=" << distributionReference.getCovariance() << std::endl;
139 RandomMixture::PointWithDescriptionCollection parameters = distribution.getParametersCollection();
140 fullprint << "parameters=" << parameters << std::endl;
141 /* distribution.setIntegrationNodesNumber(6);
142 for (UnsignedInteger i = 0; i < 6; ++i) fullprint << "standard moment n=" << i << ", value=" << distribution.getStandardMoment(i) << std::endl;*/
143 fullprint << "Standard representative=" << distribution.getStandardRepresentative().__str__() << std::endl;
144 fullprint << "blockMin=" << distribution.getBlockMin() << std::endl;
145 fullprint << "blockMax=" << distribution.getBlockMax() << std::endl;
146 fullprint << "maxSize=" << distribution.getMaxSize() << std::endl;
147 fullprint << "alpha=" << distribution.getAlpha() << std::endl;
148 fullprint << "beta=" << distribution.getBeta() << std::endl;
149 }
150 // Tests of the simplification mechanism
151 Point weights(0);
152 Collection< Distribution > coll(0);
153 // coll.add(Dirac(0.5));
154 // weights.add(1.0);
155 coll.add(Normal(1.0, 2.0));
156 weights.add(1.0);
157 coll.add(Normal(2.0, 1.0));
158 weights.add(1.0);
159 coll.add(Uniform(-2.0, 2.0));
160 weights.add(1.0);
161 coll.add(Exponential(2.0, 0.0));
162 weights.add(1.0);
163 // RandomMixture rm(coll, weights);
164 // coll.add(rm);
165 // weights.add(-2.5);
166 coll.add(Gamma(3.0, 4.0, 0.0));
167 weights.add(1.0);
168 RandomMixture distribution(coll, weights);
169 fullprint << "distribution=" << distribution << std::endl;
170 fullprint << "distribution=" << distribution.__str__() << std::endl;
171 for (UnsignedInteger i = 0; i < 30; ++i)
172 {
173 Scalar x = -12.0 + i;
174 fullprint << "pdf(" << x << ")=" << distribution.computePDF(x) << std::endl;
175 }
176 Graph graph(distribution.drawPDF());
177 Distribution ks(KernelSmoothing().build(distribution.getSample(1000000)));
178 graph.add(ks.drawPDF());
179 Description colors(2);
180 colors[0] = "red";
181 colors[1] = "green";
182 graph.setColors(colors);
183 // graph.draw("validation");
184 // Test for the projection
185 Collection<DistributionFactory> collFactories(0);
186 collFactories.add(UniformFactory());
187 collFactories.add(NormalFactory());
188 collFactories.add(TriangularFactory());
189 collFactories.add(ExponentialFactory());
190 collFactories.add(GammaFactory());
191 // collFactories.add(TrapezoidalFactory());
192 Point norms;
193 Collection< Distribution > result(distribution.project(collFactories, norms));
194 fullprint << "projections=" << result << std::endl;
195 fullprint << "norms=" << norms << std::endl;
196 //------------------------------ Multivariate tests ------------------------------//
197 // 2D RandomMixture
198 Collection< Distribution > collection(0);
199 collection.add(Normal(0.0, 1.0));
200 collection.add(Normal(0.0, 1.0));
201 collection.add(Normal(0.0, 1.0));
202
203 Matrix weightMatrix(2, 3);
204 weightMatrix(0, 0) = 1.0;
205 weightMatrix(0, 1) = -2.0;
206 weightMatrix(0, 2) = 1.0;
207 weightMatrix(1, 0) = 1.0;
208 weightMatrix(1, 1) = 1.0;
209 weightMatrix(1, 2) = -3.0;
210
211 // Build the RandomMixture
212 RandomMixture distribution2D(collection, weightMatrix);
213 fullprint << "distribution = " << distribution2D << std::endl;
214 fullprint << "range = " << distribution2D.getRange() << std::endl;
215 fullprint << "mean = " << distribution2D.getMean() << std::endl;
216 fullprint << "cov = " << distribution2D.getCovariance() << std::endl;
217 fullprint << "sigma = " << distribution2D.getStandardDeviation() << std::endl;
218 fullprint << "entropy =" << distribution.computeEntropy() << std::endl;
219 distribution2D.setBlockMin(3);
220 distribution2D.setBlockMax(10);
221
222 // Build a grid for validation
223 const Scalar xMin = distribution2D.getRange().getLowerBound()[0];
224 const Scalar xMax = distribution2D.getRange().getUpperBound()[0];
225 const Scalar yMin = distribution2D.getRange().getLowerBound()[1];
226 const Scalar yMax = distribution2D.getRange().getUpperBound()[1];
227 // Number of points of discretization
228 const UnsignedInteger nx = 4;
229 const UnsignedInteger ny = 4;
230 Point boxParameters(2);
231 boxParameters[0] = nx;
232 boxParameters[1] = ny;
233 Box boxGrid(boxParameters);
234 Sample grid(boxGrid.generate());
235 // scaling box grid
236 Point scaleFactor(2);
237 scaleFactor[0] = 0.25 * (xMax - xMin);
238 scaleFactor[1] = 0.25 * (yMax - yMin);
239 grid *= scaleFactor;
240 //translating
241 Point translateFactor(2);
242 translateFactor[0] = distribution2D.getMean()[0];
243 translateFactor[1] = distribution2D.getMean()[1];
244 grid += translateFactor;
245 // Compute PDF
246 // parameters for theoritical PDF, obtained thanks to Maple
247 const Scalar factor = sqrt(2.0) / (20 * M_PI);
248 for (UnsignedInteger index = 0; index < grid.getSize(); ++ index)
249 {
250 const Point point(grid[index]);
251 const Scalar PDF = distribution2D.computePDF(point);
252 // Very small values are not very accurate on x86, skip them
253 if (PDF < 1.e-12) continue;
254 fullprint << "pdf =" << PDF << std::endl;
255 const Scalar y = point[1];
256 const Scalar x = point[0];
257 fullprint << "pdf (ref)=" << factor * exp(-3.0 / 50.0 * y * y - 2.0 / 25 * x * y - 11.0 / 100 * x * x) << std::endl;
258 }
259 // 2D test, but too much CPU consuming
260
261 Collection< Distribution > collUniforme(0);
262 collUniforme.add(Uniform(0, 1));
263 collUniforme.add(Uniform(0, 1));
264 collUniforme.add(Uniform(0, 1));
265 // Build the RandomMixture
266 RandomMixture dist_2D(collUniforme, weightMatrix);
267 dist_2D.setBlockMin(3);
268 dist_2D.setBlockMax(8);
269
270 fullprint << "new distribution = " << dist_2D << std::endl;
271 fullprint << "range = " << dist_2D.getRange() << std::endl;
272 fullprint << "mean = " << dist_2D.getMean() << std::endl;
273 fullprint << "cov = " << dist_2D.getCovariance() << std::endl;
274 fullprint << "sigma = " << dist_2D.getStandardDeviation() << std::endl;
275
276 // Discretization on 2D grid [mu, mu+\sigma]
277 Sample newGrid(boxGrid.generate());
278 // scaling box grid
279 newGrid *= dist_2D.getStandardDeviation();
280 //translating
281 newGrid += dist_2D.getMean();
282 // Compute PDF
283 for (UnsignedInteger index = 0; index < newGrid.getSize(); ++ index)
284 {
285 const Point point(newGrid[index]);
286 const Scalar PDF = dist_2D.computePDF(point);
287 fullprint << "pdf =" << PDF << std::endl;
288 }
289 // 3D test
290 ResourceMap::SetAsUnsignedInteger("RandomMixture-DefaultMaxSize", 8290688);
291 Collection<Distribution> collectionMixture(0);
292 collectionMixture.add(Normal(2.0, 1.0));
293 collectionMixture.add(Normal(-2.0, 1.0));
294 const Mixture mixture(collectionMixture);
295 Collection<Distribution> collection3D(0);
296 collection3D.add(Normal(0.0, 1.0));
297 collection3D.add(mixture);
298 collection3D.add(Uniform(0, 1));
299 collection3D.add(Uniform(0, 1));
300 // Set weights
301 weightMatrix = Matrix(3, 4);
302 weightMatrix(0, 0) = 1.0;
303 weightMatrix(0, 1) = -0.05;
304 weightMatrix(0, 2) = 1.0;
305 weightMatrix(0, 3) = -0.5;
306
307 weightMatrix(1, 0) = 0.5;
308 weightMatrix(1, 1) = 1.0;
309 weightMatrix(1, 2) = -0.05;
310 weightMatrix(1, 3) = 0.3;
311
312 weightMatrix(2, 0) = -0.5;
313 weightMatrix(2, 1) = -0.1;
314 weightMatrix(2, 2) = 1.2;
315 weightMatrix(2, 3) = -0.8;
316
317 RandomMixture dist_3D(collection3D, weightMatrix);
318 dist_3D.setBlockMin(3);
319 dist_3D.setBlockMax(6);
320
321 fullprint << "3D distribution = " << dist_3D << std::endl;
322 fullprint << "range = " << dist_3D.getRange() << std::endl;
323 fullprint << "mean = " << dist_3D.getMean() << std::endl;
324 fullprint << "cov = " << dist_3D.getCovariance() << std::endl;
325 fullprint << "sigma = " << dist_3D.getStandardDeviation() << std::endl;
326 fullprint << "entropy =" << distribution.computeEntropy() << std::endl;
327
328 // Total number of points (is (2+2)**3)
329 // Test is CPU consuming
330 // Number of points of discretization
331 const UnsignedInteger N = 2;
332 Point box3DParameters(3, N);
333 Box box3D(box3DParameters);
334 // Grid ==> (mu, mu+sigma)
335 Sample grid3D(box3D.generate());
336 //scaling
337 grid3D *= dist_3D.getStandardDeviation();
338 // translating
339 grid3D += dist_3D.getMean();
340 for (UnsignedInteger index = 0; index < grid3D.getSize() / 4; ++ index)
341 {
342 const Point point(grid3D[4 * index]);
343 const Scalar PDF = dist_3D.computePDF(point);
344 fullprint << "pdf =" << PDF << std::endl;
345 }
346 // Test for ticket 882 (only one Dirac)
347 // The segfault was triggered during the construction...
348 RandomMixture mixture2(Collection<Distribution>(1, Dirac()));
349 // After what it was impossible to draw the PDF or the CDF due to a lack of support computation
350 Graph graphPDF(mixture2.drawPDF());
351 Graph graphCDF(mixture2.drawCDF());
352 // Test computeQuantile for the specific case of an analytical 1D mixture
353 RandomMixture case1(Collection<Distribution>(1, ChiSquare()), Point(1, 0.1));
354 Scalar q = case1.computeQuantile(0.95)[0];
355 fullprint << "case 1, q=" << q << std::endl;
356 q = case1.computeQuantile(0.95, true)[0];
357 fullprint << "case 1, q comp=" << q << std::endl;
358 RandomMixture case2(Collection<Distribution>(1, ChiSquare()), Point(1, -0.1));
359 q = case2.computeQuantile(0.95)[0];
360 fullprint << "case 2, q=" << q << std::endl;
361 q = case2.computeQuantile(0.95, true)[0];
362 fullprint << "case 2, q comp=" << q << std::endl;
363 // For ticket 953
364 {
365 TruncatedDistribution atom1(Uniform(0.0, 1.0), 0.0, 1.0);
366 Uniform atom2(0.0, 2.0);
367 Distribution sum(atom1 + atom2);
368 fullprint << "sum=" << sum << std::endl;
369 fullprint << "CDF=" << sum.computeCDF(2.0) << std::endl;
370 fullprint << "quantile=" << sum.computeQuantile(0.2) << std::endl;
371 }
372 {
373 Scalar minS = 0.2;
374 Scalar maxS = 10.0;
375 Scalar muS = (std::log(minS) + std::log(maxS)) / 2.0;
376 Scalar sigma = (std::log(maxS) - muS) / 3.0;
377 TruncatedDistribution atom1(LogNormal(muS, sigma), minS, maxS);
378 Uniform atom2(0.0, 2.0);
379 Distribution sum(atom1 + atom2);
380 fullprint << "sum=" << sum << std::endl;
381 fullprint << "CDF=" << sum.computeCDF(2.0) << std::endl;
382 fullprint << "quantile=" << sum.computeQuantile(0.2) << std::endl;
383 }
384 // For ticket 1129
385 {
386 distribution = RandomMixture(Collection<Distribution>(200, Uniform()));
387 fullprint << "CDF(0)=" << distribution.computeCDF(0.0) << std::endl;
388 }
389 }
390 catch (TestFailed & ex)
391 {
392 std::cerr << ex << std::endl;
393 return ExitCode::Error;
394 }
395
396 return ExitCode::Success;
397 }
398