bcom vs coverap: discrepancy #83
-
|
Hi all, We've been using the bcom tool for a while on a few systems and wanted to do a bit of a sanity check by making some comparisons to what the coverlap tool provides. Based on my understanding of these two papers (DOI: 10.1103/PhysRevE.65.031910, DOI:10.1021/ct2002754) if we compute the covariance overlap (using coverlap) of the 1st half of the trajectory with the full trajectory and average this with the covariance overlap of the 2nd half of the trajectory with the full trajectory, then this value should (in principle) be the same as if (using bcom) we have a block size of NFRAMES/2 (half the trajectory length) as this results in 2 blocks whose overlap is being compared with the full trajectory and then averaged together. Assuming this, we ran bcom using the following command: We then extract the line where there are 2 blocks each 1/2 of the trajectory (in this case 14950 frames as we're skipping the first 50). For a sample protein this value is 0.896. Then to test coverlap we first generate the PCA matrices using svd for 1st 1/2, 2nd 1/2, and full: Then run coverlap comparing both halves to the whole: The average of these two covariance overlap values for a sample protein is 0.828. My question is why is there this discrepancy? Is there something in the commands I'm missing? I looked through the code a bit and couldn't discern any flags/options that are being set for one that aren't being set for the other, and they both make use of the same "covarianceoverlap" command. So I assume this means that the input eigenvectors and eigenvalues are different, but I cannot figure out how or why. Or is there something inherently different in the processes that I didn't see? Cheers, Kari Gaalswyk |
Beta Was this translation helpful? Give feedback.
Replies: 1 comment 5 replies
-
|
The only things that occur to me off the top of my head are round off error (unlikely to be that large, but it’s a lot of frames and if I recall correctly we do the calculation with floats) or something to do with the alignment. Of the two, I think the latter is more likely. I think the blocked pca uses the same alignment from the full trajectory alignment, rather than redoing it. In that case, the two half alignments are slightly different. You could test this by pre-aligning the trajectory with the @tromo , do you have any thoughts? |
Beta Was this translation helpful? Give feedback.
-
|
I think you nailed it with number 2 above... |
Beta Was this translation helpful? Give feedback.
-
|
Ok. So I used the aligner tool to align the full trajectory and then ran svd and coverlap as shown below:
And while this reduces the difference, it doesn't eliminate it entirely. I ran this for 7 different trajectories (replicates of the same system) and the difference ranges from 0.0005 to 0.045 with neither always being greater than the other. For the system I mentioned earlier we now have a coverlap of 0.849 and a bcom of 0.896. I'd assume that if there were rounding issues that the difference would be a bit more regular, but perhaps not. Best, |
Beta Was this translation helpful? Give feedback.
-
|
My intuition is that what's left is numerical, but I'm not even sure how to test. |
Beta Was this translation helpful? Give feedback.
-
|
I agree. In addition, there are under-the-hood differences that will contribute to discrepancies. The svd tool uses dgesvd() to calculate the SVD, whereas bcom calls ssyev() to get an eigen decomp (and I believe the default is to normalize the eigenvalues by block size). I would only expect to get "reasonably" close numbers, not identical ones (nor a systematic difference). |
Beta Was this translation helpful? Give feedback.
-
|
for what it's worth, those differences are also well below the threshold of interpretability. |
Beta Was this translation helpful? Give feedback.
The only things that occur to me off the top of my head are round off error (unlikely to be that large, but it’s a lot of frames and if I recall correctly we do the calculation with floats) or something to do with the alignment.
Of the two, I think the latter is more likely. I think the blocked pca uses the same alignment from the full trajectory alignment, rather than redoing it. In that case, the two half alignments are slightly different. You could test this by pre-aligning the trajectory with the
alignertool, then tellingsvdnot to align. If my hypothesis is correct, that should make the discrepancy go away. In that case, it’s not a bug.@tromo , do you have any thoughts?