//****************************************************************************

//! \file kutest.hpp
//! \date 23-10-2003

//! \class vWalker
//! \brief The class representing all the little Vector Walkers from Kurchans Method

//! \class vMult
//! \brief unary function object for multiplying a vector with a number

//! \class vAddRand
//! \brief unary function object for adding a random number to a vector

//! \class vCalc
//! \brief unary function object for evaluating a func pointer vector

//! \class Tkutest
//! \brief the Main class of the Simulation

//*****************************************************************************

//     This file is part of the kutest.
//
//     The kutest is free software; you can redistribute it and/or modify
//     it under the terms of the GNU General Public License as published by
//     the Free Software Foundation; either version 2 of the License, or
//     (at your option) any later version.
//
//     The kutest 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
//     GNU General Public License for more details.
//
//     You should have received a copy of the GNU General Public License
//     along with DataView; if not, write to the Free Software
//     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//*****************************************************************************

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <vector>
#include <list>
#include <algorithm>
#include <functional>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <cmath>

#ifndef KUTEST_HPP
#define KUTEST_HPP

class Tkutest;
//! typedef for the function pointer for the ejf and egf functions
typedef double (Tkutest::*jgfunc)( const std::vector<double> x);

//*****************************************************************************
class vWalker
//-----------------------------------------------------------------------------
//! This class is a container for the vWalkers Positions and vectors providing
//! the functionality the vWalkers need
//*****************************************************************************
{

	bool cl;//!< true if this vWalker is a clone
	bool h;//!< Yuhuu an occasion for a var named h

public:
	std::vector<double> x;//!< the "Position" of the vWalker
	std::vector<double> v;//!< the vector of the vWalker

	//! \brief the copy constructor (needed for vWalker vcloning)
	//! \par in the vWalker to be cloned
	vWalker( const vWalker *in) : cl( true)
	{
		x = *( new std::vector<double>( in->x));
		v = *( new std::vector<double>( in->v));
	}

	//! \brief the constructor
	//! \par inx the "Position" of the vWalker
	//! \par inv the vector of the vWalker
	vWalker( const std::vector<double> inx, const std::vector<double> inv) : cl( false)
	{
		std::copy( inx.begin(), inx.end(), std::back_inserter( x));
		std::copy( inv.begin(), inv.end(), std::back_inserter( v));
	}

	//! find out if this vWalker is a clone an then rehabilitate him
	const bool clonetest()
	{
		h = cl;
		cl = false;
		return h;
	}
};

//*****************************************************************************
template <class Arg>
class vMult:  public std::unary_function<Arg, void>
//-----------------------------------------------------------------------------
//! A unary function object derived from unary_function to multiply a vector
//! with a number
//*****************************************************************************
{
	Arg multiplier;//!< the multiplier
public:
	//! \brief A constructor does nothing but to initialize the multiplier
	//! \param x the multiplier
	vMult(const Arg &x) : multiplier( x) { };

	//! \brief Another constructor doing nothing (if you want to initialize the multiplier later)
	vMult() : multiplier( 0){};

	//! \brief resets the multiplier
	//! \param x the new multiplier
	void mreset( const Arg &x) { multiplier = x;}

	//! \brief the function doing the real work
	//! \param x the value to be multiplied with the multiplier
	void operator()( Arg &x)
	{
		x*=multiplier;
	}
};

//*****************************************************************************
template <class Arg>
class vAddRand:  public std::unary_function<Arg, void>
//-----------------------------------------------------------------------------
//! A unary function object derived from unary_function to add a random
//! number to a vector
//*****************************************************************************
{
	double bound;//! the variance for the random number distribution
	gsl_rng *gslrand1;//! the gsl random number generator

public:
	//! \brief A constructor does nothing but to initialize the boundary
	//! \param g1 a gsl random number generator
	//! \param x the boundary
	vAddRand( gsl_rng *g1, const double &x) : gslrand1(g1), bound( x) {};

	//! \brief Another constructor doing nothing (if you want to initialize the boundary later)
	vAddRand() : bound( 0.){};

	~vAddRand() {};

	//! \brief resets the boundary
	//! \param x the new boundary
	void mreset( const double &x) { bound = x;}

	//! \brief the function doing the real work
	//! \param x the value to be added to the random number
	void operator()( Arg &x)
	{
		x+=gsl_ran_gaussian( gslrand1, bound);//rand()*2.*bound/RAND_MAX - bound;
	}
};

//*****************************************************************************
class vCalc :  public std::unary_function<jgfunc, double>
//-----------------------------------------------------------------------------
//! A unary function object derived from unary_function to evaluate a
//! function pointer value
//*****************************************************************************
{
	std::vector<double> a;//!< the argument
	Tkutest *Instance;
public:
	//! \brief The constructor does nothing but to initialize the Instance of Kutest for
	//! the function pointer evalualtion
	vCalc( Tkutest *in) : Instance( in){};

	//! \brief resets the function pointer parameter
	//! \param x the new parameter
	void mreset( const std::vector<double> &x) { a = x;}

	//! \brief the function doing the real work
	//! \param x the function pointer to evaluated
	double operator()(const jgfunc &x) const
	{
		return (Instance->*x)( a);
	}
};

