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

/*
 Copyright (C) 2006 Ferdinando Ametrano
 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb

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

/*! \file simplex.hpp
    \brief Simplex optimization method
*/

/* The implementation of the algorithm was inspired by
 * "Numerical Recipes in C", 2nd edition, Press, Teukolsky, Vetterling, Flannery
 * Chapter 10
 */

#ifndef quantlib_optimization_simplex_hpp
#define quantlib_optimization_simplex_hpp

#include <ql/math/optimization/problem.hpp>
#include <vector>

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

#if !defined(__GNUC__) || __GNUC__ > 3 || __GNUC_MINOR__ > 4
#define QL_ARRAY_EXPRESSIONS
#endif

namespace QuantLib {

    //! Multi-dimensional simplex class
    /*! This method is rather raw and requires quite a lot of
        computing resources, but it has the advantage that it does not
        need any evaluation of the cost function's gradient, and that
        it is quite easily implemented. First, we choose N+1
        starting points, given here by a starting point \f$
        \mathbf{P}_{0} \f$ and N points such that
        \f[ 
            \mathbf{P}_{\mathbf{i}}=\mathbf{P}_{0}+\lambda \mathbf{e}_{\mathbf{i}},
        \f]
        where \f$ \lambda \f$ is the problem's characteristic length scale). These 
        points will form a geometrical form called simplex.
        The principle of the downhill simplex method is, at each
        iteration, to move the worst point (highest cost function value)
        through the opposite face to a better point. When the simplex
        seems to be constrained in a valley, it will be contracted
        downhill, keeping the best point unchanged.

        \ingroup optimizers
    */
    class Simplex : public OptimizationMethod {
      public:
        /*! Constructor taking as input the characteristic length */
        Simplex(Real lambda) : lambda_(lambda) {}
        virtual EndCriteria::Type minimize(Problem& P,
                                           const EndCriteria& endCriteria);
      private:
        Real extrapolate(Problem& P,
                         Size iHighest,
                         Real& factor) const;
        Real lambda_;
        mutable std::vector<Array> vertices_;
        mutable Array values_, sum_;
    };

    // implementation

    namespace {
    // Computes the size of the simplex
    inline Real computeSimplexSize (const std::vector<Array>& vertices) {
            Array center(vertices.front().size(),0);
            for (Size i=0; i<vertices.size(); ++i)
                center += vertices[i];
            center *=1/Real(vertices.size());
            Real result = 0;
            for (Size i=0; i<vertices.size(); ++i) {
                Array temp =  vertices[i] - center;
                result += std::sqrt(DotProduct(temp,temp));
            }
            return result/Real(vertices.size());
        }
    }

    inline Real Simplex::extrapolate(Problem& P,
                              Size iHighest,
                              Real &factor) const {

        Array pTry;
        do {
            Size dimensions = values_.size() - 1;
            Real factor1 = (1.0 - factor)/dimensions;
            Real factor2 = factor1 - factor;
            #if defined(QL_ARRAY_EXPRESSIONS)
            pTry = sum_*factor1 - vertices_[iHighest]*factor2;
            #else
            // composite expressions fail to compile with gcc 3.4 on windows
            pTry = sum_*factor1;
            pTry -= vertices_[iHighest]*factor2;
            #endif
            factor *= 0.5;
        } while (!P.constraint().test(pTry) && std::fabs(factor) > QL_EPSILON);
        if (std::fabs(factor) <= QL_EPSILON) {
            return values_[iHighest];
        }
        factor *= 2.0;
        Real vTry = P.value(pTry);
        if (vTry < values_[iHighest]) {
            values_[iHighest] = vTry;
            #if defined(QL_ARRAY_EXPRESSIONS)
            sum_ += pTry - vertices_[iHighest];
            #else
            sum_ += pTry;
            sum_ -= vertices_[iHighest];
            #endif
            vertices_[iHighest] = pTry;
        }
        return vTry;

    }


