diff --git a/.gitignore b/.gitignore index 204b2f1b..be1d61b2 100644 --- a/.gitignore +++ b/.gitignore @@ -37,5 +37,7 @@ m4/ltsugar.m4 m4/ltversion.m4 m4/lt~obsolete.m4 +cscope.out +tags # travis build file .travis.yml diff --git a/AUTHORS b/AUTHORS index da1e658f..a315aebc 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,6 +1,9 @@ -Package created by Friedemann Zenke +Auryn package created by Friedemann Zenke -Contributors Lorric Ziegler, Ankur Sinha +Contributors: +Lorric Ziegler +Ankur Sinha +Emre Neftci diff --git a/ChangeLog b/ChangeLog index cf459614..a616e531 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,31 @@ -xxxx-xx-xx Friedemann Zenke +2015-07-01 Friedemann Zenke + * Added PairInteractionConnection to implement arbitrary (pairwise) + STDP windows + * Added backend to implement current based synapses with temporal dynamics + * Added code to stimulate with currents generated from a Wiener process + (NormalStimulator) + * Added the AdEx neuron model + * Added the SIFGroup neuron model + * Implemented load/save to netstate files of Checkers + * Changed backend of WeightMonitor to avoid segfaults after loading the + network state + * Implementation of moving average for progress bar moved + to System class to disambiguate from Checkers + * Fixed some problems in ComplexMatrix + * Bug fixes + +2015-02-06 Friedemann Zenke + * SpikingGroups and Connections are now serializeable and as a consequence + the network state can be saved to a single file per rank. + * SimpleMatrix has been replaced by ComplexMatrix which creates the + basis for the implementation of connection models with multiple + internal synaptic states. + * Basis for Monitors writing to binary files for increased performance was + created (e.g. BinarySpikeMonitor). + * Auryn compiles to a statically linkable library per default now, which + facilitates to decouple simulation code from simulator code. * A simplified + STDPConnection was added along with a tutorial of how to implement own new + plasticity and synapse models. * Adds integrate-and-fire model for exclusive use with current based synapses (CubaIFGroup). * Adds an example simulation for short-term-plasticity (STPConnection). diff --git a/README.md b/README.md index 8ab67871..e2723622 100644 --- a/README.md +++ b/README.md @@ -10,11 +10,23 @@ to simulate recurrent spiking neural networks with spike timing dependent plasticity (STDP). It comes with the GPLv3 (please see COPYING). +Quick start +----------- + +To download and compile the examples try: + +``` +sudo apt-get install cmake git build-essential libboost-all-dev +git clone https://github.com/fzenke/auryn.git && cd auryn +git checkout -b develop origin/develop +mkdir build && cd build +cmake ../ -DCMAKE_BUILD_TYPE=Release && make +``` + Documentation & Installation/Use -------------------------------- -Please visit http://www.fzenke.net/auryn/ - +Please visit the wiki at http://www.fzenke.net/auryn/ Requirements diff --git a/examples/sim_background.cpp b/examples/sim_background.cpp index d7b6cda6..d43dbc4c 100644 --- a/examples/sim_background.cpp +++ b/examples/sim_background.cpp @@ -615,6 +615,7 @@ int main(int ac, char* av[]) for ( int i = 0 ; i < 5 ; ++i ) { for ( int j = 0 ; j < 5 ; ++j ) { vector sublist = con_ee->get_block(i*psize,(i+1)*psize,j*psize,(j+1)*psize); + sublist.resize(50); // only record a maximum of 50 connections from each block wmon->add_to_list(sublist); } } @@ -645,6 +646,10 @@ int main(int ac, char* av[]) RateChecker * chk = new RateChecker( neurons_e , 0.1 , 20.*kappa , tau_chk); + // Use the same time constant for the online rate estimate in the progress bar + sys->set_online_rate_monitor_id(0); + sys->set_online_rate_monitor_tau(tau_chk); + if ( scaling && (errcode==0) ) { diff --git a/src/AdExGroup.cpp b/src/AdExGroup.cpp new file mode 100644 index 00000000..42d5b647 --- /dev/null +++ b/src/AdExGroup.cpp @@ -0,0 +1,224 @@ +/* +* Copyright 2014-2015 Ankur Sinha and Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +* +* If you are using Auryn or parts of it for your work please cite: +* Zenke, F. and Gerstner, W., 2014. Limits to high-speed simulations +* of spiking neural networks using general-purpose computers. +* Front Neuroinform 8, 76. doi: 10.3389/fninf.2014.00076 +*/ + +#include "AdExGroup.h" + +AdExGroup::AdExGroup(NeuronID size) : NeuronGroup(size) +{ + sys->register_spiking_group(this); + if ( evolve_locally() ) init(); +} + +void AdExGroup::calculate_scale_constants() +{ + scale_mem = dt/tau_mem; + scale_w = dt/tau_w; + scale_ampa = exp(-dt/tau_ampa); + scale_gaba = exp(-dt/tau_gaba); +} + +void AdExGroup::init() +{ + g_leak = 10e-9; + e_rest = -70e-3; + e_reset = -58e-3; + e_rev_ampa = 0; + e_rev_gaba = -80e-3; + e_thr = -50e-3; + tau_ampa = 5e-3; + tau_gaba = 10e-3; + tau_w = 30e-3; + c_mem = 200e-12; + tau_mem = c_mem/g_leak; + deltat = 2e-3; + a = 2e-9/g_leak; + b = 0./g_leak; + + w = get_state_vector("w"); + + calculate_scale_constants(); + bg_current = get_state_vector("bg_current"); + + t_g_ampa = auryn_vector_float_ptr ( g_ampa , 0 ); + t_g_gaba = auryn_vector_float_ptr ( g_gaba , 0 ); + t_bg_cur = auryn_vector_float_ptr ( bg_current , 0 ); + t_mem = auryn_vector_float_ptr ( mem , 0 ); + t_w = auryn_vector_float_ptr ( w , 0 ); + + clear(); + +} + +void AdExGroup::clear() +{ + clear_spikes(); + for (NeuronID i = 0; i < get_rank_size(); i++) { + auryn_vector_float_set (mem, i, e_rest); + auryn_vector_float_set (g_ampa, i, 0.); + auryn_vector_float_set (g_gaba, i, 0.); + auryn_vector_float_set (bg_current, i, 0.); + } +} + + +AdExGroup::~AdExGroup() +{ +} + + +void AdExGroup::evolve() +{ + + // TODO we should vectorize this code and use some fast SSE + // library such as http://gruntthepeon.free.fr/ssemath/ + // for the exponential + for (NeuronID i = 0 ; i < get_rank_size() ; ++i ) { + t_w[i] += scale_w * (a * (t_mem[i]-e_rest) - t_w[i]); + + t_mem[i] += scale_mem * ( + e_rest-t_mem[i] + + deltat * exp((t_mem[i]-e_thr)/deltat) + - t_g_ampa[i] * (t_mem[i]-e_rev_ampa) + - t_g_gaba[i] * (t_mem[i]-e_rev_gaba) + + t_bg_cur[i]-t_w[i]); + + + if (t_mem[i]>0.0) { + push_spike(i); + t_mem[i] = e_reset; + t_w[i] += b; + } + } + + auryn_vector_float_scale(scale_ampa,g_ampa); + auryn_vector_float_scale(scale_gaba,g_gaba); +} + +void AdExGroup::set_bg_current(NeuronID i, AurynFloat current) +{ + if ( localrank(i) ) + auryn_vector_float_set ( bg_current , global2rank(i) , current ) ; +} + +void AdExGroup::set_tau_w(AurynFloat tauw) +{ + tau_w = tauw; + calculate_scale_constants(); +} + +void AdExGroup::set_b(AurynFloat _b) +{ + b = _b; +} + +void AdExGroup::set_e_reset(AurynFloat ereset) +{ + e_reset = ereset; +} + +void AdExGroup::set_e_thr(AurynFloat ethr) +{ + e_thr = ethr; +} +void AdExGroup::set_e_rest(AurynFloat erest) +{ + e_rest = erest; + for (NeuronID i = 0; i < get_rank_size(); i++) + auryn_vector_float_set (mem, i, e_rest); +} + +void AdExGroup::set_a(AurynFloat _a) +{ + a = _a; +} + +void AdExGroup::set_delta_t(AurynFloat d) +{ + deltat = d; +} + +void AdExGroup::set_g_leak(AurynFloat g) +{ + g_leak = g; + tau_mem = c_mem/g_leak; + calculate_scale_constants(); +} + +void AdExGroup::set_c_mem(AurynFloat cm) +{ + c_mem = cm; + tau_mem = c_mem/g_leak; + calculate_scale_constants(); +} + +AurynFloat AdExGroup::get_bg_current(NeuronID i) { + if ( localrank(i) ) + return auryn_vector_float_get ( bg_current , global2rank(i) ) ; + else + return 0; +} + +string AdExGroup::get_output_line(NeuronID i) +{ + stringstream oss; + oss << get_mem(i) << " " << get_ampa(i) << " " << get_gaba(i) << " " + << auryn_vector_float_get (bg_current, i) <<"\n"; + return oss.str(); +} + +void AdExGroup::load_input_line(NeuronID i, const char * buf) +{ + float vmem,vampa,vgaba,vbgcur; + sscanf (buf,"%f %f %f %f",&vmem,&vampa,&vgaba,&vbgcur); + if ( localrank(i) ) { + NeuronID trans = global2rank(i); + set_mem(trans,vmem); + set_ampa(trans,vampa); + set_gaba(trans,vgaba); + auryn_vector_float_set (bg_current, trans, vbgcur); + } +} + +void AdExGroup::set_tau_ampa(AurynFloat taum) +{ + tau_ampa = taum; + calculate_scale_constants(); +} + +AurynFloat AdExGroup::get_tau_ampa() +{ + return tau_ampa; +} + +void AdExGroup::set_tau_gaba(AurynFloat taum) +{ + tau_gaba = taum; + calculate_scale_constants(); +} + +AurynFloat AdExGroup::get_tau_gaba() +{ + return tau_gaba; +} diff --git a/src/AdExGroup.h b/src/AdExGroup.h new file mode 100644 index 00000000..946a352c --- /dev/null +++ b/src/AdExGroup.h @@ -0,0 +1,104 @@ +/* +* Copyright 2014-2015 Ankur Sinha and Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +* +* If you are using Auryn or parts of it for your work please cite: +* Zenke, F. and Gerstner, W., 2014. Limits to high-speed simulations +* of spiking neural networks using general-purpose computers. +* Front Neuroinform 8, 76. doi: 10.3389/fninf.2014.00076 +*/ + +#ifndef ADEXGROUP_H_ +#define ADEXGROUP_H_ + +#include "auryn_definitions.h" +#include "NeuronGroup.h" +#include "System.h" + + +/*! \brief Conductance based Adaptive Exponential neuron model - Brette and Gerstner (2005). Default values are taken from Table 1 (4a) of Naud, Marcille, Clopath and Gerstner (2008) +*/ +class AdExGroup : public NeuronGroup +{ +private: + auryn_vector_float * bg_current; + AurynFloat e_rest, e_reset, e_rev_gaba, e_rev_ampa,e_thr, g_leak, c_mem, deltat; + AurynFloat tau_ampa, tau_gaba, tau_mem; + AurynFloat scale_ampa, scale_gaba, scale_mem, scale_w; + AurynFloat * t_w; + AurynFloat a, tau_w, b; + + /*! Stores the adaptation current. */ + auryn_vector_float * w __attribute__((aligned(16))); + + AurynFloat * t_g_ampa; + AurynFloat * t_g_gaba; + AurynFloat * t_bg_cur; + AurynFloat * t_mem; + + void init(); + void calculate_scale_constants(); + inline void integrate_state(); + inline void check_thresholds(); + virtual string get_output_line(NeuronID i); + virtual void load_input_line(NeuronID i, const char * buf); + +public: + /*! The default constructor of this NeuronGroup */ + AdExGroup(NeuronID size); + virtual ~AdExGroup(); + + /*! Controls the constant current input to neuron i in natural units of g_leak (default 500pA/10ns) */ + void set_bg_current(NeuronID i, AurynFloat current); + + /*! Set value of slope factor deltat (default 2mV) */ + void set_delta_t(AurynFloat d); + /*! Set value of a in natural units of g_leak (default 2nS/10ns) */ + void set_a(AurynFloat _a); + /*! Set value of b in natural units of g_leak (default 0nS/10ns) */ + void set_b(AurynFloat _b); + /*! Set value of V_r (default -70mV) */ + void set_e_reset(AurynFloat ereset); + /*! Set value of E_l (default -70mV) */ + void set_e_rest(AurynFloat erest); + /*! Set value of V_t (default -50mV) */ + void set_e_thr(AurynFloat ethr); + /*! Sets the w time constant (default 30ms) */ + void set_tau_w(AurynFloat tauw); + /*! Gets the current background current value for neuron i */ + AurynFloat get_bg_current(NeuronID i); + /*! Sets the leak conductance (default 10nS) */ + void set_g_leak(AurynFloat g); + /*! Sets the membrane capacitance (default 200pF) */ + void set_c_mem(AurynFloat cm); + /*! Sets the exponential time constant for the AMPA channel (default 5ms) */ + void set_tau_ampa(AurynFloat tau); + /*! Gets the exponential time constant for the AMPA channel */ + AurynFloat get_tau_ampa(); + /*! Sets the exponential time constant for the GABA channel (default 10ms) */ + void set_tau_gaba(AurynFloat tau); + /*! Gets the exponential time constant for the GABA channel */ + AurynFloat get_tau_gaba(); + /*! Resets all neurons to defined and identical initial state. */ + void clear(); + /*! The evolve method internally used by System. */ + void evolve(); +}; + +#endif /*ADEXGROUP_H_*/ + diff --git a/src/Checker.h b/src/Checker.h index ed2c111e..3c5f47ae 100644 --- a/src/Checker.h +++ b/src/Checker.h @@ -45,6 +45,18 @@ class System; */ class Checker { +private: + /* Functions necesssary for serialization */ + friend class boost::serialization::access; + template + void serialize(Archive & ar, const unsigned int version) + { + virtual_serialize(ar, version); + } + + virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) = 0; + virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) = 0; + protected: SpikingGroup * src; @@ -63,4 +75,6 @@ class Checker virtual AurynFloat get_property() = 0 ; }; +BOOST_SERIALIZATION_ASSUME_ABSTRACT(Checker) + #endif /*CHECKER_H_*/ diff --git a/src/ComplexMatrix.cpp b/src/ComplexMatrix.cpp index 68df1c43..df4998fd 100644 --- a/src/ComplexMatrix.cpp +++ b/src/ComplexMatrix.cpp @@ -1 +1,2 @@ // This file was intentionally left empty and exists only for the sake of Makefile wildcards +#include "ComplexMatrix.h" diff --git a/src/ComplexMatrix.h b/src/ComplexMatrix.h index 7d268b54..0194e894 100644 --- a/src/ComplexMatrix.h +++ b/src/ComplexMatrix.h @@ -176,11 +176,18 @@ class ComplexMatrix bool exists(NeuronID i, NeuronID j); /*! Returns the pointer to a particular element */ T * get_ptr(NeuronID i, NeuronID j); + T * get_ptr(NeuronID i, NeuronID j, NeuronID z); /*! Returns the pointer to a particular element given * its position in the data array. */ T * get_ptr(AurynLong data_index); + /*! Returns data index to a particular element specifed by i and j */ + AurynLong get_data_index(NeuronID i, NeuronID j, NeuronID z=0); + + /*! Returns data index to a particular element specifed by a data pointer */ + AurynLong get_data_index(T * ptr); + /* Methods concerning synaptic state vectors. */ void set_num_synapse_states(StateID zsize); /*! Gets pointer for the first element of a given synaptic state vector */ @@ -255,22 +262,22 @@ class ComplexMatrix }; template -T * ComplexMatrix::get_data_ptr(AurynLong i, StateID z) +T * ComplexMatrix::get_data_ptr(const AurynLong i, const StateID z) { return elementdata+z*get_datasize()+i; } template -T ComplexMatrix::get_data(AurynLong i, StateID z) +T ComplexMatrix::get_data(const AurynLong i, const StateID z) { return elementdata[z*get_datasize()+i]; } template -T * ComplexMatrix::get_data_ptr(const NeuronID * ind_ptr, StateID z) +T * ComplexMatrix::get_data_ptr(const NeuronID * ind_ptr, const StateID z) { size_t ptr_offset = ind_ptr-get_ind_begin(); - return elementdata+ptr_offset; + return elementdata+z*get_datasize()+ptr_offset; } template @@ -531,6 +538,18 @@ T * ComplexMatrix::get_ptr(NeuronID i, NeuronID j) return NULL; } +template +AurynLong ComplexMatrix::get_data_index(NeuronID i, NeuronID j, NeuronID z) +{ + return get_data_index( get_ptr(i,j) ) + z*statesize ; +} + +template +AurynLong ComplexMatrix::get_data_index(T * ptr) +{ + return ptr - get_data_begin(); +} + template T * ComplexMatrix::get_ptr(AurynLong data_index) { diff --git a/src/Connection.cpp b/src/Connection.cpp index 0552768a..066c772a 100644 --- a/src/Connection.cpp +++ b/src/Connection.cpp @@ -89,6 +89,9 @@ void Connection::set_transmitter(TransmitterType transmitter) case MEM: set_transmitter(dst->get_mem_ptr()->data); break; + case CURSYN: + set_transmitter(dst->get_cursyn_ptr()->data); + break; case NMDA: set_transmitter(dst->get_nmda_ptr()->data); break; diff --git a/src/Connection.h b/src/Connection.h index 0c0d9505..f8f33156 100644 --- a/src/Connection.h +++ b/src/Connection.h @@ -60,11 +60,17 @@ class Connection string connection_name; protected: + /*! Serialization function for saving the Connection state. Implement in derived classes to save + * additional information. + * */ virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) { ar & m_rows & n_cols & connection_name; } + /*! Serialization function for loading the Connection state. Implement in derived classes to save + * additional information. + * */ virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) { ar & m_rows & n_cols & connection_name; diff --git a/src/Monitor.cpp b/src/Monitor.cpp index d9a12a09..f75f4b99 100644 --- a/src/Monitor.cpp +++ b/src/Monitor.cpp @@ -66,3 +66,11 @@ void Monitor::free() outfile.close(); } + +void Monitor::virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) +{ +} + +void Monitor::virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) +{ +} diff --git a/src/Monitor.h b/src/Monitor.h index c979a758..5aba573e 100644 --- a/src/Monitor.h +++ b/src/Monitor.h @@ -42,6 +42,16 @@ class System; class Monitor { +private: + /*! Functions necesssary for serialization and loading saving to netstate files. */ + friend class boost::serialization::access; + template + void serialize(Archive & ar, const unsigned int version) + { + virtual_serialize(ar, version); + } + + protected: /*! Output filestream to be used in the derived classes */ ofstream outfile; @@ -53,6 +63,10 @@ class Monitor virtual void open_output_file(string filename); /*! Standard free function to be called by the destructor - closes the file stream. */ void free(); + + /*! Functions necesssary for serialization and loading saving to netstate files. */ + virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) ; + virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) ; public: /*! Toggles Monitor (in)active */ @@ -67,6 +81,8 @@ class Monitor virtual void propagate() = 0; }; +BOOST_SERIALIZATION_ASSUME_ABSTRACT(Checker) + extern System * sys; extern Logger * logger; diff --git a/src/NeuronGroup.cpp b/src/NeuronGroup.cpp index 3f4f1952..b0dfc41c 100644 --- a/src/NeuronGroup.cpp +++ b/src/NeuronGroup.cpp @@ -41,6 +41,7 @@ void NeuronGroup::init() thr = get_state_vector("thr"); g_ampa = get_state_vector("g_ampa"); g_gaba = get_state_vector("g_gaba"); + g_cursyn = get_state_vector("g_cursyn"); g_nmda = get_state_vector("g_nmda"); #ifndef CODE_ALIGNED_SSE_INSTRUCTIONS @@ -49,7 +50,8 @@ void NeuronGroup::init() || auryn_AlignOffset( thr->size, thr->data, sizeof(float), 16) || auryn_AlignOffset( g_ampa->size, g_ampa->data, sizeof(float), 16) || auryn_AlignOffset( g_nmda->size, g_nmda->data, sizeof(float), 16) - || auryn_AlignOffset( g_gaba->size, g_gaba->data, sizeof(float), 16) ) + || auryn_AlignOffset( g_gaba->size, g_gaba->data, sizeof(float), 16) + || auryn_AlignOffset( g_cursyn->size, g_cursyn->data, sizeof(float), 16) ) throw AurynMemoryAlignmentException(); #endif } @@ -80,7 +82,7 @@ void NeuronGroup::print_vec(auryn_vector_float * vec, const char * name) void NeuronGroup::print_state(NeuronID id) { - printf ("%g %g %g %g\n",auryn_vector_float_get (mem, id), auryn_vector_float_get (g_ampa, id), auryn_vector_float_get (g_gaba, id) ,auryn_vector_float_get (g_nmda, id)); + printf ("%g %g %g %g %g\n",auryn_vector_float_get (mem, id), auryn_vector_float_get (g_ampa, id), auryn_vector_float_get (g_gaba, id) ,auryn_vector_float_get (g_nmda, id), auryn_vector_float_get (g_cursyn, id)); } AurynState NeuronGroup::get_mem(NeuronID i) @@ -121,11 +123,21 @@ AurynState NeuronGroup::get_ampa(NeuronID i) return get_val(g_ampa,i); } +AurynState NeuronGroup::get_cursyn(NeuronID i) +{ + return get_val(g_cursyn,i); +} + auryn_vector_float * NeuronGroup::get_ampa_ptr() { return g_ampa; } +auryn_vector_float * NeuronGroup::get_cursyn_ptr() +{ + return g_cursyn; +} + AurynState NeuronGroup::get_gaba(NeuronID i) { return get_val(g_gaba,i); @@ -151,6 +163,11 @@ void NeuronGroup::set_ampa(NeuronID i, AurynState val) set_val(g_ampa,i,val); } +void NeuronGroup::set_cursyn(NeuronID i, AurynState val) +{ + set_val(g_cursyn,i,val); +} + void NeuronGroup::set_gaba(NeuronID i, AurynState val) { set_val(g_gaba,i,val); @@ -216,6 +233,11 @@ void NeuronGroup::print_nmda() print_vec(g_nmda,"g_nmda"); } +void NeuronGroup::print_cursyn() +{ + print_vec(g_nmda,"g_cursyn"); +} + void NeuronGroup::safe_tadd(NeuronID id, AurynWeight amount, TransmitterType t) { if (localrank(id)) diff --git a/src/NeuronGroup.h b/src/NeuronGroup.h index 844370bd..d9adaf99 100644 --- a/src/NeuronGroup.h +++ b/src/NeuronGroup.h @@ -55,6 +55,8 @@ class NeuronGroup : public SpikingGroup auryn_vector_float * g_gaba __attribute__((aligned(16))); /*! Stores the NMDA conductances of each point neuron. */ auryn_vector_float * g_nmda __attribute__((aligned(16))); + /*! Stores the CURSYN states of each point neuron. */ + auryn_vector_float * g_cursyn __attribute__((aligned(16))); /*! Stores threshold terms for moving thresholds. */ auryn_vector_float * thr __attribute__((aligned(16))); @@ -91,9 +93,15 @@ class NeuronGroup : public SpikingGroup AurynState get_ampa(NeuronID i); void set_ampa(NeuronID i,AurynState val); auryn_vector_float * get_ampa_ptr(); + + AurynState get_cursyn(NeuronID i); + void set_cursyn(NeuronID i,AurynState val); + auryn_vector_float * get_cursyn_ptr(); + AurynState get_gaba(NeuronID i); void set_gaba(NeuronID i,AurynState val); auryn_vector_float * get_gaba_ptr(); + AurynState get_nmda(NeuronID i); void set_nmda(NeuronID i,AurynState val); auryn_vector_float * get_nmda_ptr(); @@ -109,6 +117,7 @@ class NeuronGroup : public SpikingGroup void print_ampa(); void print_gaba(); void print_nmda(); + void print_cursyn(); void print_state(NeuronID id); void safe_tadd(NeuronID id, AurynWeight amount, TransmitterType t=GLUT); /*! Adds given transmitter to neuron as from a synaptic event. DEPRECATED. Moving slowly to SparseConnection transmit. */ diff --git a/src/NormalStimulator.cpp b/src/NormalStimulator.cpp new file mode 100644 index 00000000..83fccd8b --- /dev/null +++ b/src/NormalStimulator.cpp @@ -0,0 +1,94 @@ +/* +* Copyright 2014 Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#include "NormalStimulator.h" + +boost::mt19937 NormalStimulator::gen = boost::mt19937(); + +NormalStimulator::NormalStimulator(NeuronGroup * target, AurynWeight sigma, string target_state ) : Monitor( ) +{ + init(target, sigma, target_state); +} + + +void NormalStimulator::init( NeuronGroup * target, AurynWeight sigma, string target_state ) +{ + sys->register_monitor(this); + dst = target; + + set_target_state(target_state); + + normal_sigma = sigma; + + + stringstream oss; + oss << scientific << "NormalStimulator:: initializing with mean " << get_lambda(); + logger->msg(oss.str(),NOTIFICATION); + + seed(61093*communicator->rank()); + dist = new boost::normal_distribution (0.0, get_lambda()); + die = new boost::variate_generator > ( gen, *dist ); +} + +void NormalStimulator::free( ) +{ + delete dist; + delete die; +} + + +NormalStimulator::~NormalStimulator() +{ + free(); +} + +void NormalStimulator::propagate() +{ + if ( dst->evolve_locally() ) { + for ( NeuronID i = 0 ; i < dst->get_post_size() ; ++i ) { + float draw = (*die)(); + target_vector->data[i] = draw; + } + } +} + +void NormalStimulator::set_sigma(AurynWeight sigma) { + delete dist; + normal_sigma = sigma; + dist = new boost::normal_distribution (0.0, get_lambda()); +} + +AurynFloat NormalStimulator::get_sigma() { + return (AurynFloat) normal_sigma; +} + +AurynFloat NormalStimulator::get_lambda() { + return get_sigma(); +} + +void NormalStimulator::set_target_state(string state_name) { + target_vector = dst->get_state_vector(state_name); +} + +void NormalStimulator::seed(int s) +{ + gen.seed(s); +} + diff --git a/src/NormalStimulator.h b/src/NormalStimulator.h new file mode 100644 index 00000000..d4f864c5 --- /dev/null +++ b/src/NormalStimulator.h @@ -0,0 +1,107 @@ +/* +* Copyright 2014 Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#ifndef NORMALSTIMULATOR_H_ +#define NORMALSTIMULATOR_H_ + +#include "auryn_definitions.h" +#include "System.h" +#include "Logger.h" +#include "Monitor.h" +#include "NeuronGroup.h" + +#include +#include +#include + +using namespace std; + +/*! \brief Stimulator class to inject timeseries of currents to patterns (subpopulations) of neurons + * + * Instances of this class inject currents that vary over time to subpopulations of the NeuronGroup assigned. + */ + +class NormalStimulator : protected Monitor +{ +private: + + static boost::mt19937 gen; + boost::normal_distribution * dist; + boost::variate_generator > * die; + + + /*! Vector storing all the current values */ + AurynState * currents; + + /*! Vector storing all new current values */ + AurynState * newcurrents; + + /*! Target membrane */ + auryn_vector_float * target_vector; + + /*! Scale stimulus size */ + AurynWeight normal_sigma; + + /*! Default init method */ + void init(NeuronGroup * target, AurynWeight sigma, string target_state); + + void free(); + + /*! Returns the lambda parameter of the pmf for Normal. */ + AurynFloat get_lambda(); + +protected: + + /*! The target NeuronGroup */ + NeuronGroup * dst; + + +public: + /*! Default Constructor + * @param[target] The target spiking group. + * @param[rate] The firing rate of each the Normal process. + * @param[weight] The weight or unit of amount of change on the state variable + */ + NormalStimulator(NeuronGroup * target, AurynWeight sigma=1.0, string target_state="inj_current"); + + /*! Default Destructor */ + virtual ~NormalStimulator(); + + + /*! Sets the event rate of the underlying Normal generator */ + void set_sigma(AurynFloat sigma); + + /*! Returns the event rate of the underlying Normal generator */ + AurynFloat get_sigma(); + + /*! Seeds the random number generator of all NormalStimulator objects */ + void seed(int s); + + + /*! Sets the state that is stimulated with Normal input. + * This must be a valid state vector name (default = mem) */ + void set_target_state( string state_name = "inj_current" ); + + /*! Implementation of necessary propagate() function. */ + void propagate(); + +}; + +#endif /*POISSONSTIMULATOR_H_*/ diff --git a/src/PairInteractionConnection.cpp b/src/PairInteractionConnection.cpp new file mode 100644 index 00000000..bf69539c --- /dev/null +++ b/src/PairInteractionConnection.cpp @@ -0,0 +1,243 @@ +/* +* Copyright 2014-2015 Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#include "PairInteractionConnection.h" + +void PairInteractionConnection::init(AurynWeight maxw) +{ + last_spike_pre = new AurynTime[get_m_rows()]; + last_spike_post = new AurynTime[get_n_cols()]; + + for ( unsigned int i = 0 ; i < get_m_rows() ; ++i ) + last_spike_pre[i] = -1; // set this to end of range (so spike is infinitely in the future -- might cause problems without ffast-math + + for ( unsigned int i = 0 ; i < get_n_cols() ; ++i ) + last_spike_post[i] = -1; + + window_pre_post = new AurynFloat[WINDOW_MAX_SIZE]; + window_post_pre = new AurynFloat[WINDOW_MAX_SIZE]; + + // initialize window with standard exponential 20ms time constant + + set_exponential_window(); + + + // TODO write proper init for these variables + stdp_active = true; + w_max = maxw; + + set_name("PairInteractionConnection"); +} + +void PairInteractionConnection::free() +{ + delete last_spike_pre; + delete last_spike_post; + + delete window_pre_post; + delete window_post_pre; +} + +PairInteractionConnection::PairInteractionConnection(SpikingGroup * source, NeuronGroup * destination, + const char * filename, + AurynWeight maxweight , TransmitterType transmitter) +: DuplexConnection(source, destination, filename, transmitter) +{ + init(maxweight); +} + +PairInteractionConnection::PairInteractionConnection(SpikingGroup * source, NeuronGroup * destination, + AurynWeight weight, AurynFloat sparseness, + AurynWeight maxweight , TransmitterType transmitter, string name) +: DuplexConnection(source, destination, weight, sparseness, transmitter, name) +{ + init(maxweight); +} + +PairInteractionConnection::~PairInteractionConnection() +{ + free(); +} + +inline AurynWeight PairInteractionConnection::dw_fwd(NeuronID post) +{ + AurynTime diff = sys->get_clock()-last_spike_post[post]; + if ( stdp_active ) { + if ( diff >= WINDOW_MAX_SIZE ) diff = WINDOW_MAX_SIZE-1; + double dw = window_post_pre[diff]; + return dw; + } + else return 0.; +} + +inline AurynWeight PairInteractionConnection::dw_bkw(NeuronID pre) +{ + AurynTime diff = sys->get_clock()-last_spike_pre[pre]; + if ( stdp_active ) { + if ( diff >= WINDOW_MAX_SIZE ) diff = WINDOW_MAX_SIZE-1; + double dw = window_pre_post[diff]; + return dw; + } + else return 0.; +} + +inline void PairInteractionConnection::propagate_forward() +{ + NeuronID * ind = w->get_row_begin(0); // first element of index array + AurynWeight * data = w->get_data_begin(); + AurynWeight value; + TransmitterType transmitter = get_transmitter(); + SpikeContainer::const_iterator spikes_end = src->get_spikes()->end(); + // process spikes + for (SpikeContainer::const_iterator spike = src->get_spikes()->begin() ; // spike = pre_spike + spike != spikes_end ; ++spike ) { + for (NeuronID * c = w->get_row_begin(*spike) ; c != w->get_row_end(*spike) ; ++c ) { + value = data[c-ind]; + //dst->tadd( *c , value , transmitter ); + transmit( *c, value ); + //if (data[c-ind]>0 && data[c-ind]get_clock(); + } +} + +inline void PairInteractionConnection::propagate_backward() +{ + NeuronID * ind = bkw->get_row_begin(0); // first element of index array + AurynWeight ** data = bkw->get_data_begin(); + SpikeContainer::const_iterator spikes_end = dst->get_spikes_immediate()->end(); + for (SpikeContainer::const_iterator spike = dst->get_spikes_immediate()->begin() ; // spike = post_spike + spike != spikes_end ; ++spike ) { + for (NeuronID * c = bkw->get_row_begin(*spike) ; c != bkw->get_row_end(*spike) ; ++c ) { + if (*data[c-ind]get_clock(); + } +} + +void PairInteractionConnection::propagate() +{ + // propagate + propagate_forward(); + propagate_backward(); +} + +void PairInteractionConnection::load_window_from_file( const char * filename , double scale ) +{ + + stringstream oss; + oss << "PairInteractionConnection:: Loading STDP window from " << filename; + logger->msg(oss.str(),NOTIFICATION); + + // default window all zeros + for ( int i = 0 ; i < WINDOW_MAX_SIZE ; ++i ) { + window_pre_post[i] = 0; + window_post_pre[i] = 0; + } + + ifstream infile (filename); + if (!infile) { + stringstream oes; + oes << "Can't open input file " << filename; + logger->msg(oes.str(),ERROR); + return; + } + + unsigned int size; + float timebinsize; + float value; + float time; + unsigned int count = 0; + + char buffer[256]; + infile.getline (buffer,256); + sscanf (buffer,"# %u %f",&size,&timebinsize); + + if ( size > 2*WINDOW_MAX_SIZE ) + logger->msg("PairInteractionConnection:: STDP window too large ... truncating!",WARNING); + + if ( dt < timebinsize ) + logger->msg("PairInteractionConnection:: Timebinning of loaded STDP window is different from simulator timestep.",WARNING); + + double sum_pre_post = 0 ; + double sum_post_pre = 0 ; + + // read window file line-by-line + while ( infile.getline (buffer,256) ) + { + sscanf (buffer,"%f %f",&time,&value); + if ( abs(time) < WINDOW_MAX_SIZE*dt ) { + NeuronID start; + if ( time < 0 ) { + start = -(time+dt/2)/dt; // plus element is for correct rounding + window_post_pre[start] = scale*value; + sum_post_pre += scale*value; + } else { + start = (time+dt/2)/dt; + window_pre_post[start] = scale*value; + sum_pre_post += scale*value; + } + } + count++; + } + + // for ( int i = 0 ; i < WINDOW_MAX_SIZE ; ++i ) { + // cout << scientific << window_pre_post[i] << endl; + // } + // for ( int i = 0 ; i < WINDOW_MAX_SIZE ; ++i ) { + // cout << scientific << window_post_pre[i] << endl; + // } + + + oss.str(""); + oss << "PairInteractionConnection:: sum_pre_post=" + << scientific + << sum_pre_post + << " sum_post_pre=" + << sum_post_pre; + logger->msg(oss.str(),NOTIFICATION); + + infile.close(); + +} + +void PairInteractionConnection::set_exponential_window ( double Aplus, double tau_plus, double Aminus, double tau_minus) +{ + for ( int i = 0 ; i < WINDOW_MAX_SIZE ; ++i ) { + window_pre_post[i] = Aplus/tau_plus*exp(-i*dt/tau_plus); + } + + for ( int i = 0 ; i < WINDOW_MAX_SIZE ; ++i ) { + window_post_pre[i] = Aminus/tau_minus*exp(-i*dt/tau_minus); + } + + // zero floor terms + set_floor_terms(); +} + +void PairInteractionConnection::set_floor_terms( double pre_post, double post_pre ) +{ + window_pre_post[WINDOW_MAX_SIZE-1] = pre_post; + window_post_pre[WINDOW_MAX_SIZE-1] = post_pre; +} diff --git a/src/PairInteractionConnection.h b/src/PairInteractionConnection.h new file mode 100644 index 00000000..4445dc75 --- /dev/null +++ b/src/PairInteractionConnection.h @@ -0,0 +1,74 @@ +/* +* Copyright 2014-2015 Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#ifndef PAIRINTERACTIONCONNECTION_H_ +#define PAIRINTERACTIONCONNECTION_H_ + +#define WINDOW_MAX_SIZE 60000 + +#include "auryn_definitions.h" +#include "DuplexConnection.h" +#include "EulerTrace.h" + +using namespace std; + + +class PairInteractionConnection : public DuplexConnection +{ + +protected: + AurynWeight w_max; + + AurynTime * last_spike_pre; + AurynTime * last_spike_post; + + inline AurynWeight dw_fwd(NeuronID post); + inline AurynWeight dw_bkw(NeuronID pre); + + inline void propagate_forward(); + inline void propagate_backward(); + + +public: + AurynFloat * window_pre_post; + AurynFloat * window_post_pre; + + bool stdp_active; + + PairInteractionConnection(SpikingGroup * source, NeuronGroup * destination, + const char * filename, + AurynWeight maxweight=1. , TransmitterType transmitter=GLUT); + + PairInteractionConnection(SpikingGroup * source, NeuronGroup * destination, + AurynWeight weight, AurynFloat sparseness=0.05, + AurynWeight maxweight=1. , TransmitterType transmitter=GLUT, string name="PairInteractionConnection"); + virtual ~PairInteractionConnection(); + void init(AurynWeight maxw); + void free(); + + void load_window_from_file( const char * filename , double scale = 1. ); + void set_exponential_window ( double Aplus = 1e-3, double tau_plus = 20e-3, double Aminus = -1e-3, double tau_minus = 20e-3); + void set_floor_terms( double pre_post = 0.0, double post_pre = 0.0 ); + + virtual void propagate(); + +}; + +#endif /*PAIRINTERACTIONCONNECTION_H_*/ diff --git a/src/PatternMonitor.cpp b/src/PatternMonitor.cpp index a7b347c5..48e1bb26 100644 --- a/src/PatternMonitor.cpp +++ b/src/PatternMonitor.cpp @@ -180,3 +180,15 @@ void PatternMonitor::load_patterns( string filename ) oss << "PatternMonitor:: Finished loading " << patterns->size() << " patterns"; logger->msg(oss.str(),NOTIFICATION); } + +void PatternMonitor::virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) +{ + for ( NeuronID i = 0 ; i < src->get_rank_size() ; ++i ) + ar & counter[i] ; +} + +void PatternMonitor::virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) +{ + for ( NeuronID i = 0 ; i < src->get_rank_size() ; ++i ) + ar & counter[i] ; +} diff --git a/src/PatternMonitor.h b/src/PatternMonitor.h index c534fa91..b1cc72b7 100644 --- a/src/PatternMonitor.h +++ b/src/PatternMonitor.h @@ -64,6 +64,9 @@ class PatternMonitor : protected Monitor SpikingGroup * src; /*! Default init method */ void init(SpikingGroup * source, string filename, NeuronID maximum_patterns, AurynFloat binsize); + + virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ); + virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ); public: /*! Default Constructor diff --git a/src/PatternStimulator.cpp b/src/PatternStimulator.cpp index 13d6ae1a..1f88501b 100644 --- a/src/PatternStimulator.cpp +++ b/src/PatternStimulator.cpp @@ -74,11 +74,16 @@ void PatternStimulator::propagate() { if ( dst->evolve_locally() ) { - char buffer[256]; + char buffer[256]; string line; while( !timeseriesfile.eof() && filetime < sys->get_clock() ) { line.clear(); - timeseriesfile.getline (buffer,255); + timeseriesfile.getline (buffer,255); + // FIXME This buffer can quickly become to small of read lots of columns + // Instead this section should be re-written to directly read one token/column + // at a time in the loop below. Just making the buffer size large should be avoided + // not to risk buffer overflow or unecessarily large memory consumption on cluster + // machines with limited resources. line = buffer; if (line[0] == '#') continue; stringstream iss (line); diff --git a/src/PopulationRateMonitor.cpp b/src/PopulationRateMonitor.cpp index 41cec59b..8846eb14 100644 --- a/src/PopulationRateMonitor.cpp +++ b/src/PopulationRateMonitor.cpp @@ -62,3 +62,13 @@ void PopulationRateMonitor::propagate() } } } + +void PopulationRateMonitor::virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) +{ + ar & counter ; +} + +void PopulationRateMonitor::virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) +{ + ar & counter ; +} diff --git a/src/PopulationRateMonitor.h b/src/PopulationRateMonitor.h index 20f7b38c..e10546c4 100644 --- a/src/PopulationRateMonitor.h +++ b/src/PopulationRateMonitor.h @@ -58,6 +58,8 @@ class PopulationRateMonitor : protected Monitor /*! Default init method */ void init(SpikingGroup * source, string filename, AurynDouble binsize); + virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ); + virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ); public: /*! Default Constructor @param[source] The source spiking group. diff --git a/src/ProfilePoissonGroup.h b/src/ProfilePoissonGroup.h index ef23409b..20786a63 100644 --- a/src/ProfilePoissonGroup.h +++ b/src/ProfilePoissonGroup.h @@ -40,7 +40,17 @@ using namespace std; /*! \brief A SpikingGroup that creates poissonian spikes with a given rate - * and spacial profile. + * and spatial profile. + * + * This SpikingGroup is a logic extension of PoissonGroup for (relatively) stationary, + * but not uniform firing rates. It uses a similar algorithm as PoissonGroup to generate spikes + * in which each random number yields a spike, but uses a warped output array in which neurons + * can have different firing probabilities in each timestep. The resulting implementation requires + * the computation of the cumulative firing probability accross the group in every timestep. It is + * therefore substantially slower than PoissonGroup, but presumably much faster than drawing + * random numbers for each neuron in each time step. + * To set the firing rate profile use the function set_profile which needs to point to an array at + * which the firing rate profile is stored. */ class ProfilePoissonGroup : public SpikingGroup { diff --git a/src/RateChecker.cpp b/src/RateChecker.cpp index d2fd4567..428b92e8 100644 --- a/src/RateChecker.cpp +++ b/src/RateChecker.cpp @@ -70,7 +70,22 @@ AurynFloat RateChecker::get_rate() return state; } +void RateChecker::set_rate(AurynFloat r) +{ + state = r; +} + void RateChecker::reset() { - state = (popmax+popmin)/2; + set_rate((popmax+popmin)/2); +} + +void RateChecker::virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) +{ + ar & state; +} + +void RateChecker::virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) +{ + ar & state; } diff --git a/src/RateChecker.h b/src/RateChecker.h index e4da991d..2a480d8a 100644 --- a/src/RateChecker.h +++ b/src/RateChecker.h @@ -33,22 +33,26 @@ using namespace std; -/*! \brief A Checker class that tracks population firing rate and breaks - * a run if it goes out of bound. +/*! \brief A Checker class that tracks population firing rate as a moving + * average and breaks a run if it goes out of bound. * - * This class should be specified at least once with all plastic - * runs to ensure that expoding firing rates or a silent network - * does not get simulated needlessly for hours. + * This class should be specified at least once with all plastic runs to ensure + * that expoding firing rates or a silent network does not get simulated + * needlessly for hours. * - * The different constructors allow to specify different min and max - * firing rates to guard against too active or quiet networks. - * Also the timeconstant (tau) over which the moving rate average - * is computed online can be specified. Note that for highly parallel - * simulations the number of neurons in a specific group on a particular - * node can be highly reduced and therefore the rate estimate - * might become noisy. If highly parallel simulation is anticipated - * tau should be chosen longer to avoid spurious breaks caused by - * a noisy rate estimate. + * The different constructors allow to specify different min and max firing + * rates to guard against too active or quiet networks. Also the timeconstant + * (tau) over which the moving rate average is computed online can be specified. + * Allow for 3-5 x tau for the estimate to settle to its steady state value. To + * avoid accidental breaking of a run due to this effect, at initialization the + * rate estimate is assumed to be the mean of the min and max. Note further + * that this checker computes population averages over the fraction of a neuron + * group which is simulated on a particular rank. In highly parallel + * simulations when the number of neurons per rank is very the rate estimate + * might have a high variance accross ranks. If highly parallel simulation is + * anticipated tau should be chosen longer to avoid spurious breaks caused by a + * noisy rate estimate or a different checker which computes the rate of entire + * population (after a MINDELAY s minimal delay) should be used. */ class RateChecker : public Checker @@ -62,6 +66,8 @@ class RateChecker : public Checker NeuronID size; void init(AurynFloat min, AurynFloat max, AurynFloat tau); + virtual void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ); + virtual void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ); public: /*! The default constructor. * @param source the source group to monitor. @@ -82,6 +88,9 @@ class RateChecker : public Checker virtual AurynFloat get_property(); /*! Reads out the current rate estimate. */ AurynFloat get_rate(); + /*! Sets the current rate estimate -- for instance to provide a reasonable guess upon init. + * ( per default this is (max+min)/2.) */ + void set_rate(AurynFloat r); void reset(); }; diff --git a/src/SIFGroup.cpp b/src/SIFGroup.cpp new file mode 100644 index 00000000..47c06d0c --- /dev/null +++ b/src/SIFGroup.cpp @@ -0,0 +1,199 @@ +/* +* Copyright 2015 Neftci Emre and Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#include "SIFGroup.h" + +SIFGroup::SIFGroup(NeuronID size) : NeuronGroup(size) +{ + sys->register_spiking_group(this); + if ( evolve_locally() ) init(); +} + +void SIFGroup::calculate_scale_constants() +{ + scale_mem = dt/tau_mem; + scale_cursyn = exp(-dt/tau_cursyn); + scale_ampa = exp(-dt/tau_ampa); + scale_gaba = exp(-dt/tau_gaba); +} + +void SIFGroup::init() +{ + e_rest = 0e-3; + e_rev = -80e-3; + thr = .1; + tau_ampa = 4e-3; + tau_gaba = 4e-3; + tau_cursyn = 4e-3; + tau_mem = 1e-3; + set_refractory_period(4e-3); + + calculate_scale_constants(); + + ref = auryn_vector_ushort_alloc (get_vector_size()); + bg_current = get_state_vector("bg_current"); + inj_current = get_state_vector("inj_current"); + + t_g_ampa = auryn_vector_float_ptr ( g_ampa , 0 ); + t_g_gaba = auryn_vector_float_ptr ( g_gaba , 0 ); + t_g_cursyn = auryn_vector_float_ptr ( g_cursyn , 0 ); + t_bg_cur = auryn_vector_float_ptr ( bg_current , 0 ); + t_inj_cur = auryn_vector_float_ptr ( inj_current , 0 ); + t_mem = auryn_vector_float_ptr ( mem , 0 ); + t_ref = auryn_vector_ushort_ptr ( ref , 0 ); + + clear(); + +} + +void SIFGroup::clear() +{ + clear_spikes(); + for (NeuronID i = 0; i < get_rank_size(); i++) { + auryn_vector_ushort_set (ref, i, 0); + auryn_vector_float_set (mem, i, e_rest); + auryn_vector_float_set (g_ampa, i, 0.); + auryn_vector_float_set (g_gaba, i, 0.); + auryn_vector_float_set (g_cursyn, i, 0.); + auryn_vector_float_set (bg_current, i, 0.); + auryn_vector_float_set (inj_current, i, 0.); + } +} + + +SIFGroup::~SIFGroup() +{ + if ( !evolve_locally() ) return; + + auryn_vector_ushort_free (ref); +} + + +void SIFGroup::evolve() +{ + + + for (NeuronID i = 0 ; i < get_rank_size() ; ++i ) { + if (t_ref[i]==0) { + const AurynFloat dg_mem = ( + + (e_rest-t_mem[i]) + - t_g_ampa[i] * t_mem[i] + - t_g_gaba[i] * (t_mem[i]-e_rev) + + t_g_cursyn[i] + + t_inj_cur[i] + + t_bg_cur[i] ); + t_mem[i] += dg_mem*scale_mem; + + if (t_mem[i]>thr) { + push_spike(i); + t_mem[i] = e_rest ; + t_ref[i] += refractory_time ; + } + } else { + t_ref[i]-- ; + t_mem[i] = e_rest ; + } + + } + + auryn_vector_float_scale(scale_ampa,g_ampa); + auryn_vector_float_scale(scale_gaba,g_gaba); + auryn_vector_float_scale(scale_cursyn,g_cursyn); +} + +void SIFGroup::set_bg_current(NeuronID i, AurynFloat current) { + if ( localrank(i) ) + auryn_vector_float_set ( bg_current , global2rank(i) , current ) ; +} + +void SIFGroup::set_tau_mem(AurynFloat taum) +{ + tau_mem = taum; + calculate_scale_constants(); +} + +AurynFloat SIFGroup::get_bg_current(NeuronID i) { + if ( localrank(i) ) + return auryn_vector_float_get ( bg_current , global2rank(i) ) ; + else + return 0; +} + +string SIFGroup::get_output_line(NeuronID i) +{ + stringstream oss; + oss << get_mem(i) << " " << get_ampa(i) << " " << get_gaba(i) << " " << auryn_vector_ushort_get (ref, i) << "\n"; + return oss.str(); +} + +void SIFGroup::load_input_line(NeuronID i, const char * buf) +{ + float vmem,vampa,vgaba; + NeuronID vref; + sscanf (buf,"%f %f %f %u",&vmem,&vampa,&vgaba,&vref); + if ( localrank(i) ) { + NeuronID trans = global2rank(i); + set_mem(trans,vmem); + set_ampa(trans,vampa); + set_gaba(trans,vgaba); + auryn_vector_ushort_set (ref, trans, vref); + } +} + +void SIFGroup::set_tau_ampa(AurynFloat taum) +{ + tau_ampa = taum; + calculate_scale_constants(); +} + +AurynFloat SIFGroup::get_tau_ampa() +{ + return tau_ampa; +} + +void SIFGroup::set_tau_gaba(AurynFloat taum) +{ + tau_gaba = taum; + calculate_scale_constants(); +} + +AurynFloat SIFGroup::get_tau_gaba() +{ + return tau_gaba; +} + +void SIFGroup::set_refractory_period(AurynDouble t) +{ + double tmp = (unsigned short) (t/dt) - 1; + if (tmp<0) tmp = 0; + refractory_time = tmp; +} + +void SIFGroup::virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) +{ + SpikingGroup::virtual_serialize(ar,version); + ar & *ref; +} + +void SIFGroup::virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) +{ + SpikingGroup::virtual_serialize(ar,version); + ar & *ref; +} diff --git a/src/SIFGroup.h b/src/SIFGroup.h new file mode 100644 index 00000000..9eb6a72a --- /dev/null +++ b/src/SIFGroup.h @@ -0,0 +1,89 @@ +/* +* Copyright 2015 Neftci Emre and Friedemann Zenke +* +* This file is part of Auryn, a simulation package for plastic +* spiking neural networks. +* +* Auryn is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* Auryn is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with Auryn. If not, see . +*/ + +#ifndef SIFGROUP_H_ +#define SIFGROUP_H_ + +#include "auryn_definitions.h" +#include "NeuronGroup.h" +#include "System.h" + + +/*! \brief Conductance based neuron model with absolute refractoriness used for Neural Sampling + */ +class SIFGroup : public NeuronGroup +{ +private: + auryn_vector_float * bg_current; + auryn_vector_float * inj_current; + auryn_vector_ushort * ref; + unsigned short refractory_time; + AurynFloat e_rest,e_rev,thr,tau_mem; + AurynFloat tau_ampa,tau_gaba,tau_cursyn; + AurynFloat scale_ampa, scale_gaba, scale_mem, scale_cursyn; + + AurynFloat * t_g_cursyn; + AurynFloat * t_g_ampa; + AurynFloat * t_g_gaba; + AurynFloat * t_bg_cur; + AurynFloat * t_inj_cur; + AurynFloat * t_mem; + unsigned short * t_ref; + + void init(); + void calculate_scale_constants(); + inline void integrate_state(); + inline void check_thresholds(); + virtual string get_output_line(NeuronID i); + virtual void load_input_line(NeuronID i, const char * buf); + + void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ); + void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ); +public: + /*! The default constructor of this NeuronGroup */ + SIFGroup(NeuronID size); + virtual ~SIFGroup(); + + /*! Controls the constant current input (per default set so zero) to neuron i */ + void set_bg_current(NeuronID i, AurynFloat current); + + /*! Setter for refractory time [s] */ + void set_refractory_period(AurynDouble t); + + /*! Gets the current background current value for neuron i */ + AurynFloat get_bg_current(NeuronID i); + /*! Sets the membrane time constant (default 20ms) */ + void set_tau_mem(AurynFloat taum); + /*! Sets the exponential time constant for the AMPA channel (default 5ms) */ + void set_tau_ampa(AurynFloat tau); + /*! Gets the exponential time constant for the AMPA channel */ + AurynFloat get_tau_ampa(); + /*! Sets the exponential time constant for the GABA channel (default 10ms) */ + void set_tau_gaba(AurynFloat tau); + /*! Gets the exponential time constant for the GABA channel */ + AurynFloat get_tau_gaba(); + /*! Resets all neurons to defined and identical initial state. */ + void clear(); + /*! The evolve method internally used by System. */ + void evolve(); +}; + +#endif /*SIFGROUP_H_*/ + diff --git a/src/STDPConnection.cpp b/src/STDPConnection.cpp index 28a9be5f..110824f8 100644 --- a/src/STDPConnection.cpp +++ b/src/STDPConnection.cpp @@ -20,13 +20,10 @@ #include "STDPConnection.h" -void STDPConnection::init(AurynFloat eta, AurynFloat maxweight) +void STDPConnection::init(AurynFloat eta, AurynFloat tau_pre, AurynFloat tau_post, AurynFloat maxweight) { if ( dst->get_post_size() == 0 ) return; - tau_pre = 20.0e-3; - tau_post = 20.0e-3; - A = eta; // post-pre B = eta; // pre-post @@ -58,21 +55,25 @@ STDPConnection::STDPConnection(SpikingGroup * source, NeuronGroup * destination, STDPConnection::STDPConnection(SpikingGroup * source, NeuronGroup * destination, const char * filename, - AurynFloat eta, - AurynFloat maxweight , + AurynFloat eta, + AurynFloat tau_pre, + AurynFloat tau_post, + AurynFloat maxweight, TransmitterType transmitter) : DuplexConnection(source, destination, filename, transmitter) { - init(eta, maxweight); + init(eta, tau_pre, tau_post, maxweight); } STDPConnection::STDPConnection(SpikingGroup * source, NeuronGroup * destination, AurynWeight weight, AurynFloat sparseness, AurynFloat eta, - AurynFloat maxweight , + AurynFloat tau_pre, + AurynFloat tau_post, + AurynFloat maxweight, TransmitterType transmitter, string name) : DuplexConnection(source, @@ -82,7 +83,7 @@ STDPConnection::STDPConnection(SpikingGroup * source, NeuronGroup * destination, transmitter, name) { - init(eta, maxweight); + init(eta, tau_pre, tau_post, maxweight); if ( name.empty() ) set_name("STDPConnection"); } diff --git a/src/STDPConnection.h b/src/STDPConnection.h index 126236a7..3fb1dde4 100644 --- a/src/STDPConnection.h +++ b/src/STDPConnection.h @@ -40,13 +40,10 @@ class STDPConnection : public DuplexConnection { private: - void init(AurynFloat eta, AurynFloat maxweight); + void init(AurynFloat eta, AurynFloat tau_pre, AurynFloat tau_post, AurynFloat maxweight); protected: - AurynFloat tau_pre; - AurynFloat tau_post; - AurynDouble hom_fudge; PRE_TRACE_MODEL * tr_pre; @@ -70,16 +67,21 @@ class STDPConnection : public DuplexConnection STDPConnection(SpikingGroup * source, NeuronGroup * destination, const char * filename, AurynFloat eta=1, + AurynFloat tau_pre=20e-3, + AurynFloat tau_post=20e-3, AurynFloat maxweight=1. , TransmitterType transmitter=GLUT); STDPConnection(SpikingGroup * source, NeuronGroup * destination, AurynWeight weight, AurynFloat sparseness=0.05, AurynFloat eta=1, + AurynFloat tau_pre=20e-3, + AurynFloat tau_post=20e-3, AurynFloat maxweight=1. , TransmitterType transmitter=GLUT, string name = "STDPConnection" ); + virtual ~STDPConnection(); virtual void finalize(); void free(); diff --git a/src/System.cpp b/src/System.cpp index e0eb473e..87017ec9 100644 --- a/src/System.cpp +++ b/src/System.cpp @@ -30,6 +30,11 @@ void System::init() { quiet = false; set_simulation_name("default"); + set_online_rate_monitor_tau(); + // assumes that we have at least one spiking group in the sim + online_rate_monitor_id = 0; + online_rate_monitor_state = 0.0; + syncbuffer = new SyncBuffer(mpicom); stringstream oss; @@ -215,6 +220,9 @@ void System::evolve() for ( iter = spiking_groups.begin() ; iter != spiking_groups.end() ; ++iter ) (*iter)->conditional_evolve(); // evolve only if existing on rank + + // update the online rate estimate + evolve_online_rate_monitor(); } void System::evolve_independent() @@ -276,8 +284,8 @@ void System::progressbar ( double fraction, AurynTime clk ) { cout<< percent << "% "<< setiosflags(ios::fixed) << " t=" << time ; - if (checkers.size()) - cout << setprecision(1) << " f=" << checkers[0]->get_property() << " Hz "; + if ( online_rate_monitor_id >= 0 ) + cout << setprecision(1) << " f=" << online_rate_monitor_state << " Hz "; cout << std::flush; @@ -497,6 +505,8 @@ void System::set_simulation_name(string name) void System::save_network_state(string basename) { + logger->msg("Saving network state", NOTIFICATION); + string netstate_filename; { stringstream oss; @@ -509,35 +519,93 @@ void System::save_network_state(string basename) std::ofstream ofs(netstate_filename.c_str()); boost::archive::binary_oarchive oa(ofs); + logger->msg("Saving Connections ...",DEBUG); for ( unsigned int i = 0 ; i < connections.size() ; ++i ) { - // sprintf(filename, "%s.%d.%d.wmat", basename.c_str(), i, mpicom->rank()); stringstream oss; oss << "Saving connection " << i + << " \"" + << connections[i]->get_name() + << "\"" << " to stream"; - logger->msg(oss.str(),NOTIFICATION); + logger->msg(oss.str(),DEBUG); oa << *(connections[i]); } + logger->msg("Saving SpikingGroups",DEBUG); for ( unsigned int i = 0 ; i < spiking_groups.size() ; ++i ) { - // sprintf(filename, "%s.%d.%d.gstate", basename.c_str(), i, mpicom->rank()); stringstream oss; - oss << "Saving SpikingGroup " + oss << "Saving SpikingGroup ..." << i << " to stream"; - logger->msg(oss.str(),NOTIFICATION); + logger->msg(oss.str(),DEBUG); oa << *(spiking_groups[i]); } + // Save Monitors + logger->msg("Saving Monitors ...",DEBUG); + for ( unsigned int i = 0 ; i < monitors.size() ; ++i ) { + + stringstream oss; + oss << "Saving Monitor " + << i + << " to stream"; + logger->msg(oss.str(),DEBUG); + + oa << *(monitors[i]); + } + + logger->msg("Saving Checkers ...",DEBUG); + for ( unsigned int i = 0 ; i < checkers.size() ; ++i ) { + + stringstream oss; + oss << "Saving Checker " + << i + << " to stream"; + logger->msg(oss.str(),DEBUG); + + oa << *(checkers[i]); + } + ofs.close(); } +void System::save_network_state_text(string basename) +{ + logger->msg("Saving network state to textfile", NOTIFICATION); + + char filename [255]; + for ( unsigned int i = 0 ; i < connections.size() ; ++i ) { + sprintf(filename, "%s.%d.%d.wmat", basename.c_str(), i, mpicom->rank()); + + stringstream oss; + oss << "Saving connection " + << filename ; + logger->msg(oss.str(),DEBUG); + + connections[i]->write_to_file(filename); + } + + for ( unsigned int i = 0 ; i < spiking_groups.size() ; ++i ) { + sprintf(filename, "%s.%d.%d.gstate", basename.c_str(), i, mpicom->rank()); + + stringstream oss; + oss << "Saving group " + << filename ; + logger->msg(oss.str(),DEBUG); + + spiking_groups[i]->write_to_file(filename); + } +} + void System::load_network_state(string basename) { + logger->msg("Loading network state", NOTIFICATION); + string netstate_filename; { stringstream oss; @@ -548,32 +616,91 @@ void System::load_network_state(string basename) } // oss goes out of focus std::ifstream ifs(netstate_filename.c_str()); + + if ( !ifs.is_open() ) { + stringstream oss; + oss << "Error opening netstate file: " + << netstate_filename; + logger->msg(oss.str(),ERROR); + throw AurynOpenFileException(); + } + boost::archive::binary_iarchive ia(ifs); + logger->msg("Loading connections ...",DEBUG); for ( unsigned int i = 0 ; i < connections.size() ; ++i ) { stringstream oss; oss << "Loading connection " << i ; - logger->msg(oss.str(),NOTIFICATION); + logger->msg(oss.str(),DEBUG); ia >> *(connections[i]); connections[i]->finalize(); } + logger->msg("Loading SpikingGroups ...",DEBUG); for ( unsigned int i = 0 ; i < spiking_groups.size() ; ++i ) { stringstream oss; oss << "Loading group " << i ; - logger->msg(oss.str(),NOTIFICATION); + logger->msg(oss.str(),DEBUG); ia >> *(spiking_groups[i]); } + // Loading Monitors states + logger->msg("Loading Monitors ...",DEBUG); + for ( unsigned int i = 0 ; i < monitors.size() ; ++i ) { + + stringstream oss; + oss << "Loading Monitor " + << i; + logger->msg(oss.str(),DEBUG); + + ia >> *(monitors[i]); + } + + + logger->msg("Loading Checkers ...",DEBUG); + for ( unsigned int i = 0 ; i < checkers.size() ; ++i ) { + + stringstream oss; + oss << "Loading Checker " + << i; + logger->msg(oss.str(),DEBUG); + + ia >> *(checkers[i]); + } + ifs.close(); } +void System::set_online_rate_monitor_tau(AurynDouble tau) +{ + online_rate_monitor_tau = tau; + online_rate_monitor_mul = exp(-dt/tau); +} + +void System::evolve_online_rate_monitor() +{ + if ( online_rate_monitor_id >= 0 ) { + online_rate_monitor_state *= online_rate_monitor_mul; + SpikingGroup * src = spiking_groups[online_rate_monitor_id]; + online_rate_monitor_state += 1.0*src->get_spikes()->size()/online_rate_monitor_tau/src->get_size(); + } +} + +void System::set_online_rate_monitor_id( int id ) +{ + online_rate_monitor_state = 0.0; + if ( id < spiking_groups.size() ) + online_rate_monitor_id = id; + else + online_rate_monitor_id = -1; +} + #ifdef CODE_COLLECT_SYNC_TIMING_STATS AurynDouble System::get_relative_sync_time() { diff --git a/src/System.h b/src/System.h index 55a3a30f..fe784f6e 100644 --- a/src/System.h +++ b/src/System.h @@ -48,7 +48,13 @@ namespace mpi = boost::mpi; /*! \brief Class that implements system wide variables and methods to manage and run simulations. * - * This Class contains methods to manage and run sets of classes that make up the simulation. In particular it distinguishes between constituents of types SpikingGroup, Connection, Monitor and Checker. A MPI implementation should implement communicators and all that stuff in here. All the constituent object of a single simulation are stored in STL vectors. The methods evolve() and propagate() from each object in these vectors are called alternatingly from within the run procedure. + * This Class contains methods to manage and run sets of classes that make up + * the simulation. In particular it distinguishes between constituents of types + * SpikingGroup, Connection, Monitor and Checker. A MPI implementation should + * implement communicators and all that stuff in here. All the constituent + * object of a single simulation are stored in STL vectors. The methods + * evolve() and propagate() from each object in these vectors are called + * alternatingly from within the run procedure. */ class System { @@ -60,6 +66,7 @@ class System SyncBuffer * syncbuffer; + vector spiking_groups; vector connections; vector monitors; @@ -67,6 +74,14 @@ class System double simulation_time_realtime_ratio; + int online_rate_monitor_id; + double online_rate_monitor_tau; + double online_rate_monitor_mul; + double online_rate_monitor_state; + + /*! Evolves the online rate monitor for the status bar. */ + void evolve_online_rate_monitor(); + /*! Returns string with a human readable time. */ string get_nice_time ( AurynTime clk ); @@ -97,10 +112,49 @@ class System /*! Initialializes the recvs for all the MPI sync */ void sync_prepare(); + /*! Sets the SpikingGroup ID used to display the rate estimate in the + * progressbar (this typically is reflected by the order in + * which you define the SpikingGroup and NeuronGroup classes. It starts + * numbering from 0.). */ + void set_online_rate_monitor_id( int id=0 ); + + /*! Sets the timeconstant to compute the online rate average for the status bar. */ + void set_online_rate_monitor_tau( AurynDouble tau=100e-3 ); + + /*! \brief Saves network state to a netstate file + * + * This function saves the network state to one serialized file. The network + * state includes the internal state variables of all neurons and the synaptic + * connections. It currently does not save the state of any random number + * generators (v0.5) but this is planned to change in the future. Note that + * netstate files do not contain any parameters either. This was done to + * allow to run a simulation with a certain parameter set for a given amount + * of time. Save the network state and then continue the simulation from + * that point with a changed parameter set (e.g. a new stimulus set or + * similar). + * + * \param Prefix (including directory path) of the netstate file without extension + */ void save_network_state(string basename); + /*! \brief Loads network state from a netstate file + * + * \param Basename (directory and prefix of file) of the netstate file without extension + */ void load_network_state(string basename); + /*! \brief Saves the network state to human readable text files + * + * This deprecated method of saving the network state generates a large number of files + * because each Connection object or SpikingGroup creates their own respective file. This + * function might still be useful if you have code in which you analaze these files offline. + * In most cases you will want to use save_network_state and only dump a limited subset (e.g. + * all the plastic connections) in human-readable text files for analysis. + * + * \param Basename (directory and prefix of file) of the netstate file without extension + */ + void save_network_state_text(string basename); + /*! Synchronizes SpikingGroups */ void sync(); diff --git a/src/TripletConnection.cpp b/src/TripletConnection.cpp index 6b6c0b0a..80428e40 100644 --- a/src/TripletConnection.cpp +++ b/src/TripletConnection.cpp @@ -135,8 +135,6 @@ AurynWeight TripletConnection::get_hom(NeuronID i) } -/*! This function implements what happens to synapes transmitting a - * spike to neuron 'post'. */ AurynWeight TripletConnection::dw_pre(NeuronID post) { // translate post id to local id on rank: translated_spike @@ -145,8 +143,6 @@ AurynWeight TripletConnection::dw_pre(NeuronID post) return dw; } -/*! This function implements what happens to synapes experiencing a - * backpropagating action potential from neuron 'pre'. */ AurynWeight TripletConnection::dw_post(NeuronID pre, NeuronID post) { // at this point post was already translated to a local id in @@ -158,21 +154,23 @@ AurynWeight TripletConnection::dw_post(NeuronID pre, NeuronID post) void TripletConnection::propagate_forward() { - // loop over all spikes + // loop over all spikes (yields presynaptic cell ids of cells that spiked) for (SpikeContainer::const_iterator spike = src->get_spikes()->begin() ; // spike = pre_spike spike != src->get_spikes()->end() ; ++spike ) { - // loop over all postsynaptic partners + // loop over all postsynaptic partners the cells + // that are targeted by that presynaptic cell for (const NeuronID * c = w->get_row_begin(*spike) ; c != w->get_row_end(*spike) ; ++c ) { // c = post index - // transmit signal to target at postsynaptic neuron + // determines the weight of connection AurynWeight * weight = w->get_data_ptr(c); + // evokes the postsynaptic response transmit( *c , *weight ); - // handle plasticity + // handles plasticity if ( stdp_active ) { - // performs weight update + // performs weight update upon presynaptic spike *weight += dw_pre(*c); // clips too small weights diff --git a/src/TripletConnection.h b/src/TripletConnection.h index d7b05050..823cb606 100644 --- a/src/TripletConnection.h +++ b/src/TripletConnection.h @@ -48,19 +48,17 @@ class TripletConnection : public DuplexConnection { private: - friend class boost::serialization::access; - template - void save(Archive & ar, const unsigned int version) const + void virtual_serialize(boost::archive::binary_oarchive & ar, const unsigned int version ) { - ar & boost::serialization::base_object(*this); + DuplexConnection::virtual_serialize(ar,version); + ar & *w; } - template - void load(Archive & ar, const unsigned int version) + + void virtual_serialize(boost::archive::binary_iarchive & ar, const unsigned int version ) { - ar & boost::serialization::base_object(*this); - finalize(); + DuplexConnection::virtual_serialize(ar,version); + ar & *w; } - BOOST_SERIALIZATION_SPLIT_MEMBER() void init(AurynFloat tau_hom, AurynFloat eta, AurynFloat kappa, AurynFloat maxweight); void init_shortcuts(); @@ -91,22 +89,35 @@ class TripletConnection : public DuplexConnection DEFAULT_TRACE_MODEL * tr_post2; DEFAULT_TRACE_MODEL * tr_post_hom; + /*! This function propagates spikes from pre to postsynaptic cells + * and performs plasticity updates upon presynaptic spikes. */ void propagate_forward(); + + /*! This performs plasticity updates following postsynaptic spikes. To that end the postsynaptic spikes + * have to be communicated backward to the corresponding synapses connecting to presynaptic neurons. This + * is why this function is called propagate_backward ... it is remeniscent of a back-propagating action + * potential. */ void propagate_backward(); - void sort_spikes(); - /*! Action on weight upon presynaptic spike on connection with postsynaptic - * partner post. This function should be modified to define new spike based - * plasticity rules. - * @param post the postsynaptic cell from which the synaptic trace is read out*/ + + /*! This function implements the plastic update to each + * synapse at the time of a presynaptic spike. + * + * \param post the parameter specifies the postsynaptic partner for which we + * are computing the update. + * */ AurynWeight dw_pre(NeuronID post); - /*! Action on weight upon postsynaptic spike of cell post on connection - * with presynaptic partner pre. This function should be modified to define - * new spike based plasticity rules. - * @param pre the presynaptic cell in question. - * @param post the postsynaptic cell in question. - */ + /*! This function implements the plastic update to each + * synapse at the time of a postsynaptic spike. Since LTP in the minimal triplet model + * depends on the timing of the last pre and postsynaptic spike we are passing both NeuronID + * as arguments. + * + * \param pre The parameter specifies the presynaptic partner for which we + * are computing the update. + * \param post the parameter specifies the postsynaptic partner for which we + * are computing the update. + * */ AurynWeight dw_post(NeuronID pre, NeuronID post); @@ -120,6 +131,8 @@ class TripletConnection : public DuplexConnection TripletConnection(SpikingGroup * source, NeuronGroup * destination, TransmitterType transmitter=GLUT); + /*! Deprecated constructor. + */ TripletConnection(SpikingGroup * source, NeuronGroup * destination, const char * filename, AurynFloat tau_hom=10, diff --git a/src/VoltageMonitor.cpp b/src/VoltageMonitor.cpp index b3fae1b7..cbef69eb 100644 --- a/src/VoltageMonitor.cpp +++ b/src/VoltageMonitor.cpp @@ -41,8 +41,6 @@ void VoltageMonitor::init(NeuronGroup * source, NeuronID id, string filename, Au ssize = stepsize; if ( ssize < 1 ) ssize = 1; - if ( ssize < 1 ) ssize = 1; - nid = id; gid = src->rank2global(nid); paste_spikes = true; @@ -59,11 +57,16 @@ void VoltageMonitor::init(NeuronGroup * source, NeuronID id, string filename, Au void VoltageMonitor::propagate() { if (active && (sys->get_clock())%ssize==0 && sys->get_clock() < tStop && nid < src->get_size() ) { - if ( paste_spikes && - std::find(src->get_spikes_immediate()->begin(), - src->get_spikes_immediate()->end(), gid)!=src->get_spikes_immediate()->end() ) - outfile << (sys->get_time()) << " " << VOLTAGEMONITOR_PASTED_SPIKE_HEIGHT << "\n"; - else - outfile << (sys->get_time()) << " " << src->get_mem(nid) << "\n"; + double voltage = src->get_mem(nid); + if ( paste_spikes ) { + SpikeContainer * spikes = src->get_spikes_immediate(); + for ( int i = 0 ; i < spikes->size() ; ++i ) { + if ( spikes->at(i) == gid ) { + voltage = VOLTAGEMONITOR_PASTED_SPIKE_HEIGHT; + break; + } + } + } + outfile << (sys->get_time()) << " " << voltage << "\n"; } } diff --git a/src/WeightMonitor.cpp b/src/WeightMonitor.cpp index 4e4319d2..4583c686 100644 --- a/src/WeightMonitor.cpp +++ b/src/WeightMonitor.cpp @@ -48,7 +48,7 @@ WeightMonitor::WeightMonitor(SparseConnection * source, NeuronID i, NeuronID j, switch (recordingmode) { case DATARANGE : for (AurynLong c = elem_i ; c < elem_j ; ++c) - add_to_list(mat->get_data_begin()+c) ; + add_to_list(c) ; break; case SINGLE : add_to_list(i,j); @@ -80,22 +80,27 @@ void WeightMonitor::init(SparseConnection * source, NeuronID i, NeuronID j, stri // default behavior recordingmode = ELEMENTLIST; - element_list = new vector; + element_list = new vector; group_indices.push_back(0); // important for group mode elem_i = 0; elem_j = 0; } +void WeightMonitor::add_to_list(AurynLong data_index) +{ + element_list->push_back( data_index ); +} + void WeightMonitor::add_to_list(AurynWeight * ptr) { - if ( ptr != NULL ) - element_list->push_back( ptr ); + if ( ptr != NULL ) { + element_list->push_back( mat->get_data_index(ptr) ); + } } void WeightMonitor::add_to_list(NeuronID i, NeuronID j) { - AurynWeight * ptr = mat->get_ptr(i,j); - add_to_list(ptr); + add_to_list( mat->get_data_index(i,j) ); } void WeightMonitor::add_to_list( vector vec, string label ) @@ -125,6 +130,12 @@ void WeightMonitor::add_to_list( vector vec, string label ) void WeightMonitor::add_equally_spaced(NeuronID number) { + if ( number > src->get_nonzero() ) { + logger->msg("WeightMonitor:: add_equally_spaced: \ + Not enough elements in this Connection object",WARNING); + number = src->get_nonzero(); + } + for ( NeuronID i = 0 ; i < number ; ++i ) add_to_list(mat->get_data_begin()+i*mat->get_nonzero()/number); @@ -150,7 +161,7 @@ void WeightMonitor::load_data_range( NeuronID i, NeuronID j ) << j; logger->msg(oss.str(),DEBUG); for ( NeuronID a = i ; a < j ; ++a ) - element_list->push_back( mat->get_data_begin()+a ); + element_list->push_back( a ); outfile << "# Added data range " << i << "-" << j << "." << endl; } @@ -248,9 +259,10 @@ void WeightMonitor::load_pattern_connections( string filename_pre, string filena neuron_pair p; p.i = patterns_pre->at(i)[k].i; p.j = patterns_post->at(j)[l].i; - AurynWeight * ptr = mat->get_ptr(p.i,p.j); - if ( ptr != NULL ) // make sure we are counting connections that do exist + AurynWeight * ptr = mat->get_ptr(p.i,p.j); + if ( ptr != NULL ) { // make sure we are counting connections that do exist list.push_back( p ); + } if ( list.size() >= maxcon ) break; } if ( list.size() >= maxcon ) break; @@ -284,8 +296,8 @@ void WeightMonitor::load_pattern_connections( string filename_pre, string filena void WeightMonitor::record_single_synapses() { - for (vector::iterator iter = element_list->begin() ; iter != element_list->end() ; ++iter) - outfile << *(*iter) << " "; + for (vector::iterator iter = element_list->begin() ; iter != element_list->end() ; ++iter) + outfile << mat->get_data( (*iter) ) << " "; } void WeightMonitor::record_synapse_groups() @@ -295,8 +307,8 @@ void WeightMonitor::record_synapse_groups() AurynDouble sum2 = 0; for ( int k = group_indices[i-1] ; k < group_indices[i] ; ++k ) { - sum += *(element_list->at(k)); - sum2 += pow((*(element_list->at(k))),2); + sum += mat->get_data( element_list->at(k) ); + sum2 += pow( mat->get_data( element_list->at(k) ), 2); } NeuronID n = group_indices[i]-group_indices[i-1]; AurynDouble mean = sum/n; diff --git a/src/WeightMonitor.h b/src/WeightMonitor.h index 27a06c4a..aa030d35 100644 --- a/src/WeightMonitor.h +++ b/src/WeightMonitor.h @@ -75,7 +75,7 @@ class WeightMonitor : protected Monitor NeuronID elem_i; NeuronID elem_j; AurynTime ssize; - vector * element_list; + vector * element_list; vector group_indices; void init(SparseConnection * source, NeuronID i, NeuronID j, string filename, AurynTime interval); @@ -94,7 +94,9 @@ class WeightMonitor : protected Monitor void set_mat(ForwardMatrix * m); - /*! Adds a single element identified by a pointer to the recording list. */ + /*! Adds a single element to the recording list which is identified by its data index. */ + void add_to_list( AurynLong index ); + /*! Adds a single element to the recording list which is identified by a pointer. */ void add_to_list( AurynWeight * ptr ); /*! Adds a single element identified matrix coordinates (row,col) to the recording list. */ void add_to_list( NeuronID i, NeuronID j ); diff --git a/src/auryn_definitions.h b/src/auryn_definitions.h index 0975fd4d..7ce410cd 100644 --- a/src/auryn_definitions.h +++ b/src/auryn_definitions.h @@ -63,8 +63,8 @@ namespace mpi = boost::mpi; /*! The current Auryn version number */ -#define AURYNVERSION 0.5 -#define AURYNSUBVERSION 0 +#define AURYNVERSION 0.6 +#define AURYNSUBVERSION 1 /*! Toggle between memory alignment for * SIMD code. @@ -127,7 +127,8 @@ enum TransmitterType { GABA, //!< Standard Gabaergic (inhibitory) transmission. AMPA, //!< Only targets AMPA channels. NMDA, //!< Only targets NMDA. - MEM //!< Current based synapse. Adds the transmitted quantity directly to membrane voltage. + MEM, //!< Current based synapse. Adds the transmitted quantity directly to membrane voltage. + CURSYN //!< Current based synapse with dynamics. }; enum StimulusGroupModeType { MANUAL, RANDOM, SEQUENTIAL, SEQUENTIAL_REV, STIMFILE }; @@ -327,5 +328,21 @@ class AurynTimeOverFlowException: public exception } }; +class AurynStateVectorException: public exception +{ + virtual const char* what() const throw() + { + return "Auryn encountered an undefined problem when dealing with StateVectors."; + } +}; + +class AurynGenericException: public exception +{ + virtual const char* what() const throw() + { + return "Auryn encountered a problem which it deemed serious enough to break the run. \ + To debug set logger vebosity to EVERYTHING and analyze the log files."; + } +}; #endif /*AURYN_DEFINITIONS_H__*/