Adapt validated integration methods to new taylorize macro (#133)

* Adapt validated integration methods to new taylorize macro

* Fixes of broken tests for Julia v1.8

* Bump minor version and version of TaylorIntegration

* More generic evaluate methods for TMNs

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.4"
+version = "0.6.0"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -19,9 +19,9 @@ IntervalArithmetic = "^0.20"
 IntervalRootFinding = "0.5"
 RecipesBase = "1"
 Reexport = "1"
-TaylorIntegration = "0.8"
+TaylorIntegration = "0.9"
 TaylorSeries = "0.12"
-julia = "1.5, 1.6, 1.7"
+julia = "1.6, 1.7"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/evaluate.jl

@@ -78,7 +78,7 @@ end
 
 (tm::TaylorModelN{N,T,S})(a::IntervalBox{N,S}) where {N,T,S} = evaluate(tm, a)
 
-function evaluate(tm::TaylorModelN{N,T,S}, a::Array{R,1}) where {N,T,S,R}
+function evaluate(tm::TaylorModelN{N,T,S}, a::AbstractVector{R}) where {N,T,S,R}
     @assert iscontained(a, tm)
     _order = get_order(tm)
 
@@ -87,7 +87,7 @@ function evaluate(tm::TaylorModelN{N,T,S}, a::Array{R,1}) where {N,T,S,R}
     return tm.pol(a) + Δ
 end
 
-(tm::TaylorModelN{N,T,S})(a::Array{R,1}) where {N,T,S,R} = evaluate(tm, a)
+(tm::TaylorModelN{N,T,S})(a::AbstractVector{R}) where {N,T,S,R} = evaluate(tm, a)
 
 evaluate(tm::Vector{TaylorModelN{N,T,S}}, a::IntervalBox{N,S}) where {N,T,S} =
     IntervalBox( [ tm[i](a) for i in eachindex(tm) ] )

src/validatedODEs.jl

@@ -324,20 +324,19 @@ function validated_step!(f!, t::Taylor1{T}, x::Vector{Taylor1{TaylorN{T}}},
         t0::T, tmax::T, sign_tstep::Int,
         xTMN::Vector{TaylorModelN{N,T,T}}, xv::Vector{IntervalBox{N,T}},
         rem::Vector{Interval{T}}, zbox::IntervalBox{N,T}, symIbox::IntervalBox{N,T},
-        nsteps::Int, orderT::Int, abstol::T, params, parse_eqs::Bool,
+        nsteps::Int, orderT::Int, abstol::T, params,
+        tmpTaylor, arrTaylor, tmpTaylorI, arrTaylorI, parse_eqs::Bool,
         adaptive::Bool, minabstol::T, absorb::Bool,
         check_property::Function=(t, x)->true) where {N,T}
 
     # One step integration (non-validated)
-    # TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params)
-    # δt = TaylorIntegration.stepsize(x, abstol)
-    δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs)
+    δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params,
+                                        tmpTaylor, arrTaylor, parse_eqs)
     f!(dx, x, params, t)  # Update `dx[:][orderT]`
 
-    # One step integration for the initial box
-    # TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, tI, xI, dxI, xauxI, params)
-    # δtI = TaylorIntegration.stepsize(xI, abstol)
-    δtI = TaylorIntegration.taylorstep!(f!, tI, xI, dxI, xauxI, abstol, params, parse_eqs)
+    # One step integration for the *initial box*
+    δtI = TaylorIntegration.taylorstep!(f!, tI, xI, dxI, xauxI, abstol, params,
+                                        tmpTaylorI, arrTaylorI, parse_eqs)
     f!(dxI, xI, params, tI)  # Update `dxI[:][orderT+1]`
 
     # Step size
@@ -574,7 +573,11 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
     @inbounds xv[1] = evaluate(xTMN, S)
 
     # Determine if specialized jetcoeffs! method exists (built by @taylorize)
-    parse_eqs = TaylorIntegration._determine_parsing!(parse_eqs, f!, t, x, dx, params)
+    parse_eqsI, tmpTaylorI, arrTaylorI = TaylorIntegration._determine_parsing!(
+        parse_eqs, f!, tI, xI, dxI, params)
+    parse_eqs, tmpTaylor, arrTaylor = TaylorIntegration._determine_parsing!(
+        parse_eqs, f!, t, x, dx, params)
+        @assert parse_eqs == parse_eqsI
 
     local _success # if true, the validation step succeeded
     red_abstol = abstol
