Sandwich tutorial

This notebook combines Sensitivity, Covariance, and Sandwich. It demonstrates similarity without covariance data, and uncertainty plus representativity when sandy is available to load the ERRORR files.

[1]:

import os
import sys
import tempfile
from pathlib import Path

os.environ.setdefault("MPLCONFIGDIR", str(Path(tempfile.gettempdir()) / "matplotlib"))

def find_project_root(start: Path | None = None) -> Path:
    """Find the pyNDUS repository root."""
    start = (start or Path.cwd()).resolve()

    for directory in (start, *start.parents):
        if (
            (directory / "setup.py").is_file()
            and (directory / "src" / "pyNDUS").is_dir()
            and (directory / "docs" / "tutorials").is_dir()
        ):
            return directory

    raise FileNotFoundError(
        f"Could not find the pyNDUS repository root starting from {start}"
    )


PROJECT_ROOT = find_project_root()

SRC_ROOT = PROJECT_ROOT / "src"
if str(SRC_ROOT) not in sys.path:
    sys.path.insert(0, str(SRC_ROOT))

DATA_ROOT = PROJECT_ROOT / "docs" / "tutorials"
SENS_ROOT = DATA_ROOT / "MG_sensitivities"
COV_ROOT = DATA_ROOT / "MG_covariances"

assert SENS_ROOT.is_dir(), f"Sensitivity directory not found: {SENS_ROOT}"
assert COV_ROOT.is_dir(), f"Covariance directory not found: {COV_ROOT}"

[2]:

import importlib.util

import numpy as np
import pandas as pd

from pyNDUS import Sensitivity, Sandwich

HAS_SANDY = importlib.util.find_spec("sandy") is not None
if HAS_SANDY:
    from pyNDUS import Covariance
else:
    Covariance = None
    print("sandy is not installed; covariance-weighted sandwich examples are skipped.")

Sensitivity inputs

[3]:

godiva = Sensitivity(SENS_ROOT / "serpent" / "godiva_sens0.m")
jezebel = Sensitivity(SENS_ROOT / "serpent" / "jezebel_sens0.m")

assert godiva.reader == "serpent"
assert jezebel.reader == "serpent"
assert godiva.n_groups == jezebel.n_groups == 33

pd.DataFrame({
    "godiva": [", ".join(godiva.zais.keys()), godiva.n_MTs],
    "jezebel": [", ".join(jezebel.zais.keys()), jezebel.n_MTs],
}, index=["nuclides", "n_MTs"])

Reading /Users/nicoloabrate/mycodes/pyNDUS/docs/tutorials/MG_sensitivities/serpent/godiva_sens0.m
  - done
Reading /Users/nicoloabrate/mycodes/pyNDUS/docs/tutorials/MG_sensitivities/serpent/jezebel_sens0.m
  - done

[3]:

	godiva	jezebel
nuclides	U-235, U-238, U-234	Pu-239, Pu-240, Pu-241, Ga-69, Ga-71
n_MTs	55	55

Similarity

Similarity does not need covariance data. Each matrix element is normalized with exactly the two MT sensitivity vectors in its numerator; no global norm over the other reactions is used. Here the same Jezebel object is used on both sides, so every non-zero diagonal element must be one. The off-diagonal elements compare different MT profiles and may differ from one. For a real cross-system comparison, choose two sensitivity objects with a common ZA/MT set.

[4]:

similarity = Sandwich(jezebel, sens2=jezebel, similarity=True, list_resp=["keff"],
                      list_za=[942390], list_MTs=[18, 102],
                      )

assert similarity.calculation_type == "similarity"
assert not similarity.similarity.empty
for mt in (18, 102):
    value = similarity.similarity.loc[("keff", "Pu-239", mt), mt]
    assert np.isclose(value.nominal_value, 1.0)
similarity.similarity

[4]:

		MT_col	18	102
RESPONSE	ZA	MT_row
keff	Pu-239	18	0.99999999999999977796+/-0.00000000000000000004	-0.5986+/-0.0005
keff	Pu-239	102	-0.5986+/-0.0005	1.0+/-0

Covariance input for uncertainty and representativity

[5]:

LIB = "endfb_80"
EGRID = "ECCO-33"
TEMPERATURE = 300

