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

Clock.epsilon too small? #730

Closed
thesamovar opened this Issue Jul 11, 2016 · 21 comments

Comments

Projects
None yet
2 participants
@thesamovar
Member

thesamovar commented Jul 11, 2016

The following was reported on the mailing list:

When I restore a network model from last run and run the simulation, I encountered such error, ValueError: Cannot set dt from 100. us to 10. us, the time 112. s is not a multiple of 10. Us.

This error will not happen if the last simulation time is very short such as 10. S. Don’t know where is the problem.

I just checked and (112/10e-6)*10e-6 - 112 = 1.4e-14, and Clock.epsilon=1e-14. At the moment we have an absolute check for whether a time is a multiple of dt, i.e. abs(t_new-t_old)<Clock.epsilon and so this would be failing this check.

Two options, either we increase Clock.epsilon or we use a relative measure. Any thoughts?

@thesamovar thesamovar added the bug label Jul 11, 2016

@thesamovar thesamovar added this to the 2.0 milestone Jul 11, 2016

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 18, 2016

Member

We do need something more robust here (and possibly in other places), I think just lowering the absolute epsilon is not a good idea. How about using numpy.allclose with both an absolute and a relative error (http://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.allclose.html#numpy.allclose)?

Member

mstimberg commented Jul 18, 2016

We do need something more robust here (and possibly in other places), I think just lowering the absolute epsilon is not a good idea. How about using numpy.allclose with both an absolute and a relative error (http://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.allclose.html#numpy.allclose)?

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 19, 2016

Member

For background: I had a feeling we discussed this before, so I went back to check when this was added (in #318) and it seems that we didn't discuss this point in detail.

Yes, I think having an absolute and relative error check makes sense.

One option would be to look at the fractional part of (i_old*dt_old)/dt_new and raise an error if it is too far from 0. So for the example above we get a value of 1.5e-8. If we increase the time to be a day long simulation at a dt of 100us changing to 10us we get a value of 1.5e-5. Given that seems to be on the extreme side, I think I'd be happy with a cutoff of 1e-4 here, i.e. we'd tolerate an error of up to 0.01% of a timestep. That seems large enough to catch the extreme case and small enough that it wouldn't introduce much error. What do you think?

Member

thesamovar commented Jul 19, 2016

For background: I had a feeling we discussed this before, so I went back to check when this was added (in #318) and it seems that we didn't discuss this point in detail.

Yes, I think having an absolute and relative error check makes sense.

One option would be to look at the fractional part of (i_old*dt_old)/dt_new and raise an error if it is too far from 0. So for the example above we get a value of 1.5e-8. If we increase the time to be a day long simulation at a dt of 100us changing to 10us we get a value of 1.5e-5. Given that seems to be on the extreme side, I think I'd be happy with a cutoff of 1e-4 here, i.e. we'd tolerate an error of up to 0.01% of a timestep. That seems large enough to catch the extreme case and small enough that it wouldn't introduce much error. What do you think?

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 25, 2016

Member

In general I agree, but just from the point of documenting it, the "fractional part of (i_old*dt_old)/dt_new" is the same as abs(t_new - t_old)/dt_new, no? Other than that, 0.01% of a timestep sounds very reasonable.

Member

mstimberg commented Jul 25, 2016

In general I agree, but just from the point of documenting it, the "fractional part of (i_old*dt_old)/dt_new" is the same as abs(t_new - t_old)/dt_new, no? Other than that, 0.01% of a timestep sounds very reasonable.

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 25, 2016

Member

From the point of view of documenting it, yes, but I think the calculation is better without the subtraction? Maybe it's also better to document the actual calculation too? Not sure.

I'm also still wondering if there is a more principled / correct calculation we could do. Seems like there could be given that we have the integer and the two dts so in principle we could avoid the imprecision of going to a large float and then converting back to an int.

Member

thesamovar commented Jul 25, 2016

From the point of view of documenting it, yes, but I think the calculation is better without the subtraction? Maybe it's also better to document the actual calculation too? Not sure.

I'm also still wondering if there is a more principled / correct calculation we could do. Seems like there could be given that we have the integer and the two dts so in principle we could avoid the imprecision of going to a large float and then converting back to an int.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 25, 2016

Member

I'm also still wondering if there is a more principled / correct calculation we could do. Seems like there could be given that we have the integer and the two dts so in principle we could avoid the imprecision of going to a large float and then converting back to an int.

Yes, I have been wondering the same -- let's still think this over a bit before committing to a solution.

Member

mstimberg commented Jul 25, 2016

I'm also still wondering if there is a more principled / correct calculation we could do. Seems like there could be given that we have the integer and the two dts so in principle we could avoid the imprecision of going to a large float and then converting back to an int.

Yes, I have been wondering the same -- let's still think this over a bit before committing to a solution.

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 25, 2016

Member

Actually I think the calculation should be fine. With double precision, you have a 53 bit mantissa. That means you can store all integers up to 2^53 exactly. Now in order to make our check that the fractional part is less than 0.01% (or 1e-4) we'd need a resolution of approximately 1e4 times better than that, i.e. we could have the result of (i_old*dt_old)/dt_new be as large as 2**53*1e-4. If dt was as small as 1us this is about 10 days. Given that, I think it's OK. What do you think? That reasoning is a bit rough and maybe could do with some further thought, but suggests we're OK.

Member

thesamovar commented Jul 25, 2016

Actually I think the calculation should be fine. With double precision, you have a 53 bit mantissa. That means you can store all integers up to 2^53 exactly. Now in order to make our check that the fractional part is less than 0.01% (or 1e-4) we'd need a resolution of approximately 1e4 times better than that, i.e. we could have the result of (i_old*dt_old)/dt_new be as large as 2**53*1e-4. If dt was as small as 1us this is about 10 days. Given that, I think it's OK. What do you think? That reasoning is a bit rough and maybe could do with some further thought, but suggests we're OK.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 25, 2016

Member

That seems to make sense but I'll try to look into this in more detail tomorrow. A general thought: this calculation should go into a separate function that we can then easily document and -- more importantly -- test.

Member

mstimberg commented Jul 25, 2016

That seems to make sense but I'll try to look into this in more detail tomorrow. A general thought: this calculation should go into a separate function that we can then easily document and -- more importantly -- test.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 26, 2016

Member

I realized why I have a bit of an issue with our discussion on this check -- I think we are mixing two separate checks here:

  1. Is the current time a multiple of the new dt value? E.g. we don't want to allow a change to a dt of 0.5ms if the current t is 1.2ms.
  2. Is our time t too big to represent steps of dt with sufficient accuracy? E.g. for dt=1*nsecond, a t of about 18min is fine:
>>> diff((2**40 + arange(10)) * 1*nsecond - 2**40 * 1*nsecond)
array([ 0.99998942,  0.99998942,  0.99998942,  0.99998942,  0.99998942,
        0.99998942,  0.99998942,  0.99998942,  0.99998942]) * nsecond

but a t of about 13d is probably not:

>>> diff((2**50 + arange(10)) * 1*nsecond - 2**50 * 1*nsecond)
array([ 0.93132257,  1.16415322,  0.93132257,  0.93132257,  0.93132257,
        1.16415322,  0.93132257,  0.93132257,  1.16415322]) * nsecond

I think the check we have was made with 1) in mind but if I am not mistaken, your argument is also about 2). The problem is that we can (and should) check for 1) when changing dt but we this is not the point to check for 2) because it not only applies when changing dt. So I think we might need two checks, one that checks for the dt compatibility (in set_t_update_dt, the one we currently have) and one that checks for too long simulation times (in set_interval). We might consider making the latter check a warning instead of an error, because the user might not actually care about the time t being exact (and you might have a run(1e9*second) with an external stop or something like that).

