Mandelbrot Set ด้วย Python

Submitted by sugree on 6 April, 2007 - 00:41

Python in General

โจทย์ - http://www.codenone.com/node/242

ดูของ Ruby ใช้ OpenGL แล้วอยากเขียนมั่ง แต่ดันบ้าจี้ใช้ Tkinter ตามโจทย์ไปซะแล้ว เอาแค่นี้ก่อน โค้ดหลักของการหาค่าแล้วเอาไปวาดจะใช้ Tkinter

mp.py

#!/usr/bin/env python
 
from Tkinter import *
 
def mandelbrot(ul,lr,callback,step=(200,200),maxiter=255,limit=2):
    xs,ys = step
    x0,y0 = ul.real,ul.imag
    x1,y1 = lr.real,lr.imag
    xd,yd = (x1-x0)/xs,(y1-y0)/ys
    i = 0
    y = y0
    for i in range(ys):
        x = x0
        for j in range(xs):
            color = mandelbrot_point(x,y,maxiter,limit)
            callback(i,j,x,y,color)
        y += yd
 
class App:
    def __init__(self,parent):
        self.parse_args()
        self.init(parent)
 
    def parse_args(self):
        import sys
        args = sys.argv[1:]
        if not args:
            args = ['-2+1j','1-1j']
        self.upleft,self.lowright = map(complex,args[:2])
        self.width,self.height = 200,200
        self.maxiter = 30
        self.limit = 2
 
    def init(self,parent):
        self.frame = Frame(parent,width=self.width,height=self.height)
        self.frame.pack()
 
        self.canvas = Canvas(self.frame,width=self.width,height=self.height)
        self.canvas.pack(side=TOP)
 
    def palette(self,c):
        c = c*255.0/self.maxiter
        return '#%02x%02x%02x' % (c,(c+64)%256,(c+32)%256)
 
    def run(self):
        params = {
            'step': (self.width,self.height),
            'maxiter': self.maxiter,
            'limit': self.limit,
        }
        self.pixels = [0]*(self.width*self.height)
        mandelbrot(self.upleft,self.lowright,self.plot,**params)
 
    def plot(self,i,j,x,y,c):
        index = i*self.width+j
        #print i,j,c
        color = self.palette(c)
        p = self.canvas.create_line(j,i,j+1,i,fill=color)
        #print i,j,index
        self.pixels[index] = p
        self.frame.update()
 
import mandelbench
mandelbrot_point = mandelbench.is_mandel_4
 
if __name__ == '__main__':
    root = Tk()
    app = App(root)
    app.run()
    root.mainloop()

ส่วนฟังก์ชั่นจริงๆ จะใช้ mandelbrot_point ซึ่งข้างบนคือ mandelbench.is_mandel_4 หน้าตาของ mandelbench.py จะเป็นแบบนี้

mandelbench.py

def is_mandel_1(x,y,maxiter=30,limit=2):
    p = complex(x,y)
    i = 0
    z = 0+0j
    while abs(z) < limit and i < maxiter:
        z = z*z+p
        i += 1
    return i
 
def is_mandel_2(x,y,maxiter=30,limit=2):
    x0,y0 = x,y
    x2,y2 = x*x,y*y
 
    i = 0
    while x2+y2 < limit*limit and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i += 1
    return i
 
def is_mandel_3(x,y,maxiter=30,limit=2):
    x0,y0 = x,y
    x2,y2 = x*x,y*y
    lim = limit*limit
 
    i = 0
    while x2+y2 < lim and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i += 1
    return i
 
import mandelbrot
is_mandel_4 = mandelbrot.is_mandel_4

ซึ่ง is_mandel_4() ได้มาจาก Pyrex อีกที เพื่อการนี้ต้องใช้ไฟล์อีกสองไฟล์

mandelbrot.pyx

