/* -*- 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

 <https://www.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.

*/

/*! \file isotropicrandomwalk.hpp

    \brief Isotropic random walk

*/

#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 {

    //! Isotropic random walk

    /*! A variate is used to draw from a random element of a 

        probability distribution. The draw corresponds to the 

        radius of a d-dimensional sphere. The position on the

        surface of the d-dimensional sphere is randomly chosen

        with all points on the surface having the same probability,

        i.e. all directions are isotropic and the step is randomly

        drawn from the given variate.

    */

    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>

        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);

            }

        }

Unexecuted instantiation: void QuantLib::IsotropicRandomWalk<std::__1::normal_distribution<double>, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::nextReal<double*>(double*)

Unexecuted instantiation: void QuantLib::IsotropicRandomWalk<QuantLib::LevyFlightDistribution, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::nextReal<double*>(double*)

        void setDimension(Size dim) { 

            dim_ = dim;

            weights_ = Array(dim, 1.0);

        }

        void setDimension(Size dim, const Array& weights) {

            QL_REQUIRE(dim == weights.size(), "Invalid weights");

            dim_ = dim;

            weights_ = weights;

        }

Unexecuted instantiation: QuantLib::IsotropicRandomWalk<std::__1::normal_distribution<double>, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::setDimension(unsigned long, QuantLib::Array const&)

Unexecuted instantiation: QuantLib::IsotropicRandomWalk<QuantLib::LevyFlightDistribution, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::setDimension(unsigned long, QuantLib::Array const&)

        /*!

        The isotropic random walk will not adjust its draw to be within the lower and upper bounds,

        but if the limits are provided, they are used to rescale the sphere so as to make it to an

        ellipsoid, with different radius in different dimensions.

        */

        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);

        }

Unexecuted instantiation: QuantLib::IsotropicRandomWalk<std::__1::normal_distribution<double>, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::setDimension(unsigned long, QuantLib::Array const&, QuantLib::Array const&)

Unexecuted instantiation: QuantLib::IsotropicRandomWalk<QuantLib::LevyFlightDistribution, std::__1::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >::setDimension(unsigned long, QuantLib::Array const&, QuantLib::Array const&)

      protected:

        Engine engine_;

        Distribution distribution_;

        MersenneTwisterUniformRng rng_;

        Array weights_;

        Size dim_;

    };

}

#endif