Spatial evolutionary game simulator¶

S. Kolotev, A. Malyutin, E. Burovski, S. Krashakov and L. Shchur, Dynamic fractals in spatial evolutionary games, Physica A 499, 142 (2018)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import animation

%matplotlib notebook

In [2]:
import sys
import cython

print("python ", sys.version)
print("numpy", np.__version__, ", matplotlib", mpl.__version__, ", and Cython", cython.__version__)

('python ', '2.7.13 |Anaconda 4.4.0 (64-bit)| (default, Dec 20 2016, 23:09:15) \n[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]')
('numpy', '1.12.1', ', matplotlib', '2.0.2', ', and Cython', '0.25.2')


Main update function:¶

In [3]:
def evolve2(field, b, num_steps=1):
L = field.shape[0]
current = np.zeros((L, L), dtype=int)
scores = np.zeros((L, L), dtype=float)

for step in range(num_steps):
current = field.copy()
scores = np.zeros((L, L), dtype=float)

for x in range(L):
for y in range(L):
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
scores[x, y] += (1 - field[ix, jy])

if field[x, y] == 1:
scores[x, y] *= b

for x in range(L):
for y in range(L):
bestX = x
bestY = y
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
if (scores[bestX, bestY] < scores[ix, jy]):
bestX = ix
bestY = jy

field[x, y] = current[bestX, bestY]
return field


Run the simulation and animate¶

In [4]:
L = 75
field = np.zeros((L, L), dtype=int)
field[L//2, L//2] = 1

# draw the initial field
fig = plt.figure()
im = plt.imshow(field, animated=True, cmap=plt.get_cmap('Paired'))

# updater function
def updatefig(*args):
global field
field = evolve2(field, 1.81, 1)
im.set_array(field)
return im,

# animate!
anim = animation.FuncAnimation(fig, updatefig, frames=200, interval=50, blit=True)

plt.show()


Timings¶

In [5]:
L = 44
field = np.zeros((L, L), dtype=int)
field[L//2, L//2] = 1

%timeit evolve2(field, 1.81, 10)

1 loop, best of 3: 345 ms per loop


Enter Cython¶

In class, we started with the original code and did the conversion iteratively, removing yellow lines from the annotated output. Below is only the final result.

In [6]:
%load_ext cython

In [7]:
%%cython -a

import numpy as np

import cython

@cython.cdivision(True)
@cython.boundscheck(False)
def evolve2_1(long[:, ::1] field, double b, int num_steps=1):

cdef int x, y, L, i, j, ix, jy, step

L = field.shape[0]
cdef double[:, ::1] scores = np.zeros((L, L), dtype=float)

cdef double[:, ::1] _zeros = np.zeros((L, L), dtype=float)
cdef long[:, ::1] current = field.copy()

for step in range(num_steps):
current = field.copy()
scores[...] = _zeros

for x in range(L):
for y in range(L):
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
scores[x, y] += (1 - field[ix, jy])

if field[x, y] == 1:
scores[x, y] *= b

for x in range(L):
for y in range(L):
bestX = x
bestY = y
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
if (scores[bestX, bestY] < scores[ix, jy]):
bestX = ix
bestY = jy

field[x, y] = current[bestX, bestY]
return field

Out[7]:
Cython: _cython_magic_4a51e9818ff378ac421ff64c016d5995.pyx

Generated by Cython 0.25.2

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

 01:
+02: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
__pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 03:
 04: import cython
 05:
 06: @cython.cdivision(True)
 07: @cython.boundscheck(False)
+08: def evolve2_1(long[:, ::1] field, double b, int num_steps=1):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_1evolve2_1(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_1evolve2_1 = {"evolve2_1", (PyCFunction)__pyx_pw_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_1evolve2_1, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_1evolve2_1(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
__Pyx_memviewslice __pyx_v_field = { 0, 0, { 0 }, { 0 }, { 0 } };
double __pyx_v_b;
int __pyx_v_num_steps;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("evolve2_1 (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_field,&__pyx_n_s_b,&__pyx_n_s_num_steps,0};
PyObject* values[3] = {0,0,0};
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case  0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case  0:
if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_field)) != 0)) kw_args--;
else goto __pyx_L5_argtuple_error;
case  1:
if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("evolve2_1", 0, 2, 3, 1); __PYX_ERR(0, 8, __pyx_L3_error)
}
case  2:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_num_steps);
if (value) { values[2] = value; kw_args--; }
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "evolve2_1") < 0)) __PYX_ERR(0, 8, __pyx_L3_error)
}
} else {
switch (PyTuple_GET_SIZE(__pyx_args)) {
case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
break;
default: goto __pyx_L5_argtuple_error;
}
}
__pyx_v_field = __Pyx_PyObject_to_MemoryviewSlice_d_dc_long(values[0]); if (unlikely(!__pyx_v_field.memview)) __PYX_ERR(0, 8, __pyx_L3_error)
__pyx_v_b = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_b == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
if (values[2]) {
__pyx_v_num_steps = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_num_steps == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 8, __pyx_L3_error)
} else {
__pyx_v_num_steps = ((int)1);
}
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("evolve2_1", 0, 2, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 8, __pyx_L3_error)
__pyx_L3_error:;
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_evolve2_1(__pyx_self, __pyx_v_field, __pyx_v_b, __pyx_v_num_steps);

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

static PyObject *__pyx_pf_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_evolve2_1(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_field, double __pyx_v_b, int __pyx_v_num_steps) {
int __pyx_v_x;
int __pyx_v_y;
int __pyx_v_L;
int __pyx_v_i;
int __pyx_v_j;
int __pyx_v_ix;
int __pyx_v_jy;
CYTHON_UNUSED int __pyx_v_step;
__Pyx_memviewslice __pyx_v_scores = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v__zeros = { 0, 0, { 0 }, { 0 }, { 0 } };
__Pyx_memviewslice __pyx_v_current = { 0, 0, { 0 }, { 0 }, { 0 } };
int __pyx_v_bestX;
int __pyx_v_bestY;
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("evolve2_1", 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_XDEC_MEMVIEW(&__pyx_t_5, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
__pyx_r = NULL;
__pyx_L0:;
__PYX_XDEC_MEMVIEW(&__pyx_v_scores, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v__zeros, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_current, 1);
__PYX_XDEC_MEMVIEW(&__pyx_v_field, 1);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__14 = PyTuple_Pack(16, __pyx_n_s_field, __pyx_n_s_b, __pyx_n_s_num_steps, __pyx_n_s_x, __pyx_n_s_y, __pyx_n_s_L, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_ix, __pyx_n_s_jy, __pyx_n_s_step, __pyx_n_s_scores, __pyx_n_s_zeros_2, __pyx_n_s_current, __pyx_n_s_bestX, __pyx_n_s_bestY); if (unlikely(!__pyx_tuple__14)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__14);
__Pyx_GIVEREF(__pyx_tuple__14);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_4a51e9818ff378ac421ff64c016d5995_1evolve2_1, NULL, __pyx_n_s_cython_magic_4a51e9818ff378ac42); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_evolve2_1, __pyx_t_1) < 0) __PYX_ERR(0, 8, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_codeobj__15 = (PyObject*)__Pyx_PyCode_New(3, 0, 16, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__14, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_br_cache_ipython_cython__c, __pyx_n_s_evolve2_1, 8, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__15)) __PYX_ERR(0, 8, __pyx_L1_error)

 09:
 10:     cdef int x, y, L, i, j, ix, jy, step
 11:
+12:     L = field.shape[0]
  __pyx_v_L = (__pyx_v_field.shape[0]);

+13:     cdef double[:, ::1] scores = np.zeros((L, L), dtype=float)
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_L); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_L); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_3);
__pyx_t_1 = 0;
__pyx_t_3 = 0;
__pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_GIVEREF(__pyx_t_4);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_4);
__pyx_t_4 = 0;
__pyx_t_4 = PyDict_New(); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
if (PyDict_SetItem(__pyx_t_4, __pyx_n_s_dtype, ((PyObject *)(&PyFloat_Type))) < 0) __PYX_ERR(0, 13, __pyx_L1_error)
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_3, __pyx_t_4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_5 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(__pyx_t_1);
if (unlikely(!__pyx_t_5.memview)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_scores = __pyx_t_5;
__pyx_t_5.memview = NULL;
__pyx_t_5.data = NULL;

 14:
+15:     cdef double[:, ::1] _zeros = np.zeros((L, L), dtype=float)
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_L); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_L); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_3);
__pyx_t_1 = 0;
__pyx_t_3 = 0;
__pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 15, __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 = PyDict_New(); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_dtype, ((PyObject *)(&PyFloat_Type))) < 0) __PYX_ERR(0, 15, __pyx_L1_error)
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_3, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_5 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(__pyx_t_1);
if (unlikely(!__pyx_t_5.memview)) __PYX_ERR(0, 15, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v__zeros = __pyx_t_5;
__pyx_t_5.memview = NULL;
__pyx_t_5.data = NULL;

+16:     cdef long[:, ::1] current = field.copy()
  __pyx_t_6 = __pyx_memoryview_copy_slice_d_dc_long_c(__pyx_v_field); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 16, __pyx_L1_error)
__pyx_v_current = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;

 17:
+18:     for step in range(num_steps):
  __pyx_t_7 = __pyx_v_num_steps;
