Speeding up this fractal-generating code

63

32

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:

graphic

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[]).

Isaac

Posted 2012-01-18T04:39:46.967

Reputation: 2 989

Answers

69

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]

enter image description here

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]

enter image description here

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

enter image description here

References:

Fractals CDF paper

Compile to C

LibraryFunction

OpenCL

Demonstrations

Vitaliy Kaurov

Posted 2012-01-18T04:39:46.967

Reputation: 66 672

Strange, Graphics[Raster[Rescale[...]]] returns a white image in Mathematica 10 on OS X. ArrayPlot[Reverse@...] works fine. – shrx – 2014-01-30T16:08:47.203

Great answer. I'd upvote but I'm out of votes for today. – Mike Bailey – 2012-01-18T05:41:46.877

Thank you, Mike, your answer is great too! This is all compiled from bits of Documentation. I'll post links too. – Vitaliy Kaurov – 2012-01-18T06:09:37.317

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 Bailey – 2012-01-19T02:11:07.837

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 – 2012-01-19T03:18:41.137

With further experimentation, the FixedPoint[] version does yield the same results as my While[] version (I missed something that caused it not to before), but it is still about 20% slower than the While[] version. – Isaac – 2012-01-19T03:50:49.707

@Isaac My guess is that Parallelization and C would start to show advantage for the larger image sizes and more cores present. I just wanted to share different approaches documented in M. Thanks for trying these out ;-) – Vitaliy Kaurov – 2012-01-19T04:05:32.267

1@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 – 2012-01-19T11:25:10.563

33

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.

Mike Bailey

Posted 2012-01-18T04:39:46.967

Reputation: 1 875

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 – 2012-01-18T04:56:12.120

@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 Bailey – 2012-01-18T04:58:41.487

1You guys should try Raster here: Graphics[Raster[..., ColorFunction -> "TemperatureMap"]] With ColorFunction the Raster may be the fastest. – Vitaliy Kaurov – 2012-01-18T06:33:51.743

15

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.

David

Posted 2012-01-18T04:39:46.967

Reputation: 14 421

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 – 2012-01-18T09:37:41.947

11

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

Mr.Wizard

Posted 2012-01-18T04:39:46.967

Reputation: 259 163

Welcome to the party :) – Mike Bailey – 2012-01-18T05:07:50.910

1With ColorFunction the fastest will be Raster: Graphics[Raster[..., ColorFunction -> "TemperatureMap"]] – Vitaliy Kaurov – 2012-01-18T06:08:23.163