Merge pull request #525 from ChrisRackauckas-Claude/makima-via-akima-kwarg

Add modified Akima (makima) as a kwarg on AkimaInterpolation

Author: ChrisRackauckas

docs/make.jl

@@ -17,6 +17,7 @@ makedocs(
     ),
     linkcheck_ignore = [
         "https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf",
+        "https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/",
     ],
     pages = [
         "index.md",

src/interpolation_caches.jl

@@ -220,7 +220,7 @@ function LagrangeInterpolation(
 end
 
 """
-    AkimaInterpolation(u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+    AkimaInterpolation(u, t; modified = false, extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false)
 
 It is a spline interpolation built from cubic polynomials. It forms a continuously differentiable function. For more details, refer: [https://en.wikipedia.org/wiki/Akima_spline](https://en.wikipedia.org/wiki/Akima_spline).
@@ -233,6 +233,11 @@ Extrapolation extends the last cubic polynomial on each side.
 
 ## Keyword Arguments
 
+  - `modified`: if `true`, use the modified Akima (makima) formula for the slopes at the knots,
+    which adds an extra term `|m_{i+1} + m_i| / 2` to each weight. Tends to reduce overshoot
+    and oscillation on data with flat regions or repeated values. See
+    [https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/](https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/).
+    Defaults to `false`.
   - `extrapolation`: The extrapolation type applied left and right of the data. Possible options
     are `ExtrapolationType.None` (default), `ExtrapolationType.Constant`, `ExtrapolationType.Linear`
     `ExtrapolationType.Extension`, `ExtrapolationType.Periodic` and `ExtrapolationType.Reflective`.
@@ -284,7 +289,8 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
 end
 
 function AkimaInterpolation(
-        u, t; extrapolation::ExtrapolationType.T = ExtrapolationType.None,
+        u, t; modified::Bool = false,
+        extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, cache_parameters = false, assume_linear_t = 1.0e-2
     )
@@ -303,16 +309,30 @@ function AkimaInterpolation(
     m[end - 1] = 2m[end - 2] - m[end - 3]
     m[end] = 2m[end - 1] - m[end - 2]
 
-    b = (m[4:end] .+ m[1:(end - 3)]) ./ 2
-    dm = abs.(diff(m))
-    f1 = dm[3:(n + 2)]
-    f2 = dm[1:n]
-    f12 = f1 + f2
-    ind = findall(f12 .> 1.0e-9 * maximum(f12))
-    b[ind] = (
-        f1[ind] .* m[ind .+ 1] .+
-            f2[ind] .* m[ind .+ 2]
-    ) ./ f12[ind]
+    if modified
+        # Modified Akima (makima): adds |m_{i+1} + m_i| / 2 to each weight, which
+        # reduces overshoot on flat regions. The simple-average fallback still
+        # guards the case where all four neighboring slopes vanish.
+        w1 = abs.(m[4:end] .- m[3:(end - 1)]) .+
+            abs.(m[4:end] .+ m[3:(end - 1)]) ./ 2
+        w2 = abs.(m[2:(end - 2)] .- m[1:(end - 3)]) .+
+            abs.(m[2:(end - 2)] .+ m[1:(end - 3)]) ./ 2
+        w12 = w1 .+ w2
+        b = (m[2:(end - 2)] .+ m[3:(end - 1)]) ./ 2
+        ind = findall(w12 .> 1.0e-9 * maximum(w12))
+        b[ind] = (w1[ind] .* m[ind .+ 1] .+ w2[ind] .* m[ind .+ 2]) ./ w12[ind]
+    else
+        b = (m[4:end] .+ m[1:(end - 3)]) ./ 2
+        dm = abs.(diff(m))
+        f1 = dm[3:(n + 2)]
+        f2 = dm[1:n]
+        f12 = f1 + f2
+        ind = findall(f12 .> 1.0e-9 * maximum(f12))
+        b[ind] = (
+            f1[ind] .* m[ind .+ 1] .+
+                f2[ind] .* m[ind .+ 2]
+        ) ./ f12[ind]
+    end
     c = (3 .* m[3:(end - 2)] .- 2 .* b[1:(end - 1)] .- b[2:end]) ./ dt
     d = (b[1:(end - 1)] .+ b[2:end] .- 2 .* m[3:(end - 2)]) ./ dt .^ 2
 

test/interpolation_tests.jl

@@ -511,6 +511,60 @@ end
     A = @inferred(AkimaInterpolation(u, t))
     @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
     @test_throws DataInterpolations.RightExtrapolationError A(11.0)
+
+    # Modified Akima (makima) — same struct, different b computation
+    @testset "Modified Akima (makima)" begin
+        A_makima = @inferred(AkimaInterpolation(u, t; modified = true))
+        # Interpolates the data points exactly
+        for i in eachindex(t)
+            @test A_makima(t[i]) ≈ u[i]
+        end
+        # Differs from standard Akima on this data set
+        A_std = AkimaInterpolation(u, t)
+        @test any(
+            !isapprox(A_makima(ti), A_std(ti); atol = 1.0e-12)
+                for ti in 0.5:1.0:9.5
+        )
+
+        # Compare against a direct vectorized implementation of the makima formula
+        n = length(t)
+        dt_vec = diff(t)
+        m = Array{eltype(u)}(undef, n + 3)
+        m[3:(end - 2)] = diff(u) ./ dt_vec
+        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]
+        w1 = abs.(m[4:end] .- m[3:(end - 1)]) .+
+            abs.(m[4:end] .+ m[3:(end - 1)]) ./ 2
+        w2 = abs.(m[2:(end - 2)] .- m[1:(end - 3)]) .+
+            abs.(m[2:(end - 2)] .+ m[1:(end - 3)]) ./ 2
+        b_ref = (w1 .* m[2:(end - 2)] .+ w2 .* m[3:(end - 1)]) ./ (w1 .+ w2)
+        @test A_makima.b ≈ b_ref
+
+        # Makima avoids the w1 + w2 == 0 division-by-zero edge case
+        # that the original Akima formula explicitly works around: a
+        # constant-then-constant signal produces zero forward and backward
+        # slope differences at the transition.
+        u_flat = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
+        t_flat = collect(0.0:5.0)
+        A_flat = @inferred(AkimaInterpolation(u_flat, t_flat; modified = true))
+        @test all(isfinite, A_flat.b)
+        @test all(isfinite, A_flat.c)
+        @test all(isfinite, A_flat.d)
+        for i in eachindex(t_flat)
+            @test A_flat(t_flat[i]) ≈ u_flat[i]
+        end
+
+        # Extrapolation, derivatives and integrals work via the shared struct path
+        A_ext = AkimaInterpolation(
+            u, t; modified = true, extrapolation = ExtrapolationType.Extension
+        )
+        @test isfinite(A_ext(-1.0))
+        @test isfinite(A_ext(11.0))
+        @test isfinite(DataInterpolations.derivative(A_makima, 5.0))
+        @test isfinite(DataInterpolations.integral(A_makima, 0.0, 10.0))
+    end
 end
 
 @testset "ConstantInterpolation" begin