def is_mandel_4(double x,double y,int maxiter=30,int limit=2):
    cdef double x0,y0,x2,y2,lim
 
    x0,y0 = x,y
    x2,y2 = x*x,y*y
    lim = limit*limit
 
    i = 0
    while x2+y2 < lim and i < maxiter:
        x,y = x2-y2+x0,2*x*y+y0
        x2,y2 = x*x,y*y
        i = i+1
    return i

Makefile

PYINCLUDE = -I/usr/include/python$(shell python -c "import sys;print sys.version.split()[0]")
 
all:    mandelbrot.so
 
%.c:    %.pyx
        pyrexc $<
 
%.o:    %.c
        gcc -c -fPIC $(PYINCLUDE) $<
 
%.so:   %.o
        gcc -shared $< -lm -o $@
 
clean:
        @rm -f *.c *.o *.so *~ core core.*

ก่อนจะรันได้ต้องสั่ง make all ซักครั้งเพื่อให้ได้ mandelbrot.so เอาไว้ใช้ต่อไป

เฮ้อ inline ของ Ruby ดีกว่าเยอะ มีผลการรันนิดๆ ด้วย

Submitted by pphetra on 6 April, 2007 - 08:46.

Pyrex ของ python ดูเหมือนจะเอาไปใช้งานจริงๆจัง มันก็เลยหนักๆ
แต่ inline ของ ruby มันเป็น quick hack
เวลา run มันก็จะแอบไปสร้าง c code + wrapper code
แล้วก็สั่ง compile ให้เรา

Submitted by sugree on 7 April, 2007 - 12:31.

Tkinter มันช้าจริงๆ แค่ธรรมดายังแย่ zoom ยิ่งแย่กว่า ลองแบบ wxPython บ้าง

#!/usr/bin/env python
 
import wx
 
class MandelbrotSet:
    def __init__(self,x0,y0,x1,y1,w,h,limit=2,maxiter=30):
        self.x0,self.y0 = x0,y0
        self.rx,self.ry = (x1-x0)/w,(y0-y1)/h
        self.w,self.h = w,h
 
        self.limit = limit
        self.maxiter = maxiter
 
    def from_screen(self,x,y):
        return self.x0+self.rx*x,self.y0-self.ry*y
 
    def zoom_in(self,x,y):
        zx,zy = self.w/4,self.h/4
        return x-zx*self.rx,y+zy*self.ry,x+zx*self.rx,y-zy*self.ry,self.w,self.h
 
    def zoom_out(self,x,y):
        zx,zy = self.w,self.h
        return x-zx*self.rx,y+zy*self.ry,x+zx*self.rx,y-zy*self.ry,self.w,self.h
 
    def is_mandelbrot(self,x,y):
        p = complex(x,y)
        i = 0
        z = 0+0j
        while abs(z) < self.limit and i < self.maxiter:
            z = z*z+p
            i += 1
        return i
 
    def compute(self,callback):
        x = self.x0
        for xi in range(self.w):
            y = self.y0
            for yi in range(self.h):
                c = self.is_mandelbrot(x,y)
                callback(xi,yi,c)
                y -= self.ry
            x += self.rx
 
