如何加速numpy数组的分形生成？

Question

这是我写的一个小脚本，用于使用牛顿法来制作分形。

import numpy as np 
import matplotlib.pyplot as plt 

f = np.poly1d([1,0,0,-1]) # x^3 - 1 
fp = np.polyder(f) 

def newton(i, guess): 
    if abs(f(guess)) > .00001: 
     return newton(i+1, guess - f(guess)/fp(guess)) 
    else: 
     return i 

pic = [] 
for y in np.linspace(-10,10, 1000): 
    pic.append([newton(0,x+y*1j) for x in np.linspace(-10,10,1000)]) 

plt.imshow(pic) 
plt.show() 

我使用numpy的阵列，但仍然通过循环中的每个元素1000由-1000 linspaces应用newton()功能，其作用于一个单一的猜测，而不是整个阵列。

我的问题是这样的：如何改变我的方法来更好地利用numpy数组的优点？

P.S .:如果你想不用等太久就试试代码，最好使用100×100。

额外背景：
查看牛顿法找出多项式的零点。
分形的基本思想是测试复平面中的猜测并计算迭代次数以收敛到零。这就是newton()递归的结果，最终返回步数。在复平面中的猜测表示图片中的一个像素，按收敛步骤的数量着色。从一个简单的算法，你会得到这些美丽的分形。

Answer 1

我从Lauritz五Thaulow的代码工作，并能得到相当显著加速与以下代码：

import numpy as np 
import matplotlib.pyplot as plt 
from itertools import count 

def newton_fractal(xmin, xmax, ymin, ymax, xres, yres): 
    yarr, xarr = np.meshgrid(np.linspace(xmin, xmax, xres), \ 
          np.linspace(ymin, ymax, yres) * 1j) 
    arr = yarr + xarr 
    ydim, xdim = arr.shape 
    arr = arr.flatten() 
    f = np.poly1d([1,0,0,-1]) # x^3 - 1 
    fp = np.polyder(f) 
    counts = np.zeros(shape=arr.shape) 
    unconverged = np.ones(shape=arr.shape, dtype=bool) 
    indices = np.arange(len(arr)) 
    for i in count(): 
     f_g = f(arr[unconverged]) 
     new_unconverged = np.abs(f_g) > 0.00001 
     counts[indices[unconverged][~new_unconverged]] = i 
     if not np.any(new_unconverged): 
      return counts.reshape((ydim, xdim)) 
     unconverged[unconverged] = new_unconverged 
     arr[unconverged] -= f_g[new_unconverged]/fp(arr[unconverged]) 

N = 1000 
pic = newton_fractal(-10, 10, -10, 10, N, N) 

plt.imshow(pic) 
plt.show() 

对于N = 1000，我得到11.1秒的时间使用Lauritz的代码和使用这种代码1.7秒的时间。

这里有两个主要的加速。首先，我使用meshgrid来加速numpy输入值数组的创建。当N = 1000时，这实际上是加速的一个非常重要的部分。

第二次加速来自只对未收敛部分进行计算。 Lauritz在意识到他们正在放慢速度之前，正在使用屏蔽阵列。我在相当一段时间没有看过他们，但我确实记得过去掩盖的阵列是缓慢的来源。我相信这是因为它们在很大程度上是用纯Python实现的，而不是像numpy数组那样几乎完全用C写成。

Answer 2

这是我的刺伤。它快了16倍。

import numpy as np 
import matplotlib.pyplot as plt 
from itertools import count 

def newton_fractal(xmin, xmax, ymin, ymax, xres, yres): 
    arr = np.array([[x + y * 1j for x in np.linspace(xmin, xmax, xres)] \ 
     for y in np.linspace(ymin, ymax, yres)], dtype="complex") 
    f = np.poly1d([1,0,0,-1]) # x^3 - 1 
    fp = np.polyder(f) 
    counts = np.zeros(shape=arr.shape) 
    for i in count(): 
     f_g = f(arr) 
     converged = np.abs(f_g) <= 0.00001 
     counts[np.where(np.logical_and(converged, counts == 0))] = i 
     if np.all(converged): 
      return counts 
     arr -= f_g/fp(arr) 

