Smooth equivalent of ConstantInterpolation in terms of integral

[SouthEndMusic]

Is your feature request related to a problem? Please describe.

A bit of context: In my ODE problem for hydrological modelling, I want to model certain forcing events like precipitation. The precipitation is given as a time series, where constant interpolation ('forward fill') is assumed so that the exact total volume is known. Constant interpolation however does not play nice with ODE solvers because of its lack of smoothness. Therefore I want to have an interpolation type which:

• Is at least $C^1$ smooth;
• Has for each interval between data points the same integral as forward fill would have.

An additional requirement I have is that it is easy to extend this interpolation with new incoming data.

Describe the solution you’d like

The naive solution of assuming constant interpolation on the first interval and then quadratic on the others to satisfy the above conditions is not stable:

using DataInterpolations
using CairoMakie
using Random
using FindFirstFunctions: searchsortedfirst
using QuadGK

fig = Figure()

ax_itp = Axis(fig[1,1]; title = "Interpolation")
ax_int = Axis(fig[2,1]; title = "Integral")

N = 6
Random.seed!(5)
u = rand(N)
t = cumsum(rand(N))
t_eval = range(first(t), last(t), length = 250)

A = ConstantInterpolation(u, t; cache_parameters = true)
lines!(ax_itp, t_eval, A.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A), t_eval); label = "ConstantInterpolation")

a = zeros(N - 1)
b = zeros(N - 1)
c = zeros(N - 1)

# Assume constant interpolation on first interval
c[1] = u[1]

for i in 1:(N-2)
    Δt = t[i+1] - t[i]

    # value at t[i+1]
    u_ = @evalpoly Δt c[i] b[i] a[i]

    # derivative at t[i+1]
    du_ = @evalpoly Δt b[i] 2a[i]

    # integral over (t[i+1], t[i+2])
    Δt = t[i+2] - t[i+1]  
    I = Δt * u[i+1]

    # Next coefficients
    c[i+1] = u_
    b[i+1] = du_
    a[i+1] = 3 * (u[i+1] - u_ - (du_ * Δt) / 2) / Δt^2
end

function B(t_)
    i = clamp(searchsortedfirst(t, t_)-1, 1, N-1)
    @evalpoly (t_ - t[i]) c[i] b[i] a[i]
end

lines!(ax_itp, t_eval, B.(t_eval); label = "Naive smooth interpolation")
lines!(ax_int, t_eval, [quadgk(B, first(t), t_)[1] for t_ in t_eval]; label = "Naive smooth interpolation")

axislegend(ax_itp)
axislegend(ax_int)

fig

[SouthEndMusic]

Oh and of course having periods of 'negative precipitation' to compensate does not make sense. It looks like what I want is that the integral of my interpolation is a PCHIPInterpolation of the integrated values (although that is not $C^1$ smooth and can still overshoot):

using DataInterpolations
using CairoMakie
using Random
using FindFirstFunctions: searchsortedfirst

fig = Figure()

ax_itp = Axis(fig[1,1]; title = "Interpolation")
ax_int = Axis(fig[2,1]; title = "Integral")

N = 6
Random.seed!(5)
u = rand(N)
t = cumsum(rand(N))
t_eval = range(first(t), last(t), length = 250)

A1 = ConstantInterpolation(u, t; cache_parameters = true)

lines!(ax_itp, t_eval, A1.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A1), t_eval); label = "ConstantInterpolation")

u2 = copy(A1.I)
pushfirst!(u2, 0)
A2 = PCHIPInterpolation(u2, t)

lines!(ax_itp, t_eval, DataInterpolations.derivative.(Ref(A2), t_eval); label = "Naive smooth interpolation")
lines!(ax_int, t_eval, A2.(t_eval); label = "Naive smooth interpolation")

axislegend(ax_itp)
axislegend(ax_int)

fig

[SouthEndMusic]

I just realized that as posed above this problem is unsolvable. If the interpolation is nonzero at some interval and you have no precipitation in the interval afterwards, the interpolation cannot be smooth,. have an integral of $0$ in the next interval and be non-negative all at once.

