juliaintervals / taylormodels.jl Goto Github PK

Rigorous function approximation using Taylor models in Julia

License: Other

Julia 100.00%

taylormodels.jl's Introduction

TaylorModels

This package combines the IntervalArithmetic.jl and TaylorSeries.jl packages to provide Taylor models, i.e. Taylor polynomials with guaranteed error bounds to approximate functions.

An introduction is available in this video from JuliaCon 2018.

Authors

Luis Benet, Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM)
David P. Sanders, Departamento de Física, Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM)

Bibliography

Kyoko Makino, Rigorous Analysis of Nonlinear Motion in Particle Accelerators, PhD Thesis, Michigan State University (1991).
Mioara Joldes, Rigorous Polynomial Approximations and Applications, PhD Thesis, École Normale Supérieure de Lyon, ENS-Lyon (2011).

Acknowledgements

Financial support is acknowledged from DGAPA-UNAM PAPIIT grants IN-117117, IG-100616 and IG-100819. DPS acknowledges support through a Cátedra Marcos Moshinsky (2018).

taylormodels.jl's People

Contributors

Stargazers

Watchers

Forkers

lbenet mforets shreyanshtripathi barpit20 aa25desh fireofearth bibekkabi mewilhel ulinares lucaferranti standardgalactic mauricio-aguirre biomedai edtrelo

taylormodels.jl's Issues

Conflicting exports

julia> using TaylorModels

help?> derivative
WARNING: both IntervalRootFinding and TaylorIntegration export "derivative"; uses of it in module TaylorModels must be qualified
search: derivative

Couldn't find derivative
Perhaps you meant derivative or decoration
  No documentation found.

  Binding TaylorModels.derivative does not exist.

help?> ∇
WARNING: both IntervalRootFinding and TaylorIntegration export "∇"; uses of it in module TaylorModels must be qualified
"∇" can be typed by \nabla<tab>

search: ∇

Couldn't find ∇
Perhaps you meant ∇, ±, ×, ∅, ∞, ≺, ≼, ⊂, ⊃, ⊇, ⪽, !, %, &, *, +, -, /, :, <, >, \, ^, |, ~, ÷, π, ℯ, ∈ or ∉
  No documentation found.

  Binding TaylorModels.∇ does not exist.

Add a constructor for a TaylorModelN just changing remainder

TaylorModelN(u, Δ) = TaylorModelN(u.pol, Δ, u.x0, u.I)

Add linear dominated bounder method

See https://bt.pa.msu.edu/pub/papers/GOM05/GOM05.pdf.

Normalize to all variables in [-1, 1]

@mforets Do you have code for this?

Algorithm comparisons: Brusselator model

This might be an interesting problem to study the trade-off between algorithm choices. Consider the Brusselator model. Below I compared validated_integ and validated_integ2.

const A = 1.0
const B = 1.5
const B1 = B + 1

@taylorize function brusselator!(du, u, p, t)
    x, y = u
    x² = x * x
    aux = x² * y
    du[1] = A + aux - B1*x
    du[2] = B*x - aux
    return du
end

# set of initial states
U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)

# initial-value problem
bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim: 2)

prob = bruss(0.1)

# validated_integ2
@time sol_3 = solve(prob, T=25.0, alg=TMJets(orderT=6, orderQ=3, adaptive=false, absorb=false));
  2.283377 seconds (16.94 M allocations: 1.497 GiB, 44.55% gc time)

# validated_integ
@time sol_2 = solve(prob, T=30.0, alg=TMJets21a(orderT=6, orderQ=2));
  1.345497 seconds (4.99 M allocations: 418.477 MiB, 68.60% gc time)

However, the first algorithm fails if the time window is slightly larger:

@time sol_3 = solve(prob, T=28.0, alg=TMJets(orderT=6, orderQ=3, adaptive=false, absorb=false));

I'm confused as to why this happens, since you can check by plotting with plot(sol_3, vars=(1, 3)) that the flowpipe is effectively shrinking, so my gut feeling is that it would be easier to validate those final steps.

Add simplify function

Add a simplify function given a cutoff threshold such that the terms with small enough coefficients, given the threshold, are absorbed to its interval enclosure onto the remainder part.

Rename TMs domain field

Recall that:

            pol  :: Taylor1{T}    # polynomial approx (of order `ord`)
            rem  :: Interval{S}   # remainder
            x0   :: Interval{S}   # expansion point
            I    :: Interval{S}   # interval of interest

However, in several references it is the (interval) remainder that's denoted with I, and the domain is often written D. Also, the domain of interest is, in principle, not necessarily a box, while the remainder always is an interval (by definition of Taylor model).

Proposal: call the domain field dom to avoid potential sources of confusion. This is also more symmetric wrt to calling rem the remainder field.

integrate not defined for TaylorModelN

This used to be in TaylorNModel?

Better bound for Taylor Model product

Playing around with the implementation of the product of Taylor Models given in [1],
I found that his implementation gives a better bound for the resulting Taylor Model, eg, the TaylorModel1 product for f(x) = sin(x), g(x)=cos(x) + x^3, gives the following results

Green boxes correspond to the case where the diameter of the bound given by the current implementation is equal
to the diameter of the bound given by the implementation in [1].
Red boxes correspond to the case where the diameter of the bound given by [1] is less than the bound given by the
current implementation.

What do you think about this @lbenet?

[1] A Taylor model toolbox for solving ODEs implemented in MATLAB/INTLAB, Floarian Bünger

New methods for TaylorModel1 with TaylorN coefficients

There isn't an implemented version of TaylorModel1 with TaylorN (or TaylorModelN) coefficients for the following functions

Trigonometric functions (sin(x), cos(x), etc.)
Exponential and logarithmic functions.
Power functions (up to now only ^2 is supported)

These kinds of Taylor Models are required by the new ODE solver in #84 and will be needed by it in order to properly work.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Add benchmark suite

See eg. Tables 3-6 in AGK18.

Don't Understand How to Use TaylorModelN

I don't understand how to use TaylorModelN. The documentation is insufficient. Specifying the order, the approximation point, and the interval over which it should be valid makes sense, but what else should I need to specify? And why does it seem like I can't evaluate an ND Taylor model over more than 1D?

Add a contractor expanding around midpoint by default

TaylorModel1(n, X) = TaylorModel1(n, interval(mid(X)), X)

Tag v0.1.1

Revise random tests

See this paper and the discourse post. This is important since lots of tests we have use the naive implementation which is criticized.

Changes in TaylorSeries (linear_polynomial and nonlinear_polynomial)

There have been some changes in TS regarding linear_polynomial and the new functoin nonlinear_polynomial; the returned Taylor1 series have thee same order as the input one. We should update the functionality.

TaylorModels and complex evaluation

This issue is motivated by a discussion which was opened in TaylorSeries, see JuliaDiff/TaylorSeries.jl#288

Quoting @jonniedie:

I have a rational polynomial function G(s) = n(s) / d(s) where n(s) and d(s) are polynomial functions of s. The coefficients of these polynomials have uncertain parameters, which are represented by Intervals. For a given (real) value ω, I would like to evaluate the response G(ω*im) and get the argument and magnitude of the resulting complex number.

When I tried evaluating this directly using Intervals, the resulting bounds were too wide. @mforets recommended I give Taylor models a try, which is when I noticed Complex{<:TaylorModel}s wouldn't work because Complex requires <:Real parameters.

Update to TaylorSeries v0.12

You probably also want to enable CompatHelper again (here).
Sorry @lbenet for the déjà-vu 😅

Reactivate CompatHelper

Ambiguity in broadcasted(+, ::IntervalBox, ::IntervalBox)

(I was not sure whether to report in IntervalArithmetic.jl or here. Feel free to move around.)

Since a few days the following does not work anymore.

julia> using TaylorIntegration
julia> using TaylorModels: @taylorize
julia> using IntervalArithmetic: IntervalBox
julia> using TaylorModels: validated_integ

julia> @taylorize function vanderPol!(t, x, dx)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

julia> validated_integ(vanderPol!, [0., 1.], IntervalBox([0., 1.]), 0., 1., 2, 2, 1e-2)

ERROR: MethodError: Base.Broadcast.broadcasted(::typeof(+), ::IntervalBox{2,Float64}, ::IntervalBox{2,Float64}) is ambiguous. Candidates:
  broadcasted(f, X::IntervalBox, Y::IntervalBox) in IntervalArithmetic at .julia/packages/IntervalArithmetic/9q80F/src/multidim/arithmetic.jl:27
  broadcasted(f, X::IntervalBox, y) in IntervalArithmetic at .julia/packages/IntervalArithmetic/9q80F/src/multidim/arithmetic.jl:28
  broadcasted(f, x, Y::IntervalBox) in IntervalArithmetic at .julia/packages/IntervalArithmetic/9q80F/src/multidim/arithmetic.jl:29
  broadcasted(f::Union{typeof(*), typeof(+), typeof(muladd)}, arg1, arg2, args...) in DiffEqBase at .julia/packages/DiffEqBase/EFqMn/src/diffeqfastbc.jl:54
Possible fix, define
  broadcasted(::Union{typeof(*), typeof(+), typeof(muladd)}, ::IntervalBox, ::IntervalBox)
Stacktrace:
 [1] #validated_integ#31(::Int64, ::Bool, ::Function, ::Function, ::Function, ::Array{Float64,1}, ::IntervalBox{2,Float64}, ::Float64, ::Float64, ::Int64, ::Int64, ::Float64) at .julia/packages/TaylorModels/DXLdP/src/validatedODEs.jl:277
 [2] validated_integ(::Function, ::Array{Float64,1}, ::IntervalBox{2,Float64}, ::Float64, ::Float64, ::Int64, ::Int64, ::Float64) at .julia/packages/TaylorModels/DXLdP/src/validatedODEs.jl:265

It is quite mysterious because there were no changes in this package.

Add branch-and bound algorithm to bound polynomials

The idea is to have other methods to improve bounding polynomials, specially with many variables.

Add quadratic fast bounder

See https://bt.pa.msu.edu/pub/papers/GOM05/GOM05.pdf.

Update to TaylorIntegration v0.10

CompatHelper fell asleep again...

Bump TaylorModels

Commit f51a27d is not included yet in the latest tagged version but it is needed in Reachability.jl.

Range bounders for the upper bound

Currently, the linear_dominted_bounder and quadratic_fast_bounder only improves the lower bound, there isn't an implemented way to use these algorithms to improve the upper bound. I have two options in mind to solve this

Implement a modified version of the algorithms to address the upper bound case.
Make an auxiliary function whose purpose is to apply the algorithms to the negative of the Taylor Model which gives the desired results.

I think the last one is the best option since is easier to implement, but I don't know if this approach will give unexpected results for some cases.

Unexpected `evaluate` behavior with non-Float64 numbers

Is it expected that evaluate gives vastly different results when using non-Float64 numbers?

julia> using TaylorModels

julia> function eval(N)
    t = TaylorModels.Taylor1(3)
    q₁ = 1 + 2 * t + 2 * t^2
    I = interval(N(0), N(0))
    D = interval(N(-1), N(1))
    x0 = mid(D)
    TM = TaylorModels.TaylorModel1(q₁, I, x0, D)
    return evaluate(TM, domain(TM))
end;

julia> eval(Float64)
[-1, 5]

julia> eval(Float32)
[-3, 5]

julia> eval(Rational{Int})
[-3, 5]

Behaviour of the evaluate function

Dear all,
I have a question on the behavior of the evaluate function, or to be more precise, on the behavior of the iscontained
function.

At the moment it is defined as

iscontained(a, tm::TaylorModelN) = a ∈ domain(tm)-tm.x0

this leads to, in my opinion, unexpected behavior: as an example if tm is a Taylor model centered at 5, with domain [4,6],
then
tm(5)
fails the iscontained assertion, i.e., the evaluation fails on the expansion point.

