Add functionality to calculate delta values (#20)

* Processor utility contains a delta function calculator

* Finish up `integrals_calc_delta` in processor

* Add integral delta calculation per default to multi processor

* Lint

* Add docs for delta values

* Add a run through test for delta calculation

Author: trappitsch

docs/bg/deltas.rst

@@ -0,0 +1,44 @@
+============
+Delta Values
+============
+
+For all integrals, :math:`\delta`-values can be calculated automatically
+as described in the :doc:`usage documentation <../pkg/usage>`.
+Delta values are always calculated (1) with respect to the major isotope
+and (2) using the NIST database of the
+`iniabu <https://github.com/galactic-forensics/iniabu>`_ package.
+If no value can be calculated, i.e.,
+because the isotope is unknown to ``iniabu`` or because
+the natural abundance of the specified isotope is unknown,
+``np.nan`` will be returned as the :math:`\delta`-value.
+
+.. note:: All :math:`\delta`-values are reported in per mil (‰).
+
+--------------------------------
+Calculation :math:`\delta`-value
+--------------------------------
+
+Be :math:`i` and :math:`j` the number of counts
+in the nominator and denominator isotopes, respectively, and
+:math:`r` the NIST isotopic ratio of the same isotopes.
+The :math:`\delta`-value can then be calculated as:
+
+.. math::
+
+    \delta = \left( \frac{i/j}{r} - 1 \right) \times 1000 \qquad (‰)
+
+-------------
+Uncertainties
+-------------
+
+Be :math:`\sigma_i` and :math:`\sigma_j` the uncertainties of values
+:math:`i` and :math:`j`, respectively.
+The uncertainty of the :math:`\delta`-value :math:`\sigma_{\delta}`
+can be calculated as:
+
+.. math::
+
+    \sigma_{\delta} = \frac{1000}{r} \left[
+        \left(\frac{\sigma_i}{j}\right)^2 +
+        \left(\frac{i\sigma_j}{j^2}\right)^2 \right]^{1/2}
+        \qquad (‰)

docs/index.rst

@@ -71,6 +71,7 @@ Contents
     :maxdepth: 2
     :caption: Scientific Background
 
+    bg/deltas
     bg/integrals
 
 .. toctree::

docs/pkg/usage.rst

@@ -169,6 +169,9 @@ you could set the backgrounds as following:
 .. warning:: The backgrounds you define must have the same name as the peaks they are defined for.
     Multiple definitions per background can exist.
 
+.. note:: Integral and background definitions can also be performed using a ``matplotlib`` GUI.
+    For details, see :doc:`here <guis>`.
+
 To apply the integrals,
 simply run:
 
@@ -181,8 +184,23 @@ to either ``True`` or ``False``, respectively.
 Details on integral background corrections and the math behind it
 can be found :doc:`here <../bg/integrals>`.
 
-.. note:: Integral and background definitions can also be performed using a ``matplotlib`` GUI.
-    For details, see :doc:`here <guis>`.
+Finally, if your integral names follow the format used in the ``iniabu`` package,
+you can calculate :math:`\delta`-values for your individual peaks automatically.
+These values are always calculated with respect to the major isotope.
+If values are not available, ``np.nan`` will be written for that specific
+:math:`\delta`-value.
+To calculate the :math:`\delta` values, run (after calculating integrals):
+
+.. code-block:: python
+
+    crd.integrals_calc_delta()
+
+This will save the :math:`\delta`-values and associated uncertainties to
+``crd.integrals_delta``. If packages were defined, an :math:`\delta`-values
+for each package will also be calculated and stored in ``crd.integrals_delta_pkg``.
+Details on the calculation and error propagation can be found
+:doc:`here <../bg/deltas>`.
+
 
 +++++++
 Results

docs/pkg/variables.rst

@@ -55,11 +55,15 @@ Integrals
     the sum of counts in the defined area and its uncertainty.
     :math:`m` defined peaks, it will therefore be of size :math:`m \times 2`.
     The second entry of each integral is the uncertainty of the counts.
+- ``crd.integrals_delta``: :math:`\delta`-values for all integrals and uncertainties.
+    The format is the same as for ``crd.integrals``.
 - ``crd.integrals_pkg``: Integral data for packages.
     For a total of :math:`p` packages and :math:`m` defined integrals,
     this ``numpy.ndarray`` will be of size :math:`p \times m \times 2`.
     As for the integrals, the sum of counts for each peak and its uncertainty
     are given.
+- ``crd.integrals_delta_pkg``: :math:`\delta`-values for all packages.
+    The format is the same as for ``crd.integrals_pkg``.
 
 ----------------
 Single Shot Data

rimseval/multi_proc.py

@@ -103,6 +103,7 @@ def apply_to_all(self, id: int, opt_mcal: bool = False, bg_corr: bool = False):
         else:
             bg_corr = False
         crd_main.integrals_calc(bg_corr=bg_corr)
+        crd_main.integrals_calc_delta()
         self.signal_processed.emit(str(crd_main.fname.name))
         rimseval.interfacer.save_cal_file(crd_main)
 
@@ -126,6 +127,7 @@ def apply_to_all(self, id: int, opt_mcal: bool = False, bg_corr: bool = False):
 
                 # integrals
                 file.integrals_calc(bg_corr=bg_corr)
+                file.integrals_calc_delta()
 
                 # save calibration
                 rimseval.interfacer.save_cal_file(file)

rimseval/processor.py

@@ -56,7 +56,9 @@ def __init__(self, fname: Path) -> None:
 
         # Integrals
         self.integrals = None
+        self.integrals_delta = None
         self.integrals_pkg = None
+        self.integrals_delta_pkg = None
 
         # parameters for calibration and evaluation
         self._params_mcal = None  # mass calibration
@@ -600,6 +602,35 @@ def integral_windows(limits_tmp: np.array) -> List:
                 bgs_pkg,
             )
 
+    def integrals_calc_delta(self) -> None:
+        """Calculate delta values for integrals and save them in class.
+
+        This routine uses the ``iniabu`` package to calculate delta values for defined
+        integrals. It reads the peak names and calculates delta values for isotopes
+        that can be understood ``iniabu``, and calculates the delta values with
+        respect to the major isotope. These values are then saved to the class as
+        ``integrals_delta`` and ``integrals_delta_pkg``, if packages were defined.
+        Uncertainties are propagated according to Gaussian error propagation.
+        The format of the resulting arrays are identical to the ``integrals`` and
+        ``integrals_pkg`` arrays.
+
+        :raises ValueError: No integrals were calculated.
+        """
+        if self.integrals is None or self.def_integrals is None:
+            raise ValueError("No integrals were defined or calculated.")
+
+        peak_names = self.def_integrals[0]
+
+        integrals_delta = processor_utils.delta_calc(peak_names, self.integrals)
+
+        if self.integrals_pkg is not None:
+            integrals_delta_pkg = np.zeros_like(self.integrals_pkg, dtype=float)
+            for it, line in enumerate(self.integrals_pkg):
+                integrals_delta_pkg[it] = processor_utils.delta_calc(peak_names, line)
+            self.integrals_delta_pkg = integrals_delta_pkg
+
+        self.integrals_delta = integrals_delta
+
     def mass_calibration(self) -> None:
         r"""Perform a mass calibration on the data.
 

rimseval/processor_utils.py

@@ -1,12 +1,12 @@
 """Utilities for CRD processors. Mostly methods that can be jitted."""
 
-from typing import Tuple, Union
+from typing import List, Tuple, Union
 
 from numba import njit
 import numpy as np
 from scipy import optimize
 
-from .utilities import fitting, utils
+from .utilities import fitting, ini, utils
 
 
 @njit
@@ -82,6 +82,72 @@ def dead_time_correction(
     return data
 
 
+def delta_calc(names: List[str], integrals: np.ndarray) -> np.ndarray:
+    """Calculate delta values for a given set of integrals.
+
+    Use ``iniabu`` to calculate the delta values with respect to major isotope.
+    If the name of a peak is not valid or the major isotope not present, return
+    ``np.nan`` for that entry. Appropriate error propagation is done as well.
+
+    :param names: Names of the peaks as list.
+    :param integrals: Integrals, as defined in ``CRDFileProcessor.integrals``.
+
+    :return: List of delta values, same shape and format as ``integrals``.
+    """
+    # transform all names to valid ``iniabu`` names or call them ``None``
+    names_iniabu = []
+    for name in names:
+        try:
+            names_iniabu.append(ini.iso[name].name)
+        except IndexError:
+            names_iniabu.append(None)
+
+    # find major isotope names
+    major_iso_name = []
+    for name in names_iniabu:
+        if name is None:
+            major_iso_name.append(None)
+        else:
+            ele = name.split("-")[0]
+            try:
+                maj = ini._get_major_iso(ele)
+            except IndexError:
+                maj = None
+            major_iso_name.append(maj)
+
+    integrals_dict = dict(zip(names_iniabu, range(len(names_iniabu))))
+
+    integrals_delta = np.zeros_like(integrals, dtype=float)
+
+    for it, iso in enumerate(names_iniabu):
+        maj_iso = major_iso_name[it]
+
+        if iso is None or maj_iso not in names_iniabu:
+            integrals_delta[it][0] = np.nan
+            integrals_delta[it][1] = np.nan
+        else:
+            msr_nom = integrals[it][0]
+            msr_nom_unc = integrals[it][1]
+            msr_denom = integrals[integrals_dict[maj_iso]][0]
+            msr_denom_unc = integrals[integrals_dict[maj_iso]][1]
+
+            msr_ratio = msr_nom / msr_denom
+            integrals_delta[it][0] = ini.iso_delta(iso, maj_iso, msr_ratio)
+
+            # error calculation
+            std_ratio = ini.iso_ratio(iso, maj_iso)
+            integrals_delta[it][1] = (
+                1000
+                / std_ratio
+                * np.sqrt(
+                    (msr_nom_unc / msr_denom) ** 2
+                    + (msr_nom * msr_denom_unc / msr_denom**2) ** 2
+                )
+            )
+
+    return integrals_delta
+
+
 def gaussian_fit_get_max(xdata: np.ndarray, ydata: np.ndarray) -> float:
     """Fit a Gaussian to xdata and ydata and return the xvalue of the peak.
 

