From 443c2ac4c99c06866798ca14b2944d8d50e0141a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 24 May 2022 11:59:30 -0500
Subject: [PATCH 01/12] Add integral method on Discrete, Tabular, and Mixture
 classes

---
 openmc/stats/univariate.py | 39 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 5528bdb30..733b7eaa6 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -220,6 +220,16 @@ class Discrete(Univariate):
         p_arr = np.array([p_merged[x] for x in x_arr])
         return cls(x_arr, p_arr)
 
+    def integral(self):
+        """Return integral of distribution
+
+        Returns
+        -------
+        float
+            Integral of discrete distribution
+        """
+        return np.sum(self.p)
+
 class Uniform(Univariate):
     """Distribution with constant probability over a finite interval [a,b]
 
@@ -1027,6 +1037,22 @@ class Tabular(Univariate):
         p = params[len(params)//2:]
         return cls(x, p, interpolation)
 
+    def integral(self):
+        """Return integral of distribution
+
+        Returns
+        -------
+        float
+            Integral of tabular distrbution
+        """
+        if self.interpolation == 'histogram':
+            return np.sum(np.diff(self.x) * self.p[:-1])
+        elif self.interpolation == 'linear-linear':
+            return np.trapz(self.p, self.x)
+        else:
+            raise NotImplementedError(
+                f'integral() not supported for {self.inteprolation} interpolation')
+
 
 class Legendre(Univariate):
     r"""Probability density given by a Legendre polynomial expansion
@@ -1199,3 +1225,16 @@ class Mixture(Univariate):
             distribution.append(Univariate.from_xml_element(pair.find("dist")))
 
         return cls(probability, distribution)
+
+    def integral(self):
+        """Return integral of the distribution
+
+        Returns
+        -------
+        float
+            Integral of the distribution
+        """
+        return sum([
+            p*dist.integral()
+            for p, dist in zip(self.probability, self.distribution)
+        ])

From c4a9b2c8b13592e2e761985eef3f2eef3dc092f0 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 5 Apr 2022 11:34:09 -0500
Subject: [PATCH 02/12] Get gamma and xray decay sources as distributions

---
 openmc/data/decay.py | 42 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 51a594bc8..2855a9434 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -9,7 +9,9 @@ from uncertainties import ufloat, UFloat
 
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
+from openmc.stats import Discrete, Tabular
 from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
+from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
 
 
@@ -495,3 +497,43 @@ class Decay(EqualityMixin):
 
         """
         return cls(ev_or_filename)
+
+    def get_sources(self):
+        sources = {}
+        name = self.nuclide['name']
+        for particle, spectra in self.spectra.items():
+            # Only handle gammas for now
+            if particle not in ('gamma', 'xray'):
+                continue
+
+            # Create distribution for discrete
+            distributions = []
+            if spectra['continuous_flag'] in ('discrete', 'both'):
+                energies = []
+                intensities = []
+                for discrete_data in spectra['discrete']:
+                    energies.append(discrete_data['energy'].n)
+                    intensities.append(discrete_data['intensity'].n)
+                energies = np.array(energies)
+                intensities = np.array(intensities)
+                dist_discrete = Discrete(energies, intensities)  # <-- not normalized yet
+                dist_discrete._normalization = spectra['discrete_normalization'].n
+                distributions.append(dist_discrete)
+
+            # Create distribution for continuous
+            if spectra['continuous_flag'] in ('continuous', 'both'):
+                f = spectra['continuous']['probability']
+                if len(f.interpolation) > 1:
+                    raise NotImplementedError("Multiple interpolation regions: {name}, {particle}")
+                interpolation = INTERPOLATION_SCHEME[f.interpolation[0]]
+                if interpolation not in ('histogram', 'linear-linear'):
+                    raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported")
+
+                dist_continuous = Tabular(f.x, f.y, interpolation)
+                dist_continuous._intensity = spectra['continuous_normalization'].n
+                distributions.append(dist_continuous)
+
+            # Combine distribution for discrete and continuous
+            sources[particle] = distributions
+
+        return sources

From b7596245f2af6c5484c9bb73e323e39ee6fb2de9 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 5 Apr 2022 12:50:22 -0500
Subject: [PATCH 03/12] Combined decay distributions

---
 openmc/data/decay.py | 46 +++++++++++++++++++++++++++++++++++++-------
 1 file changed, 39 insertions(+), 7 deletions(-)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 2855a9434..5de07f465 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -9,7 +9,7 @@ from uncertainties import ufloat, UFloat
 
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
-from openmc.stats import Discrete, Tabular
+from openmc.stats import Discrete, Tabular, Mixture
 from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
@@ -499,15 +499,27 @@ class Decay(EqualityMixin):
         return cls(ev_or_filename)
 
     def get_sources(self):
+        """Get radioactive decay source distributions
+
+        Returns
+        -------
+        sources : dict
+            Dictionary mapping particle types (e.g., 'photon') to instances of
+            :class:`openmc.stats.Univariate`
+
+        """
         sources = {}
         name = self.nuclide['name']
         for particle, spectra in self.spectra.items():
             # Only handle gammas for now
             if particle not in ('gamma', 'xray'):
                 continue
+            # TODO: Set particle type based on 'particle' above
+            particle_type = 'photon'
+            if particle_type not in sources:
+                sources[particle_type] = []
 
             # Create distribution for discrete
-            distributions = []
             if spectra['continuous_flag'] in ('discrete', 'both'):
                 energies = []
                 intensities = []
@@ -516,9 +528,10 @@ class Decay(EqualityMixin):
                     intensities.append(discrete_data['intensity'].n)
                 energies = np.array(energies)
                 intensities = np.array(intensities)
+                intensities *= spectra['discrete_normalization'].n
                 dist_discrete = Discrete(energies, intensities)  # <-- not normalized yet
-                dist_discrete._normalization = spectra['discrete_normalization'].n
-                distributions.append(dist_discrete)
+                dist_discrete._intensity = intensities.sum()
+                sources[particle_type].append(dist_discrete)
 
             # Create distribution for continuous
             if spectra['continuous_flag'] in ('continuous', 'both'):
@@ -531,9 +544,28 @@ class Decay(EqualityMixin):
 
                 dist_continuous = Tabular(f.x, f.y, interpolation)
                 dist_continuous._intensity = spectra['continuous_normalization'].n
-                distributions.append(dist_continuous)
+                sources[particle_type].append(dist_continuous)
 
-            # Combine distribution for discrete and continuous
-            sources[particle] = distributions
+        # Combine discrete distributions
+        for dist_list in sources.values():
+            # Get list of discrete distributions
+            discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
+            dist_discrete = [dist_list[i] for i in discrete_index]
 
+            if len(dist_discrete) > 1:
+                # Create combined discrete distribution
+                probs = [1.0] * len(dist_discrete)
+                combined_dist = Discrete.merge(dist_discrete, probs)
+                combined_dist._intensity = np.sum(combined_dist.p)
+
+                for idx in reversed(discrete_index):
+                    dist_list.pop(idx)
+                dist_list.append(combined_dist)
+
+            # Combine discrete and continuous if present
+            if len(dist_list) > 1:
+                probs = [d._intensity for d in dist_list]
+                dist_list[:] = Mixture(probs, dist_list.copy())
+
+        sources = {k: (v[0] if v else None) for k, v in sources.items()}
         return sources

From 1b2b5e169c7dcf9d504ed505edc76f96ee25e87d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 23 May 2022 14:44:08 -0500
Subject: [PATCH 04/12] Map radiation type to particle type properly in
 Decay.get_sources

---
 openmc/data/decay.py | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 5de07f465..142be35e4 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -511,11 +511,21 @@ class Decay(EqualityMixin):
         sources = {}
         name = self.nuclide['name']
         for particle, spectra in self.spectra.items():
-            # Only handle gammas for now
-            if particle not in ('gamma', 'xray'):
-                continue
-            # TODO: Set particle type based on 'particle' above
-            particle_type = 'photon'
+            # Set particle type based on 'particle' above
+            particle_type = {
+                'gamma': 'photon',
+                'beta-': 'electron',
+                'ec/beta+': 'positron',
+                'alpha': 'alpha',
+                'n': 'neutron',
+                'sf': 'fragment',
+                'p': 'proton',
+                'e-': 'electron',
+                'xray': 'photon',
+                'anti-neutrino': 'anti-neutrino',
+                'neutrino': 'neutrino',
+            }[particle]
+
             if particle_type not in sources:
                 sources[particle_type] = []
 

From 6388c8f994b6587cd2d0de1e3aea2d6c3b514ec6 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 24 May 2022 12:02:31 -0500
Subject: [PATCH 05/12] Separate out a function for combine_distributions

---
 openmc/data/decay.py | 73 +++++++++++++++++++++++++++++++++-----------
 1 file changed, 55 insertions(+), 18 deletions(-)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 142be35e4..edf9c1c0b 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -1,4 +1,5 @@
 from collections.abc import Iterable
+from copy import deepcopy
 from io import StringIO
 from math import log
 import re
@@ -552,30 +553,66 @@ class Decay(EqualityMixin):
                 if interpolation not in ('histogram', 'linear-linear'):
                     raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported")
 
+                # TODO: work normalization into tabular itself
                 dist_continuous = Tabular(f.x, f.y, interpolation)
                 dist_continuous._intensity = spectra['continuous_normalization'].n
                 sources[particle_type].append(dist_continuous)
 
         # Combine discrete distributions
-        for dist_list in sources.values():
-            # Get list of discrete distributions
-            discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
-            dist_discrete = [dist_list[i] for i in discrete_index]
+        merged_sources = {}
+        for particle_type, dist_list in sources.items():
+            merged_sources[particle_type] = combine_distributions(
+                dist_list, [1.0]*len(dist_list))
 
-            if len(dist_discrete) > 1:
-                # Create combined discrete distribution
-                probs = [1.0] * len(dist_discrete)
-                combined_dist = Discrete.merge(dist_discrete, probs)
-                combined_dist._intensity = np.sum(combined_dist.p)
+        return merged_sources
 
-                for idx in reversed(discrete_index):
-                    dist_list.pop(idx)
-                dist_list.append(combined_dist)
 
-            # Combine discrete and continuous if present
-            if len(dist_list) > 1:
-                probs = [d._intensity for d in dist_list]
-                dist_list[:] = Mixture(probs, dist_list.copy())
+def combine_distributions(dists, probs):
+    """Combine distributions with specified probabilities
 
-        sources = {k: (v[0] if v else None) for k, v in sources.items()}
-        return sources
+    This function can be used to combine multiple instances of
+    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular` into a single
+    distribution. Multiple discrete distributions are merged into a single
+    distribution and the remainder of the distributions are put into a
+    :class:`~openmc.stats.Mixture` distribution.
+
+    Parameters
+    ----------
+    dists : iterable of openmc.stats.Univariate
+        Distributions to combine
+    probs : iterable of float
+        Probability (or intensity) of each distribution
+
+    """
+    # Get copy of distribution list so as not to modify the argument
+    dist_list = deepcopy(dists)
+
+    # Get list of discrete/continuous distribution indices
+    discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
+    cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
+
+    # Apply probabilites to continuous distributions
+    for i in cont_index:
+        dist = dist_list[i]
+        dist.p *= probs[i]
+        dist._intensity *= probs[i]
+
+    if discrete_index:
+        # Create combined discrete distribution
+        dist_discrete = [dist_list[i] for i in discrete_index]
+        discrete_probs = [probs[i] for i in discrete_index]
+        combined_dist = Discrete.merge(dist_discrete, discrete_probs)
+        combined_dist._intensity = np.sum(combined_dist.p)
+
+        # Replace multiple discrete distributions with merged
+        for idx in reversed(discrete_index):
+            dist_list.pop(idx)
+        dist_list.append(combined_dist)
+
+    # Combine discrete and continuous if present
+    if len(dist_list) > 1:
+        probs = [d._intensity for d in dist_list]
+        dist_list[:] = [Mixture(probs, dist_list.copy())]
+        dist_list[0]._intensity = sum(probs)
+
+    return dist_list[0]

From d4989ae64226f0e11d75220277bfb13e3bdb538f Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 24 May 2022 12:48:06 -0500
Subject: [PATCH 06/12] Use integral() methods instead of _intensity hidden
 attribute

---
 openmc/data/decay.py | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index edf9c1c0b..874981e2b 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -511,6 +511,7 @@ class Decay(EqualityMixin):
         """
         sources = {}
         name = self.nuclide['name']
+        decay_constant = self.decay_constant.n
         for particle, spectra in self.spectra.items():
             # Set particle type based on 'particle' above
             particle_type = {
@@ -538,10 +539,9 @@ class Decay(EqualityMixin):
                     energies.append(discrete_data['energy'].n)
                     intensities.append(discrete_data['intensity'].n)
                 energies = np.array(energies)
-                intensities = np.array(intensities)
-                intensities *= spectra['discrete_normalization'].n
-                dist_discrete = Discrete(energies, intensities)  # <-- not normalized yet
-                dist_discrete._intensity = intensities.sum()
+                intensity = spectra['discrete_normalization'].n
+                rates = decay_constant * intensity * np.array(intensities)
+                dist_discrete = Discrete(energies, rates)
                 sources[particle_type].append(dist_discrete)
 
             # Create distribution for continuous
@@ -553,9 +553,9 @@ class Decay(EqualityMixin):
                 if interpolation not in ('histogram', 'linear-linear'):
                     raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported")
 
-                # TODO: work normalization into tabular itself
-                dist_continuous = Tabular(f.x, f.y, interpolation)
-                dist_continuous._intensity = spectra['continuous_normalization'].n
+                intensity = spectra['continuous_normalization'].n
+                rates = decay_constant * intensity * f.y
+                dist_continuous = Tabular(f.x, rates, interpolation)
                 sources[particle_type].append(dist_continuous)
 
         # Combine discrete distributions
@@ -595,14 +595,12 @@ def combine_distributions(dists, probs):
     for i in cont_index:
         dist = dist_list[i]
         dist.p *= probs[i]
-        dist._intensity *= probs[i]
 
     if discrete_index:
         # Create combined discrete distribution
         dist_discrete = [dist_list[i] for i in discrete_index]
         discrete_probs = [probs[i] for i in discrete_index]
         combined_dist = Discrete.merge(dist_discrete, discrete_probs)
-        combined_dist._intensity = np.sum(combined_dist.p)
 
         # Replace multiple discrete distributions with merged
         for idx in reversed(discrete_index):
@@ -611,8 +609,7 @@ def combine_distributions(dists, probs):
 
     # Combine discrete and continuous if present
     if len(dist_list) > 1:
-        probs = [d._intensity for d in dist_list]
+        probs = [d.integral() for d in dist_list]
         dist_list[:] = [Mixture(probs, dist_list.copy())]
-        dist_list[0]._intensity = sum(probs)
 
     return dist_list[0]

From 7c779f73f8aa70c43283e4bbed2404ece2b2ab3a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 18 Jul 2022 07:18:23 -0500
Subject: [PATCH 07/12] Move combine_distributions to univariate.py

---
 docs/source/pythonapi/data.rst |  1 +
 openmc/data/decay.py           | 49 +---------------------------------
 openmc/stats/univariate.py     | 49 ++++++++++++++++++++++++++++++++++
 3 files changed, 51 insertions(+), 48 deletions(-)

diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index 287738774..6f56938b8 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -61,6 +61,7 @@ Core Functions
 
     atomic_mass
     atomic_weight
+    combine_distributions
     decay_constant
     dose_coefficients
     gnd_name
diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 874981e2b..f358dadf7 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -1,5 +1,4 @@
 from collections.abc import Iterable
-from copy import deepcopy
 from io import StringIO
 from math import log
 import re
@@ -10,7 +9,7 @@ from uncertainties import ufloat, UFloat
 
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
-from openmc.stats import Discrete, Tabular, Mixture
+from openmc.stats import Discrete, Tabular, combine_distributions
 from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
@@ -567,49 +566,3 @@ class Decay(EqualityMixin):
         return merged_sources
 
 
-def combine_distributions(dists, probs):
-    """Combine distributions with specified probabilities
-
-    This function can be used to combine multiple instances of
-    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular` into a single
-    distribution. Multiple discrete distributions are merged into a single
-    distribution and the remainder of the distributions are put into a
-    :class:`~openmc.stats.Mixture` distribution.
-
-    Parameters
-    ----------
-    dists : iterable of openmc.stats.Univariate
-        Distributions to combine
-    probs : iterable of float
-        Probability (or intensity) of each distribution
-
-    """
-    # Get copy of distribution list so as not to modify the argument
-    dist_list = deepcopy(dists)
-
-    # Get list of discrete/continuous distribution indices
-    discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
-    cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
-
-    # Apply probabilites to continuous distributions
-    for i in cont_index:
-        dist = dist_list[i]
-        dist.p *= probs[i]
-
-    if discrete_index:
-        # Create combined discrete distribution
-        dist_discrete = [dist_list[i] for i in discrete_index]
-        discrete_probs = [probs[i] for i in discrete_index]
-        combined_dist = Discrete.merge(dist_discrete, discrete_probs)
-
-        # Replace multiple discrete distributions with merged
-        for idx in reversed(discrete_index):
-            dist_list.pop(idx)
-        dist_list.append(combined_dist)
-
-    # Combine discrete and continuous if present
-    if len(dist_list) > 1:
-        probs = [d.integral() for d in dist_list]
-        dist_list[:] = [Mixture(probs, dist_list.copy())]
-
-    return dist_list[0]
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 733b7eaa6..f571377bf 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -1,6 +1,7 @@
 from abc import ABC, abstractmethod
 from collections import defaultdict
 from collections.abc import Iterable
+from copy import deepcopy
 from numbers import Real
 from xml.etree import ElementTree as ET
 
@@ -1238,3 +1239,51 @@ class Mixture(Univariate):
             p*dist.integral()
             for p, dist in zip(self.probability, self.distribution)
         ])
+
+
+def combine_distributions(dists, probs):
+    """Combine distributions with specified probabilities
+
+    This function can be used to combine multiple instances of
+    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular` into a single
+    distribution. Multiple discrete distributions are merged into a single
+    distribution and the remainder of the distributions are put into a
+    :class:`~openmc.stats.Mixture` distribution.
+
+    Parameters
+    ----------
+    dists : iterable of openmc.stats.Univariate
+        Distributions to combine
+    probs : iterable of float
+        Probability (or intensity) of each distribution
+
+    """
+    # Get copy of distribution list so as not to modify the argument
+    dist_list = deepcopy(dists)
+
+    # Get list of discrete/continuous distribution indices
+    discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)]
+    cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)]
+
+    # Apply probabilites to continuous distributions
+    for i in cont_index:
+        dist = dist_list[i]
+        dist.p *= probs[i]
+
+    if discrete_index:
+        # Create combined discrete distribution
+        dist_discrete = [dist_list[i] for i in discrete_index]
+        discrete_probs = [probs[i] for i in discrete_index]
+        combined_dist = Discrete.merge(dist_discrete, discrete_probs)
+
+        # Replace multiple discrete distributions with merged
+        for idx in reversed(discrete_index):
+            dist_list.pop(idx)
+        dist_list.append(combined_dist)
+
+    # Combine discrete and continuous if present
+    if len(dist_list) > 1:
+        probs = [d.integral() for d in dist_list]
+        dist_list[:] = [Mixture(probs, dist_list.copy())]
+
+    return dist_list[0]

From 5a9662aa9fbe0743a10a0f778910475dc9507a3a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 18 Jul 2022 10:54:58 -0500
Subject: [PATCH 08/12] Add test for combine_distributions, use np.array in
 Discrete/Tabular setters

---
 openmc/stats/univariate.py     | 22 +++++++++----------
 tests/unit_tests/test_stats.py | 39 ++++++++++++++++++++++++++++++----
 2 files changed, 46 insertions(+), 15 deletions(-)

diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index f571377bf..e11c9b50e 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -96,9 +96,9 @@ class Discrete(Univariate):
 
     Attributes
     ----------
-    x : Iterable of float
+    x : numpy.ndarray
         Values of the random variable
-    p : Iterable of float
+    p : numpy.ndarray
         Discrete probability for each value
 
     """
@@ -123,7 +123,7 @@ class Discrete(Univariate):
         if isinstance(x, Real):
             x = [x]
         cv.check_type('discrete values', x, Iterable, Real)
-        self._x = x
+        self._x = np.array(x, dtype=float)
 
     @p.setter
     def p(self, p):
@@ -132,7 +132,7 @@ class Discrete(Univariate):
         cv.check_type('discrete probabilities', p, Iterable, Real)
         for pk in p:
             cv.check_greater_than('discrete probability', pk, 0.0, True)
-        self._p = p
+        self._p = np.array(p, dtype=float)
 
     def cdf(self):
         return np.insert(np.cumsum(self.p), 0, 0.0)
@@ -836,9 +836,9 @@ class Tabular(Univariate):
 
     Attributes
     ----------
-    x : Iterable of float
+    x : numpy.ndarray
         Tabulated values of the random variable
-    p : Iterable of float
+    p : numpy.ndarray
         Tabulated probabilities
     interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional
         Indicate whether the density function is constant between tabulated
@@ -871,7 +871,7 @@ class Tabular(Univariate):
     @x.setter
     def x(self, x):
         cv.check_type('tabulated values', x, Iterable, Real)
-        self._x = x
+        self._x = np.array(x, dtype=float)
 
     @p.setter
     def p(self, p):
@@ -879,7 +879,7 @@ class Tabular(Univariate):
         if not self._ignore_negative:
             for pk in p:
                 cv.check_greater_than('tabulated probability', pk, 0.0, True)
-        self._p = p
+        self._p = np.array(p, dtype=float)
 
     @interpolation.setter
     def interpolation(self, interpolation):
@@ -892,8 +892,8 @@ class Tabular(Univariate):
                                       'distributions using histogram or '
                                       'linear-linear interpolation')
         c = np.zeros_like(self.x)
-        x = np.asarray(self.x)
-        p = np.asarray(self.p)
+        x = self.x
+        p = self.p
 
         if self.interpolation == 'histogram':
             c[1:] = p[:-1] * np.diff(x)
@@ -933,7 +933,7 @@ class Tabular(Univariate):
 
     def normalize(self):
         """Normalize the probabilities stored on the distribution"""
-        self.p = np.asarray(self.p) / self.cdf().max()
+        self.p /= self.cdf().max()
 
     def sample(self, n_samples=1, seed=None):
         if not self.interpolation in ('histogram', 'linear-linear'):
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index b57f9578a..404cecd13 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -13,8 +13,8 @@ def test_discrete():
     elem = d.to_xml_element('distribution')
 
     d = openmc.stats.Discrete.from_xml_element(elem)
-    assert d.x == x
-    assert d.p == p
+    np.testing.assert_array_equal(d.x, x)
+    np.testing.assert_array_equal(d.p, p)
     assert len(d) == len(x)
 
     d = openmc.stats.Univariate.from_xml_element(elem)
@@ -76,8 +76,8 @@ def test_uniform():
     assert len(d) == 2
 
     t = d.to_tabular()
-    assert t.x == [a, b]
-    assert t.p == [1/(b-a), 1/(b-a)]
+    np.testing.assert_array_equal(t.x, [a, b])
+    np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)])
     assert t.interpolation == 'histogram'
 
     exp_mean = 0.5 * (a + b)
@@ -396,3 +396,34 @@ def test_muir():
     assert within_2_sigma / n_samples >= 0.95
     within_3_sigma = np.count_nonzero(samples < 3*d.std_dev)
     assert within_3_sigma / n_samples >= 0.99
+
+
+def test_combine_distributions():
+    # Combine two discrete (same data as in test_merge_discrete)
+    x1 = [0.0, 1.0, 10.0]
+    p1 = [0.3, 0.2, 0.5]
+    d1 = openmc.stats.Discrete(x1, p1)
+    x2 = [0.5, 1.0, 5.0]
+    p2 = [0.4, 0.5, 0.1]
+    d2 = openmc.stats.Discrete(x2, p2)
+
+    # Merged distribution should have x values sorted and probabilities
+    # appropriately combined. Duplicate x values should appear once.
+    merged = openmc.stats.combine_distributions([d1, d2], [0.6, 0.4])
+    assert isinstance(merged, openmc.stats.Discrete)
+    assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0])
+    assert merged.p == pytest.approx(
+        [0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5])
+
+    # Probabilities add up but are not normalized
+    d1 = openmc.stats.Discrete([3.0], [1.0])
+    triple = openmc.stats.combine_distributions([d1, d1, d1], [1.0, 2.0, 3.0])
+    assert triple.x == pytest.approx([3.0])
+    assert triple.p == pytest.approx([6.0])
+
+    # Combine discrete and tabular
+    t1 = openmc.stats.Tabular(x2, p2)
+    mixed = openmc.stats.combine_distributions([d1, t1], [0.5, 0.5])
+    assert isinstance(mixed, openmc.stats.Mixture)
+    assert len(mixed.distribution) == 2
+    assert len(mixed.probability) == 2

From a68a3ede6a6fc85d0b8ac6eacd7b139fbbe6278e Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 18 Jul 2022 13:30:31 -0500
Subject: [PATCH 09/12] Change get_sources() -> sources property and add test

---
 openmc/data/decay.py                | 27 ++++++++++++----------
 openmc/stats/univariate.py          |  8 +++++++
 tests/unit_tests/test_data_decay.py | 35 +++++++++++++++++++++++++++++
 3 files changed, 58 insertions(+), 12 deletions(-)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index f358dadf7..96c3f562a 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -316,6 +316,12 @@ class Decay(EqualityMixin):
         'excited_state', 'mass', 'stable', 'spin', and 'parity'.
     spectra : dict
         Resulting radiation spectra for each radiation type.
+    sources : dict
+        Radioactive decay source distributions represented as a dictionary
+        mapping particle types (e.g., 'photon') to instances of
+        :class:`openmc.stats.Univariate`.
+
+        .. versionadded:: 0.13.1
 
     """
     def __init__(self, ev_or_filename):
@@ -498,16 +504,14 @@ class Decay(EqualityMixin):
         """
         return cls(ev_or_filename)
 
-    def get_sources(self):
-        """Get radioactive decay source distributions
+    @property
+    def sources(self):
+        """Radioactive decay source distributions"""
+        # If property has been computed already, return it
+        # TODO: Replace with functools.cached_property when support is Python 3.9+
+        if self._sources is not None:
+            return self._sources
 
-        Returns
-        -------
-        sources : dict
-            Dictionary mapping particle types (e.g., 'photon') to instances of
-            :class:`openmc.stats.Univariate`
-
-        """
         sources = {}
         name = self.nuclide['name']
         decay_constant = self.decay_constant.n
@@ -563,6 +567,5 @@ class Decay(EqualityMixin):
             merged_sources[particle_type] = combine_distributions(
                 dist_list, [1.0]*len(dist_list))
 
-        return merged_sources
-
-
+        self._sources = merged_sources
+        return self._sources
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index e11c9b50e..66c6e5f06 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -224,6 +224,8 @@ class Discrete(Univariate):
     def integral(self):
         """Return integral of distribution
 
+        .. versionadded:: 0.13.1
+
         Returns
         -------
         float
@@ -1041,6 +1043,8 @@ class Tabular(Univariate):
     def integral(self):
         """Return integral of distribution
 
+        .. versionadded: 0.13.1
+
         Returns
         -------
         float
@@ -1230,6 +1234,8 @@ class Mixture(Univariate):
     def integral(self):
         """Return integral of the distribution
 
+        .. versionadded:: 0.13.1
+
         Returns
         -------
         float
@@ -1250,6 +1256,8 @@ def combine_distributions(dists, probs):
     distribution and the remainder of the distributions are put into a
     :class:`~openmc.stats.Mixture` distribution.
 
+    .. versionadded:: 0.13.1
+
     Parameters
     ----------
     dists : iterable of openmc.stats.Univariate
diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py
index 06b8e6bed..f0bda1bd9 100644
--- a/tests/unit_tests/test_data_decay.py
+++ b/tests/unit_tests/test_data_decay.py
@@ -23,6 +23,14 @@ def nb90():
     return openmc.data.Decay.from_endf(filename)
 
 
+@pytest.fixture(scope='module')
+def ba137m():
+    """Ba137_m1 decay data."""
+    endf_data = os.environ['OPENMC_ENDF_DATA']
+    filename = os.path.join(endf_data, 'decay', 'dec-056_Ba_137m1.endf')
+    return openmc.data.Decay.from_endf(filename)
+
+
 @pytest.fixture(scope='module')
 def u235_yields():
     """U235 fission product yield data."""
@@ -48,6 +56,7 @@ def test_nb90_halflife(nb90):
     ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life)
     ufloat_close(nb90.decay_energy, ufloat(2265527.5, 25159.400474401213))
 
+
 def test_nb90_nuclide(nb90):
     assert nb90.nuclide['atomic_number'] == 41
     assert nb90.nuclide['mass_number'] == 90
@@ -91,3 +100,29 @@ def test_fpy(u235_yields):
     assert len(u235_yields.independent) == 3
     thermal = u235_yields.independent[0]
     ufloat_close(thermal['I135'], ufloat(0.0292737, 0.000819663))
+
+
+def test_sources(ba137m, nb90):
+    # Running .sources twice should give same objects
+    sources = ba137m.sources
+    sources2 = ba137m.sources
+    for key in sources:
+        assert sources[key] is sources2[key]
+
+    # Each source should be a univariate distribution
+    for dist in sources.values():
+        assert isinstance(dist, openmc.stats.Univariate)
+
+    # Check for presence of 662 keV gamma ray in decay of Ba137m
+    gamma_source = ba137m.sources['photon']
+    assert isinstance(gamma_source, openmc.stats.Discrete)
+    b = np.isclose(gamma_source.x, 661657.)
+    assert np.count_nonzero(b) == 1
+
+    # Check value of decay/s/atom
+    idx = np.flatnonzero(b)[0]
+    assert gamma_source.p[idx] == pytest.approx(0.004069614)
+
+    # Nb90 decays by β+ and should emit positrons, electrons, and photons
+    sources = nb90.sources
+    assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0

From 3e5afd21eb67512c05111ad49e848b52f85dc981 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 22 Jul 2022 16:42:24 -0500
Subject: [PATCH 10/12] Add missing definition of _sources in Decay.__init__

---
 openmc/data/decay.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 96c3f562a..2c2505bf5 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -337,6 +337,7 @@ class Decay(EqualityMixin):
         self.modes = []
         self.spectra = {}
         self.average_energies = {}
+        self._sources = None
 
         # Get head record
         items = get_head_record(file_obj)

From b2fed3c3cec856b85f7fb75aa8f4a20e1e974a84 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 28 Jul 2022 17:06:13 -0500
Subject: [PATCH 11/12] Fix bugs in Discrete, Tabular, and Mixture sample
 methods

---
 openmc/stats/univariate.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 66c6e5f06..c19f45365 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -139,7 +139,8 @@ class Discrete(Univariate):
 
     def sample(self, n_samples=1, seed=None):
         np.random.seed(seed)
-        return np.random.choice(self.x, n_samples, p=self.p)
+        p = self.p / self.p.sum()
+        return np.random.choice(self.x, n_samples, p=p)
 
     def normalize(self):
         """Normalize the probabilities stored on the distribution"""
@@ -944,10 +945,11 @@ class Tabular(Univariate):
                                       'linear-linear interpolation')
         np.random.seed(seed)
         xi = np.random.rand(n_samples)
-        cdf = self.cdf()
-        cdf /= cdf.max()
+
         # always use normalized probabilities when sampling
+        cdf = self.cdf()
         p = self.p / cdf.max()
+        cdf /= cdf.max()
 
         # get CDF bins that are above the
         # sampled values
@@ -962,7 +964,7 @@ class Tabular(Univariate):
         # the random number is less than the next cdf
         # entry
         x_i = self.x[cdf_idx]
-        p_i = self.p[cdf_idx]
+        p_i = p[cdf_idx]
 
         if self.interpolation == 'histogram':
             # mask where probability is greater than zero
@@ -980,7 +982,7 @@ class Tabular(Univariate):
             # get variable and probability values for the
             # next entry
             x_i1 = self.x[cdf_idx + 1]
-            p_i1 = self.p[cdf_idx + 1]
+            p_i1 = p[cdf_idx + 1]
             # compute slope between entries
             m = (p_i1 - p_i) / (x_i1 - x_i)
             # set values for zero slope
@@ -1166,9 +1168,10 @@ class Mixture(Univariate):
 
     def sample(self, n_samples=1, seed=None):
         np.random.seed(seed)
-        idx = np.random.choice(self.distribution, n_samples, p=self.probability)
+        idx = np.random.choice(range(len(self.distribution)),
+                               n_samples, p=self.probability)
 
-        out = np.zeros_like(idx)
+        out = np.empty_like(idx, dtype=float)
         for i in np.unique(idx):
             n_dist_samples = np.count_nonzero(idx == i)
             samples = self.distribution[i].sample(n_dist_samples)

From cc7aa092be7b7a9eed9c9ddd85a9e1edb0a1f14b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 28 Jul 2022 17:07:06 -0500
Subject: [PATCH 12/12] Address @eepeterson comments on #2135

---
 openmc/stats/univariate.py     |  8 ++--
 tests/unit_tests/test_stats.py | 69 ++++++++++++++++++++--------------
 2 files changed, 44 insertions(+), 33 deletions(-)

diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index c19f45365..58b7c6172 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -1254,10 +1254,10 @@ def combine_distributions(dists, probs):
     """Combine distributions with specified probabilities
 
     This function can be used to combine multiple instances of
-    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular` into a single
-    distribution. Multiple discrete distributions are merged into a single
-    distribution and the remainder of the distributions are put into a
-    :class:`~openmc.stats.Mixture` distribution.
+    :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular`. Multiple
+    discrete distributions are merged into a single distribution and the
+    remainder of the distributions are put into a :class:`~openmc.stats.Mixture`
+    distribution.
 
     .. versionadded:: 0.13.1
 
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index 404cecd13..1cfcf00c3 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -6,6 +6,12 @@ import openmc
 import openmc.stats
 
 
+def assert_sample_mean(samples, expected_mean):
+    std_dev = samples.std() / np.sqrt(samples.size)
+    assert np.abs(expected_mean - samples.mean()) < 3*std_dev
+
+
+
 def test_discrete():
     x = [0.0, 1.0, 10.0]
     p = [0.3, 0.2, 0.5]
@@ -33,13 +39,11 @@ def test_discrete():
 
     d3 = openmc.stats.Discrete(vals, probs)
 
-    # sample discrete distribution
+    # sample discrete distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d3.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_merge_discrete():
@@ -80,13 +84,12 @@ def test_uniform():
     np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)])
     assert t.interpolation == 'histogram'
 
+    # Sample distribution and check that the mean of the samples is within 3
+    # std. dev. of the expected mean
     exp_mean = 0.5 * (a + b)
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_powerlaw():
@@ -102,13 +105,11 @@ def test_powerlaw():
 
     exp_mean = 100.0 * (n+1) / (n+2)
 
-    # sample power law distribution
+    # sample power law distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_maxwell():
@@ -122,21 +123,15 @@ def test_maxwell():
 
     exp_mean = 3/2 * theta
 
-    # sample maxwell distribution
+    # sample maxwell distribution and check that the mean of the samples is
+    # within 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
     # A second sample with a different seed
     samples_2 = d.sample(n_samples, seed=200)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples_2.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples_2.mean()) < 3*std_dev
-
+    assert_sample_mean(samples_2, exp_mean)
     assert samples_2.mean() != samples.mean()
 
 
@@ -156,13 +151,11 @@ def test_watt():
     # https://doi.org/10.1016/j.physletb.2003.09.048
     exp_mean = 3/2 * a + a**2 * b / 4
 
-    # sample Watt distribution
+    # sample Watt distribution and check that the mean of the samples is within
+    # 3 std. dev. of the expected mean
     n_samples = 1_000_000
     samples = d.sample(n_samples, seed=100)
-    # check that the mean of the samples is within 3 std. dev.
-    # of the expected mean
-    std_dev = samples.std() / np.sqrt(n_samples)
-    assert np.abs(exp_mean - samples.mean()) < 3*std_dev
+    assert_sample_mean(samples, exp_mean)
 
 
 def test_tabular():
@@ -228,6 +221,11 @@ def test_mixture():
     assert mix.distribution == [d1, d2]
     assert len(mix) == 4
 
+    # Sample and make sure sample mean is close to expected mean
+    n_samples = 1_000_000
+    samples = mix.sample(n_samples)
+    assert_sample_mean(samples, (2.5 + 5.0)/2)
+
     elem = mix.to_xml_element('distribution')
 
     d = openmc.stats.Mixture.from_xml_element(elem)
@@ -427,3 +425,16 @@ def test_combine_distributions():
     assert isinstance(mixed, openmc.stats.Mixture)
     assert len(mixed.distribution) == 2
     assert len(mixed.probability) == 2
+
+    # Combine 1 discrete and 2 tabular -- the tabular distributions should
+    # combine to produce a uniform distribution with mean 0.5. The combined
+    # distribution should have a mean of 0.25.
+    t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0])
+    t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0])
+    d1 = openmc.stats.Discrete([0.0], [1.0])
+    combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5])
+
+    # Sample the combined distribution and make sure the sample mean is within
+    # uncertainty of the expected value
+    samples = combined.sample(10)
+    assert_sample_mean(samples, 0.25)