#! /usr/bin/env python #######################################################################################94x78## """ Mandelbrot program, for no particular reason. John Farrell 1999/11/02 CREDITS: Thanks to Robin Dunn for showing me how to use images and bitmaps. Thanks to Ionel Simionescu for showing me how to use the Numeric package, and consequently writing all of the computation code. Thanks to Alexander Smishlajev for suggestions on optimising loops and drawing. """ import Numeric WIDTH = 200 HEIGHT = 200 #coords = ((-1.0, -1.5), (0.0, -0.5)) coords = ((-2.0, -1.5), (1.0, 1.5)) WEIGHTS = (16, 1, 32) MAX_ITERATIONS = 255 zero = complex(0.0, 0.0) def setup(): # we use a single array for the real values corresponding to the x coordinates diff = coords[1][0] - coords[0][0] xs = 0j + (coords[0][0] + Numeric.arange(WIDTH).astype(Numeric.Float32) * diff / WIDTH ) # we use a single array for the imaginary values corresponding to the ycoordinates diff = coords[1][1] - coords[0][1] ys = 1j * (coords[0][1] + Numeric.arange(HEIGHT).astype(Numeric.Float32) * diff / HEIGHT ) # we build in direct correpondence with the pixels in the image c = Numeric.add.outer(ys, xs) z = Numeric.zeros((HEIGHT, WIDTH)).astype(Numeric.Complex32) c.shape = (HEIGHT * WIDTH,) z.shape = (HEIGHT * WIDTH,) return (c, z) def mandelbrot(progressMonitor, take=Numeric.take, repeat=Numeric.repeat, not_equal=Numeric.not_equal, greater_equal=Numeric.greater_equal, multiply=Numeric.multiply, add=Numeric.add,logical_not=Numeric.logical_not, equal=Numeric.equal, ): import time t = time.clock() (c, z) = setup() iterations = 0 i_no = Numeric.arange(WIDTH * HEIGHT) # non-overflow indices data = MAX_ITERATIONS + Numeric.zeros(WIDTH * HEIGHT) # initialize the "todo" arrays; # they will contain just the spots where we still need to iterate c_ = Numeric.array(c).astype(Numeric.Complex32) z_ = Numeric.array(z).astype(Numeric.Complex32) # replaces call to arange for each iteration bigarray = Numeric.arange(WIDTH * HEIGHT) while iterations <= MAX_ITERATIONS and len(i_no): progressMonitor.progress(iterations, MAX_ITERATIONS) # do the calculations in-place # On my box the mandelbrot calculation now takes 26 seconds # compared to 30 seconds which were needed without in-place notation. multiply(z_, z_, z_) add(z_, c_, z_) ### testing overflow = greater_equal(abs(z_), 2.0) ## not_overflow = logical_not(overflow) ## not_overflow = -(overflow - 1) # get the indices where overflow occured ### Inlined function compress... condition = not_equal(overflow, 0) oflowindicies = repeat(bigarray[:len(condition)], condition) ### Inlined function end overflowIndices = take( i_no, oflowindicies,-1) # set the pixel indices there for idx in overflowIndices: data[idx] = iterations # compute the new array of non-overflow indices ### Inlined function compress... condition = equal(overflow, 0) ## condition = not_equal(not_overflow, 0) noflowindicies = repeat(bigarray[:len(condition)], condition) ### Inlined function end i_no = take( i_no, noflowindicies, -1) #compress(not_overflow, i_no, bigarray) # update the todo arrays c_ = take( c_, noflowindicies, -1) #c_ = compress(not_overflow, c_, bigarray) z_ = take( z_, noflowindicies, -1) #z_ = compress(not_overflow, z_, bigarray) iterations = iterations + 1 progressMonitor.progress(MAX_ITERATIONS, MAX_ITERATIONS) print time.clock() - t return data try: from wxPython.wx import * class MandelTimer( wxTimer): '''Provides callback-insulated calculation''' def __init__( self, canvas, frame ): self.canvas = canvas self.frame = frame wxTimer.__init__( self ) def Notify( self, event=None ): data = mandelbrot(self.frame) import time t = time.clock() #-- John Farrell's version: Takes 13.8 seconds on my box, #-- that's about half as long as the mandelbrot calculation, #-- which takes 26 seconds on my box. ##pixels = [chr(0)] * len(data) * 3 ##for i in range(len(data)): ## di = data[i] ## i3 = i * 3 ## for j in range(3): ## pixels[i3 + j] = chr(self.colours[di, j]) ##img = wxEmptyImage(WIDTH, HEIGHT) ##d = string.join(pixels, '') ##img.SetData(d) ##self.bitmap = img.ConvertToBitmap() #-- Ionel Simionescu's version: takes 0.076 seconds on my box, #-- 180 times faster than John's version :-) data.shape = (HEIGHT, WIDTH) # build the pixel values pixels = Numeric.take( self.canvas.colours, data ) # create the image data bitmap = pixels.astype(Numeric.UnsignedInt8).tostring() # create the image itself image = wxEmptyImage(WIDTH, HEIGHT) image.SetData(bitmap) # store for canvas' use self.canvas.bitmap = image.ConvertToBitmap() print time.clock() - t # force display self.canvas.Refresh( ) class MandelCanvas(wxWindow): def __init__(self, parent, id = -1): wx.wxWindow.__init__(self, parent, id) self.parent = parent self.border = (1,1) self.SetSize(wxSize(WIDTH, HEIGHT)) self.SetBackgroundColour(wx.wxNamedColour("black")) self.bitmap = None self.colours = Numeric.zeros( (MAX_ITERATIONS + 1, 3) ) arangeMax = Numeric.arange(0, MAX_ITERATIONS + 1) self.colours[:,0] = Numeric.clip(arangeMax * WEIGHTS[0], 0, MAX_ITERATIONS) self.colours[:,1] = Numeric.clip(arangeMax * WEIGHTS[1], 0, MAX_ITERATIONS) self.colours[:,2] = Numeric.clip(arangeMax * WEIGHTS[2], 0, MAX_ITERATIONS) # start 3 millisecond timer... MandelTimer( self, parent ).Start( 10, oneShot=TRUE ) def OnPaint(self, event): dc = wxPaintDC(self) dc.BeginDrawing() if self.bitmap: dc.DrawBitmap(self.bitmap, 0, 0, false) dc.EndDrawing() class MyFrame(wxFrame): def __init__(self, parent, ID, title): wxFrame.__init__(self, parent, ID, title) self.CreateStatusBar() self.SetStatusText("") self.Centre(wxBOTH) mdb = MandelCanvas(self) self.SetAutoLayout(true) self.Layout() self.Fit() sz = self.GetClientSize() self.SetClientSize(wxSize(sz.width-7, sz.height-14)) def progress(self, done, outof): self.SetStatusText("%d/%d" % (done, outof)) class MyApp(wxApp): def OnInit(self): frame = MyFrame(NULL, -1, "Mandelbrot") frame.Show(true) self.SetTopWindow(frame) return true if __name__ == "__main__": app = MyApp(0) app.MainLoop() except ImportError: if __name__ == '__main__': print 'Starting non-gui run' class ProgMon: def progress( self, part, whole): print '.', import profile profile.run( "data = mandelbrot(ProgMon())")