Welcome to QRules!¶

QRules is a system for validating and generating particle reactions, using quantum number conservation rules. The user only has to provide a basic information of the particle reaction, such as an initial state and a final state. Helper functions provide easy ways to configure the system, but the user still has full control. QRules then constructs several hypotheses for what happens during the transition from initial to final state.

Original project: PWA Expert System

The original project was the Welcome to the PWA Expert System!. QRules originates from its ~expertsystem.reaction module.

Internal design¶

Internally, QRules consists of three major components.

1. State Transition Graphs¶

A StateTransitionGraph is a directed graph that consists of nodes and edges. In a directed graph, each edge must be connected to at least one node (in correspondence to Feynman graphs). This way, a graph describes the transition from one state to another.

The edges correspond to particles/states, in other words a collection of properties such as the quantum numbers that characterize the particle state.
Each node represents an interaction and contains all information for the transition of this specific step. Most importantly, a node contains a collection of conservation rules that have to be satisfied. An interaction node has \(M\) ingoing lines and \(N\) outgoing lines, where \(M,N \in \mathbb{Z}\), \(M > 0, N > 0\).

2. Conservation Rules¶

The central component of the expert system are the conservation rules. They belong to individual nodes and receive properties about the node itself, as well as properties of the ingoing and outgoing edges of that node. Based on those properties the conservation rules determine whether edges pass or not.

3. Solvers¶

The determination of the correct state properties in the graph is done by solvers. New properties are set for intermediate edges and interaction nodes and their validity is checked with the conservation rules.

QRules workflow¶

Preparation

1.1. Build all possible topologies. A topology is represented by a graph, in which the edges and nodes are empty (no particle information).

1.2. Fill the topology graphs with the user provided information. Typically these are the graph’s ingoing edges (initial state) and outgoing edges (final state).
Solving

2.1. Propagate quantum number information through the complete graph while respecting the specified conservation laws. Information like mass is not used in this first solving step.

2.2. Clone graphs while inserting concrete matching particles for the intermediate edges (mainly adds the mass variable).

2.3. Validate the complete graphs, so run all conservation law check that were postponed from the first step.

Table of Contents¶

Installation¶

The fastest way of installing this package is through PyPI:

python3 -m pip install qrules

This installs the latest, stable release that you can find on the stable branch. The latest version on the main branch can be installed as follows:

python3 -m pip install git+https://github.com/ComPWA/qrules@main

In that case, however, we highly recommend using the more dynamic, ‘editable installation’ instead. This goes as follows:

Get the source code (see Working with Git):

git clone https://github.com/ComPWA/qrules.git
cd qrules

[Recommended] Create a virtual environment (see here).

Install the project in ‘editable installation’, as well as additional dependencies for the developer:

# pin dependencies first!
python3 -m pip install -r reqs/PYTHON_VERSION/requirements-dev.txt
python3 -m pip install -e .

That’s all! Have a look at the Usage page to try out the package, and see Develop for tips on how to work with this ‘editable’ developer setup!

Usage¶

Main interface¶

Here’s a small example of how to use qrules!

Investigate intermediate resonances¶

import qrules as q

result = q.generate_transitions(
    initial_state="J/psi(1S)",
    final_state=["K0", "Sigma+", "p~"],
    allowed_interaction_types="strong",
    formalism_type="canonical-helicity",
)

import graphviz

dot = q.io.asdot(result, collapse_graphs=True)
graphviz.Source(dot)

Next, you use the ampform package to convert these transitions into a mathematical description that you can use to fit your data and perform Partial Wave Analysis!

Quantum number search¶

The load_pdg() function creates a ParticleCollection containing the latest PDG info. Its find() and filter() methods allows you to quickly look up the quantum numbers of a particle and, vice versa, look up particle candidates based on a set of quantum numbers.

import qrules as q

pdg = q.load_pdg()
pdg.find(22)
delta = pdg.find("Delta(1920)++")
q.io.asdict(delta)

{'name': 'Delta(1920)++',
 'pid': 22224,
 'latex': '\\Delta(1920)^{++}',
 'spin': 1.5,
 'mass': 1.92,
 'width': 0.3,
 'charge': 2,
 'isospin': {'magnitude': 1.5, 'projection': 1.5},
 'baryon_number': 1,
 'parity': {'value': 1}}

subset = pdg.filter(
    lambda p: p.spin in [2.5, 3.5, 4.5] and p.name.startswith("N")
)
subset.names

{'N(1675)+',
 'N(1675)0',
 'N(1675)~-',
 'N(1675)~0',
 'N(1680)+',
 'N(1680)0',
 'N(1680)~-',
 'N(1680)~0',
 'N(2190)+',
 'N(2190)0',
 'N(2190)~-',
 'N(2190)~0'}

Tip

See Particle database

Check allowed reactions¶

qrules can be used to check whether a transition between an initial and final state is violated by any conservation rules:

q.check_reaction_violations(
    initial_state="pi0",
    final_state=["gamma", "gamma", "gamma"],
)

{frozenset({'c_parity_conservation'})}

Advanced¶

Each of the qrules’s sub-modules offer functionality to handle more advanced reaction types. The following notebooks illustrate how use them.

Generate transitions¶

A Partial Wave Analysis starts by defining an amplitude model that describes the reaction process that is to be analyzed. Such a model is generally very complex and requires a fair amount of effort by the analyst (you). This gives a lot of room for mistakes.

The ‘expert system’ is responsible to give you advice on the form of an amplitude model, based on the problem set you define (initial state, final state, allowed interactions, intermediate states, etc.). Internally, the system propagates the quantum numbers through the reaction graph while satisfying the specified conservation rules. How to control this procedure is explained in more detail below.

Afterwards, the amplitude model of the expert system can be exported into Tensorwaves. The model can for instance be used to generate a data set (toy Monte Carlo) for this reaction and to optimize its parameters to resemble an actual data set as good as possible. For more info on that see Formulate amplitude model.

Note

Simple channels can be treated with the generate_transitions() façade function. This notebook shows how to treat more complicated cases with the StateTransitionManager.

1. Define the problem set¶

We first define the boundary conditions of our physics problem, such as initial state, final state, formalism type, etc. and pass all of that information to the StateTransitionManager. This is the main user interface class of qrules.

By default, the StateTransitionManager loads all particles from the PDG. The qrules would take a long time to check the quantum numbers of all these particles, so in this notebook, we use a smaller subset of relatively common particles.

from qrules import InteractionTypes, StateTransitionManager

stm = StateTransitionManager(
    initial_state=["J/psi(1S)"],
    final_state=["gamma", "pi0", "pi0"],
    formalism_type="helicity",
)

StateTransitionManager

The StateTransitionManager (STM) is the main user interface class of {mod}`qrules`. The boundary conditions of your physics problem, such as the initial state, final state, formalism type, etc., are defined through this interface.

create_problem_sets of the STM creates all problem sets ― using the boundary conditions of the StateTransitionManager instance. In total 3 steps are performed. The creation of reaction topologies. The creation of InitialFacts, based on a topology and the initial and final state information. And finally the solving settings such as the conservation laws and quantum number domains to use at which point of the topology.
By default, all three interaction types (EM, STRONG, and WEAK) are used in the preparation stage. However, it is also possible to choose the allowed interaction types globally via set_allowed_interaction_types.

After the preparation step, you can modify the problem sets returned by create_problem_sets to your liking. Since the output of this function contains quite a lot of information, {mod}`qrules` aids in the configuration (especially the STM).

A subset of particles that are allowed as intermediate states can also be specified: either through the STM's constructor or by setting the instance allowed_intermediate_particles.

2. Prepare Problem Sets¶

Create all ProblemSet’s using the boundary conditions of the StateTransitionManager instance. By default it uses the isobar model (tree of two-body decays) to build Topology’s. Various InitialFacts are created for each topology based on the initial and final state. Lastly some reasonable default settings for the solving process are chosen. Remember that each interaction node defines its own set of conservation laws.

The StateTransitionManager (STM) defines three interaction types:

Interaction	Strength
strong	\(60\)
electromagnetic (EM)	\(1\)
weak	\(10^{-4}\)

By default, all three are used in the preparation stage. The create_problem_sets() method of the STM generates graphs with all possible combinations of interaction nodes. An overall interaction strength is assigned to each graph and they are grouped according to this strength.

problem_sets = stm.create_problem_sets()

3. Find solutions¶

If you are happy with the default settings generated by the StateTransitionManager, just start with solving directly!

This step takes about 23 sec on an Intel(R) Core(TM) i7-6820HQ CPU of 2.70GHz running, multi-threaded.

result = stm.find_solutions(problem_sets)

The find_solutions() method returns a Result object from which you can extract the transitions. Now, you can use get_intermediate_particles() to print the names of the intermediate states that the StateTransitionManager found:

print("found", len(result.transitions), "solutions!")
result.get_intermediate_particles().names

found 414 solutions!

{'a(0)(1450)0',
 'a(0)(980)0',
 'a(1)(1260)0',
 'a(1)(1640)0',
 'a(2)(1320)0',
 'a(2)(1700)0',
 'b(1)(1235)0',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'f(0)(500)',
 'f(0)(980)',
 'f(1)(1285)',
 'f(1)(1420)',
 "f(2)'(1525)",
 'f(2)(1270)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(1420)',
 'omega(1650)',
 'omega(782)',
 'phi(1020)',
 'phi(1680)',
 'rho(1450)0',
 'rho(1700)0',
 'rho(770)0'}

About the number of solutions

The “number of transitions” is the total number of allowed StateTransitionGraph instances that the StateTransitionManager has found. This also includes all allowed spin projection combinations. In this channel, we for example consider a \(J/\psi\) with spin projection \(\pm1\) that decays into a \(\gamma\) with spin projection \(\pm1\), which already gives us four possibilities.

On the other hand, the intermediate state names that was extracted with Result.get_intermediate_particles(), is just a set of the state names on the intermediate edges of the list of transitions, regardless of spin projection.

Now we have a lot of solutions that are actually heavily suppressed (they involve two weak decays).

In general, you can modify the ProblemSets returned by create_problem_sets() directly, but the STM also comes with functionality to globally choose the allowed interaction types. So, go ahead and disable the EM and InteractionTypes.WEAK interactions:

stm.set_allowed_interaction_types([InteractionTypes.STRONG])
problem_sets = stm.create_problem_sets()
result = stm.find_solutions(problem_sets)

print("found", len(result.transitions), "solutions!")
result.get_intermediate_particles().names

found 192 solutions!

{'b(1)(1235)0',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'f(0)(500)',
 'f(0)(980)',
 "f(2)'(1525)",
 'f(2)(1270)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'rho(1450)0',
 'rho(1700)0',
 'rho(770)0'}

Now note that, since a \(\gamma\) particle appears in one of the interaction nodes, the expert system knows that this node must involve EM interactions! Because the node can be an effective interaction, the weak interaction cannot be excluded, as it contains only a subset of conservation laws.

Since only the strong interaction was supposed to be used, this results in a warning and the STM automatically corrects the mistake.

Once the EM interaction is included, this warning disappears. Be aware, however, that the EM interaction is now available globally. Hence, there now might be solutions in which both nodes are electromagnetic.

stm.set_allowed_interaction_types(
    [InteractionTypes.STRONG, InteractionTypes.EM]
)
problem_sets = stm.create_problem_sets()
result = stm.find_solutions(problem_sets)

print("found", len(result.transitions), "solutions!")
result.get_intermediate_particles().names

found 318 solutions!

{'a(0)(1450)0',
 'a(0)(980)0',
 'a(2)(1320)0',
 'a(2)(1700)0',
 'b(1)(1235)0',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'f(0)(500)',
 'f(0)(980)',
 "f(2)'(1525)",
 'f(2)(1270)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)',
 'h(1)(1170)',
 'omega(1420)',
 'omega(1650)',
 'omega(782)',
 'phi(1020)',
 'phi(1680)',
 'rho(1450)0',
 'rho(1700)0',
 'rho(770)0'}

Great! Now we selected only the strongest contributions. Be aware, though, that there are more effects that can suppress certain decays, like small branching ratios. In this example, the initial state \(J/\Psi\) can decay into \(\pi^0 + \rho^0\) or \(\pi^0 + \omega\).

decay	branching ratio
\(\omega\to\gamma+\pi^0\)	0.0828
\(\rho^0\to\gamma+\pi^0\)	0.0006

Unfortunately, the \(\rho^0\) mainly decays into \(\pi^0+\pi^0\), not \(\gamma+\pi^0\) and is therefore suppressed. This information is currently not known to the expert system, but it is possible to hand the expert system a list of allowed intermediate states.

# particles are found by name comparison,
# i.e. f2 will find all f2's and f all f's independent of their spin
stm.set_allowed_intermediate_particles(["f(0)", "f(2)"])

result = stm.find_solutions(problem_sets)

print("found", len(result.transitions), "solutions!")
result.get_intermediate_particles().names

found 138 solutions!

{'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'f(0)(500)',
 'f(0)(980)',
 "f(2)'(1525)",
 'f(2)(1270)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)'}

Now we have selected all amplitudes that involve f states. The warnings appear only to notify the user that the list of solutions is not exhaustive: for certain edges in the graph, no suitable particle was found (since only f states were allowed).

import graphviz
from qrules import io

dot = io.asdot(result, collapse_graphs=True, render_node=False)
graphviz.Source(dot)

See also

Visualize decay topologies

4. Export generated transitions¶

The Result, StateTransitionGraph, and Topology can be serialized to and from a dict with io.asdict() and io.fromdict():

from qrules import io

graph = result.transitions[0]
io.asdict(graph.topology)

{'nodes': [0, 1],
 'edges': {-1: {'ending_node_id': 0},
  0: {'originating_node_id': 0},
  1: {'originating_node_id': 1},
  2: {'originating_node_id': 1},
  3: {'originating_node_id': 0, 'ending_node_id': 1}}}

This also means that the Result can be written to JSON or YAML format with io.write() and loaded again with io.load():

io.write(result, "transitions.json")
imported_result = io.load("transitions.json")
assert imported_result == result

Handy if it takes a lot of computation time to generate the transitions!

Warning

It’s not possible to pickle a Result, because StateTransitionGraph makes use of Generic.

Tip

The ampform package can formulate amplitude models based on the state transitions created by qrules.

Particle database¶

In PWA, you usually want to search for special resonances, possibly even some not listed in the PDG. In this notebook, we go through a few ways to add or overwrite Particle instances in the database with your own particle definitions.

Loading the default database¶

In Generate transitions, we made use of the StateTransitionManager. By default, if you do not specify the particles argument, the StateTransitionManager calls the function load_default_particles(). This functions returns a ParticleCollection instance with Particle definitions from the PDG, along with additional definitions that are provided in the file additional_definitions.yml.

Here, we call this method directly to illustrate what happens (we use load_pdg(), which loads a subset):

from qrules.particle import load_pdg

particle_db = load_pdg()
print("Number of loaded particles:", len(particle_db))

Number of loaded particles: 519

In the following, we illustrate how to use the methods of the ParticleCollection class to find and ‘modify’ Particles and add() them back to the ParticleCollection.

Finding particles¶

The ParticleCollection class offers some methods to search for particles by name or by PID (see find()):

particle_db.find(333)

Particle(
  name='phi(1020)',
  pid=333,
  latex='\\phi(1020)',
  spin=1.0,
  mass=1.019461,
  width=0.004248999999999999,
  isospin=Spin(0, 0),
  parity=-1,
  c_parity=-1,
  g_parity=-1,
)

With filter(), you can perform more sophisticated searches. This is done by either passing a function or lambda.

subset = particle_db.filter(lambda p: p.name.startswith("f(2)"))
subset.names

{"f(2)'(1525)",
 'f(2)(1270)',
 'f(2)(1950)',
 'f(2)(2010)',
 'f(2)(2300)',
 'f(2)(2340)'}

subset = particle_db.filter(
    lambda p: p.strangeness == 1
    and p.spin >= 1
    and p.mass > 1.8
    and p.mass < 1.9
)
subset.names

{'K(2)(1820)+',
 'K(2)(1820)0',
 'Lambda(1820)~',
 'Lambda(1830)~',
 'Lambda(1890)~'}

subset = particle_db.filter(lambda p: p.is_lepton())
subset.names

{'e+',
 'e-',
 'mu+',
 'mu-',
 'nu(e)',
 'nu(e)~',
 'nu(mu)',
 'nu(mu)~',
 'nu(tau)',
 'nu(tau)~',
 'tau+',
 'tau-'}

Note that in each of these examples, we call the names property. This is just to only display the names, sorted alphabetically, otherwise the output becomes a bit of a mess:

particle_db.filter(lambda p: p.name.startswith("pi") and len(p.name) == 3)

ParticleCollection({
  Particle(
    name='pi0',
    pid=111,
    latex='\\pi^{0}',
    spin=0.0,
    mass=0.1349768,
    width=7.73e-09,
    isospin=Spin(1, 0),
    parity=-1,
    c_parity=+1,
    g_parity=-1,
  ),
  Particle(
    name='pi-',
    pid=-211,
    latex='\\pi^{-}',
    spin=0.0,
    mass=0.13957039000000002,
    width=2.5284e-17,
    charge=-1,
    isospin=Spin(1, -1),
    parity=-1,
    g_parity=-1,
  ),
  Particle(
    name='pi+',
    pid=211,
    latex='\\pi^{+}',
    spin=0.0,
    mass=0.13957039000000002,
    width=2.5284e-17,
    charge=1,
    isospin=Spin(1, +1),
    parity=-1,
    g_parity=-1,
  ),
})

LaTeX representation¶

Particles also contain a latex tag. Here, we use ipython to render them nicely as mathematical symbols:

from IPython.display import Math

sigmas = particle_db.filter(
    lambda p: p.name.startswith("Sigma") and p.charmness == 1
)
Math(", ".join([p.latex for p in sigmas]))

\[\displaystyle \Sigma_{c}(2520)^{+}, \Sigma_{c}(2455)^{++}, \Sigma_{c}^{0}, \Sigma_{c}(2520)^{0}, \Sigma_{c}(2455)^{+}, \Sigma_{c}(2520)^{++}\]

Adding custom particle definitions through Python¶

A quick way to modify or overwrite particles, is through your Python script or notebook. Notice that the instances in the database are Particle instances:

N1650_plus = particle_db["N(1650)+"]
N1650_plus

Particle(
  name='N(1650)+',
  pid=32212,
  latex='N(1650)^{+}',
  spin=0.5,
  mass=1.65,
  width=0.125,
  charge=1,
  isospin=Spin(1/2, +1/2),
  baryon_number=1,
  parity=-1,
)

The instances in the database are immutable. Therefore, if you want to modify, say, the width, you have to create a new Particle instance from the particle you want to modify and add() it back to the database. You can do this with create_particle():

from qrules.particle import create_particle

new_N1650_plus = create_particle(
    template_particle=N1650_plus, name="Modified N(1650)+", width=0.2
)

particle_db.add(new_N1650_plus)
particle_db["Modified N(1650)+"].width

0.2

You often also want to add the antiparticle of the particle you modified to the database. Using create_antiparticle(), it is easy to create the corresponding antiparticle object.

from qrules.particle import create_antiparticle

new_N1650_minus = create_antiparticle(
    new_N1650_plus, new_name="Modified N(1650)-"
)

particle_db.add(new_N1650_minus)
particle_db["Modified N(1650)-"]

Particle(
  name='Modified N(1650)-',
  pid=-32212,
  latex='\\overline{N(1650)^{+}}',
  spin=0.5,
  mass=1.65,
  width=0.2,
  charge=-1,
  isospin=Spin(1/2, -1/2),
  baryon_number=-1,
  parity=+1,
)

When adding additional particles you may need for your research, it is easiest to work with an existing particle as template. Let’s say we want to study \(e^+e^-\) collisions of several energies:

energies_mev = {4180, 4220, 4420, 4600}
template_particle = particle_db["J/psi(1S)"]
for energy_mev in energies_mev:
    energy_gev = energy_mev / 1e3
    new_particle = create_particle(
        template_particle, name=f"EpEm ({energy_mev} MeV)", mass=energy_gev
    )
    particle_db.add(new_particle)
len(particle_db)

particle_db.filter(lambda p: "EpEm" in p.name).names

{'EpEm (4180 MeV)', 'EpEm (4220 MeV)', 'EpEm (4420 MeV)', 'EpEm (4600 MeV)'}

Of course, it’s also possible to add any kind of custom Particle, as long as its quantum numbers comply with the gellmann_nishijima() rule:

from qrules.particle import Particle

custom = Particle(
    name="custom",
    pid=99999,
    latex=R"p_\mathrm{custom}",
    spin=1.0,
    mass=1,
    charge=1,
    isospin=(1.5, 0.5),
    charmness=1,
)
custom

Particle(
  name='custom',
  pid=99999,
  latex='p_\\mathrm{custom}',
  spin=1.0,
  mass=1.0,
  charge=1,
  isospin=Spin(3/2, +1/2),
  charmness=1,
)

particle_db += custom
len(particle_db)

Loading custom definitions from a YAML file¶

It’s also possible to add particles from a config file, with io.load(). Existing entries remain and if the imported file of particle definitions contains a particle with the same name, it is overwritten in the database.

It’s easiest to work with YAML. Here, we use the provided additional_particles.yml example file:

from qrules import io

particle_db += io.load("additional_particles.yml")

Writing to YAML¶

You can also dump the existing particle lists to YAML. You do this with the io.write() function.

io.write(instance=particle_db, filename="dumped_particle_list.yaml")

Note that the function write can dump any ParticleCollection to an output file, also a specific subset.

from qrules.particle import ParticleCollection

output = ParticleCollection()
output += particle_db["J/psi(1S)"]
output += particle_db.find(22)  # gamma
output += particle_db.filter(lambda p: p.name.startswith("f(0)"))
output += particle_db["pi0"]
output += particle_db["pi+"]
output += particle_db["pi-"]
output += particle_db["custom"]
io.write(output, "particle_list_selection.yml")
output.names

{'J/psi(1S)',
 'custom',
 'f(0)(1370)',
 'f(0)(1500)',
 'f(0)(1710)',
 'f(0)(500)',
 'f(0)(980)',
 'gamma',
 'pi+',
 'pi-',
 'pi0'}

As a side note, qrules provides JSON schemas (reaction/particle-validation.json) to validate your particle list files (see also jsonschema.validate()). If you have installed qrules as an Editable installation and use VSCode, your YAML particle list are checked automatically in the GUI.

Visualize decay topologies¶

The io module allows you to convert StateTransitionGraph and Topology instances to DOT language with asdot(). You can visualize its output with third-party libraries, such as Graphviz. This is particularly useful after running find_solutions(), which produces a Result object with a list of StateTransitionGraph instances (see Generate transitions).

Topologies¶

First of all, here are is an example of how to visualize a group of Topology instances. We use create_isobar_topologies() and create_n_body_topology() to create a few standard topologies.

import graphviz
from qrules import io
from qrules.topology import create_isobar_topologies, create_n_body_topology

topology = create_n_body_topology(2, 4)
graphviz.Source(io.asdot(topology, render_initial_state_id=True))

Note the IDs of the nodes is also rendered if there is more than node:

topologies = create_isobar_topologies(4)
graphviz.Source(io.asdot(topologies))

This can be turned on or off with the arguments of asdot():

topologies = create_isobar_topologies(3)
graphviz.Source(io.asdot(topologies, render_node=False))

asdot() provides other options as well:

topologies = create_isobar_topologies(5)
dot = io.asdot(
    topologies[0],
    render_final_state_id=False,
    render_resonance_id=True,
    render_node=False,
)
display(graphviz.Source(dot))

`StateTransitionGraph`s¶

Here, we’ll visualize the allowed transitions for the decay \(\psi' \to \gamma\eta\eta\) as an example.

import qrules as q

result = q.generate_transitions(
    initial_state="psi(2S)",
    final_state=["gamma", "eta", "eta"],
    allowed_interaction_types="EM",
)

As noted in 3. Find solutions, the transitions contain all spin projection combinations (which is necessary for the ampform package). It is possible to convert all these solutions to DOT language with asdot(). To avoid visualizing all solutions, we just take a subset of the transitions:

dot = q.io.asdot(result.transitions[::50][:3])  # just some selection

This str of DOT language for the list of StateTransitionGraph instances can then be visualized with a third-party library, for instance, with graphviz.Source:

import graphviz

dot = q.io.asdot(
    result.transitions[::50][:3], render_node=False
)  # just some selection
graphviz.Source(dot)

You can also serialize the DOT string to file with io.write(). The file extension for a DOT file is .gv:

q.io.write(result, "decay_topologies_with_spin.gv")

Collapse graphs¶

Since this list of all possible spin projections transitions is rather long, it is often useful to use strip_spin=True or collapse_graphs=True to bundle comparable graphs. First, strip_spin=True allows one collapse (ignore) the spin projections (we again show a selection only):

dot = q.io.asdot(result.transitions[:3], strip_spin=True)
graphviz.Source(dot)

Note

By default, .asdot renders edge IDs, because they represent the (final) state IDs as well. In the example above, we switched this off.

If that list is still too much, there is collapse_graphs=True, which bundles all graphs with the same final state groupings:

dot = q.io.asdot(result, collapse_graphs=True, render_node=False)
graphviz.Source(dot)

Bibliography¶

Tip

Download this bibliography as BibTeX here.

1: D. M. Asner. Dalitz Plot Analysis Formalism. In Review of Particle Physics: Volume I Reviews. January 2006. pdg.lbl.gov/2010/reviews/rpp2010-rev-dalitz-analysis-formalism.pdf.
2: S.-U. Chung et al. Partial wave analysis in 𝐾-matrix formalism. Annalen der Physik, 507(5):404–430, May 1995. doi:10.1002/andp.19955070504.

qrules¶

import qrules

A rule based system that facilitates particle reaction analysis.

QRules generates allowed particle transitions from a set of conservation rules and boundary conditions as specified by the user. The user boundary conditions for a particle reaction problem are for example the initial state, final state, and allowed interactions.

The core of qrules computes which transitions (represented by a StateTransitionGraph) are allowed between a certain initial and final state. Internally, the system propagates the quantum numbers defined by the particle module through the StateTransitionGraph, while satisfying the rules define by the conservation_rules module. See Generate transitions and Particle database.

Finally, the io module provides tools that can read and write the objects of this framework.

check_reaction_violations(initial_state: Union[str, Tuple[str, Sequence[float]], Sequence[Union[str, Tuple[str, Sequence[float]]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], mass_conservation_factor: Optional[float] = 3.0) → Set[FrozenSet[str]][source]¶

Determine violated interaction rules for a given particle reaction.

Warning

This function only guarantees to find P, C and G parity violations, if it’s a two body decay. If all initial and final states have the C/G parity defined, then these violations are also determined correctly.

Parameters

initial_state – Shortform description of the initial state w/o spin projections.
final_state – Shortform description of the final state w/o spin projections.
mass_conservation_factor – Factor with which the width is multiplied when checking for MassConservation. Set to None in order to deactivate mass conservation.

Returns

Set of least violating rules. The set can have multiple entries, as several quantum numbers can be violated. Each entry in the frozenset represents a group of rules that together violate all possible quantum number configurations.

generate_transitions(initial_state: Union[str, Tuple[str, Sequence[float]], Sequence[Union[str, Tuple[str, Sequence[float]]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], allowed_intermediate_particles: Optional[List[str]] = None, allowed_interaction_types: Optional[Union[str, List[str]]] = None, formalism_type: str = 'helicity', particles: Optional[qrules.particle.ParticleCollection] = None, mass_conservation_factor: Optional[float] = 3.0, topology_building: str = 'isobar', number_of_threads: Optional[int] = None) → qrules.transition.Result [source]¶

Generate allowed transitions between an initial and final state.

Serves as a facade to the StateTransitionManager (see Generate transitions).

Parameters

initial_state (list) – A list of particle names in the initial state. You can specify spin projections for these particles with a tuple, e.g. ("J/psi(1S)", [-1, 0, +1]). If spin projections are not specified, all projections are taken, so the example here would be equivalent to "J/psi(1S)".
final_state (list) – Same as initial_state, but for final state particles.
allowed_intermediate_particles (list, optional) – A list of particle states that you want to allow as intermediate states. This helps (1) filter out resonances and (2) speed up computation time.
allowed_interaction_types (str, optional) – Interaction types you want to consider. For instance, both "strong and EM" and ["s", "em"] results in EM and STRONG.
formalism_type (str, optional) – Formalism that you intend to use in the eventual amplitude model.
particles (ParticleCollection, optional) – The particles that you want to be involved in the reaction. Uses load_pdg by default. It’s better to use a subset for larger reactions, because of the computation times. This argument is especially useful when you want to use your own particle definitions (see Particle database).
mass_conservation_factor – Width factor that is taken into account for for the MassConservation rule.
topology_building (str) –
Technique with which to build the Topology instances. Allowed values are:
- "isobar": Isobar model (each state decays into two states)
- "nbody": Use one central node and connect initial and final states to it
number_of_threads (int) – Number of cores with which to compute the allowed transitions. Defaults to all cores on the system.

An example (where, for illustrative purposes only, we specify all arguments) would be:

>>> import qrules as q
>>> result = q.generate_transitions(
...     initial_state="D0",
...     final_state=["K~0", "K+", "K-"],
...     allowed_intermediate_particles=["a(0)(980)", "a(2)(1320)-"],
...     allowed_interaction_types="ew",
...     formalism_type="helicity",
...     particles=q.load_pdg(),
...     topology_building="isobar",
... )
>>> len(result.transitions)
4

load_default_particles() → qrules.particle.ParticleCollection [source]¶

Load the default particle list that comes with qrules.

Runs load_pdg and supplements its output definitions from the file additional_definitions.yml.

Submodules and Subpackages

io¶

import qrules.io

Serialization module for the qrules.

The io module provides tools to export or import objects from qrules to and from disk, so that they can be used by external packages, or just to store (cache) the state of the system.

asdict(instance: object) → dict [source]¶

asdot(instance: object, *, render_node: bool = True, render_final_state_id: bool = True, render_resonance_id: bool = False, render_initial_state_id: bool = False, strip_spin: bool = False, collapse_graphs: bool = False) → str [source]¶

Convert a object to a DOT language str.

Only works for objects that can be represented as a graph, particularly a StateTransitionGraph or a list of StateTransitionGraph instances.

See also

Visualize decay topologies

fromdict(definition: dict) → object [source]¶

load(filename: str) → object [source]¶

write(instance: object, filename: str) → None [source]¶

argument_handling¶

import qrules.argument_handling

Handles argument handling for rules.

Responsibilities are the check of requirements for rules and the creation of the arguments from general graph property maps. The information is extracted from the type annotations of the rules.

class RuleArgumentHandler[source]¶

Bases: object

register_rule(rule: Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule]) → Tuple[Callable, Callable][source]¶

get_required_qns(rule: Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule]) → Tuple[Set[Type[Union[qrules.quantum_numbers.EdgeQuantumNumbers.pid, qrules.quantum_numbers.EdgeQuantumNumbers.mass, qrules.quantum_numbers.EdgeQuantumNumbers.width, qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.charge, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.strangeness, qrules.quantum_numbers.EdgeQuantumNumbers.charmness, qrules.quantum_numbers.EdgeQuantumNumbers.bottomness, qrules.quantum_numbers.EdgeQuantumNumbers.topness, qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number, qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.parity, qrules.quantum_numbers.EdgeQuantumNumbers.c_parity, qrules.quantum_numbers.EdgeQuantumNumbers.g_parity]]], Set[Type[Union[qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.l_projection, qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.s_projection, qrules.quantum_numbers.NodeQuantumNumbers.parity_prefactor]]]][source]¶

combinatorics¶

import qrules.combinatorics

Perform permutations on the edges of a StateTransitionGraph.

In a StateTransitionGraph, the edges represent quantum states, while the nodes represent interactions. This module provides tools to permutate, modify or extract these edge and node properties.

class InitialFacts(edge_props: Dict[int, Tuple[qrules.particle.Particle, float]] = NOTHING, node_props: Dict[int, qrules.quantum_numbers.InteractionProperties] = NOTHING)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class InitialFacts.

edge_props: Dict[int, Tuple[qrules.particle.Particle, float]]¶

node_props: Dict[int, qrules.quantum_numbers.InteractionProperties]¶

create_initial_facts(topology: qrules.topology.Topology, particles: qrules.particle.ParticleCollection, initial_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], final_state_groupings: Optional[Union[List[List[List[str]]], List[List[str]], List[str]]] = None) → List[qrules.combinatorics.InitialFacts][source]¶

match_external_edges(graphs: List[qrules.topology.StateTransitionGraph[Tuple[qrules.particle.Particle, float]]]) → None [source]¶

perform_external_edge_identical_particle_combinatorics(graph: qrules.topology.StateTransitionGraph) → List[qrules.topology.StateTransitionGraph][source]¶

Create combinatorics clones of the StateTransitionGraph.

In case of identical particles in the initial or final state. Only identical particles, which do not enter or exit the same node allow for combinatorics!

conservation_rules¶

import qrules.conservation_rules

Collection of quantum number conservation rules for particle reactions.

This module is the place where the ‘expert’ defines the rules that verify quantum numbers of the reaction.

A rule is a function that takes quantum numbers as input and outputs a boolean. There are three different types of rules:

GraphElementRule that work on individual graph edges or nodes.
EdgeQNConservationRule that work on the interaction level, which use ingoing edges, outgoing edges as arguments. E.g.: ChargeConservation.
ConservationRule that work on the interaction level, which use ingoing edges, outgoing edges and a interaction node as arguments. E.g: parity_conservation.

The arguments can be any type of quantum number. However a rule argument resembling edges only accepts EdgeQuantumNumbers. Similarly arguments that resemble a node only accept NodeQuantumNumbers. The argument types do not have to be limited to a single quantum number, but can be a composite (see CParityEdgeInput).

Warning

Besides the rule logic itself, a rule also has the responsibility of stating its run conditions. These run conditions must be stated by the type annotations of its __call__ method. The type annotations therefore are not just there for static type checking: they also carry more information about the rule that is extracted dynamically by the solving module.

Generally, the conditions can be separated into two categories:

variable conditions
toplogical conditions

Currently, only variable conditions are being used. Topological conditions could be created in the form of Tuple instead of List.

For additive quantum numbers, the decorator additive_quantum_number_rule can be used to automatically generate the appropriate behavior.

The module is therefore strongly typed (both for the reader of the code and for type checking with mypy). An example is HelicityParityEdgeInput, which has been defined to provide type checks on parity_conservation_helicity.

class BaryonNumberConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for baryon_number conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number]) → bool ¶: Call self as a function.

class BottomnessConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for bottomness conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.bottomness], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.bottomness]) → bool ¶: Call self as a function.

class CParityEdgeInput(spin_magnitude, pid, c_parity=None)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class CParityEdgeInput.

c_parity: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.c_parity]¶

pid: qrules.quantum_numbers.EdgeQuantumNumbers.pid¶

spin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude¶

class CParityNodeInput(l_magnitude, s_magnitude)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class CParityNodeInput.

l_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude¶

s_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude¶

class ChargeConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for charge conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.charge], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.charge]) → bool ¶: Call self as a function.

class CharmConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for charmness conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.charmness], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.charmness]) → bool ¶: Call self as a function.

class ConservationRule(*args, **kwargs)[source]¶

Bases: Protocol

__call__(_ConservationRule__ingoing_edge_qns: List[Any], _ConservationRule__outgoing_edge_qns: List[Any], _ConservationRule__node_qns: Any) → bool [source]¶: Call self as a function.

class EdgeQNConservationRule(*args, **kwargs)[source]¶

Bases: Protocol

__call__(_EdgeQNConservationRule__ingoing_edge_qns: List[Any], _EdgeQNConservationRule__outgoing_edge_qns: List[Any]) → bool [source]¶: Call self as a function.

class ElectronLNConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for electron_lepton_number conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number]) → bool ¶: Call self as a function.

class GParityEdgeInput(isospin_magnitude, spin_magnitude, pid, g_parity=None)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class GParityEdgeInput.

g_parity: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.g_parity]¶

isospin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.isospin_magnitude¶

pid: qrules.quantum_numbers.EdgeQuantumNumbers.pid¶

spin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude¶

class GParityNodeInput(l_magnitude, s_magnitude)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class GParityNodeInput.

l_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude¶

s_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude¶

class GellMannNishijimaInput(charge, isospin_projection=None, strangeness=None, charmness=None, bottomness=None, topness=None, baryon_number=None, electron_lepton_number=None, muon_lepton_number=None, tau_lepton_number=None)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class GellMannNishijimaInput.

baryon_number: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number]¶

bottomness: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.bottomness]¶

charge: qrules.quantum_numbers.EdgeQuantumNumbers.charge¶

charmness: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.charmness]¶

electron_lepton_number: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number]¶

isospin_projection: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.isospin_projection]¶

muon_lepton_number: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number]¶

strangeness: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.strangeness]¶

tau_lepton_number: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number]¶

topness: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.topness]¶

class GraphElementRule(*args, **kwargs)[source]¶

Bases: Protocol

__call__(_GraphElementRule__qns: Any) → bool [source]¶: Call self as a function.

class HelicityParityEdgeInput(parity, spin_magnitude, spin_projection)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class HelicityParityEdgeInput.

parity: qrules.quantum_numbers.EdgeQuantumNumbers.parity¶

spin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude¶

spin_projection: qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection¶

class IdenticalParticleSymmetryOutEdgeInput(spin_magnitude, spin_projection, pid)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class IdenticalParticleSymmetryOutEdgeInput.

pid: qrules.quantum_numbers.EdgeQuantumNumbers.pid¶

spin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude¶

spin_projection: qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection¶

class IsoSpinEdgeInput(isospin_magnitude, isospin_projection)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class IsoSpinEdgeInput.

isospin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.isospin_magnitude¶

isospin_projection: qrules.quantum_numbers.EdgeQuantumNumbers.isospin_projection¶

class MassConservation(width_factor: float)[source]¶

Bases: object

Mass conservation rule.

__call__(ingoing_edge_qns: List[qrules.conservation_rules.MassEdgeInput], outgoing_edge_qns: List[qrules.conservation_rules.MassEdgeInput]) → bool [source]¶

Implements mass conservation.

\(M_{out} - N \cdot W_{out} < M_{in} + N \cdot W_{in}\)

It makes sure that the net mass outgoing state \(M_{out}\) is smaller than the net mass of the ingoing state \(M_{in}\). Also the width \(W\) of the states is taken into account.

class MassEdgeInput(mass, width=None)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class MassEdgeInput.

mass: qrules.quantum_numbers.EdgeQuantumNumbers.mass¶

width: Optional[qrules.quantum_numbers.EdgeQuantumNumbers.width]¶

class MuonLNConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for muon_lepton_number conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number]) → bool ¶: Call self as a function.

class SpinEdgeInput(spin_magnitude, spin_projection)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class SpinEdgeInput.

spin_magnitude: qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude¶

spin_projection: qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection¶

class SpinMagnitudeNodeInput(l_magnitude, s_magnitude)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class SpinMagnitudeNodeInput.

l_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude¶

s_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude¶

class SpinNodeInput(l_magnitude, l_projection, s_magnitude, s_projection)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class SpinNodeInput.

l_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude¶

l_projection: qrules.quantum_numbers.NodeQuantumNumbers.l_projection¶

s_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude¶

s_projection: qrules.quantum_numbers.NodeQuantumNumbers.s_projection¶

class StrangenessConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for strangeness conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.strangeness], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.strangeness]) → bool ¶: Call self as a function.

class TauLNConservation(*args, **kwargs)[source]¶

Bases: qrules.conservation_rules.EdgeQNConservationRule

Decorated via additive_quantum_number_rule.

Check for tau_lepton_number conservation.

__call__(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number]) → bool ¶: Call self as a function.

additive_quantum_number_rule(quantum_number: type) → Callable[[Any], qrules.conservation_rules.EdgeQNConservationRule][source]¶

Class decorator for creating an additive conservation rule.

Use this decorator to create a EdgeQNConservationRule for a quantum number to which an additive conservation rule applies:

\[\sum q_{in} = \sum q_{out}\]

Parameters: quantum_number – Quantum number to which you want to apply the additive conservation check. An example would be EdgeQuantumNumbers.charge.

c_parity_conservation(ingoing_edge_qns: List[qrules.conservation_rules.CParityEdgeInput], outgoing_edge_qns: List[qrules.conservation_rules.CParityEdgeInput], interaction_node_qns: qrules.conservation_rules.CParityNodeInput) → bool [source]¶

Check for \(C\)-parity conservation.

Implements \(C_{in} = C_{out}\).

clebsch_gordan_helicity_to_canonical(ingoing_spins: List[qrules.conservation_rules.SpinEdgeInput], outgoing_spins: List[qrules.conservation_rules.SpinEdgeInput], interaction_qns: qrules.conservation_rules.SpinNodeInput) → bool [source]¶

Implement Clebsch-Gordan checks.

For \(S_1, S_2\) to \(S\) and the \(L,S\) to \(J\) coupling based on the conversion of helicity to canonical amplitude sums.

Note

This rule does not check that the spin magnitudes couple correctly to L and S, as this is already performed by spin_magnitude_conservation.

g_parity_conservation(ingoing_edge_qns: List[qrules.conservation_rules.GParityEdgeInput], outgoing_edge_qns: List[qrules.conservation_rules.GParityEdgeInput], interaction_qns: qrules.conservation_rules.GParityNodeInput) → bool [source]¶

Check for \(G\)-parity conservation.

Implements for \(G_{in} = G_{out}\).

gellmann_nishijima(edge_qns: qrules.conservation_rules.GellMannNishijimaInput) → bool [source]¶

Check the Gell-Mann–Nishijima formula.

Gell-Mann–Nishijima formula:

\[Q = I_3 + \frac{1}{2}(B+S+C+B'+T)\]

where \(Q\) is charge (computed), \(I_3\) is Spin.projection of isospin, \(B\) is baryon_number, \(S\) is strangeness, \(C\) is charmness, \(B'\) is bottomness, and \(T\) is topness.

helicity_conservation(ingoing_spin_mags: List[qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude], outgoing_helicities: List[qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection]) → bool [source]¶

Implementation of helicity conservation.

Check for \(|\lambda_2-\lambda_3| \leq S_1\).

identical_particle_symmetrization(ingoing_parities: List[qrules.quantum_numbers.EdgeQuantumNumbers.parity], outgoing_edge_qns: List[qrules.conservation_rules.IdenticalParticleSymmetryOutEdgeInput]) → bool [source]¶

Verifies multi particle state symmetrization for identical particles.

In case of a multi particle state with identical particles, their exchange symmetry has to follow the spin statistic theorem.

For bosonic systems the total exchange symmetry (parity) has to be even (+1). For fermionic systems the total exchange symmetry (parity) has to be odd (-1).

In case of a particle decaying into N identical particles (N>1), the decaying particle has to have the same parity as required by the spin statistic theorem of the multi body state.

isospin_conservation(ingoing_isospins: List[qrules.conservation_rules.IsoSpinEdgeInput], outgoing_isospins: List[qrules.conservation_rules.IsoSpinEdgeInput]) → bool [source]¶

Check for isospin conservation.

Implements

\[|I_1 - I_2| \leq I \leq |I_1 + I_2|\]

Also checks \(I_{1,z} + I_{2,z} = I_z\) and if Clebsch-Gordan coefficients are all 0.

isospin_validity(isospin: qrules.conservation_rules.IsoSpinEdgeInput) → bool [source]¶: Check for valid isospin magnitude and projection.

ls_spin_validity(spin_input: qrules.conservation_rules.SpinNodeInput) → bool [source]¶: Check for valid isospin magnitude and projection.

parity_conservation(ingoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.parity], outgoing_edge_qns: List[qrules.quantum_numbers.EdgeQuantumNumbers.parity], l_magnitude: qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude) → bool [source]¶: Implement \(P_{in} = P_{out} \cdot (-1)^L\).

parity_conservation_helicity(ingoing_edge_qns: List[qrules.conservation_rules.HelicityParityEdgeInput], outgoing_edge_qns: List[qrules.conservation_rules.HelicityParityEdgeInput], parity_prefactor: qrules.quantum_numbers.NodeQuantumNumbers.parity_prefactor) → bool [source]¶

Implements parity conservation for helicity formalism.

Check the following:

\[A_{-\lambda_1-\lambda_2} = P_1 P_2 P_3 (-1)^{S_2+S_3-S_1} A_{\lambda_1\lambda_2}\]

\[\mathrm{parity\,prefactor} = P_1 P_2 P_3 (-1)^{S_2+S_3-S_1}\]

Note

Only the special case \(\lambda_1=\lambda_2=0\) may return False independent on the parity prefactor.

spin_conservation(ingoing_spins: List[qrules.conservation_rules.SpinEdgeInput], outgoing_spins: List[qrules.conservation_rules.SpinEdgeInput], interaction_qns: qrules.conservation_rules.SpinNodeInput) → bool [source]¶

Check for spin conservation.

Implements

\[|S_1 - S_2| \leq S \leq |S_1 + S_2|\]

and

\[|L - S| \leq J \leq |L + S|\]

Also checks \(M_1 + M_2 = M\) and if Clebsch-Gordan coefficients are all 0.

spin_magnitude_conservation(ingoing_spins: List[qrules.conservation_rules.SpinEdgeInput], outgoing_spins: List[qrules.conservation_rules.SpinEdgeInput], interaction_qns: qrules.conservation_rules.SpinMagnitudeNodeInput) → bool [source]¶

Check for spin conservation.

Implements

\[|S_1 - S_2| \leq S \leq |S_1 + S_2|\]

and

\[|L - S| \leq J \leq |L + S|\]

spin_validity(spin: qrules.conservation_rules.SpinEdgeInput) → bool [source]¶: Check for valid spin magnitude and projection.

default_settings¶

import qrules.default_settings

Default configuration for the qrules.

class InteractionTypes(value)[source]¶

Bases: enum.Enum

Types of interactions in the form of an enumerate.

EM = 2¶

STRONG = 1¶

WEAK = 3¶

create_default_interaction_settings(formalism_type: str, nbody_topology: bool = False, mass_conservation_factor: Optional[float] = 3.0) → Dict[qrules.default_settings.InteractionTypes, Tuple[qrules.solving.EdgeSettings, qrules.solving.NodeSettings]][source]¶

Create a container that holds the settings for the various interactions.

E.g.: strong, em and weak interaction.

particle¶

import qrules.particle

A collection of particle info containers.

The particle module is the starting point of qrules. Its main interface is the ParticleCollection, which is a collection of immutable Particle instances that are uniquely defined by their properties. As such, it can be used stand-alone as a database of quantum numbers (see Particle database).

The transition module uses the properties of Particle instances when it computes which StateTransitionGraph s are allowed between an initial state and final state.

class Particle(*, name: str, pid: int, latex: Optional[str] = None, spin, mass, width=0.0, charge: int = 0, isospin=None, strangeness: int = 0, charmness: int = 0, bottomness: int = 0, topness: int = 0, baryon_number: int = 0, electron_lepton_number: int = 0, muon_lepton_number: int = 0, tau_lepton_number: int = 0, parity=None, c_parity=None, g_parity=None)[source]¶

Bases: object

Immutable container of data defining a physical particle.

A Particle is defined by the minimum set of the quantum numbers that every possible instances of that particle have in common (the “static” quantum numbers of the particle). A “non-static” quantum number is the spin projection. Hence Particle instances do not contain spin projection information.

Particle instances are uniquely defined by their quantum numbers and properties like mass. The name and pid are therefore just labels that are not taken into account when checking if two Particle instances are equal.

Note

As opposed to classes such as EdgeQuantumNumbers and NodeQuantumNumbers, the Particle class serves as an interface to the user (see Particle database).

__eq__(other)¶: Method generated by attrs for class Particle.

baryon_number: int¶

bottomness: int¶

c_parity: Optional[qrules.quantum_numbers.Parity]¶

charge: int¶

charmness: int¶

electron_lepton_number: int¶

g_parity: Optional[qrules.quantum_numbers.Parity]¶

is_lepton() → bool [source]¶

isospin: Optional[qrules.particle.Spin]¶

latex: Optional[str]¶

mass: float¶

muon_lepton_number: int¶

name: str¶

parity: Optional[qrules.quantum_numbers.Parity]¶

pid: int¶

spin: float¶

strangeness: int¶

tau_lepton_number: int¶

topness: int¶

width: float¶

class ParticleCollection(particles: Optional[Iterable[qrules.particle.Particle]] = None)[source]¶

Bases: collections.abc.MutableSet

Searchable collection of immutable Particle instances.

__eq__(other: object) → bool [source]¶: Return self==value.

add(value: qrules.particle.Particle) → None [source]¶: Add an element.

discard(value: Union[qrules.particle.Particle, str]) → None [source]¶: Remove an element. Do not raise an exception if absent.

filter(function: Callable[[qrules.particle.Particle], bool]) → qrules.particle.ParticleCollection [source]¶

Search by Particle properties using a lambda function.

For example:

>>> from qrules.particle import load_pdg
>>> pdg = load_pdg()
>>> subset = pdg.filter(
...     lambda p: p.mass > 1.8
...     and p.mass < 2.0
...     and p.spin == 2
...     and p.strangeness == 1
... )
>>> sorted(list(subset.names))
['K(2)(1820)+', 'K(2)(1820)0']

find(search_term: Union[int, str]) → qrules.particle.Particle [source]¶: Search for a particle by either name (str) or PID (int).

property names¶

update(other: Iterable[qrules.particle.Particle]) → None [source]¶

class Spin(magnitude, projection)[source]¶

Bases: object

Safe, immutable data container for spin with projection.

__eq__(other: object) → bool [source]¶: Return self==value.

magnitude: float¶

projection: float¶

create_antiparticle(template_particle: qrules.particle.Particle, new_name: Optional[str] = None, new_latex: Optional[str] = None) → qrules.particle.Particle [source]¶

create_particle(template_particle: qrules.particle.Particle, name: Optional[str] = None, latex: Optional[str] = None, pid: Optional[int] = None, mass: Optional[float] = None, width: Optional[float] = None, charge: Optional[int] = None, spin: Optional[float] = None, isospin: Optional[qrules.particle.Spin] = None, strangeness: Optional[int] = None, charmness: Optional[int] = None, bottomness: Optional[int] = None, topness: Optional[int] = None, baryon_number: Optional[int] = None, electron_lepton_number: Optional[int] = None, muon_lepton_number: Optional[int] = None, tau_lepton_number: Optional[int] = None, parity: Optional[int] = None, c_parity: Optional[int] = None, g_parity: Optional[int] = None) → qrules.particle.Particle [source]¶

load_pdg() → qrules.particle.ParticleCollection [source]¶

Create a ParticleCollection with all entries from the PDG.

PDG info is imported from the scikit-hep/particle package.

quantum_numbers¶

import qrules.quantum_numbers

Definitions used internally for type hints and signatures.

qrules is strictly typed (enforced through mypy). This module bundles structures and definitions that don’t serve as data containers but only as type hints. EdgeQuantumNumbers and NodeQuantumNumbers are the main structures and serve as a bridge between the particle and the conservation_rules module.

class EdgeQuantumNumbers[source]¶

Bases: object

Definition of quantum numbers for edges.

This class defines the types that are used in the conservation_rules, for instance in additive_quantum_number_rule. You can also create data classes (see attr.s) with data members that are typed as the data members of EdgeQuantumNumbers (see for example HelicityParityEdgeInput) and use them in conservation rules that satisfy the appropriate rule protocol (see ConservationRule, EdgeQNConservationRule).

__eq__(other)¶: Method generated by attrs for class EdgeQuantumNumbers.

baryon_number()¶

bottomness()¶

c_parity()¶

charge()¶

charmness()¶

electron_lepton_number()¶

g_parity()¶

isospin_magnitude()¶

isospin_projection()¶

mass()¶

muon_lepton_number()¶

parity()¶

pid()¶

spin_magnitude()¶

spin_projection()¶

strangeness()¶

tau_lepton_number()¶

topness()¶

width()¶

class InteractionProperties(l_magnitude=None, l_projection=None, s_magnitude=None, s_projection=None, parity_prefactor=None)[source]¶

Bases: object

Immutable data structure containing interaction properties.

Note

As opposed to NodeQuantumNumbers, the InteractionProperties class serves as an interface to the user.

__eq__(other)¶: Method generated by attrs for class InteractionProperties.

l_magnitude: Optional[int]¶

l_projection: Optional[int]¶

parity_prefactor: Optional[float]¶

s_magnitude: Optional[float]¶

s_projection: Optional[float]¶

class NodeQuantumNumbers[source]¶

Bases: object

Definition of quantum numbers for interaction nodes.

__eq__(other)¶: Method generated by attrs for class NodeQuantumNumbers.

l_magnitude()¶

l_projection()¶

parity_prefactor()¶

s_magnitude()¶

s_projection()¶

class Parity(value: int)[source]¶

Bases: object

__eq__(other: object) → bool [source]¶: Return self==value.

value: int¶

arange(x_1: float, x_2: float, delta: float) → Generator[float, None, None][source]¶

edge_qn_type(self)¶: Method generated by attrs for class EdgeQuantumNumbers.

node_qn_type(self)¶: Method generated by attrs for class NodeQuantumNumbers.

solving¶

import qrules.solving

Functions to solve a particle reaction problem.

This module is responsible for solving a particle reaction problem stated by a StateTransitionGraph and corresponding GraphSettings. The Solver classes (e.g. CSPSolver) generate new quantum numbers (for example belonging to an intermediate state) and validate the decay processes with the rules formulated by the conservation_rules module.

Bases: qrules.solving.Solver

Solver reducing the task to a Constraint Satisfaction Problem.

Solving this done with the python-constraint module.

The variables are the quantum numbers of particles/edges, but also some composite quantum numbers which are attributed to the interaction nodes (such as angular momentum \(L\)). The conservation rules serve as the constraints and a special wrapper class serves as an adapter.

find_solutions(problem_set: qrules.solving.QNProblemSet) → qrules.solving.QNResult [source]¶

Find solutions for the given input.

It is expected that this function determines and returns all of the found solutions. In case no solutions are found a partial list of violated rules has to be given. This list of violated rules does not have to be complete.

Parameters: problem_set (QNProblemSet) – states a problem set
Returns: contains possible solutions, violated rules and not executed rules due to requirement issues.
Return type: QNResult

class EdgeSettings(conservation_rules: Set[qrules.conservation_rules.GraphElementRule] = NOTHING, rule_priorities: Dict[qrules.conservation_rules.GraphElementRule, int] = NOTHING, qn_domains: Dict[Any, Any] = NOTHING)[source]¶

Bases: object

Solver settings for a specific edge of a graph.

__eq__(other)¶: Method generated by attrs for class EdgeSettings.

conservation_rules: Set[qrules.conservation_rules.GraphElementRule]¶

qn_domains: Dict[Any, Any]¶

rule_priorities: Dict[qrules.conservation_rules.GraphElementRule, int]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class GraphElementProperties.

edge_props: Dict[int, Dict[Type[Union[qrules.quantum_numbers.EdgeQuantumNumbers.pid, qrules.quantum_numbers.EdgeQuantumNumbers.mass, qrules.quantum_numbers.EdgeQuantumNumbers.width, qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.charge, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.strangeness, qrules.quantum_numbers.EdgeQuantumNumbers.charmness, qrules.quantum_numbers.EdgeQuantumNumbers.bottomness, qrules.quantum_numbers.EdgeQuantumNumbers.topness, qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number, qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.parity, qrules.quantum_numbers.EdgeQuantumNumbers.c_parity, qrules.quantum_numbers.EdgeQuantumNumbers.g_parity]], Union[int, float]]]¶

node_props: Dict[int, Dict[Type[Union[qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.l_projection, qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.s_projection, qrules.quantum_numbers.NodeQuantumNumbers.parity_prefactor]], Union[int, float]]]¶

class GraphSettings(edge_settings: Dict[int, qrules.solving.EdgeSettings] = NOTHING, node_settings: Dict[int, qrules.solving.NodeSettings] = NOTHING)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class GraphSettings.

edge_settings: Dict[int, qrules.solving.EdgeSettings]¶

node_settings: Dict[int, qrules.solving.NodeSettings]¶

class NodeSettings(conservation_rules: Set[Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule]] = NOTHING, rule_priorities: Dict[Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule], int] = NOTHING, qn_domains: Dict[Any, Any] = NOTHING)[source]¶

Bases: object

Container class for the interaction settings.

This class can be assigned to each node of a state transition graph. Hence, these settings contain the complete configuration information which is required for the solution finding, e.g:

set of conservation rules

mapping of rules to priorities (optional)

mapping of quantum numbers to their domains

strength scale parameter (higher value means stronger force)

__eq__(other)¶: Method generated by attrs for class NodeSettings.

conservation_rules: Set[Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule]]¶

interaction_strength: float = 1.0¶

qn_domains: Dict[Any, Any]¶

rule_priorities: Dict[Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule], int]¶

class QNProblemSet(topology: qrules.topology.Topology, initial_facts: qrules.solving.GraphElementProperties, solving_settings: qrules.solving.GraphSettings)[source]¶

Bases: object

Particle reaction problem set, defined as a graph like data structure.

Parameters

topology (Topology) – a topology that represent the structure of the reaction
initial_facts (GraphElementProperties) – all of the known facts quantum numbers of the problem
solving_settings (GraphSettings) – solving specific settings such as the specific rules and variable domains for nodes and edges of the topology

__eq__(other)¶: Method generated by attrs for class QNProblemSet.

initial_facts: qrules.solving.GraphElementProperties¶

solving_settings: qrules.solving.GraphSettings¶

topology: qrules.topology.Topology¶

class QNResult(solutions: List[qrules.solving.QuantumNumberSolution] = NOTHING, not_executed_node_rules: Dict[int, Set[str]] = NOTHING, violated_node_rules: Dict[int, Set[str]] = NOTHING, not_executed_edge_rules: Dict[int, Set[str]] = NOTHING, violated_edge_rules: Dict[int, Set[str]] = NOTHING)[source]¶

Bases: object

Defines a result to a problem set processed by the solving code.

__eq__(other)¶: Method generated by attrs for class QNResult.

extend(other_result: qrules.solving.QNResult) → None [source]¶

not_executed_edge_rules: Dict[int, Set[str]]¶

not_executed_node_rules: Dict[int, Set[str]]¶

solutions: List[qrules.solving.QuantumNumberSolution]¶

violated_edge_rules: Dict[int, Set[str]]¶

violated_node_rules: Dict[int, Set[str]]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class QuantumNumberSolution.

edge_quantum_numbers: Dict[int, Dict[Type[Union[qrules.quantum_numbers.EdgeQuantumNumbers.pid, qrules.quantum_numbers.EdgeQuantumNumbers.mass, qrules.quantum_numbers.EdgeQuantumNumbers.width, qrules.quantum_numbers.EdgeQuantumNumbers.spin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.spin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.charge, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_magnitude, qrules.quantum_numbers.EdgeQuantumNumbers.isospin_projection, qrules.quantum_numbers.EdgeQuantumNumbers.strangeness, qrules.quantum_numbers.EdgeQuantumNumbers.charmness, qrules.quantum_numbers.EdgeQuantumNumbers.bottomness, qrules.quantum_numbers.EdgeQuantumNumbers.topness, qrules.quantum_numbers.EdgeQuantumNumbers.baryon_number, qrules.quantum_numbers.EdgeQuantumNumbers.electron_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.muon_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.tau_lepton_number, qrules.quantum_numbers.EdgeQuantumNumbers.parity, qrules.quantum_numbers.EdgeQuantumNumbers.c_parity, qrules.quantum_numbers.EdgeQuantumNumbers.g_parity]], Union[int, float]]]¶

node_quantum_numbers: Dict[int, Dict[Type[Union[qrules.quantum_numbers.NodeQuantumNumbers.l_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.l_projection, qrules.quantum_numbers.NodeQuantumNumbers.s_magnitude, qrules.quantum_numbers.NodeQuantumNumbers.s_projection, qrules.quantum_numbers.NodeQuantumNumbers.parity_prefactor]], Union[int, float]]]¶

class Scoresheet[source]¶

Bases: object

register_rule(graph_element_id: int, rule: Union[qrules.conservation_rules.GraphElementRule, qrules.conservation_rules.EdgeQNConservationRule, qrules.conservation_rules.ConservationRule]) → Callable[[bool], None][source]¶

property rule_calls¶

property rule_passes¶

class Solver[source]¶

Bases: abc.ABC

Interface of a Solver.

abstract find_solutions(problem_set: qrules.solving.QNProblemSet) → qrules.solving.QNResult [source]¶

Find solutions for the given input.

It is expected that this function determines and returns all of the found solutions. In case no solutions are found a partial list of violated rules has to be given. This list of violated rules does not have to be complete.

Parameters: problem_set (QNProblemSet) – states a problem set
Returns: contains possible solutions, violated rules and not executed rules due to requirement issues.
Return type: QNResult

validate_full_solution(problem_set: qrules.solving.QNProblemSet) → qrules.solving.QNResult [source]¶

topology¶

import qrules.topology

All modules related to topology building.

Responsible for building all possible topologies bases on basic user input:

number of initial state particles
number of final state particles

The main interface is the StateTransitionGraph.

class Edge(originating_node_id=None, ending_node_id=None)[source]¶

Bases: object

Struct-like definition of an edge, used in Topology.

__eq__(other)¶: Method generated by attrs for class Edge.

ending_node_id: Optional[int]¶

get_connected_nodes() → Set[int][source]¶

originating_node_id: Optional[int]¶

EdgeType¶

A TypeVar representing the type of edge properties.

alias of TypeVar(‘EdgeType’)

class FrozenDict(mapping: Optional[Mapping] = None)[source]¶

Bases: Generic[qrules.topology._K, qrules.topology._V], collections.abc.Hashable, collections.abc.Mapping

items() → a set-like object providing a view on D’s items[source]¶

keys() → a set-like object providing a view on D’s keys[source]¶

values() → an object providing a view on D’s values[source]¶

class InteractionNode(number_of_ingoing_edges: int, number_of_outgoing_edges: int)[source]¶

Bases: object

Helper class for the SimpleStateTransitionTopologyBuilder.

__eq__(other)¶: Method generated by attrs for class InteractionNode.

number_of_ingoing_edges: int¶

number_of_outgoing_edges: int¶

class SimpleStateTransitionTopologyBuilder(interaction_node_set: Iterable[qrules.topology.InteractionNode])[source]¶

Bases: object

Simple topology builder.

Recursively tries to add the interaction nodes to available open end edges/lines in all combinations until the number of open end lines matches the final state lines.

build(number_of_initial_edges: int, number_of_final_edges: int) → Tuple[qrules.topology.Topology, …][source]¶

class StateTransitionGraph(topology: qrules.topology.Topology, node_props: Mapping[int, qrules.quantum_numbers.InteractionProperties], edge_props: Mapping[int, EdgeType])[source]¶

Bases: Generic[qrules.topology.EdgeType]

Graph class that resembles a frozen Topology with properties.

This class should contain the full information of a state transition from a initial state to a final state. This information can be attached to the nodes and edges via properties. In case not all information is provided, error can be raised on property retrieval.

__eq__(other: object) → bool [source]¶: Check if two StateTransitionGraph instances are identical.

compare(other: qrules.topology.StateTransitionGraph, edge_comparator: Optional[Callable[[EdgeType, EdgeType], bool]] = None, node_comparator: Optional[Callable[[qrules.quantum_numbers.InteractionProperties, qrules.quantum_numbers.InteractionProperties], bool]] = None) → bool [source]¶

evolve(node_props: Optional[Dict[int, qrules.quantum_numbers.InteractionProperties]] = None, edge_props: Optional[Dict[int, EdgeType]] = None) → qrules.topology.StateTransitionGraph[EdgeType][source]¶

Changes the node and edge properties of a graph instance.

Since a StateTransitionGraph is frozen (cannot be modified), the evolve function will also create a shallow copy the properties.

get_edge_props(edge_id: int) → EdgeType [source]¶

get_node_props(node_id: int) → qrules.quantum_numbers.InteractionProperties [source]¶

swap_edges(edge_id1: int, edge_id2: int) → None [source]¶

class Topology(nodes, edges)[source]¶

