From 95354fb1682a529169ee5a2e951901e08de6520e Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Thu, 13 Jun 2024 15:49:09 +0200
Subject: [PATCH] Add labelled plane tree generation
---
Project.toml | 2 +
src/FeynmanDiagramGenerator.jl | 4 +
src/diagrams/QED.jl | 7 +-
src/qft/qft.jl | 4 +-
src/trees/labelled_plane_trees.jl | 134 ++++++++++++++++++++++++++++++
5 files changed, 146 insertions(+), 5 deletions(-)
create mode 100644 src/trees/labelled_plane_trees.jl
diff --git a/Project.toml b/Project.toml
index 42aae1d..b49002d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,3 +5,5 @@ version = "0.1.0"
[deps]
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
diff --git a/src/FeynmanDiagramGenerator.jl b/src/FeynmanDiagramGenerator.jl
index 38de335..95a127d 100644
--- a/src/FeynmanDiagramGenerator.jl
+++ b/src/FeynmanDiagramGenerator.jl
@@ -10,9 +10,12 @@ export external_particles, virtual_particles, process, vertices
export Forest
+export plane_trees, labelled_plane_trees
+
export can_interact
export QED
+include("trees/labelled_plane_trees.jl")
include("trees/trees.jl")
include("trees/iterator.jl")
include("trees/print.jl")
@@ -23,5 +26,6 @@ include("diagrams/interface.jl")
include("diagrams/vertex.jl")
include("diagrams/diagrams.jl")
include("diagrams/iterator.jl")
+include("diagrams/QED.jl")
end # module FeynmanDiagramGenerator
diff --git a/src/diagrams/QED.jl b/src/diagrams/QED.jl
index bc3f110..5a46acc 100644
--- a/src/diagrams/QED.jl
+++ b/src/diagrams/QED.jl
@@ -16,7 +16,9 @@ function plane_trees(n::Int64)
end
function iterate(lpt::LabelledPlaneTreeIterator)
- return ()
+ state = LabelledPlaneTreeIteratorState()
+
+
end
function iterate(lpt::LabelledPlaneTreeIterator, state)
@@ -26,5 +28,6 @@ function iterate(lpt::LabelledPlaneTreeIterator, state)
i += 1
end
end
-end
+ return (1, )
+end
diff --git a/src/qft/qft.jl b/src/qft/qft.jl
index fe6e341..cd6875c 100644
--- a/src/qft/qft.jl
+++ b/src/qft/qft.jl
@@ -32,7 +32,5 @@ 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())
+ Vertex(Electron(), Positron(), Photon())
])
diff --git a/src/trees/labelled_plane_trees.jl b/src/trees/labelled_plane_trees.jl
new file mode 100644
index 0000000..b10970d
--- /dev/null
+++ b/src/trees/labelled_plane_trees.jl
@@ -0,0 +1,134 @@
+using Combinatorics
+
+function _gen_trees(vec::Vector{Tuple{Int64,Vector{Int64}}})
+ # input is one tree structure
+ # node neighbours are partially filled from the front
+
+ n = length(vec)
+
+ # recursion ends if only two nodes are left with 1
+ sum_open_edges = sum(getindex.(vec, 1)) # total number of open edges
+ sum_open_edges -= sum(length.(getindex.(vec, 2))) # minus total number of assigned edges
+
+ @assert sum_open_edges != 0 "Cannot have empty tree"
+ @assert sum_open_edges % 2 == 0 "Cannot have uneven number of open edges"
+
+ # find first entry that has only one open edge
+ single_open_edge_nodes = Vector{Int64}()
+ i = 0
+ for node in vec
+ i += 1
+ if node[1] - length(node[2]) == 1
+ push!(single_open_edge_nodes, i)
+ end
+ end
+
+ @assert length(single_open_edge_nodes) >= 2 "There were less than two nodes with only one open edge, which cannot happen for valid trees"
+
+ results = Vector{Vector{Vector{Int64}}}()
+
+ if sum_open_edges == 2
+ # assign open edges
+ n1 = single_open_edge_nodes[1]
+ n2 = single_open_edge_nodes[2]
+ push!(vec[n1][2], n2)
+ push!(vec[n2][2], n1)
+
+ push!(results, getindex.(vec, 2))
+
+ return results
+ end
+
+ # choose first node
+ n1 = single_open_edge_nodes[1]
+
+ # iterate through second nodes: all nodes that have at least
+ for n2 in 1:n
+ if vec[n2][1] - length(vec[n2][2]) <= 1
+ # must have at least 2 open slots to be a partner for the edge
+ continue
+ end
+
+ # make new vec and connect edges
+ recurse_vec = deepcopy(vec)
+ push!(recurse_vec[n1][2], n2)
+ push!(recurse_vec[n2][2], n1)
+
+ generated_trees = _gen_trees(recurse_vec)
+ append!(results, generated_trees)
+ end
+
+ return results
+end
+
+function _labelled_plane_trees_unpermuted(N::Int64)
+ all_trees = Vector{Vector{Vector{Int64}}}()
+
+ parts = partitions(N - 2)
+ if N == 1
+ push!(all_trees, Vector{Vector{Vector{Int64}}}())
+ return all_trees
+ end
+
+ if N == 2
+ # fix empty partition when n = 2
+ parts = [Vector{Int64}()]
+ end
+
+ for origp in parts
+ p = copy(origp) # can't change original partition because of the iteration
+
+ n = length(p)
+ resize!(p, N)
+ for i in 1:N
+ if i > n
+ p[i] = 0
+ end
+ p[i] += 1
+ end
+
+ for mp in multiset_permutations(p, N)
+ # representation as vector of vectors
+ # each vector is a node
+ # each vector has the indices of the connected nodes
+ tree = [(l, Vector{Int64}()) for l in mp]
+
+ generated_trees = _gen_trees(tree)
+ append!(all_trees, generated_trees)
+ end
+ end
+
+ return all_trees
+end
+
+function _plane_tree_permutations_helper(t::Vector{Vector{Int64}}, depth::Int64)
+ if depth > length(t)
+ return [t]
+ end
+
+ all_perms = Vector{Vector{Vector{Int64}}}()
+ for permutation in permutations(t[depth])
+ permuted_t = deepcopy(t)
+ permuted_t[depth] = permutation
+ append!(all_perms, _plane_tree_permutations_helper(permuted_t, depth + 1))
+ end
+
+ return all_perms
+end
+
+# return all the possible edge permutations of a given labelled plane tree
+function _plane_tree_permutations(t::Vector{Vector{Int64}})
+ return _plane_tree_permutations_helper(t, 1)
+end
+
+function labelled_plane_trees(N::Int64)
+ all_trees = Vector{Vector{Vector{Int64}}}()
+
+ sizehint!(all_trees, factorial(3N - 3, 2N - 1))
+
+ for t in _labelled_plane_trees_unpermuted(N)
+ append!(all_trees, _plane_tree_permutations(t))
+ end
+
+ return all_trees
+end