Skip to content

Commit

Permalink
Fix dense output (#168)
Browse files Browse the repository at this point in the history 
* WIP: Fix dense ouput and update tests

* Bump minor version

* Add breaking change warning

* Rename variable

* Small fix

* Fix tests in nightly (1.10.0-DEV)

* Fix tests
  • Loading branch information
@PerezHz
PerezHz committed Jun 29, 2023
1 parent f2e0a47 commit 6086c17
Show file tree
Hide file tree
Showing 7 changed files with 53 additions and 54 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
@@ -1,7 +1,7 @@
name = "TaylorIntegration"
uuid = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
repo = "https://github.com/PerezHz/TaylorIntegration.jl.git"
version = "0.13.0"
version = "0.14.0"

[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Expand Down
11 changes: 11 additions & 0 deletions src/TaylorIntegration.jl
View file
Expand Up @@ -20,6 +20,17 @@ include("rootfinding.jl")
include("common.jl")

function __init__()

### Breaking change warning added in v0.14.0; to be deleted in next minor version (v0.15.0)
@warn("""\n\n
# Breaking changes
When dense output from `taylorinteg` is requested by the user,
the first row of the dense output polynomials matrix now represents
the Taylor expansion of the solution around the initial condition
w.r.t. the independent variable in the ODE. Previously it represented
the initial condition plus a Taylor1 perturbation.\n
""")

@static if !isdefined(Base, :get_extension)
@require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
include("../ext/TaylorIntegrationDiffEq.jl")
Expand Down
36 changes: 12 additions & 24 deletions src/explicitode.jl
View file
Expand Up @@ -364,17 +364,14 @@ for V in (:(Val{true}), :(Val{false}))
tv = Array{T}(undef, maxsteps+1)
xv = Array{U}(undef, maxsteps+1)
if $V == Val{true}
polynV = Array{Taylor1{U}}(undef, maxsteps+1)
psol = Array{Taylor1{U}}(undef, maxsteps)
end

# Initial conditions
nsteps = 1
@inbounds t[0] = t0
@inbounds tv[1] = t0
@inbounds xv[1] = x0
if $V == Val{true}
@inbounds polynV[1] = deepcopy(x)
end
sign_tstep = copysign(1, tmax-t0)

# Integration
Expand All @@ -384,7 +381,7 @@ for V in (:(Val{true}), :(Val{false}))
δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
x0 = evaluate(x, δt) # new initial condition
if $V == Val{true}
@inbounds polynV[nsteps+1] = deepcopy(x)
@inbounds psol[nsteps] = deepcopy(x)
end
@inbounds x[0] = x0
t0 += δt
Expand All @@ -401,7 +398,7 @@ for V in (:(Val{true}), :(Val{false}))
end

if $V == Val{true}
return view(tv,1:nsteps), view(xv,1:nsteps), view(polynV, 1:nsteps)
return view(tv,1:nsteps), view(xv,1:nsteps), view(psol, 1:nsteps-1)
elseif $V == Val{false}
return view(tv,1:nsteps), view(xv,1:nsteps)
end
Expand All @@ -414,17 +411,14 @@ for V in (:(Val{true}), :(Val{false}))
tv = Array{T}(undef, maxsteps+1)
xv = Array{U}(undef, maxsteps+1)
if $V == Val{true}
polynV = Array{Taylor1{U}}(undef, maxsteps+1)
psol = Array{Taylor1{U}}(undef, maxsteps)
end

# Initial conditions
nsteps = 1
@inbounds t[0] = t0
@inbounds tv[1] = t0
@inbounds xv[1] = x0
if $V == Val{true}
@inbounds polynV[1] = deepcopy(x)
end
sign_tstep = copysign(1, tmax-t0)

# Integration
Expand All @@ -434,7 +428,7 @@ for V in (:(Val{true}), :(Val{false}))
δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
x0 = evaluate(x, δt) # new initial condition
if $V == Val{true}
@inbounds polynV[nsteps+1] = deepcopy(x)
@inbounds psol[nsteps] = deepcopy(x)
end
@inbounds x[0] = x0
t0 += δt
Expand All @@ -451,7 +445,7 @@ for V in (:(Val{true}), :(Val{false}))
end

if $V == Val{true}
return view(tv,1:nsteps), view(xv,1:nsteps), view(polynV, 1:nsteps)
return view(tv,1:nsteps), view(xv,1:nsteps), view(psol, 1:nsteps-1)
elseif $V == Val{false}
return view(tv,1:nsteps), view(xv,1:nsteps)
end
Expand Down Expand Up @@ -497,7 +491,7 @@ for V in (:(Val{true}), :(Val{false}))
tv = Array{T}(undef, maxsteps+1)
xv = Array{U}(undef, dof, maxsteps+1)
if $V == Val{true}
polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1)
psol = Array{Taylor1{U}}(undef, dof, maxsteps)
end
xaux = Array{Taylor1{U}}(undef, dof)

Expand All @@ -507,9 +501,6 @@ for V in (:(Val{true}), :(Val{false}))
x0 = deepcopy(q0)
@inbounds tv[1] = t0
@inbounds xv[:,1] .= q0
if $V == Val{true}
@inbounds polynV[:,1] .= deepcopy.(x)
end
sign_tstep = copysign(1, tmax-t0)

# Integration
Expand All @@ -521,7 +512,7 @@ for V in (:(Val{true}), :(Val{false}))
evaluate!(x, δt, x0) # new initial condition
if $V == Val{true}
# Store the Taylor polynomial solution
@inbounds polynV[:,nsteps+1] .= deepcopy.(x)
@inbounds psol[:,nsteps] .= deepcopy.(x)
end
@inbounds for i in eachindex(x0)
x[i][0] = x0[i]
Expand All @@ -542,7 +533,7 @@ for V in (:(Val{true}), :(Val{false}))

if $V == Val{true}
return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:),
view(transpose(view(polynV, :, 1:nsteps)), 1:nsteps, :)
view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :)
elseif $V == Val{false}
return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:)
end
Expand All @@ -559,17 +550,14 @@ for V in (:(Val{true}), :(Val{false}))
tv = Array{T}(undef, maxsteps+1)
xv = Array{U}(undef, dof, maxsteps+1)
if $V == Val{true}
polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1)
psol = Array{Taylor1{U}}(undef, dof, maxsteps)
end

