In [1]:

    

import os, os.path
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
from scipy.special import sph_harm
import weave
from weave import converters
from numba import jit, vectorize
import numexpr
%load_ext Cython


    

In [2]:

    

from colloids import periodic


    

In [ ]:

    




    

In [4]:

    

periodic.periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,3,3])


    

    
Out[4]:





array([[ 1.,  1.,  1.]])

In [5]:

    

periodic.periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,1,2])


    

    
Out[5]:





array([[ 1.,  0.,  1.]])

In [6]:

    

periodic.periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [0.9,0.5,0.5])


    

    
Out[6]:





array([[ 0.1,  0. ,  0. ]])

In [10]:

    

reload(periodic)


    

    
Out[10]:





<module 'colloids.periodic' from '/home/mathieu/src/colloids/python/colloids/periodic.py'>

In [3]:

    

u = np.zeros([1000000, 3])
v = 3*np.random.rand(1000000, 3)
u.shape, v.shape


    

    
Out[3]:





((1000000, 3), (1000000, 3))

In [8]:

    

%timeit periodic.periodify(u, v, [0.9,0.5,0.5])


    

    




1 loop, best of 3: 388 ms per loop

In [9]:

    

def python_periodify(u, v, periods=None):
    assert u.shape == v.shape
    diff = np.array(v, float) - u
    if periods is None:
        return diff
    assert len(periods) == u.shape[-1]
    #ensures the largest coordinate is smaller than half a period
    half = 0.5*np.array(periods)
    for i in range(diff.shape[0]):
        for j in range(diff.shape[1]):
            while diff[i,j] > half[j]:
                diff[i,j] -= periods[j]
            while diff[i,j] < -half[j]:
                diff[i,j] += periods[j]
    return diff


    

In [10]:

    

python_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]))


    

    
Out[10]:





array([[ 1.,  1.,  1.]])

In [11]:

    

python_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,3,3])


    

    
Out[11]:





array([[ 1.,  1.,  1.]])

In [12]:

    

python_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,1,2])


    

    
Out[12]:





array([[ 1.,  0.,  1.]])

In [13]:

    

python_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [0.9,0.5,0.5])


    

    
Out[13]:





array([[ 0.1,  0. ,  0. ]])

In [14]:

    

%timeit python_periodify(u, v, [0.9,0.5,0.5])


    

    




1 loop, best of 3: 8.17 s per loop

In [26]:

    

def numpy_periodify(u, v, periods=None):
    """Given two arrays of points in a d-dimentional space with periodic boundary conditions, find the shortest vector between each pair"""
    assert u.shape == v.shape
    diff = np.array(v, float) - u
    if periods is None:
        return diff
    assert len(periods) == u.shape[-1]
    box_r = 1/np.array(periods)
    return diff - periods * np.floor(diff * box_r + 0.5)
    k = diff * box_r  + np.where(diff >= 0.5, 0.5, -0.5)
    diff -= k.astype(int) * periods
    return diff


    

In [27]:

    

%timeit numpy_periodify(u, v, [0.9,0.5,0.5])


    

    




10 loops, best of 3: 90.3 ms per loop

In [26]:

    

from numba import guvectorize
from math import floor
#@jit(["float64[:](float64[:], float64[:], float64)"], nopython=True)
#@guvectorize(["(float64[:], float64[:], float64[:], float64[:])"], '(n),(n),()->(n)')
def numba_periodify(u, v, period, diff):
    """Given two arrays of points in a d-dimentional space with periodic boundary conditions, find the shortest vector between each pair"""
    #assert u.shape == v.shape
    for i in range(u.shape[0]):
        diff[i] = v[i] - u[i]
    if period[0] > 0:
        for i in range(u.shape[0]):
            diff[i] -= period[0] * floor(diff[i] /period[0] + 0.5)
     #   diff -= period[0] * np.floor(diff /period[0] + 0.5)
      #  print('a')
    #k = diff * box_r  + np.where(diff >= 0.5, 0.5, -0.5)
    #diff -= k.astype(int) * periods
    #return diff


    

In [27]:

    

numba_periodify(np.array([0,0,0],float), np.array([1,1,1], float), 0.)


    

    
Out[27]:





array([ 1.,  1.,  1.])

In [28]:

    

numba_periodify(np.array([0,0,0], float), np.array([1,1,1], float), 3)


    

    
Out[28]:





array([ 1.,  1.,  1.])

In [29]:

    

numba_periodify(np.array([0,0,0], float), np.array([1,1,1], float), 1)


    

    
Out[29]:





array([ 0.,  0.,  0.])

In [30]:

    

