Add results and evaluation

Congruent in photons example (#12 )
Now targeting the correct branches Co-authored-by: Rubydragon <anton.reinhard@proton.me> Reviewed-on: #12
2024-07-10 14:21:26 +02:00 · 2024-07-10 14:17:39 +02:00 · 2024-07-04 15:31:22 +02:00 · 2024-07-03 20:24:53 +02:00 · 2024-07-02 10:50:30 +02:00 · 2024-06-24 23:31:30 +02:00
46 changed files with 1362 additions and 953 deletions
--- a/.gitattributes
+++ b/.gitattributes
@ -1,3 +1,5 @@
 input/AB->ABBBBBBBBB.txt filter=lfs diff=lfs merge=lfs -text
 input/AB->ABBBBBBB.txt filter=lfs diff=lfs merge=lfs -text
-*.zip filter=lfs diff=lfs merge=lfs -text
+*.zip filter=lfs diff=lfs merge=lfs
+*.gif filter=lfs diff=lfs merge=lfs
+*.jld2 filter=lfs diff=lfs merge=lfs
--- a/Project.toml
+++ b/Project.toml
@ -14,16 +14,22 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 NumaAllocators = "21436f30-1b4a-4f08-87af-e26101bb5379"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
 QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"

 [extras]
 CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
+QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

 [targets]
-test = ["Test"]
+test = ["SafeTestsets", "Test", "QEDbase", "QEDcore", "QEDprocesses"]
--- a/examples/congruent_in_ph.jl
+++ b/examples/congruent_in_ph.jl
@ -0,0 +1,194 @@
+ENV["UCX_ERROR_SIGNALS"] = "SIGILL,SIGBUS,SIGFPE"
+
+using MetagraphOptimization
+using QEDbase
+using QEDcore
+using QEDprocesses
+using Random
+using UUIDs
+
+using CUDA
+
+using NamedDims
+using CSV
+using JLD2
+using FlexiMaps
+
+RNG = Random.MersenneTwister(123)
+
+function mock_machine()
+    return Machine(
+        [
+            MetagraphOptimization.NumaNode(
+                0,
+                1,
+                MetagraphOptimization.default_strategy(MetagraphOptimization.NumaNode),
+                -1.0,
+                UUIDs.uuid1(),
+            ),
+        ],
+        [-1.0;;],
+    )
+end
+
+function congruent_input_momenta(processDescription::GenericQEDProcess, omega::Number)
+    # generate an input sample for given e + nk -> e' + k' process, where the nk are equal
+    inputMasses = Vector{Float64}()
+    for particle in incoming_particles(processDescription)
+        push!(inputMasses, mass(particle))
+    end
+    outputMasses = Vector{Float64}()
+    for particle in outgoing_particles(processDescription)
+        push!(outputMasses, mass(particle))
+    end
+
+    initial_momenta = [
+        i == length(inputMasses) ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for
+        i in 1:length(inputMasses)
+    ]
+
+    ss = sqrt(sum(initial_momenta) * sum(initial_momenta))
+    final_momenta = MetagraphOptimization.generate_physical_massive_moms(RNG, ss, outputMasses)
+
+    return (tuple(initial_momenta...), tuple(final_momenta...))
+end
+
+# theta ∈ [0, 2π] and phi ∈ [0, 2π]
+function congruent_input_momenta_scenario_2(
+    processDescription::GenericQEDProcess,
+    omega::Number,
+    theta::Number,
+    phi::Number,
+)
+    # -------------
+    # same as above
+
+    # generate an input sample for given e + nk -> e' + k' process, where the nk are equal
+    inputMasses = Vector{Float64}()
+    for particle in incoming_particles(processDescription)
+        push!(inputMasses, mass(particle))
+    end
+    outputMasses = Vector{Float64}()
+    for particle in outgoing_particles(processDescription)
+        push!(outputMasses, mass(particle))
+    end
+
+    initial_momenta = [
+        i == length(inputMasses) ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for
+        i in 1:length(inputMasses)
+    ]
+
+    ss = sqrt(sum(initial_momenta) * sum(initial_momenta))
+
+    # up to here
+    # ----------
+
+    # now calculate the final_momenta from omega, cos_theta and phi
+    n = number_particles(processDescription, Incoming(), Photon())
+
+    cos_theta = cos(theta)
+    omega_prime = (n * omega) / (1 + n * omega * (1 - cos_theta))
+
+    k_prime =
+        omega_prime * SFourMomentum(1, sqrt(1 - cos_theta^2) * cos(phi), sqrt(1 - cos_theta^2) * sin(phi), cos_theta)
+    p_prime = sum(initial_momenta) - k_prime
+
+    final_momenta = (k_prime, p_prime)
+
+    return (tuple(initial_momenta...), tuple(final_momenta...))
+
+end
+
+function build_psp(processDescription::GenericQEDProcess, momenta)
+    return PhaseSpacePoint(
+        processDescription,
+        PerturbativeQED(),
+        PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
+        momenta[1],
+        momenta[2],
+    )
+end
+
+# hack to fix stacksize for threading
+with_stacksize(f, n) = fetch(schedule(Task(f, n)))
+
+# scenario 2
+N = 1024    # thetas
+M = 1024    # phis
+K = 64      # omegas
+
+thetas = collect(LinRange(0, 2π, N))
+phis = collect(LinRange(0, 2π, M))
+omegas = collect(maprange(log, 2e-2, 2e-7, K))
+
+for photons in 1:5
+    # temp process to generate momenta
+    println("Generating $(K*N*M) inputs for $photons photons (Scenario 2 grid walk)...")
+    temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp())
+
+    input_momenta =
+        Array{typeof(congruent_input_momenta_scenario_2(temp_process, omegas[1], thetas[1], phis[1]))}(undef, (K, N, M))
+
+    Threads.@threads for k in 1:K
+        Threads.@threads for i in 1:N
+            Threads.@threads for j in 1:M
+                input_momenta[k, i, j] = congruent_input_momenta_scenario_2(temp_process, omegas[k], thetas[i], phis[j])
+            end
+        end
+    end
+
+
+    cu_results = CuArray{Float64}(undef, size(input_momenta))
+    fill!(cu_results, 0.0)
+
+    i = 1
+    for (in_pol, in_spin, out_pol, out_spin) in
+        Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()])
+
+        print(
+            "[$i/16] Calculating for spin/pol config: $in_pol, $in_spin -> $out_pol, $out_spin... Preparing inputs... ",
+        )
+        process = parse_process("k"^photons * "e->ke", QEDModel(), in_pol, in_spin, out_pol, out_spin)
+
+        inputs = Array{typeof(build_psp(process, input_momenta[1, 1, 1]))}(undef, (K, N, M))
+        #println("input_momenta: $input_momenta")
+        Threads.@threads for k in 1:K
+            Threads.@threads for i in 1:N
+                Threads.@threads for j in 1:M
+                    inputs[k, i, j] = build_psp(process, input_momenta[k, i, j])
+                end
+            end
+        end
+        cu_inputs = CuArray(inputs)
+
+        print("Preparing graph... ")
+        graph = gen_graph(process)
+        optimize_to_fixpoint!(ReductionOptimizer(), graph)
+        print("Preparing function... ")
+        kernel! = get_cuda_kernel(graph, process, mock_machine())
+        #func = get_compute_function(graph, process, mock_machine())
+
+        print("Calculating... ")
+        ts = 32
+        bs = Int64(length(cu_inputs) / 32)
+
+        outputs = CuArray{ComplexF64}(undef, size(cu_inputs))
+
+        @cuda threads = ts blocks = bs always_inline = true kernel!(cu_inputs, outputs, length(cu_inputs))
+        CUDA.device_synchronize()
+        cu_results += abs2.(outputs)
+
+        println("Done.")
+        i += 1
+    end
+
+    println("Writing results")
+
+    out_ph_moms = getindex.(getindex.(input_momenta, 2), 1)
+    out_el_moms = getindex.(getindex.(input_momenta, 2), 2)
+
+    results = NamedDimsArray{(:omegas, :thetas, :phis)}(Array(cu_results))
+    println("Named results array: $(typeof(results))")
+
+    @save "$(photons)_congruent_photons_grid.jld2" omegas thetas phis results
+end
--- a/images/contour_plot_congruent_in_photons.gif
+++ b/images/contour_plot_congruent_in_photons.gif
--- a/notebooks/abc_model_large.ipynb
+++ b/notebooks/abc_model_large.ipynb
@ -413,7 +413,7 @@
 ],
 "metadata": {
  "kernelspec": {
-   "display_name": "Julia 1.10.2",
+   "display_name": "Julia 1.10.4",
   "language": "julia",
   "name": "julia-1.10"
  },
@ -421,7 +421,7 @@
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
-   "version": "1.10.2"
+   "version": "1.10.4"
  }
 },
 "nbformat": 4,
--- a/notebooks/congruent_ph_eval.ipynb
+++ b/notebooks/congruent_ph_eval.ipynb
--- a/results/1_congruent_photons_grid.jld2
+++ b/results/1_congruent_photons_grid.jld2
--- a/results/2_congruent_photons_grid.jld2
+++ b/results/2_congruent_photons_grid.jld2
--- a/results/3_congruent_photons_grid.jld2
+++ b/results/3_congruent_photons_grid.jld2
--- a/results/4_congruent_photons_grid.jld2
+++ b/results/4_congruent_photons_grid.jld2
--- a/results/5_congruent_photons_grid.jld2
+++ b/results/5_congruent_photons_grid.jld2
--- a/src/MetagraphOptimization.jl
+++ b/src/MetagraphOptimization.jl
@ -6,6 +6,8 @@ A module containing tools to work on DAGs.
 module MetagraphOptimization

 using QEDbase
+using QEDcore
+using QEDprocesses

 # graph types
 export DAG
@ -64,8 +66,7 @@ export ComputeTaskABC_Sum

 # QED model
 export FeynmanDiagram, FeynmanVertex, FeynmanTie, FeynmanParticle
-export PhotonStateful, FermionStateful, AntiFermionStateful
-export QEDParticle, QEDProcessDescription, QEDProcessInput, QEDModel
+export GenericQEDProcess, QEDModel
 export ComputeTaskQED_P
 export ComputeTaskQED_S1
 export ComputeTaskQED_S2
@ -99,9 +100,6 @@ export ==, in, show, isempty, delete!, length

 export bytes_to_human_readable

-# TODO: this is probably not good
-import QEDprocesses.compute
-
 import Base.length
 import Base.show
 import Base.==
@ -114,6 +112,7 @@ import Base.delete!
 import Base.insert!
 import Base.collect

+include("QEDprocesses_patch.jl")

 include("devices/interface.jl")
 include("task/type.jl")
@ -173,22 +172,26 @@ include("optimization/split.jl")
 include("models/interface.jl")
 include("models/print.jl")

-include("models/abc/types.jl")
-include("models/abc/particle.jl")
-include("models/abc/compute.jl")
-include("models/abc/create.jl")
-include("models/abc/properties.jl")
-include("models/abc/parse.jl")
-include("models/abc/print.jl")
+include("models/physics_models/interface.jl")

-include("models/qed/types.jl")
-include("models/qed/particle.jl")
-include("models/qed/diagrams.jl")
-include("models/qed/compute.jl")
-include("models/qed/create.jl")
-include("models/qed/properties.jl")
-include("models/qed/parse.jl")
-include("models/qed/print.jl")
+#=
+include("models/physics_models/abc/types.jl")
+include("models/physics_models/abc/particle.jl")
+include("models/physics_models/abc/compute.jl")
+include("models/physics_models/abc/create.jl")
+include("models/physics_models/abc/properties.jl")
+include("models/physics_models/abc/parse.jl")
+include("models/physics_models/abc/print.jl")
+=#
+
+include("models/physics_models/qed/types.jl")
+include("models/physics_models/qed/particle.jl")
+include("models/physics_models/qed/diagrams.jl")
+include("models/physics_models/qed/compute.jl")
+include("models/physics_models/qed/create.jl")
+include("models/physics_models/qed/properties.jl")
+include("models/physics_models/qed/parse.jl")
+include("models/physics_models/qed/print.jl")

 include("devices/measure.jl")
 include("devices/detect.jl")
@ -197,7 +200,7 @@ include("devices/impl.jl")
 include("devices/numa/impl.jl")
 include("devices/cuda/impl.jl")
 include("devices/rocm/impl.jl")
-include("devices/oneapi/impl.jl")
+#include("devices/oneapi/impl.jl")

 include("scheduler/interface.jl")
 include("scheduler/greedy.jl")
--- a/src/QEDprocesses_patch.jl
+++ b/src/QEDprocesses_patch.jl
@ -0,0 +1,25 @@
+# patch QEDprocesses
+# see issue https://github.com/QEDjl-project/QEDprocesses.jl/issues/77
+@inline function QEDprocesses.number_particles(
+    proc_def::QEDbase.AbstractProcessDefinition,
+    ::Type{PS},
+) where {
+    DIR <: QEDbase.ParticleDirection,
+    PT <: QEDbase.AbstractParticleType,
+    EL <: AbstractFourMomentum,
+    PS <: ParticleStateful{DIR, PT, EL},
+}
+    return QEDprocesses.number_particles(proc_def, DIR(), PT())
+end
+
+@inline function QEDcore.ParticleStateful{DIR, SPECIES}(
+    mom::AbstractFourMomentum,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType}
+    return ParticleStateful(DIR(), SPECIES(), mom)
+end
+
+@inline function QEDcore.ParticleStateful{DIR, SPECIES, EL}(
+    mom::EL,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType, EL <: AbstractFourMomentum}
+    return ParticleStateful(DIR(), SPECIES(), mom)
+end
--- a/src/code_gen/function.jl
+++ b/src/code_gen/function.jl
@ -1,19 +1,19 @@
 """
-    get_compute_function(graph::DAG, process::AbstractProcessDescription, machine::Machine)
+    get_compute_function(graph::DAG, instance, machine::Machine)

-Return a function of signature `compute_<id>(input::AbstractProcessInput)`, which will return the result of the DAG computation on the given input.
+Return a function of signature `compute_<id>(input::input_type(instance))`, which will return the result of the DAG computation on the given input.
 """
-function get_compute_function(graph::DAG, process::AbstractProcessDescription, machine::Machine)
-    tape = gen_tape(graph, process, machine)
+function get_compute_function(graph::DAG, instance, machine::Machine)
+    tape = gen_tape(graph, instance, machine)

    initCaches = Expr(:block, tape.initCachesCode...)
-    assignInputs = Expr(:block, expr_from_fc.(tape.inputAssignCode)...)
+    assignInputs = Expr(:block, tape.inputAssignCode...)
    code = Expr(:block, expr_from_fc.(tape.computeCode)...)

    functionId = to_var_name(UUIDs.uuid1(rng[1]))
    resSym = eval(gen_access_expr(entry_device(tape.machine), tape.outputSymbol))
    expr = Meta.parse(
-        "function compute_$(functionId)(data_input::AbstractProcessInput) $(initCaches); $(assignInputs); $code; return $resSym; end",
+        "function compute_$(functionId)(data_input::$(input_type(instance))) $(initCaches); $(assignInputs); $code; return $resSym; end",
    )

    func = eval(expr)
@ -22,15 +22,15 @@ function get_compute_function(graph::DAG, process::AbstractProcessDescription, m
 end

 """
-    get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machine::Machine)
+    get_cuda_kernel(graph::DAG, instance, machine::Machine)

 Return a function of signature `compute_<id>(input::CuVector, output::CuVector, n::Int64)`, which will return the result of the DAG computation of the input on the given output variable.
 """
-function get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machine::Machine)
-    tape = gen_tape(graph, process, machine)
+function get_cuda_kernel(graph::DAG, instance, machine::Machine)
+    tape = gen_tape(graph, instance, machine)

    initCaches = Expr(:block, tape.initCachesCode...)
