Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 16 additions & 11 deletions article/figures/chris.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,25 +32,30 @@ function hinit(F, x0, t0, tend, p, reltol, abstol)
return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0
end

A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]]
fun1(t, x) = A(λ)*x

A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]]
fun1(t,x) = A(λ)*x
function fun2(t, xX)
[fun1(xX[:, 1], λ, t) fun1(
xX[:, 2:end], λ, t
)+[xX[1, 1] 0; 0 xX[2, 1]; xX[3, 1] xX[3, 1]]]
end

fun2(t,xX) = [fun1(xX[:,1],λ,t) fun1(xX[:,2:end],λ,t)+[xX[1,1] 0 ; 0 xX[2,1] ; xX[3,1] xX[3,1]] ]

t0 = 0. ; tend = 2.;
t0 = 0.0;
tend = 2.0;
p = 2
x0 = [0.,1.,1]
reltol = 1.e-4; abstol = 1.e-4
x0 = [0.0, 1.0, 1]
reltol = 1.e-4;
abstol = 1.e-4

hinit(fun1, x0, t0, tend, p, reltol, abstol)

xX0 = [funx0(λ) [0 1 ; 0 0 ; 0 0]]
xX0 = [funx0(λ) [0 1; 0 0; 0 0]]
hinit(fun2, xX0, t0, tend, p, reltol, abstol)
RelTol = reltol*ones(n,p+1)
RelTol = reltol*ones(n, p+1)
AbsTol = RelTol

hinit(fun2, xX0, t0, tend, p, RelTol, AbsTol)

RelTol = [reltol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1)
AbsTol = [abstol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1)
RelTol = [reltol*ones(n, 1) Inf*ones(n, p)]/sqrt(p+1)
AbsTol = [abstol*ones(n, 1) Inf*ones(n, p)]/sqrt(p+1)
9 changes: 4 additions & 5 deletions article/figures/test-dopri-DP5-julia.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
using OrdinaryDiffEq, Test, DiffEqDevTools,
ODEInterface, ODEInterfaceDiffEq
using OrdinaryDiffEq, Test, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq

import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

Expand All @@ -20,7 +19,7 @@ println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2])
@test sol.t[2] ≈ sol2.t[2]

prob = prob_ode_2Dlinear
sol = solve(prob, DP5(), internalnorm = (u, t) -> sqrt(sum(abs2, u)))
sol = solve(prob, DP5(); internalnorm=(u, t) -> sqrt(sum(abs2, u)))

# Change the norm due to error in dopri5.f
sol2 = solve(prob, dopri5())
Expand All @@ -40,9 +39,9 @@ println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2])
prob = deepcopy(prob_ode_2Dlinear)
#prob2 = ODEProblem(prob.f, prob.u0, (1.0, 0.0), 1.01)
prob2 = ODEProblem(prob.f, prob.u0, (0.0, 1.0), 1.01)
sol = solve(prob2, DP5(), internalnorm = (u, t) -> norm(u)/sqrt(length(u)))#sqrt(sum(abs2, u)))
sol = solve(prob2, DP5(); internalnorm=(u, t) -> norm(u)/sqrt(length(u)))#sqrt(sum(abs2, u)))

# Change the norm due to error in dopri5.f
sol2 = solve(prob2, dopri5())
println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2])
@test sol.t[2] ≈ sol2.t[2]
@test sol.t[2] ≈ sol2.t[2]
63 changes: 33 additions & 30 deletions article/figures/test_myode43.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
include("../../src/myode43/myode43.jl")


# ẋ = [λ₁x₁ ; λ₂x₂ ; (λ₁-λ₂)x₃
# x(t0) = x₀(λ) = [λ₂, 1, 1]
# x(t,t0,x_0(λ),λ) = diag(exp(λ₁t),exp(λ₂t),exp((λ₁-λ₂)t))x₀(λ)
Expand All @@ -9,48 +8,52 @@ include("../../src/myode43/myode43.jl")
#
# ∂x0_flow

A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]]
fun1(x,λ,t) = A(λ)*x

fun2(xX,λ,t) = [fun1(xX[:,1],λ,t) fun1(xX[:,2:end],λ,t)+[xX[1,1] 0 ; 0 xX[2,1] ; xX[3,1] xX[3,1]] ]



t0 = 0. ; tf = 2.; tspan = (t0,tf);
λ = [1.,2]; p = length(λ);
funx0(λ) = [λ[2],1.,1];
A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]]
fun1(x, λ, t) = A(λ)*x

function fun2(xX, λ, t)
[fun1(xX[:, 1], λ, t) fun1(
xX[:, 2:end], λ, t
)+[xX[1, 1] 0; 0 xX[2, 1]; xX[3, 1] xX[3, 1]]]
end

t0 = 0.0;
tf = 2.0;
tspan = (t0, tf);
λ = [1.0, 2];
p = length(λ);
funx0(λ) = [λ[2], 1.0, 1];
x0 = funx0(λ)
n = length(x0)
reltol = 1.e-4; abstol = 1.e-4
myode43(fun1,x0,λ,tspan,reltol,abstol)
reltol = 1.e-4;
abstol = 1.e-4
myode43(fun1, x0, λ, tspan, reltol, abstol)