class wxMandelbrotSetViewer(wx.Frame):
    def __init__(self,x0,y0,x1,y1,w,h):
        super(wxMandelbrotSetViewer,self).__init__(None,-1,'Mandelbrot Set')
        self.dc = None
 
        self.SetSize((w,h))
        self.bitmap = wx.EmptyBitmap(w,h)
        self.panel = wx.Panel(self)
        self.panel.Bind(wx.EVT_PAINT,self.on_paint)
        self.panel.Bind(wx.EVT_LEFT_UP,self.on_zoom_in)
        self.panel.Bind(wx.EVT_RIGHT_UP,self.on_zoom_out)
 
        self.w,self.h = w,h
        self.m = MandelbrotSet(x0,y0,x1,y1,w,h)
 
    def on_zoom_in(self,event):
        x,y = self.m.from_screen(event.GetX(),event.GetY())
        self.m = MandelbrotSet(*self.m.zoom_in(x,y))
        self.dc = None
        self.Refresh()
 
    def on_zoom_out(self,event):
        x,y = self.m.from_screen(event.GetX(),event.GetY())
        self.m = MandelbrotSet(*self.m.zoom_out(x,y))
        self.dc = None
        self.Refresh()
 
    def on_paint(self,event):
        if not self.dc:
            self.dc = self.draw()
        dc = wx.PaintDC(self.panel)
        dc.Blit(0,0,self.w,self.h,self.dc,0,0)
 
    def palette(self,c):
        c = c*255.0/self.m.maxiter
        return map(int,[c,(c+64)%256,(c+32)%256])
 
    def draw(self):
        dc = wx.MemoryDC()
        dc.SelectObject(self.bitmap)
        def callback(x,y,c):
            r,g,b = self.palette(c)
            dc.SetPen(wx.Pen(wx.Colour(r,g,b),1))
            dc.DrawPoint(x,y)
        self.m.compute(callback)
        return dc
 
if __name__ == '__main__':
    app = wx.PySimpleApp()
    f = wxMandelbrotSetViewer(-2.0,1.0,1.0,-1.0,200,200)
    f.Show()
    app.MainLoop()

หน้าตาเหมือนใน Ruby เกือบเปีะ เพิ่ม zoom out เข้าไปอีกหน่อยกับให้เรียก callback เพื่อวาดทันที จะได้ประหยัดหน่วยความจำอีกนิดนึง

Submitted by bow_der_kleine on 15 April, 2007 - 09:55.

#!/usr/bin/env python
 
import numpy
import pylab
 
try :
    import is_mandeln
except :
    import f2py2e
    f2py2e.compile("""
    subroutine is_mandeln(x, y, c, limit, maxiter, dimenx, dimeny)
Cf2py intent(out) :: c
 
        integer dimenx, dimeny, maxiter
        double precision x(dimenx), y(dimeny)
        double precision c(dimeny, dimenx)
        double precision limit
        double complex p,z
        integer a, color
 
        do i=1, dimenx
            do j=1, dimeny
                p = dcmplx(x(i),y(j))
                a = 0
                z = dcmplx(0.0d0,0.0d0)
                do while ( ( abs(z) < limit) .and. ( a < maxiter) )
                    z = z*z + p
                    a = a+1
                enddo
                c(j,i) = (a*256/maxiter)
            enddo
        enddo
    end subroutine
    """,'is_mandeln')
    import is_mandeln
 
def getMandeln(x0, y0, x1, y1, w, h, limit=2, maxiter=30):
    x = numpy.arange(x0,x1,float((x1-x0)/w))
    y = yachs = numpy.arange(y0,y1,float((y1-y0)/h))
    return is_mandeln.is_mandeln(x,y,2,30)
 
def plotMandeln(mandeln):
    im = pylab.imshow(mandeln, interpolation='bilinear', origin='lower', extent=[-3,3,-3,3])
    pylab.show()
 
if __name__ == '__main__':
    plotMandeln(getMandeln(-2.0,1.0,1.0,-1.0,600,400))

คล้าย ๆ inline ของ Ruby ครับ แต่เป็น FORTRAN แล้วก็เอามา plot ด้วย matplotlib จริง ๆ ตอนพลอทนี่แหละครับที่เสียเวลา

ย้าย Codenone

ประกาศย้าย Codenone ไปใช้ Forum ของ Blognone แทนครับ ตามไปตั้งกระทู้ต่อได้ที่ Codenone Forum (รายละเอียดอ่านจากกระทู้ ย้าย Codenone ไปรวมกับ Blognone)

กระทู้เก่าๆ จะย้ายตามไปในภายหลัง ตอนนี้ปิดการโพสต์กระทู้ไว้ เหลือไว้เฉพาะอ้างอิงเท่านั้น

Codenone

Book navigation

Forum List