# Adaptive sampling for slow to compute functions in 2D

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:

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


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

Question:

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.

• 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. Jan 18, 2012 at 21:24
• 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.) Jan 19, 2012 at 1:40
• 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? Jan 19, 2012 at 4:31
• 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. Jan 19, 2012 at 12:02
• This is related to this question. Oct 6, 2012 at 17:30

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_] :=
Module[
{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];

InterpolationOrder -> 0, Mesh -> All, ColorFunction -> "MintColors",
Epilog -> {PointSize[Large], Point[points], Red, Point[newpoints]}]

points = Join[points, newpoints];


The result after several iterations:

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


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


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

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):

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!

• 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.
– user21
Jan 20, 2012 at 7:48
• 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. Mar 31, 2012 at 6:27
• @Szabolcs : have you made any progress for this step by step DensityPlot, based on curvature detection ? Thank you for the current version anyway 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 *)
Do[

(* 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]];
Print["---------"];
*)

xnew = Quiet[
Reap[
NIntegrate[
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],
fake[x]],MemberQ[known,x]}];
*)
Sow[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[#] &,
xnew
, Method -> "FinestGrained"];

(* memoize those values that were computed in parallel *)
Do[
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[
NIntegrate[
mem[x],
{x, xmin, xmax},
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
MinRecursion -> r,
MaxRecursion -> r
(*,EvaluationMonitor:>Print[{"2nd: ", x,mem[x]}]*)
]
];

(* increment up to maxr *)
, {r, 0, maxr}];
res
];


Define a slow function:

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


Integrate the normal way:

AbsoluteTiming[
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:

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

{35.0655, -0.100505}
`