/***************************************************************************
 *                           HHTimeDrivenModel.cpp		                   *
 *                           -------------------                           *
 * copyright            : (C) 2015 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.                                   *
 *                                                                         *
 ***************************************************************************/

#include "../../include/neuron_model/HHTimeDrivenModel.h"
#include "../../include/neuron_model/VectorNeuronState.h"

#include <iostream>
#include <cmath>
#include <string>

#include "../../include/openmp/openmp.h"

#include "../../include/spike/EDLUTFileException.h"
#include "../../include/spike/Neuron.h"
#include "../../include/spike/InternalSpike.h"
#include "../../include/spike/PropagatedSpike.h"
#include "../../include/spike/Interconnection.h"

#include "../../include/simulation/Utils.h"

#include "../../include/openmp/openmp.h"


const float HHTimeDrivenModel::Max_V=50.0f;
const float HHTimeDrivenModel::Min_V=-90.0f;

const float HHTimeDrivenModel::aux=(HHTimeDrivenModel::TableSize-1)/( HHTimeDrivenModel::Max_V- HHTimeDrivenModel::Min_V);

float * HHTimeDrivenModel::channel_values=Generate_channel_values();


void HHTimeDrivenModel::LoadNeuronModel(string ConfigFile) throw (EDLUTFileException){
	FILE *fh;
	long Currentline = 0L;
	fh=fopen(ConfigFile.c_str(),"rt");
	if(fh){
		this->State = (VectorNeuronState *) new VectorNeuronState(N_NeuronStateVariables, true);

		//INTEGRATION METHOD
		this->integrationMethod = LoadIntegrationMethod::loadIntegrationMethod((TimeDrivenNeuronModel *)this, fh, &Currentline, N_NeuronStateVariables, N_DifferentialNeuronState, N_TimeDependentNeuronState);

		this->integrationMethod->SetBiFixedStepParameters(VT-5,VT-5,1.0);

		this->integrationMethod->Calculate_conductance_exp_values();
	}
	fclose(fh);
}

//this neuron model is implemented in a second scale.
void HHTimeDrivenModel::SynapsisEffect(int index, Interconnection * InputConnection){
	this->GetVectorNeuronState()->IncrementStateVariableAtCPU(index,N_DifferentialNeuronState+InputConnection->GetType(),InputConnection->GetWeight());
}

HHTimeDrivenModel::HHTimeDrivenModel(string NeuronTypeID, string NeuronModelID): TimeDrivenNeuronModel(NeuronTypeID, NeuronModelID, MilisecondScale){
	eexc=0.0f; //mV
	einh=-80.0f; //mV
	erest=-65.0f; //mV
	vthr=-30.0f; //mV
	cm=120.0f; //pF   
	inv_cm=1.0f/cm;
	texc=5.0f;//ms
	inv_texc=1.0f/texc;
	tinh=10.0f;//ms
	inv_tinh=1.0f/tinh;
	grest=10.0f;//nS
	gNa=20000.0;//nS
	gKd=6000.0f;//nS
	ENa=50.0f;//mV
	EK=-90.0f;//mV
	VT=-52.0f;//mV
}

HHTimeDrivenModel::~HHTimeDrivenModel(void)
{
	//if(this->channel_values){
	//	delete this->channel_values;
	//	this->channel_values=0;
	//}
}

void HHTimeDrivenModel::LoadNeuronModel() throw (EDLUTFileException){
	this->LoadNeuronModel(this->GetModelID()+".cfg");
}

VectorNeuronState * HHTimeDrivenModel::InitializeState(){
	return this->GetVectorNeuronState();
}


InternalSpike * HHTimeDrivenModel::ProcessInputSpike(Interconnection * inter, Neuron * target, double time){
	// Add the effect of the input spike
	this->SynapsisEffect(target->GetIndex_VectorNeuronState(),inter);

	return 0;
}






bool HHTimeDrivenModel::UpdateState(int index, double CurrentTime){

	
	bool * internalSpike=State->getInternalSpike();


	//float * NeuronState;
	//NeuronState[0] --> vm 
	//NeuronState[1] --> gexc 
	//NeuronState[2] --> ginh 

	for (int i = 0; i< State->GetSizeState(); i++){

		float * NeuronState=State->GetStateVariableAt(i);
	
		internalSpike[i]=this->integrationMethod->NextDifferentialEcuationValues(i, NeuronState);

		State->SetLastUpdateTime(i,CurrentTime);
	}

	CheckValidIntegration();

	return false;
}



ostream & HHTimeDrivenModel::PrintInfo(ostream & out){
	//out << "- Leaky Time-Driven Model: " << this->GetModelID() << endl;

	//out << "\tExc. Reversal Potential: " << this->eexc << "V\tInh. Reversal Potential: " << this->einh << "V\tResting potential: " << this->erest << "V" << endl;

	//out << "\tFiring threshold: " << this->vthr << "V\tMembrane capacitance: " << this->cm << "nS\tExcitatory Time Constant: " << this->texc << "s" << endl;

	//out << "\tInhibitory time constant: " << this->tinh << "s\tRefractory Period: " << this->tref << "s\tResting Conductance: " << this->grest << "nS" << endl;

	return out;
}	



void HHTimeDrivenModel::InitializeStates(int N_neurons, int OpenMPQueueIndex){
	//Initialize neural state variables.
	float * values=Get_channel_values(erest);

	//m
	float alpha_m=values[0];
	float beta_m=values[1];
	float m_inf=alpha_m/(alpha_m+beta_m);

	//h
	float alpha_h=values[2];
	float beta_h=values[3];
	float h_inf=alpha_h/(alpha_h+beta_h);


	//n
	float alpha_n=values[4];
	float beta_n=values[5];
	float n_inf=alpha_n/(alpha_n+beta_n);

	float initialization[] = {erest,m_inf, h_inf, n_inf, 0.0f, 0.0f};
	State->InitializeStates(N_neurons, initialization);

	//Initialize integration method state variables.
	this->integrationMethod->InitializeStates(N_neurons, initialization);
}


bool HHTimeDrivenModel::EvaluateSpikeCondition(float previous_V, float * NeuronState, int index, float elpasedTimeInNeuronModelScale){
	if (NeuronState[0] > this->vthr && previous_V < this->vthr){
		State->NewFiredSpike(index);
		return true;
	}
	return false;
}



void HHTimeDrivenModel::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState, int index, float elapsed_time){
	float V=NeuronState[0];
	float m=NeuronState[1];
	float h=NeuronState[2];
	float n=NeuronState[3];
	float gexc=NeuronState[4];
	float ginh=NeuronState[5];

	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[0]=(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;


	float * values=Get_channel_values(V);

	//m
	float alpha_m=values[0];
	float beta_m=values[1];
	AuxNeuronState[1]=(alpha_m*(1.0f-m)-beta_m*m);



	//h
	float alpha_h=values[2];
	float beta_h=values[3];
	AuxNeuronState[2]=(alpha_h*(1-h)-beta_h*h);

	
	//n
	float alpha_n=values[4];
	float beta_n=values[5];
	AuxNeuronState[3]=(alpha_n*(1-n)-beta_n*n);
}



void HHTimeDrivenModel::EvaluateTimeDependentEcuation(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[N_DifferentialNeuronState]<limit){
		NeuronState[N_DifferentialNeuronState]=0.0f;
	}else{
		NeuronState[N_DifferentialNeuronState]*= Conductance_values[0];
	}
	if(NeuronState[N_DifferentialNeuronState+1]<limit){
		NeuronState[N_DifferentialNeuronState+1]=0.0f;
	}else{
		NeuronState[N_DifferentialNeuronState+1]*= Conductance_values[1];
	}	
}


void HHTimeDrivenModel::Calculate_conductance_exp_values(int index, float elapsed_time){
	//excitatory synapse.
	Set_conductance_exp_values(index, 0, exp(-elapsed_time*this->inv_texc));
	//inhibitory synapse.
	Set_conductance_exp_values(index, 1, exp(-elapsed_time*this->inv_tinh));
}


int HHTimeDrivenModel::CheckSynapseTypeNumber(int Type){
	if(Type<N_TimeDependentNeuronState && Type>=0){
		return Type;
	}else{
		cout<<"Neuron model "<<this->GetTypeID()<<", "<<this->GetModelID()<<" does not support input synapses of type "<<Type<<endl;
		return 0;
	}
}