Minimum failing code:

julia> tm = TaylorModel1(10, 5, Interval(4, 6))
 5 + 1 t + [0, 0]

julia> tm(5)
ERROR: AssertionError: iscontained(a, tm)

If this is a design choice of the package I think this should be documented.

Overload evaluate for TMSol

See this comment.

Getter function for the expansion point

I think there is no such function; it can be called expansion_point.

Cleanup auxiliary interval constructors

#103 introduces some auxiliary interval constructors. This is a reminder to use IntervalArithmetic.zero and IntervalArithmetic.symmetric_box once JuliaIntervals/IntervalArithmetic.jl#446 is ready.

Upgrade to new version of TaylorIntegration

... and tag a new version

Use more efficient multiplication of two Taylor models

That does not require duplicating the order of the Taylor series.

Automatic display of bechmark results

see #46 discussion.
With use of Weave.jl or Literate.jl

Generalize validated ODEs to receive an initial Taylor model

The method

TaylorModels.jl/src/validatedODEs.jl

Line 401 in 4dd5d17

function validated_integ(f!, qq0::AbstractArray{T,1}, δq0::IntervalBox{N,T},

currently receives (qq0, δq0) representing a box. It would be nice if it can handle an initial Taylor model reach-set (TaylorModel1 of TaylorN's); the conversion from a box to a TM can be done as a pre-processing step.

links to docs [stable] gives a 404

https://juliaintervals.github.io/TaylorModels.jl/stable --> 404

https://juliaintervals.github.io/TaylorModels.jl/latest/ --> works fine

Error when changing polynomial degree for `TaylorModel1{TaylorN{T} ,T}`

The following yields an error:

julia> using TaylorModels;

julia> xN, yN = set_variables("x y", order=4, numvars=2);

julia> xtm1_order5 = exp( TaylorModel1( Taylor1(xN, 5), 0..0, 0..0, -0.5 .. 0.5) )
  1.0 + 1.0 x + 0.5 x² + 0.16666666666666666 x³ + 0.041666666666666664 x⁴ + [-4.63936, 3.17106]

julia> xtm1_order4 = TaylorModel1( Taylor1(xN, 4), 0..0, 0..0, -0.5 .. 0.5)
  1.0 x + [0, 0]

julia> new_xtm1_order4, _ = TaylorSeries.fixorder(xtm1_order5, xtm1_order4)
ERROR: MethodError: no method matching TaylorModel1(::Taylor1{TaylorN{Float64}}, ::TaylorN{Interval{Float64}}, ::Interval{Float64}, ::Interval{Float64})
...

cc: @AnderGray

New struct to contain the solution of an ODE?

Documentation link broken in readme.md

Identity for floating poing rpa

Add the method fp_rpa(tm::TaylorModel1{T, T}) where {T} = tm (since the polynomial coeffs of the input are already floating points) and similarly for the other TM types.

Remark. This identity is needed if one tries to plot a TM1 with scalar coefficients:

m = 4
p = Taylor1([0.0, 1.0, 1.0], m)
δ = 0.5
I = Interval(-δ, δ)
x₀ = Interval(0.0)
D = Interval(-1.0, 1.0)
plot(TaylorModel1(p, I, x₀, D), label="")

Gives

MethodError: no method matching fp_rpa(::TaylorModel1{Float64,Float64})
...

and if you add TaylorModels.fp_rpa(tm::TaylorModel1{Float64,Float64}) = tm then:

Rewrite to allow the use of any polynomial package

I feel like it would be useful to rewrite the package to allow the use of any polynomial package, in particular for sparse representations.

Symbolic Taylor expansions for estimating round-off error

https://soarlab.org/papers/2018_toplas_sbbjrg.pdf

Should we use absorb_remainder for validated integration?

Shortcut constructor for independent variable

julia> using TaylorModels
[ Info: Recompiling stale cache file /Users/forets/.julia/compiled/v1.1/TaylorModels/yPZW1.ji for TaylorModels [a2574a24-d24b-5d8e-9d8b-7eb85f89fcf9]

julia> TaylorModel1(2, Interval(0.5, 0.5), Interval(0.0, 1.0))
 [0.5, 0.5] + [1, 1] t + [0, 0]

julia> TaylorModel1(2, 0.5, Interval(0.0, 1.0))
ERROR: MethodError: no method matching TaylorModel1(::Taylor1{Interval{Float64}}, ::Interval{Float64}, ::Float64, ::Interval{Float64})
Closest candidates are:
  TaylorModel1(::Taylor1{T}, ::Interval{S}, ::Interval{S}, ::Interval{S}) where {T, S} at /Users/forets/.julia/dev/TaylorModels/src/constructors.jl:33
  TaylorModel1(::Integer, ::Any, ::Interval{T}) where T at /Users/forets/.julia/dev/TaylorModels/src/constructors.jl:36
  TaylorModel1(::Interval{T}, ::Integer, ::Interval{T}, ::Interval{T}) where T at /Users/forets/.julia/dev/TaylorModels/src/constructors.jl:41
  ...
Stacktrace:
 [1] TaylorModel1(::Int64, ::Float64, ::Interval{Float64}) at /Users/forets/.julia/dev/TaylorModels/src/constructors.jl:36
 [2] top-level scope at none:0

Changing the penultimate argument in this line, x0, to interval(x0) one can use the short-cut constructor TaylorModel1(2, 0.5, Interval(0.0, 1.0)).

Can not use 'taylor_var' and 'make_Taylor_model'

julia> t = taylor_var(2, 0, -0.5..0.5)
ERROR: UndefVarError: taylor_var not defined
Stacktrace:
 [1] top-level scope at REPL[2]:1

julia> t1 = make_Taylor_model(exp, 3, 0, -0.5..0.5)
ERROR: UndefVarError: make_Taylor_model not defined
Stacktrace:
 [1] top-level scope at REPL[3]:1

I tried to run the examples but failed. These two functions 'taylor_var' and 'make_Taylor_model' seem not defined in this repo. I wonder whether there's any pkg I haven't include or import?

AssertionError in validated_integ if rhs has a sine term

using Test
using TaylorModels
using TaylorModels: validated_integ

function pendulum!(dx, x, p, t)
    aux = 2 * sin(x[1])
    dx[1] = x[2]
    dx[2] = aux + 8*x[3]
    dx[3] = zero(x[1])
    return nothing
end

# Initial conditions
tini, tend = 0.0, 0.03
q0 = [1.2, 0.2, -0.78]
δq0 = IntervalBox(0 .. 0, 3)
X0 = IntervalBox(q0 .+ δq0)

# Parameters
abstol = 1e-20
orderQ = 3
orderT = 10
ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))

normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol)

