ucgamdo@ucl.ac.uk
2003-Sep-17 11:45 UTC
[R] Fractals in R and having fun! (and more persp and color)
Well, I started playing with fractals in R, and wrote a function to
generate de Mandelbrot set, which might
be of interest to some people
###########################################################################
# Mandelbrot set
###########################################################################
mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates
y = c(-1.8, 1.8), # y coordinates
b = .05, # by 'steps'
iter = 20) # maximum number of iterations
{
m = NULL # will store the results, this is the 'image' matrix
for(i in seq(x[1], x[2], by = b)) {
r = NULL # stores part of the iteration results
for(j in seq(y[1], y[2], by = b)) {
it = iter # will hold iteration at which point (i, j) breaks off
c = c(i, j) # initial point
z = c(i, j) # i: real part; j: imaginary part
for(k in 1:iter) {
# the Mandelbrot iteration formulae: z -> z*z + c
z = c(z[1]^2 - z[2]^2, 2 * z[1]*z[2]) + c
# tests if point breaks off
if((z[1] + z[2])^2 > 4) { it = k; break }
}
r = c(r, it) # stores iteration results
}
# constructs the 'image' matrix
m = rbind(m, r)
}
# the output fractal object
fractal = list(x = seq(x[1], x[2], by = b), # x coordinates
y = seq(y[1], y[2], by = b), # y coordinates
z = m) # it matrix
}
######################################################################
# here goes how it works
######################################################################
frac <- mandelbrot()
image(frac) # perhaps not very nice!
# more resolution, beware, this might take some time to calculate
frac <- mandelbrot(b = .01)
image(frac)
# now here comes the fun, lets do a persp plot of the set!
persp(log(frac$z), border = NA, theta = 225, phi = 45, col = "gray",
shade
= .4)
###END###
Well, here comes what would have been my question. How to put some nice
colors to the above perspective plot? I wanted to post this question some
weeks ago, but it was answered yesterday after all the discussion on persp
and colors. If you see my previous post, you'll the definition of
surf.colors used below:
persp(log(frac$z), border = NA, theta = 225, phi = 45, col surf.colors(frac$z,
col = c(gray(seq(0.5, 1, len = 39)), "black")), shade .4)
that was what I always wanted to do, and this was the source of my
confusion when checking the code in persp demo.
But now some new question arises, which I'll leave to the interested hacker:
Any C programmer who volunteers to write and link c code to R to calculate de
Mandelbrot set faster?
Does anyone have any idea on how to 'smooth' the valleys and ridges that
rise towards the set?, i.e. to avoid the stairs like appearence of the plot?
I hope you can have some fun with this!
Martin Maechler
2003-Sep-18 12:48 UTC
[R] Fractals in R and having fun! (and more persp and color)
>>>>> "Mario" == Mario ??? <ucgamdo at ucl.ac.uk> >>>>> on Wed, 17 Sep 2003 12:45:51 +0100 writes:Mario> Well, I started playing with fractals in R, and wrote Mario> a function to generate de Mandelbrot set, which might Mario> be of interest to some people Mario> ################################################################### Mario> # Mandelbrot set Mario> ################################################################### Mario> mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates Mario> y = c(-1.8, 1.8), # y coordinates Mario> b = .05, # by 'steps' Mario> iter = 20) # maximum number of iterations Mario> <..................> Mario> <..................> Mario> <..................> Well, only a bit more than year ago I had posted my version of mandelbrot() and a drawing function, see http://finzi.psych.upenn.edu/R/Rhelp02a/archive/5898.html [ I have a slightly updated version of the R code, that is now available as ftp://stat.ethz.ch/U/maechler/R/Mandelbrot.R ] which is an order of magnitude more efficient than yours (22 x faster for your "b = 0.01", *1), *2) and with an iterImage() function for drawing its result in several (simple) color schemes. I agree that the proper solution would *definitely* use C code! Your idea of using persp() -- while seen frequently in fractal books, is new ``for R'' however, and nice! Thank you. Regards, Martin *1) mainly because I only have the iteration loop and do the rest vectorized, but also because I have a check for symmetry and make use of it when appropriate *2) For speed comparison, note that I use 10 x more iterations and a much higher resolution by default *3) When I do the comparison, I see that your function's result slightly differs from mine -- not for the Mandelbrot set itself, but for the very ``outer points'' (the dark orange ones in the image plot). Your result is asymmetric and hence can't be quite correct.
ucgamdo@ucl.ac.uk
2003-Sep-18 16:28 UTC
[R] Fractals in R and having fun! (and more persp and color)
Hi Martin, Thanks a lot for pointing your function out to me. I tried it, and indeed it is very fast. I will have a close look at the asymmetry problem in my own code, although this sounds pretty bizarre to me, I mean, I realised that there was some asymmetry but I didn't really payed attention to that. I think that the high number of iterations that can be achieved with your function will be the solution to obtaining 'smoother' slopes in the perspective plot (i.e. the outer regions that rise progressively towards the set). If you come with any ideas on how to improve that perspective let me know. I think I'll have a go to the c code. Thanks, Mario. At 10:21 18/09/03 +0200, you wrote:> > Mario> Well, I started playing with fractals in R, and wrote > Mario> a function to generate de Mandelbrot set, which might > Mario> be of interest to some people > > Mario>###################################################################> Mario> # Mandelbrot set > Mario>###################################################################> > Mario> mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates > Mario> y = c(-1.8, 1.8), # y coordinates > Mario> b = .05, # by 'steps' > Mario> iter = 20) # maximum number ofiterations> > <..................> > <..................> > <..................> > >Well, only a bit more than year ago I had posted my version of >mandelbrot() and a drawing function, see > > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/5898.html > > [ I have a slightly updated version of the R code, that is now > available as ftp://stat.ethz.ch/U/maechler/R/Mandelbrot.R > ] > >which is an order of magnitude more efficient than yours (22 x >faster for your "b = 0.01", *1), *2) >and with an iterImage() function for drawing its result in >several (simple) color schemes. >I agree that the proper solution would *definitely* use C code! > >Your idea of using persp() -- while seen frequently in fractal books, >is new ``for R'' however, and nice! Thank you. > > >Regards, >Martin > >*1) mainly because I only have the iteration loop > and do the rest vectorized, but also because I have a check > for symmetry and make use of it when appropriate > >*2) For speed comparison, note that I use 10 x more iterations > and a much higher resolution by default > >*3) When I do the comparison, I see that your function's result > slightly differs from mine -- not for the Mandelbrot set > itself, but for the very ``outer points'' (the dark orange > ones in the image plot). Your result is > asymmetric and hence can't be quite correct. > >
I decided to take on the 'proper' solution to calculate the
Mandelbrot set in R, i.e. to do the raw calculations in C and then link
that code with R. I thought it would be a hard task, but I was pleasantly
surprised when I saw how easily was to write the bit of C code (I am not a
C programmer myself!) and how even easier was to link it with R using .C()
and R CMD SHLIB. If you are interested (or anyone else) here is the (very
fast) code to calculate the mandelbrot set, I hope you enjoy playing with
it as much as I did.
# BEGINNING OF R CODE: ################################################
#######################################################################
# Function to calculate the Mandelbrot set. This function calls a #
# C routine in order to perform the calculations faster. #
# #
# Written by Mario dos Reis. September 2003 #
#######################################################################
mandelbrot <- function(x = c(-3, 1), # x limits
y = c(-1.8, 1.8), # y limits
nx = 600, # x resolution
ny = 600, # y resolution
iter = 20) # maximun number of iterations
{
xcoo <- seq(x[1], x[2], len = nx) # x coordinates
ycoo <- seq(y[1], y[2], len = ny) # y coordinates
set = numeric(nx*ny) # this will store the output of
# the C routine
# This is the call to the C function itself
the.set = .C("mandelbrot",
xcoo = as.double(xcoo),
ycoo = as.double(ycoo),
nx = as.integer(nx),
ny = as.integer(ny),
set = as.integer(set),
iter = as.integer(iter))$set
# Create a list with elements x, y and z,
# suitable for image(), persp(), etc. and return it.
return(list(x = xcoo, y = ycoo, z = matrix(the.set, ncol = ny, byrow = T)));
}
# END OF R CODE. #######################################################
/* BEGINNING OF C CODE: ***********************************/
#include <stdio.h>
/***********************************************************
* Function to evaluate a series of complex numbers that *
* form a subset of the Argand plane, to see if *
* they belong to the Mandelbrot set. *
* This function has been written with the intention of *
* linking it to some R code in order to create a *
* practical way to play with the set. *
* *
* Written by Mario dos Reis. September 2003 *
* *
***********************************************************/
void mandelbrot(double *xcoo, double *ycoo, int *nx,
int *ny, int *set, int *iter)
/* 'xcoo' and 'ycoo' are the x and y coordinates respectively
* of all the points to be evaluated.
* 'nx' and 'ny' number of divisions along the x and y axes
* 'set' will store the practical output of the function
* 'iter' is the maximun number of iterations
*/
{
int i, j, k;
double z[2], c[2], oldz[2];
for(i = 0; i < *nx; i++) {
for(j = 0; j < *ny; j++) {
/* initialise the complex point to be tested
* c[0] (or z[0]) is the real part
* c[1] (or z[1]) is the imaginary part
*/
c[0] = xcoo[i]; c[1] = ycoo[j];
z[0] = 0; z[1] = 0;
for(k = 1; k < *iter + 1; k++) {
oldz[0] = z[0]; oldz[1] = z[1];
// the mandelbrot mapping z -> z^2 + c
z[0] = oldz[0]*oldz[0] - oldz[1]*oldz[1] + c[0];
z[1] = 2 * oldz[0]*oldz[1] + c[1];
if((z[0]*z[0] + z[1]*z[1]) > 4) {
break;
}
}
/* fills the set vector
* notice the trick 'i * (*ny) + j' to find
* the apporpiate position of the output in the
* vector set. The R function will take care of
* transforming this vector into a suitable matrix
* for plotting, etc.
*/
set[i * (*ny) + j] = k;
}
}
return;
}
/* END OF C CODE. ******************************************
To anyone interested in trying out these functions, just save the C code
into a file, say 'mandelbrot.c', to be able to use the C code with R a
shared object needs to be created. This is done simple running the
following command from the directory where 'mandelbrot.c' is stored
(note
this code was tried in Linux, for Win users I think the instructions are
similar but you need you check it out):
R CMD SHLIB mandelbrot.c
This will create the file 'mandelbrot.so' in that same directory.
Start R from this directory and paste the above R code into the console. Type
> dyn.load("mandelbrot.so")
in order to load the shared object. Now you are ready to appreaciate in
full detail the set.
> image(mandelbrot(), col = c(heat.colors(49), "black"))
This function is many orders of magnitude faster than the one I had
published previously
(see http://maths.newcastle.edu.au/~rking/R/help/03b/3390.html)
that function also had a slight error, it tested (z[1] + z[2])^2 > 4
instead of z[1]^2 + z[2]^2 > 4 which creates the asymmetric result you
pointed out in the outer regions of the set. The C code is also at
extremely fast with high zooms and iteratons, try:
xcoo = c(-1.1854735004165451, -1.1854647240344973)
ycoo = c(-0.3057632375298113, -0.3057544611477636)
frac <- mandelbrot(x = xcoo, y = ycoo, iter = 2000) # took ~ 12s in P4 2.2GHz
image(frac, col = c(heat.colors(49), "black"))
This is at least 19x faster than the code you suggested me. (same settings
took ~ 230s in the same machine).
Regards,
Mario.
At 10:21 18/09/03 +0200, you wrote:>
> Mario> Well, I started playing with fractals in R, and wrote
> Mario> a function to generate de Mandelbrot set, which might
> Mario> be of interest to some people
>
> Mario>
###################################################################> Mario> # Mandelbrot set
> Mario>
###################################################################>
> Mario> mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates
> Mario> y = c(-1.8, 1.8), # y coordinates
> Mario> b = .05, # by
'steps'
> Mario> iter = 20) # maximum number of
iterations>
> <..................>
> <..................>
> <..................>
>
>Well, only a bit more than year ago I had posted my version of
>mandelbrot() and a drawing function, see
>
> http://finzi.psych.upenn.edu/R/Rhelp02a/archive/5898.html
>
> [ I have a slightly updated version of the R code, that is now
> available as ftp://stat.ethz.ch/U/maechler/R/Mandelbrot.R
> ]
>
>which is an order of magnitude more efficient than yours (22 x
>faster for your "b = 0.01", *1), *2)
>and with an iterImage() function for drawing its result in
>several (simple) color schemes.
>I agree that the proper solution would *definitely* use C code!
>
>Your idea of using persp() -- while seen frequently in fractal books,
>is new ``for R'' however, and nice! Thank you.
>
>
>Regards,
>Martin
>
>*1) mainly because I only have the iteration loop
> and do the rest vectorized, but also because I have a check
> for symmetry and make use of it when appropriate
>
>*2) For speed comparison, note that I use 10 x more iterations
> and a much higher resolution by default
>
>*3) When I do the comparison, I see that your function's result
> slightly differs from mine -- not for the Mandelbrot set
> itself, but for the very ``outer points'' (the dark orange
> ones in the image plot). Your result is
> asymmetric and hence can't be quite correct.
>
>