/***************************************************************************
* BDFn_GPU2.h *
* ------------------- *
* copyright : (C) 2013 by Francisco Naveros *
* email : 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. *
* *
***************************************************************************/
#ifndef BDF_GPU2_H_
#define BDF_GPU2_H_
/*!
* \file BDFn_GPU2.h
*
* \author Francisco Naveros
* \date May 2013
*
* This file declares a class which implement six BDF (Backward Differentiation Formulas) integration methods (From
* first order to sixth order BDF integration method) in GPU. This method implements a progressive implementation of the
* higher order integration method using the lower order integration mehtod (BDF1->BDF2->...->BDF6). All integration
* methods in GPU are fixed step due to the parallel architecture of this one.
*/
#include "./IntegrationMethod_GPU2.h"
#include "../../include/neuron_model/TimeDrivenNeuronModel_GPU2.h"
//Library for CUDA
#include <helper_cuda.h>
/*!
* \class BDFn_GPU2
*
* \brief BDF1, BDF2, ..., BDF6 integration methods in GPU.
*
* This class abstracts the behavior of a Euler integration method for neurons in a
* time-driven spiking neural network.
* It includes internal model functions which define the behavior of integration methods
* (initialization, calculate next value, ...).
*
* \author Francisco Naveros
* \date May 2012
*/
class BDFn_GPU2 : public IntegrationMethod_GPU2 {
public:
/*!
* \brief These vectors are used as auxiliar vectors.
*/
float * AuxNeuronState;
float * AuxNeuronState_p;
float * AuxNeuronState_p1;
float * AuxNeuronState_c;
float * jacnum;
float * J;
float * inv_J;
//For Jacobian
float * AuxNeuronState2;
float * AuxNeuronState_pos;
float * AuxNeuronState_neg;
/*!
* \brief This constant matrix contains the coefficients of each order for the BDF integration mehtod.
*/
float * Coeficient;
// {{1.0,1.0,0.0,0.0,0.0,0.0,0.0},
// {1.0,1.0,0.0,0.0,0.0,0.0,0.0},
// {2.0/3.0,4.0/3.0,-1.0/3.0,0.0,0.0,0.0,0.0},
// {6.0/11.0,18.0/11.0,-9.0/11.0,2.0/11.0,0.0,0.0,0.0},
// {12.0/25.0,48.0/25.0,-36.0/25.0,16.0/25.0,-3.0/25.0,0.0,0.0},
// {60.0/137.0,300.0/137.0,-300.0/137.0,200.0/137.0,-75.0/137.0,12.0/137.0,0.0},
// {60.0/147.0,360.0/147.0,-450.0/147.0,400.0/147.0,-225.0/147.0,72.0/147.0,-10.0/147.0}};
/*!
* \brief This vector stores previous neuron state variable for all neuron. This one is used as a memory.
*/
float * PreviousNeuronState;
/*!
* \brief This vector stores previous neuron state variable for all neuron. This one is used as a memory.
*/
float * D;
/*!
* \brief This vector contains the state of each neuron (BDF order).
*/
int * state;
/*!
* \brief This value stores the order of the integration method.
*/
int BDForder;
/*!
* \brief Constructor of the class with 6 parameter.
*
* It generates a new BDF object in GPU memory indicating the order of the method.
*
* \param N_neuronStateVariables Number of state variables for each cell.
* \param N_differentialNeuronState Number of state variables witch are calculate with a differential equation for each cell.
* \param N_timeDependentNeuronState Number of state variables witch ara calculate with a time dependent equation for each cell.
* \param Total_N_thread Number of thread in GPU (in this method it is not necessary)
* \param Buffer_GPU This vector contains all the necesary GPU memory witch have been reserved in the CPU (this memory
* could be reserved directly in the GPU, but this suppose some restriction in the amount of memory witch can be reserved).
* \param BDForder BDF order (1, 2, ..., 6).
*/
__device__ BDFn_GPU2(TimeDrivenNeuronModel_GPU2* NewModel, int N_neuronStateVariables, int N_differentialNeuronState, int N_timeDependentNeuronState, void ** Buffer_GPU, int BDFOrder):IntegrationMethod_GPU2(NewModel, N_neuronStateVariables, N_differentialNeuronState, N_timeDependentNeuronState), BDForder(BDFOrder){
AuxNeuronState = ((float*)Buffer_GPU[0]);
AuxNeuronState_p = ((float*)Buffer_GPU[1]);
AuxNeuronState_p1 = ((float*)Buffer_GPU[2]);
AuxNeuronState_c = ((float*)Buffer_GPU[3]);
jacnum = ((float*)Buffer_GPU[4]);
J = ((float*)Buffer_GPU[5]);
inv_J = ((float*)Buffer_GPU[6]);
Coeficient=((float *)Buffer_GPU[7]);
PreviousNeuronState = ((float *)Buffer_GPU[8]);
D = ((float *)Buffer_GPU[9]);
state = ((int *)Buffer_GPU[10]);
AuxNeuronState2 = ((float*)Buffer_GPU[11]);
AuxNeuronState_pos = ((float*)Buffer_GPU[12]);
AuxNeuronState_neg = ((float*)Buffer_GPU[13]);
}
/*!
* \brief Class destructor.
*
* It destroys an object of this class.
*/
__device__ ~BDFn_GPU2(){
}
/*!
* \brief It calculate the next neural state varaibles of the model.
*
* It calculate the next neural state varaibles of the model.
*
* \param index Index of the cell inside the neuron model for method with memory (e.g. BDF).
* \param SizeStates Number of neurons
* \param Model The NeuronModel.
* \param NeuronState Vector of neuron state variables for all neurons.
* \param elapsed_time integration time step.
*/
__device__ void NextDifferentialEcuationValue(int index, int SizeStates, float * NeuronState, float elapsed_time){
int offset1=gridDim.x * blockDim.x;
int offset2=blockDim.x*blockIdx.x + threadIdx.x;
//If the state of the cell is 0, we use a Euler method to calculate an aproximation of the solution.
if(state[index]==0){
model->EvaluateDifferentialEcuation(index, SizeStates, NeuronState, AuxNeuronState);
for (int j=0; j<N_DifferentialNeuronState; j++){
AuxNeuronState_p[j*offset1 + offset2]= NeuronState[j*SizeStates + index] + elapsed_time*AuxNeuronState[j*offset1 + offset2];
}
}
//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*offset1 + offset2]= NeuronState[j*SizeStates + index];
for (int i=0; i<state[index]; i++){
AuxNeuronState_p[j*offset1 + offset2]+=D[i*SizeStates*N_DifferentialNeuronState + j*SizeStates + index];
}
}
}
for(int i=N_DifferentialNeuronState; i<N_NeuronStateVariables; i++){
AuxNeuronState_p[i*offset1 + offset2]=NeuronState[i*SizeStates + index];
}
model->EvaluateTimeDependentEcuation(offset2, offset1, AuxNeuronState_p, elapsed_time);
float epsi=1.0;
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){
model->EvaluateDifferentialEcuation(offset2, offset1, AuxNeuronState_p, AuxNeuronState);
for (int j=0; j<N_DifferentialNeuronState; j++){
AuxNeuronState_c[j*offset1 + offset2]=Coeficient[state[index]*7 + 0]*elapsed_time*AuxNeuronState[j*offset1 + offset2] + Coeficient[state[index]*7 + 1]*NeuronState[j*SizeStates + index];
for (int i=1; i<state[index]; i++){
AuxNeuronState_c[j*offset1 + offset2]+=Coeficient[state[index]*7 + i+1]*PreviousNeuronState[(i-1)*SizeStates*N_DifferentialNeuronState + j*SizeStates + index];
}
}
//jacobian.
Jacobian(offset2, offset1, model, AuxNeuronState_p, jacnum);
for(int z=0; z<N_DifferentialNeuronState; z++){
for(int t=0; t<N_DifferentialNeuronState; t++){
J[z*offset1*N_DifferentialNeuronState+ t*offset1 + offset2] = Coeficient[state[index]*7 + 0] * elapsed_time * jacnum[z*offset1*N_DifferentialNeuronState+ t*offset1 + offset2];
if(z==t){
J[z*offset1*N_DifferentialNeuronState+ t*offset1 + offset2]-=1;
}
}
}
this->invermat(offset2, offset1, 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*offset1*N_DifferentialNeuronState+ t*offset1 + offset2]*(AuxNeuronState_p[t*offset1 + offset2]-AuxNeuronState_c[t*offset1 + offset2]);
}
AuxNeuronState_p1[z*offset1 + offset2]=aux + AuxNeuronState_p[z*offset1 + offset2];
}
//We calculate the difference between both aproximations.
float aux=0.0;
float aux2=0.0;
for(int z=0; z<N_DifferentialNeuronState; z++){
aux=fabs(AuxNeuronState_p1[z*offset1 + offset2]-AuxNeuronState_p[z*offset1 + offset2]);
AuxNeuronState_p[z*offset1 + offset2]=AuxNeuronState_p1[z*offset1 + offset2];
if(aux>aux2){
aux2=aux;
}
}
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*SizeStates*N_DifferentialNeuronState + j*SizeStates + index]=-D[(i-1)*SizeStates*N_DifferentialNeuronState + j*SizeStates + index];
}
D[0*SizeStates*N_DifferentialNeuronState + j*SizeStates + index]=AuxNeuronState_p[j*offset1 + offset2]-NeuronState[j*SizeStates + index];
for(int i=1; i<state[index]; i++){
D[i*SizeStates*N_DifferentialNeuronState + j*SizeStates + index]+=D[(i-1)*SizeStates*N_DifferentialNeuronState + j*SizeStates + index];
}
}
if(state[index]>1){
for(int i=state[index]-2; i>0; i--){
for(int j=0; j<N_DifferentialNeuronState; j++){
PreviousNeuronState[i*SizeStates*N_DifferentialNeuronState + j*SizeStates + index] = PreviousNeuronState[(i-1)*SizeStates*N_DifferentialNeuronState + j*SizeStates + index];
}
}
for(int j=0; j<N_DifferentialNeuronState; j++){
PreviousNeuronState[0*SizeStates*N_DifferentialNeuronState + j*SizeStates + index] = NeuronState[j*SizeStates + index];
}
}
for(int z=0; z<N_DifferentialNeuronState; z++){
NeuronState[z*SizeStates + index] = AuxNeuronState_p[z*offset1 + offset2];
}
//Finaly, we evaluate the neural state variables with time dependence.
model->EvaluateTimeDependentEcuation(index, SizeStates, NeuronState, elapsed_time);
}
/*!
* \brief It calculate numerically the Jacobian.
*
* It calculate numerically the Jacobian.
*
* \param index Index of the cell inside the neuron model.
* \param SizeStates Number of neuron.
* \param Model neuron model.
* \param NeuronState Vector of neuron state variables for all neurons.
* \param jancum vector where is stored the Jacobian.
*/
__device__ void Jacobian(int index, int SizeStates, TimeDrivenNeuronModel_GPU2 * Model, float * NeuronState, float * jacnum){
float epsi=9.5367431640625e-7;
for (int j=0; j<N_DifferentialNeuronState; j++){
for (int i=0; i<N_NeuronStateVariables; i++){
AuxNeuronState2[i*SizeStates + index]=NeuronState[i*SizeStates + index];
}
AuxNeuronState2[j*SizeStates + index]+=epsi;
Model->EvaluateDifferentialEcuation(index, SizeStates, AuxNeuronState2, AuxNeuronState_pos);
AuxNeuronState2[j*SizeStates + index]-=2*epsi;
Model->EvaluateDifferentialEcuation(index, SizeStates, AuxNeuronState2, AuxNeuronState_neg);
for(int z=0; z<N_DifferentialNeuronState; z++){
jacnum[z*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]=(AuxNeuronState_pos[z*SizeStates + index]-AuxNeuronState_neg[z*SizeStates + index])/(2*epsi);
}
}
}
/*!
* \brief It calculate the inverse of a square matrix using Gauss-Jordan Method.
*
* It calculate the inverse of a square matrix using Gauss-Jordan Method.
*
* \param index Index of the cell inside the neuron model.
* \param SizeStates Number of neuron.
* \param a pointer to the square matrixs.
* \param ainv pointer to the inverse of the square matrixs.
*/
__device__ void invermat(int index, int SizeStates, float *a, float *ainv) {
if(N_DifferentialNeuronState==1){
ainv[0]=1/a[0];
}else{
float coef, inv_elemento;
int i,j, s;
for (i=0;i<N_DifferentialNeuronState;i++){
for(j=0;j<N_DifferentialNeuronState;j++){
if(i==j)
ainv[i*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]=1.0;
else
ainv[i*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]=0.0;
}
}
//Iteraciones
for (s=0;s<N_DifferentialNeuronState;s++)
{
inv_elemento=1.0/a[s*SizeStates*N_DifferentialNeuronState+ s*SizeStates + index];
for (j=0;j<N_DifferentialNeuronState;j++){
a[s*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]*=inv_elemento;
ainv[s*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]*=inv_elemento;
}
for(i=0;i<N_DifferentialNeuronState;i++)
{
if (i==s)
;
else
{
coef=a[i*SizeStates*N_DifferentialNeuronState+ s*SizeStates + index];
for(j=0;j<N_DifferentialNeuronState;j++){
a[i*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]-=a[s*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]*coef;
ainv[i*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]-=ainv[s*SizeStates*N_DifferentialNeuronState+ j*SizeStates + index]*coef;
}
}
}
}
}
}
/*!
* \brief It reset the state of the integration method for method with memory (e.g. BDF).
*
* It reset the state of the integration method for method with memory (e.g. BDF).
*
* \param index indicate witch neuron must be reseted.
*
*/
__device__ void resetState(int index){
state[index]=0;
}
};
#endif /* BDFN_GPU2_H_ */