update qcprot to cython.floating types

[richardjgowers]

use memoryviews throughout for potential performance

allows either float32 or float64 inputs

Related #3927

Changes made in this Pull Request:

• qcprot functions take either float32 or float64 inputs
• modernised qcprot to use memoryviews

PR Checklist

• Tests?
• Docs?
• CHANGELOG updated?
• Issue raised/referenced?

Developers certificate of origin

• I certify that this contribution is covered by the LGPLv2.1+ license as defined in our LICENSE and adheres to the Developer Certificate of Origin.

---

📚 Documentation preview 📚: https://mdanalysis--4273.org.readthedocs.build/en/4273/

[github-actions]

Linter Bot Results:

Hi @richardjgowers! Thanks for making this PR. We linted your code and found the following:

Some issues were found with the formatting of your code.

Code Location	Outcome
main package	✅ Passed
testsuite	⚠️ Possible failure

Please have a look at the darker-main-code and darker-test-code steps here for more details: https://github.com/MDAnalysis/mdanalysis/actions/runs/6066435327/job/16457020106

---

Please note: The black linter is purely informational, you can safely ignore these outcomes if there are no flake8 failures!

[richardjgowers]

Still need to do some tests here. rms.RMSD looks like it's converting/copying everything up to float64, but I imagine that not forcing this would be some sort of precision-regression. Instead I can probably write some tests that using lib.qcprot directly with floats/doubles work as expected and we can have some future discussion about precision and dtypes in the future.

[IAlibay]

This looks reasonable to me. Is there an associated issue this is working towards? (edit: I know there were a few out there that had been opened about fp precision).

[richardjgowers]

@IAlibay I've tagged a related issue, I think we should be moving towards numpy-like behaviour where output dtype matches input dtype and either float or double are allowed. A future discussion is the next major release of MDA will probably have a single precision option, since we're able to get 2x performance & 1/2 memory footprint for this.

The memoryview change is just a best-practice and we should gradually modernise to this as it's documented as the recommended route for best performance.

[IAlibay]

Odd, I thought you had conttributed in 2.6 for the GSD stuff, maybe that was 2.5... time flies by.

[richardjgowers]

this branch looks like it was a bug. at this point the RMSD has been calculated and the code has failed to determine a suitable rotation matrix, returning None doesn't seem to make sense

[richardjgowers]

e.g. compared to this implementation in BioPython: https://github.com/biopython/biopython/blob/master/Bio/PDB/qcprot.py#L190

[codecov]

Codecov Report

Patch and project coverage have no change.

Comparison is base (f50a097) 93.40% compared to head (a6d45dd) 93.40%.

Additional details and impacted files

@@            Coverage Diff             @@
##           develop    #4273     +/-   ##
==========================================
  Coverage    93.40%   93.40%             
==========================================
  Files          170      184     +14     
  Lines        22250    23358   +1108     
  Branches      4071     4071             
==========================================
+ Hits         20783    21818   +1035     
- Misses         951     1024     +73     
  Partials       516      516             

see 14 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[tylerjereddy]

I remember that a while back I tried to improve performance in SciPy by switching from the old NumPy buffer syntax, to memoryviews, and it actually caused a performance regression: scipy/scipy#12379 (comment)

Anyway, you may want to measure that you remain at least performance neutral, my experience with following that best practice hasn't historically been great, but hopefully the situation is better now.

[tylerjereddy]

Maybe use assert_allclose for new test code per the Note at https://numpy.org/doc/stable/reference/generated/numpy.testing.assert_array_almost_equal.html ?

[orbeckst]

The RMSD functions are quite time-critical for many applications so if there's any chance for a performance impact then please provide some benchmark data before merge. Thanks!

[orbeckst]

About float vs double. I remember vaguely that during one of the REU projects we played around with single vs double precision and arrived at the conclusion that float32 is insufficient for the algorithm to properly converge. Of course, I don't remember any details and have no idea if there are any notes anywhere. Primarily I want to say that it simply might not be as easy to say that we just run the algorithm at whatever dtype the rest of the code wants to run and that this will require careful testing.

[richardjgowers]

@orbeckst benchmarks

from MDAnalysisTests.datafiles import PSF, DCD
import MDAnalysis as mda
from MDAnalysis.lib import qcprot
import numpy as np

u = mda.Universe(PSF, DCD)

prot = u.select_atoms('protein')

w = (prot.masses / np.mean(prot.masses)).astype(np.float64)

ref = prot.positions.astype(np.float64)
u.trajectory[1]
conf = prot.positions.astype(np.float64)

rot = np.zeros(9, dtype=np.float64)
A = np.zeros(9)

N = len(prot)

then:

%%timeit

qcprot.CalcRMSDRotationalMatrix(ref, conf, len(prot), rot, w)

old code:
10.1 µs ± 46.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
new code:
9.97 µs ± 29.1 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

The number fluctuate a bit, I think there's no performance gain/loss.

I was able to reproduce the memoryview regression that @tylerjereddy flagged up, but I think it's an order of magnitude smaller than the numbers measured here, so it won't show up. It's something to keep in mind if we're writing code that has a lot of function calls (there's 3 here).

re: precision, yes I had to make the qcp function use doubles throughout to keep it working. The InnerProduct can use either float or double now. Strangely these take the same amount of time, which means the compiler isn't able to do anything clever here (where it ought to be) so probably some future digging. This does some ground work at least.