ObjSim/src/layer.cpp

#include "sys.hpp"   // for libcwd
#include "debug.hpp" // for libcwd

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// #include <QThread>

#include "layer.hpp"
#include "simloop.hpp"
#include "libcsim.hpp"

#include "spiketrain.hpp"

#include <limits>  // for detecting denormals
#include "rng_paraqueue_pool.hpp"


////////////////////////////////

layer::layer(int n)
        :SimElement(seLayer), N(n), NoiseSigma(0),
        NormPos(false), PoissonLambda(12.5), NoiseAmplitude(0),
        BalancedInhibition(0), TooManySpikes(false), mSpikeTrain(0)
{
    Name="Layer";

    int i;

    SetRandomSpikeFrequency(1.0); //Hz
    Binary = true;

    SpikeFileName = "ospikes";
    if (Binary) FileType=".bin";
    else FileType=".dat";

    RandomInputStrength=20.;

    Dmax=DMAX;
    if (MacroTimeStep <=Dmax)
    {
        Dout(dc::layer | dc::warning, "MacroTimeStep too low: " << MacroTimeStep);
    }

    Dout(dc::layer, "LayerConstructor " << "N=" << N << "TimeStep = " << dt << "ms");

    v = new  float [N];
    u = new  float [N];
    Input = new float [N];
    for (i=0;i<N;++i) Input[i] =0;
    N_firings_max=int(400*N/dt);   // upper limit on the number of fired neurons per sec

    NewArray2d(firings, N_firings_max,2);   //mybe change the dimensions?? what is more efficient? (faster?)
    // at least it would require less memory (only two pointers would have to be stored
    // indeces and timings of spikes
    // firings[N_firings_max][2]; // indeces and timings of spikes
    N_firings=1;       // spike timings
    firings[0][0]=-Dmax;   // put a dummy spike at -D for simulation efficiency
    firings[0][1]=0;   // index of the dummy spike

    mPositions.resize(N);
    Pos = &(mPositions[0]);
    NormPos=true;
    Nx=Ny=0;
}

layer::~layer()
{
    Dout(dc::layer, "Deleting Layer");
    delete[] v;
    delete[] u;
    DeleteArray2d(firings, N_firings_max);
    delete[] Input;
    if (mSpikeTrain) {
        delete mSpikeTrain;
    }
}

long layer::calcMemoryConsumption()
{
    long MemSum=0;
    MemSum += 3*N*sizeof(float); // float *v, *u;     float *Input;
    MemSum += N_firings_max*sizeof(int*);  //NewArray2d(firings, N_firings_max,2);
    MemSum += N_firings_max*2*sizeof(int); //NewArray2d(firings, N_firings_max,2);
    MemSum += N*sizeof(vector2d);  //Pos = new vector2d [N];
    return MemSum;
}


int layer::WriteSimInfo(fstream &fw)
{
    fw << "<" << seTypeString << " id=\"" << IdNumber <<  "\" Type=\"" << seType << "\" Name=\"" << Name << "\"> \n";
    fw << "<LayerDimensions N=\"" << N << "\" Nx=\"" << Nx << "\" Ny=\"" << Ny << "\"/> \n";
//      fw << "<N value=\"" << N << "\"/> \n";
//      fw << "<Nx value=\"" << Nx << "\"/> \n";
//      fw << "<Ny value=\"" << Ny << "\"/> \n";
    fw << "<NormPos value=\"" << NormPos << "\"/> \n";
    fw << "<Dmax value=\"" << Dmax << "\"/> \n";
    fw << "<Resetable value=\"" << resetable << "\"/> \n";
    fw << "</" << seTypeString << "> \n";
}

int layer::WriteSimInfo(fstream &fw, const string &ChildInfo)
{
    stringstream sstr;
    sstr << "<LayerDimensions N=\"" << N << "\" Nx=\"" << Nx << "\" Ny=\"" << Ny << "\"/> \n";
    sstr << "<NormPos value=\"" << NormPos << "\"/> \n";
    sstr << "<Dmax value=\"" << Dmax << "\"/> \n";
    sstr << "<Resetable value=\"" << resetable << "\"/> \n";
    sstr << ChildInfo;
    SimElement::WriteSimInfo(fw, sstr.str());
}

const vector<vector2d>& layer::getPositions()
{
    return mPositions;
}


int layer::SetupPositions()
{
    int i;
    // calculate neuron positions
    // default: quadratic (allmost)
    Ny = int(sqrt(float(N)));
    Nx = int(ceil(float(N)/float(Ny)));
    for (i=0;i<N;++i) Pos[i].set(i%Nx, int(i/Nx));
    NormPos=false;
}

int layer::SetupPositions(int _nx, int _ny, bool Normalize)
{
    int i;
    Nx=_nx;
    Ny=_ny;
    if ((Nx*Ny) != N) {
        Dout(dc::layer | dc::warning, "WARNING: Nx*Ny =" << Nx*Ny
             << "!= N=" << N << "  Nx=" << Nx << "  Ny=" << Ny << "  _ny=" << _ny);
    }
    if (Normalize)  {
        for (i=0; i<N; ++i) {
            Pos[i].set(float(i%Nx)/Nx, float(int(i/Nx))/Ny);
        }
    } else {
        for (i=0; i<N; ++i) {
            Pos[i].set(i%Nx, int(i/Nx));
        }
    }
    NormPos = Normalize;
}

int layer::SetupPositionsShift(int _nx, int _ny, float xshift, float yshift, bool Normalize)
{
    int i;
    Nx=_nx;
    Ny=_ny;
    if ((Nx*Ny) != N) Dout(dc::layer | dc::warning, "Nx*Ny =" << Nx*Ny << "!= N=" << N);
    if (Normalize)  {
        for (i=0;i<N;++i) Pos[i].set((float(i%Nx)/Nx)+xshift, (float(int(i/Nx))/Ny)+yshift);
    } else {
        for (i=0;i<N;++i) Pos[i].set(i%Nx+xshift, int(i/Nx)+yshift);
    }
    NormPos = Normalize;
}

int  layer::SetupPositionsLinDiscrete(int NeuronsPerPatch)
{
    int i;
    int NPatches = N/NeuronsPerPatch;
    for (i=0;i<N;++i) {
        Pos[i].set(float(i*NPatches/N)/NPatches,0);
    }
    NormPos=true;

    Ny = int(sqrt(float(N)));
    Nx = int(ceil(float(N)/float(Ny)));
}


void layer::hallo()
{
    Dout(dc::layer, "Hallo");
}

float layer::GetDt()
{
    return dt;
}

int layer::proceede(int t)
{
    int i,j,k;
    for (int i=1;i<=nris;++i)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input
    if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;      // random thalamic input

}

// prepare next sec
int layer::prepare(int sec)
{
    SimElement::prepare(sec);
    int i,j,k;
    FILE   *fs;
//     Dout(dc::layer,  Name << " firing rate=" << float(N_firings-1)/(dt*N) << "");fflush(stdout);

    Dout(dc::layer, Name << " firing rate=" << 1000*float(N_firings-1)/((MacroTimeStep+Dmax)*dt*N) << " Hz");

    // save firings

    stringstream sstr;
    sstr << sec;

    std::string tmpstring = DataDirectory+SpikeFileName+SimTag+FileType;
    const char* CurFileName = tmpstring.c_str();

    // save spike trains
    if (sec ==0) // open new file
    {
        fs = fopen(CurFileName,"w");
        if (fs) {
            if (Binary) {
                for (i=1;i<N_firings;i++)
                    if (firings[i][0] >=0)
                        fwrite(firings[i], 2*sizeof(firings[0][0]), 1, fs);

            } else {
                for (i=1;i<N_firings;i++)
                    if (firings[i][0] >=0)
                        fprintf(fs, "%d  %d\n", firings[i][0], firings[i][1]);
            }
            fclose(fs);
        } else {
            Dout(dc::layer | dc::fatal, "failed opening file " << CurFileName << " for writing");
        }

    }
    else
    { // append spikes to existing file
        fs = fopen(CurFileName,"a");
        if (fs) {
            if (Binary) {
                int **temp_firings;
                int *buffer = new int [2*N_firings];
                int count=0;
                for (i=1;i<N_firings;++i) {
                    if (firings[i][0] >= 0) {
                        buffer[2*count] = firings[i][0] + sec*MacroTimeStep;
                        buffer[2*count+1] = firings[i][1];
                        ++count;
                    }
                }
                fwrite(buffer, 2*count*sizeof(*buffer), 1, fs);
                delete[] buffer;
            } else {
                for (i=1;i<N_firings;i++)
                    if (firings[i][0] >=0)
                        fprintf(fs, "%d  %d\n", dt*firings[i][0]+sec*MacroTimeStep, firings[i][1]);
            }
            fclose(fs);
        } else {
            Dout(dc::layer | dc::fatal, "failed opening file " << CurFileName << " for writing");
        }
    }

    k=N_firings-1;
    while (MacroTimeStep-firings[k][0]<Dmax) k--;
    for (i=1;i<N_firings-k;i++)
    {
        firings[i][0]=firings[k+i][0]-MacroTimeStep;
        firings[i][1]=firings[k+i][1];
    }
    N_firings = N_firings-k;
}

int layer::reset(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    SimElement::reset(t);
    Dout(dc::layer,  Name << " RESET LAYER ");
    cerr << Name << " ERROR in layer::reset: usually the layer class has to overwrite this method\n";
    fflush(stdout);

    // delete all spikes which are still in the delay line
    // firings with firing times higher than t-MaximumDelay

    int MaximumDelay = MainSimLoop->GetMaximumDelay();
    int i,j,k;

    k=N_firings-1;
    while (t-firings[k][0]<MaximumDelay) k--;
    N_firings = k+1;

    for (int i=0;i<N;++i) {
        u[i]=0;
        v[i]=0;
        Input[i]=0;
    }

}


void layer::SetName(const char* _name)
{
    SimElement::SetName(_name);
    SpikeFileName = Name + "spikes.dat";
    MemPotFileName = Name + "MemPot.dat";
}

void layer::SetSpikeFileName(char* FileName)
{
    SpikeFileName = FileName;
    cout << "SpikeFile: " << SpikeFileName;
}

void layer::SetRandomInputStrength(float ris)
{
    RandomInputStrength = ris;
    Dout(dc::layer,  "set RandomInputStrength to " << RandomInputStrength << "");
}

void layer::SetRandomSpikeFrequency(float _rsf)
{
    RandomSpikeFrequency = _rsf;
    NumRandomInputSpikes = dt*N*RandomSpikeFrequency/1000;
    nris = (int) floor(NumRandomInputSpikes);
    RestProb = int(1000*(NumRandomInputSpikes - nris));
    Dout(dc::layer,   "  NumRandomInputSpikes= " << NumRandomInputSpikes << "  RestProb= " << RestProb << "  nris= " << nris << "");

}

void layer::SetNoiseSigma(float _sigma)
{
    NoiseSigma=_sigma;
}

void layer::SetBalancedInhibition(float value)
{
    BalancedInhibition = value;
}

void layer::SetNoiseAmplitude(float value)
{
    NoiseAmplitude = value;
}


float* layer::GetInputPointer(csimInputChannel InputNumber)
{
    return Input;
}

float* layer::GetPspPointer(csimInputChannel channel)
{
    return Input;
}


int layer::SaveSimInfo()
{
    int i;
    string FileName=Name+string(".info");
    cout << "Save SimInfo in File: " << FileName;
    FILE *fw;
    fw = fopen((DataDirectory+FileName).c_str(),"w");
    fprintf(fw, "<LayerDimensions Nx=\"%d\" Ny=\"%d\" />\n", Nx, Ny);
    for (i=0;i<N;++i)
    {
        fprintf(fw, "%d", i);
        Pos[i].fprintfcs(fw);
    }
//      fwrite(&ns, sizeof(ns), 1, fw);
//      fwrite(&nt, sizeof(nt), 1, fw);
//      fwrite(&M, sizeof(M), 1, fw);
//      fwrite(&Dmax, sizeof(Dmax), 1, fw);
//      fwrite(&maxN_pre, sizeof(maxN_pre), 1, fw);

// //     NewArray2d(delays_length,N,Dmax); // distribution of delays
// //     NewArray3d(delays,N,Dmax,M);  // arrangement of delays
//      for (i=0;i<ns;++i) fwrite(post[i], M*sizeof(post[0][0]), 1, fw);
//      for (i=0;i<ns;++i) fwrite(s[i], M*sizeof(s[0][0]), 1, fw);
//      for (i=0;i<ns;++i) fwrite(delays_length[i], Dmax*sizeof(delays_length[0][0]),1,fw);
//      for (i=0;i<ns;++i) for (j=0;j<Dmax;++j) fwrite(delays[i][j], M*sizeof(delays[0][0][0]),1,fw);


    fclose(fw);
    Dout(dc::layer,  "  saved SimInfo ");


}

int layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<NumObserve; ++i) Buffer[i] = &(v[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, NumObserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

void layer::ThrowTooManySpikes(int time)
{
    cout << Name << "Two many spikes at t=" << time << " (ignoring all)";
    N_firings=1;  // not zero to keep the dummy spike as first element in the array
    TooManySpikes=true;
}

int layer::TuneNoiseAmplitude(float DesiredFrequency)
{
    int time=0;
    NoiseAmplitude = 0.01;
    int SettlingTime = int(100./dt); // ms/dt
    int MeasureTime = int(500./dt); // ms/dt
    float Tollerance = 0.01;
    float ChangeFactor=1;
    const int MaxTuneTrials=100;
    int TuneTrials=0;
    bool Stop=false;
    int MStep=0;
    enum ChangeDirVal {up, down};
    ChangeDirVal ChDir=up;
    ChangeDirVal LastChDir=up;
    while (! Stop) {
        for (int i=0;i<SettlingTime;++i) {
            proceede(time++);
            if (( time % MacroTimeStep) == 0) {
                prepare(MStep++);
            }
        }
        int LastNSpikesStart=N_firings;
        int SpikeCounter=0;
        TooManySpikes=false;
        for (int i=0;i<MeasureTime;++i) {
            proceede(time++);
            if (( time % MacroTimeStep) == 0) {
                SpikeCounter += N_firings-LastNSpikesStart;
                prepare(MStep++);
                LastNSpikesStart=N_firings;
            }
        }
        SpikeCounter += N_firings-LastNSpikesStart;
        float SpikeFrequency = 1000*float(SpikeCounter)/(dt*MeasureTime*N);
        cout << "SC=" << SpikeCounter << " NoiseAmp=" << NoiseAmplitude
             << "  CurSpikeFrequency=" << SpikeFrequency << "Hz\n";
        if (abs(SpikeFrequency-DesiredFrequency) < Tollerance*DesiredFrequency) {
            Stop = true;
            Dout(dc::layer,  "successfully finished noise tuning");
        } else {
            if ((SpikeFrequency>DesiredFrequency) || TooManySpikes) {
                ChDir=down;
                if (LastChDir == up) {
                    ChangeFactor *= 0.5;
                }
                NoiseAmplitude /= (1+ChangeFactor);
            } else {
                ChDir=up;
                if (LastChDir == down) {
                    ChangeFactor *= 0.5;
                }
                NoiseAmplitude *= (1+ChangeFactor);
            }
            LastChDir = ChDir;
            if (TuneTrials++ > MaxTuneTrials) {
                Stop=true;
                Dout(dc::layer,  "ERROR: unsuccessfully finished noise tuning");
            }
        }
    }

    prepare(0);
    N_firings=1;
}


int layer::TuneBalancedInhibition(float DesiredFrequency)
{
    int time=0;
    BalancedInhibition = 0.1;
    int SettlingTime = int(100./dt); // ms/dt
    int MeasureTime = int(500./dt); // ms/dt
    float Tollerance = 0.01;
    float ChangeFactor=1;
    const int MaxTuneTrials=100;
    int TuneTrials=0;
    bool Stop=false;
    int MStep=0;
    enum ChangeDirVal {up, down};
    ChangeDirVal ChDir=up;
    ChangeDirVal LastChDir=up;
    while (! Stop) {
        for (int i=0;i<SettlingTime;++i) {
            proceede(time++);
            if (( time % MacroTimeStep) == 0) {
                prepare(MStep++);
            }
        }
        int LastNSpikesStart=N_firings;
        int SpikeCounter=0;
        TooManySpikes=false;
        for (int i=0;i<MeasureTime;++i) {
            proceede(time++);
            if (( time % MacroTimeStep) == 0) {
                SpikeCounter += N_firings-LastNSpikesStart;
                prepare(MStep++);
                LastNSpikesStart=N_firings;
            }
        }
        SpikeCounter += N_firings-LastNSpikesStart;
        float SpikeFrequency = 1000*float(SpikeCounter)/(dt*MeasureTime*N);
        cout << "SC=" << SpikeCounter << " BalancedInhibition=" << BalancedInhibition
             << "  CurSpikeFrequency=" << SpikeFrequency << "Hz\n";
        if (abs(SpikeFrequency-DesiredFrequency) < Tollerance*DesiredFrequency) {
            Stop = true;
            Dout(dc::layer,  "successfully finished BalancedInhibition tuning");
        } else {
            if ((SpikeFrequency<DesiredFrequency) || TooManySpikes) {
                ChDir=down;
                if (LastChDir == up) {
                    ChangeFactor *= 0.5;
                }
                BalancedInhibition /= (1+ChangeFactor);
            } else {
                ChDir=up;
                if (LastChDir == down) {
                    ChangeFactor *= 0.5;
                }
                BalancedInhibition *= (1+ChangeFactor);
            }
            LastChDir = ChDir;
            if (TuneTrials++ > MaxTuneTrials) {
                Stop=true;
                Dout(dc::layer,  "ERROR: unsuccessfully finished BalancedInhibition tuning");
            }
        }
    }

    prepare(0);
    N_firings=1;
}


SpikeTrain* layer::GetSpikeTrain()
{
    Dout(dc::layer,  "layer::GetSpikeTrain()");
    if (!mSpikeTrain) {
        Dout(dc::layer,  "created new SpikeTrain object");
        return mSpikeTrain=new SpikeTrain(firings, &N_firings, &N, dt);
    } else {
        Dout(dc::layer,  "returned old SpikeTrain object");
        return mSpikeTrain;
    }
}

/////////////////////////////
// Konstruktor der Basisklasse muss in der Initialisierungsliste aufgerufen werden
liflayer::liflayer(
    int n,  float tau,
    float thresh, float k0, float k1,
    float k2, float kffi, float kfbi,
    float taus, float _InputSat, float _NoiseSigma)
        : layer(n), TauM(tau),
        Threshold(thresh), K0(k0), K1(k1), K2(k2),
        Kffi(kffi), Kfbi(kfbi), TauS(taus), InputSaturation(_InputSat)
{
    NoiseSigma=_NoiseSigma;
    int i;
    cout << "Leaky Integrate- and Fire Neuron: K1="
         << K1 << " K2=" << K2 << " Kffi=" << Kffi << " Kfbi= " << Kfbi << "\n";
    RandomInputStrength=5;
    TauMfac = dt/TauM;
    TauSfac = exp(-dt/TauS);
    LastSpike = new int[N];
    for (i=0;i<N;++i) {
        LastSpike[i]=-2;
        v[i]=0;
        u[i]=0;
    }

    Input1 = new float[N];
    for (i=0;i<N;++i) Input1[i]=0;

    InhibitionDelay = int(1./dt);  // 1 milli second
    Dout(dc::layer,  "InhibitionDelay=" << InhibitionDelay << " TimeSteps ");
    mquer = new float [MacroTimeStep + InhibitionDelay + 1];
    squer = new float [MacroTimeStep + InhibitionDelay + 1];
    for (i=0; i<(MacroTimeStep+InhibitionDelay+1); ++i) {
        mquer[i]=0;
        squer[i]=0;
    }

    RunAvgFac = exp(-dt/2.0);
    Dout(dc::layer,  "RunAvgFac=" << RunAvgFac << "");

    RefractoryPeriod = int(2/dt);
}


liflayer::~liflayer()
{
    delete[] LastSpike;
    delete[] Input1;
    delete[] mquer;
    delete[] squer;
}


int liflayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    float CurInput;

    // Izhikevich-like noise generation (poissonian noise with strong pulses)
//      for (int i=1;i<=nris;++i)  Input[getrandom(N)]= RandomInputStrength;        // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input


    last_N_firings = N_firings;

    Iinh = K0+Kffi*squer[t]+Kfbi*mquer[t];

    for (i=0;i<N;i++)
    {
        CurInput = K1*Input[i] + K2*Input1[i] + gsl_ran_gaussian(gslr, NoiseSigma);
        u[i] += InputSaturation*CurInput/(CurInput+Iinh);
        if ((i==0) && (rec !=0)) rec->record(dt*TotalTime, squer[t], mquer[t]);  // record potentials for DEBUGing
        squer[t+InhibitionDelay] += dt*Input[i];
        v[i] +=TauMfac*(-v[i] + u[i]);

        if ((v[i] > Threshold) && (t-LastSpike[i] > 2/dt))    // did it fire?
        {
            v[i] -= Threshold;       // voltage reset
            LastSpike[i] = t;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            mquer[t+InhibitionDelay] += 1;
            // Dout(dc::layer,  "Mquer increased to " <<  mquer[t+InhibitionDelay] << "");

            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }

        u[i] *= TauSfac;
    }
    mquer[t+InhibitionDelay+1] = mquer[t+InhibitionDelay]*RunAvgFac;
    squer[t+InhibitionDelay+1] = squer[t+InhibitionDelay]*RunAvgFac;
//    Dout(dc::layer,  "Squer= " << squer[t] << "  Mquer= " << mquer[t] << "");

    for (i=0;i<N;i++)
    {
        Input[i] = 0.0;   // reset the input
        Input1[i] = 0.0;
    }
}


int liflayer::prepare(int sec)
{
    int i,j;
    // prepare for the next MacroTimeStep
    layer::prepare(sec);
    for (i=0;i<N;++i) LastSpike[i] -= MacroTimeStep;
    for (i=0;i<InhibitionDelay+1;i++)
    {
        mquer[i] = mquer[MacroTimeStep+i];
        squer[i] = squer[MacroTimeStep+i];
    }
    for (i=InhibitionDelay+1;i<MacroTimeStep+InhibitionDelay+1;++i)
    {
        mquer[i]=0;
        squer[i]=0;
    }
}

float* liflayer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA:
        Dout(dc::layer,  "EG/DG Input");
        return Input;
        break;
    case csimInputChannel_GABAa:
        Dout(dc::layer,  "recurrent input");
        return Input1;
        break;
    default:
        return Input;
        break;
    }
    return 0;
}

int liflayer::SetKfbi(float _kfbi)
{
    Kfbi=_kfbi;
}

int liflayer::InitInhibitionPot(float _mquer, float _squer)
{
    int i;
    for (i=0;i<InhibitionDelay+1;i++)
    {
        mquer[i] = _mquer;
        squer[i] = _squer;
    }
}

//////////////////////////////////

thetaliflayer::thetaliflayer(
    int n, float tau,
    float thresh, float k0, float k1,
    float k2, float kffi, float kfbi,
    float taus, float _InputSat, float _NoiseSigma, int _ThetaPeriod)

        : liflayer(n, tau, thresh, k0, k1, k2, kffi, kfbi,taus,  _InputSat, _NoiseSigma), ThetaPeriod (_ThetaPeriod)
{
    ThetaPhase = int(ThetaPeriod/dt);
    ThetaOn=true;
}

thetaliflayer::~thetaliflayer()
{

}


int thetaliflayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    float CurInput;

    // Izhikevich-like noise generation (poissonian noise with strong pulses)
//      for (int i=1;i<=nris;++i)  Input[getrandom(N)]= RandomInputStrength;        // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input


    last_N_firings = N_firings;

    Iinh = K0+Kffi*squer[t]+Kfbi*mquer[t];
    float ThetaK2;
    if (ThetaOn)
    {
        if (ThetaPhase++ >= (ThetaPeriod/dt)) ThetaPhase=0;
        ThetaK2 = K2*(dt*ThetaPhase/ThetaPeriod);
    } else ThetaK2 = K2;
//     cout << ThetaK2/K2 << " ";
    for (i=0;i<N;i++)
    {
        CurInput = K1*Input[i] + ThetaK2*Input1[i] + gsl_ran_gaussian(gslr, NoiseSigma);
        u[i] += InputSaturation*CurInput/(CurInput+Iinh);
        if ((i==0) && (rec !=0)) rec->record(dt*TotalTime, squer[t], mquer[t]);  // record potentials for DEBUGing
        squer[t+InhibitionDelay] += dt*Input[i];
        v[i] +=TauMfac*(-v[i] + u[i]);

        if ((v[i] > Threshold) && (t-LastSpike[i] > 2/dt))    // did it fire?
        {
            v[i] -= Threshold;       // voltage reset
            LastSpike[i] = t;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            mquer[t+InhibitionDelay] += 1;
            // Dout(dc::layer,  "Mquer increased to " <<  mquer[t+InhibitionDelay] << "");

            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
        mquer[t+InhibitionDelay+1] = mquer[t+InhibitionDelay]*RunAvgFac;
        squer[t+InhibitionDelay+1] = squer[t+InhibitionDelay]*RunAvgFac;
//    Dout(dc::layer,  "Squer= " << squer[t] << "  Mquer= " << mquer[t] << "");

        u[i] *= TauSfac;
    }


    for (i=0;i<N;i++)
    {
        Input[i] = 0.0;   // reset the input
        Input1[i] = 0.0;
    }
}


int thetaliflayer::prepare(int sec)
{
    int i,j;
    // prepare for the next MacroTimeStep
    layer::prepare(sec);
    for (i=0;i<N;++i) LastSpike[i] -= MacroTimeStep;
    for (i=0;i<InhibitionDelay+1;i++)
    {
        mquer[i] = mquer[MacroTimeStep+i];
        squer[i] = squer[MacroTimeStep+i];
    }
    for (i=InhibitionDelay+1;i<MacroTimeStep+InhibitionDelay+1;++i)
    {
        mquer[i]=0;
        squer[i]=0;
    }
}

void thetaliflayer::SetThetaOn(bool status)
{
    ThetaOn = status;
}
////////////////////////

int IzhParas::print()
{
    Dout(dc::layer,  "IzhParas: a=" << A << " b=" << B << " c=" << C << " d=" << D << " e=" << E << " f=" << F << "");
}


/////////////////////////////
izhlayer::izhlayer (int n, float a, float b, float c, float d, float e, float f, float _InputSat, float _InhInputSat): layer(n), A(a), B(b), C(c), D(d), E(e), F(f), InputSaturation(_InputSat), InhInputSaturation(_InhInputSat)
{
    Dout(dc::layer,  "Constructor: Quadratic Integrate- and Fire Neuron ");
    Init();
}

izhlayer::izhlayer (int n, IzhParas _paras, float _InputSat, float _InhInputSat): layer(n), A(_paras.A), B(_paras.B), C(_paras.C), D(_paras.D), E(_paras.E), F(_paras.F), InputSaturation(_InputSat), InhInputSaturation(_InhInputSat)
{
    SetNoiseSigma(_paras.sigma);
    Dout(dc::layer,  "Constructor: Quadratic Integrate- and Fire Neuron ");
    _paras.print();
    Init();
}


int izhlayer::Init()
{
    int i;
    InitEquilibriumPotentials();
    TauS = 2.0;
    TauSfac = exp(-dt/TauS);
    TauInh = 10.0;
    TauInhFac = exp(-dt/TauInh);
    ExSynPot = new float [N];
    for (i=0;i<N;++i) ExSynPot[i]=0;
    InhSynPot = new float [N];
    for (i=0;i<N;++i) InhSynPot[i]=0;
    Dout(dc::layer,  "dt = " << dt << "");
}

int izhlayer::SetTauInh(float _TauInh)
{
    TauInh = _TauInh;
    TauInhFac = exp(-dt/TauInh);
}

int izhlayer::SetTauEx(float _TauEx)
{
    TauS = _TauEx;
    TauSfac = exp(-dt/TauS);
}


izhlayer::izhlayer (IzhType type, int n,  float _InputSat, float _InhInputSat): layer(n), A(0.02), B(0.2), C(-65), D(2), E(140), F(5), InputSaturation(_InputSat), InhInputSaturation(_InhInputSat)
{
    switch (type) {
    case Integrator: {
        Dout(dc::layer,  "Initialize Izhikevich-Neuron Layer as Integrator");
        A=0.02;
        B=-0.1;
        C=-55;
        D=6;
        E=108;
        F=4.1;
    } break;
    case FastSpiking: {
        Dout(dc::layer,  "Initialize Izhikevich-Neuron Layer as FastSpiking");
        A=0.02;
        B=0.25;
        C=-65;
        D=1;
        E=140;
        F=5;
    } break;
    default:
        Dout(dc::layer,  "No correct Izhikevich-Neuron-Type");
        break;
    }

    Init();
}


izhlayer::~izhlayer()
{
    delete[] ExSynPot;
    delete[] InhSynPot;
}

int izhlayer::InitEquilibriumPotentials()
{
    int i;
    float p = (F-B)/0.04;
    float q = E/0.04;
    float FixPoint = -60;
    float FP2 = -40;
    cout << "p=" << p << " q=" << q;
    if (q<=(p*p/4.0))
    {
        FixPoint = -p/2.0 - sqrt(p*p/4 - q);
        FP2        = -p/2.0 + sqrt(p*p/4 - q);
        Dout(dc::layer,  "  Equilibrium point is v=" << FixPoint << "mV");
    } else Dout(dc::layer,  "NO Equilibrium point!! v=" << FixPoint << "mV");

    float FPu = B*FixPoint;
    float StartJitterV = 0.5; //(FP2 - FixPoint)/40;
//     float StartJitterV = 0.0; //(FP2 - FixPoint)/40;
    EquilibriumPot=FixPoint;

    for (i=0;i<N;i++)  v[i]= FixPoint + StartJitterV*(float(getrandom(1000))/1000 - 0.5);      // initial values for v  // calculate equilibrium point!!
    for (i=0;i<N;i++)  u[i]=FPu ;    // initial values for u
}


float* izhlayer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhSynPot;
    } break;
    default:
        return Input;
    }
    return 0;
}

int izhlayer::WriteSimInfo(fstream &fw)
{
    fw << "<" << seTypeString << " id=\"" << IdNumber <<  "\" Type=\"" << seType << "\" Name=\"" << Name << "\"> \n";
    fw << "<LayerType value=\"" << "Izhikevich" << "\"/> \n";
    fw << "<LayerDimensions N=\"" << N << "\" Nx=\"" << Nx << "\" Ny=\"" << Ny << "\"/> \n";
    fw << "<NormPos value=\"" << NormPos << "\"/> \n";
    fw << "<Dmax value=\"" << Dmax << "\"/> \n";
    fw << "<InputSaturation ex=\"" << InputSaturation << " \" inh=\"" << InhInputSaturation << "\"/> \n";
    fw << "<Parameter A=\"" << A << "\" B=\"" << B << "\" C=\"" << C << "\" D=\"" << D << "\" E=\"" << E << "\" F=\"" << F << "\"" << "/>\n";
    fw << "<EquilibriumPot value=\"" << EquilibriumPot << "\"/> \n";
    fw << "</" << seTypeString << "> \n";
}


int izhlayer::WriteSimInfo(fstream &fw, const string &ChildInfo)
{
    stringstream sstr;
//     sstr << "<LayerType value=\"" << "Izhikevich" << "\"/> \n";
    sstr << "<InputSaturation ex=\"" << InputSaturation << " \" inh=\"" << InhInputSaturation << "\"/> \n";
    sstr << "<Parameter A=\"" << A << "\" B=\"" << B << "\" C=\"" << C << "\" D=\"" << D << "\" E=\"" << E << "\" F=\"" << F << "\"" << "/>\n";
    sstr << "<EquilibriumPot value=\"" << EquilibriumPot << "\"/> \n";
    sstr << ChildInfo;
    layer::WriteSimInfo(fw, sstr.str());
}



int izhlayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    for (int i=1;i<=nris;++i)        Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
    if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;      // random thalamic input

    last_N_firings = N_firings;
    for (i=0;i<N;i++)
    {
        ExSynPot[i] += InputSaturation*Input[i]/(Input[i] + 1 + InhSynPot[i]);
        if ((i==0) && (rec != 0)) rec->record(dt*TotalTime, v[i], Input[i]);
        v[i] += dt*((0.04 *v[i] + F) * v[i] + E - u[i] + ExSynPot[i]);

        //    v[i] +=dt*0.5*((0.04*v[i]+F)*v[i]+E-u[i]+Input[i]); // for numerical stability // time step is 0.5 ms
        u[i] += dt * A * (B * v[i] - u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        ExSynPot[i] *= TauSfac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
    }
}


int izhlayer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 3*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(ExSynPot[i+NNumber]);
//      for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(Input[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(InhSynPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}




////////////////////////////////////////

izh2layer::izh2layer (int n, float a, float b, float c, float d, float e, float f, float _InputSat): izhlayer(n,a,b,c,d,e,f,_InputSat)
{
}

izh2layer::izh2layer (int n, IzhParas _paras, float _InputSat): izhlayer(n,_paras, _InputSat)
{
    SetNoiseSigma(_paras.sigma);
}

izh2layer::~izh2layer()
{
}

int izh2layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
//      for (int i=1;i<=nris;++i)     Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input

    last_N_firings = N_firings;
    int RecNum=52;
    float tmp=0;
//     if (rec != 0) rec->record(dt*TotalTime, u[RecNum], Input[RecNum]);
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
//    Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
//    Input[i] += gsl_ran_gaussian(gslr, NoiseSigma);
        ExSynPot[i] += gsl_ran_gaussian(gslr, NoiseSigma) + InputSaturation*Input[i]/(Input[i] + 1 + InhSynPot[i]);
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ExSynPot[i]-InhSynPot[i]);
//    v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ExSynPot[i]-InhSynPot[i] + gsl_ran_gaussian(gslr, NoiseSigma) + );
        u[i] +=dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        ExSynPot[i] *= TauSfac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
    }
}

////////////////////////////////////////

izh3layer::izh3layer(int n, float a, float b, float c, float d, float e, float f, float _InputSat, float IRPotDiff, float _InhInputSat): izhlayer(n,a,b,c,d,e,f,_InputSat, _InhInputSat)
{
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV
    InhInput = new float [N];
    for (int i=0;i<N;++i) InhInput[i] =0;

}

izh3layer::izh3layer (int n, IzhParas _paras, float _InputSat, float IRPotDiff, float _InhInputSat): izhlayer(n,_paras, _InputSat,_InhInputSat)
{
    SetNoiseSigma(_paras.sigma);
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV
    InhInput = new float [N];
    for (int i=0;i<N;++i) InhInput[i] =0;
}

izh3layer::~izh3layer()
{
    delete[] InhInput;
}

float* izh3layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    default:
        return Input;
    }
    return 0;
}


int izh3layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh3layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
//      for (int i=1;i<=nris;++i)     Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input

    last_N_firings = N_firings;
    int RecNum=52;
    float tmp=0;
//     if (rec != 0) rec->record(dt*TotalTime, u[RecNum], Input[RecNum]);
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
//    Input[i] += gsl_ran_gaussian(gslr, NoiseSigma);
        ExSynPot[i] += InputSaturation*Input[i]/(Input[i] + 1);
        InhSynPot[i] += InhInputSaturation*InhInput[i]/(InhInput[i] + 1);
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ExSynPot[i]-InhSynPot[i]*(v[i]-InhReversePot));
//    v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ExSynPot[i]-InhSynPot[i] + gsl_ran_gaussian(gslr, NoiseSigma) + );
        u[i] +=dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        InhInput[i] = 0;  // reset input
        ExSynPot[i] *= TauSfac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
    }
//     Dout(dc::layer,  "izh3layer::proceede___feddisch");fflush(stdout);
}


//////////////////////////////////////////

izh4layer::izh4layer (int n, IzhParas _paras, float _InputSat, float IRPotDiff, float _TauLink): izh3layer(n,_paras, _InputSat, IRPotDiff), TauLink(_TauLink)
{
    LinkSynPot = new float[N];
    for (int i=0;i<N;++i) LinkSynPot[i]=0;
    TauLinkFac = exp(-dt/TauLink);
}

float* izh4layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhSynPot;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Linking Input");
        return LinkSynPot;
    } break;
    default:
        return Input;
    }
    return 0;
}


izh4layer::~izh4layer()
{
    delete[] LinkSynPot;
}

int izh4layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh4layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        Input[i] *= 1+LinkSynPot[i];
        ExSynPot[i] += InputSaturation*Input[i]/(Input[i] + 1);
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ExSynPot[i]-InhSynPot[i]*(v[i]-InhReversePot));
        u[i] +=dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        ExSynPot[i] *= TauSfac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
        LinkSynPot[i] *= TauLinkFac;
    }
//     Dout(dc::layer,  "izh4layer::proceede___feddisch");fflush(stdout);
}


//////////////////////////////////////////

////////////////////////////////////////

izh5layer::izh5layer(int n, float a, float b, float c, float d, float e, float f, float _InputSat, float IRPotDiff, float _InhInputSat): izhlayer(n,a,b,c,d,e,f,_InputSat, _InhInputSat)
{
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV
}

izh5layer::izh5layer (int n, IzhParas _paras, float _InputSat, float IRPotDiff, float _InhInputSat): izhlayer(n,_paras, _InputSat,_InhInputSat)
{
    SetNoiseSigma(_paras.sigma);
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV
}

izh5layer::~izh5layer()
{
}

float* izh5layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return ExSynPot;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhSynPot;
    } break;
    default:
        return ExSynPot;
    }
    return 0;
}

int izh5layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh5layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
//      for (int i=1;i<=nris;++i)     Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input

    last_N_firings = N_firings;
    int RecNum=52;
    float tmp=0;
//     if (rec != 0) rec->record(dt*TotalTime, u[RecNum], Input[RecNum]);
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        ExSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]
                   + InputSaturation*ExSynPot[i]/(ExSynPot[i]+1)
                   - (v[i]-InhReversePot)*InhInputSaturation*InhSynPot[i]/(InhSynPot[i]+1));
        u[i] += dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }


    }
    for (i=0;i<N;i++)
    {
        ExSynPot[i] *= TauSfac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
    }
//     Dout(dc::layer,  "izh5layer::proceede___feddisch");fflush(stdout);
}


//////////////////////////////////////////

izh6layer::izh6layer (int n, IzhParas _paras, float _InputSat, float IRPotDiff, float _InhInputSat, float _NmdaSaturation): izhlayer(n,_paras, _InputSat,_InhInputSat), NmdaSaturation(_NmdaSaturation)
{
    SetNoiseSigma(_paras.sigma);
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV
    NmdaSynPot = new float [N];
    for (int i=0;i<N;++i) NmdaSynPot[i]=0;
    SetTauNmda(100);
    SetNmdaAmpaRatio(0.1);
}

izh6layer::~izh6layer()
{
    delete[] NmdaSynPot;
}

float* izh6layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return ExSynPot;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhSynPot;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return Input;
    }  break;
    default:
        return ExSynPot;
    }
    return 0;
}


int izh6layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh6layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
//      for (int i=1;i<=nris;++i)     Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input

    last_N_firings = N_firings;
    int RecNum=52;
    float tmp=0;
//     if (rec != 0) rec->record(dt*TotalTime, u[RecNum], Input[RecNum]);
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        ExSynPot[i] += Input[i];
        NmdaSynPot[i] += NmdaAmpaRatio*Input[i];
        ExSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]
                   + InputSaturation*ExSynPot[i]/(ExSynPot[i]+1)
                   + NmdaSaturation*INmda(v[i])*NmdaSynPot[i]/(NmdaSynPot[i]+1)
                   - (v[i]-InhReversePot)*InhInputSaturation*InhSynPot[i]/(InhSynPot[i]+1));
        u[i] += dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }


    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        ExSynPot[i] *= TauSfac;  // decay the input potential
        NmdaSynPot[i] *= TauNmdaFac;  // decay the input potential
        InhSynPot[i] *= TauInhFac; // decay inhibitory input potential
    }
//     Dout(dc::layer,  "izh6layer::proceede___feddisch");fflush(stdout);
}

void izh6layer::SetNmdaAmpaRatio(float ratio)
{
    NmdaAmpaRatio = ratio;
}

void izh6layer::SetTauNmda(float value)
{
    TauNmda=value;
    TauNmdaFac = exp(-dt/TauNmda);
}

int izh6layer::WriteSimInfo(fstream &fw)
{
    fw << "<" << seTypeString << " id=\"" << IdNumber <<  "\" Type=\"" << seType << "\" Name=\"" << Name << "\"> \n";
    fw << "<LayerType value=\"" << "Izhikevich6" << "\"/> \n";
    fw << "<LayerDimensions N=\"" << N << "\" Nx=\"" << Nx << "\" Ny=\"" << Ny << "\"/> \n";
    fw << "<NormPos value=\"" << NormPos << "\"/> \n";
    fw << "<Dmax value=\"" << Dmax << "\"/> \n";
    fw << "<InputSaturation ex=\"" << InputSaturation << " \" inh=\"" << InhInputSaturation << " \" nmda=\"" << NmdaSaturation << "\"/> \n";
    fw << "<EquilibriumPot value=\"" << EquilibriumPot << "\"/> \n";
    fw << "</" << seTypeString << "> \n";
}

int izh6layer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 4*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(ExSynPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(NmdaSynPot[i+NNumber]);
//      for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(Input[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(InhSynPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

//////////////////////////////////////////

izh7layer::izh7layer (int n, IzhParas _paras, float _InputSat, float IRPotDiff, float _InhInputSat, float _NmdaSaturation)
        :izhlayer(n,_paras, _InputSat,_InhInputSat), NmdaSaturation(_NmdaSaturation),
        TauAmpaRise(0.5), TauAmpaFall(2.4),
        TauGabaRise(1.0), TauGabaFall(7.0),
        TauNmdaRise(5.5), TauNmdaFall(100)
{
    SetNoiseSigma(_paras.sigma);      // 0.00233 -> 3Hz f�r IzhIntegrator
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV

    Dout(dc::layer,  "izh7layer::izh7layer");



    // use arrays from izhlayer
    AmpaFallPot = ExSynPot;
    GabaFallPot = InhSynPot;

    NmdaRisePot = new float [N];
    NmdaFallPot = new float [N];
    AmpaRisePot = new float [N];
    GabaRisePot = new float [N];
    InhInput = new float [N];
    NmdaAmpaInput = new float [N];

    // initialize arrays with zero
    for (int i=0;i<N;++i) {
        NmdaRisePot[i]=0;
        NmdaFallPot[i]=0;
        AmpaRisePot[i]=0;
        GabaRisePot[i]=0;
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;
    }

    SetTauNmda(TauNmdaRise, TauNmdaFall);
    SetTauGaba(TauGabaRise, TauGabaFall);
    SetTauAmpa(TauAmpaRise, TauAmpaFall);
    SetNmdaAmpaRatio(0.1);
}

izh7layer::~izh7layer()
{
    delete[] NmdaRisePot;
    delete[] NmdaFallPot;
    delete[] AmpaRisePot;
    delete[] GabaRisePot;
    delete[] InhInput;
    delete[] NmdaAmpaInput;
}

float* izh7layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return NmdaAmpaInput;
    }  break;
    default:
        return Input;
    }
    return 0;
}


int izh7layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh7layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
//      for (int i=1;i<=nris;++i)     Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
//      if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input

    last_N_firings = N_firings;
    int RecNum=52;
    float tmp=0;
//     if (rec != 0) rec->record(dt*TotalTime, u[RecNum], Input[RecNum]);
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
//    Input[i] += NoiseSigma*gsl_ran_poisson(gslr, PoissonLambda);
        Input[i] += NmdaAmpaInput[i];

        Input[i] *= AmpaCorr;
        InhInput[i] *= GabaCorr;
        NmdaAmpaInput[i] *= NmdaCorr;

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaRisePot[i] += InhInput[i];
        GabaFallPot[i] += InhInput[i];

//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        float Ampa = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float Gaba = GabaFallPot[i] - GabaRisePot[i];

        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]
                   + InputSaturation*Ampa/(Ampa+1)
                   + NmdaSaturation*INmda(v[i])*Nmda/(Nmda+1)
                   - (v[i]-InhReversePot)*InhInputSaturation*Gaba/(Gaba+1));
        u[i] += dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }


    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaRisePot[i] *= TauGabaRiseFac;
        GabaFallPot[i] *= TauGabaFallFac;
    }
//     Dout(dc::layer,  "izh7layer::proceede___feddisch");fflush(stdout);
}

void izh7layer::SetNmdaAmpaRatio(float ratio)
{
    NmdaAmpaRatio = ratio;
}

void izh7layer::SetTauNmda(float rise, float fall)
{
    TauNmdaRise = rise;
    TauNmdaFall = fall;
    TauNmdaRiseFac = exp(-dt/TauNmdaRise);
    TauNmdaFallFac = exp(-dt/TauNmdaFall);
    NmdaCorr = DoubleExpCorrectionFactor(TauNmdaRise, TauNmdaFall);
}

void izh7layer::SetTauAmpa(float rise, float fall)
{
    TauAmpaRise = rise;
    TauAmpaFall = fall;
    TauAmpaRiseFac = exp(-dt/TauAmpaRise);
    TauAmpaFallFac = exp(-dt/TauAmpaFall);
    AmpaCorr = DoubleExpCorrectionFactor(TauAmpaRise, TauAmpaFall);
}

void izh7layer::SetTauGaba(float rise, float fall)
{
    TauGabaRise = rise;
    TauGabaFall = fall;
    TauGabaRiseFac = exp(-dt/TauGabaRise);
    TauGabaFallFac = exp(-dt/TauGabaFall);
    GabaCorr = DoubleExpCorrectionFactor(TauGabaRise, TauGabaFall);
}

int izh7layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "Izhikevich7" << "\"/> \n";
    sstr << "<InputSaturation ex=\"" << InputSaturation << " \" inh=\"" << InhInputSaturation << " \" nmda=\"" << NmdaSaturation << "\"/> \n";
    sstr << "<EquilibriumPot value=\"" << EquilibriumPot << "\"/> \n";
    layer::WriteSimInfo(fw, sstr.str());
}

int izh7layer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 7*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(AmpaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(AmpaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(GabaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+4*nobserve] = &(GabaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+5*nobserve] = &(NmdaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+6*nobserve] = &(NmdaFallPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

/////////////////////////////////////////


izh8layer::izh8layer (int n, IzhParas _paras, float IRPotDiff)
        :izhlayer(n,_paras),
        TauAmpaRise(0.5), TauAmpaFall(2.4),
        TauGabaRise(1.0), TauGabaFall(7.0),
        TauNmdaRise(5.5), TauNmdaFall(100)
{
    SetNoiseSigma(_paras.sigma);      // 0.00233 -> 3Hz f�r IzhIntegrator
    InhReversePot = EquilibriumPot-IRPotDiff;  // default: Ruhepot - 10 mV

    Dout(dc::layer,  "izh8layer::izh8layer");



    // use arrays from izhlayer
    AmpaPot = ExSynPot;
    GabaFallPot = InhSynPot;

    NmdaRisePot = new float [N];
    NmdaFallPot = new float [N];
    AmpaRisePot = new float [N];
    AmpaFallPot = new float [N];
    GabaRisePot = new float [N];
    InhInput = new float [N];
    NmdaAmpaInput = new float [N];

    // initialize arrays with zero
    for (int i=0;i<N;++i) {
        NmdaRisePot[i]=0;
        NmdaFallPot[i]=0;
        AmpaRisePot[i]=0;
        AmpaFallPot[i]=0;
        GabaRisePot[i]=0;
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;
    }

    SetTauNmda(TauNmdaRise, TauNmdaFall);
    SetTauGaba(TauGabaRise, TauGabaFall);
    SetTauAmpa(TauAmpaRise, TauAmpaFall);
    SetNmdaAmpaRatio(0.1);
}

izh8layer::~izh8layer()
{
    delete[] NmdaRisePot;
    delete[] NmdaFallPot;
    delete[] AmpaRisePot;
    delete[] AmpaFallPot;
    delete[] GabaRisePot;
    delete[] InhInput;
    delete[] NmdaAmpaInput;
}

float* izh8layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return NmdaAmpaInput;
    }  break;
    default:
        return Input;
    }
    return 0;
}

float* izh8layer::GetPspPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case  csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return AmpaPot;
    }  break;
    default:
        return AmpaPot;
    }
    return 0;
}


int izh8layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh8layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;

    if (BinRec) BinRec->record();

    for (i=0;i<N;i++)
    {
//    Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        Input[i] += NoiseAmplitude*gsl_ran_poisson(gslr, PoissonLambda);
        Input[i] += NmdaAmpaInput[i];

        Input[i] *= AmpaCorr;
        InhInput[i] *= GabaCorr;
        NmdaAmpaInput[i] *= NmdaCorr;

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaRisePot[i] += InhInput[i];
        GabaFallPot[i] += InhInput[i];

//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        AmpaPot[i] = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float Gaba = GabaFallPot[i] - GabaRisePot[i];

        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]
                   + AmpaPot[i]
                   + INmda(v[i])*Nmda
                   - (v[i]-InhReversePot)*(Gaba+BalancedInhibition));
        u[i] += dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }


    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaRisePot[i] *= TauGabaRiseFac;
        GabaFallPot[i] *= TauGabaFallFac;
    }
//     Dout(dc::layer,  "izh8layer::proceede___feddisch");fflush(stdout);
}

void izh8layer::SetNmdaAmpaRatio(float ratio)
{
    NmdaAmpaRatio = ratio;
}

void izh8layer::SetTauNmda(float rise, float fall)
{
    TauNmdaRise = rise;
    TauNmdaFall = fall;
    TauNmdaRiseFac = exp(-dt/TauNmdaRise);
    TauNmdaFallFac = exp(-dt/TauNmdaFall);
    NmdaCorr = DoubleExpCorrectionFactor(TauNmdaRise, TauNmdaFall);
}

void izh8layer::SetTauAmpa(float rise, float fall)
{
    TauAmpaRise = rise;
    TauAmpaFall = fall;
    TauAmpaRiseFac = exp(-dt/TauAmpaRise);
    TauAmpaFallFac = exp(-dt/TauAmpaFall);
    AmpaCorr = DoubleExpCorrectionFactor(TauAmpaRise, TauAmpaFall);
}

void izh8layer::SetTauGaba(float rise, float fall)
{
    TauGabaRise = rise;
    TauGabaFall = fall;
    TauGabaRiseFac = exp(-dt/TauGabaRise);
    TauGabaFallFac = exp(-dt/TauGabaFall);
    GabaCorr = DoubleExpCorrectionFactor(TauGabaRise, TauGabaFall);
}

int izh8layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "Izhikevich8" << "\"/> \n";
    sstr << "<NoiseAmplitude value=\"" << NoiseAmplitude << "\"/> \n";
    sstr << "<BalancedInhibition value=\"" << BalancedInhibition << "\"/> \n";
    sstr << "<TauAmpa "
    << "Rise=\"" << TauAmpaRise << "\" "
    << "Fall=\"" << TauAmpaFall << "\" "
    << "/> \n";
    sstr << "<TauGaba "
    << "Rise=\"" << TauGabaRise << "\" "
    << "Fall=\"" << TauGabaFall << "\" "
    << "/> \n";
    sstr << "<TauNmda "
    << "Rise=\"" << TauNmdaRise << "\" "
    << "Fall=\"" << TauNmdaFall << "\" "
    << "/> \n";
    izhlayer::WriteSimInfo(fw, sstr.str());
}

int izh8layer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 7*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(AmpaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(AmpaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(GabaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+4*nobserve] = &(GabaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+5*nobserve] = &(NmdaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+6*nobserve] = &(NmdaFallPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

/////////////////////////////////////////


izh9layer::izh9layer(int n, IzhParas _paras, float GabaARPotDiff, float GabaBRPotDiff)
        :izhlayer(n,_paras),
        TauAmpaRise(0.5), TauAmpaFall(2.4),
        TauGabaARise(1.0), TauGabaAFall(7.0),
        TauGabaBRise(13.5), TauGabaBFall(140),
        TauNmdaRise(5.5), TauNmdaFall(100),
        GabaABRatio(0.1)
{
    SetNoiseSigma(_paras.sigma);
    InhReversePot = EquilibriumPot-GabaARPotDiff;  // GabaA
    GabaBReversePot = EquilibriumPot-GabaBRPotDiff; // GabaB

    Dout(dc::layer,  "izh9layer::izh9layer");



    // use arrays from izhlayer
    AmpaPot = ExSynPot;
    GabaAFallPot = InhSynPot;

    NmdaRisePot = new float [N];
    NmdaFallPot = new float [N];
    GabaBRisePot = new float [N];
    GabaBFallPot = new float [N];
    AmpaRisePot = new float [N];
    AmpaFallPot = new float [N];
    GabaARisePot = new float [N];
    InhInput = new float [N];
    NmdaAmpaInput = new float [N];

    // initialize arrays with zero
    for (int i=0;i<N;++i) {
        NmdaRisePot[i]=0;
        NmdaFallPot[i]=0;
        GabaBRisePot[i]=0;
        GabaBFallPot[i]=0;
        AmpaRisePot[i]=0;
        AmpaFallPot[i]=0;
        GabaARisePot[i]=0;
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;
    }

    SetTauNmda(TauNmdaRise, TauNmdaFall);
    SetTauGabaA(TauGabaARise, TauGabaAFall);
    SetTauGabaB(TauGabaBRise, TauGabaBFall);
    SetTauAmpa(TauAmpaRise, TauAmpaFall);
    SetNmdaAmpaRatio(0.1);
}

izh9layer::~izh9layer()
{
    delete[] NmdaRisePot;
    delete[] NmdaFallPot;
    delete[] AmpaRisePot;
    delete[] AmpaFallPot;
    delete[] GabaBRisePot;
    delete[] GabaBFallPot;
    delete[] GabaARisePot;
    delete[] InhInput;
    delete[] NmdaAmpaInput;
}

float* izh9layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return NmdaAmpaInput;
    }  break;
    default:
        return Input;
    }
    return 0;
}

float* izh9layer::GetPspPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case  csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return AmpaPot;
    }  break;
    default:
        return AmpaPot;
    }
    return 0;
}


int izh9layer::proceede(int TotalTime)
{
//     Dout(dc::layer,  "izh9layer::proceede");fflush(stdout);
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;

    if (BinRec) BinRec->record();

    for (i=0;i<N;i++)
    {
//    Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        Input[i] += NoiseAmplitude*gsl_ran_poisson(gslr, PoissonLambda);
        Input[i] += NmdaAmpaInput[i];

        Input[i] *= AmpaCorr;  // to normalize the amplitude of the ampa epsp
        InhInput[i] *= GabaACorr;
        NmdaAmpaInput[i] *= NmdaCorr;

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaARisePot[i] += InhInput[i];
        GabaAFallPot[i] += InhInput[i];
        GabaBRisePot[i] +=GabaABRatio*InhInput[i];
        GabaBFallPot[i] += GabaABRatio*InhInput[i];

//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        AmpaPot[i] = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float GabaA = GabaAFallPot[i] - GabaARisePot[i];
        float GabaB = GabaBFallPot[i] - GabaBRisePot[i];

        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]
                   + AmpaPot[i]
                   + INmda(v[i])*Nmda
                   - (v[i]-InhReversePot)*(GabaA+BalancedInhibition))
               - (v[i]-GabaBReversePot)*GabaB;
        u[i] += dt*A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