numba_periodify(np.array([0,0,0], float), np.array([1,1,1], float), 0.9)


    

    
Out[30]:





array([ 0.1,  0.1,  0.1])

In [35]:

    

numba_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), np.transpose([[0.9,0.5,0.5]]))


    

    
Out[35]:





array([[[ 0.1,  0.1,  0.1]],

       [[ 0. ,  0. ,  0. ]],

       [[ 0. ,  0. ,  0. ]]])

In [31]:

    

%timeit numba_periodify(u.ravel(), v.ravel(), 0.9)


    

    




10 loops, best of 3: 24.6 ms per loop

In [3]:

    

%%cython --annotate
import numpy as np
cimport numpy as np
cimport cython

@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function
@cython.cdivision(True) # turn off division by zero checking for entire function
def cython_periodify(
    np.ndarray[np.float64_t, ndim=2] u, 
    np.ndarray[np.float64_t, ndim=2] v, 
    periods=None
):
    assert u.shape[0] == v.shape[0]
    assert u.shape[1] == v.shape[1]
    #cdef np.ndarray[np.float64_t, ndim=2] diff = v - u
    if periods is None:
        return v - u
    assert len(periods) == u.shape[1]
    cdef np.ndarray[np.float64_t, ndim=2] diff = np.zeros_like(u)
    #ensures the largest coordinate is smaller than half a period
    cdef np.ndarray[np.float64_t, ndim=1] pers = np.array(periods)
    cdef np.ndarray[np.float64_t, ndim=1] half = 0.5*np.array(periods)
    cdef int i,j
    cdef float d,h,p
    for i in range(diff.shape[0]):
        for j in range(diff.shape[1]):
            d = v[i,j] - u[i,j]
            h = half[j]
            p = pers[j]
            while d > h:
                d -= p
            while d < -h:
                d += p
            diff[i,j] = d
    return diff


    

    




In file included from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1809:0,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /data/mleocmach/.cache/ipython/cython/_cython_magic_70908a5941237c0652b8ec03b5e71a99.c:257:
/home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
In file included from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:27:0,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /data/mleocmach/.cache/ipython/cython/_cython_magic_70908a5941237c0652b8ec03b5e71a99.c:257:
/home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/__multiarray_api.h:1453:1: warning: ‘_import_array’ defined but not used [-Wunused-function]
 _import_array(void)
 ^






    
Out[3]:









    
    
    
    


Generated by Cython 0.23.4

    Yellow lines hint at Python interaction.

    Click on a line that starts with a "+" to see the C code that Cython generated for it.

+01: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 02: cimport numpy as np
 03: cimport cython
 04: 
 05: @cython.boundscheck(False) # turn off bounds-checking for entire function
 06: @cython.wraparound(False)  # turn off negative index wrapping for entire function
 07: @cython.cdivision(True) # turn off division by zero checking for entire function
+08: def cython_periodify(
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_1cython_periodify(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_1cython_periodify = {"cython_periodify", (PyCFunction)__pyx_pw_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_1cython_periodify, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_1cython_periodify(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  PyArrayObject *__pyx_v_u = 0;
  PyArrayObject *__pyx_v_v = 0;
  PyObject *__pyx_v_periods = 0;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_periodify (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_u,&__pyx_n_s_v,&__pyx_n_s_periods,0};
    PyObject* values[3] = {0,0,0};
/* … */
  /* function exit code */
  goto __pyx_L0;
  __pyx_L1_error:;
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_cython_periodify(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_u, PyArrayObject *__pyx_v_v, PyObject *__pyx_v_periods) {
  PyArrayObject *__pyx_v_diff = 0;
  PyArrayObject *__pyx_v_pers = 0;
  PyArrayObject *__pyx_v_half = 0;
  int __pyx_v_i;
  int __pyx_v_j;
  float __pyx_v_d;
  float __pyx_v_h;
  float __pyx_v_p;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_diff;
  __Pyx_Buffer __pyx_pybuffer_diff;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_half;
  __Pyx_Buffer __pyx_pybuffer_half;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_pers;
  __Pyx_Buffer __pyx_pybuffer_pers;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_u;
  __Pyx_Buffer __pyx_pybuffer_u;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_v;
  __Pyx_Buffer __pyx_pybuffer_v;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_periodify", 0);
  __pyx_pybuffer_diff.pybuffer.buf = NULL;
  __pyx_pybuffer_diff.refcount = 0;
  __pyx_pybuffernd_diff.data = NULL;
  __pyx_pybuffernd_diff.rcbuffer = &__pyx_pybuffer_diff;
  __pyx_pybuffer_pers.pybuffer.buf = NULL;
  __pyx_pybuffer_pers.refcount = 0;
  __pyx_pybuffernd_pers.data = NULL;
  __pyx_pybuffernd_pers.rcbuffer = &__pyx_pybuffer_pers;
  __pyx_pybuffer_half.pybuffer.buf = NULL;
  __pyx_pybuffer_half.refcount = 0;
  __pyx_pybuffernd_half.data = NULL;
  __pyx_pybuffernd_half.rcbuffer = &__pyx_pybuffer_half;
  __pyx_pybuffer_u.pybuffer.buf = NULL;
  __pyx_pybuffer_u.refcount = 0;
  __pyx_pybuffernd_u.data = NULL;
  __pyx_pybuffernd_u.rcbuffer = &__pyx_pybuffer_u;
  __pyx_pybuffer_v.pybuffer.buf = NULL;
  __pyx_pybuffer_v.refcount = 0;
  __pyx_pybuffernd_v.data = NULL;
  __pyx_pybuffernd_v.rcbuffer = &__pyx_pybuffer_v;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_u.rcbuffer->pybuffer, (PyObject*)__pyx_v_u, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  }
  __pyx_pybuffernd_u.diminfo[0].strides = __pyx_pybuffernd_u.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_u.diminfo[0].shape = __pyx_pybuffernd_u.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_u.diminfo[1].strides = __pyx_pybuffernd_u.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_u.diminfo[1].shape = __pyx_pybuffernd_u.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_v.rcbuffer->pybuffer, (PyObject*)__pyx_v_v, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  }
  __pyx_pybuffernd_v.diminfo[0].strides = __pyx_pybuffernd_v.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_v.diminfo[0].shape = __pyx_pybuffernd_v.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_v.diminfo[1].strides = __pyx_pybuffernd_v.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_v.diminfo[1].shape = __pyx_pybuffernd_v.rcbuffer->pybuffer.shape[1];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_6);
  __Pyx_XDECREF(__pyx_t_7);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_diff.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_half.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pers.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_u.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_v.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_70908a5941237c0652b8ec03b5e71a99.cython_periodify", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_diff.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_half.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pers.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_u.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_v.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_diff);
  __Pyx_XDECREF((PyObject *)__pyx_v_pers);
  __Pyx_XDECREF((PyObject *)__pyx_v_half);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__7 = PyTuple_Pack(11, __pyx_n_s_u, __pyx_n_s_v, __pyx_n_s_periods, __pyx_n_s_diff, __pyx_n_s_pers, __pyx_n_s_half, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_d, __pyx_n_s_h, __pyx_n_s_p); if (unlikely(!__pyx_tuple__7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_tuple__7);
  __Pyx_GIVEREF(__pyx_tuple__7);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_1cython_periodify, NULL, __pyx_n_s_cython_magic_70908a5941237c0652); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_cython_periodify, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 09:     np.ndarray[np.float64_t, ndim=2] u,
 10:     np.ndarray[np.float64_t, ndim=2] v,
+11:     periods=None
    values[2] = ((PyObject *)Py_None);
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_u)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        case  1:
        if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_v)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("cython_periodify", 0, 2, 3, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
        case  2:
        if (kw_args > 0) {
          PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_periods);
          if (value) { values[2] = value; kw_args--; }
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "cython_periodify") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
      }
    } else {
      switch (PyTuple_GET_SIZE(__pyx_args)) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_u = ((PyArrayObject *)values[0]);
    __pyx_v_v = ((PyArrayObject *)values[1]);
    __pyx_v_periods = values[2];
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("cython_periodify", 0, 2, 3, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_70908a5941237c0652b8ec03b5e71a99.cython_periodify", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_u), __pyx_ptype_5numpy_ndarray, 1, "u", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_v), __pyx_ptype_5numpy_ndarray, 1, "v", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_r = __pyx_pf_46_cython_magic_70908a5941237c0652b8ec03b5e71a99_cython_periodify(__pyx_self, __pyx_v_u, __pyx_v_v, __pyx_v_periods);
 12: ):
+13:     assert u.shape[0] == v.shape[0]
  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(!Py_OptimizeFlag)) {
    if (unlikely(!(((__pyx_v_u->dimensions[0]) == (__pyx_v_v->dimensions[0])) != 0))) {
      PyErr_SetNone(PyExc_AssertionError);
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    }
  }
  #endif
+14:     assert u.shape[1] == v.shape[1]
  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(!Py_OptimizeFlag)) {
    if (unlikely(!(((__pyx_v_u->dimensions[1]) == (__pyx_v_v->dimensions[1])) != 0))) {
      PyErr_SetNone(PyExc_AssertionError);
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    }
  }
  #endif
 15:     #cdef np.ndarray[np.float64_t, ndim=2] diff = v - u
+16:     if periods is None:
  __pyx_t_1 = (__pyx_v_periods == Py_None);
  __pyx_t_2 = (__pyx_t_1 != 0);
  if (__pyx_t_2) {
/* … */
  }
+17:         return v - u
    __Pyx_XDECREF(__pyx_r);
    __pyx_t_3 = PyNumber_Subtract(((PyObject *)__pyx_v_v), ((PyObject *)__pyx_v_u)); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_r = __pyx_t_3;
    __pyx_t_3 = 0;
    goto __pyx_L0;
+18:     assert len(periods) == u.shape[1]
  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(!Py_OptimizeFlag)) {
    __pyx_t_4 = PyObject_Length(__pyx_v_periods); if (unlikely(__pyx_t_4 == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    if (unlikely(!((__pyx_t_4 == (__pyx_v_u->dimensions[1])) != 0))) {
      PyErr_SetNone(PyExc_AssertionError);
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    }
  }
  #endif
+19:     cdef np.ndarray[np.float64_t, ndim=2] diff = np.zeros_like(u)
  __pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_zeros_like); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_6))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_6);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_6, function);
    }
  }
  if (!__pyx_t_5) {
    __pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_6, ((PyObject *)__pyx_v_u)); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
  } else {
    __pyx_t_7 = PyTuple_New(1+1); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_7);
    __Pyx_GIVEREF(__pyx_t_5); PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_5); __pyx_t_5 = NULL;
    __Pyx_INCREF(((PyObject *)__pyx_v_u));
    __Pyx_GIVEREF(((PyObject *)__pyx_v_u));
    PyTuple_SET_ITEM(__pyx_t_7, 0+1, ((PyObject *)__pyx_v_u));
    __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_7, NULL); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  }
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_t_8 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_diff.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_diff = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_diff.rcbuffer->pybuffer.buf = NULL;
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    } else {__pyx_pybuffernd_diff.diminfo[0].strides = __pyx_pybuffernd_diff.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_diff.diminfo[0].shape = __pyx_pybuffernd_diff.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_diff.diminfo[1].strides = __pyx_pybuffernd_diff.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_diff.diminfo[1].shape = __pyx_pybuffernd_diff.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_8 = 0;
  __pyx_v_diff = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;
 20:     #ensures the largest coordinate is smaller than half a period
+21:     cdef np.ndarray[np.float64_t, ndim=1] pers = np.array(periods)
  __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_array); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_7))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_7);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_7);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_7, function);
    }
  }
  if (!__pyx_t_6) {
    __pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_7, __pyx_v_periods); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
  } else {
    __pyx_t_5 = PyTuple_New(1+1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_GIVEREF(__pyx_t_6); PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_6); __pyx_t_6 = NULL;
    __Pyx_INCREF(__pyx_v_periods);
    __Pyx_GIVEREF(__pyx_v_periods);
    PyTuple_SET_ITEM(__pyx_t_5, 0+1, __pyx_v_periods);
    __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_7, __pyx_t_5, NULL); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  }
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_t_9 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_pers.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_pers = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_pers.rcbuffer->pybuffer.buf = NULL;
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    } else {__pyx_pybuffernd_pers.diminfo[0].strides = __pyx_pybuffernd_pers.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_pers.diminfo[0].shape = __pyx_pybuffernd_pers.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_9 = 0;
  __pyx_v_pers = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;
+22:     cdef np.ndarray[np.float64_t, ndim=1] half = 0.5*np.array(periods)
  __pyx_t_7 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_7);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_7, __pyx_n_s_array); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  __pyx_t_7 = NULL;
  if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_5))) {
    __pyx_t_7 = PyMethod_GET_SELF(__pyx_t_5);
    if (likely(__pyx_t_7)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_7);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_5, function);
    }
  }
  if (!__pyx_t_7) {
    __pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_v_periods); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
  } else {
    __pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_6);
    __Pyx_GIVEREF(__pyx_t_7); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_7); __pyx_t_7 = NULL;
    __Pyx_INCREF(__pyx_v_periods);
    __Pyx_GIVEREF(__pyx_v_periods);
    PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_v_periods);
    __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_6, NULL); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  }
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = PyNumber_Multiply(__pyx_float_0_5, __pyx_t_3); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (!(likely(((__pyx_t_5) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_5, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_t_10 = ((PyArrayObject *)__pyx_t_5);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_half.rcbuffer->pybuffer, (PyObject*)__pyx_t_10, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_half = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_half.rcbuffer->pybuffer.buf = NULL;
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 22; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    } else {__pyx_pybuffernd_half.diminfo[0].strides = __pyx_pybuffernd_half.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_half.diminfo[0].shape = __pyx_pybuffernd_half.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_10 = 0;
  __pyx_v_half = ((PyArrayObject *)__pyx_t_5);
  __pyx_t_5 = 0;
 23:     cdef int i,j
 24:     cdef float d,h,p
+25:     for i in range(diff.shape[0]):
  __pyx_t_11 = (__pyx_v_diff->dimensions[0]);
  for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
    __pyx_v_i = __pyx_t_12;
+26:         for j in range(diff.shape[1]):
    __pyx_t_13 = (__pyx_v_diff->dimensions[1]);
    for (__pyx_t_14 = 0; __pyx_t_14 < __pyx_t_13; __pyx_t_14+=1) {
      __pyx_v_j = __pyx_t_14;
+27:             d = v[i,j] - u[i,j]
      __pyx_t_15 = __pyx_v_i;
      __pyx_t_16 = __pyx_v_j;
      __pyx_t_17 = __pyx_v_i;
      __pyx_t_18 = __pyx_v_j;
      __pyx_v_d = ((*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_v.rcbuffer->pybuffer.buf, __pyx_t_15, __pyx_pybuffernd_v.diminfo[0].strides, __pyx_t_16, __pyx_pybuffernd_v.diminfo[1].strides)) - (*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_u.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_u.diminfo[0].strides, __pyx_t_18, __pyx_pybuffernd_u.diminfo[1].strides)));
+28:             h = half[j]
      __pyx_t_19 = __pyx_v_j;
      __pyx_v_h = (*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_half.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_half.diminfo[0].strides));
+29:             p = pers[j]
      __pyx_t_20 = __pyx_v_j;
      __pyx_v_p = (*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_pers.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_pers.diminfo[0].strides));
+30:             while d > h:
      while (1) {
        __pyx_t_2 = ((__pyx_v_d > __pyx_v_h) != 0);
        if (!__pyx_t_2) break;
+31:                 d -= p
        __pyx_v_d = (__pyx_v_d - __pyx_v_p);
      }
+32:             while d < -h:
      while (1) {
        __pyx_t_2 = ((__pyx_v_d < (-__pyx_v_h)) != 0);
        if (!__pyx_t_2) break;
+33:                 d += p
        __pyx_v_d = (__pyx_v_d + __pyx_v_p);
      }
+34:             diff[i,j] = d
      __pyx_t_21 = __pyx_v_i;
      __pyx_t_22 = __pyx_v_j;
      *__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_diff.rcbuffer->pybuffer.buf, __pyx_t_21, __pyx_pybuffernd_diff.diminfo[0].strides, __pyx_t_22, __pyx_pybuffernd_diff.diminfo[1].strides) = __pyx_v_d;
    }
  }
+35:     return diff
  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_diff));
  __pyx_r = ((PyObject *)__pyx_v_diff);
  goto __pyx_L0;

In [54]:

    

%timeit cython_periodify(u, v, np.array([0.9,0.5,0.5]))


    

    




10 loops, best of 3: 55.4 ms per loop

In [47]:

    

from numba import vectorize, float64
from math import floor
@vectorize([float64(float64,float64,float64)], nopython=True)
def v_periodify(u,v,period=-1.0):
    diff = v - u
    if period <= 0:
        return diff
    return diff - period * floor(diff /period + 0.5)


    

In [48]:

    

v_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), -1)


    

    
Out[48]:





array([[ 1.,  1.,  1.]])

In [49]:

    

v_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,3,3])


    

    
Out[49]:





array([[ 1.,  1.,  1.]])

In [50]:

    

v_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [3,1,2])


    

    
Out[50]:





array([[ 1.,  0., -1.]])

In [51]:

    

v_periodify(np.array([[0,0,0]]), np.array([[1,1,1]]), [0.9,0.5,0.5])


    

    
Out[51]:





array([[ 0.1,  0. ,  0. ]])

In [52]:

    

%timeit v_periodify(u, v, 0.9)


    

    




100 loops, best of 3: 15.8 ms per loop

In [55]:

    

v_periodify(u, v, 0.9).shape


    

    
Out[55]:





(1000000, 3)

In [56]:

    

v_periodify(u, v, [0.9,0.5,0.5]).shape


    

    
Out[56]:





(1000000, 3)

In [3]:

    

pos = np.loadtxt('../../../multiscale/test_input/3d_6_0.54_0.dat')


    

In [22]:

    

pos.min(0), pos.max(0)


    

    
Out[22]:





(array([ 0.0167643 ,  0.017033  ,  0.00916006,  3.22366   ]),
 array([ 202.999  ,  202.999  ,  202.985  ,    4.85154]))

In [8]:

    

g = periodic.get_rdf(0, pos[:,:-1], 100, 203.0)


    

In [9]:

    

plt.plot(g)


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x7f433b776fd0>]






    

In [13]:

    

g = periodic.get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

In [20]:

    

plt.plot(g.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)))


    

    
Out[20]:





[<matplotlib.lines.Line2D at 0x7f4339f37150>]






    

In [70]:

    

%timeit periodic.get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

    




1 loop, best of 3: 2.87 s per loop

In [34]:

    

def get_mono_rdf(pos, Nbins, L):
    maxdist = L/2.0
    g = np.zeros([len(pos), Nbins], int)
    #conversion factor between indices and bins
    l2r = Nbins/maxdist
    for i in range(len(pos)):
        good = np.ones(len(pos), bool)
        good[i] = False
        dists = np.sqrt((periodic.periodify(pos[i], pos[good], L)**2).sum(-1))
        rs = (dists * l2r).astype(int)
        np.add.at(g[i], rs[rs<Nbins], 1)
    return g


    

In [56]:

    

g2 = get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

In [32]:

    

plt.plot(g.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)))
plt.plot(g2.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)), ls=':')


    

    
Out[32]:





[<matplotlib.lines.Line2D at 0x7f4339d88dd0>]






    

In [81]:

    

from numba import float64, int64
from math import sqrt, floor
@jit(int64[:,:](float64[:,:], int64, float64), nopython=True)
def numba_get_mono_rdf(pos, Nbins, L):
    maxdist = L/2.0
    rL = 1.0/L
    g = np.zeros((pos.shape[0], Nbins), np.int64)
    #conversion factor between indices and bins
    l2r = Nbins/maxdist
    for i in range(pos.shape[0]):
        for j in range(pos.shape[0]):
            if j==i:
                continue
            dist = 0
            for d in range(pos.shape[1]):
                diff = pos[i,d] - pos[j,d]
                diff -= L * floor(diff * rL + 0.5)
                dist += diff*diff
            dist = sqrt(dist)
            r = int(dist * l2r)
            if r<Nbins:
                g[i,r] += 1
    return g


    

In [82]:

    

g3 = numba_get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

In [83]:

    

plt.plot(g2.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)))
plt.plot(g3.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)), ls=':')


    

    
Out[83]:





[<matplotlib.lines.Line2D at 0x7f43390c1d90>]






    

In [84]:

    

%timeit numba_get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

    




1 loop, best of 3: 5.03 s per loop

In [4]:

    

%%cython --annotate
import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport sqrt, floor

@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function
@cython.cdivision(True) # turn off division by zero checking for entire function
def cython_get_mono_rdf(
    np.ndarray[np.float64_t, ndim=2] pos, 
    int Nbins, 
    float L
):
    cdef float maxdist = L/2.0
    cdef np.ndarray[np.int64_t, ndim=2] g = np.zeros((pos.shape[0], Nbins), np.int64)
    #conversion factor between indices and bins
    cdef float l2r = Nbins/maxdist
    cdef int i,j, r, d
    cdef float dist, diff
    for i in range(pos.shape[0]):
        for j in range(pos.shape[0]):
            if j==i:
                continue
            dist = 0
            for d in range(pos.shape[1]):
                diff = pos[i,d] - pos[j,d]
                diff -= L * floor(diff /L + 0.5)
                dist += diff*diff
            dist = sqrt(dist)
            r = int(dist * l2r)
            if r<Nbins:
                g[i,r] += 1
    return g


    

    




In file included from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1809:0,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /data/mleocmach/.cache/ipython/cython/_cython_magic_eab24a72aafa0e638dac17a47d127bbe.c:258:
/home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^
In file included from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:27:0,
                 from /home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /data/mleocmach/.cache/ipython/cython/_cython_magic_eab24a72aafa0e638dac17a47d127bbe.c:258:
/home/mathieu/anaconda3/envs/python2/lib/python2.7/site-packages/numpy/core/include/numpy/__multiarray_api.h:1453:1: warning: ‘_import_array’ defined but not used [-Wunused-function]
 _import_array(void)
 ^






    
Out[4]:









    
    
    
    


Generated by Cython 0.23.4

    Yellow lines hint at Python interaction.

    Click on a line that starts with a "+" to see the C code that Cython generated for it.

+01: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 02: cimport numpy as np
 03: cimport cython
 04: from libc.math cimport sqrt, floor
 05: 
 06: @cython.boundscheck(False) # turn off bounds-checking for entire function
 07: @cython.wraparound(False)  # turn off negative index wrapping for entire function
 08: @cython.cdivision(True) # turn off division by zero checking for entire function
+09: def cython_get_mono_rdf(
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_1cython_get_mono_rdf(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_1cython_get_mono_rdf = {"cython_get_mono_rdf", (PyCFunction)__pyx_pw_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_1cython_get_mono_rdf, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_1cython_get_mono_rdf(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  PyArrayObject *__pyx_v_pos = 0;
  int __pyx_v_Nbins;
  float __pyx_v_L;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_get_mono_rdf (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_pos,&__pyx_n_s_Nbins,&__pyx_n_s_L,0};
    PyObject* values[3] = {0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_pos)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        case  1:
        if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_Nbins)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("cython_get_mono_rdf", 1, 3, 3, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
        case  2:
        if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_L)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("cython_get_mono_rdf", 1, 3, 3, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "cython_get_mono_rdf") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
      }
    } else if (PyTuple_GET_SIZE(__pyx_args) != 3) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
    }
    __pyx_v_pos = ((PyArrayObject *)values[0]);
    __pyx_v_Nbins = __Pyx_PyInt_As_int(values[1]); if (unlikely((__pyx_v_Nbins == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
    __pyx_v_L = __pyx_PyFloat_AsFloat(values[2]); if (unlikely((__pyx_v_L == (float)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("cython_get_mono_rdf", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_eab24a72aafa0e638dac17a47d127bbe.cython_get_mono_rdf", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_pos), __pyx_ptype_5numpy_ndarray, 1, "pos", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 10; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_r = __pyx_pf_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_cython_get_mono_rdf(__pyx_self, __pyx_v_pos, __pyx_v_Nbins, __pyx_v_L);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  goto __pyx_L0;
  __pyx_L1_error:;
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_cython_get_mono_rdf(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_pos, int __pyx_v_Nbins, float __pyx_v_L) {
  float __pyx_v_maxdist;
  PyArrayObject *__pyx_v_g = 0;
  float __pyx_v_l2r;
  int __pyx_v_i;
  int __pyx_v_j;
  int __pyx_v_r;
  int __pyx_v_d;
  float __pyx_v_dist;
  float __pyx_v_diff;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_pos;
  __Pyx_Buffer __pyx_pybuffer_pos;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_get_mono_rdf", 0);
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  __pyx_pybuffer_pos.pybuffer.buf = NULL;
  __pyx_pybuffer_pos.refcount = 0;
  __pyx_pybuffernd_pos.data = NULL;
  __pyx_pybuffernd_pos.rcbuffer = &__pyx_pybuffer_pos;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_pos.rcbuffer->pybuffer, (PyObject*)__pyx_v_pos, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  }
  __pyx_pybuffernd_pos.diminfo[0].strides = __pyx_pybuffernd_pos.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_pos.diminfo[0].shape = __pyx_pybuffernd_pos.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_pos.diminfo[1].strides = __pyx_pybuffernd_pos.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_pos.diminfo[1].shape = __pyx_pybuffernd_pos.rcbuffer->pybuffer.shape[1];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_7);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pos.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_eab24a72aafa0e638dac17a47d127bbe.cython_get_mono_rdf", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pos.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_g);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__7 = PyTuple_Pack(12, __pyx_n_s_pos, __pyx_n_s_Nbins, __pyx_n_s_L, __pyx_n_s_maxdist, __pyx_n_s_g, __pyx_n_s_l2r, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_r, __pyx_n_s_d, __pyx_n_s_dist, __pyx_n_s_diff); if (unlikely(!__pyx_tuple__7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_tuple__7);
  __Pyx_GIVEREF(__pyx_tuple__7);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_eab24a72aafa0e638dac17a47d127bbe_1cython_get_mono_rdf, NULL, __pyx_n_s_cython_magic_eab24a72aafa0e638d); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_cython_get_mono_rdf, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 10:     np.ndarray[np.float64_t, ndim=2] pos,
 11:     int Nbins,
 12:     float L
 13: ):
+14:     cdef float maxdist = L/2.0
  __pyx_v_maxdist = (__pyx_v_L / 2.0);
+15:     cdef np.ndarray[np.int64_t, ndim=2] g = np.zeros((pos.shape[0], Nbins), np.int64)
  __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyInt_From_Py_intptr_t((__pyx_v_pos->dimensions[0])); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_Nbins); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_2);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_2);
  __Pyx_GIVEREF(__pyx_t_4);
  PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_4);
  __pyx_t_2 = 0;
  __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_int64); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = NULL;
  __pyx_t_6 = 0;
  if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_6 = 1;
    }
  }
  __pyx_t_7 = PyTuple_New(2+__pyx_t_6); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_7);
  if (__pyx_t_4) {
    __Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_4); __pyx_t_4 = NULL;
  }
  __Pyx_GIVEREF(__pyx_t_5);
  PyTuple_SET_ITEM(__pyx_t_7, 0+__pyx_t_6, __pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_2);
  PyTuple_SET_ITEM(__pyx_t_7, 1+__pyx_t_6, __pyx_t_2);
  __pyx_t_5 = 0;
  __pyx_t_2 = 0;
  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_7, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_t_8 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_5numpy_int64_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_g = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_g.rcbuffer->pybuffer.buf = NULL;
      {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    } else {__pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_g.diminfo[1].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_g.diminfo[1].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_8 = 0;
  __pyx_v_g = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;
 16:     #conversion factor between indices and bins
+17:     cdef float l2r = Nbins/maxdist
  __pyx_v_l2r = (__pyx_v_Nbins / __pyx_v_maxdist);
 18:     cdef int i,j, r, d
 19:     cdef float dist, diff
+20:     for i in range(pos.shape[0]):
  __pyx_t_9 = (__pyx_v_pos->dimensions[0]);
  for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
    __pyx_v_i = __pyx_t_10;
+21:         for j in range(pos.shape[0]):
    __pyx_t_11 = (__pyx_v_pos->dimensions[0]);
    for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
      __pyx_v_j = __pyx_t_12;
+22:             if j==i:
      __pyx_t_13 = ((__pyx_v_j == __pyx_v_i) != 0);
      if (__pyx_t_13) {
/* … */
      }
+23:                 continue
        goto __pyx_L5_continue;
+24:             dist = 0
      __pyx_v_dist = 0.0;
+25:             for d in range(pos.shape[1]):
      __pyx_t_14 = (__pyx_v_pos->dimensions[1]);
      for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
        __pyx_v_d = __pyx_t_15;
+26:                 diff = pos[i,d] - pos[j,d]
        __pyx_t_16 = __pyx_v_i;
        __pyx_t_17 = __pyx_v_d;
        __pyx_t_18 = __pyx_v_j;
        __pyx_t_19 = __pyx_v_d;
        __pyx_v_diff = ((*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_pos.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_pos.diminfo[0].strides, __pyx_t_17, __pyx_pybuffernd_pos.diminfo[1].strides)) - (*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_pos.rcbuffer->pybuffer.buf, __pyx_t_18, __pyx_pybuffernd_pos.diminfo[0].strides, __pyx_t_19, __pyx_pybuffernd_pos.diminfo[1].strides)));
+27:                 diff -= L * floor(diff /L + 0.5)
        __pyx_v_diff = (__pyx_v_diff - (__pyx_v_L * floor(((__pyx_v_diff / __pyx_v_L) + 0.5))));
+28:                 dist += diff*diff
        __pyx_v_dist = (__pyx_v_dist + (__pyx_v_diff * __pyx_v_diff));
      }
+29:             dist = sqrt(dist)
      __pyx_v_dist = sqrt(__pyx_v_dist);
+30:             r = int(dist * l2r)
      __pyx_v_r = ((int)(__pyx_v_dist * __pyx_v_l2r));
+31:             if r<Nbins:
      __pyx_t_13 = ((__pyx_v_r < __pyx_v_Nbins) != 0);
      if (__pyx_t_13) {
/* … */
      }
      __pyx_L5_continue:;
    }
  }
+32:                 g[i,r] += 1
        __pyx_t_20 = __pyx_v_i;
        __pyx_t_21 = __pyx_v_r;
        *__Pyx_BufPtrStrided2d(__pyx_t_5numpy_int64_t *, __pyx_pybuffernd_g.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_g.diminfo[0].strides, __pyx_t_21, __pyx_pybuffernd_g.diminfo[1].strides) += 1;
+33:     return g
  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_g));
  __pyx_r = ((PyObject *)__pyx_v_g);
  goto __pyx_L0;

In [65]:

    

g4 = cython_get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

In [66]:

    

plt.plot(g2.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)))
plt.plot(g4.mean(0)*203**3/len(pos)/(4/3*np.pi*np.diff((np.arange(1001)/1001.*203./2)**3)), ls=':')


    

    
Out[66]:





[<matplotlib.lines.Line2D at 0x7f4339abf210>]






    

In [69]:

    

%timeit cython_get_mono_rdf(pos[:,:-1], 1000, 203.0)


    

    




1 loop, best of 3: 7.99 s per loop

In [ ]: