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

/*

 Copyright (C) 2009 Ferdinando Ametrano

 Copyright (C) 2007 Allen Kuo

 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.

*/

#include <ql/cashflows/cashflows.hpp>

#include <ql/math/optimization/constraint.hpp>

#include <ql/math/optimization/costfunction.hpp>

#include <ql/math/optimization/simplex.hpp>

#include <ql/pricingengines/bond/bondfunctions.hpp>

#include <ql/termstructures/yield/fittedbonddiscountcurve.hpp>

#include <ql/time/daycounters/simpledaycounter.hpp>

#include <ql/utilities/dataformatters.hpp>

#include <utility>

using std::vector;

namespace QuantLib {

    class FittedBondDiscountCurve::FittingMethod::FittingCost

        : public CostFunction {

        friend class FittedBondDiscountCurve::FittingMethod;

      public:

        explicit FittingCost(

                       FittedBondDiscountCurve::FittingMethod* fittingMethod);

        Real value(const Array& x) const override;

        Array values(const Array& x) const override;

      private:

        FittedBondDiscountCurve::FittingMethod* fittingMethod_;

    };

    FittedBondDiscountCurve::FittedBondDiscountCurve(

        Natural settlementDays,

        const Calendar& calendar,

        vector<ext::shared_ptr<BondHelper> > bondHelpers,

        const DayCounter& dayCounter,

        const FittingMethod& fittingMethod,

        Real accuracy,

        Size maxEvaluations,

        Array guess,

        Real simplexLambda,

        Size maxStationaryStateIterations)

    : YieldTermStructure(settlementDays, calendar, dayCounter), accuracy_(accuracy),

      maxEvaluations_(maxEvaluations), simplexLambda_(simplexLambda),

      maxStationaryStateIterations_(maxStationaryStateIterations), guessSolution_(std::move(guess)),

      bondHelpers_(std::move(bondHelpers)), fittingMethod_(fittingMethod) {

        fittingMethod_->curve_ = this;

        setup();

    }

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(unsigned int, QuantLib::Calendar const&, std::__1::vector<boost::shared_ptr<QuantLib::BondHelper>, std::__1::allocator<boost::shared_ptr<QuantLib::BondHelper> > >, QuantLib::DayCounter const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, double, unsigned long, QuantLib::Array, double, unsigned long)

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(unsigned int, QuantLib::Calendar const&, std::__1::vector<boost::shared_ptr<QuantLib::BondHelper>, std::__1::allocator<boost::shared_ptr<QuantLib::BondHelper> > >, QuantLib::DayCounter const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, double, unsigned long, QuantLib::Array, double, unsigned long)

    FittedBondDiscountCurve::FittedBondDiscountCurve(

        const Date& referenceDate,

        vector<ext::shared_ptr<BondHelper> > bondHelpers,

        const DayCounter& dayCounter,

        const FittingMethod& fittingMethod,

        Real accuracy,

        Size maxEvaluations,

        Array guess,

        Real simplexLambda,

        Size maxStationaryStateIterations)

    : YieldTermStructure(referenceDate, Calendar(), dayCounter), accuracy_(accuracy),

      maxEvaluations_(maxEvaluations), simplexLambda_(simplexLambda),

      maxStationaryStateIterations_(maxStationaryStateIterations), guessSolution_(std::move(guess)),

      bondHelpers_(std::move(bondHelpers)), fittingMethod_(fittingMethod) {

        fittingMethod_->curve_ = this;

        setup();

    }

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(QuantLib::Date const&, std::__1::vector<boost::shared_ptr<QuantLib::BondHelper>, std::__1::allocator<boost::shared_ptr<QuantLib::BondHelper> > >, QuantLib::DayCounter const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, double, unsigned long, QuantLib::Array, double, unsigned long)

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(QuantLib::Date const&, std::__1::vector<boost::shared_ptr<QuantLib::BondHelper>, std::__1::allocator<boost::shared_ptr<QuantLib::BondHelper> > >, QuantLib::DayCounter const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, double, unsigned long, QuantLib::Array, double, unsigned long)

    FittedBondDiscountCurve::FittedBondDiscountCurve(

                            Natural settlementDays,

                            const Calendar& calendar,

                            const FittingMethod& fittingMethod,

                            Array parameters,

                            Date maxDate,

                            const DayCounter& dayCounter)

    : YieldTermStructure(settlementDays, calendar, dayCounter), accuracy_(1e-10),

      maxEvaluations_(0), guessSolution_(std::move(parameters)),

