Come velocizzare la generazione di frattali con gli array numpy?

Question

Ecco una piccola sceneggiatura che ho scritto per creare frattali usando il metodo di Newton.

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

sto usando array NumPy, ma comunque ciclicamente ogni elemento 1000-by-1000 linspaces applicare la funzione newton(), che agisce su un singolo indovinare e non un intero array.

La mia domanda è questa: Come posso modificare il mio approccio per sfruttare meglio i vantaggi degli array di numpy?

P.S .: Se si desidera provare il codice senza attendere troppo a lungo, è meglio utilizzare 100 per 100.

Sfondo aggiuntivo:
Vedere Metodo di Newton per trovare gli zeri di un polinomio.
L'idea di base del frattale è di verificare le ipotesi nel piano complesso e contare il numero di iterazioni per convergere a zero. Questo è il motivo della ricorsione in newton(), che in definitiva restituisce il numero di passaggi. Un'ipotesi nel piano complesso rappresenta un pixel nell'immagine, colorato dal numero di passaggi alla convergenza. Da un semplice algoritmo, ottieni questi splendidi frattali.

Answer 1

ho lavorato dal codice di Lauritz V. Thaulow ed era in grado di ottenere una bella significativa accelerazione con la codice seguente:

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

Per N = 1000, ho un tempo di 11,1 secondi utilizzando il codice di Lauritz e un tempo di 1,7 secondi con questo codice.

Ci sono due velocità principali qui. Innanzitutto, ho usato meshgrid per velocizzare la creazione della serie numpy dei valori di input. Questa è in realtà una parte piuttosto significativa dell'accelerazione quando N = 1000.

La seconda accelerazione deriva dal solo calcolo delle parti non convertite. Lauritz stava usando gli array mascherati per questo prima di rendersi conto che stavano rallentando le cose. Non li ho guardati da un po 'di tempo, ma ricordo che gli array mascherati sono stati una fonte di lentezza nel passato. Credo che sia perché sono in gran parte implementati in Python puro su una matrice numpy piuttosto che essere scritti quasi completamente in C come array numpy.

Answer 2

Ecco la mia pugnalata. È circa 16 volte più veloce.

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

Io non sono un esperto di NumPy, sono sicuro che quelli che sono possono optimalize ancora un po ', ma già si tratta di un enorme miglioramento in termini di velocità.

EDIT: Si è scoperto matrici mascherati non ha aiutato affatto, rimuovendoli comportato una velocità aumento del 15%, così ho rimosso gli array mascherati dalla soluzione sopra. Qualcuno può spiegare perché gli array mascherati non hanno aiutato?

Answer 3

Ho vettorializzato la funzione newton e ottenuto ca. 85 volte più veloci con 200x200 punti, 144 volte più veloce con 500x500 punti, e 148 volte più veloce con 1000x1000 punti:

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

Ok, ho risolto il ciclo infinito nel codice di Justin Peel aggiungendo una condizione di iterazioni massime nel codice, ora il codice traccia polinomi come z^4-1 e non entra in un ciclo infinito. Se qualcuno sa come migliorare questo bug, faccelo sapere.La mia soluzione, forse rende più lenta l'esecuzione del codice ma funziona. Questo è il codice:

#!/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() 

sto usando PyCharm 5 con Anaconda Python 3 e questo IDE segnala un avviso nel codice non np.any (new_unconverged)

Previsto 'tipo Union [ndarray, iterable] ', invece' bool '

Questo avviso appare anche nel codice originale di Justin Peel; e non so come risolverlo. Sono molto interessato a questo problema.