Add model, particles etc., add interaction_result and tests for it, add compute task types

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")

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
 
 """

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

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

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

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]

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

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

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

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