isotropicrandomwalk.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2015 Andres Hernandez
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/
 
#ifndef quantlib_isotropic_random_walk_hpp
#define quantlib_isotropic_random_walk_hpp
 
#include <ql/math/array.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/mathconstants.hpp>
#include <utility>
 
namespace QuantLib {
 
 
    template <class Distribution, class Engine>
    class IsotropicRandomWalk {
      public:
        IsotropicRandomWalk(Engine eng,
                            Distribution dist,
                            Size dim,
                            Array weights = Array(),
                            unsigned long seed = 0)
        : engine_(std::move(eng)), distribution_(std::move(dist)), rng_(seed), weights_(std::move(weights)), dim_(dim) {
            if (weights_.empty())
                weights_ = Array(dim, 1.0);
            else
                QL_REQUIRE(dim_ == weights_.size(), "Invalid weights");
        }
        template <class InputIterator>
        inline void nextReal(InputIterator first) {
            Real radius = distribution_(engine_);
            Array::const_iterator weight = weights_.begin();
            if (dim_ > 1) {
                //Isotropic random direction
                Real phi = M_PI*rng_.nextReal();
                for (Size i = 0; i < dim_ - 2; i++) {
                    *first++ = radius*cos(phi)*(*weight++);
                    radius *= sin(phi);
                    phi = M_PI*rng_.nextReal();
                }
                *first++ = radius*cos(2.0*phi)*(*weight++);
                *first = radius*sin(2.0*phi)*(*weight);
            }
            else {
                if (rng_.nextReal() < 0.5)
                    *first = -radius*(*weight);
                else
                    *first = radius*(*weight);
            }
        }
        inline void setDimension(Size dim) { 
            dim_ = dim;
            weights_ = Array(dim, 1.0);
        }
        inline void setDimension(Size dim, const Array& weights) {
            QL_REQUIRE(dim == weights.size(), "Invalid weights");
            dim_ = dim;
            weights_ = weights;
        }
        inline void setDimension(Size dim,
            const Array& lowerBound, const Array& upperBound) {
            QL_REQUIRE(dim == lowerBound.size(),
                "Incompatible dimension and lower bound");
            QL_REQUIRE(dim == upperBound.size(),
                "Incompatible dimension and upper bound");
            //Find largest bound
            Array bounds = upperBound - lowerBound;
            Real maxBound = bounds[0];
            for (Size j = 1; j < dim; j++) {
                if (bounds[j] > maxBound) maxBound = bounds[j];
            }
            //weights by dimension is the size of the bound
            //divided by the largest bound
            maxBound = 1.0 / maxBound;
            bounds *= maxBound;
            setDimension(dim, bounds);
        }
      protected:
        Engine engine_;
        Distribution distribution_;
        MersenneTwisterUniformRng rng_;
        Array weights_;
        Size dim_;
    };
}
#endif