basicneuron.cpp

此代码经过大量使用

/***************************************************************************
                          basicneuron.cpp  -  description
                             -------------------
    copyright            : (C) 2000, 2001, 2002 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.                                   *
 *                                                                         *
 ***************************************************************************/

using namespace std;
#include <cmath>
#include <iostream.h>
#include "basicneuron.h"
#include "functionlookup.h"


BasicNeuron::BasicNeuron(AmIdInt neuronId):
    Neuron(neuronId)
{
    // prevent risk of deleting non-existant arrays in derived classes
    pspLookup = 0;
    dPspLookup = 0;
}

BasicNeuron::~BasicNeuron()
{
}

void BasicNeuron::CalcState(float& state, float& deriv,
                        const AmTimeInt& calcTime,
                        const AmTimeInt& historySize)
{
    AmTimeInt i, funcTime, tblIndex;
    float funcWeight, stepSizeInv = 1.0/pspStepSize;
    InputHist input;

    state = 0.0;
    deriv = 0.0;
    for (i=histBeginIdx; i<historySize; i++) {
        input = inputHist[i];
        funcTime = calcTime - input.time;
        funcWeight = input.weight;

        #ifdef DEBUG_NEURON
        cout << "funcTime: " << funcTime << endl;
        cout << "funcWeight: " << funcWeight << endl;
        #endif

        tblIndex = (unsigned int)(funcTime * stepSizeInv);
        if (tblIndex < pspLSize) {  // don't go beyond the end of the tables
            state = state + (funcWeight * (pspLookup[tblIndex]));
            deriv = deriv + (funcWeight * (dPspLookup[tblIndex]));
        }
        else {
            ++histBeginIdx;
        }
    }
}