if HAS_SANDY:
    jezebel_zaids = [942390]
    jezebel_cov = {
        zaid: Covariance(zaid, TEMPERATURE, group_structure=jezebel.group_structure_as("eV"),
                         energy_unit="eV", egridname=EGRID, lib=LIB,
                         database=True, cwd=COV_ROOT,
                         )
        for zaid in jezebel_zaids
    }

    assert set(jezebel_cov) == set(jezebel_zaids)
    assert "errorr33" in jezebel_cov[942390].MFs2MTs

Uncertainty on Jezebel

The sandwich products are variance contributions. With uncertainty_output="signed_sqrt", each displayed matrix term is \(\operatorname{sign}(q_{ij})\sqrt{|q_{ij}|}\), so negative covariance contributions remain negative. These signed roots are not additive: the total standard deviation shown separately is obtained as \(\sqrt{\sum_{ij}q_{ij}}\) from the original contributions.

[6]:

if HAS_SANDY:
    uncertainty = Sandwich(jezebel, covmat=jezebel_cov, list_resp=["keff"],
                           list_mat=["total"], list_za=[942390], list_MTs=[2, 16, 18, 102],
                           list_MFs=[33], include_MF=True,
                           uncertainty_output="signed_sqrt",)

    assert uncertainty.calculation_type == "uncertainty"
    assert not uncertainty.uncertainty.empty
    # Signed square roots of the individual MT-MT contributions.
    print("Uncertainty as std. dev. in [pcm]:")
    display(uncertainty.uncertainty*1E5)
    # The physically additive variance matrix and the total standard deviation.
    print("Uncertainty as variance:")
    display(uncertainty.uncertainty_variance)
    print("Total uncertainty in pcm:")
    display(uncertainty.uncertainty_standard_deviation*1E5)

Uncertainty as std. dev. in [pcm]:

				MT_col	2	18	102
RESPONSE	MATERIAL	ZA	MF	MT_row
keff	total	Pu-239	33	2	483.6+/-0.7	-159.07+/-0.13	25.603+/-0.022
				18	-159.07+/-0.13	905.98+/-0.09	0.0+/-0
				102	25.603+/-0.022	0.0+/-0	74.40+/-0.08

Uncertainty as variance:

				MT_col	2	18	102
RESPONSE	MATERIAL	ZA	MF	MT_row
keff	total	Pu-239	33	2	(2.338+/-0.007)e-05	(-2.530+/-0.004)e-06	(6.555+/-0.011)e-08
				18	(-2.530+/-0.004)e-06	(8.2079+/-0.0016)e-05	0.0+/-0
				102	(6.555+/-0.011)e-08	0.0+/-0	(5.536+/-0.011)e-07

Total uncertainty in pcm:

		variance	standard_deviation
RESPONSE	MATERIAL
keff	total	10.109+/-0.007	1005.4+/-0.4

Representativity

The same Jezebel sensitivity is used as reference and application. This is a compact self-test: the result should be defined wherever the covariance and sensitivities overlap. For this calculation we let list_za default to all Jezebel nuclides, consistently with the normalization terms used internally.

representativity_total is the physical coefficient, obtained by summing all normalized ZA/MF/MT-pair contributions. The matrix and long table are retained only to diagnose how individual nuclear-data terms contribute to that total.

[8]:

if HAS_SANDY:
    representativity = Sandwich(jezebel, sens2=jezebel, covmat=jezebel_cov,
                                representativity=True, list_resp=["keff"],
                                list_za=None,
                                list_MTs=[18, 102], list_MFs=[33], include_MF=True,
                                )

    assert representativity.calculation_type == "representativity"
    assert not representativity.representativity.empty
    assert np.isclose(representativity.representativity_total.loc["keff"].nominal_value, 1.0)
    print("Total representativity between Jezebel and Jezebel itself:")
    display(representativity.representativity_total)
    print("Contribution to representativity:")
    display(representativity.representativity)
    display(representativity.representativity_table.head())

Total representativity between Jezebel and Jezebel itself:

RESPONSE
keff    1.00000+/-0.00020
Name: representativity_total, dtype: object

Contribution to representativity:

			MT_col	18	102
RESPONSE	ZA	MF	MT_row
keff	Pu-239	33	18	0.99330+/-0.00020	0.0+/-0
keff	Pu-239	33	102	0.0+/-0	0.006699+/-0.000010

					value
RESPONSE	ZA	MF	MT_row	MT_col
keff	Pu-239	33	18	18	0.99330+/-0.00020
			18	102	0.0+/-0
			102	18	0.0+/-0
			102	102	0.006699+/-0.000010