Interpolating Function returns unevaluated for some inputs in its domain

Questions
,
1

I ran into this while writing a script to process and plot data output from another tool. The data is expected to be a list of elements, where each element is of the form {x,y,z,intensity}, where intensity is some intensity value. The x-y-z triplets are coordinates which are expected to cover some surface. They are in no particular order. In the example the surface is a torus. What I want to do is to be able to plot the data as a surface, with the colorfunction on the surface representing the intensity value at that point, interpolated from the data.

ListSurfacePlot3D looks perfect for this, with its ColorFunction option used to paint the surface. So for that purpose, I need a function that maps {x,y,z} to an interpolated intensity value. An Interpolating Function seems perfect, with InterpolationOrder->1, since the data is not from a structured grid.

Now to the issue. I construct just such a function, but when I try to use it I find that for some inputs within the range of the function, it is being returned unevaluated. This seemed surprising, since even outside the range an interpolation function generally returns and extrapolated value. (I do admit this usage seems a bit extreme, since the 3D surface is in fact a 2D manifold.)

But what is even more peculiar in that when I apply the function to the exact list of points for which Interpolation was given values, the interpolation function returns some of them unevaluated. I would think it would be defined at least on the exact values on which it is based.

I attach the notebook.
Anyone have any thoughts on this?

Best,
David

2

There are odd things going on here. If you experiment with varying the value of your n you will see some cases where there is no problem (e.g, n=10, n=20). And I was testing where the problem was arising using

Position[applied, InterpolatingFunction[___]]

In the case n=40 another odd thing happens. The above Position expression indicates that there are no unevaluated InterpolationFunciton expressions being returned, but there are error messages generated of the form

InterpolatingFunction::dmval: "Input value {-2.6933817463773266`, 2.6933817463773266`, 0.5877852522924731`} lies outside the range of data in the interpolating function. Extrapolation will be used."

So, give this value a try using your iFunciton by executing

iFunction[-2.6933817463773266, 2.6933817463773266, 0.5877852522924731]

and indeed the same error is returned (but of course the arguments are inside the range of data in the interpolating function).
Now change the first argument of iFunction to 2.693381746377326 (i.e., drop the final 6 in the number).

iFunction[-2.693381746377326, 2.6933817463773266, 0.5877852522924731]

and it returns a value without generating an error message. Now put that 6 back in and exectute the expression once again.

iFunction[-2.6933817463773266, 2.6933817463773266, 0.5877852522924731]

And no no error message is returned even though this is the same expression that generated an error message two executions previously.

Peculiar…
This all smells bugworthy. Perhaps some internal compiled expression is not being properly generated.

3

Thanks, David. I’ll wait and see if any of the Wolfram people comment. Then if it is still a mystery I’ll send it in.

Best,

Also David

4

Hmmm. No further comments I see. I’ll send this in to tech support and see what response I get.

But, this has given rise to some interesting thoughts regarding the construction of interpolating functions from an unstructured grid. As a thought experiment, consider an unstructured grin in R3, where the location of the grid points all lie on a plane which is not parallel to x=0, y=0, or z=0. The collection of points for which values are supplied covers a range of values in all 3 axes. However, the rank of the space spanned by these points is 2, not 3. A rotation and translation could transform these points to lie in, for example, the plane z=0, in a new coordinate system. Then it is clear that there are no values spanning z axis. So although each supplied value is at a point (x,y,0), there are no additional values of z which could be used to establish a linear function representing a neighborhood about this point. The interpolation function could be told that all values interpolated on the x-y plane should be extruded in z, such that (x1,y1,z) = (x1,y1,0) for all z, but does Mathematica do this?

Best,
David

5

I have received a reply from tech support on this:

“Thanks for contacting Wolfram Technical Support and for taking the time to send in the report. This is indeed an issue of our Interpolation function and is caused by the TetGen, which fails to generate a mesh with sufficient quality. I’ve forwarded your message to our development team and I will keep you fully posted on any updates and/or workarounds of this problem.”

So this is indeed an unexpected behavior.

Best,
David

6

In Mathematica version 10 (pre-release) error messages are generated that suggest a problem with the underlying mesh.