Fix compton

Project.toml

@@ -6,8 +6,8 @@ version = "0.1.0"
 [deps]
 AccurateArithmetic = "22286c92-06ac-501d-9306-4abd417d9753"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 ComputableDAGs = "62933717-1c9d-4a3f-b06f-7ab7f17ca32d"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
 QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
 QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
@@ -17,6 +17,11 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TypeUtils = "c3b1956e-8857-4d84-9b79-890df85b1e67"
 
+[compat]
+ComputableDAGs = "0.1.1"
+QEDbase = "0.3"
+QEDcore = "0.2"
+
 [extras]
 ComputableDAGs = "62933717-1c9d-4a3f-b06f-7ab7f17ca32d"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"

src/QEDFeynman.jl

@@ -16,7 +16,6 @@ using Base.Threads
 export ParticleValue
 export ParticleA, ParticleB, ParticleC
 export ABCParticle, GenericABCProcess, ABCModel, PerturbativeABC
-export ComputeTaskABC_P
 export ComputeTaskABC_S1
 export ComputeTaskABC_S2
 export ComputeTaskABC_V
@@ -28,13 +27,12 @@ export parse_dag
 # QED model
 export FeynmanDiagram, FeynmanVertex, FeynmanTie, FeynmanParticle
 export QEDModel
-export ComputeTaskQED_P
 export ComputeTaskQED_S1
 export ComputeTaskQED_S2
 export ComputeTaskQED_V
 export ComputeTaskQED_U
 export ComputeTaskQED_Sum
-export gen_graph
+export graph
 
 export ParticleValue, ParticleValueSP
 

src/abc/compute.jl

@@ -18,19 +18,6 @@ function ComputableDAGs.input_expr(
     )
 end
 
-"""
-    compute(::ComputeTaskABC_P, data::ParticleValue)
-
-Return the particle and value as is. 
-
-0 FLOP.
-"""
-function ComputableDAGs.compute(
-    ::ComputeTaskABC_P, data::ParticleValue{P}
-)::ParticleValue{P} where {P}
-    return data
-end
-
 """
     compute(::ComputeTaskABC_U, data::ParticleValue)
 

src/abc/parse.jl

@@ -93,19 +93,21 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
     insert_edge!(graph, sum_node, global_data_out)
 
     # remember the data out nodes for connection
-    dataOutNodes = Dict()
+    data_out_nodes = Dict()
 
     if (verbose)
         println("Building graph")
     end
-    noNodes = 0
+    number_of_nodes = 0
     nodesToRead = length(nodes)
     while !isempty(nodes)
         node = popfirst!(nodes)
-        noNodes += 1
-        if (noNodes % 100 == 0)
+        number_of_nodes += 1
+        if (number_of_nodes % 100 == 0)
             if (verbose)
-                percent = string(round(100.0 * noNodes / nodesToRead; digits=2), "%")
+                percent = string(
+                    round(100.0 * number_of_nodes / nodesToRead; digits=2), "%"
+                )
                 print("\rReading Nodes... $percent")
             end
         end
@@ -114,18 +116,14 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
 
             # add nodes and edges for the state reading to u(P(Particle))
             data_in = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE), name) # read particle data node
-            compute_P = insert_node!(graph, ComputeTaskABC_P()) # compute P node
-            data_Pu = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE)) # transfer data from P to u (one ParticleValue object)
             compute_u = insert_node!(graph, ComputeTaskABC_U()) # compute U node
             data_out = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE)) # transfer data out from u (one ParticleValue object)
 
-            insert_edge!(graph, data_in, compute_P)
-            insert_edge!(graph, compute_P, data_Pu)
-            insert_edge!(graph, data_Pu, compute_u)
+            insert_edge!(graph, data_in, compute_u)
             insert_edge!(graph, compute_u, data_out)
 
             # remember the data_out node for future edges
-            dataOutNodes[node] = data_out
+            data_out_nodes[node] = data_out
         elseif occursin(regex_c, node)
             capt = match(regex_c, node)
 
@@ -140,12 +138,12 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
                 compute_S = insert_node!(graph, ComputeTaskABC_S1())
                 data_S_v = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE))
 
-                insert_edge!(graph, dataOutNodes[in1], compute_S)
+                insert_edge!(graph, data_out_nodes[in1], compute_S)
                 insert_edge!(graph, compute_S, data_S_v)
 
                 insert_edge!(graph, data_S_v, compute_v)
             else
-                insert_edge!(graph, dataOutNodes[in1], compute_v)
+                insert_edge!(graph, data_out_nodes[in1], compute_v)
             end
 
             if (occursin(regex_c, in2))
@@ -154,16 +152,16 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
                 compute_S = insert_node!(graph, ComputeTaskABC_S1())
                 data_S_v = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE))
 
-                insert_edge!(graph, dataOutNodes[in2], compute_S)
+                insert_edge!(graph, data_out_nodes[in2], compute_S)
                 insert_edge!(graph, compute_S, data_S_v)
 
                 insert_edge!(graph, data_S_v, compute_v)
             else
-                insert_edge!(graph, dataOutNodes[in2], compute_v)
+                insert_edge!(graph, data_out_nodes[in2], compute_v)
             end
 
             insert_edge!(graph, compute_v, data_out)
-            dataOutNodes[node] = data_out
+            data_out_nodes[node] = data_out
 
         elseif occursin(regex_m, node)
             # assume for now that only the first particle of the three is combined and the other two are "original" ones
@@ -176,16 +174,16 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
             compute_v = insert_node!(graph, ComputeTaskABC_V())
             data_v = insert_node!(graph, DataTask(PARTICLE_VALUE_SIZE))
 
-            insert_edge!(graph, dataOutNodes[in2], compute_v)
-            insert_edge!(graph, dataOutNodes[in3], compute_v)
+            insert_edge!(graph, data_out_nodes[in2], compute_v)
+            insert_edge!(graph, data_out_nodes[in3], compute_v)
             insert_edge!(graph, compute_v, data_v)
 
             # combine with the v of the combined other input
             compute_S2 = insert_node!(graph, ComputeTaskABC_S2())
             data_out = insert_node!(graph, DataTask(FLOAT_SIZE)) # output of a S2 task is only a float
 
             insert_edge!(graph, data_v, compute_S2)
-            insert_edge!(graph, dataOutNodes[in1], compute_S2)
+            insert_edge!(graph, data_out_nodes[in1], compute_S2)
             insert_edge!(graph, compute_S2, data_out)
 
             insert_edge!(graph, data_out, sum_node)
@@ -201,7 +199,7 @@ function parse_dag(filename::AbstractString, proc::GenericABCProcess, verbose::B
     end
 
     #put all nodes into dirty nodes set
-    graph.dirtyNodes = copy(graph.nodes)
+    graph.dirty_nodes = copy(graph.nodes)
 
     if (verbose)
         println("Generating the graph's properties")

src/abc/properties.jl

@@ -28,13 +28,6 @@ Return the compute effort of a V task.
 """
 compute_effort(t::ComputeTaskABC_V)::Float64 = 6.0
 
