From 05569a168e333ba32dd5f3d66110346f809b9ec8 Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Thu, 8 Apr 2021 09:30:09 +0200
Subject: [PATCH] Add phase parameter to current clamps. (#1474)
* `i_clamp` constructors take new optional phase argument.
* Pulse or box envelope clamps are now provided by the static method `i_clamp::box`, taking onset, duration, and amplitude parameters.
* Bump pybind11 version so we can have keyword-only arguments in the Python API.
* Make frequency and phase keyword-only arguments to Python `iclamp`.
* Change frequency units for clamps from Hz to kHz for consistency with other interfaces.
Fixes #1454.
---
arbor/backends/gpu/shared_state.cpp | 2 ++
arbor/backends/gpu/shared_state.hpp | 3 ++-
arbor/backends/gpu/stimulus.cu | 4 ++--
arbor/backends/gpu/stimulus.hpp | 1 +
arbor/backends/multicore/shared_state.cpp | 7 +++---
arbor/backends/multicore/shared_state.hpp | 3 ++-
arbor/fvm_layout.cpp | 3 +++
arbor/fvm_layout.hpp | 3 ++-
arbor/include/arbor/cable_cell_param.hpp | 19 +++++++++++------
arborio/cableio.cpp | 20 ++++++++---------
doc/concepts/decor.rst | 15 +++++++------
doc/fileformat/cable_cell.rst | 26 +++++++++++------------
doc/python/decor.rst | 2 +-
example/gap_junctions/gap_junctions.cpp | 2 +-
example/probe-demo/probe-demo.cpp | 2 +-
python/cells.cpp | 20 +++++++++--------
python/pybind11 | 2 +-
test/unit/test_fvm_lowered.cpp | 11 +++++-----
test/unit/test_mc_cell_group.cpp | 2 +-
test/unit/test_mc_cell_group_gpu.cpp | 2 +-
test/unit/test_probe.cpp | 18 ++++++++--------
test/unit/test_s_expr.cpp | 15 +++++++------
test/unit/test_spikes.cpp | 2 +-
23 files changed, 102 insertions(+), 82 deletions(-)
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index 9dfa021d..555f9a28 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -101,6 +101,7 @@ istim_state::istim_state(const fvm_stimulus_config& stim) {
std::vector<fvm_index_type> edivs;
frequency_ = make_const_view(stim.frequency);
+ phase_ = make_const_view(stim.phase);
arb_assert(n==frequency_.size());
arb_assert(n==stim.envelope_time.size());
@@ -132,6 +133,7 @@ istim_state::istim_state(const fvm_stimulus_config& stim) {
ppack_.accu_index = accu_index_.data();
ppack_.accu_to_cv = accu_to_cv_.data();
ppack_.frequency = frequency_.data();
+ ppack_.phase = phase_.data();
ppack_.envl_amplitudes = envl_amplitudes_.data();
ppack_.envl_times = envl_times_.data();
ppack_.envl_divs = envl_divs_.data();
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index 6243982e..695821fa 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -66,7 +66,8 @@ struct istim_state {
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 frequency_; // (kHz) stimulus frequency per instance.
+ array phase_; // (rad) stimulus waveform phase at t=0.
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,
diff --git a/arbor/backends/gpu/stimulus.cu b/arbor/backends/gpu/stimulus.cu
index cc3f85ca..a127128c 100644
--- a/arbor/backends/gpu/stimulus.cu
+++ b/arbor/backends/gpu/stimulus.cu
@@ -17,7 +17,7 @@ namespace kernel {
__global__
void istim_add_current_impl(int n, istim_pp pp) {
- constexpr double freq_scale = 2*pi*0.001;
+ constexpr double two_pi = 2*pi;
auto i = threadIdx.x + blockDim.x*blockIdx.x;
if (i>=n) return;
@@ -43,7 +43,7 @@ void istim_add_current_impl(int n, istim_pp pp) {
}
if (double f = pp.frequency[i]) {
- J *= std::sin(freq_scale*f*t);
+ J *= std::sin(two_pi*f*t + pp.phase[i]);
}
gpu_atomic_add(&pp.accu_stim[ai], J);
diff --git a/arbor/backends/gpu/stimulus.hpp b/arbor/backends/gpu/stimulus.hpp
index e03edd43..97447f6e 100644
--- a/arbor/backends/gpu/stimulus.hpp
+++ b/arbor/backends/gpu/stimulus.hpp
@@ -12,6 +12,7 @@ struct istim_pp {
const fvm_index_type* accu_index;
const fvm_index_type* accu_to_cv;
const fvm_value_type* frequency;
+ const fvm_value_type* phase;
const fvm_value_type* envl_amplitudes;
const fvm_value_type* envl_times;
const fvm_index_type* envl_divs;
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 0da23bca..442da337 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -92,7 +92,8 @@ void ion_state::reset() {
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))
+ frequency_(stim.frequency.begin(), stim.frequency.end(), pad(alignment)),
+ phase_(stim.phase.begin(), stim.phase.end(), pad(alignment))
{
using util::assign;
@@ -138,7 +139,7 @@ void istim_state::reset() {
}
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;
+ constexpr double two_pi = 2*math::pi<double>;
// Consider vectorizing...
for (auto i: util::count_along(accu_index_)) {
@@ -168,7 +169,7 @@ void istim_state::add_current(const array& time, const iarray& cv_to_intdom, arr
}
if (frequency_[i]) {
- J *= std::sin(freq_scale*frequency_[i]*t);
+ J *= std::sin(two_pi*frequency_[i]*t + phase_[i]);
}
accu_stim_[ai] += J;
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index 2011c668..54c0424b 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -80,7 +80,8 @@ struct istim_state {
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 frequency_; // (kHz) stimulus frequency per instance.
+ array phase_; // (rad) stimulus waveform phase at t=0.
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,
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 4f19daa7..1f917939 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -753,6 +753,7 @@ 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.phase, right.stimuli.phase);
append(left.stimuli.envelope_time, right.stimuli.envelope_time);
append(left.stimuli.envelope_amplitude, right.stimuli.envelope_amplitude);
@@ -1096,6 +1097,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
config.cv.reserve(n);
config.cv_unique.reserve(n);
config.frequency.reserve(n);
+ config.phase.reserve(n);
config.envelope_time.reserve(n);
config.envelope_amplitude.reserve(n);
@@ -1106,6 +1108,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
config.cv.push_back(cv);
config.frequency.push_back(stim.frequency);
+ config.phase.push_back(stim.phase);
std::size_t envl_n = stim.envelope.size();
std::vector<double> envl_t, envl_a;
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 403d52d4..8d8ba818 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -279,7 +279,8 @@ struct fvm_stimulus_config {
// 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<double> frequency; // [kHz]
+ std::vector<double> phase; // [rad]
std::vector<std::vector<double>> envelope_time; // [ms]
std::vector<std::vector<double>> envelope_amplitude; // [A/m²]
};
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index c21649a6..0e80bb85 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -54,27 +54,32 @@ struct i_clamp {
};
std::vector<envelope_point> envelope;
- double frequency = 0; // [Hz] 0 => constant
+ double frequency = 0; // [kHz] 0 => constant
+ double phase = 0; // [rad]
// A default constructed i_clamp, with empty envelope, describes
// a trivial stimulus, providing no current at all.
i_clamp() = default;
// The simple constructor describes a constant amplitude stimulus starting from t=0.
- explicit i_clamp(double amplitude, double frequency = 0):
+ explicit i_clamp(double amplitude, double frequency = 0, double phase = 0):
envelope({{0., amplitude}}),
- frequency(frequency)
+ frequency(frequency),
+ phase(phase)
{}
// Describe a stimulus by envelope and frequency.
- explicit i_clamp(std::vector<envelope_point> envelope, double frequency = 0):
+ explicit i_clamp(std::vector<envelope_point> envelope, double frequency = 0, double phase = 0):
envelope(std::move(envelope)),
- frequency(frequency)
+ frequency(frequency),
+ phase(phase)
{}
// 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) {}
+ static i_clamp box(double onset, double duration, double amplitude, double frequency = 0, double phase = 0) {
+ return i_clamp({{onset, amplitude}, {onset+duration, amplitude}, {onset+duration, 0.}}, frequency, phase);
+ }
+
};
// Threshold detector description.
diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp
index 9a7ff5af..d87b1413 100644
--- a/arborio/cableio.cpp
+++ b/arborio/cableio.cpp
@@ -61,7 +61,7 @@ s_expr mksexp(const i_clamp& c) {
std::transform(c.envelope.begin(), c.envelope.end(), std::back_inserter(evlps),
[](const auto& x){return slist(x.t, x.amplitude);});
auto envelope = slist("envelope"_symbol, slist_range(evlps));
- return slist("current-clamp"_symbol, envelope, c.frequency);
+ return slist("current-clamp"_symbol, envelope, c.frequency, c.phase);
}
s_expr mksexp(const threshold_detector& d) {
return slist("threshold-detector"_symbol, d.threshold);
@@ -234,14 +234,14 @@ std::vector<arb::i_clamp::envelope_point> make_envelope(const std::vector<std::v
});
return envlp;
}
-arb::i_clamp make_i_clamp(const std::vector<arb::i_clamp::envelope_point>& envlp, double freq) {
- return arb::i_clamp(envlp, freq);
+arb::i_clamp make_i_clamp(const std::vector<arb::i_clamp::envelope_point>& envlp, double freq, double phase) {
+ return arb::i_clamp(envlp, freq, phase);
}
pulse_tuple make_envelope_pulse(double delay, double duration, double amplitude) {
return pulse_tuple{delay, duration, amplitude};
}
-arb::i_clamp make_i_clamp_pulse(pulse_tuple p, double freq) {
- return arb::i_clamp(std::get<0>(p), std::get<1>(p), std::get<2>(p), freq);
+arb::i_clamp make_i_clamp_pulse(pulse_tuple p, double freq, double phase) {
+ return arb::i_clamp::box(std::get<0>(p), std::get<1>(p), std::get<2>(p), freq, phase);
}
arb::gap_junction_site make_gap_junction_site() {
return arb::gap_junction_site{};
@@ -757,10 +757,10 @@ eval_map named_evals{
"`envelope` with one or more pairs of start time and amplitude (start:real amplitude:real)")},
{"envelope-pulse", make_call<double, double, double>(make_envelope_pulse,
"'envelope-pulse' with 3 arguments (delay:real duration:real amplitude:real)")},
- {"current-clamp", make_call<std::vector<arb::i_clamp::envelope_point>, double>(make_i_clamp,
- "`current-clamp` with 2 arguments (env:envelope freq:real)")},
- {"current-clamp", make_call<pulse_tuple, double>(make_i_clamp_pulse,
- "`current-clamp` with 2 arguments (env:envelope_pulse freq:real)")},
+ {"current-clamp", make_call<std::vector<arb::i_clamp::envelope_point>, double, double>(make_i_clamp,
+ "`current-clamp` with 3 arguments (env:envelope freq:real phase:real)")},
+ {"current-clamp", make_call<pulse_tuple, double, double>(make_i_clamp_pulse,
+ "`current-clamp` with 3 arguments (env:envelope_pulse freq:real phase:real)")},
{"threshold-detector", make_call<double>(make_threshold_detector,
"'threshold-detector' with 1 argument (threshold:real)")},
{"gap-junction-site", make_call<>(make_gap_junction_site,
@@ -857,4 +857,4 @@ parse_hopefully<cable_cell_component> parse_component(const std::string& s) {
parse_hopefully<cable_cell_component> parse_component(std::istream& s) {
return parse_component(std::string(std::istreambuf_iterator<char>(s), {}));
}
-} // namespace arborio
\ No newline at end of file
+} // namespace arborio
diff --git a/doc/concepts/decor.rst b/doc/concepts/decor.rst
index 93f6c0b5..9b29d0eb 100644
--- a/doc/concepts/decor.rst
+++ b/doc/concepts/decor.rst
@@ -290,9 +290,10 @@ See :ref:`modelgapjunctions`.
A current stimulus is a DC or sinusoidal current of fixed frequency 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 is described by three parameters:
+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;
+a frequency in kilohertz, where a value of zero denotes DC; and the phase in radians at time zero.
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
@@ -307,11 +308,11 @@ constant stimuli and constant amplitude stimuli restricted to a fixed time inter
# 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))
+ # Constant amplitude 10 nA stimulus at 20 Hz, with initial phase of π/4 radians.
+ decor.place('(root)', arbor.iclamp(10, frequency=0.020, phasce=math.pi/4))
- # Stimulus at 20 Hz, amplitude 10 nA, for 40 ms starting at t = 30 ms.
- decor.place('(root)', arbor.iclamp(30, 40, 10, 20))
+ # Stimulus at 1 kHz, amplitude 10 nA, for 40 ms starting at t = 30 ms.
+ decor.place('(root)', arbor.iclamp(30, 40, 10, frequency=1))
# Piecewise linear stimulus with amplitude ranging from 0 nA to 10 nA,
# starting at t = 30 ms and stopping at t = 50 ms.
diff --git a/doc/fileformat/cable_cell.rst b/doc/fileformat/cable_cell.rst
index 5419912e..8bb9f1ff 100644
--- a/doc/fileformat/cable_cell.rst
+++ b/doc/fileformat/cable_cell.rst
@@ -148,29 +148,29 @@ The various properties and dynamics of the decor are described as follows:
(ion-reversal-potential-method "ca" (mechanism "nernst/ca"))
-.. label:: (current-clamp (envelope-pulse delay:real duration:real amplitude:real) freq:real)
+.. label:: (current-clamp (envelope-pulse delay:real duration:real amplitude:real) freq:real phase:real)
- This creates a *current clamp*. If the frequency ``freq`` (unit Hz) is zero, the current is a square
+ This creates a *current clamp*. If the frequency ``freq`` (unit kHz) is zero, the current is a square
pulse with amplitude ``amplitude`` (unit nA) starting at ``delay`` (unit ms) and lasting for ``duration``
(unit ms). If ``freq`` is non-zero, the current is sinusoidal with amplitude ``amplitude`` and frequency
- ``freq`` from time ``delay`` and lasting for ``duration``.
+ ``freq`` from time ``delay`` and lasting for ``duration``, with phase ``phase`` (unit rad) at time zero.
(More information about current clamps can be found :ref:`here <cablecell-stimuli>`).
-.. label:: (current-clamp [...(envelope time:real amplitude:real)] freq:real)
+.. label:: (current-clamp [...(envelope time:real amplitude:real)] freq:real phase:real)
This creates a *current clamp* with an amplitude governed by the given envelopes (``time`` unit ms and
- ``amplitude`` unit nA). A frequency ``freq`` (unit Hz) of zero implies that the generated current simply
+ ``amplitude`` unit nA). A frequency ``freq`` (unit kHz) of zero implies that the generated current simply
follows the envelope. A non-zero ``freq`` implies the current is sinusoidal with that frequency and amplitude
- that varies according to the envelope. (More information about current clamps can be found
- :ref:`here <cablecell-stimuli>`).
+ that varies according to the envelope. The ``phase`` (unit rad) is the phase of the sinusoidal current
+ clamp at time zero. (More information about current clamps can be found :ref:`here <cablecell-stimuli>`).
For example:
.. code::
- (current-clamp (envelope (0 10) (50 10) (50 0)) 40)
+ (current-clamp (envelope (0 10) (50 10) (50 0)) 0.04 0.15)
- This expression describes a sinusoidal current with amplitude 10nA and frequency 40Hz and that lasts
- from t = 0ms to t = 50ms, finally leaving the current at 0nA (final amplitude in the envelope).
+ This expression describes a sinusoidal current with amplitude 10Â nA and frequency 40Â Hz and that lasts
+ from t = 0Â ms to t = 50Â ms, finally leaving the current at 0Â nA (final amplitude in the envelope).
.. label:: (threshold-detector val:real).
@@ -205,7 +205,7 @@ The various properties and dynamics of the decor are described as follows:
(place (locset "mylocset") (threshold-detector 10))
- This expression places a 10 mV threshold detector on the locset labeled ``mylocset``.
+ This expression places a 10Â mV threshold detector on the locset labeled ``mylocset``.
(The definition of ``mylocset`` should be provided in a label dictionary associated
with the decor).
@@ -219,7 +219,7 @@ The various properties and dynamics of the decor are described as follows:
(default (membrane-potential -65))
- This expression sets the default membrane potential of the cell to -65 mV.
+ This expression sets the default membrane potential of the cell to -65Â mV.
Any number of paint, place and default expressions can be used to create a decor as follows:
@@ -437,4 +437,4 @@ API
---
* :ref:`Python <pycablecellformat>`
-* :ref:`C++ <cppcablecellformat>`
\ No newline at end of file
+* :ref:`C++ <cppcablecellformat>`
diff --git a/doc/python/decor.rst b/doc/python/decor.rst
index 099bfebd..7a7ca937 100644
--- a/doc/python/decor.rst
+++ b/doc/python/decor.rst
@@ -167,7 +167,7 @@ Cable cell decoration
:param str locations: description of the locset.
:param stim: the current stim.
- :type stim: :py:class:`i_clamp`
+ :type stim: :py:class:`iclamp`
:rtype: int
.. method:: place(locations, d)
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 4c16d2bc..58f37a4d 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -320,7 +320,7 @@ arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration)
// Attach a stimulus to the second cell.
if (!gid) {
- arb::i_clamp stim(0, stim_duration, 0.4);
+ auto stim = arb::i_clamp::box(0, stim_duration, 0.4);
decor.place(arb::mlocation{0, 0.5}, stim);
}
diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp
index 9e9c1600..d2391af4 100644
--- a/example/probe-demo/probe-demo.cpp
+++ b/example/probe-demo/probe-demo.cpp
@@ -122,7 +122,7 @@ struct cable_recipe: public arb::recipe {
arb::decor decor;
decor.paint(arb::reg::all(), "hh"); // HH mechanism over whole cell.
- decor.place(arb::mlocation{0, 0.}, arb::i_clamp{0., INFINITY, 1.}); // Inject a 1 nA current indefinitely.
+ decor.place(arb::mlocation{0, 0.}, arb::i_clamp{1.}); // Inject a 1 nA current indefinitely.
return arb::cable_cell(tree, {}, decor);
}
diff --git a/python/cells.cpp b/python/cells.cpp
index 76afd9e5..a9f7fabd 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -288,24 +288,25 @@ void register_cells(pybind11::module& m) {
"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, double frequency) {
- return arb::i_clamp{ts, dur, cur, frequency};
- }), "tstart"_a, "duration"_a, "current"_a, "frequency"_a=0,
+ [](double ts, double dur, double cur, double frequency, double phase) {
+ return arb::i_clamp::box(ts, dur, cur, frequency, phase);
+ }), "tstart"_a, "duration"_a, "current"_a, pybind11::kw_only(), "frequency"_a=0, "phase"_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,
+ [](double cur, double frequency, double phase) {
+ return arb::i_clamp{cur, frequency, phase};
+ }), "current"_a, pybind11::kw_only(), "frequency"_a=0, "phase"_a=0,
"Construct constant amplitude current clamp")
.def(pybind11::init(
- [](std::vector<std::pair<double, double>> envl, double frequency) {
+ [](std::vector<std::pair<double, double>> envl, double frequency, double phase) {
arb::i_clamp clamp;
for (const auto& p: envl) {
clamp.envelope.push_back({p.first, p.second});
}
clamp.frequency = frequency;
+ clamp.phase = phase;
return clamp;
- }), "envelope"_a, "frequency"_a=0,
+ }), "envelope"_a, pybind11::kw_only(), "frequency"_a=0, "phase"_a=0,
"Construct current clamp according to (time, amplitude) linear envelope")
.def_property_readonly("envelope",
[](const arb::i_clamp& obj) {
@@ -316,7 +317,8 @@ void register_cells(pybind11::module& m) {
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_readonly("frequency", &arb::i_clamp::frequency, "Oscillation frequency (kHz), zero implies DC stimulus.")
+ .def_readonly("phase", &arb::i_clamp::phase, "Oscillation initial phase (rad)")
.def("__repr__", [](const arb::i_clamp& c) {
return util::pprintf("<arbor.iclamp: frequency {} Hz>", c.frequency);})
.def("__str__", [](const arb::i_clamp& c) {
diff --git a/python/pybind11 b/python/pybind11
index 3b1dbeba..8de7772c 160000
--- a/python/pybind11
+++ b/python/pybind11
@@ -1 +1 @@
-Subproject commit 3b1dbebabc801c9cf6f0953a4c20b904d444f879
+Subproject commit 8de7772cc72daca8e947b79b83fea46214931604
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index a35e8393..0634c25f 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -345,9 +345,9 @@ TEST(fvm_lowered, stimulus) {
auto desc = make_cell_ball_and_stick(false);
// At end of stick
- desc.decorations.place(mlocation{0,1}, i_clamp{5., 80., 0.3});
+ desc.decorations.place(mlocation{0,1}, i_clamp::box(5., 80., 0.3));
// On the soma CV, which is over the approximate interval: (cable 0 0 0.1)
- desc.decorations.place(mlocation{0,0.05}, i_clamp{1., 2., 0.1});
+ desc.decorations.place(mlocation{0,0.05}, i_clamp::box(1., 2., 0.1));
std::vector<cable_cell> cells{desc};
@@ -416,12 +416,13 @@ TEST(fvm_lowered, ac_stimulus) {
segment_tree tree;
tree.append(mnpos, {0., 0., 0., 1.}, {100., 0., 0., 1.}, 1);
- const double freq = 20; // (Hz)
+ const double freq = 0.02; // (kHz)
+ const double phase = 0.1; // (radian)
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));
+ dec.place(mlocation{0, 0}, i_clamp({{0, 0}, {max_time, max_amplitude}, {max_time, 0}}, freq, phase));
std::vector<cable_cell> cells = {cable_cell(tree, {}, dec)};
cable_cell_global_properties gprop;
@@ -455,7 +456,7 @@ TEST(fvm_lowered, ac_stimulus) {
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;
+ double expected_I = t<=max_time? max_amplitude*t/max_time*std::sin(2*math::pi<double>*t*freq+phase): 0;
EXPECT_TRUE(near_relative(-expected_I, J[0]*A[0]*unit_factor, reltol));
}
}
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index efc5fbff..564b2179 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -28,7 +28,7 @@ namespace {
auto d = builder.make_cell();
d.decorations.paint("soma"_lab, "hh");
d.decorations.paint("dend"_lab, "pas");
- d.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
+ d.decorations.place(builder.location({1,1}), i_clamp::box(5, 80, 0.3));
d.decorations.place(builder.location({0, 0}), threshold_detector{0});
return d;
}
diff --git a/test/unit/test_mc_cell_group_gpu.cpp b/test/unit/test_mc_cell_group_gpu.cpp
index 7b2306e6..f7c8fbb7 100644
--- a/test/unit/test_mc_cell_group_gpu.cpp
+++ b/test/unit/test_mc_cell_group_gpu.cpp
@@ -29,7 +29,7 @@ namespace {
auto d = builder.make_cell();
d.decorations.paint("soma"_lab, "hh");
d.decorations.paint("dend"_lab, "pas");
- d.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
+ d.decorations.place(builder.location({1,1}), i_clamp::box(5, 80, 0.3));
d.decorations.place(builder.location({0, 0}), threshold_detector{0});
return d;
}
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 7e705019..1ac01253 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -110,7 +110,7 @@ void run_v_i_probe_test(const context& ctx) {
bs.decorations.set_default(cv_policy_fixed_per_branch(1));
- i_clamp stim(0, 100, 0.3);
+ auto stim = i_clamp::box(0, 100, 0.3);
bs.decorations.place(mlocation{1, 1}, stim);
cable1d_recipe rec((cable_cell(bs)));
@@ -792,7 +792,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
cv_policy policy = cv_policy_fixed_per_branch(n_cv);
d.set_default(policy);
- d.place(mlocation{0, 0}, i_clamp(0, INFINITY, 0.3));
+ d.place(mlocation{0, 0}, i_clamp(0.3));
// The time constant will be membrane capacitance / membrane conductance.
// For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance
@@ -989,8 +989,8 @@ void run_v_sampled_probe_test(const context& ctx) {
// samples at the same point on each cell will give the same value at
// 0.3 ms, but different at 0.6 ms.
- d0.place(mlocation{1, 1}, i_clamp(0, 0.5, 1.));
- d1.place(mlocation{1, 1}, i_clamp(0, 1.0, 1.));
+ d0.place(mlocation{1, 1}, i_clamp::box(0, 0.5, 1.));
+ d1.place(mlocation{1, 1}, i_clamp::box(0, 1.0, 1.));
mlocation probe_loc{1, 0.2};
std::vector<cable_cell> cells = {{bs.morph, bs.labels, d0}, {bs.morph, bs.labels, d1}};
@@ -1040,7 +1040,7 @@ void run_total_current_probe_test(const context& ctx) {
// to 0.01 F/m².
const double tau = 0.1; // [ms]
- d0.place(mlocation{0, 0}, i_clamp(0, INFINITY, 0.3));
+ d0.place(mlocation{0, 0}, i_clamp(0.3));
d0.paint(reg::all(), mechanism_desc("ca_linear").set("g", 0.01)); // [S/cm²]
d0.set_default(membrane_capacitance{0.01}); // [F/m²]
@@ -1161,13 +1161,13 @@ void run_stimulus_probe_test(const context& ctx) {
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.));
+ d0.place(mlocation{0, 0.5}, i_clamp::box(0., stim_until, 10.));
+ d0.place(mlocation{0, 0.5}, i_clamp::box(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.));
+ d1.place(mlocation{0, 1}, i_clamp::box(0., stim_until, 30.));
+ d1.place(mlocation{0, 1}, i_clamp::box(0., stim_until, -10.));
double expected_stim1 = 20;
std::vector<cable_cell> cells = {{m, {}, d0}, {m, {}, d1}};
diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp
index 543f7779..078d531f 100644
--- a/test/unit/test_s_expr.cpp
+++ b/test/unit/test_s_expr.cpp
@@ -317,7 +317,7 @@ std::ostream& operator<<(std::ostream& o, const i_clamp& c) {
for (const auto& p: c.envelope) {
o << " (" << p.t << " " << p.amplitude << ')';
}
- return o << ") " << c.frequency << ')';
+ return o << ") " << c.frequency << ' ' << c.phase << ')';
}
std::ostream& operator<<(std::ostream& o, const threshold_detector& p) {
return o << "(threshold-detector " << p.threshold << ')';
@@ -463,7 +463,7 @@ TEST(decor_literals, round_tripping) {
auto default_literals = {
"(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))"};
auto place_literals = {
- "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0)",
+ "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 10 0.25)",
"(threshold-detector -10)",
"(gap-junction-site)",
"(mechanism \"expsyn\")"};
@@ -480,8 +480,8 @@ TEST(decor_literals, round_tripping) {
for (auto l: place_literals) {
EXPECT_EQ(l, round_trip_variant<placeable>(l));
}
- auto clamp_literal = "(current-clamp (envelope-pulse 10 5 0.1) 50)";
- EXPECT_EQ("(current-clamp (envelope (10 0.1) (15 0.1) (15 0)) 50)", round_trip_variant<placeable>(clamp_literal));
+ auto clamp_literal = "(current-clamp (envelope-pulse 10 5 0.1) 50 0.5)";
+ EXPECT_EQ("(current-clamp (envelope (10 0.1) (15 0.1) (15 0)) 50 0.5)", round_trip_variant<placeable>(clamp_literal));
std::string mech_str = "(mechanism \"kamt\" (\"gbar\" 50) (\"zetam\" 0.1) (\"q10\" 5))";
auto maybe_mech = arborio::parse_expression(mech_str);
@@ -524,7 +524,7 @@ TEST(decor_expressions, round_tripping) {
"(default (ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\")))"
};
auto decorate_place_literals = {
- "(place (location 3 0.2) (current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0))",
+ "(place (location 3 0.2) (current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0.5 0.25))",
"(place (terminal) (threshold-detector -10))",
"(place (root) (gap-junction-site))",
"(place (locset \"my!ls\") (mechanism \"expsyn\"))"};
@@ -825,6 +825,7 @@ TEST(cable_cell_literals, errors) {
"(mechanism \"pas\" (\"g\" 0.5 0.1) (\"e\" 0.2))", // too many values
"(gap-junction-site 0)", // too many arguments
"(current-clamp (envelope (10 0.5) (110 0.5) (110 0)))", // too few arguments
+ "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 10)", // too few arguments
"(paint (region) (mechanism \"hh\"))", // invalid region
"(paint (tag 1) (mechanims hh))", // invalid painting
"(paint (terminal) (membrance-capacitance 0.2))", // can't paint a locset
@@ -872,7 +873,7 @@ TEST(doc_expressions, parse) {
"(locset-def \"my_locset\" (location 3 0.5))",
"(mechanism \"hh\" (\"gl\" 0.5) (\"el\" 2))",
"(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))",
- "(current-clamp (envelope (0 10) (50 10) (50 0)) 40)",
+ "(current-clamp (envelope (0 10) (50 10) (50 0)) 40 0.25)",
"(paint (tag 1) (membrane-capacitance 0.02))",
"(place (locset \"mylocset\") (threshold-detector 10))",
"(default (membrane-potential -65))",
@@ -991,4 +992,4 @@ TEST(doc_expressions, parse) {
{
EXPECT_TRUE(arborio::parse_expression(expr));
}
-}
\ No newline at end of file
+}
diff --git a/test/unit/test_spikes.cpp b/test/unit/test_spikes.cpp
index 9ededc62..b2b51586 100644
--- a/test/unit/test_spikes.cpp
+++ b/test/unit/test_spikes.cpp
@@ -222,7 +222,7 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) {
arb::decor decor;
decor.set_default(arb::cv_policy_every_segment());
decor.place("\"mid\"", arb::threshold_detector{10});
- decor.place("\"mid\"", arb::i_clamp(0.01+i*dt, duration, 0.5));
+ decor.place("\"mid\"", arb::i_clamp::box(0.01+i*dt, duration, 0.5));
decor.place("\"mid\"", arb::mechanism_desc("hh"));
arb::cable_cell cell(morpho, dict, decor);
--
GitLab