-    assignInputs = Expr(:block, expr_from_fc.(tape.inputAssignCode)...)
+    assignInputs = Expr(:block, tape.inputAssignCode...)
    code = Expr(:block, expr_from_fc.(tape.computeCode)...)

    functionId = to_var_name(UUIDs.uuid1(rng[1]))
@ -54,19 +54,19 @@ function get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machin
 end

 """
-    execute(graph::DAG, process::AbstractProcessDescription, machine::Machine, input::AbstractProcessInput)
+    execute(graph::DAG, instance, machine::Machine, input)

-Execute the code of the given `graph` on the given input particles.
+Execute the code of the given `graph` on the given input values.

 This is essentially shorthand for
 ```julia
-tape = gen_tape(graph, process, machine)
+tape = gen_tape(graph, instance, machine)
 return execute_tape(tape, input)
 ```

 See also: [`parse_dag`](@ref), [`parse_process`](@ref), [`gen_process_input`](@ref)
 """
-function execute(graph::DAG, process::AbstractProcessDescription, machine::Machine, input::AbstractProcessInput)
-    tape = gen_tape(graph, process, machine)
+function execute(graph::DAG, instance, machine::Machine, input)
+    tape = gen_tape(graph, instance, machine)
    return execute_tape(tape, input)
 end
--- a/src/code_gen/tape_machine.jl
+++ b/src/code_gen/tape_machine.jl
@ -1,10 +1,11 @@
+# TODO: do this with macros
 function call_fc(fc::FunctionCall{VectorT, 0}, cache::Dict{Symbol, Any}) where {VectorT <: SVector{1}}
    cache[fc.return_symbol] = fc.func(cache[fc.arguments[1]])
    return nothing
 end

 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT <: SVector{1}}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], cache[fc.arguments[1]])
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], cache[fc.arguments[1]])
    return nothing
 end

@ -14,12 +15,12 @@ function call_fc(fc::FunctionCall{VectorT, 0}, cache::Dict{Symbol, Any}) where {
 end

 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT <: SVector{2}}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], cache[fc.arguments[1]], cache[fc.arguments[2]])
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], cache[fc.arguments[1]], cache[fc.arguments[2]])
    return nothing
 end

 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], getindex.(Ref(cache), fc.arguments)...)
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], getindex.(Ref(cache), fc.arguments)...)
    return nothing
 end

@ -31,7 +32,7 @@ Execute the given [`FunctionCall`](@ref) on the dictionary.
 Several more specialized versions of this function exist to reduce vector unrolling work for common cases.
 """
 function call_fc(fc::FunctionCall{VectorT, M}, cache::Dict{Symbol, Any}) where {VectorT, M}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments..., getindex.(Ref(cache), fc.arguments)...)
+    cache[fc.return_symbol] = fc.func(fc.value_arguments..., getindex.(Ref(cache), fc.arguments)...)
    return nothing
 end

@ -47,12 +48,8 @@ end
 For a given function call, return an expression evaluating it.
 """
 function expr_from_fc(fc::FunctionCall{VectorT, M}) where {VectorT, M}
-    func_call = Expr(
-        :call,
-        Symbol(fc.func),
-        fc.additional_arguments...,
-        eval.(gen_access_expr.(Ref(fc.device), fc.arguments))...,
-    )
+    func_call =
+        Expr(:call, Symbol(fc.func), fc.value_arguments..., eval.(gen_access_expr.(Ref(fc.device), fc.arguments))...)

    expr = :($(eval(gen_access_expr(fc.device, fc.return_symbol))) = $func_call)
    return expr
@ -73,51 +70,32 @@ function gen_cache_init_code(machine::Machine)
    return initializeCaches
 end

-"""
-    part_from_x(type::Type, index::Int, x::AbstractProcessInput)
-
-Return the [`ParticleValue`](@ref) of the given type of particle with the given `index` from the given process input.
-
-Function is wrapped into a [`FunctionCall`](@ref) in [`gen_input_assignment_code`](@ref).
-"""
-part_from_x(type::Type, index::Int, x::AbstractProcessInput) =
-    ParticleValue{type, ComplexF64}(get_particle(x, type, index), one(ComplexF64))
-
 """
    gen_input_assignment_code(
        inputSymbols::Dict{String, Vector{Symbol}},
-        processDescription::AbstractProcessDescription,
+        instance::AbstractProblemInstance,
        machine::Machine,
-        processInputSymbol::Symbol = :input,
+        problemInputSymbol::Symbol = :data_input,
    )

-Return a `Vector{Expr}` doing the input assignments from the given `processInputSymbol` onto the `inputSymbols`.
+Return a `Vector{Expr}` doing the input assignments from the given `problemInputSymbol` onto the `inputSymbols`.
 """
 function gen_input_assignment_code(
    inputSymbols::Dict{String, Vector{Symbol}},
-    processDescription::AbstractProcessDescription,
+    instance,
    machine::Machine,
-    processInputSymbol::Symbol = :input,
+    problemInputSymbol::Symbol = :data_input,
 )
-    @assert length(inputSymbols) >=
-            sum(values(in_particles(processDescription))) + sum(values(out_particles(processDescription))) "Number of input Symbols is smaller than the number of particles in the process description"
-
-    assignInputs = Vector{FunctionCall}()
+    assignInputs = Vector{Expr}()
    for (name, symbols) in inputSymbols
-        (type, index) = type_index_from_name(model(processDescription), name)
        # make a function for this, since we can't use anonymous functions in the FunctionCall

        for symbol in symbols
            device = entry_device(machine)
            push!(
                assignInputs,
-                FunctionCall(
-                    # x is the process input
-                    part_from_x,
-                    SVector{1, Symbol}(processInputSymbol),
-                    SVector{2, Any}(type, index),
-                    symbol,
-                    device,
+                Meta.parse(
+                    "$(eval(gen_access_expr(device, symbol))) = $(input_expr(instance, name, problemInputSymbol))",
                ),
            )
        end
@ -127,14 +105,14 @@ function gen_input_assignment_code(
 end

 """
-    gen_tape(graph::DAG, process::AbstractProcessDescription, machine::Machine)
+    gen_tape(graph::DAG, instance::AbstractProblemInstance, machine::Machine, scheduler::AbstractScheduler = GreedyScheduler())

 Generate the code for a given graph. The return value is a [`Tape`](@ref).

 See also: [`execute`](@ref), [`execute_tape`](@ref)
 """
-function gen_tape(graph::DAG, process::AbstractProcessDescription, machine::Machine)
-    schedule = schedule_dag(GreedyScheduler(), graph, machine)
+function gen_tape(graph::DAG, instance, machine::Machine, scheduler::AbstractScheduler = GreedyScheduler())
+    schedule = schedule_dag(scheduler, graph, machine)

    # get inSymbols
    inputSyms = Dict{String, Vector{Symbol}}()
@ -150,23 +128,24 @@ function gen_tape(graph::DAG, process::AbstractProcessDescription, machine::Mach
    outSym = Symbol(to_var_name(get_exit_node(graph).id))

    initCaches = gen_cache_init_code(machine)
-    assignInputs = gen_input_assignment_code(inputSyms, process, machine, :input)
+    assignInputs = gen_input_assignment_code(inputSyms, instance, machine, :data_input)

-    return Tape(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), process, machine)
+    return Tape{input_type(instance)}(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine)
 end

 """
-    execute_tape(tape::Tape, input::AbstractProcessInput)
+    execute_tape(tape::Tape, input::Input) where {Input}

 Execute the given tape with the given input.

 For implementation reasons, this disregards the set [`CacheStrategy`](@ref) of the devices and always uses a dictionary.
 """
-function execute_tape(tape::Tape, input::AbstractProcessInput)
+function execute_tape(tape::Tape, input)
    cache = Dict{Symbol, Any}()
-    cache[:input] = input
+    cache[:data_input] = input
    # simply execute all the code snippets here
-    # TODO: `@assert` that process input fits the tape.process
+    @assert typeof(input) == input_type(tape.instance)
+    # TODO: `@assert` that input fits the tape.instance
    for expr in tape.initCachesCode
        @eval $expr
    end
--- a/src/code_gen/type.jl
+++ b/src/code_gen/type.jl
@ -1,19 +1,21 @@

 """
-    Tape
+    Tape{INPUT}

 TODO: update docs
+- `INPUT` the input type of the problem instance
+
 - `code::Vector{Expr}`: The julia expression containing the code for the whole graph.
 - `inputSymbols::Dict{String, Vector{Symbol}}`: A dictionary of symbols mapping the names of the input nodes of the graph to the symbols their inputs should be provided on.
 - `outputSymbol::Symbol`: The symbol of the final calculated value
 """
-struct Tape
+struct Tape{INPUT}
    initCachesCode::Vector{Expr}
-    inputAssignCode::Vector{FunctionCall}
+    inputAssignCode::Vector{Expr}
    computeCode::Vector{FunctionCall}
    inputSymbols::Dict{String, Vector{Symbol}}
    outputSymbol::Symbol
    cache::Dict{Symbol, Any}
-    process::AbstractProcessDescription
+    instance::Any
    machine::Machine
 end
--- a/src/models/interface.jl
+++ b/src/models/interface.jl
@ -1,120 +1,46 @@
-import QEDbase.mass
-import QEDbase.AbstractParticle

 """
-    AbstractPhysicsModel
+    AbstractModel

-Base type for a model, e.g. ABC-Model or QED. This is used to dispatch many functions.
+Base type for all models. From this, [`AbstractProblemInstance`](@ref)s can be constructed.
+
+See also: [`problem_instance`](@ref)
 """
-abstract type AbstractPhysicsModel end
+abstract type AbstractModel end

 """
-    ParticleValue{ParticleType <: AbstractParticle}
+    problem_instance(::AbstractModel, ::Vararg)

-A struct describing a particle during a calculation of a Feynman Diagram, together with the value that's being calculated. `AbstractParticle` is the type from the QEDbase package.
-
-`sizeof(ParticleValue())` = 48 Byte
+Interface function that must be implemented for any implementation of [`AbstractModel`](@ref). This function should return a specific [`AbstractProblemInstance`](@ref) given some parameters.
 """
-struct ParticleValue{ParticleType <: AbstractParticle, ValueType}
-    p::ParticleType
-    v::ValueType
-end
+function problem_instance end

 """
-    AbstractProcessDescription
+    AbstractProblemInstance

-Base type for process descriptions. An object of this type of a corresponding [`AbstractPhysicsModel`](@ref) should uniquely identify a process in that model.
+Base type for problem instances. An object of this type of a corresponding [`AbstractModel`](@ref) should uniquely identify a problem instance of that model.

 See also: [`parse_process`](@ref)
 """
-abstract type AbstractProcessDescription end
+abstract type AbstractProblemInstance end

 """
-    AbstractProcessInput
+    input_type(problem::AbstractProblemInstance)

-Base type for process inputs. An object of this type contains the input values (e.g. momenta) of the particles in a process.
-
-See also: [`gen_process_input`](@ref)
+Return the fully specified input type for a specific [`AbstractProblemInstance`](@ref).
 """
-abstract type AbstractProcessInput end
+function input_type end

 """
-    interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: AbstractParticle, T2 <: AbstractParticle}
+    graph(::AbstractProblemInstance)

-Interface function that must be implemented for every subtype of [`AbstractParticle`](@ref), returning the result particle type when the two given particles interact.
+Generate the [`DAG`](@ref) for the given [`AbstractProblemInstance`](@ref). Every entry node (see [`get_entry_nodes`](@ref)) to the graph must have a name set. Implement [`input_expr`](@ref) to return a valid expression for each of those names.
 """
-function interaction_result end
+function graph end

 """
-    types(::AbstractPhysicsModel)
+    input_expr(instance::AbstractProblemInstance, name::String, input_symbol::Symbol)

-Interface function that must be implemented for every subtype of [`AbstractPhysicsModel`](@ref), returning a `Vector` of the available particle types in the model.
+For the given [`AbstractProblemInstance`](@ref), the entry node name, and the symbol of the problem input (where a variable of type `input_type(...)` will exist), return an `Expr` that gets that specific input value from the input symbol.
 """
-function types end
-
-"""
-    in_particles(::AbstractProcessDescription)
-
-Interface function that must be implemented for every subtype of [`AbstractProcessDescription`](@ref).
-Returns a `<: Dict{Type{AbstractParticle}, Int}` object, representing the number of incoming particles for the process per particle type.
-
-
-    in_particles(::AbstractProcessInput)
-
-Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
-Returns a `<: Vector{AbstractParticle}` object with the values of all incoming particles for the corresponding `ProcessDescription`.
-"""
-function in_particles end
-
-"""
-    out_particles(::AbstractProcessDescription)
-
-Interface function that must be implemented for every subtype of [`AbstractProcessDescription`](@ref).
-Returns a `<: Dict{Type{AbstractParticle}, Int}` object, representing the number of outgoing particles for the process per particle type.
-
-
-    out_particles(::AbstractProcessInput)
-
-Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
-Returns a `<: Vector{AbstractParticle}` object with the values of all outgoing particles for the corresponding `ProcessDescription`.
-"""
-function out_particles end
-
-"""
-    get_particle(::AbstractProcessInput, t::Type, n::Int)
-
-Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
-Returns the `n`th particle of type `t`.
-"""
-function get_particle end
-
-"""
-    parse_process(::AbstractString, ::AbstractPhysicsModel)
-
-Interface function that must be implemented for every subtype of [`AbstractPhysicsModel`](@ref).
-Returns a `ProcessDescription` object.
-"""
-function parse_process end
-
-"""
-    gen_process_input(::AbstractProcessDescription)
-
-Interface function that must be implemented for every specific [`AbstractProcessDescription`](@ref).
-Returns a randomly generated and valid corresponding `ProcessInput`.
-"""
-function gen_process_input end
-
-"""
-    model(::AbstractProcessDescription)
-    model(::AbstarctProcessInput)
-
-Return the model of this process description or input.
-"""
-function model end
-
-"""
-    type_from_name(model::AbstractModel, name::String)
-
-For a name of a particle in the given [`AbstractModel`](@ref), return the particle's [`Type`] and index as a tuple. The input string can be expetced to be of the form \"<name><index>\".
-"""
-function type_index_from_name end
+function input_expr end
--- a/src/models/physics_models/abc/compute.jl
+++ b/src/models/physics_models/abc/compute.jl
--- a/src/models/physics_models/abc/create.jl
+++ b/src/models/physics_models/abc/create.jl
--- a/src/models/physics_models/abc/parse.jl
+++ b/src/models/physics_models/abc/parse.jl
--- a/src/models/physics_models/abc/particle.jl
+++ b/src/models/physics_models/abc/particle.jl
--- a/src/models/physics_models/abc/print.jl
+++ b/src/models/physics_models/abc/print.jl
--- a/src/models/physics_models/abc/properties.jl
+++ b/src/models/physics_models/abc/properties.jl
--- a/src/models/physics_models/abc/types.jl
+++ b/src/models/physics_models/abc/types.jl
--- a/src/models/physics_models/interface.jl
+++ b/src/models/physics_models/interface.jl
@ -0,0 +1,142 @@
+import QEDbase.AbstractParticle
+
+"""
+    AbstractPhysicsModel
+
+Base type for a model, e.g. ABC-Model or QED. This is used to dispatch many functions.
+"""
+abstract type AbstractPhysicsModel <: AbstractModel end
+
+"""
+    ParticleValue{ParticleType <: AbstractParticleStateful}
+
+A struct describing a particle during a calculation of a Feynman Diagram, together with the value that's being calculated. `AbstractParticleStateful` is the type from the QEDbase package.
+
+`sizeof(ParticleValue())` = 48 Byte
+"""
+struct ParticleValue{ParticleType <: AbstractParticleStateful, ValueType}
+    p::ParticleType
+    v::ValueType
+end
+
+"""
+TBW
+
+particle value + spin/pol info, only used on the external legs (u tasks)
+"""
+struct ParticleValueSP{ParticleType <: AbstractParticleStateful, SP <: AbstractSpinOrPolarization, ValueType}
+    p::ParticleType
+    v::ValueType
+    sp::SP
+end
+
+"""
+    AbstractProcessDescription <: AbstractProblemInstance
+
+Base type for particle scattering process descriptions. An object of this type of a corresponding [`AbstractPhysicsModel`](@ref) should uniquely identify a scattering process in that model.
+
+See also: [`parse_process`](@ref), [`AbstractProblemInstance`](@ref)
+"""
+abstract type AbstractProcessDescription end
+
+
+#TODO: i don't think giving this a base type is a good idea, the input type should just be returned of some function, allowing anything as an input type
+"""
+    AbstractProcessInput
+
+Base type for process inputs. An object of this type contains the input values (e.g. momenta) of the particles in a process.
+
+See also: [`gen_process_input`](@ref)
+"""
+abstract type AbstractProcessInput end
+
+"""
+    interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: AbstractParticle, T2 <: AbstractParticle}
+
+Interface function that must be implemented for every subtype of [`AbstractParticle`](@ref), returning the result particle type when the two given particles interact.
+"""
+function interaction_result end
+
+"""
+    types(::AbstractPhysicsModel)
+
+Interface function that must be implemented for every subtype of [`AbstractPhysicsModel`](@ref), returning a `Vector` of the available particle types in the model.
+"""
+function types end
+
+"""
+    in_particles(::AbstractProcessDescription)
+
+Interface function that must be implemented for every subtype of [`AbstractProcessDescription`](@ref).
+Returns a `<: Dict{Type{AbstractParticle}, Int}` object, representing the number of incoming particles for the process per particle type.
+
+
+    in_particles(::AbstractProcessInput)
+
+Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
+Returns a `<: Vector{AbstractParticle}` object with the values of all incoming particles for the corresponding `ProcessDescription`.
+"""
+function in_particles end
+
+"""
+    out_particles(::AbstractProcessDescription)
+
+Interface function that must be implemented for every subtype of [`AbstractProcessDescription`](@ref).
+Returns a `<: Dict{Type{AbstractParticle}, Int}` object, representing the number of outgoing particles for the process per particle type.
+
+
+    out_particles(::AbstractProcessInput)
+
+Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
+Returns a `<: Vector{AbstractParticle}` object with the values of all outgoing particles for the corresponding `ProcessDescription`.
+"""
+function out_particles end
+
+"""
+    get_particle(::AbstractProcessInput, t::Type, n::Int)
+
+Interface function that must be implemented for every subtype of [`AbstractProcessInput`](@ref).
+Returns the `n`th particle of type `t`.
+"""
+function get_particle end
+
+"""
+    parse_process(::AbstractString, ::AbstractPhysicsModel)
+
+Interface function that must be implemented for every subtype of [`AbstractPhysicsModel`](@ref).
+Returns a `ProcessDescription` object.
+"""
+function parse_process end
+
+"""
+    gen_process_input(::AbstractProcessDescription)
+
+Interface function that must be implemented for every specific [`AbstractProcessDescription`](@ref).
+Returns a randomly generated and valid corresponding `ProcessInput`.
+"""
+function gen_process_input end
+
+"""
+    model(::AbstractProcessDescription)
+    model(::AbstractProcessInput)
+
+Return the model of this process description or input.
+"""
+function model end
+
+"""
+    type_from_name(model::AbstractModel, name::String)
+
+For a name of a particle in the given [`AbstractModel`](@ref), return the particle's [`Type`] and index as a tuple. The input string can be expetced to be of the form \"<name><index>\".
+"""
+function type_index_from_name end
+
+"""
+    part_from_x(type::Type, index::Int, x::AbstractProcessInput)
+
+Return the [`ParticleValue`](@ref) of the given type of particle with the given `index` from the given process input.
+
+Function is wrapped into a [`FunctionCall`](@ref) in [`gen_input_assignment_code`](@ref).
+"""
+part_from_x(type::Type, index::Int, x::AbstractProcessInput) =
+    ParticleValue{type, ComplexF64}(get_particle(x, type, index), one(ComplexF64))
--- a/src/models/physics_models/qed/compute.jl
+++ b/src/models/physics_models/qed/compute.jl
@ -1,23 +1,40 @@
 using StaticArrays

-"""
-    compute(::ComputeTaskQED_P, data::QEDParticleValue)
+construction_string(::Incoming) = "Incoming()"
+construction_string(::Outgoing) = "Outgoing()"

-Return the particle as is and initialize the Value.
-"""
-function compute(::ComputeTaskQED_P, data::QEDParticleValue{P}) where {P <: QEDParticle}
-    # TODO do we actually need this for anything?
-    return ParticleValue{P, DiracMatrix}(data.p, one(DiracMatrix))
+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 input_expr(instance::GenericQEDProcess, name::String, psp_symbol::Symbol)
+    (type, index) = type_index_from_name(QEDModel(), name)
+
+    return Meta.parse(
+        "ParticleValueSP(
+    $type(momentum($psp_symbol, $(construction_string(particle_direction(type))), $(construction_string(particle_species(type))), Val($index))),
+    0.0im,
+    $(construction_string(spin_or_pol(instance, type, index))),
+)",
+    )
 end

 """
-    compute(::ComputeTaskQED_U, data::QEDParticleValue)
+    compute(::ComputeTaskQED_U, data::ParticleValueSP)

 Compute an outer edge. Return the particle value with the same particle and the value multiplied by an outer_edge factor.
 """
-function compute(::ComputeTaskQED_U, data::PV) where {P <: QEDParticle, PV <: QEDParticleValue{P}}
+function compute(
+    ::ComputeTaskQED_U,
+    data::ParticleValueSP{P, SP, V},
+) where {P <: ParticleStateful, V <: ValueType, SP <: AbstractSpinOrPolarization}
    part::P = data.p
-    state = base_state(particle(part), direction(part), momentum(part), spin_or_pol(part))
+    state = base_state(particle_species(part), particle_direction(part), momentum(part), SP())
    return ParticleValue{P, typeof(state)}(
        data.p,
        state, # will return a SLorentzVector{ComplexF64}, BiSpinor or AdjointBiSpinor
@ -25,15 +42,15 @@ function compute(::ComputeTaskQED_U, data::PV) where {P <: QEDParticle, PV <: QE
 end

 """
-    compute(::ComputeTaskQED_V, data1::QEDParticleValue, data2::QEDParticleValue)
+    compute(::ComputeTaskQED_V, data1::ParticleValue, data2::ParticleValue)

 Compute a vertex. Preserve momentum and particle types (e + gamma->p etc.) to create resulting particle, multiply values together and times a vertex factor.
 """
 function compute(
    ::ComputeTaskQED_V,
-    data1::PV1,
-    data2::PV2,
-) where {P1 <: QEDParticle, P2 <: QEDParticle, PV1 <: QEDParticleValue{P1}, PV2 <: QEDParticleValue{P2}}
+    data1::ParticleValue{P1, V1},
+    data2::ParticleValue{P2, V2},
+) where {P1 <: ParticleStateful, P2 <: ParticleStateful, V1 <: ValueType, V2 <: ValueType}
    p3 = QED_conserve_momentum(data1.p, data2.p)
    P3 = interaction_result(P1, P2)
    state = QED_vertex()
@ -53,7 +70,7 @@ function compute(
 end

 """
-    compute(::ComputeTaskQED_S2, data1::QEDParticleValue, data2::QEDParticleValue)
+    compute(::ComputeTaskQED_S2, data1::ParticleValue, data2::ParticleValue)

 Compute a final inner edge (2 input particles, no output particle).

@ -63,9 +80,19 @@ For valid inputs, both input particles should have the same momenta at this poin
 """
 function compute(
    ::ComputeTaskQED_S2,
-    data1::ParticleValue{P1},
-    data2::ParticleValue{P2},
-) where {P1 <: Union{AntiFermionStateful, FermionStateful}, P2 <: Union{AntiFermionStateful, FermionStateful}}
+    data1::ParticleValue{P1, V1},
+    data2::ParticleValue{P2, V2},
+) where {
+    D1 <: ParticleDirection,
+    D2 <: ParticleDirection,
+    S1 <: Union{Electron, Positron},
+    S2 <: Union{Electron, Positron},
+    V1 <: ValueType,
+    V2 <: ValueType,
+    EL <: AbstractFourMomentum,
+    P1 <: ParticleStateful{D1, S1, EL},
+    P2 <: ParticleStateful{D2, S2, EL},
+}
    #@assert isapprox(data1.p.momentum, data2.p.momentum, rtol = sqrt(eps()), atol = sqrt(eps())) "$(data1.p.momentum) vs. $(data2.p.momentum)"

    inner = QED_inner_edge(propagation_result(P1)(momentum(data1.p)))
@ -80,9 +107,9 @@ end

 function compute(
    ::ComputeTaskQED_S2,
-    data1::ParticleValue{P1},
-    data2::ParticleValue{P2},
-) where {P1 <: PhotonStateful, P2 <: PhotonStateful}
+    data1::ParticleValue{ParticleStateful{D1, Photon}, V1},
+    data2::ParticleValue{ParticleStateful{D2, Photon}, V2},
+) where {D1 <: ParticleDirection, D2 <: ParticleDirection, V1 <: ValueType, V2 <: ValueType}
    # TODO: assert that data1 and data2 are opposites
    inner = QED_inner_edge(data1.p)
    # inner edge is just a scalar, data1 and data2 are photon states that are just Complex numbers here
@ -90,11 +117,11 @@ function compute(
 end

 """
-    compute(::ComputeTaskQED_S1, data::QEDParticleValue)
+    compute(::ComputeTaskQED_S1, data::ParticleValue)

 Compute inner edge (1 input particle, 1 output particle).
 """
-function compute(::ComputeTaskQED_S1, data::QEDParticleValue{P}) where {P <: QEDParticle}
+function compute(::ComputeTaskQED_S1, data::ParticleValue{P, V}) where {P <: ParticleStateful, V <: ValueType}
    newP = propagation_result(P)
    new_p = newP(momentum(data.p))
    # inner edge is just a scalar, can multiply from either side
--- a/src/models/physics_models/qed/create.jl
+++ b/src/models/physics_models/qed/create.jl
@ -1,83 +1,58 @@

 ComputeTaskQED_Sum() = ComputeTaskQED_Sum(0)

-function _svector_from_type(processDescription::QEDProcessDescription, type, particles)
+function _svector_from_type(processDescription::GenericQEDProcess, type, particles)
    if haskey(in_particles(processDescription), type)
        return SVector{in_particles(processDescription)[type], type}(filter(x -> typeof(x) <: type, particles))
    end
    if haskey(out_particles(processDescription), type)
        return SVector{out_particles(processDescription)[type], type}(filter(x -> typeof(x) <: type, particles))
    end
-    return SVector{0, type}()
 end

 """
-    gen_process_input(processDescription::QEDProcessDescription)
+    gen_process_input(processDescription::GenericQEDProcess)

-Return a ProcessInput of randomly generated [`QEDParticle`](@ref)s from a [`QEDProcessDescription`](@ref). The process description can be created manually or parsed from a string using [`parse_process`](@ref).
+Return a ProcessInput of randomly generated [`QEDParticle`](@ref)s from a [`GenericQEDProcess`](@ref). The process description can be created manually or parsed from a string using [`parse_process`](@ref).

 Note: This uses RAMBO to create a valid process with conservation of momentum and energy.
 """
-function gen_process_input(processDescription::QEDProcessDescription)
+function gen_process_input(processDescription::GenericQEDProcess)
    massSum = 0
    inputMasses = Vector{Float64}()
-    for (particle, n) in processDescription.inParticles
-        for _ in 1:n
-            massSum += mass(particle)
-            push!(inputMasses, mass(particle))
-        end
+    for particle in incoming_particles(processDescription)
+        massSum += mass(particle)
+        push!(inputMasses, mass(particle))
    end
    outputMasses = Vector{Float64}()
-    for (particle, n) in processDescription.outParticles
-        for _ in 1:n
-            massSum += mass(particle)
-            push!(outputMasses, mass(particle))
-        end
+    for particle in outgoing_particles(processDescription)
+        massSum += mass(particle)
+        push!(outputMasses, mass(particle))
    end

    # add some extra random mass to allow for some momentum
    massSum += rand(rng[threadid()]) * (length(inputMasses) + length(outputMasses))

-
-    particles = Vector{QEDParticle}()
-    initialMomenta = generate_initial_moms(massSum, inputMasses)
-    index = 1
-    for (particle, n) in processDescription.inParticles
-        for _ in 1:n
-            mom = initialMomenta[index]
-            push!(particles, particle(mom))
-            index += 1
-        end
-    end
-
+    initial_momenta = generate_initial_moms(massSum, inputMasses)
    final_momenta = generate_physical_massive_moms(rng[threadid()], massSum, outputMasses)
-    index = 1
-    for (particle, n) in processDescription.outParticles
-        for _ in 1:n
-            push!(particles, particle(final_momenta[index]))
-            index += 1
-        end
-    end

-    inFerms = _svector_from_type(processDescription, FermionStateful{Incoming, SpinUp}, particles)
-    outFerms = _svector_from_type(processDescription, FermionStateful{Outgoing, SpinUp}, particles)
-    inAntiferms = _svector_from_type(processDescription, AntiFermionStateful{Incoming, SpinUp}, particles)
-    outAntiferms = _svector_from_type(processDescription, AntiFermionStateful{Outgoing, SpinUp}, particles)
-    inPhotons = _svector_from_type(processDescription, PhotonStateful{Incoming, PolX}, particles)
-    outPhotons = _svector_from_type(processDescription, PhotonStateful{Outgoing, PolX}, particles)
-
-    processInput =
-        QEDProcessInput(processDescription, inFerms, outFerms, inAntiferms, outAntiferms, inPhotons, outPhotons)
+    processInput = PhaseSpacePoint(
+        processDescription,
+        PerturbativeQED(),
+        PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
+        tuple(initial_momenta...),
+        tuple(final_momenta...),
+    )

    return processInput
 end

 """
-    gen_graph(process_description::QEDProcessDescription)
+    gen_graph(process_description::GenericQEDProcess)

-For a given [`QEDProcessDescription`](@ref), return the [`DAG`](@ref) that computes it.
+For a given [`GenericQEDProcess`](@ref), return the [`DAG`](@ref) that computes it.
 """
-function gen_graph(process_description::QEDProcessDescription)
+function gen_graph(process_description::GenericQEDProcess)
    initial_diagram = FeynmanDiagram(process_description)
    diagrams = gen_diagrams(initial_diagram)

@ -88,9 +63,9 @@ function gen_graph(process_description::QEDProcessDescription)

    # TODO: Not all diagram outputs should always be summed at the end, if they differ by fermion exchange they need to be diffed
    # Should not matter for n-Photon Compton processes though
-    sum_node = insert_node!(graph, make_node(ComputeTaskQED_Sum(0)), track = false, invalidate_cache = false)
-    global_data_out = insert_node!(graph, make_node(DataTask(COMPLEX_SIZE)), track = false, invalidate_cache = false)
-    insert_edge!(graph, sum_node, global_data_out, track = false, invalidate_cache = false)
+    sum_node = insert_node!(graph, make_node(ComputeTaskQED_Sum(0)); track = false, invalidate_cache = false)
+    global_data_out = insert_node!(graph, make_node(DataTask(COMPLEX_SIZE)); track = false, invalidate_cache = false)
+    insert_edge!(graph, sum_node, global_data_out; track = false, invalidate_cache = false)

    # remember the data out nodes for connection
    dataOutNodes = Dict()
@ -99,16 +74,16 @@ function gen_graph(process_description::QEDProcessDescription)
        # generate data in and U tasks
        data_in = insert_node!(
            graph,
-            make_node(DataTask(PARTICLE_VALUE_SIZE), String(particle)),
+            make_node(DataTask(PARTICLE_VALUE_SIZE), String(particle));
            track = false,
            invalidate_cache = false,
        ) # read particle data node
-        compute_u = insert_node!(graph, make_node(ComputeTaskQED_U()), track = false, invalidate_cache = false) # compute U node
+        compute_u = insert_node!(graph, make_node(ComputeTaskQED_U()); track = false, invalidate_cache = false) # compute U node
        data_out =
-            insert_node!(graph, make_node(DataTask(PARTICLE_VALUE_SIZE)), track = false, invalidate_cache = false) # transfer data out from u (one ABCParticleValue object)
+            insert_node!(graph, make_node(DataTask(PARTICLE_VALUE_SIZE)); track = false, invalidate_cache = false) # transfer data out from u (one ABCParticleValue object)

-        insert_edge!(graph, data_in, compute_u, track = false, invalidate_cache = false)
-        insert_edge!(graph, compute_u, data_out, track = false, invalidate_cache = false)
+        insert_edge!(graph, data_in, compute_u; track = false, invalidate_cache = false)
+        insert_edge!(graph, compute_u, data_out; track = false, invalidate_cache = false)

        # remember the data_out node for future edges
        dataOutNodes[String(particle)] = data_out
@ -124,19 +99,19 @@ function gen_graph(process_description::QEDProcessDescription)
                data_in1 = dataOutNodes[String(vertex.in1)]
                data_in2 = dataOutNodes[String(vertex.in2)]

-                compute_V = insert_node!(graph, make_node(ComputeTaskQED_V()), track = false, invalidate_cache = false) # compute vertex
+                compute_V = insert_node!(graph, make_node(ComputeTaskQED_V()); track = false, invalidate_cache = false) # compute vertex

-                insert_edge!(graph, data_in1, compute_V, track = false, invalidate_cache = false)
-                insert_edge!(graph, data_in2, compute_V, track = false, invalidate_cache = false)
+                insert_edge!(graph, data_in1, compute_V; track = false, invalidate_cache = false)
+                insert_edge!(graph, data_in2, compute_V; track = false, invalidate_cache = false)

                data_V_out = insert_node!(
                    graph,
-                    make_node(DataTask(PARTICLE_VALUE_SIZE)),
+                    make_node(DataTask(PARTICLE_VALUE_SIZE));
                    track = false,
                    invalidate_cache = false,
                )

-                insert_edge!(graph, compute_V, data_V_out, track = false, invalidate_cache = false)
+                insert_edge!(graph, compute_V, data_V_out; track = false, invalidate_cache = false)

                if (vertex.out == tie.in1 || vertex.out == tie.in2)
                    # out particle is part of the tie -> there will be an S2 task with it later, don't make S1 task
@ -146,18 +121,18 @@ function gen_graph(process_description::QEDProcessDescription)

                # otherwise, add S1 task
                compute_S1 =
-                    insert_node!(graph, make_node(ComputeTaskQED_S1()), track = false, invalidate_cache = false) # compute propagator
+                    insert_node!(graph, make_node(ComputeTaskQED_S1()); track = false, invalidate_cache = false) # compute propagator

-                insert_edge!(graph, data_V_out, compute_S1, track = false, invalidate_cache = false)
+                insert_edge!(graph, data_V_out, compute_S1; track = false, invalidate_cache = false)

                data_S1_out = insert_node!(
                    graph,
-                    make_node(DataTask(PARTICLE_VALUE_SIZE)),
+                    make_node(DataTask(PARTICLE_VALUE_SIZE));
                    track = false,
                    invalidate_cache = false,
                )

-                insert_edge!(graph, compute_S1, data_S1_out, track = false, invalidate_cache = false)
+                insert_edge!(graph, compute_S1, data_S1_out; track = false, invalidate_cache = false)

                # overrides potentially different nodes from previous diagrams, which is intentional
                dataOutNodes[String(vertex.out)] = data_S1_out
@ -168,16 +143,16 @@ function gen_graph(process_description::QEDProcessDescription)
        data_in1 = dataOutNodes[String(tie.in1)]
        data_in2 = dataOutNodes[String(tie.in2)]

-        compute_S2 = insert_node!(graph, make_node(ComputeTaskQED_S2()), track = false, invalidate_cache = false)
+        compute_S2 = insert_node!(graph, make_node(ComputeTaskQED_S2()); track = false, invalidate_cache = false)

-        data_S2 = insert_node!(graph, make_node(DataTask(PARTICLE_VALUE_SIZE)), track = false, invalidate_cache = false)
+        data_S2 = insert_node!(graph, make_node(DataTask(PARTICLE_VALUE_SIZE)); track = false, invalidate_cache = false)

-        insert_edge!(graph, data_in1, compute_S2, track = false, invalidate_cache = false)
-        insert_edge!(graph, data_in2, compute_S2, track = false, invalidate_cache = false)
+        insert_edge!(graph, data_in1, compute_S2; track = false, invalidate_cache = false)
+        insert_edge!(graph, data_in2, compute_S2; track = false, invalidate_cache = false)

-        insert_edge!(graph, compute_S2, data_S2, track = false, invalidate_cache = false)
+        insert_edge!(graph, compute_S2, data_S2; track = false, invalidate_cache = false)

-        insert_edge!(graph, data_S2, sum_node, track = false, invalidate_cache = false)
+        insert_edge!(graph, data_S2, sum_node; track = false, invalidate_cache = false)
        add_child!(task(sum_node))
    end

--- a/src/models/physics_models/qed/diagrams.jl
+++ b/src/models/physics_models/qed/diagrams.jl
@ -8,10 +8,10 @@ import Base.show
 """
    FeynmanParticle

-Representation of a particle for use in [`FeynmanDiagram`](@ref)s. Consist of the [`QEDParticle`](@ref) type and an id.
+Representation of a particle for use in [`FeynmanDiagram`](@ref)s. Consist of the `ParticleStateful` type and an id.
 """
 struct FeynmanParticle
-    particle::Type{<:QEDParticle}
+    particle::Type{<:ParticleStateful}
    id::Int
 end

@ -51,31 +51,21 @@ struct FeynmanDiagram
 end

 """
-    FeynmanDiagram(pd::QEDProcessDescription)
+    FeynmanDiagram(pd::GenericQEDProcess)

 Create an initial [`FeynmanDiagram`](@ref) with only its initial particles set and no vertices or ties.

 Use [`gen_diagrams`](@ref) to generate all possible diagrams from this one.
 """
-function FeynmanDiagram(pd::QEDProcessDescription)
+function FeynmanDiagram(pd::GenericQEDProcess)
    parts = Vector{FeynmanParticle}()
-    for (type, n) in pd.inParticles
-        for i in 1:n
-            push!(parts, FeynmanParticle(type, i))
-        end
-    end
-    for (type, n) in pd.outParticles
-        for i in 1:n
-            push!(parts, FeynmanParticle(type, i))
-        end
-    end
+
    ids = Dict{Type, Int64}()
-    for t in types(QEDModel())
-        if (isincoming(t))
-            ids[t] = get(pd.inParticles, t, 0)
-        else
-            ids[t] = get(pd.outParticles, t, 0)
+    for type in types(model(pd))
+        for i in 1:number_particles(pd, type)
+            push!(parts, FeynmanParticle(type, i))
        end
+        ids[type] = number_particles(pd, type)
    end

    return FeynmanDiagram([], missing, parts, ids)
@ -83,7 +73,7 @@ end

 function particle_after_tie(p::FeynmanParticle, t::FeynmanTie)
    if p == t.in1 || p == t.in2
-        return FeynmanParticle(FermionStateful{Incoming, SpinUp}, -1) # placeholder particle and id for tied particles
+        return FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, -1) # placeholder particle and id for tied particles
    end
    return p
 end
@ -114,7 +104,7 @@ end
 Return a string representation of the [`FeynmanParticle`](@ref) in a format that is readable by [`type_index_from_name`](@ref).
 """
 function String(p::FeynmanParticle)
-    return "$(String(p.particle))$(String(direction(p.particle)))$(p.id)"
+    return "$(String(p.particle))$(String(particle_direction(p.particle)))$(p.id)"
 end

 function hash(v::FeynmanVertex)
@ -162,15 +152,16 @@ function ==(d1::FeynmanDiagram, d2::FeynmanDiagram)
    )=#
 end

-copy(fd::FeynmanDiagram) =
-    FeynmanDiagram(deepcopy(fd.vertices), copy(fd.tie[]), deepcopy(fd.particles), copy(fd.type_ids))
+function copy(fd::FeynmanDiagram)
+    return FeynmanDiagram(deepcopy(fd.vertices), copy(fd.tie[]), deepcopy(fd.particles), copy(fd.type_ids))
+end

 """
-    id_for_type(d::FeynmanDiagram, t::Type{<:QEDParticle})
+    id_for_type(d::FeynmanDiagram, t::Type{<:ParticleStateful})

 Return the highest id of any particle of the given type in the diagram + 1.
 """
-function id_for_type(d::FeynmanDiagram, t::Type{<:QEDParticle})
+function id_for_type(d::FeynmanDiagram, t::Type{<:ParticleStateful})
    return d.type_ids[t] + 1
 end

@ -439,18 +430,19 @@ function remove_duplicates(compare_set::Set{FeynmanDiagram})
    return result
 end

+
 """
    is_compton(fd::FeynmanDiagram)

 Returns true iff the given feynman diagram is an (empty) diagram of a compton process like ke->k^ne
 """
 function is_compton(fd::FeynmanDiagram)
-    return fd.type_ids[FermionStateful{Incoming, SpinUp}] == 1 &&
-           fd.type_ids[FermionStateful{Outgoing, SpinUp}] == 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
+    return fd.type_ids[ParticleStateful{Incoming, Electron, SFourMomentum}] == 1 &&
+           fd.type_ids[ParticleStateful{Outgoing, Electron, SFourMomentum}] == 1 &&
+           fd.type_ids[ParticleStateful{Incoming, Positron, SFourMomentum}] == 0 &&
+           fd.type_ids[ParticleStateful{Outgoing, Positron, SFourMomentum}] == 0 &&
+           fd.type_ids[ParticleStateful{Incoming, Photon, SFourMomentum}] >= 1 &&
+           fd.type_ids[ParticleStateful{Outgoing, Photon, SFourMomentum}] >= 1
 end

 """
@ -460,8 +452,8 @@ Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram f
 """
 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],
+        [FeynmanParticle(ParticleStateful{Incoming, Photon, SFourMomentum}, i) for i in 1:n],
+        [FeynmanParticle(ParticleStateful{Outgoing, Photon, SFourMomentum}, i) for i in 1:m],
    )

    new_diagram = FeynmanDiagram(
@ -469,10 +461,10 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
        missing,
        [inFerm, outFerm, photons...],
        Dict{Type, Int64}(
-            FermionStateful{Incoming, SpinUp} => 1,
-            FermionStateful{Outgoing, SpinUp} => 1,
-            PhotonStateful{Incoming, PolX} => n,
-            PhotonStateful{Outgoing, PolX} => m,
+            ParticleStateful{Incoming, Electron, SFourMomentum} => 1,
+            ParticleStateful{Outgoing, Electron, SFourMomentum} => 1,
+            ParticleStateful{Incoming, Photon, SFourMomentum} => n,
+            ParticleStateful{Outgoing, Photon, SFourMomentum} => m,
        ),
    )

@ -484,9 +476,9 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
    while left_index <= right_index
        # left side
        v_left = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations),
            photons[order[left_index]],
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations + 1),
        )
        left_index += 1
        add_vertex!(new_diagram, v_left)
@ -497,9 +489,9 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::

        # right side
        v_right = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations),
            photons[order[right_index]],
-            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations + 1),
        )
        right_index -= 1
        add_vertex!(new_diagram, v_right)
@ -512,7 +504,6 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
    return new_diagram
 end

-
 """
    gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)

