This notebook attemps to benchmark the symbolic decayer to ORIGEN v2.2 results. The basic strategy is that at $t=0$, a single kg of a pure nuclude is decayed for 0.01, 0.1, 1.0, 10, 100 times that nuclide's half life. If the nuclide is stable then the decay calculation is run for 0.01, 0.1, 1.0, 10, 100 seconds.
Throughout this benchmark I use the mean relative error as the basis for comparison. For any points $x$ and $y$:
$e = \frac{2|x - y|}{x + y} $
This may also be given in a mass-weighted form, $w$, for mass $m$.
$w = m\cdot e$
The mass weighted form is useful to discrimnate against species that have a high relative error but are virtually absent in the decayed material. Let's start with some imports!
import os
import sys
import shutil
import warnings
import subprocess
from functools import lru_cache
import numpy as np
import matplotlib
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
from pyne.utils import QAWarning
warnings.simplefilter('ignore', QAWarning)
from pyne import nucname
from pyne import data
from pyne import origen22
from pyne.material import Material, MaterialLibrary
import sys
sys.path.insert(0, '/home/scopatz/pyne/src')
import decaygen
exp2 = np.exp2
H5NAME = 'origen-benchmark-py{0}k.h5'.format(sys.version_info[0])
mats = MaterialLibrary(H5NAME)
Testing the decayer is meaningless if the half-lives that are used do not match each other. This adds a function to compare the relative error of the pyne half life to the origen one.
t9 = origen22.parse_tape9('TAPE9.INP')
for key in list(t9.keys()):
if key not in [1, 2, 3]:
del t9[key]
def t9_half_life(nuc):
nuczz = nucname.zzaaam(nuc)
hl = None
for i in [1, 2, 3]:
hl = t9[i]['half_life'].get(nuczz, None)
if hl is not None:
break
return hl
@lru_cache()
def hl_relerr(nuc):
"""Half-life relative error."""
dhl = data.half_life(nuc) or np.inf
ohl = t9_half_life(nuc) or np.inf
if np.isinf(ohl) and np.isinf(dhl):
return 0.0
hlre = np.abs(ohl - dhl) * 2 / (ohl + dhl)
return np.inf if np.isnan(hlre) else hlre
import json
t9hls = {}
for i in [1, 2, 3]:
t9hls.update({nucname.id(int(nuc)): v for nuc, v in t9[i]['half_life'].items()})
with open('o2hls.json', 'w') as f:
json.dump(t9hls, f, sort_keys=True)
This puts the material keys into a dict mapping nuclides to the data set names (eg 'origen_922350000_3'). This filters out any data sets that contain high relative errors in the half-lives. It also filters out species for which Origen has weird metastable behavior that seems unphysical.
metastable_blacklist = {
771940001, # have 2+ metastables that ORIGEN lumps together, or
340830001, # ...
350830000, # decay to metastables without reason.
340830000, # ...
481170000, # ...
461110000, # ...
441030000, # ...
471110001, # missing branch in origen
501130001, # ...
}
nuc_keys = {}
for key in mats.keys():
_, nuc, t = key.split('_')
nuc = int(nuc)
if nuc in metastable_blacklist:
continue
chains = decaygen.genchains([(nuc,)])
maxrelerr = max([max(list(map(hl_relerr, c))) for c in chains])
if maxrelerr > 0.1:
continue
t = int(t)
if nuc not in nuc_keys:
nuc_keys[nuc] = []
nuc_keys[nuc].append(key)
sfunc = lambda x: int(x.rsplit('_', 1)[1])
for keys in nuc_keys.values():
keys.sort(key=sfunc)
Unitly functions to convert materials to pandas Series and compute the relative error of each nuclide bewteen two materials.
def mat_to_series(x):
xmass = x.mult_by_mass()
s = pd.Series(list(x.values()), list(map(int, x.keys())))
return s
def matdiff(x, y):
if not isinstance(x, pd.Series):
x = mat_to_series(x)
if not isinstance(y, pd.Series):
y = mat_to_series(y)
return (x - y).abs() * (2 / (x + y))
This constructs a new data frame that compares the origen decayed material to the pyne material decay. The comparison stores:
Note that this also filters by when Origen itself has messed up by gaining or lossing too much mass (1%). Furthermore, the maximum relative error is only computed for species that have a certain threshold mass (1e-3).
idx = []
err = []
for nuc, keys in list(nuc_keys.items()):
fresh = mats[keys[0]]
for matkey in keys[1:]:
mat = mats[matkey]
if (mat.mass < 0.99) or (mat.mass > 1.1):
# Origen lost too much mass
continue
t = mat.metadata['decay_time']
dmat = fresh.decay(t)
oser = mat_to_series(mat)
dser = mat_to_series(dmat)
md = matdiff(oser, dser)
md = md.fillna(value=0.0)
mask = (dser > 1e-3) | (oser > 1e-3)
mdmasked = md.where(mask)
maxmd = mdmasked.max()
child = mdmasked.argmax()
if np.isnan(child):
continue
childmass = max(dser[child], oser[child])
weightederr = maxmd * childmass
idx.append(matkey)
err.append((nuc, t, maxmd, child, childmass, weightederr))
df = pd.DataFrame(err, index=idx, columns=['nuc', 't', 'relerr', 'child', 'childmass', 'weightederr'])
A utility function for printing out nuclide summaries
def nucsummary(nuc, i=None):
print('Nulcide:', nucname.name(nuc))
print("origen half life:", t9_half_life(nuc))
print("data.half_life():", data.half_life(nuc))
print('relerr half life:', hl_relerr(nuc))
print()
subdf = df[df['nuc'] == nuc]
plt.figure(figsize=(10, 6))
subdf.plot(x='t', y='relerr', logx=True, label=nucname.name(int(nuc)))
plt.legend(loc=0)
plt.ylabel('t [sec]')
plt.ylabel('Relative Error')
if i is not None:
okey = 'origen_{0}_{1}'.format(nuc, i)
fkey = 'fresh_{0}_0'.format(nuc)
omat = mats[okey]
fresh = mats[fkey]
t = omat.metadata['decay_time']
dmat = fresh.decay(t)
md = matdiff(omat, dmat)
print(omat)
print()
print(dmat)
print()
print(md)
print(subdf.to_latex())
return subdf
The data frame, sorted by relative error, childmass, and nuclide.
df.sort(['relerr', 'childmass', 'nuc'], ascending=[False, False, True])
nuc | t | relerr | child | childmass | weightederr | |
---|---|---|---|---|---|---|
origen_420930001_4 | 420930001 | 2.466000e+05 | 1.994709 | 410930000 | 0.001168 | 0.002330 |
origen_420930001_5 | 420930001 | 2.466000e+06 | 1.948262 | 410930000 | 0.001181 | 0.002302 |
origen_781930001_5 | 781930001 | 3.741120e+07 | 1.629114 | 771930000 | 0.016067 | 0.026175 |
origen_781930001_4 | 781930001 | 3.741120e+06 | 1.581632 | 771930000 | 0.001406 | 0.002223 |
origen_420990000_1 | 420990000 | 2.375136e+03 | 1.557261 | 430990000 | 0.006907 | 0.010757 |
origen_491140000_4 | 491140000 | 7.190000e+02 | 1.483783 | 481140000 | 0.033710 | 0.050018 |
origen_491140000_5 | 491140000 | 7.190000e+03 | 1.483744 | 481140000 | 0.033740 | 0.050062 |
origen_491140000_3 | 491140000 | 7.190000e+01 | 1.483744 | 481140000 | 0.016870 | 0.025031 |
origen_491140000_2 | 491140000 | 7.190000e+00 | 1.483651 | 481140000 | 0.002259 | 0.003352 |
origen_380910000_1 | 380910000 | 3.474000e+02 | 0.803749 | 390910000 | 0.006907 | 0.005552 |
origen_491140001_4 | 491140001 | 4.277664e+07 | 0.674001 | 481140000 | 0.075220 | 0.050698 |
origen_491140001_5 | 491140001 | 4.277664e+08 | 0.673959 | 481140000 | 0.075290 | 0.050742 |
origen_491140001_3 | 491140001 | 4.277664e+06 | 0.673843 | 481140000 | 0.037640 | 0.025363 |
origen_491140001_2 | 491140001 | 4.277664e+05 | 0.673820 | 481140000 | 0.005041 | 0.003397 |
origen_230540000_4 | 230540000 | 4.980000e+02 | 0.633013 | 230540000 | 0.001881 | 0.001191 |
origen_360790001_4 | 360790001 | 5.000000e+02 | 0.610648 | 360790001 | 0.001835 | 0.001121 |
origen_561330001_4 | 561330001 | 1.401480e+06 | 0.587649 | 551330000 | 0.002970 | 0.001745 |
origen_842110000_4 | 842110000 | 5.160000e+00 | 0.531242 | 842110000 | 0.001683 | 0.000894 |
origen_952420001_4 | 952420001 | 4.449622e+10 | 0.519896 | 942380000 | 0.001293 | 0.000672 |
origen_80190000_4 | 80190000 | 2.688000e+02 | 0.496765 | 80190000 | 0.001622 | 0.000806 |
origen_741890000_4 | 741890000 | 6.420000e+03 | 0.472665 | 741890000 | 0.001581 | 0.000747 |
origen_320710001_4 | 320710001 | 2.041000e-01 | 0.462447 | 320710001 | 0.001564 | 0.000723 |
origen_802050000_4 | 802050000 | 3.084000e+03 | 0.448400 | 802050000 | 0.001541 | 0.000691 |
origen_461110001_1 | 461110001 | 1.980000e+02 | 0.402345 | 471110000 | 0.002104 | 0.000846 |
origen_40100000_4 | 40100000 | 4.765198e+14 | 0.387555 | 40100000 | 0.001446 | 0.000560 |
origen_781910000_4 | 781910000 | 2.445120e+06 | 0.387555 | 781910000 | 0.001446 | 0.000560 |
origen_561310000_4 | 561310000 | 9.936000e+06 | 0.346184 | 551310000 | 0.005383 | 0.001864 |
origen_280590000_4 | 280590000 | 2.398378e+13 | 0.344502 | 280590000 | 0.001383 | 0.000476 |
origen_461090001_3 | 461090001 | 2.817600e+02 | 0.341814 | 471090000 | 0.001102 | 0.000377 |
origen_812040000_2 | 812040000 | 1.193824e+07 | 0.331065 | 802040000 | 0.001955 | 0.000647 |
... | ... | ... | ... | ... | ... | ... |
origen_802040000_2 | 802040000 | 1.000000e-01 | 0.000000 | 802040000 | 1.000000 | 0.000000 |
origen_802040000_3 | 802040000 | 1.000000e+00 | 0.000000 | 802040000 | 1.000000 | 0.000000 |
origen_802040000_4 | 802040000 | 1.000000e+01 | 0.000000 | 802040000 | 1.000000 | 0.000000 |
origen_802040000_5 | 802040000 | 1.000000e+02 | 0.000000 | 802040000 | 1.000000 | 0.000000 |
origen_812030000_1 | 812030000 | 1.000000e-02 | 0.000000 | 812030000 | 1.000000 | 0.000000 |
origen_812030000_2 | 812030000 | 1.000000e-01 | 0.000000 | 812030000 | 1.000000 | 0.000000 |
origen_812030000_3 | 812030000 | 1.000000e+00 | 0.000000 | 812030000 | 1.000000 | 0.000000 |
origen_812030000_4 | 812030000 | 1.000000e+01 | 0.000000 | 812030000 | 1.000000 | 0.000000 |
origen_812030000_5 | 812030000 | 1.000000e+02 | 0.000000 | 812030000 | 1.000000 | 0.000000 |
origen_812050000_1 | 812050000 | 1.000000e-02 | 0.000000 | 812050000 | 1.000000 | 0.000000 |
origen_812050000_2 | 812050000 | 1.000000e-01 | 0.000000 | 812050000 | 1.000000 | 0.000000 |
origen_812050000_3 | 812050000 | 1.000000e+00 | 0.000000 | 812050000 | 1.000000 | 0.000000 |
origen_812050000_4 | 812050000 | 1.000000e+01 | 0.000000 | 812050000 | 1.000000 | 0.000000 |
origen_812050000_5 | 812050000 | 1.000000e+02 | 0.000000 | 812050000 | 1.000000 | 0.000000 |
origen_822060000_1 | 822060000 | 1.000000e-02 | 0.000000 | 822060000 | 1.000000 | 0.000000 |
origen_822060000_2 | 822060000 | 1.000000e-01 | 0.000000 | 822060000 | 1.000000 | 0.000000 |
origen_822060000_3 | 822060000 | 1.000000e+00 | 0.000000 | 822060000 | 1.000000 | 0.000000 |
origen_822060000_4 | 822060000 | 1.000000e+01 | 0.000000 | 822060000 | 1.000000 | 0.000000 |
origen_822060000_5 | 822060000 | 1.000000e+02 | 0.000000 | 822060000 | 1.000000 | 0.000000 |
origen_822070000_1 | 822070000 | 1.000000e-02 | 0.000000 | 822070000 | 1.000000 | 0.000000 |
origen_822070000_2 | 822070000 | 1.000000e-01 | 0.000000 | 822070000 | 1.000000 | 0.000000 |
origen_822070000_3 | 822070000 | 1.000000e+00 | 0.000000 | 822070000 | 1.000000 | 0.000000 |
origen_822070000_4 | 822070000 | 1.000000e+01 | 0.000000 | 822070000 | 1.000000 | 0.000000 |
origen_822070000_5 | 822070000 | 1.000000e+02 | 0.000000 | 822070000 | 1.000000 | 0.000000 |
origen_822080000_1 | 822080000 | 1.000000e-02 | 0.000000 | 822080000 | 1.000000 | 0.000000 |
origen_822080000_2 | 822080000 | 1.000000e-01 | 0.000000 | 822080000 | 1.000000 | 0.000000 |
origen_822080000_3 | 822080000 | 1.000000e+00 | 0.000000 | 822080000 | 1.000000 | 0.000000 |
origen_822080000_4 | 822080000 | 1.000000e+01 | 0.000000 | 822080000 | 1.000000 | 0.000000 |
origen_822080000_5 | 822080000 | 1.000000e+02 | 0.000000 | 822080000 | 1.000000 | 0.000000 |
origen_721780001_3 | 721780001 | 4.000000e+00 | 0.000000 | 721780001 | 0.500000 | 0.000000 |
2802 rows × 6 columns
The fraction of the data with a relative error greater than 50%
(df.relerr > 0.5).sum() / len(df)
0.0067808708065667384
The fraction of the data with a mass-weighted relative error greater than 50%
(df.weightederr > 0.5).sum() / len(df)
0.0
The fraction of the data with both relative and mass-weighted relative error greater than 50%
((df.relerr > 0.5) & (df.weightederr > 0.5)).sum() / len(df)
0.0
The number of decay calculations compared
len(df)
2802
The relative errors, grouped by nuclide
N = 20
plt.figure(figsize=(10, 6))
for i, (key, grp) in enumerate(df.groupby(['nuc'])):
if i > N:
break
nn = nucname.name(int(key))
p = grp.plot(x='t', y='relerr', logx=True, label=nn)
plt.xlabel('t [sec]')
plt.ylabel('Relative Error')
p.legend(p.lines, (nn,), loc=0)
plt.savefig('decay-relative-error-{0}.eps'.format(nn))
<matplotlib.figure.Figure at 0x7f34d7a2a208>
Should explain the issue with remaining high relative error
Origen seems to have the wrong branch ratios
nucsummary(420990000, 1)
Nulcide: Mo99 origen half life: 237600.0 data.half_life(): 237513.6 relerr half life: 0.000363702491362 Material: mass = 0.9939597000001054 density = -1.0 atoms per molecule = -1.0 decay_time = 2375.136 name = origen_420990000_1 ------------------------- Mo99 0.9991350755970234 Tc99 0.000864924402870568 Ru99 1.0604051653199704e-13 Material: mass = 0.9999999336058844 density = -1.0 atoms per molecule = 1.0000000354067466 --------------------------------------- Mo99 0.9930925617347804 Tc99 0.006907402611982175 Ru99 3.565323735891612e-08 420990000 0.000008 430990000 1.557261 440990000 1.999988 dtype: float64
nuc | t | relerr | child | childmass | weightederr | |
---|---|---|---|---|---|---|
origen_420990000_1 | 420990000 | 2375.136 | 1.557261 | 430990000 | 0.006907 | 0.010757 |
<matplotlib.figure.Figure at 0x7f34e1ff8780>
Origen just seems to have the wrong branch ratios
nucsummary(491140000, 3)
Nulcide: In114 origen half life: 71.9 data.half_life(): 71.9 relerr half life: 0.0 Material: mass = 0.9999699999999999 density = -1.0 atoms per molecule = -1.0 decay_time = 71.9 name = origen_491140000_3 ------------------------- Cd114 0.016870506115183457 In114 0.5000150004500135 Sn114 0.48311449343480306 Material: mass = 0.999990639989361 density = -1.0 atoms per molecule = 1.0 ------------------------ Cd114 0.002499989836490059 In114 0.500004576422662 Sn114 0.4974954337408478 481140000 1.483744 491140000 0.000000 501140000 0.029351 dtype: float64 \begin{tabular}{lrrrrrr} \toprule {} & nuc & t & relerr & child & childmass & weightederr \\ \midrule origen\_491140000\_1 & 491140000 & 0.719 & 0.029356 & 501140000 & 0.006873 & 0.000202 \\ origen\_491140000\_2 & 491140000 & 7.190 & 1.483651 & 481140000 & 0.002259 & 0.003352 \\ origen\_491140000\_3 & 491140000 & 71.900 & 1.483744 & 481140000 & 0.016870 & 0.025031 \\ origen\_491140000\_4 & 491140000 & 719.000 & 1.483783 & 481140000 & 0.033710 & 0.050018 \\ origen\_491140000\_5 & 491140000 & 7190.000 & 1.483744 & 481140000 & 0.033740 & 0.050062 \\ \bottomrule \end{tabular}
nuc | t | relerr | child | childmass | weightederr | |
---|---|---|---|---|---|---|
origen_491140000_1 | 491140000 | 0.719 | 0.029356 | 501140000 | 0.006873 | 0.000202 |
origen_491140000_2 | 491140000 | 7.190 | 1.483651 | 481140000 | 0.002259 | 0.003352 |
origen_491140000_3 | 491140000 | 71.900 | 1.483744 | 481140000 | 0.016870 | 0.025031 |
origen_491140000_4 | 491140000 | 719.000 | 1.483783 | 481140000 | 0.033710 | 0.050018 |
origen_491140000_5 | 491140000 | 7190.000 | 1.483744 | 481140000 | 0.033740 | 0.050062 |
<matplotlib.figure.Figure at 0x7f34f08eaf60>
Relative error is high for Y-91, but the mass weighted relative error is low
nucsummary(380910000, 1)
Nulcide: Sr91 origen half life: 34200.0 data.half_life(): 34740.0 relerr half life: 0.0156657963446 Material: mass = 0.9959470702700001 density = -1.0 atoms per molecule = -1.0 decay_time = 347.4 name = origen_380910000_1 ------------------------- Sr91 0.9970409368550067 Y91 0.0029589925890349494 Zr91 7.055595834118967e-08 Material: mass = 0.9999996178756235 density = -1.0 atoms per molecule = 0.9999998380652019 --------------------------------------- Sr91 0.9930928753873625 Y91 0.006907124612637552 380910000 0.000093 390910000 0.803749 400910000 NaN dtype: float64
nuc | t | relerr | child | childmass | weightederr | |
---|---|---|---|---|---|---|
origen_380910000_1 | 380910000 | 347.4 | 0.803749 | 390910000 | 0.006907 | 0.005552 |
<matplotlib.figure.Figure at 0x7f34e1fd5e10>
... and so on
df.relerr.plot(kind='hist', bins=10000, loglog=True, figsize=(10, 6))
plt.xlabel('Relative Error, on range [0, 2]')
plt.ylabel('Number of Calculations, Total = {0}'.format(len(df)))
plt.savefig('decay-relative-error-histogram.eps')