/***************************************************************************
 *                           RK4_GPU2.h                                   *
 *                           -------------------                           *
 * copyright            : (C) 2013 by Francisco Naveros                    *
 * email                : fnaveros@ugr.es                                  *
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program 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 3 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef RK4_GPU2_H_
#define RK4_GPU2_H_

/*!
 * \file RK4_GPU2.h
 *
 * \author Francisco Naveros
 * \date May 2013
 *
 * This file declares a class which implement a fixed step fourth order Runge-Kutta integration method in GPU (this class is stored
 * in GPU memory and executed in GPU. 
 */

#include "./FixedStep_GPU2.h"

#include "../../include/neuron_model/TimeDrivenNeuronModel_GPU2.h"

//Library for CUDA
#include <helper_cuda.h>

/*!
 * \class RK4_GPU2
 *
 * \brief Fourth order Runge-Kutta integration method in GPU.
 * 
 * This class abstracts the behavior of a Euler integration method for neurons in a 
 * time-driven spiking neural network.
 * It includes internal model functions which define the behavior of integration methods
 * (initialization, calculate next value, ...).
 *
 * \author Francisco Naveros
 * \date May 2012
 */
class RK4_GPU2 : public FixedStep_GPU2 {
	public:

		/*!
		 * \brief These vectors are used as auxiliar vectors.
		*/
		float * AuxNeuronState;
		float * AuxNeuronState1;
		float * AuxNeuronState2;
		float * AuxNeuronState3;
		float * AuxNeuronState4;


		/*!
		 * \brief Constructor of the class with 5 parameter.
		 *
		 * It generates a new fourth order Runge-Kutta object in GPU memory.
		 *
		 * \param N_neuronStateVariables Number of state variables for each cell.
		 * \param N_differentialNeuronState Number of state variables which are calculate with a differential equation for each cell.
		 * \param N_timeDependentNeuronState Number of state variables which ara calculate with a time dependent equation for each cell.
		 * \param Total_N_thread Number of thread in GPU (in this method it is not necessary)
		 * \param Buffer_GPU This vector contains all the necesary GPU memory which have been reserved in the CPU (this memory
		 * could be reserved directly in the GPU, but this suppose some restriction in the amount of memory which can be reserved).
		 */
		__device__ RK4_GPU2(TimeDrivenNeuronModel_GPU2* NewModel, int N_neuronStateVariables, int N_differentialNeuronState, int N_timeDependentNeuronState, void ** Buffer_GPU):FixedStep_GPU2(NewModel, N_neuronStateVariables, N_differentialNeuronState, N_timeDependentNeuronState, Buffer_GPU){
			AuxNeuronState=((float*)Buffer_GPU[1]);
			AuxNeuronState1=((float*)Buffer_GPU[2]);
			AuxNeuronState2=((float*)Buffer_GPU[3]);
			AuxNeuronState3=((float*)Buffer_GPU[4]);
			AuxNeuronState4=((float*)Buffer_GPU[5]);
		}

		/*!
		 * \brief Class destructor.
		 *
		 * It destroys an object of this class.
		 */
		__device__ ~RK4_GPU2(){
		}
		


		/*!
		 * \brief It calculate the next neural state varaibles of the model.
		 *
		 * It calculate the next neural state varaibles of the model.
		 *
		 * \param index Index of the cell inside the neuron model for method with memory (e.g. BDF).
		 * \param SizeStates Number of neurons
		 * \param Model The NeuronModel.
		 * \param NeuronState Vector of neuron state variables for all neurons.
		 */
		__device__ bool NextDifferentialEcuationValues(int index, int SizeStates, float * NeuronState){
			float previous_V=NeuronState[index];
			int offset1=gridDim.x * blockDim.x;
			int offset2=blockDim.x * blockIdx.x + threadIdx.x;

			//1st term
			model->EvaluateDifferentialEcuation(index, SizeStates, NeuronState, AuxNeuronState1, this->elapsedTimeInNeuronModelScale);
			
			//2nd term
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState[j*offset1 + offset2]= NeuronState[j*SizeStates + index] + AuxNeuronState1[j*offset1 + offset2]*this->elapsedTimeInNeuronModelScale*0.5f;
			}
			for (int j=N_DifferentialNeuronState; j<N_NeuronStateVariables; j++){
				AuxNeuronState[j*offset1 + offset2]=NeuronState[j*SizeStates + index];
			}


			model->EvaluateTimeDependentEcuation(offset2, offset1, AuxNeuronState, this->elapsedTimeInNeuronModelScale*0.5f, 0);
			model->EvaluateDifferentialEcuation(offset2, offset1, AuxNeuronState, AuxNeuronState2, this->elapsedTimeInNeuronModelScale);

			//3rd term
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState[j*offset1 + offset2]= NeuronState[j*SizeStates + index] + AuxNeuronState2[j*offset1 + offset2]*this->elapsedTimeInNeuronModelScale*0.5f;
			}

			model->EvaluateDifferentialEcuation(offset2, offset1, AuxNeuronState, AuxNeuronState3, this->elapsedTimeInNeuronModelScale);

			//4rd term
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState[j*offset1 + offset2]= NeuronState[j*SizeStates + index] + AuxNeuronState3[j*offset1 + offset2]*this->elapsedTimeInNeuronModelScale;
			}

			model->EvaluateTimeDependentEcuation(offset2, offset1, AuxNeuronState, this->elapsedTimeInNeuronModelScale*0.5f, 0);
			model->EvaluateDifferentialEcuation(offset2, offset1, AuxNeuronState, AuxNeuronState4, this->elapsedTimeInNeuronModelScale);


			for (int j=0; j<N_DifferentialNeuronState; j++){
				NeuronState[j*SizeStates + index]+=(AuxNeuronState1[j*offset1 + offset2]+2*(AuxNeuronState2[j*offset1 + offset2]+AuxNeuronState3[j*offset1 + offset2])+AuxNeuronState4[j*offset1 + offset2])*this->elapsedTimeInNeuronModelScale*0.166666666667f;
			}

			//Finaly, we evaluate the neural state variables with time dependence.
			//Model->EvaluateTimeDependentEcuation(index, SizeStates, NeuronState, this->elapsedTimeInNeuronModelScale, 1);
			for (int j=N_DifferentialNeuronState; j<N_NeuronStateVariables; j++){
				NeuronState[j*SizeStates + index]=AuxNeuronState[j*offset1 + offset2];
			}

			//Update the last spike time.
			this->model->vectorNeuronState_GPU2->LastSpikeTimeGPU[index]+=this->elapsedTimeInSeconds;

			bool spike=model->EvaluateSpikeCondition(previous_V, NeuronState, index, this->elapsedTimeInNeuronModelScale);

			return spike;
			
		}




		/*!
		 * \brief It reset the state of the integration method for method with memory (e.g. BDF).
		 *
		 * It reset the state of the integration method for method with memory (e.g. BDF).
		 *
		 * \param index indicate which neuron must be reseted.
		 *
		 */
		__device__ void resetState(int index){
		}

		/*!
		 * \brief It calculates the conductance exponential values for time driven neuron models.
		 *
		 * It calculates the conductance exponential values for time driven neuron models.
		 */
		__device__ virtual void Calculate_conductance_exp_values(){
			this->model->Initialize_conductance_exp_values(this->N_TimeDependentNeuronState,2);
			//index 0
			this->model->Calculate_conductance_exp_values(0, elapsedTimeInNeuronModelScale*0.5f);
		}

};



#endif /* RK4_GPU2_H_ */
