/***************************************************************************
 *                          IzhikevichTimeDrivenModel_GPU2.h                 *
 *                           -------------------                           *
 * copyright            : (C) 2012 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 HHTIMEDRIVENMODEL_GPU2_H_
#define HHTIMEDRIVENMODEL_GPU2_H_

/*!
 * \file LIFTimeDrivenModel_GPU.h
 *
 * \author Francisco Naveros
 * \date November 2012
 *
 * This file declares a class which abstracts a Leaky Integrate-And-Fire neuron model with one 
 * differential equation and two time dependent equations (conductances). This model is
 * implemented in GPU.
 */

#include "./TimeDrivenNeuronModel_GPU2.h"
#include "../integration_method/IntegrationMethod_GPU2.h"
#include "../integration_method/LoadIntegrationMethod_GPU2.h"

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

/*!
 * \class LIFTimeDrivenModel_1_2_GPU
 *
 *
 * \brief Leaky Integrate-And-Fire Time-Driven neuron model with a membrane potential and
 * two conductances. This model is implemented in GPU.
 *
 * This class abstracts the behavior of a neuron in a time-driven spiking neural network.
 * It includes internal model functions which define the behavior of the model
 * (initialization, update of the state, synapses effect, next firing prediction...).
 * This is only a virtual function (an interface) which defines the functions of the
 * inherited classes.
 *
 * \author Francisco Naveros
 * \date May 2013
 */

class HHTimeDrivenModel_GPU2 : public TimeDrivenNeuronModel_GPU2 {
	public:
		/*!
		 * \brief Excitatory reversal potential
		 */
		float eexc;

		/*!
		 * \brief Inhibitory reversal potential
		 */
		float einh;

		/*!
		 * \brief Resting potential
		 */
		float erest;

		/*!
		 * \brief Resting conductance
		 */
		float grest;

		/*!
		 * \brief Inverse membrane capacitance
		 */
		float inv_cm;

		/*!
		 * \brief Membrane capacitance
		 */
		float cm;

		/*!
		 * \brief Firing threshold
		 */
		float vthr;

		/*!
		 * \brief AMPA receptor time constant
		 */
		float texc;
		float inv_texc;

		/*!
		 * \brief GABA receptor time constant
		 */
		float tinh;
		float inv_tinh;


		//////////////////
		float gNa, gKd, ENa, EK, VT,area, inv_area;




		/*!
		 * \brief Number of state variables for each cell.
		*/
		static const int N_NeuronStateVariables=6;

		/*!
		 * \brief Number of state variables which are calculate with a differential equation for each cell.
		*/
		static const int N_DifferentialNeuronState=4;

		/*!
		 * \brief Number of state variables which are calculate with a time dependent equation for each cell.
		*/
		static const int N_TimeDependentNeuronState=2;

		
		/*!
		 * \brief constructor with parameters.
		 *
		 * It generates a new neuron model object.
		 *
		 * \param Eexc eexc.
		 * \param Einh einh.
		 * \param Erest erest.
		 * \param Vthr vthr.
		 * \param Cm cm.
		 * \param Texc texc.
		 * \param Tinh tinh.
		 * \param Tref tref.
		 * \param Grest grest.
		 * \param integrationName integration method type.
		 * \param N_neurons number of neurons.
		 * \param Total_N_thread total number of CUDA thread.
		 * \param Buffer_GPU Gpu auxiliar memory.
		 *
		 */
		__device__ HHTimeDrivenModel_GPU2(float new_eexc, float new_einh, float new_erest, float new_grest, float new_cm, float new_vthr, float new_texc, 
				float new_tinh, float new_gNa, float new_gKd, float new_ENa, float new_EK, float new_VT,
				char const* integrationName, int N_neurons, void ** Buffer_GPU):TimeDrivenNeuronModel_GPU2(MilisecondScale_GPU),
				eexc(new_eexc), einh(new_einh), erest(new_erest), grest(new_grest), cm(new_cm), inv_cm(1.0f/new_cm), vthr(new_vthr), texc(new_texc), 
				inv_texc(1.0f/new_texc), tinh(new_tinh), inv_tinh(1.0f/new_tinh), gNa(new_gNa),	gKd(new_gKd), ENa(new_ENa), EK(new_EK),
				VT(new_VT){
			integrationMethod_GPU2=LoadIntegrationMethod_GPU2::loadIntegrationMethod_GPU2(this, integrationName, N_NeuronStateVariables, N_DifferentialNeuronState, N_TimeDependentNeuronState, Buffer_GPU);
		
			integrationMethod_GPU2->Calculate_conductance_exp_values();	
		}

