スポンサーリンク
pythonでmandelbrotを描く¶
Mandelbrotをpythonで書くと以下の例で分かるように非常に遅い。
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import time
# color function for point at (x, y)
def mandel(x, y, max_iters):
c = complex(x, y)
z = 0.0j
for i in range(max_iters):
z = z*z + c
if z.real*z.real + z.imag*z.imag >= 4:
return i
return max_iters
def create_fractal(xmin, xmax, ymin, ymax, image, iters):
height, width = image.shape
pixel_size_x = (xmax - xmin)/width
pixel_size_y = (ymax - ymin)/height
for x in range(width):
real = xmin + x*pixel_size_x
for y in range(height):
imag = ymin + y*pixel_size_y
color = mandel(real, imag, iters)
image[y, x] = color
gimage = np.zeros((1024, 1536), dtype=np.uint8)
xmin, xmax, ymin, ymax = np.array([-2.0, 1.0, -1.0, 1.0]).astype('float32')
iters = 50
start = time.clock()
create_fractal(xmin, xmax, ymin, ymax, gimage, iters)
dt = time.clock() - start
print("Mandelbrot created on CPU in %f s" % dt)
plt.grid(False)
plt.imshow(gimage, cmap='jet')
pass
スポンサーリンク
cythonでmandelbrotを描く¶
%load_ext cython
%%cython
import matplotlib.pyplot as plt
import numpy as np
import time
cimport cython
cdef extern from "complex.h":
double cabs(double complex)
# color function for point at (x, y)
cdef unsigned char mandel_cython(double x, double y, int max_iters):
cdef double complex c, z
c = x + y*1j
z = 0.0j
for i in range(max_iters):
z = z*z + c
if cabs(z) >= 2:
return i
return max_iters
@cython.cdivision(True)
def create_fractal_cython(double xmin, double xmax, double ymin, double ymax, unsigned char[:, :] image, int iters):
cdef int x, y
cdef int height, width
cdef double pixel_size_x, pixel_size_y
cdef double real, imag
cdef unsigned char color
height = image.shape[0]
width = image.shape[1]
pixel_size_x = (xmax - xmin)/width
pixel_size_y = (ymax - ymin)/height
for x in range(width):
real = xmin + x*pixel_size_x
for y in range(height):
imag = ymin + y*pixel_size_y
color = mandel_cython(real, imag, iters)
image[y, x] = color
gimage = np.zeros((1024, 1536), dtype=np.uint8)
xmin, xmax, ymin, ymax = np.array([-2.0, 1.0, -1.0, 1.0]).astype('float32')
iters = 50
start = time.clock()
create_fractal_cython(xmin, xmax, ymin, ymax, gimage, iters)
dt = time.clock() - start
print("Mandelbrot created on CPU in %f s" % dt)
plt.grid(False)
plt.imshow(gimage, cmap='jet')
pass
224.312500/1.093750
cythonはpythonの205倍の速度でmandelbrotを描き上げることができた。
スポンサーリンク
スポンサーリンク