Cython xdrlib

[kain88-de]

See #419

This is the progress I made on the cython port of xdrlib. Most tests are passing, some are failing because currently XTC/TRR-File has less features then the old version. I'll add them later, see list below). The reimplemented Readers/Writers are now just a minimal subclass of base.Reader. I'd like to imagine this is easier to understand then the previous code

Things that have happened

short list of what I did.

• add updated xdrlib source
• wrap unmodified xdrlib with python-like file objects in cython
• add tests for file XTCFile and TRRFile classes
• use XTCFile in XDR-Reader/Writer
• use TRRFile in XDR-Reader/Writer

TODO

• clean up XTC-Reader/Writer
• clean up TRR-Reader/Writer
• use a combined base class for XTC/TRR reader|writer
• add offsets and proper seek to XTCFile
• add offsets and proper seek to TRRFile
• use baseclass for XTCFile/TRRFile
• support gaps in TRR
• re enable some of the deactivated tests
• add doc-strings for XTC-Reader/Writer
• add doc-strings for TRR-Reader/Writer
• add some more test to get ~100% coverage for the new classes
• reader.next should return second frame after initialization
• implement persistent offsets using np.savez, see benchmarks
• check docs
• add info to changelog

[kain88-de]

Was there anything really special about the old XTCTimestep? I didn't could so far to everything with the base.Timestep. If there is nothing special I'll remove these test as the class doesn't exist anymore

[richardjgowers]

https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/coordinates/xdrfile/core.py#L97

Looks like some stuff for keeping the status of the Reader and _frame which is the trajectories opinion of the frame, rather than MDA's

[richardjgowers]

Oh, and the unitcell isn't standard either. Timestep._unitcell should be the native format representation of the box, gromacs has 3 vectors

[kain88-de]

Well I already use the native format with the conversion functions defined in 'lib.mdamath'. If this is all then we don't have a reason that to keep that special Timestep class

[richardjgowers]

The problem with that is if you write out your unit cell again you might end up with rounding errors (89.999 angles) , which then isn't a cuboid any more.

[kain88-de]

According to this code that is actually what MDAnalysis is doing currently.

mdanalysis/package/MDAnalysis/coordinates/xdrfile/core.py

Line 245 in 96e3ce9

unitcell = self.convert_dimensions_to_unitcell(ts).astype(np.float32) # must be float32 (!)

in a convoluted way currently we are directly using for every save

box = triclinic_vectors(triclinic_box(box_vectos))

But I'll definitely run a check how often that can be done until the errors sum.

[richardjgowers]

I think another thing that's going to have to be checked is performance. Probably just something simple like iterating through a huge trajectory and seeing how fast it goes

[kain88-de]

I'll do a performance check once I implemented the old seek code again. For proper performance tests we should have a least two files. One with a large number of frames and another with a large number of atoms, I'm think that that memory allocations could be a problem with a large number of atoms per frame. But I'll have to test that first.

[richardjgowers]

Is this going to put the trajectory before frame 0? The Reader should read the first frame, so that calling next reads the second frame

u = mda.Universe(etc)

u.trajectory.next()  # should be second frame

[kain88-de]

Ah OK I didn't know that. I added a TODO for me and I will also add tests for all readers.

[kain88-de]

#474 definitely was a good idea. The new tests are already catching stuff that I didn't notice before.

[richardjgowers]

Yeah I was going to move a few Readers over into the new tests and I'm half expecting to find a few nice little bugs.

Don't worry about the unitcell thing too much, I'm just wary of boxes having an angle of 89.99. There's a huge difference between that and 90.0.

[kain88-de]

OK finally all the tests are passing, old and new! This means that just the optimized seeking is left to do + some cleaning up.

As far as performance goes this isn't doing to bad. Version 0.12.1 needs 24s to run the xdr tests and this needs 28s. But my hope is that this will go away once I optimized the seeking behavior

[kain88-de]

BTW is there something super special about the single frame readers? What is the expected behavior of the Readers when the trajectory only contains 1 frame?

[richardjgowers]

