[WIP] Handling of NaN and other errors in engine - #570

Questions on this (looking only at top-post, no code):

How does this affect the Monte Carlo? What @jchodera wants is not exactly an ignore but a reject: the move is rejected. The approach I thought of for this was: (1) raise a custom exception (maybe EngineNaNError); (2) the MC has a choice of either catching the error and returning NaNErrorMoveChange (which is rejected) or not catching the error.

Importantly, this is not a retry. This considers the NaN move as a valid trial, but rejects it. However, you should be able to easily check your simulation for a large number of moves rejected because of NaN.

Note that you can't just assume that a NaN-containing frame will always reject. Consider:

>>> x = paths.FunctionCV("x", lambda s: s.xyz[0][0])
>>> vol = paths.CVDefinedVolumes(x, 1.0, float("inf"))
>>> ensemble = paths.AllOutXEnsemble(vol)
>>> # okay, we don't have a `trajectory2D` function, but you can imagine it
>>> trajectory = trajectory2D([[[0.0, 0.0]], 
>>>                            [[0.1, float("nan")]]])
>>> ensemble(trajectory)
True

So if you're sampling that ensemble, you might accept that trajectory. Right?

I see. I can add that possibility, too. A problem I discovered is that OpenMM will fail with a ValueError('Particle is nan') when using CPU and run fine using CUDA but return a snapshot with NaNs. And crash once in a while (although I think that is MacOS two graphics card problem related)

Then there is still the question if you should continue until it really breaks or if you detect NaN.

Importantly, this is not a retry.

This is clear. I just remembered John wanting to auto retry. This is different from rejecting and then in the scheme doing a retry.

So, should I remove all the retry stuff? Is this unintended?

So, should I remove all the retry stuff? Is this unintended?

On that I would say ask @jchodera. I don't have a problem with it existing as an option if you've already implemented it, but if you haven't implemented it, it might not be worth implementing.

Almost done. I left the retry as an option which is disabled. So the engine passes a NaNError to the EngineMover that generates a SampleNaNError which will be caught by the SampleMover/MC an return a RejectedMoveChange. I still need to test this somehow and make the capture optional.

Okay, I implemented both. I actually like the idea that a trajectory that contains NaNs will be considered a rejected one. Not sure of the implications on the dynamics, but wouldn't this mean that under the given dynamics (including integrator) trajectories that go 'NaN' are always excluded and we do not want these anyway. Of course, dynamics that could lead to NaN are not really what we want, but for bootstrapping this should work.

I would assume that just running more steps in a PathSimulator and doing the retry will be more or less equivalent in many cases. The only difference is, that if you do not care about the dynamics so much, using retries that do not completely retry the full step, but cut off the NaN trajectory and continue could be faster.

The default is, that the Engine throws a EngineNanError (there is also now a EngineMaxLengthError but that is not the default). The Mover will catch this error and effectively create a RejectedMoveChange that contains the last valid trajectory but will always be rejected, even if it would be a valid sample. I think that cannot happen since the engine should have stopped by then. But to make sure, it will be rejected. This way NaN will not wreck a simulation.

Now that I think of it. Could it be that this will oversample trajectories in the vicinity of crashing ones?

Thanks, JHP!

Some quick comments:

Okay, I implemented both. I actually like the idea that a trajectory that
contains NaNs will be considered a rejected one. Not sure of the
implications on the dynamics, but wouldn't this mean that under the given
dynamics (including integrator) trajectories that go 'NaN' are always
excluded and we do not want these anyway.

Both "retry" and "reject" lead to excluding trajectories with NaNs, right?
"reject" means that we accumulate statistics on those rejections, "retry"
does not?

Of course, dynamics that could lead to NaN are not really what we want,
but for bootstrapping this should work.

I've tried to be very clear: There is no way to completely avoid NaNs.
Using single-precision very much increases the probability of running into
NaNs.

I would assume that just running more steps in a PathSimulator and doing
the retry will be more or less equivalent in many cases. The only
difference is, that if you do not care about the dynamics so much, using
retries that do not completely retry the full step, but cut off the NaN
trajectory and continue could be faster.

The default is, that the Engine throws a EngineNanError (there is also
now a EngineMaxLengthError but that is not the default). The Mover will
catch this error and effectively create a RejectedMoveChange that contains
the last valid trajectory but will always be rejected, even if it would be
a valid sample. I think that cannot happen since the engine should have
stopped by then. But to make sure, it will be rejected. This way NaN will
not wreck a simulation.

Now that I think of it. Could it be that this will oversample trajectories in
the vicinity of crashing ones?

Excellent point. That may lead to altered statistics in the ratio between
different classes of trajectories in the vicinity of classes that produce
NaNs. I guess that means "retry" is preferred?

Question for clarification

The Mover will catch this error and effectively create a RejectedMoveChange that contains the last valid trajectory but will always be rejected, even if it would be a valid sample

I haven't looked at the code yet, but when you say "effectively create a RejectedMoveChange," what exactly do you mean? Mainly, is there a quick and easy way to check for moves rejected because of nan? I would think that the easiest ways to do that would be to either add a reason to rejected move change (which could later include other explanations like "not in ensemble", or "failed acceptance probability test"), or to have a subclass of RejectedMoveChange specific to nan-rejected. But the main thing is that it could be really important to know that one of your replicas/ensembles is getting rejected a lot due to NaN.

NaNs in general

I've tried to be very clear: There is no way to completely avoid NaNs.

@jchodera : We heard you, and while I trust that you're right, I still don't understand why this is so much of a problem that it becomes inevitable that every practical run has a NaN (which seems to be the problem you're having). Can you point to a resource that explains this?

The only time I've encountered dynamics that had this problem was a case where the potential involved taking the square root of a number that would become negative with a large step size. In that case, NaNs were easily avoided by using an adaptive time step integrator. (I guess I had another version of that where even the adaptive time step didn't help, but that was basically because the conserved Hamiltonian had a value on the order of 10^-400, which can't be represented by a double.)

Retry vs. reject

Could it be that this will oversample trajectories in the vicinity of crashing ones?

Excellent point. That may lead to altered statistics in the ratio between
different classes of trajectories in the vicinity of classes that produce
NaNs. I guess that means "retry" is preferred?

This runs counter to my intuition. "Retrying" seems to me like it violates detailed balance. I don't see how retrying a trajectory because it went NaN differs from retrying because a move was rejected, which of course violates detailed balance.

Quick experiment, using a "configurational" sampler instead of a path sampler, but I think the argument is the same:

import numpy as np

def V(x):
    return -np.sqrt(-(x-0.5)*(x+0.5))

# no efforts at subtlety or elegance in this sampling routine
def run_monte_carlo(n_steps, beta=2.0, max_dx=0.05, initial_position=0.0, reject_nan=False):
    position = initial_position
    history = []
    n_nan_trials = 0
    for step in range(n_steps):
        old_weight = np.exp(-beta * V(position))
        new_weight = float("nan")
        while np.isnan(new_weight):
            trial = position + 2.0*max_dx*(np.random.rand() - 0.5)
            new_weight = np.exp(-beta * V(trial))
            if np.isnan(new_weight):
                n_nan_trials += 1
                if reject_nan:
                    new_weight = 0.0
        acceptance_probability = new_weight / old_weight
        if np.random.rand() <= acceptance_probability:
            position = trial
        history.append(position)
    return (history, n_nan_trials)

Run it with:

(reject, n_nan_reject) = run_monte_carlo(1000000, reject_nan=True)
(retry, n_nan_retry) = run_monte_carlo(1000000, reject_nan=False)

A quick check shows us that about 1.43% of trials are NaN in "retry" mode, and 1.38% in "reject" mode. Not really a significant difference (probably due to a 1D nature of this system).

Plotting:

# http://www.wolframalpha.com/input/?i=Integral%5BExp%5B-2.0*Sqrt%5B-(x-0.5)*(x%2B0.5)%5D%5D,+x,+-0.5,+0.5%5D
norm_expV = 0.468451

plt.hist(reject, bins=51, alpha=0.5, normed=True, label='Reject')
plt.hist(retry, bins=51, alpha=0.5, normed=True, label='Retry')
plt.plot(x, np.exp(-2.0*V(x))*norm_expV, label='Exact')
plt.plot(x, np.exp(-2.0*V(x))*norm_expV * 0.96, 'm.', label='Exact*0.96')
plt.plot(x, np.exp(-2.0*V(x))*norm_expV - 0.05, 'm-', label='Exact-0.05')
plt.xlim(-0.6, 0.8)
plt.ylim(0.0, 1.6)
plt.legend()

gives us:

You're correct that reject samples more near the NaN region, but it also samples more correctly. I'm not sure why the exact was so much higher: please feel free to check over my quickly-hacked code to make sure I didn't do anything stupid. But I tried to correct for the difference with ad-hoc shift and ad-hoc scaling of the exact result. In all cases, retry looks low.

@jchodera : We heard you, and while I trust that you're right, I still don't understand why this is so much of a problem that it becomes inevitable that every practical run has a NaN (which seems to be the problem you're having). Can you point to a resource that explains this?

There's no good reference that I know of, though this would be a fun paper to write.

The issue is with single-precision forces and weak thermostats. With double-precision forces, or with strong thermostats, this is much less of a problem. As you get to higher forces (which is often the case with forced or driven barrier crossing), the numerical stability may become poor, and more velocity or random number choices may lead to trajectories that quickly blow up. We see this all the time when using NCMC (in which we drive the potential or coordinates), and have switched to using Metropolized integrators like GHMC that can reject NaNs as a consequence.

The only time I've encountered dynamics that had this problem was a case where the potential involved taking the square root of a number that would become negative with a large step size. In that case, NaNs were easily avoided by using an adaptive time step integrator. (I guess I had another version of that where even the adaptive time step didn't help, but that was basically because the conserved Hamiltonian had a value on the order of 10^-400, which can't be represented by a double.)

Oof, there are a lot of problems with adaptive timestep integrators for molecular simulation. They are generally not symplectic or even measure-preserving, so have a whole host of problems when it comes to anything involving stationary distributions.

Thanks so much for the concrete 1D example! There are a bunch of ways to handle excluded regions in Monte Carlo: The most common is using proposals that omit the forbidden region, in which case we need to then correct for the fact that the probability of ending up in the allowed region is increased, and this correction has to go in the Metropolis-Hastings acceptance criteria in P(old | new) / P(new | old). This explains why the "retry" sampler---which effectively excludes proposals to the forbidden region without correcting the acceptance criteria---gives very erroneous results.

So @dwhswenson is totally right---"retry" will give incorrect statistics, while "reject" should give correct statistics.

I imagine the difference between the "exact" statistics and your histogram would mostly disappear if you changed

plt.plot(x, np.exp(-2.0*V(x))*norm_expV, label='Exact')

to

plt.plot(x, np.exp(-2.0*V(x)) / np.exp(-2.0*V(x)).sum(), label='Exact')

or, even better, numerically integrated over the histogram bins to get the unnormalized probability and then normalized.

Thanks, guys!

I haven't looked at the code yet, but when you say "effectively create a RejectedMoveChange," what exactly do you mean?

Okay, currently (still sorting it out) it will force to return a RejectedPathMoveChange. I added a { 'error': 'nan'} to the details in the Sample as well as the PMC. But i like the idea of using subclasses. I think it is easier to remove these from analysis than having to load all the details to check. Having an additional attribute would also be fine. i just think the first is more elegant. We can still add the concrete reason to the details.

This runs counter to my intuition. "Retrying" seems to me like it violates detailed balance. I don't see how retrying a trajectory because it went NaN differs from retrying because a move was rejected, which of course violates detailed balance.

I did not mean to say that retrying is better. Sorry for the wrong implication. I just thought that both ways will count wrong, but now I think it was actually only retrying that will do that.

So @dwhswenson is totally right---"retry" will give incorrect statistics, while "reject" should give correct statistics.

Yes, I think so. I never thought about the problem of proposing impossible moves, but ii acts as if the energy is infinite. As long as proposals are symmetric you should be fine. And into nan is zero and from nan is also zero, so we're good. And the little test shows it.

Btw. The little analysis is phantastic! Not sure how you patched that so quickly together, but it is pretty cool.

Well, this I do not understand. If I understood David correctly then he said that both will not work correctly but in his example the reject works fine. Is that correct? I would have assumed that both fail but the reject seems not to overcount. How is that possible? Retrying surely fails DB so that is out. I am confused, but I found the reason for the not fitting exact solution. I think you had a sign missing. The 0.46... is the solution, if you flip the sign in the exp. In the other case the integral is 2.24395 and so the normalization is 0.445643.

Okay, so rejecting should be the default because it is the right thing todo. We might keep the retry, becuase in cases like bootstrapping it might be faster. Does that make sense?

Okay, so rejecting should be the default because it is the right thing todo. We might keep the retry, becuase in cases like bootstrapping it might be faster. Does that make sense?

Sounds great to me!

@jchodera : Thanks for the explanation. I'm still surprised that it is such a huge problem, but it helps to at least see what the problem is.

I think you had a sign missing. The 0.46... is the solution, if you flip the sign in the exp. In the other case the integral is 2.24395 and so the normalization is 0.445643.

Nice catch! Yep -- I figured I'd done something stupid. Also, note that the correction of 0.44../0.46.. gives a number pretty close to the scaling by 0.96 that I got by eye, so that part of the plot should be pretty close to what the correct exact value would have been.

If I understood David correctly then he said that both will not work correctly but in his example the reject works fine. Is that correct? I would have assumed that both fail but the reject seems not to overcount.

I had meant that I thought reject works, from a detailed balance standpoint, though I'm a bit disquieted that we can't avoid NaN trajectories. I think this example helps make it clear why it works: as you get closer to the boundary, the slope of the potential goes to infinity. So you can sort-of pretend (not in a mathematically rigorous sense, but just conceptually) that the area with a "NaN" potential is because there's an infinite wall, and the potential is actually infinite outside the domain. In this case, the probably of accepting a move there is, of course, 0. You can try, but you have to reject.

(The other thing I was saying is that, depending which degree of freedom gets the NaN, there could be cases where the NaN-containing trajectory satisfies the ensemble, and is therefore accepted. So we need to have the engine say "this trajectory contains NaNs and can't be trusted" -- which is exactly what you're doing.)

I added a { 'error': 'nan'} to the details in the Sample as well as the PMC. But i like the idea of using subclasses. I think it is easier to remove these from analysis than having to load all the details to check.

I think either approach is fine. To me, the trade-off is:

• subclass: easier to check
• attribute: easier to generalize

The only thing is I would suggest using rejection_reason, or something like that, for the attribute name, instead of error. My thought is that many (hopefully most) rejection reasons would not be errors!

Question on this: Is there any possibility that the NaNs are being hidden somehow when running in an ipynb? I've seen this error occur in a pure-python script on alanine dipeptide TPS, but I haven't gotten it to occur in an ipynb yet.

Question on this: Is there any possibility that the NaNs are being hidden
somehow when running in an ipynb? I've seen this error occur in a
pure-python script on alanine dipeptide TPS, but I haven't gotten it to
occur in an ipynb yet.

Are you running on the same simtk.openmm.Platform within the ipynb as when
running from the command line?

Are you running on the same simtk.openmm.Platform within the ipynb as when
running from the command line?

Should be -- I'm reloading the same OPS engine from a storage file.

(I would check, but the notebook version is currently at about step 5000 of 10000)

I'm fully in support of changing that default!

Good. Will do that. And then it is ready to merge. Once it passes tests.

While I am at it. To unify the options of an engine, would it be okay to have the integrator as a regular attribute, like openmm has? We have the topology this way and it seems more consistent that way.

I think I would make a separate PR, because this affects all ToyEngines, everywhere. So not now...

While I am at it. To unify the options of an engine, would it be okay to have the integrator as a regular attribute, like openmm has? We have the topology this way and it seems more consistent that way.

DynamicsEngine shouldn't have it (and shouldn't have topology as an attribute, either). I'm not sure why ToyEngine still has a topology. When I originally did that, it was just to shoehorn it into something that looked like what you'd done for OpenMM -- several parts of the code required a topology attribute to get, e.g., n_atoms. ToyEngine should have a potential_energy_surface (or just potential) object, which should contain the information equivalent to a topology. I think the major improvements you made about the engine, snapshot features, etc., should make that really easy clean up.

I think the way it should work is this: DynamicsEngine has a minimal set of initialization parameters (which are then attributes). Each subclass should have a different set that makes sense for that subclass. integrator makes sense for both ToyEngine and OpenMMEngine, but maybe there would be some other engine where it was redundant. topology makes sense for OpenMMEngine, but does not make sense for ToyEngine (which should have a potential_energy_surface instead). It's okay to add extra parameters to the initialization of subclasses, since we don't expect them to be perfectly interchangeable.

See also #65.

Good. Will do that. And then it is ready to merge. Once it passes tests.

... you haven't pushed yet... waiting eagerly! :-) (Still at the office; would love to run the AD example overnight)

I think the way it should work is this: DynamicsEngine has a minimal set of initialization parameters (which are then attributes). Each subclass should have a different set that makes sense for that subclass. integrator makes sense for both ToyEngine and OpenMMEngine, but maybe there would be some other engine where it was redundant. topology makes sense for OpenMMEngine, but does not make sense for ToyEngine (which should have a potential_energy_surface instead). It's okay to add extra parameters to the initialization of subclasses, since we don't expect them to be perfectly interchangeable.

Alright. What I thought and wanted. So we clean this up shortly after this. DynamicsEngine now only has stuff about the number of steps, maximum frames and all the error handling which makes sense for all of them. These are options. All the really specific stuff are parameters on the constructor.

Will need a few more minutes. Need to change to named tuple .

Alright. What I thought and wanted. So we clean this up shortly after this.

Excellent. One question: will this mean I have to re-run the AD examples again? It sounds like this will break backward compatibility. (That's fine, and I absolute must start these runs ASAP so I can have the long fixed TPS run done in time for my talk next week. Just checking the situation.)

Will need a few more minutes.

No worries. I'm working on visualizations of how we use ensemble to do a lifetime calculation (rate for long MD) for that talk. Plus, the custom strategy is running its last version, and has 12000 seconds left. ;-)

Actually. The namedtuple stuff is kind of weird. Let's postpone that to the engine cleanup. It is kind of messy changing all the options['xyz'] to options.xyz

The namedtuple stuff is kind of weird. Let's postpone that to the engine cleanup.

Sure... it's also already immutable-enough, if we make the habit of accessing it as engine.foo instead of engine.options['foo']. (Because engine.foo = bar doesn't work, right?)

(Although I like the extra certainty in setup from a namedtuple; removed the need for some of the "check allowed options" code you wrote)

nosetests work. Now the real ones...

(Because engine.foo = bar doesn't work, right?)

Yes, it does not.

That's it. Ready to be merged.

Later engine cleanup will be separate PR.

Great! I've been reviewing recent changes as you made them, so I'll merge this as soon as it passes tests. Then I'll start giving it some thorough stress-tests on the AD system starting tonight, so pretty soon we should be certain that it fixed the NaN problem! @jchodera : once this is merged, it would really help if you can check that this solved your NaN problems, too!