rubydragon/FeynmanDiagrams.jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

WIP on DAG generation

Browse Source
This commit is contained in:
Rubydragon 2024-07-16 18:51:09 +02:00
parent 9859c3b928
commit a647801f12
6 changed files with 541 additions and 222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

158
notebooks/diagram_generation.ipynb
View File

@ -9,52 +9,23 @@
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Method definition number_particles(QEDbase.AbstractProcessDefinition, DIR, PT) where {DIR<:QEDbase.ParticleDirection, PT<:QEDbase.AbstractParticleType} in module QEDbase at /home/antonr/.julia/packages/QEDbase/gWHxL/src/interfaces/process.jl:208 overwritten in module FeynmanDiagramGenerator at /home/antonr/repos/FeynmanDiagramGenerator/src/QEDprocesses_patch.jl:23.\n",
"WARNING: Method definition (::Type{QEDcore.ParticleStateful{DIR, SPECIES, ELEMENT} where ELEMENT<:QEDbase.AbstractFourMomentum})(QEDbase.AbstractFourMomentum) where {DIR<:QEDbase.ParticleDirection, SPECIES<:QEDbase.AbstractParticleType} in module QEDcore at /home/antonr/.julia/packages/QEDcore/uVldP/src/phase_spaces/create.jl:7 overwritten in module MetagraphOptimization at /home/antonr/.julia/packages/MetagraphOptimization/mvCVq/src/QEDprocesses_patch.jl:15.\n",
"ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.\n"
]
}
],
"source": [
"using Revise\n",
"using FeynmanDiagramGenerator"
]
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": null,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"254-element Vector{VirtualParticle{GenericQEDProcess{Tuple{Photon, Photon, Photon, Photon, Photon, Photon, Photon, Electron}, Tuple{Photon, Electron}}, PT, 8, 2} where PT<:AbstractParticleType}:\n",
" positron: \t00000000 | 11\n",
" electron: \t00000001 | 10\n",
" positron: \t00000010 | 01\n",
" electron: \t00000011 | 00\n",
" positron: \t00000100 | 01\n",
" electron: \t00000101 | 00\n",
" positron: \t00001000 | 01\n",
" electron: \t00001001 | 00\n",
" positron: \t00010000 | 01\n",
" electron: \t00010001 | 00\n",
" ⋮\n",
" electron: \t11100001 | 10\n",
" positron: \t11100010 | 01\n",
" electron: \t11100011 | 00\n",
" positron: \t11100100 | 01\n",
" electron: \t11100101 | 00\n",
" positron: \t11101000 | 01\n",
" electron: \t11101001 | 00\n",
" positron: \t11110000 | 01\n",
" electron: \t11110001 | 00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"outputs": [],
"source": [
"proc = GenericQEDProcess(7, 1, 1, 1, 0, 0)\n",
"proc = GenericQEDProcess(3, 1, 1, 1, 0, 0)\n",
"all_particles = Set()\n",
"for fd in feynman_diagrams(proc)\n",
" push!(all_particles, virtual_particles(proc, fd)...)\n",
@ -65,125 +36,27 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": null,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{VirtualParticle, Vector{Tuple{VirtualParticle, VirtualParticle}}} with 254 entries:\n",
" positron: \t01000010 | 11 => [(photon: \t00000000 | 10, positron: \t01000010 | 01),…\n",
" positron: \t11010000 | 11 => [(photon: \t00000000 | 10, positron: \t11010000 | 01),…\n",
" electron: \t10100011 | 00 => [(photon: \t00000010 | 00, electron: \t10100001 | 00),…\n",
" positron: \t11000000 | 11 => [(photon: \t00000000 | 10, positron: \t11000000 | 01),…\n",
" positron: \t00010100 | 01 => [(photon: \t00000100 | 00, positron: \t00010000 | 01),…\n",
" positron: \t01001000 | 11 => [(photon: \t00000000 | 10, positron: \t01001000 | 01),…\n",
" positron: \t11000110 | 01 => [(photon: \t00000010 | 00, positron: \t11000100 | 01),…\n",
" electron: \t10001011 | 10 => [(photon: \t00000000 | 10, electron: \t10001011 | 00),…\n",
" positron: \t00100000 | 11 => [(photon: \t00000000 | 10, positron: \t00100000 | 01),…\n",
" positron: \t01010000 | 11 => [(photon: \t00000000 | 10, positron: \t01010000 | 01),…\n",
" electron: \t11101001 | 00 => [(photon: \t00001000 | 00, electron: \t11100001 | 00),…\n",
" electron: \t00100101 | 10 => [(photon: \t00000000 | 10, electron: \t00100101 | 00),…\n",
" electron: \t11000101 | 00 => [(photon: \t00000100 | 00, electron: \t11000001 | 00),…\n",
" electron: \t01101001 | 00 => [(photon: \t00001000 | 00, electron: \t01100001 | 00),…\n",
" positron: \t00101000 | 01 => [(photon: \t00001000 | 00, positron: \t00100000 | 01),…\n",
" positron: \t00011000 | 11 => [(photon: \t00000000 | 10, positron: \t00011000 | 01),…\n",
" positron: \t00011010 | 01 => [(photon: \t00000010 | 00, positron: \t00011000 | 01),…\n",
" electron: \t10000101 | 10 => [(photon: \t00000000 | 10, electron: \t10000101 | 00),…\n",
" electron: \t10000111 | 00 => [(photon: \t00000010 | 00, electron: \t10000101 | 00),…\n",
" ⋮ => ⋮"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"outputs": [],
"source": [
"pairs = FeynmanDiagramGenerator.particle_pairs(all_particles)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": null,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"280-element Vector{Tuple{VirtualParticle, VirtualParticle, VirtualParticle}}:\n",
" (photon: \t00000000 | 10, electron: \t01001011 | 00, positron: \t10110100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00010111 | 00, positron: \t11101000 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00110011 | 00, positron: \t11001100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t10111001 | 00, positron: \t01000110 | 01)\n",
" (photon: \t00000000 | 10, electron: \t10110101 | 00, positron: \t01001010 | 01)\n",
" (photon: \t00000000 | 10, electron: \t11110001 | 00, positron: \t00001110 | 01)\n",
" (photon: \t00000000 | 10, electron: \t11010011 | 00, positron: \t00101100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t01001101 | 00, positron: \t10110010 | 01)\n",
" (photon: \t00000000 | 10, electron: \t11101001 | 00, positron: \t00010110 | 01)\n",
" (photon: \t00000000 | 10, electron: \t10011011 | 00, positron: \t01100100 | 01)\n",
" ⋮\n",
" (photon: \t01000000 | 00, electron: \t00111001 | 00, positron: \t10000110 | 11)\n",
" (photon: \t01000000 | 00, electron: \t00000111 | 10, positron: \t10111000 | 01)\n",
" (photon: \t01000000 | 00, electron: \t00011011 | 00, positron: \t10100100 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10101011 | 00, positron: \t00010100 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10000111 | 10, positron: \t00111000 | 01)\n",
" (photon: \t01000000 | 00, electron: \t10100101 | 10, positron: \t00011010 | 01)\n",
" (photon: \t01000000 | 00, electron: \t00100111 | 00, positron: \t10011000 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10001101 | 10, positron: \t00110010 | 01)\n",
" (photon: \t01000000 | 00, electron: \t10011101 | 00, positron: \t00100010 | 11)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"outputs": [],
"source": [
"triples = FeynmanDiagramGenerator.total_particle_triples(all_particles)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"s: 40320, should be: 40320\n",
"number of triples: 280\n"
]
},
{
"data": {
"text/plain": [
"280-element Vector{Tuple{VirtualParticle, VirtualParticle, VirtualParticle}}:\n",
" (photon: \t00000000 | 10, electron: \t00001111 | 00, positron: \t11110000 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00010111 | 00, positron: \t11101000 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00011011 | 00, positron: \t11100100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00011101 | 00, positron: \t11100010 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00100111 | 00, positron: \t11011000 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00101011 | 00, positron: \t11010100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00101101 | 00, positron: \t11010010 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00110011 | 00, positron: \t11001100 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00110101 | 00, positron: \t11001010 | 01)\n",
" (photon: \t00000000 | 10, electron: \t00111001 | 00, positron: \t11000110 | 01)\n",
" ⋮\n",
" (photon: \t01000000 | 00, electron: \t10100101 | 10, positron: \t00011010 | 01)\n",
" (photon: \t01000000 | 00, electron: \t10100111 | 00, positron: \t00011000 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10101001 | 10, positron: \t00010110 | 01)\n",
" (photon: \t01000000 | 00, electron: \t10101011 | 00, positron: \t00010100 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10101101 | 00, positron: \t00010010 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10110001 | 10, positron: \t00001110 | 01)\n",
" (photon: \t01000000 | 00, electron: \t10110011 | 00, positron: \t00001100 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10110101 | 00, positron: \t00001010 | 11)\n",
" (photon: \t01000000 | 00, electron: \t10111001 | 00, positron: \t00000110 | 11)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"outputs": [],
"source": [
"function n(vp::VirtualParticle)\n",
" if !haskey(pairs, vp)\n",
@ -206,6 +79,15 @@
"\n",
"sort(triples)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"generate_DAG(proc)"
]
}
],
"metadata": {

2
src/FeynmanDiagramGenerator.jl
View File

@ -17,7 +17,7 @@ export FlatMatrix
export GenericQEDProcess, isphysical
export AbstractTreeLevelFeynmanDiagram, FeynmanVertex, FeynmanDiagram
export external_particles, virtual_particles, process
export external_particles, virtual_particles, process, generate_DAG
export VirtualParticle, number_contributions
export contributions, in_contributions, out_contributions

220
src/diagrams/diagrams.jl
View File

@ -39,12 +39,12 @@ function leaves(partition::TopologyPartition)
return partition.leaves1 + partition.leaves2 + partition.leaves3
end
#
# Feynman Diagram, tree-level, QED
#
struct FeynmanDiagram{N,E,U,T,M,FM} <: AbstractTreeLevelFeynmanDiagram where {N,E,U,T,M,FM<:FlatMatrix}
struct FeynmanDiagram{N,E,U,T,M,FM} <:
AbstractTreeLevelFeynmanDiagram where {N,E,U,T,M,FM<:FlatMatrix}
diagram_structure::FM
electron_permutation::NTuple{E,Int}
@ -59,19 +59,26 @@ struct FeynmanDiagram{N,E,U,T,M,FM} <: AbstractTreeLevelFeynmanDiagram where {N,
::Val{E},
::Val{U},
::Val{T},
::Val{M}
::Val{M},
) where {E,U,T,M}
@assert E == length(elec_perm)
@assert U == length(muon_perm)
@assert T == length(tauon_perm)
N = E + U + T
return new{N,E,U,T,M,FlatMatrix{Int64,2 * N - 2 + M,N}}(FlatMatrix(structure), NTuple{E,Int}(elec_perm), NTuple{U,Int}(muon_perm), NTuple{T,Int}(tauon_perm))
return new{N,E,U,T,M,FlatMatrix{Int64,2 * N - 2 + M,N}}(
FlatMatrix(structure),
NTuple{E,Int}(elec_perm),
NTuple{U,Int}(muon_perm),
NTuple{T,Int}(tauon_perm),
)
end
end
# representation of a virtual particle and the return type of the virtual_particles function
struct VirtualParticle{PROC<:QEDbase.AbstractProcessDefinition,PT<:QEDbase.AbstractParticleType,I,O}
struct VirtualParticle{
PROC<:QEDbase.AbstractProcessDefinition,PT<:QEDbase.AbstractParticleType,I,O
}
proc::PROC
species::PT
in_particle_contributions::NTuple{I,Bool}
@ -80,7 +87,10 @@ end
function Base.show(io::IO, vp::VirtualParticle)
pr = x -> x ? "1" : "0"
print(io, "$(vp.species): \t$(*(pr.(vp.in_particle_contributions)...)) | $(*(pr.(vp.out_particle_contributions)...))")
return print(
io,
"$(vp.species): \t$(*(pr.(vp.in_particle_contributions)...)) | $(*(pr.(vp.out_particle_contributions)...))",
)
end
function QEDbase.process(vp::VirtualParticle)
@ -96,11 +106,14 @@ out_contributions(vp::VirtualParticle) = vp.out_particle_contributions
contributions(vp::VirtualParticle) = ((in_contributions(vp), out_contributions(vp)))
is_virtual(vp::VirtualParticle) = number_contributions(vp) > 1
is_external(vp::VirtualParticle) = number_contributions(vp) == 1
import Base: +
# "addition" of the bool tuples
function +(a::Tuple{NTuple{I,Bool},NTuple{O,Bool}}, b::Tuple{NTuple{I,Bool},NTuple{O,Bool}}) where {I,O}
function +(
a::Tuple{NTuple{I,Bool},NTuple{O,Bool}}, b::Tuple{NTuple{I,Bool},NTuple{O,Bool}}
) where {I,O}
# realistically, there should never be "colliding" 1s. if there are there is probably an error and this should be asserted
#= for (i, j) in zip(a[1], b[1]) @assert !(i && j) end
for (i, j) in zip(a[2], b[2]) @assert !(i && j) end =#
@ -111,7 +124,9 @@ end
invert(::Electron) = Positron()
invert(::Positron) = Electron()
invert(::Photon) = Photon()
invert(::AbstractParticleType) = throw(InvalidInputError("unimplemented for this particle type"))
function invert(::AbstractParticleType)
throw(InvalidInputError("unimplemented for this particle type"))
end
function invert(virtual_particle::VirtualParticle)
I = length(virtual_particle.in_particle_contributions)
@ -120,7 +135,8 @@ function invert(virtual_particle::VirtualParticle)
virtual_particle.proc,
invert(virtual_particle.species),
ntuple(x -> !virtual_particle.in_particle_contributions[x], I),
ntuple(x -> !virtual_particle.out_particle_contributions[x], O))
ntuple(x -> !virtual_particle.out_particle_contributions[x], O),
)
end
# normalize the representation
@ -138,7 +154,12 @@ function normalize(virtual_particle::VirtualParticle)
end
end
function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDirection, species::AbstractParticleType, index::Int)
function _momentum_contribution(
proc::AbstractProcessDefinition,
dir::ParticleDirection,
species::AbstractParticleType,
index::Int,
)
I = number_incoming_particles(proc)
O = number_outgoing_particles(proc)
@ -151,7 +172,10 @@ function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDi
particles_seen += 1
end
if particles_seen == index
return (((is_incoming(dir) && x == c for x in 1:I)...,), ((is_outgoing(dir) && x == c for x in 1:O)...,))
return (
((is_incoming(dir) && x == c for x in 1:I)...,),
((is_outgoing(dir) && x == c for x in 1:O)...,),
)
end
end
@ -159,17 +183,25 @@ function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDi
end
function _fermion_type(proc::AbstractProcessDefinition, n::Int)
E = number_particles(proc, Incoming(), Electron()) + number_particles(proc, Outgoing(), Positron())
E =
number_particles(proc, Incoming(), Electron()) +
number_particles(proc, Outgoing(), Positron())
U = 0 # TODO add muons
T = 0 # TODO add tauons
M = number_particles(proc, Incoming(), Photon()) + number_particles(proc, Outgoing(), Photon())
M =
number_particles(proc, Incoming(), Photon()) +
number_particles(proc, Outgoing(), Photon())
N = E + U + T
# from the fermion index, get (Direction, Species, n) tuple, where n means it's the nth particle of that dir and species
if (n > 0 && n <= E)
electron_n = n
if electron_n > number_particles(proc, Incoming(), Electron())
return (Outgoing(), Positron(), electron_n - number_particles(proc, Incoming(), Electron()))
return (
Outgoing(),
Positron(),
electron_n - number_particles(proc, Incoming(), Electron()),
)
else
return (Incoming(), Electron(), electron_n)
end
@ -185,7 +217,11 @@ function _fermion_type(proc::AbstractProcessDefinition, n::Int)
# photon
photon_n = n - N
if photon_n > number_particles(proc, Incoming(), Photon())
return (Outgoing(), Photon(), photon_n - number_particles(proc, Incoming(), Photon()))
return (
Outgoing(),
Photon(),
photon_n - number_particles(proc, Incoming(), Photon()),
)
else
return (Incoming(), Photon(), photon_n)
end
@ -194,7 +230,11 @@ function _fermion_type(proc::AbstractProcessDefinition, n::Int)
electron_n = n - N - M
if electron_n > number_particles(proc, Outgoing(), Electron())
# incoming positron
return (Incoming(), Positron(), electron_n - number_particles(proc, Outgoing(), Electron()))
return (
Incoming(),
Positron(),
electron_n - number_particles(proc, Outgoing(), Electron()),
)
else
# outgoing electron
return (Outgoing(), Electron(), electron_n)
@ -245,7 +285,8 @@ end
Return true if the momenta contributions of `a` and `b` are disjunct.
"""
function disjunct(a::VirtualParticle, b::VirtualParticle)
for (a_contrib, b_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
for (a_contrib, b_contrib) in
Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
if b_contrib && a_contrib
return false
end
@ -260,7 +301,8 @@ end
Returns true if the set of particles contributing to `a` are contains the set of particles contributing to `b`.
"""
function contains(a::VirtualParticle, b::VirtualParticle)
for (a_contrib, b_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
for (a_contrib, b_contrib) in
Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
if b_contrib && !a_contrib
return false
end
@ -279,7 +321,8 @@ function make_up(a::VirtualParticle, b::VirtualParticle, c::VirtualParticle)
return false
end
# it should be unnecessary to check here that a and b can actually react. if a + b = c they must be able to if a, b and c all exist in the diagram.
for (a_contrib, b_contrib, c_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
for (a_contrib, b_contrib, c_contrib) in
Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
if c_contrib != a_contrib + b_contrib
return false
end
@ -294,7 +337,8 @@ end
Return true if a, b and c combined contain all external particles exactly once.
"""
function are_total(a::VirtualParticle, b::VirtualParticle, c::VirtualParticle)
for (a_contrib, b_contrib, c_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
for (a_contrib, b_contrib, c_contrib) in
Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
if a_contrib + b_contrib + c_contrib != 1
return false
end
@ -311,14 +355,17 @@ function particle_pairs(particles::Vector)
all_particles = sort([_pseudo_virtual_particles(proc)..., particles...])
# find pairs for every particle after the external ones (those can't have pairs)
for p_i in number_incoming_particles(proc)+number_outgoing_particles(proc)+1:length(all_particles)
for p_i in
(number_incoming_particles(proc) + number_outgoing_particles(proc) + 1):length(
all_particles
)
p = all_particles[p_i]
pairs[p] = Vector{Tuple{VirtualParticle,VirtualParticle}}()
# only need to consider external particles and virtual particles that come before p_i
for p_a_i in 1:p_i-2
for p_a_i in 1:(p_i - 2)
# and only partners between a and p_i
for p_b_i in p_a_i+1:p_i-1
for p_b_i in (p_a_i + 1):(p_i - 1)
p_a = all_particles[p_a_i]
p_b = all_particles[p_b_i]
@ -343,11 +390,20 @@ function total_particle_triples(particles::Vector)
photons = filter(p -> is_boson(particle_species(p)), working_set)
# make electrons a set for fast deletion
electrons = Set(filter(p -> is_fermion(particle_species(p)) && is_particle(particle_species(p)), working_set))
electrons = Set(
filter(
p -> is_fermion(particle_species(p)) && is_particle(particle_species(p)),
working_set,
),
)
# make positrons a set for fast lookup
positrons = Set(filter(p -> is_fermion(particle_species(p)) && is_anti_particle(particle_species(p)), working_set))
positrons = Set(
filter(
p -> is_fermion(particle_species(p)) && is_anti_particle(particle_species(p)),
working_set,
),
)
# no participant can have more than half the external particles, so every possible particle is contained here
# every photon has exactly one electron and positron partner
@ -358,7 +414,11 @@ function total_particle_triples(particles::Vector)
end
# create the only partner the ph and e could have together, then look for it in the actual positrons
expected_p = invert(VirtualParticle(proc, particle_species(e), (contributions(ph) + contributions(e))...))
expected_p = invert(
VirtualParticle(
proc, particle_species(e), (contributions(ph) + contributions(e))...
),
)
if expected_p in positrons
@assert are_total(ph, e, expected_p)
@ -376,10 +436,16 @@ end
Return a vector of `VirtualParticle` for each external particle. These are not actually virtual particles, but can be helpful as entry points.
"""
function _pseudo_virtual_particles(proc::AbstractProcessDefinition)
return sort(_external_particle.(proc, [1:number_incoming_particles(proc)+number_outgoing_particles(proc);]))
return sort(
_external_particle.(
proc, [1:(number_incoming_particles(proc) + number_outgoing_particles(proc));]
),
)
end
function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiagram{N,E,U,T,M,FM}) where {N,E,U,T,M,FM}
function virtual_particles(
proc::AbstractProcessDefinition, diagram::FeynmanDiagram{N,E,U,T,M,FM}
) where {N,E,U,T,M,FM}
fermion_lines = PriorityQueue{Int64,Int64}()
# count number of internal photons in each fermion line and make a priority queue for fermion line => number of internal photons
@ -452,11 +518,17 @@ function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiag
else # if the binding particle is unknown
# we have arrived at the "middle" of the line
# this line will be the unknown particle for the other lines
internal_photon_contributions[current_line] = cumulative_mom + unknown_photon_momentum
internal_photon_contributions[current_line] =
cumulative_mom + unknown_photon_momentum
# now we know that the fermion line that binding_particle binds to on the other end has one fewer unknown internal photons
fermion_lines[binding_particle] -= 1
# add the internal photon virtual particle
push!(result, VirtualParticle(proc, Photon(), (cumulative_mom + unknown_photon_momentum)...))
push!(
result,
VirtualParticle(
proc, Photon(), (cumulative_mom + unknown_photon_momentum)...
),
)
break
end
else # binding_particle is an external photon
@ -466,7 +538,7 @@ function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiag
end
end
return normalize.(result)[1:end-1]
return normalize.(result)[1:(end - 1)]
end
#
@ -484,14 +556,7 @@ end
function _feynman_structures(n::Int, m::Int)
f = labelled_plane_trees(n)
return ExternalPhotonIterator(
n,
m,
f,
1,
_external_photon_multiplicity(f[1], n, m),
1
)
return ExternalPhotonIterator(n, m, f, 1, _external_photon_multiplicity(f[1], n, m), 1)
end
function Base.length(it::ExternalPhotonIterator)
@ -510,7 +575,9 @@ function Base.iterate(iter::ExternalPhotonIterator, n::Nothing)
iter.fermion_structure_index += 1
iter.photon_structure_index = 1
iter.photon_structures = _external_photon_multiplicity(iter.fermion_structures[iter.fermion_structure_index], iter.N, iter.M)
iter.photon_structures = _external_photon_multiplicity(
iter.fermion_structures[iter.fermion_structure_index], iter.N, iter.M
)
else
iter.photon_structure_index += 1
end
@ -518,7 +585,9 @@ function Base.iterate(iter::ExternalPhotonIterator, n::Nothing)
return (iter.photon_structures[iter.photon_structure_index], nothing)
end
function _external_photon_multiplicity(fermion_structure::Vector{Vector{Int}}, n::Int, m::Int)
function _external_photon_multiplicity(
fermion_structure::Vector{Vector{Int}}, n::Int, m::Int
)
if m == 0
return [fermion_structure]
end
@ -529,7 +598,7 @@ function _external_photon_multiplicity(fermion_structure::Vector{Vector{Int}}, n
for line_index in eachindex(fermion_structure)
# go through indices start to end
for index in 1:length(fermion_structure[line_index])+1
for index in 1:(length(fermion_structure[line_index]) + 1)
# copy to prevent muting
new_photon_structure = deepcopy(fermion_structure)
# add new photon
@ -564,11 +633,17 @@ end
function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}) where {E,U,T,M}
N = E + U + T
f = FeynmanDiagram(iter.photon_structure, iter.e_perms[iter.e_index], iter.u_perms[iter.u_index], iter.t_perms[iter.t_index], iter.e, iter.u, iter.t, iter.m)
return (
f,
nothing
f = FeynmanDiagram(
iter.photon_structure,
iter.e_perms[iter.e_index],
iter.u_perms[iter.u_index],
iter.t_perms[iter.t_index],
iter.e,
iter.u,
iter.t,
iter.m,
)
return (f, nothing)
end
function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}, ::Nothing) where {E,U,T,M}
@ -594,11 +669,17 @@ function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}, ::Nothing) where {E
end
N = E + U + T
f = FeynmanDiagram(iter.photon_structure, iter.e_perms[iter.e_index], iter.u_perms[iter.u_index], iter.t_perms[iter.t_index], iter.e, iter.u, iter.t, iter.m)
return (
f,
nothing
f = FeynmanDiagram(
iter.photon_structure,
iter.e_perms[iter.e_index],
iter.u_perms[iter.u_index],
iter.t_perms[iter.t_index],
iter.e,
iter.u,
iter.t,
iter.m,
)
return (f, nothing)
end
function feynman_diagrams(proc::PROC) where {PROC<:GenericQEDProcess}
@ -625,10 +706,39 @@ function feynman_diagrams(in_particles::Tuple, out_particles::Tuple)
# left electrons -> left muons -> left tauons -> left photons -> right photons -> right electrons -> right muons -> right tauons
# a "left" fermion is simply an incoming fermion or outgoing antifermion of the type, while a "left" photon is an incoming photon, and the reverse for the right ones
f_iter = _feynman_structures(e + u + t, m)
e_perms = collect(permutations(Int[n+m+1:n+m+e;]))
u_perms = collect(permutations(Int[n+m+e+1:n+m+e+u;]))
t_perms = collect(permutations(Int[n+m+e+u+1:n+m+e+u+t;]))
e_perms = collect(permutations(Int[(n + m + 1):(n + m + e);]))
u_perms = collect(permutations(Int[(n + m + e + 1):(n + m + e + u);]))
t_perms = collect(permutations(Int[(n + m + e + u + 1):(n + m + e + u + t);]))
first_photon_structure, _ = iterate(f_iter)
return FeynmanDiagramIterator(Val(e), e_perms, 1, Val(u), u_perms, 1, Val(t), t_perms, 1, Val(m), f_iter, first_photon_structure)
return FeynmanDiagramIterator(
Val(e),
e_perms,
1,
Val(u),
u_perms,
1,
Val(t),
t_perms,
1,
Val(m),
f_iter,
first_photon_structure,
)
end
function virtual_particles(proc::GenericQEDProcess)
if !isempty(proc.virtual_particles_cache)
return proc.virtual_particles_cache
end
# use a set for deduplication
particles = Set{VirtualParticle}()
for fd in feynman_diagrams(proc)
push!(particles, virtual_particles(proc, fd)...)
end
# convert to vector
proc.virtual_particles_cache = [particles...]
return return proc.virtual_particles_cache
end

106
src/generic_process_def.jl
View File

@ -1,41 +1,115 @@
_assert_particle_type_tuple(::Tuple{}) = nothing
_assert_particle_type_tuple(t::Tuple{AbstractParticleType,Vararg}) = _assert_particle_type_tuple(t[2:end])
_assert_particle_type_tuple(t::Any) = throw(InvalidInputError("invalid input, provide a tuple of AbstractParticleTypes to construct a GenericQEDProcess"))
function _assert_particle_type_tuple(t::Tuple{AbstractParticleType,Vararg})
return _assert_particle_type_tuple(t[2:end])
end
function _assert_particle_type_tuple(t::Any)
throw(
InvalidInputError(
"invalid input, provide a tuple of AbstractParticleTypes to construct a GenericQEDProcess",
),
)
end
struct GenericQEDProcess{INT,OUTT} <: QEDbase.AbstractProcessDefinition where {INT<:Tuple,OUTT<:Tuple}
mutable struct GenericQEDProcess{INT,OUTT,INSP,OUTSP} <: AbstractProcessDefinition where {
INT<:Tuple,OUTT<:Tuple,INSP<:Tuple,OUTSP<:Tuple
}
incoming_particles::INT
outgoing_particles::OUTT
function GenericQEDProcess(in_particles::INT, out_particles::OUTT) where {INT<:Tuple,OUTT<:Tuple}
incoming_spins_pols::INSP
outgoing_spins_pols::OUTSP
virtual_particles_cache::Vector
function GenericQEDProcess(
in_particles::INT, out_particles::OUTT, in_sp::INSP, out_sp::OUTSP
) where {INT<:Tuple,OUTT<:Tuple,INSP<:Tuple,OUTSP<:Tuple}
_assert_particle_type_tuple(in_particles)
_assert_particle_type_tuple(out_particles)
return new{INT,OUTT}(in_particles, out_particles)
return new{INT,OUTT,INSP,OUTSP}(in_particles, out_particles, in_sp, out_sp, [])
end
"""
GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int)
Convenience constructor from numbers of input/output photons, electrons and positrons.
Uses `AllSpin()` and `AllPol()` for every particle's spin/pol by default.
"""
function GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int)
in_p = ntuple(i -> i <= in_ph ? Photon() : i <= in_ph + in_el ? Electron() : Positron(), in_ph + in_el + in_po)
out_p = ntuple(i -> i <= out_ph ? Photon() : i <= out_ph + out_el ? Electron() : Positron(), out_ph + out_el + out_po)
return GenericQEDProcess(in_p, out_p)
function GenericQEDProcess(
in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int
)
in_p = ntuple(i -> if i <= in_ph
Photon()
elseif i <= in_ph + in_el
Electron()
else
Positron()
end, in_ph + in_el + in_po)
out_p = ntuple(i -> if i <= out_ph
Photon()
elseif i <= out_ph + out_el
Electron()
else
Positron()
end, out_ph + out_el + out_po)
in_sp = tuple([i <= in_ph ? AllPol() : AllSpin() for i in 1:length(in_p)]...)
out_sp = tuple([i <= out_ph ? AllPol() : AllSpin() for i in 1:length(out_p)]...)
return GenericQEDProcess(in_p, out_p, in_sp, out_sp)
end
end
QEDprocesses.incoming_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT} = proc.incoming_particles
QEDprocesses.outgoing_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT} = proc.outgoing_particles
function spin_or_pol(
process::GenericQEDProcess,
dir::ParticleDirection,
species::AbstractParticleType,
n::Int,
)
i = 0
c = n
for p in particles(process, dir)
i += 1
if p == species
c -= 1
end
if c == 0
break
end
end
if c != 0 || n <= 0
throw(
InvalidInputError(
"could not get $n-th spin/pol of $dir $species, does not exist"
),
)
end
if is_incoming(dir)
return process.incoming_spins_pols[i]
elseif is_outgoing(dir)
return process.outgoing_spins_pols[i]
else
throw(InvalidInputError("unknown direction $(dir) given"))
end
end
function QEDprocesses.incoming_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT}
return proc.incoming_particles
end
function QEDprocesses.outgoing_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT}
return proc.outgoing_particles
end
function isphysical(proc::GenericQEDProcess)
return (number_particles(proc, Incoming(), Electron()) + number_particles(proc, Outgoing(), Positron()) ==
number_particles(proc, Incoming(), Positron()) + number_particles(proc, Outgoing(), Electron())) &&
number_particles(proc, Incoming()) + number_particles(proc, Outgoing()) >= 2
return (
number_particles(proc, Incoming(), Electron()) +
number_particles(proc, Outgoing(), Positron()) ==
number_particles(proc, Incoming(), Positron()) +
number_particles(proc, Outgoing(), Electron())
) && number_particles(proc, Incoming()) + number_particles(proc, Outgoing()) >= 2
end
function matrix_element(proc::GenericQEDProcess, psp::PhaseSpacePoint)
return nothing
end

