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