/***************************************************************************
 *                           TimeDrivenPurkinjeCell_IP.cpp                    *
 *                           -------------------                           *
 * copyright            : (C) 2016 by Richard Carrillo, Niceto Luque and   *
						  Francisco Naveros								   *
 * email                : rcarrillo@ugr.es, nluque@ugr.es and			   *
						  fnaveros@atc.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/TimeDrivenPurkinjeCell_IP.h"
#include "../../include/neuron_model/VectorNeuronState.h"

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

#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 TimeDrivenPurkinjeCell_IP::Max_V=35.0f;
const float TimeDrivenPurkinjeCell_IP::Min_V=-100.0f;

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

float * TimeDrivenPurkinjeCell_IP::channel_values=Generate_channel_values();


void TimeDrivenPurkinjeCell_IP::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->external_I) == 1){
			this->State = (VectorNeuronState *) new VectorNeuronState(N_NeuronStateVariables, true);
		}
		else{
			printf("ERROR LOADING EXTERNAL CURRIENT IN TIMEDRIVENPURKINJECELL_IP MODEL");
		}

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

		//CHECK THIS VALUES.
		this->integrationMethod->SetBiFixedStepParameters((erest+3*vthr)/4,(erest+3*vthr)/4,0.2);

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

void TimeDrivenPurkinjeCell_IP::SynapsisEffect(int index, Interconnection * InputConnection){
	this->GetVectorNeuronState()->IncrementStateVariableAtCPU(index,N_DifferentialNeuronState+InputConnection->GetType(),1e-6f*InputConnection->GetWeight());
}

TimeDrivenPurkinjeCell_IP::TimeDrivenPurkinjeCell_IP(string NeuronTypeID, string NeuronModelID): TimeDrivenNeuronModel(NeuronTypeID, NeuronModelID, MilisecondScale), 
		tau_IP(12000.0f),inv_tau_IP(1/12000.0f),g_L(0.02f),g_Ca(0.001f), g_M(0.75f), Cylinder_length_of_the_soma(0.0015f), Radius_of_the_soma(0.0008f), 
		Area(3.141592f*0.0015f*2.0f*0.0008f), inv_Area(1.0f/(3.141592f*0.0015f*2.0f*0.0008f)), Membrane_capacitance_const(1.0f), 
		inv_Membrane_capacitance_const(1.0f/1.0f)
		{
	eexc=0.0f;
	einh=-80.0f ;
	vthr= -35.0f;
	erest=-70.0f;
	texc=1.0f;
	inv_texc=1.0f/texc;
	tinh=2;
	inv_tinh=1.0f/tinh;
	tref=1.35f;
	tref_0_5=tref*0.5f;
	inv_tref_0_5=1.0f/tref_0_5;
	spkpeak=31.0f;
	ep_Capacitance=0.0475; 
	external_I=0.0f;
}


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

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

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


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

	return 0;
}




bool TimeDrivenPurkinjeCell_IP::UpdateState(int index, double CurrentTime){
	//float * NeuronState;
	//NeuronState[0] --> V
	//NeuronState[1] --> c
	//NeuronState[2] --> M
	//NeuronState[3] --> Capacitance
	//NeuronState[4] --> gexc 
	//NeuronState[5] --> ginh 


	bool * internalSpike=State->getInternalSpike();

	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 & TimeDrivenPurkinjeCell_IP::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 TimeDrivenPurkinjeCell_IP::InitializeStates(int N_neurons, int OpenMPQueueIndex){
	//Initialize neural state variables.
	float * values=Get_channel_values(erest);
	float alpha_ca=values[0];
	float inv_tau_ca=values[1];
	float alpha_M=values[2];
	float inv_tau_M=values[3];

	//c_inf
	float c_inf=alpha_ca/inv_tau_ca;

	//M_inf
	float M_inf=alpha_M/inv_tau_M;

	float initialization[] = {erest,c_inf,M_inf,Membrane_capacitance_const,0.0f,0.0f};
	State->InitializeStates(N_neurons, initialization);

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


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

	double last_spike = State->GetLastSpikeTime(index) * this->timeScale;

	if(last_spike < tref){
		if(last_spike <= tref_0_5){
			NeuronState[0]=vthr+(spkpeak-vthr)*(last_spike*inv_tref_0_5);
		}else{
			NeuronState[0]=spkpeak-(spkpeak-erest)*((last_spike-tref_0_5)*inv_tref_0_5);
		}
	}else if((last_spike - tref)<elapsedTimeInNeuronModelScale){
		NeuronState[0]=erest;
	}
	
	if(spike){
		//Capacitance
		float Inject_current=(NeuronState[4] * (this->eexc - NeuronState[0]) + NeuronState[5] * (this->einh - NeuronState[0]))*inv_Area;
		NeuronState[3]+=ep_Capacitance*inv_tau_IP;  	
	}



	return spike;
}


void TimeDrivenPurkinjeCell_IP::EvaluateDifferentialEcuation(float * NeuronState, float * AuxNeuronState, int index, float elapsed_time){
	float V=NeuronState[0];
	float ca=NeuronState[1];
	float M=NeuronState[2];
	float Membrane_capacitance=NeuronState[3];
	float g_exc=NeuronState[4];
	float g_inh=NeuronState[5];
	float last_spike=this->timeScale*State->GetLastSpikeTime(index);


	float Inject_current=(g_exc * (this->eexc - V) + g_inh * (this->einh - V) + this->external_I)*inv_Area;
	//V
	if(last_spike >= tref){
		AuxNeuronState[0]=(-g_L*(V+70.0f)-g_Ca*ca*ca*(V-125.0f)-g_M*M*(V+95.0f) + Inject_current )/Membrane_capacitance;
	
	}else if(last_spike <= tref_0_5){
		AuxNeuronState[0]=(spkpeak-vthr)*inv_tref_0_5;
	}else{
		AuxNeuronState[0]=(erest-spkpeak)*inv_tref_0_5;
	}

	float * values=Get_channel_values(V);

	//ca
	float alpha_ca=values[0];
	float inv_tau_ca=values[1];
	AuxNeuronState[1]=alpha_ca - ca*inv_tau_ca;
	
	//M	
	float alpha_M=values[2];
	float inv_tau_M=values[3];
	AuxNeuronState[2]=alpha_M - M*inv_tau_M;

	//Membrane Capacitance
	if (V<this->vthr){
		AuxNeuronState[3]=-(Membrane_capacitance-Membrane_capacitance_const*0.5)/(Membrane_capacitance_const*tau_IP*(1+Inject_current));
	}else{
		AuxNeuronState[3]=0.0f;
	}

}

void TimeDrivenPurkinjeCell_IP::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 TimeDrivenPurkinjeCell_IP::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 TimeDrivenPurkinjeCell_IP::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;
	}
}







