Windows 環境にてCython＋ OpenMPでJulia集合をガリガリ計算する．

#本日は

並列処理の典型的な素材Julia集合をCythonで並列計算していきます．
Julia集合は

![julia.png](https://qiita-image-store.s3.amazonaws.com/0/135978/e2ad5546-489b-4c31-fbf3-a8eecbff1a30.png)

といったような絵でございます．今回はループの処理の並列化に言及するので詳細などは[ジュリア集合（Julia set）](http://mt-soft.sakura.ne.jp/kyozai/excel_vba/320_vba_high/11_julia/main.html)などをご覧くださいまし．結構わかりやすいです．
もっと深めの話ですと[Julia集合の描画アルゴリズムとその精度評価](http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1793-02.pdf)
など良いかもしれません．複素関数論の知識がある方はぜひトライ！

# 環境
Windows10Pro Genuine Intel(R)CPU 0000 @2.00GHz 2.00GHz (20C 40T)
Python 3.5.3＋Cython (Miniconda3にてCythonを導入)

Cythonから　精製　もとい　生成したーコードをビルドするために cl.exe を用いているためVisual C++ 2015 Build Toolsなどが必要になるかもしれません．
# 使用コード
オライリーから出版されているCython１２章のコードを拝借します．

シングルな場合とパラレルにする場合二通りのコードを書いてそれらを比較します．
下記のディレクトリ構造に従ってコードを書いていきます．

```shell-session
.
|-- build.bat
|-- compare_res.py
|-- single
|   |-- calc_julia.pyx
|   `-- setup.py
|-- parallel
|   |-- calc_julia.pyx
|   `-- setup.py
```

## single に配置するもの
```py3:single\calc_julia.pyx
import numpy as np 

cdef inline double norm2(double complex z):
    return z.real*z.real+z.imag*z.imag

cdef int escape(double complex z,
                double complex c,
                double zmax,
                int nmax):
    cdef:
        int i = 0
        double square_zmax = zmax*zmax
    while norm2(z) < zmax and i < nmax:
        z = z*z+c
        i += 1
    return i

def calc_julia(int resolution,
                double complex c,
                double bound=1.5,
                double zmax=4.0,
                int nmax=1000):
    cdef:
        double step = 2.0*bound/resolution
        int i, j
        double complex z
        double real, imag
        int[:, ::1] counts

    counts = np.zeros((resolution+1, resolution+1), dtype=np.int32)

    for i in range(resolution+1):
        real = -bound+i*step
        for j in range(resolution+1):
            imag = -bound+j*step
            z = real+imag*1j
            counts[i, j] = escape(z, c, zmax, nmax)

    return np.asarray(counts)
```
```py3:single\setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_module = Extension(
    name="julia_single",
    sources=["calc_julia.pyx"],
)

setup(
    name = 'julia single',
    cmdclass = {'build_ext': build_ext},
    ext_modules = [ext_module],
)
```

## parallel に配置するもの

```py3:parallel\calc_julia.pyx
import numpy as np 
from cython.parallel import prange
from cython cimport boundscheck, wraparound

cdef inline double norm2(double complex z)nogil:
    return z.real*z.real+z.imag*z.imag

cdef int escape(double complex z,
                double complex c,
                double zmax,
                int nmax)nogil:
    cdef:
        int i = 0
        double square_zmax = zmax*zmax
    while norm2(z) < zmax and i < nmax:
        z = z*z+c
        i += 1
    return i


@boundscheck(False)
@wraparound(False)
def calc_julia(int resolution,
                double complex c,
                double bound=1.5,
                double zmax=4.0,
                int nmax=1000):
    cdef:
        double step = 2.0*bound/resolution
        int i, j
        double complex z
        double real, imag
        int[:, ::1] counts

    counts = np.zeros((resolution+1, resolution+1), dtype=np.int32)

    for i in prange(resolution+1,nogil=True,schedule='static',chunksize=1):
        real = -bound+i*step
        for j in range(resolution+1):
            imag = -bound+j*step
            z = real+imag*1j
            counts[i, j] = escape(z, c, zmax, nmax)

    return np.asarray(counts)
```

ポイントなのは `prange` を導入し内部で用いる関数に対してGILを解除するように明示することです．

```py3:parallel\setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_module = Extension(
    name="julia_parallel",
    sources=["calc_julia.pyx"],
    extra_compile_args=['/openmp'],
)

setup(
    name = 'julia parallel',
    cmdclass = {'build_ext': build_ext},
    ext_modules = [ext_module],
)
```

`setup.py`　はOpenMPを使ってねとコンパイラにお願いするオプションをつけます．

## その他ツール

```shell-session
> python setup.py build_ext --inplace
```

を各ディレクトリで実行するだけのバッチを作ります．

```bat:built.bat
cd single 
python setup.py build_ext --inplace
cd ..

cd parallel
python setup.py build_ext --inplace
cd ..
```

## 検証用コード

シングルとパラレルの場合で各々処理時間を計測しどれくらい効率よくなってるか確かめます．原理的には結果が同じになってしかるべきですので　`assert` で確認しています．

```py3:compare_res.py
from single import julia_single
from parallel import julia_parallel
import time

resolution = 5000
c = (0.322+0.05j)
# single
start = time.time()
jl_single = julia_single.calc_julia(resolution, c)
end = time.time()
single_time=end-start
print("single elapsed =", single_time)
# parallel
start = time.time()
jl_parallel = julia_parallel.calc_julia(resolution, c)
end = time.time()
parallel_time=end-start
print("parallel elapsed =", parallel_time)

print("single/parallel =",single_time/parallel_time)

assert (jl_single == jl_parallel).all()

from matplotlib import pyplot as plt 
import numpy as np
plt.imshow(np.log(jl_single))
plt.show()
```

## 実行例

```shell-session
> build.bat
> python compare_res.py
single elapsed = 52.052311420440674
parallel elapsed = 1.4532854557037354
single/parallel = 35.81699054108746
```

効率が35倍ほど改善されており原理的に得られる効率の倍数であるスレッド数40に近いですね．ちなみにresolution のパラメーターを5000から10000にしておおきくすると

```shell-session
single elapsed = 209.3094298839569
parallel elapsed = 5.563415050506592
single/parallel = 37.62247252519793
```
という結果で効率が＋２増えました．パラレルで動かしているときにCPUが100％回ってるので効率はこれ以上増えないと思います．

# 最後に
- Windows環境での文献がないのでとりあえず試せてよかった．
- 40TのCPUをフルに回せることができて満足．
![40T.png](https://qiita-image-store.s3.amazonaws.com/0/135978/5cf5ba29-5163-fa42-5ad4-e2257f5c12b4.png)

写真は resolution を 50000 で回してます．