From a1581182ca91bc2cd7128259be6a57c637ebaf2a Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Tue, 2 Jul 2024 10:50:30 +0200 Subject: [PATCH 1/9] WIP --- Project.toml | 2 +- src/MetagraphOptimization.jl | 32 ++--- src/models/physics_models/interface.jl | 7 +- src/models/physics_models/qed/particle.jl | 145 ++++++---------------- src/models/print.jl | 10 -- 5 files changed, 62 insertions(+), 134 deletions(-) diff --git a/Project.toml b/Project.toml index 63f9d2e..ae82e04 100644 --- a/Project.toml +++ b/Project.toml @@ -14,12 +14,12 @@ 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" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" [extras] CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" diff --git a/src/MetagraphOptimization.jl b/src/MetagraphOptimization.jl index 608b515..429b7d4 100644 --- a/src/MetagraphOptimization.jl +++ b/src/MetagraphOptimization.jl @@ -173,22 +173,24 @@ 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") diff --git a/src/models/physics_models/interface.jl b/src/models/physics_models/interface.jl index 32478cf..17b90fd 100644 --- a/src/models/physics_models/interface.jl +++ b/src/models/physics_models/interface.jl @@ -1,4 +1,3 @@ -import QEDbase.mass import QEDbase.AbstractParticle """ @@ -9,13 +8,13 @@ Base type for a model, e.g. ABC-Model or QED. This is used to dispatch many func abstract type AbstractPhysicsModel <: AbstractModel end """ - ParticleValue{ParticleType <: AbstractParticle} + ParticleValue{ParticleType <: AbstractParticleStateful} -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. +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 <: AbstractParticle, ValueType} +struct ParticleValue{ParticleType <: AbstractParticleStateful, ValueType} p::ParticleType v::ValueType end diff --git a/src/models/physics_models/qed/particle.jl b/src/models/physics_models/qed/particle.jl index d51fedd..130dc2e 100644 --- a/src/models/physics_models/qed/particle.jl +++ b/src/models/physics_models/qed/particle.jl @@ -1,8 +1,6 @@ using QEDprocesses using StaticArrays -import QEDbase.mass -# TODO check const e = sqrt(4π / 137) """ @@ -12,17 +10,6 @@ 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 @@ -35,7 +22,7 @@ struct QEDProcessDescription <: AbstractProcessDescription outParticles::Dict{Type{<:QEDParticle{Outgoing}}, Int} end -QEDParticleValue{ParticleType <: QEDParticle} = Union{ +QEDParticleValue{ParticleType <: ParticleStateful} = Union{ ParticleValue{ParticleType, BiSpinor}, ParticleValue{ParticleType, AdjointBiSpinor}, ParticleValue{ParticleType, DiracMatrix}, @@ -43,92 +30,48 @@ QEDParticleValue{ParticleType <: QEDParticle} = Union{ 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} +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}}, ::Type{ParticleStateful{Outgoing, Electron}}) = + ParticleStateful{Incoming, Photon} +interaction_result(::Type{ParticleStateful{Incoming, Electron}}, ::Type{ParticleStateful{Incoming, Positron}}) = + ParticleStateful{Incoming, Photon} 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} + ::Type{ParticleStateful{Incoming, Electron}}, + ::Type{ParticleStateful{<:ParticleDirection, Photon}}, +) = ParticleStateful{Outgoing, Electron} +interaction_result(::Type{ParticleStateful{Outgoing, Electron}}, ::Type{ParticleStateful{Incoming, Electron}}) = + ParticleStateful{Incoming, Photon} +interaction_result(::Type{ParticleStateful{Outgoing, Electron}}, ::Type{ParticleStateful{Outgoing, Positron}}) = + ParticleStateful{Incoming, Photon} 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} + ::Type{ParticleStateful{Outgoing, Electron}}, + ::Type{ParticleStateful{<:ParticleDirection, Photon}}, +) = ParticleStateful{Incoming, Electron} # 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) +interaction_result(::Type{ParticleStateful{Incoming, Positron}}, t2::Type{<:ParticleStateful}) = + interaction_result(ParticleStateful{Outgoing, Electron}, t2) +interaction_result(::Type{ParticleStateful{Outgoing, Positron}}, t2::Type{<:ParticleStateful}) = + interaction_result(ParticleStateful{Incoming, Electron}, t2) # photon commutativity -interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:QEDParticle}) = interaction_result(t2, t1) +interaction_result(t1::Type{<:ParticleStateful{<:ParticleDirection, Photon}}, t2::Type{<:ParticleStateful}) = + interaction_result(t2, t1) # but prevent stack overflow -function interaction_result(t1::Type{<:PhotonStateful}, t2::Type{<:PhotonStateful}) +function interaction_result( + t1::Type{ParticleStateful{<:ParticleDirection, Photon}}, + t2::Type{ParticleStateful{<:ParticleDirection, Photon}}, +) @assert false "Invalid interaction between particles of types $t1 and $t2" end @@ -137,18 +80,12 @@ end 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} +propagation_result(::Type{ParticleStateful{Incoming, Electron}}) = ParticleStateful{Outgoing, Electron} +propagation_result(::Type{ParticleStateful{Outgoing, Electron}}) = ParticleStateful{Incoming, Electron} +propagation_result(::Type{ParticleStateful{Incoming, Positron}}) = ParticleStateful{Outgoing, Positron} +propagation_result(::Type{ParticleStateful{Outgoing, Positron}}) = ParticleStateful{Incoming, Positron} +propagation_result(::Type{ParticleStateful{Incoming, Photon}}) = ParticleStateful{Outgoing, Photon} +propagation_result(::Type{ParticleStateful{Outgoing, Photon}}) = ParticleStateful{Incoming, Photon} """ types(::QEDModel) @@ -157,12 +94,12 @@ 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}, + ParticleStateful{Incoming, Photon}, + ParticleStateful{Outgoing, Photon}, + ParticleStateful{Incoming, Electron}, + ParticleStateful{Outgoing, Electron}, + ParticleStateful{Incoming, Positron}, + ParticleStateful{Outgoing, Positron}, ] end @@ -235,9 +172,9 @@ end @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 QEDbase.mass(::Type{<:FermionStateful}) = 1.0 +@inline QEDbase.mass(::Type{<:AntiFermionStateful}) = 1.0 +@inline QEDbase.mass(::Type{<:PhotonStateful}) = 0.0 @inline invert_momentum(p::FermionStateful{Dir, Spin}) where {Dir, Spin} = FermionStateful{Dir, Spin}(-p.momentum, p.spin) diff --git a/src/models/print.jl b/src/models/print.jl index 00b8489..e69de29 100644 --- 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 -- 2.49.0 From 6a9a7b41f18ab0df4e7a10b7d746263420fc7af1 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 3 Jul 2024 20:24:53 +0200 Subject: [PATCH 2/9] rework a lot of the QED model to use QEDcore/base/processes --- Project.toml | 8 +- notebooks/abc_model_large.ipynb | 4 +- src/MetagraphOptimization.jl | 10 +- src/QEDprocesses_patch.jl | 71 +++++ src/code_gen/function.jl | 24 +- src/code_gen/tape_machine.jl | 37 +-- src/code_gen/type.jl | 4 +- src/models/interface.jl | 4 +- src/models/physics_models/interface.jl | 1 + src/models/physics_models/qed/compute.jl | 43 +-- src/models/physics_models/qed/create.jl | 66 ++-- src/models/physics_models/qed/diagrams.jl | 111 ++++--- src/models/physics_models/qed/parse.jl | 26 +- src/models/physics_models/qed/particle.jl | 348 ++++++++++------------ src/models/physics_models/qed/print.jl | 16 +- src/task/properties.jl | 13 +- src/utility.jl | 4 +- test/Project.toml | 10 - test/runtests.jl | 4 + test/unit_tests_qed_diagrams.jl | 2 +- test/unit_tests_qedmodel.jl | 88 +++--- 21 files changed, 458 insertions(+), 436 deletions(-) create mode 100644 src/QEDprocesses_patch.jl delete mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index ae82e04..73ce957 100644 --- a/Project.toml +++ b/Project.toml @@ -19,11 +19,17 @@ 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" [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"] diff --git a/notebooks/abc_model_large.ipynb b/notebooks/abc_model_large.ipynb index e98871c..ef75af5 100644 --- 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, diff --git a/src/MetagraphOptimization.jl b/src/MetagraphOptimization.jl index 429b7d4..77a85eb 100644 --- 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 @@ -114,6 +115,7 @@ import Base.delete! import Base.insert! import Base.collect +include("QEDprocesses_patch.jl") include("devices/interface.jl") include("task/type.jl") @@ -175,6 +177,7 @@ include("models/print.jl") include("models/physics_models/interface.jl") +#= include("models/physics_models/abc/types.jl") include("models/physics_models/abc/particle.jl") include("models/physics_models/abc/compute.jl") @@ -182,6 +185,7 @@ 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") @@ -199,7 +203,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") diff --git a/src/QEDprocesses_patch.jl b/src/QEDprocesses_patch.jl new file mode 100644 index 0000000..1eba98c --- /dev/null +++ b/src/QEDprocesses_patch.jl @@ -0,0 +1,71 @@ +# 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}, +) 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 + +@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 diff --git a/src/code_gen/function.jl b/src/code_gen/function.jl index 588b350..fc9c2ba 100644 --- a/src/code_gen/function.jl +++ b/src/code_gen/function.jl @@ -1,10 +1,10 @@ """ - get_compute_function(graph::DAG, process::AbstractProcessDescription, machine::Machine) + get_compute_function(graph::DAG, instance, machine::Machine) -Return a function of signature `compute_(input::AbstractProcessInput)`, which will return the result of the DAG computation on the given input. +Return a function of signature `compute_(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)...) @@ -13,7 +13,7 @@ function get_compute_function(graph::DAG, process::AbstractProcessDescription, m 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,12 +22,12 @@ 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_(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)...) @@ -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 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 diff --git a/src/code_gen/tape_machine.jl b/src/code_gen/tape_machine.jl index e986ad4..db8ade7 100644 --- a/src/code_gen/tape_machine.jl +++ b/src/code_gen/tape_machine.jl @@ -5,7 +5,7 @@ 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{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 @@ -15,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 @@ -32,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 @@ -48,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 @@ -86,27 +82,19 @@ Return a `Vector{Expr}` doing the input assignments from the given `problemInput """ function gen_input_assignment_code( inputSymbols::Dict{String, Vector{Symbol}}, - instance::AbstractProblemInstance, + instance, machine::Machine, problemInputSymbol::Symbol = :input, ) assignInputs = Vector{FunctionCall}() for (name, symbols) in inputSymbols - (type, index) = type_index_from_name(model(instance), 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{2, Any}(type, index), - SVector{1, Symbol}(problemInputSymbol), - symbol, - device, - ), + FunctionCall(input_expr, SVector{1, Any}(name), SVector{1, Symbol}(problemInputSymbol), symbol, device), ) end end @@ -121,12 +109,7 @@ 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, - instance::AbstractProblemInstance, - machine::Machine, - scheduler::AbstractScheduler = GreedyScheduler(), -) +function gen_tape(graph::DAG, instance, machine::Machine, scheduler::AbstractScheduler = GreedyScheduler()) schedule = schedule_dag(scheduler, graph, machine) # get inSymbols @@ -145,7 +128,7 @@ function gen_tape( initCaches = gen_cache_init_code(machine) assignInputs = gen_input_assignment_code(inputSyms, instance, machine, :input) - return Tape(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine) + return Tape{input_type(instance)}(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine) end """ diff --git a/src/code_gen/type.jl b/src/code_gen/type.jl index 46e9fd1..db420fd 100644 --- a/src/code_gen/type.jl +++ b/src/code_gen/type.jl @@ -3,6 +3,8 @@ 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 @@ -14,6 +16,6 @@ struct Tape{INPUT} inputSymbols::Dict{String, Vector{Symbol}} outputSymbol::Symbol cache::Dict{Symbol, Any} - instance::AbstractProblemInstance + instance::Any machine::Machine end diff --git a/src/models/interface.jl b/src/models/interface.jl index 3ed1b34..8bf7514 100644 --- a/src/models/interface.jl +++ b/src/models/interface.jl @@ -39,8 +39,8 @@ Generate the [`DAG`](@ref) for the given [`AbstractProblemInstance`](@ref). Ever function graph end """ - input_expr(::AbstractProblemInstance, input_sym::Symbol, name::String) + input_expr(::AbstractProblemInstance, name::String, input) -For a given [`AbstractProblemInstance`](@ref), a `Symbol` on which the problem input is available (see [`input_type`](@ref)), and entry node name, return an `Expr` getting that specific input value from the +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. """ function input_expr end diff --git a/src/models/physics_models/interface.jl b/src/models/physics_models/interface.jl index 17b90fd..5654735 100644 --- a/src/models/physics_models/interface.jl +++ b/src/models/physics_models/interface.jl @@ -29,6 +29,7 @@ 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 diff --git a/src/models/physics_models/qed/compute.jl b/src/models/physics_models/qed/compute.jl index b30994f..ce94d16 100644 --- a/src/models/physics_models/qed/compute.jl +++ b/src/models/physics_models/qed/compute.jl @@ -1,13 +1,9 @@ using StaticArrays -""" - compute(::ComputeTaskQED_P, data::QEDParticleValue) +function input_expr(name::String, psp::PhaseSpacePoint) + (type, index) = type_index_from_name(QEDModel(), name) -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)) + return ParticleValue(type(momentum(psp, particle_direction(type), particle_species(type), index)), 0.0im) end """ @@ -15,9 +11,9 @@ end 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::ParticleValue{P, V}) where {P <: ParticleStateful, V <: ValueType} 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), spin_or_pol(part)) return ParticleValue{P, typeof(state)}( data.p, state, # will return a SLorentzVector{ComplexF64}, BiSpinor or AdjointBiSpinor @@ -31,9 +27,9 @@ Compute a vertex. Preserve momentum and particle types (e + gamma->p etc.) to cr """ 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() @@ -63,9 +59,16 @@ 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{ParticleStateful{D1, S1}, V1}, + data2::ParticleValue{ParticleStateful{D2, S2}, V2}, +) where { + D1 <: ParticleDirection, + D2 <: ParticleDirection, + S1 <: Union{Electron, Positron}, + S2 <: Union{Electron, Positron}, + V1 <: ValueType, + V2 <: ValueType, +} #@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 +83,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 +93,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 diff --git a/src/models/physics_models/qed/create.jl b/src/models/physics_models/qed/create.jl index 85ee3aa..b8a8a50 100644 --- a/src/models/physics_models/qed/create.jl +++ b/src/models/physics_models/qed/create.jl @@ -1,7 +1,7 @@ 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 @@ -12,72 +12,48 @@ function _svector_from_type(processDescription::QEDProcessDescription, type, par 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) diff --git a/src/models/physics_models/qed/diagrams.jl b/src/models/physics_models/qed/diagrams.jl index 9477eb6..e3f94a6 100644 --- 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) @@ -166,11 +156,11 @@ copy(fd::FeynmanDiagram) = FeynmanDiagram(deepcopy(fd.vertices), copy(fd.tie[]), deepcopy(fd.particles), copy(fd.type_ids)) """ - 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 +429,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 +451,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 +460,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 +475,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 +488,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 +503,7 @@ 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}, i) for i in 1:n], + [FeynmanParticle(ParticleStateful{Outgoing, Photon}, 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} => 1, + ParticleStateful{Outgoing, Electron} => 1, + ParticleStateful{Incoming, Photon} => n, + ParticleStateful{Outgoing, Photon} => 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}, iterations), photons[order[left_index]], - FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1), + FeynmanParticle(ParticleStateful{Incoming, Electron}, 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}, iterations), photons[order[right_index]], - FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1), + FeynmanParticle(ParticleStateful{Outgoing, Electron}, iterations + 1), ) right_index -= 1 add_vertex!(new_diagram, v_right) @@ -575,7 +566,7 @@ 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) @@ -583,8 +574,8 @@ end 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)])...] @@ -603,8 +594,8 @@ end 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 +614,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 +661,5 @@ function gen_diagrams(fd::FeynmanDiagram) end return remove_duplicates(results) + =# end diff --git a/src/models/physics_models/qed/parse.jl b/src/models/physics_models/qed/parse.jl index ffcac74..31f2e85 100644 --- a/src/models/physics_models/qed/parse.jl +++ b/src/models/physics_models/qed/parse.jl @@ -5,8 +5,16 @@ 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}() + 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\"") @@ -19,14 +27,14 @@ function parse_process(str::AbstractString, model::QEDModel) end for t in types(model) - if (isincoming(t)) + if (is_incoming(t)) inCount = count(x -> x == String(t)[1], inStr) if inCount != 0 inParticles[t] = inCount end end - if (isoutgoing(t)) + if (is_outgoing(t)) outCount = count(x -> x == String(t)[1], outStr) if outCount != 0 outParticles[t] = outCount @@ -40,5 +48,13 @@ function parse_process(str::AbstractString, model::QEDModel) throw("Encountered unknown characters in the output part of process \"$str\"") end - return QEDProcessDescription(inParticles, outParticles) + + return GenericQEDProcess( + inParticles[ParticleStateful{Incoming, Photon, SFourMomentum}], + outParticles[ParticleStateful{Outgoing, Photon, SFourMomentum}], + inParticles[ParticleStateful{Incoming, Electron, SFourMomentum}], + outParticles[ParticleStateful{Outgoing, Electron, SFourMomentum}], + inParticles[ParticleStateful{Incoming, Positron, SFourMomentum}], + outParticles[ParticleStateful{Outgoing, Positron, SFourMomentum}], + ) end diff --git a/src/models/physics_models/qed/particle.jl b/src/models/physics_models/qed/particle.jl index 130dc2e..f41fa8b 100644 --- a/src/models/physics_models/qed/particle.jl +++ b/src/models/physics_models/qed/particle.jl @@ -10,25 +10,70 @@ Singleton definition for identification of the QED-Model. """ struct QEDModel <: AbstractPhysicsModel end -""" - QEDProcessDescription <: AbstractProcessDescription +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 -A description of a process in the QED-Model. Contains the input and output particles. +QEDbase.particle_direction(::Type{<:ParticleStateful{DIR}}) where {DIR <: ParticleDirection} = DIR() +QEDbase.particle_species( + ::Type{<:ParticleStateful{DIR, SPECIES}}, +) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType} = SPECIES() -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} +_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} <: AbstractProcessDefinition where {INT <: Tuple, OUTT <: Tuple} + incoming_particles::INT + outgoing_particles::OUTT + + function GenericQEDProcess(in_particles::INT, out_particles::OUTT) where {INT <: Tuple, OUTT <: Tuple} + _assert_particle_type_tuple(in_particles) + _assert_particle_type_tuple(out_particles) + + return new{INT, OUTT}(in_particles, out_particles) + end + + """ + GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int) + + Convenience constructor from numbers of input/output photons, electrons and positrons. + """ + function GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int) + 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) + end end -QEDParticleValue{ParticleType <: ParticleStateful} = Union{ - ParticleValue{ParticleType, BiSpinor}, - ParticleValue{ParticleType, AdjointBiSpinor}, - ParticleValue{ParticleType, DiracMatrix}, - ParticleValue{ParticleType, SLorentzVector{Float64}}, - ParticleValue{ParticleType, ComplexF64}, -} +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} @@ -39,39 +84,53 @@ function interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: ParticleSta @assert false "Invalid interaction between particles of types $t1 and $t2" end -interaction_result(::Type{ParticleStateful{Incoming, Electron}}, ::Type{ParticleStateful{Outgoing, Electron}}) = - ParticleStateful{Incoming, Photon} -interaction_result(::Type{ParticleStateful{Incoming, Electron}}, ::Type{ParticleStateful{Incoming, Positron}}) = - ParticleStateful{Incoming, Photon} interaction_result( - ::Type{ParticleStateful{Incoming, Electron}}, - ::Type{ParticleStateful{<:ParticleDirection, Photon}}, -) = ParticleStateful{Outgoing, Electron} + ::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}}, ::Type{ParticleStateful{Incoming, Electron}}) = - ParticleStateful{Incoming, Photon} -interaction_result(::Type{ParticleStateful{Outgoing, Electron}}, ::Type{ParticleStateful{Outgoing, Positron}}) = - ParticleStateful{Incoming, Photon} interaction_result( - ::Type{ParticleStateful{Outgoing, Electron}}, - ::Type{ParticleStateful{<:ParticleDirection, Photon}}, -) = ParticleStateful{Incoming, Electron} + ::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}}, t2::Type{<:ParticleStateful}) = - interaction_result(ParticleStateful{Outgoing, Electron}, t2) -interaction_result(::Type{ParticleStateful{Outgoing, Positron}}, t2::Type{<:ParticleStateful}) = - interaction_result(ParticleStateful{Incoming, Electron}, t2) +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{<:ParticleDirection, Photon}}, t2::Type{<:ParticleStateful}) = - interaction_result(t2, t1) +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{<:ParticleDirection, Photon}}, - t2::Type{ParticleStateful{<:ParticleDirection, Photon}}, -) + 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 @@ -80,12 +139,18 @@ end Return the type of the inverted direction. E.g. """ -propagation_result(::Type{ParticleStateful{Incoming, Electron}}) = ParticleStateful{Outgoing, Electron} -propagation_result(::Type{ParticleStateful{Outgoing, Electron}}) = ParticleStateful{Incoming, Electron} -propagation_result(::Type{ParticleStateful{Incoming, Positron}}) = ParticleStateful{Outgoing, Positron} -propagation_result(::Type{ParticleStateful{Outgoing, Positron}}) = ParticleStateful{Incoming, Positron} -propagation_result(::Type{ParticleStateful{Incoming, Photon}}) = ParticleStateful{Outgoing, Photon} -propagation_result(::Type{ParticleStateful{Outgoing, Photon}}) = ParticleStateful{Incoming, Photon} +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) @@ -94,12 +159,12 @@ Return a Vector of the possible types of particle in the [`QEDModel`](@ref). """ function types(::QEDModel) return [ - ParticleStateful{Incoming, Photon}, - ParticleStateful{Outgoing, Photon}, - ParticleStateful{Incoming, Electron}, - ParticleStateful{Outgoing, Electron}, - ParticleStateful{Incoming, Positron}, - ParticleStateful{Outgoing, Positron}, + 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 @@ -116,94 +181,39 @@ String(::Type{SpinDown}) = "spindown" String(::Incoming) = "i" String(::Outgoing) = "o" -function String(::Type{<:PhotonStateful}) +function String(::Type{<:ParticleStateful{DIR, Photon}}) where {DIR <: ParticleDirection} return "k" end -function String(::Type{<:FermionStateful}) +function String(::Type{<:ParticleStateful{DIR, Electron}}) where {DIR <: ParticleDirection} return "e" end -function String(::Type{<:AntiFermionStateful}) +function String(::Type{<:ParticleStateful{DIR, Positron}}) where {DIR <: ParticleDirection} 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 QEDbase.mass(::Type{<:FermionStateful}) = 1.0 -@inline QEDbase.mass(::Type{<:AntiFermionStateful}) = 1.0 -@inline QEDbase.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}) + caninteract(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful}) -For two given [`QEDParticle`](@ref) types, return whether they can interact at a vertex. This is equivalent to `!issame(T1, T2)`. +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{<:QEDParticle}, T2::Type{<:QEDParticle}) +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 (T1 <: PhotonStateful && T2 <: PhotonStateful) + if (S1 == Photon && S2 == Photon) return false end for (P1, P2) in [(T1, T2), (T2, T1)] - if (P1 <: FermionStateful{Incoming} && P2 <: AntiFermionStateful{Outgoing}) + if (P1 <: ParticleStateful{Incoming, Electron} && P2 <: ParticleStateful{Outgoing, Positron}) return false end - if (P1 <: FermionStateful{Outgoing} && P2 <: AntiFermionStateful{Incoming}) + if (P1 <: ParticleStateful{Outgoing, Electron} && P2 <: ParticleStateful{Incoming, Positron}) return false end end @@ -213,30 +223,30 @@ end function type_index_from_name(::QEDModel, name::String) if startswith(name, "ki") - return (PhotonStateful{Incoming, PolX}, parse(Int, name[3:end])) + return (ParticleStateful{Incoming, Photon, SFourMomentum}, parse(Int, name[3:end])) elseif startswith(name, "ko") - return (PhotonStateful{Outgoing, PolX}, parse(Int, name[3:end])) + return (ParticleStateful{Outgoing, Photon, SFourMomentum}, parse(Int, name[3:end])) elseif startswith(name, "ei") - return (FermionStateful{Incoming, SpinUp}, parse(Int, name[3:end])) + return (ParticleStateful{Incoming, Electron, SFourMomentum}, parse(Int, name[3:end])) elseif startswith(name, "eo") - return (FermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end])) + return (ParticleStateful{Outgoing, Electron, SFourMomentum}, parse(Int, name[3:end])) elseif startswith(name, "pi") - return (AntiFermionStateful{Incoming, SpinUp}, parse(Int, name[3:end])) + return (ParticleStateful{Incoming, Positron, SFourMomentum}, parse(Int, name[3:end])) elseif startswith(name, "po") - return (AntiFermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end])) + 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{<:QEDParticle}, T2::Type{<:QEDParticle}) + issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful}) -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)`. +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{<:QEDParticle}, T2::Type{<:QEDParticle}) +function issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful}) return !caninteract(T1, T2) end @@ -250,12 +260,12 @@ Return the factor of a vertex in a QED feynman diagram. return -1im * e * gamma() end -@inline function QED_inner_edge(p::QEDParticle) - return propagator(particle(p), p.momentum) +@inline function QED_inner_edge(p::ParticleStateful) + return propagator(particle_species(p), momentum(p)) end """ - QED_conserve_momentum(p1::QEDParticle, p2::QEDParticle) + QED_conserve_momentum(p1::ParticleStateful, p2::ParticleStateful) Calculate and return a new particle from two given interacting ones at a vertex. """ @@ -265,43 +275,26 @@ function QED_conserve_momentum( ) 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}}, + Species1 <: AbstractParticleType, + Species2 <: AbstractParticleType, + P1 <: ParticleStateful{Dir1, Species1}, + P2 <: ParticleStateful{Dir2, Species2}, } P3 = interaction_result(P1, P2) - p1_mom = p1.momentum + p1_mom = momentum(p1) if (Dir1 <: Outgoing) p1_mom *= -1 end - p2_mom = p2.momentum + p2_mom = momentum(p2) if (Dir2 <: Outgoing) p2_mom *= -1 end p3_mom = p1_mom + p2_mom - if (typeof(direction(P3)) <: Incoming) - return P3(-p3_mom) + if (particle_direction(P3) isa Incoming) + return ParticleStateful(particle_direction(P3), particle_species(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}} + return ParticleStateful(particle_direction(P3), particle_species(P3), p3_mom) end """ @@ -309,37 +302,22 @@ end Return the model of this process description. """ -model(::QEDProcessDescription) = QEDModel() -model(::QEDProcessInput) = QEDModel() +model(::GenericQEDProcess) = QEDModel() +model(::PhaseSpacePoint) = 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] +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 diff --git a/src/models/physics_models/qed/print.jl b/src/models/physics_models/qed/print.jl index c5cb3f5..2ac8ac4 100644 --- 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) diff --git a/src/task/properties.jl b/src/task/properties.jl index c608843..2c4b320 100644 --- 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) diff --git a/src/utility.jl b/src/utility.jl index ac02bbf..a610663 100644 --- 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) diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 41c02e3..0000000 --- a/test/Project.toml +++ /dev/null @@ -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" diff --git a/test/runtests.jl b/test/runtests.jl index fa37b10..277f5de 100644 --- 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 +=# \ No newline at end of file diff --git a/test/unit_tests_qed_diagrams.jl b/test/unit_tests_qed_diagrams.jl index 802fb50..cd3486b 100644 --- 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) diff --git a/test/unit_tests_qedmodel.jl b/test/unit_tests_qedmodel.jl index 8126605..6a9defe 100644 --- 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)) @@ -88,8 +87,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 +102,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 +116,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 +145,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 +175,7 @@ end end end +#= @testset "Compton" begin import MetagraphOptimization.insert_node! import MetagraphOptimization.insert_edge! @@ -347,3 +338,4 @@ end @test isapprox(compute_function.(input), reduced_compute_function.(input)) end end +=# \ No newline at end of file -- 2.49.0 From b5d92b729c07b958dffb26a163aa118e3f416dcb Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 4 Jul 2024 15:31:22 +0200 Subject: [PATCH 3/9] Get compute function working --- src/models/physics_models/interface.jl | 11 ++++ src/models/physics_models/qed/compute.jl | 26 ++++++--- src/models/physics_models/qed/parse.jl | 27 ++++++---- src/models/physics_models/qed/particle.jl | 66 ++++++++++++++++++++--- 4 files changed, 107 insertions(+), 23 deletions(-) diff --git a/src/models/physics_models/interface.jl b/src/models/physics_models/interface.jl index 5654735..a345754 100644 --- a/src/models/physics_models/interface.jl +++ b/src/models/physics_models/interface.jl @@ -19,6 +19,17 @@ struct ParticleValue{ParticleType <: AbstractParticleStateful, ValueType} 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 diff --git a/src/models/physics_models/qed/compute.jl b/src/models/physics_models/qed/compute.jl index ce94d16..239c805 100644 --- a/src/models/physics_models/qed/compute.jl +++ b/src/models/physics_models/qed/compute.jl @@ -3,17 +3,24 @@ using StaticArrays function input_expr(name::String, psp::PhaseSpacePoint) (type, index) = type_index_from_name(QEDModel(), name) - return ParticleValue(type(momentum(psp, particle_direction(type), particle_species(type), index)), 0.0im) + return ParticleValueSP( + type(momentum(psp, particle_direction(type), particle_species(type), index)), + 0.0im, + spin_or_pol(process(psp), 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::ParticleValue{P, V}) where {P <: ParticleStateful, V <: ValueType} +function compute( + ::ComputeTaskQED_U, + data::ParticleValueSP{P, SP, V}, +) where {P <: ParticleStateful, V <: ValueType, SP <: AbstractSpinOrPolarization} part::P = data.p - state = base_state(particle_species(part), particle_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 @@ -21,7 +28,7 @@ function compute(::ComputeTaskQED_U, data::ParticleValue{P, V}) where {P <: Part 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. """ @@ -49,7 +56,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). @@ -59,8 +66,8 @@ For valid inputs, both input particles should have the same momenta at this poin """ function compute( ::ComputeTaskQED_S2, - data1::ParticleValue{ParticleStateful{D1, S1}, V1}, - data2::ParticleValue{ParticleStateful{D2, S2}, V2}, + data1::ParticleValue{P1, V1}, + data2::ParticleValue{P2, V2}, ) where { D1 <: ParticleDirection, D2 <: ParticleDirection, @@ -68,6 +75,9 @@ function compute( 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)" diff --git a/src/models/physics_models/qed/parse.jl b/src/models/physics_models/qed/parse.jl index 31f2e85..7b8f6e0 100644 --- a/src/models/physics_models/qed/parse.jl +++ b/src/models/physics_models/qed/parse.jl @@ -1,10 +1,18 @@ +# 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) +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, @@ -48,13 +56,14 @@ function parse_process(str::AbstractString, model::QEDModel) 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( - inParticles[ParticleStateful{Incoming, Photon, SFourMomentum}], - outParticles[ParticleStateful{Outgoing, Photon, SFourMomentum}], - inParticles[ParticleStateful{Incoming, Electron, SFourMomentum}], - outParticles[ParticleStateful{Outgoing, Electron, SFourMomentum}], - inParticles[ParticleStateful{Incoming, Positron, SFourMomentum}], - outParticles[ParticleStateful{Outgoing, Positron, SFourMomentum}], - ) + return GenericQEDProcess(in_ph, out_ph, in_el, out_el, in_pos, out_pos, in_spin_pols, out_spin_pols) end diff --git a/src/models/physics_models/qed/particle.jl b/src/models/physics_models/qed/particle.jl index f41fa8b..a21021e 100644 --- a/src/models/physics_models/qed/particle.jl +++ b/src/models/physics_models/qed/particle.jl @@ -25,29 +25,83 @@ _assert_particle_type_tuple(t::Tuple{AbstractParticleType, Vararg}) = _assert_pa _assert_particle_type_tuple(t::Any) = throw(InvalidInputError("invalid input, provide a tuple of AbstractParticleTypes to construct a GenericQEDProcess")) -struct GenericQEDProcess{INT, OUTT} <: AbstractProcessDefinition where {INT <: Tuple, OUTT <: Tuple} +struct GenericQEDProcess{INT, OUTT, INSP, OUTSP} <: + AbstractProcessDefinition where {INT <: Tuple, OUTT <: Tuple, INSP <: Tuple, OUTSP <: Tuple} incoming_particles::INT outgoing_particles::OUTT - function GenericQEDProcess(in_particles::INT, out_particles::OUTT) where {INT <: Tuple, OUTT <: Tuple} + 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) - return new{INT, OUTT}(in_particles, 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(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int) + 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) + 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) + 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 -- 2.49.0 From 55501c15c8c92f73dc90f7f8e98cf4dc19382e9d Mon Sep 17 00:00:00 2001 From: Rubydragon Date: Thu, 27 Jun 2024 13:21:53 +0200 Subject: [PATCH 4/9] WIP --- examples/congruent_in_ph.jl | 123 +++++++++++++++++++ src/models/physics_models/qed/create.jl | 47 ++++---- src/models/physics_models/qed/diagrams.jl | 10 +- test/unit_tests_qedmodel.jl | 140 +++++++++++----------- 4 files changed, 223 insertions(+), 97 deletions(-) create mode 100644 examples/congruent_in_ph.jl 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( -- 2.49.0 From 92f534f6bf061ca7024acc4bfb703b64aa94d988 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 4 Jul 2024 17:03:18 +0200 Subject: [PATCH 5/9] Fix congruent ph script --- examples/congruent_in_ph.jl | 156 ++++++++++++++++-------------------- 1 file changed, 67 insertions(+), 89 deletions(-) diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl index a0639bf..1c91757 100644 --- a/examples/congruent_in_ph.jl +++ b/examples/congruent_in_ph.jl @@ -1,8 +1,10 @@ using MetagraphOptimization using QEDbase using QEDcore +using QEDprocesses using Random using UUIDs +using CSV RNG = Random.MersenneTwister(123) @@ -21,103 +23,79 @@ function mock_machine() ) end -function congruent_input(processDescription::QEDProcessDescription, omega::Number) +function congruent_input_momenta(processDescription::GenericQEDProcess, 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 + 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 - initialMomenta = [ - i == 1 ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for - i in 1:length(inputMasses) - ] + initial_momenta = + [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)) + ss = sqrt(sum(initial_momenta) * sum(initial_momenta)) + final_momenta = MetagraphOptimization.generate_physical_massive_moms(RNG, ss, outputMasses) - 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 + 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 + +n = 1000000 +photons = 4 +omega = 2e-3 + +# n is the number of incoming photons +# omega is the number + +println("Generating $n inputs for $photons photons, omega=$omega...") + +# temp process to generate momenta +temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) +input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n] +results = [0.0im for _ in 1:n] + +for (in_pol, in_spin, out_pol, out_spin) in + Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()]) + + print("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 = build_psp.(Ref(process), input_momenta) + + print("Preparing graph... ") + # prepare function + graph = gen_graph(process) + optimize_to_fixpoint!(ReductionOptimizer(), graph) + print("Preparing function... ") + func = get_compute_function(graph, process, mock_machine()) + + print("Calculating... ") + for i in 1:n + results[i] += func(inputs[i])^2 + end + println("Done.") +end + +println("Writing results") +open("$(photons)_congruent_photons_omega_$(omega).csv", "w") do f + println(f, "out_photon_momentum;out_electron_momentum;result") + for (momentum, result) in Iterators.zip(input_momenta, results) + println(f, "$(momentum[2][1]);$(momentum[2][2]);$(result)") + end end -- 2.49.0 From 813d40cd30b8623e39f13260adbab770f4c49853 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 4 Jul 2024 20:24:39 +0200 Subject: [PATCH 6/9] Working state --- examples/congruent_in_ph.jl | 176 ++++++++++++++++++----- src/code_gen/function.jl | 4 +- src/code_gen/tape_machine.jl | 14 +- src/code_gen/type.jl | 2 +- src/models/interface.jl | 4 +- src/models/physics_models/qed/compute.jl | 24 +++- 6 files changed, 176 insertions(+), 48 deletions(-) diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl index 1c91757..1f8367d 100644 --- a/examples/congruent_in_ph.jl +++ b/examples/congruent_in_ph.jl @@ -5,6 +5,7 @@ using QEDprocesses using Random using UUIDs using CSV +using JLD2 RNG = Random.MersenneTwister(123) @@ -37,16 +38,66 @@ function congruent_input_momenta(processDescription::GenericQEDProcess, omega::N push!(outputMasses, mass(particle)) end - initial_momenta = - [i == 1 ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for i in 1:length(inputMasses)] + initial_momenta = [ + i == length(inputMasses) ? 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(initial_momenta) * sum(initial_momenta)) final_momenta = MetagraphOptimization.generate_physical_massive_moms(RNG, ss, outputMasses) return (tuple(initial_momenta...), tuple(final_momenta...)) end +# cos_theta ∈ [-1, 1] 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 + massSum = 0 + inputMasses = Vector{Float64}() + for particle in incoming_particles(processDescription) + massSum += mass(particle) + push!(inputMasses, mass(particle)) + end + outputMasses = Vector{Float64}() + for particle in outgoing_particles(processDescription) + massSum += mass(particle) + 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, ParticleStateful{Incoming, Photon, SFourMomentum}) + + 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, @@ -57,45 +108,106 @@ function build_psp(processDescription::GenericQEDProcess, momenta) ) end +return 0 + +# scenario 2 +N = 1000 +M = 1000 + +thetas = collect(LinRange(0, 2π, N)) +phis = collect(LinRange(0, 2π, M)) + +for photons in [6] + # temp process to generate momenta + for omega in [2e-3, 2e-6] + println("Generating $(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 = [ + congruent_input_momenta_scenario_2(temp_process, omega, theta, phi) for + (theta, phi) in Iterators.product(thetas, phis) + ] + results = [0.0 for _ in 1:length(input_momenta)] + + 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 = build_psp.(Ref(process), input_momenta) + + print("Preparing graph... ") + # prepare function + graph = gen_graph(process) + optimize_to_fixpoint!(ReductionOptimizer(), graph) + print("Preparing function... ") + func = get_compute_function(graph, process, mock_machine()) + + print("Calculating... ") + Threads.@threads for i in 1:length(inputs) + results[i] += abs2(func(inputs[i])) + end + 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) + + @save "$(photons)_congruent_photons_omega_$(omega)_grid.jld2" out_ph_moms out_el_moms results + end +end + + n = 1000000 -photons = 4 -omega = 2e-3 # n is the number of incoming photons # omega is the number -println("Generating $n inputs for $photons photons, omega=$omega...") +for photons in [6] + # temp process to generate momenta + for omega in [2e-3, 2e-6] + println("Generating $n inputs for $photons photons...") + temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) -# temp process to generate momenta -temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) -input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n] -results = [0.0im for _ in 1:n] + input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n] + results = [0.0 for _ in 1:length(input_momenta)] -for (in_pol, in_spin, out_pol, out_spin) in - Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()]) + i = 1 + for (in_pol, in_spin, out_pol, out_spin) in + Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()]) - print("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 = build_psp.(Ref(process), input_momenta) + 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 = build_psp.(Ref(process), input_momenta) - print("Preparing graph... ") - # prepare function - graph = gen_graph(process) - optimize_to_fixpoint!(ReductionOptimizer(), graph) - print("Preparing function... ") - func = get_compute_function(graph, process, mock_machine()) + print("Preparing graph... ") + # prepare function + graph = gen_graph(process) + optimize_to_fixpoint!(ReductionOptimizer(), graph) + print("Preparing function... ") + func = get_compute_function(graph, process, mock_machine()) - print("Calculating... ") - for i in 1:n - results[i] += func(inputs[i])^2 - end - println("Done.") -end + print("Calculating... ") + Threads.@threads for i in 1:n + results[i] += abs2(func(inputs[i])) + end + println("Done.") + i += 1 + end -println("Writing results") -open("$(photons)_congruent_photons_omega_$(omega).csv", "w") do f - println(f, "out_photon_momentum;out_electron_momentum;result") - for (momentum, result) in Iterators.zip(input_momenta, results) - println(f, "$(momentum[2][1]);$(momentum[2][2]);$(result)") + println("Writing results") + + out_ph_moms = getindex.(getindex.(input_momenta, 2), 1) + out_el_moms = getindex.(getindex.(input_momenta, 2), 2) + + @save "$(photons)_congruent_photons_omega_$(omega).jld2" out_ph_moms out_el_moms results end end diff --git a/src/code_gen/function.jl b/src/code_gen/function.jl index fc9c2ba..50b3df6 100644 --- a/src/code_gen/function.jl +++ b/src/code_gen/function.jl @@ -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])) diff --git a/src/code_gen/tape_machine.jl b/src/code_gen/tape_machine.jl index db8ade7..010fd6f 100644 --- a/src/code_gen/tape_machine.jl +++ b/src/code_gen/tape_machine.jl @@ -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 diff --git a/src/code_gen/type.jl b/src/code_gen/type.jl index db420fd..1b079dd 100644 --- a/src/code_gen/type.jl +++ b/src/code_gen/type.jl @@ -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 diff --git a/src/models/interface.jl b/src/models/interface.jl index 8bf7514..8ca05b2 100644 --- a/src/models/interface.jl +++ b/src/models/interface.jl @@ -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 diff --git a/src/models/physics_models/qed/compute.jl b/src/models/physics_models/qed/compute.jl index 239c805..01b1459 100644 --- a/src/models/physics_models/qed/compute.jl +++ b/src/models/physics_models/qed/compute.jl @@ -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))), $index)), + 0.0im, + $(construction_string(spin_or_pol(instance, type, index))), +)", ) end -- 2.49.0 From b784720859c969196fbbef5dff4b7cd03b305c55 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 4 Jul 2024 20:40:22 +0200 Subject: [PATCH 7/9] Minor fixes and commints --- examples/congruent_in_ph.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl index 1f8367d..24b430e 100644 --- a/examples/congruent_in_ph.jl +++ b/examples/congruent_in_ph.jl @@ -26,15 +26,12 @@ 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 - massSum = 0 inputMasses = Vector{Float64}() for particle in incoming_particles(processDescription) - massSum += mass(particle) push!(inputMasses, mass(particle)) end outputMasses = Vector{Float64}() for particle in outgoing_particles(processDescription) - massSum += mass(particle) push!(outputMasses, mass(particle)) end @@ -49,7 +46,7 @@ function congruent_input_momenta(processDescription::GenericQEDProcess, omega::N return (tuple(initial_momenta...), tuple(final_momenta...)) end -# cos_theta ∈ [-1, 1] and phi ∈ [0, 2π] +# theta ∈ [0, 2π] and phi ∈ [0, 2π] function congruent_input_momenta_scenario_2( processDescription::GenericQEDProcess, omega::Number, @@ -60,15 +57,12 @@ function congruent_input_momenta_scenario_2( # same as above # generate an input sample for given e + nk -> e' + k' process, where the nk are equal - massSum = 0 inputMasses = Vector{Float64}() for particle in incoming_particles(processDescription) - massSum += mass(particle) push!(inputMasses, mass(particle)) end outputMasses = Vector{Float64}() for particle in outgoing_particles(processDescription) - massSum += mass(particle) push!(outputMasses, mass(particle)) end @@ -163,7 +157,7 @@ for photons in [6] end end - +# scenario n = 1000000 # n is the number of incoming photons -- 2.49.0 From dee44dad66332603d1706f65263a20cfabd2af16 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Fri, 5 Jul 2024 02:15:03 +0200 Subject: [PATCH 8/9] Small fixes --- examples/congruent_in_ph.jl | 27 +++++++++++++++-------- src/models/physics_models/qed/diagrams.jl | 26 +++++++++------------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl index 24b430e..19acc24 100644 --- a/examples/congruent_in_ph.jl +++ b/examples/congruent_in_ph.jl @@ -1,3 +1,5 @@ +ENV["UCX_ERROR_SIGNALS"] = "SIGILL,SIGBUS,SIGFPE" + using MetagraphOptimization using QEDbase using QEDcore @@ -102,7 +104,8 @@ function build_psp(processDescription::GenericQEDProcess, momenta) ) end -return 0 +# hack to fix stacksize for threading +with_stacksize(f, n) = fetch(schedule(Task(f, n))) # scenario 2 N = 1000 @@ -111,7 +114,7 @@ M = 1000 thetas = collect(LinRange(0, 2π, N)) phis = collect(LinRange(0, 2π, M)) -for photons in [6] +for photons in 1:6 # temp process to generate momenta for omega in [2e-3, 2e-6] println("Generating $(N*M) inputs for $photons photons (Scenario 2 grid walk)...") @@ -121,7 +124,8 @@ for photons in [6] congruent_input_momenta_scenario_2(temp_process, omega, theta, phi) for (theta, phi) in Iterators.product(thetas, phis) ] - results = [0.0 for _ in 1:length(input_momenta)] + results = Array{Float64}(undef, size(input_momenta)) + fill!(results, 0.0) i = 1 for (in_pol, in_spin, out_pol, out_spin) in @@ -134,15 +138,17 @@ for photons in [6] inputs = build_psp.(Ref(process), input_momenta) print("Preparing graph... ") - # prepare function graph = gen_graph(process) optimize_to_fixpoint!(ReductionOptimizer(), graph) print("Preparing function... ") func = get_compute_function(graph, process, mock_machine()) + func(inputs[1]) print("Calculating... ") - Threads.@threads for i in 1:length(inputs) - results[i] += abs2(func(inputs[i])) + Threads.@threads for i in 1:N + Threads.@threads for j in 1:M + return results[i, j] += abs2(func(inputs[i, j])) + end end println("Done.") i += 1 @@ -157,20 +163,23 @@ for photons in [6] end end -# scenario +exit(0) + +# scenario 1 (disabled) n = 1000000 # n is the number of incoming photons # omega is the number -for photons in [6] +for photons in 1:6 # temp process to generate momenta for omega in [2e-3, 2e-6] println("Generating $n inputs for $photons photons...") temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n] - results = [0.0 for _ in 1:length(input_momenta)] + results = Array{Float64}(undef, size(input_momenta)) + fill!(results, 0.0) i = 1 for (in_pol, in_spin, out_pol, out_spin) in diff --git a/src/models/physics_models/qed/diagrams.jl b/src/models/physics_models/qed/diagrams.jl index 27082e5..0b79920 100644 --- a/src/models/physics_models/qed/diagrams.jl +++ b/src/models/physics_models/qed/diagrams.jl @@ -504,18 +504,15 @@ 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) 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], + [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( @@ -523,10 +520,10 @@ function gen_compton_diagram_from_order_one_side( 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, ), ) @@ -538,9 +535,9 @@ function gen_compton_diagram_from_order_one_side( 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) @@ -553,9 +550,9 @@ function gen_compton_diagram_from_order_one_side( 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) @@ -569,7 +566,6 @@ function gen_compton_diagram_from_order_one_side( add_tie!(new_diagram, FeynmanTie(ps[1], ps[2])) return new_diagram end -=# """ gen_compton_diagrams(n::Int, m::Int) -- 2.49.0 From 4eee23f08121fd383bae437901f5389cf4a1fa2f Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Tue, 9 Jul 2024 23:19:25 +0200 Subject: [PATCH 9/9] Run on GPU --- examples/congruent_in_ph.jl | 152 ++++++++++------------- src/MetagraphOptimization.jl | 3 - src/QEDprocesses_patch.jl | 46 ------- src/models/physics_models/qed/compute.jl | 2 +- 4 files changed, 66 insertions(+), 137 deletions(-) diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl index 19acc24..8b805b7 100644 --- a/examples/congruent_in_ph.jl +++ b/examples/congruent_in_ph.jl @@ -6,8 +6,13 @@ using QEDcore using QEDprocesses using Random using UUIDs + +using CUDA + +using NamedDims using CSV using JLD2 +using FlexiMaps RNG = Random.MersenneTwister(123) @@ -79,7 +84,7 @@ function congruent_input_momenta_scenario_2( # ---------- # now calculate the final_momenta from omega, cos_theta and phi - n = number_particles(processDescription, ParticleStateful{Incoming, Photon, SFourMomentum}) + n = number_particles(processDescription, Incoming(), Photon()) cos_theta = cos(theta) omega_prime = (n * omega) / (1 + n * omega * (1 - cos_theta)) @@ -108,109 +113,82 @@ end with_stacksize(f, n) = fetch(schedule(Task(f, n))) # scenario 2 -N = 1000 -M = 1000 +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:6 +for photons in 1:5 # temp process to generate momenta - for omega in [2e-3, 2e-6] - println("Generating $(N*M) inputs for $photons photons (Scenario 2 grid walk)...") - temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) + 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 = [ - congruent_input_momenta_scenario_2(temp_process, omega, theta, phi) for - (theta, phi) in Iterators.product(thetas, phis) - ] - results = Array{Float64}(undef, size(input_momenta)) - fill!(results, 0.0) + input_momenta = + Array{typeof(congruent_input_momenta_scenario_2(temp_process, omegas[1], thetas[1], phis[1]))}(undef, (K, N, M)) - i = 1 - for (in_pol, in_spin, out_pol, out_spin) in - Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()]) + 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 - 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 = build_psp.(Ref(process), input_momenta) - print("Preparing graph... ") - graph = gen_graph(process) - optimize_to_fixpoint!(ReductionOptimizer(), graph) - print("Preparing function... ") - func = get_compute_function(graph, process, mock_machine()) - func(inputs[1]) + cu_results = CuArray{Float64}(undef, size(input_momenta)) + fill!(cu_results, 0.0) - print("Calculating... ") + 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 - return results[i, j] += abs2(func(inputs[i, j])) + inputs[k, i, j] = build_psp(process, input_momenta[k, i, j]) end end - println("Done.") - i += 1 end + cu_inputs = CuArray(inputs) - println("Writing results") + 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()) - out_ph_moms = getindex.(getindex.(input_momenta, 2), 1) - out_el_moms = getindex.(getindex.(input_momenta, 2), 2) + print("Calculating... ") + ts = 32 + bs = Int64(length(cu_inputs) / 32) - @save "$(photons)_congruent_photons_omega_$(omega)_grid.jld2" out_ph_moms out_el_moms results - end -end - -exit(0) - -# scenario 1 (disabled) -n = 1000000 - -# n is the number of incoming photons -# omega is the number - -for photons in 1:6 - # temp process to generate momenta - for omega in [2e-3, 2e-6] - println("Generating $n inputs for $photons photons...") - temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp()) - - input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n] - results = Array{Float64}(undef, size(input_momenta)) - fill!(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 = build_psp.(Ref(process), input_momenta) - - print("Preparing graph... ") - # prepare function - graph = gen_graph(process) - optimize_to_fixpoint!(ReductionOptimizer(), graph) - print("Preparing function... ") - func = get_compute_function(graph, process, mock_machine()) - - print("Calculating... ") - Threads.@threads for i in 1:n - results[i] += abs2(func(inputs[i])) - end - 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) - - @save "$(photons)_congruent_photons_omega_$(omega).jld2" out_ph_moms out_el_moms results + 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 diff --git a/src/MetagraphOptimization.jl b/src/MetagraphOptimization.jl index 77a85eb..874de89 100644 --- a/src/MetagraphOptimization.jl +++ b/src/MetagraphOptimization.jl @@ -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.== diff --git a/src/QEDprocesses_patch.jl b/src/QEDprocesses_patch.jl index 1eba98c..282d713 100644 --- a/src/QEDprocesses_patch.jl +++ b/src/QEDprocesses_patch.jl @@ -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 diff --git a/src/models/physics_models/qed/compute.jl b/src/models/physics_models/qed/compute.jl index 01b1459..3650955 100644 --- a/src/models/physics_models/qed/compute.jl +++ b/src/models/physics_models/qed/compute.jl @@ -17,7 +17,7 @@ function input_expr(instance::GenericQEDProcess, name::String, psp_symbol::Symbo return Meta.parse( "ParticleValueSP( - $type(momentum($psp_symbol, $(construction_string(particle_direction(type))), $(construction_string(particle_species(type))), $index)), + $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))), )", -- 2.49.0