pythonparallel-processingcythonphysicsmontecarlo# How can I use 'prange' in Cython?

I'm trying to solve a two-dimensional Ising model with a Monte Carlo approach.

As it is slow, I used Cython to accelerate the code execution. I would like to push it even further and parallelize the Cython code. My idea is to split the two-dimensional lattice in two, so for any point on a lattice has its nearest neighbours on the other lattice. This way, I can randomly choose one lattice, and I can flip all the spins and this could be done in parallel since all those spins are independent.

So far this is my code:_{(inspired from http://jakevdp.github.io/blog/2017/12/11/live-coding-cython-ising-model/)}

```
%load_ext Cython
%%cython
cimport cython
cimport numpy as np
import numpy as np
from cython.parallel cimport prange
@cython.boundscheck(False)
@cython.wraparound(False)
def cy_ising_step(np.int64_t[:, :] field,float beta):
cdef int N = field.shape[0]
cdef int M = field.shape[1]
cdef int offset = np.random.randint(0,2)
cdef np.int64_t[:,] n_update = np.arange(offset,N,2,dtype=np.int64)
cdef int m,n,i,j
for m in prange(M,nogil=True):
i = m % 2
for j in range(n_update.shape[0]) :
n = n_update[j]
cy_spin_flip(field,(n+i) %N,m%M,beta)
return np.array(field,dtype=np.int64)
cdef cy_spin_flip(np.int64_t[:, :] field,int n,int m, float beta=0.4,float J=1.0):
cdef int N = field.shape[0]
cdef int M = field.shape[1]
cdef float dE = 2*J*field[n,m]*(field[(n-1)%N,m]+field[(n+1)%N,m]+field[n,(m-1)%M]+field[n,(m+1)%M])
if dE <= 0 :
field[n,m] *= -1
elif np.exp(-dE * beta) > np.random.rand():
field[n,m] *= -1
```

I tried using a `prange`

-constructor, but I'm having a lots of troubles with GIL-lock. I’m new to Cython and parallel computing so I could easily have missed something.

**The error:**

```
Discarding owned Python object not allowed without gil
Calling gil-requiring function not allowed without gil
```

Solution

From a Cython point-of-view the main problem is that `cy_spin_flip`

requires the GIL. You need to add `nogil`

to the end of its signature, and set the return type to `void`

(since by default it returns a Python object, which requires the GIL).

However, `np.exp`

and `np.random.rand`

also require the GIL, because they're Python function calls. `np.exp`

is probably easily replaced with `libc.math.exp`

. `np.random`

is a bit harder, but there's plenty of suggestions for C- and C++-based approaches: 1 2 3 4 (+ others).

A more fundamental problem is the line:

```
cdef float dE = 2*J*field[n,m]*(field[(n-1)%N,m]+field[(n+1)%N,m]+field[n,(m-1)%M]+field[n,(m+1)%M])
```

You've parallelized this with respect to `m`

(i.e. different values of `m`

are run in different threads), and each iteration changes `field`

. However in this line you are looking up several different values of `m`

. This means the whole thing is a race-condition (the result depends on which order the different threads finish) and suggests **your algorithm may be fundamentally unsuitable for parallelization**. Or that you should copy `field`

and have `field_in`

and `field_out`

. It isn't obvious to me, but this is something that you should be able to work out.

*Edit:* it does look like you've given the race condition some thought with using `i%2`

. It isn't obvious to me that this is right though. I think a working implementation of your "alternate cells" scheme would look something like:

```
for oddeven in range(2):
for m in prange(M):
for n in range(N):
# some mechanism to pick the alternate cells here.
```

i.e. you need a regular loop to pick the alternate cells outside your parallel loop.

- Python Jinja2 LaTeX Table
- Getting attributes of a class
- How can I print many significant figures in Python?
- How to allow list append() method to return the new list
- Calculate Last Friday of Month in Pandas
- Python type hint for Iterable[str] that isn't str
- How to iterate over a list in chunks
- How to exit the entire application from a Python thread?
- Running shell command and capturing the output
- How do I pass a variable by reference?
- Convert range(r) to list of strings of length 2 in python
- How can I get the start and end dates for each week?
- how to use send_message() in python-telegram-bot
- Python conditional replacement based on element type
- How can I count the number of items in an arbitrary iterable (such as a generator)?
- Find longest consecutive range of numbers in list
- Insert text in braces with asyncpg
- How does one put a link / url to the web-site's home page in Django?
- How to determine if a path is a subdirectory of another?
- Custom Keybindings for Ipython terminal
- FastAPI asynchronous background tasks blocks other requests?
- How to make sure that information from one file is duplicated into several text documents, without specific lines
- Installing a Python environment with Anaconda
- sklearn pipeline model predicting same results for all input
- Brew command not found after installing Anaconda Python
- How to get an XPath from selenium webelement or from lxml?
- Pipe PuTTY console to Python script
- How to align the axes of a figure in matplotlib?
- Persist ParentDocumentRetriever of langchain
- How to reset index in a pandas dataframe?