Skip to content

Commit b2f1851

Browse files
thorstenhaterthorstenhater
authored
Stricter checks on network construction / simulation parameters. (#2264)
Clean-up semantic checks on parameters, ensure connections have finite delay. Closes #2263
1 parent c4ff08c commit b2f1851

17 files changed

+339
-244
lines changed

‎arbor/include/arbor/cable_cell_param.hpp

+17-18
Original file line numberDiff line numberDiff line change
@@ -66,9 +66,8 @@ struct ARB_SYMBOL_VISIBLE i_clamp {
6666
const U::quantity& current):
6767
t(time.value_as(U::ms)),
6868
amplitude(current.value_as(U::nA)) {
69-
70-
if (std::isnan(t)) throw std::domain_error{"Time must be finite and convertible to ms."};
71-
if (std::isnan(amplitude)) throw std::domain_error{"Amplitude must be finite and convertible to nA."};
69+
if (!std::isfinite(t)) throw std::domain_error{"Time must be finite and convertible to ms."};
70+
if (!std::isfinite(amplitude)) throw std::domain_error{"Amplitude must be finite and convertible to nA."};
7271
}
7372
double t; // [ms]
7473
double amplitude; // [nA]
@@ -104,8 +103,8 @@ struct ARB_SYMBOL_VISIBLE i_clamp {
104103
frequency(f.value_as(U::kHz)),
105104
phase(phi.value_as(U::rad))
106105
{
107-
if (std::isnan(frequency)) throw std::domain_error{"Frequency must be finite and convertible to kHz."};
108-
if (std::isnan(phase)) throw std::domain_error{"Phase must be finite and convertible to rad."};
106+
if (!std::isfinite(frequency)) throw std::domain_error{"Frequency must be finite and convertible to kHz."};
107+
if (!std::isfinite(phase)) throw std::domain_error{"Phase must be finite and convertible to rad."};
109108
}
110109

111110
// A 'box' stimulus with fixed onset time, duration, and constant amplitude.
@@ -123,7 +122,7 @@ struct ARB_SYMBOL_VISIBLE i_clamp {
123122
// Threshold detector description.
124123
struct ARB_SYMBOL_VISIBLE threshold_detector {
125124
threshold_detector(const U::quantity& m): threshold(m.value_as(U::mV)) {
126-
if (std::isnan(threshold)) throw std::domain_error{"Threshold must be finite and in [mV]."};
125+
if (!std::isfinite(threshold)) throw std::domain_error{"Threshold must be finite and in [mV]."};
127126
}
128127
static threshold_detector from_raw_millivolts(double v) { return {v*U::mV}; }
129128
double threshold; // [mV]
@@ -139,7 +138,7 @@ struct ARB_SYMBOL_VISIBLE init_membrane_potential {
139138
init_membrane_potential() = default;
140139
init_membrane_potential(const U::quantity& m, iexpr scale=1):
141140
value(m.value_as(U::mV)), scale{scale} {
142-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mV]."};
141+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mV]."};
143142
}
144143
};
145144

@@ -151,7 +150,7 @@ struct ARB_SYMBOL_VISIBLE temperature {
151150
temperature() = default;
152151
temperature(const U::quantity& m, iexpr scale=1):
153152
value(m.value_as(U::Kelvin)), scale{scale} {
154-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [K]."};
153+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [K]."};
155154
}
156155
};
157156

@@ -162,7 +161,7 @@ struct ARB_SYMBOL_VISIBLE axial_resistivity {
162161
axial_resistivity() = default;
163162
axial_resistivity(const U::quantity& m, iexpr scale=1):
164163
value(m.value_as(U::cm*U::Ohm)), scale{scale} {
165-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [Ω·cm]."};
164+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [Ω·cm]."};
166165
}
167166
};
168167

@@ -173,7 +172,7 @@ struct ARB_SYMBOL_VISIBLE membrane_capacitance {
173172
membrane_capacitance() = default;
174173
membrane_capacitance(const U::quantity& m, iexpr scale=1):
175174
value(m.value_as(U::F/U::m2)), scale{scale} {
176-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [F/m²]."};
175+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [F/m²]."};
177176
}
178177
};
179178

@@ -185,7 +184,7 @@ struct ARB_SYMBOL_VISIBLE init_int_concentration {
185184
init_int_concentration() = default;
186185
init_int_concentration(const std::string& ion, const U::quantity& m, iexpr scale=1):
187186
ion{ion}, value(m.value_as(U::mM)), scale{scale} {
188-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mM]."};
187+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mM]."};
189188
}
190189
};
191190

