[Numpy-discussion] How to remove fortran-like loops with numpy?

[David Goldsmith]

I look forward to an instructive reply: the "Pythonic" way to do it would be to take advantage of the facts that Numpy is "pre-vectorized" and uses broadcasting, but so far I haven't been able to figure out (though I haven't yet really buckled down and tried real hard) how to broadcast a conditionally-terminated iteration where the number of iterations will vary among the array elements.  Hopefully someone else already has. :-)

DG

--- On Mon, 6/8/09, Juanjo Gomez Navarro <juanjo.gomeznavarro at gmail.com> wrote:

> From: Juanjo Gomez Navarro <juanjo.gomeznavarro at gmail.com>
> Subject: [Numpy-discussion] How to remove fortran-like loops with numpy?
> To: numpy-discussion at scipy.org
> Date: Monday, June 8, 2009, 2:52 PM
> 
> Hi all, 
> 
> I'm new in numpy. Actually, I'm new in Python. In
> order to learn a bit, I want to create a program to plot the
> Mandelbrot set. This program is quite simple, and I have
> already programmed it. The problem is that I come from
> fortran, so I use to think in "for" loops. I know
> that it is not the best way to use Python and in fact the
> performance of the program is more than poor. 
> 
> 
> 
> Here is the program:
> 
> #!/usr/bin/python
> 
> 
> 
> import numpy as np
> 
> import matplotlib.pyplot as plt
> 
> 
> 
> # Some parameters
> Xmin=-1.5
> 
> Xmax=0.5
> 
> Ymin=-1
> 
> Ymax=1
> 
> 
> Ds = 0.01
> 
> 
> # Initialization of varibles
> X = np.arange(Xmin,Xmax,Ds)
> 
> 
> Y = np.arange(Ymax,Ymin,-Ds)
> 
> 
> 
> N = np.zeros((X.shape[0],Y.shape[0]),'f')
> 
> 
> ############## Here are inefficient the calculations
> ################
> for i in range(X.shape[0]):
>   for j in range(Y.shape[0]):
>     z= complex(0.0, 0.0)
>     c = complex(X[i], Y[j])
> 
>     while N[i, j] < 30 and abs(z) < 2:
> 
>       N[i, j] += 1
>       z = z**2 + c
>     if N[i, j] == 29:
>       N[i, j]=0
> ####################################################
> 
> # And now, just for ploting...
> N = N.transpose()
> 
> fig = plt.figure()
> 
> plt.imshow(N,cmap=plt.cm.Blues)
> 
> 
> plt.title('Mandelbrot set')
> 
> plt.xticks([]); plt.yticks([])
> 
> 
> plt.show()
> 
> fig.savefig('test.png')
> 
> 
> 
> As you can see, it is very simple, but it takes several
> seconds running just to create a 200x200 plot. Fortran takes
> same time to create a 2000x2000 plot, around 100 times
> faster... So the question is, do you know how to programme
> this in a Python-like fashion in order to improve seriously
> the performance?
> 
> 
> 
> Thanks in advance
> -- 
> Juan José Gómez Navarro
> 
> Edificio CIOyN, Campus de Espinardo, 30100
> Departamento de Física 
> Universidad de Murcia
> Tfno. (+34) 968 398552
> 
> Email: juanjo.gomeznavarro at gmail.com
> 
> Web: http://ciclon.inf.um.es/Inicio.html
> 
> 
> -----Inline Attachment Follows-----
> 
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>