Skip to content

Commit

Permalink
Add sampling mf35 (#356)
Browse files Browse the repository at this point in the history 
* update

* update

---------

Co-authored-by: Luca Fiorito <[email protected]>
  • Loading branch information
@luca-fiorito-11
luca-fiorito-11 and Luca Fiorito authored Jan 10, 2025
1 parent bfa5a0a commit 8d232b0
Show file tree
Hide file tree
Showing 4 changed files with 258 additions and 12 deletions.
114 changes: 114 additions & 0 deletions notebooks/notebook_sampling_mf35_testing.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "2255de27-2a19-4cf9-9c26-122962881e06",
"metadata": {},
"source": [
"# Testing sampling of MF35 for U235"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a56e847e-3d93-446a-ac9a-e05cd7ff894f",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"import pandas as pd\n",
"import sandy"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a51b18d-3d10-435d-b737-14a8f8f3ad00",
"metadata": {},
"outputs": [],
"source": [
"endf6 = sandy.get_endf6_file(\"jeff_33\", \"xs\", 922350)\n",
"endf6.to_file(\"U235.jeff33\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a0477c60-8ed7-40f2-a0f9-b14ba763bc3d",
"metadata": {},
"outputs": [],
"source": [
"sandy.sampling.run(\"U235.jeff33 --mf 35 --samples 2 --seed35 5 --loglevel error\".split())"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1df9b6b3-6725-4a43-8f0d-647c9830b38c",
"metadata": {},
"outputs": [],
"source": [
"smp = sandy.Samples.from_excel(\"PERT_92235_MF35.xlsx\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1469dadc-3ec4-4aab-a4cf-e17a45c9c62f",
"metadata": {},
"outputs": [],
"source": [
"# Extract energy distributions from original and perturbed files\n",
"edistr = sandy.Edistr.from_endf6(endf6)\n",
"edistr0 = sandy.Edistr.from_endf6(sandy.Endf6.from_file(\"92235_0.endf6\"))\n",
"edistr1 = sandy.Edistr.from_endf6(sandy.Endf6.from_file(\"92235_1.endf6\"))\n",
"\n",
"# Calculate ratios, this should be compared to the random samples\n",
"edistr.data[\"RATIO0\"] = edistr0.data.VALUE / edistr.data.VALUE\n",
"edistr.data[\"RATIO1\"] = edistr1.data.VALUE / edistr.data.VALUE"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e6b533bf-1d3a-41fb-8095-02271774a21d",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"sns.lineplot(data=edistr.data, x=\"EOUT\", y=\"RATIO0\", errorbar=\"sd\", ax=ax, color=\"dodgerblue\")\n",
"smp.get_eright()[0].reset_index().plot(x=\"ERIGHT\", y=0, drawstyle=\"steps-pre\", ax=ax, color=\"blue\", ls=\"--\", lw=.8)\n",
"\n",
"sns.lineplot(data=edistr.data, x=\"EOUT\", y=\"RATIO1\", errorbar=\"sd\", ax=ax, color=\"tomato\")\n",
"smp.get_eright()[1].reset_index().plot(x=\"ERIGHT\", y=1, drawstyle=\"steps-pre\", ax=ax, color=\"red\", ls=\"--\", lw=.8)\n",
"\n",
"ax.set(xlim=[0,3e7], ylabel='relative perturbation', xlabel='energy / $eV$', title=\"Compare perturbed ENDF-6 files with random samples \")\n",
"ax.legend(title=\"SAMPLE\")\n",
"fig.tight_layout()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "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.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
36 changes: 28 additions & 8 deletions sandy/endf6.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3220,21 +3220,41 @@ def endf6_perturb_worker(e6, pendf, ismp,
if plpc is not None:
pass

# apply edistr perturbation
# Apply energy distribution (edistr) perturbation
if pchi is not None:
# Applies the same perturbation to all incident particle energies and K
# Applies the same perturbation to all incident particle energies (EIN) and K
edistr_pert = []

# Group data by EIN and K for processing
for (ein, k), df in sandy.Edistr.from_endf6(endf6_pert).data.groupby(['EIN', 'K']):
# Prepare dummy energy distribution data as a xs object
dummy_xs = sandy.Xs(
df.rename({"EOUT": "E"}, axis=1).set_index(["MAT","MT"])[["E","VALUE"]].pivot(columns="E").T.droplevel(level=0)
df.rename({"EOUT": "E"}, axis=1)
.set_index(["MAT","MT"])[["E","VALUE"]]
.pivot(columns="E").T.droplevel(level=0)
)

# Apply perturbation to dummy energy distribution
dummy_xs_pert = sandy.xs.xs_perturb_worker(dummy_xs, ismp, pchi, verbose=verbose)
edistr_pert.append(
dummy_xs_pert.data.stack([1, 0]).to_frame().reset_index().rename({"E": "EOUT", 0: "VALUE"}, axis=1).assign(K=k, EIN=ein)[["MAT", "MT", "K", "EIN", "EOUT", "VALUE"]] # sort columns to match Edistr.data

# Transform xs data into edistr data and append perturbed data
perturbed_data = (
dummy_xs_pert.data.stack([1, 0], future_stack=True) # Use future_stack=True to adopt the new behavior
.to_frame()
.reset_index()
.rename({"E": "EOUT", 0: "VALUE"}, axis=1)
.assign(K=k, EIN=ein)
[["MAT", "MT", "K", "EIN", "EOUT", "VALUE"]]
)
edistr_pert.append(perturbed_data)

# Combine and normalize perturbed data, then update ENDF6
endf6_pert = (
sandy.Edistr(pd.concat(edistr_pert, ignore_index=True))
.normalize()
.to_endf6(endf6_pert)
.update_intro()
)
endf6_pert = sandy.Edistr(
pd.concat(edistr_pert, ignore_index=True)
).normalize().to_endf6(endf6_pert).update_intro()

# apply xs perturbation
if pxs is not None:
Expand Down
112 changes: 112 additions & 0 deletions sandy/samples.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,10 @@ class Samples():
Return correlation matrix of samples.
get_cov
Return covariance matrix of samples.
get_eleft
Replace energy intervals with left bounds.
get_eright
Replace energy intervals with right bounds.
get_mean
Return mean vector of samples.
get_std
Expand Down Expand Up @@ -169,6 +173,114 @@ def get_corr(self):
def get_cov(self):
return self.data.T.cov()

def get_eright(self):
"""
Create a new DataFrame with an updated index, where the 'E' level in the MultiIndex
(representing energy intervals) is replaced by the right bound of each interval,
and the original 'E' level is dropped. The new level is named 'ERIGHT'.
The method performs the following steps:
1. Extracts the right bounds of the energy intervals in the 'E' level of the MultiIndex.
2. Constructs a new MultiIndex that includes the right bound values as a new level,
and removes the original 'E' level.
3. Creates a copy of the current DataFrame, replacing the index with the newly
constructed MultiIndex.
4. Returns the modified DataFrame with the updated index.
Returns
-------
pd.DataFrame
A new DataFrame with the same data as the original, but with an updated MultiIndex
where the 'E' level has been replaced by the 'ERIGHT' level.
Notes
-----
- The method assumes that the original MultiIndex has a level named 'E', which contains
intervals.
- The method will drop the 'E' level from the MultiIndex and replace it with the right bounds
of the intervals, renaming the new level to 'ERIGHT'.
Examples
--------
>>> import sandy
>>> import numpy as np
>>> endf6 = sandy.get_endf6_file('jeff_33', 'xs', 10010)
>>> smps1 = endf6.get_perturbations(1, njoy_kws=dict(err=1, chi=False, mubar=False, nubar=False, errorr33_kws=dict(mt=2, ek=[1, 2, 3])))[33]
>>> np.testing.assert_array_equal(smps1.get_eright().index.get_level_values("ERIGHT"), [2, 3])
>>> np.testing.assert_array_equal(smps1.get_eright().values, smps1.data.values)
"""
multi_index = self.data.index

# right bound of the energy intervals, which will become the new level
new_level = [x.right for x in multi_index.get_level_values("E")]

# Add the new level and drop the EnergyInterval index
new_multi_index = pd.MultiIndex.from_tuples(
[(*old, new) for old, new in zip(multi_index, new_level)],
names=[*multi_index.names, "ERIGHT"]
).droplevel("E")

# Create a copy of the samples with the new index
copy = self.data.copy()
copy.index = new_multi_index

# Return a dataframe, not the Samples object
return copy

def get_eleft(self):
"""
Create a new DataFrame with an updated index, where the 'E' level in the MultiIndex
(representing energy intervals) is replaced by the left bound of each interval,
and the original 'E' level is dropped. The new level is named 'ELEFT'.
The method performs the following steps:
1. Extracts the left bounds of the energy intervals in the 'E' level of the MultiIndex.
2. Constructs a new MultiIndex that includes the left bound values as a new level,
and removes the original 'E' level.
3. Creates a copy of the current DataFrame, replacing the index with the newly
constructed MultiIndex.
4. Returns the modified DataFrame with the updated index.
Returns
-------
pd.DataFrame
A new DataFrame with the same data as the original, but with an updated MultiIndex
where the 'E' level has been replaced by the 'ELEFT' level.
Notes
-----
- The method assumes that the original MultiIndex has a level named 'E', which contains
intervals.
- The method will drop the 'E' level from the MultiIndex and replace it with the left bounds
of the intervals, renaming the new level to 'ELEFT'.
Examples
--------
>>> import sandy
>>> import numpy as np
>>> endf6 = sandy.get_endf6_file('jeff_33', 'xs', 10010)
>>> smps1 = endf6.get_perturbations(1, njoy_kws=dict(err=1, chi=False, mubar=False, nubar=False, errorr33_kws=dict(mt=2, ek=[1, 2, 3])))[33]
>>> np.testing.assert_array_equal(smps1.get_eleft().index.get_level_values("ELEFT"), [1, 2])
>>> np.testing.assert_array_equal(smps1.get_eleft().values, smps1.data.values)
"""
multi_index = self.data.index

# right bound of the energy intervals, which will become the new level
new_level = [x.left for x in multi_index.get_level_values("E")]

# Add the new level and drop the EnergyInterval index
new_multi_index = pd.MultiIndex.from_tuples(
[(*old, new) for old, new in zip(multi_index, new_level)],
names=[*multi_index.names, "ELEFT"]
).droplevel("E")

# Create a copy of the samples with the new index
copy = self.data.copy()
copy.index = new_multi_index

# Return a dataframe, not the Samples object
return copy

def get_std(self):
return self.data.std(axis=1).rename("STD")

Expand Down
8 changes: 4 additions & 4 deletions sandy/sampling.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,10 +85,10 @@ def parse(iargs=None):

parser.add_argument('--mf',
type=int,
default=[31, 33,],
default=[31, 33, 35],
action='store',
nargs="+",
metavar="{31,33}",
metavar="{31,33,35}",
help="draw samples only from the selected MF sections "
"(default = keep all)")

Expand Down Expand Up @@ -253,7 +253,7 @@ def multi_run(foo):
>>> file = "942410.jeff33"
>>> sandy.get_endf6_file("jeff_33", "xs", 942410).to_file(file)
>>> cl = f"{file}" + " --samples 2 -O {SMP}-{ZAM} --seed33 1 --seed31 1 --mt33 2"
>>> cl = f"{file}" + " --samples 2 -O {SMP}-{ZAM} --seed33 1 --seed31 1 --seed35 1 --mt33 2"
>>> sandy.sampling.run(cl.split())
Now, let's interrupt the process after that the perturbations are
Expand Down Expand Up @@ -409,7 +409,7 @@ def run(iargs):
nubar = bool(31 in iargs.mf) and (31 in endf6.mf)
xs = bool(33 in iargs.mf) and (33 in endf6.mf or 32 in endf6.mf) # this handles together MF32 and MF33
mubar = False
chi = False
chi = bool(35 in iargs.mf) and (35 in endf6.mf)
errorr_kws = dict(
verbose=iargs.debug,
err=err_errorr,
Expand Down

0 comments on commit 8d232b0

Please sign in to comment.