import pandas as pd
import numpy as np
from numpy.random import multivariate_normal, normal
from numpy.linalg import cholesky, svd
from scipy.sparse.linalg import lsqr
from itertools import accumulate
from sklearn.linear_model import LinearRegression
from scipy.linalg import lstsq
from scipy.optimize import nnls
import pytest
import matplotlib.pyplot as plt
import seaborn as sns
import sandy
import serpentTools as sts
sens = sts.read("Godiva.i_sens0.m")
INFO: Reading Godiva.i_sens0.m INFO: - done
nzai = sens.zais[922350]
sensitivity = {}
e = sens.energies
npert = sens.perts['capture xs']
s = sens.sensitivities["keff"][0, nzai, npert].T[0]
sensitivity["capture xs"] = pd.Series(s, index=e[1:])
npert = sens.perts['fission xs']
s = sens.sensitivities["keff"][0, nzai, npert].T[0]
sensitivity["fission xs"] = pd.Series(s, index=e[1:])
npert = sens.perts['ela scatt xs']
s = sens.sensitivities["keff"][0, nzai, npert].T[0]
sensitivity["ela scatt xs"] = pd.Series(s, index=e[1:])
npert = sens.perts['inl scatt xs']
s = sens.sensitivities["keff"][0, nzai, npert].T[0]
sensitivity["inl scatt xs"] = pd.Series(s, index=e[1:])
tape = sandy.get_endf6_file("jeff_33", "xs", 922350)
mt = [2, 4, 18, 102]
err = tape.get_errorr(err=1, xs=True, nubar=False, chi=False, mubar=False, errorr33_kws=dict(mt=mt, ek=e), verbose=False)["errorr33"]
cov = err.get_cov()
INFO: Zero or no temperature was requested, NJOY processing will stop after RECONR. If you want to process 0K cross sections use `temperature=0.1`.
njoy 2016.74 12Jan24 03/03/24 21:48:34 ***************************************************************************** moder... 0.0s reconr... 0.1s ---message from rdf2bw---calculation of angular distribution not installed. moder... 5.9s errorr... 6.4s ---message from grpav---mf 3 mt 4 has threshold gt highest union energy. ---message from grpav---mf 3 mt 16 has threshold gt highest union energy. ---message from grpav---mf 3 mt 17 has threshold gt highest union energy. ---message from grpav---mf 3 mt 37 has threshold gt highest union energy. processing mat 9228 --------------------- 92-U -235 IRSN-CEA EVAL-DEC14 IRSN-CEA DAM/DEN COLLAB. ---message from rdgout---mf 3 mt 4 not found. calculation continued with sigma=0 ---message from rdgout---mf 3 mt 4 not found. calculation continued with sigma=0 covariances calculated for 4 reactions and 1006 groups 7.3s ---message from rdgout---mf 3 mt 4 not found. calculation continued with sigma=0 12.3s *****************************************************************************
/home/runner/.local/lib/python3.12/site-packages/sandy-1.1-py3.12.egg/sandy/core/endf6.py:963: DtypeWarning: Columns (1,2,3) have mixed types. Specify dtype option on import or set low_memory=False. df = pd.read_csv(
C = cov.data.copy()
#print(f"condition number original matrix: {np.linalg.cond(C):>10}")
C = pd.DataFrame(C.values + np.diag(np.diag(C.values) * 0.5 / 100 + 0.0001), index=C.index, columns=C.columns)
#print(f"condition number adjusted matrix: {np.linalg.cond(C):>10}")
cov_ = sandy.CategoryCov(C)
U, S, V = svd(C)
Lambda12 = np.diag(np.sqrt(S))
fig, ax = plt.subplots(figsize=(6, 6))
sns.heatmap(cov_.get_corr().data, ax=ax, cmap="bwr", vmin=-1, vmax=1)
fig.tight_layout()
M = cov.data.shape[0] # number of parameters
frac = S / S.sum()
acc = np.array(list(accumulate(frac)))
r = acc[acc < 0.98].size + 1
N = 5000 # number of samples
print(f"sample size: {N:>10}\nnumber of parameters: {M:>10}\nreduced number of parameters: {r:>10}")
sample size: 5000 number of parameters: 4000 reduced number of parameters: 3450
# non correlated, standardized sample
X_ = normal(loc=1., scale=1., size=N * M).reshape(N, M)
# correlated sample
X = X_ @ (U @ Lambda12).T
def model1(x): return sensitivity["capture xs"] @ x
def model2(x): return sensitivity["fission xs"] @ x
def model3(x): return sensitivity["ela scatt xs"] @ x
def model4(x): return sensitivity["inl scatt xs"] @ x
fig, ax = plt.subplots(figsize=(8, 5))
model = model1
data = {}
for fv in np.array([0.8, 0.86, 0.9, 0.95]):
n = acc[acc < fv].size + 1
Y = model(X[:n, (e.size - 1)*2:(e.size - 1)*3].T)
Z = lstsq(X[:n, (e.size - 1)*2:(e.size - 1)*3], Y)[0]
data[f"N={n}, FV={fv}"] = Z
# Y = model(X[:, (e.size - 1)*0:(e.size - 1)*1].T)
# Z = lstsq(X[:, (e.size - 1)*0:(e.size - 1)*1], Y)[0]
# data[f"N={N}"] = Z
pd.DataFrame(data, index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre")
sensitivity["capture xs"].plot(ax=ax, label="sensitivity", ls="--", drawstyle="steps-pre")
ax.legend()
ax.set(xlim=[1e-3, 20], title="fssion xs")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(8, 5))
model = model2
data = {}
for fv in np.array([0.8, 0.86, 0.9, 0.95]):
n = acc[acc < fv].size + 1
Y = model(X[:n, (e.size - 1)*3:(e.size - 1)*4].T)
Z = lstsq(X[:n, (e.size - 1)*3:(e.size - 1)*4], Y)[0]
data[f"N={n}, FV={fv}"] = Z
# Y = model(X[:, (e.size - 1)*1:(e.size - 1)*2].T)
# Z = lstsq(X[:, (e.size - 1)*1:(e.size - 1)*2], Y)[0]
# data[f"N={N}"] = Z
pd.DataFrame(data, index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre")
sensitivity["fission xs"].plot(ax=ax, label="sensitivity", ls="--", drawstyle="steps-pre")
ax.legend()
ax.set(xlim=[1e-3, 20], title="capture xs")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(8, 5))
model = model3
data = {}
for fv in np.array([0.8, 0.86, 0.9, 0.95]):
n = acc[acc < fv].size + 1
Y = model(X[:n, (e.size - 1)*0:(e.size - 1)*1].T)
Z = lstsq(X[:n, (e.size - 1)*0:(e.size - 1)*1], Y)[0]
data[f"N={n}, FV={fv}"] = Z
# Y = model(X[:, (e.size - 1)*2:(e.size - 1)*3].T)
# Z = lstsq(X[:, (e.size - 1)*2:(e.size - 1)*3], Y)[0]
# data[f"N={N}"] = Z
pd.DataFrame(data, index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre")
sensitivity["ela scatt xs"].plot(ax=ax, label="sensitivity", ls="--", drawstyle="steps-pre")
ax.legend()
ax.set(xlim=[1e-3, 20], title="ela scatt xs")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(8, 5))
model = model4
data = {}
for fv in np.array([0.8, 0.86, 0.9, 0.95]):
n = acc[acc < fv].size + 1
Y = model(X[:n, (e.size - 1)*1:(e.size - 1)*2].T)
Z = lstsq(X[:n, (e.size - 1)*1:(e.size - 1)*2], Y)[0]
data[f"N={n}, FV={fv}"] = Z
# Y = model(X[:, (e.size - 1)*3:(e.size - 1)*4].T)
# Z = lstsq(X[:, (e.size - 1)*3:(e.size - 1)*4], Y)[0]
# data[f"N={N}"] = Z
pd.DataFrame(data, index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre")
sensitivity["inl scatt xs"].plot(ax=ax, label="sensitivity", ls="--", drawstyle="steps-pre")
ax.legend()
ax.set(xlim=[1e-3, 20], title="inl scatt xs")
fig.tight_layout()
def model(x):
sens_all = pd.concat([sensitivity["ela scatt xs"], sensitivity["inl scatt xs"], sensitivity["fission xs"], sensitivity["capture xs"]], ignore_index=True).values
return sens_all @ x
fig, ax = plt.subplots(figsize=(10, 6))
model = model
fv = 0.80
n = acc[acc < fv].size + 1
Y = model(X[:n, :].T)
Z = lstsq(X[:n, :], Y)[0]
pd.Series(Z[(e.size - 1)*0:(e.size - 1)*1], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity ela scatt xs")
pd.Series(Z[(e.size - 1)*1:(e.size - 1)*2], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity inl scatt xs")
pd.Series(Z[(e.size - 1)*2:(e.size - 1)*3], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity capture xs")
pd.Series(Z[(e.size - 1)*3:(e.size - 1)*4], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity fission xs")
sensitivity["capture xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["fission xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["ela scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["inl scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
ax.legend()
ax.set(xlim=[1e-3, 20], title=f"fraction of variance: {fv} (N={n})")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(10, 6))
model = model
fv = 0.90
n = acc[acc < fv].size + 1
Y = model(X[:n, :].T)
Z = lstsq(X[:n, :], Y)[0]
pd.Series(Z[(e.size - 1)*0:(e.size - 1)*1], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity ela scatt xs")
pd.Series(Z[(e.size - 1)*1:(e.size - 1)*2], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity inl scatt xs")
pd.Series(Z[(e.size - 1)*2:(e.size - 1)*3], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity capture xs")
pd.Series(Z[(e.size - 1)*3:(e.size - 1)*4], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity fission xs")
sensitivity["capture xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["fission xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["ela scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["inl scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
ax.legend()
ax.set(xlim=[1e-3, 20], title=f"fraction of variance: {fv} (N={n})")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(10, 6))
model = model
fv = 0.95
n = acc[acc < fv].size + 1
Y = model(X[:n, :].T)
Z = lstsq(X[:n, :], Y)[0]
pd.Series(Z[(e.size - 1)*0:(e.size - 1)*1], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity ela scatt xs")
pd.Series(Z[(e.size - 1)*1:(e.size - 1)*2], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity inl scatt xs")
pd.Series(Z[(e.size - 1)*2:(e.size - 1)*3], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity capture xs")
pd.Series(Z[(e.size - 1)*3:(e.size - 1)*4], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity fission xs")
sensitivity["capture xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["fission xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["ela scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["inl scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
ax.legend()
ax.set(xlim=[1e-3, 20], title=f"fraction of variance: {fv} (N={n})")
fig.tight_layout()
fig, ax = plt.subplots(figsize=(10, 6))
model = model
fv = 0.99
n = acc[acc < fv].size + 1
Y = model(X[:n, :].T)
Z = lstsq(X[:n, :], Y)[0]
pd.Series(Z[(e.size - 1)*0:(e.size - 1)*1], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity ela scatt xs")
pd.Series(Z[(e.size - 1)*1:(e.size - 1)*2], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity inl scatt xs")
pd.Series(Z[(e.size - 1)*2:(e.size - 1)*3], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity capture xs")
pd.Series(Z[(e.size - 1)*3:(e.size - 1)*4], index=e[1:]).plot(kind="line", ax=ax, logx=True, drawstyle="steps-pre", label="sensitivity fission xs")
sensitivity["capture xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["fission xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["ela scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
sensitivity["inl scatt xs"].plot(ax=ax, ls="--", lw=2, drawstyle="steps-pre", alpha=.5)
ax.legend()
ax.set(xlim=[1e-3, 20], title=f"fraction of variance: {fv} (N={n})")
fig.tight_layout()