Discussion: algebra, and constructor

[mmikhasenko]

This issue discusses mechanisms to enable four-vector arithmetics with a customary types.

Proposal

Make a customary types (not part of this package) to explain how constructor works if given just four-vector components. That would enable basic operations like

• addition, subtraction,
• boosts and rotations
  for wide class of objects (only four-vector components are transformed)

Form messages below @grasph

1. We define an LVB.AbstractLorentzVector abtsract type for pure Lorentz vector, that is which do not have other properties than Lorentz vector components.
2. We define LVB concrete types of LVB.AbstractLorentzVector for each supported coordinate systems. These types will be used by LVB when it needs an LVB.AbstractLorentzVector (see example code below). At the same time, the concrete type ctor will acts as as accessor to the bare vector.
3. LVB overrides Base.:+(::Any, ::LVB.AbstractLorentzVector) and possibly Base.:+(::LVB.AbstractLorentzVector, ::Any), same for -. It is cleaner than in the previous definition example, because only one argument can have extra properties and therefore there is no ambiguities on how these properties are propagated.
4. LVB provides also Lorentz transformations and scalar products.

This will looks like this for Base.:+, assuming the code in the LVB module,

abstract type AbtractLorentzVector end

# Note: if we go that way, I would prefer to use the shorter name
# PxPyPzE for the vector and the longer name e.g., PxPyPzECoord
# for the coordinate system.
struct PxPyPzEVec{T} <: AbtractLorentzVector
px::T
py::T
pz::T
energy::T
end

# Constructor from any type compliant with the kinematic interface
PxPyPzEVec(q) = PxPyPzEVec(px(q), py(q), pz(q), energy(q))

const XYZTVec{T} = PxPyPzEVec{T}

px(p::PxPyPzEVec) = p.px
py(p::PxPyPzEVec) = p.py
pz(p::PxPyPzEVec) = p.pz
energy(p::PxPyPzEVec) = p.energy
coordinate_system(PxPyPzEVec) = PxPyPzE

# Similar definitions for PtEtaPhiEVec{T}, PtEtaPhiMVec{T}
# ....

function Base.:+(p::Any, q::AbstractLorentzVector)
  px = px(p) + px(q)
  py = py(p) + py(q)
  pz = pz(p) + pz(q)
  energy = energy(p) + energy(q)
  materializer(p)(PxPyPzEVec(px, py, pz, energy))
end

# And possibly we define also (not clear to me if we should):
Base.:+(q::AbstractLorentzVector, p::Any) = Base.:+(p::Any, q::AbstractLorentzVector) 

Note: I think we need to require in the kinematic interface that all component accessors (px(), py(), pz(), energy() in case of PxPyPzE coordinate) return elements of the same type. In the above code, the ctor of PxPyPzEVec(::Any) assumes this. We could relax the requirements by asking that either they have the same type or promote_rule's are implemented, but I'm not convinced we need this extra sophistication.

Earlier discussion:

• if we want to provide a default algebra implementation, then we need an abstract supertype is order to override the Base package operators. We could make inheritance from a LorentzVectorBase abstract type as optional.

• Constructor from a LorentzVector. The constructor will need to take an argument of any type. It's probably ok. The Table.jl interface has a materializer() method that returns the "constructor". Maybe there is a reason for doing this way?

[Moelf]

Let me write down three applications that can benefit from this interface (i.e. save code a downstream developer needs to write), they have different goals and thus demand different set of things from this LVBase package, we should think about what trade-off we want to make:

1. Particle gen

This is @szabo137 's original application, imagine you have:

struct GenParticle
    lv::LV
    pdgID:Int
    isSUSY::Bool
    ...
end

Clearly in this case, demanding + to be defined doesn't make sense, demanding "constructor from a LV" also doesn't make sense. So in this case, adding those two to the interface specification will burden the user. (in reality, they will probably just ignore your interface).

2. Jet clustering

This is @grasph 's original application in mind, Philippe wants to write a jet clustering algorithm that can take any jets satisfying the interface, which clearly needs + to be defined. (there's a workaround we also discussed, involved get_bare_lv()). For this application, having abstract type also helps, because otherwise Philippe's algorithm can't "restrict" the input to the function other than calling some interface function in a duck-type way (ala istable())

3. Analysis

