From 257e625dd8bdd5f27a4ff417e48e6d38ea1d365d Mon Sep 17 00:00:00 2001
From: Simon Frasch <simon.frasch@cscs.ch>
Date: Thu, 30 Jun 2022 09:47:39 +0200
Subject: [PATCH] Inhomogeneous parameters (#1887)
Adds the ability to scale parameters of (density) mechanisms based on inhomogeneous properties along a cell.
Two new types are added:
- iexpr: An expression describing the scaling factor computation.
- scaled_mechanism: A wrapper struct around a mechanisms with iexpr attached to selected parameters.
Closes #1650
---
arbor/CMakeLists.txt | 1 +
arbor/cable_cell.cpp | 36 +-
arbor/cable_cell_param.cpp | 2 +
arbor/fvm_layout.cpp | 27 +-
arbor/iexpr.cpp | 600 +++++++++++++++++++++++
arbor/include/arbor/cable_cell.hpp | 17 +-
arbor/include/arbor/cable_cell_param.hpp | 20 +-
arbor/include/arbor/iexpr.hpp | 137 ++++++
arbor/include/arbor/morph/label_dict.hpp | 7 +
arbor/include/arbor/morph/mprovider.hpp | 4 +
arbor/morph/label_dict.cpp | 25 +-
arbor/morph/mprovider.cpp | 12 +-
arbor/util/pw_over_cable.hpp | 8 +-
arborio/cableio.cpp | 92 +++-
arborio/include/arborio/label_parse.hpp | 2 +
arborio/label_parse.cpp | 75 +++
arborio/parse_helpers.hpp | 74 +++
doc/concepts/decor.rst | 24 +-
doc/concepts/labels.rst | 167 ++++++-
doc/fileformat/cable_cell.rst | 61 ++-
doc/images/iexpr_distance.svg | 45 ++
doc/images/iexpr_interp.svg | 47 ++
doc/scripts/inputs.py | 10 +
doc/scripts/make_images.py | 13 +
python/cells.cpp | 43 ++
python/proxy.hpp | 40 +-
test/unit/CMakeLists.txt | 1 +
test/unit/test_iexpr.cpp | 465 ++++++++++++++++++
test/unit/test_s_expr.cpp | 141 ++++++
29 files changed, 2126 insertions(+), 70 deletions(-)
create mode 100644 arbor/iexpr.cpp
create mode 100644 arbor/include/arbor/iexpr.hpp
create mode 100644 doc/images/iexpr_distance.svg
create mode 100644 doc/images/iexpr_interp.svg
create mode 100644 test/unit/test_iexpr.cpp
diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 569a27ec..985767ad 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -21,6 +21,7 @@ set(arbor_sources
fvm_lowered_cell_impl.cpp
hardware/memory.cpp
hardware/power.cpp
+ iexpr.cpp
io/locked_ostream.cpp
io/serialize_hex.cpp
label_resolution.cpp
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index a8a3bcf6..9f6a523f 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -1,3 +1,4 @@
+#include <memory>
#include <sstream>
#include <unordered_map>
#include <variant>
@@ -16,9 +17,6 @@
namespace arb {
-using region_map = std::unordered_map<std::string, mcable_list>;
-using locset_map = std::unordered_map<std::string, mlocation_list>;
-
using value_type = cable_cell::value_type;
using index_type = cable_cell::index_type;
using size_type = cable_cell::size_type;
@@ -138,8 +136,9 @@ struct cable_cell_impl {
return region_map.get<T>();
}
- mcable_map<density>& get_region_map(const density& desc) {
- return region_map.get<density>()[desc.mech.name()];
+ mcable_map<std::pair<density, iexpr_map>> &
+ get_region_map(const density &desc) {
+ return region_map.get<density>()[desc.mech.name()];
}
mcable_map<init_int_concentration>& get_region_map(const init_int_concentration& init) {
@@ -158,8 +157,31 @@ struct cable_cell_impl {
return region_map.get<init_reversal_potential>()[init.ion];
}
- template <typename Property>
- void paint(const region& reg, const Property& prop) {
+ void paint(const region& reg, const density& prop) {
+ this->paint(reg, scaled_mechanism<density>(prop));
+ }
+
+ void paint(const region& reg, const scaled_mechanism<density>& prop) {
+ mextent cables = thingify(reg, provider);
+ auto& mm = get_region_map(prop.t_mech);
+
+ std::unordered_map<std::string, iexpr_ptr> im;
+ for (const auto& [fst, snd]: prop.scale_expr) {
+ im.insert_or_assign(fst, thingify(snd, provider));
+ }
+
+ for (const auto& c: cables) {
+ // Skip zero-length cables in extent:
+ if (c.prox_pos == c.dist_pos) continue;
+
+ if (!mm.insert(c, {prop.t_mech, im})) {
+ throw cable_cell_error(util::pprintf("cable {} overpaints", c));
+ }
+ }
+ }
+
+ template <typename TaggedMech>
+ void paint(const region& reg, const TaggedMech& prop) {
mextent cables = thingify(reg, provider);
auto& mm = get_region_map(prop);
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index 1c2a91e4..2f48b0a7 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -1,8 +1,10 @@
#include <cfloat>
#include <cmath>
+#include <memory>
#include <numeric>
#include <vector>
#include <variant>
+#include <tuple>
#include <arbor/cable_cell.hpp>
#include <arbor/cable_cell_param.hpp>
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index e862adad..20b4aa7e 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -5,9 +5,12 @@
#include <unordered_set>
#include <unordered_map>
#include <vector>
+#include <utility>
+#include <memory>
#include <arbor/arbexcept.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/iexpr.hpp>
#include <arbor/math.hpp>
#include <arbor/morph/mcable_map.hpp>
#include <arbor/morph/mprovider.hpp>
@@ -897,9 +900,9 @@ 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());
// Density mechanisms:
-
for (const auto& [name, cables]: assignments.get<density>()) {
mechanism_info info = catalogue[name];
@@ -925,20 +928,23 @@ fvm_mechanism_data fvm_build_mechanism_data(
}
mcable_map<double> support;
- std::vector<mcable_map<double>> param_maps;
+ std::vector<mcable_map<std::pair<double, iexpr_ptr>>> param_maps;
param_maps.resize(n_param);
- for (auto& on_cable: cables) {
- const auto& mech = on_cable.second.mech;
+ for (const auto& [cable, density_iexpr]: cables) {
+ const auto& mech = density_iexpr.first.mech;
+ const auto& scale_expr = density_iexpr.second;
+
verify_mechanism(info, mech);
- mcable cable = on_cable.first;
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);
+ auto scale_it = scale_expr.find(param_names[i]);
+ param_maps[i].insert(cable, {value, scale_it == scale_expr.end()
+ ? unit_scale : scale_it->second});
}
}
@@ -954,7 +960,12 @@ fvm_mechanism_data fvm_build_mechanism_data(
area += area_on_cable;
for (std::size_t i = 0; i<n_param; ++i) {
- param_on_cv[i] += embedding.integrate_area(c.branch, pw_over_cable(param_maps[i], c, 0.));
+ param_on_cv[i] += embedding.integrate_area(
+ 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);
+ }));
}
}
@@ -1261,7 +1272,7 @@ fvm_mechanism_data fvm_build_mechanism_data(
for (auto i: cv_order) {
const i_clamp& stim = stimuli[i].item;
auto cv = stimuli_cv[i];
- double cv_area_scale = 1000./D.cv_area[cv]; // constant scales from nA/µm² to A/m².
+ double cv_area_scale = 1000./D.cv_area[cv]; // constant scale from nA/µm² to A/m².
config.cv.push_back(cv);
config.frequency.push_back(stim.frequency);
diff --git a/arbor/iexpr.cpp b/arbor/iexpr.cpp
new file mode 100644
index 00000000..f2fe81ac
--- /dev/null
+++ b/arbor/iexpr.cpp
@@ -0,0 +1,600 @@
+// Implementations for inhomogeneous expressions.
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <optional>
+#include <stdexcept>
+#include <variant>
+#include <sstream>
+
+#include <arbor/arbexcept.hpp>
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/iexpr.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/util/any_visitor.hpp>
+#include <arbor/math.hpp>
+
+namespace arb {
+
+namespace iexpr_impl {
+namespace {
+
+msize_t common_parent_branch(msize_t branch_a, msize_t branch_b, const morphology& m) {
+ // Locations on different branches.
+ // Find first common parent branch. Branch id of parent is
+ // always smaller.
+ while (branch_a != branch_b) {
+ if (branch_b == mnpos || (branch_a != mnpos && branch_a > branch_b))
+ branch_a = m.branch_parent(branch_a);
+ else
+ branch_b = m.branch_parent(branch_b);
+ }
+
+ return branch_a;
+}
+
+// compute the distance between any two points on the same morphology
+double compute_distance(const mlocation& loc_a, const mlocation& loc_b, const mprovider& p) {
+ if (loc_a.branch == loc_b.branch) return std::abs(p.embedding().integrate_length(loc_a, loc_b));
+
+ // If mnpos, locations are on different sides of root. Take
+ // distance to root in this case. Otherwise, take distance to
+ // end of parent branch
+ const auto base_branch = common_parent_branch(loc_a.branch, loc_b.branch, p.morphology());
+ const auto base_loc = base_branch == mnpos ? mlocation{0, 0.0} : mlocation{base_branch, 1.0};
+
+ // compute distance to distal end of parent branch and add
+ // together
+ return std::abs(p.embedding().integrate_length(loc_a, base_loc)) +
+ std::abs(p.embedding().integrate_length(loc_b, base_loc));
+};
+
+// Compute the distance in proximal direction. Will return nullopt if loc_prox is not between origin
+// and loc_dist
+std::optional<double> compute_proximal_distance(const mlocation& loc_prox,
+ const mlocation& loc_dist,
+ const mprovider& p) {
+
+ // check order if on same branch
+ if (loc_prox.branch == loc_dist.branch && loc_prox.pos > loc_dist.pos) return std::nullopt;
+
+ // Special case root, for which no direction can be assumed. Always return the actual distance
+ // in this case.
+ if (loc_prox.pos == 0.0 && p.morphology().branch_parent(loc_prox.branch) == mnpos)
+ return p.embedding().integrate_length(loc_prox, loc_dist);
+
+ // check if loc_prox branch is in proximal direction from loc_dist
+ auto b = loc_dist.branch;
+ while (b > loc_prox.branch) {
+ b = p.morphology().branch_parent(b);
+ if (b == mnpos) return std::nullopt;
+ }
+ if (b != loc_prox.branch) return std::nullopt;
+
+ return p.embedding().integrate_length(loc_prox, loc_dist);
+};
+
+enum class direction { any, proximal, distal };
+
+// compute the minimum distance in the given direction from the given locations towards loc_eval.
+// Returs nullopt if loc_eval cannot be found in the direction.
+template <direction Dir>
+std::optional<double> distance_from_locations(
+ const std::variant<mlocation_list, mextent>& locations,
+ const mlocation& loc_eval,
+ const mprovider& p) {
+ return std::visit(
+ arb::util::overload(
+ [&](const mlocation_list& arg) -> std::optional<double> {
+ std::optional<double> min_dist, dist;
+ for (const auto& loc: arg) {
+ if constexpr (Dir == direction::proximal) {
+ dist = compute_proximal_distance(loc_eval, loc, p);
+ }
+ else if constexpr (Dir == direction::distal) {
+ dist = compute_proximal_distance(loc, loc_eval, p);
+ }
+ else {
+ dist = compute_distance(loc, loc_eval, p);
+ }
+ if (dist)
+ min_dist = std::min(
+ min_dist.value_or(std::numeric_limits<double>::max()), dist.value());
+ }
+ return min_dist;
+ },
+ [&](const mextent& arg) -> std::optional<double> {
+ std::optional<double> min_dist, dist;
+ for (const auto& c: arg) {
+ if (c.branch == loc_eval.branch && c.prox_pos < loc_eval.pos &&
+ c.dist_pos > loc_eval.pos)
+ return std::nullopt;
+ if constexpr (Dir == direction::proximal) {
+ dist = compute_proximal_distance(loc_eval, {c.branch, c.prox_pos}, p);
+ }
+ else if constexpr (Dir == direction::distal) {
+ dist = compute_proximal_distance({c.branch, c.dist_pos}, loc_eval, p);
+ }
+ else {
+ dist = std::min(compute_distance({c.branch, c.dist_pos}, loc_eval, p),
+ compute_distance({c.branch, c.prox_pos}, loc_eval, p));
+ }
+ if (dist)
+ min_dist = std::min(
+ min_dist.value_or(std::numeric_limits<double>::max()), dist.value());
+ }
+ return min_dist;
+ }),
+ locations);
+}
+
+struct scalar: public iexpr_interface {
+ scalar(double v): value(v) {}
+
+ double eval(const mprovider&, const mcable&) const override { return value; }
+
+ double value;
+};
+
+struct radius: public iexpr_interface {
+ radius(double s): scale(s) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ auto loc_eval = mlocation{c.branch, (c.dist_pos + c.prox_pos) / 2};
+ return scale * p.embedding().radius(loc_eval);
+ }
+
+ double scale;
+};
+
+struct distance: public iexpr_interface {
+ distance(double s, std::variant<mlocation_list, mextent> l):
+ scale(s),
+ locations(std::move(l)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ auto loc_eval = mlocation{c.branch, (c.dist_pos + c.prox_pos) / 2};
+
+ return scale *
+ distance_from_locations<direction::any>(locations, loc_eval, p).value_or(0.0);
+ }
+
+ double scale;
+ std::variant<mlocation_list, mextent> locations;
+};
+
+struct proximal_distance: public iexpr_interface {
+ proximal_distance(double s, std::variant<mlocation_list, mextent> l):
+ scale(s),
+ locations(std::move(l)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ auto loc_eval = mlocation{c.branch, (c.dist_pos + c.prox_pos) / 2};
+
+ return scale *
+ distance_from_locations<direction::proximal>(locations, loc_eval, p).value_or(0.0);
+ }
+
+ double scale;
+ std::variant<mlocation_list, mextent> locations;
+};
+
+struct distal_distance: public iexpr_interface {
+ distal_distance(double s, std::variant<mlocation_list, mextent> l):
+ scale(s),
+ locations(std::move(l)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ auto loc_eval = mlocation{c.branch, (c.dist_pos + c.prox_pos) / 2};
+
+ return scale *
+ distance_from_locations<direction::distal>(locations, loc_eval, p).value_or(0.0);
+ }
+
+ double scale;
+ std::variant<mlocation_list, mextent> locations;
+};
+
+struct interpolation: public iexpr_interface {
+ interpolation(double prox_value,
+ std::variant<mlocation_list, mextent> prox_list,
+ double dist_value,
+ std::variant<mlocation_list, mextent> dist_list):
+ prox_v(prox_value),
+ dist_v(dist_value),
+ prox_l(std::move(prox_list)),
+ dist_l(std::move(dist_list)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ auto loc_eval = mlocation{c.branch, (c.dist_pos + c.prox_pos) / 2};
+
+ const auto d1 = distance_from_locations<direction::distal>(prox_l, loc_eval, p);
+ if (!d1) return 0.0;
+
+ const auto d2 = distance_from_locations<direction::proximal>(dist_l, loc_eval, p);
+ if (!d2) return 0.0;
+
+ const auto sum = d1.value() + d2.value();
+ if (!sum) return (prox_v + dist_v) * 0.5;
+
+ return prox_v * (d2.value() / sum) + dist_v * (d1.value() / sum);
+ }
+
+ double prox_v, dist_v;
+ std::variant<mlocation_list, mextent> prox_l, dist_l;
+};
+
+struct add: public iexpr_interface {
+ add(iexpr_ptr l, iexpr_ptr r): left(std::move(l)), right(std::move(r)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return left->eval(p, c) + right->eval(p, c);
+ }
+
+ iexpr_ptr left;
+ iexpr_ptr right;
+};
+
+struct sub: public iexpr_interface {
+ sub(iexpr_ptr l, iexpr_ptr r): left(std::move(l)), right(std::move(r)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return left->eval(p, c) - right->eval(p, c);
+ }
+
+ iexpr_ptr left;
+ iexpr_ptr right;
+};
+
+struct mul: public iexpr_interface {
+ mul(iexpr_ptr l, iexpr_ptr r): left(std::move(l)), right(std::move(r)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return left->eval(p, c) * right->eval(p, c);
+ }
+
+ iexpr_ptr left;
+ iexpr_ptr right;
+};
+
+struct div: public iexpr_interface {
+ div(iexpr_ptr l, iexpr_ptr r): left(std::move(l)), right(std::move(r)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return left->eval(p, c) / right->eval(p, c);
+ }
+
+ iexpr_ptr left;
+ iexpr_ptr right;
+};
+
+struct exp: public iexpr_interface {
+ exp(iexpr_ptr v): value(std::move(v)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return std::exp(value->eval(p, c));
+ }
+
+ iexpr_ptr value;
+};
+
+struct log: public iexpr_interface {
+ log(iexpr_ptr v): value(std::move(v)) {}
+
+ double eval(const mprovider& p, const mcable& c) const override {
+ return std::log(value->eval(p, c));
+ }
+
+ iexpr_ptr value;
+};
+
+} // namespace
+} // namespace iexpr_impl
+
+iexpr::iexpr(double value) { *this = iexpr::scalar(value); }
+
+iexpr iexpr::scalar(double value) { return iexpr(iexpr_type::scalar, std::make_tuple(value)); }
+
+iexpr iexpr::pi() { return iexpr::scalar(math::pi<double>); }
+
+iexpr iexpr::distance(double scale, locset loc) {
+ return iexpr(
+ iexpr_type::distance, std::make_tuple(scale, std::variant<locset, region>(std::move(loc))));
+}
+
+iexpr iexpr::distance(locset loc) {
+ return iexpr::distance(1.0, std::move(loc));
+}
+
+iexpr iexpr::distance(double scale, region reg) {
+ return iexpr(
+ iexpr_type::distance, std::make_tuple(scale, std::variant<locset, region>(std::move(reg))));
+}
+
+iexpr iexpr::distance(region reg) {
+ return iexpr::distance(1.0, std::move(reg));
+}
+
+iexpr iexpr::proximal_distance(double scale, locset loc) {
+ return iexpr(iexpr_type::proximal_distance,
+ std::make_tuple(scale, std::variant<locset, region>(std::move(loc))));
+}
+
+iexpr iexpr::proximal_distance(locset loc) {
+ return iexpr::proximal_distance(1.0, std::move(loc));
+}
+
+iexpr iexpr::proximal_distance(double scale, region reg) {
+ return iexpr(iexpr_type::proximal_distance,
+ std::make_tuple(scale, std::variant<locset, region>(std::move(reg))));
+}
+
+iexpr iexpr::proximal_distance(region reg) {
+ return iexpr::proximal_distance(1.0, std::move(reg));
+}
+
+iexpr iexpr::distal_distance(double scale, locset loc) {
+ return iexpr(iexpr_type::distal_distance,
+ std::make_tuple(scale, std::variant<locset, region>(std::move(loc))));
+}
+
+iexpr iexpr::distal_distance(locset loc) {
+ return iexpr::distal_distance(1.0, std::move(loc));
+}
+
+iexpr iexpr::distal_distance(double scale, region reg) {
+ return iexpr(iexpr_type::distal_distance,
+ std::make_tuple(scale, std::variant<locset, region>(std::move(reg))));
+}
+
+iexpr iexpr::distal_distance(region reg) {
+ return iexpr::distal_distance(1.0, std::move(reg));
+}
+
+iexpr iexpr::interpolation(double prox_value,
+ locset prox_list,
+ double dist_value,
+ locset dist_list) {
+ return iexpr(iexpr_type::interpolation,
+ std::make_tuple(prox_value,
+ std::variant<locset, region>(std::move(prox_list)),
+ dist_value,
+ std::variant<locset, region>(std::move(dist_list))));
+}
+
+iexpr iexpr::interpolation(double prox_value,
+ region prox_list,
+ double dist_value,
+ region dist_list) {
+ return iexpr(iexpr_type::interpolation,
+ std::make_tuple(prox_value,
+ std::variant<locset, region>(std::move(prox_list)),
+ dist_value,
+ std::variant<locset, region>(std::move(dist_list))));
+}
+
+iexpr iexpr::radius(double scale) { return iexpr(iexpr_type::radius, std::make_tuple(scale)); }
+
+iexpr iexpr::radius() { return iexpr::radius(1.0); }
+
+iexpr iexpr::diameter(double scale) { return iexpr(iexpr_type::diameter, std::make_tuple(scale)); }
+
+iexpr iexpr::diameter() { return iexpr::diameter(1.0); }
+
+iexpr iexpr::add(iexpr left, iexpr right) {
+ return iexpr(iexpr_type::add, std::make_tuple(std::move(left), std::move(right)));
+}
+
+iexpr iexpr::sub(iexpr left, iexpr right) {
+ return iexpr(iexpr_type::sub, std::make_tuple(std::move(left), std::move(right)));
+}
+
+iexpr iexpr::mul(iexpr left, iexpr right) {
+ return iexpr(iexpr_type::mul, std::make_tuple(std::move(left), std::move(right)));
+}
+
+iexpr iexpr::div(iexpr left, iexpr right) {
+ return iexpr(iexpr_type::div, std::make_tuple(std::move(left), std::move(right)));
+}
+
+iexpr iexpr::exp(iexpr value) { return iexpr(iexpr_type::exp, std::make_tuple(std::move(value))); }
+
+iexpr iexpr::log(iexpr value) { return iexpr(iexpr_type::log, std::make_tuple(std::move(value))); }
+
+iexpr iexpr::named(std::string name) {
+ return iexpr(iexpr_type::named, std::make_tuple(std::move(name)));
+}
+
+iexpr_ptr thingify(const iexpr& expr, const mprovider& m) {
+ switch (expr.type()) {
+ case iexpr_type::scalar:
+ return iexpr_ptr(new iexpr_impl::scalar(
+ std::get<0>(std::any_cast<const std::tuple<double>&>(expr.args()))));
+ case iexpr_type::distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+
+ return std::visit(
+ [&](auto&& arg) {
+ return iexpr_ptr(new iexpr_impl::distance(scale, thingify(arg, m)));
+ },
+ var);
+ }
+ case iexpr_type::proximal_distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+
+ return std::visit(
+ [&](auto&& arg) {
+ return iexpr_ptr(new iexpr_impl::proximal_distance(scale, thingify(arg, m)));
+ },
+ var);
+ }
+ case iexpr_type::distal_distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(expr.args()));
+
+ return std::visit(
+ [&](auto&& arg) {
+ return iexpr_ptr(new iexpr_impl::distal_distance(scale, thingify(arg, m)));
+ },
+ var);
+ }
+ case iexpr_type::interpolation: {
+ const auto& t = std::any_cast<const std::
+ tuple<double, std::variant<locset, region>, double, std::variant<locset, region>>&>(
+ expr.args());
+ auto prox_list = std::visit(
+ [&](auto&& arg) -> std::variant<mlocation_list, mextent> { return thingify(arg, m); },
+ std::get<1>(t));
+
+ auto dist_list = std::visit(
+ [&](auto&& arg) -> std::variant<mlocation_list, mextent> { return thingify(arg, m); },
+ std::get<3>(t));
+ return iexpr_ptr(new iexpr_impl::interpolation(
+ std::get<0>(t), std::move(prox_list), std::get<2>(t), std::move(dist_list)));
+ }
+ case iexpr_type::radius:
+ return iexpr_ptr(new iexpr_impl::radius(
+ std::get<0>(std::any_cast<const std::tuple<double>&>(expr.args()))));
+ case iexpr_type::diameter:
+ return iexpr_ptr(new iexpr_impl::radius(
+ 2.0 * std::get<0>(std::any_cast<const std::tuple<double>&>(expr.args()))));
+ case iexpr_type::add:
+ return iexpr_ptr(new iexpr_impl::add(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m),
+ thingify(std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m)));
+ case iexpr_type::sub:
+ return iexpr_ptr(new iexpr_impl::sub(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m),
+ thingify(std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m)));
+ case iexpr_type::mul:
+ return iexpr_ptr(new iexpr_impl::mul(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m),
+ thingify(std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m)));
+ case iexpr_type::div:
+ return iexpr_ptr(new iexpr_impl::div(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m),
+ thingify(std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(expr.args())), m)));
+ case iexpr_type::exp:
+ return iexpr_ptr(new iexpr_impl::exp(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr>&>(expr.args())), m)));
+ case iexpr_type::log:
+ return iexpr_ptr(new iexpr_impl::log(
+ thingify(std::get<0>(std::any_cast<const std::tuple<iexpr>&>(expr.args())), m)));
+ case iexpr_type::named:
+ return m.iexpr(std::get<0>(std::any_cast<const std::tuple<std::string>&>(expr.args())));
+ }
+
+ throw std::runtime_error("thingify iexpr: Unknown iexpr type");
+ return nullptr;
+}
+
+std::ostream& operator<<(std::ostream& o, const iexpr& e) {
+ o << "(";
+
+ switch (e.type()) {
+ case iexpr_type::scalar: {
+ o << "scalar " << std::get<0>(std::any_cast<const std::tuple<double>&>(e.args()));
+ break;
+ }
+ case iexpr_type::distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ o << "distance " << scale << " ";
+
+ std::visit([&](auto&& arg) { o << arg; }, var);
+ break;
+ }
+ case iexpr_type::proximal_distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ o << "proximal-distance " << scale << " ";
+
+ std::visit([&](auto&& arg) { o << arg; }, var);
+ break;
+ }
+ case iexpr_type::distal_distance: {
+ const auto& scale = std::get<0>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ const auto& var = std::get<1>(
+ std::any_cast<const std::tuple<double, std::variant<locset, region>>&>(e.args()));
+ o << "distal-distance " << scale << " ";
+
+ std::visit([&](auto&& arg) { o << arg; }, var);
+ break;
+ }
+ case iexpr_type::interpolation: {
+ using arg_type =
+ std::tuple<double, std::variant<locset, region>, double, std::variant<locset, region>>;
+
+ o << "interpolation " << std::get<0>(std::any_cast<const arg_type&>(e.args())) << " ";
+ std::visit(
+ [&](auto&& arg) { o << arg; }, std::get<1>(std::any_cast<const arg_type&>(e.args())));
+ o << " " << std::get<2>(std::any_cast<const arg_type&>(e.args())) << " ";
+ std::visit(
+ [&](auto&& arg) { o << arg; }, std::get<3>(std::any_cast<const arg_type&>(e.args())));
+ break;
+ }
+ case iexpr_type::radius: {
+ o << "radius " << std::get<0>(std::any_cast<const std::tuple<double>&>(e.args()));
+ break;
+ }
+ case iexpr_type::diameter: {
+ o << "diameter " << std::get<0>(std::any_cast<const std::tuple<double>&>(e.args()));
+ break;
+ }
+ case iexpr_type::add: {
+ o << "add " << std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args())) << " "
+ << std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::sub: {
+ o << "sub " << std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args())) << " "
+ << std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::mul: {
+ o << "mul " << std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args())) << " "
+ << std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::div: {
+ o << "div " << std::get<0>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args())) << " "
+ << std::get<1>(std::any_cast<const std::tuple<iexpr, iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::exp: {
+ o << "exp " << std::get<0>(std::any_cast<const std::tuple<iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::log: {
+ o << "log " << std::get<0>(std::any_cast<const std::tuple<iexpr>&>(e.args()));
+ break;
+ }
+ case iexpr_type::named: {
+ o << "iexpr \"" << std::get<0>(std::any_cast<const std::tuple<std::string>&>(e.args()))
+ << "\"";
+ break;
+ }
+ default: throw std::runtime_error("print iexpr: Unknown iexpr type");
+ }
+
+ o << ")";
+ return o;
+}
+
+} // namespace arb
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index e2927c8f..6556deb0 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <memory>
#include <string>
#include <unordered_map>
#include <utility>
@@ -11,6 +12,7 @@
#include <arbor/cable_cell_param.hpp>
#include <arbor/common_types.hpp>
#include <arbor/constants.hpp>
+#include <arbor/iexpr.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/mcable_map.hpp>
@@ -195,20 +197,27 @@ struct ARB_SYMBOL_VISIBLE cable_probe_ion_ext_concentration_cell {
// Forward declare the implementation, for PIMPL.
struct cable_cell_impl;
-
// Typed maps for access to painted and placed assignments:
//
// Mechanisms and initial ion data are further keyed by
// mechanism name and ion name respectively.
+using iexpr_map = std::unordered_map<std::string, iexpr_ptr>;
+
template <typename T>
-using region_assignment =
+using region_assignment = std::conditional_t<
+ std::is_same_v<T, init_int_concentration> ||
+ std::is_same_v<T, init_ext_concentration> ||
+ std::is_same_v<T, init_reversal_potential>,
+ std::unordered_map<std::string, mcable_map<T>>,
std::conditional_t<
- std::is_same<T, density>::value || std::is_same<T, init_int_concentration>::value ||
+ 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::unordered_map<std::string, mcable_map<T>>,
- mcable_map<T>>;
+ std::conditional_t<std::is_same_v<T, density>,
+ std::unordered_map<std::string, mcable_map<std::pair<T, iexpr_map>>>,
+ mcable_map<T>>>>;
template <typename T>
struct placed {
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 0e8a78ef..2c6cca39 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -6,12 +6,15 @@
#include <unordered_map>
#include <string>
#include <variant>
+#include <any>
#include <arbor/export.hpp>
#include <arbor/arbexcept.hpp>
#include <arbor/cv_policy.hpp>
+#include <arbor/iexpr.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morph/locset.hpp>
+#include <arbor/morph/primitives.hpp>
namespace arb {
@@ -192,6 +195,7 @@ private:
std::unordered_map<std::string, double> param_;
};
+
// Tagged mechanism types for dispatching decor::place() and decor::paint() calls
struct ARB_SYMBOL_VISIBLE junction {
mechanism_desc mech;
@@ -228,6 +232,19 @@ struct ARB_SYMBOL_VISIBLE ion_reversal_potential_method {
mechanism_desc method;
};
+template <typename TaggedMech>
+struct ARB_SYMBOL_VISIBLE scaled_mechanism {
+ TaggedMech t_mech;
+ std::unordered_map<std::string, iexpr> scale_expr;
+
+ explicit scaled_mechanism(TaggedMech m) : t_mech(std::move(m)) {}
+
+ scaled_mechanism& scale(std::string name, iexpr expr) {
+ scale_expr.insert_or_assign(name, expr);
+ return *this;
+ }
+};
+
using paintable =
std::variant<init_membrane_potential,
axial_resistivity,
@@ -237,7 +254,8 @@ using paintable =
init_int_concentration,
init_ext_concentration,
init_reversal_potential,
- density>;
+ density,
+ scaled_mechanism<density>>;
using placeable =
std::variant<i_clamp,
diff --git a/arbor/include/arbor/iexpr.hpp b/arbor/include/arbor/iexpr.hpp
new file mode 100644
index 00000000..6921b777
--- /dev/null
+++ b/arbor/include/arbor/iexpr.hpp
@@ -0,0 +1,137 @@
+#pragma once
+
+// Implementations for inhomogeneous expressions.
+
+#include <any>
+#include <memory>
+#include <string>
+#include <ostream>
+
+#include <arbor/export.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/region.hpp>
+
+namespace arb {
+
+struct mprovider;
+
+enum class iexpr_type {
+ scalar,
+ distance,
+ proximal_distance,
+ distal_distance,
+ interpolation,
+ radius,
+ diameter,
+ add,
+ sub,
+ mul,
+ div,
+ exp,
+ log,
+ named
+};
+
+struct ARB_SYMBOL_VISIBLE iexpr {
+ // Convert to scalar expr type
+ iexpr(double value);
+
+ iexpr_type type() const { return type_; }
+
+ const std::any& args() const { return args_; }
+
+ static iexpr scalar(double value);
+
+ static iexpr pi();
+
+ static iexpr distance(double scale, locset loc);
+
+ static iexpr distance(locset loc);
+
+ static iexpr distance(double scale, region reg);
+
+ static iexpr distance(region reg);
+
+ static iexpr proximal_distance(double scale, locset loc);
+
+ static iexpr proximal_distance(locset loc);
+
+ static iexpr proximal_distance(double scale, region reg);
+
+ static iexpr proximal_distance(region reg);
+
+ static iexpr distal_distance(double scale, locset loc);
+
+ static iexpr distal_distance(locset loc);
+
+ static iexpr distal_distance(double scale, region reg);
+
+ static iexpr distal_distance(region reg);
+
+ static iexpr interpolation(double prox_value,
+ locset prox_list,
+ double dist_value,
+ locset dist_list);
+
+ static iexpr interpolation(double prox_value,
+ region prox_list,
+ double dist_value,
+ region dist_list);
+
+ static iexpr radius(double scale);
+
+ static iexpr radius();
+
+ static iexpr diameter(double scale);
+
+ static iexpr diameter();
+
+ static iexpr add(iexpr left, iexpr right);
+
+ static iexpr sub(iexpr left, iexpr right);
+
+ static iexpr mul(iexpr left, iexpr right);
+
+ static iexpr div(iexpr left, iexpr right);
+
+ static iexpr exp(iexpr value);
+
+ static iexpr log(iexpr value);
+
+ static iexpr named(std::string name);
+
+
+private:
+ iexpr(iexpr_type type, std::any args): type_(type), args_(std::move(args)) {}
+
+ iexpr_type type_;
+ std::any args_;
+};
+
+ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const iexpr& e);
+
+ARB_ARBOR_API inline iexpr operator+(iexpr a, iexpr b) { return iexpr::add(std::move(a), std::move(b)); }
+
+ARB_ARBOR_API inline iexpr operator-(iexpr a, iexpr b) { return iexpr::sub(std::move(a), std::move(b)); }
+
+ARB_ARBOR_API inline iexpr operator*(iexpr a, iexpr b) { return iexpr::mul(std::move(a), std::move(b)); }
+
+ARB_ARBOR_API inline iexpr operator/(iexpr a, iexpr b) { return iexpr::div(std::move(a), std::move(b)); }
+
+ARB_ARBOR_API inline iexpr operator+(iexpr a) { return a; }
+
+ARB_ARBOR_API inline iexpr operator-(iexpr a) { return iexpr::mul(-1.0, std::move(a)); }
+
+struct ARB_SYMBOL_VISIBLE iexpr_interface {
+
+ virtual double eval(const mprovider& p, const mcable& c) const = 0;
+
+ virtual ~iexpr_interface() = default;
+};
+
+using iexpr_ptr = std::shared_ptr<iexpr_interface>;
+
+ARB_ARBOR_API iexpr_ptr thingify(const iexpr& expr, const mprovider& m);
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/label_dict.hpp b/arbor/include/arbor/morph/label_dict.hpp
index b9d415e8..f7c43ec6 100644
--- a/arbor/include/arbor/morph/label_dict.hpp
+++ b/arbor/include/arbor/morph/label_dict.hpp
@@ -7,14 +7,18 @@
#include <arbor/export.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
+#include <arbor/iexpr.hpp>
namespace arb {
class ARB_ARBOR_API label_dict {
using ps_map = std::unordered_map<std::string, arb::locset>;
using reg_map = std::unordered_map<std::string, arb::region>;
+ using iexpr_map = std::unordered_map<std::string, arb::iexpr>;
+
ps_map locsets_;
reg_map regions_;
+ iexpr_map iexpressions_;
public:
// construct a label dict with SWC tags predefined
@@ -24,12 +28,15 @@ public:
label_dict& set(const std::string& name, locset ls);
label_dict& set(const std::string& name, region reg);
+ label_dict& set(const std::string& name, iexpr e);
std::optional<arb::region> region(const std::string& name) const;
std::optional<arb::locset> locset(const std::string& name) const;
+ std::optional<arb::iexpr> iexpr(const std::string& name) const;
const ps_map& locsets() const;
const reg_map& regions() const;
+ const iexpr_map& iexpressions() const;
std::size_t size() const;
};
diff --git a/arbor/include/arbor/morph/mprovider.hpp b/arbor/include/arbor/morph/mprovider.hpp
index 441ca2b5..40196c71 100644
--- a/arbor/include/arbor/morph/mprovider.hpp
+++ b/arbor/include/arbor/morph/mprovider.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <memory>
#include <string>
#include <unordered_map>
@@ -8,6 +9,7 @@
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/label_dict.hpp>
#include <arbor/util/expected.hpp>
+#include <arbor/iexpr.hpp>
namespace arb {
@@ -20,6 +22,7 @@ struct ARB_ARBOR_API mprovider {
// Throw exception on missing or recursive definition.
const mextent& region(const std::string& name) const;
const mlocation_list& locset(const std::string& name) const;
+ const iexpr_ptr& iexpr(const std::string& name) const;
// Read-only access to morphology and constructed embedding.
const auto& morphology() const { return morphology_; }
@@ -37,6 +40,7 @@ private:
// Maps are mutated only during initialization phase of mprovider.
mutable std::unordered_map<std::string, util::expected<mextent, circular_def>> regions_;
mutable std::unordered_map<std::string, util::expected<mlocation_list, circular_def>> locsets_;
+ mutable std::unordered_map<std::string, util::expected<iexpr_ptr, circular_def>> iexpressions_;
// Non-null only during initialization phase.
mutable const label_dict* label_dict_ptr;
diff --git a/arbor/morph/label_dict.cpp b/arbor/morph/label_dict.cpp
index bc6bca89..ba0c4d83 100644
--- a/arbor/morph/label_dict.cpp
+++ b/arbor/morph/label_dict.cpp
@@ -22,7 +22,7 @@ size_t label_dict::size() const {
}
label_dict& label_dict::set(const std::string& name, arb::locset ls) {
- if (regions_.count(name)) {
+ if (regions_.count(name) || iexpressions_.count(name)) {
throw label_type_mismatch(name);
}
locsets_[name] = std::move(ls);
@@ -30,13 +30,21 @@ label_dict& label_dict::set(const std::string& name, arb::locset ls) {
}
label_dict& label_dict::set(const std::string& name, arb::region reg) {
- if (locsets_.count(name)) {
+ if (locsets_.count(name) || iexpressions_.count(name)) {
throw label_type_mismatch(name);
}
regions_[name] = std::move(reg);
return *this;
}
+label_dict& label_dict::set(const std::string& name, arb::iexpr e) {
+ if (locsets_.count(name) || regions_.count(name)) {
+ throw label_type_mismatch(name);
+ }
+ iexpressions_.insert_or_assign(name, std::move(e));
+ return *this;
+}
+
void label_dict::import(const label_dict& other, const std::string& prefix) {
for (const auto& entry: other.locsets()) {
set(prefix+entry.first, entry.second);
@@ -44,6 +52,9 @@ void label_dict::import(const label_dict& other, const std::string& prefix) {
for (const auto& entry: other.regions()) {
set(prefix+entry.first, entry.second);
}
+ for (const auto& entry: other.iexpressions()) {
+ set(prefix+entry.first, entry.second);
+ }
}
std::optional<region> label_dict::region(const std::string& name) const {
@@ -58,6 +69,12 @@ std::optional<locset> label_dict::locset(const std::string& name) const {
return it->second;
}
+std::optional<iexpr> label_dict::iexpr(const std::string& name) const {
+ auto it = iexpressions_.find(name);
+ if (it==iexpressions_.end()) return std::nullopt;
+ return it->second;
+}
+
const std::unordered_map<std::string, locset>& label_dict::locsets() const {
return locsets_;
}
@@ -66,4 +83,8 @@ const std::unordered_map<std::string, region>& label_dict::regions() const {
return regions_;
}
+const std::unordered_map<std::string, iexpr>& label_dict::iexpressions() const {
+ return iexpressions_;
+}
+
} // namespace arb
diff --git a/arbor/morph/mprovider.cpp b/arbor/morph/mprovider.cpp
index 330fbf09..13082587 100644
--- a/arbor/morph/mprovider.cpp
+++ b/arbor/morph/mprovider.cpp
@@ -18,11 +18,15 @@ void mprovider::init() {
if (!label_dict_ptr) return;
for (const auto& pair: label_dict_ptr->regions()) {
- (void)region(pair.first);
+ (void)(this->region(pair.first));
}
for (const auto& pair: label_dict_ptr->locsets()) {
- (void)locset(pair.first);
+ (void)(this->locset(pair.first));
+ }
+
+ for (const auto& pair: label_dict_ptr->iexpressions()) {
+ (void)(this->iexpr(pair.first));
}
label_dict_ptr = nullptr;
@@ -72,5 +76,9 @@ const mlocation_list& mprovider::locset(const std::string& name) const {
return try_lookup(*this, name, locsets_, locsets_ptr);
}
+const iexpr_ptr& mprovider::iexpr(const std::string& name) const {
+ const auto* locsets_ptr = label_dict_ptr ? &(label_dict_ptr->iexpressions()) : nullptr;
+ return try_lookup(*this, name, iexpressions_, locsets_ptr);
+}
} // namespace arb
diff --git a/arbor/util/pw_over_cable.hpp b/arbor/util/pw_over_cable.hpp
index 9345b86b..cfab6db1 100644
--- a/arbor/util/pw_over_cable.hpp
+++ b/arbor/util/pw_over_cable.hpp
@@ -8,8 +8,8 @@ namespace util {
namespace impl {
struct get_value {
template <typename X>
- double operator()(const X& x) const { return x.value; }
- double operator()(double x) const { return x; }
+ double operator()(const mcable&, const X& x) const { return x.value; }
+ double operator()(const mcable&, double x) const { return x; }
};
} // namespace impl
@@ -39,7 +39,7 @@ util::pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt
if (el.first.prox_pos>pw_right) {
pw.push_back(pw_right, el.first.prox_pos, dflt_value);
}
- pw.push_back(el.first.prox_pos, el.first.dist_pos, projection(el.second));
+ pw.push_back(el.first.prox_pos, el.first.dist_pos, projection(el.first, el.second));
}
double pw_right = pw.empty()? 0: pw.bounds().second;
@@ -53,4 +53,4 @@ util::pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt
return pw;
}
} // namespace util
-} // namespace arb
\ No newline at end of file
+} // namespace arb
diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp
index 2727788b..103246eb 100644
--- a/arborio/cableio.cpp
+++ b/arborio/cableio.cpp
@@ -5,6 +5,8 @@
#include <arbor/s_expr.hpp>
#include <arbor/util/pp_util.hpp>
#include <arbor/util/any_visitor.hpp>
+#include <arbor/iexpr.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arborio/label_parse.hpp>
#include <arborio/cableio.hpp>
@@ -88,6 +90,20 @@ s_expr mksexp(const synapse& j) {
s_expr mksexp(const density& j) {
return slist("density"_symbol, mksexp(j.mech));
}
+s_expr mksexp(const iexpr& j) {
+ std::stringstream s;
+ s << j;
+ return parse_s_expr(s.str());
+}
+template <typename TaggedMech>
+s_expr mksexp(const scaled_mechanism<TaggedMech>& j) {
+ std::vector<s_expr> expressions;
+ std::transform(j.scale_expr.begin(),
+ j.scale_expr.end(),
+ std::back_inserter(expressions),
+ [](const auto& x) { return slist(s_expr(x.first), mksexp(x.second)); });
+ return s_expr{"scaled-mechanism"_symbol, s_expr{mksexp(j.t_mech), slist_range(expressions)}};
+}
s_expr mksexp(const mpoint& p) {
return slist("point"_symbol, p.x, p.y, p.z, p.radius);
}
@@ -134,6 +150,9 @@ s_expr mksexp(const label_dict& dict) {
for (auto& r: dict.regions()) {
defs = s_expr(slist("region-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs));
}
+ for (auto& r: dict.iexpressions()) {
+ defs = s_expr(slist("iexpr-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs));
+ }
return {"label-dict"_symbol, std::move(defs)};
}
s_expr mksexp(const morphology& morph) {
@@ -262,13 +281,17 @@ decor make_decor(const std::vector<std::variant<place_tuple, paint_pair, default
// Define maker for locset pairs, region pairs and label_dicts
using locset_pair = std::pair<std::string, locset>;
using region_pair = std::pair<std::string, region>;
-locset_pair make_locset_pair(const std::string& name, const locset& desc) {
- return locset_pair{name, desc};
+using iexpr_pair = std::pair<std::string, iexpr>;
+locset_pair make_locset_pair(std::string name, locset desc) {
+ return locset_pair{std::move(name), std::move(desc)};
+}
+region_pair make_region_pair(std::string name, region desc) {
+ return region_pair{std::move(name), std::move(desc)};
}
-region_pair make_region_pair(const std::string name, const region& desc) {
- return region_pair{name, desc};
+iexpr_pair make_iexpr_pair(std::string name, iexpr e) {
+ return iexpr_pair{std::move(name), std::move(e)};
}
-label_dict make_label_dict(const std::vector<std::variant<locset_pair, region_pair>>& args) {
+label_dict make_label_dict(const std::vector<std::variant<locset_pair, region_pair, iexpr_pair>>& args) {
label_dict d;
for(const auto& a: args) {
std::visit([&](auto&& x){d.set(x.first, x.second);}, a);
@@ -380,6 +403,48 @@ struct make_mech_call {
}
};
+
+// Test whether a list of arguments passed as a std::vector<std::any> can be converted
+// to a string followed by any number of std::pair<std::string, arb::iexpr>
+template <typename TaggedMech>
+struct scaled_mechanism_match {
+ bool operator()(const std::vector<std::any>& args) const {
+ // First argument is the mech name
+ if (args.size() < 1) return false;
+ if (!match<TaggedMech>(args.front().type())) return false;
+
+ // The rest of the arguments should be tuples
+ for (auto it = args.begin()+1; it != args.end(); ++it) {
+ if (!match<std::tuple<std::string, arb::iexpr>>(it->type())) return false;
+ }
+ return true;
+ }
+};
+// Create a mechanism_desc from a std::vector<std::any>.
+template <typename TaggedMech>
+struct scaled_mechanism_eval {
+ arb::scaled_mechanism<TaggedMech> operator()(const std::vector<std::any>& args) {
+ auto d = scaled_mechanism(eval_cast<TaggedMech>(args.front()));
+
+ for (auto it = args.begin() + 1; it != args.end(); ++it) {
+ auto p = eval_cast<std::tuple<std::string, arb::iexpr>>(*it);
+ d.scale(std::move(std::get<0>(p)), std::move(std::get<1>(p)));
+ }
+ return d;
+ }
+};
+// Wrap mech_match and mech_eval in an evaluator
+template <typename TaggedMech>
+struct make_scaled_mechanism_call {
+ evaluator state;
+ make_scaled_mechanism_call(const char* msg="scaled-mechanism"):
+ state(scaled_mechanism_eval<TaggedMech>(), scaled_mechanism_match<TaggedMech>(), msg)
+ {}
+ operator evaluator() const {
+ return state;
+ }
+};
+
// Test whether a list of arguments passed as a std::vector<std::any> can be converted
// to 2 integers followed by at least 1 msegment
struct branch_match {
@@ -554,6 +619,7 @@ parse_hopefully<std::any> eval(const s_expr& e, const eval_map& map, const eval_
if (auto l = parse_label_expression(e)) {
if (match<region>(l->type())) return eval_cast<region>(l.value());
if (match<locset>(l->type())) return eval_cast<locset>(l.value());
+ if (match<iexpr>(l->type())) return eval_cast<iexpr>(l.value());
}
// Or it could be a cv-policy expression
@@ -608,6 +674,8 @@ 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)")},
+ {"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)")},
{"place", make_call<locset, threshold_detector, std::string>(make_place, "'place' with 3 arguments (ls:locset t:threshold-detector name:string)")},
{"place", make_call<locset, junction, std::string>(make_place, "'place' with 3 arguments (ls:locset gj:junction name:string)")},
@@ -622,6 +690,7 @@ eval_map named_evals{
{"paint", make_call<region, ion_diffusivity>(make_paint, "'paint' with 2 arguments (reg:region v:ion-diffusivity)")},
{"paint", make_call<region, init_reversal_potential>(make_paint, "'paint' with 2 arguments (reg:region v:ion-reversal-potential)")},
{"paint", make_call<region, density>(make_paint, "'paint' with 2 arguments (reg:region v:density)")},
+ {"paint", make_call<region, scaled_mechanism<density>>(make_paint, "'paint' with 2 arguments (reg:region v:scaled-mechanism)")},
{"default", make_call<init_membrane_potential>(make_default, "'default' with 1 argument (v:membrane-potential)")},
{"default", make_call<temperature_K>(make_default, "'default' with 1 argument (v:temperature-kelvin)")},
@@ -638,6 +707,8 @@ eval_map named_evals{
"'locset-def' with 2 arguments (name:string ls:locset)")},
{"region-def", make_call<std::string, region>(make_region_pair,
"'region-def' with 2 arguments (name:string reg:region)")},
+ {"iexpr-def", make_call<std::string, iexpr>(make_iexpr_pair,
+ "'iexpr-def' with 2 arguments (name:string e:iexpr)")},
{"point", make_call<double, double, double, double>(make_point,
"'point' with 4 arguments (x:real y:real z:real radius:real)")},
@@ -648,8 +719,8 @@ eval_map named_evals{
{"decor", make_arg_vec_call<place_tuple, paint_pair, defaultable>(make_decor,
"'decor' with 1 or more `paint`, `place` or `default` arguments")},
- {"label-dict", make_arg_vec_call<locset_pair, region_pair>(make_label_dict,
- "'label-dict' with 1 or more `locset-def` or `region-def` arguments")},
+ {"label-dict", make_arg_vec_call<locset_pair, region_pair, iexpr_pair>(make_label_dict,
+ "'label-dict' with 1 or more `locset-def` or `region-def` or `iexpr-def` arguments")},
{"morphology", make_arg_vec_call<branch_tuple>(make_morphology,
"'morphology' 1 or more `branch` arguments")},
@@ -667,11 +738,12 @@ eval_map named_evals{
eval_vec unnamed_evals{
make_call<std::string, double>(std::make_tuple<std::string, double>, "tuple<std::string, double>"),
- make_call<double, double>(std::make_tuple<double, double>, "tuple<double, double>")
+ make_call<double, double>(std::make_tuple<double, double>, "tuple<double, double>"),
+ make_call<std::string, arb::iexpr>(std::make_tuple<std::string, arb::iexpr>, "tuple<std::string, arb::iexpr>"),
};
-inline parse_hopefully<std::any> parse(arb::s_expr s) {
- return eval(std::move(s), named_evals, unnamed_evals);
+inline parse_hopefully<std::any> parse(const arb::s_expr& s) {
+ return eval(s, named_evals, unnamed_evals);
}
ARB_ARBORIO_API parse_hopefully<std::any> parse_expression(const std::string& s) {
diff --git a/arborio/include/arborio/label_parse.hpp b/arborio/include/arborio/label_parse.hpp
index 5f96ede8..d937cad1 100644
--- a/arborio/include/arborio/label_parse.hpp
+++ b/arborio/include/arborio/label_parse.hpp
@@ -7,6 +7,7 @@
#include <arbor/morph/region.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/util/expected.hpp>
+#include <arbor/iexpr.hpp>
#include <arbor/s_expr.hpp>
#include <arborio/export.hpp>
@@ -26,6 +27,7 @@ ARB_ARBORIO_API parse_label_hopefully<std::any> parse_label_expression(const arb
ARB_ARBORIO_API parse_label_hopefully<arb::region> parse_region_expression(const std::string& s);
ARB_ARBORIO_API parse_label_hopefully<arb::locset> parse_locset_expression(const std::string& s);
+ARB_ARBORIO_API parse_label_hopefully<arb::iexpr> parse_iexpr_expression(const std::string& s);
namespace literals {
diff --git a/arborio/label_parse.cpp b/arborio/label_parse.cpp
index 961754cc..6700295c 100644
--- a/arborio/label_parse.cpp
+++ b/arborio/label_parse.cpp
@@ -6,6 +6,7 @@
#include <arbor/arbexcept.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
+#include <arbor/iexpr.hpp>
#include <arbor/util/expected.hpp>
@@ -112,6 +113,66 @@ std::unordered_multimap<std::string, evaluator> eval_map {
"'join' with at least 2 arguments: (locset locset [...locset])")},
{"sum", make_fold<arb::locset>(static_cast<arb::locset(*)(arb::locset, arb::locset)>(arb::sum),
"'sum' with at least 2 arguments: (locset locset [...locset])")},
+
+
+ // iexpr
+ {"iexpr", make_call<std::string>(arb::iexpr::named, "iexpr with 1 argument: (value:string)")},
+
+ {"scalar", make_call<double>(arb::iexpr::scalar, "iexpr with 1 argument: (value:double)")},
+
+ {"pi", make_call<>(arb::iexpr::pi, "iexpr with no argument")},
+
+ {"distance", make_call<double, arb::locset>(static_cast<arb::iexpr(*)(double, arb::locset)>(arb::iexpr::distance),
+ "iexpr with 2 arguments: (scale:double, loc:locset)")},
+ {"distance", make_call<arb::locset>(static_cast<arb::iexpr(*)(arb::locset)>(arb::iexpr::distance),
+ "iexpr with 1 argument: (loc:locset)")},
+ {"distance", make_call<double, arb::region>(static_cast<arb::iexpr(*)(double, arb::region)>(arb::iexpr::distance),
+ "iexpr with 2 arguments: (scale:double, reg:region)")},
+ {"distance", make_call<arb::region>(static_cast<arb::iexpr(*)(arb::region)>(arb::iexpr::distance),
+ "iexpr with 1 argument: (reg:region)")},
+
+ {"proximal-distance", make_call<double, arb::locset>(static_cast<arb::iexpr(*)(double, arb::locset)>(arb::iexpr::proximal_distance),
+ "iexpr with 2 arguments: (scale:double, loc:locset)")},
+ {"proximal-distance", make_call<arb::locset>(static_cast<arb::iexpr(*)(arb::locset)>(arb::iexpr::proximal_distance),
+ "iexpr with 1 argument: (loc:locset)")},
+ {"proximal-distance", make_call<double, arb::region>(static_cast<arb::iexpr(*)(double, arb::region)>(arb::iexpr::proximal_distance),
+ "iexpr with 2 arguments: (scale:double, reg:region)")},
+ {"proximal-distance", make_call<arb::region>(static_cast<arb::iexpr(*)(arb::region)>(arb::iexpr::proximal_distance),
+ "iexpr with 1 arguments: (reg:region)")},
+
+ {"distal-distance", make_call<double, arb::locset>(static_cast<arb::iexpr(*)(double, arb::locset)>(arb::iexpr::distal_distance),
+ "iexpr with 2 arguments: (scale:double, loc:locset)")},
+ {"distal-distance", make_call<arb::locset>(static_cast<arb::iexpr(*)(arb::locset)>(arb::iexpr::distal_distance),
+ "iexpr with 1 argument: (loc:locset)")},
+ {"distal-distance", make_call<double, arb::region>(static_cast<arb::iexpr(*)(double, arb::region)>(arb::iexpr::distal_distance),
+ "iexpr with 2 arguments: (scale:double, reg:region)")},
+ {"distal-distance", make_call<arb::region>(static_cast<arb::iexpr(*)(arb::region)>(arb::iexpr::distal_distance),
+ "iexpr with 1 argument: (reg:region)")},
+
+ {"interpolation", make_call<double, arb::locset, double, locset>(static_cast<arb::iexpr(*)(double, arb::locset, double, arb::locset)>(arb::iexpr::interpolation),
+ "iexpr with 4 arguments: (prox_value:double, prox_list:locset, dist_value:double, dist_list:locset)")},
+ {"interpolation", make_call<double, arb::region, double, region>(static_cast<arb::iexpr(*)(double, arb::region, double, arb::region)>(arb::iexpr::interpolation),
+ "iexpr with 4 arguments: (prox_value:double, prox_list:region, dist_value:double, dist_list:region)")},
+
+ {"radius", make_call<double>(static_cast<arb::iexpr(*)(double)>(arb::iexpr::radius), "iexpr with 1 argument: (value:double)")},
+ {"radius", make_call<>(static_cast<arb::iexpr(*)()>(arb::iexpr::radius), "iexpr with no argument")},
+
+ {"diameter", make_call<double>(static_cast<arb::iexpr(*)(double)>(arb::iexpr::diameter), "iexpr with 1 argument: (value:double)")},
+ {"diameter", make_call<>(static_cast<arb::iexpr(*)()>(arb::iexpr::diameter), "iexpr with no argument")},
+
+ {"exp", make_call<arb::iexpr>(arb::iexpr::exp, "iexpr with 1 argument: (value:iexpr)")},
+ {"exp", make_call<double>(arb::iexpr::exp, "iexpr with 1 argument: (value:double)")},
+
+ {"log", make_call<arb::iexpr>(arb::iexpr::log, "iexpr with 1 argument: (value:iexpr)")},
+ {"log", make_call<double>(arb::iexpr::log, "iexpr with 1 argument: (value:double)")},
+
+ {"add", make_conversion_fold<arb::iexpr, arb::iexpr, double>(arb::iexpr::add, "iexpr with at least 2 arguments: ((iexpr | double) (iexpr | double) [...(iexpr | double)])")},
+
+ {"sub", make_conversion_fold<arb::iexpr, arb::iexpr, double>(arb::iexpr::sub, "iexpr with at least 2 arguments: ((iexpr | double) (iexpr | double) [...(iexpr | double)])")},
+
+ {"mul", make_conversion_fold<arb::iexpr, arb::iexpr, double>(arb::iexpr::mul, "iexpr with at least 2 arguments: ((iexpr | double) (iexpr | double) [...(iexpr | double)])")},
+
+ {"div", make_conversion_fold<arb::iexpr, arb::iexpr, double>(arb::iexpr::div, "iexpr with at least 2 arguments: ((iexpr | double) (iexpr | double) [...(iexpr | double)])")},
};
parse_label_hopefully<std::any> eval(const s_expr& e);
@@ -254,4 +315,18 @@ ARB_ARBORIO_API parse_label_hopefully<arb::locset> parse_locset_expression(const
}
}
+parse_label_hopefully<arb::iexpr> parse_iexpr_expression(const std::string& s) {
+ if (auto e = eval(parse_s_expr(s))) {
+ if (e->type() == typeid(iexpr)) {
+ return {std::move(std::any_cast<iexpr&>(*e))};
+ }
+ return util::unexpected(
+ label_parse_error(
+ concat("Invalid iexpr description: '", s)));
+ }
+ else {
+ return util::unexpected(label_parse_error(std::string()+e.error().what()));
+ }
+}
+
} // namespace arborio
diff --git a/arborio/parse_helpers.hpp b/arborio/parse_helpers.hpp
index 7a611ae1..5b0c37f9 100644
--- a/arborio/parse_helpers.hpp
+++ b/arborio/parse_helpers.hpp
@@ -4,12 +4,14 @@
#include <string>
#include <sstream>
#include <iostream>
+#include <typeinfo>
#include <arbor/assert.hpp>
#include <arbor/arbexcept.hpp>
#include <arbor/util/expected.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
+#include <arbor/iexpr.hpp>
namespace arborio {
using namespace arb;
@@ -23,18 +25,33 @@ template <> inline
bool match<arb::locset>(const std::type_info& info) { return info == typeid(arb::locset); }
template <> inline
bool match<arb::region>(const std::type_info& info) { return info == typeid(arb::region); }
+template <> inline
+bool match<arb::iexpr>(const std::type_info& info) { return info == typeid(arb::iexpr); }
// Convert a value wrapped in a std::any to target type.
template <typename T>
T eval_cast(std::any arg) {
return std::move(std::any_cast<T&>(arg));
}
+
template <> inline
double eval_cast<double>(std::any arg) {
if (arg.type()==typeid(int)) return std::any_cast<int>(arg);
return std::any_cast<double>(arg);
}
+template <typename T>
+T conversion_cast(std::any arg) {
+ return eval_cast<T>(std::move(arg));
+}
+
+template <typename T, typename Q, typename... Types>
+T conversion_cast(std::any arg) {
+ if (match<Q>(arg.type())) return T(eval_cast<Q>(arg));
+
+ return conversion_cast<T, Types...>(arg);
+}
+
// Test whether a list of arguments passed as a std::vector<std::any> can be converted
// to the types in Args.
//
@@ -99,6 +116,47 @@ struct fold_match {
}
};
+
+// Fold with first converting from any of the "Types" to type T.
+template <typename T, typename... Types>
+struct fold_conversion_eval {
+ using fold_fn = std::function<T(T, T)>;
+ fold_fn f;
+
+ using anyvec = std::vector<std::any>;
+ using iterator = anyvec::iterator;
+
+ fold_conversion_eval(fold_fn f): f(std::move(f)) {}
+
+ T fold_impl(iterator left, iterator right) {
+ if (std::distance(left,right)==1u) {
+ return conversion_cast<T, Types...>(std::move(*left));
+ }
+ // Compute fold. Order is important for left-associative operations like division and
+ // subtraction
+ auto back = right - 1;
+ return f(fold_impl(left, back), conversion_cast<T, Types...>(std::move(*back)));
+ }
+
+ std::any operator()(anyvec args) {
+ return fold_impl(args.begin(), args.end());
+ }
+};
+
+// Test if all args match at least one of the "Types" types
+template <typename... Types>
+struct fold_conversion_match {
+ bool operator()(const std::vector<std::any>& args) const {
+ if (args.size() < 2) return false;
+
+ bool m = true;
+ for (const auto& a : args) {
+ m = m && (match<Types>(a.type()) || ...);
+ }
+ return m;
+ }
+};
+
// Evaluator: member of make_call, make_arg_vec_call, make_mech_call, make_branch_call, make_unordered_call
struct evaluator {
using any_vec = std::vector<std::any>;
@@ -163,6 +221,22 @@ struct make_fold {
}
};
+// Fold with first converting from any of the "Types" to type T.
+template <typename T, typename... Types>
+struct make_conversion_fold {
+ evaluator state;
+
+ template <typename F>
+ make_conversion_fold(F&& f, const char* msg = "fold_conversion"):
+ state(fold_conversion_eval<T, Types...>(std::forward<F>(f)),
+ fold_conversion_match<T, Types...>(),
+ msg) {}
+
+ operator evaluator() const {
+ return state;
+ }
+};
+
// Test whether a list of arguments passed as a std::vector<std::any> can be converted
// to a std::vector<std::variant<Args...>>.
//
diff --git a/doc/concepts/decor.rst b/doc/concepts/decor.rst
index f817193a..03a255fc 100644
--- a/doc/concepts/decor.rst
+++ b/doc/concepts/decor.rst
@@ -62,6 +62,7 @@ The types of dynamics, and where they can be defined, are
cable properties, ✓, ✓, ✓
ion initial conditions, ✓, ✓, ✓
density mechanism, ✓, --, --
+ scaled-mechanism (density), ✓, --, --
ion rev pot mechanism, --, ✓, ✓
ion valence, --, --, ✓
@@ -166,7 +167,28 @@ Take for example the built-in mechanism for passive leaky dynamics:
.. _cablecell-ions:
-4. Ion species
+.. _cablecell-scaled-mechs:
+
+4. Scaled mechanisms
+~~~~~~~~~~~~~~~~~~~~~
+Mechanism parameters are usually homogeneous along a cell. However, sometimes it is useful to scale parameters based on inhomogeneous properties.
+:ref:`Inhomogeneous expressions <labels-iexpr>` provide a way to describe a desired scaling formula, which for example can include the cell radius or the distance to a given set of locations.
+The name is inspired by NeuroML's https://docs.neuroml.org/Userdocs/Schemas/Cells.html#schema-inhomogeneousparameter.
+Such an expression is evaluated along the cell and yields a scaling factor, which is multiplied with the base value of the selected parameter.
+Internally, this evaluation and scaling is done at mid-points of the cable partition of the cell.
+Currently, only parameters of :ref:`density mechanisms <cablecell-density-mechs>` can be scaled.
+
+
+.. code-block:: Python
+
+ # Create mechanism with custom conductance (range)
+ m = arbor.mechanism('pas', {'g': 0.1})
+
+ decor = arbor.decor()
+ # paint a scaled density mechanism, where 'g' is scaled with the distance from the root.
+ decor.paint('"dend"', arbor.scaled_mechanism(arbor.density(m), {'g': '(distance 1.0 (root))'}))
+
+5. Ion species
~~~~~~~~~~~~~~
Arbor allows arbitrary ion species to be defined, to extend the default
diff --git a/doc/concepts/labels.rst b/doc/concepts/labels.rst
index 845674e4..a103e0e5 100644
--- a/doc/concepts/labels.rst
+++ b/doc/concepts/labels.rst
@@ -98,6 +98,10 @@ the region of cables that have radius less than 0.5 μm
.. _labels-expressions:
+.. glossary::
+ iexpr
+ An iexpr is an inhomogeneous expression, that can be evaluated at any point on a cell.
+
Expressions
-----------
@@ -608,6 +612,162 @@ Region expressions
Two regions (left and middle) and their intersection (right).
+
+.. _labels-iexpr:
+
+Inhomogeneous Expressions
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. label:: (scalar value:real)
+
+ A scalar of given value.
+
+.. label:: (pi)
+
+ A scalar expression representing the pi constant.
+
+.. label:: (distance scale:real loc:locset)
+
+ The minimum distance to points within the locset ``loc``. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+ .. figure:: ../images/iexpr_distance.svg
+ :width: 600
+ :align: center
+
+ The distance between any two points (the evaluation location and a location within the locset), is calculated **along** the entire tree, even across the root.
+ Therefore, a distance expression is defined on the entire cell and only zero if evaluated at a location within the locset (or the scale parameter is set to zero).
+
+.. label:: (distance loc:locset)
+
+ A distance expression with a default scaling factor of 1.0.
+
+.. label:: (distance scale:real reg:region)
+
+ The minimum distance to the region ``reg``. Evaluates to zero within the region. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+.. label:: (distance reg:region)
+
+ A distance expression with a default scaling factor of 1.0.
+
+.. label:: (proximal-distance scale:real loc:locset)
+
+ The minimum distance in proximal direction from the points within the locset ``loc``. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+ .. figure:: ../gen-images/iexpr_prox_dis.svg
+ :width: 600
+ :align: center
+
+ Example of a proximal-distance expression with a single input location. **Left**: Input location. **Right**: Area where the expression evaluates to non-zero values.
+
+.. label:: (proximal-distance loc:locset)
+
+ A proximal-distance expression with a default scaling factor of 1.0.
+
+.. label:: (proximal-distance scale:real reg:region)
+
+ The minimum distance in proximal direction from the region ``reg``. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+.. label:: (proximal-distance reg:region)
+
+ A proximal-distance expression with a default scaling factor of 1.0.
+
+.. label:: (distal-distance scale:real loc:locset)
+
+ The minimum distance in distal direction from the points within the locset ``loc``. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+ .. figure:: ../gen-images/iexpr_dist_dis.svg
+ :width: 600
+ :align: center
+
+ Example of a distal-distance expression with a single input location. **Left**: Input location. **Right**: Area, where the expression evaluates to non-zero values.
+
+.. label:: (distal-distance loc:locset)
+
+ A distal-distance expression with a default scaling factor of 1.0.
+
+.. label:: (distal-distance scale:real reg:region)
+
+ The minimum distance in distal direction from the region ``reg``. The scaling parameter ``scale`` has unit :math:`{\mu m}^{-1}`
+ and is multiplied with the distance, such that the result is unitless.
+
+.. label:: (distal-distance reg:region)
+
+ A distal-distance expression with a default scaling factor of 1.0.
+
+.. label:: (interpolation prox_value:real prox_loc:locset dist_value:real dist_loc:locset)
+
+ Interpolates between the closest point in proximal direction in locset ``prox_loc`` and the closest point in
+ distal direction ``dist_loc`` with the assosiated unitless values ``prox_value`` and ``dist_value``.
+ Evaluates to zero, if no point is located in each required direction.
+
+ **Note**: At any fork, an interpolation expression may be discontinuous, if the distance to the closest location within the distal locset differs along each attached branch.
+
+ .. figure:: ../images/iexpr_interp.svg
+ :width: 600
+ :align: center
+
+ Example of an interpolation expression. **Red**: The root of the morphology. **Blue**: The proximal locset consisting of a single location.
+ **Green**: The distal locset, consisting of four locations. Given these locsets, an interpolation expression only evaluates to non-zero in the highlighted area.
+ For locations 3 and 4 of the distal locset, there is no location within the proximal locset, that is between them and the root (in proximal direction),
+ and therefore an interpolation expression cannot be evaluated and defaults to zero.
+ Contrary, for locations 1 and 2 of the distal locset, there is a location within the proximal locset in proximal direction.
+
+
+.. label:: (interpolation prox_value:real prox_reg:region dist_value:real dist_reg:region)
+
+ Interpolates between the region ``prox_reg`` in proximal diretion and the region ``dist_reg`` in distal direction
+ with the associated unitless values ``prox_value`` and ``dist_value``. If evaluated inside either region, returns the corresponding value.
+ Evaluates to zero, if no region is located in each required direction.
+
+.. label:: (radius scale:real)
+
+ The radius of the cell at a given point multiplied with the ``scale`` parameter with unit :math:`{\mu m}^{-1}`.
+
+.. label:: (radius)
+
+ A radius expression with a default scaling factor of 1.0.
+
+.. label:: (diameter scale:real)
+
+ The diameter of the cell at a given point multiplied with the ``scale`` parameter with unit :math:`{\mu m}^{-1}`.
+
+.. label:: (diameter)
+
+ A diameter expression with a default scaling factor of 1.0.
+
+.. label:: (add (iexpr | real) (iexpr | real) [... (iexpr | real)])
+
+ Addition of at least two inhomogeneous expressions or real numbers.
+
+.. label:: (sub (iexpr | real) (iexpr | real) [... (iexpr | real)])
+
+ Subtraction of at least two inhomogeneous expressions or real numbers.
+ The expression is evaluated from the left to right, subtracting each element from the first one in turn.
+
+.. label:: (mul (iexpr | real) (iexpr | real) [... (iexpr | real)])
+
+ Multiplication of at least two inhomogeneous expressions or real numbers.
+
+.. label:: (div (iexpr | real) (iexpr | real) [... (iexpr | real)])
+
+ Division of at least two inhomogeneous expressions or real numbers.
+ The expression is evaluated from the left to right, dividing the first element by each divisor in turn.
+
+.. label:: (exp value:(iexpr | real))
+
+ The exponential function of the inhomogeneous expression or real ``value``.
+
+.. label:: (log value:(iexpr | real))
+
+ The logarithm of the inhomogeneous expression or real ``value``.
+
+
+
.. _labels-thingify:
Thingification
@@ -674,7 +834,7 @@ Label Dictionaries
.. glossary::
label
A label is a string assigned to an :ref:`expression <labels-expressions>`, and used to refer to the expression or the
- concrete :term:`region` or :term:`locset` generated when the expression is applied to a morphology.
+ concrete :term:`region` or :term:`locset` or :term:`iexpr` generated when the expression is applied to a morphology.
Although any string is a valid label, it is a good idea to avoid labels that would
also be valid expressions in the region DSL; creating a label ``"(tag 1)"`` will only
@@ -686,7 +846,7 @@ lead to confusion.
Label dictionaries are used to create a cable-cell along with the :ref:`morphology <morph>`
and a :ref:`decor <cablecell-decoration>`. The decorations can be painted or placed on
-the regions or locsets defined in the label dictionary by referring to their labels.
+the regions, locsets or iexpr defined in the label dictionary by referring to their labels.
.. code-block:: python
:caption: Example of a label dictionary in python:
@@ -697,7 +857,8 @@ the regions or locsets defined in the label dictionary by referring to their lab
'dend': '(tag 3)', # dend is every cable with tab 3 in the morphology
'root': '(root)', # typically the start of the soma is at the root of the cell.
'stim_site': '(location 0 0.5)', # site for the stimulus, in the middle of branch 0.
- 'axon_end': '(restrict (terminal) (region "axon"))'} # end of the axon.
+ 'axon_end': '(restrict (terminal) (region "axon"))', # end of the axon.
+ 'rad_expr': '(radius 0.5)'} # iexpr evaluating the radius scaled by 0.5
})
diff --git a/doc/fileformat/cable_cell.rst b/doc/fileformat/cable_cell.rst
index 745f35b6..70385202 100644
--- a/doc/fileformat/cable_cell.rst
+++ b/doc/fileformat/cable_cell.rst
@@ -54,11 +54,23 @@ The components of the label dictionary are the following:
This expression identifies the midpoint of the branch with id 3 as "my_locset".
-Any number of locset and region definitions can be grouped in a label dictionary as follows:
+.. label:: (iexpr-def label:string e:iexpr)
-.. label:: (label-dict [...def:region-def/locset-def])
+ This defines a ``label`` which can be used to refer to the iexpr ``e``.
+ For example:
+
+ .. code:: lisp
+
+ (iexpr-def "my_iexpr" (radius 0.5))
+
+ This expression identifies the radius iexpr with a scaling factor 0.5.
- This describes a label dictionary of zero or more region and locset definitons.
+
+Any number of locset, region an iexpr definitions can be grouped in a label dictionary as follows:
+
+.. label:: (label-dict [...def:region-def/locset-def/iexpr-def])
+
+ This describes a label dictionary of zero or more region, locset and iexpr definitons.
For example:
.. code:: lisp
@@ -68,7 +80,8 @@ Any number of locset and region definitions can be grouped in a label dictionary
(locset-def "root" (root))
(region-def "all" (all))
(region-def "my_region" (radius-ge (region "my_soma") 1.5))
- (locset-def "terminal" (terminal)))
+ (locset-def "terminal" (terminal))
+ (iexpr-def "my_iexpr" (radius 0.5)))
Decor
-----
@@ -84,20 +97,21 @@ of the cell):
.. csv-table:: Property applicability.
:widths: 20, 10, 10, 10
- , **placeable**, **paintable**, **defaultable**
- initial membrane potential, --, ✓, ✓
- axial resistivity, --, ✓, ✓
- temperature, --, ✓, ✓
- membrane capacitance, --, ✓, ✓
- ion initial internal concentration, --, ✓, ✓
- ion initial external concentration, --, ✓, ✓
- ion initial reversal potential, --, ✓, ✓
- ion reversal potential method, --, --, ✓
- density mechanism, --, ✓, --
- point mechanism, ✓, --, --
- junction mechanism, ✓, --, --
- current clamp, ✓, --, --
- threshold detector, ✓, --, --
+ , **placeable**, **paintable**, **defaultable**
+ initial membrane potential, --, ✓, ✓
+ axial resistivity, --, ✓, ✓
+ temperature, --, ✓, ✓
+ membrane capacitance, --, ✓, ✓
+ ion initial internal concentration, --, ✓, ✓
+ ion initial external concentration, --, ✓, ✓
+ ion initial reversal potential, --, ✓, ✓
+ ion reversal potential method, --, --, ✓
+ density mechanism, --, ✓, --
+ scaled-mechanism (density), --, ✓, --
+ point mechanism, ✓, --, --
+ junction mechanism, ✓, --, --
+ current clamp, ✓, --, --
+ threshold detector, ✓, --, --
The various properties and dynamics of the decor are described as follows:
@@ -155,15 +169,20 @@ The various properties and dynamics of the decor are described as follows:
.. label:: (density method:mechanism)
- This describes a *density* mechanism whose behavior is is defined by ``mechanism``.
+ This describes a *density* mechanism whose behavior is defined by ``mechanism``.
+
+.. label:: (scaled-mechanism p:density [...(param:string e:iexpr)])
+
+ This describes a *density* mechanism, which is modified by scaling of individual parameters with
+ inhomogeneous scaling expressions.
.. label:: (synapse method:mechanism)
- This describes a *synapse* (point) mechanism whose behavior is is defined by ``mechanism``.
+ This describes a *synapse* (point) mechanism whose behavior is defined by ``mechanism``.
.. label:: (junction method:mechanism)
- This describes a *gap-junction* mechanism whose behavior is is defined by ``mechanism``.
+ This describes a *gap-junction* mechanism whose behavior is defined by ``mechanism``.
.. label:: (current-clamp (envelope-pulse delay:real duration:real amplitude:real) freq:real phase:real)
diff --git a/doc/images/iexpr_distance.svg b/doc/images/iexpr_distance.svg
new file mode 100644
index 00000000..0782b49c
--- /dev/null
+++ b/doc/images/iexpr_distance.svg
@@ -0,0 +1,45 @@
+<?xml version="1.0" encoding="utf-8"?>
+<!-- Generator: Adobe Illustrator 25.2.1, SVG Export Plug-In . SVG Version: 6.00 Build 0) -->
+<svg version="1.1" id="Layer_1" xmlns:ev="http://www.w3.org/2001/xml-events"
+ xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px" viewBox="0 0 1920 1080"
+ style="enable-background:new 0 0 1920 1080;" xml:space="preserve">
+<style type="text/css">
+ .st0{fill:#D3D3D3;stroke:#D3D3D3;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st1{fill:#FFFFFF;stroke:#FF0000;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st2{fill:#FFFFFF;stroke:#000000;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st3{fill:none;stroke:#000000;stroke-width:4;stroke-miterlimit:10;}
+</style>
+<g id="lines">
+ <polygon class="st0" points="577,612.7 770.9,612.7 770.9,554.5 964.8,554.5 960,554.2 1154,578.5 1163.6,501.5 969.7,477.3
+ 964.8,477 770.9,477 770.9,418.8 577,418.8 "/>
+ <polygon class="st0" points="1177.8,573.8 1556.2,338.7 1552.8,340.2 1840.6,234 1834.5,215.6 1540.5,303.4 1537.2,304.9
+ 1139.8,506.2 "/>
+ <polygon class="st0" points="1150.6,562.8 1490,684 1506.3,638.4 1166.9,517.2 "/>
+ <polygon class="st0" points="1483,680.1 1734.5,862.7 1746.7,847.6 1513.3,642.3 "/>
+ <polygon class="st0" points="1509,682.9 1696.5,572.9 1689.1,573.4 1834.5,621.9 1840.6,603.5 1695.2,555 1687.8,555.6
+ 1487.3,639.5 "/>
+ <polygon class="st0" points="577,418.8 237.6,496.4 237.6,496.4 92.1,496.4 92.1,535.2 237.6,535.2 237.6,535.2 577,612.7 "/>
+</g>
+<g id="points">
+ <circle class="st1" cx="577" cy="515.8" r="38.8"/>
+ <circle class="st2" cx="1486.4" cy="355.7" r="32.3"/>
+ <circle class="st2" cx="1663.7" cy="578.4" r="32.3"/>
+</g>
+<g>
+ <g>
+ <line class="st3" x1="1420.2" y1="392.4" x2="1171.9" y2="532.5"/>
+ <g>
+ <polygon points="1444.8,378.5 1402.2,383.9 1418.4,393.4 1418.2,412.2 "/>
+ </g>
+ </g>
+</g>
+<line class="st3" x1="1506.3" y1="661.5" x2="1171.9" y2="532.5"/>
+<g>
+ <g>
+ <line class="st3" x1="1597.3" y1="613.7" x2="1506.3" y2="661.5"/>
+ <g>
+ <polygon points="1622.3,600.6 1579.5,604.7 1595.4,614.7 1594.6,633.5 "/>
+ </g>
+ </g>
+</g>
+</svg>
diff --git a/doc/images/iexpr_interp.svg b/doc/images/iexpr_interp.svg
new file mode 100644
index 00000000..88c77836
--- /dev/null
+++ b/doc/images/iexpr_interp.svg
@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="utf-8"?>
+<!-- Generator: Adobe Illustrator 25.2.1, SVG Export Plug-In . SVG Version: 6.00 Build 0) -->
+<svg version="1.1" id="Layer_1" xmlns:ev="http://www.w3.org/2001/xml-events"
+ xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px" viewBox="0 0 1920 1080"
+ style="enable-background:new 0 0 1920 1080;" xml:space="preserve">
+<style type="text/css">
+ .st0{fill:#D3D3D3;stroke:#D3D3D3;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st1{stroke:#D3D3D3;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st2{fill:#FFFFFF;stroke:#FF0000;stroke-width:9.697;stroke-miterlimit:9.697;}
+ .st3{fill:none;stroke:#00BF0A;stroke-width:12;stroke-miterlimit:10;}
+ .st4{fill:none;stroke:#006FE9;stroke-width:12;stroke-miterlimit:10;}
+ .st5{fill:none;}
+ .st6{font-family:'MyriadPro-Regular';}
+ .st7{font-size:60px;}
+</style>
+<g id="lines">
+ <polygon class="st0" points="577,612.7 770.9,612.7 770.9,554.5 964.8,554.5 960,554.2 1154,578.5 1163.6,501.5 969.7,477.3
+ 964.8,477 770.9,477 770.9,418.8 577,418.8 "/>
+ <polygon class="st0" points="1177.8,573.8 1556.2,338.7 1552.8,340.2 1840.6,234 1834.5,215.6 1540.5,303.4 1537.2,304.9
+ 1139.8,506.2 "/>
+ <polygon class="st0" points="1150.6,562.8 1490,684 1506.3,638.4 1166.9,517.2 "/>
+ <polygon class="st0" points="1483,680.1 1734.5,862.7 1746.7,847.6 1513.3,642.3 "/>
+ <polygon class="st0" points="1509,682.9 1696.5,572.9 1689.1,573.4 1834.5,621.9 1840.6,603.5 1695.2,555 1687.8,555.6
+ 1487.3,639.5 "/>
+ <polygon class="st0" points="577,418.8 237.6,496.4 237.6,496.4 92.1,496.4 92.1,535.2 237.6,535.2 237.6,535.2 577,612.7 "/>
+ <polygon class="st1" points="1320.3,623.4 1490,684 1506.3,638.4 1336.6,577.8 "/>
+ <polygon class="st1" points="1483,680.1 1734.5,862.7 1746.7,847.6 1513.3,642.3 "/>
+ <polygon class="st1" points="1509,682.9 1696.5,572.9 1689.1,573.4 1834.5,621.9 1840.6,603.5 1695.2,555 1687.8,555.6
+ 1487.3,639.5 "/>
+</g>
+<g id="points">
+ <circle class="st2" cx="577" cy="515.8" r="38.8"/>
+</g>
+<circle class="st3" cx="1840.6" cy="612.7" r="24.8"/>
+<circle class="st4" cx="1336.2" cy="600.6" r="24.8"/>
+<circle class="st3" cx="92.1" cy="517.2" r="24.8"/>
+<circle class="st3" cx="1739.7" cy="855.1" r="24.8"/>
+<circle class="st3" cx="1840.6" cy="224.1" r="24.8"/>
+<rect x="1692.5" y="884.4" class="st5" width="72" height="82"/>
+<text transform="matrix(1 0 0 1 1692.4697 927.0187)" class="st6 st7">1</text>
+<rect x="1793.4" y="643" class="st5" width="72" height="82"/>
+<text transform="matrix(1 0 0 1 1793.4165 685.6184)" class="st6 st7">2</text>
+<rect x="1804.6" y="267.2" class="st5" width="72" height="82"/>
+<text transform="matrix(1 0 0 1 1804.6423 309.8443)" class="st6 st7">3</text>
+<rect x="44.9" y="559.6" class="st5" width="72" height="82"/>
+<text transform="matrix(1 0 0 1 44.8954 602.1823)" class="st6 st7">4</text>
+</svg>
diff --git a/doc/scripts/inputs.py b/doc/scripts/inputs.py
index a59dd0ae..de2ea3cd 100644
--- a/doc/scripts/inputs.py
+++ b/doc/scripts/inputs.py
@@ -521,6 +521,16 @@ reg_rhs = {"type": "region", "value": [(1, 0.0, 1.0)]}
reg_and = {"type": "region", "value": [(1, 0.0, 0.5)]}
reg_or = {"type": "region", "value": [(0, 0.5, 1.0), (1, 0.0, 1.0)]}
+
+############# iexpr (label_morph)
+
+iexpr_directional_loc = {"type": "locset", "value": [(0, 1.0)]}
+iexpr_dist_dis = {
+ "type": "region",
+ "value": [(1, 0.0, 1.0), (2, 0.0, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0)],
+}
+iexpr_prox_dis = {"type": "region", "value": [(0, 0.0, 1.0)]}
+
############# locsets (tutorial_morph)
tut_ls_root = {"type": "locset", "value": [(0, 0.0)]}
diff --git a/doc/scripts/make_images.py b/doc/scripts/make_images.py
index 0c742a37..2e6ee0c9 100644
--- a/doc/scripts/make_images.py
+++ b/doc/scripts/make_images.py
@@ -511,6 +511,19 @@ def generate(path=""):
label_image(inputs.label_morph, [inputs.reg_radgt5], path + "/radiusgt_label.svg")
label_image(inputs.label_morph, [inputs.reg_radge5], path + "/radiusge_label.svg")
+ ####################### iexpr
+
+ label_image(
+ inputs.label_morph,
+ [inputs.iexpr_directional_loc, inputs.iexpr_prox_dis],
+ path + "/iexpr_prox_dis.svg",
+ )
+ label_image(
+ inputs.label_morph,
+ [inputs.iexpr_directional_loc, inputs.iexpr_dist_dis],
+ path + "/iexpr_dist_dis.svg",
+ )
+
####################### Tutorial examples
morph_image([inputs.tutorial_morph], ["segments"], path + "/tutorial_morph.svg")
diff --git a/python/cells.cpp b/python/cells.cpp
index 73384031..1020143d 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -131,6 +131,11 @@ std::string mechanism_desc_str(const arb::mechanism_desc& md) {
md.name(), util::dictionary_csv(md.values()));
}
+std::string scaled_density_desc_str(const arb::scaled_mechanism<arb::density>& p) {
+ return util::pprintf("({}, {})",
+ mechanism_desc_str(p.t_mech.mech), util::dictionary_csv(p.scale_expr));
+}
+
void register_cells(pybind11::module& m) {
using namespace pybind11::literals;
using std::optional;
@@ -424,6 +429,38 @@ 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) + ">";});
+ // arb::scaled_mechanism<arb::density>
+
+ pybind11::class_<arb::scaled_mechanism<arb::density>> scaled_mechanism(
+ m, "scaled_mechanism", "For painting a scaled density mechanism on a region.");
+ scaled_mechanism
+ .def(pybind11::init(
+ [](arb::density dens) { return arb::scaled_mechanism<arb::density>(std::move(dens)); }))
+ .def(pybind11::init(
+ [](arb::density dens, const std::unordered_map<std::string, std::string>& scales) {
+ auto s = arb::scaled_mechanism<arb::density>(std::move(dens));
+ for (const auto& it: scales) {
+ s.scale(it.first, arborio::parse_iexpr_expression(it.second).unwrap());
+ }
+ return s;
+ }))
+ .def(
+ "scale",
+ [](arb::scaled_mechanism<arb::density>& s, std::string name, const std::string& ex) {
+ s.scale(std::move(name), arborio::parse_iexpr_expression(ex).unwrap());
+ return s;
+ },
+ pybind11::arg_v("name", "name of parameter to scale."),
+ pybind11::arg_v("ex", "iexpr for scaling"),
+ "Add a scaling expression to a parameter.")
+ .def("__repr__",
+ [](const arb::scaled_mechanism<arb::density>& d) {
+ return "<arbor.scaled_mechanism<density> " + scaled_density_desc_str(d) + ">";
+ })
+ .def("__str__", [](const arb::scaled_mechanism<arb::density>& d) {
+ return "<arbor.scaled_mechanism<density> " + scaled_density_desc_str(d) + ">";
+ });
+
// arb::synapse
pybind11::class_<arb::synapse> synapse(m, "synapse", "For placing a synaptic mechanism on a locset.");
@@ -659,6 +696,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::scaled_mechanism<arb::density>& mechanism) {
+ dec.paint(arborio::parse_region_expression(region).unwrap(), mechanism);
+ },
+ "region"_a, "mechanism"_a,
+ "Associate a scaled density mechanism with a region.")
// Paint membrane/static properties.
.def("paint",
[](arb::decor& dec,
diff --git a/python/proxy.hpp b/python/proxy.hpp
index 4d79f7c8..31333e13 100644
--- a/python/proxy.hpp
+++ b/python/proxy.hpp
@@ -1,6 +1,7 @@
#pragma once
#include <any>
+#include <string>
#include <arborio/label_parse.hpp>
@@ -9,12 +10,19 @@
#include "strprintf.hpp"
namespace pyarb {
+
+struct iexpr_proxy {
+ ::arb::iexpr iexpr;
+ std::string label;
+};
+
struct label_dict_proxy {
using str_map = std::unordered_map<std::string, std::string>;
arb::label_dict dict;
str_map cache;
std::vector<std::string> locsets;
std::vector<std::string> regions;
+ std::vector<std::string> iexpressions;
label_dict_proxy() = default;
@@ -39,7 +47,7 @@ struct label_dict_proxy {
label_dict_proxy(const label_dict_proxy&) = default;
std::size_t size() const {
- return locsets.size() + regions.size();
+ return locsets.size() + regions.size() + iexpressions.size();
}
void import(const label_dict_proxy& other, std::string prefix = "") {
@@ -51,20 +59,20 @@ struct label_dict_proxy {
void set(const std::string& name, const std::string& desc) {
using namespace std::string_literals;
- // The following code takes an input name and a region or locset
+ // The following code takes an input name and a region or locset or iexpr
// description, e.g.:
// name='reg', desc='(tag 4)'
// name='loc', desc='(terminal)'
// name='foo', desc='(join (tag 2) (tag 3))'
// Then it parses the description, and tests whether the description
- // is a region or locset, and updates the label dictionary appropriately.
+ // is a region or locset or iexpr, and updates the label dictionary appropriately.
// Errors occur when:
- // * a region is described with a name that matches an existing locset
+ // * a region is described with a name that matches an existing locset or iexpr
// (and vice versa.)
// * the description is not well formed, e.g. it contains a syntax error.
- // * the description is well-formed, but describes neither a region or locset.
+ // * the description is well-formed, but describes neither a region or locset or iexpr.
try {
- // Evaluate the s-expression to build a region/locset.
+ // Evaluate the s-expression to build a region/locset/iexpr.
auto result = arborio::parse_label_expression(desc);
if (!result) { // an error parsing / evaluating description.
throw result.error();
@@ -79,6 +87,11 @@ struct label_dict_proxy {
auto it = std::lower_bound(locsets.begin(), locsets.end(), name);
if (it==locsets.end() || *it!=name) locsets.insert(it, name);
}
+ else if (result->type()==typeid(arb::iexpr)) { // describes a iexpr.
+ dict.set(name, std::move(std::any_cast<arb::iexpr&>(*result)));
+ auto it = std::lower_bound(iexpressions.begin(), iexpressions.end(), name);
+ if (it==iexpressions.end() || *it!=name) iexpressions.insert(it, name);
+ }
else {
// Successfully parsed an expression that is neither region nor locset.
throw util::pprintf("The definition of '{} = {}' does not define a valid region or locset.", name, desc);
@@ -115,10 +128,13 @@ struct label_dict_proxy {
for (auto& x: dict.locsets()) {
s += util::pprintf(" (locset \"{}\" {})", x.first, x.second);
}
+ for (auto& x: dict.iexpressions()) {
+ s += util::pprintf(" (iexpr \"{}\" {})", x.first, x.second);
+ }
s += ")";
return s;
}
-
+
bool contains(const std::string& name) const {
return cache.find(name) != cache.end();
}
@@ -135,6 +151,7 @@ struct label_dict_proxy {
void clear_cache() {
regions.clear();
locsets.clear();
+ iexpressions.clear();
cache.clear();
}
@@ -155,9 +172,18 @@ struct label_dict_proxy {
cache[lab] = s.str();
}
}
+ for (const auto& [lab, ls]: dict.iexpressions()) {
+ if (!cache.count(lab)) {
+ std::stringstream s;
+ s << ls;
+ iexpressions.push_back(lab);
+ cache[lab] = s.str();
+ }
+ }
// Sort the region and locset names
std::sort(regions.begin(), regions.end());
std::sort(locsets.begin(), locsets.end());
+ std::sort(iexpressions.begin(), iexpressions.end());
}
};
}
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 2fe2ff90..db708766 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -113,6 +113,7 @@ set(unit_sources
test_fvm_layout.cpp
test_fvm_lowered.cpp
test_diffusion.cpp
+ test_iexpr.cpp
test_index.cpp
test_kinetic_linear.cpp
test_lexcmp.cpp
diff --git a/test/unit/test_iexpr.cpp b/test/unit/test_iexpr.cpp
new file mode 100644
index 00000000..9043381a
--- /dev/null
+++ b/test/unit/test_iexpr.cpp
@@ -0,0 +1,465 @@
+#include "../gtest.h"
+#include "../common_cells.hpp"
+#include "fvm_layout.hpp"
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/iexpr.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arborio/label_parse.hpp>
+
+#include <cmath>
+
+using namespace arb;
+using namespace arborio::literals;
+
+TEST(iexpr, distance_locset) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 2.0;
+
+ // test distance to single point
+ auto ex = thingify(arb::iexpr::distance(scale, arb::mlocation{1, 0.2}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), scale * 3.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 12.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 22.0);
+
+ // test distance to multiple points
+ ex = thingify(arb::iexpr::distance(
+ scale, arb::mlocation_list({arb::mlocation{1, 1.0}, arb::mlocation{2, 1.0}})),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 15.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), scale * 5.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 10.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 30.0);
+}
+
+TEST(iexpr, distance_region) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 3.0;
+
+ // test distance to single cable region
+ auto ex = thingify(arb::iexpr::distance(scale, arb::mcable{0, 0.2, 0.8}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 12.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 12.0);
+
+ // test distance to multi cable region
+ ex = thingify(arb::iexpr::distance(
+ scale, arb::mcable_list{arb::mcable{1, 0.2, 0.8}, arb::mcable{2, 0.6, 1.0}}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 2.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 22.0);
+}
+
+TEST(iexpr, proximal_distance_locset) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 2.0;
+
+ // test distance to single point
+ auto ex = thingify(arb::iexpr::proximal_distance(scale, arb::mlocation{1, 0.2}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multiple points
+ ex = thingify(arb::iexpr::proximal_distance(
+ scale, arb::mlocation_list({arb::mlocation{1, 0.2}, arb::mlocation{2, 1.0}})),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 10.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+}
+
+TEST(iexpr, proximal_distance_region) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 3.0;
+
+ // test distance to single cable region
+ auto ex = thingify(arb::iexpr::proximal_distance(scale, arb::mcable{1, 0.2, 0.8}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multi cable region
+ ex = thingify(
+ arb::iexpr::proximal_distance(scale,
+ arb::mcable_list{
+ arb::mcable{1, 0.2, 0.8}, arb::mcable{2, 0.6, 1.0}, arb::mcable{3, 0.1, 0.3}}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 7.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 2.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+}
+
+TEST(iexpr, distal_distance_locset) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 2.0;
+
+ // test distance to single point
+ auto ex = thingify(arb::iexpr::distal_distance(scale, arb::mlocation{1, 0.2}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), scale * 3.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multiple points
+ ex = thingify(arb::iexpr::distal_distance(
+ scale, arb::mlocation_list({arb::mlocation{1, 0.2}, arb::mlocation{3, 0.4}})),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), scale * 3.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 2.0);
+
+ // test distance to root
+ ex = thingify(arb::iexpr::distal_distance(scale, arb::ls::root()), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), scale * 5.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 10.0);
+}
+
+TEST(iexpr, distal_distance_region) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 3.0;
+
+ // test distance to single cable region
+ auto ex = thingify(arb::iexpr::distal_distance(scale, arb::mcable{1, 0.2, 0.8}), prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multi cable region
+ ex = thingify(
+ arb::iexpr::distal_distance(scale,
+ arb::mcable_list{
+ arb::mcable{1, 0.2, 0.8}, arb::mcable{2, 0.3, 0.4}, arb::mcable{3, 0.1, 0.2}}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), scale * 2.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), scale * 6.0);
+}
+
+TEST(iexpr, interpolation_locset) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double prox_value = 2.0;
+ const double dist_value = 3.0;
+
+ // test single point
+ auto ex = thingify(arb::iexpr::interpolation(
+ prox_value, arb::mlocation{1, 0.2}, dist_value, arb::mlocation{1, 0.2}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test evaluation on ends of interval
+ ex = thingify(arb::iexpr::interpolation(
+ prox_value, arb::mlocation{1, 0.2}, dist_value, arb::mlocation{1, 0.8}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.2, 0.2}), prox_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.8, 0.8}), dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multiple points
+ ex =
+ thingify(arb::iexpr::interpolation(prox_value,
+ arb::mlocation{0, 0.3},
+ dist_value,
+ arb::mlocation_list(
+ {arb::mlocation{1, 0.2}, arb::mlocation{2, 0.8}, arb::mlocation{3, 1.0}})),
+ prov);
+ EXPECT_DOUBLE_EQ(
+ ex->eval(prov, {0, 0.0, 1.0}), 7.0 / 9.0 * prox_value + 2.0 / 9.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(
+ ex->eval(prov, {2, 0.0, 1.0}), 6.0 / 23.0 * prox_value + 17.0 / 23.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test root
+ ex = thingify(
+ arb::iexpr::interpolation(prox_value, arb::ls::root(), dist_value, arb::ls::terminal()),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 15.0 / 20.0 * prox_value + 5.0 / 20.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 5.0 / 20.0 * prox_value + 15.0 / 20.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 10.0 / 30.0 * prox_value + 20.0 / 30.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 10.0 / 20.0 * prox_value + 10.0 / 20.0 * dist_value);
+}
+
+TEST(iexpr, interpolation_region) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 1}, 3);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 30, 1}, 4);
+ tree.append(mnpos, {0, 0, 0, 2}, {0, 0, -20, 2}, 2);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double prox_value = 2.0;
+ const double dist_value = 3.0;
+
+ // test distance to single point
+ auto ex =
+ thingify(arb::iexpr::interpolation(
+ prox_value, arb::mcable{1, 0.2, 0.7}, dist_value, arb::mcable{1, 0.2, 0.7}),
+ prov);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {0, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+
+ // test distance to multiple points
+ ex = thingify(
+ arb::iexpr::interpolation(prox_value,
+ arb::mcable{0, 0.1, 0.3},
+ dist_value,
+ arb::mcable_list(
+ {arb::mcable{1, 0.2, 0.4}, arb::mcable{2, 0.2, 1.0}, arb::mcable{3, 0.1, 1.0}})),
+ prov);
+ EXPECT_DOUBLE_EQ(
+ ex->eval(prov, {0, 0.0, 1.0}), 7.0 / 9.0 * prox_value + 2.0 / 9.0 * dist_value);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {1, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {2, 0.0, 1.0}), 0.0);
+ EXPECT_DOUBLE_EQ(ex->eval(prov, {3, 0.0, 1.0}), 0.0);
+}
+
+TEST(iexpr, scalar) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::scalar(2.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0);
+}
+
+TEST(iexpr, radius) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 2}, 2);
+ tree.append(0, {0, 0, 10, 10}, {0, 0, 30, 5}, 3);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 2.0;
+ auto e = thingify(arb::iexpr::radius(scale), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), scale * 10.0);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {1, 0.0, 1.0}), scale * 1.5);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {2, 0.0, 1.0}), scale * 7.5);
+}
+
+TEST(iexpr, diameter) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+ tree.append(0, {0, 0, 10, 1}, {0, 0, 20, 2}, 2);
+ tree.append(0, {0, 0, 10, 10}, {0, 0, 30, 5}, 3);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ const double scale = 2.0;
+ auto e = thingify(arb::iexpr::diameter(scale), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), scale * 20.0);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {1, 0.0, 1.0}), scale * 3.0);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {2, 0.0, 1.0}), scale * 15.0);
+}
+
+TEST(iexpr, add) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::add(arb::iexpr::scalar(2.0), arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 + 3.0 * 10.0);
+
+ // check operator
+ e = thingify(2.0 + arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 + 3.0 * 10.0);
+
+ // check unary operator
+ e = thingify(+arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 3.0 * 10.0);
+}
+
+TEST(iexpr, sub) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::sub(arb::iexpr::scalar(2.0), arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 - 3.0 * 10.0);
+
+ // check operator
+ e = thingify(2.0 - arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 - 3.0 * 10.0);
+
+ // check unary operator
+ e = thingify(-arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), -3.0 * 10.0);
+}
+
+TEST(iexpr, mul) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::mul(arb::iexpr::scalar(2.0), arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 * 3.0 * 10.0);
+
+ // check operator
+ e = thingify(2.0 * arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 * 3.0 * 10.0);
+}
+
+TEST(iexpr, div) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::div(arb::iexpr::scalar(2.0), arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 / (3.0 * 10.0));
+
+ // check operator
+ e = thingify(2.0 / arb::iexpr::radius(3.0), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), 2.0 / (3.0 * 10.0));
+}
+
+TEST(iexpr, exp) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::exp(arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), std::exp(3.0 * 10.0));
+}
+
+TEST(iexpr, log) {
+ segment_tree tree;
+ tree.append(mnpos, {0, 0, 0, 10}, {0, 0, 10, 10}, 1);
+
+ arb::mprovider prov(arb::morphology(std::move(tree)));
+
+ auto e = thingify(arb::iexpr::log(arb::iexpr::radius(3.0)), prov);
+ EXPECT_DOUBLE_EQ(e->eval(prov, {0, 0.0, 1.0}), std::log(3.0 * 10.0));
+}
+
+TEST(iexpr, fvm_layout) {
+ const double radius = 2.0;
+ const double radius_scale = 3.0;
+ const std::string scaled_param = "gnabar";
+ auto iexpr_radius = arb::iexpr::radius(radius_scale);
+
+ soma_cell_builder builder(12.6157 / 2.0);
+ builder.add_branch(0, 100, radius, radius, 4, "dend");
+
+ auto desc_ref = builder.make_cell();
+ desc_ref.decorations.paint("soma"_lab, density("pas"));
+ desc_ref.decorations.paint("dend"_lab, density("hh"));
+
+ auto desc_scaled = builder.make_cell();
+ desc_scaled.decorations.paint("soma"_lab, density("pas"));
+ desc_scaled.decorations.paint("dend"_lab, scaled_mechanism(density("hh")).scale(scaled_param, iexpr_radius));
+
+ cable_cell_global_properties gprop_coalesce;
+ gprop_coalesce.default_parameters = neuron_parameter_defaults;
+ gprop_coalesce.coalesce_synapses = true;
+
+ cable_cell cell_ref(desc_ref);
+ fvm_mechanism_data data_ref = fvm_build_mechanism_data(gprop_coalesce,
+ {cell_ref},
+ {0},
+ {{0, {}}},
+ fvm_cv_discretize({cell_ref}, neuron_parameter_defaults));
+
+ cable_cell cell_scaled(desc_scaled);
+ fvm_mechanism_data data_scaled = fvm_build_mechanism_data(gprop_coalesce,
+ {cell_scaled},
+ {0},
+ {{0, {}}},
+ fvm_cv_discretize({cell_scaled}, neuron_parameter_defaults));
+
+ // compare parameter values between reference and scaled fvm data
+ for(const auto& m_ref : data_ref.mechanisms) {
+ auto it_scaled = data_scaled.mechanisms.find(m_ref.first);
+ ASSERT_NE(it_scaled, data_scaled.mechanisms.end()); // make sure mech exists in both
+
+ const auto& param_ref = m_ref.second.param_values;
+ const auto& param_scaled = it_scaled->second.param_values;
+ ASSERT_EQ(param_ref.size(), param_scaled.size());
+
+ for (auto i: util::count_along(m_ref.second.param_values)) {
+ ASSERT_STREQ(param_ref[i].first.c_str(), param_scaled[i].first.c_str());
+ ASSERT_EQ(param_ref[i].second.size(), param_ref[i].second.size());
+
+ const double scale = param_scaled[i].first == scaled_param ? radius * radius_scale : 1.0;
+ for (auto j: util::count_along(param_ref[i].second)) {
+ EXPECT_DOUBLE_EQ(scale * param_ref[i].second[j], param_scaled[i].second[j]);
+ }
+ }
+ }
+}
diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp
index be436e78..9920ee86 100644
--- a/test/unit/test_s_expr.cpp
+++ b/test/unit/test_s_expr.cpp
@@ -206,6 +206,15 @@ std::string round_trip_locset(const char* in) {
}
}
+std::string round_trip_iexpr(const char* in) {
+ if (auto x = parse_iexpr_expression(in)) {
+ return util::pprintf("{}", std::any_cast<arb::iexpr>(*x));
+ }
+ else {
+ return x.error().what();
+ }
+}
+
TEST(cv_policies, round_tripping) {
auto literals = {"(every-segment (tag 42))",
@@ -244,6 +253,88 @@ TEST(cv_policies, bad) {
EXPECT_THROW(check("(tag 42)"), cv_policy_parse_error); // not a cv_policy
}
+TEST(iexpr, round_tripping) {
+ EXPECT_EQ("(cable 3 0 1)", round_trip_label<arb::region>("(branch 3)"));
+ EXPECT_EQ("(intersect (tag 1) (intersect (tag 2) (tag 3)))", round_trip_label<arb::region>("(intersect (tag 1) (tag 2) (tag 3))"));
+ auto iexpr_literals = {
+ "(scalar 2.1)",
+ "(distance 3.2 (region \"foo\"))",
+ "(distance 3.2 (location 3 0.2))",
+ "(proximal-distance 3.2 (region \"foo\"))",
+ "(proximal-distance 3.2 (location 3 0.2))",
+ "(distal-distance 3.2 (region \"foo\"))",
+ "(distal-distance 3.2 (location 3 0.2))",
+ "(interpolation 3.2 (region \"foo\") 4.3 (radius-gt (tag 3) 1))",
+ "(interpolation 3.2 (location 3 0.2) 4.3 (distal (tag 2)))",
+ "(radius 2.1)",
+ "(diameter 2.1)",
+ "(exp (scalar 2.1))",
+ "(log (scalar 2.1))",
+ "(add (scalar 2.1) (radius 3.2))",
+ "(sub (scalar 2.1) (radius 3.2))",
+ "(mul (scalar 2.1) (radius 3.2))",
+ "(div (scalar 2.1) (radius 3.2))",
+ };
+ for (auto l: iexpr_literals) {
+ EXPECT_EQ(l, round_trip_label<arb::iexpr>(l));
+ EXPECT_EQ(l, round_trip_iexpr(l));
+ }
+
+ // check double values for input instead of explicit scalar iexpr
+ auto mono_iexpr = {std::string("exp"), std::string("log")};
+ auto duo_iexpr = {std::string("add"), std::string("sub"), std::string("mul"), std::string("div")};
+ constexpr auto v1 = "1.2";
+ constexpr auto v2 = "1.2";
+ for(const auto& l : mono_iexpr) {
+ auto l_input = "(" + l + " " + v1 + ")";
+ auto l_output = "(" + l + " (scalar " + v1 + "))";
+ EXPECT_EQ(l_output.c_str(), round_trip_label<arb::iexpr>(l_input.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_iexpr(l_input.c_str()));
+ }
+ for(const auto& l : duo_iexpr) {
+ auto l_input_dd = "(" + l + " " + v1 + + " " + v2 +")";
+ auto l_input_sd = "(" + l + " (scalar " + v1 + + ") " + v2 +")";
+ auto l_input_ds = "(" + l + " " + v1 + + " (scalar " + v2 +"))";
+ auto l_output = "(" + l + " (scalar " + v1 + ") (scalar " + v2 +"))";
+ EXPECT_EQ(l_output.c_str(), round_trip_label<arb::iexpr>(l_input_dd.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_iexpr(l_input_dd.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_label<arb::iexpr>(l_input_sd.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_iexpr(l_input_sd.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_label<arb::iexpr>(l_input_ds.c_str()));
+ EXPECT_EQ(l_output.c_str(), round_trip_iexpr(l_input_ds.c_str()));
+ }
+
+ // test order for more than two arguments
+ EXPECT_EQ("(add (add (add (scalar 1.1) (scalar 2.2)) (scalar 3.3)) (scalar 4.4))",
+ round_trip_label<arb::iexpr>("(add 1.1 2.2 3.3 4.4)"));
+ EXPECT_EQ("(sub (sub (sub (scalar 1.1) (scalar 2.2)) (scalar 3.3)) (scalar 4.4))",
+ round_trip_label<arb::iexpr>("(sub 1.1 2.2 3.3 4.4)"));
+ EXPECT_EQ("(mul (mul (mul (scalar 1.1) (scalar 2.2)) (scalar 3.3)) (scalar 4.4))",
+ round_trip_label<arb::iexpr>("(mul 1.1 2.2 3.3 4.4)"));
+ EXPECT_EQ("(div (div (div (scalar 1.1) (scalar 2.2)) (scalar 3.3)) (scalar 4.4))",
+ round_trip_label<arb::iexpr>("(div 1.1 2.2 3.3 4.4)"));
+
+ // test default constructors
+ EXPECT_EQ("(distance 1 (location 3 0.2))",
+ round_trip_label<arb::iexpr>("(distance (location 3 0.2))"));
+ EXPECT_EQ("(distance 1 (region \"foo\"))",
+ round_trip_label<arb::iexpr>("(distance (region \"foo\"))"));
+ EXPECT_EQ("(distal-distance 1 (location 3 0.2))",
+ round_trip_label<arb::iexpr>("(distal-distance (location 3 0.2))"));
+ EXPECT_EQ("(distal-distance 1 (region \"foo\"))",
+ round_trip_label<arb::iexpr>("(distal-distance (region \"foo\"))"));
+ EXPECT_EQ("(proximal-distance 1 (location 3 0.2))",
+ round_trip_label<arb::iexpr>("(proximal-distance (location 3 0.2))"));
+ EXPECT_EQ("(proximal-distance 1 (region \"foo\"))",
+ round_trip_label<arb::iexpr>("(proximal-distance (region \"foo\"))"));
+ EXPECT_EQ("(radius 1)",
+ round_trip_label<arb::iexpr>("(radius)"));
+ EXPECT_EQ("(diameter 1)",
+ round_trip_label<arb::iexpr>("(diameter)"));
+ EXPECT_EQ("(scalar 3.14159)",
+ round_trip_label<arb::iexpr>("(pi)"));
+}
+
TEST(regloc, round_tripping) {
EXPECT_EQ("(cable 3 0 1)", round_trip_label<arb::region>("(branch 3)"));
EXPECT_EQ("(intersect (tag 1) (intersect (tag 2) (tag 3)))", round_trip_label<arb::region>("(intersect (tag 1) (tag 2) (tag 3))"));
@@ -453,6 +544,7 @@ using place_tuple = std::tuple<arb::locset, arb::placeable, std::string>;
using paint_pair = std::pair<arb::region, arb::paintable>;
using locset_pair = std::pair<std::string, locset>;
using region_pair = std::pair<std::string, region>;
+using iexpr_pair = std::pair<std::string, iexpr>;
std::ostream& operator<<(std::ostream& o, const i_clamp& c) {
o << "(current-clamp (envelope";
@@ -504,6 +596,15 @@ std::ostream& operator<<(std::ostream& o, const synapse& p) {
std::ostream& operator<<(std::ostream& o, const density& p) {
return o << "(density " << p.mech << ')';
}
+template <typename TaggedMech>
+std::ostream& operator<<(std::ostream& o, const scaled_mechanism<TaggedMech>& p) {
+ o << "(scaled-mechanism " << p.t_mech;
+ for (const auto& it: p.scale_expr) {
+ o << " (\"" << it.first << "\" " << it.second << ')';
+ }
+ o << ")";
+ return o;
+}
std::ostream& operator<<(std::ostream& o, const ion_reversal_potential_method& p) {
return o << "(ion-reversal-potential-method \"" << p.ion << "\" " << p.method << ')';
}
@@ -538,6 +639,9 @@ std::ostream& operator<<(std::ostream& o, const locset_pair& p) {
std::ostream& operator<<(std::ostream& o, const region_pair& p) {
return o << "(region-def \"" << p.first << "\" " << p.second << ")";
}
+std::ostream& operator<<(std::ostream& o, const iexpr_pair& p) {
+ return o << "(iexpr-def \"" << p.first << "\" " << p.second << ")";
+}
template <typename T>
std::string to_string(const T& obj) {
@@ -615,6 +719,8 @@ TEST(decor_literals, round_tripping) {
auto paint_literals = {
"(density (mechanism \"hh\"))",
"(density (mechanism \"pas\" (\"g\" 0.02)))",
+ "(scaled-mechanism (density (mechanism \"pas\" (\"g\" 0.02))))",
+ "(scaled-mechanism (density (mechanism \"pas\" (\"g\" 0.02))) (\"g\" (exp (add (radius 2.1) (scalar 3.2)))))",
};
auto default_literals = {
"(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))",
@@ -669,6 +775,7 @@ TEST(decor_expressions, round_tripping) {
"(paint (intersect (cable 2 0 0.5) (region \"axon\")) (ion-external-concentration \"h\" -50.1))",
"(paint (region \"my_region\") (ion-reversal-potential \"na\" 30))",
"(paint (cable 2 0.1 0.4) (density (mechanism \"hh\")))",
+ "(paint (cable 2 0.1 0.4) (scaled-mechanism (density (mechanism \"pas\" (\"g\" 0.02))) (\"g\" (exp (add (distance 2.1 (region \"my_region\")) (scalar 3.2))))))",
"(paint (all) (density (mechanism \"pas\" (\"g\" 0.02))))"
};
auto decorate_default_literals = {
@@ -709,12 +816,20 @@ TEST(label_dict_expressions, round_tripping) {
"(region-def \"my region\" (all))",
"(region-def \"reg42\" (cable 0 0.1 0.9))"
};
+ auto iexpr_def_literals = {
+ "(iexpr-def \"my_iexpr\" (radius 1.2))",
+ "(iexpr-def \"my_iexpr_2\" (iexpr \"my_iexpr\"))",
+ };
+
for (auto l: locset_def_literals) {
EXPECT_EQ(l, round_trip<locset_pair>(l));
}
for (auto l: region_def_literals) {
EXPECT_EQ(l, round_trip<region_pair>(l));
}
+ for (auto l: iexpr_def_literals) {
+ EXPECT_EQ(l, round_trip<iexpr_pair>(l));
+ }
}
TEST(morphology_literals, round_tripping) {
@@ -756,6 +871,18 @@ TEST(decor, round_tripping) {
" (density \n"
" (mechanism \"pas\")))\n"
" (paint \n"
+ " (region \"dend\")\n"
+ " (scaled-mechanism \n"
+ " (density \n"
+ " (mechanism \"pas\"))))\n"
+ " (paint \n"
+ " (region \"soma\")\n"
+ " (scaled-mechanism \n"
+ " (density \n"
+ " (mechanism \"pas\"))\n"
+ " (\"g\" \n"
+ " (radius 2.1))))\n"
+ " (paint \n"
" (region \"soma\")\n"
" (density \n"
" (mechanism \"hh\")))\n"
@@ -788,6 +915,11 @@ TEST(label_dict, round_tripping) {
" (meta-data \n"
" (version \"" + arborio::acc_version() + "\"))\n"
" (label-dict \n"
+ " (iexpr-def \"my_iexpr\" \n"
+ " (log \n"
+ " (mul \n"
+ " (scalar 3.5)\n"
+ " (diameter 4.3))))\n"
" (region-def \"soma\" \n"
" (tag 1))\n"
" (region-def \"dend\" \n"
@@ -963,6 +1095,8 @@ TEST(cable_cell, round_tripping) {
" (point 206.300000 0.000000 0.000000 0.200000)\n"
" 3)))\n"
" (label-dict \n"
+ " (iexpr-def \"my_iexpr\" \n"
+ " (radius 2.1))\n"
" (region-def \"soma\" \n"
" (tag 1))\n"
" (region-def \"dend\" \n"
@@ -981,6 +1115,13 @@ TEST(cable_cell, round_tripping) {
" (density \n"
" (mechanism \"hh\" \n"
" (\"el\" 0.500000))))\n"
+ " (paint \n"
+ " (region \"soma\")\n"
+ " (scaled-mechanism \n"
+ " (density \n"
+ " (mechanism \"pas\"))\n"
+ " (\"g\" \n"
+ " (iexpr \"my_iexpr\"))))\n"
" (place \n"
" (location 0 1)\n"
" (current-clamp \n"
--
GitLab