From 6caf252a1d5f6267b32d3d3f60d2363b9a7824f0 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 11:14:58 +0100 Subject: [PATCH 01/12] add importance_reweighting method to NestedSamples --- anesthetic/samples.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 491b194f..a40b3866 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -611,6 +611,38 @@ def dlogX(self, nsamples=None): else: return WeightedDataFrame(dlogX, self.index, weights=self.weights) + def importance_reweighting(self, logL_new, replace=True, add=False): + """Perform importance re-weighting on the log-likelihood. + + Parameters + ---------- + logL_new: np.array + New log-likelihood values. Should have the same shape as `logL`. + + replace: bool + Replace the current `logL` with the new `logL_new`. + + add: bool + Add the new `logL_new` to the current `logL`. + + Returns + ------- + samples: NestedSamples + Importance re-weighted samples. + """ + samples = merge_nested_samples((self, )) + if not replace ^ add: + raise ValueError("One and only one of args 'replace' or 'add' " + "should be set to True.") + elif replace: + samples.logL = logL_new + elif add: + samples.logL += logL_new + samples = merge_nested_samples( + (samples[samples.logL > samples.logL_birth], ) + ) + return samples + def _compute_nlive(self, logL_birth): if is_int(logL_birth): nlive = logL_birth From 902526f3e0b4df7cc751d1cf33396d049b57d068 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 11:24:18 +0100 Subject: [PATCH 02/12] version bump to 2.0.0-beta.3 --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 03e296b5..c48edf1e 100644 --- a/README.rst +++ b/README.rst @@ -3,7 +3,7 @@ anesthetic: nested sampling visualisation ========================================= :anesthetic: nested sampling visualisation :Author: Will Handley and Lukas Hergt -:Version: 2.0.0-beta.2 +:Version: 2.0.0-beta.3 :Homepage: https://github.com/williamjameshandley/anesthetic :Documentation: http://anesthetic.readthedocs.io/ From bb7a013b6ecf386f9e40e96849f30ef9f2e6ad34 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 12:40:11 +0100 Subject: [PATCH 03/12] change add and replace args to action arg taking a string --- anesthetic/samples.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index a40b3866..71a464b5 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -611,7 +611,7 @@ def dlogX(self, nsamples=None): else: return WeightedDataFrame(dlogX, self.index, weights=self.weights) - def importance_reweighting(self, logL_new, replace=True, add=False): + def importance_reweighting(self, logL_new, action='replace'): """Perform importance re-weighting on the log-likelihood. Parameters @@ -619,11 +619,11 @@ def importance_reweighting(self, logL_new, replace=True, add=False): logL_new: np.array New log-likelihood values. Should have the same shape as `logL`. - replace: bool - Replace the current `logL` with the new `logL_new`. - - add: bool - Add the new `logL_new` to the current `logL`. + action: str + Can be any of {'replace', 'add'}. + * replace: Replace the current `logL` with the new `logL_new`. + * add: Add the new `logL_new` to the current `logL`. + default: 'replace' Returns ------- @@ -631,13 +631,14 @@ def importance_reweighting(self, logL_new, replace=True, add=False): Importance re-weighted samples. """ samples = merge_nested_samples((self, )) - if not replace ^ add: - raise ValueError("One and only one of args 'replace' or 'add' " - "should be set to True.") - elif replace: + if action == 'replace': samples.logL = logL_new - elif add: + elif action == 'add': samples.logL += logL_new + else: + raise NotImplementedError("`action` needs to be one of " + "{'replace', 'add'}, but '%s' was " + "requested." % action) samples = merge_nested_samples( (samples[samples.logL > samples.logL_birth], ) ) From bfd51830f904a35ac81984698231f2c620f68c92 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 12:41:19 +0100 Subject: [PATCH 04/12] rename importance_reweighting to importance_sample --- anesthetic/samples.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 71a464b5..98abd57f 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -611,7 +611,7 @@ def dlogX(self, nsamples=None): else: return WeightedDataFrame(dlogX, self.index, weights=self.weights) - def importance_reweighting(self, logL_new, action='replace'): + def importance_sample(self, logL_new, action='replace'): """Perform importance re-weighting on the log-likelihood. Parameters From fc529951dc7550cc5bf2fdc3de54d9eeb908ce03 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 13:01:44 +0100 Subject: [PATCH 05/12] append merge_nested_samples docstring to reflect additional utility for passing a single run --- anesthetic/samples.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 98abd57f..ad1d4a14 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -675,7 +675,9 @@ def merge_nested_samples(runs): Parameters ---------- runs: list(NestedSamples) - list or array-like of nested sampling runs. + List or array-like of one or more nested sampling runs. + If only a single run is provided, this recalculates the live points and + as such can be used for masked runs. Returns ------- From d129ba80a721424ab0a12244b8bd3e9700130145 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 13:06:14 +0100 Subject: [PATCH 06/12] add test for importance_samples --- tests/test_samples.py | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/test_samples.py b/tests/test_samples.py index 0d52b66c..a5b3ecbd 100644 --- a/tests/test_samples.py +++ b/tests/test_samples.py @@ -648,6 +648,40 @@ def test_posterior_points(): assert_array_equal(ns.posterior_points(0.5), ns.posterior_points(0.5)) +def test_importance_samples(): + np.random.seed(3) + ns0 = NestedSamples(root='./tests/example_data/pc') + + with pytest.raises(NotImplementedError): + ns0.importance_sample(ns0.logL, action='spam') + + ns1 = ns0.importance_sample(ns0.logL, action='replace') + assert_array_equal(ns0.logL, ns1.logL) + assert_array_equal(ns0.logL_birth, ns1.logL_birth) + assert_array_equal(ns0.weights, ns1.weights) + + ns1 = ns0.importance_sample(np.zeros_like(ns0.logL), action='add') + assert_array_equal(ns0.logL, ns1.logL) + assert_array_equal(ns0.logL_birth, ns1.logL_birth) + assert_array_equal(ns0.weights, ns1.weights) + + mask = ((ns0.x0 > -0.1) & (ns0.x2 > 0.3) & (ns0.x4 < 3)).to_numpy() + ns1 = merge_nested_samples((ns0[mask], )) + logL_new = np.where(mask, ns0.logL, -np.inf) + ns2 = ns0.importance_sample(logL_new=logL_new) + assert_array_equal(ns1.logL, ns2.logL) + assert_array_equal(ns1.logL_birth, ns2.logL_birth) + assert_array_equal(ns1.weights, ns2.weights) + NS1 = ns1.ns_output(nsamples=5000) + NS2 = ns2.ns_output(nsamples=5000) + assert NS2.logZ.mean() == pytest.approx(NS1.logZ.mean(), rel=1e-2) + assert NS2.logZ.std() == pytest.approx(NS1.logZ.std(), abs=1e-2) + assert NS2.D.mean() == pytest.approx(NS1.D.mean(), rel=1e-2) + assert NS2.D.std() == pytest.approx(NS1.D.std(), abs=1e-2) + assert NS2.d.mean() == pytest.approx(NS1.d.mean(), rel=1e-2) + assert NS2.d.std() == pytest.approx(NS1.d.std(), abs=1e-2) + + def test_wedding_cake(): np.random.seed(3) wc = WeddingCake(4, 0.5, 0.01) From 45562a086fd397d029c87e7bc467d5f06ffd8d7d Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 13:42:57 +0100 Subject: [PATCH 07/12] add action='mask' to importance_sample and corresponding test --- anesthetic/samples.py | 4 ++++ tests/test_samples.py | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index ad1d4a14..8a9db3f4 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -623,6 +623,8 @@ def importance_sample(self, logL_new, action='replace'): Can be any of {'replace', 'add'}. * replace: Replace the current `logL` with the new `logL_new`. * add: Add the new `logL_new` to the current `logL`. + * mask: treat `logL_new` as a boolean mask and only keep the + corresponding (True) samples. default: 'replace' Returns @@ -635,6 +637,8 @@ def importance_sample(self, logL_new, action='replace'): samples.logL = logL_new elif action == 'add': samples.logL += logL_new + elif action == 'mask': + samples = samples[logL_new] else: raise NotImplementedError("`action` needs to be one of " "{'replace', 'add'}, but '%s' was " diff --git a/tests/test_samples.py b/tests/test_samples.py index a5b3ecbd..4f9b310a 100644 --- a/tests/test_samples.py +++ b/tests/test_samples.py @@ -681,6 +681,11 @@ def test_importance_samples(): assert NS2.d.mean() == pytest.approx(NS1.d.mean(), rel=1e-2) assert NS2.d.std() == pytest.approx(NS1.d.std(), abs=1e-2) + ns2 = ns0.importance_sample(mask, action='mask') + assert_array_equal(ns1.logL, ns2.logL) + assert_array_equal(ns1.logL_birth, ns2.logL_birth) + assert_array_equal(ns1.weights, ns2.weights) + def test_wedding_cake(): np.random.seed(3) From c56a4290b57bf82e2d2a17b18e8f0d1decf32994 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 13:47:00 +0100 Subject: [PATCH 08/12] add 'mask' to the sets in docstring and NotImplementedError message --- anesthetic/samples.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 8a9db3f4..6f5379e1 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -620,7 +620,7 @@ def importance_sample(self, logL_new, action='replace'): New log-likelihood values. Should have the same shape as `logL`. action: str - Can be any of {'replace', 'add'}. + Can be any of {'replace', 'add', 'mask'}. * replace: Replace the current `logL` with the new `logL_new`. * add: Add the new `logL_new` to the current `logL`. * mask: treat `logL_new` as a boolean mask and only keep the @@ -641,8 +641,8 @@ def importance_sample(self, logL_new, action='replace'): samples = samples[logL_new] else: raise NotImplementedError("`action` needs to be one of " - "{'replace', 'add'}, but '%s' was " - "requested." % action) + "{'replace', 'add', 'mask'}, but '%s' " + "was requested." % action) samples = merge_nested_samples( (samples[samples.logL > samples.logL_birth], ) ) From 611a8689e747497ffaf5eca36cef73331a4427d7 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 14:15:27 +0100 Subject: [PATCH 09/12] change action default from 'replace' to 'add' --- anesthetic/samples.py | 16 ++++++++-------- tests/test_samples.py | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 6f5379e1..0ccb43b1 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -611,7 +611,7 @@ def dlogX(self, nsamples=None): else: return WeightedDataFrame(dlogX, self.index, weights=self.weights) - def importance_sample(self, logL_new, action='replace'): + def importance_sample(self, logL_new, action='add'): """Perform importance re-weighting on the log-likelihood. Parameters @@ -620,12 +620,12 @@ def importance_sample(self, logL_new, action='replace'): New log-likelihood values. Should have the same shape as `logL`. action: str - Can be any of {'replace', 'add', 'mask'}. - * replace: Replace the current `logL` with the new `logL_new`. + Can be any of {'add', 'replace', 'mask'}. * add: Add the new `logL_new` to the current `logL`. + * replace: Replace the current `logL` with the new `logL_new`. * mask: treat `logL_new` as a boolean mask and only keep the corresponding (True) samples. - default: 'replace' + default: 'add' Returns ------- @@ -633,15 +633,15 @@ def importance_sample(self, logL_new, action='replace'): Importance re-weighted samples. """ samples = merge_nested_samples((self, )) - if action == 'replace': - samples.logL = logL_new - elif action == 'add': + if action == 'add': samples.logL += logL_new + elif action == 'replace': + samples.logL = logL_new elif action == 'mask': samples = samples[logL_new] else: raise NotImplementedError("`action` needs to be one of " - "{'replace', 'add', 'mask'}, but '%s' " + "{'add', 'replace', 'mask'}, but '%s' " "was requested." % action) samples = merge_nested_samples( (samples[samples.logL > samples.logL_birth], ) diff --git a/tests/test_samples.py b/tests/test_samples.py index 4f9b310a..50b4e71e 100644 --- a/tests/test_samples.py +++ b/tests/test_samples.py @@ -667,7 +667,7 @@ def test_importance_samples(): mask = ((ns0.x0 > -0.1) & (ns0.x2 > 0.3) & (ns0.x4 < 3)).to_numpy() ns1 = merge_nested_samples((ns0[mask], )) - logL_new = np.where(mask, ns0.logL, -np.inf) + logL_new = np.where(mask, 0, -np.inf) ns2 = ns0.importance_sample(logL_new=logL_new) assert_array_equal(ns1.logL, ns2.logL) assert_array_equal(ns1.logL_birth, ns2.logL_birth) From 4bafa5ff450529889ec70b28637f22996e45b9ac Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 14:43:58 +0100 Subject: [PATCH 10/12] fix importance_sample tests --- tests/test_samples.py | 58 +++++++++++++++++++++++-------------------- 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/tests/test_samples.py b/tests/test_samples.py index 50b4e71e..b1d625b9 100644 --- a/tests/test_samples.py +++ b/tests/test_samples.py @@ -651,40 +651,44 @@ def test_posterior_points(): def test_importance_samples(): np.random.seed(3) ns0 = NestedSamples(root='./tests/example_data/pc') + pi0 = ns0.set_beta(0) + NS0 = ns0.ns_output(nsamples=2000) with pytest.raises(NotImplementedError): ns0.importance_sample(ns0.logL, action='spam') - ns1 = ns0.importance_sample(ns0.logL, action='replace') - assert_array_equal(ns0.logL, ns1.logL) - assert_array_equal(ns0.logL_birth, ns1.logL_birth) - assert_array_equal(ns0.weights, ns1.weights) + ns_masked = ns0.importance_sample(ns0.logL, action='replace') + assert_array_equal(ns0.logL, ns_masked.logL) + assert_array_equal(ns0.logL_birth, ns_masked.logL_birth) + assert_array_equal(ns0.weights, ns_masked.weights) + + ns_masked = ns0.importance_sample(np.zeros_like(ns0.logL), action='add') + assert_array_equal(ns0.logL, ns_masked.logL) + assert_array_equal(ns0.logL_birth, ns_masked.logL_birth) + assert_array_equal(ns0.weights, ns_masked.weights) + + mask = ((ns0.x0 > -0.3) & (ns0.x2 > 0.2) & (ns0.x4 < 3.5)).to_numpy() + ns_masked = merge_nested_samples((ns0[mask], )) + V_prior = pi0[mask].weights.sum() / pi0.weights.sum() + V_posterior = ns0[mask].weights.sum() / ns0.weights.sum() - ns1 = ns0.importance_sample(np.zeros_like(ns0.logL), action='add') - assert_array_equal(ns0.logL, ns1.logL) - assert_array_equal(ns0.logL_birth, ns1.logL_birth) - assert_array_equal(ns0.weights, ns1.weights) + ns1 = ns0.importance_sample(mask, action='mask') + assert_array_equal(ns_masked.logL, ns1.logL) + assert_array_equal(ns_masked.logL_birth, ns1.logL_birth) + assert_array_equal(ns_masked.weights, ns1.weights) - mask = ((ns0.x0 > -0.1) & (ns0.x2 > 0.3) & (ns0.x4 < 3)).to_numpy() - ns1 = merge_nested_samples((ns0[mask], )) logL_new = np.where(mask, 0, -np.inf) - ns2 = ns0.importance_sample(logL_new=logL_new) - assert_array_equal(ns1.logL, ns2.logL) - assert_array_equal(ns1.logL_birth, ns2.logL_birth) - assert_array_equal(ns1.weights, ns2.weights) - NS1 = ns1.ns_output(nsamples=5000) - NS2 = ns2.ns_output(nsamples=5000) - assert NS2.logZ.mean() == pytest.approx(NS1.logZ.mean(), rel=1e-2) - assert NS2.logZ.std() == pytest.approx(NS1.logZ.std(), abs=1e-2) - assert NS2.D.mean() == pytest.approx(NS1.D.mean(), rel=1e-2) - assert NS2.D.std() == pytest.approx(NS1.D.std(), abs=1e-2) - assert NS2.d.mean() == pytest.approx(NS1.d.mean(), rel=1e-2) - assert NS2.d.std() == pytest.approx(NS1.d.std(), abs=1e-2) - - ns2 = ns0.importance_sample(mask, action='mask') - assert_array_equal(ns1.logL, ns2.logL) - assert_array_equal(ns1.logL_birth, ns2.logL_birth) - assert_array_equal(ns1.weights, ns2.weights) + ns1 = ns0.importance_sample(logL_new=logL_new) + NS1 = ns1.ns_output(nsamples=2000) + assert_array_equal(ns1, ns_masked) + logZ_V = NS0.logZ.mean() + np.log(V_posterior) - np.log(V_prior) + assert abs(NS1.logZ.mean() - logZ_V) < 1.5 * NS1.logZ.std() + + logL_new = np.where(mask, 0, -1e30) + ns1 = ns0.importance_sample(logL_new=logL_new) + NS1 = ns1.ns_output(nsamples=2000) + logZ_V = NS0.logZ.mean() + np.log(V_posterior) + assert abs(NS1.logZ.mean() - logZ_V) < 1.5 * NS1.logZ.std() def test_wedding_cake(): From f4d315a11c9113b55806c13744c612d709efdc4d Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 15:11:49 +0100 Subject: [PATCH 11/12] move importance_sample to MCMCSamples and inherit in NestedSamples --- anesthetic/samples.py | 47 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 0ccb43b1..3cc8c417 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -341,6 +341,40 @@ def plot_2d(self, axes, *args, **kwargs): return fig, axes + def importance_sample(self, logL_new, action='add'): + """Perform importance re-weighting on the log-likelihood. + + Parameters + ---------- + logL_new: np.array + New log-likelihood values. Should have the same shape as `logL`. + + action: str + Can be any of {'add', 'replace', 'mask'}. + * add: Add the new `logL_new` to the current `logL`. + * replace: Replace the current `logL` with the new `logL_new`. + * mask: treat `logL_new` as a boolean mask and only keep the + corresponding (True) samples. + default: 'add' + + Returns + ------- + samples: MCMCSamples + Importance re-weighted samples. + """ + samples = self.copy() + if action == 'add': + samples.logL += logL_new + elif action == 'replace': + samples.logL = logL_new + elif action == 'mask': + samples = samples[logL_new] + else: + raise NotImplementedError("`action` needs to be one of " + "{'add', 'replace', 'mask'}, but '%s' " + "was requested." % action) + return samples + def _limits(self, paramname): limits = self.limits.get(paramname, (None, None)) if limits[0] == limits[1]: @@ -633,16 +667,9 @@ def importance_sample(self, logL_new, action='add'): Importance re-weighted samples. """ samples = merge_nested_samples((self, )) - if action == 'add': - samples.logL += logL_new - elif action == 'replace': - samples.logL = logL_new - elif action == 'mask': - samples = samples[logL_new] - else: - raise NotImplementedError("`action` needs to be one of " - "{'add', 'replace', 'mask'}, but '%s' " - "was requested." % action) + samples = super(NestedSamples, samples).importance_sample( + logL_new=logL_new, action=action + ) samples = merge_nested_samples( (samples[samples.logL > samples.logL_birth], ) ) From 955776a27072b5f4542186b41fbd602c36fdb7a8 Mon Sep 17 00:00:00 2001 From: lukashergt Date: Thu, 20 Aug 2020 15:15:00 +0100 Subject: [PATCH 12/12] fix merge_nested_samples docstring: two->one --- anesthetic/samples.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/anesthetic/samples.py b/anesthetic/samples.py index 3cc8c417..c4c4dec3 100644 --- a/anesthetic/samples.py +++ b/anesthetic/samples.py @@ -701,7 +701,7 @@ def _constructor(self): def merge_nested_samples(runs): - """Merge two or more nested sampling runs. + """Merge one or more nested sampling runs. Parameters ----------