		/*!
		 * \brief Class destructor.
		 *
		 * It destroys an object of this class.
		 */
		__device__ virtual ~HHTimeDrivenModel_GPU2(){
			delete integrationMethod_GPU2;
		}


		/*!
		 * \brief Update the neuron state variables.
		 *
		 * It updates the neuron state variables.
		 *
		 * \param index The cell index inside the StateGPU. 
		 * \param AuxStateGPU Auxiliary incremental conductance vector.
		 * \param StateGPU Neural state variables.
		 * \param LastUpdateGPU Last update time
		 * \param LastSpikeTimeGPU Last spike time
		 * \param InternalSpikeGPU In this vector is stored if a neuron must generate an output spike.
		 * \param SizeStates Number of neurons
		 * \param CurrentTime Current time.
		 *
		 * \return True if an output spike have been fired. False in other case.
		 */
		__device__ void UpdateState(double CurrentTime)
		{
			int index = blockIdx.x * blockDim.x + threadIdx.x;
			while (index<vectorNeuronState_GPU2->SizeStates){
				
				vectorNeuronState_GPU2->VectorNeuronStates_GPU[(this->N_DifferentialNeuronState)*vectorNeuronState_GPU2->SizeStates + index]+=vectorNeuronState_GPU2->AuxStateGPU[index];
				vectorNeuronState_GPU2->VectorNeuronStates_GPU[(this->N_DifferentialNeuronState+1)*vectorNeuronState_GPU2->SizeStates + index]+=vectorNeuronState_GPU2->AuxStateGPU[vectorNeuronState_GPU2->SizeStates + index];

				vectorNeuronState_GPU2->InternalSpikeGPU[index]=this->integrationMethod_GPU2->NextDifferentialEcuationValues(index, vectorNeuronState_GPU2->SizeStates, vectorNeuronState_GPU2->VectorNeuronStates_GPU);

				vectorNeuronState_GPU2->LastUpdateGPU[index]=CurrentTime;

				this->CheckValidIntegration(index);

				index+=blockDim.x*gridDim.x;
			}
		} 


		/*!
		 * \brief It evaluates if a neuron must spike.
		 *
		 * It evaluates if a neuron must spike.
		 *
		 * \param previous_V previous membrane potential
		 * \param NeuronState neuron state variables.
		 * \param index Neuron index inside the neuron model.
 		 * \param elapsedTimeInNeuronModelScale integration method step.
		 * \return It returns if a neuron must spike.
		 */
		__device__ virtual bool EvaluateSpikeCondition(float previous_V, float * NeuronState, int index, float elapsedTimeInNeuronModelScale){
			if (NeuronState[index] > this->vthr && previous_V < this->vthr){
				vectorNeuronState_GPU2->LastSpikeTimeGPU[index]=0.0;
				return true;
			}
			return false;
		}


		__device__ void EvaluateDifferentialEcuation(int index, int SizeStates, float * NeuronState, float * AuxNeuronState, float elapsed_time){
			float V=NeuronState[index];
			float m=NeuronState[SizeStates + index];
			float h=NeuronState[2*SizeStates + index];
			float n=NeuronState[3*SizeStates + index];
			float gexc=NeuronState[4*SizeStates + index];
			float ginh=NeuronState[5*SizeStates + index];

			if(V<EK){
				V=EK;
			}
			if(V>ENa){
				V=ENa;
			}

			if(m<0){
				m=0.0f;
			}
			if(m>1){
				m=1.0f;
			}
			if(h<0){
				h=0.0f;
			}
			if(h>1){
				h=1.0f;
			}
			if(n<0){
				n=0.0f;
			}
			if(n>0.72){
				n=0.72f;
			}

			//v			
			AuxNeuronState[blockDim.x*blockIdx.x + threadIdx.x]=(gexc * (this->eexc - V) + ginh * (this->einh - V) + grest * (this->erest - V) + gNa*m*m*m*h*(ENa - V) + gKd*n*n*n*n*(EK - V))*this->inv_cm;
			
			//m
			float alpha_m, beta_m;
			if((13.0f-V+VT)!=0.0f){
				alpha_m=0.32f*(13.0f-V+VT)/(__expf((13.0f-V+VT)*0.25f)-1.0f);
			}else{
				alpha_m=1.279999894092334f;
			}

			if((V-VT-40.0f)!=0.0f){
				beta_m=0.28f*(V-VT-40.0f)/(__expf((V-VT-40.0f)*0.2f)-1.0f);
			}else{
				beta_m=1.400030971093304f;
			}
			AuxNeuronState[blockDim.x*gridDim.x + blockDim.x*blockIdx.x + threadIdx.x]=(alpha_m*(1.0f-m)-beta_m*m);
			
			//h
			float alpha_h=0.128f*__expf((17.0f-V+VT)*0.0555555555555555555f);
			float beta_h=4.0f/(1.0f+__expf((40.0f-V+VT)*0.2f));
			AuxNeuronState[2*blockDim.x*gridDim.x + blockDim.x*blockIdx.x + threadIdx.x]=(alpha_h*(1-h)-beta_h*h);
			
			//n
			float alpha_n;
			if((15.0f-V+VT)!=0.0f){
				alpha_n=0.032f*(15.0f-V+VT)/(__expf((15.0f-V+VT)*0.2f)-1.0f);
			}else{
				alpha_n=0.1599999867615418f;
			}
			float beta_n=0.5f*__expf((10.0f-V+VT)*0.025f);
			AuxNeuronState[3*blockDim.x*gridDim.x + blockDim.x*blockIdx.x + threadIdx.x]=(alpha_n*(1-n)-beta_n*n);
		}


		/*!
		 * \brief It evaluates the time depedendent ecuation in NeuronState for elapsed_time and it stores the results in NeuronState.
		 *
		 * It evaluates the time depedendent ecuation in NeuronState for elapsed_time and it stores the results in NeuronState.
		 *
		 * \param index index inside the NeuronState vector.
		 * \param SizeStates number of element in NeuronState vector.
		 * \param NeuronState value of the neuron state variables where time dependent equations are evaluated.
		 * \param elapsed_time integration time step.
		 * \param elapsed_time_index index inside the conductance_exp_values array.
		 */
		__device__ void EvaluateTimeDependentEcuation(int index, int SizeStates, float * NeuronState, float elapsed_time, int elapsed_time_index){
			float limit=1e-20;

			float * Conductance_values=this->Get_conductance_exponential_values(elapsed_time_index);
			
			if(NeuronState[this->N_DifferentialNeuronState*SizeStates + index]<limit){
				NeuronState[this->N_DifferentialNeuronState*SizeStates + index]=0.0f;
			}else{
				NeuronState[this->N_DifferentialNeuronState*SizeStates + index]*=  Conductance_values[0];
			}
			if(NeuronState[(this->N_DifferentialNeuronState+1)*SizeStates + index]<limit){
				NeuronState[(this->N_DifferentialNeuronState+1)*SizeStates + index]=0.0f;
			}else{
				NeuronState[(this->N_DifferentialNeuronState+1)*SizeStates + index]*=  Conductance_values[1];
			}
		}

		/*!
		 * \brief It calculates the conductace exponential value for an elapsed time.
		 *
		 * It calculates the conductace exponential value for an elapsed time.
		 *
		 * \param index elapsed time index .
		 * \param elapses_time elapsed time.
		 */
		__device__ void Calculate_conductance_exp_values(int index, float elapsed_time){
			//excitatory synapse.
			Set_conductance_exp_values(index, 0, __expf(-elapsed_time*this->inv_texc));
			//inhibitory synapse.
			Set_conductance_exp_values(index, 1, __expf(-elapsed_time*this->inv_tinh));
		}

};


#endif /* LIFTIMEDRIVENMODEL_1_2_GPU2_H_ */