pic = newton_fractal(-10, 10, -10, 10, 100, 100) 

plt.imshow(pic) 
plt.show() 

我不是numpy的专家，我相信那些能够optimalize它多一些，但已经是在速度方面的一个巨大的进步。

编辑：原来蒙面阵列并没有在所有帮助，去除它们导致了15％的速度增长，所以我已经移除上述溶液蒙面阵列。任何人都可以解释为什么蒙面阵列没有帮助？

Answer 3

我矢量化牛顿函数，并获得约。 85倍的速度与200×200点，500×500点，144倍的速度，并与1000×1000点的148倍的速度：

import numpy as np 
import matplotlib.pyplot as plt 

f = np.poly1d([1,0,0,-1]) # x^3 - 1 
fp = np.polyder(f) 
def newton(i, guess):    
    a = np.empty(guess.shape,dtype=int) 
    a[:] = i 
    j = np.abs(f(guess))>.00001 
    if np.any(j):   
     a[j] = newton(i+1, guess[j] - np.divide(f(guess[j]),fp(guess[j])))   
    return a 

npts = 1000 
x = np.linspace(-10,10,npts) 
y = np.linspace(-10,10,npts) 
xx, yy = np.meshgrid(x, y) 
pic = np.reshape(newton(0,np.ravel(xx+yy*1j)),[npts,npts]) 
plt.imshow(pic) 
plt.show() 

Answer 4

好的，我已经解决了Justin Peel的代码中的无限循环，在代码中添加了最大迭代条件，现在代码会绘制像z^4-1这样的多项式，并且它不会进入无限循环。如果有人知道如何改善这个错误，请告诉我们。我的解决方案，可能会使代码的执行速度变慢，但它起作用。 这是代码：

#!/usr/bin/python 
    # -*- coding: utf-8 -*- 

    import numpy as np 
    import itertools 
    import matplotlib.pyplot as plt 

    __author__ = 'Tobal' 
    __version__ = 1.0 


    def newton_fractal(f, xmin, xmax, ymin, ymax, xres, yres, tolerance, maxiters): 
     yarr, xarr = np.meshgrid(np.linspace(xmin, xmax, xres), np.linspace(ymin, ymax, yres) * 1j) 
     arr = yarr + xarr 
     ydim, xdim = arr.shape 
     arr = arr.flatten() 
     fp = np.polyder(f, m=1) 
     counts = np.zeros(shape=arr.shape) 
     unconverged = np.ones(shape=arr.shape, dtype=bool) 
     indices = np.arange(len(arr)) 
     iters = 0 
     for i in itertools.count(): 
      f_g = f(arr[unconverged]) 
      new_unconverged = np.abs(f_g) > tolerance 
      counts[indices[unconverged][~new_unconverged]] = i 
      if not np.any(new_unconverged) or iters >= maxiters: 
       return counts.reshape((ydim, xdim)) 
      iters += 1 
      unconverged[unconverged] = new_unconverged 
      arr[unconverged] -= f_g[new_unconverged]/fp(arr[unconverged]) 


    pic = newton_fractal(np.poly1d([1., 0., 0., 0., -1.]), -10, 10, -10, 10, 1000, 1000, 0.00001, 1000) 
    plt.imshow(pic, cmap=plt.get_cmap('gnuplot')) 
    plt.title(r'$Newton^{\prime} s\;\;Fractal\;\;Of\;\;P\,(z) =z^4-1$', fontsize=18, y=1.03) 
    plt.tight_layout() 
    plt.show() 

我使用Pycharm 5阿纳康达Python 3中和此IDE报告中的代码不np.any（new_unconverged）

警告应为“类型Union [ndarray，iterable]'，得到'bool'而不是

此警告也出现在原始Justin Peel代码中;我不知道如何解决它。我对这个问题很感兴趣。