Pythonでマンデルブロ集合の計算を高速化する

Last updated at 2024-12-26Posted at 2024-12-25

マンデルブロ集合とは

マンデルブロ集合は、複素数平面上の特定の点の集合で、数学的にも視覚的にも非常に興味深い性質を持つフラクタル図形です。これはフランスの数学者ブノワ・マンデルブロ（Benoît B. Mandelbrot）の名前にちなんで名付けられました。

定義

マンデルブロ集合は、次の反復関数に基づいて定義されます：

z_{n+1} = z_n^2 + c

$z$ : 複素数（初期値は $z_0 = 0$)
$c$ : 複素平面上の任意の点

マンデルブロ集合は、以下の条件を満たす複素数 $ c $ の集合です：

反復を無限に続けても $ |z_n| \leq 2 $ を満たす場合、$ c $ はマンデルブロ集合に含まれます。
一方、反復によって $ |z_n| $ が $ 2 $ を超えると、$ c $ はマンデルブロ集合に含まれません（発散するとみなします）。

直感的な説明

初期条件：
- 各点 $ c $ について、$ z_0 = 0 $ を初期値として計算を始めます。
反復の結果：
- $ z $ が無限に反復しても発散しない（つまり、ある範囲内に留まる）場合、その点 $ c $ はマンデルブロ集合に属します。
視覚化：
- $ c $ を複素数平面上の点と考え、それぞれの点が集合に含まれるかどうかを計算します。
- 含まれる点を黒で描き、発散する点を発散の速さに応じて異なる色で塗ると、独特な図形が現れます。

マンデルブロ集合の特徴

フラクタル性：
- 拡大しても似た構造が繰り返し現れる自己相似性を持っています。
複雑性と単純性の融合：
- 簡単な数学式 $ z_{n+1} = z_n^2 + c $ に基づいていますが、生成される図形は非常に複雑です。
境界の無限の複雑さ：
- 境界部分を拡大すると、無限の細部が現れます。これがフラクタル図形の典型的な特徴です。

応用と意義

数学と芸術：
- マンデルブロ集合の形状は、数学的興味を超え、芸術的な美しさを持つことで知られています。
物理学と自然現象：
- フラクタルは自然界のパターン（例えば、海岸線の形状や雪の結晶）をモデル化する際にも用いられます。
カオス理論：
- マンデルブロ集合は、カオスと秩序の境界を示す例としても研究されています。

視覚的な例

マンデルブロ集合を描画すると、中心に独特な形状が現れ、周囲には複雑で美しい模様が広がります。拡大するほど新しいパターンが現れ、数学と芸術の交差点としても広く知られています。

マンデルブロ集合を描くコード

基本的な実装を以下に示します。計算に17.12秒かかりました。

import matplotlib.pyplot as plt
import numpy as np


def mandelbrot(c, max_iter):
    """
    マンデルブロ集合の点が発散するかどうかを判定します。
    発散までにかかる反復回数を返します。
    """
    z = 0
    for n in range(max_iter):
        if abs(z) > 2:
            return n
        z = z * z + c
    return max_iter


def mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter):
    """
    指定した範囲のマンデルブロ集合を計算します。
    """
    # グリッドを作成
    x = np.linspace(xmin, xmax, width)
    y = np.linspace(ymin, ymax, height)
    mandelbrot_image = np.zeros((height, width))

    for i in range(height):
        for j in range(width):
            c = x[j] + 1j * y[i]
            mandelbrot_image[i, j] = mandelbrot(c, max_iter)
    return mandelbrot_image


# パラメータ設定
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width, height = 800, 800
max_iter = 100

# マンデルブロ集合を計算
image = mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter)

# 描画
plt.figure(figsize=(10, 10))
plt.imshow(image, extent=(xmin, xmax, ymin, ymax), cmap="hot", interpolation="nearest")
plt.colorbar(label="Iterations")
plt.title("Mandelbrot Set")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

並列化による高速化

joblibを使用してマルチプロセスで高速化します。CPU coreは16です。計算に3.14秒かかりました。

import time

