From 14a27dcb7668bffc99c23e12e48eae6c554e611c Mon Sep 17 00:00:00 2001
From: Giacomo Magni <giac.magni@gmail.com>
Date: Fri, 11 Nov 2022 14:59:35 +0100
Subject: [PATCH] init check on cramer rao bounds, not working

---
 src/smefit/analyze/fisher.py | 50 +++++++++++++++++++++++++++++++++++-
 src/smefit/analyze/report.py | 11 ++++++++
 2 files changed, 60 insertions(+), 1 deletion(-)

diff --git a/src/smefit/analyze/fisher.py b/src/smefit/analyze/fisher.py
index 42783944..0377b5a3 100644
--- a/src/smefit/analyze/fisher.py
+++ b/src/smefit/analyze/fisher.py
@@ -6,9 +6,12 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from rich.progress import track
 
+from ..log import logging
 from .latex_tools import latex_packages
 from .pca import impose_constrain
 
+_logger = logging.getLogger(__name__)
+
 
 class FisherCalculator:
 
@@ -66,7 +69,15 @@ def compute_linear(self):
         )
 
     def compute_quadratic(self, posterior_df, smeft_predictions):
-        """Compute quadratic Fisher information."""
+        """Compute quadratic Fisher information.
+
+        Parameters
+        ----------
+            posterior_df : pd.DataFrame
+                fit results
+            smeft_predictions: np.ndarray
+                array with all the predictions for each replica
+        """
         quad_fisher = []
 
         # compute some average values over the replicas
@@ -167,6 +178,43 @@ def normalize(table, norm, log):
             table = np.log(table[table > 0.0])
         return table.replace(np.nan, 0.0)
 
+    def test_cramer_rao_bound(self, posterior_df):
+        r"""Test Cramer Rao bound, asserting if:
+
+        .. math ::
+            I(c_{i}) - Var(c_{i}) \le 0
+
+        Parameters
+        ----------
+            posterior_df : pd.DataFrame
+                fit results
+
+        """
+        fish_mat = (
+            self.new_LinearCorrections
+            @ self.datasets.InvCovMat
+            @ self.new_LinearCorrections.T
+        )
+        fish_mat = (fish_mat + fish_mat.T) / 2
+
+        posterior_df = posterior_df[self.free_parameters]
+        covariance = posterior_df.cov().values
+        v, u = np.linalg.eig(fish_mat)
+        inv_fish = u @ np.diag(1 / v) @ u.T
+        cr_limit, _ = np.linalg.eig(inv_fish - covariance)
+        cr_limit = pd.Series(cr_limit.real, posterior_df.columns)
+
+        # cr_limit, _ = np.linalg.eig(np.eye(covariance.shape[0]) - fish_mat @ covariance)
+        # cr_limit = pd.Series(cr_limit.real, posterior_df.columns)
+        try:
+            np.testing.assert_array_less(cr_limit, np.zeros_like(cr_limit))
+            _logger.info(f"Cramer Rao bounds are satisfied!")
+        except AssertionError:
+            _logger.warning(
+                "Following coefficients violate Cramer Rao bound:  I(c)^-1 - Var(c) < 0"
+            )
+            _logger.warning(cr_limit[cr_limit > 0])
+
     def groupby_data(self, table, data_groups, norm, log):
         """Merge fisher per data group."""
         summary_table = pd.merge(
diff --git a/src/smefit/analyze/report.py b/src/smefit/analyze/report.py
index 9c520662..de605423 100644
--- a/src/smefit/analyze/report.py
+++ b/src/smefit/analyze/report.py
@@ -375,6 +375,7 @@ def fisher(
         plot=None,
         fit_list=None,
         log=False,
+        test_cramer_rao_bound=False,
     ):
         """Fisher information table and plots runner.
 
@@ -394,6 +395,8 @@ def fisher(
             By default all the fits included in the report
         log: bool, optional
             if True shows the log of the Fisher informaltion
+        test_cramer_rao_bound: bool, optinal
+            if True test the Cramer Rao bound
 
         """
         figs_list, links_list = None, None
@@ -401,9 +404,14 @@ def fisher(
             fit_list = self.fits[self.fits == fit_list]
 
         for fit in fit_list:
+            _logger.info(f"Computing Fisher for fit {fit.name}")
+
             compute_quad = fit.config["use_quad"]
             fisher_cal = FisherCalculator(fit.coefficients, fit.datasets, compute_quad)
             fisher_cal.compute_linear()
+            if test_cramer_rao_bound and not compute_quad:
+                fisher_cal.test_cramer_rao_bound(fit.results)
+
             fisher_cal.lin_fisher = fisher_cal.normalize(
                 fisher_cal.lin_fisher, norm=norm, log=log
             )
@@ -414,6 +422,9 @@ def fisher(
             # if necessary compute the quadratic Fisher
             if compute_quad:
                 fisher_cal.compute_quadratic(fit.results, fit.smeft_predictions)
+                # if test_cramer_rao_bound:
+                #     fisher_cal.test_cramer_rao_bound(fit.results)
+
                 fisher_cal.quad_fisher = fisher_cal.normalize(
                     fisher_cal.quad_fisher, norm=norm, log=log
                 )