boost/random/uniform_on_sphere.hpp

0001 /* boost random/uniform_on_sphere.hpp header file
0002  *
0003  * Copyright Jens Maurer 2000-2001
0004  * Copyright Steven Watanabe 2011
0005  * Distributed under the Boost Software License, Version 1.0. (See
0006  * accompanying file LICENSE_1_0.txt or copy at
0007  * http://www.boost.org/LICENSE_1_0.txt)
0008  *
0009  * See http://www.boost.org for most recent version including documentation.
0010  *
0011  * $Id$
0012  *
0013  * Revision history
0014  *  2001-02-18  moved to individual header files
0015  */
0016
0017 #ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
0018 #define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
0019
0020 #include <vector>
0021 #include <algorithm>     // std::transform
0022 #include <functional>    // std::bind2nd, std::divides
0023 #include <boost/assert.hpp>
0024 #include <boost/random/detail/config.hpp>
0025 #include <boost/random/detail/operators.hpp>
0026 #include <boost/random/normal_distribution.hpp>
0027
0028 namespace boost {
0029 namespace random {
0030
0031 /**
0032  * Instantiations of class template uniform_on_sphere model a
0033  * \random_distribution. Such a distribution produces random
0034  * numbers uniformly distributed on the unit sphere of arbitrary
0035  * dimension @c dim. The @c Cont template parameter must be a STL-like
0036  * container type with begin and end operations returning non-const
0037  * ForwardIterators of type @c Cont::iterator.
0038  */
0039 template<class RealType = double, class Cont = std::vector<RealType> >
0040 class uniform_on_sphere
0041 {
0042 public:
0043     typedef RealType input_type;
0044     typedef Cont result_type;
0045
0046     class param_type
0047     {
0048     public:
0049
0050         typedef uniform_on_sphere distribution_type;
0051
0052         /**
0053          * Constructs the parameters of a uniform_on_sphere
0054          * distribution, given the dimension of the sphere.
0055          */
0056         explicit param_type(int dim_arg = 2) : _dim(dim_arg)
0057         {
0058             BOOST_ASSERT(_dim >= 0);
0059         }
0060
0061         /** Returns the dimension of the sphere. */
0062         int dim() const { return _dim; }
0063
0064         /** Writes the parameters to a @c std::ostream. */
0065         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0066         {
0067             os << parm._dim;
0068             return os;
0069         }
0070
0071         /** Reads the parameters from a @c std::istream. */
0072         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0073         {
0074             is >> parm._dim;
0075             return is;
0076         }
0077
0078         /** Returns true if the two sets of parameters are equal. */
0079         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0080         { return lhs._dim == rhs._dim; }
0081
0082         /** Returns true if the two sets of parameters are different. */
0083         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0084
0085     private:
0086         int _dim;
0087     };
0088
0089     /**
0090      * Constructs a @c uniform_on_sphere distribution.
0091      * @c dim is the dimension of the sphere.
0092      *
0093      * Requires: dim >= 0
0094      */
0095     explicit uniform_on_sphere(int dim_arg = 2)
0096       : _container(dim_arg), _dim(dim_arg) { }
0097
0098     /**
0099      * Constructs a @c uniform_on_sphere distribution from its parameters.
0100      */
0101     explicit uniform_on_sphere(const param_type& parm)
0102       : _container(parm.dim()), _dim(parm.dim()) { }
0103
0104     // compiler-generated copy ctor and assignment operator are fine
0105
0106     /** Returns the dimension of the sphere. */
0107     int dim() const { return _dim; }
0108
0109     /** Returns the parameters of the distribution. */
0110     param_type param() const { return param_type(_dim); }
0111     /** Sets the parameters of the distribution. */
0112     void param(const param_type& parm)
0113     {
0114         _dim = parm.dim();
0115         _container.resize(_dim);
0116     }
0117
0118     /**
0119      * Returns the smallest value that the distribution can produce.
0120      * Note that this is required to approximate the standard library's
0121      * requirements.  The behavior is defined according to lexicographical
0122      * comparison so that for a container type of std::vector,
0123      * dist.min() <= x <= dist.max() where x is any value produced
0124      * by the distribution.
0125      */
0126     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0127     {
0128         result_type result(_dim);
0129         if(_dim != 0) {
0130             result.front() = RealType(-1.0);
0131         }
0132         return result;
0133     }
0134     /**
0135      * Returns the largest value that the distribution can produce.
0136      * Note that this is required to approximate the standard library's
0137      * requirements.  The behavior is defined according to lexicographical
0138      * comparison so that for a container type of std::vector,
0139      * dist.min() <= x <= dist.max() where x is any value produced
0140      * by the distribution.
0141      */
0142     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0143     {
0144         result_type result(_dim);
0145         if(_dim != 0) {
0146             result.front() = RealType(1.0);
0147         }
0148         return result;
0149     }
0150
0151     /**
0152      * Effects: Subsequent uses of the distribution do not depend
0153      * on values produced by any engine prior to invoking reset.
0154      */
0155     void reset() {}
0156
0157     /**
0158      * Returns a point uniformly distributed over the surface of
0159      * a sphere of dimension dim().
0160      */
0161     template<class Engine>
0162     const result_type & operator()(Engine& eng)
0163     {
0164         using std::sqrt;
0165         switch(_dim)
0166         {
0167         case 0: break;
0168         case 1:
0169             {
0170                 if(uniform_01<RealType>()(eng) < 0.5) {
0171                     *_container.begin() = -1;
0172                 } else {
0173                     *_container.begin() = 1;
0174                 }
0175                 break;
0176             }
0177         case 2:
0178             {
0179                 uniform_01<RealType> uniform;
0180                 RealType sqsum;
0181                 RealType x, y;
0182                 do {
0183                     x = uniform(eng) * 2 - 1;
0184                     y = uniform(eng) * 2 - 1;
0185                     sqsum = x*x + y*y;
0186                 } while(sqsum == 0 || sqsum > 1);
0187                 RealType mult = 1/sqrt(sqsum);
0188                 typename Cont::iterator iter = _container.begin();
0189                 *iter = x * mult;
0190                 iter++;
0191                 *iter = y * mult;
0192                 break;
0193             }
0194         case 3:
0195             {
0196                 uniform_01<RealType> uniform;
0197                 RealType sqsum;
0198                 RealType x, y;
0199                 do {
0200                     x = uniform(eng) * 2 - 1;
0201                     y = uniform(eng) * 2 - 1;
0202                     sqsum = x*x + y*y;
0203                 } while(sqsum > 1);
0204                 RealType mult = 2 * sqrt(1 - sqsum);
0205                 typename Cont::iterator iter = _container.begin();
0206                 *iter = x * mult;
0207                 ++iter;
0208                 *iter = y * mult;
0209                 ++iter;
0210                 *iter = 2 * sqsum - 1;
0211                 break;
0212             }
0213         default:
0214             {
0215                 detail::unit_normal_distribution<RealType> normal;
0216                 RealType sqsum;
0217                 do {
0218                     sqsum = 0;
0219                     for(typename Cont::iterator it = _container.begin();
0220                         it != _container.end();
0221                         ++it) {
0222                         RealType val = normal(eng);
0223                         *it = val;
0224                         sqsum += val * val;
0225                     }
0226                 } while(sqsum == 0);
0227                 // for all i: result[i] /= sqrt(sqsum)
0228                 RealType inverse_distance = 1 / sqrt(sqsum);
0229                 for(typename Cont::iterator it = _container.begin();
0230                     it != _container.end();
0231                     ++it) {
0232                     *it *= inverse_distance;
0233                 }
0234             }
0235         }
0236         return _container;
0237     }
0238
0239     /**
0240      * Returns a point uniformly distributed over the surface of
0241      * a sphere of dimension param.dim().
0242      */
0243     template<class Engine>
0244     result_type operator()(Engine& eng, const param_type& parm) const
0245     {
0246         return uniform_on_sphere(parm)(eng);
0247     }
0248
0249     /** Writes the distribution to a @c std::ostream. */
0250     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd)
0251     {
0252         os << sd._dim;
0253         return os;
0254     }
0255
0256     /** Reads the distribution from a @c std::istream. */
0257     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd)
0258     {
0259         is >> sd._dim;
0260         sd._container.resize(sd._dim);
0261         return is;
0262     }
0263
0264     /**
0265      * Returns true if the two distributions will produce identical
0266      * sequences of values, given equal generators.
0267      */
0268     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
0269     { return lhs._dim == rhs._dim; }
0270
0271     /**
0272      * Returns true if the two distributions may produce different
0273      * sequences of values, given equal generators.
0274      */
0275     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)
0276
0277 private:
0278     result_type _container;
0279     int _dim;
0280 };
0281
0282 } // namespace random
0283
0284 using random::uniform_on_sphere;
0285
0286 } // namespace boost
0287
0288 #endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP