8. Speed up 2D particles with Cython

8.1. Python code

import random

import math

import matplotlib.pyplot as plt

import numpy as np

import time

import timeit


# main function
def simulate_py(n, time_step, show_plot=False):

    # random particles
    np.random.seed(0)
    
    p = np.random.rand(n, 2)
    
# Reducing 2D to 1D    p = set_second_column_to_zero(p)

    print("P(0):\n", p)

    # time step loop
    for i in range(time_step):
        
#        print("time step", i)
        
        # calculate total force for each particle
        total_force = combined_force(p, n)
        
#        print("total_force \n",total_force)
        
        # calculate displacement for each particle
        x_det = displacement(total_force, delta_t=0.0001)
        
        #print("x_det \n",x_det)

        # update the particles
        p = update_position(p, x_det)
        
        pos = p
        
        #print("pos \n",pos)
        
        # plot particles
        colors = ['red', 'green', 'blue', 'orange'] 
        if show_plot:

            if i % 2 == 0:

                update_plot(pos,colors)

                print("in iteration", i)
            
    # plot finally result
    print("P({}): ".format(time_step), p)
    
    return p



# Reducing 2 dimensions to 1 dimension for model checking
def set_second_column_to_zero(matrix):

    matrix[:, 1] = 0
    
    return matrix


# Calculate the strength of the repulsion
def f(r, c1=1, c2=1):

    # Calculate the force

    mag = c1 / np.linalg.norm(r)

    #print("mag",mag)

    return mag


# Calculate the total force for each particle
def combined_force(p, n):
    
    total_force = np.zeros_like(p)
    
    for i in range(n):
        
        fn_sum = np.zeros(2)
        
        for j in range(n):
            
            if j != i:
                
                r = p[j] - p[i]
                
                #print("r",r)
                #print("sign",np.sign(r))
                
                fn =  -1 * f(r) * np.sign(r) 

                fn_sum += fn 
                
            total_force[i] = fn_sum
            
    return total_force



# Calculate the displacement between two particles
def displacement(total_force, eta=1, delta_t=1):

    displacement = total_force / eta * delta_t

    return displacement



# Update the position of particles
def update_position(p, delta_r, min_x=0, max_x=1):
    
    new_pos = p + delta_r
    
    x_out_of_bounds = np.logical_or(new_pos[:,0] > max_x, new_pos[:,0] < min_x)
    
    y_out_of_bounds = np.logical_or(new_pos[:,1] > max_x, new_pos[:,1] < min_x)
    
    new_pos[x_out_of_bounds, 0] = np.clip(new_pos[x_out_of_bounds, 0], min_x, max_x)
    
    new_pos[y_out_of_bounds, 1] = np.clip(new_pos[y_out_of_bounds, 1], min_x, max_x)
    
    return new_pos


# Plot
def update_plot(pos,color):

    plt.clf()

    xpos = pos[:, 0]
    
    ypos = pos[:, 1]

    N = len(pos)
    N_color = len(color)
    for i in range(N):
        plt.plot(xpos[i], ypos[i], "o", color=color[i % N_color])

    plt.xlim(left=-0.1, right=1.1)
    
    plt.ylim(bottom=-0.1, top=1.1)

    plt.grid()

    plt.draw()

    plt.pause(0.0001)



# Example usage:
# colors = [random.choice(['r', 'g', 'b', 'y', 'm']) for _ in range(n)]
 
    
#p = simulate_py(n=10, time_step=1000, show_plot=True)

# %prun simulate_py(n=10, time_step=1000, show_plot=False)
compute_time_py = timeit.timeit(lambda: simulate_py(n=10, time_step=1000, show_plot=False), number=1)

print("simulate_py execution time:", compute_time_py)

P(0):
 [[0.5488135  0.71518937]
 [0.60276338 0.54488318]
 [0.4236548  0.64589411]
 [0.43758721 0.891773  ]
 [0.96366276 0.38344152]
 [0.79172504 0.52889492]
 [0.56804456 0.92559664]
 [0.07103606 0.0871293 ]
 [0.0202184  0.83261985]
 [0.77815675 0.87001215]]
P(1000):  [[0.41916695 0.66633421]
 [0.7845857  0.29828616]
 [0.         0.36688618]
 [0.17622127 1.        ]
 [1.         0.        ]
 [1.         0.12007927]
 [0.72441852 1.        ]
 [0.         0.        ]
 [0.         0.93898401]
 [1.         1.        ]]
simulate_py execution time: 1.06936289999976

8.2. Cython code

%load_ext cython

%%cython --force -c=/openmp -a
#%%cython -a

import numpy as np

cimport numpy as np

from cpython cimport array

import array

import cython

import random

import math

import matplotlib.pyplot as plt

import time

import timeit

from cython.parallel import prange


# random particles
cdef initial_particles(int n):
    
    np.random.seed(0)
    
    cdef np.ndarray[np.float64_t, ndim=2] p
    p = np.random.rand(n, 2)
    
    return p


cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
cdef simulate_cy(int n, int time_step, show_plot=False):
    
    cdef np.ndarray[np.float64_t, ndim=2] p

    p = initial_particles(n)
    print("initial_P:\n",p)

    cdef int i
    cdef np.ndarray[np.float64_t, ndim=2] total_force, x_det, pos
    
    for i in range(time_step):
        
        total_force = combined_force(p, n)
    
        x_det = displacement(total_force, delta_t=0.0001)
        
        # update the particles
        p = update_position(p, x_det)
        
        pos = p
        
        # plot particles
        colors = ['red', 'green', 'blue', 'orange'] 
        if show_plot:

            if i % 2 == 0:

                update_plot(pos,colors)

    print("P({}): \n".format(time_step), p)
    
    return p


# Calculate the strength of the repulsion
cdef f(np.ndarray[np.float64_t, ndim=1]r, double c1=1, double c2=1):

    # Calculate the force

    cdef double mag
    
    mag = c1 / math.sqrt(square_sum(r)) 
    
    return mag



cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)

#calculate square
cdef square_sum(np.ndarray[np.float64_t, ndim=1] arr):
    cdef Py_ssize_t i, n
    cdef double s = 0
    n = arr.shape[0]
    for i in prange(n, nogil=True, num_threads=1):
        s += arr[i] * arr[i]
#    print(type(s))
    return s


@cython.boundscheck(False)
@cython.wraparound(False)
# Calculate the total force for each particle
cdef combined_force(np.ndarray[np.float64_t, ndim=2] p, int n):
    
    cdef np.ndarray[np.float64_t, ndim=2] total_force
    
    total_force= np.zeros_like(p,dtype=np.double)
    
    cdef Py_ssize_t i,j
    
    cdef np.ndarray[np.float64_t, ndim=1] r, fn_sum, fn
    
    for i in range(n):
        
        fn_sum = np.zeros(2,dtype=np.double)
        
        for j in range(n):
            
            if j != i:
                
                r = p[j] - p[i]
                
                fn =  -1 * f(r) * np.sign(r) 

                fn_sum += fn 
                
            total_force[i] = fn_sum
            
    return total_force



# Calculate the displacement between two particles
cpdef displacement(np.ndarray[np.float64_t, ndim=2] total_force, double eta=1, double delta_t=1):
    cdef np.ndarray[np.float64_t, ndim=2] displacement

    displacement = total_force / eta * delta_t

    return displacement


cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
# Update the position of particles
cdef update_position(np.ndarray[np.float64_t, ndim=2] p, np.ndarray[np.float64_t, ndim=2] delta_r, double min_x=0, double max_x=1):
    
    cdef np.ndarray[np.float64_t, ndim=2] new_pos
    
    new_pos = p + delta_r
    
    x_out_of_bounds = np.logical_or(new_pos[:,0] > max_x, new_pos[:,0] < min_x)
    
    y_out_of_bounds = np.logical_or(new_pos[:,1] > max_x, new_pos[:,1] < min_x)
    

    for i in range(len(new_pos)):
        if x_out_of_bounds[i]:
            new_pos[i, 0] = clip(new_pos[i, 0], min_x, max_x)
        if y_out_of_bounds[i]:
            new_pos[i, 1] = clip(new_pos[i, 1], min_x, max_x)
    
    return new_pos


cdef double clip(a, min_value, max_value):
    return min(max(a, min_value), max_value)


# Plot
cdef update_plot(pos,color):

    plt.clf()

    xpos = pos[:, 0]
    
    ypos = pos[:, 1]

    N = len(pos)
    N_color = len(color)
    for i in range(N):
        plt.plot(xpos[i], ypos[i], "o", color=color[i % N_color])

    plt.xlim(left=-0.1, right=1.1)
    
    plt.ylim(bottom=-0.1, top=1.1)

    plt.grid()

    plt.draw()

    plt.pause(0.0001)

compute_time_cy = timeit.timeit(lambda: simulate_cy(n=10, time_step=1000, show_plot=False), number=1)

print("simulate_cy execution time:", compute_time_cy)

initial_P:
 [[0.5488135  0.71518937]
 [0.60276338 0.54488318]
 [0.4236548  0.64589411]
 [0.43758721 0.891773  ]
 [0.96366276 0.38344152]
 [0.79172504 0.52889492]
 [0.56804456 0.92559664]
 [0.07103606 0.0871293 ]
 [0.0202184  0.83261985]
 [0.77815675 0.87001215]]
P(1000): 
 [[0.41916695 0.66633421]
 [0.7845857  0.29828616]
 [0.         0.36688618]
 [0.17622127 1.        ]
 [1.         0.        ]
 [1.         0.12007927]
 [0.72441852 1.        ]
 [0.         0.        ]
 [0.         0.93898401]
 [1.         1.        ]]
simulate_cy execution time: 0.42031910000150674

Generated by Cython 0.29.33

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

+001: #%%cython -a

  __pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_4) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 002: 

+003: import numpy as np

  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 3, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 004: 

 005: cimport numpy as np

 006: 

 007: from cpython cimport array

 008: 

