From b7623d13e4d49146006cf4e2024e4433cb936350 Mon Sep 17 00:00:00 2001
From: Sam Yates <yates@cscs.ch>
Date: Thu, 2 Nov 2017 17:18:51 +0100
Subject: [PATCH] Add mechanism parameter setting/new implementation. (#377)
Fixes #350
* Replace parameter_list with mechanism_spec.
* Add prototype for mechanism parameter schema checking.
* Allow mechanism weights to be set after construction.
* Combine range parameters on density mechanisms by linear contribution in CVs.
* Cable segment electrical parameters are now member variables.
* Publish mechanism parameter information through new method `mechanism::field_info`; note this will be replaced/improved in upcoming dynamic mechanism catalog work.
* Access mechanism parameter scalars and range data via `mechanism::field_view_ptr` and `mechanism::field_value_ptr` methods.
* Allow mechanism 'global' parameters to be set via a method of specializing mechanisms (and giving them corresponding aliases).
* Extend recipe interface to allow querying of per-cell-kind global information for use by cell group implementations.
* Add unit tests for above - note that linear density mechanism parameter test is tightly coupled with the FVM discretization scheme.
---
miniapp/miniapp_recipes.cpp | 8 +-
modcc/blocks.hpp | 19 +-
modcc/cprinter.cpp | 117 +++++++-
modcc/cudaprinter.cpp | 77 ++++++
modcc/parser.cpp | 72 ++++-
modcc/parser.hpp | 3 +
src/CMakeLists.txt | 1 -
src/backends/fvm.hpp | 2 +-
src/backends/gpu/fvm.hpp | 12 +-
src/backends/multicore/fvm.hpp | 14 +-
src/backends/multicore/stimulus.hpp | 3 -
src/cell.hpp | 20 +-
src/fvm_multicell.hpp | 319 +++++++++++++++------
src/mechanism.hpp | 57 +++-
src/mechinfo.hpp | 108 ++++++++
src/memory/array.hpp | 19 --
src/parameter_list.cpp | 108 --------
src/parameter_list.hpp | 229 ---------------
src/recipe.hpp | 8 +-
src/segment.hpp | 161 +++--------
src/util/simple_table.hpp | 44 +++
tests/common_cells.hpp | 69 ++---
tests/simple_recipes.hpp | 14 +
tests/unit/CMakeLists.txt | 2 +-
tests/unit/common.hpp | 12 +
tests/unit/test.cpp | 4 +
tests/unit/test_cell.cpp | 3 +-
tests/unit/test_fvm_multi.cpp | 369 +++++++++++++++++++++++--
tests/unit/test_matrix.cpp | 8 +-
tests/unit/test_matrix.cu | 2 +-
tests/unit/test_mechinfo.cpp | 45 +++
tests/unit/test_parameters.cpp | 41 ---
tests/unit/test_synapses.cpp | 17 +-
tests/validation/validate_kinetic.cpp | 13 +-
tests/validation/validate_synapses.cpp | 2 +-
35 files changed, 1281 insertions(+), 721 deletions(-)
create mode 100644 src/mechinfo.hpp
delete mode 100644 src/parameter_list.cpp
delete mode 100644 src/parameter_list.hpp
create mode 100644 src/util/simple_table.hpp
create mode 100644 tests/unit/test_mechinfo.cpp
delete mode 100644 tests/unit/test_parameters.cpp
diff --git a/miniapp/miniapp_recipes.cpp b/miniapp/miniapp_recipes.cpp
index 91384d47..cae791fc 100644
--- a/miniapp/miniapp_recipes.cpp
+++ b/miniapp/miniapp_recipes.cpp
@@ -37,12 +37,12 @@ cell make_basic_cell(
}
if (segment->is_dendrite()) {
- segment->add_mechanism(pas_parameters());
- segment->mechanism("membrane").set("r_L", 100);
+ segment->add_mechanism("pas");
+ segment->rL = 100;
}
}
- cell.soma()->add_mechanism(hh_parameters());
+ cell.soma()->add_mechanism("hh");
cell.add_detector({0,0}, 20);
auto distribution = std::uniform_real_distribution<float>(0.f, 1.0f);
@@ -62,7 +62,7 @@ cell make_basic_cell(
EXPECTS(!terminals.empty());
- arb::parameter_list syn_default(syn_type);
+ arb::mechanism_spec syn_default(syn_type);
for (unsigned i=0; i<num_synapses; ++i) {
unsigned id = terminals[i%terminals.size()];
cell.add_synapse({id, distribution(rng)}, syn_default);
diff --git a/modcc/blocks.hpp b/modcc/blocks.hpp
index 3ea27608..fc26c91d 100644
--- a/modcc/blocks.hpp
+++ b/modcc/blocks.hpp
@@ -33,6 +33,8 @@ struct Id {
// string == no value
unit_tokens units;
+ std::pair<Token, Token> range; // empty component => no range set
+
Id(Token const& t, std::string const& v, unit_tokens const& u)
: token(t), value(v), units(u)
{}
@@ -40,7 +42,22 @@ struct Id {
Id() {}
bool has_value() const {
- return value.size()>0;
+ return !value.empty();
+ }
+
+ bool has_range() const {
+ return !range.first.spelling.empty();
+ }
+
+ std::string unit_string() const {
+ std::string u;
+ for (auto& t: units) {
+ if (!u.empty()) {
+ u += ' ';
+ }
+ u += t.spelling;
+ }
+ return u;
}
std::string const& name() const {
diff --git a/modcc/cprinter.cpp b/modcc/cprinter.cpp
index 8e7ad9b0..7c935f61 100644
--- a/modcc/cprinter.cpp
+++ b/modcc/cprinter.cpp
@@ -1,5 +1,6 @@
#include <algorithm>
#include <string>
+#include <unordered_set>
#include "cprinter.hpp"
#include "lexer.hpp"
@@ -19,6 +20,7 @@ std::string CPrinter::emit_source() {
// and a list of all scalar types
std::vector<VariableExpression*> scalar_variables;
std::vector<VariableExpression*> array_variables;
+
for(auto& sym: module_->symbols()) {
if(auto var = sym.second->is_variable()) {
if(var->is_range()) {
@@ -132,12 +134,21 @@ std::string CPrinter::emit_source() {
if (module_->kind() == moduleKind::density) {
text_.add_line("// add the user-supplied weights for converting from current density");
text_.add_line("// to per-compartment current in nA");
+ text_.add_line("if (weights.size()) {");
+ text_.increase_indentation();
if(optimize_) {
text_.add_line("memory::copy(weights, view(weights_, size()));");
}
else {
text_.add_line("memory::copy(weights, weights_(0, size()));");
}
+ text_.decrease_indentation();
+ text_.add_line("}");
+ text_.add_line("else {");
+ text_.increase_indentation();
+ text_.add_line("memory::fill(weights_, 1.0);");
+ text_.decrease_indentation();
+ text_.add_line("}");
text_.add_line();
}
@@ -179,10 +190,6 @@ std::string CPrinter::emit_source() {
text_.add_line("}");
text_.add_line();
- text_.add_line("void set_params() override {");
- text_.add_line("}");
- text_.add_line();
-
text_.add_line("std::string name() const override {");
text_.increase_indentation();
text_.add_line("return \"" + module_name + "\";");
@@ -200,6 +207,16 @@ std::string CPrinter::emit_source() {
text_.add_line("}");
text_.add_line();
+ // Override `set_weights` method only for density mechanisms.
+ if (module_->kind() == moduleKind::density) {
+ text_.add_line("void set_weights(array&& weights) override {");
+ text_.increase_indentation();
+ text_.add_line("memory::copy(weights, weights_(0, size()));");
+ text_.decrease_indentation();
+ text_.add_line("}");
+ text_.add_line();
+ }
+
// return true/false indicating if cell has dependency on k
auto const& ions = module_->neuron_block().ions;
auto find_ion = [&ions] (ionKind k) {
@@ -351,6 +368,97 @@ std::string CPrinter::emit_source() {
text_.add_line();
}
+ // TODO: replace field_info() generation implemenation with separate schema info generation
+ // as per #349.
+ auto field_info_string = [](const std::string& kind, const Id& id) {
+ return "field_spec{field_spec::" + kind + ", " +
+ "\"" + id.unit_string() + "\", " +
+ (id.has_value()? id.value: "0") +
+ (id.has_range()? ", " + id.range.first.spelling + "," + id.range.second.spelling: "") +
+ "}";
+ };
+
+ std::unordered_set<std::string> scalar_set;
+ for (auto& v: scalar_variables) {
+ scalar_set.insert(v->name());
+ }
+
+ std::vector<Id> global_param_ids;
+ std::vector<Id> instance_param_ids;
+
+ for (const Id& id: module_->parameter_block().parameters) {
+ auto var = id.token.spelling;
+ (scalar_set.count(var)? global_param_ids: instance_param_ids).push_back(id);
+ }
+ const std::vector<Id>& state_ids = module_->state_block().state_variables;
+
+ text_.add_line("util::optional<field_spec> field_info(const char* id) const /* override */ {");
+ text_.increase_indentation();
+ text_.add_line("static const std::pair<const char*, field_spec> field_tbl[] = {");
+ text_.increase_indentation();
+ for (const auto& id: global_param_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", "+field_info_string("global",id )+"},");
+ }
+ for (const auto& id: instance_param_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", "+field_info_string("parameter", id)+"},");
+ }
+ for (const auto& id: state_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", "+field_info_string("state", id)+"},");
+ }
+ text_.decrease_indentation();
+ text_.add_line("};");
+ text_.add_line();
+ text_.add_line("auto* info = util::table_lookup(field_tbl, id);");
+ text_.add_line("return info? util::just(*info): util::nothing;");
+ text_.decrease_indentation();
+ text_.add_line("}");
+ text_.add_line();
+
+ if (!instance_param_ids.empty() || !state_ids.empty()) {
+ text_.add_line("view base::* field_view_ptr(const char* id) const override {");
+ text_.increase_indentation();
+ text_.add_line("static const std::pair<const char*, view "+class_name+"::*> field_tbl[] = {");
+ text_.increase_indentation();
+ for (const auto& id: instance_param_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ for (const auto& id: state_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ text_.decrease_indentation();
+ text_.add_line("};");
+ text_.add_line();
+ text_.add_line("auto* pptr = util::table_lookup(field_tbl, id);");
+ text_.add_line("return pptr? static_cast<view base::*>(*pptr): nullptr;");
+ text_.decrease_indentation();
+ text_.add_line("}");
+ text_.add_line();
+ }
+
+ if (!global_param_ids.empty()) {
+ text_.add_line("value_type base::* field_value_ptr(const char* id) const override {");
+ text_.increase_indentation();
+ text_.add_line("static const std::pair<const char*, value_type "+class_name+"::*> field_tbl[] = {");
+ text_.increase_indentation();
+ for (const auto& id: global_param_ids) {
+ auto var = id.token.spelling;
+ text_.add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ text_.decrease_indentation();
+ text_.add_line("};");
+ text_.add_line();
+ text_.add_line("auto* pptr = util::table_lookup(field_tbl, id);");
+ text_.add_line("return pptr? static_cast<value_type base::*>(*pptr): nullptr;");
+ text_.decrease_indentation();
+ text_.add_line("}");
+ text_.add_line();
+ }
+
//////////////////////////////////////////////
//////////////////////////////////////////////
@@ -409,6 +517,7 @@ void CPrinter::emit_headers() {
text_.add_line("#include <backends/event.hpp>");
text_.add_line("#include <backends/multi_event_stream_state.hpp>");
text_.add_line("#include <util/pprintf.hpp>");
+ text_.add_line("#include <util/simple_table.hpp>");
text_.add_line();
}
diff --git a/modcc/cudaprinter.cpp b/modcc/cudaprinter.cpp
index 3eeecf84..1d9be13e 100644
--- a/modcc/cudaprinter.cpp
+++ b/modcc/cudaprinter.cpp
@@ -1,4 +1,6 @@
#include <algorithm>
+#include <string>
+#include <unordered_set>
#include "cudaprinter.hpp"
#include "lexer.hpp"
@@ -46,6 +48,7 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
buffer().add_line("#include <backends/fvm_types.hpp>");
buffer().add_line("#include <backends/multi_event_stream_state.hpp>");
buffer().add_line("#include <backends/gpu/kernels/detail.hpp>");
+ buffer().add_line("#include <util/simple_table.hpp>");
buffer().add_line();
buffer().add_line("namespace arb { namespace gpu{");
@@ -317,7 +320,16 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
if (m.kind() == moduleKind::density) {
buffer().add_line("// add the user-supplied weights for converting from current density");
buffer().add_line("// to per-compartment current in nA");
+ buffer().add_line("if (weights.size()) {");
+ buffer().increase_indentation();
buffer().add_line("memory::copy(weights, weights_(0, size()));");
+ buffer().decrease_indentation();
+ buffer().add_line("}");
+ buffer().add_line("else {");
+ buffer().increase_indentation();
+ buffer().add_line("memory::fill(weights_, 1.0);");
+ buffer().decrease_indentation();
+ buffer().add_line("}");
buffer().add_line();
}
@@ -378,6 +390,16 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
buffer().add_line("}");
buffer().add_line();
+ // Override `set_weights` method only for density mechanisms.
+ if (module_->kind() == moduleKind::density) {
+ buffer().add_line("void set_weights(array&& weights) override {");
+ buffer().increase_indentation();
+ buffer().add_line("memory::copy(weights, weights_(0, size()));");
+ buffer().decrease_indentation();
+ buffer().add_line("}");
+ buffer().add_line();
+ }
+
//////////////////////////////////////////////
// print ion channel interface
//////////////////////////////////////////////
@@ -525,6 +547,61 @@ CUDAPrinter::CUDAPrinter(Module &m, bool o)
}
}
+ std::unordered_set<std::string> scalar_set;
+ for (auto& v: scalar_variables) {
+ scalar_set.insert(v->name());
+ }
+
+ std::vector<Id> global_param_ids;
+ std::vector<Id> instance_param_ids;
+
+ for (const Id& id: module_->parameter_block().parameters) {
+ auto var = id.token.spelling;
+ (scalar_set.count(var)? global_param_ids: instance_param_ids).push_back(id);
+ }
+ const std::vector<Id>& state_ids = module_->state_block().state_variables;
+
+ if (!instance_param_ids.empty() || !state_ids.empty()) {
+ buffer().add_line("view base::* field_view_ptr(const char* id) const override {");
+ buffer().increase_indentation();
+ buffer().add_line("static const std::pair<const char*, view "+class_name+"::*> field_tbl[] = {");
+ buffer().increase_indentation();
+ for (const auto& id: instance_param_ids) {
+ auto var = id.token.spelling;
+ buffer().add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ for (const auto& id: state_ids) {
+ auto var = id.token.spelling;
+ buffer().add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ buffer().decrease_indentation();
+ buffer().add_line("};");
+ buffer().add_line();
+ buffer().add_line("auto* pptr = util::table_lookup(field_tbl, id);");
+ buffer().add_line("return pptr? static_cast<view base::*>(*pptr): nullptr;");
+ buffer().decrease_indentation();
+ buffer().add_line("}");
+ buffer().add_line();
+ }
+
+ if (!global_param_ids.empty()) {
+ buffer().add_line("value_type base::* field_value_ptr(const char* id) const override {");
+ buffer().increase_indentation();
+ buffer().add_line("static const std::pair<const char*, value_type "+class_name+"::*> field_tbl[] = {");
+ buffer().increase_indentation();
+ for (const auto& id: global_param_ids) {
+ auto var = id.token.spelling;
+ buffer().add_line("{\""+var+"\", &"+class_name+"::"+var+"},");
+ }
+ buffer().decrease_indentation();
+ buffer().add_line("};");
+ buffer().add_line();
+ buffer().add_line("auto* pptr = util::table_lookup(field_tbl, id);");
+ buffer().add_line("return pptr? static_cast<value_type base::*>(*pptr): nullptr;");
+ buffer().decrease_indentation();
+ buffer().add_line("}");
+ buffer().add_line();
+ }
//////////////////////////////////////////////
//////////////////////////////////////////////
diff --git a/modcc/parser.cpp b/modcc/parser.cpp
index 84f79413..4c55321b 100644
--- a/modcc/parser.cpp
+++ b/modcc/parser.cpp
@@ -487,19 +487,14 @@ void Parser::parse_parameter_block() {
// look for equality
if(token_.type==tok::eq) {
get_token(); // consume '='
- if(token_.type==tok::minus) {
- parm.value = "-";
- get_token();
- }
- if(token_.type != tok::integer && token_.type != tok::real) {
+ parm.value = value_literal();
+ if(status_ == lexerStatus::error) {
success = 0;
goto parm_exit;
}
- parm.value += token_.spelling; // store value as a string
- get_token();
}
- // get the parameters
+ // get the units
if(line==location_.line && token_.type == tok::lparen) {
parm.units = unit_description();
if(status_ == lexerStatus::error) {
@@ -508,6 +503,14 @@ void Parser::parse_parameter_block() {
}
}
+ // get the range
+ if(line==location_.line && token_.type == tok::lt) {
+ parm.range = range_description();
+ if(status_ == lexerStatus::error) {
+ success = 0;
+ goto parm_exit;
+ }
+ }
block.parameters.push_back(parm);
}
@@ -595,6 +598,26 @@ ass_exit:
return;
}
+// Parse a value (integral or real) with possible preceding unary minus,
+// and return as a string.
+std::string Parser::value_literal() {
+ std::string value;
+
+ if(token_.type==tok::minus) {
+ value = "-";
+ get_token();
+ }
+ if(token_.type != tok::integer && token_.type != tok::real) {
+ error(pprintf("numeric constant not an integer or real number '%'", token_));
+ return "";
+ }
+ else {
+ value += token_.spelling;
+ get_token();
+ return value;
+ }
+}
+
std::vector<Token> Parser::unit_description() {
static const tok legal_tokens[] = {tok::identifier, tok::divide, tok::real, tok::integer};
int startline = location_.line;
@@ -602,7 +625,7 @@ std::vector<Token> Parser::unit_description() {
// check that we start with a left parenthesis
if(token_.type != tok::lparen) {
- error(pprintf("unit description must start with a parenthesis '%'", tokens));
+ error(pprintf("unit description must start with a parenthesis '%'", token_));
goto unit_exit;
}
@@ -611,7 +634,7 @@ std::vector<Token> Parser::unit_description() {
while(token_.type != tok::rparen) {
// check for illegal tokens or a new line
if( !is_in(token_.type,legal_tokens) || startline < location_.line ) {
- error(pprintf("incorrect unit description '%'", tokens));
+ error(pprintf("incorrect unit description '%'", token_));
goto unit_exit;
}
@@ -626,6 +649,35 @@ unit_exit:
return tokens;
}
+std::pair<Token, Token> Parser::range_description() {
+ int startline = location_.line;
+ Token lb, ub;
+
+ if(token_.type != tok::lt) {
+ error(pprintf("range description must start with a left angle bracket '%'", token_));
+ return {};
+ }
+
+ get_token();
+ lb = token_;
+
+ if(token_.type != tok::comma) {
+ error(pprintf("range description must separate lower and upper bound with a comma '%'", token_));
+ return {};
+ }
+
+ get_token();
+ ub = token_;
+
+ if(token_.type != tok::gt) {
+ error(pprintf("range description must end with a right angle bracket '%'", token_));
+ return {};
+ }
+
+ get_token();
+ return {lb, ub};
+}
+
// Returns a prototype expression for a function or procedure call
// Takes an optional argument that allows the user to specify the
// name of the prototype, which is used for prototypes where the name
diff --git a/modcc/parser.hpp b/modcc/parser.hpp
index 10a33e59..43fb323b 100644
--- a/modcc/parser.hpp
+++ b/modcc/parser.hpp
@@ -2,6 +2,7 @@
#include <memory>
#include <string>
+#include <utility>
#include "expression.hpp"
#include "lexer.hpp"
@@ -60,6 +61,8 @@ private:
std::vector<Token> comma_separated_identifiers();
std::vector<Token> unit_description();
+ std::string value_literal();
+ std::pair<Token, Token> range_description();
/// build the identifier list
void add_variables_to_symbols();
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 44488b94..6712d047 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -11,7 +11,6 @@ set(BASE_SOURCES
hardware/power.cpp
model.cpp
morphology.cpp
- parameter_list.cpp
partition_load_balance.cpp
profiling/memory_meter.cpp
profiling/meter_manager.cpp
diff --git a/src/backends/fvm.hpp b/src/backends/fvm.hpp
index be7a1153..42ba619c 100644
--- a/src/backends/fvm.hpp
+++ b/src/backends/fvm.hpp
@@ -11,7 +11,7 @@ struct null_backend: public multicore::backend {
return false;
}
- static mechanism make_mechanism(
+ static mechanism_ptr make_mechanism(
const std::string&,
size_type,
const_iview,
diff --git a/src/backends/gpu/fvm.hpp b/src/backends/gpu/fvm.hpp
index ebdf9961..95b03c3f 100644
--- a/src/backends/gpu/fvm.hpp
+++ b/src/backends/gpu/fvm.hpp
@@ -65,12 +65,12 @@ struct backend {
// mechanism infrastructure
using ion_type = ion<backend>;
-
- using mechanism = mechanism_ptr<backend>;
-
using stimulus = gpu::stimulus<backend>;
- static mechanism make_mechanism(
+ using mechanism = arb::mechanism<backend>;
+ using mechanism_ptr = std::unique_ptr<mechanism>;
+
+ static mechanism_ptr make_mechanism(
const std::string& name,
size_type mech_id,
const_iview vec_ci,
@@ -137,11 +137,11 @@ struct backend {
}
private:
- using maker_type = mechanism (*)(size_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
+ using maker_type = mechanism_ptr (*)(size_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
template <template <typename> class Mech>
- static mechanism maker(size_type mech_id, const_iview vec_ci, const_view vec_t, const_view vec_t_to, const_view vec_dt, view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
+ static mechanism_ptr maker(size_type mech_id, const_iview vec_ci, const_view vec_t, const_view vec_t_to, const_view vec_dt, view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
return arb::make_mechanism<Mech<backend>>
(mech_id, vec_ci, vec_t, vec_t_to, vec_dt, vec_v, vec_i, std::move(weights), std::move(node_indices));
}
diff --git a/src/backends/multicore/fvm.hpp b/src/backends/multicore/fvm.hpp
index 801f0744..af90096b 100644
--- a/src/backends/multicore/fvm.hpp
+++ b/src/backends/multicore/fvm.hpp
@@ -57,12 +57,12 @@ struct backend {
// mechanism infrastructure
//
using ion_type = ion<backend>;
-
- using mechanism = mechanism_ptr<backend>;
-
using stimulus = multicore::stimulus<backend>;
- static mechanism make_mechanism(
+ using mechanism = arb::mechanism<backend>;
+ using mechanism_ptr = std::unique_ptr<mechanism>;
+
+ static mechanism_ptr make_mechanism(
const std::string& name,
size_type mech_id,
const_iview vec_ci,
@@ -75,7 +75,7 @@ struct backend {
throw std::out_of_range("no mechanism in database : " + name);
}
- return mech_map_.find(name)->second(mech_id, vec_ci, vec_t, vec_t_to, vec_dt, vec_v, vec_i, array(weights), iarray(node_indices));
+ return mech_map_.find(name)->second(mech_id, vec_ci, vec_t, vec_t_to, vec_dt, vec_v, vec_i, memory::make_const_view(weights), memory::make_const_view(node_indices));
}
static bool has_mechanism(const std::string& name) {
@@ -149,11 +149,11 @@ struct backend {
}
private:
- using maker_type = mechanism (*)(value_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
+ using maker_type = mechanism_ptr (*)(value_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
template <template <typename> class Mech>
- static mechanism maker(value_type mech_id, const_iview vec_ci, const_view vec_t, const_view vec_t_to, const_view vec_dt, view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
+ static mechanism_ptr maker(value_type mech_id, const_iview vec_ci, const_view vec_t, const_view vec_t_to, const_view vec_dt, view vec_v, view vec_i, array&& weights, iarray&& node_indices) {
return arb::make_mechanism<Mech<backend>>
(mech_id, vec_ci, vec_t, vec_t_to, vec_dt, vec_v, vec_i, std::move(weights), std::move(node_indices));
}
diff --git a/src/backends/multicore/stimulus.hpp b/src/backends/multicore/stimulus.hpp
index d43b3d7f..2318ff9f 100644
--- a/src/backends/multicore/stimulus.hpp
+++ b/src/backends/multicore/stimulus.hpp
@@ -38,9 +38,6 @@ public:
return 0;
}
- void set_params() override {
- }
-
std::string name() const override {
return "stimulus";
}
diff --git a/src/cell.hpp b/src/cell.hpp
index 7879ab6a..28149edf 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <map>
#include <mutex>
#include <stdexcept>
#include <thread>
@@ -29,6 +30,7 @@ int find_compartment_index(
compartment_model const& graph
);
+// Probe type for cell descriptions.
struct cell_probe_address {
enum probe_kind {
membrane_voltage, membrane_current
@@ -38,6 +40,20 @@ struct cell_probe_address {
probe_kind kind;
};
+// Global parameter type for cell descriptions.
+
+struct specialized_mechanism {
+ std::string mech_name; // underlying mechanism
+
+ // parameters specify global constants for the specialized mechanism
+ std::vector<std::pair<std::string, double>> parameters;
+};
+
+struct cell_global_properties {
+ // Mechanisms specialized by mechanism-global parameter settings.
+ std::map<std::string, specialized_mechanism> special_mechs;
+};
+
// used in constructor below
struct clone_cell_t {};
constexpr clone_cell_t clone_cell{};
@@ -52,7 +68,7 @@ public:
struct synapse_instance {
segment_location location;
- parameter_list mechanism;
+ mechanism_spec mechanism;
};
struct stimulus_instance {
@@ -158,7 +174,7 @@ public:
//////////////////
// synapses
//////////////////
- void add_synapse(segment_location loc, parameter_list p)
+ void add_synapse(segment_location loc, mechanism_spec p)
{
synapses_.push_back(synapse_instance{loc, std::move(p)});
}
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index c7cbe89d..0fd6b4ee 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -191,6 +191,7 @@ public:
/// mechanism type
using mechanism = typename backend::mechanism;
+ using mechanism_ptr = typename backend::mechanism_ptr;
/// stimulus type
using stimulus = typename backend::stimulus;
@@ -221,7 +222,7 @@ public:
std::size_t size() const { return matrix_.size(); }
/// return reference to in iterable container of the mechanisms
- std::vector<mechanism>& mechanisms() { return mechanisms_; }
+ std::vector<mechanism_ptr>& mechanisms() { return mechanisms_; }
/// return reference to list of ions
std::map<ionKind, ion_type>& ions() { return ions_; }
@@ -249,12 +250,12 @@ public:
}
/// Return reference to the mechanism that matches name.
- /// The reference is const, because it this information should not be
+ /// The reference is const, because this information should not be
/// modified by the caller, however it is needed for unit testing.
- util::optional<const mechanism&> find_mechanism(const std::string& name) const {
+ util::optional<const mechanism_ptr&> find_mechanism(const std::string& name) const {
auto it = std::find_if(
std::begin(mechanisms_), std::end(mechanisms_),
- [&name](const mechanism& m) {return m->name()==name;});
+ [&name](const mechanism_ptr& m) {return m->name()==name;});
return it==mechanisms_.end() ? util::nothing: util::just(*it);
}
@@ -376,7 +377,7 @@ private:
array voltage_;
/// the set of mechanisms present in the cell
- std::vector<mechanism> mechanisms_;
+ std::vector<mechanism_ptr> mechanisms_;
/// the ion species
std::map<ionKind, ion_type> ions_;
@@ -420,6 +421,52 @@ private:
std::vector<value_type>& tmp_cv_areas,
std::vector<value_type>& cv_capacitance
);
+
+ // TODO: This process should be simpler when we can deal with mechanism prototypes and have
+ // separate initialization.
+ //
+ // Create possibly-specialized mechanism and add to mechanism set.
+ // Weights are unset, and should be set specifically with mechanism::set_weights().
+ mechanism& make_mechanism(
+ const std::string& name,
+ const std::map<std::string, specialized_mechanism>& special_mechs,
+ const std::vector<size_type>& node_indices)
+ {
+ std::string impl_name = name;
+ std::vector<std::pair<std::string, double>> global_params;
+
+ if (special_mechs.count(name)) {
+ const auto& spec_mech = special_mechs.at(name);
+ impl_name = spec_mech.mech_name;
+ global_params = spec_mech.parameters;
+ }
+
+ size_type mech_id = mechanisms_.size();
+ auto m = backend::make_mechanism(impl_name, mech_id, cv_to_cell_, time_, time_to_, dt_comp_, voltage_, current_, {}, node_indices);
+ if (impl_name!=name) {
+ m->set_alias(name);
+ }
+
+ for (const auto& pv: global_params) {
+ auto field = m->field_value_ptr(pv.first);
+ if (!field) {
+ throw std::invalid_argument("no scalar parameter "+pv.first+" in mechanism "+m->name());
+ }
+ m.get()->*field = pv.second;
+ }
+
+ mechanisms_.push_back(std::move(m));
+ return *mechanisms_.back();
+ }
+
+ // Throwing-wrapper around mechanism (range) parameter look up.
+ static view mech_field(mechanism& m, const std::string& param_name) {
+ auto p = m.field_view_ptr(param_name);
+ if (!p) {
+ throw std::invalid_argument("no parameter "+param_name+" in mechanism "+m.name());
+ }
+ return m.*p;
+ }
};
////////////////////////////////////////////////////////////////////////////////
@@ -442,6 +489,9 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
segment_cv_range cv_range;
+ auto cm = seg->cm;
+ auto rL = seg->rL;
+
if (auto soma = seg->as_soma()) {
// confirm assumption that there is one compartment in soma
if (ncomp!=1) {
@@ -449,10 +499,9 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
}
auto i = comp_ival.first;
auto area = math::area_sphere(soma->radius());
- auto c_m = soma->mechanism("membrane").get("c_m").value;
tmp_cv_areas[i] += area;
- cv_capacitance[i] += area*c_m;
+ cv_capacitance[i] += area*cm;
cv_range.segment_cvs = {comp_ival.first, comp_ival.first+1};
cv_range.areas = {0.0, area};
@@ -488,9 +537,6 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
// The first control volume of any cell corresponds to the soma
// and the first half of the first cable compartment of that cell.
- auto c_m = cable->mechanism("membrane").get("c_m").value;
- auto r_L = cable->mechanism("membrane").get("r_L").value;
-
auto divs = div_compartments<div_compartment_integrator>(cable, ncomp);
// assume that this segment has a parent, which is the case so long
@@ -535,7 +581,7 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
auto V2 = div.right.volume;
auto h = h1+h2;
- auto conductance = 1/r_L*h*V1*V2/(h2*h2*V1+h1*h1*V2);
+ auto conductance = 1/rL*h*V1*V2/(h2*h2*V1+h1*h1*V2);
// the scaling factor of 10^2 is to convert the quantity
// to micro Siemens [μS]
face_conductance[i] = 1e2 * conductance / h;
@@ -545,8 +591,8 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
tmp_cv_areas[j] += al;
tmp_cv_areas[i] += ar;
- cv_capacitance[j] += al * c_m;
- cv_capacitance[i] += ar * c_m;
+ cv_capacitance[j] += al * cm;
+ cv_capacitance[i] += ar * cm;
}
}
else {
@@ -576,6 +622,14 @@ void fvm_multicell<Backend>::initialize(
ncell_ = size(gids);
std::size_t targets_count = 0u;
+ // Handle any global parameters for these cell groups.
+ // (Currently: just specialized mechanisms).
+ std::map<std::string, specialized_mechanism> special_mechs;
+ util::any gprops = rec.get_global_properties(cable1d_neuron);
+ if (gprops.has_value()) {
+ special_mechs = util::any_cast<cell_global_properties&>(gprops).special_mechs;
+ }
+
// Take cell descriptions from recipe. These are used initially
// to count compartments for allocation of data structures, and
// then interrogated again for the details for each cell in turn.
@@ -611,15 +665,36 @@ void fvm_multicell<Backend>::initialize(
}
memory::copy(cv_to_cell_tmp, cv_to_cell_);
- // look up table: mechanism name -> list of cv_range objects
- std::map<std::string, std::vector<segment_cv_range>> mech_to_cv_range;
+ // TODO: mechanism parameters are currently indexed by string keys; more efficient
+ // to use the mechanism member pointers, when these become easily accessible via
+ // the mechanism catalogue interface.
+
+ // look up table: mechanism name -> list of cv_range objects and parameter settings.
+
+ struct mech_area_contrib {
+ size_type index;
+ value_type area;
+ };
+
+ struct mech_info {
+ segment_cv_range cv_range;
+ // Note: owing to linearity constraints, the only parameters for which it is
+ // sensible to modify are those which make a linear contribution to currents
+ // (or ion fluxes, etc.)
+ std::map<std::string, value_type> param_map;
+ std::vector<mech_area_contrib> contributions;
+ };
+
+ std::map<std::string, std::vector<mech_info>> mech_map;
- // look up table: point mechanism (synapse) name -> CV indices and target numbers.
- struct syn_cv_and_target {
+ // look up table: point mechanism (synapse) name -> list of CV indices, target numbers, parameters.
+ struct syn_info {
cell_lid_type cv;
cell_lid_type target;
+ std::map<std::string, value_type> param_map;
};
- std::map<std::string, std::vector<syn_cv_and_target>> syn_mech_map;
+
+ std::map<std::string, std::vector<syn_info>> syn_mech_map;
// initialize vector used for matrix creation.
std::vector<size_type> group_parent_index(ncomp);
@@ -631,9 +706,9 @@ void fvm_multicell<Backend>::initialize(
// Create each cell:
// Allocate scratch storage for calculating quantities used to build the
- // linear system: these will later be copied into target-specific storage
+ // linear system: these will later be copied into target-specific storage.
- // face_conductance_[i] = area_face / (r_L * delta_x);
+ // face_conductance_[i] = area_face / (rL * delta_x);
std::vector<value_type> face_conductance(ncomp); // [µS]
/// cv_capacitance_[i] is the capacitance of CV membrane
std::vector<value_type> cv_capacitance(ncomp); // [µm^2*F*m^-2 = pF]
@@ -681,26 +756,24 @@ void fvm_multicell<Backend>::initialize(
face_conductance, tmp_cv_areas, cv_capacitance);
for (const auto& mech: seg->mechanisms()) {
- if (mech.name()!="membrane") {
- mech_to_cv_range[mech.name()].push_back(cv_range);
- }
+ mech_map[mech.name()].push_back({cv_range, mech.values()});
}
}
for (const auto& syn: c.synapses()) {
const auto& name = syn.mechanism.name();
- auto& map_entry = syn_mech_map[name];
cell_lid_type syn_cv = comp_ival.first + find_cv_index(syn.location, graph);
cell_lid_type target_index = targets_count++;
- map_entry.push_back(syn_cv_and_target{syn_cv, target_index});
+ syn_mech_map[name].push_back({syn_cv, target_index, syn.mechanism.values()});
}
//
// add the stimuli
//
+ // TODO: use same process as for synapses!
// step 1: pack the index and parameter information into flat vectors
std::vector<size_type> stim_index;
std::vector<value_type> stim_durations;
@@ -726,7 +799,7 @@ void fvm_multicell<Backend>::initialize(
cv_to_cell_, time_, time_to_, dt_comp_,
voltage_, current_, memory::make_const_view(stim_index));
stim->set_parameters(stim_amplitudes, stim_durations, stim_delays);
- mechanisms_.push_back(mechanism(stim));
+ mechanisms_.push_back(mechanism_ptr(stim));
}
// calculate spike detector handles are their corresponding compartment indices
@@ -770,85 +843,177 @@ void fvm_multicell<Backend>::initialize(
matrix_ = matrix_type(
group_parent_index, cell_comp_bounds, cv_capacitance, face_conductance);
- // For each density mechanism build the full node index, i.e the list of
- // compartments with that mechanism, then build the mechanism instance.
- std::vector<size_type> mech_cv_index(ncomp);
- std::vector<value_type> mech_cv_weight(ncomp);
+ // Keep cv index list for each mechanism for ion set up below.
std::map<std::string, std::vector<size_type>> mech_to_cv_index;
- for (auto const& entry: mech_to_cv_range) {
+
+ // Working vectors (re-used per mechanism).
+ std::vector<size_type> mech_cv(ncomp);
+ std::vector<value_type> mech_weight(ncomp);
+
+ for (auto& entry: mech_map) {
const auto& mech_name = entry.first;
- const auto& seg_cv_ranges = entry.second;
+ auto& segments = entry.second;
+
+ mech_cv.clear();
+ mech_weight.clear();
+
+ // Three passes are performed over the segment list:
+ // 1. Compute the CVs and area contributions where the mechanism is instanced.
+ // 2. Build table of modified parameters together with default values.
+ // 3. Compute weights and parameters.
+ // The mechanism is instantiated after the first pass, in order to gain
+ // access to default mechanism parameter values.
- // Clear the pre-allocated storage for mechanism indexes and weights.
- // Reuse the same vectors each time to have only one malloc and free
- // outside of the loop for each
- mech_cv_index.clear();
- mech_cv_weight.clear();
+ for (auto& seg: segments) {
+ const auto& rng = seg.cv_range;
+ seg.contributions.reserve(rng.size());
- for (auto& rng: seg_cv_ranges) {
if (rng.has_parent()) {
- // locate the parent cv in the partially constructed list of cv indexes
- auto it = algorithms::binary_find(mech_cv_index, rng.parent_cv);
- if (it == mech_cv_index.end()) {
- mech_cv_index.push_back(rng.parent_cv);
- mech_cv_weight.push_back(0);
+ auto cv = rng.parent_cv;
+
+ auto it = algorithms::binary_find(mech_cv, cv);
+ size_type pos = it - mech_cv.begin();
+
+ if (it == mech_cv.end()) {
+ mech_cv.push_back(cv);
}
- auto pos = std::distance(std::begin(mech_cv_index), it);
- // add area contribution to the parent cv for the segment
- mech_cv_weight[pos] += rng.areas.first;
+ seg.contributions.push_back({pos, rng.areas.first});
}
- util::append(mech_cv_index, make_span(rng.segment_cvs));
- util::append(mech_cv_weight, subrange_view(tmp_cv_areas, rng.segment_cvs));
- // adjust the last CV
- mech_cv_weight.back() = rng.areas.second;
+ for (auto cv: make_span(rng.segment_cvs)) {
+ size_type pos = mech_cv.size();
+ mech_cv.push_back(cv);
+ seg.contributions.push_back({pos, tmp_cv_areas[cv]});
+ }
- EXPECTS(mech_cv_weight.size()==mech_cv_index.size());
+ // Last CV contribution may be only partial, so adjust.
+ seg.contributions.back().area = rng.areas.second;
+ }
+
+ auto nindex = mech_cv.size();
+
+ EXPECTS(std::is_sorted(mech_cv.begin(), mech_cv.end()));
+ EXPECTS(nindex>0);
+
+ auto& mech = make_mechanism(mech_name, special_mechs, mech_cv);
+
+ // Save the indices for ion set up below.
+
+ mech_to_cv_index[mech_name] = mech_cv;
+
+ // Build modified (non-global) parameter table.
+
+ struct param_tbl_entry {
+ std::vector<value_type> values; // staged for writing to mechanism
+ view data; // view to corresponding data in mechanism
+ value_type dflt; // default value for parameter
+ };
+
+ std::map<std::string, param_tbl_entry> param_tbl;
+
+ for (const auto& seg: segments) {
+ for (const auto& pv: seg.param_map) {
+ if (param_tbl.count(pv.first)) {
+ continue;
+ }
+
+ // Grab default value from mechanism data.
+ auto& entry = param_tbl[pv.first];
+ entry.data = mech_field(mech, pv.first);
+ entry.dflt = entry.data[0];
+ entry.values.assign(nindex, 0.);
+ }
+ }
+
+ // Perform another pass of segment list to compute weights and (non-global) parameters.
+
+ mech_weight.assign(nindex, 0.);
+
+ for (const auto& seg: segments) {
+ for (auto cw: seg.contributions) {
+ mech_weight[cw.index] += cw.area;
+
+ for (auto& entry: param_tbl) {
+ value_type v = entry.second.dflt;
+ const auto& name = entry.first;
+
+ auto it = seg.param_map.find(name);
+ if (it != seg.param_map.end()) {
+ v = it->second;
+ }
+
+ entry.second.values[cw.index] += cw.area*v;
+ }
+ }
+ }
+
+ for (auto& entry: param_tbl) {
+ for (size_type i = 0; i<nindex; ++i) {
+ entry.second.values[i] /= mech_weight[i];
+ }
+ memory::copy(entry.second.values, entry.second.data);
}
// Scale the weights to get correct units (see w_i^d in formulation docs)
// The units for the density channel weights are [10^2 μm^2 = 10^-10 m^2],
// which requires that we scale the areas [μm^2] by 10^-2
- for (auto& w: mech_cv_weight) {
+
+ for (auto& w: mech_weight) {
w *= 1e-2;
}
-
- size_type mech_id = mechanisms_.size();
- mechanisms_.push_back(
- backend::make_mechanism(mech_name, mech_id, cv_to_cell_, time_, time_to_, dt_comp_, voltage_, current_, mech_cv_weight, mech_cv_index));
-
- // Save the indices for ion set up below.
- mech_to_cv_index[mech_name] = mech_cv_index;
+ mech.set_weights(memory::make_const_view(mech_weight));
}
target_handles.resize(targets_count);
// Create point (synapse) mechanisms.
- for (auto& syni: syn_mech_map) {
+ for (auto& map_entry: syn_mech_map) {
size_type mech_id = mechanisms_.size();
- const auto& mech_name = syni.first;
- auto& cv_assoc = syni.second;
+ const auto& mech_name = map_entry.first;
+ auto& syn_data = map_entry.second;
+ auto n_instance = syn_data.size();
- // Sort CV indices but keep track of their corresponding targets.
- auto cv_index = [](syn_cv_and_target x) { return x.cv; };
- util::stable_sort_by(cv_assoc, cv_index);
- std::vector<cell_lid_type> cv_indices = assign_from(transform_view(cv_assoc, cv_index));
+ // Build permutation p such that p[j] is the index into
+ // syn_data for the jth synapse of this mechanism type as ordered by cv index.
- // Create the mechanism.
- // An empty weight vector is supplied, because there are no weights applied to point
- // processes, because their currents are calculated with the target units of [nA]
- mechanisms_.push_back(
- backend::make_mechanism(mech_name, mech_id, cv_to_cell_, time_, time_to_, dt_comp_, voltage_, current_, {}, cv_indices));
+ auto cv_of = [&](cell_lid_type i) { return syn_data[i].cv; };
+
+ std::vector<cell_lid_type> p(n_instance);
+ std::iota(p.begin(), p.end(), 0u);
+ util::sort_by(p, cv_of);
+
+ std::vector<cell_lid_type> mech_cv;
+ mech_cv.reserve(n_instance);
+
+ // Build mechanism cv index vector and targets.
+ for (auto i: make_span(0u, n_instance)) {
+ const auto& syn = syn_data[p[i]];
+ mech_cv.push_back(syn.cv);
+ target_handles[syn.target] = target_handle(mech_id, i, cv_to_cell_tmp[syn.cv]);
+ }
+
+ auto& mech = make_mechanism(mech_name, special_mechs, mech_cv);
// Save the indices for ion set up below.
- mech_to_cv_index[mech_name] = cv_indices;
+ mech_to_cv_index[mech_name] = mech_cv;
- // Make the target handles.
- cell_lid_type instance = 0;
- for (auto entry: cv_assoc) {
- target_handles[entry.target] = target_handle(mech_id, instance++, cv_to_cell_tmp[entry.cv]);
+ // Update the mechanism parameters.
+ std::map<std::string, std::vector<std::pair<cell_lid_type, fvm_value_type>>> param_assigns;
+ for (auto i: make_span(0u, n_instance)) {
+ for (const auto& pv: syn_data[p[i]].param_map) {
+ param_assigns[pv.first].push_back({i, pv.second});
+ }
+ }
+
+ for (const auto& pa: param_assigns) {
+ view field_data = mech_field(mech, pa.first);
+ host_array field_values = field_data;
+ for (const auto &iv: pa.second) {
+ field_values[iv.first] = iv.second;
+ }
+ memory::copy(field_values, field_data);
}
}
@@ -866,7 +1031,7 @@ void fvm_multicell<Backend>::initialize(
std::vector<size_type> indexes(index_set.begin(), index_set.end());
// create the ion state
- if(indexes.size()) {
+ if (indexes.size()) {
ions_[ion] = indexes;
}
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index 2711c679..33d53872 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -4,16 +4,44 @@
#include <memory>
#include <string>
+#include <backends/fvm_types.hpp>
#include <backends/event.hpp>
#include <backends/multi_event_stream_state.hpp>
#include <ion.hpp>
-#include <parameter_list.hpp>
#include <util/indirect.hpp>
#include <util/meta.hpp>
#include <util/make_unique.hpp>
namespace arb {
+struct field_spec {
+ enum field_kind {
+ parameter, // defined in 'PARAMETER' block and a 'RANGE' variable.
+ global, // defined in 'PARAMETER' block and a 'GLOBAL' variable.
+ state, // defined in 'STATE' block; run-time, read only values.
+ };
+ enum field_kind kind = parameter;
+
+ std::string units;
+
+ fvm_value_type default_value = 0;
+ fvm_value_type lower_bound = std::numeric_limits<fvm_value_type>::lowest();
+ fvm_value_type upper_bound = std::numeric_limits<fvm_value_type>::max();
+
+ // Until C++14, we need a ctor to provide default values instead of using
+ // default member initializers and aggregate initialization.
+ field_spec(
+ enum field_kind kind = parameter,
+ std::string units = "",
+ fvm_value_type default_value = 0.,
+ fvm_value_type lower_bound = std::numeric_limits<fvm_value_type>::lowest(),
+ fvm_value_type upper_bound = std::numeric_limits<fvm_value_type>::max()
+ ):
+ kind(kind), units(units), default_value(default_value), lower_bound(lower_bound), upper_bound(upper_bound)
+ {}
+};
+
+
enum class mechanismKind {point, density};
/// The mechanism type is templated on a memory policy type.
@@ -60,7 +88,12 @@ public:
return node_index_;
}
- virtual void set_params() = 0;
+ // Save pointers to data for use with GPU-side mechanisms;
+ // TODO: might be able to remove this method if we separate instantiation
+ // from initialization.
+ virtual void set_params() {}
+ virtual void set_weights(array&& weights) {} // override for density mechanisms
+
virtual std::string name() const = 0;
virtual std::size_t memory() const = 0;
virtual void nrn_init() = 0;
@@ -71,6 +104,24 @@ public:
virtual void set_ion(ionKind k, ion_type& i, const std::vector<size_type>& index) = 0;
virtual mechanismKind kind() const = 0;
+ // Mechanism instances with different global parameter settings can be distinguished by alias.
+ std::string alias() const {
+ return alias_.empty()? name(): alias_;
+ }
+
+ void set_alias(std::string alias) {
+ alias_ = std::move(alias);
+ }
+
+ // For non-global fields:
+ virtual view mechanism::* field_view_ptr(const char* id) const { return nullptr; }
+ // For global fields:
+ virtual value_type mechanism::* field_value_ptr(const char* id) const { return nullptr; }
+
+ // Convenience wrappers for field access methods with string parameter.
+ view mechanism::* field_view_ptr(const std::string& id) const { return field_view_ptr(id.c_str()); }
+ value_type mechanism::* field_value_ptr(const std::string& id) const { return field_value_ptr(id.c_str()); }
+
// net_receive() is used internally by deliver_events(), but
// is exposed primarily for unit testing.
virtual void net_receive(int, value_type) {};
@@ -94,6 +145,8 @@ public:
view vec_i_;
// Maps mechanism instance index to compartment index.
iarray node_index_;
+
+ std::string alias_;
};
template <class Backend>
diff --git a/src/mechinfo.hpp b/src/mechinfo.hpp
new file mode 100644
index 00000000..cd40c387
--- /dev/null
+++ b/src/mechinfo.hpp
@@ -0,0 +1,108 @@
+#pragma once
+
+/* Mechanism schema classes, catalogue and parameter specification.
+ *
+ * Catalogue and schemata have placeholder implementations, to be
+ * completed in future work.
+ *
+ * The `mechanism_spec` class is the public interface for describing
+ * a mechanism and its (non-global) parameters in cable1d_cell
+ * recipes. It presents a map-like interface for accessing and querying
+ * parameter values, and parameter assignments will be validated
+ * against a corresponding schema when that infrastructure is in
+ * place.
+ */
+
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace arb {
+
+struct mechanism_schema_field {
+ void validate(double) const {}
+ double default_value = 0;
+};
+
+struct mechanism_schema {
+ const mechanism_schema_field* field(const std::string& key) const {
+ static mechanism_schema_field dummy_field;
+ return &dummy_field;
+ }
+
+ static const mechanism_schema* dummy_schema() {
+ static mechanism_schema d;
+ return &d;
+ }
+};
+
+class mechanism_spec {
+public:
+ struct field_proxy {
+ mechanism_spec* m;
+ std::string key;
+
+ field_proxy& operator=(double v) {
+ m->set(key, v);
+ return *this;
+ }
+
+ operator double() const {
+ return m->get(key);
+ }
+ };
+
+ // implicit
+ mechanism_spec(std::string name): name_(std::move(name)) {
+ // get schema pointer from global catalogue, or throw
+ schema_ = mechanism_schema::dummy_schema();
+ if (!schema_) {
+ throw std::runtime_error("no mechanism "+name_);
+ }
+ }
+
+ // implicit
+ mechanism_spec(const char* name): mechanism_spec(std::string(name)) {}
+
+ mechanism_spec& set(std::string key, double value) {
+ auto field_schema = schema_->field(key);
+ if (!field_schema) {
+ throw std::runtime_error("no field "+key+" in mechanism "+name_);
+ }
+
+ field_schema->validate(value);
+ param_[key] = value;
+ return *this;
+ }
+
+ double operator[](const std::string& key) const {
+ return get(key);
+ }
+
+ double get(const std::string& key) const {
+ auto field_schema = schema_->field(key);
+ if (!field_schema) {
+ throw std::runtime_error("no field "+key+" in mechanism "+name_);
+ }
+
+ auto it = param_.find(key);
+ return it==param_.end()? field_schema->default_value: it->second;
+ }
+
+ field_proxy operator[](const std::string& key) {
+ return {this, key};
+ }
+
+ const std::map<std::string, double>& values() const {
+ return param_;
+ }
+
+ const std::string& name() const { return name_; }
+
+private:
+ std::string name_;
+ std::map<std::string, double> param_;
+ const mechanism_schema* schema_; // non-owning; schema must have longer lifetime
+};
+
+} // namespace arb
diff --git a/src/memory/array.hpp b/src/memory/array.hpp
index cdc97fe4..c5e358e6 100644
--- a/src/memory/array.hpp
+++ b/src/memory/array.hpp
@@ -223,25 +223,6 @@ public:
std::copy(b, e, this->begin());
}
- template <typename Seq>
- array(
- const Seq& seq,
- arb::util::enable_if_t<
- !std::is_convertible<Seq, std::size_t>::value
- && !impl::is_array_t<Seq>::value >* = nullptr
- ):
- base(coordinator_type().allocate(arb::util::size(seq)))
- {
-#ifdef VERBOSE
- std::cerr << util::green("array(iterator, iterator)")
- << " " << util::type_printer<array>::print()
- << "\n this " << util::pretty_printer<array>::print(*this) << "\n";
- //<< "\n other " << util::pretty_printer<Other>::print(other) << std::endl;
-#endif
- auto canon = arb::util::canonical_view(seq);
- std::copy(std::begin(canon), std::end(canon), this->begin());
- }
-
// use the accessors provided by array_view
// this enforces the requirement that accessing all of or a sub-array of an
// array should return a view, not a new array.
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
deleted file mode 100644
index 5b5d974c..00000000
--- a/src/parameter_list.cpp
+++ /dev/null
@@ -1,108 +0,0 @@
-#include <algorithm>
-#include <iostream>
-
-#include "parameter_list.hpp"
-
-namespace arb {
-
-bool parameter_list::add_parameter(parameter p) {
- if (has_parameter(p.name)) {
- return false;
- }
-
- parameters_.push_back(std::move(p));
-
- return true;
-}
-
-bool parameter_list::has_parameter(std::string const& n) const {
- return find_by_name(n) != parameters_.end();
-}
-
-int parameter_list::num_parameters() const {
- return parameters_.size();
-}
-
-// returns true if parameter was succesfully updated
-// returns false if parameter was not updated, i.e. if
-// - no parameter with name n exists
-// - value is not in the valid range
-bool parameter_list::set(std::string const& n, parameter_list::value_type v) {
- auto p = find_by_name(n);
- if(p!=parameters_.end()) {
- if(p->is_in_range(v)) {
- p->value = v;
- return true;
- }
- }
- return false;
-}
-
-parameter& parameter_list::get(std::string const& n) {
- auto it = find_by_name(n);
- if (it==parameters_.end()) {
- throw std::domain_error(
- "parameter list does not contain parameter"
- );
- }
- return *it;
-}
-
-const parameter& parameter_list::get(std::string const& n) const {
- auto it = find_by_name(n);
- if (it==parameters_.end()) {
- throw std::domain_error(
- "parameter list does not contain parameter"
- );
- }
- return *it;
-}
-
-std::string const& parameter_list::name() const {
- return mechanism_name_;
-}
-
-std::vector<parameter> const& parameter_list::parameters() const {
- return parameters_;
-}
-
-auto parameter_list::find_by_name(std::string const& n)
- -> decltype(parameters_.begin())
-{
- return
- std::find_if(
- parameters_.begin(), parameters_.end(),
- [&n](parameter const& p) {return p.name == n;}
- );
-}
-
-auto parameter_list::find_by_name(std::string const& n) const
- -> decltype(parameters_.begin())
-{
- return
- std::find_if(
- parameters_.begin(), parameters_.end(),
- [&n](parameter const& p) {return p.name == n;}
- );
-}
-
-} // namespace arb
-
-std::ostream&
-operator<<(std::ostream& o, arb::parameter const& p) {
- return o
- << "parameter("
- << "name " << p.name
- << " : value " << p.value
- << " : range " << p.range
- << ")";
-}
-
-std::ostream&
-operator<<(std::ostream& o, arb::parameter_list const& l) {
- o << "parameters \"" << l.name() << "\" :\n";
- for(arb::parameter const& p : l.parameters()) {
- o << " " << p << "\n";
- }
- return o;
-}
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
deleted file mode 100644
index bc47ef0a..00000000
--- a/src/parameter_list.hpp
+++ /dev/null
@@ -1,229 +0,0 @@
-#pragma once
-
-#include <limits>
-#include <ostream>
-#include <stdexcept>
-#include <string>
-#include <vector>
-
-namespace arb {
-
- template <typename T>
- struct value_range {
- using value_type = T;
-
- value_range()
- : min(std::numeric_limits<value_type>::quiet_NaN()),
- max(std::numeric_limits<value_type>::quiet_NaN())
- { }
-
- value_range(value_type const& left, value_type const& right)
- : min(left),
- max(right)
- {
- if(left>right) {
- throw std::out_of_range(
- "parameter range must have left <= right"
- );
- }
- }
-
- bool has_lower_bound() const
- {
- return min==min;
- }
-
- bool has_upper_bound() const
- {
- return max==max;
- }
-
- bool is_in_range(value_type const& v) const
- {
- if(has_lower_bound()) {
- if(v<min) return false;
- }
- if(has_upper_bound()) {
- if(v>max) return false;
- }
- return true;
- }
-
- value_type min;
- value_type max;
- };
-
- struct parameter {
- using value_type = double;
- using range_type = value_range<value_type>;
-
- parameter() = delete;
-
- parameter(std::string n, value_type v, range_type r=range_type())
- : name(std::move(n)),
- value(v),
- range(r)
- {
- if(!is_in_range(v)) {
- throw std::out_of_range(
- "parameter value is out of permitted value range"
- );
- }
- }
-
- bool is_in_range(value_type v) const
- {
- return range.is_in_range(v);
- }
-
- std::string name;
- value_type value;
- range_type range;
- };
-
- // Use a dumb container class for now
- // might have to use a more sophisticated interface in the future if need be
- class parameter_list {
- public :
-
- using value_type = double;
-
- parameter_list(std::string n)
- : mechanism_name_(std::move(n))
- { }
-
- bool has_parameter(std::string const& n) const;
- bool add_parameter(parameter p);
-
- // returns true if parameter was succesfully updated
- // returns false if parameter was not updated, i.e. if
- // - no parameter with name n exists
- // - value is not in the valid range
- bool set(std::string const& n, value_type v);
- parameter& get(std::string const& n);
- const parameter& get(std::string const& n) const;
-
- std::string const& name() const;
-
- std::vector<parameter> const& parameters() const;
-
- int num_parameters() const;
-
- private:
-
- std::vector<parameter> parameters_;
- std::string mechanism_name_;
-
- // need two versions for const and non-const applications
- auto find_by_name(std::string const& n)
- -> decltype(parameters_.begin());
-
- auto find_by_name(std::string const& n) const
- -> decltype(parameters_.begin());
-
- };
-
- ///////////////////////////////////////////////////////////////////////////
- // predefined parameter sets
- ///////////////////////////////////////////////////////////////////////////
-
- /// default set of parameters for the cell membrane that are added to every
- /// segment when it is created.
- class membrane_parameters
- : public parameter_list
- {
- public:
-
- using base = parameter_list;
-
- using base::value_type;
-
- using base::set;
- using base::get;
- using base::parameters;
- using base::has_parameter;
-
- membrane_parameters()
- : base("membrane")
- {
- // from nrn/src/nrnoc/membdef.h
- //#define DEF_cm 1. // uF/cm^2
- //#define DEF_Ra 35.4 // ohm-cm
- //#define DEF_celsius 6.3 // deg-C
- //#define DEF_vrest -65. // mV
- // r_L is called Ra in Neuron
- //base::add_parameter({"c_m", 10e-6, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
- base::add_parameter({"c_m", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
- base::add_parameter({"r_L", 100.00, {0., 1e9}}); // equivalent to Ra in Neuron : Ohm.cm
- }
- };
-
- /// parameters for the classic Hodgkin & Huxley model (1952)
- class hh_parameters
- : public parameter_list
- {
- public:
-
- using base = parameter_list;
-
- using base::value_type;
-
- using base::set;
- using base::get;
- using base::parameters;
- using base::has_parameter;
-
- hh_parameters()
- : base("hh")
- {
- base::add_parameter({"gnabar", 0.12, {0, 1e9}});
- base::add_parameter({"gkbar", 0.036, {0, 1e9}});
- base::add_parameter({"gl", 0.0003,{0, 1e9}});
- base::add_parameter({"el", -54.3});
- }
- };
-
- /// parameters for passive channel
- class pas_parameters
- : public parameter_list
- {
- public:
-
- using base = parameter_list;
-
- using base::value_type;
-
- using base::set;
- using base::get;
- using base::parameters;
- using base::has_parameter;
-
- pas_parameters()
- : base("pas")
- {
- base::add_parameter({"g", 0.001, {0, 1e9}});
- base::add_parameter({"e", -70});
- }
- };
-
-} // namespace arb
-
-template <typename T>
-std::ostream& operator<<(std::ostream& o, arb::value_range<T> const& r)
-{
- o << "[";
- if(r.has_lower_bound())
- o << r.min;
- else
- o<< "-inf";
- o << ", ";
- if(r.has_upper_bound())
- o << r.max;
- else
- o<< "inf";
- return o << "]";
-}
-
-std::ostream& operator<<(std::ostream& o, arb::parameter const& p);
-std::ostream& operator<<(std::ostream& o, arb::parameter_list const& l);
-
diff --git a/src/recipe.hpp b/src/recipe.hpp
index 52c8c377..106fac4b 100644
--- a/src/recipe.hpp
+++ b/src/recipe.hpp
@@ -14,6 +14,8 @@ namespace arb {
struct probe_info {
cell_member_type id;
probe_tag tag;
+
+ // Address type will be specific to cell kind of cell `id.gid`.
util::any address;
};
@@ -53,7 +55,8 @@ class recipe {
public:
virtual cell_size_type num_cells() const = 0;
- virtual util::unique_any get_cell_description(cell_gid_type) const = 0;
+ // Cell description type will be specific to cell kind of cell with given gid.
+ virtual util::unique_any get_cell_description(cell_gid_type gid) const = 0;
virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
virtual cell_size_type num_sources(cell_gid_type) const = 0;
@@ -62,6 +65,9 @@ public:
virtual std::vector<cell_connection> connections_on(cell_gid_type) const = 0;
virtual probe_info get_probe(cell_member_type probe_id) const = 0;
+
+ // Global property type will be specific to given cell kind.
+ virtual util::any get_global_properties(cell_kind) const { return util::any{}; };
};
} // namespace arb
diff --git a/src/segment.hpp b/src/segment.hpp
index bfdfa1f0..66534669 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -8,19 +8,12 @@
#include <compartment.hpp>
#include <math.hpp>
#include <morphology.hpp>
-#include <parameter_list.hpp>
+#include <mechinfo.hpp>
#include <point.hpp>
#include <util/make_unique.hpp>
namespace arb {
-template <typename T,
- typename valid = typename std::is_floating_point<T>::type>
-struct segment_properties {
- T rL = 180.0; // resistivity [Ohm.cm]
- T cm = 0.01; // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
-};
-
// forward declarations of segment specializations
class soma_segment;
class cable_segment;
@@ -87,79 +80,44 @@ public:
return false;
}
- segment_properties<value_type> properties;
-
- void add_mechanism(parameter_list p)
- {
- auto it = std::find_if(
- mechanisms_.begin(), mechanisms_.end(),
- [&p](parameter_list const& l){return l.name() == p.name();}
- );
- if(it!=mechanisms_.end()) {
- throw std::out_of_range(
- "attempt to add a mechanism parameter set to a segment which has an existing mechanism with the same name"
- );
- }
-
- mechanisms_.push_back(std::move(p));
+ util::optional<mechanism_spec&> mechanism(const std::string& name) {
+ auto it = std::find_if(mechanisms_.begin(), mechanisms_.end(),
+ [&](mechanism_spec& m) { return m.name()==name; });
+ return it==mechanisms_.end()? util::nothing: util::just(*it);
}
- parameter_list& mechanism(std::string const& n)
- {
- auto it = std::find_if(
- mechanisms_.begin(), mechanisms_.end(),
- [&n](parameter_list const& l){return l.name() == n;}
- );
- if(it==mechanisms_.end()) {
- throw std::out_of_range(
- "attempt to access a parameter that is not defined in a segment"
- );
+ void add_mechanism(mechanism_spec mech) {
+ auto m = mechanism(mech.name());
+ if (m) {
+ m.get() = std::move(mech);
}
-
- return *it;
- }
-
- const parameter_list& mechanism(std::string const& n) const
- {
- auto it = std::find_if(
- mechanisms_.begin(), mechanisms_.end(),
- [&n](parameter_list const& l){return l.name() == n;}
- );
- if(it==mechanisms_.end()) {
- throw std::out_of_range(
- "attempt to access a parameter that is not defined in a segment"
- );
+ else {
+ mechanisms_.push_back(std::move(mech));
}
-
- return *it;
}
- std::vector<parameter_list>& mechanisms()
- {
+ const std::vector<mechanism_spec>& mechanisms() {
return mechanisms_;
}
- const std::vector<parameter_list>& mechanisms() const
- {
+ const std::vector<mechanism_spec>& mechanisms() const {
return mechanisms_;
}
+ // common electrical properties
+ value_type rL = 100.0; // resistivity [Ohm.cm]
+ value_type cm = 0.01; // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
+
protected:
+ segment(section_kind kind): kind_(kind) {}
+
section_kind kind_;
- std::vector<parameter_list> mechanisms_;
+ std::vector<mechanism_spec> mechanisms_;
};
class placeholder_segment : public segment {
public:
- using base = segment;
- using base::kind_;
- using base::value_type;
- using base::size_type;
-
- placeholder_segment()
- {
- kind_ = section_kind::none;
- }
+ placeholder_segment(): segment(section_kind::none) {}
std::unique_ptr<segment> clone() const override {
// use default copy constructor
@@ -191,26 +149,10 @@ public:
class soma_segment : public segment {
public:
- using base = segment;
- using base::kind_;
- using base::value_type;
- using base::point_type;
- using base::size_type;
-
soma_segment() = delete;
- soma_segment(value_type r):
- radius_{r}
- {
- kind_ = section_kind::soma;
- mechanisms_.push_back(membrane_parameters());
- }
-
- soma_segment(value_type r, point_type const& c):
- soma_segment(r)
- {
- center_ = c;
- }
+ explicit soma_segment(value_type r, point_type c = point_type{}):
+ segment(section_kind::soma), radius_{r}, center_(c) {}
std::unique_ptr<segment> clone() const override {
// use default copy constructor
@@ -273,46 +215,24 @@ public:
cable_segment() = delete;
// constructors for a cable with no location information
- cable_segment(
- section_kind k,
- std::vector<value_type> r,
- std::vector<value_type> lens
- ) {
- kind_ = k;
- assert(k==section_kind::dendrite || k==section_kind::axon);
-
- radii_ = std::move(r);
- lengths_ = std::move(lens);
-
- // add default membrane parameters
- mechanisms_.push_back(membrane_parameters());
+ cable_segment(section_kind k, std::vector<value_type> r, std::vector<value_type> lens):
+ segment(k), radii_(std::move(r)), lengths_(std::move(lens))
+ {
+ assert(kind_==section_kind::dendrite || kind_==section_kind::axon);
}
- cable_segment(
- section_kind k,
- value_type r1,
- value_type r2,
- value_type len
- )
- : cable_segment{k, std::vector<value_type>{r1, r2}, std::vector<value_type>{len}}
- { }
+ cable_segment(section_kind k, value_type r1, value_type r2, value_type len):
+ //cable_segment{k, std::vector<value_type>{r1, r2}, std::vector<value_type>{len}}
+ cable_segment{k, {r1, r2}, decltype(lengths_){len}}
+ {}
// constructor that lets the user describe the cable as a
// seriew of radii and locations
- cable_segment(
- section_kind k,
- std::vector<value_type> r,
- std::vector<point_type> p
- ) {
- kind_ = k;
- assert(k==section_kind::dendrite || k==section_kind::axon);
-
- radii_ = std::move(r);
- locations_ = std::move(p);
+ cable_segment(section_kind k, std::vector<value_type> r, std::vector<point_type> p):
+ segment(k), radii_(std::move(r)), locations_(std::move(p))
+ {
+ assert(kind_==section_kind::dendrite || kind_==section_kind::axon);
update_lengths();
-
- // add default membrane parameters
- mechanisms_.push_back(membrane_parameters());
}
// helper that lets user specify with the radius and location
@@ -324,9 +244,9 @@ public:
value_type r2,
point_type const& p1,
point_type const& p2
- )
- : cable_segment(k, {r1, r2}, {p1, p2})
- { }
+ ):
+ cable_segment(k, {r1, r2}, {p1, p2})
+ {}
std::unique_ptr<segment> clone() const override {
// use default copy constructor
@@ -435,8 +355,7 @@ public:
}
private:
- void update_lengths()
- {
+ void update_lengths() {
if (locations_.size()) {
lengths_.resize(num_sub_segments());
for (size_type i=0; i<num_sub_segments(); ++i) {
@@ -446,8 +365,8 @@ private:
}
size_type num_compartments_ = 1;
- std::vector<value_type> lengths_;
std::vector<value_type> radii_;
+ std::vector<value_type> lengths_;
std::vector<point_type> locations_;
};
diff --git a/src/util/simple_table.hpp b/src/util/simple_table.hpp
new file mode 100644
index 00000000..32d6ee42
--- /dev/null
+++ b/src/util/simple_table.hpp
@@ -0,0 +1,44 @@
+#pragma once
+
+#include <iterator>
+#include <cstring>
+#include <utility>
+
+// Linear lookup of values in a table indexed by a key.
+//
+// Tables are any sequence of pairs or tuples, with the key as the first element.
+
+namespace arb {
+namespace util {
+
+namespace impl {
+ struct key_equal {
+ template <typename U, typename V>
+ bool operator()(U&& u, V&& v) const {
+ return std::forward<U>(u)==std::forward<V>(v);
+ }
+
+ // special case for C strings:
+ bool operator()(const char* u, const char* v) const {
+ return !std::strcmp(u, v);
+ }
+ };
+};
+
+// Return pointer to value (second element in entry) in table if key found,
+// otherwise nullptr.
+
+template <typename PairSeq, typename Key, typename Eq = impl::key_equal>
+auto table_lookup(PairSeq&& seq, const Key& key, Eq eq = Eq{})
+ -> decltype(&std::get<1>(*std::begin(seq)))
+{
+ for (auto&& entry: seq) {
+ if (eq(std::get<0>(entry), key)) {
+ return &std::get<1>(entry);
+ }
+ }
+ return nullptr;
+}
+
+} // namespace util
+} // namespace arb
diff --git a/tests/common_cells.hpp b/tests/common_cells.hpp
index dc6d0d19..088967db 100644
--- a/tests/common_cells.hpp
+++ b/tests/common_cells.hpp
@@ -4,7 +4,7 @@
#include <recipe.hpp>
#include <segment.hpp>
#include <math.hpp>
-#include <parameter_list.hpp>
+#include <mechinfo.hpp>
namespace arb {
@@ -14,7 +14,7 @@ namespace arb {
* Soma:
* diameter: 18.8 µm
* mechanisms: HH (default params)
- * bulk resistivitiy: 100 Ω·cm
+ * bulk resistivitiy: 100 Ω·cm [default]
* capacitance: 0.01 F/m² [default]
*
* Stimuli:
@@ -25,8 +25,7 @@ inline cell make_cell_soma_only(bool with_stim = true) {
cell c;
auto soma = c.add_soma(18.8/2.0);
- soma->mechanism("membrane").set("r_L", 100);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
if (with_stim) {
c.add_stimulus({0,0.5}, {10., 100., 0.1});
@@ -39,7 +38,7 @@ inline cell make_cell_soma_only(bool with_stim = true) {
* Create cell with a soma and unbranched dendrite:
*
* Common properties:
- * bulk resistivity: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm [default]
* capacitance: 0.01 F/m² [default]
*
* Soma:
@@ -50,7 +49,6 @@ inline cell make_cell_soma_only(bool with_stim = true) {
* mechanisms: passive (default params)
* diameter: 1 µm
* length: 200 µm
- * bulk resistivity: 100 Ω·cm
* compartments: 4
*
* Stimulus:
@@ -61,14 +59,13 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
cell c;
auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
c.add_cable(0, section_kind::dendrite, 1.0/2, 1.0/2, 200.0);
for (auto& seg: c.segments()) {
- seg->mechanism("membrane").set("r_L", 100);
if (seg->is_dendrite()) {
- seg->add_mechanism(pas_parameters());
+ seg->add_mechanism("pas");
seg->set_compartments(4);
}
}
@@ -83,7 +80,7 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
* Create cell with a soma and unbranched tapered dendrite:
*
* Common properties:
- * bulk resistivity: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm [default]
* capacitance: 0.01 F/m² [default]
*
* Soma:
@@ -106,14 +103,13 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
cell c;
auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
c.add_cable(0, section_kind::dendrite, 1.0/2, 0.4/2, 200.0);
for (auto& seg: c.segments()) {
- seg->mechanism("membrane").set("r_L", 100);
if (seg->is_dendrite()) {
- seg->add_mechanism(pas_parameters());
+ seg->add_mechanism("pas");
seg->set_compartments(4);
}
}
@@ -127,12 +123,16 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
/*
* Create cell with a soma and unbranched dendrite with varying diameter:
*
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm [default]
+ * capacitance: 0.01 F/m² [default]
+ *
* Soma:
* mechanisms: HH
* diameter: 12.6157 µm
*
* Dendrite:
- * mechanisms: none
+ * mechanisms: passive (default params)
* length: 100 µm
* membrane resistance: 100 Ω·cm
* compartments: 4
@@ -145,7 +145,7 @@ inline cell make_cell_ball_and_squiggle(bool with_stim = true) {
cell c;
auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
std::vector<cell::value_type> radii;
std::vector<cell::point_type> points;
@@ -166,9 +166,8 @@ inline cell make_cell_ball_and_squiggle(bool with_stim = true) {
c.add_cable(0, std::move(dendrite));
for (auto& seg: c.segments()) {
- seg->mechanism("membrane").set("r_L", 100);
if (seg->is_dendrite()) {
- seg->add_mechanism(pas_parameters());
+ seg->add_mechanism("pas");
seg->set_compartments(4);
}
}
@@ -185,7 +184,7 @@ inline cell make_cell_ball_and_squiggle(bool with_stim = true) {
* O----======
*
* Common properties:
- * bulk resistivity: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm [default]
* capacitance: 0.01 F/m² [default]
*
* Soma:
@@ -207,16 +206,15 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
cell c;
auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
for (auto& seg: c.segments()) {
- seg->mechanism("membrane").set("r_L", 100);
if (seg->is_dendrite()) {
- seg->add_mechanism(pas_parameters());
+ seg->add_mechanism("pas");
seg->set_compartments(4);
}
}
@@ -238,7 +236,7 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
* membrane reversal potential: -65 mV (default)
* diameter: 1 µm
* length: 1000 µm
- * bulk resistivity: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm [default]
* capacitance: 0.01 F/m² [default]
* compartments: 4
*
@@ -247,13 +245,6 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
*
* Note: zero-volume soma added with same mechanisms, as
* work-around for some existing fvm modelling issues.
- *
- * TODO: Set the correct values when parameters are generally
- * settable!
- *
- * We can't currently change leak parameters
- * from defaults, so we scale other electrical parameters
- * proportionally.
*/
inline cell make_cell_simple_cable(bool with_stim = true) {
@@ -267,22 +258,13 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
double gbar = 0.000025;
double I = 0.1;
- // fudge factor! can't change passive membrane
- // conductance from gbar0 = 0.001
-
- double gbar0 = 0.001;
- double f = gbar/gbar0;
-
- // scale everything else
- r_L *= f;
- c_m /= f;
- I /= f;
+ mechanism_spec pas("pas");
+ pas["g"] = gbar;
for (auto& seg: c.segments()) {
- seg->add_mechanism(pas_parameters());
- seg->mechanism("membrane").set("r_L", r_L);
- seg->mechanism("membrane").set("c_m", c_m);
- // seg->mechanism("pas").set("g", gbar);
+ seg->add_mechanism(pas);
+ seg->rL = r_L;
+ seg->cm = c_m;
if (seg->is_dendrite()) {
seg->set_compartments(4);
@@ -296,5 +278,4 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
}
return c;
}
-
} // namespace arb
diff --git a/tests/simple_recipes.hpp b/tests/simple_recipes.hpp
index 338b0a55..5417c0fb 100644
--- a/tests/simple_recipes.hpp
+++ b/tests/simple_recipes.hpp
@@ -34,8 +34,22 @@ public:
pvec_.push_back({probe_id, tag, std::move(address)});
}
+ void add_specialized_mechanism(std::string name, specialized_mechanism m) {
+ cell_gprop.special_mechs[name] = std::move(m);
+ }
+
+ util::any get_global_properties(cell_kind k) const override {
+ switch (k) {
+ case cable1d_neuron:
+ return cell_gprop;
+ default:
+ return util::any{};
+ }
+ }
+
protected:
std::unordered_map<cell_gid_type, std::vector<probe_info>> probes_;
+ cell_global_properties cell_gprop;
};
// Convenience derived recipe class for wrapping n copies of a single
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 11808a13..62eb4695 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -57,7 +57,7 @@ set(TEST_SOURCES
test_multi_event_stream.cpp
test_nop.cpp
test_optional.cpp
- test_parameters.cpp
+ test_mechinfo.cpp
test_partition.cpp
test_path.cpp
test_point.cpp
diff --git a/tests/unit/common.hpp b/tests/unit/common.hpp
index af9833f5..445f5992 100644
--- a/tests/unit/common.hpp
+++ b/tests/unit/common.hpp
@@ -5,6 +5,9 @@
* more than one unit test.
*/
+#include <utility>
+#include <cmath>
+
#include "../gtest.h"
namespace testing {
@@ -143,4 +146,13 @@ template <typename FPType, typename Seq1, typename Seq2>
return ::testing::AssertionSuccess();
}
+// Assert two floating point values are within a relative tolerance.
+inline ::testing::AssertionResult near_relative(double a, double b, double relerr) {
+ double tol = relerr*std::max(std::abs(a), std::abs(b));
+ if (std::abs(a-b)>tol) {
+ return ::testing::AssertionFailure() << "relative error between floating point numbers " << a << " and " << b << " exceeds tolerance " << relerr;
+ }
+ return ::testing::AssertionSuccess();
+}
+
} // namespace testing
diff --git a/tests/unit/test.cpp b/tests/unit/test.cpp
index d0a2d43a..d971b2f9 100644
--- a/tests/unit/test.cpp
+++ b/tests/unit/test.cpp
@@ -3,9 +3,13 @@
#include <numeric>
#include <vector>
+#include <communication/global_policy.hpp>
+
#include "../gtest.h"
int main(int argc, char **argv) {
+ arb::communication::global_policy_guard g(argc, argv);
+
::testing::InitGoogleTest(&argc, argv);
return RUN_ALL_TESTS();
}
diff --git a/tests/unit/test_cell.cpp b/tests/unit/test_cell.cpp
index 8fbcabd4..5c197809 100644
--- a/tests/unit/test_cell.cpp
+++ b/tests/unit/test_cell.cpp
@@ -178,8 +178,7 @@ TEST(cell_type, clone)
c.add_cable(1, section_kind::dendrite, 0.2, 0.15, 20);
c.segment(2)->set_compartments(5);
- parameter_list exp_default("expsyn");
- c.add_synapse({1, 0.3}, exp_default);
+ c.add_synapse({1, 0.3}, mechanism_spec("expsyn"));
c.add_detector({0, 0.5}, 10.0);
diff --git a/tests/unit/test_fvm_multi.cpp b/tests/unit/test_fvm_multi.cpp
index eb7cafce..a9c57ec0 100644
--- a/tests/unit/test_fvm_multi.cpp
+++ b/tests/unit/test_fvm_multi.cpp
@@ -1,5 +1,7 @@
#include <cstddef>
#include <fstream>
+#include <set>
+#include <vector>
#include "../gtest.h"
@@ -7,12 +9,17 @@
#include <cell.hpp>
#include <common_types.hpp>
#include <fvm_multicell.hpp>
+#include <load_balance.hpp>
+#include <model.hpp>
#include <recipe.hpp>
#include <sampler_map.hpp>
+#include <sampling.hpp>
+#include <schedule.hpp>
+#include <segment.hpp>
#include <util/meta.hpp>
#include <util/rangeutil.hpp>
-#include "../test_util.hpp"
+#include "common.hpp"
#include "../common_cells.hpp"
#include "../simple_recipes.hpp"
@@ -53,19 +60,13 @@ TEST(fvm_multi, init)
const auto m = cell.model();
EXPECT_EQ(m.tree.num_segments(), 2u);
- auto& soma_hh = cell.soma()->mechanism("hh");
+ auto& soma_hh = (cell.soma()->mechanism("hh")).get();
soma_hh.set("gnabar", 0.12);
soma_hh.set("gkbar", 0.036);
soma_hh.set("gl", 0.0003);
soma_hh.set("el", -54.3);
- // check that parameter values were set correctly
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gnabar").value, 0.12);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gkbar").value, 0.036);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
-
cell.segment(1)->set_compartments(10);
std::vector<fvm_cell::target_handle> targets;
@@ -124,10 +125,10 @@ TEST(fvm_multi, multi_init)
EXPECT_EQ(cells[1].segment(2)->num_compartments(), 4u);
EXPECT_EQ(cells[1].segment(3)->num_compartments(), 4u);
- cells[0].add_synapse({1, 0.4}, parameter_list("expsyn"));
- cells[0].add_synapse({1, 0.4}, parameter_list("expsyn"));
- cells[1].add_synapse({2, 0.4}, parameter_list("exp2syn"));
- cells[1].add_synapse({3, 0.4}, parameter_list("expsyn"));
+ cells[0].add_synapse({1, 0.4}, "expsyn");
+ cells[0].add_synapse({1, 0.4}, "expsyn");
+ cells[1].add_synapse({2, 0.4}, "exp2syn");
+ cells[1].add_synapse({3, 0.4}, "expsyn");
cells[1].add_detector({0, 0}, 3.3);
@@ -265,7 +266,7 @@ TEST(fvm_multi, mechanism_indexes)
cell c;
auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
// add dendrite of length 200 um and diameter 1 um with passive channel
c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
@@ -273,13 +274,13 @@ TEST(fvm_multi, mechanism_indexes)
c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
auto& segs = c.segments();
- segs[1]->add_mechanism(pas_parameters());
- segs[2]->add_mechanism(hh_parameters());
- segs[3]->add_mechanism(pas_parameters());
+ segs[1]->add_mechanism("pas");
+ segs[2]->add_mechanism("hh");
+ segs[3]->add_mechanism("pas");
for (auto& seg: segs) {
if (seg->is_dendrite()) {
- seg->mechanism("membrane").set("r_L", 100);
+ seg->rL = 100;
seg->set_compartments(4);
}
}
@@ -328,6 +329,338 @@ TEST(fvm_multi, mechanism_indexes)
}
}
+namespace {
+ double wm_impl(double wa, double xa) {
+ return wa? xa/wa: 0;
+ }
+
+ template <typename... R>
+ double wm_impl(double wa, double xa, double w, double x, R... rest) {
+ return wm_impl(wa+w, xa+w*x, rest...);
+ }
+
+ // Computed weighted mean (w*x + ...) / (w + ...).
+ template <typename... R>
+ double wmean(double w, double x, R... rest) {
+ return wm_impl(w, w*x, rest...);
+ }
+}
+
+// Test area-weighted linear combination of density mechanism parameters.
+
+TEST(fvm_multi, density_weights) {
+ using namespace arb;
+
+ // Create a cell with 4 segments:
+ // - Soma (segment 0) plus three dendrites (1, 2, 3) meeting at a branch point.
+ // - HH mechanism on all segments.
+ // - Dendritic segments are given 3 compartments each.
+ //
+ // The CV corresponding to the branch point should comprise the terminal
+ // 1/6 of segment 1 and the initial 1/6 of segments 2 and 3.
+ //
+ // The HH mechanism current density parameters ('gnabar', 'gkbar' and 'gl') are set
+ // differently for each segment:
+ //
+ // soma: all default values (gnabar = 0.12, gkbar = .036, gl = .0003)
+ // segment 1: gl = .0002
+ // segment 2: gkbar = .05
+ // segment 3: gkbar = .07, gl = .0004
+ //
+ // Geometry:
+ // segment 1: 100 µm long, 1 µm diameter cylinder.
+ // segment 2: 200 µm long, diameter linear taper from 1 µm to 0.2 µm.
+ // segment 3: 150 µm long, 0.8 µm diameter cylinder.
+ //
+ // Use divided compartment view on segments to compute area contributions.
+
+ cell c;
+ auto soma = c.add_soma(12.6157/2.0);
+
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.1, 200);
+ c.add_cable(1, section_kind::dendrite, 0.4, 0.4, 150);
+
+ auto& segs = c.segments();
+
+ double dflt_gkbar = .036;
+ double dflt_gnabar = 0.12;
+ double dflt_gl = 0.0003;
+
+ double seg1_gl = .0002;
+ double seg2_gkbar = .05;
+ double seg3_gkbar = .0004;
+ double seg3_gl = .0004;
+
+ for (int i = 0; i<4; ++i) {
+ segment& seg = *segs[i];
+ seg.set_compartments(3);
+
+ mechanism_spec hh("hh");
+ switch (i) {
+ case 1:
+ hh["gl"] = seg1_gl;
+ break;
+ case 2:
+ hh["gkbar"] = seg2_gkbar;
+ break;
+ case 3:
+ hh["gkbar"] = seg3_gkbar;
+ hh["gl"] = seg3_gl;
+ break;
+ default: ;
+ }
+ seg.add_mechanism(hh);
+ }
+
+ int ncv = 10;
+ std::vector<double> expected_gkbar(ncv, dflt_gkbar);
+ std::vector<double> expected_gnabar(ncv, dflt_gnabar);
+ std::vector<double> expected_gl(ncv, dflt_gl);
+
+ double soma_area = soma->area();
+ auto seg1_divs = div_compartments<div_compartment_by_ends>(segs[1]->as_cable());
+ auto seg2_divs = div_compartments<div_compartment_by_ends>(segs[2]->as_cable());
+ auto seg3_divs = div_compartments<div_compartment_by_ends>(segs[3]->as_cable());
+
+ // CV 0: mix of soma and left of segment 1
+ expected_gl[0] = wmean(soma_area, dflt_gl, seg1_divs(0).left.area, seg1_gl);
+
+ expected_gl[1] = seg1_gl;
+ expected_gl[2] = seg1_gl;
+
+ // CV 3: mix of right of segment 1 and left of segments 2 and 3.
+ expected_gkbar[3] = wmean(seg1_divs(2).right.area, dflt_gkbar, seg2_divs(0).left.area, seg2_gkbar, seg3_divs(0).left.area, seg3_gkbar);
+ expected_gl[3] = wmean(seg1_divs(2).right.area, seg1_gl, seg2_divs(0).left.area, dflt_gl, seg3_divs(0).left.area, seg3_gl);
+
+ // CV 4-6: just segment 2
+ expected_gkbar[4] = seg2_gkbar;
+ expected_gkbar[5] = seg2_gkbar;
+ expected_gkbar[6] = seg2_gkbar;
+
+ // CV 7-9: just segment 3
+ expected_gkbar[7] = seg3_gkbar;
+ expected_gkbar[8] = seg3_gkbar;
+ expected_gkbar[9] = seg3_gkbar;
+ expected_gl[7] = seg3_gl;
+ expected_gl[8] = seg3_gl;
+ expected_gl[9] = seg3_gl;
+
+ // Generate the lowered fvm cell.
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
+
+ fvm_cell fvcell;
+ fvcell.initialize({0}, cable1d_recipe(c), targets, probe_map);
+
+ // Check CV area assumptions.
+ // Note: area integrator used here and in `fvm_multicell` may differ, and so areas computed may
+ // differ some due to rounding area, even given that we're dealing with simple truncated cones
+ // for segments. Check relative error within a tolerance of (say) 10 epsilon.
+ auto cv_areas = fvcell.cv_areas();
+ double area_relerr = 10*std::numeric_limits<double>::epsilon();
+ EXPECT_TRUE(testing::near_relative(cv_areas[0],
+ soma_area+seg1_divs(0).left.area, area_relerr));
+ EXPECT_TRUE(testing::near_relative(cv_areas[1],
+ seg1_divs(0).right.area+seg1_divs(1).left.area, area_relerr));
+ EXPECT_TRUE(testing::near_relative(cv_areas[3],
+ seg1_divs(2).right.area+seg2_divs(0).left.area+seg3_divs(0).left.area, area_relerr));
+ EXPECT_TRUE(testing::near_relative(cv_areas[6],
+ seg2_divs(2).right.area, area_relerr));
+
+ // Grab the HH parameters from the mechanism.
+ EXPECT_EQ(1u, fvcell.mechanisms().size());
+ auto& hh_mech = *fvcell.mechanisms().front();
+
+ auto gnabar_field = hh_mech.field_view_ptr("gnabar");
+ auto gkbar_field = hh_mech.field_view_ptr("gkbar");
+ auto gl_field = hh_mech.field_view_ptr("gl");
+
+ EXPECT_TRUE(testing::seq_almost_eq<double>(expected_gnabar, hh_mech.*gnabar_field));
+ EXPECT_TRUE(testing::seq_almost_eq<double>(expected_gkbar, hh_mech.*gkbar_field));
+ EXPECT_TRUE(testing::seq_almost_eq<double>(expected_gl, hh_mech.*gl_field));
+}
+
+// Test specialized mechanism behaviour.
+
+TEST(fvm_multi, specialized_mechs) {
+ using namespace arb;
+
+ // Create ball and stick cells with the 'test_kin1' mechanism, which produces
+ // a voltage-independent current density of the form a + exp(-t/tau) as a function
+ // of time t.
+ //
+ // 1. Default 'test_kin1': tau = 10 [ms].
+ //
+ // 2. Specialized version 'custom_kin1' with tau = 20 [ms].
+ //
+ // 3. Cell with both test_kin1 and custom_kin1.
+
+ specialized_mechanism custom_kin1 = {"test_kin1", {{"tau", 20.0}}};
+
+ cell cells[3];
+
+ for (int i = 0; i<3; ++i) {
+ cell& c = cells[i];
+ c.add_soma(6.0);
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
+
+ c.segment(1)->set_compartments(4);
+ for (auto& seg: c.segments()) {
+ if (!seg->is_soma()) {
+ seg->as_cable()->set_compartments(4);
+ }
+
+ switch (i) {
+ case 0:
+ seg->add_mechanism("test_kin1");
+ break;
+ case 1:
+ seg->add_mechanism("custom_kin1");
+ break;
+ case 2:
+ seg->add_mechanism("test_kin1");
+ seg->add_mechanism("custom_kin1");
+ break;
+ }
+ }
+ }
+
+ cable1d_recipe rec(cells);
+ rec.add_specialized_mechanism("custom_kin1", custom_kin1);
+
+ cell_probe_address where{{1, 0.3}, cell_probe_address::membrane_current};
+ rec.add_probe(0, 0, where);
+ rec.add_probe(1, 0, where);
+ rec.add_probe(2, 0, where);
+
+ {
+ // Test initialization and global parameter values.
+
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
+
+ fvm_cell fvcell;
+ fvcell.initialize({0, 1, 2}, rec, targets, probe_map);
+
+ std::map<std::string, fvm_cell::mechanism*> mechmap;
+ for (auto& m: fvcell.mechanisms()) {
+ // (names of mechanisms should _all_ be 'test_kin1', but aliases will differ)
+ EXPECT_EQ("test_kin1", m->name());
+ mechmap[m->alias()] = m.get();
+ }
+
+ ASSERT_EQ(2u, mechmap.size());
+ EXPECT_NE(0u, mechmap.count("test_kin1"));
+ EXPECT_NE(0u, mechmap.count("custom_kin1"));
+
+ // Both mechanisms are of the same type, so we can use the
+ // same member pointer.
+ auto fptr = mechmap.begin()->second->field_value_ptr("tau");
+
+ ASSERT_NE(nullptr, fptr);
+
+ EXPECT_EQ(10.0, mechmap["test_kin1"]->*fptr);
+ EXPECT_EQ(20.0, mechmap["custom_kin1"]->*fptr);
+ }
+
+ {
+ // Test dynamics:
+ // 1. Current at same point on cell 0 at time 10 ms should equal that
+ // on cell 1 at time 20 ms.
+ // 2. Current for cell 2 should be sum of currents for cells 0 and 1 at any given time.
+
+ std::vector<double> samples[3];
+
+ sampler_function sampler = [&](cell_member_type pid, probe_tag, std::size_t n, const sample_record* records) {
+ for (std::size_t i = 0; i<n; ++i) {
+ double v = *util::any_cast<const double*>(records[i].data);
+ samples[pid.gid].push_back(v);
+ }
+ };
+
+ float times[] = {10.f, 20.f};
+
+ auto decomp = partition_load_balance(rec, hw::node_info{1u, 0u});
+ model sim(rec, decomp);
+ sim.add_sampler(all_probes, explicit_schedule(times), sampler);
+ sim.run(30.0, 1.f/1024);
+
+ ASSERT_EQ(2u, samples[0].size());
+ ASSERT_EQ(2u, samples[1].size());
+ ASSERT_EQ(2u, samples[2].size());
+
+ // Integration isn't exact: let's aim for one part in 10'000.
+ double relerr = 0.0001;
+ EXPECT_TRUE(testing::near_relative(samples[0][0], samples[1][1], relerr));
+ EXPECT_TRUE(testing::near_relative(samples[0][0]+samples[1][0], samples[2][0], relerr));
+ EXPECT_TRUE(testing::near_relative(samples[0][1]+samples[1][1], samples[2][1], relerr));
+ }
+}
+
+// Test synapses with differing parameter settings.
+
+TEST(fvm_multi, synapse_parameters) {
+ using namespace arb;
+
+ cell c;
+ c.add_soma(6.0);
+ c.add_cable(0, section_kind::dendrite, 0.4, 0.4, 100.0);
+ c.segment(1)->set_compartments(4);
+
+ // Add synapses out-of-order, with a parameter value a function of position
+ // on the segment, so we can test that parameters are properly associated
+ // after re-ordering.
+
+ struct pset {
+ double x; // segment position
+ double tau1;
+ double tau2;
+ };
+
+ pset settings[] = {
+ {0.8, 1.5, 2.5},
+ {0.1, 1.6, 3.7},
+ {0.5, 1.7, 3.6},
+ {0.6, 0.8, 2.5},
+ {0.4, 0.9, 3.4},
+ {0.9, 1.1, 2.3}
+ };
+
+ for (auto s: settings) {
+ mechanism_spec m("exp2syn");
+ c.add_synapse({1, s.x}, m.set("tau1", s.tau1).set("tau2", s.tau2));
+ }
+
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
+
+ fvm_cell fvcell;
+ fvcell.initialize({0}, cable1d_recipe(c), targets, probe_map);
+
+ EXPECT_EQ(1u, fvcell.mechanisms().size());
+ auto& exp2syn_mech = *fvcell.mechanisms().front();
+
+ auto tau1_ptr = exp2syn_mech.field_view_ptr("tau1");
+ auto tau2_ptr = exp2syn_mech.field_view_ptr("tau2");
+
+ // Compare tau1, tau2 values from settings and from mechanism, ignoring order.
+ std::set<std::pair<double, double>> expected;
+ for (auto s: settings) {
+ expected.insert({s.tau1, s.tau2});
+ }
+
+ unsigned n = exp2syn_mech.size();
+ ASSERT_EQ(util::size(settings), n);
+
+ std::set<std::pair<double, double>> values;
+ for (unsigned i = 0; i<n; ++i) {
+ values.insert({(exp2syn_mech.*tau1_ptr)[i], (exp2syn_mech.*tau2_ptr)[i]});
+ }
+
+ EXPECT_EQ(expected, values);
+}
+
struct handle_info {
unsigned cell;
std::string mech;
@@ -378,7 +711,7 @@ void run_target_handle_test(std::vector<handle_info> all_handles) {
x.cv += 5; // offset for cell 1
}
- cells[x.cell].add_synapse({seg_id, pos}, parameter_list(x.mech));
+ cells[x.cell].add_synapse({seg_id, pos}, x.mech);
handles[x.cell].push_back(x);
}
diff --git a/tests/unit/test_matrix.cpp b/tests/unit/test_matrix.cpp
index 22f66f02..3bbab835 100644
--- a/tests/unit/test_matrix.cpp
+++ b/tests/unit/test_matrix.cpp
@@ -87,6 +87,8 @@ TEST(matrix, zero_diagonal)
// elements should be ignored).
// These submatrices should leave the rhs as-is when solved.
+ using memory::make_const_view;
+
// Three matrices, sizes 3, 3 and 2, with no branching.
std::vector<size_type> p = {0, 0, 1, 3, 3, 5, 5};
std::vector<size_type> c = {0, 3, 5, 7};
@@ -96,9 +98,9 @@ TEST(matrix, zero_diagonal)
EXPECT_EQ(3u, m.num_cells());
auto& A = m.state_;
- A.d = vvec({2, 3, 2, 0, 0, 4, 5});
- A.u = vvec({0, -1, -1, 0, -1, 0, -2});
- A.rhs = vvec({3, 5, 7, 7, 8, 16, 32});
+ A.d = make_const_view(vvec({2, 3, 2, 0, 0, 4, 5}));
+ A.u = make_const_view(vvec({0, -1, -1, 0, -1, 0, -2}));
+ A.rhs = make_const_view(vvec({3, 5, 7, 7, 8, 16, 32}));
// Expected solution:
std::vector<value_type> expected = {4, 5, 6, 7, 8, 9, 10};
diff --git a/tests/unit/test_matrix.cu b/tests/unit/test_matrix.cu
index 195f82ea..c3156b62 100644
--- a/tests/unit/test_matrix.cu
+++ b/tests/unit/test_matrix.cu
@@ -291,7 +291,7 @@ TEST(matrix, assemble)
std::generate(dt.begin(), dt.end(), [&](){return dt_dist(gen);});
// Voltage and current values
- m_mc.assemble(host_array(dt), host_array(group_size, -64), host_array(group_size, 10));
+ m_mc.assemble(on_host(dt), host_array(group_size, -64), host_array(group_size, 10));
m_mc.solve();
m_gpu.assemble(on_gpu(dt), gpu_array(group_size, -64), gpu_array(group_size, 10));
m_gpu.solve();
diff --git a/tests/unit/test_mechinfo.cpp b/tests/unit/test_mechinfo.cpp
new file mode 100644
index 00000000..1e92ba08
--- /dev/null
+++ b/tests/unit/test_mechinfo.cpp
@@ -0,0 +1,45 @@
+#include <map>
+#include <string>
+#include <vector>
+
+#include "mechinfo.hpp"
+
+#include "../gtest.h"
+#include "../test_util.hpp"
+
+// TODO: expand tests when we have exported mechanism schemata
+// from modcc.
+
+using namespace arb;
+
+TEST(mechanism_spec, setting) {
+ mechanism_spec m("foo");
+
+ m.set("a", 3.2);
+ m.set("b", 4.3);
+
+ auto dflt = m["c"];
+ EXPECT_EQ(0., dflt); // note: 0 default is artefact of dummy schema
+
+ EXPECT_EQ(3.2, m["a"]);
+ EXPECT_EQ(4.3, m["b"]);
+
+ m["b"] = 5.4;
+ m["d"] = 6.5;
+
+ EXPECT_EQ(3.2, m["a"]);
+ EXPECT_EQ(5.4, m["b"]);
+ EXPECT_EQ(6.5, m["d"]);
+
+ // Check values() method is faithful:
+ auto p = m.values();
+ EXPECT_EQ(3u, p.size());
+
+ EXPECT_TRUE(p.count("a"));
+ EXPECT_TRUE(p.count("b"));
+ EXPECT_TRUE(p.count("d"));
+
+ EXPECT_EQ(p["a"], m["a"]);
+ EXPECT_EQ(p["b"], m["b"]);
+ EXPECT_EQ(p["d"], m["d"]);
+}
diff --git a/tests/unit/test_parameters.cpp b/tests/unit/test_parameters.cpp
deleted file mode 100644
index f8e29ca3..00000000
--- a/tests/unit/test_parameters.cpp
+++ /dev/null
@@ -1,41 +0,0 @@
-#include <fstream>
-
-#include "../gtest.h"
-#include "../test_util.hpp"
-
-#include "../src/parameter_list.hpp"
-
-// test out the parameter infrastructure
-TEST(parameters, setting)
-{
- arb::parameter_list list("test");
- EXPECT_EQ(list.name(), "test");
- EXPECT_EQ(list.num_parameters(), 0);
-
- arb::parameter p("a", 0.12, {0, 10});
-
- // add_parameter() returns a bool that indicates whether
- // it was able to successfull add the parameter
- EXPECT_TRUE(list.add_parameter(std::move(p)));
- EXPECT_EQ(list.num_parameters(), 1);
-
- // test in place construction of a parameter
- EXPECT_TRUE(list.add_parameter({"b", -3.0}));
- EXPECT_EQ(list.num_parameters(), 2);
-
- // check that adding a parameter that already exists returns false
- // and does not increase the number of parameters
- EXPECT_FALSE(list.add_parameter({"b", -3.0}));
- EXPECT_EQ(list.num_parameters(), 2);
-
- auto& parms = list.parameters();
- EXPECT_EQ(parms[0].name, "a");
- EXPECT_EQ(parms[0].value, 0.12);
- EXPECT_EQ(parms[0].range.min, 0);
- EXPECT_EQ(parms[0].range.max, 10);
-
- EXPECT_EQ(parms[1].name, "b");
- EXPECT_EQ(parms[1].value, -3);
- EXPECT_FALSE(parms[1].range.has_lower_bound());
- EXPECT_FALSE(parms[1].range.has_upper_bound());
-}
diff --git a/tests/unit/test_synapses.cpp b/tests/unit/test_synapses.cpp
index 804218f3..7951b4fb 100644
--- a/tests/unit/test_synapses.cpp
+++ b/tests/unit/test_synapses.cpp
@@ -16,14 +16,11 @@ TEST(synapses, add_to_cell)
// Soma with diameter 12.6157 um and HH channel
auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
+ soma->add_mechanism("hh");
- parameter_list exp_default("expsyn");
- parameter_list exp2_default("exp2syn");
-
- cell.add_synapse({0, 0.1}, exp_default);
- cell.add_synapse({1, 0.2}, exp2_default);
- cell.add_synapse({0, 0.3}, exp_default);
+ cell.add_synapse({0, 0.1}, "expsyn");
+ cell.add_synapse({1, 0.2}, "exp2syn");
+ cell.add_synapse({0, 0.3}, "expsyn");
EXPECT_EQ(3u, cell.synapses().size());
const auto& syns = cell.synapses();
@@ -44,6 +41,7 @@ TEST(synapses, add_to_cell)
TEST(synapses, expsyn_basic_state)
{
using namespace arb;
+ using memory::make_const_view;
using size_type = multicore::backend::size_type;
using value_type = multicore::backend::value_type;
@@ -62,7 +60,7 @@ TEST(synapses, expsyn_basic_state)
synapse_type::array voltage(num_comp, -65.0);
synapse_type::array current(num_comp, 1.0);
- auto mech = make_mechanism<synapse_type>(0, cell_index, time, time_to, dt, voltage, current, weights, node_index);
+ auto mech = make_mechanism<synapse_type>(0, cell_index, time, time_to, dt, voltage, current, make_const_view(weights), make_const_view(node_index));
auto ptr = dynamic_cast<synapse_type*>(mech.get());
auto n = ptr->size();
@@ -105,6 +103,7 @@ TEST(synapses, expsyn_basic_state)
TEST(synapses, exp2syn_basic_state)
{
using namespace arb;
+ using memory::make_const_view;
using size_type = multicore::backend::size_type;
using value_type = multicore::backend::value_type;
@@ -123,7 +122,7 @@ TEST(synapses, exp2syn_basic_state)
synapse_type::array voltage(num_comp, -65.0);
synapse_type::array current(num_comp, 1.0);
- auto mech = make_mechanism<synapse_type>(0, cell_index, time, time_to, dt, voltage, current, weights, node_index);
+ auto mech = make_mechanism<synapse_type>(0, cell_index, time, time_to, dt, voltage, current, make_const_view(weights), make_const_view(node_index));
auto ptr = dynamic_cast<synapse_type*>(mech.get());
auto n = ptr->size();
diff --git a/tests/validation/validate_kinetic.cpp b/tests/validation/validate_kinetic.cpp
index 986d8cbc..b0077df1 100644
--- a/tests/validation/validate_kinetic.cpp
+++ b/tests/validation/validate_kinetic.cpp
@@ -23,6 +23,7 @@
void run_kinetic_dt(
arb::backend_kind backend,
arb::cell& c,
+ arb::cell_probe_address probe,
float t_end,
nlohmann::json meta,
const std::string& ref_file)
@@ -32,7 +33,7 @@ void run_kinetic_dt(
float sample_dt = g_trace_io.sample_dt();
cable1d_recipe rec{c};
- rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
+ rec.add_probe(0, 0, probe);
probe_label plabels[1] = {"soma.mid", {0u, 0u}};
meta["sim"] = "arbor";
@@ -71,7 +72,8 @@ void validate_kinetic_kin1(arb::backend_kind backend) {
// 20 µm diameter soma with single mechanism, current probe
cell c;
auto soma = c.add_soma(10);
- soma->add_mechanism(std::string("test_kin1"));
+ soma->add_mechanism("test_kin1");
+ cell_probe_address probe{{0, 0.5}, cell_probe_address::membrane_current};
nlohmann::json meta = {
{"model", "test_kin1"},
@@ -79,7 +81,7 @@ void validate_kinetic_kin1(arb::backend_kind backend) {
{"units", "nA"}
};
- run_kinetic_dt(backend, c, 100.f, meta, "numeric_kin1.json");
+ run_kinetic_dt(backend, c, probe, 100.f, meta, "numeric_kin1.json");
}
void validate_kinetic_kinlva(arb::backend_kind backend) {
@@ -89,7 +91,8 @@ void validate_kinetic_kinlva(arb::backend_kind backend) {
cell c;
auto soma = c.add_soma(10);
c.add_stimulus({0,0.5}, {20., 130., -0.025});
- soma->add_mechanism(std::string("test_kinlva"));
+ soma->add_mechanism("test_kinlva");
+ cell_probe_address probe{{0, 0.5}, cell_probe_address::membrane_voltage};
nlohmann::json meta = {
{"model", "test_kinlva"},
@@ -97,7 +100,7 @@ void validate_kinetic_kinlva(arb::backend_kind backend) {
{"units", "mV"}
};
- run_kinetic_dt(backend, c, 300.f, meta, "numeric_kinlva.json");
+ run_kinetic_dt(backend, c, probe, 300.f, meta, "numeric_kinlva.json");
}
diff --git a/tests/validation/validate_synapses.cpp b/tests/validation/validate_synapses.cpp
index 03c2ba70..55fc0e62 100644
--- a/tests/validation/validate_synapses.cpp
+++ b/tests/validation/validate_synapses.cpp
@@ -38,7 +38,7 @@ void run_synapse_test(
};
cell c = make_cell_ball_and_stick(false); // no stimuli
- parameter_list syn_default(syn_type);
+ mechanism_spec syn_default(syn_type);
c.add_synapse({1, 0.5}, syn_default);
// injected spike events
--
GitLab