Visual Servoing Platform  version 3.6.1 under development (2024-02-13)
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 
46 namespace
47 {
48 class 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 
57 private:
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 
66 protected:
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 
105 public:
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 
120 TEST_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 
143 TEST_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 
178 TEST_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 
259 TEST_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 
304 TEST_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 
368 int 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
383 int 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.
Definition: vpGaussRand.h:116
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:123
int uniform(int a, int b)
Definition: vpUniRand.cpp:160
uint32_t next()
Definition: vpUniRand.cpp:146