Merge 11291df into 7010459

Project.toml

@@ -8,6 +8,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
@@ -24,7 +25,8 @@ julia = "1.3, 1.4, 1.5, 1.6"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "Test"]
+test = ["LinearAlgebra", "Random", "Test"]

src/TaylorModels.jl

@@ -21,7 +21,7 @@ import Base: setindex!, getindex, copy, firstindex, lastindex,
 using TaylorSeries: derivative, ∇
 
 import TaylorSeries: integrate, get_order, evaluate, pretty_print,
-    constant_term, linear_polynomial, fixorder
+    constant_term, linear_polynomial, fixorder, get_numvars
 
 import IntervalArithmetic: showfull
 
@@ -31,9 +31,9 @@ import LinearAlgebra: norm
 #         taylor1_var, integrate, degree,
 #         calculate_set, Taylor_step
 
-export TaylorModel1, RTaylorModel1, TaylorModelN
+export TaylorModel1, RTaylorModel1, TaylorModelN, TMSol
 
-export remainder, polynomial, domain, expansion_point,
+export remainder, polynomial, domain, expansion_point, flowpipe, get_xTM,
     rpa, fp_rpa, bound_remainder,
     validated_integ, validated_integ2
 
@@ -47,9 +47,9 @@ include("evaluate.jl")
 include("rpa_functions.jl")
 include("arithmetic.jl")
 include("integration.jl")
-include("recipe.jl")
 include("show.jl")
 include("validatedODEs.jl")
+include("recipe.jl")
 
 # include("Taylor1/Taylor1.jl")
 # include("TaylorN/TaylorN.jl")

src/auxiliary.jl

@@ -104,3 +104,18 @@ end
 @inline zero_interval(T) = zero(Interval{T})
 @inline zero_box(N, T) = IntervalBox(zero_interval(T), Val(N))
 @inline symmetric_box(N, T) = IntervalBox(Interval{T}(-1, 1), Val(N))
+
+
+# TMSol utilities
+@inline firstindex(a::TMSol) = firstindex(a.time)
+@inline lastindex(a::TMSol)  = lastindex(a.time)
+@inline Base.length(a::TMSol) = length(a.time)
+@inline Base.iterate(a::TMSol, state=0) = state ≥ lastindex(a) ? nothing : (a[state+1], state+1)
+@inline Base.eachindex(a::TMSol) = firstindex(a):lastindex(a)
+
+getindex(a::TMSol, n::Integer) = a.xTM[:,n]
+getindex(a::TMSol, u::UnitRange) = a.xTM[:,u]
+getindex(a::TMSol, c::Colon) = a.xTM[:,c]
+getindex(a::TMSol, n::Integer, m::Integer) = a.xTM[m,n]
+getindex(a::TMSol, c::Colon, m::Integer) = a.xTM[m,c]
+getindex(a::TMSol, u::UnitRange, m::Integer) = a.xTM[m,u]

src/constructors.jl

@@ -146,3 +146,41 @@ TaylorModelN(a::T, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) wh
 @inline polynomial(tm::TaylorModelN) = tm.pol
 @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
+
+"""
+    TMSol{N,T,V1,V2,M}
+
+Structure containing the solution of a validated integration.
+
+# Fields
+    `time :: AbstractVector{T}`  Vector containing the expansion time of the `TaylorModel1` solutions
+
+    `fp   :: AbstractVector{IntervalBox{N,T}}`  IntervalBox vector representing the flowpipe
+
+    `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}},
+        M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
+    time :: V1
+    fp   :: V2
+    xTM  :: M
+
+    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)
+        return new{N,T,V1,V2,M}(time, fp, xTM)
+    end
+end
+
+@inline expansion_point(a::TMSol) = getfield(a,:time)
+@inline expansion_point(a::TMSol, n::Int) = getindex(getfield(a,:time),n)
+@inline flowpipe(a::TMSol) = getfield(a,:fp)
+@inline flowpipe(a::TMSol, n::Int) = getindex(getfield(a,:fp),n)
+@inline get_xTM(a::TMSol) = getfield(a,:xTM)
+@inline get_xTM(a::TMSol, n::Int) = getindex(getfield(a,:xTM),:,n)
+@inline domain(a::TMSol) = domain.(getindex(getfield(a, :xTM), 1, :)) # vector!
+@inline domain(a::TMSol, n::Int) = domain(getindex(getfield(a, :xTM), 1, n))
+@inline get_numvars(::TMSol{N,T,V1,V2,M}) where {N,T,V1,V2,M} = N

src/recipe.jl

@@ -44,3 +44,95 @@ end
 
     xs, ys
 end
+
+@recipe function g(sol::TMSol; vars=(0,1), timediv=1)
+    @assert length(vars) == length(unique(vars)) == 2 "`vars` must have 2 distinct integer indices"
+
+    return _plot_intvbox(sol, vars=vars, timediv=timediv)
+end
+
+function _plot_intvbox(sol::TMSol; vars=(0,1), timediv=1)
+    domT = mince_in_time(sol, var=0, timediv=timediv)
+
+    if 0 ∈ vars
+        tup0 = findfirst(0 ∈ vars)
+        var1 = vars[3-tup0]
+        v1 = _mince_in_time(sol, domT, var1, timediv)
+        if tup0 == 1
+            return @. IntervalBox(domT, v1)
+        else
+            return @. IntervalBox(v1, domT)
+        end
+    end
+
+    var1, var2 = vars
+    v1 = _mince_in_time(sol, domT, var1, timediv)
+    v2 = _mince_in_time(sol, domT, var2, timediv)
+    return @. IntervalBox(v1, v2)
+end
+
+"""
+    mince_in_time(sol::TMSol; var=0, timediv=1) --> ::Vector{Interval}
+
+For `var=0`, this function divides the time domain of each entry of `sol` in
+`timediv` parts (`timediv==1` is the initial domain), and returns the time
+intervals where the solution holds. This is useful for plotting or finding 
+specific events.
+For `var ≥ 1`, this function evaluates the flowpipes at the split domain 
+intervals, which is useful to decrease the overapproximations associated
+to the whole time domain.
+"""
+function mince_in_time(sol::TMSol; var::Int=0, timediv::Int=1)
+    @assert timediv > 0 "`timediv must be 1 or larger"
+    @assert 0 ≤ var ≤ get_numvars(sol)
+
+    domT = _mince_in_time(sol, Val(true), timediv)
+    var == 0 && return domT
+    return _mince_in_time(sol, domT, var, timediv)
+end
+
+# Mince in time var (var == 0)
+function _mince_in_time(sol::TMSol, ::Val{true}, timediv::Int=1)
+    # Case timediv == 1
+    if timediv == 1
+        return expansion_point(sol) .+ domain(sol)
+    end
+
+    # Case timediv > 1
+    domT = Vector{typeof(domain(sol,1))}(undef, timediv*length(sol))
+    i0 = 1
+    @inbounds for indT in eachindex(sol)
+        i1 = indT*timediv
+        domT[i0:i1] .= expansion_point(sol,indT) .+ mince(domain(sol,indT), timediv)
+        i0 = i1 + 1
+    end
+
+    return domT
+end
+
+# Mince other var (var > 0)
+function _mince_in_time(sol::TMSol, domT::Vector{Interval{T}}, var::Int, 
+        timediv::Int=1) where {T}
+    N = get_numvars(sol)
+    @assert 1 ≤ var ≤ N
+
+    # Case timediv == 1
+    if timediv == 1
+        return @. getindex(getfield(flowpipe(sol), :v), var)
+    end
+
+    # Case timediv > 1
+    vv = similar(domT)
+    normalized_box = symmetric_box(N, Float64)
+    δt = mince(domain(sol,1), timediv)
+
+    i0 = 1
+    @inbounds for indT in eachindex(sol)
+        i1 = indT*timediv
+        δt .= mince(domain(sol,indT), timediv)
+        @. vv[i0:i1] = evaluate(evaluate(sol[indT,var], δt), (normalized_box,))
+        i0 = i1 + 1
+    end
+
+    return vv    
+end

src/validatedODEs.jl

@@ -2,6 +2,7 @@
 
 const _DEF_MINABSTOL = 1.0e-50
 
+
 """
     remainder_taylorstep!(f!, t, x, dx, xI, dxI, δI, δt, params)
 
@@ -628,7 +629,7 @@ function validated_integ(f!, X0, t0::T, tmax::T, orderQ::Int, orderT::Int, absto
 
     end
 
-    return view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v, :, 1:nsteps)
+    return TMSol(view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v,:,1:nsteps))
 end
 
 """
@@ -904,5 +905,5 @@ function validated_integ2(f!, X0, t0::T, tf::T, orderQ::Int, orderT::Int,
         end
     end
 
-    return view(tv, 1:nsteps), view(xv, 1:nsteps), view(xTM1v, :, 1:nsteps)
+    return TMSol(view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v,:,1:nsteps))
 end