// Copyright 2005, Ondra Kamenik // Quasi Monte Carlo quadrature. /* This defines quasi Monte Carlo quadratures for cube and for a function multiplied by normal density. The quadrature for a cube is named |QMCarloCubeQuadrature| and integrates: $$\int_{\langle 0,1\rangle^n}f(x){\rm d}x$$ The quadrature for a function of normally distributed parameters is named |QMCarloNormalQuadrature| and integrates: $${1\over\sqrt{(2\pi)^n}}\int_{(-\infty,\infty)^n}f(x)e^{-{1\over 2}x^Tx}{\rm d}x$$ For a cube we define |qmcpit| as iterator of |QMCarloCubeQuadrature|, and for the normal density multiplied function we define |qmcnpit| as iterator of |QMCarloNormalQuadrature|. The quasi Monte Carlo method generates low discrepancy points with equal weights. The one dimensional low discrepancy sequences are generated by |RadicalInverse| class, the sequences are combined for higher dimensions by |HaltonSequence| class. The Halton sequence can use a permutation scheme; |PermutattionScheme| is an abstract class for all permutaton schemes. We have three implementations: |WarnockPerScheme|, |ReversePerScheme|, and |IdentityPerScheme|. */ #ifndef QUASI_MCARLO_H #define QUASI_MCARLO_H #include "int_sequence.hh" #include "quadrature.hh" #include "Vector.hh" #include /* This abstract class declares |permute| method which permutes coefficient |c| having index of |i| fro the base |base| and returns the permuted coefficient which must be in $0,\ldots,base-1$. */ class PermutationScheme { public: PermutationScheme() { } virtual ~PermutationScheme() { } virtual int permute(int i, int base, int c) const = 0; }; /* This class represents an integer number |num| as $c_0+c_1b+c_2b^2+\ldots+c_jb^j$, where $b$ is |base| and $c_0,\ldots,c_j$ is stored in |coeff|. The size of |IntSequence| coeff does not grow with growing |num|, but is fixed from the very beginning and is set according to supplied maximum |maxn|. The basic method is |eval| which evaluates the |RadicalInverse| with a given permutation scheme and returns the point, and |increase| which increases |num| and recalculates the coefficients. */ class RadicalInverse { int num; int base; int maxn; int j; IntSequence coeff; public: RadicalInverse(int n, int b, int mxn); RadicalInverse(const RadicalInverse &ri) : num(ri.num), base(ri.base), maxn(ri.maxn), j(ri.j), coeff(ri.coeff) { } const RadicalInverse & operator=(const RadicalInverse &radi) { num = radi.num; base = radi.base; maxn = radi.maxn; j = radi.j; coeff = radi.coeff; return *this; } double eval(const PermutationScheme &p) const; void increase(); void print() const; }; /* This is a vector of |RadicalInverse|s, each |RadicalInverse| has a different prime as its base. The static members |primes| and |num_primes| define a precalculated array of primes. The |increase| method of the class increases indices in all |RadicalInverse|s and sets point |pt| to contain the points in each dimension. */ class HaltonSequence { private: static int primes[]; static int num_primes; protected: int num; int maxn; vector ri; const PermutationScheme &per; Vector pt; public: HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p); HaltonSequence(const HaltonSequence &hs) : num(hs.num), maxn(hs.maxn), ri(hs.ri), per(hs.per), pt(hs.pt) { } const HaltonSequence &operator=(const HaltonSequence &hs); void increase(); const Vector & point() const { return pt; } const int getNum() const { return num; } void print() const; protected: void eval(); }; /* This is a specification of quasi Monte Carlo quadrature. It consists of dimension |dim|, number of points (or level) |lev|, and the permutation scheme. This class is common to all quasi Monte Carlo classes. */ class QMCSpecification { protected: int dim; int lev; const PermutationScheme &per_scheme; public: QMCSpecification(int d, int l, const PermutationScheme &p) : dim(d), lev(l), per_scheme(p) { } virtual ~QMCSpecification() { } int dimen() const { return dim; } int level() const { return lev; } const PermutationScheme & getPerScheme() const { return per_scheme; } }; /* This is an iterator for quasi Monte Carlo over a cube |QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of the same dimension as given by the specification. An iterator can be constructed from a given number |n|, or by a copy constructor. For technical reasons, there is also an empty constructor; for that reason, every member is a pointer. */ class qmcpit { protected: const QMCSpecification *spec; HaltonSequence *halton; ParameterSignal *sig; public: qmcpit(); qmcpit(const QMCSpecification &s, int n); qmcpit(const qmcpit &qpit); ~qmcpit(); bool operator==(const qmcpit &qpit) const; bool operator!=(const qmcpit &qpit) const { return !operator==(qpit); } const qmcpit &operator=(const qmcpit &qpit); qmcpit &operator++(); const ParameterSignal & signal() const { return *sig; } const Vector & point() const { return halton->point(); } double weight() const; void print() const { halton->print(); } }; /* This is an easy declaration of quasi Monte Carlo quadrature for a cube. Everything important has been done in its iterator |qmcpit|, so we only inherit from general |Quadrature| and reimplement |begin| and |numEvals|. */ class QMCarloCubeQuadrature : public QuadratureImpl, public QMCSpecification { public: QMCarloCubeQuadrature(int d, int l, const PermutationScheme &p) : QuadratureImpl(d), QMCSpecification(d, l, p) { } virtual ~QMCarloCubeQuadrature() { } int numEvals(int l) const { return l; } protected: qmcpit begin(int ti, int tn, int lev) const { return qmcpit(*this, ti*level()/tn + 1); } }; /* This is an iterator for |QMCarloNormalQuadrature|. It is equivalent to an iterator for quasi Monte Carlo cube quadrature but a point. The point is obtained from a point of |QMCarloCubeQuadrature| by a transformation $\langle 0,1\rangle\rightarrow\langle-\infty,\infty\rangle$ aplied to all dimensions. The transformation yields a normal distribution from a uniform distribution on $\langle0,1\rangle$. It is in fact |NormalICDF|. */ class qmcnpit : public qmcpit { protected: Vector *pnt; public: qmcnpit(); qmcnpit(const QMCSpecification &spec, int n); qmcnpit(const qmcnpit &qpit); ~qmcnpit(); bool operator==(const qmcnpit &qpit) const { return qmcpit::operator==(qpit); } bool operator!=(const qmcnpit &qpit) const { return !operator==(qpit); } const qmcnpit &operator=(const qmcnpit &qpit); qmcnpit &operator++(); const ParameterSignal & signal() const { return *sig; } const Vector & point() const { return *pnt; } void print() const { halton->print(); pnt->print(); } }; /* This is an easy declaration of quasi Monte Carlo quadrature for a cube. Everything important has been done in its iterator |qmcnpit|, so we only inherit from general |Quadrature| and reimplement |begin| and |numEvals|. */ class QMCarloNormalQuadrature : public QuadratureImpl, public QMCSpecification { public: QMCarloNormalQuadrature(int d, int l, const PermutationScheme &p) : QuadratureImpl(d), QMCSpecification(d, l, p) { } virtual ~QMCarloNormalQuadrature() { } int numEvals(int l) const { return l; } protected: qmcnpit begin(int ti, int tn, int lev) const { return qmcnpit(*this, ti*level()/tn + 1); } }; /* Declares Warnock permutation scheme. */ class WarnockPerScheme : public PermutationScheme { public: int permute(int i, int base, int c) const; }; /* Declares reverse permutation scheme. */ class ReversePerScheme : public PermutationScheme { public: int permute(int i, int base, int c) const; }; /* Declares no permutation (identity) scheme. */ class IdentityPerScheme : public PermutationScheme { public: int permute(int i, int base, int c) const { return c; } }; #endif