Interpolations of a function that has NaN

roi.holtzman · March 7, 2023, 11:26pm

I am trying to interpolate a 2d function that is defined only on a partial set of my parameters.
Where it is not defined, my 2d array has NaNs. Note that in my example the range of where the function is defined is not a rectangle, but it is a function of x, y.

When I perform the interpolation, I get a function that has all NaNs. Is there a way to overcome this issue?

I am using Interpolations.jl, and my code is like the example in the documentation.
Here is a minimal example which includes NaNs:

A_x1 = 1:.1:10
A_x2 = 1:.5:20
f(x1, x2) = log(x1+x2)
A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
A[1:2, :] .= NaN
A[:, 1:2] .= NaN

itp = interpolate(A, BSpline(Cubic(Line(OnGrid()))))

The result of itp is all NaN.

rafael.guerra · March 8, 2023, 12:12am

Instead of using a matrix with NaNs, one possibility is to take the remaining valid entries as scattered points. See the example below using ScatteredInterpolation.jl.

ScatteredInterpolation.jl code

using ScatteredInterpolation, StaticArrays, Plots

A_x1 = 1:.1:10
A_x2 = 1:.5:20
f(x1, x2) = log(x1+x2)
A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
A[50, :] .= NaN
A[:, 20] .= NaN

CI = findall(!isnan, A)
points = SVector{2,Float64}[]
samples = Float64[]
for ci in CI
    push!(points, SA[A_x1[ci[1]], A_x2[ci[2]]])
    push!(samples, A[ci])
end
points = reduce(hcat, points)

nx1, nx2 = length(A_x1), length(A_x2)
X1, X2 = vec(repeat(A_x1, nx2)), vec(repeat(A_x2', nx1))
gridPoints = [X1 X2]'

itp = ScatteredInterpolation.interpolate(Multiquadratic(), points, samples)
interpolated = evaluate(itp, gridPoints)
gridded = reshape(interpolated, nx1, nx2)

plot(heatmap(A), heatmap(gridded))

The packages GMT.jl and DIVAnd.jl are also very good for this task: see this other post.

roi.holtzman · March 9, 2023, 3:34pm

Thanks for your response!
It has been useful, but I am still not sure how to deal with my problem.

My function f(x, y) is defined on some non-rectangular grid. Namely, there is a line y_{edge}(x) (known only numerically) that gives the boundary of where my function f(x, y) is defined. So when I calculate the values f(x,y) on some 2d grid I put NaN where the function is not defined.

The best scenario for me would be to simply extrapolate the function f(x,y) where it is not defined.

I tried to fill all the NaN values with flat values outside the region where the function is defined and to do the interpolations, but that does not work very well.

If you have other ideas on how to do that, I would love to hear!

JM_Beckers · March 9, 2023, 4:00pm

With DIVAnd you can define a mask for your grid onto which you want to interpolate and the data points can be scattered. So basically you defined a 2D grid and mask it for the regions where you know your function is not defined. Only THEN you ask DIVAnd to interpolate and the interpolation will taken into account the mask. So you already created that mask in your tests and only need to use the (scattered) data points you actually have for the interpolation.

Also is your edge function known at the resolution of the interpolation grid ? Otherwise you probably need to interpolate it first to the right resolution too.

Hope it makes sense, otherwise I will provide a MWE.

JM_Beckers · March 9, 2023, 4:34pm

using DIVAnd
using PyPlot
using Statistics



NX=300
NY=300
len=0.1

ND=NX
# Fine grid
xg=collect(range(0,stop=1,length=NX))
yg=collect(range(0,stop=1,length=NY))


# Boundary on coarse grid
yedge=[0.1,0.2,0.3,0.2,0.2]
xedge=[0.1,0.3,0.5,0.7,0.9]


# Interpolate fine boundary
yedgefine,s=DIVAndrun(trues(NX),(ones(size(xg)) / (xg[2]-xg[1]),),(xg,),(xedge,),yedge.-mean(yedge),len,0.01)
yedgefine=yedgefine.+mean(yedge)


figure()

plot(xg,yedgefine,"-",xedge,yedge,"o")

figure()


# Sample a function

fun(x,y) = 2*(sin.(6x) * cos.(6y)) .+ (0 .- 1) .* x .* y
x = 0.5.+0.25.*randn(ND);
y = 0.5  .+ 0.25 .* randn(ND);
f = fun.(x,y)+0.1*randn(ND)


# create the grid for the interpolation
xi,yi= DIVAnd.ndgrid(xg,yg)


# all points are valid points except those below the edge
mask = trues(size(xi));

for j=1:NY
    for i=1:NX
        if yi[i,j]<yedgefine[i]
            mask[i,j]=false
        end
    end
end
            
# metrics
pm = ones(size(xi)) / (xi[2,1]-xi[1,1]);
pn = ones(size(xi)) / (yi[1,2]-yi[1,1]);

# Analysis
epsilon2=1.0
fi,s=DIVAndrun(mask,(pm,pn),(xi,yi),(x,y),f,(len,len),epsilon2)


# Plot it with the edge function
pcolor(xi,yi,fi,shading="nearest"),colorbar()
plot(xg,yedgefine,"-")

rafael.guerra · March 10, 2023, 7:14pm

@JM_Beckers, could you please verify your code as it is giving errors?

After correcting the espilon2 typo, the subsequent command:

fi, s = DIVAndrun(mask, (pm,pn), (xi,yi), (x,y), f, (len,len), epsilon2)

throws:

ERROR: CanonicalIndexError: setindex! not defined for LazyGrids.GridAV{Float64, 1, 2}     
Stacktrace:
  [1] error_if_canonical_setindex(::IndexCartesian, ::LazyGrids.GridAV{Float64, 1, 2}, ::Int64, ::Int64)

Using Julia 1.8.5 on Windows and DIVAnd v2.7.9.

JM_Beckers · March 10, 2023, 9:25pm

I think there is a clash for function ndgrid which is also defined in Methods · LazyGrids.jl.
Probably you loaded that one ? I changed to DIVAnd.ndgrid to make sure you use the right version. Let me know if it works.
(I also corrected the typo, thanks)

rafael.guerra · March 10, 2023, 10:19pm

@JM_Beckers, thank you, you were spot on.

roi.holtzman · March 10, 2023, 10:47pm

Thanks for this wonderful example!

I tried to implement it in my case, but it did not work for me.
I get the following error:

MethodError: no method matching DIVAndrun(::BitMatrix, ::Tuple{Matrix{Float64}, Matrix{Float64}}, ::Tuple{Matrix{Float64}, Matrix{Float64}}, ::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::typeof(Sλq_2d_grid_interpolate_function), ::Tuple{Float64, Float64}, ::Float64)

For me the function f in your code is a function that is defined using Interpolations.jl.
To have a function and not an interpolations type, I defined
f(x, y) = interpolations_f(x,y).
Still, it does not work.
Do you have any idea of how to make it work?

JM_Beckers · March 11, 2023, 1:29am

No, f is an array containing the scattered data values located in (x,y), where x and y are the arrays containing the two coordinate values of each value of f.

In the MWE they are created by a random set of coordinates and a know function,

x = 0.5.+0.25.*randn(ND);
y = 0.5  .+ 0.25 .* randn(ND);
f = fun.(x,y)+0.1*randn(ND)

but you can provide the arrays you want. Basically you tell DIVAnd in which points you know the function value and which is the value.