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() 

Iは、単一の推測ではなく、アレイ全体に作用newton()関数を適用する1000によって-1000 linspacesの各要素を介してnumpyのアレイを使用して、それにもかかわらず、ループしています。

私の質問はこれですnumpy配列の利点をより良く活用するために私のアプローチを変更するにはどうすればよいですか？

P.S.：あまりにも長く待たずにコードを試したい場合は、100×100の方が良いでしょう。

追加の背景：
多項式のゼロを見つける方法を参照してください。
フラクタルの基本的な考え方は、複素平面の推測をテストし、反復の回数を数えてゼロに収束させることです。これは、最終的にステップ数を返すnewton()の再帰についてです。複素平面内の推測は、収束までのステップ数で表される画像内のピクセルを表す。簡単なアルゴリズムから、これらの美しいフラクタルを得ることができます。

Answer 1

私はLauritz V.タヴロウのコードから働いて、とかなり大幅なスピードアップを取得することができました次のコード：N = 1000の場合

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() 

、私はこのコードを使用してLauritzのコードと1.7秒の時間を使って11.1秒の時間を取得します。

ここには2つのスピードアップがあります。まず、numgpyの入力値の配列の作成を高速化するためにmeshgridを使用しました。 N = 1000のとき、これは実際にスピードアップのかなりの部分です。

2番目のスピードアップは、収束していない部分の計算だけです。 Lauritzは、マスクされた配列を使用して、遅くなっていることを認識しました。私はかなりの時間でそれらを見ていないが、私はマスクされた配列が過去の遅さの源であることを覚えている。純粋なPythonではnumpy配列のようにCで書かれているのではなく、numpy配列で実装されているからです。

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ことができると確信していますが、すでにそれは速度の点で大きな改善です。

EDIT：それはそれらを取り除く、マスクされた配列はまったく助けにはならなかったが判明した15％の速度増加になったので、私は上記の溶液からマスクされた配列を削除しました。マスクされた配列がなぜ役に立たなかったのか誰でも説明できますか？

Answer 3

私はニュートン関数をベクトル化し、 85倍の速さ144倍の速さ200×200ポイントと500×500ポイント、148倍の速さで1000×1000ポイントを持つ：

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

私はジャスティン・ピールのコードで無限ループを解いて、コードで最大反復条件を追加しました。コードは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() 

私はアナコンダのPython 3でPycharm 5を使用していますし、このIDEは「タイプ予想np.any（new_unconverged）ないコード

に警告を報告します

この警告は、元のジャスティン・ピールのコードにも表示されています。私はそれを解決する方法を知らない。私はこの問題に非常に興味があります。