Hi everyone,
I am have run into a problem with interpolating values across large arrays using Julia. I have read these Discourse discussions pertaining to my problem but have not been able to use them to create a solution:
- Scattered Interpolations - good information that exposed me to a lot of packages that do interpolation well. Honestly, this problem is quite similar to mine but I am not trying to related it back to Julia.
- Matrix Interpolations - great info on keeping my types clean.
- Using Interpolations - good theory and how to convert from Python to Julia but not what I am looking for.
- Trouble Making Topographical EEG Interpolation - cool visualization ideas but not what I am trying to accomplish.
Here is the information on what I have been attempting:
Problem Background
I have a grid with dimensions m x n where the grid is rectangular in shape but could be square. It is represented as a matrix constructed like grid = ones(m, n). From there, I have an array of observations that I would like to use as a basis to interpolate values across the grid. These observations can be thought of as observations = ([xs], [ys], [values]) where xs and ys represent lists of x and y points on the grid where each x and y pair have an associated value that should be used to interpolate values surrounding each x and y location for use in a 2D interpolation.
What I run into is that in my specific implementations, I have large arrays that cause memory errors. Further, I am curious if I am fundamentally misunderstanding something about interpolation or how Julia handles interpolations as my numerous efforts have failed. Here are my attempts at finding a solution for a specific implementation of my problem:
Attempts on a Specific Problem
In this case, what I was attempting is to define my grid like grid = ones(500, 500) with 31 observations of interest. My observations look like this:
julia> observations
Tuple{Array{Float64,1},Array{Float64,1},Array{Float32,1}}
([187.0, 313.0, 157.0, 343.0, 132.0, 368.0, 157.0, 343.0, 187.0, 313.0, 85.0, 415.0, 46.0, 454.0, 85.0, 415.0, 250.0, 250.0, 57.0, 443.0, 129.0, 371.0, 164.0, 336.0, 250.0, 250.0, 250.0, 180.0, 320.0, 118.0, 382.0],
[57.0, 57.0, 134.0, 134.0, 250.0, 250.0, 366.0, 366.0, 443.0, 443.0, 129.0, 129.0, 250.0, 250.0, 371.0, 371.0, 132.0, 368.0, 313.0, 313.0, 415.0, 415.0, 403.0, 403.0, 417.0, 454.0, 475.0, 466.0, 466.0, 433.0, 433.0],
Float32[-19.995152, -46.38638, -28.329832, 81.12224, -67.002556, 21.798035, 32.538525, 38.022156, -14.682513, -56.849937, -20.680933, 95.41259, 22.357513, 36.404854, -37.628937, -3.4352758, -25.271269, 36.05853, 26.666668, -29.088951, -7.4670925, -47.769047, -29.178823, -43.958958, 29.638771, -10.448129, 5.71033, -25.612432, 37.268364, -22.657557, -2.0058606])
Attempt 1: Using ScatteredInterpolation.jl
My first thought was to used ScatteredInterpolation.jl since I am basically working with scatter data, at least I thought. Trying to read the sparse documentation provided by the package, my first thought was to convert my observations to an array of length, 250000, with every value that I don’t care about set to one. Further, I modified the xs and ys to contain each x value and y value possible resulting in 250000 ordered pairs.
When I tried to interpolate like this interpolate(MultiQuadratic(), [xs, ys], samples), I consistently got OutOfMemoryError() messages. I couldn’t figure out how to remedy this and moved on to other approaches.
Attempt 2: Using Interpolations.jl
I found Interpolations.jl’s documentation a bit more friendly and informative. Using the docs, my initial thoughts was to do the following:
julia> nodes = (values, )
(Any[22.357513f0, 26.666668f0, -20.680933f0, -37.628937f0, -22.657557f0, -7.4670925f0, -67.002556f0, -28.329832f0, 32.538525f0, -29.178823f0 … -43.958958f0, 81.12224f0, 38.022156f0, 21.798035f0, -47.769047f0, -2.0058606f0, 95.41259f0, -3.4352758f0, -29.088951f0, 36.404854f0],)
julia> grid = ones(500, 500);
julia> using Interpolations
julia> itp = interpolate(nodes, grid, Gridded(Linear()))
ERROR: MethodError: no method matching interpolate(::Tuple{Array{Any,1}}, ::Array{Float64,2}, ::Gridded{Linear})
I consistently struggled with getting the right form it was expecting. When I did get my data into an acceptable form, I wound up with these errors:
julia> itp = interpolate(nodes, grid, Gridded(Linear()))
ERROR: DimensionMismatch("knot vectors must have the same axes as the corresponding dimension of the array")
Which I did not understand as I thought I had dimensions correctly matched. Further, what is a knot? From here, I got frustrated and turned to my friend who is an expert in MATLAB for help.
Attempt 3: Reimplementation of MATLAB’s interp2 Function
Thanks to this handy blog, I came across a Julia implementation of interp2! It built upon Interpolations.jl so I was able to understand what it was doing based on Attempt 2 and reading the docs on MATLAB’s interp2 function. I created this function:
function interp2(X, Y, V, Xq, Yq)
knots = (X,Y)
itp = interpolate(knots, V, Gridded(Linear()))
itp[Xq, Yq]
end
I redefined xs and ys to only be 31 ordered pairs (each x and y were in their own array) and values as the only 31 values I cared about. I then invoked this reimplementation to try to find the interpolated value at 300, 300 and got similar errors from Attempt 1 and 2:
julia> interp2(xs, ys, values, 300, 300)
ERROR: MethodError: no method matching interpolate(::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Gridded{Linear})
Closest candidates are:
interpolate(::Tuple{Vararg{Union{AbstractArray{T,1}, Tuple} where T,N}}, ::AbstractArray{Tel,N}, ::IT) where {Tel, N, IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, Gridded},N} where N}, Gridded}} at /home/cedarprince/.julia/packages/Interpolations/TBuvH/src/gridded/gridded.jl:83
interpolate(::Type{TWeights}, ::Type{TC}, ::Any, ::IT) where {TWeights, TC, IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, BSpline},N} where N}, BSpline}} at /home/cedarprince/.julia/packages/Interpolations/TBuvH/src/b-splines/b-splines.jl:159
interpolate(::AbstractArray{var"#s67",1} where var"#s67"<:Number, ::AbstractArray{TEl,1}, ::TInterpolationType) where {TEl, TInterpolationType<:Interpolations.MonotonicInterpolationType} at /home/cedarprince/.julia/packages/Interpolations/TBuvH/src/monotonic/monotonic.jl:177
...
Stacktrace:
[1] interp2(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Int64) at ./REPL[95]:3
[2] top-level scope at REPL[122]:1
This is so confusing to me because based on the MATLAB docs, I was fully expecting this to work out of all my attempts.
Conclusion
As I have not been able to solve this problem, I have therefore come to the following conclusions:
- I am missing something fundamental about how Julia handles interpolation - both for gridded and non-gridded situations.
- I somewhere have a fundamental misunderstanding about interpolations.
- I discovered bugs that I have been running into consistently.
- Either I lack intuition about how interpolations work or there is something strange about how Julia does it as I come from Python and MATLAB.
Would anyone be able to please help me? I have been struggling all day and would greatly appreciate the help. Also, shoutout to the interpolations folks in the Julia community! I am sure these are great tools - I am just struggling with them!
Take care and have a great rest of the day.
~ TCP 
P.S. If more information is needed, please let me know and I would be happy to include code snippets, data, or what you may need.
However, the package name itself is still not informative enough for me. As someone who is looking for a proper solution for plotting, it is really unintuitive to think at first that the desired function lies underneath the GeoStats.jl package.