Chris@16: /* boost random/uniform_on_sphere.hpp header file
Chris@16:  *
Chris@16:  * Copyright Jens Maurer 2000-2001
Chris@16:  * Copyright Steven Watanabe 2011
Chris@16:  * Distributed under the Boost Software License, Version 1.0. (See
Chris@16:  * accompanying file LICENSE_1_0.txt or copy at
Chris@16:  * http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  *
Chris@16:  * See http://www.boost.org for most recent version including documentation.
Chris@16:  *
Chris@101:  * $Id$
Chris@16:  *
Chris@16:  * Revision history
Chris@16:  *  2001-02-18  moved to individual header files
Chris@16:  */
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
Chris@16: #define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
Chris@16: 
Chris@16: #include <vector>
Chris@16: #include <algorithm>     // std::transform
Chris@16: #include <functional>    // std::bind2nd, std::divides
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/random/detail/config.hpp>
Chris@16: #include <boost/random/detail/operators.hpp>
Chris@16: #include <boost/random/normal_distribution.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace random {
Chris@16: 
Chris@16: /**
Chris@16:  * Instantiations of class template uniform_on_sphere model a
Chris@16:  * \random_distribution. Such a distribution produces random
Chris@16:  * numbers uniformly distributed on the unit sphere of arbitrary
Chris@16:  * dimension @c dim. The @c Cont template parameter must be a STL-like
Chris@16:  * container type with begin and end operations returning non-const
Chris@16:  * ForwardIterators of type @c Cont::iterator. 
Chris@16:  */
Chris@16: template<class RealType = double, class Cont = std::vector<RealType> >
Chris@16: class uniform_on_sphere
Chris@16: {
Chris@16: public:
Chris@16:     typedef RealType input_type;
Chris@16:     typedef Cont result_type;
Chris@16: 
Chris@16:     class param_type
Chris@16:     {
Chris@16:     public:
Chris@16: 
Chris@16:         typedef uniform_on_sphere distribution_type;
Chris@16: 
Chris@16:         /**
Chris@16:          * Constructs the parameters of a uniform_on_sphere
Chris@16:          * distribution, given the dimension of the sphere.
Chris@16:          */
Chris@16:         explicit param_type(int dim_arg = 2) : _dim(dim_arg)
Chris@16:         {
Chris@16:             BOOST_ASSERT(_dim >= 0);
Chris@16:         }
Chris@16: 
Chris@16:         /** Returns the dimension of the sphere. */
Chris@16:         int dim() const { return _dim; }
Chris@16: 
Chris@16:         /** Writes the parameters to a @c std::ostream. */
Chris@16:         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
Chris@16:         {
Chris@16:             os << parm._dim;
Chris@16:             return os;
Chris@16:         }
Chris@16: 
Chris@16:         /** Reads the parameters from a @c std::istream. */
Chris@16:         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
Chris@16:         {
Chris@16:             is >> parm._dim;
Chris@16:             return is;
Chris@16:         }
Chris@16: 
Chris@16:         /** Returns true if the two sets of parameters are equal. */
Chris@16:         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
Chris@16:         { return lhs._dim == rhs._dim; }
Chris@16: 
Chris@16:         /** Returns true if the two sets of parameters are different. */
Chris@16:         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
Chris@16: 
Chris@16:     private:
Chris@16:         int _dim;
Chris@16:     };
Chris@16: 
Chris@16:     /**
Chris@16:      * Constructs a @c uniform_on_sphere distribution.
Chris@16:      * @c dim is the dimension of the sphere.
Chris@16:      *
Chris@16:      * Requires: dim >= 0
Chris@16:      */
Chris@16:     explicit uniform_on_sphere(int dim_arg = 2)
Chris@16:       : _container(dim_arg), _dim(dim_arg) { }
Chris@16: 
Chris@16:     /**
Chris@16:      * Constructs a @c uniform_on_sphere distribution from its parameters.
Chris@16:      */
Chris@16:     explicit uniform_on_sphere(const param_type& parm)
Chris@16:       : _container(parm.dim()), _dim(parm.dim()) { }
Chris@16: 
Chris@16:     // compiler-generated copy ctor and assignment operator are fine
Chris@16: 
Chris@16:     /** Returns the dimension of the sphere. */
Chris@16:     int dim() const { return _dim; }
Chris@16: 
Chris@16:     /** Returns the parameters of the distribution. */
Chris@16:     param_type param() const { return param_type(_dim); }
Chris@16:     /** Sets the parameters of the distribution. */
Chris@16:     void param(const param_type& parm)
Chris@16:     {
Chris@16:         _dim = parm.dim();
Chris@16:         _container.resize(_dim);
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns the smallest value that the distribution can produce.
Chris@16:      * Note that this is required to approximate the standard library's
Chris@16:      * requirements.  The behavior is defined according to lexicographical
Chris@16:      * comparison so that for a container type of std::vector,
Chris@16:      * dist.min() <= x <= dist.max() where x is any value produced
Chris@16:      * by the distribution.
Chris@16:      */
Chris@16:     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
Chris@16:     {
Chris@16:         result_type result(_dim);
Chris@16:         if(_dim != 0) {
Chris@16:             result.front() = RealType(-1.0);
Chris@16:         }
Chris@16:         return result;
Chris@16:     }
Chris@16:     /**
Chris@16:      * Returns the largest value that the distribution can produce.
Chris@16:      * Note that this is required to approximate the standard library's
Chris@16:      * requirements.  The behavior is defined according to lexicographical
Chris@16:      * comparison so that for a container type of std::vector,
Chris@16:      * dist.min() <= x <= dist.max() where x is any value produced
Chris@16:      * by the distribution.
Chris@16:      */
Chris@16:     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
Chris@16:     {
Chris@16:         result_type result(_dim);
Chris@16:         if(_dim != 0) {
Chris@16:             result.front() = RealType(1.0);
Chris@16:         }
Chris@16:         return result;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Effects: Subsequent uses of the distribution do not depend
Chris@16:      * on values produced by any engine prior to invoking reset.
Chris@16:      */
Chris@101:     void reset() {}
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns a point uniformly distributed over the surface of
Chris@16:      * a sphere of dimension dim().
Chris@16:      */
Chris@16:     template<class Engine>
Chris@16:     const result_type & operator()(Engine& eng)
Chris@16:     {
Chris@101:         using std::sqrt;
Chris@101:         switch(_dim)
Chris@101:         {
Chris@101:         case 0: break;
Chris@101:         case 1:
Chris@101:             {
Chris@101:                 if(uniform_01<RealType>()(eng) < 0.5) {
Chris@101:                     *_container.begin() = -1;
Chris@101:                 } else {
Chris@101:                     *_container.begin() = 1;
Chris@101:                 }
Chris@101:             }
Chris@101:         case 2:
Chris@101:             {
Chris@101:                 uniform_01<RealType> uniform;
Chris@101:                 RealType sqsum;
Chris@101:                 RealType x, y;
Chris@101:                 do {
Chris@101:                     x = uniform(eng) * 2 - 1;
Chris@101:                     y = uniform(eng) * 2 - 1;
Chris@101:                     sqsum = x*x + y*y;
Chris@101:                 } while(sqsum == 0 || sqsum > 1);
Chris@101:                 RealType mult = 1/sqrt(sqsum);
Chris@101:                 typename Cont::iterator iter = _container.begin();
Chris@101:                 *iter = x * mult;
Chris@101:                 iter++;
Chris@101:                 *iter = y * mult;
Chris@101:                 break;
Chris@101:             }
Chris@101:         case 3:
Chris@101:             {
Chris@101:                 uniform_01<RealType> uniform;
Chris@101:                 RealType sqsum;
Chris@101:                 RealType x, y;
Chris@101:                 do {
Chris@101:                     x = uniform(eng) * 2 - 1;
Chris@101:                     y = uniform(eng) * 2 - 1;
Chris@101:                     sqsum = x*x + y*y;
Chris@101:                 } while(sqsum > 1);
Chris@101:                 RealType mult = 2 * sqrt(1 - sqsum);
Chris@101:                 typename Cont::iterator iter = _container.begin();
Chris@101:                 *iter = x * mult;
Chris@101:                 ++iter;
Chris@101:                 *iter = y * mult;
Chris@101:                 ++iter;
Chris@101:                 *iter = 2 * sqsum - 1;
Chris@101:                 break;
Chris@101:             }
Chris@101:         default:
Chris@101:             {
Chris@101:                 detail::unit_normal_distribution<RealType> normal;
Chris@101:                 RealType sqsum;
Chris@101:                 do {
Chris@101:                     sqsum = 0;
Chris@101:                     for(typename Cont::iterator it = _container.begin();
Chris@101:                         it != _container.end();
Chris@101:                         ++it) {
Chris@101:                         RealType val = normal(eng);
Chris@101:                         *it = val;
Chris@101:                         sqsum += val * val;
Chris@101:                     }
Chris@101:                 } while(sqsum == 0);
Chris@101:                 // for all i: result[i] /= sqrt(sqsum)
Chris@101:                 std::transform(_container.begin(), _container.end(), _container.begin(),
Chris@101:                                std::bind2nd(std::multiplies<RealType>(), 1/sqrt(sqsum)));
Chris@101:             }
Chris@16:         }
Chris@16:         return _container;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns a point uniformly distributed over the surface of
Chris@16:      * a sphere of dimension param.dim().
Chris@16:      */
Chris@16:     template<class Engine>
Chris@16:     result_type operator()(Engine& eng, const param_type& parm) const
Chris@16:     {
Chris@16:         return uniform_on_sphere(parm)(eng);
Chris@16:     }
Chris@16: 
Chris@16:     /** Writes the distribution to a @c std::ostream. */
Chris@16:     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd)
Chris@16:     {
Chris@16:         os << sd._dim;
Chris@16:         return os;
Chris@16:     }
Chris@16: 
Chris@16:     /** Reads the distribution from a @c std::istream. */
Chris@16:     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd)
Chris@16:     {
Chris@16:         is >> sd._dim;
Chris@16:         sd._container.resize(sd._dim);
Chris@16:         return is;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns true if the two distributions will produce identical
Chris@16:      * sequences of values, given equal generators.
Chris@16:      */
Chris@16:     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
Chris@101:     { return lhs._dim == rhs._dim; }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns true if the two distributions may produce different
Chris@16:      * sequences of values, given equal generators.
Chris@16:      */
Chris@16:     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)
Chris@16: 
Chris@16: private:
Chris@16:     result_type _container;
Chris@16:     int _dim;
Chris@16: };
Chris@16: 
Chris@16: } // namespace random
Chris@16: 
Chris@16: using random::uniform_on_sphere;
Chris@16: 
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP