XDR file seeks and tells working again on large files (closes #677)

[mnmelo]

Did not add tests, as they'd require a >2GB test file (a reason why this bug cropped up again after having been addressed in the past).

Also tested in a 32 bit vagrant box. For the record, the only failing test there was:

======================================================================
FAIL: MDAnalysisTests.coordinates.test_lammps.TestLammpsDataTriclinic.test_triclinicness
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/usr/local/lib/python2.7/dist-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/vagrant/mdanalysis/testsuite/build/lib.linux-i686-2.7/MDAnalysisTests/coordinates/test_lammps.py", line 249, in test_triclinicness
    assert_(self.u.dimensions[5] == 120.)
  File "/usr/local/lib/python2.7/dist-packages/numpy/testing/utils.py", line 53, in assert_
    raise AssertionError(smsg)
AssertionError:

As discussed in #677, I'd vote for restoring the old behavior of _xdr.tell(), mostly for debugging convenience.

[kain88-de]

unused import

[jbarnoud]

I am very sorry to point out that the fix does not seem to work on my test case. EDIT: That was my fault, see comment bellow.

In [10]: u = mda.Universe(top, traj)

#Until frame 25138, it works
In [11]: target_frame = 25138
In [12]: u.trajectory[target_frame]
Out[12]: < Timestep 25138 with unit cell dimensions [ 241.50689697  243.38421631   87.10479736   90.           90.           90.        ] >
In [13]: print(u.trajectory._xdr.offsets[target_frame], len(u.atoms))
(4294835832, 45904)

# From frame 25139, it fails
In [14]: target_frame = 25139
In [15]: u.trajectory[target_frame]
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-15-fa427d781897> in <module>()
----> 1 u.trajectory[target_frame]

/home/jon/dev/mdanalysis/package/MDAnalysis/coordinates/base.pyc in __getitem__(self, frame)
   1162         if isinstance(frame, int):
   1163             frame = apply_limits(frame)
-> 1164             return self._read_frame(frame)
   1165         elif isinstance(frame, (list, np.ndarray)):
   1166             def listiter(frames):

/home/jon/dev/mdanalysis/package/MDAnalysis/coordinates/XDR.pyc in _read_frame(self, i)
    140         self._xdr.seek(i)
    141         self._frame = i - 1
--> 142         return self._read_next_timestep()
    143 
    144     def _read_next_timestep(self, ts=None):

/home/jon/dev/mdanalysis/package/MDAnalysis/coordinates/XDR.pyc in _read_next_timestep(self, ts)
    148         if ts is None:
    149             ts = self.ts
--> 150         frame = self._xdr.read()
    151         self._frame += 1
    152         self._frame_to_ts(frame, ts)

/home/jon/dev/mdanalysis/package/MDAnalysis/lib/formats/xdrlib.pyx in MDAnalysis.lib.formats.xdrlib.XTCFile.read (MDAnalysis/lib/formats/xdrlib.c:7672)()
    590                                       <rvec*>xyz.data, <float*> &prec)
    591         if return_code != EOK and return_code != EENDOFFILE:
--> 592             raise RuntimeError('XTC Read Error occured: {}'.format(
    593                 error_message[return_code]))
    594 

RuntimeError: XTC Read Error occured: magic

[jbarnoud]

Nevermind, I did not clean the offset... (Isn't the validity of the offset supposed to be checked by trying to seek to the last frame? If so, it should not have used the old offset.)

[kain88-de]

Isn't the validity of the offset supposed to be checked by trying to seek to the last frame?

See #656 and #631. Currently the offsets aren't validated because I didn't find the old approach robust enough. In #656 I'm working on a new approach that the offsets are reloaded once we notice that they fail.

[mnmelo]

I left the offset types as off_t, as per #677.

I added the low-level _bytes_seek and _bytes_tell methods.

As had already been discussed in #441, xtc_seek was returning an errorcode larger than the error message list. This caused a cryptic IndexError when trying to raise a RuntimeError for a failed seek. The larger error code was returned because the original xdrfile implementation doesn't have a specific error for it. Since the fseeko call reports its error via errno, I now pass that value instead. It is printed as is with the RuntimeError exception, and no lookup is done in the XDR error list.

For reference I got errno=22 when seeking beyond the 16TB mark (well before the 64-bit overflow). This reflects an overflow of the filesystem allocation units (ext3 with 4K clustersize, in my case), but the errno itself is not very informative, since it is the same one would get if attempting to go to negative offsets, or passing a whence argument other than SEEK_SET, SEEK_CUR or SEEK_END. (Let's hope that whoever gets bitten by this error in the future finds this message)

[kain88-de]

Use a dictionary

[kain88-de]

the same goes here. We don't need to create variables with the correct type for return values of functions. Cython takes care of that. But for arguments to C functions it is good to avoid ambiguity.

[kain88-de]

Looks good to me besides the comments.

[mnmelo]

@kain88-de, for consistency, I was going through the code and there are many other initializations as int for variables that'll receive return codes. If there's an overhead, even small, I'd say let's make them all cdef, no?
The code is a bit inconsistent in that return codes are sometimes compared to zero, sometimes to EOK (which turns out to be equivalent...). Should we settle on comparing to EOK? Also, since error codes are cdef enum does it mean there'll be further overhead when comparing a Python object to them (more case for keeping return code variables cdef).
What do you say?

[kain88-de]

If there's an overhead, even small, I'd say let's make them all cdef, no?

Nope I just checked cython generates these variables with the correct type for us. So no overhead.

The code is a bit inconsistent in that return codes are sometimes compared to zero, sometimes to EOK (which turns out to be equivalent...). Should we settle on comparing to EOK?

Yes I have been inconsistent there. It should be compared to EOK for clarity. You can do that if you want. Otherwise I'll do it myself.

Also, since error codes are cdef enum does it mean there'll be further overhead when comparing a Python object to them (more case for keeping return code variables cdef).

I doubt that. You can check the produced C-code yourself by compiling the pyx with cython --annotate xdrlib.pyx. That will produce a html showing the cython-code. The more yellow the code the more expensive it is. The more yellowish a line is the more python specific stuff is called. If you click on a line it will show the C-code produced for it.

[mnmelo]

Ok, then I clean all these cdefs only for return codes, even if they're later compared to error codes.

[kain88-de]

Ok, then I clean all these cdefs only for return codes, even if they're later compared to error codes.

Thank you.

[kain88-de]

@mnmelo. Looks good. You can merge it if you are done.