diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl new file mode 100644 index 0000000..a0639bf --- /dev/null +++ b/examples/congruent_in_ph.jl @@ -0,0 +1,123 @@ +using MetagraphOptimization +using QEDbase +using QEDcore +using Random +using UUIDs + +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(processDescription::QEDProcessDescription, omega::Number) + # generate an input sample for given e + nk -> e' + k' process, where the nk are equal + massSum = 0 + inputMasses = Vector{Float64}() + for (particle, n) in processDescription.inParticles + for _ in 1:n + massSum += mass(particle) + push!(inputMasses, mass(particle)) + end + end + outputMasses = Vector{Float64}() + for (particle, n) in processDescription.outParticles + for _ in 1:n + massSum += mass(particle) + push!(outputMasses, mass(particle)) + end + end + + initialMomenta = [ + i == 1 ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for + i in 1:length(inputMasses) + ] + + # add some extra random mass to allow for some momentum + ss = sqrt(sum(initialMomenta) * sum(initialMomenta)) + + result = Vector{QEDProcessInput}() + sizehint!(result, 16) + + spin_pol_combinations = Iterators.product( + [SpinUp, SpinDown], [SpinUp, SpinDown], [PolX, PolY], [PolX, PolY] + ) + for (in_spin, out_spin, in_pol, out_pol) in spin_pol_combinations + + # get the electron first, then the n photons + particles = Vector{QEDParticle}() + + for (particle, n) in processDescription.inParticles + if particle <: FermionStateful + mom = initialMomenta[1] + push!(particles, particle(mom, in_spin())) + elseif particle <: PhotonStateful + for i in 1:n + mom = initialMomenta[i + 1] + push!(particles, particle(mom, in_pol())) + end + else + @assert false + end + end + + final_momenta = MetagraphOptimization.generate_physical_massive_moms( + RNG, ss, outputMasses + ) + index = 1 + for (particle, n) in processDescription.outParticles + for _ in 1:n + if particle <: FermionStateful + push!(particles, particle(final_momenta[index], out_spin())) + elseif particle <: PhotonStateful + push!(particles, particle(final_momenta[index], out_pol())) + end + index += 1 + end + end + + inFerms = MetagraphOptimization._svector_from_type( + processDescription, FermionStateful{Incoming,in_spin}, particles + ) + outFerms = MetagraphOptimization._svector_from_type( + processDescription, FermionStateful{Outgoing,out_spin}, particles + ) + inAntiferms = MetagraphOptimization._svector_from_type( + processDescription, AntiFermionStateful{Incoming,in_spin}, particles + ) + outAntiferms = MetagraphOptimization._svector_from_type( + processDescription, AntiFermionStateful{Outgoing,out_spin}, particles + ) + inPhotons = MetagraphOptimization._svector_from_type( + processDescription, PhotonStateful{Incoming,in_pol}, particles + ) + outPhotons = MetagraphOptimization._svector_from_type( + processDescription, PhotonStateful{Outgoing,out_pol}, particles + ) + + processInput = QEDProcessInput( + processDescription, + inFerms, + outFerms, + inAntiferms, + outAntiferms, + inPhotons, + outPhotons, + ) + + push!(result, processInput) + end + + return result +end diff --git a/src/models/physics_models/qed/create.jl b/src/models/physics_models/qed/create.jl index b8a8a50..94454d8 100644 --- a/src/models/physics_models/qed/create.jl +++ b/src/models/physics_models/qed/create.jl @@ -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 diff --git a/src/models/physics_models/qed/diagrams.jl b/src/models/physics_models/qed/diagrams.jl index e3f94a6..27082e5 100644 --- a/src/models/physics_models/qed/diagrams.jl +++ b/src/models/physics_models/qed/diagrams.jl @@ -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}) @@ -509,7 +510,9 @@ end Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram for the given order and n input and m output photons. """ -function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int) +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], @@ -587,7 +590,6 @@ function gen_compton_diagrams(n::Int, m::Int) return vcat(diagrams...) end - """ gen_compton_diagrams_one_side(n::Int, m::Int) diff --git a/test/unit_tests_qedmodel.jl b/test/unit_tests_qedmodel.jl index 6a9defe..9a48050 100644 --- a/test/unit_tests_qedmodel.jl +++ b/test/unit_tests_qedmodel.jl @@ -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(