PetaVision  Alpha
LIFGap.cpp
1 /*
2  * LIFGap.cpp
3  *
4  * Created on: Jul 29, 2011
5  * Author: garkenyon
6  */
7 
8 #include "LIFGap.hpp"
9 #include "../connections/HyPerConn.hpp"
10 #include "../include/default_params.h"
11 #include "../include/pv_common.h"
12 #include "../io/fileio.hpp"
13 #include "utils/cl_random.h"
14 
15 #include <assert.h>
16 #include <cmath>
17 #include <float.h>
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 
22 void LIFGap_update_state_original(
23  const int nbatch,
24  const int numNeurons,
25  const float time,
26  const float dt,
27 
28  const int nx,
29  const int ny,
30  const int nf,
31  const int lt,
32  const int rt,
33  const int dn,
34  const int up,
35 
36  LIF_params *params,
37  taus_uint4 *rnd,
38 
39  float *V,
40  float *Vth,
41  float *G_E,
42  float *G_I,
43  float *G_IB,
44  float *GSynHead,
45  float *activity,
46 
47  const float *gapStrength);
48 
49 void LIFGap_update_state_beginning(
50  const int nbatch,
51  const int numNeurons,
52  const float time,
53  const float dt,
54 
55  const int nx,
56  const int ny,
57  const int nf,
58  const int lt,
59  const int rt,
60  const int dn,
61  const int up,
62 
63  LIF_params *params,
64  taus_uint4 *rnd,
65 
66  float *V,
67  float *Vth,
68  float *G_E,
69  float *G_I,
70  float *G_IB,
71  float *GSynHead,
72  float *activity,
73 
74  const float *gapStrength);
75 
76 void LIFGap_update_state_arma(
77  const int nbatch,
78  const int numNeurons,
79  const float time,
80  const float dt,
81 
82  const int nx,
83  const int ny,
84  const int nf,
85  const int lt,
86  const int rt,
87  const int dn,
88  const int up,
89 
90  LIF_params *params,
91  taus_uint4 *rnd,
92 
93  float *V,
94  float *Vth,
95  float *G_E,
96  float *G_I,
97  float *G_IB,
98  float *GSynHead,
99  float *activity,
100 
101  const float *gapStrength);
102 
103 namespace PV {
104 
105 LIFGap::LIFGap() { initialize_base(); }
106 
107 LIFGap::LIFGap(const char *name, HyPerCol *hc) {
108  initialize_base();
109  initialize(name, hc, "LIFGap_update_state");
110 }
111 
112 LIFGap::~LIFGap() { free(gapStrength); }
113 
114 int LIFGap::initialize_base() {
115  numChannels = 4;
116  gapStrength = NULL;
117  gapStrengthInitialized = false;
118  return PV_SUCCESS;
119 }
120 
121 // Initialize this class
122 /*
123  *
124  */
125 int LIFGap::initialize(const char *name, HyPerCol *hc, const char *kernel_name) {
126  int status = LIF::initialize(name, hc, kernel_name);
127  return status;
128 }
129 
130 void LIFGap::allocateConductances(int num_channels) {
131  LIF::allocateConductances(num_channels - 1); // CHANNEL_GAP doesn't have a conductance per se.
132  gapStrength = (float *)calloc((size_t)getNumNeuronsAllBatches(), sizeof(*gapStrength));
133  if (gapStrength == nullptr) {
134  Fatal().printf(
135  "%s: rank %d process unable to allocate memory for gapStrength: %s\n",
136  getDescription_c(),
137  parent->columnId(),
138  strerror(errno));
139  }
140 }
141 
142 void LIFGap::calcGapStrength() {
143  bool needsNewCalc = !gapStrengthInitialized;
144  if (!needsNewCalc) {
145  for (auto &c : recvConns) {
146  HyPerConn *conn = dynamic_cast<HyPerConn *>(c);
147  if (conn != nullptr) {
148  continue;
149  }
150  if (conn->getChannelCode() == CHANNEL_GAP && mLastUpdateTime < conn->getLastUpdateTime()) {
151  needsNewCalc = true;
152  break;
153  }
154  }
155  }
156  if (!needsNewCalc) {
157  return;
158  }
159 
160  for (int k = 0; k < getNumNeuronsAllBatches(); k++) {
161  gapStrength[k] = (float)0;
162  }
163  for (auto &c : recvConns) {
164  if (c == nullptr or c->getChannelCode() != CHANNEL_GAP) {
165  continue;
166  }
167  pvAssert(c->getPost() == this);
168  auto *weightUpdater = c->getComponentByType<BaseWeightUpdater>();
169  if (weightUpdater and weightUpdater->getPlasticityFlag() and parent->columnId() == 0) {
170  WarnLog().printf(
171  "%s: %s on CHANNEL_GAP has plasticity flag set to true\n",
172  getDescription_c(),
173  c->getDescription_c());
174  }
175  c->deliverUnitInput(gapStrength);
176  }
177  gapStrengthInitialized = true;
178 }
179 
180 Response::Status LIFGap::registerData(Checkpointer *checkpointer) {
181  auto status = LIF::registerData(checkpointer);
182  if (!Response::completed(status)) {
183  return status;
184  }
185  checkpointPvpActivityFloat(checkpointer, "gapStrength", gapStrength, false /*not extended*/);
186  return Response::SUCCESS;
187 }
188 
189 Response::Status LIFGap::readStateFromCheckpoint(Checkpointer *checkpointer) {
190  if (initializeFromCheckpointFlag) {
191  auto status = LIF::readStateFromCheckpoint(checkpointer);
192  if (!Response::completed(status)) {
193  return status;
194  }
195  readGapStrengthFromCheckpoint(checkpointer);
196  return Response::SUCCESS;
197  }
198  else {
199  return Response::NO_ACTION;
200  }
201 }
202 
203 void LIFGap::readGapStrengthFromCheckpoint(Checkpointer *checkpointer) {
204  checkpointer->readNamedCheckpointEntry(
205  std::string(name), std::string("gapStrength"), false /*not constant*/);
206 }
207 
208 Response::Status LIFGap::updateState(double time, double dt) {
209  calcGapStrength();
210 
211  const int nx = clayer->loc.nx;
212  const int ny = clayer->loc.ny;
213  const int nf = clayer->loc.nf;
214  const PVHalo *halo = &clayer->loc.halo;
215  const int nbatch = clayer->loc.nbatch;
216 
217  float *GSynHead = GSyn[0];
218  float *activity = clayer->activity->data;
219 
220  switch (method) {
221  case 'a':
222  LIFGap_update_state_arma(
223  nbatch,
224  getNumNeurons(),
225  time,
226  dt,
227  nx,
228  ny,
229  nf,
230  halo->lt,
231  halo->rt,
232  halo->dn,
233  halo->up,
234  &lParams,
235  randState->getRNG(0),
236  clayer->V,
237  Vth,
238  G_E,
239  G_I,
240  G_IB,
241  GSynHead,
242  activity,
243  gapStrength);
244  break;
245  case 'b':
246  LIFGap_update_state_beginning(
247  nbatch,
248  getNumNeurons(),
249  time,
250  dt,
251  nx,
252  ny,
253  nf,
254  halo->lt,
255  halo->rt,
256  halo->dn,
257  halo->up,
258  &lParams,
259  randState->getRNG(0),
260  clayer->V,
261  Vth,
262  G_E,
263  G_I,
264  G_IB,
265  GSynHead,
266  activity,
267  gapStrength);
268  break;
269  case 'o':
270  LIFGap_update_state_original(
271  nbatch,
272  getNumNeurons(),
273  time,
274  dt,
275  nx,
276  ny,
277  nf,
278  halo->lt,
279  halo->rt,
280  halo->dn,
281  halo->up,
282  &lParams,
283  randState->getRNG(0),
284  clayer->V,
285  Vth,
286  G_E,
287  G_I,
288  G_IB,
289  GSynHead,
290  activity,
291  gapStrength);
292  break;
293  default: break;
294  }
295  return Response::SUCCESS;
296 }
297 
298 } // namespace PV
299 
301 //
302 // implementation of LIF kernels
303 //
304 
305 inline float LIFGap_Vmem_derivative(
306  const float Vmem,
307  const float G_E,
308  const float G_I,
309  const float G_IB,
310  const float G_Gap,
311  const float V_E,
312  const float V_I,
313  const float V_IB,
314  const float sum_gap,
315  const float Vrest,
316  const float tau) {
317  float totalconductance = 1.0f + G_E + G_I + G_IB + sum_gap;
318  float Vmeminf = (Vrest + V_E * G_E + V_I * G_I + V_IB * G_IB + G_Gap) / totalconductance;
319  return totalconductance * (Vmeminf - Vmem) / tau;
320 }
321 
322 //
323 // update the state of a LIFGap layer (spiking)
324 //
325 // assume called with 1D kernel
326 //
327 // LIFGap_update_state_original uses an Euler scheme for V where the conductances over the entire
328 // timestep are taken to be the values calculated at the end of the timestep
329 // LIFGap_update_state_beginning uses a Heun scheme for V, using values of the conductances at both
330 // the beginning and end of the timestep. Spikes in the input are applied at the beginning of the
331 // timestep.
332 //
333 
334 void LIFGap_update_state_original(
335  const int nbatch,
336  const int numNeurons,
337  const float time,
338  const float dt,
339 
340  const int nx,
341  const int ny,
342  const int nf,
343  const int lt,
344  const int rt,
345  const int dn,
346  const int up,
347 
348  LIF_params *params,
349  taus_uint4 *rnd,
350  float *V,
351  float *Vth,
352  float *G_E,
353  float *G_I,
354  float *G_IB,
355  float *GSynHead,
356  float *activity,
357 
358  const float *gapStrength) {
359  int k;
360 
361  const float exp_tauE = expf(-dt / params->tauE);
362  const float exp_tauI = expf(-dt / params->tauI);
363  const float exp_tauIB = expf(-dt / params->tauIB);
364  const float exp_tauVth = expf(-dt / params->tauVth);
365 
366  const float dt_sec = 0.001f * dt; // convert to seconds
367 
368  for (k = 0; k < nx * ny * nf * nbatch; k++) {
369  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
370 
371  //
372  // kernel (nonheader part) begins here
373  //
374 
375  // local param variables
376  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
377 
378  // local variables
379  float l_activ;
380 
381  taus_uint4 l_rnd = rnd[k];
382 
383  float l_V = V[k];
384  float l_Vth = Vth[k];
385 
386  float l_G_E = G_E[k];
387  float l_G_I = G_I[k];
388  float l_G_IB = G_IB[k];
389  float l_gapStrength = gapStrength[k];
390 
391  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
392  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
393  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
394  float *GSynGap = &GSynHead[CHANNEL_GAP * nbatch * numNeurons];
395  float l_GSynExc = GSynExc[k];
396  float l_GSynInh = GSynInh[k];
397  float l_GSynInhB = GSynInhB[k];
398  float l_GSynGap = GSynGap[k];
399 
400  // define local param variables
401  //
402  tau = params->tau;
403  Vexc = params->Vexc;
404  Vinh = params->Vinh;
405  VinhB = params->VinhB;
406  Vrest = params->Vrest;
407 
408  VthRest = params->VthRest;
409  deltaVth = params->deltaVth;
410  deltaGIB = params->deltaGIB;
411 
412  // add noise
413  //
414 
415  l_rnd = cl_random_get(l_rnd);
416  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
417  l_rnd = cl_random_get(l_rnd);
418  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
419  }
420 
421  l_rnd = cl_random_get(l_rnd);
422  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
423  l_rnd = cl_random_get(l_rnd);
424  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
425  }
426 
427  l_rnd = cl_random_get(l_rnd);
428  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
429  l_rnd = cl_random_get(l_rnd);
430  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
431  }
432 
433  const float GMAX = 10.0f;
434  float tauInf, VmemInf;
435 
436  // The portion of code below uses the original method of calculating l_V.
437  l_G_E = l_GSynExc + l_G_E * exp_tauE;
438  l_G_I = l_GSynInh + l_G_I * exp_tauI;
439  l_G_IB = l_GSynInhB + l_G_IB * exp_tauIB;
440 
441  l_G_E = (l_G_E > GMAX) ? GMAX : l_G_E;
442  l_G_I = (l_G_I > GMAX) ? GMAX : l_G_I;
443  l_G_IB = (l_G_IB > GMAX) ? GMAX : l_G_IB;
444 
445  tauInf = (dt / tau) * (1.0f + l_G_E + l_G_I + l_G_IB + l_gapStrength);
446  VmemInf = (Vrest + l_G_E * Vexc + l_G_I * Vinh + l_G_IB * VinhB + l_GSynGap)
447  / (1.0f + l_G_E + l_G_I + l_G_IB + l_gapStrength);
448 
449  l_V = VmemInf + (l_V - VmemInf) * expf(-tauInf);
450 
451  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
452  // End of code unique to original method
453 
454  bool fired_flag = (l_V > l_Vth);
455 
456  l_activ = fired_flag ? 1.0f : 0.0f;
457  l_V = fired_flag ? Vrest : l_V;
458  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
459  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
460 
461  //
462  // These actions must be done outside of kernel
463  // 1. set activity to 0 in boundary (if needed)
464  // 2. update active indices
465  //
466 
467  // store local variables back to global memory
468  //
469  rnd[k] = l_rnd;
470 
471  activity[kex] = l_activ;
472 
473  V[k] = l_V;
474  Vth[k] = l_Vth;
475 
476  G_E[k] = l_G_E; // G_E_final;
477  G_I[k] = l_G_I; // G_I_final;
478  G_IB[k] = l_G_IB; // G_IB_final;
479  // gapStrength[k] doesn't change;
480 
481  // We blank GSyn here in original, but not in beginning or arma. Why?
482  GSynExc[k] = 0.0f;
483  GSynInh[k] = 0.0f;
484  GSynInhB[k] = 0.0f;
485  GSynGap[k] = 0.0f;
486 
487  } // loop over k
488 }
489 
490 void LIFGap_update_state_beginning(
491  const int nbatch,
492  const int numNeurons,
493  const float time,
494  const float dt,
495 
496  const int nx,
497  const int ny,
498  const int nf,
499  const int lt,
500  const int rt,
501  const int dn,
502  const int up,
503 
504  LIF_params *params,
505  taus_uint4 *rnd,
506  float *V,
507  float *Vth,
508  float *G_E,
509  float *G_I,
510  float *G_IB,
511  float *GSynHead,
512  float *activity,
513 
514  const float *gapStrength) {
515  int k;
516 
517  const float exp_tauE = expf(-dt / params->tauE);
518  const float exp_tauI = expf(-dt / params->tauI);
519  const float exp_tauIB = expf(-dt / params->tauIB);
520  const float exp_tauVth = expf(-dt / params->tauVth);
521 
522  const float dt_sec = 0.001f * dt; // convert to seconds
523 
524  for (k = 0; k < nx * ny * nf * nbatch; k++) {
525  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
526 
527  //
528  // kernel (nonheader part) begins here
529  //
530 
531  // local param variables
532  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
533 
534  // local variables
535  float l_activ;
536 
537  taus_uint4 l_rnd = rnd[k];
538 
539  float l_V = V[k];
540  float l_Vth = Vth[k];
541 
542  // The correction factors to the conductances are so that if l_GSyn_* is the same every
543  // timestep,
544  // then the asymptotic value of l_G_* will be l_GSyn_*
545  float l_G_E = G_E[k];
546  float l_G_I = G_I[k];
547  float l_G_IB = G_IB[k];
548  float l_gapStrength = gapStrength[k];
549 
550  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
551  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
552  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
553  float *GSynGap = &GSynHead[CHANNEL_GAP * nbatch * numNeurons];
554  float l_GSynExc = GSynExc[k];
555  float l_GSynInh = GSynInh[k];
556  float l_GSynInhB = GSynInhB[k];
557  float l_GSynGap = GSynGap[k];
558 
559  // define local param variables
560  //
561  tau = params->tau;
562  Vexc = params->Vexc;
563  Vinh = params->Vinh;
564  VinhB = params->VinhB;
565  Vrest = params->Vrest;
566 
567  VthRest = params->VthRest;
568  deltaVth = params->deltaVth;
569  deltaGIB = params->deltaGIB;
570 
571  // add noise
572  //
573 
574  l_rnd = cl_random_get(l_rnd);
575  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
576  l_rnd = cl_random_get(l_rnd);
577  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
578  }
579 
580  l_rnd = cl_random_get(l_rnd);
581  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
582  l_rnd = cl_random_get(l_rnd);
583  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
584  }
585 
586  l_rnd = cl_random_get(l_rnd);
587  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
588  l_rnd = cl_random_get(l_rnd);
589  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
590  }
591 
592  const float GMAX = 10.0f;
593 
594  // The portion of code below uses the newer method of calculating l_V.
595  float G_E_initial, G_I_initial, G_IB_initial, G_E_final, G_I_final, G_IB_final;
596  float dV1, dV2, dV;
597 
598  G_E_initial = l_G_E + l_GSynExc;
599  G_I_initial = l_G_I + l_GSynInh;
600  G_IB_initial = l_G_IB + l_GSynInhB;
601 
602  G_E_initial = (G_E_initial > GMAX) ? GMAX : G_E_initial;
603  G_I_initial = (G_I_initial > GMAX) ? GMAX : G_I_initial;
604  G_IB_initial = (G_IB_initial > GMAX) ? GMAX : G_IB_initial;
605 
606  G_E_final = G_E_initial * exp_tauE;
607  G_I_final = G_I_initial * exp_tauI;
608  G_IB_final = G_IB_initial * exp_tauIB;
609 
610  dV1 = LIFGap_Vmem_derivative(
611  l_V,
612  G_E_initial,
613  G_I_initial,
614  G_IB_initial,
615  l_GSynGap,
616  Vexc,
617  Vinh,
618  VinhB,
619  l_gapStrength,
620  Vrest,
621  tau);
622  dV2 = LIFGap_Vmem_derivative(
623  l_V + dt * dV1,
624  G_E_final,
625  G_I_final,
626  G_IB_final,
627  l_GSynGap,
628  Vexc,
629  Vinh,
630  VinhB,
631  l_gapStrength,
632  Vrest,
633  tau);
634  dV = (dV1 + dV2) * 0.5f;
635  l_V = l_V + dt * dV;
636 
637  l_G_E = G_E_final;
638  l_G_I = G_I_final;
639  l_G_IB = G_IB_final;
640 
641  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
642  // End of code unique to newer method.
643 
644  bool fired_flag = (l_V > l_Vth);
645 
646  l_activ = fired_flag ? 1.0f : 0.0f;
647  l_V = fired_flag ? Vrest : l_V;
648  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
649  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
650 
651  //
652  // These actions must be done outside of kernel
653  // 1. set activity to 0 in boundary (if needed)
654  // 2. update active indices
655  //
656 
657  // store local variables back to global memory
658  //
659  rnd[k] = l_rnd;
660 
661  activity[kex] = l_activ;
662 
663  V[k] = l_V;
664  Vth[k] = l_Vth;
665 
666  G_E[k] = l_G_E; // G_E_final;
667  G_I[k] = l_G_I; // G_I_final;
668  G_IB[k] = l_G_IB; // G_IB_final;
669  // gapStrength[k] doesn't change
670 
671  // We blank GSyn here in original, but not in beginning or arma. Why?
672 
673  } // loop over k
674 }
675 
676 void LIFGap_update_state_arma(
677  const int nbatch,
678  const int numNeurons,
679  const float time,
680  const float dt,
681 
682  const int nx,
683  const int ny,
684  const int nf,
685  const int lt,
686  const int rt,
687  const int dn,
688  const int up,
689 
690  LIF_params *params,
691  taus_uint4 *rnd,
692  float *V,
693  float *Vth,
694  float *G_E,
695  float *G_I,
696  float *G_IB,
697  float *GSynHead,
698  float *activity,
699 
700  const float *gapStrength) {
701  int k;
702 
703  const float exp_tauE = expf(-dt / params->tauE);
704  const float exp_tauI = expf(-dt / params->tauI);
705  const float exp_tauIB = expf(-dt / params->tauIB);
706  const float exp_tauVth = expf(-dt / params->tauVth);
707 
708  const float dt_sec = 0.001f * dt; // convert to seconds
709 
710  for (k = 0; k < nx * ny * nf * nbatch; k++) {
711  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
712 
713  //
714  // kernel (nonheader part) begins here
715  //
716 
717  // local param variables
718  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
719 
720  const float GMAX = 10.0f;
721 
722  // local variables
723  float l_activ;
724 
725  taus_uint4 l_rnd = rnd[k];
726 
727  float l_V = V[k];
728  float l_Vth = Vth[k];
729 
730  // The correction factors to the conductances are so that if l_GSyn_* is the same every
731  // timestep,
732  // then the asymptotic value of l_G_* will be l_GSyn_*
733  float l_G_E = G_E[k];
734  float l_G_I = G_I[k];
735  float l_G_IB = G_IB[k];
736  float l_gapStrength = gapStrength[k];
737 
738  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
739  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
740  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
741  float *GSynGap = &GSynHead[CHANNEL_GAP * nbatch * numNeurons];
742  float l_GSynExc = GSynExc[k];
743  float l_GSynInh = GSynInh[k];
744  float l_GSynInhB = GSynInhB[k];
745  float l_GSynGap = GSynGap[k];
746 
747  //
748  // start of LIF2_update_exact_linear
749  //
750 
751  // define local param variables
752  //
753  tau = params->tau;
754  Vexc = params->Vexc;
755  Vinh = params->Vinh;
756  VinhB = params->VinhB;
757  Vrest = params->Vrest;
758 
759  VthRest = params->VthRest;
760  deltaVth = params->deltaVth;
761  deltaGIB = params->deltaGIB;
762 
763  // add noise
764  //
765 
766  l_rnd = cl_random_get(l_rnd);
767  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
768  l_rnd = cl_random_get(l_rnd);
769  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
770  }
771 
772  l_rnd = cl_random_get(l_rnd);
773  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
774  l_rnd = cl_random_get(l_rnd);
775  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
776  }
777 
778  l_rnd = cl_random_get(l_rnd);
779  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
780  l_rnd = cl_random_get(l_rnd);
781  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
782  }
783 
784  // The portion of code below uses the newer method of calculating l_V.
785  float G_E_initial, G_I_initial, G_IB_initial, G_E_final, G_I_final, G_IB_final;
786  float tau_inf_initial, tau_inf_final, V_inf_initial, V_inf_final;
787 
788  G_E_initial = l_G_E + l_GSynExc;
789  G_I_initial = l_G_I + l_GSynInh;
790  G_IB_initial = l_G_IB + l_GSynInhB;
791  tau_inf_initial = tau / (1.0f + G_E_initial + G_I_initial + G_IB_initial + l_gapStrength);
792  V_inf_initial =
793  (Vrest + Vexc * G_E_initial + Vinh * G_I_initial + VinhB * G_IB_initial + l_GSynGap)
794  / (1.0f + G_E_initial + G_I_initial + G_IB_initial + l_gapStrength);
795 
796  G_E_initial = (G_E_initial > GMAX) ? GMAX : G_E_initial;
797  G_I_initial = (G_I_initial > GMAX) ? GMAX : G_I_initial;
798  G_IB_initial = (G_IB_initial > GMAX) ? GMAX : G_IB_initial;
799 
800  G_E_final = G_E_initial * exp_tauE;
801  G_I_final = G_I_initial * exp_tauI;
802  G_IB_final = G_IB_initial * exp_tauIB;
803  tau_inf_final = tau / (1.0f + G_E_final + G_I_final + G_IB_final + l_gapStrength);
804  V_inf_final = (Vrest + Vexc * G_E_final + Vinh * G_I_final + VinhB * G_IB_final + l_GSynGap)
805  / (1.0f + G_E_final + G_I_final + G_IB_final + l_gapStrength);
806 
807  float tau_slope = (tau_inf_final - tau_inf_initial) / dt;
808  float f1 = tau_slope == 0.0f ? expf(-dt / tau_inf_initial)
809  : powf(tau_inf_final / tau_inf_initial, -1 / tau_slope);
810  float f2 = tau_slope == -1.0f
811  ? tau_inf_initial / dt * logf(tau_inf_final / tau_inf_initial + 1.0f)
812  : (1 - tau_inf_initial / dt * (1 - f1)) / (1 + tau_slope);
813  float f3 = 1.0f - f1 - f2;
814  l_V = f1 * l_V + f2 * V_inf_initial + f3 * V_inf_final;
815 
816  l_G_E = G_E_final;
817  l_G_I = G_I_final;
818  l_G_IB = G_IB_final;
819 
820  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
821  // End of code unique to newer method.
822 
823  //
824  // start of update_f
825  //
826 
827  bool fired_flag = (l_V > l_Vth);
828 
829  l_activ = fired_flag ? 1.0f : 0.0f;
830  l_V = fired_flag ? Vrest : l_V;
831  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
832  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
833 
834  //
835  // These actions must be done outside of kernel
836  // 1. set activity to 0 in boundary (if needed)
837  // 2. update active indices
838  //
839 
840  // store local variables back to global memory
841  //
842  rnd[k] = l_rnd;
843 
844  activity[kex] = l_activ;
845 
846  V[k] = l_V;
847  Vth[k] = l_Vth;
848 
849  G_E[k] = l_G_E;
850  G_I[k] = l_G_I;
851  G_IB[k] = l_G_IB;
852  // gapStrength[k] doesn't change
853 
854  // We blank GSyn here in original, but not in beginning or arma. Why?
855  }
856 }
static bool completed(Status &a)
Definition: Response.hpp:49