diff --git a/docs/make.jl b/docs/make.jl index 8a2835d7..94546761 100644 --- a/docs/make.jl +++ b/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) diff --git a/docs/src/inverting_integrals.md b/docs/src/inverting_integrals.md new file mode 100644 index 00000000..2d10f848 --- /dev/null +++ b/docs/src/inverting_integrals.md @@ -0,0 +1,53 @@ +# 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 +``` + diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index 4ff881d9..d5b69e77 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -90,7 +90,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 diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl index a6219b35..e8ecface 100644 --- a/src/integral_inverses.jl +++ b/src/integral_inverses.jl @@ -1,8 +1,23 @@ abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T} end _integral(A::AbstractIntegralInverseInterpolation, idx, t) = throw(IntegralNotFoundError()) + +""" + 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()) +""" + some stuff +""" struct LinearInterpolationIntInv{uType, tType, itpType, T} <: AbstractIntegralInverseInterpolation{T} u::uType @@ -36,6 +51,9 @@ function _interpolate( u, idx end +""" + some stuff +""" struct ConstantInterpolationIntInv{uType, tType, itpType, T} <: AbstractIntegralInverseInterpolation{T} u::uType