Fix broken WHAMing if an InterfaceSet has more than 10 interfaces - #967
Conversation
Codecov Report
@@ Coverage Diff @@
## master #967 +/- ##
=======================================
Coverage 80.25% 80.26%
=======================================
Files 136 136
Lines 14449 14449
=======================================
+ Hits 11596 11597 +1
+ Misses 2853 2852 -1
Continue to review full report at Codecov.
|
dwhswenson
left a comment
There was a problem hiding this comment.
Huh... no idea how that happened, but that code dates back to 2015. My guess is that there was some weirdness with pandas that forced me to name things as strings, and then that fact got forgotten in the new TIS analysis (written in 2017, IIRC).
|
@dwhswenson Quick question; nothing got published with OPS and more than 10 interfaces in an interface set, right? As this wrong ordering can lead to wrong answers (for me it was 7 orders of magnitude) if (by sheer luck) there is data from the inner interface that is valid for outer interfaces. |
|
@arjunwadhawan You might want to check on this. Note that this should only be relevant for the new TIS analysis. This will only happen if you set the cutoff very small (or have interfaces that overlap way more than is reasonable). If you use the rule of thumb that successive interfaces should have 20% overlap, then fewer than 1 in a million trajectories from interface 1 should overlap with interface 10, and that would be ignored by cutoff. |
|
@sroet : With the error, how did plots of the total crossing probability look? I would hope that this problem would have been evident in the TCP, compared to the per-interface crossing probability. |
The reason why it happened for me was because I started sampling with an A->B path in all interfaces, assuming it wouldn't matter for convergence if I ran 100_000 MC steps (not unreasonable for real systems), followed by the error described in #920 which was why I set the cutoff to (essentially) 0 (which would could have prevented the "wrong" result. @dwhswenson so for me it looked reasonable as my toy system had a reasonable cliff (10 KbT triangle barrier between -1 and 1, with stable states at -2 and 2.5), but the TCP should have (ofcourse) been -10 instead of the -3 found.
And from the tcp this was not entirely clear (to me) why is was 7 orders of magnitude off (the BA one is should be mirrored for this plot). Looking at it with what I know now, you indeed see the giant cliff that is way earlier than you would expect and then recovers again.
The current plots looks like this (1/10 th of newly sampled data (that came from a run that decorrelated properly before running production)):
|
|
I'm having trouble interpreting the plots from the old data. Is the first one the total crossing prob? The standard thing is to plot TCP on the same plot as the ensemble crossing probabilities, as you do at the end for the current data. TCP should look exactly like the innermost crossing probability until noise in the innermost crossing prob starts to matter. It's hard to tell, but it doesn't look like that would be the case if you put the first plot on the first per-ensemble plot. |
Yeah that is the TCP, and I started getting better at plotting these graphs. It is indeed better do plot them together with the ensemble crossings. Looking at them side by side (on my monitor) it is not the case that the total tcp follow the innermost probability, it drops more rapid than the innermost probability |
Okay, I'll hope that this has therefore been noticed by anyone who has tried it before. At least @arjunwadhawan and anyone else in the @bolhuis group would have made such plots, and noticed if something was way off. But this is a critical bug. Cutting 1.4.2 release ASAP. |
This fixes an error that took a long time to track down and manifests itself with the same error as #920 (or with a wrong answer in more extreme cases).
TLDR: WHAMing with more than 10 interfaces in an interfaceset was horribly broken, should now work
14bcc2f fixes the issue and adds a regression test
bf8241c updates the test file to use
pytestinstead (incremental update for #756)How this PR came to be:
This manifested itself as the same error as #920 (partial stacktrace):
From that issue we know that that happens if everything before the interface is set to
0and everything after that interface is also0. So the first thing I did was setting cutoff to0to circumvent the same issue as in #920.This did not solve the error so I decided to add print functions to the
prep_reverse_cumulativefunction ofnumerics.wham.py(sidenote:%autoreloadis great and makes it so that I can debug this without waiting an hour to reload my data).This is is the first block of print(
df), before any cleaning, notice anything peculiar about the column ordering?Now here is where things go horribly wrong: line
137-142ofwham.pyare:Because this is how
self.interfaceslooks before converting it to a series:and this is after
So this now maps the 3rd interface to the data of the 10th, but even worse is the 13th interface is mapped to the data of the second! The data from my 2nd interface does not have any data that is past the 13th, so this results in a column of
0s after cleaning (everything up to the 13th interface is set to 0).Now, this PR does not actually alter
wham.pybecause the real issue is the sorting of the columns of the dataframe.This is done at line
124-127ofopenpathsampling/analysis/tis/crossing_probability.pySo this sorts axis=1, and what we see is that th columns are sorted, but as if they are
strs, as those are sorted character by character, which is how you end up with0, 1, 10, 11, ..., 19, 2, 20, 3, .... Now where do these strings come from?That is what is actually fixed in this PR
If we look at the called function (
histograms_to_pandas_dataframe) on line 432-433 the follwing code exists:and these
strindices end up in the dataframe. Hence this PR to convert them to anintinstead.