import matplotlib.pyplot as plt
import numpy as np
from joblib import Parallel, cpu_count, delayed


def mandelbrot(c, max_iter):
    """
    マンデルブロ集合の点が発散するかどうかを判定します。
    発散までにかかる反復回数を返します。
    """
    z = 0
    for n in range(max_iter):
        if abs(z) > 2:
            return n
        z = z * z + c
    return max_iter


def mandelbrot_row(y, x_vals, max_iter):
    """
    1行分のマンデルブロ集合を計算します。
    """
    return [mandelbrot(x + 1j * y, max_iter) for x in x_vals]


def mandelbrot_set_parallel(xmin, xmax, ymin, ymax, width, height, max_iter, n_jobs=-1):
    """
    指定した範囲のマンデルブロ集合を並列計算します。
    """
    # グリッドを作成
    x_vals = np.linspace(xmin, xmax, width)
    y_vals = np.linspace(ymin, ymax, height)

    # 並列計算
    mandelbrot_image = Parallel(n_jobs=n_jobs)(
        delayed(mandelbrot_row)(y, x_vals, max_iter) for y in y_vals
    )
    return np.array(mandelbrot_image)


# パラメータ設定
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width, height = 2000, 2000  # 高解像度
max_iter = 100
n_jobs = -1  # 全ての利用可能なコアを使用

# 使用するコア数を取得
num_cores = cpu_count() if n_jobs == -1 else n_jobs
print(f"Number of CPU cores available: {num_cores}")

# 計算時間を計測
start_time = time.time()
image = mandelbrot_set_parallel(
    xmin, xmax, ymin, ymax, width, height, max_iter, n_jobs=n_jobs
)
end_time = time.time()

# 計算時間を表示
print(f"Calculation time with parallelization: {end_time - start_time:.2f} seconds")
print(f"Number of CPU cores used: {num_cores}")

# 描画
plt.figure(figsize=(12, 12))
plt.imshow(image, extent=(xmin, xmax, ymin, ymax), cmap="hot", interpolation="nearest")
plt.colorbar(label="Iterations")
plt.title("Mandelbrot Set (High Resolution, Parallelized)")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

Numbaによる高速化

Numbaをインストールしてください。

pip install numba

一回目の実行時間は2.01秒、二回目以降は1.74秒程度でした。一回目はコンパイルが行われるため、遅いです。

import time

import matplotlib.pyplot as plt
import numpy as np

from numba import jit


@jit(nopython=True)
def mandelbrot(c, max_iter):
    """
    マンデルブロ集合の点が発散するかどうかを判定します。
    発散までにかかる反復回数を返します。
    """
    z = 0
    for n in range(max_iter):
        if abs(z) > 2:
            return n
        z = z * z + c
    return max_iter


@jit(nopython=True, parallel=True)
def mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter):
    """
    指定した範囲のマンデルブロ集合を計算します。
    """
    # グリッドを作成
    x = np.linspace(xmin, xmax, width)
    y = np.linspace(ymin, ymax, height)
    mandelbrot_image = np.zeros((height, width))

    for i in range(height):
        for j in range(width):
            c = x[j] + 1j * y[i]
            mandelbrot_image[i, j] = mandelbrot(c, max_iter)
    return mandelbrot_image


# パラメータ設定
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width, height = 2000, 2000  # 高解像度
max_iter = 100

# 計算時間を計測
start_time = time.time()
image = mandelbrot_set(xmin, xmax, ymin, ymax, width, height, max_iter)
end_time = time.time()

# 計算時間を表示
print(f"Calculation time (with numba): {end_time - start_time:.2f} seconds")

# 描画
plt.figure(figsize=(12, 12))
plt.imshow(image, extent=(xmin, xmax, ymin, ymax), cmap="hot", interpolation="nearest")
plt.colorbar(label="Iterations")
plt.title("Mandelbrot Set (High Resolution, Numba Optimized)")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

Cythonによる高速化

Cythonをインストールしてください。

pip install cython

マンデルブロ集合を計算する関数をCythonコードとして書きます。

mandelbrot.pyx

cimport cython
import numpy as np
cimport numpy as np