-"""
-    compute_effort(t::ComputeTaskABC_P)
-
-Return the compute effort of a P task.
-"""
-compute_effort(t::ComputeTaskABC_P)::Float64 = 0.0
-
 """
     compute_effort(t::ComputeTaskABC_Sum)
 
@@ -59,13 +52,6 @@ Return the number of children of a ComputeTaskABC_S2 (always 2).
 """
 children(::ComputeTaskABC_S2) = 2
 
-"""
-    children(::ComputeTaskABC_P)
-
-Return the number of children of a ComputeTaskABC_P (always 1).
-"""
-children(::ComputeTaskABC_P) = 1
-
 """
     children(::ComputeTaskABC_U)
 

src/abc/types.jl

@@ -6,11 +6,11 @@ Singleton definition for identification of the ABC-Model.
 struct ABCModel <: AbstractPhysicsModel end
 
 """
-    PerturbativeABC <: AbstractModel
+    PerturbativeABC <: AbstractModelDefinition
 
 The model being used for the ABC model.
 """
-struct PerturbativeABC <: AbstractModelDefinition end
+struct PerturbativeABC <: QEDbase.AbstractModelDefinition end
 
 """
     ABCParticle
@@ -54,13 +54,6 @@ S task with two children.
 """
 struct ComputeTaskABC_S2 <: AbstractComputeTask end
 
-"""
-    ComputeTaskABC_P <: AbstractComputeTask
-
-P task with no children.
-"""
-struct ComputeTaskABC_P <: AbstractComputeTask end
-
 """
     ComputeTaskABC_V <: AbstractComputeTask
 
@@ -92,7 +85,6 @@ Constant vector of all tasks of the ABC-Model.
 ABC_TASKS = [
     ComputeTaskABC_S1,
     ComputeTaskABC_S2,
-    ComputeTaskABC_P,
     ComputeTaskABC_V,
     ComputeTaskABC_U,
     ComputeTaskABC_Sum,

src/interface.jl

@@ -5,7 +5,7 @@ import QEDbase.AbstractParticle
 
 Base type for a model, e.g. ABC-Model or QED. This is used to dispatch many functions.
 """
-abstract type AbstractPhysicsModel <: AbstractModel end
+abstract type AbstractPhysicsModel end
 
 """
     ParticleValue{ParticleType <: AbstractParticleStateful}
@@ -99,9 +99,9 @@ Return the model of this process description or input.
 function model end
 
 """
-    type_from_name(model::AbstractModel, name::String)
+    type_from_name(model::AbstractPhysicsModel, name::String)
 
-For a name of a particle in the given `AbstractModel`, return the particle's `Type` and index as a tuple. The input string can be expetced to be of the form \"<name><index>\".
+For a name of a particle in the given `AbstractPhysicsModel`, return the particle's `Type` and index as a tuple. The input string can be expetced to be of the form \"<name><index>\".
 """
 function type_index_from_name end
 

src/qed/compute.jl

@@ -51,15 +51,20 @@ function ComputableDAGs.compute(
 ) where {P1<:ParticleStateful,P2<:ParticleStateful,V1<:ValueType,V2<:ValueType}
     p3 = QED_conserve_momentum(data1.p, data2.p)
     state = QED_vertex()
-    if (typeof(data1.v) <: AdjointBiSpinor)
-        state = data1.v * state
+    if data1.v isa AdjointBiSpinor
+        @assert !(data2.v isa AdjointBiSpinor)
+        state = data1.v * state * data2.v
+    elseif data1.v isa BiSpinor
+        @assert !(data2.v isa BiSpinor)
+        state = data2.v * state * data1.v
+    elseif data2.v isa AdjointBiSpinor
+        @assert !(data1.v isa AdjointBiSpinor)
+        state = data2.v * state * data1.v
+    elseif data2.v isa BiSpinor
+        @assert !(data1.v isa BiSpinor)
+        state = data1.v * state * data2.v
     else
-        state = state * data1.v
-    end
-    if (typeof(data2.v) <: AdjointBiSpinor)
-        state = data2.v * state
-    else
-        state = state * data2.v
+        @assert "invalid V task inputs"
     end
 
     dataOut = ParticleValue(p3, state)
@@ -88,9 +93,12 @@ function ComputableDAGs.compute(
     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)"
+    inner1 = QED_inner_edge(data1.p)
+    inner2 = QED_inner_edge(data2.p)
 
-    inner = QED_inner_edge(propagated_particle(data1.p))
+    # TODO: This is broken. It's currently not possible (i think) to find out which of the two is the correct
+    # inner edge value to take. Likely the graph building has to be changed to provide the correct one first or similar.
+    inner = is_outgoing(data1.p) ? inner1 : inner2
 
     # inner edge is just a "scalar", data1 and data2 are bispinor/adjointbispinnor, need to keep correct order
     if typeof(data1.v) <: BiSpinor
@@ -106,8 +114,12 @@ function ComputableDAGs.compute(
     data2::ParticleValue{ParticleStateful{D2,Photon},V2},
 ) where {D1<:ParticleDirection,D2<:ParticleDirection,V1<:ValueType,V2<:ValueType}
     # TODO: assert that data1 and data2 are opposites
+    @assert isapprox(
+        momentum(data1.p), momentum(data2.p), rtol=sqrt(eps()), atol=sqrt(eps())
+    ) "$(momentum(data1.p)) vs. $(momentum(data2.p))"
+
     inner = QED_inner_edge(data1.p)
-    # inner edge is just a scalar, data1 and data2 are photon states that are just Complex numbers here
+
     return data1.v * inner * data2.v
 end
 
@@ -119,37 +131,38 @@ Compute inner edge (1 input particle, 1 output particle).
 function ComputableDAGs.compute(
     ::ComputeTaskQED_S1, data::ParticleValue{P,V}
 ) where {P<:ParticleStateful,V<:ValueType}
+    inner = QED_inner_edge(data.p)
     new_p = propagated_particle(data.p)
     # inner edge is just a scalar, can multiply from either side
     if typeof(data.v) <: BiSpinor
-        return ParticleValue(new_p, QED_inner_edge(new_p) * data.v)
+        return ParticleValue(new_p, inner * data.v)
     else
-        return ParticleValue(new_p, data.v * QED_inner_edge(new_p))
+        return ParticleValue(new_p, data.v * inner)
     end
 end
 
 """
     compute(::ComputeTaskQED_Sum, data...)
     compute(::ComputeTaskQED_Sum, data::AbstractArray)
 
-Compute a sum over the vector. Use an algorithm that accounts for accumulated errors in long sums with potentially large differences in magnitude of the summands.
+Compute a sum over the vector and return the `abs2()` of it.
 
 Linearly many FLOP with growing data.
 """
-function ComputableDAGs.compute(::ComputeTaskQED_Sum, data...)::ComplexF64
+function ComputableDAGs.compute(::ComputeTaskQED_Sum, data...)
     # TODO: want to use sum_kbn here but it doesn't seem to support ComplexF64, do it element-wise?
     s = 0.0im
     for d in data
         s += d
     end
-    return s
+    return abs2(s)
 end
 
-function ComputableDAGs.compute(::ComputeTaskQED_Sum, data::AbstractArray)::ComplexF64
+function ComputableDAGs.compute(::ComputeTaskQED_Sum, data::AbstractArray)
     # TODO: want to use sum_kbn here but it doesn't seem to support ComplexF64, do it element-wise?
     s = 0.0im
     for d in data
         s += d
     end
-    return s
+    return abs2(s)
 end