/***************************************************************************
 *                           LIFTimeDrivenModel_1_4.cpp                    *
 *                           -------------------                           *
 * copyright            : (C) 2013 by Jesus Garrido and Francisco Naveros  *
 * email                : jgarrido@atc.ugr.es, 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/LIFTimeDrivenModel_1_4.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"

const float LIFTimeDrivenModel_1_4::Max_V=0.0f;
const float LIFTimeDrivenModel_1_4::Min_V=-100.0f;

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

float * LIFTimeDrivenModel_1_4::channel_values=Generate_channel_values();


void LIFTimeDrivenModel_1_4::LoadNeuronModel(string ConfigFile) throw (EDLUTFileException){
	FILE *fh;
	long Currentline = 0L;
	fh=fopen(ConfigFile.c_str(),"rt");
	if(fh){
		Currentline=1L;
		skip_comments(fh,Currentline);
		if(fscanf(fh,"%f",&this->eexc)==1){
			skip_comments(fh,Currentline);

			if (fscanf(fh,"%f",&this->einh)==1){
				skip_comments(fh,Currentline);

				if(fscanf(fh,"%f",&this->erest)==1){
					skip_comments(fh,Currentline);

					if(fscanf(fh,"%f",&this->vthr)==1){
						skip_comments(fh,Currentline);

						if(fscanf(fh,"%f",&this->cm)==1){
							inv_cm = 1.0f / cm;
							skip_comments(fh,Currentline);

							if(fscanf(fh,"%f",&this->tampa)==1){
								inv_tampa=1.0f/tampa;
								skip_comments(fh,Currentline);

								if(fscanf(fh,"%f",&this->tnmda)==1){
									inv_tnmda=1.0f/tnmda;
									skip_comments(fh,Currentline);
									
									if(fscanf(fh,"%f",&this->tinh)==1){
										inv_tinh=1.0f/tinh;
										skip_comments(fh,Currentline);

										if(fscanf(fh,"%f",&this->tgj)==1){
											inv_tgj=1.0f/tgj;
											skip_comments(fh,Currentline);

											if(fscanf(fh,"%f",&this->tref)==1){
												skip_comments(fh,Currentline);

												if(fscanf(fh,"%f",&this->grest)==1){
													skip_comments(fh,Currentline);

													if(fscanf(fh,"%f",&this->fgj)==1){
														skip_comments(fh,Currentline);

														this->State = (VectorNeuronState *) new VectorNeuronState(N_NeuronStateVariables, true);

													} else {
														throw EDLUTFileException(13,60,3,1,Currentline);
													}
												} else {
													throw EDLUTFileException(13,61,3,1,Currentline);
												}
											} else {
												throw EDLUTFileException(13,62,3,1,Currentline);
											}
										} else {
											throw EDLUTFileException(13,63,3,1,Currentline);
										}
									} else {
										throw EDLUTFileException(13,64,3,1,Currentline);
									}
								} else {
									throw EDLUTFileException(13,65,3,1,Currentline);
								}
							} else {
								throw EDLUTFileException(13,66,3,1,Currentline);
							}
						} else {
							throw EDLUTFileException(13,67,3,1,Currentline);
						}
					} else {
						throw EDLUTFileException(13,68,3,1,Currentline);
					}
				} else {
					throw EDLUTFileException(13,69,3,1,Currentline);
				}
			} else {
				throw EDLUTFileException(13,70,3,1,Currentline);
			}
		} else {
			throw EDLUTFileException(13,71,3,1,Currentline);
		}

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

		this->integrationMethod->SetBiFixedStepParameters((erest+vthr)/2,(2*erest+vthr)/3,0);

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

void LIFTimeDrivenModel_1_4::SynapsisEffect(int index, Interconnection * InputConnection){
	this->GetVectorNeuronState()->IncrementStateVariableAtCPU(index,N_DifferentialNeuronState+InputConnection->GetType(),InputConnection->GetWeight());
}

//this neuron model is implemented in a second scale.
LIFTimeDrivenModel_1_4::LIFTimeDrivenModel_1_4(string NeuronTypeID, string NeuronModelID): TimeDrivenNeuronModel(NeuronTypeID, NeuronModelID, SecondScale), eexc(0), einh(0), erest(0), vthr(0), cm(0), tampa(0), tnmda(0), tinh(0), tgj(0),
		tref(0), grest(0){
}

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

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

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


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

	return 0;
}


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

	
	bool * internalSpike=State->getInternalSpike();
	int Size=State->GetSizeState();

	//NeuronState[0] --> vm 
	//NeuronState[1] --> gampa 
	//NeuronState[2] --> gnmda 
	//NeuronState[3] --> ginh 
	//NeuronState[4] --> ggj 


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

		float * NeuronState=State->GetStateVariableAt(i);

		bool spike=false;
		internalSpike[i]=this->integrationMethod->NextDifferentialEcuationValues(i, NeuronState);

		State->SetLastUpdateTime(i,CurrentTime);
	}

	CheckValidIntegration();

	return false;
}



ostream & LIFTimeDrivenModel_1_4::PrintInfo(ostream & out){
	out << "- Leaky Time-Driven Model 1_4: " << 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\tAMPA Time Constant: " << this->tampa << "sNMDA Time Constant: " << this->tnmda << "s" << endl;

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

	return out;
}	



void LIFTimeDrivenModel_1_4::InitializeStates(int N_neurons, int OpenMPQueueIndex){
	//Initialize neural state variables.
	float initialization[] = {erest,0.0f,0.0f,0.0f,0.0f};
	State->InitializeStates(N_neurons, initialization);

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


bool LIFTimeDrivenModel_1_4::EvaluateSpikeCondition(float previous_V, float * NeuronState, int index, float elapsedTimeInNeuronModelScale){
	float vm_cou = NeuronState[0] + this->fgj * NeuronState[4];
	if (vm_cou > this->vthr  && this->GetVectorNeuronState()->GetLastSpikeTime(index)>this->tref){
		NeuronState[0] = this->erest;
		State->NewFiredSpike(index);
		this->integrationMethod->resetState(index);
		return true;
	}
	return false;
}


//void LIFTimeDrivenModel_1_4::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState, int index, float elapsed_time){
//	if(this->GetVectorNeuronState()->GetLastSpikeTime(index)>this->tref){
//		float iampa = NeuronState[1]*(this->eexc-NeuronState[0]);
//		//float gnmdainf = 1.0f/(1.0f + ExponentialTable::GetResult(-62.0f*NeuronState[0])*(1.2f/3.57f));
//		float gnmdainf = *Get_channel_values(NeuronState[0]);
//		float inmda = NeuronState[2]*gnmdainf*(this->eexc-NeuronState[0]);
//		float iinh = NeuronState[3]*(this->einh-NeuronState[0]);
//		AuxNeuronState[0]=(iampa + inmda + iinh + this->grest* (this->erest-NeuronState[0]))*1.e-9f/this->cm;
//	}else if((this->GetVectorNeuronState()->GetLastSpikeTime(index)+elapsed_time)>this->tref){
//		float fraction=(this->GetVectorNeuronState()->GetLastSpikeTime(index)+elapsed_time-this->tref)/elapsed_time;
//		float iampa = NeuronState[1]*(this->eexc-NeuronState[0]);
//		//float gnmdainf = 1.0f/(1.0f + ExponentialTable::GetResult(-62.0f*NeuronState[0])*(1.2f/3.57f));
//		float gnmdainf = *Get_channel_values(NeuronState[0]);
//		float inmda = NeuronState[2]*gnmdainf*(this->eexc-NeuronState[0]);
//		float iinh = NeuronState[3]*(this->einh-NeuronState[0]);
//		AuxNeuronState[0]=fraction*(iampa + inmda + iinh + this->grest* (this->erest-NeuronState[0]))*1.e-9f/this->cm;
//	}else{
//		AuxNeuronState[0]=0;
//	}
//}

void LIFTimeDrivenModel_1_4::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState, int index, float elapsed_time){
	if (this->GetVectorNeuronState()->GetLastSpikeTime(index)>this->tref){
		float iampa = NeuronState[1] * (this->eexc - NeuronState[0]);
		//float gnmdainf = 1.0f/(1.0f + ExponentialTable::GetResult(-62.0f*NeuronState[0])*(1.2f/3.57f));
		float gnmdainf = *Get_channel_values(NeuronState[0]);
		float inmda = NeuronState[2] * gnmdainf*(this->eexc - NeuronState[0]);
		float iinh = NeuronState[3] * (this->einh - NeuronState[0]);
		AuxNeuronState[0] = (iampa + inmda + iinh + this->grest* (this->erest - NeuronState[0]))*1.e-9f * this->inv_cm;
	}
	else if ((this->GetVectorNeuronState()->GetLastSpikeTime(index) + elapsed_time)>this->tref){
		float fraction = (this->GetVectorNeuronState()->GetLastSpikeTime(index) + elapsed_time - this->tref) / elapsed_time;
		float iampa = NeuronState[1] * (this->eexc - NeuronState[0]);
		//float gnmdainf = 1.0f/(1.0f + ExponentialTable::GetResult(-62.0f*NeuronState[0])*(1.2f/3.57f));
		float gnmdainf = *Get_channel_values(NeuronState[0]);
		float inmda = NeuronState[2] * gnmdainf*(this->eexc - NeuronState[0]);
		float iinh = NeuronState[3] * (this->einh - NeuronState[0]);
		AuxNeuronState[0] = fraction*(iampa + inmda + iinh + this->grest* (this->erest - NeuronState[0]))*1.e-9f * this->inv_cm;
	}
	else{
		AuxNeuronState[0] = 0;
	}
}

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

}

void LIFTimeDrivenModel_1_4::Calculate_conductance_exp_values(int index, float elapsed_time){
	//ampa synapse.
	Set_conductance_exp_values(index, 0, exp(-elapsed_time*this->inv_tampa));
	//nmda synapse.
	Set_conductance_exp_values(index, 1, exp(-elapsed_time*this->inv_tnmda));
	//inhibitory synapse.
	Set_conductance_exp_values(index, 2, exp(-elapsed_time*this->inv_tinh));
	//gap junction synapse.
	Set_conductance_exp_values(index, 3, exp(-elapsed_time*this->inv_tgj));
}


int LIFTimeDrivenModel_1_4::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;
	}
}