[SouthEndMusic]

For now I'm letting go of the condition of being able to easily add new data to the interpolation, as that poses quite some problems. My new idea is to replace each block of precipitation by a hill that has support over the neighboring intervals, using degree 2 B-spline basis functions:

using DataInterpolations
using CairoMakie
using Random
using QuadGK

N = 25
Random.seed!(8)
u = rand(N)
t = collect(1:N)
t_eval = range(first(t), last(t), length = 250)
A1 = ConstantInterpolation(u, t; cache_parameters = true)

k = t[3:end]
pushfirst!(k, t[1])
pushfirst!(k, first(t))
pushfirst!(k, first(t))
push!(k, last(t))
push!(k, last(t))
sc = zeros(N)
out = zeros(length(sc), length(t_eval))

for (i, t_) in enumerate(t_eval)
    sc .= 0
    DataInterpolations.spline_coefficients!(sc, 2, k, t_)
    out[:, i] .= sc
end

I = sum(out, dims=2) * (t_eval[2] - t_eval[1])
c = zeros(N)
c[1:(N-1)] .= diff(t) .* u[1:(N-1)] ./ I[1:(N-1)]

f = Figure()
ax_itp = Axis(f[1,1]; title = "Interpolation")
ax_int = Axis(f[2,1]; title = "Integral")

lines!(ax_itp, t_eval, A1.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A1), t_eval); label = "ConstantInterpolation")

function A2(t_)
    DataInterpolations.spline_coefficients!(sc, 2, k, t_)
    return sum(sc .* c)
end

lines!(ax_itp, t_eval, A2.(t_eval); label = "QuadraticSpline")
lines!(ax_int, t_eval, [quadgk(A2, first(t), t_)[1] for t_ in t_eval]; label = "QuadraticSpline")

axislegend(ax_itp)
axislegend(ax_int)

f

It must be noted that this only seems to work well with equally spaced data.

[SouthEndMusic]

Better fit and handling of boundaries:

using DataInterpolations
using CairoMakie
using Random
using QuadGK

N = 25
Random.seed!(8)
u = rand(N)
t = collect(1:N)

function refine(u,t)
    N = 2*length(t)-1
    t_new = zeros(N)
    u_new = zeros(N)
    Δt = diff(t)
    t_new[1:2:end] .= t
    t_new[2:2:end] .= t[1:end-1] + Δt/2
    u_new[1:2:end] .= u
    u_new[2:2:end] .= u[1:end-1]
    N = 
    u_new, t_new, N
end

u, t, N = refine(u, t)
t_eval = range(first(t), last(t), length = 500)
A1 = ConstantInterpolation(u, t; cache_parameters = true)

k = copy(t)
pushfirst!(k, first(t))
pushfirst!(k, first(t))
push!(k, last(t))
push!(k, last(t))
sc = zeros(N+1)
out = zeros(N-1, length(t_eval))

for (i, t_) in enumerate(t_eval)
    sc .= 0
    DataInterpolations.spline_coefficients!(sc, 2, k, t_)
    out[1, i] += sc[1]
    out[1, i] += sc[2]
    out[2:(N-2), i] .= sc[3:(N-1)]
    out[N-1, i] += sc[N]
    out[N-1, i] += sc[N+1]
end

Δt = t_eval[2] - t_eval[1]
I = sum(out, dims=2) * Δt
I = I[:,1]
c = diff(t) .* u[1:(N-1)] ./ I
c_extended = copy(c)
pushfirst!(c_extended, first(c))
push!(c_extended, last(c))

f = Figure()
ax_itp = Axis(f[1,1]; title = "Interpolation")
ax_int = Axis(f[2,1]; title = "Integral")
ax_b = Axis(f[3,1]; title = "Basis functions")
scatter!(ax_b, t, zero(t))

for i in 1:(N-1)
    lines!(ax_b, t_eval, out[i, :])
end

