Test if argument is inside domain of InterpolatingFunction

Question

The documentation says In standard output format, only the domain element of an InterpolatingFunction object is printed explicitly. The remaining elements are indicated by <>

But this is not the same as the actual domain where it doesn't use extrapolated values, see this example with a $\mathbb{R}^2 \rightarrow \mathbb{R}$ function:

SeedRandom[1];
dat = RandomReal[1, {9, 3}];
paramdat = Map[{{#[[1]], #[[2]]}, #[[3]]} &, dat];

ip = Interpolation[paramdat, InterpolationOrder -> 1]
(* InterpolatingFunction[{{0.128821,0.825163},{0.11142,0.925275}},<>] *)

ip["Domain"]
(* {{0.128821, 0.825163}, {0.11142, 0.925275}} *)

{x, y} = {0.12887062309893396`, 0.11147780180592484`};

ip[[1, 1, 1]] < x < ip[[1, 1, 2]] && ip[[1, 2, 1]] < y < ip[[1, 2, 2]]
(* True *)

ip[x, y];
(* InterpolatingFunction::dmval:
   "Input value {0.128871,0.111478} lies outside the range of data in the
    interpolating function. Extrapolation will be used."
*)

Is there a way to get the actual domain of an InterpolatingFunction?

(and can I see all possible ip["something"], I only found "Domain" by chance)

Answer 1

One way is to create a function that uses Check to see if the message was printed.

InDomainQ[ip_,x_]:=Quiet[Check[ip@@x;True,False,InterpolatingFunction::dmval]]

For the function in the question:

region = RegionPlot[InDomainQ[ip,{x,y}],
    {x,ip[[1,1,1]],ip[[1,1,2]]},
    {y,ip[[1,2,1]],ip[[1,2,2]]},
    MaxRecursion->4,
    PlotPoints->50,
    Epilog->{Red,Point[{x,y}]}]

This seems to be the convex hull of the points:

<<ComputationalGeometry`
hullpts = ip["Grid"][[ConvexHull[ip["Grid"]]]]
hull = Graphics[{
   {Red, Point[{x, y}]},
   {Opacity[0.2], Green, Polygon[hullpts]}
   }];
Show[region, hull]

But how often this is the case I don't know.

Edit: For functions $\mathbb{R}^3 \rightarrow \mathbb{R}$ it actually seems the idea of using the convex hull to determine the domain is more reliable than checking for the message

SeedRandom[1];
dat = RandomReal[1, {9, 4}];
paramdat = Map[{Most[#], Last[#]} &, dat];
ip = Interpolation[paramdat, InterpolationOrder -> 1];

ip[10000,0,0] is clearly outside the domain but does not trigger a message, and checking that the argument is inside the rectangle given by ip["Domain"] would only be partially helpful.

InDomainQ[ip_InterpolatingFunction, x_List] := 
 Quiet[Check[ip @@ x; True, False, InterpolatingFunction::dmval], 
  InterpolatingFunction::dmval]

(* Plot InDomainQ region together with convex hull *)
d = 0.05;
pts = Select[
   Flatten[Table[{x, y, z}, {x, 0, 1, d}, {y, 0, 1, d}, {z, 0, 1, d}],2],
   (InDomainQ[ip, #] &)
   ];

(* ConvexHull3D is in MathWorldPackages.zip
available at http://library.wolfram.com/infocenter/MathSource/4775 *)
hull = ConvexHull3D[ip["Grid"]];

Graphics3D[{
  {Green, Opacity[0.2], EdgeForm[Opacity[0.3]], 
   Cuboid[# - d/2, # + d/2] & /@ pts},
  {Blue, Sphere[ip["Grid"], 0.05]},
  {Opacity[0.3], hull}
  }]

Green is where InDomainQ is True. The blue spheres is the domain part of the interpolated points, and the polygons connecting them form their convex hull.

Since ComplexHull3D returns such nicely oriented Polygons it's easy to construct a new InDomainQ function simply by checking if the signed distance from a given point to each polygon has the same sign:

(* pt's side of the plane passing through the p1,p2,p3 *)
SideOfPlane[{p1_,p2_,p3_},pt_]:=Module[
  {v1=p2-p1,v2=p3-p1,n},
  n=Cross[v1,v2];
  Sign[n.(pt-p1)]
]

GetDomainFunction[ip_]:=Module[
  { hull=ConvexHull3D[ip["Grid"]], planes},
  planes=hull/.Polygon[a_]:>a;
  (* Here is where the orientation from ConvexHull3D comes in
     if the sign from SideOfPlane is same for all the polygons
     it is in the convex hull *)
  Function[{x,y,z},Evaluate[
   Equal@@Function[plane,SideOfPlane[plane,{x,y,z}]]/@planes
  ]]
];

(* InDomainQ[x,y,z], only for this instance of ip *)
InDomainQ=GetDomainFunction[ip];

hull = ConvexHull3D[ip["Grid"]];
Show[{
  Graphics3D[{Opacity[0.4],Green,hull}],
  RegionPlot3D[InDomainQ[x,y,z],{x,0,1},{y,0,1},{z,0,1},
    MaxRecursion->4,
    PlotPoints->30,
    PlotStyle->Directive[Opacity[0.4],Blue]]
}]