Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
testRand.cpp
1/****************************************************************************
2 *
3 * ViSP, open source Visual Servoing Platform software.
4 * Copyright (C) 2005 - 2023 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 https://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/vpMath.h>
44#include <visp3/core/vpTime.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} // namespace
119
120TEST_CASE("Check Gaussian draw", "[visp_rand]")
121{
122 std::vector<double> vec(100000);
123 const double sigma = 5.0, mean = -7.5;
124 vpGaussRand rng(sigma, mean);
125
126 vpChrono chrono;
127 chrono.start();
128 for (size_t i = 0; i < vec.size(); i++) {
129 vec[i] = rng();
130 }
131 chrono.stop();
132
133 std::cout << vec.size() << " Gaussian draw performed in " << chrono.getDurationMs() << " ms" << std::endl;
134 double calculated_sigma = vpMath::getStdev(vec);
135 double calculated_mean = vpMath::getMean(vec);
136 std::cout << "Calculated sigma: " << calculated_sigma << std::endl;
137 std::cout << "Calculated mean: " << calculated_mean << std::endl;
138
139 CHECK(calculated_sigma == Approx(sigma).epsilon(0.01));
140 CHECK(calculated_mean == Approx(mean).epsilon(0.01));
141}
142
143TEST_CASE("Check Gaussian draw independance", "[visp_rand]")
144{
145 const double sigma = 5.0, mean = -7.5;
146
147 SECTION("Two simultaneous vpGaussRand instances with the same seed should produce the same results")
148 {
149 vpGaussRand rng1(sigma, mean), rng2(sigma, mean);
150
151 for (int i = 0; i < 10; i++) {
152 CHECK(rng1() == Approx(rng2()).margin(1e-6));
153 }
154 }
155 SECTION("Two vpGaussRand instances with the same seed should produce the same results")
156 {
157 std::vector<double> vec1, vec2;
158 {
159 vpGaussRand rng(sigma, mean);
160 for (int i = 0; i < 10; i++) {
161 vec1.push_back(rng());
162 }
163 }
164 {
165 vpGaussRand rng(sigma, mean);
166 for (int i = 0; i < 10; i++) {
167 vec2.push_back(rng());
168 }
169 }
170 REQUIRE(vec1.size() == vec2.size());
171
172 for (size_t i = 0; i < vec1.size(); i++) {
173 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
174 }
175 }
176}
177
178TEST_CASE("Check uniform draw", "[visp_rand]")
179{
180 const int niters = 500000;
181
182 SECTION("vpUniRand")
183 {
184 vpUniRand rng;
185 int inside = 0;
186
187 vpChrono chrono;
188 chrono.start();
189 for (int i = 0; i < niters; i++) {
190 double x = rng();
191 double y = rng();
192
193 if (sqrt(x * x + y * y) <= 1.0) {
194 inside++;
195 }
196 }
197 double pi = 4.0 * inside / niters;
198 chrono.stop();
199
200 double pi_error = pi - M_PI;
201 std::cout << "vpUniRand calculated pi: " << pi << " in " << chrono.getDurationMs() << " ms" << std::endl;
202 std::cout << "pi error: " << pi_error << std::endl;
203
204 CHECK(pi == Approx(M_PI).margin(0.005));
205 }
206
207 SECTION("C++ rand()")
208 {
209 srand(0);
210 int inside = 0;
211
212 vpChrono chrono;
213 chrono.start();
214 for (int i = 0; i < niters; i++) {
215 double x = static_cast<double>(rand()) / RAND_MAX;
216 double y = static_cast<double>(rand()) / RAND_MAX;
217
218 if (sqrt(x * x + y * y) <= 1.0) {
219 inside++;
220 }
221 }
222 double pi = 4.0 * inside / niters;
223 chrono.stop();
224
225 double pi_error = pi - M_PI;
226 std::cout << "C++ rand() calculated pi: " << pi << " in " << chrono.getDurationMs() << " ms" << std::endl;
227 std::cout << "pi error: " << pi_error << std::endl;
228
229 CHECK(pi == Approx(M_PI).margin(0.01));
230 }
231
232 SECTION("Old ViSP vpUniRand implementation")
233 {
234 vpUniRand rng;
235 int inside = 0;
236
237 vpChrono chrono;
238 chrono.start();
239 for (int i = 0; i < niters; i++) {
240 double x = rng();
241 double y = rng();
242
243 if (sqrt(x * x + y * y) <= 1.0) {
244 inside++;
245 }
246 }
247 double pi = 4.0 * inside / niters;
248 chrono.stop();
249
250 double pi_error = pi - M_PI;
251 std::cout << "Old ViSP vpUniRand implementation calculated pi: " << pi << " in " << chrono.getDurationMs() << " ms"
252 << std::endl;
253 std::cout << "pi error: " << pi_error << std::endl;
254
255 CHECK(pi == Approx(M_PI).margin(0.005));
256 }
257}
258
259TEST_CASE("Check uniform draw range", "[visp_rand]")
260{
261 const int niters = 1000;
262 vpUniRand rng;
263
264 SECTION("Check[0.0, 1.0) range")
265 {
266 for (int i = 0; i < niters; i++) {
267 double r = rng();
268 CHECK(r >= 0.0);
269 CHECK(r < 1.0);
270 }
271 }
272
273 SECTION("Check[-7, 10) range")
274 {
275 const int a = -7, b = 10;
276 for (int i = 0; i < niters; i++) {
277 int r = rng.uniform(a, b);
278 CHECK(r >= a);
279 CHECK(r < b);
280 }
281 }
282
283 SECTION("Check[-4.5f, 105.3f) range")
284 {
285 const float a = -4.5f, b = 105.3f;
286 for (int i = 0; i < niters; i++) {
287 float r = rng.uniform(a, b);
288 CHECK(r >= a);
289 CHECK(r < b);
290 }
291 }
292
293 SECTION("Check[14.6, 56.78) range")
294 {
295 const double a = 14.6, b = 56.78;
296 for (int i = 0; i < niters; i++) {
297 double r = rng.uniform(a, b);
298 CHECK(r >= a);
299 CHECK(r < b);
300 }
301 }
302}
303
304TEST_CASE("Check uniform draw independance", "[visp_rand]")
305{
306 SECTION("Two simultaneous vpUniRand instances with the same seed should produce the same results")
307 {
308 {
309 vpUniRand rng1, rng2;
310
311 for (int i = 0; i < 10; i++) {
312 CHECK(rng1.next() == rng2.next());
313 }
314 }
315 {
316 vpUniRand rng1, rng2;
317
318 for (int i = 0; i < 10; i++) {
319 CHECK(rng1.uniform(-1.0, 1.0) == Approx(rng2.uniform(-1.0, 1.0)).margin(1e-6));
320 }
321 }
322 }
323 SECTION("Two vpUniRand instances with the same seed should produce the same results")
324 {
325 {
326 std::vector<uint32_t> vec1, vec2;
327 {
328 vpUniRand rng;
329 for (int i = 0; i < 10; i++) {
330 vec1.push_back(rng.next());
331 }
332 }
333 {
334 vpUniRand rng;
335 for (int i = 0; i < 10; i++) {
336 vec2.push_back(rng.next());
337 }
338 }
339 REQUIRE(vec1.size() == vec2.size());
340
341 for (size_t i = 0; i < vec1.size(); i++) {
342 CHECK(vec1[i] == vec2[i]);
343 }
344 }
345 {
346 std::vector<double> vec1, vec2;
347 {
348 vpUniRand rng;
349 for (int i = 0; i < 10; i++) {
350 vec1.push_back(rng.uniform(-1.0, 2.0));
351 }
352 }
353 {
354 vpUniRand rng;
355 for (int i = 0; i < 10; i++) {
356 vec2.push_back(rng.uniform(-1.0, 2.0));
357 }
358 }
359 REQUIRE(vec1.size() == vec2.size());
360
361 for (size_t i = 0; i < vec1.size(); i++) {
362 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
363 }
364 }
365 }
366}
367
368int main(int argc, char *argv[])
369{
370 Catch::Session session; // There must be exactly one instance
371
372 // Let Catch (using Clara) parse the command line
373 session.applyCommandLine(argc, argv);
374
375 int numFailed = session.run();
376
377 // numFailed is clamped to 255 as some unices only use the lower 8 bits.
378 // This clamping has already been applied, so just return it here
379 // You can also do any post run clean-up here
380 return numFailed;
381}
382#else
383int main() { return EXIT_SUCCESS; }
384#endif
void start(bool reset=true)
Definition vpTime.cpp:397
void stop()
Definition vpTime.cpp:412
double getDurationMs()
Definition vpTime.cpp:386
Class for generating random number with normal probability density.
static double getStdev(const std::vector< double > &v, bool useBesselCorrection=false)
Definition vpMath.cpp:345
static double getMean(const std::vector< double > &v)
Definition vpMath.cpp:294
Class for generating random numbers with uniform probability density.
Definition vpUniRand.h:122
int uniform(int a, int b)
uint32_t next()