From eb84fd399b1d3b1f61345f8a7d5d6aa46f5be75a Mon Sep 17 00:00:00 2001 From: CosmoMatt Date: Wed, 1 May 2024 15:33:44 +0100 Subject: [PATCH] remove non-log class variables and their use in examples --- examples/legacy_gaussian_diagcov.py | 8 +-- examples/legacy_gaussian_nondiagcov.py | 18 +++---- examples/legacy_rastrigin.py | 18 +++---- examples/legacy_rosenbrock.py | 18 +++---- harmonic/evidence.py | 70 +++++++++++++++++--------- tests/test_evidence.py | 67 +++++++++--------------- 6 files changed, 101 insertions(+), 98 deletions(-) diff --git a/examples/legacy_gaussian_diagcov.py b/examples/legacy_gaussian_diagcov.py index 3f328262..bd639119 100644 --- 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 index d80a5654..3d9254d3 100644 --- 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 index 69ddc58a..6ee86e1f 100644 --- 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 index b5b63629..7e6fbc74 100644 --- 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 index 20552d19..13cb794b 100644 --- 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) diff --git a/tests/test_evidence.py b/tests/test_evidence.py index ad8a21c0..dd8c5b88 100644 --- a/tests/test_evidence.py +++ b/tests/test_evidence.py @@ -32,9 +32,9 @@ def test_constructor(model): rho = cbe.Evidence(nchains, model) assert rho.nchains == nchains - assert rho.evidence_inv == pytest.approx(0.0) - assert rho.evidence_inv_var == pytest.approx(0.0) - assert rho.evidence_inv_var_var == pytest.approx(0.0) + assert rho.ln_evidence_inv == pytest.approx(0.0) + assert rho.ln_evidence_inv_var == pytest.approx(0.0) + assert rho.ln_evidence_inv_var_var == pytest.approx(0.0) assert rho.running_sum.size == nchains assert rho.nsamples_per_chain.size == nchains assert rho.shift_value == pytest.approx(0.0) @@ -82,9 +82,9 @@ def test_process_run_with_shift(model): / nchains**3 ) - assert rho.evidence_inv == pytest.approx(evidence_inv, abs=1e-7) - assert rho.evidence_inv_var == pytest.approx(evidence_inv_var) - assert rho.evidence_inv_var_var == pytest.approx(evidence_inv_var_var) + assert np.exp(rho.ln_evidence_inv) == pytest.approx(evidence_inv, abs=1e-7) + assert np.exp(rho.ln_evidence_inv_var) == pytest.approx(evidence_inv_var) + assert np.exp(rho.ln_evidence_inv_var_var) == pytest.approx(evidence_inv_var_var) rho = cbe.Evidence(nchains, model, cbe.Shifting.MEAN_SHIFT) np.random.seed(1) @@ -106,9 +106,9 @@ def test_process_run_with_shift(model): / nchains**3 ) - assert rho.evidence_inv == pytest.approx(evidence_inv, abs=1e-7) - assert rho.evidence_inv_var == pytest.approx(evidence_inv_var) - assert rho.evidence_inv_var_var == pytest.approx(evidence_inv_var_var) + assert np.exp(rho.ln_evidence_inv) == pytest.approx(evidence_inv, abs=1e-7) + assert np.exp(rho.ln_evidence_inv_var) == pytest.approx(evidence_inv_var) + assert np.exp(rho.ln_evidence_inv_var_var) == pytest.approx(evidence_inv_var_var) def test_add_chains(): @@ -134,11 +134,11 @@ def test_add_chains(): cal_ev = cbe.Evidence(nchains, sphere, cbe.Shifting.MEAN_SHIFT) cal_ev.add_chains(chain) - print("cal_ev.evidence_inv = {}".format(cal_ev.evidence_inv)) + print("cal_ev.evidence_inv = {}".format(np.exp(cal_ev.ln_evidence_inv))) - assert cal_ev.evidence_inv == pytest.approx(0.159438606) - assert cal_ev.evidence_inv_var == pytest.approx(1.164628268e-07) - assert cal_ev.evidence_inv_var_var**0.5 == pytest.approx(1.142786462e-08) + assert np.exp(cal_ev.ln_evidence_inv) == pytest.approx(0.159438606) + assert np.exp(cal_ev.ln_evidence_inv_var) == pytest.approx(1.164628268e-07) + assert np.exp(cal_ev.ln_evidence_inv_var_var)**0.5 == pytest.approx(1.142786462e-08) nsamples1 = 300 chains1 = ch.Chains(ndim) @@ -153,9 +153,9 @@ def test_add_chains(): ev.add_chains(chains1) ev.add_chains(chains2) - assert ev.evidence_inv == pytest.approx(0.159438606) - assert ev.evidence_inv_var == pytest.approx(1.164628268e-07) - assert ev.evidence_inv_var_var**0.5 == pytest.approx(1.142786462e-08) + assert np.exp(ev.ln_evidence_inv) == pytest.approx(0.159438606) + assert np.exp(ev.ln_evidence_inv_var) == pytest.approx(1.164628268e-07) + assert np.exp(ev.ln_evidence_inv_var_var)**0.5 == pytest.approx(1.142786462e-08) return @@ -223,8 +223,6 @@ def test_compute_evidence(model): ev_inv = 1e10 ev_inv_var = 2e10 ev = cbe.Evidence(nchains, model) - ev.evidence_inv = ev_inv - ev.evidence_inv_var = ev_inv_var ev.ln_evidence_inv = np.log(ev_inv) ev.ln_evidence_inv_var = np.log(ev_inv_var) @@ -248,8 +246,6 @@ def test_compute_ln_inv_evidence_errors(model): ln_ev_inv = 10 ln_ev_inv_var = 2 * ln_ev_inv ev = cbe.Evidence(nchains, model) - ev.evidence_inv = np.exp(ln_ev_inv) - ev.evidence_inv_var = np.exp(ln_ev_inv_var) ev.ln_evidence_inv = ln_ev_inv ev.ln_evidence_inv_var = ln_ev_inv_var @@ -261,8 +257,6 @@ def test_compute_ln_inv_evidence_errors(model): ln_ev_inv = 10 ln_ev_inv_var = ln_ev_inv ev = cbe.Evidence(nchains, model) - ev.evidence_inv = np.exp(ln_ev_inv) - ev.evidence_inv_var = np.exp(ln_ev_inv_var) ev.ln_evidence_inv = ln_ev_inv ev.ln_evidence_inv_var = ln_ev_inv_var @@ -278,8 +272,6 @@ def test_compute_ln_inv_evidence_errors(model): ln_ev_inv = 10 ln_ev_inv_var = 0.5 * ln_ev_inv ev = cbe.Evidence(nchains, model) - ev.evidence_inv = np.exp(ln_ev_inv) - ev.evidence_inv_var = np.exp(ln_ev_inv_var) ev.ln_evidence_inv = ln_ev_inv ev.ln_evidence_inv_var = ln_ev_inv_var @@ -303,8 +295,8 @@ def test_compute_bayes_factors(): ev1_inv = 1e10 ev1_inv_var = 2e10 ev1 = cbe.Evidence(nchains, sphere) - ev1.evidence_inv = ev1_inv - ev1.evidence_inv_var = ev1_inv_var + ev1.ln_evidence_inv = np.log(ev1_inv) + ev1.ln_evidence_inv_var = np.log(ev1_inv_var) ev1.chains_added = True ndim = 4 @@ -313,8 +305,8 @@ def test_compute_bayes_factors(): ev2_inv = 3e10 ev2_inv_var = 4e10 ev2 = cbe.Evidence(nchains, sphere) - ev2.evidence_inv = ev2_inv - ev2.evidence_inv_var = ev2_inv_var + ev2.ln_evidence_inv = np.log(ev2_inv) + ev2.ln_evidence_inv_var = np.log(ev2_inv_var) ev2.chains_added = True bf12_check = ev2_inv / ev1_inv * (1.0 + ev2_inv_var / ev2_inv**2) @@ -332,19 +324,6 @@ def test_compute_bayes_factors(): assert bf12 == pytest.approx(np.exp(ln_bf12)) assert bf12_std == pytest.approx(np.exp(ln_bf12_std)) - # Test bayes factor reduces to single evidence calculation. - ev2_inv = 1.0 - ev2_inv_var = 0.0 - ev2 = cbe.Evidence(nchains, sphere) - ev2.evidence_inv = ev2_inv - ev2.evidence_inv_var = ev2_inv_var - ev2.chains_added = True - (bf12, bf12_std) = cbe.compute_bayes_factor(ev1, ev2) - - (evidence, evidence_std) = ev1.compute_evidence() - assert bf12 == pytest.approx(evidence) - assert bf12_std == pytest.approx(evidence_std) - @pytest.mark.parametrize("model", models_to_test_2) def test_serialization(model): @@ -379,9 +358,9 @@ def test_serialization(model): # Test evidence objects the same assert ev1.nchains == ev2.nchains - assert ev1.evidence_inv == ev2.evidence_inv - assert ev1.evidence_inv_var == ev2.evidence_inv_var - assert ev1.evidence_inv_var_var == ev2.evidence_inv_var_var + assert ev1.ln_evidence_inv == ev2.ln_evidence_inv + assert ev1.ln_evidence_inv_var == ev2.ln_evidence_inv_var + assert ev1.ln_evidence_inv_var_var == ev2.ln_evidence_inv_var_var assert ev1.running_sum.size == ev2.running_sum.size assert ev1.nsamples_per_chain.size == ev2.nsamples_per_chain.size assert ev1.shift_value == ev2.shift_value