NumPy vs. iterative vs. ???

[Tim Lavoie]

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()