# 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. – Daniel Lichtblau Jan 18 '12 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.) – David Jan 19 '12 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? – Vitaliy Kaurov Jan 19 '12 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. – David Z Jan 19 '12 at 12:02
• This is related to this question. – adavid Oct 6 '12 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 '12 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. – Sjoerd C. de Vries Mar 31 '12 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 – Pierre Jul 15 '16 at 14:37