Cython手抜き実装のモンテカルロ法のスピードを速くする．

本日はCython

以前書いた
Pythonでモンテカルロ法の実装を手抜きで速くする．

の　Cython　を用いた実装があまりにも手抜きすぎたので補足しておきます．

元々のPythonのコードはこれ

import random
NUM=100000000
def monte():
    counter = 0
    for i in range(NUM):
        x = random.random()
        y = random.random()
        if x*x+y*y < 1.0:
            counter += 1
    pi = 4.0*counter/NUM
    print(pi)


def main():
    monte()

if __name__ == '__main__':
    main()

　手抜きCythonはこれ

#monte.pyx
import random

cdef int NUM = 100000000

cdef cmonte():
    cdef :
        int counter = 0
        int i=0
        double x
        double y
    for i in range(NUM):
        x = random.random()
        y = random.random()
        if x*x + y*y < 1.0:
            counter += 1

    cdef double pi = 4.0*counter/NUM
    return pi

def monte():
    pi=cmonte()
    print(pi)

100秒ほどから17秒（12inch MacBookで計測）に縮んだのはいいけれどもっと頑張れるはず．
randomの部分はPythonのコードなのでここをCの実装を使えば良さそう．

Cython 改良版

from libc.stdlib cimport rand, RAND_MAX

cdef int NUM = 100000000

def monte():
    cdef :
        int counter = 0
        int i=0
        double x
        double y
    for i in range(NUM):
        x = (rand()+1.0)/(RAND_MAX+2.0)
        y = (rand()+1.0)/(RAND_MAX+2.0)
        if x*x + y*y < 1.0:
            counter += 1

    pi = 4.0*counter/NUM
    print(pi)

setup.py も書きましょう．

#setup.py
from setuptools import setup, Extension

ext_modules = [
    Extension(
        name='monte',
        sources=['monte.pyx']
        )
]

setup(
    name = 'cymonte',
    ext_modules = ext_modules
)

次をターミナルで実行するとモジュールができます．Windowsでもいけると思います．

$ python setup.py build_ext --inplace

あとは動かすためのメインスクリプトを書く．

#main.py
import monte
monte.monte()

さあ，動かしましょう．

$ time python main.py
real    0m2.081s
user    0m1.935s
sys     0m0.070s

チューニングをしてかなり速くなりました！
Numbaの実装，　C++の実装に匹敵するほどまで近けることができましたね．

Reference:

Monte Carlo Simulation with Cython