Pass sanity checks

src/computable_dags/compute.jl

@@ -19,7 +19,7 @@ end
 import ComputableDAGs: compute, compute_effort, children
 
 const e = sqrt(4π / 137.035999177)
-const VERTEX = 1#-1im * e * gamma()
+const VERTEX = -1im * e * gamma()
 
 compute_effort(::ComputeTask_BaseState) = 0
 compute_effort(::ComputeTask_Propagator) = 0
@@ -63,8 +63,7 @@ function compute(
     )
     return Propagated( # "propagated" because it goes directly into the next pair
         species,
-        1,
-        #state,
+        state,
         # bispinor, adjointbispinor, or lorentzvector
     )
 end
@@ -101,8 +100,7 @@ function compute(
     vp_mom = _vp_momentum(input.vp, input.psp)
     vp_species = particle_species(input.vp)
     inner = QEDbase.propagator(vp_species, vp_mom)
-    return 1
-    #return inner
+    return inner
 end
 
 struct Unpropagated{PARTICLE_T<:AbstractParticleType,VALUE_T}
@@ -190,8 +188,8 @@ function compute(::ComputeTask_CollectPairs, args::Vararg{N,T}) where {N,T}
     return sum(args)
 end
 function compute(::ComputeTask_CollectTriples, args::Vararg{N,T}) where {N,T}
-    println("$([args...])")
-    println("$(sum(args))")
+    #println("$([args...])")
+    #println("$(sum(args))")
     return sum(args)
 end
 function compute(::ComputeTask_SpinPolCumulation, args::Vararg{N,T}) where {N,T}

src/computable_dags/generation.jl

@@ -414,17 +414,37 @@ function relative_sign_pair(p1::VirtualParticle, p2::VirtualParticle)
         return false
     end
 
-    (s1, n1) = _canonical_index(p1)
-    (s2, n2) = _canonical_index(p2)
+    local vec::Vector{OPEN_FERMION_CYCLE_T}
 
-    new_cycle = if (s1 == :left && s2 == :right)
-        (n1, n2)
-    elseif (s1 == :right && s2 == :left)
-        (n2, n1)
-    else
-        @assert false
+    (l, el) = if p1.species == Electron
+        _canonical_index(p1)
+    elseif p2.species == Electron
+        _canonical_index(p2)
+    end
+    (r, po) = if p1.species == Positron
+        _canonical_index(p1)
+    elseif p2.species == Positron
+        _canonical_index(p2)
     end
 
+    vec = [p1.open_cycles..., p2.open_cycles..., (el, po)]
+
+    n = _count_closed_cycles(vec)
+
+    return n % 2 == 1
+end
+
+"""
+    relative_sign_triple(p2::VirtualParticle, p3::VirtualParticle, p3::VirtualParticle)
+
+Returns true if the Triple task combining electron p1, and positron p2 into a full diagram family needs a relative sign in the DAG.
+"""
+function relative_sign_triple(p1::VirtualParticle, p2::VirtualParticle, p3::VirtualParticle)
+    # p1 + p2 is the photon
+    n1 = relative_sign_pair(p1, p2)
+
+    new_cycle = (_canonical_index(p1)[2], _canonical_index(p2)[2])
+
     res = VirtualParticle(
         p1.proc,
         Photon(),
@@ -436,57 +456,14 @@ function relative_sign_pair(p1::VirtualParticle, p2::VirtualParticle)
         ),
     )
 
-    @assert make_up(p1, p2, res) "$p1 + $p2 - $res"
-
-    n_before = length(p1.open_cycles) + length(p2.open_cycles)
-    n_after = length(res.open_cycles)
-
-    if sort([p1.open_cycles..., p2.open_cycles...]) == res.open_cycles
-        # no change in the cycles, so the two new fermion/antifermions match up and close their cycle -> negate
-        @info "1 P [TRUE] $p1 + $p2 - $res"
-        return true
-    elseif n_before - 1 == n_after
-        # one cycle removed -> negate
-        @info "2 P [TRUE] $p1 + $p2 - $res"
-        return true
-    elseif n_before - 2 == n_after
-        # one cycle removed -> negate
-        @info "2 P [TRUE] $p1 + $p2 - $res"
-        return true
-    elseif n_before + 1 == n_after
-        # one cycle opened -> no negation
-        @info "3 P [FALS] $p1 + $p2 - $res"
-        return false
-    else
-        # number of cycles not changed but cycle itself changed -> no negation
-        @info "4 P [FALS] $p1 + $p2 - $res"
-        @assert n_before == n_after
-        return false
-    end
-end
-
-"""
-    relative_sign_triple(p1::VirtualParticle, p2::VirtualParticle, p3::VirtualParticle)
-
-Returns true if the Triple task combining p1, p2 and p3 into a full diagram family needs a relative sign in the DAG.
-"""
-function relative_sign_triple(p1::VirtualParticle, p2::VirtualParticle, p3::VirtualParticle)
-    # p1 is the photon, p2 + p3 is the "inverse" photon, question is whether p2 and p3 close a cycle
-    n1 = relative_sign_pair(p2, p3)
+    # cycles closed?
+    closed_cycle_count = _count_closed_cycles(
+        OPEN_FERMION_CYCLE_T[p3.open_cycles..., res.open_cycles...]
+    )
 
-    # cycle closed? -> negate
-    n2 = length(p1.open_cycles) == 1
-    if n2
-        #@info "  T [TRUE] $p1 + $p2 + $p3"
-    else
-        #@info "  T [FALS] $p1 + $p2 + $p3"
-    end
+    n2 = closed_cycle_count % 2 == 1
 
-    if n1 != n2 # xor
-        return true
-    else
-        return false
-    end
+    return n1 != n2 # xor
 end
 
 """
@@ -654,7 +631,7 @@ function ComputableDAGs.graph(proc::PROC) where {PROC<:AbstractProcessDefinition
                 continue
             end
 
-            negate = relative_sign_triple(ph, el, po)
+            negate = relative_sign_triple(el, po, ph)
 
             compute_triples = if negate
                 insert_node!(g, ComputeTask_TripleNegated())

src/diagrams/diagrams.jl

@@ -991,6 +991,68 @@ function reduce_cycles(vec::Vector{OPEN_FERMION_CYCLE_T})
     return vec
 end
 
+function _count_closed_cycles(vec::Vector{OPEN_FERMION_CYCLE_T})
+    closed_count = 0
+    # Adapted from reduce_cycles
+    while true
+        # Flag to check if changes occur
+        changed = false
+        new_vec = OPEN_FERMION_CYCLE_T[]
+        skip_indices = Set{Int}()
+
+        for i in 1:length(vec)
+            if i in skip_indices
+                continue
+            end
+
+            fused = false
+            for j in (i + 1):length(vec)
+                if j in skip_indices
+                    continue
+                end
+
+                a, b = vec[i]
+                c, d = vec[j]
+
+                # Fuse if b == c
+                if b == c
+                    push!(new_vec, (a, d))
+                    push!(skip_indices, j)
+                    changed = true
+                    fused = true
+                    break
+                    # Fuse if d == a
+                elseif d == a
+                    push!(new_vec, (c, b))
+                    push!(skip_indices, j)
+                    changed = true
+                    fused = true
+                    break
+                end
+            end
+
+            # Add the original tuple if it wasn't fused and numbers aren't equal
+            if !fused
+                if vec[i][1] == vec[i][2]
+                    closed_count += 1   # only here do we actually close a cycle
+                else
+                    push!(new_vec, vec[i])
+                end
+            end
+        end
+
+        # Update the vector
+        vec = new_vec
+
+        # Break if no changes
+        if !changed
+            break
+        end
+    end
+
+    return closed_count
+end
+
 function _find_cycles(left_ferms::Vector{Int}, right_ferms::Vector{Int})
     all_cycles = OPEN_FERMION_CYCLE_T[]
     for (l, r) in Iterators.zip(left_ferms, right_ferms)