AssertionError: get_numvars() == length(vals)

Stacktrace:
  [1] _evaluate(a::TaylorN{Interval{Float64}}, vals::Tuple{Interval{Float64}})
    @ TaylorSeries ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:229
  [2] _evaluate
    @ ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:240 [inlined]
  [3] evaluate(a::TaylorN{Interval{Float64}}, vals::Interval{Float64}; sorting::Bool)
    @ TaylorSeries ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:223
  [4] evaluate
    @ ~/.julia/packages/TaylorSeries/tveWm/src/evaluate.jl:223 [inlined]
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [7] getindex
    @ ./broadcast.jl:575 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:983 [inlined]
 [11] copyto!
    @ ./broadcast.jl:936 [inlined]
 [12] copy
    @ ./broadcast.jl:908 [inlined]
 [13] materialize
    @ ./broadcast.jl:883 [inlined]
 [14] broadcasted
    @ ~/.julia/packages/IntervalArithmetic/dhnJB/src/multidim/arithmetic.jl:29 [inlined]
 [15] picard_remainder!(f!::typeof(pendulum!), t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xxI::Vector{Taylor1{TaylorN{Interval{Float64}}}}, dxxI::Vector{Taylor1{TaylorN{Interval{Float64}}}}, δI::IntervalBox{3, Float64}, δt::Interval{Float64}, Δx::IntervalBox{3, Float64}, Δ0::IntervalBox{3, Float64}, params::Nothing)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:108
 [16] remainder_taylorstep!(f!::Function, t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xI::Vector{Taylor1{Interval{Float64}}}, dxI::Vector{Taylor1{Interval{Float64}}}, δI::IntervalBox{3, Float64}, δt::Interval{Float64}, params::Nothing)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:32
 [17] validated_step!(f!::typeof(pendulum!), t::Taylor1{Float64}, x::Vector{Taylor1{TaylorN{Float64}}}, dx::Vector{Taylor1{TaylorN{Float64}}}, xaux::Vector{Taylor1{TaylorN{Float64}}}, tI::Taylor1{Float64}, xI::Vector{Taylor1{Interval{Float64}}}, dxI::Vector{Taylor1{Interval{Float64}}}, xauxI::Vector{Taylor1{Interval{Float64}}}, t0::Float64, tmax::Float64, sign_tstep::Int64, xTMN::Vector{TaylorModelN{3, Float64, Float64}}, xv::Vector{IntervalBox{3, Float64}}, rem::Vector{Interval{Float64}}, zbox::IntervalBox{3, Float64}, symIbox::IntervalBox{3, Float64}, nsteps::Int64, orderT::Int64, abstol::Float64, params::Nothing, parse_eqs::Bool, check_property::TaylorModels.var"#36#38")
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:361
 [18] validated_integ(f!::Function, X0::IntervalBox{3, Float64}, t0::Float64, tmax::Float64, orderQ::Int64, orderT::Int64, abstol::Float64, params::Nothing; maxsteps::Int64, parse_eqs::Bool, check_property::Function)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:533
 [19] validated_integ (repeats 2 times)
    @ ~/.julia/dev/TaylorModels/src/validatedODEs.jl:481 [inlined]
 [20] top-level scope
    @ In[15]:27
 [21] eval
    @ ./boot.jl:360 [inlined]
 [22] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