This is an example I learned from @pviscone by looking at what scikit-hep has, when people read nanoAOD from CMS and loads up muon/electron etc., it's possible to defined the "behavior" of the collection. And one of the possible behavior involves adding muons and muons, that roughly looks like:

struct Muon <: CMSParticle
    lv::LV
    charge::Int8
end

in this case, it's reasonable to defined + to also add the electric charges of the particles (it's non-ambiguous and useful). But, constructor from arbitrary LV doesn't make sense, since a bare LV doesn't have charge. Also in this case, asking for sub-typing is too much, since collaboration / big project has their own primary type hierarchy logic.

Summary

One argument is to say "LV is LV, particle is particle", but I can sympathize with people who scratch their heads and think "it's just LV, why it's so damn complicated". After all we're doing some practical physics, not trying to use these LVs for mathematical proof...

Interface\Application	Generator	Jet clustering	Analysis
arithmetic	👎	👍	overload
constructor	👎	N/A	👎
sub typing	👎	👍	👎

[mmikhasenko]

thinking a bit more, I'd be in favor of keeping this package minimal: no arithmetic, no transformations, no constraints to constructor.
Defining +,-,*is not a big deal, and it will differ case to case.
We will have other packages, which know how to transform, and how to do algebra.

[szabo137]

I agree with you, therefore I proposed a compromise here. With this, the package stays light weight, yet keeps the possibility for library developer to require arithmetics via the AbstractLorentzVector.

[grasph]

About the constructor

The materializer() method mechanism of Tables.jl allows a user to provide the "constructor" without modifying the package that provide the table-like type they use even if the package does not provide it. I think we should use the same approach.

We need to add to the list of functions to implements for compliance with the interface (I will use LVB as shortcut of LorentzVectorBase):

LVB.materializer(p::MyLorentzVector)

and this method must return a Callable defined as:

LVB.materializer(p::MyLorentzVector)(q) ≡ r

with r a copy of p whose Lorentz vector was set to the Lorentz vector value of q.

This method will allow the development of packages that provides transformation of Lorentz-vector-like objects e.g. boost and rotation.

❓ Can the dependency of the callable returned by LVB.materializer(p::MyLorentzVector) on the actual instance of MyLorentzVector alter the performance ? If it is the case, we can require instead providing a method LVB.materializer(p::Type{MyLorentzVector}) that returns a callable taking two arguments, the MyLorentzVector p to modify and the object q containing the new Lorentz vector components.

About the arithmetic

The behaviour of a function acting on AbstractLorentzVectors and using the arithmetic defined by the subtypes will depend on their implementation, unless:

• arithmetic properties (commutativity, distributivity, etc.) is respected for each property of the subtypes. This is very restrictive e.g., it cannot be achieved for the PDG id of a particle.

• or we better define the behaviour of the operators in the other properties than kinematic and the implementation of the function use these definitions to achieve a well-defined behaviour.

An example of definition can be defining the sum p+q (resp. the difference p-q) as the object p with the Lorentz vector of q added (resp. subtracted) to the Lorentz vector of p.

In addition to the + and - operations we discussed at the hackatron, we need matrix multiplication to make it useful. This totals to 6 methods: +(::Particle, ::ALV), -(::Particle, ::ALV), -(::Particle), *(::Number, ::Particle), *(::Particle, ::Number), *(::AbstractLorentzMatrix, ::Particle), with ALV a shortcut for AbstractLorentzVector. Our specs will also need to define the AbstactLorentzMatrix type.

If I define a Particle type and desire to use the boost function of an algorithm package relying on the AbstractLorentzVector interface, the number of methods to implement to use the package will definitely discourage me and I will implement my own boost function.

Therefore, I would prefer other developers to not rely on me to implement the arithmetic. I don't have examples of algorithms, when the component accessors defined by our interface, together with the materializer() method, will not be sufficient.

I am not convinced by the "compromise". It comes from an attempt to reconcile two orthogonal views: subtyping AbstactLorentzVector to benefit from default overrides of Base operators (the original idea), and subtyping to define an interface.

@oliver was advocating that subtyping is absolutely needed. I would need an example where the inheritance-free interface does not work, in order to understand if this package can and should solve the case. It might be that what we actual need is an abstract type for pure Lorentz vectors.

I'm also adding @graeme-a-stewart to the discussion.

[mmikhasenko]

I wonder, conceptually, do we to want to implement transformations (boosts, rotations)? - naively, if I have an object that is a Lorentz vector, I know how to boost, and rotate. How often, different applications would need to do something different (additional) when boost_z, or rotate_y?

