SciML · ChrisRackauckas · Jun 23, 2022
diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style = "sciml"
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
@@ -0,0 +1,42 @@
+name: format-check
+
+on:
+  push:
+    branches:
+      - 'master'
+      - 'release-'
+    tags: '*'
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+
+      - uses: actions/checkout@v1
+      - name: Install JuliaFormatter and format
+        # This will use the latest version by default but you can set the version like so:
+        #
+        # julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))'
+        run: |
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
+          julia  -e 'using JuliaFormatter; format(".", verbose=true)'
+      - name: Format check
+        run: |
+          julia -e '
+          out = Cmd(`git diff --name-only`) |> read |> String
+          if out == ""
+              exit(0)
+          else
+              @error "Some files have not been formatted !!!"
+              write(stdout, out)
+              exit(1)
+          end'
diff --git a/examples/double_pendulum.jl b/examples/double_pendulum.jl
@@ -2,113 +2,122 @@ using DifferentialEquations, Plots
 
 #Reference: https://diego.assencio.com/?index=e5ac36fcb129ce95a61f8e8ce0572dbf
 
-
 #Solving the double pendulum with a traditional ODE method
 #==========================================================#
 function doublependulum(du, u, params, t)
     l1 = params[1]
     l2 = params[2]
     m1 = params[3]
     m2 = params[4]
-    g  = params[5]
+    g = params[5]
 
     p1 = u[1]
     p2 = u[2]
     θ1 = u[3]
     θ2 = u[4]
 
-    h1 = p1*p2*sin(θ1-θ2)/(l1*l2*(m1+m2*sin(θ1-θ2)^2))
-    h2 = (m2*l2^2*p1^2 + (m1 + m2) * l1^2 * p2^2 - 2*m2*l1*l2*p1*p2*cos(θ1-θ2))/(2 * l1^2 * l2^2 *(m1 + m2*sin(θ1-θ2)^2)^2)
-
-    d_p1 = -(m1 + m2)* g*l1*sin(θ1) - h1 + h2*sin(2*(θ1-θ2))
-    d_θ1 = (l2*p1 - l1*p2*cos(θ1-θ2))/(l1^2 * l2*(m1 + m2*sin(θ1 - θ2)^2))
-    d_p2 = -m2*g*l2*sin(θ2) + h1 - h2*sin(2*(θ1-θ2))
-    d_θ2 = (-m2*l2*p1*cos(θ1 - θ2) + (m1 + m2)*l1*p2) / (m2*l1*l2^2*(m1 + m2*sin(θ1-θ2)^2))
+    h1 = p1 * p2 * sin(θ1 - θ2) / (l1 * l2 * (m1 + m2 * sin(θ1 - θ2)^2))
+    h2 = (m2 * l2^2 * p1^2 + (m1 + m2) * l1^2 * p2^2 -
+          2 * m2 * l1 * l2 * p1 * p2 * cos(θ1 - θ2)) /
+         (2 * l1^2 * l2^2 * (m1 + m2 * sin(θ1 - θ2)^2)^2)
 
+    d_p1 = -(m1 + m2) * g * l1 * sin(θ1) - h1 + h2 * sin(2 * (θ1 - θ2))
+    d_θ1 = (l2 * p1 - l1 * p2 * cos(θ1 - θ2)) / (l1^2 * l2 * (m1 + m2 * sin(θ1 - θ2)^2))
+    d_p2 = -m2 * g * l2 * sin(θ2) + h1 - h2 * sin(2 * (θ1 - θ2))
+    d_θ2 = (-m2 * l2 * p1 * cos(θ1 - θ2) + (m1 + m2) * l1 * p2) /
+           (m2 * l1 * l2^2 * (m1 + m2 * sin(θ1 - θ2)^2))
 
     du .= [d_p1, d_p2, d_θ1, d_θ2]
     return nothing
 end
 
-l1 = 1. #length of pendulum1
-l2 = 2. #length of pendulum2
-m1 = 1. #mass of pendulum1
-m2 = 1. #mass of pendulum2
+l1 = 1.0 #length of pendulum1
+l2 = 2.0 #length of pendulum2
+m1 = 1.0 #mass of pendulum1
+m2 = 1.0 #mass of pendulum2
 g = 9.81 #gravity
 
 params = [l1, l2, m1, m2, g]
-times = (0., 25.)
-u0 = [1.,1.,1.,1.]
+times = (0.0, 25.0)
+u0 = [1.0, 1.0, 1.0, 1.0]
 prob = ODEProblem(doublependulum, u0, times, params)
-sol1 = solve(prob, AutoVern7(Rodas5()), dt = .005)
+sol1 = solve(prob, AutoVern7(Rodas5()), dt = 0.005)
 
 #plot solution
-plot(sol1, vars=1, xlim=(0,20), label="Momentum1")
-plot!(sol1, vars=2, xlim=(0,20), label="Momentum2")
-plot!(sol1, vars=3, xlim=(0,20), label="theta1")
-plot!(sol1, vars=4, xlim=(0,20), label="theta2")
+plot(sol1, vars = 1, xlim = (0, 20), label = "Momentum1")
+plot!(sol1, vars = 2, xlim = (0, 20), label = "Momentum2")
+plot!(sol1, vars = 3, xlim = (0, 20), label = "theta1")
+plot!(sol1, vars = 4, xlim = (0, 20), label = "theta2")
 #==========================================================#
 
-
 #Now with HamiltonianProblem()
 #==========================================================#
 function H(p, θ, params)
     l1 = params[1]
     l2 = params[2]
     m1 = params[3]
     m2 = params[4]
-    g  = params[5]
+    g = params[5]
 
-    return  (m2*l2^2*p[1]^2 + (m1+m2)*l1^2*p[2]^2 - 2*m2*l1*l2*p[1]*p[2]*cos(θ[1]-θ[2]) ) /
-        (2*m2*l1^2*l2^2*(m1+m2*sin(θ[1]-θ[2])^2)) -
-        (m1+m2)*g*l1*cos(θ[1]) - m2*g*l2*cos(θ[2])
+    return (m2 * l2^2 * p[1]^2 + (m1 + m2) * l1^2 * p[2]^2 -
+            2 * m2 * l1 * l2 * p[1] * p[2] * cos(θ[1] - θ[2])) /
+           (2 * m2 * l1^2 * l2^2 * (m1 + m2 * sin(θ[1] - θ[2])^2)) -
+           (m1 + m2) * g * l1 * cos(θ[1]) - m2 * g * l2 * cos(θ[2])
 end
 
-l1 = 1. #length of pendulum1
-l2 = 2. #length of pendulum2
-m1 = 1. #mass of pendulum1
-m2 = 1. #mass of pendulum2
+l1 = 1.0 #length of pendulum1
+l2 = 2.0 #length of pendulum2
+m1 = 1.0 #mass of pendulum1
+m2 = 1.0 #mass of pendulum2
 g = 9.81 #gravity
 
 params = [l1, l2, m1, m2, g]
-q0 = [1.0,1.0]
-p0 = [1.0,1.0]
-times = (0.,25.)
+q0 = [1.0, 1.0]
+p0 = [1.0, 1.0]
+times = (0.0, 25.0)
 prob = HamiltonianProblem(H, q0, p0, times, params)
-sol2 = solve(prob, SofSpa10(), dt = .05)
+sol2 = solve(prob, SofSpa10(), dt = 0.05)
 
-plot(sol2, vars=1, xlim=(0,20), label="Momentum1")
-plot!(sol2, vars=2, xlim=(0,20), label="Momentum2")
-plot!(sol2, vars=3, xlim=(0,20), label="theta1")
-plot!(sol2, vars=4, xlim=(0,20), label="theta2")
+plot(sol2, vars = 1, xlim = (0, 20), label = "Momentum1")
+plot!(sol2, vars = 2, xlim = (0, 20), label = "Momentum2")
+plot!(sol2, vars = 3, xlim = (0, 20), label = "theta1")
+plot!(sol2, vars = 4, xlim = (0, 20), label = "theta2")
 #==========================================================#
 
-
 #Animation
 #==========================================================#
 function make_pretty_gif(sol)
     timepoints = sol.t
 
-    x1 = l1*sin.(sol[3,:])
-    y1 = -l1*cos.(sol[3,:])
-    x2 = x1 + l2*sin.(sol[4,:])
-    y2 = y1 - l2*cos.(sol[4,:])
+    x1 = l1 * sin.(sol[3, :])
+    y1 = -l1 * cos.(sol[3, :])
+    x2 = x1 + l2 * sin.(sol[4, :])
+    y2 = y1 - l2 * cos.(sol[4, :])
 
-    axis_lim = (l1+l2)*1.2
+    axis_lim = (l1 + l2) * 1.2
 
     anim = Animation()
-    for i =1:length(timepoints)
-        str = string("Time = ", round(timepoints[i],digits=1), " sec")
-        plot([0,x1[i]], [0,y1[i]], size=(400,300), xlim=(-axis_lim,axis_lim), ylim=(-axis_lim,1), markersize = 10, markershape = :circle,label ="",axis = [])
-        plot!([x1[i],x2[i]], [y1[i],y2[i]], markersize = 10, markershape = :circle,label ="",title = str, title_location = :left)
+    for i in 1:length(timepoints)
+        str = string("Time = ", round(timepoints[i], digits = 1), " sec")
+        plot([0, x1[i]], [0, y1[i]], size = (400, 300), xlim = (-axis_lim, axis_lim),
+             ylim = (-axis_lim, 1), markersize = 10, markershape = :circle, label = "",
+             axis = [])
+        plot!([x1[i], x2[i]], [y1[i], y2[i]], markersize = 10, markershape = :circle,
+              label = "", title = str, title_location = :left)
 
         if i > 8 #rainbow trail
-            plot!([x2[i-2:i]],   [y2[i-2:i]],  alpha = 0.15, linewidth = 2, color = :red, label=nothing)
-            plot!([x2[i-3:i-2]], [y2[i-3:i-2]],alpha = 0.15, linewidth = 2, color = :orange, label=nothing)
-            plot!([x2[i-4:i-3]], [y2[i-4:i-3]],alpha = 0.15, linewidth = 2, color = :yellow, label=nothing)
-            plot!([x2[i-6:i-4]], [y2[i-6:i-4]],alpha = 0.15, linewidth = 2, color = :green, label=nothing)
-            plot!([x2[i-7:i-6]], [y2[i-7:i-6]],alpha = 0.15, linewidth = 2, color = :blue, label=nothing)
-            plot!([x2[i-8:i-7]], [y2[i-8:i-7]],alpha = 0.15, linewidth = 2, color = :purple, label=nothing)
+            plot!([x2[(i - 2):i]], [y2[(i - 2):i]], alpha = 0.15, linewidth = 2,
+                  color = :red, label = nothing)
+            plot!([x2[(i - 3):(i - 2)]], [y2[(i - 3):(i - 2)]], alpha = 0.15, linewidth = 2,
+                  color = :orange, label = nothing)
+            plot!([x2[(i - 4):(i - 3)]], [y2[(i - 4):(i - 3)]], alpha = 0.15, linewidth = 2,
+                  color = :yellow, label = nothing)
+            plot!([x2[(i - 6):(i - 4)]], [y2[(i - 6):(i - 4)]], alpha = 0.15, linewidth = 2,
+                  color = :green, label = nothing)
+            plot!([x2[(i - 7):(i - 6)]], [y2[(i - 7):(i - 6)]], alpha = 0.15, linewidth = 2,
+                  color = :blue, label = nothing)
+            plot!([x2[(i - 8):(i - 7)]], [y2[(i - 8):(i - 7)]], alpha = 0.15, linewidth = 2,
+                  color = :purple, label = nothing)
         end
         frame(anim)
     end

diff --git a/examples/pendulum.jl b/examples/pendulum.jl
@@ -3,7 +3,7 @@ using DifferentialEquations, Plots
 #Solving the simple pendulum with a traditional ODE method
 #==========================================================#
 function pendulum(du, u, params, t)
-#reference: http://www.pgccphy.net/ref/advmech.pdf page 6
+    #reference: http://www.pgccphy.net/ref/advmech.pdf page 6
 
     g = params[1] #gravitational acceleration
     m = params[2] #mass
@@ -12,50 +12,47 @@ function pendulum(du, u, params, t)
     θ = u[1]
     ℒ = u[2]
 
-    d_θ = ℒ/(m*l^2)
-    d_ℒ = -m*g*l*sin(θ)
+    d_θ = ℒ / (m * l^2)
+    d_ℒ = -m * g * l * sin(θ)
 
     du .= [d_θ, d_ℒ]
     return nothing
 end
 
 g = 9.81
-m = 2.
-l = 1.
-params = [g,m,l]
-u0 = [1.0,1.0]
+m = 2.0
+l = 1.0
+params = [g, m, l]
+u0 = [1.0, 1.0]
 
-prob1 = ODEProblem(pendulum, u0, (0., 100.), params)
-sol1 = solve(prob1, AutoVern7(Rodas5()), dt = .05)
+prob1 = ODEProblem(pendulum, u0, (0.0, 100.0), params)
+sol1 = solve(prob1, AutoVern7(Rodas5()), dt = 0.05)
 
 #==========================================================#
 
-
-
 #Solving the simple pendulum with the DiffEqPhysics.jl HamiltonianProblem()
 #==========================================================#
 function H(ℒ, θ, params, t)
     g = params[1] #gravitational acceleration
     m = params[2] #mass
     l = params[3] #length
 
-    return ℒ^2/(2*m*l^2) + m*g*l*(1-cos(θ))
+    return ℒ^2 / (2 * m * l^2) + m * g * l * (1 - cos(θ))
 end
 
 g = 9.81
-m = 2.
-l = 1.
-params = [g,m,l]
-θ₀ = 1.
-ℒ₀ = 1.
-
-prob2 = HamiltonianProblem(H, ℒ₀, θ₀, (0., 100.), params)
-sol2 = solve(prob2, SofSpa10(), dt = .05);
+m = 2.0
+l = 1.0
+params = [g, m, l]
+θ₀ = 1.0
+ℒ₀ = 1.0
+
+prob2 = HamiltonianProblem(H, ℒ₀, θ₀, (0.0, 100.0), params)
+sol2 = solve(prob2, SofSpa10(), dt = 0.05);
 #==========================================================#
 
-
 #Plotting
 #==========================================================#
-plot(sol1, vars=1, tspan=(0,20), label="ODE Method")
-plot!(sol2.t, sol2[2,:], xlim=(0,20), label="Hamiltonian Method")
+plot(sol1, vars = 1, tspan = (0, 20), label = "ODE Method")
+plot!(sol2.t, sol2[2, :], xlim = (0, 20), label = "Hamiltonian Method")
 #They produce the same solution!
diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
@@ -25,9 +25,9 @@ AD (`ForwardDiff`).
 `H` may be defined with or without time as fourth argument. If both methods are defined,
 that with 4 arguments is used.
 """
-function HamiltonianProblem(H, p0::S, q0::T, tspan, param=nothing; kwargs...) where {S,T}
-  iip = T <: AbstractArray && !(T <: SArray) && S <: AbstractArray && !(S <: SArray)
-  HamiltonianProblem{iip}(H, p0, q0, tspan, param; kwargs...)
+function HamiltonianProblem(H, p0::S, q0::T, tspan, param = nothing; kwargs...) where {S, T}
+    iip = T <: AbstractArray && !(T <: SArray) && S <: AbstractArray && !(S <: SArray)
+    HamiltonianProblem{iip}(H, p0, q0, tspan, param; kwargs...)
 end
 
 struct PhysicsTag end
@@ -40,14 +40,18 @@ function generic_derivative(q0::Number, hami, x)
     ForwardDiff.derivative(hami, x)
 end
 
-function HamiltonianProblem{false}((dp, dq)::Tuple{Any,Any}, p0, q0, tspan, param=nothing; kwargs...)
-    return ODEProblem(DynamicalODEFunction{false}(dp, dq), ArrayPartition(p0, q0), tspan, param; kwargs...)
+function HamiltonianProblem{false}((dp, dq)::Tuple{Any, Any}, p0, q0, tspan,
+                                   param = nothing; kwargs...)
+    return ODEProblem(DynamicalODEFunction{false}(dp, dq), ArrayPartition(p0, q0), tspan,
+                      param; kwargs...)
 end
-function HamiltonianProblem{true}((dp, dq)::Tuple{Any,Any}, p0, q0, tspan, param=nothing; kwargs...)
-    return ODEProblem(DynamicalODEFunction{true}(dp, dq), ArrayPartition(p0, q0), tspan, param; kwargs...)
+function HamiltonianProblem{true}((dp, dq)::Tuple{Any, Any}, p0, q0, tspan, param = nothing;
+                                  kwargs...)
+    return ODEProblem(DynamicalODEFunction{true}(dp, dq), ArrayPartition(p0, q0), tspan,
+                      param; kwargs...)
 end
 
-function HamiltonianProblem{false}(H, p0, q0, tspan, param=nothing; kwargs...)
+function HamiltonianProblem{false}(H, p0, q0, tspan, param = nothing; kwargs...)
     if DiffEqBase.numargs(H) == 4
         dp = (p, q, param, t) -> generic_derivative(q0, q -> -H(p, q, param, t), q)
         dq = (p, q, param, t) -> generic_derivative(q0, p -> H(p, q, param, t), p)
@@ -59,18 +63,22 @@ function HamiltonianProblem{false}(H, p0, q0, tspan, param=nothing; kwargs...)
     return HamiltonianProblem{false}((dp, dq), p0, q0, tspan, param; kwargs...)
 end
 
-function HamiltonianProblem{true}(H, p0, q0, tspan, param=nothing; kwargs...)
+function HamiltonianProblem{true}(H, p0, q0, tspan, param = nothing; kwargs...)
     let cp = ForwardDiff.GradientConfig(PhysicsTag(), p0),
         cq = ForwardDiff.GradientConfig(PhysicsTag(), q0),
         vfalse = Val(false)
-        
+
         if DiffEqBase.numargs(H) == 4
-            dp = (Δp, p, q, param, t) -> ForwardDiff.gradient!(Δp, q->-H(p, q, param, t), q, cq, vfalse)
-            dq = (Δq, p, q, param, t) -> ForwardDiff.gradient!(Δq, p-> H(p, q, param, t), p, cp, vfalse)
+            dp = (Δp, p, q, param, t) -> ForwardDiff.gradient!(Δp, q -> -H(p, q, param, t),
+                                                               q, cq, vfalse)
+            dq = (Δq, p, q, param, t) -> ForwardDiff.gradient!(Δq, p -> H(p, q, param, t),
+                                                               p, cp, vfalse)
         else
             issue_depwarn()
-            dp = (Δp, p, q, param, t) -> ForwardDiff.gradient!(Δp, q->-H(p, q, param), q, cq, vfalse)
-            dq = (Δq, p, q, param, t) -> ForwardDiff.gradient!(Δq, p-> H(p, q, param), p, cp, vfalse)
+            dp = (Δp, p, q, param, t) -> ForwardDiff.gradient!(Δp, q -> -H(p, q, param), q,
+                                                               cq, vfalse)
+            dq = (Δq, p, q, param, t) -> ForwardDiff.gradient!(Δq, p -> H(p, q, param), p,
+                                                               cp, vfalse)
         end
         return HamiltonianProblem{true}((dp, dq), p0, q0, tspan, param; kwargs...)
     end