Let me remark that for aux = 2 there is no error.

MethodError in time-dependent ODE (Duffing oscillator)

The problem seems to be due to the presence of the time variable in the right-hand side of the ODEs using validated_integ2 (the example works with validated_integ).

using ReachabilityAnalysis, Plots

@taylorize function duffing!(du, u, p, t)
    local α = -1.0
    local β = 1.0
    local δ = 0.3
    local γ = 0.37

    x, v = u
    f = γ * cos(ω * t)

    du[1] = u[2]
    du[2] = - α*x - δ*v - β*x^3 + f
end

ω = 1.2
T = 2*pi / ω
X0 = Singleton([1.0, 0.0]) ⊕ BallInf(zeros(2), 0.1)
prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim=2);

# uses validated_integ2
sol = solve(prob, tspan=(0.0, 20*T), alg=TMJets());

MethodError: no method matching Float64(::Taylor1{Float64})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  (::Type{T})(::Base.TwicePrecision) where T<:Number at twiceprecision.jl:243
  ...

Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::Taylor1{Float64})
    @ Base ./number.jl:7
  [2] promote(a::TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}, b::Taylor1{Float64})
    @ TaylorModels ~/.julia/dev/TaylorModels/src/promotion.jl:85
  [3] +(b::TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}, a::Taylor1{Float64})
    @ TaylorSeries ~/.julia/packages/TaylorSeries/tveWm/src/arithmetic.jl:114
  [4] duffing!
    @ ./In[71]:13 [inlined]
  [5] duffing!(du::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, u::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, p::Nothing, t::Taylor1{Float64})
    @ Main ./none:0
  [6] picard_iteration(f!::typeof(duffing!), dx::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, xTM1K::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, params::Nothing, t::Taylor1{Float64}, x0::Vector{TaylorModels.TaylorModelN{2, Float64, Float64}}, box::IntervalBox{2, Float64}, #unused#::Val{true})
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:642
  [7] _validate_step!(xTM1K::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, f!::Function, dx::Vector{TaylorModels.TaylorModel1{TaylorN{Float64}, Float64}}, x0::Vector{TaylorModels.TaylorModelN{2, Float64, Float64}}, params::Nothing, x::Vector{Taylor1{TaylorN{Float64}}}, t::Taylor1{Float64}, box::IntervalBox{2, Float64}, dof::Int64, rem::Vector{IntervalArithmetic.Interval{Float64}}, abstol::Float64, δt::Float64, sign_tstep::Int64, E::Vector{IntervalArithmetic.Interval{Float64}}, E′::Vector{IntervalArithmetic.Interval{Float64}}, polv::Vector{Taylor1{TaylorN{Float64}}}, low_ratiov::Vector{Float64}, hi_ratiov::Vector{Float64}, adaptive::Bool, minabstol::Float64; ε::Float64, δ::Float64, validatesteps::Int64, extrasteps::Int64)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:696
  [8] validated_integ2(f!::typeof(duffing!), X0::IntervalBox{2, Float64}, t0::Float64, tf::Float64, orderQ::Int64, orderT::Int64, abstol::Float64, params::Nothing; parse_eqs::Bool, maxsteps::Int64, absorb::Bool, adaptive::Bool, minabstol::Float64, validatesteps::Int64, ε::Float64, δ::Float64, absorb_steps::Int64)
    @ TaylorModels ~/.julia/dev/TaylorModels/src/validatedODEs.jl:826
  [9] post(alg::TMJets21b{Float64, ZonotopeEnclosure}, ivp::InitialValueProblem{BlackBoxContinuousSystem{typeof(duffing!)}, MinkowskiSum{Float64, Singleton{Float64, Vector{Float64}}, BallInf{Float64, Vector{Float64}}}}, timespan::IntervalArithmetic.Interval{Float64}; Δt0::IntervalArithmetic.Interval{Float64}, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21b{Float64, ZonotopeEnclosure}}}})
    @ ReachabilityAnalysis ~/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/TMJets21b/post.jl:36
 [10] solve(::InitialValueProblem{BlackBoxContinuousSystem{typeof(duffing!)}, MinkowskiSum{Float64, Singleton{Float64, Vector{Float64}}, BallInf{Float64, Vector{Float64}}}}; kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21b{Float64, ZonotopeEnclosure}}}})
    @ ReachabilityAnalysis ~/.julia/dev/ReachabilityAnalysis/src/Continuous/solve.jl:61
 [11] top-level scope
    @ In[71]:23
 [12] eval
    @ ./boot.jl:360 [inlined]
 [13] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

