From efb612dcc0d2934065cef9aa2b9a434b97f9ef81 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 28 Nov 2023 16:31:12 +0100 Subject: [PATCH] Tests compile; poisson offset/reset broken. --- arbor/include/arbor/cable_cell_param.hpp | 53 ++-- arbor/include/arbor/event_generator.hpp | 48 ++-- arbor/include/arbor/schedule.hpp | 229 +++--------------- arbor/include/arbor/simulation.hpp | 4 +- arbor/include/arbor/units.hpp | 2 + arbor/schedule.cpp | 254 +++++++++++++++++--- arbor/simulation.cpp | 10 +- python/schedule.cpp | 170 +++++++------ python/schedule.hpp | 59 ++--- python/simulation.cpp | 2 +- python/single_cell_model.cpp | 15 +- test/common_cells.cpp | 17 +- test/common_cells.hpp | 2 - test/unit-distributed/test_communicator.cpp | 10 +- test/unit/test_cable_cell.cpp | 5 +- test/unit/test_cable_cell_group.cpp | 7 +- test/unit/test_diffusion.cpp | 17 +- test/unit/test_domain_decomposition.cpp | 2 - test/unit/test_event_delivery.cpp | 4 +- test/unit/test_event_generators.cpp | 14 +- test/unit/test_fvm_layout.cpp | 15 +- test/unit/test_fvm_lowered.cpp | 29 +-- test/unit/test_lif_cell_group.cpp | 24 +- test/unit/test_merge_events.cpp | 24 +- test/unit/test_probe.cpp | 57 ++--- test/unit/test_recipe.cpp | 143 +++++------ test/unit/test_schedule.cpp | 136 +++-------- test/unit/test_sde.cpp | 12 +- test/unit/test_serdes.cpp | 14 +- test/unit/test_simulation.cpp | 26 +- test/unit/test_spike_source.cpp | 43 ++-- test/unit/test_spikes.cpp | 10 +- test/unit/test_synapses.cpp | 2 - test/unit/test_v_clamp.cpp | 14 +- 34 files changed, 713 insertions(+), 760 deletions(-) diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp index 333db33510..5f33303a5e 100644 --- a/arbor/include/arbor/cable_cell_param.hpp +++ b/arbor/include/arbor/cable_cell_param.hpp @@ -37,22 +37,29 @@ struct cable_cell_ion_data { std::optional diffusivity; }; -// 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 -// amplitude 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. - +/** + * Current clamp; 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 + * amplitude 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 ARB_SYMBOL_VISIBLE i_clamp { struct envelope_point { + /** + * Current at point in time + * + * @param t, must be convertible to time + * @param amplitude must be convertible to current + */ envelope_point(const units::quantity& t, const units::quantity& amplitude): t(t.value_as(units::ms)), @@ -70,7 +77,14 @@ struct ARB_SYMBOL_VISIBLE i_clamp { // a trivial stimulus, providing no current at all. i_clamp() = default; - // The simple constructor describes a constant amplitude stimulus starting from t=0. + /** + * Constant amplitude stimulus starting at t = 0. + * + * @param amplitude must be convertible to current + * @param frequency, must be convertible to frequency; gives a sine current if not zero + * @param frequency, must be convertible to radians, phase shift of sine. + */ + explicit i_clamp(const units::quantity& amplitude, const units::quantity& frequency = 0*units::kHz, const units::quantity& phase = 0*units::rad): @@ -112,6 +126,15 @@ struct ARB_SYMBOL_VISIBLE init_membrane_potential { iexpr value = NAN; // [mV] }; +struct ARB_SYMBOL_VISIBLE init_membrane_potential_ { + double value = NAN; // [mV] + iexpr scale = 1; // [1] + + init_membrane_potential_(const units::quantity& m, iexpr scale=1): value(m.value_as(units::mV)), scale{scale} {} + +}; + + struct ARB_SYMBOL_VISIBLE temperature_K { iexpr value = NAN; // [K] }; diff --git a/arbor/include/arbor/event_generator.hpp b/arbor/include/arbor/event_generator.hpp index c378ef7d03..cab7e5fe2b 100644 --- a/arbor/include/arbor/event_generator.hpp +++ b/arbor/include/arbor/event_generator.hpp @@ -1,10 +1,5 @@ #pragma once -#include -#include -#include -#include -#include #include #include @@ -101,34 +96,29 @@ struct event_generator { inline event_generator empty_generator( cell_local_label_type target, - float weight) -{ + float weight) { return event_generator(std::move(target), weight, schedule()); } // Generate events at integer multiples of dt that lie between tstart and tstop. -inline event_generator regular_generator( - cell_local_label_type target, - float weight, - time_type tstart, - time_type dt, - time_type tstop=terminal_time) -{ +inline event_generator regular_generator(cell_local_label_type target, + float weight, + const units::quantity& tstart, + const units::quantity& dt, + const units::quantity& tstop=terminal_time*units::ms) { return event_generator(std::move(target), weight, regular_schedule(tstart, dt, tstop)); } -template -inline event_generator poisson_generator( - cell_local_label_type target, - float weight, - time_type tstart, - time_type rate_kHz, - const RNG& rng, - time_type tstop=terminal_time) -{ - return event_generator(std::move(target), weight, poisson_schedule(tstart, rate_kHz, rng, tstop)); +inline event_generator poisson_generator(cell_local_label_type target, + float weight, + const units::quantity& tstart, + const units::quantity& rate_kHz, + seed_type seed = default_seed, + const units::quantity& tstop=terminal_time*units::ms) { + // TODO(TH) handle seed + return event_generator(std::move(target), weight, poisson_schedule(tstart, rate_kHz, seed, tstop)); } @@ -137,10 +127,16 @@ inline event_generator poisson_generator( template inline event_generator explicit_generator(cell_local_label_type target, float weight, - const S& s) -{ + const S& s) { return event_generator(std::move(target), weight, explicit_schedule(s)); } +template inline +event_generator explicit_generator_from_milliseconds(cell_local_label_type target, + float weight, + const S& s) { + return event_generator(std::move(target), weight, explicit_schedule_from_milliseconds(s)); +} + } // namespace arb diff --git a/arbor/include/arbor/schedule.hpp b/arbor/include/arbor/schedule.hpp index da3921bc36..8856cceac2 100644 --- a/arbor/include/arbor/schedule.hpp +++ b/arbor/include/arbor/schedule.hpp @@ -1,190 +1,39 @@ #pragma once #include -#include #include -#include #include #include #include +#include #include #include #include #include #include +#include // Time schedules for probe–sampler associations. namespace arb { +using engine_type = std::mt19937_64; +using seed_type = std::remove_cv_t; + +constexpr static auto default_seed = engine_type::default_seed; + using time_event_span = std::pair; inline time_event_span as_time_event_span(const std::vector& v) { - return {v.data(), v.data()+v.size()}; + return {v.data(), v.data() + v.size()}; } -// Common schedules - -// Default schedule is empty. - -struct ARB_ARBOR_API empty_schedule { - void reset() {} - time_event_span events(time_type t0, time_type t1) { - static time_type no_time; - return {&no_time, &no_time}; - } -}; - -// Schedule at k·dt for integral k≥0 within the interval [t0, t1). -struct ARB_ARBOR_API regular_schedule_impl { - explicit regular_schedule_impl(time_type t0, time_type dt, time_type t1): - t0_(t0), t1_(t1), dt_(dt), oodt_(1./dt) - { - if (t0_<0) t0_ = 0; - }; - - void reset() {} - time_event_span events(time_type t0, time_type t1); - - ARB_SERDES_ENABLE(regular_schedule_impl, t0_, t1_, dt_, oodt_); - -private: - time_type t0_, t1_, dt_; - time_type oodt_; - - std::vector times_; -}; - -// Schedule at times given explicitly via a provided sorted sequence. -struct ARB_ARBOR_API explicit_schedule_impl { - explicit_schedule_impl(const explicit_schedule_impl&) = default; - explicit_schedule_impl(explicit_schedule_impl&&) = default; - - template - explicit explicit_schedule_impl(const Seq& seq): - start_index_(0) - { - using std::begin; - using std::end; - - times_.assign(begin(seq), end(seq)); - arb_assert(std::is_sorted(times_.begin(), times_.end())); - } - - void reset() { - start_index_ = 0; - } - - time_event_span events(time_type t0, time_type t1); - - ARB_SERDES_ENABLE(explicit_schedule_impl, start_index_, times_); - -private: - std::ptrdiff_t start_index_; - std::vector times_; -}; - -// Schedule at Poisson point process with rate 1/mean_dt, -// restricted to non-negative times. -template -class poisson_schedule_impl { -public: - poisson_schedule_impl(time_type tstart, time_type rate_kHz, const RandomNumberEngine& rng, time_type tstop): - tstart_(tstart), exp_(rate_kHz), rng_(rng), reset_state_(rng), next_(tstart), tstop_(tstop) - { - arb_assert(tstart_>=0); - arb_assert(tstart_ <= tstop_); - step(); - } - - void reset() { - rng_ = reset_state_; - next_ = tstart_; - step(); - } - - time_event_span events(time_type t0, time_type t1) { - // if we start after the maximal allowed time, we have nothing to do - if (t0 >= tstop_) return {}; - - // restrict by maximal allowed time - t1 = std::min(t1, tstop_); - - times_.clear(); - - while (next_ exp_; - RandomNumberEngine rng_; - RandomNumberEngine reset_state_; - time_type next_; - std::vector times_; - time_type tstop_; -}; - // Type erased wrapper // A schedule describes a sequence of time values used for sampling. Schedules // are queried monotonically in time: if two method calls `events(t0, t1)` // and `events(t2, t3)` are made without an intervening call to `reset()`, // then 0 ≤ _t0_ ≤ _t1_ ≤ _t2_ ≤ _t3_. - -template -void serialize(arb::serializer& s, const K& k, const arb::explicit_schedule_impl&); -template -void deserialize(arb::serializer& s, const K& k, arb::explicit_schedule_impl&); - -template -void serialize(arb::serializer& s, const K& k, const arb::regular_schedule_impl&); -template -void deserialize(arb::serializer& s, const K& k, arb::regular_schedule_impl&); - -template -ARB_ARBOR_API void serialize(::arb::serializer& ser, - const K& k, - const arb::empty_schedule& t) { - ser.begin_write_map(::arb::to_serdes_key(k)); - ser.end_write_map(); -} - -template -ARB_ARBOR_API void deserialize(::arb::serializer& ser, - const K& k, - arb::empty_schedule& t) { - ser.begin_read_map(::arb::to_serdes_key(k)); - ser.end_read_map(); -} - -// These are custom to get the reset in. -template -void serialize(::arb::serializer& ser, const K& k, const ::arb::poisson_schedule_impl& t) { - ser.begin_write_map(arb::to_serdes_key(k)); - ARB_SERDES_WRITE(tstart_); - ARB_SERDES_WRITE(tstop_); - ser.end_write_map(); -} - -template -void deserialize(::arb::serializer& ser, const K& k, ::arb::poisson_schedule_impl& t) { - ser.begin_read_map(arb::to_serdes_key(k)); - ARB_SERDES_READ(tstart_); - ARB_SERDES_READ(tstop_); - ser.end_read_map(); - t.reset(); -} - -class schedule { -public: +struct ARB_ARBOR_API schedule { schedule(); template , schedule>>> @@ -206,12 +55,13 @@ class schedule { return *this; } - time_event_span events(time_type t0, time_type t1) { - return impl_->events(t0, t1); - } + time_event_span events(time_type t0, time_type t1) { return impl_->events(t0, t1); } void reset() { impl_->reset(); } + // Discard the next n events. Used in tests. + auto discard(std::size_t n) { return impl_->discard(n); } + template friend ARB_ARBOR_API void serialize(serializer& s, const K& k, const schedule& v) { v.impl_->t_serialize(s, to_serdes_key(k)); } template @@ -219,10 +69,10 @@ class schedule { )); } private: - struct interface { virtual time_event_span events(time_type t0, time_type t1) = 0; virtual void reset() = 0; + virtual void discard(std::size_t n) = 0; virtual std::unique_ptr clone() = 0; virtual ~interface() {} virtual void t_serialize(serializer&, const std::string&k) const = 0; @@ -237,46 +87,37 @@ class schedule { struct wrap: interface { explicit wrap(const Impl& impl): wrapped(impl) {} explicit wrap(Impl&& impl): wrapped(std::move(impl)) {} - virtual time_event_span events(time_type t0, time_type t1) override { return wrapped.events(t0, t1); } - virtual void reset() override { wrapped.reset(); } - virtual iface_ptr clone() override { return std::make_unique>(wrapped); } - virtual void t_serialize(serializer& s, const std::string& k) const override { serialize(s, k, wrapped); } - virtual void t_deserialize(serializer& s, const std::string& k) override { deserialize(s, k, wrapped); } + time_event_span events(time_type t0, time_type t1) override { return wrapped.events(t0, t1); } + void reset() override { wrapped.reset(); } + void discard(std::size_t n) override { return wrapped.discard(n); } + iface_ptr clone() override { return std::make_unique>(wrapped); } + void t_serialize(serializer& s, const std::string& k) const override { wrapped.t_serialize(s, k); } + void t_deserialize(serializer& s, const std::string& k) override { wrapped.t_deserialize(s, k); } Impl wrapped; }; }; // Constructors -inline schedule::schedule(): schedule(empty_schedule{}) {} -inline schedule regular_schedule(time_type t0, - time_type dt, - time_type t1 = std::numeric_limits::max()) { - return schedule(regular_schedule_impl(t0, dt, t1)); -} +/// Regular schedule with start `t0`, interval `dt`, and optional end `t1`. +schedule ARB_ARBOR_API regular_schedule(const units::quantity& t0, + const units::quantity& dt, + const units::quantity& t1 = std::numeric_limits::max()*units::ms); -inline schedule regular_schedule(time_type dt) { - return regular_schedule(0, dt); -} +/// Regular schedule with interval `dt`. +schedule ARB_ARBOR_API regular_schedule(const units::quantity& dt); -template -inline schedule explicit_schedule(const Seq& seq) { - return schedule(explicit_schedule_impl(seq)); -} +schedule ARB_ARBOR_API explicit_schedule(const std::vector& seq); +schedule ARB_ARBOR_API explicit_schedule_from_milliseconds(const std::vector& seq); -inline schedule explicit_schedule(const std::initializer_list& seq) { - return schedule(explicit_schedule_impl(seq)); -} +schedule ARB_ARBOR_API poisson_schedule(const units::quantity& tstart, + const units::quantity& rate, + seed_type seed = default_seed, + const units::quantity& tstop=terminal_time*units::ms); -template -inline schedule poisson_schedule(time_type rate_kHz, const RandomNumberEngine& rng, time_type tstop=terminal_time) { - return schedule(poisson_schedule_impl(0., rate_kHz, rng, tstop)); -} - -template -inline schedule poisson_schedule(time_type tstart, time_type rate_kHz, const RandomNumberEngine& rng, time_type tstop=terminal_time) { - return schedule(poisson_schedule_impl(tstart, rate_kHz, rng, tstop)); -} +schedule ARB_ARBOR_API poisson_schedule(const units::quantity& rate, + seed_type seed = default_seed, + const units::quantity& tstop=terminal_time*units::ms); } // namespace arb diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp index c6cf845d9c..b5de7b4a5c 100644 --- a/arbor/include/arbor/simulation.hpp +++ b/arbor/include/arbor/simulation.hpp @@ -1,8 +1,6 @@ #pragma once -#include #include -#include #include #include @@ -50,7 +48,7 @@ class ARB_ARBOR_API simulation { void reset(); - time_type run(time_type tfinal, time_type dt); + time_type run(const units::quantity& tfinal, const units::quantity& dt); // Minimum delay in network τ time_type min_delay(); diff --git a/arbor/include/arbor/units.hpp b/arbor/include/arbor/units.hpp index e837324b9e..659aa90388 100644 --- a/arbor/include/arbor/units.hpp +++ b/arbor/include/arbor/units.hpp @@ -25,6 +25,7 @@ constexpr auto Kelvin = ::units::Kelvin; using ::units::s; using ::units::ms; +using ::units::ns; using ::units::m; constexpr auto cm = centi * m; @@ -63,4 +64,5 @@ using ::units::mol; constexpr auto M = mol / m.pow(3); constexpr auto mM = milli * M; +using ::units::is_valid; } diff --git a/arbor/schedule.cpp b/arbor/schedule.cpp index 2a0d6befd6..2de0b30241 100644 --- a/arbor/schedule.cpp +++ b/arbor/schedule.cpp @@ -1,53 +1,253 @@ #include -#include -#include -#include #include #include #include -// Implementations for specific schedules. - namespace arb { -// Regular schedule implementation. +// Schedule at Poisson point process with rate 1/mean_dt, +// restricted to non-negative times. +struct poisson_schedule_impl { + poisson_schedule_impl(time_type tstart, time_type rate_kHz, seed_type seed, time_type tstop): + tstart_(tstart), exp_(rate_kHz), rng_(seed), seed_(seed), next_(tstart), tstop_(tstop) { + arb_assert(tstart_>=0); + arb_assert(tstart_ <= tstop_); + step(); + } + + void reset() { + rng_ = engine_type{seed_}; + next_ = tstart_; + step(); + } + + void discard(std::size_t n) { rng_.discard(n); } -time_event_span regular_schedule_impl::events(time_type t0, time_type t1) { - times_.clear(); + time_event_span events(time_type t0, time_type t1) { + // if we start after the maximal allowed time, we have nothing to do + if (t0 >= tstop_) return {}; - t0 = std::max(t0, t0_); - t1 = std::min(t1, t1_); + // restrict by maximal allowed time + t1 = std::min(t1, tstop_); - if (t1>t0) { - times_.reserve(1+std::size_t((t1-t0)*oodt_)); + times_.clear(); - long long n = t0*oodt_; - time_type t = n*dt_; + while (next_ + void t_serialize(::arb::serializer& ser, const K& k) const { + const auto& t = *this; + ser.begin_write_map(arb::to_serdes_key(k)); + ARB_SERDES_WRITE(tstart_); + ARB_SERDES_WRITE(tstop_); + ser.end_write_map(); + } + + template + void t_deserialize(::arb::serializer& ser, const K& k) { + auto& t = *this; + ser.begin_read_map(arb::to_serdes_key(k)); + ARB_SERDES_READ(tstart_); + ARB_SERDES_READ(tstop_); + ser.end_read_map(); + t.reset(); + } + + time_type tstart_; + std::exponential_distribution exp_; + engine_type rng_; + seed_type seed_; + time_type next_; + std::vector times_; + time_type tstop_; +}; + +schedule poisson_schedule(const units::quantity& tstart, + const units::quantity& rate, + seed_type seed, + const units::quantity& tstop) { + return schedule(poisson_schedule_impl(tstart.value_as(units::ms), + rate.value_as(units::kHz), + seed, + tstop.value_as(units::ms))); +} + +schedule poisson_schedule(const units::quantity& rate, + seed_type seed, + const units::quantity& tstop) { + return poisson_schedule(0.*units::ms, rate, seed, tstop); +} + + +struct empty_schedule { + void reset() {} + time_event_span events(time_type t0, time_type t1) { + static time_type no_time; + return {&no_time, &no_time}; + } + + + void t_serialize(::arb::serializer& ser, + const std::string& k) const { + ser.begin_write_map(::arb::to_serdes_key(k)); + ser.end_write_map(); + } + + void t_deserialize(::arb::serializer& ser, + const std::string& k) { + ser.begin_read_map(::arb::to_serdes_key(k)); + ser.end_read_map(); + } + + void discard(std::size_t) {} +}; + +schedule::schedule(): schedule(empty_schedule{}) {} + +// Schedule at k·dt for integral k≥0 within the interval [t0, t1). +struct ARB_ARBOR_API regular_schedule_impl { + explicit regular_schedule_impl(time_type t0, time_type dt, time_type t1): + t0_(t0), t1_(t1), dt_(dt), oodt_(1./dt) + { + if (t0_<0) t0_ = 0; + }; + + void reset() {} + + ARB_SERDES_ENABLE(regular_schedule_impl, t0_, t1_, dt_, oodt_); + + template + void t_serialize(::arb::serializer& ser, const K& k) const { + const auto& t = *this; + ser.begin_write_map(arb::to_serdes_key(k)); + ARB_SERDES_WRITE(t0_); + ARB_SERDES_WRITE(t1_); + ARB_SERDES_WRITE(dt_); + ser.end_write_map(); + } + + template + void t_deserialize(::arb::serializer& ser, const K& k) { + auto& t = *this; + ser.begin_read_map(arb::to_serdes_key(k)); + ARB_SERDES_READ(t0_); + ARB_SERDES_READ(t1_); + ARB_SERDES_READ(dt_); + oodt_ = 1.0/dt_; + ser.end_read_map(); + } + + time_event_span events(time_type t0, time_type t1) { + times_.clear(); + + t0 = std::max(t0, t0_); + t1 = std::min(t1, t1_); + + if (t1>t0) { + times_.reserve(1+std::size_t((t1-t0)*oodt_)); + + long long n = t0*oodt_; + time_type t = n*dt_; + + while (t times_; +}; + +schedule regular_schedule(const units::quantity& t0, + const units::quantity& dt, + const units::quantity& t1) { + return schedule(regular_schedule_impl(t0.value_as(units::ms), + dt.value_as(units::ms), + t1.value_as(units::ms))); } -// Explicit schedule implementation. +schedule regular_schedule(const units::quantity& dt) { + return regular_schedule(0*units::ms, dt); +} -time_event_span explicit_schedule_impl::events(time_type t0, time_type t1) { - time_event_span view = as_time_event_span(times_); +// Schedule at times given explicitly via a provided sorted sequence. +struct explicit_schedule_impl { + explicit_schedule_impl(const explicit_schedule_impl&) = default; + explicit_schedule_impl(explicit_schedule_impl&&) = default; - const time_type* lb = std::lower_bound(view.first+start_index_, view.second, t0); - const time_type* ub = std::lower_bound(lb, view.second, t1); + explicit explicit_schedule_impl(std::vector seq): + start_index_(0), times_(std::move(seq)) { + arb_assert(std::is_sorted(times_.begin(), times_.end())); + } - start_index_ = ub-view.first; - return {lb, ub}; + void reset() { start_index_ = 0; } + + time_event_span events(time_type t0, time_type t1) { + time_event_span view = as_time_event_span(times_); + + const time_type* lb = std::lower_bound(view.first+start_index_, view.second, t0); + const time_type* ub = std::lower_bound(lb, view.second, t1); + + start_index_ = ub-view.first; + return {lb, ub}; + } + + template + void t_serialize(::arb::serializer& ser, const K& k) const { + const auto& t = *this; + ser.begin_write_map(arb::to_serdes_key(k)); + ARB_SERDES_WRITE(start_index_); + ARB_SERDES_WRITE(times_); + ser.end_write_map(); + } + + template + void t_deserialize(::arb::serializer& ser, const K& k) { + auto& t = *this; + ser.begin_read_map(arb::to_serdes_key(k)); + ARB_SERDES_READ(start_index_); + ARB_SERDES_READ(times_); + ser.end_read_map(); + } + + void discard(std::size_t) {} + + std::ptrdiff_t start_index_; + std::vector times_; +}; + +schedule explicit_schedule(const std::vector& seq) { + std::vector res; + for (const auto& t: seq) res.push_back(t.value_as(units::ms)); + return schedule(explicit_schedule_impl(res)); } +schedule explicit_schedule_from_milliseconds(const std::vector& seq) { + return schedule(explicit_schedule_impl(seq)); +} + + } // namespace arb diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index dac25d6fef..2ceda88c21 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -1,5 +1,4 @@ #include -#include #include #include @@ -15,13 +14,10 @@ #include "cell_group.hpp" #include "cell_group_factory.hpp" #include "communication/communicator.hpp" -#include "execution_context.hpp" #include "merge_events.hpp" #include "thread_private_spike_store.hpp" #include "threading/threading.hpp" -#include "util/filter.hpp" #include "util/maputil.hpp" -#include "util/partition.hpp" #include "util/span.hpp" #include "profile/profiler_macro.hpp" @@ -596,11 +592,11 @@ void simulation::reset() { void simulation::update(const connectivity& rec) { impl_->update(rec); } -time_type simulation::run(time_type tfinal, time_type dt) { - if (dt <= 0.0) { +time_type simulation::run(const units::quantity& tfinal, const units::quantity& dt) { + if (dt.value() <= 0.0) { throw domain_error("Finite time-step must be supplied."); } - return impl_->run(tfinal, dt); + return impl_->run(tfinal.value_as(units::ms), dt.value_as(units::ms)); } sampler_association_handle simulation::add_sampler( diff --git a/python/schedule.cpp b/python/schedule.cpp index 70950a4a60..e9ec42a1dc 100644 --- a/python/schedule.cpp +++ b/python/schedule.cpp @@ -13,57 +13,70 @@ namespace py = pybind11; namespace pyarb { std::ostream& operator<<(std::ostream& o, const regular_schedule_shim& x) { - return o << ""; + if (x.tstop.has_value()) { + return o << ""; + } + else { + return o << ""; + } } std::ostream& operator<<(std::ostream& o, const explicit_schedule_shim& e) { - o << ""; -}; + return o << ""; +} std::ostream& operator<<(std::ostream& o, const poisson_schedule_shim& p) { - return o << ""; -}; + if (p.tstop) { + return o << ""; + } + else { + return o << ""; + } +} static std::vector as_vector(std::pair ts) { return std::vector(ts.first, ts.second); } -// -// regular_schedule shim -// - -regular_schedule_shim::regular_schedule_shim(arb::time_type t0, arb::time_type delta_t, py::object t1) { +regular_schedule_shim::regular_schedule_shim(const arb::units::quantity& t0, + const arb::units::quantity& delta_t, + std::optional t1) { set_tstart(t0); set_dt(delta_t); set_tstop(t1); } -regular_schedule_shim::regular_schedule_shim(arb::time_type delta_t) { - set_tstart(0.); +regular_schedule_shim::regular_schedule_shim(const arb::units::quantity& delta_t) { + set_tstart(0.*arb::units::ms); set_dt(delta_t); } -void regular_schedule_shim::set_tstart(arb::time_type t) { - pyarb::assert_throw(is_nonneg()(t), "tstart must be a non-negative number"); +void regular_schedule_shim::set_tstart(const arb::units::quantity& t) { + pyarb::assert_throw(is_nonneg()(t.value()), "tstart must be a non-negative number"); + pyarb::assert_throw(arb::units::is_valid(t.convert_to(arb::units::ms)), "must be convertible to time"); tstart = t; -}; +} -void regular_schedule_shim::set_tstop(py::object t) { - tstop = py2optional( - t, "tstop must be a non-negative number, or None", is_nonneg()); -}; +void regular_schedule_shim::set_tstop(std::optional t) { + if (t.has_value()) { + pyarb::assert_throw(arb::units::is_valid(t.value().convert_to(arb::units::ms)), "must be convertible to time"); + } + tstop = t; +} -void regular_schedule_shim::set_dt(arb::time_type delta_t) { - pyarb::assert_throw(is_positive()(delta_t), "dt must be a positive number"); - dt = delta_t; -}; +void regular_schedule_shim::set_dt(const arb::units::quantity& t) { + pyarb::assert_throw(is_positive()(t.value()), "dt must be a positive number"); + pyarb::assert_throw(arb::units::is_valid(t.convert_to(arb::units::ms)), "must be convertible to time"); + dt = t; +} regular_schedule_shim::time_type regular_schedule_shim::get_tstart() const { return tstart; @@ -78,10 +91,9 @@ regular_schedule_shim::opt_time_type regular_schedule_shim::get_tstop() const { } arb::schedule regular_schedule_shim::schedule() const { - return arb::regular_schedule( - tstart, - dt, - tstop.value_or(arb::terminal_time)); + return arb::regular_schedule(tstart, + dt, + tstop.value_or(arb::terminal_time*arb::units::ms)); } std::vector regular_schedule_shim::events(arb::time_type t0, arb::time_type t1) { @@ -93,17 +105,14 @@ std::vector regular_schedule_shim::events(arb::time_type t0, arb return as_vector(sched.events(t0, t1)); } -// -// explicit_schedule shim -// - -//struct explicit_schedule_shim { -explicit_schedule_shim::explicit_schedule_shim(std::vector t) { - set_times(t); +explicit_schedule_shim::explicit_schedule_shim(const std::vector& seq) { + std::vector ts; + for (const auto t: seq) ts.push_back(t.value_as(arb::units::ms)); + set_times_ms(std::move(ts)); } // getter and setter (in order to assert when being set) -void explicit_schedule_shim::set_times(std::vector t) { +void explicit_schedule_shim::set_times_ms(std::vector t) { times = std::move(t); // Sort the times in ascending order if necessary @@ -118,12 +127,12 @@ void explicit_schedule_shim::set_times(std::vector t) { } }; -std::vector explicit_schedule_shim::get_times() const { +std::vector explicit_schedule_shim::get_times_ms() const { return times; } arb::schedule explicit_schedule_shim::schedule() const { - return arb::explicit_schedule(times); + return arb::explicit_schedule_from_milliseconds(times); } std::vector explicit_schedule_shim::events(arb::time_type t0, arb::time_type t1) { @@ -135,71 +144,56 @@ std::vector explicit_schedule_shim::events(arb::time_type t0, ar return as_vector(sched.events(t0, t1)); } -// -// poisson_schedule shim -// - -poisson_schedule_shim::poisson_schedule_shim( - arb::time_type ts, - arb::time_type f, - rng_type::result_type s, - py::object tstop) -{ +poisson_schedule_shim::poisson_schedule_shim(const arb::units::quantity& ts, + const arb::units::quantity& f, + arb::seed_type s, + std::optional tstop) { set_tstart(ts); set_freq(f); seed = s; set_tstop(tstop); } -poisson_schedule_shim::poisson_schedule_shim(arb::time_type f) { - set_tstart(0.); +poisson_schedule_shim::poisson_schedule_shim(const arb::units::quantity& f) { + set_tstart(0.*arb::units::ms); set_freq(f); - seed = 0; } -void poisson_schedule_shim::set_tstart(arb::time_type t) { - pyarb::assert_throw(is_nonneg()(t), "tstart must be a non-negative number"); +void poisson_schedule_shim::set_tstart(const arb::units::quantity& t) { + pyarb::assert_throw(is_nonneg()(t.value()), "tstart must be a non-negative number"); tstart = t; }; -void poisson_schedule_shim::set_freq(arb::time_type f) { - pyarb::assert_throw(is_nonneg()(f), "frequency must be a non-negative number"); +void poisson_schedule_shim::set_freq(const arb::units::quantity& f) { + pyarb::assert_throw(is_nonneg()(f.value()), "frequency must be a non-negative number"); freq = f; }; -void poisson_schedule_shim::set_tstop(py::object t) { - tstop = py2optional( - t, "tstop must be a non-negative number, or None", is_nonneg()); +void poisson_schedule_shim::set_tstop(std::optional t) { + // TODO(TH) + // if (t.has_value()) pyarb::assert_throw(is_nonneg()(t.value()), "frequency must be a non-negative number"); + tstop = t; }; -arb::time_type poisson_schedule_shim::get_tstart() const { - return tstart; -} - -arb::time_type poisson_schedule_shim::get_freq() const { - return freq; -} - -poisson_schedule_shim::opt_time_type poisson_schedule_shim::get_tstop() const { - return tstop; -} - arb::schedule poisson_schedule_shim::schedule() const { - return arb::poisson_schedule(tstart, freq, rng_type(seed), tstop.value_or(arb::terminal_time)); + return arb::poisson_schedule(tstart, freq, seed, tstop.value_or(arb::terminal_time*arb::units::ms)); } -std::vector poisson_schedule_shim::events(arb::time_type t0, arb::time_type t1) { - pyarb::assert_throw(is_nonneg()(t0), "t0 must be a non-negative number"); - pyarb::assert_throw(is_nonneg()(t1), "t1 must be a non-negative number"); +std::vector poisson_schedule_shim::events(const arb::units::quantity& t0, + const arb::units::quantity& t1) { + auto beg = t0.value_as(arb::units::ms); + auto end = t1.value_as(arb::units::ms); + pyarb::assert_throw(is_nonneg()(beg), "t0 must be a non-negative number"); + pyarb::assert_throw(is_nonneg()(end), "t1 must be a non-negative number"); arb::schedule sched = poisson_schedule_shim::schedule(); - return as_vector(sched.events(t0, t1)); + return as_vector(sched.events(beg, end)); } void register_schedules(py::module& m) { using namespace py::literals; - using time_type = arb::time_type; + using time_type = arb::units::quantity; py::class_ schedule_base(m, "schedule_base", "Schedule abstract base class."); @@ -208,13 +202,13 @@ void register_schedules(py::module& m) { "Describes a regular schedule with multiples of dt within the interval [tstart, tstop)."); regular_schedule - .def(py::init(), + .def(py::init>(), "tstart"_a, "dt"_a, "tstop"_a = py::none(), "Construct a regular schedule with arguments:\n" " tstart: The delivery time of the first event in the sequence [ms].\n" " dt: The interval between time points [ms].\n" " tstop: No events delivered after this time [ms], None by default.") - .def(py::init(), + .def(py::init(), "dt"_a, "Construct a regular schedule, starting from t = 0 and never terminating, with arguments:\n" " dt: The interval between time points [ms].\n") @@ -240,10 +234,10 @@ void register_schedules(py::module& m) { "times"_a, "Construct an explicit schedule with argument:\n" " times: A list of times [ms], [] by default.") - .def_property("times", &explicit_schedule_shim::get_times, &explicit_schedule_shim::set_times, + .def_property("times_ms", &explicit_schedule_shim::get_times_ms, &explicit_schedule_shim::set_times_ms, "A list of times [ms].") .def("events", &explicit_schedule_shim::events, - "A view of monotonically increasing time values in the half-open interval [t0, t1).") + "A view of monotonically increasing time values in the half-open interval [t0, t1) in [ms].") .def("__str__", util::to_string) .def("__repr__", util::to_string); @@ -252,8 +246,8 @@ void register_schedules(py::module& m) { "Describes a schedule according to a Poisson process within the interval [tstart, tstop)."); poisson_schedule - .def(py::init(), - "tstart"_a = 0., "freq"_a, "seed"_a = 0, "tstop"_a = py::none(), + .def(py::init>(), + "tstart"_a = 0.*arb::units::ms, "freq"_a, "seed"_a = 0, "tstop"_a, "Construct a Poisson schedule with arguments:\n" " tstart: The delivery time of the first event in the sequence [ms], 0 by default.\n" " freq: The expected frequency [kHz].\n" diff --git a/python/schedule.hpp b/python/schedule.hpp index 619f77236a..7fade42297 100644 --- a/python/schedule.hpp +++ b/python/schedule.hpp @@ -28,20 +28,23 @@ struct schedule_shim_base { // a regular_schedule in python are manipulating this type. This is converted to // an arb::regular_schedule when a C++ recipe is created from a Python recipe. struct regular_schedule_shim: schedule_shim_base { - using time_type = arb::time_type; + using time_type = arb::units::quantity; using opt_time_type = std::optional; - time_type tstart = {}; - time_type dt = 0; - opt_time_type tstop = {}; + time_type tstart; + time_type dt = 0*arb::units::ms; + opt_time_type tstop; - regular_schedule_shim(time_type t0, time_type delta_t, pybind11::object t1); - explicit regular_schedule_shim(time_type delta_t); + regular_schedule_shim(const time_type& t0, + const time_type& delta_t, + opt_time_type t1); + + explicit regular_schedule_shim(const time_type& delta_t); // getter and setter (in order to assert when being set) - void set_tstart(time_type t); - void set_dt(time_type delta_t); - void set_tstop(pybind11::object t); + void set_tstart(const time_type& t); + void set_dt(const time_type& delta_t); + void set_tstop(opt_time_type t); time_type get_tstart() const; time_type get_dt() const; @@ -60,11 +63,11 @@ struct explicit_schedule_shim: schedule_shim_base { std::vector times; explicit_schedule_shim() = default; - explicit_schedule_shim(std::vector t); + explicit_schedule_shim(const std::vector& t); // getter and setter (in order to assert when being set) - void set_times(std::vector t); - std::vector get_times() const; + void set_times_ms(std::vector t); + std::vector get_times_ms() const; arb::schedule schedule() const override; @@ -76,28 +79,30 @@ struct explicit_schedule_shim: schedule_shim_base { // Python are manipulating this type. This is converted to an // arb::poisson_schedule when a C++ recipe is created from a Python recipe. struct poisson_schedule_shim: schedule_shim_base { - using rng_type = std::mt19937_64; - using opt_time_type = std::optional; + arb::units::quantity tstart; // ms + arb::units::quantity freq; // kHz + std::optional tstop; // ms + arb::seed_type seed = arb::default_seed; - arb::time_type tstart; // ms - arb::time_type freq; // kHz - opt_time_type tstop; // ms - rng_type::result_type seed; + poisson_schedule_shim(const arb::units::quantity& ts, + const arb::units::quantity& f, + arb::seed_type s, + std::optional tstop); + poisson_schedule_shim(const arb::units::quantity& f); - poisson_schedule_shim(arb::time_type ts, arb::time_type f, rng_type::result_type s, pybind11::object tstop); - poisson_schedule_shim(arb::time_type f); + void set_tstart(const arb::units::quantity& t); + void set_freq(const arb::units::quantity& f); + void set_tstop(std::optional f); - void set_tstart(arb::time_type t); - void set_freq(arb::time_type f); - void set_tstop(pybind11::object t); + const auto& get_tstop() const { return tstop; } + const auto& get_tstart() const { return tstart; } + const auto& get_freq() const { return freq; } - arb::time_type get_tstart() const; - arb::time_type get_freq() const; - opt_time_type get_tstop() const; arb::schedule schedule() const override; - std::vector events(arb::time_type t0, arb::time_type t1); + // TODO(TH) this should be symmetrical... + std::vector events(const arb::units::quantity& t0, const arb::units::quantity& t1); }; } diff --git a/python/simulation.cpp b/python/simulation.cpp index 59ee2747ba..5f7edadd4b 100644 --- a/python/simulation.cpp +++ b/python/simulation.cpp @@ -116,7 +116,7 @@ class simulation_shim { } } - arb::time_type run(arb::time_type tfinal, arb::time_type dt) { + arb::time_type run(const arb::units::quantity& tfinal, const arb::units::quantity& dt) { return sim_->run(tfinal, dt); } diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp index 020a637373..2b3af09930 100644 --- a/python/single_cell_model.cpp +++ b/python/single_cell_model.cpp @@ -34,7 +34,7 @@ namespace pyarb { // Stores the location and sampling frequency for a probe in a single cell model. struct probe_site { arb::locset locset; // Location of sample on morphology. - double frequency; // Sampling frequency [kHz]. + arb::units::quantity frequency; // Sampling frequency [kHz]. arb::cell_tag_type tag; // Tag = unique name }; @@ -161,12 +161,15 @@ class single_cell_model { // m.probe('voltage', '(location 2 0.5)') // m.probe('voltage', 'term') - void probe(const std::string& what, const arb::locset& where, const arb::cell_tag_type& tag, double frequency) { + void probe(const std::string& what, + const arb::locset& where, + const arb::cell_tag_type& tag, + const arb::units::quantity& frequency) { if (what != "voltage") { throw pyarb_error( util::pprintf("{} does not name a valid variable to trace (currently only 'voltage' is supported)", what)); } - if (frequency<=0) { + if (frequency.value() <= 0) { throw pyarb_error( util::pprintf("sampling frequency is not greater than zero", what)); } @@ -177,7 +180,7 @@ class single_cell_model { event_generators_.push_back(event_generator); } - void run(double tfinal, double dt) { + void run(const arb::units::quantity& tfinal, const arb::units::quantity& dt) { single_cell_recipe rec(cell_, probes_, gprop, event_generators_); auto domdec = arb::partition_load_balance(rec, ctx_); @@ -255,7 +258,7 @@ void register_single_cell(pybind11::module& m) { "dt"_a = 0.025, "Run model from t=0 to t=tfinal ms.") .def("probe", - [](single_cell_model& m, const char* what, const char* where, const char* tag, double frequency) { + [](single_cell_model& m, const char* what, const char* where, const char* tag, const arb::units::quantity& frequency) { m.probe(what, arborio::parse_locset_expression(where).unwrap(), tag, frequency);}, "what"_a, "where"_a, "tag"_a, "frequency"_a, "Sample a variable on the cell.\n" @@ -264,7 +267,7 @@ void register_single_cell(pybind11::module& m) { " tag: Unique name for this probe.\n" " frequency: The target frequency at which to sample [kHz].") .def("probe", - [](single_cell_model& m, const char* what, const arb::mlocation& where, const char* tag, double frequency) { + [](single_cell_model& m, const char* what, const arb::mlocation& where, const char* tag, const arb::units::quantity& frequency) { m.probe(what, where, tag, frequency);}, "what"_a, "where"_a, "tag"_a, "frequency"_a, "Sample a variable on the cell.\n" diff --git a/test/common_cells.cpp b/test/common_cells.cpp index 5d1f3aa2bc..903e4e0a43 100644 --- a/test/common_cells.cpp +++ b/test/common_cells.cpp @@ -1,5 +1,4 @@ #include -#include "arbor/morph/morphology.hpp" #include "common_cells.hpp" namespace arb { @@ -181,7 +180,9 @@ cable_cell_description make_cell_soma_only(bool with_stim) { auto c = builder.make_cell(); c.decorations.paint("soma"_lab, density("hh")); if (with_stim) { - c.decorations.place(builder.location({0,0.5}), i_clamp{10., 100., 0.1}, "cc"); + c.decorations.place(builder.location({0,0.5}), + i_clamp::box(10.*arb::units::ms, 100.*arb::units::ms, 0.1*arb::units::nA), + "cc"); } return {c.morph, c.labels, c.decorations}; @@ -216,7 +217,9 @@ cable_cell_description make_cell_ball_and_stick(bool with_stim) { c.decorations.paint("soma"_lab, density("hh")); c.decorations.paint("dend"_lab, density("pas")); if (with_stim) { - c.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3}, "cc"); + c.decorations.place(builder.location({1,1}), + i_clamp::box(5*arb::units::ms, 80*arb::units::ms, 0.3*arb::units::nA), + "cc"); } return {c.morph, c.labels, c.decorations}; @@ -254,8 +257,12 @@ cable_cell_description make_cell_ball_and_3stick(bool with_stim) { c.decorations.paint("soma"_lab, density("hh")); c.decorations.paint("dend"_lab, density("pas")); if (with_stim) { - c.decorations.place(builder.location({2,1}), i_clamp{5., 80., 0.45}, "cc0"); - c.decorations.place(builder.location({3,1}), i_clamp{40., 10.,-0.2}, "cc1"); + c.decorations.place(builder.location({2,1}), + i_clamp::box(5.*arb::units::ms, 80.*arb::units::ms, 0.45*arb::units::nA), + "cc0"); + c.decorations.place(builder.location({3,1}), + i_clamp::box(40.*arb::units::ms, 10.*arb::units::ms,-0.2*arb::units::nA), + "cc1"); } return {c.morph, c.labels, c.decorations}; diff --git a/test/common_cells.hpp b/test/common_cells.hpp index 9b22af70e7..6e9cda85a7 100644 --- a/test/common_cells.hpp +++ b/test/common_cells.hpp @@ -1,5 +1,3 @@ -#include - #include #include #include diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp index 9f06a3dda9..80b3274476 100644 --- a/test/unit-distributed/test_communicator.cpp +++ b/test/unit-distributed/test_communicator.cpp @@ -1,7 +1,6 @@ #include #include "test.hpp" -#include #include #include @@ -11,7 +10,6 @@ #include #include "communication/communicator.hpp" -#include "execution_context.hpp" #include "fvm_lowered_cell.hpp" #include "lif_cell_group.hpp" #include "cable_cell_group.hpp" @@ -204,7 +202,7 @@ namespace { tree.append(arb::mnpos, {0, 0, 0.0, 1.0}, {0, 0, 200, 1.0}, 1); arb::decor decor; decor.set_default(arb::cv_policy_fixed_per_branch(10)); - decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10}, "src"); + decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10*arb::units::mV}, "src"); decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "tgt"); return arb::cable_cell(arb::morphology(tree), decor); } @@ -277,7 +275,7 @@ namespace { tree.append(arb::mnpos, {0, 0, 0.0, 1.0}, {0, 0, 200, 1.0}, 1); arb::decor decor; decor.set_default(arb::cv_policy_fixed_per_branch(10)); - decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10}, "src"); + decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10*arb::units::mV}, "src"); decor.place(arb::ls::uniform(arb::reg::all(), 0, size_, gid), arb::synapse("expsyn"), "tgt"); return arb::cable_cell(arb::morphology(tree), decor); } @@ -354,8 +352,8 @@ namespace { decor.place(arb::ls::uniform(arb::reg::all(), 2, 2, gid), arb::synapse("expsyn"), "synapses_1"); } else { - decor.place(arb::ls::uniform(arb::reg::all(), 0, 2, gid), arb::threshold_detector{10}, "detectors_0"); - decor.place(arb::ls::uniform(arb::reg::all(), 3, 3, gid), arb::threshold_detector{10}, "detectors_1"); + decor.place(arb::ls::uniform(arb::reg::all(), 0, 2, gid), arb::threshold_detector{10*arb::units::mV}, "detectors_0"); + decor.place(arb::ls::uniform(arb::reg::all(), 3, 3, gid), arb::threshold_detector{10*arb::units::mV}, "detectors_1"); } return arb::cable_cell(arb::morphology(tree), decor); } diff --git a/test/unit/test_cable_cell.cpp b/test/unit/test_cable_cell.cpp index 1e59691bd2..15042c88a9 100644 --- a/test/unit/test_cable_cell.cpp +++ b/test/unit/test_cable_cell.cpp @@ -2,6 +2,7 @@ #include "../common_cells.hpp" #include +#include #include #include @@ -33,9 +34,9 @@ TEST(cable_cell, lid_ranges) { // Note: there are 2 terminal points. decorations.place("term"_lab, synapse("expsyn"), "t0"); decorations.place("term"_lab, synapse("expsyn"), "t1"); - decorations.place("term"_lab, threshold_detector{-10}, "s0"); + decorations.place("term"_lab, threshold_detector{-10*arb::units::mV}, "s0"); decorations.place(empty_sites, synapse("expsyn"), "t2"); - decorations.place("term"_lab, threshold_detector{-20}, "s1"); + decorations.place("term"_lab, threshold_detector{-20*arb::units::mV}, "s1"); decorations.place(three_sites, synapse("expsyn"), "t3"); decorations.place("term"_lab, synapse("exp2syn"), "t3"); diff --git a/test/unit/test_cable_cell_group.cpp b/test/unit/test_cable_cell_group.cpp index f143df0a25..cbb3871454 100644 --- a/test/unit/test_cable_cell_group.cpp +++ b/test/unit/test_cable_cell_group.cpp @@ -7,7 +7,6 @@ #include "epoch.hpp" #include "fvm_lowered_cell.hpp" #include "cable_cell_group.hpp" -#include "util/rangeutil.hpp" #include "common.hpp" #include "../common_cells.hpp" @@ -29,8 +28,8 @@ namespace { auto d = builder.make_cell(); d.decorations.paint("soma"_lab, density("hh")); d.decorations.paint("dend"_lab, density("pas")); - d.decorations.place(builder.location({1,1}), i_clamp::box(5, 80, 0.3), "clamp0"); - d.decorations.place(builder.location({0, 0}), threshold_detector{0}, "detector0"); + d.decorations.place(builder.location({1,1}), i_clamp::box(5*arb::units::ms, 80*arb::units::ms, 0.3*arb::units::nA), "clamp0"); + d.decorations.place(builder.location({0, 0}), threshold_detector{0*arb::units::mV}, "detector0"); return d; } } @@ -72,7 +71,7 @@ TEST(cable_cell_group, sources) { for (int i=0; i<20; ++i) { auto desc = make_cell(); if (i==0 || i==3 || i==17) { - desc.decorations.place(mlocation{0, 0.3}, threshold_detector{2.3}, "detector1"); + desc.decorations.place(mlocation{0, 0.3}, threshold_detector{2.3*arb::units::mV}, "detector1"); } cells.emplace_back(desc); diff --git a/test/unit/test_diffusion.cpp b/test/unit/test_diffusion.cpp index 348b135ca8..7bcc173837 100644 --- a/test/unit/test_diffusion.cpp +++ b/test/unit/test_diffusion.cpp @@ -1,5 +1,4 @@ #include -#include #include #include @@ -55,7 +54,7 @@ struct linear: public recipe { std::vector event_generators(arb::cell_gid_type gid) const override { std::vector result; for (const auto& [t, w]: inject_at) { - result.push_back(arb::explicit_generator({"Zap"}, w, std::vector{t})); + result.push_back(arb::explicit_generator_from_milliseconds({"Zap"}, w, std::vector{t})); } return result; } @@ -117,8 +116,8 @@ testing::AssertionResult run(const linear& rec, const result_t exp) { }; auto ctx = make_context({arbenv::default_concurrency(), with_gpu}); auto sim = simulation{rec, ctx, partition_load_balance(rec, ctx)}; - sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1), sampler); - sim.run(0.11, 0.01); + sim.add_sampler(arb::all_probes, arb::regular_schedule(0.1*arb::units::ms), sampler); + sim.run(0.11*arb::units::ms, 0.01*arb::units::ms); return all_near(sample_values, exp, epsilon); } @@ -316,21 +315,21 @@ TEST(diffusion, setting_diffusivity) { { R r; r.gprop.add_ion("bla", 1, 23, 42, 0, 0); - EXPECT_THROW(simulation(r).run(1, 1), illegal_diffusive_mechanism); + EXPECT_THROW(simulation(r).run(1*arb::units::ms, 1*arb::units::ms), illegal_diffusive_mechanism); } // BAD: Trying to use a partially diffusive ion { R r; r.gprop.add_ion("bla", 1, 23, 42, 0, 0); r.dec.paint("(tag 1)"_reg, ion_diffusivity{"bla", 13}); - EXPECT_THROW(simulation(r).run(1, 1), cable_cell_error); + EXPECT_THROW(simulation(r).run(1*arb::units::ms, 1*arb::units::ms), cable_cell_error); } // OK: Using the global default { R r; r.gprop.add_ion("bla", 1, 23, 42, 0, 8); r.dec.paint("(tag 1)"_reg, ion_diffusivity{"bla", 13}); - EXPECT_NO_THROW(simulation(r).run(1, 1)); + EXPECT_NO_THROW(simulation(r).run(1*arb::units::ms, 1*arb::units::ms)); } // OK: Using the cell default { @@ -338,14 +337,14 @@ TEST(diffusion, setting_diffusivity) { r.gprop.add_ion("bla", 1, 23, 42, 0, 0); r.dec.set_default(ion_diffusivity{"bla", 8}); r.dec.paint("(tag 1)"_reg, ion_diffusivity{"bla", 13}); - EXPECT_NO_THROW(simulation(r).run(1, 1)); + EXPECT_NO_THROW(simulation(r).run(1*arb::units::ms, 1*arb::units::ms)); } // BAD: Using an unknown species { R r; r.dec.set_default(ion_diffusivity{"bla", 8}); r.dec.paint("(tag 1)"_reg, ion_diffusivity{"bla", 13}); - EXPECT_THROW(simulation(r).run(1, 1), cable_cell_error); + EXPECT_THROW(simulation(r).run(1*arb::units::ms, 1*arb::units::ms), cable_cell_error); } } diff --git a/test/unit/test_domain_decomposition.cpp b/test/unit/test_domain_decomposition.cpp index ab8a7c3e4a..10e0b754c6 100644 --- a/test/unit/test_domain_decomposition.cpp +++ b/test/unit/test_domain_decomposition.cpp @@ -1,7 +1,5 @@ #include -#include - #include #include #include diff --git a/test/unit/test_event_delivery.cpp b/test/unit/test_event_delivery.cpp index 89718a1958..8800224e24 100644 --- a/test/unit/test_event_delivery.cpp +++ b/test/unit/test_event_delivery.cpp @@ -37,7 +37,7 @@ struct test_recipe: public n_cable_cell_recipe { decor decorations; decorations.place(mlocation{0, 0.5}, synapse("expsyn"), "synapse"); - decorations.place(mlocation{0, 0.5}, threshold_detector{-64}, "detector"); + decorations.place(mlocation{0, 0.5}, threshold_detector{-64*arb::units::mV}, "detector"); decorations.place(mlocation{0, 0.5}, junction("gj"), "gapjunction"); cable_cell c(st, decorations, labels); @@ -74,7 +74,7 @@ std::vector run_test_sim(const recipe& R, const group_gids_type& } sim.inject_events(cell_events); - sim.run((n+1)*ev_delta_t, 0.01); + sim.run((n+1)*ev_delta_t*arb::units::ms, 0.01*arb::units::ms); std::vector spike_gids; util::sort_by(spikes, [](auto s) { return s.time; }); diff --git a/test/unit/test_event_generators.cpp b/test/unit/test_event_generators.cpp index 024fbc1a08..84d9082669 100644 --- a/test/unit/test_event_generators.cpp +++ b/test/unit/test_event_generators.cpp @@ -5,10 +5,6 @@ #include #include -#include "util/rangeutil.hpp" - -#include "common.hpp" - using namespace arb; namespace{ @@ -18,7 +14,7 @@ namespace{ } TEST(event_generators, assign_and_copy) { - event_generator gen = regular_generator({"l2"}, 5., 0.5, 0.75); + event_generator gen = regular_generator({"l2"}, 5., 0.5*arb::units::ms, 0.75*arb::units::ms); gen.resolve_label([](const cell_local_label_type&) {return 2;}); spike_event expected{2, 0.75, 5.}; @@ -57,7 +53,7 @@ TEST(event_generators, regular) { cell_lid_type lid = 3; float weight = 3.14; - event_generator gen = regular_generator(label, weight, t0, dt); + event_generator gen = regular_generator(label, weight, t0*arb::units::ms, dt*arb::units::ms); gen.resolve_label([lid](const cell_local_label_type&) {return lid;}); // Helper for building a set of expected events. @@ -93,7 +89,7 @@ TEST(event_generators, seq) { expected.push_back({0, time, weight}); } - event_generator gen = explicit_generator(l0, weight, times); + event_generator gen = explicit_generator_from_milliseconds(l0, weight, times); gen.resolve_label([](const cell_local_label_type&) {return 0;}); EXPECT_EQ(expected, as_vector(gen.events(0, 100.))); gen.reset(); @@ -126,8 +122,6 @@ TEST(event_generators, seq) { } TEST(event_generators, poisson) { - std::mt19937_64 G; - time_type t0 = 0; time_type t1 = 10; time_type lambda = 10; // expect 10 events per ms @@ -135,7 +129,7 @@ TEST(event_generators, poisson) { cell_lid_type lid = 2; float weight = 42; - event_generator gen = poisson_generator(label, weight, t0, lambda, G); + event_generator gen = poisson_generator(label, weight, t0*arb::units::ms, lambda*arb::units::Hz); gen.resolve_label([lid](const cell_local_label_type&) {return lid;}); pse_vector int1 = as_vector(gen.events(0, t1)); diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp index 5569171d7d..7175c99f6a 100644 --- a/test/unit/test_fvm_layout.cpp +++ b/test/unit/test_fvm_layout.cpp @@ -21,7 +21,6 @@ #include "util/maputil.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" -#include "io/sepval.hpp" #include "common.hpp" #include "common_morphologies.hpp" @@ -71,10 +70,12 @@ namespace { auto description = builder.make_cell(); description.decorations.paint("soma"_lab, density("hh")); description.decorations.paint("dend"_lab, density("pas")); - description.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3}, "clamp"); - + description.decorations.place(builder.location({1,1}), + i_clamp{5*arb::units::nA, 80*arb::units::kHz, 0.3*arb::units::rad}, + "clamp"); s.builders.push_back(std::move(builder)); descriptions.push_back(description); + } // Cell 1: ball and 3-stick, but with uneven dendrite @@ -123,8 +124,12 @@ namespace { desc.decorations.paint(c2, membrane_capacitance{0.013}); desc.decorations.paint(c3, membrane_capacitance{0.018}); - desc.decorations.place(b.location({2,1}), i_clamp{5., 80., 0.45}, "clamo0"); - desc.decorations.place(b.location({3,1}), i_clamp{40., 10.,-0.2}, "clamp1"); + desc.decorations.place(b.location({2,1}), + i_clamp::box( 5.*arb::units::ms, 80.*arb::units::ms, 0.45*arb::units::nA), + "clamp0"); + desc.decorations.place(b.location({3,1}), + i_clamp::box(40.*arb::units::ms, 10.*arb::units::ms, -0.2*arb::units::nA), + "clamp1"); desc.decorations.set_default(axial_resistivity{90}); diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp index 19b064f60e..8254f2db5c 100644 --- a/test/unit/test_fvm_lowered.cpp +++ b/test/unit/test_fvm_lowered.cpp @@ -1,5 +1,4 @@ #include -#include #include #include @@ -25,11 +24,8 @@ #include "backends/multicore/fvm.hpp" #include "fvm_lowered_cell.hpp" #include "fvm_lowered_cell_impl.hpp" -#include "util/meta.hpp" -#include "util/maputil.hpp" #include "util/rangeutil.hpp" #include "util/span.hpp" -#include "util/transform.hpp" #include "common.hpp" #include "mech_private_field_access.hpp" @@ -233,7 +229,7 @@ TEST(fvm_lowered, target_handles) { descriptions[1].decorations.place(mlocation{2, 0.2}, synapse("exp2syn"), "syn2"); descriptions[1].decorations.place(mlocation{2, 0.8}, synapse("expsyn"), "syn3"); - descriptions[1].decorations.place(mlocation{0, 0}, threshold_detector{3.3}, "detector"); + descriptions[1].decorations.place(mlocation{0, 0}, threshold_detector{3.3*arb::units::mV}, "detector"); cable_cell cells[] = {descriptions[0], descriptions[1]}; @@ -289,9 +285,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::box(5., 80., 0.3), "clamp0"); + desc.decorations.place(mlocation{0,1}, i_clamp::box(5.*arb::units::ms, 80.*arb::units::ms, 0.3*arb::units::nA), "clamp0"); // On the soma CV, which is over the approximate interval: (cable 0 0 0.1) - desc.decorations.place(mlocation{0,0.05}, i_clamp::box(1., 2., 0.1), "clamp1"); + desc.decorations.place(mlocation{0,0.05}, i_clamp::box(1.*arb::units::ms, 2.*arb::units::ms, 0.1*arb::units::nA), "clamp1"); std::vector cells{desc}; @@ -360,7 +356,8 @@ TEST(fvm_lowered, ac_stimulus) { 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, phase), "clamp"); + dec.place(mlocation{0, 0}, + i_clamp({{0*arb::units::ms, 0*arb::units::nA}, {max_time*arb::units::ms, max_amplitude*arb::units::nA}, {max_time*arb::units::ms, 0*arb::units::nA}}, freq*arb::units::kHz, phase*arb::units::rad), "clamp"); std::vector cells = {cable_cell(tree, dec)}; cable_cell_global_properties gprop; @@ -474,12 +471,12 @@ TEST(fvm_lowered, derived_mechs) { } }; - float times[] = {10.f, 20.f}; + std::vector times{10.f, 20.f}; auto decomp = partition_load_balance(rec, context); simulation sim(rec, context, decomp); - sim.add_sampler(all_probes, explicit_schedule(times), sampler); - sim.run(30.0, 1.f/1024); + sim.add_sampler(all_probes, explicit_schedule_from_milliseconds(times), sampler); + sim.run(30.0*arb::units::ms, 1.f/1024*arb::units::ms); ASSERT_EQ(2u, samples[0].size()); ASSERT_EQ(2u, samples[1].size()); @@ -509,7 +506,7 @@ TEST(fvm_lowered, null_region) { auto decomp = partition_load_balance(rec, context); simulation sim(rec, context, decomp); - EXPECT_NO_THROW(sim.run(30.0, 1.f/1024)); + EXPECT_NO_THROW(sim.run(30.0*arb::units::ms, 1.f/1024*arb::units::ms)); } @@ -830,7 +827,7 @@ TEST(fvm_lowered, post_events_shared_state) { auto ndetectors = detectors_per_cell_[gid]; auto offset = 1.0 / ndetectors; for (unsigned i = 0; i < ndetectors; ++i) { - decor.place(arb::mlocation{0, offset * i}, arb::threshold_detector{10}, "detector"+std::to_string(i)); + decor.place(arb::mlocation{0, offset * i}, arb::threshold_detector{10*arb::units::mV}, "detector"+std::to_string(i)); } decor.place(arb::mlocation{0, 0.5}, synapse_, "syanpse"); @@ -924,15 +921,15 @@ TEST(fvm_lowered, label_data) { decor.set_default(arb::cv_policy_fixed_per_branch(10)); decor.place(uniform(all(), 0, 3, 42), arb::synapse("expsyn"), "4_synapses"); decor.place(uniform(all(), 4, 4, 42), arb::synapse("expsyn"), "1_synapse"); - decor.place(uniform(all(), 5, 5, 42), arb::threshold_detector{10}, "1_detector"); + decor.place(uniform(all(), 5, 5, 42), arb::threshold_detector{10*arb::units::mV}, "1_detector"); cells_.push_back(arb::cable_cell(arb::morphology(tree), decor)); } { arb::decor decor; decor.set_default(arb::cv_policy_fixed_per_branch(10)); - decor.place(uniform(all(), 0, 2, 24), arb::threshold_detector{10}, "3_detectors"); - decor.place(uniform(all(), 3, 4, 24), arb::threshold_detector{10}, "2_detectors"); + decor.place(uniform(all(), 0, 2, 24), arb::threshold_detector{10*arb::units::mV}, "3_detectors"); + decor.place(uniform(all(), 3, 4, 24), arb::threshold_detector{10*arb::units::mV}, "2_detectors"); decor.place(uniform(all(), 5, 6, 24), arb::junction("gj"), "2_gap_junctions"); decor.place(uniform(all(), 7, 7, 24), arb::junction("gj"), "1_gap_junction"); diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp index 7a03cd8022..77d4f3c70c 100644 --- a/test/unit/test_lif_cell_group.cpp +++ b/test/unit/test_lif_cell_group.cpp @@ -12,8 +12,6 @@ #include #include -#include "lif_cell_group.hpp" - using namespace arb; // Simple ring network of LIF neurons. // with one regularly spiking cell (fake cell) connected to the first cell in the ring. @@ -62,7 +60,7 @@ class ring_recipe: public arb::recipe { // regularly spiking cell. if (gid == 0) { // Produces just a single spike at time 0ms. - return spike_source_cell("src", explicit_schedule({0.f})); + return spike_source_cell("src", explicit_schedule_from_milliseconds({0.})); } // LIF cell. auto cell = lif_cell("src", "tgt"); @@ -149,7 +147,7 @@ class probe_recipe: public arb::recipe { } } std::vector event_generators(cell_gid_type) const override { - return {regular_generator({"tgt"}, 200.0, 2.0, 1.0, 6.0)}; + return {regular_generator({"tgt"}, 200.0, 2.0*arb::units::ms, 1.0*arb::units::ms, 6.0*arb::units::ms)}; } size_t n_conn_ = 0; @@ -195,10 +193,7 @@ TEST(lif_cell_group, spikes) { events.push_back({0, {{0, 50, 1000}}}); sim.inject_events(events); - - time_type tfinal = 100; - time_type dt = 0.01; - sim.run(tfinal, dt); + sim.run(100*arb::units::ms, 0.01*arb::units::ms); // we expect 4 spikes: 2 by both neurons EXPECT_EQ(4u, sim.num_spikes()); @@ -211,9 +206,6 @@ TEST(lif_cell_group, ring) double weight = 1000; double delay = 1; - // Total simulation time. - time_type simulation_time = 100; - auto recipe = ring_recipe(num_lif_cells, weight, delay); // Creates a simulation with a ring recipe of lif neurons simulation sim(recipe); @@ -227,7 +219,7 @@ TEST(lif_cell_group, ring) ); // Runs the simulation for simulation_time with given timestep - sim.run(simulation_time, 0.01); + sim.run(100*arb::units::ms, 0.01*arb::units::ms); // The total number of cells in all the cell groups. // There is one additional fake cell (regularly spiking cell). EXPECT_EQ(num_lif_cells + 1u, recipe.num_cells()); @@ -274,7 +266,7 @@ TEST(lif_cell_group, probe) { auto rec = probe_recipe{}; auto sim = simulation(rec); - sim.add_sampler(all_probes, regular_schedule(0.025), fun); + sim.add_sampler(all_probes, regular_schedule(0.025*arb::units::ms), fun); std::vector spikes; @@ -282,7 +274,7 @@ TEST(lif_cell_group, probe) { [&spikes](const std::vector& spk) { for (const auto& s: spk) spikes.push_back(s.time); } ); - sim.run(10, 0.005); + sim.run(10*arb::units::ms, 0.005*arb::units::ms); std::vector exp = {{ 0, -18 }, { 0.025, -17.9750624 }, { 0.05, -17.9502492 }, @@ -710,7 +702,7 @@ TEST(lif_cell_group, probe_with_connections) { auto rec = probe_recipe{5}; auto sim = simulation(rec); - sim.add_sampler(all_probes, regular_schedule(0.025), fun); + sim.add_sampler(all_probes, regular_schedule(0.025*arb::units::ms), fun); std::vector spikes; @@ -718,7 +710,7 @@ TEST(lif_cell_group, probe_with_connections) { [&spikes](const std::vector& spk) { for (const auto& s: spk) spikes.push_back(s.time); } ); - sim.run(10, 0.005); + sim.run(10*arb::units::ms, 0.005*arb::units::ms); std::vector exp = {{ 0, -18 }, { 0.025, -17.9750624 }, { 0.05, -17.9502492 }, diff --git a/test/unit/test_merge_events.cpp b/test/unit/test_merge_events.cpp index 7b17c4094b..98060b610e 100644 --- a/test/unit/test_merge_events.cpp +++ b/test/unit/test_merge_events.cpp @@ -28,13 +28,15 @@ static void merge_events( const pse_vector& old_events, pse_vector& pending, std::vector& generators, - pse_vector& new_events) -{ + pse_vector& new_events) { util::sort(pending); - merge_cell_events(t_from, t_to, util::range_pointer_view(old_events), util::range_pointer_view(pending), generators, new_events); + merge_cell_events(t_from, t_to, + util::range_pointer_view(old_events), + util::range_pointer_view(pending), + generators, + new_events); } - std::vector empty_gens; // Test the trivial case of merging empty sets @@ -157,7 +159,7 @@ TEST(merge_events, X) {0, 26, 4}, }; - auto gen = regular_generator({"l0"}, 42.f, t0, 5); + auto gen = regular_generator({"l0"}, 42.f, t0*arb::units::ms, 5*arb::units::ms); gen.resolve_label([](const cell_local_label_type&) {return 2;}); std::vector generators = {gen}; @@ -203,8 +205,8 @@ TEST(merge_events, tourney_seq) util::sort(expected); auto - g1 = explicit_generator(l0, w1, times), - g2 = explicit_generator(l0, w2, times); + g1 = explicit_generator_from_milliseconds(l0, w1, times), + g2 = explicit_generator_from_milliseconds(l0, w2, times); g1.resolve_label([](const cell_local_label_type&) {return 0;}); g2.resolve_label([](const cell_local_label_type&) {return 0;}); @@ -223,9 +225,7 @@ TEST(merge_events, tourney_seq) } // Test the tournament tree on a large set of Poisson generators. -TEST(merge_events, tourney_poisson) -{ - using rndgen = std::mt19937_64; +TEST(merge_events, tourney_poisson) { // Number of poisson generators. // Not a power of 2, so that there will be "null" leaf nodes in the // tournament tree. @@ -241,8 +241,8 @@ TEST(merge_events, tourney_poisson) float weight = i; // the first and last generators have the same seed to test that sorting // of events with the same time but different weights works properly. - rndgen G(i%(ngen-1)); - auto gen = poisson_generator(label, weight, t0, lambda, G); + auto G = i%(ngen-1); + auto gen = poisson_generator(label, weight, t0*arb::units::ms, lambda*arb::units::kHz, G); gen.resolve_label([lid](const cell_local_label_type&) {return lid;}); generators.push_back(std::move(gen)); } diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp index 90771c9800..9e2aa6038c 100644 --- a/test/unit/test_probe.cpp +++ b/test/unit/test_probe.cpp @@ -1,7 +1,6 @@ #include #include -#include #include #include @@ -110,7 +109,7 @@ void run_v_i_probe_test(context ctx) { bs.decorations.set_default(cv_policy_fixed_per_branch(1)); - auto stim = i_clamp::box(0, 100, 0.3); + auto stim = i_clamp::box(0.*arb::units::ms, 100*arb::units::ms, 0.3*arb::units::nA); bs.decorations.place(mlocation{1, 1}, stim, "clamp"); cable1d_recipe rec((cable_cell(bs))); @@ -778,7 +777,7 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { cv_policy policy = cv_policy_fixed_per_branch(n_cv); d.set_default(policy); - d.place(mlocation{0, 0}, i_clamp(0.3), "clamp"); + d.place(mlocation{0, 0}, i_clamp(0.3*arb::units::nA), "clamp"); // The time constant will be membrane capacitance / membrane conductance. // For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance @@ -786,7 +785,7 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { d.paint(reg::all(), density("ca_linear", {{"g", 0.01}})); // [S/cm²] d.set_default(membrane_capacitance{0.01}); // [F/m²] - const double tau = 0.1; // [ms] + auto tau = 0.1*arb::units::ms; cable1d_recipe rec(cable_cell(m, d)); rec.catalogue() = cat; @@ -826,7 +825,7 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { std::vector i_axial(n_axial_probe); std::vector i_memb(n_cv), i_stim(n_cv); - sim.add_sampler(all_probes, explicit_schedule({20*tau}), + sim.add_sampler(all_probes, explicit_schedule(std::vector{20*tau}), [&](probe_metadata pm, std::size_t n_sample, const sample_record* samples) { @@ -869,8 +868,8 @@ void run_axial_and_ion_current_sampled_probe_test(context ctx) { } }); - const double dt = 0.025; // [ms] - sim.run(20*tau+dt, dt); + auto dt = 0.025*arb::units::ms; // [ms] + sim.run(20*tau + dt, dt); ASSERT_EQ(n_cv, i_memb.size()); ASSERT_EQ(n_cv, i_stim.size()); @@ -913,11 +912,11 @@ auto run_simple_samplers(const arb::context& ctx, std::vector> traces(n_probe); for (unsigned i = 0; i cells = {{bs.morph, d0, bs.labels}, {bs.morph, d1, bs.labels}}; @@ -1038,8 +1037,8 @@ void run_total_current_probe_test(context ctx) { // For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance // to 0.01 F/m². - const double tau = 0.1; // [ms] - d0.place(mlocation{0, 0}, i_clamp(0.3), "clamp0"); + auto tau = 0.1*arb::units::ms; // [ms] + d0.place(mlocation{0, 0}, i_clamp(0.3*arb::units::nA), "clamp0"); d0.paint(reg::all(), density("ca_linear", {{"g", 0.01}})); // [S/cm²] d0.set_default(membrane_capacitance{0.01}); // [F/m²] @@ -1069,19 +1068,21 @@ void run_total_current_probe_test(context ctx) { for (unsigned i = 0; i<2; ++i) { SCOPED_TRACE(i); - const double t_end = 21*tau; // [ms] + auto t_end = 21*tau; // [ms] - traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, + traces[i] = run_simple_sampler, mcable_list>(ctx, t_end.value(), cells, {i, "Itotal"}, - cable_probe_total_current_cell{}, {tau, 20*tau}).at(0); + cable_probe_total_current_cell{}, + {tau.value(), 20*tau.value()}).at(0); - ion_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, + ion_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end.value(), cells, {i, "Iion"}, - cable_probe_total_ion_current_cell{}, {tau, 20*tau}).at(0); + cable_probe_total_ion_current_cell{}, + {tau.value(), 20*tau.value()}).at(0); - stim_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end, cells, + stim_traces[i] = run_simple_sampler, mcable_list>(ctx, t_end.value(), cells, {i, "Istim"}, - cable_probe_stimulus_current_cell{}, {tau, 20*tau}).at(0); + cable_probe_stimulus_current_cell{}, {tau.value(), 20*tau.value()}).at(0); ASSERT_EQ(2u, traces[i].size()); ASSERT_EQ(2u, ion_traces[i].size()); @@ -1157,19 +1158,21 @@ void run_stimulus_probe_test(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 stim_from = 0.*arb::units::ms; + auto stim_until = 1.*arb::units::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::box(0., stim_until, 10.), "clamp0"); - d0.place(mlocation{0, 0.5}, i_clamp::box(0., stim_until, 20.), "clamp1"); + d0.place(mlocation{0, 0.5}, i_clamp::box(stim_from, stim_until, 10.*arb::units::nA), "clamp0"); + d0.place(mlocation{0, 0.5}, i_clamp::box(stim_from, stim_until, 20.*arb::units::nA), "clamp1"); double expected_stim0 = 30; d1.set_default(policy); - d1.place(mlocation{0, 1}, i_clamp::box(0., stim_until, 30.), "clamp0"); - d1.place(mlocation{0, 1}, i_clamp::box(0., stim_until, -10.), "clamp1"); + d1.place(mlocation{0, 1}, i_clamp::box(stim_from, stim_until, 30.*arb::units::nA), "clamp0"); + d1.place(mlocation{0, 1}, i_clamp::box(stim_from, stim_until, -10.*arb::units::nA), "clamp1"); double expected_stim1 = 20; std::vector cells = {{m, d0}, {m, d1}}; @@ -1179,9 +1182,9 @@ void run_stimulus_probe_test(context ctx) { trace_data, mcable_list> traces[2]; for (unsigned i: {0u, 1u}) { - traces[i] = run_simple_sampler, mcable_list>(ctx, 2.5*stim_until, cells, {i, "Istim"}, + traces[i] = run_simple_sampler, mcable_list>(ctx, 2.5*stim_until.value(), cells, {i, "Istim"}, cable_probe_stimulus_current_cell{}, - {stim_until/2, 2*stim_until}).at(0); + {stim_until.value()/2, 2*stim_until.value()}).at(0); ASSERT_EQ(3u, traces[i].meta.size()); for ([[maybe_unused]] unsigned cv: {0u, 1u, 2u}) { diff --git a/test/unit/test_recipe.cpp b/test/unit/test_recipe.cpp index 27b04cb749..7c67d3c760 100644 --- a/test/unit/test_recipe.cpp +++ b/test/unit/test_recipe.cpp @@ -17,83 +17,72 @@ #include "../common_cells.hpp" using namespace arb; -using arb::util::make_span; namespace { - class custom_recipe: public recipe { - public: - custom_recipe(std::vector cells, - std::vector> conns, - std::vector> gjs, - std::vector> gens): - num_cells_(cells.size()), - connections_(conns), - gap_junctions_(gjs), - event_generators_(gens), - cells_(cells) {} - - cell_size_type num_cells() const override { - return num_cells_; - } - arb::util::unique_any get_cell_description(cell_gid_type gid) const override { - return cells_.at(gid); - } - cell_kind get_cell_kind(cell_gid_type gid) const override { - return cell_kind::cable; - } - std::vector gap_junctions_on(cell_gid_type gid) const override { - return gap_junctions_.at(gid); - } - std::vector connections_on(cell_gid_type gid) const override { - return connections_.at(gid); - } - std::vector event_generators(cell_gid_type gid) const override { - return event_generators_.at(gid); - } - std::any get_global_properties(cell_kind) const override { - arb::cable_cell_global_properties a; - a.default_parameters = arb::neuron_parameter_defaults; - return a; - } - - private: - cell_size_type num_cells_; - std::vector> connections_; - std::vector> gap_junctions_; - std::vector> event_generators_; - std::vector cells_; - }; - - cable_cell custom_cell(cell_size_type num_detectors, cell_size_type num_synapses, cell_size_type num_gj) { - arb::segment_tree tree; - tree.append(arb::mnpos, {0,0,0,10}, {0,0,20,10}, 1); // soma - tree.append(0, {0,0, 20, 2}, {0,0, 320, 2}, 3); // dendrite - - arb::cable_cell cell(tree, {}); - - arb::decor decorations; - - // Add a num_detectors detectors to the cell. - for (auto i: util::make_span(num_detectors)) { - decorations.place(arb::mlocation{0,(double)i/num_detectors}, arb::threshold_detector{10}, "detector"+std::to_string(i)); - } - - // Add a num_synapses synapses to the cell. - for (auto i: util::make_span(num_synapses)) { - decorations.place(arb::mlocation{0,(double)i/num_synapses}, arb::synapse("expsyn"), "synapse"+std::to_string(i)); - } - - // Add a num_gj gap_junctions to the cell. - for (auto i: util::make_span(num_gj)) { - decorations.place(arb::mlocation{0,(double)i/num_gj}, arb::junction("gj"), "gapjunction"+std::to_string(i)); - } - - return arb::cable_cell(tree, decorations); +struct custom_recipe: public recipe { + custom_recipe(std::vector cells, + std::vector> conns, + std::vector> gjs, + std::vector> gens): + num_cells_(cells.size()), + connections_(std::move(conns)), + gap_junctions_(std::move(gjs)), + event_generators_(std::move(gens)), + cells_(std::move(cells)) { + gprop.default_parameters = arb::neuron_parameter_defaults; } + + cell_size_type num_cells() const override { return num_cells_; } + arb::util::unique_any get_cell_description(cell_gid_type gid) const override { return cells_.at(gid); } + cell_kind get_cell_kind(cell_gid_type gid) const override { return cell_kind::cable; } + std::vector gap_junctions_on(cell_gid_type gid) const override { return gap_junctions_.at(gid); } + std::vector connections_on(cell_gid_type gid) const override { return connections_.at(gid); } + std::vector event_generators(cell_gid_type gid) const override { return event_generators_.at(gid); } + std::any get_global_properties(cell_kind) const override { return gprop; } + +private: + cell_size_type num_cells_; + std::vector> connections_; + std::vector> gap_junctions_; + std::vector> event_generators_; + std::vector cells_; + arb::cable_cell_global_properties gprop; +}; + +cable_cell custom_cell(cell_size_type num_detectors, cell_size_type num_synapses, cell_size_type num_gj) { + arb::segment_tree tree; + tree.append(arb::mnpos, {0,0,0,10}, {0,0,20,10}, 1); // soma + tree.append(0, {0,0, 20, 2}, {0,0, 320, 2}, 3); // dendrite + + arb::cable_cell cell(tree, {}); + + arb::decor decorations; + + // Add a num_detectors detectors to the cell. + for (auto i: util::make_span(num_detectors)) { + decorations.place(arb::mlocation{0,(double)i/num_detectors}, + arb::threshold_detector{10*arb::units::ms}, "detector"+std::to_string(i)); + } + + // Add a num_synapses synapses to the cell. + for (auto i: util::make_span(num_synapses)) { + decorations.place(arb::mlocation{0,(double)i/num_synapses}, + arb::synapse("expsyn"), + "synapse"+std::to_string(i)); + } + + // Add a num_gj gap_junctions to the cell. + for (auto i: util::make_span(num_gj)) { + decorations.place(arb::mlocation{0,(double)i/num_gj}, + arb::junction("gj"), + "gapjunction"+std::to_string(i)); + } + + return arb::cable_cell(tree, decorations); } +} // namespace -TEST(recipe, gap_junctions) -{ +TEST(recipe, gap_junctions) { auto context = make_context({arbenv::default_concurrency(), -1}); auto cell_0 = custom_cell(0, 0, 3); @@ -203,18 +192,18 @@ TEST(recipe, event_generators) { auto cell_1 = custom_cell(2, 1, 0); { std::vector - gens_0 = {arb::explicit_generator({"synapse0"}, 0.1, std::vector{1.0}), - arb::explicit_generator({"synapse1"}, 0.1, std::vector{2.0})}, - gens_1 = {arb::explicit_generator({"synapse0"}, 0.1, std::vector{1.0})}; + gens_0 = {arb::explicit_generator_from_milliseconds({"synapse0"}, 0.1, std::vector{1.0}), + arb::explicit_generator_from_milliseconds({"synapse1"}, 0.1, std::vector{2.0})}, + gens_1 = {arb::explicit_generator_from_milliseconds({"synapse0"}, 0.1, std::vector{1.0})}; auto recipe_0 = custom_recipe({cell_0, cell_1}, {{}, {}}, {{}, {}}, {gens_0, gens_1}); auto decomp_0 = partition_load_balance(recipe_0, context); - EXPECT_NO_THROW(simulation(recipe_0, context, decomp_0).run(1, 0.1)); + EXPECT_NO_THROW(simulation(recipe_0, context, decomp_0).run(1*arb::units::ms, 0.1*arb::units::ms)); } { std::vector - gens_0 = {arb::regular_generator({"totally-not-a-synapse-42"}, 0.1, 0, 0.001)}, + gens_0 = {arb::regular_generator({"totally-not-a-synapse-42"}, 0.1, 0*arb::units::ms, 0.001*arb::units::ms)}, gens_1 = {}; auto recipe_0 = custom_recipe({cell_0, cell_1}, {{}, {}}, {{}, {}}, {gens_0, gens_1}); diff --git a/test/unit/test_schedule.cpp b/test/unit/test_schedule.cpp index 8452385160..59fd06220c 100644 --- a/test/unit/test_schedule.cpp +++ b/test/unit/test_schedule.cpp @@ -87,7 +87,7 @@ TEST(schedule, regular) { // Use exact fp representations for strict equality testing. std::vector expected = {0, 0.25, 0.5, 0.75, 1.0}; - schedule S = regular_schedule(0.25); + schedule S = regular_schedule(0.25*arb::units::ms); EXPECT_EQ(expected, as_vector(S.events(0, 1.25))); S.reset(); @@ -100,12 +100,12 @@ TEST(schedule, regular) { TEST(schedule, regular_invariants) { SCOPED_TRACE("regular_invariants"); - run_invariant_checks(regular_schedule(0.3), 3, 12, 7); + run_invariant_checks(regular_schedule(0.3*arb::units::ms), 3, 12, 7); } TEST(schedule, regular_reset) { SCOPED_TRACE("regular_reset"); - run_reset_check(regular_schedule(0.3), 3, 12, 7); + run_reset_check(regular_schedule(0.3*arb::units::ms), 3, 12, 7); } TEST(schedule, regular_rounding) { @@ -120,7 +120,7 @@ TEST(schedule, regular_rounding) { time_type t0 = t1-10*dt; time_type t2 = t1+10*dt; - schedule S = regular_schedule(t0, dt); + schedule S = regular_schedule(t0*arb::units::ms, dt*arb::units::ms); auto int_l = as_vector(S.events(t0, t1)); auto int_r = as_vector(S.events(t1, t2)); @@ -144,10 +144,10 @@ TEST(schedule, regular_rounding) { } TEST(schedule, explicit_schedule) { - time_type times[] = {0.1, 0.3, 1.0, 1.25, 1.7, 2.2}; - std::vector expected = {0.1, 0.3, 1.0}; + std::vector times{0.1, 0.3, 1.0, 1.25, 1.7, 2.2}; + std::vector expected{0.1, 0.3, 1.0}; - schedule S = explicit_schedule(times); + schedule S = explicit_schedule_from_milliseconds(times); EXPECT_EQ(expected, as_vector(S.events(0, 1.25))); S.reset(); @@ -161,47 +161,19 @@ TEST(schedule, explicit_schedule) { TEST(schedule, explicit_invariants) { SCOPED_TRACE("explicit_invariants"); - time_type times[] = {0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13}; - run_invariant_checks(explicit_schedule(times), 0.4, 10.2, 5); + std::vector times{0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13}; + run_invariant_checks(explicit_schedule_from_milliseconds(times), 0.4, 10.2, 5); } TEST(schedule, explicit_reset) { SCOPED_TRACE("explicit_reset"); - time_type times[] = {0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13}; - run_reset_check(explicit_schedule(times), 0.4, 10.2, 5); + std::vector times{0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13}; + run_reset_check(explicit_schedule_from_milliseconds(times), 0.4, 10.2, 5); } -// A Uniform Random Bit Generator[*] adaptor that deliberately -// skews the generated numbers by raising their quantile to -// the given power. -// -// [*] Not actually uniform. - -template -struct skew_adaptor { - using result_type = typename RNG::result_type; - static constexpr result_type min() { return RNG::min(); } - static constexpr result_type max() { return RNG::max(); } - - explicit skew_adaptor(double power): power_(power) {} - result_type operator()() { - constexpr double scale = (double)(max()-min()); - constexpr double ooscale = 1./scale; - - double x = ooscale*(G_()-min()); - x = std::pow(x, power_); - return min()+scale*x; - } - -private: - RNG G_; - double power_; -}; - -template -double poisson_schedule_dispersion(int nbin, double rate_kHz, RNG& G) { - schedule S = poisson_schedule(rate_kHz, G); +double poisson_schedule_dispersion(int nbin, double rate_kHz) { + schedule S = poisson_schedule(rate_kHz*arb::units::kHz); std::vector bin(nbin); for (auto t: time_range(S.events(0, nbin))) { @@ -238,8 +210,7 @@ TEST(schedule, poisson_uniformity) { constexpr double chi2_lb = 888.56352318146696; constexpr double chi2_ub = 1118.9480663231843; - std::mt19937_64 G; - double dispersion = poisson_schedule_dispersion(N, .813, G); + double dispersion = poisson_schedule_dispersion(N, .813); double test_value = N*dispersion; EXPECT_GT(test_value, chi2_lb); EXPECT_LT(test_value, chi2_ub); @@ -247,34 +218,12 @@ TEST(schedule, poisson_uniformity) { // Run one sample K-S test for uniformity, with critical // value for the finite K-S statistic Dn of α=0.01. - schedule S = poisson_schedule(100., G); + schedule S = poisson_schedule(100.*arb::units::kHz); auto events = as_vector(S.events(0,1)); int n = (int)events.size(); double dn = ks::dn_statistic(events); EXPECT_LT(ks::dn_cdf(dn, n), 0.99); - - // Check that these tests fail for a non-Poisson - // source. - - skew_adaptor W(1.5); - dispersion = poisson_schedule_dispersion(N, .813, W); - test_value = N*dispersion; - - EXPECT_FALSE(test_value>=chi2_lb && test_value<=chi2_ub); - - S = poisson_schedule(100., W); - events = as_vector(S.events(0,1)); - n = (int)events.size(); - dn = ks::dn_statistic(events); - - // This test is currently failing, because we can't - // use a sufficiently high `n` in the `dn_cdf` function - // to get enough discrimination from the K-S test at - // 1%. TODO: Fix this by implementing n>140 case in - // `dn_cdf`. - - // EXPECT_GT(ks::dn_cdf(dn, n), 0.99); } TEST(schedule, poisson_rate) { @@ -284,37 +233,26 @@ TEST(schedule, poisson_rate) { constexpr double alpha = 0.01; constexpr double lambda = 123.4; - std::mt19937_64 G; - schedule S = poisson_schedule(lambda, G); + schedule S = poisson_schedule(lambda*arb::units::kHz); int n = (int)time_range(S.events(0, 1)).size(); double cdf = poisson::poisson_cdf_approx(n, lambda); EXPECT_GT(cdf, alpha/2); EXPECT_LT(cdf, 1-alpha/2); - - // Check that the test fails for a non-Poisson - // source. - - skew_adaptor W(1.5); - S = poisson_schedule(lambda, W); - n = (int)time_range(S.events(0, 1)).size(); - cdf = poisson::poisson_cdf_approx(n, lambda); - - EXPECT_FALSE(cdf>=alpha/2 && cdf<=1-alpha/2); } TEST(schedule, poisson_invariants) { SCOPED_TRACE("poisson_invariants"); - std::mt19937_64 G; - G.discard(100); - run_invariant_checks(poisson_schedule(0.81, G), 5.1, 15.3, 7); + auto sched = poisson_schedule(0.81*arb::units::kHz); + sched.discard(100); + run_invariant_checks(sched, 5.1, 15.3, 7); } TEST(schedule, poisson_reset) { SCOPED_TRACE("poisson_reset"); - std::mt19937_64 G; - G.discard(200); - run_reset_check(poisson_schedule(.11, G), 1, 10, 7); + auto sched = poisson_schedule(0.11*arb::units::kHz); + sched.discard(200); + run_reset_check(sched, 1, 10, 7); } TEST(schedule, poisson_offset) { @@ -323,42 +261,38 @@ TEST(schedule, poisson_offset) { const double offset = 3.3; - std::mt19937_64 G1; - G1.discard(300); - + auto sched1 = poisson_schedule(.234*arb::units::kHz); + sched1.discard(300); std::vector expected; - for (auto t: as_vector(poisson_schedule(.234, G1).events(0., 100.))) { + for (auto t: as_vector(sched1.events(0., 100.))) { t += offset; if (t<100.) { expected.push_back(t); } } - std::mt19937_64 G2; - G2.discard(300); - + auto sched2 = poisson_schedule(offset*arb::units::ms, .234*arb::units::kHz); + sched2.discard(300); EXPECT_TRUE(seq_almost_eq(expected, - as_vector(poisson_schedule(offset, .234, G2).events(0., 100.)))); + as_vector(sched2.events(0., 100.)))); } TEST(schedule, poisson_offset_reset) { SCOPED_TRACE("poisson_reset"); - std::mt19937_64 G; - G.discard(400); - run_reset_check(poisson_schedule(3.3, 9.1, G), 1, 10, 7); + auto sched = poisson_schedule(3.3*arb::units::ms, 0.81*arb::units::kHz); + sched.discard(400); + run_reset_check(sched, 1, 10, 7); } TEST(schedule, poisson_tstop) { SCOPED_TRACE("poisson_tstop"); - std::mt19937_64 G; - G.discard(500); - - const double tstop = 50; - auto const times = as_vector(poisson_schedule(0, .234, G, tstop).events(0., 100.)); + auto sched = poisson_schedule(0*arb::units::ms, 0.234*arb::units::kHz, default_seed, 50*arb::units::ms); + sched.discard(500); + auto const times = as_vector(sched.events(0., 100.)); auto const max = std::max_element(begin(times), end(times)); EXPECT_TRUE(max != end(times)); - EXPECT_TRUE(*max <= tstop); + EXPECT_TRUE(*max <= 50); } diff --git a/test/unit/test_sde.cpp b/test/unit/test_sde.cpp index 7383b7ef02..e877a33713 100644 --- a/test/unit/test_sde.cpp +++ b/test/unit/test_sde.cpp @@ -409,7 +409,7 @@ TEST(sde, reproducibility) { // simulation parameters unsigned ncells = 4; unsigned ncvs = 2; - double const dt = 0.5; + auto const dt = 0.5*arb::units::ms; unsigned nsteps = 6; // Decorations with a bunch of stochastic processes @@ -490,7 +490,7 @@ TEST(sde, normality) { unsigned ncells = 4; unsigned nsynapses = 100; unsigned ncvs = 100; - double const dt = 0.5; + auto dt = 0.5*arb::units::ms;; unsigned nsteps = 50; // make labels (and locations for synapses) @@ -647,7 +647,7 @@ TEST(sde, solver) { unsigned ncells = 4; unsigned nsynapses = 2000; unsigned ncvs = 1; - double const dt = 1.0/512; // need relatively small time steps due to low accuracy + auto dt = 1.0/512*arb::units::ms; // need relatively small time steps due to low accuracy unsigned nsteps = 100; unsigned nsims = 4; @@ -773,7 +773,7 @@ TEST(sde, solver) { auto test = [&] (auto func, const auto& stats) { for (unsigned int i=1; i event_generators(arb::cell_gid_type gid) const override { std::vector res; - if (!gid) res.push_back(arb::regular_generator({"synapse"}, 1, 0.5, 0.73)); + if (!gid) res.push_back(arb::regular_generator({"synapse"}, 1, 0.5*arb::units::ms, 0.73*arb::units::ms)); return {}; } @@ -190,8 +190,8 @@ void sampler(arb::probe_metadata pm, } TEST(serdes, single_cell) { - double dt = 0.5; - double T = 5; + auto dt = 0.5*arb::units::ms; + auto T = 5*arb::units::ms; // Result std::vector result_pre; @@ -231,8 +231,8 @@ TEST(serdes, single_cell) { } TEST(serdes, network) { - double dt = 0.5; - double T = 5; + auto dt = 0.5*arb::units::ms; + auto T = 5*arb::units::ms; // Result std::vector result_pre; @@ -248,7 +248,7 @@ TEST(serdes, network) { model.num = 10; auto simulation = arb::simulation{model}; simulation.add_sampler(arb::all_probes, - arb::regular_schedule(dt), + arb::regular_schedule(dt*arb::units::ms), sampler); // Run simulation forward && snapshot diff --git a/test/unit/test_simulation.cpp b/test/unit/test_simulation.cpp index 88bdec4072..b11ac2d804 100644 --- a/test/unit/test_simulation.cpp +++ b/test/unit/test_simulation.cpp @@ -1,6 +1,5 @@ #include -#include #include #include @@ -17,7 +16,6 @@ #include "util/rangeutil.hpp" #include "util/transform.hpp" -#include "common.hpp" using namespace arb; struct play_spikes: public recipe { @@ -25,10 +23,7 @@ struct play_spikes: public recipe { cell_size_type num_cells() const override { return spike_times_.size(); } cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::spike_source; } - util::unique_any get_cell_description(cell_gid_type gid) const override { - return spike_source_cell("src", spike_times_.at(gid)); - } - + util::unique_any get_cell_description(cell_gid_type gid) const override { return spike_source_cell("src", spike_times_.at(gid)); } std::vector spike_times_; }; @@ -52,7 +47,7 @@ TEST(simulation, null) { auto c = arb::make_context(); auto d = arb::partition_load_balance(r, c); auto s = arb::simulation(r, c, d); - s.run(0.05, 0.01); + s.run(0.05*arb::units::ms, 0.01*arb::units::ms); } // Test with simulation builder @@ -60,22 +55,23 @@ TEST(simulation, null_builder) { auto r = null_recipe{}; { arb::simulation s = arb::simulation::create(r); - s.run(0.05, 0.01); + s.run(0.05*arb::units::ms, + 0.01*arb::units::ms); } { arb::simulation s = arb::simulation::create(r).set_seed(42); - s.run(0.05, 0.01); + s.run(0.05*arb::units::ms, 0.01*arb::units::ms); } { auto c = arb::make_context(); arb::simulation s = arb::simulation::create(r).set_context(c); - s.run(0.05, 0.01); + s.run(0.05*arb::units::ms, 0.01*arb::units::ms); } { auto c = arb::make_context(); auto d = arb::partition_load_balance(r, c); arb::simulation s = arb::simulation::create(r).set_context(c).set_decomposition(d); - s.run(0.05, 0.01); + s.run(0.05*arb::units::ms, 0.01*arb::units::ms); } } @@ -85,7 +81,7 @@ TEST(simulation, spike_global_callback) { std::vector spike_times; for (unsigned i = 0; i expected_spikes; @@ -107,7 +103,7 @@ TEST(simulation, spike_global_callback) { double tfinal = 0.7*t_max; constexpr double dt = 0.01; - sim.run(tfinal, dt); + sim.run(tfinal*arb::units::ms, dt*arb::units::ms); auto spike_lt = [](spike a, spike b) { return a.time trigger_times = {1., 2., 3.}; double delay = 10; unsigned n = 5; - lif_chain rec(n, delay, explicit_schedule(trigger_times)); + lif_chain rec(n, delay, explicit_schedule_from_milliseconds(trigger_times)); // Expect spike times to be almost exactly according to trigger times, // plus delays along the chain of cells. @@ -203,7 +199,7 @@ TEST(simulation, restart) { double t = 0; do { double run_to = std::min(tfinal, t + run_time); - t = sim.run(run_to, dt); + t = sim.run(run_to*arb::units::ms, dt*arb::units::ms); ASSERT_EQ(t, run_to); } while (t // Test that a spike_source_cell_group identifies itself with the correct // cell_kind enum value. -TEST(spike_source, cell_kind) -{ +TEST(spike_source, cell_kind) { ss_recipe rec(1u, spike_source_cell("src", explicit_schedule({}))); cell_label_range srcs, tgts; spike_source_cell_group group({0}, rec, srcs, tgts); @@ -37,8 +36,7 @@ static std::vector spike_times(const std::vector& evs) { // Test that a spike_source_cell_group produces a sequence of spikes with spike // times corresponding to the underlying time_seq. -TEST(spike_source, matches_time_seq) -{ +TEST(spike_source, matches_time_seq) { auto test_seq = [](schedule seq) { ss_recipe rec(1u, spike_source_cell("src", seq)); cell_label_range srcs, tgts; @@ -57,16 +55,14 @@ TEST(spike_source, matches_time_seq) EXPECT_EQ(spike_times(group.spikes()), as_vector(seq.events(10, 20))); }; - std::mt19937_64 G; - test_seq(regular_schedule(0, 1)); - test_seq(poisson_schedule(10, G)); // produce many spikes in each interval - test_seq(poisson_schedule(1e-6, G)); // very unlikely to produce any spikes in either interval + test_seq(regular_schedule(1*arb::units::ms)); + test_seq(poisson_schedule(10*arb::units::kHz)); // produce many spikes in each interval + test_seq(poisson_schedule(1e-6*arb::units::kHz)); // very unlikely to produce any spikes in either interval } // Test that a spike_source_cell_group will produce the same sequence of spikes // after being reset. -TEST(spike_source, reset) -{ +TEST(spike_source, reset) { auto test_seq = [](schedule seq) { ss_recipe rec(1u, spike_source_cell("src", seq)); cell_label_range srcs, tgts; @@ -87,16 +83,14 @@ TEST(spike_source, reset) EXPECT_EQ(spikes1, spikes2); }; - std::mt19937_64 G; - test_seq(regular_schedule(0, 1)); - test_seq(poisson_schedule(10, G)); // produce many spikes in each interval - test_seq(poisson_schedule(1e-6, G)); // very unlikely to produce any spikes in either interval + test_seq(regular_schedule(10*arb::units::ms)); + test_seq(poisson_schedule(100*arb::units::kHz)); // produce many spikes in each interval + test_seq(poisson_schedule(1e-6*arb::units::kHz)); // very unlikely to produce any spikes in either interval } // Test that a spike_source_cell_group will produce the expected // output when the underlying time_seq is finite. -TEST(spike_source, exhaust) -{ +TEST(spike_source, exhaust) { // This test assumes that seq will exhaust itself before t=10 ms. auto test_seq = [](schedule seq) { ss_recipe rec(1u, spike_source_cell("src", seq)); @@ -112,12 +106,11 @@ TEST(spike_source, exhaust) EXPECT_LT(group.spikes().back().time, time_type(10)); }; - test_seq(regular_schedule(0, 1, 5)); - test_seq(explicit_schedule({0.3, 2.3, 4.7})); + test_seq(regular_schedule(0*arb::units::ms, 1*arb::units::ms, 5*arb::units::ms)); + test_seq(explicit_schedule_from_milliseconds(std::vector{0.3, 2.3, 4.7})); } -TEST(spike_source, multiple) -{ +TEST(spike_source, multiple) { // This test assumes that seq will exhaust itself before t=10 ms. auto test_seq = [](auto&&... seqs) { std::vector schedules{seqs...}; @@ -145,13 +138,7 @@ TEST(spike_source, multiple) EXPECT_LT(group.spikes().back().time, time_type(10)); }; - auto seqs = std::vector{regular_schedule(0, 1, 5), - explicit_schedule({0.3, 2.3, 4.7})}; + std::vector seqs{regular_schedule(0*arb::units::ms, 1*arb::units::ms, 5*arb::units::ms), + explicit_schedule_from_milliseconds(std::vector{0.3, 2.3, 4.7})}; test_seq(seqs); - test_seq(std::vector{regular_schedule(0, 1, 5), - explicit_schedule({0.3, 2.3, 4.7})}); - test_seq(regular_schedule(0, 1, 5), - explicit_schedule({0.3, 2.3, 4.7})); - auto reg_sched = regular_schedule(0, 1, 5); - test_seq(reg_sched, explicit_schedule({0.3, 2.3, 4.7})); } diff --git a/test/unit/test_spikes.cpp b/test/unit/test_spikes.cpp index c27b6d1990..031b9bce89 100644 --- a/test/unit/test_spikes.cpp +++ b/test/unit/test_spikes.cpp @@ -203,8 +203,8 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher) { } TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) { - double dt = 0.025; - double duration = 1; + auto dt = 0.025*arb::units::ms; + auto duration = 1*arb::units::ms; arb::segment_tree tree; tree.append(arb::mnpos, { -6.3, 0.0, 0.0, 6.3}, { 6.3, 0.0, 0.0, 6.3}, 1); @@ -223,8 +223,8 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) { for (unsigned i = 0; i < 8; i++) { arb::decor decor; decor.set_default(arb::cv_policy_every_segment()); - decor.place("mid"_lab, arb::threshold_detector{10}, "detector"); - decor.place("mid"_lab, arb::i_clamp::box(0.01+i*dt, duration, 0.5), "clamp"); + decor.place("mid"_lab, arb::threshold_detector{10*arb::units::mV}, "detector"); + decor.place("mid"_lab, arb::i_clamp::box(0.01*arb::units::ms + i*dt, duration, 0.5*arb::units::nA), "clamp"); decor.place("mid"_lab, arb::synapse("expsyn"), "synapse"); arb::cable_cell cell(morpho, decor, dict); @@ -243,7 +243,7 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) { } for (unsigned i = 1; i < spikes.size(); ++i) { - EXPECT_NEAR(dt, spikes[i].time - spikes[i-1].time, 1e-4); + EXPECT_NEAR(dt.value(), spikes[i].time - spikes[i-1].time, 1e-4); } } diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp index 2c00ce0382..8e56151faa 100644 --- a/test/unit/test_synapses.cpp +++ b/test/unit/test_synapses.cpp @@ -1,7 +1,6 @@ #include #include -#include #include #include @@ -10,7 +9,6 @@ #include #include "backends/multicore/fvm.hpp" -#include "util/maputil.hpp" #include "util/range.hpp" #include "../common_cells.hpp" diff --git a/test/unit/test_v_clamp.cpp b/test/unit/test_v_clamp.cpp index 2142d45ce6..cb5dd6d460 100644 --- a/test/unit/test_v_clamp.cpp +++ b/test/unit/test_v_clamp.cpp @@ -41,7 +41,7 @@ struct v_proc_recipe: public arb::recipe { {arb::cable_probe_membrane_voltage{"(location 0 0.625)"_ls}, "Um-(0, 0.625)"}}; // dend center: 0.75/2 + 0.25 } std::vector event_generators(arb::cell_gid_type) const override { - return {arb::regular_generator({"tgt"}, 5.0, 0.2, 0.05)}; + return {arb::regular_generator({"tgt"}, 5.0, 0.2*arb::units::ms, 0.05*arb::units::ms)}; } std::any get_global_properties(arb::cell_kind) const override { return gprop; } @@ -90,8 +90,8 @@ TEST(v_process, clamp) { } }; auto sim = arb::simulation(v_proc_recipe{true, false}); - sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun); - sim.run(1.0, 0.005); + sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05*arb::units::ms), fun); + sim.run(1.0*arb::units::ms, 0.005*arb::units::ms); um_s_type exp_soma{{ 0, -65 }, { 0.05, -42 }, @@ -158,8 +158,8 @@ TEST(v_process, limit) { } }; auto sim = arb::simulation(v_proc_recipe{false, true}); - sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun); - sim.run(1.0, 0.005); + sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05*arb::units::ms), fun); + sim.run(1.0*arb::units::ms, 0.005*arb::units::ms); um_s_type exp_soma{{ 0, -65 }, { 0.05, -60 }, @@ -228,8 +228,8 @@ TEST(v_process, clamp_fine) { auto rec = v_proc_recipe{true, false}; rec.gprop.default_parameters.discretization = arb::cv_policy_max_extent(0.5); auto sim = arb::simulation(rec); - sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05), fun); - sim.run(1.0, 0.005); + sim.add_sampler(arb::all_probes, arb::regular_schedule(0.05*arb::units::ms), fun); + sim.run(1.0*arb::units::ms, 0.005*arb::units::ms); um_s_type exp_soma{{ 0, -65 }, { 0.05, -42 },