for (__pyx_t_8 = 0; __pyx_t_8 < __pyx_t_7; __pyx_t_8+=1) {
__pyx_v_step = __pyx_t_8;

+19:         current = field.copy()
    __pyx_t_6 = __pyx_memoryview_copy_slice_d_dc_long_c(__pyx_v_field); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 19, __pyx_L1_error)
__PYX_XDEC_MEMVIEW(&__pyx_v_current, 1);
__pyx_v_current = __pyx_t_6;
__pyx_t_6.memview = NULL;
__pyx_t_6.data = NULL;

+20:         scores[...] = _zeros
    if (unlikely(__pyx_memoryview_copy_contents(__pyx_v__zeros, __pyx_v_scores, 2, 2, 0) < 0)) __PYX_ERR(0, 20, __pyx_L1_error)

 21:
+22:         for x in range(L):
    __pyx_t_9 = __pyx_v_L;
for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
__pyx_v_x = __pyx_t_10;

+23:             for y in range(L):
      __pyx_t_11 = __pyx_v_L;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_y = __pyx_t_12;

+24:                 for i in range(-1, 2):
        for (__pyx_t_13 = -1L; __pyx_t_13 < 2; __pyx_t_13+=1) {
__pyx_v_i = __pyx_t_13;

+25:                     for j in range(-1, 2):
          for (__pyx_t_14 = -1L; __pyx_t_14 < 2; __pyx_t_14+=1) {
__pyx_v_j = __pyx_t_14;

+26:                         ix = (x + i) % L
            __pyx_v_ix = ((__pyx_v_x + __pyx_v_i) % __pyx_v_L);

+27:                         jy = (y + j) % L
            __pyx_v_jy = ((__pyx_v_y + __pyx_v_j) % __pyx_v_L);

+28:                         scores[x, y] += (1 - field[ix, jy])
            __pyx_t_15 = __pyx_v_ix;
__pyx_t_16 = __pyx_v_jy;
if (__pyx_t_15 < 0) __pyx_t_15 += __pyx_v_field.shape[0];
if (__pyx_t_16 < 0) __pyx_t_16 += __pyx_v_field.shape[1];
__pyx_t_17 = __pyx_v_x;
__pyx_t_18 = __pyx_v_y;
if (__pyx_t_17 < 0) __pyx_t_17 += __pyx_v_scores.shape[0];
if (__pyx_t_18 < 0) __pyx_t_18 += __pyx_v_scores.shape[1];
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_17 * __pyx_v_scores.strides[0]) )) + __pyx_t_18)) )) += (1 - (*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_15 * __pyx_v_field.strides[0]) )) + __pyx_t_16)) ))));
}
}

 29:
+30:                 if field[x, y] == 1:
        __pyx_t_19 = __pyx_v_x;
__pyx_t_20 = __pyx_v_y;
if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_v_field.shape[0];
if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_v_field.shape[1];
__pyx_t_21 = (((*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_19 * __pyx_v_field.strides[0]) )) + __pyx_t_20)) ))) == 1) != 0);
if (__pyx_t_21) {
/* … */
}
}
}

+31:                     scores[x, y] *= b
          __pyx_t_22 = __pyx_v_x;
__pyx_t_23 = __pyx_v_y;
if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_v_scores.shape[0];
if (__pyx_t_23 < 0) __pyx_t_23 += __pyx_v_scores.shape[1];
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_22 * __pyx_v_scores.strides[0]) )) + __pyx_t_23)) )) *= __pyx_v_b;

 32:
+33:         for x in range(L):
    __pyx_t_9 = __pyx_v_L;
for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
__pyx_v_x = __pyx_t_10;

+34:             for y in range(L):
      __pyx_t_11 = __pyx_v_L;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_y = __pyx_t_12;

+35:                 bestX = x
        __pyx_v_bestX = __pyx_v_x;

+36:                 bestY = y
        __pyx_v_bestY = __pyx_v_y;

+37:                 for i in range(-1, 2):
        for (__pyx_t_13 = -1L; __pyx_t_13 < 2; __pyx_t_13+=1) {
__pyx_v_i = __pyx_t_13;

+38:                     for j in range(-1, 2):
          for (__pyx_t_14 = -1L; __pyx_t_14 < 2; __pyx_t_14+=1) {
__pyx_v_j = __pyx_t_14;

+39:                         ix = (x + i) % L
            __pyx_v_ix = ((__pyx_v_x + __pyx_v_i) % __pyx_v_L);

+40:                         jy = (y + j) % L
            __pyx_v_jy = ((__pyx_v_y + __pyx_v_j) % __pyx_v_L);

+41:                         if (scores[bestX, bestY] < scores[ix, jy]):
            __pyx_t_24 = __pyx_v_bestX;
__pyx_t_25 = __pyx_v_bestY;
if (__pyx_t_24 < 0) __pyx_t_24 += __pyx_v_scores.shape[0];
if (__pyx_t_25 < 0) __pyx_t_25 += __pyx_v_scores.shape[1];
__pyx_t_26 = __pyx_v_ix;
__pyx_t_27 = __pyx_v_jy;
if (__pyx_t_26 < 0) __pyx_t_26 += __pyx_v_scores.shape[0];
if (__pyx_t_27 < 0) __pyx_t_27 += __pyx_v_scores.shape[1];
__pyx_t_21 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_24 * __pyx_v_scores.strides[0]) )) + __pyx_t_25)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_26 * __pyx_v_scores.strides[0]) )) + __pyx_t_27)) )))) != 0);
if (__pyx_t_21) {
/* … */
}
}
}

+42:                             bestX = ix
              __pyx_v_bestX = __pyx_v_ix;

+43:                             bestY = jy
              __pyx_v_bestY = __pyx_v_jy;

 44:
+45:                 field[x, y] = current[bestX, bestY]
        __pyx_t_28 = __pyx_v_bestX;
__pyx_t_29 = __pyx_v_bestY;
if (__pyx_t_28 < 0) __pyx_t_28 += __pyx_v_current.shape[0];
if (__pyx_t_29 < 0) __pyx_t_29 += __pyx_v_current.shape[1];
__pyx_t_30 = __pyx_v_x;
__pyx_t_31 = __pyx_v_y;
if (__pyx_t_30 < 0) __pyx_t_30 += __pyx_v_field.shape[0];
if (__pyx_t_31 < 0) __pyx_t_31 += __pyx_v_field.shape[1];
*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_30 * __pyx_v_field.strides[0]) )) + __pyx_t_31)) )) = (*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_current.data + __pyx_t_28 * __pyx_v_current.strides[0]) )) + __pyx_t_29)) )));
}
}
}

+46:     return field
  __Pyx_XDECREF(__pyx_r);
__pyx_t_1 = __pyx_memoryview_fromslice(__pyx_v_field, 2, (PyObject *(*)(char *)) __pyx_memview_get_long, (int (*)(char *, PyObject *)) __pyx_memview_set_long, 0);; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_r = __pyx_t_1;
__pyx_t_1 = 0;
goto __pyx_L0;

In [10]:
L = 44
field = np.zeros((L, L), dtype=int)
field[L//2, L//2] = 1

%timeit evolve2_1(field, 1.81, 10)

100 loops, best of 3: 1.98 ms per loop

In [11]:
345 / 1.98

Out[11]:
174.24242424242425

The recommended code structure is to have a python-facing def function which takes care of the input and output, allocated memory etc, and delegates all heavy lifting to an internal cdef function.

In [ ]:
%%cython -a

import numpy as np
import cython

def evolve2_2(field, b, num_steps=1):
# validate input here
field = np.atleast_2d(field)
# XXX check that field is integer dtype etc

# allocate memory here
L = field.shape[0]
scores = np.zeros((L, L), dtype=float)
current = field.copy()
_zeros = np.zeros((L, L), dtype=float)

cdef:
int _num_steps = num_steps
double _b = b

# do the work (field array is updated in-place)
_evolve2_2_impl(field, _b, _num_steps, scores, _zeros, current)

# convert memoryviews to numpy arrays
return np.asarray(field)

@cython.cdivision(True)
@cython.boundscheck(False)
cdef void _evolve2_2_impl(long[:, ::1] field, double b, int num_steps,
double[:, ::1] scores, double[:, ::1] _zeros,
long[:, ::1] current) nogil:

cdef int x, y, L, i, j, ix, jy, step

L = field.shape[0]

for step in range(num_steps):
current[...] = field
scores[...] = _zeros

for x in range(L):
for y in range(L):
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
scores[x, y] += (1 - field[ix, jy])

if field[x, y] == 1:
scores[x, y] *= b

for x in range(L):
for y in range(L):
bestX = x
bestY = y
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
if (scores[bestX, bestY] < scores[ix, jy]):
bestX = ix
bestY = jy

field[x, y] = current[bestX, bestY]

Out[ ]:
Cython: _cython_magic_451bc548913a9407f1c439613b3672aa.pyx

Generated by Cython 0.25.2

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

 01:
+02: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
__pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 03: import cython
 04:
+05: def evolve2_2(field, b, num_steps=1):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_451bc548913a9407f1c439613b3672aa_1evolve2_2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_451bc548913a9407f1c439613b3672aa_1evolve2_2 = {"evolve2_2", (PyCFunction)__pyx_pw_46_cython_magic_451bc548913a9407f1c439613b3672aa_1evolve2_2, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_451bc548913a9407f1c439613b3672aa_1evolve2_2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
PyObject *__pyx_v_field = 0;
PyObject *__pyx_v_b = 0;
PyObject *__pyx_v_num_steps = 0;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("evolve2_2 (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_field,&__pyx_n_s_b,&__pyx_n_s_num_steps,0};
PyObject* values[3] = {0,0,0};
values[2] = ((PyObject *)__pyx_int_1);
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case  0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case  0:
if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_field)) != 0)) kw_args--;
else goto __pyx_L5_argtuple_error;
case  1:
if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
else {
__Pyx_RaiseArgtupleInvalid("evolve2_2", 0, 2, 3, 1); __PYX_ERR(0, 5, __pyx_L3_error)
}
case  2:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_num_steps);
if (value) { values[2] = value; kw_args--; }
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "evolve2_2") < 0)) __PYX_ERR(0, 5, __pyx_L3_error)
}
} else {
switch (PyTuple_GET_SIZE(__pyx_args)) {
case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
break;
default: goto __pyx_L5_argtuple_error;
}
}
__pyx_v_field = values[0];
__pyx_v_b = values[1];
__pyx_v_num_steps = values[2];
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("evolve2_2", 0, 2, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 5, __pyx_L3_error)
__pyx_L3_error:;
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_46_cython_magic_451bc548913a9407f1c439613b3672aa_evolve2_2(__pyx_self, __pyx_v_field, __pyx_v_b, __pyx_v_num_steps);

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

static PyObject *__pyx_pf_46_cython_magic_451bc548913a9407f1c439613b3672aa_evolve2_2(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_field, PyObject *__pyx_v_b, PyObject *__pyx_v_num_steps) {
PyObject *__pyx_v_L = NULL;
PyObject *__pyx_v_scores = NULL;
PyObject *__pyx_v_current = NULL;
PyObject *__pyx_v__zeros = NULL;
int __pyx_v__num_steps;
double __pyx_v__b;
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("evolve2_2", 0);
__Pyx_INCREF(__pyx_v_field);
/* … */
/* 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_XDEC_MEMVIEW(&__pyx_t_7, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_8, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_9, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_10, 1);
__pyx_r = NULL;
__pyx_L0:;
__Pyx_XDECREF(__pyx_v_L);
__Pyx_XDECREF(__pyx_v_scores);
__Pyx_XDECREF(__pyx_v_current);
__Pyx_XDECREF(__pyx_v__zeros);
__Pyx_XDECREF(__pyx_v_field);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__14 = PyTuple_Pack(9, __pyx_n_s_field, __pyx_n_s_b, __pyx_n_s_num_steps, __pyx_n_s_L, __pyx_n_s_scores, __pyx_n_s_current, __pyx_n_s_zeros_2, __pyx_n_s_num_steps_2, __pyx_n_s_b_2); if (unlikely(!__pyx_tuple__14)) __PYX_ERR(0, 5, __pyx_L1_error)
__Pyx_GOTREF(__pyx_tuple__14);
__Pyx_GIVEREF(__pyx_tuple__14);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_451bc548913a9407f1c439613b3672aa_1evolve2_2, NULL, __pyx_n_s_cython_magic_451bc548913a9407f1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 5, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_evolve2_2, __pyx_t_1) < 0) __PYX_ERR(0, 5, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_codeobj__15 = (PyObject*)__Pyx_PyCode_New(3, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__14, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_br_cache_ipython_cython__c, __pyx_n_s_evolve2_2, 5, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__15)) __PYX_ERR(0, 5, __pyx_L1_error)

 06:     # validate input here
+07:     field = np.atleast_2d(field)
  __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_atleast_2d); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 7, __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);
}
}
if (!__pyx_t_2) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_v_field); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_v_field};
__pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_3)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_v_field};
__pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else
#endif
{
__pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __pyx_t_2 = NULL;
__Pyx_INCREF(__pyx_v_field);
__Pyx_GIVEREF(__pyx_v_field);
PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_v_field);
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 7, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
}
}
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_DECREF_SET(__pyx_v_field, __pyx_t_1);
__pyx_t_1 = 0;

 08:     # XXX check that field is integer dtype etc
 09:
 10:     # allocate memory here
+11:     L = field.shape[0]
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_v_field, __pyx_n_s_shape); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_3 = __Pyx_GetItemInt(__pyx_t_1, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 11, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_L = __pyx_t_3;
__pyx_t_3 = 0;

+12:     scores = np.zeros((L, L), dtype=float)
  __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_zeros); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_INCREF(__pyx_v_L);
__Pyx_GIVEREF(__pyx_v_L);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_v_L);
__Pyx_INCREF(__pyx_v_L);
__Pyx_GIVEREF(__pyx_v_L);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_L);
__pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_3);
__pyx_t_3 = 0;
__pyx_t_3 = PyDict_New(); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, ((PyObject *)(&PyFloat_Type))) < 0) __PYX_ERR(0, 12, __pyx_L1_error)
__pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 12, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_v_scores = __pyx_t_2;
__pyx_t_2 = 0;

+13:     current = field.copy()
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_v_field, __pyx_n_s_copy); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__pyx_t_4 = NULL;
if (CYTHON_UNPACK_METHODS && likely(PyMethod_Check(__pyx_t_3))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
if (likely(__pyx_t_4)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_3, function);
}
}
if (__pyx_t_4) {
__pyx_t_2 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_4); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 13, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
} else {
__pyx_t_2 = __Pyx_PyObject_CallNoArg(__pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 13, __pyx_L1_error)
}
__Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_v_current = __pyx_t_2;
__pyx_t_2 = 0;

+14:     _zeros = np.zeros((L, L), dtype=float)
  __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__Pyx_INCREF(__pyx_v_L);
__Pyx_GIVEREF(__pyx_v_L);
PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_v_L);
__Pyx_INCREF(__pyx_v_L);
__Pyx_GIVEREF(__pyx_v_L);
PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_v_L);
__pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_GIVEREF(__pyx_t_2);
PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2);
__pyx_t_2 = 0;
__pyx_t_2 = PyDict_New(); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_dtype, ((PyObject *)(&PyFloat_Type))) < 0) __PYX_ERR(0, 14, __pyx_L1_error)
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_v__zeros = __pyx_t_1;
__pyx_t_1 = 0;

 15:
 16:     cdef:
+17:         int _num_steps = num_steps
  __pyx_t_5 = __Pyx_PyInt_As_int(__pyx_v_num_steps); if (unlikely((__pyx_t_5 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 17, __pyx_L1_error)
__pyx_v__num_steps = __pyx_t_5;

+18:         double _b = b
  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_v_b); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 18, __pyx_L1_error)
__pyx_v__b = __pyx_t_6;

 19:
 20:     # do the work (field array is updated in-place)
+21:     _evolve2_2_impl(field, _b, _num_steps, scores, _zeros, current)
  __pyx_t_7 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_long(__pyx_v_field);
if (unlikely(!__pyx_t_7.memview)) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_t_8 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(__pyx_v_scores);
if (unlikely(!__pyx_t_8.memview)) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_t_9 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(__pyx_v__zeros);
if (unlikely(!__pyx_t_9.memview)) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_t_10 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_long(__pyx_v_current);
if (unlikely(!__pyx_t_10.memview)) __PYX_ERR(0, 21, __pyx_L1_error)
__pyx_f_46_cython_magic_451bc548913a9407f1c439613b3672aa__evolve2_2_impl(__pyx_t_7, __pyx_v__b, __pyx_v__num_steps, __pyx_t_8, __pyx_t_9, __pyx_t_10);
__PYX_XDEC_MEMVIEW(&__pyx_t_7, 1);
__pyx_t_7.memview = NULL;
__pyx_t_7.data = NULL;
__PYX_XDEC_MEMVIEW(&__pyx_t_8, 1);
__pyx_t_8.memview = NULL;
__pyx_t_8.data = NULL;
__PYX_XDEC_MEMVIEW(&__pyx_t_9, 1);
__pyx_t_9.memview = NULL;
__pyx_t_9.data = NULL;
__PYX_XDEC_MEMVIEW(&__pyx_t_10, 1);
__pyx_t_10.memview = NULL;
__pyx_t_10.data = NULL;

 22:
 23:     # convert memoryviews to numpy arrays
+24:     return np.asarray(field)
  __Pyx_XDECREF(__pyx_r);
__pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2);
__pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_asarray); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_2 = NULL;
if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
__pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
if (likely(__pyx_t_2)) {
PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
__Pyx_INCREF(__pyx_t_2);
__Pyx_INCREF(function);
__Pyx_DECREF_SET(__pyx_t_4, function);
}
}
if (!__pyx_t_2) {
__pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_v_field); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
} else {
#if CYTHON_FAST_PYCALL
if (PyFunction_Check(__pyx_t_4)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_v_field};
__pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else
#endif
#if CYTHON_FAST_PYCCALL
if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_v_field};
__pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_GOTREF(__pyx_t_1);
} else
#endif
{
__pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 24, __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 = NULL;
__Pyx_INCREF(__pyx_v_field);
__Pyx_GIVEREF(__pyx_v_field);
PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_v_field);
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_3, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 24, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
}
}
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_r = __pyx_t_1;
__pyx_t_1 = 0;
goto __pyx_L0;

 25:
 26:
 27: @cython.cdivision(True)
 28: @cython.boundscheck(False)
+29: cdef void _evolve2_2_impl(long[:, ::1] field, double b, int num_steps,
static void __pyx_f_46_cython_magic_451bc548913a9407f1c439613b3672aa__evolve2_2_impl(__Pyx_memviewslice __pyx_v_field, double __pyx_v_b, int __pyx_v_num_steps, __Pyx_memviewslice __pyx_v_scores, __Pyx_memviewslice __pyx_v__zeros, __Pyx_memviewslice __pyx_v_current) {
int __pyx_v_x;
int __pyx_v_y;
int __pyx_v_L;
int __pyx_v_i;
int __pyx_v_j;
int __pyx_v_ix;
int __pyx_v_jy;
CYTHON_UNUSED int __pyx_v_step;
int __pyx_v_bestX;
int __pyx_v_bestY;
/* … */
/* function exit code */
goto __pyx_L0;
__pyx_L1_error:;
__Pyx_WriteUnraisable("_cython_magic_451bc548913a9407f1c439613b3672aa._evolve2_2_impl", __pyx_clineno, __pyx_lineno, __pyx_filename, 0, 1);
__pyx_L0:;
}

 30:                                   double[:, ::1] scores, double[:, ::1] _zeros,
 31:                                   long[:, ::1] current) nogil:
 32:
 33:     cdef int x, y, L, i, j, ix, jy, step
 34:
+35:     L = field.shape[0]
  __pyx_v_L = (__pyx_v_field.shape[0]);

 36:
+37:     for step in range(num_steps):
  __pyx_t_1 = __pyx_v_num_steps;
for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
__pyx_v_step = __pyx_t_2;

+38:         current[...] = field
    if (unlikely(__pyx_memoryview_copy_contents(__pyx_v_field, __pyx_v_current, 2, 2, 0) < 0)) __PYX_ERR(0, 38, __pyx_L1_error)

+39:         scores[...] = _zeros
    if (unlikely(__pyx_memoryview_copy_contents(__pyx_v__zeros, __pyx_v_scores, 2, 2, 0) < 0)) __PYX_ERR(0, 39, __pyx_L1_error)

 40:
+41:         for x in range(L):
    __pyx_t_3 = __pyx_v_L;
for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
__pyx_v_x = __pyx_t_4;

+42:             for y in range(L):
      __pyx_t_5 = __pyx_v_L;
for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) {
__pyx_v_y = __pyx_t_6;

+43:                 for i in range(-1, 2):
        for (__pyx_t_7 = -1L; __pyx_t_7 < 2; __pyx_t_7+=1) {
__pyx_v_i = __pyx_t_7;

+44:                     for j in range(-1, 2):
          for (__pyx_t_8 = -1L; __pyx_t_8 < 2; __pyx_t_8+=1) {
__pyx_v_j = __pyx_t_8;

+45:                         ix = (x + i) % L
            __pyx_v_ix = ((__pyx_v_x + __pyx_v_i) % __pyx_v_L);

+46:                         jy = (y + j) % L
            __pyx_v_jy = ((__pyx_v_y + __pyx_v_j) % __pyx_v_L);

+47:                         scores[x, y] += (1 - field[ix, jy])
            __pyx_t_9 = __pyx_v_ix;
__pyx_t_10 = __pyx_v_jy;
if (__pyx_t_9 < 0) __pyx_t_9 += __pyx_v_field.shape[0];
if (__pyx_t_10 < 0) __pyx_t_10 += __pyx_v_field.shape[1];
__pyx_t_11 = __pyx_v_x;
__pyx_t_12 = __pyx_v_y;
if (__pyx_t_11 < 0) __pyx_t_11 += __pyx_v_scores.shape[0];
if (__pyx_t_12 < 0) __pyx_t_12 += __pyx_v_scores.shape[1];
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_11 * __pyx_v_scores.strides[0]) )) + __pyx_t_12)) )) += (1 - (*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_9 * __pyx_v_field.strides[0]) )) + __pyx_t_10)) ))));
}
}

 48:
