Problem with loading XTC offsets in parallel

[mkhoshle]

Hello all,

I have an MPI job where I have 48 trajectory segments in the XTC format (newtraj_0.xtc,...newtraj_47.xtc). In order to have 48 segments one can just use a small size xtc file and copy this single segment 48 times to get a list of segments. I did not put all the segments in the repository because the whole size was large (30 GB). I repeat my job several times to collect statistics.
In my MPI job I use 48 processes and each process tries to copy one of the segments (e.g. newtraj_6.xtc for rank ID 6) to another folder. This copying is to avoid cach effects. After each rank finish copying, it creates an index file for its own segment by doing print(len(mda.Universe(PSF, newtraj_rankID.xtc).trajectory)). I have a barrier to synchronize all processes.

Also I tested to see if the error is due to the corrupted offset file or not. In my code before creating the universe I do the following and all the ranks pass it successfuly,

offsets_ = []
offsets = glob.glob(os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/.*.xtc_offsets.npz'))))
for offset in offsets:
    offsets_.append(mda.coordinates.XDR.read_numpy_offsets(offset))

but some of them fail when creating the universe. i.e. here:

u = mda.Universe(PSF, XTCs)

The following helps to avoid the error, but some ranks are taking much much longer than others as a result of failure to load the offset files.

try:
    u = mda.Universe(PSF, XTC)
except IOError:
    u = mda.Universe(PSF, XTC, refresh_offsets=True)

Expected behavior
What I expect to see from my code is that each rank should create its index file and they should be able to successfully use the chain-reader to create the universe.
Note:
This error happens when I test my code on more than one node. On a single workstation (node), it works fine. The index files are being written independently by each rank. So there should not be any overlap and I also tested the code with rank0 only creating the indices. However, it does not work either. The offset files cannot be corrupted except having a race condition. Is there any method in MDAnalysis like creating the universe, getting n_frames, etc that lead to rewriting the index file?

Actual behavior
Some ranks fail due to the following error:

Traceback (most recent call last):
  File "/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/main1-ga4py-updated-rmsd-chain-reader.py", line 73, in <module>
    u = mda.Universe(PSF, XTCs)
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/core/universe.py", line 305, in __init__
    self.load_new(coordinatefile, **kwargs)
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/core/universe.py", line 535, in load_new
    self.trajectory = reader(filename, **kwargs)
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/coordinates/chain.py", line 123, in __init__
    for filename in self.filenames]
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/coordinates/core.py", line 83, in reader
    return Reader(filename, **kwargs)
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.py", line 144, in __init__
    self._load_offsets()
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.py", line 183, in _load_offsets
    data = read_numpy_offsets(fname)
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/MDAnalysis/coordinates/XDR.py", line 84, in read_numpy_offsets
    return {k: v for k, v in six.iteritems(np.load(filename))}
  File "/home/mkhoshle/miniconda2/envs/ga4py/lib/python2.7/site-packages/numpy/lib/npyio.py", line 429, in load
    "Failed to interpret file %s as a pickle" % repr(file))
IOError: Failed to interpret file '/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/files/.newtraj_41_1.xtc_offsets.npz' as a pickle

Code to reproduce the behavior
Here is the simple code to reproduce the problem:

import numpy as np
import MDAnalysis as mda
from shutil import copyfile
import glob, os
import mpi4py
from mpi4py import MPI

MPI.Init
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()

j = 1    # j is the repeat number
PSF = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/adk4AKE.psf')))
longXTC = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'traj_data_{}'.format(size),'newtraj_{}.xtc'.format(rank))))
longXTC1 = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/newtraj_{}_{}.xtc'.format(rank,j))))

copyfile(longXTC, longXTC1)
print(len(mda.Universe(PSF, longXTC1).trajectory))

