enabled selecting between serial and OpenMP distance functions

[orbeckst]

• all functions in lib.distances acquired the mode= keyword;
  • "serial" selects the standard serial C/Cython code from lib._distances;
  • "OpenMP" selects the OpenMP parallelized C/Cython code from
    lib._distances_openmp
• on systems without an OpenMP capable compiler, this will likely just mean
  that the OpenMP versions are actually serial (the OpenMP pragma should
  get ignored) --- need to be confirmed/tested
• tests in test_distances were refactored so that both serial and OpenMP
  versions are tested
• the current implementation is not very pretty (see lib.distances.select())
  but will make it very easy to add additional acceleration schemes such as
  CUDA or numba as long as the call signatures of the calc_* functions remain
  the same.

[orbeckst]

Btw, I am very open for any comments on how to make this nicer, and I am also not very happy with the name mode for the keyword argument as used to determine if running in serial or OpenMP, e.g

d = MDAnalysis.lib.distances.distance_array(g.positions, ow.positions, box=u.dimensions, 
                                            mode="OpenMP")

[kain88-de]

@orbeckst instead of mode you could use a boolean switch parallel=True/False. Or if the openMP code supports a variable number of jobs/CPUs some packages have adopted a jobs argument, stating the number of CPU's used. -1 then means use all available CPU's.

[jbarnoud]

I agree with @kain88-de : mode is a rather vague keyword.

[richardjgowers]

Could we make this _select so people don't try and use directly?

[orbeckst]

I changed it to _run().

[richardjgowers]

Maybe backend instead of mode?

[jbarnoud]

backend looks clearer already.  I would go for parallel.
Le 11 nov. 2015 18:54, Richard Gowers [email protected] a écrit :Maybe backend instead of mode?

—Reply to this email directly or view it on GitHub.

[orbeckst]

On 11 Nov, 2015, at 09:54, Richard Gowers wrote:

Maybe backend instead of mode?

Funnily enough, that's exactly what I just changed it to. It's a bit technical but the best description I can think of. parallel is not good because parallel="serial" sounds really strange.

With backend we can also imagine backend="CUDA", backend="OpenCL", backend="SPIDAL" etc.

[hainm]

I am +1 for backend if you have more than 2 choices.

[orbeckst]

Rebased to include

• mode --> backend
• select() --> _run()

Also removed cimport c_numpy (see mailing list for additional discussion).

[hainm]

I am not sure why duplicate a lot of code here

$ diff distances_openmp.pyx distances.pyx 
19,20c19,20
< Parallel distance calculation library --- :mod:`MDAnalysis.lib._distances_openmp`
< =================================================================================
---
> Distance calculation library --- :mod:`MDAnalysis.lib._distances`
> =================================================================
22,23c22
< 
< Contains OpenMP versions of the contents of "calc_distances.h"
---
> Serial versions of all distance calculations
25a25,28
> cimport cython
> cimport c_numpy
> c_numpy.import_array()
> 
49a53
>     void minimum_image(double *x, float *box, float *inverse_box)
83c87
<                                    <double*>result.data) 
---
>                                    <double*>result.data)
255a260,307
> 
> 
> @cython.boundscheck(False)
> def contact_matrix_no_pbc(coord, sparse_contacts, cutoff):
>     cdef int rows = len(coord)
>     cdef double cutoff2 = cutoff ** 2
>     cdef float[:, ::1] coord_view = coord
> 
>     cdef int i, j
>     cdef double[3] rr;
>     cdef double dist
>     for i in range(rows):
>         sparse_contacts[i, i] = True
>         for j in range(i+1, rows):
>             rr[0] = coord_view[i, 0] - coord_view[j, 0]
>             rr[1] = coord_view[i, 1] - coord_view[j, 1]
>             rr[2] = coord_view[i, 2] - coord_view[j, 2]
>             dist = rr[0]*rr[0] + rr[1]*rr[1] + rr[2]*rr[2]
>             if dist < cutoff2:
>                 sparse_contacts[i, j] = True
>                 sparse_contacts[j, i] = True
> 
> 
> @cython.boundscheck(False)
> def contact_matrix_pbc(coord, sparse_contacts, box, cutoff):
>     cdef int rows = len(coord)
>     cdef double cutoff2 = cutoff ** 2
>     cdef float[:, ::1] coord_view = coord
>     cdef float[::1] box_view = box
>     cdef float[::1] box_inv = 1. / box
> 
>     cdef int i, j
>     cdef double[3] rr;
>     cdef double dist
>     for i in range(rows):
>         sparse_contacts[i, i] = True
>         for j in range(i+1, rows):
>             rr[0] = coord_view[i, 0] - coord_view[j, 0]
>             rr[1] = coord_view[i, 1] - coord_view[j, 1]
>             rr[2] = coord_view[i, 2] - coord_view[j, 2]
> 
>             minimum_image(rr, &box_view[0], &box_inv[0])
> 
>             dist = rr[0]*rr[0] + rr[1]*rr[1] + rr[2]*rr[2]
> 
>             if dist < cutoff2:
>                 sparse_contacts[i, j] = True
>                 sparse_contacts[j, i] = True

[richardjgowers]

This is pretty minor, but can we make the backend selection case insensitive? I think you could just call this 'openmp' and then call .lower() on whatever string is passed.

[orbeckst]

With OpenMP-versions available I can now implement #530 and remove lib.parallel.distances.

(This comment mainly serves to link to #530.)