data/evaluate.jl

@@ -21,18 +21,10 @@ processes = [
     "ABC Process: 'AB->ABBBBB'",
 ]
 
-
 function proc_to_n(str::AbstractString)
     parts = split(str, "'")
-
-    infix = parts[2]
-
-    parts = split(infix, "->")
-
-    suffix = parts[2]
-
-    k_count = count(c -> c == 'k', suffix)
-
+    parts = split(parts[2], "->")
+    k_count = count(c -> c == 'k', parts[2])
     return k_count
 end
 
@@ -43,7 +35,6 @@ function beautify_title(str::AbstractString)
     infix = parts[2]
     sufsuffix = parts[3]
 
-
     parts = split(infix, "->")
 
     prefix = parts[1]

data/evaluate_gen.jl

@@ -0,0 +1,148 @@
+using CSV
+using DataFrames
+using Plots
+using StatsPlots
+using LaTeXStrings
+
+if (length(ARGS) < 1)
+    println("Please use with \"input_file.csv\"")
+end
+
+function proc_to_n(str::AbstractString)
+    parts = split(str, "->")
+    k_count = count(c -> c == 'k', parts[2])
+    return k_count
+end
+
+input_file = ARGS[1]
+df = CSV.read(input_file, DataFrame)
+
+# plotting with process size as x axis
+THREADS = [4]
+
+for threads in THREADS
+    title_string = "n-photon Compton diagram generation, $threads threads"
+
+    df_filt = filter(:cpu_threads => x -> x == threads, df)
+    df_filt = filter(:process_name => x -> proc_to_n(x) >= 1, df_filt)
+    df_filt.graph_gen_mean = @. df_filt.graph_gen_mean / 1e9
+    df_filt.graph_gen_std = @. df_filt.graph_gen_std / 1e9
+
+    df_filt.process_size = @. proc_to_n(df_filt.process_name)
+
+    @df df_filt scatter(
+        :process_size,
+        :graph_gen_mean,
+        yerror = :graph_gen_std,
+        label = "graph generation time",
+        markersize = 7,
+    )
+
+    plot!(
+        title = title_string,
+        yscale = :log10,
+        legend = :outerbottom,
+        minorgrid = true,
+        xticks = :process_size,
+        yticks = [1e-3, 1e-2, 1e-1, 1e-0, 1e1],
+        xgrid = false,
+        xminorticks = false,
+        legendcolumns = 1,
+        legend_font_pointsize = 10,
+        size = (800, 600),
+        ylabel = "time (s)",
+        xlabel = "process size (#)",
+    )
+
+    savefig("gen_times_$(threads)_threads.pdf")
+
+    # graph size
+    title_string = "n-photon Compton unreduced graph size"
+
+    @df df_filt scatter(:process_size, :graph_nodes, label = "nodes", markershape = :circle)
+    @df df_filt scatter!(:process_size, :graph_edges, label = "edges", markershape = :square)
+    @df df_filt scatter!(:process_size, :graph_u_nodes, label = "U-nodes", markershape = :star)
+    @df df_filt scatter!(:process_size, :graph_v_nodes, label = "V-nodes", markershape = :utriangle)
+    @df df_filt scatter!(:process_size[2:end], :graph_s1_nodes[2:end], label = "S1-nodes", markershape = :x)
+    @df df_filt scatter!(:process_size, :graph_s2_nodes, label = "S2-nodes", markershape = :diamond)
+
+    plot!(
+        title = title_string,
+        yscale = :log10,
+        legend = :outerbottom,
+        yminorgrid = true,
+        xticks = :process_size,
+        yticks = [1e1, 1e3, 1e5, 1e7],
+        xgrid = false,
+        xminorticks = false,
+        legendcolumns = 2,
+        legend_font_pointsize = 10,
+        size = (800, 600),
+        ylabel = "(#)",
+        xlabel = "process size (#)",
+    )
+
+    savefig("compton_graph_size_unreduced.pdf")
+
+
+    # graph size
+    title_string = "n-photon Compton reduced graph size"
+
+    @df df_filt scatter(:process_size, :graph_nodes_reduced, label = "nodes", markershape = :circle)
+    @df df_filt scatter!(:process_size, :graph_edges_reduced, label = "edges", markershape = :square)
+    @df df_filt scatter!(:process_size, :graph_u_nodes_reduced, label = "U-nodes", markershape = :star)
+    @df df_filt scatter!(:process_size, :graph_v_nodes_reduced, label = "V-nodes", markershape = :utriangle)
+    @df df_filt scatter!(:process_size[2:end], :graph_s1_nodes_reduced[2:end], label = "S1-nodes", markershape = :x)
+    @df df_filt scatter!(:process_size, :graph_s2_nodes_reduced, label = "S2-nodes", markershape = :diamond)
+
+    plot!(
+        title = title_string,
+        yscale = :log10,
+        legend = :outerbottom,
+        yminorgrid = true,
+        xticks = :process_size,
+        yticks = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6],
+        xgrid = false,
+        xminorticks = false,
+        legendcolumns = 2,
+        legend_font_pointsize = 10,
+        size = (800, 600),
+        ylabel = "(#)",
+        xlabel = "process size (#)",
+    )
+
+    savefig("compton_graph_size_reduced.pdf")
+
+
+    # graph size versus
+    title_string = "n-photon Compton graph sizes"
+
+    @df df_filt scatter(:process_size, :graph_nodes, label = "nodes", markershape = :circle)
+    @df df_filt scatter!(:process_size, :graph_edges, label = "edges", markershape = :square)
+    @df df_filt scatter!(:process_size, :graph_nodes_reduced, label = "nodes (after reduction)", markershape = :star)
+    @df df_filt scatter!(
+        :process_size,
+        :graph_edges_reduced,
+        label = "edges (after reduction)",
+        markershape = :utriangle,
+    )
+
+    plot!(
+        title = title_string,
+        yscale = :log10,
+        legend = :outerbottom,
+        yminorgrid = true,
+        xticks = :process_size,
+        yticks = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6],
+        xgrid = false,
+        xminorticks = false,
+        legendcolumns = 2,
+        legend_font_pointsize = 10,
+        size = (800, 600),
+        ylabel = "(#)",
+        xlabel = "process size (#)",
+    )
+
+    savefig("compton_graph_size_versus.pdf")
+
+end

