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

Question

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?

Answer 1

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.

Adaptive sampling is what I asked about (and solved) here: Adaptive sampling for slow to compute functions in 2D

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

The following version is a mixture between Sjoerd's answer and Szabolcs answer. I'm using the CountorPlot example, and not the RegionPlot example, since I believe there's a bug on the EvaluationMonitor option of the RegionPlot (and also, because there's no such thing as a ListRegionPlot).

First, I start by making the function a little slower, just to be on a "real" scenario.

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

The method is straightforward (and if someone is willing to make my code simpler, please go ahead). We just divide the 2D space.

Let's start by the serial case for speed benchmarking:

AbsoluteTiming@ContourPlot[fun[{x, y}], {x, 0, 5}, {y, 0, 5},
 Contours -> Range[0, 1, .1], ColorFunctionScaling -> False, 
 MaxRecursion -> 3, PlotPoints -> {6, 6}, Frame -> False, 
 ImageSize -> 200, ClippingStyle -> Automatic
]

Now, we parallelize it:

coreMesh = {10, 10};
xRange = {0, 5};
yRange = {0, 5};
totalPlotPoints = {6, 6};
GraphicsGrid[allPlots = ParallelTable[
 ContourPlot[
  fun[{x, y}],
   {x, (xI - 1) (xRange[[2]] - xRange[[1]])/coreMesh[[1]],
    xI (xRange[[2]] - xRange[[1]])/coreMesh[[1]]},
   {y, (yI - 1) (yRange[[2]] - yRange[[1]])/coreMesh[[2]],
    yI (yRange[[2]] - yRange[[1]])/coreMesh[[2]]}, 
 Contours -> Range[0, 1, .1], ColorFunctionScaling -> False, 
 MaxRecursion -> 3, 
 PlotPoints -> {Max[2, Floor[totalPlotPoints[[1]]/coreMesh[[1]]]], 
  Max[2, Floor[totalPlotPoints[[2]]/coreMesh[[2]]]]}, 
 Frame -> False, ImageSize -> Floor[200/coreMesh[[1]]], 
 ClippingStyle -> Automatic], {yI, coreMesh[[2]], 1, -1}, {xI, 1, 
 coreMesh[[1]]}], Spacings -> 0]

Show[allPlots, PlotRange -> All, ImageSize -> 200]

I get similar artefacts as the ones identified by Sjoerd for the RegionPlot.

On this speed testing case, I'm going to divide on a coarser mesh (and I slightly improved the code readability).

coreMesh = {4, 4};
xRange = {0, 5};
yRange = {0, 5};
{stepX, stepY} = {(xRange[[2]] - xRange[[1]])/coreMesh[[1]],
 (yRange[[2]] - yRange[[1]])/coreMesh[[2]]};
totalPlotPoints = {6, 6};
{time, allPoints} =AbsoluteTiming@Flatten[ParallelTable[Reap[
 ContourPlot[z = fun[{x, y}],
  {x, (xI - 1) stepX, xI*stepX}, {y, (yI - 1) stepY, yI*stepY}, 
  Contours -> Range[0, 1, .1], MaxRecursion -> 3, 
  PlotPoints -> {Max[2,Floor[totalPlotPoints[[1]]/coreMesh[[1]]]], 
    Max[2, Floor[totalPlotPoints[[2]]/coreMesh[[2]]]]}, 
  EvaluationMonitor :> Sow[{x, y, z}]]][[2, 1]],
{yI, coreMesh[[2]], 1, -1}, {xI, 1, coreMesh[[1]]}], 2];



{time, Length[allPoints]}
(*{1.82858, 4472}*)

A little more points, and yet, less than half the time (on 8 threads...). I got the best performance on a 3 by 3 mesh, but since the amount of points was smaller, the comparison wouldn't be honest.

Graphics[Point[allPoints[[All, {1, 2}]]], ImageSize -> Small]

Putting all points together:

AbsoluteTiming@ListContourPlot[allPoints, Contours -> Range[0, 1, .1], 
 ColorFunctionScaling -> False, Frame -> False, ImageSize -> 200
]

So, on total, 2.1 seconds, instead of 5.4 seconds. Not huge... but this method is very susceptible to non optimized work balances, since all the complexity may end up being on just one kernel (which, on this example, is not very far from the truth).