From c2687cdc01d51a3e8f614d10580730472961dc7a Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 24 Nov 2023 01:26:10 +0100
Subject: [PATCH] Add model, particles etc., add interaction_result and tests
for it, add compute task types
---
src/MetagraphOptimization.jl | 20 ++++-
src/models/abc/compute.jl | 12 +--
src/models/abc/particle.jl | 16 ++--
src/models/qed/compute.jl | 80 +++++++++++++++++
src/models/qed/particle.jl | 168 +++++++++++++++++++++++++++++++++++
src/models/qed/types.jl | 51 +++++++++++
test/runtests.jl | 3 +
test/unit_tests_abcmodel.jl | 2 +-
test/unit_tests_execution.jl | 4 +-
test/unit_tests_qedmodel.jl | 51 +++++++++++
10 files changed, 389 insertions(+), 18 deletions(-)
create mode 100644 src/models/qed/compute.jl
create mode 100644 src/models/qed/particle.jl
create mode 100644 src/models/qed/types.jl
create mode 100644 test/unit_tests_qedmodel.jl
diff --git a/src/MetagraphOptimization.jl b/src/MetagraphOptimization.jl
index daea574..3df3ac1 100644
--- a/src/MetagraphOptimization.jl
+++ b/src/MetagraphOptimization.jl
@@ -54,7 +54,7 @@ export get_operations
# ABC model
export ParticleValue
export ParticleA, ParticleB, ParticleC
-export ABCProcessDescription, ABCProcessInput, ABCModel
+export ABCParticle, ABCProcessDescription, ABCProcessInput, ABCModel
export ComputeTaskABC_P
export ComputeTaskABC_S1
export ComputeTaskABC_S2
@@ -62,6 +62,16 @@ export ComputeTaskABC_V
export ComputeTaskABC_U
export ComputeTaskABC_Sum
+# QED model
+export PhotonStateful, FermionStateful, AntiFermionStateful
+export QEDParticle, QEDProcessDescription, QEDProcessInput, QEDModel
+export ComputeTaskQED_P
+export ComputeTaskQED_S1
+export ComputeTaskQED_S2
+export ComputeTaskQED_V
+export ComputeTaskQED_U
+export ComputeTaskQED_Sum
+
# code generation related
export execute
export parse_dag, parse_process
@@ -162,6 +172,14 @@ include("models/abc/properties.jl")
include("models/abc/parse.jl")
include("models/abc/print.jl")
+include("models/qed/types.jl")
+include("models/qed/particle.jl")
+include("models/qed/compute.jl")
+#include("models/qed/create.jl")
+#include("models/qed/properties.jl")
+#include("models/qed/parse.jl")
+#include("models/qed/print.jl")
+
include("devices/measure.jl")
include("devices/detect.jl")
include("devices/impl.jl")
diff --git a/src/models/abc/compute.jl b/src/models/abc/compute.jl
index e0d1110..13bcffa 100644
--- a/src/models/abc/compute.jl
+++ b/src/models/abc/compute.jl
@@ -14,12 +14,12 @@ end
"""
compute(::ComputeTaskABC_U, data::ABCParticleValue)
-Compute an outer edge. Return the particle value with the same particle and the value multiplied by an outer_edge factor.
+Compute an outer edge. Return the particle value with the same particle and the value multiplied by an ABC_outer_edge factor.
1 FLOP.
"""
function compute(::ComputeTaskABC_U, data::ABCParticleValue{P})::ABCParticleValue{P} where {P <: ABCParticle}
- return ABCParticleValue{P}(data.p, data.v * outer_edge(data.p))
+ return ABCParticleValue{P}(data.p, data.v * ABC_outer_edge(data.p))
end
"""
@@ -34,8 +34,8 @@ function compute(
data1::ABCParticleValue{P1},
data2::ABCParticleValue{P2},
)::ABCParticleValue where {P1 <: ABCParticle, P2 <: ABCParticle}
- p3 = preserve_momentum(data1.p, data2.p)
- dataOut = ABCParticleValue{typeof(p3)}(p3, data1.v * vertex() * data2.v)
+ p3 = ABC_preserve_momentum(data1.p, data2.p)
+ dataOut = ABCParticleValue{typeof(p3)}(p3, data1.v * ABC_vertex() * data2.v)
return dataOut
end
@@ -59,7 +59,7 @@ function compute(
@assert isapprox(data1.p.momentum.py, -data2.p.momentum.py, rtol = 0.001, atol = sqrt(eps())) "py: $(data1.p.momentum.py) vs. $(data2.p.momentum.py)"
@assert isapprox(data1.p.momentum.pz, -data2.p.momentum.pz, rtol = 0.001, atol = sqrt(eps())) "pz: $(data1.p.momentum.pz) vs. $(data2.p.momentum.pz)"
=#
- inner = inner_edge(data1.p)
+ inner = ABC_inner_edge(data1.p)
return data1.v * inner * data2.v
end
@@ -71,7 +71,7 @@ Compute inner edge (1 input particle, 1 output particle).
11 FLOP.
"""
function compute(::ComputeTaskABC_S1, data::ABCParticleValue{P})::ABCParticleValue{P} where {P <: ABCParticle}
- return ABCParticleValue{P}(data.p, data.v * inner_edge(data.p))
+ return ABCParticleValue{P}(data.p, data.v * ABC_inner_edge(data.p))
end
"""
diff --git a/src/models/abc/particle.jl b/src/models/abc/particle.jl
index fbaef32..f374ed0 100644
--- a/src/models/abc/particle.jl
+++ b/src/models/abc/particle.jl
@@ -119,48 +119,48 @@ function square(p::ABCParticle)
end
"""
- inner_edge(p::ABCParticle)
+ ABC_inner_edge(p::ABCParticle)
Return the factor of the inner edge with the given (virtual) particle.
Takes 10 effective FLOP. (3 here + 7 in square(p))
"""
-function inner_edge(p::ABCParticle)
+function ABC_inner_edge(p::ABCParticle)
return 1.0 / (square(p) - mass(typeof(p)) * mass(typeof(p)))
end
"""
- outer_edge(p::ABCParticle)
+ ABC_outer_edge(p::ABCParticle)
Return the factor of the outer edge with the given (real) particle.
Takes 0 effective FLOP.
"""
-function outer_edge(p::ABCParticle)
+function ABC_outer_edge(p::ABCParticle)
return 1.0
end
"""
- vertex()
+ ABC_vertex()
Return the factor of a vertex.
Takes 0 effective FLOP since it's constant.
"""
-function vertex()
+function ABC_vertex()
i = 1.0
lambda = 1.0 / 137.0
return i * lambda
end
"""
- preserve_momentum(p1::ABCParticle, p2::ABCParticle)
+ ABC_preserve_momentum(p1::ABCParticle, p2::ABCParticle)
Calculate and return a new particle from two given interacting ones at a vertex.
Takes 4 effective FLOP.
"""
-function preserve_momentum(p1::ABCParticle, p2::ABCParticle)
+function ABC_preserve_momentum(p1::ABCParticle, p2::ABCParticle)
t3 = interaction_result(typeof(p1), typeof(p2))
p3 = t3(p1.momentum + p2.momentum)
return p3
diff --git a/src/models/qed/compute.jl b/src/models/qed/compute.jl
new file mode 100644
index 0000000..844e3f1
--- /dev/null
+++ b/src/models/qed/compute.jl
@@ -0,0 +1,80 @@
+
+"""
+ compute(::ComputeTaskQED_P, data::QEDParticleValue)
+
+Return the particle as is and initialize the Value.
+"""
+function compute(::ComputeTaskQED_P, data::QEDParticleValue{P})::QEDParticleValue{P} where {P <: QEDParticle}
+ # TODO do we actually need this for anything?
+ return QEDParticleValue{P}(data.p, one(DiracMatrix))
+end
+
+"""
+ compute(::ComputeTaskQED_U, data::QEDParticleValue)
+
+Compute an outer edge. Return the particle value with the same particle and the value multiplied by an outer_edge factor.
+"""
+function compute(::ComputeTaskQED_U, data::QEDParticleValue{P})::QEDParticleValue{P} where {P <: QEDParticle}
+ return QEDParticleValue{P}(
+ data.p,
+ base_state(particle(data.p), direction(data.p), momentum(data.p), spin_or_pol(data.p)), # will return a SLorentzVector{ComplexF64}, BiSpinor or AdjointBiSpinor
+ )
+end
+
+"""
+ compute(::ComputeTaskQED_V, data1::QEDParticleValue, data2::QEDParticleValue)
+
+Compute a vertex. Preserve momentum and particle types (e + gamma->p etc.) to create resulting particle, multiply values together and times a vertex factor.
+"""
+function compute(
+ ::ComputeTaskQED_V,
+ data1::QEDParticleValue{P1},
+ data2::QEDParticleValue{P2},
+)::QEDParticleValue where {P1 <: QEDParticle, P2 <: QEDParticle}
+ p3 = QED_preserve_momentum(data1.p, data2.p)
+ dataOut = QEDParticleValue{typeof(p3)}(p3, data1.v * ABC_vertex() * data2.v)
+ return dataOut
+end
+
+"""
+ compute(::ComputeTaskQED_S2, data1::QEDParticleValue, data2::QEDParticleValue)
+
+Compute a final inner edge (2 input particles, no output particle).
+
+For valid inputs, both input particles should have the same momenta at this point.
+
+12 FLOP.
+"""
+function compute(::ComputeTaskQED_S2, data1::ParticleValue{P}, data2::ParticleValue{P})::ComplexF64
+ #=
+ @assert isapprox(abs(data1.p.momentum.E), abs(data2.p.momentum.E), rtol = 0.001, atol = sqrt(eps())) "E: $(data1.p.momentum.E) vs. $(data2.p.momentum.E)"
+ @assert isapprox(data1.p.momentum.px, -data2.p.momentum.px, rtol = 0.001, atol = sqrt(eps())) "px: $(data1.p.momentum.px) vs. $(data2.p.momentum.px)"
+ @assert isapprox(data1.p.momentum.py, -data2.p.momentum.py, rtol = 0.001, atol = sqrt(eps())) "py: $(data1.p.momentum.py) vs. $(data2.p.momentum.py)"
+ @assert isapprox(data1.p.momentum.pz, -data2.p.momentum.pz, rtol = 0.001, atol = sqrt(eps())) "pz: $(data1.p.momentum.pz) vs. $(data2.p.momentum.pz)"
+ =#
+ inner = QED_inner_edge(data1.p)
+ return data1.v * inner * data2.v
+end
+
+"""
+ compute(::ComputeTaskQED_S1, data::QEDParticleValue)
+
+Compute inner edge (1 input particle, 1 output particle).
+
+11 FLOP.
+"""
+function compute(::ComputeTaskQED_S1, data::QEDParticleValue{P})::QEDParticleValue where {P <: QEDParticle}
+ # TODO invert P for result (incoming becomes outgoing, outgoing becomes incoming)
+ return QEDParticleValue{P}(data.p, data.v * QED_inner_edge(data.p))
+end
+
+"""
+ compute(::ComputeTaskQED_Sum, data::Vector{Float64})
+
+Compute a sum over the vector. Use an algorithm that accounts for accumulated errors in long sums with potentially large differences in magnitude of the summands.
+
+Linearly many FLOP with growing data.
+"""
+function compute(::ComputeTaskQED_Sum, data::Vector{Float64})::Float64
+ return sum_kbn(data)
+end
diff --git a/src/models/qed/particle.jl b/src/models/qed/particle.jl
new file mode 100644
index 0000000..673e649
--- /dev/null
+++ b/src/models/qed/particle.jl
@@ -0,0 +1,168 @@
+import QEDbase.mass
+
+"""
+ QEDModel <: AbstractPhysicsModel
+
+Singleton definition for identification of the QED-Model.
+"""
+struct QEDModel <: AbstractPhysicsModel end
+
+"""
+ QEDParticle
+
+Base type for all particles in the [`QEDModel`](@ref).
+
+Its template parameter specifies the particle's direction.
+
+The concrete types contain singletons of the types that they are, like `Photon` and `Electron` from QEDbase, and their state descriptions.
+"""
+abstract type QEDParticle{Direction <: ParticleDirection} <: AbstractParticle end
+
+"""
+ QEDProcessDescription <: AbstractProcessDescription
+
+A description of a process in the QED-Model. Contains the input and output particles.
+
+See also: [`in_particles`](@ref), [`out_particles`](@ref), [`parse_process`](@ref)
+"""
+struct QEDProcessDescription <: AbstractProcessDescription
+ inParticles::Dict{Type{<:QEDParticle{Incoming}}, Int}
+ outParticles::Dict{Type{<:QEDParticle{Outgoing}}, Int}
+end
+
+"""
+ QEDProcessInput <: AbstractProcessInput
+
+Input for a QED Process. Contains the [`QEDProcessDescription`](@ref) of the process it is an input for, and the values of the in and out particles.
+
+See also: [`gen_process_input`](@ref)
+"""
+struct QEDProcessInput <: AbstractProcessInput
+ process::QEDProcessDescription
+ inParticles::Vector{QEDParticle}
+ outParticles::Vector{QEDParticle}
+end
+
+QEDParticleValue{ParticleType <: QEDParticle} =
+ ParticleValue{ParticleType, <:Union{BiSpinor, AdjointBiSpinor, DiracMatrix, SLorentzVector{ComplexF64}, ComplexF64}}
+
+"""
+ PhotonStateful <: QEDParticle
+
+A photon of the [`QEDModel`](@ref) with its state.
+"""
+struct PhotonStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
+ momentum::SFourMomentum
+ # this will maybe change to the full polarization vector? or do i need both
+ polarization::AbstractPolarization
+end
+
+"""
+ FermionStateful <: QEDParticle
+
+A fermion of the [`QEDModel`](@ref) with its state.
+"""
+struct FermionStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
+ momentum::SFourMomentum
+ spin::AbstractSpin
+end
+
+"""
+ AntiFermionStateful <: QEDParticle
+
+An anti-fermion of the [`QEDModel`](@ref) with its state.
+"""
+struct AntiFermionStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
+ momentum::SFourMomentum
+ spin::AbstractSpin
+end
+
+"""
+ interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
+
+For 2 given (non-equal) particle types, return the third.
+"""
+function interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
+ @assert false "Invalid interaction between particles of types $t1 and $t2"
+end
+
+interaction_result(::Type{FermionStateful{Incoming}}, ::Type{FermionStateful{Outgoing}}) = PhotonStateful{Incoming}
+interaction_result(::Type{FermionStateful{Incoming}}, ::Type{AntiFermionStateful{Incoming}}) = PhotonStateful{Incoming}
+interaction_result(::Type{FermionStateful{Incoming}}, ::Type{<:PhotonStateful}) = FermionStateful{Outgoing}
+
+interaction_result(::Type{FermionStateful{Outgoing}}, ::Type{FermionStateful{Incoming}}) = PhotonStateful{Incoming}
+interaction_result(::Type{FermionStateful{Outgoing}}, ::Type{AntiFermionStateful{Outgoing}}) = PhotonStateful{Incoming}
+interaction_result(::Type{FermionStateful{Outgoing}}, ::Type{<:PhotonStateful}) = FermionStateful{Incoming}
+
+# antifermion mirror
+interaction_result(::Type{AntiFermionStateful{Incoming}}, t2::Type{<:QEDParticle}) =
+ interaction_result(FermionStateful{Outgoing}, t2)
+interaction_result(::Type{AntiFermionStateful{Outgoing}}, t2::Type{<:QEDParticle}) =
+ interaction_result(FermionStateful{Incoming}, t2)
+
+# photon commutativity
+interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:QEDParticle}) = interaction_result(t2, t1)
+
+# but prevent stack overflow
+function interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:PhotonStateful})
+ @assert false "Invalid interaction between particles of types $t1 and $t2"
+end
+
+"""
+ types(::QEDModel)
+
+Return a Vector of the possible types of particle in the [`QEDModel`](@ref).
+"""
+function types(::QEDModel)
+ return [PhotonStateful, FermionStateful, AntiFermionStateful]
+end
+
+# type piracy?
+String(::Type{Incoming}) = "Incoming"
+String(::Type{Outgoing}) = "Outgoing"
+
+function String(::Type{PhotonStateful{Dir}}) where {Dir <: ParticleDirection}
+ return "Photon" * String(Dir)
+end
+function String(::Type{FermionStateful{Dir}}) where {Dir <: ParticleDirection}
+ return "Fermion" * String(Dir)
+end
+function String(::Type{AntiFermionStateful{Dir}}) where {Dir <: ParticleDirection}
+ return "Anti-Fermion" * String(Dir)
+end
+
+@inline particle(::PhotonStateful) = Photon()
+@inline particle(::FermionStateful) = Fermion()
+@inline particle(::AntiFermionStateful) = AntiFermion()
+
+@inline momentum(p::PhotonStateful)::SFourMomentum{ComplexF64} = p.momentum
+@inline momentum(p::FermionStateful)::SFourMomentum{ComplexF64} = p.momentum
+@inline momentum(p::AntiFermionStateful)::SFourMomentum{ComplexF64} = p.momentum
+
+@inline spin_or_pol(p::PhotonStateful)::AbstractPolarization = p.polarization
+@inline spin_or_pol(p::FermionStateful)::AbstractSpin = p.spin
+@inline spin_or_pol(p::AntiFermionStateful)::AbstractSpin = p.spin
+
+@inline direction(::PhotonStateful{Dir})::ParticleDirection = Dir()
+@inline direction(::FermionStateful{Dir})::ParticleDirection = Dir()
+@inline direction(::AntiFermionStateful{Dir})::ParticleDirection = Dir()
+
+"""
+ QED_vertex()
+
+Return the factor of a vertex in a QED feynman diagram.
+"""
+@inline function QED_vertex()::SLorentzVector{DiracMatrix}
+ return ComplexF64(0, 1) * gamma()
+end
+
+"""
+ QED_preserve_momentum(p1::QEDParticle, p2::QEDParticle)
+
+Calculate and return a new particle from two given interacting ones at a vertex.
+"""
+function QED_preserve_momentum(p1::QEDParticle, p2::QEDParticle)
+ t3 = interaction_result(typeof(p1), typeof(p2))
+ p3 = t3(p1.momentum + p2.momentum)
+ return p3
+end
diff --git a/src/models/qed/types.jl b/src/models/qed/types.jl
new file mode 100644
index 0000000..9923014
--- /dev/null
+++ b/src/models/qed/types.jl
@@ -0,0 +1,51 @@
+"""
+ ComputeTaskQED_S1 <: AbstractComputeTask
+
+S task with a single child.
+"""
+struct ComputeTaskQED_S1 <: AbstractComputeTask end
+
+"""
+ ComputeTaskQED_S2 <: AbstractComputeTask
+
+S task with two children.
+"""
+struct ComputeTaskQED_S2 <: AbstractComputeTask end
+
+"""
+ ComputeTaskQED_P <: AbstractComputeTask
+
+P task with no children.
+"""
+struct ComputeTaskQED_P <: AbstractComputeTask end
+
+"""
+ ComputeTaskQED_V <: AbstractComputeTask
+
+v task with two children.
+"""
+struct ComputeTaskQED_V <: AbstractComputeTask end
+
+"""
+ ComputeTaskQED_U <: AbstractComputeTask
+
+u task with a single child.
+"""
+struct ComputeTaskQED_U <: AbstractComputeTask end
+
+"""
+ ComputeTaskQED_Sum <: AbstractComputeTask
+
+Task that sums all its inputs, n children.
+"""
+mutable struct ComputeTaskQED_Sum <: AbstractComputeTask
+ children_number::Int
+end
+
+"""
+ QED_TASKS
+
+Constant vector of all tasks of the QED-Model.
+"""
+QED_TASKS =
+ [ComputeTaskQED_S1, ComputeTaskQED_S2, ComputeTaskQED_P, ComputeTaskQED_V, ComputeTaskQED_U, ComputeTaskQED_Sum]
diff --git a/test/runtests.jl b/test/runtests.jl
index 83c0c8e..3930b57 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,6 +18,9 @@ end
@safetestset "ABC-Model Unit Tests" begin
include("unit_tests_abcmodel.jl")
end
+@safetestset "QED-Model Unit Tests" begin
+ include("unit_tests_qedmodel.jl")
+end
@safetestset "Node Reduction Unit Tests" begin
include("node_reduction.jl")
end
diff --git a/test/unit_tests_abcmodel.jl b/test/unit_tests_abcmodel.jl
index 0814ba0..cff7ecc 100644
--- a/test/unit_tests_abcmodel.jl
+++ b/test/unit_tests_abcmodel.jl
@@ -19,5 +19,5 @@ testparticles = [ParticleA(def_momentum), ParticleB(def_momentum), ParticleC(def
end
@testset "Vertex" begin
- @test isapprox(MetagraphOptimization.vertex(), 1 / 137.0)
+ @test isapprox(MetagraphOptimization.ABC_vertex(), 1 / 137.0)
end
diff --git a/test/unit_tests_execution.jl b/test/unit_tests_execution.jl
index f2b2b81..c2046d1 100644
--- a/test/unit_tests_execution.jl
+++ b/test/unit_tests_execution.jl
@@ -38,8 +38,8 @@ function ground_truth_graph_result(input::ABCProcessInput)
@test isapprox(getMass2(diagram1_C.momentum), getMass2(diagram1_Cp.momentum))
@test isapprox(getMass2(diagram2_C.momentum), getMass2(diagram2_Cp.momentum))
- inner1 = MetagraphOptimization.inner_edge(diagram1_C)
- inner2 = MetagraphOptimization.inner_edge(diagram2_C)
+ inner1 = MetagraphOptimization.ABC_inner_edge(diagram1_C)
+ inner2 = MetagraphOptimization.ABC_inner_edge(diagram2_C)
diagram1_result = inner1 * constant
diagram2_result = inner2 * constant
diff --git a/test/unit_tests_qedmodel.jl b/test/unit_tests_qedmodel.jl
new file mode 100644
index 0000000..ead6388
--- /dev/null
+++ b/test/unit_tests_qedmodel.jl
@@ -0,0 +1,51 @@
+using MetagraphOptimization
+using QEDbase
+
+import MetagraphOptimization.interaction_result
+
+def_momentum = SFourMomentum(1.0, 0.0, 0.0, 0.0)
+
+testparticleTypes = [
+ PhotonStateful{Incoming},
+ PhotonStateful{Outgoing},
+ FermionStateful{Incoming},
+ FermionStateful{Outgoing},
+ AntiFermionStateful{Incoming},
+ AntiFermionStateful{Outgoing},
+]
+#testparticles = [ParticleA(def_momentum), ParticleB(def_momentum), ParticleC(def_momentum)]
+
+function caninteract(t1::Type{<:QEDParticle}, t2::Type{<:QEDParticle})
+ if (t1 == t2)
+ return false
+ end
+ if (t1 <: PhotonStateful && t2 <: PhotonStateful)
+ return false
+ end
+
+ for (p1, p2) in [(t1, t2), (t2, t1)]
+ if (p1 == FermionStateful{Incoming} && p2 == AntiFermionStateful{Outgoing})
+ return false
+ end
+ if (p1 == FermionStateful{Outgoing} && p2 == AntiFermionStateful{Incoming})
+ return false
+ end
+ end
+
+ return true
+end
+
+function issame(t1::Type{<:QEDParticle}, t2::Type{<:QEDParticle})
+ return !caninteract(t1, t2)
+end
+
+@testset "Interaction Result" begin
+ for p1 in testparticleTypes, p2 in testparticleTypes
+ if !caninteract(p1, p2)
+ @test_throws AssertionError interaction_result(p1, p2)
+ else
+ @test interaction_result(p1, p2) in setdiff(testparticleTypes, [p1, p2])
+ @test issame(interaction_result(p1, p2), interaction_result(p2, p1))
+ end
+ end
+end