QEDInput for differential cross-section interface

[AntonReinhard]

After some talking with @szabo137 we came up with the following idea for the public differential cross-section interface:

The current interfaces have to check that the input/output phasespaces have the correct dimensionalities to fit to the process definition. We suggest adding a QEDProcessInput type that includes all the necessary information of the process in its template parameters and also one point in the phasespace. Something like

struct QEDProcessInput{Process, N1, N2, N3, N4, N5, N6}
    proc::Process
    inFerms::SVector{N1, FermionStateful{Incoming}}
    outFerms::SVector{N2, FermionStateful{Outgoing}}
    inAntiferms::SVector{N3, AntiFermionStateful{Incoming}}
    outAntiferms::SVector{N4, AntiFermionStateful{Outgoing}}
    inPhotons::SVector{N5, PhotonStateful{Incoming}}
    outPhotons::SVector{N6, PhotonStateful{Outgoing}}
end

Then, a diffCS implementation taking just a QEDProcessInput and a model would be provided, assembling the process, in-phase-space-definition, in-phase-space-point, out-phase-space-definition, and out-phase-space-definition from them and dispatching to the internal implementation. The internal implementations could then be hidden (add _) and would no longer need to provide a safety check for the dimensionalities of the phase spaces.

It would make the usage very simple, by allowing calls such as:

julia> output = diffCS(model, input)
julia> outputs = diffCS.(model, inputs)
julia> cuOutputs = diffCS.(model, CuArray(inputs))

[szabo137]

I like the idea! However, there is at least two open questions on my side:

Type stability

I assume that FermionStateful is the combination of the respective fermion momentum, direction, and spin. However, this means, for instance, that the SVector in the field inFerm contains different element types, if, e.g., the spins are different: assume we have something like this

using QEDbase
using StaticArrays

struct FermionStateful{D<:ParticleDirection, S<:AbstractSpin}
  mom::SFourMomentum
end

Then we quickly run into performance issues:

setup

mom = rand(SFourMomentum)
s1 = FermionStateful{Incoming,SpinUp}(mom)
s2 = FermionStateful{Incoming,SpinDown}(mom)

v_stable = SVector(s1,s1)
v_unstable = SVector(s1,s2)

benchmark: indexing

using BenchmarkTools

julia> @benchmark $v_stable[1]

BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  2.000 ns … 11.292 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.042 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.065 ns ±  0.119 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁                  █                  █                    ▂
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ █
  2 ns         Histogram: log(frequency) by time     2.12 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark $v_unstable[1]

BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  17.667 ns … 593.917 ns  ┊ GC (min … max): 0.00% … 92.83%
 Time  (median):     17.916 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   18.984 ns ±  17.794 ns  ┊ GC (mean ± σ):  2.84% ±  2.95%

  ▆█▆▄▂                                                        ▂
  ██████▆▆██▇▄▃▃▁▁▁▄▅▅▅▅▄▆▆▇▇█████████▇▅▆▆▇▆▆▆▆▇▇██████▇▇▆▆▆▅▆ █
  17.7 ns       Histogram: log(frequency) by time      25.4 ns <

 Memory estimate: 48 bytes, allocs estimate: 1.

Notice in the unstable version, the indexing allocates memory.

One could just remove the spin/polarization from the stateful particles, and add them to the process definition.

Suggested solution

The particle state (generalized to arbitrary subtypes of AbstractParticleType)

struct ParticleState{D<:ParticleDirection,P<:AbstractParticleType}
  mom::SFourMomentum
end

The input struct should be renamed, may to just PhaseSpace:

struct PhaseSpace{
  PROC<:AbstractProcessDefinition,
  MODEL<:AbstractModelDefinition,
  PSDEF<:AbstractPhasespaceDefinition,
  N_IN_ELECTRON, 
  N_OUT_ELECTRON,
  N_IN_POSITRON, 
  N_OUT_POSITRON,
  N_IN_PHOTON,
  N_OUT_PHOTON,
  }

  proc::PROC
  model::MODEL
  PS_DEF::PSDEF

  in_electron::SVector{N_IN_ELECTRON,ParticleState{Incoming,Electron}}
  out_electron::SVector{N_OUT_ELECTRON,ParticleState{Outgoing,Electron}}
  in_positron::SVector{N_IN_POSITRON,ParticleState{Incoming,Positron}}
  out_positron::SVector{N_OUT_POSITRON,ParticleState{Outgoing,Positron}}
  in_photon::SVector{N_IN_PHOTON,ParticleState{Incoming,Photon}}
  out_photon::SVector{N_OUT_PHOTON,ParticleState{Outgoing,Photon}}
end

Here the spins of all particles are stored as type parameters in the definition of PROC, which makes the fields storing the particle states stable.

Construction

Having a field for each possible type of input and output particle makes the construction of an object of the QEDinput somewhat cumbersome. I think the unused fields (e.g. if there is no anti-particle involved in the process) could be just empty SVector. However, these empty fields need to be of the right type, which would look like

state_type = ParticleState{Incoming,Positron}
empty_field = SVector{0,state_type}()

This seems not very user-friendly. As far as I know, there ways to set default values for fields are very limited in Julia. From the top of my head, one could use Base.@kwdef, but this could be rather slow, because keywords tend to be slow.

[AntonReinhard]

Is it a good idea to include the Spin/Pol information in the process definition? At least in the diagram/formula generation algorithms I have seen, the spin/pol is disregarded at this stage. If the spins/pols are known beforehand, specific optimizations should be possible because, if I understand correctly, some generated diagrams would suddenly become impossible. That could be nice, but also not all samples would then be applicable, namely only samples with the same spin/pol distributions as the initial process.

I see two problems with your proposed solution:

1. From a usability perspective, I think it's very unintuitive that the spins/pols are not immediately part of a given particle. Related to what I wrote above.
2. Extending this definition to Muons and Tauons would add 8 more SVectors to the struct and 8 more type parameters. That seems infeasible.

I'm wondering if it might not be more useful to only differentiate between incoming versus outgoing particles in the phase space, and keep the info about particle type in the particle definition like this:

struct ParticleState
    mom::SFourMomentum
    type::Type{<:AbstractParticle}
    spinorpol::Type{<:AbstractSpinOrPolarization}
    dir::Type{<:ParticleDirection}
end

struct PhaseSpace{
  PROC<:AbstractProcessDefinition,
  MODEL<:AbstractModelDefinition,
  PSDEF<:AbstractPhasespaceDefinition,
  N_IN_PARTICLES,
  N_OUT_PARTICLES,
  }

  proc::PROC
  model::MODEL
  PS_DEF::PSDEF

  in_particles::SVector{N_IN_PARTICLES,ParticleState}
  out_particles::SVector{N_OUT_PARTICLES,ParticleState}
end

This would include all the necessary information and also be easily extensible or reusable in a similar way for non-QED particles. It would keep type stability since the actual type is always just ParticleState. The only issues I can see are that the particle direction is now duplicated (available here and as implied information through placement in in_particles or out_particles) and that the type information now takes up space in each ParticleState object. It's only an extra 24 Bytes though.

Dispatch can be used as before, and even has been used very similarly by QEDprocesses already in definitions like

propagator(particle_type::FermionLike, P::QEDbase.AbstractFourMomentum)

[szabo137]

This is solved by #51 and #57. Thanks a lot, @AntonReinhard for your help with this.