0
$\begingroup$

I have function that is costly to evaluate, fun[x,y], that takes two real arguments and returns True or False. I am trying to visualize in the x,y-plane the region where the function returns True.

I took a sample of points on a 2-dimensional grid on the x,y-plane, so I generated a list of the form:

{{x1,y1,True}, {x2,y2,True}, {x3,y3,False}, ...}

consisting of the coordinates of each point and the return value of the function.

Using ArrayPlot I can turn this into a visualization of the region where the function is true, but it is ugly because it has very straight edges. Is there a way to obtain a smoothed boundary?

Here is an example

fun[x_, y_] := (Boole[x^2 <= y])
data = Flatten[
   Table[{x, y, fun[x, y]}, {x, 0, 1, 0.1}, {y, 0, 1, 0.1}], {1, 2}];
ListPlot[{Cases[data, {_, _, 1}], Cases[data, {_, _, 0}]}[[All, 
  All, {1, 2}]]]

enter image description here

Of course, fun[x,y] is a dummy example, in reality it is expensive to evaluate. It would be nice if could get a smooth region boundary.

$\endgroup$
  • $\begingroup$ Can you post an example? $\endgroup$ – Henrik Schumacher Oct 22 '17 at 14:05
  • $\begingroup$ @HenrikSchumacher I added an example. $\endgroup$ – becko Oct 22 '17 at 14:20
  • $\begingroup$ Sigh As code? $\endgroup$ – Henrik Schumacher Oct 22 '17 at 14:21
  • $\begingroup$ @HenrikSchumacher Done $\endgroup$ – becko Oct 22 '17 at 14:22
  • $\begingroup$ Hm. Can we assume that the value of fun depends on the fulfillment of few rather smooth inequalities? $\endgroup$ – Henrik Schumacher Oct 22 '17 at 14:27
5
$\begingroup$

Then what about this: I interpolate with a BezierFunction and plot it's superlevel set:

dx = 0.1;
fun[x_, y_] := (Boole[x^2 <= y])
data = Table[{x, y, fun[x, y]}, {x, 0, 1, dx}, {y, 0, 1, dx}];
f = BezierFunction[data];
Show[
 Quiet@RegionPlot[f[x, y][[3]] >= .5, {x, 0, 1}, {y, 0, 1}],
 ListPlot[{
    Cases[Flatten[data, 1], {_, _, 1}],
    Cases[Flatten[data, 1], {_, _, 0}]}[[All, All, {1, 2}]]
  ]
 ]

enter image description here

BezierFunction only accepts arguments in the unit interval/square/cube. At least for rectangular domains, we can rescale and shift coordinates to obtain a slightly more flexible tool:

maxx = 3;
maxy = 2;
minx = miny = -1;
dx = 0.1;
dy = 0.2;
fun[x_, y_] := (Boole[x^2 <= y])
data = Table[{x, y, fun[x, y]}, {x, minx, maxx, dx}, {y, miny, maxy, dy}];
f = BezierFunction[data.DiagonalMatrix[{1./(maxx - minx), 1./(maxy - minx), 1.}]];
Show[
 Quiet@RegionPlot[ f[(x - minx)/(maxx - minx) , (y - miny)/(maxy - miny) ][[3]] >= .5, {x, minx, maxx}, {y, miny, maxy}], 
 ListPlot[{Cases[Flatten[data, 1], {_, _, 1}], Cases[Flatten[data, 1], {_, _, 0}]}[[All, All, {1, 2}]]]
]

enter image description here

However, I do not know how to get rid of the requirement of even spacings of data points... Maybe one can also experiment with Interpolation to obtain suitable functions f. My experience with Interpolation was that the solutions were a bit too tight to data. Using smoothing splines instead may also help.

$\endgroup$
  • $\begingroup$ That looks nice. I'll try it on my data and let you know. Thanks! $\endgroup$ – becko Oct 22 '17 at 15:13
  • $\begingroup$ How do the $u,v$ parameters on the Bezier relate to my original $x,y$? In this example the $x,y$ go from 0 to 1, but in general this is not the case. $\endgroup$ – becko Oct 26 '17 at 7:49
  • $\begingroup$ Thanks for the edit. Yes, the even spacing is too tight a requirement. $\endgroup$ – becko Oct 26 '17 at 12:13
  • $\begingroup$ Okay, I posted a new approach... $\endgroup$ – Henrik Schumacher Oct 26 '17 at 12:30
  • $\begingroup$ Yes! I was looking at it. Massive! Thanks a lot for your help! $\endgroup$ – becko Oct 26 '17 at 12:32
4
$\begingroup$

One way is to use the function as RegionFunction in ContourPlot. In the example below, I desined expensivefunction such that it keeps track about the evaluation points. You may use the options PlotPoints and MaxRecursion to control the amount of function evaluations and the accuracy of the sampling of the domain.

expensivefunction[x_?NumericQ, y_?NumericQ] := (
   Internal`StuffBag[pts, {x, y}];
   x^2 + 0.1 Sin[3 Pi x] + y^2 <= 1.
   );
f = {x, y} \[Function] expensivefunction[x, y];

pts = Internal`Bag[{}];
g = ContourPlot[1, {x, -1.2, 1.2}, {y, -1.2, 1.2},
   RegionFunction -> f,
   Mesh -> All,
   BoundaryStyle -> Black,
   PlotPoints -> 2,
   MaxRecursion -> 2
   ];
Show[{g, Graphics[Point[Internal`BagPart[pts, All]]]}]

enter image description here

With higher values for MaxRecursion you get better resolution, but of course, also mor points to evaluate your function on. The following is the same as above but with MaxRecursion -> 5:

enter image description here

For your example, one could try this

expensivefunction[x_?NumericQ, y_?NumericQ] := (Internal`StuffBag[pts, {x, y}];x^2 - y <= 0.);
f = {x, y} \[Function] expensivefunction[x, y];
pts = Internal`Bag[{}];
g = ContourPlot[1., {x, 0, 1}, {y, 0, 1},
   Mesh -> All,
   PlotPoints -> {10, 5},
   MaxRecursion -> 0,
   RegionFunction -> f
   ];
Show[{g, Graphics[Point[Internal`BagPart[pts, All]]]}]
Internal`BagLength[pts]

Unfortunately, it seems to be impossible to tell ContourPlot exactly from which points to sample...

enter image description here

$\endgroup$
  • $\begingroup$ This is a nice approach. But I would like to control the sampled points manually, if possible. So I generate the points myself, then the question is how I can plot a nice region out of it. $\endgroup$ – becko Oct 22 '17 at 14:22
3
$\begingroup$

And yet another approach. First, I try to detect the region of the interface between 0- and 1-valued points. Afterwards I use the minimimizer of the Dirichlet energy with respect 0- and 1-valued boundary conditions in this region; the final interface will be the 1/2-levelset.

First, the data points:

maxx = 3;
maxy = 2;
minx = miny = -1;
fun[x_, y_] := (Boole[x^2 <= y]);
nn = 600;
data0 = Transpose[{RandomReal[{minx, maxx}, nn], RandomReal[{miny, maxy}, nn]}];
data1 = Join[data0, Transpose[{fun @@@ data0}], 2];
listplot = Show[Graphics[], ListPlot[{Cases[data1, {_, _, 1}], Cases[data1, {_, _, 0}]}[[All, All, {1, 2}]], 
   PlotStyle -> {ColorData[97][3], ColorData[97][4]}]
  ]

enter image description here

Next, we coarsly detext the interface region:

Needs["TriangleLink`"]
R = TriangleCreate[];
pts = data1[[All, 1 ;; 2]];
TriangleSetPoints[R, pts];
S = TriangleTriangulate[R, "a1.0"];
pts = TriangleGetPoints[S];
faces = TriangleGetElements[S];
pos = Flatten[Position[Equal @@@ Partition[data1[[Flatten[faces], 3]], 3], False, 1]];
nfaces = faces[[pos]];
plist = Union @@ nfaces;
npts = pts[[plist]];
lookuptable = AssociationThread[plist -> Range[Length[plist]]];
nfaces = Partition[Lookup[lookuptable, Flatten[nfaces]], 3];
nullpos = Flatten[Position[data1[[plist, 3]], 0, 1]];
onepos = Flatten[Position[data1[[plist, 3]], 1, 1]];
R = MeshRegion[npts, Polygon[nfaces]];
With[{edges = Developer`ToPackedArray[MeshCells[R, 1][[All, 1]]]},
  R0 = MeshRegion[MeshCoordinates[R], Line[Select[edges, SubsetQ[nullpos, #] &]]];
  R1 = MeshRegion[MeshCoordinates[R], Line[Select[edges, SubsetQ[onepos, #] &]]];
  ];
GraphicsRow[{Show[listplot, R], Show[listplot, R0, R1]}, ImageSize -> Large]

enter image description here

Finally, we solve a boundary value problem. This can be done easier with the NDSolve facilities, but as I am not that much into these details, I do it a bit more verbatim.

S = DiscretizeRegion[R, MaxCellMeasure -> 0.0001];
n = Length[MeshCoordinates[S]];
nullpos = Flatten[Position[RegionMember[R0]@MeshCoordinates[S], True, 1]];
onepos = Flatten[Position[RegionMember[R1]@MeshCoordinates[S], True, 1]];
plist = Join[onepos, nullpos];
Module[{vd, sd, cdata, mdata, dpde, dbc, load, damping, bcdata, y, x, u},
  Needs["NDSolve`FEM`"];
  Rdiscr = ToElementMesh[
    "Coordinates" -> MeshCoordinates[S],
    "MeshElements" -> {TriangleElement[MeshCells[S, 2][[All, 1]]]},
    "MeshOrder" -> 1,
    "NodeReordering" -> False];
  vd = NDSolve`VariableData[{"DependentVariables","Space"} -> {{u}, {x, y}}]; 
  sd = NDSolve`SolutionData[{"Space"} -> {Rdiscr}]; 
  cdata = InitializePDECoefficients[vd, sd,
    "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
    "MassCoefficients" -> {{1}}, "LoadCoefficients" -> {{0}}
    ];
  mdata = InitializePDEMethodData[vd, sd];
  dpde = DiscretizePDE[cdata, mdata, sd];
  dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd]; {load, 
    stiffness, damping, mass} = dpde["All"];
  ];
costraints = SparseArray[Transpose[{Range[Length[plist]], plist}] -> 1., {Length[plist], n}, 0.];
L = ArrayFlatten[{{stiffness, Transpose[costraints]}, {costraints, 0.}}];
b = Join[ConstantArray[0., n], ConstantArray[1., Length[onepos]], ConstantArray[0., Length[nullpos]]];
x = LinearSolve[L, b][[1 ;; n]];
f = ElementMeshInterpolation[{Rdiscr}, x];
GraphicsGrid[{
  {HighlightMesh[S, {0, nullpos}], HighlightMesh[S, {0, onepos}]},
  {Plot3D[f[u, v], Element[{u, v}, S]],
   Show[ContourPlot[f[u, v] == 1/2, Element[{u, v}, S]], listplot]
   }}
 , ImageSize -> Large]

enter image description here

This should also work for quite general domains.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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