diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 24e73acd1a904f874c61252fe01c46735c95871d..f0d0924be86d7c48b4e15fea08c81d80fb51bcf4 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -3,6 +3,7 @@
set(arbor_sources
arbexcept.cpp
assert.cpp
+ backends/multicore/fvm.cpp
backends/multicore/mechanism.cpp
backends/multicore/shared_state.cpp
backends/multicore/stimulus.cpp
@@ -60,6 +61,7 @@ set(arbor_sources
if(ARB_WITH_CUDA)
list(APPEND arbor_sources
+ backends/gpu/fvm.cpp
backends/gpu/mechanism.cpp
backends/gpu/mechanism.cu
backends/gpu/shared_state.cpp
diff --git a/arbor/backends/gpu/fvm.cpp b/arbor/backends/gpu/fvm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c07e4b0922f353ef438fc10302c00bf43a4987b
--- /dev/null
+++ b/arbor/backends/gpu/fvm.cpp
@@ -0,0 +1,18 @@
+#include <string>
+
+#include <arbor/mechanism.hpp>
+#include "fvm.hpp"
+#include "mechanism.hpp"
+
+// Provides implementation of backend::mechanism_field_data.
+
+namespace arb {
+namespace gpu {
+
+fvm_value_type* backend::mechanism_field_data(arb::mechanism* mptr, const std::string& field) {
+ arb::gpu::mechanism* m = dynamic_cast<arb::gpu::mechanism*>(mptr);
+ return m? m->field_data(field): nullptr;
+}
+
+} // namespace gpu
+} // namespace arb
diff --git a/arbor/backends/gpu/fvm.hpp b/arbor/backends/gpu/fvm.hpp
index f19a71bd1f8a400028a289c25ebdc0ebc4f2dc0c..b40b0a63a77c3f9d4b744148a9174708a303a8d3 100644
--- a/arbor/backends/gpu/fvm.hpp
+++ b/arbor/backends/gpu/fvm.hpp
@@ -72,6 +72,8 @@ struct backend {
thresholds,
context);
}
+
+ static value_type* mechanism_field_data(arb::mechanism* mptr, const std::string& field);
};
} // namespace gpu
diff --git a/arbor/backends/gpu/mechanism.cpp b/arbor/backends/gpu/mechanism.cpp
index 81d5b29d845b54ce9c9f774c72e92ceb88ec0241..e6c0b46e71f29605ebd06ec0d97b858ba0d0e4b1 100644
--- a/arbor/backends/gpu/mechanism.cpp
+++ b/arbor/backends/gpu/mechanism.cpp
@@ -195,6 +195,14 @@ void mechanism::set_parameter(const std::string& key, const std::vector<fvm_valu
}
}
+fvm_value_type* mechanism::field_data(const std::string& field_var) {
+ if (auto opt_ptr = value_by_key(field_table(), field_var)) {
+ return *opt_ptr.value();
+ }
+
+ return nullptr;
+}
+
void multiply_in_place(fvm_value_type* s, const fvm_index_type* p, int n);
void mechanism::initialize() {
diff --git a/arbor/backends/gpu/mechanism.hpp b/arbor/backends/gpu/mechanism.hpp
index 42c7e57a8085cb02761e7b44d7bfeeaa17788061..402507d690afe8af7638eb9035ef5e4d3464d749 100644
--- a/arbor/backends/gpu/mechanism.hpp
+++ b/arbor/backends/gpu/mechanism.hpp
@@ -55,6 +55,10 @@ public:
void set_parameter(const std::string& key, const std::vector<fvm_value_type>& values) override;
+ // Peek into mechanism state variable; implements arb::gpu::backend::mechanism_field_data.
+ // Returns pointer to GPU memory corresponding to state variable data.
+ fvm_value_type* field_data(const std::string& state_var);
+
void initialize() override;
protected:
diff --git a/arbor/backends/multicore/fvm.cpp b/arbor/backends/multicore/fvm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2432d43d62be2d8a9bba346152d595020a1cb635
--- /dev/null
+++ b/arbor/backends/multicore/fvm.cpp
@@ -0,0 +1,18 @@
+#include <string>
+
+#include <arbor/mechanism.hpp>
+#include "fvm.hpp"
+#include "mechanism.hpp"
+
+// Provides implementation of backend::mechanism_field_data.
+
+namespace arb {
+namespace multicore {
+
+fvm_value_type* backend::mechanism_field_data(arb::mechanism* mptr, const std::string& field) {
+ arb::multicore::mechanism* m = dynamic_cast<arb::multicore::mechanism*>(mptr);
+ return m? m->field_data(field): nullptr;
+}
+
+} // namespace multicore
+} // namespace arb
diff --git a/arbor/backends/multicore/fvm.hpp b/arbor/backends/multicore/fvm.hpp
index 1f9c019ddeff696d121ad5bf156d1fe7ffc5a61c..46767b3490e9bb287a4bc9965b15356614284350 100644
--- a/arbor/backends/multicore/fvm.hpp
+++ b/arbor/backends/multicore/fvm.hpp
@@ -3,6 +3,8 @@
#include <string>
#include <vector>
+#include <arbor/mechanism.hpp>
+
#include "backends/event.hpp"
#include "backends/multicore/matrix_state.hpp"
#include "backends/multicore/multi_event_stream.hpp"
@@ -59,6 +61,8 @@ struct backend {
thresholds,
context);
}
+
+ static fvm_value_type* mechanism_field_data(arb::mechanism* mptr, const std::string& field);
};
} // namespace multicore
diff --git a/arbor/backends/multicore/mechanism.cpp b/arbor/backends/multicore/mechanism.cpp
index 22b9e3c521258b62fdc69656d79974bb313a0112..87597207a70e29dcc0591e06792ff7b46255a352 100644
--- a/arbor/backends/multicore/mechanism.cpp
+++ b/arbor/backends/multicore/mechanism.cpp
@@ -15,6 +15,7 @@
#include "util/maputil.hpp"
#include "util/padded_alloc.hpp"
#include "util/range.hpp"
+#include "util/rangeutil.hpp"
#include "backends/multicore/mechanism.hpp"
#include "backends/multicore/multicore_common.hpp"
@@ -212,5 +213,14 @@ void mechanism::initialize() {
}
}
+fvm_value_type* mechanism::field_data(const std::string& field_var) {
+ if (auto opt_ptr = value_by_key(field_table(), field_var)) {
+ return *opt_ptr.value();
+ }
+
+ return nullptr;
+}
+
+
} // namespace multicore
} // namespace arb
diff --git a/arbor/backends/multicore/mechanism.hpp b/arbor/backends/multicore/mechanism.hpp
index bf7970b01e85841a056b40e9fda181fc8ea486b2..0d1357304e357f9b56f7f1c5736a6957b1f4b7d1 100644
--- a/arbor/backends/multicore/mechanism.hpp
+++ b/arbor/backends/multicore/mechanism.hpp
@@ -65,6 +65,9 @@ public:
void set_parameter(const std::string& key, const std::vector<fvm_value_type>& values) override;
+ // Peek into mechanism state variable; implements arb::multicore::backend::mechanism_field_data.
+ fvm_value_type* field_data(const std::string& state_var);
+
protected:
size_type width_ = 0; // Instance width (number of CVs/sites)
size_type width_padded_ = 0; // Width rounded up to multiple of pad/alignment.
diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp
index 5f77ca1e4edda257e218289cce7b3161c64061da..99d3000ed777180d007c765e781a029b0f851ad2 100644
--- a/arbor/fvm_lowered_cell.hpp
+++ b/arbor/fvm_lowered_cell.hpp
@@ -21,6 +21,18 @@ struct fvm_integration_result {
util::range<const fvm_value_type*> sample_value;
};
+// A sample for a probe may be derived from multiple 'raw' sampled
+// values from the backend.
+//
+// While supported probes are at this point all simple scalar values,
+// fvm_probe_info will be the class that represents the mapping
+// between a single sample result and the back-end raw probe handles.
+
+struct fvm_probe_info {
+ // nullptr => nothing to probe
+ probe_handle raw_handle = nullptr;
+};
+
// Common base class for FVM implementation on host or gpu back-end.
struct fvm_lowered_cell {
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 0f54936bb7220964aa3c8d1ef0b1d577ff52ea95..f35f446eea592b4397049c8d32538418a0b7ab9a 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -136,6 +136,15 @@ private:
void set_gpu() {
if (context_.gpu->has_gpu()) context_.gpu->set_gpu();
}
+
+ // Translate cell probe descriptions into probe handles etc.
+ fvm_probe_info resolve_probe_address(
+ std::size_t cell_idx,
+ const util::any& paddr,
+ const fvm_cv_discretization& D,
+ const fvm_mechanism_data& M,
+ const std::vector<target_handle>& handles,
+ const std::unordered_map<std::string, mechanism*>& mech_instance_by_name);
};
template <typename Backend>
@@ -316,16 +325,16 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
};
}
-template <typename B>
-void fvm_lowered_cell_impl<B>::update_ion_state() {
+template <typename Backend>
+void fvm_lowered_cell_impl<Backend>::update_ion_state() {
state_->ions_init_concentration();
for (auto& m: mechanisms_) {
m->write_ions();
}
}
-template <typename B>
-void fvm_lowered_cell_impl<B>::assert_voltage_bounded(fvm_value_type bound) {
+template <typename Backend>
+void fvm_lowered_cell_impl<Backend>::assert_voltage_bounded(fvm_value_type bound) {
auto v_minmax = state_->voltage_bounds();
if (v_minmax.first>=-bound && v_minmax.second<=bound) {
return;
@@ -337,8 +346,8 @@ void fvm_lowered_cell_impl<B>::assert_voltage_bounded(fvm_value_type bound) {
v_minmax.first<-bound? v_minmax.first: v_minmax.second);
}
-template <typename B>
-void fvm_lowered_cell_impl<B>::initialize(
+template <typename Backend>
+void fvm_lowered_cell_impl<Backend>::initialize(
const std::vector<cell_gid_type>& gids,
const recipe& rec,
std::vector<fvm_index_type>& cell_to_intdom,
@@ -444,6 +453,9 @@ void fvm_lowered_cell_impl<B>::initialize(
target_handles.resize(mech_data.n_target);
+ // Keep track of mechanisms by name for probe lookup.
+ std::unordered_map<std::string, mechanism*> mechptr_by_name;
+
unsigned mech_id = 0;
for (auto& m: mech_data.mechanisms) {
auto& name = m.first;
@@ -498,6 +510,7 @@ void fvm_lowered_cell_impl<B>::initialize(
auto minst = mech_instance(name);
minst.mech->instantiate(mech_id++, *state_, minst.overrides, layout);
+ mechptr_by_name[name] = minst.mech.get();
for (auto& pv: config.param_values) {
minst.mech->set_parameter(pv.first, pv.second);
@@ -526,25 +539,11 @@ void fvm_lowered_cell_impl<B>::initialize(
for (cell_lid_type j: make_span(rec.num_probes(gid))) {
probe_info pi = rec.get_probe({gid, j});
- auto where = any_cast<cell_probe_address>(pi.address);
-
- fvm_size_type cv;
- probe_handle handle;
-
- switch (where.kind) {
- case cell_probe_address::membrane_voltage:
- cv = D.geometry.location_cv(cell_idx, where.location, cv_prefer::cv_empty);
- handle = state_->voltage.data()+cv;
- break;
- case cell_probe_address::membrane_current:
- cv = D.geometry.location_cv(cell_idx, where.location, cv_prefer::cv_nonempty);
- handle = state_->current_density.data()+cv;
- break;
- default:
- throw arbor_internal_error("fvm_lowered_cell: unrecognized probeKind");
- }
+ fvm_probe_info info = resolve_probe_address(cell_idx, pi.address, D, mech_data, target_handles, mechptr_by_name);
- probe_map.insert({pi.id, {handle, pi.tag}});
+ if (info.raw_handle) {
+ probe_map.insert({pi.id, {info.raw_handle, pi.tag}});
+ }
}
}
@@ -554,8 +553,8 @@ void fvm_lowered_cell_impl<B>::initialize(
}
// Get vector of gap_junctions
-template <typename B>
-std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
+template <typename Backend>
+std::vector<fvm_gap_junction> fvm_lowered_cell_impl<Backend>::fvm_gap_junctions(
const std::vector<cable_cell>& cells,
const std::vector<cell_gid_type>& gids,
const recipe& rec, const fvm_cv_discretization& D) {
@@ -599,8 +598,8 @@ std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
return v;
}
-template <typename B>
-fvm_size_type fvm_lowered_cell_impl<B>::fvm_intdom(
+template <typename Backend>
+fvm_size_type fvm_lowered_cell_impl<Backend>::fvm_intdom(
const recipe& rec,
const std::vector<cell_gid_type>& gids,
std::vector<fvm_index_type>& cell_to_intdom) {
@@ -649,4 +648,111 @@ fvm_size_type fvm_lowered_cell_impl<B>::fvm_intdom(
return intdom_id;
}
+template <typename Backend>
+fvm_probe_info fvm_lowered_cell_impl<Backend>::resolve_probe_address(
+ std::size_t cell_idx,
+ const util::any& paddr,
+ const fvm_cv_discretization& D,
+ const fvm_mechanism_data& M,
+ const std::vector<target_handle>& handles,
+ const std::unordered_map<std::string, mechanism*>& mech_instance_by_name)
+{
+ using util::any_cast;
+ using util::value_by_key;
+ using util::binary_search_index;
+ using namespace cv_prefer;
+
+ // Probe address can be one of a number of cable cell-specific
+ // probe types; dispatch on type and compute probe info accordingly.
+ //
+ // Mechanisms or ions not being instantiated at a particular address
+ // is not treated as an error; a request for a state variable that
+ // does not exist on a mechanism is, however.
+
+ auto location_cv = [&](auto* p, cv_prefer::type prefer) -> fvm_index_type {
+ return D.geometry.location_cv(cell_idx, p->location, prefer);
+ };
+
+ struct mechanism_data {
+ const fvm_value_type* state_data = nullptr;
+ unsigned id;
+ explicit operator bool() const { return state_data; }
+ };
+
+ auto lookup_mechanism_data = [&](auto* p) -> mechanism_data {
+ if (mechanism* m = value_by_key(mech_instance_by_name, p->mechanism).value_or(nullptr)) {
+ const fvm_value_type* data = Backend::mechanism_field_data(m, p->state);
+ if (!data) {
+ throw cable_cell_error("no state variable '"+p->state+"' in mechanism '"+p->mechanism+"'");
+ }
+ return {data, m->mechanism_id()};
+ }
+ return {};
+ };
+
+ auto ion_location_index = [&](auto* p) -> util::optional<fvm_index_type> {
+ if (state_->ion_data.count(p->ion)) {
+ auto cv = location_cv(p, cv_nonempty);
+ return util::binary_search_index(M.ions.at(p->ion).cv, cv);
+ }
+ return util::nullopt;
+ };
+
+ if (auto *p = any_cast<cell_probe_membrane_voltage>(&paddr)) {
+ const fvm_value_type* src = state_->voltage.data() + location_cv(p, cv_empty);
+ return fvm_probe_info{src};
+ }
+ else if (auto* p = any_cast<cell_probe_total_ionic_current_density>(&paddr)) {
+ const fvm_value_type* src = state_->current_density.data() + location_cv(p, cv_nonempty);
+ return fvm_probe_info{src};
+ }
+ else if (auto* p = any_cast<cell_probe_density_state>(&paddr)) {
+ if (auto data = lookup_mechanism_data(p)) {
+ auto cv = location_cv(p, cv_nonempty);
+ if (auto opt_i = binary_search_index(M.mechanisms.at(p->mechanism).cv, cv)) {
+ return fvm_probe_info{data.state_data+opt_i.value()};
+ }
+ }
+ return fvm_probe_info{};
+ }
+ else if (auto* p = any_cast<cell_probe_point_state>(&paddr)) {
+ if (p->target>=handles.size()) {
+ return fvm_probe_info{};
+ }
+
+ if (auto data = lookup_mechanism_data(p)) {
+ // Confirm mechanism is actually at this target.
+ const auto& th = handles[p->target];
+ if (th.mech_id == data.id) {
+ return fvm_probe_info{data.state_data+th.mech_index};
+ }
+ }
+ return fvm_probe_info{};
+ }
+ else if (auto* p = any_cast<cell_probe_ion_current_density>(&paddr)) {
+ if (auto opt_i = ion_location_index(p)) {
+ const fvm_value_type* data = state_->ion_data.at(p->ion).iX_.data();
+ return fvm_probe_info{data+opt_i.value()};
+ }
+ return fvm_probe_info{};
+ }
+ else if (auto* p = any_cast<cell_probe_ion_int_concentration>(&paddr)) {
+ if (auto opt_i = ion_location_index(p)) {
+ const fvm_value_type* data = state_->ion_data.at(p->ion).Xi_.data();
+ return fvm_probe_info{data+opt_i.value()};
+ }
+ return fvm_probe_info{};
+ }
+ else if (auto* p = any_cast<cell_probe_ion_ext_concentration>(&paddr)) {
+ if (auto opt_i = ion_location_index(p)) {
+ const fvm_value_type* data = state_->ion_data.at(p->ion).Xo_.data();
+ return fvm_probe_info{data+opt_i.value()};
+ }
+ return fvm_probe_info{};
+ }
+ else {
+ throw cable_cell_error("unrecognized probe address type");
+ }
+}
+
} // namespace arb
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index fd876e1e2182f406f2057588fda10d6a26677128..ed9d3a927a8bd2681dd1e4e7e12b8f00eb3738c0 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -29,14 +29,49 @@ struct lid_range {
begin(b), end(e) {}
};
-// Probe type for cell descriptions.
-struct cell_probe_address {
- enum probe_kind {
- membrane_voltage, membrane_current
- };
+// Each kind of probe has its own type for representing its address,
+// as below:
+// Voltage estimate [mV] at `location`, possibly interpolated.
+struct cell_probe_membrane_voltage {
mlocation location;
- probe_kind kind;
+};
+
+// Total current density [A/m²] across membrane excluding capacitive current at `location`.
+struct cell_probe_total_ionic_current_density {
+ mlocation location;
+};
+
+// Value of state variable `state` in density mechanism `mechanism` in CV at `location`.
+struct cell_probe_density_state {
+ mlocation location;
+ std::string mechanism;
+ std::string state;
+};
+
+// Value of state variable `key` in point mechanism `source` at target `target`.
+struct cell_probe_point_state {
+ cell_lid_type target;
+ std::string mechanism;
+ std::string state;
+};
+
+// Current density [A/m²] across membrane attributed to the ion `source` at `location`.
+struct cell_probe_ion_current_density {
+ mlocation location;
+ std::string ion;
+};
+
+// Ionic internal concentration [mmol/L] of ion `source` at `location`, possibly interpolated.
+struct cell_probe_ion_int_concentration {
+ mlocation location;
+ std::string ion;
+};
+
+// Ionic external concentration [mmol/L] of ion `source` at `location`, possibly interpolated.
+struct cell_probe_ion_ext_concentration {
+ mlocation location;
+ std::string ion;
};
// Forward declare the implementation, for PIMPL.
diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp
index ec9413bc59097e39e59405e55b57995c488f7b11..c8354c26c5a3d6da8fd991a9c6cab0b062567a98 100644
--- a/example/dryrun/dryrun.cpp
+++ b/example/dryrun/dryrun.cpp
@@ -52,7 +52,6 @@ using arb::cell_size_type;
using arb::cell_member_type;
using arb::cell_kind;
using arb::time_type;
-using arb::cell_probe_address;
// Generate a cell.
arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params);
@@ -127,12 +126,10 @@ public:
}
arb::probe_info get_probe(cell_member_type id) const override {
- // Get the appropriate kind for measuring voltage.
- cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
arb::mlocation loc{0, 0.0};
- return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
}
private:
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 8a39dd23bca8ccb49ec2b591e6f86b00d7d27c98..4168aa769a9468a37809041eeacf38c547781b07 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -53,7 +53,6 @@ using arb::cell_size_type;
using arb::cell_member_type;
using arb::cell_kind;
using arb::time_type;
-using arb::cell_probe_address;
// Writes voltage trace as a json file.
void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigned rank);
@@ -100,12 +99,9 @@ public:
}
arb::probe_info get_probe(cell_member_type id) const override {
- // Get the appropriate kind for measuring voltage.
- cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
arb::mlocation loc{0, 1.};
-
- return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
}
arb::util::any get_global_properties(cell_kind k) const override {
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index e5234ec96980dae656103e1697d7ed7389ba58d5..bc392b43919086c6e4cc08338d6009c905b35d97 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -28,7 +28,6 @@ using arb::cell_size_type;
using arb::cell_member_type;
using arb::cell_kind;
using arb::time_type;
-using arb::cell_probe_address;
// Writes voltage trace as a json file.
void write_trace_json(const arb::trace_data<double>& trace);
@@ -125,12 +124,9 @@ public:
arb_assert(id.gid==0); // There is one cell,
arb_assert(id.index==0); // with one probe.
- // Get the appropriate kind for measuring voltage
- cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma
arb::mlocation loc{0, 0.0};
-
- return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
}
};
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index 46b222cc11824717bfef78456451add30988f996..561b09df2dc296c3cf73072d2134d58296f56989 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -53,7 +53,6 @@ using arb::cell_size_type;
using arb::cell_member_type;
using arb::cell_kind;
using arb::time_type;
-using arb::cell_probe_address;
// Writes voltage trace as a json file.
void write_trace_json(const arb::trace_data<double>& trace);
@@ -115,12 +114,9 @@ public:
}
arb::probe_info get_probe(cell_member_type id) const override {
- // Get the appropriate kind for measuring voltage.
- cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
arb::mlocation loc{0, 0.0};
-
- return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
+ return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
}
arb::util::any get_global_properties(arb::cell_kind) const override {
diff --git a/example/single/single.cpp b/example/single/single.cpp
index 3946d9c759357488820b6697e86c8763869ed4ca..4e62819a99477ad9dc6fe2d594732bc8a67f6017 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -39,7 +39,7 @@ struct single_recipe: public arb::recipe {
arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
arb::mlocation mid_soma = {0, 0.5};
- arb::cell_probe_address probe = {mid_soma, arb::cell_probe_address::membrane_voltage};
+ arb::cell_probe_membrane_voltage probe = {mid_soma};
// Probe info consists of: the probe id, a tag value to distinguish this probe
// from others for any attached sampler (unused), and the cell probe address.
diff --git a/python/recipe.cpp b/python/recipe.cpp
index 9d9e342323af35f12b0a0c5641d672b237da87b4..0ce022d900ad799aafa0e5f5ecdcf5278ebd81ac 100644
--- a/python/recipe.cpp
+++ b/python/recipe.cpp
@@ -29,20 +29,14 @@ arb::util::unique_any py_recipe_shim::get_cell_description(arb::cell_gid_type gi
"Python error already thrown");
}
-arb::cell_probe_address::probe_kind probe_kind_from_string(const std::string& name) {
- if (name == "voltage") {
- return arb::cell_probe_address::probe_kind::membrane_voltage;
+arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) {
+ if (kind == "voltage") {
+ return arb::probe_info{id, 0, arb::cell_probe_membrane_voltage{loc}};
}
- else if (name == "current") {
- return arb::cell_probe_address::probe_kind::membrane_current;
+ else if (kind == "ionic current density") {
+ return arb::probe_info{id, 0, arb::cell_probe_total_ionic_current_density{loc}};
}
- else throw pyarb_error(util::pprintf("invalid probe kind: {}, neither voltage nor current", name));
-}
-
-arb::probe_info cable_probe(std::string kind, arb::cell_member_type id, arb::mlocation loc) {
- auto pkind = probe_kind_from_string(kind);
- arb::cell_probe_address probe{loc, pkind};
- return arb::probe_info{id, pkind, probe};
+ else throw pyarb_error(util::pprintf("unrecognized probe kind: {}", kind));
};
std::vector<arb::event_generator> convert_gen(std::vector<pybind11::object> pygens, arb::cell_gid_type gid) {
@@ -180,7 +174,7 @@ void register_recipe(pybind11::module& m) {
// Probes
m.def("cable_probe", &cable_probe,
"Description of a probe at a location on a cable cell with id available for monitoring data of kind "\
- "where kind is one of 'voltage' or 'current'.",
+ "where kind is one of 'voltage' or 'ionic current density'.",
"kind"_a, "id"_a, "location"_a);
pybind11::class_<arb::probe_info> probe(m, "probe");
diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp
index 249e0325e855b5d9253adb80579239ba70f3e190..036477652662c9d4f3e6cf38136d87ddbcaa51d1 100644
--- a/python/single_cell_model.cpp
+++ b/python/single_cell_model.cpp
@@ -116,9 +116,8 @@ struct single_cell_recipe: arb::recipe {
}
// For now only voltage can be selected for measurement.
- auto kind = arb::cell_probe_address::membrane_voltage;
const auto& loc = probes_[probe_id.index].site;
- return arb::probe_info{probe_id, kind, arb::cell_probe_address{loc, kind}};
+ return arb::probe_info{probe_id, 0, arb::cell_probe_membrane_voltage{loc}};
}
// gap junctions
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index d867d294b15d406f4f933ca16e181084b1507887..9f9d6455931cf9e1eb0040cd2bdca5fa1a392805 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -23,6 +23,7 @@ set(test_mechanisms
test_cl_valence
test_ca_read_valence
read_eX
+ write_Xi_Xo
write_multiple_eX
write_eX
read_cai_init
diff --git a/test/unit/mod/fixed_ica_current.mod b/test/unit/mod/fixed_ica_current.mod
index f3887c09d3ee87bd1339f765cc337a9cb67782e8..4756a70b6ad7e2beefe7f4baa2a47cdd64ece5d4 100644
--- a/test/unit/mod/fixed_ica_current.mod
+++ b/test/unit/mod/fixed_ica_current.mod
@@ -2,23 +2,23 @@
NEURON {
SUFFIX fixed_ica_current
- USEION ca WRITE ica VALENCE 2
- RANGE ica_density
+ USEION ca WRITE ica
+ RANGE current_density
}
PARAMETER {
- ica_density = 0
+ current_density = 0
}
ASSIGNED {}
INITIAL {
- ica = ica_density
+ ica = current_density
}
STATE {}
BREAKPOINT {
- ica = ica_density
+ ica = current_density
}
diff --git a/test/unit/mod/write_Xi_Xo.mod b/test/unit/mod/write_Xi_Xo.mod
new file mode 100644
index 0000000000000000000000000000000000000000..fa916d21f901b8340e0021881e8a7ddca94b0282
--- /dev/null
+++ b/test/unit/mod/write_Xi_Xo.mod
@@ -0,0 +1,33 @@
+: Ion concentration writer, with state variable.
+: Used for testing mechanism and ion probes.
+
+NEURON {
+ SUFFIX write_Xi_Xo
+ USEION x WRITE xi, xo
+ GLOBAL xi0, xo0, s0
+}
+
+PARAMETER {
+ xi0 = 1
+ xo0 = 2
+ s0 = 3
+}
+
+ASSIGNED {}
+
+STATE {
+ s
+}
+
+INITIAL {
+ s = s0
+ xi = xi0
+ xo = xo0
+}
+
+BREAKPOINT {
+ s = s0
+ xi = xi0
+ xo = xo0
+}
+
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index bbdc25fd98ff962bb85476bc3a7b68b36e7f8463..c09f25eef696e26127fb0da910b989f197796028 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -515,7 +515,7 @@ TEST(fvm_layout, density_norm_area) {
std::vector<cable_cell> cells{std::move(cell)};
- int ncv = 11; //ncomp + 1
+ int ncv = 11;
std::vector<double> expected_gkbar(ncv, dflt_gkbar);
std::vector<double> expected_gl(ncv, dflt_gl);
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 884d2df185a5c3fc8e720fe2b26c758e74e8bf89..69a8e30d8d90ea697261474b461f7489f4aae6fa 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -449,7 +449,7 @@ TEST(fvm_lowered, derived_mechs) {
cable1d_recipe rec(cells);
rec.catalogue().derive("custom_kin1", "test_kin1", {{"tau", 20.0}});
- cell_probe_address where{{1, 0.3}, cell_probe_address::membrane_current};
+ cell_probe_total_ionic_current_density where{{1, 0.3}};
rec.add_probe(0, 0, where);
rec.add_probe(1, 0, where);
rec.add_probe(2, 0, where);
@@ -662,7 +662,7 @@ TEST(fvm_lowered, ionic_currents) {
const double jca = 1.5;
mechanism_desc m1("fixed_ica_current");
- m1["ica_density"] = jca;
+ m1["current_density"] = jca;
// Mechanism models a well-mixed fixed-depth volume without replenishment,
// giving a linear response to ica over time.
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 7ac92d3c395a403d72ae9ef15f76535b148ac4ad..9b8582e102771117103cdf8cc1bbd4a2a3285337 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -1,25 +1,71 @@
#include "../gtest.h"
-#include <arbor/common_types.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/common_types.hpp>
+#include <arbor/mechcat.hpp>
+#include <arbor/mechinfo.hpp>
+#include <arbor/version.hpp>
+#include <arborenv/gpu_env.hpp>
#include "backends/event.hpp"
#include "backends/multicore/fvm.hpp"
+#ifdef ARB_GPU_ENABLED
+#include "backends/gpu/fvm.hpp"
+#endif
#include "fvm_lowered_cell_impl.hpp"
+#include "memory/cuda_wrappers.hpp"
#include "util/rangeutil.hpp"
#include "common.hpp"
+#include "unit_test_catalogue.hpp"
#include "../common_cells.hpp"
#include "../simple_recipes.hpp"
using namespace arb;
-using fvm_cell = fvm_lowered_cell_impl<multicore::backend>;
-using shared_state = multicore::backend::shared_state;
-ACCESS_BIND(std::unique_ptr<shared_state> fvm_cell::*, fvm_state_ptr, &fvm_cell::state_);
+using multicore_fvm_cell = fvm_lowered_cell_impl<multicore::backend>;
+using multicore_shared_state = multicore::backend::shared_state;
+ACCESS_BIND(std::unique_ptr<multicore_shared_state> multicore_fvm_cell::*, multicore_fvm_state_ptr, &multicore_fvm_cell::state_);
+
+
+template <typename Backend>
+struct backend_access {
+ using fvm_cell = multicore_fvm_cell;
+
+ static multicore_shared_state& state(fvm_cell& cell) {
+ return *(cell.*multicore_fvm_state_ptr).get();
+ }
+
+ static fvm_value_type deref(const fvm_value_type* p) { return *p; }
+};
+
+#ifdef ARB_GPU_ENABLED
+
+using gpu_fvm_cell = fvm_lowered_cell_impl<gpu::backend>;
+using gpu_shared_state = gpu::backend::shared_state;
+ACCESS_BIND(std::unique_ptr<gpu_shared_state> gpu_fvm_cell::*, gpu_fvm_state_ptr, &gpu_fvm_cell::state_);
-TEST(probe, fvm_lowered_cell) {
- execution_context context;
+template <>
+struct backend_access<gpu::backend> {
+ using fvm_cell = gpu_fvm_cell;
+
+ static gpu_shared_state& state(fvm_cell& cell) {
+ return *(cell.*gpu_fvm_state_ptr).get();
+ }
+
+ static fvm_value_type deref(const fvm_value_type* p) {
+ fvm_value_type r;
+ memory::cuda_memcpy_d2h(&r, p, sizeof(r));
+ return r;
+ }
+};
+
+#endif
+
+template <typename Backend>
+void run_v_i_probe_test(const context& ctx) {
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+ auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
cable_cell bs = make_cell_ball_and_stick(false);
@@ -32,15 +78,15 @@ TEST(probe, fvm_lowered_cell) {
mlocation loc1{1, 1};
mlocation loc2{1, 0.3};
- rec.add_probe(0, 10, cell_probe_address{loc0, cell_probe_address::membrane_voltage});
- rec.add_probe(0, 20, cell_probe_address{loc1, cell_probe_address::membrane_voltage});
- rec.add_probe(0, 30, cell_probe_address{loc2, cell_probe_address::membrane_current});
+ rec.add_probe(0, 10, cell_probe_membrane_voltage{loc0});
+ rec.add_probe(0, 20, cell_probe_membrane_voltage{loc1});
+ rec.add_probe(0, 30, cell_probe_total_ionic_current_density{loc2});
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
- fvm_cell lcell(context);
+ fvm_cell lcell(*ctx);
lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
EXPECT_EQ(3u, rec.num_probes(0));
@@ -58,16 +104,16 @@ TEST(probe, fvm_lowered_cell) {
// for the voltage probes (cell membrane potential should
// be constant), and zero for the current probe.
- auto& state = *(lcell.*fvm_state_ptr).get();
+ auto& state = backend_access<Backend>::state(lcell);
auto& voltage = state.voltage;
- auto resting = voltage[0];
+ fvm_value_type resting = voltage[0];
EXPECT_NE(0.0, resting);
// (Probe handles are just pointers in this implementation).
- EXPECT_EQ(resting, *p0);
- EXPECT_EQ(resting, *p1);
- EXPECT_EQ(0.0, *p2);
+ EXPECT_EQ(resting, deref(p0));
+ EXPECT_EQ(resting, deref(p1));
+ EXPECT_EQ(0.0, deref(p2));
// After an integration step, expect voltage probe values
// to differ from resting, and between each other, and
@@ -78,11 +124,376 @@ TEST(probe, fvm_lowered_cell) {
lcell.integrate(0.01, 0.0025, {}, {});
- EXPECT_NE(resting, *p0);
- EXPECT_NE(resting, *p1);
- EXPECT_NE(*p0, *p1);
- EXPECT_NE(0.0, *p2);
+ EXPECT_NE(resting, deref(p0));
+ EXPECT_NE(resting, deref(p1));
+ EXPECT_NE(deref(p0), deref(p1));
+ EXPECT_NE(0.0, deref(p2));
+
+ fvm_value_type v = voltage[0];
+ EXPECT_EQ(v, deref(p0));
+}
+
+template <typename Backend>
+void run_expsyn_g_probe_test(const context& ctx, bool coalesce_synapses = false) {
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+ auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
+
+ const double tau = 2.0;
+ EXPECT_EQ(tau, global_default_catalogue()["expsyn"].parameters.at("tau").default_value);
+
+ // Ball-and-stick cell, two synapses, both in same CV.
+ mlocation loc0{1, 0.8};
+ mlocation loc1{1, 1.0};
+
+ cable_cell bs = make_cell_ball_and_stick(false);
+ bs.place(loc0, "expsyn");
+ bs.place(loc1, "expsyn");
+ bs.default_parameters.discretization = cv_policy_fixed_per_branch(2);
+
+ cable1d_recipe rec(bs, coalesce_synapses);
+ rec.add_probe(0, 10, cell_probe_point_state{0u, "expsyn", "g"});
+ rec.add_probe(0, 20, cell_probe_point_state{1u, "expsyn", "g"});
+
+ std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
+ probe_association_map<probe_handle> probe_map;
+
+ fvm_cell lcell(*ctx);
+ lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
+
+ EXPECT_EQ(2u, rec.num_probes(0));
+ EXPECT_EQ(2u, probe_map.size());
+
+ EXPECT_EQ(10, probe_map.at({0, 0}).tag);
+ EXPECT_EQ(20, probe_map.at({0, 1}).tag);
+
+ probe_handle p0 = probe_map.at({0, 0}).handle;
+ probe_handle p1 = probe_map.at({0, 1}).handle;
+
+ // Expect initial probe values to be intial synapse g == 0.
+
+ EXPECT_EQ(0.0, deref(p0));
+ EXPECT_EQ(0.0, deref(p1));
+
+ if (coalesce_synapses) {
+ // Should be the same raw pointer!
+ EXPECT_EQ(p0, p1);
+ }
+
+ // Integrate to 3 ms, with one event at 1ms to first expsyn weight 0.5,
+ // and another at 2ms to second, weight 1.
+
+ std::vector<deliverable_event> evs = {
+ {1.0, targets[0], 0.5},
+ {2.0, targets[1], 1.0}
+ };
+ const double tfinal = 3.;
+ const double dt = 0.001;
+ lcell.integrate(tfinal, dt, evs, {});
+
+ fvm_value_type g0 = deref(p0);
+ fvm_value_type g1 = deref(p1);
+
+ // Expected value: weight*exp(-(t_final-t_event)/tau).
+ double expected_g0 = 0.5*std::exp(-(tfinal-1.0)/tau);
+ double expected_g1 = 1.0*std::exp(-(tfinal-2.0)/tau);
+
+ const double rtol = 1e-6;
+ if (coalesce_synapses) {
+ EXPECT_TRUE(testing::near_relative(expected_g0+expected_g1, g0, rtol));
+ EXPECT_TRUE(testing::near_relative(expected_g0+expected_g1, g1, rtol));
+ }
+ else {
+ EXPECT_TRUE(testing::near_relative(expected_g0, g0, rtol));
+ EXPECT_TRUE(testing::near_relative(expected_g1, g1, rtol));
+ }
+}
+
+template <typename Backend>
+void run_ion_density_probe_test(const context& ctx) {
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+ auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
+
+ // Use test mechanism write_Xi_Xo to check ion concentration probes and
+ // density mechanism state probes.
+
+ auto cat = make_unit_test_catalogue();
+ cat.derive("write_ca1", "write_Xi_Xo", {{"xi0", 1.25}, {"xo0", 1.5}, {"s0", 1.75}}, {{"x", "ca"}});
+ cat.derive("write_ca2", "write_Xi_Xo", {{"xi0", 2.25}, {"xo0", 2.5}, {"s0", 2.75}}, {{"x", "ca"}});
+ cat.derive("write_na3", "write_Xi_Xo", {{"xi0", 3.25}, {"xo0", 3.5}, {"s0", 3.75}}, {{"x", "na"}});
+
+ // Simple constant diameter cable, 3 CVs.
+
+ cable_cell cable(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cable.default_parameters.discretization = cv_policy_fixed_per_branch(3);
+
+ // Calcium ions everywhere, half written by write_ca1, half by write_ca2.
+ // Sodium ions only on distal half.
+
+ cable.paint(mcable{0, 0., 0.5}, "write_ca1");
+ cable.paint(mcable{0, 0.5, 1.}, "write_ca2");
+ cable.paint(mcable{0, 0.5, 1.}, "write_na3");
+
+ // Place probes in each CV.
+
+ mlocation loc0{0, 0.1};
+ mlocation loc1{0, 0.5};
+ mlocation loc2{0, 0.9};
+
+ cable1d_recipe rec(cable);
+ rec.catalogue() = cat;
+
+ // Probe (0, 0): ca internal on CV 0.
+ rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc0, "ca"});
+ // Probe (0, 1): ca internal on CV 1.
+ rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc1, "ca"});
+ // Probe (0, 2): ca internal on CV 2.
+ rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc2, "ca"});
+
+ // Probe (0, 3): ca external on CV 0.
+ rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc0, "ca"});
+ // Probe (0, 4): ca external on CV 1.
+ rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc1, "ca"});
+ // Probe (0, 5): ca external on CV 2.
+ rec.add_probe(0, 0, cell_probe_ion_ext_concentration{loc2, "ca"});
+
+ // Probe (0, 6): na internal on CV 0.
+ rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc0, "na"});
+ // Probe (0, 7): na internal on CV 2.
+ rec.add_probe(0, 0, cell_probe_ion_int_concentration{loc2, "na"});
+
+ // Probe (0, 8): write_ca2 state 's' in CV 0.
+ rec.add_probe(0, 0, cell_probe_density_state{loc0, "write_ca2", "s"});
+ // Probe (0, 9): write_ca2 state 's' in CV 1.
+ rec.add_probe(0, 0, cell_probe_density_state{loc1, "write_ca2", "s"});
+ // Probe (0, 10): write_ca2 state 's' in CV 2.
+ rec.add_probe(0, 0, cell_probe_density_state{loc2, "write_ca2", "s"});
+
+ std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
+ probe_association_map<probe_handle> probe_map;
+
+ fvm_cell lcell(*ctx);
+ lcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
+
+ // Should be no sodium ion instantiated on CV 0, so probe (0, 6) should
+ // have been silently discared. Similarly, write_ca2 is not instantiated on
+ // CV 0, and so probe (0, 8) should have been discarded. All other probes
+ // should be in the map.
+
+ EXPECT_EQ(11u, rec.num_probes(0));
+ EXPECT_EQ(9u, probe_map.size());
+
+ probe_handle ca_int_cv0 = probe_map.at({0, 0}).handle;
+ probe_handle ca_int_cv1 = probe_map.at({0, 1}).handle;
+ probe_handle ca_int_cv2 = probe_map.at({0, 2}).handle;
+ probe_handle ca_ext_cv0 = probe_map.at({0, 3}).handle;
+ probe_handle ca_ext_cv1 = probe_map.at({0, 4}).handle;
+ probe_handle ca_ext_cv2 = probe_map.at({0, 5}).handle;
+ EXPECT_EQ(0u, probe_map.count({0, 6}));
+ probe_handle na_int_cv2 = probe_map.at({0, 7}).handle;
+ EXPECT_EQ(0u, probe_map.count({0, 8}));
+ probe_handle write_ca2_s_cv1 = probe_map.at({0, 9}).handle;
+ probe_handle write_ca2_s_cv2 = probe_map.at({0, 10}).handle;
+
+ // Ion concentrations should have been written in initialization.
+ // For CV 1, calcium concentration should be mean of the two values
+ // from write_ca1 and write_ca2.
+
+ EXPECT_EQ(1.25, deref(ca_int_cv0));
+ EXPECT_DOUBLE_EQ((1.25+2.25)/2., deref(ca_int_cv1));
+ EXPECT_EQ(2.25, deref(ca_int_cv2));
+
+ EXPECT_EQ(1.5, deref(ca_ext_cv0));
+ EXPECT_DOUBLE_EQ((1.5+2.5)/2., deref(ca_ext_cv1));
+ EXPECT_EQ(2.5, deref(ca_ext_cv2));
+
+ EXPECT_EQ(3.25, deref(na_int_cv2));
+
+ // State variable in write_ca2 should be the same in both CV 1 and 2.
+ // The raw handles should be different addresses, however.
+
+ EXPECT_EQ(2.75, deref(write_ca2_s_cv1));
+ EXPECT_EQ(2.75, deref(write_ca2_s_cv2));
+ EXPECT_NE(write_ca2_s_cv1, write_ca2_s_cv2);
+}
+
+template <typename Backend>
+void run_ion_current_probe_test(const context& ctx) {
+ using fvm_cell = typename backend_access<Backend>::fvm_cell;
+ auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
+
+ // Use test mechanism fixed_ica_current, and a derived mechanism for sodium, to
+ // write to specific ion currents.
+
+ auto cat = make_unit_test_catalogue();
+ cat.derive("fixed_ina_current", "fixed_ica_current", {}, {{"ca", "na"}});
+
+ cable_cell cells[2];
+
+ // Simple constant diameter cable, 3 CVs.
+
+ cells[0] = cable_cell(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cells[0].default_parameters.discretization = cv_policy_fixed_per_branch(3);
+
+ // Calcium ions everywhere, half with current density jca0, half with jca1.
+ // Sodium ions only on distal half, with current densitry jna1.
+
+ const double jca0 = 1.5; // [A/m²]
+ const double jca1 = 2.0;
+ const double jna1 = 2.5;
+
+ // Scaling factor 0.1 is to convert our current densities in [A/m²] to NMODL units [mA/cm²].
+
+ cells[0].paint(mcable{0, 0., 0.5}, mechanism_desc("fixed_ica_current").set("current_density", 0.1*jca0));
+ cells[0].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ica_current").set("current_density", 0.1*jca1));
+ cells[0].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ina_current").set("current_density", 0.1*jna1));
+
+ // Make a second cable cell, with same layout but 3 times the current.
+
+ cells[1] = cable_cell(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cells[1].default_parameters.discretization = cv_policy_fixed_per_branch(3);
+
+ cells[1].paint(mcable{0, 0., 0.5}, mechanism_desc("fixed_ica_current").set("current_density", 0.3*jca0));
+ cells[1].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ica_current").set("current_density", 0.3*jca1));
+ cells[1].paint(mcable{0, 0.5, 1.}, mechanism_desc("fixed_ina_current").set("current_density", 0.3*jna1));
+
+ // Place probes in each CV on cell 0, plus one in the last CV on cell 1.
+
+ mlocation loc0{0, 0.1};
+ mlocation loc1{0, 0.5};
+ mlocation loc2{0, 0.9};
+
+ cable1d_recipe rec(cells);
+ rec.catalogue() = cat;
+
+ // Probe (0, 0): ica on CV 0.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc0, "ca"});
+ // Probe (0, 1): ica on CV 1.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc1, "ca"});
+ // Probe (0, 2): ica on CV 2.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc2, "ca"});
+
+ // Probe (0, 3): ina on CV 0.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc0, "na"});
+ // Probe (0, 4): ina on CV 1.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc1, "na"});
+ // Probe (0, 5): ina on CV 2.
+ rec.add_probe(0, 0, cell_probe_ion_current_density{loc2, "na"});
+
+ // Probe (0, 6): total ion current density on CV 0.
+ rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc0});
+ // Probe (0, 7): total ion current density on CV 1.
+ rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc1});
+ // Probe (0, 8): total ion current density on CV 2.
+ rec.add_probe(0, 0, cell_probe_total_ionic_current_density{loc2});
+
+ // Probe (1, 0): ica on CV 5 (CV 2 of cell 1).
+ rec.add_probe(1, 0, cell_probe_ion_current_density{loc2, "ca"});
+ // Probe (1, 1): total ion current density on CV 5 (CV 2 of cell 1).
+ rec.add_probe(1, 0, cell_probe_total_ionic_current_density{loc2});
+
+ std::vector<target_handle> targets;
+ std::vector<fvm_index_type> cell_to_intdom;
+ probe_association_map<probe_handle> probe_map;
+
+ fvm_cell lcell(*ctx);
+ lcell.initialize({0, 1}, rec, cell_to_intdom, targets, probe_map);
+
+ // Should be no sodium ion instantiated on CV 0, so probe (0, 3) should
+ // have been silently discared.
+
+ EXPECT_EQ(9u, rec.num_probes(0));
+ EXPECT_EQ(2u, rec.num_probes(1));
+ EXPECT_EQ(10u, probe_map.size());
+
+ probe_handle ica_cv0 = probe_map.at({0, 0}).handle;
+ probe_handle ica_cv1 = probe_map.at({0, 1}).handle;
+ probe_handle ica_cv2 = probe_map.at({0, 2}).handle;
+ EXPECT_EQ(0u, probe_map.count({0, 3}));
+ probe_handle ina_cv1 = probe_map.at({0, 4}).handle;
+ probe_handle ina_cv2 = probe_map.at({0, 5}).handle;
+ probe_handle i_cv0 = probe_map.at({0, 6}).handle;
+ probe_handle i_cv1 = probe_map.at({0, 7}).handle;
+ probe_handle i_cv2 = probe_map.at({0, 8}).handle;
+
+ probe_handle ica_cv5 = probe_map.at({1, 0}).handle;
+ probe_handle i_cv5 = probe_map.at({1, 1}).handle;
+
+ // Integrate cell for a bit, and check that currents add up as we expect.
+
+ lcell.integrate(0.01, 0.0025, {}, {});
+
+ EXPECT_DOUBLE_EQ(jca0, deref(ica_cv0));
+ EXPECT_DOUBLE_EQ((jca0+jca1)/2, deref(ica_cv1));
+ EXPECT_DOUBLE_EQ(jca1, deref(ica_cv2));
+
+ EXPECT_DOUBLE_EQ(jna1/2, deref(ina_cv1));
+ EXPECT_DOUBLE_EQ(jna1, deref(ina_cv2));
+
+ EXPECT_DOUBLE_EQ(jca0, deref(i_cv0));
+ EXPECT_DOUBLE_EQ(jna1/2+jca0/2+jca1/2, deref(i_cv1));
+ EXPECT_DOUBLE_EQ(jna1+jca1, deref(i_cv2));
+
+ // Currents on cell 1 should be 3 times those on cell 0.
+
+ EXPECT_DOUBLE_EQ(jca1*3, deref(ica_cv5));
+ EXPECT_DOUBLE_EQ((jna1+jca1)*3, deref(i_cv5));
+}
+
+TEST(probe, multicore_v_i) {
+ context ctx = make_context();
+ run_v_i_probe_test<multicore::backend>(ctx);
+}
+
+TEST(probe, multicore_expsyn_g) {
+ context ctx = make_context();
+ SCOPED_TRACE("uncoalesced synapses");
+ run_expsyn_g_probe_test<multicore::backend>(ctx, false);
+ SCOPED_TRACE("coalesced synapses");
+ run_expsyn_g_probe_test<multicore::backend>(ctx, true);
+}
+
+TEST(probe, multicore_ion_conc) {
+ context ctx = make_context();
+ run_ion_density_probe_test<multicore::backend>(ctx);
+}
+
+TEST(probe, multicore_ion_currents) {
+ context ctx = make_context();
+ run_ion_current_probe_test<multicore::backend>(ctx);
+}
+
+#ifdef ARB_GPU_ENABLED
+TEST(probe, gpu_v_i) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_v_i_probe_test<gpu::backend>(ctx);
+ }
+}
+
+TEST(probe, gpu_expsyn_g) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ SCOPED_TRACE("uncoalesced synapses");
+ run_expsyn_g_probe_test<gpu::backend>(ctx, false);
+ SCOPED_TRACE("coalesced synapses");
+ run_expsyn_g_probe_test<gpu::backend>(ctx, true);
+ }
+}
+
+TEST(probe, multicore_ion_conc) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_ion_density_probe_test<gpu::backend>(ctx);
+ }
+}
- EXPECT_EQ(voltage[0], *p0);
+TEST(probe, multicore_ion_currents) {
+ context ctx = make_context(proc_allocation{1, arbenv::default_gpu()});
+ if (has_gpu(ctx)) {
+ run_ion_current_probe_test<gpu::backend>(ctx);
+ }
}
+#endif
diff --git a/test/unit/unit_test_catalogue.cpp b/test/unit/unit_test_catalogue.cpp
index 0d893f288939b4ec65796e413c1f134429b57a0d..8d01b374326be96af215daa0ccc5856d96dea500 100644
--- a/test/unit/unit_test_catalogue.cpp
+++ b/test/unit/unit_test_catalogue.cpp
@@ -29,6 +29,7 @@
#include "mechanisms/test_cl_valence.hpp"
#include "mechanisms/test_ca_read_valence.hpp"
#include "mechanisms/read_eX.hpp"
+#include "mechanisms/write_Xi_Xo.hpp"
#include "mechanisms/write_multiple_eX.hpp"
#include "mechanisms/write_eX.hpp"
#include "mechanisms/read_cai_init.hpp"
@@ -74,6 +75,7 @@ mechanism_catalogue make_unit_test_catalogue() {
ADD_MECH(cat, test_cl_valence)
ADD_MECH(cat, test_ca_read_valence)
ADD_MECH(cat, read_eX)
+ ADD_MECH(cat, write_Xi_Xo)
ADD_MECH(cat, write_multiple_eX)
ADD_MECH(cat, write_eX)
ADD_MECH(cat, read_cai_init)