Mandelbrot - Optimierung

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Hallo,

Ich habe vor einigen Wochen die Aufgabe bekommen, ein Programm zu programmieren, welches dazu in der Lage ist, die Mandelbrot-Menge darzustellen. Da ich mit Python nicht eine extra Klasse für komplexe Zahlen oder mühsame Umrechnereien brauchte, entschied ich mich für diese Sprache. Nun bin ich soweit, dass die Mandelbrot-Menge bereits ziemlich sauber dargestellt wird, nur dauert die Rechnerei noch unglaublich lange und das Programm selber stockt auch ziemlich stark.

Code: Alles auswählen

#!/usr/bin/python
# Filename: mandelbrot.py

import matplotlib.pyplot as plt

def mandelbrot():
    global c, c0
    c = c ** 2 + c0

plt.xlabel(r'$re(z)$')
plt.ylabel(r'$i \cdot im(z)$')
plt.title('The Mandelbrot Set')

start = complex(-3,-3) # starting point

re = start.real # store real and imaginary part in c (part to calculate with) and c0 (fix during calculations)
im = start.imag
c = complex(re,im)
c0 = complex(re,im)

stepx = complex(0.01,0)
stepy = complex(0,0.01)
steps = 601
iterations = 100

for h in range(0,steps):
    for k in range(0,steps): # do this for all points in the selected area
        
        if abs(c) >= 2: # check divergence
            start += stepx
            
        else:
            
            for i in range(0,iterations):
                
                if abs(c ** 2 + c0) >= 2: # check divergence of value
                    break

                elif i < iterations - 1:
                    mandelbrot()
                
                elif i == iterations - 1:
                    plt.plot(round(c0.real,2), round(c0.imag,2), 'k,') # print(c0)
                    
                
            start += stepx
            
        re = start.real
        im = start.imag
        c = complex(re,im)
        c0 = complex(re,im)
    
    start = start + stepy - complex(steps * stepx.real,0)

plt.axis([-3,2,-2,2])
plt.grid(True)
plt.show()
Nun frage ich mich, wie ich das Script bzw. die Rechnereien optimieren kann, damit es nicht mehr so stockt und die Rechenzeit deutlich verringert wird? Ich habe vor Python noch nie irgendetwas programmiert und kenne mich demzufolge ziemlich schlecht mit Programmieren aus... Gibt es vllt. einen "besseren" Plotter für Python als matplot, der schneller arbeitet? Wie würdet ihr eine Zoomfunktion einbauen, die mir erlaubt, unendlich weit hineinzuzoomen? Dabei müsste ich halt einfach immer wieder die Rechnung für alle Punkte im ausgewählten Bereich durchführen, aber wie mache ich sowas?

Danke im Voraus für Hilfe...

MfG
Zuletzt geändert von Huy am Montag 26. April 2010, 16:17, insgesamt 1-mal geändert.
Gabelmensch
User
Beiträge: 79
Registriert: Montag 12. Oktober 2009, 11:50

Schneller Plotten geht auf jedenfall mit Gnuplot, vermutlich kann Gnuplot von sich selbst aus auch Mandelbrote darstellen.
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Vorweg: Ich habe matplotlib noch nie benutzt, habe es auch nicht installiert und kann deinen Code darum nicht testen.

Dennoch ein paar Gedanken:
Deine innere Schleife wird 36 Mio. mal durchlaufen. Falls du die Anzahl dieser Iterationen reduzieren kannst, solltet du das tun.
Auf jeden Fall solltest du bei der inneren Schleife ansetzen:
Der Funktionsaufruf ist - so wie die Funktion gemacht ist - nicht nur hässlich, sondern auch total überflüssig und zudem teurer, als wenn du die eine Zeile direkt in die Schleife setzt.
Die letzte Abfrage mit elif ist teuer und überflüssig. Es genügt doch, diesen Schritt durchzuführen, wenn die Iterationen durchgelaufen sind.
Und: Falls du mit Python 2.x arbeitest, solltest du xrange statt range verwenden, um nicht unnötigerweise immer wieder Listen zu erzeugen.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

@Gabelmensch: Gibt es irgendwo eine richtige Dokumentation von gnuplot.py? Habe es jetzt installiert, aber keine Ahnung, wie man es benutzt...

@numerix: Die Iterationen sollten schon um die 100 liegen, um saubere Resultate liefern zu können. Die Zeile setze ich also direkt in die Schleife. Bezüglich der letzten elif-Abfrage: Wenn ich diesen Schritt durchführe, wenn die Iterationen durchgelaufen sind, wird dann nicht auch im Falle "abs(c ** 2 + c0) >= 2" geplottet? Das möchte ich dann ja nicht... Oder stehe ich gerade aufm Schlauch?

Und ja, ich arbeite mit Python 2.6. Danke für die Tipps bisher.

MfG
Benutzeravatar
gkuhl
User
Beiträge: 600
Registriert: Dienstag 25. November 2008, 18:03
Wohnort: Hong Kong

Hier gibt's ein Beispiel mit Matplotlib: http://www.scipy.org/Tentative_NumPy_Tu ... et_Example
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Huy hat geschrieben:Bezüglich der letzten elif-Abfrage: Wenn ich diesen Schritt durchführe, wenn die Iterationen durchgelaufen sind, wird dann nicht auch im Falle "abs(c ** 2 + c0) >= 2" geplottet?
Nicht, wenn du die for-Schleife um einen else-Teil ergänzt. Der wird nur ausgeführt, wenn die Schleife nicht mittels break vorzeitig verlassen wurde.
Benutzeravatar
cofi
Python-Forum Veteran
Beiträge: 4432
Registriert: Sonntag 30. März 2008, 04:16
Wohnort: RGFybXN0YWR0

numerix hat geschrieben:Und: Falls du mit Python 2.x arbeitest, solltest du xrange statt range verwenden, um nicht unnötigerweise immer wieder Listen zu erzeugen.

Code: Alles auswählen

Python 2.6.5 (r265:79063, Mar 18 2010, 23:38:15) 
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> xrange(10**10)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to int
In Python 2.x muss man da leider vorsichtig sein, wenn man mit sehr grossen Zahlen arbeitet. Aber ein einfacher eigener Generator ist da auch kein Problem.
BlackJack

@cofi: Was soll die überflüssige Bemerkung denn jetzt? Der OP verwendet 2.x und `range()` da würde eine Länge von ``10**10`` auch nicht funktionieren -- alleine die Zeiger auf die Objekte würden da 37 Gigabyte benötigen. Das ist ein "Problem" was hier überhaupt nicht existiert.

@Huy: Pack die ganze Berechnung mal in eine Funktion. Lokale Namen werden schneller nachgeschlagen.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Was meinst du mit der ganzen Berechnung konkret? Alle drei for-Schleifen? Die innerste for-Schleife?

MfG
Gabelmensch
User
Beiträge: 79
Registriert: Montag 12. Oktober 2009, 11:50

Schau mal hier: http://rosettacode.org/wiki/Mandelbrot_set

gnuplot.py habe ich nur mal zum testen verwendet, wenn ich mal gnuplot in Py nutze, nehme ich meist subprocess.
BlackJack

@Huy: Ich meine "alles". Auf Modulebene sollten nach Möglichkeit sowieso nur ``import``\e und Definitionen von Funktionen, Klassen, und Konstanten stehen. Das macht Module einfacher wiederverwendbar und testbar.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Sorry, aber ich glaube, ich stehe gerade aufm Schlauch.

Code: Alles auswählen

#!/usr/bin/python
# Filename: mandelbrot.py

import matplotlib.pyplot as plt

plt.xlabel(r'$re(z)$')
plt.ylabel(r'$i \cdot im(z)$')
plt.title('The Mandelbrot Set')

start = complex(-3,-3) # starting point

re = start.real # store real and imaginary part in c (part to calculate with) and c0 (fix during calculations)
im = start.imag
c = complex(re,im)
c0 = complex(re,im)

stepx = complex(0.01,0)
stepy = complex(0,0.01)
steps = 601
iterations = 100

def mandelbrot(): 

    global start, re, im, c, c0, stepx, stepy, steps, iterations

    for h in xrange(0,steps):
        for k in xrange(0,steps): # do this for all points in the selected area
        
            if abs(c) >= 2: # check divergence
                start += stepx
            
            else:
            
                for i in xrange(0,iterations):

                    if abs(c ** 2 + c0) >= 2: # check divergence of value
                        break

                    else:
                        c = c ** 2 + c0
                
                else:
                    plt.plot(round(c0.real,2), round(c0.imag,2), 'k,') # print(c0)
                
                start += stepx
            
            re = start.real
            im = start.imag
            c = complex(re,im)
            c0 = complex(re,im)
    
        start = start + stepy - complex(steps * stepx.real,0)

mandelbrot()

plt.axis([-3,2,-2,2])
plt.grid(True)
plt.show()
Warum sollte das hier denn schneller berechnet werden so? Oder wie soll ich die Funktion gestalten? Sorry für die dummen Fragen, aber habe echt kein Plan von...

MfG
EyDu
User
Beiträge: 4881
Registriert: Donnerstag 20. Juli 2006, 23:06
Wohnort: Berlin

Weil der Zugriff auf lokale Namen schneller geht als auf globale.

Dein "global" in Zeile 24 solltest du übrigens loswerden. Funktionen sollten so wenige Seiteneffekte haben wie möglich. Du möchtest Parameter und Rückgabewerte verwenden.
Das Leben ist wie ein Tennisball.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Code: Alles auswählen

