# mandel_weave.py David MacQuigg ece175 3/17/08 '''Functions for computing the Mandlebrot set, and testing the speed of various implementations. Note [1] ''' from numpy import zeros # to initialize an array Note [2] import weave # to weave C code inline with Python Note [3] def getpts(cx0, cy0, dx): '''Compute color index values for an array of points in a Mandelbrot image. This version is simple Python with maximum clarity, and no concern about speed. See the next function getptsOpt() for a faster version. Given a starting point cx0, cy0 in the complex plane, and an increment dx for the space between points, return a color index value at each of 100 points. For each point, iterate the Mandelbrot generator z = z**2 + c, and save the number of iterations before the "orbit" of z reaches "escape altitude" |z| > 2. This number is an indication of how close c is to the "Mandelbrot Set", and will be used later to index a color from a palette. If the escape condition is not reached by iteration 999, save that number as an indicator that the iteration is most likely convergent (i.e. the point c is in the Mandelbrot Set). Doing 100 points at a time will minimize the overhead of calling this function. >>> getpts(-.222, 0.8, 1e-5) [279, 314, 266, 311, 273, 321, 302, 329, 496, 999,\ 482, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 635, 564, 604, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999,\ 999, 999, 999, 999, 999, 999, 999, 999, 999, 999] ''' vals = [] # an empty list for p in range(100): # 100 points cx = cx0 + p * dx # moving right x . . . . . . c = complex(cx, cy0) z = complex(0,0) for iternum in range(1000): # 999 max iteration if abs(z) > 2.0: break # point has diverged z = z**2 + c colorindex = iternum vals.append(colorindex) return vals def getptsOpt(cx0, cy0, dx): '''Optimized version of getpts() 1) Handle real and imag parts explicitly, no complex number conversions. z = zx + j*zy # a point in the complex plane c = cx0 + j*cy0 # starting point (z = 0) 2) Avoid square-root computation in abs(z). 3) Use x*x instead of x**2. >>> getptsOpt(-.222, 0.8, 1e-6) [279, 280, 283, 284, 292, 297, 342, 417, 388, 336,\ 314, 335, 286, 286, 341, 315, 320, 319, 271, 267,\ 266, 267, 268, 270, 273, 278, 296, 361, 347, 300,\ 311, 313, 294, 286, 283, 278, 276, 419, 288, 277,\ 273, 272, 272, 275, 436, 389, 282, 283, 285, 293,\ 321, 311, 317, 318, 852, 317, 314, 339, 320, 304,\ 302, 304, 308, 354, 331, 350, 347, 336, 330, 328,\ 329, 496, 334, 332, 337, 373, 363, 386, 570, 846,\ 496, 431, 539, 427, 433, 457, 484, 556, 734, 643,\ 999, 820, 576, 569, 503, 558, 715, 544, 509, 483] ''' vals = [] # an empty list for p in range(100): # 100 points cx = cx0 + p * dx # moving right x . . . . . . zx = zy = 0.0 for i in range(1000): # 999 max iteration zx2 = zx*zx zy2 = zy*zy if (zx2 + zy2) > 4.0 : break # escape zy = 2*zx*zy + cy0 zx = zx2 - zy2 + cx vals.append(i) return vals def getptsCstub(cx0, cy0, dx): '''stub for C-coded getpts function''' # Initialize a numpy array of 100 C-style ints vals = zeros(100, int) # Set up your C code as a raw multiline string mycode = r''' for (int p=0; p<100; p++) { // compute vals from cx0, cy0 and dx vals[p] = p; } ''' # List all Python variables your C code interacts with varlist = ['cx0', 'cy0', 'dx', 'vals'] # Run your C code inline weave.inline(mycode, varlist) return list(vals) # convert the numpy array back to a Python list def getptsC(cx0, cy0, dx): '''C-coded equivalent of getptsOpt function''' vals = zeros(100, int) # initialize an array of 100 C-style ints # Set up the inline C code as a raw multiline string mycode = r''' double cx = cx0 - dx; for (int p=0; p<100; p++) { double zx = 0.0; double zy = 0.0; int i; cx += dx; for (i=0; i<999; i++) { double zx2 = zx*zx; double zy2 = zy*zy; if ((zx2 + zy2) > 4.0) break; zy = 2*zx*zy + cy0; zx = zx2 - zy2 + cx; } vals[p] = i; } ''' weave.inline(mycode, ['cx0', 'cy0', 'dx', 'vals']) return list(vals) if __name__ == '__main__': '''Run tests if this is the main module.''' import time # for speed tests ## from mandelbrotC import getpts as getptsC # alternate C-coded version def speedtest(func=getptsOpt, cx0=0.3, cy0=0.5, dx=1e-3): '''Test specified function and return speed in points per second. >>> speedtest(getptsOpt, -.222, 0.8, 0.0001) 730 ''' t0 = t = time.clock() tmax = t0 + 1.0 # 1 second total time Npoints = 0 while(t < tmax): func(cx0, cy0, dx) Npoints += 100 t = time.clock() return Npoints/(t - t0) def compare(func1, func2): '''Compare the speed of two functions on a variety of problems. >>> compare(getpts, getptsOpt) Compare 1) getpts to 2) getptsOpt Region speed1 speed2 ratio convergent 547 649 1.19 edge 3096 3976 1.28 divergent 24200 29064 1.20 ''' print "Compare 1) %s to 2) %s" % (func1.__name__, func2.__name__) print "Region speed1 speed2 ratio" spd1 = speedtest(func1, cx0=0.1) spd2 = speedtest(func2, cx0=0.1) print " convergent %10d %10d %8.2f" % (spd1, spd2, spd2/spd1) spd1 = speedtest(func1, cx0=0.3) spd2 = speedtest(func2, cx0=0.3) print " edge %10d %10d %8.2f" % (spd1, spd2, spd2/spd1) spd1 = speedtest(func1, cx0=0.5) spd2 = speedtest(func2, cx0=0.5) print " divergent %10d %10d %8.2f" % (spd1, spd2, spd2/spd1) from doctest import testmod testmod(verbose=True) ''' Notes [1] References on Mandelbrot set and fractals. http://en.wikipedia.org/wiki/Mandelbrot_set - good overview, then lots of mathematical detail - image gallery of a zoom sequence (d://downloads/mandelbrots/) http://4dsolutions.net/ocn/fractals.html - nice images, good short explanation http://www.math.utah.edu/~pa/math/mandelbrot/mandelbrot.html - Java applet including control panel with color editor, etc, [2] http://numpy.scipy.org/ - package for scientific computing with Python [3] http://scipy.org/Weave - package for acceleration of Python by inclusion of C. '''