// SPNET: Spiking neural network with axonal conduction delays and STDP
// Created by Eugene M. Izhikevich, May 17, 2004, San Diego, CA
// Saves spiking data each second in file spikes.dat
// To plot spikes, use MATLAB code: load spikes.dat;plot(spikes(:,1),spikes(:,2),'.');
#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h> 

#define getrandom(max1) ((rand()%(int)((max1)))) // random integer between 0 and max-1
												
const	int		Ne = 800;		// excitatory neurons			
const	int		Ni = 200;		// inhibitory neurons				 
const	int		N  = Ne+Ni;		// total number of neurons	
const	int		M  = 100;		// the number of synapses per neuron 
const	int		D  = 20;		// maximal axonal conduction delay
		float	sm = 10.0;		// maximal synaptic strength		
int		post[N][M];				// indices of postsynaptic neurons
float	s[N][M], sd[N][M];		// matrix of synaptic weights and their derivatives
short	delays_length[N][D];	// distribution of delays
short	delays[N][D][M];		// arrangement of delays   
int		N_pre[N], I_pre[N][3*M], D_pre[N][3*M];	// presynaptic information
float	*s_pre[N][3*M], *sd_pre[N][3*M];		// presynaptic weights
float	LTP[N][1001+D], LTD[N];	// STDP functions 
float	a[N], d[N];				// neuronal dynamics parameters
float	v[N], u[N];				// activity variables
int		N_firings;				// the number of fired neurons 
const int N_firings_max=100*N;	// upper limit on the number of fired neurons per sec
int		firings[N_firings_max][2]; // indices and timings of spikes

const int group_size = 10;		// size of group for stimulation
const int words = 5;			// number of words in song
int 	groups[group_size][words];	// groups of stimulatory neurons
int		delay_w[words];			// delays between onset and words
int		onset = 0;			// ms from beginning of second for onset of song

float stability[N][1000];		// stability of neurons firing across seconds

void initialize() {	int i,j,k,jj,dd, exists, r;
	for (i=0;i<Ne;i++) a[i]=0.02;// RS type
	for (i=Ne;i<N;i++) a[i]=0.1; // FS type

	for (i=0;i<Ne;i++) d[i]=8.0; // RS type
	for (i=Ne;i<N;i++) d[i]=2.0; // FS type

	for (i=0;i<N;i++) for (j=0;j<M;j++) 
	{
		do{
			exists = 0;		// avoid multiple synapses
			if (i<Ne) r = getrandom(N);
			else	  r = getrandom(Ne);// inh -> exc only
			if (r==i) exists=1;									// no self-synapses 
			for (k=0;k<j;k++) if (post[i][k]==r) exists = 1;	// synapse already exists  
		}while (exists == 1);
		post[i][j]=r;
	}
	for (i=0;i<Ne;i++)	for (j=0;j<M;j++) s[i][j]=6.0;  // initial exc. synaptic weights
	for (i=Ne;i<N;i++)	for (j=0;j<M;j++) s[i][j]=-5.0; // inhibitory synaptic weights
  	for (i=0;i<N;i++)	for (j=0;j<M;j++) sd[i][j]=0.0; // synaptic derivatives 
  	for (i=0;i<N;i++)
	{
		short ind=0;
		if (i<Ne)
		{
			for (j=0;j<D;j++) 
			{	delays_length[i][j]=M/D;	// uniform distribution of exc. synaptic delays
				for (k=0;k<delays_length[i][j];k++)
					delays[i][j][k]=ind++;
			}
		}
		else
		{
			for (j=0;j<D;j++) delays_length[i][j]=0;
			delays_length[i][0]=M;			// all inhibitory delays are 1 ms
			for (k=0;k<delays_length[i][0];k++)
					delays[i][0][k]=ind++;
		}
	}
	
  	for (i=0;i<N;i++)
	{
		N_pre[i]=0;
		for (j=0;j<Ne;j++)
		for (k=0;k<M;k++)
		if (post[j][k] == i)		// find all presynaptic neurons 
		{
			I_pre[i][N_pre[i]]=j;	// add this neuron to the list
			for (dd=0;dd<D;dd++)	// find the delay
				for (jj=0;jj<delays_length[j][dd];jj++)
					if (post[j][delays[j][dd][jj]]==i) D_pre[i][N_pre[i]]=dd;
			s_pre[i][N_pre[i]]=&s[j][k];	// pointer to the synaptic weight	
			sd_pre[i][N_pre[i]++]=&sd[j][k];// pointer to the derivative
		}
	}

	for (i=0;i<N;i++)	for (j=0;j<1+D;j++) LTP[i][j]=0.0;
	for (i=0;i<N;i++)	LTD[i]=0.0;
	for (i=0;i<N;i++)	v[i]=-65.0;		// initial values for v
	for (i=0;i<N;i++)	u[i]=0.2*v[i];	// initial values for u

	N_firings=1;		// spike timings
	firings[0][0]=-D;	// put a dummy spike at -D for simulation efficiency 
	firings[0][1]=0;	// index of the dummy spike
	
	for (i=0;i<N;i++)	for (j=0;j<1000;j++) stability[i][j]=0.0;
}

