# Interpolate.jl a vector of `Point`

Suppose a have a custom type `Point`:

``````struct Point{T}
x::T
y::T
end
``````

I would like to find a way to use Interpolations.jl to interpolate a `Vector{Point{T}}`, but I can’t get my head around it. Say I have this

``````julia> p=[Point(1,1),Point(2,2),Point(3,3)]
3-element Array{Point{Int64},1}:
Point{Int64}(1, 1)
Point{Int64}(2, 2)
Point{Int64}(3, 3)
``````

Is there an alternative to manually converting to `Array`

``````using Interpolations
x = [i.x for i in p]
y = [i.y for i in p]
interpolate((x,),y,Gridded(Linear()))
``````

?

I tried to cast this as a StaticArray

``````julia> rp = reinterpret(SVector{2,Int},p,(3,))
3-element Array{StaticArrays.SArray{Tuple{2},Int64,1,2},1}:
[1, 1]
[2, 2]
[3, 3]
``````

but of course that returns both `x` and `y` as interpolated values, which is not what I want:

``````julia> i = interpolate(rp,BSpline(Linear()),OnGrid())
3-element interpolate(::Array{StaticArrays.SArray{Tuple{2},Int64,1,2},1}, BSpline(Linear()), OnGrid()) with element type StaticArrays.SArray{Tuple{2},Float64,1,2}:
[1, 1]
[2, 2]
[3, 3]

julia> i[1.5]
2-element StaticArrays.SArray{Tuple{2},Float64,1,2}:
1.5
1.5
``````

@sglyon maybe?

I don’t know that there is a way to tell Interpolations one struct field should be the “grid” and the other struct field should be interpolated.

So, I think, for now, you are stuck with the manual version you posted.

Perhaps @tim.holy has a different alternative that would work with the forthcoming re-write of Interpolations.jl

You’ll have to create a special type that performs arithmetic operations just on `y`. Check out the multivalued test in Interpolations’ test suite for info on exactly which operations you need to define.

2 Likes

I knew there must be a way! Awesome!

When you say

you mean to modify the arithmetic methods from this

``````struct Point{T}
x::T
y::T
end

(+)(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y)
(-)(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y)
(*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) = Point(p.x/n, p.y/n)
``````

to that?

``````(+)(p1::Point, p2::Point) = p1.y+p2.y
(-)(p1::Point, p2::Point) = p1.y-p2.y
(*)(n::Number, p::Point) = n*p.y
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) =  p.y/n
``````

Doing so yields an error (sample below) in all interpolations but `Linear()`. For me `Linear()` is ok for now, but I just want to check I’m doing this correctly. Notice that I find that interpolating a `Vector{Point{T}}` takes three times as long it takes to interpolate a `Vector{T}`, created from all `y` fields of `Vector{Point{T}}`. (that result is net of this conversion). Am I missing something or is this to be expected?

``````# very similar to https://github.com/JuliaMath/Interpolations.jl/blob/master/test/b-splines/multivalued.jl

module Points
using Interpolations
using BenchmarkTools

import Base: +, -, *, /

struct Point{T}
x::T
y::T
end

# (+)(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y)
# (-)(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y)
# (*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
# (*)(p::Point, n::Number) = n*p
# (/)(p::Point, n::Number) = Point(p.x/n, p.y/n)
(+)(p1::Point, p2::Point) = p1.y+p2.y
(-)(p1::Point, p2::Point) = p1.y-p2.y
(*)(n::Number, p::Point) = n*p.y
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) =  p.y/n
Base.zero(::Type{Point{T}}) where {T} = Point(zero(T),zero(T))
Base.promote_rule(::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{T1}, ::Type{Point{T2}}) where {T1<:Number,T2} = Point{promote_type(T1,T2)}

n = 1_000_000
m = 10_000
# cast data as vector of Point
data = reinterpret(Point{Float64},hcat(collect(1:n),rand(n))',(n,))
# data to test interpolation on
rdata = rand(m)
# convert vec of points to separate x and y vectors
x = [i.x for i in data]
y = [i.y for i in data]

# my favourite version
function itp_points(d,rd)
itp = scale(interpolate(d, BSpline(Linear()),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
out = similar(rd)
for i in 1:length(rd)
out[i] = itp[rd[i]]
end
out
end

# only interpolate the y vector
function itp_x_y(x,y,rd)
itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
out = similar(rd)
for i in 1:length(rd)
out[i] = itp[rd[i]]
end
out
end

# need to first convert Point to vector!
# an in fact that conversion would have to happen FOR EACH interpolation
# in the below loop... so much more costly in the end!
function itp_pointsxy(d,rd)
x = [i.x for i in data]
y = [i.y for i in data]
itp = scale(interpolate(y, BSpline(Linear()),OnGrid()),linspace(x[1],x[end],length(x)))
out = similar(rd)
for i in 1:length(rd)
out[i] = itp[rd[i]]
end
out
end

# compile
itp_points(data[1:10],rdata[1:5])
itp_x_y(x[1:10],y[1:10],rdata[1:5])
itp_pointsxy(data[1:10],rdata[1:5])

# # BenchmarkTools
t1 = @benchmark itp_points(\$data,\$rdata)
t2 = @benchmark itp_x_y(\$x,\$y,\$rdata)
t3 = @benchmark itp_pointsxy(\$data,\$rdata)

end

julia v0.6.2> Points.t1
BenchmarkTools.Trial:
memory estimate:  15.34 MiB
allocs estimate:  8
--------------
minimum time:     3.784 ms (0.00% GC)
median time:      6.310 ms (14.08% GC)
mean time:        6.688 ms (20.41% GC)
maximum time:     23.664 ms (0.00% GC)
--------------
samples:          746
evals/sample:     1

julia v0.6.2> Points.t2
BenchmarkTools.Trial:
memory estimate:  7.71 MiB
allocs estimate:  8
--------------
minimum time:     1.184 ms (0.00% GC)
median time:      2.036 ms (0.00% GC)
mean time:        2.413 ms (23.11% GC)
maximum time:     7.950 ms (0.00% GC)
--------------
samples:          2064
evals/sample:     1

julia v0.6.2> Points.t3
BenchmarkTools.Trial:
memory estimate:  23.42 MiB
allocs estimate:  30018
--------------
minimum time:     10.395 ms (3.52% GC)
median time:      12.448 ms (21.00% GC)
mean time:        14.273 ms (16.24% GC)
maximum time:     42.232 ms (11.52% GC)
--------------
samples:          350
evals/sample:     1

# checking output
julia v0.6.2> Points.itp_points(Points.data[1:10],Points.rdata[1:5])
5-element Array{Float64,1}:
1.2433
0.97767
1.14297
1.32135
1.18468

WARNING: replacing module Points
ERROR: LoadError: MethodError: Cannot `convert` an object of type Float64 to an object of type Points.Point{Float64}
This may have arisen from a call to the constructor Points.Point{Float64}(...),
since type constructors fall back to convert methods.
Stacktrace:
[1] setindex! at ./multidimensional.jl:300 [inlined]
[2] _A_ldiv_B_md!(::Array{Points.Point{Float64},1}, ::Base.LinAlg.LU{Float64,Tridiagonal{Float64}}, ::Array{Points.Point{Float64},1}, ::CartesianRange{CartesianIndex{0}}, ::CartesianRange{CartesianIndex{0}}, ::Array{Points.Point{Float64},1}) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/../filter1d.jl:30
[3] _A_ldiv_B_md!(::Array{Points.Point{Float64},1}, ::WoodburyMatrices.Woodbury{Float64,Base.LinAlg.LU{Float64,Tridiagonal{Float64}},SparseMatrixCSC{Float64,Int64},SparseMatrixCSC{Float64,Int64},Array{Float64,2},Array{Float64,2}}, ::Array{Points.Point{Float64},1}, ::CartesianRange{CartesianIndex{0}}, ::CartesianRange{CartesianIndex{0}}, ::Array{Points.Point{Float64},1}) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/../filter1d.jl:76
[4] A_ldiv_B_md!(::Array{Points.Point{Float64},1}, ::WoodburyMatrices.Woodbury{Float64,Base.LinAlg.LU{Float64,Tridiagonal{Float64}},SparseMatrixCSC{Float64,Int64},SparseMatrixCSC{Float64,Int64},Array{Float64,2},Array{Float64,2}}, ::Array{Points.Point{Float64},1}, ::Int64, ::Array{Points.Point{Float64},1}) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/../filter1d.jl:15
[5] prefilter!(::Type{Float64}, ::Array{Points.Point{Float64},1}, ::Type{Interpolations.BSpline{Interpolations.Quadratic{Interpolations.Flat}}}, ::Type{Interpolations.OnGrid}) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/prefiltering.jl:58
[6] prefilter(::Type{Float64}, ::Type{Points.Point{Float64}}, ::Array{Points.Point{Float64},1}, ::Type{Interpolations.BSpline{Interpolations.Quadratic{Interpolations.Flat}}}, ::Type{Interpolations.OnGrid}) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/prefiltering.jl:40
[7] interpolate(::Type{Float64}, ::Type{Points.Point{Float64}}, ::Array{Points.Point{Float64},1}, ::Interpolations.BSpline{Interpolations.Quadratic{Interpolations.Flat}}, ::Interpolations.OnGrid) at /Users/florian.oswald/.julia/v0.6/Interpolations/src/b-splines/b-splines.jl:73
[8] itp_points(::Array{Points.Point{Float64},1}, ::Array{Float64,1}) at /Users/florian.oswald/.julia/v0.6/DCEGM/src/int2.jl:41
You still need to return a `Point` object, so it’s something like
``````function (+)(p1::Point, p2::Point)