void BasicNeuron::InputSpike(SynapseItr& inSynapse, AmTimeInt inTime, unsigned int numSyn = 0)
{
    // If the neuron is within a refractory period, either ignore the
    // input (training off) or call Train(inTime)
    if ( (inTime - spikeTime) <= refPeriod ) {

        // Don't do anything if spikeTime == 0 because this neuron
        // has not yet fired and the simulation time could be less
        // than the refractory period.
        if (spikeTime) {
            if ( IsTraining() ) {
                Train(inTime);
            }
            return;
        }
    }

    unsigned int i, iterate, converged, histSize, tblIndex;
    AmTimeInt calcTime, funcTime;
    float currState = 0.0;
    float currDeriv = 0.0;
    float stateDelta = 0.0;
    float threshCrs = 0.0;
    float lstThreshCrs = 0.0;
    InputHist tmpInput;

    NEvent eventTp = NOACTION;
    iterate = 1;
    converged = 0;
    calcTime = 0;
    funcTime = 0;
    histSize = 0;
    tblIndex = 0;

    if (!maxThreshCrs) {
        maxThreshCrs = pspLSize * pspStepSize;

        // find the convergence resolution (the resolution at which two values
        // of threshCrs are considered to be identical) -- must be <= simStepSize
        if (simStepSize > pspStepSize) {
            if (pspStepSize > (simStepSize / 2.0)) {
                convergeRes = pspStepSize;
            }
            else {
                convergeRes = simStepSize / 2.0;
            }
        }
        else {
            convergeRes = simStepSize;
        }
    }

    calcTime = inTime;
    inputTime = inTime;
    currTime = inTime;

    if (numSyn) {
        tmpInput.weight = 0.0;
        float synWeight = 0.0;
        for (unsigned int j=0; j<numSyn; ++j) {
            synWeight = inSynapse[j]->GetWeight();
            tmpInput.weight = tmpInput.weight + synWeight;
        }
    }
    else {
        tmpInput.weight = inSynapse[0]->GetWeight();
    }
    tmpInput.time = inTime;
    inputHist.push_back(tmpInput);

    histSize = inputHist.size();

    if (trainingMode) {
        SynapseHist synHist;
        synHist.time = inTime;

        if (numSyn) {
            for (unsigned int i=0; i<numSyn; ++i) {
                synHist.syn = inSynapse[i];
                synapseHist.push_back(synHist);
            }
        }
        else {
            synHist.syn = inSynapse[0];
            synapseHist.push_back(synHist);
        }
    }


    // Return if the neuron is going to spike during the current time
    // step and the input weight is positive
    if (inTime == schedSpikeTime) {
        if (tmpInput.weight > 0.0) {
            return;
        }
    }

    #ifdef DEBUG_NEURON
    cout << "\nNeuron " << nId << " receiving a spike from " << inNeuron << endl;
    for (i=0; i<histSize; i++) {
        tmpInput = inputHist[i];
        cout << "inTimeHist[" << i << "] " << tmpInput.time << endl;
        cout << "inWeightHist[" << i << "] " << tmpInput.weight << endl;
    }

    cout << "calcTime: " << calcTime << endl;
    cout << "inTime: " << inTime << endl;
    cout << "currTime: " << currTime << endl;
    #endif

    i = 0;

    // Use Newton's method to determine if and when the spike will occur
    /**************************************************************************
    *    1) Determine the membrane potential (currState) at time (calcTime -
    *       inTimeHist[i]) by summing the state of each inTimeHist[]
    *       (use pspLookup).
    *    2) Find the derivative of the function for calcTime (dPspLookup).
    *    3) Calculate intercept with thresholdPtnl.
    *    4) Set new calcTime to time of intercept.
    *    5) Repeat until:
    *        a) Two successive iterations result in a change in calcTime
    *           smaller than convergeRes (the convergence resolution).
    *            (Converges)
    *        b) The derivative of the function becomes negative.
    *            (Does not converge)
    **************************************************************************/
    // Take care of a special case first:
    // If this is the first input (starting from a rest potential)
    // then the PSP can only cross the threshold at t = synTimeConst
    // because of constraints on the maximum value of weights.
    #ifdef DEBUG_NEURON
    cout << "Starting main loop...\n";
    #endif

    if (histSize == 1) {
        iterate = 0;
        funcTime = int( synTimeConst * 1000.0 );
        tblIndex = (funcTime / pspStepSize);
        currState = (inputHist[0].weight) * (pspLookup[tblIndex]);
        if (currState >= thresholdPtnl) {
            calcTime += funcTime;
            converged = 1;
        }
        else {
            converged = 0;
        }
    }
    else {
        // find a good starting point for N's method
//        calcTime = inTime + int( memTimeConst * 1000.0 );
//
//        while (iterate) {
//            CalcState(currState, currDeriv, calcTime, histSize);
//
//            if ( (currDeriv < 0.0) || (currState > thresholdPtnl) ) {
//                // went too far -- choose a smaller starting point
//                calcTime -= ( ( calcTime - inTime ) / 2 );
//            }
//            else {
//                // good starting point
//                iterate = 0;
//            }
//        }
//
//        iterate = 1;
    }

    lstThreshCrs = float(calcTime);
    while (iterate) {
        currState = 0.0;
        currDeriv = 0.0;
        // debug msg

        // calcTime has to be rounded to the nearest simStepSize
        RoundTime(calcTime);

        #ifdef DEBUG_NEURON
        cout << "calcTime: " << calcTime << endl;
        #endif

        // should be inline
        CalcState(currState, currDeriv, calcTime, histSize);
//        for (i=0; i<histSize; i++) {
//            funcTime = calcTime - inTimeHist[i];
//            funcWeight = inWeightHist[i];
//
//            // debugging msg
//            cout << "funcTime: " << funcTime << endl;
//            cout << "funcWeight: " << funcWeight << endl;
//            // end debugging
//
//            tblIndex = (funcTime / pspStepSize);
//            currState = currState + (funcWeight * (pspLookup[tblIndex]));
//            currDeriv = currDeriv + (funcWeight * (dPspLookup[tblIndex]));
//        }
        #ifdef DEBUG_NEURON
        cout << "currState: " << currState << endl;
        cout << "currDeriv: " << currDeriv << endl;
        #endif

        if ( (currDeriv < 0.0) && (currState < thresholdPtnl) ) {
            converged = 0;
            iterate = 0;

            // set the scheduled spike time to zero to keep
            // the neuron from spiking when an inhibitory input
            // spike has canceled a previously scheduled spike.
            schedSpikeTime = 0;
        }
        else if (currState > thresholdPtnl) {
            converged = 1;
            iterate = 0;
        }
        else {
            stateDelta = fabs(thresholdPtnl - currState);
            threshCrs = (stateDelta / currDeriv) + lstThreshCrs;

            #ifdef DEBUG_NEURON
            cout << "stateDelta: " << stateDelta << endl;
            cout << "threshCrs: " << threshCrs << endl;
            cout << "lstThreshCrs: " << lstThreshCrs << endl;
            #endif

            if ( (threshCrs - lstThreshCrs) < convergeRes) {
                converged = 1;
                iterate = 0;

                #ifdef DEBUG_NEURON
                cout << "threshCrs - lstThreshCrs < convergeRes\n";
                cout << "Converged\n";
                #endif
            }
            else {
                converged = 0;
                if (threshCrs > (maxThreshCrs + lstThreshCrs)) {
                    iterate = 0;
                }
                calcTime = int(threshCrs);
                lstThreshCrs = threshCrs;

                // set the scheduled spike time to zero to keep
                // the neuron from spiking when an inhibitory input
                // spike has canceled a previously scheduled spike.
                schedSpikeTime = 0;
            }

        }
    }
    if (converged) {
        if (schedSpikeTime) {
            eventTp = RESPIKE;
        }
        else {
            eventTp = SPIKE;
        }

        // Make sure we wait at least one cycle before spiking
        if (calcTime == inTime) {
            calcTime += simStepSize;
        }
        schedSpikeTime = calcTime;

        #ifdef DEBUG_NEURON
        cout << "Scheduling spike at " << calcTime << endl;
        #endif

        SNet->ScheduleNEvent( eventTp, calcTime, this );
    }
}

