plot: add pull plots, "extra" option for wrapper

examples/002_model_functions/model_functions.py

@@ -40,13 +40,15 @@ def exponential_model(x, A_0=1., x_0=5.):
 
 x_data = [0.38, 0.83, 1.96, 2.82, 4.28, 4.69, 5.97, 7.60, 7.62, 8.81, 9.87, 10.91]
 y_data = [1.60, 1.66, 2.12, 3.05, 3.57, 4.65, 6.21, 7.57, 8.27, 10.79, 14.27, 18.48]
-x_error = 0.3
-y_error_rel = 0.05
+x_error = 0.2  # 0.2 in absolute units of x
+y_error = 0.5  # 0.5 in absolute units of y
+y_error_rel = 0.03  # 3% of the y model value
 
 # kafe2.xy_fit needs to be called twice to do two fits:
-kafe2.xy_fit(linear_model, x_data, y_data, x_error=x_error, y_error_rel=y_error_rel)
+kafe2.xy_fit(linear_model, x_data, y_data,
+             x_error=x_error, y_error=y_error, y_error_rel=y_error_rel)
 kafe2.xy_fit(exponential_model, x_data, y_data,
-             x_error=x_error, y_error_rel=y_error_rel, profile=True)
+             x_error=x_error, y_error=y_error, y_error_rel=y_error_rel)
 # Make sure to specify profile=True whenever you use a nonlinear model function.
 # A model function is linear if it is a linear function of each of its parameters.
 # The model function does not need to be a linear function of the independent variable x.
@@ -62,7 +64,11 @@ def exponential_model(x, A_0=1., x_0=5.):
 
     # When Python functions are used as custom model functions kafe2 does not know
     # how to express them as LaTeX. The LaTeX can be manually defined like this:
-    model_expression=["{a}{x} + {b}", "{A_0} e^{{{x}/{x_0}}}"]
+    model_expression=["{a}{x} + {b}", "{A_0} e^{{{x}/{x_0}}}"],
     # Parameter names have to be put between {}. To get {} for LaTex double them like {{ or }}.
     # When using SymPy to define model function the LaTeX expression can be derived automatically.
+
+    # Add a so-called pull plot in order to compare the influence of individual data points on the models.
+    extra="pull"
+    # Alternative extra plots: "residual" or "ratio".
 )

kafe2/config/kafe2.yaml

@@ -30,6 +30,7 @@ fit:
       y: '$y$'
     ratio_label: 'Ratio'
     residual_label: 'Residual'
+    pull_label: 'Pull'
     error_label: "%(model_label)s $\\pm 1\\sigma$"
 
     style: !include plot_style_color.yaml

kafe2/fit/_base/plot.py

@@ -163,6 +163,11 @@ class PlotAdapterBase:
             plot_adapter_method="plot_residual",
             target_axes="residual",
         ),
+        pull=dict(
+            plot_style_as="data",
+            plot_adapter_method="plot_pull",
+            target_axes="pull",
+        ),
     )
 
     AVAILABLE_X_SCALES = ("linear",)
@@ -626,6 +631,27 @@ def plot_residual(self, target_axes, error_contributions=("data",), **kwargs):
             **kwargs,
         )
 
+    def plot_pull(self, target_axes, error_contributions=("data",), **kwargs):
+        """Plot the pull of the data values on the model to a specified :py:obj:`matplotlib.axes.Axes` object.
+
+        :param matplotlib.axes.Axes target_axes: The :py:obj:`matplotlib` axes used for plotting.
+        :param error_contributions: Which error contributions to include when plotting the data.
+            Can either be ``data``, ``'model'`` or both.
+        :type error_contributions: str or Tuple[str]
+        :param dict kwargs: Keyword arguments accepted by :py:obj:`matplotlib.pyplot.errorbar`.
+        :return: plot handle(s)
+        """
+
+        _xmin, _xmax = self.x_range
+        target_axes.fill_between([_xmin, _xmax], -2, 2, color=[0.0, 0.0, 0.0, 0.1])
+        target_axes.fill_between([_xmin, _xmax], -1, 1, color=[0.0, 0.0, 0.0, 0.1])
+        target_axes.hlines(y=0, xmin=_xmin, xmax=_xmax, colors="black", linestyles=":")
+
+        _yerr = self._get_total_error(error_contributions)
+        _pull = (self.data_y - self.model_y) / _yerr
+
+        return target_axes.errorbar(self.data_x, _pull, xerr=0, yerr=0, **kwargs)
+
     # Overridden by multi plot adapters
     def get_formatted_model_function(self, **kwargs):
         """return model function string"""
