thetaNeuron.cpp

// Summary: Theta neuron class function definitions, see thetaNeuron header for more details
// Author: Sam McKennoch
// Date: 2007-2009

#include "stdafx.h"
#include "thetaNeuron.h"

using namespace std;

//Display to Console the Theta Neuron Properties
void thetaNeuron::displayThetaNeuron(void) {
        cout << endl;
        cout << "Neuron Properties: " << endl;
        cout << "\tBeta: " << beta << endl;
        cout << "\tIo: " << io << endl;
        cout << "\tPositive Fixed Point: " << positive_fixed_point << endl;
        cout << "\tPhase: " << phase << endl;
        cout << "\tNumber of Inputs: " << num_inputs << endl;                       
        cout << "\tWeights: ";  
        showArray<double>(weights,num_inputs);                     
        cout << endl << endl;                       
}                                              

//Simulate Theta Neuron
int thetaNeuron::runThetaNeuron(double ti[], double ts[], int max_num_ts=MAX_NUM_OUTPUT_SPIKES) {  
        double thm; //The Phase Right Before a Specific Input Spike is Received
        double thp; //The Phase Right After a Specific Input Spike is Received
        double t_bl=(1/beta)*atanh(beta/tan((positive_fixed_point+OFFSET)/2)); //Baseline Firing Time
        int num_output_spikes=0; //Number of Output Spikes
        int* indices = NULL;
        indices = new int[num_inputs]; //Indices of Input Spikes

        //Copy ti array so original is not overwritten
        double* ti_copy = NULL;
        ti_copy = new double[num_inputs];
        for (int j=0;j<num_inputs;j++) {
                ti_copy[j]=ti[j];
        }

        //Sort Input Spike Times To Speed Up Processing
        sortArray(ti_copy,indices,num_inputs);

        if (verbose==true) {
                cout << "Input Spike Times: ";
                showArray<double>(ti_copy,num_inputs);
                cout << "Input Spike Indices: ";
                showArray<int>(indices,num_inputs);
                cout << "Baseline Firing Time: " << t_bl << endl;
        }

        if (method==EventDriven) {

                double prev_time=0; //The Previous Time at which an Input Spike Arrived
                double prev_thp=phase; //The Phase Directly after the Previous Input Arrived
                double c1, c2; //Temporary Coefficients

                for (int j=0; j<num_inputs; j++) {
                        c1=(1/beta)*tan(prev_thp/2);
                        if (j==0) {
                                c2=beta*(c1+tanh(-beta*ti_copy[j]))/(1+c1*tanh(-beta*ti_copy[j]));
                        }
                        else {
                                c2=beta*(c1+tanh(-beta*(ti_copy[j]-ti_copy[j-1])))/(1+c1*tanh(-beta*(ti_copy[j]-ti_copy[j-1])));     
                        }
                        thm=2*atan(c2);
                        thp=2*atan(ALPHA*weights[indices[j]]+c2);
                        if (verbose==true) {
                                cout << "--->Input Spike Time " << j << " (Neuron " << indices[j] << "): " << ti_copy[j] << endl;                                            
                                cout << "--->Pre -Spike Phase: " << thm << endl;
                                cout << "--->Post-Spike Phase: " << thp << endl;
                                cout << "--->Weight: " << weights[indices[j]] << endl;
                                if (thp<positive_fixed_point && thp>-positive_fixed_point) {
                                        cout << "--->Operating Below Positive Fixed Point (" << positive_fixed_point << ") -  Neuron may not fire again" << endl;
                                }                                
                        }

                        //Check for firing on spike
                        if (thp>PI && 2*atan(c2)<=PI) {
                                ts[num_output_spikes]=ti_copy[j];
                                num_output_spikes++;
                                thp=thp-2*PI;
                                if (verbose==true) {
                                        cout << "---> Neuron Output Spike Time " << num_output_spikes << ": " << ts[num_output_spikes-1] << endl;
                                }                          
                        }

                        //Check for firing on trajectory; If Io<0 the previous thp must be
                        //between the positive fixed point and pi and thm must be between 
                        //-pi and the negative fixed point in order for the neuron to have 
                        //fired on trajectory, regardless of the weight.
                        else if (prev_thp<PI && prev_thp>positive_fixed_point && 2*atan(c2)>-PI && 2*atan(c2)<=(-positive_fixed_point+10*EPS)) {
                                ts[num_output_spikes]=prev_time+(1/beta)*atanh(1/c1);
                                num_output_spikes++;
                                if (verbose==true) {
                                        cout << "---> Neuron Output Spike Time " << num_output_spikes << ": " << ts[num_output_spikes-1] << endl;
                                }                            
                        }
                        prev_time=ti_copy[j];
                        prev_thp=thp;                    

                        if (num_output_spikes==max_num_ts) {
                                return num_output_spikes;
                        }
                }

                //Calculate thm later to save time, unless needed to spike checking
                if (num_inputs>0)
                        c1=tan(thp/2)/beta;
                else
                        c1=tan(phase/2)/beta;

                //If phase is between the positive fixed point and pi, calculate the final ts
                //with or without previous input spikes
                if (num_inputs>0 && (thp>(positive_fixed_point+1e-5))) {
                        ts[num_output_spikes]=ti_copy[num_inputs-1]+(1/beta)*atanh(1/c1);
                        num_output_spikes++;                                
                        if (verbose==true) {
                                cout << "---> Neuron Output Spike Time " << num_output_spikes << ": " << ts[num_output_spikes-1] << endl;
                        }                                                        
                }
                else if (num_inputs==0 && (phase>(positive_fixed_point+1e-5))) {
                        ts[num_output_spikes]=(1/beta)*atanh(1/c1);
                        num_output_spikes++;
                        if (verbose==true) {
                                cout << "---> Neuron Output Spike Time " << num_output_spikes << ": " << ts[num_output_spikes-1] << endl;
                        }                                                        
                }
        }

        else if (method==Numerical) {
                int f_count=0; //The number of input spikes that have been processed
                double theta_now=phase; //The current value of neuron phase
                double theta_past=phase; //The previous value of neuron phase
                double int_limit; //The maximum number of numerical integration iterations

                //Integrate long enough to include effects of all input spikes
                if (num_inputs>0) {int_limit=(ti_copy[num_inputs-1]+1.5*t_bl)/TIMESTEP;}
                else {int_limit=1.1*t_bl/TIMESTEP;}

                for (int j=0;j<int_limit;j++) {
                        theta_now=theta_past + TIMESTEP * ((1-cos(theta_past))+ALPHA*io*(1+cos(theta_past)));
                        for (int k=f_count;k<num_inputs;k++) {
                                if ((ti_copy[k]>=j*TIMESTEP) && (ti_copy[k]<((j+1)*TIMESTEP))) {
                                        if (verbose==true) {
                                                cout << "--->Input Spike Time " << k << " (Neuron " << indices[k] << "): " << ti_copy[k] << endl;
                                        }  
                                        thm=theta_now;
                                        if (verbose==true) {
                                                cout << "--->Pre -Spike Phase: " << thm << endl;
                                        }
                                        theta_now=2*atan(ALPHA*weights[indices[k]]+tan(theta_now/2));
                                        thp=theta_now;
                                        if (verbose==true) {
                                                cout << "--->Post-Spike Phase: " << thp << endl;
                                                cout << "--->Weight: " << weights[indices[k]] << endl;
                                        }
                                        f_count++;
                                }
                                else
                                        break;
                        }

                        //Check to see if no further spikes are possible/desired
                        if (f_count==num_inputs && theta_now<positive_fixed_point) {
                                break;
                        }

                        //Check for output spike
                        if (theta_now>PI && theta_past<=PI) {
                                num_output_spikes++;
                                ts[num_output_spikes-1]=j*TIMESTEP;
                                if (verbose==true) {
                                        cout << "---> Neuron Output Spike Time " << num_output_spikes << ": " << ts[num_output_spikes-1] << endl;
                                }
                                theta_now=theta_now-2*PI;
                                if (num_output_spikes==max_num_ts) {
                                        return num_output_spikes;
                                }
                        }

                        //Special Case of Refractory Firing
                        if (theta_now>-PI && theta_past<=-PI) {
                                theta_now=theta_now+2*PI;
                        }

                        //Prepare for Next Iteration
                        theta_past=theta_now;
                }
        }

        if (verbose==true) {
                cout << endl;
        }

        //Delete Dynamically Allocated Variables
        delete [] indices;
        indices = NULL;
        delete [] ti_copy;
        ti_copy = NULL;

        return num_output_spikes;


}

---