remove non-log class variables and their use in examples

astro-informatics · May 1, 2024 · eb84fd3 · eb84fd3
1 parent 6f674a7
commit eb84fd3
Show file tree

Hide file tree

Showing 6 changed files with 101 additions and 98 deletions.
diff --git a/examples/legacy_gaussian_diagcov.py b/examples/legacy_gaussian_diagcov.py
@@ -216,18 +216,18 @@ def run_example(
     hm.logs.debug_log("---------------------------------")
     hm.logs.debug_log(
         "Inv Evidence: analytic = {}, estimate = {}".format(
-            np.exp(ln_rho), cal_ev.evidence_inv
+            np.exp(ln_rho), np.exp(cal_ev.ln_evidence_inv)
         )
     )
     hm.logs.debug_log(
         "Inv Evidence: std = {}, std / estimate = {}".format(
-            np.sqrt(cal_ev.evidence_inv_var),
-            np.sqrt(cal_ev.evidence_inv_var) / cal_ev.evidence_inv,
+            np.sqrt(np.exp(cal_ev.ln_evidence_inv_var)),
+            np.sqrt(np.exp(cal_ev.ln_evidence_inv_var)) / np.exp(cal_ev.ln_evidence_inv),
         )
     )
     hm.logs.info_log(
         "Inv Evidence: 100 * |analytic - estimate| / estimate = {}%".format(
-            100.0 * np.abs(np.exp(ln_rho) - cal_ev.evidence_inv) / cal_ev.evidence_inv
+            100.0 * np.abs(np.exp(ln_rho) - np.exp(cal_ev.ln_evidence_inv)) / np.exp(cal_ev.ln_evidence_inv)
         )
     )
 

diff --git a/examples/legacy_gaussian_nondiagcov.py b/examples/legacy_gaussian_nondiagcov.py
@@ -186,13 +186,13 @@ def run_example(
         hm.logs.debug_log("---------------------------------")
         hm.logs.debug_log(
             "Inv Evidence: analytic = {}, estimate = {}".format(
-                np.exp(-ln_evidence_analytic), ev.evidence_inv
+                np.exp(-ln_evidence_analytic), np.exp(ev.ln_evidence_inv)
             )
         )
         hm.logs.debug_log(
             "Inv Evidence: std = {}, std / estimate = {}".format(
-                np.sqrt(ev.evidence_inv_var),
-                np.sqrt(ev.evidence_inv_var) / ev.evidence_inv,
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)),
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)) / np.exp(ev.ln_evidence_inv),
             )
         )
         hm.logs.debug_log(
@@ -202,13 +202,13 @@ def run_example(
         )
         hm.logs.debug_log(
             "Inv Evidence: sqrt( var(var) ) / var = {}".format(
-                np.sqrt(ev.evidence_inv_var_var) / ev.evidence_inv_var
+                np.sqrt(np.exp(ev.ln_evidence_inv_var_var)) / np.exp(ev.ln_evidence_inv_var)
             )
         )
         hm.logs.info_log(
             "Inv Evidence: |analytic - estimate| / estimate = {}".format(
-                np.abs(np.exp(-ln_evidence_analytic) - ev.evidence_inv)
-                / ev.evidence_inv
+                np.abs(np.exp(-ln_evidence_analytic) - np.exp(ev.ln_evidence_inv))
+                / np.exp(ev.ln_evidence_inv)
             )
         )
 
@@ -375,9 +375,9 @@ def run_example(
             plt.show(block=False)
             # ==================================================================
 
-        evidence_inv_summary[i_realisation, 0] = ev.evidence_inv
-        evidence_inv_summary[i_realisation, 1] = ev.evidence_inv_var
-        evidence_inv_summary[i_realisation, 2] = ev.evidence_inv_var_var
+        evidence_inv_summary[i_realisation, 0] = np.exp(ev.ln_evidence_inv)
+        evidence_inv_summary[i_realisation, 1] = np.exp(ev.ln_evidence_inv_var)
+        evidence_inv_summary[i_realisation, 2] = np.exp(ev.ln_evidence_inv_var_var)
 
     clock = time.process_time() - clock
     hm.logs.info_log("Execution_time = {}s".format(clock))

diff --git a/examples/legacy_rastrigin.py b/examples/legacy_rastrigin.py
@@ -270,13 +270,13 @@ def run_example(
         # ======================================================================
         hm.logs.debug_log(
             "Inv Evidence: numerical = {}, estimate = {}".format(
-                1.0 / evidence_numerical_integration, ev.evidence_inv
+                1.0 / evidence_numerical_integration, np.exp(ev.ln_evidence_inv)
             )
         )
         hm.logs.debug_log(
             "Inv Evidence: std = {}, std / estimate = {}".format(
-                np.sqrt(ev.evidence_inv_var),
-                np.sqrt(ev.evidence_inv_var) / ev.evidence_inv,
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)),
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)) / np.exp(ev.ln_evidence_inv),
             )
         )
         hm.logs.debug_log(
@@ -286,13 +286,13 @@ def run_example(
         )
         hm.logs.debug_log(
             "Inv Evidence: sqrt( var(var) )/ var = {}".format(
-                np.sqrt(ev.evidence_inv_var_var) / ev.evidence_inv_var
+                np.sqrt(np.exp(ev.ln_evidence_inv_var_var)) / np.exp(ev.ln_evidence_inv_var)
             )
         )
         hm.logs.info_log(
             "Inv Evidence: |numerical - estimate| / estimate = {}".format(
-                np.abs(1.0 / evidence_numerical_integration - ev.evidence_inv)
-                / ev.evidence_inv
+                np.abs(1.0 / evidence_numerical_integration - np.exp(ev.ln_evidence_inv))
+                / np.exp(ev.ln_evidence_inv)
             )
         )
 
@@ -409,9 +409,9 @@ def run_example(
                 created_plots = True
 
         # Save out realisations for voilin plot.
-        evidence_inv_summary[i_realisation, 0] = ev.evidence_inv
-        evidence_inv_summary[i_realisation, 1] = ev.evidence_inv_var
-        evidence_inv_summary[i_realisation, 2] = ev.evidence_inv_var_var
+        evidence_inv_summary[i_realisation, 0] = np.exp(ev.ln_evidence_inv)
+        evidence_inv_summary[i_realisation, 1] = np.exp(ev.ln_evidence_inv_var)
+        evidence_inv_summary[i_realisation, 2] = np.exp(ev.ln_evidence_inv_var_var)
 
     # ===========================================================================
     # End Timer.

diff --git a/examples/legacy_rosenbrock.py b/examples/legacy_rosenbrock.py
@@ -308,13 +308,13 @@ def run_example(
         # ======================================================================
         hm.logs.debug_log(
             "Inv Evidence: numerical = {}, estimate = {}".format(
-                1.0 / evidence_numerical_integration, ev.evidence_inv
+                1.0 / evidence_numerical_integration, np.exp(ev.ln_evidence_inv)
             )
         )
         hm.logs.debug_log(
             "Inv Evidence: std = {}, std / estimate = {}".format(
-                np.sqrt(ev.evidence_inv_var),
-                np.sqrt(ev.evidence_inv_var) / ev.evidence_inv,
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)),
+                np.sqrt(np.exp(ev.ln_evidence_inv_var)) / np.exp(ev.ln_evidence_inv),
             )
         )
         hm.logs.debug_log(
@@ -324,13 +324,13 @@ def run_example(
         )
         hm.logs.debug_log(
             "Inv Evidence: sqrt( var(var) )/ var = {}".format(
-                np.sqrt(ev.evidence_inv_var_var) / ev.evidence_inv_var
+                np.sqrt(np.exp(ev.ln_evidence_inv_var_var)) / np.exp(ev.ln_evidence_inv_var)
             )
         )
         hm.logs.info_log(
             "Inv Evidence: |numerical - estimate| / estimate = {}".format(
-                np.abs(1.0 / evidence_numerical_integration - ev.evidence_inv)
-                / ev.evidence_inv
+                np.abs(1.0 / evidence_numerical_integration - np.exp(ev.ln_evidence_inv))
+                / np.exp(ev.ln_evidence_inv)
             )
         )
 
@@ -449,9 +449,9 @@ def run_example(
             created_plots = True
 
         # Save out realisations for voilin plot.
-        evidence_inv_summary[i_realisation, 0] = ev.evidence_inv
-        evidence_inv_summary[i_realisation, 1] = ev.evidence_inv_var
-        evidence_inv_summary[i_realisation, 2] = ev.evidence_inv_var_var
+        evidence_inv_summary[i_realisation, 0] = np.exp(ev.ln_evidence_inv)
+        evidence_inv_summary[i_realisation, 1] = np.exp(ev.ln_evidence_inv_var)
+        evidence_inv_summary[i_realisation, 2] = np.exp(ev.ln_evidence_inv_var_var)
 
     # ===========================================================================
     # End Timer.

diff --git a/harmonic/evidence.py b/harmonic/evidence.py
@@ -69,9 +69,6 @@ def __init__(self, nchains: int, model, shift=Shifting.MEAN_SHIFT):
         # Chain parameters and realspace statistics
         self.nchains = nchains
         self.ndim = model.ndim
-        self.evidence_inv = 0.0
-        self.evidence_inv_var = 0.0
-        self.evidence_inv_var_var = 0.0
         self.kurtosis = 0.0
         self.n_eff = 0
 
@@ -184,13 +181,6 @@ def process_run(self):
             + np.log((kur - 1) + 2.0 / (n_eff - 1))
         )
 
-        # Compute inverse evidence statistics in real-space. In certain settings
-        # these values may return nan due to float overflow: in these cases one
-        # should use the log-space values.
-        self.evidence_inv = np.exp(self.ln_evidence_inv)
-        self.evidence_inv_var = np.exp(self.ln_evidence_inv_var)
-        self.evidence_inv_var_var = np.exp(self.ln_evidence_inv_var_var)
-
         return
 
     def get_masks(self, chain_start_ixs: jnp.ndarray) -> jnp.ndarray:
@@ -424,15 +414,24 @@ def compute_evidence(self):
                 - evidence_std (double): Estimate of standard deviation of
                   evidence.
 
+        Raises:
+
+            ValueError: if inverse evidence or its variance overflows.
         """
 
         self.check_basic_diagnostic()
 
-        common_factor = 1.0 + self.evidence_inv_var / (self.evidence_inv**2)
+        evidence_inv = np.exp(self.ln_evidence_inv)
+        evidence_inv_var = np.exp(self.ln_evidence_inv_var)
 
-        evidence = common_factor / self.evidence_inv
+        if np.isnan(evidence_inv) or np.isnan(evidence_inv_var):
+            raise ValueError("Evidence is too large to represent in non-log space. Use log-space values instead.")
 
-        evidence_std = np.sqrt(self.evidence_inv_var) / (self.evidence_inv**2)
+        common_factor = 1.0 + evidence_inv_var / (evidence_inv**2)
+
+        evidence = common_factor / evidence_inv
+
+        evidence_std = np.sqrt(evidence_inv_var) / (evidence_inv**2)
 
         return (evidence, evidence_std)
 
@@ -555,6 +554,9 @@ def compute_bayes_factor(ev1, ev2):
         ValueError: Raised if model 1 does not have chains added.
 
         ValueError: Raised if model 2 does not have chains added.
+        
+        ValueError: If inverse evidence or its variance for model 1 or model 2 too large 
+            to store in non-log space.
 
     """
 
@@ -566,14 +568,25 @@ def compute_bayes_factor(ev1, ev2):
     ev1.check_basic_diagnostic()
     ev2.check_basic_diagnostic()
 
-    common_factor = 1.0 + ev1.evidence_inv_var / (ev1.evidence_inv**2)
+    evidence_inv_ev1 = np.exp(ev1.ln_evidence_inv)
+    evidence_inv_var_ev1 = np.exp(ev1.ln_evidence_inv_var)
+
+    evidence_inv_ev2 = np.exp(ev2.ln_evidence_inv)
+    evidence_inv_var_ev2 = np.exp(ev2.ln_evidence_inv_var)
 
-    bf12 = ev2.evidence_inv / ev1.evidence_inv * common_factor
+    if np.isnan(evidence_inv_ev1) or np.isnan(evidence_inv_var_ev1):
+        raise ValueError("Evidence for model 1 is too large to represent in non-log space. Use log-space values instead.")
+    if np.isnan(evidence_inv_ev2) or np.isnan(evidence_inv_var_ev2):
+        raise ValueError("Evidence for model 2 is too large to represent in non-log space. Use log-space values instead.")
+
+    common_factor = 1.0 + evidence_inv_var_ev1 / (evidence_inv_ev1**2)
+
+    bf12 = evidence_inv_ev2 / evidence_inv_ev1 * common_factor
 
     bf12_std = np.sqrt(
-        ev1.evidence_inv**2 * ev2.evidence_inv_var
-        + ev2.evidence_inv**2 * ev1.evidence_inv_var
-    ) / (ev1.evidence_inv**2)
+        evidence_inv_ev1**2 * evidence_inv_var_ev2
+        + evidence_inv_ev2**2 * evidence_inv_var_ev1
+    ) / (evidence_inv_ev1**2)
 
     return (bf12, bf12_std)
 
@@ -612,17 +625,28 @@ def compute_ln_bayes_factor(ev1, ev2):
     ev1.check_basic_diagnostic()
     ev2.check_basic_diagnostic()
 
-    common_factor = 1.0 + ev1.evidence_inv_var / (ev1.evidence_inv**2)
+    evidence_inv_ev1 = np.exp(ev1.ln_evidence_inv)
+    evidence_inv_var_ev1 = np.exp(ev1.ln_evidence_inv_var)
+
+    evidence_inv_ev2 = np.exp(ev2.ln_evidence_inv)
+    evidence_inv_var_ev2 = np.exp(ev2.ln_evidence_inv_var)
+
+    if np.isnan(evidence_inv_ev1) or np.isnan(evidence_inv_var_ev1):
+        raise ValueError("Evidence for model 1 is too large to represent in non-log space. Use log-space values instead.")
+    if np.isnan(evidence_inv_ev2) or np.isnan(evidence_inv_var_ev2):
+        raise ValueError("Evidence for model 2 is too large to represent in non-log space. Use log-space values instead.")
+
+    common_factor = 1.0 + evidence_inv_var_ev1 / (evidence_inv_ev1**2)
 
     ln_bf12 = (
-        np.log(ev2.evidence_inv) - np.log(ev1.evidence_inv) + np.log(common_factor)
+        np.log(evidence_inv_ev2) - np.log(evidence_inv_ev1) + np.log(common_factor)
     )
 
     factor = (
-        ev1.evidence_inv**2 * ev2.evidence_inv_var
-        + ev2.evidence_inv**2 * ev1.evidence_inv_var
+        evidence_inv_ev1**2 * evidence_inv_var_ev2
+        + evidence_inv_ev2**2 * evidence_inv_var_ev1
     )
 
-    ln_bf12_std = 0.5 * np.log(factor) - 2.0 * np.log(ev1.evidence_inv)
+    ln_bf12_std = 0.5 * np.log(factor) - 2.0 * np.log(evidence_inv_ev1)
 
     return (ln_bf12, ln_bf12_std)