From 6a5bdc4fb471d7b547c1fad84a38cb0ac88b1aac Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Tue, 2 Apr 2024 15:14:51 +0200 Subject: [PATCH 1/6] Tiny fix for |epoch| < dt. --- arbor/simulation.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index 95963e85fd..b806a1a547 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -355,7 +355,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: // @@ -399,6 +399,8 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { if (tfinal<=epoch_.t1) return epoch_.t1; + auto epoch_length = std::max(t_interval_, dt); + // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { epoch next = e; @@ -485,8 +487,8 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { }; epoch prev = epoch_; - epoch current = next_epoch(prev, t_interval_); - epoch next = next_epoch(current, t_interval_); + epoch current = next_epoch(prev, epoch_length); + epoch next = next_epoch(current, epoch_length); if (epoch_callback_) epoch_callback_(current.t0, tfinal); @@ -507,7 +509,7 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { for (;;) { prev = current; current = next; - next = next_epoch(next, t_interval_); + next = next_epoch(next, epoch_length); if (next.empty()) break; g.run([&]() { exchange(prev); enqueue(next); }); From 4fb547691925034b0798e4ba46b2ddcec69137f1 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 3 Apr 2024 10:47:04 +0200 Subject: [PATCH 2/6] Tighten the checks on simulation::run. --- arbor/connection.hpp | 2 -- arbor/include/arbor/recipe.hpp | 6 +++--- arbor/simulation.cpp | 19 ++++++++++--------- 3 files changed, 13 insertions(+), 14 deletions(-) diff --git a/arbor/connection.hpp b/arbor/connection.hpp index e83fa28722..b002ee5194 100644 --- a/arbor/connection.hpp +++ b/arbor/connection.hpp @@ -1,7 +1,5 @@ #pragma once -#include - #include #include #include diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index dbcba5fcde..84f832acf6 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -54,8 +54,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::isnan(weight)) throw std::domain_error("Connection weight must be finite."); + if (std::isnan(delay) || delay <= 0) throw std::domain_error("Connection delay must be positive, finite, and given in units of [ms]."); } }; @@ -69,7 +69,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::isnan(weight)) throw std::domain_error("Gap junction weight must be finite."); } }; diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index b806a1a547..a905980058 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -19,6 +19,7 @@ #include "threading/threading.hpp" #include "util/maputil.hpp" #include "util/span.hpp" +#include "util/strprintf.hpp" #include "profile/profiler_macro.hpp" namespace arb { @@ -394,12 +395,12 @@ 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(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; - - auto epoch_length = std::max(t_interval_, dt); + // Compare up to picoseconds + time_type eps = 1e-9; + if (!std::isfinite(tfinal) || tfinal < eps) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); + if (epoch_.t1 - tfinal > eps) throw std::domain_error(util::pprintf("simulation: tfinal={}ms is in the past, current time of simulation is {}ms", tfinal, epoch_.t1)); + if (!std::isfinite(dt) || dt < eps) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); + if (dt - t_interval_ > eps) throw std::domain_error(util::pprintf("simulation: dt={}ms is larger than epoch length by {}, chose at most half the minimal connection delay {}ms.", dt, dt - t_interval_, t_interval_)); // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { @@ -487,8 +488,8 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { }; epoch prev = epoch_; - epoch current = next_epoch(prev, epoch_length); - epoch next = next_epoch(current, epoch_length); + epoch current = next_epoch(prev, t_interval_); + epoch next = next_epoch(current, t_interval_); if (epoch_callback_) epoch_callback_(current.t0, tfinal); @@ -509,7 +510,7 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { for (;;) { prev = current; current = next; - next = next_epoch(next, epoch_length); + next = next_epoch(next, t_interval_); if (next.empty()) break; g.run([&]() { exchange(prev); enqueue(next); }); From 12261a242d76d269bf6993ff935a2cf00b75c458 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 3 Apr 2024 11:34:02 +0200 Subject: [PATCH 3/6] Stricter checks, adjust tests. --- arbor/simulation.cpp | 4 ++-- test/unit/test_lif_cell_group.cpp | 2 +- test/unit/test_serdes.cpp | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index a905980058..da842c5119 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -397,10 +397,10 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { // Compare up to picoseconds time_type eps = 1e-9; - if (!std::isfinite(tfinal) || tfinal < eps) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); - if (epoch_.t1 - tfinal > eps) throw std::domain_error(util::pprintf("simulation: tfinal={}ms is in the past, current time of simulation is {}ms", tfinal, epoch_.t1)); if (!std::isfinite(dt) || dt < eps) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); if (dt - t_interval_ > eps) throw std::domain_error(util::pprintf("simulation: dt={}ms is larger than epoch length by {}, chose at most half the minimal connection delay {}ms.", dt, dt - t_interval_, t_interval_)); + if (!std::isfinite(tfinal) || tfinal < eps) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); + if (tfinal - epoch_.t1 < dt) throw std::domain_error(util::pprintf("simulation: tfinal={}ms doesn't make progress of least one dt; current time of simulation is {}ms and dt {}ms", tfinal, epoch_.t1, dt)); // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp index 8e13b222fa..03b55fb742 100644 --- a/test/unit/test_lif_cell_group.cpp +++ b/test/unit/test_lif_cell_group.cpp @@ -708,7 +708,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 From 2739392c0a2e7627268273c40e336adaa983d188 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Wed, 3 Apr 2024 13:31:19 +0200 Subject: [PATCH 4/6] Loosen checks a tiny bit, ensure min_delay is >0. --- arbor/include/arbor/recipe.hpp | 6 +-- arbor/simulation.cpp | 12 ++++-- python/test/unit/test_simulation.py | 66 +++++++++++++++++++++++++++++ 3 files changed, 77 insertions(+), 7 deletions(-) create mode 100644 python/test/unit/test_simulation.py diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp index 84f832acf6..bd0b495abb 100644 --- a/arbor/include/arbor/recipe.hpp +++ b/arbor/include/arbor/recipe.hpp @@ -54,8 +54,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::domain_error("Connection weight must be finite."); - if (std::isnan(delay) || delay <= 0) throw std::domain_error("Connection delay must be positive, finite, and given 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]."); } }; @@ -69,7 +69,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::domain_error("Gap junction weight must be finite."); + if (!std::isfinite(weight)) throw std::domain_error("Gap junction weight must be finite."); } }; diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index da842c5119..90ee9f3beb 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -111,7 +111,11 @@ class simulation_state { void inject_events(const cse_vector& events); - 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_; @@ -398,9 +402,9 @@ time_type simulation_state::run(time_type tfinal, time_type dt) { // Compare up to picoseconds time_type eps = 1e-9; if (!std::isfinite(dt) || dt < eps) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); - if (dt - t_interval_ > eps) throw std::domain_error(util::pprintf("simulation: dt={}ms is larger than epoch length by {}, chose at most half the minimal connection delay {}ms.", dt, dt - t_interval_, t_interval_)); - if (!std::isfinite(tfinal) || tfinal < eps) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); - if (tfinal - epoch_.t1 < dt) throw std::domain_error(util::pprintf("simulation: tfinal={}ms doesn't make progress of least one dt; current time of simulation is {}ms and dt {}ms", tfinal, epoch_.t1, dt)); + if (!std::isfinite(tfinal) || tfinal - epoch_.t1 < eps) { + throw std::domain_error("simulation: tfinal must be finite, positive, larger than the current time, and in [ms]"); + } // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { diff --git a/python/test/unit/test_simulation.py b/python/test/unit/test_simulation.py new file mode 100644 index 0000000000..50b3d3e2f4 --- /dev/null +++ b/python/test/unit/test_simulation.py @@ -0,0 +1,66 @@ +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 test_cannot_relive_the_past(self): + T = 1 * U.ms + dt = 0.01 * U.ms + rec = DelayRecipe(2*dt) + sim = A.simulation(rec) + sim.run(T + dt, dt) + self.assertRaises(ValueError, sim.run, T, dt) From 4797b9ffaf6a3ca875388866a1262c592f7d99d8 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Fri, 20 Sep 2024 07:54:43 +0200 Subject: [PATCH 5/6] Fix tests. --- arbor/simulation.cpp | 6 ++---- python/test/unit/test_simulation.py | 8 +++++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp index f41bbcae96..359de8fa1f 100644 --- a/arbor/simulation.cpp +++ b/arbor/simulation.cpp @@ -398,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. - // Compare up to picoseconds - time_type eps = 1e-9; - if (!std::isfinite(dt) || dt < eps) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]"); - if (!std::isfinite(tfinal) || tfinal < eps) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]"); + 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]"); // Compute following epoch, with max time tfinal. auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch { diff --git a/python/test/unit/test_simulation.py b/python/test/unit/test_simulation.py index 50b3d3e2f4..73755ba7c5 100644 --- a/python/test/unit/test_simulation.py +++ b/python/test/unit/test_simulation.py @@ -45,6 +45,7 @@ def probes(self, _): def global_properties(self, _): return A.neuron_cable_properties() + class TestDelayNetwork(unittest.TestCase): def test_zero_delay(self): rec = DelayRecipe(0.0 * U.ms) @@ -57,10 +58,11 @@ def test_dt_half_delay(self): sim = A.simulation(rec) sim.run(T, dt) - def test_cannot_relive_the_past(self): + def dt_must_be_finite(self): T = 1 * U.ms dt = 0.01 * U.ms rec = DelayRecipe(2*dt) sim = A.simulation(rec) - sim.run(T + dt, dt) - self.assertRaises(ValueError, sim.run, T, dt) + 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) From 73acec6685f5be5c468e308ff293108e51733b91 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Fri, 20 Sep 2024 07:55:39 +0200 Subject: [PATCH 6/6] Format. --- python/test/unit/test_simulation.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/python/test/unit/test_simulation.py b/python/test/unit/test_simulation.py index 73755ba7c5..764c28e5d7 100644 --- a/python/test/unit/test_simulation.py +++ b/python/test/unit/test_simulation.py @@ -6,10 +6,12 @@ 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 @@ -23,24 +25,23 @@ def cell_kind(self, _): 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) + 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'))) + 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)] + return [A.connection((src, "det"), "syn", 0.42, self.delay)] def probes(self, _): - return [A.cable_probe_membrane_voltage('(location 0 0.5)', 'Um')] + return [A.cable_probe_membrane_voltage("(location 0 0.5)", "Um")] def global_properties(self, _): return A.neuron_cable_properties() @@ -54,15 +55,15 @@ def test_zero_delay(self): def test_dt_half_delay(self): T = 1 * U.ms dt = 0.01 * U.ms - rec = DelayRecipe(2*dt) + 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) + 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) + self.assertRaises(ValueError, sim.run, T, 1.0 / 0.0 * U.ms) + self.assertRaises(ValueError, sim.run, 1.0 / 0.0 * U.ms, dt)