lines!(ax_itp, t_eval, A1.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A1), t_eval); label = "ConstantInterpolation")

function A2(t_)
    DataInterpolations.spline_coefficients!(sc, 2, k, t_)
    return sum(sc .* c_extended)
end

lines!(ax_itp, t_eval, A2.(t_eval); label = "QuadraticSpline")
lines!(ax_int, t_eval, [quadgk(A2, first(t), t_)[1] for t_ in t_eval]; label = "QuadraticSpline")

f

It actually looks like this curve has a pretty simple sigmoid unit between the middles of consecutive intervals.

[SouthEndMusic]

The sigmoid section can be used with constant sections in between to generalize to non-equally spaced data:

[SouthEndMusic]

🙌

using DataInterpolations
using FindFirstFunctions: searchsortedfirstcorrelated
using CairoMakie
using Random
using QuadGK

N = 25
Random.seed!(8)
u = rand(N)
t = cumsum(rand(N))

t_eval = range(first(t), last(t), length = 500)
A1 = ConstantInterpolation(u, t; cache_parameters = true)

f = Figure()
ax_itp = Axis(f[1,1]; title = "Interpolation")
ax_int = Axis(f[2,1]; title = "Integral")

lines!(ax_itp, t_eval, A1.(t_eval); label = "ConstantInterpolation")
lines!(ax_int, t_eval, DataInterpolations.integral.(Ref(A1), t_eval); label = "ConstantInterpolation")

Δt = diff(t)
d = [min(Δt[i], Δt[i+1])/2 for i in 1:(N-2)]
Δu = diff(u)

function A2(t_)
    idx = clamp(searchsortedfirstcorrelated(t, t_, 1) - 1, 1, length(t) - 1)
    d_lower = (idx == 1) ? 0 : d[idx-1]
    d_upper = (idx == N-1) ? 0 : d[idx]
    return if (t_ - t[idx]) < d_lower
        u[idx] - 0.5 * Δu[idx-1] * (((t_ - t[idx])/d_lower - 1))^2
    elseif (t[idx+1] - t_) < d_upper
        u[idx] + 0.5 * Δu[idx] * ((1 - (t[idx+1] - t_)/d_upper))^2
    else
        u[idx]
    end
end

lines!(ax_itp, t_eval, A2.(t_eval))
lines!(ax_int, t_eval, [quadgk(A2, first(t), t_)[1] for t_ in t_eval])

f

[SouthEndMusic]

Setting a global maximum on $d$:

using DataInterpolations
using FindFirstFunctions: searchsortedfirstcorrelated
using CairoMakie
using Random
using QuadGK

N = 10
Random.seed!(8)
u = rand(N)
t = cumsum(rand(N))

f = Figure()
ax_itp = Axis(f[1,1]; title = "Interpolation")
ax_int = Axis(f[2,1]; title = "Integral")

Δt = diff(t)
Δu = diff(u)

for d_max in 0.00:0.05:0.3
    d = [min(Δt[i], Δt[i+1], 2d_max)/2 for i in 1:(N-2)]

    function A2(t_)
        idx = clamp(searchsortedfirstcorrelated(t, t_, 1) - 1, 1, length(t) - 1)
        d_lower = (idx == 1) ? 0 : d[idx-1]
        d_upper = (idx == N-1) ? 0 : d[idx]
        return if (t_ - t[idx]) < d_lower
            u[idx] - 0.5 * Δu[idx-1] * (((t_ - t[idx])/d_lower - 1))^2
        elseif (t[idx+1] - t_) < d_upper
            u[idx] + 0.5 * Δu[idx] * ((1 - (t[idx+1] - t_)/d_upper))^2
        else
            u[idx]
        end
    end

    label = "d_max = $d_max"

    lines!(ax_itp, t_eval, A2.(t_eval); label)
    lines!(ax_int, t_eval, [quadgk(A2, first(t), t_)[1] for t_ in t_eval]; label)
end

axislegend(ax_itp)
axislegend(ax_int)

f