ObjSim/src/spiketrain.cpp

#include "sys.hpp"   // for libcwd
#include "debug.hpp" // for libcwd

#include <iostream>
#include "spiketrain.hpp"
#include "ringbuffer.hpp"
#include "libcsim.hpp"


/** \note only use this constructor if firings don't change during the 
 lifetime of the SpikeTrain object 
*/
SpikeTrain::SpikeTrain(int** _firings, int* _N_firings, T_NNeurons* _NNeurons, const float & _dt):
	firings(_firings), 
	N_firings(_N_firings), 
	NNeurons(_NNeurons), 
	dt(_dt), 
	firings_loaded(0), 
	firings_fromLayer(_firings) 
{
	std::cout << "SpikeTrain constructor for running simulation\n";
}



/**  @brief return maximum spike frequency, calculated from sub spike trains of length NSpikes 
 * 
 * algorithm: 
 * 	for each spiking neuron store last NSpikes spikes in a ring buffer
 * 	calculate difference between last spike time and earlyest spike in ring buffer
 *	frequency is (NSpikes-a)/t_diff (NSpikes spikes have NSpikes-1 inter spike intervals) 
 * 
 * @param NSpikes 
 * @return 
 */
float SpikeTrain::MaxSpikeFreq(const int& NSpikes)
{
	if (NSpikes<2) {
		return 0;
	}
	int** firings_data=0;
	if (firings_loaded) {
		firings_data=firings_loaded;
	} else {
		firings_data=firings;
	}
	int tdiff;
	int MaxTDiff=1000;
	vector<RingBuffer<int> > vecrbuff(*NNeurons, RingBuffer<int>(NSpikes));
	//vector<int> TimeDiff(*NNeurons, MaxTDiff);
	int MinTDiff = MaxTDiff;
	RingBuffer<int>* PRngBuff;
	for (int spike=0;spike<*N_firings;++spike) {
		int CurNeuronNumber = firings_data[spike][1];
		if (CurNeuronNumber>=*NNeurons)  {
			*NNeurons=CurNeuronNumber+1;
			vecrbuff.resize(*NNeurons, RingBuffer<int>(NSpikes));
			// TimeDiff.resize(*NNeurons, MaxTDiff);
		}
		PRngBuff = &vecrbuff[CurNeuronNumber];
		PRngBuff->write(firings_data[spike][0]);	
		// cout << firings_data[spike][0] << " ";
		if (PRngBuff->isFull()) {
			tdiff = PRngBuff->last()-PRngBuff->first();
			if (tdiff<MinTDiff) {
				// 				cout << "tdiff=" << tdiff << "\n";
				// PRngBuff->print();
				MinTDiff=tdiff;
			}
		}
	}
	return 1000*(NSpikes-1)/(MinTDiff*dt);
}

SpikeTrain::SpikeTrain(const char* filename, const float& _dt): dt(_dt)
{
	N_firings=&mN_firings;
	NNeurons=&mNNeurons;
	LoadSpikeFile(filename);
}

void SpikeTrain::LoadSpikeFile(const char* filename)
{
	size_t NBytes=FileSize(filename);
	cout << "FileSize=" << NBytes << "\n";
	*N_firings=NBytes/(2*sizeof(int));
    cout << "N_firings=" << *N_firings << "\n";
	NewArray2d(firings_loaded, *N_firings, 2);
	FILE* fs = fopen(filename, "r");
	int bcount=0;
	*NNeurons=0;
	for (int spike=0;spike<*N_firings;++spike) {
		bcount=fread(&firings_loaded[spike][0], sizeof(int), 1, fs);
		if (bcount==0) {
			cout << "ERROR, spike=" << spike << "\n";
			break;
		}
		bcount=fread(&firings_loaded[spike][1], sizeof(int), 1, fs);
		if (bcount==0) {
			cout << "ERROR2, spike=" << spike << "\n";
			break;
		}
		if (firings_loaded[spike][1]>=*NNeurons) {
			*NNeurons=firings_loaded[spike][1]+1;
		}
	}
	// cout << "SpikeTrain::LoadSpikeFile: *NNeurons=" << *NNeurons << "\n";
	fclose(fs);
	if (firings_loaded) {
		firings=firings_loaded;
	};	
}

SpikeTrain::~SpikeTrain()  
{
	if (firings_loaded) {
		DeleteArray2d(firings_loaded, *N_firings);
	}
}

/**  returns the spike data
 * 
 * @param SpikeTimes spike times 
 * @param NeuronNumbers neuron numbers
 */
void SpikeTrain::SpikeData(vector<double> * SpikeTimes, vector<double> * NeuronNumbers)
{
	NeuronNumbers->resize(*N_firings);
	SpikeTimes->resize(*N_firings);
	for (int SpikeNr=0;SpikeNr<*N_firings;++SpikeNr) {
		(*NeuronNumbers)[SpikeNr]=firings[SpikeNr][1];		
		(*SpikeTimes)[SpikeNr]=firings[SpikeNr][0]*dt;
	}
}

void SpikeTrain::MeanSpikeRate(const double& WindowSize, vector<double> *timeScale, vector<double> * spikeRate)
{
	spikeRate->clear();
	int TimeStart = min(0, firings[0][0]);
	int TimeStop = firings[*N_firings-1][0]+1;
	int NTimeSteps=TimeStop-TimeStart;
	spikeRate->resize(NTimeSteps);
	timeScale->resize(NTimeSteps);
	int t=TimeStart;
	for (vector<double>::iterator it=timeScale->begin(); it!=timeScale->end(); ++it)
	{
		(*it) = t*dt;
		++t;
	}
	int WinSizeSteps = static_cast<int>(round(WindowSize/dt));
	double FreqFactor = 1000.0/(WindowSize*(*NNeurons));
	for (int spike=0;spike<*N_firings;++spike) {
		int CurSpikeTime=firings[spike][0]-TimeStart;
		int WinEnd=min(CurSpikeTime+WinSizeSteps, NTimeSteps);
		for (int step=CurSpikeTime;step<WinEnd;++step) {
			(*spikeRate)[step]+=FreqFactor;
		}
	}
}

/** calculate statistical parameters of mean neural activity (psth)
 * 
 * @param timeScale 
 * @param spikeRate 
 */
void MeanActivityParameters(vector<double> &timeScale, vector<double> & spikeRate)
{
	int nTimeSteps = spikeRate.size();
	if ((nTimeSteps<=0) || nTimeSteps!= spikeRate.size()) {
		return;
	}
	double totalTime = timeScale.back()-timeScale.front();
	double dt = totalTime/nTimeSteps;
	int totalNoSpikeTime=0;
	int maxNoSpikeTime=0;
	int curNoSpikeTime=0;
	for (int tStep=0;tStep<nTimeSteps;++tStep) {
		if (spikeRate[tStep]==0) {
			++curNoSpikeTime;
			++totalNoSpikeTime;
			if (maxNoSpikeTime<curNoSpikeTime) {
				++maxNoSpikeTime;
			}
		}
	}
}