diff --git a/notebooks/diagram_gen.ipynb b/notebooks/diagram_gen.ipynb
index f2c3dad..ede40b3 100644
--- a/notebooks/diagram_gen.ipynb
+++ b/notebooks/diagram_gen.ipynb
@@ -6,7 +6,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "using Revise; using QEDbase; using QEDprocesses; using MetagraphOptimization; using BenchmarkTools"
+    "using Revise; using QEDbase; using QEDprocesses; using MetagraphOptimization; using BenchmarkTools; using DataStructures"
    ]
   },
   {
@@ -23,15 +23,24 @@
     "AntifermionTypes = Union{Type{FermionStateful{Outgoing}},Type{AntiFermionStateful{Incoming}}};\n"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Feynman Diagram Generation\n",
+    "Data structure: Ve\n",
+    "\n"
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 83,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "gen_diagram_helper (generic function with 4 methods)"
+       "gen_diagrams (generic function with 1 method)"
       ]
      },
      "metadata": {},
@@ -39,156 +48,975 @@
     }
    ],
    "source": [
-    "function gen_diagram_helper(photons::Int, fermionsIn::Int, fermionsOut::Int, antifermionsIn::Int, antifermionsOut::Int, trace::String = \"\", fermionsN::Int = 1, photonsN::Int = 1)\n",
-    "    subtrace = \"ph\"\n",
-    "    if (fermionsIn == 1)\n",
-    "        subtrace *= \" + f(in)\"\n",
+    "import MetagraphOptimization.QEDParticle\n",
+    "import MetagraphOptimization.caninteract\n",
+    "import MetagraphOptimization.interaction_result\n",
+    "import MetagraphOptimization.propagation_result\n",
+    "import MetagraphOptimization.direction\n",
+    "\n",
+    "import Base.copy\n",
+    "import Base.hash\n",
+    "import Base.==\n",
+    "import Base.show\n",
+    "\n",
+    "struct FeynmanParticle\n",
+    "    particle::Type{<:QEDParticle}\n",
+    "    id::Int\n",
+    "end\n",
+    "\n",
+    "show(io::IO, p::FeynmanParticle) = print(io, \"$(String(p.particle))_$(String(direction(p.particle)))_$(p.id)\")\n",
+    "\n",
+    "struct FeynmanVertex\n",
+    "    in1::FeynmanParticle\n",
+    "    in2::FeynmanParticle\n",
+    "    out::FeynmanParticle\n",
+    "end\n",
+    "\n",
+    "show(io::IO, v::FeynmanVertex) = print(io, \"$(v.in1) + $(v.in2) -> $(v.out)\")\n",
+    "\n",
+    "function Base.hash(v::FeynmanVertex)\n",
+    "    return hash(Set{FeynmanParticle}([v.in1, v.in2]))\n",
+    "end\n",
+    "\n",
+    "function ==(v1::FeynmanVertex, v2::FeynmanVertex)\n",
+    "    return Set{FeynmanParticle}([v1.in1, v1.in2]) == Set{FeynmanParticle}([v2.in1, v2.in2])\n",
+    "end\n",
+    "\n",
+    "struct FeynmanTie\n",
+    "    in1::FeynmanParticle\n",
+    "    in2::FeynmanParticle\n",
+    "end\n",
+    "\n",
+    "show(io::IO, t::FeynmanTie) = print(io, \"$(t.in1) -- $(t.in2)\")\n",
+    "\n",
+    "function Base.hash(t::FeynmanTie)\n",
+    "    return hash(Set{FeynmanParticle}([t.in1, t.in2]))\n",
+    "end\n",
+    "\n",
+    "function ==(t1::FeynmanTie, t2::FeynmanTie)\n",
+    "    return Set{FeynmanParticle}([t1.in1, t1.in2]) == Set{FeynmanParticle}([t2.in1, t2.in2])\n",
+    "end\n",
+    "\n",
+    "struct FeynmanDiagram\n",
+    "    vertices::Vector{Set{FeynmanVertex}}\n",
+    "    ties::Vector{Set{FeynmanTie}}\n",
+    "    particles::Vector{FeynmanParticle}\n",
+    "end\n",
+    "\n",
+    "function FeynmanDiagram(pd::QEDProcessDescription)\n",
+    "    parts = Vector{FeynmanParticle}()\n",
+    "    for (type, n) in pd.inParticles\n",
+    "        for i in 1:n\n",
+    "            push!(parts, FeynmanParticle(type, i))\n",
+    "        end\n",
     "    end\n",
-    "    if (antifermionsIn == 1)\n",
-    "        subtrace *= \" + af(in)\"\n",
+    "    for (type, n) in pd.outParticles\n",
+    "        for i in 1:n\n",
+    "            push!(parts, FeynmanParticle(type, i))\n",
+    "        end\n",
     "    end\n",
-    "    if (fermionsOut == 1)\n",
-    "        subtrace *= \" + f(out)\"\n",
+    "    return FeynmanDiagram([], [], parts)\n",
+    "end\n",
+    "\n",
+    "function ==(d1::FeynmanDiagram, d2::FeynmanDiagram)\n",
+    "    return d1.vertices == d2.vertices && d1.ties == d2.ties && d1.particles == d2.particles\n",
+    "end\n",
+    "\n",
+    "function Base.hash(d::FeynmanDiagram)\n",
+    "    return hash((d.vertices, d.ties, d.particles))\n",
+    "end\n",
+    "\n",
+    "function show(io::IO, d::FeynmanDiagram)\n",
+    "    print(io, \"Initial Particles: [\")\n",
+    "    first = true\n",
+    "    for p in d.particles\n",
+    "        if first\n",
+    "            first = false\n",
+    "            print(io, \"$p\")\n",
+    "        else\n",
+    "            print(io, \", $p\")\n",
+    "        end\n",
     "    end\n",
-    "    if (antifermionsOut == 1)\n",
-    "        subtrace *= \" + af(out)\"\n",
+    "    print(io, \"]\\n\")\n",
+    "    for l in eachindex(d.vertices)\n",
+    "        print(io, \"Virtuality Level $l:\\n\")\n",
+    "        print(io, \"  Vertices: [\")\n",
+    "        first = true\n",
+    "        for v in d.vertices[l]\n",
+    "            if first\n",
+    "                first = false\n",
+    "                print(io, \"$v\")\n",
+    "            else\n",
+    "                print(io, \", $v\")\n",
+    "            end\n",
+    "        end\n",
+    "        print(io, \"]\\n  Ties: [\")\n",
+    "        first = true\n",
+    "        for t in d.ties[l]\n",
+    "            if first\n",
+    "                first = false\n",
+    "                print(io, \"$t\")\n",
+    "            else\n",
+    "                print(io, \", $t\")\n",
+    "            end\n",
+    "        end\n",
+    "        print(io, \"]\\n\")\n",
     "    end\n",
-    "    if photons == 1 && (fermionsIn + antifermionsOut == 1) && (antifermionsIn + fermionsOut == 1)\n",
-    "        println(\"Average over $fermionsN fermion perturbations and sum over $photonsN photon perturbations for diagrams of form: $(trace)$(subtrace)\")\n",
-    "        return nothing\n",
+    "end\n",
+    "\n",
+    "copy(fd::FeynmanDiagram) = FeynmanDiagram(deepcopy(fd.vertices), deepcopy(fd.ties), deepcopy(fd.particles))\n",
+    "\n",
+    "function id_for_type(d::FeynmanDiagram, t::Type{<:QEDParticle})\n",
+    "    id = 1\n",
+    "    for l in 0:length(d.vertices)\n",
+    "        particles = get_particles(d, l)\n",
+    "        for p in particles\n",
+    "            if (p.particle <: t)\n",
+    "                id = max(id, p.id + 1)\n",
+    "            end\n",
+    "        end\n",
+    "    end\n",
+    "    return id\n",
+    "end\n",
+    "\n",
+    "function can_apply_vertex(particles::Vector{FeynmanParticle}, vertex::FeynmanVertex)\n",
+    "    return vertex.in1 in particles && vertex.in2 in particles && !(vertex.out in particles)\n",
+    "end\n",
+    "\n",
+    "function apply_vertex!(particles::Vector{FeynmanParticle}, vertex::FeynmanVertex)\n",
+    "    @assert can_apply_vertex(particles, vertex)\n",
+    "    length_before = length(particles)\n",
+    "    filter!(x -> x != vertex.in1 && x != vertex.in2, particles)\n",
+    "    push!(particles, vertex.out)\n",
+    "    @assert length(particles) == length_before - 1\n",
+    "    return nothing\n",
+    "end\n",
+    "\n",
+    "function can_apply_tie(particles::Vector{FeynmanParticle}, tie::FeynmanTie)\n",
+    "    return tie.in1 in particles && tie.in2 in particles\n",
+    "end\n",
+    "\n",
+    "function apply_tie!(particles::Vector{FeynmanParticle}, tie::FeynmanTie)\n",
+    "    @assert can_apply_tie(particles, tie)\n",
+    "    @assert can_tie(tie.in1.particle, tie.in2.particle)\n",
+    "    length_before = length(particles)\n",
+    "    filter!(x -> x != tie.in1 && x != tie.in2, particles)\n",
+    "    @assert length(particles) == length_before - 2\n",
+    "    return nothing\n",
+    "end\n",
+    "\n",
+    "function get_particles(fd::FeynmanDiagram, level::Int = -1)\n",
+    "    if level == -1\n",
+    "        level = length(fd.vertices)\n",
+    "    end\n",
+    "\n",
+    "    working_particles = copy(fd.particles)\n",
+    "    for l in 1:length(fd.vertices)\n",
+    "        if l > level\n",
+    "            break\n",
+    "        end\n",
+    "        for v in fd.vertices[l]\n",
+    "            apply_vertex!(working_particles, v)\n",
+    "        end\n",
+    "        for t in fd.ties[l]\n",
+    "            apply_tie!(working_particles, t)\n",
+    "        end\n",
+    "    end\n",
+    "\n",
+    "    return working_particles\n",
+    "end\n",
+    "\n",
+    "function add_vertex!(fd::FeynmanDiagram, vertex::FeynmanVertex)\n",
+    "    for i in eachindex(fd.vertices)\n",
+    "        if (can_apply_vertex(get_particles(fd, i - 1), vertex))\n",
+    "            push!(fd.vertices[i], vertex)\n",
+    "            return nothing\n",
+    "        end\n",
+    "    end\n",
+    "\n",
+    "    if !can_apply_vertex(get_particles(fd), vertex)\n",
+    "        @assert false \"Can't add vertex $vertex to diagram\"\n",
+    "    end\n",
+    "\n",
+    "    push!(fd.vertices, Set{FeynmanVertex}())\n",
+    "    push!(fd.vertices[end], vertex)\n",
+    "    push!(fd.ties, Set{FeynmanTie}())\n",
+    "    return nothing\n",
+    "end\n",
+    "\n",
+    "function add_vertex(fd::FeynmanDiagram, vertex::FeynmanVertex)\n",
+    "    newfd = copy(fd)\n",
+    "    add_vertex!(newfd, vertex)\n",
+    "    return newfd\n",
+    "end\n",
+    "\n",
+    "function add_tie!(fd::FeynmanDiagram, tie::FeynmanTie)\n",
+    "    for i in eachindex(fd.ties)\n",
+    "        if (can_apply_tie(get_particles(fd, i - 1), tie))\n",
+    "            push!(fd.ties[i], tie)\n",
+    "            return nothing\n",
+    "        end\n",
+    "    end\n",
+    "\n",
+    "    if !can_apply_tie(get_particles(fd), tie)\n",
+    "        @assert false \"Can't add tie $tie to diagram\"\n",
+    "    end\n",
+    "\n",
+    "    push!(fd.ties, Set{FeynmanTie}())\n",
+    "    push!(fd.ties[end], tie)\n",
+    "    push!(fd.vertices, Set{FeynmanVertex}())\n",
+    "    return nothing\n",
+    "end\n",
+    "\n",
+    "function add_tie(fd::FeynmanDiagram, tie::FeynmanTie)\n",
+    "    newfd = copy(fd)\n",
+    "    add_tie!(newfd, tie)\n",
+    "    return newfd\n",
+    "end\n",
+    "\n",
+    "# whether the diagram is connected and fully constructed\n",
+    "function isvalid(fd::FeynmanDiagram)\n",
+    "    @assert length(fd.vertices) == length(fd.ties)\n",
+    "\n",
+    "    if get_particles(fd) != []\n",
+    "        return false\n",
+    "    end\n",
+    "\n",
+    "    return true\n",
+    "end\n",
+    "\n",
+    "\n",
+    "function possible_vertices(fd::FeynmanDiagram)\n",
+    "    possibilities = Vector{FeynmanVertex}()\n",
+    "    \n",
+    "    particles = get_particles(fd)\n",
+    "    for i in 1:length(particles)\n",
+    "        for j in i+1:length(particles)\n",
+    "            p1 = particles[i]\n",
+    "            p2 = particles[j]\n",
+    "            if (caninteract(p1.particle, p2.particle))\n",
+    "                interaction_res = propagation_result(interaction_result(p1.particle, p2.particle))\n",
+    "                v = FeynmanVertex(p1, p2, FeynmanParticle(interaction_res, id_for_type(fd, interaction_res)))\n",
+    "                @assert !(v.out in particles) \"$v is in $fd\"\n",
+    "                @assert can_apply_vertex(particles, v)\n",
+    "                push!(possibilities, v)\n",
+    "            end\n",
+    "        end\n",
+    "    end\n",
+    "\n",
+    "    if (length(possibilities) == 0)\n",
+    "        #println(\"No vertices found for $particles\")\n",
+    "    end\n",
+    "\n",
+    "    return possibilities\n",
+    "end\n",
+    "\n",
+    "function can_tie(p1::Type, p2::Type)\n",
+    "    if p1 == propagation_result(p2)\n",
+    "        return true\n",
+    "    end\n",
+    "    if (p1 <: PhotonStateful && p2 <: PhotonStateful)\n",
+    "        return true\n",
+    "    end\n",
+    "    return false\n",
+    "end\n",
+    "\n",
+    "function possible_ties(fd::FeynmanDiagram)\n",
+    "    possibilities = Vector{FeynmanTie}()\n",
+    "\n",
+    "    particles = get_particles(fd)\n",
+    "\n",
+    "    for i in 1:length(particles)\n",
+    "        p1 = particles[i]\n",
+    "        if p1 in fd.particles\n",
+    "            continue\n",
+    "        end\n",
+    "        for j in i+1:length(particles)\n",
+    "            p2 = particles[j]\n",
+    "            if p2 in fd.particles\n",
+    "                continue\n",
+    "            end\n",
+    "            if (can_tie(p1.particle, p2.particle))\n",
+    "                t = FeynmanTie(p1, p2)\n",
+    "                @assert can_apply_tie(particles, t)\n",
+    "                push!(possibilities, t)\n",
+    "            end\n",
+    "        end\n",
     "    end\n",
     "    \n",
-    "    #println(\"trace: $trace  photons: $(photons)  fermions: $(fermionsIn) / $fermionsOut  antifermions: $antifermionsIn / $antifermionsOut\")\n",
+    "    return possibilities\n",
+    "end\n",
     "\n",
-    "    if (photons != 0)\n",
-    "        if (fermionsIn != 0)\n",
-    "            # photon + fermionin\n",
-    "            gen_diagram_helper(photons - 1, fermionsIn, fermionsOut, antifermionsIn, antifermionsOut, trace * \"ph + f(in) = f(in) -- \", fermionsN * fermionsIn, photonsN * photons)\n",
-    "        end\n",
-    "        if (fermionsOut != 0)\n",
-    "            # photon + fermionout\n",
-    "            gen_diagram_helper(photons - 1, fermionsIn + 1, fermionsOut, antifermionsIn, antifermionsOut, trace * \"ph + f(out) = f(out) -- \", fermionsN * fermionsOut, photonsN * photons)\n",
+    "function gen_diagrams(fd::FeynmanDiagram)\n",
+    "    working = Deque{FeynmanDiagram}()\n",
+    "    results = Set{FeynmanDiagram}()\n",
+    "\n",
+    "    push!(working, fd)\n",
+    "\n",
+    "    while !isempty(working)\n",
+    "        d = pop!(working)\n",
+    "        if (isvalid(d))\n",
+    "            push!(results, d)\n",
+    "            continue\n",
     "        end\n",
     "\n",
-    "        if (antifermionsIn != 0)\n",
-    "            # photon + antifermionin\n",
-    "            gen_diagram_helper(photons - 1, fermionsIn, fermionsOut, antifermionsIn, antifermionsOut, trace * \"ph + af(in) = af(in) -- \", fermionsN * antifermionsIn, photonsN * photons)\n",
+    "        possibilities = possible_vertices(d)\n",
+    "        for v in possibilities\n",
+    "            newfd = add_vertex(d, v)\n",
+    "            push!(working, newfd)\n",
     "        end\n",
-    "        if (antifermionsOut != 0)\n",
-    "            # photon + antifermionOut\n",
-    "            gen_diagram_helper(photons - 1, fermionsIn, fermionsOut, antifermionsIn, antifermionsOut, trace * \"ph + af(out) = af(out) -- \", fermionsN * antifermionsOut, photonsN * photons)\n",
-    "        end\n",
-    "\n",
-    "        if fermionsIn + fermionsOut + antifermionsIn + antifermionsOut != 0\n",
-    "            return nothing\n",
+    "        possibilities = possible_ties(d)\n",
+    "        for t in possibilities\n",
+    "            push!(working, add_tie(d, t))\n",
     "        end\n",
     "    end\n",
     "\n",
-    "    if (fermionsIn != 0)\n",
-    "        if (fermionsOut != 0)\n",
-    "            # fermionIn + fermionOut\n",
-    "            gen_diagram_helper(photons + 1, fermionsIn - 1, fermionsOut - 1, antifermionsIn, antifermionsOut, trace * \"f(in) + f(out) = ph -- \", fermionsN * fermionsIn * fermionsOut, photonsN)\n",
-    "            return nothing\n",
-    "        end\n",
-    "        if (antifermionsIn != 0)\n",
-    "            #fermionIn + antfermionIn\n",
-    "            gen_diagram_helper(photons + 1, fermionsIn - 1, fermionsOut, antifermionsIn - 1, antifermionsOut, trace * \"f(in) + af(in) = ph -- \", fermionsN * fermionsIn * antifermionsIn, photonsN)\n",
-    "            return nothing\n",
-    "        end\n",
-    "        #fermionIn + antifermionOut is impossible\n",
-    "    end\n",
+    "    return results\n",
+    "end\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 84,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Diagram 1: Initial Particles: [k_in_1, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [e_out_1 + k_out_1 -> e_out_2, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_2]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [k_in_1, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_1 -> e_out_2, e_in_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_2]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Compton Scattering\n",
+    "fd = FeynmanDiagram(parse_process(\"ke->ke\", QEDModel()))\n",
     "\n",
-    "    if (fermionsOut != 0)\n",
-    "        if (antifermionsOut != 0)\n",
-    "            #fermionOut + antifermionOut\n",
-    "            gen_diagram_helper(photons + 1, fermionsIn, fermionsOut - 1, antifermionsIn, antifermionsOut - 1, trace * \"f(out) + af(out) = ph -- \", fermionsN * fermionsOut * antifermionsOut, photonsN)\n",
-    "            return nothing\n",
-    "        end\n",
-    "        #fermionOut + antifermionIn is impossible\n",
-    "    end\n",
+    "diagrams = gen_diagrams(fd)\n",
     "\n",
-    "    if (antifermionsIn != 0)\n",
-    "        if (antifermionsOut != 0)\n",
-    "            #antifermionIn + antifermionOut\n",
-    "            gen_diagram_helper(photons + 1, fermionsIn, fermionsOut, antifermionsIn - 1, antifermionsOut - 1, trace * \"af(in) + af(out) = ph -- \", fermionsN * antifermionsIn * antifermionsOut, photonsN)\n",
-    "            return nothing\n",
-    "        end\n",
-    "    end\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 101,
    "metadata": {},
    "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 1354 samples with 1 evaluation.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m3.042 ms\u001b[22m\u001b[39m … \u001b[35m  6.234 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 34.23%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m3.667 ms               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m3.692 ms\u001b[22m\u001b[39m ± \u001b[32m577.114 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m2.87% ±  7.97%\n",
+       "\n",
+       "  \u001b[39m▂\u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[34m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\u001b[39m▅\u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m▅\u001b[39m▄\u001b[39m▆\u001b[39m▇\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▅\u001b[39m▄\u001b[34m▄\u001b[39m\u001b[39m▅\u001b[39m▄\u001b[39m▅\u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▅\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▃\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▃\u001b[39m▃\u001b[39m▂\u001b[39m▂\u001b[39m \u001b[39m▃\n",
+       "  3.04 ms\u001b[90m         Histogram: frequency by time\u001b[39m        5.72 ms \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m2.54 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m24047\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Average over 1 fermion perturbations and sum over 2 photon perturbations for diagrams of form: ph + f(in) = f(in) -- ph + f(in) + f(out)\n"
+      "Initial Diagram: Initial Particles: [k_in_1, k_in_2, e_in_1, e_out_1, k_out_1]\n",
+      "\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Found 12 Diagrams for 2-Photon Compton\n",
+      "Diagram 1: Initial Particles: [k_in_1, k_in_2, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_2 + e_in_1 -> e_in_2, e_out_1 + k_out_1 -> e_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [k_in_1 + e_in_2 -> e_in_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_3]\n",
+      "\n"
      ]
     }
    ],
    "source": [
-    "# Compton:\n",
-    "gen_diagram_helper(2, 1, 1, 0, 0)"
+    "# 2-Photon Compton Scattering\n",
+    "two_k_compton = FeynmanDiagram(parse_process(\"kke->ke\", QEDModel()))\n",
+    "\n",
+    "display(@benchmark gen_diagrams(two_k_compton))\n",
+    "diagrams = gen_diagrams(two_k_compton)\n",
+    "\n",
+    "println(\"Found $(length(diagrams)) Diagrams for 2-Photon Compton\")\n",
+    "println(\"Diagram 1: $(first(diagrams))\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 100,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 124 samples with 1 evaluation.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m35.054 ms\u001b[22m\u001b[39m … \u001b[35m48.433 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 5.92%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m40.638 ms              \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m3.66%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m40.457 ms\u001b[22m\u001b[39m ± \u001b[32m 2.120 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m3.49% ± 2.97%\n",
+       "\n",
+       "  \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▄\u001b[39m \u001b[39m▅\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▁\u001b[32m▁\u001b[39m\u001b[34m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m▄\u001b[39m█\u001b[39m█\u001b[39m▄\u001b[39m▅\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n",
+       "  \u001b[39m▃\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▅\u001b[39m▆\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▆\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m█\u001b[39m▃\u001b[39m█\u001b[39m▆\u001b[39m▁\u001b[39m▅\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[34m▅\u001b[39m\u001b[39m▆\u001b[39m▅\u001b[39m▃\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▃\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▃\u001b[39m \u001b[39m▃\n",
+       "  35.1 ms\u001b[90m         Histogram: frequency by time\u001b[39m        45.9 ms \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m26.97 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m262392\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Initial Diagram: Initial Particles: [k_in_1, k_in_2, e_in_1, e_out_1, k_out_1]\n",
+      "\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Found 72 Diagrams for 3-Photon Compton\n",
+      "Diagram 1: Initial Particles: [k_in_1, k_in_2, k_in_3, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_2 + e_out_1 -> e_out_2, k_in_3 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [k_out_1 + e_in_2 -> e_in_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: [k_in_1 + e_in_3 -> e_in_4]\n",
+      "  Ties: []\n",
+      "Virtuality Level 4:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_4]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# 3-Photon Compton Scattering\n",
+    "three_k_compton = FeynmanDiagram(parse_process(\"kkke->ke\", QEDModel()))\n",
+    "\n",
+    "display(@benchmark gen_diagrams(three_k_compton))\n",
+    "diagrams = gen_diagrams(three_k_compton)\n",
+    "\n",
+    "println(\"Found $(length(diagrams)) Diagrams for 3-Photon Compton\")\n",
+    "println(\"Diagram 1: $(first(diagrams))\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 103,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10 samples with 1 evaluation.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m520.818 ms\u001b[22m\u001b[39m … \u001b[35m588.348 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m3.87% … 6.33%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m553.318 ms               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m4.65%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m556.689 ms\u001b[22m\u001b[39m ± \u001b[32m 19.426 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m4.75% ± 1.01%\n",
+       "\n",
+       "  \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m█\u001b[39m \u001b[39m█\u001b[34m█\u001b[39m\u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[39m \n",
+       "  \u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m█\u001b[39m▁\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n",
+       "  521 ms\u001b[90m           Histogram: frequency by time\u001b[39m          588 ms \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m336.89 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m3305766\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Found 480 Diagrams for 4-Photon Compton\n",
+      "Diagram 1: Initial Particles: [k_in_1, k_in_2, k_in_3, k_in_4, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_1 -> e_out_2, e_in_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [k_in_2 + e_out_2 -> e_out_3, k_in_3 + e_in_2 -> e_in_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: [k_in_4 + e_out_3 -> e_out_4]\n",
+      "  Ties: []\n",
+      "Virtuality Level 4:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_3 -- e_out_4]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# 4-Photon Compton Scattering\n",
+    "four_k_compton = FeynmanDiagram(parse_process(\"kkkke->ke\", QEDModel()))\n",
+    "\n",
+    "display(@benchmark gen_diagrams(four_k_compton))\n",
+    "diagrams = gen_diagrams(four_k_compton)\n",
+    "\n",
+    "println(\"Found $(length(diagrams)) Diagrams for 4-Photon Compton\")\n",
+    "println(\"Diagram 1: $(first(diagrams))\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 104,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 1 sample with 1 evaluation.\n",
+       " Single result which took \u001b[34m8.170 s\u001b[39m (7.15% GC) to evaluate,\n",
+       " with a memory estimate of \u001b[33m4.69 GiB\u001b[39m, over \u001b[33m47449329\u001b[39m allocations."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Found 3600 Diagrams for 5-Photon Compton\n",
+      "Diagram 1: Initial Particles: [k_in_1, k_in_2, k_in_3, k_in_4, k_in_5, e_in_1, e_out_1, k_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_2 + e_out_1 -> e_out_2, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [k_in_4 + e_in_2 -> e_in_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: [k_in_3 + e_in_3 -> e_in_4]\n",
+      "  Ties: []\n",
+      "Virtuality Level 4:\n",
+      "  Vertices: [k_out_1 + e_in_4 -> e_in_5]\n",
+      "  Ties: []\n",
+      "Virtuality Level 5:\n",
+      "  Vertices: [k_in_5 + e_in_5 -> e_in_6]\n",
+      "  Ties: []\n",
+      "Virtuality Level 6:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_6]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# 5-Photon Compton Scattering\n",
+    "five_k_compton = FeynmanDiagram(parse_process(\"kkkkke->ke\", QEDModel()))\n",
+    "\n",
+    "display(@benchmark gen_diagrams(five_k_compton))\n",
+    "diagrams = gen_diagrams(five_k_compton)\n",
+    "\n",
+    "println(\"Found $(length(diagrams)) Diagrams for 5-Photon Compton\")\n",
+    "println(\"Diagram 1: $(first(diagrams))\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 87,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Average over 1 fermion perturbations and sum over 2 photon perturbations for diagrams of form: ph + f(in) = f(in) -- ph + f(in) + af(in)\n",
-      "Average over 1 fermion perturbations and sum over 2 photon perturbations for diagrams of form: ph + af(in) = af(in) -- ph + f(in) + af(in)\n"
+      "Diagram 1: Initial Particles: [p_in_1, e_in_1, p_out_1, e_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_in_1 + p_out_1 -> k_out_1, e_in_1 + e_out_1 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [p_in_1, e_in_1, p_out_1, e_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_in_1 + e_in_1 -> k_out_1, p_out_1 + e_out_1 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Bhabha Scattering\n",
+    "fd = FeynmanDiagram(parse_process(\"ep->ep\", QEDModel()))\n",
+    "\n",
+    "diagrams = gen_diagrams(fd)\n",
+    "\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 88,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Diagram 1: Initial Particles: [e_in_1, e_in_2, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [e_in_1 + e_out_2 -> k_out_1, e_in_2 + e_out_1 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [e_in_1, e_in_2, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [e_in_2 + e_out_2 -> k_out_2, e_in_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_2 -- k_out_1]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Moller Scattering\n",
+    "fd = FeynmanDiagram(parse_process(\"ee->ee\", QEDModel()))\n",
+    "\n",
+    "diagrams = gen_diagrams(fd)\n",
+    "\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Diagram 1: Initial Particles: [p_in_1, e_in_1, k_out_1, k_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_in_1 + k_out_2 -> e_out_1, e_in_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_1 -- e_in_2]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [p_in_1, e_in_1, k_out_1, k_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_in_1 + k_out_1 -> e_out_1, e_in_1 + k_out_2 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_1 -- e_in_2]\n",
+      "\n"
      ]
     }
    ],
    "source": [
     "# Pair annihilation\n",
-    "gen_diagram_helper(2, 1, 0, 1, 0)"
+    "fd = FeynmanDiagram(parse_process(\"ep->kk\", QEDModel()))\n",
+    "\n",
+    "diagrams = gen_diagrams(fd)\n",
+    "\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
+    "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 90,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Average over 1 fermion perturbations and sum over 2 photon perturbations for diagrams of form: ph + af(out) = af(out) -- ph + f(out) + af(out)\n"
+      "Diagram 1: Initial Particles: [k_in_1, k_in_2, p_out_1, e_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + p_out_1 -> e_in_1, k_in_2 + e_out_1 -> e_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_1 -- e_out_2]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [k_in_1, k_in_2, p_out_1, e_out_1]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_1 -> e_out_2, k_in_2 + p_out_1 -> e_in_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_2 -- e_in_1]\n",
+      "\n"
      ]
     }
    ],
    "source": [
     "# Pair production\n",
-    "gen_diagram_helper(2, 0, 1, 0, 1)"
+    "fd = FeynmanDiagram(parse_process(\"kk->pe\", QEDModel()))\n",
+    "\n",
+    "diagrams = gen_diagrams(fd)\n",
+    "\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
+    "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 91,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Average over 2 fermion perturbations and sum over 1 photon perturbations for diagrams of form: ph + f(in) = f(in) -- f(in) + f(out) = ph -- ph + f(out) + af(out)\n",
-      "Average over 2 fermion perturbations and sum over 1 photon perturbations for diagrams of form: ph + af(out) = af(out) -- f(in) + f(out) = ph -- ph + f(out) + af(out)\n"
+      "Found 16 diagrams:\n",
+      "Diagram 1: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + p_out_1 -> e_in_2, e_in_1 + e_out_2 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_1 + k_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_2 -- e_out_3]\n",
+      "\n",
+      "Diagram 2: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_2 -> e_out_3, e_in_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [p_out_1 + e_out_3 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 3: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + p_out_1 -> e_in_2, e_in_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_2 + e_in_2 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 4: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [e_in_1 + e_out_2 -> k_out_1, k_in_1 + e_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [p_out_1 + e_out_3 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 5: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + p_out_1 -> e_in_2, e_in_1 + e_out_2 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_1 + e_in_2 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 6: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + p_out_1 -> e_in_2, e_in_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_2 + k_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_2 -- e_out_3]\n",
+      "\n",
+      "Diagram 7: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_out_1 + e_out_2 -> k_out_1, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_1 + k_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_2 -- e_out_3]\n",
+      "\n",
+      "Diagram 8: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_out_1 + e_out_1 -> k_out_1, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_2 + e_in_2 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 9: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_out_1 + e_out_2 -> k_out_1, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_1 + e_in_2 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_1 -- k_out_2]\n",
+      "\n",
+      "Diagram 10: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_1 -> e_out_3, p_out_1 + e_out_2 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_in_1 + e_out_3 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_2 -- k_out_1]\n",
+      "\n",
+      "Diagram 11: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_2 -> e_out_3, p_out_1 + e_out_1 -> k_out_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_in_1 + e_out_3 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [k_out_2 -- k_out_1]\n",
+      "\n",
+      "Diagram 12: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_2 -> e_out_3, e_in_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [p_out_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_3 -- e_in_2]\n",
+      "\n",
+      "Diagram 13: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [p_out_1 + e_out_1 -> k_out_1, k_in_1 + e_in_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_out_2 + k_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_in_2 -- e_out_3]\n",
+      "\n",
+      "Diagram 14: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_1 -> e_out_3, p_out_1 + e_out_2 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_in_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_3 -- e_in_2]\n",
+      "\n",
+      "Diagram 15: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [e_in_1 + e_out_2 -> k_out_1, k_in_1 + e_out_1 -> e_out_3]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [p_out_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_3 -- e_in_2]\n",
+      "\n",
+      "Diagram 16: Initial Particles: [k_in_1, e_in_1, p_out_1, e_out_1, e_out_2]\n",
+      "Virtuality Level 1:\n",
+      "  Vertices: [k_in_1 + e_out_2 -> e_out_3, p_out_1 + e_out_1 -> k_out_1]\n",
+      "  Ties: []\n",
+      "Virtuality Level 2:\n",
+      "  Vertices: [e_in_1 + k_out_1 -> e_in_2]\n",
+      "  Ties: []\n",
+      "Virtuality Level 3:\n",
+      "  Vertices: []\n",
+      "  Ties: [e_out_3 -- e_in_2]\n",
+      "\n"
      ]
     }
    ],
    "source": [
     "# Trident\n",
-    "gen_diagram_helper(1, 1, 2, 0, 1)"
+    "fd = FeynmanDiagram(parse_process(\"ke->epe\", QEDModel()))\n",
+    "\n",
+    "diagrams = gen_diagrams(fd)\n",
+    "\n",
+    "println(\"Found $(length(diagrams)) diagrams:\")\n",
+    "c = 1\n",
+    "for d in diagrams\n",
+    "    println(\"Diagram $c: $d\")\n",
+    "    c += 1\n",
+    "end"
    ]
   }
  ],
diff --git a/src/models/qed/particle.jl b/src/models/qed/particle.jl
index 1b56175..749d98a 100644
--- a/src/models/qed/particle.jl
+++ b/src/models/qed/particle.jl
@@ -169,6 +169,9 @@ end
 String(::Type{Incoming}) = "Incoming"
 String(::Type{Outgoing}) = "Outgoing"
 
+String(::Incoming) = "in"
+String(::Outgoing) = "out"
+
 function String(::Type{<:PhotonStateful})
     return "k"
 end