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 #include "ompl/util/RandomNumbers.h"
38 #include "ompl/util/Exception.h"
39 #include "ompl/util/Console.h"
40 #include <chrono>
41 #include <mutex>
42 #include <memory>
43 #include <boost/math/constants/constants.hpp>
44 #include <boost/scoped_ptr.hpp>
45 #include <boost/random/uniform_on_sphere.hpp>
46 #include <boost/random/variate_generator.hpp>
47 
49 namespace
50 {
54  class RNGSeedGenerator
55  {
56  public:
57  RNGSeedGenerator()
58  : firstSeed_(std::chrono::duration_cast<std::chrono::microseconds>(
59  std::chrono::high_resolution_clock::now().time_since_epoch())
60  .count())
61  , sGen_(firstSeed_)
62  , sDist_(1, 1000000000)
63  {
64  }
65 
66  std::uint_fast32_t firstSeed()
67  {
68  std::lock_guard<std::mutex> slock(rngMutex_);
69  return firstSeed_;
70  }
71 
72  void setSeed(std::uint_fast32_t seed)
73  {
74  std::lock_guard<std::mutex> slock(rngMutex_);
75  if (seed > 0)
76  {
77  if (someSeedsGenerated_)
78  {
79  OMPL_ERROR("Random number generation already started. Changing seed now will not lead to "
80  "deterministic sampling.");
81  }
82  else
83  {
84  // In this case, since no seeds have been generated yet, so we remember this seed as the first one.
85  firstSeed_ = seed;
86  }
87  }
88  else
89  {
90  if (someSeedsGenerated_)
91  {
92  OMPL_WARN("Random generator seed cannot be 0. Ignoring seed.");
93  return;
94  }
95  OMPL_WARN("Random generator seed cannot be 0. Using 1 instead.");
96  seed = 1;
97  }
98  sGen_.seed(seed);
99  }
100 
101  std::uint_fast32_t nextSeed()
102  {
103  std::lock_guard<std::mutex> slock(rngMutex_);
104  someSeedsGenerated_ = true;
105  return sDist_(sGen_);
106  }
107 
108  private:
109  bool someSeedsGenerated_{false};
110  std::uint_fast32_t firstSeed_;
111  std::mutex rngMutex_;
112  std::ranlux24_base sGen_;
113  std::uniform_int_distribution<> sDist_;
114  };
115 
116  std::once_flag g_once;
117  boost::scoped_ptr<RNGSeedGenerator> g_RNGSeedGenerator;
118 
119  void initRNGSeedGenerator()
120  {
121  g_RNGSeedGenerator.reset(new RNGSeedGenerator());
122  }
123 
124  RNGSeedGenerator &getRNGSeedGenerator()
125  {
126  std::call_once(g_once, &initRNGSeedGenerator);
127  return *g_RNGSeedGenerator;
128  }
129 } // namespace
131 
133 class ompl::RNG::SphericalData
134 {
135 public:
137  using container_type_t = std::vector<double>;
138 
140  using spherical_dist_t = boost::uniform_on_sphere<double, container_type_t>;
141 
143  using variate_generator_t = boost::variate_generator<std::mt19937 *, spherical_dist_t>;
144 
146  SphericalData(std::mt19937 *generatorPtr) : generatorPtr_(generatorPtr){};
147 
149  container_type_t generate(unsigned int dim)
150  {
151  // Assure that the dimension is in the range of the vector.
152  growVector(dim);
153 
154  // Assure that the dimension is allocated:
155  allocateDimension(dim);
156 
157  // Return the generator
158  return (*dimVector_.at(dim).second)();
159  };
160 
162  void reset()
163  {
164  // Iterate over each dimension
165  for (auto &i : dimVector_)
166  // Check if the variate_generator is allocated
167  if (bool(i.first))
168  // It is, reset THE DATA (not the pointer)
169  i.first->reset();
170  // No else, this is an uninitialized dimension.
171  };
172 
173 private:
175  using dist_gen_pair_t = std::pair<std::shared_ptr<spherical_dist_t>, std::shared_ptr<variate_generator_t>>;
176 
178  std::vector<dist_gen_pair_t> dimVector_;
179 
181  std::mt19937 *generatorPtr_;
182 
184  void growVector(unsigned int dim)
185  {
186  // Iterate until the index associated with this dimension is in the vector
187  while (dim >= dimVector_.size())
188  // Create a pair of empty pointers:
189  dimVector_.emplace_back();
190  };
191 
193  void allocateDimension(unsigned int dim)
194  {
195  // Only do this if unallocated, so check that:
196  if (dimVector_.at(dim).first == nullptr)
197  {
198  // It is not allocated, so....
199  // First construct the distribution
200  dimVector_.at(dim).first = std::make_shared<spherical_dist_t>(dim);
201  // Then the variate generator
202  dimVector_.at(dim).second = std::make_shared<variate_generator_t>(generatorPtr_, *dimVector_.at(dim).first);
203  }
204  // No else, the pointer is already allocated.
205  };
206 };
208 
209 std::uint_fast32_t ompl::RNG::getSeed()
210 {
211  return getRNGSeedGenerator().firstSeed();
212 }
213 
214 void ompl::RNG::setSeed(std::uint_fast32_t seed)
215 {
216  getRNGSeedGenerator().setSeed(seed);
217 }
218 
220  : localSeed_(getRNGSeedGenerator().nextSeed())
221  , generator_(localSeed_)
222  , sphericalDataPtr_(std::make_shared<SphericalData>(&generator_))
223 {
224 }
225 
226 ompl::RNG::RNG(std::uint_fast32_t localSeed)
227  : localSeed_(localSeed), generator_(localSeed_), sphericalDataPtr_(std::make_shared<SphericalData>(&generator_))
228 {
229 }
230 
231 void ompl::RNG::setLocalSeed(std::uint_fast32_t localSeed)
232 {
233  // Store the seed
234  localSeed_ = localSeed;
235 
236  // Change the generator's seed
237  generator_.seed(localSeed_);
238 
239  // Reset the distributions used by the variate generators, as they can cache values
240  uniDist_.reset();
241  normalDist_.reset();
242  sphericalDataPtr_->reset();
243 }
244 
245 double ompl::RNG::halfNormalReal(double r_min, double r_max, double focus)
246 {
247  assert(r_min <= r_max);
248 
249  const double mean = r_max - r_min;
250  double v = gaussian(mean, mean / focus);
251 
252  if (v > mean)
253  v = 2.0 * mean - v;
254  double r = v >= 0.0 ? v + r_min : r_min;
255  return r > r_max ? r_max : r;
256 }
257 
258 int ompl::RNG::halfNormalInt(int r_min, int r_max, double focus)
259 {
260  auto r = (int)floor(halfNormalReal((double)r_min, (double)(r_max) + 1.0, focus));
261  return (r > r_max) ? r_max : r;
262 }
263 
264 // From: "Uniform Random Rotations", Ken Shoemake, Graphics Gems III,
265 // pg. 124-132
266 void ompl::RNG::quaternion(double value[4])
267 {
268  double x0 = uniDist_(generator_);
269  double r1 = sqrt(1.0 - x0), r2 = sqrt(x0);
270  double t1 = 2.0 * boost::math::constants::pi<double>() * uniDist_(generator_),
271  t2 = 2.0 * boost::math::constants::pi<double>() * uniDist_(generator_);
272  double c1 = cos(t1), s1 = sin(t1);
273  double c2 = cos(t2), s2 = sin(t2);
274  value[0] = s1 * r1;
275  value[1] = c1 * r1;
276  value[2] = s2 * r2;
277  value[3] = c2 * r2;
278 }
279 
280 // From Effective Sampling and Distance Metrics for 3D Rigid Body Path Planning, by James Kuffner, ICRA 2004
281 void ompl::RNG::eulerRPY(double value[3])
282 {
283  value[0] = boost::math::constants::pi<double>() * (-2.0 * uniDist_(generator_) + 1.0);
284  value[1] = acos(1.0 - 2.0 * uniDist_(generator_)) - boost::math::constants::pi<double>() / 2.0;
285  value[2] = boost::math::constants::pi<double>() * (-2.0 * uniDist_(generator_) + 1.0);
286 }
287 
288 void ompl::RNG::uniformNormalVector(std::vector<double> &v)
289 {
290  // Generate a random value, the variate_generator is returning a shallow_array_adaptor, which will modify the value
291  // array:
292  v = sphericalDataPtr_->generate(v.size());
293 }
294 
295 // See: http://math.stackexchange.com/a/87238
296 void ompl::RNG::uniformInBall(double r, std::vector<double> &v)
297 {
298  // Draw a random point on the unit sphere
299  uniformNormalVector(v);
300 
301  // Draw a random radius scale
302  double radiusScale = r * std::pow(uniformReal(0.0, 1.0), 1.0 / static_cast<double>(v.size()));
303 
304  // Scale the point on the unit sphere
305  std::transform(v.begin(), v.end(), v.begin(), [radiusScale](double x) { return radiusScale * x; });
306 }
307 
308 void ompl::RNG::uniformProlateHyperspheroidSurface(const std::shared_ptr<const ProlateHyperspheroid> &phsPtr,
309  double value[])
310 {
311  // Variables
312  // The spherical point as a std::vector
313  std::vector<double> sphere(phsPtr->getDimension());
314 
315  // Get a random point on the sphere
316  uniformNormalVector(sphere);
317 
318  // Transform to the PHS
319  phsPtr->transform(&sphere[0], value);
320 }
321 
322 void ompl::RNG::uniformProlateHyperspheroid(const std::shared_ptr<const ProlateHyperspheroid> &phsPtr, double value[])
323 {
324  // Variables
325  // The spherical point as a std::vector
326  std::vector<double> sphere(phsPtr->getDimension());
327 
328  // Get a random point in the sphere
329  uniformInBall(1.0, sphere);
330 
331  // Transform to the PHS
332  phsPtr->transform(&sphere[0], value);
333 }
void uniformNormalVector(std::vector< double > &v)
Uniform random sampling of a unit-length vector. I.e., the surface of an n-ball. The return variable ...
void uniformProlateHyperspheroidSurface(const std::shared_ptr< const ProlateHyperspheroid > &phsPtr, double value[])
Uniform random sampling of the surface of a prolate hyperspheroid, a special symmetric type of n-dime...
point now()
Get the current time point.
Definition: Time.h:122
void transform(const double sphere[], double phs[]) const
Transform a point from a sphere to PHS. The return variable phs is expected to already exist.
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,...
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...
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 quaternion(double value[4])
Uniform random unit quaternion sampling. The computed value has the order (x,y,z,w)....
void uniformInBall(double r, std::vector< double > &v)
Uniform random sampling of the content of an n-ball, with a radius appropriately distributed between ...
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
unsigned int getDimension() const
The state dimension of the PHS.
void uniformProlateHyperspheroid(const std::shared_ptr< const ProlateHyperspheroid > &phsPtr, double value[])
Uniform random sampling of a prolate hyperspheroid, a special symmetric type of n-dimensional ellipse...
#define OMPL_ERROR(fmt,...)
Log a formatted error string.
Definition: Console.h:64
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...
void eulerRPY(double value[3])
Uniform random sampling of Euler roll-pitch-yaw angles, each in the range (-pi, pi]....
RNG()
Constructor. Always sets a different random seed.