Bases: object

Directed Feynman-like graph without edge or node properties.

Forms the underlying topology of StateTransitionGraph. The graphs are directed, meaning the edges are ingoing and outgoing to specific nodes (since feynman graphs also have a time axis). Note that a Topology is not strictly speaking a graph from graph theory, because it allows open edges, like a Feynman-diagram.

__eq__(other)¶: Method generated by attrs for class Topology.

edges: qrules.topology.FrozenDict[int, qrules.topology.Edge]¶

get_edge_ids_ingoing_to_node(node_id: int) → Set[int][source]¶

get_edge_ids_outgoing_from_node(node_id: int) → Set[int][source]¶

get_originating_final_state_edge_ids(node_id: int) → Set[int][source]¶

get_originating_initial_state_edge_ids(node_id: int) → Set[int][source]¶

incoming_edge_ids: FrozenSet[int]¶

intermediate_edge_ids: FrozenSet[int]¶

is_isomorphic(other: qrules.topology.Topology) → bool [source]¶

Check if two graphs are isomorphic.

Returns: True if the two graphs have a one-to-one mapping of the node IDs and edge IDs.
Return type: bool

nodes: FrozenSet[int]¶

organize_edge_ids() → qrules.topology.Topology [source]¶

Create a new topology with edge IDs in range [-m, n+i].

where m is the number of incoming_edge_ids, n is the number of outgoing_edge_ids, and i is the number of intermediate_edge_ids.

In other words, relabel the edges so that:

incoming_edge_ids lies in the range [-1, -2, ...]
outgoing_edge_ids lies in the range [0, 1, ..., n]
intermediate_edge_ids lies in the range [n+1, n+2, ...]

outgoing_edge_ids: FrozenSet[int]¶

swap_edges(edge_id1: int, edge_id2: int) → qrules.topology.Topology [source]¶

create_isobar_topologies(number_of_final_states: int) → Tuple[qrules.topology.Topology, …][source]¶

create_n_body_topology(number_of_initial_states: int, number_of_final_states: int) → qrules.topology.Topology [source]¶

get_originating_node_list(topology: qrules.topology.Topology, edge_ids: Iterable[int]) → List[int][source]¶

Get list of node ids from which the supplied edges originate from.

Parameters: edge_ids ([int]) – list of edge ids for which the origin node is searched for
Returns: a list of node ids
Return type: [int]

transition¶

import qrules.transition

Find allowed transitions between an initial and final state.

class ExecutionInfo(not_executed_node_rules: Dict[int, Set[str]] = NOTHING, violated_node_rules: Dict[int, Set[str]] = NOTHING, not_executed_edge_rules: Dict[int, Set[str]] = NOTHING, violated_edge_rules: Dict[int, Set[str]] = NOTHING)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class ExecutionInfo.

clear() → None [source]¶

extend(other_result: qrules.transition.ExecutionInfo, intersect_violations: bool = False) → None [source]¶

not_executed_edge_rules: Dict[int, Set[str]]¶

not_executed_node_rules: Dict[int, Set[str]]¶

violated_edge_rules: Dict[int, Set[str]]¶

violated_node_rules: Dict[int, Set[str]]¶

class ProblemSet(topology: qrules.topology.Topology, initial_facts: qrules.combinatorics.InitialFacts, solving_settings: qrules.solving.GraphSettings)[source]¶

Bases: object

Particle reaction problem set, defined as a graph like data structure.

Parameters

topology – Topology that contains the structure of the reaction.
initial_facts – InitialFacts that contain the info of initial and final state in connection with the topology.
solving_settings – Solving related settings such as the conservation rules and the quantum number domains.

__eq__(other)¶: Method generated by attrs for class ProblemSet.

initial_facts: qrules.combinatorics.InitialFacts¶

solving_settings: qrules.solving.GraphSettings¶

to_qn_problem_set() → qrules.solving.QNProblemSet [source]¶

topology: qrules.topology.Topology¶

class Result(transitions: List[qrules.topology.StateTransitionGraph[Tuple[qrules.particle.Particle, float]]] = NOTHING, formalism_type: Optional[str] = None)[source]¶

Bases: object

__eq__(other)¶: Method generated by attrs for class Result.

formalism_type: Optional[str]¶

get_final_state() → List[qrules.particle.Particle][source]¶

get_initial_state() → List[qrules.particle.Particle][source]¶

get_intermediate_particles() → qrules.particle.ParticleCollection [source]¶: Extract the names of the intermediate state particles.

transitions: List[qrules.topology.StateTransitionGraph[Tuple[qrules.particle.Particle, float]]]¶

class SolvingMode(value)[source]¶

Bases: enum.Enum

Types of modes for solving.

FAST = 1¶: Find “likeliest” solutions only.

FULL = 2¶: Find all possible solutions.

class StateTransitionManager(initial_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], final_state: Sequence[Union[str, Tuple[str, Sequence[float]]]], particles: Optional[qrules.particle.ParticleCollection] = None, allowed_intermediate_particles: Optional[List[str]] = None, interaction_type_settings: Optional[Dict[qrules.default_settings.InteractionTypes, Tuple[qrules.solving.EdgeSettings, qrules.solving.NodeSettings]]] = None, formalism_type: str = 'helicity', topology_building: str = 'isobar', number_of_threads: Optional[int] = None, solving_mode: qrules.transition.SolvingMode = <SolvingMode.FAST: 1>, reload_pdg: bool = False, mass_conservation_factor: Optional[float] = 3.0)[source]¶

Bases: object

Main handler for decay topologies.

add_final_state_grouping(fs_group: List[Union[str, List[str]]]) → None [source]¶

create_problem_sets() → Dict[float, List[qrules.transition.ProblemSet]][source]¶

find_solutions(problem_sets: Dict[float, List[qrules.transition.ProblemSet]]) → qrules.transition.Result [source]¶: Check for solutions for a specific set of interaction settings.

property formalism_type¶

set_allowed_interaction_types(allowed_interaction_types: List[qrules.default_settings.InteractionTypes]) → None [source]¶

set_allowed_intermediate_particles(particle_names: List[str]) → None [source]¶

Architectural Decision Records¶

This log lists the architectural decisions for the QRules repository:

[ADR-000] Use ADRs¶

Status: accepted

Deciders: @redeboer, @spflueger

