Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Endf get cov #189

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 112 additions & 0 deletions notebooks/notebook_Xs_cov_from_enf6.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sandy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import seaborn as sns\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"endf6 = sandy.get_endf6_file('jeff_33','xs', 922380)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plot_corr(cov_obj, ax, **kwargs):\n",
" cov_obj_ = sandy.CategoryCov(cov_obj) if not isinstance(cov_obj, sandy.CategoryCov) else cov_obj\n",
" add = {\"cbar\": True, \"vmin\": -1, \"vmax\": 1, \"cmap\": \"RdBu\"}\n",
" for k, v in kwargs.items():\n",
" add[k] = v\n",
" ax = sns.heatmap(cov_obj_.get_corr().data, ax=ax, **add)\n",
" return ax"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot cross section correlation matrix from Endf6 files (MF=33)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cov = endf6._get_xs_cov()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(12, 12), dpi=100)\n",
"plot_corr(cov, ax=ax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot only capture (mt=102) reaction"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(6, 6), dpi=100)\n",
"plot_corr(cov.data.loc[(9237, 102), (9237, 102)], ax=ax)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python3 (sandy-devel)",
"language": "python",
"name": "sandy-devel"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.10"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
159 changes: 159 additions & 0 deletions sandy/core/endf6.py
Original file line number Diff line number Diff line change
Expand Up @@ -2475,3 +2475,162 @@ def covariance_info(self, nubar=True, xs=True, mubar=True, chi=True):
'chi': run_chi if chi else False,
}
return cov_info

def _get_xs_cov(self, mf=[31, 33], mt=range(1, 10000),
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
mat=range(1, 10000)):
"""
Extract from endf6 file the cross section/nubar covariance matrices.

Parameters
----------
mf : `int` or `list`, optional
MF number. The default is [31, 33].
mt : `int` or `list`, optional
MT number. The default is range(1, 10000).
mat : `int` or `list`, optional
MAT number. The default is range(1, 10000).

Returns
-------
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
`sandy.CategoryCov` or None if no covariance information is found.
Covariance matrix.

Examples
--------
>>> endf6 = sandy.get_endf6_file('jeff_33','xs', 922380)
>>> cov = endf6._get_xs_cov()
>>> xs_cov = sandy.XsCov.from_endf6(endf6)
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
>>> assert (cov.data == xs_cov).all().all()
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
"""
listmf_ = [mf] if isinstance(mf, int) else mf
listmt_ = [mt] if isinstance(mt, int) else mt
listmat_ = [mat] if isinstance(mat, int) else mat
tape = self.filter_by(listmf=listmf_, listmt=listmt_, listmat=listmat_)
data = []
for mat, mf, mt in tape.data:
sec = tape.read_section(mat, mf, mt)
for sub in sec["SUB"].values():
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
mat1 = sub['MAT1'] if sub['MAT1'] != 0 else mat
mt1 = sub['MT1']
covs = []
# Loop NI-type covariances
for nisec in sub["NI"].values():
lb = nisec["LB"]
if lb == 5:
if nisec["LS"] == 0:
foo = sandy.EnergyCov.from_lb5_asym
else:
foo = sandy.EnergyCov.from_lb5_sym
cov = foo(nisec["EK"], nisec["FKK"])
elif lb == 1:
foo = sandy.EnergyCov.from_lb1
cov = foo(nisec["EK"], nisec["FK"])
elif lb == 2:
foo = sandy.EnergyCov.from_lb2
cov = foo(nisec["EK"], nisec["FK"])
elif lb == 6:
foo = sandy.EnergyCov.from_lb6
cov = foo(nisec["EK"], nisec["EL"], nisec["FKL"])
else:
logging.warning(f"skip 'LB={lb}' covariance for"
f" [({mat}/{mt}), ({mat1}/{mt1})]")
continue
covs.append(cov)
if not covs:
continue
cov = sandy.EnergyCov.sum_covs(*covs)
idx = pd.MultiIndex.from_product(
[[mat], [mt], list(cov.data.index)],
names=["MAT", "MT", "E"],
)
idx1 = pd.MultiIndex.from_product(
[[mat], [mt1], list(cov.data.columns)],
names=["MAT1", "MT1", "E1"],
)
df = pd.DataFrame(cov.data.values, index=idx, columns=idx1) \
.stack(level=["MAT1", "MT1", "E1"]) \
.rename("VAL") \
.reset_index()
data.append(df)
if not data:
return None
df = pd.concat(data)
return sandy.CategoryCov.from_stack(df, index=["MAT", "MT", "E"],
columns=["MAT1", "MT1", "E1"],
values='VAL')

def _get_lpc_cov(self):
return

def _get_edistr_cov(self):
return

def get_cov(self, process_mf=[], mf=None, mt=None, **kwds_njoy):
"""
Get covariance matrices for nubar/xs/energy distribution/Legendre
polynomial coefficients from Endf6 or process with NJOY.

Parameters
----------
process_mf : `int` or `list`, optional
DESCRIPTION. The default is [].
mf : `int` or `list`, optional
MF number. The default is None.
**kwds_njoy : `dict`
Extra arguments for proccessing the covariance matrix with NJOY.

Returns
-------
cov : `dict` or `sandy.CategoryCov`
Dictionary containing the covariance matrix of different mf data.
If only one mf is introduced, the method returns the
`sandy.CategoryCov` object directly.

Examples
--------
Nubar covariance matrix from endf6 file:
>>> endf6 = sandy.get_endf6_file('jeff_33','xs', 922380)
>>> xs_cov = sandy.XsCov.from_endf6(endf6)
>>> nubar_cov = endf6.get_cov(mf=31)
>>> assert (nubar_cov.data.values == xs_cov.loc[(9237, 456), (9237, 456)].values).all().all()
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved

Xs covariance matrix from endf6 file:
>>> endf_xs_cov = endf6.get_cov(mf=33)
>>> size = len(endf_xs_cov.data)
>>> assert (endf_xs_cov.data.values == xs_cov.iloc[:size, :size].values).all().all()

Both together:
>>> endf_xs_cov = endf6.get_cov(mf=[31, 33])
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
>>> size = len(endf_xs_cov[33].data)
>>> assert (endf_xs_cov[33].data.values == xs_cov.iloc[:size, :size].values).all().all()
"""
# Select the mf:
if mf is not None:
listmf_ = [mf] if isinstance(mf, int) else mf
else:
listmf_ = list(self.to_series().index.get_level_values("MF")
.intersection([31, 33, 34, 35]))
list_process = [process_mf] if isinstance(process_mf, int) else process_mf
# Dividing the list information:
listmf_ = list(pd.Index(listmf_).difference(pd.Index(list_process)))
if len(list_process) != 0:
kwds_njoy["nubar"] = True if 31 in list_process else False
kwds_njoy["xs"] = True if 33 in list_process else False
kwds_njoy["chi"] = True if 35 in list_process else False
kwds_njoy["mubar"] = True if 34 in list_process else False
out = self.get_errorr(**kwds_njoy)
cov = out.get_cov(mf=list_process)
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
if len(listmf_) != 0:
cov = {} if len(list_process) == 0 else cov
if 31 in listmf_:
cov[31] = self._get_xs_cov(mf=[31])
if 33 in listmf_:
cov[33] = self._get_xs_cov(mf=[33])
AitorBengoechea marked this conversation as resolved.
Show resolved Hide resolved
if 34 in listmf_:
cov[34] = self._get_lpc_cov()
if 35 in listmf_:
cov[35] = self._get_edistr_cov()
# If only one mf is calculated, the return is directly the `CategoryCov` object
if len(cov) == 1:
[(key, cov)] = cov.items()
return cov