123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186 |
- #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;
- }
- }
- }
- }
|