////////////////////////////////////////////////////////////////////////
/// \class RAT::Optimisers::SimulatedAnnealing
///
/// \brief  Optimisation using simulated annealing
///
/// \author Matt Mottram <m.mottram@qmul.ac.uk> -- contact person
///
/// REVISION HISTORY:\n
///     26/03/2014 : M Mottram - first revision \n
///
/// \details Optimisation using simulated annealing (from Numerical Recipes)
///         Annealing schedule parameters are not necessarily obvious to
///         choose; the defaults here are from the SNOMAN version of this method.
///
////////////////////////////////////////////////////////////////////////

#ifndef __RAT_Optimiser_SimulatedAnnealing__
#define __RAT_Optimiser_SimulatedAnnealing__

#include <string>
#include <vector>
#include <math.h>

#include <RAT/Optimiser.hh>

namespace RAT
{

namespace Optimisers
{

class SimulatedAnnealing : public Optimiser
{
public:
  virtual std::string
  GetName() const { return SimulatedAnnealing::Name(); }

  static std::string
  Name() { return std::string( "simulatedAnnealing" ); }

  void Initialise(const std::string& init );

  void BeginOfRun( DS::Run& ) { }

  void EndOfRun( DS::Run& ) { }

  virtual void SetI( const std::string& param, const int value );

  virtual void SetD( const std::string& param, const double value );

  virtual double Minimise();

  virtual double Maximise();

protected:

  /// Runs the SimulatedAnnealing optimisation according to the global annealing schedule parameters
  ///
  /// @return Result from the PDF for the optimised best fit vertex parameters
  double AnnealOptimise();

  /// Simulated annealing routine from Numerical Recipes: amebsa
  ///
  /// @param[in,out] p vector of (N+1)*N points defining N+1 points in N dimensional space
  /// @param[in,out] y vector of N+1 values of the calculated PDF values for the N+1 points
  /// @param[in] ndim number of dimensions (N)
  /// @param[in,out] pb vector of N points describing the current best fit point
  /// @param[in,out] yb calculated PDF value for the best fit points
  /// @param[in] ftol fractional tolerance (difference) between best and worst PDF values for the given simplex
  /// @param[in,out] iter number of iterations to step through, forced return when iter becomes negative
  /// @param[in] temptr temperature for the given stage in the annealing schedule
  void Amebsa(std::vector<double>& p, std::vector<double>& y, int ndim, std::vector<double>& pb,
              double& yb, double ftol, int& iter, double temptr);

  /// Simulated annealing routine from Numerical Recipes: amotsa
  ///
  /// Extrapolates from across the face of the simplex from a high point, saves the new point
  /// if better than the original high point
  ///
  /// @param[in] p vector of (N+1)*N points defining N+1 points in N dimensional space
  /// @param[in,out] y vector of N+1 values of the calculated PDF values for the N+1 points
  /// @param[in,out] psum the sums (over N+1 values) of the values for each dimension in the simplex
  /// @param[in] ndim number of dimensions (N)
  /// @param[in,out] pb vector of N points describing the current best fit point
  /// @param[in,out] yb calculated PDF value for the best fit point
  /// @param[in] ihi index of the worst point in the simplex
  /// @param[in,out] yhi calculated PDF value for the worst fit point
  /// @param[in] fac Extrapolation factor through the simplex
  /// @param[in,out] tt thermal fluctuation
  /// @return PDF value for the new (or old) point with an additional thermal fluctuation added, then removed
  double Amotsa(std::vector<double>& p, std::vector<double>& y, std::vector<double>& psum, int ndim,
                std::vector<double>& pb, double& yb, int ihi, double& yhi, double fac, const double& tt);

  /// Sums the values of each dimension for N+1 points in the N-dimensional parameter space
  ///
  /// @param[in,out] psum the sums (over N+1 values) of the values for each dimension in the simplex
  /// @param[in] ndim the number of dimensions (N)
  /// @param[in] p vector of (N+1)*N points defining N+1 points in N dimensional space
  void GetPSum(std::vector<double>& psum, const int& ndim, const std::vector<double>& p);

  int fIterations; /// <Number of iterations (simplex resize) at a given temperature
  double fTolerance; /// <fractional tolerance at each temperature between highest and lowest point on simplex
  int fScheduleSteps; /// <Number of steps in method (i.e. number of temperatures)
  int fSchedulePower; /// <Power, a, in T = T0 * (1 - i/N)^a schedule
  double fStartTemperature; /// <Starting temperature, defaults to Ybest - Yworst from starting simplex
  double fMinFactor; /// <Minimise or maximise

};
}
}

#endif // __RAT_Optimiser_SimulatedAnnealing__