Get lpc cov

notebooks/notebook_Lpc_cov_from_enf6.ipynb

@@ -0,0 +1,112 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sandy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import seaborn as sns\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "endf6 = sandy.get_endf6_file('jeff_33','xs', 922380)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_corr(cov_obj, ax, **kwargs):\n",
+    "    cov_obj_ = sandy.CategoryCov(cov_obj) if not isinstance(cov_obj, sandy.CategoryCov) else cov_obj\n",
+    "    add = {\"cbar\": True, \"vmin\": -1, \"vmax\": 1, \"cmap\": \"RdBu\"}\n",
+    "    for k, v in kwargs.items():\n",
+    "        add[k] = v\n",
+    "    ax = sns.heatmap(cov_obj_.get_corr().data, ax=ax, **add)\n",
+    "    return ax"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot lpc correlation matrix from Endf6 files (MF=34)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cov = endf6._get_lpc_cov()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(12, 12), dpi=100)\n",
+    "plot_corr(cov, ax=ax)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot only the first polynomial for elastic cross section (mt=2) reaction"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(6, 6), dpi=100)\n",
+    "plot_corr(cov.data.loc[(9237, 2, 1), (9237, 2, 1)], ax=ax)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python3 (sandy-devel)",
+   "language": "python",
+   "name": "sandy-devel"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.10"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

sandy/core/endf6.py

@@ -2475,3 +2475,102 @@ def covariance_info(self, nubar=True, xs=True, mubar=True, chi=True):
             'chi': run_chi if chi else False,
             }
         return cov_info
+
+    def _get_lpc_cov(self, mt=range(1, 10000), mat=range(1, 10000), p=None):
+        """
+        Extract global Legendre Polynomials coefficients covariance matrix
+        from `sandy.Endf6` instance.
+
+        Parameters
+        ----------
+        mat : `int` or `list`, optional
+            MAT number. The default is range(1, 10000).
+        mt : `int` or `list`, optional
+            MT number. The default is range(1, 1000).
+        p : `int`, optional
+            maximum order of Legendre polynomial coefficients
+
+        Returns
+        -------
+        `sandy.CategoryCov`
+            Covariance matrix divided according to mat, mt, index of the
+            Legendre coefficient and energy.
+
+        Examples
+        --------
+        >>> tape = sandy.get_endf6_file("jeff_33",'xs',922380)
+        >>> out = tape._get_lpc_cov()
+        >>> out.data.index.get_level_values("MT").unique()
+        Int64Index([2], dtype='int64', name='MT')
+
+        >>> out.data.index.get_level_values("L").unique()
+        Int64Index([1, 2, 3, 4, 5, 6], dtype='int64', name='L')
+
+        Filter p:
+        >>> out = tape._get_lpc_cov(p=4)
+        >>> out.data.index.get_level_values("L").unique()
+        Int64Index([1, 2, 3, 4], dtype='int64', name='L')
+        """
+        listmt_ = [mt] if isinstance(mt, int) else mt
+        listmat_ = [mat] if isinstance(mat, int) else mat
+        tape = self.filter_by(listmf=[34],
+                              listmt=listmt_,
+                              listmat=listmat_)
+        data = []
+        for mat, mf, mt in tape.data:
+            sec = tape.read_section(mat, mf, mt)
+            for (mat1, mt1), rsec in sec["REAC"].items():
+                if mat1 == 0:
+                    mat1 = mat
+                for (l, l1), psec in rsec["P"].items():
+                    covs = []
+                    for i, nisec in psec["NI"].items():
+                        if nisec["LB"] == 5:
+                            foo = sandy.EnergyCov.from_lb5_asym if nisec["LS"] == 0 else sandy.EnergyCov.from_lb5_sym
+                            cov = foo(nisec["EK"], nisec["FKK"])
+                            covs.append(cov)
+                        elif nisec["LB"] == 1:
+                            cov = sandy.EnergyCov.from_lb1(nisec["EK"],
+                                                           nisec["FK"])
+                            covs.append(cov)
+                        elif nisec["LB"] == 2:
+                            cov = sandy.EnergyCov.from_lb2(nisec["EK"],
+                                                           nisec["FK"])
+                            covs.append(cov)
+                        elif nisec["LB"] == 6:
+                            cov = sandy.EnergyCov.from_lb6(nisec["EK"],
+                                                           nisec["EL"],
+                                                           nisec["FKL"])
+                            covs.append(cov)
+                        else:
+                            lb = nisec["LB"]
+                            logging.warning("skip LB={lb} covariance for [({mat}/{mt}), ({nat1}/{mt1})]")
+                            continue
+                    if len(covs) == 0:
+                        logging.debug("\tsubsection MAT1={}/MT1={} did not provide accetable covariances"
+                                      .format(mat1, mt1))
+                        continue
+                    cov = sandy.EnergyCov.sum_covs(*covs)
+                    index = pd.MultiIndex.from_product(
+                        [[mat], [mt], [l], cov.data.index],
+                        names=["MAT", "MT", "L", "E"],
+                        )
+                    columns = pd.MultiIndex.from_product(
+                        [[mat1], [mt1], [l1], cov.data.index],
+                        names=["MAT1", "MT1", "L1", "E1"],
+                        )
+                    df = pd.DataFrame(cov.data.values,
+                                      index=index,
+                                      columns=columns)\
+                        .stack(level=["MAT1", "MT1", "L1", "E1"])\
+                        .rename("VAL").reset_index()
+                    data.append(df)
+        if not data:
+            return pd.DataFrame()
+        data = pd.concat(data)
+        if p is not None:
+            data = data.query(f"L <= {p} & L1 <= {p}")
+        return sandy.CategoryCov.from_stack(data,
+                                            index=["MAT", "MT", "L", "E"],
+                                            columns=["MAT1", "MT1", "L1", "E1"],
+                                            values='VAL')