Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Stricter checks on network construction / simulation parameters. #2264

Draft
wants to merge 8 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 0 additions & 2 deletions arbor/connection.hpp
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
#pragma once

#include <cstdint>

#include <arbor/common_types.hpp>
#include <arbor/spike.hpp>
#include <arbor/spike_event.hpp>
Expand Down
6 changes: 3 additions & 3 deletions arbor/include/arbor/recipe.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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.");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Message doesn't match test: should this test be !std::isfinite(weight)?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's better, yes.

if (std::isnan(delay) || delay <= 0) throw std::domain_error("Connection delay must be positive, finite, and given in units of [ms].");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inasmuch as we're presuming NaNs work as expected, we could just have the test !(delay > 0) instead of std::isnan(delay) || delay <= 0.

Do we really need to enforce that delay is finite? If so, then the test should include that.

Also, not being familiar (yet) with how LLNL units works, why do we need to specify that the quantity is in milliseconds? Can't we just convert as required or else assert in the type?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, the two issues go hand in hand: (42 * U.ms).value_as(U.mV) == NaN. So, receiving a nan can mean either we got nan * U.ms or an erroneous unit. This is why the message allows for both.

}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have quite specific exceptions such as arb::bad_connection_source_gid defined in <arbexcept.hpp>; it would be consistent to define some exceptions to throw here that derive from arb::arbor_exception.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am a bit torn here, std::domain_error fits the basic checks for isnan etc well and has the advantage that it plays nice with pybind11 (it'll automatically convert to ValueError).

};

Expand All @@ -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.");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar comments as above: testing finite or testing Nan? We should use an arbor exception.

}
};

Expand Down
13 changes: 8 additions & 5 deletions arbor/simulation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -355,7 +356,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:
//
Expand Down Expand Up @@ -394,10 +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;
// 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));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems all a bit fiddly; as mentioned in the general comments, I think we should just leave dt interpretation up to the integrators and they can make a sensible choice, e.g. clip it by epoch duration or interpret very small dt values as being larger, based on whatever they have to do. Having eps here is untidy because it splits the responsibility for sane dt interpretation between the main loop and the integrators.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, if someone wants to run arbor for ever, do we really need to stop them? We already have a check for tfinal being less than the end time of the preceding epoch, where we return the correct 'simulated up to' time, so we don't really need tests for it being zero, or negative; in fact a tfinal of zero should be a valid no-op in my opinion.

In short, a isnan test for tfinal remains sufficient, I believe.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In endeffect you are suggesting that instead of throwing an error, the cell groups' individual advance methods should decide. Do enable that, we'd have to pass all possible options, currently min_delay and dt, down the stack. I'd rather make the choice uniformly, especially since #2053.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No we don't have to pass it down the stack: we know epochs are at most min_delay/2 long. The epoch already has all the info they need, and if there is more info it will be in the cell kind global data, which they also have.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#2053 is just cable cells, and I still want to partially revert it so that we can use flexible time steps with a different cable cell integrator, even if it's just to make them line up correctly with epoch intervals.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As an example, if we have the last epoch (ending on tfinal) having duration 3.1 dt, wouldn't it be best to set the fixed dt for that epoch to be e.g. the duration/4 and then we would know that all cell states across cell groups were actually at the same integration time?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, if someone wants to run arbor for ever, do we really need to stop them? We already have a check for tfinal being less than the end time of the preceding epoch, where we return the correct 'simulated up to' time, so we don't really need tests for it being zero, or negative; in fact a tfinal of zero should be a valid no-op in my opinion.

Technically I agree. Practically I know that I'd spent way too much time looking for the reason why my simulation did nothing in that situation. Especially since our time_type doesn't distinguish time points from durations, so run(5 *ms) really could mean run until t=5ms or run for 5ms. Thus I think it's friendlier towards the user to tell them 'this is a no-op and you likely didn't mean that'. Especially, as I cannot currently imagine situations where this is a semantically meaningful request.

Copy link
Contributor

@halfflat halfflat Apr 3, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As a general principle, I think we should keep things general unless there is a good reason not to. I mean, the parameter is literally called tfinal, and we return the simulation state time - it's a simple interface with simple semantics (from the outside) and there's nothing stopping us adding checks on the outside of that if we want to provide more hand holds in the Python interface (though I don't think we should there either).

We can always change the name to run_to.

The semantics of running a simulation from t = 0 for 0 seconds should be that the state reflects the initial conditions. It's not what we'd expect someone to do in normal circumstances, but it's the consistent result, and may arise in circumstances where the simulator is being driven by another process or co-simulation.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Damn you and your arguments :D


// Compute following epoch, with max time tfinal.
auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch {
Expand Down
2 changes: 1 addition & 1 deletion test/unit/test_lif_cell_group.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -708,7 +708,7 @@ TEST(lif_cell_group, probe_with_connections) {
[&spikes](const std::vector<spike>& 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<Um_type> exp = {{ 0, -18 },
{ 0.025, -17.9750624 },
{ 0.05, -17.9502492 },
Expand Down
4 changes: 2 additions & 2 deletions test/unit/test_serdes.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading