search

dimchris / mdanalysis

Automatically exported from code.google.com/p/mdanalysis
0 stars 0 forks source link

while calculating rmsd, it returns nan in some places where is it supposed to be 0 #111

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
The following code produces the bug:

import MDAnalysis
from MDAnalysis.tests.datafiles import PRMncdf, NCDF
import MDAnalysis.analysis.align
u = MDAnalysis.Universe(PRMncdf, NCDF)
ca = u.selectAtoms('name CA')
x = ca.get_positions(u.trajectory[2])
rmsd = MDAnalysis.analysis.align.rmsd(x, x)
print 'rmsd(t[2], t[2]) = ', rmsd
assert rmsd == 0

Produces 'nan' instead of 0.

(Doing this for all atoms gives 1.9621127973121917e-06, which is close enough 
to 0.)

Original issue reported on code.google.com by bala.bio...@gmail.com on 4 Jun 2012 at 7:10

GoogleCodeExporter commented 9 years ago 
I'm having trouble finding the files associated with PRMncdf and NCDF, and 
didn't see them in the git data repo. Is there another branch, or am I just 
missing something? Also, can you confirm that this problem arises with other 
test data besides these two files?

Original comment by joshua.a...@gmail.com on 9 Jun 2012 at 6:06

GoogleCodeExporter commented 9 years ago 
THese files are in the DevelopmentBranch and after installing the 
MDAnalysisTests package you can access them in the way describe above. Or find 
them in the MDAnalysisTests/data directory. 

I have only seen this nan in this specific test case.

Original comment by orbeckst on 9 Jun 2012 at 7:43

GoogleCodeExporter commented 9 years ago 
Hi,
I have attached a sample trajectory and top file. I used the following code and 
i found the value nan instead zero at few places.

import MDAnalysis
from MDAnalysis.analysis.align import rmsd
import numpy as np

traj=MDAnalysis.Universe('TEST.prmtop','TEST1.NCDF')

T1=traj.selectAtoms('name CA')
NO_FR=len(traj.trajectory)
RMSD=np.zeros((NO_FR,NO_FR),dtype=float)

for i in range(NO_FR):
    for j in range(NO_FR):
        cor1=T1.coordinates(traj.trajectory[i])
        cor2=T1.coordinates(traj.trajectory[j])
        RMSD[i][j]=rmsd(cor1,cor2)
print RMSD

Original comment by bala.bio...@gmail.com on 11 Jun 2012 at 6:13

Attachments:

GoogleCodeExporter commented 9 years ago 
Hi Bala,

Just so I can do a quick test external of MDAnalysis, can you post the indices 
of the CA atoms either as a text file or a .npy file (zeros based)? Since I can 
read the netcdf file outside of MDAnalysis, I want to see if I can reproduce 
the bug with my raw pyqcprot code. 

Also just as a numpy performance tip, `RMSD[i,j]` is about 3x faster than 
`RMSD[i][j]` for an array of that size, since the latter amounts to first 
grabbing the entire row and then the column, instead of going to the data 
element directly.

Original comment by joshua.a...@gmail.com on 11 Jun 2012 at 12:27

GoogleCodeExporter commented 9 years ago 
I cannot replicate this error on my machine:

In [12]: import MDAnalysis
In [13]: from MDAnalysis.tests.datafiles import PRMncdf, NCDF
In [14]: import MDAnalysis.analysis.align
In [15]: u = MDAnalysis.Universe(PRMncdf, NCDF)
In [16]: ca = u.selectAtoms('name CA')
In [17]: x = ca.get_positions(u.trajectory[2])
In [18]: rmsd = MDAnalysis.analysis.align.rmsd(x, x)
In [19]: print 'rmsd(t[2], t[2]) = ', rmsd
rmsd(t[2], t[2]) =  3.06563401572e-06

I am however getting nan values in the first two diagonal entries of your 
second example. I can also confirm that using the bare coordinates with my 
qcpprot code (outside of the context of MDAnalysis) produces the same error. 
I'm looking into this now and will get back to you on this soon.

Josh

Original comment by joshua.a...@gmail.com on 16 Jun 2012 at 5:38

GoogleCodeExporter commented 9 years ago 
Ok, figured out the problem. It was fixed in the most recent version of Doug 
Theobald's qcp code, and had to do with small floating point erros that 
sometimes resulted in negative numbers where there shouldn't have been any. 

I'm attaching a fixed version of pyqcprot.pyx since I'm having trouble wrapping 
my head around git and forcing the build system to regenerate the c code right 
now. Hopefully one of the developers can patch the code in and cythonize it for 
me.

Original comment by joshua.a...@gmail.com on 16 Jun 2012 at 8:45

Attachments:

GoogleCodeExporter commented 9 years ago 
This issue was closed by revision 435430e64caf.

Original comment by orbeckst on 17 Jun 2012 at 1:14

GoogleCodeExporter commented 9 years ago 
Thanks Josh. I updated the code and regenerated the c code.

On my machine, I could reproduce the NAN with the old code whereas the new one 
gives a small number. Because of this machine dependency I did not add a 
specific test case.

Original comment by orbeckst on 17 Jun 2012 at 1:20