diff --git a/arbor/cell_group.hpp b/arbor/cell_group.hpp
index 4c73352c48810daf5ed32ab2b853a8dd073c2387..c624946b473ee24b45d4113556476ddb32804fb0 100644
--- a/arbor/cell_group.hpp
+++ b/arbor/cell_group.hpp
@@ -37,6 +37,13 @@ public:
virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function, sampling_policy) = 0;
virtual void remove_sampler(sampler_association_handle) = 0;
virtual void remove_all_samplers() = 0;
+
+ // Probe metadata queries might also be called while a simulation is running, and so should
+ // also be thread-safe.
+
+ virtual std::vector<probe_metadata> get_probe_metadata(cell_member_type) const {
+ return {};
+ }
};
using cell_group_ptr = std::unique_ptr<cell_group>;
diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp
index 83f934586ae4e7a6a11a16390d9fd7548be220ff..9a4f397219f172928d05e01f53410e694cde8e00 100644
--- a/arbor/fvm_lowered_cell.hpp
+++ b/arbor/fvm_lowered_cell.hpp
@@ -11,6 +11,7 @@
#include <arbor/fvm_types.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/util/any_ptr.hpp>
#include <arbor/util/variant.hpp>
#include "backends/event.hpp"
@@ -40,12 +41,18 @@ struct fvm_integration_result {
struct fvm_probe_scalar {
probe_handle raw_handles[1] = {nullptr};
util::variant<mlocation, cable_probe_point_info> metadata;
+
+ util::any_ptr get_metadata_ptr() const {
+ return util::visit([](const auto& x) -> util::any_ptr { return &x; }, metadata);
+ }
};
struct fvm_probe_interpolated {
probe_handle raw_handles[2] = {nullptr, nullptr};
double coef[2];
mlocation metadata;
+
+ util::any_ptr get_metadata_ptr() const { return &metadata; }
};
struct fvm_probe_multi {
@@ -56,6 +63,10 @@ struct fvm_probe_multi {
raw_handles.shrink_to_fit();
util::visit([](auto& v) { v.shrink_to_fit(); }, metadata);
}
+
+ util::any_ptr get_metadata_ptr() const {
+ return util::visit([](const auto& x) -> util::any_ptr { return &x; }, metadata);
+ }
};
struct fvm_probe_weighted_multi {
@@ -68,6 +79,8 @@ struct fvm_probe_weighted_multi {
weight.shrink_to_fit();
metadata.shrink_to_fit();
}
+
+ util::any_ptr get_metadata_ptr() const { return &metadata; }
};
// Trans-membrane currents require special handling!
@@ -88,11 +101,16 @@ struct fvm_probe_membrane_currents {
weight.shrink_to_fit();
cv_cables_divs.shrink_to_fit();
}
+
+ util::any_ptr get_metadata_ptr() const { return &metadata; }
};
struct missing_probe_info {
// dummy data...
std::array<probe_handle, 0> raw_handles;
+ void* metadata = nullptr;
+
+ util::any_ptr get_metadata_ptr() const { return util::any_ptr{}; }
};
struct fvm_probe_data {
@@ -123,6 +141,10 @@ struct fvm_probe_data {
info));
}
+ util::any_ptr get_metadata_ptr() const {
+ return util::visit([](const auto& i) -> util::any_ptr { return i.get_metadata_ptr(); }, info);
+ }
+
sample_size_type n_raw() const { return raw_handle_range().size(); }
explicit operator bool() const { return !util::get_if<missing_probe_info>(info); }
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index 5d7d35e09512b0b13589f83f111e827ddc22dae1..170e1fbb15191c5c5c0594eca11556fe303768fd 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -159,6 +159,11 @@ locset uniform(region reg, unsigned left, unsigned right, uint64_t seed);
// Proportional location on every branch.
locset on_branches(double pos);
+// Proportional locations on each component:
+// For each component C of the region, find locations L
+// s.t. dist(h, L) = r * max {dist(h, t) | t is a distal point in C}.
+locset on_components(double relpos, region reg);
+
// Support of a locset (x s.t. x in locset).
locset support(locset);
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index 5e60d945a372d1763ba28feda44709ceb4d68871..ebdedff41db03c70244a1b7ef84fc323fcfa3e01 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -133,10 +133,10 @@ region branch(msize_t);
// Region with all segments with segment tag id.
region tagged(int id);
-// Region with all segments distal from another region
+// Region up to `distance` distal from points in `start`.
region distal_interval(locset start, double distance);
-// Region with all segments proximal from another region
+// Region up to `distance` proximal from points in `start`.
region proximal_interval(locset end, double distance);
// Region with all segments with radius less than/less than or equal to r
diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp
index 3dc177a4467455fb8d1c7ee2c7b2cad6afdda4e0..c26d022577f47837538686ff68953c79c8b3d74c 100644
--- a/arbor/include/arbor/simulation.hpp
+++ b/arbor/include/arbor/simulation.hpp
@@ -39,6 +39,10 @@ public:
void remove_all_samplers();
+ // Return probe metadata, one entry per probe associated with supplied probe id,
+ // or an empty vector if no local match for probe id.
+ std::vector<probe_metadata> get_probe_metadata(cell_member_type probe_id) const;
+
std::size_t num_spikes() const;
// Set event binning policy on all our groups.
diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp
index f5d85ff572bac32d9c3672726828d560b74ffe23..2b6950d2e0bbce24be3205e988d2ba22f1e9e1fb 100644
--- a/arbor/mc_cell_group.cpp
+++ b/arbor/mc_cell_group.cpp
@@ -62,8 +62,8 @@ void mc_cell_group::reset() {
spikes_.clear();
sample_events_.clear();
- for (auto &assoc: sampler_map_) {
- assoc.sched.reset();
+ for (auto &entry: sampler_map_) {
+ entry.second.sched.reset();
}
for (auto& b: binners_) {
@@ -133,10 +133,7 @@ void run_samples(
sample_records.push_back(sample_record{time_type(raw_times[i]), &raw_samples[i]});
}
- // Metadata may be an mlocation or cell_lid_type (for point mechanism state probes).
- util::visit(
- [&](auto& metadata) { sc.sampler({sc.probe_id, sc.tag, sc.index, &metadata}, n_sample, sample_records.data()); },
- p.metadata);
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
void run_samples(
@@ -166,7 +163,7 @@ void run_samples(
sample_records.push_back(sample_record{time_type(raw_times[offset]), &ctmp[j]});
}
- sc.sampler({sc.probe_id, sc.tag, sc.index, &p.metadata}, n_sample, sample_records.data());
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
void run_samples(
@@ -196,10 +193,7 @@ void run_samples(
sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
}
- // Metadata may be an mlocation_list or mcable_list.
- util::visit(
- [&](auto& metadata) { sc.sampler({sc.probe_id, sc.tag, sc.index, &metadata}, n_sample, sample_records.data()); },
- p.metadata);
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
void run_samples(
@@ -242,8 +236,7 @@ void run_samples(
sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
}
- // (Unlike fvm_probe_multi, we only have mcable_list metadata.)
- sc.sampler({sc.probe_id, sc.tag, sc.index, &p.metadata}, n_sample, sample_records.data());
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
void run_samples(
@@ -305,7 +298,7 @@ void run_samples(
sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
}
- sc.sampler({sc.probe_id, sc.tag, sc.index, &p.metadata}, n_sample, sample_records.data());
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
// Generic run_samples dispatches on probe info variant type.
@@ -391,37 +384,44 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
std::vector<deliverable_event> exact_sampling_events;
- for (auto& sa: sampler_map_) {
- auto sample_times = util::make_range(sa.sched.events(tstart, ep.tfinal));
- if (sample_times.empty()) {
- continue;
- }
-
- sample_size_type n_times = sample_times.size();
- max_samples_per_call = std::max(max_samples_per_call, n_times);
+ {
+ std::lock_guard<std::mutex> guard(sampler_mex_);
- for (cell_member_type pid: sa.probe_ids) {
- auto cell_index = gid_index_map_.at(pid.gid);
+ for (auto& sm_entry: sampler_map_) {
+ // Ignore sampler_association_handle, just need the association itself.
+ sampler_association& sa = sm_entry.second;
- probe_tag tag = probe_map_.tag.at(pid);
- unsigned index = 0;
- for (const fvm_probe_data& pdata: probe_map_.data_on(pid)) {
- call_info.push_back({sa.sampler, pid, tag, index++, &pdata, n_samples, n_samples + n_times*pdata.n_raw()});
- auto intdom = cell_to_intdom_[cell_index];
+ auto sample_times = util::make_range(sa.sched.events(tstart, ep.tfinal));
+ if (sample_times.empty()) {
+ continue;
+ }
- for (auto t: sample_times) {
- for (probe_handle h: pdata.raw_handle_range()) {
- sample_event ev{t, (cell_gid_type)intdom, {h, n_samples++}};
- sample_events.push_back(ev);
- }
- if (sa.policy==sampling_policy::exact) {
- target_handle h(-1, 0, intdom);
- exact_sampling_events.push_back({t, h, 0.f});
+ sample_size_type n_times = sample_times.size();
+ max_samples_per_call = std::max(max_samples_per_call, n_times);
+
+ for (cell_member_type pid: sa.probe_ids) {
+ auto cell_index = gid_index_map_.at(pid.gid);
+
+ probe_tag tag = probe_map_.tag.at(pid);
+ unsigned index = 0;
+ for (const fvm_probe_data& pdata: probe_map_.data_on(pid)) {
+ call_info.push_back({sa.sampler, pid, tag, index++, &pdata, n_samples, n_samples + n_times*pdata.n_raw()});
+ auto intdom = cell_to_intdom_[cell_index];
+
+ for (auto t: sample_times) {
+ for (probe_handle h: pdata.raw_handle_range()) {
+ sample_event ev{t, (cell_gid_type)intdom, {h, n_samples++}};
+ sample_events.push_back(ev);
+ }
+ if (sa.policy==sampling_policy::exact) {
+ target_handle h(-1, 0, intdom);
+ exact_sampling_events.push_back({t, h, 0.f});
+ }
}
}
}
+ arb_assert(n_samples==call_info.back().end_offset);
}
- arb_assert(n_samples==call_info.back().end_offset);
}
// Sort exact sampling events into staged events for delivery.
@@ -481,20 +481,45 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
void mc_cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probe_ids,
schedule sched, sampler_function fn, sampling_policy policy)
{
+ std::lock_guard<std::mutex> guard(sampler_mex_);
+
std::vector<cell_member_type> probeset =
util::assign_from(util::filter(util::keys(probe_map_.tag), probe_ids));
if (!probeset.empty()) {
- sampler_map_.add(h, sampler_association{std::move(sched), std::move(fn), std::move(probeset), policy});
+ auto result = sampler_map_.insert({h, sampler_association{std::move(sched), std::move(fn), std::move(probeset), policy}});
+ arb_assert(result.second);
}
}
void mc_cell_group::remove_sampler(sampler_association_handle h) {
- sampler_map_.remove(h);
+ std::lock_guard<std::mutex> guard(sampler_mex_);
+ sampler_map_.erase(h);
}
void mc_cell_group::remove_all_samplers() {
+ std::lock_guard<std::mutex> guard(sampler_mex_);
sampler_map_.clear();
}
+std::vector<probe_metadata> mc_cell_group::get_probe_metadata(cell_member_type probe_id) const {
+ // Probe associations are fixed after construction, so we do not need to grab the mutex.
+
+ util::optional<probe_tag> maybe_tag = util::value_by_key(probe_map_.tag, probe_id);
+ if (!maybe_tag) {
+ return {};
+ }
+
+ auto data = probe_map_.data_on(probe_id);
+
+ std::vector<probe_metadata> result;
+ result.reserve(data.size());
+ unsigned index = 0;
+ for (const fvm_probe_data& item: data) {
+ result.push_back(probe_metadata{probe_id, *maybe_tag, index++, item.get_metadata_ptr()});
+ }
+
+ return result;
+}
+
} // namespace arb
diff --git a/arbor/mc_cell_group.hpp b/arbor/mc_cell_group.hpp
index 8285a4f0648f1c5239b98cf12ae32e9f9df00821..e342dbf13633d981677f5a06e9c07b5ad3d22145 100644
--- a/arbor/mc_cell_group.hpp
+++ b/arbor/mc_cell_group.hpp
@@ -3,6 +3,7 @@
#include <cstdint>
#include <functional>
#include <iterator>
+#include <mutex>
#include <unordered_map>
#include <vector>
@@ -57,6 +58,8 @@ public:
void remove_all_samplers() override;
+ std::vector<probe_metadata> get_probe_metadata(cell_member_type probe_id) const override;
+
private:
// List of the gids of the cells in the group.
std::vector<cell_gid_type> gids_;
@@ -94,6 +97,9 @@ private:
// Collection of samplers to be run against probes in this group.
sampler_association_map sampler_map_;
+ // Mutex for thread-safe access to sampler associations.
+ std::mutex sampler_mex_;
+
// Lookup table for target ids -> local target handle indices.
std::vector<std::size_t> target_handle_divisions_;
};
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index dca7c313a770e23c39b0e86dd6ae56ee57892f84..365c66a66c97f558a0069c2a6176cb3c4b1226bf 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -260,14 +260,13 @@ mlocation_list thingify_(const boundary_& n, const mprovider& p) {
mlocation_list L;
for (const mextent& comp: comps) {
- mlocation_list proximal_set, distal_set;
+ arb_assert(!comp.empty());
+ arb_assert(thingify_(most_proximal_{region{comp}}, p).size()==1u);
- for (const mcable& c: comp) {
- proximal_set.push_back({c.branch, c.prox_pos});
- distal_set.push_back({c.branch, c.dist_pos});
- }
+ mlocation_list distal_set;
+ util::assign(distal_set, util::transform_view(comp, [](auto c) { return dist_loc(c); }));
- L = sum(L, minset(p.morphology(), proximal_set));
+ L = sum(L, {prox_loc(comp.front())});
L = sum(L, maxset(p.morphology(), distal_set));
}
return support(std::move(L));
@@ -297,13 +296,16 @@ mlocation_list thingify_(const cboundary_& n, const mprovider& p) {
mlocation_list L;
for (const mextent& comp: comps) {
- mlocation_list proximal_set, distal_set;
-
mextent ccomp = thingify(reg::complete(comp), p);
- for (const mcable& c: ccomp.cables()) {
- proximal_set.push_back({c.branch, c.prox_pos});
- distal_set.push_back({c.branch, c.dist_pos});
- }
+
+ // Note: if component contains the head of a top-level cable,
+ // the completion might not be connected (!).
+
+ mlocation_list proximal_set;
+ util::assign(proximal_set, util::transform_view(ccomp, [](auto c) { return prox_loc(c); }));
+
+ mlocation_list distal_set;
+ util::assign(distal_set, util::transform_view(ccomp, [](auto c) { return dist_loc(c); }));
L = sum(L, minset(p.morphology(), proximal_set));
L = sum(L, maxset(p.morphology(), distal_set));
@@ -315,6 +317,83 @@ std::ostream& operator<<(std::ostream& o, const cboundary_& x) {
return o << "(cboundary " << x.reg << ")";
}
+// Proportional on components of a region.
+
+struct on_components_: locset_tag {
+ explicit on_components_(double relpos, region reg):
+ relpos(relpos), reg(std::move(reg)) {}
+ double relpos;
+ region reg;
+};
+
+locset on_components(double relpos, region reg) {
+ return locset(on_components_(relpos, std::move(reg)));
+}
+
+mlocation_list thingify_(const on_components_& n, const mprovider& p) {
+ if (n.relpos<0 || n.relpos>1) {
+ return {};
+ }
+
+ std::vector<mextent> comps = components(p.morphology(), thingify(n.reg, p));
+ std::vector<mlocation> L;
+
+ for (const mextent& comp: comps) {
+ arb_assert(!comp.empty());
+ arb_assert(thingify_(most_proximal_{region{comp}}, p).size()==1u);
+
+ mlocation prox = prox_loc(comp.front());
+ auto d_from_prox = [&](mlocation x) { return p.embedding().integrate_length(prox, x); };
+
+ if (n.relpos==0) {
+ L.push_back(prox);
+ }
+ else if (n.relpos==1) {
+ double diameter = 0;
+ mlocation_list most_distal = {prox};
+
+ for (mcable c: comp) {
+ mlocation x = dist_loc(c);
+ double d = d_from_prox(x);
+
+ if (d>diameter) {
+ most_distal = {x};
+ diameter = d;
+ }
+ else if (d==diameter) {
+ most_distal.push_back(x);
+ }
+ }
+
+ util::append(L, most_distal);
+ }
+ else {
+ double diameter = util::max_value(util::transform_view(comp,
+ [&](auto c) { return d_from_prox(dist_loc(c)); }));
+
+ double d = n.relpos*diameter;
+ for (mcable c: comp) {
+ double d0 = d_from_prox(prox_loc(c));
+ double d1 = d_from_prox(dist_loc(c));
+
+ if (d0<=d && d<=d1) {
+ double s = d0==d1? 0: (d-d0)/(d1-d0);
+ s = std::min(1.0, std::fma(s, c.dist_pos-c.prox_pos, c.prox_pos));
+ L.push_back(mlocation{c.branch, s});
+ }
+ }
+
+ }
+ }
+
+ util::sort(L);
+ return L;
+}
+
+std::ostream& operator<<(std::ostream& o, const on_components_& x) {
+ return o << "(on_components " << x.relpos << " " << x.reg << ")";
+}
+
// Uniform locset.
struct uniform_ {
diff --git a/arbor/sampler_map.hpp b/arbor/sampler_map.hpp
index 50d0ff9b4a00df551752a65d2810a206e5fb2659..495a676d4b3030586f4b4f8c01e651e0c726f2c6 100644
--- a/arbor/sampler_map.hpp
+++ b/arbor/sampler_map.hpp
@@ -1,20 +1,15 @@
#pragma once
-/*
- * Helper classes for managing sampler/schedule associations in
- * cell group classes (see sampling_api doc).
- */
+// Helper classes for managing sampler/schedule associations in
+// cell group classes (see sampling_api doc).
-#include <functional>
-#include <mutex>
#include <unordered_map>
+#include <vector>
#include <arbor/common_types.hpp>
#include <arbor/sampling.hpp>
#include <arbor/schedule.hpp>
-#include "util/transform.hpp"
-
namespace arb {
// An association between a samplers, schedule, and set of probe ids, as provided
@@ -27,38 +22,6 @@ struct sampler_association {
sampling_policy policy;
};
-// Maintain a set of associations paired with handles used for deletion.
-
-class sampler_association_map {
-public:
- void add(sampler_association_handle h, sampler_association assoc) {
- std::lock_guard<std::mutex> lock(m_);
- map_.insert({h, std::move(assoc)});
- }
-
- void remove(sampler_association_handle h) {
- std::lock_guard<std::mutex> lock(m_);
- map_.erase(h);
- }
-
- void clear() {
- std::lock_guard<std::mutex> lock(m_);
- map_.clear();
- }
-
-private:
- using assoc_map = std::unordered_map<sampler_association_handle, sampler_association>;
- assoc_map map_;
- std::mutex m_;
-
- static sampler_association& second(assoc_map::value_type& p) { return p.second; }
- auto assoc_view() { return util::transform_view(map_, &sampler_association_map::second); }
-
-public:
- // Range-like view presents just the associations, omitting the handles.
-
- auto begin() { return assoc_view().begin(); }
- auto end() { return assoc_view().end(); }
-};
+using sampler_association_map = std::unordered_map<sampler_association_handle, sampler_association>;
} // namespace arb
diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp
index 7e32ec774cc7218f03c68849fd714d5e94debcf3..726e798ed516be5e8d9b6a5903fad0f6bb8031c6 100644
--- a/arbor/simulation.cpp
+++ b/arbor/simulation.cpp
@@ -65,6 +65,8 @@ public:
void remove_all_samplers();
+ std::vector<probe_metadata> get_probe_metadata(cell_member_type) const;
+
std::size_t num_spikes() const {
return communicator_.num_spikes();
}
@@ -98,7 +100,11 @@ private:
std::unique_ptr<spike_double_buffer> local_spikes_;
// Hash table for looking up the the local index of a cell with a given gid
- std::unordered_map<cell_gid_type, cell_size_type> gid_to_local_;
+ struct gid_local_info {
+ cell_size_type cell_index;
+ cell_size_type group_index;
+ };
+ std::unordered_map<cell_gid_type, gid_local_info> gid_to_local_;
communicator communicator_;
@@ -146,16 +152,18 @@ simulation_state::simulation_state(
pending_events_.resize(num_local_cells);
event_generators_.resize(num_local_cells);
- cell_local_size_type lidx = 0;
+ cell_size_type lidx = 0;
+ cell_size_type grpidx = 0;
for (const auto& group_info: decomp.groups) {
for (auto gid: group_info.gids) {
// Store mapping of gid to local cell index.
- gid_to_local_[gid] = lidx;
+ gid_to_local_[gid] = gid_local_info{lidx, grpidx};
// Set up the event generators for cell gid.
event_generators_[lidx] = rec.event_generators(gid);
++lidx;
}
+ ++grpidx;
}
// Generate the cell groups in parallel, with one task per cell group.
@@ -398,6 +406,15 @@ void simulation_state::remove_all_samplers() {
sassoc_handles_.clear();
}
+std::vector<probe_metadata> simulation_state::get_probe_metadata(cell_member_type probe_id) const {
+ if (auto linfo = util::value_by_key(gid_to_local_, probe_id.gid)) {
+ return cell_groups_.at(linfo->group_index)->get_probe_metadata(probe_id);
+ }
+ else {
+ return {};
+ }
+}
+
void simulation_state::set_binning_policy(binning_kind policy, time_type bin_interval) {
foreach_group(
[&](cell_group_ptr& group) { group->set_binning_policy(policy, bin_interval); });
@@ -410,9 +427,9 @@ void simulation_state::inject_events(const pse_vector& events) {
if (e.time<t_) {
throw bad_event_time(e.time, t_);
}
- // gid_to_local_ maps gid to index into local set of cells.
+ // gid_to_local_ maps gid to index in local cells and of corresponding cell group.
if (auto lidx = util::value_by_key(gid_to_local_, e.target.gid)) {
- pending_events_[*lidx].push_back(e);
+ pending_events_[lidx->cell_index].push_back(e);
}
}
}
@@ -452,6 +469,10 @@ void simulation::remove_all_samplers() {
impl_->remove_all_samplers();
}
+std::vector<probe_metadata> simulation::get_probe_metadata(cell_member_type probe_id) const {
+ return impl_->get_probe_metadata(probe_id);
+}
+
std::size_t simulation::num_spikes() const {
return impl_->num_spikes();
}
diff --git a/example/lfp/lfp.cpp b/example/lfp/lfp.cpp
index 6a0066b641d4baf7f75d8f5ced9cc583b77143cd..38217129f7f8556d17846f4441634a8424ee46a3 100644
--- a/example/lfp/lfp.cpp
+++ b/example/lfp/lfp.cpp
@@ -18,6 +18,7 @@ using arb::util::any;
using arb::util::any_cast;
using arb::util::any_ptr;
using arb::cell_gid_type;
+using arb::cell_member_type;
// Recipe represents one cable cell with one synapse, together with probes for total trans-membrane current, membrane voltage,
// ionic current density, and synaptic conductance. A sequence of spikes are presented to the one synapse on the cell.
@@ -70,12 +71,11 @@ private:
void make_cell() {
using namespace arb;
- // Set up morphology as two branches:
segment_tree tree;
- // * soma, length 20 μm radius 10 μm, with SWC tag 1.
- tree.append(arb::mnpos, {0, 0, 10, 10}, {0, 0, -10, 10}, 1);
- // * apical dendrite, length 490 μm, radius 1 μm, with SWC tag 4.
- tree.append(arb::mnpos, {0, 0, 10, 1}, {0, 0, 500, 1}, 4);
+ // Soma, length 20 μm radius 10 μm, with SWC tag 1.
+ auto soma_apex = tree.append(arb::mnpos, {0, 0, -10, 10}, {0, 0, 10, 10}, 1);
+ // Apical dendrite, length 490 μm, radius 1 μm, with SWC tag 4.
+ tree.append(soma_apex, {0, 0, 10, 1}, {0, 0, 500, 1}, 4);
cell_ = cable_cell(tree);
@@ -90,8 +90,8 @@ private:
cell_.paint(reg::tagged(1), "hh"); // (default parameters)
cell_.paint(reg::tagged(4), mechanism_desc("pas").set("e", -70));
- // Add exponential synapse at centre of soma (0.5 along branch 0).
- synapse_location_ = mlocation{0, 0.5};
+ // Add exponential synapse at centre of soma.
+ synapse_location_ = ls::on_components(0.5, reg::tagged(1));
cell_.place(synapse_location_, mechanism_desc("expsyn").set("e", 0).set("tau", 2));
}
};
@@ -99,17 +99,14 @@ private:
struct position { double x, y, z; };
struct lfp_sampler {
- lfp_sampler(const arb::place_pwlin& p, std::vector<position> electrodes, double sigma):
- placement(p), electrodes(std::move(electrodes)), sigma(sigma) {}
-
- // Compute response coefficients for each electrode, given a set of cable-like sources.
- void initialize(const arb::mcable_list& cables) {
+ lfp_sampler(const arb::place_pwlin& placement, const arb::mcable_list& cables, const std::vector<position>& electrodes, double sigma) {
+ // Compute response coefficients for each electrode, given a set of cable-like sources.
const unsigned n_electrode = electrodes.size();
response.assign(n_electrode, std::vector<double>(cables.size()));
std::vector<arb::mpoint> midpoints;
std::transform(cables.begin(), cables.end(), std::back_inserter(midpoints),
- [this](const auto& c) { return placement.at({c.branch, 0.5*(c.prox_pos+c.dist_pos)}); });
+ [&placement](const auto& c) { return placement.at({c.branch, 0.5*(c.prox_pos+c.dist_pos)}); });
const double coef = 1/(4*M_PI*sigma); // [Ω·m]
for (unsigned i = 0; i<n_electrode; ++i) {
@@ -126,27 +123,9 @@ struct lfp_sampler {
}
}
- void reset() {
- response.clear();
- lfp_time.clear();
- lfp_voltage.clear();
- }
-
- bool is_initialized() const {
- return !response.empty();
- }
-
// On receipt of a sequence of cell-wide current samples, apply response matrix and save results to lfp_voltage.
arb::sampler_function callback() {
return [this](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
- auto cables_ptr = any_cast<const arb::mcable_list*>(pm.meta);
- assert(cables_ptr);
-
- if (!is_initialized()) {
- // The first time we get metadata, build the response matrix.
- initialize(*cables_ptr);
- }
-
std::vector<double> currents;
lfp_voltage.resize(response.size());
@@ -167,9 +146,6 @@ struct lfp_sampler {
std::vector<std::vector<double>> lfp_voltage; // [mV] (one vector per electrode)
private:
- const arb::place_pwlin placement; // Represents cell morphology in space.
- const std::vector<position> electrodes; // [μm]
- const double sigma; // [S/m]
std::vector<std::vector<double>> response; // [MΩ]
};
@@ -220,11 +196,16 @@ int main(int argc, char** argv) {
{30, 0, 100}
};
- auto sample_schedule = arb::regular_schedule(sample_dt);
-
arb::morphology cell_morphology = any_cast<arb::cable_cell>(R.get_cell_description(0)).morphology();
arb::place_pwlin placed_cell(cell_morphology);
- lfp_sampler lfp(placed_cell, electrodes, 3.0);
+
+ auto probe0_metadata = sim.get_probe_metadata(cell_member_type{0, 0});
+ assert(probe0_metadata.size()==1); // Should only be one probe associated with this id.
+ arb::mcable_list current_cables = *any_cast<const arb::mcable_list*>(probe0_metadata.at(0).meta);
+
+ lfp_sampler lfp(placed_cell, current_cables, electrodes, 3.0);
+
+ auto sample_schedule = arb::regular_schedule(sample_dt);
sim.add_sampler(arb::one_probe({0, 0}), sample_schedule, lfp.callback(), arb::sampling_policy::exact);
arb::trace_vector<double, arb::mlocation> membrane_voltage;
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 1771eda6085b34217a4a19f9b3922286cf2a854e..5075872077fa91ff7ddbb80cf3571edae9b869f2 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -23,7 +23,6 @@
#include "fvm_lowered_cell.hpp"
#include "fvm_lowered_cell_impl.hpp"
#include "mech_private_field_access.hpp"
-#include "sampler_map.hpp"
#include "util/meta.hpp"
#include "util/maputil.hpp"
#include "util/rangeutil.hpp"
diff --git a/test/unit/test_kinetic_linear.cpp b/test/unit/test_kinetic_linear.cpp
index b9d5707f82cdad9eb14f07dc9fa404c70d9bb311..7ed540cc6f18616ff24215fe11d01031f0de5f96 100644
--- a/test/unit/test_kinetic_linear.cpp
+++ b/test/unit/test_kinetic_linear.cpp
@@ -13,7 +13,6 @@
#include "mech_private_field_access.hpp"
#include "fvm_lowered_cell.hpp"
#include "fvm_lowered_cell_impl.hpp"
-#include "sampler_map.hpp"
#include "simple_recipes.hpp"
#include "unit_test_catalogue.hpp"
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index fd5ddbafa14f9b92eaeab240ee1bc39711bca54f..0b42a5417618492c63d79f7ea009216c6f9ed298 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -218,9 +218,24 @@ TEST(locset, thingify) {
EXPECT_EQ(thingify(begb1, mp), (ll{{1,0}}));
EXPECT_EQ(thingify(begb2, mp), (ll{{2,0}}));
EXPECT_EQ(thingify(begb3, mp), (ll{{3,0}}));
-
- // In the absence of a spherical root, there is no branch 4.
EXPECT_THROW(thingify(begb4, mp), no_such_branch);
+
+ auto oc_beg_b0 = ls::on_components(0, reg::branch(0));
+ auto oc_mid_b1 = ls::on_components(0.5, reg::branch(1));
+ auto oc_mid_b123 = ls::on_components(0.5, join(reg::branch(1), reg::branch(2), reg::branch(3)));
+ auto oc_end_b123 = ls::on_components(1, join(reg::branch(1), reg::branch(2), reg::branch(3)));
+ auto oc_end_b123h = ls::on_components(1, join(reg::branch(1), reg::branch(2), reg::cable(3, 0, 0.5)));
+ auto oc_mid_b02 = ls::on_components(0.5, join(reg::branch(0), reg::branch(2)));
+ auto oc_end_b02 = ls::on_components(1, join(reg::branch(0), reg::branch(2)));
+
+ EXPECT_EQ(thingify(oc_beg_b0, mp), (ll{{0, 0}}));
+ EXPECT_EQ(thingify(oc_mid_b1, mp), (ll{{1, 0.5}}));
+ EXPECT_EQ(thingify(oc_mid_b123, mp), (ll{{2, 0.5}, {3, 0.25}}));
+ EXPECT_EQ(thingify(oc_end_b123, mp), (ll{{3, 1}}));
+ EXPECT_EQ(thingify(oc_end_b123h, mp), (ll{{2, 1}, {3, 0.5}}));
+ EXPECT_EQ(thingify(oc_mid_b02, mp), (ll{{0, 0.5}, {2, 0.5}}));
+ EXPECT_EQ(thingify(oc_end_b02, mp), (ll{{0, 1}, {2, 1}}));
+ EXPECT_EQ(thingify(ls::on_components(0.25, reg::cable(1, 0.5, 1)), mp), (ll{{1, 0.625}}));
}
{
auto mp = mprovider(morphology(sm));
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 0d6a48c2826a09d2a99cdb2858968925da8db4ec..7b363f2adbb6c81be5aa6fd6d1847d1b74da8723 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -13,6 +13,7 @@
#include <arbor/schedule.hpp>
#include <arbor/simple_sampler.hpp>
#include <arbor/simulation.hpp>
+#include <arbor/util/any_ptr.hpp>
#include <arbor/util/pp_util.hpp>
#include <arbor/version.hpp>
#include <arborenv/gpu_env.hpp>
@@ -1275,3 +1276,57 @@ TEST(probe, gpu_##x) { \
ARB_PP_FOREACH(RUN_GPU, PROBE_TESTS)
#endif // def ARB_GPU_ENABLED
+
+// Test simulator `get_probe_metadata` interface.
+// (No need to run this on GPU back-end as well.)
+
+TEST(probe, get_probe_metadata) {
+ // Reuse multiprobe test set-up to confirm simulator::get_probe_metadata returns
+ // correct vector of metadata.
+
+ cable_cell cell(common_morphology::m_mlt_b6);
+
+ // Paint mechanism on branches 1, 2, and 5, omitting branch 4.
+ cell.paint(reg::branch(1), mechanism_desc("param_as_state").set("p", 10.));
+ cell.paint(reg::branch(2), mechanism_desc("param_as_state").set("p", 20.));
+ cell.paint(reg::branch(5), mechanism_desc("param_as_state").set("p", 50.));
+
+ cable1d_recipe rec({cell}, false);
+ rec.catalogue() = make_unit_test_catalogue(global_default_catalogue());
+ rec.add_probe(0, 7, cable_probe_density_state{ls::terminal(), "param_as_state", "s"});
+
+ context ctx = make_context();
+ partition_hint_map phints = {
+ {cell_kind::cable, {partition_hint::max_size, partition_hint::max_size, true}}
+ };
+ simulation sim(rec, partition_load_balance(rec, ctx, phints), ctx);
+
+ std::vector<probe_metadata> mm = sim.get_probe_metadata({0, 0});
+ ASSERT_EQ(3u, mm.size());
+
+ EXPECT_EQ((cell_member_type{0, 0}), mm[0].id);
+ EXPECT_EQ((cell_member_type{0, 0}), mm[1].id);
+ EXPECT_EQ((cell_member_type{0, 0}), mm[2].id);
+
+ EXPECT_EQ(0u, mm[0].index);
+ EXPECT_EQ(1u, mm[1].index);
+ EXPECT_EQ(2u, mm[2].index);
+
+ EXPECT_EQ(7, mm[0].tag);
+ EXPECT_EQ(7, mm[1].tag);
+ EXPECT_EQ(7, mm[2].tag);
+
+ const mlocation* l0 = util::any_cast<const mlocation*>(mm[0].meta);
+ const mlocation* l1 = util::any_cast<const mlocation*>(mm[1].meta);
+ const mlocation* l2 = util::any_cast<const mlocation*>(mm[2].meta);
+
+ ASSERT_TRUE(l0);
+ ASSERT_TRUE(l1);
+ ASSERT_TRUE(l2);
+
+ std::vector<mlocation> locs = {*l0, *l1, *l2};
+ util::sort(locs);
+ EXPECT_EQ((mlocation{1, 1.}), locs[0]);
+ EXPECT_EQ((mlocation{2, 1.}), locs[1]);
+ EXPECT_EQ((mlocation{5, 1.}), locs[2]);
+}