diff --git a/.gitignore b/.gitignore
index 036ba49182d024c4de5daf45f4026bf7184de43d..4b1dccf5bf06dd3dfcdfda2c86255f5167315b43 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 74dca037f3b2e655bf65c2517c0ee4a2e81c65fd..21ccedc18d8c3421701fe56b42bf656a118ce1fb 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 d6c66065a54b7064f346b6c366a6194858e83d4f..ab6ae51665e8e0dd8140a1cffde8b2d2bf212fe1 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 c09cdd12989275aaf53659ee8652d0f955919411..e08de3aece8b56e5c959f548f03c1aac885bb496 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 e7c2e8f8effcae528331ed5620cb137837c191d4..1ae389365cbe006603f225a9daf7e29a8f9c90bb 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 7d306ce7cf68ab36190fdd888b589129df0a127b..dda56f51eef8381e60d3849c40a27d945415963a 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 8b7d03823f13e0006d7ec6ade1682f8612efc24a..34c4f39b77e9d896064f6a3fb01778a005d719cb 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 94a8e1071f2e8c4f3eacd5089975aea4a3d2fb35..7862419dcce192a8dc31763abf61a979361c23cc 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 09902f5c27701e6da3b0d6cf8fe674431763a80a..b3d09765e2c343655d1eb7ba32ff237bd12c856a 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 deb578e43cbb7b5db61b95877488acae4864cffb..efa1255d0a6e95abb2595629b715a1e02776d42b 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 cb0570bc1a05cf2f728eba8b8dbda29cc9986dc7..221c967d149d8bd1b15067d0f34acd16e920e90f 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 b4307d0a843b040d4e6abac56b4c752c149cfd95..42748a57a3f4e35231431ebd2a006e0decff2013 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 17ab8d8eb4067492711d9a74e0e98816d2a7ffbb..52fb6eff4e116b67b86165b38f9e6a72dd5c2a06 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 a5050ed84c32994c7c2d57a83c403abe9a4df703..bebe1aaad38b51bca727bf4aaa23a828b9ae4d01 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 966c472f2e50e734c43e9ae3278a724c43168331..3c3ba9e7f95321b06d768530d6e10abf415e3a25 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 8be039c645c83c64183d3b945069b75099cd1570..0afbad2317cf5bde241a96c04a6bba3cf9217372 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 a125cc82d8d393ce7e48993ecd5aae95224bd7b7..f8015e40f41acee4216da1a77d3d40e47ca7e920 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 b341b25736cca5247a63fc58ba1d291a1a52bfb7..e97d41acab3e53459a55cca65ca466e1f325757c 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 d7f21181b8fb8ac45a76676cd447ff7de4216ae0..388832d7613707293c12bcfad72bf852a123c315 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 0000000000000000000000000000000000000000..43de3438c3407bc7d14331ff52878e345751c69d
--- /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 0000000000000000000000000000000000000000..d1a69fb92cb864d0657906825f8995175ad720c8
--- /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 0000000000000000000000000000000000000000..2d00a481121929fd66520668fbd59460e8adb7fb
--- /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 56c0e20836ca8c8ccbd3542067518925567e53d9..320fd85236a2d2ca89df94c9a8accda2eb601480 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 0000000000000000000000000000000000000000..33e3eb245f6bf03009f6f02b6604cec9d24feedf
--- /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 0000000000000000000000000000000000000000..276f625f7c2cea61e58fdb115e0749a150a34f88
--- /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 2bc80734030ec6b4e4f2434e0bc582f07e62b493..70005455bb8a7bc35c340027539a4069f0902f3a 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 47c835f0f977a4c7201ec6039779e28ae09fc79a..1da953023f45151a3391c8a18e4e99c685c88b38 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 8151816504cf82e429ce02ed66e424d24ca2c948..0000000000000000000000000000000000000000
--- 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 6b76116267a4f47c19f8fc27b67b01aa08fc5796..dcde114f668634d4e23bf43b32e738dd78770ef3 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 b75897d5ff1fd83b7906d591d151e486050b9f53..81facaa19ed8ca027053bc3d53fb6c8997e7ee3f 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 cc86179f5f63fc912c8d4ac3bdc09b28e70a6f98..dfd0bf4d9186fb845abdbc21229c620edbcccb3e 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 432c1187bf0967c1ff68ea868654a77007b7d0ca..708ad24fe124eba3893870cf2d2922d24c214e03 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 b53c318bd9ac73f16fd1299fe55486ec7136d5da..010d1398b00474ffca8732fbf723d0a34a3f81b4 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 7b326fea93b61d110efd547c89c3e52d1ccd376d..2150bf0291101802ecbb157851eb39bb5408da27 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 f0a7ff8ed748d8314cfb3356d01a220b4d68ff26..7fc8cfa34cfde32fb71348cb03ffa9b30cdb6b14 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 e8b70a5ae3e7c0afc18e78a9c8e998f03f323a7c..34081c48db414460ce4ec819270c1b9a1afa0d66 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 c3e6d3168539c4e65524b087ae1d3ba26f7af9b9..109d38a0624fc826350adb5928b9cb56a94d6def 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 b3d542c437df837eeb245cb480b1899dab018a95..bc881d3ae81fbd724f79fd6a6ed57bfa451e168c 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 0000000000000000000000000000000000000000..274691777b4489c7c462badaedd6adb84bdb1b9d
--- /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 5af6591ce71cee0722a280d1f0b2b84881940ea7..f332e60b574e166344fe7e5cfd1de67d11364339 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 4277d7b71f8a63a91aa1be7b8d2fa40def3ad6cc..3bdc5a7d4c19e4efca9f2ac1a40bfda1b07b21b2 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 0000000000000000000000000000000000000000..a770c55bbf3f64f4c4d6efd6ae10f06497350aaf
--- /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 67e4b32f100ad8a23993ba971834ed9c0ea31325..daf0d1b8504b1f963c3cf494f811d5a515b122d7 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 c666a3b49890623daf7d9bafaefcbe04e22ff9be..a5800ce7e32d62d981b5256945a865107f619bfe 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 36892930db35fba373fdad991e068e76721febe3..f42737eb5467dcb635e6bcbcbacb3f155eb63ee1 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 1c6c68b59ecc2d9c17d8440165f6b283b1d017c5..7c700ac90a09500803919e3b88163a89164e856d 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 a11be575dd7ad71ebdc6d6f7369e7f172f56715e..351504fce1ad2871953893c1343a9de1177a9755 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 e61d59a9cea9cfa04a616551396db0a6ae7d9ec5..bccdf44e2834dae58c755c4bfc621606fe5b61de 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 d3adf7ad0c3cc16b798922c68a39792a00ba89dd..dbaf1b0a124e3a2807c3fa512f6c56d13038c656 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 98fe3e7ce81d1acf5fe4878a4ea9271d03b1d3a8..f1e9237833353ba3b38117a20db08d07205a47f1 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 571f332f001353191c766704e277bd362de74590..fb05e9abe5247a4314178e810d06d29288078d5d 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 3fb0d713f2f217558752a5423b5a569a327c2a4f..d13aa718b00583b850ef74967fd67dd7f85c97d3 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 0000000000000000000000000000000000000000..a2948c7cd1f045e02c1f62f356d3f74a884a9eb5
--- /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));
+}