    inline EndCriteria::Type Simplex::minimize(Problem& P,
                                        const EndCriteria& endCriteria) {
        // set up of the problem
        //Real ftol = endCriteria.functionEpsilon();    // end criteria on f(x) (see Numerical Recipes in C++, p.410)
        Real xtol = endCriteria.rootEpsilon();          // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/)
        Size maxStationaryStateIterations_
            = endCriteria.maxStationaryStateIterations();
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        Array x_ = P.currentValue();
        Integer iterationNumber_=0;

        // Initialize vertices of the simplex
        bool end = false;
        Size n = x_.size(), i;
        vertices_ = std::vector<Array>(n+1, x_);
        for (i=0; i<n; i++) {
            Array direction(n, 0.0);
            direction[i] = 1.0;
            P.constraint().update(vertices_[i+1], direction, lambda_);
        }
        // Initialize function values at the vertices of the simplex
        values_ = Array(n+1, 0.0);
        for (i=0; i<=n; i++)
            values_[i] = P.value(vertices_[i]);
        // Loop looking for minimum
        do {
            sum_ = Array(n, 0.0);
            Size i;
            for (i=0; i<=n; i++)
                sum_ += vertices_[i];
            // Determine the best (iLowest), worst (iHighest)
            // and 2nd worst (iNextHighest) vertices
            Size iLowest = 0;
            Size iHighest, iNextHighest;
            if (values_[0]<values_[1]) {
                iHighest = 1;
                iNextHighest = 0;
            } else {
                iHighest = 0;
                iNextHighest = 1;
            }
            for (i=1;i<=n; i++) {
                if (values_[i]>values_[iHighest]) {
                    iNextHighest = iHighest;
                    iHighest = i;
                } else {
                    if ((values_[i]>values_[iNextHighest]) && i!=iHighest)
                        iNextHighest = i;
                }
                if (values_[i]<values_[iLowest])
                    iLowest = i;
            }
            // Now compute accuracy, update iteration number and check end criteria
            //// Numerical Recipes exit strategy on fx (see NR in C++, p.410)
            //Real low = values_[iLowest];
            //Real high = values_[iHighest];
            //Real rtol = 2.0*std::fabs(high - low)/
            //    (std::fabs(high) + std::fabs(low) + QL_EPSILON);
            //++iterationNumber_;
            //if (rtol < ftol ||
            //    endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
            // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl
            Real simplexSize = computeSimplexSize(vertices_);
            ++iterationNumber_;
            if (simplexSize < xtol ||
                endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
                endCriteria.checkStationaryPoint(0.0, 0.0,
                    maxStationaryStateIterations_, ecType);
                endCriteria.checkMaxIterations(iterationNumber_, ecType);
                x_ = vertices_[iLowest];
                Real low = values_[iLowest];
                P.setFunctionValue(low);
                P.setCurrentValue(x_);
                return ecType;
            }
            // If end criteria is not met, continue
            Real factor = -1.0;
            Real vTry = extrapolate(P, iHighest, factor);
            if ((vTry <= values_[iLowest]) && (factor == -1.0)) {
                factor = 2.0;
                extrapolate(P, iHighest, factor);
            } else if (std::fabs(factor) > QL_EPSILON) {
                if (vTry >= values_[iNextHighest]) {
                    Real vSave = values_[iHighest];
                    factor = 0.5;
                    vTry = extrapolate(P, iHighest, factor);
                    if (vTry >= vSave && std::fabs(factor) > QL_EPSILON) {
                        for (Size i=0; i<=n; i++) {
                            if (i!=iLowest) {
                                #if defined(QL_ARRAY_EXPRESSIONS)
                                vertices_[i] =
                                    0.5*(vertices_[i] + vertices_[iLowest]);
                                #else
                                vertices_[i] += vertices_[iLowest];
                                vertices_[i] *= 0.5;
                                #endif
                                values_[i] = P.value(vertices_[i]);
                            }
                        }
                    }
                }
            }
            // If can't extrapolate given the constraints, exit
            if (std::fabs(factor) <= QL_EPSILON) {
                x_ = vertices_[iLowest];
                Real low = values_[iLowest];
                P.setFunctionValue(low);
                P.setCurrentValue(x_);
                return EndCriteria::StationaryFunctionValue;
            }
        } while (end == false);
        QL_FAIL("optimization failed: unexpected behaviour");
    }


}

#endif
