The switch to PyTNG has been fun but also challenging, as switching to a new language at short notice has required some adjustment in both mindset and practicality. In good news, my experience with TNG itself has allowed me to rapidly make strides in changing the interface and adding test cases. I have developed a very rough prototype that reads coordinates, velocities and forces if present. This is validated by a range of new test cases and test files that match those in the C++ TNG library. The result of my handiwork can be seen in #29.
This gets the job done but is far from satisfactory for reasons I will outline below. One upshot is that there are a whole bunch of new test cases for later modifications.
The counting of “frames” by MD steps means that you have to do a whole bunch of useless blank reads. A TNG file can have as many potentially arbitrary blocks as it wants and it is difficult/wasteful to keep track of all of them. Instead we should try and spool off the blocks sequentially! That makes much much more sense as it will give us data when present only and exactly the types we want, all we need is a matching cdef cython class to hold the information and a way to pull them out of the file.
This sounds easy right? A draft PR of the start of this can be found in #32.
How exactly to spool the blocks off the file is something I still need to work on. There are lots of possible functions that seem to do similar things. At the core of most of them are twin calls to
tng_block_read_next, but I am a bit unclear about what exact state the
tng_trajectory_t struct needs to be in to call this correctly. This is because the blocks do not contain (obvious) information linking them together (I would love to be wrong here). Instead as far as I can tell the position is just kept track of with an
int64_t counter of the file position and a struct that marks the current
Watch this space.
Once we have the data off the block we will need to cast it to a numpy array. This can be done with the following code which is a bit ugly but gets the job done.
#declare mem wrapper and array cdef MemoryWrapper wrap cdef float* array = NULL #wrap python memory and bind it to the values at C level wrap = MemoryWrapper(dim1 * dim2 * sizeof(float)) # deals with ref counts, PyMemMalloc and PyMemFree values = <float*> wrap.ptr function_that_binds_ptr_2_ref_counted_float_array(&values) #declare numpy things cdef np.ndarray reshaped_npy_arr cdef npy_intp_dims dims # dimension holder cdef int err cdef int ndim = 2 #set the dims dims = dim1 dims = dim2 #bind the values to the numpy array reshaped_npy_arr = PyArray_SimpleNewFromData(nd, dims, NPY_FLOAT, wrap.ptr) Py_INCREF(wrap) #increment reference counter to take new binding into account
Keep working on the way to spool blocks off with pyTNG, may need some ideas from richard on this. Magnus has also kindly given me some feedback on my TNG++ code which I will incorporate and move forward with.