+49:                 if field[x, y] == 1:
        __pyx_t_13 = __pyx_v_x;
__pyx_t_14 = __pyx_v_y;
if (__pyx_t_13 < 0) __pyx_t_13 += __pyx_v_field.shape[0];
if (__pyx_t_14 < 0) __pyx_t_14 += __pyx_v_field.shape[1];
__pyx_t_15 = (((*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_13 * __pyx_v_field.strides[0]) )) + __pyx_t_14)) ))) == 1) != 0);
if (__pyx_t_15) {
/* … */
}
}
}

+50:                     scores[x, y] *= b
          __pyx_t_16 = __pyx_v_x;
__pyx_t_17 = __pyx_v_y;
if (__pyx_t_16 < 0) __pyx_t_16 += __pyx_v_scores.shape[0];
if (__pyx_t_17 < 0) __pyx_t_17 += __pyx_v_scores.shape[1];
*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_16 * __pyx_v_scores.strides[0]) )) + __pyx_t_17)) )) *= __pyx_v_b;

 51:
+52:         for x in range(L):
    __pyx_t_3 = __pyx_v_L;
for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
__pyx_v_x = __pyx_t_4;

+53:             for y in range(L):
      __pyx_t_5 = __pyx_v_L;
for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) {
__pyx_v_y = __pyx_t_6;

+54:                 bestX = x
        __pyx_v_bestX = __pyx_v_x;

+55:                 bestY = y
        __pyx_v_bestY = __pyx_v_y;

+56:                 for i in range(-1, 2):
        for (__pyx_t_7 = -1L; __pyx_t_7 < 2; __pyx_t_7+=1) {
__pyx_v_i = __pyx_t_7;

+57:                     for j in range(-1, 2):
          for (__pyx_t_8 = -1L; __pyx_t_8 < 2; __pyx_t_8+=1) {
__pyx_v_j = __pyx_t_8;

+58:                         ix = (x + i) % L
            __pyx_v_ix = ((__pyx_v_x + __pyx_v_i) % __pyx_v_L);

+59:                         jy = (y + j) % L
            __pyx_v_jy = ((__pyx_v_y + __pyx_v_j) % __pyx_v_L);

+60:                         if (scores[bestX, bestY] < scores[ix, jy]):
            __pyx_t_18 = __pyx_v_bestX;
__pyx_t_19 = __pyx_v_bestY;
if (__pyx_t_18 < 0) __pyx_t_18 += __pyx_v_scores.shape[0];
if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_v_scores.shape[1];
__pyx_t_20 = __pyx_v_ix;
__pyx_t_21 = __pyx_v_jy;
if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_v_scores.shape[0];
if (__pyx_t_21 < 0) __pyx_t_21 += __pyx_v_scores.shape[1];
__pyx_t_15 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_18 * __pyx_v_scores.strides[0]) )) + __pyx_t_19)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_scores.data + __pyx_t_20 * __pyx_v_scores.strides[0]) )) + __pyx_t_21)) )))) != 0);
if (__pyx_t_15) {
/* … */
}
}
}

