Mandelbrot Numba/Numpy 矢量化? [英] Mandelbrot Numba/Numpy vectorization?
问题描述
我使用 kivy 在 Python 中编写了一个交互式 mandelbrot 渲染器，您可以在其中使用鼠标指针进行缩放，并正在尽我所能对其进行优化.我目前使用这个实现来渲染设置/缩放(这是一个小片段，只是用来渲染它的两个函数):
将 numba 导入为 nb将 numpy 导入为 np@nb.njit(缓存=真，并行=真)def mandelbrot(c_r, c_i,maxIt): #mandelbrot 函数z_r = 0z_i = 0z_r2 = 0z_i2=0对于 nb.prange(maxIt) 中的 x:z_i = 2 * z_r * z_i + c_iz_r = z_r2 - z_i2 + c_rz_r2 = z_r * z_rz_i2 = z_i * z_i如果 z_r2 + z_i2 >4:返回 x返回最大值@nb.njit(缓存=真，并行=真)def DrawSet(W, H, xStart, xDist, yStart, yDist, maxIt):array = np.zeros((H, W, 3), dtype=np.uint8) #array 为每个像素保存 'hsv' 元组对于 nb.prange(0, W) 中的 x:c_r = (x/W)* xDist + xStart #一些计算实部的数学对于范围内的 y (0, H):c_i = -((y/H) * yDist + yStart) #计算虚部的更多数学运算cIt = mandelbrot(c_r, c_i, maxIt)颜色 = int((255 * cIt)/maxIt)array[y,x] = (color, 255, 255) #添加色调值return array #returns hsv 数组，稍后使用 PIL 显示我目前的表现相当不错.它可以在大约 0.08 - 0.09 秒内渲染一个 500x500 的区域，其中每个点都有界(所以基本上是黑色图片，最坏的情况)，迭代 300 次.我将 Numba JIT 与并行范围函数prange()"一起使用，这有很大帮助.
但是，我听说矢量化通常是渲染此类分形的最快方法.经过大量研究(我对矢量化很陌生)，我设法将这个实现放在一起:
将 numba 导入为 nb将 numpy 导入为 npdef DrawSet(W, H, xStart, xEnd, yStart, yEnd, maxIt):array = np.zeros((H,W,3), dtype = np.uint8) # 3D 数组包含每个像素的 'hsv' 元组 (hue,saturation,value)x = np.linspace(xStart, xEnd, W).reshape((1, W)) #将水平像素缩放到 x 轴y = np.linspace(yStart, yEnd, H).reshape((H, 1)) #将垂直像素缩放到y轴c = x + 1j * y #从x轴(实数)和y轴(虚数)创建复平面z = np.zeros(c.shape, dtype= np.complex128)div_time = np.zeros(z.shape, dtype= int)m = np.full(c.shape, True, dtype= bool)div_time = loop(z, c, div_time, m, maxIt)array[:,:,0] = (div_time/maxIt) * 255 -20 #adding 'hue' 值array[:,:,1] = 255 #添加饱和度"值数组[:,:,2] = 255 #添加值"返回数组@nb.vectorize(nb.int64[:,:](nb.complex128[:,:], nb.complex128[:,:], nb.int64[:,:], nb.boolean[:,:],nb.int64))def循环(z，c，div_time，m，maxIt):对于我在范围内(maxIt):z[m] = z[m]**2 + c[m]发散 = np.greater(np.abs(z), 2, out=np.full(c.shape, False), where=m)div_time[发散] = im[np.abs(z) >2] = 错误返回 div_time没有@nb.vectorize 装饰器，它运行得非常慢.(500x500 的最坏情况为 4 秒，300 It.).使用 @nb.vectorize 装饰器，我收到此错误:
<预><代码>回溯(最近一次调用最后一次):文件MandelBrot.py"，第 13 行，在 <module> 中.从测试导入 DrawSet文件C:\Users\User\Documents\Code\Python\Mandelbrot-GUI\test.py"，第 25 行，在 <module> 中.def循环(z，c，div_time，m，maxIt):文件C:\Users\User\AppData\Local\Programs\Python\Python38\lib\site-packages\numba\np\ufunc\decorators.py"，第 119 行，包装中在 ftylist 中签名:类型错误:签名"对象不可迭代我做错了什么?我是否以正确的方式定义了所有的 numba 签名?这种矢量化方法甚至会超过我当前的实现吗?
我会感谢每一个建议！提前致谢.
您的实现已经矢量化了！
矢量化的想法是创建通用函数，以元素方式对其进行操作数组.您只需定义对单个元素执行的操作，向量化机制将允许使用数组调用您的函数.
该函数计算单个点 c:
def mandelbrot_point(c, max_it):z = 0j对于我在范围内(max_it):z = z**2 + c如果 abs(z) >2:返回我返回 0您可以使用 Numpy 对其进行矢量化:
@np.vectorizedef mandelbrot_numpy(c, max_it):z = 0j对于我在范围内(max_it):z = z**2 + c如果 abs(z) >2:返回我返回 0或者您可以使用 Numba 对其进行矢量化.请注意，函数的签名描述了如何处理单个点:
@nb.vectorize([nb.int64(nb.complex128, nb.int64)])def mandelbrot_numba(c, max_it):z = 0j对于我在范围内(max_it):z = z**2 + c如果 abs(z) >2:返回我返回 0然后您可以使用任意维数的标量或数组调用向量化函数:
<预><代码>>>>p = 0.4+0.4j>>>mandelbrot_point(p, 99)8>>>mandelbrot_numpy(p, 99)数组(8)>>>mandelbrot_numba(p, 99)8>>>x = np.linspace(-2, 2, 11)>>>mandelbrot_numpy(x, 99)数组([0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1])>>>mandelbrot_numba(x, 99)数组([0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1])>>>x = np.atleast_2d(x)>>>y = x.T>>>c = x + 1j * y>>>mandelbrot_numpy(c, 99)数组([[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],[ 0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1],[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])>>>mandelbrot_numba(c, 99)数组([[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],[ 0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1],[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])Numpy 的 vectorize 极大地简化了您的代码，但作为文档说它主要是为了方便，而不是为了性能.实现本质上是一个 for 循环.
根据我的测量，Numpy 向量化版本比您的原始实现略快，而 Numba 向量化版本快一个数量级.
I wrote an interactive mandelbrot renderer in Python using kivy where you can zoom with the mouse pointer and am trying to optimize it as best as I can. I currently use this implementation to render the set/zoom (this is a small snippet, just the two functions used to render it):
import numba as nb
import numpy as np
@nb.njit(cache= True, parallel = True)
def mandelbrot(c_r, c_i,maxIt): #mandelbrot function
z_r = 0
z_i = 0
z_r2 = 0
z_i2= 0
for x in nb.prange(maxIt):
z_i = 2 * z_r * z_i + c_i
z_r = z_r2 - z_i2 + c_r
z_r2 = z_r * z_r
z_i2 = z_i * z_i
if z_r2 + z_i2 > 4:
return x
return maxIt
@nb.njit(cache= True, parallel = True)
def DrawSet(W, H, xStart, xDist, yStart, yDist, maxIt):
array = np.zeros((H, W, 3), dtype=np.uint8) #array that holds 'hsv' tuple for every pixel
for x in nb.prange(0, W):
c_r = (x/W)* xDist + xStart #some math to calculate real part
for y in range (0, H):
c_i = -((y/H) * yDist + yStart) #some more math to calculate imaginary part
cIt = mandelbrot(c_r, c_i, maxIt)
color = int((255 * cIt) / maxIt)
array[y,x] = (color, 255, 255) #adds hue value
return array #returns hsv array, gets later displayed using PIL
I currently get pretty good performance. It can render a 500x500 area where every point is bounded (so basically a black picture, the worst-case) in about 0.08 - 0.09 seconds, with 300 iterations. I'm using Numba JIT along with the parallelized range function 'prange()', which helps quite a lot.
However, I've heard that vectorization is usually the fastest for rendering these kind of fractals. After a lot of research (I'm pretty new to vectorization), I managed to put this implementation together:
import numba as nb
import numpy as np
def DrawSet(W, H, xStart, xEnd, yStart, yEnd, maxIt):
array = np.zeros((H,W,3), dtype = np.uint8) # 3D array containing 'hsv' tuple (hue,saturation,value) of each pixel
x = np.linspace(xStart, xEnd, W).reshape((1, W)) #scaling horizontal pixels to x-axis
y = np.linspace(yStart, yEnd, H).reshape((H, 1)) #scaling vertical pixels to y-axis
c = x + 1j * y #creating complex plane out of x axis (real) and y axis (imaginary)
z = np.zeros(c.shape, dtype= np.complex128)
div_time = np.zeros(z.shape, dtype= int)
m = np.full(c.shape, True, dtype= bool)
div_time = loop(z, c, div_time, m, maxIt)
array[:,:,0] = (div_time/maxIt) * 255 -20 #adding 'hue' value
array[:,:,1] = 255 #adding 'saturation' value
array[:,:,2] = 255 #adding 'value'
return array
@nb.vectorize(nb.int64[:,:](nb.complex128[:,:], nb.complex128[:,:], nb.int64[:,:], nb.boolean[:,:], nb.int64))
def loop(z, c, div_time, m, maxIt):
for i in range(maxIt):
z[m] = z[m]**2 + c[m]
diverged = np.greater(np.abs(z), 2, out=np.full(c.shape, False), where=m)
div_time[diverged] = i
m[np.abs(z) > 2] = False
return div_time
Without, the @nb.vectorize decorator, it runs horribly slow. (4 seconds worst-case for 500x500, 300 It.). With the @nb.vectorize decorator, I'm getting this error:
Traceback (most recent call last):
File "MandelBrot.py", line 13, in <module>
from test import DrawSet
File "C:\Users\User\Documents\Code\Python\Mandelbrot-GUI\test.py", line 25, in <module>
def loop(z, c, div_time, m, maxIt):
File "C:\Users\User\AppData\Local\Programs\Python\Python38\lib\site-packages\numba\np\ufunc\decorators.py", line 119, in wrap
for sig in ftylist:
TypeError: 'Signature' object is not iterable
What am I doing wrong? Am I defining all the numba signatures the correct way? Would this vectorization method even top my current implementation?
I would be thankful for every advice! Thank you in advance.
Your implementation is already vectorized!
The idea of vectorization is to create universal functions that operate elementwise on arrays. You only need to define what to do on a single element and the vectorization mechanism will allow your function to be called using arrays.
This function does the calculation for a single point c:
def mandelbrot_point(c, max_it):
z = 0j
for i in range(max_it):
z = z**2 + c
if abs(z) > 2:
return i
return 0
You can vectorize it using Numpy:
@np.vectorize
def mandelbrot_numpy(c, max_it):
z = 0j
for i in range(max_it):
z = z**2 + c
if abs(z) > 2:
return i
return 0
Or you can vectorize it using Numba. Notice that the function's signature describes what to do with a single point:
@nb.vectorize([nb.int64(nb.complex128, nb.int64)])
def mandelbrot_numba(c, max_it):
z = 0j
for i in range(max_it):
z = z**2 + c
if abs(z) > 2:
return i
return 0
Then you can call the vectorized functions using scalars or arrays of any number of dimensions:
>>> p = 0.4+0.4j
>>> mandelbrot_point(p, 99)
8
>>> mandelbrot_numpy(p, 99)
array(8)
>>> mandelbrot_numba(p, 99)
8
>>> x = np.linspace(-2, 2, 11)
>>> mandelbrot_numpy(x, 99)
array([0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1])
>>> mandelbrot_numba(x, 99)
array([0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1])
>>> x = np.atleast_2d(x)
>>> y = x.T
>>> c = x + 1j * y
>>> mandelbrot_numpy(c, 99)
array([[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],
[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],
[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],
[ 0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1],
[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],
[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],
[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],
[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])
>>> mandelbrot_numba(c, 99)
array([[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],
[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],
[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],
[ 0, 0, 0, 0, 0, 0, 6, 2, 1, 1, 1],
[ 0, 2, 6, 6, 0, 0, 8, 2, 1, 1, 0],
[ 0, 2, 2, 3, 5, 17, 3, 1, 1, 1, 0],
[ 0, 1, 1, 2, 2, 2, 1, 1, 1, 1, 0],
[ 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])
Numpy's vectorize greatly simplifies your code, but as the docs say it is provided primarily for convenience, not for performance. The implementation is essentially a for loop.
According to my measurements, the Numpy-vectorized version is slightly faster than your original implementation, and the Numba-vectorized version is an order of magnitude faster.
这篇关于Mandelbrot Numba/Numpy 矢量化?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
