Merge pull request #282 from SouthEndMusic/invert_integrals

Invert interpolation integrals

docs/make.jl

@@ -13,6 +13,6 @@ makedocs(modules = [DataInterpolations],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/DataInterpolations/stable/"),
     pages = ["index.md", "Methods" => "methods.md",
-        "Interface" => "interface.md", "Manual" => "manual.md"])
+        "Interface" => "interface.md", "Manual" => "manual.md", "Inverting Integrals" => "inverting_integrals.md"])
 
 deploydocs(repo = "github.com/SciML/DataInterpolations.jl"; push_preview = true)

docs/src/inverting_integrals.md

@@ -0,0 +1,52 @@
+# Inverting integrals
+
+Solving implicit integral problems of the form:
+
+```math
+\begin{equation}
+    \text{find $t$ such that } \int_{t_1}^t f(\tau)\text{d}\tau = V \ge 0
+\end{equation}
+```
+
+is supported for interpolations $f$ that are strictly positive and of one of these types:
+
+  - `ConstantInterpolation`
+  - `LinearInterpolation`
+
+This is done by creating an 'integral inverse' interpolation object which can efficiently compute $t$ for a given value of $V$, see the example below.
+
+```@example inverting_integrals
+using Random #hide
+Random.seed!(1234) # hide
+using DataInterpolations
+using Plots
+
+# Create LinearInterpolation object from the
+u = sqrt.(1:25) + (2.0 * rand(25) .- 1.0) / 3
+t = cumsum(rand(25))
+A = LinearInterpolation(u, t)
+
+# Create LinearInterpolationIntInv object
+# from the LinearInterpolation object
+A_intinv = DataInterpolations.invert_integral(A)
+
+# Get the t values up to and including the
+# solution to the integral problem
+V = 25.0
+t_ = A_intinv(V)
+ts = t[t .<= t_]
+push!(ts, t_)
+
+# Plot results
+plot(A; label = "Linear Interpolation")
+plot!(ts, A.(ts), fillrange = 0.0, fillalpha = 0.75,
+    fc = :blues, lw = 0, label = "Area of $V")
+```
+
+## Docstrings
+
+```@docs
+DataInterpolations.invert_integral
+ConstantInterpolationIntInv
+LinearInterpolationIntInv
+```

src/DataInterpolations.jl

@@ -18,6 +18,7 @@ include("interpolation_methods.jl")
 include("plot_rec.jl")
 include("derivatives.jl")
 include("integrals.jl")
+include("integral_inverses.jl")
 include("online.jl")
 include("show.jl")
 
@@ -75,6 +76,18 @@ function Base.showerror(io::IO, e::DerivativeNotFoundError)
     print(io, DERIVATIVE_NOT_FOUND_ERROR)
 end
 
+const INTEGRAL_INVERSE_NOT_FOUND_ERROR = "Cannot invert the integral analytically. Please use Numerical methods."
+struct IntegralInverseNotFoundError <: Exception end
+function Base.showerror(io::IO, e::IntegralInverseNotFoundError)
+    print(io, INTEGRAL_INVERSE_NOT_FOUND_ERROR)
+end
+
+const INTEGRAL_NOT_INVERTIBLE_ERROR = "The Interpolation is not positive everywhere so its integral is not invertible."
+struct IntegralNotInvertibleError <: Exception end
+function Base.showerror(io::IO, e::IntegralNotInvertibleError)
+    print(io, INTEGRAL_NOT_INVERTIBLE_ERROR)
+end
+
 const MUST_COPY_ERROR = "A copy must be made of u, t to filter missing data"
 struct MustCopyError <: Exception end
 function Base.showerror(io::IO, e::MustCopyError)
@@ -84,7 +97,7 @@ end
 export LinearInterpolation, QuadraticInterpolation, LagrangeInterpolation,
        AkimaInterpolation, ConstantInterpolation, QuadraticSpline, CubicSpline,
        BSplineInterpolation, BSplineApprox, CubicHermiteSpline,
-       QuinticHermiteSpline
+       QuinticHermiteSpline, LinearInterpolationIntInv, ConstantInterpolationIntInv
 
 # added for RegularizationSmooth, JJS 11/27/21
 ### Regularization data smoothing and interpolation

src/integral_inverses.jl

