CI badges and formatting (#134)

.JuliaFormatter.toml

@@ -0,0 +1 @@
+style = "sciml"

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/codequality.yml

@@ -0,0 +1,10 @@
+name: Code Quality Check
+
+on: [pull_request]
+
+jobs:
+  code-style:
+    name: Format Suggestions
+    runs-on: ubuntu-latest
+    steps:
+      - uses: julia-actions/julia-format@v3

README.md

@@ -1,15 +1,19 @@
+# Differentiable programming for Differential equations: a review
+
+[![CI](https://github.com/ODINN-SciML/DiffEqSensitivity-Review/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/ODINN-SciML/DiffEqSensitivity-Review/actions/workflows/CI.yml)
 ![example workflow](https://github.com/ODINN-SciML/DiffEqSensitivity-Review/actions/workflows/latex.yml/badge.svg)
 ![example workflow](https://github.com/ODINN-SciML/DiffEqSensitivity-Review/actions/workflows/biblatex.yml/badge.svg)
-[![All Contributors](https://img.shields.io/github/all-contributors/ODINN-SciML/DiffEqSensitivity-Review?color=ee8449&style=flat-square)](#contributors)
 
-# Differentiable programming for Differential equations: a review
+[![All Contributors](https://img.shields.io/github/all-contributors/ODINN-SciML/DiffEqSensitivity-Review?color=ee8449&style=flat-square)](#contributors)
+[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 
 ### ⚠️ New preprint available! 📖 ⚠️
 
 The review paper is now available as a preprint on arXiv: https://arxiv.org/abs/2406.09699 
 
 If you want to cite this work, please use this BibTex citation:
-```
+
+```bibtex
 @misc{sapienza2024differentiable,
       title={Differentiable Programming for Differential Equations: A Review}, 
       author={Facundo Sapienza and Jordi Bolibar and Frank Schäfer and Brian Groenke and Avik Pal and Victor Boussange and Patrick Heimbach and Giles Hooker and Fernando Pérez and Per-Olof Persson and Christopher Rackauckas},

code/DirectMethods/Comparison/direct-comparision.jl

@@ -20,7 +20,7 @@ abstol = 1e-5
 function oscilatior!(du, u, p, t)
     ω = p[1]
     du[1] = u[2]
-    du[2] = - ω^2 * u[1]
+    du[2] = -ω^2 * u[1]
     nothing
 end
 
@@ -36,14 +36,14 @@ function solution_derivative(t, u0, p)
     ω = p[1]
     A₀ = u0[2] / ω
     B₀ = u0[1]
-    return A₀ * ( t * cos(ω * t) - sin(ω * t)/ω ) - B₀ * t * sin(ω * t)
+    return A₀ * (t * cos(ω * t) - sin(ω * t) / ω) - B₀ * t * sin(ω * t)
 end
 
 ######### Simple example of how to run the dynamcics ###########
 
 # Solve numerical problem
 prob = ODEProblem(oscilatior!, u0, tspan, p)
-sol = solve(prob, Tsit5(), reltol=reltol, abstol=abstol)
+sol = solve(prob, Tsit5(), reltol = reltol, abstol = abstol)
 
 u_final = sol.u[end][1]
 
@@ -65,17 +65,17 @@ function finitediff_solver(h, t, u0, p, reltol, abstol)
     tspan = (0.0, t)
     # Forward model with -h
     prob₋ = ODEProblem(oscilatior!, u0, tspan, p₋)
-    sol₋ = solve(prob₋, Tsit5(), reltol=reltol, abstol=abstol)
+    sol₋ = solve(prob₋, Tsit5(), reltol = reltol, abstol = abstol)
     # Forward model with +h
     prob₊ = ODEProblem(oscilatior!, u0, tspan, p₊)
-    sol₊ = solve(prob₊, Tsit5(), reltol=reltol, abstol=abstol)
+    sol₊ = solve(prob₊, Tsit5(), reltol = reltol, abstol = abstol)
 
-    return (sol₊.u[end][1] - sol₋.u[end][1]) /(2h)
+    return (sol₊.u[end][1] - sol₋.u[end][1]) / (2h)
 end
 
 ######### Simulation with differerent stepsizes ###########
 
-stepsizes = 2.0.^collect(round(log2(eps(Float64))):1:0)
+stepsizes = 2.0 .^ collect(round(log2(eps(Float64))):1:0)
 times = collect(t₀:1.0:t₁)
 
 # True derivative computend analytially
@@ -85,71 +85,85 @@ derivative_true = solution_derivative(t₁, u0, p)
 # derivative_numerical = finitediff_numerical.(stepsizes, Ref(t₁), Ref(u0), Ref(p))
 derivative_numerical = finitediff_numerical.(stepsizes, Ref(t₁), Ref(u0), Ref(p))
 derivative_finitediff_exact = finitediff_numerical.(stepsizes, Ref(t₁), Ref(u0), Ref(p))
-error_finitediff_exact = abs.((derivative_numerical .- derivative_true)./derivative_true)
+error_finitediff_exact = abs.((derivative_numerical .- derivative_true) ./ derivative_true)
 
 # Finite differences with solution from solver and low tolerance
-derivative_solver_low = finitediff_solver.(stepsizes, Ref(t₁), Ref(u0), Ref(p), Ref(1e-6), Ref(1e-6))
-error_finitediff_low = abs.((derivative_solver_low .- derivative_true)./derivative_true)
+derivative_solver_low = finitediff_solver.(
+    stepsizes, Ref(t₁), Ref(u0), Ref(p), Ref(1e-6), Ref(1e-6))
+error_finitediff_low = abs.((derivative_solver_low .- derivative_true) ./ derivative_true)
 
 # Finite differences with solution from solver and high tolerance
-derivative_solver_high = finitediff_solver.(stepsizes, Ref(t₁), Ref(u0), Ref(p), Ref(1e-12), Ref(1e-12))
-error_finitediff_high = abs.((derivative_solver_high .- derivative_true)./derivative_true)
+derivative_solver_high = finitediff_solver.(
+    stepsizes, Ref(t₁), Ref(u0), Ref(p), Ref(1e-12), Ref(1e-12))
+error_finitediff_high = abs.((derivative_solver_high .- derivative_true) ./ derivative_true)
 
 # Complex step differentiation with solution from solver and high tolerance
 u0_complex = ComplexF64.(u0)
 
-derivative_complex_low = complexstep_differentiation.(Ref(x -> solve(ODEProblem(oscilatior!, u0_complex, tspan, [x]), Tsit5(), reltol=1e-6, abstol=1e-6).u[end][1]), Ref(p[1]), stepsizes)
+derivative_complex_low = complexstep_differentiation.(
+    Ref(x -> solve(ODEProblem(oscilatior!, u0_complex, tspan, [x]), Tsit5(), reltol = 1e-6, abstol = 1e-6).u[end][1]),
+    Ref(p[1]),
+    stepsizes)
 error_complex_low = abs.((derivative_true .- derivative_complex_low) ./ derivative_true)
-derivative_complex_high = complexstep_differentiation.(Ref(x -> solve(ODEProblem(oscilatior!, u0_complex, tspan, [x]), Tsit5(), reltol=1e-12, abstol=1e-12).u[end][1]), Ref(p[1]), stepsizes)
+derivative_complex_high = complexstep_differentiation.(
+    Ref(x -> solve(ODEProblem(oscilatior!, u0_complex, tspan, [x]), Tsit5(), reltol = 1e-12, abstol = 1e-12).u[end][1]),
+    Ref(p[1]),
+    stepsizes)
 error_complex_high = abs.((derivative_true .- derivative_complex_high) ./ derivative_true)
 
 # Complex step Differentiation
-derivative_complex_exact = ComplexDiff.derivative.(ω -> solution(t₁, u0, [ω]), p[1], stepsizes)
-error_complex_exact = abs.((derivative_complex_exact .- derivative_true)./derivative_true)
+derivative_complex_exact = ComplexDiff.derivative.(
+    ω -> solution(t₁, u0, [ω]), p[1], stepsizes)
+error_complex_exact = abs.((derivative_complex_exact .- derivative_true) ./ derivative_true)
 
 # Forward AD applied to numerical solver
-derivative_AD_low = Zygote.gradient(p->solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol=1e-6, abstol=1e-6).u[end][1], p)[1][1]
+derivative_AD_low = Zygote.gradient(
+    p -> solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol = 1e-6, abstol = 1e-6).u[end][1],
+    p)[1][1]
 error_AD_low = abs((derivative_true - derivative_AD_low) / derivative_true)
 
-derivative_AD_high = Zygote.gradient(p->solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol=1e-12, abstol=1e-12).u[end][1], p)[1][1]
+derivative_AD_high = Zygote.gradient(
+    p -> solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol = 1e-12, abstol = 1e-12).u[end][1],
+    p)[1][1]
 error_AD_high = abs((derivative_true - derivative_AD_high) / derivative_true)
 
-
 ######### Figure ###########
 
-color_finitediff = RGBf(192/255, 57/255, 43/255)
-color_finitediff_low = RGBf(230/255, 126/255, 34/255)
-color_complex = RGBf(41/255, 128/255, 185/255)
-color_complex_low = RGBf(52/255, 152/255, 219/255)
-color_AD = RGBf(142/255, 68/255, 173/255)
-color_AD_low = RGBf(155/255, 89/255, 182/255)
+color_finitediff = RGBf(192 / 255, 57 / 255, 43 / 255)
+color_finitediff_low = RGBf(230 / 255, 126 / 255, 34 / 255)
+color_complex = RGBf(41 / 255, 128 / 255, 185 / 255)
+color_complex_low = RGBf(52 / 255, 152 / 255, 219 / 255)
+color_AD = RGBf(142 / 255, 68 / 255, 173 / 255)
+color_AD_low = RGBf(155 / 255, 89 / 255, 182 / 255)
 
-fig = Figure(resolution=(1000, 400))
-ax = Axis(fig[1, 1], xlabel = L"Stepsize ($\varepsilon$)", ylabel = L"\text{Absolute relative error}",
-          xscale = log10, yscale=log10, xlabelsize=24, ylabelsize=24, xticklabelsize=18, yticklabelsize=18)
+fig = Figure(resolution = (1000, 400))
+ax = Axis(fig[1, 1], xlabel = L"Stepsize ($\varepsilon$)",
+    ylabel = L"\text{Absolute relative error}", xscale = log10, yscale = log10,
+    xlabelsize = 24, ylabelsize = 24, xticklabelsize = 18, yticklabelsize = 18)
 
 # Plot derivatived of true solution (no numerical solver)
-lines!(ax, stepsizes, error_finitediff_exact, label=L"\text{Finite differences (exact solution)}",
-    color=color_finitediff, linewidth=2, linestyle = :dash)
+lines!(ax, stepsizes, error_finitediff_exact,
+    label = L"\text{Finite differences (exact solution)}",
+    color = color_finitediff, linewidth = 2, linestyle = :dash)
 lines!(ax, stepsizes, error_complex_exact,
-    label=L"\text{Complex step differentiation (exact solution)}",
-    color=color_complex, linewidth=2, linestyle = :dash)
-    
+    label = L"\text{Complex step differentiation (exact solution)}",
+    color = color_complex, linewidth = 2, linestyle = :dash)
+
 # Plot derivatives computed on top of numerical solver with finite differences
-scatter!(ax, stepsizes, error_finitediff_low,
-    label=L"Finite differences (tol=$10^{-6}$)", color=color_finitediff_low,
-    marker ='•', markersize=20)
-scatter!(ax, stepsizes, error_finitediff_high,
-    label=L"Finite differences (tol=$10^{-12}$)", color=color_finitediff,
-    marker ='•', markersize=30)
+scatter!(
+    ax, stepsizes, error_finitediff_low, label = L"Finite differences (tol=$10^{-6}$)",
+    color = color_finitediff_low, marker = '•', markersize = 20)
+scatter!(
+    ax, stepsizes, error_finitediff_high, label = L"Finite differences (tol=$10^{-12}$)",
+    color = color_finitediff, marker = '•', markersize = 30)
 
 # Plot derivatives computed on top of numerical solver with complex step method
 scatter!(ax, stepsizes, error_complex_low,
-    label=L"Complex step differentiation (tol=$10^{-6}$)",
-    color=color_complex_low, marker ='∘', markersize=20)
+    label = L"Complex step differentiation (tol=$10^{-6}$)",
+    color = color_complex_low, marker = '∘', markersize = 20)
 scatter!(ax, stepsizes, error_complex_high,
-    label=L"Complex step differentiation (tol=$10^{-12}$)", color=color_complex,
-    marker ='∘', markersize=30)
+    label = L"Complex step differentiation (tol=$10^{-12}$)",
+    color = color_complex, marker = '∘', markersize = 30)
 
 # AD
 # hlines!(ax, [error_AD_low, error_AD_high], color=color_AD, linewidth=1.5)
@@ -158,20 +172,19 @@ scatter!(ax, stepsizes, error_complex_high,
 # plot!(ax, [stepsizes[begin], stepsizes[end]],[error_AD_high, error_AD_high],
 #     color=color_AD, label=L"Forward AD (tol=$10^{-12}$)", marker='•', markersize=25)
 
-lines!(ax, stepsizes, repeat([error_AD_low], length(stepsizes)), 
-    color=color_AD_low, label=L"Forward AD (tol=$10^{-6}$)", linewidth=2)
-lines!(ax, stepsizes, repeat([error_AD_high], length(stepsizes)), 
-    color=color_AD, label=L"Forward AD (tol=$10^{-12}$)", linewidth=3)
+lines!(ax, stepsizes, repeat([error_AD_low], length(stepsizes)),
+    color = color_AD_low, label = L"Forward AD (tol=$10^{-6}$)", linewidth = 2)
+lines!(ax, stepsizes, repeat([error_AD_high], length(stepsizes)),
+    color = color_AD, label = L"Forward AD (tol=$10^{-12}$)", linewidth = 3)
 
 # Add legend
 fig[1, 2] = Legend(fig, ax)
 
 !ispath("Figures") && mkpath("Figures")
 save("Figures/DirectMethods_comparison.pdf", fig)
 
-
 ######### Benchmark ###########
 
 # It looks like complex step has better performance... both in speed and momory allocation.
 # @benchmark derivative_complex_low = complexstep_differentiation.(Ref(x -> solve(ODEProblem(oscilatior!, u0_complex, tspan, [x]), Tsit5(), reltol=1e-6, abstol=1e-6).u[end][1]), Ref(p[1]), [1e-5])
-# @benchmark derivative_AD_low = Zygote.gradient(p->solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol=1e-6, abstol=1e-6).u[end][1], p)[1][1]
+# @benchmark derivative_AD_low = Zygote.gradient(p->solve(ODEProblem(oscilatior!, u0, tspan, p), Tsit5(), reltol=1e-6, abstol=1e-6).u[end][1], p)[1][1]

code/DirectMethods/ComplexStep/complex_solver.jl

@@ -3,7 +3,7 @@ using OrdinaryDiffEq
 function dyn!(du::Array{Complex{Float64}}, u::Array{Complex{Float64}}, p, t)
     ω = p[1]
     du[1] = u[2]
-    du[2] = - ω^2 * u[1]
+    du[2] = -ω^2 * u[1]
 end
 
 tspan = [0.0, 10.0]
@@ -15,4 +15,5 @@ function complexstep_differentiation(f::Function, p::Float64, ε::Float64)
     return imag(f(p_complex)) / ε
 end
 
-complexstep_differentiation(x -> solve(ODEProblem(dyn!, u0, tspan, [x]), Tsit5()).u[end][1], 20., 1e-3)
+complexstep_differentiation(
+    x -> solve(ODEProblem(dyn!, u0, tspan, [x]), Tsit5()).u[end][1], 20.0, 1e-3)

code/DirectMethods/DualNumbers/dualnumber_definition.jl

@@ -11,38 +11,41 @@ end
 # Outer constructor
 function DualNumber(value::F) where {F <: AbstractFloat}
     DualNumber(value, F(0.0))
-end 
+end
 
 # Chain rules for binary opperators
 
 # Binary sum
-Base.:(+)(a::DualNumber, b::DualNumber) = DualNumber(value      = a.value + b.value,
-                                                     derivative = a.derivative + b.derivative)
+function Base.:(+)(a::DualNumber, b::DualNumber)
+    DualNumber(value = a.value + b.value, derivative = a.derivative + b.derivative)
+end
 
 # Binary product
-Base.:(*)(a::DualNumber, b::DualNumber) = DualNumber(value      = a.value * b.value,
-                                                     derivative = a.value*b.derivative + a.derivative*b.value)
+function Base.:(*)(a::DualNumber, b::DualNumber)
+    DualNumber(value = a.value * b.value,
+        derivative = a.value * b.derivative + a.derivative * b.value)
+end
 
 # Power
-Base.:(^)(a::DualNumber, b::AbstractFloat) = DualNumber(value      = a.value ^ b, 
-                                                        derivative = b * a.value^(b-1) * a.derivative)
-
+function Base.:(^)(a::DualNumber, b::AbstractFloat)
+    DualNumber(value = a.value^b, derivative = b * a.value^(b - 1) * a.derivative)
+end
 
 # Special functions
 
 function Base.:(sin)(a::DualNumber)
     value = sin(a.value)
     derivative = a.derivative * cos(a.value)
-    return DualNumber(value=value, derivative=derivative)
+    return DualNumber(value = value, derivative = derivative)
 end
 
 # Now we define a series of variables. We are interested in computing the derivative with respect to the variable "a":
 
-a = DualNumber(value=1.0, derivative=1.0)
+a = DualNumber(value = 1.0, derivative = 1.0)
 
-b = DualNumber(value=2.0, derivative=0.0)
-c = DualNumber(value=3.0, derivative=0.0)
+b = DualNumber(value = 2.0, derivative = 0.0)
+c = DualNumber(value = 3.0, derivative = 0.0)
 
 # Now, we can evaluate a new DualNumber
 result = a * b * c
-# println("The derivative of a*b*c with respect to a is: ", result.derivative)
+# println("The derivative of a*b*c with respect to a is: ", result.derivative)

code/DirectMethods/DualNumbers/dualnumber_tolerances.jl

@@ -15,20 +15,20 @@ This generates solutions u(t) = (t-θ)^5/5 that can be solved exactly with a 5th
 """
 function dyn!(du, u, p, t)
     θ = p[1]
-    du .= (t .- θ).^4.0
+    du .= (t .- θ) .^ 4.0
 end
 
 p = [1.0]
 
 prob = ODEProblem(dyn!, u0, tspan, p)
-sol  = solve(prob, Tsit5(), reltol=reltol, abstol=abstol)
+sol = solve(prob, Tsit5(), reltol = reltol, abstol = abstol)
 
 # We can see that the time steps increase with non-stop
 # @show diff(sol.t)
 
 function loss(p, sensealg)
     prob = ODEProblem(dyn!, u0, tspan, p)
-    sol = solve(prob, Tsit5(), sensealg=sensealg, reltol=reltol, abstol=abstol)
+    sol = solve(prob, Tsit5(), sensealg = sensealg, reltol = reltol, abstol = abstol)
     @show "Number of time steps: ", length(sol.t)
     sol.u[end][1]
 end
@@ -52,21 +52,23 @@ condition(u, t, integrator) = true
 function printstepsize!(integrator)
     if length(integrator.sol.t) > 1
         println("Stepsize at step ", length(integrator.sol.t), ":   ",
-            integrator.sol.t[end] - integrator.sol.t[end-1])
+            integrator.sol.t[end] - integrator.sol.t[end - 1])
     end
 end
 
 cb = DiscreteCallback(condition, printstepsize!)
 
 # g1 = Zygote.gradient(p -> loss(p, ForwardDiffSensitivity()), internalnorm = (u,t) -> sum(abs2,u/length(u)), p)
-g1 = Zygote.gradient(p -> solve(ODEProblem(dyn!, u0, tspan, p),
-                                Tsit5(),
-                                sensealg = ForwardDiffSensitivity(),
-                                saveat = 0.1,
-                                internalnorm = (u,t) -> sum(abs2, u/length(u)),
-                                callback = cb,
-                                reltol=1e-6,
-                                abstol=1e-6).u[end][1], p)
+g1 = Zygote.gradient(
+    p -> solve(ODEProblem(dyn!, u0, tspan, p),
+    Tsit5(),
+    sensealg = ForwardDiffSensitivity(),
+    saveat = 0.1,
+    internalnorm = (u, t) -> sum(abs2, u / length(u)),
+    callback = cb,
+    reltol = 1e-6,
+    abstol = 1e-6).u[end][1],
+    p)
 @show g1
 
 # Forward Sensitivity
@@ -82,15 +84,15 @@ g1 = Zygote.gradient(p -> solve(ODEProblem(dyn!, u0, tspan, p),
 
 # Corrected AD
 # g3 = ForwardDiff.gradient(p -> loss(p, nothing), p)
-g3 = Zygote.gradient(p -> solve(ODEProblem(dyn!, u0, tspan, p),
-                                Tsit5(),
-                                sensealg = ForwardDiffSensitivity(),
-                                # saveat = 0.1,
-                                # callback = cb,
-                                reltol=1e-6,
-                                abstol=1e-6).u[end][1], p)
+g3 = Zygote.gradient(
+    p -> solve(ODEProblem(dyn!, u0, tspan, p),
+    Tsit5(),
+    sensealg = ForwardDiffSensitivity(),     # saveat = 0.1,     # callback = cb,
+    reltol = 1e-6,
+    abstol = 1e-6).u[end][1],
+    p)
 @show g3
 
 @show grad_true(p)
 
-# Define customized RK(4) solver with given timesteps to show the divergence of forward sensitivities
+# Define customized RK(4) solver with given timesteps to show the divergence of forward sensitivities