/***************************************************************************                          statisticsoutput.cpp  -  description                             -------------------    begin                : Wed Oct 24 2001    copyright            : (C) 2001 by Matt Grover    email                : mpgrover@sourceforge.net ***************************************************************************//*************************************************************************** *                                                                         * *   This program is free software; you can redistribute it and/or modify  * *   it under the terms of the GNU General Public License as published by  * *   the Free Software Foundation; either version 2 of the License, or     * *   (at your option) any later version.                                   * *                                                                         * ***************************************************************************/#include "config.h"#include <map>#include <vector>#include <stdexcept>#include <amygdala/network.h>#include <amygdala/neuron.h>#include <amygdala/statisticsoutput.h>using namespace std;using namespace Amygdala;StatisticsOutput::StatisticsOutput():SpikeOutput(){//    combinedHistogram = 0;    maxPeakId = 0;    maxPeakTime = 0;    meanRate = 0.0;    mostActiveCount = 0;    mostActiveId = 0;    totalSpikeCount = 0;    beginTime = 0;    stepSize = 1;    logging  = false;    logFd = NULL;    traceGroups = 0;}StatisticsOutput::~StatisticsOutput(){    // TODO: delete the elements of histogram and clear the    // map and do the same with outputHistory    combinedHistogram.clear();    if(logFd && logging) fclose(logFd);}void StatisticsOutput::OutputEvent(Neuron* nrn, AmTimeInt eventTime){    if(logging) Log(nrn, eventTime);    vector<AmTimeInt>& hist = outputHistory[nrn->GetId()];    hist.push_back(eventTime);    calcTime = eventTime;}void StatisticsOutput::AddTrace(unsigned int groupId){	if (groupId >= (sizeof(AmGroupInt)*8)) {		throw runtime_error("GroupId is too large");	}	unsigned int state=1;	state <<= groupId;	traceGroups |= state;}void StatisticsOutput::ClearHistory(){    map< AmIdInt, vector<unsigned int> >::iterator histItr;    histItr = histogram.begin();    while ( histItr != histogram.end() ) {        histItr->second.clear();        histItr++;    }    histogram.clear();    map< AmIdInt, vector<AmTimeInt> >::iterator outItr;    outItr = outputHistory.begin();    while ( outItr != outputHistory.end() ) {        outItr->second.clear();        outItr++;    }    outputHistory.clear();    combinedHistogram.clear();    beginTime = Network::GetNetworkRef()->SimTime();}vector<unsigned int>& StatisticsOutput::Histogram(){    // TODO: Check this function for correctness.    unsigned int i, eventCount = 0, eventIdx = 0, numElem = 0;    AmTimeInt endStep, binSize;    if (calcTime != Network::GetNetworkRef()->SimTime()) {        calcTime = Network::GetNetworkRef()->SimTime();        lastCalcTime = calcTime;        combinedHistogram.clear();    }    if (!combinedHistogram.size()) {        // find the bin size and number of elements in combinedHistogram        binSize = stepSize * 1000;        numElem = (calcTime - beginTime) / binSize;        // initialize combinedHistogram        for (i=0; i<numElem; i++) {            combinedHistogram.push_back(0);        }        maxPeakId = 0;        maxPeakTime = 0;        // Fill the histogram based on the raw data in outputHistory.        // This will be for all neurons, so use an iterator        map< AmIdInt, vector<AmTimeInt> >::iterator histItr;        histItr = outputHistory.begin();        while ( histItr != outputHistory.end() ) {            vector<AmTimeInt>& events = histItr->second;            unsigned int peakCount = 0;            eventCount = 0;            eventIdx = 0;            for (i=0; i<combinedHistogram.size(); i++) {                endStep = (i + 1) * binSize + beginTime;                if (eventIdx >= events.size()) {                    break;                }                while (events[eventIdx] < endStep) {                    eventCount++;                    if (++eventIdx >= events.size()) {                        break;                    }                }                combinedHistogram[i] += eventCount;                if (eventCount > peakCount) {                    peakCount = eventCount;                    maxPeakId = histItr->first;                    maxPeakTime = (i * binSize + beginTime) / 1000;                }                if (events.size() > mostActiveCount) {                    mostActiveCount = events.size();                    mostActiveId = histItr->first;                }                eventCount = 0;            }            histItr++;        }    }    return combinedHistogram;}vector<unsigned int>& StatisticsOutput::Histogram(AmIdInt neuronId){    // TODO: Check this function for correctness.    unsigned int i, eventCount = 0, eventIdx = 0, numElem = 0;    AmTimeInt endStep, binSize;    // clear out the vectors in histogram if they are old    if (calcTime != Network::GetNetworkRef()->SimTime()) {        calcTime = Network::GetNetworkRef()->SimTime();        lastCalcTime = calcTime;        map< AmIdInt, vector<unsigned int> >::iterator itr;        itr = histogram.begin();        while (itr != histogram.end()) {            itr->second.clear();            itr++;        }    }    vector<unsigned int>& nidHistogram = histogram[neuronId];    if (!nidHistogram.size()) {        // find the bin size and number of elements in combinedHistogram        binSize = stepSize * 1000;        numElem = (calcTime - beginTime) / binSize;        // initialize combinedHistogram        for (i=0; i<numElem; i++) {            nidHistogram.push_back(0);        }        vector<AmTimeInt>& events = outputHistory[neuronId];        eventCount = 0;        eventIdx = 0;        for (i=0; i<nidHistogram.size(); i++) {            endStep = (i + 1) * binSize + beginTime;            if (eventIdx >= events.size()) {                break;            }            while (events[eventIdx] < endStep) {                eventCount++;                if (++eventIdx >= events.size()) {                    break;                }            }            nidHistogram[i] += eventCount;            eventCount = 0;        }    }    return nidHistogram;}vector<unsigned int> StatisticsOutput::Histogram(AmTimeInt start, AmTimeInt end){    // TODO: Check this function for correctness.    unsigned int i, eventCount = 0, eventIdx = 0, numElem = 0;    AmTimeInt endStep, binSize;    vector<unsigned int> intervalHist;    if (start >= end) {        return intervalHist;    }    // find the bin size and number of elements in combinedHistogram    binSize = stepSize * 1000;    numElem = (end - start) / binSize;    // initialize combinedHistogram    for (i=0; i<numElem; i++) {        intervalHist.push_back(0);    }    // Fill the histogram based on the raw data in outputHistory.    // This will be for all neurons, so use an iterator    map< AmIdInt, vector<AmTimeInt> >::iterator histItr;    histItr = outputHistory.begin();    while ( histItr != outputHistory.end() ) {        vector<AmTimeInt>& events = histItr->second;        eventCount = 0;        eventIdx = 0;        for (i=0; i<intervalHist.size(); i++) {            endStep = (i + 1) * binSize + start;            if (eventIdx >= events.size()) {                break;            }            while (events[eventIdx] < endStep) {                eventCount++;                if (++eventIdx >= events.size()) {                    break;                }            }            intervalHist[i] += eventCount;            eventCount = 0;        }        histItr++;    }    return intervalHist;}unsigned int StatisticsOutput::TotalOutputSpikes(){    unsigned int numEvents = 0;    if (!combinedHistogram.size()) {        Histogram();    }    map< AmIdInt, vector<AmTimeInt> >::iterator outItr;    outItr = outputHistory.begin();    while ( outItr != outputHistory.end() ) {        numEvents += outItr->second.size();        outItr++;    }    return numEvents;}unsigned int StatisticsOutput::TotalOutputSpikes(AmIdInt neuronId){    if (!combinedHistogram.size()) {        Histogram();    }    unsigned int numEvents = outputHistory[neuronId].size();    return numEvents;}float StatisticsOutput::MeanSpikeRate(){    unsigned int numEvents = 0;    AmTimeInt endTime = Network::GetNetworkRef()->SimTime();    float avg;    map< AmIdInt, vector<AmTimeInt> >::iterator outItr;    outItr = outputHistory.begin();    while ( outItr != outputHistory.end() ) {        numEvents += outItr->second.size();        outItr++;    }    if ( (endTime - beginTime) > 0 ) {        avg = ( float(numEvents) / float(endTime - beginTime) ) * 1000000.0;    }    else {        avg = 0.0;    }    return avg;}float StatisticsOutput::MeanSpikeRate(AmIdInt neuronId){    float avg;    AmTimeInt endTime = Network::GetNetworkRef()->SimTime();    vector<AmTimeInt>& hist = outputHistory[neuronId];    unsigned int numEvents = hist.size();    if ( (endTime - beginTime) > 0 ) {        avg = ( float(numEvents) / float(endTime - beginTime) ) * 1000000.0;    }    else {        avg = 0.0;    }    return avg;}float StatisticsOutput::PeakSpikeRate(){    unsigned int i;    // TODO: This algorithm is good enough for temporary use,    // but a better one needs to be produced.  If a small    // stepSize is in use, this algorithm will give wildly    // inacurate results (1 spike/1 ms step = 1000 spikes/second)    if (!combinedHistogram.size()) {        Histogram();    }    for (i=0; i<combinedHistogram.size(); i++) {        // find combined peak rate and combined peak time    }    return combinedPeakRate;}	float StatisticsOutput::PeakSpikeRate(AmIdInt neuronId){    // TODO: This algorithm is good enough for temporary use,    // but a better one needs to be produced.  If a small    // stepSize is in use, this algorithm will give wildly    // inacurate results (1 spike/1 ms step = 1000 spikes/second)    vector<unsigned int> nHist = Histogram(neuronId);    unsigned int i, maxCount = 0;    for (i=0; i<nHist.size(); i++) {        if (nHist[i] > maxCount) {            maxCount = nHist[i];        }    }    float rate = (float(maxCount) / float(stepSize)) * 1000.0;    return rate;}	AmTimeInt StatisticsOutput::PeakRateTime(){    if (!combinedHistogram.size()) {        Histogram();    }    if (combinedPeakRate == 0.0) {        PeakSpikeRate();    }    return combinedPeakTime;}	AmTimeInt StatisticsOutput::PeakRateTime(AmIdInt neuronId){    vector<unsigned int> nHist = Histogram(neuronId);    unsigned int i, maxIndex = 0, maxCount = 0;    for (i=0; i<nHist.size(); i++) {        if (nHist[i] > maxCount) {            maxCount = nHist[i];            maxIndex = i;        }    }    return maxIndex * stepSize;}	AmIdInt StatisticsOutput::PeakNeuron(){    if (!combinedHistogram.size()) {        Histogram();    }    return maxPeakId;}	AmIdInt StatisticsOutput::MostActiveNeuron(){    if (!combinedHistogram.size()) {        Histogram();    }    return mostActiveId;}	void StatisticsOutput::SetStepSize(AmTimeInt step){    if (step != stepSize) {        combinedHistogram.clear();        map< AmIdInt, vector<unsigned int> >::iterator histItr;        histItr = histogram.begin();        while ( histItr != histogram.end() ) {            histItr->second.clear();            histItr++;        }        histogram.clear();    }    stepSize = step;}void StatisticsOutput::LogSpikeTimes(string filename, AmTimeInt start, AmTimeInt end){    logFd = fopen(filename.c_str(), "w");    if(!logFd) throw "Canot open file: " + filename;    logging = true;    logStart = start;    logEnd   = end;}void StatisticsOutput::Log(Neuron* nrn, AmTimeInt eventTime){    if(eventTime > logEnd){        fclose(logFd);        logging = false;        return;    }    	bool traceNeuron = (traceGroups & nrn->GetOutputGroupIndex());    if(traceNeuron){        fprintf(logFd, "%ld %ld\n", eventTime/1000, nrn->GetId());    }}void StatisticsOutput::CloseLog(){    if(logFd && logging) fclose(logFd);    logging = false;    logFd = NULL;}