# Test if argument is inside domain of InterpolatingFunction

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)

-
Take a look at the package << DifferentialEquationsInterpolatingFunctionAnatomy . It's documented in the advanced NDSolve docs, and if you look at the package source, you'll find all these ip["something"] things. The package is simply a wrapper over these. – Szabolcs Dec 8 '12 at 0:35
@Szabolcs Thanks, that was an easy read :) – ssch Dec 8 '12 at 0:44
InterpolatingFunction itself did not seem very well documented. That documentation wasn't linked either. So unless I were using NDSolve or looking here, that would have been hard to find. – ThomasH Jun 6 '14 at 16:33

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

-