创建在Python中设置Mandelbrot的性能提高

import numpy as np
from PIL import Image

z = 0

real_axis = np.linspace(-2,1,num=3000)
imaginary_axis = np.linspace(1,-1,num=2000)

complex_grid = [[complex(np.float64(a),np.float64(b)) for a in real_axis] for b in imaginary_axis]

pixel_grid = np.zeros((2000,3000,3),dtype=np.uint8)

for complex_list in complex_grid:
    for complex_number in complex_list:
        for iteration in range(255):
            z = z**2 + complex_number
            if (z.real**2+z.imag**2)**0.5 > 2:
                pixel_grid[complex_grid.index(complex_list),complex_list.index(complex_number)]=[iteration,iteration,iteration]
                break
            else:
                continue
        z = 0

#1 楼

我将重用最近在代码审查上发布在此处的答案的某些部分。蟒蛇。特别是多个嵌套循环。

NumPy可以帮助向量化您的代码，即在这种情况下，更多
循环在C后端而不是在Python
解释器中完成。我强烈建议听一听演讲
失去循环：Jake的NumPy进行快速数值计算。
VanderPlas。复杂的网格，后跟用于迭代网格的嵌套循环，当交给Python解释器时，图像很慢。幸运的是，NumPy可以减轻您的很多负担。

例如

 real_axis = np.linspace(-2, 1, num=3000)
imaginary_axis = np.linspace(1, -1, num=2000)
complex_grid = [[complex(np.float64(a),np.float64(b)) for a in real_axis] for b in imaginary_axis]
  
 
可以成为
 
 
 n_rows, n_cols = 2000, 3000
complex_grid_np = np.zeros((n_rows, n_cols), dtype=np.complex)
real, imag = np.meshgrid(real_axis, imaginary_axis)
complex_grid_np.real = real
complex_grid_np.imag = imag
 
 
 
没有循环，只是普通的简单的NumPy。
 
 
也一样
 
 
 for complex_list in complex_grid:
    for complex_number in complex_list:
        for iteration in range(255):
            z = z**2 + complex_number
            if (z.real**2+z.imag**2)**0.5 > 2:
                pixel_grid[complex_grid.index(complex_list),complex_list.index(complex_number)]=[iteration,iteration,iteration]
                break
            else:
                continue
        z = 0
 
 
 
可以转换为
 
 
 z_grid_np = np.zeros_like(complex_grid_np)
elements_todo = np.ones((n_rows, n_cols), dtype=bool)
for iteration in range(255):
    z_grid_np[elements_todo] = \
        z_grid_np[elements_todo]**2 + complex_grid_np[elements_todo]
    mask = np.logical_and(np.absolute(z_grid_np) > 2, elements_todo)
    pixel_grid_np[mask, :] = (iteration, iteration, iteration)
    elements_todo = np.logical_and(elements_todo, np.logical_not(mask))
 
 
 
，这只是一个循环，而不是三个嵌套循环。在这里，与您处理break情况一样，还需要更多技巧。如果elements_todo值未标记为已完成，则仅用于对其进行更新。没有这个，也许还有更好的解决方案。
 
我添加了以下几行
 
 
 z 
 
 
可以根据您的参考实现验证我的结果。
 
在我测试的几种complex_grid_close = np.allclose(np.array(complex_grid), complex_grid_np)
pixel_grid_close = np.allclose(pixel_grid, pixel_grid_np)
print("Results were similar: {}".format(all((complex_grid_close, pixel_grid_close))))
组合上，矢量化代码在我的机器上的速度快9-10倍。例如。适用于n_rows/n_cols：
 
 
 n_rows, n_cols = 1000, 1500 
 
 
减小尺寸
 
在查看代码时，我以某种方式略微忽略的一个方面是，由于所有颜色通道值都相同，因此您实际上是在创建灰度图像。考虑到这一点，您可以轻松地将程序处理的数据大小从Looped generation took 61.989842s
Vectorized generation took 6.656926s
Results were similar: True
减少到3000x2000x3。尽管我不是该领域的专家，但这可能会帮助您的程序提高缓存效率。
 
 
编辑/附录：更多时间安排
 
包括trichoplax在他的答案和Peter Cordes在这样的评论中建议的“无平方根”优化
 
 
 3000x2000 
 /> 
将为您提供mask = np.logical_and((z_grid_np.real**2+z_grid_np.imag**2) > 4, elements_todo)
的另一个半秒，即大约是原始解决方案速度的12倍。 / pre> 
 
 Reinderien关于Mandelbrot集对称性的提示的快速实现将再次为此增加大约2倍。
 
 
 n_rows, n_cols = 1000, 1500 
 
 
但是，与最初的相比，我快速的黑客方法并没有导致输出完全在10 loops, best of 5: 4.98 s per loop
10 loops, best of 5: 4.28 s per loop (in grayscale, 14x)
的容差范围内。有趣的是，它似乎在单个像素处相差一个，但在视觉上仍然相同。由于这篇文章已经很长了，我将把重新实现作为练习留给读者。
 
根据您的需求，除所有前面提到的优化。
 
 
 10 loops, best of 5: 2.54 s per loop (~24x)
10 loops, best of 5: 2.07 s per loop (in grayscale, ~30x)
 
 

评论
\ $ \ begingroup \ $

我很难理解您发布的基于numpy的代码，可能是因为有一些新命令，但是我将尝试矢量化编程！

 \ $ \ endgroup \ $

 –伊恩

 19 Mar 27 '19在6:24 
\ $ \ begingroup \ $

您能否缩小代码的哪一部分/什么命令使您头痛？然后，我将在需要的地方添加更多详细信息。

 \ $ \ endgroup \ $

 – AlexV

 19 Mar 27 '19在7:02 
\ $ \ begingroup \ $

 @Peter Cordes：不需要为实部和虚部分别设置阵列。如果要手动执行，可以使用arr.real和arr.imag来访问这两个部分。

 \ $ \ endgroup \ $

 – AlexV

 19 Mar 27 '19在8:05 
\ $ \ begingroup \ $

如果w * h乘以每个复数双16字节大于32k（L1d大小）或256k（典型的L2高速缓存大小），也值得一看缓存块：重复遍历整个数组的一部分。 （例如1500 * 1000 * 16 = 24MB，仅适合大型Xeon上的L3缓存，或者完全不适合普通台式机CPU。）1500 * 16B = 24kB，因此重复循环超过1行可能是一个胜利。 （或正如其他答案所指出的那样，问题的不同区域具有不同的典型迭代次数，因此当所有像素均达到| m |> 4时，在正方形拼贴中进行操作可能会让您在几次迭代之后停止

 \ $ \ endgroup \ $

 – Peter Cordes

 19 Mar 27 '19在8:07 
\ $ \ begingroup \ $

这相当复杂，恕不恕我直言，因为这种代码审查和OP似乎具有的经验显然不突出。也许我们可以在其他地方继续讨论？

 \ $ \ endgroup \ $

 – AlexV

 19 Mar 27 '19在8:19 

#2 楼

这将涵盖性能以及Python样式。

将常量保存在一个地方

您当前拥有神奇的数字2000和3000，即图像的分辨率。将它们保存到名为X，Y或W，H的变量中。

提及您的需求

您不仅依赖于Python 3和Jupyter-您还依赖numpy和枕头。如果您还没有的话，这些文件应该放在requirements.txt文件中。

根本不保存您的复杂网格

。 complex_number应该基于range表达式在循环中动态形成。免责声明：如果您要向量化（您应该这样做），那么情况恰好相反-您将保留复杂的网格，并且会丢失一些网格

不要使用index查找

您正在使用index来获取坐标。不要这样做-也要在循环中形成坐标。

Mandelbrot是对称的

请注意，它是镜像的。这意味着您可以将计算时间减半，并将每个像素保存到顶部和底部一半。

我将展示一些示例代码，其中包含上述所有建议。只需（几乎）执行@Alex所说的，我就可以在实现过程中完成一半的工作，但有一个区别：容纳我描述的对称性优化。

#3 楼

Mandelbrot特定的优化

这些可以与其他答案中的Python特定的优化结合。

避免冗余平方根
br />
等同于

if (z.real**2+z.imag**2)**0.5 > 2:

（只需将原始比较的两边都平方即可得到优化的比较）除非使用平方，否则请避免平方。

从原点超过2的任何点将继续向无穷远移动。因此，检查点是否已超出半径2的圆或是否已超出完全包含该圆的其他有限形状并不重要。例如，检查点是否在正方形而不是圆形之外，可以避免对实部和虚部进行平方。这也意味着您将需要更多的迭代，但是数量很少，因此每次迭代都要更快，这会被忽略。例如， />可以替换为

if z.real ** 2 + z.imag ** 2 > 4:

此建议的例外情况是，如果圆对您的输出很重要。如果仅将集合内的点绘制为黑色，将集合外的点绘制为白色，则两种方法的图像都相同。但是，如果计算一个点要逃逸的迭代次数，并用它来确定集合外点的颜色，则彩色条纹的形状在正方形边界处会不同于圆形边界处。布景的内部是相同的，但外面的颜色将以不同的形状排列。

在您的示例图像中，几乎看不到彩色条纹，大多数外部和内部都是黑色的。在这种情况下，我怀疑使用这种优化方法会在外观上产生重大差异。但是，如果将来更改为显示更宽的条纹，则可能需要取消此优化（取决于您想要的外观）。

尽可能对内部进行硬编码

集合的内部比外部花费更长的时间进行计算。内部的每个像素保证进行255次迭代（如果为更高质量的图像增加最大迭代次数，则保证迭代次数更多），而外部的每个像素花费的迭代次数少于此次数。绝大多数外部像素仅需进行几次迭代。

如果您希望代码可适应于放大到任意位置，那么您将不知道图像的哪些部分将成为内部要点。但是，如果只希望此代码生成整个图像集，则可以避免计算已知的内部像素，从而显着提高速度。例如，如果您检查像素是在主心形中还是在大圆圈中，则可以为所有这些像素分配255的迭代计数，而无需实际进行任何迭代。您增加的最大迭代次数越多，就值得将其排除在外的圆圈越多，因为平均外部像素与平均内部像素之间的计算时间差将继续显着不同。

I不知道这些圆的确切中心和半径，也不知道心形的确切方程式，但是选择不与外部重叠的粗略估计仍然会对速度产生很大的影响。甚至排除一些眼睛选择的完全在内部的正方形也会有所帮助。

#4 楼

不要使用矢量化的numpy，而应使用numba jit

使用numpy计算Mandelbrot集并不是一个很好的选择，因为相同的数据将被反复存储和从内存中加载和加载到内存中，从而破坏了缓存。更好的选择是使用jit编译器来加速关键代码路径，例如numba jit。在这种情况下，4个字符可使函数运行速度提高200倍。使用功能注释@jit时，此代码在3000x2000分辨率下可在2秒内运行，而不是400秒，而没有任何特殊技巧。距离估算值可以使外观更平滑：

#5 楼

我不是python专家。我对Mandlebrot的生成非常满意（我花了很多时间在自定义的Julia Set生成器上。）迭代。忘记干净代码或好的OOP原则。对于类似这样的重复项，您希望尽可能的细腻。 >
仅从第一行开始，想象一下内存中幕后发生的事情。您有一个复数的实例。您想对其平方...因此它必须创建另一个复杂对象实例以保存平方值。然后，向其添加另一个复数-这意味着您将创建另一个Complex实例以保存加法的结果。

您将左右创建对象实例，然后正在按3000 x 2000 x 255倍的顺序进行操作。创建多个类实例听起来并不多，但是当您进行十亿次操作时，它会把事情拖累。

z = z**2 + complex_number
if (z.real**2+z.imag**2)**0.5 > 2:
    pixel_grid[complex_grid.index(complex_list),complex_list.index(complex_number)]=[iteration,iteration,iteration]
    break
else:
    continue

没有对象被创建和销毁。归结为尽可能高效。看起来不那么好看……但是当您做某事十亿次时，即使节省一百万分之一秒也可以节省15分钟。简化逻辑，以便您不必进行平方根运算-如果您可以在渐变方式上进行一些细微的改动，请在边界框中使用“是x或y”更改“幅度”检查“

#6 楼

您可以使用一些技巧来使Mandelbrot渲染器真正飞行。进入一个周期。我发现，最经济的检测方法是进行x次迭代，测试是否与以前相同，然后递增x并重复。

先绘制一个半分辨率版本

您的情况就是一张1000x1500的图片。计算每个像素代表真实图像中的一个像素。然后，如果一个像素完全被具有相同迭代计数的其他像素包围，则可以假定它也具有该迭代计数，并跳过计算。大量的时间。每当您计算不可跳过的像素以查找以前可能被视为可跳过但不可跳过的其他像素时，也应使用泛洪填充样式算法。这应该可以解决大多数问题。

还要注意，这是递归的。在计算1000x1500版本之前，您应该先计算500x750版本，然后再计算250x375版本等。精度很高，这可能会浪费大量的计算时间。但是，严格来说，您只需要对一个像素使用高精度。

我们从位置\ $ p_0 \ $开始，并遵循通常的迭代公式：

\ $ p_ {x + 1} = {p_x} ^ 2 + p_0 \ $

我们将\ $ p_x \ $的所有值记录为常规的双精度复数。现在我们计算\ $ q \ $，但是通过计算\ $ d \ $来完成，其中\ $ d_x = q_x-p_x \ $：

\ $ d_ {x + 1} = 2d_xp_x + {d_x} ^ 2 +（q_0-p_0）\ $

这有点复杂，但是我们只需要使用双精度数，因此深度缩放时它快得多。

一个问题是\ $ p \ $序列必须至少与\ $ q \ $序列一样长，而且我们无法提前告知最佳的\ $ p \ $序列。实际上，我们经常需要使用高精度算术来计算新的\ $ p \ $序列，因为我们发现了具有更长转义时间的像素。在它周围，Python很慢。 NumPy可以完成繁重的工作，可以大大加快速度，但是与用C编写的相同代码相比，这很尴尬。我建议学习使用Ctypes并编写一个小的C库以遵循迭代公式。

#7 楼

要进行更多优化，您可以深入研究Fractint的源代码。它是在80年代末/ 90年代初编写的，其硬件比现代CPU慢数千倍（但可以在不到一分钟的时间内生成640x480瓦片的图像）。 FractInt的“重大交易”之一是，大多数实现都使用整数数学来实现定点算术（当浮点数由慢速库或通过可选的昂贵的外部芯片来模拟时，则要大得多。 （请参阅Intel 8087至80387）。

另一个改进之处：将图像分成正方形。由于Mandelbrot集已连接，因此如果正方形在其边界上没有该集的点，则其内部也没有该集的点。这导致边缘跟随作为一种策略，可以大大减少必须实际计算的像素数。

源以及MS-DOS和Win 3.x可执行文件仍然可用。 >