1 files changed, 164 insertions, 0 deletions
diff --git a/code/core/neuron.cpp b/code/core/neuron.cpp
new file mode 100644
index 0000000..c35afcc
--- /dev/null
+++ b/code/core/neuron.cpp
@@ -0,0 +1,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]);
+}