Add support for synced spins/pols (#31)

* Improve accuracy of finestructure constant

* Remove QEDprocesses dependency in main Project.toml

* Remove some unnecessary type parameters

* Implement Synced Spins/Pols support

* Add tests for synced spin and pol, update literate.jl example

* Add test for synced Spins and multiple independent syncs

* Add QED dev dependencies for CI

.github/workflows/BuildDeployDoc.yml

@@ -22,8 +22,7 @@ jobs:
       - name: Install dependencies # TODO: currently needs dev versions of QEDbase and QEDprocesses, remove once they are released again
         run: |
           julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
-          julia --project=docs/ -e 'using Pkg; Pkg.add(url="https://github.com/QEDjl-project/QEDbase.jl/")'
-          julia --project=docs/ -e 'using Pkg; Pkg.add(url="https://github.com/QEDjl-project/QEDprocesses.jl/")'
+          julia --project=docs/ ./add_QED_dev.jl
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token

.github/workflows/unit_tests.yml

@@ -26,9 +26,10 @@ jobs:
         with:
           version: "1.10"
 
-      - name: Instantiate
+      - name: Instantiate # TODO: currently needs dev versions of QEDbase and QEDprocesses, remove once they are released again
         run: |
           julia --project=./ -e 'using Pkg; Pkg.instantiate()'
+          julia --project=./ ./add_QED_dev.jl
 
       - name: Run tests
         run: julia --project=./ -t 4 -e 'using Pkg; Pkg.test()' -O0

.gitignore

@@ -35,3 +35,6 @@ Manifest.toml
 .julia
 **/.ipynb_checkpoints/
 *.bkp
+
+*.gif
+*.jld2

Project.toml

@@ -11,7 +11,6 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
 QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
-QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
@@ -26,8 +25,11 @@ Reexport = "1"
 julia = "1.10"
 
 [extras]
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["SafeTestsets", "Test"]
+test = ["SafeTestsets", "Test", "Random", "QEDprocesses", "RuntimeGeneratedFunctions"]

add_QED_dev.jl

@@ -0,0 +1,7 @@
+using Pkg
+
+Pkg.add(; url="https://github.com/QEDjl-project/QEDbase.jl", rev="dev")
+@warn "This repository depends on the dev branch of QEDbase.jl\n It is NOT ready for release!"
+
+Pkg.add(; url="https://github.com/QEDjl-project/QEDprocesses.jl.git", rev="dev")
+@warn "This repository depends on the dev branch of QEDprocesses.jl\n It is NOT ready for release!"

docs/src/examples/compton.jl

@@ -16,17 +16,13 @@ using QEDprocesses
 n = 4;
 
 # Now we setup the scattering process accordingly. We consider all spin/polarization
-# combinations of the particles except for the incoming photons.
-#
-# !!! note
-#     The functionality is not quite ready yet, but in the future, the incoming photons
-#     can be synced using `SyncedPol` from QEDbase. This would emulate all incoming
-#     photons having the same polarization, for example from a laser.
+# combinations of the particles except for the incoming photons, where the polarizations are synced using [`QEDbase.SyncedPolarization`](@extref). 
+# This emulates all synced photons having the same, but still indefinite, polarization, for example from a laser.
 proc = ScatteringProcess(
-    (Electron(), ntuple(_ -> Photon(), n)...),  # incoming particles
-    (Electron(), Photon()),                     # outgoing particles
-    (AllSpin(), ntuple(_ -> PolX(), n)...),     # incoming particle spin/pols
-    (AllSpin(), AllPol()),                      # outgoing particle spin/pols
+    (Electron(), ntuple(_ -> Photon(), n)...),     # incoming particles
+    (Electron(), Photon()),                        # outgoing particles
+    (AllSpin(), ntuple(_ -> SyncedPol(1), n)...),  # incoming particle spin/pols
+    (AllSpin(), AllPol()),                         # outgoing particle spin/pols
 )
 
 # The [`feynman_diagrams`](@ref) function returns an iterator for all possible Feynman diagrams

docs/src/lib/internal.md

@@ -3,6 +3,9 @@
 ```@meta
 CurrentModule = QEDFeynmanDiagrams
 ```
+```@index
+Pages = ["internal.md"]
+```
 
 ## Types
 
@@ -13,9 +16,10 @@ AbstractTreeLevelFeynmanDiagram
 ## Functions
 
 ```@docs
-disjunct
+_is_index_valid_combination
+_pseudo_virtual_particles
+are_total
 contains
+disjunct
 make_up
-are_total
-_pseudo_virtual_particles
 ```

docs/src/lib/public.md

@@ -1,5 +1,9 @@
 # Public Interface
 
+```@index
+Pages = ["public.md"]
+```
+
 ## Types
 
 ```@docs
@@ -9,9 +13,9 @@ VirtualParticle
 
 ## Functions
 ```@docs
-feynman_diagrams
 external_particles
-virtual_particles
-process
+feynman_diagrams
 generate_DAG
+process
+virtual_particles
 ```

src/computable_dags/compute.jl

@@ -16,7 +16,7 @@ end
 # import so we don't have to repeat it all the time
 import ComputableDAGs: compute, compute_effort, children
 
-const e = sqrt(4π / 137)
+const e = sqrt(4π / 137.035999177)
 const VERTEX = -1im * e * gamma()
 
 compute_effort(::ComputeTask_BaseState) = 0
@@ -42,9 +42,7 @@ struct BaseStateInput{PS_T<:AbstractParticleStateful,SPIN_POL_T<:AbstractSpinOrP
     spin_pol::SPIN_POL_T
 end
 
-function compute(
-    ::ComputeTask_BaseState, input::BaseStateInput{PS,SPIN_POL}
-) where {PS,SPIN_POL}
+function compute(::ComputeTask_BaseState, input::BaseStateInput)
     species = particle_species(input.particle)
     if is_outgoing(input.particle)
         species = _invert(species)
@@ -82,7 +80,7 @@ end
 end
 
 function _vp_momentum(
-    vp::VirtualParticle{PROC,SPECIES,I,O}, psp::PhaseSpacePoint
+    vp::VirtualParticle{PROC,SPECIES,I,O}, psp::AbstractPhaseSpacePoint
 ) where {PROC,SPECIES,I,O}
     return _masked_sum(momenta(psp, Incoming()), _in_contributions(vp)) -
            _masked_sum(momenta(psp, Outgoing()), _out_contributions(vp))
@@ -113,48 +111,46 @@ end
 
 @inline function compute( # photon, electron
     ::ComputeTask_Pair,
-    photon::Propagated{Photon,V1},
-    electron::Propagated{Electron,V2},
-) where {V1,V2}
+    photon::Propagated{Photon},
+    electron::Propagated{Electron},
+)
     return Unpropagated(Electron(), photon.value * VERTEX * electron.value) # photon - electron -> electron
 end
 @inline function compute( # photon, positron
     ::ComputeTask_Pair,
-    photon::Propagated{Photon,V1},
-    positron::Propagated{Positron,V2},
-) where {V1,V2}
+    photon::Propagated{Photon},
+    positron::Propagated{Positron},
+)
     return Unpropagated(Positron(), positron.value * VERTEX * photon.value) # photon - positron -> positron
 end
 @inline function compute( # electron, positron
     ::ComputeTask_Pair,
-    electron::Propagated{Electron,V1},
-    positron::Propagated{Positron,V2},
-) where {V1,V2}
+    electron::Propagated{Electron},
+    positron::Propagated{Positron},
+)
     return Unpropagated(Photon(), positron.value * VERTEX * electron.value)  # electron - positron -> photon
 end
 
-@inline function compute(
-    ::ComputeTask_PropagatePairs, prop::PROP_V, photon::Unpropagated{Photon,VAL}
-) where {PROP_V,VAL}
+@inline function compute(::ComputeTask_PropagatePairs, prop, photon::Unpropagated{Photon})
     return Propagated(Photon(), photon.value * prop)
 end
 @inline function compute(
-    ::ComputeTask_PropagatePairs, prop::PROP_V, electron::Unpropagated{Electron,VAL}
-) where {PROP_V,VAL}
+    ::ComputeTask_PropagatePairs, prop, electron::Unpropagated{Electron}
+)
     return Propagated(Electron(), prop * electron.value)
 end
 @inline function compute(
-    ::ComputeTask_PropagatePairs, prop::PROP_V, positron::Unpropagated{Positron,VAL}
-) where {PROP_V,VAL}
+    ::ComputeTask_PropagatePairs, prop, positron::Unpropagated{Positron}
+)
     return Propagated(Positron(), positron.value * prop)
 end
 
 @inline function compute(
     ::ComputeTask_Triple,
-    photon::Propagated{Photon,V1},
-    electron::Propagated{Electron,V2},
-    positron::Propagated{Positron,V3},
-) where {V1,V2,V3}
+    photon::Propagated{Photon},
+    electron::Propagated{Electron},
+    positron::Propagated{Positron},
+)
     return positron.value * (VERTEX * photon.value) * electron.value
 end
 

src/computable_dags/generation.jl

@@ -122,6 +122,8 @@ function ComputableDAGs.input_expr(
 end
 
 function ComputableDAGs.input_type(p::AbstractProcessDefinition)
+    # TODO correctly assemble abstract types here
+    # See https://github.com/QEDjl-project/QEDFeynmanDiagrams.jl/issues/29
     return Any
     in_t = QEDcore._assemble_tuple_type(incoming_particles(p), Incoming(), SFourMomentum)
     out_t = QEDcore._assemble_tuple_type(outgoing_particles(p), Outgoing(), SFourMomentum)
@@ -163,10 +165,14 @@ end
 _dir_str(::Incoming) = "inc"
 _dir_str(::Outgoing) = "out"
 
+# the possible spins or pols for generating base state tasks
 _spin_pols(::AllSpin) = (SpinUp(), SpinDown())
+_spin_pols(::SyncedSpin) = (SpinUp(), SpinDown())
 _spin_pols(::SpinUp) = (SpinUp(),)
 _spin_pols(::SpinDown) = (SpinDown(),)
+
 _spin_pols(::AllPol) = (PolX(), PolY())
+_spin_pols(::SyncedPol) = (PolX(), PolY())
 _spin_pols(::PolX) = (PolX(),)
 _spin_pols(::PolY) = (PolY(),)
 
@@ -307,19 +313,71 @@ _edge_index_from_species(::Electron) = 2
 _edge_index_from_species(::Positron) = 3
 _edge_index_from_vp(vp::VirtualParticle) = _edge_index_from_species(particle_species(vp))
 
+"""
+    _is_index_valid_combination(proc::AbstractProcessDefinition, index::Tuple)
+
+Internal function for DAG generation. Checks for a given process and a spin/pol combination whether the spin/pol combination is
+part of the process, including checking for [`QEDbase.SyncedPolarization`](@extref) and [`QEDbase.SyncedSpin`](@extref).
+"""
+function _is_index_valid_combination(proc::AbstractProcessDefinition, index::Tuple)
+    proc_spin_pols = (incoming_spin_pols(proc)..., outgoing_spin_pols(proc)...)
+
+    # for synced spins/pols, remember the first occurrence and its definite spin/pol, then check that later ones are the same
+    synced_pols = Dict{SyncedPol,AbstractDefinitePolarization}()
+    synced_spins = Dict{SyncedSpin,AbstractDefiniteSpin}()
+
+    for (i, sp) in index
+        if proc_spin_pols[i] isa AllSpin || proc_spin_pols[i] isa AllPol
+            # for allspin and allpol, everything is allowed
+            continue
+        end
+        if proc_spin_pols[i] == sp
+            # sp is always definite, so if they're equal the combination is always allowed
+            continue
+        end
+        if proc_spin_pols[i] isa SyncedSpin
+            if !haskey(synced_spins, proc_spin_pols[i]) # insert first occurrence
+                synced_spins[proc_spin_pols[i]] = sp
+                continue
+            end
+            if synced_spins[proc_spin_pols[i]] == sp # otherwise, check if sp is synced
+                continue
+            end
+            # the spin is not synced
+            return false
+        end
+
+        if proc_spin_pols[i] isa SyncedPol
+            if !haskey(synced_pols, proc_spin_pols[i])
+                synced_pols[proc_spin_pols[i]] = sp
+                continue
+            end
+            if synced_pols[proc_spin_pols[i]] == sp
+                continue
+            end
+            # the pol is not synced
+            return false
+        end
+
+        error("encountered unknown spin or polarization type")
+    end
+
+    return true
+end
+
 """
     generate_DAG(proc::AbstractProcessDefinition)
 
 Generate and return a [`ComputableDAGs.DAG`](@extref), representing the computation for the squared matrix element of this scattering process, summed over spin and polarization combinations allowed by the process.
 """
 function generate_DAG(proc::AbstractProcessDefinition)
     particles = virtual_particles(proc)                  # virtual particles that will be input to propagator tasks
+    # TODO apparently this sort is deprecated, change it
     pairs = sort(particle_pairs(particles))              # pairs to generate the pair tasks
     triples = sort(total_particle_triples(particles))    # triples to generate the triple tasks
 
     graph = DAG()
 
-    # TODO: use the spin/pol iterator here once it has been implemented
     # -- Base State Tasks --
     propagated_outputs = Dict{VirtualParticle,Vector{Node}}()
     for dir in (Incoming(), Outgoing())
@@ -394,7 +452,16 @@ function generate_DAG(proc::AbstractProcessDefinition)
             )
 
             for in_nodes in Iterators.product(particles_data_out_nodes...)
+                # get the spin/pol config of the input particles from the data_out names
+                index = _parse_node_names(in_nodes[1].name, in_nodes[2].name)
+                # index is a tuple of tuples, containing the particle index and their definite spin/pol
+                if !_is_index_valid_combination(proc, index)
+                    # skip this pair creation if the spin/pol combination doesn't exist
+                    continue
+                end
+
                 # make the compute pair nodes for every combination of the found input_particle_nodes to get all spin/pol combinations
+
                 compute_pair = insert_node!(graph, ComputeTask_Pair())
                 pair_data_out = insert_node!(graph, DataTask(0))
 
@@ -412,9 +479,6 @@ function generate_DAG(proc::AbstractProcessDefinition)
                 )
                 insert_edge!(graph, compute_pair, pair_data_out)
 
-                # get the spin/pol config of the input particles from the data_out names
-                index = _parse_node_names(in_nodes[1].name, in_nodes[2].name)
-
                 if !haskey(pair_output_nodes_by_spin_pol, index)
                     pair_output_nodes_by_spin_pol[index] = Vector()
                 end
@@ -458,6 +522,12 @@ function generate_DAG(proc::AbstractProcessDefinition)
         positrons = propagated_outputs[po]
 
         for (a, b, c) in Iterators.product(photons, electrons, positrons) # for each spin/pol config of each part
+            index = _parse_node_names(a.name, b.name, c.name)
+            if !_is_index_valid_combination(proc, index)
+                # skip this triple creation if the spin/pol combination doesn't exist, same as for pairs
+                continue
+            end
+
             compute_triples = insert_node!(graph, ComputeTask_Triple())
             data_triples = insert_node!(graph, DataTask(0))
 
@@ -467,7 +537,6 @@ function generate_DAG(proc::AbstractProcessDefinition)
 
             insert_edge!(graph, compute_triples, data_triples)
 
-            index = _parse_node_names(a.name, b.name, c.name)
             if !haskey(triples_results, index)
                 triples_results[index] = Vector{DataTaskNode}()
             end

test/runtests.jl

@@ -1,6 +1,5 @@
 using SafeTestsets
 
-@safetestset "Dummy" begin
-    # TODO: There are no tests yet
-    @test 1 == 1
+@safetestset "Synced Spins and Polarizations" begin
+    include("synced_spin_pol.jl")
 end