diff options
Diffstat (limited to 'code')
58 files changed, 3453 insertions, 0 deletions
diff --git a/code/Makefile b/code/Makefile new file mode 100644 index 0000000..804a0ec --- /dev/null +++ b/code/Makefile @@ -0,0 +1,9 @@ +.PHONY: all clean + +SUBDIRS=core trainer matlab + +all: + for I in $(SUBDIRS); do cd $$I; make; cd ..; done + +clean: + for I in $(SUBDIRS); do cd $$I; make clean; cd ..; done diff --git a/code/anatool/plot_spikes b/code/anatool/plot_spikes new file mode 100755 index 0000000..58fe722 --- /dev/null +++ b/code/anatool/plot_spikes @@ -0,0 +1,12 @@ +#!/bin/bash + +( + if [ $# -eq 1 ]; then + echo "set xrange [0:$1]" + fi + if [ $# -eq 2 ]; then + echo "set xrange [$1:$2]" + fi + echo "plot 'spikes.out' using 1:2 with dots" +) | gnuplot -persist + 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 ¤t) { + 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 ¤t); // 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 diff --git a/code/glue/Makefile b/code/glue/Makefile new file mode 100644 index 0000000..35bd16b --- /dev/null +++ b/code/glue/Makefile @@ -0,0 +1,8 @@ +.PHONY: all clean +all: repeat-trace-cmd + +clean: + rm *~ repeat-trace-cmd + +repeat-trace-cmd: repeat-trace-cmd.c + gcc repeat-trace-cmd.c -o repeat-trace-cmd -O23
\ No newline at end of file diff --git a/code/glue/da-controlled-sim-wrapper b/code/glue/da-controlled-sim-wrapper new file mode 100755 index 0000000..918de4a --- /dev/null +++ b/code/glue/da-controlled-sim-wrapper @@ -0,0 +1,62 @@ +#!/bin/sh + +# check param count +if [ ! $# -eq 7 ]; then + echo 'wrong parameter count (see the source for parameter order)' >&2 + # 1. model name (e.g. the "if" from "sim-if") + # 2. controller name (relative to trainer-dir) + # 3. input_neuron_file + # 4. input_synapse_file + # 5. output_neuron_file + # 6. output_synapse_file + # 7. performance output + exit 1 +fi + +# determine the path of the simulaton program +SIM=`dirname $0`/../core/sim-$1 +if [ ! -x $SIM ]; then + echo "executable ($SIM) does not exist" >&2 + exit 1 +fi + +# determine the path of the controller programm +CTL=`dirname $0`/../trainer/$2-$1 +if [ ! -x $CTL ]; then + echo "executable ($CTL) does not exist" >&2 + exit 1 +fi + + +# create tmp dir +FIFODIR=`mktemp -td fasimu.XXXXXXXXXX` + +# create fifos +mkfifo $FIFODIR/spike_in +mkfifo $FIFODIR/spike_out +mkfifo $FIFODIR/trace_in +mkfifo $FIFODIR/global_in +mkfifo $FIFODIR/global_out + +# TODO: check if an additional i/o file is an executable + +# launch controller and simulator +#echo $CTL - $FIFODIR/trace_in $FIFODIR/global_in $FIFODIR/global_out $FIFODIR/spike_in $FIFODIR/spike_out "2> trainer.err &" +$CTL $7 $FIFODIR/trace_in $FIFODIR/global_in $FIFODIR/global_out $FIFODIR/spike_in $FIFODIR/spike_out 2> trainer.err & + +#echo $SIM "2> sim.err" $3 $4 $FIFODIR/spike_in $FIFODIR/global_in $5 $6 $FIFODIR/spike_out $FIFODIR/global_out $FIFODIR/trace_in +$SIM 2> sim.err $3 $4 $FIFODIR/spike_in $FIFODIR/global_in $5 $6 $FIFODIR/spike_out $FIFODIR/global_out $FIFODIR/trace_in + +# hint: simulator params are + # input_neuron_file + # input_synapse_file + # input_spike_file + # input_global_file + # output_neuron_file + # output_synapse_file + # output_spike_file + # output_global_file + # trace_commando_file + +# delete tmp dir +rm -R $FIFODIR diff --git a/code/glue/distill-performance b/code/glue/distill-performance new file mode 100755 index 0000000..af90f2f --- /dev/null +++ b/code/glue/distill-performance @@ -0,0 +1,5 @@ +#!/bin/bash + +if [ -f performance.out.raw ]; then + cat performance.out.raw | tr "\r" "\n" | grep ^PERF | cut -d" " -f2 > performance.out +fi diff --git a/code/glue/exec-matlab b/code/glue/exec-matlab new file mode 100755 index 0000000..ccb0d23 --- /dev/null +++ b/code/glue/exec-matlab @@ -0,0 +1,8 @@ +#!/bin/sh + +if [ -f $0.m ]; then + env - bash -c "(echo \"$1\"; cat $0.m) | matlab -nodesktop" +# env - bash -c "cat $0.m | octave" +else + echo file $0.m not found +fi diff --git a/code/glue/extract-matlab-matrix b/code/glue/extract-matlab-matrix new file mode 100755 index 0000000..b5021a4 --- /dev/null +++ b/code/glue/extract-matlab-matrix @@ -0,0 +1,43 @@ +#!/bin/bash + +# TODO: use tail instead of cat where only the top of the file is needed + +# extract the multiplier of the matrix +MUL=`cat $1 \ +| egrep -o 'e\+[0-9]* \*' \ +| tr -d 'e+ *'` + +if [ -z "$MUL" ]; then + MUL=1 +fi + +# get the number of cols +COLS=`cat $1 \ +| egrep '([:space:]*([01](\.[0-9]*){0,1})[:space:]*)+$' \ +| tail -n1 \ +| wc -w` + +# read the matrix, multiply to correct value and put it out +cat $1 \ +| egrep '([:space:]*([01](\.[0-9]*){0,1})[:space:]*)+$' \ +| tr " " "\n" \ +| egrep -v '^$' \ +| while read; do + echo $MUL '*' $REPLY +done \ +| bc \ +| sed 's/\([0-9]\+\)\.0*/\1/' \ +| ( I=1 + while read; do + if [ $I -eq $COLS ]; then + echo "$REPLY" + I=1 + else + echo -n "$REPLY," + I=$(( $I + 1 )) + fi + done ) + + +# old debug stuff +#echo $MUL $COLS
\ No newline at end of file diff --git a/code/glue/plot_sliding_perf b/code/glue/plot_sliding_perf new file mode 100755 index 0000000..c71d42b --- /dev/null +++ b/code/glue/plot_sliding_perf @@ -0,0 +1,24 @@ +#!/bin/bash + +if [ -f "performance.out" ]; then + echo "x=load('performance.out'); + x=(x > 1); + y=length(x); + p=zeros(y-100,1); + for i=1:(y-100) + p(i)=mean(x(i:(i+100),1)); + end + p" \ + | octave -q \ + | tr -dc "0123456789.\n" \ + |grep -v "^$" \ + > performance.out.sliding-avg + + PWD=`pwd` + + echo "set title 'performance (sliding avg, window size 100) $PWD' + set terminal postscript + set output 'performance.out.ps' + plot 'performance.out.sliding-avg' using 1 with lines, 0.5" \ + | gnuplot +fi diff --git a/code/glue/plot_spike_time_hist b/code/glue/plot_spike_time_hist new file mode 100755 index 0000000..ad45305 --- /dev/null +++ b/code/glue/plot_spike_time_hist @@ -0,0 +1,19 @@ +#!/bin/bash + +if [ -f "spikes.out" ]; then + cat spikes.out | cut -d, -f1 > spikes.out.timing + + echo "x = load('spikes.out.timing'); + y = hist(x, ceil(max(x)/10)); + y'" \ + | octave -q \ + | tr -dc "0123456789.\n" \ + | grep -v "^$" \ + > spikes.out.binned-timing + + echo "set title 'population frequency $PWD' + set terminal postscript + set output 'spikes.out.binned-timing.ps' + plot 'spikes.out.binned-timing' using 1 with lines" \ + | gnuplot +fi
\ No newline at end of file diff --git a/code/glue/print-params b/code/glue/print-params new file mode 100755 index 0000000..3df1e19 --- /dev/null +++ b/code/glue/print-params @@ -0,0 +1,3 @@ +#!/bin/bash + +echo $@ diff --git a/code/glue/repeat-trace-cmd b/code/glue/repeat-trace-cmd Binary files differnew file mode 100755 index 0000000..b889e0b --- /dev/null +++ b/code/glue/repeat-trace-cmd diff --git a/code/glue/repeat-trace-cmd.c b/code/glue/repeat-trace-cmd.c new file mode 100644 index 0000000..ac35733 --- /dev/null +++ b/code/glue/repeat-trace-cmd.c @@ -0,0 +1,27 @@ +#include <stdio.h> + +int main(int argc, char **argv) { + double t, dt; + long n; + + if (argc != 4) { + fprintf(stderr, "ERROR: wrong argument count\nUse %s total_time time_per_trace \"trace command(s) \"\n", argv[0]); + return -1; + } + + if ((sscanf(argv[1], "%lf", &t) != 1) || + (sscanf(argv[2], "%lf", &dt) != 1)) { + fprintf(stderr, "failed to read arg 1/2\n"); + return -1; + } + printf("%f, %f\n", t, dt); + // print the full command once + + // now print enough newline (= command repetitions) + // TODO: be faster than lame-duck-speed + n = (long) (t / dt); // on step already passed because of above printf statement + while (n>1) { + printf("\n"); + n--; + } +} diff --git a/code/glue/sim-wrapper b/code/glue/sim-wrapper new file mode 100755 index 0000000..9825409 --- /dev/null +++ b/code/glue/sim-wrapper @@ -0,0 +1,24 @@ +#!/bin/sh + +if [ ! $# -eq 10 ]; then + echo 'wrong parameter count (see ./sim-current for parameter order and add the model (current/if/...) as the first param)' >&2 + exit 1 +fi + +# determine the path of the simulaton program +SIM=`dirname $0`/../core/sim-$1 +if [ ! -x $SIM ]; then + echo "executable ($SIM) does not exist" >&2 + exit 1 +fi + +# check if one of the input files is executable +if [ -x $2 -o -x $3 -o -x $4 -o -x $5 -o -x $5 -o -x $6 -o -x $6 -o -x $7 -o -x $8 -o -x $9 ]; then + # yes -> interactive simulation + # create the FIFOs to communicate + echo Interactive spike program is not NOT IMPLEMENTED + exit 1 +else + # no -> static simulation + $SIM $2 $3 $4 $5 $6 $7 $8 $9 $10 +fi diff --git a/code/matlab/Makefile b/code/matlab/Makefile new file mode 100644 index 0000000..3546e44 --- /dev/null +++ b/code/matlab/Makefile @@ -0,0 +1,10 @@ +.PHONY: all symlinks clean + +all: symlinks + +symlinks: + ls | grep '.m$$' | sed 's/\.m$$//' | xargs -n1 ln -f -s ../glue/exec-matlab + +clean: + ls | grep '.m$$' | sed 's/\.m$$//' | xargs rm || true + rm *~ || true
\ No newline at end of file diff --git a/code/matlab/analye-perfomance.m b/code/matlab/analye-perfomance.m new file mode 100644 index 0000000..75b399b --- /dev/null +++ b/code/matlab/analye-perfomance.m @@ -0,0 +1,14 @@ +%% load data +perf = load('performance.out'); +perf = perf > 1.0; +s = length(perf); + +%% compute sliding avergae +ws = 250; % window size +perf_avg = zeros(s-ws,1); +for i = 1:(s-ws) + perf_avg(i,1) = mean(perf(i:(i+ws),1)); +end + +%% plot +plot(perf_avg);
\ No newline at end of file diff --git a/code/matlab/analye-stdp-freq-dep.m b/code/matlab/analye-stdp-freq-dep.m new file mode 100644 index 0000000..85c4030 --- /dev/null +++ b/code/matlab/analye-stdp-freq-dep.m @@ -0,0 +1,12 @@ + +raw = load('synapse.destilled'); + +res = []; +for i=1:(floor(length(raw) / 200)), + row = raw(((i*200)-199):((i*200)-99),2)'; + res(i,:) = [mean(row), var(row)]; +end + +for i=1:length(res), + fprintf(2,'%d, %f, %f\n', int32(i), res(i,1), res(i,2)) +end
\ No newline at end of file diff --git a/code/matlab/analyze_weight_development.m b/code/matlab/analyze_weight_development.m new file mode 100644 index 0000000..da34cd4 --- /dev/null +++ b/code/matlab/analyze_weight_development.m @@ -0,0 +1,30 @@ +%% -- load the synapse file + +syn_raw = load('synapse.out'); + +%% -- generate mean weights (and some config stuff) + +num_syn = 96227; +num_steps = length(syn_raw) / num_syn; +types = [ 1*ones(100,1); 2*ones(600,1); 3*ones(100,1); 4*ones(200,1) ]; + +syn_mean = zeros(num_steps, 16); +syn_count = zeros(num_steps, 16); +for i = 1:num_steps + ba = (i-1)*num_syn+1; + for j = ba:(ba+num_syn-1) + l = syn_raw(j,:); + src = l(1,2) + 1; + dst = l(1,3) + 1; + w = l(1,5); + + syn_mean(i, 4*types(src)+types(dst)-4) = syn_mean(i, 4*types(src)+types(dst)-4) + w; + syn_count(i, 4*types(src)+types(dst)-4) = syn_count(i, 4*types(src)+types(dst)-4) + 1; + end +end +syn_mean = syn_mean ./ syn_count; + +%% plot it + +plot(syn_mean); +legend('II', 'IB', 'IO', 'IX', 'BI', 'BB', 'BO', 'BX', 'OI', 'OB', 'OO', 'OX', 'XI', 'XB', 'XO', 'XX');
\ No newline at end of file diff --git a/code/matlab/plot_stdp_param_scout.m b/code/matlab/plot_stdp_param_scout.m new file mode 100644 index 0000000..2342fd7 --- /dev/null +++ b/code/matlab/plot_stdp_param_scout.m @@ -0,0 +1,51 @@ +%% load the raw data from synapses output file +raw = load('synapse.out'); + + +%% get mean and variance from the 1000 neurons +ns = 999; +l = floor(length(raw)/ns); +raw2 = zeros(l,1); +raw2_var = zeros(l,1); +for i=1:l, + % hint: adapted to read only the first 100 out of 1000 neurons + raw2(i,1) = mean(raw((i*ns-ns+1):(i*ns-900), 4), 1); + raw2_var(i,1) = var(raw((i*ns-ns+1):(i*ns-900), 4), 1); +end + +%% erase duplicate lines (the silence period of simulation) +l2 = floor(l/2) - 1; +raw3 = zeros(l2,1); +raw3_var = zeros(l2,1); +for i=0:l2, + raw3(i+1,1) = raw2((2*i+1),1); + raw3_var(i+1,1) = raw2_var((2*i+1),1); +end + +%% display graphs +nx = 4; +nt = 100; +res = zeros(nx,nx); +k =1; +for i=0:(nx-1), + for j=0:(nx-1), + cur = raw3( ((i*nx + j)*nt + 1):((i*nx + j)*nt + nt - 1), 1); + cur_var = raw3_var( ((i*nx + j)*nt + 1):((i*nx + j)*nt + nt - 1), 1); + tr = (cur(nt - 1) > 0) + 2 * (sum(cur < 0) > 0); + res(i+1,j+1) = tr; + + %if ((mod(i,5) == 0) && (mod(j,5) == 0)), + %if (tr == 2), + %if (i<=10 && j>10), + subplot(nx, nx, k); + k = k + 1; + %if (tr == 3), + plot( [ 1:0.5:(nt/2) ]', [ cur, zeros(nt-1, 1), cur-cur_var, cur+cur_var ]); + title( sprintf('i=%d, j=%d', i, j)); + %end; + %axis([-1 (nt+1) -1 1]) + set(gca,'xtick',0:10:(nt/2)) + %set(gca,'ytick',[]) + %end; + end; +end;
\ No newline at end of file diff --git a/code/matlab/random_spikes.m b/code/matlab/random_spikes.m new file mode 100644 index 0000000..89baf02 --- /dev/null +++ b/code/matlab/random_spikes.m @@ -0,0 +1,15 @@ +% fill the extenrally +%num_neurons = 2; % wich should receive a spike +%spike_freq = 1; % local per neuron per second +%duration = 2; % [s] + +current = 100; % this should drive my neurons crazy + +format long + +num = num_neurons * spike_freq * duration; +res = [ sort(rand(num, 1) .* duration), floor(rand(num, 1) .* num_neurons), ones(num, 1) * current ]; + +for i=1:length(res), + fprintf(2,'%f,%d,%f\n', res(i,1), int32(floor(res(i,2))), res(i,3)) +end diff --git a/code/matlab/random_topo.m b/code/matlab/random_topo.m new file mode 100644 index 0000000..73d94ac --- /dev/null +++ b/code/matlab/random_topo.m @@ -0,0 +1,60 @@ +%% config +% the values below should be set externally +% num_neurons = 10; +% connection_density = 0.5; +% inhibitory_fraction = 0.2; + +min_weight = 0; +max_weight = 0.004; +fie = - 1; +fei = 15 / 20; +fii = 0; + +min_delay = 0.001; % >= 0.1 ms +max_delay = 0.005; % < 2 ms + + +%% init weights + +weights = min_weight + (max_weight - min_weight) * rand(num_neurons); +weights = weights .* (rand(num_neurons) < connection_density); + +l = num_neurons - inhibitory_fraction * num_neurons + 1; +h = num_neurons; + +% make incoming and outgoing weight of each inhibitory neuron proportional +weights(1:(l-1), l:h) = weights(l:h, 1:(l-1))'; + +% scaling weights for inhibitory connections +weights(l:h, 1:(l-1)) = fie * weights(l:h, 1:(l-1)); +weights(l:h, l:h) = fii * weights(l:h, l:h); +weights(1:(l-1), l:h) = fei * weights(1:(l-1), l:h); + + +%% init delays + +delay = min_delay + (max_delay - min_delay) * rand(num_neurons); + +% make inhibitory delays shorter +delay(1:(l-1), l:h) = 0.1 * delay(1:(l-1), l:h); +delay(l:h, 1:(l-1)) = 0.1 * delay(l:h, 1:(l-1)); + + +%% print resulting topology config + +[ a, b ] = find(weights ~= 0); +index = find(weights ~= 0); + +format long + +res = [ a, b, delay(index), weights(index) ]; + +for i=1:length(res), + constness = (max(res(i,1:2)) / num_neurons >= 1 - inhibitory_fraction) % is src or dst and inhibitory neuron? + if (constness) + fprintf(2,'%d,%d,1,%f,%f\n', int32(res(i,1))-1, int32(res(i,2))-1, res(i,3), res(i,4)); + else + fprintf(2,'%d,%d,0,%f,%f\n', int32(res(i,1))-1, int32(res(i,2))-1, res(i,3), res(i,4)); + end +end + diff --git a/code/trainer/Makefile b/code/trainer/Makefile new file mode 100644 index 0000000..1fc48d1 --- /dev/null +++ b/code/trainer/Makefile @@ -0,0 +1,29 @@ +# HINT: the paradigm is not to create object files to allow all otpimizations +# take place even though the code is scattered across many files + +CC=g++ +CCFLAGS=-O3 -ggdb +LDFLAGS=-lpthread -lm +INCLUDE=-I/home/huwald/src/boost -I../core + +BASE_SRC_FILES=../core/model_switch.h ../core/fileutils.h ../core/fileutils.cpp + +.PHONY: all clean wordcount + +all: reinforce_synapse check_stdp_freq-dep + +clean: + rm -f *~ trainer massif.*.* reinforce_synapse-* check_stdp_freq-dep-* + +.PHONY: reinforce_synapse +reinforce_synapse: reinforce_synapse-dalif +reinforce_synapse-%: reinforce_synapse.cpp reinforce_synapse.h ../core/models/%.h $(BASE_SRC_FILES) + $(CC) -o $@ $(CCFLAGS) reinforce_synapse.cpp $(INCLUDE) $(LDFLAGS) -DMODEL_`echo $* | tr '[:lower:]' '[:upper:]'` + +.PHONY: check_stdp_freq-dep +check_stdp_freq-dep: check_stdp_freq-dep-dalif +check_stdp_freq-dep-%: check_stdp_freq-dep.cpp check_stdp_freq-dep.h ../core/models/%.h $(BASE_SRC_FILES) + $(CC) -o $@ $(CCFLAGS) check_stdp_freq-dep.cpp $(INCLUDE) $(LDFLAGS) -DMODEL_`echo $* | tr '[:lower:]' '[:upper:]'` + +wordcount: + wc *h *cpp Makefile diff --git a/code/trainer/check_stdp_freq-dep.cpp b/code/trainer/check_stdp_freq-dep.cpp new file mode 100644 index 0000000..f606d81 --- /dev/null +++ b/code/trainer/check_stdp_freq-dep.cpp @@ -0,0 +1,160 @@ +#include <stdlib.h> +#include "fileutils.h" +#include "math.h" +#include "unistd.h" + +#include "check_stdp_freq-dep.h" +#include "fileutils.cpp" +#include "model_switch.h" + +using namespace std; + +int main(int argc, char **argv) { + // check cmd line sanity + if (argc != 5) { + fprintf(stderr, "Wrong argument count\n\n" + "Call format:\n" + "%s\n\t" + "performance out\n\t" + "trace cmd out\n\t" + "global out\n\t" + "spike out\n\t" + "\n" + "Special names allowed:\n\t- (standart input)\n\t0 (/dev/null)\n", argv[0]); + return -1; + } + + Trainer *t = new Trainer(argc, argv); + t->run(); + + pthread_join(t->thread_write, NULL); +} + +Trainer::Trainer(int argc, char** argv) { + // init vars + currentEpoch = 0; + epochDuration = 10.0; // [s] + neurons = 1000; // number of neurons to send noise to + voltage = 0.1; // [V] + md = 1.2; // 10.0; + mss = 1.1; //1.2; + fs = 100; // number of frequencies to try + frd = 1.0; // relative difference between to frequencies (f_i+1 = frd * f_i) + fad = 0.5; // absolute difference between to frequencies (f_i+1 = fad + f_i) + + // open all file descriptors in an order complementary to the simulators one + // to avoid deadlocks + fd_spike_out = fd_magic(argv[4], true); + fd_global_out = fd_magic(argv[3], true); + fd_performance_out = fd_magic(argv[1], true); + fd_trace_out = fd_magic(argv[2], true); + + // create read and write threads + pthread_create(&thread_write, NULL, (void* (*)(void*)) &write_spikes, this); +} + +void Trainer::run() { + char *str_trace = "%f; synapse\n"; + + // init global sim variables + MS_Global msg; + msg_init(msg); + msg.dopamin_level = 0.0; + + double ta = 0.009821, //0.0088541, + la = 0.140249; // 0.126445; + + /* + // loop over both vars to examine + for (msg.stdp_tau_plus = msg.stdp_tau_minus / md; + msg.stdp_tau_plus <= msg.stdp_tau_minus * md; + msg.stdp_tau_plus *= mss) { + for (msg.stdp_lambda_plus = msg.stdp_lambda_minus / md; + msg.stdp_lambda_plus <= msg.stdp_lambda_minus * md; + msg.stdp_lambda_plus *= mss) {*/ + // loop over both vars to examine + for (msg.stdp_tau_plus = ta / md; + msg.stdp_tau_plus <= ta * md; + msg.stdp_tau_plus *= mss) { + for (msg.stdp_lambda_plus = la / md; + msg.stdp_lambda_plus <= la * md; + msg.stdp_lambda_plus *= mss) { + + // print the parameters to the performance output + msg_print(msg, fd_performance_out); + fprintf(fd_performance_out, "\n"); + + // print the global params + fprintf(fd_global_out, "%f, ", currentEpoch * epochDuration); + msg_print(msg, fd_global_out); + fprintf(fd_global_out, "\n"); + + // let the simulation proceed + fprintf(fd_trace_out, str_trace, epochDuration); + currentEpoch++; + + // repeat this 2*n-1 times (n=number of different frequency trials) + for (int i=0; i < 2*fs-1; i++) { + fprintf(fd_trace_out, "\n"); + currentEpoch++; + } + } + } + + fclose(fd_trace_out); + fclose(fd_global_out); +} + +// ---- send indepenent poisson noise w/ increasing fequency---- +void *write_spikes(Trainer *t) { + // calculate how often we have to try all frequencies (=outer loop) + // WARN: ignore minor numerical instabilities + int max = (int) floor(2.0 * log(t->md) / log(t->mss) ) + 1; + max *= max; // there are two nested loops of the same size + + double time = 0.0; // global time (that one send to the simulator) + + // for each paramter config (set in the main routine) + for (int i=0; i<max; i++) { + + double freq = 1.0; + + // examine a set of frequencies + for (int j=0; j < t->fs; j++) { + // send out the spikes + double localtime = 0.0; + double nextRefSpike = 0.0; + double refFreq = 10.0; // [Hz] + int dst = -1; + while (localtime < t->epochDuration) { + // starting with the second call ... + if (dst != -1) { + // check if we have to send a spike to the ref neuron + if (localtime > nextRefSpike) { + fprintf(t->fd_spike_out, "%f, %d, %f\n", time + nextRefSpike, 0, t->voltage); + nextRefSpike += 1.0 / refFreq; + } + + // send spike to the simulator + fprintf(t->fd_spike_out, "%f, %d, %f\n", time + localtime, dst, t->voltage); + }else{ + } + + localtime -= log(1.0 - drand48()) / (freq * t->neurons); // possion distributed spike timing + dst = 1 + rand() % (t->neurons - 1); // random neuron (except reference neuron 0) + } + + // increase time (twice because of the silence period after each noise period) + time = (i * t->fs + j) * 2.0 * t->epochDuration; + + // increase frequency + freq *= t->frd; + freq += t->fad; + } + } + + // close fd because fscanf sucks + fclose(t->fd_spike_out); + + return NULL; +} diff --git a/code/trainer/check_stdp_freq-dep.h b/code/trainer/check_stdp_freq-dep.h new file mode 100644 index 0000000..cf429b6 --- /dev/null +++ b/code/trainer/check_stdp_freq-dep.h @@ -0,0 +1,46 @@ +#ifndef TRAINER_H +#define TRAINER_H + +#include <stdio.h> +#include <pthread.h> +#include <map> +#include <queue> +#include "boost/tuple/tuple.hpp" + + +using namespace std; + +class Trainer { + public: + FILE *fd_spike_out, + *fd_global_out, + *fd_trace_out, + *fd_performance_out; + + // init stuff + Trainer(int argc, char** argv); + + // main routine + void run(); + + // state vars + long currentEpoch; + + // thread related + pthread_t thread_write; + + // configuration + double md; // multiplicative difference (>1) + double mss; // multiplicative step size (>1) + double fs, frd, fad; // number of frequency steps and relative step size + double epochDuration; + double voltage; // per outgoing (random) spike + long neurons; +}; + +// seperate thread to read all spikes neccessary because reading and +// writing to these descriptors could block and thus cause a deadlock +void *read_spikes(Trainer *t); +void *write_spikes(Trainer *t); + +#endif // TRAINER_H diff --git a/code/trainer/mem1.cpp b/code/trainer/mem1.cpp new file mode 100644 index 0000000..3b522b4 --- /dev/null +++ b/code/trainer/mem1.cpp @@ -0,0 +1,412 @@ +#include <stdlib.h> +#include "fileutils.h" +#include "math.h" + +#include "reinforce_synapse.h" +#include "fileutils.cpp" +#include "model_switch.h" + +using namespace std; + +int main(int argc, char **argv) { + // check cmd line sanity + if (argc != 7) { + fprintf(stderr, "Wrong argument count\n\n" + "Call format:\n" + "%s\n\t" + "performance out\n\t" + "trace cmd out\n\t" + "global out\n\t" + "global in\n\t" + "spike out\n\t" + "spike in\n\t" + "\n" + "Special names allowed:\n\t- (standart input)\n\t0 (/dev/null)\n", argv[0]); + return -1; + } + + Trainer *t = new Trainer(argc, argv); + t->run(); +} + +//===== Initialisation ===================================== + +Trainer::Trainer(int argc, char** argv) { + initConfig(); + initState(); + initGroups(); // determine input and output neurons + + initFiles(); + initThreads(); +} + +void Trainer::initConfig() { + neurons = 1000; + neuronsPerSymbol = 200; + noiseFreq = 10.0; // [Hz] + noiseVoltage = 0.03; // [V] + reward = 0.1; + + epochDuration = 1.0; // [s] + numTrials = 1000; + numSymbols = 2; + //readoutDelay = 1; + refractoryPeriods = 3; +} + +void Trainer::initState() { + dopamin_level = 0.0; + currentTrial = 0; + state = 0; + + msg_init(msg); + msg.dopamin_level = dopamin_level; + + groupFreq.resize(numSymbols); + for (int i=0; i<numSymbols; i++) { + groupFreq[i] = 0; + } +} + +void Trainer::initGroups() { + ioNeurons[i] = new set<int>(); + for (int j=0; j<neuronsPerSymbol;) { + int n = rand() % numNeurons; + if (!isNeurons[i]->count(n)) { + ioNeurons[i]->insert(n); + j++; + } + } +} + +void Trainer::initFiles() { + // open all file descriptors in an order complementary to the simulators one + // to avoid deadlocks + fd_spike_in = fd_magic(argv[6], false); + fd_global_in = fd_magic(argv[4], false); + fd_spike_out = fd_magic(argv[5], true); + fd_global_out = fd_magic(argv[3], true); + fd_performance_out = fd_magic(argv[1], true); + fd_trace_out = fd_magic(argv[2], true); +} + +void Trainer::initThreads() { + // init locks + pthread_mutex_init(&incomingSpikeLock, NULL); + pthread_mutex_init(&writerLock, NULL); + + // create read and write threads + pthread_create(&thread_read, NULL, (void* (*)(void*)) &read_spikes, this); + pthread_create(&thread_write, NULL, (void* (*)(void*)) &write_spikes, this); +} + +//===== Core trainer ==================================== + +void Trainer::pushGlobal(double time) { + fprintf(fd_global_out, "%f, ", time); + msg_print(msg, fd_global_out); + fprintf(fd_global_out, "\n"); + fflush(fd_global_out); +} + +// hint time is delta time! +void Trainer::pushTrace(double time) { + const char *str_trace = "%f; spikes (0; 1); global; neuron (0; 1); synapse (0; 1)\n"; + fprintf(fd_trace_out, str_trace, time); + fflush(fd_trace_out); +} + +bool Trainer::readGlobal() { + double _foo_dbl; + char str_raw[128], + str_msg[128]; + str_raw[0] = 0; + + // read a single line + if (fgets((char*) str_raw, 128, fd_global_in) == NULL) { + fprintf(stderr, "ERROR: global status file descriptor from simulator closed unexpectedly\n"); + return false; + } + + // parse it + if ((sscanf((char*) str_raw, "%lf, %[^\n]\n", &_foo_dbl, (char*) str_msg) != 2) + || (!msg_parse(msg, (char*) str_msg))) { + fprintf(stderr, "ERROR: reading global status from simulator failed\n\t\"%s\"\n", (char*) str_raw); + return false; + } + + return true; +} + +void binIncomingSpikes() { + // reset bins + for (int i=0; i<groupFreq.size(); i++) + groupFreq[i] = 0; + + // lock spike queue + pthread_yield(); // give the spike reading thread chance to finish ... this is not more than ugly semifix wrong par! + pthread_mutex_lock(&incomingSpikeLock); + + // read all spikes in the correct time window + while ((!incomingSpikes.empty()) && (incomingSpikes.front().get<0>() <= currentEpoch * epochDuration)) { + // drop event out of queue + SpikeEvent se = incomingSpikes.front(); + double time = se.get<0>(); + int neuron = se.get<1>(); + incomingSpikes.pop(); + + // check if it belongs to the previous bin (and ignore it if this is the case) + if (time < (currentEpoch - 1) * epochDuration) { + fprintf(stderr, "WARN: spike reading thread to slow; unprocessed spike of the past discovered\n%f\t%f\t%d\t%f\n", + time, (double) (currentEpoch - 1) * epochDuration, currentEpoch, epochDuration); + continue; + } + + // check membership in each group and increase group frequency + for (int i=0; i < ioNeurons.size(); i++) + if (ioNeuros[i]->count(neuron)) + groupFreq[i]++; + } + + pthread_mutex_unlock(&incomingSpikeLock); +} + +void Trainer::addBaselineSpikes() { +} + +void Trainer::addSymbolSpikes() { + +} + + +double Trainer::calcSignalStrength() { + if (symbolHist.empty()) { + fprintf(stderr, "Writer thread is too slow; missed the current symbol\n"); + exit(-1); + } + + int fs, fn = 0; // freq signal, freq noise + for (int i=0; i<numSymbols; i++) { + if (i == symbolHist.front()) { + fs = groupFreq[i]; + }else{ + fm = fmax(fm, groupFreq[i]); + } + } + + if (fn == 0) { + return fs * INFINITY; + }else{ + return ((double) fs) / fn; + } +} + +void Trainer::run() { + // rough description of this function + // . start an epoch + // . wait for it's end + // . process incomig spikes (binning) + // . select if a reward takes place + // . print reward value + // . send out the reward signal + + // send out the full trace command once (later it will be repeated by sending newline) + pushTrace(epochDuration); + + // send the first two global states (at t=0 and t=1.5 [bintime] to allow the simulation to + // be initialized (before the causality of the loop below is met) + pushGlobal(0.0); + msg_process(msg, 1.5 * epochDuration); + dopamin_level = msg.dopamin_level; + pushGlobal(1.5 * epochDuration); + + // loop until the experiment is done + for (; currentEpoch * epochDuration < entireDuration; currentEpoch++) { + + // send a new trace command (do it as early as possible although it is + // only executed after the new global is send out at the bottom of this loop) + if ((currentEpoch + 2) * epochDuration < entireDuration) { + // repeat the previous trace command + fprintf(fd_trace_out, "\n"); + fflush(fd_trace_out); + }else{ + pushTrace(entireDuration - (currentEpoch + 1) * epochDuration); + } + + // send new spikes + pthread_mutex_lock(&outgoingSpikeLock); + addBaselineSpikes(); + if (state == 0) addSymbolSpikes(); + pthread_cond_signal(&outgoingSpikeCond); + pthread_mutex_unlock(&outgoingSpikeLock); + + // wait for the end of the epoch (by reading the global state resulting from it) + if (!readGlobal()) + break; + + // process incomig spikes (binning) of the previous epoch + if (currentEpoch > 0) + binIncomingSpikes(); + + // proceed the global state to keep it in sync with the simulator's global state + // the local dopamin level is kept seperately and aged only one epochDuration to + // avoid oscillation effects in dopamin level + msg_process(msg, 1.5 * epochDuration); + dopamin_level *= exp( - epochDuration / msg.dopamin_tau ); + + // do various actions depeding on state (thus lock mutex of the writer thread) + + + switch (state) { + case 0: // a signal is sent + state++; + + case 1: // we are waiting for the signal to be reproduce + // get fraction of the current symbol's freq compared to the strongest wrong symbol + double ss = calcSignalStrength(); + + // check if the reward condition is met + if (ss > 1) { + dopmain_level += reward; + }else{ + state++; // lost signal -> next state (and finally a new trial) + currentSymbol = rand() % numSymbols; // determine new symbol to display + } + + break; + + default: // the signal has been lost (in the last round); refractory time + ++state %= refractoryPeriods; + } + + /*if ((currentEpoch > 1) && ((*neuronFreq[0])[0] > 0) && ((*neuronFreq[1])[1] > 0)) { + dopamin_level += da_single_reward; + fprintf(fd_performance_out, "+"); + }else{ + fprintf(fd_performance_out, "-"); + }*/ + + // performance and "debug" output + if (currentEpoch > 1) { + //fprintf(fd_performance_out, "\n"); + fprintf(fd_performance_out, "\t%f\t%d\t%d\n", dopamin_level, (*neuronFreq[0])[0], (*neuronFreq[1])[1]); + }else{ + // fake output as acutal data is not available, yet + fprintf(fd_performance_out, "\t%f\t%d\t%d\n", dopamin_level, (int) 0, (int) 0); + } + + // set the new DA level + msg.dopamin_level = dopamin_level; + + // print new global state + // (do this even if there has been no evaluation of the performance yet, + // because it is neccessary for the simulator to proceed) + pushGlobal(((double) currentEpoch + 2.5) * epochDuration); + } + + fclose(fd_trace_out); + + // terminate child threads + pthread_cancel(thread_read); + pthread_cancel(thread_write); +} + +void *read_spikes(Trainer *t) { + double lastSpike = -INFINITY; // used to check if the spikes are coming in order + + // read spikes until eternity + while (!feof(t->fd_spike_in)) { + // read one line from stdin (blocking) + char buf[128]; + if (fgets((char*) buf, 128, t->fd_spike_in) == NULL) continue; // this should stop the loop because of EOF + + // parse the input + double time, current; + int neuron; + switch (sscanf((char*) buf, "%lf, %d, %lf\n", &time, &neuron, ¤t)) { + case 3: + // format is ok, continue + break; + default: + // format is wrong, stop + fprintf(stderr, "ERROR: malformatted incoming spike:\n\t%s\n", &buf); + return NULL; + } + + if (lastSpike > time) { + fprintf(stderr, "WARN: out of order spike detected (coming from simulator)\n\t%f\t%d\n", time, neuron); + continue; + } + + lastSpike = time; + + // add the spike to the queue of spikes + pthread_mutex_lock(&(t->incomingSpikeLock)); + t->incomingSpikes.push(boost::make_tuple(time, neuron, current)); + pthread_mutex_unlock(&(t->incomingSpikeLock)); + } + + // we shouldn't reach this point in a non-error case + fprintf(stderr, "ERROR: EOF in incoming spike stream\n"); + // TODO: kill entire programm + return NULL; +} + +void *write_spikes(Trainer *t) { + // at the moment: generate noise until the file descriptor blocks + double time = 0.0; + + // PAR HINT: + // loop until exactly one spike after the entire duration is send out + // this will block on full buffer on the file descriptor and thus keep + // the thread busy early enough + + + /* // ---- send 100% dependent spike train --- + time = 0.005; + while (time <= t->entireDuration) { + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, 0, 1.0); + time += 0.012; + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, 1, 1.0); + time += 1.0; + }*/ + + + /* // ---- send indepenent poisson noise ---- + while (time <= t->entireDuration) { + // calc timing, intensity and destination of the spike + // HINT: + // * log(...) is negative + // * drand48() returns something in [0,1), to avoid log(0) we transform it to (0,1] + time -= log(1.0 - drand48()) / (t->freq * t->neurons); + int dst = rand() % t->neurons; + double current = t->voltage; + + // send it to the simulator + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, dst, current); + }*/ + + // ---- send indepenent poisson noise w7 increasing fequency---- + double blafoo = 0; + t->freq = 1.0; + while (time <= t->entireDuration) { + if (time - blafoo > 100.0) { + blafoo += 200.0; + t->freq += 1.0; + time += 100.0; // time jump to let ET recover to zero + } + // calc timing, intensity and destination of the spike + // HINT: + // * log(...) is negative + // * drand48() returns something in [0,1), to avoid log(0) we transform it to (0,1] + time -= log(1.0 - drand48()) / (t->freq * t->neurons); + int dst = rand() % t->neurons; + double current = t->voltage; + + // send it to the simulator + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, dst, current); + } + + // close fd because fscanf sucks + fclose(t->fd_spike_out); +} diff --git a/code/trainer/mem1.h b/code/trainer/mem1.h new file mode 100644 index 0000000..31fc2b0 --- /dev/null +++ b/code/trainer/mem1.h @@ -0,0 +1,72 @@ +#ifndef TRAINER_H +#define TRAINER_H + +#include <stdio.h> +#include <pthread.h> +#include <map> +#include <queue> +#include "boost/tuple/tuple.hpp" + + +using namespace std; + +class Trainer { + public: + FILE *fd_spike_in, + *fd_spike_out, + *fd_global_out, + *fd_global_in, + *fd_trace_out, + *fd_performance_out; + + // init stuff + Trainer(int argc, char** argv); + void initConfig(); + void initState(); + void initGroups(); + void initFiles(); + void initThreads(); + + // main routine + void run(); + + // state vars + long currentTrial; + vector<int> groupFreq; // spikes fired during epoch for each symbol (=group of neurons) + vector< set<int>* > ioNeurons; // a set of neurons for each symbol to which the symbol is written and from wich it is read later + //queue< int > symbolHist; // stores the symbols displayed; written by thread_write, read by main thread + int currentSymbol; + double dopamin_level; + int state; + MS_Global msg; + + // synchronisation vars + typedef boost::tuple<double, int, double> SpikeEvent; // <what time, wich neuron, current> + queue<SpikeEvent> incomingSpikes; + // TODO: outgoingSpikes; + pthread_mutex_t incomingSpikeLock, writerLock; + // TODO: , outgoingSpikeLock; (including a condition for the writer to wakeup upon if a previously empyt queue has been filled) + pthread_t thread_read, thread_write; + + // configuration + // network + long neurons; // total number of neurons + long neuronsPerSymbol; + double noiseFreq; // of poisson noise (per neuron) + double noiseVoltage; // per noise spike + double reward; // per succesful trial + + // learning task + double epochDuration; // a trial consists of several epochs + long numTrials; + int numSymbols; // number of different things to remember + //int readoutDelay; // number of epochs between symbol-write and symbol-read-epoch + int refractoryPeriods; // how many epochs to wait after a trial finished until we start with a new trial +}; + +// seperate thread to read all spikes neccessary because reading and +// writing to these descriptors could block and thus cause a deadlock +void *read_spikes(Trainer *t); +void *write_spikes(Trainer *t); + +#endif // TRAINER_H diff --git a/code/trainer/reinforce_synapse.cpp b/code/trainer/reinforce_synapse.cpp new file mode 100644 index 0000000..bf6fc7f --- /dev/null +++ b/code/trainer/reinforce_synapse.cpp @@ -0,0 +1,302 @@ +#include <stdlib.h> +#include "fileutils.h" +#include "math.h" + +#include "reinforce_synapse.h" +#include "fileutils.cpp" +#include "model_switch.h" + +using namespace std; + +int main(int argc, char **argv) { + // check cmd line sanity + if (argc != 7) { + fprintf(stderr, "Wrong argument count\n\n" + "Call format:\n" + "%s\n\t" + "performance out\n\t" + "trace cmd out\n\t" + "global out\n\t" + "global in\n\t" + "spike out\n\t" + "spike in\n\t" + "\n" + "Special names allowed:\n\t- (standart input)\n\t0 (/dev/null)\n", argv[0]); + return -1; + } + + Trainer *t = new Trainer(argc, argv); + t->run(); + // TODO: finalize +} + +Trainer::Trainer(int argc, char** argv) { + // init vars + currentEpoch = 0; + dopamin_level = 0.0; + + epochDuration = 0.01; // [s] + //epochDuration = 1.0; // [s] + entireDuration = 20000.0; // [s] + neurons = 2; // number of neurons to send noise to + freq = 1.0; // [Hz] per Neuron + voltage = 0.1; // [V] + da_single_reward = 0.01; + + neuronFreq[0] = (map<int, int>*) NULL; + neuronFreq[1] = (map<int, int>*) NULL; + + // open all file descriptors in an order complementary to the simulators one + // to avoid deadlocks + fd_spike_in = fd_magic(argv[6], false); + fd_global_in = fd_magic(argv[4], false); + fd_spike_out = fd_magic(argv[5], true); + fd_global_out = fd_magic(argv[3], true); + fd_performance_out = fd_magic(argv[1], true); + fd_trace_out = fd_magic(argv[2], true); + + // init locks + pthread_mutex_init(&incomingSpikeLock, NULL); + + // create read and write threads + pthread_create(&thread_read, NULL, (void* (*)(void*)) &read_spikes, this); + pthread_create(&thread_write, NULL, (void* (*)(void*)) &write_spikes, this); +} + +void Trainer::run() { + // start an epoch + // wait for it's end + // process incomig spikes (binning) + // select if a reward takes place + // print reward value (TODO: into a seperate, externally given file descriptor) + // send out the reward signal + + char *str_trace = "%f; spikes (0; 1); global; neuron (0; 1); synapse (0; 1)\n"; + + // send out the full trace commande once (later it will be repeated by sending newline) + fprintf(fd_trace_out, str_trace, epochDuration); + fflush(fd_trace_out); + + // send the first two global states (at t=0 and t=1.5 [bintime] to allow the simulation to + // be initialized (before the causality of the loop below is met) + MS_Global msg; + msg_init(msg); + msg.dopamin_level = dopamin_level; + + // set the tau-levels like in Izhi's network + msg.stdp_tau_minus = 1.5 * msg.stdp_tau_plus; + msg.stdp_lambda_plus = msg.stdp_lambda_minus; + + fprintf(fd_global_out, "0.0, "); + msg_print(msg, fd_global_out); + fprintf(fd_global_out, "\n"); + + msg_process(msg, 1.5 * epochDuration); + dopamin_level = msg.dopamin_level; + + + fprintf(fd_global_out, "%f, ", 1.5 * epochDuration); + msg_print(msg, fd_global_out); + fprintf(fd_global_out, "\n"); + + fflush(fd_global_out); + + // loop until the experiment is done + for (; currentEpoch * epochDuration < entireDuration; currentEpoch++) { + // send a new trace command (do it as early as possible although it is + // only executed after the new global is send out at the bottom of this loop) + if ((currentEpoch + 2) * epochDuration < entireDuration) { + // repeat the previous trace command + fprintf(fd_trace_out, "\n"); + }else{ + fprintf(fd_trace_out, str_trace, entireDuration - (currentEpoch + 1) * epochDuration); + } + fflush(fd_trace_out); + + // wait for the end of the epoch (by reading the global state resulting from it) + char str_raw[128], str_msg[128]; str_raw[0] = 0; + double _foo_dbl; + if (fgets((char*) str_raw, 128, fd_global_in) == NULL) { + fprintf(stderr, "ERROR: global status file descriptor from simulator closed unexpectedly\n"); + break; + } + if ((sscanf((char*) str_raw, "%lf, %[^\n]\n", &_foo_dbl, (char*) str_msg) != 2) + || (!msg_parse(msg, (char*) str_msg))) { + fprintf(stderr, "ERROR: reading global status from simulator failed\n\t\"%s\"\n", (char*) str_raw); + break; + } + + // process incomig spikes (binning) of the previous epoch + if (currentEpoch > 0) { + // shift the bins + if (neuronFreq[0]) { + delete neuronFreq[0]; + neuronFreq[0] = neuronFreq[1]; + }else{ + neuronFreq[0] = new map<int, int>(); + } + neuronFreq[1] = new map<int, int>(); + + // read all spikes in the correct time window + pthread_mutex_lock(&incomingSpikeLock); + while ((!incomingSpikes.empty()) && (incomingSpikes.front().get<0>() <= currentEpoch * epochDuration)) { + // drop event out of queue + SpikeEvent se = incomingSpikes.front(); + double time = se.get<0>(); + int neuron = se.get<1>(); + incomingSpikes.pop(); + + // check if it belongs to the previous bin (and ignore it if this is the case) + if (time < (currentEpoch - 1) * epochDuration) { + fprintf(stderr, "WARN: spike reading thread to slow; unprocessed spike of the past discovered\n%f\t%f\t%d\t%f\n", time, (double) (currentEpoch - 1) * epochDuration, currentEpoch, epochDuration); + continue; + } + + // increment the frequency counter (relies on int being default constructable to value 0) + (*neuronFreq[1])[neuron]++; + } + + pthread_mutex_unlock(&incomingSpikeLock); + } + + // proceed the global state to keep it in sync with the simulator's global state + // the local dopamin level is kept seperately and aged only one epochDuration to + // avoid oscillation effects in dopamin level + msg_process(msg, 1.5 * epochDuration); + dopamin_level *= exp( - epochDuration / msg.dopamin_tau ); + + // select if the reward takes place + if ((currentEpoch > 1) && ((*neuronFreq[0])[0] > 0) && ((*neuronFreq[1])[1] > 0)) { + dopamin_level += da_single_reward; + fprintf(fd_performance_out, "+"); + }else{ + fprintf(fd_performance_out, "-"); + } + + if (currentEpoch > 1) { + //fprintf(fd_performance_out, "\n"); + fprintf(fd_performance_out, "\t%f\t%d\t%d\n", dopamin_level, (*neuronFreq[0])[0], (*neuronFreq[1])[1]); + }else{ + // fake output as acutal data i not available, yet + fprintf(fd_performance_out, "\t%f\t%d\t%d\n", dopamin_level, (int) 0, (int) 0); + } + + // set the new DA level + msg.dopamin_level = dopamin_level; + + // print new global state + // (do this even if there has been no evaluation of the performance yet, + // because it is neccessary for the simulator to proceed) + + fprintf(fd_global_out, "%f, ", ((double) currentEpoch + 2.5) * epochDuration); + msg_print(msg, fd_global_out); + fprintf(fd_global_out, "\n"); + fflush(fd_global_out); + } + + fclose(fd_trace_out); + + // terminate child threads + pthread_cancel(thread_read); + pthread_cancel(thread_write); +} + +void *read_spikes(Trainer *t) { + double lastSpike = -INFINITY; // used to check if the spikes are coming in order + + // read spikes until eternity + while (!feof(t->fd_spike_in)) { + // read one line from stdin (blocking) + char buf[128]; + if (fgets((char*) buf, 128, t->fd_spike_in) == NULL) continue; // this should stop the loop because of EOF + + // parse the input + double time, current; + int neuron; + switch (sscanf((char*) buf, "%lf, %d, %lf\n", &time, &neuron, ¤t)) { + case 3: + // format is ok, continue + break; + default: + // format is wrong, stop + fprintf(stderr, "ERROR: malformatted incoming spike:\n\t%s\n", &buf); + return NULL; + } + + if (lastSpike > time) { + fprintf(stderr, "WARN: out of order spike detected (coming from simulator)\n\t%f\t%d\n", time, neuron); + continue; + } + + lastSpike = time; + + // add the spike to the queue of spikes + pthread_mutex_lock(&(t->incomingSpikeLock)); + t->incomingSpikes.push(boost::make_tuple(time, neuron, current)); + pthread_mutex_unlock(&(t->incomingSpikeLock)); + } + + // we shouldn't reach this point in a non-error case + fprintf(stderr, "ERROR: EOF in incoming spike stream\n"); + // TODO: kill entire programm + return NULL; +} + +void *write_spikes(Trainer *t) { + // at the moment: generate noise until the file descriptor blocks + double time = 0.0; + + // PAR HINT: + // loop until exactly one spike after the entire duration is send out + // this will block on full buffer on the file descriptor and thus keep + // the thread busy early enough + + + /* // ---- send 100% dependent spike train --- + time = 0.005; + while (time <= t->entireDuration) { + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, 0, 1.0); + time += 0.012; + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, 1, 1.0); + time += 1.0; + }*/ + + + // ---- send indepenent poisson noise ---- + while (time <= t->entireDuration) { + // calc timing, intensity and destination of the spike + // HINT: + // * log(...) is negative + // * drand48() returns something in [0,1), to avoid log(0) we transform it to (0,1] + time -= log(1.0 - drand48()) / (t->freq * t->neurons); + int dst = rand() % t->neurons; + double current = t->voltage; + + // send it to the simulator + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, dst, current); + } + + /*// ---- send indepenent poisson noise w7 increasing fequency---- + double blafoo = 0; + t->freq = 1.0; + while (time <= t->entireDuration) { + if (time - blafoo > 100.0) { + blafoo += 200.0; + t->freq += 1.0; + time += 100.0; // time jump to let ET recover to zero + } + // calc timing, intensity and destination of the spike + // HINT: + // * log(...) is negative + // * drand48() returns something in [0,1), to avoid log(0) we transform it to (0,1] + time -= log(1.0 - drand48()) / (t->freq * t->neurons); + int dst = rand() % t->neurons; + double current = t->voltage; + + // send it to the simulator + fprintf(t->fd_spike_out, "%f, %d, %f\n", time, dst, current); + }*/ + + // close fd because fscanf sucks + fclose(t->fd_spike_out); +} diff --git a/code/trainer/reinforce_synapse.h b/code/trainer/reinforce_synapse.h new file mode 100644 index 0000000..46b0083 --- /dev/null +++ b/code/trainer/reinforce_synapse.h @@ -0,0 +1,55 @@ +#ifndef TRAINER_H +#define TRAINER_H + +#include <stdio.h> +#include <pthread.h> +#include <map> +#include <queue> +#include "boost/tuple/tuple.hpp" + + +using namespace std; + +class Trainer { + public: + FILE *fd_spike_in, + *fd_spike_out, + *fd_global_out, + *fd_global_in, + *fd_trace_out, + *fd_performance_out; + + // init stuff + Trainer(int argc, char** argv); + + // main routine + void run(); + + // state vars + long currentEpoch; + map<int, int> *neuronFreq[2]; // stores if a surveilled neuron fired during the current or last epoch + double dopamin_level; + + // synchronisation vars + typedef boost::tuple<double, int, double> SpikeEvent; // <what time, wich neuron, current> + queue<SpikeEvent> incomingSpikes; + // TODO: outgoingSpikes; + pthread_mutex_t incomingSpikeLock; + // TODO: , outgoingSpikeLock; (including a condition for the writer to wakeup upon if a previously empyt queue has been filled) + pthread_t thread_read, thread_write; + + // configuration + double epochDuration; + double entireDuration; + double freq; // of outgoing noise per neuron + double voltage; // per outgoing (random) spike + long neurons; + double da_single_reward; +}; + +// seperate thread to read all spikes neccessary because reading and +// writing to these descriptors could block and thus cause a deadlock +void *read_spikes(Trainer *t); +void *write_spikes(Trainer *t); + +#endif // TRAINER_H diff --git a/code/trainer/test.cpp b/code/trainer/test.cpp new file mode 100644 index 0000000..ceb2484 --- /dev/null +++ b/code/trainer/test.cpp @@ -0,0 +1,13 @@ +#include <stdio.h> +#include <math.h> + +int main() { + fprintf("aaa\n"); + while (true);/* + double d=0; + printf("aa %lf\n", d); + d = 0; + scanf("%lf\n", &d); + printf("%lf\n", d);*/ + return 0; +} |