@@ -0,0 +1,95 @@
+abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end
+
+"""
+    invert_integral(A::AbstractInterpolation)::AbstractIntegralInverseInterpolation
+
+Creates the inverted integral interpolation object from the given interpolation. Conditions:
+
+  - The range of `A` must be strictly positive
+  - There must be an ordering defined on the data type of `A.u`
+  - This is currently only supported for ConstantInterpolation and LinearInterpolation
+
+## Arguments
+
+  - `A`: interpolation object satisfying the above requirements
+"""
+invert_integral(A::AbstractInterpolation) = throw(IntegralInverseNotFoundError())
+
+_integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError())
+
+function _derivative(A::AbstractIntegralInverseInterpolation, t::Number, iguess)
+    inv(A.itp(A(t))), A.idx_prev[]
+end
+
+"""
+    some stuff
+"""
+struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
+    u::uType
+    t::tType
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    itp::itpType
+    function LinearInterpolationIntInv(u, t, A)
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+            u, t, A.extrapolate, Ref(1), A)
+    end
+end
+
+function invertible_integral(A::LinearInterpolation{<:AbstractVector{<:Number}})
+    return all(A.u .> 0)
+end
+
+function invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}})
+    !invertible_integral(A) && throw(IntegralNotInvertibleError())
+    return LinearInterpolationIntInv(A.t, A.I, A)
+end
+
+function _interpolate(
+        A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess)
+    Δt = t - A.t[idx]
+    x = A.itp.u[idx]
+    u = A.u[idx] + 2Δt / (x + sqrt(x^2 + A.itp.p.slope[idx] * 2Δt))
+    u, idx
+end
+
+"""
+    some stuff
+"""
+struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
+       AbstractIntegralInverseInterpolation{T}
+    u::uType
+    t::tType
+    extrapolate::Bool
+    idx_prev::Base.RefValue{Int}
+    itp::itpType
+    function ConstantInterpolationIntInv(u, t, A)
+        new{typeof(u), typeof(t), typeof(A), eltype(u)}(
+            u, t, A.extrapolate, Ref(1), A
+        )
+    end
+end
+
+function invertible_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}})
+    return all(A.u .> 0)
+end
+
+function invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}})
+    !invertible_integral(A) && throw(IntegralNotInvertibleError())
+    return ConstantInterpolationIntInv(A.t, A.I, A)
+end
+
+function _interpolate(
+        A::ConstantInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
+    idx = get_idx(A.t, t, iguess; ub_shift = 0)
+    if A.itp.dir === :left
+        # :left means that value to the left is used for interpolation
+        idx_ = get_idx(A.t, t, idx; lb = 1, ub_shift = 0)
+    else
+        # :right means that value to the right is used for interpolation
+        idx_ = get_idx(A.t, t, idx; side = :first, lb = 1, ub_shift = 0)
+    end
+    A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx
+end

src/integrals.jl

@@ -6,20 +6,22 @@ end
 function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     ((t1 < A.t[1] || t1 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
+    !hasfield(typeof(A), :I) && throw(IntegralNotFoundError())
     # the index less than or equal to t1
-    idx1 = max(1, min(searchsortedlast(A.t, t1), length(A.t) - 1))
+    idx1 = get_idx(A.t, t1, 0)
     # the index less than t2
-    idx2 = max(1, min(searchsortedlast(A.t, t2), length(A.t) - 1))
-    if A.t[idx2] == t2
-        idx2 -= 1
-    end
-    total = zero(eltype(A.u))
-    for idx in idx1:idx2
-        lt1 = idx == idx1 ? t1 : A.t[idx]
-        lt2 = idx == idx2 ? t2 : A.t[idx + 1]
-        total += _integral(A, idx, lt2) - _integral(A, idx, lt1)
+    idx2 = get_idx(A.t, t2, 0; idx_shift = -1, side = :first)
+
+    total = A.I[idx2] - A.I[idx1]
+    return if t1 == t2
+        zero(total)
+    else
+        total += _integral(A, idx1, A.t[idx1])
+        total -= _integral(A, idx1, t1)
+        total += _integral(A, idx2, t2)
+        total -= _integral(A, idx2, A.t[idx2])
+        total
     end
-    total
 end
 
 function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
@@ -29,7 +31,8 @@ function _integral(A::LinearInterpolation{<:AbstractVector{<:Number}},
     Δt * (A.u[idx] + A.p.slope[idx] * Δt / 2)
 end
 
-function _integral(A::ConstantInterpolation{<:AbstractVector}, idx::Number, t::Number)
+function _integral(
+        A::ConstantInterpolation{<:AbstractVector{<:Number}}, idx::Number, t::Number)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
         return A.u[idx] * t