+61:                             bestX = ix
              __pyx_v_bestX = __pyx_v_ix;

+62:                             bestY = jy
              __pyx_v_bestY = __pyx_v_jy;

 63:
+64:                 field[x, y] = current[bestX, bestY]
        __pyx_t_22 = __pyx_v_bestX;
__pyx_t_23 = __pyx_v_bestY;
if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_v_current.shape[0];
if (__pyx_t_23 < 0) __pyx_t_23 += __pyx_v_current.shape[1];
__pyx_t_24 = __pyx_v_x;
__pyx_t_25 = __pyx_v_y;
if (__pyx_t_24 < 0) __pyx_t_24 += __pyx_v_field.shape[0];
if (__pyx_t_25 < 0) __pyx_t_25 += __pyx_v_field.shape[1];
*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_field.data + __pyx_t_24 * __pyx_v_field.strides[0]) )) + __pyx_t_25)) )) = (*((long *) ( /* dim=1 */ ((char *) (((long *) ( /* dim=0 */ (__pyx_v_current.data + __pyx_t_22 * __pyx_v_current.strides[0]) )) + __pyx_t_23)) )));
}
}
}

In [ ]:
L = 44
field = np.zeros((L, L), dtype=int)
field[L//2, L//2] = 1

%timeit evolve2_2(field, 1.81, 10)


Now, for something completely different: numba just-in-time compilation¶

In [1]:
import numba

In [21]:
def evolve3(field, b, num_steps=1):
L = field.shape[0]
_zeros = np.zeros((L, L), dtype=float)
_int_zeros = np.zeros((L, L), dtype=int)
current = _int_zeros.copy()
scores = _zeros.copy()
evolve3_impl(field, b, num_steps, _zeros, _int_zeros)
return field

# Note a single decorator: @jit

@numba.jit(nopython=True)
def evolve3_impl(field, b, num_steps, _zeros, _int_zeros):
L = field.shape[0]
current = _int_zeros.copy()

for step in range(num_steps):
current = field.copy()
scores = _zeros.copy()

for x in range(L):
for y in range(L):
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
scores[x, y] += (1 - field[ix, jy])

if field[x, y] == 1:
scores[x, y] *= b

for x in range(L):
for y in range(L):
bestX = x
bestY = y
for i in range(-1, 2):
for j in range(-1, 2):
ix = (x + i) % L
jy = (y + j) % L
if (scores[bestX, bestY] < scores[ix, jy]):
bestX = ix
bestY = jy

field[x, y] = current[bestX, bestY]
return field

In [23]:
L = 44
field = np.zeros((L, L), dtype=int)
field[L//2, L//2] = 1

%timeit evolve3(field, 1.81, 10)

100 loops, best of 3: 2.12 ms per loop

In [ ]: