Improve TUD W post-processing and try fixing the PR templates

.github/PULL_REQUEST_TEMPLATE/new_benchmark.md

@@ -24,8 +24,9 @@ Please select the transport code for which the benchmark is added.
 
 ## Checklist:
 
-- [ ] I have added the required YAML files.
+- [ ] I have added the required YAML files for post-processing.
 - [ ] I have updated the documentation to include the benchmark description.
 - [ ] I have updated the documentation summary table which gives an overview of the implemented benchmark. 
-- [ ] Have the benchmark inputs been added to the appropriate respository. 
+- [ ] The benchmark is available in the default ``run_cfg.yaml`` file
+- [ ] The benchmark inputs have been added to the appropriate repository. 
 - [ ] My changes generate no new warnings.

docs/source/dev/insertbenchmarks.rst

@@ -248,9 +248,9 @@ The currently supported modifiers are:
   If no *column* argument is provided, the cumulative sum is computed on the 'Value' column by default. The argument *norm* is True by default.
 
 * ``gaussian_broadening``: applies Gaussian broadening to the 'Value' column. The optional keyarg to provide is *fwhm_frac*, 
-which specifies the fraction of the FWHM (Full Width at Half Maximum) to use for Gaussian broadening. This can be provided either 
-as a single float value, which will be applied uniformly to all energy bins, or as a list of float values with the same length as the 
-'Energy' column, allowing for a different broadening parameter for each energy bin. If not specified, the default value is 0.1 (10%).
+  which specifies the fraction of the FWHM (Full Width at Half Maximum) to use for Gaussian broadening. This can be provided either 
+  as a single float value, which will be applied uniformly to all energy bins, or as a list of float values with the same length as the 
+  'Energy' column, allowing for a different broadening parameter for each energy bin. If not specified, the default value is 0.1 (10%).
 
 More than one modifiers can be applied in series to a single tally.
 If your benchmark requires a new modifier, please refer to :ref:`add_tally_mod`.

src/jade/app/app.py

@@ -18,7 +18,7 @@
 from jade.config.raw_config import ConfigRawProcessor
 from jade.config.run_config import RunConfig, RunMode
 from jade.config.status import GlobalStatus
-from jade.helper.__optionals__ import TKINTER_AVAIL
+from jade.helper.__optionals__ import TKINTER_AVAIL, OMC_AVAIL
 
 if TKINTER_AVAIL:
     from jade.gui.post_config_gui import PostConfigGUI
@@ -223,6 +223,11 @@ def get_config(
 
         to_process = {}
         for code, lib, bench in successful:
+            # if openmc is not available in system, skip processing
+            if code == CODE.OPENMC and not OMC_AVAIL:
+                logging.warning(f"OpenMC not installed. Skipping {bench} processing.")
+                continue
+
             if force:
                 # process all the successful simulations, force override
                 raw_cfg = get_config(root_cfg, code, bench)
@@ -235,6 +240,7 @@ def get_config(
                 if (code, lib, bench) not in self.status.raw_data:
                     if subset is not None and bench not in subset:
                         continue
+
                     raw_cfg = get_config(root_cfg, code, bench)
                     if raw_cfg is None:
                         continue

src/jade/helper/aux_functions.py

@@ -6,6 +6,8 @@
 from pathlib import Path
 from typing import TYPE_CHECKING, Union
 
+import numpy as np
+import pandas as pd
 import yaml
 
 if TYPE_CHECKING:
@@ -135,3 +137,26 @@ class VerboseSafeDumper(yaml.SafeDumper):
 
     def ignore_aliases(self, data):
         return True
+
+
+def same_index(
+    index1: pd.Index | pd.MultiIndex, index2: pd.Index | pd.MultiIndex
+) -> bool:
+    """Check if two pandas indices are the same, allowing for small numerical differences."""
+    if index1.nlevels != index2.nlevels:
+        return False
+    for level in range(index1.nlevels):
+        subindex1 = index1.get_level_values(level)
+        subindex2 = index2.get_level_values(level)
+        if not subindex1.equals(subindex2):
+            # Check if both indices are numeric and match within a tolerance
+            try:
+                idx1 = np.array(subindex1, dtype=float)
+                idx2 = np.array(subindex2, dtype=float)
+                if np.allclose(idx1, idx2, rtol=1e-3):
+                    continue
+                else:
+                    return False
+            except Exception:
+                return False
+    return True

src/jade/post/plotter.py

@@ -1,14 +1,16 @@
 from __future__ import annotations
 
+import logging
 import math
+import warnings
 from abc import ABC, abstractmethod
 
+import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import seaborn as sns
 from f4enix.input.libmanager import LibManager
-import matplotlib
 from matplotlib.axes import Axes
 from matplotlib.figure import Figure
 from matplotlib.lines import Line2D
@@ -17,6 +19,8 @@
 from matplotlib.ticker import AutoLocator, AutoMinorLocator, LogLocator, MultipleLocator
 
 from jade.config.atlas_config import PlotConfig, PlotType
+from jade.helper.aux_functions import same_index
+from jade.post.manipulate_tally import ComparisonType, compare_data
 
 matplotlib.use("Agg")  # use a non-interactive backend
 LM = LibManager()
@@ -523,41 +527,57 @@ def _get_figure(self) -> tuple[Figure, list[Axes]]:
             raise ValueError(f"Style {style} not recognized")
 
         # compute the ratios
+        to_plot = []
         if subcases:
-            to_plot = []
-            for codelib, df in self.data[1:]:
-                to_plot.append(
-                    (
-                        codelib,
-                        df.set_index([subcases[0], self.cfg.x])[self.cfg.y]
-                        / self.data[0][1].set_index([subcases[0], self.cfg.x])[
-                            self.cfg.y
-                        ],
-                    )
-                )
+            ref = self.data[0][1].set_index([subcases[0], self.cfg.x])
         else:
-            to_plot = [
-                (
-                    codelib,
-                    df.set_index(self.cfg.x)[self.cfg.y]
-                    / self.data[0][1].set_index(self.cfg.x)[self.cfg.y],
+            ref = self.data[0][1].set_index(self.cfg.x)
+        val1 = ref[self.cfg.y].sort_index()
+        err1 = ref["Error"].sort_index()
+
+        for codelib, df in self.data[1:]:
+            if subcases:
+                target = df.set_index([subcases[0], self.cfg.x])
+            else:
+                target = df.set_index(self.cfg.x)
+            val2 = target[self.cfg.y].sort_index()
+            err2 = target["Error"].sort_index()
+            # sometimes there are index which are numerical and may have slight
+            # differences due to rounding. In reality the two must be the same
+            # in a C/E plot
+            if same_index(val1.index, val2.index) is False:
+                logging.error(
+                    f"Indices do not match between reference and {codelib}: "
+                    f"{val1.index}, {val2.index}"
                 )
-                for (codelib, df) in self.data[1:]
-            ]
+                raise RuntimeError("Indices do not match.")
+            else:
+                val2.index = val1.index
+                err2.index = err1.index
+
+            values, errors = compare_data(
+                val1,
+                val2,
+                err1,
+                err2,
+                comparison_type=ComparisonType.RATIO,
+            )
+            to_plot.append((codelib, values, errors))
 
         # Plot the data
-        for idx, (codelib, df) in enumerate(to_plot):
+        for idx, (codelib, df_vals, df_errors) in enumerate(to_plot):
             # Split the dfs into the subcases if needed
             if subcases:
                 dfs = []
                 for value in subcases[1]:
                     try:
-                        subset = df.loc[value]
+                        subset_val = df_vals.loc[value]
+                        subset_err = df_errors.loc[value]
                     except KeyError:
                         continue
-                    dfs.append((value, subset))
+                    dfs.append((value, subset_val, subset_err))
             else:
-                dfs = [(None, df)]
+                dfs = [(None, df_vals, df_errors)]
 
             # If this is the first lib, create the plot
             if idx == 0:
@@ -571,7 +591,7 @@ def _get_figure(self) -> tuple[Figure, list[Axes]]:
                     axes = ax
 
             # plot all subcases
-            for i, (case, df1) in enumerate(dfs):
+            for i, (case, dfv, dfe) in enumerate(dfs):
                 if i == 0:
                     label = codelib
                 else:
@@ -593,8 +613,8 @@ def _get_figure(self) -> tuple[Figure, list[Axes]]:
 
                 if style == "step":
                     axes[i].step(
-                        df1.index,
-                        df1.values,
+                        dfv.index,
+                        dfv.values,
                         label=label,
                         color=COLORS[idx],
                         linestyle=LINESTYLES[idx],
@@ -606,22 +626,40 @@ def _get_figure(self) -> tuple[Figure, list[Axes]]:
                         _apply_CE_limits(
                             ce_limits[0],
                             ce_limits[1],
-                            df1.values,
-                            df1.index,
+                            dfv.values,
+                            dfv.index,
                             axes[i],
                             idx,
                             label,
                         )
                     else:
-                        ax.scatter(
-                            df1.index,
-                            df1.values,
-                            label=label,
-                            color=COLORS[idx],
-                            marker=MARKERS[idx],
-                            # the marker should be not filled
-                            facecolors="none",
-                        )
+                        with warnings.catch_warnings():
+                            warnings.simplefilter("ignore", category=UserWarning)
+                            axes[i].scatter(
+                                dfv.index,
+                                dfv.values,
+                                label=label,
+                                color=COLORS[idx],
+                                marker=MARKERS[idx],
+                                # the marker should be not filled
+                                facecolors="none",
+                            )
+                    # add error bars
+                    axes[i].errorbar(
+                        dfv.index,
+                        dfv.values,
+                        yerr=dfe.values * dfv.values,
+                        # fmt="none",
+                        ecolor=COLORS[idx],
+                        # elinewidth=0.5,
+                        # capsize=2,
+                        # label=None,
+                    )
+                    # if it is a scatter plot we must check for categorical X axis
+                    # as this is not automatically detected by matplotlib
+                    if dfv.index.dtype == str or dfv.index.dtype == object:
+                        axes[i].set_xticks(range(len(dfv.index)))
+                        axes[i].set_xticklabels(dfv.index)
 
         # put the legend in the top right corner if it was not already placed
         if not axes[0].get_legend():
@@ -1003,21 +1041,25 @@ def _apply_CE_limits(
         alpha=0,
     )
 
-    # normal points
-    ax.scatter(
-        norm[0],
-        norm[1],
-        label=label,
-        color=COLORS[idx],
-        marker=MARKERS[idx],
-        # the marker should be not filled
-        facecolors="none",
-    )
-    # upper and lower limits
-    ax.scatter(upper[0], upper[1], marker=CARETUPBASE, c=COLORS[idx], facecolors="none")
-    ax.scatter(
-        lower[0], lower[1], marker=CARETDOWNBASE, c=COLORS[idx], facecolors="none"
-    )
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", category=UserWarning)
+        # normal points
+        ax.scatter(
+            norm[0],
+            norm[1],
+            label=label,
+            color=COLORS[idx],
+            marker=MARKERS[idx],
+            # the marker should be not filled
+            facecolors="none",
+        )
+        # upper and lower limits
+        ax.scatter(
+            upper[0], upper[1], marker=CARETUPBASE, c=COLORS[idx], facecolors="none"
+        )
+        ax.scatter(
+            lower[0], lower[1], marker=CARETDOWNBASE, c=COLORS[idx], facecolors="none"
+        )
     # additional legend
     leg = [
         Line2D(
@@ -1045,7 +1087,7 @@ def _apply_CE_limits(
     combined = handles + leg
 
     if label is not None:
-        ax.legend(handles=combined, loc="best")
+        ax.legend(handles=combined, bbox_to_anchor=(1, 1))
 
 
 def _rotate_ticks(ax: Axes) -> None: