rsmd atom_indices checks fix

[PolyachenkoYA]

As I understand, atom_indices are indices of atoms in the target, and ref_atom_indices are indices in the reference.

However here atom_indices are compared to the reference size, but they shouldn't be since they index atoms in a different trajectory (topology).

The same applies here where ref_atom_indices are compared to the target size.

[rmcgibbo]

Could you also add a test that fails with the current version put passes given this change?

[PolyachenkoYA]

Sure, I got one - that's why I went to find this bug in the 1st place. 2 pdb files and the code are in the email - git doesn't want to attach them here. I want to get rmsd between the 0-th chain in the uR pdb and the 1-th chain in the C pdb. They are actually the same chain (physically), so they have the same number of atoms, and atoms are identical. But they have different indices.

…

On Sun, Aug 2, 2020, 20:20 Robert T. McGibbon ***@***.***> wrote: Could you also add a test that fails with the current version put passes given this change? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUSJSEKB6HQT5ZZLWRDR6WN5DANCNFSM4PSRG7BA> .

[rmcgibbo]

Maybe just a little synthetic example that could be added to the test suite?

[PolyachenkoYA]

What do you mean by synthetic? You can do it with pdb-s wìth only 1 and 2 atoms and not 5000, but other from that - I'm not sure.

…

On Mon, Aug 3, 2020, 00:34 Robert T. McGibbon ***@***.***> wrote: Maybe just a little synthetic example that could be added to the test suite? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUR5X3ABLUIA7TQHRU3R6XLXDANCNFSM4PSRG7BA> .

[rmcgibbo]

I was thinking something tiny and/or using only the files that are already included in the repo. It shouldn't be necessary to add another 2 PDB files to this repository and force everyone to download them when they clone the repo.

[PolyachenkoYA]

well, there is not much text in these pdb-s, so I can just hardcode them and then write them to the /tmp, load from there, and then delete them. Is that ok? I mean I should build a trajectory somehow. It can be done with add_smth() many times, but I think writing real pdb-s to /tmp would be shorter and easier.

…

On Mon, Aug 3, 2020 at 5:18 PM Robert T. McGibbon ***@***.***> wrote: I was thinking something tiny and/or using only the files that are already included in the repo. It shouldn't be necessary to add another 2 PDB files to this repository and force everyone to download them when they clone the repo. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUQQNH2FN5XE4UCKYLLR63BKPANCNFSM4PSRG7BA> .

[rmcgibbo]

What I really have in mind is something like this:

def test_gh_1571():
  n_atoms_1 = 100
  n_atoms_2 = 10
  n_frames = 1
  trj_1 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_1, 3)
  trj_2 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_2, 3)
  md.rmsd(trj_1, trj_2, ...)  # something that shows the bug

[PolyachenkoYA]

def test_gh_1571(): n_frames = 1 n_atoms_1 = 1 n_atoms_2 = 2 top_1 = md.Topology() top_1.add_chain() top_1.add_residue('RS2', top_1.chain(0)) top_1.add_atom('A2', 'H', top_1.residue(0)) top_2 = md.Topology() top_2.add_chain() top_2.add_residue('RS1', top_2.chain(0)) top_2.add_atom('A1', 'H', top_2.residue(0)) top_2.add_chain() top_2.add_residue('RS2', top_2.chain(1)) top_2.add_atom('A2', 'H', top_2.residue(1)) # here the 2nd chain in the top_2 is rmsd-compatible to the one in the top_1 so we should be able to compute rsmd between them. trj_1 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_1, 3), top_1) trj_2 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_2, 3), top_2)

…

On Mon, Aug 3, 2020 at 8:24 PM Robert T. McGibbon ***@***.***> wrote: What I really have in mind is something like this: def test_gh_1571(): n_atoms_1 = 100 n_atoms_2 = 10 n_frames = 1 trj_1 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_1, 3) trj_2 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_2, 3) md.rmsd(trj_1, trj_2, ...) # something that shows the bug — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUU2342DT4DZI2T35MDR63XGTANCNFSM4PSRG7BA> .