Assert (or warn) if evaluation is computed in the proper domain

Currently, evaluate does not check if the evaluation point (or interval) is contained within the domain defining the Taylor model.

julia> using TaylorModels

julia> tm = TaylorModel1(5, 0.0, -0.5 .. 0.5)
 1.0 t + [0, 0]

julia> em = exp(tm)
 1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + [-0, 2.33541e-05]

julia> domain(em)
[-0.5, 0.5]

julia> evaluate(em, 1.5)
[4.46171, 4.46175]

There is no assertion nor warning about evaluating outside the domain; the result in this case is actually wrong. The problem is that this happens silently.

We should use @assert, or at least @warn when the point (or interval) used in the evaluation is not completely contained within the domain of the Taylor model.

cc: @dpsanders, @UzielLinares

TaylorSeries `numtype` and failing tests

In JuliaDiff/TaylorSeries.jl#289, the behavior of eltype is changed for Taylor1 and TaylorN; see JuliaDiff/TaylorSeries.jl#224. The old behavior is obtained by using numtype (which in TaylorSeries.jl is not exported); something similar was done in JuliaIntervals/IntervalArithmetic.jl#496.

That change in TaylorSeries makes some tests to fail in TayloModels.jl.

Return TaylorModel1 from validated_integ

As we discussed elsewhere, since we have that info for each time step covered, it is a good idea to return it.