#!/usr/bin/python
# Filename: mandelbrot.py

import matplotlib.pyplot as plt

def mandelbrot(steps, start, stepx, stepy, iterations): 

    c = start
    c0 = start

    for h in xrange(0,steps):
        for k in xrange(0,steps): # do this for all points in the selected area
        
            if abs(c) >= 2: # check divergence
                start += stepx
            
            else:
            
                for i in xrange(0,iterations):

                    if abs(c ** 2 + c0) >= 2: # check divergence of value
                        break

                    else:
                        c = c ** 2 + c0
                
                else:
                    plt.plot(round(c0.real,2), round(c0.imag,2), 'k,') # print(c0)
                
                start += stepx
         
            c = start
            c0 = start
    
        start = start + stepy - complex(steps * stepx.real,0)

mandelbrot(601, complex(-3,-3), complex(0.01,0), complex(0,0.01), 100)

plt.xlabel(r'$re(z)$')
plt.ylabel(r'$i \cdot im(z)$')
plt.title('The Mandelbrot Set')

plt.axis([-3,2,-2,2])
plt.grid(True)
plt.show()
So? Oder habe ich das immer noch falsch verstanden?

MfG
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Huy hat geschrieben:So? Oder habe ich das immer noch falsch verstanden?
Hast du richtig verstanden. Die innere Schleife lässt sich allerdings noch weiter optimieren. Statt if .. else kannst du die Berechnung von c als ersten Schritt in der Schleife durchführen und danach nur mittels abs(c) das Abbruchkriterium prüfen. Wenn du dann noch multiplizierst statt zu potenzieren wird die Berechnung fast doppelt so schnell. Wenn du dann noch psyco hinzunimmst, ist die Berechnung gegenüber dem von dir gezeigten Code etwa um den Faktor 6 schneller. Das eigentliche Plotten ist da nicht enthalten, weil ich - wie gesagt - matplotlib nicht installiert habe.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Das Potenzieren und die if-else-Abfrage habe ich jetzt entsprechend verändert. Habe noch versucht psyco zu installieren, allerdings rausgefunden, dass das auf 64bit-System nicht läuft, bzw. es nicht unbedingt deutlich schneller ist als sonst, da das Script auf nem 64bit-System ohnehin ähnlich schnell läuft wie auf nem 32bit-System + psyco.

Gibt es noch irgendetwas zu optimieren? Wie kann man jetzt die Zoomfunktion verändern, bzw. wie verändere ich, was matplotlib.pyplot "macht"? Oder ist das im Endeffekt zu kompliziert?

MfG
Benutzeravatar
HerrHagen
User
Beiträge: 430
Registriert: Freitag 6. Juni 2008, 19:07

Matplotlib ist schnell genug. Das es so langsam läuft, hat IMO andere Ursachen.
Dein Hauptproblem ist das du matplotlib mit plt.plot jeden Punkt einzeln übergibst. So ist die Verwendung dieser Funktion einfach nicht gedacht. Die Punkte werden in vielfacher Ausführung in einer Datenstruktur gehalten die für diese Aufgabe nicht geeignet ist. Du erzeugst bei jedem Aufruf von plt.plot ein neues matplotlib-Objekt. Schau dir noch mal das folgende Beispiel an:
http://www.scipy.org/Tentative_NumPy_Tu ... et_Example
Hier wird das Ergebnis komplett in einer sehr effizienten Datenstruktur gehalten, nämlich einem numpy-array (entspricht c-Array + Zusatzinformationen) und dann als ganzes matplotlib übergeben. Entsprechend schnell liegt das Ergebnis vor.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Ich verstehe leider das Beispiel auf scipy nicht wirklich; also was genau gerechnet wird etc. Habe mal versucht, die Daten in einem numpy-Array festzuhalten, allerdings bin ich nicht ganz sicher, ob ich das richtig mache. Ich definiere zuerst ein Array mit empty() und füge dann immer Daten hinzu. Allerdings muss ich beim empty() Befehl die Dimensionen des Arrays festlegen, wobei ich ja keine Ahnung habe, wie lange es wird... Ist das schon die falsche Richtung? Wie muss ich das angehen?

MfG
BlackJack

@Huy: Natürlich weisst Du wie gross das Array werden muss -- Du weist doch wieviele Punkte Du in X- und Y-Richtung berechnest.
Huy
User
Beiträge: 18
Registriert: Freitag 23. April 2010, 20:25

Ehm, nein. Bzw. ich weiss wieviele Punkte ich berechne, aber ich weiss ja nicht, wieviele Punkte nicht divergieren. Ob ein Punkt divergiert oder nicht, finde ich ja erst mittels for-Schleife heraus. Oder?

MfG
Antworten