# Initial conditions
@inbounds t[0] = t0
x0 = deepcopy(q0)
@inbounds tv[1] = t0
@inbounds xv[:,1] .= q0
if $V == Val{true}
@inbounds polynV[:,1] .= deepcopy.(x)
end
sign_tstep = copysign(1, tmax-t0)

# Integration
Expand All @@ -581,7 +569,7 @@ for V in (:(Val{true}), :(Val{false}))
evaluate!(x, δt, x0) # new initial condition
if $V == Val{true}
# Store the Taylor polynomial solution
@inbounds polynV[:,nsteps+1] .= deepcopy.(x)
@inbounds psol[:,nsteps] .= deepcopy.(x)
end
@inbounds for i in eachindex(x0)
x[i][0] = x0[i]
Expand All @@ -602,7 +590,7 @@ for V in (:(Val{true}), :(Val{false}))

if $V == Val{true}
return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:),
view(transpose(view(polynV, :, 1:nsteps)), 1:nsteps, :)
view(transpose(view(psol, :, 1:nsteps-1)), 1:nsteps-1, :)
elseif $V == Val{false}
return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:)
end
Expand Down
16 changes: 8 additions & 8 deletions src/rootfinding.jl
View file
Expand Up @@ -180,11 +180,11 @@ function taylorinteg(f!, g, q0::Array{U,1}, t0::T, tmax::T,
# Re-initialize the Taylor1 expansions
t = t0 + Taylor1( T, order )
x .= Taylor1.( q0, order )
return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, rv, params,
maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol, rv, params;
maxsteps, eventorder, newtoniter, nrabstol)
else
return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol,params,
maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
return _taylorinteg!(f!, g, t, x, dx, q0, t0, tmax, abstol,params;
maxsteps, eventorder, newtoniter, nrabstol)
end
end

Expand Down Expand Up @@ -355,11 +355,11 @@ function taylorinteg(f!, g, q0::Array{U,1}, trange::AbstractVector{T},
# Re-initialize the Taylor1 expansions
t = t0 + Taylor1( T, order )
x .= Taylor1.(q0, order)
return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, params,
maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol, rv, params;
maxsteps, eventorder, newtoniter, nrabstol)
else
return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol,params,
maxsteps=maxsteps, eventorder=eventorder, newtoniter=newtoniter, nrabstol=nrabstol)
return _taylorinteg!(f!, g, t, x, dx, q0, trange, abstol,params;
maxsteps, eventorder, newtoniter, nrabstol)
end
end

Expand Down
20 changes: 10 additions & 10 deletions test/many_ode.jl
View file
Expand Up @@ -68,13 +68,13 @@ import Logging: Warn
@test abs(xv2[2,1] - 4.8) ≤ eps(4.8)

# Output includes Taylor polynomial solution
tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg(
tv, xv, psol = (@test_logs (Warn, max_iters_reached()) taylorinteg(
eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2))
@test size(polynV) == (3, 2)
@test size(psol) == (2, 2)
@test xv[1,:] == q0
@test polynV[1,:] == Taylor1.(q0, _order)
@test xv[2,:] == evaluate.(polynV[2, :], tv[2]-tv[1])
@test polynV[2,2] == Taylor1(ones(_order+1))
@test xv[2,:] == evaluate.(psol[1, :], tv[2]-tv[1])
@test psol[1,1] == Taylor1([3.0^i for i in 1:_order+1])
@test psol[1,2] == Taylor1(ones(_order+1))
end

@testset "Tests: dot{x}=x.^2, x(0) = [3, 3]" begin
Expand Down Expand Up @@ -120,11 +120,11 @@ import Logging: Warn
@test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 1.0e-11

# Output includes Taylor polynomial solution
tv, xv, polynV = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, Val(true), nothing)
@test size(polynV) == size(xv)
@test polynV[1,:] == Taylor1.(q0, _order)
@test xv[end,1] == evaluate.(polynV[end, 1], tv[end]-tv[end-1])
@test polynV[:,1] == polynV[:,2]
tv, xv, psol = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, Val(true), nothing)
@test size(psol) == ( size(xv, 1)-1, size(xv, 2) )
@test psol[1,:] == fill(Taylor1([3.0^i for i in 1:_order+1]), 2)
@test xv[end,1] == evaluate.(psol[end, 1], tv[end]-tv[end-1])
@test psol[:,1] == psol[:,2]
end

@testset "Test non-autonomous ODE (2): dot{x}=cos(t)" begin
Expand Down
9 changes: 4 additions & 5 deletions test/one_ode.jl
View file
Expand Up @@ -81,13 +81,12 @@ import Logging: Warn
@test abs(xv2[5] - 2.0) ≤ eps(2.0)

# Output includes Taylor polynomial solution
tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg(
tv, xv, psol = (@test_logs (Warn, max_iters_reached()) taylorinteg(
eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2))
@test length(polynV) == 3
@test length(psol) == 2
@test xv[1] == x0
@test polynV[1] == Taylor1(x0, _order)
@test xv[2] == evaluate(polynV[2], tv[2]-tv[1])
@test polynV[2] == Taylor1(ones(_order+1))
@test psol[1] == Taylor1(ones(_order+1))
@test xv[2] == evaluate(psol[1], tv[2]-tv[1])
end

@testset "Tests: dot{x}=x^2, x(0) = 3; nsteps <= maxsteps" begin
Expand Down
13 changes: 7 additions & 6 deletions test/taylorize.jl
View file
@@ -1,6 +1,6 @@
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration, DiffEqBase
using TaylorIntegration, OrdinaryDiffEq
using Test
using LinearAlgebra: norm
using InteractiveUtils: methodswith
Expand Down Expand Up @@ -134,11 +134,10 @@ import Logging: Warn
# Output includes Taylor polynomial solution
tv3t, xv3t, polynV3t = (@test_logs (Warn, max_iters_reached()) taylorinteg(
xdot3, x0, t0, tf, _order, _abstol, Val(true), 0.0, maxsteps=2))
@test length(polynV3t) == 3
@test length(polynV3t) == 2
@test xv3t[1] == x0
@test polynV3t[1] == Taylor1(x0, _order)
@test xv3t[2] == evaluate(polynV3t[2], tv3t[2]-tv3t[1])
@test polynV3t[2] == Taylor1([(-1.0)^i for i=0:_order])
@test polynV3t[1] == Taylor1([(-1.0)^i for i=0:_order])
@test xv3t[2] == evaluate(polynV3t[1], tv3t[2]-tv3t[1])
end


Expand Down Expand Up @@ -1280,7 +1279,7 @@ import Logging: Warn
[Array{Taylor1{_S},4}(undef, 0, 0, 0, 0)]))

# Issue 96: deal with `elseif`s, `continue` and `break`
@test newex1.args[2].args[6].args[2].args[3] == :(
ex = :(
for i = 1:length(q)
if i == 1
TaylorSeries.mul!(dq[i], 2, q[i], ord)
Expand All @@ -1299,6 +1298,8 @@ import Logging: Warn
end
end)

@test newex1.args[2].args[6].args[2].args[3] == Base.remove_linenums!(ex)

# Throws no error
ex = :(
function err_arr_indx!(Dz, z, p, t)
Expand Down

2 comments on commit 6086c17

@PerezHz
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/86536

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.14.0 -m "<description of version>" 6086c170565681afb4d2ec761fc712749cebd91c
git push origin v0.14.0

Please sign in to comment.