Как ускорить фрактальную генерацию с помощью массивов 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.

Второе ускорение происходит только при выполнении вычислений на несвязанных участках. Лориц использовал маскированные массивы для этого, прежде чем осознать, что они замедляют работу. Я не смотрел на них довольно долгое время, но я помню маскированные массивы, которые были источником медлительности в прошлом. Я считаю, что это происходит потому, что они в значительной степени реализованы в чистом 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 раз быстрее, с 200х200 точек, 144 раз быстрее с 500x500 точек, и 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() 

Я использую PyCharm 5 с Anaconda Python 3, и это IDE сообщает предупреждение в коде не np.any (new_unconverged)

ожидаемый тип» Union [ndarray, iterable] ', вместо этого получил «bool»

Это предупреждение появляется в исходном коде Justin Peel; и я не знаю, как его решить. Я очень заинтересован в этой проблеме.