[SciPy-user] threading and scipy

[Xavier Gnata]

Damian Eads wrote:
> Hi there,
>
> I am running some of my code through a very large data set on quad-core 
> cluster nodes. A simple grep confirms that most parts of Scipy (e.g. 
> linalg) do not use the Py_BEGIN_ALLOW_THREADS and Py_END_ALLOW_THREADS 
> macros (or the numpy equivalents).
>
> [eads at pumpkin scipy]$ grep ALLOW_THREADS `find ~/work/repo/scipy -name 
> "*.[ch]*"` | grep -v "build" | sed 's/:.*//g' | sort | uniq
> /home/eads/work/repo/scipy/scipy/sandbox/netcdf/_netcdf.c
> /home/eads/work/repo/scipy/scipy/sandbox/netcdf/.svn/text-base/_netcdf.c.svn-base
> /home/eads/work/repo/scipy/scipy/stsci/convolve/src/_lineshapemodule.c
> /home/eads/work/repo/scipy/scipy/stsci/convolve/src/.svn/text-base/_lineshapemodule.c.svn-base
>
> Numpy seems to have a lot more coverage, though.
>
> [eads at pumpkin scipy]$ grep ALLOW_THREADS `find ~/work/numpy-1.0.4/numpy 
> -name "*.c"` | sed 's/:.*//g' | sort | uniq
> /home/eads/work/numpy-1.0.4/numpy/core/blasdot/_dotblas.c
> /home/eads/work/numpy-1.0.4/numpy/core/src/arrayobject.c
> /home/eads/work/numpy-1.0.4/numpy/core/src/multiarraymodule.c
>
> [eads at pumpkin scipy]$
>
> Is it true if my code is heavily dependent on Scipy (I do image 
> processing on large images with ndimage) and I use the %bg command in 
> IPython, most of the time there will be only one thread running with the 
> others blocked?
>
> I anticipate others will insist that I read up on the caveats of 
> multi-threaded programming (mutual exclusion, locking, critical regions, 
> etc.) so I should mention that I am a pretty seasoned with it, having 
> done quite a bit of work with pthreads. However, I am new to threading 
> in python and I heard there are issues, specifically only one thread is 
> allowed access to the global interpreter lock at a time.
>
> I would like to run some filters on 300 images. These filters change 
> from one iteration to the next of the program. When all the filtering is 
> finished, a single thread needs to see the result of all the computation 
> (all the result images) to compute so inter-image statistics. Then, I 
> start the process over. I'd really like to spawn 4+ threads, one each 
> working on a different image.
>
> Being that I don't see any code in ndimage that releases the global 
> interpreter lock, is it true that if I wrote code to spawn separate 
> filter threads, only one would execute at a time?
>
> Please advise.
>
> Thank you!
>
> Damian
>   
Well if you have 300 images you should run let say 6 times the same 
python *script* (not thread) on 6 sets of 50 images. It is what we call 
" embarrassingly parallel". 
http://en.wikipedia.org/wiki/Embarrassingly_parallel
No threads, no locks, no mutex :)

Now, if now want to perform some explicit multi-threading computations 
in the numpy/scipy framework,
you could try to insert C code + openmp when the performances matter.

One very stupid example is :
import scipy
from scipy.weave import converters
from scipy import zeros
import pylab

max = 100000

code = """
#include <math.h>
long i;
double s=1;
long k=10000;
#pragma omp parallel for
for (i=1; i<max; i++) {
for (k=1; k<max; k++) {
s +=i;
}
}
return_val = s;
        """
val = scipy.weave.inline(code,
             ['max'],
             type_converters=converters.blitz,
             compiler = 'gcc',extra_compile_args 
=['-fopenmp'],libraries=['gomp'])


It works at least with gcc-4.3 bug not with gcc-4.2 on my ubuntu.
Using scipy.weave and openmp, you can write threaded 
performance-critical sections of code.
You do not have to learn the details of the pthreads API and that is why 
openmp is cool.

The next step is either MPI or explicit threading using a thread lib 
(but I do hope you don(t need this because it is another story).

"I'd really like to spawn 4+ threads, one each working on a different image.": No. You should really run 4 times the same script at the same time. There are several solutions to know when these 4 exe are done but if someone can provide us with a *clean and easy* way...

Xavier