Initial commitHEADmaster

30 files changed, 1928 insertions, 0 deletions

diff --git a/code/core/Makefile b/code/core/Makefile
new file mode 100644
index 0000000..df0dd28
--- /dev/null
+++ b/code/core/Makefile
@@ -0,0 +1,57 @@
+# HINT: the paradigma is not to create object files to allow all otpimizations
+#	take place even though the code is scattere across many files
+
+CC=g++
+CCFLAGS=-O3 -ggdb
+LDFLAGS=
+INCLUDE=-I/home/huwald/src/boost
+
+# flags to test
+# -ffast-math -fno-finite-math-only
+
+# debug options
+#   normal program:
+#DEBUGOPTS=
+#   verbose status line for every (!) spike:
+DEBUGOPTS=-DDEBUG_STATUSLINE
+#   enale profiling:
+#DEBUGOPTS=-pg
+
+# list of files every build depends on
+#BASE_SRC_FILES=neuron.cpp simulate.cpp synapse.cpp
+BASE_SRC_FILES=\
+	model_switch.h \
+	reward.h \
+	reward.cpp \
+	event.h \
+	event.cpp \
+	interface.h \
+	interface.cpp \
+	min.h \
+	max.h \
+	neuron.cpp \
+	neuron.h \
+	simulate.cpp \
+	simulate.h \
+	synapse.cpp \
+	synapse.h \
+	tracepoints.h \
+	tracepoints.cpp \
+	fileutils.cpp
+
+.PHONY: all clean wordcount
+
+all: sim print_defaults
+
+clean:
+	rm -f *~ massif.*.* sim print_defaults
+
+
+sim: $(BASE_SRC_FILES) Makefile
+	$(CC) -o $@ $(CCFLAGS) simulate.cpp $(INCLUDE) $(LDFLAGS) $(DEBUGOPTS) -DMODEL_`echo $* | tr '[:lower:]' '[:upper:]'`
+
+print_defaults: $(BASE_SRC_FILES) print_defaults.cpp
+	$(CC) -o $@ $(CCFLAGS) print_defaults.cpp $(LDFLAGS) $(DEBUGOPTS) -DMODEL_`echo $* | tr '[:lower:]' '[:upper:]'`
+
+wordcount:
+	wc *h *cpp Makefile
diff --git a/code/core/bin.cpp b/code/core/bin.cpp
new file mode 100644
index 0000000..e93e4c8
--- /dev/null
+++ b/code/core/bin.cpp
@@ -0,0 +1,6 @@
+#include "bin.h"
+
+void Bin::bin(int neuron) {
+  if (neurons->count(neuron))
+    count++;
+}
diff --git a/code/core/bin.h b/code/core/bin.h
new file mode 100644
index 0000000..302df56
--- /dev/null
+++ b/code/core/bin.h
@@ -0,0 +1,15 @@
+#ifndef BIN_H
+#define BIN_H
+
+#include <set>
+
+class Bin {
+public:
+  Bin(std::set<int> * neurons) : neurons(neurons), count(0) {}
+  void bin(int neuron);
+
+  std::set<int> * neurons;
+  int count;
+};
+
+#endif // BIN_H
diff --git a/code/core/event.cpp b/code/core/event.cpp
new file mode 100644
index 0000000..cb6f5ff
--- /dev/null
+++ b/code/core/event.cpp
@@ -0,0 +1,128 @@
+#include "simulate.h"
+
+#include "event.h"
+
+void Event::execute() {
+  switch (type) {
+  case 0:
+    ((InternalSpike*) this)->execute();
+    break;
+
+  case 1:
+    ((ExternalSpike*) this)->execute();
+    break;
+
+  case 2:
+    ((ExternalNoise*) this)->execute();
+    break;
+
+  case 3:
+    ((IntrinsicNeuron*) this)->execute();
+    break;
+
+  case 4:
+    ((VirtualEvent*) this)->vexecute();
+    break;
+  }
+}
+
+void Event::free() {
+  switch (type) {
+  case 0:
+    delete (InternalSpike*) this;
+    break;
+
+  case 1:
+    delete (ExternalSpike*) this;
+    break;
+
+  case 2:
+    delete (ExternalNoise*) this;
+    break;
+
+  case 3:
+    delete (IntrinsicNeuron*) this;
+    break;
+
+  case 4:
+    delete (VirtualEvent*) this;
+    break;
+  }
+}
+
+// comparison of two events regarding their order in the event list; spike events are processed prior to non-spike events (as spikes are incoming and non-spike are intrinsic and depending on all events up to the present)
+bool Event::operator>(Event &e) {
+  return (time > e.time) || ((time == e.time) && (type > e.type));
+}
+
+void Spike::execute() {
+    double nextEventTime = s.neurons[dst].processCurrent(time, current);
+
+    // check if there is an event to occur
+    if (nextEventTime != INFINITY)
+      s.addEvent(new IntrinsicNeuron(nextEventTime, dst));
+}
+
+void ExternalSpike::execute() {
+  // exec common spike code
+  ((Spike *) this)->execute();
+}
+
+void InternalSpike::execute() {
+  // the spike has to be registered as delieverd in the sending synapse
+  synapse->lastSpike = time;
+  if (synapse->firstSpike == -INFINITY)
+    synapse->firstSpike = time;
+
+  // exec common spike code
+  ((Spike *) this)->execute();
+
+  // synaptic scaling
+  s.neurons[dst].normalizeWeight(g.sumWeight);
+}
+
+void IntrinsicNeuron::execute() {
+  // check if a spike does occur
+  bool spikeDoesOccur;
+  double nextEventTime = s.neurons[dst].generateSpike(time, spikeDoesOccur);
+
+  // check if a spike should be generated
+  if (!spikeDoesOccur)
+    return;
+  
+  // add next intrinsic event
+  if (nextEventTime != INFINITY)
+    s.addEvent(new IntrinsicNeuron(nextEventTime, dst));
+  
+  // create the spike for every neuron where it will be recieved
+  for (SynapseDstList::iterator i = s.neurons[dst].sout.begin(); i != s.neurons[dst].sout.end(); i++) {
+    double ntime = time, ncurrent;
+    i->computePostsynapticPulse(ntime, ncurrent);
+    s.addEvent(new InternalSpike(ntime, i->dst, ncurrent, &(*i)));
+  }
+
+  // check if we should output a spike of this neuron
+  for (s.spikeOIfList.iterator i = s.spikeOIfList.begin(); i != s.spikeOIfList.end(); i++)
+    if (i->isContained(dst)) {
+      SpikeMUX s(dst, time);
+      i->pushObject(s);
+    }
+
+  // bin spikes
+  for (std::set<Bin*>::iterator bin = s.binSets.begin(); bin != s.to.traceBinSets.end(); bin++)
+    (*bin)->bin(dst);
+}
+
+void ExternalNoise::execute() {
+  // pick a random neuron to send noise to (using a temporary event)
+  s.addEvent(new ExternalSpike(time, rand() % (s.numNeurons * 4 / 5), g.en_current));
+
+  // select the next timepoint from a poisson process
+  s.addEvent(new ExternalNoise(time - (log(1.0 - drand48()) / (g.en_freq * s.numNeurons))));
+}
+
+virtual void GlobalUpdate::vexecute() {
+  static const double td = 0.001;
+  g.evolve(td);
+  s.addEvent(new GlobalUpdate(time + td));
+}
diff --git a/code/core/event.h b/code/core/event.h
new file mode 100644
index 0000000..20d348a
--- /dev/null
+++ b/code/core/event.h
@@ -0,0 +1,74 @@
+#ifndef EVENT_H
+#define EVENT_H
+
+class Event {
+public:
+  double time;
+  int type;
+
+  Event(double ti, int ty) : time(ti), type(ty) {}
+  void execute();
+  void free();
+  bool operator>(Event &e1); // comparison of two events regarding their order in the event list; spike events are processed prior to non-spike events (as spikes are incoming and non-spike are intrinsic and depending on all events up to the present)
+
+private:
+  Event() {};
+};
+
+class Spike : public Event {
+public:
+  int dst;
+  double current;
+
+  Spike(double time, int type, int dst, double current) : Event(time, type), dst(dst), current(current) {}
+  void execute();
+};
+
+
+class ExternalSpike : public Spike {
+public:
+  ExternalSpike(double time, int dst, double current) : Spike(time, 1, dst, current) {}
+  void execute();
+};
+
+class InternalSpike : public Spike {
+public:
+  Synapse *synapse;
+
+  InternalSpike(double time, int dst, double current, Synapse *synapse) : Spike(time, 0, dst, current), synapse(synapse) {}
+  void execute();
+};
+
+class IntrinsicNeuron : public Event {
+public:
+  int dst;
+
+  IntrinsicNeuron(double time, int dst) : Event(time, 3), dst(dst) {}
+  void execute();
+};
+
+class ExternalNoise : public Event {
+public:
+  ExternalNoise(double time) : Event(time, 2) {}
+  void execute();
+};
+
+class VirtualEvent : public Event {
+public:
+  VirtualEvent(double time) : Event(time, 4) {}
+  virtual void vexecute() {}
+  virtual ~VirtualEvent() {}
+};
+
+class GlobalUpdate : public VirtualEvent {
+public:
+  GlobalUpdate(double time) : VirtualEvent(time) {}
+  virtual void vexecute();
+};
+
+class PEventGreater {
+public:
+  bool operator() (Event *e1, Event *e2) { return *e1 > *e2; }
+};
+
+#endif // EVENT_H
diff --git a/code/core/fileutils.cpp b/code/core/fileutils.cpp
new file mode 100644
index 0000000..9d080ce
--- /dev/null
+++ b/code/core/fileutils.cpp
@@ -0,0 +1,37 @@
+#include <string.h>
+#include <stdlib.h>
+
+#include "fileutils.h"
+
+// direction is either false=input or true=output and used to overload parameter '-'
+FILE *fd_magic(char *str, bool direction) {
+  static bool usedDir[2] = { false, false };
+
+  // check if stdin/-out is wanted
+  if (strcmp(str, "-") == 0) {
+    // check if we already used this
+    if (usedDir[direction]) {
+      fprintf(stderr, "stdin/-out cannot be used twice\n");
+      exit(-1);
+    }
+
+    // mark that we used it and return in
+    usedDir[direction] = true;
+    return direction ? stdout : stdin;
+  }
+
+  // check if stdout is wanted
+  if (strcmp(str, "0") == 0) {
+    // replace filename and proceed
+    str = "/dev/null";
+  }
+
+  // open file traditionally
+  FILE *fd = fopen(str, direction ? "w" : "r" );
+  if (fd == NULL) {
+    fprintf(stderr, "Failed to open %s\n", str);
+    exit(-1);
+  }
+
+  return fd;
+}
diff --git a/code/core/fileutils.h b/code/core/fileutils.h
new file mode 100644
index 0000000..324dab7
--- /dev/null
+++ b/code/core/fileutils.h
@@ -0,0 +1,10 @@
+#ifndef FILEUTILS_H
+#define FILEUTILS_H
+
+#include <stdio.h>
+
+// direction is either false=input or true=output and used to overload parameter '-'
+FILE *fd_magic(char *str, bool direction);
+
+
+#endif // FILEUTILS_H
diff --git a/code/core/global.cpp b/code/core/global.cpp
new file mode 100644
index 0000000..83ec666
--- /dev/null
+++ b/code/core/global.cpp
@@ -0,0 +1,115 @@
+#include <math.h>
+
+#include "interface.h"
+
+#include "global.h"
+
+Global::Global() :
+    // neuron related
+  voltage_tau(0.05),  // [s]
+  voltage_base(-0.065), // [V] rest potential 
+  threshold(-0.05),    // [V] above which a spike is emitted
+
+  Wmax(0.004), // [V] 
+  Wmin(0.0),   // [V]
+  sumWeight(0.16),    // [V]
+  absoluteRefractoryTime(0.001), // [V]
+
+  // dopamin
+  dopamin_level(0.0), // [?]
+  dopamin_tau(0.005),   // [s]
+
+  // STDP
+  stdp_tau_plus(0.014),     // [s]
+  stdp_tau_minus(0.034),    // [s]
+  stdp_lambda_plus(0.066),  // [1]
+  stdp_lambda_minus(0.066), // [1]
+  stdp_et_tau(0.1),       // [s] eligibility trace tau
+
+  // IP
+  ip_sliding_avg_tau(), // TODO
+  ip_lambda_R(),
+  ip_lambda_C(),
+  ip_dst_mom1(),
+  ip_dst_mom2(),
+
+  // external noise
+  en_current(0.1) // [V]
+  en_freq(1.0) // [Hz/Neuron]
+
+  // trainer
+  trainer_eval_delay(0.001), // [s]
+  trainer_numSymbols(2),            // [1]
+  trainer_refractoryTime(0.2),     // [s]
+  trainer_rewardAmount(0.1),       // [?]
+  trainer_rd_c1(0.05), // readout delay with
+  trainer_rd_c2(0.0001), //   delay = c1 + t*c2 + r*c3 + t*r'*c4
+  trainer_rd_c3(0.0), // where
+  trainer_rd_c4(0.0001); //   r and r' are two distinct random numbers (equ. dist in [0,1))
+{}
+
+void Global::evolve(double td) {
+  dopamin_level = decay_dopamin(dopamin_level, td);
+}
+
+double decay_dopamin(double level, double td) {
+  return level * exp( - td / dopamin_tau );
+}
+
+template<class Action>
+bool Global::reflect<Action>(Action &a) {
+return
+  // neuron related
+  reflect(voltage_tau, "voltage_tau")
+  && reflect(voltage_base, "voltage_base")
+  && reflect(threshold, "threshold")
+
+  && reflect(Wmin, "Wmin")
+  && reflect(Wmax, "Wmax")
+  && reflect(sumWeight, "sumWeight")
+  && reflect(absoluteRefractoryTime, "absoluteRefractoryTime")
+
+  // dopamin
+  && reflect(dopamin_level, "dopamin_level")
+  && reflect(dopamin_tau, "dopamin_tau")
+
+  // STDP
+  && reflect(stdp_tau_plus, "stdp_tau_plus")
+  && reflect(stdp_tau_minus, "stdp_tau_minus")
+  && reflect(stdp_lambda_plus, "stdp_lambda_plus")
+  && reflect(stdp_lambda_minus, "stdp_lambda_minus")
+  && reflect(stdp_et_tau, "stdp_et_tau")
+
+  // IP
+  && reflect(ip_sliding_avg_tau, "ip_sliding_avg_tau")
+  && reflect(ip_lambda_R, "ip_lambda_R")
+  && reflect(ip_lambda_C, "ip_lambda_C")
+  && reflect(ip_dst_mom1, "ip_dst_mom1")
+  && reflect(ip_dst_mom2, "ip_dst_mom2")
+
+  // external noise
+  && reflect(en_current, "en_current")
+  && reflect(en_freq, "en_freq")
+
+  // trainer
+  && reflect(trainer_eval_delay, "trainer_eval_delay")
+  && reflect(trainer_numSymbols, "trainer_numSymbols")
+  && reflect(trainer_refractoryTime, "trainer_refractoryTime")
+  && reflect(trainer_rewardAmount, "trainer_rewardAmount")
+  && reflect(trainer_rd_c1, "trainer_rd_c1")
+  && reflect(trainer_rd_c1, "trainer_rd_c2")
+  && reflect(trainer_rd_c1, "trainer_rd_c3")
+  && reflect(trainer_rd_c1, "trainer_rd_c4");
+}
+
+static id_type Global::numElements() {
+  return 1;
+}
+
+static Global * singleton(int num) {
+  if (num==0) {
+    return &global;
+  }else{
+    DIE("requested global object with id != 0");
+  }
+}
diff --git a/code/core/global.h b/code/core/global.h
new file mode 100644
index 0000000..8fec3d4
--- /dev/null
+++ b/code/core/global.h
@@ -0,0 +1,63 @@
+#ifndef GLOBAL_H
+#define GLOBAL_H
+
+class Global {
+public:
+  Global();
+
+  // methods
+  void evolve(double td);
+  double decay_dopamin(double level, double td); // return how much the level _would_ have been decayed
+
+  // variables (and quasi constants)
+  //   * neuron related
+  double voltage_tau;  // [s]
+  double voltage_base; // [V] rest potential 
+  double threshold;    // [V] above which a spike is emitted
+
+  double Wmax, Wmin;   // [V]
+  double sumWeight;    // [V]
+  double absoluteRefractoryTime; // [V]
+
+  //   * dopamin
+  double dopamin_level; // [?]
+  double dopamin_tau;   // [s]
+
+  //   * STDP
+  double stdp_tau_plus;     // [s]
+  double stdp_tau_minus;    // [s]
+  double stdp_lambda_plus;  // [1]
+  double stdp_lambda_minus; // [1]
+  double stdp_et_tau;       // [s] eligibility trace tau
+
+  //   * IP
+  double ip_sliding_avg_tau;
+  double ip_lambda_R;
+  double ip_lambda_C;
+  double ip_dst_mom1;
+  double ip_dst_mom2;
+
+  //   * external noise
+  double en_current; // [V]
+  double en_freq; // [Hz/Neuron]
+
+  //   * trainer
+  double trainer_eval_delay; // [s]
+  int trainer_numSymbols;            // [1]
+  double trainer_refractoryTime;     // [s]
+  double trainer_rewardAmount;       // [?]
+  double trainer_rd_c1; // readout delay with
+  double trainer_rd_c2; //   delay = c1 + t*c2 + r*c3 + t*r'*c4
+  double trainer_rd_c3; // where
+  double trainer_rd_c4; //   r and r' are two distinct random numbers (equ. dist in [0,1))
+
+  // reflection
+  template<class Action> bool reflect(Action &a);
+  typedef int id_type;
+  static id_type numElements();
+  static Neuron * singleton(int num); // return neuron # num
+} global;
+
+Global &g = global;
+
+#endif //  GLOBAL_H
diff --git a/code/core/interface.cpp b/code/core/interface.cpp
new file mode 100644
index 0000000..c2d3a3b
--- /dev/null
+++ b/code/core/interface.cpp
@@ -0,0 +1,262 @@
+#include <unistd.h>
+#include <asm-generic/errno.h>
+#include <boost/lexical_cast.hpp>
+
+#include "type2name.h"
+#include "log.h"
+#include "min.h"
+
+#include "interface.h"
+
+/* action definitions (the stuff implementing the different reflection primtives) */
+
+// (abstract) action (to implement filtering stuff)
+template<class interface>
+class Action {
+public:
+  Action<interface>(interface *ifc) : ifc(ifc) {}
+
+  bool isFiltered(char *desc) { // true -> don't use this element
+    return ! (ifc->ife.empty() || ifc->ifo.count(desc));
+  }
+
+  interface *ifc;
+};
+
+// count the length of a record assuming binary encoding
+template<class interface>
+class ActionCountRecordSize : public Action<interface> {
+public:
+  ActionCountRecordSize<interface>(interface ifc) : Action<interface>(ifc), result(0) {}
+
+  template<class T>
+  bool operator() (T &val, char *desc) {
+    if (this->isFiltered(desc)) return true;
+    result += sizeof(val);
+    return true;
+  }
+
+  int result;
+};
+
+// put the chosen titles in a string
+template<class interface>
+class ActionWriteTitle : public Action<interface> {
+public:
+  ActionWriteTitle<interface>(interface ifc) : Action<interface>(ifc) {}
+
+  template<class T>
+  bool operator() (T &val, char * desc) {
+    if (this->isFiltered(desc)) return true;
+    this->ifc->bufWrite(desc);
+    this->ifc->bufWrite(" ");
+    return true;
+  }
+};
+
+// put the choosen types in a string
+template<class interface>
+class ActionWriteTypes : public Action<interface> {
+public:
+  ActionWriteTypes<interface>(interface ifc) : Action<interface>(ifc) {}
+
+  template<class T>
+  bool operator() (T &val, char * desc) {
+    if (this->isFiltered(desc)) return true;
+    this->ifc->buf.append(type2name<T>());
+    this->ifc->buf.append(" ");
+    return true;
+  }
+};
+
+// put the serialized values to the string (space-delimited)
+// put the binary values to the string (no delimiter)
+template<class interface>
+class ActionWriteValues : public Action<interface> {
+public:
+  ActionWriteValues<interface>(interface ifc) : Action<interface>(ifc) {}
+
+  template<class T>
+  bool operator() (T &val, char *desc) {
+    if (this->isFiltered(desc)) return true;
+    return this->ifc->bufWrite(val);
+  }
+};
+
+// read serialized values
+template<class interface>
+class ActionRead : public Action<interface> {
+public:
+  ActionRead<interface>(interface ifc) : Action<interface>(ifc) {}
+
+  template<class T>
+  bool operator() (T &val, char * desc) {
+    if (this->isFiltered(desc)) return true;
+    return this->ifc->bufRead(val);
+  }
+};
+
+/* implementation of InputInterface */
+
+template<class T>
+InputInterface<T>::InputInterface(int fd, bool binary, IfObjectFilter *ifo, IfElementFilter *ife)
+  : binary(binary), fd(fd), ifo(ifo), ifo(ife), buf(), rpos(0), npos(0), wpos(0), eof(false) {
+  buf.reserve(buf_size);
+}
+
+template<class T>
+bool InputInterface<T>::garantueeBuf(size_t size) {
+  // no chance?
+  if (eof) return (wpos - rpos >= size);
+  if (size >= buf_size/2) return false;
+
+  // relocate buffer?
+  if (rpos > buf_size/2) {
+    buf = buf.substr(rpos);
+    buf.reserve(buf_size);
+    npos -= rpos;
+    wpos -= rpos;
+    rpos = 0;
+  }
+
+  // read more?
+  do {
+    if (wpos < buf_size) {
+      ssize_t res = read(fd, buf.data(), buf_size - wpos);
+      switch (res) {
+      case 0: // EOF
+	eof = true;
+	break;
+	
+      case EAGAIN:
+	// TO DECIDE: wait?
+	// read again (if neccessary)
+	break;
+	
+      case EIO: case EBADF: case EINVAL: case EINTR: case EFAULT:
+	return false;
+
+      default:
+	// read res bytes
+	wpos += res;
+      }
+    }
+  } while ((wpos - rpos < size) || eof); // loop in case of EAGAIN
+
+  return wpos - rpos >= size;
+}
+
+// make sure that a "\n"-terminated, non-empty string is in the buffer (starting at rpos)
+template<class T>
+bool InputInterface<T>::garantueeLine() {
+  if (npos <= rpos) npos = rpos + 1;
+  while (true) {
+    if ((npos > wpos) && !garantueeBuf(npos - rpos)) return false; // increase amount of data in buffer
+    if (buf[npos] == '\n') return true;
+    npos++;
+  }
+}
+
+template<class T>
+double InputInterface<T>::peekNextTime() {
+  double res;
+  if (!bufRead(res, false)) {
+    if (eof) {
+      return INFINITY;
+    }else{
+      DIE("peeking next time failed");
+    }
+  }
+  return res;
+}
+
+
+template<class Tif> template<class Tval>
+bool InputInterface<Tif>::bufRead(Tval &val, bool proceedRPos=true) {
+  if (binary) {
+    if (!garantueeBuf(sizeof(Tval))) return false;
+    val = *((Tval*) &(buf[npos]));
+    if (proceedRPos) rpos += sizeof(Tval);
+  }else{
+    if (!garantueeLine()) return false;
+
+    // find next delimiter
+    size_t p = buf.find_first_of(" \n", rpos);
+    if (!p || (p == buf.npos)) return false;
+    
+    // parse substring
+    try {
+      val = boost::lexical_cast<Tval>(buf.substr(rpos, p-1));
+    } catch(...) {
+      return false;
+    }
+
+    if (proceedRPos) rpos = p + 1;
+  }
+  return true;
+}
+
+template<class T>
+bool InputInterface<T>::readEntireFile() {
+  return readFileUntil(INFINITY);
+}
+
+template<typename T>
+bool InputInterface<T>::readFileUntil(double time) {
+  // iterate over events
+  while (peekNextTime() <= min(time, INFINITY)) {
+    // read (and forget) time to proceed rpos
+    double _foo; bufRead(_foo);
+
+    // get element id
+    typename T::id_type id;
+    if (!bufRead(id)) return false;
+    // TODO: check if the id refers to a valid object
+
+    // ignore it?
+    if (!ifo->empty() && !ifo->count(id)) return false;
+
+    // construct interface action and apply it
+    ActionRead<InputInterface<T> > ia(this);
+    if (!(T::singleton(id)->interface(ia))) return false;
+  }
+  return true;
+}
+
+template<class T>
+OutputInterface<T>::OutputInterface(int fd, bool binary, IfObjectFilter *ifo, IfElementFilter *ife)
+  : fd(fd), binary(binary), ifo(ifo), ife(ife), buf() {}
+
+template<class T>
+bool OutputInterface<T>::pushObject(T *o) {
+  ActionWriteValues<OutputInterface<T> > ia(this);
+  return o->interface(ia);
+}
+
+template<class T>
+bool OutputInterface<T>::pushClass() {
+  if (ifo->empty()) {
+    // for every object of the class
+    for (int i=0; i<T::numElements(); i++)
+      if (!pushObject(T::singleton(i))) return false;
+  }else{
+    // for all objects specified by ifo
+    for (IfObjectFilter::iterator i = ifo->begin(); i != ifo->end(); i++)
+      if (!pushObject(T::singleton(*i))) return false;
+  }
+}
+
+template<class T>
+bool OutputInterface<T>::isContained(typename T::id_type id) {
+  return ifo->empty() || ifo->count(id);
+}
+
+template<class Tif> template<class Tval>
+bool OutputInterface<Tif>::bufWrite(Tval &val) {
+  try {
+    buf.append(boost::lexical_cast<std::string>(val));
+  }catch (...) {
+    return false;
+  }
+  return true;
+}
diff --git a/code/core/interface.h b/code/core/interface.h
new file mode 100644
index 0000000..b53b741
--- /dev/null
+++ b/code/core/interface.h
@@ -0,0 +1,62 @@
+#ifndef INTERFACE_H
+#define INTERFACE_H
+
+#include <string>
+#include <set>
+#include <boost/lexical_cast.hpp>
+
+#include "math.h"
+
+// filters to select which objects/variables (not) to select
+typedef std::set<int> IfObjectFilter;
+typedef std::set<std::string> IfElementFilter;
+
+// interface-class .. holds long-time stuff for reflecting on a per class (template), per file (instance) and per direction (class) basis 
+template<class T>
+class InputInterface {
+public:
+  // A. the methods you want to call
+  // creation
+  InputInterface<T>(int fd, bool binary, IfObjectFilter ifo, IfElementFilter ife);
+  bool readEntireFile();
+  bool readFileUntil(double time);
+  double peekNextTime(); // read a line ahead (if possible) and check at what time it occurs
+
+  // B. internal state
+  bool binary;  
+  int fd;
+  std::string buf;
+  size_t rpos, npos, wpos;
+  static const size_t buf_size = 32768;
+  bool eof;
+  IfObjectFilter ifo;
+  IfElementFilter ife;
+
+  // C. internal functions
+  bool garantueeBuf(size_t size); // returns only after wpos - rpos >= size, but reads as much as possible without blocking
+  bool garantueeLine();
+
+  template<class Tval> bool bufRead(Tval &val, bool proccedRPos);
+};
+
+template<class T>
+class OutputInterface {
+public:
+  OutputInterface<T>(int fd, bool binary, IfObjectFilter *ifo, IfElementFilter *ife);
+
+  bool pushObject(T *o); // serialize one object
+  bool pushClass(); // serialize the entire class (selective by filter)
+  bool isContained(typename T::id_type id);
+
+  // internal state
+  int fd;
+  std::string buf;
+  bool binary;
+  IfObjectFilter *ifo;
+  IfElementFilter *ife;
+
+  // internal functions
+  template<class Tval> bool bufWrite(Tval &val);
+};
+
+#endif // INTERFACE_H
diff --git a/code/core/log.h b/code/core/log.h
new file mode 100644
index 0000000..57fb21b
--- /dev/null
+++ b/code/core/log.h
@@ -0,0 +1,7 @@
+#ifndef LOG_H
+#define LOG_H
+
+#define WARN(str) fprintf(stderr, "WARN at " __FILE__ ":%d: " str "\n", __LINE__, str);
+#define DIE(str) { fprintf(stderr, "FATAL at " __FILE__ ":%d: " str "\n", __LINE__, str); exit(-1); }
+
+#endif
diff --git a/code/core/max.h b/code/core/max.h
new file mode 100644
index 0000000..490b338
--- /dev/null
+++ b/code/core/max.h
@@ -0,0 +1,13 @@
+#ifndef MAX_H
+#define MAX_H
+
+template<class L, class R>
+inline L max(L l, R r) {
+  if (l > r) {
+    return l;
+  }else{
+    return r;
+  }
+}
+
+#endif // MAX_H
diff --git a/code/core/min.h b/code/core/min.h
new file mode 100644
index 0000000..78472ab
--- /dev/null
+++ b/code/core/min.h
@@ -0,0 +1,13 @@
+#ifndef MIN_H
+#define MIN_H
+
+template<class L, class R>
+inline L min(L l, R r) {
+  if (l < r) {
+    return l;
+  }else{
+    return r;
+  }
+}
+
+#endif // MIN_H
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]);
+}
diff --git a/code/core/neuron.h b/code/core/neuron.h
new file mode 100644
index 0000000..03c0dbb
--- /dev/null
+++ b/code/core/neuron.h
@@ -0,0 +1,64 @@
+/*************************************************************
+neuron.h
+  implementation the axon hillock
+
+  STORES
+  * current tuning of neuron body (parameters)
+  * list of synapses (-> connections to other neurons)
+**************************************************************/
+
+#ifndef NEURON_H
+#define NEURON_H
+
+#include <boost/intrusive/slist.hpp>
+
+#include "synapse.h"
+
+using namespace boost::intrusive;
+
+typedef slist<Synapse, member_hook<Synapse, slist_member_hook<>, &Synapse::hook_dst> > SynapseDstList;
+typedef slist<Synapse, member_hook<Synapse, slist_member_hook<>, &Synapse::hook_src> > SynapseSrcList;
+
+class Neuron {
+public:
+  Neuron();
+  void init();
+
+  // functions operating on neuronal state
+  double evolve(double time);
+  double processCurrent(double time, double current);
+  double generateSpike(double time, bool &doesOccur);
+  double predictSpike(double time);
+  void intrinsicPlasticity(double td);
+
+  // reflection
+  template<class Action> bool reflect(Action &a);
+  typedef int id_type;
+  static id_type numElements();
+  static Neuron * singleton(int num); // return neuron # num
+
+  // list of outgoing and incoming synapses
+  SynapseSrcList sin;
+  SynapseDstList sout;
+  //Synapse **sin, *sout;// list of _in_coming and _out_going _s_ynapses; HINT: sin is a list of pointers to synapses in the corresponding outgoing neuron
+
+  // basic neuron properties
+  double voltage; // [V]
+  double refractoryTime; // [s]
+
+  // IP related
+  bool ip; // if this neuron has IP at all
+  double ip_est_mom1, ip_est_mom2; // estimated moments (using exponential sliding average (exact integration))
+  double ip_dst_mom1, ip_dst_mom2; // targeted moments
+  double ip_R, ip_C; // internal correspondences to christinas model
+  double fac_voltage_tau, fac_current; // resulting coefficients for neuron behaviour
+
+  // local clocks
+  double lastEvent; // timepoint when the last event occured
+  double lastSpike;
+};
+
+const int numNeurons = 1000;
+Neuron n[numNeurons];
+
+#endif // NEURON_H
diff --git a/code/core/print_defaults.cpp b/code/core/print_defaults.cpp
new file mode 100644
index 0000000..7fbc7c2
--- /dev/null
+++ b/code/core/print_defaults.cpp
@@ -0,0 +1,10 @@
+#include <stdio.h>
+
+
+
+int main () {
+  TO BE IMPLEMENTED AGAIN
+  using the relflection interface
+
+  return 0;
+}
diff --git a/code/core/regex.cpp b/code/core/regex.cpp
new file mode 100644
index 0000000..ef98209
--- /dev/null
+++ b/code/core/regex.cpp
@@ -0,0 +1,134 @@
+// file access
+#include <fcntl.h>
+#include <errno.h>
+#include <list>
+#include <boost/xpressive/xpressive_static.hpp>
+#include <boost/xpressive/regex_actions.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/foreach.hpp>
+
+#include "interface.h"
+
+#include "log.h"
+#include "simulate.h"
+#include "tracepoints.h"
+#include "global.h"
+#include "neuron.h"
+#include "synapse.h"
+#include "spike.h"
+
+#include "regex.h"
+
+// pollute seperate namespace with regex expressions
+namespace tracepoints {
+
+  using namespace boost::xpressive;
+
+  /* parser result variables */
+  std::set<int> intSet;
+  std::set<std::string> stringSet;
+  double targetTimeOffset = 0.0;
+  std::string filename;
+  bool ioDirection; // tue -> write; false -> read
+  bool ioBinary;
+
+  /* parser callback implementations */
+  struct push_impl
+  {
+    typedef void result_type;     // Result type, needed for tr1::result_of
+    
+    template<typename Set, class Value>
+    void operator()(Set &s, Value const &val) const { s.insert(val); }
+  };
+  function<push_impl>::type const push = {{}};
+
+  /* I/O interface launcher */
+  template<class T>
+  struct start_interface_impl
+  {
+    typedef void result_type;     // Result type, needed for tr1::result_of
+    
+    void operator() () const {
+      // open src/dst file
+      int fd = open(filename.c_str, write ? (O_WRONLY | O_APPEND | O_CREAT) : (O_RDONLY | O_NONBLOCK));
+      if (fd == -1) DIE("could not open file");
+
+      // execute task
+      if (ioDirection) { // write
+	if (boost::is_same<T, SpikeMUX>::value) {
+	  // create a persistent interface to write out spikes as they occur
+	  s.spikeOIfList.push_back(new OutputInterface<SpikeMUX>(fd, ioBinary, stringSet, intSet));
+	}else{
+	  // create a temporary interface to write out any state
+	  OutputInterface<T> oif(fd, ioBinary, stringSet, intSet);
+	  oif.pushClass();
+	  close(fd);
+	}
+      }else{ // read
+	// create long living input interface (terminates itself when EOF is discovered)
+	FileInputEvent<T>::createInputStream(new InputInterface<T>(fd, ioBinary, stringSet, intSet));
+      }
+      
+      // clean up vars for next command
+      ioBinary = false;
+      stringSet.empty();
+      intSet.empty();
+    }
+  };
+
+  function<start_interface_impl<Global> >::type const start_interface_global = {{}};
+  function<start_interface_impl<Neuron> >::type const start_interface_neuron = {{}};
+  function<start_interface_impl<Synapse> >::type const start_interface_synapse = {{}};
+  function<start_interface_impl<Spike> >::type const start_interface_spike = {{}};
+
+  /* regex pattern definitions */
+
+  // character sets
+  cregex Delim = as_xpr(';');
+  cregex VarNameChar = alnum | as_xpr('_');
+
+  // basic elements
+  cregex Time = (+_d >> !('.' >> +_d))[ ref(targetOffseTime)=as<double>(_) ];
+
+  cregex Index = (+_d)[ push(as<int>(_)) ];
+
+  cregex Element = (+varnamechar)[ push(as<std::string>(_)) ];
+
+  cregex Filename = (+varnamechar)[ ref(filename)=<std::string>(_) ];
+
+  cregex IODirWrite = as_xpr('>') [ ref(ioDirection)=true ];
+  cregex IODirRead = as_xpr('>') [ ref(ioDirection)=false ];
+  cregex IODirection = IODirWrite | IODirRead;
+
+  cregex IOBinaryTrue = as_xpr('#') [ ref(ioBinary)=true ];
+  cregex IOBinary = !IOBinaryTrue;
+
+  cregex IOCmd = IOBinary >> IODirection >> Filename;
+
+  // lists
+  cregex ObjectListInner = *space >> (s1= index) >> *space;
+  cregex ObjectList = *space >> '(' >> object_list_inner >> *( delim >> object_list_inner ) >> ')';
+
+  cregex ElementListInner = *space >> (s1= element) >> *space;
+  cregex ElementList = *space >> '{' >> element_list_inner >> *( delim >> element_list_inner ) >> '}';
+
+  // trace commands
+  cregex TraceCmdTail = *space >> !ElementList >> *space >> !ObjectList >> *space >> IOCmd;
+
+  cregex TraceGlobal  = (icase( "global" )  >> TraceCmdTail [ start_interface_global() ];
+  cregex TraceNeuron  = (icase( "neuron" )  >> TraceCmdTail [ start_interface_neuron() ];
+  cregex TraceSynapse = (icase( "synapse" ) >> TraceCmdTail [ start_interface_synapse() ];
+  cregex TraceSpike   = (icase( "spike" )   >> TraceCmdTail [ start_interface_spike() ];
+
+  // whole line
+  // target format  what
+  //  cregex channel = *space >> +delim >> *space >> ( neuron | synapse | global | spikes );
+  //  cregex traceline = time  >> *channel >> *( delim | space );
+  cregex traceLine = !icase("proceed") >> +space >> time;
+  cregex line = *space >> *(inputElem | outputElem | traceElem) >> *( delim | space ));
+
+  bool parseRequest(char *str) {
+    return regex_match(str, line);
+  }
+  
+};
diff --git a/code/core/regex.h b/code/core/regex.h
new file mode 100644
index 0000000..0d1f3d3
--- /dev/null
+++ b/code/core/regex.h
@@ -0,0 +1,8 @@
+#ifndef REGEX_H
+#define REGEX_H
+
+namespace regex {
+  bool parseRequest(char *str);
+};
+
+#endif // REGEX_H
diff --git a/code/core/reward.cpp b/code/core/reward.cpp
new file mode 100644
index 0000000..0738c19
--- /dev/null
+++ b/code/core/reward.cpp
@@ -0,0 +1,105 @@
+#include <vector>
+
+#include "simulate.h"
+#include "bin.h"
+
+#include "reward.h"
+
+// only called once to init the da-reward system; later the more specific events are called
+void Reward::vexecute() {
+  // for ease of viewing ascendingly select 50 neurons per symbol
+  inputNeurons  = new IOPop(g.trainer_numSymbols);
+  outputNeurons = new IOPop(g.trainer_numSymbols);
+  for (int i = 0; i<g.trainer_numSymbols; i++)
+    for (int j=i*50; j<(i+1)*50; j++) {
+      (*inputNeurons)[i].insert(j);
+      (*outputNeurons)[i].insert(799 - j);
+    }
+
+
+  // present the first symbol after 0.5 second
+  s.addEvent(new Reward_Input(time + 10, inputNeurons, outputNeurons));
+}
+
+// choose and present a random input, then wait a random time
+void Reward_Input::vexecute() {
+  int symbol = rand() % g.trainer_numSymbols;
+  fprintf(stderr, "Pin: %d                              \n", symbol);
+
+  // excite symbol specific input neurons with 0.1 V
+  for (std::set<int>::iterator i = (*inputNeurons)[symbol].begin(); i != (*inputNeurons)[symbol].end(); i++)
+    s.addEvent(new ExternalSpike(time + drand48() * 0.001, *i, 0.15));
+
+  // wait a (not yet) random time until evaluation
+  s.addEvent(new Reward_EnableBinning(time + g.trainer_eval_delay, inputNeurons, outputNeurons, symbol));
+}
+
+// start the binning of relevant neurons
+void Reward_EnableBinning::vexecute() {
+  // add bins
+  vector<Bin*> *bins = new vector<Bin*>(g.trainer_numSymbols);
+  for (int i=0; i < (*outputNeurons).size(); i++) {
+    Bin *b = new Bin(&((*outputNeurons)[i]));
+    s.to.traceBinSets.insert(b);
+    (*bins)[i] = b;
+  }
+
+  // wait a small amount until reading the result of this binning
+  double delay = g.trainer_rd_c1
+               + g.trainer_rd_c2 * time
+               + g.trainer_rd_c3 * drand48()
+               + g.trainer_rd_c4 * time * drand48();
+  s.addEvent(new Reward_Readout(time + delay, symbol, inputNeurons, outputNeurons, bins));
+}
+
+// read the output frequencies and give reward
+void Reward_Readout::vexecute() {
+  if (estimatePerformance() > 1.0) {
+    deployReward(g.trainer_rewardAmount);
+  }else{
+    deployReward(-g.trainer_rewardAmount);
+  }
+
+  // delete bin trace commands
+  for (int i=0; i<bins->size(); i++) {
+    s.to.traceBinSets.erase((*bins)[i]);
+    delete (*bins)[i];
+  }
+  delete bins;
+
+  // wait a refractory time
+  s.addEvent(new Reward_Input(time + g.trainer_refractoryTime, inputNeurons, outputNeurons));
+}
+
+double Reward_Readout::estimatePerformance() {
+  int max_freq = 0, // max. freq of all populations _except_ target population
+    target_freq;
+
+  for (int i=0; i<bins->size(); i++) {
+    fprintf(stderr, "Pout: %d -> %d              \n", i, (*bins)[i]->count);
+    if (i == symbol) {
+      target_freq = (*bins)[i]->count;
+    }else{
+      max_freq = max(max_freq, (*bins)[i]->count);
+    }
+  }
+
+  double res;
+  if (max_freq == 0) {
+    if (target_freq == 0) {
+      res = 0.0;
+    }else{
+      res = ((double) target_freq) * INFINITY;
+    }
+  }else{
+    res = ((double) target_freq) / max_freq;
+  }
+
+  fprintf(stderr, "PERF: %f\n", res);
+  return res;
+}
+
+void Reward_Readout::deployReward(double reward) {
+  g.dopamin_level += reward;
+  s.da_history.push_front(pair<double,double>(time, g.dopamin_level));
+}
diff --git a/code/core/reward.h b/code/core/reward.h
new file mode 100644
index 0000000..1bc13f4
--- /dev/null
+++ b/code/core/reward.h
@@ -0,0 +1,53 @@
+#ifndef REWARD_H
+#define REWARD_H
+
+#include <vector>
+#include <set>
+
+#include <queue>
+
+#include "event.h"
+#include "bin.h"
+
+using namespace std;
+
+class Reward : public VirtualEvent {
+public:
+  typedef vector< std::set<int> > IOPop;
+
+  Reward(double time) : VirtualEvent(time) {}
+  Reward(double time, IOPop *inputNeurons, IOPop *outputNeurons) : VirtualEvent(time), inputNeurons(inputNeurons), outputNeurons(outputNeurons) {}
+  virtual void vexecute();
+
+  // state
+  IOPop *inputNeurons, *outputNeurons;
+};
+
+class Reward_Input : public Reward {
+public:
+  Reward_Input(double time, IOPop *inputNeurons, IOPop *outputNeurons) : Reward(time, inputNeurons, outputNeurons) {}
+  virtual void vexecute();  
+};
+
+class Reward_EnableBinning : public Reward {
+public:
+  Reward_EnableBinning(double time, IOPop *inputNeurons, IOPop *outputNeurons, int symbol) : Reward(time, inputNeurons,  outputNeurons), symbol(symbol) {}
+  virtual void vexecute();
+
+  int symbol;
+};
+
+class Reward_Readout : public Reward {
+public:
+  Reward_Readout(double time, int symbol, IOPop *inputNeurons, IOPop *outputNeurons, vector<Bin*> *bins) : Reward(time, inputNeurons, outputNeurons), symbol(symbol), bins(bins) {}
+  virtual void vexecute();
+
+  double estimatePerformance();
+  void deployReward(double reward);
+
+  int symbol;
+  vector<Bin*> * bins;
+};
+
+
+#endif // REWARD_H
diff --git a/code/core/simulate.cpp b/code/core/simulate.cpp
new file mode 100644
index 0000000..69252bb
--- /dev/null
+++ b/code/core/simulate.cpp
@@ -0,0 +1,134 @@
+#include <utility>
+#include <map>
+#include <signal.h>
+
+#include "model_switch.h"
+#include "synapse.h"
+#include "neuron.h"
+#include "topology.h"
+#include "reward.h"
+#include "tracepoints.h"
+#include "fileutils.h"
+#include "log.h"
+
+#include "simulate.h"
+
+// include cpp files to allow all otpimizations to take place across all files
+#include "synapse.cpp"
+#include "neuron.cpp"
+#include "topology.cpp"
+#include "fileutils.cpp"
+#include "event.cpp"
+#include "reward.cpp"
+#include "bin.cpp"
+#include "tracepoints.cpp" // should be excluded to speed up compilation
+
+
+Simulation::Simulation() :
+  // debug related stuff
+#ifdef DEBUG_STATUSLINE
+  numSpikes(0),
+  charCount(0),
+#endif
+
+  // init globals
+  currentTime(0.0),
+{
+  // set the initial dopamin level
+  da_history.push_front(pair<double, double>(0.0, g.dopamin_level));
+}
+
+bool Simulation::Step() {
+  // check if there are pending events
+  if (pendingSpikes.empty())
+    return false;
+
+  // retrieve and check next event
+  Event *e = pendingSpikes.top();
+  pendingSpikes.pop();
+
+  // proceed to the new time
+  if (currentTime > e->time)
+    DIE("tried to execute event of the past");
+  currentTime = e->time;
+
+#ifdef DEBUG_STATUSLINE
+  if (numSpikes % 16384 == 0) {
+    // remove old line
+    while (charCount-- > 0)
+      fprintf(stderr, " ");
+    fprintf(stderr, "\r");
+
+    // print new line
+    charCount += fprintf(stderr, "%f s, %d, %lld events (%d\tpending:", e->type, currentTime, numSpikes, pendingSpikes.size());
+    for (map<int, int>::iterator i = eventCount.begin(); i != eventCount.end(); i++) {
+      charCount += fprintf(stderr, "\t%d: %d", i->first, i->second);
+      charCount += 7; // tab
+    }
+    charCount += fprintf(stderr, ")\r");
+    charCount += 7; // tab
+  }
+  fflush(stderr);
+
+  numSpikes++;
+  eventCount[e->type]--;
+#endif
+
+  // execute event
+  e->execute();
+  e->free();
+
+  return true;
+}
+
+bool Simulation::proceedTime(double targetTime) {
+  // insert an intrinsic event for each neuron to have it computed up to exactly targetTime
+  for (int i=0; i < numNeurons; i++) {
+    addEvent(new IntrinsicNeuron(targetTime, i));
+  }
+
+  // proceed until no event is left before the time-point to be reached
+  // then also currentTime == targetTime holds
+  while (!pendingSpikes.empty() && (pendingSpikes.top()-> time <= targetTime)) {
+    Step();
+  }
+
+  // check if there are events left ... with the new event system this should always be the case
+  // and does not anymore indicate that the network is alive
+  return !pendingSpikes.empty();
+}
+
+void Simulation::addEvent(Event *e) {
+  pendingSpikes.push(e);
+#ifdef DEBUG_STATUSLINE
+  eventCount[e->type]++;
+#endif
+}
+
+void initStaticVars() {
+  Synapse::numSynapses = 0;
+}
+
+int main(int argc, char **argv) {
+  // ignore broken pipe signals (they should be handled by stdio file handlers and our usual error functions)
+  signal(SIGPIPE, SIG_IGN);
+
+  initStaticVars();
+
+  // load the network
+  // tp_load(s.fd_isynapse);
+  // neuron_load(s.fd_ineuron);
+
+  // init services implemented via events by adding one of them to the event population
+  s.addEvent(new ExternalNoise(0.0));
+  s.addEvent(new Reward(0.0));
+
+  // check for autoexciting events (and init neurons btw)
+  s.proceedTime(0.0);
+
+  // pass control to the tracepoint parser
+  // this allows to run the simulation interactively
+  executeTracepoints(s.fd_trace, s.fd_ispike, s.fd_iglobal);
+
+  return 0;
+}
diff --git a/code/core/simulate.h b/code/core/simulate.h
new file mode 100644
index 0000000..3d84b6e
--- /dev/null
+++ b/code/core/simulate.h
@@ -0,0 +1,70 @@
+#ifndef SIMULATE_H
+#define SIMULATE_H
+
+#include <queue>
+#include <list>
+
+#include "neuron.h"
+#include "synapse.h"
+#include "tracepoints.h"
+#include "event.h"
+#include "interface.h"
+#include "spike.h"
+
+using namespace std;
+ 
+class Simulation {
+public:
+  Simulation();
+
+  // ----- FUNCTIONS -----
+
+  // controlling functions
+  bool Step(); // do a single step (fire one neuron); return if the net is dead
+  bool proceedTime(double td);
+  bool proceedSpikes(long count);
+
+  void addEvent(Event*);
+
+  // ----- SIMULATION STATE -----
+
+  // simulation global variables
+  double currentTime;
+
+  // list of spikes and events
+  priority_queue<Event, vector<Event*>, PEventGreater> pendingSpikes;
+
+  // list of discontinous dopamin changes (time, amount), sorted by time (newest is front)
+  std::list<pair<double, double> > da_history;
+
+  // ----- HELPER VARS -----
+
+  // reflection & traces
+  std::list<OutputInterface<SpikeMUX> *> spikeOIfList;
+
+  std::set<Bin*> binSets;
+
+  // file descriptors
+  FILE *fd_ineuron,
+    *fd_isynapse,
+    *fd_ispike,
+    *fd_iglobal,
+    *fd_oneuron,
+    *fd_osynapse,
+    *fd_ospike,
+    *fd_oglobal,
+    *fd_trace; // command file
+
+  // debug vars
+#ifdef DEBUG_STATUSLINE
+  long long numSpikes; // number of spikes (counted after axon hillock) processed so far
+  map<int, int> eventCount;
+  int charCount;
+#endif
+};
+
+// init neural network
+Simulation s;
+
+
+#endif // SIMULATE_H
diff --git a/code/core/spike.cpp b/code/core/spike.cpp
new file mode 100644
index 0000000..7b12c18
--- /dev/null
+++ b/code/core/spike.cpp
@@ -0,0 +1,28 @@
+#include "spike.h"
+
+static SpikeMUX * SpikeMUX::singleton(id_type num) {
+  staticSpikeMUX.dst = num;
+  return &staticSpikeMUX;  
+}
+
+// special case for reading a spike from a file
+template<>
+bool SpikeMUX::reflect<ActionRead>(ActionRead &a) {
+  bool res = reflect(dst, "dst")
+    && reflect(time, "time")
+    && reflect(current, "current")
+    && (time > s.currentTime);
+
+  if (res)
+    s.addEvent(new ExternalSpike(time, dst, current));
+
+  return res;
+}
+
+// general reflection case
+template<class Action>
+bool SpikeMUX::reflect<Action>(Action &a) {
+  return reflect(dst, "dst")
+    && reflect(time, "time")
+    && reflect(current, "current");
+}
diff --git a/code/core/spike.h b/code/core/spike.h
new file mode 100644
index 0000000..2b2e4da
--- /dev/null
+++ b/code/core/spike.h
@@ -0,0 +1,28 @@
+#ifndef SPIKE_H
+#define SPIKE_H
+
+#include <math.h>
+
+#include "neuron.h"
+#include "log.h"
+
+class SpikeMUX {
+public:
+  typedef Neuron::id_type id_type; // WARN: this id refers to the spiking neuron, not the spike itself wich can not be addressed directly
+  SpikeMUX(id_type dst) : dst(dst), time(-INFINITY), current(0.0) {}
+  SpikeMUX(id_type dst, double time) : dst(dst), time(time), current(0.0) {}
+
+  // reflection
+  template<class Action> bool reflect(Action &a);
+  static id_type numElements() { DIE("Spike::numElements() called"); }
+  static SpikeMUX * singleton(id_type num); // do some reflection magic to multiplex between IntrinsicNeuron and ExternalNoise
+
+  // properties
+  id_type dst;
+  double time;
+  double current;
+};
+
+SpikeMUX staticSpikeMUX(-1); // single copy to use for reading spikes
+
+#endif // SPIKE_H
diff --git a/code/core/synapse.cpp b/code/core/synapse.cpp
new file mode 100644
index 0000000..cc7399a
--- /dev/null
+++ b/code/core/synapse.cpp
@@ -0,0 +1,118 @@
+#include <math.h>
+
+#include "synapse.h"
+
+Synapse::Synapse() :
+  src(-1), dst(-1), // seeing this value indicates an unused synapse
+  weight(0.0),
+  eligibility(0.0),
+  delay(0.0),
+  constant(false),
+  firstSpike(- INFINITY),
+  lastSpike(- INFINITY),
+  evolvedUntil(0.0)
+{}
+
+// given an input time and current calculate their output values
+void Synapse::computePostsynapticPulse(double &time, double &current) {
+  time += delay;
+  current = weight;
+}
+
+void Synapse::evolve(double time) {
+  if (time < evolvedUntil)
+    return;
+
+  // get the dopamin history starting from a dopamin event earlier than the synpase's last event
+  std::list<pair<double, double> >::iterator i = s.da_history.begin();
+  while (i->first > evolvedUntil)
+    i++;
+
+  // evolve synapse state for each interval given by the dopamin history
+  while (evolvedUntil < time) {
+    double da_level = i->second;
+    if (i->first < evolvedUntil) // the first dopamin event partially lies outside our time window
+      da_level = g.decay_dopamin(da_level, evolvedUntil - i->first);
+
+    double td = - evolvedUntil;
+    if (i == s.da_history.begin()) {
+      // this is the last (=newest) element of the history
+      td += time;
+      evolvedUntil = time;
+    }else{
+      i--;
+      if (i->first > time) {
+	// dopamin history goes furher (in time) than we have to simulate
+	td += time;
+	evolvedUntil = time;
+      }else{
+	td += i->first;
+	evolvedUntil += td;
+      }
+    }
+
+    /* below is the exact solution of the equation system
+         dc/dt = -Cc          c - eligibility trace, C - cf. tau
+         dd/dt = -Dd          d - dopamin level, D - cf. tau
+         dw/dt = cd           w - weight
+    */
+    
+    if (!constant) {
+      double tau = g.dopamin_tau * g.stdp_et_tau / (g.dopamin_tau + g.stdp_et_tau);
+      double dw = da_level * eligibility * tau * (1 - exp( -td / tau ));
+      weight += dw;
+      weight = fmax( msg.Wmin, fmin( msg.Wmax, weight ));
+    }
+    eligibility *= exp( -td / g.stdp_et_tau );
+  }
+}
+
+// modify eligibility trace according to a single STDP event given by
+// a pre-post time difference (= t_post - t_pre)
+void Synapse::stdp(double td) {
+  if (td >= 0.0) {// WARN: this doesn't work with negative zero
+    eligibility = g.stdp_lambda_plus * exp( - td / g.stdp_tau_plus );
+  }else{
+    eligibility = - g.stdp_lambda_minus * exp( td / g.stdp_tau_minus );
+  }
+}
+
+static void Synapse::updateTopology() {
+  // clear the sin/sout lists of all neurons
+  for (int i=0; i<numNeurons; i++) {
+    n[i].sout.clear();
+    n[i].sin.clear();
+  }
+
+  // iterate over all synapses wich have a valid source and destination (neuron) and add them to sin/sout lists
+  int i;
+  for (i=0; (i<maxSynapses) && (syn[i].src != -1) && (syn[i].dst != -1); i++) {
+    n[syn[i].src].sout.push_front(syn[i]);
+    n[syn[i].dst].sin.push_front(syn[i]);
+  }
+
+  // store number of (used) synapses
+  numSynapses = i;
+}
+
+template<class Action>
+bool Synapse::reflect<Action>(Action &a) {
+  return
+    a(src, "src")
+    && a(dst, "dst")
+    && a(weight, "weight")
+    && a(eligibility, "eligibility")
+    && a(delay, "delay")
+    && a(constant, "constant")
+    && a(firstSpike, "firstSpike")
+    && a(lastSpike, "lastSpike")
+    && a(evolvedUntil, "evolvedUntil");
+}
+
+static id_type Synapse::numElements() {
+  return numSynapsesDec;
+}
+
+static Synapse * Synapse::singleton(id_type id) {
+  return &(synapses[id]);
+}
diff --git a/code/core/synapse.h b/code/core/synapse.h
new file mode 100644
index 0000000..3eef12d
--- /dev/null
+++ b/code/core/synapse.h
@@ -0,0 +1,48 @@
+#ifndef SYNAPSE_H
+#define SYNAPSE_H
+
+#include <boost/intrusive/slist.hpp>
+
+using namespace boost::intrusive;
+
+class Neuron;
+
+class Synapse {
+public:
+  Synapse();
+
+  // functions
+  void computePostsynapticPulse(double &time, double &current); // given an input time and current calculate their output value
+  void evolve(double time);
+  void stdp(double dt);
+  static void updateTopology(); // recreate all sin/sout lists in the neurons
+
+  // reflection
+  template<class Action> bool reflect(Action &a);
+  typedef int id_type;
+  static id_type numElements();
+  static Synapse * singleton(id_type num); // return neuron # num
+  static id_type numSynapses;
+
+  // model specifics
+  int src, dst;       // pre- and postsynaptic neuron
+  double weight;      // [V]
+  double eligibility; // [?]
+  double delay;       // total delay from presynaptic axon hillock to postsynaptic axon hillock
+  bool constant;      // if the weight can change (HINT: weight < 0 implies constness, too)
+
+  // sim helper vars
+  double firstSpike,   // time the first/last spike arrived at it's dst neuron in a given interval
+         lastSpike;    // HINT: firstSpike is set to -INFINITY every time after it is processed
+                       // by the STDP rule
+  double evolvedUntil; // up to wich time has the model been evolved
+
+  // hock to use intrusive lists (defined in neuron.h)
+  slist_member_hook<> hook_src, hook_dst;  
+};
+
+const int maxSynapses = 200000;
+Synapse syn[maxSynapses];
+
+
+#endif // SYNAPSE_H
diff --git a/code/core/tracepoints.cpp b/code/core/tracepoints.cpp
new file mode 100644
index 0000000..7243f1c
--- /dev/null
+++ b/code/core/tracepoints.cpp
@@ -0,0 +1,57 @@
+#include "event.h"
+#include "regex.h"
+
+#include "simulate.h"
+
+/* input event streams */
+
+template<class T>
+static void FileInputEvent<T>::CreateInputStream(InputInterface<T> *iface) {
+  double time = iface->peekNextTime();
+  if (time == INFINITY) {
+    delete iface;
+    return;
+  }
+  if (time < s.currentTime) {
+    DIE("tried to include a file with events of the past");
+  }
+  s.addEvent(new FileInputEvent(iface, time));
+}
+
+template<class T>
+void FileInputEvent<T>::vexecute() {
+  iface->readFileUntil(time);
+  CreateInputStream(iface);
+}
+
+/* command loop */
+
+bool executeTracepoints() {
+  // the main loop of the program
+  char *str = NULL;
+  size_t _foo;
+
+  // loop over trace commands
+  while (getline(&str, &_foo, stdin) > 0) {
+    // parse request
+    if (!regex::parseRequest(str))
+      DIE("Invalid tracepoint format");
+
+    // proceed time if a run command was given
+    if (regex::targetTimeOffset != 0.0) {
+      if (!s.proceedTime(s.currentTime + tracepoints::targetTimeOffset))
+	fprintf(stderr, "Warning: network is dead\n");
+      // reset target time
+      targetTimeOffset = 0.0;
+      // reset spike output list
+      BOOST_FOREARCH(SpikeMUX *i, s.spikeOIfList) { delete i; }
+      s.spikeOIfList.empty();
+    }
+
+
+    free(str);
+    str = NULL;
+  }
+  
+  return true;
+}
diff --git a/code/core/tracepoints.h b/code/core/tracepoints.h
new file mode 100644
index 0000000..d0ea9fe
--- /dev/null
+++ b/code/core/tracepoints.h
@@ -0,0 +1,25 @@
+#ifndef TRACEPOINTS_H
+#define TRACEPOINTS_H
+
+#include <list>
+#include <set>
+
+#include "interface.h"
+#include "bin.h"
+#include "event.h"
+
+// read tracepoints from fd and execute them
+bool executeTracepoints(FILE *fd_trace, FILE *fd_spike, FILE *fd_global);
+
+template<class T>
+class FileInputEvent : public VirtualEvent {
+public:
+  FileInputEvent<T>(InputInterface<T> *iface, double time) : VirtualEvent(time), iface(iface) {}
+  static void createInputStream(InputInterface<T> *iface);
+  void vexecute();
+
+  InputInterface<T> *iface;
+};
+
+
+#endif // TRACEPOINTS_H
diff --git a/code/core/type2name.h b/code/core/type2name.h
new file mode 100644
index 0000000..1917a6a
--- /dev/null
+++ b/code/core/type2name.h
@@ -0,0 +1,20 @@
+#ifndef TYPE2NAME_H
+#define TYPE2NAME_H
+
+// if the type is not known
+template<class T>
+char * type2name() { return "type_not_stringified"; }
+
+template<>
+char * type2name<bool>() { return "bool"; }
+
+template<>
+char * type2name<char>() { return "char"; }
+
+template<>
+char * type2name<int>() { return "int"; }
+
+template<>
+char * type2name<double>() { return "double"; }
+
+#endif // TYPE2NAME_H