/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*

 Copyright (C) 2014 Jose Aparicio

 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.

*/

#include <ql/experimental/math/convolvedstudentt.hpp>

#include <ql/errors.hpp>

#include <ql/math/factorial.hpp>

#include <ql/math/distributions/normaldistribution.hpp>

#include <ql/math/solvers1d/brent.hpp>

#include <ql/math/functional.hpp>

#include <boost/math/distributions/students_t.hpp>

namespace QuantLib {

    CumulativeBehrensFisher::CumulativeBehrensFisher(const std::vector<Integer>& degreesFreedom,

                                                     const std::vector<Real>& factors)

    : degreesFreedom_(degreesFreedom), factors_(factors), polyConvolved_(std::vector<Real>(1, 1.))

    {

        QL_REQUIRE(degreesFreedom.size() == factors.size(),

            "Incompatible sizes in convolution.");

        for (int i : degreesFreedom) {

            QL_REQUIRE(i % 2 != 0, "Even degree of freedom not allowed");

            QL_REQUIRE(i >= 0, "Negative degree of freedom not allowed");

        }

        for(Size i=0; i<degreesFreedom_.size(); i++)

            polynCharFnc_.push_back(polynCharactT((degreesFreedom[i]-1)/2));

        // adjust the polynomial coefficients by the factors in the linear

        //   combination:

        for(Size i=0; i<degreesFreedom_.size(); i++) {

            Real multiplier = 1.;

            for(Size k=1; k<polynCharFnc_[i].size(); k++) {

                multiplier *= std::abs(factors_[i]);

                polynCharFnc_[i][k] *= multiplier;

            }

        }

        //convolution, here it is a product of polynomials and exponentials

        for (auto& i : polynCharFnc_)

            polyConvolved_ = convolveVectorPolynomials(polyConvolved_, i);

        // trim possible zeros that might have arised:

        auto it = polyConvolved_.rbegin();

        while (it != polyConvolved_.rend()) {

            if (*it == 0.) {

                polyConvolved_.pop_back();

                it = polyConvolved_.rbegin();

              }else{

                  break;

              }

          }

          // cache 'a' value (the exponent)

          for(Size i=0; i<degreesFreedom_.size(); i++)

              a_ += std::sqrt(static_cast<Real>(degreesFreedom_[i]))

                * std::abs(factors_[i]);

          a2_ = a_ * a_;

    }

    std::vector<Real> CumulativeBehrensFisher::polynCharactT(Natural n) const {

        Natural nu = 2 * n +1;

        std::vector<Real> low(1,1.), high(1,1.);

        high.push_back(std::sqrt(static_cast<Real>(nu)));

        if(n==0) return low;

        if(n==1) return high;

        for(Size k=1; k<n; k++) {

            std::vector<Real> recursionFactor(1,0.); // 0 coef

            recursionFactor.push_back(0.); // 1 coef

            recursionFactor.push_back(nu/((2.*k+1.)*(2.*k-1.))); // 2 coef

            std::vector<Real> lowUp =

                convolveVectorPolynomials(recursionFactor, low);

            //add them up:

            for(Size i=0; i<high.size(); i++)

                lowUp[i] += high[i];

            low = high;

            high = lowUp;

        }

        return high;

    }

    std::vector<Real> CumulativeBehrensFisher::convolveVectorPolynomials(

        const std::vector<Real>& v1,

        const std::vector<Real>& v2) const {

    #if defined(QL_EXTRA_SAFETY_CHECKS)

        QL_REQUIRE(!v1.empty() && !v2.empty(),

            "Incorrect vectors in polynomial.");

    #endif

        const std::vector<Real>& shorter = v1.size() < v2.size() ? v1 : v2;

        const std::vector<Real>& longer = (v1 == shorter) ? v2 : v1;

        Size newDegree = v1.size()+v2.size()-2;

        std::vector<Real> resultB(newDegree+1, 0.);

        for(Size polyOrdr=0; polyOrdr<resultB.size(); polyOrdr++) {

            for(Size i=std::max<Integer>(0, polyOrdr-longer.size()+1);

                i<=std::min(polyOrdr, shorter.size()-1); i++)

                resultB[polyOrdr] += shorter[i]*longer[polyOrdr-i];

        }

        return resultB;

    }

    Probability CumulativeBehrensFisher::operator()(const Real x) const {

        // 1st & 0th terms with the table integration

        Real integral = polyConvolved_[0] * std::atan(x/a_);

        Real squared = a2_ + x*x;

        Real rootsqr = std::sqrt(squared);

        Real atan2xa = std::atan2(-x,a_);

        if(polyConvolved_.size()>1)

            integral += polyConvolved_[1] * x/squared;

        for(Size exponent = 2; exponent <polyConvolved_.size(); exponent++) {

            integral -= polyConvolved_[exponent] *

                Factorial::get(exponent-1) * std::sin((exponent)*atan2xa)

                    /std::pow(rootsqr, static_cast<Real>(exponent));

         }

        return .5 + integral / M_PI;

    }

    Probability

    CumulativeBehrensFisher::density(const Real x) const {

        Real squared = a2_ + x*x;

        Real integral = polyConvolved_[0] * a_ / squared;

        Real rootsqr = std::sqrt(squared);

        Real atan2xa = std::atan2(-x,a_);

        for(Size exponent=1; exponent <polyConvolved_.size(); exponent++) {

            integral += polyConvolved_[exponent] *

                Factorial::get(exponent) * std::cos((exponent+1)*atan2xa)

                    /std::pow(rootsqr, static_cast<Real>(exponent+1) );

        }

        return integral / M_PI;

    }

    InverseCumulativeBehrensFisher::InverseCumulativeBehrensFisher(

        const std::vector<Integer>& degreesFreedom,

        const std::vector<Real>& factors,

        Real accuracy)

    : normSqr_(std::inner_product(factors.begin(), factors.end(),

        factors.begin(), Real(0.))),

      accuracy_(accuracy), distrib_(degreesFreedom, factors) { }

    Real InverseCumulativeBehrensFisher::operator()(const Probability q) const {

        Probability effectiveq;

        Real sign;

        // since the distrib is symmetric solve only on the right side:

        if(q==0.5) {

            return 0.;

        }else if(q < 0.5) {

            sign = -1.;

            effectiveq = 1.-q;

        }else{

            sign = 1.;

            effectiveq = q;

        }

        Real xMin =

            InverseCumulativeNormal::standard_value(effectiveq) * normSqr_;

        // inversion will fail at the Brent's bounds-check if this is not enough

        // (q is very close to 1.), in a bad combination fails around 1.-1.e-7

        Real xMax = 1.e6;

        return sign *

            Brent().solve([&](Real x) -> Real { return distrib_(x) - effectiveq; },

                          accuracy_, (xMin+xMax)/2., xMin, xMax);

    }

}