#ifndef __JFIT__JGRADIENT__
#define __JFIT__JGRADIENT__

#include <limits>
#include <vector>
#include <cmath>
#include <memory>
#include <ostream>
#include <iomanip>

#include "JLang/JManip.hh"
#include "JLang/JException.hh"
#include "Jeep/JMessage.hh"


/**
 * \author mdejong
 */

namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {

  /**
   * Auxiliary data structure for fit parameter.
   */
  struct JParameter_t {
    /**
     * Virtual destructor.
     */
    virtual ~JParameter_t()
    {}

    /**
     * Apply step.
     *
     * \param  step              step
     */
    virtual void apply(const double step) = 0;
  };


  /**
   * Auxiliary data structure for editable parameter.
   */
  struct JModifier_t :
    public std::shared_ptr<JParameter_t>
  {
    /**
     * Constructor.
     *
     * \param  name              name
     * \param  parameter         parameter
     * \param  value             value
     */
    JModifier_t(const std::string& name,
		JParameter_t*      parameter,
		const double       value) :
      std::shared_ptr<JParameter_t>(parameter),
      name (name),
      value(value)
    {}

    std::string name;
    double      value;
  };
  

  /**
   * Conjugate gradient fit.
   */
  struct JGradient :
    public std::vector<JModifier_t>
  {
    /**
     * Constructor.
     *
     * The number of iterations and epsilon refer to the number of steps and 
     * the distance to the minimum, respectively.\n
     * The number of extra steps can be used to overcome a possible hurdle on the way.
     *
     * \param  Nmax              maximum number of iterations
     * \param  Nextra            maximum number of extra steps
     * \param  epsilon           epsilon
     * \param  debug             debug
     */
    JGradient(const size_t Nmax     = std::numeric_limits<size_t>::max(),
	      const size_t Nextra   = 0,
	      const double epsilon  = 1.0e-4,
	      const int    debug    = 3) :
      Nmax   (Nmax),
      Nextra (Nextra),
      epsilon(epsilon),
      debug  (debug)
    {}


    /**
     * Fit.
     *
     * The template parameter should provide for the following function operator.
     * <pre>
     *    double operator()(int option);
     * </pre>
     * The value of the option corresponds to the following cases.
     *  - 0 => step wise improvement of the chi2;
     *  - 1 => evaluation of the chi2 before the determination of the gradient of the chi2; and
     *  - 2 => evaluation of the derivative of the chi2 to each fit parameter.
     *
     * \param  getChi2           chi2 function
     * \return                   chi2
     */
    template<class T>
    double operator()(const T& getChi2)
    {
      using namespace std;
      using namespace JPP;

      vector<double> chi2(5, numeric_limits<double>::max());

      if (this->empty()) {
	return numeric_limits<double>::max();
      }

      chi2[0] = this->evaluate(getChi2);
      
      const size_t N = this->size();

      vector<double> H(N);
      vector<double> G(N);

      for (size_t i = 0; i != N; ++i) {
	G[i] = -1.0 * gradient[i];
	H[i] = G[i];
	gradient[i] = H[i];
      }

      numberOfIterations = 0;

      for ( ; numberOfIterations != Nmax; ++numberOfIterations) {

	DEBUG("chi2[0]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);

	// minimise chi2 in direction of gradient

	chi2[1] = chi2[0];
	chi2[2] = chi2[1];

	size_t m = 0;
	
	for (double ds = 1.0; ds > 1.0e-3; ) {

	  this->move(+1.0 * ds);

	  chi2[3] = getChi2(0);

	  DEBUG("chi2[3]  " << setw(4) << m << ' ' << FIXED(12,5) << chi2[3] << ' ' << FIXED(12,5) << ds << endl);

	  if (chi2[3] < chi2[2]) {

	    chi2[1] = chi2[2];
	    chi2[2] = chi2[3];

	    m = 0;
	    
	    continue;
	  }

	  if (ds == 1.0) {

	    if (m == 0) {
	      chi2[4] = chi2[3];
	    }

	    if (m != Nextra) {

	      ++m;

	      continue;

	    } else {

	      for ( ; m != 0; --m) {
		this->move(-1.0 * ds);
	      }

	      chi2[3] = chi2[4];
	    }
	  }

	  this->move(-1.0 * ds);

	  if (chi2[2] < chi2[3]) {

	    // final step based on parabolic interpolation through following points
	    //
	    // x1 = -1 * ds  ->  chi2[1]
	    // x2 =  0 * ds  ->  chi2[2]
	    // x3 = +1 * ds  ->  chi2[3]

	    const double f21 = chi2[2] - chi2[1];   // f(x2) - f(x1)
	    const double f23 = chi2[2] - chi2[3];   // f(x2) - f(x3)

	    const double xs  =  0.5 * (f21 - f23) / (f23 + f21);

	    this->move(+1.0 * xs * ds);

	    chi2[3] = getChi2(0);
	    
	    if (chi2[3] < chi2[2]) {

	      chi2[2] = chi2[3];

	    } else {

	      this->move(-1.0 * xs * ds);

	      chi2[2] = getChi2(0);
	    }

	    DEBUG("chi2[2]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[2] << ' ' << SCIENTIFIC(12,5) << ds << endl);

	    break;

	  } else {

	    ds *= 0.5;
	  }
	}

	if (fabs(chi2[2] - chi2[0]) < epsilon * 0.5 * (fabs(chi2[0]) + fabs(chi2[2]))) {

	  chi2[0] = chi2[2];

	  break;
        }

        chi2[0] = this->evaluate(getChi2);

        double gg  = 0.0;
        double dgg = 0.0;

        for (size_t i = 0; i != N; ++i){
	  gg  +=  G[i]*G[i];
	  dgg += (gradient[i] + G[i]) * gradient[i];
        }

	if (gg == 0.0) {
	  break;
        }

        dgg /= gg;

        for (size_t i = 0; i != N; ++i){
	  G[i] = -1.0 * gradient[i];
	  H[i] =  G[i] + dgg * H[i];
	  gradient[i] = H[i];
        }
      }

      DEBUG("chi2[0]  " << setw(4) << numberOfIterations << ' ' << FIXED(12,5) << chi2[0] << endl);

      return chi2[0];
    }


    size_t  Nmax;              //!< maximum number of iterations
    size_t  Nextra;            //!< maximum number of extra steps
    double  epsilon;           //!< epsilon
    int     debug;             //!< debug

    size_t  numberOfIterations;

  private:
    /**
     * Evaluate gradient.
     *
     * \return                   chi2
     */
    template<class T>
    double evaluate(const T& getChi2)
    {
      using namespace std;
      using namespace JPP;

      const size_t N = this->size();

      gradient.resize(N);

      for (std::vector<double>::iterator i = gradient.begin(); i != gradient.end(); ++i) {
	*i = 0.0;
      }
      
      const double x0 = getChi2(1);

      size_t width = 1;

      for (size_t i = 0; i != N; ++i) {
	if ((*this)[i].name.size() > width) {
	  width = (*this)[i].name.size();
	}
      }

      double V = 0.0;

      for (size_t i = 0; i != N; ++i) {

	if ((*this)[i].value != 0.0) {

	  (*this)[i]->apply(+0.5 * (*this)[i].value);

	  const double x1 = getChi2(2);

	  (*this)[i]->apply(-0.5 * (*this)[i].value);
	  (*this)[i]->apply(-0.5 * (*this)[i].value);

	  const double x2 = getChi2(2);

	  gradient[i] = x1 - x2;
	
	  (*this)[i]->apply(+0.5 * (*this)[i].value);

	  DEBUG(setw(width) << left << (*this)[i].name << right << ' ' << FIXED(12,5) << (*this)[i].value << ' ' << FIXED(12,5) << gradient[i] << endl);
	  
	} else {

	  gradient[i] = 0.0;
	}

	V += gradient[i] * gradient[i];
      }

      DEBUG(setw(width) << left << "|gradient|" << right << ' ' << FIXED(12,5) << sqrt(V) << endl);

      return x0;
    }
    

    /**
     * Move.
     *
     * \param  factor            factor
     */
    void move(const double factor)
    {
      if        (factor > 0.0) {
	for (size_t i = 0; i != this->size(); ++i) {
	  (*this)[ i ]->apply((*this)[ i ].value * gradient[ i ] * factor);
	}
      }	else if (factor < 0.0) {
	for (size_t i = this->size(); i != 0; --i) {
	  (*this)[i-1]->apply((*this)[i-1].value * gradient[i-1] * factor);
	}
      }
    }
    
    std::vector<double> gradient;
  };
}

#endif