PetaVision  Alpha
LIF.cpp
1 /*
2  * LIF.cpp
3  *
4  * Created on: Dec 21, 2010
5  * Author: Craig Rasmussen
6  */
7 
8 #include "LIF.hpp"
9 #include "HyPerLayer.hpp"
10 
11 #include "checkpointing/CheckpointEntryRandState.hpp"
12 #include "include/default_params.h"
13 #include "include/pv_common.h"
14 #include "io/fileio.hpp"
15 #include "io/randomstateio.hpp"
16 
17 #include <cassert>
18 #include <cfloat>
19 #include <cmath>
20 #include <cstdio>
21 #include <cstdlib>
22 #include <cstring>
23 #include <memory>
24 
25 void LIF_update_state_arma(
26  const int nbatch,
27  const int numNeurons,
28  const float time,
29  const float dt,
30 
31  const int nx,
32  const int ny,
33  const int nf,
34  const int lt,
35  const int rt,
36  const int dn,
37  const int up,
38 
39  LIF_params *params,
40  taus_uint4 *rnd,
41 
42  float *V,
43  float *Vth,
44  float *G_E,
45  float *G_I,
46  float *G_IB,
47  float *GSynHead,
48  float *activity);
49 
50 void LIF_update_state_beginning(
51  const int nbatch,
52  const int numNeurons,
53  const float time,
54  const float dt,
55 
56  const int nx,
57  const int ny,
58  const int nf,
59  const int lt,
60  const int rt,
61  const int dn,
62  const int up,
63 
64  LIF_params *params,
65  taus_uint4 *rnd,
66 
67  float *V,
68  float *Vth,
69  float *G_E,
70  float *G_I,
71  float *G_IB,
72  float *GSynHead,
73  float *activity);
74 
75 void LIF_update_state_original(
76  const int nbatch,
77  const int numNeurons,
78  const float time,
79  const float dt,
80 
81  const int nx,
82  const int ny,
83  const int nf,
84  const int lt,
85  const int rt,
86  const int dn,
87  const int up,
88 
89  LIF_params *params,
90  taus_uint4 *rnd,
91 
92  float *V,
93  float *Vth,
94  float *G_E,
95  float *G_I,
96  float *G_IB,
97  float *GSynHead,
98  float *activity);
99 
100 namespace PV {
101 
102 LIF::LIF() { initialize_base(); }
103 
104 LIF::LIF(const char *name, HyPerCol *hc) {
105  initialize_base();
106  initialize(name, hc, "LIF_update_state");
107 }
108 
109 LIF::~LIF() {
110  if (numChannels > 0) {
111  // conductances allocated contiguously so this frees all
112  free(G_E);
113  }
114  free(Vth);
115  delete randState;
116  free(methodString);
117 }
118 
119 int LIF::initialize_base() {
120  numChannels = 3;
121  randState = NULL;
122  Vth = NULL;
123  G_E = NULL;
124  G_I = NULL;
125  G_IB = NULL;
126  methodString = NULL;
127 
128  return PV_SUCCESS;
129 }
130 
131 int LIF::initialize(const char *name, HyPerCol *hc, const char *kernel_name) {
132  HyPerLayer::initialize(name, hc);
133  return PV_SUCCESS;
134 }
135 
136 int LIF::ioParamsFillGroup(enum ParamsIOFlag ioFlag) {
137  HyPerLayer::ioParamsFillGroup(ioFlag);
138 
139  ioParam_Vrest(ioFlag);
140  ioParam_Vexc(ioFlag);
141  ioParam_Vinh(ioFlag);
142  ioParam_VinhB(ioFlag);
143  ioParam_VthRest(ioFlag);
144  ioParam_tau(ioFlag);
145  ioParam_tauE(ioFlag);
146  ioParam_tauI(ioFlag);
147  ioParam_tauIB(ioFlag);
148  ioParam_tauVth(ioFlag);
149  ioParam_deltaVth(ioFlag);
150  ioParam_deltaGIB(ioFlag);
151 
152  // NOTE: in LIFDefaultParams, noise ampE, ampI, ampIB were
153  // ampE=0*NOISE_AMP*( 1.0/TAU_EXC )
154  // *(( TAU_INH * (V_REST-V_INH) + TAU_INHB * (V_REST-V_INHB) ) / (V_EXC-V_REST))
155  // ampI=0*NOISE_AMP*1.0
156  // ampIB=0*NOISE_AMP*1.0
157  //
158 
159  ioParam_noiseAmpE(ioFlag);
160  ioParam_noiseAmpI(ioFlag);
161  ioParam_noiseAmpIB(ioFlag);
162  ioParam_noiseFreqE(ioFlag);
163  ioParam_noiseFreqI(ioFlag);
164  ioParam_noiseFreqIB(ioFlag);
165 
166  ioParam_method(ioFlag);
167  return 0;
168 }
169 void LIF::ioParam_Vrest(enum ParamsIOFlag ioFlag) {
170  parent->parameters()->ioParamValue(ioFlag, name, "Vrest", &lParams.Vrest, (float)V_REST);
171 }
172 void LIF::ioParam_Vexc(enum ParamsIOFlag ioFlag) {
173  parent->parameters()->ioParamValue(ioFlag, name, "Vexc", &lParams.Vexc, (float)V_EXC);
174 }
175 void LIF::ioParam_Vinh(enum ParamsIOFlag ioFlag) {
176  parent->parameters()->ioParamValue(ioFlag, name, "Vinh", &lParams.Vinh, (float)V_INH);
177 }
178 void LIF::ioParam_VinhB(enum ParamsIOFlag ioFlag) {
179  parent->parameters()->ioParamValue(ioFlag, name, "VinhB", &lParams.VinhB, (float)V_INHB);
180 }
181 void LIF::ioParam_VthRest(enum ParamsIOFlag ioFlag) {
182  parent->parameters()->ioParamValue(ioFlag, name, "VthRest", &lParams.VthRest, (float)VTH_REST);
183 }
184 void LIF::ioParam_tau(enum ParamsIOFlag ioFlag) {
185  parent->parameters()->ioParamValue(ioFlag, name, "tau", &lParams.tau, (float)TAU_VMEM);
186 }
187 void LIF::ioParam_tauE(enum ParamsIOFlag ioFlag) {
188  parent->parameters()->ioParamValue(ioFlag, name, "tauE", &lParams.tauE, (float)TAU_EXC);
189 }
190 void LIF::ioParam_tauI(enum ParamsIOFlag ioFlag) {
191  parent->parameters()->ioParamValue(ioFlag, name, "tauI", &lParams.tauI, (float)TAU_INH);
192 }
193 void LIF::ioParam_tauIB(enum ParamsIOFlag ioFlag) {
194  parent->parameters()->ioParamValue(ioFlag, name, "tauIB", &lParams.tauIB, (float)TAU_INHB);
195 }
196 void LIF::ioParam_tauVth(enum ParamsIOFlag ioFlag) {
197  parent->parameters()->ioParamValue(ioFlag, name, "tauVth", &lParams.tauVth, (float)TAU_VTH);
198 }
199 void LIF::ioParam_deltaVth(enum ParamsIOFlag ioFlag) {
200  parent->parameters()->ioParamValue(
201  ioFlag, name, "deltaVth", &lParams.deltaVth, (float)DELTA_VTH);
202 }
203 void LIF::ioParam_deltaGIB(enum ParamsIOFlag ioFlag) {
204  parent->parameters()->ioParamValue(
205  ioFlag, name, "deltaGIB", &lParams.deltaGIB, (float)DELTA_G_INHB);
206 }
207 void LIF::ioParam_noiseAmpE(enum ParamsIOFlag ioFlag) {
208  parent->parameters()->ioParamValue(ioFlag, name, "noiseAmpE", &lParams.noiseAmpE, 0.0f);
209 }
210 void LIF::ioParam_noiseAmpI(enum ParamsIOFlag ioFlag) {
211  parent->parameters()->ioParamValue(ioFlag, name, "noiseAmpI", &lParams.noiseAmpI, 0.0f);
212 }
213 void LIF::ioParam_noiseAmpIB(enum ParamsIOFlag ioFlag) {
214  parent->parameters()->ioParamValue(ioFlag, name, "noiseAmpIB", &lParams.noiseAmpIB, 0.0f);
215 }
216 
217 void LIF::ioParam_noiseFreqE(enum ParamsIOFlag ioFlag) {
218  parent->parameters()->ioParamValue(ioFlag, name, "noiseFreqE", &lParams.noiseFreqE, 250.0f);
219  if (ioFlag == PARAMS_IO_READ) {
220  float dt_sec = 0.001f * (float)parent->getDeltaTime(); // seconds
221  if (dt_sec * lParams.noiseFreqE > 1.0f) {
222  lParams.noiseFreqE = 1.0f / dt_sec;
223  }
224  }
225 }
226 
227 void LIF::ioParam_noiseFreqI(enum ParamsIOFlag ioFlag) {
228  parent->parameters()->ioParamValue(ioFlag, name, "noiseFreqI", &lParams.noiseFreqI, 250.0f);
229  if (ioFlag == PARAMS_IO_READ) {
230  float dt_sec = 0.001f * (float)parent->getDeltaTime(); // seconds
231  if (dt_sec * lParams.noiseFreqI > 1.0f) {
232  lParams.noiseFreqI = 1.0f / dt_sec;
233  }
234  }
235 }
236 
237 void LIF::ioParam_noiseFreqIB(enum ParamsIOFlag ioFlag) {
238  parent->parameters()->ioParamValue(ioFlag, name, "noiseFreqIB", &lParams.noiseFreqIB, 250.0f);
239  if (ioFlag == PARAMS_IO_READ) {
240  float dt_sec = 0.001f * (float)parent->getDeltaTime(); // seconds
241  if (dt_sec * lParams.noiseFreqIB > 1.0f) {
242  lParams.noiseFreqIB = 1.0f / dt_sec;
243  }
244  }
245 }
246 
247 void LIF::ioParam_method(enum ParamsIOFlag ioFlag) {
248  // Read the integration method: one of 'arma' (preferred), 'beginning' (deprecated), or
249  // 'original' (deprecated).
250  const char *default_method = "arma";
251  parent->parameters()->ioParamString(
252  ioFlag, name, "method", &methodString, default_method, true /*warnIfAbsent*/);
253  if (ioFlag != PARAMS_IO_READ) {
254  return;
255  }
256 
257  assert(methodString);
258  if (methodString[0] == '\0') {
259  free(methodString);
260  methodString = strdup(default_method);
261  if (methodString == NULL) {
262  Fatal().printf(
263  "%s: unable to set method string: %s\n", getDescription_c(), strerror(errno));
264  }
265  }
266  method = methodString ? methodString[0]
267  : 'a'; // Default is ARMA; 'beginning' and 'original' are deprecated.
268  if (method != 'o' && method != 'b' && method != 'a') {
269  if (parent->columnId() == 0) {
270  ErrorLog().printf(
271  "LIF::setLIFParams error. Layer \"%s\" has method \"%s\". Allowable values are "
272  "\"arma\", \"beginning\" and \"original\".",
273  name,
274  methodString);
275  }
276  MPI_Barrier(parent->getCommunicator()->communicator());
277  exit(EXIT_FAILURE);
278  }
279  if (method != 'a') {
280  if (parent->columnId() == 0) {
281  WarnLog().printf(
282  "LIF layer \"%s\" integration method \"%s\" is deprecated. Method \"arma\" is "
283  "preferred.\n",
284  name,
285  methodString);
286  }
287  }
288 }
289 
290 int LIF::setActivity() {
291  float *activity = clayer->activity->data;
292  memset(activity, 0, sizeof(float) * clayer->numExtendedAllBatches);
293  return 0;
294 }
295 
296 Response::Status
297 LIF::communicateInitInfo(std::shared_ptr<CommunicateInitInfoMessage const> message) {
298  return HyPerLayer::communicateInitInfo(message);
299 }
300 
301 Response::Status LIF::allocateDataStructures() {
302  auto status = HyPerLayer::allocateDataStructures();
303  if (!Response::completed(status)) {
304  return status;
305  }
306 
307  // // a random state variable is needed for every neuron/clthread
308  randState = new Random(getLayerLoc(), false /*isExtended*/);
309  if (randState == nullptr) {
310  Fatal().printf(
311  "LIF::initialize: %s unable to create object of Random class.\n", getDescription_c());
312  }
313 
314  int numNeurons = getNumNeuronsAllBatches();
315  assert(Vth); // Allocated when HyPerLayer::allocateDataStructures() called allocateBuffers().
316  for (size_t k = 0; k < numNeurons; k++) {
317  Vth[k] = lParams.VthRest; // lParams.VthRest is set in setLIFParams
318  }
319  return Response::SUCCESS;
320 }
321 
322 void LIF::allocateBuffers() {
323  allocateConductances(numChannels);
324  Vth = (float *)calloc((size_t)getNumNeuronsAllBatches(), sizeof(float));
325  if (Vth == NULL) {
326  Fatal().printf(
327  "%s: rank %d process unable to allocate memory for Vth: %s\n",
328  getDescription_c(),
329  parent->columnId(),
330  strerror(errno));
331  }
332  HyPerLayer::allocateBuffers();
333 }
334 
335 void LIF::allocateConductances(int num_channels) {
336  assert(num_channels >= 3); // Need exc, inh, and inhb at a minimum.
337  const int numNeurons = getNumNeuronsAllBatches();
338  G_E = (float *)calloc((size_t)(getNumNeuronsAllBatches() * numChannels), sizeof(float));
339  if (G_E == NULL) {
340  Fatal().printf(
341  "%s: rank %d process unable to allocate memory for %d conductances: %s\n",
342  getDescription_c(),
343  parent->columnId(),
344  num_channels,
345  strerror(errno));
346  }
347 
348  G_I = G_E + 1 * numNeurons;
349  G_IB = G_E + 2 * numNeurons;
350 }
351 
352 Response::Status LIF::readStateFromCheckpoint(Checkpointer *checkpointer) {
353  if (initializeFromCheckpointFlag) {
354  auto status = HyPerLayer::readStateFromCheckpoint(checkpointer);
355  if (!Response::completed(status)) {
356  return status;
357  }
358  readVthFromCheckpoint(checkpointer);
359  readG_EFromCheckpoint(checkpointer);
360  readG_IFromCheckpoint(checkpointer);
361  readG_IBFromCheckpoint(checkpointer);
362  readRandStateFromCheckpoint(checkpointer);
363  return Response::SUCCESS;
364  }
365  else {
366  return Response::NO_ACTION;
367  }
368 }
369 
370 void LIF::readVthFromCheckpoint(Checkpointer *checkpointer) {
371  checkpointer->readNamedCheckpointEntry(std::string(name), "Vth", false /*not constant*/);
372 }
373 
374 void LIF::readG_EFromCheckpoint(Checkpointer *checkpointer) {
375  checkpointer->readNamedCheckpointEntry(std::string(name), "G_E", false /*not constant*/);
376 }
377 
378 void LIF::readG_IFromCheckpoint(Checkpointer *checkpointer) {
379  checkpointer->readNamedCheckpointEntry(std::string(name), "G_I", false /*not constant*/);
380 }
381 
382 void LIF::readG_IBFromCheckpoint(Checkpointer *checkpointer) {
383  checkpointer->readNamedCheckpointEntry(std::string(name), "G_IB", false /*not constant*/);
384 }
385 
386 void LIF::readRandStateFromCheckpoint(Checkpointer *checkpointer) {
387  checkpointer->readNamedCheckpointEntry(std::string(name), "rand_state", false /*not constant*/);
388 }
389 
390 Response::Status LIF::registerData(Checkpointer *checkpointer) {
391  auto status = HyPerLayer::registerData(checkpointer);
392  if (!Response::completed(status)) {
393  return status;
394  }
395  checkpointPvpActivityFloat(checkpointer, "Vth", Vth, false /*not extended*/);
396  checkpointPvpActivityFloat(checkpointer, "G_E", G_E, false /*not extended*/);
397  checkpointPvpActivityFloat(checkpointer, "G_I", G_I, false /*not extended*/);
398  checkpointPvpActivityFloat(checkpointer, "G_IB", G_IB, false /*not extended*/);
399  checkpointRandState(checkpointer, "rand_state", randState, false /*not extended*/);
400  return Response::SUCCESS;
401 }
402 
403 Response::Status LIF::updateState(double time, double dt) {
404  update_timer->start();
405 
406  const int nx = clayer->loc.nx;
407  const int ny = clayer->loc.ny;
408  const int nf = clayer->loc.nf;
409  const PVHalo *halo = &clayer->loc.halo;
410  const int nbatch = clayer->loc.nbatch;
411 
412  float *GSynHead = GSyn[0];
413  float *activity = clayer->activity->data;
414 
415  switch (method) {
416  case 'a':
417  LIF_update_state_arma(
418  nbatch,
419  getNumNeurons(),
420  time,
421  dt,
422  nx,
423  ny,
424  nf,
425  halo->lt,
426  halo->rt,
427  halo->dn,
428  halo->up,
429  &lParams,
430  randState->getRNG(0),
431  clayer->V,
432  Vth,
433  G_E,
434  G_I,
435  G_IB,
436  GSynHead,
437  activity);
438  break;
439  case 'b':
440  LIF_update_state_beginning(
441  nbatch,
442  getNumNeurons(),
443  time,
444  dt,
445  nx,
446  ny,
447  nf,
448  halo->lt,
449  halo->rt,
450  halo->dn,
451  halo->up,
452  &lParams,
453  randState->getRNG(0),
454  clayer->V,
455  Vth,
456  G_E,
457  G_I,
458  G_IB,
459  GSynHead,
460  activity);
461  break;
462  case 'o':
463  LIF_update_state_original(
464  nbatch,
465  getNumNeurons(),
466  time,
467  dt,
468  nx,
469  ny,
470  nf,
471  halo->lt,
472  halo->rt,
473  halo->dn,
474  halo->up,
475  &lParams,
476  randState->getRNG(0),
477  clayer->V,
478  Vth,
479  G_E,
480  G_I,
481  G_IB,
482  GSynHead,
483  activity);
484  break;
485  default: assert(0); break;
486  }
487  update_timer->stop();
488  return Response::SUCCESS;
489 }
490 
491 float LIF::getChannelTimeConst(enum ChannelType channel_type) {
492  float channel_time_const = 0.0f;
493  switch (channel_type) {
494  case CHANNEL_EXC: channel_time_const = lParams.tauE; break;
495  case CHANNEL_INH: channel_time_const = lParams.tauI; break;
496  case CHANNEL_INHB: channel_time_const = lParams.tauIB; break;
497  default: channel_time_const = 0.0f; break;
498  }
499  return channel_time_const;
500 }
501 
502 int LIF::findPostSynaptic(
503  int dim,
504  int maxSize,
505  int col,
506  // input: which layer, which neuron
507  HyPerLayer *lSource,
508  float pos[],
509 
510  // output: how many of our neurons are connected.
511  // an array with their indices.
512  // an array with their feature vectors.
513  int *nNeurons,
514  int nConnectedNeurons[],
515  float *vPos) {
516  return 0;
517 }
518 
519 } // namespace PV
520 
522 //
523 // implementation of LIF kernels
524 //
525 
526 inline float LIF_Vmem_derivative(
527  const float Vmem,
528  const float G_E,
529  const float G_I,
530  const float G_IB,
531  const float V_E,
532  const float V_I,
533  const float V_IB,
534  const float Vrest,
535  const float tau) {
536  float totalconductance = 1.0f + G_E + G_I + G_IB;
537  float Vmeminf = (Vrest + V_E * G_E + V_I * G_I + V_IB * G_IB) / totalconductance;
538  return totalconductance * (Vmeminf - Vmem) / tau;
539 }
540 //
541 // update the state of a retinal layer (spiking)
542 //
543 // assume called with 1D kernel
544 //
545 // LIF_update_state_original uses an Euler scheme for V where the conductances over the entire
546 // timestep are taken to be the values calculated at the end of the timestep
547 // LIF_update_state_beginning uses a Heun scheme for V, using values of the conductances at both the
548 // beginning and end of the timestep. Spikes in the input are applied at the beginning of the
549 // timestep.
550 // LIF_update_state_arma uses an auto-regressive moving average filter for V, applying the GSyn at
551 // the start of the timestep and assuming that tau_inf and V_inf vary linearly over the timestep.
552 // See van Hateren, Journal of Vision (2005), p. 331.
553 //
554 void LIF_update_state_original(
555  const int nbatch,
556  const int numNeurons,
557  const float time,
558  const float dt,
559 
560  const int nx,
561  const int ny,
562  const int nf,
563  const int lt,
564  const int rt,
565  const int dn,
566  const int up,
567 
568  LIF_params *params,
569  taus_uint4 *rnd,
570  float *V,
571  float *Vth,
572  float *G_E,
573  float *G_I,
574  float *G_IB,
575  float *GSynHead,
576  float *activity) {
577  int k;
578 
579  const float exp_tauE = expf(-dt / params->tauE);
580  const float exp_tauI = expf(-dt / params->tauI);
581  const float exp_tauIB = expf(-dt / params->tauIB);
582  const float exp_tauVth = expf(-dt / params->tauVth);
583 
584  const float dt_sec = 0.001f * dt; // convert to seconds
585 
586  for (k = 0; k < nx * ny * nf * nbatch; k++) {
587  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
588 
589  //
590  // kernel (nonheader part) begins here
591  //
592 
593  // local param variables
594  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
595 
596  const float GMAX = 10.0f;
597 
598  // local variables
599  float l_activ;
600 
601  taus_uint4 l_rnd = rnd[k];
602 
603  float l_V = V[k];
604  float l_Vth = Vth[k];
605 
606  float l_G_E = G_E[k];
607  float l_G_I = G_I[k];
608  float l_G_IB = G_IB[k];
609 
610  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
611  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
612  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
613  float l_GSynExc = GSynExc[k];
614  float l_GSynInh = GSynInh[k];
615  float l_GSynInhB = GSynInhB[k];
616 
617  // temporary arrays
618  float tauInf, VmemInf;
619 
620  //
621  // start of LIF2_update_exact_linear
622  //
623 
624  // define local param variables
625  //
626  tau = params->tau;
627  Vexc = params->Vexc;
628  Vinh = params->Vinh;
629  VinhB = params->VinhB;
630  Vrest = params->Vrest;
631 
632  VthRest = params->VthRest;
633  deltaVth = params->deltaVth;
634  deltaGIB = params->deltaGIB;
635 
636  // add noise
637  //
638 
639  l_rnd = cl_random_get(l_rnd);
640  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
641  l_rnd = cl_random_get(l_rnd);
642  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
643  }
644 
645  l_rnd = cl_random_get(l_rnd);
646  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
647  l_rnd = cl_random_get(l_rnd);
648  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
649  }
650 
651  l_rnd = cl_random_get(l_rnd);
652  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
653  l_rnd = cl_random_get(l_rnd);
654  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
655  }
656 
657  l_G_E = l_GSynExc + l_G_E * exp_tauE;
658  l_G_I = l_GSynInh + l_G_I * exp_tauI;
659  l_G_IB = l_GSynInhB + l_G_IB * exp_tauIB;
660 
661  l_G_E = (l_G_E > GMAX) ? GMAX : l_G_E;
662  l_G_I = (l_G_I > GMAX) ? GMAX : l_G_I;
663  l_G_IB = (l_G_IB > GMAX) ? GMAX : l_G_IB;
664 
665  tauInf = (dt / tau) * (1.0f + l_G_E + l_G_I + l_G_IB);
666  VmemInf = (Vrest + l_G_E * Vexc + l_G_I * Vinh + l_G_IB * VinhB)
667  / (1.0f + l_G_E + l_G_I + l_G_IB);
668 
669  l_V = VmemInf + (l_V - VmemInf) * expf(-tauInf);
670 
671  //
672  // start of LIF2_update_finish
673  //
674 
675  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
676 
677  //
678  // start of update_f
679  //
680 
681  bool fired_flag = (l_V > l_Vth);
682 
683  l_activ = fired_flag ? 1.0f : 0.0f;
684  l_V = fired_flag ? Vrest : l_V;
685  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
686  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
687 
688  //
689  // These actions must be done outside of kernel
690  // 1. set activity to 0 in boundary (if needed)
691  // 2. update active indices
692  //
693 
694  // store local variables back to global memory
695  //
696  rnd[k] = l_rnd;
697 
698  activity[kex] = l_activ;
699 
700  V[k] = l_V;
701  Vth[k] = l_Vth;
702 
703  G_E[k] = l_G_E;
704  G_I[k] = l_G_I;
705  G_IB[k] = l_G_IB;
706 
707  GSynExc[k] = 0.0f;
708  GSynInh[k] = 0.0f;
709  GSynInhB[k] = 0.0f;
710 
711  } // loop over k
712 }
713 
714 void LIF_update_state_beginning(
715  const int nbatch,
716  const int numNeurons,
717  const float time,
718  const float dt,
719 
720  const int nx,
721  const int ny,
722  const int nf,
723  const int lt,
724  const int rt,
725  const int dn,
726  const int up,
727 
728  LIF_params *params,
729  taus_uint4 *rnd,
730  float *V,
731  float *Vth,
732  float *G_E,
733  float *G_I,
734  float *G_IB,
735  float *GSynHead,
736  // float * GSynExc,
737  // float * GSynInh,
738  // float * GSynInhB,
739  float *activity) {
740  int k;
741 
742  const float exp_tauE = expf(-dt / params->tauE);
743  const float exp_tauI = expf(-dt / params->tauI);
744  const float exp_tauIB = expf(-dt / params->tauIB);
745  const float exp_tauVth = expf(-dt / params->tauVth);
746 
747  const float dt_sec = 0.001f * dt; // convert to seconds
748 
749  for (k = 0; k < nx * ny * nf * nbatch; k++) {
750 
751  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
752 
753  //
754  // kernel (nonheader part) begins here
755  //
756 
757  // local param variables
758  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
759 
760  const float GMAX = 10.0f;
761 
762  // local variables
763  float l_activ;
764 
765  taus_uint4 l_rnd = rnd[k];
766 
767  float l_V = V[k];
768  float l_Vth = Vth[k];
769 
770  // The correction factors to the conductances are so that if l_GSyn_* is the same every
771  // timestep,
772  // then the asymptotic value of l_G_* will be l_GSyn_*
773  float l_G_E = G_E[k];
774  float l_G_I = G_I[k];
775  float l_G_IB = G_IB[k];
776 
777  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
778  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
779  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
780  float l_GSynExc = GSynExc[k];
781  float l_GSynInh = GSynInh[k];
782  float l_GSynInhB = GSynInhB[k];
783 
784  //
785  // start of LIF2_update_exact_linear
786  //
787 
788  // define local param variables
789  //
790  tau = params->tau;
791  Vexc = params->Vexc;
792  Vinh = params->Vinh;
793  VinhB = params->VinhB;
794  Vrest = params->Vrest;
795 
796  VthRest = params->VthRest;
797  deltaVth = params->deltaVth;
798  deltaGIB = params->deltaGIB;
799 
800  // add noise
801  //
802 
803  l_rnd = cl_random_get(l_rnd);
804  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
805  l_rnd = cl_random_get(l_rnd);
806  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
807  }
808 
809  l_rnd = cl_random_get(l_rnd);
810  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
811  l_rnd = cl_random_get(l_rnd);
812  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
813  }
814 
815  l_rnd = cl_random_get(l_rnd);
816  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
817  l_rnd = cl_random_get(l_rnd);
818  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
819  }
820 
821  // The portion of code below uses the newer method of calculating l_V.
822  float G_E_initial, G_I_initial, G_IB_initial, G_E_final, G_I_final, G_IB_final;
823  float dV1, dV2, dV;
824 
825  G_E_initial = l_G_E + l_GSynExc;
826  G_I_initial = l_G_I + l_GSynInh;
827  G_IB_initial = l_G_IB + l_GSynInhB;
828 
829  G_E_initial = (G_E_initial > GMAX) ? GMAX : G_E_initial;
830  G_I_initial = (G_I_initial > GMAX) ? GMAX : G_I_initial;
831  G_IB_initial = (G_IB_initial > GMAX) ? GMAX : G_IB_initial;
832 
833  G_E_final = G_E_initial * exp_tauE;
834  G_I_final = G_I_initial * exp_tauI;
835  G_IB_final = G_IB_initial * exp_tauIB;
836 
837  dV1 = LIF_Vmem_derivative(
838  l_V, G_E_initial, G_I_initial, G_IB_initial, Vexc, Vinh, VinhB, Vrest, tau);
839  dV2 = LIF_Vmem_derivative(
840  l_V + dt * dV1, G_E_final, G_I_final, G_IB_final, Vexc, Vinh, VinhB, Vrest, tau);
841  dV = (dV1 + dV2) * 0.5f;
842  l_V = l_V + dt * dV;
843 
844  l_G_E = G_E_final;
845  l_G_I = G_I_final;
846  l_G_IB = G_IB_final;
847 
848  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
849  // End of code unique to newer method.
850 
851  //
852  // start of update_f
853  //
854 
855  bool fired_flag = (l_V > l_Vth);
856 
857  l_activ = fired_flag ? 1.0f : 0.0f;
858  l_V = fired_flag ? Vrest : l_V;
859  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
860  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
861 
862  //
863  // These actions must be done outside of kernel
864  // 1. set activity to 0 in boundary (if needed)
865  // 2. update active indices
866  //
867 
868  // store local variables back to global memory
869  //
870  rnd[k] = l_rnd;
871 
872  activity[kex] = l_activ;
873 
874  V[k] = l_V;
875  Vth[k] = l_Vth;
876 
877  G_E[k] = l_G_E;
878  G_I[k] = l_G_I;
879  G_IB[k] = l_G_IB;
880 
881  } // loop over k
882 }
883 
884 void LIF_update_state_arma(
885  const int nbatch,
886  const int numNeurons,
887  const float time,
888  const float dt,
889 
890  const int nx,
891  const int ny,
892  const int nf,
893  const int lt,
894  const int rt,
895  const int dn,
896  const int up,
897 
898  LIF_params *params,
899  taus_uint4 *rnd,
900  float *V,
901  float *Vth,
902  float *G_E,
903  float *G_I,
904  float *G_IB,
905  float *GSynHead,
906  float *activity) {
907  int k;
908 
909  const float exp_tauE = expf(-dt / params->tauE);
910  const float exp_tauI = expf(-dt / params->tauI);
911  const float exp_tauIB = expf(-dt / params->tauIB);
912  const float exp_tauVth = expf(-dt / params->tauVth);
913 
914  const float dt_sec = 0.001f * dt; // convert to seconds
915 
916  for (k = 0; k < nx * ny * nf * nbatch; k++) {
917  int kex = kIndexExtendedBatch(k, nbatch, nx, ny, nf, lt, rt, dn, up);
918 
919  //
920  // kernel (nonheader part) begins here
921  //
922 
923  // local param variables
924  float tau, Vrest, VthRest, Vexc, Vinh, VinhB, deltaVth, deltaGIB;
925 
926  const float GMAX = 10.0;
927 
928  // local variables
929  float l_activ;
930 
931  taus_uint4 l_rnd = rnd[k];
932 
933  float l_V = V[k];
934  float l_Vth = Vth[k];
935 
936  // The correction factors to the conductances are so that if l_GSyn_* is the same every
937  // timestep,
938  // then the asymptotic value of l_G_* will be l_GSyn_*
939  float l_G_E = G_E[k];
940  float l_G_I = G_I[k];
941  float l_G_IB = G_IB[k];
942 
943  float *GSynExc = &GSynHead[CHANNEL_EXC * nbatch * numNeurons];
944  float *GSynInh = &GSynHead[CHANNEL_INH * nbatch * numNeurons];
945  float *GSynInhB = &GSynHead[CHANNEL_INHB * nbatch * numNeurons];
946  float l_GSynExc = GSynExc[k];
947  float l_GSynInh = GSynInh[k];
948  float l_GSynInhB = GSynInhB[k];
949 
950  //
951  // start of LIF2_update_exact_linear
952  //
953 
954  // define local param variables
955  //
956  tau = params->tau;
957  Vexc = params->Vexc;
958  Vinh = params->Vinh;
959  VinhB = params->VinhB;
960  Vrest = params->Vrest;
961 
962  VthRest = params->VthRest;
963  deltaVth = params->deltaVth;
964  deltaGIB = params->deltaGIB;
965 
966  // add noise
967  //
968 
969  l_rnd = cl_random_get(l_rnd);
970  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqE) {
971  l_rnd = cl_random_get(l_rnd);
972  l_GSynExc = l_GSynExc + params->noiseAmpE * cl_random_prob(l_rnd);
973  }
974 
975  l_rnd = cl_random_get(l_rnd);
976  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqI) {
977  l_rnd = cl_random_get(l_rnd);
978  l_GSynInh = l_GSynInh + params->noiseAmpI * cl_random_prob(l_rnd);
979  }
980 
981  l_rnd = cl_random_get(l_rnd);
982  if (cl_random_prob(l_rnd) < dt_sec * params->noiseFreqIB) {
983  l_rnd = cl_random_get(l_rnd);
984  l_GSynInhB = l_GSynInhB + params->noiseAmpIB * cl_random_prob(l_rnd);
985  }
986 
987  // The portion of code below uses the newer method of calculating l_V.
988  float G_E_initial, G_I_initial, G_IB_initial, G_E_final, G_I_final, G_IB_final;
989  float tau_inf_initial, tau_inf_final, V_inf_initial, V_inf_final;
990 
991  G_E_initial = l_G_E + l_GSynExc;
992  G_I_initial = l_G_I + l_GSynInh;
993  G_IB_initial = l_G_IB + l_GSynInhB;
994  tau_inf_initial = tau / (1 + G_E_initial + G_I_initial + G_IB_initial);
995  V_inf_initial = (Vrest + Vexc * G_E_initial + Vinh * G_I_initial + VinhB * G_IB_initial)
996  / (1 + G_E_initial + G_I_initial + G_IB_initial);
997 
998  G_E_initial = (G_E_initial > GMAX) ? GMAX : G_E_initial;
999  G_I_initial = (G_I_initial > GMAX) ? GMAX : G_I_initial;
1000  G_IB_initial = (G_IB_initial > GMAX) ? GMAX : G_IB_initial;
1001 
1002  G_E_final = G_E_initial * exp_tauE;
1003  G_I_final = G_I_initial * exp_tauI;
1004  G_IB_final = G_IB_initial * exp_tauIB;
1005 
1006  tau_inf_final = tau / (1 + G_E_final + G_I_final + G_IB_initial);
1007  V_inf_final = (Vrest + Vexc * G_E_final + Vinh * G_I_final + VinhB * G_IB_final)
1008  / (1 + G_E_final + G_I_final + G_IB_final);
1009 
1010  float tau_slope = (tau_inf_final - tau_inf_initial) / dt;
1011  float f1 = tau_slope == 0.0f ? expf(-dt / tau_inf_initial)
1012  : powf(tau_inf_final / tau_inf_initial, -1 / tau_slope);
1013  float f2 = tau_slope == -1.0f
1014  ? tau_inf_initial / dt * logf(tau_inf_final / tau_inf_initial + 1.0f)
1015  : (1 - tau_inf_initial / dt * (1 - f1)) / (1 + tau_slope);
1016  float f3 = 1.0f - f1 - f2;
1017  l_V = f1 * l_V + f2 * V_inf_initial + f3 * V_inf_final;
1018 
1019  l_G_E = G_E_final;
1020  l_G_I = G_I_final;
1021  l_G_IB = G_IB_final;
1022 
1023  l_Vth = VthRest + (l_Vth - VthRest) * exp_tauVth;
1024  // End of code unique to newer method.
1025 
1026  //
1027  // start of update_f
1028  //
1029 
1030  bool fired_flag = (l_V > l_Vth);
1031 
1032  l_activ = fired_flag ? 1.0f : 0.0f;
1033  l_V = fired_flag ? Vrest : l_V;
1034  l_Vth = fired_flag ? l_Vth + deltaVth : l_Vth;
1035  l_G_IB = fired_flag ? l_G_IB + deltaGIB : l_G_IB;
1036 
1037  //
1038  // These actions must be done outside of kernel
1039  // 1. set activity to 0 in boundary (if needed)
1040  // 2. update active indices
1041  //
1042 
1043  // store local variables back to global memory
1044  //
1045  rnd[k] = l_rnd;
1046 
1047  activity[kex] = l_activ;
1048 
1049  V[k] = l_V;
1050  Vth[k] = l_Vth;
1051 
1052  G_E[k] = l_G_E;
1053  G_I[k] = l_G_I;
1054  G_IB[k] = l_G_IB;
1055  }
1056 }
LIF_params
Definition: LIF.hpp:23
PV::HyPerLayer::ioParamsFillGroup
virtual int ioParamsFillGroup(enum ParamsIOFlag ioFlag) override
Definition: HyPerLayer.cpp:571
PV::LIF::ioParamsFillGroup
virtual int ioParamsFillGroup(enum ParamsIOFlag ioFlag) override
Definition: LIF.cpp:136
PVHalo_
Definition: PVLayerLoc.h:14
PV::Response::completed
static bool completed(Status &a)
Definition: Response.hpp:49
PV::HyPerLayer::initialize
int initialize(const char *name, HyPerCol *hc)
Definition: HyPerLayer.cpp:129
taus_uint4_
Definition: pv_types.h:103