@@ -197,7 +196,7 @@ struct ARB_SYMBOL_VISIBLE ion_diffusivity {
197196
ion_diffusivity() = default;
198197
ion_diffusivity(const std::string& ion, const U::quantity& m, iexpr scale=1):
199198
ion{ion}, value(m.value_as(U::m2/U::s)), scale{scale} {
200-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [m²/s]."};
199+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [m²/s]."};
201200
}
202201
};
203202

@@ -209,7 +208,7 @@ struct ARB_SYMBOL_VISIBLE init_ext_concentration {
209208
init_ext_concentration() = default;
210209
init_ext_concentration(const std::string& ion, const U::quantity& m, iexpr scale=1):
211210
ion{ion}, value(m.value_as(U::mM)), scale{scale} {
212-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mM]."};
211+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mM]."};
213212
}
214213
};
215214

@@ -221,7 +220,7 @@ struct ARB_SYMBOL_VISIBLE init_reversal_potential {
221220
init_reversal_potential() = default;
222221
init_reversal_potential(const std::string& ion, const U::quantity& m, iexpr scale=1):
223222
ion{ion}, value(m.value_as(U::mV)), scale{scale} {
224-
if (std::isnan(value)) throw std::domain_error{"Value must be finite and in [mV]."};
223+
if (!std::isfinite(value)) throw std::domain_error{"Value must be finite and in [mV]."};
225224
}
226225
};
227226

@@ -460,13 +459,13 @@ struct ARB_SYMBOL_VISIBLE cable_cell_global_properties {
460459

461460
auto &ion_data = default_parameters.ion_data[ion_name];
462461
ion_data.init_int_concentration = init_iconc.value_as(U::mM);
463-
if (std::isnan(*ion_data.init_int_concentration)) throw std::domain_error("init_int_concentration must be finite and convertible to mM");
462+
if (!std::isfinite(*ion_data.init_int_concentration)) throw std::domain_error("init_int_concentration must be finite and convertible to mM");
464463
ion_data.init_ext_concentration = init_econc.value_as(U::mM);
465-
if (std::isnan(*ion_data.init_ext_concentration)) throw std::domain_error("init_ext_concentration must be finite and convertible to mM");
464+
if (!std::isfinite(*ion_data.init_ext_concentration)) throw std::domain_error("init_ext_concentration must be finite and convertible to mM");
466465
ion_data.init_reversal_potential = init_revpot.value_as(U::mV);
467-
if (std::isnan(*ion_data.init_reversal_potential)) throw std::domain_error("init_reversal_potential must be finite and convertible to mV");
466+
if (!std::isfinite(*ion_data.init_reversal_potential)) throw std::domain_error("init_reversal_potential must be finite and convertible to mV");
468467
ion_data.diffusivity = diffusivity.value_as(U::m2/U::s);
469-
if (std::isnan(*ion_data.diffusivity) || *ion_data.diffusivity < 0) throw std::domain_error("diffusivity must be positive, finite, and convertible to m2/s");
468+
if (!std::isfinite(*ion_data.diffusivity) || *ion_data.diffusivity < 0) throw std::domain_error("diffusivity must be positive, finite, and convertible to m2/s");
470469
}
471470

472471
void add_ion(const std::string& ion_name,

‎arbor/include/arbor/profile/profiler.hpp

+1-2
Original file line numberDiff line numberDiff line change
@@ -31,8 +31,7 @@ struct profile {
3131
double wall_time;
3232
};
3333

34-
// TODO: remove declaration and update the docs
35-
void profiler_clear();
34+
ARB_ARBOR_API void profiler_clear();
3635
ARB_ARBOR_API void profiler_initialize(context ctx);
3736
ARB_ARBOR_API void profiler_enter(std::size_t region_id);
3837
ARB_ARBOR_API void profiler_leave();

‎arbor/include/arbor/recipe.hpp

+3-3
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,8 @@ struct cell_connection_base {
5858

5959
cell_connection_base(L src, cell_local_label_type dst, float w, const U::quantity& d):
6060
source(std::move(src)), target(std::move(dst)), weight(w), delay(d.value_as(U::ms)) {
61-
if (std::isnan(weight)) throw std::out_of_range("Connection weight must be finite.");
62-
if (std::isnan(delay) || delay < 0) throw std::out_of_range("Connection delay must be non-negative and infinite in units of [ms].");
61+
if (!std::isfinite(weight)) throw std::domain_error("Connection weight must be finite.");
62+
if (!std::isfinite(delay) || delay <= 0) throw std::domain_error("Connection delay must be positive, finite, and given in units of [ms].");
6363
}
6464
};
6565

@@ -73,7 +73,7 @@ struct gap_junction_connection {
7373

7474
gap_junction_connection(cell_global_label_type peer, cell_local_label_type local, double g):
7575
peer(std::move(peer)), local(std::move(local)), weight(g) {
76-
if (std::isnan(weight)) throw std::out_of_range("Gap junction weight must be finite.");
76+
if (!std::isfinite(weight)) throw std::domain_error("Gap junction weight must be finite.");
7777
}
7878
};
7979

‎arbor/lif_cell_group.hpp

+7-7
Original file line numberDiff line numberDiff line change
@@ -43,13 +43,13 @@ struct ARB_SYMBOL_VISIBLE lif_lowered_cell {
4343
V_m = lif.V_m.value_as(U::mV);
4444
t_ref = lif.t_ref.value_as(U::ms);
4545

46-
if (std::isnan(V_th)) throw std::out_of_range("V_th must be finite and in [mV]");
47-
if (std::isnan(tau_m) || tau_m < 0) throw std::out_of_range("tau_m must be positive, finite, and in [ms]");
48-
if (std::isnan(C_m) || C_m < 0) throw std::out_of_range("C_m must be positive, finite, and in [pF]");
49-
if (std::isnan(E_L)) throw std::out_of_range("E_L must be finite and in [mV]");
50-
if (std::isnan(E_R)) throw std::out_of_range("E_R must be finite and in [mV]");
51-
if (std::isnan(V_m)) throw std::out_of_range("V_m must be finite and in [mV]");
52-
if (std::isnan(t_ref) || t_ref < 0) throw std::out_of_range("t_ref must be positive, finite, and in [ms]");
46+
if (!std::isfinite(V_th)) throw std::domain_error("V_th must be finite and in [mV]");
47+
if (!std::isfinite(tau_m) || tau_m < 0) throw std::domain_error("tau_m must be positive, finite, and in [ms]");
48+
if (!std::isfinite(C_m) || C_m < 0) throw std::domain_error("C_m must be positive, finite, and in [pF]");
49+
if (!std::isfinite(E_L)) throw std::domain_error("E_L must be finite and in [mV]");
50+
if (!std::isfinite(E_R)) throw std::domain_error("E_R must be finite and in [mV]");
51+
if (!std::isfinite(V_m)) throw std::domain_error("V_m must be finite and in [mV]");
52+
if (!std::isfinite(t_ref) || t_ref < 0) throw std::domain_error("t_ref must be positive, finite, and in [ms]");
5353
}
5454

5555
ARB_SERDES_ENABLE(lif_lowered_cell, source, target, tau_m, V_th, C_m, E_L, E_R, V_m, t_ref);

‎arbor/profile/profiler.cpp

+12
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,12 @@ class profiler {
118118
static profiler p;
119119
return p;
120120
}
121+
122+
void clear() {
123+
for (auto& r: recorders_) r.clear();
124+
name_index_.clear();
125+
region_names_.clear();
126+
}
121127
};
122128

123129

@@ -384,6 +390,11 @@ ARB_ARBOR_API void profiler_leave() {
384390
profiler::get_global_profiler().leave();
385391
}
386392

