NumPy vs. iterative vs. ???

[Tim Hochberg]

tim.lavoie at mts.net (Tim Lavoie) writes:

Fun!

I played around with your code and came up with the attached. It's
about twice as fast by pruning the set of points to look at each
iteration, but no longer returns the full arrays of c and z. (I think
that could be added back in if nescessary). I agree that NumPy may not
be the most elegant way to attack this kind of problem, but it's
probably as good as your going to get from within Python.

Enjoy.


-tim

#!/usr/bin/env python
#
# Mandelbrot ASCII-art using Numeric Python 1.0beta1
#
# Rob Hooft, 1996. Distribute freely.
#
# Additional kluges and hackery, Tim Lavoie, 2000
# Yet more by Tim Hochberg, 2000

from Numeric import *
from NumTut import *
import profile


def fraccalc(frac='z**2+c',LowX=-2.1, HighX=0.7, LowY=-1.2, HighY=1.2, \
                                          stepx=250, stepy=250, maxiter=30):
	xx=arange(LowX,HighX,(HighX-LowX)/stepx)
	yy=arange(HighY,LowY,(LowY-HighY)/stepy)*1j

	c=ravel(xx+yy[:,NewAxis])
	z=zeros(c.shape,Complex)

	fraccode = compile(frac,'<string>','eval')
	indices = arange(len(z))

	iter_count = zeros([stepx*stepy])
	for iter in range(maxiter):
		z = eval(fraccode)
		diverged = greater(abs(z),2.0)
		divergedIndices = compress(diverged, indices)
		for index in divergedIndices:
			iter_count[index] = (iter % 2)*255
		# This is the equivalent of three compresses only faster
		nonzeroNotDiverged = nonzero(logical_not(diverged))
		z = take(z, nonzeroNotDiverged)
		c = take(c, nonzeroNotDiverged)
		indices = take(indices, nonzeroNotDiverged)
		
	iter_count.shape = (stepy,stepx)
	return iter_count
		

if __name__ == "__main__":
	profile.run('results = fraccalc()')

	view(results)
	raw_input()



> I've been tinkering (no, not TkInter-ing, at least yet) with fractals in
> Python, but am not sure if my approach is sensible. By the way, I mean
> iterative-by-nature, Mandelbrot set sort of fractals.
> 
> I like the idea of NumPy, in that it does a lot for me and makes for more
> elegant code. However, I'm not sure how well it fits this application, since
> it prefers to apply each operation to every element in the array. Also,
> keeping several arrays (some complex-type) for the whole works is very
> RAM-hungry as well, at an element per pixel per array.
> 
> I do like the flexibility of doing things like this in Python, for instance
> I pass the function to iterate as a parameter, and eval() it.
> 
> Doing things the old, nested-loop and iteration way, I can more easily
> evaluate the individual elements, and use a *TINY* fraction of the memory.
> The problem is, I get whacked with function-call overhead, so it all goes
> from slow to worse. I do like the idea of passing back the raw result data
> though, so that I can apply different color parameters without re-doing all
> the basic calculation.
> 
> I've considered breaking up the data into chunks that I can process in
> smaller pieces, say a row at a time; would this make more sense than either
> previous approach?
> 
> The code which follows started out as Rob Hooft's example from the NumPy
> demos, further mangled by myself. Ideas, suggestions and outright derision
> are all ok...
> 
>     Thanks,
>     Tim
> 
> -- 
> "Now let us peel back the foreskin of misconception and apply
> the wire brush of enlightenment" -- Geoff Miller
> 
> 
> 
> #!/usr/bin/env python
> #
> # Mandelbrot ASCII-art using Numeric Python 1.0beta1
> #
> # Rob Hooft, 1996. Distribute freely.
> #
> # Additional kluges and hackery, Tim Lavoie, 2000
> 
> from Numeric import *
> from NumTut import *
> import profile
> 
> 
> def fraccalc(frac='z**2+c',LowX=-2.1, HighX=0.7, LowY=-1.2, HighY=1.2, \
>                                           stepx=250, stepy=250, maxiter=30):
> 	xx=arange(LowX,HighX,(HighX-LowX)/stepx)
> 	yy=arange(HighY,LowY,(LowY-HighY)/stepy)*1j
> 
> 	c=ravel(xx+yy[:,NewAxis])
> 	z=zeros(c.shape,Complex)
> 	
> 	iter_count=resize(array([128]),c.shape)
> 	fin = zeros(c.shape)
> 	print iter_count.shape, z.shape, c.shape
> 
> 	
> 	fraccode = compile(frac,'<string>','eval')
> 
> 	for iter in range(maxiter):
> 		z = where(fin,z,eval(fraccode))
> 		fin = where(greater(abs(z),2.0),1,0)
> 			# c=where(finished,0+0j,c)
> 			# z=where(finished,0+0j,z)
> 		iter_count=where(equal(fin,1), iter_count, iter)
> 
> 	# return iter_count.tostring()
> 	iter_count = where(equal(fmod (iter_count, 2), 1),0,255)
> 	iter_count.shape = (stepy,stepx)
> 	z.shape = (stepy,stepx)
> 	c.shape = (stepy,stepx)
> 	print iter_count.shape, z.shape, c.shape
> 	return (iter_count, z, c)
> 		
> 
> if __name__ == "__main__":
> 	profile.run('results = fraccalc()')
> 
> 	view(results[0])
> 	raw_input()
> 
>