Increase functionality of RTaylorModel1 (#131)

* Check monotonicity bounds and exploit dispatch

* Functionality for RTaylorModel1{TaylorN} and tests

* Fixes and add test coverage

* Avoid duplicating code for integrate, fixes, docstrings and tests

* Add centered_dom

* Bump patch version

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorModels"
 uuid = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 repo = "https://github.com/JuliaIntervals/TaylorModels.jl.git"
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"

src/TaylorModels.jl

@@ -31,7 +31,7 @@ import LinearAlgebra: norm
 export TaylorModel1, RTaylorModel1, TaylorModelN, TMSol
 
 export remainder, polynomial, domain, expansion_point, flowpipe, get_xTM,
-    rpa, fp_rpa, bound_remainder,
+    rpa, fp_rpa, bound_remainder, centered_dom,
     validated_integ, validated_integ2
 
 export linear_dominated_bounder, quadratic_fast_bounder

src/arithmetic.jl

@@ -61,7 +61,7 @@ for TM in tupleTMs
             a_order = get_order(a)
             b_order = get_order(b)
             rnegl_order = a_order + b_order
-            aux = a.dom - a.x0
+            aux = centered_dom(a)
 
             # Returned polynomial
             res = a.pol * b.pol
@@ -70,7 +70,7 @@ for TM in tupleTMs
             # Remainder of the product
             if $TM == TaylorModel1
 
-                # Remaing terms of the product as reduced Taylor1 (factored polynomial)
+                # Remaining terms of the product as reduced Taylor1 (factored polynomial)
                 rnegl = Taylor1(zero(res[0]), rnegl_order)
                 for k in order+1:rnegl_order
                     @inbounds for i = 0:k
@@ -84,7 +84,7 @@ for TM in tupleTMs
                 Δ = remainder_product(a, b, aux, Δnegl)
             else
 
-                # Remaing terms of the product as reduced Taylor1 (factored polynomial)
+                # Remaining terms of the product as reduced Taylor1 (factored polynomial)
                 rnegl = Taylor1(zero(res[0]), rnegl_order-order)
                 for k in order+1:rnegl_order
                     @inbounds for i = 0:k
@@ -122,6 +122,7 @@ for TM in tupleTMs
 end
 
 # Remainder of the product
+# TaylorModel1
 function remainder_product(a, b, aux, Δnegl)
     Δa = a.pol(aux)
     Δb = b.pol(aux)
@@ -135,12 +136,6 @@ function remainder_product(a::TaylorModel1{TaylorN{T}, S},
     # An N-dimensional symmetrical IntervalBox is assumed
     # to bound the TaylorN part
     auxQ = IntervalBox(-1 .. 1, Val(N))
-    Δ = remainder_product(a, b, auxT, auxQ, Δnegl)
-    return Δ
-end
-function remainder_product(a::TaylorModel1{TaylorN{T}, S},
-                           b::TaylorModel1{TaylorN{T}, S},
-                           auxT, auxQ, Δnegl) where {T, S}
     Δa = a.pol(auxT)(auxQ)
     Δb = b.pol(auxT)(auxQ)
     Δ = Δnegl(auxQ) + Δb * a.rem + Δa * b.rem + a.rem * b.rem
@@ -153,17 +148,32 @@ function remainder_product(a::TaylorModel1{TaylorModelN{N,T,S},S},
     Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem
 
     # Evaluate at the TMN centered domain
-    auxN = a[0].dom - a[0].x0
+    auxN = centered_dom(a[0])
     ΔN = Δ(auxN)
     return ΔN
 end
-function remainder_product(a::RTaylorModel1, b::RTaylorModel1, aux, Δnegl, order)
+# RTaylorModel1
+function remainder_product(a, b, aux, Δnegl, order)
     Δa = a.pol(aux)
     Δb = b.pol(aux)
     V = aux^(order+1)
     Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem * V
     return Δ
 end
+function remainder_product(a::RTaylorModel1{TaylorN{T},S},
+                           b::RTaylorModel1{TaylorN{T},S},
+                           aux, Δnegl, order) where {T, S}
+    N = get_numvars()
+    # An N-dimensional symmetrical IntervalBox is assumed
+    # to bound the TaylorN part
+    auxQ = symmetric_box(N, T)
+    Δa = a.pol(aux)(auxQ)
+    Δb = b.pol(aux)(auxQ)
+    V = aux^(order+1)
+    Δ = Δnegl(auxQ) + Δb * a.rem + Δa * b.rem + a.rem * b.rem * V
+    return Δ
+end
+
 
 # Division
 function /(a::TaylorModel1, b::TaylorModel1)
@@ -207,7 +217,7 @@ function truncate_taylormodel(a::RTaylorModel1, m::Integer)
 
     apol = Taylor1(copy(a.pol.coeffs[1:m+1]))
     bpol = Taylor1(copy(a.pol.coeffs))
-    aux = a.dom - a.x0
+    aux = centered_dom(a)
     Δnegl = bound_truncation(RTaylorModel1, bpol, aux, m)
     Δ = Δnegl + a.rem * (aux)^(order-m)
     return RTaylorModel1( apol, Δ, a.x0, a.dom )
@@ -253,7 +263,7 @@ function *(a::TaylorModelN, b::TaylorModelN)
     b_order = get_order(b)
     rnegl_order = a_order + b_order
     @assert rnegl_order ≤ get_order()
-    aux = a.dom - a.x0
+    aux = centered_dom(a)
 
     # Returned polynomial
     res = a.pol * b.pol

src/auxiliary.jl

@@ -55,7 +55,7 @@ for TM in tupleTMs
             apol, bpol = TaylorSeries.fixorder(apol0, bpol0)
 
             # Bound for the neglected part of the polynomial
-            dom = a.dom - a.x0
+            dom = centered_dom(a)
             Δa = bound_truncation($TM, apol0, dom, order) + remainder(a)
             Δb = bound_truncation($TM, bpol0, dom, order) + remainder(b)
 
@@ -86,7 +86,7 @@ function fixorder(a::TaylorModelN, b::TaylorModelN)
     apol, bpol = TaylorSeries.fixorder(apol0, bpol0)
 
     # Bound for the neglected part of the polynomial
-    dom = a.dom - a.x0
+    dom = centered_dom(a)
     Δa = bound_truncation(TaylorModelN, apol0, dom, order) + remainder(a)
     Δb = bound_truncation(TaylorModelN, bpol0, dom, order) + remainder(b)
 

src/bounds.jl

@@ -18,43 +18,16 @@ function bound_remainder(::Type{TaylorModel1}, f::Function, polf::Taylor1, polfI
 
     _order = get_order(polf) + 1
     fTIend = polfI[_order]
-    if (sup(fTIend) < 0 || inf(fTIend) > 0)
-        # Absolute remainder is monotonic
-        a = Interval(inf(I))
-        b = Interval(sup(I))
-        Δlo = f(a) - polf(a-x0)
-        # Δlo = f(a) - bound_taylor1(polf, a-x0)
-        Δhi = f(b) - polf(b-x0)
-        # Δhi = f(b) - bound_taylor1(polf, b-x0)
-        Δx0 = f(x0) - polf[0]
-        Δ = hull(Δlo, Δx0, Δhi)
-    else
-        # Lagrange bound
-        Δ = fTIend * (I-x0)^_order
-    end
-    return Δ
+    bb = sup(fTIend) < 0 || inf(fTIend) > 0
+    return _monot_bound_remainder(TaylorModel1, Val(bb), f, polf, polfI, x0, I)
 end
 function bound_remainder(::Type{TaylorModel1}, f::Function, polf::Taylor1{TaylorN{T}}, polfI::Taylor1, x0, I::Interval) where {T}
     # The domain of the TaylorN part is assumed to be
     # the normalized symmetric box
     _order = get_order(polf) + 1
     fTIend = polfI[_order]
-    if (sup(fTIend) < 0 || inf(fTIend) > 0)
-        # Absolute remainder is monotonic
-        a = Interval(inf(I))
-        b = Interval(sup(I))
-        symIbox = IntervalBox(-1 .. 1, get_numvars())
-        Δlo = (f(a) - polf(a-x0))(symIbox)
-        # Δlo = f(a) - bound_taylor1(polf, a-x0)
-        Δhi = (f(b) - polf(b-x0))(symIbox)
-        # Δhi = f(b) - bound_taylor1(polf, b-x0)
-        Δx0 = (f(x0) - polf[0])(symIbox)
-        Δ = hull(Δlo, Δx0, Δhi)
-    else
-        # Lagrange bound
-        Δ = fTIend * (I-x0)^_order
-    end
-    return Δ
+    bb = sup(fTIend) < 0 || inf(fTIend) > 0
+    return _monot_bound_remainder(TaylorModel1, Val(bb), f, polf, polfI, x0, I)
 end
 
 
@@ -77,22 +50,98 @@ function bound_remainder(::Type{RTaylorModel1}, f::Function, polf::Taylor1, polf
     fTIend = polfI[_order+1]
     a = Interval(inf(I))
     b = Interval(sup(I))
-    if (sup(fTIend) < 0 || inf(fTIend) > 0) && isempty(a ∩ x0) && isempty(b ∩ x0)
-        # Error is monotonic
-        denom_lo = (a-x0)^_order
-        Δlo = f(a) - polf(a-x0)
-        # Δlo = f(a) - bound_taylor1(polf, a-x0)
-        Δlo = Δlo / denom_lo
-        denom_hi = (b-x0)^_order
-        Δhi = f(b) - polf(b-x0)
-        # Δhi = f(b) - bound_taylor1(polf, b-x0)
-        Δhi = Δhi / denom_hi
-        Δ = hull(Δlo, Δhi)
-    else
-        # Lagrange coefficient
-        Δ = polfI[_order]
-    end
-    return Δ
+    bb = (sup(fTIend) < 0 || inf(fTIend) > 0) && isempty(a ∩ x0) && isempty(b ∩ x0)
+    return _monot_bound_remainder(RTaylorModel1, Val(bb), f, polf, polfI, x0, I)
+end
+
+
+"""
+    _monot_bound_remainder(::Type{TaylorModel1}, ::Val{true}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+
+Computes the remainder exploiting monotonicity; see Prop 2.2.1 in Mioara Joldes' PhD thesis (pp 52).
+"""
+@inline function _monot_bound_remainder(::Type{TaylorModel1}, ::Val{true}, f::Function,
+        polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+    # Absolute remainder is monotonic
+    a = Interval(inf(I))
+    b = Interval(sup(I))
+    Δlo = f(a) - polf(a-x0)
+    # Δlo = f(a) - bound_taylor1(polf, a-x0)
+    Δhi = f(b) - polf(b-x0)
+    # Δhi = f(b) - bound_taylor1(polf, b-x0)
+    Δx0 = f(x0) - polf[0]
+    return hull(Δlo, Δx0, Δhi)
+end
+@inline function _monot_bound_remainder(::Type{TaylorModel1}, ::Val{true}, f::Function,
+        polf::Taylor1{TaylorN{T}}, polfI::Taylor1, x0, I::Interval) where {T}
+    # Absolute remainder is monotonic
+    a = Interval(inf(I))
+    b = Interval(sup(I))
+    symIbox = IntervalBox(-1 .. 1, get_numvars())
+    Δlo = (f(a) - polf(a-x0))(symIbox)
+    # Δlo = f(a) - bound_taylor1(polf, a-x0)
+    Δhi = (f(b) - polf(b-x0))(symIbox)
+    # Δhi = f(b) - bound_taylor1(polf, b-x0)
+    Δx0 = (f(x0) - polf[0])(symIbox)
+    return hull(Δlo, Δx0, Δhi)
+end
+"""
+    _monot_bound_remainder(::Type{TaylorModel1}, ::Val{false}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+
+Computes the remainder using Lagrange bound
+"""
+@inline function _monot_bound_remainder(::Type{TaylorModel1}, ::Val{false}, f::Function,
+        polf, polfI, x0, I)
+    _order = get_order(polf) + 1
+    fTIend = polfI[_order]
+    # Lagrange bound
+    return fTIend * (I-x0)^_order
+end
+
+"""
+    _monot_bound_remainder(::Type{RTaylorModel1}, ::Val{true}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+
+Computes the remainder exploiting monotonicity; see Prop 2.3.7 in Mioara Joldes' PhD thesis (pp 67).
+"""
+@inline function _monot_bound_remainder(::Type{RTaylorModel1}, ::Val{true}, f::Function,
+        polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+
+    _order = get_order(polf) + 1
+    a = Interval(inf(I))
+    b = Interval(sup(I))
+    # Error is monotonic
+    denom_lo = (a-x0)^_order
+    Δlo = f(a) - polf(a-x0)
+    # Δlo = f(a) - bound_taylor1(polf, a-x0)
+    Δlo = Δlo / denom_lo
+    denom_hi = (b-x0)^_order
+    Δhi = f(b) - polf(b-x0)
+    # Δhi = f(b) - bound_taylor1(polf, b-x0)
+    Δhi = Δhi / denom_hi
+    return hull(Δlo, Δhi)
+end
+@inline function _monot_bound_remainder(::Type{RTaylorModel1}, ::Val{true}, f::Function,
+        polf::Taylor1{TaylorN{T}}, polfI::Taylor1, x0, I::Interval) where {T}
+
+    _order = get_order(polf) + 1
+    a = Interval(inf(I))
+    b = Interval(sup(I))
+    symIbox = IntervalBox(-1 .. 1, get_numvars())
+    # Error is monotonic
+    denom_lo = (a-x0)^_order
+    Δlo = (f(a) - polf(a-x0))(symIbox)
+    # Δlo = f(a) - bound_taylor1(polf, a-x0)
+    Δlo = Δlo / denom_lo
+    denom_hi = (b-x0)^_order
+    Δhi = (f(b) - polf(b-x0))(symIbox)
+    # Δhi = f(b) - bound_taylor1(polf, b-x0)
+    Δhi = Δhi / denom_hi
+    return hull(Δlo, Δhi)
+end
+@inline function _monot_bound_remainder(::Type{RTaylorModel1}, ::Val{false}, f::Function,
+        polf::Taylor1, polfI::Taylor1, x0, I::Interval)
+    _order = get_order(polf) + 1
+    return polfI[_order]
 end
 
 

src/constructors.jl

@@ -53,6 +53,8 @@ for TM in tupleTMs
         @inline polynomial(tm::$TM) = tm.pol
         @inline domain(tm::$TM) = tm.dom
         @inline expansion_point(tm::$TM) = tm.x0
+        # Centered domain
+        @inline centered_dom(tm::$TM) = domain(tm) - expansion_point(tm)
     end
 end
 
@@ -147,6 +149,8 @@ TaylorModelN(a::T, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) wh
 @inline domain(tm::TaylorModelN) = tm.dom
 @inline expansion_point(tm::TaylorModelN) = tm.x0
 @inline get_numvars(::TaylorModelN{N,T,S}) where {N,T,S} = N
+# Centered domain
+@inline centered_dom(tm::TaylorModelN) = domain(tm) .- expansion_point(tm)
 
 """
     TMSol{N,T,V1,V2,M}
@@ -158,7 +162,7 @@ Structure containing the solution of a validated integration.
 
     `fp   :: AbstractVector{IntervalBox{N,T}}`  IntervalBox vector representing the flowpipe
 
-    `xTMv :: AbstractMatrix{TaylorModel1{TaylorN{T},T}}`  Matrix whose entry `xTMv[i,t]` represents 
+    `xTMv :: AbstractMatrix{TaylorModel1{TaylorN{T},T}}`  Matrix whose entry `xTMv[i,t]` represents
     the `TaylorModel1` of the i-th dependent variable, obtained at time time[t].
 """
 struct TMSol{N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
@@ -167,7 +171,7 @@ struct TMSol{N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}
     fp   :: V2
     xTM  :: M
 
-    function TMSol(time::V1, fp::V2, xTM::M) where 
+    function TMSol(time::V1, fp::V2, xTM::M) where
             {N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
             M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
         @assert length(time) == length(fp) == size(xTM,2) && N == size(xTM,1)