diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp index d68c0c73e2..3ff0dd0ef8 100644 --- a/arbor/include/arbor/cable_cell_param.hpp +++ b/arbor/include/arbor/cable_cell_param.hpp @@ -66,9 +66,8 @@ struct ARB_SYMBOL_VISIBLE i_clamp { const U::quantity& current): t(time.value_as(U::ms)), amplitude(current.value_as(U::nA)) { - - if (std::isnan(t)) throw std::domain_error{"Time must be finite and convertible to ms."}; - if (std::isnan(amplitude)) throw std::domain_error{"Amplitude must be finite and convertible to nA."}; + if (!std::isfinite(t)) throw std::domain_error{"Time must be finite and convertible to ms."}; + if (!std::isfinite(amplitude)) throw std::domain_error{"Amplitude must be finite and convertible to nA."}; } double t; // [ms] double amplitude; // [nA] @@ -104,8 +103,8 @@ struct ARB_SYMBOL_VISIBLE i_clamp { frequency(f.value_as(U::kHz)), phase(phi.value_as(U::rad)) { - if (std::isnan(frequency)) throw std::domain_error{"Frequency must be finite and convertible to kHz."}; - if (std::isnan(phase)) throw std::domain_error{"Phase must be finite and convertible to rad."}; + if (!std::isfinite(frequency)) throw std::domain_error{"Frequency must be finite and convertible to kHz."}; + if (!std::isfinite(phase)) throw std::domain_error{"Phase must be finite and convertible to rad."}; } // A 'box' stimulus with fixed onset time, duration, and constant amplitude. @@ -123,7 +122,7 @@ struct ARB_SYMBOL_VISIBLE i_clamp { // Threshold detector description. struct ARB_SYMBOL_VISIBLE threshold_detector { threshold_detector(const U::quantity& m): threshold(m.value_as(U::mV)) { - if (std::isnan(threshold)) throw std::domain_error{"Threshold must be finite and in [mV]."}; + if (!std::isfinite(threshold)) throw std::domain_error{"Threshold must be finite and in [mV]."}; } static threshold_detector from_raw_millivolts(double v) { return {v*U::mV}; } double threshold; // [mV] @@ -139,7 +138,7 @@ struct ARB_SYMBOL_VISIBLE init_membrane_potential { init_membrane_potential() = default; init_membrane_potential(const U::quantity& m, iexpr scale=1): value(m.value_as(U::mV)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mV]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mV]."}; } }; @@ -151,7 +150,7 @@ struct ARB_SYMBOL_VISIBLE temperature { temperature() = default; temperature(const U::quantity& m, iexpr scale=1): value(m.value_as(U::Kelvin)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [K]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [K]."}; } }; @@ -162,7 +161,7 @@ struct ARB_SYMBOL_VISIBLE axial_resistivity { axial_resistivity() = default; axial_resistivity(const U::quantity& m, iexpr scale=1): value(m.value_as(U::cm*U::Ohm)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [Ω·cm]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [Ω·cm]."}; } }; @@ -173,7 +172,7 @@ struct ARB_SYMBOL_VISIBLE membrane_capacitance { membrane_capacitance() = default; membrane_capacitance(const U::quantity& m, iexpr scale=1): value(m.value_as(U::F/U::m2)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [F/m²]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [F/m²]."}; } }; @@ -185,7 +184,7 @@ struct ARB_SYMBOL_VISIBLE init_int_concentration { init_int_concentration() = default; init_int_concentration(const std::string& ion, const U::quantity& m, iexpr scale=1): ion{ion}, value(m.value_as(U::mM)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mM]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mM]."}; } }; @@ -197,7 +196,7 @@ struct ARB_SYMBOL_VISIBLE ion_diffusivity { ion_diffusivity() = default; ion_diffusivity(const std::string& ion, const U::quantity& m, iexpr scale=1): ion{ion}, value(m.value_as(U::m2/U::s)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [m²/s]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [m²/s]."}; } }; @@ -209,7 +208,7 @@ struct ARB_SYMBOL_VISIBLE init_ext_concentration { init_ext_concentration() = default; init_ext_concentration(const std::string& ion, const U::quantity& m, iexpr scale=1): ion{ion}, value(m.value_as(U::mM)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mM]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mM]."}; } }; @@ -221,7 +220,7 @@ struct ARB_SYMBOL_VISIBLE init_reversal_potential { init_reversal_potential() = default; init_reversal_potential(const std::string& ion, const U::quantity& m, iexpr scale=1): ion{ion}, value(m.value_as(U::mV)), scale{scale} { - if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mV]."}; + if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mV]."}; } }; @@ -461,13 +460,13 @@ struct ARB_SYMBOL_VISIBLE cable_cell_global_properties { auto &ion_data = default_parameters.ion_data[ion_name]; ion_data.init_int_concentration = init_iconc.value_as(U::mM); - if (std::isnan(*ion_data.init_int_concentration)) throw std::domain_error("init_int_concentration must be finite and convertible to mM"); + if (!std::isfinite(*ion_data.init_int_concentration)) throw std::domain_error("init_int_concentration must be finite and convertible to mM"); ion_data.init_ext_concentration = init_econc.value_as(U::mM); - if (std::isnan(*ion_data.init_ext_concentration)) throw std::domain_error("init_ext_concentration must be finite and convertible to mM"); + if (!std::isfinite(*ion_data.init_ext_concentration)) throw std::domain_error("init_ext_concentration must be finite and convertible to mM"); ion_data.init_reversal_potential = init_revpot.value_as(U::mV); - if (std::isnan(*ion_data.init_reversal_potential)) throw std::domain_error("init_reversal_potential must be finite and convertible to mV"); + if (!std::isfinite(*ion_data.init_reversal_potential)) throw std::domain_error("init_reversal_potential must be finite and convertible to mV"); ion_data.diffusivity = diffusivity.value_as(U::m2/U::s); - if (std::isnan(*ion_data.diffusivity) || *ion_data.diffusivity < 0) throw std::domain_error("diffusivity must be positive, finite, and convertible to m2/s"); + if (!std::isfinite(*ion_data.diffusivity) || *ion_data.diffusivity < 0) throw std::domain_error("diffusivity must be positive, finite, and convertible to m2/s"); } void add_ion(const std::string& ion_name, diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index a8bd8eb17e..39daadb55a 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -58,8 +58,8 @@ struct cell_connection_base { cell_connection_base(L src, cell_local_label_type dst, float w, const U::quantity& d): source(std::move(src)), target(std::move(dst)), weight(w), delay(d.value_as(U::ms)) { - if (std::isnan(weight)) throw std::out_of_range("Connection weight must be finite."); - if (std::isnan(delay) || delay < 0) throw std::out_of_range("Connection delay must be non-negative and infinite in units of [ms]."); + if (!std::isfinite(weight)) throw std::domain_error("Connection weight must be finite."); + if (!std::isfinite(delay) || delay <= 0) throw std::domain_error("Connection delay must be positive, finite, and given in units of [ms]."); } }; @@ -73,7 +73,7 @@ struct gap_junction_connection { gap_junction_connection(cell_global_label_type peer, cell_local_label_type local, double g): peer(std::move(peer)), local(std::move(local)), weight(g) { - if (std::isnan(weight)) throw std::out_of_range("Gap junction weight must be finite."); + if (!std::isfinite(weight)) throw std::domain_error("Gap junction weight must be finite."); } }; diff --git a/arbor/lif_cell_group.hpp b/arbor/lif_cell_group.hpp index ff922334ad..4cb1bc1546 100644 --- a/arbor/lif_cell_group.hpp +++ b/arbor/lif_cell_group.hpp @@ -43,13 +43,13 @@ struct ARB_SYMBOL_VISIBLE lif_lowered_cell { V_m = lif.V_m.value_as(U::mV); t_ref = lif.t_ref.value_as(U::ms); - if (std::isnan(V_th)) throw std::out_of_range("V_th must be finite and in [mV]"); - if (std::isnan(tau_m) || tau_m < 0) throw std::out_of_range("tau_m must be positive, finite, and in [ms]"); - if (std::isnan(C_m) || C_m < 0) throw std::out_of_range("C_m must be positive, finite, and in [pF]"); - if (std::isnan(E_L)) throw std::out_of_range("E_L must be finite and in [mV]"); - if (std::isnan(E_R)) throw std::out_of_range("E_R must be finite and in [mV]"); - if (std::isnan(V_m)) throw std::out_of_range("V_m must be finite and in [mV]"); - if (std::isnan(t_ref) || t_ref < 0) throw std::out_of_range("t_ref must be positive, finite, and in [ms]"); + if (!std::isfinite(V_th)) throw std::domain_error("V_th must be finite and in [mV]"); + if (!std::isfinite(tau_m) || tau_m < 0) throw std::domain_error("tau_m must be positive, finite, and in [ms]"); + if (!std::isfinite(C_m) || C_m < 0) throw std::domain_error("C_m must be positive, finite, and in [pF]"); + if (!std::isfinite(E_L)) throw std::domain_error("E_L must be finite and in [mV]"); + if (!std::isfinite(E_R)) throw std::domain_error("E_R must be finite and in [mV]"); + if (!std::isfinite(V_m)) throw std::domain_error("V_m must be finite and in [mV]"); + if (!std::isfinite(t_ref) || t_ref < 0) throw std::domain_error("t_ref must be positive, finite, and in [ms]"); } ARB_SERDES_ENABLE(lif_lowered_cell, source, target, tau_m, V_th, C_m, E_L, E_R, V_m, t_ref); diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 998963790b..359de8fa1f 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -20,6 +20,7 @@ #include "threading/threading.hpp" #include "util/maputil.hpp" #include "util/span.hpp" +#include "util/strprintf.hpp" #include "profile/profiler_macro.hpp" namespace arb { @@ -109,7 +110,11 @@ class simulation_state { void set_remote_spike_filter(const spike_predicate& p) { return communicator_.set_remote_spike_filter(p); } - time_type min_delay() { return communicator_.min_delay(); } + time_type min_delay() { + auto tau = communicator_.min_delay(); + if (tau <= 0.0 || !std::isfinite(tau)) throw std::domain_error("Minimum connection delay must be strictly positive and finite."); + return tau; + } spike_export_function global_export_callback_; spike_export_function local_export_callback_; @@ -354,7 +359,7 @@ void simulation_state::reset() { time_type simulation_state::run(time_type tfinal, time_type dt) { // Progress simulation to time tfinal, through a series of integration epochs // of length at most t_interval_. t_interval_ is chosen to be no more than - // than half the network minimum delay. + // than half the network minimum delay and minimally the timestep `dt`. // // There are three simulation tasks that can be run partially in parallel: // @@ -393,10 +398,8 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { // Requires state at end of run(), with epoch_.id==k: // * U(k) and D(k) have completed. + if (!std::isfinite(dt) || dt <= 0) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); if (!std::isfinite(tfinal) || tfinal < 0) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); - if (!std::isfinite(dt) || tfinal < 0) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); - - if (tfinal<=epoch_.t1) return epoch_.t1; // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { @@ -580,9 +583,9 @@ void simulation::update(const recipe& rec) { impl_->update(rec); } time_type simulation::run(const units::quantity& tfinal, const units::quantity& dt) { auto dt_ms = dt.value_as(units::ms); - if (dt_ms <= 0.0 || std::isnan(dt_ms)) throw domain_error("Finite time-step must be supplied."); + if (dt_ms <= 0.0 || !std::isfinite(dt_ms)) throw domain_error("Finite time-step must be supplied."); auto tfinal_ms = tfinal.value_as(units::ms); - if (tfinal_ms <= 0.0 || std::isnan(tfinal_ms)) throw domain_error("Finite time-step must be supplied."); + if (tfinal_ms <= 0.0 || !std::isfinite(tfinal_ms)) throw domain_error("Finite time-step must be supplied."); return impl_->run(tfinal_ms, dt_ms); } diff --git a/example/brunel/brunel.cpp b/example/brunel/brunel.cpp index bfcfca9d0b..15f0f89d55 100644 --- a/example/brunel/brunel.cpp +++ b/example/brunel/brunel.cpp @@ -96,7 +96,7 @@ class brunel_recipe: public recipe { ncells_exc_(nexc), ncells_inh_(ninh), delay_(delay), seed_(seed) { // Make sure that in_degree_prop in the interval (0, 1] if (in_degree_prop <= 0.0 || in_degree_prop > 1.0) { - throw std::out_of_range("The proportion of incoming connections should be in the interval (0, 1]."); + throw std::domain_error("The proportion of incoming connections should be in the interval (0, 1]."); } // Set up the parameters. diff --git a/python/test/unit/test_simulation.py b/python/test/unit/test_simulation.py new file mode 100644 index 0000000000..764c28e5d7 --- /dev/null +++ b/python/test/unit/test_simulation.py @@ -0,0 +1,69 @@ +import arbor as A +from arbor import units as U +import unittest + +""" +End to end tests at simulation level. +""" + + +class DelayRecipe(A.recipe): + """ + Construct a simple network with configurable axonal delay. + """ + + def __init__(self, delay): + A.recipe.__init__(self) + self.delay = delay + + def num_cells(self): + return 2 + + def cell_kind(self, _): + return A.cell_kind.cable + + def cell_description(self, _): + # morphology + tree = A.segment_tree() + tree.append(A.mnpos, A.mpoint(-1, 0, 0, 1), A.mpoint(1, 0, 0, 1), tag=1) + + decor = ( + A.decor() + .place("(location 0 0.5)", A.synapse("expsyn"), "syn") + .place("(location 0 0.5)", A.threshold_detector(-15 * U.mV), "det") + .paint("(all)", A.density("hh")) + ) + + return A.cable_cell(tree, decor, A.label_dict()) + + def connections_on(self, gid): + src = (gid + 1) % self.num_cells() + return [A.connection((src, "det"), "syn", 0.42, self.delay)] + + def probes(self, _): + return [A.cable_probe_membrane_voltage("(location 0 0.5)", "Um")] + + def global_properties(self, _): + return A.neuron_cable_properties() + + +class TestDelayNetwork(unittest.TestCase): + def test_zero_delay(self): + rec = DelayRecipe(0.0 * U.ms) + self.assertRaises(ValueError, A.simulation, rec) + + def test_dt_half_delay(self): + T = 1 * U.ms + dt = 0.01 * U.ms + rec = DelayRecipe(2 * dt) + sim = A.simulation(rec) + sim.run(T, dt) + + def dt_must_be_finite(self): + T = 1 * U.ms + dt = 0.01 * U.ms + rec = DelayRecipe(2 * dt) + sim = A.simulation(rec) + self.assertRaises(ValueError, sim.run, T, 0 * U.ms) + self.assertRaises(ValueError, sim.run, T, 1.0 / 0.0 * U.ms) + self.assertRaises(ValueError, sim.run, 1.0 / 0.0 * U.ms, dt) diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp index 72d79d858a..b8aba8e5af 100644 --- a/test/unit/test_lif_cell_group.cpp +++ b/test/unit/test_lif_cell_group.cpp @@ -704,7 +704,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*U::ms, 0.005*U::ms); + sim.run(10*U::ms, 0.0025*U::ms); std::vector exp = {{ 0, -18 }, { 0.025, -17.9750624 }, { 0.05, -17.9502492 }, diff --git a/test/unit/test_serdes.cpp b/test/unit/test_serdes.cpp index 546362daac..1855720934 100644 --- a/test/unit/test_serdes.cpp +++ b/test/unit/test_serdes.cpp @@ -233,7 +233,7 @@ TEST(serdes, single_cell) { } TEST(serdes, network) { - auto dt = 0.5*arb::units::ms; + auto dt = 0.05*arb::units::ms; auto T = 5*arb::units::ms; // Result @@ -302,7 +302,7 @@ TEST(serdes, host_device_arrays) { } TEST(serdes, single_cell_gpu) { - auto dt = 0.5*arb::units::ms; + auto dt = 0.05*arb::units::ms; auto T = 5*arb::units::ms; // Result