Merge c5ade92937 into d3bc1669d3

openmc/data/mf33_njoy.py

@@ -0,0 +1,385 @@
+#!/usr/bin/env python3
+"""
+Minimal NJOY/ERRORR driver to produce MF=33 (multigroup covariance) as text.
+Run a short NJOY chain and return dict with keys: mat, ek, tape33.
+
+Assumptions
+-----------
+- `energy_grid_ev` is a strictly increasing sequence of group boundaries in **eV** (length G+1).
+- Output is relative covariances (IRELCO=1) on an explicit group structure (IGN=1).
+"""
+
+from __future__ import annotations
+
+import os
+import shutil
+import subprocess as sp
+from pathlib import Path
+from tempfile import TemporaryDirectory
+from typing import Any, Dict, List, Optional, Sequence, Union
+
+
+NJOY_tolerance = 0.01
+NJOY_sig0 = [1e10]
+NJOY_thermr_emax = 10.0
+
+# ------------------------- small helpers -------------------------
+
+def _read_mat_from_endf(endf_path: Path) -> int:
+    """Infer MAT from an ENDF-6 evaluation (text)."""
+    first_nonzero = None
+    with endf_path.open("r", errors="ignore") as f:
+        for line in f:
+            if len(line) < 75:
+                continue
+            try:
+                mat = int(line[66:70])
+                mf = int(line[70:72])
+                mt = int(line[72:75])
+            except ValueError:
+                continue
+            if mat > 0 and first_nonzero is None:
+                first_nonzero = mat
+            if mat > 0 and mf == 1 and mt == 451:
+                return mat
+    if first_nonzero is None:
+        raise ValueError(f"Could not infer MAT from ENDF file: {endf_path}")
+    return first_nonzero
+
+def _endf_has_mf33(endf_path: Path, mat: int) -> bool:
+    """Return True if the ENDF tape contains any MF=33 section for this MAT."""
+    with endf_path.open("r", errors="ignore") as f:
+        for line in f:
+            if len(line) < 75:
+                continue
+            try:
+                m  = int(line[66:70])
+                mf = int(line[70:72])
+            except ValueError:
+                continue
+            if m == mat and mf == 33:
+                return True
+    return False
+
+def _validate_energy_grid_ev(ek: Sequence[float]) -> List[float]:
+    ek = [float(x) for x in ek]
+    if len(ek) < 2:
+        raise ValueError("Energy grid must have at least 2 boundaries (G+1).")
+    if ek[0] <= 0.0:
+        raise ValueError(
+            f"Energy grid lower boundary must be positive (got {ek[0]:g} eV). "
+            f"ERRORR cannot integrate from zero energy. Use a small positive "
+            f"value like 1e-5 eV instead."
+        )
+    for i in range(1, len(ek)):
+        if not (ek[i] > ek[i-1]):
+            raise ValueError("Energy grid boundaries must be strictly increasing (in eV).")
+    return ek
+
+# ------------------------- NJOY deck builder -------------------------
+
+def _moder_input(nin: int, nout: int) -> str:
+    return f"moder\n{nin:d} {nout:d} /\n"
+
+def _reconr_input(endfin, pendfout, mat, err=NJOY_tolerance):
+    return (
+        "reconr\n"
+        f"{endfin:d} {pendfout:d} /\n"
+        f"'mf33'/\n"
+        f"{mat:d} 0 0 /\n"
+        f"{err:g} 0. /\n"
+        "0/\n"
+    )
+
+def _thermr_input(endfin, pendfin, pendfout, mat, temperature,
+                   err=NJOY_tolerance, emax=NJOY_thermr_emax):
+    return (
+        "thermr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"0 {mat:d} 20 1 1 0 0 1 221 0 /\n"
+        f"{temperature:.1f} /\n"
+        f"{err:g} {emax:g} /\n"
+    )
+
+def _broadr_input(endfin, pendfin, pendfout, mat, temperature, err=NJOY_tolerance):
+    return (
+        "broadr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 0 0 0. /\n"
+        f"{err:g} /\n"
+        f"{temperature:.1f} /\n"
+        "0 /\n"
+    )
+
+def _unresr_input(endfin, pendfin, pendfout, mat, temperature,
+                   sig0=NJOY_sig0):
+    nsig0 = len(sig0)
+    return (
+        "unresr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 {nsig0:d} 0 /\n"
+        f"{temperature:.1f} /\n"
+        + " ".join(f"{s:.2E}" for s in sig0) + " /\n"
+        "0 /\n"
+    )
+
+def _purr_input(endfin, pendfin, pendfout, mat, temperature,
+                 sig0=NJOY_sig0, bins=20, ladders=32):
+    nsig0 = len(sig0)
+    return (
+        "purr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 {nsig0:d} {bins:d} {ladders:d} 0 /\n"
+        f"{temperature:.1f} /\n"
+        + " ".join(f"{s:.2E}" for s in sig0) + " /\n"
+        "0 /\n"
+    )
+
+
+def _errorr_mf33_input(
+    *,
+    endfin: int,
+    pendfin: int,
+    errorrout: int,
+    mat: int,
+    temperature: float,
+    ek: Sequence[float],
+    iwt: int = 6,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> str:
+    irelco = 1 if relative else 0
+    iread  = 1 if mt is not None else 0
+    iwt_   = 1 if spectrum is not None else iwt
+    pf     = int(iprint)
+    nk     = len(ek) - 1
+
+    lines = [
+        "errorr",
+        f"{endfin:d} {pendfin:d} 0 {errorrout:d} 0 /",
+        f"{mat:d} 1 {iwt_:d} {pf:d} {irelco} /",   # ign=1 (explicit grid)
+        f"{pf:d} {temperature:.1f} /",
+        f"{iread:d} 33 {irespr:d} /",
+    ]
+
+    # MT list
+    if iread == 1:
+        mtlist = [mt] if isinstance(mt, int) else list(mt)
+        lines.append(f"{len(mtlist):d} 0 /")
+        lines.append(" ".join(str(m) for m in mtlist) + " /")
+
+    # Explicit energy grid
+    lines.append(f"{nk} /")
+    lines.append(" ".join(f"{x:.5e}" for x in ek) + " /")
+
+    # User weight spectrum (TAB1, iwt=1)
+    if iwt_ == 1 and spectrum is not None:
+        n_pairs = len(spectrum) // 2
+        lines.append(f" 0.0 0.0 0 0 1 {n_pairs:d} /")
+        lines.append(f"{n_pairs:d} 1 /")
+        for k in range(n_pairs):
+            lines.append(f"{spectrum[2*k]:.6e} {spectrum[2*k+1]:.6e} /")
+        lines.append("/")
+
+    return "\n".join(lines) + "\n"
+
+
+# ---- deck assembly ---------------------------------------------------------
+
+def _build_deck(
+    mat: int,
+    ek: Sequence[float],
+    temperature: float,
+    *,
+    err: float = NJOY_tolerance,
+    thermr: bool = False,
+    unresr: bool = False,
+    purr: bool = False,
+    iwt: int = 2,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> str:
+    e = 21      # internal formatted ENDF
+    p = e + 1   # running PENDF tape number
+    deck = _moder_input(20, -e)
+
+    # RECONR
+    deck += _reconr_input(-e, -p, mat=mat, err=err)
+
+    # BROADR
+    if temperature > 0:
+        o = p + 1
+        deck += _broadr_input(-e, -p, -o, mat=mat,
+                               temperature=temperature, err=err)
+        p = o
+
+    # THERMR (free-gas)
+    if thermr:
+        o = p + 1
+        deck += _thermr_input(0, -p, -o, mat=mat,
+                               temperature=temperature, err=err)
+        p = o
+
+    # UNRESR
+    if unresr:
+        o = p + 1
+        deck += _unresr_input(-e, -p, -o, mat=mat,
+                               temperature=temperature)
+        p = o
+
+    # PURR
+    if purr:
+        o = p + 1
+        deck += _purr_input(-e, -p, -o, mat=mat,
+                             temperature=temperature)
+        p = o
+
+    # ERRORR (MF=33)
+    deck += _errorr_mf33_input(
+        endfin=-e, pendfin=-p, errorrout=33,
+        mat=mat, temperature=temperature, ek=ek,
+        iwt=iwt, spectrum=spectrum, relative=relative,
+        irespr=irespr, mt=mt, iprint=iprint,
+    )
+    deck += "stop\n"
+    return deck
+
+
+# ---- NJOY runner -----------------------------------------------------------
+
+def _run_njoy(deck: str, endf: Path, exe: Optional[str]) -> Dict[int, str]:
+    if exe is None:
+        exe = os.environ.get("NJOY")
+        if not exe:
+            raise ValueError(
+                "NJOY executable not provided and $NJOY env var is unset."
+            )
+
+    with TemporaryDirectory() as td:
+        tmpdir = Path(td)
+        shutil.copy(endf, tmpdir / "tape20")
+        (tmpdir / "input").write_text(deck)
+
+        proc = sp.Popen(
+            exe, shell=True, cwd=str(tmpdir),
+            stdin=sp.PIPE, stdout=sp.PIPE, stderr=sp.PIPE,
+        )
+        stdout, stderr = proc.communicate(input=deck.encode())
+
+        if proc.returncode != 0:
+            # Save work dir for debugging
+            dest = Path("njoy_failed_outputs")
+            if dest.exists():
+                shutil.rmtree(dest)
+            shutil.copytree(tmpdir, dest)
+            raise RuntimeError(
+                f"NJOY failed (rc={proc.returncode}). "
+                f"Work dir saved to '{dest}'.\n"
+                f"stderr: {stderr.decode(errors='replace')[:2000]}"
+            )
+
+        tp33 = tmpdir / "tape33"
+        if not tp33.exists():
+            raise RuntimeError("NJOY ran but did not produce tape33.")
+
+        return {33: tp33.read_text()}
+
+
+# ---- public API ------------------------------------------------------------
+
+def generate_errorr_mf33(
+    endf_path: str | Path,
+    energy_grid_ev: Sequence[float],
+    *,
+    njoy_exec: Optional[str] = None,
+    mat: Optional[int] = None,
+    temperature: float = 293.6,
+    # processing chain
+    err: float = NJOY_tolerance,
+    minimal_processing: bool = False,
+    thermr: bool = False,
+    unresr: bool = True,
+    purr: bool = True,
+    # ERRORR options
+    iwt: int = 6,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> Dict[str, Any]:
+    """Run NJOY/ERRORR and return MF=33 tape33 as text.
+
+    Parameters
+    ----------
+    endf_path
+        Path to ENDF-6 evaluation (text).
+    energy_grid_ev
+        Explicit group boundaries (G+1 values) in eV.
+    njoy_exec
+        NJOY executable/command.  Falls back to $NJOY.
+    mat
+        MAT number.  Inferred from the file when omitted.
+    temperature
+        Processing temperature in K.
+    err
+        Reconstruction tolerance for RECONR / BROADR.
+    minimal_processing
+        If True (default), force thermr/unresr/purr off.
+    thermr, unresr, purr
+        Individual module toggles.  All ignored when
+        ``minimal_processing=True``.
+    iwt
+        Weight function option (default 2 = constant).
+        Ignored when *spectrum* is given.
+    spectrum
+        User weight function as flat [E0, W0, E1, W1, …].
+        Forces iwt=1.
+    relative
+        Produce relative covariances (default True).
+    irespr
+        Resonance-parameter covariance processing
+        (0 = area sensitivity, 1 = 1%% sensitivity).
+    mt
+        Restrict ERRORR to specific reaction numbers.
+    iprint
+        NJOY print flag.
+
+    Returns
+    -------
+    dict
+        Keys: "mat", "ek" (validated eV grid), "tape33" (text).
+    """
+    endf_path = Path(endf_path).resolve()
+    if not endf_path.exists():
+        raise FileNotFoundError(endf_path)
+
+    ek = _validate_energy_grid_ev(energy_grid_ev)
+    if mat is None:
+        mat = _read_mat_from_endf(endf_path)
+
+    if not _endf_has_mf33(endf_path, mat):
+        raise ValueError(
+            f"ENDF file {endf_path.name} (MAT={mat}) contains no MF=33 "
+            f"covariance data; nothing for ERRORR to process."
+        )
+
+    if minimal_processing:
+        thermr = unresr = purr = False
+
+    deck = _build_deck(
+        mat=mat, ek=ek, temperature=float(temperature),
+        err=err, thermr=thermr, unresr=unresr, purr=purr,
+        iwt=iwt, spectrum=spectrum, relative=relative,
+        irespr=irespr, mt=mt, iprint=iprint,
+    )
+    print("Running NJOY to produce MF33 covariance data")
+    tapes = _run_njoy(deck, endf_path, exe=njoy_exec)
+
+    return {"mat": int(mat), "ek": ek, "tape33": tapes[33]}

openmc/data/neutron.py

@@ -109,6 +109,7 @@ class IncidentNeutron(EqualityMixin):
         self.reactions = {}
         self._urr = {}
         self._resonances = None
+        self._mg_covariance = None
 
     def __contains__(self, mt):
         return mt in self.reactions
@@ -229,6 +230,30 @@ class IncidentNeutron(EqualityMixin):
             cv.check_type('probability tables', value, ProbabilityTables)
         self._urr = urr
 
+    @property
+    def mg_covariance(self):
+        return self._mg_covariance
+
+    @mg_covariance.setter
+    def mg_covariance(self, cov):
+        if cov is None:
+            self._mg_covariance = None
+            return
+        if not hasattr(cov, 'write_mf33_group'):
+            raise TypeError(
+                "mg_covariance must provide write_mf33_group "
+                "(e.g. NeutronXSCovariances)")
+        self._mg_covariance = cov
+
+    # Optional alias if you already use `data.covariance = ...` somewhere
+    @property
+    def covariance(self):
+        return self._mg_covariance
+    
+    @covariance.setter
+    def covariance(self, cov):
+        self.mg_covariance = cov
+
     @property
     def temperatures(self):
         return [f"{int(round(kT / K_BOLTZMANN))}K" for kT in self.kTs]
@@ -429,7 +454,16 @@ class IncidentNeutron(EqualityMixin):
             if self.fission_energy is not None:
                 fer_group = g.create_group('fission_energy_release')
                 self.fission_energy.to_hdf5(fer_group)
+            
+            # Write covariance data (per-reaction MF=33 schema)
+            if self._mg_covariance is not None:
+                cov = self._mg_covariance
+                cov_root = g.require_group("covariance")
 
+                # Prefer the shared writer if available (NeutronXSCovariances)
+                if hasattr(cov, "write_mf33_group"):
+                    cov.write_mf33_group(cov_root)
+        
     @classmethod
     def from_hdf5(cls, group_or_filename):
         """Generate continuous-energy neutron interaction data from HDF5 group
@@ -504,7 +538,13 @@ class IncidentNeutron(EqualityMixin):
         if 'fission_energy_release' in group:
             fer_group = group['fission_energy_release']
             data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
-
+        
+        # Read covariance data (per-reaction MF=33 schema)
+        if 'covariance' in group and 'mf33' in group['covariance']:
+            from .xs_covariance_njoy import NeutronXSCovariances
+            data._mg_covariance = NeutronXSCovariances._read_mf33_group(
+                group['covariance']['mf33'], name=name,)
+            
         return data
 
     @classmethod

openmc/data/xs_covariance_njoy.py

@@ -0,0 +1,873 @@
+"""openmc.data.xs_covariance
+
+Tiny model + HDF5 I/O for multigroup **cross section** covariances (ERRORR MF=33).
+
+This is intentionally narrow:
+- neutron MF=33 only
+- explicit energy grid (group edges in eV)
+- **relative** covariances (what ERRORR typically produces with IRELCO=1)
+
+Public API
+----------
+NeutronXSCovariances.from_endf(...)
+NeutronXSCovariances.to_hdf5(...)
+NeutronXSCovariances.from_hdf5(...)
+
+Everything else is internal parsing.
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Sequence, Tuple
+
+import numpy as np
+import scipy.linalg as la
+from .mf33_njoy import generate_errorr_mf33
+
+log = logging.getLogger(__name__)
+
+# -----------------------------------------------------------------------------
+# ERRORR tape33 (ENDF-6 text) parser for MF=33
+# -----------------------------------------------------------------------------
+
+from collections import namedtuple
+
+_CONT = namedtuple("CONT", "C1 C2 L1 L2 N1 N2")
+_LIST = namedtuple("LIST", "C1 C2 L1 L2 NPL N2 B")
+
+
+def _endf_int(field: str) -> int:
+    field = field.strip()
+    if not field:
+        return 0
+    try:
+        return int(field)
+    except ValueError:
+        return int(float(field))
+
+
+def _endf_float(field: str) -> float:
+    s = field.strip()
+    if not s:
+        return 0.0
+    if "e" in s.lower():
+        return float(s)
+    # ENDF-style: mantissa + exponent sign+digits, no 'e'
+    for i in range(len(s) - 1, 0, -1):
+        if s[i] in "+-" and s[i - 1].isdigit():
+            mant = s[:i]
+            exp = s[i:]
+            return float(f"{mant}e{exp}")
+    return float(s)
+
+
+def _parse_endf_ids(line: str) -> Tuple[int, int, int]:
+    if len(line) < 75:
+        return (0, 0, 0)
+    return (int(line[66:70]), int(line[70:72]), int(line[72:75]))
+
+
+def _read_cont_line(line: str) -> _CONT:
+    return _CONT(
+        _endf_float(line[0:11]),
+        _endf_float(line[11:22]),
+        _endf_int(line[22:33]),
+        _endf_int(line[33:44]),
+        _endf_int(line[44:55]),
+        _endf_int(line[55:66]),
+    )
+
+
+def _read_list_record(lines: List[str], ipos: int) -> Tuple[_LIST, int]:
+    C = _read_cont_line(lines[ipos])
+    npl = int(C.N1)
+    ipos += 1
+    vals: List[float] = []
+    nlines = int(np.ceil(npl / 6.0))
+    for _ in range(nlines):
+        ln = lines[ipos]
+        for j in range(6):
+            a = 11 * j
+            b = a + 11
+            if a >= 66:
+                break
+            vals.append(_endf_float(ln[a:b]))
+        ipos += 1
+    vals = vals[:npl]
+    return _LIST(C.C1, C.C2, C.L1, C.L2, npl, C.N2, vals), ipos
+
+
+def _parse_mf33_section_lines(section_lines: List[str], mat: int, mt: int) -> Dict[str, Any]:
+    i = 0
+    C0 = _read_cont_line(section_lines[i])
+    i += 1
+    out: Dict[str, Any] = {"MAT": mat, "MF": 33, "MT": mt, "ZA": C0.C1, "AWR": C0.C2}
+
+    reaction_pairs: Dict[int, np.ndarray] = {}
+    for _ in range(int(C0.N2)):  # number of reaction pairs
+        C = _read_cont_line(section_lines[i])
+        i += 1
+        mt1 = int(C.L2)
+        ng = int(C.N2)
+        M = np.zeros((ng, ng))
+        while True:
+            L, i = _read_list_record(section_lines, i)
+            ngcol = int(L.L1)
+            grow = int(L.N2)
+            gcol = int(L.L2)
+            # LIST rows are 1-based
+            M[grow - 1, gcol - 1 : gcol - 1 + ngcol] = np.asarray(L.B, dtype=float)
+            if grow >= ng:
+                break
+        reaction_pairs[mt1] = M
+
+    out["COVS"] = reaction_pairs
+    return out
+
+
+def parse_errorr_mf33_text(tape33_text: str, mat: Optional[int] = None) -> Dict[int, Dict[str, Any]]:
+    lines = tape33_text.splitlines(True)
+    reactions: Dict[int, Dict[str, Any]] = {}
+
+    current: List[str] = []
+    current_key: Optional[Tuple[int, int, int]] = None
+
+    def flush():
+        nonlocal current, current_key
+        if current_key is None or not current:
+            current = []
+            current_key = None
+            return
+        m, mf, mt = current_key
+        if mf == 33 and mt > 0 and (mat is None or m == mat):
+            reactions[mt] = _parse_mf33_section_lines(current, m, mt)
+        current = []
+        current_key = None
+
+    for ln in lines:
+        m, mf, mt = _parse_endf_ids(ln)
+
+        if m == 0 and mf == 0 and mt == 0:
+            flush()
+            continue
+        if mt == 0:
+            flush()
+            continue
+
+        key = (m, mf, mt)
+        if current_key is None:
+            current_key = key
+        elif key != current_key:
+            flush()
+            current_key = key
+
+        current.append(ln)
+
+    flush()
+    return reactions
+
+# -----------------------------------------------------------------------------
+# Validating raw covariance data from NJOY
+# -----------------------------------------------------------------------------
+@dataclass
+class RawBlockDiagnostics:
+    """Stage-1 QA diagnostics for one raw MF=33 block."""
+    mt: int
+    mt1: int
+    is_self_covariance: bool
+    status: str = "pass"   # "pass", "warn", "fail"
+    messages: List[str] = field(default_factory=list)
+
+    # --- eigenvalue spectrum diagnostics (self-blocks only) ---
+    n_negative_eigenvalues: Optional[int] = None
+    negativity_ratio: Optional[float] = None      # |λ_min_neg| / λ_max_pos
+    negative_spectral_fraction: Optional[float] = None  # Σ|λ_neg| / Σ|λ_all|
+
+    # --- correlation matrix bounds (self-blocks only) ---
+    max_abs_correlation: Optional[float] = None
+
+    # --- relative uncertainty magnitude (self-blocks only) ---
+    max_rel_std: Optional[float] = None
+    n_groups_rel_std_above_1: Optional[int] = None
+    n_groups_rel_std_above_10: Optional[int] = None
+ 
+    def _flag(self, level: str, msg: str):
+        """Set status to *level* (unless already 'fail') and append *msg*."""
+        if level == "fail" or self.status != "fail":
+            self.status = level
+        self.messages.append(msg)
+ 
+ 
+def _rel_frob(A: np.ndarray, B: np.ndarray) -> float:
+    """||A-B||_F / max(||A||_F, ||B||_F, 1)."""
+    nA, nB = np.linalg.norm(A, "fro"), np.linalg.norm(B, "fro")
+    return float(np.linalg.norm(A - B, "fro") / max(nA, nB, 1.0))
+ 
+ 
+def validate_raw_mf33_reactions(
+    reactions: Dict[int, Dict[str, Any]],
+    energy_grid_ev: Sequence[float],
+    *,
+    atol: float = 1e-14,
+    rtol: float = 1e-10,
+    warn_only_missing_partner: bool = True,
+    negativity_ratio_warn: float = 1e-3,
+    negativity_ratio_fail: float = 1e-1,
+    rel_std_warn: float = 1.0,
+    rel_std_fail: float = 10.0,
+) -> Dict[int, Dict[int, RawBlockDiagnostics]]:
+    """Stage-1 integrity checks on raw MF=33 matrices from ERRORR tape33.
+ 
+    Checks per block: shape == (G,G), finiteness, and—for self blocks—
+    symmetry, non-negative diagonal, eigenvalue spectrum, correlation
+    matrix bounds, and relative uncertainty magnitude; for cross blocks,
+    transpose-partner consistency C(mt,mt1) ≈ C(mt1,mt)^T.
+
+    Parameters
+    ----------
+    negativity_ratio_warn, negativity_ratio_fail
+        Thresholds on |λ_min_neg| / λ_max_pos for the eigenvalue spectrum
+        check.  Following the convention that 1e-3 is a warning and 1e-1
+        signals the approximate (repaired) matrix will differ substantially.
+    rel_std_warn, rel_std_fail
+        Thresholds on group-wise relative standard deviations
+        (sqrt of diagonal, since IRELCO=1).  Values above rel_std_warn
+        (default 1.0 = 100%) are flagged as warnings; above rel_std_fail
+        (default 10.0 = 1000%) as likely processing artifacts.
+ 
+    Returns nested dict ``diagnostics[mt][mt1] = RawBlockDiagnostics``.
+    """
+    G = len(energy_grid_ev) - 1
+    shape_exp = (G, G)
+    diags: Dict[int, Dict[int, RawBlockDiagnostics]] = {}
+ 
+    for mt, sec in reactions.items():
+        diags[mt] = {}
+        for mt1, M in sec.get("COVS", {}).items():
+            A = np.asarray(M, dtype=np.float64)
+            is_self = int(mt) == int(mt1)
+            dx = RawBlockDiagnostics(int(mt), int(mt1), is_self)
+ 
+            # 1) Shape / finiteness (hard failures) -------------------------
+            if A.shape != shape_exp:
+                dx._flag("fail", f"Wrong shape {A.shape}, expected {shape_exp}.")
+            if not np.all(np.isfinite(A)):
+                dx._flag("fail", "Matrix contains NaN or Inf values.")
+            if dx.status == "fail":
+                diags[mt][mt1] = dx
+                continue
+ 
+            # 2) Self-covariance: symmetry + diagonal -----------------------
+            if is_self:
+                res = _rel_frob(A, A.T)
+                if not np.allclose(A, A.T, atol=atol, rtol=rtol):
+                    dx._flag("fail", f"Not symmetric (residual={res:.3e}).")
+                d = np.diag(A)
+                n_neg = int(np.count_nonzero(d < -atol))
+                n_zero = int(np.count_nonzero(np.isclose(d, 0.0, atol=atol)))
+                if n_neg > 0:
+                    dx._flag("fail", f"{n_neg} negative diagonal entries; "
+                             f"min(diag)={float(d.min()):.6e}.")
+                elif n_zero > 0:
+                    dx._flag("warn", f"{n_zero} zero-variance diagonal groups.")
+
+                # 2a) Eigenvalue spectrum diagnostics -----------------------
+                #     Symmetrise before computing eigenvalues so that eigh
+                #     is applicable even when ERRORR left tiny asymmetries.
+                A_sym = 0.5 * (A + A.T)
+                eigenvalues = la.eigvalsh(A_sym)
+
+                neg_mask = eigenvalues < -atol
+                pos_mask = eigenvalues > atol
+                n_neg_eig = int(neg_mask.sum())
+                dx.n_negative_eigenvalues = n_neg_eig
+
+                lam_max_pos = float(eigenvalues[pos_mask].max()) if pos_mask.any() else 0.0
+                abs_neg = np.abs(eigenvalues[neg_mask]) if neg_mask.any() else np.empty(0)
+                lam_min_neg_abs = float(abs_neg.max()) if abs_neg.size > 0 else 0.0
+
+                if lam_max_pos > 0.0 and lam_min_neg_abs > 0.0:
+                    dx.negativity_ratio = lam_min_neg_abs / lam_max_pos
+                else:
+                    dx.negativity_ratio = 0.0
+
+                total_spectral_weight = float(np.abs(eigenvalues).sum())
+                if total_spectral_weight > 0.0:
+                    dx.negative_spectral_fraction = float(abs_neg.sum()) / total_spectral_weight
+                else:
+                    dx.negative_spectral_fraction = 0.0
+
+                if n_neg_eig > 0:
+                    if dx.negativity_ratio >= negativity_ratio_fail:
+                        dx._flag(
+                            "fail",
+                            f"Eigenvalue spectrum: {n_neg_eig} negative eigenvalue(s), "
+                            f"negativity ratio={dx.negativity_ratio:.3e} "
+                            f"(>= {negativity_ratio_fail:.0e} threshold); "
+                            f"repaired matrix will differ substantially from original."
+                        )
+                    elif dx.negativity_ratio >= negativity_ratio_warn:
+                        dx._flag(
+                            "warn",
+                            f"Eigenvalue spectrum: {n_neg_eig} negative eigenvalue(s), "
+                            f"negativity ratio={dx.negativity_ratio:.3e} "
+                            f"(>= {negativity_ratio_warn:.0e} threshold)."
+                        )
+                    else:
+                        # Small negative eigenvalues — informational only
+                        dx.messages.append(
+                            f"Eigenvalue spectrum: {n_neg_eig} small negative "
+                            f"eigenvalue(s), negativity ratio={dx.negativity_ratio:.3e}."
+                        )
+
+                # 2b) Correlation matrix bounds check -----------------------
+                #     Convert self-covariance to correlation and verify |ρ|≤1.
+                d_safe = d.copy()
+                d_safe[d_safe <= 0.0] = np.inf  # skip zero-variance groups
+                inv_std = 1.0 / np.sqrt(d_safe)
+                corr = A_sym * np.outer(inv_std, inv_std)
+                # Zero out rows/cols that had zero variance (they are undefined)
+                zero_var_mask = d <= 0.0
+                corr[zero_var_mask, :] = 0.0
+                corr[:, zero_var_mask] = 0.0
+                np.fill_diagonal(corr, 1.0)  # diagonal is ρ=1 by definition
+
+                max_abs_rho = float(np.max(np.abs(
+                    corr[np.triu_indices_from(corr, k=1)]
+                ))) if G > 1 else 0.0
+                dx.max_abs_correlation = max_abs_rho
+
+                if max_abs_rho > 1.0 + rtol:
+                    n_violating = int(np.count_nonzero(
+                        np.abs(corr[np.triu_indices_from(corr, k=1)]) > 1.0 + rtol
+                    ))
+                    dx._flag(
+                        "warn",
+                        f"Correlation matrix has {n_violating} off-diagonal "
+                        f"element(s) with |rho| > 1 (max |rho|={max_abs_rho:.6f}); "
+                        f"unphysical covariance entries."
+                    )
+
+                # 2c) Relative uncertainty magnitude check ------------------
+                #     For IRELCO=1, sqrt(diag) gives group-wise relative sigma.
+                rel_std = np.sqrt(np.maximum(d, 0.0))
+                dx.max_rel_std = float(rel_std.max()) if rel_std.size > 0 else 0.0
+                dx.n_groups_rel_std_above_1 = int(np.count_nonzero(rel_std > rel_std_warn))
+                dx.n_groups_rel_std_above_10 = int(np.count_nonzero(rel_std > rel_std_fail))
+
+                if dx.n_groups_rel_std_above_10 > 0:
+                    dx._flag(
+                        "warn",
+                        f"Relative uncertainty: {dx.n_groups_rel_std_above_10} group(s) "
+                        f"with sigma_rel > {rel_std_fail} "
+                        f"(max={dx.max_rel_std:.3f}); "
+                        f"likely ERRORR processing artifact."
+                    )
+                elif dx.n_groups_rel_std_above_1 > 0:
+                    dx._flag(
+                        "warn",
+                        f"Relative uncertainty: {dx.n_groups_rel_std_above_1} group(s) "
+                        f"with sigma_rel > {rel_std_warn} "
+                        f"(max={dx.max_rel_std:.3f})."
+                    )
+ 
+            # 3) Cross-covariance: transpose partner ------------------------
+            else:
+                B_raw = reactions.get(int(mt1), {}).get("COVS", {}).get(int(mt))
+                if B_raw is None:
+                    pass  # Expected: ERRORR stores only the upper-triangle block
+                else:
+                    B = np.asarray(B_raw, dtype=np.float64)
+                    if B.shape != shape_exp:
+                        dx._flag("fail", f"Partner wrong shape {B.shape}.")
+                    elif not np.all(np.isfinite(B)):
+                        dx._flag("fail", "Partner contains NaN/Inf.")
+                    elif not np.allclose(A, B.T, atol=atol, rtol=rtol):
+                        dx._flag("fail",
+                                f"C({mt},{mt1}) != C({mt1},{mt})^T "
+                                f"(residual={_rel_frob(A, B.T):.3e}).")
+ 
+            diags[mt][mt1] = dx
+ 
+    return diags
+
+# -----------------------------------------------------------------------------
+# Covariance factor computation (eigendecomposition + QR)
+# -----------------------------------------------------------------------------
+
+@dataclass
+class CovFactorResult:
+    """Container for a covariance factorization after diagnostics are done.
+
+    Attributes
+    ----------
+    L : np.ndarray
+        Lower-triangular factor, shape (G, r) where r = effective_rank.
+        Satisfies L @ L.T ≈ original covariance (after any regularization).
+    effective_rank : int
+        Number of positive eigenvalues retained.
+    full_size : int
+        Original matrix dimension G.
+    method : str
+        Which path was taken: "cholesky", "eigen_qr", or "zero_matrix".
+    """
+    L: np.ndarray
+    effective_rank: int
+    full_size: int
+    method: str = "eigen_qr"
+  
+
+def _enforce_symmetry(A: np.ndarray) -> np.ndarray:
+    """Force exact symmetry: A_sym = (A + A^T) / 2."""
+    return 0.5 * (A + A.T)
+
+
+def _strip_zero_variance(A: np.ndarray):
+    """
+    Remove rows/columns with zero diagonal (zero variance).
+    """
+    diag = np.diag(A)
+    nz = np.flatnonzero(diag > 0.0)
+    if len(nz) == 0:
+        return np.empty((0, 0), dtype=A.dtype), nz
+    return A[np.ix_(nz, nz)], nz
+
+def compute_covariance_factor(
+    cov_matrix: np.ndarray,
+    tol: float = 1e-10,
+
+) -> CovFactorResult:
+    """
+    Compute a lower-triangular factor L such that L @ L.T ≈ Sigma.
+    """
+    G = cov_matrix.shape[0]
+
+    # --- 0. Enforce exact symmetry ---
+    A = _enforce_symmetry(np.asarray(cov_matrix, dtype=np.float64))
+
+    # --- 1. Strip zero-variance rows/columns ---
+    Ar, nz = _strip_zero_variance(A)
+
+    if Ar.size == 0:
+        return CovFactorResult(
+            L=np.zeros((G, 0), dtype=np.float64),
+            effective_rank=0,
+            full_size=G,
+            method="zero_matrix",
+        )
+
+    Gr = Ar.shape[0]
+
+    # --- 2a. Fast path: standard Cholesky ---
+    try:
+        Lr = la.cholesky(Ar, lower=True)
+
+        L_full = np.zeros((G, Gr), dtype=np.float64)
+        L_full[nz, :] = Lr
+
+        return CovFactorResult(
+            L=L_full,
+            effective_rank=Gr,
+            full_size=G,
+            method="cholesky"
+        )
+    except la.LinAlgError:
+        pass
+
+    # --- 2b. Fallback: eigendecomposition + QR ---
+    eigenvalues, V = la.eigh(Ar)
+
+    max_eig = eigenvalues.max() if eigenvalues.max() > 0 else 1.0
+    eigenvalues[eigenvalues / max_eig < tol] = 0.0
+    eigenvalues[eigenvalues < 0.0] = 0.0
+
+    pos_mask = eigenvalues > 0.0
+    r = int(pos_mask.sum())
+
+    if r == 0:
+        return CovFactorResult(
+            L=np.zeros((G, 0), dtype=np.float64),
+            effective_rank=0,
+            full_size=G,
+            method="eigen_qr",
+        )
+
+    V_pos = V[:, pos_mask]
+    sqrt_lam = np.sqrt(eigenvalues[pos_mask])
+    Athin = V_pos * sqrt_lam[np.newaxis, :]
+
+    _, R = la.qr(Athin.T, mode="economic")
+    Lr = R.T
+
+    L_full = np.zeros((G, r), dtype=np.float64)
+    L_full[nz, :] = Lr
+
+    return CovFactorResult(
+        L=L_full,
+        effective_rank=r,
+        full_size=G,
+        method="eigen_qr",
+    )
+
+# -----------------------------------------------------------------------------
+# Data model + HDF5 I/O
+# -----------------------------------------------------------------------------
+
+@dataclass
+class NeutronXSCovariances:
+    """Neutron multigroup **relative** covariance data for cross sections (MF=33)."""
+
+    name: str
+    energy_grid_ev: np.ndarray              # shape (G+1,)
+    reactions: Dict[int, Dict[str, Any]]    # MT -> parsed dict from ERRORR
+    mat: Optional[int] = None
+    temperature_k: Optional[float] = None   # what we processed at (single T)
+    factor_results: Optional[Dict[int, Dict[int, CovFactorResult]]] = None 
+
+    @classmethod
+    def from_endf(
+        cls,
+        endf_path: str | Path,
+        energy_grid_ev: Sequence[float],
+        *,
+        njoy_exec: Optional[str] = None,
+        mat: Optional[int] = None,
+        temperature: float = 293.6,
+        name: Optional[str] = None,
+        compute_factors: bool = True,
+        eig_tol: float = 1e-10,
+    ) -> "NeutronXSCovariances":
+        res = generate_errorr_mf33(
+            endf_path,
+            energy_grid_ev,
+            njoy_exec=njoy_exec,
+            mat=mat,
+            temperature=temperature,
+        )
+        mat_used = int(res["mat"])
+        ek = np.asarray(res["ek"], dtype=float)
+
+        reactions = parse_errorr_mf33_text(res["tape33"], mat=mat_used)
+
+        # --- Parsed block inventory ---
+        n_self, n_cross = 0, 0
+        cross_pairs: List[Tuple[int, int]] = []
+        for mt_key, sec in reactions.items():
+            for mt1_key in sec.get("COVS", {}):
+                if int(mt_key) == int(mt1_key):
+                    n_self += 1
+                else:
+                    n_cross += 1
+                    cross_pairs.append((int(mt_key), int(mt1_key)))
+        log.info(
+            "ERRORR tape33 parsed: %d MT section(s), "
+            "%d self-covariance block(s), %d cross-covariance block(s)",
+            len(reactions), n_self, n_cross,
+        )
+        if cross_pairs:
+            for mt_a, mt1_a in cross_pairs:
+                log.info("  cross-covariance: MT %d <-> MT %d", mt_a, mt1_a)
+        else:
+            log.info(
+                "  No explicit cross-covariance blocks in ERRORR output.  "
+                "Cross-reaction correlations (if any) are handled implicitly "
+                "via NC-type derivation relations in the evaluation."
+            )
+
+        raw_validation = validate_raw_mf33_reactions(reactions, ek, atol=1e-14, rtol=1e-10)
+        
+        n_warnings = 0
+        for mt, blocks in raw_validation.items():
+            for mt1, qa in blocks.items():
+                if qa.status == "fail":
+                    raise ValueError(
+                        f"Stage-1 MF=33 validation failed for MT {mt} -> MT1 {mt1}: "
+                        + "; ".join(qa.messages)
+                    )
+                elif qa.status == "warn":
+                    n_warnings += 1
+                    log.info(
+                        "Stage-1 MF=33 warning for MT %s -> MT1 %s: %s",
+                        mt, mt1, "; ".join(qa.messages)
+                    )
+        # Pre-compute covariance factors for all sub-blocks
+        n_eigen_qr = 0
+        fcache: Optional[Dict[int, Dict[int, CovFactorResult]]] = None
+        if compute_factors:
+            fcache = {}
+            for mt_key, sec in reactions.items():
+                fcache[mt_key] = {}
+                for mt1, M in sec.get("COVS", {}).items():
+                    if int(mt_key) != int(mt1):
+                        log.info(
+                            "  MT %s -> MT1 %s: cross-block (%d x %d), "
+                            "stored raw (no factorization)",
+                            mt_key, mt1,
+                            M.shape[0], M.shape[1] if M.ndim > 1 else 0,
+                        )
+                        continue  # cross-block: store raw only
+                    result = compute_covariance_factor(
+                        np.asarray(M, dtype=np.float64),
+                        tol=eig_tol,
+                    )
+                    fcache[mt_key][mt1] = result
+                    if result.method == "eigen_qr":
+                        n_eigen_qr += 1
+                    log.info(
+                        "  MT %s -> MT1 %s: method=%-9s  rank=%d/%d  ",
+                        mt_key, mt1,
+                        result.method,
+                        result.effective_rank,
+                        result.full_size
+                    )
+
+        # --- One-line summary ---
+        log.info(
+            "%s: %d self-block(s), %d cross-block(s), "
+            "%d eigen_qr fallback(s), %d warning(s)",
+            name if name is not None else Path(endf_path).stem,
+            n_self, n_cross, n_eigen_qr, n_warnings,
+        )
+
+        if name is None:
+            name = Path(endf_path).stem
+
+        return cls(
+            name=str(name),
+            energy_grid_ev=ek,
+            reactions=reactions,
+            mat=mat_used,
+            temperature_k=float(temperature),
+            factor_results=fcache,
+        )
+    
+    # -----------------------------------------------------------------
+    # HDF5 writing
+    # -----------------------------------------------------------------
+
+    def to_hdf5(self, filename: str | Path, store_raw_covariance: bool = True) -> None:
+        """
+        Write covariances to a standalone HDF5 file.
+        """
+        import h5py
+
+        filename = Path(filename)
+        filename.parent.mkdir(parents=True, exist_ok=True)
+
+        with h5py.File(filename, "w") as f:
+            self.write_mf33_group(f, store_raw_covariance=store_raw_covariance)
+
+    def write_mf33_group(self, h5_group, store_raw_covariance: bool = True,
+                         eig_tol: float = 1e-10) -> None:
+        """
+        Write the "mf33" sub-tree into an already-open HDF5 group.
+        """
+        if "mf33" in h5_group:
+            del h5_group["mf33"]
+        mf33 = h5_group.create_group("mf33")
+
+        mf33.attrs["format"] = np.bytes_("openmc.mf33.v1")
+        mf33.attrs["source"] = np.bytes_("njoy errorr")
+        mf33.attrs["relative"] = 1  # int flag – portable across languages
+        mf33.attrs["store_raw_covariance"] = int(store_raw_covariance)
+        #mf33.attrs["factorization"] = np.bytes_("eigen_qr_thin")
+        if self.mat is not None:
+            mf33.attrs["mat"] = int(self.mat)
+        if self.temperature_k is not None:
+            mf33.attrs["temperature_k"] = float(self.temperature_k)
+
+        mf33.create_dataset(
+            "energy_grid_ev",
+            data=np.asarray(self.energy_grid_ev, dtype=np.float64),
+        )
+        nonempty_mts = sorted(mt for mt, sec in self.reactions.items() if sec.get("COVS", {}))
+        
+        mf33.create_dataset(
+            "mts",
+            data=np.asarray(nonempty_mts, dtype=np.int32),
+        )
+
+        if store_raw_covariance:
+            greact = mf33.create_group("reactions")
+        gfact = mf33.create_group("factors")
+
+        cached = self.factor_results or {}
+
+        for mt, sec in self.reactions.items():
+            if not sec.get("COVS", {}):
+                continue 
+            mt_str = str(int(mt))
+
+            if store_raw_covariance:
+                gmt = greact.create_group(mt_str)
+                for attr_name in ("ZA", "AWR"):
+                    if attr_name in sec:
+                        gmt.attrs[attr_name] = float(sec[attr_name])
+
+            covs: Dict[int, np.ndarray] = sec.get("COVS", {})
+            for mt1, M in covs.items():
+                M_arr = np.asarray(M, dtype=np.float64)
+                ds_name = str(int(mt1))
+                is_self = int(mt) == int(mt1)
+
+                # ---- raw covariance ----
+                # Self-blocks: optional. Cross-blocks: always stored (no factor).
+                if store_raw_covariance or not is_self:
+                    if not store_raw_covariance:
+                        # Ensure the reactions group exists for cross-blocks
+                        if "reactions" not in mf33:
+                            greact = mf33.create_group("reactions")
+                        if mt_str not in greact:
+                            gmt = greact.create_group(mt_str)
+                            for attr_name in ("ZA", "AWR"):
+                                if attr_name in sec:
+                                    gmt.attrs[attr_name] = float(sec[attr_name])
+                        else:
+                            gmt = greact[mt_str]
+                    gmt.create_dataset(
+                        ds_name,
+                        data=M_arr,
+                        compression="gzip",
+                        shuffle=True,
+                    )
+
+                # ---- Triangular factor (self-blocks only) ----
+                if is_self:
+                    result = cached.get(mt, {}).get(mt1, None)
+                    if result is None:
+                        result = compute_covariance_factor(
+                            M_arr, tol=eig_tol,)
+
+                    if mt_str not in gfact:
+                        gmt_fact = gfact.create_group(mt_str)
+                        for attr_name in ("ZA", "AWR"):
+                            if attr_name in sec:
+                                gmt_fact.attrs[attr_name] = float(sec[attr_name])
+                    else:
+                        gmt_fact = gfact[mt_str]
+
+                    gmt_fact.create_dataset(
+                        ds_name,
+                        data=result.L,
+                        compression="gzip",
+                        shuffle=True,
+                    )
+
+    # -----------------------------------------------------------------
+    # HDF5 reading
+    # -----------------------------------------------------------------
+ 
+    @classmethod
+    def _read_mf33_group(cls, mf33_group, name: str = "unknown") -> "NeutronXSCovariances":
+        """Reconstruct a :class:`NeutronXSCovariances` from an open HDF5
+        ``mf33`` group (the inverse of :meth:`write_mf33_group`).
+ 
+        Parameters
+        ----------
+        mf33_group : h5py.Group
+            The ``mf33`` group inside the HDF5 file.
+        name : str
+            Nuclide name to attach to the returned object.
+ 
+        Returns
+        -------
+        NeutronXSCovariances
+        """
+        energy_grid_ev = mf33_group["energy_grid_ev"][()]
+        mt_list = mf33_group["mts"][()]
+ 
+        mat = int(mf33_group.attrs["mat"]) if "mat" in mf33_group.attrs else None
+        temperature_k = (
+            float(mf33_group.attrs["temperature_k"])
+            if "temperature_k" in mf33_group.attrs
+            else None
+        )
+ 
+        has_raw = "reactions" in mf33_group
+        has_factors = "factors" in mf33_group
+ 
+        reactions: Dict[int, Dict[str, Any]] = {}
+        factor_results: Dict[int, Dict[int, CovFactorResult]] = {}
+ 
+        for mt in mt_list:
+            mt_str = str(int(mt))
+            sec: Dict[str, Any] = {"MAT": mat, "MF": 33, "MT": int(mt)}
+ 
+            # --- Read raw covariance matrices (optional) ---
+            covs: Dict[int, np.ndarray] = {}
+            if has_raw and mt_str in mf33_group["reactions"]:
+                gmt = mf33_group["reactions"][mt_str]
+                for attr_name in ("ZA", "AWR"):
+                    if attr_name in gmt.attrs:
+                        sec[attr_name] = float(gmt.attrs[attr_name])
+                for ds_name in gmt:
+                    mt1 = int(ds_name)
+                    covs[mt1] = gmt[ds_name][()]
+            sec["COVS"] = covs
+ 
+            # --- Read triangular factors ---
+            mt_factors: Dict[int, CovFactorResult] = {}
+            if has_factors and mt_str in mf33_group["factors"]:
+                gmt_fact = mf33_group["factors"][mt_str]
+                for attr_name in ("ZA", "AWR"):
+                    if attr_name in gmt_fact.attrs:
+                        sec.setdefault(attr_name, float(gmt_fact.attrs[attr_name]))
+                for ds_name in gmt_fact:
+                    mt1 = int(ds_name)
+                    L = gmt_fact[ds_name][()]
+                    G = L.shape[0]
+                    r = L.shape[1] if L.ndim == 2 else 0
+                    mt_factors[mt1] = CovFactorResult(
+                        L=L,
+                        effective_rank=r,
+                        full_size=G,
+                        method="from_hdf5",
+                    )
+ 
+                    # If raw covariance was not stored, synthesise a
+                    # minimal COVS entry so downstream code that iterates
+                    # over reactions["COVS"] still works.
+                    if mt1 not in covs:
+                        covs[mt1] = L @ L.T
+ 
+            reactions[int(mt)] = sec
+            factor_results[int(mt)] = mt_factors
+ 
+        return cls(
+            name=name,
+            energy_grid_ev=energy_grid_ev,
+            reactions=reactions,
+            mat=mat,
+            temperature_k=temperature_k,
+            factor_results=factor_results if factor_results else None,
+        )
+ 
+    @classmethod
+    def from_hdf5(cls, filename: "str | Path") -> "NeutronXSCovariances":
+        """Read covariances from a standalone HDF5 file written by
+        :meth:`to_hdf5`.
+ 
+        Parameters
+        ----------
+        filename : str or Path
+            Path to the HDF5 file.
+ 
+        Returns
+        -------
+        NeutronXSCovariances
+        """
+        import h5py
+ 
+        filename = Path(filename)
+        with h5py.File(filename, "r") as f:
+            if "mf33" in f:
+                mf33 = f["mf33"]
+            else:
+                raise KeyError(
+                    f"HDF5 file '{filename}' does not contain an 'mf33' group."
+                )
+            name = filename.stem
+            return cls._read_mf33_group(mf33, name=name)

tests/unit_tests/test_xs_covariance.py

@@ -0,0 +1,578 @@
+"""Unit tests for openmc.data.xs_covariance_njoy
+
+These tests use synthetic data only (no NJOY executable or ENDF files
+are required).  They test compute_covariance_factor on positive-definite, semi-definite, 
+and zero matrices. The HDF5 round-trip for NeutronXSCovariances (with and without raw 
+covariance storage). Finally, its testing the integration with export_to_hdf5
+ and from_hdf5 via the mg_covariance property.
+"""
+
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+
+from openmc.data.xs_covariance_njoy import (
+    CovFactorResult,
+    NeutronXSCovariances,
+    RawBlockDiagnostics,
+    compute_covariance_factor,
+    validate_raw_mf33_reactions,
+)
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+ 
+# Default group count and energy grid used by validate_raw_mf33_reactions tests
+G = 4
+EK = np.logspace(-5, 7, G + 1).tolist()  # G+1 boundaries in eV
+ 
+ 
+def _posdef(n, seed=42):
+    """Return a random (n x n) symmetric positive-definite matrix."""
+    rng = np.random.default_rng(seed)
+    A = rng.standard_normal((n, n))
+    return A @ A.T + 0.1 * np.eye(n)
+ 
+ 
+def _semidef(n, rank, seed=42):
+    """Return a random (n x n) symmetric PSD matrix with given rank."""
+    rng = np.random.default_rng(seed)
+    B = rng.standard_normal((n, rank))
+    return B @ B.T
+ 
+ 
+def _make_mock_covariances(G=10, seed=42):
+    """Build a NeutronXSCovariances object with synthetic data."""
+    energy = np.logspace(-5, 7, G + 1)
+    cov_mat = _posdef(G, seed=seed)
+    result = compute_covariance_factor(cov_mat)
+ 
+    reactions = {
+        2: {
+            "MAT": 2631, "MF": 33, "MT": 2,
+            "ZA": 26056.0, "AWR": 55.45,
+            "COVS": {2: cov_mat},
+        },
+    }
+    factors = {2: {2: result}}
+ 
+    return NeutronXSCovariances(
+        name="Fe56",
+        energy_grid_ev=energy,
+        reactions=reactions,
+        mat=2631,
+        temperature_k=293.6,
+        factor_results=factors,
+    )
+ 
+ 
+def _wrap_self(mt, M):
+    """Wrap a single self-block into the reactions dict structure."""
+    return {mt: {"COVS": {mt: M}}}
+ 
+ 
+def _wrap_pair(mt, mt1, M_self, M_cross,
+               M_self1=None, M_partner=None):
+    """Wrap a (mt, mt1) pair with cross-block(s)."""
+    reactions = {
+        mt:  {"COVS": {mt: M_self, mt1: M_cross}},
+        mt1: {"COVS": {}},
+    }
+    if M_self1 is not None:
+        reactions[mt1]["COVS"][mt1] = M_self1
+    if M_partner is not None:
+        reactions[mt1]["COVS"][mt] = M_partner
+    return reactions
+ 
+ 
+# ===========================================================================
+# compute_covariance_factor
+# ===========================================================================
+ 
+def test_factor_posdef_cholesky():
+    """Positive-definite input should use the Cholesky path."""
+    Sigma = _posdef(5)
+    result = compute_covariance_factor(Sigma)
+ 
+    assert result.method == "cholesky"
+    assert result.effective_rank == 5
+    assert result.full_size == 5
+    assert result.L.shape == (5, 5)
+    np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-10)
+ 
+ 
+def test_factor_semidef_eigen_qr():
+    """Rank-deficient input should fall back to eigen_qr."""
+    Sigma = _semidef(8, rank=3)
+    result = compute_covariance_factor(Sigma)
+ 
+    assert result.method == "eigen_qr"
+    assert result.effective_rank == 3
+    assert result.full_size == 8
+    assert result.L.shape == (8, 3)
+    np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-8)
+ 
+ 
+def test_factor_zero_matrix():
+    """All-zero input should return rank-0 with the zero_matrix method."""
+    result = compute_covariance_factor(np.zeros((6, 6)))
+ 
+    assert result.method == "zero_matrix"
+    assert result.effective_rank == 0
+    assert result.full_size == 6
+    assert result.L.shape == (6, 0)
+ 
+ 
+def test_factor_partial_zero_variance():
+    """Matrix with some zero-diagonal rows should still factorise the
+    non-zero sub-block correctly."""
+    Sigma = np.zeros((5, 5))
+    small = _posdef(3, seed=99)
+    Sigma[1:4, 1:4] = small
+ 
+    result = compute_covariance_factor(Sigma)
+ 
+    assert result.effective_rank == 3
+    assert result.L.shape[0] == 5
+    np.testing.assert_allclose(result.L @ result.L.T, Sigma, atol=1e-10)
+ 
+ 
+def test_factor_symmetry_enforced():
+    """Slightly asymmetric input should not crash — symmetry is forced."""
+    Sigma = _posdef(4)
+    Sigma[0, 1] += 1e-12
+ 
+    result = compute_covariance_factor(Sigma)
+    assert result.effective_rank == 4
+ 
+# ===========================================================================
+# HDF5 round-trip (standalone file)
+# ===========================================================================
+ 
+def test_hdf5_roundtrip_with_raw_covariance(tmp_path):
+    """Full round-trip storing both raw covariance and L factors."""
+    cov = _make_mock_covariances(G=10)
+    h5_path = tmp_path / "cov_test.h5"
+ 
+    cov.to_hdf5(h5_path, store_raw_covariance=True)
+    cov_read = NeutronXSCovariances.from_hdf5(h5_path)
+ 
+    # Metadata
+    assert cov_read.mat == cov.mat
+    assert cov_read.temperature_k == cov.temperature_k
+    np.testing.assert_allclose(cov_read.energy_grid_ev, cov.energy_grid_ev)
+ 
+    # Reaction keys preserved
+    assert set(cov_read.reactions.keys()) == set(cov.reactions.keys())
+ 
+    # Raw covariance matrices match
+    for mt in cov.reactions:
+        for mt1 in cov.reactions[mt]["COVS"]:
+            original = cov.reactions[mt]["COVS"][mt1]
+            loaded = cov_read.reactions[mt]["COVS"][mt1]
+            np.testing.assert_allclose(loaded, original, atol=1e-12)
+ 
+    # L factors match
+    for mt in cov.factor_results:
+        for mt1 in cov.factor_results[mt]:
+            L_orig = cov.factor_results[mt][mt1].L
+            L_read = cov_read.factor_results[mt][mt1].L
+            np.testing.assert_allclose(L_read, L_orig, atol=1e-12)
+ 
+ 
+def test_hdf5_roundtrip_factors_only(tmp_path):
+    """Round-trip with store_raw_covariance=False.
+ 
+    The raw covariance should be reconstructed as L @ L.T on read.
+    """
+    cov = _make_mock_covariances(G=8)
+    h5_path = tmp_path / "cov_factors_only.h5"
+ 
+    cov.to_hdf5(h5_path, store_raw_covariance=False)
+    cov_read = NeutronXSCovariances.from_hdf5(h5_path)
+ 
+    # Verify "reactions" group is absent in the HDF5
+    with h5py.File(h5_path, "r") as f:
+        assert "reactions" not in f["mf33"]
+ 
+    # COVS should be reconstructed from L factors
+    for mt in cov.reactions:
+        for mt1 in cov.reactions[mt]["COVS"]:
+            original = cov.reactions[mt]["COVS"][mt1]
+            reconstructed = cov_read.reactions[mt]["COVS"][mt1]
+            np.testing.assert_allclose(reconstructed, original, atol=1e-8)
+ 
+ 
+def test_hdf5_schema_attributes(tmp_path):
+    """Verify the expected HDF5 attributes and datasets exist."""
+    cov = _make_mock_covariances(G=5)
+    h5_path = tmp_path / "schema_check.h5"
+    cov.to_hdf5(h5_path)
+ 
+    with h5py.File(h5_path, "r") as f:
+        mf33 = f["mf33"]
+        assert mf33.attrs["format"] == b"openmc.mf33.v1"
+        assert mf33.attrs["relative"] == 1
+        assert mf33.attrs["mat"] == 2631
+        assert "energy_grid_ev" in mf33
+        assert "mts" in mf33
+        assert "factors" in mf33
+ 
+        # Energy grid has G+1 entries
+        assert mf33["energy_grid_ev"].shape == (6,)
+ 
+ 
+def test_hdf5_multiple_reactions(tmp_path):
+    """Round-trip with multiple MT values."""
+    G_local = 6
+    energy = np.logspace(-5, 7, G_local + 1)
+    cov_2 = _posdef(G_local, seed=10)
+    cov_102 = _posdef(G_local, seed=20)
+ 
+    reactions = {
+        2: {"MAT": 2631, "MF": 33, "MT": 2,
+            "ZA": 26056.0, "AWR": 55.45,
+            "COVS": {2: cov_2}},
+        102: {"MAT": 2631, "MF": 33, "MT": 102,
+              "ZA": 26056.0, "AWR": 55.45,
+              "COVS": {102: cov_102}},
+    }
+    cov_obj = NeutronXSCovariances(
+        name="Fe56", energy_grid_ev=energy,
+        reactions=reactions, mat=2631, temperature_k=293.6,
+    )
+    h5_path = tmp_path / "multi_mt.h5"
+    cov_obj.to_hdf5(h5_path)
+    cov_read = NeutronXSCovariances.from_hdf5(h5_path)
+ 
+    assert set(cov_read.reactions.keys()) == {2, 102}
+    np.testing.assert_allclose(
+        cov_read.reactions[2]["COVS"][2], cov_2, atol=1e-12,
+    )
+    np.testing.assert_allclose(
+        cov_read.reactions[102]["COVS"][102], cov_102, atol=1e-12,
+    )
+ 
+ 
+# ===========================================================================
+# write_mf33_group / _read_mf33_group (embedded in another HDF5)
+# ===========================================================================
+ 
+def test_write_read_mf33_into_existing_group(tmp_path):
+    """Write mf33 into a pre-existing HDF5 group, then read it back."""
+    cov = _make_mock_covariances(G=5)
+    h5_path = tmp_path / "embedded.h5"
+ 
+    with h5py.File(h5_path, "w") as f:
+        nuc = f.create_group("Fe56")
+        cov_root = nuc.create_group("covariance")
+        cov.write_mf33_group(cov_root, store_raw_covariance=True)
+ 
+    with h5py.File(h5_path, "r") as f:
+        mf33_group = f["Fe56"]["covariance"]["mf33"]
+        cov_read = NeutronXSCovariances._read_mf33_group(
+            mf33_group, name="Fe56",
+        )
+ 
+    assert cov_read.name == "Fe56"
+    assert cov_read.mat == 2631
+    assert 2 in cov_read.reactions
+ 
+ 
+def test_overwrite_existing_mf33(tmp_path):
+    """Calling write_mf33_group twice should replace, not duplicate."""
+    cov = _make_mock_covariances(G=4)
+    h5_path = tmp_path / "overwrite.h5"
+ 
+    with h5py.File(h5_path, "w") as f:
+        root = f.create_group("root")
+        cov.write_mf33_group(root)
+        cov.write_mf33_group(root)
+ 
+    with h5py.File(h5_path, "r") as f:
+        assert "mf33" in f["root"]
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — hard failures (shape & finiteness)
+# ===========================================================================
+ 
+def test_validate_wrong_shape_fails():
+    M = np.eye(G + 1)  # deliberately wrong
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "fail"
+    assert any("Wrong shape" in m for m in dx.messages)
+    # Downstream diagnostics untouched (None) because the check `continue`s
+    assert dx.n_negative_eigenvalues is None
+    assert dx.max_abs_correlation is None
+ 
+ 
+def test_validate_nan_in_matrix_fails():
+    M = np.eye(G)
+    M[1, 1] = np.nan
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "fail"
+    assert any("NaN or Inf" in m for m in dx.messages)
+    assert dx.n_negative_eigenvalues is None
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — clean self-block baseline
+# ===========================================================================
+ 
+def test_validate_clean_selfblock_passes():
+    # Small variances (~0.01) so rel_std ~ 0.1 → well below warn threshold
+    M = 0.01 * _posdef(G, seed=1)
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "pass"
+    assert dx.mt == 2 and dx.mt1 == 2
+    assert dx.is_self_covariance is True
+    # Diagnostic fields should be populated for self-blocks
+    assert dx.n_negative_eigenvalues is not None
+    assert dx.max_abs_correlation is not None
+    assert dx.max_rel_std is not None
+    assert dx.n_negative_eigenvalues == 0
+    assert dx.max_abs_correlation <= 1.0 + 1e-10
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — symmetry
+# ===========================================================================
+ 
+def test_validate_asymmetric_selfblock_fails():
+    M = _posdef(G, seed=2)
+    M[0, 1] += 0.5  # break symmetry noticeably
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "fail"
+    assert any("Not symmetric" in m for m in dx.messages)
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — diagonal: negative / zero
+# ===========================================================================
+ 
+def test_validate_negative_diagonal_fails():
+    M = np.eye(G) * 0.01
+    M[2, 2] = -0.01
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "fail"
+    assert any("negative diagonal" in m for m in dx.messages)
+ 
+ 
+def test_validate_zero_diagonal_warns():
+    # One group has zero variance; rest are fine.
+    M = np.eye(G) * 0.01
+    M[1, 1] = 0.0
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "warn"
+    assert any("zero-variance" in m for m in dx.messages)
+ 
+ 
+def test_validate_negative_and_zero_diagonal_is_fail_not_warn():
+    """If both negative and zero diagonals are present, fail wins."""
+    M = np.eye(G) * 0.01
+    M[0, 0] = 0.0
+    M[1, 1] = -0.01
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    assert diags[2][2].status == "fail"
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — eigenvalue spectrum
+# ===========================================================================
+ 
+def test_validate_small_negative_eigenvalue_is_informational_only():
+    """neg_ratio well below warn threshold (1e-3) → status stays pass, but
+    an informational message is appended."""
+    rng = np.random.default_rng(5)
+    Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
+    eigs = np.array([1.0, 0.5, 0.2, -1e-6])
+    M = (Q * eigs) @ Q.T
+    M = 0.5 * (M + M.T)  # numerical symmetry
+ 
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    # Diagonal might not all be positive due to random rotation;
+    # check only the eigenvalue-specific behavior.
+    if dx.status == "pass":
+        assert dx.n_negative_eigenvalues >= 1
+        assert dx.negativity_ratio is not None
+        assert dx.negativity_ratio < 1e-3
+        assert any("small negative" in m for m in dx.messages)
+ 
+ 
+def test_validate_moderate_negative_eigenvalue_warns():
+    """neg_ratio between warn (1e-3) and fail (1e-1) → warn."""
+    eigs = np.array([1.0, 0.5, 0.2, -1e-2])  # ratio = 1e-2
+    rng = np.random.default_rng(7)
+    Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
+    M = (Q * eigs) @ Q.T
+    M = 0.5 * (M + M.T)
+ 
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    # Filter out the (possibly flaky) diagonal/correlation side-effects:
+    assert dx.n_negative_eigenvalues == 1
+    assert 1e-3 <= dx.negativity_ratio < 1e-1
+    assert any("negativity ratio" in m and "threshold" in m for m in dx.messages)
+ 
+ 
+def test_validate_large_negative_eigenvalue_fails():
+    """neg_ratio >= fail threshold (1e-1) → fail."""
+    eigs = np.array([1.0, 0.5, 0.2, -0.5])  # ratio = 0.5
+    rng = np.random.default_rng(11)
+    Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
+    M = (Q * eigs) @ Q.T
+    M = 0.5 * (M + M.T)
+ 
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "fail"
+    assert dx.n_negative_eigenvalues == 1
+    assert dx.negativity_ratio >= 1e-1
+    assert any("repaired matrix will differ" in m for m in dx.messages)
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — correlation bound |rho| <= 1
+# ===========================================================================
+ 
+def test_validate_correlation_overshoot_warns():
+    """Construct a symmetric matrix with an off-diagonal that exceeds the
+    geometric-mean bound sqrt(d_ii * d_jj)."""
+    M = np.eye(G) * 1.0
+    M[0, 1] = M[1, 0] = 1.5  # rho = 1.5 > 1
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    # Status is at least warn; could be fail because this also produces a
+    # negative eigenvalue (2x2 block has det < 0).
+    assert dx.status in ("warn", "fail")
+    assert dx.max_abs_correlation > 1.0
+    assert any("Correlation matrix" in m and "rho" in m for m in dx.messages)
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — relative uncertainty magnitude
+# ===========================================================================
+ 
+def test_validate_rel_std_above_warn_warns():
+    """sqrt(diag) > rel_std_warn (1.0) but < rel_std_fail (10.0) → warn."""
+    M = np.eye(G) * 0.01
+    M[0, 0] = 4.0  # sigma_rel = 2.0 → between 1 and 10
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "warn"
+    assert dx.n_groups_rel_std_above_1 == 1
+    assert dx.n_groups_rel_std_above_10 == 0
+    assert any("Relative uncertainty" in m for m in dx.messages)
+ 
+ 
+def test_validate_rel_std_above_fail_threshold_flagged_as_artifact():
+    """sqrt(diag) > rel_std_fail (10.0) → warn with 'artifact' message."""
+    M = np.eye(G) * 0.01
+    M[0, 0] = 200.0  # sigma_rel ~ 14 > 10
+    diags = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    dx = diags[2][2]
+    assert dx.status == "warn"
+    assert dx.n_groups_rel_std_above_10 >= 1
+    assert any("processing artifact" in m for m in dx.messages)
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — threshold customization
+# ===========================================================================
+ 
+def test_validate_custom_negativity_thresholds():
+    """Verify that relaxing thresholds turns a fail into a non-fail."""
+    eigs = np.array([1.0, 0.5, 0.2, -0.5])
+    rng = np.random.default_rng(11)
+    Q, _ = np.linalg.qr(rng.standard_normal((G, G)))
+    M = (Q * eigs) @ Q.T
+    M = 0.5 * (M + M.T)
+ 
+    # With defaults: fail.
+    strict = validate_raw_mf33_reactions(_wrap_self(2, M), EK)
+    assert strict[2][2].status == "fail"
+ 
+    # With extremely loose thresholds, the eigenvalue-spectrum branch no
+    # longer fires. Diagonal / correlation checks may still warn depending
+    # on the random rotation, but no "repaired matrix" message.
+    loose = validate_raw_mf33_reactions(
+        _wrap_self(2, M), EK,
+        negativity_ratio_warn=10.0,
+        negativity_ratio_fail=100.0,
+    )
+    assert not any(
+        "repaired matrix will differ" in m
+        for m in loose[2][2].messages
+    )
+ 
+ 
+# ===========================================================================
+# validate_raw_mf33_reactions — cross-block partner consistency
+# ===========================================================================
+ 
+def test_validate_cross_block_partner_absent_passes():
+    """ERRORR stores only the upper triangle; missing partner is expected."""
+    Mself = 0.01 * _posdef(G, seed=21)
+    Mcross = np.full((G, G), 0.001)
+    reactions = _wrap_pair(2, 102, Mself, Mcross)  # no partner at (102, 2)
+    diags = validate_raw_mf33_reactions(reactions, EK)
+    dx = diags[2][102]
+    assert dx.is_self_covariance is False
+    assert dx.status == "pass"
+    # Cross-blocks shouldn't populate self-only fields
+    assert dx.n_negative_eigenvalues is None
+    assert dx.max_abs_correlation is None
+ 
+ 
+def test_validate_cross_block_partner_consistent_passes():
+    Mself = 0.01 * _posdef(G, seed=22)
+    Mself1 = 0.01 * _posdef(G, seed=23)
+    Mcross = np.full((G, G), 0.001)
+    reactions = _wrap_pair(
+        2, 102, Mself, Mcross,
+        M_self1=Mself1,
+        M_partner=Mcross.T,   # correct partner
+    )
+    diags = validate_raw_mf33_reactions(reactions, EK)
+    assert diags[2][102].status == "pass"
+    assert diags[102][2].status == "pass"
+ 
+ 
+def test_validate_cross_block_partner_shape_mismatch_fails():
+    Mself = 0.01 * _posdef(G, seed=24)
+    Mcross = np.full((G, G), 0.001)
+    bad_partner = np.full((G, G + 1), 0.001)  # wrong shape
+    reactions = _wrap_pair(
+        2, 102, Mself, Mcross,
+        M_self1=0.01 * _posdef(G, seed=25),
+        M_partner=bad_partner,
+    )
+    diags = validate_raw_mf33_reactions(reactions, EK)
+    assert diags[2][102].status == "fail"
+    assert any("Partner wrong shape" in m for m in diags[2][102].messages)
+ 
+ 
+def test_validate_cross_block_partner_not_transpose_fails():
+    Mself = 0.01 * _posdef(G, seed=28)
+    Mcross = np.full((G, G), 0.001)
+    wrong_partner = np.full((G, G), 0.002)  # finite, right shape, wrong values
+    reactions = _wrap_pair(
+        2, 102, Mself, Mcross,
+        M_self1=0.01 * _posdef(G, seed=29),
+        M_partner=wrong_partner,
+    )
+    diags = validate_raw_mf33_reactions(reactions, EK)
+    assert diags[2][102].status == "fail"
+    msg = " ".join(diags[2][102].messages)
+    assert "!=" in msg or "C(" in msg
+