I now see that I have had to write a reason why I want a plain sage
script, which generates Mandelbrot set. I want to show for students,
who haven't coding before, how easy is to generate Mandelbrot set
with sage. After adding some additional functions and classes to
sage init everything should be okay.
Thank you, for help.
> On Sep 1, 9:52 pm, Robert Bradshaw <[EMAIL PROTECTED]> wrote:
> On Sep 1, 2008, at 3:18 AM, Simon King wrote:
>
> > Hi,
>
> > one hint for doing *fast* computations:
> > Use Cython, rather than pure Python!
>
> Certainly!
>
>
>
> > From my point of view, Cython is one of the main strengths of Sage:
> > You can do your computations with C data types, and you should do so,
> > because it is faster than using "complicated" data types such as Sage-
> > Integer. But you still have the nice feel of Python.
>
> > A few months ago, I made up a Mandelbrot creator myself, just for fun
> > and for practicing. It consists of a class "Mandelbrot". In the
> > initiation, you can pass various parameters (the x- and y-range, the
> > number of pixel on both axes, the maximal number of iterations).
>
> > If i understand the meaning of your parameters correctly, the
> > following produces the image that you want:
> > sage: attach Mandelbrot.pyx
> > # the code is below
> > sage: MB=Mandelbrot(xmin=-2,xmax=1,xsteps=200,
> > ymin=-1,ymax=1,ysteps=200, maxIter=200)
> > # i think this is your example
> > sage: time MB.make()
> > CPU times: user 0.34 s, sys: 0.00 s, total: 0.34 s
> > Wall time: 0.34 s
> > sage: show(MB.img)
> > # Objects of class <Mandelbrot> have an attribute img, which is
> > computed with the method make()
> > # and shown with show()
>
> > Main reasons why the code is fast:
> > 1. The computations are done with data of type <double>. For example,
> > cdef double cr
> > cdef double ci
> > describes real and imaginary part of the complex number for which we
> > iterate
> > 2. All loop variables are of type <int>. You do
> > for i in range(n):
> > ...
> > which is slow. Faster is
> > for i in xrange(n):
> > ...
> > Using C types is still faster:
> > cdef int i
> > cdef int N=n # hence, the <Integer> n is converted into a C
> > <int>
> > for i from 0 <= i < N:
> > ...
>
> Note that if i is a cdef variable, one only can do
>
> for i in range(N)
>
> and it will be converted to the for-from style loop.
>
>
>
>
>
> > One very useful feature of Sage for developing fast code is
> > "annotation". By typing
> >> sage -cython -a Mandelbrot.pyx
> > from the command line, i get a file Mandelbrot.html, which i can view
> > with a browser. That file shows me my own code -- with hints on how
> > efficient the code is (the more yellow a line of code the more slow it
> > is executed).
>
> > The code below also relies on an experimental package for image
> > manipulation, that you can get (from the command line) with
> >> sage -i PIL-1.1.5
>
> > Cheers
> > Simon
>
> > ###################################################
> > ## Here is the code of Mandelbrot.pyx
> > import sage
> > import sage.all
> > import Image
> > import colorsys
> > # requires PIL
>
> > class Mandelbrot:
> > def __init__(self,xmin=-2.4,xmax=1.3,xsteps=600,
> > ymin=-1.7,ymax=1.7,ysteps=600, maxIter=120):
> > """
> > Class for creating a picture of the Mandelbrot set.
> > MB = Mandelbrot()
> > The picture can be created with MB.make() and shown with
> > show(MB.img)
> > """
> > self.colors = {} # used for cacheing
> > self.img = Image.new("RGB",(xsteps,ysteps),(0,0,0))
> > self.d2 = float(4) # sequence is unbounded if it comes outside
> > a circle of radius 2
> > self.N = int(maxIter) # maximal number of iterations
> > if xmin>=xmax:
> > raise ValueError, "xmax<=xmin"
> > self.xmin=float(xmin)
> > self.xmax=float(xmax)
> > if xsteps==0:
> > raise ValueError, "xsteps==0"
> > self.dx=(self.xmax-self.xmin)/float(xsteps)
> > self.xsteps = int(xsteps)
> > if ymin>=ymax:
> > raise ValueError, "ymax<=ymin"
> > self.ymin=float(ymin)
> > self.ymax=float(ymax)
> > if ysteps==0:
> > raise ValueError, "ysteps==0"
> > self.dy=(self.ymax-self.ymin)/float(ysteps)
> > self.ysteps = int(ysteps)
>
> > def color(self,i):
> > "computes the color of a pixel if the iteration leaves the
> > circle of radius 2 after i steps"
> > # the result is cached in self.colors
> > if i==self.N:
> > return (0,0,0)
> > CL = self.colors.get(i)
> > if CL:
> > return CL
> > (r,g,b) = colorsys.hsv_to_rgb(1.0-float(i)/self.N,1.0,1.0)
> > CL = (int(256*r)-1,int(256*g)-1,int(256*b)-1)
> > self.colors[i] = CL
> > return CL
>
> > def point(self,x,y):
> > """
> > return the number of iterations after which point (x,y)
> > leaves the circle of radius 2.
> > Maximum number of iterations is self.N
> > """
> > # Morally, this is C code!
> > cdef double cr = self.xmin+x*self.dx
> > cdef double ci = self.ymin+y*self.dy
> > cdef double zr = 0
> > cdef double zi = 0
> > cdef double zrneu,zineu
> > cdef double zr2=0
> > cdef double zi2=0
> > cdef double d2 = self.d2
> > cdef int N = self.N
> > cdef int i
> > for i from 0 <= i <= N:
> > zrneu = zr2-zi2 + cr
> > zineu = 2*zr*zi + ci
> > zr=zrneu
> > zi=zineu
> > zr2=zr**2
> > zi2=zi**2
> > if (zr2+zi2)>d2:
> > self.img.putpixel((x,y),self.color(i))
> > return i
> > return self.N
>
> > def make (self):
> > "Create a picture of the Mandelbrot set"
> > cdef int i
> > cdef int j
> > cdef int xm = self.xsteps
> > cdef int ym = self.ysteps
> > for i from 0 <= i < xm:
> > for j from 0 <= j < ym:
> > self.point(i,j)
>
> Another optimization, in you're very inner loop it would be faster to
> use zr*zr rather than zr**2 (as the latter uses pow() which is much
> slower than a multiplication). Still, I bet the bulk of your time is
> spent on the line "self.img.putpixel((x,y),self.color(i))" rather
> than computing the set.
>
> Another point (though this is just style, but I would like to show
> the flexibility) is that one can declare multiple variables on the
> same line (e.g. cdef "double zr, zi, zr2, zi2") and do parallel
> assignments just as in Python "zr, zi = zr2-zi2 + cr, 2*zr*zi + ci"
>
> - Robert
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
