1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2004 Ferdinando Ametrano
5 
6  This file is part of QuantLib, a free-software/open-source library
7  for financial quantitative analysts and developers - http://quantlib.org/
8 
9  QuantLib is free software: you can redistribute it and/or modify it
10  under the terms of the QuantLib license.  You should have received a
11  copy of the license along with this program; if not, please email
12  <quantlib-dev@lists.sf.net>. The license is also available online at
13  <http://quantlib.org/license.shtml>.
14 
15  This program is distributed in the hope that it will be useful, but WITHOUT
16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17  FOR A PARTICULAR PURPOSE.  See the license for more details.
18 */
19 
20 /*! \file randomizedlds.hpp
21     \brief Randomized low-discrepancy sequence
22 */
23 
24 #ifndef quantlib_randomized_lds_hpp
25 #define quantlib_randomized_lds_hpp
26 
27 #include <ql/math/randomnumbers/randomsequencegenerator.hpp>
28 #include <ql/math/randomnumbers/mt19937uniformrng.hpp>
29 
30 namespace QuantLib {
31 
32     //! Randomized (random shift) low-discrepancy sequence
33     /*! Random-shifts a uniform low-discrepancy sequence of dimension
34         \f$ N \f$ by adding (modulo 1 for each coordinate) a pseudo-random
35         uniform deviate in \f$ (0, 1)^N. \f$
36         It is used for implementing Randomized Quasi Monte Carlo.
37 
38         The uniform low discrepancy sequence is supplied by LDS; the
39         uniform pseudo-random sequence is supplied by PRS.
40 
41         Both class LDS and PRS must implement the following interface:
42         \code
43             LDS::sample_type LDS::nextSequence() const;
44             Size LDS::dimension() const;
45         \endcode
46 
47         \pre LDS and PRS must have the same dimension \f$ N \f$
48 
49         \warning Inverting LDS and PRS is possible, but it doesn't
50                  make sense.
51 
52         \todo implement the other randomization algorithms
53 
54         \test correct initialization is tested.
55     */
56     template <class LDS,
57               class PRS = RandomSequenceGenerator<MersenneTwisterUniformRng> >
58     class RandomizedLDS {
59       public:
60         typedef Sample<std::vector<Real> > sample_type;
61         RandomizedLDS(const LDS& ldsg,
62                       const PRS& prsg);
63         RandomizedLDS(const LDS& ldsg);
64         RandomizedLDS(Size dimensionality,
65                       BigNatural ldsSeed = 0,
66                       BigNatural prsSeed = 0);
67         //! returns next sample using a given randomizing vector
68         const sample_type& nextSequence() const;
lastSequence() const69         const sample_type& lastSequence() const {
70             return x;
71         }
72         /*! update the randomizing vector and re-initialize
73             the low discrepancy generator */
nextRandomizer()74         void nextRandomizer() {
75             randomizer_ = prsg_.nextSequence();
76             ldsg_ = pristineldsg_;
77         }
dimension() const78         Size dimension() const {return dimension_;}
79       private:
80         mutable LDS ldsg_, pristineldsg_; // mutable because nextSequence is const
81         PRS prsg_;
82         Size dimension_;
83         mutable sample_type x, randomizer_;
84     };
85 
86     template <class LDS, class PRS>
RandomizedLDS(const LDS & ldsg,const PRS & prsg)87     RandomizedLDS<LDS, PRS>::RandomizedLDS(const LDS& ldsg, const PRS& prsg)
88     : ldsg_(ldsg), pristineldsg_(ldsg),
89       prsg_(prsg), dimension_(ldsg_.dimension()),
90       x(std::vector<Real> (dimension_), 1.0), randomizer_(std::vector<Real> (dimension_), 1.0) {
91 
92         QL_REQUIRE(prsg_.dimension()==dimension_,
93                    "generator mismatch: "
94                    << dimension_ << "-dim low discrepancy "
95                    << "and " << prsg_.dimension() << "-dim pseudo random")
96 
97         randomizer_ = prsg_.nextSequence();
98 
99     }
100 
101     template <class LDS, class PRS>
RandomizedLDS(const LDS & ldsg)102     RandomizedLDS<LDS, PRS>::RandomizedLDS(const LDS& ldsg)
103     : ldsg_(ldsg), pristineldsg_(ldsg),
104       prsg_(ldsg_.dimension()), dimension_(ldsg_.dimension()),
105       x(std::vector<Real> (dimension_), 1.0), randomizer_(std::vector<Real> (dimension_), 1.0) {
106 
107         randomizer_ = prsg_.nextSequence();
108 
109     }
110 
111     template <class LDS, class PRS>
RandomizedLDS(Size dimensionality,BigNatural ldsSeed,BigNatural prsSeed)112     RandomizedLDS<LDS, PRS>::RandomizedLDS(Size dimensionality,
113                                            BigNatural ldsSeed,
114                                            BigNatural prsSeed)
115     : ldsg_(dimensionality, ldsSeed), pristineldsg_(dimensionality, ldsSeed),
116       prsg_(dimensionality, prsSeed), dimension_(dimensionality),
117       x(std::vector<Real> (dimensionality), 1.0), randomizer_(std::vector<Real> (dimensionality), 1.0) {
118 
119         randomizer_ = prsg_.nextSequence();
120     }
121 
122     template <class LDS, class PRS>
123     inline const typename RandomizedLDS<LDS, PRS>::sample_type&
nextSequence() const124     RandomizedLDS<LDS, PRS>::nextSequence() const {
125     typename LDS::sample_type sample =
126         ldsg_.nextSequence();
127     x.weight = randomizer_.weight * sample.weight;
128     for (Size i = 0; i < dimension_; i++) {
129         x.value[i] =  randomizer_.value[i] + sample.value[i];
130         if (x.value[i]>1.0)
131             x.value[i] -= 1.0;
132     }
133     return x;
134     }
135 
136 }
137 
138 
139 #endif
140