void BasicNeuron::SetMaxScaledWeight()
{
    unsigned int i;
    float lastPtnl = 0.0;
    float maxPtnl = 0.0;

    for (i=0; i<pspLSize; i++) {
        maxPtnl = pspLookup[i];

        if (maxPtnl < lastPtnl) {
            maxPtnl = lastPtnl;
            break;
        }
        else {
            lastPtnl = maxPtnl;
        }
    }

    if (maxPtnl == 0.0)
        return;

    // maxScaledWeight is increased a little bit beyond its
    // true value to make sure that a weight of 1.0 will always
    // exceed the threshold at the peak of the psp curve.
    maxScaledWeight = ( thresholdPtnl / maxPtnl ) + 0.000001;
}

float* BasicNeuron::InitializeLookupTable(int index)
{
    unsigned int i;

    switch (index) {
        case 0:
            pspLookup = new float[pspLSize];
            for (i=0; i<pspLSize; i++) {
                pspLookup[i] = 0.0;
            }
            if ( FillLookupTables(index) ) {
                return pspLookup;
            }
            else {
                cerr << "FillLookupTables() failed in Neuron.\n";
                return 0;
            }
            break;
        case 1:
            dPspLookup = new float[pspLSize];
            for (i=0; i<pspLSize; i++) {
                dPspLookup[i] = 0.0;
            }
            if ( FillLookupTables(index) ) {
                return dPspLookup;
            }
            else {
                cerr << "FillLookupTables() failed in Neuron.\n";
                return 0;
            }
            break;
        default:
            return 0;
    }
}

float* BasicNeuron::GetTableParams(int index, int& numParams)
{
    float* paramList;
    switch (index) {
        case 0:
            numParams = 1;
            paramList = new float[1];
            paramList[0] = synTimeConst;
            break;
        case 1:
            numParams = 1;
            paramList = new float[1];
            paramList[0] = synTimeConst;
            break;
        default:
            return 0;
    }
    return paramList;
}

int BasicNeuron::SetLookupTables(FunctionLookup* funcRef)
{
    pspLookup = funcRef->GetLookupTable(this, 0, pspLSize, pspStepSize);
    dPspLookup = funcRef->GetLookupTable(this, 1, pspLSize, pspStepSize);

    // test for table pointers
    if ( pspLookup == 0 || dPspLookup == 0 ) {
        return 0;
    }

    SetMaxScaledWeight();
    return 1;
}

int BasicNeuron::FillLookupTables(int index)
{

    if (pspLSize == 0 || pspStepSize == 0) {
        return 0;
    }

    float calcTm = 0.0;
    unsigned int i;

    switch (index) {
        case 0:
            //cout << "Filling PSP Table..." << endl;
            for (i=0; i<pspLSize; i++) {
                calcTm = (float(i) * pspStepSize);
                PspKernel(calcTm, &pspLookup[i]);
                /*if ( (modff( float(i) / 10, &temp )) == 0.0 ) {
                    cout << "pspLookup[" << i << "] = " << pspLookup[i] << endl;
                }*/
            }
            break;
        case 1:
            //cout << "Filling dPSP Table..." << endl;
            for (i=0; i<pspLSize; i++) {
                calcTm = (float(i) * pspStepSize);
                DPspKernel(calcTm, &dPspLookup[i]);
                /*if ( (modff( float(i) / 10, &temp )) == 0.0 ) {
                    cout << "dPspLookup[" << i << "] = " << dPspLookup[i] << endl;
                }*/
            }
            break;
        default:
            return 0;
    }

    usePspLookup = true;
    return 1;
}

inline void BasicNeuron::PspKernel(float calcTime, float* pspElement)
{
    // Model the PSP as Ep = ( t/( Ts^2 ) ) * exp[ -t / Ts ]
    float term1 = 0.0, term2 = 0.0, expTerm = 0.0, convCalcTime;

    // convert time from microseconds to miliseconds
    convCalcTime = calcTime / 1000.0;

    term1 = (convCalcTime/(pow(synTimeConst, 2)));
    expTerm = -convCalcTime/synTimeConst;
    term2 = exp(expTerm);
    *pspElement = term1 * term2;
    //*pspElement = (calcTime/(pow(memTimeConst, 2)))*exp(-calcTime/memTimeConst);
}

inline void BasicNeuron::DPspKernel(float calcTime, float* dPspElement)
{
    float tempVar, convCalcTime;

    // convert time from microseconds to miliseconds
    convCalcTime = calcTime / 1000.0;

    tempVar = exp(-convCalcTime/synTimeConst)/(pow(synTimeConst, 2));
    tempVar = tempVar - (convCalcTime/(pow(synTimeConst, 3))*exp(-convCalcTime/synTimeConst));

    // Divide the result by 1000 so that it has units of mV/us (converted from
    // mV/ms).
    *dPspElement = tempVar / 1000.0;
}