data/qed_gen_results.csv

@@ -0,0 +1,9 @@
+process_name,cpu_threads,graph_gen_samples,graph_gen_mean,graph_gen_std,graph_gen_median,graph_nodes,graph_data_nodes,graph_u_nodes,graph_v_nodes,graph_s1_nodes,graph_s2_nodes,graph_edges,graph_nodes_reduced,graph_data_nodes_reduced,graph_u_nodes_reduced,graph_v_nodes_reduced,graph_s1_nodes_reduced,graph_s2_nodes_reduced,graph_edges_reduced,graph_mem,graph_mem_reduced
+ke->ke,4,621,108190.55233494364,15933.10900718239,108119.0,26,15,4,4,0,2,29,26,15,4,4,0,2,29,23980.0,25292.0
+ke->kke,4,622,153201.41318327974,23413.42529125677,151041.5,77,41,5,18,6,6,101,59,32,5,12,3,6,77,35423.0,47607.0
+ke->kkke,4,598,404339.652173913,35709.9373272436,398624.0,356,181,6,96,48,24,493,148,77,6,32,8,24,221,95474.0,156874.0
+ke->kkkke,4,613,1.9826555334420882e6,120582.30781418575,1.972286e6,2183,1095,7,600,360,120,3015,543,275,7,110,30,120,885,470405.0,788269.0
+ke->kkkkke,4,565,1.3094592946902655e7,626078.9227119224,1.3065283e7,15866,7937,8,4320,2880,720,21617,2234,1121,8,312,72,720,3977,3.619736e6,5.464312e6
+ke->kkkkkke,4,386,1.46437015619171e8,3.4493988193742386e6,1.46310623e8,131069,65539,9,35280,25200,5040,176419,13441,6725,9,1358,308,5040,24869,2.7220139e7,4.1642163e7
+ke->kkkkkkke,4,95,1.4519048337473683e9,3.4428550325367525e7,1.44881422e9,1209632,604821,10,322560,241920,40320,1612821,90592,45301,10,4160,800,40320,175381,2.47916222e8,3.6695295e8
+ke->kkkkkkkke,4,10,2.09393503407e10,1.2609409886873345e9,2.1513879476e10,12337955,6168983,11,3265920,2540160,362880,16329623,778859,389435,11,22338,4194,362880,1526945,2.636058833e9,3.753461609e9

examples/qed_gen_bench.jl

@@ -0,0 +1,106 @@
+using MetagraphOptimization
+using DataFrames
+using CSV
+using BenchmarkTools
+using StatsBase
+
+results_filename = "qed_gen_results.csv"
+
+df = DataFrame(
+    process_name = String[],
+    cpu_threads = Int[],
+    graph_gen_samples = Int[],
+    graph_gen_mean = Float64[],
+    graph_gen_std = Float64[],
+    graph_gen_median = Float64[],
+    graph_nodes = Int[],
+    graph_data_nodes = Int[],
+    graph_u_nodes = Int[],
+    graph_v_nodes = Int[],
+    graph_s1_nodes = Int[],
+    graph_s2_nodes = Int[],
+    graph_edges = Int[],
+    graph_nodes_reduced = Int[],
+    graph_data_nodes_reduced = Int[],
+    graph_u_nodes_reduced = Int[],
+    graph_v_nodes_reduced = Int[],
+    graph_s1_nodes_reduced = Int[],
+    graph_s2_nodes_reduced = Int[],
+    graph_edges_reduced = Int[],
+    graph_mem = Float64[],
+    graph_mem_reduced = Float64[],
+)
+
+# if they exist, read existing results and append new ones
+if isfile(results_filename)
+    df = CSV.read(results_filename, DataFrame)
+end
+
+processes = [
+    "ke->ke",
+    "ke->kke",
+    "ke->kkke",
+    "ke->kkkke",
+    "ke->kkkkke",
+    "ke->kkkkkke",
+    "ke->kkkkkkke",
+    "ke->kkkkkkkke",
+    #"ke->kkkkkkkkke",
+]
+
+function bench_process(process::AbstractString)
+    println("Benchmarking $process...")
+    model = QEDModel()
+
+    proc = parse_process(process, model)
+
+    gen_bench = @benchmark gen_graph($proc) gcsample = true seconds = 300
+
+    graph = gen_graph(proc)
+
+    props = GraphProperties(graph)
+    node_dict = countmap(typeof.(graph.nodes))
+    graph_size = Base.summarysize(graph)
+
+    optim = ReductionOptimizer()
+    optimize_to_fixpoint!(optim, graph)
+
+    props_reduced = GraphProperties(graph)
+    node_dict_reduced = countmap(typeof.(graph.nodes))
+    graph_size_reduced = Base.summarysize(graph)
+
+    push!(
+        df,
+        Dict(
+            :process_name => process,
+            :cpu_threads => Threads.nthreads(),
+            :graph_gen_samples => length(gen_bench.times),
+            :graph_gen_mean => mean(gen_bench.times),
+            :graph_gen_std => std(gen_bench.times),
+            :graph_gen_median => median(gen_bench.times),
+            :graph_nodes => props.noNodes,
+            :graph_data_nodes => get(node_dict, DataTaskNode{DataTask}, 0),
+            :graph_u_nodes => get(node_dict, ComputeTaskNode{ComputeTaskQED_U}, 0),
+            :graph_v_nodes => get(node_dict, ComputeTaskNode{ComputeTaskQED_V}, 0),
+            :graph_s1_nodes => get(node_dict, ComputeTaskNode{ComputeTaskQED_S1}, 0),
+            :graph_s2_nodes => get(node_dict, ComputeTaskNode{ComputeTaskQED_S2}, 0),
+            :graph_edges => props.noEdges,
+            :graph_nodes_reduced => props_reduced.noNodes,
+            :graph_data_nodes_reduced => get(node_dict_reduced, DataTaskNode{DataTask}, 0),
+            :graph_u_nodes_reduced => get(node_dict_reduced, ComputeTaskNode{ComputeTaskQED_U}, 0),
+            :graph_v_nodes_reduced => get(node_dict_reduced, ComputeTaskNode{ComputeTaskQED_V}, 0),
+            :graph_s1_nodes_reduced => get(node_dict_reduced, ComputeTaskNode{ComputeTaskQED_S1}, 0),
+            :graph_s2_nodes_reduced => get(node_dict_reduced, ComputeTaskNode{ComputeTaskQED_S2}, 0),
+            :graph_edges_reduced => props_reduced.noEdges,
+            :graph_mem => graph_size,
+            :graph_mem_reduced => graph_size_reduced,
+        ),
+    )
+    return nothing
+end
+
+for process in processes
+    bench_process(process)
+end
+
+CSV.write(results_filename, df)

src/models/qed/create.jl

@@ -114,12 +114,8 @@ function gen_graph(process_description::QEDProcessDescription)
         dataOutNodes[String(particle)] = data_out
     end
 
-    #dataOutBackup = copy(dataOutNodes)
-
+    # TODO: this should be parallelizable somewhat easily
     for diagram in diagrams
-        # the intermediate (virtual) particles change across
-        #dataOutNodes = copy(dataOutBackup)
-
         tie = diagram.tie[]
 
         # handle the vertices

src/models/qed/diagrams.jl

@@ -447,75 +447,88 @@ Returns true iff the given feynman diagram is an (empty) diagram of a compton pr
 function is_compton(fd::FeynmanDiagram)
     return fd.type_ids[FermionStateful{Incoming, SpinUp}] == 1 &&
            fd.type_ids[FermionStateful{Outgoing, SpinUp}] == 1 &&
-           fd.type_ids[PhotonStateful{Incoming, PolX}] == 1 &&
+           fd.type_ids[AntiFermionStateful{Incoming, SpinUp}] == 0 &&
+           fd.type_ids[AntiFermionStateful{Outgoing, SpinUp}] == 0 &&
+           fd.type_ids[PhotonStateful{Incoming, PolX}] >= 1 &&
            fd.type_ids[PhotonStateful{Outgoing, PolX}] >= 1
 end
 
 """
-    n_photon_compton_diagrams(n::Int)
+    gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)
 
-Special case diagram generation for n-Photon Compton processes, i.e., processes of the form ke->k^ne
+Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram for the given order and n input and m output photons.
 """
-function n_photon_compton_diagram(n::Int)
-    inEl = FeynmanParticle(FermionStateful{Incoming, SpinUp}, 1)
-    outEl = FeynmanParticle(FermionStateful{Outgoing, SpinUp}, 1)
+function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)
+    photons = vcat(
+        [FeynmanParticle(PhotonStateful{Incoming, PolX}, i) for i in 1:n],
+        [FeynmanParticle(PhotonStateful{Outgoing, PolX}, i) for i in 1:m],
+    )
 
-    photons = [FeynmanParticle(PhotonStateful{Outgoing, PolX}, i) for i in 1:n]
-    push!(photons, FeynmanParticle(PhotonStateful{Incoming, PolX}, 1))
+    new_diagram = FeynmanDiagram(
+        [],
+        missing,
+        [inFerm, outFerm, photons...],
+        Dict{Type, Int64}(
+            FermionStateful{Incoming, SpinUp} => 1,
+            FermionStateful{Outgoing, SpinUp} => 1,
+            PhotonStateful{Incoming, PolX} => n,
+            PhotonStateful{Outgoing, PolX} => m,
+        ),
+    )
 
-    perms = permutations([i for i in 1:length(photons)])
+    left_index = 1
+    right_index = length(order)
 
-    diagrams = Vector{FeynmanDiagram}()
-    for order in perms
-        new_diagram = FeynmanDiagram(
-            [],
-            missing,
-            [inEl, outEl, photons...],
-            Dict{Type, Int64}(
-                FermionStateful{Incoming, SpinUp} => 1,
-                FermionStateful{Outgoing, SpinUp} => 1,
-                PhotonStateful{Incoming, PolX} => 1,
-                PhotonStateful{Outgoing, PolX} => n,
-            ),
+    iterations = 1
+
+    while left_index <= right_index
+        # left side
+        v_left = FeynmanVertex(
+            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
+            photons[order[left_index]],
+            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
         )
+        left_index += 1
+        add_vertex!(new_diagram, v_left)
 
-        left_index = 1
-        right_index = length(order)
-
-        iterations = 1
-
-        while left_index <= right_index
-            # left side
-            v_left = FeynmanVertex(
-                FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
-                photons[order[left_index]],
-                FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
-            )
-            left_index += 1
-            add_vertex!(new_diagram, v_left)
-
-            if (left_index > right_index)
-                break
-            end
-
-            # right side
-            v_right = FeynmanVertex(
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
-                photons[order[right_index]],
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
-            )
-            right_index -= 1
-            add_vertex!(new_diagram, v_right)
-
-            iterations += 1
+        if (left_index > right_index)
+            break
         end
 
-        @assert possible_tie(new_diagram) !== missing
-        add_tie!(new_diagram, possible_tie(new_diagram))
-        push!(diagrams, new_diagram)
+        # right side
+        v_right = FeynmanVertex(
+            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
+            photons[order[right_index]],
+            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
+        )
+        right_index -= 1
+        add_vertex!(new_diagram, v_right)
+
+        iterations += 1
     end
 
-    return diagrams
+    @assert possible_tie(new_diagram) !== missing
+    add_tie!(new_diagram, possible_tie(new_diagram))
+    return return new_diagram
+end
+
+"""
+    gen_compton_diagrams(n::Int, m::Int)
+
+Special case diagram generation for Compton processes, i.e., processes of the form k^ne->k^me
+"""
+function gen_compton_diagrams(n::Int, m::Int)
+    inFerm = FeynmanParticle(FermionStateful{Incoming, SpinUp}, 1)
+    outFerm = FeynmanParticle(FermionStateful{Outgoing, SpinUp}, 1)
+
+    perms = [permutations([i for i in 1:(n + m)])...]
+
+    diagrams = [Vector{FeynmanDiagram}() for i in 1:nthreads()]
+    @threads for order in perms
+        push!(diagrams[threadid()], gen_compton_diagram_from_order(order, inFerm, outFerm, n, m))
+    end
+
+    return vcat(diagrams...)
 end
 
 """
@@ -525,7 +538,10 @@ From a given feynman diagram in its initial state, e.g. when created through the
 """
 function gen_diagrams(fd::FeynmanDiagram)
     if is_compton(fd)
-        return n_photon_compton_diagram(fd.type_ids[PhotonStateful{Outgoing, PolX}])
+        return gen_compton_diagrams(
+            fd.type_ids[PhotonStateful{Incoming, PolX}],
+            fd.type_ids[PhotonStateful{Outgoing, PolX}],
+        )
     end
 
     working = Set{FeynmanDiagram}()