Proposal: Benchmarks for phase space points

.gitignore

@@ -1,7 +1,15 @@
-/docs/build/
+docs/build/
 
 # File generated by Pkg, the package manager, based on a corresponding Project.toml
 # It records a fixed state of all packages used by the project. As such, it should not be
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+
+
+####
+# benchmarks
+####
+benchmark/*.json
+benchmark/plots
+

benchmark/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+BenchmarkPlots = "ab8c0f59-4072-4e0d-8f91-a91e1495eb26"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

benchmark/README.md

@@ -0,0 +1,27 @@
+# Benchmarks for QEDprocesses.jl
+
+This folder contains some micro benchmarks for the core functionality of
+`QEDprocesses.jl`.
+
+## Run benchmarks
+
+To run the benchmarks locally on your machine, go to the subfolder `benchmark` and execute
+the script `run_local.jl`:
+
+```console
+cd benchmark
+julia run_local.jl
+```
+
+This builds the benchmark suite, tunes the benchmarks, runs them and saves the results
+locally in a file `bench.json`.
+
+## Plot results
+
+For plotting the results, one can just run
+
+```console
+julia plot_bench.jl
+```
+
+which creates all benchmark plots and stores them in the folder `plots`.

benchmark/benchmark.jl

@@ -0,0 +1,108 @@
+using QEDbase
+using QEDprocesses
+using BenchmarkTools
+using StaticArrays
+using Random
+
+include("test_implementation.jl")
+
+const SUITE = BenchmarkGroup()
+
+SUITE["stateful"] = BenchmarkGroup()
+
+MOM_IN_1 = rand(SFourMomentum)
+MOM_IN_1K = rand(SFourMomentum, 1000)
+MOM_IN_1M = rand(SFourMomentum, 1000000)
+SUITE["stateful"]["construction"] = BenchmarkGroup()
+SUITE["stateful"]["construction"]["single"] = @benchmarkable ParticleStateful(
+    Incoming(), TestFermion(), mom, AllSpin()
+) setup = (mom = $MOM_IN_1)
+SUITE["stateful"]["construction"]["broadcast 1e3"] = @benchmarkable ParticleStateful.(
+    Incoming(), TestFermion(), mom, AllSpin()
+) setup = (mom = $MOM_IN_1K)
+SUITE["stateful"]["construction"]["broadcast 1e6"] = @benchmarkable ParticleStateful.(
+    Incoming(), TestFermion(), mom, AllSpin()
+) setup = (mom = $MOM_IN_1M)
+
+PARTICLE_STATEFUL = ParticleStateful(Incoming(), TestFermion(), MOM_IN_1, AllSpin())
+SUITE["stateful"]["accessor"] = BenchmarkGroup()
+SUITE["stateful"]["accessor"]["is_incoming"] = @benchmarkable QEDprocesses.is_incoming(
+    particle
+) setup = (particle = $PARTICLE_STATEFUL)
+SUITE["stateful"]["accessor"]["is_outgoing"] = @benchmarkable QEDprocesses.is_outgoing(
+    particle
+) setup = (particle = $PARTICLE_STATEFUL)
+SUITE["stateful"]["accessor"]["is_fermion"] = @benchmarkable is_fermion(particle) setup = (
+    particle = $PARTICLE_STATEFUL
+)
+SUITE["stateful"]["accessor"]["is_boson"] = @benchmarkable is_boson(particle) setup = (
+    particle = $PARTICLE_STATEFUL
+)
+SUITE["stateful"]["accessor"]["is_particle"] = @benchmarkable is_particle(particle) setup = (
+    particle = $PARTICLE_STATEFUL
+)
+SUITE["stateful"]["accessor"]["is_anti_particle"] = @benchmarkable is_anti_particle(
+    particle
+) setup = (particle = $PARTICLE_STATEFUL)
+SUITE["stateful"]["accessor"]["mass"] = @benchmarkable mass(particle) setup = (
+    particle = $PARTICLE_STATEFUL
+)
+SUITE["stateful"]["accessor"]["charge"] = @benchmarkable charge(particle) setup = (
+    particle = $PARTICLE_STATEFUL
+)
+
+MODEL = TestModel()
+PSDEF = TestPhasespaceDef()
+
+PROCESS_STRINGS = Dict(
+    "BF->BF" => TestProcess(
+        SVector(TestBoson(), TestFermion()), SVector(TestBoson(), TestFermion())
+    ),
+    "BB->BB" =>
+        TestProcess(SVector(TestBoson(), TestBoson()), SVector(TestBoson(), TestBoson())),
+    "FF->FF" => TestProcess(
+        SVector(TestFermion(), TestFermion()), SVector(TestFermion(), TestFermion())
+    ),
+    "FF->BF" => TestProcess(
+        SVector(TestFermion(), TestFermion()), SVector(TestBoson(), TestFermion())
+    ),
+    "FB->FF" => TestProcess(
+        SVector(TestBoson(), TestFermion()), SVector(TestFermion(), TestFermion())
+    ),
+)
+
+SUITE["phase space point"] = BenchmarkGroup()
+SUITE["phase space point"]["generate"] = BenchmarkGroup()
+for (proc_string, proc) in PROCESS_STRINGS
+    in_moms = rand(SFourMomentum, number_incoming_particles(proc))
+    out_moms = rand(SFourMomentum, number_outgoing_particles(proc))
+
+    in_particles = incoming_particles(proc)
+    out_particles = outgoing_particles(proc)
+
+    in_particle_states = SVector(
+        collect(
+            ParticleStateful(Incoming(), in_particles[i], in_moms[i]) for
+            i in 1:number_incoming_particles(proc)
+        )...,
+    )
+    out_particle_states = SVector(
+        collect(
+            ParticleStateful(Outgoing(), out_particles[i], out_moms[i]) for
+            i in 1:number_outgoing_particles(proc)
+        )...,
+    )
+
+    SUITE["phase space point"]["construction"][proc_string] = @benchmarkable PhaseSpacePoint(
+        $proc, $MODEL, $PSDEF, $in_particle_states, $out_particle_states
+    )
+    SUITE["phase space point"]["generate"][proc_string] = @benchmarkable generate_phase_space(
+        $proc, $MODEL, $PSDEF, $in_moms, $out_moms
+    )
+
+    psp = generate_phase_space(proc, MODEL, PSDEF, in_moms, out_moms)
+    SUITE["phase space point"]["momentum"][proc_string] = @benchmarkable momentum(
+        $psp, Incoming(), 1
+    )
+    SUITE["phase space point"]["indexing"][proc_string] = @benchmarkable $psp[Incoming(), 1]
+end

benchmark/plot_bench.jl

@@ -0,0 +1,44 @@
+using Pkg: Pkg
+Pkg.activate(".")
+
+using BenchmarkTools
+using BenchmarkPlots, StatsPlots
+
+include("utils.jl")
+
+plotpath = "plots"
+if !isdir(plotpath)
+    mkdir(plotpath)
+end
+
+result = BenchmarkTools.load("bench.json")[1]
+
+data = result["stateful"]["accessor"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: particle stateful accessor")
+savefig(joinpath(plotpath, "stateful_accessor.pdf"))
+
+data = result["stateful"]["construction"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: construction particle stateful")
+savefig(joinpath(plotpath, "stateful_construction.pdf"))
+
+data = result["phase space point"]["generate"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: generate phase space points")
+savefig(joinpath(plotpath, "phase_space_generation.pdf"))
+
+data = result["phase space point"]["construction"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: construction phase space points")
+savefig(joinpath(plotpath, "phase_space_construction.pdf"))
+
+data = result["phase space point"]["momentum"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: phase space points momentum access")
+savefig(joinpath(plotpath, "phase_space_momentum_access.pdf"))
+
+data = result["phase space point"]["indexing"]
+P = plot(data; yscale=:log10, ylim=_find_y_lims(data))
+plot!(P; title="benchmark: phase space points indexingaccess")
+savefig(joinpath(plotpath, "phase_space_momentum_indexing.pdf"))

benchmark/run_local.jl

@@ -0,0 +1,11 @@
+using Pkg: Pkg
+
+Pkg.activate(".")
+Pkg.develop(; path="./..")
+
+include("benchmark.jl")
+
+tune!(SUITE)
+result = run(SUITE; verbose=true)
+
+BenchmarkTools.save("bench.json", result)

benchmark/test_implementation.jl

@@ -0,0 +1,34 @@
+
+# dummy particles
+struct TestFermion <: FermionLike end
+QEDbase.mass(::TestFermion) = 1.0
+QEDbase.charge(::TestFermion) = 2.0
+struct TestBoson <: BosonLike end
+QEDbase.mass(::TestBoson) = 0.0
+QEDbase.charge(::TestBoson) = -2.0
+
+const PARTICLE_SET = [TestFermion(), TestBoson()]
+
+"""
+
+    TestProcess(rng,incoming_particles,outgoing_particles)
+
+"""
+struct TestProcess{IP<:AbstractVector,OP<:AbstractVector} <: AbstractProcessDefinition
+    incoming_particles::IP
+    outgoing_particles::OP
+end
+
+function TestProcess(rng::AbstractRNG, N_in::Int, N_out::Int)
+    in_particles = rand(rng, PARTICLE_SET, N_in)
+    out_particles = rand(rng, PARTICLE_SET, N_out)
+    return TestProcess(in_particles, out_particles)
+end
+
+QEDprocesses.incoming_particles(proc::TestProcess) = proc.incoming_particles
+QEDprocesses.outgoing_particles(proc::TestProcess) = proc.outgoing_particles
+
+struct TestPhasespaceDef <: AbstractPhasespaceDefinition end
+
+struct TestModel <: AbstractModelDefinition end
+QEDprocesses.fundamental_interaction_type(::TestModel) = :test_interaction

benchmark/utils.jl

@@ -0,0 +1,18 @@
+"""
+    _find_y_lims(data)
+
+Return the limits which encloses the data, i.e. the supremum and infimum of the data in powers of 10.
+This is useful for fining limit of logarithmic data.
+
+"""
+function _find_y_lims(data::BenchmarkTools.BenchmarkGroup)
+    t_min = time(minimum(data))
+    data_min = minimum(values(t_min))
+    ymin = 10^floor(log10(data_min))
+
+    t_max = time(maximum(data))
+    data_max = maximum(values(t_max))
+    ymax = 10^(ceil(log10(data_max)))
+
+    return (ymin, ymax)
+end

src/QEDprocesses.jl

@@ -45,4 +45,6 @@ include("momentum_generation.jl")
 include("propagators.jl")
 include("probabilities.jl")
 include("cross_sections.jl")
+
+include("patch_QEDbase.jl")
 end

src/patch_QEDbase.jl

@@ -0,0 +1,9 @@
+#############
+# Patches for `QEDbase.jl`
+# remove if this went into `QEDbase.jl`
+#############
+
+Broadcast.broadcastable(dir::Incoming) = Ref(dir)
+Broadcast.broadcastable(dir::Outgoing) = Ref(dir)
+Broadcast.broadcastable(part::AbstractParticleType) = Ref(part)
+Broadcast.broadcastable(spin_or_pol::AbstractSpinOrPolarization) = Ref(spin_or_pol)