EDIT: Although I have posted an answer based on my current progress, this in incomplete. Please see the "open issues" section in the answer.

Most plotting functions in Mathematica adjust the sampling density dynamically, based on the slope of the function:

 1 + Exp[10 (Sqrt[x^2 + y^2] - 3)]), {x, 0, 5}, {y, 0, 5}, 
 Mesh -> All, PlotRange -> All]

Mathematica graphics

Unfortunately the internal algorithm that does this is not directly callable.


I have a slow to compute 2D function (takes up to 10-40 seconds for a single point even though it's in optimized C++ called through LibraryLink) that I need to sample. How can I sample it adaptively, in a convenient and controlled way?

Since the function is so slow to compute,

  • I'd like to be able to take the existing points, and refine them more if needed (i.e. continue the computation using the existing results)

  • I can't use DensityPlot directly because I can't control how many points it is going to compute, and I need an upper limit on that (i.e. an upper limit on computation time). It also can't be interrupted and continued later.

I am looking for

  1. effective methods to do this

  2. implementations in Mathematica (either as an answer or pointers to libraries)

The messy details:

I am trying to compute a phase diagram and map the boundaries between the phases precisely. So I don't need the function value everywhere, only where it very suddenly drops. The function is either of magnitude 1 (say, between $0.1 \div 1$), or very small (close to zero).

The function is computed using Monte Carlo methods, so at a small scale it doesn't appear smooth, and I might get inconsistent results close to the phase boundaries on subsequent runs of the function.

This should give you an idea of what sort of function I'd like to apply this to, which might be important when choosing a method.

  • 1
    $\begingroup$ Why not threshold to get 0-1 values? Could do a bit map, say. Then do something more extensive in (typically thin) regions that contain both 0s and 1s. $\endgroup$ Commented Jan 18, 2012 at 21:24
  • $\begingroup$ I'd say you have to code this manually, i.e. generate a crude grid, define points between the existing function values, convolve the thing with some flattener and find out whether there's still a lot of fluctuation going on, and then deciding to calculate the value for those fluctuating intermediate points. I mean in the end what you've got is more of a data list generation than a plotting problem. (I don't think you can do more with Plot than changing MaxRecursion and PlotPoints.) $\endgroup$
    – David
    Commented Jan 19, 2012 at 1:40
  • $\begingroup$ If the function is computed with Monte Carlo, do you actually use ListDensityPlot for a set of points then? Or somehow the function has symbolic form? $\endgroup$ Commented Jan 19, 2012 at 4:31
  • 1
    $\begingroup$ This SciComp question may be related. It's not about Mathematica code, but some of the links might be useful resources if you have to code it yourself. $\endgroup$
    – David Z
    Commented Jan 19, 2012 at 12:02
  • 1
    $\begingroup$ This is related to this question. $\endgroup$
    – adavid
    Commented Oct 6, 2012 at 17:30

2 Answers 2


Update: I described an alternative approach based on built in plotting functions in this answer. That approach is not very practical here though because I need to be able to handle points at arbitrary positions while built in functions work with a rectangle-based mesh. I am still looking for improvements.

I came up with this very naive approach and implementation (I know that the implementation is not optimal at all):

First let's define a test function (same one as in the question):

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

These functions will subdivide lines in the Delaunay triangulation of the points if 1. the points are further apart than a threshold (i.e. the resolution is controlled) and 2. the function values in the two points also differ by more than another threshold.

<< ComputationalGeometry`

makeLines[tri_] := Union[Sort /@ Flatten[Thread /@ tri, 1]]

subdivision[points_, values_, valueThreshold_, distanceThreshold_] :=
  {tri, lines, linesToDivide},
  tri = DelaunayTriangulation[points];
  lines = makeLines[tri];
  linesToDivide = 
   Pick[lines, (Abs[values[[#1]] - values[[#2]]] > valueThreshold && 
        Norm[points[[#1]] - points[[#2]]] > distanceThreshold ) & @@@ lines];
  Mean /@ (linesToDivide /. n_Integer :> points[[n]])

Let's define an initial point grid to compute the function in:

points = Tuples[Range[0, 5, 1], 2];

We can iterate this function to add more and more points and recursively subdivide the grid (evaluate the following commands together repreatedly):

values = fun /@ N[points];
newpoints = subdivision[points, values, 0.1, 0.1];

ListDensityPlot[Flatten /@ Thread[{points, values}], 
 InterpolationOrder -> 0, Mesh -> All, ColorFunction -> "MintColors", 
 Epilog -> {PointSize[Large], Point[points], Red, Point[newpoints]}]

points = Join[points, newpoints];

Mathematica graphics Mathematica graphics Mathematica graphics Mathematica graphics

The result after several iterations:

values = fun /@ N[points];
ListDensityPlot[Flatten /@ Thread[{points, values}], 
 InterpolationOrder -> 0, Mesh -> All, ColorFunction -> "MintColors"]

Mathematica graphics

Open question: My aim is to minimize the number of points I need to compute while getting a precise approximation. This is probably not the best subdivision method for it. What are some easy-to-implement better methods?

I think ideally the decision for refining the grid should be made based on some sort of curvature. Take for example the following function:

ContourPlot[Erf[1/(1 + 20 x^2) - y], {x, -3, 3}, {y, -3, 3}]

Mathematica graphics

Using a valueThreshold of 0.3 and distanceThreshold of 0.1, and a starting grid with a spacing of 0.5 produces this:

Mathematica graphics

Let's turn on interpolation (because I can't turn interpolation off in DensityPlot) and compare it with a DensityPlot made using similar options (PlotPoints -> 12, MaxRecursion -> 15):

Mathematica graphics

The curvature-based DensityPlot (right) is clearly much better. Furthermore, my method won't properly detect "fjord-like" structures (similar to the one in this example). It tend to "jump" over them, this is why some artefacts are visible in the middle of the plot.

Thanks to @ruebenko for the hints and ideas he sent me!

  • $\begingroup$ very nice. Here is an idea for an improvement. Compute an initial mesh. Subdivide everything and compute the new points. Then it should be possible estimate if and where the next subdivision is necessary. You could also look at the gradient of change. $\endgroup$
    – user21
    Commented Jan 20, 2012 at 7:48
  • $\begingroup$ If your going to mesh a symmetric function I'd expect the result to be asymmetric as well. In your last mesh plot the result is clearly asymetric, so perhaps there's another angle of attack. $\endgroup$ Commented Mar 31, 2012 at 6:27
  • 1
    $\begingroup$ @Szabolcs : have you made any progress for this step by step DensityPlot, based on curvature detection ? Thank you for the current version anyway $\endgroup$
    – sekisushai
    Commented Jul 15, 2016 at 14:37

I have cooked up a routine that performs adaptive sampling in parallel for integration. With a little tweaking, it should work for multidimensional plotting routines such as this--I think they are related problems as "parallel adaptive sampling" which Mathematica uses in several areas such as plotting, polynomial fitting, and integration.

My use case is very expensive integrand. The idea is to trick NIntegrate into giving yielding where it would like to integrate, but not actually doing the integration right away, but rather deferring it to be done in parallel batch. These values are then memoized, and provided to a second invocation of NIntegrate, which merely looks up the (pre-computed) values that it needs.

(It was made in desperation and probably is a horrible implementation it but I hope the idea/solution are sound.)

ParallelNInt[f_, {xmin_, xmax_}, maxr_] := Module[
   {memoize, mtemp, res, fake, mem, xnew, known, minr, integrand},
   (* keep a list of known x values *)
   known = {};
   (* ensure numeric evaluation inside NIntegrate *)
   mem[x_?NumericQ] := memoize[x];
   fake[x_?NumericQ] := 0;
   (* increment MaxRecursion from 0 to maxr *)
    (* First Invocation: "trick" NIntegrate 
    so that it either looks up values from previous
    lower recursive levels, or quickly yields a fake value *)
    Print["Recursion r = " <> ToString[r]];
    xnew = Quiet[
         If[MemberQ[known, x], mem[x], fake[x]],
         {x, xmin, xmax},
         Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
         MinRecursion -> r,
         MaxRecursion -> r,
         EvaluationMonitor :> (
           Print[{"1st: ",x,If[MemberQ[known,x],mem[x],
        ][[2, 1]]
    (* return the previous Second Invocation result if no new data *)

        If[Length[xnew] == 0, Return[res]];
    (* for values that we have not seen, evaluate in parallel *)
    xnew = Complement[xnew, known];
    mtemp = ParallelMap[
      f[#] &,
      , Method -> "FinestGrained"];
    (* memoize those values that were computed in parallel *)
     memoize[xnew[[i]]] = mtemp[[i]]
     , {i, Length[mtemp]}];
    (* append those x-values to the known list *)
    known = Join[known, xnew];
    (* Second Invocation: NIntegrate with memoized values *)
    res = Quiet[
       {x, xmin, xmax},
       Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
       MinRecursion -> r,
       MaxRecursion -> r
       (*,EvaluationMonitor:>Print[{"2nd: ", x,mem[x]}]*)
    (* increment up to maxr *)
    , {r, 0, maxr}];

Define a slow function:

f[x_?NumericQ] := (Pause[5]; 1/x Cos[Log[x]/x]);

Integrate the normal way:

 NIntegrate[f[x], {x, 0.1, 0.2}, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  MinRecursion -> 1, MaxRecursion -> 1]

{95.0225, -0.100505}

Integrate with this parallel function:

 ParallelNInt[f, {0.1, 0.2}, 1]

{35.0655, -0.100505}

Addressing a few points in your question:

  • I use "LocalAdaptive" because I believe that better-reuses data from previous recursions.
  • One may use an integration technique to generate data and then use a List___Plot command. Or, adapt this code to make graphics instead of integrations. Anything that can Reap/Sow is amenable.
  • Uses native Mathematica routines and no 3rd party packages (such as Cubature)
  • Perhaps can be used with NIntegrate and Monte Carlo methods, I am not sure if they obey MaxRecursion or just how terribly sequential they are.
  • Speed-up is merely through parallelization and it is easy to check against serial versions since this solution is supposed to be a order-of-operations difference tops.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.