summaryrefslogtreecommitdiff
path: root/code/core/neuron.cpp
blob: c35afcc7d68dddd598b67c1543a1b29af9bec07c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
//#include "model_switch.h"

#include <stdio.h>

#include "neuron.h"

Neuron::Neuron() {
  // neurons are initialized before the Global objects is; therefore this call to init() is preliminary
  // and has to be repeated after the Global object has been loaded
  init();
}

Neuron::init() {
  // general neuron params
  voltage = global.base_voltage;
  refractoryTime = 0.0;

  // IP
  fac_voltage_tau = 1.0;
  fax_current = 1.0;

  ip_R = 1.0;
  ip_C = 1.0;

  ip_est_mom1 = global.ip_dst_mom1;
  ip_est_mom2 = global.ip_dst_mom2;
  ip_dst_mom1 = global.ip_dst_mom1;
  ip_dst_mom2 = global.ip_dst_mom2;

  // clock
  lastEvent = 0.0;
  lastSpike = - INFINITY;
 
  // WARN: synapse lists not initialized (!) (well .. they have a constructor)
}

// evolve the neuron's state time [seconds]
//   RETURN: the time difference between the this and the last update
double Neuron::evolve(double time) {
  // update internal clock
  double dt = time - lastEvent;
  lastEvent = time;

  // voltage decay
  voltage -= (voltage - global.voltage_base) * (1.0 - exp( - td / (global.voltage_tau * fac_voltage_tau)));
  
  return dt;
}

// apply an incoming spike of current to this neurons state giving it the current time [s] of the simulation; returns the time of the next spike to occur
double Neuron::processCurrent(double time, double current) {
  // process the model until the current time
  evolve(time);

  // add spike
  if (time > refractoryTime)
    voltage += fac_current * current;

  // return when the neuron should fire the next time (0.0 is an option, too)
  return predictSpike(time);
}

// generate a spike wich has been predicted earlier to occur at time [s]
// RETURN when the next spike is expected to occur and the current of the spike
double Neuron::generateSpike(double time, bool &doesOccur) {
  // update the model (to the after-firing-state)
  evolve(time);

  // check if a spike occurs (the spike date might be outdated)
  doesOccur = (voltage >= global.fire_threshold) && (time >= refractoryTime);

  if (doesOccur) {
    // 0. update internal state (instantaneously)
    voltage = global.voltage_base;
    refractoryTime = time + global.absoluteRefractoryTime;

    // 1. do intrinsic plasticity
    intrinsicPlasticity(time - lastSpike);

    // 2. do stdp related accounting in the incoming synapses
    //    (check each synapse for event(s))
    for (SynapseSrcList::iterator i = sin.begin(); i != sin.end(); i++) {
      // do not touch constant synapses (calc'ing the values below is a waste of time)
      if (i->constant)
	continue;

      // check if pre-neuron has fired since we have fired the last time
      // this info is stored in the synapse between pre- and this neuron
      if ((lastSpike != -INFINITY) && (i->firstSpike != -INFINITY)) {
	// depression
	i->evolve(i->firstSpike);
	i->stdp(lastSpike - i->firstSpike);

	// facilitation
	i->evolve(time);
	i->stdp(time - i->lastSpike);

	// reset firstFired to allow the next max calculation
	i->firstSpike = -INFINITY;
      }
    }

    // 3. update the last time this neuron has fired
    lastSpike = time;
  }
  
  return predictSpike(time);
}

void Neuron::predictSpike(double time) {
  if (voltage >= global.voltage_threshold) {
    return fmax(time, refractoryTime);
  }else{
    return INFINITY;
  }
}

void Neuron::intrinsicPlasticity(double td) {
  if (!ip)
    return;

  // update rate estimation
  double edt = exp( td / g.ip_sliding_avg_tau );
  double freq = 1.0 / td; // HINT: td=0 is not possible (because of an absolute refractory period)
  ip_est_mom1 = edt * ip_est_mom1 + (1 - edt) * freq;
  ip_est_mom2 = edt * ip_est_mom2 + (1 - edt) * freq * freq;

  // modify internal representation
  ip_R += g.ip_lambda_R * (ip_est_mom1 - ip_dst_mom1) * td;

  double m2d = g.ip_lambda_C * (ip_est_mom2 - ip_dst_mom2);
  ip_C += pow(ip_C, 2) * m2d * td
        / ( 1.0 + ip_C * m2d * td );

  // update coefficients
  fac_voltage_tau = ip_R * ip_C;
  fac_current = ip_R;
}

template<class Action>
bool Neuron::reflect<Action>(Action &a) {
  return 
       a(voltage, "voltage")
    && a(refractoryTime, "refractoryTime")
    && a(ip, "ip")
    && a(ip_est_mom1, "ip_est_mom1")
    && a(ip_est_mom2, "ip_est_mom2")
    && a(ip_dst_mom1, "ip_dst_mom1")
    && a(ip_dst_mom2, "ip_dst_mom2")
    && a(ip_R, "ip_R")
    && a(ip_C, "ip_C")
    && a(fac_voltage_tau, "fac_voltage_tau")
    && a(fac_current, "fac_current")
    && a(lastEvent, "lastEvent")
    && a(lastSpike, "lastSpike");
}

static id_type Neuron::numElements() {
  return s.numNeurons;
}

static Neuron * Neuron::singleton(id_type id) {
  return &(s.neurons[id]);
}
contact: Jan Huwald // Impressum