# Speeding up this fractal-generating code

I used the code below (which is a sample from this gist containing more similar code) in my answer to my own question about Mandelbrot-like sets for functions other than the simple quadratic on Math.SE to generate this image:

cosineEscapeTime =
Compile[{{c, _Complex}},
Block[{z = c, n = 2, escapeRadius = 10 \[Pi], maxIterations = 100},
While[And[Abs[z] <= escapeRadius, n < maxIterations],
z = Cos[z] + c; n++]; n]]

Block[{center = {0.5527, 0.9435}, radius = 0.1},
DensityPlot[
cosineEscapeTime[x + I y], {x, center[[1]] - radius,
center[[2]] + radius}, PlotPoints -> 250, AspectRatio -> 1,
ColorFunction -> "TemperatureMap"]]


What could I do to improve the speed/time-efficiency of this code? Is there any reasonable way to parallelize it? (I'm running Mathematica 8 on an 8-core machine.)

edit Thanks all for the help so far. I wanted to post an update with what I'm seeing based on the answers so far and see if I get any further refinements before I accept an answer. Without going to hand-written C code and/or OpenCL/CUDA stuff, the best so far seems to be to use cosineEscapeTime as defined above, but replace the Block[...DensityPlot[]] with:

Block[{center = {0.5527, 0.9435}, radius = 0.1, n = 500},
Graphics[
Raster[Rescale@
ParallelTable[
cosineEscapeTime[x + I y],
ColorFunction -> "TemperatureMap"], ImageSize -> n]
]


Probably in large part because it parallelizes over my 8 cores, this runs in a little under 1 second versus about 27 seconds for my original code (based on AbsoluteTiming[]).

-

Use these 3 components: compile, C, parallel computing.

Also to speed up coloring instead of ArrayPlot use

Graphics[Raster[Rescale[...], ColorFunction -> "TemperatureMap"]]


In such cases Compile is essential. Compile to C with parallelization will speed it up even more, but you need to have a C compiler installed. Note difference for usage of C and parallelization may show for rather greater image resolution and more cores.

mandelComp =
Compile[{{c, _Complex}},
Module[{num = 1},
FixedPoint[(num++; #^2 + c) &, 0, 99,
SameTest -> (Re[#]^2 + Im[#]^2 >= 4 &)]; num],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True];

data = ParallelTable[
a + I b, {a, -.715, -.61, .0001}, {b, -.5, -.4, .0001}];

Graphics[Raster[Rescale[mandelComp[data]],
ColorFunction -> "TemperatureMap"], ImageSize -> 800, PlotRangePadding -> 0]


This is just a prototype - you can figure out a better coloring. Another way is to use LibraryFunction - we have Mandelbrot built in:

mlf = LibraryFunctionLoad["demo_numerical", "mandelbrot", {Complex},
Integer];
n = 501; samples =
Table[mlf[x + I y], {y, -1.25, 1.25, 2.5/(n - 1)}, {x, -2., .5,
2.5/(n - 1)}];
colormap =
Function[If[# == 0, {0., 0., 0.},  Part[r, #]]] /.
r -> RandomReal[1, {1000, 3}];
Graphics[Raster[Map[colormap, samples, {2}]], ImageSize -> 512]


Now, if you have a proper NVIDIA graphics card you can do some GPU computing with CUDA or OpenCL. I use OpenCL here because I got the source (from documentation btw):

Needs["OpenCLLink"]

src = "
__kernel void mandelbrot_kernel(__global mint * set, float zoom, \
float bailout, mint width, mint height) {
int xIndex = get_global_id(0);
int yIndex = get_global_id(1);
int ii;

float x0 = zoom*(width/3 - xIndex);
float y0 = zoom*(height/2 - yIndex);
float tmp, x = 0, y = 0;
float c;

if (xIndex < width && yIndex < height) {
for (ii = 0; (x*x+y*y <= bailout) && (ii < MAX_ITERATIONS); \
ii++) {
tmp = x*x - y*y +x0;
y = 2*x*y + y0;
x = tmp;
}
c = ii - log(log(sqrt(x*x + y*y)))/log(2.0);
if (ii == MAX_ITERATIONS) {
set[3*(xIndex + yIndex*width)] = 0;
set[3*(xIndex + yIndex*width) + 1] = 0;
set[3*(xIndex + yIndex*width) + 2] = 0;
} else {
set[3*(xIndex + yIndex*width)] = ii*c/4 + 20;
set[3*(xIndex + yIndex*width) + 1] = ii*c/4;
set[3*(xIndex + yIndex*width) + 2] = ii*c/4 + 5;
}
}
}
";

MandelbrotSet =
"mandelbrot_kernel", {{_Integer, _, "Output"}, "Float",
"Float", _Integer, _Integer}, {16, 16},
"Defines" -> {"MAX_ITERATIONS" -> 100}];

width = 2048;
height = 1024;
mem = OpenCLMemoryAllocate[Integer, {height, width, 3}];

res = MandelbrotSet[mem, 0.0017, 8.0, width, height, {width, height}];

Image[OpenCLMemoryGet[First[res]], "Byte"]


References:

Fractals CDF paper

Compile to C

Demonstrations

-
Great answer. I'd upvote but I'm out of votes for today. –  Mike Bantegui Jan 18 '12 at 5:41
Thank you, Mike, your answer is great too! This is all compiled from bits of Documentation. I'll post links too. –  Vitaliy Kaurov Jan 18 '12 at 6:09
If any answer is the "right" one (of which all of these are), this would be the one. Beautifully demonstrates usage of every concept mentioned, plus it has references to the documentation. –  Mike Bantegui Jan 19 '12 at 2:11
I played with this a little using AbsoluteTiming[Graphics[Raster[Rescale@ParallelTable[...],ColorFunction -> "TemperatureMap"]]] and it seems like your FixedPoint[] formulation of the escape-time function is slower and produces different results than my original While[] formulation, and adding , CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True inside the Compile[] of my While[] version doesn't make any significant difference. As you suggested, the Graphics[Raster[Rescale@...]] formulation is much faster than DensityPlot[]. –  Isaac Jan 19 '12 at 3:18
@Isaac at one point I wasted a day or so trying to find the fastest way to obtain the mandelbrot set. The fastest I could do was While[(iters < maxiter) && (Abs@z < 2),iters++;z = z^2 + c] compiled to C, and this was a bit faster than FixedPoint, Nest and other similar approaches. ie, I found the same as you. –  acl Jan 19 '12 at 11:25

Many plots can be speeded up by pre-generating the data set you want and then plotting the resulting list. In any case, it's not coincidence that Table and Plot have similar syntax, only that Plot does additional things like finding out the range to be displayed, the interpolation strength, and so forth. If you're already sure what kind of picture you want to get out, you may want to generate the data set by hand, like mentioned by Mike.

In your case, you can add just a few characters to the block part,

Block[{center = {0.5527, 0.9435}, radius = 0.1, data},
data = ParallelTable[
cosineEscapeTime[x + I y],
];
ListDensityPlot[
data,
AspectRatio -> 1,
ColorFunction -> "TemperatureMap"
]
]


This speeds up the process by more than an order of magnitude for me, and requires significantly less memory than the automatic plot. However, there's no interpolation in tricky corners of the plot here, so if you want a larger image you may want to change the parameters, i.e. decrease step size etc.

-
I'd like to second this. To explain a bit more: DensityPlot assumed a more or less smooth function (which this isn't), and tries to figure out where it need to sample more points for sufficient detail. Relevant options are PlotPoints and MaxRecursion. Both DensityPlot and ListDensityPlot will interpolate between sampled points, which can again be slow. I'd suggest using an Image instead of ListDensityPlot. –  Szabolcs Jan 18 '12 at 9:37

On version 7 you will need to use DistributeDefinitions:

{xc, yc} = {0.5527, 0.9435};
r = 0.1;
nx = 350;
ny = 350;

DistributeDefinitions[nx, ny, xc, yc, r];

pts = ParallelTable[
cosineEscapeTime[x + I y],
{x, xc - r, xc + r, (2 r)/nx},
{y, yc - r, yc + r, (2 r)/ny}
];


For faster plotting try using ArrayPlot:

ArrayPlot[Reverse@pts, ColorFunction -> "TemperatureMap"]


Or better still as Vitaliy recommends:

Graphics[Raster[Rescale@pts, ColorFunction -> "TemperatureMap"]]

-
Welcome to the party :) –  Mike Bantegui Jan 18 '12 at 5:07
With ColorFunction the fastest will be Raster: Graphics[Raster[..., ColorFunction -> "TemperatureMap"]] –  Vitaliy Kaurov Jan 18 '12 at 6:08

This might be an excellent candidate for ParallelTable;

MakeFractal[f_, nx_, ny_, {cx_, cy_}, {rx_, ry_}] := Module[{pts},
DistributeDefinitions[nx, ny, cx, cy, rx, ry, f];
pts = ParallelTable[f[x + I y],
{x, cx - rx, cx + rx, (2 rx)/nx},
{y, cy - ry, cy + ry, (2 ry)/ny}];
ArrayPlot[Reverse@pts, ColorFunction -> "TemperatureMap"]
]


Note that evaluation time is very fast, but plotting isn't.

You may wish to alternatively adjust the PlotPoints and MaxRecursion parameters of DensityPlot. MaxRecursion controls how far deep Mathematica goes for each plot point to determine the function value to use.

Too high of a value for MaxRecursion can lead to very long evaluation, especially with fractals.

Another way to help speed is to use CompilationTarget->"C" and RuntimeOptions->Speed. These sometimes provide a (small) speedup.

Here's some timing values (Note your value of PlotPoints->250 is a bit excessive):

(* 4.05 seconds on my machine *)
AbsoluteTiming[
Block[{center = {0.5527, 0.9435}, radius = 0.1},
DensityPlot[cosineEscapeTime[x + I y],
PlotPoints -> 120, MaxRecursion -> 2,
AspectRatio -> 1, ColorFunction -> "TemperatureMap"]]
]

(* 2.07 seconds on my machine *)
AbsoluteTiming[
pts = ParallelTable[
cosineEscapeTime[x + I y], {x, xc - r, xc + r, (2 r)/nx}, {y,
yc - r, yc + r, ( 2 r)/ny}];
ListDensityPlot[pts, AspectRatio -> 1,
ColorFunction -> "TemperatureMap"]
]


And individually for the ParallelTable approach:

(* 0.753 seconds on my machine *)
AbsoluteTiming[
pts = ParallelTable[
cosineEscapeTime[x + I y], {x, xc - r, xc + r, (2 r)/nx}, {y,
yc - r, yc + r, ( 2 r)/ny}];
]

(* 1.397 seconds on my machine *)
AbsoluteTiming[
ListDensityPlot[pts, AspectRatio -> 1,
ColorFunction -> "TemperatureMap"]
]


Using ArrayPlot as Mr.Wizard suggests is much faster:

(* 0.233 seconds on my machine *)
AbsoluteTiming[
ArrayPlot[Reverse@pts,
ColorFunction -> "TemperatureMap"
]
]


As you can see, plotting takes up a rather large amount of time.

-
Am I giving up anything significant in plotting an explicit table of points rather than having Mathematica/DensityPlot figure out for me dynamically which points to plot? –  Isaac Jan 18 '12 at 4:56
@Isaac: Well, you have to play around with the number of points which is a downside. In the other version you can just adjust PlotPoints and MaxRecursion and Mathematica "knows" where to evaluate next. Otherwise, you have the same "formatting" options that DensityPlot` has –  Mike Bantegui Jan 18 '12 at 4:58
You guys should try Raster here: Graphics[Raster[..., ColorFunction -> "TemperatureMap"]] With ColorFunction the Raster may be the fastest. –  Vitaliy Kaurov Jan 18 '12 at 6:33