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_AddTraceback("_cython_magic_4a51e9818ff378ac421ff64c016d5995.evolve2_1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __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_AddTraceback("_cython_magic_4a51e9818ff378ac421ff64c016d5995.evolve2_1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __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_AddTraceback("_cython_magic_451bc548913a9407f1c439613b3672aa.evolve2_2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __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_AddTraceback("_cython_magic_451bc548913a9407f1c439613b3672aa.evolve2_2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __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