Merge pull request #5 from control-toolbox/temp

Temp

Author: joseph-gergaud

Project.toml

@@ -4,20 +4,30 @@ version = "0.1.0"
 authors = ["Olivier Cots <olivier.cots@toulouse-inp.fr>"]
 
 [deps]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
+ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+DiffEqDevTools = "2.48.0"
 DifferentiationInterface = "0.7.9"
+DualNumbers = "0.6.9"
 ForwardDiff = "1.2.2"
 LaTeXStrings = "1.4.0"
 LinearAlgebra = "1.12.0"
+ODEInterface = "0.5.0"
+ODEInterfaceDiffEq = "3.13.5"
+ODEProblemLibrary = "1.2.2"
 OrdinaryDiffEq = "6.102.1"
 Plots = "1.41.1"
 SciMLSensitivity = "7.90.0"

article/figures/chris.jl

@@ -0,0 +1,56 @@
+###############################################################################
+# From https://github.com/SciML/ODE.jl/blob/master/src/ODE.jl
+# estimator for initial step based on book
+# "Solving Ordinary Differential Equations I" by Hairer et al., p.169
+#
+#function hinit(F, x0, t0::T, tend, p, reltol, abstol) where T
+function hinit(F, x0, t0, tend, p, reltol, abstol)
+    # Returns first step, direction of integration and F evaluated at t0
+    tdir = sign(tend-t0)
+    tdir==0 && error("Zero time span")
+    tau = max(reltol*norm(x0, Inf), abstol)
+    d0 = norm(x0, Inf)/tau
+    f0 = F(t0, x0)
+    d1 = norm(f0, Inf)/tau
+    if d0 < 1e-5 || d1 < 1e-5
+        h0 = 1e-6
+    else
+        h0 = 0.01*(d0/d1)
+    end
+    # perform Euler step
+    x1 = x0 + tdir*h0*f0
+    f1 = F(t0 + tdir*h0, x1)
+    # estimate second derivative
+    d2 = norm(f1 - f0, Inf)/(tau*h0)
+    if max(d1, d2) <= 1e-15
+        h1 = max((10)^(-6), (10)^(-3)*h0)
+    else
+        pow = -(2 + log10(max(d1, d2)))/(p + 1)
+        h1 = 10^pow
+    end
+
+    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
+
+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.;
+p = 2
+x0 = [0.,1.,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]]
+hinit(fun2, xX0, t0, tend, p, reltol, abstol)
+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)

article/figures/test-dopri-DP5-julia.jl

@@ -0,0 +1,48 @@
+using OrdinaryDiffEq,  Test, DiffEqDevTools,
+      ODEInterface, ODEInterfaceDiffEq
+
+import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear
+
+prob = prob_ode_linear
+#=
+prob.tspan
+prob.f
+prob.u0
+prob.p
+prob.kwargs
+prob.problem_type
+=#
+
+sol = solve(prob, DP5())
+
+sol2 = solve(prob, dopri5())
+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)))
+
+# Change the norm due to error in dopri5.f
+sol2 = solve(prob, dopri5())
+println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2])
+@test sol.t[2] ≈ sol2.t[2]
+
+prob = deepcopy(prob_ode_linear)
+prob
+#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())
+
+sol2 = solve(prob2, dopri5())
+println("sol.t[2] - sol2.t[2] = ", sol.t[2] - sol2.t[2])
+@test 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)))
+
+# 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]

article/figures/test_myode43.jl

@@ -0,0 +1,56 @@
+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₀(λ)
+# (∂x(t,t0,x_0(λ),λ)/∂x0) = diag(exp(λ₁t),exp(λ₂t),exp((λ₁-λ₂)t))
+# ∂x(t,t0,x_0(λ),λ)/∂λ = [λ₂texp(λ₁t) exp(λ₁t) ; 0 texp(λ₂t) ; texp((λ₁-λ₂)t)) -texp((λ₁-λ₂)t))]
+# 
+# ∂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];  
+x0 = funx0(λ)
+n = length(x0)
+reltol = 1.e-4; abstol = 1.e-4
+myode43(fun1,x0,λ,tspan,reltol,abstol)
+
+xX0 = [x0 [0 1 ; 0 0 ; 0 0]]
+
+myode43(fun2,xX0,λ,tspan,reltol,abstol)
+
+RelTol = reltol*ones(n,p+1) 
+AbsTol = RelTol
+
+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)
+
+inith(fun1,x0,λ,t0,abstol,reltol)
+
+inith(fun2,xX0,λ,t0,abstol,reltol)
+RelTol = reltol*ones(n,p+1) 
+AbsTol = RelTol
+inith(fun2,xX0,λ,t0,AbsTol,RelTol)
+
+
+epsi = eps()
+epsi = 0
+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)
+

article/figures/test_zygote.jl

@@ -0,0 +1,93 @@
+using OrdinaryDiffEq
+using LinearAlgebra
+using Plots
+using LaTeXStrings
+using ForwardDiff: ForwardDiff
+using Zygote: Zygote
+#using Mooncake # plante à l'instalation
+using SciMLSensitivity
+using DifferentiationInterface
+using DataFrames
+
+include("../../src/CTDiffFlow.jl")
+using .CTDiffFlow
+include("../../src/myode43/myode43.jl")
+
+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)
+        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]; 
+N = 10
+t0 = 0. ; tf = 1.; tspan = (t0,tf);
+λ = [1.,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)]
+
+
+#dict_Zygote = Dict(:∂λ_sol => sol_∂λ_flow)
+#dict_ForwardDiff = Dict(:∂λ_sol => sol_∂λ_flow)
+
+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)
+#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(),λ)))
+# 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]
+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(),λ)))
+# Zygote with ∂λ_flow
+∂λ_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)))
+
+# ForwardDiff with ∂λ_flow
+∂λ_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)))
+
+function _flow(λ)
+    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(),λ)))
+
+jacobian(_flow,AutoForwardDiff(),λ)
+
+latexify(df,arraystyle=:pmatrix,env=:tabular)
+
+sol_Zygote = df[in.(df.AutoDiff, Ref(["solution","Zygote"])),:]
+sol_ForwardDiff = df[in.(df.AutoDiff, Ref(["solution", "ForwardDiff"])),:]

article/figures/time-steps-bruss.jl

@@ -0,0 +1,102 @@
+
+using Pkg
+Pkg.activate(".")
+using OrdinaryDiffEq
+using LinearAlgebra
+using Plots
+using LaTeXStrings
+using ForwardDiff: ForwardDiff
+using Zygote: Zygote
+#using Mooncake # plante à l'instalation
+using SciMLSensitivity
+using DifferentiationInterface
+
+include("../../src/CTDiffFlow.jl")
+using .CTDiffFlow
+#
+# Definition of the second member
+function bruss(x,par,t)
+  λ = par[1]
+  x₁ = x[1]
+  x₂ = x[2]
+  return [1+x₁^2*x₂-(λ+1)*x₁ , λ*x₁-x₁^2*x₂]
+end
+
+t0 = 0.; tf = 1.
+tspan = (t0,tf)
+λ = [3.]
+p = length(λ)
+x01 = 1.3
+funx0(λ) = [x01, λ[1]]
+tol = 1.e-4
+n = 2
+
+plt1 = plot(); plt2 = plot()
+algo = Tsit5()
+reltol = 1.e-4; abstol = 1.e-4
+
+# Flow
+ivp = ODEProblem(bruss, funx0(λ), (t0,tf), λ)
+sol_ivp = solve(ivp, alg=algo; reltol=reltol,abstol=abstol)
+println("Flow")
+println("Times for the initial flow T = ")
+println(sol_ivp.t)
+dict_T = Dict(:ivp => sol_ivp.t)
+dt = sol_ivp.t[2]
+
+# -------------------------------
+println("Vatiational equation")
+# ------------------------------
+
+∂λ_flow_var = CTDiffFlow.build_∂λ_flow_var(bruss,t0,funx0,tf,λ)
+
+#println(∂λ_flow_var(t0,funx0,tf,λ;reltol=reltol,abstol=abstol))
+sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol)
+sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u))
+println("Times default values= ", sol_var.t)
+dict_T[:var_reltol] = sol_var.t
+my_Inf = prevfloat(typemax(Float64))
+RelTol = [reltol*ones(n,1) Inf*ones(n,1)]/sqrt(p+1)
+AbsTol = [abstol*ones(n,1) Inf*ones(n,1)]/sqrt(p+1)
+sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol)
+sol_var = ∂λ_flow_var((t0,tf),funx0,λ;reltol=RelTol,abstol=AbsTol, dt=dt)#,internalnorm = (u,t)->norm(u) )
+println("RelTol = ", RelTol)
+println("Times = ", sol_var.t)
+println(sol_ivp.t - sol_var.t)
+dict_T[:var_RelTol] = sol_var.t
+
+# ne fonctionne pas 
+#@btime ∂λ_flow_var(t0,funx0,tf,λ;reltol=reltol,abstol=abstol)
+
+println("Automatic differentiation")
+println("with ForwardDiff")
+∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ)
+sol_diff_auto_flow, T = ∂λ_flow(t0,funx0,tf,λ;reltol=RelTol,abstol=AbsTol,print_times=true)
+println("Times : ", T)
+dict_T[:diff_auto_ForwardDiff] = T
+#sol_diff_auto_flow, T = ∂λ_flow(tspan,funx0,λ;reltol=RelTol,abstol=AbsTol,print_times=true)
+#= 
+println("internalnorm = (u,t)->norm(u)")
+println(∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,internalnorm = (u,t)->norm(u),print_times=true))
+#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol)
+=#
+println("With Zygote")
+# with Zygote
+∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ; backend=AutoZygote())
+sol_diff_auto_flow_tf, T = ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,print_times=true)
+dict_T[:diff_auto_Zygote] = T
+println(sol_diff_auto_flow_tf-sol_var)
+#sol_diff_auto_flow,T = ∂λ_flow(tspan,funx0,λ;reltol=reltol,abstol=abstol,print_times=true)
+#println(sol_diff_auto_flow)
+#reshape(sol_diff_auto_flow,2,7)
+
+#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol)
+
+println("With Mooncake")
+# plante
+#∂λ_flow = CTDiffFlow.build_∂λ_flow(bruss,t0,funx0,tf,λ; backend=AutoMooncake())
+#println(∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol,print_times=true))
+#@btime ∂λ_flow(t0,funx0,tf,λ;reltol=reltol,abstol=abstol)
+
+#plot(plt1,plt2)
+