Does that make sense?

Member

mstimberg commented Jul 26, 2016

I realized why I have a bit of an issue with our discussion on this check -- I think we are mixing two separate checks here:

  1. Is the current time a multiple of the new dt value? E.g. we don't want to allow a change to a dt of 0.5ms if the current t is 1.2ms.
  2. Is our time t too big to represent steps of dt with sufficient accuracy? E.g. for dt=1*nsecond, a t of about 18min is fine:
>>> diff((2**40 + arange(10)) * 1*nsecond - 2**40 * 1*nsecond)
array([ 0.99998942,  0.99998942,  0.99998942,  0.99998942,  0.99998942,
        0.99998942,  0.99998942,  0.99998942,  0.99998942]) * nsecond

but a t of about 13d is probably not:

>>> diff((2**50 + arange(10)) * 1*nsecond - 2**50 * 1*nsecond)
array([ 0.93132257,  1.16415322,  0.93132257,  0.93132257,  0.93132257,
        1.16415322,  0.93132257,  0.93132257,  1.16415322]) * nsecond

I think the check we have was made with 1) in mind but if I am not mistaken, your argument is also about 2). The problem is that we can (and should) check for 1) when changing dt but we this is not the point to check for 2) because it not only applies when changing dt. So I think we might need two checks, one that checks for the dt compatibility (in set_t_update_dt, the one we currently have) and one that checks for too long simulation times (in set_interval). We might consider making the latter check a warning instead of an error, because the user might not actually care about the time t being exact (and you might have a run(1e9*second) with an external stop or something like that).

