1 | #include <boost/random/variate_generator.hpp>
|
---|
2 | #include <boost/random/uniform_on_sphere.hpp>
|
---|
3 | #include <boost/random.hpp>
|
---|
4 | #include <iostream>
|
---|
5 | #include <vector>
|
---|
6 |
|
---|
7 | int main(int argc, const char *argv[])
|
---|
8 | {
|
---|
9 | // Typedefs to the generator, container, distribution and variate generator
|
---|
10 | typedef boost::mt19937 generator_t;
|
---|
11 | typedef std::vector<double> container_t;
|
---|
12 | typedef boost::uniform_on_sphere<double, container_t> spherical_dist_t;
|
---|
13 | typedef boost::variate_generator<generator_t*, spherical_dist_t > variate_generator_t;
|
---|
14 |
|
---|
15 | //Specify the dimensions to try and the the number of samples per dimension
|
---|
16 | const unsigned int minDim = 1u;
|
---|
17 | const unsigned int maxDim = 25u;
|
---|
18 | const unsigned int numSamples = 10000u;
|
---|
19 |
|
---|
20 | // Track the largest error found:
|
---|
21 | double maxError = 0.0;
|
---|
22 |
|
---|
23 | // Create the generator and specify the maximum dimension to try
|
---|
24 | generator_t gen;
|
---|
25 |
|
---|
26 | // Iterate over the dimensions drawing samples and announcing if they're not of magnitude 1.0
|
---|
27 | for (unsigned int n = minDim; n <= maxDim; ++n)
|
---|
28 | {
|
---|
29 | // Create the distribution and the variate generator
|
---|
30 | spherical_dist_t dist(n);
|
---|
31 | variate_generator_t variate(&gen, dist);
|
---|
32 |
|
---|
33 | // Draw some samples
|
---|
34 | for (unsigned int i = 0u; i < numSamples; ++i)
|
---|
35 | {
|
---|
36 | // Draw a sample
|
---|
37 | container_t sample = variate();
|
---|
38 |
|
---|
39 | // Calculate it's magnitude
|
---|
40 | double magnitude = 0.0;
|
---|
41 | for (unsigned int j = 0u; j < n; ++j)
|
---|
42 | {
|
---|
43 | magnitude = magnitude + sample.at(j)*sample.at(j);
|
---|
44 | }
|
---|
45 | magnitude = std::sqrt(magnitude);
|
---|
46 |
|
---|
47 | // Calculate the error in the magnitude:
|
---|
48 | double magnitudeError = std::abs(1.0 - magnitude);
|
---|
49 |
|
---|
50 | // Update the maximum error:
|
---|
51 | maxError = std::max(maxError, magnitudeError);
|
---|
52 |
|
---|
53 | // Check if this is above an arbitrary tolerance and is an example of the bug:
|
---|
54 | if ( magnitudeError > 1E-12)
|
---|
55 | {
|
---|
56 | std::cout << std::scientific << "A " << n << " dimensional vector has a magnitude of " << magnitude << " which is " << magnitudeError << " away from 1.0" << std::endl;
|
---|
57 | }
|
---|
58 | }
|
---|
59 | }
|
---|
60 |
|
---|
61 | // Statistics:
|
---|
62 | std::cout << "Drawing " << (maxDim - minDim + 1u)*numSamples << " samples from " << minDim << "D to " << maxDim << "D found a maximum absolute magnitude error of " << maxError << std::endl;
|
---|
63 |
|
---|
64 | return 0;
|
---|
65 | }
|
---|
66 |
|
---|