SobolGenerator.hh

Go to the documentation of this file.

#ifndef NPSTAT_SOBOLGENERATOR_HH_
#define NPSTAT_SOBOLGENERATOR_HH_
 
/*!
// \file SobolGenerator.hh
//
// \brief Generator of Sobol low-discrepancy sequences
//
// Author: I. Volobouev
//
// October 2010
*/
 
#include "npstat/rng/AbsRandomGenerator.hh"
 
namespace npstat {
    /**
    // Generator of Sobol sequences. Based on the "i8_sobol" function
    // written by John Burkardt:
    // http://people.sc.fsu.edu/~jburkardt/cpp_src/sobol/sobol.html
    */
    class SobolGenerator : public AbsRandomGenerator
    {
    public:
        /**
        // Maximum number of dimensions and maximum number of points which
        // can be generated by this code (the latter equals pow(2,LOG_MAX))
        */
        enum {
            DIM_MAX = 1111,
            LOG_MAX = 62
        };
 
        /**
        // Constructor arguments are as follows:
        //
        // dim           -- dimensionality of the point set
        //
        // maxPowerOfTwo -- maximum number of points which can be made
        //                  by this generator will be 2^maxPowerOfTwo.
        //                  maxPowerOfTwo must not be less than 2 and
        //                  must not exceed LOG_MAX.
        //
        // nSkip         -- number of events to skip at the beginning
        //                  in order to "warm up" the generator.
        //                  Reasonable values of this argument are
        //                  2^(maxPowerOfTwo-k), with some small k > 0.
        */
        explicit SobolGenerator(unsigned dim, unsigned maxPowerOfTwo=LOG_MAX,
                                unsigned nSkip=0U);
        inline virtual ~SobolGenerator() {}
 
        inline virtual unsigned dim() const {return dim_num;}
        virtual void run(double* buf, unsigned bufSize, unsigned nPoints);
 
        double operator()();
        inline unsigned long long maxPoints() const {return 1ULL << maxcol;}
 
    protected:
        void nextBitSet(long long* buf, unsigned bufSize);
        inline unsigned maxPowerOfTwo() const {return maxcol;}
 
    private:
        void update();
 
        double recipd_;
        long long count;
        long long lastq[DIM_MAX];
        long long v[DIM_MAX][LOG_MAX];
        unsigned dim_num;
        unsigned maxcol;
    };
}
 
#endif // NPSTAT_SOBOLGENERATOR_HH_