@ -520,8 +511,8 @@ Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram f
 """
 function gen_compton_diagram_from_order_one_side(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],
+        [FeynmanParticle(ParticleStateful{Incoming, Photon, SFourMomentum}, i) for i in 1:n],
+        [FeynmanParticle(ParticleStateful{Outgoing, Photon, SFourMomentum}, i) for i in 1:m],
    )

    new_diagram = FeynmanDiagram(
@ -529,10 +520,10 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
        missing,
        [inFerm, outFerm, photons...],
        Dict{Type, Int64}(
-            FermionStateful{Incoming, SpinUp} => 1,
-            FermionStateful{Outgoing, SpinUp} => 1,
-            PhotonStateful{Incoming, PolX} => n,
-            PhotonStateful{Outgoing, PolX} => m,
+            ParticleStateful{Incoming, Electron, SFourMomentum} => 1,
+            ParticleStateful{Outgoing, Electron, SFourMomentum} => 1,
+            ParticleStateful{Incoming, Photon, SFourMomentum} => n,
+            ParticleStateful{Outgoing, Photon, SFourMomentum} => m,
        ),
    )

@ -544,9 +535,9 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
    while left_index <= right_index
        # left side
        v_left = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations),
            photons[order[left_index]],
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations + 1),
        )
        left_index += 1
        add_vertex!(new_diagram, v_left)
@ -559,9 +550,9 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
        if (iterations == 1)
            # right side
            v_right = FeynmanVertex(
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
+                FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations),
                photons[order[right_index]],
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
+                FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations + 1),
            )
            right_index -= 1
            add_vertex!(new_diagram, v_right)
@ -576,15 +567,14 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
    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)
+    inFerm = FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, 1)
+    outFerm = FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, 1)

    perms = [permutations([i for i in 1:(n + m)])...]

@ -596,15 +586,14 @@ function gen_compton_diagrams(n::Int, m::Int)
    return vcat(diagrams...)
 end

-
 """
    gen_compton_diagrams_one_side(n::Int, m::Int)

 Special case diagram generation for Compton processes, i.e., processes of the form k^ne->k^me, but generating from one end, yielding larger diagrams
 """
 function gen_compton_diagrams_one_side(n::Int, m::Int)
-    inFerm = FeynmanParticle(FermionStateful{Incoming, SpinUp}, 1)
-    outFerm = FeynmanParticle(FermionStateful{Outgoing, SpinUp}, 1)
+    inFerm = FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, 1)
+    outFerm = FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, 1)

    perms = [permutations([i for i in 1:(n + m)])...]

@ -623,12 +612,15 @@ 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 gen_compton_diagrams_one_side(
-            fd.type_ids[PhotonStateful{Incoming, PolX}],
-            fd.type_ids[PhotonStateful{Outgoing, PolX}],
+        return gen_compton_diagrams(
+            fd.type_ids[ParticleStateful{Incoming, Photon, SFourMomentum}],
+            fd.type_ids[ParticleStateful{Outgoing, Photon, SFourMomentum}],
        )
    end

+    throw(error("Unimplemented for non-compton!"))
+
+    #=
    working = Set{FeynmanDiagram}()
    results = Set{FeynmanDiagram}()

@ -667,4 +659,5 @@ function gen_diagrams(fd::FeynmanDiagram)
    end

    return remove_duplicates(results)
+    =#
 end
--- a/src/models/physics_models/qed/parse.jl
+++ b/src/models/physics_models/qed/parse.jl
@ -0,0 +1,69 @@
+
+# TODO generalize this somehow, probably using AllSpin/AllPol as defaults
+"""
+    parse_process(string::AbstractString, model::QEDModel)
+
+Parse a string representation of a process, such as "ke->ke" into the corresponding [`QEDProcessDescription`](@ref).
+"""
+function parse_process(
+    str::AbstractString,
+    model::QEDModel,
+    inphpol::AbstractDefinitePolarization,
+    inelspin::AbstractDefiniteSpin,
+    outphpol::AbstractDefinitePolarization,
+    outelspin::AbstractDefiniteSpin,
+)
+    inParticles = Dict(
+        ParticleStateful{Incoming, Photon, SFourMomentum} => 0,
+        ParticleStateful{Incoming, Electron, SFourMomentum} => 0,
+        ParticleStateful{Incoming, Positron, SFourMomentum} => 0,
+    )
+    outParticles = Dict(
+        ParticleStateful{Outgoing, Photon, SFourMomentum} => 0,
+        ParticleStateful{Outgoing, Electron, SFourMomentum} => 0,
+        ParticleStateful{Outgoing, Positron, SFourMomentum} => 0,
+    )
+
+    if !(contains(str, "->"))
+        throw("Did not find -> while parsing process \"$str\"")
+    end
+
+    (inStr, outStr) = split(str, "->")
+
+    if (isempty(inStr) || isempty(outStr))
+        throw("Process (\"$str\") input or output part is empty!")
+    end
+
+    for t in types(model)
+        if (is_incoming(t))
+            inCount = count(x -> x == String(t)[1], inStr)
+
+            if inCount != 0
+                inParticles[t] = inCount
+            end
+        end
+        if (is_outgoing(t))
+            outCount = count(x -> x == String(t)[1], outStr)
+            if outCount != 0
+                outParticles[t] = outCount
+            end
+        end
+    end
+
+    if length(inStr) != sum(values(inParticles))
+        throw("Encountered unknown characters in the input part of process \"$str\"")
+    elseif length(outStr) != sum(values(outParticles))
+        throw("Encountered unknown characters in the output part of process \"$str\"")
+    end
+
+    in_ph = inParticles[ParticleStateful{Incoming, Photon, SFourMomentum}]
+    in_el = inParticles[ParticleStateful{Incoming, Electron, SFourMomentum}]
+    in_pos = inParticles[ParticleStateful{Incoming, Positron, SFourMomentum}]
+    out_ph = outParticles[ParticleStateful{Outgoing, Photon, SFourMomentum}]
+    out_el = outParticles[ParticleStateful{Outgoing, Electron, SFourMomentum}]
+    out_pos = outParticles[ParticleStateful{Outgoing, Positron, SFourMomentum}]
+    in_spin_pols = tuple([i <= in_ph ? inphpol : inelspin for i in 1:(in_ph + in_el)]...)
+    out_spin_pols = tuple([i <= out_ph ? outphpol : outelspin for i in 1:(out_ph + out_el)]...)
+
+    return GenericQEDProcess(in_ph, out_ph, in_el, out_el, in_pos, out_pos, in_spin_pols, out_spin_pols)
+end
--- a/src/models/physics_models/qed/particle.jl
+++ b/src/models/physics_models/qed/particle.jl
@ -0,0 +1,377 @@
+using QEDprocesses
+using StaticArrays
+
+const e = sqrt(4π / 137)
+
+"""
+    QEDModel <: AbstractPhysicsModel
+
+Singleton definition for identification of the QED-Model.
+"""
+struct QEDModel <: AbstractPhysicsModel end
+
+QEDbase.is_incoming(::Type{<:ParticleStateful{Incoming}}) = true
+QEDbase.is_outgoing(::Type{<:ParticleStateful{Outgoing}}) = true
+QEDbase.is_incoming(::Type{<:ParticleStateful{Outgoing}}) = false
+QEDbase.is_outgoing(::Type{<:ParticleStateful{Incoming}}) = false
+
+QEDbase.particle_direction(::Type{<:ParticleStateful{DIR}}) where {DIR <: ParticleDirection} = DIR()
+QEDbase.particle_species(
+    ::Type{<:ParticleStateful{DIR, SPECIES}},
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType} = SPECIES()
+
+_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"))
+
+struct GenericQEDProcess{INT, OUTT, INSP, OUTSP} <:
+       AbstractProcessDefinition where {INT <: Tuple, OUTT <: Tuple, INSP <: Tuple, OUTSP <: Tuple}
+    incoming_particles::INT
+    outgoing_particles::OUTT
+
+    incoming_spins_pols::INSP
+    outgoing_spins_pols::OUTSP
+
+    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)
+
+        # TODO: check in_sp/out_sp fits to particles (spin->fermion, pol->photon)
+
+        return new{INT, OUTT, INSP, OUTSP}(in_particles, out_particles, in_sp, out_sp)
+    end
+
+    """
+        GenericQEDProcess{INSP, OUTSP}(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int, in_sp::INSP, out_sp::OUTSP)
+
+    Convenience constructor from numbers of input/output photons, electrons and positrons.
+
+    `in_sp` and `out_sp` are tuples of the spin and polarization configurations of the given particles. Their order is Photons, then Electrons, then Positrons for both incoming and outgoing particles.
+
+    # TODO: make default for in_sp and out_sp available for convenience
+    """
+    function GenericQEDProcess(
+        in_ph::Int,
+        out_ph::Int,
+        in_el::Int,
+        out_el::Int,
+        in_po::Int,
+        out_po::Int,
+        in_sp::INSP,
+        out_sp::OUTSP,
+    ) where {INSP <: Tuple, OUTSP <: Tuple}
+        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, in_sp, out_sp)
+    end
+end
+
+function spin_or_pol(
+    process::GenericQEDProcess,
+    type::Type{ParticleStateful{DIR, SPECIES, EL}},
+    n::Int,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType, EL <: AbstractFourMomentum}
+    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 DIR <: Incoming
+        return process.incoming_spins_pols[i]
+    elseif DIR <: Outgoing
+        return process.outgoing_spins_pols[i]
+    else
+        throw(InvalidInputError("unknown direction $(DIR()) given"))
+    end
+end
+
+QEDprocesses.incoming_particles(proc::GenericQEDProcess) = proc.incoming_particles
+QEDprocesses.outgoing_particles(proc::GenericQEDProcess) = proc.outgoing_particles
+
+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
+end
+
+function input_type(p::GenericQEDProcess)
+    in_t = QEDcore._assemble_tuple_type(incoming_particles(p), Incoming(), SFourMomentum)
+    out_t = QEDcore._assemble_tuple_type(outgoing_particles(p), Outgoing(), SFourMomentum)
+    return PhaseSpacePoint{
+        typeof(p),
+        PerturbativeQED,
+        PhasespaceDefinition{SphericalCoordinateSystem, ElectronRestFrame},
+        Tuple{in_t...},
+        Tuple{out_t...},
+    }
+end
+
+ValueType = Union{BiSpinor, AdjointBiSpinor, DiracMatrix, SLorentzVector{Float64}, ComplexF64}
+
+"""
+    interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
+
+For two given particle types that can interact, return the third.
+"""
+function interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: ParticleStateful, T2 <: ParticleStateful}
+    @assert false "Invalid interaction between particles of types $t1 and $t2"
+end
+
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{Incoming, Positron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{DIR, Photon, EL}},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = ParticleStateful{Outgoing, Electron, EL}
+
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{Outgoing, Positron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{DIR, Photon, EL}},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = ParticleStateful{Incoming, Electron, EL}
+
+# antifermion mirror
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Positron, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum} = interaction_result(ParticleStateful{Outgoing, Electron, EL}, t2)
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Positron, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum} = interaction_result(ParticleStateful{Incoming, Electron, EL}, t2)
+
+# photon commutativity
+interaction_result(
+    t1::Type{ParticleStateful{DIR, Photon, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = interaction_result(t2, t1)
+
+# but prevent stack overflow
+function interaction_result(
+    t1::Type{ParticleStateful{DIR1, Photon, EL}},
+    t2::Type{ParticleStateful{DIR2, Photon, EL}},
+) where {DIR1 <: ParticleDirection, DIR2 <: ParticleDirection, EL <: AbstractFourMomentum}
+    @assert false "Invalid interaction between particles of types $t1 and $t2"
+end
+
+"""
+    propagation_result(t1::Type{T}) where {T <: QEDParticle}
+
+Return the type of the inverted direction. E.g.
+"""
+propagation_result(::Type{ParticleStateful{Incoming, Electron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Electron, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Electron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Electron, EL}
+propagation_result(::Type{ParticleStateful{Incoming, Positron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Positron, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Positron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Positron, EL}
+propagation_result(::Type{ParticleStateful{Incoming, Photon, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Photon, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Photon, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Photon, EL}
+
+"""
+    types(::QEDModel)
+
+Return a Vector of the possible types of particle in the [`QEDModel`](@ref).
+"""
+function types(::QEDModel)
+    return [
+        ParticleStateful{Incoming, Photon, SFourMomentum},
+        ParticleStateful{Outgoing, Photon, SFourMomentum},
+        ParticleStateful{Incoming, Electron, SFourMomentum},
+        ParticleStateful{Outgoing, Electron, SFourMomentum},
+        ParticleStateful{Incoming, Positron, SFourMomentum},
+        ParticleStateful{Outgoing, Positron, SFourMomentum},
+    ]
+end
+
+# type piracy?
+String(::Type{Incoming}) = "Incoming"
+String(::Type{Outgoing}) = "Outgoing"
+
+String(::Type{PolX}) = "polx"
+String(::Type{PolY}) = "poly"
+
+String(::Type{SpinUp}) = "spinup"
+String(::Type{SpinDown}) = "spindown"
+
+String(::Incoming) = "i"
+String(::Outgoing) = "o"
+
+function String(::Type{<:ParticleStateful{DIR, Photon}}) where {DIR <: ParticleDirection}
+    return "k"
+end
+function String(::Type{<:ParticleStateful{DIR, Electron}}) where {DIR <: ParticleDirection}
+    return "e"
+end
+function String(::Type{<:ParticleStateful{DIR, Positron}}) where {DIR <: ParticleDirection}
+    return "p"
+end
+
+"""
+    caninteract(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
+
+For two given `ParticleStateful` types, return whether they can interact at a vertex. This is equivalent to `!issame(T1, T2)`.
+
+See also: [`issame`](@ref) and [`interaction_result`](@ref)
+"""
+function caninteract(
+    T1::Type{<:ParticleStateful{D1, S1}},
+    T2::Type{<:ParticleStateful{D2, S2}},
+) where {D1 <: ParticleDirection, S1 <: AbstractParticleType, D2 <: ParticleDirection, S2 <: AbstractParticleType}
+    if (T1 == T2)
+        return false
+    end
+    if (S1 == Photon && S2 == Photon)
+        return false
+    end
+
+    for (P1, P2) in [(T1, T2), (T2, T1)]
+        if (P1 <: ParticleStateful{Incoming, Electron} && P2 <: ParticleStateful{Outgoing, Positron})
+            return false
+        end
+        if (P1 <: ParticleStateful{Outgoing, Electron} && P2 <: ParticleStateful{Incoming, Positron})
+            return false
+        end
+    end
+
+    return true
+end
+
+function type_index_from_name(::QEDModel, name::String)
+    if startswith(name, "ki")
+        return (ParticleStateful{Incoming, Photon, SFourMomentum}, parse(Int, name[3:end]))
+    elseif startswith(name, "ko")
+        return (ParticleStateful{Outgoing, Photon, SFourMomentum}, parse(Int, name[3:end]))
+    elseif startswith(name, "ei")
+        return (ParticleStateful{Incoming, Electron, SFourMomentum}, parse(Int, name[3:end]))
+    elseif startswith(name, "eo")
+        return (ParticleStateful{Outgoing, Electron, SFourMomentum}, parse(Int, name[3:end]))
+    elseif startswith(name, "pi")
+        return (ParticleStateful{Incoming, Positron, SFourMomentum}, parse(Int, name[3:end]))
+    elseif startswith(name, "po")
+        return (ParticleStateful{Outgoing, Positron, SFourMomentum}, parse(Int, name[3:end]))
+    else
+        throw("Invalid name for a particle in the QED model")
+    end
+end
+
+"""
+    issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
+
+For two given `ParticleStateful` types, return whether they are equivalent for the purpose of a Feynman Diagram. That means e.g. an `Incoming` `AntiFermion` is the same as an `Outgoing` `Fermion`. This is equivalent to `!caninteract(T1, T2)`.
+
+See also: [`caninteract`](@ref) and [`interaction_result`](@ref)
+"""
+function issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
+    return !caninteract(T1, T2)
+end
+
+"""
+    QED_vertex()
+
+Return the factor of a vertex in a QED feynman diagram.
+"""
+@inline function QED_vertex()::SLorentzVector{DiracMatrix}
+    # Peskin-Schroeder notation
+    return -1im * e * gamma()
+end
+
+@inline function QED_inner_edge(p::ParticleStateful)
+    return propagator(particle_species(p), momentum(p))
+end
+
+"""
+    QED_conserve_momentum(p1::ParticleStateful, p2::ParticleStateful)
+
+Calculate and return a new particle from two given interacting ones at a vertex.
+"""
+function QED_conserve_momentum(
+    p1::P1,
+    p2::P2,
+) where {
+    Dir1 <: ParticleDirection,
+    Dir2 <: ParticleDirection,
+    Species1 <: AbstractParticleType,
+    Species2 <: AbstractParticleType,
+    P1 <: ParticleStateful{Dir1, Species1},
+    P2 <: ParticleStateful{Dir2, Species2},
+}
+    P3 = interaction_result(P1, P2)
+    p1_mom = momentum(p1)
+    if (Dir1 <: Outgoing)
+        p1_mom *= -1
+    end
+    p2_mom = momentum(p2)
+    if (Dir2 <: Outgoing)
+        p2_mom *= -1
+    end
+
+    p3_mom = p1_mom + p2_mom
+    if (particle_direction(P3) isa Incoming)
+        return ParticleStateful(particle_direction(P3), particle_species(P3), -p3_mom)
+    end
+    return ParticleStateful(particle_direction(P3), particle_species(P3), p3_mom)
+end
+
+"""
+    model(::AbstractProcessDescription)
+
+Return the model of this process description.
+"""
+model(::GenericQEDProcess) = QEDModel()
+model(::PhaseSpacePoint) = QEDModel()
+
+function get_particle(
+    input::PhaseSpacePoint,
+    t::Type{ParticleStateful{DIR, SPECIES}},
+    n::Int,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType}
+    i = 0
+    for p in particles(input, DIR())
+        if p isa t
+            i += 1
+            if i == n
+                return p
+            end
+        end
+    end
+    @assert false "Invalid type given"
+end
--- a/src/models/physics_models/qed/print.jl
+++ b/src/models/physics_models/qed/print.jl
@ -1,8 +1,9 @@

+#=
 """
-    show(io::IO, process::QEDProcessDescription)
+    show(io::IO, process::GenericQEDProcess)

-Pretty print an [`QEDProcessDescription`](@ref) (no newlines).
+Pretty print an [`GenericQEDProcess`](@ref) (no newlines).

 ```jldoctest
 julia> using MetagraphOptimization
@ -14,7 +15,7 @@ julia> print(parse_process("kk->ep", QEDModel()))
 QED Process: 'kk->ep'
 ```
 """
-function show(io::IO, process::QEDProcessDescription)
+function show(io::IO, process::GenericQEDProcess)
    # types() gives the types in order (QED) instead of random like keys() would
    print(io, "QED Process: \'")
    for type in types(QEDModel())
@ -34,7 +35,7 @@ end


 """
-    String(process::QEDProcessDescription)
+    String(process::GenericQEDProcess)

 Create a short string suitable as a filename or similar, describing the given process.

@ -64,7 +65,9 @@ function String(process::QEDProcessDescription)
    end
    return str
 end
+=#

+#=
 """
    show(io::IO, processInput::QEDProcessInput)

@ -92,7 +95,9 @@ function show(io::IO, processInput::QEDProcessInput)
    end
    return nothing
 end
+=#

+#=
 """
    show(io::IO, particle::T) where {T <: QEDParticle}

@ -102,13 +107,14 @@ function show(io::IO, particle::T) where {T <: QEDParticle}
    print(io, "$(String(typeof(particle))): $(particle.momentum)")
    return nothing
 end
+=#

 """
    show(io::IO, particle::FeynmanParticle)

 Pretty print a [`FeynmanParticle`](@ref) (no newlines).
 """
-show(io::IO, p::FeynmanParticle) = print(io, "$(String(p.particle))_$(String(direction(p.particle)))_$(p.id)")
+show(io::IO, p::FeynmanParticle) = print(io, "$(String(p.particle))_$(String(particle_direction(p.particle)))_$(p.id)")

 """
    show(io::IO, particle::FeynmanVertex)
--- a/src/models/physics_models/qed/properties.jl
+++ b/src/models/physics_models/qed/properties.jl
--- a/src/models/physics_models/qed/types.jl
+++ b/src/models/physics_models/qed/types.jl
--- a/src/models/print.jl
+++ b/src/models/print.jl
@ -1,10 +0,0 @@
-
-"""
-    show(io::IO, particleValue::ParticleValue)
-
-Pretty print a [`ParticleValue`](@ref), no newlines.
-"""
-function show(io::IO, particleValue::ParticleValue)
-    print(io, "($(particleValue.p), value: $(particleValue.v))")
-    return nothing
-end
--- a/src/models/qed/parse.jl
+++ b/src/models/qed/parse.jl
@ -1,44 +0,0 @@
-
-"""
-    parse_process(string::AbstractString, model::QEDModel)
-
-Parse a string representation of a process, such as "ke->ke" into the corresponding [`QEDProcessDescription`](@ref).
-"""
-function parse_process(str::AbstractString, model::QEDModel)
-    inParticles = Dict{Type, Int}()
-    outParticles = Dict{Type, Int}()
-
-    if !(contains(str, "->"))
-        throw("Did not find -> while parsing process \"$str\"")
-    end
-
-    (inStr, outStr) = split(str, "->")
-
-    if (isempty(inStr) || isempty(outStr))
-        throw("Process (\"$str\") input or output part is empty!")
-    end
-
-    for t in types(model)
-        if (isincoming(t))
-            inCount = count(x -> x == String(t)[1], inStr)
-
-            if inCount != 0
-                inParticles[t] = inCount
-            end
-        end
-        if (isoutgoing(t))
-            outCount = count(x -> x == String(t)[1], outStr)
-            if outCount != 0
-                outParticles[t] = outCount
-            end
-        end
-    end
-
-    if length(inStr) != sum(values(inParticles))
-        throw("Encountered unknown characters in the input part of process \"$str\"")
-    elseif length(outStr) != sum(values(outParticles))
-        throw("Encountered unknown characters in the output part of process \"$str\"")
-    end
-
-    return QEDProcessDescription(inParticles, outParticles)
-end
--- a/src/models/qed/particle.jl
+++ b/src/models/qed/particle.jl
@ -1,408 +0,0 @@
-using QEDprocesses
-using StaticArrays
-import QEDbase.mass
-
-# TODO check
-const e = sqrt(4π / 137)
-
-"""
-    QEDModel <: AbstractPhysicsModel
-
-Singleton definition for identification of the QED-Model.
-"""
-struct QEDModel <: AbstractPhysicsModel end
-
-"""
-    QEDParticle
-
-Base type for all particles in the [`QEDModel`](@ref).
-
-Its template parameter specifies the particle's direction.
-
-The concrete types contain singletons of the types that they are, like `Photon` and `Electron` from QEDbase, and their state descriptions.
-"""
-abstract type QEDParticle{Direction <: ParticleDirection} <: AbstractParticle end
-
-"""
-    QEDProcessDescription <: AbstractProcessDescription
-
-A description of a process in the QED-Model. Contains the input and output particles.
-
-See also: [`in_particles`](@ref), [`out_particles`](@ref), [`parse_process`](@ref)
-"""
-struct QEDProcessDescription <: AbstractProcessDescription
-    inParticles::Dict{Type{<:QEDParticle{Incoming}}, Int}
-    outParticles::Dict{Type{<:QEDParticle{Outgoing}}, Int}
-end
-
-QEDParticleValue{ParticleType <: QEDParticle} = Union{
-    ParticleValue{ParticleType, BiSpinor},
-    ParticleValue{ParticleType, AdjointBiSpinor},
-    ParticleValue{ParticleType, DiracMatrix},
-    ParticleValue{ParticleType, SLorentzVector{Float64}},
-    ParticleValue{ParticleType, ComplexF64},
-}
-
-"""
-    PhotonStateful <: QEDParticle
-
-A photon of the [`QEDModel`](@ref) with its state.
-"""
-struct PhotonStateful{Direction <: ParticleDirection, Pol <: AbstractDefinitePolarization} <: QEDParticle{Direction}
-    momentum::SFourMomentum
-end
-
-PhotonStateful{Direction}(mom::SFourMomentum) where {Direction <: ParticleDirection} =
-    PhotonStateful{Direction, PolX}(mom)
-
-PhotonStateful{Dir, Pol}(ph::PhotonStateful) where {Dir, Pol} = PhotonStateful{Dir, Pol}(ph.momentum)
-
-"""
-    FermionStateful <: QEDParticle
-
-A fermion of the [`QEDModel`](@ref) with its state.
-"""
-struct FermionStateful{Direction <: ParticleDirection, Spin <: AbstractDefiniteSpin} <: QEDParticle{Direction}
-    momentum::SFourMomentum
-    # TODO: mass for electron/muon/tauon representation?
-end
-
-FermionStateful{Direction}(mom::SFourMomentum) where {Direction <: ParticleDirection} =
-    FermionStateful{Direction, SpinUp}(mom)
-
-FermionStateful{Dir, Spin}(f::FermionStateful) where {Dir, Spin} = FermionStateful{Dir, Spin}(f.momentum)
-
-"""
-    AntiFermionStateful <: QEDParticle
-
-An anti-fermion of the [`QEDModel`](@ref) with its state.
-"""
-struct AntiFermionStateful{Direction <: ParticleDirection, Spin <: AbstractDefiniteSpin} <: QEDParticle{Direction}
-    momentum::SFourMomentum
-    # TODO: mass for electron/muon/tauon representation?
-end
-
-AntiFermionStateful{Direction}(mom::SFourMomentum) where {Direction <: ParticleDirection} =
-    AntiFermionStateful{Direction, SpinUp}(mom)
-
-AntiFermionStateful{Dir, Spin}(f::AntiFermionStateful) where {Dir, Spin} = AntiFermionStateful{Dir, Spin}(f.momentum)
-
-"""
-    interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
-
-For two given particle types that can interact, return the third.
-"""
-function interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
-    @assert false "Invalid interaction between particles of types $t1 and $t2"
-end
-
-interaction_result(
-    ::Type{FermionStateful{Incoming, Spin1}},
-    ::Type{FermionStateful{Outgoing, Spin2}},
-) where {Spin1, Spin2} = PhotonStateful{Incoming, PolX}
-interaction_result(
-    ::Type{FermionStateful{Incoming, Spin1}},
-    ::Type{AntiFermionStateful{Incoming, Spin2}},
-) where {Spin1, Spin2} = PhotonStateful{Incoming, PolX}
-interaction_result(::Type{FermionStateful{Incoming, Spin1}}, ::Type{<:PhotonStateful}) where {Spin1} =
-    FermionStateful{Outgoing, SpinUp}
-
-interaction_result(
-    ::Type{FermionStateful{Outgoing, Spin1}},
-    ::Type{FermionStateful{Incoming, Spin2}},
-) where {Spin1, Spin2} = PhotonStateful{Incoming, PolX}
-interaction_result(
-    ::Type{FermionStateful{Outgoing, Spin1}},
-    ::Type{AntiFermionStateful{Outgoing, Spin2}},
-) where {Spin1, Spin2} = PhotonStateful{Incoming, PolX}
-interaction_result(::Type{FermionStateful{Outgoing, Spin1}}, ::Type{<:PhotonStateful}) where {Spin1} =
-    FermionStateful{Incoming, SpinUp}
-
-# antifermion mirror
-interaction_result(::Type{AntiFermionStateful{Incoming, Spin}}, t2::Type{<:QEDParticle}) where {Spin} =
-    interaction_result(FermionStateful{Outgoing, Spin}, t2)
-interaction_result(::Type{AntiFermionStateful{Outgoing, Spin}}, t2::Type{<:QEDParticle}) where {Spin} =
-    interaction_result(FermionStateful{Incoming, Spin}, t2)
-
-# photon commutativity
-interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:QEDParticle}) = interaction_result(t2, t1)
-
-# but prevent stack overflow
-function interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:PhotonStateful})
-    @assert false "Invalid interaction between particles of types $t1 and $t2"
-end
-
-"""
-    propagation_result(t1::Type{T}) where {T <: QEDParticle}
-
-Return the type of the inverted direction. E.g.
-"""
-propagation_result(::Type{FermionStateful{Incoming, Spin}}) where {Spin <: AbstractDefiniteSpin} =
-    FermionStateful{Outgoing, Spin}
-propagation_result(::Type{FermionStateful{Outgoing, Spin}}) where {Spin <: AbstractDefiniteSpin} =
-    FermionStateful{Incoming, Spin}
-propagation_result(::Type{AntiFermionStateful{Incoming, Spin}}) where {Spin <: AbstractDefiniteSpin} =
-    AntiFermionStateful{Outgoing, Spin}
-propagation_result(::Type{AntiFermionStateful{Outgoing, Spin}}) where {Spin <: AbstractDefiniteSpin} =
-    AntiFermionStateful{Incoming, Spin}
-propagation_result(::Type{PhotonStateful{Incoming, Pol}}) where {Pol <: AbstractDefinitePolarization} =
-    PhotonStateful{Outgoing, Pol}
-propagation_result(::Type{PhotonStateful{Outgoing, Pol}}) where {Pol <: AbstractDefinitePolarization} =
-    PhotonStateful{Incoming, Pol}
-
-"""
-    types(::QEDModel)
-
-Return a Vector of the possible types of particle in the [`QEDModel`](@ref).
-"""
-function types(::QEDModel)
-    return [
-        PhotonStateful{Incoming, PolX},
-        PhotonStateful{Outgoing, PolX},
-        FermionStateful{Incoming, SpinUp},
-        FermionStateful{Outgoing, SpinUp},
-        AntiFermionStateful{Incoming, SpinUp},
-        AntiFermionStateful{Outgoing, SpinUp},
-    ]
-end
-
-# type piracy?
-String(::Type{Incoming}) = "Incoming"
-String(::Type{Outgoing}) = "Outgoing"
-
-String(::Type{PolX}) = "polx"
-String(::Type{PolY}) = "poly"
-
-String(::Type{SpinUp}) = "spinup"
-String(::Type{SpinDown}) = "spindown"
-
-String(::Incoming) = "i"
-String(::Outgoing) = "o"
-
-function String(::Type{<:PhotonStateful})
-    return "k"
-end
-function String(::Type{<:FermionStateful})
-    return "e"
-end
-function String(::Type{<:AntiFermionStateful})
-    return "p"
-end
-
-function unique_name(::Type{PhotonStateful{Dir, Pol}}) where {Dir, Pol}
-    return String(PhotonStateful) * String(Dir) * String(Pol)
-end
-function unique_name(::Type{FermionStateful{Dir, Spin}}) where {Dir, Spin}
-    return String(FermionStateful) * String(Dir) * String(Spin)
-end
-function unique_name(::Type{AntiFermionStateful{Dir, Spin}}) where {Dir, Spin}
-    return String(AntiFermionStateful) * String(Dir) * String(Spin)
-end
-
-@inline particle(::PhotonStateful) = Photon()
-@inline particle(::FermionStateful) = Electron()
-@inline particle(::AntiFermionStateful) = Positron()
-
-@inline momentum(p::PhotonStateful)::SFourMomentum = p.momentum
-@inline momentum(p::FermionStateful)::SFourMomentum = p.momentum
-@inline momentum(p::AntiFermionStateful)::SFourMomentum = p.momentum
-
-@inline spin_or_pol(p::PhotonStateful{Dir, Pol}) where {Dir, Pol <: AbstractDefinitePolarization} = Pol()
-@inline spin_or_pol(p::FermionStateful{Dir, Spin}) where {Dir, Spin <: AbstractDefiniteSpin} = Spin()
-@inline spin_or_pol(p::AntiFermionStateful{Dir, Spin}) where {Dir, Spin <: AbstractDefiniteSpin} = Spin()
-
-@inline direction(
-    ::Type{P},
-) where {P <: Union{FermionStateful{Incoming}, AntiFermionStateful{Incoming}, PhotonStateful{Incoming}}} = Incoming()
-@inline direction(
-    ::Type{P},
-) where {P <: Union{FermionStateful{Outgoing}, AntiFermionStateful{Outgoing}, PhotonStateful{Outgoing}}} = Outgoing()
-
-@inline direction(
-    ::P,
-) where {P <: Union{FermionStateful{Incoming}, AntiFermionStateful{Incoming}, PhotonStateful{Incoming}}} = Incoming()
-@inline direction(
-    ::P,
-) where {P <: Union{FermionStateful{Outgoing}, AntiFermionStateful{Outgoing}, PhotonStateful{Outgoing}}} = Outgoing()
-
-@inline isincoming(::QEDParticle{Incoming}) = true
-@inline isincoming(::QEDParticle{Outgoing}) = false
-@inline isoutgoing(::QEDParticle{Incoming}) = false
-@inline isoutgoing(::QEDParticle{Outgoing}) = true
-
-@inline isincoming(::Type{<:QEDParticle{Incoming}}) = true
-@inline isincoming(::Type{<:QEDParticle{Outgoing}}) = false
-@inline isoutgoing(::Type{<:QEDParticle{Incoming}}) = false
-@inline isoutgoing(::Type{<:QEDParticle{Outgoing}}) = true
-
-@inline mass(::Type{<:FermionStateful}) = 1.0
-@inline mass(::Type{<:AntiFermionStateful}) = 1.0
-@inline mass(::Type{<:PhotonStateful}) = 0.0
-
-@inline invert_momentum(p::FermionStateful{Dir, Spin}) where {Dir, Spin} =
-    FermionStateful{Dir, Spin}(-p.momentum, p.spin)
-@inline invert_momentum(p::AntiFermionStateful{Dir, Spin}) where {Dir, Spin} =
-    AntiFermionStateful{Dir, Spin}(-p.momentum, p.spin)
-@inline invert_momentum(k::PhotonStateful{Dir, Spin}) where {Dir, Spin} =
-    PhotonStateful{Dir, Spin}(-k.momentum, k.polarization)
-
-
-"""
-    caninteract(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
-
-For two given [`QEDParticle`](@ref) types, return whether they can interact at a vertex. This is equivalent to `!issame(T1, T2)`.
-
-See also: [`issame`](@ref) and [`interaction_result`](@ref)
-"""
-function caninteract(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
-    if (T1 == T2)
-        return false
-    end
-    if (T1 <: PhotonStateful && T2 <: PhotonStateful)
-        return false
-    end
-
-    for (P1, P2) in [(T1, T2), (T2, T1)]
-        if (P1 <: FermionStateful{Incoming} && P2 <: AntiFermionStateful{Outgoing})
-            return false
-        end
-        if (P1 <: FermionStateful{Outgoing} && P2 <: AntiFermionStateful{Incoming})
-            return false
-        end
-    end
-
-    return true
-end
-
-function type_index_from_name(::QEDModel, name::String)
-    if startswith(name, "ki")
-        return (PhotonStateful{Incoming, PolX}, parse(Int, name[3:end]))
-    elseif startswith(name, "ko")
-        return (PhotonStateful{Outgoing, PolX}, parse(Int, name[3:end]))
-    elseif startswith(name, "ei")
-        return (FermionStateful{Incoming, SpinUp}, parse(Int, name[3:end]))
-    elseif startswith(name, "eo")
-        return (FermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end]))
-    elseif startswith(name, "pi")
-        return (AntiFermionStateful{Incoming, SpinUp}, parse(Int, name[3:end]))
-    elseif startswith(name, "po")
-        return (AntiFermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end]))
-    else
-        throw("Invalid name for a particle in the QED model")
-    end
-end
-
-"""
-    issame(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
-
-For two given [`QEDParticle`](@ref) types, return whether they are equivalent for the purpose of a Feynman Diagram. That means e.g. an `Incoming` `AntiFermion` is the same as an `Outgoing` `Fermion`. This is equivalent to `!caninteract(T1, T2)`.
-
-See also: [`caninteract`](@ref) and [`interaction_result`](@ref)
-"""
-function issame(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
-    return !caninteract(T1, T2)
-end
-
-"""
-    QED_vertex()
-
-Return the factor of a vertex in a QED feynman diagram.
-"""
-@inline function QED_vertex()::SLorentzVector{DiracMatrix}
-    # Peskin-Schroeder notation
-    return -1im * e * gamma()
-end
-
-@inline function QED_inner_edge(p::QEDParticle)
-    return propagator(particle(p), p.momentum)
-end
-
-"""
-    QED_conserve_momentum(p1::QEDParticle, p2::QEDParticle)
-
-Calculate and return a new particle from two given interacting ones at a vertex.
-"""
-function QED_conserve_momentum(
-    p1::P1,
-    p2::P2,
-) where {
-    Dir1 <: ParticleDirection,
-    Dir2 <: ParticleDirection,
-    SpinPol1 <: AbstractSpinOrPolarization,
-    SpinPol2 <: AbstractSpinOrPolarization,
-    P1 <: Union{FermionStateful{Dir1, SpinPol1}, AntiFermionStateful{Dir1, SpinPol1}, PhotonStateful{Dir1, SpinPol1}},
-    P2 <: Union{FermionStateful{Dir2, SpinPol2}, AntiFermionStateful{Dir2, SpinPol2}, PhotonStateful{Dir2, SpinPol2}},
-}
-    P3 = interaction_result(P1, P2)
-    p1_mom = p1.momentum
-    if (Dir1 <: Outgoing)
-        p1_mom *= -1
-    end
-    p2_mom = p2.momentum
-    if (Dir2 <: Outgoing)
-        p2_mom *= -1
-    end
-
-    p3_mom = p1_mom + p2_mom
-    if (typeof(direction(P3)) <: Incoming)
-        return P3(-p3_mom)
-    end
-    return P3(p3_mom)
-end
-
-"""
-    QEDProcessInput <: AbstractProcessInput
-
-Input for a QED Process. Contains the [`QEDProcessDescription`](@ref) of the process it is an input for, and the values of the in and out particles.
-
-See also: [`gen_process_input`](@ref)
-"""
-struct QEDProcessInput{N1, N2, N3, N4, N5, N6} <: AbstractProcessInput
-    process::QEDProcessDescription
-    inFerms::SVector{N1, FermionStateful{Incoming, SpinUp}}
-    outFerms::SVector{N2, FermionStateful{Outgoing, SpinUp}}
-    inAntiferms::SVector{N3, AntiFermionStateful{Incoming, SpinUp}}
-    outAntiferms::SVector{N4, AntiFermionStateful{Outgoing, SpinUp}}
-    inPhotons::SVector{N5, PhotonStateful{Incoming, PolX}}
-    outPhotons::SVector{N6, PhotonStateful{Outgoing, PolX}}
-end
-
-"""
-    model(::AbstractProcessDescription)
-
-Return the model of this process description.
-"""
-model(::QEDProcessDescription) = QEDModel()
-model(::QEDProcessInput) = QEDModel()
-
-function copy(process::QEDProcessDescription)
-    return QEDProcessDescription(copy(process.inParticles), copy(process.outParticles))
-end
-
-==(p1::QEDProcessDescription, p2::QEDProcessDescription) =
-    p1.inParticles == p2.inParticles && p1.outParticles == p2.outParticles
-
-function in_particles(process::QEDProcessDescription)
-    return process.inParticles
-end
-
-function out_particles(process::QEDProcessDescription)
-    return process.outParticles
-end
-
-function get_particle(input::QEDProcessInput, t::Type{Particle}, n::Int)::Particle where {Particle}
-    if (t <: FermionStateful{Incoming})
-        return input.inFerms[n]
-    elseif (t <: FermionStateful{Outgoing})
-        return input.outFerms[n]
-    elseif (t <: AntiFermionStateful{Incoming})
-        return input.inAntiferms[n]
-    elseif (t <: AntiFermionStateful{Outgoing})
-        return input.outAntiferms[n]
-    elseif (t <: PhotonStateful{Incoming})
-        return input.inPhotons[n]
-    elseif (t <: PhotonStateful{Outgoing})
-        return input.outPhotons[n]
-    end
-    @assert false "Invalid type given"
-end
--- a/src/scheduler/greedy.jl
+++ b/src/scheduler/greedy.jl
@ -4,7 +4,7 @@

 A greedy implementation of a scheduler, creating a topological ordering of nodes and naively balancing them onto the different devices.
 """
-struct GreedyScheduler end
+struct GreedyScheduler <: AbstractScheduler end

 function schedule_dag(::GreedyScheduler, graph::DAG, machine::Machine)
    nodeQueue = PriorityQueue{Node, Int}()
--- a/src/scheduler/interface.jl
+++ b/src/scheduler/interface.jl
@ -1,10 +1,10 @@

 """
-    Scheduler
+    AbstractScheduler

 Abstract base type for scheduler implementations. The scheduler is used to assign each node to a device and create a topological ordering of tasks.
 """
-abstract type Scheduler end
+abstract type AbstractScheduler end

 """
    schedule_dag(::Scheduler, ::DAG, ::Machine)
--- a/src/scheduler/type.jl
+++ b/src/scheduler/type.jl
@ -5,10 +5,11 @@ using StaticArrays

 Type representing a function call with `N` parameters. Contains the function to call, argument symbols, the return symbol and the device to execute on.
 """
-struct FunctionCall{VectorType <: AbstractVector, M}
+struct FunctionCall{VectorType <: AbstractVector, N}
    func::Function
-    arguments::VectorType
-    additional_arguments::SVector{M, Any} # additional arguments (as values) for the function call, will be prepended to the other arguments
+    # TODO: this should be a tuple
+    value_arguments::SVector{N, Any}    # value arguments for the function call, will be prepended to the other arguments
+    arguments::VectorType               # symbols of the inputs to the function call
    return_symbol::Symbol
    device::AbstractDevice
 end
--- a/src/task/compute.jl
+++ b/src/task/compute.jl
@ -32,7 +32,7 @@ function get_function_call(
    inSymbols::AbstractVector,
    outSymbol::Symbol,
 ) where {CompTask <: AbstractComputeTask}
-    return [FunctionCall(compute, inSymbols, SVector{1, Any}(t), outSymbol, device)]
+    return [FunctionCall(compute, SVector{1, Any}(t), inSymbols, outSymbol, device)]
 end

 function get_function_call(node::ComputeTaskNode)
@ -64,8 +64,8 @@ function get_function_call(node::DataTaskNode)
    return [
        FunctionCall(
            unpack_identity,
-            SVector{1, Symbol}(Symbol(to_var_name(first(children(node)).id))),
            SVector{0, Any}(),
+            SVector{1, Symbol}(Symbol(to_var_name(first(children(node)).id))),
            Symbol(to_var_name(node.id)),
            first(children(node)).device,
        ),
@ -77,8 +77,8 @@ function get_init_function_call(node::DataTaskNode, device::AbstractDevice)

    return FunctionCall(
        unpack_identity,
-        SVector{1, Symbol}(Symbol("$(to_var_name(node.id))_in")),
        SVector{0, Any}(),
+        SVector{1, Symbol}(Symbol("$(to_var_name(node.id))_in")),
        Symbol(to_var_name(node.id)),
        device,
    )
--- a/src/task/properties.jl
+++ b/src/task/properties.jl
@ -3,28 +3,21 @@

 Fallback implementation of the compute function of a compute task, throwing an error.
 """
-function compute(t::AbstractTask, data...)
-    return error("Need to implement compute()")
-end
+function compute end

 """
    compute_effort(t::AbstractTask)

 Fallback implementation of the compute effort of a task, throwing an error.
 """
-function compute_effort(t::AbstractTask)::Float64
-    # default implementation using compute
-    return error("Need to implement compute_effort()")
-end
+function compute_effort end

 """
    data(t::AbstractTask)

 Fallback implementation of the data of a task, throwing an error.
 """
-function data(t::AbstractTask)::Float64
-    return error("Need to implement data()")
-end
+function data end

 """
    compute_effort(t::AbstractDataTask)
--- a/src/utility.jl
+++ b/src/utility.jl
@ -1,3 +1,6 @@
+using Roots
+using ForwardDiff
+
 """
    noop()

@ -249,7 +252,6 @@ end

 first_derivative(func) = x -> ForwardDiff.derivative(func, float(x))

-
 function generate_physical_massive_moms(rng, ss, masses; x0 = 0.1)
    n = length(masses)
    massless_moms = generate_physical_massless_moms(rng, ss, n)
--- a/test/Project.toml
+++ b/test/Project.toml
@ -1,10 +0,0 @@
-[deps]
-AccurateArithmetic = "22286c92-06ac-501d-9306-4abd417d9753"
-QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
-QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
--- a/test/runtests.jl
+++ b/test/runtests.jl
@ -1,5 +1,6 @@
 using SafeTestsets

+#=
@safetestset "Utility Unit Tests                  " begin
    include("unit_tests_utility.jl")
 end
@ -18,9 +19,11 @@ end
@safetestset "ABC-Model Unit Tests                " begin
    include("unit_tests_abcmodel.jl")
 end
+=#
@safetestset "QED-Model Unit Tests                " begin
    include("unit_tests_qedmodel.jl")
 end
+#=
@safetestset "QED Feynman Diagram Generation Tests" begin
    include("unit_tests_qed_diagrams.jl")
 end
@ -39,3 +42,4 @@ end
@safetestset "Known Graph Tests                   " begin
    include("known_graphs.jl")
 end
+=#
--- a/test/unit_tests_qed_diagrams.jl
+++ b/test/unit_tests_qed_diagrams.jl
@ -22,7 +22,7 @@ trident = ("Trident", parse_process("ke->epe", model), 8)

        n_particles = 0
        for type in types(model)
-            if (isincoming(type))
+            if (is_incoming(type))
                n_particles += get(process.inParticles, type, 0)
            else
                n_particles += get(process.outParticles, type, 0)
--- a/test/unit_tests_qedmodel.jl
+++ b/test/unit_tests_qedmodel.jl
@ -1,5 +1,6 @@
 using MetagraphOptimization
 using QEDbase
+using QEDcore
 using QEDprocesses
 using StatsBase # for countmap
 using Random
@ -9,8 +10,6 @@ import MetagraphOptimization.caninteract
 import MetagraphOptimization.issame
 import MetagraphOptimization.interaction_result
 import MetagraphOptimization.propagation_result
-import MetagraphOptimization.direction
-import MetagraphOptimization.spin_or_pol
 import MetagraphOptimization.QED_vertex

 def_momentum = SFourMomentum(1.0, 0.0, 0.0, 0.0)
@ -18,32 +17,32 @@ def_momentum = SFourMomentum(1.0, 0.0, 0.0, 0.0)
 RNG = Random.MersenneTwister(0)

 testparticleTypes = [
-    PhotonStateful{Incoming, PolX},
-    PhotonStateful{Outgoing, PolX},
-    FermionStateful{Incoming, SpinUp},
-    FermionStateful{Outgoing, SpinUp},
-    AntiFermionStateful{Incoming, SpinUp},
-    AntiFermionStateful{Outgoing, SpinUp},
+    ParticleStateful{Incoming, Photon, SFourMomentum},
+    ParticleStateful{Outgoing, Photon, SFourMomentum},
+    ParticleStateful{Incoming, Electron, SFourMomentum},
+    ParticleStateful{Outgoing, Electron, SFourMomentum},
+    ParticleStateful{Incoming, Positron, SFourMomentum},
+    ParticleStateful{Outgoing, Positron, SFourMomentum},
 ]

 testparticleTypesPropagated = [
-    PhotonStateful{Outgoing, PolX},
-    PhotonStateful{Incoming, PolX},
-    FermionStateful{Outgoing, SpinUp},
-    FermionStateful{Incoming, SpinUp},
-    AntiFermionStateful{Outgoing, SpinUp},
-    AntiFermionStateful{Incoming, SpinUp},
+    ParticleStateful{Outgoing, Photon, SFourMomentum},
+    ParticleStateful{Incoming, Photon, SFourMomentum},
+    ParticleStateful{Outgoing, Electron, SFourMomentum},
+    ParticleStateful{Incoming, Electron, SFourMomentum},
+    ParticleStateful{Outgoing, Positron, SFourMomentum},
+    ParticleStateful{Incoming, Positron, SFourMomentum},
 ]

-function compton_groundtruth(input::QEDProcessInput)
+function compton_groundtruth(input::PhaseSpacePoint)
    # p1k1 -> p2k2
    # formula: −(ie)^2 (u(p2) slashed(ε1) S(p2 − k1) slashed(ε2) u(p1) + u(p2) slashed(ε2) S(p1 + k1) slashed(ε1) u(p1))

-    p1 = input.inFerms[1]
-    p2 = input.outFerms[1]
+    p1 = momentum(psp, Incoming(), 2)
+    p2 = momentum(psp, Outgoing(), 2)

-    k1 = input.inPhotons[1]
-    k2 = input.outPhotons[1]
+    k1 = momentum(psp, Incoming(), 1)
+    k2 = momentum(psp, Outgoing(), 1)

    u_p1 = base_state(Electron(), Incoming(), p1.momentum, spin_or_pol(p1))
    u_p2 = base_state(Electron(), Outgoing(), p2.momentum, spin_or_pol(p2))
@ -57,8 +56,8 @@ function compton_groundtruth(input::QEDProcessInput)
    virt2_mom = p1.momentum + k1.momentum
    @test isapprox(p2.momentum + k2.momentum, virt2_mom)

-    s_p2_k1 = propagator(Electron(), virt1_mom)
-    s_p1_k1 = propagator(Electron(), virt2_mom)
+    s_p2_k1 = QEDbase.propagator(Electron(), virt1_mom)
+    s_p1_k1 = QEDbase.propagator(Electron(), virt2_mom)

    diagram1 = u_p2 * (eps_1 * QED_vertex()) * s_p2_k1 * (eps_2 * QED_vertex()) * u_p1
    diagram2 = u_p2 * (eps_2 * QED_vertex()) * s_p1_k1 * (eps_1 * QED_vertex()) * u_p1
@ -66,7 +65,6 @@ function compton_groundtruth(input::QEDProcessInput)
    return diagram1 + diagram2
 end

-
@testset "Interaction Result" begin
    import MetagraphOptimization.QED_conserve_momentum

@ -88,8 +86,8 @@ end
        @test issame(typeof(resultParticle), interaction_result(p1, p2))

        totalMom = zero(SFourMomentum)
-        for (p, mom) in [(p1, testParticle1.momentum), (p2, testParticle2.momentum), (p3, resultParticle.momentum)]
-            if (typeof(direction(p)) <: Incoming)
+        for (p, mom) in [(p1, momentum(testParticle1)), (p2, momentum(testParticle2)), (p3, momentum(resultParticle))]
+            if (typeof(particle_direction(p)) <: Incoming)
                totalMom += mom
            else
                totalMom -= mom
@ -103,7 +101,7 @@ end
@testset "Propagation Result" begin
    for (p, propResult) in zip(testparticleTypes, testparticleTypesPropagated)
        @test issame(propagation_result(p), propResult)
-        @test direction(propagation_result(p)(def_momentum)) != direction(p(def_momentum))
+        @test particle_direction(propagation_result(p)(def_momentum)) != particle_direction(p(def_momentum))
    end
 end

@ -117,38 +115,23 @@ end
    end

    @testset "Known processes" begin
-        compton_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, FermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 1, FermionStateful{Outgoing, SpinUp} => 1),
-        )
+        compton_process = GenericQEDProcess(1, 1, 1, 1, 0, 0)

        @test parse_process("ke->ke", QEDModel()) == compton_process

-        positron_compton_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, AntiFermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 1, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        positron_compton_process = GenericQEDProcess(1, 1, 0, 0, 1, 1)

        @test parse_process("kp->kp", QEDModel()) == positron_compton_process

-        trident_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, FermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(FermionStateful{Outgoing, SpinUp} => 2, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        trident_process = GenericQEDProcess(1, 0, 1, 2, 0, 1)

        @test parse_process("ke->eep", QEDModel()) == trident_process

-        pair_production_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 2),
-            Dict{Type, Int}(FermionStateful{Outgoing, SpinUp} => 1, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        pair_production_process = GenericQEDProcess(2, 0, 0, 1, 0, 1)

        @test parse_process("kk->pe", QEDModel()) == pair_production_process

-        pair_annihilation_process = QEDProcessDescription(
-            Dict{Type, Int}(FermionStateful{Incoming, SpinUp} => 1, AntiFermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 2),
-        )
+        pair_annihilation_process = GenericQEDProcess(0, 2, 1, 0, 1, 0)

        @test parse_process("pe->kk", QEDModel()) == pair_annihilation_process
    end
@ -161,12 +144,18 @@ end

        for i in 1:100
            input = gen_process_input(process)
-            @test length(input.inFerms) == get(process.inParticles, FermionStateful{Incoming, SpinUp}, 0)
-            @test length(input.inAntiferms) == get(process.inParticles, AntiFermionStateful{Incoming, SpinUp}, 0)
-            @test length(input.inPhotons) == get(process.inParticles, PhotonStateful{Incoming, PolX}, 0)
-            @test length(input.outFerms) == get(process.outParticles, FermionStateful{Outgoing, SpinUp}, 0)
-            @test length(input.outAntiferms) == get(process.outParticles, AntiFermionStateful{Outgoing, SpinUp}, 0)
-            @test length(input.outPhotons) == get(process.outParticles, PhotonStateful{Outgoing, PolX}, 0)
+            @test length(input.inFerms) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Electron, SFourMomentum}, 0)
+            @test length(input.inAntiferms) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Positron, SFourMomentum}, 0)
+            @test length(input.inPhotons) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Photon, SFourMomentum}, 0)
+            @test length(input.outFerms) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Electron, SFourMomentum}, 0)
+            @test length(input.outAntiferms) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Positron, SFourMomentum}, 0)
+            @test length(input.outPhotons) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Photon, SFourMomentum}, 0)

            @test isapprox(
                sum([
@ -185,6 +174,7 @@ end
    end
 end

+#=
@testset "Compton" begin
    import MetagraphOptimization.insert_node!
    import MetagraphOptimization.insert_edge!
@ -211,97 +201,97 @@ end
    graph = DAG()

    # s to output (exit node)
-    d_exit = insert_node!(graph, make_node(DataTask(16)), track = false)
+    d_exit = insert_node!(graph, make_node(DataTask(16)); track=false)

-    sum_node = insert_node!(graph, make_node(ComputeTaskQED_Sum(2)), track = false)
+    sum_node = insert_node!(graph, make_node(ComputeTaskQED_Sum(2)); track=false)

-    d_s0_sum = insert_node!(graph, make_node(DataTask(16)), track = false)
-    d_s1_sum = insert_node!(graph, make_node(DataTask(16)), track = false)
+    d_s0_sum = insert_node!(graph, make_node(DataTask(16)); track=false)
+    d_s1_sum = insert_node!(graph, make_node(DataTask(16)); track=false)

    # final s compute
-    s0 = insert_node!(graph, make_node(ComputeTaskQED_S2()), track = false)
-    s1 = insert_node!(graph, make_node(ComputeTaskQED_S2()), track = false)
+    s0 = insert_node!(graph, make_node(ComputeTaskQED_S2()); track=false)
+    s1 = insert_node!(graph, make_node(ComputeTaskQED_S2()); track=false)

    # data from v0 and v1 to s0
-    d_v0_s0 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_v1_s0 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_v2_s1 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_v3_s1 = insert_node!(graph, make_node(DataTask(96)), track = false)
+    d_v0_s0 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_v1_s0 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_v2_s1 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_v3_s1 = insert_node!(graph, make_node(DataTask(96)); track=false)

    # v0 and v1 compute
-    v0 = insert_node!(graph, make_node(ComputeTaskQED_V()), track = false)
-    v1 = insert_node!(graph, make_node(ComputeTaskQED_V()), track = false)
-    v2 = insert_node!(graph, make_node(ComputeTaskQED_V()), track = false)
-    v3 = insert_node!(graph, make_node(ComputeTaskQED_V()), track = false)
+    v0 = insert_node!(graph, make_node(ComputeTaskQED_V()); track=false)
+    v1 = insert_node!(graph, make_node(ComputeTaskQED_V()); track=false)
+    v2 = insert_node!(graph, make_node(ComputeTaskQED_V()); track=false)
+    v3 = insert_node!(graph, make_node(ComputeTaskQED_V()); track=false)

    # data from uPhIn, uPhOut, uElIn, uElOut to v0 and v1
-    d_uPhIn_v0 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uElIn_v0 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uPhOut_v1 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uElOut_v1 = insert_node!(graph, make_node(DataTask(96)), track = false)
+    d_uPhIn_v0 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uElIn_v0 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uPhOut_v1 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uElOut_v1 = insert_node!(graph, make_node(DataTask(96)); track=false)

    # data from uPhIn, uPhOut, uElIn, uElOut to v2 and v3
-    d_uPhOut_v2 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uElIn_v2 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uPhIn_v3 = insert_node!(graph, make_node(DataTask(96)), track = false)
-    d_uElOut_v3 = insert_node!(graph, make_node(DataTask(96)), track = false)
+    d_uPhOut_v2 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uElIn_v2 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uPhIn_v3 = insert_node!(graph, make_node(DataTask(96)); track=false)
+    d_uElOut_v3 = insert_node!(graph, make_node(DataTask(96)); track=false)

    # uPhIn, uPhOut, uElIn and uElOut computes
-    uPhIn = insert_node!(graph, make_node(ComputeTaskQED_U()), track = false)
-    uPhOut = insert_node!(graph, make_node(ComputeTaskQED_U()), track = false)
-    uElIn = insert_node!(graph, make_node(ComputeTaskQED_U()), track = false)
-    uElOut = insert_node!(graph, make_node(ComputeTaskQED_U()), track = false)
+    uPhIn = insert_node!(graph, make_node(ComputeTaskQED_U()); track=false)
+    uPhOut = insert_node!(graph, make_node(ComputeTaskQED_U()); track=false)
+    uElIn = insert_node!(graph, make_node(ComputeTaskQED_U()); track=false)
+    uElOut = insert_node!(graph, make_node(ComputeTaskQED_U()); track=false)

    # data into U
-    d_uPhIn = insert_node!(graph, make_node(DataTask(16), "ki1"), track = false)
-    d_uPhOut = insert_node!(graph, make_node(DataTask(16), "ko1"), track = false)
-    d_uElIn = insert_node!(graph, make_node(DataTask(16), "ei1"), track = false)
-    d_uElOut = insert_node!(graph, make_node(DataTask(16), "eo1"), track = false)
+    d_uPhIn = insert_node!(graph, make_node(DataTask(16), "ki1"); track=false)
+    d_uPhOut = insert_node!(graph, make_node(DataTask(16), "ko1"); track=false)
+    d_uElIn = insert_node!(graph, make_node(DataTask(16), "ei1"); track=false)
+    d_uElOut = insert_node!(graph, make_node(DataTask(16), "eo1"); track=false)

    # now for all the edges
-    insert_edge!(graph, d_uPhIn, uPhIn, track = false)
-    insert_edge!(graph, d_uPhOut, uPhOut, track = false)
-    insert_edge!(graph, d_uElIn, uElIn, track = false)
-    insert_edge!(graph, d_uElOut, uElOut, track = false)
+    insert_edge!(graph, d_uPhIn, uPhIn; track=false)
+    insert_edge!(graph, d_uPhOut, uPhOut; track=false)
+    insert_edge!(graph, d_uElIn, uElIn; track=false)
+    insert_edge!(graph, d_uElOut, uElOut; track=false)

-    insert_edge!(graph, uPhIn, d_uPhIn_v0, track = false)
-    insert_edge!(graph, uPhOut, d_uPhOut_v1, track = false)
-    insert_edge!(graph, uElIn, d_uElIn_v0, track = false)
-    insert_edge!(graph, uElOut, d_uElOut_v1, track = false)
+    insert_edge!(graph, uPhIn, d_uPhIn_v0; track=false)
+    insert_edge!(graph, uPhOut, d_uPhOut_v1; track=false)
+    insert_edge!(graph, uElIn, d_uElIn_v0; track=false)
+    insert_edge!(graph, uElOut, d_uElOut_v1; track=false)

-    insert_edge!(graph, uPhIn, d_uPhIn_v3, track = false)
-    insert_edge!(graph, uPhOut, d_uPhOut_v2, track = false)
-    insert_edge!(graph, uElIn, d_uElIn_v2, track = false)
-    insert_edge!(graph, uElOut, d_uElOut_v3, track = false)
+    insert_edge!(graph, uPhIn, d_uPhIn_v3; track=false)
+    insert_edge!(graph, uPhOut, d_uPhOut_v2; track=false)
+    insert_edge!(graph, uElIn, d_uElIn_v2; track=false)
+    insert_edge!(graph, uElOut, d_uElOut_v3; track=false)

-    insert_edge!(graph, d_uPhIn_v0, v0, track = false)
-    insert_edge!(graph, d_uPhOut_v1, v1, track = false)
-    insert_edge!(graph, d_uElIn_v0, v0, track = false)
-    insert_edge!(graph, d_uElOut_v1, v1, track = false)
+    insert_edge!(graph, d_uPhIn_v0, v0; track=false)
+    insert_edge!(graph, d_uPhOut_v1, v1; track=false)
+    insert_edge!(graph, d_uElIn_v0, v0; track=false)
+    insert_edge!(graph, d_uElOut_v1, v1; track=false)

-    insert_edge!(graph, d_uPhIn_v3, v3, track = false)
-    insert_edge!(graph, d_uPhOut_v2, v2, track = false)
-    insert_edge!(graph, d_uElIn_v2, v2, track = false)
-    insert_edge!(graph, d_uElOut_v3, v3, track = false)
+    insert_edge!(graph, d_uPhIn_v3, v3; track=false)
+    insert_edge!(graph, d_uPhOut_v2, v2; track=false)
+    insert_edge!(graph, d_uElIn_v2, v2; track=false)
+    insert_edge!(graph, d_uElOut_v3, v3; track=false)

-    insert_edge!(graph, v0, d_v0_s0, track = false)
-    insert_edge!(graph, v1, d_v1_s0, track = false)
-    insert_edge!(graph, v2, d_v2_s1, track = false)
-    insert_edge!(graph, v3, d_v3_s1, track = false)
+    insert_edge!(graph, v0, d_v0_s0; track=false)
+    insert_edge!(graph, v1, d_v1_s0; track=false)
+    insert_edge!(graph, v2, d_v2_s1; track=false)
+    insert_edge!(graph, v3, d_v3_s1; track=false)

-    insert_edge!(graph, d_v0_s0, s0, track = false)
-    insert_edge!(graph, d_v1_s0, s0, track = false)
+    insert_edge!(graph, d_v0_s0, s0; track=false)
+    insert_edge!(graph, d_v1_s0, s0; track=false)

-    insert_edge!(graph, d_v2_s1, s1, track = false)
-    insert_edge!(graph, d_v3_s1, s1, track = false)
+    insert_edge!(graph, d_v2_s1, s1; track=false)
+    insert_edge!(graph, d_v3_s1, s1; track=false)

-    insert_edge!(graph, s0, d_s0_sum, track = false)
-    insert_edge!(graph, s1, d_s1_sum, track = false)
+    insert_edge!(graph, s0, d_s0_sum; track=false)
+    insert_edge!(graph, s1, d_s1_sum; track=false)

-    insert_edge!(graph, d_s0_sum, sum_node, track = false)
-    insert_edge!(graph, d_s1_sum, sum_node, track = false)
+    insert_edge!(graph, d_s0_sum, sum_node; track=false)
+    insert_edge!(graph, d_s1_sum, sum_node; track=false)

-    insert_edge!(graph, sum_node, d_exit, track = false)
+    insert_edge!(graph, sum_node, d_exit; track=false)

    input = [gen_process_input(process) for _ in 1:1000]

@ -314,9 +304,12 @@ end
    @test isapprox(compton_function.(input), compton_groundtruth.(input))
 end

-@testset "Equal results after optimization" for optimizer in
-                                                [ReductionOptimizer(), RandomWalkOptimizer(MersenneTwister(0))]
-    @testset "Process $proc_str" for proc_str in ["ke->ke", "kp->kp", "kk->ep", "ep->kk", "ke->kke", "ke->kkke"]
+@testset "Equal results after optimization" for optimizer in [
+    ReductionOptimizer(), RandomWalkOptimizer(MersenneTwister(0))
+]
+    @testset "Process $proc_str" for proc_str in [
+        "ke->ke", "kp->kp", "kk->ep", "ep->kk", "ke->kke", "ke->kkke"
+    ]
        model = QEDModel()
        process = parse_process(proc_str, model)
        machine = Machine(
@ -347,3 +340,4 @@ end
        @test isapprox(compute_function.(input), reduced_compute_function.(input))
    end
 end
+=#
Author	SHA1	Message	Date
Anton Reinhard	5be7ca99e7	Add results and evaluation Some checks failed MetagraphOptimization_CI / test (push) Failing after 1m33s Details MetagraphOptimization_CI / docs (push) Failing after 1m33s Details	2024-07-10 14:21:26 +02:00
rubydragon	1ae39a8caa	Congruent in photons example (#12 ) Some checks failed MetagraphOptimization_CI / docs (push) Failing after 1m50s Details MetagraphOptimization_CI / test (push) Failing after 1m53s Details Now targeting the correct branches Co-authored-by: Rubydragon <anton.reinhard@proton.me> Reviewed-on: #12	2024-07-10 14:17:39 +02:00
Anton Reinhard	b5d92b729c	Get compute function working Some checks failed MetagraphOptimization_CI / docs (push) Failing after 1m32s Details MetagraphOptimization_CI / test (push) Failing after 1m49s Details	2024-07-04 15:31:22 +02:00
Anton Reinhard	6a9a7b41f1	rework a lot of the QED model to use QEDcore/base/processes	2024-07-03 20:24:53 +02:00
Anton Reinhard	a1581182ca	WIP	2024-07-02 10:50:30 +02:00
Anton Reinhard	1b4ba285c3	WIP refactor Some checks failed MetagraphOptimization_CI / test (push) Failing after 7m59s Details MetagraphOptimization_CI / docs (push) Failing after 8m6s Details	2024-06-24 23:31:30 +02:00
Anton Reinhard	2921882fd4	EOD Some checks failed MetagraphOptimization_CI / docs (push) Failing after 8m32s Details MetagraphOptimization_CI / test (push) Failing after 11m24s Details	2024-05-24 19:20:59 +02:00