I think with the SingleFrameReader, what was happening previously was each Reader had an implementation of __getitem__ etc so I wrote the class to unify that. It might be possible to use Reader for them too now that's also been cleaned up

[kain88-de]

So I'm finally trying to read the offsets from the file in the same way we used to before. But I can't quite seem to make if possible. This code compiles but doesn't interface with the correctly with the C xdrlib. offsets is still a NULL-ptr after calling the read_xtc_n_frames. This shouldn't have happened. Anyone has an idea of why?

[kain88-de]

I checked the memory is allocated in the xdrlib but the pointer at the cython level is never updated.

[kain88-de]

This [SO}(http://stackoverflow.com/questions/1398307/how-can-i-allocate-memory-and-return-it-via-a-pointer-parameter-to-the-calling#1398321) posts answers it. If I want to change the address a pointer points to in a function I have to pass a pointer to that pointer. Ok so offset reading works I just need to include it in the API (YEAH)

[kain88-de]

I'm hitting a wall again here. Every time I'm calling this function from python I get SystemError: error return without exception set. I have no clue where that comes from. The xdrlib functions seem to be called correctly if I add printf statements there. As far as I can make it out there are no errors occuring. @mnmelo did you see similar things implementing the seeking in SWIG?

[mnmelo]

Hmm... didn't bump into that one. My first reaction was to blame the frame type, as they always must be 64bit (large file access, and the such).

Another thing I noticed is that, in xdrfile.c/xdrfile.h the exdr_message is being built with exdrNR-1 members. This might be related, since exdrNR is out-of-bounds, and that's what seek returns on error. Maybe you can change xdrstdio_setpos or xdr_seek to return another error (though none of the available ones really fit it).

[kain88-de]

I overlooked exdr_message so far but that will be useful to get better errors messages. So i'm sure that is not the problem.

[mnmelo]

Yea... I couldn't also tell it would be the cause, but do flag it as problematic because exdrNR shouldn't be used for an error code and it will break when we attempt to retrieve the error message from it.
Anyway, a Google search led me to http://pythonextensionpatterns.readthedocs.org/en/latest/exceptions.html, where it is pointed out that somewhere a function is returning NULL without setting an exception. None of the seek functions can return a NULL...

[kain88-de]

So I think I'm done with the largest part of the implementation. The offsets are now also calculated directly in the xdrfile-library.

I've started to do some benchmarks. The test-suite runs in the same time. But when I try to benchmark bigger files I don't get conclusive results. I have a gut feeling that this version is faster but I can't really prove it right now. It would be nice if others could also benchmark and test this code. When I have time this week I'll run some more tests on different systems.

1.7 GB XTC

MDA v0.12.1

In [1]: from MDAnalysis.coordinates.XTC import XTCReader

In [2]: %timeit -n1 -r1 XTCReader('all_bb.xtc')
1 loops, best of 1: 1.91 s per loop

MDA this PR

In [1]: from MDAnalysis.coordinates.XTC import XTCReader

In [2]: %timeit -n1 -r1 XTCReader('all_bb.xtc')
1 loops, best of 1: 597 ms per loop

150 MB XTC

MDA v0.12.1

In [1]: from MDAnalysis.coordinates.XTC import XTCReader

In [2]: %timeit  XTCReader('test.xtc')
1 loops, best of 3: 206 ms per loop

MDA this PR

In [1]: from MDAnalysis.coordinates.XTC import XTCReader

In [2]: %timeit  XTCReader('test.xtc')
The slowest run took 4.08 times longer than the fastest. This could mean that an intermediate result is being cached 
1 loops, best of 3: 50.2 ms per loop

[kain88-de]

Also I'm not sure why the travis tests are failing. The files exists in the repository and when I do a fresh checkout everything works on my laptop.

[kain88-de]

@richardjgowers I added the changelog information. You can merge this if you like.

[richardjgowers]

Sped up (exception in english)

[kain88-de]

fixed.

[richardjgowers]

Awesome. So afaik it's just DCD and lots of six things that need doing now?

[kain88-de]

@richardjgowers hopefully yes. But the dcd will be a big one again.