@cython.boundscheck(False)  # 境界チェックを無効化
@cython.wraparound(False)   # 配列の負のインデックスを無効化
def mandelbrot_set_cython(double xmin, double xmax, double ymin, double ymax,
                          int width, int height, int max_iter):
    cdef:
        double[:] x = np.linspace(xmin, xmax, width)
        double[:] y = np.linspace(ymin, ymax, height)
        int[:, :] mandelbrot_image = np.zeros((height, width), dtype=np.int32)
        int i, j, n
        double zx, zy, zx2, zy2, cx, cy

    for i in range(height):
        for j in range(width):
            cx = x[j]
            cy = y[i]
            zx = zy = 0.0
            zx2 = zy2 = 0.0
            for n in range(max_iter):
                if zx2 + zy2 > 4.0:
                    break
                zy = 2.0 * zx * zy + cy
                zx = zx2 - zy2 + cx
                zx2 = zx * zx
                zy2 = zy * zy
            mandelbrot_image[i, j] = n

    return np.array(mandelbrot_image)

Cythonコードをコンパイルするためのコードを書きます。

setup.py

import numpy as np
from setuptools import setup

from Cython.Build import cythonize

setup(ext_modules=cythonize("mandelbrot.pyx"), include_dirs=[np.get_include()])

Cythonコードをコンパイルします。

python setup.py build_ext --inplace

メインコードを書きます。

main.py

import time

import matplotlib.pyplot as plt
from mandelbrot import mandelbrot_set_cython

# パラメータ設定
xmin, xmax, ymin, ymax = -2.0, 1.0, -1.5, 1.5
width, height = 2000, 2000  # 高解像度
max_iter = 100

# 計算時間を計測
start_time = time.time()
image = mandelbrot_set_cython(xmin, xmax, ymin, ymax, width, height, max_iter)
end_time = time.time()

# 計算時間を表示
print(f"Calculation time with Cython: {end_time - start_time:.2f} seconds")

# 描画
plt.figure(figsize=(12, 12))
plt.imshow(image, extent=(xmin, xmax, ymin, ymax), cmap="hot", interpolation="nearest")
plt.colorbar(label="Iterations")
plt.title("Mandelbrot Set (High Resolution, Cython Optimized)")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

メインコードを実行するとマンデルブロ集合を計算できます。計算に0.17秒かかりました。

実行環境に依存しますが、ファイル構成は以下のような構造になると思います。

>TREE /F
C:.
│  main.py
│  mandelbrot.c
│  mandelbrot.cp310-win_amd64.pyd
│  mandelbrot.pyx
│  setup.py
│
└─build
    ├─lib.win-amd64-cpython-310
    │      mandelbrot.cp310-win_amd64.pyd
    │
    └─temp.win-amd64-cpython-310
        └─Release
                mandelbrot.cp310-win_amd64.exp
                mandelbrot.cp310-win_amd64.lib
                mandelbrot.obj

ファイル構成の解説

以下は、Cythonを使用してマンデルブロ集合を高速化したプロジェクトのファイル構成と、それぞれの役割についての説明です。

ルートディレクトリ

1. `main.py`

役割: メインのPythonスクリプト。
内容:
- Cythonで高速化されたマンデルブロ集合関数を呼び出し、結果を描画します。
- 高解像度画像を生成し、計算時間を計測します。
実行方法: python main.py でスクリプトを実行します。

2. `mandelbrot.c`

役割: Cythonによって生成されたC言語のソースコード。
内容:
- mandelbrot.pyx を基に、自動生成されたCコード。
- PythonとCの間でやり取りを行うためのコードや、高速化されたマンデルブロ集合の計算ロジックが含まれています。
生成方法: setup.py を使って Cython 拡張をビルドする際に生成されました。
注意: 通常、このファイルを直接編集する必要はありません。

3. `mandelbrot.cp310-win_amd64.pyd`

役割: Windows用のPython拡張モジュール（Cライブラリ）。
内容:
- コンパイルされたバイナリコード。
- .pyd ファイルは、Windows環境で .so（共有ライブラリ）に相当する形式。
- このモジュールを main.py からインポートすることで、Cythonで高速化されたマンデルブロ関数を利用できます。

4. `mandelbrot.pyx`

役割: Cythonソースコード。
内容:
- PythonとC言語の中間的な形式で記述されたコード。
- マンデルブロ集合の計算ロジックをC言語に近い形で記述し、速度を向上させています。
使用方法:
- setup.py を使ってコンパイルします。
- コンパイル結果として .c ファイルと .pyd ファイルが生成されます。

5. `setup.py`

役割: Cythonソースコードをコンパイルするためのスクリプト。
内容:
- setuptools と Cython を利用して、.pyx ファイルをCコードに変換し、さらにそれをコンパイルしてPython拡張モジュール（.pyd）を生成します。
使用方法:
```
python setup.py build_ext --inplace
```
を実行することで、コンパイルが行われます。

`build` ディレクトリ

1. `lib.win-amd64-cpython-310/mandelbrot.cp310-win_amd64.pyd`

役割: 同じディレクトリ構成でインストール可能なPythonモジュール。
内容:
- setup.py の build_ext コマンドによって生成されたバイナリ拡張。
- 同じ名前の .pyd ファイルですが、モジュールがビルドディレクトリに配置されています。

2. `temp.win-amd64-cpython-310/Release/`

役割: ビルド時に生成される一時ファイルを格納するディレクトリ。
内容:
- mandelbrot.cp310-win_amd64.exp: .pyd ファイルを生成する際にリンク情報を含むエクスポートファイル。
- mandelbrot.cp310-win_amd64.lib: .pyd ファイルに対応するライブラリファイル。
- mandelbrot.obj: .pyx ファイルから生成されたオブジェクトファイル。

ファイル構成の流れ

mandelbrot.pyx:
- ソースコードをCython形式で記述。
setup.py:
- mandelbrot.pyx を C コード（mandelbrot.c）に変換し、さらにコンパイルして .pyd モジュールを生成。
mandelbrot.cp310-win_amd64.pyd:
- Pythonからインポート可能な拡張モジュール。
main.py:
- mandelbrot.cp310-win_amd64.pyd をインポートし、高速なマンデルブロ集合の計算を実行。

注意事項

.pyd ファイルは、Pythonのバージョン（例: cp310）やOS環境（例: win_amd64）に依存します。他の環境では再コンパイルが必要です。
build ディレクトリ内のファイルは、中間生成物であり、通常は直接操作する必要はありません。

長々と書きましたが、コンパイルした後に使うものは以下のファイルだけです。

>TREE /F
C:.
    main.py
    mandelbrot.cp310-win_amd64.pyd

これでmain.pyを実行できます。

計算環境

Microsoft Windows [Version 10.0.26100.2605]
13th Gen Intel(R) Core(TM) i5-13400F
TotalPhysicalMemory 34055766016
Python 3.10.15
numpy 1.26.3
joblib 1.4.2
Cython 3.0.11
numba 0.60.0

補足：Cythonコードの解説

以下に、Cythonコードの各行について詳細な解説を行います。

ヘッダー部

1. `cimport cython`

役割: Cythonの特別な機能（例: 型指定や最適化オプション）をインポートします。
詳細: 通常のPython import ではなく、Cython専用のキーワード cimport を使用。

2. `import numpy as np`

役割: NumPyライブラリをインポートします。
詳細: 通常のPythonでのNumPyのインポート。

3. `cimport numpy as np`

役割: NumPyのCython専用インターフェイスをインポートします。
詳細: 配列の型指定やCレベルの効率的な操作を可能にします。

デコレーター

4. `@cython.boundscheck(False)`

役割: 配列アクセス時の境界チェックを無効化します。
詳細: パフォーマンス向上のため、範囲外アクセスがないことを前提にチェックを省略。

5. `@cython.wraparound(False)`

役割: 配列インデックスでの負の値（例: -1）の対応を無効化します。
詳細: Cythonで高速化する際、負のインデックスが不要な場合に使います。

関数ヘッダー

6. `def mandelbrot_set_cython(double xmin, double xmax, double ymin, double ymax, int width, int height, int max_iter):`

役割: マンデルブロ集合を計算する関数の定義。
引数:
- xmin, xmax, ymin, ymax（double）: 計算領域の範囲（実部と虚部の範囲）。
- width, height（int）: 計算領域の解像度。
- max_iter（int）: 各点の最大反復回数。

型指定と変数定義

7. `cdef:`

役割: Cythonで使用するC言語型の変数を宣言。

8. `double[:] x = np.linspace(xmin, xmax, width)`

役割: 実数軸の範囲を均等に分割してNumPy配列として生成。
型指定: double[:] は1次元配列（C配列のような形式）。

9. `double[:] y = np.linspace(ymin, ymax, height)`

役割: 虚数軸の範囲を均等に分割してNumPy配列として生成。

10. `int[:, :] mandelbrot_image = np.zeros((height, width), dtype=np.int32)`

役割: 計算結果を格納する2次元配列を初期化。
型指定: int[:, :] は2次元配列。

11. `int i, j, n`

役割: ループ用の整数型変数を宣言。

12. `double zx, zy, zx2, zy2, cx, cy`

役割: マンデルブロ集合計算用の実数型変数を宣言。

ループと計算

13. `for i in range(height):`

役割: 虚数軸方向（高さ）を反復。

14. `for j in range(width):`

役割: 実数軸方向（幅）を反復。

15. `cx = x[j]`

役割: 現在の点の実部を取得。

16. `cy = y[i]`

役割: 現在の点の虚部を取得。

17. `zx = zy = 0.0`

役割: z（現在の値）の初期化（実部と虚部を0に設定）。

18. `zx2 = zy2 = 0.0`

役割: zx と zy の2乗値を保持する変数を初期化。

19. `for n in range(max_iter):`

役割: マンデルブロ集合の反復計算を実行。

20. `if zx2 + zy2 > 4.0:`

役割: 現在の点が発散するかどうかを判定。
詳細: 発散条件は (|z|^2 = zx^2 + zy^2 > 4)。

21. `break`

役割: 発散を確認した場合、ループを終了。

22. `zy = 2.0 * zx * zy + cy`

役割: 次の反復での虚部を計算。

23. `zx = zx2 - zy2 + cx`

役割: 次の反復での実部を計算。

24. `zx2 = zx * zx`

役割: 実部の2乗を計算。

25. `zy2 = zy * zy`

役割: 虚部の2乗を計算。

26. `mandelbrot_image[i, j] = n`

役割: 現在の点の反復回数を結果配列に格納。

結果の返却

27. `return np.array(mandelbrot_image)`

役割: 計算結果の2次元配列をNumPy配列として返却。

補足まとめ

このコードは、マンデルブロ集合の各点を効率的に計算するために、Cythonの型指定や最適化機能を活用しています。特にboundscheckやwraparoundを無効化することで、計算速度が向上しています。また、np.linspaceとNumPy配列の型指定を組み合わせることで、PythonとC言語のパフォーマンスのギャップを埋めています。

おまけ：OpenMPによるマルチスレッド処理

CythonコードをOpenMPを使用して並列化したバージョンです。並列化することで、複数のコアを利用してマンデルブロ集合の計算を高速化します。

cimport cython
from cython.parallel import prange
import numpy as np
cimport numpy as np
import os

@cython.boundscheck(False)  # 境界チェックを無効化
@cython.wraparound(False)   # 配列の負のインデックスを無効化
@cython.cdivision(True)     # 除算のエラーを無効化
def mandelbrot_set_cython(double xmin, double xmax, double ymin, double ymax,
                          int width, int height, int max_iter):
    """
    OpenMPを使用してマンデルブロ集合の計算を並列化
    """
    cdef:
        double[:] x = np.linspace(xmin, xmax, width)
        double[:] y = np.linspace(ymin, ymax, height)
        int[:, :] mandelbrot_image = np.zeros((height, width), dtype=np.int32)
        int i, j, n
        double zx, zy, zx2, zy2, cx, cy
        int num_threads = os.cpu_count()  # 使用するスレッド数を動的に取得

    # 並列化されたループ
    for i in prange(height, nogil=True, num_threads=num_threads):  # OpenMPを使用して並列化
        for j in range(width):
            cx = x[j]
            cy = y[i]
            zx = zy = 0.0
            zx2 = zy2 = 0.0
            for n in range(max_iter):
                if zx2 + zy2 > 4.0:
                    break
                zy = 2.0 * zx * zy + cy
                zx = zx2 - zy2 + cx
                zx2 = zx * zx
                zy2 = zy * zy
            mandelbrot_image[i, j] = n

    return np.array(mandelbrot_image)