Does that make sense?

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 26, 2016

Member

It's not quite so simple, because we could allow a change of dt from 0.5ms to 0.75ms if we can run for t=1.5ms (old i=3, new i=2). We could disallow this but it seems a bit arbitrary to do so. But you're definitely right that we ought to have explicit warnings for the case where t is too big relative to dt (e.g. could warn when i is more than 2**49 on the assumption that we want it to be accurate to within 0.01% of a timestep, as I suggested earlier).

Member

thesamovar commented Jul 26, 2016

It's not quite so simple, because we could allow a change of dt from 0.5ms to 0.75ms if we can run for t=1.5ms (old i=3, new i=2). We could disallow this but it seems a bit arbitrary to do so. But you're definitely right that we ought to have explicit warnings for the case where t is too big relative to dt (e.g. could warn when i is more than 2**49 on the assumption that we want it to be accurate to within 0.01% of a timestep, as I suggested earlier).

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 26, 2016

Member

It's not quite so simple, because we could allow a change of dt from 0.5ms to 0.75ms if we can run for t=1.5ms (old i=3, new i=2). We could disallow this but it seems a bit arbitrary to do so.

Not sure I understand, this passes test number 1 ("Is the current time a multiple of the new dt value?"), doesn't it?

Member

mstimberg commented Jul 26, 2016

It's not quite so simple, because we could allow a change of dt from 0.5ms to 0.75ms if we can run for t=1.5ms (old i=3, new i=2). We could disallow this but it seems a bit arbitrary to do so.

Not sure I understand, this passes test number 1 ("Is the current time a multiple of the new dt value?"), doesn't it?

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 26, 2016

Member

Oops, sorry I misunderstood: I thought you meant timestep not time!

Member

thesamovar commented Jul 26, 2016

Oops, sorry I misunderstood: I thought you meant timestep not time!

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 26, 2016

Member

OK so yes, I agree. We want two checks indeed, one for the validity of dt, and the other for the size of i. We could optionally allow both to be warnings instead of errors? Maybe a preference whether or not an incompatible dt raises an error or just a warning?

Member

thesamovar commented Jul 26, 2016

OK so yes, I agree. We want two checks indeed, one for the validity of dt, and the other for the size of i. We could optionally allow both to be warnings instead of errors? Maybe a preference whether or not an incompatible dt raises an error or just a warning?

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 26, 2016

Member

Default being error of course.

Member

thesamovar commented Jul 26, 2016

Default being error of course.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 26, 2016

Member

I think that the first check (validity of dt) should be an error, because it gives wrong results and the user might have made a mistake (e.g. they wanted to start at t 0s again but instead they continued a previous simulation). For the second check I think it could well be a warning, it is not really a user mistake but a limitation of the time representation in Brian. If for whatever reason they want to do this crazily long simulation with small dt's than they don't have any alternative. The warning would just make sure they are aware of the numerical problems.

Now what would be our criterion for the dt validity check? abs(t_new - t_old)/dt_new < 1e-4? I think it would be straightforward to explain and should never happen for normal use cases (in the report triggering this issue, that value would be 1.4e-9).

In general, I'd prefer not to use a preference for this, I am worried it would be the kind of preference that is never used and does more harm than good.

Member

mstimberg commented Jul 26, 2016

I think that the first check (validity of dt) should be an error, because it gives wrong results and the user might have made a mistake (e.g. they wanted to start at t 0s again but instead they continued a previous simulation). For the second check I think it could well be a warning, it is not really a user mistake but a limitation of the time representation in Brian. If for whatever reason they want to do this crazily long simulation with small dt's than they don't have any alternative. The warning would just make sure they are aware of the numerical problems.

Now what would be our criterion for the dt validity check? abs(t_new - t_old)/dt_new < 1e-4? I think it would be straightforward to explain and should never happen for normal use cases (in the report triggering this issue, that value would be 1.4e-9).

In general, I'd prefer not to use a preference for this, I am worried it would be the kind of preference that is never used and does more harm than good.

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 26, 2016

Member

I'm still not sure what the best expression is. I tried out two:

f1 = lambda i,dt_old,dt_new: abs(int(round((i*dt_old)/dt_new))*dt_new-i*dt_old)/dt_new
f2 = lambda i,dt_old,dt_new: abs(i*dt_old/dt_new-int(round((i*dt_old)/dt_new)))

using the following test code:

from numpy import *
from random import randint
N = 1000000
imax = 2**40
err1 = 0
err2 = 0
dt_old = 100e-6
dt_new = 10e-6
for _ in xrange(N):
    i = randint(0, imax)
    err1 += f1(i, dt_old, dt_new)
    err2 += f2(i, dt_old, dt_new)
print err1/N, err2/N

For most values of imax, the average error is lower for f2 than f1 (if imax=2**50 the average errors are .14 and .03 for f1 and f2). However, for dt_new=50e-6 f1 always gets an error of 0 whereas f2 gets nonzero errors, so it's not entirely straightforward. In any case, for most reasonable values of i (e.g. less than 2**40) the errors are typically tiny and similar for both f1 and f2. For example, for imax=2**30 the average errors are 1.8e-7 and 1.3e-7 so it doesn't make much difference, and it's only when i starts to get close to 2**40 that they start producing errors that are larger than the 1e-4 cutoff we've been talking about. So now we'd be talking about simulation with 2**40 timesteps, which is pretty implausible. (A single Python pass takes about 40ns and so even doing that 2**40 times would take half a day to run.)

In conclusion, I don't think it makes much difference which validity check expression we use for most reasonable ranges. The fractional part one (f2) might be slightly more accurate at the extremes, but there are also some relative values of dt where your expression (f1) gives no error at all and the fractional part one gives some error. So, go with whichever you prefer I think. :)

OK for not having a preference btw.

Member

thesamovar commented Jul 26, 2016

I'm still not sure what the best expression is. I tried out two:

f1 = lambda i,dt_old,dt_new: abs(int(round((i*dt_old)/dt_new))*dt_new-i*dt_old)/dt_new
f2 = lambda i,dt_old,dt_new: abs(i*dt_old/dt_new-int(round((i*dt_old)/dt_new)))

using the following test code:

from numpy import *
from random import randint
N = 1000000
imax = 2**40
err1 = 0
err2 = 0
dt_old = 100e-6
dt_new = 10e-6
for _ in xrange(N):
    i = randint(0, imax)
    err1 += f1(i, dt_old, dt_new)
    err2 += f2(i, dt_old, dt_new)
print err1/N, err2/N

For most values of imax, the average error is lower for f2 than f1 (if imax=2**50 the average errors are .14 and .03 for f1 and f2). However, for dt_new=50e-6 f1 always gets an error of 0 whereas f2 gets nonzero errors, so it's not entirely straightforward. In any case, for most reasonable values of i (e.g. less than 2**40) the errors are typically tiny and similar for both f1 and f2. For example, for imax=2**30 the average errors are 1.8e-7 and 1.3e-7 so it doesn't make much difference, and it's only when i starts to get close to 2**40 that they start producing errors that are larger than the 1e-4 cutoff we've been talking about. So now we'd be talking about simulation with 2**40 timesteps, which is pretty implausible. (A single Python pass takes about 40ns and so even doing that 2**40 times would take half a day to run.)

In conclusion, I don't think it makes much difference which validity check expression we use for most reasonable ranges. The fractional part one (f2) might be slightly more accurate at the extremes, but there are also some relative values of dt where your expression (f1) gives no error at all and the fractional part one gives some error. So, go with whichever you prefer I think. :)

OK for not having a preference btw.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 27, 2016

Member

I'd say the first measure, because it measures the quantity relevant for practical purposes: the difference between the time t in the old and the new representation. If we use your example of dt_old=100e-6 and dt_new=50e-6, then you can see the difference between the two measures at time step 13. Since 100e-6 does not have a finite binary representation, even 13*dt_old/dt_old will not return 13 exactly. The measure that f2 uses will give the difference between 13*dt_old/dt_new ( = 26.000000000000004) and the value 26. The difference is tiny of course but there is a difference. However, 13*dt_old and 26*dt_new, i.e. the time where the simulation stopped and where it will start are exactly the same, so everything is fine and this is why f1 returns 0. I'd therefore go for this measure, I feel it is more straightforward (and as you said, it will probably not make any difference in practice).

So to summarize the current proposal:

  • Let's use f1 and a cutoff of 1e-4 (raising an error if the check fails)
  • In set_interval, check that the end time can still correctly represent dt steps with no more than 1e-4 error (and raise a warning if not). Do we need a clever check for that or is abs((i_end + 1)*dt - i_end*dt - dt) < 1e-4 enough?
Member

mstimberg commented Jul 27, 2016

I'd say the first measure, because it measures the quantity relevant for practical purposes: the difference between the time t in the old and the new representation. If we use your example of dt_old=100e-6 and dt_new=50e-6, then you can see the difference between the two measures at time step 13. Since 100e-6 does not have a finite binary representation, even 13*dt_old/dt_old will not return 13 exactly. The measure that f2 uses will give the difference between 13*dt_old/dt_new ( = 26.000000000000004) and the value 26. The difference is tiny of course but there is a difference. However, 13*dt_old and 26*dt_new, i.e. the time where the simulation stopped and where it will start are exactly the same, so everything is fine and this is why f1 returns 0. I'd therefore go for this measure, I feel it is more straightforward (and as you said, it will probably not make any difference in practice).

So to summarize the current proposal:

  • Let's use f1 and a cutoff of 1e-4 (raising an error if the check fails)
  • In set_interval, check that the end time can still correctly represent dt steps with no more than 1e-4 error (and raise a warning if not). Do we need a clever check for that or is abs((i_end + 1)*dt - i_end*dt - dt) < 1e-4 enough?
@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 27, 2016

Member

OK!

Is it enough to just check that i_end is smaller than 2**40 or something like that? I'm not sure that your expression works because of effects like, e.g., 1e50-1e-50==1e50. Also, shouldn't it be abs((i_end + 1)*dt - i_end*dt - dt)/dt < 1e-4? Otherwise if dt<1e-4 we'd be allowing an error of more than dt.

Just quickly checking abs((i_end + 1)*dt - i_end*dt - dt)/dt for large values of i_end it does seem that it reaches 1e-4 at around 2**43 to 2**44 for a wide range of values of dt (which makes sense because the whole expression is things multiplied by dt, and float accuracy is scale invariant across a wide range). So I'd say just checking i_end<2**40 is probably enough. It's rare that it'll be larger than that, and probably if it is people need a warning anyway!

Member

thesamovar commented Jul 27, 2016

OK!

Is it enough to just check that i_end is smaller than 2**40 or something like that? I'm not sure that your expression works because of effects like, e.g., 1e50-1e-50==1e50. Also, shouldn't it be abs((i_end + 1)*dt - i_end*dt - dt)/dt < 1e-4? Otherwise if dt<1e-4 we'd be allowing an error of more than dt.

Just quickly checking abs((i_end + 1)*dt - i_end*dt - dt)/dt for large values of i_end it does seem that it reaches 1e-4 at around 2**43 to 2**44 for a wide range of values of dt (which makes sense because the whole expression is things multiplied by dt, and float accuracy is scale invariant across a wide range). So I'd say just checking i_end<2**40 is probably enough. It's rare that it'll be larger than that, and probably if it is people need a warning anyway!

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 27, 2016

Member

Is it enough to just check that i_end is smaller than 2**40 or something like that? I'm not sure that your expression works because of effects like, e.g., 1e50-1e-50==1e50.

The /dt (or 1e-4*dt) was indeed a mistake. But otherwise than that the expression should work, I think.

So I'd say just checking i_end<2**40 is probably enough. It's rare that it'll be larger than that, and probably if it is people need a warning anyway!

Right, I'd be fine with a simple check like that, it is certainly not something where it is worth complicating things.

That said, there are two more situations where we are using Clock.epsilon and I think it would be good to deal with them in a more principled way now (or maybe get rid of them). We have two occurrences of the following code:

        new_i = np.uint64(np.round(target_t/new_dt))
        new_t = new_i*new_dt
        if (new_t == target_t or
                    np.abs(new_t-target_t) <= self.epsilon*np.abs(new_t)):
            new_timestep = new_i
        else:
            new_timestep = np.uint64(np.ceil(target_t/new_dt))

Once when setting a new t (should only happen when restoring a network) and once for calculating the end time step. Do we have examples where we need the "approximately equal" comparison? For sure there is no such example in our test suite, everything still passes when I remove it...

The second situation where we use Clock.epsilon is when we calculate the "current clocks":

        minclock, min_time = min(clocks_times, key=lambda k: k[1])
        curclocks = set(clock for clock, time in clocks_times if
                        (time == min_time or
                         abs(time - min_time) < Clock.epsilon))

Again, do we really need the "approximately equal" check here? I think it is important that we deal with multiple clocks that have the same dt here (which can happen if you define dt=... manually for objects), but in this case their times will be exactly the same. I think we put this in to deal with multiple clocks that have dt's that are multiples of each others to make sure that they are recognized as such. Practically, though, I don't think e.g. dt * i and ((dt/2) * (i*2)) ever differ until we get to those crazily high timestep values we talked about earlier.

Member

mstimberg commented Jul 27, 2016

Is it enough to just check that i_end is smaller than 2**40 or something like that? I'm not sure that your expression works because of effects like, e.g., 1e50-1e-50==1e50.

The /dt (or 1e-4*dt) was indeed a mistake. But otherwise than that the expression should work, I think.

So I'd say just checking i_end<2**40 is probably enough. It's rare that it'll be larger than that, and probably if it is people need a warning anyway!

Right, I'd be fine with a simple check like that, it is certainly not something where it is worth complicating things.

That said, there are two more situations where we are using Clock.epsilon and I think it would be good to deal with them in a more principled way now (or maybe get rid of them). We have two occurrences of the following code:

        new_i = np.uint64(np.round(target_t/new_dt))
        new_t = new_i*new_dt
        if (new_t == target_t or
                    np.abs(new_t-target_t) <= self.epsilon*np.abs(new_t)):
            new_timestep = new_i
        else:
            new_timestep = np.uint64(np.ceil(target_t/new_dt))

Once when setting a new t (should only happen when restoring a network) and once for calculating the end time step. Do we have examples where we need the "approximately equal" comparison? For sure there is no such example in our test suite, everything still passes when I remove it...

The second situation where we use Clock.epsilon is when we calculate the "current clocks":

        minclock, min_time = min(clocks_times, key=lambda k: k[1])
        curclocks = set(clock for clock, time in clocks_times if
                        (time == min_time or
                         abs(time - min_time) < Clock.epsilon))

Again, do we really need the "approximately equal" check here? I think it is important that we deal with multiple clocks that have the same dt here (which can happen if you define dt=... manually for objects), but in this case their times will be exactly the same. I think we put this in to deal with multiple clocks that have dt's that are multiples of each others to make sure that they are recognized as such. Practically, though, I don't think e.g. dt * i and ((dt/2) * (i*2)) ever differ until we get to those crazily high timestep values we talked about earlier.

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 27, 2016

Member

For the first one, I agree that if the only case we change t is when restoring networks an exact equality check is fine. I guess that's a bit of old code from when we allowed changing t.

For the second one, I wonder what's the best solution. I agree that halving timesteps will be exact, but we know it's not exact for other multiples such as the original issue: 11200000*10e-6-1120000*100e-6 is approximately 1.4e-14. Of course, using Clock.epsilon doesn't work here either. What about we use abs(i_1*dt_1-i_2*dt_2)/min(dt_1, dt_2)<1e-4 as the condition? This means that the times are within 0.01% of a dt from the point of view of both clocks. For this example the LHS of the expression is 1e-9.

Member

thesamovar commented Jul 27, 2016

For the first one, I agree that if the only case we change t is when restoring networks an exact equality check is fine. I guess that's a bit of old code from when we allowed changing t.

For the second one, I wonder what's the best solution. I agree that halving timesteps will be exact, but we know it's not exact for other multiples such as the original issue: 11200000*10e-6-1120000*100e-6 is approximately 1.4e-14. Of course, using Clock.epsilon doesn't work here either. What about we use abs(i_1*dt_1-i_2*dt_2)/min(dt_1, dt_2)<1e-4 as the condition? This means that the times are within 0.01% of a dt from the point of view of both clocks. For this example the LHS of the expression is 1e-9.

@mstimberg

This comment has been minimized.

Show comment
Hide comment
@mstimberg

mstimberg Jul 28, 2016

Member

For the first one, I agree that if the only case we change t is when restoring networks an exact equality check is fine. I guess that's a bit of old code from when we allowed changing t.

Right, I'll remove that. There's identical code for the end of a simulation as well, though. I don't really see what could go wrong there, the comment states "Sets the start and end value of the clock precisely if possible (using epsilon) or rounding up if not. This assures that multiple calls to Network.run will not re-run the same time step.". We have tests that do several runs with and without changing dt's but none fails when I remove this comparison. I'd say unless we can come up with a test case, let's remove it.

What about we use abs(i_1_dt_1-i_2_dt_2)/min(dt_1, dt_2)<1e-4 as the condition? This means that the times are within 0.01% of a dt from the point of view of both clocks. For this example the LHS of the expression is 1e-9.

Ok, let's do that. It's doing a bit more computation per time step but the more important optimisation would be to special case the "one clock" case where we don't have to do any of this, anyway.

I'll try to address all these point in a branch and then we can discuss the result (or any problems that came up during the implementation).

Member

mstimberg commented Jul 28, 2016

For the first one, I agree that if the only case we change t is when restoring networks an exact equality check is fine. I guess that's a bit of old code from when we allowed changing t.

Right, I'll remove that. There's identical code for the end of a simulation as well, though. I don't really see what could go wrong there, the comment states "Sets the start and end value of the clock precisely if possible (using epsilon) or rounding up if not. This assures that multiple calls to Network.run will not re-run the same time step.". We have tests that do several runs with and without changing dt's but none fails when I remove this comparison. I'd say unless we can come up with a test case, let's remove it.

What about we use abs(i_1_dt_1-i_2_dt_2)/min(dt_1, dt_2)<1e-4 as the condition? This means that the times are within 0.01% of a dt from the point of view of both clocks. For this example the LHS of the expression is 1e-9.

Ok, let's do that. It's doing a bit more computation per time step but the more important optimisation would be to special case the "one clock" case where we don't have to do any of this, anyway.

I'll try to address all these point in a branch and then we can discuss the result (or any problems that came up during the implementation).

@thesamovar

This comment has been minimized.

Show comment
Hide comment
@thesamovar

thesamovar Jul 28, 2016

Member

Ah I missed the bit about end of simulation. I think in that case we should continue to use the approximate check, but instead use abs(target_t-new_t)/new_dt<1e-4. Also, let's make all these 1e-4 be the new Clock.epsilon (or Clock.rel_epsilon if you prefer?). Example case where target_t!=new_t is the original example. In this case, the two paths give the same result, but that's just chance (50% of the time it would be different).

I think the reason it doesn't come up in our tests is all our tests are short runs with everything being nice, small integer multiples of dt.

Member

thesamovar commented Jul 28, 2016

Ah I missed the bit about end of simulation. I think in that case we should continue to use the approximate check, but instead use abs(target_t-new_t)/new_dt<1e-4. Also, let's make all these 1e-4 be the new Clock.epsilon (or Clock.rel_epsilon if you prefer?). Example case where target_t!=new_t is the original example. In this case, the two paths give the same result, but that's just chance (50% of the time it would be different).

I think the reason it doesn't come up in our tests is all our tests are short runs with everything being nice, small integer multiples of dt.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment