/***************************************************************************
 *                           Bifixed_RK2.cpp                                       *
 *                           -------------------                           *
 * 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.                                   *
 *                                                                         *
 ***************************************************************************/

#include "../../include/integration_method/Bifixed_RK2.h"
#include "../../include/neuron_model/TimeDrivenNeuronModel.h"
#include "../../include/neuron_model/VectorNeuronState.h"


Bifixed_RK2::Bifixed_RK2(TimeDrivenNeuronModel * NewModel, int N_neuronStateVariables, int N_differentialNeuronState, int N_timeDependentNeuronState):BiFixedStep(NewModel,"Bifixed_RK2", N_neuronStateVariables, N_differentialNeuronState, N_timeDependentNeuronState, false, false){

}

Bifixed_RK2::~Bifixed_RK2(){
}
		

bool Bifixed_RK2::NextDifferentialEcuationValues(int index, float * NeuronState){

	float AuxNeuronState[MAX_VARIABLES];
	float AuxNeuronState1[MAX_VARIABLES];
	float AuxNeuronState2[MAX_VARIABLES];
	float PreviousNeuronState[MAX_VARIABLES];
	bool spike=false;


	if(integrationMethodState[index]==0){
		memcpy(PreviousNeuronState,NeuronState,sizeof(float)*N_NeuronStateVariables);

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

		this->model->EvaluateTimeDependentEcuation(AuxNeuronState, elapsedTimeInNeuronModelScale, 0);
		this->model->EvaluateDifferentialEcuation(AuxNeuronState, AuxNeuronState2, index, elapsedTimeInNeuronModelScale);


		for (int j=0; j<N_DifferentialNeuronState; j++){
			NeuronState[j]+=(AuxNeuronState1[j]+AuxNeuronState2[j])*elapsedTimeInNeuronModelScale*0.5f;
		}

		this->model->EvaluateTimeDependentEcuation(NeuronState, elapsedTimeInNeuronModelScale, 0);

		if(NeuronState[0]>startVoltageThreshold){
			integrationMethodState[index]=1;
			//Restore the neuron model state to a previous state.
			memcpy(NeuronState,PreviousNeuronState,sizeof(float)*N_NeuronStateVariables);
		}else{
			//Update the last spike time.
			this->model->GetVectorNeuronState()->AddElapsedTime(index,elapsedTimeInSeconds);
		}
	}

	if(integrationMethodState[index]>0){
		for(int iteration=0; iteration<N_BiFixedSteps; iteration++){
			float previous_V=NeuronState[0];

			//1st term
			this->model->EvaluateDifferentialEcuation(NeuronState, AuxNeuronState1, index, bifixedElapsedTimeInNeuronModelScale);
			
			//2nd term
			memcpy(AuxNeuronState, NeuronState,sizeof(float)*N_NeuronStateVariables);
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState[j]= NeuronState[j] + AuxNeuronState1[j]*bifixedElapsedTimeInNeuronModelScale;
			}


			this->model->EvaluateTimeDependentEcuation(AuxNeuronState, bifixedElapsedTimeInNeuronModelScale, 1);
			this->model->EvaluateDifferentialEcuation(AuxNeuronState, AuxNeuronState2, index, bifixedElapsedTimeInNeuronModelScale);


			for (int j=0; j<N_DifferentialNeuronState; j++){
				NeuronState[j]+=(AuxNeuronState1[j]+AuxNeuronState2[j])*bifixedElapsedTimeInNeuronModelScale*0.5f;
			}

			this->model->EvaluateTimeDependentEcuation(NeuronState, bifixedElapsedTimeInNeuronModelScale, 1);

			//Update the last spike time.
			this->model->GetVectorNeuronState()->AddElapsedTime(index,bifixedElapsedTimeInSeconds);

			bool aux_spike=this->model->EvaluateSpikeCondition(previous_V, NeuronState, index, bifixedElapsedTimeInNeuronModelScale);
			if(aux_spike){
				spike=true;
			}

			if(NeuronState[0]>startVoltageThreshold && integrationMethodState[index]==1){
				integrationMethodState[index]=2;
			}else if(NeuronState[0]<endVoltageThreshold && integrationMethodState[index]==2){
				integrationMethodState[index]=3;
				integrationMethodCounter[index]=N_postBiFixedSteps;
			}
			if(integrationMethodCounter[index]>0 && integrationMethodState[index]==3){
				integrationMethodCounter[index]--;
				if(integrationMethodCounter[index]==0){
					integrationMethodState[index]=0;
				}
			}
		}
		if(integrationMethodState[index]==1){
			integrationMethodState[index]=0;
		}
	}

	return spike;
}


ostream & Bifixed_RK2::PrintInfo(ostream & out){
	out << "Integration Method Type: " << this->GetType() << endl;

	return out;
}	

void Bifixed_RK2::InitializeStates(int N_neurons, float * initialization){
	integrationMethodCounter=new int [N_neurons]();
	integrationMethodState=new int [N_neurons]();
}


void Bifixed_RK2::Calculate_conductance_exp_values(){
	this->model->Initialize_conductance_exp_values(this->N_TimeDependentNeuronState,2);
	//index 0
	this->model->Calculate_conductance_exp_values(0, elapsedTimeInNeuronModelScale);
	//index 1
	this->model->Calculate_conductance_exp_values(1, bifixedElapsedTimeInNeuronModelScale);
}