Visual Servoing Platform version 3.5.0
testRand.cpp
1/****************************************************************************
2 *
3 * ViSP, open source Visual Servoing Platform software.
4 * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5 *
6 * This software is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 * See the file LICENSE.txt at the root directory of this source
11 * distribution for additional information about the GNU GPL.
12 *
13 * For using ViSP with software that can not be combined with the GNU
14 * GPL, please contact Inria about acquiring a ViSP Professional
15 * Edition License.
16 *
17 * See http://visp.inria.fr for more information.
18 *
19 * This software was developed at:
20 * Inria Rennes - Bretagne Atlantique
21 * Campus Universitaire de Beaulieu
22 * 35042 Rennes Cedex
23 * France
24 *
25 * If you have questions regarding the use of this file, please contact
26 * Inria at visp@inria.fr
27 *
28 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30 *
31 * Description:
32 * Test pseudo random number generator.
33 *
34 *****************************************************************************/
35
36#include <visp3/core/vpConfig.h>
37
38#ifdef VISP_HAVE_CATCH2
39#define CATCH_CONFIG_RUNNER
40#include <catch.hpp>
41
42#include <visp3/core/vpGaussRand.h>
43#include <visp3/core/vpTime.h>
44#include <visp3/core/vpMath.h>
45
46namespace
47{
48class vpUniRandOld
49{
50 long a;
51 long m; // 2^31-1
52 long q; // integer part of m/a
53 long r; // r=m mod a
54 double normalizer; // we use a normalizer > m to ensure ans will never be 1
55 // (it is the case if x = 739806647)
56
57private:
58 inline void draw0()
59 {
60 long k = x / q; // temp value for computing without overflow
61 x = a * (x - k * q) - k * r;
62 if (x < 0)
63 x += m; // compute x without overflow
64 }
65
66protected:
67 long x;
68 double draw1()
69 {
70 const long ntab = 33; // we work on a 32 elements array.
71 // the 33rd one is actually the first value of y.
72 const long modulo = ntab - 2;
73
74 static long y = 0;
75 static long T[ntab];
76
77 long j; // index of T
78
79 // step 0
80 if (!y) { // first time
81 for (j = 0; j < ntab; j++) {
82 draw0();
83 T[j] = x;
84 } // compute table T
85 y = T[ntab - 1];
86 }
87
88 // step 1
89 j = y & modulo; // compute modulo ntab+1 (the first element is the 0th)
90
91 // step 3
92 y = T[j];
93 double ans = (double)y / normalizer;
94
95 // step 4
96 // generate x(k+i) and set y=x(k+i)
97 draw0();
98
99 // refresh T[j];
100 T[j] = x;
101
102 return ans;
103 }
104
105public:
107 explicit vpUniRandOld(const long seed = 0)
108 : a(16807), m(2147483647), q(127773), r(2836), normalizer(2147484721.0), x((seed) ? seed : 739806647)
109 {
110 }
111
113 virtual ~vpUniRandOld() {};
114
116 double operator()() { return draw1(); }
117};
118}
119
120TEST_CASE("Check Gaussian draw", "[visp_rand]") {
121 std::vector<double> vec(100000);
122 const double sigma = 5.0, mean = -7.5;
123 vpGaussRand rng(sigma, mean);
124
125 vpChrono chrono;
126 chrono.start();
127 for (size_t i = 0; i < vec.size(); i++) {
128 vec[i] = rng();
129 }
130 chrono.stop();
131
132 std::cout << vec.size() << " Gaussian draw performed in " << chrono.getDurationMs() << " ms" << std::endl;
133 double calculated_sigma = vpMath::getStdev(vec);
134 double calculated_mean = vpMath::getMean(vec);
135 std::cout << "Calculated sigma: " << calculated_sigma << std::endl;
136 std::cout << "Calculated mean: " << calculated_mean << std::endl;
137
138 CHECK(calculated_sigma == Approx(sigma).epsilon(0.01));
139 CHECK(calculated_mean == Approx(mean).epsilon(0.01));
140}
141
142TEST_CASE("Check Gaussian draw independance", "[visp_rand]") {
143 const double sigma = 5.0, mean = -7.5;
144
145 SECTION("Two simultaneous vpGaussRand instances with the same seed should produce the same results")
146 {
147 vpGaussRand rng1(sigma, mean), rng2(sigma, mean);
148
149 for (int i = 0; i < 10; i++) {
150 CHECK(rng1() == Approx(rng2()).margin(1e-6));
151 }
152 }
153 SECTION("Two vpGaussRand instances with the same seed should produce the same results")
154 {
155 std::vector<double> vec1, vec2;
156 {
157 vpGaussRand rng(sigma, mean);
158 for (int i = 0; i < 10; i++) {
159 vec1.push_back(rng());
160 }
161 }
162 {
163 vpGaussRand rng(sigma, mean);
164 for (int i = 0; i < 10; i++) {
165 vec2.push_back(rng());
166 }
167 }
168 REQUIRE(vec1.size() == vec2.size());
169
170 for (size_t i = 0; i < vec1.size(); i++) {
171 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
172 }
173 }
174}
175
176TEST_CASE("Check uniform draw", "[visp_rand]") {
177 const int niters = 500000;
178
179 SECTION("vpUniRand")
180 {
181 vpUniRand rng;
182 int inside = 0;
183
184 vpChrono chrono;
185 chrono.start();
186 for (int i = 0; i < niters; i++) {
187 double x = rng();
188 double y = rng();
189
190 if (sqrt(x*x + y * y) <= 1.0) {
191 inside++;
192 }
193 }
194 double pi = 4.0 * inside / niters;
195 chrono.stop();
196
197 double pi_error = pi - M_PI;
198 std::cout << "vpUniRand calculated pi: " << pi << " in " << chrono.getDurationMs() << " ms" << std::endl;
199 std::cout << "pi error: " << pi_error << std::endl;
200
201 CHECK(pi == Approx(M_PI).margin(0.005));
202 }
203
204 SECTION("C++ rand()")
205 {
206 srand(0);
207 int inside = 0;
208
209 vpChrono chrono;
210 chrono.start();
211 for (int i = 0; i < niters; i++) {
212 double x = static_cast<double>(rand()) / RAND_MAX;
213 double y = static_cast<double>(rand()) / RAND_MAX;
214
215 if (sqrt(x*x + y * y) <= 1.0) {
216 inside++;
217 }
218 }
219 double pi = 4.0 * inside / niters;
220 chrono.stop();
221
222 double pi_error = pi - M_PI;
223 std::cout << "C++ rand() calculated pi: " << pi << " in " << chrono.getDurationMs() << " ms" << std::endl;
224 std::cout << "pi error: " << pi_error << std::endl;
225
226 CHECK(pi == Approx(M_PI).margin(0.01));
227 }
228
229 SECTION("Old ViSP vpUniRand implementation")
230 {
231 vpUniRand rng;
232 int inside = 0;
233
234 vpChrono chrono;
235 chrono.start();
236 for (int i = 0; i < niters; i++) {
237 double x = rng();
238 double y = rng();
239
240 if (sqrt(x*x + y * y) <= 1.0) {
241 inside++;
242 }
243 }
244 double pi = 4.0 * inside / niters;
245 chrono.stop();
246
247 double pi_error = pi - M_PI;
248 std::cout << "Old ViSP vpUniRand implementation calculated pi: " << pi << " in "
249 << chrono.getDurationMs() << " ms" << std::endl;
250 std::cout << "pi error: " << pi_error << std::endl;
251
252 CHECK(pi == Approx(M_PI).margin(0.005));
253 }
254}
255
256TEST_CASE("Check uniform draw range", "[visp_rand]") {
257 const int niters = 1000;
258 vpUniRand rng;
259
260 SECTION("Check[0.0, 1.0) range")
261 {
262 for (int i = 0; i < niters; i++) {
263 double r = rng();
264 CHECK(r >= 0.0);
265 CHECK(r < 1.0);
266 }
267 }
268
269 SECTION("Check[-7, 10) range")
270 {
271 const int a = -7, b = 10;
272 for (int i = 0; i < niters; i++) {
273 int r = rng.uniform(a, b);
274 CHECK(r >= a);
275 CHECK(r < b);
276 }
277 }
278
279 SECTION("Check[-4.5f, 105.3f) range")
280 {
281 const float a = -4.5f, b = 105.3f;
282 for (int i = 0; i < niters; i++) {
283 float r = rng.uniform(a, b);
284 CHECK(r >= a);
285 CHECK(r < b);
286 }
287 }
288
289 SECTION("Check[14.6, 56.78) range")
290 {
291 const double a = 14.6, b = 56.78;
292 for (int i = 0; i < niters; i++) {
293 double r = rng.uniform(a, b);
294 CHECK(r >= a);
295 CHECK(r < b);
296 }
297 }
298}
299
300TEST_CASE("Check uniform draw independance", "[visp_rand]") {
301 SECTION("Two simultaneous vpUniRand instances with the same seed should produce the same results")
302 {
303 {
304 vpUniRand rng1, rng2;
305
306 for (int i = 0; i < 10; i++) {
307 CHECK(rng1.next() == rng2.next());
308 }
309 }
310 {
311 vpUniRand rng1, rng2;
312
313 for (int i = 0; i < 10; i++) {
314 CHECK(rng1.uniform(-1.0, 1.0) == Approx(rng2.uniform(-1.0, 1.0)).margin(1e-6));
315 }
316 }
317 }
318 SECTION("Two vpUniRand instances with the same seed should produce the same results")
319 {
320 {
321 std::vector<uint32_t> vec1, vec2;
322 {
323 vpUniRand rng;
324 for (int i = 0; i < 10; i++) {
325 vec1.push_back(rng.next());
326 }
327 }
328 {
329 vpUniRand rng;
330 for (int i = 0; i < 10; i++) {
331 vec2.push_back(rng.next());
332 }
333 }
334 REQUIRE(vec1.size() == vec2.size());
335
336 for (size_t i = 0; i < vec1.size(); i++) {
337 CHECK(vec1[i] == vec2[i]);
338 }
339 }
340 {
341 std::vector<double> vec1, vec2;
342 {
343 vpUniRand rng;
344 for (int i = 0; i < 10; i++) {
345 vec1.push_back(rng.uniform(-1.0, 2.0));
346 }
347 }
348 {
349 vpUniRand rng;
350 for (int i = 0; i < 10; i++) {
351 vec2.push_back(rng.uniform(-1.0, 2.0));
352 }
353 }
354 REQUIRE(vec1.size() == vec2.size());
355
356 for (size_t i = 0; i < vec1.size(); i++) {
357 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
358 }
359 }
360 }
361}
362
363int main(int argc, char* argv[])
364{
365 Catch::Session session; // There must be exactly one instance
366
367 // Let Catch (using Clara) parse the command line
368 session.applyCommandLine(argc, argv);
369
370 int numFailed = session.run();
371
372 // numFailed is clamped to 255 as some unices only use the lower 8 bits.
373 // This clamping has already been applied, so just return it here
374 // You can also do any post run clean-up here
375 return numFailed;
376}
377#else
378int main()
379{
380 return 0;
381}
382#endif
void start(bool reset=true)
Definition: vpTime.cpp:409
void stop()
Definition: vpTime.cpp:424
double getDurationMs()
Definition: vpTime.cpp:392
Class for generating random number with normal probability density.
Definition: vpGaussRand.h:121
static double getStdev(const std::vector< double > &v, bool useBesselCorrection=false)
Definition: vpMath.cpp:291
static double getMean(const std::vector< double > &v)
Definition: vpMath.cpp:241
Class for generating random numbers with uniform probability density.
Definition: vpUniRand.h:101
int uniform(int a, int b)
Definition: vpUniRand.cpp:163
uint32_t next()
Definition: vpUniRand.cpp:149