euler.cpp

此代码经过大量使用

/***************************************************************************
                          euler.cpp  -  description
                             -------------------
    copyright            : (C) 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.                                   *
 *                                                                         *
 ***************************************************************************/

#include <iostream.h>
#include <math.h>
#include "euler.h"

Euler::Euler(unsigned int step,
            unsigned int min,
            unsigned int max,
            float synapticTimeConst,
            float membraneTimeConst):
    stepSize(step),
    plotMin(min),
    plotMax(max)
{
    if (synapticTimeConst == 0) {
        tauS = 2000.0;
    }
    else {
        tauS = synapticTimeConst * 1000;
    }

    if (membraneTimeConst == 0) {
        tauM = 10000.0;
    }
    else {
        tauM = membraneTimeConst * 1000;
    }

    if (plotMax <= plotMin) {
        cerr << "Plot min and max invalid.\n";
        return;
    }

    plotSize = (plotMax - plotMin)/stepSize;
    plotData = new float[plotSize];
    plotDeriv = new float[plotSize];
    	
    unsigned int i = 0;
    for (i = 0; i < plotSize; i++) {
        plotData[i] = 0.0;
        plotDeriv[i] = 0.0;
    }
}

Euler::~Euler()
{
}


/*void Euler::GenPlot()
{
    unsigned int i;

    lastV = 0.0;	// initial condition
    lastSlope = 0.0;
    time = 0;
    // Get the initial slope
    //AlphaExcite();
    AlphaInhibit();
    lastSlope = currSlope;

    for (i=0; i< plotSize; i++) {
        //AlphaExcite();
        AlphaInhibit();
        plotData[i] = currV;

        // display value for every 500 iterations
        if (drem(i, 1000) == 0.0) {
            cout << currV << endl;
        }
        lastV = currV;
        lastSlope = currSlope;
        time = time + stepSize;
    }
}*/

void Euler::GenExcitatoryTable(int resolution, float* tbl, float* derivTbl)
{
    unsigned int i;
    int tblStep = 0, counter = 0;

    lastV = 0.0;	// initial condition
    lastSlope = 0.0;
    time = 0;

    tblStep = resolution / stepSize;

    // Make sure the table is empty
    for (i=0; i < plotSize; i++) {
        plotData[i] = 0.0;
        plotDeriv[i] = 0.0;
    }

    // Get the initial slope
    AlphaExcite();
    lastSlope = currSlope;

    for (i=0; i< plotSize; i++) {
        AlphaExcite();
        plotData[i] = currV;
        plotDeriv[i] = currSlope;

        lastV = currV;
        lastSlope = currSlope;
        time = time + stepSize;
    }

    // Copy the tables
    for (i=0; i < plotSize; i+=tblStep) {
        tbl[counter] = plotData[i];
        //cout << "plotData[" << i << "] = " << plotData[i] << endl;
        derivTbl[counter] = plotDeriv[i];
        counter++;
    }
}

void Euler::GenInhibitoryTable(int resolution, float* tbl, float* derivTbl)
{
    unsigned int i;
    int tblStep = 0, counter = 0;

    lastV = 0.0;	// initial condition
    lastSlope = 0.0;
    time = 0;

    tblStep = resolution / stepSize;

    // Make sure the table is empty
    for (i=0; i < plotSize; i++) {
        plotData[i] = 0.0;
        plotDeriv[i] = 0.0;
    }

    // Get the initial slope
    AlphaInhibit();
    lastSlope = currSlope;

    for (i=0; i< plotSize; i++) {
        AlphaInhibit();
        plotData[i] = currV;
        plotDeriv[i] = currSlope;

        lastV = currV;
        lastSlope = currSlope;
        time = time + stepSize;
    }

    // Copy the tables
    for (i=0; i < plotSize; i+=tblStep) {
        tbl[counter] = plotData[i];
        derivTbl[counter] = plotDeriv[i];
        counter++;
    }
}

void Euler::AlphaExcite()
{
	float dV = 0.0, I = 0.0;
	
	// find new V value
	currV = lastV + float(stepSize)*(lastSlope);

	// calculate new derivative
	I = ( exp( -(float(time) + float(stepSize) )/tauS)/tauS );
	dV = I - ( currV/tauM );
	currSlope = dV;
}

void Euler::AlphaInhibit()
{
	float dV = 0.0, I = 0.0;

	// Increase the value of the synaptic time constant
	// for inhibitory currents.  This is being done
	// as a quick hack to make AlphaNeuron more useful
	// for synchronizing populations (the ipsp should
	// be slower than the epsp -- see Mus Silicium by
	// Hopfield and Brody).  Ideally, the neuron model
	// should allow time constants to be set for both
	// i and e psps, but it does not at this time.  When
	// that feature is added, this should be removed to
	// allow users more control over AlphaNeuron behavior.
	// TODO: Remove this when inhibitory synaptic TCs can
	// be set separately from excitatory TCs.
	static float itauS = tauS+4000.0;

	// find new V value
	currV = lastV + float(stepSize)*(lastSlope);

	// calculate new derivative
	I = ( (float(time) + float(stepSize)) * exp( -(float(time) +
	        float(stepSize) )/itauS)/(itauS * itauS) );

	dV = I - ( currV/tauM );
	currSlope = dV;
}