@@ -1240,6 +1266,9 @@ def plot(
         residual=False,
         residual_range=None,
         residual_height_share=0.25,
+        pull=False,
+        pull_range=None,
+        pull_height_share=0.25,
         plot_width_share=0.5,
         font_scale=1.0,
         figsize=None,
@@ -1268,8 +1297,15 @@ def plot(
         :return: dictionary containing information about the plotted objects
         :rtype: dict
         """
-        if ratio and residual:
-            raise NotImplementedError("Cannot plot ratio and residual at the same time.")
+        _num_extra_plots = 0
+        if ratio:
+            _num_extra_plots += 1
+        if residual:
+            _num_extra_plots += 1
+        if pull:
+            _num_extra_plots += 1
+        if _num_extra_plots > 1:
+            raise NotImplementedError("Only one out of ratio, residual, and pull can be used at the same time.")
 
         with rc_context(kafe2_rc):
             rcParams["font.size"] *= font_scale
@@ -1281,10 +1317,14 @@ def plot(
                 _axes_keys += ("ratio",)
                 _height_ratios[0] -= ratio_height_share
                 _height_ratios.append(ratio_height_share)
-            elif residual:
+            if residual:
                 _axes_keys += ("residual",)
                 _height_ratios[0] -= residual_height_share
                 _height_ratios.append(residual_height_share)
+            if pull:
+                _axes_keys += ("pull",)
+                _height_ratios[0] -= pull_height_share
+                _height_ratios.append(pull_height_share)
 
             _all_plot_results = []
             for i in range(len(self._fits) if self._separate_figs else 1):
@@ -1356,6 +1396,24 @@ def plot(
                         _axis.set_ylim((_low, _high))
                     else:
                         _axis.set_ylim(residual_range)
+                if pull:
+                    _axis = self._current_axes["pull"]
+                    _pull_label = kc("fit", "plot", "pull_label")
+                    _axis.set_ylabel(_pull_label)
+                    if pull_range is None:
+                        _plot_adapters = self._get_plot_adapters()[i : i + 1] if self._separate_figs else self._get_plot_adapters()
+                        _max_abs_deviation = 2 # Ensure 1 sigma/2 sigma bands are visible
+                        for _plot_adapter in _plot_adapters:
+                            _max_abs_deviation = max(
+                                _max_abs_deviation,
+                                np.max(np.abs((_plot_adapter.data_y - _plot_adapter.model_y) / _plot_adapter.data_yerr)),
+                            )
+                        # Small gap between highest error bar and plot border:
+                        _low = -_max_abs_deviation * 1.2
+                        _high = _max_abs_deviation * 1.2
+                        _axis.set_ylim((_low, _high))
+                    else:
+                        _axis.set_ylim(pull_range)
 
                 _all_plot_results.append(_plot_results)
 

kafe2/fit/util/wrapper.py

@@ -477,6 +477,7 @@ def plot(
     legend=True,
     fit_info=True,
     error_band=True,
+    extra=None,
     profile=None,
     plot_profile=None,
     show=True,
@@ -527,6 +528,8 @@ def plot(
     :type fit_info: bool
     :param error_band: whether the model error band should be shown.
     :type error_band: bool
+    :param extra: additional, supplementary plots to show below the main plot.
+    :type error_band: "ratio", "residual", or "pull".
     :param profile: whether the profile likelihood method should be used for asymmetric parameter
         errors and profile/contour plots.
     :type profile: bool
@@ -648,12 +651,17 @@ def plot(
     if y_ticks is not None:
         _plot.y_ticks = y_ticks
 
+    if extra not in ["ratio", "residual", "pull"]:
+        raise ValueError(f"Unknown extra plot: '{extra}'. Available: 'ratio', 'residual', 'pull'.")
     if len(fits) > 0:  # Do not plot if only CustomFit.
         _plot.plot(
             legend=legend,
             fit_info=fit_info,
             asymmetric_parameter_errors=profile,
             font_scale=font_scale,
+            ratio=extra == "ratio",
+            residual=extra == "residual",
+            pull=extra == "pull",
         )
 
         if save: