Pymbar4 bootstrapping

[mrshirts]

Adding support for bootstrapped errors to all expectations.

Computes bootstrap uncertainties and covariances for:

• compute_free_energy_differences
• compute_expectations
• compute_multiple_expectations
• compute_perturbed_free_energies
• compute_entropy_and_enthalpy

These are accessed by using the uncertainty_method = "bootstrap". Note that MBAR needs to be initialized with "n_bootstrap=(some integer)" samples in order to do use this method - an error is omitted otherwise. Possibly could be done differently, but the logic is simpler this way for now.

I could use a little help in rewriting the pytests. Currently, to get bootstrap uncertainties from various functions, as one will now want to check both the default uncertainties and the bootstrap uncertainties, and MBAR needs to be initialized differently each way, which the pytests don't seem to handle.

Note this is built off of #444, so many of the differences are from that branch.

[codecov]

Codecov Report

Merging #448 (ed4749d) into pymbar4 (ee92147) will decrease coverage by 0.46%.
The diff coverage is 73.79%.

[zhang-ivy]

@mrshirts : I should be able to test this out in the next couple days.

What's your expected timeline for getting this merged and pymbar 4 released? The reason I ask is because I am generating data for a paper and ideally I would be able to use a stable version of pymbar to do the bootstrapping.

@mikemhenry is planning on fixing the pymbar 3 unit tests (which will take about a day) so that we can get #302 merged in and then put out a minor release, but if you're planning on having pymbar 4 released very soon (this week), then maybe its not worth his effort.

[mrshirts]

I would like to do pymbar4 in the next ~10 days or so - depends mostly on other stakeholders signing on. I would not like to extend pymbar3 much longer. I have already identified a problem in pymbar3 bootstrapping in #302 (it bootstraps among ALL samples, rather than bootstrapping among samples at each K), so I don't know how much to trust the #302 version exactly.

Additional testing of this branch by people would be useful and ensure that it gets out soon!

[zhang-ivy]

@mrshirts : I had a discussion today with @jchodera and the perses / openmmtools developers -- since our tools depend on pymbar3, we are willing to take over maintaining pymbar3, so I think @mikemhenry will work on cutting a bugfix release that fixes the CI and contains the bootstrap branch you wrote awhile ago (#302).

I'll make sure to carefully compare the code in this branch with that of the older branch and will fix any bugs in the old branch. I will also test both branches to make sure they yield similar results. Will post an update once I've done this!

[mrshirts]

I have concerns about maintaining multiple branches. I think that it's fine to keep pymbar3 alive for a while to make sure that toolchains have time to adapt, and port any bugfixes back, but I think it's better to keep new development linear. The output formats of pymbar3 are really restrictive in addition.

Note that #302 is incorrect in at least one way (randomly selects from all samples, instead of ones from one point), and it doesn't support bootstrapping of expectations, etc.

@zhang-ivy, shoot me an email at Colorado address and we can set up a time to talk over things with @jchodera / perses developers? I'd rather allocate time in helping people to move forward rather than supporting old, out-of-date code.

[mikemhenry]

@zhang-ivy are you going to make a new PR with just the bootstrapping stuff, or take overt this one? I just want to make sure either on this PR or a future PR, the correct branch pymbar3 is selected.

[zhang-ivy]

@mikemhenry : This branch is for pymbar4, the one that we want to merge into the pymbar3 long term maintenance branch is in #302

[mrshirts]

Yes, so what we want is once #302 is merged in to master, and a 3.0.8 version is released, then afterwards master should point to pymbar4, and what is currently master should be pymbar3_ltm or something like that. @zhang-ivy , @jchodera , is that your understanding?

[mrshirts]

The failing histogram test just seems to be that a histogram with random samples is not getting data in that histogram. I will check over that problem again - maybe the way it's handled needs to be tweaked so this is less likely to happen.

[zhang-ivy]

Yes, so what we want is once #302 is merged in to master, and a 3.0.8 version is released, then afterwards master should point to pymbar4, and what is currently master should be pymbar3_ltm or something like that. @zhang-ivy , @jchodera , is that your understanding?

Yes!

[mrshirts]

@Lnaden, thoughts about adjusting the pymbar tests relating to uncertainties? Specifically, what would you look at for the fact that we can't just add 'bootstrap' as an uncertainty_method variable, since MBAR needs to be initialized with bootstraps for this to happen? So therefore, when uncertainty_method = bootstrap is called, we need to those MBAR objects to be initialized with bootstrapping - but we don't want ALL of them to be initialized with bootstraps, because that takes a lot of extra computational time.

[mrshirts]

I think the function is supposed to be inside of the loop, as the function (if I recall!) need to be redefined with the loop. I remember struggling with the right way to do this. Now, there's probably a better way to do it than I did, but I will need to go back and make sure the behavior hasn't changed. Before release at least.

[mrshirts]

Ah, I see what you did. OK, I think that will work, I will double check to make sure. It probably improves performance . . .

[Lnaden]

thoughts about adjusting the pymbar tests relating to uncertainties?

A couple of options. The simplest is to do nothing special and just initialize the n_bootstraps on the mbar_and_test fixture. Because the mbar_and_test fixture is a scope=module fixture, it only ever gets initialized once for the whole module of tests per parameter, a total of twice. At n_bootstrap=1000, this adds about 15 seconds to the total test suite time on my local machine.

There is a separate option where we can rebuild the tests to pass bootstrap parameters, but I think that adds far more complexity than is needed here.

Let me do a quick tinkering to test bootstrap uncertainty method and to add bad initialization options to make sure things fail correctly as well.

[mrshirts]

it only ever gets initialized once for the whole module of tests per parameter, a total of twice.

Ah, I didn't realize that!

And 1000 bootstraps are not needed. 200 should be plenty for decent uncertainties (so 3 seconds?) - maybe even 50, but then you might get a few statistically off samples.

[mrshirts]

Note that when running pymbar4 and calling compute_free_energy_differences() I see the following warning:
Failed to reach a solution to within tolerance with hybr: trying next method
Not sure if this is a problem, but just wanted to point out this behavior in case its not what you intended it to be.

Is it returning a successful result in the end? Then it is working as intended (it tried 'hybr', then something else). Maybe there should be an additional message saying what method succeeded?

[mrshirts]

Any idea why this is happening? I would've thought it'd be resolved by #443 but then I realized this error occurs at line 730 in mbar.py, and the PR fixes a later line.

The fix was in #441, #443 was a bugfix to #441, so maybe it's still #441?

[mrshirts]

One other question that would be great to know on a hard system - how does the pymbar4 bootstrapping work for entropy and enthalply errors? (n_bootstraps in initializations and uncertainty_method = "bootstrap" in compute_entropy_and_enthalpy).

[mrshirts]

Thanks for checking these!

[zhang-ivy]

Is it returning a successful result in the end? Then it is working as intended (it tried 'hybr', then something else). Maybe there should be an additional message saying what method succeeded?

Yes, its returning a successful result. I think adding an additional message with the method that succeeded would be nice!

The fix was in #441, #443 was a bugfix to #441, so maybe it's still #441?

Just saw this comment: https://github.com/choderalab/pymbar/pull/441/files#diff-39154a13f665187718c90bf1b7a14a9edb282028dd0fbd94a8dfb02168e42580R648-R652
In case the link doesn't work, its this line:

logfactor = 0  # make sure all results are larger than this number.
                       # We tried 1 before, but expecations that are all very small (like
                       # fraction folded when it is low) cannot be computed accurately. 
                       # it's possible  that something really small but > 0 might avoid
                       # errors, but no errors have occured yet. 

Do you think we should change logfactor to small number like 0.001?

One other question that would be great to know on a hard system - how does the pymbar4 bootstrapping work for entropy and enthalply errors? (n_bootstraps in initializations and uncertainty_method = "bootstrap" in compute_entropy_and_enthalpy).

Sure I can check this later tonight!

[mrshirts]

Do you think we should change logfactor to small number like 0.001?

In the current code, it's something like 4x(the smallest mantissa accessible to numpy, like 1e-15). That has worked so far. 0.001 is not so good, since then numbers near in magnitude 0.001 can't be averaged.

[mrshirts]

Sure I can check this later tonight!

Great!

[zhang-ivy]

@mrshirts : I just tested bootstrapping for compute_free_energy_differences() and compute_entropy_and_enthalpy() with this code (with pymbar4 only, since pymbar3 doesn't do bootstrapping for compute_entropy_and_enthalpy()):

    # Test compute_free_energy_differences()
    decorrelated_u_ln, decorrelated_N_l = analyzer._compute_mbar_decorrelated_energies()
    mbar = MBAR(decorrelated_u_ln, decorrelated_N_l)
    results = mbar.compute_free_energy_differences(compute_uncertainty=True)
    print("compute_free_energy_differences()", results['Delta_f'][0, -1], results['dDelta_f'][0, -1])
    
    # Test compute_entropy_and_enthalpy()
    results = mbar.compute_entropy_and_enthalpy()
    print("compute_entropy_and_enthalpy()", results['Delta_f'][0, -1], results['dDelta_f'][0, -1], results['Delta_u'][0, -1], results['dDelta_u'][0, -1], results['Delta_s'][0, -1], results['dDelta_s'][0, -1])
    
    # Test compute_free_energy_differences() with bootstrapping
    mbar = MBAR(decorrelated_u_ln, decorrelated_N_l, n_bootstraps=200)
    results = mbar.compute_free_energy_differences(compute_uncertainty=True, uncertainty_method='bootstrap')
    print("compute_free_energy_differences() bootstrapped", results['Delta_f'][0, -1], results['dDelta_f'][0, -1])
    
    # Test compute_entropy_and_enthalpy() with bootstrapping
    results = mbar.compute_entropy_and_enthalpy(uncertainty_method='bootstrap')
    print("compute_entropy_and_enthalpy() bootstrapped", results['Delta_f'][0, -1], results['dDelta_f'][0, -1], results['Delta_u'][0, -1], results['dDelta_u'][0, -1], results['Delta_s'][0, -1], results['dDelta_s'][0, -1])
    

Here are the results.

For case 1 (barnase:barstar mutation):

compute_free_energy_differences() -52.62531138448979 0.04448291789808753
compute_entropy_and_enthalpy() -52.62531138449739 0.044482917898087226 -58.26693039969541 9.695425671134956 -5.64161901519401 9.693498249036294
Failed to reach a solution to within tolerance with hybr: trying next method
compute_free_energy_differences() bootstrapped -52.62531138448979 0.10343092206869589
compute_entropy_and_enthalpy() bootstrapped -52.62531138449739 0.10343092206869589 -58.26693039969541 69.68171363279014 -5.64161901519401 69.67825680667228

For case 2 (ala dipeptide mutation):

compute_free_energy_differences() -52.27795853969293 0.056782925355530894
compute_entropy_and_enthalpy() -52.277958539692946 0.05678292535553127 -57.77404736274184 1.8697858250872512 -5.496088823048922 1.8701373837035669
Failed to reach a solution to within tolerance with hybr: trying next method
compute_free_energy_differences() bootstrapped -52.27795853969293 0.16964559913815944
compute_entropy_and_enthalpy() bootstrapped -52.277958539692946 0.16964559913815944 -57.77404736274184 26.242413812942953 -5.496088823048922 26.239604799661123

Bootstrapping for free energy differences looks good (as I already found in my tests last week), but the bootstrapped uncertainty for the entropy and enthalpy looks much larger than expected...

[zhang-ivy]

In the current code, it's something like 4x(the smallest mantissa accessible to numpy, like 1e-15). That has worked so far. 0.001 is not so good, since then numbers near in magnitude 0.001 can't be averaged.

I just tried copying this line from pymbar4 mbar.py:

logfactor = 4.0 * np.finfo(np.float64).eps
        # make sure all results are larger than this number.
        # We tried 1 before, but expecations that are all very small (like
        # fraction folded when it is low) cannot be computed accurately.
        # 0 causes warnings in the test with divide by zero, as does 1*eps (though fewer),
        # and even occasionally 2*eps, so we chooose 4*eps

to replace logfactor here.

But I'm still seeing the same warning with pymbar3:

/home/zhangi/miniconda3/envs/perses-paper/lib/python3.9/site-packages/pymbar-3.0.5-py3.9.egg/pymbar/mbar.py:777: RuntimeWarning: divide by zero encountered in log
  Log_W_nk[:, sa] = np.log(A_n[i, :]) + Log_W_nk[:, l]

I'm not seeing this error when running pymbar4. I can take a closer look at the diff for computeExpectationsInner next week (will be on vacation the rest of the week), but if you have any pointers for what change in pymbar4 will be necessary to port over to pymbar3, that would be super helpful!

[mrshirts]

but the bootstrapped uncertainty for the entropy and enthalpy looks much larger than expected...

Interesting, can you post the data set somewhere I can look at it? (like google drive?) It was looking OK for my tests, but maybe there's an error that didn't show up with that relatively easy test.

[mrshirts]

But I'm still seeing the same warning with pymbar3:

Yeah, I'd have to take a look at the data set to see.

[zhang-ivy]

The data for these experiments isn't too big, so i'm just going to attach them here. I included the u_kn and n_k for barnase:barstar mutation as well as for ala dipeptide mutation.
pymbar3_divide_by_zero.zip

[zhang-ivy]

Note that barnase:barstar and ala dipeptide have 24 and 12 thermodynamic states, respectively, but I've added unsampled endstates for lambda = 0 and lambda = 1 such that the states are ordered like [unsampled endstate for lambda 0, sampled endstate for lambda 0, ... sampled endstate for lambda = 1, unsampled endstate for lambda = 1], so the matrices have K = 26 and 14, respectively. Not sure this matters, but just noted it in case it does.

[mrshirts]

Hmm. I can't seem to load them. When I unzip and run I get:

with open("u_kn.npy", "rb") as f:
u_kn = np.load(f)['arr_0']

I get: ValueError: Cannot load file containing pickled data when allow_pickle=False

And when I try
u_kn = np.load(f,allow_pickle=True)['arr_0']

I get:
*** _pickle.UnpicklingError: Failed to interpret file <_io.BufferedReader name='u_kn.npy'> as a pickle

Any thoughts?

The size of u_kn.npy is 10932608.

[zhang-ivy]

Hmm, this is working for me:

with open("bnbs/u_kn.npy", "rb") as f:
    u_kn = np.load(f, allow_pickle=True)

with open("bnbs/n_k.npy", "rb") as f:
    n_k = np.load(f, allow_pickle=True)

[mrshirts]

Hah, figured out what was doing wrong. Back to work!

[mrshirts]

Ah, I figured out what went wrong with the expectations bootstrapping. A couple of arrays that needed to be bootstrapped were not. It appears to give very close results now.

[zhang-ivy]

@mrshirts : Nice!

Did you have a chance to look into why this warning is appearing for pymbar3 and not pymbar4?

/home/zhangi/miniconda3/envs/perses-paper/lib/python3.9/site-packages/pymbar-3.0.5-py3.9.egg/pymbar/mbar.py:777: RuntimeWarning: divide by zero encountered in log
  Log_W_nk[:, sa] = np.log(A_n[i, :]) + Log_W_nk[:, l]

If not, no worries, I will look into that this week.

[mrshirts]

I'll try to take a look later today.

[mrshirts]

There's one really weird point - 24090 is much lower than the other points, and is giving problems. Not sure if that affects things.

[mrshirts]

So it fails because with the logic, there always be a number that is zero. So the idea is when subtracting the minimum value for each row, you actually subtract (minimum_value - small amount), which will result of all entries in that line to be greater than zero. So setting a slightly larger value for logfactor should solve the problem.

[zhang-ivy]

There's one really weird point - 24090 is much lower than the other points, and is giving problems. Not sure if that affects things.

Can you clarify what you're referring to when you say 24090? And what specifically is lower? the free energy?

[zhang-ivy]

So it fails because with the logic, there always be a number that is zero. So the idea is when subtracting the minimum value for each row, you actually subtract (minimum_value - small amount), which will result of all entries in that line to be greater than zero. So setting a slightly larger value for logfactor should solve the problem.

Hmm I set logfactor to logfactor = 4.0 * np.finfo(np.float64).eps in my local version of pymbar3's bootstrap branch and I'm still seeing the same warning. Same thing happens if I bump it up a bit and use logfactor = 10.0 * np.finfo(np.float64).eps

Update: if I set logfactor to 0.001, the warning goes away! But earlier in this thread, you said "0.001 is not so good, since then numbers near in magnitude 0.001 can't be averaged."

Also, I should note that this warning only shows up when I call computeEntropyAndEnthalpy(), which is not actually something I'm interested in computing (I usually only compute getFreeEnergyDifferences()), so its not imperative for my work to fix this warning, but for long term maintenance purposes, its probably something Mike Henry and I should fix.

[mrshirts]

Ah, OK, I understand better now, and why it was solved in pymbar4. When we subtract out the minimum, in order for it not to be zero, we need to subtract np.min(A) + some factor such that (A-np.min(A))/np.min(A) is greater than eps. This logic is fixed in pymbar 4, but not pymbar 3 - i.e. the factor we subtract is not too small in absolute terms, but not in relative terms; so technically A-(A_min-eps) = eps for the lowest value, but that is not enough, since if A is too big, then A-(A_min-eps) will round to 0. See improved logic in pymbar 4 on

pymbar/pymbar/mbar.py

Line 872 in e8a68ff

logfactors = np.zeros(len(A_list))

and the next few lines.

I don't fully understand why the code was not working for very small expectations (expectations of occupation probabilities that were much lower than 1), but it clearly was not, and the new changes seem to fix this. I will continue to test with the new code to see if there are any issues we didn't get to.

[mrshirts]

Note that the particular data set is not great, because there's one data point that is 1400 units below the rest of the data points, leading to highly likely overflow. So I guess it's a hard case! But at least you see why the pymbar 4 one does work, and can decide how much to port over.

[zhang-ivy]

Thanks for the explanation! I ported the pymbar4 changes you mentioned to my local branch of pymbar3 and the divide by 0 warning goes away.

But at least you see why the pymbar 4 one does work, and can decide how much to port over.

It seems like its a pretty easy fix to port the changes over (only involves changing a few lines), but I'm curious how bad is it to see the divide by 0 warning? In other words, are the returned expectations wrong when we see the warning? I would assume they aren't, but i'm not sure what the implications of dividing by 0 are on the resulting expectations.

[mrshirts]

but I'm curious how bad is it to see the divide by 0 warning? In other words, are the returned expectations wrong when we see the warning? I would assume they aren't, but i'm not sure what the implications of dividing by 0 are on the resulting expectations.

My understanding is that give the error, the point that is zero (which is the minimum) is left out, which may or many not affect the results. If the lowest point is very non-representative it would affect it a lot, if it was relatively representative then it would not affect things much.

[zhang-ivy]

Ah I see. It sounds like it would definitely be good to port the changes to pymbar3 then. I'll create a new branch for that fix.