Merge pull request SciML#287 from SouthEndMusic/invert_integrals

Cleanup of Invert integrals

docs/src/inverting_integrals.md

@@ -1,19 +1,19 @@
 # Inverting integrals
 
-Solving implicit integral problems of the form:
+Solving implicit integral problems of the following form is supported:
 
 ```math
 \begin{equation}
-    \text{find $t$ such that } \int_{t_1}^t f(\tau)\text{d}\tau = V \ge 0
+    \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:
+where $t_1$ is given by `first(A.t)`. This 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.
+This is achieved 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

src/integral_inverses.jl

@@ -6,8 +6,8 @@ abstract type AbstractIntegralInverseInterpolation{T} <: AbstractInterpolation{T
 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
+  - `A.u` must be a number type (on which an ordering is defined)
+  - This is currently only supported for `ConstantInterpolation` and `LinearInterpolation`
 
 ## Arguments
 
@@ -22,7 +22,16 @@ function _derivative(A::AbstractIntegralInverseInterpolation, t::Number, iguess)
 end
 
 """
-    some stuff
+    LinearInterpolationIntInv(u, t, A)
+
+It is the interpolation of the inverse of the integral of a `LinearInterpolation`.
+Can be easily constructed with `invert_integral(A::LinearInterpolation{<:AbstractVector{<:Number}})`
+
+## Arguments
+
+  - `u` : Given by `A.t`
+  - `t` : Given by `A.I` (the cumulative integral of `A`)
+  - `A` : The `LinearInterpolation` object
 """
 struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
        AbstractIntegralInverseInterpolation{T}
@@ -31,9 +40,10 @@ struct LinearInterpolationIntInv{uType, tType, itpType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     itp::itpType
+    safetycopy::Bool
     function LinearInterpolationIntInv(u, t, A)
         new{typeof(u), typeof(t), typeof(A), eltype(u)}(
-            u, t, A.extrapolate, Ref(1), A)
+            u, t, A.extrapolate, Ref(1), A, A.safetycopy)
     end
 end
 
@@ -56,7 +66,16 @@ function _interpolate(
 end
 
 """
-    some stuff
+    ConstantInterpolationIntInv(u, t, A)
+
+It is the interpolation of the inverse of the integral of a `ConstantInterpolation`.
+Can be easily constructed with `invert_integral(A::ConstantInterpolation{<:AbstractVector{<:Number}})`
+
+## Arguments
+
+  - `u` : Given by `A.t`
+  - `t` : Given by `A.I` (the cumulative integral of `A`)
+  - `A` : The `ConstantInterpolation` object
 """
 struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
        AbstractIntegralInverseInterpolation{T}
@@ -65,9 +84,10 @@ struct ConstantInterpolationIntInv{uType, tType, itpType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     itp::itpType
+    safetycopy::Bool
     function ConstantInterpolationIntInv(u, t, A)
         new{typeof(u), typeof(t), typeof(A), eltype(u)}(
-            u, t, A.extrapolate, Ref(1), A
+            u, t, A.extrapolate, Ref(1), A, A.safetycopy
         )
     end
 end

src/online.jl

@@ -1,10 +1,19 @@
 import Base: append!, push!
 
+function add_integral_values!(A)
+    integral_values = cumsum([_integral(A, idx, A.t[idx + 1]) - _integral(A, idx, A.t[idx])
+                              for idx in (length(A.I) - 1):(length(A.t) - 1)])
+    pop!(A.I)
+    integral_values .+= last(A.I)
+    append!(A.I, integral_values)
+end
+
 function push!(A::LinearInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
     push!(A.u.parent, u)
     push!(A.t.parent, t)
     slope = linear_interpolation_parameters(A.u, A.t, length(A.t) - 1)
     push!(A.p.slope, slope)
+    add_integral_values!(A)
     A
 end
 
@@ -15,12 +24,14 @@ function push!(A::QuadraticInterpolation{U, T}, u::eltype(U), t::eltype(T)) wher
     push!(A.p.l₀, l₀)
     push!(A.p.l₁, l₁)
     push!(A.p.l₂, l₂)
+    add_integral_values!(A)
     A
 end
 
 function push!(A::ConstantInterpolation{U, T}, u::eltype(U), t::eltype(T)) where {U, T}
     push!(A.u.parent, u)
     push!(A.t.parent, t)
+    add_integral_values!(A)
     A
 end
 
@@ -34,6 +45,7 @@ function append!(
     slope = linear_interpolation_parameters.(
         Ref(A.u), Ref(A.t), length_old:(length(A.t) - 1))
     append!(A.p.slope, slope)
+    add_integral_values!(A)
     A
 end
 
@@ -43,6 +55,7 @@ function append!(
     u, t = munge_data(u, t, true)
     append!(A.u.parent, u)
     append!(A.t.parent, t)
+    add_integral_values!(A)
     A
 end
 
@@ -59,5 +72,6 @@ function append!(
     append!(A.p.l₀, l₀)
     append!(A.p.l₁, l₁)
     append!(A.p.l₂, l₂)
+    add_integral_values!(A)
     A
 end

test/online_tests.jl

@@ -18,6 +18,7 @@ for method in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolatio
         @test getfield(func1.p, name) == getfield(func2.p, name)
     end
     @test func1(ts) == func2(ts)
+    @test func1.I == func2.I
 
     func1 = method(u1, t1)
     push!(func1, 1.0, 4.0)
@@ -28,4 +29,5 @@ for method in [LinearInterpolation, QuadraticInterpolation, ConstantInterpolatio
         @test getfield(func1.p, name) == getfield(func2.p, name)
     end
     @test func1(ts) == func2(ts)
+    @test func1.I == func2.I
 end