@@ -587,6 +590,7 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
         (_success, δt, red_abstol) = validated_step!(f!, t, x, dx, xaux, tI, xI, dxI, xauxI,
                                 t0, tmax, sign_tstep, xTMN, xv, rem, zB, S,
                                 nsteps, orderT, red_abstol, params,
+                                tmpTaylor, arrTaylor, tmpTaylorI, arrTaylorI,
                                 parse_eqs, adaptive, minabstol,
                                 absorb, check_property)
         δtI = sign_tstep * Interval(zt, sign_tstep*δt)
@@ -598,7 +602,7 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
         @inbounds t[0] = t0
         @inbounds tI[0] = t0
         @. begin
-            xTM1v[:, nsteps] = TaylorModel1(x, rem, zI, δtI)  # deepcopy?
+            xTM1v[:, nsteps] = TaylorModel1(deepcopy(x), rem, zI, δtI) # deepcopy is needed!
             x = Taylor1(evaluate(x, δt), orderT)
             # dx = Taylor1(zero(constant_term(x)), orderT)
             xI = Taylor1(evaluate(xTMN, (S,)), orderT+1)
@@ -829,11 +833,13 @@ function validated_integ2(f!, X0, t0::T, tf::T, orderQ::Int, orderT::Int,
     @inbounds tv[1] = t0
     @inbounds xv[1] = evaluate(xTMN, S)
 
-    parse_eqs = TaylorIntegration._determine_parsing!(parse_eqs, f!, t, x, dx, params)
+    parse_eqs, tmpTaylor, arrTaylor = TaylorIntegration._determine_parsing!(parse_eqs,
+        f!, t, x, dx, params)
     red_abstol = abstol
 
     while t0 * sign_tstep < tf * sign_tstep
-        δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params, parse_eqs)
+        δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux, abstol, params,
+            tmpTaylor, arrTaylor, parse_eqs)
         f!(dx, x, params, t)
 
         δt = min(δt, sign_tstep*(tf-t0))

test/TM1.jl

@@ -418,7 +418,7 @@ end
         tma = acos(tv)
         tmb = cos(tma)
         @test tmb == cos(acos(tv))
-        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-16
+        @test sup(norm(tmb.pol - tv.pol, Inf)) < 5.0e-16
 
         tma = tan(tv)
         tmb = atan(tma)
@@ -959,7 +959,7 @@ end
         tma = acos(tv)
         tmb = cos(tma)
         @test tmb == cos(acos(tv))
-        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-16
+        @test sup(norm(tmb.pol - tv.pol, Inf)) < 5.0e-16
 
         tma = tan(tv)
         tmb = atan(tma)

test/TMN.jl

@@ -230,7 +230,7 @@ set_variables(Interval{Float64}, [:x, :y], order=_order_max)
         tma = acos(ym)
         tmb = cos(tma)
         @test tmb == cos(acos(ym))
-        @test sup(norm(tmb.pol - ym.pol, Inf)) < 1.0e-16
+        @test sup(norm(tmb.pol - ym.pol, Inf)) < 5.0e-16
 
         tma = tan(xm)
         tmb = atan(tma)

test/validated_integ.jl

@@ -30,7 +30,7 @@ function test_integ(fexact, t0, qTM, q0, δq0)
     # Box computed from the exact solution must be within q
     bb = all(fexact(t0+δtI, q0 .+ q0ξB) .⊆ q)
     # Display details if bb is false
-    bb || @show(t0, domt, remainder.(qTM), 
+    bb || @show(t0, domt, remainder.(qTM),
             δt, δtI, q0ξ, q0ξB, q,
             fexact(t0+δtI, q0 .+ q0ξB))
     return bb
@@ -76,6 +76,7 @@ end
                 @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
             end
 
+            # Check equality of solutions using `parse_eqs=false` or `parse_eqs=true`
             solf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                 adaptive=false)
             qvf, qTMf = getfield.((solf,), 2:3)
@@ -118,6 +119,13 @@ end
                 @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
             end
 
+            # Check equality of solutions using `parse_eqs=false` or `parse_eqs=true`
+            solf = validated_integ2(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
+                adaptive=false)
+            qvf, qTMf = getfield.((solf,), 2:3)
+            @test length(qvf) == length(qv)
+            @test qTM == qTMf
+
             # initializaton with a Taylor model
             X0tm = get_xTM(sol,1)
             sol2 = validated_integ2(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)