//*****************************************************************************
class Tkutest
//-----------------------------------------------------------------------------
//! The Simulation class; overload when you implement a different Energy function
//*****************************************************************************
{
	double fos;//!< helper variable used in the EJF and EGF calculations
	double bx;//!< buffer for x used in the EJF and EGF calculations
	double by;//! <buffer for y used in the EJF and EGF calculations

	// The Energy Jacobi Matrix Functions

	//! \brief EJF Matrix function component 1 1
	//! \param x the position of the walker in Phase Space
	inline double ejf11( const std::vector<double> x);
	//! \brief EJF Matrix function component 1 2
	//! \param x the position of the walker in Phase Space
	inline double ejf12( const std::vector<double> x);
	//! \brief EJF Matrix function component 2 1
	//! \param x the position of the walker in Phase Space
	inline double ejf21( const std::vector<double> x);
	//! \brief EJF Matrix function component 2 2
	//! \param x the position of the walker in Phase Space
	inline double ejf22( const std::vector<double> x);

	// The Energy Gradient Functions

	//! \brief EGF Vector function component 1
	//! \param x the position of the walker in Phase Space
	inline double egf1( const std::vector<double> x);
	//! \brief EGF Vector function component 2
	//! \param x the position of the walker in Phase Space
	inline double egf2( const std::vector<double> x);

protected:
	int dim;//!< the Dimension of the Simulation
	double dT;//!< the timestep
	double bound;//!< the boundaries for the Sim
	bool inited;//!< if init has allready been called
	vCalc *calc;//!< instance of the unary_function for calculating the EG and the EJ
	vMult<double> *mult;//!< instance of the unary_function to multiply a vector with dT (in this case)
	vMult<double> *amult;//!< instance of the unary_function to multiply a vector with A (in this case)
	vAddRand<double> *ar;//!< instance of the unary_function to add a random number to a vector

	gsl_rng *gslrand1; //!< the gsl random number generator (currently used just for the Gauss Dist)

	std::vector<jgfunc> egf;//!< the Energy Gradient function vector
	std::vector<double> eg;//!< the Energy Gradient vector
	std::vector<std::vector<jgfunc> > ejf;//!< the Energy Jacobi function Matrix
	std::vector<std::vector<double> > ej;//!< the Energy Jacobi Matrix

	std::vector<jgfunc>::iterator egfit;//!< the Energy Gradient function vector iterator
	std::vector<double>::iterator egit;//!< the Energy Gradient vector Iterator
	std::vector<std::vector<jgfunc> >::iterator ejfit;//!< the Energy Jacobi function Matrix iterator
	std::vector<std::vector<double> >::iterator ejit;//!< the Energy Jacobi Matrix iterator

	std::list<vWalker> vvec;//!< vector containing the vWalkers of the sim
	std::list<vWalker>::iterator vvecit;//!< iterator for the vWalkers

	std::vector<double>::iterator vit;//!< simple vector iterators allways can come in handy

	std::ofstream ofile;//!< the output file

	//! From this ejf matrix later the numerical matrix will be calculated for every walker!
	//! \attention Overwrite when you implement a different Energy function
	//! \brief initializes the The Energy Jacobi Function (ejf) matrix
	virtual void SetJFMatrix();

	//! From this egf vector later the numerical vector will be calculated for every walker!
	//! \attention Overwrite when you implement a different Energy function
	//! \brief initializes the The Energy Gradient Function (egf) vector
	virtual void SetGFVector();

	//! \brief this sets the initial vWalker distribution
	virtual void SetStartCfg( int numvw = 500);

	//! \brief move of a kurchan step for one vWalker ( Position Displacement)
	//! \param walk the walker to be treated
	inline void move( vWalker &walk);

	//! \brief rotation of a kurchan step for one vWalker (Vector Rotation)
	//! \param walk the walker to be treated
	//! \param ev the directional vector times the EJ Matrix
	//! \param a the a value
	inline void rotate( vWalker &walk, const std::vector<double> &ev, const double a);

	//! \brief calculate the Energy Jacoby Matrix from the function Matrix
	//! \param x the Position of the vWalker
	inline void calcEMatrix( std::vector<double> x);

public:

	//! \brief the constructor
	Tkutest();
	//! \brief the destructor
	~Tkutest();

	//! \brief the main working function doing all the necessary stuff
	//! \attention Use Init before starting the Sim
	//! \param max number of steps to be carried out
	//! \return true if it worked false if any kind of error appeared
	bool Execute( int max = 1);

	//! \brief Initializes everything necessary for the Sim
	//! \attention Use this function before starting the Sim
	//! \param inbound the boundaries for the Phase Space
	//! \param numvw number of initial vWalkers
	//! \param indT the Timestep
	//! \param indis the random displacement variance (= 2 times the Temperature)
	void Init( double inbound, int numvw = 500, double indT = 0.001, double indis = 0.25);

	//! \brief write out the current phase Space
	void write();
};

#endif //KUTEST_HPP