Congruent in photons example (#12)
Now targeting the correct branches Co-authored-by: Rubydragon <anton.reinhard@proton.me> Reviewed-on: #12
This commit is contained in:
parent
b5d92b729c
commit
1ae39a8caa
194
examples/congruent_in_ph.jl
Normal file
194
examples/congruent_in_ph.jl
Normal file
@ -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
|
||||
@ -100,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.==
|
||||
|
||||
@ -1,25 +1,5 @@
|
||||
# patch QEDprocesses
|
||||
# see issue https://github.com/QEDjl-project/QEDprocesses.jl/issues/77
|
||||
@inline function QEDprocesses.number_particles(
|
||||
proc_def::QEDbase.AbstractProcessDefinition,
|
||||
dir::DIR,
|
||||
::PT,
|
||||
) where {DIR <: QEDbase.ParticleDirection, PT <: QEDbase.AbstractParticleType}
|
||||
return count(x -> x isa PT, particles(proc_def, dir))
|
||||
end
|
||||
|
||||
@inline function QEDprocesses.number_particles(
|
||||
proc_def::QEDbase.AbstractProcessDefinition,
|
||||
::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 QEDprocesses.number_particles(
|
||||
proc_def::QEDbase.AbstractProcessDefinition,
|
||||
::Type{PS},
|
||||
@ -43,29 +23,3 @@ end
|
||||
) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType, EL <: AbstractFourMomentum}
|
||||
return ParticleStateful(DIR(), SPECIES(), mom)
|
||||
end
|
||||
|
||||
@inline function QEDbase.momentum(
|
||||
psp::AbstractPhaseSpacePoint{MODEL, PROC, PS_DEF, INT, OUTT},
|
||||
dir::ParticleDirection,
|
||||
species::AbstractParticleType,
|
||||
n::Int,
|
||||
) where {MODEL, PROC, PS_DEF, INT, OUTT}
|
||||
# TODO: can be done through fancy template recursion too with 0 overhead
|
||||
i = 0
|
||||
c = n
|
||||
for p in particles(psp, dir)
|
||||
i += 1
|
||||
if particle_species(p) isa typeof(species)
|
||||
c -= 1
|
||||
end
|
||||
if c == 0
|
||||
break
|
||||
end
|
||||
end
|
||||
|
||||
if c != 0 || n <= 0
|
||||
throw(InvalidInputError("could not get $n-th momentum of $dir $species, does not exist"))
|
||||
end
|
||||
|
||||
return momenta(psp, dir)[i]
|
||||
end
|
||||
|
||||
@ -7,7 +7,7 @@ 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]))
|
||||
@ -30,7 +30,7 @@ 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]))
|
||||
|
||||
@ -75,7 +75,7 @@ end
|
||||
inputSymbols::Dict{String, Vector{Symbol}},
|
||||
instance::AbstractProblemInstance,
|
||||
machine::Machine,
|
||||
problemInputSymbol::Symbol = :input,
|
||||
problemInputSymbol::Symbol = :data_input,
|
||||
)
|
||||
|
||||
Return a `Vector{Expr}` doing the input assignments from the given `problemInputSymbol` onto the `inputSymbols`.
|
||||
@ -84,9 +84,9 @@ function gen_input_assignment_code(
|
||||
inputSymbols::Dict{String, Vector{Symbol}},
|
||||
instance,
|
||||
machine::Machine,
|
||||
problemInputSymbol::Symbol = :input,
|
||||
problemInputSymbol::Symbol = :data_input,
|
||||
)
|
||||
assignInputs = Vector{FunctionCall}()
|
||||
assignInputs = Vector{Expr}()
|
||||
for (name, symbols) in inputSymbols
|
||||
# make a function for this, since we can't use anonymous functions in the FunctionCall
|
||||
|
||||
@ -94,7 +94,9 @@ function gen_input_assignment_code(
|
||||
device = entry_device(machine)
|
||||
push!(
|
||||
assignInputs,
|
||||
FunctionCall(input_expr, SVector{1, Any}(name), SVector{1, Symbol}(problemInputSymbol), symbol, device),
|
||||
Meta.parse(
|
||||
"$(eval(gen_access_expr(device, symbol))) = $(input_expr(instance, name, problemInputSymbol))",
|
||||
),
|
||||
)
|
||||
end
|
||||
end
|
||||
@ -126,7 +128,7 @@ function gen_tape(graph::DAG, instance, machine::Machine, scheduler::AbstractSch
|
||||
outSym = Symbol(to_var_name(get_exit_node(graph).id))
|
||||
|
||||
initCaches = gen_cache_init_code(machine)
|
||||
assignInputs = gen_input_assignment_code(inputSyms, instance, machine, :input)
|
||||
assignInputs = gen_input_assignment_code(inputSyms, instance, machine, :data_input)
|
||||
|
||||
return Tape{input_type(instance)}(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine)
|
||||
end
|
||||
@ -140,7 +142,7 @@ For implementation reasons, this disregards the set [`CacheStrategy`](@ref) of t
|
||||
"""
|
||||
function execute_tape(tape::Tape, input)
|
||||
cache = Dict{Symbol, Any}()
|
||||
cache[:input] = input
|
||||
cache[:data_input] = input
|
||||
# simply execute all the code snippets here
|
||||
@assert typeof(input) == input_type(tape.instance)
|
||||
# TODO: `@assert` that input fits the tape.instance
|
||||
|
||||
@ -11,7 +11,7 @@ TODO: update docs
|
||||
"""
|
||||
struct Tape{INPUT}
|
||||
initCachesCode::Vector{Expr}
|
||||
inputAssignCode::Vector{FunctionCall}
|
||||
inputAssignCode::Vector{Expr}
|
||||
computeCode::Vector{FunctionCall}
|
||||
inputSymbols::Dict{String, Vector{Symbol}}
|
||||
outputSymbol::Symbol
|
||||
|
||||
@ -39,8 +39,8 @@ Generate the [`DAG`](@ref) for the given [`AbstractProblemInstance`](@ref). Ever
|
||||
function graph end
|
||||
|
||||
"""
|
||||
input_expr(::AbstractProblemInstance, name::String, input)
|
||||
input_expr(instance::AbstractProblemInstance, name::String, input_symbol::Symbol)
|
||||
|
||||
For a given [`AbstractProblemInstance`](@ref), the problem input (of type `input_type(...)`) and an entry node name, return a that specific input value from the input symbol.
|
||||
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 input_expr end
|
||||
|
||||
@ -1,12 +1,26 @@
|
||||
using StaticArrays
|
||||
|
||||
function input_expr(name::String, psp::PhaseSpacePoint)
|
||||
construction_string(::Incoming) = "Incoming()"
|
||||
construction_string(::Outgoing) = "Outgoing()"
|
||||
|
||||
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 ParticleValueSP(
|
||||
type(momentum(psp, particle_direction(type), particle_species(type), index)),
|
||||
0.0im,
|
||||
spin_or_pol(process(psp), type, index),
|
||||
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
|
||||
|
||||
|
||||
@ -8,7 +8,6 @@ function _svector_from_type(processDescription::GenericQEDProcess, type, particl
|
||||
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
|
||||
|
||||
"""
|
||||
@ -64,9 +63,9 @@ function gen_graph(process_description::GenericQEDProcess)
|
||||
|
||||
# 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()
|
||||
@ -75,16 +74,16 @@ function gen_graph(process_description::GenericQEDProcess)
|
||||
# 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
|
||||
@ -100,19 +99,19 @@ function gen_graph(process_description::GenericQEDProcess)
|
||||
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
|
||||
@ -122,18 +121,18 @@ function gen_graph(process_description::GenericQEDProcess)
|
||||
|
||||
# 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
|
||||
@ -144,16 +143,16 @@ function gen_graph(process_description::GenericQEDProcess)
|
||||
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
|
||||
|
||||
|
||||
@ -152,8 +152,9 @@ 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{<:ParticleStateful})
|
||||
@ -503,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)
|
||||
|
||||
@ -511,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(ParticleStateful{Incoming, Photon}, i) for i in 1:n],
|
||||
[FeynmanParticle(ParticleStateful{Outgoing, Photon}, 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(
|
||||
@ -520,10 +520,10 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
|
||||
missing,
|
||||
[inFerm, outFerm, photons...],
|
||||
Dict{Type, Int64}(
|
||||
ParticleStateful{Incoming, Electron} => 1,
|
||||
ParticleStateful{Outgoing, Electron} => 1,
|
||||
ParticleStateful{Incoming, Photon} => n,
|
||||
ParticleStateful{Outgoing, Photon} => m,
|
||||
ParticleStateful{Incoming, Electron, SFourMomentum} => 1,
|
||||
ParticleStateful{Outgoing, Electron, SFourMomentum} => 1,
|
||||
ParticleStateful{Incoming, Photon, SFourMomentum} => n,
|
||||
ParticleStateful{Outgoing, Photon, SFourMomentum} => m,
|
||||
),
|
||||
)
|
||||
|
||||
@ -535,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(ParticleStateful{Incoming, Electron}, iterations),
|
||||
FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations),
|
||||
photons[order[left_index]],
|
||||
FeynmanParticle(ParticleStateful{Incoming, Electron}, iterations + 1),
|
||||
FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations + 1),
|
||||
)
|
||||
left_index += 1
|
||||
add_vertex!(new_diagram, v_left)
|
||||
@ -550,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(ParticleStateful{Outgoing, Electron}, iterations),
|
||||
FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations),
|
||||
photons[order[right_index]],
|
||||
FeynmanParticle(ParticleStateful{Outgoing, Electron}, iterations + 1),
|
||||
FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations + 1),
|
||||
)
|
||||
right_index -= 1
|
||||
add_vertex!(new_diagram, v_right)
|
||||
@ -566,7 +566,6 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
|
||||
add_tie!(new_diagram, FeynmanTie(ps[1], ps[2]))
|
||||
return new_diagram
|
||||
end
|
||||
=#
|
||||
|
||||
"""
|
||||
gen_compton_diagrams(n::Int, m::Int)
|
||||
@ -587,7 +586,6 @@ function gen_compton_diagrams(n::Int, m::Int)
|
||||
return vcat(diagrams...)
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
gen_compton_diagrams_one_side(n::Int, m::Int)
|
||||
|
||||
|
||||
@ -56,8 +56,8 @@ function compton_groundtruth(input::PhaseSpacePoint)
|
||||
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
|
||||
@ -65,7 +65,6 @@ function compton_groundtruth(input::PhaseSpacePoint)
|
||||
return diagram1 + diagram2
|
||||
end
|
||||
|
||||
|
||||
@testset "Interaction Result" begin
|
||||
import MetagraphOptimization.QED_conserve_momentum
|
||||
|
||||
@ -202,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]
|
||||
|
||||
@ -305,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(
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user