Speeding up this fractal-generating code

Question

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[[1]] + radius}, {y, center[[2]] - 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],
      {y, center[[2]] - radius, center[[2]] + radius, 2 radius/(n - 1)}, 
      {x, center[[1]] - radius, center[[1]] + radius, 2 radius/(n - 1)}], 
   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[]).

Answer 1

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 = 
  OpenCLFunctionLoad[src, 
   "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

LibraryFunction

OpenCL

Demonstrations

Answer 2

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],
   {x, center[[1]] - radius, center[[1]] + radius},
   {y, center[[2]] - radius, center[[2]] + radius},
   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.

Answer 3

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],
        {x, center[[1]] - radius, center[[1]] + radius, 2 radius/250},
        {y, center[[2]] - radius, center[[2]] + radius, 2 radius/250}
    ];
    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.

Answer 4

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"]]