diff --git a/.gitignore b/.gitignore
index fdb1db50d0eddc2e8190301f5a785944b261b246..57389ad9a432bf12e26fee407bf7c969923765c4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,3 +66,6 @@ commit.msg
# eclipse remote sync folders
.ptp-sync-folder
+
+# Mac default files
+.DS_Store
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 4730d3d18565649a951cad791dc007cf76f10f5f..31e7122445088a9258d64958c05a763d31c6b19c 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(miniapp)
add_subdirectory(generators)
+add_subdirectory(brunel)
diff --git a/example/brunel/CMakeLists.txt b/example/brunel/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a7ed3fd070661657090b2316941f6ebbd4591a21
--- /dev/null
+++ b/example/brunel/CMakeLists.txt
@@ -0,0 +1,23 @@
+set(HEADERS
+ io.hpp
+ partitioner.hpp
+)
+set(MINIAPP_SOURCES
+ brunel_miniapp.cpp
+ io.cpp
+)
+
+add_executable(brunel_miniapp.exe ${MINIAPP_SOURCES} ${HEADERS})
+
+target_link_libraries(brunel_miniapp.exe LINK_PUBLIC ${ARB_LIBRARIES})
+target_link_libraries(brunel_miniapp.exe LINK_PUBLIC ${EXTERNAL_LIBRARIES})
+
+if(ARB_WITH_MPI)
+ target_link_libraries(brunel_miniapp.exe LINK_PUBLIC ${MPI_C_LIBRARIES})
+ set_property(TARGET brunel_miniapp.exe APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}")
+endif()
+
+set_target_properties(brunel_miniapp.exe
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/example/miniapp/brunel"
+)
diff --git a/example/brunel/brunel_miniapp.cpp b/example/brunel/brunel_miniapp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea6ca86d0821b016a77aaacc01aeeadf148c6143
--- /dev/null
+++ b/example/brunel/brunel_miniapp.cpp
@@ -0,0 +1,319 @@
+#include <cmath>
+#include <exception>
+#include <event_generator.hpp>
+#include <iostream>
+#include <fstream>
+#include <memory>
+#include <vector>
+#include <json/json.hpp>
+#include <common_types.hpp>
+#include <communication/communicator.hpp>
+#include <communication/global_policy.hpp>
+#include <io/exporter_spike_file.hpp>
+#include <lif_cell_description.hpp>
+#include <model.hpp>
+#include "partitioner.hpp"
+#include <profiling/profiler.hpp>
+#include <profiling/meter_manager.hpp>
+#include <recipe.hpp>
+#include <set>
+#include <threading/threading.hpp>
+#include <util/config.hpp>
+#include <util/debug.hpp>
+#include <util/ioutil.hpp>
+#include <util/nop.hpp>
+#include <vector>
+#include "io.hpp"
+#include <hardware/gpu.hpp>
+#include <hardware/node_info.hpp>
+
+using namespace arb;
+
+using global_policy = communication::global_policy;
+using file_export_type = io::exporter_spike_file<global_policy>;
+void banner(hw::node_info);
+using communicator_type = communication::communicator<communication::global_policy>;
+
+// Samples m unique values in interval [start, end) - gid.
+// We exclude gid because we don't want self-loops.
+std::vector<cell_gid_type> sample_subset(cell_gid_type gid, cell_gid_type start, cell_gid_type end, unsigned m) {
+ std::set<cell_gid_type> s;
+ std::mt19937 gen(gid + 42);
+ std::uniform_int_distribution<cell_gid_type> dis(start, end - 1);
+ while (s.size() < m) {
+ auto val = dis(gen);
+ if (val != gid) {
+ s.insert(val);
+ }
+ }
+ return {s.begin(), s.end()};
+}
+
+/*
+ A Brunel network consists of nexc excitatory LIF neurons and ninh inhibitory LIF neurons.
+ Each neuron in the network receives in_degree_prop * nexc excitatory connections
+ chosen randomly, in_degree_prop * ninh inhibitory connections and next (external) Poisson connections.
+ All the connections have the same delay. The strenght of excitatory and Poisson connections is given by
+ parameter weight, whereas the strength of inhibitory connections is rel_inh_strength * weight.
+ Poisson neurons all spike independently with expected number of spikes given by parameter poiss_lambda.
+ Because of the refractory period, the activity is mostly driven by Poisson neurons and
+ recurrent connections have a small effect.
+ */
+class brunel_recipe: public recipe {
+public:
+ brunel_recipe(cell_size_type nexc, cell_size_type ninh, cell_size_type next, double in_degree_prop,
+ float weight, float delay, float rel_inh_strength, double poiss_lambda, int seed = 42):
+ ncells_exc_(nexc), ncells_inh_(ninh), ncells_ext_(next), delay_(delay), seed_(seed) {
+ // Make sure that in_degree_prop in the interval (0, 1]
+ if (in_degree_prop <= 0.0 || in_degree_prop > 1.0) {
+ std::out_of_range("The proportion of incoming connections should be in the interval (0, 1].");
+ }
+
+ // Set up the parameters.
+ weight_exc_ = weight;
+ weight_inh_ = -rel_inh_strength * weight_exc_;
+ weight_ext_ = weight;
+ in_degree_exc_ = std::round(in_degree_prop * nexc);
+ in_degree_inh_ = std::round(in_degree_prop * ninh);
+ // each cell receives next incoming Poisson sources with mean rate poiss_lambda, which is equivalent
+ // to a single Poisson source with mean rate next*poiss_lambda
+ lambda_ = next * poiss_lambda;
+ }
+
+ cell_size_type num_cells() const override {
+ return ncells_exc_ + ncells_inh_;
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::lif_neuron;
+ }
+
+ std::vector<cell_connection> connections_on(cell_gid_type gid) const override {
+ std::vector<cell_connection> connections;
+ // Add incoming excitatory connections.
+ for (auto i: sample_subset(gid, 0, ncells_exc_, in_degree_exc_)) {
+ cell_member_type source{cell_gid_type(i), 0};
+ cell_member_type target{gid, 0};
+ cell_connection conn(source, target, weight_exc_, delay_);
+ connections.push_back(conn);
+ }
+
+ // Add incoming inhibitory connections.
+ for (auto i: sample_subset(gid, ncells_exc_, ncells_exc_ + ncells_inh_, in_degree_inh_)) {
+ cell_member_type source{cell_gid_type(i), 0};
+ cell_member_type target{gid, 0};
+ cell_connection conn(source, target, weight_inh_, delay_);
+ connections.push_back(conn);
+ }
+ return connections;
+ }
+
+ util::unique_any get_cell_description(cell_gid_type gid) const override {
+ auto cell = lif_cell_description();
+ cell.tau_m = 10;
+ cell.V_th = 10;
+ cell.C_m = 20;
+ cell.E_L = 0;
+ cell.V_m = 0;
+ cell.V_reset = 0;
+ cell.t_ref = 2;
+ return cell;
+ }
+
+ std::vector<event_generator_ptr> event_generators(cell_gid_type gid) const override {
+ std::vector<arb::event_generator_ptr> gens;
+
+ std::mt19937_64 G;
+ G.seed(gid + seed_);
+
+ using pgen = poisson_generator<std::mt19937_64>;
+
+ time_type t0 = 0;
+ cell_member_type target{gid, 0};
+
+ gens.push_back(arb::make_event_generator<pgen>(
+ target,
+ weight_ext_,
+ G,
+ t0,
+ lambda_));
+ return gens;
+ }
+
+ cell_size_type num_sources(cell_gid_type) const override {
+ return 1;
+ }
+
+ cell_size_type num_targets(cell_gid_type) const override {
+ return 1;
+ }
+
+ cell_size_type num_probes(cell_gid_type) const override {
+ return 0;
+ }
+
+ probe_info get_probe(cell_member_type probe_id) const override {
+ return {};
+ }
+
+private:
+ // Number of excitatory cells.
+ cell_size_type ncells_exc_;
+
+ // Number of inhibitory cells.
+ cell_size_type ncells_inh_;
+
+ // Number of Poisson connections each neuron receives
+ cell_size_type ncells_ext_;
+
+ // Weight of excitatory synapses.
+ float weight_exc_;
+
+ // Weight of inhibitory synapses.
+ float weight_inh_;
+
+ // Weight of external Poisson cell synapses.
+ float weight_ext_;
+
+ // Delay of all synapses.
+ float delay_;
+
+ // Number of connections that each neuron receives from excitatory population.
+ int in_degree_exc_;
+
+ // Number of connections that each neuron receives from inhibitory population.
+ int in_degree_inh_;
+
+ // Expected number of poisson spikes.
+ double lambda_;
+
+ // Seed used for the Poisson spikes generation.
+ int seed_;
+};
+
+using util::any_cast;
+using util::make_span;
+using global_policy = communication::global_policy;
+using file_export_type = io::exporter_spike_file<global_policy>;
+
+int main(int argc, char** argv) {
+ arb::communication::global_policy_guard global_guard(argc, argv);
+ try {
+ arb::util::meter_manager meters;
+ meters.start();
+ std::cout << util::mask_stream(global_policy::id()==0);
+ // read parameters
+ io::cl_options options = io::read_options(argc, argv, global_policy::id()==0);
+ hw::node_info nd;
+ nd.num_cpu_cores = threading::num_threads();
+ nd.num_gpus = hw::num_gpus()>0? 1: 0;
+ banner(nd);
+
+ meters.checkpoint("setup");
+
+ // The size of excitatory population.
+ cell_size_type nexc = options.nexc;
+
+ // The size of inhibitory population.
+ cell_size_type ninh = options.ninh;
+
+ // The size of Poisson (external) population.
+ cell_size_type next = options.next;
+
+ // Fraction of connections each neuron receives from each of the 3 populations.
+ double in_degree_prop = options.syn_per_cell_prop;
+
+ // Weight of excitatory and poisson connections.
+ float w = options.weight;
+
+ // Delay of all the connections.
+ float d = options.delay;
+
+ // Relative strength of inhibitory connections with respect to excitatory connections.
+ float rel_inh_strength = options.rel_inh_strength;
+
+ // Expected number of spikes from a single poisson cell per ms.
+ int poiss_lambda = options.poiss_lambda;
+
+ // The number of cells in a single cell group.
+ cell_size_type group_size = options.group_size;
+
+ unsigned seed = options.seed;
+
+ brunel_recipe recipe(nexc, ninh, next, in_degree_prop, w, d, rel_inh_strength, poiss_lambda, seed);
+
+ auto register_exporter = [] (const io::cl_options& options) {
+ return util::make_unique<file_export_type>
+ (options.file_name, options.output_path,
+ options.file_extension, options.over_write);
+ };
+
+ auto decomp = decompose(recipe, group_size);
+ model m(recipe, decomp);
+
+ // Initialize the spike exporting interface
+ std::unique_ptr<file_export_type> file_exporter;
+ if (options.spike_file_output) {
+ if (options.single_file_per_rank) {
+ file_exporter = register_exporter(options);
+
+ m.set_local_spike_callback(
+ [&](const std::vector<spike>& spikes) {
+ file_exporter->output(spikes);
+ }
+ );
+ }
+ else if(communication::global_policy::id()==0) {
+ file_exporter = register_exporter(options);
+
+ m.set_global_spike_callback(
+ [&](const std::vector<spike>& spikes) {
+ file_exporter->output(spikes);
+ }
+ );
+ }
+ }
+ meters.checkpoint("model-init");
+
+ // run model
+ m.run(options.tfinal, options.dt);
+
+ meters.checkpoint("model-simulate");
+
+ // output profile and diagnostic feedback
+ util::profiler_output(0.001, options.profile_only_zero);
+ std::cout << "There were " << m.num_spikes() << " spikes\n";
+
+ auto report = util::make_meter_report(meters);
+ std::cout << report;
+ if (global_policy::id()==0) {
+ std::ofstream fid;
+ fid.exceptions(std::ios_base::badbit | std::ios_base::failbit);
+ fid.open("meters.json");
+ fid << std::setw(1) << util::to_json(report) << "\n";
+ }
+ }
+ catch (io::usage_error& e) {
+ // only print usage/startup errors on master
+ std::cerr << util::mask_stream(global_policy::id()==0);
+ std::cerr << e.what() << "\n";
+ return 1;
+ }
+ catch (std::exception& e) {
+ std::cerr << e.what() << "\n";
+ return 2;
+ }
+ return 0;
+}
+
+void banner(hw::node_info nd) {
+ std::cout << "==========================================\n";
+ std::cout << " Arbor miniapp\n";
+ std::cout << " - distributed : " << global_policy::size()
+ << " (" << std::to_string(global_policy::kind()) << ")\n";
+ std::cout << " - threads : " << nd.num_cpu_cores
+ << " (" << threading::description() << ")\n";
+ std::cout << " - gpus : " << nd.num_gpus << "\n";
+ std::cout << "==========================================\n";
+}
+
diff --git a/example/brunel/io.cpp b/example/brunel/io.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..afb406761a6f35ed76aa62bfdbeff1f4d1d8afba
--- /dev/null
+++ b/example/brunel/io.cpp
@@ -0,0 +1,206 @@
+#include <algorithm>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <sstream>
+#include <string>
+#include <type_traits>
+#include <tclap/CmdLine.h>
+#include <util/meta.hpp>
+#include <util/optional.hpp>
+#include "io.hpp"
+
+// Let TCLAP understand value arguments that are of an optional type.
+namespace TCLAP {
+ template <typename V>
+ struct ArgTraits<arb::util::optional<V>> {
+ using ValueCategory = ValueLike;
+ };
+} // namespace TCLAP
+
+namespace arb {
+ namespace util {
+ // Using static here because we do not want external linkage for this operator.
+ template <typename V>
+ static std::istream& operator>>(std::istream& I, optional<V>& v) {
+ V u;
+ if (I >> u) {
+ v = u;
+ }
+ return I;
+ }
+ }
+
+ namespace io {
+ // Override annoying parameters listed back-to-front behaviour.
+ //
+ // TCLAP argument creation _prepends_ its arguments to the internal
+ // list (_argList), where standard options --help etc. are already
+ // pre-inserted.
+ //
+ // reorder_arguments() reverses the arguments to restore ordering,
+ // and moves the standard options to the end.
+ class CustomCmdLine: public TCLAP::CmdLine {
+ public:
+ CustomCmdLine(const std::string &message, const std::string &version = "none"):
+ TCLAP::CmdLine(message, ' ', version, true)
+ {}
+
+ void reorder_arguments() {
+ _argList.reverse();
+ for (auto opt: {"help", "version", "ignore_rest"}) {
+ auto i = std::find_if(
+ _argList.begin(), _argList.end(),
+ [&opt](TCLAP::Arg* a) { return a->getName()==opt; });
+
+ if (i!=_argList.end()) {
+ auto a = *i;
+ _argList.erase(i);
+ _argList.push_back(a);
+ }
+ }
+ }
+ };
+
+ // Update an option value from command line argument if set.
+ template <
+ typename T,
+ typename Arg,
+ typename = util::enable_if_t<std::is_base_of<TCLAP::Arg, Arg>::value>
+ >
+ static void update_option(T& opt, Arg& arg) {
+ if (arg.isSet()) {
+ opt = arg.getValue();
+ }
+ }
+
+ // Read options from (optional) json file and command line arguments.
+ cl_options read_options(int argc, char** argv, bool allow_write) {
+ cl_options options;
+ std::string save_file = "";
+
+ // Parse command line arguments.
+ try {
+ cl_options defopts;
+
+ CustomCmdLine cmd("nest brunel miniapp harness", "0.1");
+
+ TCLAP::ValueArg<uint32_t> nexc_arg
+ ("n", "n-excitatory", "total number of cells in the excitatory population",
+ false, defopts.nexc, "integer", cmd);
+
+ TCLAP::ValueArg<uint32_t> ninh_arg
+ ("m", "n-inhibitory", "total number of cells in the inhibitory population",
+ false, defopts.ninh, "integer", cmd);
+
+ TCLAP::ValueArg<uint32_t> next_arg
+ ("e", "n-external", "total number of incoming Poisson (external) connections per cell.",
+ false, defopts.ninh, "integer", cmd);
+
+ TCLAP::ValueArg<double> syn_prop_arg
+ ("p", "in-degree-prop", "the proportion of connections both the excitatory and inhibitory populations that each neuron receives",
+ false, defopts.syn_per_cell_prop, "double", cmd);
+
+ TCLAP::ValueArg<float> weight_arg
+ ("w", "weight", "the weight of all excitatory connections",
+ false, defopts.weight, "float", cmd);
+
+ TCLAP::ValueArg<float> delay_arg
+ ("d", "delay", "the delay of all connections",
+ false, defopts.delay, "float", cmd);
+
+ TCLAP::ValueArg<float> rel_inh_strength_arg
+ ("g", "rel-inh-w", "relative strength of inhibitory synapses with respect to the excitatory ones",
+ false, defopts.rel_inh_strength, "float", cmd);
+
+ TCLAP::ValueArg<double> poiss_lambda_arg
+ ("l", "lambda", "Expected number of spikes from a single poisson cell per ms",
+ false, defopts.poiss_lambda, "double", cmd);
+
+ TCLAP::ValueArg<double> tfinal_arg
+ ("t", "tfinal", "length of the simulation period [ms]",
+ false, defopts.tfinal, "time", cmd);
+
+ TCLAP::ValueArg<double> dt_arg
+ ("s", "delta-t", "simulation time step [ms] (this parameter is ignored)",
+ false, defopts.dt, "time", cmd);
+
+ TCLAP::ValueArg<uint32_t> group_size_arg
+ ("G", "group-size", "number of cells per cell group",
+ false, defopts.group_size, "integer", cmd);
+
+ TCLAP::ValueArg<uint32_t> seed_arg
+ ("S", "seed", "seed for poisson spike generators",
+ false, defopts.seed, "integer", cmd);
+
+ TCLAP::SwitchArg spike_output_arg
+ ("f","spike-file-output","save spikes to file", cmd, false);
+
+ TCLAP::SwitchArg profile_only_zero_arg
+ ("z", "profile-only-zero", "Only output profile information for rank 0",
+ cmd, false);
+
+ TCLAP::SwitchArg verbose_arg
+ ("v", "verbose", "Present more verbose information to stdout", cmd, false);
+
+ cmd.reorder_arguments();
+ cmd.parse(argc, argv);
+
+ // Handle verbosity separately from other options: it is not considered part
+ // of the saved option state.
+ options.verbose = verbose_arg.getValue();
+ update_option(options.nexc, nexc_arg);
+ update_option(options.ninh, ninh_arg);
+ update_option(options.next, next_arg);
+ update_option(options.syn_per_cell_prop, syn_prop_arg);
+ update_option(options.weight, weight_arg);
+ update_option(options.delay, delay_arg);
+ update_option(options.rel_inh_strength, rel_inh_strength_arg);
+ update_option(options.poiss_lambda, poiss_lambda_arg);
+ update_option(options.tfinal, tfinal_arg);
+ update_option(options.dt, dt_arg);
+ update_option(options.group_size, group_size_arg);
+ update_option(options.seed, seed_arg);
+ update_option(options.spike_file_output, spike_output_arg);
+ update_option(options.profile_only_zero, profile_only_zero_arg);
+
+ if (options.group_size < 1) {
+ throw usage_error("minimum of one cell per group");
+ }
+
+ if (options.rel_inh_strength <= 0 || options.rel_inh_strength > 1) {
+ throw usage_error("relative strength of inhibitory connections must be in the interval (0, 1].");
+ }
+ }
+ catch (TCLAP::ArgException& e) {
+ throw usage_error("error parsing command line argument "+e.argId()+": "+e.error());
+ }
+
+ // If verbose output requested, emit option summary.
+ if (options.verbose) {
+ std::cout << options << "\n";
+ }
+
+ return options;
+ }
+
+ std::ostream& operator<<(std::ostream& o, const cl_options& options) {
+ o << "simulation options:\n";
+ o << " excitatory cells : " << options.nexc << "\n";
+ o << " inhibitory cells : " << options.ninh << "\n";
+ o << " Poisson connections per cell : " << options.next << "\n";
+ o << " proportion of synapses/cell from each population : " << options.syn_per_cell_prop << "\n";
+ o << " weight of excitatory synapses : " << options.weight << "\n";
+ o << " relative strength of inhibitory synapses : " << options.rel_inh_strength << "\n";
+ o << " delay of all synapses : " << options.delay << "\n";
+ o << " expected number of spikes from a single poisson cell per ms: " << options.poiss_lambda << "\n";
+ o << "\n";
+ o << " simulation time : " << options.tfinal << "\n";
+ o << " dt : " << options.dt << "\n";
+ o << " group size : " << options.group_size << "\n";
+ o << " seed : " << options.seed << "\n";
+ return o;
+ }
+ } // namespace io
+} // namespace arbor
diff --git a/example/brunel/io.hpp b/example/brunel/io.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..45c0c65c7d1b90c4b7201866a82281d0dc3d535e
--- /dev/null
+++ b/example/brunel/io.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+#include <common_types.hpp>
+#include <util/optional.hpp>
+#include <util/path.hpp>
+
+namespace arb {
+namespace io {
+ // Holds the options for a simulation run.
+ // Default constructor gives default options.
+ struct cl_options {
+ // Cell parameters:
+ uint32_t nexc = 400;
+ uint32_t ninh = 100;
+ uint32_t next = 40;
+ double syn_per_cell_prop = 0.05;
+ float weight = 1.2;
+ float delay = 0.1;
+ float rel_inh_strength = 1;
+ double poiss_lambda = 1;
+
+ // Simulation running parameters:
+ double tfinal = 100.;
+ double dt = 1;
+ uint32_t group_size = 10;
+ uint32_t seed = 42;
+
+ // Parameters for spike output.
+ bool spike_file_output = false;
+ bool single_file_per_rank = false;
+ bool over_write = true;
+ std::string output_path = "./";
+ std::string file_name = "spikes";
+ std::string file_extension = "gdf";
+
+ // Turn on/off profiling output for all ranks.
+ bool profile_only_zero = false;
+
+ // Be more verbose with informational messages.
+ bool verbose = false;
+ };
+
+ class usage_error: public std::runtime_error {
+ public:
+ template <typename S>
+ usage_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
+ };
+
+ class model_description_error: public std::runtime_error {
+ public:
+ template <typename S>
+ model_description_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
+ };
+
+ std::ostream& operator<<(std::ostream& o, const cl_options& opt);
+
+ cl_options read_options(int argc, char** argv, bool allow_write = true);
+} // namespace io
+} // namespace arbor
diff --git a/example/brunel/partitioner.hpp b/example/brunel/partitioner.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f56c0915b443a06284695669c5fcddf7fb94904
--- /dev/null
+++ b/example/brunel/partitioner.hpp
@@ -0,0 +1,66 @@
+#include <communication/global_policy.hpp>
+#include <domain_decomposition.hpp>
+#include <hardware/node_info.hpp>
+#include <recipe.hpp>
+
+namespace arb {
+ domain_decomposition decompose(const recipe& rec, const unsigned group_size) {
+ struct partition_gid_domain {
+ partition_gid_domain(std::vector<cell_gid_type> divs):
+ gid_divisions(std::move(divs))
+ {}
+
+ int operator()(cell_gid_type gid) const {
+ auto gid_part = util::partition_view(gid_divisions);
+ return gid_part.index(gid);
+ }
+
+ const std::vector<cell_gid_type> gid_divisions;
+ };
+
+ cell_size_type num_global_cells = rec.num_cells();
+ unsigned num_domains = communication::global_policy::size();
+ int domain_id = communication::global_policy::id();
+
+ auto dom_size = [&](unsigned dom) -> cell_gid_type {
+ const cell_gid_type B = num_global_cells/num_domains;
+ const cell_gid_type R = num_global_cells - num_domains*B;
+ return B + (dom<R);
+ };
+
+ // Global load balance
+ std::vector<cell_gid_type> gid_divisions;
+ auto gid_part = make_partition(
+ gid_divisions, util::transform_view(util::make_span(0, num_domains), dom_size));
+
+ auto range = gid_part[domain_id];
+ cell_size_type num_local_cells = range.second - range.first;
+
+ unsigned num_groups = num_local_cells / group_size + (num_local_cells%group_size== 0 ? 0 : 1);
+ std::vector<group_description> groups;
+
+ // Local load balance
+ // i.e. all the groups that the current rank (domain) owns
+ for (unsigned i = 0; i < num_groups; ++i) {
+ unsigned start = i * group_size;
+ unsigned end = std::min(start + group_size, num_local_cells);
+ std::vector<cell_gid_type> group_elements;
+
+ for (unsigned j = start; j < end; ++j) {
+ group_elements.push_back(j);
+ }
+
+ groups.push_back({cell_kind::lif_neuron, std::move(group_elements), backend_kind::multicore});
+ }
+
+ domain_decomposition d;
+ d.num_domains = num_domains;
+ d.domain_id = domain_id;
+ d.num_local_cells = num_local_cells;
+ d.num_global_cells = num_global_cells;
+ d.groups = std::move(groups);
+ d.gid_domain = partition_gid_domain(std::move(gid_divisions));
+
+ return d;
+ }
+}
diff --git a/example/brunel/readme.md b/example/brunel/readme.md
new file mode 100644
index 0000000000000000000000000000000000000000..eb1e2f77d0be33332d22db01d99d5cec56a00e93
--- /dev/null
+++ b/example/brunel/readme.md
@@ -0,0 +1,35 @@
+## Miniapp for simulating the Brunel network of LIF neurons.
+
+The network consists of 2 main populations: excitatory and inhibitory populations.
+Each neuron from the network receives a fixed number (proportional to the size
+of the population) of incoming connections from both of these groups.
+In addition to the input from excitatory and inhibitory populations,
+each neuron receives a fixed number of connections from the Poisson neurons
+producing the Poisson-like input that is integrated into the LIF cell group so
+that the communication of Poisson events is bypassed.
+
+### Parameters
+
+The parameters that can be passed as command-line arguments are the following:
+
+* `-n` (`--n_excitatory`): total number of cells in the excitatory population.
+* `-m` (`--n_inhibitory`): total number of cells in the inhibitory population.
+* `-e` (`--n_external`): number of incoming Poisson sources each excitatory and inhibitory neuron receives.
+* `-p` (`--in_degree_prop`): the proportions of connections from both populations that each neurons receives.
+* `-w` (`--weight`): weight of all excitatory connections.
+* `-g` (`--rel_inh_w`): relative strength of inhibitory synapses with respect to the excitatory ones.
+* `-d` (`--delay`): the delay of all connections.
+* `-l` (`--lambda`): rate of Poisson cells (kHz).
+* `-t` (`--tfinal`): length of the simulation period (ms).
+* `-s` (`--delta_t`): simulation time step (ms). (this parameter is ignored)
+* `-G` (`--group-size`): number of cells per cell group
+* `-S` (`--seed`): seed of the Poisson sources attached to cells.
+* `-f` (`--spike-file-output`): save spikes to file (Bool).
+* `-z` (`--profile-only-zero`): only output profile information for rank 0.
+* `-v` (`--verbose`): present more verbose information to stdout.
+
+For example, we could run the miniapp as follows:
+
+```
+./example/miniapp/brunel/brunel_miniapp.exe -n 400 -m 100 -e 20 -p 0.1 -w 1.2 -d 1 -g 0.5 -l 5 -t 100 -s 1 -G 50 -S 123 -f
+```
diff --git a/example/miniapp/CMakeLists.txt b/example/miniapp/CMakeLists.txt
index 70bfd6176610ab77f4255fb75d9ef951e7420256..b309737a8b253fe8927189a88da1ba979edbd23c 100644
--- a/example/miniapp/CMakeLists.txt
+++ b/example/miniapp/CMakeLists.txt
@@ -21,3 +21,4 @@ set_target_properties(
PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/example"
)
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fd29c43e070fbac8d35fe9e1bbc45934a1e72838..348f47a57e8c9f5707dfbe9a998f42a20da2cce7 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,7 @@ set(BASE_SOURCES
common_types_io.cpp
cell.cpp
event_binner.cpp
+ lif_cell_group.cpp
hardware/affinity.cpp
hardware/gpu.cpp
hardware/memory.cpp
diff --git a/src/cell_group_factory.cpp b/src/cell_group_factory.cpp
index a1e5fc1e9e13a8d7595878e5ffb0bb84e31820d7..f1f9711fc8c0b07a0cf0c5a0ac4f860d368d57ba 100644
--- a/src/cell_group_factory.cpp
+++ b/src/cell_group_factory.cpp
@@ -5,6 +5,7 @@
#include <domain_decomposition.hpp>
#include <dss_cell_group.hpp>
#include <fvm_multicell.hpp>
+#include <lif_cell_group.hpp>
#include <mc_cell_group.hpp>
#include <recipe.hpp>
#include <rss_cell_group.hpp>
@@ -28,6 +29,9 @@ cell_group_ptr cell_group_factory(const recipe& rec, const group_description& gr
case cell_kind::regular_spike_source:
return make_cell_group<rss_cell_group>(group.gids, rec);
+ case cell_kind::lif_neuron:
+ return make_cell_group<lif_cell_group>(group.gids, rec);
+
case cell_kind::data_spike_source:
return make_cell_group<dss_cell_group>(group.gids, rec);
diff --git a/src/common_types.hpp b/src/common_types.hpp
index 19a5aacf1c7dfb91a10a1949cf3dd012897c8fda..41284522a6c12286a0dbc25024f6531d0d505642 100644
--- a/src/common_types.hpp
+++ b/src/common_types.hpp
@@ -69,6 +69,7 @@ using sample_size_type = std::int32_t;
enum cell_kind {
cable1d_neuron, // Our own special mc neuron
+ lif_neuron, // Leaky-integrate and fire neuron
regular_spike_source, // Regular spiking source
data_spike_source, // Spike source from values inserted via description
};
diff --git a/src/common_types_io.cpp b/src/common_types_io.cpp
index 9e478115f053a4774c0a26b86bda20182966e986..3a3fa34e7605d3cb2c7da8c984341d31c2907d41 100644
--- a/src/common_types_io.cpp
+++ b/src/common_types_io.cpp
@@ -15,6 +15,8 @@ std::ostream& operator<<(std::ostream& o, arb::cell_kind k) {
return o << "cable1d_neuron";
case arb::cell_kind::data_spike_source:
return o << "data_spike_source";
+ case arb::cell_kind::lif_neuron:
+ return o << "lif_neuron";
}
return o;
}
diff --git a/src/lif_cell_description.hpp b/src/lif_cell_description.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..24c0973feb7a8b8e15ad753d2144f11511f863d2
--- /dev/null
+++ b/src/lif_cell_description.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+// Model parameteres of leaky integrate and fire neuron model.
+struct lif_cell_description {
+ // Neuronal parameters.
+ double tau_m = 10; // Membrane potential decaying constant [ms].
+ double V_th = 10; // Firing threshold [mV].
+ double C_m = 20; // Membrane capacitance [pF].
+ double E_L = 0; // Resting potential [mV].
+ double V_m = E_L; // Initial value of the Membrane potential [mV].
+ double V_reset = E_L; // Reset potential [mV].
+ double t_ref = 2; // Refractory period [ms].
+};
diff --git a/src/lif_cell_group.cpp b/src/lif_cell_group.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e8b18852547924cf5671e1d3c5b1f9cf8cdea9a9
--- /dev/null
+++ b/src/lif_cell_group.cpp
@@ -0,0 +1,104 @@
+#include <lif_cell_group.hpp>
+
+using namespace arb;
+
+// Constructor containing gid of first cell in a group and a container of all cells.
+lif_cell_group::lif_cell_group(std::vector<cell_gid_type> gids, const recipe& rec):
+gids_(std::move(gids))
+{
+ // Default to no binning of events
+ set_binning_policy(binning_kind::none, 0);
+
+ cells_.reserve(gids_.size());
+ last_time_updated_.resize(gids_.size());
+
+ for (auto lid: util::make_span(0, gids_.size())) {
+ cells_.push_back(util::any_cast<lif_cell_description>(rec.get_cell_description(gids_[lid])));
+ }
+}
+
+cell_kind lif_cell_group::get_cell_kind() const {
+ return cell_kind::lif_neuron;
+}
+
+void lif_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& event_lanes) {
+ PE("lif");
+ if (event_lanes.size() > 0) {
+ for (auto lid: util::make_span(0, gids_.size())) {
+ // Advance each cell independently.
+ advance_cell(ep.tfinal, dt, lid, event_lanes[lid]);
+ }
+ }
+ PL();
+}
+
+const std::vector<spike>& lif_cell_group::spikes() const {
+ return spikes_;
+}
+
+void lif_cell_group::clear_spikes() {
+ spikes_.clear();
+}
+
+// TODO: implement sampler
+void lif_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probe_ids,
+ schedule sched, sampler_function fn, sampling_policy policy) {}
+void lif_cell_group::remove_sampler(sampler_association_handle h) {}
+void lif_cell_group::remove_all_samplers() {}
+
+// TODO: implement binner_
+void lif_cell_group::set_binning_policy(binning_kind policy, time_type bin_interval) {
+}
+
+void lif_cell_group::reset() {
+ spikes_.clear();
+ last_time_updated_.clear();
+}
+
+// Advances a single cell (lid) with the exact solution (jumps can be arbitrary).
+// Parameter dt is ignored, since we make jumps between two consecutive spikes.
+void lif_cell_group::advance_cell(time_type tfinal, time_type dt, cell_gid_type lid, pse_vector& event_lane) {
+ // Current time of last update.
+ auto t = last_time_updated_[lid];
+ auto& cell = cells_[lid];
+ const auto n_events = event_lane.size();
+
+ // Integrate until tfinal using the exact solution of membrane voltage differential equation.
+ for (unsigned i=0; i<n_events; ++i ) {
+ auto& ev = event_lane[i];
+ const auto time = ev.time;
+ auto weight = ev.weight;
+
+ if (time < t) continue; // skip event if a neuron is in refactory period
+ if (time >= tfinal) break; // end of integration interval
+
+ // if there are events that happened at the same time as this event, process them as well
+ while (i + 1 < n_events && event_lane[i+1].time <= time) {
+ weight += event_lane[i+1].weight;
+ i++;
+ }
+
+ // Let the membrane potential decay.
+ auto decay = exp(-(time - t) / cell.tau_m);
+ cell.V_m *= decay;
+ auto update = weight / cell.C_m;
+ // Add jump due to spike.
+ cell.V_m += update;
+ t = time;
+ // If crossing threshold occurred
+ if (cell.V_m >= cell.V_th) {
+ cell_member_type spike_neuron_gid = {gids_[lid], 0};
+ spike s = {spike_neuron_gid, t};
+ spikes_.push_back(s);
+
+ // Advance the last_time_updated to account for the refractory period.
+ t += cell.t_ref;
+
+ // Reset the voltage to resting potential.
+ cell.V_m = cell.E_L;
+ }
+ }
+
+ // This is the last time a cell was updated.
+ last_time_updated_[lid] = t;
+}
diff --git a/src/lif_cell_group.hpp b/src/lif_cell_group.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a1067381f9b940eb794343e794e349b9c82a652
--- /dev/null
+++ b/src/lif_cell_group.hpp
@@ -0,0 +1,53 @@
+#pragma once
+#include <algorithm>
+#include <threading/timer.hpp>
+#include <cell_group.hpp>
+#include <event_queue.hpp>
+#include <lif_cell_description.hpp>
+#include <profiling/profiler.hpp>
+#include <recipe.hpp>
+#include <util/unique_any.hpp>
+#include <vector>
+
+namespace arb {
+class lif_cell_group: public cell_group {
+public:
+ using value_type = double;
+
+ lif_cell_group() = default;
+
+ // Constructor containing gid of first cell in a group and a container of all cells.
+ lif_cell_group(std::vector<cell_gid_type> gids, const recipe& rec);
+
+ virtual cell_kind get_cell_kind() const override;
+ virtual void reset() override;
+ virtual void set_binning_policy(binning_kind policy, time_type bin_interval) override;
+ virtual void advance(epoch epoch, time_type dt, const event_lane_subrange& events) override;
+
+ virtual const std::vector<spike>& spikes() const override;
+ virtual void clear_spikes() override;
+
+ // Sampler association methods below should be thread-safe, as they might be invoked
+ // from a sampler call back called from a different cell group running on a different thread.
+ virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function, sampling_policy) override;
+ virtual void remove_sampler(sampler_association_handle) override;
+ virtual void remove_all_samplers() override;
+
+private:
+ // Advances a single cell (lid) with the exact solution (jumps can be arbitrary).
+ // Parameter dt is ignored, since we make jumps between two consecutive spikes.
+ void advance_cell(time_type tfinal, time_type dt, cell_gid_type lid, pse_vector& event_lane);
+
+ // List of the gids of the cells in the group.
+ std::vector<cell_gid_type> gids_;
+
+ // Cells that belong to this group.
+ std::vector<lif_cell_description> cells_;
+
+ // Spikes that are generated (not necessarily sorted).
+ std::vector<spike> spikes_;
+
+ // Time when the cell was last updated.
+ std::vector<time_type> last_time_updated_;
+};
+} // namespace arb
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index d913b07b58700a540f2db93d95f72300caa2917d..f2efc57058477160ec83cce1d79a9aa1c54c8a94 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -49,6 +49,7 @@ set(TEST_SOURCES
test_event_queue.cpp
test_filter.cpp
test_fvm_multi.cpp
+ test_lif_cell_group.cpp
test_mc_cell_group.cpp
test_lexcmp.cpp
test_mask_stream.cpp
diff --git a/tests/unit/test_lif_cell_group.cpp b/tests/unit/test_lif_cell_group.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8cc7f921a22892defacbbf555b34c29eaadf7401
--- /dev/null
+++ b/tests/unit/test_lif_cell_group.cpp
@@ -0,0 +1,233 @@
+#include "../gtest.h"
+#include <cell_group_factory.hpp>
+#include <fstream>
+#include <lif_cell_description.hpp>
+#include <lif_cell_group.hpp>
+#include <load_balance.hpp>
+#include <model.hpp>
+#include <recipe.hpp>
+#include <rss_cell.hpp>
+#include <rss_cell_group.hpp>
+
+using namespace arb;
+// Simple ring network of LIF neurons.
+// with one regularly spiking cell (fake cell) connected to the first cell in the ring.
+class ring_recipe: public arb::recipe {
+public:
+ ring_recipe(cell_size_type n_lif_cells, float weight, float delay):
+ n_lif_cells_(n_lif_cells), weight_(weight), delay_(delay)
+ {}
+
+ cell_size_type num_cells() const override {
+ return n_lif_cells_ + 1;
+ }
+
+ // LIF neurons have gid in range [1..n_lif_cells_] whereas fake cell is numbered with 0.
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ if (gid == 0) {
+ return cell_kind::regular_spike_source;
+ }
+ return cell_kind::lif_neuron;
+ }
+
+ std::vector<cell_connection> connections_on(cell_gid_type gid) const override {
+ if (gid == 0) {
+ return {};
+ }
+
+ // In a ring, each cell has just one incoming connection.
+ std::vector<cell_connection> connections;
+ // gid-1 >= 0 since gid != 0
+ cell_member_type source{(gid - 1) % n_lif_cells_, 0};
+ cell_member_type target{gid, 0};
+ cell_connection conn(source, target, weight_, delay_);
+ connections.push_back(conn);
+
+ // If first LIF cell, then add
+ // the connection from the last LIF cell as well
+ if (gid == 1) {
+ cell_member_type source{n_lif_cells_, 0};
+ cell_member_type target{gid, 0};
+ cell_connection conn(source, target, weight_, delay_);
+ connections.push_back(conn);
+ }
+
+ return connections;
+ }
+
+ util::unique_any get_cell_description(cell_gid_type gid) const override {
+ // regularly spiking cell.
+ if (gid == 0) {
+ // Produces just a single spike at time 0ms.
+ auto rs = arb::rss_cell();
+ rs.start_time = 0;
+ rs.period = 1;
+ rs.stop_time = 0.5;
+ return rs;
+ }
+ // LIF cell.
+ return lif_cell_description();
+ }
+
+ cell_size_type num_sources(cell_gid_type) const override {
+ return 1;
+ }
+ cell_size_type num_targets(cell_gid_type) const override {
+ return 1;
+ }
+ cell_size_type num_probes(cell_gid_type) const override {
+ return 0;
+ }
+ probe_info get_probe(cell_member_type probe_id) const override {
+ return {};
+ }
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
+private:
+ cell_size_type n_lif_cells_;
+ float weight_, delay_;
+};
+
+// LIF cells connected in the manner of a path 0->1->...->n-1.
+class path_recipe: public arb::recipe {
+public:
+ path_recipe(cell_size_type n, float weight, float delay):
+ ncells_(n), weight_(weight), delay_(delay)
+ {}
+
+ cell_size_type num_cells() const override {
+ return ncells_;
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::lif_neuron;
+ }
+
+ std::vector<cell_connection> connections_on(cell_gid_type gid) const override {
+ if (gid == 0) {
+ return {};
+ }
+ std::vector<cell_connection> connections;
+ cell_member_type source{gid - 1, 0};
+ cell_member_type target{gid, 0};
+ cell_connection conn(source, target, weight_, delay_);
+ connections.push_back(conn);
+
+ return connections;
+ }
+
+ util::unique_any get_cell_description(cell_gid_type gid) const override {
+ return lif_cell_description();
+ }
+
+ cell_size_type num_sources(cell_gid_type) const override {
+ return 1;
+ }
+ cell_size_type num_targets(cell_gid_type) const override {
+ return 1;
+ }
+ cell_size_type num_probes(cell_gid_type) const override {
+ return 0;
+ }
+ probe_info get_probe(cell_member_type probe_id) const override {
+ return {};
+ }
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
+private:
+ cell_size_type ncells_;
+ float weight_, delay_;
+};
+
+TEST(lif_cell_group, recipe)
+{
+ ring_recipe rr(100, 1, 0.1);
+ EXPECT_EQ(101u, rr.num_cells());
+ EXPECT_EQ(2u, rr.connections_on(1u).size());
+ EXPECT_EQ(1u, rr.connections_on(55u).size());
+ EXPECT_EQ(0u, rr.connections_on(1u)[0].source.gid);
+ EXPECT_EQ(100u, rr.connections_on(1u)[1].source.gid);
+}
+
+TEST(lif_cell_group, spikes) {
+ // make two lif cells
+ path_recipe recipe(2, 1000, 0.1);
+
+ hw::node_info nd;
+ nd.num_cpu_cores = threading::num_threads();
+
+ auto decomp = partition_load_balance(recipe, nd);
+ model m(recipe, decomp);
+
+ std::vector<postsynaptic_spike_event> events;
+
+ // First event to trigger the spike (first neuron).
+ events.push_back({{0, 0}, 1, 1000});
+
+ // This event happens inside the refractory period of the previous
+ // event, thus, should be ignored (first neuron)
+ events.push_back({{0, 0}, 1.1, 1000});
+
+ // This event happens long after the refractory period of the previous
+ // event, should thus trigger new spike (first neuron).
+ events.push_back({{0, 0}, 50, 1000});
+
+ m.inject_events(events);
+
+ time_type tfinal = 100;
+ time_type dt = 0.01;
+ m.run(tfinal, dt);
+
+ // we expect 4 spikes: 2 by both neurons
+ EXPECT_EQ(4u, m.num_spikes());
+}
+
+TEST(lif_cell_group, ring)
+{
+ // Total number of LIF cells.
+ cell_size_type num_lif_cells = 99;
+ double weight = 1000;
+ double delay = 1;
+
+ hw::node_info nd;
+ nd.num_cpu_cores = threading::num_threads();
+
+ // Total simulation time.
+ time_type simulation_time = 100;
+
+ auto recipe = ring_recipe(num_lif_cells, weight, delay);
+ auto decomp = partition_load_balance(recipe, nd);
+
+ // Creates a model with a ring recipe of lif neurons
+ model mod(recipe, decomp);
+
+ std::vector<spike> spike_buffer;
+
+ mod.set_global_spike_callback(
+ [&spike_buffer](const std::vector<spike>& spikes) {
+ spike_buffer.insert(spike_buffer.end(), spikes.begin(), spikes.end());
+ }
+ );
+
+ // Runs the simulation for simulation_time with given timestep
+ mod.run(simulation_time, 0.01);
+ // The total number of cells in all the cell groups.
+ // There is one additional fake cell (regularly spiking cell).
+ EXPECT_EQ(num_lif_cells + 1u, recipe.num_cells());
+
+ for (auto& spike : spike_buffer) {
+ // Assumes that delay = 1
+ // We expect that Regular Spiking Cell spiked at time 0s.
+ if (spike.source.gid == 0) {
+ EXPECT_EQ(0, spike.time);
+ // Other LIF cell should spike consecutively.
+ } else {
+ EXPECT_EQ(spike.source.gid, spike.time);
+ }
+ }
+}
+