+009: import array

  __pyx_t_1 = __Pyx_Import(__pyx_n_s_array, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_array, __pyx_t_1) < 0) __PYX_ERR(0, 9, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 010: 

 011: import cython

 012: 

+013: import random

  __pyx_t_1 = __Pyx_Import(__pyx_n_s_random, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_random, __pyx_t_1) < 0) __PYX_ERR(0, 13, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 014: 

+015: import math

  __pyx_t_1 = __Pyx_Import(__pyx_n_s_math, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_math, __pyx_t_1) < 0) __PYX_ERR(0, 15, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 016: 

+017: import matplotlib.pyplot as plt

  __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 17, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF(__pyx_n_s__7);
  __Pyx_GIVEREF(__pyx_n_s__7);
  PyList_SET_ITEM(__pyx_t_1, 0, __pyx_n_s__7);
  __pyx_t_2 = __Pyx_Import(__pyx_n_s_matplotlib_pyplot, __pyx_t_1, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 17, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_plt, __pyx_t_2) < 0) __PYX_ERR(0, 17, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 018: 

+019: import time

  __pyx_t_2 = __Pyx_Import(__pyx_n_s_time, 0, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_time, __pyx_t_2) < 0) __PYX_ERR(0, 19, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 020: 

+021: import timeit

  __pyx_t_2 = __Pyx_Import(__pyx_n_s_timeit, 0, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 21, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_timeit, __pyx_t_2) < 0) __PYX_ERR(0, 21, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 022: 

 023: from cython.parallel import prange

 024: 

 025: 

 026: # random particles

+027: cdef initial_particles(int n):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_initial_particles(int __pyx_v_n) {
  PyArrayObject *__pyx_v_p = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_p;
  __Pyx_Buffer __pyx_pybuffer_p;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("initial_particles", 0);
  __pyx_pybuffer_p.pybuffer.buf = NULL;
  __pyx_pybuffer_p.refcount = 0;
  __pyx_pybuffernd_p.data = NULL;
  __pyx_pybuffernd_p.rcbuffer = &__pyx_pybuffer_p;
/* … */
  /* 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_6);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.initial_particles", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_p);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 028: 

+029:     np.random.seed(0)

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 29, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_random); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 29, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_seed); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 29, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && likely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
    }
  }
  __pyx_t_1 = (__pyx_t_3) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_3, __pyx_int_0) : __Pyx_PyObject_CallOneArg(__pyx_t_2, __pyx_int_0);
  __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 29, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 030: 

 031:     cdef np.ndarray[np.float64_t, ndim=2] p

+032:     p = np.random.rand(n, 2)

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 32, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_random); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 32, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_rand); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 32, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 32, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = NULL;
  __pyx_t_5 = 0;
  if (CYTHON_UNPACK_METHODS && likely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
      __pyx_t_5 = 1;
    }
  }
  #if CYTHON_FAST_PYCALL
  if (PyFunction_Check(__pyx_t_2)) {
    PyObject *__pyx_temp[3] = {__pyx_t_4, __pyx_t_3, __pyx_int_2};
    __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_2, __pyx_temp+1-__pyx_t_5, 2+__pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 32, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  } else
  #endif
  #if CYTHON_FAST_PYCCALL
  if (__Pyx_PyFastCFunction_Check(__pyx_t_2)) {
    PyObject *__pyx_temp[3] = {__pyx_t_4, __pyx_t_3, __pyx_int_2};
    __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_2, __pyx_temp+1-__pyx_t_5, 2+__pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 32, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  } else
  #endif
  {
    __pyx_t_6 = PyTuple_New(2+__pyx_t_5); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 32, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    if (__pyx_t_4) {
      __Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_4); __pyx_t_4 = NULL;
    }
    __Pyx_GIVEREF(__pyx_t_3);
    PyTuple_SET_ITEM(__pyx_t_6, 0+__pyx_t_5, __pyx_t_3);
    __Pyx_INCREF(__pyx_int_2);
    __Pyx_GIVEREF(__pyx_int_2);
    PyTuple_SET_ITEM(__pyx_t_6, 1+__pyx_t_5, __pyx_int_2);
    __pyx_t_3 = 0;
    __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_6, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 32, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  }
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 32, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
    __pyx_t_5 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_5 < 0)) {
      PyErr_Fetch(&__pyx_t_8, &__pyx_t_9, &__pyx_t_10);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_v_p, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_8); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_10);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_8, __pyx_t_9, __pyx_t_10);
      }
      __pyx_t_8 = __pyx_t_9 = __pyx_t_10 = 0;
    }
    __pyx_pybuffernd_p.diminfo[0].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_p.diminfo[0].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_p.diminfo[1].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_p.diminfo[1].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[1];
    if (unlikely(__pyx_t_5 < 0)) __PYX_ERR(0, 32, __pyx_L1_error)
  }
  __pyx_t_7 = 0;
  __pyx_v_p = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 033: 

+034:     return p

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_p));
  __pyx_r = ((PyObject *)__pyx_v_p);
  goto __pyx_L0;

 035: 

 036: 

 037: cimport cython

 038: @cython.boundscheck(False)  # Deactivate bounds checking

 039: @cython.wraparound(False)   # Deactivate negative indexing.

+040: cdef simulate_cy(int n, int time_step, show_plot=False):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_simulate_cy(int __pyx_v_n, int __pyx_v_time_step, struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_simulate_cy *__pyx_optional_args) {
  PyObject *__pyx_v_show_plot = ((PyObject *)Py_False);
  PyArrayObject *__pyx_v_p = 0;
  int __pyx_v_i;
  PyArrayObject *__pyx_v_total_force = 0;
  PyArrayObject *__pyx_v_x_det = 0;
  PyArrayObject *__pyx_v_pos = 0;
  PyObject *__pyx_v_colors = NULL;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_p;
  __Pyx_Buffer __pyx_pybuffer_p;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_pos;
  __Pyx_Buffer __pyx_pybuffer_pos;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_total_force;
  __Pyx_Buffer __pyx_pybuffer_total_force;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x_det;
  __Pyx_Buffer __pyx_pybuffer_x_det;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("simulate_cy", 0);
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_show_plot = __pyx_optional_args->show_plot;
    }
  }
  __pyx_pybuffer_p.pybuffer.buf = NULL;
  __pyx_pybuffer_p.refcount = 0;
  __pyx_pybuffernd_p.data = NULL;
  __pyx_pybuffernd_p.rcbuffer = &__pyx_pybuffer_p;
  __pyx_pybuffer_total_force.pybuffer.buf = NULL;
  __pyx_pybuffer_total_force.refcount = 0;
  __pyx_pybuffernd_total_force.data = NULL;
  __pyx_pybuffernd_total_force.rcbuffer = &__pyx_pybuffer_total_force;
  __pyx_pybuffer_x_det.pybuffer.buf = NULL;
  __pyx_pybuffer_x_det.refcount = 0;
  __pyx_pybuffernd_x_det.data = NULL;
  __pyx_pybuffernd_x_det.rcbuffer = &__pyx_pybuffer_x_det;
  __pyx_pybuffer_pos.pybuffer.buf = NULL;
  __pyx_pybuffer_pos.refcount = 0;
  __pyx_pybuffernd_pos.data = NULL;
  __pyx_pybuffernd_pos.rcbuffer = &__pyx_pybuffer_pos;
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_12);
  __Pyx_XDECREF(__pyx_t_13);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pos.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x_det.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.simulate_cy", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pos.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x_det.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_p);
  __Pyx_XDECREF((PyObject *)__pyx_v_total_force);
  __Pyx_XDECREF((PyObject *)__pyx_v_x_det);
  __Pyx_XDECREF((PyObject *)__pyx_v_pos);
  __Pyx_XDECREF(__pyx_v_colors);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_simulate_cy {
  int __pyx_n;
  PyObject *show_plot;
};

 041: 

 042:     cdef np.ndarray[np.float64_t, ndim=2] p

 043: 

+044:     p = initial_particles(n)

  __pyx_t_1 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_initial_particles(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 44, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 44, __pyx_L1_error)
  __pyx_t_2 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
    __pyx_t_3 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_t_2, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_3 < 0)) {
      PyErr_Fetch(&__pyx_t_4, &__pyx_t_5, &__pyx_t_6);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_v_p, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_4); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_6);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_4, __pyx_t_5, __pyx_t_6);
      }
      __pyx_t_4 = __pyx_t_5 = __pyx_t_6 = 0;
    }
    __pyx_pybuffernd_p.diminfo[0].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_p.diminfo[0].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_p.diminfo[1].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_p.diminfo[1].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[1];
    if (unlikely(__pyx_t_3 < 0)) __PYX_ERR(0, 44, __pyx_L1_error)
  }
  __pyx_t_2 = 0;
  __pyx_v_p = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+045:     print("initial_P:\n",p)

  __pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 45, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF(__pyx_kp_u_initial_P);
  __Pyx_GIVEREF(__pyx_kp_u_initial_P);
  PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_kp_u_initial_P);
  __Pyx_INCREF(((PyObject *)__pyx_v_p));
  __Pyx_GIVEREF(((PyObject *)__pyx_v_p));
  PyTuple_SET_ITEM(__pyx_t_1, 1, ((PyObject *)__pyx_v_p));
  __pyx_t_7 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_t_1, NULL); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 45, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;

 046: 

 047:     cdef int i

 048:     cdef np.ndarray[np.float64_t, ndim=2] total_force, x_det, pos

 049: 

+050:     for i in range(time_step):

  __pyx_t_3 = __pyx_v_time_step;
  __pyx_t_8 = __pyx_t_3;
  for (__pyx_t_9 = 0; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
    __pyx_v_i = __pyx_t_9;

 051: 

+052:         total_force = combined_force(p, n)

    __pyx_t_7 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_combined_force(((PyArrayObject *)__pyx_v_p), __pyx_v_n); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 52, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    if (!(likely(((__pyx_t_7) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_7, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 52, __pyx_L1_error)
    __pyx_t_10 = ((PyArrayObject *)__pyx_t_7);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
      __pyx_t_11 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_t_10, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_11 < 0)) {
        PyErr_Fetch(&__pyx_t_6, &__pyx_t_5, &__pyx_t_4);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_v_total_force, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_6); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_4);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_6, __pyx_t_5, __pyx_t_4);
        }
        __pyx_t_6 = __pyx_t_5 = __pyx_t_4 = 0;
      }
      __pyx_pybuffernd_total_force.diminfo[0].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_total_force.diminfo[0].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_total_force.diminfo[1].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_total_force.diminfo[1].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[1];
      if (unlikely(__pyx_t_11 < 0)) __PYX_ERR(0, 52, __pyx_L1_error)
    }
    __pyx_t_10 = 0;
    __Pyx_XDECREF_SET(__pyx_v_total_force, ((PyArrayObject *)__pyx_t_7));
    __pyx_t_7 = 0;

 053: 

+054:         x_det = displacement(total_force, delta_t=0.0001)

    __Pyx_GetModuleGlobalName(__pyx_t_7, __pyx_n_s_displacement); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 54, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 54, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_INCREF(((PyObject *)__pyx_v_total_force));
    __Pyx_GIVEREF(((PyObject *)__pyx_v_total_force));
    PyTuple_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_v_total_force));
    __pyx_t_12 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 54, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_12);
    if (PyDict_SetItem(__pyx_t_12, __pyx_n_s_delta_t, __pyx_float_0_0001) < 0) __PYX_ERR(0, 54, __pyx_L1_error)
    __pyx_t_13 = __Pyx_PyObject_Call(__pyx_t_7, __pyx_t_1, __pyx_t_12); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 54, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_13);
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
    if (!(likely(((__pyx_t_13) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_13, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 54, __pyx_L1_error)
    __pyx_t_10 = ((PyArrayObject *)__pyx_t_13);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x_det.rcbuffer->pybuffer);
      __pyx_t_11 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x_det.rcbuffer->pybuffer, (PyObject*)__pyx_t_10, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_11 < 0)) {
        PyErr_Fetch(&__pyx_t_4, &__pyx_t_5, &__pyx_t_6);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x_det.rcbuffer->pybuffer, (PyObject*)__pyx_v_x_det, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_4); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_6);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_4, __pyx_t_5, __pyx_t_6);
        }
        __pyx_t_4 = __pyx_t_5 = __pyx_t_6 = 0;
      }
      __pyx_pybuffernd_x_det.diminfo[0].strides = __pyx_pybuffernd_x_det.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x_det.diminfo[0].shape = __pyx_pybuffernd_x_det.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x_det.diminfo[1].strides = __pyx_pybuffernd_x_det.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x_det.diminfo[1].shape = __pyx_pybuffernd_x_det.rcbuffer->pybuffer.shape[1];
      if (unlikely(__pyx_t_11 < 0)) __PYX_ERR(0, 54, __pyx_L1_error)
    }
    __pyx_t_10 = 0;
    __Pyx_XDECREF_SET(__pyx_v_x_det, ((PyArrayObject *)__pyx_t_13));
    __pyx_t_13 = 0;

 055: 

 056:         # update the particles

+057:         p = update_position(p, x_det)

    __pyx_t_13 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_update_position(((PyArrayObject *)__pyx_v_p), ((PyArrayObject *)__pyx_v_x_det), NULL); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 57, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_13);
    if (!(likely(((__pyx_t_13) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_13, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 57, __pyx_L1_error)
    __pyx_t_2 = ((PyArrayObject *)__pyx_t_13);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
      __pyx_t_11 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_t_2, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_11 < 0)) {
        PyErr_Fetch(&__pyx_t_6, &__pyx_t_5, &__pyx_t_4);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_v_p, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_6); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_4);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_6, __pyx_t_5, __pyx_t_4);
        }
        __pyx_t_6 = __pyx_t_5 = __pyx_t_4 = 0;
      }
      __pyx_pybuffernd_p.diminfo[0].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_p.diminfo[0].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_p.diminfo[1].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_p.diminfo[1].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[1];
      if (unlikely(__pyx_t_11 < 0)) __PYX_ERR(0, 57, __pyx_L1_error)
    }
    __pyx_t_2 = 0;
    __Pyx_DECREF_SET(__pyx_v_p, ((PyArrayObject *)__pyx_t_13));
    __pyx_t_13 = 0;

 058: 

+059:         pos = p

    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_pos.rcbuffer->pybuffer);
      __pyx_t_11 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_pos.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_v_p), &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_11 < 0)) {
        PyErr_Fetch(&__pyx_t_4, &__pyx_t_5, &__pyx_t_6);
        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)) {
          Py_XDECREF(__pyx_t_4); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_6);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_4, __pyx_t_5, __pyx_t_6);
        }
        __pyx_t_4 = __pyx_t_5 = __pyx_t_6 = 0;
      }
      __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];
      if (unlikely(__pyx_t_11 < 0)) __PYX_ERR(0, 59, __pyx_L1_error)
    }
    __Pyx_INCREF(((PyObject *)__pyx_v_p));
    __Pyx_XDECREF_SET(__pyx_v_pos, ((PyArrayObject *)__pyx_v_p));

 060: 

 061:         # plot particles

+062:         colors = ['red', 'green', 'blue', 'orange']

    __pyx_t_13 = PyList_New(4); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 62, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_13);
    __Pyx_INCREF(__pyx_n_u_red);
    __Pyx_GIVEREF(__pyx_n_u_red);
    PyList_SET_ITEM(__pyx_t_13, 0, __pyx_n_u_red);
    __Pyx_INCREF(__pyx_n_u_green);
    __Pyx_GIVEREF(__pyx_n_u_green);
    PyList_SET_ITEM(__pyx_t_13, 1, __pyx_n_u_green);
    __Pyx_INCREF(__pyx_n_u_blue);
    __Pyx_GIVEREF(__pyx_n_u_blue);
    PyList_SET_ITEM(__pyx_t_13, 2, __pyx_n_u_blue);
    __Pyx_INCREF(__pyx_n_u_orange);
    __Pyx_GIVEREF(__pyx_n_u_orange);
    PyList_SET_ITEM(__pyx_t_13, 3, __pyx_n_u_orange);
    __Pyx_XDECREF_SET(__pyx_v_colors, ((PyObject*)__pyx_t_13));
    __pyx_t_13 = 0;

+063:         if show_plot:

    __pyx_t_14 = __Pyx_PyObject_IsTrue(__pyx_v_show_plot); if (unlikely(__pyx_t_14 < 0)) __PYX_ERR(0, 63, __pyx_L1_error)
    if (__pyx_t_14) {
/* … */
    }
  }

 064: 

+065:             if i % 2 == 0:

      __pyx_t_14 = ((__Pyx_mod_long(__pyx_v_i, 2) == 0) != 0);
      if (__pyx_t_14) {
/* … */
      }

 066: 

+067:                 update_plot(pos,colors)

        __pyx_t_13 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_update_plot(((PyObject *)__pyx_v_pos), __pyx_v_colors); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 67, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_13);
        __Pyx_DECREF(__pyx_t_13); __pyx_t_13 = 0;

 068: 

+069:     print("P({}): \n".format(time_step), p)

  __pyx_t_12 = __Pyx_PyObject_GetAttrStr(__pyx_kp_u_P, __pyx_n_s_format); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 69, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_12);
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_time_step); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 69, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_7 = NULL;
  if (CYTHON_UNPACK_METHODS && likely(PyMethod_Check(__pyx_t_12))) {
    __pyx_t_7 = PyMethod_GET_SELF(__pyx_t_12);
    if (likely(__pyx_t_7)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_12);
      __Pyx_INCREF(__pyx_t_7);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_12, function);
    }
  }
  __pyx_t_13 = (__pyx_t_7) ? __Pyx_PyObject_Call2Args(__pyx_t_12, __pyx_t_7, __pyx_t_1) : __Pyx_PyObject_CallOneArg(__pyx_t_12, __pyx_t_1);
  __Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 69, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_13);
  __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
  __pyx_t_12 = PyTuple_New(2); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 69, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_12);
  __Pyx_GIVEREF(__pyx_t_13);
  PyTuple_SET_ITEM(__pyx_t_12, 0, __pyx_t_13);
  __Pyx_INCREF(((PyObject *)__pyx_v_p));
  __Pyx_GIVEREF(((PyObject *)__pyx_v_p));
  PyTuple_SET_ITEM(__pyx_t_12, 1, ((PyObject *)__pyx_v_p));
  __pyx_t_13 = 0;
  __pyx_t_13 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_t_12, NULL); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 69, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_13);
  __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
  __Pyx_DECREF(__pyx_t_13); __pyx_t_13 = 0;

 070: 

+071:     return p

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_p));
  __pyx_r = ((PyObject *)__pyx_v_p);
  goto __pyx_L0;

 072: 

 073: 

 074: # Calculate the strength of the repulsion

+075: cdef f(np.ndarray[np.float64_t, ndim=1]r, double c1=1, double c2=1):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_f(PyArrayObject *__pyx_v_r, struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_f *__pyx_optional_args) {
  double __pyx_v_c1 = ((double)1.0);
  double __pyx_v_mag;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_r;
  __Pyx_Buffer __pyx_pybuffer_r;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("f", 0);
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_c1 = __pyx_optional_args->c1;
    }
  }
  __pyx_pybuffer_r.pybuffer.buf = NULL;
  __pyx_pybuffer_r.refcount = 0;
  __pyx_pybuffernd_r.data = NULL;
  __pyx_pybuffernd_r.rcbuffer = &__pyx_pybuffer_r;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_r.rcbuffer->pybuffer, (PyObject*)__pyx_v_r, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 75, __pyx_L1_error)
  }
  __pyx_pybuffernd_r.diminfo[0].strides = __pyx_pybuffernd_r.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_r.diminfo[0].shape = __pyx_pybuffernd_r.rcbuffer->pybuffer.shape[0];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_r.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.f", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_r.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_f {
  int __pyx_n;
  double c1;
  double c2;
};

 076: 

 077:     # Calculate the force

 078: 

 079:     cdef double mag

 080: 

+081:     mag = c1 / math.sqrt(square_sum(r))

  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_c1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_math); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_square_sum(((PyArrayObject *)__pyx_v_r)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_5, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3);
  __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyNumber_Divide(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 81, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_mag = __pyx_t_6;

 082: 

+083:     return mag

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_4 = PyFloat_FromDouble(__pyx_v_mag); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 83, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_r = __pyx_t_4;
  __pyx_t_4 = 0;
  goto __pyx_L0;

 084: 

 085: 

 086: 

 087: cimport cython

 088: @cython.boundscheck(False)

 089: @cython.wraparound(False)

 090: 

 091: #calculate square

+092: cdef square_sum(np.ndarray[np.float64_t, ndim=1] arr):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_square_sum(PyArrayObject *__pyx_v_arr) {
  Py_ssize_t __pyx_v_i;
  CYTHON_UNUSED Py_ssize_t __pyx_v_n;
  double __pyx_v_s;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_arr;
  __Pyx_Buffer __pyx_pybuffer_arr;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("square_sum", 0);
  __pyx_pybuffer_arr.pybuffer.buf = NULL;
  __pyx_pybuffer_arr.refcount = 0;
  __pyx_pybuffernd_arr.data = NULL;
  __pyx_pybuffernd_arr.rcbuffer = &__pyx_pybuffer_arr;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_arr.rcbuffer->pybuffer, (PyObject*)__pyx_v_arr, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 92, __pyx_L1_error)
  }
  __pyx_pybuffernd_arr.diminfo[0].strides = __pyx_pybuffernd_arr.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_arr.diminfo[0].shape = __pyx_pybuffernd_arr.rcbuffer->pybuffer.shape[0];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_6);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_arr.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.square_sum", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_arr.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 093:     cdef Py_ssize_t i, n

+094:     cdef double s = 0

  __pyx_v_s = 0.0;

+095:     n = arr.shape[0]

  __pyx_v_n = (__pyx_v_arr->dimensions[0]);

+096:     for i in prange(n, nogil=True, num_threads=1):

  {
      #ifdef WITH_THREAD
      PyThreadState *_save;
      Py_UNBLOCK_THREADS
      __Pyx_FastGIL_Remember();
      #endif
      /*try:*/ {
        __pyx_t_1 = __pyx_v_n;
        if ((1 == 0)) abort();
        {
            #if ((defined(__APPLE__) || defined(__OSX__)) && (defined(__GNUC__) && (__GNUC__ > 2 || (__GNUC__ == 2 && (__GNUC_MINOR__ > 95)))))
                #undef likely
                #undef unlikely
                #define likely(x)   (x)
                #define unlikely(x) (x)
            #endif
            __pyx_t_3 = (__pyx_t_1 - 0 + 1 - 1/abs(1)) / 1;
            if (__pyx_t_3 > 0)
            {
                #ifdef _OPENMP
                #pragma omp parallel
                #endif /* _OPENMP */
                {
                    #ifdef _OPENMP
                    #pragma omp for firstprivate(__pyx_v_i) lastprivate(__pyx_v_i) reduction(+:__pyx_v_s) num_threads(1)
                    #endif /* _OPENMP */
                    for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_3; __pyx_t_2++){
                        {
                            __pyx_v_i = (Py_ssize_t)(0 + 1 * __pyx_t_2);
/* … */
      /*finally:*/ {
        /*normal exit:*/{
          #ifdef WITH_THREAD
          __Pyx_FastGIL_Forget();
          Py_BLOCK_THREADS
          #endif
          goto __pyx_L5;
        }
        __pyx_L5:;
      }
  }

+097:         s += arr[i] * arr[i]

                            __pyx_t_4 = __pyx_v_i;
                            __pyx_t_5 = __pyx_v_i;
                            __pyx_v_s = (__pyx_v_s + ((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_arr.rcbuffer->pybuffer.buf, __pyx_t_4, __pyx_pybuffernd_arr.diminfo[0].strides)) * (*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_arr.rcbuffer->pybuffer.buf, __pyx_t_5, __pyx_pybuffernd_arr.diminfo[0].strides))));
                        }
                    }
                }
            }
        }
        #if ((defined(__APPLE__) || defined(__OSX__)) && (defined(__GNUC__) && (__GNUC__ > 2 || (__GNUC__ == 2 && (__GNUC_MINOR__ > 95)))))
            #undef likely
            #undef unlikely
            #define likely(x)   __builtin_expect(!!(x), 1)
            #define unlikely(x) __builtin_expect(!!(x), 0)
        #endif
      }

 098: #    print(type(s))

+099:     return s

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_6 = PyFloat_FromDouble(__pyx_v_s); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 99, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_r = __pyx_t_6;
  __pyx_t_6 = 0;
  goto __pyx_L0;

 100: 

 101: 

 102: @cython.boundscheck(False)

 103: @cython.wraparound(False)

 104: # Calculate the total force for each particle

+105: cdef combined_force(np.ndarray[np.float64_t, ndim=2] p, int n):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_combined_force(PyArrayObject *__pyx_v_p, int __pyx_v_n) {
  PyArrayObject *__pyx_v_total_force = 0;
  Py_ssize_t __pyx_v_i;
  Py_ssize_t __pyx_v_j;
  PyArrayObject *__pyx_v_r = 0;
  PyArrayObject *__pyx_v_fn_sum = 0;
  PyArrayObject *__pyx_v_fn = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_fn;
  __Pyx_Buffer __pyx_pybuffer_fn;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_fn_sum;
  __Pyx_Buffer __pyx_pybuffer_fn_sum;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_p;
  __Pyx_Buffer __pyx_pybuffer_p;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_r;
  __Pyx_Buffer __pyx_pybuffer_r;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_total_force;
  __Pyx_Buffer __pyx_pybuffer_total_force;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("combined_force", 0);
  __pyx_pybuffer_total_force.pybuffer.buf = NULL;
  __pyx_pybuffer_total_force.refcount = 0;
  __pyx_pybuffernd_total_force.data = NULL;
  __pyx_pybuffernd_total_force.rcbuffer = &__pyx_pybuffer_total_force;
  __pyx_pybuffer_r.pybuffer.buf = NULL;
  __pyx_pybuffer_r.refcount = 0;
  __pyx_pybuffernd_r.data = NULL;
  __pyx_pybuffernd_r.rcbuffer = &__pyx_pybuffer_r;
  __pyx_pybuffer_fn_sum.pybuffer.buf = NULL;
  __pyx_pybuffer_fn_sum.refcount = 0;
  __pyx_pybuffernd_fn_sum.data = NULL;
  __pyx_pybuffernd_fn_sum.rcbuffer = &__pyx_pybuffer_fn_sum;
  __pyx_pybuffer_fn.pybuffer.buf = NULL;
  __pyx_pybuffer_fn.refcount = 0;
  __pyx_pybuffernd_fn.data = NULL;
  __pyx_pybuffernd_fn.rcbuffer = &__pyx_pybuffer_fn;
  __pyx_pybuffer_p.pybuffer.buf = NULL;
  __pyx_pybuffer_p.refcount = 0;
  __pyx_pybuffernd_p.data = NULL;
  __pyx_pybuffernd_p.rcbuffer = &__pyx_pybuffer_p;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_v_p, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 105, __pyx_L1_error)
  }
  __pyx_pybuffernd_p.diminfo[0].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_p.diminfo[0].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_p.diminfo[1].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_p.diminfo[1].shape = __pyx_pybuffernd_p.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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_r.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.combined_force", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_r.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_total_force);
  __Pyx_XDECREF((PyObject *)__pyx_v_r);
  __Pyx_XDECREF((PyObject *)__pyx_v_fn_sum);
  __Pyx_XDECREF((PyObject *)__pyx_v_fn);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 106: 

 107:     cdef np.ndarray[np.float64_t, ndim=2] total_force

 108: 

+109:     total_force= np.zeros_like(p,dtype=np.double)

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros_like); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF(((PyObject *)__pyx_v_p));
  __Pyx_GIVEREF(((PyObject *)__pyx_v_p));
  PyTuple_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_v_p));
  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, __pyx_t_5) < 0) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 109, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __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_ERR(0, 109, __pyx_L1_error)
  __pyx_t_6 = ((PyArrayObject *)__pyx_t_5);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
    __pyx_t_7 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_t_6, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_7 < 0)) {
      PyErr_Fetch(&__pyx_t_8, &__pyx_t_9, &__pyx_t_10);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_v_total_force, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_8); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_10);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_8, __pyx_t_9, __pyx_t_10);
      }
      __pyx_t_8 = __pyx_t_9 = __pyx_t_10 = 0;
    }
    __pyx_pybuffernd_total_force.diminfo[0].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_total_force.diminfo[0].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_total_force.diminfo[1].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_total_force.diminfo[1].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[1];
    if (unlikely(__pyx_t_7 < 0)) __PYX_ERR(0, 109, __pyx_L1_error)
  }
  __pyx_t_6 = 0;
  __pyx_v_total_force = ((PyArrayObject *)__pyx_t_5);
  __pyx_t_5 = 0;

 110: 

 111:     cdef Py_ssize_t i,j

 112: 

 113:     cdef np.ndarray[np.float64_t, ndim=1] r, fn_sum, fn

 114: 

+115:     for i in range(n):

  __pyx_t_7 = __pyx_v_n;
  __pyx_t_11 = __pyx_t_7;
  for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
    __pyx_v_i = __pyx_t_12;

 116: 

+117:         fn_sum = np.zeros(2,dtype=np.double)

    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_5 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_double); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    if (PyDict_SetItem(__pyx_t_5, __pyx_n_s_dtype, __pyx_t_2) < 0) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_tuple_, __pyx_t_5); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 117, __pyx_L1_error)
    __pyx_t_13 = ((PyArrayObject *)__pyx_t_2);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer);
      __pyx_t_14 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer, (PyObject*)__pyx_t_13, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_14 < 0)) {
        PyErr_Fetch(&__pyx_t_10, &__pyx_t_9, &__pyx_t_8);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer, (PyObject*)__pyx_v_fn_sum, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_10); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_8);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_10, __pyx_t_9, __pyx_t_8);
        }
        __pyx_t_10 = __pyx_t_9 = __pyx_t_8 = 0;
      }
      __pyx_pybuffernd_fn_sum.diminfo[0].strides = __pyx_pybuffernd_fn_sum.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_fn_sum.diminfo[0].shape = __pyx_pybuffernd_fn_sum.rcbuffer->pybuffer.shape[0];
      if (unlikely(__pyx_t_14 < 0)) __PYX_ERR(0, 117, __pyx_L1_error)
    }
    __pyx_t_13 = 0;
    __Pyx_XDECREF_SET(__pyx_v_fn_sum, ((PyArrayObject *)__pyx_t_2));
    __pyx_t_2 = 0;
/* … */
  __pyx_tuple_ = PyTuple_Pack(1, __pyx_int_2); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 117, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple_);
  __Pyx_GIVEREF(__pyx_tuple_);

 118: 

+119:         for j in range(n):

    __pyx_t_14 = __pyx_v_n;
    __pyx_t_15 = __pyx_t_14;
    for (__pyx_t_16 = 0; __pyx_t_16 < __pyx_t_15; __pyx_t_16+=1) {
      __pyx_v_j = __pyx_t_16;

 120: 

+121:             if j != i:

      __pyx_t_17 = ((__pyx_v_j != __pyx_v_i) != 0);
      if (__pyx_t_17) {
/* … */
      }

 122: 

+123:                 r = p[j] - p[i]

        __pyx_t_2 = __Pyx_GetItemInt(((PyObject *)__pyx_v_p), __pyx_v_j, Py_ssize_t, 1, PyInt_FromSsize_t, 0, 0, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 123, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_2);
        __pyx_t_5 = __Pyx_GetItemInt(((PyObject *)__pyx_v_p), __pyx_v_i, Py_ssize_t, 1, PyInt_FromSsize_t, 0, 0, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 123, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __pyx_t_3 = PyNumber_Subtract(__pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 123, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
        __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
        if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 123, __pyx_L1_error)
        __pyx_t_13 = ((PyArrayObject *)__pyx_t_3);
        {
          __Pyx_BufFmt_StackElem __pyx_stack[1];
          __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_r.rcbuffer->pybuffer);
          __pyx_t_18 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_r.rcbuffer->pybuffer, (PyObject*)__pyx_t_13, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
          if (unlikely(__pyx_t_18 < 0)) {
            PyErr_Fetch(&__pyx_t_8, &__pyx_t_9, &__pyx_t_10);
            if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_r.rcbuffer->pybuffer, (PyObject*)__pyx_v_r, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
              Py_XDECREF(__pyx_t_8); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_10);
              __Pyx_RaiseBufferFallbackError();
            } else {
              PyErr_Restore(__pyx_t_8, __pyx_t_9, __pyx_t_10);
            }
            __pyx_t_8 = __pyx_t_9 = __pyx_t_10 = 0;
          }
          __pyx_pybuffernd_r.diminfo[0].strides = __pyx_pybuffernd_r.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_r.diminfo[0].shape = __pyx_pybuffernd_r.rcbuffer->pybuffer.shape[0];
          if (unlikely(__pyx_t_18 < 0)) __PYX_ERR(0, 123, __pyx_L1_error)
        }
        __pyx_t_13 = 0;
        __Pyx_XDECREF_SET(__pyx_v_r, ((PyArrayObject *)__pyx_t_3));
        __pyx_t_3 = 0;

 124: 

+125:                 fn =  -1 * f(r) * np.sign(r)

        __pyx_t_3 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_f(((PyArrayObject *)__pyx_v_r), NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        __pyx_t_5 = PyNumber_Multiply(__pyx_int_neg_1, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
        __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_2);
        __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_sign); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_1);
        __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
        __pyx_t_2 = NULL;
        if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
          __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_1);
          if (likely(__pyx_t_2)) {
            PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
            __Pyx_INCREF(__pyx_t_2);
            __Pyx_INCREF(function);
            __Pyx_DECREF_SET(__pyx_t_1, function);
          }
        }
        __pyx_t_3 = (__pyx_t_2) ? __Pyx_PyObject_Call2Args(__pyx_t_1, __pyx_t_2, ((PyObject *)__pyx_v_r)) : __Pyx_PyObject_CallOneArg(__pyx_t_1, ((PyObject *)__pyx_v_r));
        __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
        if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
        __pyx_t_1 = PyNumber_Multiply(__pyx_t_5, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 125, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_1);
        __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 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_ERR(0, 125, __pyx_L1_error)
        __pyx_t_13 = ((PyArrayObject *)__pyx_t_1);
        {
          __Pyx_BufFmt_StackElem __pyx_stack[1];
          __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn.rcbuffer->pybuffer);
          __pyx_t_18 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn.rcbuffer->pybuffer, (PyObject*)__pyx_t_13, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
          if (unlikely(__pyx_t_18 < 0)) {
            PyErr_Fetch(&__pyx_t_10, &__pyx_t_9, &__pyx_t_8);
            if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn.rcbuffer->pybuffer, (PyObject*)__pyx_v_fn, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
              Py_XDECREF(__pyx_t_10); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_8);
              __Pyx_RaiseBufferFallbackError();
            } else {
              PyErr_Restore(__pyx_t_10, __pyx_t_9, __pyx_t_8);
            }
            __pyx_t_10 = __pyx_t_9 = __pyx_t_8 = 0;
          }
          __pyx_pybuffernd_fn.diminfo[0].strides = __pyx_pybuffernd_fn.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_fn.diminfo[0].shape = __pyx_pybuffernd_fn.rcbuffer->pybuffer.shape[0];
          if (unlikely(__pyx_t_18 < 0)) __PYX_ERR(0, 125, __pyx_L1_error)
        }
        __pyx_t_13 = 0;
        __Pyx_XDECREF_SET(__pyx_v_fn, ((PyArrayObject *)__pyx_t_1));
        __pyx_t_1 = 0;

 126: 

+127:                 fn_sum += fn

        __pyx_t_1 = PyNumber_InPlaceAdd(((PyObject *)__pyx_v_fn_sum), ((PyObject *)__pyx_v_fn)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 127, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_1);
        if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 127, __pyx_L1_error)
        __pyx_t_13 = ((PyArrayObject *)__pyx_t_1);
        {
          __Pyx_BufFmt_StackElem __pyx_stack[1];
          __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer);
          __pyx_t_18 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer, (PyObject*)__pyx_t_13, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
          if (unlikely(__pyx_t_18 < 0)) {
            PyErr_Fetch(&__pyx_t_8, &__pyx_t_9, &__pyx_t_10);
            if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_fn_sum.rcbuffer->pybuffer, (PyObject*)__pyx_v_fn_sum, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
              Py_XDECREF(__pyx_t_8); Py_XDECREF(__pyx_t_9); Py_XDECREF(__pyx_t_10);
              __Pyx_RaiseBufferFallbackError();
            } else {
              PyErr_Restore(__pyx_t_8, __pyx_t_9, __pyx_t_10);
            }
            __pyx_t_8 = __pyx_t_9 = __pyx_t_10 = 0;
          }
          __pyx_pybuffernd_fn_sum.diminfo[0].strides = __pyx_pybuffernd_fn_sum.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_fn_sum.diminfo[0].shape = __pyx_pybuffernd_fn_sum.rcbuffer->pybuffer.shape[0];
          if (unlikely(__pyx_t_18 < 0)) __PYX_ERR(0, 127, __pyx_L1_error)
        }
        __pyx_t_13 = 0;
        __Pyx_DECREF_SET(__pyx_v_fn_sum, ((PyArrayObject *)__pyx_t_1));
        __pyx_t_1 = 0;

 128: 

+129:             total_force[i] = fn_sum

      if (unlikely(__Pyx_SetItemInt(((PyObject *)__pyx_v_total_force), __pyx_v_i, ((PyObject *)__pyx_v_fn_sum), Py_ssize_t, 1, PyInt_FromSsize_t, 0, 0, 0) < 0)) __PYX_ERR(0, 129, __pyx_L1_error)
    }
  }

 130: 

+131:     return total_force

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_total_force));
  __pyx_r = ((PyObject *)__pyx_v_total_force);
  goto __pyx_L0;

 132: 

 133: 

 134: 

 135: # Calculate the displacement between two particles

+136: cpdef displacement(np.ndarray[np.float64_t, ndim=2] total_force, double eta=1, double delta_t=1):

static PyObject *__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_1displacement(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_displacement(PyArrayObject *__pyx_v_total_force, CYTHON_UNUSED int __pyx_skip_dispatch, struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_displacement *__pyx_optional_args) {
  double __pyx_v_eta = ((double)1.0);
  double __pyx_v_delta_t = ((double)1.0);
  PyArrayObject *__pyx_v_displacement = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_displacement;
  __Pyx_Buffer __pyx_pybuffer_displacement;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_total_force;
  __Pyx_Buffer __pyx_pybuffer_total_force;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("displacement", 0);
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_eta = __pyx_optional_args->eta;
      if (__pyx_optional_args->__pyx_n > 1) {
        __pyx_v_delta_t = __pyx_optional_args->delta_t;
      }
    }
  }
  __pyx_pybuffer_displacement.pybuffer.buf = NULL;
  __pyx_pybuffer_displacement.refcount = 0;
  __pyx_pybuffernd_displacement.data = NULL;
  __pyx_pybuffernd_displacement.rcbuffer = &__pyx_pybuffer_displacement;
  __pyx_pybuffer_total_force.pybuffer.buf = NULL;
  __pyx_pybuffer_total_force.refcount = 0;
  __pyx_pybuffernd_total_force.data = NULL;
  __pyx_pybuffernd_total_force.rcbuffer = &__pyx_pybuffer_total_force;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_v_total_force, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 136, __pyx_L1_error)
  }
  __pyx_pybuffernd_total_force.diminfo[0].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_total_force.diminfo[0].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_total_force.diminfo[1].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_total_force.diminfo[1].shape = __pyx_pybuffernd_total_force.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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_displacement.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.displacement", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_displacement.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_displacement);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_1displacement(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyObject *__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_1displacement(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  PyArrayObject *__pyx_v_total_force = 0;
  double __pyx_v_eta;
  double __pyx_v_delta_t;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("displacement (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_total_force,&__pyx_n_s_eta,&__pyx_n_s_delta_t,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);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_total_force)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (kw_args > 0) {
          PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_eta);
          if (value) { values[1] = value; kw_args--; }
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (kw_args > 0) {
          PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_delta_t);
          if (value) { values[2] = value; kw_args--; }
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "displacement") < 0)) __PYX_ERR(0, 136, __pyx_L3_error)
      }
    } else {
      switch (PyTuple_GET_SIZE(__pyx_args)) {
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_total_force = ((PyArrayObject *)values[0]);
    if (values[1]) {
      __pyx_v_eta = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_eta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 136, __pyx_L3_error)
    } else {
      __pyx_v_eta = ((double)1.0);
    }
    if (values[2]) {
      __pyx_v_delta_t = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_delta_t == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 136, __pyx_L3_error)
    } else {
      __pyx_v_delta_t = ((double)1.0);
    }
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("displacement", 0, 1, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 136, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.displacement", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_total_force), __pyx_ptype_5numpy_ndarray, 1, "total_force", 0))) __PYX_ERR(0, 136, __pyx_L1_error)
  __pyx_r = __pyx_pf_46_cython_magic_5971c10418fa3492d0bfce342614b144_displacement(__pyx_self, __pyx_v_total_force, __pyx_v_eta, __pyx_v_delta_t);
  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_5971c10418fa3492d0bfce342614b144_displacement(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_total_force, double __pyx_v_eta, double __pyx_v_delta_t) {
  __Pyx_LocalBuf_ND __pyx_pybuffernd_total_force;
  __Pyx_Buffer __pyx_pybuffer_total_force;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("displacement", 0);
  __pyx_pybuffer_total_force.pybuffer.buf = NULL;
  __pyx_pybuffer_total_force.refcount = 0;
  __pyx_pybuffernd_total_force.data = NULL;
  __pyx_pybuffernd_total_force.rcbuffer = &__pyx_pybuffer_total_force;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer, (PyObject*)__pyx_v_total_force, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 136, __pyx_L1_error)
  }
  __pyx_pybuffernd_total_force.diminfo[0].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_total_force.diminfo[0].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_total_force.diminfo[1].strides = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_total_force.diminfo[1].shape = __pyx_pybuffernd_total_force.rcbuffer->pybuffer.shape[1];
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2.__pyx_n = 2;
  __pyx_t_2.eta = __pyx_v_eta;
  __pyx_t_2.delta_t = __pyx_v_delta_t;
  __pyx_t_1 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_displacement(__pyx_v_total_force, 0, &__pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 136, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.displacement", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_total_force.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_displacement {
  int __pyx_n;
  double eta;
  double delta_t;
};

 137:     cdef np.ndarray[np.float64_t, ndim=2] displacement

 138: 

+139:     displacement = total_force / eta * delta_t

  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_eta); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyNumber_Divide(((PyObject *)__pyx_v_total_force), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_delta_t); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = PyNumber_Multiply(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 139, __pyx_L1_error)
  __pyx_t_4 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_displacement.rcbuffer->pybuffer);
    __pyx_t_5 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_displacement.rcbuffer->pybuffer, (PyObject*)__pyx_t_4, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_5 < 0)) {
      PyErr_Fetch(&__pyx_t_6, &__pyx_t_7, &__pyx_t_8);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_displacement.rcbuffer->pybuffer, (PyObject*)__pyx_v_displacement, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_6); Py_XDECREF(__pyx_t_7); Py_XDECREF(__pyx_t_8);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_6, __pyx_t_7, __pyx_t_8);
      }
      __pyx_t_6 = __pyx_t_7 = __pyx_t_8 = 0;
    }
    __pyx_pybuffernd_displacement.diminfo[0].strides = __pyx_pybuffernd_displacement.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_displacement.diminfo[0].shape = __pyx_pybuffernd_displacement.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_displacement.diminfo[1].strides = __pyx_pybuffernd_displacement.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_displacement.diminfo[1].shape = __pyx_pybuffernd_displacement.rcbuffer->pybuffer.shape[1];
    if (unlikely(__pyx_t_5 < 0)) __PYX_ERR(0, 139, __pyx_L1_error)
  }
  __pyx_t_4 = 0;
  __pyx_v_displacement = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;

 140: 

+141:     return displacement

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_displacement));
  __pyx_r = ((PyObject *)__pyx_v_displacement);
  goto __pyx_L0;

 142: 

 143: 

 144: cimport cython

 145: @cython.boundscheck(False)  # Deactivate bounds checking

 146: @cython.wraparound(False)   # Deactivate negative indexing.

 147: # Update the position of particles

+148: cdef update_position(np.ndarray[np.float64_t, ndim=2] p, np.ndarray[np.float64_t, ndim=2] delta_r, double min_x=0, double max_x=1):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_update_position(PyArrayObject *__pyx_v_p, PyArrayObject *__pyx_v_delta_r, struct __pyx_opt_args_46_cython_magic_5971c10418fa3492d0bfce342614b144_update_position *__pyx_optional_args) {
  double __pyx_v_min_x = ((double)0.0);
  double __pyx_v_max_x = ((double)1.0);
  PyArrayObject *__pyx_v_new_pos = 0;
  PyObject *__pyx_v_x_out_of_bounds = NULL;
  PyObject *__pyx_v_y_out_of_bounds = NULL;
  Py_ssize_t __pyx_v_i;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_delta_r;
  __Pyx_Buffer __pyx_pybuffer_delta_r;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_new_pos;
  __Pyx_Buffer __pyx_pybuffer_new_pos;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_p;
  __Pyx_Buffer __pyx_pybuffer_p;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("update_position", 0);
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_min_x = __pyx_optional_args->min_x;
      if (__pyx_optional_args->__pyx_n > 1) {
        __pyx_v_max_x = __pyx_optional_args->max_x;
      }
    }
  }
  __pyx_pybuffer_new_pos.pybuffer.buf = NULL;
  __pyx_pybuffer_new_pos.refcount = 0;
  __pyx_pybuffernd_new_pos.data = NULL;
  __pyx_pybuffernd_new_pos.rcbuffer = &__pyx_pybuffer_new_pos;
  __pyx_pybuffer_p.pybuffer.buf = NULL;
  __pyx_pybuffer_p.refcount = 0;
  __pyx_pybuffernd_p.data = NULL;
  __pyx_pybuffernd_p.rcbuffer = &__pyx_pybuffer_p;
  __pyx_pybuffer_delta_r.pybuffer.buf = NULL;
  __pyx_pybuffer_delta_r.refcount = 0;
  __pyx_pybuffernd_delta_r.data = NULL;
  __pyx_pybuffernd_delta_r.rcbuffer = &__pyx_pybuffer_delta_r;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_p.rcbuffer->pybuffer, (PyObject*)__pyx_v_p, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 148, __pyx_L1_error)
  }
  __pyx_pybuffernd_p.diminfo[0].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_p.diminfo[0].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_p.diminfo[1].strides = __pyx_pybuffernd_p.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_p.diminfo[1].shape = __pyx_pybuffernd_p.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_delta_r.rcbuffer->pybuffer, (PyObject*)__pyx_v_delta_r, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 148, __pyx_L1_error)
  }
  __pyx_pybuffernd_delta_r.diminfo[0].strides = __pyx_pybuffernd_delta_r.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_delta_r.diminfo[0].shape = __pyx_pybuffernd_delta_r.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_delta_r.diminfo[1].strides = __pyx_pybuffernd_delta_r.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_delta_r.diminfo[1].shape = __pyx_pybuffernd_delta_r.rcbuffer->pybuffer.shape[1];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_XDECREF(__pyx_t_9);
  __Pyx_XDECREF(__pyx_t_10);
  __Pyx_XDECREF(__pyx_t_11);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_delta_r.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_new_pos.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.update_position", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_delta_r.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_new_pos.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_p.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_new_pos);
  __Pyx_XDECREF(__pyx_v_x_out_of_bounds);
  __Pyx_XDECREF(__pyx_v_y_out_of_bounds);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 149: 

 150:     cdef np.ndarray[np.float64_t, ndim=2] new_pos

 151: 

+152:     new_pos = p + delta_r

  __pyx_t_1 = PyNumber_Add(((PyObject *)__pyx_v_p), ((PyObject *)__pyx_v_delta_r)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 152, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 152, __pyx_L1_error)
  __pyx_t_2 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_new_pos.rcbuffer->pybuffer);
    __pyx_t_3 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_new_pos.rcbuffer->pybuffer, (PyObject*)__pyx_t_2, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_3 < 0)) {
      PyErr_Fetch(&__pyx_t_4, &__pyx_t_5, &__pyx_t_6);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_new_pos.rcbuffer->pybuffer, (PyObject*)__pyx_v_new_pos, &__Pyx_TypeInfo_nn___pyx_t_5numpy_float64_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_4); Py_XDECREF(__pyx_t_5); Py_XDECREF(__pyx_t_6);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_4, __pyx_t_5, __pyx_t_6);
      }
      __pyx_t_4 = __pyx_t_5 = __pyx_t_6 = 0;
    }
    __pyx_pybuffernd_new_pos.diminfo[0].strides = __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_new_pos.diminfo[0].shape = __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_new_pos.diminfo[1].strides = __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_new_pos.diminfo[1].shape = __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.shape[1];
    if (unlikely(__pyx_t_3 < 0)) __PYX_ERR(0, 152, __pyx_L1_error)
  }
  __pyx_t_2 = 0;
  __pyx_v_new_pos = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 153: 

+154:     x_out_of_bounds = np.logical_or(new_pos[:,0] > max_x, new_pos[:,0] < min_x)

  __Pyx_GetModuleGlobalName(__pyx_t_7, __pyx_n_s_np); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __pyx_t_8 = __Pyx_PyObject_GetAttrStr(__pyx_t_7, __pyx_n_s_logical_or); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
/* … */
  __pyx_slice__2 = PySlice_New(Py_None, Py_None, Py_None); if (unlikely(!__pyx_slice__2)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_slice__2);
  __Pyx_GIVEREF(__pyx_slice__2);
  __pyx_t_7 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_new_pos), __pyx_tuple__3); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __pyx_t_9 = PyFloat_FromDouble(__pyx_v_max_x); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_9);
  __pyx_t_10 = PyObject_RichCompare(__pyx_t_7, __pyx_t_9, Py_GT); __Pyx_XGOTREF(__pyx_t_10); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
  __pyx_t_9 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_new_pos), __pyx_tuple__3); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_9);
  __pyx_t_7 = PyFloat_FromDouble(__pyx_v_min_x); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __pyx_t_11 = PyObject_RichCompare(__pyx_t_9, __pyx_t_7, Py_LT); __Pyx_XGOTREF(__pyx_t_11); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  __pyx_t_7 = NULL;
  __pyx_t_3 = 0;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_8))) {
    __pyx_t_7 = PyMethod_GET_SELF(__pyx_t_8);
    if (likely(__pyx_t_7)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_8);
      __Pyx_INCREF(__pyx_t_7);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_8, function);
      __pyx_t_3 = 1;
    }
  }
  #if CYTHON_FAST_PYCALL
  if (PyFunction_Check(__pyx_t_8)) {
    PyObject *__pyx_temp[3] = {__pyx_t_7, __pyx_t_10, __pyx_t_11};
    __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_8, __pyx_temp+1-__pyx_t_3, 2+__pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 154, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
  } else
  #endif
  #if CYTHON_FAST_PYCCALL
  if (__Pyx_PyFastCFunction_Check(__pyx_t_8)) {
    PyObject *__pyx_temp[3] = {__pyx_t_7, __pyx_t_10, __pyx_t_11};
    __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_8, __pyx_temp+1-__pyx_t_3, 2+__pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 154, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
  } else
  #endif
  {
    __pyx_t_9 = PyTuple_New(2+__pyx_t_3); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 154, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_9);
    if (__pyx_t_7) {
      __Pyx_GIVEREF(__pyx_t_7); PyTuple_SET_ITEM(__pyx_t_9, 0, __pyx_t_7); __pyx_t_7 = NULL;
    }
    __Pyx_GIVEREF(__pyx_t_10);
    PyTuple_SET_ITEM(__pyx_t_9, 0+__pyx_t_3, __pyx_t_10);
    __Pyx_GIVEREF(__pyx_t_11);
    PyTuple_SET_ITEM(__pyx_t_9, 1+__pyx_t_3, __pyx_t_11);
    __pyx_t_10 = 0;
    __pyx_t_11 = 0;
    __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_8, __pyx_t_9, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 154, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
  }
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
  __pyx_v_x_out_of_bounds = __pyx_t_1;
  __pyx_t_1 = 0;
  __pyx_tuple__3 = PyTuple_Pack(2, __pyx_slice__2, __pyx_int_0); if (unlikely(!__pyx_tuple__3)) __PYX_ERR(0, 154, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__3);
  __Pyx_GIVEREF(__pyx_tuple__3);

 155: 

+156:     y_out_of_bounds = np.logical_or(new_pos[:,1] > max_x, new_pos[:,1] < min_x)

  __Pyx_GetModuleGlobalName(__pyx_t_8, __pyx_n_s_np); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __pyx_t_9 = __Pyx_PyObject_GetAttrStr(__pyx_t_8, __pyx_n_s_logical_or); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_9);
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
  __pyx_t_8 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_new_pos), __pyx_tuple__4); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __pyx_t_11 = PyFloat_FromDouble(__pyx_v_max_x); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_11);
  __pyx_t_10 = PyObject_RichCompare(__pyx_t_8, __pyx_t_11, Py_GT); __Pyx_XGOTREF(__pyx_t_10); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
  __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
  __pyx_t_11 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_new_pos), __pyx_tuple__4); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_11);
  __pyx_t_8 = PyFloat_FromDouble(__pyx_v_min_x); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __pyx_t_7 = PyObject_RichCompare(__pyx_t_11, __pyx_t_8, Py_LT); __Pyx_XGOTREF(__pyx_t_7); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
  __pyx_t_8 = NULL;
  __pyx_t_3 = 0;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_9))) {
    __pyx_t_8 = PyMethod_GET_SELF(__pyx_t_9);
    if (likely(__pyx_t_8)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_9);
      __Pyx_INCREF(__pyx_t_8);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_9, function);
      __pyx_t_3 = 1;
    }
  }
  #if CYTHON_FAST_PYCALL
  if (PyFunction_Check(__pyx_t_9)) {
    PyObject *__pyx_temp[3] = {__pyx_t_8, __pyx_t_10, __pyx_t_7};
    __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_9, __pyx_temp+1-__pyx_t_3, 2+__pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_8); __pyx_t_8 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  } else
  #endif
  #if CYTHON_FAST_PYCCALL
  if (__Pyx_PyFastCFunction_Check(__pyx_t_9)) {
    PyObject *__pyx_temp[3] = {__pyx_t_8, __pyx_t_10, __pyx_t_7};
    __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_9, __pyx_temp+1-__pyx_t_3, 2+__pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_8); __pyx_t_8 = 0;
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
  } else
  #endif
  {
    __pyx_t_11 = PyTuple_New(2+__pyx_t_3); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 156, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_11);
    if (__pyx_t_8) {
      __Pyx_GIVEREF(__pyx_t_8); PyTuple_SET_ITEM(__pyx_t_11, 0, __pyx_t_8); __pyx_t_8 = NULL;
    }
    __Pyx_GIVEREF(__pyx_t_10);
    PyTuple_SET_ITEM(__pyx_t_11, 0+__pyx_t_3, __pyx_t_10);
    __Pyx_GIVEREF(__pyx_t_7);
    PyTuple_SET_ITEM(__pyx_t_11, 1+__pyx_t_3, __pyx_t_7);
    __pyx_t_10 = 0;
    __pyx_t_7 = 0;
    __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_9, __pyx_t_11, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
  }
  __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
  __pyx_v_y_out_of_bounds = __pyx_t_1;
  __pyx_t_1 = 0;
/* … */
  __pyx_tuple__4 = PyTuple_Pack(2, __pyx_slice__2, __pyx_int_1); if (unlikely(!__pyx_tuple__4)) __PYX_ERR(0, 156, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__4);
  __Pyx_GIVEREF(__pyx_tuple__4);

 157: 

 158: 

+159:     for i in range(len(new_pos)):

  __pyx_t_12 = PyObject_Length(((PyObject *)__pyx_v_new_pos)); if (unlikely(__pyx_t_12 == ((Py_ssize_t)-1))) __PYX_ERR(0, 159, __pyx_L1_error)
  __pyx_t_13 = __pyx_t_12;
  for (__pyx_t_14 = 0; __pyx_t_14 < __pyx_t_13; __pyx_t_14+=1) {
    __pyx_v_i = __pyx_t_14;

+160:         if x_out_of_bounds[i]:

    __pyx_t_1 = __Pyx_GetItemInt(__pyx_v_x_out_of_bounds, __pyx_v_i, Py_ssize_t, 1, PyInt_FromSsize_t, 0, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 160, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __pyx_t_15 = __Pyx_PyObject_IsTrue(__pyx_t_1); if (unlikely(__pyx_t_15 < 0)) __PYX_ERR(0, 160, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    if (__pyx_t_15) {
/* … */
    }

+161:             new_pos[i, 0] = clip(new_pos[i, 0], min_x, max_x)

      __pyx_t_16 = __pyx_v_i;
      __pyx_t_17 = 0;
      __pyx_t_1 = PyFloat_FromDouble((*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_new_pos.diminfo[0].strides, __pyx_t_17, __pyx_pybuffernd_new_pos.diminfo[1].strides))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 161, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __pyx_t_9 = PyFloat_FromDouble(__pyx_v_min_x); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 161, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_9);
      __pyx_t_11 = PyFloat_FromDouble(__pyx_v_max_x); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 161, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_11);
      __pyx_t_17 = __pyx_v_i;
      __pyx_t_16 = 0;
      *__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_new_pos.diminfo[0].strides, __pyx_t_16, __pyx_pybuffernd_new_pos.diminfo[1].strides) = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_clip(__pyx_t_1, __pyx_t_9, __pyx_t_11);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
      __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;

+162:         if y_out_of_bounds[i]:

    __pyx_t_11 = __Pyx_GetItemInt(__pyx_v_y_out_of_bounds, __pyx_v_i, Py_ssize_t, 1, PyInt_FromSsize_t, 0, 0, 0); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 162, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_11);
    __pyx_t_15 = __Pyx_PyObject_IsTrue(__pyx_t_11); if (unlikely(__pyx_t_15 < 0)) __PYX_ERR(0, 162, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
    if (__pyx_t_15) {
/* … */
    }
  }

+163:             new_pos[i, 1] = clip(new_pos[i, 1], min_x, max_x)

      __pyx_t_16 = __pyx_v_i;
      __pyx_t_17 = 1;
      __pyx_t_11 = PyFloat_FromDouble((*__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_new_pos.diminfo[0].strides, __pyx_t_17, __pyx_pybuffernd_new_pos.diminfo[1].strides))); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 163, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_11);
      __pyx_t_9 = PyFloat_FromDouble(__pyx_v_min_x); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 163, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_9);
      __pyx_t_1 = PyFloat_FromDouble(__pyx_v_max_x); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 163, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __pyx_t_17 = __pyx_v_i;
      __pyx_t_16 = 1;
      *__Pyx_BufPtrStrided2d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_new_pos.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_new_pos.diminfo[0].strides, __pyx_t_16, __pyx_pybuffernd_new_pos.diminfo[1].strides) = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_clip(__pyx_t_11, __pyx_t_9, __pyx_t_1);
      __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
      __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 164: 

+165:     return new_pos

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(((PyObject *)__pyx_v_new_pos));
  __pyx_r = ((PyObject *)__pyx_v_new_pos);
  goto __pyx_L0;

 166: 

 167: 

+168: cdef double clip(a, min_value, max_value):

static double __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_clip(PyObject *__pyx_v_a, PyObject *__pyx_v_min_value, PyObject *__pyx_v_max_value) {
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("clip", 0);
/* … */
  /* 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_WriteUnraisable("_cython_magic_5971c10418fa3492d0bfce342614b144.clip", __pyx_clineno, __pyx_lineno, __pyx_filename, 1, 0);
  __pyx_r = 0;
  __pyx_L0:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

+169:     return min(max(a, min_value), max_value)

  __Pyx_INCREF(__pyx_v_max_value);
  __pyx_t_1 = __pyx_v_max_value;
  __Pyx_INCREF(__pyx_v_min_value);
  __pyx_t_2 = __pyx_v_min_value;
  __Pyx_INCREF(__pyx_v_a);
  __pyx_t_3 = __pyx_v_a;
  __pyx_t_5 = PyObject_RichCompare(__pyx_t_2, __pyx_t_3, Py_GT); __Pyx_XGOTREF(__pyx_t_5); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 169, __pyx_L1_error)
  __pyx_t_6 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_6 < 0)) __PYX_ERR(0, 169, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  if (__pyx_t_6) {
    __Pyx_INCREF(__pyx_t_2);
    __pyx_t_4 = __pyx_t_2;
  } else {
    __Pyx_INCREF(__pyx_t_3);
    __pyx_t_4 = __pyx_t_3;
  }
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_INCREF(__pyx_t_4);
  __pyx_t_2 = __pyx_t_4;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_3 = PyObject_RichCompare(__pyx_t_1, __pyx_t_2, Py_LT); __Pyx_XGOTREF(__pyx_t_3); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 169, __pyx_L1_error)
  __pyx_t_6 = __Pyx_PyObject_IsTrue(__pyx_t_3); if (unlikely(__pyx_t_6 < 0)) __PYX_ERR(0, 169, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (__pyx_t_6) {
    __Pyx_INCREF(__pyx_t_1);
    __pyx_t_4 = __pyx_t_1;
  } else {
    __Pyx_INCREF(__pyx_t_2);
    __pyx_t_4 = __pyx_t_2;
  }
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_7 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_7 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 169, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_r = __pyx_t_7;
  goto __pyx_L0;

 170: 

 171: 

 172: # Plot

+173: cdef update_plot(pos,color):

static PyObject *__pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_update_plot(PyObject *__pyx_v_pos, PyObject *__pyx_v_color) {
  PyObject *__pyx_v_xpos = NULL;
  PyObject *__pyx_v_ypos = NULL;
  Py_ssize_t __pyx_v_N;
  PyObject *__pyx_v_N_color = NULL;
  PyObject *__pyx_v_i = NULL;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("update_plot", 0);
/* … */
  /* function exit code */
  __pyx_r = Py_None; __Pyx_INCREF(Py_None);
  goto __pyx_L0;
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_6);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.update_plot", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  __pyx_L0:;
  __Pyx_XDECREF(__pyx_v_xpos);
  __Pyx_XDECREF(__pyx_v_ypos);
  __Pyx_XDECREF(__pyx_v_N_color);
  __Pyx_XDECREF(__pyx_v_i);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 174: 

+175:     plt.clf()

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_plt); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 175, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_clf); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 175, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
    }
  }
  __pyx_t_1 = (__pyx_t_2) ? __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_2) : __Pyx_PyObject_CallNoArg(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 175, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 176: 

+177:     xpos = pos[:, 0]

  __pyx_t_1 = __Pyx_PyObject_GetItem(__pyx_v_pos, __pyx_tuple__3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 177, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_xpos = __pyx_t_1;
  __pyx_t_1 = 0;

 178: 

+179:     ypos = pos[:, 1]

  __pyx_t_1 = __Pyx_PyObject_GetItem(__pyx_v_pos, __pyx_tuple__4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 179, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_ypos = __pyx_t_1;
  __pyx_t_1 = 0;

 180: 

+181:     N = len(pos)

  __pyx_t_4 = PyObject_Length(__pyx_v_pos); if (unlikely(__pyx_t_4 == ((Py_ssize_t)-1))) __PYX_ERR(0, 181, __pyx_L1_error)
  __pyx_v_N = __pyx_t_4;

+182:     N_color = len(color)

  __pyx_t_4 = PyObject_Length(__pyx_v_color); if (unlikely(__pyx_t_4 == ((Py_ssize_t)-1))) __PYX_ERR(0, 182, __pyx_L1_error)
  __pyx_t_1 = PyInt_FromSsize_t(__pyx_t_4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 182, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_N_color = __pyx_t_1;
  __pyx_t_1 = 0;

+183:     for i in range(N):

  __pyx_t_1 = PyInt_FromSsize_t(__pyx_v_N); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 183, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_PyObject_CallOneArg(__pyx_builtin_range, __pyx_t_1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 183, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (likely(PyList_CheckExact(__pyx_t_3)) || PyTuple_CheckExact(__pyx_t_3)) {
    __pyx_t_1 = __pyx_t_3; __Pyx_INCREF(__pyx_t_1); __pyx_t_4 = 0;
    __pyx_t_5 = NULL;
  } else {
    __pyx_t_4 = -1; __pyx_t_1 = PyObject_GetIter(__pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 183, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __pyx_t_5 = Py_TYPE(__pyx_t_1)->tp_iternext; if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 183, __pyx_L1_error)
  }
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  for (;;) {
    if (likely(!__pyx_t_5)) {
      if (likely(PyList_CheckExact(__pyx_t_1))) {
        if (__pyx_t_4 >= PyList_GET_SIZE(__pyx_t_1)) break;
        #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
        __pyx_t_3 = PyList_GET_ITEM(__pyx_t_1, __pyx_t_4); __Pyx_INCREF(__pyx_t_3); __pyx_t_4++; if (unlikely(0 < 0)) __PYX_ERR(0, 183, __pyx_L1_error)
        #else
        __pyx_t_3 = PySequence_ITEM(__pyx_t_1, __pyx_t_4); __pyx_t_4++; if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 183, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        #endif
      } else {
        if (__pyx_t_4 >= PyTuple_GET_SIZE(__pyx_t_1)) break;
        #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
        __pyx_t_3 = PyTuple_GET_ITEM(__pyx_t_1, __pyx_t_4); __Pyx_INCREF(__pyx_t_3); __pyx_t_4++; if (unlikely(0 < 0)) __PYX_ERR(0, 183, __pyx_L1_error)
        #else
        __pyx_t_3 = PySequence_ITEM(__pyx_t_1, __pyx_t_4); __pyx_t_4++; if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 183, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        #endif
      }
    } else {
      __pyx_t_3 = __pyx_t_5(__pyx_t_1);
      if (unlikely(!__pyx_t_3)) {
        PyObject* exc_type = PyErr_Occurred();
        if (exc_type) {
          if (likely(__Pyx_PyErr_GivenExceptionMatches(exc_type, PyExc_StopIteration))) PyErr_Clear();
          else __PYX_ERR(0, 183, __pyx_L1_error)
        }
        break;
      }
      __Pyx_GOTREF(__pyx_t_3);
    }
    __Pyx_XDECREF_SET(__pyx_v_i, __pyx_t_3);
    __pyx_t_3 = 0;
/* … */
  }
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

+184:         plt.plot(xpos[i], ypos[i], "o", color=color[i % N_color])

    __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_plt); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_plot); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_3 = __Pyx_PyObject_GetItem(__pyx_v_xpos, __pyx_v_i); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_6 = __Pyx_PyObject_GetItem(__pyx_v_ypos, __pyx_v_i); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    __pyx_t_7 = PyTuple_New(3); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __Pyx_GIVEREF(__pyx_t_3);
    PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_3);
    __Pyx_GIVEREF(__pyx_t_6);
    PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_6);
    __Pyx_INCREF(__pyx_n_u_o);
    __Pyx_GIVEREF(__pyx_n_u_o);
    PyTuple_SET_ITEM(__pyx_t_7, 2, __pyx_n_u_o);
    __pyx_t_3 = 0;
    __pyx_t_6 = 0;
    __pyx_t_6 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    __pyx_t_3 = PyNumber_Remainder(__pyx_v_i, __pyx_v_N_color); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_8 = __Pyx_PyObject_GetItem(__pyx_v_color, __pyx_t_3); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_8);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    if (PyDict_SetItem(__pyx_t_6, __pyx_n_s_color, __pyx_t_8) < 0) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
    __pyx_t_8 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_7, __pyx_t_6); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_8);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;

 185: 

+186:     plt.xlim(left=-0.1, right=1.1)

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_plt); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 186, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_8 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_xlim); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 186, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyDict_NewPresized(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 186, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_left, __pyx_float_neg_0_1) < 0) __PYX_ERR(0, 186, __pyx_L1_error)
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_right, __pyx_float_1_1) < 0) __PYX_ERR(0, 186, __pyx_L1_error)
  __pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_8, __pyx_empty_tuple, __pyx_t_1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 186, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;

 187: 

+188:     plt.ylim(bottom=-0.1, top=1.1)

  __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_n_s_plt); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 188, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_ylim); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 188, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = __Pyx_PyDict_NewPresized(2); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 188, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  if (PyDict_SetItem(__pyx_t_6, __pyx_n_s_bottom, __pyx_float_neg_0_1) < 0) __PYX_ERR(0, 188, __pyx_L1_error)
  if (PyDict_SetItem(__pyx_t_6, __pyx_n_s_top, __pyx_float_1_1) < 0) __PYX_ERR(0, 188, __pyx_L1_error)
  __pyx_t_8 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_empty_tuple, __pyx_t_6); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 188, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;

 189: 

+190:     plt.grid()

  __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_n_s_plt); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 190, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_grid); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 190, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
    }
  }
  __pyx_t_8 = (__pyx_t_6) ? __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_6) : __Pyx_PyObject_CallNoArg(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 190, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;

 191: 

+192:     plt.draw()

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_plt); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 192, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_draw); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 192, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
    __pyx_t_1 = PyMethod_GET_SELF(__pyx_t_6);
    if (likely(__pyx_t_1)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
      __Pyx_INCREF(__pyx_t_1);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_6, function);
    }
  }
  __pyx_t_8 = (__pyx_t_1) ? __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_1) : __Pyx_PyObject_CallNoArg(__pyx_t_6);
  __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 192, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;

 193: 

+194:     plt.pause(0.0001)

  __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_n_s_plt); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 194, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_pause); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 194, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
    }
  }
  __pyx_t_8 = (__pyx_t_6) ? __Pyx_PyObject_Call2Args(__pyx_t_1, __pyx_t_6, __pyx_float_0_0001) : __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_float_0_0001);
  __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 194, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;

 195: 

+196: compute_time_cy = timeit.timeit(lambda: simulate_cy(n=10, time_step=1000, show_plot=False), number=1)

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_2lambda(PyObject *__pyx_self, CYTHON_UNUSED PyObject *unused); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_5971c10418fa3492d0bfce342614b144_2lambda = {"lambda", (PyCFunction)__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_2lambda, METH_NOARGS, 0};
static PyObject *__pyx_pw_46_cython_magic_5971c10418fa3492d0bfce342614b144_2lambda(PyObject *__pyx_self, CYTHON_UNUSED PyObject *unused) {
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("lambda (wrapper)", 0);
  __pyx_r = __pyx_lambda_funcdef_46_cython_magic_5971c10418fa3492d0bfce342614b144_lambda(__pyx_self);

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_lambda_funcdef_46_cython_magic_5971c10418fa3492d0bfce342614b144_lambda(CYTHON_UNUSED PyObject *__pyx_self) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("lambda", 0);
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2.__pyx_n = 1;
  __pyx_t_2.show_plot = Py_False;
  __pyx_t_1 = __pyx_f_46_cython_magic_5971c10418fa3492d0bfce342614b144_simulate_cy(10, 0x3E8, &__pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_5971c10418fa3492d0bfce342614b144.lambda", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_timeit); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_timeit); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_CyFunction_New(&__pyx_mdef_46_cython_magic_5971c10418fa3492d0bfce342614b144_2lambda, 0, __pyx_n_s_lambda, NULL, __pyx_n_s_cython_magic_5971c10418fa3492d0, __pyx_d, NULL); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_GIVEREF(__pyx_t_2);
  PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_2);
  __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_number, __pyx_int_1) < 0) __PYX_ERR(0, 196, __pyx_L1_error)
  __pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_compute_time_cy, __pyx_t_4) < 0) __PYX_ERR(0, 196, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 197: 

+198: print("simulate_cy execution time:", compute_time_cy)

  __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_compute_time_cy); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 198, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 198, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_INCREF(__pyx_kp_u_simulate_cy_execution_time);
  __Pyx_GIVEREF(__pyx_kp_u_simulate_cy_execution_time);
  PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_kp_u_simulate_cy_execution_time);
  __Pyx_GIVEREF(__pyx_t_4);
  PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_4);
  __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_t_2, NULL); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 198, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

8.3. Compare

from IPython.display import HTML
import pandas as pd

data = {
    'Methods': [ 'Python code','Cython code'],
    'Speed(s)': [compute_time_py, compute_time_cy],
    'Percentage(%)': [100, compute_time_py/compute_time_cy*100]
}
df = pd.DataFrame(data)

# Creating style functions
def add_border(val):
    return 'border: 1px solid black'

# Applying style functions to data boxes
styled_df = df.style.applymap(add_border)

# Defining CSS styles
table_style = [
    {'selector': 'table', 'props': [('border-collapse', 'collapse')]},
    {'selector': 'th, td', 'props': [('border', '1px solid black')]}
]

# Adding styles to stylised data boxes
styled_df.set_table_styles(table_style)

# Displaying stylised data boxes in Jupyter Notebook
HTML(styled_df.to_html())

	Methods	Speed(s)	Percentage(%)
0	Python code	1.069363	100.000000
1	Cython code	0.426124	250.951049

8.4. random particles

def random(n):
    
    np.random.seed(0)
    
    p = np.random.rand(n, 2)
    
    return p

8.4.1. main function

8.4.1.1. add types

Since all the functions introduced in the function are of numpy type, the arrays here are defined as np.ndarray[np.float64_t, ndim=2], which represents a 2-dimensional array with elements of np.float64_t.

8.4.1.2. boundscheck and wraparound

There is a for loop, so close boundscheck and wraparound.

8.4.1.3. Note:

np.ndarray[np.float64_t, ndim=2] is a Cython type for numpy arrays. Its definition contains information about the data type and the dimension of the array, allowing functions and methods from the numpy library to be conveniently called. Arrays of type np.ndarray can be manipulated using a variety of numpy methods and functions, providing a high degree of flexibility and versatility.

In contrast, double[:, :] is a Cython-specific type that represents only a two-dimensional array, with no information about the specified data type. Using arrays of type double[:, :] may not allow you to use all numpy functions and methods due to the lack of information about the data type, but its access is faster and can improve the execution of your code to some extent.

Therefore, arrays of type np.ndarray should be used if functions and methods from the numpy library need to be called, or if numerical calculations need to be performed, while arrays of type double[:, :] may be considered if higher execution efficiency is sought.

However,

In Cython, for multi-dimensional arrays, you must use the NumPy array type np.ndarray to perform operations such as division, otherwise an error will occur

cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
cpdef simulate_cy(int n, int time_step, show_plot=False):
    
    cdef np.ndarray[np.float64_t, ndim=2] p

    p = random(n)
    print("initial_P:\n",p)

    cdef int i
    cdef np.ndarray[np.float64_t, ndim=2] total_force, x_det, pos
    
    for i in range(time_step):
        
        total_force = combined_force(p, n)
    
        x_det = displacement(total_force, delta_t=0.0001)
        
        # update the particles
        p = update_position(p, x_det)
        
        pos = p
        
        # plot particles
        colors = ['red', 'green', 'blue', 'orange'] 
        if show_plot:

            if i % 2 == 0:

                update_plot(pos,colors)

#                print("in iteration", i)
            
    # plot finally result
    print("P({}): \n".format(time_step), p)
    
    return p

8.4.2. Calculate the strength of the repulsion

8.4.2.1. add types

Since all the functions introduced in the function are of numpy type, the arrays here are defined as np.ndarray[np.float64_t, ndim=2], which represents a 2-dimensional array with elements of np.float64_t.

8.4.2.2. change np.linalg.norm(r)

Replace the original np.linalg.norm(r) with math.sqrt(square_sum(r)) because np.linalg.norm(r) contains the for loop

8.4.2.3. boundscheck and wraparound

There is a for loop, so close boundscheck and wraparound.

8.4.2.4. Using multiple threads

from cython.parallel import prange

# Calculate the strength of the repulsion
cpdef f(np.ndarray[np.float64_t, ndim=1]r, double c1=1, double c2=1):

    # Calculate the force
    
    cdef np.ndarray[np.float64_t, ndim=1] abs_r
    
    abs_r = abs(r)
    
    cdef double mag
    
    mag = c1 / math.sqrt(square_sum(r)) 

    return mag


@cython.boundscheck(False)
@cython.wraparound(False)
#calculate square
cpdef square_sum(np.ndarray[np.float64_t, ndim=1] arr):
    cdef Py_ssize_t i, n
    cdef double s = 0
    n = arr.shape[0]
    for i in prange(n, nogil=True, num_threads=1):
        s += arr[i] * arr[i]
#    print(type(s))
    return s

8.4.3. Calculate the total force for each particle

8.4.3.1. add types

Since all the functions introduced in the function are of numpy type, the arrays here are defined as np.ndarray[np.float64_t, ndim=2], which represents a 2-dimensional array with elements of np.float64_t.

For i, j …… in a for loop, we define them to be of type Py_ssize_t

8.4.3.2. change np.zeros_like()

8.4.3.3. boundscheck and wraparound

There is a for loop, so close boundscheck and wraparound.

8.4.3.4. can not use multiple threads

Because of the repeated calls to numpy in the for loop, gil is needed here and you cannot use multiple threads here

cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
# Calculate the total force for each particle
cdef combined_force(np.ndarray[np.float64_t, ndim=2] p, n):
    
    cdef np.ndarray[np.float64_t, ndim=2] total_force
    
#    total_force = np.zeros_like(p)
    total_force = np.empty_like(p)

    cdef Py_ssize_t k, l
    for k in prange(total_force.shape[0], nogil=True, num_threads=4):
        for l in range(total_force.shape[1]):
            total_force[k, l] = 0
    
    cdef np.ndarray[np.float64_t, ndim=1] fn_sum
    
    cdef np.ndarray[np.float64_t, ndim=1] r
    
    cdef np.ndarray[np.float64_t, ndim=1] fn
 
    
    cdef Py_ssize_t i, j
    
    for i in range(n):
        
        fn_sum = np.zeros(2)
        
        for j in range(n):
            
            if j != i:
                
                r = p[j] - p[i]

                fn =  -1 * f(r) * np.sign(r) 

                fn_sum += fn
           
            total_force[i] = fn_sum
            
    return total_force

8.4.4. Calculate the displacement between two particles

8.4.4.1. add types

Since all the functions introduced in the function are of numpy type, the arrays here are defined as np.ndarray[np.float64_t, ndim=2], which represents a 2-dimensional array with elements of np.float64_t.

8.4.4.2. can not use memoryviews

Because there is division here, it is not possible to compute double[:, :] / double[:, :] in cython, but only as a numpy type

# Calculate the displacement between two particles
cpdef displacement(np.ndarray[np.float64_t, ndim=2] total_force, double eta=1, double delta_t=1):
    cdef np.ndarray[np.float64_t, ndim=2] displacement

    displacement = total_force / eta * delta_t

    return displacement

8.5. Update the position of particles

8.5.1. add types

8.5.2. np.clip

Changr ‘new_pos[x_out_of_bounds, 0] = np.clip(new_pos[x_out_of_bounds, 0], min_x, max_x)’

8.5.3. can not use memoryviews

because np.logical_or()

cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
# Update the position of particles
cdef update_position(np.ndarray[np.float64_t, ndim=2] p, np.ndarray[np.float64_t, ndim=2] delta_r, double min_x=0, double max_x=1):
    
    cdef np.ndarray[np.float64_t, ndim=2] new_pos
    
    new_pos = p + delta_r
    
    x_out_of_bounds = np.logical_or(new_pos[:,0] > max_x, new_pos[:,0] < min_x)
    
    y_out_of_bounds = np.logical_or(new_pos[:,1] > max_x, new_pos[:,1] < min_x)
    

    for i in range(len(new_pos)):
        if x_out_of_bounds[i]:
            new_pos[i, 0] = clip(new_pos[i, 0], min_x, max_x)
        if y_out_of_bounds[i]:
            new_pos[i, 1] = clip(new_pos[i, 1], min_x, max_x)
    
    return new_pos


cdef double clip(a, min_value, max_value):
    return min(max(a, min_value), max_value)