Skip to content
Snippets Groups Projects
Commit df3bc45d authored by Sam Yates's avatar Sam Yates Committed by Benjamin Cumming
Browse files

Add single cell example. (#703)

Adds a small self contained example demonstrating the use of SWC and morphology specifications, to fill a void that will be left when the miniapp is removed.

* Adds `single-cell` example, together with short README and example SWC file.
* Removes redundant constructors for context objects.
* Corrects and simplifies some of the context comments.
* Add a method to `mc_segment` for calculating (approximately) a lower bound on the length constant of that segment for a given frequency, for use with a NEURON 'd-lambda' style discretization rule.
parent 187d6a24
No related branches found
No related tags found
No related merge requests found
...@@ -14,10 +14,6 @@ ...@@ -14,10 +14,6 @@
namespace arb { namespace arb {
execution_context::execution_context():
execution_context(proc_allocation())
{}
execution_context::execution_context(const proc_allocation& resources): execution_context::execution_context(const proc_allocation& resources):
distributed(make_local_context()), distributed(make_local_context()),
thread_pool(std::make_shared<threading::task_system>(resources.num_threads)), thread_pool(std::make_shared<threading::task_system>(resources.num_threads)),
...@@ -25,10 +21,6 @@ execution_context::execution_context(const proc_allocation& resources): ...@@ -25,10 +21,6 @@ execution_context::execution_context(const proc_allocation& resources):
: std::make_shared<gpu_context>()) : std::make_shared<gpu_context>())
{} {}
context make_context() {
return context(new execution_context(), [](execution_context* p){delete p;});
}
context make_context(const proc_allocation& p) { context make_context(const proc_allocation& p) {
return context(new execution_context(p), [](execution_context* p){delete p;}); return context(new execution_context(p), [](execution_context* p){delete p;});
} }
......
...@@ -24,8 +24,7 @@ struct execution_context { ...@@ -24,8 +24,7 @@ struct execution_context {
task_system_handle thread_pool; task_system_handle thread_pool;
gpu_context_handle gpu; gpu_context_handle gpu;
execution_context(); execution_context(const proc_allocation& resources = proc_allocation{});
execution_context(const proc_allocation& resources);
// Use a template for constructing with a specific distributed context. // Use a template for constructing with a specific distributed context.
// Specialised implementations are implemented in execution_context.cpp. // Specialised implementations are implemented in execution_context.cpp.
......
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
namespace arb { namespace arb {
/// Requested dry-run parameters // Requested dry-run parameters.
struct dry_run_info { struct dry_run_info {
unsigned num_ranks; unsigned num_ranks;
unsigned num_cells_per_rank; unsigned num_cells_per_rank;
...@@ -13,17 +13,18 @@ struct dry_run_info { ...@@ -13,17 +13,18 @@ struct dry_run_info {
num_cells_per_rank(cells_per_rank) {} num_cells_per_rank(cells_per_rank) {}
}; };
/// A subset of local computation resources to use in a computation. // A description of local computation resources to use in a computation.
// By default, a proc_allocation will comprise one thread and no GPU.
struct proc_allocation { struct proc_allocation {
unsigned num_threads; unsigned num_threads;
// The gpu id corresponds to the `int device` parameter used by CUDA API calls // The gpu id corresponds to the `int device` parameter used by
// to identify gpu devices. // CUDA API calls to identify gpu devices.
// Set to -1 to indicate that no GPU device is to be used. // A gpud id of -1 indicates no GPU device is to be used.
// see CUDA documenation for cudaSetDevice and cudaDeviceGetAttribute // See CUDA documenation for cudaSetDevice and cudaDeviceGetAttribute.
int gpu_id; int gpu_id;
// By default use one thread and no GPU.
proc_allocation(): proc_allocation(1, -1) {} proc_allocation(): proc_allocation(1, -1) {}
proc_allocation(unsigned threads, int gpu): proc_allocation(unsigned threads, int gpu):
...@@ -36,30 +37,24 @@ struct proc_allocation { ...@@ -36,30 +37,24 @@ struct proc_allocation {
} }
}; };
// arb::execution_context is a container defined in the implementation for state // arb::execution_context encapsulates the execution resources used in
// related to execution resources, specifically thread pools, gpus and MPI // a simulation, namely the task system thread pools, GPU handle, and
// communicators. // MPI communicator if applicable.
// Forward declare execution_context. // Forward declare execution_context.
struct execution_context; struct execution_context;
// arb::context is an opaque handle for this container presented in the // arb::context is an opaque handle for the execution context for use
// public API. // in the public API, implemented as a unique pointer.
// It doesn't make sense to copy contexts, so we use a std::unique_ptr to
// implement the handle with lifetime management.
// //
// Because execution_context is an incomplete type, a destructor prototype must // As execution_context is an incomplete type, an explicit deleter must be
// be provided. // provided.
using context = std::unique_ptr<execution_context, void (*)(execution_context*)>; using context = std::unique_ptr<execution_context, void (*)(execution_context*)>;
// Helpers for creating contexts. These are implemented in the back end. // Helpers for creating contexts. These are implemented in the back end.
// Non-distributed context that uses all detected threads and one GPU if available. // Non-distributed context using the requested resources.
context make_context(); context make_context(const proc_allocation& resources = proc_allocation{});
// Non-distributed context that uses resources described by resources
context make_context(const proc_allocation& resources);
// Distributed context that uses MPI communicator comm, and local resources // Distributed context that uses MPI communicator comm, and local resources
// described by resources. Or dry run context that uses dry_run_info. // described by resources. Or dry run context that uses dry_run_info.
......
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#include <arbor/assert.hpp> #include <arbor/assert.hpp>
#include <arbor/common_types.hpp> #include <arbor/common_types.hpp>
#include <arbor/math.hpp>
#include <arbor/morphology.hpp> #include <arbor/morphology.hpp>
#include <arbor/mechinfo.hpp> #include <arbor/mechinfo.hpp>
#include <arbor/point.hpp> #include <arbor/point.hpp>
...@@ -134,6 +135,11 @@ public: ...@@ -134,6 +135,11 @@ public:
return false; return false;
} }
// Approximate frequency-dependent length constant lower bound.
virtual value_type length_constant(value_type freq_Hz) const {
return 0;
}
util::optional<mechanism_desc&> mechanism(const std::string& name) { util::optional<mechanism_desc&> mechanism(const std::string& name) {
auto it = std::find_if(mechanisms_.begin(), mechanisms_.end(), auto it = std::find_if(mechanisms_.begin(), mechanisms_.end(),
[&](mechanism_desc& m) { return m.name()==name; }); [&](mechanism_desc& m) { return m.name()==name; });
...@@ -291,6 +297,25 @@ public: ...@@ -291,6 +297,25 @@ public:
return std::accumulate(lengths_.begin(), lengths_.end(), value_type{}); return std::accumulate(lengths_.begin(), lengths_.end(), value_type{});
} }
value_type length_constant(value_type freq_Hz) const override {
// Following Hine and Carnevale (2001), "NEURON: A Tool for Neuroscientists",
// Neuroscientist 7, pp. 123-135.
//
// λ(f) = approx. sqrt(diameter/(pi*f*rL*cm))/2.
//
// Pick smallest non-zero diameter in the segment.
value_type r_min = 0;
for (auto r: radii_) {
if (r>0 && (r_min==0 || r<r_min)) r_min = r;
}
value_type d_min = r_min*2e-6; // [m]
value_type rc = rL*0.01*cm; // [s/m]
value_type lambda = std::sqrt(d_min/(math::pi<double>*freq_Hz*rc))/2.; // [m]
return lambda*1e6; // [µm]
}
bool has_locations() const bool has_locations() const
{ {
return locations_.size() > 0; return locations_.size() > 0;
......
...@@ -9,3 +9,4 @@ add_subdirectory(brunel) ...@@ -9,3 +9,4 @@ add_subdirectory(brunel)
add_subdirectory(bench) add_subdirectory(bench)
add_subdirectory(ring) add_subdirectory(ring)
add_subdirectory(gap_junctions) add_subdirectory(gap_junctions)
add_subdirectory(single)
add_executable(single-cell EXCLUDE_FROM_ALL single.cpp)
add_dependencies(examples single-cell)
target_link_libraries(single-cell PRIVATE arbor arborenv arbor-sup)
# 'single' example.
Example of simulating a single neuron with morphology described by an SWC file.
A cell is constructed from a supplied morphology with H–H channels
on the soma and passive channels on the dendrites. A simple exponential
synapse is added at the end of the last dendrite in the morphology,
and is triggered at time t = 1 ms.
The simulation outputs a trace of the soma membrane voltage in a simple CSV
format.
## Features
The example demonstrates the use of:
* Generating a morphology from an SWC file.
* Using a morphology to construct a cable cell.
* Injecting an artificial spike event into the simulation.
* Adding a voltage probe to a cell and running a sampler on the simulation.
## Running the example
By default, `single-cell` will simulate a 'ball-and-stick' neuron for 20 ms,
with a maxium dt of 0.025 ms and samples taken every 0.1 ms. The default
synaptic weight is 0.01 µS.
### Command line options
| Option | Effect |
|-----------------------|--------|
| -m, --morphology FILE | Load the morphology from FILE in SWC format |
| -d, --dt TIME | Set the maximum integration time step [ms] |
| -t, --t-end TIME | Set the simulation duration [ms] |
| -w, --weight WEIGHT | Set the synaptic weight [µS] |
This diff is collapsed.
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
#include <arbor/load_balance.hpp>
#include <arbor/mc_cell.hpp>
#include <arbor/morphology.hpp>
#include <arbor/swcio.hpp>
#include <arbor/simulation.hpp>
#include <arbor/simple_sampler.hpp>
#include <sup/tinyopt.hpp>
struct options {
std::string swc_file;
double t_end = 20;
double dt = 0.025;
float syn_weight = 0.01;
};
options parse_options(int argc, char** argv);
arb::morphology default_morphology();
arb::morphology read_swc(const std::string& path);
struct single_recipe: public arb::recipe {
explicit single_recipe(arb::morphology m): morpho(m) {}
arb::cell_size_type num_cells() const override { return 1; }
arb::cell_size_type num_probes(arb::cell_gid_type) const override { return 1; }
arb::cell_size_type num_targets(arb::cell_gid_type) const override { return 1; }
arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
arb::segment_location mid_soma = {0, 0.5};
arb::cell_probe_address probe = {mid_soma, arb::cell_probe_address::membrane_voltage};
// Probe info consists of: the probe id, a tag value to distinguish this probe
// from others for any attached sampler (unused), and the cell probe address.
return {probe_id, 0, probe};
}
arb::cell_kind get_cell_kind(arb::cell_gid_type) const override {
return arb::cell_kind::cable1d_neuron;
}
arb::util::unique_any get_cell_description(arb::cell_gid_type) const override {
arb::mc_cell c = make_mc_cell(morpho);
// Add HH mechanism to soma, passive channels to dendrites.
// Discretize dendrites according to the NEURON d-lambda rule.
for (auto& segment: c.segments()) {
if (segment->is_soma()) {
segment->add_mechanism("hh");
}
else {
segment->add_mechanism("pas");
double dx = segment->length_constant(100.)*0.3;
unsigned n = std::ceil(segment->as_cable()->length()/dx);
segment->set_compartments(n);
}
}
// Add synapse to last segment.
arb::cell_lid_type last_segment = morpho.components()-1;
arb::segment_location end_last_segment = { last_segment, 1. };
c.add_synapse(end_last_segment, "exp2syn");
return c;
}
arb::morphology morpho;
};
int main(int argc, char** argv) {
try {
options opt = parse_options(argc, argv);
single_recipe R(opt.swc_file.empty()? default_morphology(): read_swc(opt.swc_file));
auto context = arb::make_context();
arb::simulation sim(R, arb::partition_load_balance(R, context), context);
// Attach a sampler to the probe described in the recipe, sampling every 0.1 ms.
arb::trace_data<double> trace;
sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1), arb::make_simple_sampler(trace));
// Trigger the single synapse (target is gid 0, index 0) at t = 1 ms with
// the given weight.
arb::spike_event spike = {{0, 0}, 1., opt.syn_weight};
sim.inject_events({spike});
sim.run(opt.t_end, opt.dt);
std::cout << std::fixed << std::setprecision(4);
for (auto entry: trace) {
std::cout << entry.t << ", " << entry.v << "\n";
}
}
catch (std::exception& e) {
std::cerr << "caught exception: " << e.what() << "\n";
return 2;
}
}
options parse_options(int argc, char** argv) {
using namespace to;
options opt;
char** arg = argv+1;
while (*arg) {
if (auto dt = parse_opt<double>(arg, 'd', "dt")) {
opt.dt = dt.value();
}
else if (auto t_end = parse_opt<double>(arg, 't', "t-end")) {
opt.t_end = t_end.value();
}
else if (auto weight = parse_opt<float>(arg, 'w', "weight")) {
opt.syn_weight = weight.value();
}
else if (auto swc = parse_opt<std::string>(arg, 'm', "morphology")) {
opt.swc_file = swc.value();
}
else {
usage(argv[0], "[-m|--morphology SWCFILE] [-d|--dt TIME] [-t|--t-end TIME] [-w|--weight WEIGHT]");
std::exit(1);
}
}
return opt;
}
// If no SWC file is given, the default morphology consists
// of a soma of radius 6.3 µm and a single unbranched dendrite
// of length 200 µm and radius decreasing linearly from 0.5 µm
// to 0.2 µm.
arb::morphology default_morphology() {
arb::morphology m;
// Points in a morphology are expressed as (x, y, z, r),
// with r being the radius. Units are µm.
m.soma = {0., 0., 0., 6.3};
std::vector<arb::section_point> dendrite = {
{6.3, 0., 0., 0.5},
{206.3, 0., 0., 0.2}
};
m.add_section(dendrite, 0);
return m;
}
arb::morphology read_swc(const std::string& path) {
std::ifstream f(path);
if (!f) throw std::runtime_error("unable to open SWC file: "+path);
return arb::swc_as_morphology(arb::parse_swc_file(f));
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment