[WIP] Use an epsilon independent of dt when comparing times

brian2/core/clocks.py

@@ -36,8 +36,7 @@ def check_dt(new_dt, old_dt, target_t):
     ------
     ValueError
         If using the new dt value would lead to a difference in the target
-        time of more than `Clock.epsilon_dt` times ``new_dt`` (by default,
-        0.01% of the new dt).
+        time of more than `Clock.epsilon_dt` (by default, 10ns).
 
     Examples
     --------
@@ -52,7 +51,7 @@ def check_dt(new_dt, old_dt, target_t):
     old_t = np.int64(np.round(target_t / old_dt)) * old_dt
     new_t = np.int64(np.round(target_t / new_dt)) * new_dt
     error_t = target_t
-    if abs(new_t - old_t)/new_dt > Clock.epsilon_dt:
+    if abs(new_t - old_t) > Clock.epsilon_dt:
         raise ValueError(('Cannot set dt from {old} to {new}, the '
                           'time {t} is not a multiple of '
                           '{new}').format(old=str(old_dt * second),
@@ -74,10 +73,8 @@ class Clock(VariableOwner):
     Notes
     -----
     Clocks are run in the same `Network.run` iteration if `~Clock.t` is the
-    same. The condition for two
-    clocks to be considered as having the same time is
-    ``abs(t1-t2)<epsilon*abs(t1)``, a standard test for equality of floating
-    point values. The value of ``epsilon`` is ``1e-14``.
+    same. The condition for two clocks to be considered as having the same time
+    is ``abs(t1-t2)<10*ns``.
     '''
 
     def __init__(self, dt, name='clock*'):
@@ -121,7 +118,7 @@ def _set_t_update_dt(self, target_t=0*second):
     def _calc_timestep(self, target_t):
         '''
         Calculate the integer time step for the target time. If it cannot be
-        exactly represented (up to 0.01% of dt), round up.
+        exactly represented (up to 10ns), round up.
 
         Parameters
         ----------
@@ -136,7 +133,7 @@ def _calc_timestep(self, target_t):
         new_i = np.int64(np.round(target_t / self.dt_))
         new_t = new_i * self.dt_
         if (new_t == target_t or
-                        np.abs(new_t - target_t)/self.dt_ <= Clock.epsilon_dt):
+                        np.abs(new_t - target_t) <= Clock.epsilon_dt):
             new_timestep = new_i
         else:
             new_timestep = np.int64(np.ceil(target_t / self.dt_))
@@ -189,9 +186,9 @@ def set_interval(self, start, end):
                                      self._i_end),
                         'many_timesteps')
 
-    #: The relative difference for times (in terms of dt) so that they are
+    #: The absolute difference for times (in seconds) so that clocks are
     #: considered identical.
-    epsilon_dt = 1e-4
+    epsilon_dt = 1e-8
 
 
 class DefaultClockProxy(object):

brian2/core/network.py

@@ -877,7 +877,7 @@ def _nextclocks(self):
         minclock, min_time, minclock_dt = min(clocks_times_dt, key=lambda k: k[1])
         curclocks = set(clock for clock, time, dt in clocks_times_dt if
                         (time == min_time or
-                         abs(time - min_time)/min(minclock_dt, dt) < Clock.epsilon_dt))
+                         abs(time - min_time) < Clock.epsilon_dt))
         return minclock, curclocks
 
     @device_override('network_run')

brian2/devices/cpp_standalone/brianlib/clocks.h

@@ -19,7 +19,7 @@ class Clock
 	double *dt;
 	int64_t *timestep;
 	double *t;
-	Clock(double _epsilon=1e-14) : epsilon(_epsilon) { i_end = 0;};
+	Clock(double _epsilon=1e-8) : epsilon(_epsilon) { i_end = 0;};
     inline void tick()
     {
         timestep[0] += 1;
@@ -30,18 +30,18 @@ class Clock
 	{
         int i_start = fround(start/dt[0]);
         double t_start = i_start*dt[0];
-        if(t_start==start || fabs(t_start-start)<=epsilon*fabs(t_start))
+        if(t_start==start || fabs(t_start-start)<=epsilon)
         {
             timestep[0] = i_start;
         } else
         {
-            timestep[0] = (int)ceil(start/dt[0]);
+            timestep[0] = (int64_t)ceil(start/dt[0]);
         }
         i_end = fround(end/dt[0]);
         double t_end = i_end*dt[0];
-        if(!(t_end==end || fabs(t_end-end)<=epsilon*fabs(t_end)))
+        if(!(t_end==end || fabs(t_end-end)<=epsilon))
         {
-            i_end = (int)ceil(end/dt[0]);
+            i_end = (int64_t)ceil(end/dt[0]);
         }
 	}
 private:

brian2/devices/cpp_standalone/templates/network.cpp

@@ -7,8 +7,6 @@
 #include<utility>
 {{ openmp_pragma('include') }}
 
-#define Clock_epsilon 1e-14
-
 double Network::_last_run_time = 0.0;
 double Network::_last_run_completed_fraction = 0.0;
 
@@ -158,7 +156,7 @@ Clock* Network::next_clocks()
     {
         Clock *clock = *i;
         double s = clock->t[0];
-        if(s==t || fabs(s-t)<=Clock_epsilon)
+        if(s==t || fabs(s-t) <= clock->epsilon)
             curclocks.insert(clock);
     }
     return minclock;

brian2/devices/cpp_standalone/templates/objects.cpp

@@ -59,7 +59,7 @@ SynapticPathway {{path.name}}(
 
 //////////////// clocks ///////////////////
 {% for clock in clocks | sort(attribute='name') %}
-Clock {{clock.name}};  // attributes will be set in run.cpp
+Clock {{clock.name}}({{clock.epsilon_dt}});  // attributes will be set in run.cpp
 {% endfor %}
 
 {% if profiled_codeobjects is defined %}

brian2/tests/test_clocks.py

@@ -59,10 +59,36 @@ def test_set_interval_warning():
     assert logs[0][1].endswith('many_timesteps')
 
 
+@attr('codegen-independent')
+def test_very_long_dt():
+    # See github issue #1054
+    clock1 = Clock(dt=0.1*ms)
+    clock2 = Clock(dt=100000*second)
+    clock1.set_interval(0*ms, 1*ms)
+    clock2.set_interval(0*ms, 1*ms)  # The clock should advance
+    assert clock1.timestep[:] == 0
+    assert clock2.timestep[:] == 0
+    assert clock1._i_end == 10
+    assert clock2._i_end == 1
+    # Simulate advancing the clock
+    clock1.variables['timestep'].set_value(clock1._i_end)
+    clock1.variables['t'].set_value(clock1._i_end * clock1.dt)
+    clock2.variables['timestep'].set_value(clock2._i_end)
+    clock2.variables['t'].set_value(clock2._i_end * clock2.dt)
+
+    clock1.set_interval(1*ms, 2*ms)
+    clock2.set_interval(1*ms, 2*ms)  # The clock should not advance
+    assert clock1.timestep[:] == 10
+    assert clock2.timestep[:] == 1
+    assert clock1._i_end == 20
+    assert clock2._i_end == 1
+
+
 if __name__ == '__main__':
     test_clock_attributes()
     restore_initial_state()
     test_clock_dt_change()
     restore_initial_state()
     test_defaultclock()
     test_set_interval_warning()
+    test_very_long_dt()

brian2/tests/test_network.py

@@ -1372,6 +1372,24 @@ def test_small_runs():
     assert_allclose(mon_1.t_[:], mon_2.t_[:])
     assert_allclose(mon_1.v_[:], mon_2.v_[:])
 
+@attr('standalone-compatible')
+@with_setup(teardown=reinit_devices)
+def test_long_run():
+    defaultclock.dt = 0.1 * ms
+    group = NeuronGroup(1, 'x : 1')
+    group.run_regularly('x += 1')
+    # Timesteps are internally stored as 64bit integers, but previous versions
+    # converted them into 32bit integers along the way. We'll make sure that
+    # this is not the case and everything runs fine. To not actually run such a
+    # long simulation we manually advance the clock.
+    net = Network(group, name='network')
+    start_step = 2**31-5
+    start_time = start_step*defaultclock.dt_
+    with catch_logs() as l:
+        net.t_ = start_time  # for runtime
+        device.insert_code('main', 'network.t = {};'.format(start_time))  # for standalone
+        net.run(6 * defaultclock.dt)
+    assert group.x == 6
 
 @attr('codegen-independent')
 @with_setup(teardown=reinit_devices)
@@ -1409,6 +1427,7 @@ def test_multiple_runs_function_change():
     device.build(direct_call=False, **device.build_options)
     assert_equal(mon.v[0], [1, 2, 3, 4])
 
+
 if __name__ == '__main__':
     BrianLogger.log_level_warn()
     for t in [
@@ -1470,6 +1489,7 @@ def test_multiple_runs_function_change():
             test_magic_scope,
             test_runtime_rounding,
             test_small_runs,
+            test_long_run,
             test_long_run_dt_change,
             test_multiple_runs_constant_change,
             test_multiple_runs_function_change