diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index a184c0e9ad2da0bfceaee4a1287a35ed44b5e8ae..e9779a4af66a19878c352a3f3a2d829f30cfe02c 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -6,11 +6,9 @@ set(arbor_sources
backends/multicore/fvm.cpp
backends/multicore/mechanism.cpp
backends/multicore/shared_state.cpp
- backends/multicore/stimulus.cpp
communication/communicator.cpp
communication/dry_run_context.cpp
benchmark_cell_group.cpp
- builtin_mechanisms.cpp
cable_cell.cpp
cable_cell_param.cpp
cell_group_factory.cpp
@@ -69,7 +67,6 @@ if(ARB_WITH_GPU)
backends/gpu/mechanism.cpp
backends/gpu/mechanism.cu
backends/gpu/shared_state.cpp
- backends/gpu/stimulus.cpp
backends/gpu/stimulus.cu
backends/gpu/threshold_watcher.cu
backends/gpu/matrix_assemble.cu
diff --git a/arbor/backends/builtin_mech_proto.hpp b/arbor/backends/builtin_mech_proto.hpp
deleted file mode 100644
index 892f122f8165f38fe8dd7e32062c4f85e9b8904c..0000000000000000000000000000000000000000
--- a/arbor/backends/builtin_mech_proto.hpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#pragma once
-
-#include <arbor/mechanism.hpp>
-#include <arbor/mechinfo.hpp>
-
-namespace arb {
-
-// Stimulus
-
-inline const mechanism_info& builtin_stimulus_info() {
- using spec = mechanism_field_spec;
- static mechanism_info info = {
- // globals
- {},
- // parameters
- {
- {"delay", {spec::parameter, "ms", 0, 0}},
- {"duration", {spec::parameter, "ms", 0, 0}},
- {"amplitude", {spec::parameter, "nA", 0, 0}}
- },
- // state
- {},
- // ions
- {},
- // fingerprint
- "##builtin_stimulus"
- };
-
- return info;
-};
-
-template <typename B>
-concrete_mech_ptr<B> make_builtin_stimulus();
-
-} // namespace arb
diff --git a/arbor/backends/gpu/hip_api.hpp b/arbor/backends/gpu/hip_api.hpp
index 9e5e98fc6f1e684e4f9d98bd68ce34316446539f..283f68bf4a90ff16e87f99dbbbd3dc92a0850161 100644
--- a/arbor/backends/gpu/hip_api.hpp
+++ b/arbor/backends/gpu/hip_api.hpp
@@ -83,7 +83,7 @@ inline api_error_type device_mem_get_info(ARGS&&... args) {
__device__
inline double gpu_atomic_add(double* address, double val) {
-return atomicAdd(address, val);
+ return atomicAdd(address, val);
}
__device__
diff --git a/arbor/backends/gpu/math_cu.hpp b/arbor/backends/gpu/math_cu.hpp
index d7658e6d81ab3c566de1b0c3a3a03967b38e2a0b..e2692b9d9a6106276083f2503a25de433a2055e3 100644
--- a/arbor/backends/gpu/math_cu.hpp
+++ b/arbor/backends/gpu/math_cu.hpp
@@ -1,9 +1,11 @@
#pragma once
+
+// Implementations of mathematical operations required by generated CUDA mechanisms and back-end methods.
+
#include <cfloat>
+#include <cmath>
#include "gpu_api.hpp"
-// Implementations of mathematical operations required
-// by generated CUDA mechanisms.
namespace arb {
namespace gpu {
@@ -38,5 +40,13 @@ inline T max(T lhs, T rhs) {
return lhs<rhs? rhs: lhs;
}
+template <typename T>
+__device__
+inline T lerp(T a, T b, T u) {
+ return std::fma(u, b, std::fma(-u, a, a));
+}
+
+constexpr double pi = 3.1415926535897932384626433832795l;
+
} // namespace gpu
} // namespace arb
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index ee5648457a74e335610e5087f9da18854df335e8..9dfa021d89131a221792cb7429b75721574527a1 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -10,6 +10,7 @@
#include "backends/multi_event_stream_state.hpp"
#include "memory/copy.hpp"
#include "memory/wrappers.hpp"
+#include "util/index_into.hpp"
#include "util/rangeutil.hpp"
using arb::memory::make_const_view;
@@ -86,6 +87,82 @@ void ion_state::reset() {
memory::copy(init_eX_, eX_);
}
+// istim_state methods:
+
+istim_state::istim_state(const fvm_stimulus_config& stim) {
+ using util::assign;
+
+ // Translate instance-to-CV index from stim to istim_state index vectors.
+ std::vector<fvm_index_type> accu_index_stage;
+ assign(accu_index_stage, util::index_into(stim.cv, stim.cv_unique));
+
+ std::size_t n = accu_index_stage.size();
+ std::vector<fvm_value_type> envl_a, envl_t;
+ std::vector<fvm_index_type> edivs;
+
+ frequency_ = make_const_view(stim.frequency);
+
+ arb_assert(n==frequency_.size());
+ arb_assert(n==stim.envelope_time.size());
+ arb_assert(n==stim.envelope_amplitude.size());
+
+ edivs.reserve(n+1);
+ edivs.push_back(0);
+
+ for (auto i: util::make_span(n)) {
+ arb_assert(stim.envelope_time[i].size()==stim.envelope_amplitude[i].size());
+ arb_assert(util::is_sorted(stim.envelope_time[i]));
+
+ util::append(envl_a, stim.envelope_amplitude[i]);
+ util::append(envl_t, stim.envelope_time[i]);
+ edivs.push_back(fvm_index_type(envl_t.size()));
+ }
+
+ accu_index_ = make_const_view(accu_index_stage);
+ accu_to_cv_ = make_const_view(stim.cv_unique);
+ accu_stim_ = array(accu_index_.size());
+ envl_amplitudes_ = make_const_view(envl_a);
+ envl_times_ = make_const_view(envl_t);
+ envl_divs_ = make_const_view(edivs);
+
+ // Initial indices into envelope match partition divisions; ignore last (index n) element.
+ envl_index_ = envl_divs_;
+
+ // Initialize ppack pointers.
+ ppack_.accu_index = accu_index_.data();
+ ppack_.accu_to_cv = accu_to_cv_.data();
+ ppack_.frequency = frequency_.data();
+ ppack_.envl_amplitudes = envl_amplitudes_.data();
+ ppack_.envl_times = envl_times_.data();
+ ppack_.envl_divs = envl_divs_.data();
+ ppack_.accu_stim = accu_stim_.data();
+ ppack_.envl_index = envl_index_.data();
+ // The following ppack fields must be set in add_current() before queuing kernel.
+ ppack_.time = nullptr;
+ ppack_.cv_to_intdom = nullptr;
+ ppack_.current_density = nullptr;
+}
+
+std::size_t istim_state::size() const {
+ return frequency_.size();
+}
+
+void istim_state::zero_current() {
+ memory::fill(accu_stim_, 0.);
+}
+
+void istim_state::reset() {
+ zero_current();
+ memory::copy(envl_divs_, envl_index_);
+}
+
+void istim_state::add_current(const array& time, const iarray& cv_to_intdom, array& current_density) {
+ ppack_.time = time.data();
+ ppack_.cv_to_intdom = cv_to_intdom.data();
+ ppack_.current_density = current_density.data();
+ istim_add_current_impl((int)size(), ppack_);
+}
+
// Shared state methods:
shared_state::shared_state(
@@ -136,6 +213,10 @@ void shared_state::add_ion(
std::forward_as_tuple(charge, ion_info, 1u));
}
+void shared_state::configure_stimulus(const fvm_stimulus_config& stims) {
+ stim_data = istim_state(stims);
+}
+
void shared_state::reset() {
memory::copy(init_voltage, voltage);
memory::fill(current_density, 0);
@@ -147,6 +228,7 @@ void shared_state::reset() {
for (auto& i: ion_data) {
i.second.reset();
}
+ stim_data.reset();
}
void shared_state::zero_currents() {
@@ -155,6 +237,7 @@ void shared_state::zero_currents() {
for (auto& i: ion_data) {
i.second.zero_current();
}
+ stim_data.zero_current();
}
void shared_state::ions_init_concentration() {
@@ -175,6 +258,10 @@ void shared_state::add_gj_current() {
add_gj_current_impl(n_gj, gap_junctions.data(), voltage.data(), current_density.data());
}
+void shared_state::add_stimulus_current() {
+ stim_data.add_current(time, cv_to_intdom, current_density);
+}
+
std::pair<fvm_value_type, fvm_value_type> shared_state::time_bounds() const {
return minmax_value_impl(n_intdom, time.data());
}
diff --git a/arbor/backends/gpu/shared_state.cu b/arbor/backends/gpu/shared_state.cu
index 62e7c4d450c3a34c868076e9bbe3cdd52ea0d418..c01376ff506ecc5fef06986de246dfb1282e8427 100644
--- a/arbor/backends/gpu/shared_state.cu
+++ b/arbor/backends/gpu/shared_state.cu
@@ -79,7 +79,7 @@ __global__ void take_samples_impl(
auto end = s.ev_data+s.end_offset[i];
for (auto p = begin; p!=end; ++p) {
sample_time[p->offset] = time[i];
- sample_value[p->offset] = *p->handle;
+ sample_value[p->offset] = p->handle? *p->handle: 0;
}
}
}
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index 8d0254b3e95feeb9734a1bac6b3261d95d95d926..6243982e6b45428b81a20a2d53cbf60e055d4690 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -10,6 +10,7 @@
#include "fvm_layout.hpp"
#include "backends/gpu/gpu_store_types.hpp"
+#include "backends/gpu/stimulus.hpp"
namespace arb {
namespace gpu {
@@ -60,6 +61,42 @@ struct ion_state {
void reset();
};
+struct istim_state {
+ // Immutable data (post construction/initialization):
+ iarray accu_index_; // Instance to accumulator index (accu_stim_ index) map.
+ iarray accu_to_cv_; // Accumulator index to CV map.
+
+ array frequency_; // (Hz) stimulus frequency per instance.
+ array envl_amplitudes_; // (A/m²) stimulus envelope amplitudes, partitioned by instance.
+ array envl_times_; // (A/m²) stimulus envelope timepoints, partitioned by instance.
+ iarray envl_divs_; // Partition divisions for envl_ arrays,
+
+ // Mutable data:
+ array accu_stim_; // (A/m²) accumulated stim current / CV area, one per CV with a stimulus.
+ iarray envl_index_; // Per instance index into envl_ arrays, corresponding to last sample time.
+
+ // Parameter pack presents pointers to state arrays, relevant shared state to GPU kernels.
+ // Initialized at state construction.
+ istim_pp ppack_;
+
+ // Zero stim current.
+ void zero_current();
+
+ // Zero stim current, reset indices.
+ void reset();
+
+ // Contribute to current density:
+ void add_current(const array& time, const iarray& cv_to_intdom, array& current_density);
+
+ // Number of stimuli:
+ std::size_t size() const;
+
+ // Construct state from i_clamp data; references to shared state vectors are used to initialize ppack.
+ istim_state(const fvm_stimulus_config& stim_data);
+
+ istim_state() = default;
+};
+
struct shared_state {
fvm_size_type n_intdom = 0; // Number of distinct integration domains.
fvm_size_type n_detector = 0; // Max number of detectors on all cells.
@@ -84,8 +121,8 @@ struct shared_state {
array time_since_spike; // Stores time since last spike on any detector, organized by cell.
iarray src_to_spike; // Maps spike source index to spike index
+ istim_state stim_data;
std::unordered_map<std::string, ion_state> ion_data;
-
deliverable_event_stream deliverable_events;
shared_state() = default;
@@ -109,6 +146,8 @@ struct shared_state {
int charge,
const fvm_ion_config& ion_data);
+ void configure_stimulus(const fvm_stimulus_config&);
+
void zero_currents();
void ions_init_concentration();
@@ -122,6 +161,9 @@ struct shared_state {
// Update gap_junction state
void add_gj_current();
+ // Update stimulus state and add current contributions.
+ void add_stimulus_current();
+
// Return minimum and maximum time value [ms] across cells.
std::pair<fvm_value_type, fvm_value_type> time_bounds() const;
diff --git a/arbor/backends/gpu/stimulus.cpp b/arbor/backends/gpu/stimulus.cpp
deleted file mode 100644
index ba05951e9c1ade209ab4bd37fbdf65a9b00301ac..0000000000000000000000000000000000000000
--- a/arbor/backends/gpu/stimulus.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-#include "backends/builtin_mech_proto.hpp"
-#include "backends/gpu/mechanism.hpp"
-#include "backends/gpu/mechanism_ppack_base.hpp"
-
-#include "stimulus.hpp"
-
-namespace arb {
-namespace gpu {
-
-class stimulus: public arb::gpu::mechanism {
-public:
- const mechanism_fingerprint& fingerprint() const override {
- static mechanism_fingerprint hash = "##builtin_stimulus";
- return hash;
- }
- std::string internal_name() const override { return "_builtin_stimulus"; }
- mechanismKind kind() const override { return ::arb::mechanismKind::point; }
- mechanism_ptr clone() const override { return mechanism_ptr(new stimulus()); }
-
- void init() override {}
- void advance_state() override {}
- void compute_currents() override {
- stimulus_current_impl(size(), pp_);
- }
-
- void write_ions() override {}
- void apply_events(deliverable_event_stream::state events) override {}
-
- mechanism_ppack_base* ppack_ptr() override {
- return &pp_;
- }
-
-protected:
- std::size_t object_sizeof() const override { return sizeof(*this); }
-
- mechanism_field_table field_table() override {
- return {
- {"delay", &pp_.delay},
- {"duration", &pp_.duration},
- {"amplitude", &pp_.amplitude}
- };
- }
-
- mechanism_field_default_table field_default_table() override {
- return {
- {"delay", 0},
- {"duration", 0},
- {"amplitude", 0}
- };
- }
-
-private:
- stimulus_pp pp_;
-};
-
-} // namespace gpu
-
-template <>
-concrete_mech_ptr<gpu::backend> make_builtin_stimulus() {
- return concrete_mech_ptr<gpu::backend>(new arb::gpu::stimulus());
-}
-
-} // namespace arb
diff --git a/arbor/backends/gpu/stimulus.cu b/arbor/backends/gpu/stimulus.cu
index a207b376614d2874127c4349aad7dc63a1b597b5..cc3f85cad05d6f2df435d33e98a7a333ce8e8958 100644
--- a/arbor/backends/gpu/stimulus.cu
+++ b/arbor/backends/gpu/stimulus.cu
@@ -1,33 +1,62 @@
+#include <cmath>
+
#include <arbor/fvm_types.hpp>
-#include "gpu_api.hpp"
-#include "gpu_common.hpp"
-#include "stimulus.hpp"
+#include "backends/gpu/gpu_api.hpp"
+#include "backends/gpu/gpu_common.hpp"
+#include "backends/gpu/math_cu.hpp"
+#include "backends/gpu/stimulus.hpp"
namespace arb {
namespace gpu {
+// TODO: Implement version with reproducibility-friendly accumulations.
+// See Arbor issue #1059.
+
namespace kernel {
- __global__
- void stimulus_current_impl(int n, stimulus_pp pp) {
- auto i = threadIdx.x + blockDim.x*blockIdx.x;
- if (i<n) {
- auto t = pp.vec_t_[pp.vec_di_[i]];
- if (t>=pp.delay[i] && t<pp.delay[i]+pp.duration[i]) {
- // use subtraction because the electrode currents are specified
- // in terms of current into the compartment
- gpu_atomic_add(pp.vec_i_+pp.node_index_[i], -pp.weight_[i]*pp.amplitude[i]);
- }
- }
+
+__global__
+void istim_add_current_impl(int n, istim_pp pp) {
+ constexpr double freq_scale = 2*pi*0.001;
+
+ auto i = threadIdx.x + blockDim.x*blockIdx.x;
+ if (i>=n) return;
+
+ fvm_index_type ei_left = pp.envl_divs[i];
+ fvm_index_type ei_right = pp.envl_divs[i+1];
+
+ fvm_index_type ai = pp.accu_index[i];
+ fvm_index_type cv = pp.accu_to_cv[ai];
+ double t = pp.time[pp.cv_to_intdom[cv]];
+
+ if (ei_left==ei_right || t<pp.envl_times[ei_left]) return;
+
+ fvm_index_type& ei = pp.envl_index[i];
+ while (ei+1<ei_right && pp.envl_times[ei+1]<=t) ++ei;
+
+ double J = pp.envl_amplitudes[ei]; // current density (A/m²)
+ if (ei+1<ei_right) {
+ // linearly interpolate:
+ double J1 = pp.envl_amplitudes[ei+1];
+ double u = (t-pp.envl_times[ei])/(pp.envl_times[ei+1]-pp.envl_times[ei]);
+ J = lerp(J, J1, u);
+ }
+
+ if (double f = pp.frequency[i]) {
+ J *= std::sin(freq_scale*f*t);
}
-} // namespace kernel
+ gpu_atomic_add(&pp.accu_stim[ai], J);
+ gpu_atomic_sub(&pp.current_density[cv], J);
+}
+
+} // namespace kernel
-void stimulus_current_impl(int n, const stimulus_pp& pp) {
+void istim_add_current_impl(int n, const istim_pp& pp) {
constexpr unsigned block_dim = 128;
const unsigned grid_dim = impl::block_count(n, block_dim);
- kernel::stimulus_current_impl<<<grid_dim, block_dim>>>(n, pp);
+ kernel::istim_add_current_impl<<<grid_dim, block_dim>>>(n, pp);
}
} // namespace gpu
diff --git a/arbor/backends/gpu/stimulus.hpp b/arbor/backends/gpu/stimulus.hpp
index b5d41dada90543e967f14750f166949f0fa6c33f..e03edd4322568af25cb9f40ba47c7c01f830251a 100644
--- a/arbor/backends/gpu/stimulus.hpp
+++ b/arbor/backends/gpu/stimulus.hpp
@@ -1,17 +1,30 @@
#pragma once
-#include <backends/gpu/mechanism_ppack_base.hpp>
+#include <arbor/fvm_types.hpp>
namespace arb {
namespace gpu {
-struct stimulus_pp: mechanism_ppack_base {
- fvm_value_type* delay;
- fvm_value_type* duration;
- fvm_value_type* amplitude;
+// Pointer representation of state is passed to GPU kernels.
+
+struct istim_pp {
+ // Stimulus constant and mutable data:
+ const fvm_index_type* accu_index;
+ const fvm_index_type* accu_to_cv;
+ const fvm_value_type* frequency;
+ const fvm_value_type* envl_amplitudes;
+ const fvm_value_type* envl_times;
+ const fvm_index_type* envl_divs;
+ fvm_value_type* accu_stim;
+ fvm_index_type* envl_index;
+
+ // Pointers to shared state data:
+ const fvm_value_type* time;
+ const fvm_index_type* cv_to_intdom;
+ fvm_value_type* current_density;
};
-void stimulus_current_impl(int n, const stimulus_pp&);
+void istim_add_current_impl(int n, const istim_pp& pp);
} // namespace gpu
} // namespace arb
diff --git a/arbor/backends/gpu/threshold_watcher.cu b/arbor/backends/gpu/threshold_watcher.cu
index 61706f9ee4c5ea215a110a34f0e90e998a30b73f..553f85e4cf3725edbc71dc6ab36378a392805e0b 100644
--- a/arbor/backends/gpu/threshold_watcher.cu
+++ b/arbor/backends/gpu/threshold_watcher.cu
@@ -3,7 +3,7 @@
#include <arbor/fvm_types.hpp>
#include "backends/threshold_crossing.hpp"
-#include "gpu_common.hpp"
+#include "math_cu.hpp"
#include "stack_cu.hpp"
namespace arb {
@@ -11,12 +11,6 @@ namespace gpu {
namespace kernel {
-template <typename T>
-__device__
-inline T lerp(T a, T b, T u) {
- return std::fma(u, b, std::fma(-u, a, a));
-}
-
/// kernel used to test for threshold crossing test code.
/// params:
/// t : current time (ms)
@@ -67,7 +61,7 @@ void test_thresholds_impl(
// The threshold has been passed, so estimate the time using
// linear interpolation
auto pos = (thresh - v_prev)/(v - v_prev);
- crossing_time = lerp(t_before[intdom], t_after[intdom], pos);
+ crossing_time = gpu::lerp(t_before[intdom], t_after[intdom], pos);
if (record_time_since_spike) {
time_since_spike[spike_idx] = t_after[intdom] - crossing_time;
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 02f4e9f880559a49bd66317249519d8f01e757e9..0da23bcaaa7cfa4e184782fdd30bcaebce0c9859 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -16,6 +16,7 @@
#include "backends/event.hpp"
#include "io/sepval.hpp"
+#include "util/index_into.hpp"
#include "util/padded_alloc.hpp"
#include "util/rangeutil.hpp"
@@ -44,7 +45,6 @@ inline unsigned min_alignment(unsigned align) {
using pad = util::padded_allocator<>;
-
// ion_state methods:
ion_state::ion_state(
@@ -87,6 +87,95 @@ void ion_state::reset() {
std::copy(init_eX_.begin(), init_eX_.end(), eX_.begin());
}
+// istim_state methods:
+
+istim_state::istim_state(const fvm_stimulus_config& stim, unsigned align):
+ alignment(min_alignment(align)),
+ accu_to_cv_(stim.cv_unique.begin(), stim.cv_unique.end(), pad(alignment)),
+ frequency_(stim.frequency.begin(), stim.frequency.end(), pad(alignment))
+{
+ using util::assign;
+
+ // Translate instance-to-CV index from stim to istim_state index vectors.
+ assign(accu_index_, util::index_into(stim.cv, accu_to_cv_));
+ accu_stim_.resize(accu_to_cv_.size());
+
+ std::size_t n = accu_index_.size();
+ std::vector<fvm_value_type> envl_a, envl_t;
+ std::vector<fvm_index_type> edivs;
+
+ arb_assert(n==frequency_.size());
+ arb_assert(n==stim.envelope_time.size());
+ arb_assert(n==stim.envelope_amplitude.size());
+
+ edivs.reserve(n+1);
+ edivs.push_back(0);
+
+ for (auto i: util::make_span(n)) {
+ arb_assert(stim.envelope_time[i].size()==stim.envelope_amplitude[i].size());
+ arb_assert(util::is_sorted(stim.envelope_time[i]));
+
+ util::append(envl_a, stim.envelope_amplitude[i]);
+ util::append(envl_t, stim.envelope_time[i]);
+ edivs.push_back(fvm_index_type(envl_t.size()));
+ }
+
+ assign(envl_amplitudes_, envl_a);
+ assign(envl_times_, envl_t);
+ assign(envl_divs_, edivs);
+ envl_index_.assign(edivs.data(), edivs.data()+n);
+}
+
+void istim_state::zero_current() {
+ util::fill(accu_stim_, 0);
+}
+
+void istim_state::reset() {
+ zero_current();
+
+ std::size_t n = envl_index_.size();
+ std::copy(envl_divs_.data(), envl_divs_.data()+n, envl_index_.begin());
+}
+
+void istim_state::add_current(const array& time, const iarray& cv_to_intdom, array& current_density) {
+ constexpr double freq_scale = 2*math::pi<double>*0.001;
+
+ // Consider vectorizing...
+ for (auto i: util::count_along(accu_index_)) {
+ // Advance index into envelope until either
+ // - the next envelope time is greater than simulation time, or
+ // - it is the last valid index for the envelope.
+
+ fvm_index_type ei_left = envl_divs_[i];
+ fvm_index_type ei_right = envl_divs_[i+1];
+
+ fvm_index_type ai = accu_index_[i];
+ fvm_index_type cv = accu_to_cv_[ai];
+ double t = time[cv_to_intdom[cv]];
+
+ if (ei_left==ei_right || t<envl_times_[ei_left]) continue;
+
+ fvm_index_type& ei = envl_index_[i];
+ while (ei+1<ei_right && envl_times_[ei+1]<=t) ++ei;
+
+ double J = envl_amplitudes_[ei]; // current density (A/m²)
+ if (ei+1<ei_right) {
+ // linearly interpolate:
+ arb_assert(envl_times_[ei]<=t && envl_times_[ei+1]>t);
+ double J1 = envl_amplitudes_[ei+1];
+ double u = (t-envl_times_[ei])/(envl_times_[ei+1]-envl_times_[ei]);
+ J = math::lerp(J, J1, u);
+ }
+
+ if (frequency_[i]) {
+ J *= std::sin(freq_scale*frequency_[i]*t);
+ }
+
+ accu_stim_[ai] += J;
+ current_density[cv] -= J;
+ }
+}
+
// shared_state methods:
shared_state::shared_state(
@@ -155,6 +244,10 @@ void shared_state::add_ion(
std::forward_as_tuple(charge, ion_info, alignment));
}
+void shared_state::configure_stimulus(const fvm_stimulus_config& stims) {
+ stim_data = istim_state(stims, alignment);
+}
+
void shared_state::reset() {
std::copy(init_voltage.begin(), init_voltage.end(), voltage.begin());
util::fill(current_density, 0);
@@ -166,6 +259,8 @@ void shared_state::reset() {
for (auto& i: ion_data) {
i.second.reset();
}
+
+ stim_data.reset();
}
void shared_state::zero_currents() {
@@ -174,6 +269,7 @@ void shared_state::zero_currents() {
for (auto& i: ion_data) {
i.second.zero_current();
}
+ stim_data.zero_current();
}
void shared_state::ions_init_concentration() {
@@ -228,6 +324,10 @@ void shared_state::add_gj_current() {
}
}
+void shared_state::add_stimulus_current() {
+ stim_data.add_current(time, cv_to_intdom, current_density);
+}
+
std::pair<fvm_value_type, fvm_value_type> shared_state::time_bounds() const {
return util::minmax_value(time);
}
@@ -245,10 +345,11 @@ void shared_state::take_samples(
auto begin = s.begin_marked(i);
auto end = s.end_marked(i);
+ // Null handles are explicitly permitted, and always give a sample of zero.
// (Note: probably not worth explicitly vectorizing this.)
for (auto p = begin; p<end; ++p) {
sample_time[p->offset] = time[i];
- sample_value[p->offset] = *p->handle;
+ sample_value[p->offset] = p->handle? *p->handle: 0;
}
}
}
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index c36908fd32335e166318e7416c979ac3f654848d..2011c66861dd87b5fc829422421b2ea2388245d5 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -69,11 +69,41 @@ struct ion_state {
// Set ionic current density to zero.
void zero_current();
- // Zero currents, reset concentrations, and reset reversal potential from
- // initial values.
+ // Zero currents, reset concentrations, and reset reversal potential from initial values.
void reset();
};
+struct istim_state {
+ unsigned alignment = 1; // Alignment and padding multiple.
+
+ // Immutable data (post initialization):
+ iarray accu_index_; // Instance to accumulator index (accu_stim_ index) map.
+ iarray accu_to_cv_; // Accumulator index to CV map.
+
+ array frequency_; // (Hz) stimulus frequency per instance.
+ array envl_amplitudes_; // (A/m²) stimulus envelope amplitudes, partitioned by instance.
+ array envl_times_; // (A/m²) stimulus envelope timepoints, partitioned by instance.
+ iarray envl_divs_; // Partition divisions for envl_ arrays,
+
+ // Mutable data:
+ array accu_stim_; // (A/m²) accumulated stim current / CV area, one per CV with a stimulus.
+ iarray envl_index_; // Per instance index into envl_ arrays, corresponding to last sample time.
+
+ // Zero stim current.
+ void zero_current();
+
+ // Zero stim current, reset indices.
+ void reset();
+
+ // Contribute to current density:
+ void add_current(const array& time, const iarray& cv_to_intdom, array& current_density);
+
+ // Construct state from i_clamp data:
+ istim_state(const fvm_stimulus_config& stim_data, unsigned align);
+
+ istim_state() = default;
+};
+
struct shared_state {
unsigned alignment = 1; // Alignment and padding multiple.
util::padded_allocator<> alloc; // Allocator with corresponging alignment/padding.
@@ -101,8 +131,8 @@ struct shared_state {
array time_since_spike; // Stores time since last spike on any detector, organized by cell.
iarray src_to_spike; // Maps spike source index to spike index
+ istim_state stim_data;
std::unordered_map<std::string, ion_state> ion_data;
-
deliverable_event_stream deliverable_events;
shared_state() = default;
@@ -126,6 +156,8 @@ struct shared_state {
int charge,
const fvm_ion_config& ion_data);
+ void configure_stimulus(const fvm_stimulus_config&);
+
void zero_currents();
void ions_init_concentration();
@@ -141,6 +173,9 @@ struct shared_state {
// Update gap_junction state
void add_gj_current();
+ // Update stimulus state and add current contributions.
+ void add_stimulus_current();
+
// Return minimum and maximum time value [ms] across cells.
std::pair<fvm_value_type, fvm_value_type> time_bounds() const;
diff --git a/arbor/backends/multicore/stimulus.cpp b/arbor/backends/multicore/stimulus.cpp
deleted file mode 100644
index d6a8c139a816c1c97739a52e62415dc9dd66d7ca..0000000000000000000000000000000000000000
--- a/arbor/backends/multicore/stimulus.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-#include <cmath>
-
-#include <arbor/fvm_types.hpp>
-
-#include "backends/builtin_mech_proto.hpp"
-#include "backends/multicore/mechanism.hpp"
-
-namespace arb {
-
-namespace multicore {
-class stimulus: public arb::multicore::mechanism {
-public:
- const mechanism_fingerprint& fingerprint() const override {
- static mechanism_fingerprint hash = "##builtin_stimulus";
- return hash;
- }
- std::string internal_name() const override { return "_builtin_stimulus"; }
- mechanismKind kind() const override { return ::arb::mechanismKind::point; }
- mechanism_ptr clone() const override { return mechanism_ptr(new stimulus()); }
-
- void init() override {}
- void advance_state() override {}
- void compute_currents() override {
- size_type n = size();
- for (size_type i=0; i<n; ++i) {
- auto cv = node_index_[i];
- auto t = vec_t_[vec_di_[cv]];
-
- if (t>=delay[i] && t<delay[i]+duration[i]) {
- // Amplitudes are given as a current into a compartment, so subtract.
- vec_i_[cv] -= weight_[i]*amplitude[i];
- }
- }
- }
- void write_ions() override {}
- void apply_events(deliverable_event_stream::state events) override {}
-
-protected:
- std::size_t object_sizeof() const override { return sizeof(*this); }
-
- mechanism_field_table field_table() override {
- return {
- {"delay", &delay},
- {"duration", &duration},
- {"amplitude", &litude}
- };
- }
-
- mechanism_field_default_table field_default_table() override {
- return {
- {"delay", 0},
- {"duration", 0},
- {"amplitude", 0}
- };
- }
-private:
- fvm_value_type* delay;
- fvm_value_type* duration;
- fvm_value_type* amplitude;
-};
-} // namespace multicore
-
-template <>
-concrete_mech_ptr<multicore::backend> make_builtin_stimulus() {
- return concrete_mech_ptr<multicore::backend>(new arb::multicore::stimulus());
-}
-
-} // namespace arb
diff --git a/arbor/builtin_mechanisms.cpp b/arbor/builtin_mechanisms.cpp
deleted file mode 100644
index 77a288593964e54f1e6403f2b959f133e0a50a4d..0000000000000000000000000000000000000000
--- a/arbor/builtin_mechanisms.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <arbor/mechcat.hpp>
-
-#include "backends/builtin_mech_proto.hpp"
-
-#include "backends/multicore/fvm.hpp"
-#if ARB_HAVE_GPU
-#include "backends/gpu/fvm.hpp"
-#endif
-
-namespace arb {
-
-template <typename B>
-concrete_mech_ptr<B> make_builtin_stimulus();
-
-mechanism_catalogue build_builtin_mechanisms() {
- mechanism_catalogue cat;
-
- cat.add("_builtin_stimulus", builtin_stimulus_info());
-
- cat.register_implementation("_builtin_stimulus", make_builtin_stimulus<multicore::backend>());
-
-#if ARB_HAVE_GPU
- cat.register_implementation("_builtin_stimulus", make_builtin_stimulus<gpu::backend>());
-#endif
-
- return cat;
-}
-
-const mechanism_catalogue& builtin_mechanisms() {
- static mechanism_catalogue cat = build_builtin_mechanisms();
- return cat;
-}
-
-} // namespace arb
-
diff --git a/arbor/builtin_mechanisms.hpp b/arbor/builtin_mechanisms.hpp
deleted file mode 100644
index 1a6688e3821ca39f896c29cf7ddfc48d49c3470e..0000000000000000000000000000000000000000
--- a/arbor/builtin_mechanisms.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#pragma once
-
-#include <arbor/mechcat.hpp>
-
-namespace arb {
-
-const mechanism_catalogue& builtin_mechanisms();
-
-} // namespace arb
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index e97a6e4d251b9ab667028d261f688be8450d17a5..4f19daa7a6e6bbdbd6568a33ee40c87106dd5fa8 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -750,6 +750,12 @@ fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& r
}
}
+ append(left.stimuli.cv, right.stimuli.cv);
+ append(left.stimuli.cv_unique, right.stimuli.cv_unique);
+ append(left.stimuli.frequency, right.stimuli.frequency);
+ append(left.stimuli.envelope_time, right.stimuli.envelope_time);
+ append(left.stimuli.envelope_amplitude, right.stimuli.envelope_amplitude);
+
left.n_target += right.n_target;
left.post_events |= right.post_events;
@@ -1076,6 +1082,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
if (!cell.stimuli().empty()) {
const auto& stimuli = cell.stimuli();
+ fvm_stimulus_config config;
std::vector<size_type> stimuli_cv;
assign_by(stimuli_cv, stimuli, [&D, cell_idx](auto& p) {
@@ -1085,19 +1092,38 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
assign(cv_order, count_along(stimuli));
sort_by(cv_order, [&](size_type i) { return stimuli_cv[i]; });
- fvm_mechanism_config config;
- config.kind = mechanismKind::point;
- // (param_values entries must be ordered by parameter name)
- config.param_values = {{"amplitude", {}}, {"delay", {}}, {"duration", {}}};
+ std::size_t n = stimuli.size();
+ config.cv.reserve(n);
+ config.cv_unique.reserve(n);
+ config.frequency.reserve(n);
+ config.envelope_time.reserve(n);
+ config.envelope_amplitude.reserve(n);
for (auto i: cv_order) {
- config.cv.push_back(stimuli_cv[i]);
- config.param_values[0].second.push_back(stimuli[i].item.amplitude);
- config.param_values[1].second.push_back(stimuli[i].item.delay);
- config.param_values[2].second.push_back(stimuli[i].item.duration);
+ 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².
+
+ config.cv.push_back(cv);
+ config.frequency.push_back(stim.frequency);
+
+ std::size_t envl_n = stim.envelope.size();
+ std::vector<double> envl_t, envl_a;
+ envl_t.reserve(envl_n);
+ envl_a.reserve(envl_n);
+
+ for (auto [t, a]: stim.envelope) {
+ envl_t.push_back(t);
+ envl_a.push_back(a*cv_area_scale);
+ }
+ config.envelope_time.push_back(std::move(envl_t));
+ config.envelope_amplitude.push_back(std::move(envl_a));
}
- M.mechanisms["_builtin_stimulus"] = std::move(config);
+ std::unique_copy(config.cv.begin(), config.cv.end(), std::back_inserter(config.cv_unique));
+ config.cv_unique.shrink_to_fit();
+
+ M.stimuli = std::move(config);
}
// Ions:
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index cdc7c4dae5c8f84dd0d085b5e9f101a559876653..403d52d4b28eea81e2b536c70d2e92ef2a5c4146 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -265,7 +265,23 @@ struct fvm_ion_config {
// Ion-specific (initial) reversal potential per CV.
std::vector<value_type> init_revpot;
+};
+
+struct fvm_stimulus_config {
+ using value_type = fvm_value_type;
+ using index_type = fvm_index_type;
+ // CV index for each stimulus instance; monotonically increasing.
+ std::vector<index_type> cv;
+
+ // CV indices where one or more stimuli are present; strictly monotonically increasing.
+ std::vector<index_type> cv_unique;
+
+ // Frequency, amplitude info, per instance.
+ // Note that amplitudes have been scaled by 1/CV area so that they are represent as current densities, not currents.
+ std::vector<double> frequency; // [Hz]
+ std::vector<std::vector<double>> envelope_time; // [ms]
+ std::vector<std::vector<double>> envelope_amplitude; // [A/m²]
};
struct fvm_mechanism_data {
@@ -275,6 +291,9 @@ struct fvm_mechanism_data {
// Ion config, indexed by ion name.
std::unordered_map<std::string, fvm_ion_config> ions;
+ // Stimulus config.
+ fvm_stimulus_config stimuli;
+
// Total number of targets (point-mechanism points).
std::size_t n_target = 0;
diff --git a/arbor/fvm_lowered_cell.hpp b/arbor/fvm_lowered_cell.hpp
index cd480f5917e859ec661fcf5fd298e55c74524aac..3e857f8a4d98e044abd096259ca5dacbc7035c4f 100644
--- a/arbor/fvm_lowered_cell.hpp
+++ b/arbor/fvm_lowered_cell.hpp
@@ -84,9 +84,24 @@ struct fvm_probe_weighted_multi {
util::any_ptr get_metadata_ptr() const { return &metadata; }
};
+struct fvm_probe_interpolated_multi {
+ std::vector<probe_handle> raw_handles; // First half take coef[0], second half coef[1].
+ std::vector<double> coef[2];
+ mcable_list metadata;
+
+ void shrink_to_fit() {
+ raw_handles.shrink_to_fit();
+ coef[0].shrink_to_fit();
+ coef[1].shrink_to_fit();
+ metadata.shrink_to_fit();
+ }
+
+ util::any_ptr get_metadata_ptr() const { return &metadata; }
+};
+
// Trans-membrane currents require special handling!
struct fvm_probe_membrane_currents {
- std::vector<probe_handle> raw_handles; // Voltage per CV.
+ std::vector<probe_handle> raw_handles; // Voltage per CV, followed by stim current densities.
std::vector<mcable> metadata; // Cables from each CV, in CV order.
std::vector<unsigned> cv_parent; // Parent CV index for each CV.
@@ -94,6 +109,9 @@ struct fvm_probe_membrane_currents {
std::vector<double> weight; // Area of cable : area of CV.
std::vector<unsigned> cv_cables_divs; // Partitions metadata by CV index.
+ std::vector<double> stim_scale; // CV area for scaling raw stim current densities.
+ std::vector<unsigned> stim_cv; // CV index corresponding to each stim raw handle.
+
void shrink_to_fit() {
raw_handles.shrink_to_fit();
metadata.shrink_to_fit();
@@ -101,6 +119,8 @@ struct fvm_probe_membrane_currents {
cv_parent_cond.shrink_to_fit();
weight.shrink_to_fit();
cv_cables_divs.shrink_to_fit();
+ stim_scale.shrink_to_fit();
+ stim_cv.shrink_to_fit();
}
util::any_ptr get_metadata_ptr() const { return &metadata; }
@@ -120,6 +140,7 @@ struct fvm_probe_data {
fvm_probe_data(fvm_probe_interpolated p): info(std::move(p)) {}
fvm_probe_data(fvm_probe_multi p): info(std::move(p)) {}
fvm_probe_data(fvm_probe_weighted_multi p): info(std::move(p)) {}
+ fvm_probe_data(fvm_probe_interpolated_multi p): info(std::move(p)) {}
fvm_probe_data(fvm_probe_membrane_currents p): info(std::move(p)) {}
std::variant<
@@ -128,6 +149,7 @@ struct fvm_probe_data {
fvm_probe_interpolated,
fvm_probe_multi,
fvm_probe_weighted_multi,
+ fvm_probe_interpolated_multi,
fvm_probe_membrane_currents
> info = missing_probe_info{};
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index d5126cf0cb10888fd9d198befa53ece58c386339..6086383fcf2483e62d512fc5a27dfabda90aafc9 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -23,7 +23,6 @@
#include <arbor/recipe.hpp>
#include <arbor/util/any_visitor.hpp>
-#include "builtin_mechanisms.hpp"
#include "execution_context.hpp"
#include "fvm_layout.hpp"
#include "fvm_lowered_cell.hpp"
@@ -261,6 +260,15 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
state_->set_dt();
PL();
+ // Add stimulus current contributions.
+ // (Note: performed after dt, time_to calculation, in case we
+ // want to use mean current contributions as opposed to point
+ // sample.)
+
+ PE(advance_integrate_stimuli)
+ state_->add_stimulus_current();
+ PL();
+
// Take samples at cell time if sample time in this step interval.
PE(advance_integrate_samples);
@@ -428,8 +436,7 @@ void fvm_lowered_cell_impl<Backend>::initialize(
// Mechanism instantiator helper.
auto mech_instance = [&catalogue](const std::string& name) {
- auto cat = builtin_mechanisms().has(name)? &builtin_mechanisms(): catalogue;
- return cat->instance<backend>(name);
+ return catalogue->instance<backend>(name);
};
// Check for physically reasonable membrane volages?
@@ -486,7 +493,7 @@ void fvm_lowered_cell_impl<Backend>::initialize(
D.init_membrane_potential, D.temperature_K, D.diam_um, std::move(src_to_spike),
data_alignment? data_alignment: 1u);
- // Instantiate mechanisms and ions.
+ // Instantiate mechanisms, ions, and stimuli.
for (auto& i: mech_data.ions) {
const std::string& ion_name = i.first;
@@ -499,6 +506,10 @@ void fvm_lowered_cell_impl<Backend>::initialize(
}
}
+ if (!mech_data.stimuli.cv.empty()) {
+ state_->configure_stimulus(mech_data.stimuli);
+ }
+
target_handles.resize(mech_data.n_target);
// Keep track of mechanisms by name for probe lookup.
@@ -530,7 +541,6 @@ void fvm_lowered_cell_impl<Backend>::initialize(
auto cv = layout.cv[i];
layout.weight[i] = 1000/D.cv_area[cv];
- // (builtin stimulus, for example, has no targets)
if (config.target.empty()) continue;
target_handle handle(mech_id, i, cv_to_intdom[cv]);
@@ -775,6 +785,7 @@ void fvm_lowered_cell_impl<Backend>::resolve_probe_address(
cable_probe_total_ion_current_density,
cable_probe_total_ion_current_cell,
cable_probe_total_current_cell,
+ cable_probe_stimulus_current_cell,
cable_probe_density_state,
cable_probe_density_state_cell,
cable_probe_point_state,
@@ -841,28 +852,44 @@ void resolve_probe(const cable_probe_axial_current& p, probe_resolution_data<B>&
template <typename B>
void resolve_probe(const cable_probe_total_ion_current_density& p, probe_resolution_data<B>& R) {
+ // Use interpolated probe with coeffs 1, -1 to represent difference between accumulated current density and stimulus.
for (mlocation loc: thingify(p.locations, R.cell.provider())) {
- R.result.push_back(fvm_probe_scalar{
- {R.state->current_density.data() + R.D.geometry.location_cv(R.cell_idx, loc, cv_prefer::cv_nonempty)},
+ fvm_index_type cv = R.D.geometry.location_cv(R.cell_idx, loc, cv_prefer::cv_nonempty);
+ const double* current_cv_ptr = R.state->current_density.data() + cv;
+
+ auto opt_i = util::binary_search_index(R.M.stimuli.cv_unique, cv);
+ const double* stim_cv_ptr = opt_i? R.state->stim_data.accu_stim_.data()+*opt_i: nullptr;
+
+ R.result.push_back(fvm_probe_interpolated{
+ {current_cv_ptr, stim_cv_ptr},
+ {1., -1.},
loc});
}
}
template <typename B>
void resolve_probe(const cable_probe_total_ion_current_cell& p, probe_resolution_data<B>& R) {
- fvm_probe_weighted_multi r;
+ fvm_probe_interpolated_multi r;
+ std::vector<const double*> stim_handles;
for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
- const double* ptr = R.state->current_density.data()+cv;
+ const double* current_cv_ptr = R.state->current_density.data()+cv;
+ auto opt_i = util::binary_search_index(R.M.stimuli.cv_unique, cv);
+ const double* stim_cv_ptr = opt_i? R.state->stim_data.accu_stim_.data()+*opt_i: nullptr;
+
for (auto cable: R.D.geometry.cables(cv)) {
double area = R.cell.embedding().integrate_area(cable); // [µm²]
if (area>0) {
- r.raw_handles.push_back(ptr);
- r.weight.push_back(0.001*area); // Scale from [µm²·A/m²] to [nA].
+ r.raw_handles.push_back(current_cv_ptr);
+ stim_handles.push_back(stim_cv_ptr);
+ r.coef[0].push_back(0.001*area); // Scale from [µm²·A/m²] to [nA].
+ r.coef[1].push_back(-r.coef[0].back());
r.metadata.push_back(cable);
}
}
}
+
+ util::append(r.raw_handles, stim_handles);
r.shrink_to_fit();
R.result.push_back(std::move(r));
}
@@ -878,6 +905,9 @@ void resolve_probe(const cable_probe_total_current_cell& p, probe_resolution_dat
[cv0](auto cv) { return cv+1==0? cv: cv-cv0; }));
util::assign(r.cv_parent_cond, util::subrange_view(R.D.face_conductance, cell_cv_ival));
+ const auto& stim_cvs = R.M.stimuli.cv_unique;
+ const fvm_value_type* stim_src = R.state->stim_data.accu_stim_.data();
+
r.cv_cables_divs = {0};
for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
r.raw_handles.push_back(R.state->voltage.data()+cv);
@@ -892,6 +922,39 @@ void resolve_probe(const cable_probe_total_current_cell& p, probe_resolution_dat
}
r.cv_cables_divs.push_back(r.metadata.size());
}
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ auto opt_i = util::binary_search_index(stim_cvs, cv);
+ if (!opt_i) continue;
+
+ r.raw_handles.push_back(stim_src+*opt_i);
+ r.stim_cv.push_back(cv-cv0);
+ r.stim_scale.push_back(0.001*R.D.cv_area[cv]); // Scale from [µm²·A/m²] to [nA].
+ }
+ r.shrink_to_fit();
+ R.result.push_back(std::move(r));
+}
+
+template <typename B>
+void resolve_probe(const cable_probe_stimulus_current_cell& p, probe_resolution_data<B>& R) {
+ fvm_probe_weighted_multi r;
+
+ const auto& stim_cvs = R.M.stimuli.cv_unique;
+ const fvm_value_type* src = R.state->stim_data.accu_stim_.data();
+
+ for (auto cv: R.D.geometry.cell_cvs(R.cell_idx)) {
+ auto opt_i = util::binary_search_index(stim_cvs, cv);
+ const double* ptr = opt_i? src+*opt_i: nullptr;
+
+ for (auto cable: R.D.geometry.cables(cv)) {
+ double area = R.cell.embedding().integrate_area(cable); // [µm²]
+ if (area>0) {
+ r.raw_handles.push_back(ptr);
+ r.weight.push_back(0.001*area); // Scale from [µm²·A/m²] to [nA].
+ r.metadata.push_back(cable);
+ }
+ }
+ }
+
r.shrink_to_fit();
R.result.push_back(std::move(r));
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 46468c903e4c0a86c4dc71777263b98ce83fee95..3185cf0f240f86cba756432da3af1b6e68555273 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -84,7 +84,7 @@ struct cable_probe_axial_current {
locset locations;
};
-// Total current density [A/m²] across membrane _excluding_ capacitive current at `location`.
+// Total current density [A/m²] across membrane _excluding_ capacitive and stimulus current at `location`.
// Sample value type: `cable_sample_range`
// Sample metadata type: `mlocation`
struct cable_probe_total_ion_current_density {
@@ -96,11 +96,16 @@ struct cable_probe_total_ion_current_density {
// Sample metadata type: `mcable_list`
struct cable_probe_total_ion_current_cell {};
-// Total membrane current [nA] across components of the cell.
+// Total membrane current [nA] across components of the cell _excluding_ stimulus currents.
// Sample value type: `cable_sample_range`
// Sample metadata type: `mcable_list`
struct cable_probe_total_current_cell {};
+// Stimulus currents [nA] across components of the cell.
+// Sample value type: `cable_sample_range`
+// Sample metadata type: `mcable_list`
+struct cable_probe_stimulus_current_cell {};
+
// Value of state variable `state` in density mechanism `mechanism` in CV at `location`.
// Sample value type: `double`
// Sample metadata type: `mlocation`
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index edd9a0427418d3d3049068d522473559465c7c39..e7f5de4cf1cd3a69e83f430d5e9b8b67232eda0b 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -33,19 +33,48 @@ struct cable_cell_ion_data {
std::optional<double> init_reversal_potential;
};
-// Current clamp description for stimulus specification.
+// Clamp current is described by a sine wave with amplitude governed by a
+// piecewise linear envelope. A frequency of zero indicates that the current is
+// simply that given by the envelope.
+//
+// The envelope is given by a series of envelope_point values:
+// * The time points must be monotonically increasing.
+// * Onset and initial amplitude is given by the first point.
+// * The amplitude for time after the last time point is that of the last
+// amplitgude point; an explicit zero amplitude point must be provided if the
+// envelope is intended to have finite support.
+//
+// Periodic envelopes are not supported, but may well be a feature worth
+// considering in the future.
+
struct i_clamp {
- using value_type = double;
+ struct envelope_point {
+ double t; // [ms]
+ double amplitude; // [nA]
+ };
- value_type delay = 0; // [ms]
- value_type duration = 0; // [ms]
- value_type amplitude = 0; // [nA]
+ std::vector<envelope_point> envelope;
+ double frequency = 0; // [Hz] 0 => constant
+ // A default constructed i_clamp, with empty envelope, describes
+ // a trivial stimulus, providing no current at all.
i_clamp() = default;
- i_clamp(value_type delay, value_type duration, value_type amplitude):
- delay(delay), duration(duration), amplitude(amplitude)
+ // The simple constructor describes a constant amplitude stimulus starting from t=0.
+ explicit i_clamp(double amplitude, double frequency = 0):
+ envelope({{0., amplitude}}),
+ frequency(frequency)
{}
+
+ // Describe a stimulus by envelope and frequency.
+ explicit i_clamp(std::vector<envelope_point> envelope, double frequency = 0):
+ envelope(std::move(envelope)),
+ frequency(frequency)
+ {}
+
+ // A 'box' stimulus with fixed onset time, duration, and constant amplitude.
+ i_clamp(double onset, double duration, double amplitude, double frequency = 0):
+ i_clamp({{onset, amplitude}, {onset+duration, amplitude}, {onset+duration, 0.}}, frequency) {}
};
// Threshold detector description.
diff --git a/arbor/include/arbor/mechcat.hpp b/arbor/include/arbor/mechcat.hpp
index c10b6f4bd96fc1ada8ad2830e35cf88ad323c309..876e548243ec4f88120ec59eb0c501d7fc7d4cda 100644
--- a/arbor/include/arbor/mechcat.hpp
+++ b/arbor/include/arbor/mechcat.hpp
@@ -27,8 +27,7 @@
// after any modification to the catalogue.
//
// There is in addition a global default mechanism catalogue object that is
-// populated with any builtin mechanisms and mechanisms generated from
-// module files included with arbor.
+// populated with any mechanisms generated from module files included with arbor.
//
// When a mechanism name of the form "mech/param=value,..." is requested, if the
// mechanism of that name does not already exist in the catalogue, it will be
diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp
index ba77e032e51673045b3bb972d00aee2e5979b62e..2a013afe696e8f1668900a84a8bce36631f365a5 100644
--- a/arbor/mc_cell_group.cpp
+++ b/arbor/mc_cell_group.cpp
@@ -241,6 +241,53 @@ void run_samples(
sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
}
+void run_samples(
+ const fvm_probe_interpolated_multi& p,
+ const sampler_call_info& sc,
+ const fvm_value_type* raw_times,
+ const fvm_value_type* raw_samples,
+ std::vector<sample_record>& sample_records,
+ fvm_probe_scratch& scratch)
+{
+ const sample_size_type n_raw_per_sample = p.raw_handles.size();
+ const sample_size_type n_interp_per_sample = n_raw_per_sample/2;
+ sample_size_type n_sample = (sc.end_offset-sc.begin_offset)/n_raw_per_sample;
+ arb_assert((sc.end_offset-sc.begin_offset)==n_sample*n_raw_per_sample);
+ arb_assert((unsigned)n_interp_per_sample==p.coef[0].size());
+ arb_assert((unsigned)n_interp_per_sample==p.coef[1].size());
+
+ auto& sample_ranges = std::get<std::vector<cable_sample_range>>(scratch);
+ sample_ranges.clear();
+ sample_records.clear();
+
+ auto& tmp = std::get<std::vector<double>>(scratch);
+ tmp.clear();
+ tmp.reserve(n_interp_per_sample*n_sample);
+
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_raw_per_sample+sc.begin_offset;
+ const auto* raw_a = raw_samples + offset;
+ const auto* raw_b = raw_a + n_interp_per_sample;
+ for (sample_size_type i = 0; i<n_interp_per_sample; ++i) {
+ tmp.push_back(raw_a[i]*p.coef[0][i]+raw_b[i]*p.coef[1][i]);
+ }
+ }
+
+ const double* tmp_ptr = tmp.data();
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ sample_ranges.push_back({tmp_ptr, tmp_ptr+n_interp_per_sample});
+ tmp_ptr += n_interp_per_sample;
+ }
+
+ const auto& csample_ranges = sample_ranges;
+ for (sample_size_type j = 0; j<n_sample; ++j) {
+ auto offset = j*n_interp_per_sample+sc.begin_offset;
+ sample_records.push_back(sample_record{time_type(raw_times[offset]), &csample_ranges[j]});
+ }
+
+ sc.sampler({sc.probe_id, sc.tag, sc.index, p.get_metadata_ptr()}, n_sample, sample_records.data());
+}
+
void run_samples(
const fvm_probe_membrane_currents& p,
const sampler_call_info& sc,
@@ -256,6 +303,8 @@ void run_samples(
const auto n_cable = p.metadata.size();
const auto n_cv = p.cv_parent_cond.size();
const auto cables_by_cv = util::partition_view(p.cv_cables_divs);
+ const auto n_stim = p.stim_scale.size();
+ arb_assert(n_stim+n_cv==(unsigned)n_raw_per_sample);
auto& sample_ranges = std::get<std::vector<cable_sample_range>>(scratch);
sample_ranges.clear();
@@ -291,6 +340,16 @@ void run_samples(
}
}
+ const double* stim = raw_samples+offset+n_cv;
+ for (auto i: util::make_span(n_stim)) {
+ double cv_stim_I = stim[i]*p.stim_scale[i];
+ unsigned cv = p.stim_cv[i];
+ arb_assert(cv<n_cv);
+
+ for (auto cable_i: util::make_span(cables_by_cv[cv])) {
+ tmp_base[cable_i] -= cv_stim_I*p.weight[cable_i];
+ }
+ }
sample_ranges.push_back({tmp_base, tmp_base+n_cable});
}
diff --git a/doc/concepts/decor.rst b/doc/concepts/decor.rst
index 2fd773000fe7fc98d7c3e2ca3ee9e91e350f7de1..2177775dca3e2aea59e0b2693879584769576bcc 100644
--- a/doc/concepts/decor.rst
+++ b/doc/concepts/decor.rst
@@ -287,6 +287,37 @@ See :ref:`modelgapjunctions`.
4. Stimuli
~~~~~~~~~~
+A current stimulus is a DC or sinusoidal current of fixed frequemcy with a time-varying amplitude
+governed by a piecewise-linear envelope.
+
+The stimulus is described by two parameters: a frequency in Hertz, where a value of zero denotes DC;
+and a sequence of points (*t*\ :sub:`i`\ , *a*\ :sub:`i`\ ) describing the envelope, where the times
+*t*\ :sub:`i` are in milliseconds and the amplitudes *a*\ :sub:`i` are in nanoamperes.
+
+The stimulus starts at the first timepoint *t*\ :sub:`0` with amplitude *a*\ :sub:`0`, and the amplitude
+is then interpolated linearly between successive points. The last envelope point
+(*t*\ :sub:`n`\ , *a*\ :sub:`n`\ ) describes a constant amplitude *a*\ :sub:`n` from
+the time *t*\ :sub:`n` onwards.
+
+Stimulus objects in the C++ and Python interfaces have simple constructors for describing
+constant stimuli and constant amplitude stimuli restricted to a fixed time interval.
+
+.. code-block:: Python
+
+ # Constant stimulus, amplitude 10 nA.
+ decor.place('(root)', arbor.iclamp(10))
+
+ # Constant amplitude 10 nA stimulus at 20 Hz.
+ decor.place('(root)', arbor.iclamp(10, 20))
+
+ # Stimulus at 20 Hz, amplitude 10 nA, for 40 ms starting at t = 30 ms.
+ decor.place('(root)', arbor.iclamp(30, 40, 10, 20))
+
+ # Piecewise linear stimulus with amplitude ranging from 0 nA to 10 nA,
+ # starting at t = 30 ms and stopping at t = 50 ms.
+ decor.place('(root)', arbor.iclamp([(30, 0), (37, 10), (43, 8), (50, 0)])
+
+
.. _cablecell-probes:
5. Probes
diff --git a/doc/cpp/probe_sample.rst b/doc/cpp/probe_sample.rst
index 7ca2d0b822e35bb54d0cdcf54755801dfa23180f..76e16bfac2599945a6f253363eb8fa86b36eb810 100644
--- a/doc/cpp/probe_sample.rst
+++ b/doc/cpp/probe_sample.rst
@@ -159,7 +159,7 @@ Membrane current density at given locations _excluding_ capacitive currents.
struct cable_probe_total_ion_current_cell {};
-Membrane current _excluding_ capacitive currents across components of the cell.
+Membrane current _excluding_ capacitive currents and stimuli across components of the cell.
* Sample value: ``cable_sample_range``. Each value is the current in
nanoamperes across an unbranched component of the cell, as determined
@@ -173,7 +173,7 @@ Membrane current _excluding_ capacitive currents across components of the cell.
struct cable_probe_total_current_cell {};
-Total membrance current across components of the cell.
+Total membrance current excluding current stimuli across components of the cell.
* Sample value: ``cable_sample_range``. Each value is the current in
nanoamperes across an unbranched component of the cell, as determined
@@ -182,6 +182,19 @@ Total membrance current across components of the cell.
* Metadata: ``mcable_list``. Each cable in the cable list describes
the unbranched component for the corresponding sample value.
+.. code::
+
+ struct cable_probe_stimulus_current_cell {};
+
+Total stimulus currents applied across components of the cell.
+
+* Sample value: ``cable_sample_range``. Each value is the current in
+ nanoamperes across an unbranched component of the cell, as determined
+ by the discretisation. Components of CVs where no stimulus is present
+ will report a corresponding stimulus value of zero.
+
+* Metadata: ``mcable_list``. Each cable in the cable list describes
+ the unbranched component for the corresponding sample value.
Ion concentration
^^^^^^^^^^^^^^^^^
diff --git a/doc/python/probe_sample.rst b/doc/python/probe_sample.rst
index 663682660ee266e060daf76b60a3550e683c6009..a0459d3e33603f112afc74534a0f95e546ba9dbb 100644
--- a/doc/python/probe_sample.rst
+++ b/doc/python/probe_sample.rst
@@ -69,7 +69,7 @@ Total ionic current
.. py:function:: cable_probe_total_ion_current_cell()
Transmembrane current (nA) _excluding_ capacitive currents across each
- cable in each CV of the cell discretization.
+ cable in each CV of the cell discretization. Stimulus currents are not included.
Metadata: the list of corresponding :class:`cable` objects.
@@ -77,7 +77,14 @@ Total transmembrane current
.. py:function:: cable_probe_total_current_cell()
Transmembrane current (nA) _including_ capacitive currents across each
- cable in each CV of the cell discretization.
+ cable in each CV of the cell discretization. Stimulus currents are not included.
+
+ Metadata: the list of corresponding :class:`cable` objects.
+
+Total stimulus current
+ .. py:function:: cable_probe_stimulus_current_cell()
+
+ Total stimulus current (nA) across each cable in each CV of the cell discretization.
Metadata: the list of corresponding :class:`cable` objects.
diff --git a/python/cable_probes.cpp b/python/cable_probes.cpp
index fb7e1be401523bb6dc95a08a2e5c7188a5aa0683..a8543a07e4ca570bc02792b432a0ebc88ab340b6 100644
--- a/python/cable_probes.cpp
+++ b/python/cable_probes.cpp
@@ -157,6 +157,10 @@ arb::probe_info cable_probe_total_current_cell() {
return arb::cable_probe_total_current_cell{};
}
+arb::probe_info cable_probe_stimulus_current_cell() {
+ return arb::cable_probe_stimulus_current_cell{};
+}
+
arb::probe_info cable_probe_density_state(const char* where, const char* mechanism, const char* state) {
return arb::cable_probe_density_state{arb::locset(where), mechanism, state};
};
@@ -243,12 +247,15 @@ void register_cable_probes(pybind11::module& m, pyarb_global_ptr global_ptr) {
m.def("cable_probe_total_current_cell", &cable_probe_total_current_cell,
"Probe specification for cable cell total transmembrane current for each cable in each CV.");
+ m.def("cable_probe_stimulus_current_cell", &cable_probe_stimulus_current_cell,
+ "Probe specification for cable cell stimulus current across each cable in each CV.");
+
m.def("cable_probe_density_state", &cable_probe_density_state,
"Probe specification for a cable cell density mechanism state variable at points in a location set.",
"where"_a, "mechanism"_a, "state"_a);
m.def("cable_probe_density_state_cell", &cable_probe_density_state_cell,
- "Probe specification for a cable cell density mechanism state variable on each cable in each CV.",
+ "Probe specification for a cable cell density mechanism state variable on each cable in each CV where defined.",
"mechanism"_a, "state"_a);
m.def("cable_probe_point_state", &cable_probe_point_state,
diff --git a/python/cells.cpp b/python/cells.cpp
index 5788bb86b437befab6ecfd356aa051965d64e922..76afd9e5fe18fc2969968ccd58c7f9e72ef27fa8 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -285,20 +285,42 @@ void register_cells(pybind11::module& m) {
// arb::i_clamp
pybind11::class_<arb::i_clamp> i_clamp(m, "iclamp",
- "A current clamp, for injecting a single pulse of current with fixed duration and current.");
+ "A current clamp for injecting a DC or fixed frequency current governed by a piecewise linear envelope.");
i_clamp
.def(pybind11::init(
- [](double ts, double dur, double cur) {
- return arb::i_clamp{ts, dur, cur};
- }), "tstart"_a=0, "duration"_a=0, "current"_a=0)
- .def_readonly("tstart", &arb::i_clamp::delay, "Time at which current starts [ms]")
- .def_readonly("duration", &arb::i_clamp::duration, "Duration of the current injection [ms]")
- .def_readonly("current", &arb::i_clamp::amplitude, "Amplitude of the injected current [nA]")
- .def("__repr__", [](const arb::i_clamp& c){
- return util::pprintf("(iclamp {} {} {})", c.delay, c.duration, c.amplitude);})
- .def("__str__", [](const arb::i_clamp& c){
- return util::pprintf("<arbor.iclamp: tstart {} ms, duration {} ms, current {} nA>",
- c.delay, c.duration, c.amplitude);});
+ [](double ts, double dur, double cur, double frequency) {
+ return arb::i_clamp{ts, dur, cur, frequency};
+ }), "tstart"_a, "duration"_a, "current"_a, "frequency"_a=0,
+ "Construct finite duration current clamp, constant amplitude")
+ .def(pybind11::init(
+ [](double cur, double frequency) {
+ return arb::i_clamp{cur, frequency};
+ }), "current"_a, "frequency"_a=0,
+ "Construct constant amplitude current clamp")
+ .def(pybind11::init(
+ [](std::vector<std::pair<double, double>> envl, double frequency) {
+ arb::i_clamp clamp;
+ for (const auto& p: envl) {
+ clamp.envelope.push_back({p.first, p.second});
+ }
+ clamp.frequency = frequency;
+ return clamp;
+ }), "envelope"_a, "frequency"_a=0,
+ "Construct current clamp according to (time, amplitude) linear envelope")
+ .def_property_readonly("envelope",
+ [](const arb::i_clamp& obj) {
+ std::vector<std::pair<double, double>> envl;
+ for (const auto& p: obj.envelope) {
+ envl.push_back({p.t, p.amplitude});
+ }
+ return envl;
+ },
+ "List of (time [ms], amplitude [nA]) points comprising the piecewise linear envelope")
+ .def_readonly("frequency", &arb::i_clamp::frequency, "Oscillation frequency [Hz], zero implies DC stimulus.")
+ .def("__repr__", [](const arb::i_clamp& c) {
+ return util::pprintf("<arbor.iclamp: frequency {} Hz>", c.frequency);})
+ .def("__str__", [](const arb::i_clamp& c) {
+ return util::pprintf("<arbor.iclamp: frequency {} Hz>", c.frequency);});
// arb::threshold_detector
diff --git a/python/test/unit/test_cable_probes.py b/python/test/unit/test_cable_probes.py
index 7d8c663dda8e4cdc2c74afa45b0ba9f31da1562b..8a39d3cfff1ad3a45224e39663b424b5dc2c01db 100644
--- a/python/test/unit/test_cable_probes.py
+++ b/python/test/unit/test_cable_probes.py
@@ -29,6 +29,7 @@ class cc_recipe(A.recipe):
dec.place('(location 0 0.08)', "expsyn")
dec.place('(location 0 0.09)', "exp2syn")
+ dec.place('(location 0 0.1)', A.iclamp(20.))
dec.paint('(all)', "hh")
self.cell = A.cable_cell(st, A.label_dict(), dec)
@@ -89,6 +90,8 @@ class cc_recipe(A.recipe):
A.cable_probe_ion_ext_concentration(where='(location 0 0.14)', ion='na'),
# probe id (0, 15)
A.cable_probe_ion_ext_concentration_cell(ion='na'),
+ # probe id (0, 15)
+ A.cable_probe_stimulus_current_cell()
]
def cell_description(self, gid):
@@ -172,6 +175,10 @@ class CableProbes(unittest.TestCase):
self.assertEqual(1, len(m))
self.assertEqual(all_cv_cables, m[0])
+ m = sim.probe_metadata((0, 16))
+ self.assertEqual(1, len(m))
+ self.assertEqual(all_cv_cables, m[0])
+
def suite():
# specify class and test functions in tuple (here: all tests starting with 'test' from class Contexts
suite = unittest.makeSuite(CableProbes, ('test'))
diff --git a/python/test/unit/test_decor.py b/python/test/unit/test_decor.py
new file mode 100644
index 0000000000000000000000000000000000000000..81fe7e497d530f8138fa70ef23e3ea64a1ef9947
--- /dev/null
+++ b/python/test/unit/test_decor.py
@@ -0,0 +1,61 @@
+# -*- coding: utf-8 -*-
+
+import unittest
+import arbor as A
+
+# to be able to run .py file from child directory
+import sys, os
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../../')))
+
+try:
+ import options
+except ModuleNotFoundError:
+ from test import options
+
+"""
+Tests for decor and decoration wrappers.
+TODO: Coverage for more than just iclamp.
+"""
+
+class DecorClasses(unittest.TestCase):
+ def test_iclamp(self):
+ # Constant amplitude iclamp:
+ clamp = A.iclamp(10);
+ self.assertEqual(0, clamp.frequency)
+ self.assertEqual([(0, 10)], clamp.envelope)
+
+ clamp = A.iclamp(10, 20);
+ self.assertEqual(20, clamp.frequency)
+ self.assertEqual([(0, 10)], clamp.envelope)
+
+ # Square pulse:
+ clamp = A.iclamp(100, 20, 3);
+ self.assertEqual(0, clamp.frequency)
+ self.assertEqual([(100, 3), (120, 3), (120, 0)], clamp.envelope)
+
+ clamp = A.iclamp(100, 20, 3, 7);
+ self.assertEqual(7, clamp.frequency)
+ self.assertEqual([(100, 3), (120, 3), (120, 0)], clamp.envelope)
+
+ # Explicit envelope:
+ envelope = [(1, 10), (3, 30), (5, 50), (7, 0)]
+ clamp = A.iclamp(envelope);
+ self.assertEqual(0, clamp.frequency)
+ self.assertEqual(envelope, clamp.envelope)
+
+ clamp = A.iclamp(envelope, 7);
+ self.assertEqual(7, clamp.frequency)
+ self.assertEqual(envelope, clamp.envelope)
+
+def suite():
+ # specify class and test functions in tuple (here: all tests starting with 'test' from class Contexts
+ suite = unittest.makeSuite(DecorClasses, ('test'))
+ return suite
+
+def run():
+ v = options.parse_arguments().verbosity
+ runner = unittest.TextTestRunner(verbosity = v)
+ runner.run(suite())
+
+if __name__ == "__main__":
+ run()
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index e4a5635c38e8c2fcc42e3896efa43f4984238df1..5d11cff324a497ea27a99cfd47a81350eda2e1c1 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -1,3 +1,4 @@
+#include <cmath>
#include <string>
#include <vector>
@@ -353,9 +354,6 @@ TEST(fvm_lowered, stimulus) {
const fvm_size_type soma_cv = 0u;
const fvm_size_type tip_cv = 5u;
- // now we have two stims :
- //
- //
// The implementation of the stimulus is tested by creating a lowered cell, then
// testing that the correct currents are injected at the correct control volumes
// as during the stimulus windows.
@@ -373,10 +371,6 @@ TEST(fvm_lowered, stimulus) {
fvm_cell fvcell(context);
fvcell.initialize({0}, cable1d_recipe(cells), cell_to_intdom, targets, probe_map);
- mechanism* stim = find_mechanism(fvcell, "_builtin_stimulus");
- ASSERT_TRUE(stim);
- EXPECT_EQ(2u, stim->size());
-
auto& state = *(fvcell.*private_state_ptr).get();
auto& J = state.current_density;
auto& T = state.time;
@@ -384,30 +378,86 @@ TEST(fvm_lowered, stimulus) {
// Test that no current is injected at t=0.
memory::fill(J, 0.);
memory::fill(T, 0.);
- stim->update_current();
+ state.add_stimulus_current();
for (auto j: J) {
EXPECT_EQ(j, 0.);
}
+ constexpr double reltol = 1e-10;
+ using testing::near_relative;
+
// Test that 0.1 nA current is injected at soma at t=1.
memory::fill(J, 0.);
memory::fill(T, 1.);
- stim->update_current();
+ state.add_stimulus_current();
constexpr double unit_factor = 1e-3; // scale A/m²·µm² to nA
- EXPECT_DOUBLE_EQ(-0.1, J[soma_cv]*A[soma_cv]*unit_factor);
+ EXPECT_TRUE(near_relative(-0.1, J[soma_cv]*A[soma_cv]*unit_factor, reltol));
// Test that 0.1 nA is again injected at t=1.5, for a total of 0.2 nA.
memory::fill(T, 1.);
- stim->update_current();
- EXPECT_DOUBLE_EQ(-0.2, J[soma_cv]*A[soma_cv]*unit_factor);
+ state.add_stimulus_current();
+ EXPECT_TRUE(near_relative(-0.2, J[soma_cv]*A[soma_cv]*unit_factor, reltol));
// Test that at t=10, no more current is injected at soma, and that
// that 0.3 nA is injected at dendrite tip.
memory::fill(T, 10.);
- stim->update_current();
- EXPECT_DOUBLE_EQ(-0.2, J[soma_cv]*A[soma_cv]*unit_factor);
- EXPECT_DOUBLE_EQ(-0.3, J[tip_cv]*A[tip_cv]*unit_factor);
+ state.add_stimulus_current();
+ EXPECT_TRUE(near_relative(-0.2, J[soma_cv]*A[soma_cv]*unit_factor, reltol));
+ EXPECT_TRUE(near_relative(-0.3, J[tip_cv]*A[tip_cv]*unit_factor, reltol));
+}
+
+TEST(fvm_lowered, ac_stimulus) {
+ // Simple cell (one CV) with oscillating stimulus.
+
+ arb::execution_context context; // Just use default context for this one!
+
+ decor dec;
+ segment_tree tree;
+ tree.append(mnpos, {0., 0., 0., 1.}, {100., 0., 0., 1.}, 1);
+
+ const double freq = 20; // (Hz)
+ const double max_amplitude = 30; // (nA)
+ const double max_time = 8; // (ms)
+
+ // Envelope is linear ramp from 0 to max_time.
+ dec.place(mlocation{0, 0}, i_clamp({{0, 0}, {max_time, max_amplitude}, {max_time, 0}}, freq));
+ std::vector<cable_cell> cells = {cable_cell(tree, {}, dec)};
+
+ cable_cell_global_properties gprop;
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters, context);
+ const auto& A = D.cv_area;
+
+ std::vector<target_handle> targets;
+ probe_association_map probe_map;
+ std::vector<fvm_index_type> cell_to_intdom(cells.size(), 0);
+
+ fvm_cell fvcell(context);
+ fvcell.initialize({0}, cable1d_recipe(cells), cell_to_intdom, targets, probe_map);
+
+ auto& state = *(fvcell.*private_state_ptr).get();
+ auto& J = state.current_density;
+ auto& T = state.time;
+
+ // Current from t = 0 to max_time should be given by
+ // I = max_amplitude * t/max_time * sin(2Ï€ f t)
+ // where t is in ms and f = freq/1000 is the frequency in kHz.
+
+ constexpr double reltol = 1e-10;
+ using testing::near_relative;
+
+ for (double t: {0., 0.1*max_time, 0.7*max_time, 1.1*max_time}) {
+ constexpr double unit_factor = 1e-3; // scale A/m²·µm² to nA
+
+ memory::fill(J, 0.);
+ memory::fill(T, t);
+ state.add_stimulus_current();
+
+ double expected_I = t<=max_time? max_amplitude*t/max_time*std::sin(2*math::pi<double>*t*0.001*freq): 0;
+ EXPECT_TRUE(near_relative(-expected_I, J[0]*A[0]*unit_factor, reltol));
+ }
}
// Test derived mechanism behaviour.
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 867821992a765ad4b050f4dd2e4604096aee66a5..f08a4cf0f34f3a6a82f421a1b2fc07e4944f6926 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -1,6 +1,7 @@
#include "../gtest.h"
#include <cmath>
+#include <map>
#include <vector>
#include <arbor/cable_cell.hpp>
@@ -142,17 +143,19 @@ void run_v_i_probe_test(const context& ctx) {
// Voltage probes are interpolated, so expect fvm_probe_info
// to wrap an fvm_probe_interpolated; ion current density is
- // a scalar, so should wrap fvm_probe_scalar.
+ // implemented as an interpolated probe too, in order to account
+ // for stimulus contributions.
ASSERT_TRUE(std::get_if<fvm_probe_interpolated>(&probe_map.data_on({0, 0}).front().info));
- ASSERT_TRUE(std::get_if<fvm_probe_interpolated>(&probe_map.data_on({0, 0}).front().info));
- ASSERT_TRUE(std::get_if<fvm_probe_scalar>(&probe_map.data_on({0, 2}).front().info));
+ ASSERT_TRUE(std::get_if<fvm_probe_interpolated>(&probe_map.data_on({0, 1}).front().info));
+ ASSERT_TRUE(std::get_if<fvm_probe_interpolated>(&probe_map.data_on({0, 2}).front().info));
probe_handle p0a = get_probe_raw_handle({0, 0}, 0);
probe_handle p0b = get_probe_raw_handle({0, 0}, 1);
probe_handle p1a = get_probe_raw_handle({0, 1}, 0);
probe_handle p1b = get_probe_raw_handle({0, 1}, 1);
- probe_handle p2 = get_probe_raw_handle({0, 2});
+ probe_handle p2a = get_probe_raw_handle({0, 2}, 0);
+ probe_handle p2b = get_probe_raw_handle({0, 2}, 1);
// Ball-and-stick cell with default discretization policy should
// have three CVs, one for branch 0, one trivial one covering the
@@ -172,7 +175,7 @@ void run_v_i_probe_test(const context& ctx) {
// Expect initial raw probe handle values to be the resting potential for
// the voltage probes (cell membrane potential should be constant), and
- // zero for the current probe.
+ // zero for the current probe (including stimulus component).
fvm_value_type resting = voltage[0];
EXPECT_NE(0.0, resting);
@@ -181,7 +184,8 @@ void run_v_i_probe_test(const context& ctx) {
EXPECT_EQ(resting, deref(p0b));
EXPECT_EQ(resting, deref(p1a));
EXPECT_EQ(resting, deref(p1a));
- EXPECT_EQ(0.0, deref(p2));
+ EXPECT_EQ(0.0, deref(p2a));
+ EXPECT_EQ(0.0, deref(p2b));
// After an integration step, expect voltage probe values
// to differ from resting, and for there to be a non-zero current.
@@ -192,7 +196,8 @@ void run_v_i_probe_test(const context& ctx) {
EXPECT_NE(resting, deref(p0b));
EXPECT_NE(resting, deref(p1a));
EXPECT_NE(resting, deref(p1b));
- EXPECT_NE(0.0, deref(p2));
+ EXPECT_NE(0.0, deref(p2a));
+ EXPECT_NE(0.0, deref(p2b));
}
template <typename Backend>
@@ -770,7 +775,7 @@ void run_partial_density_probe_test(const context& ctx) {
template <typename Backend>
void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
// On a passive cable in steady-state, the capacitive membrane current will be zero,
- // and the axial currents should balance the ionic membrane currents in any CV.
+ // and the axial currents should balance the stimulus and ionic membrane currents in any CV.
//
// Membrane currents not associated with an ion are not visible, so we use a custom
// mechanism 'ca_linear' that presents a passive, constant conductance membrane current
@@ -802,8 +807,8 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
cable_cell cell(m, {}, d);
- // Place axial current probes at CV boundaries and make a cell-wide probe for
- // total ionic membrane current.
+ // Place axial current probes at CV boundaries and make cell-wide probes for
+ // total ionic membrane current and stimulus currents.
mlocation_list cv_boundaries;
util::assign(cv_boundaries,
@@ -818,8 +823,9 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
rec.add_probe(0, 0, cable_probe_axial_current{loc});
}
- // Use tag 1 to indicate it's a whole-cell probe.
+ // Use tags 1 and 2 for the whole-cell probes.
rec.add_probe(0, 1, cable_probe_total_ion_current_cell{});
+ rec.add_probe(0, 2, cable_probe_stimulus_current_cell{});
partition_hint_map phints = {
{cell_kind::cable, {partition_hint::max_size, partition_hint::max_size, true}}
@@ -829,7 +835,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
// Take a sample at 20 tau, and run sim for just a bit longer.
std::vector<double> i_axial(n_axial_probe);
- std::vector<double> i_memb;
+ std::vector<double> i_memb, i_stim;
sim.add_sampler(all_probes, explicit_schedule({20*tau}),
[&](probe_metadata pm, std::size_t n_sample, const sample_record* samples)
@@ -837,7 +843,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
// Expect exactly one sample.
ASSERT_EQ(1u, n_sample);
- if (pm.tag==1) { // (whole cell probe)
+ if (pm.tag!=0) { // (whole cell probe)
const mcable_list* m = any_cast<const mcable_list*>(pm.meta);
ASSERT_NE(nullptr, m);
// Metadata should comprise one cable per CV.
@@ -847,8 +853,9 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
ASSERT_NE(nullptr, s);
ASSERT_EQ(s->first+n_cv, s->second);
+ auto* i_var = pm.tag==1? &i_memb: &i_stim;
for (const double* p = s->first; p!=s->second; ++p) {
- i_memb.push_back(*p);
+ i_var->push_back(*p);
}
}
else { // axial current probe
@@ -881,7 +888,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
double net_axial_flux = i<n_axial_probe? i_axial[i]: 0;
net_axial_flux -= i>0? i_axial[i-1]: 0;
- EXPECT_TRUE(testing::near_relative(net_axial_flux, -i_memb[i], 1e-6));
+ EXPECT_TRUE(testing::near_relative(net_axial_flux, -i_memb[i]-i_stim[i], 1e-6));
}
}
@@ -1010,6 +1017,7 @@ void run_v_sampled_probe_test(const context& ctx) {
EXPECT_NE(trace0[1].v, trace1[1].v);
}
+
template <typename Backend>
void run_total_current_probe_test(const context& ctx) {
// Model two passive Y-shaped cells with a similar but not identical
@@ -1046,6 +1054,7 @@ void run_total_current_probe_test(const context& ctx) {
trace_data<std::vector<double>, mcable_list> traces[2];
trace_data<std::vector<double>, mcable_list> ion_traces[2];
+ trace_data<std::vector<double>, mcable_list> stim_traces[2];
// Run the cells sampling at Ï„ and 20Ï„ for both total membrane
// current and total membrane ionic current.
@@ -1053,7 +1062,6 @@ void run_total_current_probe_test(const context& ctx) {
auto run_cells = [&](bool interior_forks) {
auto flags = interior_forks? cv_policy_flag::interior_forks: cv_policy_flag::none;
cv_policy policy = cv_policy_fixed_per_branch(n_cv_per_branch, flags);
- //for (auto& c: cells) { c.discretization() = policy; }
d0.set_default(policy);
d1.set_default(policy);
std::vector<cable_cell> cells = {{m, {}, d0}, {m, {}, d1}};
@@ -1070,6 +1078,9 @@ void run_total_current_probe_test(const context& ctx) {
ion_traces[i] = run_simple_sampler<std::vector<double>, mcable_list>(ctx, t_end, cells, i,
cable_probe_total_ion_current_cell{}, {tau, 20*tau}).at(0);
+ stim_traces[i] = run_simple_sampler<std::vector<double>, mcable_list>(ctx, t_end, cells, i,
+ cable_probe_stimulus_current_cell{}, {tau, 20*tau}).at(0);
+
ASSERT_EQ(2u, traces[i].size());
ASSERT_EQ(2u, ion_traces[i].size());
@@ -1083,26 +1094,30 @@ void run_total_current_probe_test(const context& ctx) {
// same support and same metadata.
ASSERT_EQ((n_cv_per_branch+(int)interior_forks)*n_branch, traces[i].meta.size());
- EXPECT_EQ(ion_traces[i].meta, traces[i].meta);
+ ASSERT_EQ(ion_traces[i].meta, traces[i].meta);
EXPECT_EQ(ion_traces[i][0].v.size(), traces[i][0].v.size());
EXPECT_EQ(ion_traces[i][1].v.size(), traces[i][1].v.size());
- // Check total membrane currents are individually non-zero, but sum is, both
- // at t=Ï„ (j=0) and t=20Ï„ (j=1).
+ // Check total membrane currents + stimulus currents are individually non-zero, but sum is,
+ // both at t=Ï„ (j=0) and t=20Ï„ (j=1).
+
+ ASSERT_EQ(ion_traces[i].meta, stim_traces[i].meta);
for (unsigned j: {0u, 1u}) {
- double max_abs_i_memb = 0;
- double sum_i_memb = 0;
- for (auto i_memb: traces[i][j].v) {
- EXPECT_NE(0.0, i_memb);
- max_abs_i_memb = std::max(max_abs_i_memb, std::abs(i_memb));
- sum_i_memb += i_memb;
+ double max_abs_current = 0;
+ double sum_current = 0;
+ for (auto k: util::count_along(traces[i].meta)) {
+ double current = traces[i][j].v[k] + stim_traces[i][j].v[k];
+
+ EXPECT_NE(0.0, current);
+ max_abs_current = std::max(max_abs_current, std::abs(current));
+ sum_current += current;
}
- EXPECT_NEAR(0.0, sum_i_memb, 1e-6*max_abs_i_memb);
+ ASSERT_NEAR(0.0, sum_current, 1e-6*max_abs_current);
}
- // Confirm that total and ion currents differ at Ï„ but are close at 20Ï„.
+ // Confirm that total (non-stim) and ion currents differ at Ï„ but are close at 20Ï„.
for (unsigned k = 0; k<traces[i].size(); ++k) {
const double rtol_large = 1e-3;
@@ -1134,6 +1149,52 @@ void run_total_current_probe_test(const context& ctx) {
}
}
+
+template <typename Backend>
+void run_stimulus_probe_test(const context& ctx) {
+ // Model two simple stick cable cells, 3 CVs each, and stimuli on cell 0, cv 1
+ // and cell 1, cv 2. Run both cells in the same cell group.
+
+ const double stim_until = 1.; // [ms]
+ auto m = make_stick_morphology();
+ cv_policy policy = cv_policy_fixed_per_branch(3);
+
+ decor d0, d1;
+ d0.set_default(policy);
+ d0.place(mlocation{0, 0.5}, i_clamp(0., stim_until, 10.));
+ d0.place(mlocation{0, 0.5}, i_clamp(0., stim_until, 20.));
+ double expected_stim0 = 30;
+
+ d1.set_default(policy);
+ d1.place(mlocation{0, 1}, i_clamp(0., stim_until, 30.));
+ d1.place(mlocation{0, 1}, i_clamp(0., stim_until, -10.));
+ double expected_stim1 = 20;
+
+ std::vector<cable_cell> cells = {{m, {}, d0}, {m, {}, d1}};
+
+ // Sample the cells during the stimulus, and after.
+
+ trace_data<std::vector<double>, mcable_list> traces[2];
+
+ for (unsigned i: {0u, 1u}) {
+ traces[i] = run_simple_sampler<std::vector<double>, mcable_list>(ctx, 2.5*stim_until, cells, i,
+ cable_probe_stimulus_current_cell{}, {stim_until/2, 2*stim_until}).at(0);
+
+ ASSERT_EQ(3u, traces[i].meta.size());
+ for (unsigned cv: {0u, 1u, 2u}) {
+ ASSERT_EQ(2u, traces[i].size());
+ }
+ }
+
+ // Every sample in each trace should be zero _except_ the first sample for cell 0, cv 1
+ // and the first sample for cell 1, cv 2.
+
+ EXPECT_EQ((std::vector<double>{0, expected_stim0, 0}), traces[0][0]);
+ EXPECT_EQ((std::vector<double>{0, 0, expected_stim1}), traces[1][0]);
+ EXPECT_EQ((std::vector<double>(3)), traces[0][1]);
+ EXPECT_EQ((std::vector<double>(3)), traces[1][1]);
+}
+
template <typename Backend>
void run_exact_sampling_probe_test(const context& ctx) {
// As the exact sampling implementation interacts with the event delivery
diff --git a/test/unit/test_spikes.cpp b/test/unit/test_spikes.cpp
index d207f586c7adff0180f456ea70126cb83ad57f9e..9ededc628810435abcfb5bc1f591300c22b7f3dc 100644
--- a/test/unit/test_spikes.cpp
+++ b/test/unit/test_spikes.cpp
@@ -22,9 +22,11 @@ using namespace arb;
#ifndef USE_BACKEND
using backend = multicore::backend;
#define SPIKES_TEST_CLASS spikes
+#define GPU_ID -1
#else
using backend = USE_BACKEND;
#define SPIKES_TEST_CLASS spikes_gpu
+#define GPU_ID 0
#endif
TEST(SPIKES_TEST_CLASS, threshold_watcher) {
@@ -211,7 +213,7 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) {
dict.set("mid", arb::ls::on_branches(0.5));
arb::proc_allocation resources;
- resources.gpu_id = arbenv::default_gpu();
+ resources.gpu_id = GPU_ID;
auto context = arb::make_context(resources);
std::vector<arb::spike> spikes;