/***************************************************************************
 *                           Bifixed_BDFn.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_BDFn.h"
#include "../../include/neuron_model/TimeDrivenNeuronModel.h"
#include "../../include/neuron_model/VectorNeuronState.h"

#include <math.h>

//const float Bifixed_BDFn::Coeficient [7][7]={{1.0f,1.0f,0.0f,0.0f,0.0f,0.0f,0.0f},
//{1.0f,1.0f,0.0f,0.0f,0.0f,0.0f,0.0f},
//{2.0f/3.0f,4.0f/3.0f,-1.0f/3.0f,0.0f,0.0f,0.0f,0.0f},
//{6.0f/11.0f,18.0f/11.0f,-9.0f/11.0f,2.0f/11.0f,0.0f,0.0f,0.0f},
//{12.0f/25.0f,48.0f/25.0f,-36.0f/25.0f,16.0f/25.0f,-3.0f/25.0f,0.0f,0.0f},
//{60.0f/137.0f,300.0f/137.0f,-300.0f/137.0f,200.0f/137.0f,-75.0f/137.0f,12.0f/137.0f,0.0f},
//{60.0f/147.0f,360.0f/147.0f,-450.0f/147.0f,400.0f/147.0f,-225.0f/147.0f,72.0f/147.0f,-10.0f/147.0f}};


Bifixed_BDFn::Bifixed_BDFn(TimeDrivenNeuronModel * NewModel, int N_neuronStateVariables, int N_differentialNeuronState, int N_timeDependentNeuronState, int BDForder):BiFixedStep(NewModel,"Bifixed_BDFn", N_neuronStateVariables, N_differentialNeuronState, N_timeDependentNeuronState, true, true),BDForder(BDForder){
	Coeficient[0][0][0]=0;
}

Bifixed_BDFn::~Bifixed_BDFn(){
	if(BDForder>1){
		for(int i=0; i<BDForder-1; i++){
			delete [] PreviousNeuronState[i];
		}
		delete [] PreviousNeuronState;
	}
	
	for(int i=0; i<BDForder; i++){
		delete [] D[i];
	}
	delete [] D;
	delete [] state;

	delete PreviousIntegrationStep;
	delete CoeficientSelector;
}
		
bool Bifixed_BDFn::NextDifferentialEcuationValues(int index, float * NeuronState){
	
	if(Coeficient[0][0][0]==0){
		float decrement=bifixedElapsedTimeInNeuronModelScale/elapsedTimeInNeuronModelScale;
		float increment=elapsedTimeInNeuronModelScale/bifixedElapsedTimeInNeuronModelScale;

		//decrement of h
		Coeficient[0][0][0]=1.0f;
		Coeficient[0][0][1]=0.0f;
		Coeficient[0][0][2]=1.0f;
		Coeficient[0][0][3]=0.0f;
		Coeficient[0][1][0]=1.0f;
		Coeficient[0][1][1]=0.0f;
		Coeficient[0][1][2]=1.0f;
		Coeficient[0][1][3]=0.0f;
		Coeficient[0][2][0]=2.0f/3.0f; //beta 0
		Coeficient[0][2][1]=(decrement-1)/3.0f; //beta 1
		Coeficient[0][2][2]=1+decrement*decrement/3.0f;//alpha 0
		Coeficient[0][2][3]=-decrement*decrement/3.0f;//alpha 1

		//Equal h;
		Coeficient[1][0][0]=1.0f;
		Coeficient[1][0][1]=0.0f;
		Coeficient[1][0][2]=1.0f;
		Coeficient[1][0][3]=0.0f;
		Coeficient[1][1][0]=1.0f;
		Coeficient[1][1][1]=0.0f;
		Coeficient[1][1][2]=1.0f;
		Coeficient[1][1][3]=0.0f;
		Coeficient[1][2][0]=2.0f/3.0f;
		Coeficient[1][2][1]=0.0f;
		Coeficient[1][2][2]=4.0f/3.0f;
		Coeficient[1][2][3]=-1.0f/3.0f;

		//increment of h
		Coeficient[2][0][0]=1.0f;
		Coeficient[2][0][1]=0.0f;
		Coeficient[2][0][2]=1.0f;
		Coeficient[2][0][3]=0.0f;
		Coeficient[2][1][0]=1.0f;
		Coeficient[2][1][1]=0.0f;
		Coeficient[2][1][2]=1.0f;
		Coeficient[2][1][3]=0.0f;
		Coeficient[2][2][0]=2.0f/3.0f;
		Coeficient[2][2][1]=(increment-1)/3.0f;
		Coeficient[2][2][2]=1+increment*increment/3.0f;
		Coeficient[2][2][3]=-increment*increment/3.0f;
	}

	bool spike=false;



	float AuxNeuronState[MAX_VARIABLES];
	float AuxNeuronState1[MAX_VARIABLES];
	float AuxNeuronState_p[MAX_VARIABLES];
	float AuxNeuronState_p1[MAX_VARIABLES];
	float AuxNeuronState_c[MAX_VARIABLES];
	float jacnum[MAX_VARIABLES*MAX_VARIABLES];
	float J[MAX_VARIABLES*MAX_VARIABLES];
	float inv_J[MAX_VARIABLES*MAX_VARIABLES];

	if(integrationMethodState[index]==0){
		//If the state of the cell is 0, we use a Euler method to calculate an aproximation of the solution.
		if(state[index]==0){
			this->model->EvaluateDifferentialEcuation(NeuronState, AuxNeuronState, index, elapsedTimeInNeuronModelScale);
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState_p[j]= NeuronState[j] + elapsedTimeInNeuronModelScale*AuxNeuronState[j];
			}
		}
		//In this case we use the value of previous states to calculate an aproximation of the solution.
		else{
			float inverseIntegrationStepFactor=PreviousIntegrationStep[index]/elapsedTimeInNeuronModelScale;
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState_p[j]= NeuronState[j];
				for (int i=0; i<state[index]; i++){
					AuxNeuronState_p[j]+=D[i][index*N_DifferentialNeuronState+j]/*/inverseIntegrationStepFactor*/;
				}
			}
		}

		for(int i=N_DifferentialNeuronState; i<N_NeuronStateVariables; i++){
			AuxNeuronState_p[i]=NeuronState[i];
		}


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



		float epsi=1.0f;
		int k=0;

		this->model->EvaluateDifferentialEcuation(NeuronState, AuxNeuronState1, index, elapsedTimeInNeuronModelScale);

		//This integration method is an implicit method. We use this loop to iteratively calculate the implicit value.
		//epsi is the difference between two consecutive aproximation of the implicit method. 
		while (epsi>1e-16 && k<5){
			this->model->EvaluateDifferentialEcuation(AuxNeuronState_p, AuxNeuronState, index, elapsedTimeInNeuronModelScale);
	
			for (int j=0; j<N_DifferentialNeuronState; j++){
				AuxNeuronState_c[j]=Coeficient[CoeficientSelector[index]][state[index]][0]*elapsedTimeInNeuronModelScale*AuxNeuronState[j] + Coeficient[CoeficientSelector[index]][state[index]][1]*elapsedTimeInNeuronModelScale*AuxNeuronState1[j] + Coeficient[CoeficientSelector[index]][state[index]][2]*NeuronState[j];
				for (int i=1; i<state[index]; i++){
					AuxNeuronState_c[j]+=Coeficient[CoeficientSelector[index]][state[index]][i+2]*PreviousNeuronState[i-1][index*N_DifferentialNeuronState + j];
				}
			}

			//jacobian.
			Jacobian(AuxNeuronState_p, jacnum, index);
		
			for(int z=0; z<N_DifferentialNeuronState; z++){
				for(int t=0; t<N_DifferentialNeuronState; t++){
					J[z*N_DifferentialNeuronState + t] = Coeficient[CoeficientSelector[index]][state[index]][0] * elapsedTimeInNeuronModelScale * jacnum[z*N_DifferentialNeuronState + t];
					if(z==t){
						J[z*N_DifferentialNeuronState + t]-=1;
					}
				}
			}

			this->invermat(J,inv_J);

			for(int z=0; z<N_DifferentialNeuronState; z++){
				float aux=0.0;
				for (int t=0; t<N_DifferentialNeuronState; t++){
					aux+=inv_J[z*N_DifferentialNeuronState+t]*(AuxNeuronState_p[t]-AuxNeuronState_c[t]);
				}
				AuxNeuronState_p1[z]=aux + AuxNeuronState_p[z];
			}

			//We calculate the difference between both aproximations.
			float aux=0.0f;
			float aux2=0.0f;
			for(int z=0; z<N_DifferentialNeuronState; z++){
				aux=fabs(AuxNeuronState_p1[z]-AuxNeuronState_p[z]);
				if(aux>aux2){
					aux2=aux;
				}
			}

			memcpy(AuxNeuronState_p , AuxNeuronState_p1 ,sizeof(float)* N_DifferentialNeuronState);

			epsi=aux2;
			k++;
		}

		if(NeuronState[0]>startVoltageThreshold){
			integrationMethodState[index]=1;

			//Restore the neuron model state to a previous state.
			this->model->GetVectorNeuronState()->AddElapsedTime(index,-elapsedTimeInSeconds);

			//Comes form a small step and goes to a small step
			if(	CoeficientSelector[index]==2){
				CoeficientSelector[index]=1;
			}else{//goes to a smaller step.
				CoeficientSelector[index]=0;
			}
		}else{
			//We increase the state of the integration method.
			if(state[index]<BDForder){
				state[index]++;
			}

			//We acumulate these new values for the next step.
			for (int j=0; j<N_DifferentialNeuronState; j++){

				for(int i=(state[index]-1); i>0; i--){ 
					D[i][index*N_DifferentialNeuronState + j]=-D[i-1][index*N_DifferentialNeuronState + j];
				}
				D[0][index*N_DifferentialNeuronState + j]=AuxNeuronState_p[j]-NeuronState[j];
				for(int i=1; i<state[index]; i++){ 
					D[i][index*N_DifferentialNeuronState + j]+=D[i-1][index*N_DifferentialNeuronState + j];
				}
			}

			if(state[index]>1){
				for(int i=state[index]-2; i>0; i--){
					memcpy(PreviousNeuronState[i] + (index*N_DifferentialNeuronState), PreviousNeuronState[i-1] + (index*N_DifferentialNeuronState) ,sizeof(float)* N_DifferentialNeuronState);
				}
				
				memcpy(PreviousNeuronState[0] + (index*N_DifferentialNeuronState), NeuronState ,sizeof(float)* N_DifferentialNeuronState);
			}
			memcpy(NeuronState, AuxNeuronState_p ,sizeof(float)* N_DifferentialNeuronState);



			//Finaly, we evaluate the neural state variables with time dependence.
			this->model->EvaluateTimeDependentEcuation(NeuronState, elapsedTimeInNeuronModelScale, 0);

			//update the integration step size.
			PreviousIntegrationStep[index]=elapsedTimeInNeuronModelScale;

			//Set the coeficient selector to 1 for the next iteration.
			CoeficientSelector[index]=1;

			//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];

			//If the state of the cell is 0, we use a Euler method to calculate an aproximation of the solution.
			if(state[index]==0){
				this->model->EvaluateDifferentialEcuation(NeuronState, AuxNeuronState, index, bifixedElapsedTimeInNeuronModelScale);
				for (int j=0; j<N_DifferentialNeuronState; j++){
					AuxNeuronState_p[j]= NeuronState[j] + bifixedElapsedTimeInNeuronModelScale*AuxNeuronState[j];
				}
			}
			//In this case we use the value of previous states to calculate an aproximation of the solution.
			else{
				float inverseIntegrationStepFactor=PreviousIntegrationStep[index]/bifixedElapsedTimeInNeuronModelScale;
				for (int j=0; j<N_DifferentialNeuronState; j++){
					AuxNeuronState_p[j]= NeuronState[j];
					for (int i=0; i<state[index]; i++){
						AuxNeuronState_p[j]+=D[i][index*N_DifferentialNeuronState+j]/*/inverseIntegrationStepFactor*/;
					}
				}
			}

			for(int i=N_DifferentialNeuronState; i<N_NeuronStateVariables; i++){
				AuxNeuronState_p[i]=NeuronState[i];
			}

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

			float epsi=1.0f;
			int k=0;


			this->model->EvaluateDifferentialEcuation(NeuronState, AuxNeuronState1, index, bifixedElapsedTimeInNeuronModelScale);

			//This integration method is an implicit method. We use this loop to iteratively calculate the implicit value.
			//epsi is the difference between two consecutive aproximation of the implicit method. 
			while (epsi>1e-16 && k<5){
				this->model->EvaluateDifferentialEcuation(AuxNeuronState_p, AuxNeuronState, index, bifixedElapsedTimeInNeuronModelScale);

				for (int j=0; j<N_DifferentialNeuronState; j++){
					AuxNeuronState_c[j]=Coeficient[CoeficientSelector[index]][state[index]][0]*bifixedElapsedTimeInNeuronModelScale*AuxNeuronState[j] + Coeficient[CoeficientSelector[index]][state[index]][1]*bifixedElapsedTimeInNeuronModelScale*AuxNeuronState1[j] + Coeficient[CoeficientSelector[index]][state[index]][2]*NeuronState[j];
					for (int i=1; i<state[index]; i++){
						AuxNeuronState_c[j]+=Coeficient[CoeficientSelector[index]][state[index]][i+2]*PreviousNeuronState[i-1][index*N_DifferentialNeuronState + j];
					}
				}

				//jacobian.
				Jacobian(AuxNeuronState_p, jacnum, index);
			
				for(int z=0; z<N_DifferentialNeuronState; z++){
					for(int t=0; t<N_DifferentialNeuronState; t++){
						J[z*N_DifferentialNeuronState + t] = Coeficient[CoeficientSelector[index]][state[index]][0] * bifixedElapsedTimeInNeuronModelScale * jacnum[z*N_DifferentialNeuronState + t];
						if(z==t){
							J[z*N_DifferentialNeuronState + t]-=1;
						}
					}
				}

				this->invermat(J,inv_J);

				for(int z=0; z<N_DifferentialNeuronState; z++){
					float aux=0.0;
					for (int t=0; t<N_DifferentialNeuronState; t++){
						aux+=inv_J[z*N_DifferentialNeuronState+t]*(AuxNeuronState_p[t]-AuxNeuronState_c[t]);
					}
					AuxNeuronState_p1[z]=aux + AuxNeuronState_p[z];
				}

				//We calculate the difference between both aproximations.
				float aux=0.0f;
				float aux2=0.0f;
				for(int z=0; z<N_DifferentialNeuronState; z++){
					aux=fabs(AuxNeuronState_p1[z]-AuxNeuronState_p[z]);
					if(aux>aux2){
						aux2=aux;
					}
				}

				memcpy(AuxNeuronState_p , AuxNeuronState_p1 ,sizeof(float)* N_DifferentialNeuronState);

				epsi=aux2;
				k++;
			}

			//We increase the state of the integration method.
			if(state[index]<BDForder){
				state[index]++;
			}


			//We acumulate these new values for the next step.
			for (int j=0; j<N_DifferentialNeuronState; j++){

				for(int i=(state[index]-1); i>0; i--){ 
					D[i][index*N_DifferentialNeuronState + j]=-D[i-1][index*N_DifferentialNeuronState + j];
				}
				D[0][index*N_DifferentialNeuronState + j]=AuxNeuronState_p[j]-NeuronState[j];
				for(int i=1; i<state[index]; i++){ 
					D[i][index*N_DifferentialNeuronState + j]+=D[i-1][index*N_DifferentialNeuronState + j];
				}
			}

			if(state[index]>1){
				for(int i=state[index]-2; i>0; i--){
					memcpy(PreviousNeuronState[i] + (index*N_DifferentialNeuronState), PreviousNeuronState[i-1] + (index*N_DifferentialNeuronState) ,sizeof(float)* N_DifferentialNeuronState);
				}
				
				memcpy(PreviousNeuronState[0] + (index*N_DifferentialNeuronState), NeuronState ,sizeof(float)* N_DifferentialNeuronState);
			}
			memcpy(NeuronState, AuxNeuronState_p ,sizeof(float)* N_DifferentialNeuronState);



			//Finaly, we evaluate the neural state variables with time dependence.
			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;
			}

			//Set the CoeficientSelector to 1.
			CoeficientSelector[index]=1;
			PreviousIntegrationStep[index]=bifixedElapsedTimeInNeuronModelScale;

			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;
					CoeficientSelector[index]=2;
				}
			}
		}
		if(integrationMethodState[index]==1){
			integrationMethodState[index]=0;
			CoeficientSelector[index]=2;
		}

	}

	return spike;
}










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

	return out;
}	

void Bifixed_BDFn::InitializeStates(int N_neurons, float * initialization){
	if(BDForder>1){
		PreviousNeuronState = (float **)new float* [BDForder-1];
		for(int i=0; i<(BDForder-1); i++){
			PreviousNeuronState[i] = new float [N_neurons*N_DifferentialNeuronState];
		}
	}
	D = (float **)new float* [BDForder];
	for(int i=0; i<BDForder; i++){
		D[i] = new float [N_neurons*N_DifferentialNeuronState];
	}

	state = new int [N_neurons]();


	PreviousIntegrationStep=new float [N_neurons]();
	CoeficientSelector=new int [N_neurons];
	
	for(int i=0; i<N_neurons; i++){
		CoeficientSelector[i]=1;
	}


	integrationMethodCounter=new int [N_neurons]();
	integrationMethodState=new int [N_neurons]();
}


void Bifixed_BDFn::resetState(int index){
	state[index]=0;
}


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