alphaneuron.cpp

此代码经过大量使用

/***************************************************************************
                          alphaneuron.cpp  -  description
                             -------------------
    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 <math.h>

#include "alphaneuron.h"
#include "euler.h"
#include "types.h"
#include "network.h"
#include "functionlookup.h"

// Uncomment for debugging messages
//#define DEBUG_ALPHANEURON

AlphaNeuron::AlphaNeuron(AmIdInt neuronId):
    Neuron(neuronId)
{
    epspLookup = 0;
    edPspLookup = 0;
    ipspLookup = 0;
    idPspLookup = 0;
    refPeriod = 20;
    eTblFilled = false;
    iTblFilled = false;
}

AlphaNeuron::~AlphaNeuron()
{
}

void AlphaNeuron::InputSpike(SynapseItr& inSynapse, AmTimeInt inTime, unsigned int numSyn = 0)
{
/*******************************************************************
* MODIFY FOR AlphaNeuron!!!!!
*******************************************************************/

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

    #ifdef DEBUG_ALPHANEURON
    cout << "\nNeuron " << nId << " receiving a spike from " << inNeuron << endl;
    #endif

    // don't do anything if still in refractory period
    if ( spikeTime && (inTime - spikeTime) <= refPeriod ) {
        return;
    }
    NEvent eventTp = NOACTION;
    iterate = 1;
    converged = 0;
    calcTime = 0;
    funcTime = 0;
    histSize = 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;
        for (unsigned int i=0; i<numSyn; ++i) {
            tmpInput.weight += inSynapse[i]->GetWeight();
        }
    }
    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);
        }
    }


    #ifdef DEBUG_ALPHANEURON
    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 no change in calcTime.
    *            (Converges)
    *        b) The derivative of the function becomes negative.
    *            (Does not converge)
    **************************************************************************/
    #ifdef DEBUG_ALPHANEURON
    cout << "Starting main loop...\n";
    #endif

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

        RoundTime(calcTime);

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

        for (i=histBeginIdx; i<histSize; i++) {
            tmpInput = inputHist[i];
            funcTime = calcTime - tmpInput.time;
            funcWeight = tmpInput.weight;

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

            tblIndex = (funcTime / pspStepSize);
            if (tblIndex < pspLSize) {
                if ( funcWeight > 0.0 ) {
                    currState = currState + (funcWeight * (epspLookup[tblIndex]));
                    currDeriv = currDeriv + (funcWeight * (edPspLookup[tblIndex]));
                }
                else {
                    currState = currState + (funcWeight * (ipspLookup[tblIndex]));
                    currDeriv = currDeriv + (funcWeight * (idPspLookup[tblIndex]));
                }
            }
            else {
                ++histBeginIdx;
            }
        }

        #ifdef DEBUG_ALPHANEURON
        cout << "currState: " << currState << endl;
        cout << "currDerive: " << 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_ALPHANEURON
            cout << "stateDelta: " << stateDelta << endl;
            cout << "threshCrs: " << threshCrs << endl;
            cout << "lstThreshCrs: " << lstThreshCrs << endl;
            #endif
            if ((threshCrs - lstThreshCrs) < convergeRes) {
                converged = 1;
                iterate = 0;
                //schedSpikeTime = int(threshCrs);
                #ifdef DEBUG_ALPHANEURON
                cout << "threshCrs - lstThreshCrs < convergeRes\n";
                cout << "Converged\n";
                #endif
            }
            else {
                converged = 0;
                if (threshCrs > maxThreshCrs) {
                    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_ALPHANEURON
        cout << "Scheduling spike at " << calcTime << endl;
        #endif
        
        SNet->ScheduleNEvent( eventTp, calcTime, this );
    }
}

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

    for (i=0; i<pspLSize; i++) {
        maxPtnl = epspLookup[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;

    // Scale the inhibitory psp to match the amplitude of the epsp.
    // This will ensure consistent behavior for all time constants
    // (the maximum of the ipsp gets lower with increasing synaptic
    // time constants).
    float imaxPtnl = 0.0;
    lastPtnl = 0.0;
    for (i=0; i<pspLSize; i++) {
        imaxPtnl = ipspLookup[i];

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

    float scaleFactor = 0.9*maxPtnl/imaxPtnl;
    for (i=0; i<pspLSize; i++) {
        ipspLookup[i] = scaleFactor*ipspLookup[i];
        idPspLookup[i] = scaleFactor*idPspLookup[i];
    }
}

int AlphaNeuron::SetLookupTables(FunctionLookup* funcRef)
{
    epspLookup = funcRef->GetLookupTable(this, 0, pspLSize, pspStepSize);
    edPspLookup = funcRef->GetLookupTable(this, 1, pspLSize, pspStepSize);
    ipspLookup = funcRef->GetLookupTable(this, 2, pspLSize, pspStepSize);
    idPspLookup = funcRef->GetLookupTable(this, 3, pspLSize, pspStepSize);

    // test for table pointers
    if ( epspLookup == 0 || edPspLookup == 0 || ipspLookup == 0 || idPspLookup == 0 ) {
        return 0;
    }

    SetMaxScaledWeight();
    return 1;
}

float* AlphaNeuron::InitializeLookupTable(int index)
{
    unsigned int res = pspStepSize, i;
    float* returnTbl;

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

    // initialize arrays, if needed
    if (!eTblFilled) {
        epspLookup = new float[pspLSize];
        edPspLookup = new float[pspLSize];
        for (i=0;i<pspLSize;i++) {
            epspLookup[i] = 0.0;
            edPspLookup[i] = 0.0;
        }
    }

    if (!iTblFilled) {
        ipspLookup = new float[pspLSize];
        idPspLookup = new float[pspLSize];
        for (i=0;i<pspLSize;i++) {
            ipspLookup[i] = 0.0;
            idPspLookup[i] = 0.0;
        }
    }

    Euler* eul = new Euler(1, 0, pspLSize * 100, synTimeConst, memTimeConst);

    switch (index) {
        case 0:
            if (!eTblFilled)
                eul->GenExcitatoryTable(res, epspLookup, edPspLookup);
            delete eul;
            eTblFilled = true;
            returnTbl = epspLookup;
            break;
        case 1:
            if (!eTblFilled)
                eul->GenExcitatoryTable(res, epspLookup, edPspLookup);
            delete eul;
            eTblFilled = true;
            returnTbl = edPspLookup;
            break;
        case 2:
            if (!iTblFilled)
                eul->GenInhibitoryTable(res, ipspLookup, idPspLookup);
            delete eul;
            iTblFilled = true;
            returnTbl = ipspLookup;
            break;
        case 3:
            if (!iTblFilled)
                eul->GenInhibitoryTable(res, ipspLookup, idPspLookup);
            delete eul;
            iTblFilled = true;
            returnTbl = idPspLookup;
            break;
        default:
            delete eul;
            return 0;
    }
    return returnTbl;
}

float* AlphaNeuron::GetTableParams(int index, int& numParams)
{
    // this will be the same for all tables, so ignore the index
    numParams = 2;
    float* paramList = new float[2];
    paramList[0] = synTimeConst;
    paramList[1] = memTimeConst;
    return paramList;
}