🤖 Format .jl files

Author: ocots

article/figures/chris.jl

@@ -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)

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

@@ -1,5 +1,4 @@
-using OrdinaryDiffEq,  Test, DiffEqDevTools,
-      ODEInterface, ODEInterfaceDiffEq
+using OrdinaryDiffEq, Test, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq
 
 import ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear
 
@@ -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())
@@ -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]

article/figures/test_myode43.jl

@@ -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₀(λ)
@@ -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)

article/figures/test_zygote.jl

@@ -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"])), :]