diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 4c01f7ca1fd5104f735be0433ad17193f3425e19..0e73398aeeb3ef5da3dc86ccda9a67063a7cf16a 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -27,16 +27,17 @@
# OTHER DEALINGS IN THE SOFTWARE.
#
# For more information, please refer to <http://unlicense.org/>
-
+
import os
import ycm_core
-
+
# These are the compilation flags that will be used in case there's no
# compilation database set (by default, one is not set).
# CHANGE THIS LIST OF FLAGS. YES, THIS IS THE DROID YOU HAVE BEEN LOOKING FOR.
flags = [
'-DNDEBUG',
- '-std=c++11',
+ '-DARB_HAVE_CTHREAD',
+ '-std=c++14',
'-x',
'c++',
'-I',
@@ -44,18 +45,17 @@ flags = [
'-I',
'.',
'-I',
- 'include',
+ 'modcc',
'-I',
- 'external',
+ 'include',
'-I',
- 'miniapp',
+ 'arbor',
'-I',
- 'modcc',
+ 'ext/json/single_include',
'-I',
- 'tests/ubench/google-benchmark/include',
+ 'build/include',
'-I',
- '/cm/shared/apps/cuda/8.0.44/include',
- '-DARB_HAVE_GPU'
+ 'aux/include',
]
# Set this to the absolute path to the folder (NOT the file!) containing the
@@ -65,18 +65,18 @@ flags = [
# Most projects will NOT need to set this to anything; you can just change the
# 'flags' list of compilation flags. Notice that YCM itself uses that approach.
compilation_database_folder = ''
-
+
if os.path.exists( compilation_database_folder ):
database = ycm_core.CompilationDatabase( compilation_database_folder )
else:
database = None
-
+
SOURCE_EXTENSIONS = [ '.cpp', '.cxx', '.cc', '.c', '.m', '.mm' ]
-
+
def DirectoryOfThisScript():
return os.path.dirname( os.path.abspath( __file__ ) )
-
-
+
+
def MakeRelativePathsInFlagsAbsolute( flags, working_directory ):
if not working_directory:
return list( flags )
@@ -85,32 +85,32 @@ def MakeRelativePathsInFlagsAbsolute( flags, working_directory ):
path_flags = [ '-isystem', '-I', '-iquote', '--sysroot=' ]
for flag in flags:
new_flag = flag
-
+
if make_next_absolute:
make_next_absolute = False
if not flag.startswith( '/' ):
new_flag = os.path.join( working_directory, flag )
-
+
for path_flag in path_flags:
if flag == path_flag:
make_next_absolute = True
break
-
+
if flag.startswith( path_flag ):
path = flag[ len( path_flag ): ]
new_flag = path_flag + os.path.join( working_directory, path )
break
-
+
if new_flag:
new_flags.append( new_flag )
return new_flags
-
-
+
+
def IsHeaderFile( filename ):
extension = os.path.splitext( filename )[ 1 ]
return extension in [ '.h', '.hxx', '.hpp', '.hh' ]
-
-
+
+
def GetCompilationInfoForFile( filename ):
# The compilation_commands.json file generated by CMake does not have entries
# for header files. So we do our best by asking the db for flags for a
@@ -127,8 +127,8 @@ def GetCompilationInfoForFile( filename ):
return compilation_info
return None
return database.GetCompilationInfoForFile( filename )
-
-
+
+
def FlagsForFile( filename, **kwargs ):
if database:
# Bear in mind that compilation_info.compiler_flags_ does NOT return a
@@ -136,11 +136,11 @@ def FlagsForFile( filename, **kwargs ):
compilation_info = GetCompilationInfoForFile( filename )
if not compilation_info:
return None
-
+
final_flags = MakeRelativePathsInFlagsAbsolute(
compilation_info.compiler_flags_,
compilation_info.compiler_working_dir_ )
-
+
# NOTE: This is just for YouCompleteMe; it's highly likely that your project
# does NOT need to remove the stdlib flag. DO NOT USE THIS IN YOUR
# ycm_extra_conf IF YOU'RE NOT 100% SURE YOU NEED IT.
@@ -151,7 +151,7 @@ def FlagsForFile( filename, **kwargs ):
else:
relative_to = DirectoryOfThisScript()
final_flags = MakeRelativePathsInFlagsAbsolute( flags, relative_to )
-
+
return {
'flags': final_flags,
'do_cache': True
diff --git a/aux/include/aux/json_params.hpp b/aux/include/aux/json_params.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f7a848e2c54d8473ea5eec715aa54fd9adc47628
--- /dev/null
+++ b/aux/include/aux/json_params.hpp
@@ -0,0 +1,42 @@
+#include <array>
+#include <exception>
+
+#include <arbor/util/optional.hpp>
+
+#include <nlohmann/json.hpp>
+
+namespace aux {
+
+// Search a json object for an entry with a given name.
+// If found, return the value and remove from json object.
+template <typename T>
+arb::util::optional<T> find_and_remove_json(const char* name, nlohmann::json& j) {
+ auto it = j.find(name);
+ if (it==j.end()) {
+ return arb::util::nullopt;
+ }
+ T value = std::move(*it);
+ j.erase(name);
+ return std::move(value);
+}
+
+template <typename T>
+void param_from_json(T& x, const char* name, nlohmann::json& j) {
+ if (auto o = find_and_remove_json<T>(name, j)) {
+ x = *o;
+ }
+}
+
+template <typename T, size_t N>
+void param_from_json(std::array<T, N>& x, const char* name, nlohmann::json& j) {
+ std::vector<T> y;
+ if (auto o = find_and_remove_json<std::vector<T>>(name, j)) {
+ y = *o;
+ if (y.size()!=N) {
+ throw std::runtime_error("parameter "+std::string(name)+" requires "+std::to_string(N)+" values");
+ }
+ std::copy(y.begin(), y.end(), x.begin());
+ }
+}
+
+} // namespace aux
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 0238e1f150f4ec95c35659accc211e89cf0f4601..56f2755605e752c76be0a221db7ebfe5d92fc16a 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -2,3 +2,4 @@ add_subdirectory(miniapp)
add_subdirectory(generators)
add_subdirectory(brunel)
add_subdirectory(bench)
+add_subdirectory(ring)
diff --git a/example/bench/parameters.cpp b/example/bench/parameters.cpp
index 74328feae748023d7ac9773050fc138b5774ce6c..b9ac4b2f8ffb0819865b3e2ba7a21b694a9e08cd 100644
--- a/example/bench/parameters.cpp
+++ b/example/bench/parameters.cpp
@@ -1,10 +1,9 @@
#include <exception>
#include <fstream>
+#include <iostream>
#include <string>
-#include <arbor/util/optional.hpp>
-
-#include <nlohmann/json.hpp>
+#include <aux/json_params.hpp>
#include "parameters.hpp"
@@ -42,23 +41,9 @@ std::ostream& operator<<(std::ostream& o, const bench_params& p) {
return o;
}
-using arb::util::optional;
-using nlohmann::json;
-
-// Search a json object for an entry with a given name.
-// If found, return the value and remove from json object.
-template<typename T>
-optional<T> extract(const char* name, json& j) {
- auto it = j.find(name);
- if (it==j.end()) {
- return arb::util::nullopt;
- }
- T value = std::move(*it);
- j.erase(name);
- return std::move(value);
-}
-
bench_params read_options(int argc, char** argv) {
+ using aux::param_from_json;
+
bench_params params;
if (argc<2) {
std::cout << "Using default parameters.\n";
@@ -79,27 +64,13 @@ bench_params read_options(int argc, char** argv) {
nlohmann::json json;
json << f;
- if (auto o = extract<std::string>("name", json)) {
- params.name = *o;
- }
- if (auto o = extract<unsigned>("num-cells", json)) {
- params.num_cells = *o;
- }
- if (auto o = extract<double>("duration", json)) {
- params.duration = *o;
- }
- if (auto o = extract<double>("min-delay", json)) {
- params.network.min_delay = *o;
- }
- if (auto o = extract<unsigned>("fan-in", json)) {
- params.network.fan_in = *o;
- }
- if (auto o = extract<double>("realtime-ratio", json)) {
- params.cell.realtime_ratio = *o;
- }
- if (auto o = extract<double>("spike-frequency", json)) {
- params.cell.spike_freq_hz = *o;
- }
+ param_from_json(params.name, "name", json);
+ param_from_json(params.num_cells, "num-cells", json);
+ param_from_json(params.duration, "duration", json);
+ param_from_json(params.network.min_delay, "min-delay", json);
+ param_from_json(params.network.fan_in, "fan-in", json);
+ param_from_json(params.cell.realtime_ratio, "realtime-ratio", json);
+ param_from_json(params.cell.spike_freq_hz, "spike-frequency", json);
for (auto it=json.begin(); it!=json.end(); ++it) {
std::cout << " Warning: unused input parameter: \"" << it.key() << "\"\n";
@@ -108,4 +79,3 @@ bench_params read_options(int argc, char** argv) {
return params;
}
-
diff --git a/example/ring/CMakeLists.txt b/example/ring/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2fd6c79f86162c50bed0fc1023c8944cad9ac3b9
--- /dev/null
+++ b/example/ring/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(ring ring.cpp)
+
+target_link_libraries(ring PRIVATE arbor arbor-aux ext-tclap ext-json)
diff --git a/example/ring/parameters.hpp b/example/ring/parameters.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f1d685735fb01fd03f55e4f1388bcc7ce1392d7
--- /dev/null
+++ b/example/ring/parameters.hpp
@@ -0,0 +1,77 @@
+#include <iostream>
+
+#include <array>
+#include <cmath>
+#include <fstream>
+#include <random>
+
+#include <arbor/mc_cell.hpp>
+
+#include <aux/json_params.hpp>
+
+// Parameters used to generate the random cell morphologies.
+struct cell_parameters {
+ cell_parameters() = default;
+
+ // Maximum number of levels in the cell (not including the soma)
+ unsigned max_depth = 5;
+
+ // The following parameters are described as ranges.
+ // The first value is at the soma, and the last value is used on the last level.
+ // Values at levels in between are found by linear interpolation.
+ std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
+ std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
+ std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
+};
+
+struct ring_params {
+ ring_params() = default;
+
+ std::string name = "default";
+ unsigned num_cells = 10;
+ double min_delay = 10;
+ double duration = 100;
+ cell_parameters cell;
+};
+
+ring_params read_options(int argc, char** argv) {
+ using aux::param_from_json;
+
+ ring_params params;
+ if (argc<2) {
+ std::cout << "Using default parameters.\n";
+ return params;
+ }
+ if (argc>2) {
+ throw std::runtime_error("More than command line one option not permitted.");
+ }
+
+ std::string fname = argv[1];
+ std::cout << "Loading parameters from file: " << fname << "\n";
+ std::ifstream f(fname);
+
+ if (!f.good()) {
+ throw std::runtime_error("Unable to open input parameter file: "+fname);
+ }
+
+ nlohmann::json json;
+ json << f;
+
+ param_from_json(params.name, "name", json);
+ param_from_json(params.num_cells, "num-cells", json);
+ param_from_json(params.duration, "duration", json);
+ param_from_json(params.min_delay, "min-delay", json);
+ param_from_json(params.cell.max_depth, "depth", json);
+ param_from_json(params.cell.branch_probs, "branch-probs", json);
+ param_from_json(params.cell.compartments, "compartments", json);
+ param_from_json(params.cell.lengths, "lengths", json);
+
+ if (!json.empty()) {
+ for (auto it=json.begin(); it!=json.end(); ++it) {
+ std::cout << " Warning: unused input parameter: \"" << it.key() << "\"\n";
+ }
+ std::cout << "\n";
+ }
+
+ return params;
+}
diff --git a/example/ring/readme.md b/example/ring/readme.md
new file mode 100644
index 0000000000000000000000000000000000000000..fb9950a54d6e738d7c929147c54b0b0cae624090
--- /dev/null
+++ b/example/ring/readme.md
@@ -0,0 +1,3 @@
+# Ring Example
+
+A miniapp that demonstrates how to describe how to build a simple ring network.
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11caa577bdd361491f8804d3982c03fbc6e85ff4
--- /dev/null
+++ b/example/ring/ring.cpp
@@ -0,0 +1,318 @@
+/*
+ * A miniapp that demonstrates how to make a ring model
+ *
+ */
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include <nlohmann/json.hpp>
+
+#include <arbor/assert_macro.hpp>
+#include <arbor/common_types.hpp>
+#include <arbor/distributed_context.hpp>
+#include <arbor/execution_context.hpp>
+#include <arbor/load_balance.hpp>
+#include <arbor/mc_cell.hpp>
+#include <arbor/profile/meter_manager.hpp>
+#include <arbor/simple_sampler.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/recipe.hpp>
+
+#include <aux/json_meter.hpp>
+
+#include "parameters.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);
+
+// Generate a cell.
+arb::mc_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params);
+
+class ring_recipe: public arb::recipe {
+public:
+ ring_recipe(unsigned num_cells, cell_parameters params, unsigned min_delay):
+ num_cells_(num_cells),
+ cell_params_(params),
+ min_delay_(min_delay)
+ {}
+
+ cell_size_type num_cells() const override {
+ return num_cells_;
+ }
+
+ arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
+ return branch_cell(gid, cell_params_);
+ }
+
+ cell_kind get_cell_kind(cell_gid_type gid) const override {
+ return cell_kind::cable1d_neuron;
+ }
+
+ // Each cell has one spike detector (at the soma).
+ cell_size_type num_sources(cell_gid_type gid) const override {
+ return 1;
+ }
+
+ // The cell has one target synapse, which will be connected to cell gid-1.
+ cell_size_type num_targets(cell_gid_type gid) const override {
+ return 1;
+ }
+
+ // Each cell has one incoming connection, from cell with gid-1.
+ std::vector<arb::cell_connection> connections_on(cell_gid_type gid) const override {
+ cell_gid_type src = gid? gid-1: num_cells_-1;
+ return {arb::cell_connection({src, 0}, {gid, 0}, event_weight_, min_delay_)};
+ }
+
+ // Return one event generator on gid 0. This generates a single event that will
+ // kick start the spiking.
+ std::vector<arb::event_generator> event_generators(cell_gid_type gid) const override {
+ std::vector<arb::event_generator> gens;
+ if (!gid) {
+ gens.push_back(arb::explicit_generator(arb::pse_vector{{{0, 0}, 0.1, 1.0}}));
+ }
+ 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 {
+ 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;
+ // Measure at the soma.
+ arb::segment_location loc(0, 0.0);
+
+ return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ }
+
+private:
+ cell_size_type num_cells_;
+ cell_parameters cell_params_;
+ double min_delay_;
+ float event_weight_ = 0.01;
+};
+
+struct cell_stats {
+ using size_type = std::uint64_t;
+ size_type ncells = 0;
+ size_type nsegs = 0;
+ size_type ncomp = 0;
+
+ cell_stats(arb::distributed_context_handle ctx, arb::recipe& r) {
+ size_type nranks = ctx->size();
+ size_type rank = ctx->id();
+ ncells = r.num_cells();
+ size_type cells_per_rank = ncells/nranks;
+ size_type b = rank*cells_per_rank;
+ size_type e = (rank==nranks-1)? ncells: (rank+1)*cells_per_rank;
+ for (size_type i=b; i<e; ++i) {
+ auto c = arb::util::any_cast<arb::mc_cell>(r.get_cell_description(i));
+ nsegs += c.num_segments();
+ ncomp += c.num_compartments();
+ }
+ nsegs = ctx->sum(nsegs);
+ ncomp = ctx->sum(ncomp);
+ }
+
+ friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) {
+ return o << "cell stats: "
+ << s.ncells << " cells; "
+ << s.nsegs << " segments; "
+ << s.ncomp << " compartments.";
+ }
+};
+
+
+int main(int argc, char** argv) {
+ // default serial context
+ arb::execution_context context;
+
+ try {
+#ifdef ARB_MPI_ENABLED
+ aux::with_mpi guard(argc, argv, false);
+ context.distributed = mpi_context(MPI_COMM_WORLD);
+#endif
+#ifdef ARB_PROFILE_ENABLED
+ profile::profiler_initialize(context.thread_pool);
+#endif
+
+ const bool root = context.distributed->id()==0;
+
+ auto params = read_options(argc, argv);
+
+ arb::profile::meter_manager meters(context.distributed);
+ meters.start();
+
+ // Create an instance of our recipe.
+ ring_recipe recipe(params.num_cells, params.cell, params.min_delay);
+ cell_stats stats(context.distributed, recipe);
+ std::cout << stats << "\n";
+
+ // Use a node description that uses the number of threads used by the
+ // threading back end, and 1 gpu if available.
+ arb::proc_allocation nd = arb::local_allocation(context);
+ nd.num_gpus = nd.num_gpus>=1? 1: 0;
+
+ auto decomp = arb::partition_load_balance(recipe, nd, context);
+
+ // Construct the model.
+ arb::simulation sim(recipe, decomp, context);
+
+ // 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
+ sim.add_sampler(arb::one_probe(probe_id), sched, arb::make_simple_sampler(voltage));
+
+ // Set up recording of spikes to a vector on the root process.
+ std::vector<arb::spike> recorded_spikes;
+ if (root) {
+ sim.set_global_spike_callback(
+ [&recorded_spikes](const std::vector<arb::spike>& spikes) {
+ recorded_spikes.insert(recorded_spikes.end(), spikes.begin(), spikes.end());
+ });
+ }
+
+ meters.checkpoint("model-init");
+
+ // Run the simulation for 100 ms, with time steps of 0.025 ms.
+ sim.run(params.duration, 0.025);
+
+ meters.checkpoint("model-run");
+
+ auto ns = sim.num_spikes();
+ std::cout << "\n" << ns << " spikes generated at rate of "
+ << params.duration/ns << " ms between spikes\n";
+
+ // Write spikes to file
+ if (root) {
+ std::ofstream fid("spikes.gdf");
+ if (!fid.good()) {
+ std::cerr << "Warning: unable to open file spikes.gdf for spike output\n";
+ }
+ else {
+ char linebuf[45];
+ for (auto spike: recorded_spikes) {
+ auto n = std::snprintf(
+ linebuf, sizeof(linebuf), "%u %.4f\n",
+ unsigned{spike.source.gid}, float(spike.time));
+ fid.write(linebuf, n);
+ }
+ }
+ }
+
+ // Write the samples to a json file.
+ write_trace_json(voltage);
+
+ auto report = arb::profile::make_meter_report(meters);
+ std::cout << report;
+ }
+ catch (std::exception& e) {
+ std::cerr << "exception caught in ring miniapp:\n" << e.what() << "\n";
+ return 1;
+ }
+
+ return 0;
+}
+
+void write_trace_json(const arb::trace_data<double>& trace) {
+ std::string path = "./voltages.json";
+
+ nlohmann::json json;
+ json["name"] = "ring 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";
+}
+
+// Helper used to interpolate in branch_cell.
+template <typename T>
+double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
+ double p = i * 1./(n-1);
+ double r0 = r[0];
+ double r1 = r[1];
+ return r[0] + p*(r1-r0);
+}
+
+arb::mc_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
+ arb::mc_cell cell;
+
+ // Add soma.
+ auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm².
+ soma->rL = 100;
+ soma->add_mechanism("hh");
+
+ std::vector<std::vector<unsigned>> levels;
+ levels.push_back({0});
+
+ // Standard mersenne_twister_engine seeded with gid.
+ std::mt19937 gen(gid);
+ std::uniform_real_distribution<double> dis(0, 1);
+
+ double dend_radius = 0.5; // Diameter of 1 μm for each cable.
+
+ unsigned nsec = 1;
+ for (unsigned i=0; i<params.max_depth; ++i) {
+ // Branch prob at this level.
+ double bp = interp(params.branch_probs, i, params.max_depth);
+ // Length at this level.
+ double l = interp(params.lengths, i, params.max_depth);
+ // Number of compartments at this level.
+ unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
+
+ std::vector<unsigned> sec_ids;
+ for (unsigned sec: levels[i]) {
+ for (unsigned j=0; j<2; ++j) {
+ if (dis(gen)<bp) {
+ sec_ids.push_back(nsec++);
+ auto dend = cell.add_cable(sec, arb::section_kind::dendrite, dend_radius, dend_radius, l);
+ dend->set_compartments(nc);
+ dend->add_mechanism("pas");
+ dend->rL = 100;
+ }
+ }
+ }
+ if (sec_ids.empty()) {
+ break;
+ }
+ levels.push_back(sec_ids);
+ }
+
+ // Add spike threshold detector at the soma.
+ cell.add_detector({0,0}, 10);
+
+ // Add a synapse to the mid point of the first dendrite.
+ cell.add_synapse({1, 0.5}, "expsyn");
+
+ return cell;
+}
+