Technical story: A large number of issues in the expertsystem are to be correlated (e.g. #40, #44, #22) so that resulting PRs (in this case, #42) lacked direction. This led us to consider ADRs.

Context and Problem Statement¶

We want to record architectural decisions made in this project. Which format and structure should these records follow?

Considered Options¶

MADR 2.1.2 – The Markdown Architectural Decision Records
Michael Nygard’s template – The first incarnation of the term “ADR”
Sustainable Architectural Decisions – The Y-Statements
Other templates listed at github.com/joelparkerhenderson/architecture_decision_record.
Formless – No conventions for file format and structure

Decision Outcome¶

Chosen option: “MADR 2.1.2”, because

Implicit assumptions should be made explicit. Design documentation is important to enable people understanding the decisions later on. See also A rational design process: How and why to fake it.
The MADR format is lean and fits our development style.
The MADR structure is comprehensible and facilitates usage & maintenance.
The MADR project is vivid.
Version 2.1.2 is the latest one available when starting to document ADRs.

[ADR-001] Amplitude model¶

Status: accepted
Deciders: @redeboer @spflueger

Context and problem statement¶

From the perspective of a PWA fitter package, the responsibility of the expertsystem is to construct a AmplitudeModel that serves as blueprint for a function that can be evaluated. Such a function has the following requirements:

It should be able to compute a list of real-valued intensities \(\mathbb{R}^m\) from a dataset of four-momenta \(\mathbb{R}^{m\times n\times4}\), where \(m\) is the number of events and \(n\) is the number of final state particles.
It should contain parameters that can be tweaked, so that they can be optimized with regard to a certain estimator.

Technical story¶

#382: Coupling parameters in the AmplitudeModel is difficult (has to be done through the place where they are used in the dynamics or intensity section) and counter-intuitive (cannot be done through the parameters section)
#440: when overwriting existing dynamics, old parameters are not cleaned up from the parameters section
#441: parameters contain a name that can be changed, but that results in a mismatch between the key that is used in the parameters section and the name of the parameter to which that entry points.
ComPWA#226: Use a math language for the blueprint of the function. This was also discussed early to mid 2020, but dropped in favor of custom python code + amplitf. The reasoning was that the effort of writing some new math language plus generators converting a mathematical expression into a function (using various back-ends) requires too much manpower.

Decision drivers¶

Solution requirements¶

The AmplitudeModel has to be convertible to a function which can be evaluated using various computation back-ends (numpy, tensorflow, theano, jax, …)
Ideally, the model should be complete in the sense that it contains all information to construct the complete model. This means that some “common” functions like a Breit-Wigner and Blatt-Weisskopf form factors should also be contained inside the AmplitudeModel. This guarantees reproducibility!
Adding new operators/models should not trigger many code modifications (open-closed principle), for instance adding new dynamics or formalisms.
Extendible:
- Add or replace current parts of an existing model. For example replace the dynamics part of some decay.
- Change a function plus a dataset to an estimator function. This is a subtle but important point. The function should hide its details (which backend and its mathematical expression) and yet be extendable to an estimator.
Definition and easy extraction of components. Components are certain sub-parts of the complete mathematical expression. This is at least needed for the calculation of fit fractions, or plotting individual parts of the intensity.

Considered solutions¶

Customized Python classes¶

Currently (v0.6.8), the AmplitudeModel contains five sections (instances of specific classes):

kinematics: defines initial and final state
particles: particle definitions (spin, etc.)
dynamics: a mapping that defines which dynamics type to apply to which particle
intensity: the actual amplitude model that is to be converted by a fitter package into a function as described above
parameters: an inventory of parameters that are used in intensity and dynamics

This structure can be represented in YAML, see an example here.

A fitter package converts intensity together with dynamics into a function. Any references to parameters that intensity or dynamics contain are converted into a parameter of the function. The parameters are initialized with the value as listed in the parameters section of the AmplitudeModel.

Alternative solutions¶

SymPy¶

Some useful SymPy pages:

Defining the amplitude model in terms of SymPY¶

Parameters would become a mapping of Symbols to initial, suggested values and dynamics would be a mapping of Symbols to ‘suggested’ expressions. Intensity will be the eventual combined expression.

# !pip install matplotlib==3.3.3 sympy==1.7.1

from typing import Dict

import attr
import sympy as sp


@attr.s
class AmplitudeModel:
    initial_values: Dict[sp.Symbol, float] = attr.ib(default=dict())
    dynamics: Dict[sp.Symbol, sp.Function] = attr.ib(default=dict())
    intensity: sp.Expr = attr.ib(default=None)

There needs to be one symbol \(x\) that represents the four-momentum input:

x = sp.Symbol("x")

As an example, let’s create an AmplitudeModel with an intensity that is a sum of Gaussians. Each Gaussian here takes the rôle of a dynamics function:

model = AmplitudeModel()

N_COMPONENTS = 3
for i in range(1, N_COMPONENTS + 1):
    mu = sp.Symbol(fR"\mu_{i}")
    sigma = sp.Symbol(fR"\sigma_{i}")
    model.initial_values.update(
        {
            mu: float(i),
            sigma: 1 / (2 * i),
        }
    )
    gauss = sp.exp(-((x - mu) ** 2) / (sigma ** 2)) / (
        sigma * sp.sqrt(2 * sp.pi)
    )
    dyn_symbol = sp.Symbol(fR"\mathrm{{dyn}}_{i}")
    model.dynamics[dyn_symbol] = gauss

coherent_sum = sum(model.dynamics)
model.intensity = coherent_sum

model.initial_values

{\mu_1: 1.0,
 \sigma_1: 0.5,
 \mu_2: 2.0,
 \sigma_2: 0.25,
 \mu_3: 3.0,
 \sigma_3: 0.16666666666666666}

model.intensity

\[\displaystyle \mathrm{dyn}_1 + \mathrm{dyn}_2 + \mathrm{dyn}_3\]

Dynamics are inserted into the intensity expression of the model:

model.intensity.subs(model.dynamics)

\[\displaystyle \frac{\sqrt{2} e^{- \frac{\left(- \mu_{3} + x\right)^{2}}{\sigma_{3}^{2}}}}{2 \sqrt{\pi} \sigma_{3}} + \frac{\sqrt{2} e^{- \frac{\left(- \mu_{2} + x\right)^{2}}{\sigma_{2}^{2}}}}{2 \sqrt{\pi} \sigma_{2}} + \frac{\sqrt{2} e^{- \frac{\left(- \mu_{1} + x\right)^{2}}{\sigma_{1}^{2}}}}{2 \sqrt{\pi} \sigma_{1}}\]

And, for evaluating, the ‘suggested’ initial parameter values are inserted:

model.intensity.subs(model.dynamics).subs(model.initial_values)

\[\displaystyle \frac{1.0 \sqrt{2} e^{- 4.0 \left(x - 1.0\right)^{2}}}{\sqrt{\pi}} + \frac{2.0 \sqrt{2} e^{- 64.0 \left(0.5 x - 1\right)^{2}}}{\sqrt{\pi}} + \frac{3.0 \sqrt{2} e^{- 324.0 \left(0.333333333333333 x - 1\right)^{2}}}{\sqrt{\pi}}\]

Here’s a small helper function to plot this model:

def plot_model(model: AmplitudeModel) -> None:
    total_plot = sp.plotting.plot(
        model.intensity.subs(model.dynamics).subs(model.initial_values),
        (x, 0, 4),
        show=False,
        line_color="black",
    )
    p1 = sp.plotting.plot(
        model.dynamics[sp.Symbol(R"\mathrm{dyn}_1")].subs(
            model.initial_values
        ),
        (x, 0, 4),
        line_color="red",
        show=False,
    )
    p2 = sp.plotting.plot(
        model.dynamics[sp.Symbol(R"\mathrm{dyn}_2")].subs(
            model.initial_values
        ),
        (x, 0, 4),
        line_color="blue",
        show=False,
    )
    p3 = sp.plotting.plot(
        model.dynamics[sp.Symbol(R"\mathrm{dyn}_3")].subs(
            model.initial_values
        ),
        (x, 0, 4),
        line_color="green",
        show=False,
    )
    total_plot.extend(p1)
    total_plot.extend(p2)
    total_plot.extend(p3)
    total_plot.show()

plot_model(model)

Now we can couple parameters like this:

model.initial_values[sp.Symbol(R"\sigma_1")] = sp.Symbol(R"\sigma_3")
plot_model(model)
model.initial_values[sp.Symbol(R"\sigma_3")] = 1
plot_model(model)

And it’s also possible to insert custom dynamics:

model.dynamics[sp.Symbol(R"\mathrm{dyn}_3")] = sp.sqrt(x)
plot_model(model)

Implementation in TensorWaves¶

Credits @spflueger

# !pip install jax==0.2.8 jaxlib==0.1.59 numpy==1.19.5 tensorflow==2.4.0

1. Create a double gaussian amp with SymPy¶

When building the model, we should be careful to pass the parameters as arguments as well, otherwise frameworks like jax can’t determine the gradient.

import math

import sympy as sp

x, A1, mu1, sigma1, A2, mu2, sigma2 = sp.symbols(
    "x, A1, mu1, sigma1, A2, mu2, sigma2"
)
gaussian1 = (
    A1
    / (sigma1 * sp.sqrt(2.0 * math.pi))
    * sp.exp(-((x - mu1) ** 2) / (2 * sigma1))
)
gaussian2 = (
    A2
    / (sigma2 * sp.sqrt(2.0 * math.pi))
    * sp.exp(-((x - mu2) ** 2) / (2 * sigma2))
)

gauss_sum = gaussian1 + gaussian2
gauss_sum

\[\displaystyle \frac{0.398942280401433 A_{1} e^{- \frac{\left(- \mu_{1} + x\right)^{2}}{2 \sigma_{1}}}}{\sigma_{1}} + \frac{0.398942280401433 A_{2} e^{- \frac{\left(- \mu_{2} + x\right)^{2}}{2 \sigma_{2}}}}{\sigma_{2}}\]

2. Convert this expression to a function using lambdify¶

TensorFlow as backend:

import inspect

tf_gauss_sum = sp.lambdify(
    (x, A1, mu1, sigma1, A2, mu2, sigma2), gauss_sum, "tensorflow"
)
print(inspect.getsource(tf_gauss_sum))

def _lambdifygenerated(x, A1, mu1, sigma1, A2, mu2, sigma2):
    return (0.398942280401433*A1*exp(-1/2*pow(-mu1 + x, 2)/sigma1)/sigma1 + 0.398942280401433*A2*exp(-1/2*pow(-mu2 + x, 2)/sigma2)/sigma2)

NumPy as backend:

numpy_gauss_sum = sp.lambdify(
    (x, A1, mu1, sigma1, A2, mu2, sigma2), gauss_sum, "numpy"
)
print(inspect.getsource(numpy_gauss_sum))

def _lambdifygenerated(x, A1, mu1, sigma1, A2, mu2, sigma2):
    return (0.398942280401433*A1*exp(-1/2*(-mu1 + x)**2/sigma1)/sigma1 + 0.398942280401433*A2*exp(-1/2*(-mu2 + x)**2/sigma2)/sigma2)

Jax as backend:

from jax import numpy as jnp
from jax import scipy as jsp

jax_gauss_sum = sp.lambdify(
    (x, A1, mu1, sigma1, A2, mu2, sigma2),
    gauss_sum,
    modules=(jnp, jsp.special),
)
print(inspect.getsource(jax_gauss_sum))

def _lambdifygenerated(x, A1, mu1, sigma1, A2, mu2, sigma2):
    return (0.398942280401433*A1*exp(-1/2*(-mu1 + x)**2/sigma1)/sigma1 + 0.398942280401433*A2*exp(-1/2*(-mu2 + x)**2/sigma2)/sigma2)

3. Natively create the respective packages¶

import math

import tensorflow as tf


def gaussian(x, A, mu, sigma):
    return (
        A
        / (sigma * tf.sqrt(tf.constant(2.0, dtype=tf.float64) * math.pi))
        * tf.exp(
            -tf.pow(-tf.constant(0.5, dtype=tf.float64) * (x - mu) / sigma, 2)
        )
    )


def native_tf_gauss_sum(x_, A1_, mu1_, sigma1_, A2_, mu2_, sigma2_):
    return gaussian(x_, A1_, mu1_, sigma1_) + gaussian(x_, A2_, mu2_, sigma2_)


# @jx.pmap
def jax_gaussian(x, A, mu, sigma):
    return (
        A
        / (sigma * jnp.sqrt(2.0 * math.pi))
        * jnp.exp(-((-0.5 * (x - mu) / sigma) ** 2))
    )


def native_jax_gauss_sum(x_, A1_, mu1_, sigma1_, A2_, mu2_, sigma2_):
    return jax_gaussian(x_, A1_, mu1_, sigma1_) + jax_gaussian(
        x_, A2_, mu2_, sigma2_
    )

4. Compare performance¶

import numpy as np

parameter_values = (1.0, 0.0, 0.1, 2.0, 2.0, 0.2)
np_x = np.random.uniform(-1, 3, 10000)
tf_x = tf.constant(np_x)


def evaluate_with_parameters(function):
    def wrapper():
        return function(np_x, *(parameter_values))

    return wrapper


def call_native_tf():
    func = native_tf_gauss_sum
    params = tuple(tf.Variable(v, dtype=tf.float64) for v in parameter_values)

    def wrapper():
        return func(tf_x, *params)

    return wrapper

import timeit

from jax.config import config

config.update("jax_enable_x64", True)

print(
    "sympy tf lambdify",
    timeit.timeit(evaluate_with_parameters(tf_gauss_sum), number=100),
)
print(
    "sympy numpy lambdify",
    timeit.timeit(evaluate_with_parameters(numpy_gauss_sum), number=100),
)
print(
    "sympy jax lambdify",
    timeit.timeit(evaluate_with_parameters(jax_gauss_sum), number=100),
)
print("native tf", timeit.timeit(call_native_tf(), number=100))

print(
    "native jax",
    timeit.timeit(evaluate_with_parameters(native_jax_gauss_sum), number=100),
)

sympy tf lambdify 0.0923923499999546
sympy numpy lambdify 0.021591910000097414
sympy jax lambdify 0.15498641800013502
native tf 0.1038757030000852
native jax 0.19923505700012356

5. Handling parameters¶

Some options:

5.1 Changing parameter values¶

Can be done in the model itself…

But how can the values be propagated to the AmplitudeModel?

Well, if an amplitude model only defines parameters with a name and the values are supplied in the function evaluation, then everything is decoupled and there are no problems.

5.2 Changing parameter names¶

Names can be changed in the sympy AmplitudeModel. Since this sympy model serves as the source of truth for the Function, all things generated from this model will reflect the name changes as well.

But going even further, since the Parameters are passed into the functions as arguments, the whole naming becomes irrelevant anyways.

tf_var_A1 = tf.Variable(1.0, dtype=tf.float64) <- does not carry a name!!

5.3 Coupling parameters¶

This means that one parameter is just assigned to another one?

result = evaluate_with_parameters(jax_gauss_sum)()
result

DeviceArray([2.8856688 , 1.81293941, 0.8762776 , ..., 1.57889689,
             3.45055759, 3.30175027], dtype=float64)

np_x

array([1.64006008, 1.43829976, 2.77864469, ..., 1.39104059, 2.24092382,
       1.72490419])

import matplotlib.pyplot as plt

plt.hist(np_x, bins=100, weights=result);

parameter_values = (1.0, 0.0, 0.1, 2.0, 2.0, 0.1)
result = evaluate_with_parameters(jax_gauss_sum)()
plt.hist(np_x, bins=100, weights=result);

6. Exchange a gaussian with some other function¶

This should be easy if you know the exact expression that you want to replace:

from sympy.abc import C, a, b, x

expr = sp.sin(a * x) + sp.cos(b * x)
expr

\[\displaystyle \sin{\left(a x \right)} + \cos{\left(b x \right)}\]

expr.subs(sp.sin(a * x), C)

\[\displaystyle C + \cos{\left(b x \right)}\]

7. Matrix operations?¶

from sympy.physics.quantum.dagger import Dagger

spin_density = sp.MatrixSymbol("rho", 3, 3)
amplitudes = sp.Matrix([[1 + sp.I], [2 + sp.I], [3 + sp.I]])

dummy_intensity = sp.re(
    Dagger(amplitudes) * spin_density * amplitudes,
    evaluate=False
    # evaluate=False is important otherwise it generates some function that cant
    # be lambdified anymore
)
dummy_intensity

\[\begin{split}\displaystyle \operatorname{re}{\left(\left[\begin{matrix}1 - i & 2 - i & 3 - i\end{matrix}\right] \rho \left[\begin{matrix}1 + i\\2 + i\\3 + i\end{matrix}\right]\right)}\end{split}\]

tf_intensity = sp.lambdify(
    spin_density,
    dummy_intensity,
    modules=(tf,),
)
print(inspect.getsource(tf_intensity))

def _lambdifygenerated(rho):
    return (real(matmul(matmul(constant([[1 - 1j, 2 - 1j, 3 - 1j]]), rho), constant([[1 + 1j], [2 + 1j], [3 + 1j]]))))

real0 = tf.constant(0, dtype=tf.float64)
real1 = tf.constant(1, dtype=tf.float64)
intensity_result = tf_intensity(
    np.array(
        [
            [
                tf.complex(real1, real0),
                tf.complex(real0, real0),
                -tf.complex(real0, real1),
            ],
            [
                tf.complex(real0, real0),
                tf.complex(real1, real0),
                tf.complex(real0, real0),
            ],
            [
                tf.complex(real0, real1),
                tf.complex(real0, real0),
                tf.complex(real1, real0),
            ],
        ]
    ),
)
intensity_result

<tf.Tensor: shape=(1, 1), dtype=float64, numpy=array([[13.]])>

Python operator library¶

See Python’s built-in operator library

What we have now¶

Build test model

import qrules as q

result = q.generate_transitions(
    initial_state=[("J/psi(1S)", [-1, 1])],
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=["N(1440)"],
    allowed_interaction_types="strong",
)
model = q.generate_amplitudes(result)
for particle in result.get_intermediate_particles():
    model.dynamics.set_breit_wigner(particle.name)
q.io.write(model, "recipe.yml")

Visualize the decay:

import graphviz

graphs = result.collapse_graphs()
dot = q.io.asdot(graphs)
graphviz.Source(dot)

model.parameters

FitParameters([
    FitParameter(name='Magnitude_J/psi(1S)_to_N(1440)+_0.5+p~_-0.5;N(1440)+_to_eta_0+p_0.5;', value=1.0, fix=False),
    FitParameter(name='Magnitude_J/psi(1S)_to_N(1440)+_0.5+p~_0.5;N(1440)+_to_eta_0+p_0.5;', value=1.0, fix=False),
    FitParameter(name='Magnitude_J/psi(1S)_to_N(1440)~-_0.5+p_-0.5;N(1440)~-_to_eta_0+p~_0.5;', value=1.0, fix=False),
    FitParameter(name='Magnitude_J/psi(1S)_to_N(1440)~-_0.5+p_0.5;N(1440)~-_to_eta_0+p~_0.5;', value=1.0, fix=False),
    FitParameter(name='MesonRadius_J/psi(1S)', value=1.0, fix=True),
    FitParameter(name='MesonRadius_N(1440)+', value=1.0, fix=True),
    FitParameter(name='MesonRadius_N(1440)~-', value=1.0, fix=True),
    FitParameter(name='Phase_J/psi(1S)_to_N(1440)+_0.5+p~_-0.5;N(1440)+_to_eta_0+p_0.5;', value=0.0, fix=False),
    FitParameter(name='Phase_J/psi(1S)_to_N(1440)+_0.5+p~_0.5;N(1440)+_to_eta_0+p_0.5;', value=0.0, fix=False),
    FitParameter(name='Phase_J/psi(1S)_to_N(1440)~-_0.5+p_-0.5;N(1440)~-_to_eta_0+p~_0.5;', value=0.0, fix=False),
    FitParameter(name='Phase_J/psi(1S)_to_N(1440)~-_0.5+p_0.5;N(1440)~-_to_eta_0+p~_0.5;', value=0.0, fix=False),
    FitParameter(name='Position_N(1440)+', value=1.44, fix=False),
    FitParameter(name='Position_N(1440)~-', value=1.44, fix=False),
    FitParameter(name='Width_N(1440)+', value=0.35, fix=False),
    FitParameter(name='Width_N(1440)~-', value=0.35, fix=False),
])

Implementation with operators¶

See this answer on Stack Overflow:

import operator

MAKE_BINARY = lambda opfn: lambda self, other: BinaryOp(  # noqa: E731
    self, asMagicNumber(other), opfn
)
MAKE_RBINARY = lambda opfn: lambda self, other: BinaryOp(  # noqa: E731
    asMagicNumber(other), self, opfn
)


class MagicNumber(object):
    __add__ = MAKE_BINARY(operator.add)
    __sub__ = MAKE_BINARY(operator.sub)
    __mul__ = MAKE_BINARY(operator.mul)
    __radd__ = MAKE_RBINARY(operator.add)
    __rsub__ = MAKE_RBINARY(operator.sub)
    __rmul__ = MAKE_RBINARY(operator.mul)
    # __div__  = MAKE_BINARY(operator.div)
    # __rdiv__ = MAKE_RBINARY(operator.div)
    __truediv__ = MAKE_BINARY(operator.truediv)
    __rtruediv__ = MAKE_RBINARY(operator.truediv)
    __floordiv__ = MAKE_BINARY(operator.floordiv)
    __rfloordiv__ = MAKE_RBINARY(operator.floordiv)

    def __neg__(self, other):
        return UnaryOp(self, lambda x: -x)

    @property
    def value(self):
        return self.eval()


class Constant(MagicNumber):
    def __init__(self, value):
        self.value_ = value

    def eval(self):
        return self.value_


class Parameter(Constant):
    def __init__(self):
        super(Parameter, self).__init__(0.0)

    def setValue(self, v):
        self.value_ = v

    value = property(fset=setValue, fget=lambda self: self.value_)


class BinaryOp(MagicNumber):
    def __init__(self, op1, op2, operation):
        self.op1 = op1
        self.op2 = op2
        self.opn = operation

    def eval(self):
        return self.opn(self.op1.eval(), self.op2.eval())


class UnaryOp(MagicNumber):
    def __init__(self, op1, operation):
        self.op1 = op1
        self.operation = operation

    def eval(self):
        return self.opn(self.op1.eval())


asMagicNumber = (
    lambda x: x if isinstance(x, MagicNumber) else Constant(x)  # noqa: E731
)
asMagicNumber(2).eval()

Other ideas¶

Option 1: parameter container

Remove name from the FitParameter class and give the FitParameters collection class the responsibility to keep track of ‘names’ of the FitParameters as keys in a dict. In the AmplitudeModel, locations where a FitParameter should be inserted are indicated by an immutable (!) str that should exist as a key in the FitParameters.

Such a setup best reflects the structure of the AmplitudeModel that we have now (best illustrated by expected_recipe, note in particular YAML anchors like &par1/*par1). It also allows one to couple FitParameters. See following snippet:

import attr


# the new FitParameter class would have this structure
@attr.s
class Parameter:
    value: float = attr.ib()
    fix: bool = attr.ib(default=False)


# the new FitParameters collection would have such a structure
mapping = {
    "par1": Parameter(1.0),
    "par2": Parameter(2.0, fix=False),
}


# intensity nodes and dynamics classes contain immutable strings
class Dynamics:
    pass


@attr.s
class CustomDynamics(Dynamics):
    par: str = attr.ib(on_setattr=attr.setters.frozen, kw_only=True)


dyn1 = CustomDynamics(par="par1")
dyn2 = CustomDynamics(par="par2")

# Parameters would be coupled like this
mapping["par1"] = mapping["par2"]
assert mapping["par2"] is mapping["par1"]
assert mapping["par1"] == {
    "par1": Parameter(1.0),
    "par2": Parameter(1.0),
}

Option 2: read-only parameter manager

Remove the FitParameters collection class altogether and use something like immutable InitialParameter instances in the dynamics and intensity section of the AmplitudeModel. The AmplitudeModel then starts to serve as a read-only’ template. A fitter package like tensorwaves can then loop over the AmplitudeModel structure to extract the InitialParameter instances and convert them to something like an FitParameter.

Here’s a rough sketch with tensorwaves in mind.

from typing import Dict, Generator, List

import attr

from expertsystem.amplitude.model import (
    AmplitudeModel,
    Dynamics,
    Node,
    ParticleDynamics,
)
from qrules.particle import Particle


@attr.s
class InitialParameter:
    name: str = attr.ib()
    value: float = attr.ib()
    # fix: bool = attr.ib(default=False)


@attr.s
class FitParameter:
    name: str = attr.ib(on_setattr=attr.setters.frozen)
    value: float = attr.ib()
    fix: bool = attr.ib(default=False)


class FitParameterManager:
    """Manages all fit parameters of the model"""

    def __init__(self, model: AmplitudeModel) -> None:
        self.__model: AmplitudeModel
        self.__parameter_couplings: Dict[str, str]

    @property
    def parameters(self) -> List[FitParameter]:
        initial_parameters = list(__yield_parameter(self.__model))
        self.__apply_couplings()
        return self.__convert(initial_parameters)

    def couple_parameters(self, parameter1: str, parameter2: str) -> None:
        pass

    def __convert(self, params: List[InitialParameter]) -> List[FitParameter]:
        pass


@attr.s
class CustomDynamics(Dynamics):
    parameter: InitialParameter = attr.ib(kw_only=True)

    @staticmethod
    def from_particle(particle: Particle):
        pass


def __yield_parameter(
    instance: object,
) -> Generator[InitialParameter, None, None]:
    if isinstance(instance, InitialParameter):
        yield instance
    elif isinstance(instance, (ParticleDynamics, Node)):
        for item in instance.values():
            yield from __yield_parameter(item)
    elif isinstance(instance, (list,)):
        for item in instance:
            yield from __yield_parameter(item)
    elif attr.has(instance.__class__):
        for field in attr.fields(instance.__class__):
            field_value = getattr(instance, field.name)
            yield from __yield_parameter(field_value)


# usage in tensorwaves
amp_model = AmplitudeModel()
kinematics: HelicityKinematics = ...
builder = IntensityBuilder(kinematics)

intensity = builder.create(amp_model)  # this would call amp_model.parameters
parameters: Dict[str, float] = intensity.parameters
# PROBLEM?: fix status is lost at this point

data_sample = generate_data(...)
dataset = kinematics.convert(data_sample)

parameters["Width_f(0)(980)"] = 0.2  # name is immutable at this point

# name of a parameter can be changed in the AmplitudeModel though
# and then call builder again
intensity(dataset, parameters)

Evaluation¶

Pros and Cons¶

Customized Python classes (current state)¶

Positive
- “Faster” implementation / prototyping possible compared to python operators
- No additional dependencies
Negative
- Not open-closed to new models
- Conversion to various back-ends not DRY
- Function replacement or extension feature becomes very difficult to handle.
- Model is not complete, since no complete mathematical description is used. For example Breit-Wigner functions are referred to directly and their implementations is not defined in the amplitude model.

SymPy ¶

Positive
- Easy to render amplitude model as LaTeX
- Model description is complete! Absolutely all information about the model is included. (reproducibility)
- Follows open-closed principle. New models and formalism can be added without any changes to other interfacing components (here: tensorwaves)
- Use lambdify to convert the expression to any back-end
- Use Expr.subs (substitute) to couple parameters or replace components of the model, for instance to set custom dynamics
Negative
- lambdify becomes a core dependency while its behavior cannot be modified, but is defined by sympy.
- Need to keep track of components in the expression tree with symbol mappings

Python’s operator library ¶

Positive
- More control over different components of in the expression tree
- More control over convert functionality to functions
- No additional dependencies
Negative
- Essentially re-inventing SymPy

Decision outcome¶

Use SymPy. Initially, we leave the existing amplitude builders (modules helicity_decay and canonical_decay) alongside a SymPy implementation, so that it’s possible to compare the results. Once it turns out the this set-up results in the same results and a comparable performance, we replace the old amplitude builders with the new SymPy implementation.

[ADR-002] Inserting dynamics¶

Status: proposed
Deciders: @redeboer @spflueger

Context and problem statement¶

Physics models usually include assumptions that simplify the structure of the model. For example, splitting a model into a product of independent parts, in which every part contains a certain responsibility. In case of partial wave amplitude models, we can make a separation into a spin part and a dynamical part. While the spin part controls the probability w.r.t angular kinematic variables, the dynamics controls the probability on variable like the invariant mass of states.

Generally, a dynamics part is simply a function, which is defined in complex space, and consists of:

a mathematical expression (sympy.Expr)
a set of parameters in that expression that can be tweaked (optimized)
a set of (kinematic) variables to which the expression applies

Technical story¶

#124: see Issues with existing set-up.
#440: no way to supply custom dynamics. Or at least, tensorwaves does not understand those custom dynamics.
ADR-001: parameters and variables are to be expressed as sympy.Symbols.
#454: dynamics are specified as a mapping of sympy.Function to a sympy.Expr, but now there is no way to supply those sympy.Exprs with expected sympy.Symbols (parameters and variables).

Issues with existing set-up¶

There is no clear way to apply dynamics functions to a specific decaying particle, that is, to a specific edge of the StateTransitionGraphs (STG). Currently, we just work with a mapping of Particles to some dynamics expression, but this is not feasible when there there are identical particles on the edges.
The set of variables to which a dynamics expression applies, is determined by the position within the STG that it is applied to. For instance, a relativistic Breit-Wigner that is applied to the resonance in some 1-to-3 isobar decay (described by an STG with final state edges 2, 3, 4 and intermediate edge 1) would work on the invariant mass of edge 3 and 4 (mSq_3+4).
Just like variables, parameters need to be identifiable from their position within the STG (take a relativistic Breit-Wigner with form factors, which would require break-up momentum as a parameter), but also require some suggested starting value (e.g. expected pole position). These starting values are usually taken from the edge and node properties within the STG.

Decision drivers¶

The following points are nice to have or can influence the decision but are not essential and can be part of the users responsibility.

The parameters that a dynamics functions requires, are registered automatically and linked together.
Kinematic variables used in dynamics functions are also linked appropriately.
It is easy to define custom dynamics (no boilerplate code).

Solution requirements¶

It is easy to apply dynamics to specific components of the STGs. Note: it’s important that the dynamics can be applied to resonances of some selected graphs and not generally all graphs in which the resonance appears.
Where possible, suggested (initial) parameter values are provided as well.
It is possible to use and inspect the dynamics expression itself independently from the expertsystem.
Follow open-closed principle. Probably the most important decision driver. The solution should be flexible enough to handle any possible scenario, without having to change the interface defined in requirement 1!

Considered solutions¶

Group 1: expression builder¶

To satisfy requirement 1, we propose the following syntax:

# model: ModelInfo
# graph: StateTransitionGraph
model.set_dynamics(graph, edge_id=1, expression_builder)

Another style would be to have ModelInfo contain a reference to the list of StateTransitionGraphs. The user then needs some other way to express which edges to apply the dynamics function to:

model.set_dynamics(
    filter=lambda p: p.name.startswith("f(0)"),
    edge_id=1,
    expression_builder,
)

Here, expression_builder is some function or method that can create a dynamics expression. It can also be a class that contains both the implementation of the expression and a static method to build itself from a StateTransitionGraph.

The dynamics expression needs to be formulated in such a way that it satisfies the rest of the requirements. The following options illustrate three different ways of formulating a dynamics expression, each taking a relativistic Breit-Wigner and a relativistic Breit-Wigner with form factor as example.

Using composition¶

from typing import Tuple

import attr
import sympy as sp
from helpers import (
    StateTransitionGraph,
    blatt_weisskopf,
    determine_attached_final_state,
    two_body_momentum_squared,
)

try:
    from typing import Protocol
except ImportError:
    from typing_extensions import Protocol  # type: ignore

A frozen DynamicExpression class keeps track of variables, parameters, and the dynamics expression in which they should appear:

@attr.s(frozen=True)
class DynamicsExpression:
    variables: Tuple[sp.Symbol, ...] = attr.ib()
    parameters: Tuple[sp.Symbol, ...] = attr.ib()
    expression: sp.Expr = attr.ib()

    def substitute(self) -> sp.Expr:
        return self.expression(*self.variables, *self.parameters)

The expression attribute can be formulated as a simple Python function that takes sympy.Symbols as arguments and returns a sympy.Expr:

def relativistic_breit_wigner(
    mass: sp.Symbol, mass0: sp.Symbol, gamma0: sp.Symbol
) -> sp.Expr:
    return gamma0 * mass0 / (mass0 ** 2 - mass ** 2 - gamma0 * mass0 * sp.I)

def relativistic_breit_wigner_with_form_factor(
    mass: sp.Symbol,
    mass0: sp.Symbol,
    gamma0: sp.Symbol,
    m_a: sp.Symbol,
    m_b: sp.Symbol,
    angular_momentum: sp.Symbol,
    meson_radius: sp.Symbol,
) -> sp.Expr:
    q_squared = two_body_momentum_squared(mass, m_a, m_b)
    q0_squared = two_body_momentum_squared(mass0, m_a, m_b)
    ff2 = blatt_weisskopf(q_squared, meson_radius, angular_momentum)
    ff02 = blatt_weisskopf(q0_squared, meson_radius, angular_momentum)
    width = gamma0 * (mass0 / mass) * (ff2 / ff02)
    width = width * sp.sqrt(q_squared / q0_squared)
    return (
        relativistic_breit_wigner(mass, mass0, width)
        * mass0
        * gamma0
        * sp.sqrt(ff2)
    )

The DynamicsExpression container class enables us to provide the expression with correctly named Symbols for the decay that is being described. Here, we use some naming scheme for an \(f_0(980)\) decaying to final state edges 3 and 4 (say \(\pi^0\pi^0\)):

bw_decay_f0 = DynamicsExpression(
    variables=sp.symbols("m_3+4", seq=True),
    parameters=sp.symbols(R"m_f(0)(980) \Gamma_f(0)(980)"),
    expression=relativistic_breit_wigner,
)
bw_decay_f0.substitute()

\[\displaystyle \frac{\Gamma_{f(0)(980)} m_{f(0)(980)}}{- i \Gamma_{f(0)(980)} m_{f(0)(980)} - m_{3+4}^{2} + m_{f(0)(980)}^{2}}\]

For each dynamics expression, we have to provide a ‘builder’ function that can create a DynamicsExpression for a specific edge within the StateTransitionGraph:

def relativistic_breit_wigner_from_graph(
    graph: StateTransitionGraph, edge_id: int
) -> DynamicsExpression:
    edge_ids = determine_attached_final_state(graph, edge_id)
    final_state_ids = map(str, edge_ids)
    mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
    particle, _ = graph.get_edge_props(edge_id)
    mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
    gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
    return DynamicsExpression(
        variables=(mass),
        parameters=(mass0, gamma0),
        expression=relativistic_breit_wigner(mass, mass0, gamma0),
    )

def relativistic_breit_wigner_with_form_factor_from_graph(
    graph: StateTransitionGraph, edge_id: int
) -> DynamicsExpression:
    edge_ids = determine_attached_final_state(graph, edge_id)
    final_state_ids = map(str, edge_ids)
    mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
    particle, _ = graph.get_edge_props(edge_id)
    mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
    gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
    m_a = sp.Symbol(f"m_{edge_ids[0]}")
    m_b = sp.Symbol(f"m_{edge_ids[1]}")
    angular_momentum = particle.spin  # helicity formalism only!
    meson_radius = sp.Symbol(fR"R_{{{particle.latex}}}")
    return DynamicsExpression(
        variables=(mass),
        parameters=(mass0, gamma0, m_a, m_b, angular_momentum, meson_radius),
        expression=relativistic_breit_wigner_with_form_factor(
            mass, mass0, gamma0, m_a, m_b, angular_momentum, meson_radius
        ),
    )

The fact that DynamicsExpression.expression is just a Python function, allows one to inspect the dynamics formulation of these functions independently, purely in terms of SymPy:

m, m0, w0 = sp.symbols(R"m m_0 \Gamma")
evaluated_bw = relativistic_breit_wigner(m, 1.0, 0.3)
sp.plot(sp.Abs(evaluated_bw), (m, 0, 2), axis_center=(0, 0), ylim=(0, 1))
sp.plot(sp.arg(evaluated_bw), (m, 0, 2), axis_center=(0, 0), ylim=(0, sp.pi))
relativistic_breit_wigner(m, m0, w0)

\[\displaystyle \frac{\Gamma m_{0}}{- i \Gamma m_{0} - m^{2} + m_{0}^{2}}\]

This closes the gap between the code and the theory that is being implemented.

Alternative signature¶

An alternative way to specify the expression is:

def expression(
    variables: Tuple[sp.Symbol, ...], parameters: Tuple[sp.Symbol, ...]
) -> sp.Expr:
    pass

Here, one would however need to unpack the variables and parameters. The advantage is that the signature becomes more general.

Type checking¶

There is no way to enforce the appropriate signature on the builder function, other than following a Protocol:

class DynamicsBuilder(Protocol):
    def __call__(
        self, graph: StateTransitionGraph, edge_id: int
    ) -> DynamicsExpression:
        ...

This DynamicsBuilder protocol would be used in the syntax proposed at Considered solutions.

It carries another subtle problem, though: a Protocol is only used in static type checking, while potential problems with the implementation of the builder and its respective expression only arrise at runtime.

Subclassing sympy.Function¶

from abc import abstractmethod

import sympy as sp
from helpers import (
    StateTransitionGraph,
    blatt_weisskopf,
    determine_attached_final_state,
    two_body_momentum_squared,
)

One way to address the Cons of Using composition, is to sub-class Function. The expression is implemented by overwriting the inherited eval() method and the builder is provided through the class through an additional from_graph class method. The interface would look like this:

class DynamicsFunction(sp.Function):
    @classmethod
    @abstractmethod
    def eval(cls, *args: sp.Symbol) -> sp.Expr:
        """Implementation of the dynamics function."""

    @classmethod
    @abstractmethod
    def from_graph(cls, graph: StateTransitionGraph, edge_id: int) -> sp.Basic:
        pass

As can be seen from the implementation of a relativistic Breit-Wigner, the implementation of the expression is nicely kept together with the implementation of the expression:

class RelativisticBreitWigner(DynamicsFunction):
    @classmethod
    def eval(cls, *args: sp.Symbol) -> sp.Expr:
        mass = args[0]
        mass0 = args[1]
        gamma0 = args[2]
        return (
            gamma0 * mass0 / (mass0 ** 2 - mass ** 2 - gamma0 * mass0 * sp.I)
        )

    @classmethod
    def from_graph(
        cls, graph: StateTransitionGraph, edge_id: int
    ) -> "RelativisticBreitWigner":
        edge_ids = determine_attached_final_state(graph, edge_id)
        final_state_ids = map(str, edge_ids)
        mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
        particle, _ = graph.get_edge_props(edge_id)
        mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
        gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
        return cls(mass, mass0, gamma0)

It becomes a bit less clear when using a form factor, but the DynamicsFunction base class enforces a correct interfaces:

class RelativisticBreitWignerWithFF(DynamicsFunction):
    @classmethod
    def eval(cls, *args: sp.Symbol) -> sp.Expr:
        # Arguments
        mass = args[0]
        mass0 = args[1]
        gamma0 = args[2]
        m_a = args[3]
        m_b = args[4]
        angular_momentum = args[5]
        meson_radius = args[6]
        # Computed variables
        q_squared = two_body_momentum_squared(mass, m_a, m_b)
        q0_squared = two_body_momentum_squared(mass0, m_a, m_b)
        ff2 = blatt_weisskopf(q_squared, meson_radius, angular_momentum)
        ff02 = blatt_weisskopf(q0_squared, meson_radius, angular_momentum)
        width = gamma0 * (mass0 / mass) * (ff2 / ff02)
        width = width * sp.sqrt(q_squared / q0_squared)
        # Expression
        return (
            RelativisticBreitWigner(mass, mass0, width)
            * mass0
            * gamma0
            * sp.sqrt(ff2)
        )

    @classmethod
    def from_graph(
        cls, graph: StateTransitionGraph, edge_id: int
    ) -> "RelativisticBreitWignerWithFF":
        edge_ids = determine_attached_final_state(graph, edge_id)
        final_state_ids = map(str, edge_ids)
        mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
        particle, _ = graph.get_edge_props(edge_id)
        mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
        gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
        m_a = sp.Symbol(f"m_{edge_ids[0]}")
        m_b = sp.Symbol(f"m_{edge_ids[1]}")
        angular_momentum = particle.spin  # helicity formalism only!
        meson_radius = sp.Symbol(fR"R_{{{particle.latex}}}")
        return cls(
            mass, mass0, gamma0, m_a, m_b, angular_momentum, meson_radius
        )

The expression_builder used in the syntax proposed at Considered solutions, would now just be a class that is derived of DynamicsFunction.

The sympy.Function class provides mixin methods, so that the derived class behaves as a sympy expression. So the expression can be inspected with the usual sympy tools (compare the Pros of Using composition):

m, m0, w0 = sp.symbols(R"m m_0 \Gamma")
evaluated_bw = RelativisticBreitWigner(m, 1.0, 0.3)
sp.plot(sp.Abs(evaluated_bw), (m, 0, 2), axis_center=(0, 0), ylim=(0, 1))
sp.plot(sp.arg(evaluated_bw), (m, 0, 2), axis_center=(0, 0), ylim=(0, sp.pi))
RelativisticBreitWigner(m, m0, w0)

\[\displaystyle \frac{\Gamma m_{0}}{- i \Gamma m_{0} - m^{2} + m_{0}^{2}}\]

Subclassing sympy.Expr¶

from abc import abstractmethod
from typing import Any, Callable, Set, Type

import sympy as sp
from helpers import (
    StateTransitionGraph,
    blatt_weisskopf,
    determine_attached_final_state,
    two_body_momentum_squared,
)

The major disadvantage of Subclassing sympy.Function, is that there is no way to identify which Symbols are variables and which are parameters. This can be solved by sub-classing from sympy.core.expr.Expr.

An example of a class that does this is WignerD. There, the implementation of the dynamics expression can be evaluated through a doit() call. This method can call anything, but sympy seems to follow the convention that it returns an ‘evaluated’ version of the class itself, where ‘evaluated’ means that any randomly named method of the class has been called on the *args that are implemented through the __new__ method (the examples below make this clearer).

For our purposes, the follow DynamicsExpr base class illustrates the interface that we expect. Here, evaluate is where expression is implemented and (just as in Subclassing sympy.Function) from_graph is the builder method.

class DynamicsExpr(sp.Expr):
    @classmethod
    @abstractmethod
    def __new__(cls, *args: sp.Symbol, **hints: Any) -> sp.Expr:
        pass

    @abstractmethod
    def doit(self, **hints: Any) -> sp.Expr:
        pass

    @abstractmethod
    def evaluate(self) -> sp.Expr:
        pass

    @classmethod
    @abstractmethod
    def from_graph(cls, graph: StateTransitionGraph, edge_id: int) -> sp.Basic:
        pass

The __new__ and doit methods split the construction from the evaluation of the expression. This allows one to distinguish variables and parameters and present them as properties:

class RelativisticBreitWigner(DynamicsExpr):
    def __new__(cls, *args: sp.Symbol, **hints: Any) -> sp.Expr:
        if len(args) != 3:
            raise ValueError(f"3 parameters expected, got {len(args)}")
        args = sp.sympify(args)
        evaluate = hints.get("evaluate", False)
        if evaluate:
            return sp.Expr.__new__(cls, *args).evaluate()  # type: ignore  # pylint: disable=no-member
        return sp.Expr.__new__(cls, *args)

    @property
    def mass(self) -> sp.Symbol:
        return self.args[0]

    @property
    def mass0(self) -> sp.Symbol:
        return self.args[1]

    @property
    def gamma0(self) -> sp.Symbol:
        return self.args[2]

    @property
    def variables(self) -> Set[sp.Symbol]:
        return {self.mass}

    @property
    def parameters(self) -> Set[sp.Symbol]:
        return {self.mass0, self.gamma0}

    def doit(self, **hints: Any) -> sp.Expr:
        return RelativisticBreitWigner(*self.args, **hints, evaluate=True)

    def evaluate(self) -> sp.Expr:
        return (
            self.gamma0
            * self.mass0
            / (
                self.mass0 ** 2
                - self.mass ** 2
                - self.gamma0 * self.mass0 * sp.I
            )
        )

    @classmethod
    def from_graph(
        cls, graph: StateTransitionGraph, edge_id: int
    ) -> "RelativisticBreitWigner":
        edge_ids = determine_attached_final_state(graph, edge_id)
        final_state_ids = map(str, edge_ids)
        mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
        particle, _ = graph.get_edge_props(edge_id)
        mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
        gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
        return cls(mass, mass0, gamma0)

class RelativisticBreitWignerWithFF(DynamicsExpr):
    def __new__(cls, *args: sp.Symbol, **hints: Any) -> sp.Expr:
        if len(args) != 7:
            raise ValueError(f"7 parameters expected, got {len(args)}")
        args = sp.sympify(args)
        evaluate = hints.get("evaluate", False)
        if evaluate:
            return sp.Expr.__new__(cls, *args).evaluate()  # type: ignore  # pylint: disable=no-member
        return sp.Expr.__new__(cls, *args)

    def doit(self, **hints: Any) -> sp.Expr:
        return RelativisticBreitWignerWithFF(
            *self.args, **hints, evaluate=True
        )

    @property
    def mass(self) -> sp.Symbol:
        return self.args[0]

    @property
    def mass0(self) -> sp.Symbol:
        return self.args[1]

    @property
    def gamma0(self) -> sp.Symbol:
        return self.args[2]

    @property
    def m_a(self) -> sp.Symbol:
        return self.args[3]

    @property
    def m_b(self) -> sp.Symbol:
        return self.args[4]

    @property
    def angular_momentum(self) -> sp.Symbol:
        return self.args[5]

    @property
    def meson_radius(self) -> sp.Symbol:
        return self.args[6]

    def evaluate(self) -> sp.Expr:
        # Computed variables
        q_squared = two_body_momentum_squared(self.mass, self.m_a, self.m_b)
        q0_squared = two_body_momentum_squared(self.mass0, self.m_a, self.m_b)
        ff2 = blatt_weisskopf(
            q_squared, self.meson_radius, self.angular_momentum
        )
        ff02 = blatt_weisskopf(
            q0_squared, self.meson_radius, self.angular_momentum
        )
        width = self.gamma0 * (self.mass0 / self.mass) * (ff2 / ff02)
        width = width * sp.sqrt(q_squared / q0_squared)
        # Expression
        return (
            RelativisticBreitWigner(self.mass, self.mass0, width)
            * self.mass0
            * self.gamma0
            * sp.sqrt(ff2)
        )

    @classmethod
    def from_graph(
        cls, graph: StateTransitionGraph, edge_id: int
    ) -> "RelativisticBreitWignerWithFF":
        edge_ids = determine_attached_final_state(graph, edge_id)
        final_state_ids = map(str, edge_ids)
        mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
        particle, _ = graph.get_edge_props(edge_id)
        mass0 = sp.Symbol(f"m_{{{particle.latex}}}")
        gamma0 = sp.Symbol(fR"\Gamma_{{{particle.latex}}}")
        m_a = sp.Symbol(f"m_{edge_ids[0]}")
        m_b = sp.Symbol(f"m_{edge_ids[1]}")
        angular_momentum = particle.spin  # helicity formalism only!
        meson_radius = sp.Symbol(fR"R_{{{particle.latex}}}")
        return cls(
            mass, mass0, gamma0, m_a, m_b, angular_momentum, meson_radius
        )

The following illustrates the difference with Subclassing sympy.Function. First, notice that a class derived from DynamicsExpr is still identifiable upon construction:

m, m0, w0 = sp.symbols(R"m m_0 \Gamma")
rel_bw = RelativisticBreitWigner(m, m0, w0)
rel_bw

\[\displaystyle RelativisticBreitWigner\left(m, m_{0}, \Gamma\right)\]

The way in which this expression is rendered in a Jupyter notebook can be changed by overwriting the _pretty and/or _latex methods.

Only once doit() is called, is the DynamicsExpr converted into a mathematical expression:

evaluated_bw = rel_bw.doit()
evaluated_bw

\[\displaystyle \frac{\Gamma m_{0}}{- i \Gamma m_{0} - m^{2} + m_{0}^{2}}\]

sp.plot(
    sp.Abs(evaluated_bw.subs({m0: 1, w0: 0.2})),
    (m, 0, 2),
    axis_center=(0, 0),
    ylim=(0, 1),
);

Decorator¶

There are a lot of implicit conventions that need to be followed to provide a correct implementation of a DynamicsExpr. Some of this may be mitigated by proving some class decorator that can easily construct the __new__() and doit() methods for you.

def dynamics_expression(n_args: int) -> Callable[[type], Type[DynamicsExpr]]:
    def decorator(decorated_class: type) -> Type[DynamicsExpr]:
        def __new__(cls, *args: sp.Symbol, **hints: Any) -> sp.Expr:
            if len(args) != n_args:
                raise ValueError(
                    f"{n_args} parameters expected, got {len(args)}"
                )
            args = sp.sympify(args)
            evaluate = hints.get("evaluate", False)
            if evaluate:
                return sp.Expr.__new__(cls, *args).evaluate()
            return sp.Expr.__new__(cls, *args)

        def doit(self, **hints: Any) -> sp.Expr:
            return decorated_class(*self.args, **hints, evaluate=True)

        decorated_class.__new__ = __new__
        decorated_class.doit = doit
        return decorated_class

    return decorator

This saves some lines of code:

@dynamics_expression(n_args=3)
class Gauss(DynamicsExpr):
    @property
    def mass(self) -> sp.Symbol:
        return self.args[0]

    @property
    def mu(self) -> sp.Symbol:
        return self.args[1]

    @property
    def sigma(self) -> sp.Symbol:
        return self.args[2]

    @property
    def variables(self) -> Set[sp.Symbol]:
        return {self.mass}

    @property
    def parameters(self) -> Set[sp.Symbol]:
        return {self.mu, self.sigma}

    def evaluate(self) -> sp.Expr:
        return sp.exp(-((self.mass - self.mu) ** 2) / self.sigma ** 2)

    @classmethod
    def from_graph(
        cls, graph: StateTransitionGraph, edge_id: int
    ) -> "RelativisticBreitWigner":
        edge_ids = determine_attached_final_state(graph, edge_id)
        final_state_ids = map(str, edge_ids)
        mass = sp.Symbol(f"m_{{{'+'.join(final_state_ids)}}}")
        particle, _ = graph.get_edge_props(edge_id)
        mass0 = sp.Symbol(fR"\mu_{{{particle.latex}}}")
        gamma0 = sp.Symbol(fR"\sigma_{{{particle.latex}}}")
        return cls(mass, mass0, gamma0)

x, mu, sigma = sp.symbols(R"x \mu \sigma")
Gauss(x, mu, w0)

\[\displaystyle Gauss\left(x, \mu, \Gamma\right)\]

Gauss(x, mu, sigma).doit()

\[\displaystyle e^{- \frac{\left(- \mu + x\right)^{2}}{\sigma^{2}}}\]

Issue with lambdify¶

It’s not possible to plot a DynamicsExpr directly as long as no lambdify hook has been provided: doit() has to be executed first.

sp.plot(sp.Abs(rel_bw.subs({m0: 1, w0: 0.2})).doit(), (m, 0, 2));

Group 2: expression-only¶

A second branch of solutions would propose the following interface:

# model: ModelInfo
# graph: StateTransitionGraph
model.set_dynamics(graph, edge_id=1, expression)

The key difference is the usage of general sympy expression sympy.Expr as an argument instead of constructing this through some builder object.

Solution evaluation¶

1: Expression builder¶

All of the solutions have the drawback arising from the choice of interface using a expression_builder. This enforces the logic of correctly coupling variables and parameters into these builders. This is extremely hard to get right, since the code has to be able to handle arbitrarily complex models. And always knowing what the user would like to do is more or less impossible. Therefore it is much better to use a already built expression that is assumed to be correctly built (see solution group 2).

All of the solutions in this group also have the following additional drawbacks. These are however more related to the correct building of the dynamics expression:

There is an implicit assumption on the signature of the expression: the first arguments are assumed to be the (kinematic) variables and the remaining arguments are parameters. In addition, the arguments cannot be keywords, but have to be positional.
The number of variables and parameters is only verified at runtime (no static typing, other than a check that each of the elements is sympy.Symbol).

Composition is the cleanest design, but is less in tune with the design of sympy. Subclassing sympy.Function and Subclassing sympy.Expr follow sympy implementations, but result in an obscure inheritance hierarchy with implicit conventions. This can result in some nasty bugs, for instance if one were to __call__ method in either the sympy.Function or sympy.Expr implementation.

Pros and Cons that are specific to each of the implementations are listed below.

Using composition ¶

Positive
- Implementation of the expression is transparent
Negative
- Alternative signature.
- The only way to see that relativistic_breit_wigner_from_graph is the builder for relativistic_breit_wigner, is from its name. This makes it implementing custom dynamics inconvenient and error-prone.
- Signature of the builder can only be checked with a Protocol, see Type checking.

Subclassing sympy.Function ¶

Positive
- DynamicsFunction behaves as a Function
- Implementation of the builder is kept together with the implementation of the expression.
Negative
- It’s not possible to identify variables and parameters

Subclassing sympy.Expr ¶

Positive
- When recursing through the amplitude model, it is still possible to identify instances of DynamicsExpr (before doit() has been called).
- Additional properties and methods can be added and carried around by the class.
Negative
- Boilerplate code required when implementing custom dynamics
- Issue with lambdify

2: Expression-only¶

Positive: This choice of interface follows the principle of SOLID more than solution group 1. By handing a complete expression of the dynamics to the setter, its sole responsibility is to insert this expression at the correct place in the full model expression.

Negative: There are no direct negative aspects to this solution, as it just splits up responsibilities. The construction of the expression with the correct linking of parameters and initial values etc has to be performed by some other code. This code is subject to the same issues mentioned in the individual solutions of group 1.

Decision outcome¶

Use a composition based solution from group 2.

Important is the definition of the interface following solution group 2. This ensures to be open-closed and keep the responsibilities separated.

The expertsystem favors composition over inheritance: we intend to use inheritance only to define interfaces, not to insert behavior. As such, the design of the expertsystem is fundamentally different than that of SymPy. That’s another reason to favor composition here: the interfaces are not determined by the dependency and instead remain contained within the dynamics class.

We decide to keep responsibilities as separated as possible. This means that:

The only responsibility of set_dynamics method is to attribute some expression (sympy.Expr) the correct symbol within the complete amplitude model. For now, this position is specified using some StateTransitionGraph and an edge ID, but this syntax may be improved later (see #458):
```
def set_dynamics(
        self,
        graph: StateTransitionGraph,
        edge_id: int,
        expression: sp.Expr,
    ) -> None:
        # dynamics_symbol = graph + edge_id
        # self.dynamics[dynamics_symbol] = expression
```
It is assumed that the expression is correct.
The user has the responsibility of formulating the expression with the correct symbols. To aid the user in the construction of such expressions some building code can handle some of the common tasks, such as
- A VariablePool can facilitate finding the correct symbol names (to avoid typos).
```
mass = variable_pool.get_invariant_mass(graph, edge_id)
```
- A dynamics module provides descriptions of common line-shapes as well as some helper functions. An example would be:
```
inv_mass, mass0, gamma0 = build_relativistic_breit_wigner(
    graph, edge_id, particle
)
rel_bw: sympy.Expr = relativistic_breit_wigner(inv_mass, mass0, gamma0)
model.set_dynamics(graph, edge, rel_bw)
```

The SympyModel has the responsibility of defining a the full model in terms of an expression and keeping track of variables and parameters, for instance:

@attr.s(kw_only=True)
class SympyModel:
    top: sp.Expr = attr.ib()
    # intensities: Dict[sp.Symbol, sp.Expr] = attr.ib(default=dict())
    # amplitudes: Dict[sp.Symbol, sp.Expr] = attr.ib(default=dict())
    dynamics: Dict[sp.Symbol, sp.Expr] = attr.ib(default=dict())
    parameters: Set[sp.Symbol]
    variables: Set[sp.Symbol]  # or: VariablePool

    def full_expression(self) -> sp.Expr:
        ...

For new ADRs, please use adr/template.md as basis. This template was inspired by MADR. General information about architectural decision records is available at adr.github.io.

Welcome to QRules!¶

Internal design¶

1. State Transition Graphs¶

2. Conservation Rules¶

3. Solvers¶

QRules workflow¶

Table of Contents¶

Installation¶

Usage¶

Main interface¶

Investigate intermediate resonances¶

Quantum number search¶

Check allowed reactions¶

Advanced¶

Generate transitions¶

1. Define the problem set¶

2. Prepare Problem Sets¶

3. Find solutions¶

4. Export generated transitions¶

Particle database¶

Loading the default database¶

Finding particles¶

LaTeX representation¶

Adding custom particle definitions through Python¶

Loading custom definitions from a YAML file¶

Writing to YAML¶

Visualize decay topologies¶

Topologies¶

StateTransitionGraphs¶

Bibliography¶

qrules¶

io¶

argument_handling¶

combinatorics¶

conservation_rules¶

default_settings¶

particle¶

quantum_numbers¶

solving¶

topology¶

transition¶

Architectural Decision Records¶

[ADR-000] Use ADRs¶

Context and Problem Statement¶

Considered Options¶

Decision Outcome¶

[ADR-001] Amplitude model¶

Context and problem statement¶

Technical story¶

Decision drivers¶

Solution requirements¶

Considered solutions¶

Customized Python classes¶

Alternative solutions¶

Evaluation¶

Pros and Cons¶

Decision outcome¶

[ADR-002] Inserting dynamics¶

Context and problem statement¶

Technical story¶

Issues with existing set-up¶

Decision drivers¶

Solution requirements¶

Considered solutions¶

Group 1: expression builder¶

Group 2: expression-only¶

Solution evaluation¶

1: Expression builder¶

2: Expression-only¶

Decision outcome¶

`StateTransitionGraph`s¶