ENH: use Fraction for spin values

docs/conf.py

@@ -49,6 +49,7 @@ def pick_newtype_attrs(some_type: type) -> list:
 api_target_substitutions: dict[str, str | tuple[str, str]] = {
     "EdgeQuantumNumberTypes": ("obj", "qrules.quantum_numbers.EdgeQuantumNumberTypes"),
     "EdgeType": "typing.TypeVar",
+    "Fraction": ("obj", "fractions.Fraction"),
     "GraphEdgePropertyMap": ("obj", "qrules.argument_handling.GraphEdgePropertyMap"),
     "GraphElementProperties": ("obj", "qrules.solving.GraphElementProperties"),
     "GraphNodePropertyMap": ("obj", "qrules.argument_handling.GraphNodePropertyMap"),
@@ -66,6 +67,7 @@ def pick_newtype_attrs(some_type: type) -> list:
     "Rule": ("obj", "qrules.argument_handling.Rule"),
     "SpinFormalism": ("obj", "qrules.transition.SpinFormalism"),
     "StateDefinition": ("obj", "qrules.combinatorics.StateDefinition"),
+    "StateDefinitionInput": ("obj", "qrules.combinatorics.StateDefinitionInput"),
     "StateTransition": ("obj", "qrules.transition.StateTransition"),
     "typing.Literal[-1, 1]": "typing.Literal",
 }

src/qrules/__init__.py

@@ -23,7 +23,12 @@
 import attrs
 
 from qrules import io
-from qrules.combinatorics import InitialFacts, StateDefinition, create_initial_facts
+from qrules.combinatorics import (
+    InitialFacts,
+    StateDefinition,
+    StateDefinitionInput,
+    create_initial_facts,
+)
 from qrules.conservation_rules import (
     BaryonNumberConservation,
     BottomnessConservation,
@@ -73,7 +78,7 @@ def check_reaction_violations(  # noqa: C901, PLR0917
     mass_conservation_factor: float | None = 3.0,
     particle_db: ParticleCollection | None = None,
     max_angular_momentum: int = 1,
-    max_spin_magnitude: float = 2.0,
+    max_spin_magnitude: float = 2,
 ) -> set[frozenset[str]]:
     """Determine violated interaction rules for a given particle reaction.
 
@@ -264,15 +269,15 @@ def check_edge_qn_conservation() -> set[frozenset[str]]:
 
 
 def generate_transitions(  # noqa: PLR0917
-    initial_state: StateDefinition | Sequence[StateDefinition],
-    final_state: Sequence[StateDefinition],
+    initial_state: StateDefinitionInput | Sequence[StateDefinitionInput],
+    final_state: Sequence[StateDefinitionInput],
     allowed_intermediate_particles: list[str] | None = None,
     allowed_interaction_types: str | Iterable[str] | None = None,
     formalism: SpinFormalism = "canonical-helicity",
     particle_db: ParticleCollection | None = None,
     mass_conservation_factor: float | None = 3.0,
     max_angular_momentum: int = 2,
-    max_spin_magnitude: float = 2.0,
+    max_spin_magnitude: float = 2,
     topology_building: str = "isobar",
     number_of_threads: int | None = None,
 ) -> ReactionInfo:

src/qrules/argument_handling.py

@@ -8,6 +8,7 @@
 from __future__ import annotations
 
 import inspect
+from fractions import Fraction
 from typing import TYPE_CHECKING, Any, Callable, Generic, TypeVar, Union
 
 import attrs
@@ -22,7 +23,7 @@
 if TYPE_CHECKING:
     from collections.abc import Sequence
 
-Scalar = Union[int, float]
+Scalar = Union[int, float, Fraction]
 
 Rule = Union[GraphElementRule, EdgeQNConservationRule, ConservationRule]
 """Any type of rule"""

src/qrules/combinatorics.py

@@ -11,8 +11,10 @@
 from collections import OrderedDict
 from collections.abc import Iterable, Mapping, Sequence
 from copy import deepcopy
+from fractions import Fraction
 from typing import TYPE_CHECKING, Any, Callable, Union
 
+from qrules.argument_handling import Scalar
 from qrules.particle import ParticleWithSpin
 from qrules.quantum_numbers import InteractionProperties, arange
 from qrules.topology import MutableTransition, Topology, get_originating_node_list
@@ -21,13 +23,27 @@
     from qrules.particle import ParticleCollection
 
 
-StateWithSpins = tuple[str, Sequence[float]]
+StateWithSpins = tuple[str, Sequence[Fraction]]
 StateDefinition = Union[str, StateWithSpins]
 """Particle name, optionally with a list of spin projections."""
+StateDefinitionInput = Union[str, tuple[str, Sequence[Scalar]]]
+"""Input type for `StateDefinition` permitting also `int` and `float`"""
 InitialFacts = MutableTransition[ParticleWithSpin, InteractionProperties]
 """A `.Transition` with only initial and final state information."""
 
 
+def as_state_definition(
+    definition: StateDefinitionInput,
+) -> StateDefinition:
+    if type(definition) is str:
+        return definition
+    if type(definition) is tuple:
+        name, state = definition
+        return name, list(map(Fraction, state))  # type: ignore  # noqa: PGH003
+    msg = f"value has to be of type {StateDefinitionInput}, got {type(definition)}"
+    raise ValueError(msg)
+
+
 class _KinematicRepresentation:  # noqa: PLW1641
     def __init__(
         self,
@@ -182,13 +198,14 @@ def fill_groupings(
 
 def create_initial_facts(
     topology: Topology,
-    initial_state: Sequence[StateDefinition],
-    final_state: Sequence[StateDefinition],
+    initial_state: Sequence[StateDefinitionInput],
+    final_state: Sequence[StateDefinitionInput],
     particle_db: ParticleCollection,
 ) -> list[InitialFacts]:
     states = __create_states_with_spin_projections(
         list(topology.incoming_edge_ids) + list(topology.outgoing_edge_ids),
-        list(initial_state) + list(final_state),
+        list(map(as_state_definition, initial_state))
+        + list(map(as_state_definition, final_state)),
         particle_db,
     )
     spin_states = __generate_spin_combinations(states, particle_db)
@@ -215,9 +232,9 @@ def fill_spin_projections(state: StateDefinition) -> StateWithSpins:
         if isinstance(state, str):
             particle_name = state
             particle = particle_db[particle_name]
-            spin_projections = set(arange(-particle.spin, particle.spin + 1, 1.0))
-            if particle.mass == 0.0 and 0.0 in spin_projections:
-                spin_projections.remove(0.0)
+            spin_projections = set(arange(-particle.spin, particle.spin + 1))
+            if particle.mass == 0.0 and Fraction(0) in spin_projections:
+                spin_projections.remove(Fraction(0))
             return particle_name, sorted(spin_projections)
         return state
 
@@ -256,20 +273,20 @@ def populate_edge_with_spin_projections(
 
 def permutate_topology_kinematically(
     topology: Topology,
-    initial_state: list[StateDefinition],
-    final_state: list[StateDefinition],
+    initial_state: Sequence[StateDefinitionInput] | Sequence[StateDefinition],
+    final_state: Sequence[StateDefinitionInput] | Sequence[StateDefinition],
     final_state_groupings: list[list[list[str]]]
     | list[list[str]]
     | list[str]
     | None = None,
 ) -> list[Topology]:
-    def strip_spin(state: StateDefinition) -> str:
+    def strip_spin(state: StateDefinitionInput) -> str:
         if isinstance(state, tuple):
             return state[0]
         return state
 
     edge_ids = sorted(topology.incoming_edge_ids) + sorted(topology.outgoing_edge_ids)
-    states = initial_state + final_state
+    states = list(initial_state) + list(final_state)
     return _generate_kinematic_permutations(
         topology,
         particle_names={i: strip_spin(s) for i, s in zip(edge_ids, states)},

src/qrules/conservation_rules.py

@@ -47,6 +47,7 @@
 
 import operator
 from copy import deepcopy
+from fractions import Fraction
 from functools import reduce
 from textwrap import dedent
 from typing import Any, Callable, Optional, Protocol, Union
@@ -59,7 +60,7 @@
 from qrules.quantum_numbers import arange
 
 
-def _is_boson(spin_magnitude: float) -> bool:
+def _is_boson(spin_magnitude: Fraction) -> bool:
     return abs(spin_magnitude % 1) < 0.01
 
 
@@ -235,6 +236,7 @@ class CParityEdgeInput:
 
 @frozen
 class CParityNodeInput:
+    # These converters currently do not do anything, as "NewType"s do not have constructors
     l_magnitude: NodeQN.l_magnitude = field(converter=NodeQN.l_magnitude)
     s_magnitude: NodeQN.s_magnitude = field(converter=NodeQN.s_magnitude)
 
@@ -264,8 +266,8 @@ def _get_c_parity_multiparticle(
                 # if boson
                 if _is_boson(part_qns[0].spin_magnitude):
                     return (-1) ** int(ang_mom)
-                coupled_spin = interaction_qns.s_magnitude
-                if isinstance(coupled_spin, int) or coupled_spin.is_integer():
+                coupled_spin = Fraction(interaction_qns.s_magnitude)
+                if isinstance(coupled_spin, int) or coupled_spin.denominator == 1:
                     return (-1) ** int(ang_mom + coupled_spin)
         return None
 
@@ -314,12 +316,12 @@ def check_multistate_g_parity(
             double_state_qns[0].pid, double_state_qns[1].pid
         ):
             ang_mom = interaction_qns.l_magnitude
-            if isinstance(isospin, int) or isospin.is_integer():
+            if isinstance(isospin, int) or isospin.denominator == 1:
                 # if boson
                 if _is_boson(double_state_qns[0].spin_magnitude):
                     return (-1) ** int(ang_mom + isospin)
                 coupled_spin = interaction_qns.s_magnitude
-                if isinstance(coupled_spin, int) or coupled_spin.is_integer():
+                if isinstance(coupled_spin, int) or coupled_spin.denominator == 1:
                     return (-1) ** int(ang_mom + coupled_spin + isospin)
         return None
 
@@ -429,8 +431,8 @@ def _check_particles_identical(
 
 @frozen
 class _Spin:
-    magnitude: float
-    projection: float
+    magnitude: Union[Fraction, NodeQN.s_magnitude, NodeQN.l_magnitude]
+    projection: Union[Fraction, NodeQN.s_projection, NodeQN.l_projection]
 
 
 def _is_clebsch_gordan_coefficient_zero(
@@ -468,26 +470,25 @@ class SpinMagnitudeNodeInput:
 def ls_spin_validity(spin_input: SpinNodeInput) -> bool:
     r"""Check for valid isospin magnitude and projection."""
     return _check_spin_valid(
-        float(spin_input.l_magnitude), float(spin_input.l_projection)
-    ) and _check_spin_valid(
-        float(spin_input.s_magnitude), float(spin_input.s_projection)
-    )
+        spin_input.l_magnitude, spin_input.l_projection
+    ) and _check_spin_valid(spin_input.s_magnitude, spin_input.s_projection)
 
 
 def _check_magnitude(
-    in_part: list[float],
-    out_part: list[float],
+    in_part: list[Fraction],
+    out_part: list[Fraction],
     interaction_qns: Optional[Union[SpinMagnitudeNodeInput, SpinNodeInput]],
 ) -> bool:
-    def couple_mags(j_1: float, j_2: float) -> list[float]:
+    def couple_mags(j_1: Fraction, j_2: Fraction) -> list[Fraction]:
         return [
-            x / 2.0 for x in range(int(2 * abs(j_1 - j_2)), int(2 * (j_1 + j_2 + 1)), 2)
+            Fraction(x, 2)
+            for x in range(int(2 * abs(j_1 - j_2)), int(2 * (j_1 + j_2 + 1)), 2)
         ]
 
     def couple_magnitudes(
-        magnitudes: list[float],
+        magnitudes: list[Fraction],
         interaction_qns: Optional[Union[SpinMagnitudeNodeInput, SpinNodeInput]],
-    ) -> set[float]:
+    ) -> set[Fraction]:
         if len(magnitudes) == 1:
             return set(magnitudes)
 
@@ -570,13 +571,13 @@ def __spin_couplings(spin1: _Spin, spin2: _Spin) -> set[_Spin]:
 
     :math:`|S_1 - S_2| \leq S \leq |S_1 + S_2|` and :math:`M_1 + M_2 = M`
     """
-    s_1 = spin1.magnitude
-    s_2 = spin2.magnitude
-
     sum_proj = spin1.projection + spin2.projection
     return {
-        _Spin(x, sum_proj)
-        for x in arange(abs(s_1 - s_2), s_1 + s_2 + 1, 1.0)
+        _Spin(Fraction(x), Fraction(sum_proj))
+        for x in arange(
+            abs(spin1.magnitude - spin2.magnitude),
+            spin1.magnitude + spin2.magnitude + 1,
+        )
         if x >= abs(sum_proj)
         and not _is_clebsch_gordan_coefficient_zero(spin1, spin2, _Spin(x, sum_proj))
     }
@@ -592,19 +593,17 @@ class IsoSpinEdgeInput:
     )
 
 
