Interpolations WeightedIndex trouble

General Usage
,
1

Have a bit of an advanced Interpolations.jl question here… I guess I just don’t understand how to use the weightedindices or its output.

Based on the docs, I assumed that I pass it the point to which I would want to interpolate, and then I can just use the output to interpolate rather than using the interpolation object itself. It seems that isn’t quite how it works?

julia> using Interpolations

julia> y = sort(rand(15))
15-element Array{Float64,1}:
 0.047341383036748175
 0.05128111599355356
 0.10902355783965456
 0.16089341732308426
 0.3088983112692745
 0.5533832682360411
 0.5771335063852496
 0.6641445249519693
 0.6662254205787557
 0.7066663941780806
 0.758685974796028
 0.7757959797477969
 0.8656767947261979
 0.9239102725958703
 0.9757619516714087

julia> x = collect(1:10.0)
10-element Array{Float64,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0

julia> F =  sin.(x.*y')
10×15 Array{Float64,2}:
 0.0473237  0.0512586  0.108808  0.1602    0.304009    0.525569   0.545624   0.616386   0.618023   0.649302    0.687968   0.700285    0.761534   0.797964   0.828129
 0.0945414  0.102383   0.216323  0.316262  0.57924     0.894256   0.914499   0.970739   0.97173    0.987628    0.998573   0.999816    0.987138   0.961874   0.928395
 0.141547   0.153237   0.32127   0.464155  0.79964     0.996011   0.987131   0.91242    0.909847   0.852941    0.761442   0.72718     0.518044   0.361487   0.21267
 0.188236   0.203689   0.422402  0.600058  0.944344    0.800458   0.739992   0.466221   0.458842   0.309747    0.106646   0.0383993  -0.315624  -0.526134  -0.689975
 0.234503   0.253605   0.518519  0.720461  0.999654    0.365971   0.253141  -0.178174  -0.188402  -0.381797   -0.606648  -0.672356   -0.927172  -0.995694  -0.986184
 0.280244   0.302855   0.608478  0.822254  0.960335   -0.177757  -0.315713  -0.746825  -0.755069  -0.890483   -0.987184  -0.998341   -0.886222  -0.674084  -0.41561
 0.325357   0.351308   0.691212  0.902808  0.83011    -0.668426  -0.782295  -0.997992  -0.998809  -0.972683   -0.826231  -0.753003   -0.221593   0.183147   0.520254
 0.369742   0.398838   0.765738  0.960041  0.621304   -0.959571  -0.99546   -0.824902  -0.815378  -0.589027   -0.212075  -0.0767419   0.598982   0.894851   0.998854
 0.413298   0.445319   0.831172  0.992475  0.353685   -0.964286  -0.886156  -0.301137  -0.283226   0.0767367   0.518408   0.643436    0.998023   0.895514   0.599535
 0.455927   0.490629   0.886736  0.999273  0.0525853  -0.681165  -0.489791   0.350645   0.370056   0.705748    0.964535   0.995393    0.694705   0.18461   -0.32673

julia> itp = LinearInterpolation((x,y),F,extrapolation_bc=Line())
10×15 extrapolate(interpolate((::Array{Float64,1},::Array{Float64,1}), ::Array{Float64,2}, Gridded(Linear())), Line()) with element type Float64:
 0.0473237  0.0512586  0.108808  0.1602    0.304009    0.525569   0.545624   0.616386   0.618023   0.649302    0.687968   0.700285    0.761534   0.797964   0.828129
 0.0945414  0.102383   0.216323  0.316262  0.57924     0.894256   0.914499   0.970739   0.97173    0.987628    0.998573   0.999816    0.987138   0.961874   0.928395
 0.141547   0.153237   0.32127   0.464155  0.79964     0.996011   0.987131   0.91242    0.909847   0.852941    0.761442   0.72718     0.518044   0.361487   0.21267
 0.188236   0.203689   0.422402  0.600058  0.944344    0.800458   0.739992   0.466221   0.458842   0.309747    0.106646   0.0383993  -0.315624  -0.526134  -0.689975
 0.234503   0.253605   0.518519  0.720461  0.999654    0.365971   0.253141  -0.178174  -0.188402  -0.381797   -0.606648  -0.672356   -0.927172  -0.995694  -0.986184
 0.280244   0.302855   0.608478  0.822254  0.960335   -0.177757  -0.315713  -0.746825  -0.755069  -0.890483   -0.987184  -0.998341   -0.886222  -0.674084  -0.41561
 0.325357   0.351308   0.691212  0.902808  0.83011    -0.668426  -0.782295  -0.997992  -0.998809  -0.972683   -0.826231  -0.753003   -0.221593   0.183147   0.520254
 0.369742   0.398838   0.765738  0.960041  0.621304   -0.959571  -0.99546   -0.824902  -0.815378  -0.589027   -0.212075  -0.0767419   0.598982   0.894851   0.998854
 0.413298   0.445319   0.831172  0.992475  0.353685   -0.964286  -0.886156  -0.301137  -0.283226   0.0767367   0.518408   0.643436    0.998023   0.895514   0.599535
 0.455927   0.490629   0.886736  0.999273  0.0525853  -0.681165  -0.489791   0.350645   0.370056   0.705748    0.964535   0.995393    0.694705   0.18461   -0.32673

julia> wi = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp)...,(5.0,0.5))
(Interpolations.WeightedAdjIndex{2,Float64}(4, (0.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(1, (1.5, -0.5)))

julia> F[wi...]
0.2249513239238678

julia> itp(5,0.5)
0.5043358836315593
2

WeightedIndexes refer to the integral indices of the underlying array, not the spatial location.

3

Ok I can see that now from the result. Presumably there is an intermediate function that Interpolations.jl uses to go from the actual point (in this case (5,0.5)), down to the indices. Could you point me in that direction? I am looking through where I think it should be and I do not see anything.

https://github.com/JuliaMath/Interpolations.jl/blob/master/src/gridded/indexing.jl

4

Figured it out I think! Not sure if this would always be true or not because my example uses one of the convenience constructors, but this piece in the call to weightedindexes

Interpolations.itpinfo(itp)

should be

Interpolations.itpinfo(itp.itp)

Ouput

julia> wi = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp.itp)...,(5,0.5))
(Interpolations.WeightedAdjIndex{2,Float64}(4, (0.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(5, (0.21834990953368805, 0.781650090466312)))

julia> F[wi...]
0.5043358836315593 # now same as itp(5,0.5)

I should, and this builds on what Tim was saying (and already did say in the Dev Docs), that this specific way only works for Linear interpolation. Otherwise the weights should be used to index the interpolant coefficients.

5

Just a further warning on this. The extrapolation that takes place by doing this is linear but not the same as the Line() boundary condition given by the package.