xX0 = [x0 [0 1 ; 0 0 ; 0 0]]
xX0 = [x0 [0 1; 0 0; 0 0]]

myode43(fun2,xX0,λ,tspan,reltol,abstol)
myode43(fun2, xX0, λ, tspan, reltol, abstol)

RelTol = reltol*ones(n,p+1)
RelTol = reltol*ones(n, p+1)
AbsTol = RelTol

myode43(fun2,xX0,λ,tspan,RelTol,AbsTol)
myode43(fun2, xX0, λ, tspan, RelTol, AbsTol)

RelTol = [reltol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1)
AbsTol = [abstol*ones(n,1) Inf*ones(n,p)]/sqrt(p+1)
RelTol = [reltol*ones(n, 1) Inf*ones(n, p)]/sqrt(p+1)
AbsTol = [abstol*ones(n, 1) Inf*ones(n, p)]/sqrt(p+1)

inith(fun1,x0,λ,t0,abstol,reltol)
inith(fun1, x0, λ, t0, abstol, reltol)

inith(fun2,xX0,λ,t0,abstol,reltol)
RelTol = reltol*ones(n,p+1)
inith(fun2, xX0, λ, t0, abstol, reltol)
RelTol = reltol*ones(n, p+1)
AbsTol = RelTol
inith(fun2,xX0,λ,t0,AbsTol,RelTol)

inith(fun2, xX0, λ, t0, AbsTol, RelTol)

epsi = eps()
epsi = 0
xX0 = [x0 [epsi 1 ; epsi epsi ; epsi epsi]]
xX0 = [x0 [epsi 1; epsi epsi; epsi epsi]]
my_Inf = prevfloat(typemax(Float64))
#my_Inf = Inf
RelTol = [reltol*ones(n,1) my_Inf*ones(n,p)]/sqrt(p+1)
AbsTol = [abstol*ones(n,1) my_Inf*ones(n,p)]/sqrt(p+1)
inith(fun2,xX0,λ,t0,AbsTol,RelTol)
myode43(fun2,xX0,λ,tspan,RelTol,AbsTol)

RelTol = [reltol*ones(n, 1) my_Inf*ones(n, p)]/sqrt(p+1)
AbsTol = [abstol*ones(n, 1) my_Inf*ones(n, p)]/sqrt(p+1)
inith(fun2, xX0, λ, t0, AbsTol, RelTol)
myode43(fun2, xX0, λ, tspan, RelTol, AbsTol)
109 changes: 74 additions & 35 deletions article/figures/test_zygote.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,81 +13,120 @@ include("../../src/CTDiffFlow.jl")
using .CTDiffFlow
include("../../src/myode43/myode43.jl")

function my_euler(fun,t0,x0,tf,λ,N)
ti = t0; xi = x0
function my_euler(fun, t0, x0, tf, λ, N)
ti = t0;
xi = x0
h = (tf-t0)/N
for i in 1:N
xi += h*fun(xi,λ,ti)
xi += h*fun(xi, λ, ti)
ti = ti+h
end
return xi
end


# Definition of the edo
A(λ) = [λ[1] 0 0 ; 0 λ[2] 0 ; 0 0 λ[1]-λ[2]]
fun1(x,λ,t) = A(λ)*x
funx0(λ) = [λ[2],1.,1];
A(λ) = [λ[1] 0 0; 0 λ[2] 0; 0 0 λ[1]-λ[2]]
fun1(x, λ, t) = A(λ)*x
funx0(λ) = [λ[2], 1.0, 1];
N = 10
t0 = 0. ; tf = 1.; tspan = (t0,tf);
λ = [1.,2];
t0 = 0.0;
tf = 1.0;
tspan = (t0, tf);
λ = [1.0, 2];
reltol = 1.e-8
abstol = reltol
# ∂λ_flow(λ)
sol_∂λ_flow = [λ[2]*exp(λ[1]*tf) exp(λ[1]*tf)
0. tf*exp(λ[2]*tf)
tf*exp((λ[1]-λ[2])*tf) -tf*exp((λ[1]-λ[2])*tf)]

sol_∂λ_flow = [
λ[2]*exp(λ[1]*tf) exp(λ[1]*tf)
0.0 tf*exp(λ[2]*tf)
tf*exp((λ[1]-λ[2])*tf) -tf*exp((λ[1]-λ[2])*tf)
]

#dict_Zygote = Dict(:∂λ_sol => sol_∂λ_flow)
#dict_ForwardDiff = Dict(:∂λ_sol => sol_∂λ_flow)

df = DataFrame(AutoDiff = String[], algo=String[], Jacobian = Matrix{<:Real}[])
df = DataFrame(; AutoDiff=String[], algo=String[], Jacobian=Matrix{<:Real}[])
push!(df, ("solution", "", sol_∂λ_flow))

# Zygote with my_euler
_flow(λ) = my_euler(fun1,t0,funx0(λ),tf,λ,N)
_flow(λ) = my_euler(fun1, t0, funx0(λ), tf, λ, N)
#dict_Zygote[:my_euler] = jacobian(_flow,AutoZygote(),λ)
#dict_ForwardDiff[:my_euler] = jacobian(_flow,AutoForwardDiff(),λ)
push!(df, ("Zygote", "my-euler", jacobian(_flow,AutoZygote(),λ)))
push!(df, ("ForwardDiff", "my-euler", jacobian(_flow,AutoForwardDiff(),λ)))
push!(df, ("Zygote", "my-euler", jacobian(_flow, AutoZygote(), λ)))
push!(df, ("ForwardDiff", "my-euler", jacobian(_flow, AutoForwardDiff(), λ)))
# Zygote with numerical integration
function _flow_int(λ)
ivp = ODEProblem(fun1, funx0(λ), (t0,tf), λ)
#algo = get(ode_kwargs, :alg, Tsit5())
#println("algo = ", algo)

sol = solve(ivp, alg = Euler(),reltol=reltol,abstol=abstol,adaptive=false, dt = (tf-t0)/N)
return sol.u[end]
ivp = ODEProblem(fun1, funx0(λ), (t0, tf), λ)
#algo = get(ode_kwargs, :alg, Tsit5())
#println("algo = ", algo)

sol = solve(
ivp; alg=Euler(), reltol=reltol, abstol=abstol, adaptive=false, dt=(tf-t0)/N
)
return sol.u[end]
end

#dict_Zygote[:Euler] = jacobian(_flow_int,AutoZygote(),λ)
#dict_ForwardDiff[:Euler] = jacobian(_flow_int,AutoForwardDiff(),λ)
push!(df, ("Zygote", "Euler", jacobian(_flow_int,AutoZygote(),λ)))
push!(df, ("ForwardDiff", "Euler", jacobian(_flow_int,AutoForwardDiff(),λ)))
push!(df, ("Zygote", "Euler", jacobian(_flow_int, AutoZygote(), λ)))
push!(df, ("ForwardDiff", "Euler", jacobian(_flow_int, AutoForwardDiff(), λ)))
# Zygote with ∂λ_flow
∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ; backend=AutoZygote())
∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1, t0, funx0, tf, λ; backend=AutoZygote())
#dict_Zygote[:CTDiffFlowZygoteEuler] = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)#,print_times=false)
push!(df, ("Zygote", "CTDiffFlow-Euler", ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)))
push!(
df,
(
"Zygote",
"CTDiffFlow-Euler",
∂λ_flow(
t0,
funx0,
tf,
λ;
reltol=reltol,
abstol=abstol,
alg=Euler(),
adaptive=false,
dt=(tf-t0)/N,
),
),
)

# ForwardDiff with ∂λ_flow
∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1,t0,funx0,tf,λ; backend=AutoForwardDiff())
∂λ_flow = CTDiffFlow.build_∂λ_flow(fun1, t0, funx0, tf, λ; backend=AutoForwardDiff())
#dict_ForwardDiff[:CTDiffFlowForwardDiffEuler] = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)#,print_times=false)
push!(df, ("ForwardDiff", "my-euler", ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol, alg = Euler(),adaptive=false, dt = (tf-t0)/N)))
push!(
df,
(
"ForwardDiff",
"my-euler",
∂λ_flow(
t0,
funx0,
tf,
λ;
reltol=reltol,
abstol=abstol,
alg=Euler(),
adaptive=false,
dt=(tf-t0)/N,
),
),
)

function _flow(λ)
T, X = myode43(fun1,funx0(λ),λ,(t0,tf),abstol,reltol)
T, X = myode43(fun1, funx0(λ), λ, (t0, tf), abstol, reltol)
return X[end]
end

#dict_Zygote[:myode43Zygote] = jacobian(_flow,AutoZygote(),λ)
#dict_ForwardDiff[:myode43ForwardDiff] = jacobian(_flow,AutoForwardDiff(),λ)
push!(df, ("Zygote", "my-ode43", [NaN;;]))
push!(df, ("ForwardDiff", "my-ode43", jacobian(_flow,AutoForwardDiff(),λ)))
push!(df, ("ForwardDiff", "my-ode43", jacobian(_flow, AutoForwardDiff(), λ)))

jacobian(_flow,AutoForwardDiff(),λ)
jacobian(_flow, AutoForwardDiff(), λ)

latexify(df,arraystyle=:pmatrix,env=:tabular)
latexify(df; arraystyle=:pmatrix, env=:tabular)

sol_Zygote = df[in.(df.AutoDiff, Ref(["solution","Zygote"])),:]
sol_ForwardDiff = df[in.(df.AutoDiff, Ref(["solution", "ForwardDiff"])),:]
sol_Zygote = df[in.(df.AutoDiff, Ref(["solution", "Zygote"])), :]
sol_ForwardDiff = df[in.(df.AutoDiff, Ref(["solution", "ForwardDiff"])), :]
Loading