CSpikeStats.cpp

Go to the documentation of this file.

#include "CSpikeStats.h"
#include "tsfuncs.h"
#include <numeric>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;

void CSpikeStats::AnalyseSpikes()
{
        
        double timediff;
        double prediff, postdiff;
        double currentfraction;
        titdata currentdata;
        titdata currenttime;
        titdata lasttime;
        titdata next,last;
        vector<int>::iterator currentoccbin;
        vector<double>::iterator currentvalbin;
        
        currentdata = Data->begin()+1;
        currenttime = Time->begin()+1;
        lasttime = currenttime - 1;
        next = currentdata + 1;
        last = currentdata -1;
        CreateBins();
        DerivTotalOcc = 0;
        AbsoluteTotalOcc = 0;
        cout << " Analyzing spikes "  << flush;
        while (distance(next,Data->end()) > 1)
        {
                timediff = *currenttime - *lasttime;
                prediff = *currentdata - *last;
                postdiff = *next -  *currentdata;
                
                if ((timediff < 2 * samplerate) && IsSpikeAbsolute(prediff,postdiff,*currentdata,-1e15) )
                {
                        currentoccbin = AbsoluteValueOcc.begin();
                        currentvalbin = AbsoluteValueBin.begin();
                        while ( *currentdata > *currentvalbin)
                        {
                                currentoccbin++;
                                currentvalbin++;
                        }
                        *currentoccbin += 1;
                        AbsoluteTotalOcc += 1;
                        currentoccbin = DerivativeOcc.begin();
                        currentvalbin = DerivativeBin.begin();
                        while ( (fabs(prediff) >= *currentvalbin) && (currentoccbin != DerivativeOcc.end()))
                        {
                                currentoccbin++;
                                currentvalbin++;
                        }
                        *currentoccbin += 1;
                        DerivTotalOcc += 1;
                }
                lasttime ++;
                currentdata++;
                currenttime++;
                next++;
                last++;
        }
        currentfraction = 0;
        currentoccbin = DerivativeOcc.end() -1;
        currentvalbin = DerivativeBin.end()-1;
        while (currentfraction < HeightFraction)
        {
                currentfraction += *currentoccbin/DerivTotalOcc;
                currentoccbin--;
                currentvalbin--;
        }
        HeightThreshold = *currentvalbin;
        cout << HeightFraction * 100 << " percent threshold: " << HeightThreshold << endl;
        cout << "   .... done " << endl << flush;
        cout << "Total Number Deriv: " << DerivTotalOcc << " Absolute: " << AbsoluteTotalOcc << endl;
}

void CSpikeStats::CreateBins()
{
        double timediff;
        double prediff, postdiff;
        
        int nbins;
        titdata currentdata;
        titdata currenttime;
        titdata lasttime;
        titdata next,last;
        
        currentdata = Data->begin()+1;
        currenttime = Time->begin()+1;
        lasttime = currenttime - 1;
        next = currentdata +1;
        last = currentdata -1;
        absolutemax = 0;
        derivmax = 0;
        absolutemin=1e19;
        derivmin=1e19;
        if (AbsoluteValueBin.size() != 0)
        {
                AbsoluteValueBin.clear();
                DerivativeBin.clear();
        }
        while ( next != Data->end() )
        {
                timediff = *currenttime - *lasttime;
                prediff = *currentdata - *last;
                postdiff = *next -  *currentdata;
                if ( (timediff < 2 * samplerate) && IsSpikeAbsolute(prediff,postdiff,*currentdata,-1e15) )
                {
                        if ( *currentdata > absolutemax)
                                absolutemax = *currentdata;
                        //if ( fabs(*currentdata) < absolutemin)
                        //      absolutemin = fabs(*currentdata);
                        if ( *currentdata < absolutemin)
                                absolutemin = *currentdata;
                        if (prediff > derivmax)
                                derivmax = prediff;
                        if (prediff < derivmin)
                                derivmin = prediff; 
                }
                lasttime ++;
                currentdata++;
                currenttime++;
                next++;
                last++;
        }
        derivstep = (derivmax - derivmin)/nDerivBins + 0.01;
        absolutestep = (absolutemax-absolutemin)/nValueBins + 0.01;
        cout << "AbsMax: " << absolutemax << endl;
        cout << "Absmin: " << absolutemin << endl;
        cout << "DerivMax: " << derivmax << endl;
        cout << "DerivMin: " << derivmin << endl;
        cout << "Step  Absolute: " << absolutestep << " Deriv: " << derivstep << endl;
        cout << "nBins Absolute: " << nValueBins   << " Deriv: " << nDerivBins << endl; 
        for (int i = 0; i < nValueBins; ++i)
                AbsoluteValueBin.push_back(absolutemin + (i+1) * absolutestep);
        for (int i = 0; i < nDerivBins; ++i)
                DerivativeBin.push_back(derivmin + (i+1) * derivstep);
        AbsoluteValueOcc.assign(nValueBins,0);
        DerivativeOcc.assign(nValueBins,0);
        cout << "First Bin Deriv: "<< DerivativeBin.front() << " Last Bin Deriv : " << DerivativeBin.back() << endl << flush;
}

int CSpikeStats::WriteStats(std::string filename)
{
        if (AbsoluteValueBin.empty() || DerivativeBin.empty())
        {
                cerr << "No statistics gathered yet !";
                return 100;
        }
        
        ofstream outfile;
        
        outfile.open((filename+".deriv_stats").c_str());
        for (int i = 0; i < DerivativeBin.size(); ++i)
        {
                outfile.precision(3);
                outfile << setw(15) << setfill(' ') << DerivativeBin.at(i);
                outfile <<   setw(15) << setfill(' ') << DerivativeOcc.at(i)/DerivTotalOcc << endl;
        }
        outfile.close();
        outfile.open((filename+".val_stats").c_str());
        for (int i = 0; i < AbsoluteValueBin.size(); ++i)
        {
                outfile.precision(3);
                outfile << setw(15) << setfill(' ')<< AbsoluteValueBin.at(i);
                outfile << setw(15) << setfill(' ')<< AbsoluteValueOcc.at(i)/AbsoluteTotalOcc << endl;
        }
        outfile.close();
        return 0;
}
CSpikeStats::CSpikeStats()
{
        HeightFraction = 0;
        HeightThreshold = 0;
        Data = NULL;
        Time = NULL;
}
CSpikeStats::~CSpikeStats()
{}