10
src/metagraph_impl/compute.jl
View File

@ -7,7 +7,15 @@ struct ComputeTask_Triple <: AbstractComputeTask end # from a triple
struct ComputeTask_CollectTriples <: AbstractComputeTask end # sum over triples results and
# import compute so we don't have to repeat it all the time
import MetagraphOptimization.compute
import MetagraphOptimization: compute, compute_effort
compute_effort(::ComputeTask_BaseState) = 0
compute_effort(::ComputeTask_Propagator) = 0
compute_effort(::ComputeTask_Pair) = 0
compute_effort(::ComputeTask_CollectPairs) = 0
compute_effort(::ComputeTask_PropagatePairs) = 0
compute_effort(::ComputeTask_Triple) = 0
compute_effort(::ComputeTask_CollectTriples) = 0
struct BaseStateInput{PS_T<:AbstractParticleStateful,SPIN_POL_T<:AbstractSpinOrPolarization}
particle::PS_T

265
src/metagraph_impl/generation.jl
View File

@ -1,19 +1,264 @@
using MetagraphOptimization: insert_node!, insert_edge!, make_node, make_edge
function generate_DAG(proc::GenericQEDProcess)
# use a set for deduplication
particles = Set()
for fd in feynman_diagrams(proc)
push!(all_particles, virtual_particles(proc, fd)...)
_construction_string(::Incoming) = "Incoming()"
_construction_string(::Outgoing) = "Outgoing()"
_construction_string(::Electron) = "Electron()"
_construction_string(::Positron) = "Positron()"
_construction_string(::Photon) = "Photon()"
_construction_string(::PolX) = "PolX()"
_construction_string(::PolY) = "PolY()"
_construction_string(::SpinUp) = "SpinUp()"
_construction_string(::SpinDown) = "SpinDown()"
function _parse_particle(name::String)
local dir
if startswith(name, "inc_")
dir = Incoming()
elseif startswith(name, "out_")
dir = Outgoing()
else
throw(InvalidInputError("failed to parse particle direction from \"$name\""))
end
# convert to vector
particles = [particles...]
name = name[4:end]
external_particles = _pseudo_virtual_particles(proc)
pairs = particle_pairs(particles)
triples = total_particle_triples(particles)
local species
if startswith(name, "el")
species = Electron()
elseif startswith(name, "ph")
species = Photon()
elseif startswith(name, "po")
species = Positron()
else
throw(InvalidInputError("failed to parse particle species from name \"$name\""))
end
name = name[3:end]
local spin_pol
if startswith(name, "su")
spin_pol = SpinUp()
elseif startswith(name, "sd")
spin_pol = SpinDown()
elseif startswith(name, "px")
spin_pol = PolX()
elseif startswith(name, "py")
spin_pol = PolY()
else
throw(
InvalidInputError(
"failed to parse particle spin or polarization from \"$name\""
),
)
end
name = name[3:end]
index = parse(Int, name)
return (dir, species, spin_pol, index)
end
function input_expr(instance::GenericQEDProcess, name::String, psp_symbol::Symbol)
if startswith(name, "bs_")
(dir, species, spin_pol, index) = _parse_particle(name[4:end])
dir_str = _construction_string(dir)
species_str = _construction_string(species)
sp_str = _construction_string(spin_pol)
return Meta.parse(
"BaseStateInput(
ParticleStateful($dir_str, $species_str, momentum($psp_symbol, $dir_str, $species_str, $index)),
$sp_str,
)",
)
elseif startswith(name, "pr_")
index = parse(Int, name[4:end]) # get index of the virtual particle in the process
vp = virtual_particles(proc)[index]
return Meta.parse("PropagatorInput(
VirtualParticle(
process($psp_symbol),
$species_str,
$(vp.in_particle_contributions),
$(vp.out_particle_contributions)
),
$psp_symbol
)")
else
throw(InvalidInputError("failed to parse node name \"$name\""))
end
end
_species_str(::Photon) = "ph"
_species_str(::Electron) = "el"
_species_str(::Positron) = "po"
_spin_pol_str(::SpinUp) = "su"
_spin_pol_str(::SpinDown) = "sd"
_spin_pol_str(::PolX) = "px"
_spin_pol_str(::PolY) = "py"
_dir_str(::Incoming) = "inc"
_dir_str(::Outgoing) = "out"
_spin_pols(::AllSpin) = (SpinUp(), SpinDown())
_spin_pols(::SpinUp) = (SpinUp(),)
_spin_pols(::SpinDown) = (SpinDown(),)
_spin_pols(::AllPol) = (PolX(), PolY())
_spin_pols(::PolX) = (PolX(),)
_spin_pols(::PolY) = (PolY(),)
_is_external(p::VirtualParticle) = number_contributions(p) == 1
function _base_state_name(p::VirtualParticle)
proc = process(p)
dir = sum(in_contributions(p)) == 1 ? Incoming() : Outgoing()
# find particle in the contributions
index = 0
contribs = is_incoming(dir) ? in_contributions(p) : out_contributions(p)
for p in particles(proc, dir)
index += 1
if contribs[index]
break
end
end
species = particles(proc, dir)[index]
# find particle index of *this species*
species_index = 0
for i in 1:index
if particles(proc, dir)[i] == species
species_index += 1
end
end
spin_pol = spin_or_pol(proc, dir, species, species_index)
return string.(
"bs_$(_dir_str(dir))_$(_species_str(species))_",
_spin_pol_str.(_spin_pols(spin_pol)),
"_$(species_index)",
)
end
function generate_DAG(proc::GenericQEDProcess)
external_particles = _pseudo_virtual_particles(proc) # external particles that will be input to base_state tasks
particles = virtual_particles(proc) # virtual particles that will be input to propagator tasks
pairs = particle_pairs(particles) # pairs to generate the pair tasks
triples = total_particle_triples(particles) # triples to generate the triple tasks
graph = DAG()
# -- Base State Tasks --
base_state_task_outputs = Dict()
for dir in (Incoming(), Outgoing())
for species in (Electron(), Positron(), Photon())
for index in 1:number_particles(proc, dir, species)
for spin_pol in _spin_pols(spin_or_pol(proc, dir, species, index))
# gen entry nodes
# names are "bs_<dir>_<species>_<spin/pol>_<index>"
data_node_name = "bs_$(_dir_str(dir))_$(_species_str(species))_$(_spin_pol_str(spin_pol))_$(index)"
data_in = insert_node!(
graph,
make_node(DataTask(0), data_node_name);
track=false,
invalidate_cache=false,
)
# generate initial base_state tasks
compute_base_state = insert_node!(
graph,
make_node(ComputeTask_BaseState());
track=false,
invalidate_cache=false,
)
data_out = insert_node!(
graph, make_node(DataTask(0)); track=false, invalidate_cache=false
)
insert_edge!(graph, data_in, compute_base_state)
insert_edge!(graph, compute_base_state, data_out)
base_state_task_outputs[data_node_name] = data_out
end
end
end
end
# -- Propagator Tasks --
propagator_task_outputs = Dict()
vp_index = 0
for vp in virtual_particles(proc)
vp_index += 1
data_node_name = "pr_$vp_index"
data_in = insert_node!(
graph, make_node(DataTask(0)); track=false, invalidate_cache=false
)
compute_vp_propagator = insert_node!(
graph, make_node(ComputeTask_Propagator()); track=false, invalidate_cache=false
)
data_out = insert_node!(
graph, make_node(DataTask(0)); track=false, invalidate_cache=false
)
insert_edge!(graph, data_in, compute_vp_propagator)
insert_edge!(graph, compute_vp_propagator, data_out)
propagator_task_outputs[data_node_name] = data_out
end
# -- Pair Tasks --
pair_task_outputs = Dict()
for (product_particle, input_particle_vector) in pairs
# for all spins/pols of particles in product_particles do ...
pair_task_outputs[product_particle] = Vector()
for input_particles in input_particle_vector
particles_data_out_nodes = (Vector(), Vector())
c = 0
for p in input_particles
c += 1
if (is_external(p))
# grab from base_states (broadcast over _base_state_name because it is a tuple for different spin_pols)
push!.(
Ref(particles_data_out_nodes[c]),
getindex.(Ref(base_state_task_outputs), _base_state_name(p)),
)
else
# grab from propagated particles
push!(particles_date_out_nodes[c], pair_task_outputs[p])
end
end
for in_nodes in Iterators.product(input_particles...)
# make the compute pair nodes for every combination of the found input_particle_nodes to get all spin/pol combinations
#insert_node!(graph, )
end
# make the collect pair and propagate nodes
end
data_out_propagated = insert_node!(
graph, make_node(DataTask(0)); track=false, invalidate_caches=false
)
pair_task_outputs[p] = data_out_propagated
# ... end do
end
return graph
end