diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 7e0a828d7895e79fafa34a82285f0c023cc8d2ef..4730d3d18565649a951cad791dc007cf76f10f5f 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(miniapp) +add_subdirectory(generators) diff --git a/example/generators/CMakeLists.txt b/example/generators/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..1e743f11c286f926f294fd0fa68bc48cca67e91e --- /dev/null +++ b/example/generators/CMakeLists.txt @@ -0,0 +1,15 @@ +add_executable(event_gen.exe event_gen.cpp) + +target_link_libraries(event_gen.exe LINK_PUBLIC ${ARB_LIBRARIES}) +target_link_libraries(event_gen.exe LINK_PUBLIC ${EXTERNAL_LIBRARIES}) + +if(ARB_WITH_MPI) + target_link_libraries(event_gen.exe LINK_PUBLIC ${MPI_C_LIBRARIES}) + set_property(TARGET event_gen.exe APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}") +endif() + +set_target_properties( + event_gen.exe + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/example" +) diff --git a/example/generators/event_gen.cpp b/example/generators/event_gen.cpp new file mode 100644 index 0000000000000000000000000000000000000000..df69cd2e8599d735fe2d94b32fcd05d5dac27a9f --- /dev/null +++ b/example/generators/event_gen.cpp @@ -0,0 +1,178 @@ +/* + * A miniapp that demonstrates how to use event generators. + * + * The miniapp builds a simple model of a single cell, with one compartment + * corresponding to the soma. The soma has a single synapse, to which two + * event generators, one inhibitory, and one excitatory, are attached. + */ + +#include <fstream> +#include <iomanip> +#include <iostream> + +#include <json/json.hpp> + +#include <cell.hpp> +#include <common_types.hpp> +#include <event_generator.hpp> +#include <hardware/node_info.hpp> +#include <load_balance.hpp> +#include <model.hpp> +#include <recipe.hpp> +#include <simple_sampler.hpp> + +using arb::cell_gid_type; +using arb::cell_lid_type; +using arb::cell_size_type; +using arb::cell_member_type; +using arb::cell_kind; +using arb::time_type; +using arb::cell_probe_address; + +// Writes voltage trace as a json file. +void write_trace_json(const arb::trace_data<double>& trace); + +class generator_recipe: public arb::recipe { +public: + // There is just the one cell in the model + cell_size_type num_cells() const override { + return 1; + } + + // Create cell with just a single compartment for the soma: + // soma diameter: 18.8 µm + // mechanisms: pas [default params] + // bulk resistivitiy: 100 Ω·cm [default] + // capacitance: 0.01 F/m² [default] + // synapses: 1 * expsyn + arb::util::unique_any get_cell_description(cell_gid_type gid) const override { + arb::cell c; + + c.add_soma(18.8/2.0); // convert 18.8 μm diameter to radius + c.soma()->add_mechanism("pas"); + + // Add one synapse at the soma. + // This synapse will be the target for all events, from both + // event_generators. + auto syn_spec = arb::mechanism_spec("expsyn"); + c.add_synapse({0, 0.5}, syn_spec); + + return std::move(c); + } + + cell_kind get_cell_kind(cell_gid_type gid) const override { + EXPECTS(gid==0); // There is only one cell in the model + return cell_kind::cable1d_neuron; + } + + // The cell has one target synapse, which receives both inhibitory and exchitatory inputs. + cell_size_type num_targets(cell_gid_type gid) const override { + EXPECTS(gid==0); // There is only one cell in the model + return 1; + } + + // Return two generators attached to the one cell. + std::vector<arb::event_generator_ptr> event_generators(cell_gid_type gid) const override { + EXPECTS(gid==0); // There is only one cell in the model + + using RNG = std::mt19937_64; + using pgen = arb::poisson_generator<RNG>; + + auto hz_to_freq = [](double hz) { return hz*1e-3; }; + time_type t0 = 0; + + // Define frequencies and weights for the excitatory and inhibitory generators. + double lambda_e = hz_to_freq(500); + double lambda_i = hz_to_freq(20); + double w_e = 0.001; + double w_i = -0.005; + + // Make two event generators. + std::vector<arb::event_generator_ptr> gens; + + // Add excitatory generator + gens.push_back( + arb::make_event_generator<pgen>( + cell_member_type{0,0}, // Target synapse (gid, local_id). + w_e, // Weight of events to deliver + RNG(29562872), // Random number generator to use + t0, // Events start being delivered from this time + lambda_e)); // Expected frequency (events per ms) + + // Add inhibitory generator + gens.push_back( + arb::make_event_generator<pgen>( + cell_member_type{0,0}, w_i, RNG(86543891), t0, lambda_i)); + + return gens; + } + + // There is one probe (for measuring voltage at the soma) on the cell + cell_size_type num_probes(cell_gid_type gid) const override { + EXPECTS(gid==0); // There is only one cell in the model + return 1; + } + + arb::probe_info get_probe(cell_member_type id) const override { + EXPECTS(id.gid==0); // There is one cell, + EXPECTS(id.index==0); // with one probe. + + // Get the appropriate kind for measuring voltage + cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage; + // Measure at the soma + arb::segment_location loc(0, 0.0); + + return arb::probe_info{id, kind, cell_probe_address{loc, kind}}; + } +}; + +int main() { + // Create an instance of our recipe. + generator_recipe recipe; + + // Make the domain decomposition for the model + auto node = arb::hw::get_node_info(); + auto decomp = arb::partition_load_balance(recipe, node); + + // Construct the model. + arb::model model(recipe, decomp); + + // Set up the probe that will measure voltage in the cell. + + // The id of the only probe on the cell: the cell_member type points to (cell 0, probe 0) + auto probe_id = cell_member_type{0, 0}; + // The schedule for sampling is 10 samples every 1 ms. + auto sched = arb::regular_schedule(0.1); + // This is where the voltage samples will be stored as (time, value) pairs + arb::trace_data<double> voltage; + // Now attach the sampler at probe_id, with sampling schedule sched, writing to voltage + model.add_sampler(arb::one_probe(probe_id), sched, arb::make_simple_sampler(voltage)); + + // Run the model for 1 s (1000 ms), with time steps of 0.01 ms. + model.run(50, 0.01); + + // Write the samples to a json file. + write_trace_json(voltage); +} + +void write_trace_json(const arb::trace_data<double>& trace) { + std::string path = "./voltages.json"; + + nlohmann::json json; + json["name"] = "event_gen_demo"; + json["units"] = "mV"; + json["cell"] = "0.0"; + json["probe"] = "0"; + + auto& jt = json["data"]["time"]; + auto& jy = json["data"]["voltage"]; + + for (const auto& sample: trace) { + jt.push_back(sample.t); + jy.push_back(sample.v); + } + + std::ofstream file(path); + file << std::setw(1) << json << "\n"; +} + diff --git a/example/generators/readme.md b/example/generators/readme.md new file mode 100644 index 0000000000000000000000000000000000000000..d910aadcea678e192046a49acd6a0f61a2f1c299 --- /dev/null +++ b/example/generators/readme.md @@ -0,0 +1,232 @@ +# Event Generator Example + +A miniapp that demonstrates how to describe event generators in a `arb::recipe`. + +This miniapp uses a simple model of a single cell, with one compartment corresponding to the soma. +The soma has one synapse, to which two event generators, one inhibitory and one excitatory, are attached. + +The following topics are covered: +* How to describe a simple single cell model with a `recipe`. +* How to connect Poisson event generators to a synapse in a `recipe`. +* How to sample the voltage on the soma of a cell. + +## The Recipe + +We need to create a recipe that describes the model. +All models derive from the `arb::recipe` base class. + +We call the recipe for this example `generator_recipe`, and declare it as follows: + +```C++ +class generator_recipe: public arb::recipe { +public: + // There is just the one cell in the model + cell_size_type num_cells() const override { + return 1; + } + + // ... +}; +``` + +### Describing The Network + +Every user-defined recipe must provide implementations for the following three methods: +* `recipe::num_cells()`: return the total number of cells in the model. +* `recipe::get_cell_description(gid)`: return a description of the cell with `gid`. +* `recipe::get_cell_kind(gid)`: return an `arb::cell_kind` enum value for the cell type. + +In addition to these three, we also need to implement +`recipe::num_targets(gid)` to return 1, to indicate that there is one target +(i.e. synapse) on the cell. + +This single cell model has no connections for spike communication, so the +default `recipe::connections_on(gid)` method can use the default implementation, +which returns an empty list of connections for a cell. + +Above you can see the definition of `recipe::num_cells()`, which is trivial for this model, which has only once cell. +The `recipe::get_cell_description(gid)` and `recipe::get_cell_kind(gid)` methods are defined to return a single compartment cell with one synapse, and a `arb::cell_kind::cable1d_neuron` respectively: + +```C++ + // Return an arb::cell that describes a single compartment cell with + // passive dynamics, and a single expsyn synapse. + arb::util::unique_any get_cell_description(cell_gid_type gid) const override { + arb::cell c; + + c.add_soma(18.8/2.0); + c.soma()->add_mechanism("pas"); + + // Add one synapse at the soma. + // This synapse will be the target for all events, from both event_generators. + auto syn_spec = arb::mechanism_spec("expsyn"); + c.add_synapse({0, 0.5}, syn_spec); + + return std::move(c); // move into unique_any wrapper + } + + cell_kind get_cell_kind(cell_gid_type gid) const override { + return cell_kind::cable1d_neuron; + } + + // There is one synapse, i.e. target, on the cell. + cell_size_type num_targets(cell_gid_type gid) const override { + EXPECTS(gid==0); // There is only one cell in the model + return 1; + } +``` + +### Event Generators + +For our demo, we want to attach two Poisson event generators to the synapse on our cell. + +1. An excitatory synapse with spiking frequency λ\_e and weight w\_e. +2. An inhibitory synapse with spiking frequency λ\_i and weight w\_i. + +To add the event generators to the synapse, we implement the `recipe::event_generators(gid)` method. + +The implementation of this with hard-coded frequencies and weights is: + +```C++ + std::vector<arb::event_generator_ptr> event_generators(cell_gid_type gid) const override { + // The type of random number generator to use. + using RNG = std::mt19937_64; + + // The poisson_generator is templated on the random number generator type. + using pgen = arb::poisson_generator<RNG>; + + auto hz_to_freq = [](double hz) { return hz*1e-3; }; + time_type t0 = 0; + + // Define frequencies and weights for the excitatory and inhibitory generators. + double lambda_e = hz_to_freq(500); + double lambda_i = hz_to_freq(20); + double w_e = 0.001; + double w_i = -0.005; + + // Make two event generators. + std::vector<arb::event_generator_ptr> gens; + + // Add excitatory generator + gens.push_back( + arb::make_event_generator<pgen>( + cell_member_type{0,0}, // Target synapse (gid, local_id). + w_e, // Weight of events to deliver + RNG(29562872), // Random number generator to use + t0, // Events start being delivered from this time + lambda_e)); // Expected frequency (events per ms) + + // Add inhibitory generator + gens.push_back( + arb::make_event_generator<pgen>( + cell_member_type{0,0}, w_i, RNG(86543891), t0, lambda_i)); + + return gens; + } +``` + +The `recipe::event_generators(gid)` method returns a vector of `event_generator_ptr` that are attached to the cell with `gid`. + +The `event_generator_ptr` is an alias for a `unique_ptr` wrapped around an `event_generator`. + +```C++ +using event_generator_ptr = std::unique_ptr<event_generator>; +``` + +In the implementation, an empty vector is created, and the generators are created and `push_back`ed into the vector one after the other. + +The helper function `make_event_generator` is used to simplify the process of creating an event generator and wrapping it in a `unique_ptr`. +It has the same semantics as [http://en.cppreference.com/w/cpp/memory/unique_ptr/make_unique](`std::make_unique`): + +* It takes as a template parameter the specialized type of `event_generator`, in this case `pgen`. +* Then it takes as arguments the arguments to pass to the constructor of `pgen`. + +Of the arguments passed to the Poisson event generators, the random number generator state require further explanation. +Each Poisson generator has its own private random number generator state. +The initial random number state is provided on construction. +For a real world model, the state should have a seed that is some reproducable hash of `gid` and the generator id, to ensure reproducable random streams. +For this simple example, we use hard coded seeds to initialize the random number state. + +### Sampling Voltages + +To visualise the result of our experiment, we want to sample the voltage in our cell at regular time points and save the resulting sequence to a JSON file. + +There are three parts to this process. + +1. Define the a `probe` in the recipe, which describes the location and variable to be sampled, i.e. the soma and voltage respectively for this example. +2. Attach a sampler to this probe once the model has been created. +3. Write the sampled values to a JSON file once the model has been run. + +#### 1. Define probe in recipe + +The `recipe::num_probes(gid)` and `recipe::get_probe(id)` have to be defined for sampling. +In our case, the cell has one probe, which refers to the voltage at the soma. + +```C++ + // There is one probe (for measuring voltage at the soma) on the cell + cell_size_type num_probes(cell_gid_type gid) const override { + return 1; + } + + arb::probe_info get_probe(cell_member_type id) const override { + // Get the appropriate kind for measuring voltage + cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage; + + // The location at which we measure: position 0 on segment 0. + // The cell has only one segment, segment 0, which is the soma. + arb::segment_location loc(0, 0.0); + + // Put this together into a `probe_info` + return arb::probe_info{id, kind, cell_probe_address{loc, kind}}; + } +``` + +#### 2. Attach sampler to model + +Once the model has been created, the model has internal state that tracks the value of the voltage at the probe location described by the recipe. +We must attach a sampler to this probe to get sample values. + +The sampling interface is rich, and can be extended in many ways. +For our simple use case there are three bits of information that need to be provided when creating a sampler + +1. The `probe_id` that identifies the probe (generated in the recipe). +2. The schedule to use for sampling (in our case 10 samples every ms). +3. The location where we want to save the samples for outputing later. + +```C++ + // The id of the only probe (cell 0, probe 0) + auto probe_id = cell_member_type{0, 0}; + // The schedule for sampling is 10 samples every 1 ms. + auto sched = arb::regular_schedule(0.1); + // Where the voltage samples will be stored as (time, value) pairs + arb::trace_data<double> voltage; + // Now attach the sampler: + model.add_sampler(arb::one_probe(probe_id), sched, arb::make_simple_sampler(voltage)); +``` + +When the model is run, the `simple_sampler` that we attached to the probe will store the sample values as (time, voltage) pairs in the `voltage` vector. + +#### 3. Output to JSON + +The function `write_trace_json()` uses the [JSON library](https://github.com/nlohmann/json) that is included in the Arbor library to write the output. +We don't document its implementation here, instead we note that the format of the output was chosen to be compatible with the `tsplot` script provided by Arbor for plotting traces. + +## Visualising the Results. + +The voltage samples are saved to `voltages.json`, and can be visualised using the `tsplot` Python script distributed with Arbor. +Here is an example set of steps used to build, run and plot the voltage trace: + +```bash +# assuming that build is a subdirectory of the main arbor project path +cd build + +# build the event generator demo +make -j event_gen.exe + +# run the event generator demo +./example/event_gen.exe + +# draw a plot of the results. +# uses the matplotlib in Python. +../scripts/tsplot voltages.json +``` + diff --git a/src/model.cpp b/src/model.cpp index 3afb9f4df9a40690d7d7db4ac6faea7494769ab5..36e8c8e1c033ec9bb6f2f69905427a797e6b2ba2 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -284,10 +284,11 @@ void model::inject_events(const pse_vector& events) { for (auto& e: events) { if (e.time<t_) { throw std::runtime_error( - "model::inject_events(): attempt to inject an event at time " + "model::inject_events(): attempt to inject an event at time: " + std::to_string(e.time) - + ", when model state is at time " - + std::to_string(t_)); + + " ms, which is earlier than the current model time: " + + std::to_string(t_) + + " ms. Events must be injected on or after the current model time."); } // local_cell_index returns an optional type that evaluates // to true iff the gid is a local cell.