# Is it possible to speed up ContourPlot on multi-core machines?

It is not very difficult to face a function for which ContourPlot works too slow. And it seems natural that this function can be parallelized well. Anyway, naive Parallelize@ContourPlot produces "ContourPlot[...] cannot be parallelized; "...

So, is it possible to parallelize ContourPlot?

-
Have you tried generating a discrete mesh and using ListContourPlot instead? The mesh does not have to be regular. – Verbeia Jan 30 '12 at 20:53
@Verbeia: actually I have not. You suggest generate points using ParallelMap and then plot them? – faleichik Jan 30 '12 at 20:59
Yes, that is what I had in mind. I wasn't sure if that was worth posting as an answer. – Verbeia Jan 30 '12 at 21:01
@verbeia I assume that would work, but it would mean you work on a regular grid. It doesn't use adaptive gridding that ContourPlot applies to increase sample density in rapidly changing spots. So the result can be rather different. – Sjoerd C. de Vries Jan 30 '12 at 21:02
@Sjoerd one could use Szabolcs's approach here: mathematica.stackexchange.com/questions/216/… – acl Jan 30 '12 at 21:04

I second @Verbeia's suggestion: compute the function on a mesh of points and use ListContourPlot. The disadvantage is that ListContourPlot has no adaptive sampling, so it'd be preferable if we could do our own adaptive sampling somehow. Adaptive sampling can give you a much better result while needing to compute the function in far less points---and the problem here is indeed computation time. So ContourPlot with its adaptive sampling might give a better result in less time on a single CPU than ListContourPlot will with a high resolution mesh computed on many CPUs.

The method I implemented there is usable (I am using it for something very similar to what you describe) but it is not nearly as good as ContourPlot's own. So one might still try to somehow make use of it. I'm quoting one suggestion I received from Leonid Shifrin there (in a comment):

You probably can control the DensityPlot, although not directly. Since it calls your function, you can simply Sow the values until some criteria (which you define) is violated (or satisfied). Then, you stop via throwing an exception, and catching it in the outer function, but still inside Reap. Alternatively, you could just start fooling DensityPlot by supplying faked values (perhaps, interpolated, or whatever), and it will stop by itself, I guess. Not sure this will work for you, but it may be worth trying.

I have not tried to implement this before, but I think it could work if your function is sufficiently smooth (which mine is definitely not, but yours may be).

Here's a quick sample implementation of how it could work:

First, let's define a sample function to plot:

fun[{x_, y_}] := 1/(1 + Exp[10 (Norm[{x, y}] - 3)])


Let's divide both the $x$ and $y$ axes into 5 parts on the interval $[0,5]$ and generate a mesh of points:

initialDivision = Range[0, 5];

points = N@Tuples[initialDivision, {2}];


Calculate function values on the intial mesh. This can be parallelized (just use ParallelMap)

values = fun /@ points;


This counter i will be used to control the maximal subdivisions in ContourPlot:

i = 0;


Now put the following code into a single cell, and evaluate it several times. Each time a finer and finer approximation will be computed. The points where function values have been computed will also be visualized. Note that I fixed the plot points in ContourPlot to force it to use the same initial mesh that I used, and I also fixed the number of contours.

if = Interpolation@ArrayFlatten[{{points, List /@ values}}]

{plot, {newpoints}} = Reap[
ContourPlot[if[x, y], {x, 0, 5}, {y, 0, 5},
Contours -> Range[0, 1, .1], MaxRecursion -> (++i),
PlotPoints -> Length[initialDivision],
EvaluationMonitor :> Sow[{x, y}]]
];
plot

newpoints = Complement[newpoints, points];
newvalues = fun /@ newpoints;  (* <-- this can be parallelized *)
points = Join[points, newpoints];
values = Join[values, newvalues];

Graphics[Point[points]]


After a few iterations the contour plot and the point mesh will look like this (note that the code above only plots the contours for the previous step, not the current results):

After 3 iterations, this method has computed the function value in 3809 points for this particular function.

Let's compare this with a plain ContourPlot using the same parameters:

ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5},
PlotPoints -> 6, MaxRecursion -> 3]


The quality of the plot is about the same with a plain ContourPlot as well.

How many points did the plain CountoutPlot use?

Reap[ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5}, PlotPoints -> 6,
MaxRecursion -> 3, EvaluationMonitor :> Sow[{x, y}]]][[2, 1]] // Length

(* ==> 3790 *)


It uses almost the same number of points, so if the bottleneck is computing f, the method I described is going to be almost as fast as ContourPlot on a single core, with the advantage that it is parallelizable for multiple cores.

The next step would be packaging this up into a self-contained function, but seeing how the quality improves step by step is also valuable as you can make decisions about when to stop calculating (and avoid excessive computation times).

I find it quite disappointing that all those nice and fast algorithms that plotting functions use (fast Voronoi cells, Delaunay trinagulation, adaptive sampling) are not directly accessible by users. We either have to use hacks to access these algorithms or reimplement them.

-
Very nice! I thought it could work, but now you've implemented this, which is much better than just an idea. +1. – Leonid Shifrin Jan 30 '12 at 22:21
I liked the idea. I'll try to make it neater and post it in a little while – Rojolalalalalalalalalalalalala Jan 31 '12 at 1:30
The first time I read this, I missed the fact that you were pulling the point info from ContourPlot instead of calculating it on your own. Clever. – rcollyer Jan 31 '12 at 16:51

One idea, which was originally appealing to me, was to simply split the plot domain in four quarter size pieces and have these pieces calculated in parallel using ParallelTable and then combined using Show. It appears though that the overhead of getting the graphics data back to Show is pretty large, so this only yields some extra speed if the function is computationally slow.

a = 36;
(parPlot = Show[
ParallelTable[
RegionPlot[
Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] +Sin[x] +Cos[y], \[Pi]] < \[Pi]/2,
{x, -a + i a, a + i a}, {y, -a + i a, a + i a},
PlotRange -> {{-a, a}, {-a, a}}, PlotPoints -> 50],
{i, 0, 1}, {j, 0, 1}]
]); // AbsoluteTiming

{25.5344605, Null}

parPlot


The original serial calculation:

a = 36;
g = RegionPlot[
Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] + Sin[x] +
Cos[y], \[Pi]] < \[Pi]/2, {x, -a, a}, {y, -a, a},
PlotPoints -> 100]; // AbsoluteTiming

(*
==> {51.2839332, Null}
*)

g


It's twice as slow, but it looks slightly better in places. RegionPlot is worse than ContourPlot in this sense.

-

Try generating a mesh of points using ParallelMap and then plot them using ListContourPlot.

You might find the techniques in this question on adaptive sampling for hard-to-compute functions useful.

-

The Presentations Application (which I am the author of) has a few routines to aid in the parallel processing of graphics.

Here is a serial contour plot of a complex function. (Presentations has routines for drawing complex functions directly with complex variables.)

<<Presentations

Module[{f = z \[Function] Sin[z/2], z, zmin = 0, zmax = 4 (1 + I)},
serial0 =
Draw2D[
{ComplexCartesianContour[f[z], {z, zmin, zmax}, Abs,
Contours -> Range[0, 4, 0.1],
ColorFunctionScaling -> False,
ColorFunction -> (ColorData["Rainbow"][Rescale[#, {0, 4}]] &),
PlotPoints -> 5,
MaxRecursion -> 4]},
Frame -> True,
ImageSize -> 300]
]; // AbsoluteTiming
serial0

{0.298017, Null}


The following statement parallel loads all the Presentations package routines to the processors.

ParallelLoadPresentations[]


The following routine will subdivide the complex plane into rectangular regions.

ComplexParallelPartitionDomain[{var, nRe, nIm, start, end, overlap}] will generate a list of nRe x nIm complex subiterators in the form {{var, start1, end1},...} for partitioning complex drawing object domains. The parameter overlap has a default value of 0 but otherwise may be used to overlap each subdomain except on the right and top boundaries.

Here, since I have six processors, we will use

ComplexParallelPartitionDomain[{z, 2, 3, 0, 4 (1 + I), 0}]
{{z, 0, 2 + (4 I)/3}, {z, (4 I)/3, 2 + (8 I)/3},
{z, (8 I)/3, 2 + 4 I}, {z, 2, 4 + (4 I)/3},
{z, 2 + (4 I)/3, 4 + (8 I)/3}, {z, 2 + (8 I)/3, 4 + 4 I}}


I find that a convenient way to do parallel processing is to use the ParallelSubmit and WaitAll mechanism. Here is the code for the parallel processing version.

Module[{f = z \[Function] Sin[z/2], z, zmin = 0, zmax = 4 (1 + I)},
par1 =
Draw2D[
{WaitAll[
ParallelSubmit[{f},
ComplexCartesianContour[f[z], #, Abs,
Contours -> Range[0, 4, 0.1],
ColorFunctionScaling -> False,

ColorFunction -> (ColorData["Rainbow"][
Rescale[#, {0, 3.7}]] &),
PlotPoints -> 6,
MaxRecursion -> 3]] & /@
ComplexParallelPartitionDomain[{z, 2, 3, 0, 4 (1 + I), 0}]
]},
Frame -> True,
ImageSize -> 300]
]; // AbsoluteTiming
par1
{0.101006, Null}


Next I will try Sjoerd's example above.

region = {x, y} \[Function] (Mod[Sqrt[x^2 + y^2] - 7/2 ArcTan[x, y] + Sin[x] +
Cos[y], \[Pi]] < \[Pi]/2)

(serial1 = RegionPlot[region[x, y], {x, -35, 35}, {y, -35, 35},
PlotPoints -> 20,
MaxRecursion -> 5]); // AbsoluteTiming
serial1
{44.546548, Null}


There is a similar partitioning function for x-y domains.

ParallelPartitionDomain[{x, 2, -35, 35, 0}, {y, 3, -35, 35, 0}]
{{{x, -35, 0}, {y, -35, -(35/3)}}, {{x, -35, 0}, {y, -(35/3), 35/3}},
{{x, -35, 0}, {y, 35/3, 35}}, {{x, 0, 35}, {y, -35, -(35/3)}},
{{x, 0, 35}, {y, -(35/3), 35/3}}, {{x, 0, 35}, {y, 35/3, 35}}}


There is one problem with partitioning RegionPlots in this way. If we want the automatic outlining of the region, then with partitioning we obtain some boundary lines on the edges of the partition. These are also on the edges of the serial plot but much less noticable. Drawing the boundary lines definitely improves the image. One way to handle this is to use a darker color for the region and the same color for the boundary lines.

par2 =
Draw2D[
{WaitAll[
ParallelSubmit[{region},
RegionDraw[region[x, y], Evaluate[Sequence @@ #],
MaxRecursion -> 3, PlotPoints -> 20,
PlotStyle -> Darker@Orange,
BoundaryStyle -> Darker@Orange]] & /@
ParallelPartitionDomain[{x, 2, -35, 35, 0}, {y, 3, -35, 35, 0}]
]},
Frame -> True,
ImageSize -> 350]; // AbsoluteTiming
par2
{2.196126, Null}


Notice that we obtained a speedup of a factor of 20, despite using only six processors, and the Parallel Kernel Status showing a speedup of 5.45. This is because we were able to use much less recursion to obtain the same quality of plot. These comparisions are a bit tricky because we have to compare the two plots by eye and we have two variables, PlotPoints and MaxRecursion` to adjust. Nevertheless, a smaller piece of a plot will generally be topologically simpler and we can often achieve speedups that exceed the number of processors.

-