393+
ARB_ARBOR_API void profiler_clear() {
394+
profiler::get_global_profiler().clear();
395+
}
396+
397+
387398
ARB_ARBOR_API region_id_type profiler_region_id(const std::string& name) {
388399
if (!is_valid_region_string(name)) {
389400
throw std::runtime_error(std::string("'")+name+"' is not a valid profiler region name.");
@@ -419,6 +430,7 @@ ARB_ARBOR_API std::ostream& print_profiler_summary(std::ostream& os, double limi
419430

420431
#else
421432

433+
ARB_ARBOR_API void profiler_clear() {}
422434
ARB_ARBOR_API void profiler_leave() {}
423435
ARB_ARBOR_API void profiler_enter(region_id_type) {}
424436
ARB_ARBOR_API profile profiler_summary();

‎arbor/simulation.cpp

+22-7
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#include "threading/threading.hpp"
2020
#include "util/maputil.hpp"
2121
#include "util/span.hpp"
22+
#include "util/strprintf.hpp"
2223
#include "profile/profiler_macro.hpp"
2324

2425
namespace arb {
@@ -112,7 +113,11 @@ class simulation_state {
112113

113114
void set_remote_spike_filter(const spike_predicate& p) { return communicator_.set_remote_spike_filter(p); }
114115

115-
time_type min_delay() { return communicator_.min_delay(); }
116+
time_type min_delay() {
117+
auto tau = communicator_.min_delay();
118+
if (tau <= 0.0 || !std::isfinite(tau)) throw std::domain_error("Minimum connection delay must be strictly positive and finite.");
119+
return tau;
120+
}
116121

117122
spike_export_function global_export_callback_;
118123
spike_export_function local_export_callback_;
@@ -357,7 +362,7 @@ void simulation_state::reset() {
357362
time_type simulation_state::run(time_type tfinal, time_type dt) {
358363
// Progress simulation to time tfinal, through a series of integration epochs
359364
// of length at most t_interval_. t_interval_ is chosen to be no more than
360-
// than half the network minimum delay.
365+
// than half the network minimum delay and minimally the timestep `dt`.
361366
//
362367
// There are three simulation tasks that can be run partially in parallel:
363368
//
@@ -396,10 +401,20 @@ time_type simulation_state::run(time_type tfinal, time_type dt) {
396401
// Requires state at end of run(), with epoch_.id==k:
397402
// * U(k) and D(k) have completed.
398403

404+
if (!std::isfinite(dt) || dt <= 0) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]");
399405
if (!std::isfinite(tfinal) || tfinal < 0) throw std::domain_error("simulation: tfinal must be finite, positive, and in [ms]");
400-
if (!std::isfinite(dt) || tfinal < 0) throw std::domain_error("simulation: dt must be finite, positive, and in [ms]");
401-
402-
if (tfinal<=epoch_.t1) return epoch_.t1;
406+
// NOTE It would be preferable to have the cell groups to make smart
407+
// decisions instead. However, since there is a tradeoff since silently
408+
// installing a smaller timestep might result in inacceptable runtimes.
409+
// NOTE we are using tau/2 since that is the actual integration time due to
410+
// the double buffering
411+
// NOTE we are _also_ using float epsilon and some fudging as connections store
412+
// delay as float. This shakes out be roughly 1e-6, equivalent to 1ps.
413+
if (dt - t_interval_ > 10*std::numeric_limits<float>::epsilon()) {
414+
throw arbor_exception(
415+
util::pprintf("simulation: Timestep {}ms is too large, must be less than half the minimum network delay ({}ms)",
416+
dt, t_interval_));
417+
}
403418

404419
// Compute following epoch, with max time tfinal.
405420
auto next_epoch = [tfinal](epoch e, time_type interval) -> epoch {
@@ -583,9 +598,9 @@ void simulation::update(const recipe& rec) { impl_->update(rec); }
583598

584599
time_type simulation::run(const units::quantity& tfinal, const units::quantity& dt) {
585600
auto dt_ms = dt.value_as(units::ms);
586-
if (dt_ms <= 0.0 || std::isnan(dt_ms)) throw domain_error("Finite time-step must be supplied.");
601+
if (dt_ms <= 0.0 || !std::isfinite(dt_ms)) throw domain_error("Finite time-step must be supplied.");
587602
auto tfinal_ms = tfinal.value_as(units::ms);
588-
if (tfinal_ms <= 0.0 || std::isnan(tfinal_ms)) throw domain_error("Finite time-step must be supplied.");
603+
if (tfinal_ms <= 0.0 || !std::isfinite(tfinal_ms)) throw domain_error("Finite time-step must be supplied.");
589604
return impl_->run(tfinal_ms, dt_ms);
590605
}
591606

‎doc/cpp/profiler.rst

+4-2
Original file line numberDiff line numberDiff line change
@@ -193,8 +193,10 @@ and the profiler regions can be reset.
193193
profile::profiler_clear();
194194
}
195195
196-
After a call to ``util::profiler_clear``, all counters and timers are set to zero.
197-
This could be used, for example, to generate separate profiler reports for model building and model execution phases.
196+
After a call to ``util::profiler_clear``, all counters and timers are set to
197+
zero. This could be used, for example, to generate separate profiler reports for
198+
model building and model execution phases. It is also required to restart a
199+
simulation terminated with a caught exception.
198200

199201
Profiler output
200202
~~~~~~~~~~~~~~~

‎example/brunel/brunel.cpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ struct cl_options {
5252
bool use_cc = false;
5353
// Simulation running parameters:
5454
double tfinal = 100.;
55-
double dt = 1;
55+
double dt = 0.05;
5656
uint32_t group_size = 10;
5757
uint32_t seed = 42;
5858
// Parameters for spike output.
@@ -102,7 +102,7 @@ class brunel_recipe: public recipe {
102102
ncells_exc_(nexc), ncells_inh_(ninh), delay_(delay), seed_(seed), use_cable_cells(use_cc) {
103103
// Make sure that in_degree_prop in the interval (0, 1]
104104
if (in_degree_prop <= 0.0 || in_degree_prop > 1.0) {
105-
throw std::out_of_range("The proportion of incoming connections should be in the interval (0, 1].");
105+
throw std::domain_error("The proportion of incoming connections should be in the interval (0, 1].");
106106
}
107107

108108
// Set up the parameters.

‎example/plasticity/plasticity.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,7 @@ void print_header(double from, double to) {
105105
<< "|---------+-------------+----------|\n";
106106
}
107107

108-
const double dt = 0.05;
108+
const double dt = 0.025;
109109

110110
int main(int argc, char** argv) {
111111
auto rec = recipe(3);

‎python/profiler.cpp

+15-11
Original file line numberDiff line numberDiff line change
@@ -54,17 +54,21 @@ void register_profiler(pybind11::module& m) {
5454
.def("__repr__", [](arb::profile::meter_report& r){return "<arbor.meter_report>";});
5555

5656
#ifdef ARB_PROFILE_ENABLED
57-
m.def("profiler_initialize", [](context_shim& ctx) {
58-
arb::profile::profiler_initialize(ctx.context);
59-
});
60-
m.def("profiler_summary",
61-
[](double limit){
62-
std::stringstream stream;
63-
arb::profile::print_profiler_summary(stream, limit);
64-
return stream.str();
65-
},
66-
"limit"_a=0.0,
67-
"Show summary of the profile; printing contributions above `limit` percent. Defaults to showing all.");
57+
m
58+
.def("profiler_initialize", [](context_shim& ctx) {
59+
arb::profile::profiler_initialize(ctx.context);
60+
})
61+
.def("profiler_summary",
62+
[](double limit){
63+
std::stringstream stream;
64+
arb::profile::print_profiler_summary(stream, limit);
65+
return stream.str();
66+
},
67+
"limit"_a=0.0,
68+
"Show summary of the profile; printing contributions above `limit` percent. Defaults to showing all.")
69+
.def("profiler_clear",
70+
[] { arb::profile::profiler_clear(); },
71+
"Reset the profiler.");
6872
#endif
6973
}
7074

0 commit comments

Comments
 (0)