      maxDate_(maxDate), fittingMethod_(fittingMethod) {

        fittingMethod_->curve_ = this;

        setup();

    }

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(unsigned int, QuantLib::Calendar const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, QuantLib::Array, QuantLib::Date, QuantLib::DayCounter const&)

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(unsigned int, QuantLib::Calendar const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, QuantLib::Array, QuantLib::Date, QuantLib::DayCounter const&)

    FittedBondDiscountCurve::FittedBondDiscountCurve(

                                const Date& referenceDate,

                                const FittingMethod& fittingMethod,

                                Array parameters,

                                Date maxDate,

                                const DayCounter& dayCounter)

    : YieldTermStructure(referenceDate, Calendar(), dayCounter), accuracy_(1e-10),

      maxEvaluations_(0), guessSolution_(std::move(parameters)),

      maxDate_(maxDate), fittingMethod_(fittingMethod) {

        fittingMethod_->curve_ = this;

        setup();

    }

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(QuantLib::Date const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, QuantLib::Array, QuantLib::Date, QuantLib::DayCounter const&)

Unexecuted instantiation: QuantLib::FittedBondDiscountCurve::FittedBondDiscountCurve(QuantLib::Date const&, QuantLib::FittedBondDiscountCurve::FittingMethod const&, QuantLib::Array, QuantLib::Date, QuantLib::DayCounter const&)

    void FittedBondDiscountCurve::resetGuess(const Array& guess) {

        QL_REQUIRE(guess.empty() || guess.size() == fittingMethod_->size(), "guess is of wrong size");

        guessSolution_ = guess;

        update();

    }

    void FittedBondDiscountCurve::performCalculations() const {

        if (maxEvaluations_!= 0) {

            // we need to fit, so we require helpers

            QL_REQUIRE(!bondHelpers_.empty(), "no bond helpers given");

        }

        if (maxEvaluations_ == 0) {

            // no fit, but we need either an explicit max date or

            // helpers from which to deduce it

            QL_REQUIRE(maxDate_ != Date() || !bondHelpers_.empty(),

                       "no bond helpers or max date given");

        }

        if (!bondHelpers_.empty()) {

            maxDate_ = Date::minDate();

            Date refDate = referenceDate();

            // double check bond quotes still valid and/or instruments not expired

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

                ext::shared_ptr<Bond> bond = bondHelpers_[i]->bond();

                QL_REQUIRE(bondHelpers_[i]->quote()->isValid(),

                           io::ordinal(i+1) << " bond (maturity: " <<

                           bond->maturityDate() << ") has an invalid price quote");

                Date bondSettlement = bond->settlementDate();

                QL_REQUIRE(bondSettlement>=refDate,

                           io::ordinal(i+1) << " bond settlemente date (" <<

                           bondSettlement << ") before curve reference date (" <<

                           refDate << ")");

                QL_REQUIRE(BondFunctions::isTradable(*bond, bondSettlement),

                           io::ordinal(i+1) << " bond non tradable at " <<

                           bondSettlement << " settlement date (maturity"

                           " being " << bond->maturityDate() << ")");

                maxDate_ = std::max(maxDate_, bondHelpers_[i]->pillarDate());

                bondHelpers_[i]->setTermStructure(

                                                  const_cast<FittedBondDiscountCurve*>(this));

            }

        }

        fittingMethod_->init();

        fittingMethod_->calculate();

    }

    FittedBondDiscountCurve::FittingMethod::FittingMethod(

        bool constrainAtZero,

        const Array& weights,

        ext::shared_ptr<OptimizationMethod> optimizationMethod,

        Array l2,

        const Real minCutoffTime,

        const Real maxCutoffTime,

        Constraint constraint)

    : constrainAtZero_(constrainAtZero), weights_(weights), l2_(std::move(l2)),

      calculateWeights_(weights.empty()), optimizationMethod_(std::move(optimizationMethod)),

      constraint_(std::move(constraint)),

      minCutoffTime_(minCutoffTime), maxCutoffTime_(maxCutoffTime) {

        if (constraint_.empty())

            constraint_ = NoConstraint();

    }

    void FittedBondDiscountCurve::FittingMethod::init() {

        if (curve_->maxEvaluations_ == 0)

            return; // we can skip the rest

        // yield conventions

        DayCounter yieldDC = curve_->dayCounter();

        Compounding yieldComp = Compounded;

        Frequency yieldFreq = Annual;

        Size n = curve_->bondHelpers_.size();

        costFunction_ = ext::make_shared<FittingCost>(this);

        for (auto& bondHelper : curve_->bondHelpers_) {

            bondHelper->setTermStructure(curve_);

        }

        if (calculateWeights_) {

            if (weights_.empty())

                weights_ = Array(n);

            Real squaredSum = 0.0;

            for (Size i=0; i<curve_->bondHelpers_.size(); ++i) {

                ext::shared_ptr<Bond> bond = curve_->bondHelpers_[i]->bond();

                Real amount = curve_->bondHelpers_[i]->quote()->value();

                Bond::Price price(amount, curve_->bondHelpers_[i]->priceType());

                Date bondSettlement = bond->settlementDate();

                Rate ytm = BondFunctions::yield(*bond, price,

                                                yieldDC, yieldComp, yieldFreq,

                                                bondSettlement);

                Time dur = BondFunctions::duration(*bond, ytm,

                                                   yieldDC, yieldComp, yieldFreq,

                                                   Duration::Modified,

                                                   bondSettlement);

                weights_[i] = 1.0/dur;

                squaredSum += weights_[i]*weights_[i];

            }

            weights_ /= std::sqrt(squaredSum);

        }

        QL_REQUIRE(weights_.size() == n,

                   "Given weights do not cover all boostrapping helpers");

        if (!l2_.empty()) {

            QL_REQUIRE(l2_.size() == size(),

                       "Given penalty factors do not cover all parameters");

            QL_REQUIRE(!curve_->guessSolution_.empty(), "L2 penalty requires a guess");

        }

    }

    void FittedBondDiscountCurve::FittingMethod::calculate() {

        if (curve_->maxEvaluations_ == 0)

        {

            // Don't calculate, simply use the given parameters to

            // provide a fitted curve.  This turns the instance into

            // an evaluator of the parametric curve, for example

            // allowing to use the parameters for a credit spread

            // curve calculated with bonds in one currency to be

            // coupled to a discount curve in another currency.

            QL_REQUIRE(curve_->guessSolution_.size() == size(),

                       "wrong number of parameters");

            solution_ = curve_->guessSolution_;

            numberOfIterations_ = 0;

            costValue_ = Null<Real>();

            errorCode_ = EndCriteria::None;

            return;

        }

        FittingCost& costFunction = *costFunction_;

        // start with the guess solution, if it exists

        Array x(size(), 0.0);

        if (!curve_->guessSolution_.empty()) {

            QL_REQUIRE(curve_->guessSolution_.size() == size(), "wrong size for guess");

            x = curve_->guessSolution_;

        }

        // workaround for backwards compatibility

        ext::shared_ptr<OptimizationMethod> optimization = optimizationMethod_;

        if (!optimization) {

            optimization = ext::make_shared<Simplex>(curve_->simplexLambda_);

        }

        Problem problem(costFunction, constraint_, x);

        Real rootEpsilon = curve_->accuracy_;

        Real functionEpsilon =  curve_->accuracy_;

        Real gradientNormEpsilon = curve_->accuracy_;

        EndCriteria endCriteria(curve_->maxEvaluations_,

                                curve_->maxStationaryStateIterations_,

                                rootEpsilon,

                                functionEpsilon,

                                gradientNormEpsilon);

        errorCode_ = optimization->minimize(problem,endCriteria);

        solution_ = problem.currentValue();

        numberOfIterations_ = problem.functionEvaluation();

        costValue_ = problem.functionValue();

        // save the results as the guess solution, in case of recalculation

        curve_->guessSolution_ = solution_;

    }

    FittedBondDiscountCurve::FittingMethod::FittingCost::FittingCost(

                        FittedBondDiscountCurve::FittingMethod* fittingMethod)

    : fittingMethod_(fittingMethod) {}

    Real FittedBondDiscountCurve::FittingMethod::FittingCost::value(

                                                       const Array& x) const {

        Real squaredError = 0.0;

        Array vals = values(x);

        for (Real val : vals) {

            squaredError += val;

        }

        return squaredError;

    }

    Array FittedBondDiscountCurve::FittingMethod::FittingCost::values(const Array &x) const {

        Size n = fittingMethod_->curve_->bondHelpers_.size();

        Size N = fittingMethod_->l2_.size();

        // set solution so that fittingMethod_->curve_ represents the current trial

        // the final solution will be set in FittingMethod::calculate() later on

        fittingMethod_->solution_ = x;

        Array values(n + N);

        for (Size i=0; i<n; ++i) {

            ext::shared_ptr<BondHelper> helper = fittingMethod_->curve_->bondHelpers_[i];

            Real error = helper->impliedQuote() - helper->quote()->value();

            Real weightedError = fittingMethod_->weights_[i] * error;

            values[i] = weightedError * weightedError;

        }

        if (N != 0) {

            for (Size i = 0; i < N; ++i) {

                Real error = x[i] - fittingMethod_->curve_->guessSolution_[i];

                values[i + n] = fittingMethod_->l2_[i] * error * error;

            }

        }

        return values;

    }

}