LIF neurons: CPU backend with Brunel Network miniapp (#441)

Two main contributions:

1) Implementation of LIF neuron model with no kernel and no external input (I_e=0)

The input current to each neuron is therefore just the sum of all the weights of incoming spikes. We integrate in jumps dt = min(t_final - t, t_event - t), since we know the exact solution of the differential equal describing the membrane potential.

2) 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. 

.gitignore

@@ -66,3 +66,6 @@ commit.msg
 
 # eclipse remote sync folders
 .ptp-sync-folder
+
+# Mac default files
+.DS_Store

example/CMakeLists.txt

 add_subdirectory(miniapp)
 add_subdirectory(generators)
+add_subdirectory(brunel)

example/brunel/CMakeLists.txt

+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"
+)

example/brunel/brunel_miniapp.cpp

+#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";
+}
+

example/brunel/io.cpp

+#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

example/brunel/io.hpp

+#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

example/brunel/partitioner.hpp

+#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;
+    }
+}

example/brunel/readme.md

+## 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
+```

example/miniapp/CMakeLists.txt

@@ -21,3 +21,4 @@ set_target_properties(
     PROPERTIES
     RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/example"
 )
+

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

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);
 

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
 };

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;
 }

src/lif_cell_description.hpp

+#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].
+};

src/lif_cell_group.cpp

+#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;
+}

src/lif_cell_group.hpp

+#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

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

tests/unit/test_lif_cell_group.cpp

+#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);
+        }
+    }
+}
+