void createSong() {
	int i, j;
	for (j=0;j<words;++j) {
		for (i=0;i<group_size;++i) {
			groups[i][j] = getrandom(Ne);	// only excitatory neurons
		}
	}
	delay_w[0] = 0;		// delay since onset for word expression
	delay_w[1] = 12;
	delay_w[2] = 27;
	delay_w[3] = 46;
	delay_w[4] = 64;
}

void run() { // runs the simulation
	int		i, j, k, sec, t, word;
	float	I[N];
	FILE	*fs;
	
	//for (sec=0; sec<60*60*24; sec++)		// simulation of 1 day
	for (sec=0; sec<=60*10+7; ++sec) //simulation for 10 minutes
	{
		word = 0;
		for (t=0;t<1000;t++)				// simulation of 1 sec
		{
			for (i=0;i<N;i++) I[i] = 0.0;	// reset the input 
			
			if (sec < 60*10) {			 
				for (k=0;k<N/1000;k++)
					I[getrandom(N)]=20.0;		// random thalamic input
			
			
				/* Instead of random thalamic input, activate sets of neurons
				 * associated with a word in a song. 
				 */
				if (t>=onset && word<words) {
					if (t-onset == delay_w[word]) {
						for (k=0;k<group_size;++k) {
							I[groups[k][word]]=30.0;
						}
						word++;
					}
				}
			} else if (sec==603 || sec==606) {	// post-training testing: 2 instances of song presentation
				if (t-onset == delay_w[word]) {
					for (k=0;k<group_size;++k) {
						I[groups[k][word]]=30.0;
					}
					word++;
				}
			} else if (sec==607) {				// testing: only first word of song
				if (t-onset == delay_w[word]) {
					for (k=0;k<group_size;++k) {
						I[groups[k][word]]=30.0;
					}
				}
			}
			
			for (i=0;i<N;i++) 
			if (v[i]>=30)					// did it fire?
			{
				v[i] = -65.0;					// voltage reset
				u[i]+=d[i];					// recovery variable reset
				LTP[i][t+D]= 0.1;		
				LTD[i]=0.12;
				for (j=0;j<N_pre[i];j++) *sd_pre[i][j]+=LTP[I_pre[i][j]][t+D-D_pre[i][j]-1];// this spike was after pre-synaptic spikes
				firings[N_firings  ][0]=t;
				firings[N_firings++][1]=i;
				if (N_firings == N_firings_max) {
					printf("Too many spikes at t=%d (ignoring all)\n", t);
					N_firings=1;
				}
			}
			k=N_firings;
			while (t-firings[--k][0] <D)
			{
				for (j=0; j< delays_length[firings[k][1]][t-firings[k][0]]; j++)
				{
					i=post[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]]; 
					I[i]+=s[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]];
					if (firings[k][1] <Ne) // this spike is before postsynaptic spikes
						sd[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]]-=LTD[i];
				}
			}
			for (i=0;i<N;i++)
			{
				v[i]+=0.5*((0.04*v[i]+5)*v[i]+140-u[i]+I[i]); // for numerical stability
				v[i]+=0.5*((0.04*v[i]+5)*v[i]+140-u[i]+I[i]); // time step is 0.5 ms
				u[i]+=a[i]*(0.2*v[i]-u[i]);
				LTP[i][t+D+1]=0.95*LTP[i][t+D];
				LTD[i]*=0.95;
			}
		}
		cout << "sec=" << sec << ", firing rate=" << float(N_firings)/N << endl;
		
		if (sec<600) {
			for (i=0;i<N_firings;i++)
				if (firings[i][0] >=0)
					stability[firings[i][0]][firings[i][1]] += 1.0/600;
			/* each second, all timing-neuron pairs that fired within the second are
			 * incremented 1/600 so that if a neuron consistently fires at a certain
			 * time through all 600 seconds of training, it will have a stability 
			 * rating of 1.0
			 */ 
		}
		
		// sec from 601 to 602 is a cool-down period
		
		// analysis of data
		if (sec==603) {
			// spikes1: just the activation of second 603, presentation of noiseless song
			fs = fopen("spikes1.dat","w");
			for (i=0;i<N_firings;i++)
				if (firings[i][0] >=0)
					fprintf(fs, "%d  %d \n", firings[i][0], firings[i][1]);
			fclose(fs);
			// spikes2: after accumulation of stability rating, those that reach threshold are recorded
			fs = fopen("spikes2.dat","w");
			for (i=0;i<1000;i++) {
				for (j=0;j<N;j++) {
					if (stability[i][j] > ((j>=Ne)?.1:.06))
						fprintf(fs, "%d %d \n", i, j);
				}
			}
			// spikes3: the firings that are directly due to input
			fclose(fs);
			fs = fopen("spikes3.dat","w");
			for (word=0;word<words;++word) {
				for (k=0;k<group_size;++k) {
					fprintf(fs, "%d %d \n", delay_w[word]+3, groups[k][word]);
				}
			}
			fclose(fs);
			// stability: the raw accumulated stability rating of all 1000 neurons at every ms across training
			fs = fopen("stability.dat","w");
			for (i=0;i<N;i++) {
				for (j=0;j<1000;j++) {
						fprintf(fs, "%f ", stability[i][j]);
				}
				fprintf(fs, "\n");
			}
			fclose(fs);
			
			// the stability matrix is no longer needed, so the memory is cleared for use
			for (i=0;i<N;i++) for (j=0;j<1000;j++) stability[i][j] = 0;
			
			/* to determine which neurons are most likely indicative of having seen the song,
			 * take the timing-neuron firing pairs that occurred very often from training, match
			 * it against a presentation of the song without noise.  That eliminates the possible
			 * synfire activations that are randomly hit.  Since those would include the neurons
			 * that fired due to input, subtract away the pairs that are due to input.
			 * Calculates are worked off a point system: being above threshold gains the timing-neuron
			 * pair 1 point; being active in noiseless song presentation gains 1 point; being active due
			 * to input loses 2 points.  Hence by the end of the process, whichever pairs have 2 points
			 * are our predictor neurons.
			 */
			// noiseless song contribution
			for (i=0;i<N_firings;i++)
				if (firings[i][0] >=0) stability[firings[i][0]][firings[i][1]] += 1;
			printf("Noiseless contribution added.\n");
			// thresholded pairs contribution
			fs = fopen("spikes2.dat","r");
			while (fscanf(fs,"%d%d", &i, &j)!=EOF) {
				stability[i][j]+=1;
			}
			fclose(fs);
			printf("Stability statistics added.\n");
			// input activity subtraction
			for (word=0;word<words;++word) {
				for (k=0;k<group_size;++k) {
					stability[delay_w[word]+3][groups[k][word]]-=2;
				}
			}
			printf("Original words removed.\n");
			// spikes4: the predictor neurons
			k = 0;
			fs = fopen("spikes4.dat","w");
			for (i=0;i<1000;i++) {
				for (j=0;j<N;j++) {
					if (stability[i][j] == 2) {
						fprintf(fs, "%d %d \n", i, j);
						k++;
					}
				}
			}
			fclose(fs);
			printf("%d Predictor neurons\n",k);
		}
		
		// sec from 604 to 605 is a cool-down period
		
		if (sec==606) {
			// a full song was just shown to the network, what fired?
			int pfire = 0;
			fs = fopen("spikes5.dat","w");
			for (i=0;i<N_firings;i++) {
				if (firings[i][0] >=0) {
					fprintf(fs,"%d %d \n", firings[i][0], firings[i][1]);
					// predictor neuron positions are still in stability matrix with value 2
					// if the current firing pair would make that value 3, then a predictor 
					// neuron has fired
					if (stability[firings[i][0]][firings[i][1]] + 1 == 3) {
						pfire++;
						printf("Time: %d Neuron: %d \n", firings[i][0], firings[i][1]);
					}
				}
			}
			printf("%d fired.\n", pfire);
			fclose(fs);
		}
		if (sec==607) {
			// a single word was just shown, what fired?
			int pfire = 0;
			fs = fopen("spikes6.dat","w");
			for (i=0;i<N_firings;i++) {
				if (firings[i][0] >=0) {
					fprintf(fs,"%d %d \n", firings[i][0], firings[i][1]);
					if (stability[firings[i][0]][firings[i][1]] + 1 == 3) {
						pfire++;
						printf("Time: %d Neuron: %d \n", firings[i][0], firings[i][1]);
					}
				}
			}
			printf("%d fired.\n", pfire);
			fclose(fs);
		}

		for (i=0;i<N;i++)		// prepare for the next sec
			for (j=0;j<D+1;j++)
			LTP[i][j]=LTP[i][1000+j];
		k=N_firings-1;
		while (1000-firings[k][0]<D) k--;
		for (i=1;i<N_firings-k;i++)
		{
			firings[i][0]=firings[k+i][0]-1000;
			firings[i][1]=firings[k+i][1];
		}
		N_firings = N_firings-k;

		for (i=0;i<Ne;i++)	// modify only exc connections
		for (j=0;j<M;j++)
		{
			sd[i][j]*=0.9;
			s[i][j]+=0.01+sd[i][j];
			if (s[i][j]>sm) s[i][j]=sm;
			if (s[i][j]<0) s[i][j]=0.0;
		}
	}
}

int main() {
	initialize();	// assign connections, weights, etc. 
	createSong();	// assigns random neurons to words in song
	run();
	return 0;
}