Merge remote-tracking branch 'upstream/master' into throw_extrapolati…

…on_error

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/DataInterpolations.jl

@@ -23,7 +23,7 @@ include("online.jl")
 include("show.jl")
 
 (interp::AbstractInterpolation)(t::Number) = _interpolate(interp, t)
-function (interp::AbstractInterpolation)(t::Number, i::Integer) 
+function (interp::AbstractInterpolation)(t::Number, i::Integer)
     interp.idx_prev[] = i
     _interpolate(interp, t)
 end

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/interpolation_caches.jl

@@ -481,6 +481,7 @@ function BSplineInterpolation(
         u, t, d, pVecType, knotVecType; extrapolate = false, safetycopy = true)
     u, t = munge_data(u, t, safetycopy)
     n = length(t)
+    n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), n + d + 1)
@@ -615,6 +616,7 @@ function BSplineApprox(
         u, t, d, h, pVecType, knotVecType; extrapolate = false, safetycopy = true)
     u, t = munge_data(u, t, safetycopy)
     n = length(t)
+    h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), h + d + 1)

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/interpolation_tests.jl

@@ -616,6 +616,12 @@ end
         @test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]
         @test [A(t[1]), A(t[end])] == [u[1], u[end]]
 
+        @test_throws ErrorException("BSplineInterpolation needs at least d + 1, i.e. 4 points.") BSplineInterpolation(
+            u[1:3], t[1:3], 3, :Uniform, :Uniform)
+        @test_throws ErrorException("BSplineInterpolation needs at least d + 1, i.e. 5 points.") BSplineInterpolation(
+            u[1:4], t[1:4], 4, :ArcLen, :Average)
+        @test_nowarn BSplineInterpolation(u[1:3], t[1:3], 2, :Uniform, :Uniform)
+
         # Test extrapolation
         A = BSplineInterpolation(u, t, 2, :ArcLen, :Average; extrapolate = true)
         @test A(-1.0) == u[1]
@@ -634,6 +640,12 @@ end
         @test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]
         test_cached_index(A)
 
+        @test_throws ErrorException("BSplineApprox needs at least d + 1, i.e. 3 control points.") BSplineApprox(
+            u, t, 2, 2, :Uniform, :Uniform)
+        @test_throws ErrorException("BSplineApprox needs at least d + 1, i.e. 4 control points.") BSplineApprox(
+            u, t, 3, 3, :ArcLen, :Average)
+        @test_nowarn BSplineApprox(u, t, 2, 3, :Uniform, :Uniform)
+
         # Test extrapolation
         A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform; extrapolate = true)
         @test A(-1.0) == u[1]

test/online_tests.jl

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