Fast grid-based neighbor search implementation.#1996
Fast grid-based neighbor search implementation.#1996seb-buch wants to merge 35 commits intoMDAnalysis:developfrom
Conversation
…be added again for conditional compilation
…be added again for conditional compilation
Pre-converting np.resize to PyMem stuff
…d after merge/rebase
orbeckst
left a comment
There was a problem hiding this comment.
Lots of good work.
I am not sure that I am 100% qualified to review so take my comments with a grain of salt.
I hope that some other @MDAnalysis/coredevs will have a look, too.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| Neighbor search library --- :mod:`MDAnalysis.lib.grid` | ||
| ====================================================== | ||
|
|
||
| This Neighbor search library is a serialized Cython port of the NS grid search implemented in GROMACS. |
There was a problem hiding this comment.
I want to make sure we properly attribute:
- link to Gromacs website
- cite Gromacs paper in which they describe the neighbor search
- add "and published under the LGPL)"
There was a problem hiding this comment.
A bit more details on what the idea is and why it is needed would be good.
Usage example would be great, especially as this is a versatile tool that people might want to use in own analysis.
There was a problem hiding this comment.
Add functions and classes for which docs should be generated.
There was a problem hiding this comment.
If we cite then please add a reference to 'Understanding Molecular Dynamics Simulations' by Frenkel. The algorithm used is called a cell-list and the book is honestly the only reference I have found for this. It's described in the appendix together with verlet-lists
There was a problem hiding this comment.
Thanks for your comments,
The documentation is clearly lacking and needs to be written. As well as the reference to the original code/algorithm. I must admit that I wanted to push the code so @ayushsuhane can compare with his solution.
From what I have tested grid NS is faster when you perform NS on all atoms but the tree-base approach may still be faster for small calculation (eg you only neighbors for just a few residues)
There was a problem hiding this comment.
Getting it out for @ayushsuhane is great, I very much agree with you. (I wouldn't be doing my job if I weren't pointing out anything that should, in my opinion, be addressed at some point down the line.)
|
|
||
| This Neighbor search library is a serialized Cython port of the NS grid search implemented in GROMACS. | ||
| """ | ||
|
|
There was a problem hiding this comment.
Can you add a comment to the code as to which Gromacs source code files inspired this code? Which version of Gromacs? Add a note on the license of the file (probably LGPL) – best is to copy the license header from the Gromacs file if you took a lot of the code from there.
|
|
||
| # Useful Functions | ||
| cdef real rvec_norm2(const rvec a) nogil: | ||
| return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ] |
There was a problem hiding this comment.
PEP8 will likely not like the formatting here and elsewhere.
Run it locally and clean up...
There was a problem hiding this comment.
(Well, I guess pep8 does not check cython files, does it?)
There was a problem hiding this comment.
Indeed, AFAIK, pep8 just cares about pure python files.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| if use: # Accept this shift vector. | ||
| if self.c_pbcbox.ntric_vec >= MAX_NTRICVEC: | ||
| with gil: | ||
| print("\nWARNING: Found more than %d triclinic " |
There was a problem hiding this comment.
Can we get rid of print and use regular logging/warning? You already grabbed the GIL so it shouldn't be an issue to use warnings.warn().
There was a problem hiding this comment.
Good point. Specially here since the only time I see this warning, it was triggered by a bug.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| "your box.") | ||
| print(np.array(box)) | ||
|
|
||
| for i in range(self.c_pbcbox.ntric_vec): |
There was a problem hiding this comment.
All this output for just a warning?? Shouldn't this rather be for a true exception?
There was a problem hiding this comment.
To be honest, I don't know. In the original code, the warning is just printed out with more consequences. (See: https://github.com/aar2163/GROMACS/blob/05ff1bdbc33224ec6b7bb1d328abd0a4b3caf156/src/gmxlib/pbc.c#L454)
I never managed to trigger it anyway.
| return searcher.search(ids) | ||
|
|
||
|
|
||
| def test_pbc_badbox(): |
There was a problem hiding this comment.
better write this as a parametrized test with pytest.mark.parametrize to catch each individual case.
There was a problem hiding this comment.
I didn't know this one! It will make the test code more efficient
There was a problem hiding this comment.
check the docs. This is one feature of pytest I really like.
| pbcbox.dx(bad, a) | ||
| pbcbox.dx(a, bad) | ||
|
|
||
| assert_equal(pbcbox.dx(a, b), dx) |
There was a problem hiding this comment.
use assert_almost_equal for floats.
| pbcbox.dx(a, bad) | ||
|
|
||
| assert_equal(pbcbox.dx(a, b), dx) | ||
| assert_allclose(pbcbox.distance(a, b), np.sqrt(np.sum(dx*dx)), atol=1e-5) |
There was a problem hiding this comment.
I think we have also been using assert_almost_equal in these cases, you can set decimals.
|
|
||
| pbcbox = nsgrid.PBCBox(box) | ||
|
|
||
| assert_allclose(pbcbox.put_atoms_in_bbox(coords), results, atol=1e-5) |
There was a problem hiding this comment.
One bad habit of mine... You're right assert_almost_equal is better suited
| with pytest.raises(TypeError): | ||
| nsgrid.FastNS(None, 1) | ||
|
|
||
| def test_nsgrid_badcutoff(universe): |
kain88-de
left a comment
There was a problem hiding this comment.
I didn't check to understand the algorithms used yet. I did notice a few issues with memory management.
|
|
||
| # Useful Functions | ||
| cdef real rvec_norm2(const rvec a) nogil: | ||
| return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ] |
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| if self.debug: | ||
| print("Total number of pairs={}".format(npairs)) | ||
|
|
||
| # ref_bead = 13937 |
There was a problem hiding this comment.
this can be removed. I don't see how this can be useful
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| cellindex_probe = self.grid.coord2cellid(probe) | ||
|
|
||
| if cellindex == cellindex_probe and xi != 1 and yi != 1 and zi != 1: | ||
| # if self.debug and debug: |
There was a problem hiding this comment.
either remove all debug code or uncomment it.
There was a problem hiding this comment.
These comments are remains from some bug hunts... I will remove them
| cdef bint prepared | ||
| cdef NSGrid grid | ||
|
|
||
| def __init__(self, u, cutoff, coords=None, prepare=True, debug=False, max_gridsize=5000): |
There was a problem hiding this comment.
it doesn't need a universe. Just use a coords and box argument like the other functions in lib
|
|
||
| with nogil: | ||
| for i in range(size_search): | ||
| if memory_error: |
There was a problem hiding this comment.
where is this value ever changed? I see all the checks and they sound nice but I don't see where you set this value
There was a problem hiding this comment.
What do you mean? when memory error is set to True? it is if the NSResults.resizefails to reallocate memory.
There was a problem hiding this comment.
I found it later. It's a bit hidden away in the nested loops.
There was a problem hiding this comment.
Yes, unfortunately, I do not see any cleaner way to propagate the memory allocation failure without requiring the GIL (and thus making the code slower)
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| with gil: | ||
| self.beadids = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.size * self.nbeads_per_cell) #np.empty((self.size, nbeads_max), dtype=np.int) | ||
| if not self.beadids: | ||
| raise MemoryError("Could not allocate memory for NSGrid.beadids ({} bits requested)".format(sizeof(ns_int) * self.size * self.nbeads_per_cell)) |
There was a problem hiding this comment.
This is a memory leak! The array allocated to beadcounts is not released if this malloc call fails
There was a problem hiding this comment.
Good call! I did not catch this one when testing as PyMem_Malloc does not fail without a bug or a humoungous system
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| cdef ns_int *beadcounts = NULL | ||
|
|
||
| # Allocate memory | ||
| beadcounts = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.size) |
There was a problem hiding this comment.
I would replace this with cdef ns_int_t[:] beadcounts = np.empty(self.size, dtype=ns_int). Then python does proper reference counting and you avoid the memory error mentioned earlier.
There was a problem hiding this comment.
My example code assumes that you have something like
ns_int_t = np.int64_t
ns_int = np.int64
at the beginning of the file to distinguish between the int type and the numpy variable. It's common convention also used in the cython docs and examples.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| self.nbeads_per_cell = self.nbeads[cellindex] | ||
|
|
||
| # Allocate memory | ||
| with gil: |
There was a problem hiding this comment.
I would remove the indentation here and open another nogil block below.
| self.beadids = NULL | ||
| self.cellids = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.ncoords) | ||
| if not self.cellids: | ||
| raise MemoryError("Could not allocate memory from NSGrid.cellids ({} bits requested)".format(sizeof(ns_int) * self.ncoords)) |
There was a problem hiding this comment.
another memory leak. If allocating previous arrays works their memory will not be released again in the python process. I assume there are more memory leaks like this in the code.
There was a problem hiding this comment.
I don't think there is a memory leak as other "self.array_variables" will be deallocated by dealloc when NSGrid is thrown away
There was a problem hiding this comment.
Is this also true when the call to __init__ fails? I believe you for other functions but for __init__ I do not know what the interpreter will do. The question is also if some of the allocations should go into a __cinit__ function.
There was a problem hiding this comment.
I think that dealloc is always called, so even when init failed... but I will check that
There was a problem hiding this comment.
According to an old mailing list entry this dealloc should be called and work fine. So no memory leak. But according to the docs the __init__ method can be called twice.
Under some circumstances it is possible for init() to be called more than once or not to be called at all, so your other methods should be designed to be robust in such situations.
If we call malloc twice that would be a memory leak. More troubling is that the memory might not be initialized at all.
| self.pairs_buffer = None | ||
| self.pair_coordinates_buffer = None | ||
|
|
||
| def __dealloc__(self): |
There was a problem hiding this comment.
by the garbage collection when the class is deleted from memory
|
Hi Seb, I also made a benchmark with the current state of gridsearch except I changed the dependency from MDAnalysis.universe to box dimensions. While grid search is fastest in most of the cases, tree structure becomes advantageous for smaller cutoff distances. I think it is mainly because of huge number of cells for smaller distances. One trick could be to limit the maximum number of cells and then search for the desired radius. Also, correct me if I am wrong, the code in its current state cannot deal with non periodic systems, which I believe also remains to be included. |
|
Thanks Ayush for the benchmarks, I glad to see that they reflect the quick and dirty ones I did while testing the improvements of the code. |
|
I haven’t found any information for init and deallocation Problems in the
cython docs. But they say that init can in some circumstances be called
twice and implementations should be good with that. Do we get a memory leak
when we call PyMalloc twice?
Unfortunately they don’t mention when init is called twice. I think we also
have to check our other cython code for this.
…On Thu 19. Jul 2018 at 09:44, Sébastien Buchoux ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In package/MDAnalysis/lib/nsgrid.pyx
<#1996 (comment)>
:
> + self.cutoff, self.size, <ns_int> (ncoords / self.size)
+ ))
+ print("NSGrid: Size={}x{}x{}={}".format(self.ncells[XX], self.ncells[YY], self.ncells[ZZ], self.size))
+
+ self.cell_offsets[XX] = 0
+ self.cell_offsets[YY] = self.ncells[XX]
+ self.cell_offsets[ZZ] = self.ncells[XX] * self.ncells[YY]
+
+ # Allocate memory
+ self.nbeads = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.size)
+ if not self.nbeads:
+ raise MemoryError("Could not allocate memory from NSGrid.nbeads ({} bits requested)".format(sizeof(ns_int) * self.size))
+ self.beadids = NULL
+ self.cellids = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.ncoords)
+ if not self.cellids:
+ raise MemoryError("Could not allocate memory from NSGrid.cellids ({} bits requested)".format(sizeof(ns_int) * self.ncoords))
I think that *dealloc* is always called, so even when *init* failed...
but I will check that
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1996 (comment)>,
or mute the thread
<https://github.com/notifications/unsubscribe-auth/AEGnVmg_I7ipP9Xf2s-NqqjuwG6XyEzeks5uIDjigaJpZM4VU11B>
.
|
| def __init__(self, real[:,::1] box): | ||
| self.update(box) | ||
|
|
||
| cdef void fast_update(self, real[:,::1] box) nogil: |
There was a problem hiding this comment.
In general, because this isn't a parallel module, all references to gil can be removed right?
There was a problem hiding this comment.
I would leave the GIL code in – you never know under which circumstances someone is trying to run the code and having the GIL/NOGIL here seems like valuable information, i.e., someone has thought about the code. Seems a shame to throw this knowledge away.
| self.fast_update(box) | ||
|
|
||
|
|
||
| cdef void fast_pbc_dx(self, rvec ref, rvec other, rvec dx) nogil: |
There was a problem hiding this comment.
So this is a different algorithm to what we use for all other minimum image calculations. It looks like it doesn't have to search all images and pick the minimum (which is what our triclinic pbc has to do) in which case if this works just as well then it's an improvement.
Either way, I don't like having a different minimum image solution in this class, we should use the same throughout MDA. So either this is better, in which case implement it everywhere, or it's wrong, in which case use the one in calc_distances.h
There was a problem hiding this comment.
It would be good to test this against the calc_distances implementation for triclinic cells with small angles and random data. I do understand how this works and that it can be faster then the full image search It would also work independent of box type. As @richardjgowers says this should use the code in calc_distances.h and replace it if this turns out to be faster.
| raise ValueError("Not 3 D coordinates") | ||
| return self.fast_distance(&a[XX], &b[XX]) | ||
|
|
||
| cdef real[:, ::1]fast_put_atoms_in_bbox(self, real[:,::1] coords) nogil: |
There was a problem hiding this comment.
I think this is identical to lib.distances.pack_into_box
package/MDAnalysis/lib/nsgrid.pyx
Outdated
|
|
||
| # Preallocate memory | ||
| self.allocation_size = search_ids.shape[0] + 1 | ||
| self.pairs = <ipair *> PyMem_Malloc(sizeof(ipair) * self.allocation_size) |
There was a problem hiding this comment.
I thought __cinit__ was for allocating C memory
| return self.npairs | ||
|
|
||
| cdef int resize(self, ns_int new_size) nogil: | ||
| # Important: If this function returns 0, it means that memory allocation failed |
There was a problem hiding this comment.
This is counterintuitive, most other things returning 0 is no errors
There was a problem hiding this comment.
we could use some constants in the beginngin like OK and ERROR.
| cdef real[:] coord_i, coord_j | ||
| from collections import defaultdict | ||
|
|
||
| indices_buffer = defaultdict(list) |
There was a problem hiding this comment.
We can use c++ vectors instead here
| self.coordinates_buffer.append(np.array(coords_buffer[elm])[sorted_indices]) | ||
| self.distances_buffer.append(np.sqrt(dists_buffer[elm])[sorted_indices]) | ||
|
|
||
| def get_indices(self): |
There was a problem hiding this comment.
get_X methods aren't very Python, could just use a property
|
I found this library for memory allocation: https://github.com/explosion/cymem It seems like it takes care of a lot of book keeping and error checks which would simplify a lot of this code. I'll try and find some time to see if it works inside these objects |
kain88-de
left a comment
There was a problem hiding this comment.
I gave the code a more thorough read today. The code is generally of good quality. There are still some issues left we should address now.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| rvec mhbox_diag | ||
| real max_cutoff2 | ||
| ns_int ntric_vec | ||
| ns_int[DIM] tric_shift[MAX_NTRICVEC] |
There was a problem hiding this comment.
Is this a leftover from somewhere else? I see that you set those values but I can't find where you use them. Could you please show me.
package/MDAnalysis/lib/nsgrid.pyx
Outdated
| PyMem_Free(self.pair_distances2) | ||
|
|
||
| cdef int add_neighbors(self, ns_int beadid_i, ns_int beadid_j, real distance2) nogil: | ||
| # Important: If this function returns 0, it means that memory allocation failed |
There was a problem hiding this comment.
I'm OK with this now. But a return value should either be a useful variable or an exit code. Not both. The better option is to later query self.npairs or to update on out pointer.
| continue | ||
|
|
||
|
|
||
| for zi in range(DIM): |
There was a problem hiding this comment.
this loop is missing a memory error check
| current_beadid = search_ids_view[i] | ||
|
|
||
| cellindex = self.grid.cellids[current_beadid] | ||
| self.grid.cellid2cellxyz(cellindex, cellxyz) |
|
|
||
| self.coords_bbox = self.box.fast_put_atoms_in_bbox(coords) | ||
|
|
||
| if cutoff < 0: |
There was a problem hiding this comment.
I prefer if these failures are at the beginning. This avoids unnecessary work when we know that the computation can't be done anyway.
|
|
||
| return 1 | ||
|
|
||
| def get_pairs(self): |
There was a problem hiding this comment.
this function can get stale
results.add_neighbors(...)
results.get_pairs(...)
results.add_neighbors(...)
# will not update pairs buffer values
results.get_pairs(...) This is true for all get_* functions. The solution would be that add_neighbor is invalidating the buffers.
I know this is not used incorrectly at the moment and cannot be once the result is passed back into the python stack. For future proofing it would still be nice if get_pairs will always works. Who knows what we use this for in the future.
Codecov Report
@@ Coverage Diff @@
## develop #1996 +/- ##
========================================
Coverage 88.48% 88.48%
========================================
Files 142 142
Lines 17203 17203
Branches 2635 2635
========================================
Hits 15222 15222
Misses 1385 1385
Partials 596 596
Continue to review full report at Codecov.
|
|
@ayushsuhane you currently have a good overview of the cell list algorithms. Could you also give this a careful read and say if it matches your needs. |
|
|
||
| for j in range(self.grid.nbeads[cellindex_adjacent]): | ||
| bid = self.grid.beadids[cellindex_adjacent * self.grid.nbeads_per_cell + j] | ||
| if checked[bid] != 0: |
There was a problem hiding this comment.
I think it could be better to have the checked flag on the cell rather than every coordinate in a cell. Since every cell is an entity as opposed to ever atom in cell-lists. AFAIK, every time we need to do the evaluation, we need to check atleast one cell (all coordinates in the cell).
Here, cell_index probe might become useful to reduce the the computations since the cutoff radius is always smaller than or equal to the cellsize. To prove my point, consider a box with size 10, cellsize 2 in all direction. For a cutoff distance of 0.4 and at the centre of any cell, all the shifts with the probe will still land the coordinates in the same cell, and if the checkcell is present, we just need to compute the distances once and not check every coordinate individually.
|
|
||
| # Get the cell index corresponding to the coord | ||
| cellindex_adjacent = self.grid.coord2cellid(shifted_coord) | ||
| cellindex_probe = self.grid.coord2cellid(probe) |
There was a problem hiding this comment.
I am not sure if I understand the use of probe here?
| probe[ZZ] = self.coords[current_beadid, ZZ] + (zi - 1) * self.cutoff | ||
| # Make sure the shifted coordinates is inside the brick-shaped box | ||
| for m in range(DIM - 1, -1, -1): | ||
| while shifted_coord[m] < 0: |
There was a problem hiding this comment.
I think it would be better if order can be used here in any way. It would be helpful for searching all the pairs such as in guess_bonds i.e we just need to evaluate half the pairs i.e. 27/2 cells to get all the pairs (as every mirror pair will be included exactly once in those half-searches).
|
|
||
| self.prepared = True | ||
|
|
||
| def search(self, search_ids=None): |
There was a problem hiding this comment.
I would have preferred if it takes input coordinates i.e. u.atoms.positions as argument, but this also works (specially considering its use with universe instances.
|
|
||
| # Allocate memory | ||
| self.beadids = <ns_int *> PyMem_Malloc(sizeof(ns_int) * self.size * self.nbeads_per_cell) #np.empty((self.size, nbeads_max), dtype=np.int) | ||
| if not self.beadids: |
There was a problem hiding this comment.
the name nbeads_per_cell is little confusing, its the maximum number of beads in any cell right?
| return <ns_int> (coord[ZZ] / self.cellsize[ZZ]) * (self.ncells[XX] * self.ncells[YY]) +\ | ||
| <ns_int> (coord[YY] / self.cellsize[YY]) * self.ncells[XX] + \ | ||
| <ns_int> (coord[XX] / self.cellsize[XX]) | ||
|
|
There was a problem hiding this comment.
Can use pre-calculated cell_offsets, as it is a highly used function.
| if (box[XX, XX] == 0) or (box[YY, YY] == 0) or (box[ZZ, ZZ] == 0): | ||
| raise ValueError("Box does not correspond to PBC=xyz") | ||
| self.fast_update(box) | ||
|
|
There was a problem hiding this comment.
I think a functionality to handle no pbc can be implemented here. Probably something like
def fast_update_nopbc(box):
min_dim = coords.min(axis=0)
max_dim = coords.max(axis=0)
box[:3] = max_dim.max(axis=0) - min_dim.min(axis=0)
box[3:] = 90. # making it an orthogonal box
and handle the same in search function.
|
@ayushsuhane can you take the docs (& any other changes made since you forked from this) and add it to your PR |
|
Yes, All of the desired changes are already there. I will recheck once again though. |
Changes made in this Pull Request:
PR Checklist