tests/func/conftest.py

@@ -0,0 +1,49 @@
+"""Fixtures for functional tests."""
+
+from pathlib import Path
+
+import pytest
+import numpy as np
+
+from rimseval.processor import CRDFileProcessor
+
+
+@pytest.fixture
+def crd_int_delta(crd_file) -> CRDFileProcessor:
+    """Initialize a CRDFileProcessor file, set some integrals, and return it.
+
+    The integrals set are fake, no actual calculation takes place. This fixture simply
+    serves to test functionality (i.e., delta calculation) after processing of the CRD
+    file.
+
+    :return: CRDFileProcessor instance with integrals defined and set.
+    """
+    _, _, _, fname = crd_file
+    crd = CRDFileProcessor(Path(fname))
+
+    peak_names = ["54Fe", "56Fe", "57Fe", "58Fe", "244Pu", "238Pu", "bg"]
+    peak_limits = np.array(
+        [
+            [53.8, 54.2],
+            [55.8, 56.2],
+            [56.8, 57.2],
+            [57.8, 58.2],
+            [243.8, 244.2],
+            [237.8, 238.2],
+            [238.5, 238.7],
+        ]
+    )
+    crd.def_integrals = peak_names, peak_limits
+    crd.integrals = np.array(
+        [
+            [5.84500000e04, 2.41764348e02],
+            [9.17540000e05, 9.57883083e02],
+            [2.11900000e04, 1.45567854e02],
+            [2.82000000e03, 5.31036722e01],
+            [10000, 100],
+            [144, 12],
+            [34212, 185],
+        ]
+    )
+    crd.integrals_pkg = np.array([crd.integrals, crd.integrals])
+    return crd

tests/func/test_processor.py

@@ -186,6 +186,16 @@ def test_integrals_definition_delete_undefined_background(crd_file):
     assert "Int2" not in crd.def_backgrounds[0]
 
 
+def test_integrals_delta_calc(crd_int_delta):
+    """Assure that the delta calculation on the integral returns correct values."""
+    crd_int_delta.integrals_calc_delta()
+    deltas_exp = np.array([0, 0, 0, 0, np.nan, np.nan, np.nan])
+    deltas_rec = crd_int_delta.integrals_delta.transpose()[0]  # test only deltas
+    np.testing.assert_almost_equal(deltas_exp, deltas_rec)
+    for line in crd_int_delta.integrals_pkg:
+        np.testing.assert_almost_equal(crd_int_delta.integrals, line)
+
+
 def test_packages(crd_file):
     """Simple test to ensure packages are made in two ways correctly."""
     _, ions_per_shot, _, fname = crd_file

tests/func/test_processor_utils.py

@@ -43,6 +43,17 @@ def test_create_packages():
     np.testing.assert_equal(pkg_data_rec, pkg_data_exp)
 
 
+def test_delta_calc():
+    """Take an integrals like array and calculate delta values.
+
+    Details of delta calculation is not tested.
+    """
+    names = ["Fe54", "Fe56", "244Pu", "bg"]
+    integrals = np.array([[10000, 100], [100000, 240], [100, 10], [2001, 21]])
+    deltas = pu.delta_calc(names, integrals)
+    assert np.isnan(deltas[2:3]).all()  # last two must be nans
+
+
 def test_integrals_bg_corr():
     """Background correction for defined integrals."""
     integrals = np.array([[10, np.sqrt(10)], [40, np.sqrt(40)]])