PCA test sometimes fails...

[richardjgowers]

I get this test failing locally when I run in parallel:

=============================================================== FAILURES ===============================================================
___________________________________________________________ test_given_mean ____________________________________________________________
[gw0] linux -- Python 3.6.10 /home/richard/miniconda3/envs/mda/bin/python

pca_aligned = <MDAnalysis.analysis.pca.PCA object at 0x7f7cda848588>, u = <Universe with 3341 atoms>

    def test_given_mean(pca_aligned, u):
        pca = PCA(u, select=SELECTION, align=False,
                  mean=pca_aligned._mean).run()
>       assert_almost_equal(pca_aligned.cov, pca.cov, decimal=5)
E       AssertionError: 
E       Arrays are not almost equal to 5 decimals
E       
E       Mismatched elements: 900 / 900 (100%)
E       Max absolute difference: 1.0380358
E       Max relative difference: 106.43334522
E        x: array([[ 6.85943e-02,  6.14962e-03, -2.85235e-02, -2.46964e-03,
E                4.93089e-03, -1.45766e-02,  6.77579e-03,  6.73560e-03,
E               -6.87083e-03, -1.39547e-02,  4.34639e-03, -1.46411e-02,...
E        y: array([[ 1.10663e+00, -4.76013e-01,  2.59657e-01,  8.86885e-01,
E               -4.62178e-01,  3.13851e-01,  7.73986e-01, -4.43201e-01,
E                2.00352e-01,  4.87414e-01, -3.32454e-01,  1.01439e-01,...

analysis/test_pca.py:169: AssertionError

This is when I run pytest -n 12 from the MDAnalysisTests/ directory. If I run pytest -n 12 analysis/ it doesn't fail any more. Varying core count seems to make it fail/not fail, so it feels like something unthreadsafe is happening?

I've tried reducing the scope of the fixture and this still fails.

[IAlibay]

So looking into this more, I think this test might actually be broken. If you comment out every other test except test_given_mean, you'll have a test that fails even if you run it in serial. I suspect some kind of transfer_to_memory issue is happening here.

Trying to isolate the test (might have messed up somewhere), it looks like the main issue here is that pca_aligned._mean just returns None, so no mean value is actually passed?

import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis import align
from MDAnalysis.analysis.pca import PCA
from MDAnalysisTests.datafiles import PSF, DCD

SELECTION = 'backbone and name CA and resid 1-10'

u_aligned = mda.Universe(PSF, DCD, in_memory=True)
pca_aligned = PCA(u_aligned, select=SELECTION, align=True).run()
print(pca_aligned.cov[0][0]) # 0.0685

u = mda.Universe(PSF, DCD, in_memory=True)
pca_mean_aligned = PCA(u, select=SELECTION, align=False, mean=pca_aligned._mean).run()
print(pca_mean_aligned.cov[0][0]) # 1.106

[Balasubra]

PCA routine fails when one passes the AtomGroup to mean. The reason being the self.mean not flattened, causing a numpy broadcast error at lines 266.

u = mda.Universe(PSF, DCD)
sel_grp = u.select_atoms('name CA')
sel_grp.write('Ref.pdb')
Ref = mda.Universe('Ref.pdb')
PSF_pca = pca.PCA( u, select='name CA', mean=Ref.select_atoms('all') ).run()

To fix the bug, i flattened the self.mean in line 225 as self.mean = self._mean.positions.ravel() and now the broadcast error is resolved.

I am also thinking that it would be useful to raise an error if self._atoms.n_atoms != self._mean.n_atoms @ line 224. It may happen that the user, by mistake, input different AtomGroup to select and mean parameters.

[fiona-naughton]

So having had a look, it seems like this is a combination of things:

1. As @IAlibay pointed out, pca_aligned._mean is None (since no mean is provided to pca_algined), so the pca in this test isn't actually running with an input mean like it's supposed to (it'll just recalculate it from the trajectory again).
   The calculated mean is actually stored in pca_aligned.mean, but as an array, while the input format for reference mean is as an atomgroup, so it's a little more annoying to fix. Is there any particular reason for requiring the reference as a atomgroup rather than a position array?

2. As @Balasubra pointed out: if you did actually provide a reference mean, it breaks because the position array is assumed to be flattened; since the test is broke we never caught this, but as above should be an easy fix

3. Even aside from the issue with the mean: the test compares PCA run (on a supposedly unaligned trajectory) with align=False (pca) to one with align=True (pca_aligned) - which should be different and fail. But since they use the same universe and pca_aligned runs first, that alignment carries over to pca so the test passes unless (I presume) you have just the right setup that causes the two to run on different cores; hence the only sometimes failing. Presumably adjusting so they don't work off the same universe will fix.

[IAlibay]

Thanks for chasing this one up @fiona-naughton

I'm going to put it down as a 2.0 milestone, it's a pretty glaring issue (at least if you try to use that code path) so we should try to fix it.

@fiona-naughton if it's ok I'm going to assign you as being in charge of making sure this gets addressed?