//         Dout(dc::layer,  N_firings << "spike"); fflush(stdout);
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }


    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaARisePot[i] *= TauGabaARiseFac;
        GabaAFallPot[i] *= TauGabaAFallFac;
        GabaBRisePot[i] *= TauGabaBRiseFac;
        GabaBFallPot[i] *= TauGabaBFallFac;
    }
//     Dout(dc::layer,  "izh9layer::proceede___feddisch");fflush(stdout);
}

void izh9layer::SetNmdaAmpaRatio(float ratio)
{
    NmdaAmpaRatio = ratio;
}


void izh9layer::SetGabaABRatio(float ratio)
{
    NmdaAmpaRatio = ratio;
}

void izh9layer::SetTauNmda(float rise, float fall)
{
    TauNmdaRise = rise;
    TauNmdaFall = fall;
    TauNmdaRiseFac = exp(-dt/TauNmdaRise);
    TauNmdaFallFac = exp(-dt/TauNmdaFall);
    NmdaCorr = DoubleExpCorrectionFactor(TauNmdaRise, TauNmdaFall);
}

void izh9layer::SetTauAmpa(float rise, float fall)
{
    TauAmpaRise = rise;
    TauAmpaFall = fall;
    TauAmpaRiseFac = exp(-dt/TauAmpaRise);
    TauAmpaFallFac = exp(-dt/TauAmpaFall);
    AmpaCorr = DoubleExpCorrectionFactor(TauAmpaRise, TauAmpaFall);
}

void izh9layer::SetTauGabaA(float rise, float fall)
{
    TauGabaARise = rise;
    TauGabaAFall = fall;
    TauGabaARiseFac = exp(-dt/TauGabaARise);
    TauGabaAFallFac = exp(-dt/TauGabaAFall);
    GabaACorr = DoubleExpCorrectionFactor(TauGabaARise, TauGabaAFall);
}

void izh9layer::SetTauGabaB(float rise, float fall)
{
    TauGabaBRise = rise;
    TauGabaBFall = fall;
    TauGabaBRiseFac = exp(-dt/TauGabaBRise);
    TauGabaBFallFac = exp(-dt/TauGabaBFall);
    GabaBCorr = DoubleExpCorrectionFactor(TauGabaBRise, TauGabaBFall);
}


int izh9layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "Izhikevich9" << "\"/> \n";
    sstr << "<NoiseAmplitude value=\"" << NoiseAmplitude << "\"/> \n";
    sstr << "<BalancedInhibition value=\"" << BalancedInhibition << "\"/> \n";
    sstr << "<Ampa "
    << "TauRise=\"" << TauAmpaRise << "\" "
    << "TauFall=\"" << TauAmpaFall << "\" "
    << "/> \n";
    sstr << "<GabaA "
    << "TauRise=\"" << TauGabaARise << "\" "
    << "TauFall=\"" << TauGabaAFall << "\" "
    << "ReversePot=\"" << InhReversePot << "\" "
    << "/> \n";
    sstr << "<GabaB "
    << "TauRise=\"" << TauGabaBRise << "\" "
    << "TauFall=\"" << TauGabaBFall << "\" "
    << "ReversePot=\"" << GabaBReversePot << "\" "
    << "/> \n";
    sstr << "<Nmda "
    << "TauRise=\"" << TauNmdaRise << "\" "
    << "TauFall=\"" << TauNmdaFall << "\" "
    << "/> \n";
    izhlayer::WriteSimInfo(fw, sstr.str());
}

int izh9layer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 9*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(AmpaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(AmpaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(GabaARisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+4*nobserve] = &(GabaAFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+5*nobserve] = &(NmdaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+6*nobserve] = &(NmdaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+7*nobserve] = &(GabaBRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+8*nobserve] = &(GabaBFallPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}



/////////////////////////////////////////
izhshuntlayer::izhshuntlayer(
    int n, IzhParas _izparas, LevyShuntingParas _shuntparas,
    float _InputSat, float _NoiseSigma)

        :izhlayer(n, _izparas, _InputSat),
        NoiseSigma(_NoiseSigma), inh(_shuntparas)
{
    int i;

    Dout(dc::layer,  "IzhShuntLayer Constructor");
    Dout(dc::layer,  "InputSaturation = " << InputSaturation << "");
    Input1 = new float[N];
    for (i=0;i<N;++i) Input1[i]=0;

    InhibitionDelay = int(1./dt);  // 1 milli second
    Dout(dc::layer,  "InhibitionDelay=" << InhibitionDelay << " TimeSteps ");
    mquer = new float [MacroTimeStep + InhibitionDelay + 1];
    squer = new float [MacroTimeStep + InhibitionDelay + 1];
    for (i=0; i<(MacroTimeStep+InhibitionDelay+1); ++i) {
        mquer[i]=0;
        squer[i]=0;
    }
    RunAvgFac = exp(-dt/2.0);
    Dout(dc::layer,  "RunAvgFac=" << RunAvgFac << "");
}


izhshuntlayer::~izhshuntlayer()
{
    delete[] Input1;
    delete[] mquer;
    delete[] squer;
}

int izhshuntlayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    float CurInput;
    float tmpinp;
    for (int i=1;i<=nris;++i)        Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
    if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;      // random thalamic input


    last_N_firings = N_firings;

    Iinh = inh.K0+inh.Kffi*squer[t]+inh.Kfbi*mquer[t];


    for (i=0;i<N;i++)
    {
//    ExSynPot[i] += InputSaturation*Input[i]/(Input[i] + 1 + InhSynPot[i]);
        CurInput = inh.K1*Input[i] + inh.K2*Input1[i] + abs(gsl_ran_gaussian(gslr, NoiseSigma));
        tmpinp = InputSaturation*CurInput/(CurInput + Iinh);
        if ((i==10) && (rec != 0)) rec->record(dt*TotalTime, v[i], Iinh);
        squer[t+InhibitionDelay] += dt*Input[i];   // running average of external input
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ tmpinp);
        u[i] +=A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            mquer[t+InhibitionDelay] += 1;              // running average of layer activity
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

    mquer[t+InhibitionDelay+1] = mquer[t+InhibitionDelay]*RunAvgFac;
    squer[t+InhibitionDelay+1] = squer[t+InhibitionDelay]*RunAvgFac;

    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        Input1[i]=0;
//    ExSynPot[i] *= TauSfac;  // decay the input potential
//    InhSynPot[i] *= TauSfac; // decay inhibitory input potential
    }

}

int izhshuntlayer::prepare(int step)
{
    int i,j;
    // prepare for the next MacroTimeStep
    layer::prepare(step);
    for (i=0;i<InhibitionDelay+1;i++)
    {
        mquer[i] = mquer[MacroTimeStep+i];
        squer[i] = squer[MacroTimeStep+i];
    }
    for (i=InhibitionDelay+1;i<MacroTimeStep+InhibitionDelay+1;++i)
    {
        mquer[i]=0;
        squer[i]=0;
    }
}

int izhshuntlayer::SetKfbi(float _kfbi)
{
    inh.Kfbi = _kfbi;
}


float* izhshuntlayer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA:
        Dout(dc::layer,  "EG/DG Input");
        return Input;
        break;
    case csimInputChannel_GABAa:
        Dout(dc::layer,  "recurrent input");
        return Input1;
        break;
    default:
        return Input;
        break;
    }
    return 0;
}
///////////////////////////////////////////


PresynIzhShuntLayer::PresynIzhShuntLayer(
    int n, IzhParas _izparas, LevyShuntingParas _shuntparas,
    float _InputSat, float _NoiseSigma)
        :izhshuntlayer(n, _izparas, _shuntparas, _InputSat, _NoiseSigma)
{
    PresynInhInput = new float[N];
    for (int i=0;i<N;++i) PresynInhInput[i]=0;
}

PresynIzhShuntLayer::~PresynIzhShuntLayer()
{
    delete[] PresynInhInput;
}

float* PresynIzhShuntLayer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA:
        Dout(dc::layer,  "EG/DG Input");
        return Input;
        break;
    case csimInputChannel_GABAa:
        Dout(dc::layer,  "recurrent input");
        return Input1;
        break;
    case csimInputChannel_NMDA_AMPA:
        Dout(dc::layer,  "Presynaptic Inhibition Input");
        return PresynInhInput;
    default:
        Dout(dc::layer,  "EG/DG Input");
        return Input;
        break;
    }
    return 0;
}

int PresynIzhShuntLayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    float CurInput;
    float tmpinp;
    for (int i=1;i<=nris;++i)        Input[getrandom(N)]= RandomInputStrength;     // random thalamic input
    if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;      // random thalamic input


    last_N_firings = N_firings;

    Iinh = inh.K0+inh.Kffi*squer[t]+inh.Kfbi*mquer[t];


    for (i=0;i<N;i++)
    {
//    ExSynPot[i] += InputSaturation*Input[i]/(Input[i] + 1 + InhSynPot[i]);
        CurInput = inh.K1*Input[i] + inh.K2*Input1[i]/(1+PresynInhInput[i]) + abs(gsl_ran_gaussian(gslr, NoiseSigma));
        tmpinp = InputSaturation*CurInput/(CurInput + Iinh);
        if ((i==10) && (rec != 0)) rec->record(dt*TotalTime, v[i], PresynInhInput[i]);
        squer[t+InhibitionDelay] += dt*Input[i];   // running average of external input
        v[i] +=dt*((0.04*v[i]+F)*v[i]+E-u[i]+ tmpinp);
        u[i] +=A*(B*v[i]-u[i]);

        if (v[i]>=30)                 // did it fire?
        {
            v[i] = C;                    // voltage reset
            u[i]+= D;                    // recovery variable reset
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            mquer[t+InhibitionDelay] += 1;              // running average of layer activity
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

    mquer[t+InhibitionDelay+1] = mquer[t+InhibitionDelay]*RunAvgFac;
    squer[t+InhibitionDelay+1] = squer[t+InhibitionDelay]*RunAvgFac;

    for (i=0;i<N;i++)
    {
        Input[i] = 0; // reset input
        Input1[i]=0;
        PresynInhInput[i] *= RunAvgFac;  // decay PresynapticInputPotential
//    ExSynPot[i] *= TauSfac;  // decay the input potential
//    InhSynPot[i] *= TauSfac; // decay inhibitory input potential
    }

}



//////////////////////////////////

lif2layer::lif2layer(int n, float tau, float thresh, float k0, float k1, float k2, float kffi, float kfbi, float taus, float _DesiredActivity, float _pylearnrate):
        liflayer(n, tau, thresh, k0, k1, k2,kffi, kfbi, taus)
{
    int i;
    DesiredActivity = N*_DesiredActivity;
    PyrToInhLearningRate = _pylearnrate/N;
    // initialize weights to inhibitory interneuron
    Dinh = new float[N];
    for (i=0;i<N;++i) Dinh[i]=1.0;

    InhibitoryActivity = new float [MacroTimeStep + InhibitionDelay + 1];
    for (i=0; i<(MacroTimeStep+InhibitionDelay+1); ++i) {
        InhibitoryActivity[i]=0;
    }

    learnpti=true;
    TauInhLearn=2.0;
    InhLearnFac = exp(-dt/TauInhLearn);

}

lif2layer::~lif2layer()
{
    delete[] Dinh;
    delete[] InhibitoryActivity;
}

int lif2layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;
    float CurInput;
    for (int i=1;i<=nris;++i)  Input[getrandom(N)]= RandomInputStrength;       // random thalamic input
    if (getrandom(1000)< RestProb)  Input[getrandom(N)]= RandomInputStrength;      // random thalamic input

    last_N_firings = N_firings;

    Iinh = K0+Kffi*squer[t]+Kfbi*InhibitoryActivity[t];

    for (i=0;i<N;i++)
    {
        CurInput = K1*Input[i] + K2*Input1[i];
        u[i] += CurInput/(CurInput+Iinh);
        if ((i==0) && (rec != 0)) rec->record(dt*TotalTime, mquer[t], Dinh[i]);
        squer[t+InhibitionDelay] += dt*Input[i];
        v[i] +=TauMfac*(-v[i] + u[i]);

        if ((v[i] > Threshold) && (t-LastSpike[i] > 2/dt))    // did it fire?
        {
            v[i] -= Threshold;       // voltage reset
            LastSpike[i] = t;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            mquer[t+InhibitionDelay] += 1;
            InhibitoryActivity[t+InhibitionDelay] += Dinh[i];
            if (learnpti) Dinh[i] += PyrToInhLearningRate*(mquer[t+InhibitionDelay-1]-DesiredActivity);
//         cout << PyrToInhLearningRate*(mquer[t+InhibitionDelay-1]-DesiredActivity) << " ";
            // Dout(dc::layer,  "Mquer increased to " <<  mquer[t+InhibitionDelay] << "");

            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
        mquer[t+InhibitionDelay+1] = mquer[t+InhibitionDelay]*RunAvgFac;
        squer[t+InhibitionDelay+1] = squer[t+InhibitionDelay]*RunAvgFac;
        InhibitoryActivity[t+InhibitionDelay+1] = InhibitoryActivity[t+InhibitionDelay]*InhLearnFac;
//    Dout(dc::layer,  "Squer= " << squer[t] << "  Mquer= " << mquer[t] << "");

        u[i] *= TauSfac;
    }

    for (i=0;i<N;i++)
    {
        Input[i] = 0.0;   // reset the input
        Input1[i] = 0.0;
    }
}


int lif2layer::prepare(int sec)
{
    int i,j;
    // prepare for the next MacroTimeStep
    layer::prepare(sec);
    for (i=0;i<N;++i) LastSpike[i] -= MacroTimeStep;
    for (i=0;i<InhibitionDelay+1;i++)
    {
        mquer[i] = mquer[MacroTimeStep+i];
        squer[i] = squer[MacroTimeStep+i];
        InhibitoryActivity[i] = InhibitoryActivity[MacroTimeStep+i];
    }
    for (i=InhibitionDelay+1;i<MacroTimeStep+InhibitionDelay+1;++i)
    {
        mquer[i]=0;
        squer[i]=0;
        InhibitoryActivity[i]=0;
    }
    float DinhCount=0;
    for (i=0;i<N;++i) DinhCount += Dinh[i];
    cout << "  DinhCount=" << DinhCount << "  DesAct=" << DesiredActivity << "Act=" << mquer[0] << " ";
}

int lif2layer::SetLearnPti(bool set)
{
    learnpti=set;
}

/////////////////////////////////

MMN01Layer::MMN01Layer(int n, float _TauThres, float _RestingThres, float _ThresInc, float _TauFeeding): layer(n), RestingThreshold(_RestingThres), ThresInc(_ThresInc), inh(0)
{
    // from class layer:
    // v: FeedingPot
    // u: Threshold
    ThresholdFac = exp(-dt/_TauThres);
    FeedingFac = exp(-dt/_TauFeeding);
    for (int i=0;i<N;++i) {
        v[i]=0;
        u[i]=0;
    }
}

MMN01Layer::~MMN01Layer()
{
    delete[] inh;
}

int MMN01Layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = 2*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(u[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}


float* MMN01Layer::GetInputPointer(csimInputChannel)
{
    return v;
}


int MMN01Layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        v[i] *= FeedingFac;
        u[i] *= ThresholdFac;
//    if ((i==0) && (rec != 0)) rec->record(dt*TotalTime, v[i], u[i]);

        if (v[i]>=RestingThreshold + u[i])                    // did it fire?
        {
            u[i] += ThresInc;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }
//      for (i=0;i<N;i++)
//      {
//    Input[i] = 0; // reset input
//      }

}



/////////////////////////////////

MMN02Layer::MMN02Layer(int n, float _TauThres, float _RestingThres, float _ThresInc, float _TauFeeding, float _TauInhibition): MMN01Layer(n, _TauThres, _RestingThres, _ThresInc, _TauFeeding)
{
    InhibitionFac = exp(-dt/_TauInhibition);
    inh = new float [N];
    for (int i=0;i<N;++i) {
        inh[i]=0;
    }
}

MMN02Layer::~MMN02Layer()
{
}

int MMN02Layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = 3*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(inh[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(u[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}


float* MMN02Layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return v;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return inh;
    } break;
    default:
        return v;
    }
    return 0;
}



int MMN02Layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        v[i] *= FeedingFac;
        v[i] += gsl_ran_gaussian(gslr, NoiseSigma);
        u[i] *= ThresholdFac;
        inh[i] *= InhibitionFac;
//    if ((i==0) && (rec != 0)) rec->record(dt*TotalTime, v[i], u[i]);

        if ((v[i]-inh[i])>=RestingThreshold + u[i])                   // did it fire?
        {
            u[i] += ThresInc;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

}


///////////////////////////////////////
/////////////////////////////////

MMN03Layer::MMN03Layer(int n, float _TauThres, float _RestingThres, float _ThresInc, float _TauFeeding, float _TauInhibition, float _TauLinking): MMN02Layer(n, _TauThres, _RestingThres, _ThresInc, _TauFeeding, _TauInhibition), TauLinking(_TauLinking)
{
    Dout(dc::layer,  "MMN03Layer::MMN03Layer, n=" << n << "");
    fflush(stdout);
    LinkingFac = exp(-dt/_TauLinking);
    Linking = new float [N];
    for (int i=0;i<N;++i) {
        Linking[i]=0;
    }
}

MMN03Layer::~MMN03Layer()
{
}

int MMN03Layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "MMN03" << "\"/> \n";
    sstr << "<TauLinking value=\"" << TauLinking << "\"/> \n";
    layer::WriteSimInfo(fw, sstr.str());
}


int MMN03Layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = 4*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(inh[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(u[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(Linking[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}


float* MMN03Layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return v;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return inh;
    } break;
    case csimInputChannel_Linking: {
        Dout(dc::layer,  "Linking Input");
        return Linking;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "NMDA/AMPA input requested, MMN03Layer has no NMDA current, using Linking Input instead");
        return Linking;
    } break;
    default:
        return v;
    }
    return 0;
}


int MMN03Layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        v[i] *= FeedingFac;                       // Feeding potential
        Linking[i] *= LinkingFac;              // linking potential
        // v[i] += gsl_ran_gaussian(gslr, NoiseSigma);    // Gausssches Rauschen
        u[i] *= ThresholdFac;                          // dynamische Schwelle
        inh[i] *= InhibitionFac;                  //inhibitorisches Potential

        if ((((Linking[i]+1)*v[i]-inh[i])+gsl_ran_gaussian(gslr,NoiseSigma)) >=RestingThreshold + u[i]) //did it fire?
        {
            u[i] += ThresInc;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

}


int MMN03Layer::reset(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    for (int i=0;i<N;++i) {
        v[i]=0;
        u[i]=0;
        inh[i]=0;
        Linking[i]=0;
    }

    // delete all spikes which are still in the delay line
    // firings with firing times higher than t-MaximumDelay
    int MaximumDelay = MainSimLoop->GetMaximumDelay();
    int k=N_firings-1;
    while (t-firings[k][0]<MaximumDelay) k--;
    N_firings = k+1;
}

/////////////////////////////////////////


MMN04Layer::MMN04Layer(int n, float _TauThres, float _RestingThres, float _ThresInc, float _TauFeeding, float _TauInhibition, float _TauLinking,float _TauEnduFeeding,float _TauEnduLinking): MMN03Layer(n, _TauThres, _RestingThres, _ThresInc, _TauFeeding, _TauInhibition, _TauLinking), TauEnduFeeding(_TauEnduFeeding), TauEnduLinking(_TauEnduLinking)
{
    Dout(dc::layer,  "MMN04Layer::MMN04Layer, n=" << n << "");
    fflush(stdout);
    EnduFeedingFac = exp(-dt/_TauEnduFeeding);
    EnduLinkingFac = exp(-dt/_TauEnduLinking);
    endu = new float [N];
    enduL = new float [N];
    v_Self = new float [N];
    for (int i=0;i<N;++i) {
        endu[i]=0;
        // enduL[i]=0;
        v_Self[i]=0;
    }
}

MMN04Layer::~MMN04Layer()
{
}

int MMN04Layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "MMN04" << "\"/> \n";
    sstr << "<TauEnduFeeding value=\"" << TauEnduFeeding << "\"/> \n";
    sstr << "<TauEnduLinking value=\"" << TauEnduLinking << "\"/> \n";
    layer::WriteSimInfo(fw, sstr.str());
}


int MMN04Layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = 4*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(inh[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(u[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(Linking[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(endu[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(enduL[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}


float* MMN04Layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return v;
    }  break;
    case csimInputChannel_AMPA2: {
        Dout(dc::layer,  "Excitatory Input");
        return v_Self;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return inh;
    } break;
    case csimInputChannel_Linking: {
        Dout(dc::layer,  "Linking Input");
        return Linking;
    } break;
    case csimInputChannel_Endu: {
        Dout(dc::layer,  "Enduring Input");
        return endu;
    } break;
    case csimInputChannel_EnduLinking: {
        Dout(dc::layer,  "Enduring Linking-Input");
        return enduL;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "NMDA/AMPA input requested, MMN03Layer has no NMDA current, using Linking Input instead");
        return Linking;
    } break;
    default:
        return v;
    }
    return 0;
}


int MMN04Layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        v[i] *= FeedingFac;// Feeding potential
        v_Self[i] *= FeedingFac;
        Linking[i] *= LinkingFac;           // linking potential
        // v[i] += gsl_ran_gaussian(gslr, NoiseSigma);    // Gausssches Rauschen
        u[i] *= ThresholdFac;                          // dynamische Schwelle
        inh[i] *= InhibitionFac;                  //inhibitorisches Potential
        endu[i] *= EnduFeedingFac;            // enduring Potential
        //enduL[i] *= EnduLinkingFac;   // Enduring Linking Potential

        //Wenn paramater von 0908 betrachtet werden sollen: if (((( (2*Linking[i])*(2*v[i]))+endu[i]+v_Self[i]-inh[i])+gsl_ran_gaussian(gslr,NoiseSigma)) >=RestingThreshold + u[i]) //did it fire?

        if (((3*Linking[i]*2*v[i])+endu[i]+v_Self[i]-inh[i]+gsl_ran_gaussian(gslr,NoiseSigma)) >=RestingThreshold + u[i]) //did it fire?
        {
            u[i] += ThresInc;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

}


int MMN04Layer::reset(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    for (int i=0;i<N;++i) {
        v[i]=0;
        u[i]=0;
        inh[i]=0;
        Linking[i]=0;
        endu[i]=0;
        // enduL[i]=0;
    }

    // delete all spikes which are still in the delay line
    // firings with firing times higher than t-MaximumDelay
    int MaximumDelay = MainSimLoop->GetMaximumDelay();
    int k=N_firings-1;
    while (t-firings[k][0]<MaximumDelay) k--;
    N_firings = k+1;
}












MMN05Layer::MMN05Layer(int n, float _TauThres, float _RestingThres, float _ThresInc, float _TauFeeding, float _TauInhibition, float _TauLinking,float _TauEnduFeeding,float _TauEnduLinking): MMN03Layer(n, _TauThres, _RestingThres, _ThresInc, _TauFeeding, _TauInhibition, _TauLinking), TauEnduFeeding(_TauEnduFeeding), TauEnduLinking(_TauEnduLinking)
{
    Dout(dc::layer,  "MMN05Layer::MMN05Layer, n=" << n << "");
    fflush(stdout);
    EnduFeedingFac = exp(-dt/_TauEnduFeeding);
    EnduLinkingFac = exp(-dt/_TauEnduLinking);
    endu = new float [N];
    enduL = new float [N];
    v_Self = new float [N];
    for (int i=0;i<N;++i) {
        endu[i]=0;
        // enduL[i]=0;
        v_Self[i]=0;
    }
}

MMN05Layer::~MMN05Layer()
{
}

int MMN05Layer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "MMN05" << "\"/> \n";
    sstr << "<TauEnduFeeding value=\"" << TauEnduFeeding << "\"/> \n";
    sstr << "<TauEnduLinking value=\"" << TauEnduLinking << "\"/> \n";
    layer::WriteSimInfo(fw, sstr.str());
}


int MMN05Layer::StartBinRec(int nobserve, int StartNumber)
{
    int NumObserve = 4*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(inh[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(u[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(Linking[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(endu[i]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(enduL[i]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}


float* MMN05Layer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return v;
    }  break;
    case csimInputChannel_AMPA2: {
        Dout(dc::layer,  "Excitatory Input");
        return v_Self;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return inh;
    } break;
    case csimInputChannel_Linking: {
        Dout(dc::layer,  "Linking Input");
        return Linking;
    } break;
    case csimInputChannel_Endu: {
        Dout(dc::layer,  "Enduring Input");
        return endu;
    } break;
    case csimInputChannel_EnduLinking: {
        Dout(dc::layer,  "Enduring Linking-Input");
        return enduL;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "NMDA/AMPA input requested, MMN03Layer has no NMDA current, using Linking Input instead");
        return Linking;
    } break;
    default:
        return v;
    }
    return 0;
}


int MMN05Layer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;
    if (BinRec) BinRec->record();
    for (i=0;i<N;i++)
    {
        v[i] *= FeedingFac;// Feeding potential
        v_Self[i] *= FeedingFac;
        Linking[i] *= LinkingFac;           // linking potential
        // v[i] += gsl_ran_gaussian(gslr, NoiseSigma);    // Gausssches Rauschen
        u[i] *= ThresholdFac;                          // dynamische Schwelle
        inh[i] *= InhibitionFac;                  //inhibitorisches Potential
        endu[i] *= EnduFeedingFac;            // enduring Potential
        //enduL[i] *= EnduLinkingFac;   // Enduring Linking Potential

        //Wenn paramater von 0908 betrachtet werden sollen: if (((( (2*Linking[i])*(2*v[i]))+endu[i]+v_Self[i]-inh[i])+gsl_ran_gaussian(gslr,NoiseSigma)) >=RestingThreshold + u[i]) //did it fire?

        if ((1+Linking[i])*(v[i]+endu[i]+v_Self[i]-inh[i])+gsl_ran_gaussian(gslr,NoiseSigma) >=RestingThreshold + u[i]) //did it fire?
        {
            u[i] += ThresInc;
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            if (N_firings == N_firings_max) {
                ThrowTooManySpikes(t);
            }
        }
    }

}


int MMN05Layer::reset(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    for (int i=0;i<N;++i) {
        v[i]=0;
        u[i]=0;
        inh[i]=0;
        Linking[i]=0;
        endu[i]=0;
        enduL[i]=0;
    }

    // delete all spikes which are still in the delay line
    // firings with firing times higher than t-MaximumDelay
    int MaximumDelay = MainSimLoop->GetMaximumDelay();
    int k=N_firings-1;
    while (t-firings[k][0]<MaximumDelay) k--;
    N_firings = k+1;
}


///////////////////////////////////////
AmpaNmdaGabaChannels::AmpaNmdaGabaChannels(
    int N, float _TauAmpaRise, float _TauAmpaFall, float _TauNmdaRise, float _TauNmdaFall
    ,float _TauGabaARise, float _TauGabaAFall, float _TauGabaBRise, float _TauGabaBFall
    ,float _GabaABRatio, float _NmdaAmpaRatio):
        TauAmpaRise(_TauAmpaRise)
        ,TauAmpaFall(_TauAmpaFall)
        ,TauNmdaRise(_TauNmdaRise)
        ,TauNmdaFall(_TauNmdaFall)
        ,TauGabaARise(_TauGabaARise)
        ,TauGabaAFall(_TauGabaAFall)
        ,TauGabaBRise(_TauGabaBRise)
        ,TauGabaBFall(_TauGabaBFall)
        ,GabaABRatio(_GabaABRatio)
        ,NmdaAmpaRatio(_NmdaAmpaRatio)
        ,GabaAReversePot(-70)
        ,GabaBReversePot(-90)
        ,NaReversePot(0)
        ,NNeurons(N)
{

    MainSimLoop = GetGlobalSimLoop();
    DeltaT = MainSimLoop->GetDeltaT();

    SetTauNmda(TauNmdaRise, TauNmdaFall);
    SetTauAmpa(TauAmpaRise, TauAmpaFall);
    SetTauGabaA(TauGabaARise, TauGabaAFall);
    SetTauGabaB(TauGabaBRise, TauGabaBFall);


// initialize arrays
    AmpaRisePot = new float [N];
    AmpaFallPot = new float [N];
    AmpaPot = new float [N];
    NmdaRisePot = new float [N];
    NmdaFallPot = new float [N];
    GabaARisePot = new float [N];
    GabaAFallPot = new float [N];
    GabaBRisePot = new float [N];
    GabaBFallPot = new float [N];
    InhInput = new float [N];
    NmdaAmpaInput = new float [N];

    for (int i=0;i<N;++i) {
        AmpaRisePot[i]=0;
        AmpaFallPot[i]=0;
        AmpaPot[i]=0;
        NmdaRisePot[i]=0;
        NmdaFallPot[i]=0;
        GabaARisePot[i]=0;
        GabaAFallPot[i]=0;
        GabaBRisePot[i]=0;
        GabaBFallPot[i]=0;
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;
    }
}

AmpaNmdaGabaChannels::~AmpaNmdaGabaChannels()
{
    delete [] AmpaRisePot;
    delete [] AmpaFallPot;
    delete [] AmpaPot;
    delete [] NmdaRisePot;
    delete [] NmdaFallPot;
    delete [] GabaARisePot;
    delete [] GabaAFallPot;
    delete [] GabaBRisePot;
    delete [] GabaBFallPot;
    delete [] InhInput;
    delete [] NmdaAmpaInput;
}

long AmpaNmdaGabaChannels::calcMemConsumption()
{
    long MemSum=0;
    MemSum += sizeof(float)*NNeurons*11;
    return MemSum;
}

void AmpaNmdaGabaChannels::reset()
{
    for (int i=0;i<NNeurons;++i) {
        AmpaRisePot[i]=0;
        AmpaFallPot[i]=0;
        AmpaPot[i]=0;
        NmdaRisePot[i]=0;
        NmdaFallPot[i]=0;
        GabaARisePot[i]=0;
        GabaAFallPot[i]=0;
        GabaBRisePot[i]=0;
        GabaBFallPot[i]=0;
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;
    }
}

void AmpaNmdaGabaChannels::SetTauNmda(float rise, float fall)
{
    TauNmdaRise = rise;
    TauNmdaFall = fall;
    TauNmdaRiseFac = exp(-DeltaT/TauNmdaRise);
    TauNmdaFallFac = exp(-DeltaT/TauNmdaFall);
    NmdaCorr = DoubleExpCorrectionFactor(TauNmdaRise, TauNmdaFall);
}

void AmpaNmdaGabaChannels::SetTauAmpa(float rise, float fall)
{
    TauAmpaRise = rise;
    TauAmpaFall = fall;
    TauAmpaRiseFac = exp(-DeltaT/TauAmpaRise);
    TauAmpaFallFac = exp(-DeltaT/TauAmpaFall);
    AmpaCorr = DoubleExpCorrectionFactor(TauAmpaRise, TauAmpaFall);
}

void AmpaNmdaGabaChannels::SetTauGabaA(float rise, float fall)
{
    TauGabaARise = rise;
    TauGabaAFall = fall;
    TauGabaARiseFac = exp(-DeltaT/TauGabaARise);
    TauGabaAFallFac = exp(-DeltaT/TauGabaAFall);
    GabaACorr = DoubleExpCorrectionFactor(TauGabaARise, TauGabaAFall);
}

void AmpaNmdaGabaChannels::SetTauGabaB(float rise, float fall)
{
    TauGabaBRise = rise;
    TauGabaBFall = fall;
    TauGabaBRiseFac = exp(-DeltaT/TauGabaBRise);
    TauGabaBFallFac = exp(-DeltaT/TauGabaBFall);
    GabaBCorr = DoubleExpCorrectionFactor(TauGabaBRise, TauGabaBFall);
}

void AmpaNmdaGabaChannels::SetNmdaAmpaRatio(float ratio)
{
    // ratio defined as psp peak ratio at xx mV
    // R*PotNmda*INmda(voltage) == PotAmpa*IAmpa(voltage)*ratio
    // R = PotAmpa*IAmpa(voltage)*ratio / (PotNmda*INmda(voltage))
    // PotNmda == PotAmpa ==>
    // R = IAmpa(voltage)*ratio / INmda(voltage)
    NmdaAmpaEPSPRatio=ratio;
    float voltage = 60; // mV
    float NaRevPot = 0; //mV
    float AmpaCurrent = NaRevPot-voltage; // if conductance gAmpa == 1
    float NmdaCurrent = INmda(voltage);
    NmdaAmpaRatio = AmpaCurrent*NmdaAmpaEPSPRatio/NmdaCurrent;
}

void AmpaNmdaGabaChannels::WriteSimInfo(stringstream &sstr)
{
    sstr << "<Ampa "
    << "TauRise=\"" << TauAmpaRise << "\" "
    << "TauFall=\"" << TauAmpaFall << "\" "
    << "AmpaCorr=\"" << AmpaCorr << "\" "
    << "/> \n";
    sstr << "<GabaA "
    << "TauRise=\"" << TauGabaARise << "\" "
    << "TauFall=\"" << TauGabaAFall << "\" "
    << "ReversePot=\"" << GabaAReversePot << "\" "
    << "/> \n";
    sstr << "<GabaB "
    << "TauRise=\"" << TauGabaBRise << "\" "
    << "TauFall=\"" << TauGabaBFall << "\" "
    << "ReversePot=\"" << GabaBReversePot << "\" "
    << "/> \n";
    sstr << "<Nmda "
    << "TauRise=\"" << TauNmdaRise << "\" "
    << "TauFall=\"" << TauNmdaFall << "\" "
    << "NmdaAmpaRatio=\"" << NmdaAmpaRatio << "\" "
    << "NmdaAmpaEPSPRatio=\"" << NmdaAmpaEPSPRatio << "\" "
    << "/> \n";
}

/////////////////////////////

DecoLifLayer::DecoLifLayer(int n, DecoParas Paras): layer(n), AmpaNmdaGabaChannels(n), RestingPot(Paras.RestingPot), ResetPot(Paras.ResetPot), Threshold(Paras.Threshold), CMembrane(Paras.CMembrane), GLeak(Paras.GLeak), AbsoluteRefractoryPeriod(Paras.AbsoluteRefractoryPeriod), SpikePot(Paras.SpikePot)
{
    Refractory = new int [N];
    for (int i=0;i<N;++i) {
        v[i] = RestingPot;
        Refractory[i]=0;
    }
    CMembraneFac = 0.001*dt/CMembrane;
    Dout(dc::layer,  "cmemfac=" << CMembraneFac << "");
    arf = int(round(AbsoluteRefractoryPeriod/dt));

#ifdef USE_RNG_THREAD
    // initialize random number generator thread queue
    Dout(dc::layer,  "initialize random number generator thread");
    mRndNumQueue = RngParaQueuePool::getRngQueuePool()->getRngQueue(kRngPoisson, PoissonLambda,10,2000);
#endif //USE_RNG_THREAD      
}


DecoLifLayer::~DecoLifLayer()
{
    delete [] Refractory;
}

int DecoLifLayer::WriteSimInfo(fstream &fw)
{
    stringstream sstr;
    sstr << "<LayerType value=\"" << "DecoLif" << "\"/> \n";
    sstr << "<NoiseAmplitude value=\"" << NoiseAmplitude << "\"/> \n";
    sstr << "<BalancedInhibition value=\"" << BalancedInhibition << "\"/> \n";
    sstr << "<RestingPot value=\"" << RestingPot << "\"/> \n";
    sstr << "<ResetPot value=\"" << ResetPot << "\"/> \n";
    sstr << "<Threshold value=\"" << Threshold << "\"/> \n";
    sstr << "<CMembrane value=\"" << CMembrane << "\"/> \n";
    sstr << "<GLeak value=\"" << GLeak << "\"/> \n";
    sstr << "<AbsoluteRefractoryPeriod value=\"" << AbsoluteRefractoryPeriod << "\"/> \n";
    AmpaNmdaGabaChannels::WriteSimInfo(sstr);
    layer::WriteSimInfo(fw, sstr.str());
}

int DecoLifLayer::proceede(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    int i,j,k;

    last_N_firings = N_firings;

    if (BinRec) BinRec->record();

    for (i=0;i<N;i++)
    {
//    Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));

#ifdef USE_RNG_THREAD
        // use random number generator thread
        Input[i] += NoiseAmplitude*mRndNumQueue->getRandomNumber();
#else
        // use no thread for rng
        Input[i] += NoiseAmplitude*gsl_ran_poisson(gslr, PoissonLambda);
#endif //USE_RNG_THREAD      

        Input[i] += NmdaAmpaInput[i];

        Input[i] *= AmpaCorr;
        InhInput[i] *= GabaACorr;
        NmdaAmpaInput[i] *= NmdaCorr;

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaARisePot[i] += InhInput[i];
        GabaAFallPot[i] += InhInput[i];
        GabaBRisePot[i] +=GabaABRatio*InhInput[i];
        GabaBFallPot[i] += GabaABRatio*InhInput[i];

        // if (i == 0) cerr << "Nmda=" <<  NmdaRisePot[i] << " " << NmdaFallPot[i] << "\n";

//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        AmpaPot[i] = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float GabaA = GabaAFallPot[i] - GabaARisePot[i];
        float GabaB = GabaBFallPot[i] - GabaBRisePot[i];

        v[i] += CMembraneFac*(-GLeak*(v[i]-RestingPot)
                              - AmpaPot[i]*(v[i]-NaReversePot)
                              + INmda(v[i])*Nmda
                              - (v[i]-GabaAReversePot)*(GabaA+BalancedInhibition)
                              - (v[i]-GabaBReversePot)*GabaB
                             );

        if (Refractory[i]>0) {
            --Refractory[i];
            v[i] = ResetPot;                 // voltage reset
        } else {
            if (v[i]>=Threshold) {   // did it fire?
//          if (i == 0) cerr << N_firings << "spike\n"; fflush(stdout);
//          v[i] = ResetPot;                    // voltage reset
                v[i] = SpikePot;                    // fake spike for LIF Neuron
                Refractory[i] = arf;
                firings[N_firings  ][0]=t;
                firings[N_firings++][1]=i;
                if (N_firings == N_firings_max) {
                    ThrowTooManySpikes(t);
                }
            }
        }

    }
    for (i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaARisePot[i] *= TauGabaARiseFac;
        GabaAFallPot[i] *= TauGabaAFallFac;
        GabaBRisePot[i] *= TauGabaBRiseFac;
        GabaBFallPot[i] *= TauGabaBFallFac;
        //if (GabaARisePot[i] != 0 && GabaARisePot[i] < numeric_limits<float>::min( )) cerr << "PERFORMANCE-WARNING: denormal floating point number\n";
    }
}

int DecoLifLayer::proceedeInThread(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;


}

int DecoLifLayer::calculateInThread(int Start_Index, int Stop_Index)
{
    int i,j,k;

    last_N_firings = N_firings;

    if (BinRec) BinRec->record();

    for (i=Start_Index;i<Stop_Index;i++)
    {
//    Input[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        Input[i] += NoiseAmplitude*gsl_ran_poisson(gslr, PoissonLambda);
        Input[i] += NmdaAmpaInput[i];

        Input[i] *= AmpaCorr;
        InhInput[i] *= GabaACorr;
        NmdaAmpaInput[i] *= NmdaCorr;

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaARisePot[i] += InhInput[i];
        GabaAFallPot[i] += InhInput[i];
        GabaBRisePot[i] +=GabaABRatio*InhInput[i];
        GabaBFallPot[i] += GabaABRatio*InhInput[i];

        // if (i == 0) cerr << "Nmda=" <<  NmdaRisePot[i] << " " << NmdaFallPot[i] << "\n";

//    InhSynPot[i] += abs(gsl_ran_gaussian(gslr, NoiseSigma));
        AmpaPot[i] = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float GabaA = GabaAFallPot[i] - GabaARisePot[i];
        float GabaB = GabaBFallPot[i] - GabaBRisePot[i];

        v[i] += CMembraneFac*(-GLeak*(v[i]-RestingPot)
                              - AmpaPot[i]*(v[i]-NaReversePot)
                              + INmda(v[i])*Nmda
                              - (v[i]-GabaAReversePot)*(GabaA+BalancedInhibition)
                              - (v[i]-GabaBReversePot)*GabaB
                             );

        if (Refractory[i]>0) {
            --Refractory[i];
            v[i] = ResetPot;                 // voltage reset
        } else {
            if (v[i]>=Threshold) {   // did it fire?
//          if (i == 0) cerr << N_firings << "spike\n"; fflush(stdout);
//          v[i] = ResetPot;                    // voltage reset
                v[i] = SpikePot;                    // fake spike for LIF Neuron
                Refractory[i] = arf;
                // MUTEX for protecting firings!!
                firings[N_firings  ][0]=t;
                firings[N_firings++][1]=i;
                if (N_firings == N_firings_max) {
                    ThrowTooManySpikes(t);
                }
            }
        }

    }
    for (i=Start_Index;i<Stop_Index;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaARisePot[i] *= TauGabaARiseFac;
        GabaAFallPot[i] *= TauGabaAFallFac;
        GabaBRisePot[i] *= TauGabaBRiseFac;
        GabaBFallPot[i] *= TauGabaBFallFac;
        if (GabaARisePot[i] != 0 && GabaARisePot[i] < 0.00000000001) cerr << "PERFORMANCE-WARNING: small floating point number\n";
    }
}


int DecoLifLayer::reset(int TotalTime)
{
    int t = TotalTime % MacroTimeStep;
    AmpaNmdaGabaChannels::reset();
    for (int i=0;i<N;++i) {
        Refractory[i]=0;
        v[i] = RestingPot;
    }

    // delete all spikes which are still in the delay line
    // firings with firing times higher than t-MaximumDelay
    int MaximumDelay = SimElement::MainSimLoop->GetMaximumDelay();
    int k=N_firings-1;
    while (t-firings[k][0]<MaximumDelay) k--;
    N_firings = k+1;
}


float* DecoLifLayer::GetInputPointer(csimInputChannel InputNumber)
{
    Dout(dc::layer,  "InputNumber=" << int(InputNumber) << "");
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return NmdaAmpaInput;
    }  break;
    default:
        return Input;
    }
    return 0;
}

float* DecoLifLayer::GetPspPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case  csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return AmpaPot;
    }  break;
    default:
        return AmpaPot;
    }
    return 0;
}

long DecoLifLayer::calcMemoryConsumption()
{
    long MemSum=0;
    MemSum += calcMemConsumption();
    MemSum += layer::calcMemoryConsumption();
    MemSum += sizeof(int)*N; //int* Refractory;

}


/** Starts binary recording of membrane potential and synaptic potentials
 *
 * @param nobserve number of neurons to observe
 * @param NNumber starting number. Neurons in the range [NNumber,NNumber+nobserve-1] are recorded
 * @return
 */
int DecoLifLayer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 9*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(AmpaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(AmpaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(GabaARisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+4*nobserve] = &(GabaAFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+5*nobserve] = &(NmdaRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+6*nobserve] = &(NmdaFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+7*nobserve] = &(GabaBRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+8*nobserve] = &(GabaBFallPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

//Hodgkin-Huxley-Model-Neuron

HodgHuxLayer::HodgHuxLayer(int n,float vl,float vna,float vk,float gl,float gna,float gk,float cm,float vshift,float spikedetectionthreshold,float mmin,float mmax,float mstart, int mn) : layer(n), AmpaNmdaGabaChannels(n),gL(gl),gNa(gna),gK(gk),CM(cm),VShift(vshift),Mstart(mstart),MN(mn)
{
    SpikeDetectionThreshold=spikedetectionthreshold-VShift;
    VL=vl-VShift;
    VNa=vna-VShift;
    VK=vk-VShift;
    Mmin=mmin-VShift;
    Mmax=mmax-VShift;
    hParas[0]=0.07;
    hParas[1]=20;
    hParas[2]=1;
    hParas[3]=30-VShift;
    hParas[4]=0.1;
    hParas[5]=1;
    mParas[0]=0.1;
    mParas[1]=24-VShift;
    mParas[2]=0.1;
    mParas[3]=1;
    mParas[4]=4;
    mParas[5]=18;
    nParas[0]=0.01;
    nParas[1]=10-VShift;
    nParas[2]=0.1;
    nParas[3]=1;
    nParas[4]=0.125;
    nParas[5]=80;
    hAlphaArr=new float[MN];
    mAlphaArr=new float[MN];
    nAlphaArr=new float[MN];
    hBetaArr=new float[MN];
    mBetaArr=new float[MN];
    nBetaArr=new float[MN];
    setxAlphaBeta(HodgHux_h,hParas);
    setxAlphaBeta(HodgHux_m,mParas);
    setxAlphaBeta(HodgHux_n,nParas);

    hNa=new float[N];
    mNa=new float[N];
    nK=u;
    M=v;
    int index=M2index(Mstart);
    for (int i=0;i<N;i++)
    {
        hNa[i]=hAlphaArr[index]/(hAlphaArr[index]+hBetaArr[index]);
        mNa[i]=mAlphaArr[index]/(mAlphaArr[index]+mBetaArr[index]);
        nK[i]=nAlphaArr[index]/(nAlphaArr[index]+nBetaArr[index]);
        M[i]=Mstart;
    }

    alreadySpiked=new bool[N];
    for (int i=0;i<N;i++) alreadySpiked[i]=false;

    NaReversePot=VNa;
    GabaAReversePot=VK;
    GabaBReversePot=VL;
}

HodgHuxLayer::~HodgHuxLayer()
{
    delete[] hAlphaArr;
    delete[] mAlphaArr;
    delete[] nAlphaArr;
    delete[] hBetaArr;
    delete[] mBetaArr;
    delete[] nBetaArr;

    delete[] hNa;
    delete[] mNa;

    delete[] alreadySpiked;
}

int HodgHuxLayer::proceede(int t)
{
    //Dout(dc::layer, "Anfang von proceede Layer"); fflush(stdout);
    //Dout(dc::layer, M[0]);
    if (BinRec) BinRec->record();
    for (int i=0;i<N;i++)
    {
//    float ma,ha,na,mb,hb,nb;
//    na = nParas[0]*(nParas[1]-M[i])/(exp((nParas[1]-M[i])*nParas[2])-nParas[3]);
//    nb = nParas[4]*exp(-(M[i]+VShift)/nParas[5]);
//    ma = mParas[0]*(mParas[1]-M[i])/(exp((mParas[1]-M[i])*mParas[2])-mParas[3]);
//    mb = mParas[4]*exp(-(M[i]+VShift)/mParas[5]);
//    ha = hParas[0]*exp(-(M[i]+VShift)/hParas[1]);
//    hb = hParas[2]/(exp((hParas[3]-M[i])*hParas[4])+hParas[5]);
//    mNa[i]+=dt*(ma*(1-mNa[i])-mb*mNa[i]);
//    hNa[i]+=dt*(ha*(1-hNa[i])-hb*hNa[i]);
//    nK[i]+=dt*(na*(1-nK[i])-nb*nK[i]);

        int index=M2index(M[i]);
        mNa[i]+=dt*(mAlphaArr[index]*(1-mNa[i])-mBetaArr[index]*mNa[i]);
        hNa[i]+=dt*(hAlphaArr[index]*(1-hNa[i])-hBetaArr[index]*hNa[i]);
        nK[i]+=dt*(nAlphaArr[index]*(1-nK[i])-nBetaArr[index]*nK[i]);

        M[i]-=dt/CM*(gL*(M[i]-VL)
                     +gNa*(mNa[i])*(mNa[i])*(mNa[i])*(hNa[i])*(M[i]-VNa)
                     +gK*(nK[i])*(nK[i])*(nK[i])*(nK[i])*(M[i]-VK)
                     //       -Input[i]
                    );

        //Input[i] += NoiseAmplitude*gsl_ran_poisson(gslr, PoissonLambda);
        Input[i] += NmdaAmpaInput[i];

        //Dout(dc::layer, "Vor ...Corr"); fflush(stdout);
        Input[i] *= AmpaCorr;
        InhInput[i] *= GabaACorr;
        NmdaAmpaInput[i] *= NmdaCorr;
        //Dout(dc::layer, "...Corr klappt"); fflush(stdout);

        AmpaRisePot[i] += Input[i];
        AmpaFallPot[i] += Input[i];
        NmdaRisePot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        NmdaFallPot[i] += NmdaAmpaRatio*NmdaAmpaInput[i];
        GabaARisePot[i] += InhInput[i];
        GabaAFallPot[i] += InhInput[i];
        GabaBRisePot[i] +=GabaABRatio*InhInput[i];
        GabaBFallPot[i] += GabaABRatio*InhInput[i];
        //Dout(dc::layer, "...Input[i] klappt"); fflush(stdout);

        AmpaPot[i] = AmpaFallPot[i] - AmpaRisePot[i];
        float Nmda = NmdaFallPot[i] - NmdaRisePot[i];
        float GabaA = GabaAFallPot[i] - GabaARisePot[i];
        float GabaB = GabaBFallPot[i] - GabaBRisePot[i];
        M[i]+=dt/CM*(- AmpaPot[i]*(M[i]-NaReversePot)
                     + INmda(M[i])*Nmda
                     - (M[i]-GabaAReversePot)*(GabaA+BalancedInhibition)
                     - (M[i]-GabaBReversePot)*GabaB
                    );
        //Dout(dc::layer, "...Pot klappt"); fflush(stdout);

        if ((M[i]>SpikeDetectionThreshold) && (!alreadySpiked[i]))
        {
            firings[N_firings  ][0]=t;
            firings[N_firings++][1]=i;
            alreadySpiked[i]=true;
        }
        if (alreadySpiked[i] && M[i]<SpikeDetectionThreshold) alreadySpiked[i]=false;
    }
    //Dout(dc::layer, "Ende von proceede Layer"); fflush(stdout);
    for (int i=0;i<N;i++)
    {
        Input[i] = 0;     // reset input
        InhInput[i]=0;
        NmdaAmpaInput[i]=0;

        // decay synaptic potentials
        NmdaRisePot[i] *= TauNmdaRiseFac;
        NmdaFallPot[i] *= TauNmdaFallFac;
        AmpaRisePot[i] *= TauAmpaRiseFac;
        AmpaFallPot[i] *= TauAmpaFallFac;
        GabaARisePot[i] *= TauGabaARiseFac;
        GabaAFallPot[i] *= TauGabaAFallFac;
        GabaBRisePot[i] *= TauGabaBRiseFac;
        GabaBFallPot[i] *= TauGabaBFallFac;
    }
    //Dout(dc::layer, "wirklich Ende von proceede"); fflush(stdout);
}

void HodgHuxLayer::setxAlphaBeta(HodgHuxGate gate,float paras[6])
{
    int i;
    float vi;
    float (HodgHuxLayer::*parapointer)[6];
    switch (gate) {
    case HodgHux_h:
        for (i=0;i<MN;i++)
        {
            vi=Mmin+(Mmax-Mmin)*(i/(float)MN);
            hAlphaArr[i]=paras[0]*exp(-(vi+VShift)/paras[1]);
            hBetaArr[i]=paras[2]/(exp((paras[3]-vi)*paras[4])+paras[5]);
            parapointer=&HodgHuxLayer::hParas;
        }
        break;
    case HodgHux_m:
        for (i=0;i<MN;i++)
        {
            vi=Mmin+(Mmax-Mmin)*(i/(float)MN);
            mAlphaArr[i]=paras[0]*(paras[1]-vi)/(exp((paras[1]-vi)*paras[2])-paras[3]);
            mBetaArr[i]=paras[4]*exp(-(vi+VShift)/paras[5]);
            parapointer=&HodgHuxLayer::mParas;
        }
        break;
    case HodgHux_n:
        for (i=0;i<MN;i++)
        {
            vi=Mmin+(Mmax-Mmin)*(i/(float)MN);
            nAlphaArr[i]=paras[0]*(paras[1]-vi)/(exp((paras[1]-vi)*paras[2])-paras[3]);
            nBetaArr[i]=paras[4]*exp(-(vi+VShift)/paras[5]);
            parapointer=&HodgHuxLayer::nParas;
        }
        break;
    }
    for (i=0;i<6;i++) (this->*parapointer)[i]=paras[i];

    if (parapointer==&HodgHuxLayer::nParas)   cout<<"n ";
    if (parapointer==&HodgHuxLayer::mParas)   cout<<"m ";
    if (parapointer==&HodgHuxLayer::hParas)   cout<<"h ";
    Dout(dc::layer, paras[0]<<" "<<paras[1]<<" "<<paras[2]<<" "<<paras[3]<<" "<<paras[4]<<" "<<paras[5]);
    Dout(dc::layer, (this->*parapointer)[0]<<" "<<(this->*parapointer)[1]<<" "<<(this->*parapointer)[2]<<" "<<(this->*parapointer)[3]<<" "<<(this->*parapointer)[4]<<" "<<(this->*parapointer)[5]);
}

int HodgHuxLayer::M2index(float m)
{
    int zwischen=(int)( (m-Mmin)/(Mmax-Mmin)*MN );
    //Dout(dc::layer, "Angang von M2index. m="<<m<<" out="<<zwischen);
    if (m<Mmin)
    {
        //Dout(dc::layer, "Ende von M2index. Returned 0."); fflush(stdout);
        return 0;
    }
    else if (m>Mmax)
    {
        //Dout(dc::layer, "Ende von M2index. Returned MN-1."); fflush(stdout);
        return MN-1;
    }
    else
    {
        //Dout(dc::layer, "Ende von M2index. Returned zwischen."); fflush(stdout);
        return (int)( (m-Mmin)/(Mmax-Mmin)*MN );
    }
}

float* HodgHuxLayer::GetInputPointer(csimInputChannel InputNumber)
{
    switch (InputNumber) {
    case csimInputChannel_AMPA: {
        Dout(dc::layer,  "Excitatory Input");
        return Input;
    }  break;
    case csimInputChannel_GABAa: {
        Dout(dc::layer,  "Inhibitory Input");
        return InhInput;
    } break;
    case csimInputChannel_NMDA_AMPA: {
        Dout(dc::layer,  "Mixed: Excitatory and NMDA-Input");
        return NmdaAmpaInput;
    }  break;
    default:
        return Input;
    }
    return 0;
}

int HodgHuxLayer::WriteSimInfo(fstream &fw)
{
    fw << "<" << seTypeString << " id=\"" << IdNumber <<  "\" Type=\"" << seType << "\" Name=\"" << Name << "\"> \n";
    fw << "     <LayerType value=\"" << "Hodgkin-Huxley" << "\"/> \n";
    fw << "     <LayerDimensions N=\"" << N << "\" Nx=\"" << Nx << "\" Ny=\"" << Ny << "\"/> \n";
    fw << "</" << seTypeString << "> \n";
}

int HodgHuxLayer::StartBinRec(int nobserve, int NNumber)
{
    int NumObserve = 6*nobserve;
    float** Buffer = new float* [NumObserve];
    for (int i=0;i<nobserve; ++i) Buffer[i] = &(v[i+NNumber]);
//      for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(hNa[i+NNumber]);
//      for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(mNa[i+NNumber]);
//      for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(nK[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+nobserve] = &(AmpaPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+2*nobserve] = &(GabaARisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+3*nobserve] = &(GabaAFallPot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+4*nobserve] = &(GabaBRisePot[i+NNumber]);
    for (int i=0;i<nobserve; ++i) Buffer[i+5*nobserve] = &(GabaBFallPot[i+NNumber]);
    string FileName(MemPotFileName);
    BinRec = new BinRecorder(MacroTimeStep, NumObserve, nobserve, Buffer, dt, (DataDirectory+FileName).c_str());
}

/////////////////////////////////////////