From b63cca1c50a2a4d0d7d24f1b8de1fdd5581806cf Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Mon, 21 Nov 2022 11:15:24 +0100
Subject: [PATCH] =?UTF-8?q?=E2=9A=A1=20Voltage=20Processes=20(#2033)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Add the VOLTAGE_PROCESS mechanism kind to modcc, allowing for direct writing to the membrane voltage
Obviously these are extremely sharp tools and can break the cable model, so we add some constraints to their
use (see docs). Still, user discretion is required.
As a sneak peek for what this might be useful for
* implementing voltage clamps
* limiting membrane voltage (irritatingly also a kind of clamping)
* faking ABSTRACT_CELL like behaviour (although that might require a bit more work)
Closes #1343
---
.gitignore | 10 +-
arbor/arbexcept.cpp | 5 +
arbor/backends/gpu/shared_state.cpp | 41 ++-
arbor/backends/gpu/shared_state.hpp | 8 +-
arbor/backends/multicore/shared_state.cpp | 51 ++--
arbor/backends/multicore/shared_state.hpp | 8 +-
arbor/cable_cell.cpp | 7 +
arbor/fvm_layout.cpp | 226 +++++++++++----
arbor/fvm_lowered_cell_impl.hpp | 63 +++--
arbor/include/arbor/arbexcept.hpp | 6 +
arbor/include/arbor/cable_cell.hpp | 4 +-
arbor/include/arbor/cable_cell_param.hpp | 11 +
arbor/include/arbor/lif_cell.hpp | 2 +-
arbor/include/arbor/mechanism_abi.h | 5 +-
arbor/include/arbor/simd/simd.hpp | 6 +
arborio/cableio.cpp | 4 +
doc/fileformat/nmodl.rst | 41 ++-
example/CMakeLists.txt | 1 +
example/probe-demo/probe-demo.cpp | 16 +-
example/v_clamp/CMakeLists.txt | 3 +
example/v_clamp/README.md | 38 +++
example/v_clamp/v-clamp.cpp | 167 +++++++++++
mechanisms/CMakeLists.txt | 2 +-
mechanisms/default/v_clamp.mod | 12 +
mechanisms/default/v_limit.mod | 13 +
modcc/blocks.cpp | 10 +-
modcc/identifier.hpp | 2 +-
modcc/mechanism.hpp | 43 ---
modcc/module.cpp | 27 +-
modcc/module.hpp | 3 +
modcc/parser.cpp | 6 +-
modcc/printer/cprinter.cpp | 64 ++++-
modcc/printer/cprinter.hpp | 6 +-
modcc/printer/gpuprinter.cpp | 34 ++-
modcc/printer/printerutil.hpp | 1 +
modcc/token.cpp | 2 +
modcc/token.hpp | 2 +-
python/cells.cpp | 18 +-
python/example/v-clamp.py | 55 ++++
scripts/run_cpp_examples.sh | 2 +-
scripts/run_python_examples.sh | 1 +
test/unit-modcc/mod_files/test_v_process.mod | 13 +
test/unit-modcc/test_module.cpp | 9 +
test/unit-modcc/test_parser.cpp | 10 +
test/unit/CMakeLists.txt | 1 +
test/unit/test_abi.cpp | 8 +-
test/unit/test_fvm_lowered.cpp | 4 +-
test/unit/test_kinetic_linear.cpp | 2 +-
test/unit/test_mech_temp_diam.cpp | 4 +-
test/unit/test_s_expr.cpp | 4 +
test/unit/test_segment_tree.cpp | 2 +-
test/unit/test_synapses.cpp | 4 +-
test/unit/test_v_clamp.cpp | 276 +++++++++++++++++++
53 files changed, 1136 insertions(+), 227 deletions(-)
create mode 100644 example/v_clamp/CMakeLists.txt
create mode 100644 example/v_clamp/README.md
create mode 100644 example/v_clamp/v-clamp.cpp
create mode 100644 mechanisms/default/v_clamp.mod
create mode 100644 mechanisms/default/v_limit.mod
delete mode 100644 modcc/mechanism.hpp
create mode 100755 python/example/v-clamp.py
create mode 100644 test/unit-modcc/mod_files/test_v_process.mod
create mode 100644 test/unit/test_v_clamp.cpp
diff --git a/.gitignore b/.gitignore
index 036ba491..4b1dccf5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,7 +66,7 @@ mechanisms/gpu/*.hpp
mechanisms/gpu/*.cu
# build path
-build*
+build/*
commit.msg
@@ -89,3 +89,11 @@ dist
# generated image files by Python examples
python/example/*.svg
python/example/*.csv
+
+# made by direnv
+.direnv
+.envrc
+
+# build/compile caches
+_skbuild
+.ccls-cache
diff --git a/arbor/arbexcept.cpp b/arbor/arbexcept.cpp
index 74dca037..21ccedc1 100644
--- a/arbor/arbexcept.cpp
+++ b/arbor/arbexcept.cpp
@@ -37,6 +37,11 @@ bad_cell_description::bad_cell_description(cell_kind kind, cell_gid_type gid):
kind(kind)
{}
+invalid_mechanism_kind::invalid_mechanism_kind(arb_mechanism_kind kind):
+ arbor_exception(pprintf("Invalid mechanism kind: {})", kind)),
+ kind(kind)
+{}
+
bad_connection_source_gid::bad_connection_source_gid(cell_gid_type gid, cell_gid_type src_gid, cell_size_type num_cells):
arbor_exception(pprintf("Model building error on cell {}: connection source gid {} is out of range: there are only {} cells in the model, in the range [{}:{}].", gid, src_gid, num_cells, 0, num_cells-1)),
gid(gid), src_gid(src_gid), num_cells(num_cells)
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index d6c66065..ab6ae516 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -216,24 +216,6 @@ shared_state::shared_state(
add_scalar(temperature_degC.size(), temperature_degC.data(), -273.15);
}
-void shared_state::set_parameter(mechanism& m, const std::string& key, const std::vector<arb_value_type>& values) {
- if (values.size()!=m.ppack_.width) throw arbor_internal_error("mechanism parameter size mismatch");
- const auto& store = storage.at(m.mechanism_id());
-
- arb_value_type* data = nullptr;
- for (arb_size_type i = 0; i<m.mech_.n_parameters; ++i) {
- if (key==m.mech_.parameters[i].name) {
- data = store.parameters_[i];
- break;
- }
- }
-
- if (!data) throw arbor_internal_error(util::pprintf("no such mechanism parameter '{}'", key));
-
- if (!m.ppack_.width) return;
- memory::copy(memory::make_const_view(values), memory::device_view<arb_value_type>(data, m.ppack_.width));
-}
-
void shared_state::update_prng_state(mechanism& m) {
if (!m.mech_.n_random_variables) return;
auto const mech_id = m.mechanism_id();
@@ -252,7 +234,11 @@ const arb_value_type* shared_state::mechanism_state_data(const mechanism& m, con
return nullptr;
}
-void shared_state::instantiate(mechanism& m, unsigned id, const mechanism_overrides& overrides, const mechanism_layout& pos_data) {
+void shared_state::instantiate(mechanism& m,
+ unsigned id,
+ const mechanism_overrides& overrides,
+ const mechanism_layout& pos_data,
+ const std::vector<std::pair<std::string, std::vector<arb_value_type>>>& params) {
assert(m.iface_.backend == arb_backend_kind_gpu);
using util::make_range;
using util::make_span;
@@ -313,19 +299,28 @@ void shared_state::instantiate(mechanism& m, unsigned id, const mechanism_overri
{
// Allocate bulk storage
std::size_t count = (m.mech_.n_state_vars + m.mech_.n_parameters + 1)*width_padded + m.mech_.n_globals;
- store.data_ = array(count, NAN);
+ store.data_ = array(count, NAN);
chunk_writer writer(store.data_.data(), width_padded);
// First sub-array of data_ is used for weight_
m.ppack_.weight = writer.append_with_padding(pos_data.weight, 0);
- // Set fields
+ // Set parameters to either default or explicit setting
for (auto idx: make_span(m.mech_.n_parameters)) {
- store.parameters_[idx] = writer.fill(m.mech_.parameters[idx].default_value);
+ const auto& param = m.mech_.parameters[idx];
+ const auto& it = std::find_if(params.begin(), params.end(),
+ [&](const auto& k) { return k.first == param.name; });
+ if (it != params.end()) {
+ if (it->second.size() != width) throw arbor_internal_error("mechanism field size mismatch");
+ store.parameters_[idx] = writer.append_with_padding(it->second, param.default_value);
+ }
+ else {
+ store.parameters_[idx] = writer.fill(param.default_value);
+ }
}
+ // Make STATE var the default
for (auto idx: make_span(m.mech_.n_state_vars)) {
store.state_vars_[idx] = writer.fill(m.mech_.state_vars[idx].default_value);
}
-
// Assign global scalar parameters. NB: Last chunk, since it breaks the width striding.
for (auto idx: make_span(m.mech_.n_globals)) store.globals_[idx] = m.mech_.globals[idx].default_value;
for (auto& [k, v]: overrides.globals) {
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index c09cdd12..e08de3ae 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -177,9 +177,11 @@ struct ARB_ARBOR_API shared_state {
);
// Setup a mechanism and tie its backing store to this object
- void instantiate(arb::mechanism&, unsigned, const mechanism_overrides&, const mechanism_layout&);
-
- void set_parameter(mechanism&, const std::string&, const std::vector<arb_value_type>&);
+ void instantiate(mechanism&,
+ unsigned,
+ const mechanism_overrides&,
+ const mechanism_layout&,
+ const std::vector<std::pair<std::string, std::vector<arb_value_type>>>&);
void update_prng_state(mechanism&);
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index e7c2e8f8..1ae38936 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -437,27 +437,6 @@ std::size_t extend_width(const arb::mechanism& mech, std::size_t width) {
}
} // anonymous namespace
-void shared_state::set_parameter(mechanism& m, const std::string& key, const std::vector<arb_value_type>& values) {
- if (values.size()!=m.ppack_.width) throw arbor_internal_error("mechanism field size mismatch");
-
- bool found = false;
- arb_value_type* data = nullptr;
- for (arb_size_type i = 0; i<m.mech_.n_parameters; ++i) {
- if (key==m.mech_.parameters[i].name) {
- data = m.ppack_.parameters[i];
- found = true;
- break;
- }
- }
-
- if (!found) throw arbor_internal_error(util::pprintf("no such parameter '{}'", key));
-
- if (!m.ppack_.width) return;
-
- auto width_padded = extend_width<arb_value_type>(m, m.ppack_.width);
- copy_extend(values, util::range_n(data, width_padded), values.back());
-}
-
const arb_value_type* shared_state::mechanism_state_data(const mechanism& m, const std::string& key) {
for (arb_size_type i = 0; i<m.mech_.n_state_vars; ++i) {
if (key==m.mech_.state_vars[i].name) {
@@ -508,7 +487,11 @@ void shared_state::update_prng_state(mechanism& m) {
// * For indices in the padded tail of node_index_, set index to last valid CV index.
// * For indices in the padded tail of ion index maps, set index to last valid ion index.
-void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_overrides& overrides, const mechanism_layout& pos_data) {
+void shared_state::instantiate(arb::mechanism& m,
+ unsigned id,
+ const mechanism_overrides& overrides,
+ const mechanism_layout& pos_data,
+ const std::vector<std::pair<std::string, std::vector<arb_value_type>>>& params) {
// Mechanism indices and data require:
// * an alignment that is a multiple of the mechansim data_alignment();
// * a size which is a multiple of partition_width() for SIMD access.
@@ -518,9 +501,14 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
util::padded_allocator<> pad(m.data_alignment());
+ if (storage.find(id) != storage.end()) {
+ throw arbor_internal_error("Duplicate mechanism id in MC shared state.");
+ }
+ auto& store = storage[id];
+ auto width = pos_data.cv.size();
// Assign non-owning views onto shared state:
m.ppack_ = {0};
- m.ppack_.width = pos_data.cv.size();
+ m.ppack_.width = width;
m.ppack_.mechanism_id = id;
m.ppack_.vec_ci = cv_to_cell.data();
m.ppack_.vec_di = cv_to_intdom.data();
@@ -537,9 +525,6 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
bool mult_in_place = !pos_data.multiplicity.empty();
bool peer_indices = !pos_data.peer_cv.empty();
- if (storage.find(id) != storage.end()) throw arb::arbor_internal_error("Duplicate mech id in shared state");
- auto& store = storage[id];
-
// store indices for random number generation
store.gid_ = pos_data.gid;
store.idx_ = pos_data.idx;
@@ -585,10 +570,20 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
// First sub-array of data_ is used for weight_
m.ppack_.weight = writer.append(pos_data.weight, 0);
- // Set fields
+ // Set parameters: either default, or explicit override
for (auto idx: make_span(m.mech_.n_parameters)) {
- m.ppack_.parameters[idx] = writer.fill(m.mech_.parameters[idx].default_value);
+ const auto& param = m.mech_.parameters[idx];
+ const auto& it = std::find_if(params.begin(), params.end(),
+ [&](const auto& k) { return k.first == param.name; });
+ if (it != params.end()) {
+ if (it->second.size() != width) throw arbor_internal_error("mechanism field size mismatch");
+ m.ppack_.parameters[idx] = writer.append(it->second, param.default_value);
+ }
+ else {
+ m.ppack_.parameters[idx] = writer.fill(param.default_value);
+ }
}
+ // Set initial state values
for (auto idx: make_span(m.mech_.n_state_vars)) {
m.ppack_.state_vars[idx] = writer.fill(m.mech_.state_vars[idx].default_value);
}
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index 7d306ce7..dda56f51 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -185,9 +185,11 @@ struct ARB_ARBOR_API shared_state {
arb_seed_type cbprng_seed_ = 0u
);
- void instantiate(mechanism&, unsigned, const mechanism_overrides&, const mechanism_layout&);
-
- void set_parameter(mechanism&, const std::string&, const std::vector<arb_value_type>&);
+ void instantiate(mechanism&,
+ unsigned,
+ const mechanism_overrides&,
+ const mechanism_layout&,
+ const std::vector<std::pair<std::string, std::vector<arb_value_type>>>&);
void update_prng_state(mechanism&);
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 8b7d0382..34c4f39b 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -55,6 +55,9 @@ std::string show(const paintable& item) {
else if constexpr (std::is_same_v<density, T>) {
os << "density:" << p.mech.name();
}
+ else if constexpr (std::is_same_v<voltage_process, T>) {
+ os << "voltage-process:" << p.mech.name();
+ }
},
item);
return os.str();
@@ -136,6 +139,10 @@ struct cable_cell_impl {
return region_map.get<T>();
}
+ mcable_map<voltage_process>& get_region_map(const voltage_process& v) {
+ return region_map.get<voltage_process>()[v.mech.name()];
+ }
+
mcable_map<std::pair<density, iexpr_map>> &
get_region_map(const density &desc) {
return region_map.get<density>()[desc.mech.name()];
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 94a8e107..7862419d 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -81,12 +81,13 @@ std::vector<V> unique_union(const std::vector<V>& a, const std::vector<V>& b) {
}
return u;
}
+
} // anonymous namespace
// Building CV geometry
// --------------------
-
+//
// CV geometry
cv_geometry::cv_geometry(const cable_cell& cell, const locset& ls):
@@ -794,6 +795,16 @@ ARB_ARBOR_API std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> f
return gj_conns;
}
+std::unordered_map<std::string, fvm_mechanism_config>
+make_voltage_mechanism_config(const cable_cell_global_properties& gprop,
+ const region_assignment<voltage_process> assignments,
+ const mechanism_catalogue& catalogue,
+ iexpr_ptr unit_scale,
+ unsigned cell_idx,
+ const fvm_cv_discretization& D,
+ const concrete_embedding& embedding,
+ const mprovider& provider);
+
fvm_mechanism_data fvm_build_mechanism_data(
const cable_cell_global_properties& gprop,
const cable_cell& cell,
@@ -821,6 +832,47 @@ ARB_ARBOR_API fvm_mechanism_data fvm_build_mechanism_data(
return combined;
}
+
+// Verify mechanism ion usage, parameter values.
+void verify_mechanism(const cable_cell_global_properties& gprop,
+ const fvm_cv_discretization& D,
+ const mechanism_info& info,
+ const mechanism_desc& desc) {
+ const auto& global_ions = gprop.ion_species;
+
+ for (const auto& pv: desc.values()) {
+ if (!info.parameters.count(pv.first)) {
+ throw no_such_parameter(desc.name(), pv.first);
+ }
+ if (!info.parameters.at(pv.first).valid(pv.second)) {
+ throw invalid_parameter_value(desc.name(), pv.first, pv.second);
+ }
+ }
+
+ for (const auto& [ion, dep]: info.ions) {
+ if (!global_ions.count(ion)) {
+ throw cable_cell_error(
+ "mechanism "+desc.name()+" uses ion "+ion+ " which is missing in global properties");
+ }
+
+ if (dep.verify_ion_charge) {
+ if (dep.expected_ion_charge!=global_ions.at(ion)) {
+ throw cable_cell_error(
+ "mechanism "+desc.name()+" uses ion "+ion+ " expecting a different valence");
+ }
+ }
+
+ if (dep.write_reversal_potential && (dep.write_concentration_int || dep.write_concentration_ext)) {
+ throw cable_cell_error("mechanism "+desc.name()+" writes both reversal potential and concentration");
+ }
+
+ auto is_diffusive = D.diffusive_ions.count(ion);
+ if (dep.access_concentration_diff && !is_diffusive) {
+ throw illegal_diffusive_mechanism(desc.name(), ion);
+ }
+ }
+}
+
// Construct FVM mechanism data for a single cell.
fvm_mechanism_data fvm_build_mechanism_data(
@@ -834,8 +886,9 @@ fvm_mechanism_data fvm_build_mechanism_data(
using index_type = arb_index_type;
using value_type = arb_value_type;
- const mechanism_catalogue& catalogue = gprop.catalogue;
+ const auto& catalogue = gprop.catalogue;
const auto& embedding = cell.embedding();
+ const auto& provider = cell.provider();
const auto& global_dflt = gprop.default_parameters;
const auto& dflt = cell.default_parameters();
@@ -847,43 +900,6 @@ fvm_mechanism_data fvm_build_mechanism_data(
const auto& assignments = cell.region_assignments();
- // Verify mechanism ion usage, parameter values.
- auto verify_mechanism = [&gprop, &D](const mechanism_info& info, const mechanism_desc& desc) {
- const auto& global_ions = gprop.ion_species;
-
- for (const auto& pv: desc.values()) {
- if (!info.parameters.count(pv.first)) {
- throw no_such_parameter(desc.name(), pv.first);
- }
- if (!info.parameters.at(pv.first).valid(pv.second)) {
- throw invalid_parameter_value(desc.name(), pv.first, pv.second);
- }
- }
-
- for (const auto& [ion, dep]: info.ions) {
- if (!global_ions.count(ion)) {
- throw cable_cell_error(
- "mechanism "+desc.name()+" uses ion "+ion+ " which is missing in global properties");
- }
-
- if (dep.verify_ion_charge) {
- if (dep.expected_ion_charge!=global_ions.at(ion)) {
- throw cable_cell_error(
- "mechanism "+desc.name()+" uses ion "+ion+ " expecting a different valence");
- }
- }
-
- if (dep.write_reversal_potential && (dep.write_concentration_int || dep.write_concentration_ext)) {
- throw cable_cell_error("mechanism "+desc.name()+" writes both reversal potential and concentration");
- }
-
- auto is_diffusive = D.diffusive_ions.count(ion);
- if (dep.access_concentration_diff && !is_diffusive) {
- throw illegal_diffusive_mechanism(desc.name(), ion);
- }
- }
- };
-
// Track ion usage of mechanisms so that ions are only instantiated where required.
std::unordered_map<std::string, std::vector<index_type>> ion_support;
@@ -905,8 +921,22 @@ fvm_mechanism_data fvm_build_mechanism_data(
std::unordered_map<std::string, mcable_map<double>> init_iconc_mask;
std::unordered_map<std::string, mcable_map<double>> init_econc_mask;
- iexpr_ptr unit_scale = thingify(iexpr::scalar(1.0), cell.provider());
-
+ iexpr_ptr unit_scale = thingify(iexpr::scalar(1.0), provider);
+
+ // Voltage mechanisms
+
+ {
+ auto assigns = assignments.get<voltage_process>();
+ auto configs = make_voltage_mechanism_config(gprop,
+ assigns,
+ catalogue,
+ unit_scale,
+ cell_idx,
+ D,
+ embedding,
+ provider);
+ M.mechanisms.insert(configs.begin(), configs.end());
+ }
// Density mechanisms:
for (const auto& [name, cables]: assignments.get<density>()) {
mechanism_info info = catalogue[name];
@@ -941,7 +971,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
const auto& mech = density_iexpr.first.mech;
const auto& scale_expr = density_iexpr.second;
- verify_mechanism(info, mech);
+ verify_mechanism(gprop, D, info, mech);
const auto& set_params = mech.values();
support.insert(cable, 1.);
@@ -969,7 +999,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
c.branch,
pw_over_cable(
param_maps[i], c, 0., [&](const mcable &c, const auto& x) {
- return x.first * x.second->eval(cell.provider(), c);
+ return x.first * x.second->eval(provider, c);
}));
}
}
@@ -1061,7 +1091,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
std::size_t offset = 0;
for (const placed<synapse>& pm: entry.second) {
const auto& mech = pm.item.mech;
- verify_mechanism(info, mech);
+ verify_mechanism(gprop, D, info, mech);
synapse_instance in{};
@@ -1210,7 +1240,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
for (const placed<junction>& pm: placements) {
const auto& mech = pm.item.mech;
- verify_mechanism(info, mech);
+ verify_mechanism(gprop, D, info, mech);
const auto& set_params = mech.values();
junction_desc per_lid;
@@ -1369,12 +1399,12 @@ fvm_mechanism_data fvm_build_mechanism_data(
}
if (auto di = D.diffusive_ions.find(ion); di != D.diffusive_ions.end()) {
- config.is_diffusive = true;
- config.face_diffusivity = di->second.face_diffusivity;
+ config.is_diffusive = true;
+ config.face_diffusivity = di->second.face_diffusivity;
}
- config.econc_written = write_xo.count(ion);
- config.iconc_written = write_xi.count(ion);
+ config.econc_written = write_xo.count(ion);
+ config.iconc_written = write_xi.count(ion);
if (!config.cv.empty()) M.ions[ion] = std::move(config);
}
@@ -1391,7 +1421,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
throw cable_cell_error("expected reversal potential mechanism for ion " +ion +", got "+ revpot.name() +" which has " +arb_mechanism_kind_str(info.kind));
}
- verify_mechanism(info, revpot);
+ verify_mechanism(gprop, D, info, revpot);
revpot_specified.insert(ion);
bool writes_this_revpot = false;
@@ -1464,4 +1494,100 @@ fvm_mechanism_data fvm_build_mechanism_data(
return M;
}
+std::unordered_map<std::string, fvm_mechanism_config>
+make_voltage_mechanism_config(const cable_cell_global_properties& gprop,
+ const region_assignment<voltage_process> assignments,
+ const mechanism_catalogue& catalogue,
+ iexpr_ptr unit_scale,
+ unsigned cell_idx,
+ const fvm_cv_discretization& D,
+ const concrete_embedding& embedding,
+ const mprovider& provider) {
+ std::unordered_map<std::string, fvm_mechanism_config> result;
+ std::unordered_set<mcable> voltage_support;
+ for (const auto& [name, cables]: assignments) {
+ auto info = catalogue[name];
+ if (info.kind != arb_mechanism_kind_voltage) {
+ throw cable_cell_error("expected voltage mechanism, got " +name +" which has " +arb_mechanism_kind_str(info.kind));
+ }
+
+ fvm_mechanism_config config;
+ config.kind = arb_mechanism_kind_voltage;
+
+ std::vector<std::string> param_names;
+ assign(param_names, util::keys(info.parameters));
+ sort(param_names);
+
+ std::vector<double> param_dflt;
+ std::size_t n_param = param_names.size();
+ param_dflt.reserve(n_param);
+ config.param_values.reserve(n_param);
+
+ for (std::size_t i = 0; i<n_param; ++i) {
+ const auto& p = param_names[i];
+ config.param_values.emplace_back(p, std::vector<arb_value_type>{});
+ param_dflt.push_back(info.parameters.at(p).default_value);
+ }
+
+ mcable_map<double> support;
+ std::vector<mcable_map<std::pair<double, iexpr_ptr>>> param_maps;
+
+ param_maps.resize(n_param);
+
+ for (const auto& [cable, density_iexpr]: cables) {
+ const auto& mech = density_iexpr.mech;
+
+ verify_mechanism(gprop, D, info, mech);
+ const auto& set_params = mech.values();
+
+ support.insert(cable, 1.);
+ for (std::size_t i = 0; i<n_param; ++i) {
+ double value = value_by_key(set_params, param_names[i]).value_or(param_dflt[i]);
+ param_maps[i].insert(cable, {value, unit_scale});
+ }
+ }
+
+ std::vector<double> param_on_cv(n_param);
+
+ for (auto cv: D.geometry.cell_cvs(cell_idx)) {
+ double area = 0;
+ util::fill(param_on_cv, 0.);
+
+ for (mcable c: D.geometry.cables(cv)) {
+ double area_on_cable = embedding.integrate_area(c.branch, pw_over_cable(support, c, 0.));
+ if (!area_on_cable) continue;
+ area += area_on_cable;
+ for (std::size_t i = 0; i<n_param; ++i) {
+ auto map = param_maps[i];
+ auto fun = [&](const auto& c, const auto& x) {
+ return x.first * x.second->eval(provider, c);
+ };
+ auto pw = pw_over_cable(map, c, 0., fun);
+ param_on_cv[i] += embedding.integrate_area(c.branch, pw);
+ }
+ }
+
+ if (area>0) {
+ config.cv.push_back(cv);
+ config.norm_area.push_back(area/D.cv_area[cv]);
+
+ double oo_area = 1./area;
+ for (auto i: count_along(param_on_cv)) {
+ config.param_values[i].second.push_back(param_on_cv[i]*oo_area);
+ }
+ }
+ }
+
+ for (const auto& [cable, _]: support) {
+ if (voltage_support.count(cable)) {
+ throw cable_cell_error("Multiple voltage processes on a single cable");
+ }
+ voltage_support.insert(cable);
+ }
+ if (!config.cv.empty()) result.emplace(name, std::move(config));
+ }
+ return result;
+}
+
+
} // namespace arb
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 09902f5c..b3d09765 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -98,6 +98,7 @@ private:
value_type tmin_ = 0;
std::vector<mechanism_ptr> mechanisms_; // excludes reversal potential calculators.
std::vector<mechanism_ptr> revpot_mechanisms_;
+ std::vector<mechanism_ptr> voltage_mechanisms_;
// Optional non-physical voltage check threshold
std::optional<double> check_voltage_mV_;
@@ -166,6 +167,10 @@ void fvm_lowered_cell_impl<Backend>::reset() {
state_->reset();
set_tmin(0);
+ for (auto& m: voltage_mechanisms_) {
+ m->initialize();
+ }
+
for (auto& m: revpot_mechanisms_) {
m->initialize();
}
@@ -188,6 +193,10 @@ void fvm_lowered_cell_impl<Backend>::reset() {
m->initialize();
}
+ for (auto& m: voltage_mechanisms_) {
+ m->initialize();
+ }
+
// NOTE: Threshold watcher reset must come after the voltage values are set,
// as voltage is implicitly read by watcher to set initial state.
threshold_watcher_.reset(state_->voltage);
@@ -286,21 +295,28 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
state_->integrate_diffusion();
PL();
- // Integrate mechanism state.
-
+ // Integrate mechanism state for density
for (auto& m: mechanisms_) {
state_->update_prng_state(*m);
m->update_state();
}
// Update ion concentrations.
-
PE(advance:integrate:ionupdate);
update_ion_state();
PL();
- // Update time and test for spike threshold crossings.
+ // voltage mechs run now; after the cable_solver, but before the
+ // threshold test
+ for (auto& m: voltage_mechanisms_) {
+ m->update_current();
+ }
+ for (auto& m: voltage_mechanisms_) {
+ state_->update_prng_state(*m);
+ m->update_state();
+ }
+ // Update time and test for spike threshold crossings.
PE(advance:integrate:threshold);
threshold_watcher_.test(&state_->time_since_spike);
PL();
@@ -545,10 +561,7 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
std::unordered_map<std::string, mechanism*> mechptr_by_name;
unsigned mech_id = 0;
- for (auto& m: mech_data.mechanisms) {
- auto& name = m.first;
- auto& config = m.second;
-
+ for (const auto& [name, config]: mech_data.mechanisms) {
mechanism_layout layout;
layout.cv = config.cv;
layout.multiplicity = config.multiplicity;
@@ -599,6 +612,7 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
layout.weight[i] = config.local_weight[i] * 1000/D.cv_area[cv];
}
break;
+ case arb_mechanism_kind_voltage:
case arb_mechanism_kind_density:
// Current density contributions from mechanism are already in [A/m²].
@@ -617,19 +631,28 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
break;
}
- auto minst = mech_instance(name);
- state_->instantiate(*minst.mech, mech_id++, minst.overrides, layout);
- mechptr_by_name[name] = minst.mech.get();
-
- for (auto& pv: config.param_values) {
- state_->set_parameter(*minst.mech, pv.first, pv.second);
- }
+ auto [mech, over] = mech_instance(name);
+ state_->instantiate(*mech, mech_id, over, layout, config.param_values);
+ mechptr_by_name[name] = mech.get();
+ ++mech_id;
- if (config.kind==arb_mechanism_kind_reversal_potential) {
- revpot_mechanisms_.emplace_back(minst.mech.release());
- }
- else {
- mechanisms_.emplace_back(minst.mech.release());
+ switch (config.kind) {
+ case arb_mechanism_kind_gap_junction:
+ case arb_mechanism_kind_point:
+ case arb_mechanism_kind_density: {
+ mechanisms_.emplace_back(mech.release());
+ break;
+ }
+ case arb_mechanism_kind_reversal_potential: {
+ revpot_mechanisms_.emplace_back(mech.release());
+ break;
+ }
+ case arb_mechanism_kind_voltage: {
+ voltage_mechanisms_.emplace_back(mech.release());
+ break;
+ }
+ default:;
+ throw invalid_mechanism_kind(config.kind);
}
}
diff --git a/arbor/include/arbor/arbexcept.hpp b/arbor/include/arbor/arbexcept.hpp
index deb578e4..efa1255d 100644
--- a/arbor/include/arbor/arbexcept.hpp
+++ b/arbor/include/arbor/arbexcept.hpp
@@ -6,6 +6,7 @@
#include <arbor/common_types.hpp>
#include <arbor/export.hpp>
+#include <arbor/mechanism_abi.h>
// Arbor-specific exception hierarchy.
@@ -41,6 +42,11 @@ struct ARB_SYMBOL_VISIBLE bad_cell_probe: arbor_exception {
cell_kind kind;
};
+struct ARB_SYMBOL_VISIBLE invalid_mechanism_kind: arbor_exception {
+ invalid_mechanism_kind(arb_mechanism_kind);
+ arb_mechanism_kind kind;
+};
+
struct ARB_SYMBOL_VISIBLE bad_cell_description: arbor_exception {
bad_cell_description(cell_kind kind, cell_gid_type gid);
cell_gid_type gid;
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index cb0570bc..221c967d 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -213,7 +213,7 @@ using region_assignment = std::conditional_t<
std::conditional_t<
std::is_same<T, init_int_concentration>::value ||
std::is_same<T, init_ext_concentration>::value || std::is_same<T, init_reversal_potential>::value ||
- std::is_same<T, ion_diffusivity>::value,
+ std::is_same<T, ion_diffusivity>::value || std::is_same<T, voltage_process>::value,
std::unordered_map<std::string, mcable_map<T>>,
std::conditional_t<std::is_same_v<T, density>,
std::unordered_map<std::string, mcable_map<std::pair<T, iexpr_map>>>,
@@ -238,7 +238,7 @@ using location_assignment =
mlocation_map<T>>;
using cable_cell_region_map = static_typed_map<region_assignment,
- density, init_membrane_potential, axial_resistivity,
+ density, voltage_process, init_membrane_potential, axial_resistivity,
temperature_K, membrane_capacitance, init_int_concentration,
ion_diffusivity, init_ext_concentration, init_reversal_potential>;
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index b4307d0a..42748a57 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -227,6 +227,16 @@ struct ARB_SYMBOL_VISIBLE density {
}
};
+struct ARB_SYMBOL_VISIBLE voltage_process {
+ mechanism_desc mech;
+ explicit voltage_process(mechanism_desc m): mech(std::move(m)) {}
+ voltage_process(mechanism_desc m, const std::unordered_map<std::string, double>& params): mech(std::move(m)) {
+ for (const auto& [param, value]: params) {
+ mech.set(param, value);
+ }
+ }
+};
+
struct ARB_SYMBOL_VISIBLE ion_reversal_potential_method {
std::string ion;
mechanism_desc method;
@@ -255,6 +265,7 @@ using paintable =
init_ext_concentration,
init_reversal_potential,
density,
+ voltage_process,
scaled_mechanism<density>>;
using placeable =
diff --git a/arbor/include/arbor/lif_cell.hpp b/arbor/include/arbor/lif_cell.hpp
index 17ab8d8e..52fb6eff 100644
--- a/arbor/include/arbor/lif_cell.hpp
+++ b/arbor/include/arbor/lif_cell.hpp
@@ -20,7 +20,7 @@ struct ARB_SYMBOL_VISIBLE lif_cell {
double t_ref = 2; // Refractory period [ms].
lif_cell() = delete;
- lif_cell(cell_tag_type source, cell_tag_type target): source(std::move(source)), target(std::move(target)) {}
+ lif_cell(cell_tag_type source, cell_tag_type target): source(std::move(source)), target(std::move(target)) {}
};
// LIF probe metadata, to be passed to sampler callbacks. Intentionally left blank.
diff --git a/arbor/include/arbor/mechanism_abi.h b/arbor/include/arbor/mechanism_abi.h
index a5050ed8..bebe1aaa 100644
--- a/arbor/include/arbor/mechanism_abi.h
+++ b/arbor/include/arbor/mechanism_abi.h
@@ -18,8 +18,6 @@ extern "C" {
#define ARB_NO_ALIAS restrict
#endif
-
-
// Version
#define ARB_MECH_ABI_VERSION_MAJOR 0
#define ARB_MECH_ABI_VERSION_MINOR 3
@@ -35,6 +33,7 @@ typedef uint32_t arb_mechanism_kind;
#define arb_mechanism_kind_density 2
#define arb_mechanism_kind_reversal_potential 3
#define arb_mechanism_kind_gap_junction 4
+#define arb_mechanism_kind_voltage 5
typedef uint32_t arb_backend_kind;
#define arb_backend_kind_nil 0
@@ -46,6 +45,8 @@ inline const char* arb_mechanism_kind_str(const arb_mechanism_kind& mech) {
case arb_mechanism_kind_density: return "density mechanism kind";
case arb_mechanism_kind_point: return "point mechanism kind";
case arb_mechanism_kind_reversal_potential: return "reversal potential mechanism kind";
+ case arb_mechanism_kind_gap_junction: return "gap junction mechanism kind";
+ case arb_mechanism_kind_voltage: return "voltage mechanism kind";
default: return "unknown mechanism kind";
}
}
diff --git a/arbor/include/arbor/simd/simd.hpp b/arbor/include/arbor/simd/simd.hpp
index 966c472f..3c3ba9e7 100644
--- a/arbor/include/arbor/simd/simd.hpp
+++ b/arbor/include/arbor/simd/simd.hpp
@@ -206,6 +206,12 @@ namespace detail {
return *this;
}
+ template <typename Other>
+ indirect_indexed_expression& operator*=(const Other& s) {
+ compound_indexed_mul(s, p, index, width, constraint);
+ return *this;
+ }
+
template <typename Other>
indirect_indexed_expression& operator-=(const Other& s) {
compound_indexed_add(neg(s), p, index, width, constraint);
diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp
index 8be039c6..0afbad23 100644
--- a/arborio/cableio.cpp
+++ b/arborio/cableio.cpp
@@ -90,6 +90,9 @@ s_expr mksexp(const synapse& j) {
s_expr mksexp(const density& j) {
return slist("density"_symbol, mksexp(j.mech));
}
+s_expr mksexp(const voltage_process& j) {
+ return slist("voltage-process"_symbol, mksexp(j.mech));
+}
s_expr mksexp(const iexpr& j) {
std::stringstream s;
s << j;
@@ -674,6 +677,7 @@ eval_map named_evals{
{"junction", make_call<arb::mechanism_desc>(make_wrapped_mechanism<junction>, "'junction' with 1 argumnet (m: mechanism)")},
{"synapse", make_call<arb::mechanism_desc>(make_wrapped_mechanism<synapse>, "'synapse' with 1 argumnet (m: mechanism)")},
{"density", make_call<arb::mechanism_desc>(make_wrapped_mechanism<density>, "'density' with 1 argumnet (m: mechanism)")},
+ {"voltage-process", make_call<arb::mechanism_desc>(make_wrapped_mechanism<voltage_process>, "'voltage-process' with 1 argumnet (m: mechanism)")},
{"scaled-mechanism", make_scaled_mechanism_call<arb::density>("'scaled_mechanism' with a density argument, and 0 or more parameter scaling expressions"
"(d:density (param:string val:iexpr))")},
{"place", make_call<locset, i_clamp, std::string>(make_place, "'place' with 3 arguments (ls:locset c:current-clamp name:string)")},
diff --git a/doc/fileformat/nmodl.rst b/doc/fileformat/nmodl.rst
index a125cc82..f8015e40 100644
--- a/doc/fileformat/nmodl.rst
+++ b/doc/fileformat/nmodl.rst
@@ -159,7 +159,6 @@ Arbor-specific features
of a gap-junction connection as well as the local site. The peer membrane potential is
made available through the ``v_peer`` variable while the local membrane potential
is available through ``v``, as usual.
-
* Arbor offers a number of additional unary math functions which may offer improved performance
compared to hand-rolled solutions (especially with the vectorized and GPU backends).
All of the following functions take a single argument `x` and return a
@@ -175,6 +174,46 @@ Arbor-specific features
signum(x) sign of argument :math:`\begin{align*} +1 & ~~ \text{if} ~x \gt 0, \\ -1 & ~~ \text{if} ~x \lt 0, \\ 0 & ~~ \text{otherwise}. \end{align*}`
exprelr(x) guarded exponential :math:`x e^{1-x}`
================== ===================================== =========
+
+Voltage Processes
+-----------------
+
+Some cases require direct manipulation of the membrane voltage ``v``; which is
+normally prohibited and for good reason so. For these limited application,
+however, we offer mechanisms that are similar to ``density`` mechanism, but are
+tagged with ``VOLTAGE_PROCESS`` where normally ``SUFFIX`` would be used.
+
+This is both a very sharp tool and a somewhat experimental feature. Depending on
+our experience, it might be changed or removed. Using a ``VOLTAGE_PROCESS``,
+voltage clamping and limiting can be implemented, c.f. relevant examples in the
+``default`` catalogue. Example: limiting membrane voltage from above and below
+
+.. code:: none
+
+ NEURON {
+ VOLTAGE_PROCESS v_limit
+ GLOBAL v_low, v_high
+ }
+
+ PARAMETER {
+ v_high = 20 (mV)
+ v_low = -70 (mV)
+ }
+
+ BREAKPOINT {
+ v = max(min(v, v_high), v_low)
+ }
+
+As of the current implementation, we note the following details and constraints
+
+* only the ``INITIAL`` and ``BREAKPOINT`` procedures are called.
+* no ``WRITE`` access to ionic quantities is allowed.
+* only one ``VOLTAGE_PROCESS`` maybe present on a single location, adding more
+ results in an exception.
+* the ``BREAKPOINT`` callback will execute _after_ the cable solver. A
+ consequence of this is that if the initial membrane potential :math:`V_0` is
+ unequal to that of a potentially applied voltage clamp :math:`V_c`, the first
+ timestep will observe :math:`V_0`.
.. _format-sde:
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index b341b257..e97d41ac 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -14,4 +14,5 @@ add_subdirectory(probe-demo)
add_subdirectory(plasticity)
add_subdirectory(lfp)
add_subdirectory(diffusion)
+add_subdirectory(v_clamp)
add_subdirectory(ornstein_uhlenbeck)
diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp
index d7f21181..388832d7 100644
--- a/example/probe-demo/probe-demo.cpp
+++ b/example/probe-demo/probe-demo.cpp
@@ -199,14 +199,14 @@ void vector_sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_rec
for (unsigned j = 0; j<n_entities; ++j) {
std::cout << samples[i].time << ", ";
- if (cables_ptr) {
- arb::mcable where = (*cables_ptr)[j];
- std::cout << where.prox_pos << ", " << where.dist_pos << ", ";
- } else {
- arb::mlocation loc = (*point_info_ptr)[j].loc;
- std::cout << loc.pos << ", ";
- }
- std::cout << lo[j] << '\n';
+ if (cables_ptr) {
+ arb::mcable where = (*cables_ptr)[j];
+ std::cout << where.prox_pos << ", " << where.dist_pos << ", ";
+ } else {
+ arb::mlocation loc = (*point_info_ptr)[j].loc;
+ std::cout << loc.pos << ", ";
+ }
+ std::cout << lo[j] << '\n';
}
}
}
diff --git a/example/v_clamp/CMakeLists.txt b/example/v_clamp/CMakeLists.txt
new file mode 100644
index 00000000..43de3438
--- /dev/null
+++ b/example/v_clamp/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(voltage-clamp EXCLUDE_FROM_ALL v-clamp.cpp)
+add_dependencies(examples voltage-clamp)
+target_link_libraries(voltage-clamp PRIVATE arbor arborio arborenv arbor-sup ext-tinyopt)
diff --git a/example/v_clamp/README.md b/example/v_clamp/README.md
new file mode 100644
index 00000000..d1a69fb9
--- /dev/null
+++ b/example/v_clamp/README.md
@@ -0,0 +1,38 @@
+# Voltage Clamp Example.
+
+Modified example of single neuron with morphology described by an SWC file,
+a voltage clamp is added to certain regions.
+
+A cell is constructed from a supplied morphology with H–H channels
+on the soma and passive channels on the dendrites. A simple exponential
+synapse is added at the end of the last dendrite in the morphology,
+and is triggered at time t = 1 ms.
+
+The simulation outputs a trace of the soma membrane voltage in a simple CSV
+format.
+
+## Features
+
+The example demonstrates the use of:
+
+* Generating a morphology from an SWC file.
+* Using a morphology to construct a cable cell.
+* Injecting an artificial spike event into the simulation.
+* Adding a voltage probe to a cell and running a sampler on the simulation.
+
+## Running the example
+
+By default, `single-cell` will simulate a 'ball-and-stick' neuron for 20 ms,
+with a maxium dt of 0.025 ms and samples taken every 0.1 ms. The default
+synaptic weight is 0.01 µS.
+
+### Command line options
+
+| Option | Effect |
+|-----------------------|--------|
+| -m, --morphology FILE | Load the morphology from FILE in SWC format |
+| -d, --dt TIME | Set the maximum integration time step [ms] |
+| -t, --t-end TIME | Set the simulation duration [ms] |
+| -w, --weight WEIGHT | Set the synaptic weight [µS] |
+| -u, --voltage VALUE | Clamp soma to voltage [mV] |
+
diff --git a/example/v_clamp/v-clamp.cpp b/example/v_clamp/v-clamp.cpp
new file mode 100644
index 00000000..2d00a481
--- /dev/null
+++ b/example/v_clamp/v-clamp.cpp
@@ -0,0 +1,167 @@
+#include <any>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <arborio/label_parse.hpp>
+
+#include <arbor/load_balance.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/segment_tree.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/simple_sampler.hpp>
+
+#include <arborio/swcio.hpp>
+
+#include <tinyopt/tinyopt.h>
+
+using namespace arborio::literals;
+
+struct options {
+ std::string swc_file;
+ double t_end = 20;
+ double dt = 0.025;
+ float syn_weight = 0.01;
+ std::optional<double> voltage;
+ arb::cv_policy policy = arb::default_cv_policy();
+};
+
+options parse_options(int argc, char** argv);
+arb::morphology default_morphology();
+arb::morphology read_swc(const std::string& path);
+
+struct single_recipe: public arb::recipe {
+ explicit single_recipe(arb::morphology m, arb::cv_policy pol): morpho(std::move(m)) {
+ gprop.default_parameters = arb::neuron_parameter_defaults;
+ gprop.default_parameters.discretization = pol;
+ }
+
+ arb::cell_size_type num_cells() const override { return 1; }
+
+ std::vector<arb::probe_info> get_probes(arb::cell_gid_type) const override {
+ arb::mlocation mid_soma = {0, 0.5};
+ arb::cable_probe_membrane_voltage probe = {mid_soma};
+
+ // Probe info consists of a probe address and an optional tag, for use
+ // by the sampler callback.
+
+ return {arb::probe_info{probe, 0}};
+ }
+
+ arb::cell_kind get_cell_kind(arb::cell_gid_type) const override {
+ return arb::cell_kind::cable;
+ }
+
+ std::any get_global_properties(arb::cell_kind) const override {
+ return gprop;
+ }
+
+ arb::util::unique_any get_cell_description(arb::cell_gid_type) const override {
+ arb::label_dict dict;
+ using arb::reg::tagged;
+ dict.set("soma", tagged(1));
+ dict.set("dend", join(tagged(3), tagged(4), tagged(42)));
+
+ auto decor = arb::decor{}
+ // Add HH mechanism to soma, passive channels to dendrites.
+ .paint("soma"_lab, arb::density("hh"))
+ .paint("dend"_lab, arb::density("pas"))
+ // Add synapse to last branch.
+ .place(arb::mlocation{ morpho.num_branches()-1, 1. }, arb::synapse("exp2syn"), "synapse");
+ if (voltage_clamp) decor.paint("soma"_lab, arb::voltage_process("v_clamp"));
+
+ return arb::cable_cell(morpho, decor, dict);
+ }
+
+ std::optional<double> voltage_clamp;
+ arb::morphology morpho;
+ arb::cable_cell_global_properties gprop;
+};
+
+int main(int argc, char** argv) {
+ try {
+ options opt = parse_options(argc, argv);
+ single_recipe R(opt.swc_file.empty()? default_morphology(): read_swc(opt.swc_file), opt.policy);
+ R.voltage_clamp = opt.voltage;
+ arb::simulation sim(R);
+
+ // Attach a sampler to the probe described in the recipe, sampling every 0.1 ms.
+
+ arb::trace_vector<double> traces;
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1), arb::make_simple_sampler(traces));
+
+ // Trigger the single synapse (target is gid 0, index 0) at t = 1 ms with
+ // the given weight.
+
+ arb::spike_event spike = {0, 1., opt.syn_weight};
+ arb::cell_spike_events cell_spikes = {0, {spike}};
+ sim.inject_events({cell_spikes});
+
+ sim.run(opt.t_end, opt.dt);
+
+ for (auto entry: traces.at(0)) {
+ std::cout << entry.t << ", " << entry.v << "\n";
+ }
+ }
+ catch (std::exception& e) {
+ std::cerr << "caught exception: " << e.what() << "\n";
+ return 2;
+ }
+}
+
+options parse_options(int argc, char** argv) {
+ using namespace to;
+ options opt;
+
+ char** arg = argv+1;
+ while (*arg) {
+ if (auto dt = parse<double>(arg, "-d", "--dt")) {
+ opt.dt = dt.value();
+ }
+ else if (auto t_end = parse<double>(arg, "-t", "--t-end")) {
+ opt.t_end = t_end.value();
+ }
+ else if (auto weight = parse<float>(arg, "-w", "--weight")) {
+ opt.syn_weight = weight.value();
+ }
+ else if (auto voltage = parse<double>(arg, "-u", "--voltage")) {
+ opt.voltage = voltage.value();
+ }
+ else if (auto swc = parse<std::string>(arg, "-m", "--morphology")) {
+ opt.swc_file = swc.value();
+ }
+ else if (auto nseg = parse<unsigned>(arg, "-n", "--cv-per-branch")) {
+ opt.policy = arb::cv_policy_fixed_per_branch(nseg.value());
+ }
+ else {
+ usage(argv[0], "[-m|--morphology SWCFILE] [-d|--dt TIME] [-t|--t-end TIME] [-w|--weight WEIGHT] [-n|--cv-per-branch N]");
+ std::exit(1);
+ }
+ }
+ return opt;
+}
+
+// If no SWC file is given, the default morphology consists
+// of a soma of radius 6.3 µm and a single unbranched dendrite
+// of length 200 µm and radius decreasing linearly from 0.5 µm
+// to 0.2 µm.
+
+arb::morphology default_morphology() {
+ arb::segment_tree tree;
+
+ tree.append(arb::mnpos, { -6.3, 0.0, 0.0, 6.3}, { 6.3, 0.0, 0.0, 6.3}, 1);
+ tree.append( 0, { 6.3, 0.0, 0.0, 0.5}, {206.3, 0.0, 0.0, 0.2}, 3);
+
+ return arb::morphology(tree);
+}
+
+arb::morphology read_swc(const std::string& path) {
+ std::ifstream f(path);
+ if (!f) throw std::runtime_error("unable to open SWC file: "+path);
+
+ return arborio::load_swc_arbor(arborio::parse_swc(f));
+}
diff --git a/mechanisms/CMakeLists.txt b/mechanisms/CMakeLists.txt
index 56c0e208..320fd852 100644
--- a/mechanisms/CMakeLists.txt
+++ b/mechanisms/CMakeLists.txt
@@ -17,7 +17,7 @@ make_catalogue(
make_catalogue(
NAME default
- MOD exp2syn expsyn expsyn_curr expsyn_stdp hh kamt kdrmt nax nernst pas gj decay inject
+ MOD exp2syn expsyn expsyn_curr expsyn_stdp hh kamt kdrmt nax nernst pas gj decay inject v_clamp v_limit
VERBOSE ${ARB_CAT_VERBOSE}
ADD_DEPS ON)
diff --git a/mechanisms/default/v_clamp.mod b/mechanisms/default/v_clamp.mod
new file mode 100644
index 00000000..33e3eb24
--- /dev/null
+++ b/mechanisms/default/v_clamp.mod
@@ -0,0 +1,12 @@
+NEURON {
+ VOLTAGE_PROCESS v_clamp
+ GLOBAL v0
+}
+
+PARAMETER {
+ v0 = -60 (mV)
+}
+
+BREAKPOINT {
+ v = v0
+}
diff --git a/mechanisms/default/v_limit.mod b/mechanisms/default/v_limit.mod
new file mode 100644
index 00000000..276f625f
--- /dev/null
+++ b/mechanisms/default/v_limit.mod
@@ -0,0 +1,13 @@
+NEURON {
+ VOLTAGE_PROCESS v_limit
+ GLOBAL v_low, v_high
+}
+
+PARAMETER {
+ v_high = 20 (mV)
+ v_low = -70 (mV)
+}
+
+BREAKPOINT {
+ v = max(min(v, v_high), v_low)
+}
diff --git a/modcc/blocks.cpp b/modcc/blocks.cpp
index 2bc80734..70005455 100644
--- a/modcc/blocks.cpp
+++ b/modcc/blocks.cpp
@@ -36,7 +36,15 @@ ARB_LIBMODCC_API std::ostream& operator<<(std::ostream& os, IonDep const& I) {
}
ARB_LIBMODCC_API std::ostream& operator<<(std::ostream& os, moduleKind const& k) {
- return os << (k==moduleKind::density ? "density" : "point process");
+ switch (k) {
+ case moduleKind::density: os << "density"; break;
+ case moduleKind::voltage: os << "voltage"; break;
+ case moduleKind::point: os << "point"; break;
+ case moduleKind::revpot: os << "revpot"; break;
+ case moduleKind::junction: os << "junction"; break;
+ default: os << "???"; break;
+ }
+ return os;
}
ARB_LIBMODCC_API std::ostream& operator<<(std::ostream& os, NeuronBlock const& N) {
diff --git a/modcc/identifier.hpp b/modcc/identifier.hpp
index 47c835f0..1da95302 100644
--- a/modcc/identifier.hpp
+++ b/modcc/identifier.hpp
@@ -5,7 +5,7 @@
#include <stdexcept>
enum class moduleKind {
- point, density, junction, revpot
+ point, density, junction, revpot, voltage,
};
/// indicate how a variable is accessed
diff --git a/modcc/mechanism.hpp b/modcc/mechanism.hpp
deleted file mode 100644
index 81518165..00000000
--- a/modcc/mechanism.hpp
+++ /dev/null
@@ -1,43 +0,0 @@
-#pragma once
-
-#include "lib/vector/include/Vector.h"
-
-/*
- Abstract base class for all mechanisms
-
- This works well for the standard interface that is exported by all mechanisms,
- i.e. compute_currents(), etc. The overhead of using virtual dispatch
- for such functions is negligable compared to the cost of the operations themselves.
- However, the friction between compile time and run time dispatch has to be considered
- carefully:
- - we want to dispatch on template parameters like vector width, target
- hardware etc. Maybe these could be template parameters to the base
- class?
- - how to expose mechanism functionality, e.g. for computing a mechanism-
- specific quantity for visualization?
-*/
-
-class Mechanism {
-public:
- // typedefs for storage
- using value_type = double;
- using array = memory::HostVector<value_type>;
- using view_type = memory::HostView<value_type>;
-
- Mechanism(std::string const& name)
- : name_(name)
- {}
-
- //virtual void state() = 0;
- //virtual void jacobi() = 0;
- virtual void current() = 0;
- virtual void init() = 0;
-
- std::string const& name() {
- return name_;
- }
-
-protected:
- std::string name_;
-};
-
diff --git a/modcc/module.cpp b/modcc/module.cpp
index 6b761162..dcde114f 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -604,6 +604,9 @@ bool Module::semantic() {
post_events_api.first->body(post_events_api.second->body()->clone());
}
+ // check voltage mechanisms before rev pot ... otherwise we are in trouble
+ check_voltage_mechanism();
+
// Are we writing an ionic reversal potential? If so, change the moduleKind to
// `revpot` and assert that the mechanism is 'pure': it has no state variables;
// it contributes to no currents, ionic or otherwise; it isn't a point mechanism;
@@ -651,12 +654,13 @@ void Module::add_variables_to_symbols() {
return symbols_[token.spelling] = make_symbol<WhiteNoise>(token.location, token.spelling);
};
- sourceKind current_kind = kind_==moduleKind::density? sourceKind::current_density: sourceKind::current;
- sourceKind conductance_kind = kind_==moduleKind::density? sourceKind::conductivity: sourceKind::conductance;
+ auto current_kind = kind_ == moduleKind::density ? sourceKind::current_density : sourceKind::current;
+ auto conductance_kind = kind_ == moduleKind::density ? sourceKind::conductivity : sourceKind::conductance;
+ auto v_access = kind_ == moduleKind::voltage ? accessKind::write : accessKind::read;
create_indexed_variable("current_", current_kind, accessKind::write, "", Location());
create_indexed_variable("conductivity_", conductance_kind, accessKind::write, "", Location());
- create_indexed_variable("v", sourceKind::voltage, accessKind::read, "", Location());
+ create_indexed_variable("v", sourceKind::voltage, v_access, "", Location());
create_indexed_variable("v_peer", sourceKind::peer_voltage, accessKind::read, "", Location());
create_indexed_variable("dt", sourceKind::dt, accessKind::read, "", Location());
@@ -1010,3 +1014,20 @@ void Module::check_revpot_mechanism() {
kind_ = moduleKind::revpot;
}
+
+void Module::check_voltage_mechanism() {
+ if (kind_ != moduleKind::voltage) return;
+
+ auto impure = neuron_block_.has_nonspecific_current();
+ for (const auto& ion: neuron_block_.ions) {
+ impure |= ion.writes_concentration_int();
+ impure |= ion.writes_concentration_ext();
+ impure |= ion.writes_current();
+ impure |= ion.writes_rev_potential();
+ }
+
+ if (impure) {
+ error("Voltage mechanisms may not write ionic quantities..");
+ return;
+ }
+}
diff --git a/modcc/module.hpp b/modcc/module.hpp
index b75897d5..81facaa1 100644
--- a/modcc/module.hpp
+++ b/modcc/module.hpp
@@ -150,6 +150,9 @@ private:
return s == symbols_.end() ? false : s->second->kind() == kind;
}
+ // Check requirements for membrane potential manipulation.
+ void check_voltage_mechanism();
+
// Check requirements for reversal potential setters.
void check_revpot_mechanism();
diff --git a/modcc/parser.cpp b/modcc/parser.cpp
index cc86179f..dfd0bf4d 100644
--- a/modcc/parser.cpp
+++ b/modcc/parser.cpp
@@ -236,17 +236,19 @@ void Parser::parse_neuron_block() {
case tok::suffix:
case tok::point_process:
+ case tok::voltage_process:
case tok::junction_process:
neuron_block.kind = (token_.type == tok::suffix) ? moduleKind::density :
+ (token_.type == tok::voltage_process) ? moduleKind::voltage :
(token_.type == tok::point_process) ? moduleKind::point : moduleKind::junction;
- // set the modul kind
+ // set the module kind
module_->kind(neuron_block.kind);
get_token(); // consume SUFFIX / POINT_PROCESS
// assert that a valid name for the Neuron has been specified
if (token_.type != tok::identifier) {
- error(pprintf("invalid name for SUFFIX, found '%'", token_.spelling));
+ error(pprintf("invalid name for mechanism, found '%'", token_.spelling));
return;
}
neuron_block.name = token_.spelling;
diff --git a/modcc/printer/cprinter.cpp b/modcc/printer/cprinter.cpp
index 432c1187..708ad24f 100644
--- a/modcc/printer/cprinter.cpp
+++ b/modcc/printer/cprinter.cpp
@@ -232,7 +232,8 @@ ARB_LIBMODCC_API std::string emit_cpp_source(const Module& module_, const printe
auto emit_body = [&](APIMethod *p, bool add=false) {
auto flags = ApiFlags{}
.additive(add)
- .point(moduleKind::point == module_.kind());
+ .point(moduleKind::point == module_.kind())
+ .voltage(moduleKind::voltage == module_.kind());
if (with_simd) {
emit_simd_api_body(out, p, vars.scalars, flags);
} else {
@@ -536,19 +537,26 @@ void emit_api_body(std::ostream& out, APIMethod* method, const ApiFlags& flags)
for (auto& sym: indexed_vars) {
auto d = decode_indexed_variable(sym->external_variable());
+ auto write_voltage = sym->external_variable()->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
out << "arb_value_type " << cprint(sym) << " = ";
- if (sym->is_read() || (sym->is_write() && d.additive)) {
+ if (sym->is_read() || (sym->is_write() && d.additive) || write_voltage) {
out << scaled(d.scale) << deref(d) << ";\n";
}
else {
out << "0;\n";
}
}
+
out << cprint(body);
for (auto& sym: indexed_vars) {
if (!sym->external_variable()->is_write()) continue;
auto d = decode_indexed_variable(sym->external_variable());
+ auto write_voltage = sym->external_variable()->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
+ if (write_voltage) d.readonly = false;
+
bool use_weight = d.always_use_weight || !flags.is_point;
if (d.readonly) throw compiler_exception("Cannot assign to read-only external state: "+sym->to_string());
std::string
@@ -564,6 +572,10 @@ void emit_api_body(std::ostream& out, APIMethod* method, const ApiFlags& flags)
"{0} = fma({1}{2}, {3}, {0});\n",
var, scale, weight, name);
}
+ else if (write_voltage) {
+ // SAFETY: we only ever allow *one* V-PROCESS per CV, so this is OK.
+ out << fmt::format("{0} = {1};\n", var, name);
+ }
else if (d.accumulate) {
out << fmt::format("{} = fma({}{}, {}, {});\n",
var, scale, weight, name, var);
@@ -705,11 +717,15 @@ void SimdPrinter::visit(BlockExpression* block) {
EXITM(out_, "block");
}
-void emit_simd_state_read(std::ostream& out, LocalVariable* local, simd_expr_constraint constraint) {
+void emit_simd_state_read(std::ostream& out, LocalVariable* local, simd_expr_constraint constraint, const ApiFlags& flags) {
ENTER(out);
out << "simd_value " << local->name();
- if (local->is_read() || (local->is_write() && decode_indexed_variable(local->external_variable()).additive)) {
+ auto write_voltage = local->external_variable()->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
+ auto is_additive = local->is_write() && decode_indexed_variable(local->external_variable()).additive;
+
+ if (local->is_read() || is_additive || write_voltage) {
auto d = decode_indexed_variable(local->external_variable());
if (d.scalar()) {
out << " = simd_cast<simd_value>(" << pp_var_pfx << d.data_var
@@ -760,9 +776,12 @@ void emit_simd_state_update(std::ostream& out,
const ApiFlags& flags) {
if (!external->is_write()) return;
ENTER(out);
-
auto d = decode_indexed_variable(external);
+ if (external->data_source() == sourceKind::voltage && flags.can_write_voltage) {
+ d.readonly = false;
+ }
+
if (d.readonly) {
throw compiler_exception("Cannot assign to read-only external state: "+external->to_string());
}
@@ -779,6 +798,9 @@ void emit_simd_state_update(std::ostream& out,
ss << "S::mul(" << name << ", simd_cast<simd_value>(" << as_c_double(1/d.scale) << "))";
scaled = ss.str();
}
+ auto write_voltage = external->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
+ if (write_voltage) d.readonly = false;
std::string weight = (d.always_use_weight || !flags.is_point) ? "w_" : "simd_cast<simd_value>(1.0)";
@@ -807,6 +829,34 @@ void emit_simd_state_update(std::ostream& out,
data, index, weight, scaled);
}
}
+ else if (write_voltage) {
+ /* For voltage processes we *assign* to the potential field.
+ ** SAFETY: only one V-PROCESS per CV allowed
+ */
+ if (d.index_var_kind == index_kind::node) {
+ if (constraint == simd_expr_constraint::contiguous) {
+ out << fmt::format("indirect({} + {}, simd_width_) = {};\n",
+ data, node, name);
+ }
+ else {
+ // We need this instead of simple assignment!
+ out << fmt::format("{{\n"
+ " simd_value t_{}0_ = simd_cast<simd_value>(0.0);\n"
+ " assign(t_{}0_, indirect({}, simd_cast<simd_index>({}), simd_width_, constraint_category_));\n"
+ " {} -= t_{}0_;\n"
+ " indirect({}, simd_cast<simd_index>({}), simd_width_, constraint_category_) += S::mul({}, {});\n"
+ "}}\n",
+ name,
+ name, data, node,
+ scaled, name,
+ data, node, weight, scaled);
+ }
+ }
+ else {
+ out << fmt::format("indirect({}, {}, simd_width_, index_constraint::none) = {};\n",
+ data, index, name);
+ }
+ }
else if (d.accumulate) {
if (d.index_var_kind == index_kind::node) {
std::string tempvar = "t_" + external->name();
@@ -903,7 +953,7 @@ void emit_simd_body_for_loop(
emit_simd_index_initialize(out, indices, constraint);
for (auto& sym: indexed_vars) {
- emit_simd_state_read(out, sym, constraint);
+ emit_simd_state_read(out, sym, constraint, flags);
}
simdprint printer(body, scalars);
@@ -988,7 +1038,7 @@ void emit_simd_api_body(std::ostream& out, APIMethod* method,
else {
// We may nonetheless need to read a global scalar indexed variable.
for (auto& sym: scalar_indexed_vars) {
- emit_simd_state_read(out, sym, simd_expr_constraint::other);
+ emit_simd_state_read(out, sym, simd_expr_constraint::other, flags);
}
out << fmt::format("for (arb_size_type i_ = 0; i_ < {}width; i_ += simd_width_) {{\n",
diff --git a/modcc/printer/cprinter.hpp b/modcc/printer/cprinter.hpp
index b53c318b..010d1398 100644
--- a/modcc/printer/cprinter.hpp
+++ b/modcc/printer/cprinter.hpp
@@ -51,15 +51,17 @@ struct ApiFlags {
bool ppack_iface=true;
bool use_additive=false;
bool is_point=false;
+ bool can_write_voltage=false;
ApiFlags& loop(bool v) { cv_loop = v; return *this; }
ApiFlags& iface(bool v) { ppack_iface = v; return *this; }
ApiFlags& additive(bool v) { use_additive = v; return *this; }
ApiFlags& point(bool v) { is_point = v; return *this; }
+ ApiFlags& voltage(bool v) { can_write_voltage = v; return *this; }
};
-const ApiFlags net_recv_flags = {false, false, true, false};
-const ApiFlags post_evt_flags = {false, false, false, false};
+const ApiFlags net_recv_flags = {false, false, true}; // No CV loop, no PPACK, use additive
+const ApiFlags post_evt_flags = {false, false}; // No CV loop, no PPACK
class ARB_LIBMODCC_API SimdPrinter: public Visitor {
public:
diff --git a/modcc/printer/gpuprinter.cpp b/modcc/printer/gpuprinter.cpp
index 7b326fea..2150bf02 100644
--- a/modcc/printer/gpuprinter.cpp
+++ b/modcc/printer/gpuprinter.cpp
@@ -29,7 +29,7 @@ static std::string scaled(double coeff) {
void emit_api_body_cu(std::ostream& out, APIMethod* method, const ApiFlags&);
-void emit_state_read_cu(std::ostream& out, LocalVariable* local);
+void emit_state_read_cu(std::ostream& out, LocalVariable* local, const ApiFlags&);
void emit_state_update_cu(std::ostream& out, Symbol* from, IndexedVariable* external, const ApiFlags&);
const char* index_id(Symbol *s);
@@ -200,7 +200,7 @@ ARB_LIBMODCC_API std::string emit_gpu_cu_source(const Module& module_, const pri
<< "void " << e->name() << "(arb_mechanism_ppack params_) {\n" << indent
<< "int n_ = params_.width;\n"
<< "int tid_ = threadIdx.x + blockDim.x*blockIdx.x;\n";
- emit_api_body_cu(out, e, ApiFlags{}.point(is_point_proc).additive(additive));
+ emit_api_body_cu(out, e, ApiFlags{}.point(is_point_proc).additive(additive).voltage(moduleKind::voltage == module_.kind()));
out << popindent << "}\n\n";
}
};
@@ -411,7 +411,7 @@ void emit_api_body_cu(std::ostream& out, APIMethod* e, const ApiFlags& flags) {
}
for (auto& sym: indexed_vars) {
- emit_state_read_cu(out, sym);
+ emit_state_read_cu(out, sym, flags);
}
out << cuprint(body);
@@ -438,10 +438,12 @@ namespace {
};
}
-void emit_state_read_cu(std::ostream& out, LocalVariable* local) {
+void emit_state_read_cu(std::ostream& out, LocalVariable* local, const ApiFlags& flags) {
+ auto write_voltage = local->external_variable()->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
out << "arb_value_type " << cuprint(local) << " = ";
auto d = decode_indexed_variable(local->external_variable());
- if (local->is_read() || (local->is_write() && d.additive)) {
+ if (local->is_read() || (local->is_write() && d.additive) || write_voltage) {
if (d.scale != 1) {
out << as_c_double(d.scale) << "*";
}
@@ -453,10 +455,15 @@ void emit_state_read_cu(std::ostream& out, LocalVariable* local) {
}
-void emit_state_update_cu(std::ostream& out, Symbol* from,
- IndexedVariable* external, const ApiFlags& flags) {
+void emit_state_update_cu(std::ostream& out,
+ Symbol* from,
+ IndexedVariable* external,
+ const ApiFlags& flags) {
if (!external->is_write()) return;
auto d = decode_indexed_variable(external);
+ auto write_voltage = external->data_source() == sourceKind::voltage
+ && flags.can_write_voltage;
+ if (write_voltage) d.readonly = false;
if (d.readonly) {
throw compiler_exception("Cannot assign to read-only external state: "+external->to_string());
}
@@ -466,8 +473,8 @@ void emit_state_update_cu(std::ostream& out, Symbol* from,
auto data = pp_var_pfx + d.data_var;
auto index = index_i_name(d.outer_index_var());
auto var = deref(d);
- std::string weight = (d.always_use_weight || !flags.is_point) ? pp_var_pfx + "weight[tid_]" : "1.0";
- weight = scale + weight;
+ auto use_weight = d.always_use_weight || !flags.is_point;
+ std::string weight = scale + (use_weight ? pp_var_pfx + "weight[tid_]" : "1.0");
if (d.additive && flags.use_additive) {
out << name << " -= " << var << ";\n";
@@ -478,6 +485,15 @@ void emit_state_update_cu(std::ostream& out, Symbol* from,
out << var << " = fma(" << weight << ", " << name << ", " << var << ");\n";
}
}
+ else if (write_voltage) {
+ /* SAFETY:
+ ** - Only one V-PROCESS per CV
+ ** - these can never be point mechs
+ ** - they run separatly from density/point mechs
+ */
+ out << name << " -= " << var << ";\n"
+ << var << " = fma(" << weight << ", " << name << ", " << var << ");\n";
+ }
else if (d.accumulate) {
if (flags.is_point) {
out << "::arb::gpu::reduce_by_key(" << weight << "*" << name << ',' << data << ", " << index << ", lane_mask_);\n";
diff --git a/modcc/printer/printerutil.hpp b/modcc/printer/printerutil.hpp
index f0a7ff8e..7fc8cfa3 100644
--- a/modcc/printer/printerutil.hpp
+++ b/modcc/printer/printerutil.hpp
@@ -60,6 +60,7 @@ struct namespace_declaration_close {
inline const char* module_kind_str(const Module& m) {
switch (m.kind()) {
case moduleKind::density: return "arb_mechanism_kind_density";
+ case moduleKind::voltage: return "arb_mechanism_kind_voltage";
case moduleKind::point: return "arb_mechanism_kind_point";
case moduleKind::revpot: return "arb_mechanism_kind_reversal_potential";
case moduleKind::junction: return "arb_mechanism_kind_gap_junction";
diff --git a/modcc/token.cpp b/modcc/token.cpp
index e8b70a5a..34081c48 100644
--- a/modcc/token.cpp
+++ b/modcc/token.cpp
@@ -41,6 +41,7 @@ static Keyword keywords[] = {
{"UNITSOFF", tok::unitsoff},
{"UNITSON", tok::unitson},
{"SUFFIX", tok::suffix},
+ {"VOLTAGE_PROCESS", tok::voltage_process},
{"NONSPECIFIC_CURRENT", tok::nonspecific_current},
{"USEION", tok::useion},
{"READ", tok::read},
@@ -132,6 +133,7 @@ static TokenString token_strings[] = {
{"UNITSOFF", tok::unitsoff},
{"UNITSON", tok::unitson},
{"SUFFIX", tok::suffix},
+ {"VOLTAGE_PROCESS", tok::voltage_process},
{"NONSPECIFIC_CURRENT", tok::nonspecific_current},
{"USEION", tok::useion},
{"READ", tok::read},
diff --git a/modcc/token.hpp b/modcc/token.hpp
index c3e6d316..109d38a0 100644
--- a/modcc/token.hpp
+++ b/modcc/token.hpp
@@ -64,7 +64,7 @@ enum class tok {
range, local, conserve, compartment,
solve, method, steadystate,
threadsafe, global,
- point_process, junction_process,
+ point_process, junction_process, voltage_process,
from, to,
// prefix binary operators
diff --git a/python/cells.cpp b/python/cells.cpp
index b3d542c4..bc881d3a 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -466,8 +466,6 @@ void register_cells(pybind11::module& m) {
.def(pybind11::init([](const std::string& i, double v) -> arb::ion_diffusivity { return {i, v}; }))
.def("__repr__", [](const arb::ion_diffusivity& d){return "D" + d.ion + "=" + std::to_string(d.value);});
- // arb::density
-
pybind11::class_<arb::density> density(m, "density", "For painting a density mechanism on a region.");
density
.def(pybind11::init([](const std::string& name) {return arb::density(name);}))
@@ -478,6 +476,16 @@ void register_cells(pybind11::module& m) {
.def("__repr__", [](const arb::density& d){return "<arbor.density " + mechanism_desc_str(d.mech) + ">";})
.def("__str__", [](const arb::density& d){return "<arbor.density " + mechanism_desc_str(d.mech) + ">";});
+ pybind11::class_<arb::voltage_process> voltage_process(m, "voltage_process", "For painting a voltage_process mechanism on a region.");
+ voltage_process
+ .def(pybind11::init([](const std::string& name) {return arb::voltage_process(name);}))
+ .def(pybind11::init([](arb::mechanism_desc mech) {return arb::voltage_process(mech);}))
+ .def(pybind11::init([](const std::string& name, const std::unordered_map<std::string, double>& params) {return arb::voltage_process(name, params);}))
+ .def(pybind11::init([](arb::mechanism_desc mech, const std::unordered_map<std::string, double>& params) {return arb::voltage_process(mech, params);}))
+ .def_readonly("mech", &arb::voltage_process::mech, "The underlying mechanism.")
+ .def("__repr__", [](const arb::voltage_process& d){return "<arbor.voltage_process " + mechanism_desc_str(d.mech) + ">";})
+ .def("__str__", [](const arb::voltage_process& d){return "<arbor.voltage_process " + mechanism_desc_str(d.mech) + ">";});
+
// arb::scaled_mechanism<arb::density>
pybind11::class_<arb::scaled_mechanism<arb::density>> scaled_mechanism(
@@ -828,6 +836,12 @@ void register_cells(pybind11::module& m) {
},
"region"_a, "mechanism"_a,
"Associate a density mechanism with a region.")
+ .def("paint",
+ [](arb::decor& dec, const char* region, const arb::voltage_process& mechanism) {
+ return dec.paint(arborio::parse_region_expression(region).unwrap(), mechanism);
+ },
+ "region"_a, "mechanism"_a,
+ "Associate a voltage process mechanism with a region.")
.def("paint",
[](arb::decor& dec, const char* region, const arb::scaled_mechanism<arb::density>& mechanism) {
dec.paint(arborio::parse_region_expression(region).unwrap(), mechanism);
diff --git a/python/example/v-clamp.py b/python/example/v-clamp.py
new file mode 100755
index 00000000..27469177
--- /dev/null
+++ b/python/example/v-clamp.py
@@ -0,0 +1,55 @@
+#!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
+
+import arbor
+import pandas # You may have to pip install these.
+import seaborn # You may have to pip install these.
+
+# (1) Create a morphology with a single (cylindrical) segment of length=diameter=6 μm
+tree = arbor.segment_tree()
+tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)
+
+# (2) Define the soma and its midpoint
+labels = arbor.label_dict({"soma": "(tag 1)", "midpoint": "(location 0 0.5)"})
+
+# (3) Create and set up a decor object
+
+decor = (
+ arbor.decor()
+ .set_property(Vm=-40)
+ .paint('"soma"', arbor.density("hh"))
+ .paint('"soma"', arbor.voltage_process("v_clamp/v0=-42"))
+ .place('"midpoint"', arbor.iclamp(10, 2, 0.8), "iclamp")
+ .place('"midpoint"', arbor.threshold_detector(-10), "detector")
+)
+
+# (4) Create cell and the single cell model based on it
+cell = arbor.cable_cell(tree, decor, labels)
+
+# (5) Make single cell model.
+m = arbor.single_cell_model(cell)
+
+# (6) Attach voltage probe sampling at 10 kHz (every 0.1 ms).
+m.probe("voltage", '"midpoint"', frequency=10)
+
+# (7) Run simulation for 30 ms of simulated activity.
+m.run(tfinal=30)
+
+# (8) Print spike times.
+if len(m.spikes) > 0:
+ print("{} spikes:".format(len(m.spikes)))
+ for s in m.spikes:
+ print("{:3.3f}".format(s))
+else:
+ print("no spikes")
+
+# (9) Plot the recorded voltages over time.
+print("Plotting results ...")
+seaborn.set_theme() # Apply some styling to the plot
+df = pandas.DataFrame({"t/ms": m.traces[0].time, "U/mV": m.traces[0].value})
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV", ci=None).savefig(
+ "single_cell_model_result.svg"
+)
+
+# (10) Optionally, you can store your results for later processing.
+df.to_csv("single_cell_model_result.csv", float_format="%g")
diff --git a/scripts/run_cpp_examples.sh b/scripts/run_cpp_examples.sh
index 5af6591c..f332e60b 100755
--- a/scripts/run_cpp_examples.sh
+++ b/scripts/run_cpp_examples.sh
@@ -26,7 +26,7 @@ check () {
fi
}
-for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity ou
+for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity ou voltage-clamp
do
echo " - $ex"
dir=`echo $ex | tr ' ' '_'`
diff --git a/scripts/run_python_examples.sh b/scripts/run_python_examples.sh
index 4277d7b7..3bdc5a7d 100755
--- a/scripts/run_python_examples.sh
+++ b/scripts/run_python_examples.sh
@@ -31,3 +31,4 @@ $PREFIX python python/example/network_ring_gpu.py # by default, gpu_id=None
$PREFIX python python/example/network_two_cells_gap_junctions.py
$PREFIX python python/example/diffusion.py
$PREFIX python python/example/plasticity.py
+$PREFIX python python/example/v-clamp.py
diff --git a/test/unit-modcc/mod_files/test_v_process.mod b/test/unit-modcc/mod_files/test_v_process.mod
new file mode 100644
index 00000000..a770c55b
--- /dev/null
+++ b/test/unit-modcc/mod_files/test_v_process.mod
@@ -0,0 +1,13 @@
+NEURON {
+ VOLTAGE_PROCESS vmech
+}
+
+INITIAL { c = 0 }
+
+STATE { c }
+
+PARAMETER {
+ tau_c = 150 (ms)
+}
+
+BREAKPOINT { v = 42 }
diff --git a/test/unit-modcc/test_module.cpp b/test/unit-modcc/test_module.cpp
index 67e4b32f..daf0d1b8 100644
--- a/test/unit-modcc/test_module.cpp
+++ b/test/unit-modcc/test_module.cpp
@@ -123,6 +123,15 @@ TEST(Module, read_write_ion) {
EXPECT_TRUE(m.semantic());
}
+TEST(Module, v_process) {
+ Module m(io::read_all(DATADIR "/mod_files/test_v_process.mod"), "test_v_process.mod");
+ EXPECT_NE(m.buffer().size(), 0u);
+
+ Parser p(m, false);
+ EXPECT_TRUE(p.parse());
+ EXPECT_TRUE(m.semantic());
+}
+
// Regression test in #1893 we found that the solver segfaults when handed a
// naked comparison statement.
TEST(Module, solver_bug_1893) {
diff --git a/test/unit-modcc/test_parser.cpp b/test/unit-modcc/test_parser.cpp
index c666a3b4..a5800ce7 100644
--- a/test/unit-modcc/test_parser.cpp
+++ b/test/unit-modcc/test_parser.cpp
@@ -79,6 +79,16 @@ TEST(Parser, full_file) {
EXPECT_EQ(p.status(), lexerStatus::happy);
}
+TEST(Parser, v_proc) {
+ Module m(io::read_all(DATADIR "/mod_files/test_v_process.mod"), "v_process.mod");
+ if (m.buffer().size() == 0) {
+ std::cout << "skipping Parser.full_file test because unable to open input file" << std::endl;
+ return;
+ }
+ Parser p(m);
+ EXPECT_EQ(p.status(), lexerStatus::happy);
+}
+
TEST(Parser, procedure) {
std::vector<const char*> calls = {
"PROCEDURE foo(x, y) {\n"
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 36892930..f42737eb 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -145,6 +145,7 @@ set(unit_sources
test_unique_any.cpp
test_vector.cpp
test_version.cpp
+ test_v_clamp.cpp
# unit test driver
test.cpp
diff --git a/test/unit/test_abi.cpp b/test/unit/test_abi.cpp
index 1c6c68b5..7c700ac9 100644
--- a/test/unit/test_abi.cpp
+++ b/test/unit/test_abi.cpp
@@ -59,7 +59,7 @@ TEST(abi, multicore_initialisation) {
layout.weight.assign(ncv, 1.);
for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
- EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+ EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout, {}));
{
ASSERT_EQ(globals.size(), mech.mech_.n_globals);
@@ -137,7 +137,7 @@ TEST(abi, multicore_null) {
layout.weight.assign(ncv, 1.);
for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
- EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+ EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout, {}));
}
#ifdef ARB_GPU_ENABLED
@@ -202,7 +202,7 @@ TEST(abi, gpu_initialisation) {
layout.weight.assign(ncv, 1.);
for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
- EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+ EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout, {}));
{
ASSERT_EQ(globals.size(), mech.mech_.n_globals);
@@ -279,7 +279,7 @@ TEST(abi, gpu_null) {
layout.weight.assign(ncv, 1.);
for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
- EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+ EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout, {}));
}
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index a11be575..351504fc 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -608,8 +608,8 @@ TEST(fvm_lowered, ionic_concentrations) {
ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, read_cai_mech->data_alignment());
shared_state->add_ion("ca", 2, ion_config);
- shared_state->instantiate(*read_cai_mech, 0, overrides, layout);
- shared_state->instantiate(*write_cai_mech, 1, overrides, layout);
+ shared_state->instantiate(*read_cai_mech, 0, overrides, layout, {});
+ shared_state->instantiate(*write_cai_mech, 1, overrides, layout, {});
shared_state->reset();
diff --git a/test/unit/test_kinetic_linear.cpp b/test/unit/test_kinetic_linear.cpp
index e61d59a9..bccdf44e 100644
--- a/test/unit/test_kinetic_linear.cpp
+++ b/test/unit/test_kinetic_linear.cpp
@@ -56,7 +56,7 @@ void run_test(std::string mech_name,
layout.cv.push_back(i);
}
- shared_state->instantiate(*test, 0, overrides, layout);
+ shared_state->instantiate(*test, 0, overrides, layout, {});
shared_state->reset();
test->initialize();
diff --git a/test/unit/test_mech_temp_diam.cpp b/test/unit/test_mech_temp_diam.cpp
index d3adf7ad..dbaf1b0a 100644
--- a/test/unit/test_mech_temp_diam.cpp
+++ b/test/unit/test_mech_temp_diam.cpp
@@ -46,7 +46,7 @@ void run_celsius_test() {
layout.cv.push_back(i);
}
- shared_state->instantiate(*celsius_test, 0, overrides, layout);
+ shared_state->instantiate(*celsius_test, 0, overrides, layout, {});
shared_state->reset();
// expect 0 value in state 'c' after init:
@@ -96,7 +96,7 @@ void run_diam_test() {
ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, celsius_test->data_alignment());
- shared_state->instantiate(*celsius_test, 0, overrides, layout);
+ shared_state->instantiate(*celsius_test, 0, overrides, layout, {});
shared_state->reset();
// expect 0 value in state 'd' after init:
diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp
index 98fe3e7c..f1e92378 100644
--- a/test/unit/test_s_expr.cpp
+++ b/test/unit/test_s_expr.cpp
@@ -597,6 +597,9 @@ std::ostream& operator<<(std::ostream& o, const synapse& p) {
std::ostream& operator<<(std::ostream& o, const density& p) {
return o << "(density " << p.mech << ')';
}
+std::ostream& operator<<(std::ostream& o, const voltage_process& p) {
+ return o << "(voltage-process " << p.mech << ')';
+}
template <typename TaggedMech>
std::ostream& operator<<(std::ostream& o, const scaled_mechanism<TaggedMech>& p) {
o << "(scaled-mechanism " << p.t_mech;
@@ -718,6 +721,7 @@ TEST(decor_literals, round_tripping) {
"(ion-external-concentration \"h\" -50.1)",
"(ion-reversal-potential \"na\" 30)"};
auto paint_literals = {
+ "(voltage-process (mechanism \"hh\"))",
"(density (mechanism \"hh\"))",
"(density (mechanism \"pas\" (\"g\" 0.02)))",
"(scaled-mechanism (density (mechanism \"pas\" (\"g\" 0.02))))",
diff --git a/test/unit/test_segment_tree.cpp b/test/unit/test_segment_tree.cpp
index 571f332f..fb05e9ab 100644
--- a/test/unit/test_segment_tree.cpp
+++ b/test/unit/test_segment_tree.cpp
@@ -358,7 +358,7 @@ TEST(segment_tree, tag_roots) {
** \
** 2 - 3
** \
- ** 5
+ ** 5
*/
{
arb::segment_tree tree;
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 3fb0d713..d13aa718 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -107,8 +107,8 @@ TEST(synapses, syn_basic_state) {
std::vector<index_type> syn_mult(num_syn, 1);
std::vector<value_type> syn_weight(num_syn, 1.0);
- state.instantiate(*expsyn, 0, {}, {syn_cv, {}, syn_weight, syn_mult});
- state.instantiate(*exp2syn, 1, {}, {syn_cv, {}, syn_weight, syn_mult});
+ state.instantiate(*expsyn, 0, {}, {syn_cv, {}, syn_weight, syn_mult}, {});
+ state.instantiate(*exp2syn, 1, {}, {syn_cv, {}, syn_weight, syn_mult}, {});
// Parameters initialized to default values?
diff --git a/test/unit/test_v_clamp.cpp b/test/unit/test_v_clamp.cpp
new file mode 100644
index 00000000..a2948c7c
--- /dev/null
+++ b/test/unit/test_v_clamp.cpp
@@ -0,0 +1,276 @@
+#include <gtest/gtest.h>
+
+#include "common.hpp"
+
+#include <arbor/arbexcept.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/domain_decomposition.hpp>
+#include <arbor/load_balance.hpp>
+#include <arbor/recipe.hpp>
+#include <arbor/schedule.hpp>
+#include <arbor/simulation.hpp>
+
+#include <arborio/label_parse.hpp>
+using namespace arborio::literals;
+
+struct recipe: public arb::recipe {
+ recipe(bool clamp, bool limit): limit{limit}, clamp{clamp} {
+ gprop.default_parameters = arb::neuron_parameter_defaults;
+ gprop.default_parameters.discretization = arb::cv_policy_max_extent(1.0);
+ }
+
+ arb::cell_size_type num_cells() const override { return 1; }
+ arb::cell_kind get_cell_kind(arb::cell_gid_type gid) const override { return arb::cell_kind::cable; }
+ std::vector<arb::cell_connection> connections_on(arb::cell_gid_type gid) const override { return {}; }
+ arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override {
+ auto tree = arb::segment_tree{};
+ auto p = arb::mnpos;
+ p = tree.append(p, {0, 0, 0, 1 }, {0, 0, 1, 1 }, 1); // soma 0-1
+ p = tree.append(p, {0, 0, 1, 0.1}, {0, 0, 4, 0.1}, 2); // dend 1-4
+ auto decor = arb::decor{}
+ .paint("(tag 1)"_reg, arb::density("hh"))
+ .paint("(tag 2)"_reg, arb::density("pas"))
+ .place("(location 0 0.125)"_ls, arb::synapse("expsyn"), "tgt");
+
+ if (clamp) decor.paint("(tag 1)"_reg, arb::voltage_process("v_clamp/v0=-42"));
+ if (limit) decor.paint("(tag 1)"_reg, arb::voltage_process("v_limit/v_high=0,v_low=-60"));
+ return arb::cable_cell(arb::morphology{tree}, decor);
+ }
+ std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override {
+ return { arb::cable_probe_membrane_voltage{"(location 0 0.125)"_ls}, // soma center: 0.25/2
+ arb::cable_probe_membrane_voltage{"(location 0 0.625)"_ls}}; // dend center: 0.75/2 + 0.25
+ }
+ std::vector<arb::event_generator> event_generators(arb::cell_gid_type) const override {
+ return {arb::regular_generator({"tgt"}, 5.0, 0.2, 0.05)};
+ }
+ std::any get_global_properties(arb::cell_kind) const override { return gprop; }
+
+ arb::cable_cell_global_properties gprop;
+ bool limit = false;
+ bool clamp = false;
+};
+
+struct Um_type {
+ constexpr static double delta = 1e-6;
+
+ double t;
+ double u;
+
+ friend std::ostream& operator<<(std::ostream& os, const Um_type& um) {
+ os << "{ " << um.t << ", " << um.u << " }";
+ return os;
+ }
+
+ friend bool operator==(const Um_type& lhs, const Um_type& rhs) {
+ return (std::abs(lhs.t - rhs.t) <= delta)
+ && (std::abs(lhs.u - rhs.u) <= delta);
+ }
+};
+
+using um_s_type = std::vector<Um_type>;
+
+TEST(v_process, clamp) {
+ auto u_soma = um_s_type{};
+ auto u_dend = um_s_type{};
+ auto fun = [&u_soma, &u_dend](arb::probe_metadata pm,
+ std::size_t n,
+ const arb::sample_record* samples) {
+ for (std::size_t ix = 0ul; ix < n; ++ix) {
+ const auto& [t, v] = samples[ix];
+ double u = *arb::util::any_cast<const double*>(v);
+ if (pm.id.index == 0) {
+ u_soma.push_back({t, u});
+ }
+ else if (pm.id.index == 1) {
+ u_dend.push_back({t, u});
+ }
+ else {
+ throw std::runtime_error{"Unexpected probe id"};
+ }
+ }
+ };
+ auto sim = arb::simulation(recipe{true, false});
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun);
+ sim.run(1.0, 0.005);
+
+ um_s_type exp_soma{{ 0, -65 },
+ { 0.05, -42 },
+ { 0.095, -42 },
+ { 0.145, -42 },
+ { 0.2, -42 },
+ { 0.25, -42 },
+ { 0.3, -42 },
+ { 0.35, -42 },
+ { 0.4, -42 },
+ { 0.45, -42 },
+ { 0.5, -42 },
+ { 0.55, -42 },
+ { 0.6, -42 },
+ { 0.65, -42 },
+ { 0.7, -42 },
+ { 0.75, -42 },
+ { 0.8, -42 },
+ { 0.85, -42 },
+ { 0.9, -42 },
+ { 0.95, -42 },};
+ um_s_type exp_dend{{ 0, -65 },
+ { 0.05, -42.1167152 },
+ { 0.095, -42.0917893 },
+ { 0.145, -42.0464582 },
+ { 0.2, -41.9766685 },
+ { 0.25, -0.124773816 },
+ { 0.3, -0.0708536802 },
+ { 0.35, -0.0526604489 },
+ { 0.4, -0.043518448 },
+ { 0.45, -0.0380607556 },
+ { 0.5, -0.0344769712 },
+ { 0.55, -0.0319770579 },
+ { 0.6, -0.0301575597 },
+ { 0.65, -0.0287897554 },
+ { 0.7, -0.0277342043 },
+ { 0.75, -0.0269013197 },
+ { 0.8, -0.0262312483 },
+ { 0.85, -0.0256827719 },
+ { 0.9, -0.0252268065 },
+ { 0.95, -0.0248424087 }};
+ ASSERT_TRUE(testing::seq_eq(u_soma, exp_soma));
+ ASSERT_TRUE(testing::seq_eq(u_dend, exp_dend));
+}
+
+TEST(v_process, limit) {
+ auto u_soma = um_s_type{};
+ auto u_dend = um_s_type{};
+ auto fun = [&u_soma, &u_dend](arb::probe_metadata pm,
+ std::size_t n,
+ const arb::sample_record* samples) {
+ for (std::size_t ix = 0ul; ix < n; ++ix) {
+ const auto& [t, v] = samples[ix];
+ double u = *arb::util::any_cast<const double*>(v);
+ if (pm.id.index == 0) {
+ u_soma.push_back({t, u});
+ }
+ else if (pm.id.index == 1) {
+ u_dend.push_back({t, u});
+ }
+ else {
+ throw std::runtime_error{"Unexpected probe id"};
+ }
+ }
+ };
+ auto sim = arb::simulation(recipe{false, true});
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun);
+ sim.run(1.0, 0.005);
+
+ um_s_type exp_soma{{ 0, -65 },
+ { 0.05, -60 },
+ { 0.095, -60 },
+ { 0.145, -60 },
+ { 0.2, -60 },
+ { 0.25, -0.000425679283 },
+ { 0.3, 0 },
+ { 0.35, 0 },
+ { 0.4, 0 },
+ { 0.45, 0 },
+ { 0.5, 0 },
+ { 0.55, 0 },
+ { 0.6, 0 },
+ { 0.65, 0 },
+ { 0.7, 0 },
+ { 0.75, 0 },
+ { 0.8, 0 },
+ { 0.85, 0 },
+ { 0.9, 0 },
+ { 0.95, 0 }};
+ um_s_type exp_dend{{ 0, -65 },
+ { 0.05, -60.032677 },
+ { 0.095, -60.0308171 },
+ { 0.145, -60.028791 },
+ { 0.2, -60.0266703 },
+ { 0.25, -0.0178456839 },
+ { 0.3, -0.0168966758 },
+ { 0.35, -0.0162935748 },
+ { 0.4, -0.0158862702 },
+ { 0.45, -0.0156348508 },
+ { 0.5, -0.0155055671 },
+ { 0.55, -0.0154673288 },
+ { 0.6, -0.0154936751 },
+ { 0.65, -0.0155635334 },
+ { 0.7, -0.0156609133 },
+ { 0.75, -0.015774124 },
+ { 0.8, -0.0158948866 },
+ { 0.85, -0.0160175185 },
+ { 0.9, -0.0161382535 },
+ { 0.95, -0.0162547064 },};
+ ASSERT_TRUE(testing::seq_eq(u_soma, exp_soma));
+ ASSERT_TRUE(testing::seq_eq(u_dend, exp_dend));
+}
+
+TEST(v_process, clamp_fine) {
+ auto u_soma = um_s_type{};
+ auto u_dend = um_s_type{};
+ auto fun = [&u_soma, &u_dend](arb::probe_metadata pm,
+ std::size_t n,
+ const arb::sample_record* samples) {
+ for (std::size_t ix = 0ul; ix < n; ++ix) {
+ const auto& [t, v] = samples[ix];
+ double u = *arb::util::any_cast<const double*>(v);
+ if (pm.id.index == 0) {
+ u_soma.push_back({t, u});
+ }
+ else if (pm.id.index == 1) {
+ u_dend.push_back({t, u});
+ }
+ else {
+ throw std::runtime_error{"Unexpected probe id"};
+ }
+ }
+ };
+ auto rec = recipe{true, false};
+ rec.gprop.default_parameters.discretization = arb::cv_policy_max_extent(0.5);
+ auto sim = arb::simulation(rec);
+ sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun);
+ sim.run(1.0, 0.005);
+
+ um_s_type exp_soma{{ 0, -65 },
+ { 0.05, -42 },
+ { 0.095, -42 },
+ { 0.145, -42 },
+ { 0.2, -42 },
+ { 0.25, -42 },
+ { 0.3, -42 },
+ { 0.35, -42 },
+ { 0.4, -42 },
+ { 0.45, -42 },
+ { 0.5, -42 },
+ { 0.55, -42 },
+ { 0.6, -42 },
+ { 0.65, -42 },
+ { 0.7, -42 },
+ { 0.75, -42 },
+ { 0.8, -42 },
+ { 0.85, -42 },
+ { 0.9, -42 },
+ { 0.95, -42 },};
+ um_s_type exp_dend{{ 0, -65 },
+ { 0.05, -42.1164544 },
+ { 0.095, -42.0915289 },
+ { 0.145, -42.0461939 },
+ { 0.2, -41.9764002 },
+ { 0.25, -0.124099713 },
+ { 0.3, -0.0701885048 },
+ { 0.35, -0.0519983288 },
+ { 0.4, -0.0428578593 },
+ { 0.45, -0.0374010627 },
+ { 0.5, -0.0338178467 },
+ { 0.55, -0.0313183138 },
+ { 0.6, -0.0294990809 },
+ { 0.65, -0.0281314686 },
+ { 0.7, -0.0270760611 },
+ { 0.75, -0.0262432877 },
+ { 0.8, -0.025573305 },
+ { 0.85, -0.0250249014 },
+ { 0.9, -0.0245689972 },
+ { 0.95, -0.0241846519 },};
+ ASSERT_TRUE(testing::seq_eq(u_soma, exp_soma));
+ ASSERT_TRUE(testing::seq_eq(u_dend, exp_dend));
+}
--
GitLab