Merge pull request #49 from helgee/akima

Add Akima interpolation method

src/DataInterpolations.jl

@@ -20,7 +20,7 @@ include("interpolation_utils.jl")
 include("interpolation_alg/interpolation_methods.jl")
 include("plot_rec.jl")
 
-export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
+export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation, AkimaInterpolation,
        ZeroSpline, QuadraticSpline, CubicSpline, BSplineInterpolation, BSplineApprox,
        Loess, GPInterpolation, Curvefit
 end # module

src/caches/interpolation_caches.jl

@@ -41,6 +41,43 @@ function LagrangeInterpolation(u,t,n=nothing)
   LagrangeInterpolation{true}(u,t,n)
 end
 
+### Akima Interpolation
+struct AkimaInterpolation{uType,tType,bType,cType,dType,FT,T} <: AbstractInterpolation{FT,T}
+  u::uType
+  t::tType
+  b::bType
+  c::cType
+  d::dType
+  AkimaInterpolation{FT}(u,t,b,c,d) where FT = new{typeof(u),typeof(t),typeof(b),typeof(c),
+                                                     typeof(d),FT,eltype(u)}(u,t,b,c,d)
+end
+
+function AkimaInterpolation(u, t)
+  u, t = munge_data(u, t)
+  n = length(t)
+  dt = diff(t)
+  m = Array{eltype(u)}(undef, n + 3)
+  m[3:end-2] = diff(u) ./ dt
+  m[2] = 2m[3] - m[4]
+  m[1] = 2m[2] - m[3]
+  m[end-1] = 2m[end-2] - m[end-3]
+  m[end] = 2m[end-1] - m[end-2]
+
+  b = 0.5 .* (m[4:end] .+ m[1:end-3])
+  dm = abs.(diff(m))
+  f1 = dm[3:n+2]
+  f2 = dm[1:n]
+  f12 = f1 + f2
+  ind = findall(f12 .> 1e-9 * maximum(f12))
+  b[ind] = (f1[ind] .* m[ind .+ 1] .+
+            f2[ind] .* m[ind .+ 2]) ./ f12[ind]
+  c = (3.0 .* m[3:end-2] .- 2.0 .* b[1:end-1] .- b[2:end]) ./ dt
+  d = (b[1:end-1] .+ b[2:end] .- 2.0 .* m[3:end-2]) ./ dt.^2
+
+  AkimaInterpolation{true}(u,t,b,c,d)
+end
+
+
 ### ZeroSpline Interpolation
 struct ZeroSpline{uType,tType,dirType,FT,T} <: AbstractInterpolation{FT,T}
   u::uType

src/interpolation_alg/interpolation_methods.jl

@@ -85,6 +85,14 @@ function (A::LagrangeInterpolation{<:AbstractMatrix{<:Number}})(t::Number)
   N/D
 end
 
+function (A::AkimaInterpolation{<:AbstractVector{<:Number}})(t::Number)
+  i = searchsortedlast(A.t, t)
+  i == 0 && return A.u[1]
+  i == length(A.t) && return A.u[end]
+  wj = t - A.t[i]
+  @evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]
+end
+
 # ZeroSpline Interpolation
 function (A::ZeroSpline{<:AbstractVector{<:Number}})(t::Number)
   if A.dir === :left

test/interpolation_tests.jl

@@ -58,6 +58,24 @@ A = LagrangeInterpolation(u,t)
 @test A(1.5) ≈ [2.25,2.25]
 @test A(3.5) ≈ [12.25,12.25]
 
+# Akima Interpolation
+u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
+t = collect(0.0:10.0)
+A = AkimaInterpolation(u, t)
+
+@test A(0.0) ≈ 0.0
+@test A(0.5) ≈ 1.375
+@test A(1.0) ≈ 2.0
+@test A(1.5) ≈ 1.5
+@test A(2.5) ≈ 1.953125
+@test A(3.5) ≈ 2.484375
+@test A(4.5) ≈ 4.1363636363636366866103344
+@test A(5.1) ≈ 5.9803623910336236590978842
+@test A(6.5) ≈ 5.5067291516462386624652936
+@test A(7.2) ≈ 5.2031367459745245795943447
+@test A(8.6) ≈ 4.1796554159017080820603951
+@test A(9.9) ≈ 3.4110386597938129327189927
+@test A(10.0) ≈ 3.0
 
 # ZeroSpline Interpolation
 u = [1.0, 2.0, 0.0, 1.0]