Skip to content

Mass COM calculation#686

Closed
richardjgowers wants to merge 2 commits intodevelopfrom
RGpositions
Closed

Mass COM calculation#686
richardjgowers wants to merge 2 commits intodevelopfrom
RGpositions

Conversation

@richardjgowers
Copy link
Member

So over in #670 we mooted having a mass COM method, so I mocked this up quickly in Cython. (it's roughly 30x faster than previous)

Basically, given a list of 0 based resids, returns the center of geometry for each resid. So is the equivalent of:
np.array([r.center_of_geometry() for r in u.residues])

So there's lots of rough edges as to how the function works, including adding a mass (and charge?) weighting option so it returns CoM.

So my question is, is this as good as cython as we can write? (maybe @kain88-de or anyone better at cython than me can say)

I did a numpy typed version:
1206092

And a memory view version:
62fdf99

And couldn't see any difference in the performance

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can tell the compiler how the array is stored (Row or Column Major). That can be done by giving it the strides directly float[:, ::1 is row major for example the C/Python standard. That should give a small boost but requires that the given arrays follow that memory layout.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would keep your current version as this is less pain for someone who uses the function.

@kain88-de
Copy link
Member

What was your benchmark? I don't see anything greatly wrong here so it should be fine and generally faster then a python version.

Fix the doc string and it should be good to go.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can turn this off globally in setup.py file. check my pytraj's one

https://github.com/Amber-MD/pytraj/blob/master/setup.py#L327-L347

@richardjgowers
Copy link
Member Author

@kain88-de I'm loading the PSF/DCD system and doing

In [1]: import MDAnalysis as mda

In [2]: u = mda.Universe('adk.psf','adk_dims.dcd')

In [3]: from MDAnalysis.lib.trans import residue_positions

In [4]: indices = u.atoms.resids - 1

In [5]: %timeit out = residue_positions(u.atoms.positions, indices)
10000 loops, best of 3: 140 µs per loop

In [6]: %timeit long_way = np.array([r.center_of_geometry() for r in u.residues])
100 loops, best of 3: 3.03 ms per loop

I'll look at defining the layout and see..

@kain88-de
Copy link
Member

You don't need to define the layout. That was more as an FYI. It should be absolutely fine here as it is.

@kain88-de kain88-de deleted the RGpositions branch May 6, 2016 12:05
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants