Implement inplace GriddedInterpolation (#496)

* Implement inplace GriddedInterpolation

* Fix io test

* Add inplace gridded example to doc

* Bump version to 0.14.0

Author: jaemolihm

Project.toml

@@ -1,6 +1,6 @@
 name = "Interpolations"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.13.7"
+version = "0.14.0"
 
 [deps]
 AxisAlgorithms = "13072b0f-2c55-5437-9ae7-d433b7a33950"

docs/src/control.md

@@ -123,6 +123,15 @@ Af = [a for a in A if !ismissing(a)]
 itp = interpolate((xf, ), Af, Gridded(Linear()))
 ```
 
+In-place gridded interpolation is also possible:
+```julia
+x = 1:4
+y = view(rand(4), :)
+itp = interpolate!((x,), y, Gridded(Linear()))
+y .= 0
+@show itp(2.5) # 0
+```
+
 ## Parametric splines
 
 Given a set a knots with coordinates `x(t)` and `y(t)`, a parametric spline `S(t) = (x(t),y(t))` parametrized by `t in [0,1]` can be constructed with the following code adapted from a [post](http://julia-programming-language.2336112.n4.nabble.com/Parametric-splines-td37794.html#a37818) by Tomas Lycken:

src/gridded/gridded.jl

@@ -16,7 +16,7 @@ const GridIndex{T} = Union{AbstractVector{T}, Tuple}
 
 struct GriddedInterpolation{T,N,TCoefs,IT<:DimSpec{Gridded},K<:Tuple{Vararg{AbstractVector}}} <: AbstractInterpolation{T,N,IT}
     knots::K
-    coefs::Array{TCoefs,N}
+    coefs::TCoefs
     it::IT
 end
 
@@ -29,12 +29,12 @@ end
 """
     GriddedInterpolation(typeOfWeights::Type{<:Real},
                          knots::NTuple{N, Union{ AbstractVector{T}, Tuple } },
-                         array::AbstractArray{TCoefs,N},
+                         array::AbstractArray{Tel,N},
                          interpolationType::DimSpec{<:Gridded})
 
     Construct a GriddedInterpolation for generic knots from an AbstractArray
 """
-function GriddedInterpolation(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, it::IT) where {N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}
+function GriddedInterpolation(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {N,Tel,TWeights<:Real,IT<:DimSpec{Gridded},pad}
     isconcretetype(IT) || error("The b-spline type must be a leaf type (was $IT)")
 
     check_gridded(it, knots, axes(A))
@@ -44,7 +44,7 @@ function GriddedInterpolation(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::A
     else
         T = typeof(c * first(A))
     end
-    GriddedInterpolation{T,N,TCoefs,IT,typeof(knots)}(knots, A, it)
+    GriddedInterpolation{T,N,typeof(A),IT,typeof(knots)}(knots, A, it)
 end
 
 """
@@ -146,7 +146,7 @@ axes(A::GriddedInterpolation) = axes(A.coefs)
 itpflag(A::GriddedInterpolation) = A.it
 
 function interpolate(::Type{TWeights}, ::Type{TCoefs}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {TWeights,TCoefs,Tel,N,IT<:DimSpec{Gridded}}
-    GriddedInterpolation(TWeights, knots, A, it)
+    GriddedInterpolation(TWeights, knots, copy(A), it)
 end
 """
     itp = interpolate((nodes1, nodes2, ...), A, interpmode)
@@ -169,7 +169,7 @@ end
 interpolate!(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {TWeights,Tel,N,IT<:DimSpec{Gridded}} =
     GriddedInterpolation(TWeights, knots, A, it)
 function interpolate!(knots::NTuple{N,GridIndex}, A::AbstractArray{Tel,N}, it::IT) where {Tel,N,IT<:DimSpec{Gridded}}
-    interpolate!(tweight(A), tcoef(A), knots, A, it)
+    interpolate!(tweight(A), knots, A, it)
 end
 
 lbounds(itp::GriddedInterpolation) = first.(itp.knots)

test/gridded/gridded.jl

@@ -5,75 +5,78 @@ using Interpolations, Test
     front(r::AbstractRange) = range(first(r), step=step(r), length=length(r)-1)
 
     for D in (Constant, Linear)
-        ## 1D
-        a = rand(5)
-        knots = (range(1, stop=length(a), length=length(a)),)
-        itp = @inferred(interpolate(knots, a, Gridded(D())))
-        @inferred(itp(2))
-        @inferred(itp(CartesianIndex(2)))
-        for i = 1:length(a)
-            @test itp(i) ≈ a[i]
-            @test itp(CartesianIndex(i)) ≈ a[i]
-        end
-        @inferred(itp(knots...))
-        @test itp(knots...) ≈ a
-        # compare scalar indexing and vector indexing
-        x = front(knots[1] .+ 0.1)
-        v = itp(x)
-        for i = 1:length(x)
-            @test v[i] ≈ itp(x[i])
-        end
-        x = [2.3,2.2]   # non-increasing order
-        v = itp(x)
-        for i = 1:length(x)
-            @test v[i] ≈ itp(x[i])
-        end
-        # compare against BSpline
-        itpb = @inferred(interpolate(a, BSpline(D())))
-        for x in range(1.1, stop=4.9, length=101)
-            @test itp(x) ≈ itpb(x)
-        end
+        for constructor in (interpolate, interpolate!)
+            isinplace = constructor == interpolate!
+            ## 1D
+            a = rand(5)
+            knots = (range(1, stop=length(a), length=length(a)),)
+            itp = @inferred(constructor(knots, a, Gridded(D())))
+            @inferred(itp(2))
+            @inferred(itp(CartesianIndex(2)))
+            for i = 1:length(a)
+                @test itp(i) ≈ a[i]
+                @test itp(CartesianIndex(i)) ≈ a[i]
+            end
+            @inferred(itp(knots...))
+            @test itp(knots...) ≈ a
+            # compare scalar indexing and vector indexing
+            x = front(knots[1] .+ 0.1)
+            v = itp(x)
+            for i = 1:length(x)
+                @test v[i] ≈ itp(x[i])
+            end
+            x = [2.3,2.2]   # non-increasing order
+            v = itp(x)
+            for i = 1:length(x)
+                @test v[i] ≈ itp(x[i])
+            end
+            # compare against BSpline
+            itpb = @inferred(constructor(a, BSpline(D())))
+            for x in range(1.1, stop=4.9, length=101)
+                @test itp(x) ≈ itpb(x)
+            end
 
-        knots = (range(0.0, stop=1.0, length=length(a)),)
-        itp = @inferred(interpolate(knots, a, Gridded(D())))
-        @test itp([0.1, 0.2, 0.3]) == [itp(0.1), itp(0.2), itp(0.3)]
+            knots = (range(0.0, stop=1.0, length=length(a)),)
+            itp = @inferred(constructor(knots, a, Gridded(D())))
+            @test itp([0.1, 0.2, 0.3]) == [itp(0.1), itp(0.2), itp(0.3)]
 
-        ## 2D
-        A = rand(6,5)
-        knots = (range(1, stop=size(A,1), length=size(A,1)), range(1, stop=size(A,2), length=size(A,2)))
-        itp = @inferred(interpolate(knots, A, Gridded(D())))
-        @test parent(itp) === A
-        @inferred(itp(2, 2))
-        @inferred(itp(CartesianIndex((2,2))))
-        for j = 2:size(A,2)-1, i = 2:size(A,1)-1
-            @test itp(i,j) ≈ A[i,j]
-            @test itp(CartesianIndex((i,j))) ≈ A[i,j]
-        end
-        @test itp(knots...) ≈ A
-        @inferred(itp(knots...))
-        # compare scalar indexing and vector indexing
-        x, y = front(knots[1] .+ 0.1), front(knots[2] .+ 0.6)
-        v = itp(x,y)
-        for j = 1:length(y), i = 1:length(x)
-            @test v[i,j] ≈ itp(x[i],y[j])
-        end
-        # check the fallback vector indexing
-        x = [2.3,2.2]   # non-increasing order
-        y = [3.5,2.8]
-        v = itp(x,y)
-        for j = 1:length(y), i = 1:length(x)
-            @test v[i,j] ≈ itp(x[i],y[j])
-        end
-        # compare against BSpline
-        itpb = @inferred(interpolate(A, BSpline(D())))
-        for x in range(1.1, stop=5.9, length=101), y in range(1.1, stop=4.9, length=101)
-            @test itp(x,y) ≈ itpb(x,y)
-        end
+            ## 2D
+            A = rand(6,5)
+            knots = (range(1, stop=size(A,1), length=size(A,1)), range(1, stop=size(A,2), length=size(A,2)))
+            itp = @inferred(constructor(knots, A, Gridded(D())))
+            isinplace && @test parent(itp) === A
+            @inferred(itp(2, 2))
+            @inferred(itp(CartesianIndex((2,2))))
+            for j = 2:size(A,2)-1, i = 2:size(A,1)-1
+                @test itp(i,j) ≈ A[i,j]
+                @test itp(CartesianIndex((i,j))) ≈ A[i,j]
+            end
+            @test itp(knots...) ≈ A
+            @inferred(itp(knots...))
+            # compare scalar indexing and vector indexing
+            x, y = front(knots[1] .+ 0.1), front(knots[2] .+ 0.6)
+            v = itp(x,y)
+            for j = 1:length(y), i = 1:length(x)
+                @test v[i,j] ≈ itp(x[i],y[j])
+            end
+            # check the fallback vector indexing
+            x = [2.3,2.2]   # non-increasing order
+            y = [3.5,2.8]
+            v = itp(x,y)
+            for j = 1:length(y), i = 1:length(x)
+                @test v[i,j] ≈ itp(x[i],y[j])
+            end
+            # compare against BSpline
+            itpb = @inferred(constructor(A, BSpline(D())))
+            for x in range(1.1, stop=5.9, length=101), y in range(1.1, stop=4.9, length=101)
+                @test itp(x,y) ≈ itpb(x,y)
+            end
 
-        A = rand(8,20)
-        knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
-        itp = interpolate(knots, A, Gridded(D()))
-        @test itp(4,1.2) ≈ A[2,6]
+            A = rand(8,20)
+            knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
+            itp = constructor(knots, A, Gridded(D()))
+            @test itp(4,1.2) ≈ A[2,6]
+        end
     end
 
     # issue #248
@@ -101,4 +104,13 @@ using Interpolations, Test
     successive_knots_warning = "Successive repeated knots detected. Consider using `move_knots` keyword to Interpolations.deduplicate_knots!"
     @test_logs (:warn, successive_knots_warning) Interpolations.deduplicate_knots!(duplicated_knots)
     @test allunique(duplicated_knots)
+
+    # inplace gridded interpolation, issue #495
+    knots = ([0., 1.],)
+    y = view([1., 2.], :)
+    f1 = interpolate(knots, y, Gridded(Linear()))
+    f2 = interpolate!(knots, y, Gridded(Linear()))
+    y .= 0
+    @test f1(0.5) ≈ 1.5
+    @test f2(0.5) ≈ 0
 end

test/io.jl

@@ -51,7 +51,7 @@ using Test
               summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, (Gridded(Linear()), Gridded(Constant{Nearest}()))) with element type Float64"
 
         # issue #260
-        A = (1:4)/4
+        A = collect((1:4)/4)
         itp = interpolate((range(0.0, stop=0.3, length=4),), A, Gridded(Linear()))
         io = IOBuffer()
         show(io, MIME("text/plain"), itp)