Add benchmarking infrastructure (#151)

* Add PkgBenchmark and BenchmarkTools as deps

* Add files for benchmarks (Kepler problem)

* Add run_benchmarks.jl

Running this file, creates a comparison from the branch `test/v0.8.12` to this one.

* Add new case for Kepler benchmarks

* Add many-spin benchmarks

* Add date to output filenames, and Dates.jl to Project.toml

* Add benchmark/README.md and adapt .gitignore

* Small fixes for run_benchmark.jl

* Small fixes

* Move run_benchmarks.jl to benchmark/

... a cleaner move.

* Update benchmark/README.md

.gitignore

@@ -4,3 +4,4 @@ examples/JuliaCon2017/.ipynb_checkpoints/
 docs/build/
 docs/site/
 .vscode/
+benchmark/20*

Project.toml

@@ -4,24 +4,29 @@ repo = "https://github.com/PerezHz/TaylorIntegration.jl.git"
 version = "0.8.11"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
 Espresso = "6912e4f1-e036-58b0-9138-08d1e6358ea9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
+BenchmarkTools = "1.3"
 DiffEqBase = "6"
 Elliptic = "0.5, 1.0"
 Espresso = "0.6.1"
 OrdinaryDiffEq = "5, 6"
+PkgBenchmark = "0.2"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
 Requires = "0.5.2, 1"

benchmark/README.md

@@ -0,0 +1,14 @@
+# Benchmarks
+
+The simplest way to run the benchmarks is to use the command-line
+and, being at the root-directory of the package, run the command
+```bash
+$ julia --project=. benchmark/run_benchmarks.jl
+```
+
+The results are stored in three files in the `benchmark/` directory:
+`yyyy-mm-dd_comparacion.md` contains a markdown file comparing the
+target and base versions, and further information of the parameters
+used for the benchmarks. The actual benchmarks are stored in
+`yyyy-mm-dd_resultadosBase.json` and `yyyy-mm-dd_resultadosTarg.json`,
+for the base and target branches.

benchmark/benchmarks.jl

@@ -0,0 +1,38 @@
+using BenchmarkTools, TaylorIntegration
+
+const directorypath = dirname(@__FILE__)
+
+const BT = BenchmarkTools
+
+SUITE = BenchmarkGroup()
+
+# ==========
+#  Include files to benchmarks with parameters
+# ==========
+#
+BT.DEFAULT_PARAMETERS.samples = 25
+# BT.DEFAULT_PARAMETERS.seconds = 12
+
+include("kepler_benchmarks.jl")
+include("manyspin_benchmarks.jl")
+
+# ToDo:
+# pendulum
+# Lyapunov spectrum
+# jet transport
+# common interface
+# Poincare maps
+# @taylorize with @threads
+
+# ==========
+# Tune
+# ==========
+println("Load `tune.json`")
+paramspath = joinpath(directorypath, "tune.json")
+if isfile(paramspath)
+    loadparams!(SUITE, BT.load(paramspath)[1], :evals, :samples);
+else
+    tune!(SUITE)
+    BT.save(paramspath, params(SUITE));
+end
+

benchmark/kepler_benchmarks.jl

@@ -0,0 +1,137 @@
+# ==========
+# Kepler problem
+# ==========
+
+# ==========
+# Define equations of motion
+# ==========
+# Using `^`
+@taylorize function kepler1!(dq, q, p, t)
+    local μ = p
+    r_p3d2 = (q[1]^2+q[2]^2)^1.5
+    #
+    dq[1] = q[3]
+    dq[2] = q[4]
+    dq[3] = μ * q[1]/r_p3d2
+    dq[4] = μ * q[2]/r_p3d2
+    #
+    return nothing
+end
+
+@taylorize function kepler2!(dq, q, p, t)
+    nn, μ = p
+    r2 = zero(q[1])
+    for i = 1:nn
+        r2_aux = r2 + q[i]^2
+        r2 = r2_aux
+    end
+    r_p3d2 = r2^(3/2)
+    for j = 1:nn
+        dq[j] = q[nn+j]
+        dq[nn+j] = μ*q[j]/r_p3d2
+    end
+    nothing
+end
+
+# Using `sqrt`
+@taylorize function kepler5!(dq, q, p, t)
+    local μ = -1.0
+    r = sqrt(q[1]^2+q[2]^2)
+    r_p3d2 = r^3
+    dq[1] = q[3]
+    dq[2] = q[4]
+    dq[3] = μ * q[1] / r_p3d2
+    dq[4] = μ * q[2] / r_p3d2
+    return nothing
+end
+
+@taylorize function kepler6!(dq, q, p, t)
+    local NN = 2
+    local μ = p
+    r2 = zero(q[1])
+    for i = 1:NN
+        r2_aux = r2 + q[i]^2
+        r2 = r2_aux
+    end
+    r = sqrt(r2)
+    r_p3d2 = r^3
+    for j = 1:NN
+        dq[j] = q[NN+j]
+        dq[NN+j] =  μ * q[j] / r_p3d2
+    end
+    nothing
+end
+
+# ==========
+# Run benchmarks
+# ==========
+let
+    # const TI = TaylorIntegration
+    local _abstol = 1.0e-20
+    local _order = 28
+    local t0 = 0.0
+    local tf1 = 2π*1000.0
+    local maxsteps1 = 500_000
+    local tf2 = 2π*10_000.0
+    local maxsteps2 = 50_000_000
+    local tf3 = 2π*100_000.0
+    local maxsteps3 = 500_000_000
+
+    pars = (2, -1.0)
+    q0 = [0.2, 0.0, 0.0, 3.0]
+
+    # ==========
+    # Kepler_parsed
+    # ==========
+    SUITE["Kepler"] = BenchmarkGroup()
+
+    SUITE["Kepler"]["kepler1-1"] = @benchmarkable taylorinteg(
+        kepler1!, $q0, $t0, $tf1, $_order, $_abstol, -1.0, maxsteps=$maxsteps1)
+    SUITE["Kepler"]["kepler2-1"] = @benchmarkable taylorinteg(
+        kepler2!, $q0, $t0, $tf1, $_order, $_abstol, $pars, maxsteps=$maxsteps1)
+
+    SUITE["Kepler"]["kepler5-1"] = @benchmarkable taylorinteg(
+        kepler5!, $q0, $t0, $tf1, $_order, $_abstol, maxsteps=$maxsteps1)
+    SUITE["Kepler"]["kepler6-1"] = @benchmarkable taylorinteg(
+        kepler6!, $q0, $t0, $tf1, $_order, $_abstol, -1.0, maxsteps=$maxsteps1)
+
+    SUITE["Kepler"]["kepler1-2"] = @benchmarkable taylorinteg(
+        kepler1!, $q0, $t0, $tf2, $_order, $_abstol, -1.0, maxsteps=$maxsteps2)
+    SUITE["Kepler"]["kepler2-2"] = @benchmarkable taylorinteg(
+        kepler2!, $q0, $t0, $tf2, $_order, $_abstol, $pars, maxsteps=$maxsteps2)
+
+    SUITE["Kepler"]["kepler5-2"] = @benchmarkable taylorinteg(
+        kepler5!, $q0, $t0, $tf2, $_order, $_abstol, maxsteps=$maxsteps2)
+    SUITE["Kepler"]["kepler6-2"] = @benchmarkable taylorinteg(
+        kepler6!, $q0, $t0, $tf2, $_order, $_abstol, -1.0, maxsteps=$maxsteps2)
+
+    SUITE["Kepler"]["kepler1-3"] = @benchmarkable taylorinteg(
+        kepler1!, $q0, $t0, $tf3, $_order, $_abstol, -1.0, maxsteps=$maxsteps3)
+    SUITE["Kepler"]["kepler2-3"] = @benchmarkable taylorinteg(
+        kepler2!, $q0, $t0, $tf3, $_order, $_abstol, $pars, maxsteps=$maxsteps3)
+
+    SUITE["Kepler"]["kepler5-3"] = @benchmarkable taylorinteg(
+        kepler5!, $q0, $t0, $tf3, $_order, $_abstol, maxsteps=$maxsteps3)
+    SUITE["Kepler"]["kepler6-3"] = @benchmarkable taylorinteg(
+        kepler6!, $q0, $t0, $tf3, $_order, $_abstol, -1.0, maxsteps=$maxsteps3)
+
+    # # ==========
+    # # KeplerNotParsed
+    # # ==========
+    # SUITE["KeplerNotParsed"] = BenchmarkGroup()
+
+    # SUITE["KeplerNotParsed"]["kepler1"] = @benchmarkable taylorinteg(
+    #     kepler1!, $q0, $t0, $tf, $_order, $_abstol, -1.0, maxsteps=$maxsteps,
+    #     parse_eqs=false)
+    # SUITE["KeplerNotParsed"]["kepler2"] = @benchmarkable taylorinteg(
+    #     kepler2!, $q0, $t0, $tf, $_order, $_abstol, $pars, maxsteps=$maxsteps,
+    #     parse_eqs=false)
+    #
+    # SUITE["KeplerNotParsed"]["kepler5"] = @benchmarkable taylorinteg(
+    #     kepler5!, $q0, $t0, $tf, $_order, $_abstol, maxsteps=$maxsteps,
+    #     parse_eqs=false)
+    # SUITE["KeplerNotParsed"]["kepler6"] = @benchmarkable taylorinteg(
+    #     kepler6!, $q0, $t0, $tf, $_order, $_abstol, -1.0, maxsteps=$maxsteps,
+    #     parse_eqs=false)
+    #
+end

benchmark/manyspin_benchmarks.jl

@@ -0,0 +1,113 @@
+# ==========
+# Classical many-spin system
+# ==========
+
+using Random
+
+# ==========
+# Equations of motion
+# ==========
+@taylorize function spin_odes!(dq, q, params, t)
+    local nn = params[1]
+    local B = params[2]
+    local J = params[3]
+
+    # Initilize vars
+    x = Array{eltype(q),1}(undef, nn)
+    y = Array{eltype(q),1}(undef, nn)
+    z = Array{eltype(q),1}(undef, nn)
+    for i = 1:nn
+        x[i] = q[3i-2]
+        y[i] = q[3i-1]
+        z[i] = q[3i  ]
+    end
+
+    local zz = zero(x[1])
+    sum_s = Array{eltype(q),1}(undef, nn)
+    aux   = Array{eltype(q),1}(undef, nn)
+    for i = 1:nn
+        sum_s[i] = zz
+        for a = 1:nn
+            aux[i] = sum_s[i]
+            xJ = J[a,i] * x[a]
+            sum_s[i] = aux[i] + xJ
+        end
+        dq[3i-2] = -B[i] * y[i]
+        dq[3i-1] =  B[i] * x[i] + z[i] * sum_s[i]
+        dq[3i  ] =              - y[i] * sum_s[i]
+    end
+
+    return dq
+end
+
+# ==========
+# Some functions
+# ==========
+function fix_B(nn)
+    @assert nn < 21
+    local B20 = [0.80000305175781250, 0.85261507034301753, 1.1022420883178712, 0.98346004486083982,
+        1.0131068229675293, 0.88758363723754885, 0.81881780624389644, 1.0715457916259765,
+        1.0717185020446778, 1.1738771438598632, 0.95340080261230464, 1.0077665328979493,
+        1.1323861122131347, 0.81382875442504887, 0.82138462066650386, 1.0118800163269044,
+        1.0684597015380859, 0.80307922363281248, 0.95336618423461916, 0.82673683166503908]
+        return B20[1:nn]
+end
+function fix_J(nn; α=1.4)
+    JJ = zeros(nn, nn)
+    @inbounds for d = 1:nn-1
+        val = 1.0 / d^α
+        @inbounds for i = 1:nn-d
+            JJ[i, i+d] = val
+            JJ[i+d, i] = val
+        end
+    end
+    return JJ
+end
+function ini_conditions!(q, sz)
+    # sz = similar(q, length(q)//3)
+    @. sz = 2.0 * (rand() - 0.5)
+
+    @inbounds for i in eachindex(sz)
+        as = sqrt(1.0 - sz[i]*sz[i])
+        pha = 2.0 * pi * rand()
+        q[3i-2] = as * cos(pha)
+        q[3i-1] = as * sin(pha)
+        q[3i  ] = sz[i]
+    end
+
+    return nothing
+end
+
+# ==========
+# Run benchmarks
+# ==========
+let
+    local _abstol = 1.0e-20
+    local _order = 26
+    local L1 = 3
+    local pars1 = (L1, fix_B(L1), fix_J(L1))
+    local L2 = 7
+    local pars2 = (L2, fix_B(L2), fix_J(L2))
+    local t0 = 0.0
+    local tf1 = 100.0
+    local tf2 = 1000.0
+    local maxsteps = 100_000
+    Random.seed!(1023)
+    local q1 = zeros(3*L1)
+    ini_conditions!(q1, similar(q1, L1))
+    Random.seed!(1023)
+    local q2 = zeros(3*L2)
+    ini_conditions!(q2, similar(q2, L2))
+
+    SUITE["ManySpin"] = BenchmarkGroup()
+
+    SUITE["ManySpin"]["manyspin1-1"] = @benchmarkable taylorinteg(
+        spin_odes!, $q1, $t0, $tf1, $_order, $_abstol, $pars1, maxsteps=$maxsteps)
+    SUITE["ManySpin"]["manyspin1-2"] = @benchmarkable taylorinteg(
+        spin_odes!, $q1, $t0, $tf2, $_order, $_abstol, $pars1, maxsteps=$maxsteps)
+
+    SUITE["ManySpin"]["manyspin2-1"] = @benchmarkable taylorinteg(
+        spin_odes!, $q2, $t0, $tf1, $_order, $_abstol, $pars2, maxsteps=$maxsteps)
+    SUITE["ManySpin"]["manyspin2-2"] = @benchmarkable taylorinteg(
+        spin_odes!, $q2, $t0, $tf2, $_order, $_abstol, $pars2, maxsteps=$maxsteps)
+end

benchmark/pendulum_benchmarks.jl

@@ -0,0 +1,41 @@
+# ==========
+# Multiple pendula
+# ==========
+
+# ==========
+# Define equations of motion
+# ==========
+@taylorize function multpendula1!(dx, x, p, t)
+    for i in p[2]
+        dx[i] = x[p[1]+i]
+        dx[i+p[1]] = -sin( x[i] )
+    end
+    return nothing
+end
+
+# ==========
+# Constants
+# ==========
+const _abstol = 1.0e-20
+const _order = 20
+const pars = (3, 1:3)
+const t0 = 0.0
+const tf1 = 1000.0
+const tf2 = 10_000.0
+const maxsteps = 10_000
+const q0 = [pi-0.001, 0.0, pi-0.001, 0.0,  pi-0.001, 0.0]
+
+# ==========
+# Run benchmarks
+# ==========
+SUITE["Pendumum"] = BenchmarkGroup()
+
+SUITE["Pendumum"]["pendulum1-1"] = @benchmarkable taylorinteg(
+    multpendula1!, $q0, $t0, $tf1, $_order, $_abstol, $pars, maxsteps=$maxsteps)
+SUITE["Pendumum"]["pendulum1-2"] = @benchmarkable taylorinteg(
+    multpendula1!, $q0, $t0, $tf2, $_order, $_abstol, $pars, maxsteps=$maxsteps)
+
+SUITE["Pendumum"]["pendulum1-1"] = @benchmarkable taylorinteg(
+    multpendula1!, $q0, $t0, $tf1, $_order, $_abstol, $pars, maxsteps=$maxsteps)
+SUITE["Pendumum"]["pendulum1-2"] = @benchmarkable taylorinteg(
+    multpendula1!, $q0, $t0, $tf2, $_order, $_abstol, $pars, maxsteps=$maxsteps)