statisticsoutput.cpp

amygdata的神经网络算法源代码

/***************************************************************************
                          statisticsoutput.cpp  -  description
                             -------------------
    begin                : Wed Oct 24 2001
    copyright            : (C) 2001 by Matt Grover
    email                : [e-mail]
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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;
}