mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge pull request #2135 from paulromano/decay-data
Add ability to get sources from `Decay` class
This commit is contained in:
Ethan Peterson
GitHub
commit
3f3b0282f7
5 changed files with 303 additions and 51 deletions
|
|
@ -61,6 +61,7 @@ Core Functions
|
|||
|
||||
atomic_mass
|
||||
atomic_weight
|
||||
combine_distributions
|
||||
decay_constant
|
||||
dose_coefficients
|
||||
gnd_name
|
||||
|
|
|
|||
|
|
@ -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, 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
|
||||
|
||||
|
||||
|
|
@ -314,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):
|
||||
|
|
@ -329,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)
|
||||
|
|
@ -495,3 +504,69 @@ class Decay(EqualityMixin):
|
|||
|
||||
"""
|
||||
return cls(ev_or_filename)
|
||||
|
||||
@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
|
||||
|
||||
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 = {
|
||||
'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] = []
|
||||
|
||||
# Create distribution for discrete
|
||||
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)
|
||||
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
|
||||
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")
|
||||
|
||||
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
|
||||
merged_sources = {}
|
||||
for particle_type, dist_list in sources.items():
|
||||
merged_sources[particle_type] = combine_distributions(
|
||||
dist_list, [1.0]*len(dist_list))
|
||||
|
||||
self._sources = merged_sources
|
||||
return self._sources
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -95,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
|
||||
|
||||
"""
|
||||
|
|
@ -122,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):
|
||||
|
|
@ -131,14 +132,15 @@ 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)
|
||||
|
||||
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"""
|
||||
|
|
@ -220,6 +222,18 @@ 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
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Integral of discrete distribution
|
||||
"""
|
||||
return np.sum(self.p)
|
||||
|
||||
class Uniform(Univariate):
|
||||
"""Distribution with constant probability over a finite interval [a,b]
|
||||
|
||||
|
|
@ -825,9 +839,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
|
||||
|
|
@ -860,7 +874,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):
|
||||
|
|
@ -868,7 +882,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):
|
||||
|
|
@ -881,8 +895,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)
|
||||
|
|
@ -922,7 +936,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'):
|
||||
|
|
@ -931,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
|
||||
|
|
@ -949,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
|
||||
|
|
@ -967,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
|
||||
|
|
@ -1027,6 +1042,24 @@ class Tabular(Univariate):
|
|||
p = params[len(params)//2:]
|
||||
return cls(x, p, interpolation)
|
||||
|
||||
def integral(self):
|
||||
"""Return integral of distribution
|
||||
|
||||
.. versionadded: 0.13.1
|
||||
|
||||
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
|
||||
|
|
@ -1135,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)
|
||||
|
|
@ -1199,3 +1233,68 @@ class Mixture(Univariate):
|
|||
distribution.append(Univariate.from_xml_element(pair.find("dist")))
|
||||
|
||||
return cls(probability, distribution)
|
||||
|
||||
def integral(self):
|
||||
"""Return integral of the distribution
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Integral of the distribution
|
||||
"""
|
||||
return sum([
|
||||
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`. 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
|
||||
|
||||
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]
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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]
|
||||
|
|
@ -13,8 +19,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)
|
||||
|
|
@ -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():
|
||||
|
|
@ -76,17 +80,16 @@ 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'
|
||||
|
||||
# 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)
|
||||
|
|
@ -396,3 +394,47 @@ 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
|
||||
|
||||
# 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)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue