End of Day; Add topology generation and iterators

Project.toml

@@ -2,3 +2,6 @@ name = "FeynmanDiagramGenerator"
 uuid = "8f7ec768-9801-4863-9e03-72391fd2204f"
 authors = ["antonr <s1509337@msx.tu-dresden.de>"]
 version = "0.1.0"
+
+[deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/FeynmanDiagramGenerator.jl

@@ -2,8 +2,18 @@ module FeynmanDiagramGenerator
 
 export Forest
 
+export Electron, Positron, Muon, AntiMuon, Tauon, AntiTauon, Photon, ZBoson, WNBoson, WPBoson, Higgs
+export can_interact
+export QED
+
 include("trees/trees.jl")
 include("trees/iterator.jl")
 include("trees/print.jl")
 
+include("qft/particles.jl")
+include("qft/qft.jl")
+
+include("diagrams/diagrams.jl")
+include("diagrams/iterator.jl")
+
 end # module FeynmanDiagramGenerator

src/diagrams/diagrams.jl

@@ -0,0 +1,37 @@
+using Combinatorics
+
+struct FeynmanDiagramDefinition
+    n::Int
+end
+
+struct FeynmanDiagramTopology
+    # largest subtree
+    subtree1::AbstractTree
+    # second largest subtree
+    subtree2::AbstractTree
+    # third largest subtree
+    subtree3::AbstractTree
+end
+
+struct TopologyPartition
+    leaves1::Int
+    leaves2::Int
+    leaves3::Int
+end
+
+function is_valid(partition::TopologyPartition)
+    if partition.leaves2 > partition.leaves1 || partition.leaves3 > partition.leaves2
+        return false
+    end
+    if partition.leaves1 > partition.leaves2 + partition.leaves3
+        return false
+    end
+    if partition.leaves1 <= 0 || partition.leaves2 <= 0 || partition.leaves3 <= 0
+        return false
+    end
+    return true
+end
+
+function leaves(partition::TopologyPartition)
+    return partition.leaves1 + partition.leaves2 + partition.leaves3
+end

src/diagrams/iterator.jl

@@ -0,0 +1,92 @@
+
+"""
+    iterates inside one partition
+"""
+mutable struct FeynmanPartitionIterator
+    topology::TopologyPartition
+    forest1_iterator::Union{ForestIterator,Missing}
+    forest2_iterator::Union{ForestIterator,Missing}
+    forest3_iterator::Union{ForestIterator,Missing}
+end
+
+mutable struct FeynmanDiagramTopologyIterator
+    leaves::Int
+    partitions_state::Vector{Int64}
+    partition::FeynmanPartitionIterator
+end
+
+function FeynmanPartitionIterator(topology::TopologyPartition)
+    return FeynmanPartitionIterator(topology, ForestIterator(topology.leaves1), ForestIterator(topology.leaves2), ForestIterator(topology.leaves3))
+end
+
+function FeynmanDiagramTopologyIterator(def::FeynmanDiagramDefinition)
+    @assert def.n >= 4 "A Feynman diagram must have at least 4 legs"
+    (p, _) = iterate(partitions(def.n))
+    while length(p) != 3 || !is_valid(TopologyPartition(p[1], p[2], p[3]))
+        (p, _) = iterate(partitions(def.n), p)
+    end
+    return FeynmanDiagramTopologyIterator(def.n, p, FeynmanPartitionIterator(TopologyPartition(p[1], p[2], p[3])))
+end
+
+function FeynmanDiagramTopology(iterator::FeynmanPartitionIterator)
+    return FeynmanDiagramTopology(tree(iterator.forest1_iterator), tree(iterator.forest2_iterator), tree(iterator.forest3_iterator))
+end
+
+function FeynmanDiagramTopology(iterator::FeynmanDiagramTopologyIterator)
+    return FeynmanDiagramTopology(iterator.partition)
+end
+
+function _iterate(iterator::FeynmanPartitionIterator)
+    # todo: need symmetry breaking here
+    if !is_end(iterator.forest3_iterator)
+        _iterate(iterator.forest3_iterator)
+    elseif !is_end(iterator.forest2_iterator)
+        iterator.forest3_iterator = ForestIterator(iterator.topology.leaves3)
+        _iterate(iterator.forest2_iterator)
+    else
+        iterator.forest3_iterator = ForestIterator(iterator.topology.leaves3)
+        iterator.forest2_iterator = ForestIterator(iterator.topology.leaves2)
+        _iterate(iterator.forest1_iterator)
+    end
+    return nothing
+end
+
+function _iterate(iterator::FeynmanDiagramTopologyIterator)
+    if !is_end(iterator.partition)
+        _iterate(iterator.partition)
+    else
+        while true
+            ret = iterate(partitions(iterator.leaves), iterator.partitions_state)
+
+            if (isnothing(ret))
+                return false
+            end
+
+            (p, _) = ret
+            iterator.partitions_state = p
+            if length(p) == 3 && is_valid(TopologyPartition(p[1], p[2], p[3]))
+                iterator.partition = FeynmanPartitionIterator(TopologyPartition(p[1], p[2], p[3]))
+                break
+            end
+        end
+    end
+    return true
+end
+
+function is_end(iterator::FeynmanPartitionIterator)
+    return is_end(iterator.forest1_iterator) && is_end(iterator.forest2_iterator) && is_end(iterator.forest3_iterator)
+end
+
+function iterate(def::FeynmanDiagramDefinition)
+    state = FeynmanDiagramTopologyIterator(def)
+
+    return (FeynmanDiagramTopology(state), state)
+end
+
+function iterate(def::FeynmanDiagramDefinition, state::FeynmanDiagramTopologyIterator)
+    if !_iterate(state)
+        return nothing
+    end
+
+    return (FeynmanDiagramTopology(state), state)
+end

src/qft/particles.jl

@@ -0,0 +1,21 @@
+
+abstract type AbstractParticle end
+
+abstract type Fermion <: AbstractParticle end
+abstract type AntiFermion <: AbstractParticle end
+abstract type Boson <: AbstractParticle end
+
+struct Electron <: Fermion end
+struct Positron <: AntiFermion end
+struct Muon <: Fermion end
+struct AntiMuon <: AntiFermion end
+struct Tauon <: Fermion end
+struct AntiTauon <: AntiFermion end
+
+struct Photon <: Boson end
+struct ZBoson <: Boson end
+struct WNBoson <: Boson end
+struct WPBoson <: Boson end
+struct Higgs <: Boson end
+
+

src/qft/qft.jl

@@ -0,0 +1,38 @@
+
+"""
+    AbstractQFT
+
+Base type for any QFT. 
+"""
+abstract type AbstractQFT end
+
+struct Vertex
+    p1::AbstractParticle
+    p2::AbstractParticle
+    p3::AbstractParticle
+end
+
+Base.isequal(v1::Vertex, v2::Vertex) = Set([v1.p1, v1.p2, v1.p3]) == Set([v2.p1, v2.p2, v2.p3])
+
+is_valid(v::Vertex) = v.p1 != v.p2 && v.p2 != v.p3
+
+
+struct QFT <: AbstractQFT
+    allowed_vertices::Vector{Vertex}
+end
+
+function can_interact(theory::AbstractQFT, v::Vertex)
+    for allowed_vertex in theory.allowed_vertices
+        if isequal(allowed_vertex, v)
+            return true
+        end
+    end
+    return false
+end
+can_interact(theory::AbstractQFT, p1::AbstractParticle, p2::AbstractParticle, p3::AbstractParticle) = can_interact(theory, Vertex(p1, p2, p3))
+
+QED = QFT([
+    Vertex(Electron(), Positron(), Photon()),
+    Vertex(Muon(), AntiMuon(), Photon()),
+    Vertex(Tauon(), AntiTauon(), Photon())
+])

src/trees/iterator.jl

@@ -39,6 +39,8 @@ Return whether the given [`ForestIterator`](@ref) is exhausted.
 is_end(::Missing) = true
 is_end(iterator::ForestIterator) = iterator.left_leaves + 1 >= iterator.right_leaves && is_end(iterator.left_subiterator) && is_end(iterator.right_subiterator)
 
+leaves(iterator::ForestIterator) = iterator.left_leaves + iterator.right_leaves
+
 """
     tree(iterator::ForestIterator)
 

test/qft.jl

@@ -0,0 +1,23 @@
+using FeynmanDiagramGenerator
+
+import FeynmanDiagramGenerator.QFT
+
+@testset "QED" begin
+    theory = QED
+
+    for (p1, p2, p3, b) in [
+        (Electron(), Positron(), Photon(), true),
+        (Muon(), AntiMuon(), Photon(), true),
+        (Tauon(), AntiTauon(), Photon(), true),
+        (Electron(), Positron(), Positron(), false),
+        (Electron(), AntiMuon(), Photon(), false),
+        (Electron(), Positron(), ZBoson(), false)
+    ]
+        @test can_interact(theory, p1, p2, p3) == b
+        @test can_interact(theory, p1, p2, p3) == b
+        @test can_interact(theory, p2, p1, p3) == b
+        @test can_interact(theory, p2, p3, p1) == b
+        @test can_interact(theory, p3, p1, p2) == b
+        @test can_interact(theory, p3, p2, p1) == b
+    end
+end

test/runtests.jl

@@ -3,3 +3,7 @@ using SafeTestsets
 @safetestset "Trees" begin
     include("trees.jl")
 end
+
+@safetestset "QFT" begin
+    include("qft.jl")
+end