修正点の説明

from cython.parallel import prange:
- prange を使用して外側のループ（for i in range(height)）を並列化します。
- prange はOpenMPの並列化を可能にするCythonの特殊な構文です。
nogil=True:
- Pythonのグローバルインタプリタロック（GIL）を解放して並列処理を許可します。
int num_threads = os.cpu_count():
- 使用するスレッド数を指定します。スレッド数は環境に合わせて変更可能です。
内側ループの非並列化:
- OpenMPのルール上、内側のループをそのまま保持することで競合を避けます。
ゼロ除算のエラー無効化
- @cython.cdivision(True)（除算のエラー無効化）は、Cythonで速度最適化のために使用されるオプションです。ゼロ除算が発生しないことが分かっているので、速度向上のために使用します。

OpenMPを使用するため、setup.pyを修正してOpenMPサポートを有効にします。

from setuptools import setup, Extension
from Cython.Build import cythonize
import numpy as np

extensions = [
    Extension(
        "mandelbrot",
        ["mandelbrot.pyx"],
        extra_compile_args=["-openmp"],  # OpenMPを有効化
        include_dirs=[np.get_include()],
    )
]

setup(
    ext_modules=cythonize(extensions),
)

MSVC専用の OpenMP オプションである -openmp を使用しています。GCCやClangでは -fopenmp を使用します。

LLVMとBuild Tools for Visual Studioをインストールします。

LLVM

githubのサイトからWindows用バイナリを入手。今回はLLVM-19.1.6-win64.exeを選択。exeを実行してそのままインストール。途中のオプションでパスを通しておくと便利。

openmpが入っているかは"C:\Program Files\LLVM\bin\libomp.dll"があるかを確認する。

Build Tools for Visual Studio

MicrosoftのサイトからBuild Tools for Visual Studio 2022を選択。exeを実行して、C++によるデスクトップ開発をインストールする。

コンパイルします。

python setup.py build_ext --inplace

メインコードは同じです。実行するとマンデルブロ集合を計算できます。計算に0.05秒かかりました。

PythonAdvent Calendar 2024

Day 14

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Pythonでマンデルブロ集合の計算を高速化する

マンデルブロ集合とは

定義

直感的な説明

マンデルブロ集合の特徴

応用と意義

視覚的な例

マンデルブロ集合を描くコード

並列化による高速化

Numbaによる高速化

Cythonによる高速化

ファイル構成の解説

ルートディレクトリ

1. main.py

2. mandelbrot.c

3. mandelbrot.cp310-win_amd64.pyd

4. mandelbrot.pyx

5. setup.py

build ディレクトリ

1. lib.win-amd64-cpython-310/mandelbrot.cp310-win_amd64.pyd

2. temp.win-amd64-cpython-310/Release/

ファイル構成の流れ

注意事項

計算環境

補足：Cythonコードの解説

ヘッダー部

1. cimport cython

2. import numpy as np

3. cimport numpy as np

デコレーター

4. @cython.boundscheck(False)

5. @cython.wraparound(False)

関数ヘッダー

6. def mandelbrot_set_cython(double xmin, double xmax, double ymin, double ymax, int width, int height, int max_iter):

型指定と変数定義

7. cdef:

8. double[:] x = np.linspace(xmin, xmax, width)

9. double[:] y = np.linspace(ymin, ymax, height)

10. int[:, :] mandelbrot_image = np.zeros((height, width), dtype=np.int32)

11. int i, j, n

12. double zx, zy, zx2, zy2, cx, cy