Interpolation over "flat" dimensions

[nilshg]

Not sure whether this is an issue or not (there may well be theoretical reasons for why the current state is correct!), but I'm currently interpolating a multidimensional object linearly within a loop, with grids changing in each iteration. In some iterations, the grid can in certain dimensions become "flat", i.e. just a vector containing repeats of the same value. In this case, Interpolations.jl returns NaNs, which at least to me came as a surprise.

When collapsing the grid down to one point (I inserted a line saying (std(gridvals)>0.1) || (gridpoints=1)), an error is raised that dimensions of length 1 are not (yet) supported - is this something that will happen sometime soon? And might it be worthwhile to raise an error if the grid is flat as well, rather than returning NaNs silently, which might be harder to debug?

[tomasaschan]

Hm. That is definietely not the intended behavior. Can you post a minimal example that reproduces the problem?

[nilshg]

xgrid = collect(linspace(1,10,20))
ygrid = collect(linspace(1,10,20))
func(x,y) = x^2 + y^2
vals = Array(Float64, 20,20)
for x = 1:20, y = 1:20
  vals[x,y] = func(xgrid[x],ygrid[y])
end
using Interpolations
interpolate((xgrid,ygrid),vals, Gridded(Linear()))  # produces expected results

but then:

ygrid = collect(linspace(1,1,20))
for x = 1:20, y = 1:20
  vals[x,y] = func(xgrid[x],ygrid[y])
end
interpolate((xgrid,ygrid), vals, Gridded(Linear()))  # gives only NaNs

This is on Julia 0.4.3 with the latest Interpolations master. Can you reproduce this on your machine?

[tomasaschan]

Yes, I can reproduce.

The simple answer to why this happens, is that since you have 0 distance between your y-points, the interpolation algorithm (which basically does Δz / Δy to interpolate in the y-direction) divides by zero, producing a NaN, which then poisons the rest of the calculations.

Regarding support for length-1 vectors I refer to @timholy, who implemented the Gridded interpolations; I don't know how difficult it would be, but I suppose it might not be too hard, in which case a PR is definitely welcome if you'd like to give it a try. (I can promise to give it a try now that I know it's a feature that someone wants, but I can't give you a timeline - don't hold your breath. And Tim is usually super busy...).

Regarding raising errors, I think the most sensible thing is to throw as soon as two grid coordinates are the same (i.e. if length(unique(grid)) != length(grid)), since otherwise we'd need a lot of special-case logic to make sure that we don't get into trouble if there are multiple data points with different values corresponding to the same coordinate. (If they have the same value, it's not as big a problem - we can just discard one).

As a workaround, maybe you can do something like this?

itp = if length(unique(ygrid)) == 1
    interpolate((xgrid,), vals[:,1])
else
    interpolate((xgrid,ygrid), vals)
end

[nilshg]

That's quite a painful workaround in my case, as it would involve changing all instances in the code where I'm indexing into the extra dimensions of the Lininterp object - I went for the slightly more ecologically unfriedly option of restricting the grid to consist of three points, with the constant value in the middle and some small "fake" variation around it for the interpolation to work.

While doing that I noticed something else I hadn't realised before: with the code above, the interpolant does not equal the underlying function at the knots - may I ask why? I was under the impression that when using the Gridded methods intp[xgrid[i], ygrid[j]] would always equal func(xgrid[i], ygrid[j])?

[tomasaschan]

Yes, itp[xgrid[i], ygrid[j]] == func(xgrid[i], ygrid[j]) should always hold for Gridded, but since (I assume) you're using floating-point stuff, could you check if you get different results with isapprox? Regardless, could you post a complete code sample for which equality doesn't hold?

[nilshg]

Sorry, scratch that - I think there was a mistake in a bit of code I wrote earlier, trying to reproduce this now I find exact equality (even doing == rather than isapprox)

[tomasaschan]

Phew! 😅

[mkitti]

Updating the minimum working example as follows:

using Interpolations

const linspace = LinRange

xgrid = collect(linspace(1,10,20))
ygrid = collect(linspace(1,10,20))
func(x,y) = x^2 + y^2
vals = Array{Float64}(undef, 20,20)
for x = 1:20, y = 1:20
  vals[x,y] = func(xgrid[x],ygrid[y])
end
interpolate((xgrid,ygrid),vals, Gridded(Linear()))  # produces expected results

ygrid = collect(linspace(1,1,20))
for x = 1:20, y = 1:20
  vals[x,y] = func(xgrid[x],ygrid[y])
end
interpolate((xgrid,ygrid), vals, Gridded(Linear()))  # gives only NaNs

From the last line, I currently I get:

┌ Warning: Duplicated knots were deduplicated. Use Interpolations.deduplicate_knots!(knots) explicitly to avoid this warning.
│   k1 =
│    20-element Vector{Float64}:
│     1.0
│     1.0
│     1.0
│     1.0
│     ⋮
│     1.0
│     1.0
│     1.0
└ @ Interpolations ~/.julia/packages/Interpolations/y4lLj/src/gridded/gridded.jl:77
20×20 interpolate((::Vector{Float64},::Vector{Float64}), ::Matrix{Float64}, Gridded(Linear())) with element type Float64:
   2.0        2.0        2.0      …    2.0        2.0        2.0
   3.17175    3.17175    3.17175       3.17175    3.17175    3.17175
   4.79224    4.79224    4.79224       4.79224    4.79224    4.79224
   6.8615     6.8615     6.8615        6.8615     6.8615     6.8615
   9.3795     9.3795     9.3795        9.3795     9.3795     9.3795
  12.3463    12.3463    12.3463   …   12.3463    12.3463    12.3463
  15.7618    15.7618    15.7618       15.7618    15.7618    15.7618
  19.626     19.626     19.626        19.626     19.626     19.626
  23.9391    23.9391    23.9391       23.9391    23.9391    23.9391
  28.7008    28.7008    28.7008       28.7008    28.7008    28.7008
  33.9114    33.9114    33.9114   …   33.9114    33.9114    33.9114
  39.5706    39.5706    39.5706       39.5706    39.5706    39.5706
  45.6787    45.6787    45.6787       45.6787    45.6787    45.6787
  52.2355    52.2355    52.2355       52.2355    52.2355    52.2355
  59.241     59.241     59.241        59.241     59.241     59.241
  66.6953    66.6953    66.6953   …   66.6953    66.6953    66.6953
  74.5983    74.5983    74.5983       74.5983    74.5983    74.5983
  82.9501    82.9501    82.9501       82.9501    82.9501    82.9501
  91.7507    91.7507    91.7507       91.7507    91.7507    91.7507
 101.0      101.0      101.0         101.0      101.0      101.0

[mkitti]

Closing at @jaemolihm's suggestion.