RandomNumbers.cpp
1 /*********************************************************************
2 * Software License Agreement (BSD License)
3 *
4 * Copyright (c) 2008, Willow Garage, Inc.
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 *
11 * * Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 * * Redistributions in binary form must reproduce the above
14 * copyright notice, this list of conditions and the following
15 * disclaimer in the documentation and/or other materials provided
16 * with the distribution.
17 * * Neither the name of the Willow Garage nor the names of its
18 * contributors may be used to endorse or promote products derived
19 * from this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
28 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
29 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
31 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
32 * POSSIBILITY OF SUCH DAMAGE.
33 *********************************************************************/
34 
35 /* Author: Ioan Sucan, Jonathan Gammell*/
36 
37 // Enable the use of shallow_array_adaptor to create a uBLAS-vector-view of C-style array without copying data
38 #define BOOST_UBLAS_SHALLOW_ARRAY_ADAPTOR
39 
40 #include "ompl/util/RandomNumbers.h"
41 #include "ompl/util/Exception.h"
42 #include "ompl/util/Console.h"
43 #include <mutex>
44 #include <memory>
45 #include <boost/math/constants/constants.hpp>
46 #include <boost/scoped_ptr.hpp>
47 #include <boost/random/uniform_on_sphere.hpp>
48 #include <boost/random/variate_generator.hpp>
49 // For boost::numeric::ublas::shallow_array_adaptor:
50 #include <boost/numeric/ublas/vector.hpp>
51 
53 namespace
54 {
58  class RNGSeedGenerator
59  {
60  public:
61  RNGSeedGenerator()
62  : firstSeed_(std::chrono::duration_cast<std::chrono::microseconds>(
63  std::chrono::system_clock::now() - std::chrono::system_clock::time_point::min()).count())
64  , sGen_(firstSeed_)
65  , sDist_(1, 1000000000)
66  {
67  }
68 
69  std::uint_fast32_t firstSeed()
70  {
71  std::lock_guard<std::mutex> slock(rngMutex_);
72  return firstSeed_;
73  }
74 
75  void setSeed(std::uint_fast32_t seed)
76  {
77  std::lock_guard<std::mutex> slock(rngMutex_);
78  if (seed > 0)
79  {
80  if (someSeedsGenerated_)
81  {
82  OMPL_ERROR("Random number generation already started. Changing seed now will not lead to "
83  "deterministic sampling.");
84  }
85  else
86  {
87  // In this case, since no seeds have been generated yet, so we remember this seed as the first one.
88  firstSeed_ = seed;
89  }
90  }
91  else
92  {
93  if (someSeedsGenerated_)
94  {
95  OMPL_WARN("Random generator seed cannot be 0. Ignoring seed.");
96  return;
97  }
98  OMPL_WARN("Random generator seed cannot be 0. Using 1 instead.");
99  seed = 1;
100  }
101  sGen_.seed(seed);
102  }
103 
104  std::uint_fast32_t nextSeed()
105  {
106  std::lock_guard<std::mutex> slock(rngMutex_);
107  someSeedsGenerated_ = true;
108  return sDist_(sGen_);
109  }
110 
111  private:
112  bool someSeedsGenerated_{false};
113  std::uint_fast32_t firstSeed_;
114  std::mutex rngMutex_;
115  std::ranlux24_base sGen_;
116  std::uniform_int_distribution<> sDist_;
117  };
118 
119  std::once_flag g_once;
120  boost::scoped_ptr<RNGSeedGenerator> g_RNGSeedGenerator;
121 
122  void initRNGSeedGenerator()
123  {
124  g_RNGSeedGenerator.reset(new RNGSeedGenerator());
125  }
126 
127  RNGSeedGenerator &getRNGSeedGenerator()
128  {
129  std::call_once(g_once, &initRNGSeedGenerator);
130  return *g_RNGSeedGenerator;
131  }
132 } // namespace
134 
136 class ompl::RNG::SphericalData
137 {
138 public:
140  using container_type_t = boost::numeric::ublas::shallow_array_adaptor<double>;
141 
143  using spherical_dist_t = boost::uniform_on_sphere<double, container_type_t>;
144 
146  using variate_generator_t = boost::variate_generator<std::mt19937 *, spherical_dist_t>;
147 
149  SphericalData(std::mt19937 *generatorPtr) : generatorPtr_(generatorPtr){};
150 
152  container_type_t generate(unsigned int dim)
153  {
154  // Assure that the dimension is in the range of the vector.
155  growVector(dim);
156 
157  // Assure that the dimension is allocated:
158  allocateDimension(dim);
159 
160  // Return the generator
161  return (*dimVector_.at(dim).second)();
162  };
163 
165  void reset()
166  {
167  // Iterate over each dimension
168  for (auto &i : dimVector_)
169  {
170  // Check if the variate_generator is allocated
171  if (bool(i.first))
172  {
173  // It is, reset THE DATA (not the pointer)
174  i.first->reset();
175  }
176  // No else, this is an uninitialized dimension.
177  }
178  };
179 
180 private:
182  using dist_gen_pair_t = std::pair<std::shared_ptr<spherical_dist_t>, std::shared_ptr<variate_generator_t>>;
183 
185  std::vector<dist_gen_pair_t> dimVector_;
186 
188  std::mt19937 *generatorPtr_;
189 
191  void growVector(unsigned int dim)
192  {
193  // Iterate until the index associated with this dimension is in the vector
194  while (dim >= dimVector_.size())
195  {
196  // Create a pair of empty pointers:
197  dimVector_.emplace_back();
198  }
199  };
200 
202  void allocateDimension(unsigned int dim)
203  {
204  // Only do this if unallocated, so check that:
205  if (dimVector_.at(dim).first == nullptr)
206  {
207  // It is not allocated, so....
208  // First construct the distribution
209  dimVector_.at(dim).first = std::make_shared<spherical_dist_t>(dim);
210  // Then the variate generator
211  dimVector_.at(dim).second = std::make_shared<variate_generator_t>(generatorPtr_, *dimVector_.at(dim).first);
212  }
213  // No else, the pointer is already allocated.
214  };
215 };
217 
218 std::uint_fast32_t ompl::RNG::getSeed()
219 {
220  return getRNGSeedGenerator().firstSeed();
221 }
222 
223 void ompl::RNG::setSeed(std::uint_fast32_t seed)
224 {
225  getRNGSeedGenerator().setSeed(seed);
226 }
227 
229  : localSeed_(getRNGSeedGenerator().nextSeed())
230  , generator_(localSeed_)
231  , sphericalDataPtr_(std::make_shared<SphericalData>(&generator_))
232 {
233 }
234 
235 ompl::RNG::RNG(std::uint_fast32_t localSeed)
236  : localSeed_(localSeed)
237  , generator_(localSeed_)
238  , sphericalDataPtr_(std::make_shared<SphericalData>(&generator_))
239 {
240 }
241 
242 void ompl::RNG::setLocalSeed(std::uint_fast32_t localSeed)
243 {
244  // Store the seed
245  localSeed_ = localSeed;
246 
247  // Change the generator's seed
248  generator_.seed(localSeed_);
249 
250  // Reset the distributions used by the variate generators, as they can cache values
251  uniDist_.reset();
252  normalDist_.reset();
253  sphericalDataPtr_->reset();
254 }
255 
256 double ompl::RNG::halfNormalReal(double r_min, double r_max, double focus)
257 {
258  assert(r_min <= r_max);
259 
260  const double mean = r_max - r_min;
261  double v = gaussian(mean, mean / focus);
262 
263  if (v > mean)
264  v = 2.0 * mean - v;
265  double r = v >= 0.0 ? v + r_min : r_min;
266  return r > r_max ? r_max : r;
267 }
268 
269 int ompl::RNG::halfNormalInt(int r_min, int r_max, double focus)
270 {
271  auto r = (int)floor(halfNormalReal((double)r_min, (double)(r_max) + 1.0, focus));
272  return (r > r_max) ? r_max : r;
273 }
274 
275 // From: "Uniform Random Rotations", Ken Shoemake, Graphics Gems III,
276 // pg. 124-132
277 void ompl::RNG::quaternion(double value[4])
278 {
279  double x0 = uniDist_(generator_);
280  double r1 = sqrt(1.0 - x0), r2 = sqrt(x0);
281  double t1 = 2.0 * boost::math::constants::pi<double>() * uniDist_(generator_),
282  t2 = 2.0 * boost::math::constants::pi<double>() * uniDist_(generator_);
283  double c1 = cos(t1), s1 = sin(t1);
284  double c2 = cos(t2), s2 = sin(t2);
285  value[0] = s1 * r1;
286  value[1] = c1 * r1;
287  value[2] = s2 * r2;
288  value[3] = c2 * r2;
289 }
290 
291 // From Effective Sampling and Distance Metrics for 3D Rigid Body Path Planning, by James Kuffner, ICRA 2004
292 void ompl::RNG::eulerRPY(double value[3])
293 {
294  value[0] = boost::math::constants::pi<double>() * (-2.0 * uniDist_(generator_) + 1.0);
295  value[1] = acos(1.0 - 2.0 * uniDist_(generator_)) - boost::math::constants::pi<double>() / 2.0;
296  value[2] = boost::math::constants::pi<double>() * (-2.0 * uniDist_(generator_) + 1.0);
297 }
298 
299 void ompl::RNG::uniformNormalVector(unsigned int n, double value[])
300 {
301  // Create a uBLAS-vector-view of the C-style array without copying data
302  SphericalData::container_type_t rVector(n, value);
303 
304  // Generate a random value, the variate_generator is returning a shallow_array_adaptor, which will modify the value
305  // array:
306  rVector = sphericalDataPtr_->generate(n);
307 }
308 
309 // See: http://math.stackexchange.com/a/87238
310 void ompl::RNG::uniformInBall(double r, unsigned int n, double value[])
311 {
312  // Draw a random point on the unit sphere
313  uniformNormalVector(n, value);
314 
315  // Draw a random radius scale
316  double radiusScale = r * std::pow(uniformReal(0.0, 1.0), 1.0 / static_cast<double>(n));
317 
318  // Scale the point on the unit sphere
319  for (unsigned int i = 0u; i < n; ++i)
320  {
321  value[i] *= radiusScale;
322  }
323 }
324 
325 #if OMPL_HAVE_EIGEN3
326 void ompl::RNG::uniformProlateHyperspheroidSurface(const std::shared_ptr<const ProlateHyperspheroid> &phsPtr,
327  double value[])
328 {
329  // Variables
330  // The spherical point as a std::vector
331  std::vector<double> sphere(phsPtr->getDimension());
332 
333  // Get a random point on the sphere
334  uniformNormalVector(phsPtr->getDimension(), &sphere[0]);
335 
336  // Transform to the PHS
337  phsPtr->transform(&sphere[0], value);
338 }
339 
340 void ompl::RNG::uniformProlateHyperspheroid(const std::shared_ptr<const ProlateHyperspheroid> &phsPtr, double value[])
341 {
342  // Variables
343  // The spherical point as a std::vector
344  std::vector<double> sphere(phsPtr->getDimension());
345 
346  // Get a random point in the sphere
347  uniformInBall(1.0, phsPtr->getDimension(), &sphere[0]);
348 
349  // Transform to the PHS
350  phsPtr->transform(&sphere[0], value);
351 }
352 #endif
STL namespace.
void quaternion(double value[4])
Uniform random unit quaternion sampling. The computed value has the order (x,y,z,w). The return variable value is expected to already exist.
RNG()
Constructor. Always sets a different random seed.
void eulerRPY(double value[3])
Uniform random sampling of Euler roll-pitch-yaw angles, each in the range (-pi, pi]. The computed value has the order (roll, pitch, yaw). The return variable value is expected to already exist.
int halfNormalInt(int r_min, int r_max, double focus=3.0)
Generate a random integer using a half-normal distribution. The value is within specified bounds ([r_...
void uniformNormalVector(unsigned int n, double value[])
Uniform random sampling of a unit-length vector. I.e., the surface of an n-ball. The return variable ...
#define OMPL_ERROR(fmt,...)
Log a formatted error string.
Definition: Console.h:64
double uniformReal(double lower_bound, double upper_bound)
Generate a random real within given bounds: [lower_bound, upper_bound)
Definition: RandomNumbers.h:74
static std::uint_fast32_t getSeed()
Get the seed used to generate the seeds of each RNG instance. Passing the returned value to setSeed()...
#define OMPL_WARN(fmt,...)
Log a formatted warning string.
Definition: Console.h:66
double halfNormalReal(double r_min, double r_max, double focus=3.0)
Generate a random real using a half-normal distribution. The value is within specified bounds [r_min...
point now()
Get the current time point.
Definition: Time.h:70
void setLocalSeed(std::uint_fast32_t localSeed)
Set the seed used for the instance of a RNG. Use this function to ensure that an instance of an RNG g...
static void setSeed(std::uint_fast32_t seed)
Set the seed used to generate the seeds of each RNG instance. Use this function to ensure the same se...
void uniformInBall(double r, unsigned int n, double value[])
Uniform random sampling of the content of an n-ball, with a radius appropriately distributed between ...
double gaussian(double mean, double stddev)
Generate a random real using a normal distribution with given mean and variance.