diff --git a/.gitignore b/.gitignore
index e9d54cb796031fdcbba7935bec51511be3f6de6b..b710bef2b72c6b9f6eea87de726bf180c4f3a107 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,3 +66,5 @@ mechanisms/*.hpp
# build path
build
+
+commit.msg
diff --git a/.gitmodules b/.gitmodules
index ae780729fd3c2e22a3a6f5029fb0b611977ccfeb..86e0c20fe68579e7b0f1ccdd1dba498cf16b13fe 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -7,3 +7,6 @@
[submodule "modparser"]
path = external/modparser
url = git@github.com:eth-cscs/modparser.git
+[submodule "external/tclap"]
+ path = external/tclap
+ url = https://github.com/eile/tclap.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc45593eb2451be9818e2b75671d91341647e4e4..5d082f12be8d6c555372be33a8737959f4e2a08b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -65,6 +65,7 @@ externalproject_add(modparser
include_directories(${CMAKE_SOURCE_DIR}/external)
+include_directories(${CMAKE_SOURCE_DIR}/external/tclap/include)
include_directories(${CMAKE_SOURCE_DIR}/include)
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${CMAKE_SOURCE_DIR}/miniapp)
diff --git a/commit.msg b/commit.msg
deleted file mode 100644
index 3abb1935698a6efe63d5ea651d8c7d5c017541fe..0000000000000000000000000000000000000000
--- a/commit.msg
+++ /dev/null
@@ -1,9 +0,0 @@
-* change `probe_sort` enum to scoped enum
- - renamed to `probeKind`
- - refactored out of class
- - updated coding guidelines wiki with enum rules
-* refactor `std::pair` into structs with meaningfull name types in `cell.hpp`
- - not just probes: stimulii and detectors too.
-* add profiling region around sampling in `cell_group`
-* change output data format for traces to json
-* remove white space at end of lines (looking at you Sam)
diff --git a/external/tclap b/external/tclap
new file mode 160000
index 0000000000000000000000000000000000000000..f41dcb5ce3d063c9fe95623193bba693338f3edb
--- /dev/null
+++ b/external/tclap
@@ -0,0 +1 @@
+Subproject commit f41dcb5ce3d063c9fe95623193bba693338f3edb
diff --git a/miniapp/io.cpp b/miniapp/io.cpp
index 887848b549b7bad46574b92ac9fcde74faae5d5c..8367a1403a1ff58ab657239deafc0fd024731295 100644
--- a/miniapp/io.cpp
+++ b/miniapp/io.cpp
@@ -1,21 +1,95 @@
+#include <fstream>
+#include <exception>
+
+#include <tclap/CmdLine.h>
+
#include "io.hpp"
namespace nest {
namespace mc {
namespace io {
-// read simulation options from json file with name fname
-// for now this is just a placeholder
-options read_options(std::string fname) {
- // 10 cells, 1 synapses per cell, 10 compartments per segment
- return {200, 1, 100};
+/// read simulation options from json file with name fname
+/// if file name is empty or if file is not a valid json file a default
+/// set of parameters is returned :
+/// 1000 cells, 500 synapses per cell, 100 compartments per segment
+cl_options read_options(int argc, char** argv) {
+
+ // set default options
+ const cl_options default_options{"", 1000, 500, 100};
+
+ cl_options options;
+ // parse command line arguments
+ try {
+ TCLAP::CmdLine cmd("mod2c performance benchmark harness", ' ', "0.1");
+
+ TCLAP::ValueArg<uint32_t> ncells_arg(
+ "n", "ncells", "total number of cells in the model",
+ false, 1000,"non negative integer");
+ TCLAP::ValueArg<uint32_t> nsynapses_arg(
+ "s", "nsynapses", "number of synapses per cell",
+ false, 500,"non negative integer");
+ TCLAP::ValueArg<uint32_t> ncompartments_arg(
+ "c", "ncompartments", "number of compartments per segment",
+ false, 100,"non negative integer");
+ TCLAP::ValueArg<std::string> ifile_arg(
+ "i", "ifile", "number of compartments per segment",
+ false, "","file name string");
+
+ cmd.add(ncells_arg);
+ cmd.add(nsynapses_arg);
+ cmd.add(ncompartments_arg);
+ cmd.add(ifile_arg);
+
+ cmd.parse(argc, argv);
+
+ options.cells = ncells_arg.getValue();
+ options.synapses_per_cell = nsynapses_arg.getValue();
+ options.compartments_per_segment = ncompartments_arg.getValue();
+ options.ifname = ifile_arg.getValue();
+ }
+ // catch any exceptions in command line handling
+ catch(TCLAP::ArgException &e) {
+ std::cerr << "error: parsing command line arguments:\n "
+ << e.error() << " for arg " << e.argId() << "\n";
+ exit(1);
+ }
+
+ if(options.ifname == "") {
+ return options;
+ }
+ else {
+ std::ifstream fid(options.ifname, std::ifstream::in);
+
+ if(fid.is_open()) {
+ // read json data in input file
+ nlohmann::json fopts;
+ fid >> fopts;
+
+ try {
+ options.cells = fopts["cells"];
+ options.synapses_per_cell = fopts["synapses"];
+ options.compartments_per_segment = fopts["compartments"];
+ }
+ catch(std::exception e) {
+ std::cerr << "error: unable to open parameters in "
+ << options.ifname << ", using defaults instead\n"
+ << "... " << e.what() << "\n";
+ exit(1);
+ }
+ return options;
+ }
+ }
+
+ return default_options;
}
-std::ostream& operator<<(std::ostream& o, const options& opt) {
+std::ostream& operator<<(std::ostream& o, const cl_options& options) {
o << "simultion options:\n";
- o << " cells : " << opt.cells << "\n";
- o << " compartments/segment : " << opt.compartments_per_segment << "\n";
- o << " synapses/cell : " << opt.synapses_per_cell << "\n";
+ o << " cells : " << options.cells << "\n";
+ o << " compartments/segment : " << options.compartments_per_segment << "\n";
+ o << " synapses/cell : " << options.synapses_per_cell << "\n";
+ o << " input file name : " << options.ifname << "\n";
return o;
}
diff --git a/miniapp/io.hpp b/miniapp/io.hpp
index e8a8762ceb34b7d7148dfb7eb831217ac633d3d3..fea4dbd911e4f7214500349c07611d1e8f93cf9d 100644
--- a/miniapp/io.hpp
+++ b/miniapp/io.hpp
@@ -7,15 +7,16 @@ namespace mc {
namespace io {
// holds the options for a simulation run
-struct options {
- int cells;
- int synapses_per_cell;
- int compartments_per_segment;
+struct cl_options {
+ std::string ifname;
+ uint32_t cells;
+ uint32_t synapses_per_cell;
+ uint32_t compartments_per_segment;
};
-std::ostream& operator<<(std::ostream& o, const options& opt);
+std::ostream& operator<<(std::ostream& o, const cl_options& opt);
-options read_options(std::string fname);
+cl_options read_options(int argc, char** argv);
} // namespace io
} // namespace mc
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index 74b551ff4f07fe795514caafde901de1721108dc..97824cfe4dbec20272c46a9db0939f096ab01dc5 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -32,14 +32,14 @@ struct model {
communicator_type communicator;
std::vector<cell_group> cell_groups;
- int num_groups() const {
+ unsigned num_groups() const {
return cell_groups.size();
}
void run(double tfinal, double dt) {
auto t = 0.;
auto delta = communicator.min_delay();
- while(t<tfinal) {
+ while (t<tfinal) {
mc::threading::parallel_for::apply(
0, num_groups(),
[&](int i) {
@@ -68,7 +68,7 @@ struct model {
// calculate the source and synapse distribution serially
std::vector<id_type> target_counts(num_groups());
std::vector<id_type> source_counts(num_groups());
- for (auto i=0; i<num_groups(); ++i) {
+ for (auto i=0u; i<num_groups(); ++i) {
target_counts[i] = cell_groups[i].cell().synapses()->size();
source_counts[i] = cell_groups[i].spike_sources().size();
}
@@ -87,9 +87,13 @@ struct model {
auto com_policy = communicator.communication_policy();
auto global_source_map = com_policy.make_map(source_map.back());
auto domain_idx = communicator.domain_id();
- for (auto i=0; i<num_groups(); ++i) {
- cell_groups[i].set_source_gids(source_map[i]+global_source_map[domain_idx]);
- cell_groups[i].set_target_gids(target_map[i]+communicator.target_gid_from_group_lid(0));
+ for (auto i=0u; i<num_groups(); ++i) {
+ cell_groups[i].set_source_gids(
+ source_map[i]+global_source_map[domain_idx]
+ );
+ cell_groups[i].set_target_gids(
+ target_map[i]+communicator.target_gid_from_group_lid(0)
+ );
}
mc::util::profiler_leave(2);
}
@@ -166,10 +170,10 @@ struct model {
namespace parameters {
namespace synapses {
// synapse delay
- constexpr double delay = 5.0; // ms
+ constexpr float delay = 20.0; // ms
// connection weight
- constexpr double weight = 0.005; // uS
+ constexpr double weight_per_cell = 0.3; // uS
}
}
@@ -187,9 +191,7 @@ void setup();
cell_group make_lowered_cell(int cell_index, const mc::cell& c);
/// models
-void ring_model(nest::mc::io::options& opt, model& m);
-void all_to_all_model(nest::mc::io::options& opt, model& m);
-
+void all_to_all_model(nest::mc::io::cl_options& opt, model& m);
///////////////////////////////////////
// main
@@ -200,11 +202,11 @@ int main(int argc, char** argv) {
setup();
// read parameters
- mc::io::options opt;
+ mc::io::cl_options options;
try {
- opt = mc::io::read_options("");
+ options = mc::io::read_options(argc, argv);
if (!global_policy::id()) {
- std::cout << opt << "\n";
+ std::cout << options << "\n";
}
}
catch (std::exception& e) {
@@ -213,18 +215,29 @@ int main(int argc, char** argv) {
}
model m;
- all_to_all_model(opt, m);
+ all_to_all_model(options, m);
//
// time stepping
//
- auto tfinal = 20.;
- auto dt = 0.01;
+ auto tfinal = 50.;
+ auto dt = 0.025;
auto id = m.communicator.domain_id();
- if (!id) {
- m.communicator.add_spike({0, 5});
+ if (id==0) {
+ std::cout << "\n";
+ std::cout << ":: simulation to " << tfinal << " ms in "
+ << std::ceil(tfinal / dt) << " steps of "
+ << dt << " ms\n";
+ }
+
+ // add some spikes to the system to start it
+ auto first = m.communicator.group_gid_first(id);
+ first += 20 - (first%20); // round up to multiple of 20
+ auto last = m.communicator.group_gid_first(id+1);
+ for (auto i=first; i<last; i+=20) {
+ m.communicator.add_spike({i, 0});
}
m.run(tfinal, dt);
@@ -234,44 +247,48 @@ int main(int argc, char** argv) {
if (!id) {
std::cout << "there were " << m.communicator.num_spikes() << " spikes\n";
}
- m.dump_traces();
-#ifdef SPLAT
- if (!global_policy::id()) {
- //for (auto i=0u; i<m.cell_groups.size(); ++i) {
- m.cell_groups[0].splat("cell0.txt");
- m.cell_groups[1].splat("cell1.txt");
- m.cell_groups[2].splat("cell2.txt");
- //}
- }
-#endif
+ m.dump_traces();
}
///////////////////////////////////////
// models
///////////////////////////////////////
-void ring_model(nest::mc::io::options& opt, model& m) {
+void all_to_all_model(nest::mc::io::cl_options& options, model& m) {
//
// make cells
//
+ // calculate how many synapses per cell
+ auto synapses_per_cell =
+ std::min(options.synapses_per_cell, options.cells-1);
+ auto is_all_to_all = synapses_per_cell == (options.cells-1);
+
// make a basic cell
- auto basic_cell = make_cell(opt.compartments_per_segment, 1);
+ auto basic_cell =
+ make_cell(options.compartments_per_segment, synapses_per_cell);
+
+ auto num_domains = global_policy::size();
+ auto domain_id = global_policy::id();
// make a vector for storing all of the cells
- m.cell_groups = std::vector<cell_group>(opt.cells);
+ id_type ncell_global = options.cells;
+ id_type ncell_local = ncell_global / num_domains;
+ int remainder = ncell_global - (ncell_local*num_domains);
+ if (domain_id<remainder) {
+ ncell_local++;
+ }
+
+ m.cell_groups = std::vector<cell_group>(ncell_local);
// initialize the cells in parallel
mc::threading::parallel_for::apply(
- 0, opt.cells,
+ 0, ncell_local,
[&](int i) {
- // initialize cell
- mc::util::profiler_enter("setup");
- mc::util::profiler_enter("make cell");
+ mc::util::profiler_enter("setup", "cells");
m.cell_groups[i] = make_lowered_cell(i, basic_cell);
- mc::util::profiler_leave();
- mc::util::profiler_leave();
+ mc::util::profiler_leave(2);
}
);
@@ -280,48 +297,52 @@ void ring_model(nest::mc::io::options& opt, model& m) {
//
m.init_communicator();
- for (auto i=0u; i<(id_type)opt.cells; ++i) {
- m.communicator.add_connection({
- i, (i+1)%opt.cells,
- parameters::synapses::weight, parameters::synapses::delay
- });
- }
-
- m.update_gids();
-}
+ mc::util::profiler_enter("setup", "connections");
-void all_to_all_model(nest::mc::io::options& opt, model& m) {
- //
- // make cells
- //
+ // RNG distributions for connection delays and source cell ids
+ auto weight_distribution = std::exponential_distribution<float>(0.75);
+ auto source_distribution =
+ std::uniform_int_distribution<uint32_t>(0u, options.cells-2);
- // make a basic cell
- auto basic_cell = make_cell(opt.compartments_per_segment, opt.cells-1);
+ // calculate the weight of synaptic connections, which is chosen so that
+ // the sum of all synaptic weights on a cell is
+ // parameters::synapses::weight_per_cell
+ float weight = parameters::synapses::weight_per_cell / synapses_per_cell;
- // make a vector for storing all of the cells
- id_type ncell_global = opt.cells;
- id_type ncell_local = ncell_global / m.communicator.num_domains();
- int remainder = ncell_global - (ncell_local*m.communicator.num_domains());
- if (m.communicator.domain_id()<remainder) {
- ncell_local++;
+ // loop over each local cell and build the list of synapse connections
+ // that terminate on the cell
+ for (auto lid=0u; lid<ncell_local; ++lid) {
+ auto target = m.communicator.target_gid_from_group_lid(lid);
+ auto gid = m.communicator.group_gid_from_group_lid(lid);
+ auto gen = std::mt19937(gid); // seed with cell gid for reproducability
+ // add synapses to cell
+ auto i = 0u;
+ auto cells_added = 0u;
+ while (cells_added < synapses_per_cell) {
+ auto source = is_all_to_all ? i : source_distribution(gen);
+ if (gid!=source) {
+ m.communicator.add_connection({
+ source, target++, weight,
+ parameters::synapses::delay + weight_distribution(gen)
+ });
+ cells_added++;
+ }
+ ++i;
+ }
}
- m.cell_groups = std::vector<cell_group>(ncell_local);
+ m.communicator.construct();
- // initialize the cells in parallel
- mc::threading::parallel_for::apply(
- 0, ncell_local,
- [&](int i) {
- mc::util::profiler_enter("setup", "cells");
- m.cell_groups[i] = make_lowered_cell(i, basic_cell);
- mc::util::profiler_leave(2);
- }
- );
+ //for (auto con : m.communicator.connections()) std::cout << con << "\n";
+
+ m.update_gids();
//
- // network creation
+ // setup probes
//
- m.init_communicator();
+
+ mc::util::profiler_leave();
+ mc::util::profiler_enter("probes");
// monitor soma and dendrite on a few cells
float sample_dt = 0.1;
@@ -335,30 +356,15 @@ void all_to_all_model(nest::mc::io::options& opt, model& m) {
auto probe_soma = m.cell_groups[lid].probe_gid_range().first;
auto probe_dend = probe_soma+1;
- m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_soma, "vsoma", gid, sample_dt));
- m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_dend, "vdend", gid, sample_dt));
- }
-
- mc::util::profiler_enter("setup", "connections");
- // lid is local cell/group id
- for (auto lid=0u; lid<ncell_local; ++lid) {
- auto target = m.communicator.target_gid_from_group_lid(lid);
- auto gid = m.communicator.group_gid_from_group_lid(lid);
- // tid is global cell/group id
- for (auto tid=0u; tid<ncell_global; ++tid) {
- if (gid!=tid) {
- m.communicator.add_connection({
- tid, target++,
- parameters::synapses::weight, parameters::synapses::delay
- });
- }
- }
+ m.cell_groups[lid].add_sampler(
+ m.make_simple_sampler(probe_soma, "vsoma", gid, sample_dt)
+ );
+ m.cell_groups[lid].add_sampler(
+ m.make_simple_sampler(probe_dend, "vdend", gid, sample_dt)
+ );
}
- m.communicator.construct();
- mc::util::profiler_leave(2);
-
- m.update_gids();
+ mc::util::profiler_leave(2);
}
///////////////////////////////////////
@@ -390,8 +396,8 @@ mc::cell make_cell(int compartments_per_segment, int num_synapses) {
// add dendrite of length 200 um and diameter 1 um with passive channel
std::vector<mc::cable_segment*> dendrites;
dendrites.push_back(cell.add_cable(0, mc::segmentKind::dendrite, 0.5, 0.5, 200));
- //dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25, 100));
- //dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25, 100));
+ dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+ dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
for (auto d : dendrites) {
d->add_mechanism(mc::pas_parameters());
@@ -399,13 +405,14 @@ mc::cell make_cell(int compartments_per_segment, int num_synapses) {
d->mechanism("membrane").set("r_L", 100);
}
- // add stimulus
- //cell.add_stimulus({1,1}, {5., 80., 0.3});
-
- cell.add_detector({0,0}, 30);
+ cell.add_detector({0,0}, 20);
+ auto gen = std::mt19937();
+ auto distribution = std::uniform_real_distribution<float>(0.f, 1.0f);
+ // distribute the synapses at random locations the terminal dendrites in a
+ // round robin manner
for (auto i=0; i<num_synapses; ++i) {
- cell.add_synapse({1, 0.5});
+ cell.add_synapse({2+(i%2), distribution(gen)});
}
// add probes:
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 70bfdeb13123c780744c5875ae125bf631bf6142..e3582b2362c8fe48984a41c21373f3755781f522 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -102,6 +102,9 @@ public:
// integrate cell state
cell_.advance(tnext - cell_.time());
+ if(!cell_.is_physical_solution()) {
+ std::cerr << "warning: solution out of bounds\n";
+ }
nest::mc::util::profiler_enter("events");
// check for new spikes
diff --git a/src/fvm_cell.hpp b/src/fvm_cell.hpp
index 90e7f7be0d373274aba0b4ff28b908465e674cb5..dde7637cea2221aa2e7d519d4a7c6bccdcc9599d 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -96,7 +96,6 @@ public:
std::vector<mechanism_type>& mechanisms() { return mechanisms_; }
/// return reference to list of ions
- //std::map<mechanisms::ionKind, ion_type> ions_;
std::map<mechanisms::ionKind, ion_type>& ions() { return ions_; }
std::map<mechanisms::ionKind, ion_type> const& ions() const { return ions_; }
@@ -133,6 +132,15 @@ public:
/// returns voltage at a segment location
value_type voltage(segment_location loc) const;
+ /// flags if solution is physically realistic.
+ /// here we define physically realistic as the voltage being within reasonable bounds.
+ /// use a simple test of the voltage at the soma is reasonable, i.e. in the range
+ /// v_soma \in (-1000mv, 1000mv)
+ bool is_physical_solution() const {
+ auto v = voltage_[0];
+ return (v>-1000.) && (v<1000.);
+ }
+
value_type time() const { return t_; }
value_type probe(uint32_t i) const {