I'd benefit from having this functions inside the LVB package. Some mechanism like materializer would be needed indeed

@grasph, I do not fully follow, what is a type of q in

LVB.materializer(p::MyLorentzVector) = q->build_p_from_q(p,q)

i.e. how to pass the "minimal" Lorentz vector information to p. Could you explain, please?

For arithmetic, I do not see how we can define it, but I also do not fully follow the proposal.

[grasph]

@grasph, I do not fully follow, what is a type of q in

LVB.materializer(p::MyLorentzVector) = q->build_p_from_q(p,q)

q is an instance of a type compliant with the LVB kinematic interface. The implementation of build_p_from_q will use the LVB Lorentz vector component accessors (px(q), py(q),... or pt(q), eta(q),...).

For arithmetic, I do not see how we can define it, but I also do not fully follow the proposal.

In the definition example I gave, only the Lorentz vector properties are summed, and other properties, if any, are copied from the 1st argument, which takes a special place.

I would like to make a proposal which is for me more satisfactory.

1. We define an LVB.AbstractLorentzVector abtsract type for pure Lorentz vector, that is which do not have other properties than Lorentz vector components.
2. We define LVB concrete types of LVB.AbstractLorentzVector for each supported coordinate systems. These types will be used by LVB when it needs an LVB.AbstractLorentzVector (see example code below). At the same time, the concrete type ctor will acts as as accessor to the bare vector.
3. LVB overrides Base.:+(::Any, ::LVB.AbstractLorentzVector) and possibly Base.:+(::LVB.AbstractLorentzVector, ::Any), same for -. It is cleaner than in the previous definition example, because only one argument can have extra properties and therefore there is no ambiguities on how these properties are propagated.
4. LVB provides also Lorentz transformations and scalar products.

This will looks like this for Base.:+, assuming the code in the LVB module,

abstract type AbtractLorentzVector end

# Note: if we go that way, I would prefer to use the shorter name
# PxPyPzE for the vector and the longer name e.g., PxPyPzECoord
# for the coordinate system.
struct PxPyPzEVec{T} <: AbtractLorentzVector
px::T
py::T
pz::T
energy::T
end

# Constructor from any type compliant with the kinematic interface
PxPyPzEVec(q) = PxPyPzEVec(px(q), py(q), pz(q), energy(q))

const XYZTVec{T} = PxPyPzEVec{T}

px(p::PxPyPzEVec) = p.px
py(p::PxPyPzEVec) = p.py
pz(p::PxPyPzEVec) = p.pz
energy(p::PxPyPzEVec) = p.energy
coordinate_system(PxPyPzEVec) = PxPyPzE

# Similar definitions for PtEtaPhiEVec{T}, PtEtaPhiMVec{T}
# ....

function Base.:+(p::Any, q::AbstractLorentzVector)
  px = px(p) + px(q)
  py = py(p) + py(q)
  pz = pz(p) + pz(q)
  energy = energy(p) + energy(q)
  materializer(p)(PxPyPzEVec(px, py, pz, energy))
end

# And possibly we define also (not clear to me if we should):
Base.:+(q::AbstractLorentzVector, p::Any) = Base.:+(p::Any, q::AbstractLorentzVector) 

Note: I think we need to require in the kinematic interface that all component accessors (px(), py(), pz(), energy() in case of PxPyPzE coordinate) return elements of the same type. In the above code, the ctor of PxPyPzEVec(::Any) assumes this. We could relax the requirements by asking that either they have the same type or promote_rule's are implemented, but I'm not convinced we need this extra sophistication.

[mmikhasenko]

Can we resurrect the discussion? When is our next meeting?

[graeme-a-stewart]

Hello @mmikhasenko! Yes. I am about to make a 0.4.4 release of JetReconstruction.jl and then this item was the one that I wanted to tackle next.

[mmikhasenko]

I do not see objections to the approach proposed by @grasph
It's an extension to the current package that would enable algebraic operations.

I'm going to move the proposal to the header of the issue, and hope we find time to tackle it in the coming weeks.
Meanwhile,

1. I release the current version that is already functional (Streamlining and cleanup #24 Update readme and docs #26 was a great improvement, thanks @szabo137 ) as 0.1.0
2. the next release with the abstract type and arithmetics will have some breaking changes