Visual Servoing Platform  version 3.4.1 under development (2021-10-17)
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 
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 }
119 
120 TEST_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 
142 TEST_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 
176 TEST_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 
256 TEST_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 
300 TEST_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 
363 int 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
378 int main()
379 {
380  return 0;
381 }
382 #endif
double getDurationMs()
Definition: vpTime.cpp:392
uint32_t next()
Definition: vpUniRand.cpp:149
static double getStdev(const std::vector< double > &v, bool useBesselCorrection=false)
Definition: vpMath.cpp:250
int uniform(int a, int b)
Definition: vpUniRand.cpp:163
void start(bool reset=true)
Definition: vpTime.cpp:409
static double getMean(const std::vector< double > &v)
Definition: vpMath.cpp:200
void stop()
Definition: vpTime.cpp:424
Class for generating random number with normal probability density.
Definition: vpGaussRand.h:120
Class for generating random numbers with uniform probability density.
Definition: vpUniRand.h:100