-def _check_spin_valid(magnitude: float, projection: float) -> bool:
-    if magnitude % 0.5 != 0.0:
+def _check_spin_valid(magnitude: Fraction, projection: Fraction) -> bool:
+    if magnitude.denominator not in {1, 2}:
         return False
     if abs(projection) > magnitude:
         return False
-    return float(projection - magnitude).is_integer()
+    return (projection - magnitude).denominator == 1
 
 
 def isospin_validity(isospin: IsoSpinEdgeInput) -> bool:
     r"""Check for valid isospin magnitude and projection."""
-    return _check_spin_valid(
-        float(isospin.isospin_magnitude), float(isospin.isospin_projection)
-    )
+    return _check_spin_valid(isospin.isospin_magnitude, isospin.isospin_projection)
 
 
 def isospin_conservation(
@@ -642,7 +641,7 @@ class SpinEdgeInput:
 
 def spin_validity(spin: SpinEdgeInput) -> bool:
     r"""Check for valid spin magnitude and projection."""
-    return _check_spin_valid(float(spin.spin_magnitude), float(spin.spin_projection))
+    return _check_spin_valid(spin.spin_magnitude, spin.spin_projection)
 
 
 def spin_conservation(
@@ -710,8 +709,8 @@ def spin_magnitude_conservation(
         len(ingoing_spin_magnitudes) == 2 and len(outgoing_spin_magnitudes) == 1
     ):
         return _check_magnitude(
-            [float(x) for x in ingoing_spin_magnitudes],
-            [float(x) for x in outgoing_spin_magnitudes],
+            [Fraction(x) for x in ingoing_spin_magnitudes],
+            [Fraction(x) for x in outgoing_spin_magnitudes],
             interaction_qns,
         )
 
@@ -856,7 +855,7 @@ def calculate_hypercharge(
         or edge_qns.tau_lepton_number
     ):
         return True
-    isospin_3 = 0.0
+    isospin_3 = Fraction(0)
     if edge_qns.isospin_projection:
         isospin_3 = edge_qns.isospin_projection
     return float(edge_qns.charge) == isospin_3 + 0.5 * calculate_hypercharge(edge_qns)

src/qrules/io/_dict.py

@@ -4,6 +4,7 @@
 
 import json
 from collections import abc
+from fractions import Fraction
 from os.path import dirname, realpath
 from typing import Any
 
@@ -28,7 +29,7 @@ def from_attrs_decorated(inst: Any) -> dict:
     )
 
 
-def _value_serializer(inst: type, field: attrs.Attribute, value: Any) -> Any:  # noqa: ARG001
+def _value_serializer(inst: type, field: attrs.Attribute, value: Any) -> Any:  # noqa: ARG001, PLR0911
     if isinstance(value, abc.Mapping):
         if all(isinstance(p, Particle) for p in value.values()):
             return {k: v.name for k, v in value.items()}
@@ -43,6 +44,8 @@ def _value_serializer(inst: type, field: attrs.Attribute, value: Any) -> Any:  #
             "magnitude": value.magnitude,
             "projection": value.projection,
         }
+    if isinstance(value, Fraction):
+        return float(value)
     return value