[PolyachenkoYA]

sorry, forgot the main line def test_gh_1571(): n_frames = 1 n_atoms_1 = 1 n_atoms_2 = 2 top_1 = md.Topology() top_1.add_chain() top_1.add_residue('RS2', top_1.chain(0)) top_1.add_atom('A2', 'H', top_1.residue(0)) top_2 = md.Topology() top_2.add_chain() top_2.add_residue('RS1', top_2.chain(0)) top_2.add_atom('A1', 'H', top_2.residue(0)) top_2.add_chain() top_2.add_residue('RS2', top_2.chain(1)) top_2.add_atom('A2', 'H', top_2.residue(1)) # here the 2nd chain in the top_2 is rmsd-compatible to the one in the top_1 so we should be able to compute rsmd between them. trj_1 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_1, 3), top_1) trj_2 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_2, 3), top_2) md.rmsd(trj_2, trj_1, atom_indices=[0], ref_atom_indices=[1]) On Mon, Aug 3, 2020 at 10:39 PM Поляченко, Юрий Анатольевич < polyachenko.yua@phystech.edu> wrote:

…

def test_gh_1571(): n_frames = 1 n_atoms_1 = 1 n_atoms_2 = 2 top_1 = md.Topology() top_1.add_chain() top_1.add_residue('RS2', top_1.chain(0)) top_1.add_atom('A2', 'H', top_1.residue(0)) top_2 = md.Topology() top_2.add_chain() top_2.add_residue('RS1', top_2.chain(0)) top_2.add_atom('A1', 'H', top_2.residue(0)) top_2.add_chain() top_2.add_residue('RS2', top_2.chain(1)) top_2.add_atom('A2', 'H', top_2.residue(1)) # here the 2nd chain in the top_2 is rmsd-compatible to the one in the top_1 so we should be able to compute rsmd between them. trj_1 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_1, 3), top_1) trj_2 = md.Trajectory(np.random.RandomState(0).randn(n_frames, n_atoms_2, 3), top_2) On Mon, Aug 3, 2020 at 8:24 PM Robert T. McGibbon < ***@***.***> wrote: > What I really have in mind is something like this: > > def test_gh_1571(): > n_atoms_1 = 100 > n_atoms_2 = 10 > n_frames = 1 > trj_1 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_1, 3) > trj_2 = md.Trajectory(xyz=np.random.RandomState(0).randn(n_frames, n_atoms_2, 3) > md.rmsd(trj_1, trj_2, ...) # something that shows the bug > > — > You are receiving this because you authored the thread. > Reply to this email directly, view it on GitHub > <#1571 (comment)>, or > unsubscribe > <https://github.com/notifications/unsubscribe-auth/AJB7VUU2342DT4DZI2T35MDR63XGTANCNFSM4PSRG7BA> > . >

[rmcgibbo]

Could you add that into the pull request?

[PolyachenkoYA]

What do you mean? Should I somehow add it in the tests? If yes, how?

…

On Mon, Aug 3, 2020 at 10:43 PM Robert T. McGibbon ***@***.***> wrote: Could you add that into the pull request? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUUJYUZORSUVUMFTFN3R64HMPANCNFSM4PSRG7BA> .

[rmcgibbo]

Your pull request currently modifies the file mdtraj/rmsd/_rmsd.pyx, but could you also add that test to the file tests/test_rmsd.py? That way it will be run by the the continuous integration system.

[PolyachenkoYA]

done

…

On Mon, Aug 3, 2020 at 11:00 PM Robert T. McGibbon ***@***.***> wrote: Your pull request currently modifies the file mdtraj/rmsd/_rmsd.pyx, but could you also add that test to the file tests/test_rmsd.py? That way it will be run by the the continuous integration system. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1571 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJB7VUV44Z2ZEX4F3IRHMHTR64JMXANCNFSM4PSRG7BA> .

[rmcgibbo]

Thanks, @PolyachenkoYA