123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765 |
- #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());
- }
- /////////////////////////////////////////
|