1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2015 Andres Hernandez 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 isotropicrandomwalk.hpp 21 \brief Isotropic random walk 22 */ 23 24 #ifndef quantlib_isotropic_random_walk_hpp 25 #define quantlib_isotropic_random_walk_hpp 26 27 #include <ql/mathconstants.hpp> 28 #include <ql/math/randomnumbers/mt19937uniformrng.hpp> 29 #include <ql/math/array.hpp> 30 #include <boost/random/variate_generator.hpp> 31 32 namespace QuantLib { 33 34 //! Isotropic random walk 35 /*! A variate is used to draw from a random element of a 36 probability distribution. The draw corresponds to the 37 radius of a d-dimensional sphere. The position on the 38 surface of the d-dimensional sphere is randomly chosen 39 with all points on the surface having the same probability, 40 i.e. all directions are isotropic and the step is randomly 41 drawn from the given variate. 42 */ 43 template <class Distribution, class Engine> 44 class IsotropicRandomWalk { 45 public: 46 typedef boost::variate_generator<Engine, Distribution> VariateGenerator; IsotropicRandomWalk(const Engine & eng,Distribution dist,Size dim,const Array & weights=Array (),unsigned long seed=0)47 IsotropicRandomWalk(const Engine& eng, Distribution dist, Size dim, 48 const Array& weights = Array(), 49 unsigned long seed = 0) : 50 variate_(eng, dist), rng_(seed), 51 weights_(weights), dim_(dim) { 52 if (weights_.empty()) 53 weights_ = Array(dim, 1.0); 54 else 55 QL_REQUIRE(dim_ == weights_.size(), "Invalid weights"); 56 } 57 template <class InputIterator> nextReal(InputIterator first) const58 inline void nextReal(InputIterator first) const { 59 Real radius = variate_(); 60 Array::const_iterator weight = weights_.begin(); 61 if (dim_ > 1) { 62 //Isotropic random direction 63 Real phi = M_PI*rng_.nextReal(); 64 for (Size i = 0; i < dim_ - 2; i++) { 65 *first++ = radius*cos(phi)*(*weight++); 66 radius *= sin(phi); 67 phi = M_PI*rng_.nextReal(); 68 } 69 *first++ = radius*cos(2.0*phi)*(*weight++); 70 *first = radius*sin(2.0*phi)*(*weight); 71 } 72 else { 73 if (rng_.nextReal() < 0.5) 74 *first = -radius*(*weight); 75 else 76 *first = radius*(*weight); 77 } 78 } setDimension(Size dim)79 inline void setDimension(Size dim) { 80 dim_ = dim; 81 weights_ = Array(dim, 1.0); 82 } setDimension(Size dim,const Array & weights)83 inline void setDimension(Size dim, const Array& weights) { 84 QL_REQUIRE(dim == weights.size(), "Invalid weights"); 85 dim_ = dim; 86 weights_ = weights; 87 } 88 /*! 89 The isotropic random walk will not adjust its draw to be within the lower and upper bounds, 90 but if the limits are provided, they are used to rescale the sphere so as to make it to an 91 ellipsoid, with different radius in different dimensions. 92 */ setDimension(Size dim,const Array & lowerBound,const Array & upperBound)93 inline void setDimension(Size dim, 94 const Array& lowerBound, const Array& upperBound) { 95 QL_REQUIRE(dim == lowerBound.size(), 96 "Incompatible dimension and lower bound"); 97 QL_REQUIRE(dim == upperBound.size(), 98 "Incompatible dimension and upper bound"); 99 //Find largest bound 100 Array bounds = upperBound - lowerBound; 101 Real maxBound = bounds[0]; 102 for (Size j = 1; j < dim; j++) { 103 if (bounds[j] > maxBound) maxBound = bounds[j]; 104 } 105 //weights by dimension is the size of the bound 106 //divided by the largest bound 107 maxBound = 1.0 / maxBound; 108 bounds *= maxBound; 109 setDimension(dim, bounds); 110 } 111 protected: 112 mutable VariateGenerator variate_; 113 MersenneTwisterUniformRng rng_; 114 Array weights_; 115 Size dim_; 116 }; 117 } 118 #endif 119