Interpolate.jl a vector of `Point`

General Usage
1

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
2

@sglyon maybe? :slight_smile:

3

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

4

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
5

I knew there must be a way! Awesome!

6

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(Quadratic(Flat())),OnGrid()),linspace(d[1].x,d[end].x,length(d)))
	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

# traditional version
# 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

# actual cost of traditional:
# 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


# error for Quadratic()

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
 [9] include_from_node1(::String) at ./loading.jl:576
 [10] include(::String) at ./sysimg.jl:14
while loading /Users/florian.oswald/.julia/v0.6/DCEGM/src/int2.jl, in expression starting on line 74
7

You still need to return a Point object, so it’s something like

function (+)(p1::Point, p2::Point)
    @noinline checksame(p1, p2) = p1.x == p2.x || throw(ArgumentError("arithmetic supported only for constant x, got $p1 and $p2"))

    checksame(p1, p2)
    Point(p1.x, p1.y + p2.y)
end