comm.Barrier()
XTCs = glob.glob(os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/*.xtc'))))
u = mda.Universe(PSF, XTCs)

MPI.Finalize

Any feedback and comment will be appreciated.

Currently version of MDAnalysis

• Which version are you using? (0.17.0)
• Which version of Python (Happens with both Python 2 & 3)
• Which operating system? (CentOS release 6.7)

[jbarnoud]

Thank you for reporting the issue. Out of curiosity, can you open the offset files directly with numpy.load without MPI. If you can, can you open the trajectories in MDAnalysis without MPI? Things should not fail like that, but we need to know how much is due to MPI, and at what point do things stop working.

[richardjgowers]

FWIW we also see the above error intermittently on Travis (1 in 50? runs) which also is also a (slightly) parallel environment.

This are the relevant lines in numpy for loading the npz file:
https://github.com/numpy/numpy/blob/master/numpy/lib/npyio.py#L412

It looks like it's first failing to recognise the file as a zip file then failing to use it as a pickle.

[mkhoshle]

@jbarnoud Yes I am able to load the offsets even in parallel. As I mentioned in my issue report all the processes are able to load all the the offsets. Therefore, all of them are able to do the following:

offsets_ = []
offsets = glob.glob(os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/.*.xtc_offsets.npz'))))
for offset in offsets:
    offsets_.append(mda.coordinates.XDR.read_numpy_offsets(offset))

However, I am not sure what happens when they try to create the universe (i.e. u = mda.Universe(PSF, XTCs)) and some of the processes fail with the error above. Yes, serially one process can do u = mda.Universe(PSF, XTCs) and it does not fail. The error happens exactly at return {k: v for k, v in six.iteritems(np.load(filename))} at read_numpy_offsets in XDR.py. I tried to run this in a simple test case without using MDAnalysis at all. The error happens when I do it in parallel and as a result I think it is not MDAnalysis error. It does not seem to be numpy error either. Do you know why the file name start with .? It might be due to the way these files are being handled.
But I have no idea what should I do to avoid this error!!!

[mkhoshle]

@richardjgowers Does not seem to be a numpy problem. It should be something with XTC I think.

[orbeckst]

Did I understand you correctly in #1988 (comment) that the line

mdanalysis/package/MDAnalysis/coordinates/XDR.py

Line 84 in 5f151b4

return {k: v for k, v in six.iteritems(np.load(filename))}

fails when using MDAnalysis and that if you use only the line

offset_dict = {k: v for k, v in six.iteritems(np.load(filename))}

in parallel (under mpi4py, without using any MDAnalysis) then it fails, too? In this case, is filename the same or different for the different MPI ranks?

[orbeckst]

Do you know why the file name start with .?

Offset files are named with a leading period to make them hidden files under Unix-like operating systems. I don't think that these files are treated in any other way by the OS except that ls does not show them by default.

[mkhoshle]

Yes, not just that line. Here is the code I use outside MDAnalysis.

def read_numpy_offsets(filename):
    return {k: v for k, v in six.iteritems(np.load(filename))}

def _load_offsets(filename):
    fname = filename

    if not isfile(fname):
       ctime = getctime(filename)
       size = getsize(filename)
       try:
          np.savez(filename,
          offsets=offsets, size=size, ctime=ctime,
                         n_atoms=data["n_atoms"])
       except Exception as e:
           print('111')
           print("Couldn't save offsets because: {}".format(e))

    data = read_numpy_offsets(fname)
    ctime_ok = size_ok = n_atoms_ok = False

    try:
       ctime_ok = getctime(filename) == data['ctime']
       size_ok = getsize(filename) == data['size']
       n_atoms_ok = data['n_atoms']
    except KeyError:
       pass

    if not (ctime_ok and size_ok and n_atoms_ok):
       print("Reload offsets from trajectory\n "
                          "ctime or size or n_atoms did not match")
       ctime = getctime(filename)
       size = getsize(filename)
       try:
          print('rank {}'.format(rank))
          np.savez(filename,
          offsets=offsets, size=size, ctime=ctime,
                         n_atoms=data["n_atoms"])
       except Exception as e:
             print("Couldn't save offsets because: {}".format(e))
    
PSF = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/adk4AKE.psf')))
longXTC = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'traj_data_{}'.format(size),'newtraj_{}.xtc'.format(rank))))
longXTC1 = os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/newtraj_{}_{}.xtc'.format(rank,j))))

copyfile(longXTC, longXTC1)
print(len(mda.Universe(PSF, longXTC1).trajectory))

comm.Barrier()
offsets = glob.glob(os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/.*.xtc_offsets.npz'))))

for i in range(len(offsets)):
    fname = offsets[i]
    _load_offsets(fname)

MPI.Finalize

The error happens here:

Traceback (most recent call last):
  File "/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/test.py", line 87, in <module>
    _load_offsets(fname)
  File "/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/test.py", line 40, in _load_offsets
    data = read_numpy_offsets(fname)
  File "/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/test.py", line 24, in read_numpy_offsets
    return {k: v for k, v in six.iteritems(np.load(filename))}
  File "/home/mkhoshle/miniconda3/envs/py3/lib/python3.6/site-packages/numpy/lib/npyio.py", line 431, in load
    "Failed to interpret file %s as a pickle" % repr(file))
OSError: Failed to interpret file '/oasis/scratch/comet/mkhoshle/temp_project/MPI-traj-split/files/.newtraj_10_1.xtc_offsets.npz' as a pickle

and when I comment this part in my code, the error does not happen at all, but it happens when I do not comment it out.

try:
       ctime_ok = getctime(filename) == data['ctime']
       size_ok = getsize(filename) == data['size']
       n_atoms_ok = data['n_atoms']
    except KeyError:
       pass

    if not (ctime_ok and size_ok and n_atoms_ok):
       print("Reload offsets from trajectory\n "
                          "ctime or size or n_atoms did not match")
       ctime = getctime(filename)
       size = getsize(filename)
       try:
          print('rank {}'.format(rank))
          np.savez(filename,
          offsets=offsets, size=size, ctime=ctime,
                         n_atoms=data["n_atoms"])
       except Exception as e:
             print("Couldn't save offsets because: {}".format(e))

I believe this specific line is responsible for the error:
np.savez(filename,offsets=offsets, size=size, ctime=ctime,n_atoms=data["n_atoms"])

What happens if we remove this line? What does it do?

[orbeckst]

Yes, you are right, that line is a problem. It seems to me that in your _load_offsets() function you do not define offsets before you use it. Thus, when your code reaches

       try:
          np.savez(filename,
          offsets=offsets, size=size, ctime=ctime,
                         n_atoms=data["n_atoms"])

it will take offsets from the enclosing (global) scope where it is a list of file names (and not an array)

offsets = glob.glob(os.path.abspath(os.path.normpath(os.path.join(os.getcwd(),'files/.*.xtc_offsets.npz'))))

It then writes the filename list into the "offsets" file. Apparently, savez() is happy to do this – or it fails but the exception block just lets you move past it and then there's no offset file on disk (or an empty one with 0 size?).

(I actually don't understand why read_numpy_offsets() would not read a file that was saved successfully – I'd assume you get the filename list back instead of an array of integers.)

In any case, I'd first fix the problem with offsets not being defined properly and then see what happens. (Btw. data is also not defined inside the function and is taken from the global scope.)

[mkhoshle]

I tested the piece of code I showed in my initial post with the following line commented out in the MDAnalysis package and everything works fine with no error.

    def _load_offsets(self):
        """load frame offsets from file, reread them from the trajectory if that
        fails"""
        fname = offsets_filename(self.filename)

        if not isfile(fname):
            self._read_offsets(store=True)
            return

        data = read_numpy_offsets(fname)
        ctime_ok = size_ok = n_atoms_ok = False

        try:
            ctime_ok = getctime(self.filename) == data['ctime']
            size_ok = getsize(self.filename) == data['size']
            n_atoms_ok = self._xdr.n_atoms == data['n_atoms']
        except KeyError:
            # we tripped over some old offset formated file
            pass

        if not (ctime_ok and size_ok and n_atoms_ok):
            warnings.warn("Reload offsets from trajectory\n "
                          "ctime or size or n_atoms did not match")
#            self._read_offsets(store=True)
        else:
            self._xdr.set_offsets(data['offsets'])

The problem is that when the condition if not (ctime_ok and size_ok and n_atoms_ok): is true, then the self._read_offsets(store=True) will be called and then it rewrites the offset fille because
np.savez(offsets_filename(self.filename), offsets=offsets, size=size, ctime=ctime, n_atoms=self._xdr.n_atoms)
will be executed. This means that you will have race condition in parallel. When each process is working on a specific segment everything is fine but when each process wants to access all segments then race condition will occur if several processes happen to rewrite it at the same time.

[kain88-de]

Oh yeah if multiple processes want to access the same file you are prone to race conditions. I thought your code is written so that each process has it's own file.

[mkhoshle]

Each process creates the offset file for one segment and there is a barrier to avoid race condition. But after the offset is written when I use chain-reader then each process wants to access all segments and if the condition if not (ctime_ok and size_ok and n_atoms_ok): is true then it means that several processes can write to the offset file at the same time.

[mkhoshle]

When I do not use chain-reader and each segment is assigned to one process only, then there is no problem and no race condition.