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

#include <math.h>

const float 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}};


BDFn::BDFn(TimeDrivenNeuronModel * NewModel, int N_neuronStateVariables, int N_differentialNeuronState, int N_timeDependentNeuronState, int BDForder):FixedStep(NewModel,"BDFn", N_neuronStateVariables, N_differentialNeuronState, N_timeDependentNeuronState, true, true),BDForder(BDForder){
}

BDFn::~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;

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

	float previous_V=NeuronState[0];

	float AuxNeuronState[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 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{
		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];
			}
		}
	}

	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 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[state[index]][0]*elapsedTimeInNeuronModelScale*AuxNeuronState[j] + Coeficient[state[index]][1]*NeuronState[j];
			for (int i=1; i<state[index]; i++){
				AuxNeuronState_c[j]+=Coeficient[state[index]][i+1]*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[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++;
	}

	//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 last spike time.
	this->model->GetVectorNeuronState()->AddElapsedTime(index,elapsedTimeInSeconds);

	bool spike=this->model->EvaluateSpikeCondition(previous_V, NeuronState, index, elapsedTimeInNeuronModelScale);

	return spike;
}

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

	return out;
}	

void 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]();
}


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


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

