dtype of array returned by *Group.center() is not consistent

[zemanj]

Expected behaviour

ag.center(weights=None).dtype should equal ag.center(weights=np.ones(len(ag))).dtype

Actual behaviour

If weights.dtype is of higher precision than ag.positions.dtype, the resulting dtype is that of weights.

Code to reproduce the behaviour

import numpy as np
import MDAnalysis as mda
from MDAnalysisTests.datafiles import PSF, DCD
u = mda.Universe(PSF, DCD)
ag = u.atoms

print(ag.center(weights=None).dtype)  # prints dtype('float32')
print(ag.center(weights=np.ones(len(ag))).dtype)  # prints dtype('float64')

Current version of MDAnalysis

0.18.1-dev

Root of the problem

The line causing the problem is this one:

return np.average(coords, weights=weights, axis=0)

According to the numpy docs:

The return type is Float if a is of integer type, otherwise it is of the same type as a.

The problem here is that np.average() doesn't behave the way the numpy docs claim it should. I already opened an issue for this.

Other affected methods

• center_of_geometry(): returns array of dtype=np.float32
• center_of_mass(): returns array of dtype=np.float64

I see no reason why those methods should return arrays with different precisions.

Proposed solution

One could argue that the dtype of ag.center(...) should always be that of ag.positions since it computes a coordinate.
However, if we enforce that, tests for methods such as center_of_mass() fail. In fact, I fear that methods using center() might loose too much precision if the result of center() is in single precision.

I'd personally opt for always returning double precision, which could be enforced by promoting the weights if required. The part of the code computing the results if the compound parameter is not 'group' could also easily be promoted to double precision.

[orbeckst]

I think I am ok with methods that do computation return float64.

I think that as soon as you use such a result to, e.g., translate coordinates, it will down-cast to float32

c = ag.center_of_mass()
ag.positions -= c

so coordinates are not suddenly promoted to float64.

[zemanj]

tl;dr:
I suggest enforcing all operations in AtomGroup.center() to be executed in double precision, and therefore also returning resuts in double precision.

Long version:
@orbeckst Sure, assigning values to an existing ndarray won't change its precision.
As you said, doing anything such as

a = np.ones(10, dtype=np.float32)
b = np.ones(10, dtype=np.float64)
a += b

won't change a.dtype, so it will remain np.float32. However, when assigning objects, for example

a = a + b

a.dtype will be np.float64.

This does not apply to AtomGroup.positions in MDA because positions is a @property and internally, the corresponding setter only ever uses value assignments:

@positions.setter
    def positions(self, values):
        ts = self.universe.trajectory.ts
        ts.positions[self.ix, :] = values

So, even though one could naively expect type promotion, (apparent!) object assignment such as

shift = np.ones(ag.positions.shape, dtype=np.float64)
u.atoms.positions = u.atoms.positions + shift

is actually value assignment and therefore won't change u.atoms.positions.dtype.

My concern is actually a bit different, it is about the precision which is used during computation of observables such as center_of_geometry(). Currently, the culprit is the entire calculation being performed in single precision. It's not just the result that is down-cast to float32. I think we should always enforce operations involving (potentially) many floating point numbers to be executed in double precision.
For center_of_mass(), this is done automatically since masses are stored in double precision, and np.average() performs its calculations in the precision corresponding the highest one of its operands.

When executing computations in double precision, I see no reason why the results should be cast to single precision afterwards. In fact, downcasting would mean to intentionally throw away potentially useful information.

[richardjgowers]

One sticking point here is that distance_array and friends all expect float32, and fail loudly about it. This makes it annoying to chain together different steps if you keep having to up/down cast. Maybe we should just use Cython's cdef floating and have float32/64 paths.

[zemanj]

I tested some approaches to enforce double precision in *Group.center() computations.
Casting self.atoms.positions to np.float64 for weights=None seems most applicable.
If weights != None casting weights to double precision should be faster than casting coordinates.

[zemanj]

@richardjgowers This can be fixed easily by automatic type conversion: Simply replace lines such as

ref = reference.copy('C')

by

ref = reference.astype(np.float32, order='C', copy=True)

See also PR #1936

[zemanj]

BTW there are these notes in some of the lib.distances functions, stating that arrays have to be copied internally. IMHO that's only true if coords.flags['C_CONTIGUOUS'] is False, since input arrays are never modified by these functions.

[richardjgowers]

@zemanj yeah there's no reason why we have to copy.... if you can get it working without copies that would be fine too. The type conversions look good

[kain88-de]

One sticking point here is that distance_array and friends all expect float32, and fail loudly about it. This makes it annoying to chain together different steps if you keep having to up/down cast. Maybe we should just use Cython's cdef floating and have float32/64 paths.

This will be difficult with our C headers that actually implement the distance calculations (a reason why pure cython implementations are cleaner). The C functions are not parametrized on different floating point types. The easiest way without copying all functions is to switch the files to c++ mode and use templates. This would allow the functions to accept single and double precision floats.

[richardjgowers]

@kain88-de or just go pure Cython?

[kain88-de]

possible but we have to figure out how to do the optional openmp backend then. It's possible to use openmp in cython with the prange loop but I never tried to make it optional.

[zemanj]

I'd opt for C++ templating here:
I have an update for lib.distances in the pipeline that speeds up RDF calculations by a factor of 3 to 20, depending on the number of OpenMP threads. It uses a more efficient memory layout for coordinates which allows to maximize cache locality and furthermore helps the compiler to auto-vectorize a lot of loops without requiring explicit SIMD instructions in the code. It also circumvents the use of distance arrays for RDF computations, which currently impose a rather strict limit on the number of particles one can compute RDFs for (one simply runs out of RAM).
Most part of this would be pretty much impossible with a pure Cython implementation, and porting my current state of work to C++ should be rather straightforward. If you want to have a preview, you can check out the lib_distances branch from my MDA fork.