From 516b0a285b074d47b2a30cf0462af5676ad0102d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 16:56:35 -0500 Subject: [PATCH 01/30] implement wideband likelihood --- src/pint/bayesian.py | 65 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 51 insertions(+), 14 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index afdf1ff80..3ba070813 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -6,7 +6,7 @@ from scipy.stats import norm, uniform from pint.models.priors import Prior, UniformUnboundedRV -from pint.residuals import Residuals +from pint.residuals import Residuals, WidebandTOAResiduals class BayesianTiming: @@ -54,8 +54,10 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self.model = deepcopy(model) self.toas = toas - if toas.is_wideband(): - raise NotImplementedError("Wideband TOAs are not yet supported.") + self.is_wideband = toas.is_wideband() + + # if toas.is_wideband(): + # raise NotImplementedError("Wideband TOAs are not yet supported.") self.param_labels = self.model.free_params self.params = [getattr(self.model, par) for par in self.param_labels] @@ -91,11 +93,12 @@ def _validate_priors(self): ) def _decide_likelihood_method(self): - """Weighted least squares with normalization term (wls), or Generalized - least-squares with normalization term (gls).""" + """Weighted least squares with normalization term (wls), or Generalized leasts + quares with normalization term (gls), for narrow-band (nb) or wide-band (wb) + dataset.""" if "NoiseComponent" not in self.model.component_types: - return "wls" + method = "wls" else: correlated_errors_present = np.any( [ @@ -104,13 +107,17 @@ def _decide_likelihood_method(self): ] ) if not correlated_errors_present: - return "wls" + method = "wls" else: raise NotImplementedError( "GLS likelihood for correlated noise is not yet implemented." ) # return "gls" + toa_type = "wb" if self.is_wideband else "nb" + + return f"{method}-{toa_type}" + def lnprior(self, params): """Basic implementation of a factorized log prior. More complex priors must be separately implemented. @@ -165,9 +172,11 @@ def lnlikelihood(self, params): Returns: float: The value of the log-likelihood at params """ - if self.likelihood_method == "wls": - return self._wls_lnlikelihood(params) - elif self.likelihood_method == "gls": + if self.likelihood_method == "wls-nb": + return self._wls_nb_lnlikelihood(params) + elif self.likelihood_method == "wls-wb": + return self._wls_wb_lnlikelihood(params) + elif self.likelihood_method in ["gls-nb", "gls-wb"]: raise NotImplementedError( "GLS likelihood for correlated noise is not yet implemented." ) @@ -190,10 +199,11 @@ def lnposterior(self, params): else: return lnpr + self.lnlikelihood(params) - def _wls_lnlikelihood(self, params): - """Implementation of Log-Likelihood function for uncorrelated noise only. - `wls' stands for weighted least squares. Also includes the normalization - term to enable sampling over white noise parameters (EFAC and EQUAD). + def _wls_nb_lnlikelihood(self, params): + """Implementation of Log-Likelihood function for uncorrelated noise only for + narrow-band TOAs. `wls' stands for weighted least squares. Also includes the + normalization term to enable sampling over white noise parameters (EFAC and + EQUAD). Args: params (array-like): Parameters @@ -207,3 +217,30 @@ def _wls_lnlikelihood(self, params): chi2 = res.calc_chi2() sigmas = self.model.scaled_toa_uncertainty(self.toas).si.value return -chi2 / 2 - np.sum(np.log(sigmas)) + + def _wls_wb_lnlikelihood(self, params): + """Implementation of Log-Likelihood function for uncorrelated noise only for + wide-band TOAs. `wls' stands for weighted least squares. Also includes the + normalization terms to enable sampling over white noise parameters (EFAC, EQUAD, + DMEFAC and DMEQUAD). + + Args: + params (array-like): Parameters + + Returns: + float: The value of the log-likelihood at params + """ + params_dict = dict(zip(self.param_labels, params)) + self.model.set_param_values(params_dict) + + res = WidebandTOAResiduals(self.toas, self.model, track_mode=self.track_mode) + + chi2_toa = res.toa.calc_chi2() + sigmas_toa = self.model.scaled_toa_uncertainty(self.toas).si.value + lnL_toa = -chi2_toa / 2 - np.sum(np.log(sigmas_toa)) + + chi2_dm = res.dm.calc_chi2() + sigmas_dm = self.model.scaled_dm_uncertainty(self.toas).si.value + lnL_dm = -chi2_dm / 2 - np.sum(np.log(sigmas_dm)) + + return lnL_toa + lnL_dm From 7151559049224006a4c02d2b2896996ca4032efd Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 16:57:07 -0500 Subject: [PATCH 02/30] fix test_bayesian.py --- tests/test_bayesian.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_bayesian.py b/tests/test_bayesian.py index 2b127a84b..2499b017c 100644 --- a/tests/test_bayesian.py +++ b/tests/test_bayesian.py @@ -68,7 +68,7 @@ def test_no_noise(data_NGC6440E): bt = BayesianTiming(model, toas) maxlike_params = np.array([param.value for param in bt.params], dtype=float) lnl = bt.lnlikelihood(maxlike_params) - assert bt.likelihood_method == "wls" and not np.isnan(lnl) + assert bt.likelihood_method == "wls-nb" and not np.isnan(lnl) def test_white_noise(data_NGC6440E_efac): @@ -76,7 +76,7 @@ def test_white_noise(data_NGC6440E_efac): bt = BayesianTiming(model, toas) maxlike_params = np.array([param.value for param in bt.params], dtype=float) lnl = bt.lnlikelihood(maxlike_params) - assert bt.likelihood_method == "wls" and not np.isnan(lnl) + assert bt.likelihood_method == "wls-nb" and not np.isnan(lnl) def test_lnlikelihood_unit_efac(data_NGC6440E, data_NGC6440E_efac): @@ -138,7 +138,7 @@ def test_prior_dict(data_NGC6440E_efac): assert np.isfinite(lnpr) -def test_wideband_exception(data_J0740p6620_wb): - model, toas = data_J0740p6620_wb - with pytest.raises(NotImplementedError): - bt = BayesianTiming(model, toas) +# def test_wideband_exception(data_J0740p6620_wb): +# model, toas = data_J0740p6620_wb +# with pytest.raises(NotImplementedError): +# bt = BayesianTiming(model, toas) From 8b769594d4347e3b6ff0677730cfa4a0f4fe530f Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 17:01:46 -0500 Subject: [PATCH 03/30] cleanup --- src/pint/bayesian.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 3ba070813..96576dedc 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -56,9 +56,6 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self.is_wideband = toas.is_wideband() - # if toas.is_wideband(): - # raise NotImplementedError("Wideband TOAs are not yet supported.") - self.param_labels = self.model.free_params self.params = [getattr(self.model, par) for par in self.param_labels] self.nparams = len(self.param_labels) From a767d929f47ab985d078ae44833531d83a381cdc Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 17:11:23 -0500 Subject: [PATCH 04/30] test_wideband_data --- tests/test_bayesian.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/tests/test_bayesian.py b/tests/test_bayesian.py index 2499b017c..4998869b2 100644 --- a/tests/test_bayesian.py +++ b/tests/test_bayesian.py @@ -138,7 +138,21 @@ def test_prior_dict(data_NGC6440E_efac): assert np.isfinite(lnpr) -# def test_wideband_exception(data_J0740p6620_wb): -# model, toas = data_J0740p6620_wb -# with pytest.raises(NotImplementedError): -# bt = BayesianTiming(model, toas) +def test_wideband_data(data_J0740p6620_wb): + model, toas = data_J0740p6620_wb + bt = BayesianTiming(model, toas) + + assert bt.is_wideband and bt.likelihood_method == "wls-wb" + + test_cube = 0.5 * np.ones(bt.nparams) + test_params = bt.prior_transform(test_cube) + assert np.all(np.isfinite(test_params)) + + lnpr = bt.lnprior(test_params) + assert np.isfinite(lnpr) + + lnl = bt.lnlikelihood(test_params) + assert np.isfinite(lnl) + + lnp = bt.lnposterior(test_params) + assert np.isfinite(lnp) and np.isclose(lnp, lnpr + lnl) From 53d8aa74d739bf4efa7c7fdbc229d45fab1bf9bd Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 17:13:36 -0500 Subject: [PATCH 05/30] WidebandTOAResiduals doesn't support track_mode for some reason. --- src/pint/bayesian.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 96576dedc..405a4990f 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -230,7 +230,11 @@ def _wls_wb_lnlikelihood(self, params): params_dict = dict(zip(self.param_labels, params)) self.model.set_param_values(params_dict) - res = WidebandTOAResiduals(self.toas, self.model, track_mode=self.track_mode) + res = WidebandTOAResiduals( + self.toas, + self.model, + # track_mode=self.track_mode + ) chi2_toa = res.toa.calc_chi2() sigmas_toa = self.model.scaled_toa_uncertainty(self.toas).si.value From 61235091f4f67633f87732d4128045014ae908a0 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 19 Oct 2022 17:53:34 -0500 Subject: [PATCH 06/30] ok, track_mode works now. --- src/pint/bayesian.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 405a4990f..d931df3d7 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -231,9 +231,7 @@ def _wls_wb_lnlikelihood(self, params): self.model.set_param_values(params_dict) res = WidebandTOAResiduals( - self.toas, - self.model, - # track_mode=self.track_mode + self.toas, self.model, toa_resid_args={"track_mode": self.track_mode} ) chi2_toa = res.toa.calc_chi2() From 89032192fd04ef1f3923702ead5e23b2da5a4f1d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 20 Oct 2022 15:28:25 -0500 Subject: [PATCH 07/30] Wideband Bayesian example --- .../bayesian-example-J1614-2230-wb.py | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 docs/examples/bayesian-example-J1614-2230-wb.py diff --git a/docs/examples/bayesian-example-J1614-2230-wb.py b/docs/examples/bayesian-example-J1614-2230-wb.py new file mode 100644 index 000000000..e31261446 --- /dev/null +++ b/docs/examples/bayesian-example-J1614-2230-wb.py @@ -0,0 +1,69 @@ +from pint.models import get_model_and_toas +from pint.bayesian import BayesianTiming +from pint.models.priors import Prior +from pint.config import examplefile +from scipy.stats import uniform + +import numpy as np +import matplotlib.pyplot as plt +import nestle +import corner + +# Wideband par and tim files +parfile = examplefile("J1614-2230_NANOGrav_12yv3.wb.gls.par") +timfile = examplefile("J1614-2230_NANOGrav_12yv3.wb.tim") +model, toas = get_model_and_toas(parfile, timfile) + +# The par files have a lot of free parameters. +# Let us keep only a few of them free so that the sampling finishes in a reasonable time. +for param_label, param in model.get_params_dict().items(): + getattr(model, param_label).frozen = True + +model.ELAT.frozen = False +model.ELONG.frozen = False +for par in model.free_params: + param = getattr(model, par) + param_min = float(param.value - 3 * param.uncertainty_value) + param_span = float(6 * param.uncertainty_value) + param.prior = Prior(uniform(param_min, param_span)) + +model.DM.frozen = False +dm = float(model.DM.value) +model.DM.prior = Prior(uniform(dm - 0.05, dm + 0.05)) + +model.EFAC1.frozen = False +model.EFAC1.prior = Prior(uniform(0.1, 2)) + +model.DMEFAC1.frozen = False +model.DMEFAC1.prior = Prior(uniform(0.1, 5)) + +# Make sure the free parameters are what we expect. +print(model.free_params) + +# Create the BayesianTiming object +bt = BayesianTiming(model, toas) + +# Test out if the likelihood function works. +test_cube = 0.5 * np.ones(bt.nparams) +test_params = bt.prior_transform(test_cube) +test_lnl = bt.lnlikelihood(test_params) + +# Do the sampling using nestle +result_nestle = nestle.sample( + bt.lnlikelihood, + bt.prior_transform, + bt.nparams, + method="multi", + npoints=200, + dlogz=0.5, + callback=nestle.print_progress, +) + +# Visualize the posterior distribution +fig = corner.corner( + result_nestle.samples, + weights=result_nestle.weights, + labels=bt.param_labels, + range=[0.999] * bt.nparams, +) +plt.show() From a45550d59faa5bd2906f43b9c1f7a3bf9352eb38 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 20 Oct 2022 15:29:29 -0500 Subject: [PATCH 08/30] comment about run time --- docs/examples/bayesian-example-J1614-2230-wb.py | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/examples/bayesian-example-J1614-2230-wb.py b/docs/examples/bayesian-example-J1614-2230-wb.py index e31261446..91f482f87 100644 --- a/docs/examples/bayesian-example-J1614-2230-wb.py +++ b/docs/examples/bayesian-example-J1614-2230-wb.py @@ -49,6 +49,7 @@ test_lnl = bt.lnlikelihood(test_params) # Do the sampling using nestle +# This takes about over 30 minutes to finish. result_nestle = nestle.sample( bt.lnlikelihood, bt.prior_transform, From e9f3bd4addd488bd1c1bc1e0d6d02de4e7c4943c Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 20 Oct 2022 15:32:35 -0500 Subject: [PATCH 09/30] print statement --- docs/examples/bayesian-example-J1614-2230-wb.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/examples/bayesian-example-J1614-2230-wb.py b/docs/examples/bayesian-example-J1614-2230-wb.py index 91f482f87..35eb5492e 100644 --- a/docs/examples/bayesian-example-J1614-2230-wb.py +++ b/docs/examples/bayesian-example-J1614-2230-wb.py @@ -14,7 +14,7 @@ timfile = examplefile("J1614-2230_NANOGrav_12yv3.wb.tim") model, toas = get_model_and_toas(parfile, timfile) -# The par files have a lot of free parameters. +# The par file has a lot of free parameters. # Let us keep only a few of them free so that the sampling finishes in a reasonable time. for param_label, param in model.get_params_dict().items(): getattr(model, param_label).frozen = True @@ -47,6 +47,7 @@ test_cube = 0.5 * np.ones(bt.nparams) test_params = bt.prior_transform(test_cube) test_lnl = bt.lnlikelihood(test_params) +print(test_lnl) # Do the sampling using nestle # This takes about over 30 minutes to finish. From caf8d3288523dd4a9ba6c8a233d433b5727d0275 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 20 Oct 2022 16:37:26 -0500 Subject: [PATCH 10/30] docstring formatting fix --- src/pint/bayesian.py | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index d931df3d7..8540438ae 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -120,10 +120,12 @@ def lnprior(self, params): More complex priors must be separately implemented. Args: - params (array-like): Parameters + params : (array-like) + Parameters Returns: - float: Value of the log-prior at params + float : + Value of the log-prior at params """ if len(params) != self.nparams: raise IndexError( @@ -145,11 +147,12 @@ def prior_transform(self, cube): More complex prior transforms must be separately implemented. Args: - cube (array-like): Sample drawn from a uniform distribution defined in an - nparams-dimensional unit hypercube. + cube : (array-like) + Sample drawn from a uniform distribution defined in an nparams-dimensional unit hypercube. Returns: - ndarray : Sample drawn from the prior distribution + ndarray : + Sample drawn from the prior distribution """ result = np.array( [param.prior._rv.ppf(x) for x, param in zip(cube, self.params)] @@ -164,10 +167,12 @@ def lnlikelihood(self, params): is the generalized least-squares metric. For reference, see, e.g., Lentati+ 2013. Args: - params (array-like): Parameters + params : (array-like) + Parameters Returns: - float: The value of the log-likelihood at params + float : + The value of the log-likelihood at params """ if self.likelihood_method == "wls-nb": return self._wls_nb_lnlikelihood(params) @@ -185,10 +190,12 @@ def lnposterior(self, params): is not evaluated. Args: - params (array-like): Parameters + params : (array-like) + Parameters Returns: - float: The value of the log-posterior at params + float : + The value of the log-posterior at params """ lnpr = self.lnprior(params) if np.isnan(lnpr): @@ -203,10 +210,12 @@ def _wls_nb_lnlikelihood(self, params): EQUAD). Args: - params (array-like): Parameters + params : (array-like) + Parameters Returns: - float: The value of the log-likelihood at params + float : + The value of the log-likelihood at params """ params_dict = dict(zip(self.param_labels, params)) self.model.set_param_values(params_dict) @@ -222,10 +231,12 @@ def _wls_wb_lnlikelihood(self, params): DMEFAC and DMEQUAD). Args: - params (array-like): Parameters + params : (array-like) + Parameters Returns: - float: The value of the log-likelihood at params + float : + The value of the log-likelihood at params """ params_dict = dict(zip(self.param_labels, params)) self.model.set_param_values(params_dict) From b4eef52ad7bc2a761f7ffdffbd06a65ba4986a9e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 20 Oct 2022 17:04:38 -0500 Subject: [PATCH 11/30] _decide_likelihood_method should directly set the likelihood method. --- src/pint/bayesian.py | 43 ++++++++++++++----------------------------- 1 file changed, 14 insertions(+), 29 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 8540438ae..a72ba9a2c 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -76,7 +76,7 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self._validate_priors() - self.likelihood_method = self._decide_likelihood_method() + self._decide_likelihood_method() self.track_mode = "use_pulse_numbers" if use_pulse_numbers else "nearest" @@ -94,26 +94,20 @@ def _decide_likelihood_method(self): quares with normalization term (gls), for narrow-band (nb) or wide-band (wb) dataset.""" - if "NoiseComponent" not in self.model.component_types: - method = "wls" + if ( + "NoiseComponent" not in self.model.component_types + or not self.model.has_correlated_errors + ): + if self.is_wideband: + self.likelihood_method = "wls-wb" + self._lnlikelihood = self._wls_wb_lnlikelihood + else: + self.likelihood_method = "wls-nb" + self._lnlikelihood = self._wls_nb_lnlikelihood else: - correlated_errors_present = np.any( - [ - nc.introduces_correlated_errors - for nc in self.model.NoiseComponent_list - ] + raise NotImplementedError( + "GLS likelihood for correlated noise is not yet implemented." ) - if not correlated_errors_present: - method = "wls" - else: - raise NotImplementedError( - "GLS likelihood for correlated noise is not yet implemented." - ) - # return "gls" - - toa_type = "wb" if self.is_wideband else "nb" - - return f"{method}-{toa_type}" def lnprior(self, params): """Basic implementation of a factorized log prior. @@ -174,16 +168,7 @@ def lnlikelihood(self, params): float : The value of the log-likelihood at params """ - if self.likelihood_method == "wls-nb": - return self._wls_nb_lnlikelihood(params) - elif self.likelihood_method == "wls-wb": - return self._wls_wb_lnlikelihood(params) - elif self.likelihood_method in ["gls-nb", "gls-wb"]: - raise NotImplementedError( - "GLS likelihood for correlated noise is not yet implemented." - ) - else: - raise ValueError(f"Unknown likelihood method '{self.likelihood_method}'.") + return self._lnlikelihood(params) def lnposterior(self, params): """Log-posterior function. If the prior evaluates to zero, the likelihood From 6995d42f1ba2a7a1db1b2b4fae7f3c5b1759a2b7 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 2 Dec 2022 14:28:26 -0600 Subject: [PATCH 12/30] cleanup --- src/pint/bayesian.py | 32 +++++++++++++------------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index a72ba9a2c..2ec755079 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -90,25 +90,25 @@ def _validate_priors(self): ) def _decide_likelihood_method(self): - """Weighted least squares with normalization term (wls), or Generalized leasts - quares with normalization term (gls), for narrow-band (nb) or wide-band (wb) + """Weighted least squares with normalization term (wls), or Generalized least + squares with normalization term (gls), for narrow-band (nb) or wide-band (wb) dataset.""" if ( - "NoiseComponent" not in self.model.component_types - or not self.model.has_correlated_errors + "NoiseComponent" in self.model.component_types + and self.model.has_correlated_errors ): - if self.is_wideband: - self.likelihood_method = "wls-wb" - self._lnlikelihood = self._wls_wb_lnlikelihood - else: - self.likelihood_method = "wls-nb" - self._lnlikelihood = self._wls_nb_lnlikelihood - else: raise NotImplementedError( "GLS likelihood for correlated noise is not yet implemented." ) + if self.is_wideband: + self.likelihood_method = "wls-wb" + self._lnlikelihood = self._wls_wb_lnlikelihood + else: + self.likelihood_method = "wls-nb" + self._lnlikelihood = self._wls_nb_lnlikelihood + def lnprior(self, params): """Basic implementation of a factorized log prior. More complex priors must be separately implemented. @@ -148,10 +148,7 @@ def prior_transform(self, cube): ndarray : Sample drawn from the prior distribution """ - result = np.array( - [param.prior._rv.ppf(x) for x, param in zip(cube, self.params)] - ) - return result + return np.array([param.prior._rv.ppf(x) for x, param in zip(cube, self.params)]) def lnlikelihood(self, params): """The Log-likelihood function. If the model does not contain any noise components or @@ -183,10 +180,7 @@ def lnposterior(self, params): The value of the log-posterior at params """ lnpr = self.lnprior(params) - if np.isnan(lnpr): - return -np.inf - else: - return lnpr + self.lnlikelihood(params) + return -np.inf if np.isnan(lnpr) else lnpr + self.lnlikelihood(params) def _wls_nb_lnlikelihood(self, params): """Implementation of Log-Likelihood function for uncorrelated noise only for From e9fb87efdaca3edd9c438172e40727f578005db7 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 9 Dec 2022 17:14:23 -0600 Subject: [PATCH 13/30] fix doc build? --- docs/conf.py | 1 + src/pint/utils.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 19997984f..3eaf67578 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -119,6 +119,7 @@ "examples/compare_tempo2_B1855.py", "examples/example_dmx_ranges.py", "examples/example_pulse_numbers.py", + "examples/bayesian-example-NGC6440E.py", "conf.py", "_ext", ] diff --git a/src/pint/utils.py b/src/pint/utils.py index 885dc118b..ff698e8e5 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1758,8 +1758,8 @@ def compute_hash(filename): cryptographically robust. It uses the SHA256 algorithm, which is known to be vulnerable to a length-extension attack. - Parameter - --------- + Parameters + ---------- f : str or Path or file-like The source of input. If file-like, it should return ``bytes`` not ``str`` - that is, the file should be opened in binary mode. From 314c53ef641ad1829a006679409fd5250504c93d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 13 Feb 2023 15:55:10 -0600 Subject: [PATCH 14/30] Update conf.py --- docs/conf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/conf.py b/docs/conf.py index 399b82c31..f0061ff76 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -120,6 +120,7 @@ "examples/example_dmx_ranges.py", "examples/example_pulse_numbers.py", "examples/bayesian-example-NGC6440E.py", + "examples/bayesian-example-J1614-2230-wb.py", "conf.py", "_ext", ] From 22716c49b9cabc2a953edf55e9fe2ad0901c9685 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 09:51:45 -0500 Subject: [PATCH 15/30] fix broken merge and tests --- src/pint/bayesian.py | 41 ++++++++++++++++++++++------------------- tests/test_bayesian.py | 6 +++--- 2 files changed, 25 insertions(+), 22 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 8bb484b4b..96f86cfa2 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -76,7 +76,7 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self._validate_priors() - self._decide_likelihood_method() + self.likelihood_method = self._decide_likelihood_method() self.track_mode = "use_pulse_numbers" if use_pulse_numbers else "nearest" @@ -104,19 +104,16 @@ def _decide_likelihood_method(self): ) else: return "wls" - # return "gls" def lnprior(self, params): """Basic implementation of a factorized log prior. More complex priors must be separately implemented. Args: - params : (array-like) - Parameters + params (array-like): Parameters Returns: - float : - Value of the log-prior at params + float: Value of the log-prior at params """ if len(params) != self.nparams: raise IndexError( @@ -138,12 +135,11 @@ def prior_transform(self, cube): More complex prior transforms must be separately implemented. Args: - cube : (array-like) - Sample drawn from a uniform distribution defined in an nparams-dimensional unit hypercube. + cube (array-like): Sample drawn from a uniform distribution defined in an + nparams-dimensional unit hypercube. Returns: - ndarray : - Sample drawn from the prior distribution + ndarray : Sample drawn from the prior distribution """ return np.array([param.prior._rv.ppf(x) for x, param in zip(cube, self.params)]) @@ -155,26 +151,33 @@ def lnlikelihood(self, params): is the generalized least-squares metric. For reference, see, e.g., Lentati+ 2013. Args: - params : (array-like) - Parameters + params (array-like): Parameters Returns: - float : - The value of the log-likelihood at params + float: The value of the log-likelihood at params """ - return self._lnlikelihood(params) + if self.likelihood_method == "wls": + return ( + self._wls_wb_lnlikelihood(params) + if self.is_wideband + else self._wls_nb_lnlikelihood(params) + ) + elif self.likelihood_method == "gls": + raise NotImplementedError( + "GLS likelihood for correlated noise is not yet implemented." + ) + else: + raise ValueError(f"Unknown likelihood method '{self.likelihood_method}'.") def lnposterior(self, params): """Log-posterior function. If the prior evaluates to zero, the likelihood is not evaluated. Args: - params : (array-like) - Parameters + params (array-like): Parameters Returns: - float : - The value of the log-posterior at params + float: The value of the log-posterior at params """ lnpr = self.lnprior(params) return lnpr + self.lnlikelihood(params) if np.isfinite(lnpr) else -np.inf diff --git a/tests/test_bayesian.py b/tests/test_bayesian.py index fdeff5e7c..4487984ed 100644 --- a/tests/test_bayesian.py +++ b/tests/test_bayesian.py @@ -94,7 +94,7 @@ def test_no_noise(data_NGC6440E): bt = BayesianTiming(model, toas) maxlike_params = np.array([param.value for param in bt.params], dtype=float) lnl = bt.lnlikelihood(maxlike_params) - assert bt.likelihood_method == "wls-nb" and not np.isnan(lnl) + assert bt.likelihood_method == "wls" and not np.isnan(lnl) def test_white_noise(data_NGC6440E_efac): @@ -102,7 +102,7 @@ def test_white_noise(data_NGC6440E_efac): bt = BayesianTiming(model, toas) maxlike_params = np.array([param.value for param in bt.params], dtype=float) lnl = bt.lnlikelihood(maxlike_params) - assert bt.likelihood_method == "wls-nb" and not np.isnan(lnl) + assert bt.likelihood_method == "wls" and not np.isnan(lnl) def test_lnlikelihood_unit_efac(data_NGC6440E, data_NGC6440E_efac): @@ -182,7 +182,7 @@ def test_wideband_data(data_J0740p6620_wb): model, toas = data_J0740p6620_wb bt = BayesianTiming(model, toas) - assert bt.is_wideband and bt.likelihood_method == "wls-wb" + assert bt.is_wideband and bt.likelihood_method == "wls" test_cube = 0.5 * np.ones(bt.nparams) test_params = bt.prior_transform(test_cube) From 24fa3b084245b17d92cc02aef02b0f1fa9b8e220 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 12:57:22 -0500 Subject: [PATCH 16/30] compute_pulse_numbers --- src/pint/bayesian.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 96f86cfa2..7a1b9e9a3 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -54,6 +54,11 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self.model = deepcopy(model) self.toas = toas + if use_pulse_numbers: + self.toas.compute_pulse_numbers(self.model) + + self.track_mode = "use_pulse_numbers" if use_pulse_numbers else "nearest" + self.is_wideband = toas.is_wideband() self.param_labels = self.model.free_params @@ -78,8 +83,6 @@ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): self.likelihood_method = self._decide_likelihood_method() - self.track_mode = "use_pulse_numbers" if use_pulse_numbers else "nearest" - def _validate_priors(self): for param in self.params: if not hasattr(param, "prior") or param.prior is None: From 0a041c7adcdc5998fa57aa5340c63a3dbc89233e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 14:35:05 -0500 Subject: [PATCH 17/30] example --- docs/examples/bayesian-wideband-example.py | 71 +++++++ src/pint/bayesian.py | 3 +- src/pint/data/examples/test-wb-0.par | 31 +++ src/pint/data/examples/test-wb-0.tim | 207 +++++++++++++++++++++ 4 files changed, 311 insertions(+), 1 deletion(-) create mode 100644 docs/examples/bayesian-wideband-example.py create mode 100644 src/pint/data/examples/test-wb-0.par create mode 100644 src/pint/data/examples/test-wb-0.tim diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py new file mode 100644 index 000000000..c521381ac --- /dev/null +++ b/docs/examples/bayesian-wideband-example.py @@ -0,0 +1,71 @@ +from pint.models import get_model_and_toas +from pint.fitter import WidebandDownhillFitter +from pint.models.priors import Prior +from scipy.stats import uniform +from pint.bayesian import BayesianTiming +import emcee +import numpy as np +from pint.logging import setup as setup_log +import corner +import matplotlib.pyplot as plt +from pint.config import examplefile + +setup_log(level="WARNING") + +m, t = get_model_and_toas(examplefile("test-wb-0.par"), examplefile("test-wb-0.tim")) + +ftr = WidebandDownhillFitter(t, m) +ftr.fit_toas() +m = ftr.model + +prior_info = {} +for par in m.free_params: + param = getattr(m, par) + param_min = float(param.value - 10 * param.uncertainty_value) + param_max = float(param.value + 10 * param.uncertainty_value) + prior_info[par] = {"distr": "uniform", "pmin": param_min, "pmax": param_max} + +prior_info["EFAC1"] = {"distr": "normal", "mu": 1, "sigma": 0.1} +prior_info["DMEFAC1"] = {"distr": "normal", "mu": 1, "sigma": 0.1} + +m.EFAC1.frozen = False +m.EFAC1.uncertainty_value = 0.01 +m.DMEFAC1.frozen = False +m.DMEFAC1.uncertainty_value = 0.01 + +bt = BayesianTiming(m, t, use_pulse_numbers=True, prior_info=prior_info) + +print("Number of parameters = ", bt.nparams) +print("Likelihood method = ", bt.likelihood_method) + +nwalkers = 25 +sampler = emcee.EnsembleSampler(nwalkers, bt.nparams, bt.lnposterior) + +maxlike_params = np.array([param.value for param in bt.params], dtype=float) +maxlike_errors = [param.uncertainty_value for param in bt.params] +start_points = ( + np.repeat([maxlike_params], nwalkers).reshape(bt.nparams, nwalkers).T + + np.random.randn(nwalkers, bt.nparams) * maxlike_errors +) + +print("Running emcee...") +chain_length = 1000 +sampler.run_mcmc( + start_points, + chain_length, + progress=True, +) + +samples_emcee = sampler.get_chain(flat=True, discard=100) + +for idx, param_chain in enumerate(samples_emcee.T): + plt.subplot(bt.nparams, 1, idx + 1) + plt.plot(param_chain) + plt.ylabel(bt.param_labels[idx]) + plt.autoscale() +plt.show() + +fig = corner.corner( + samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params +) +plt.show() diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 7a1b9e9a3..48e0c2d97 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -120,7 +120,8 @@ def lnprior(self, params): """ if len(params) != self.nparams: raise IndexError( - f"The number of input parameters ({len(params)}) should be the same as the number of free parameters ({self.nparams})." + f"The number of input parameters ({len(params)}) should be the same " + f"as the number of free parameters ({self.nparams})." ) lnsum = 0.0 diff --git a/src/pint/data/examples/test-wb-0.par b/src/pint/data/examples/test-wb-0.par new file mode 100644 index 000000000..568e9b5b1 --- /dev/null +++ b/src/pint/data/examples/test-wb-0.par @@ -0,0 +1,31 @@ +# Created: 2023-08-18T12:08:17.544471 +# PINT_version: 0.9.6+119.g22716c49 +# User: Abhimanyu Susobhanan (abhimanyu) +# Host: abhimanyu-VirtualBox +# OS: Linux-6.2.0-26-generic-x86_64-with-glibc2.35 +# Python: 3.10.12 (main, Jul 5 2023, 18:54:27) [GCC 11.2.0] +# Format: pint +PSRJ WBTEST +EPHEM DE440 +DILATEFREQ N +DMDATA N +NTOA 0 +CHI2 0.0 +RAJ 4:00:00.00000000 1 0.00000000000000000000 +DECJ 15:00:00.00000000 1 0.00000000000000000000 +PMRA 0.0 +PMDEC 0.0 +PX 0.0 +POSEPOCH 55000.0000000000000000 +F0 100.0 1 0.0 +F1 1e-15 1 0.0 +PEPOCH 55000.0000000000000000 +PLANET_SHAPIRO N +DM 10.0 1 0.0 +DM1 0.0001 1 0.0 +DMEPOCH 55000.0000000000000000 +TZRMJD 55005.0251256281408097 +TZRSITE ssb +TZRFRQ inf +EFAC tel gmrt 1 +DMEFAC tel gmrt 1 \ No newline at end of file diff --git a/src/pint/data/examples/test-wb-0.tim b/src/pint/data/examples/test-wb-0.tim new file mode 100644 index 000000000..a7d3d52c1 --- /dev/null +++ b/src/pint/data/examples/test-wb-0.tim @@ -0,0 +1,207 @@ +FORMAT 1 +C Created: 2023-08-18T12:08:45.036835 +C PINT_version: 0.9.6+119.g22716c49 +C User: Abhimanyu Susobhanan (abhimanyu) +C Host: abhimanyu-VirtualBox +C OS: Linux-6.2.0-26-generic-x86_64-with-glibc2.35 +C Python: 3.10.12 (main, Jul 5 2023, 18:54:27) [GCC 11.2.0] +fake 400.000000 54000.0000000338171065 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999689348312622012 -pn -8683387108.0 +fake 800.000000 54010.0502512459059028 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999817394902294414 -pn -8596545796.0 +fake 1200.000000 54020.1005025362641782 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999891183572727027 -pn -8509705405.0 +fake 400.000000 54030.1507537448405324 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999631939252680306 -pn -8422866134.0 +fake 800.000000 54040.2010049845390741 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9995842407658120776 -pn -8336028086.0 +fake 1200.000000 54050.2512562274484143 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999814299821878541 -pn -8249191442.0 +fake 400.000000 54060.3015075897839931 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9996197994101480066 -pn -8162356296.0 +fake 800.000000 54070.3517587435774421 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999725184871209322 -pn -8075522630.0 +fake 1200.000000 54080.4020101060183797 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999767790770363772 -pn -7988690482.0 +fake 400.000000 54090.4522613409235186 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999686831665311229 -pn -7901859813.0 +fake 800.000000 54100.5025125235372685 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999841888231202267 -pn -7815030452.0 +fake 1200.000000 54110.5527637856061109 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999855240286045288 -pn -7728202298.0 +fake 400.000000 54120.6030150412243056 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999827872380426933 -pn -7641375182.0 +fake 800.000000 54130.6532663832630324 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999810559785563071 -pn -7554548810.0 +fake 1200.000000 54140.7035175807754050 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999944845353083878 -pn -7467722991.0 +fake 400.000000 54150.7537688665358102 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999636388273509543 -pn -7380897481.0 +fake 800.000000 54160.8040201505214583 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999806942922156322 -pn -7294071936.0 +fake 1200.000000 54170.8542713282116436 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999794849102732609 -pn -7207246151.0 +fake 400.000000 54180.9045226137458797 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999755115559108437 -pn -7120419888.0 +fake 800.000000 54190.9547738557419675 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999909233122702556 -pn -7033592836.0 +fake 1200.000000 54201.0050251104141204 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999546290756704471 -pn -6946764852.0 +fake 400.000000 54211.0552763242439815 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999844631460155727 -pn -6859935773.0 +fake 800.000000 54221.1055276941583796 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999672173798245581 -pn -6773105381.0 +fake 1200.000000 54231.1557789516169097 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999734640037218725 -pn -6686273645.0 +fake 400.000000 54241.2060301601370255 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999652091265824406 -pn -6599440510.0 +fake 800.000000 54251.2562814636363426 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999855249592944891 -pn -6512605883.0 +fake 1200.000000 54261.3065327059755903 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9996220393044169695 -pn -6425769852.0 +fake 400.000000 54271.3567838821766782 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999763661348765612 -pn -6338932478.0 +fake 800.000000 54281.4070351975145023 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999840281158982123 -pn -6252093789.0 +fake 1200.000000 54291.4572863853987268 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999763278947545926 -pn -6165253979.0 +fake 400.000000 54301.5075376360473032 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999901200977491657 -pn -6078413208.0 +fake 800.000000 54311.5577889072436342 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999719461962584742 -pn -5991571602.0 +fake 1200.000000 54321.6080401637861922 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999673262590097732 -pn -5904729427.0 +fake 400.000000 54331.6582915012324653 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999966387002557102 -pn -5817886914.0 +fake 800.000000 54341.7085426773096991 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999787029487756375 -pn -5731044244.0 +fake 1200.000000 54351.7587939677584142 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999932454784072644 -pn -5644201711.0 +fake 400.000000 54361.8090452441937731 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999987919231125403 -pn -5557359572.0 +fake 800.000000 54371.8592964550520139 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999755049704442288 -pn -5470518003.0 +fake 1200.000000 54381.9095477115447685 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999640428002742632 -pn -5383677276.0 +fake 400.000000 54391.9597989894186458 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998219673397774525 -pn -5296837611.0 +fake 800.000000 54402.0100502489624421 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999871983091432774 -pn -5209999116.0 +fake 1200.000000 54412.0603015108758565 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.00000703329270763 -pn -5123161982.0 +fake 400.000000 54422.1105527431482870 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999884527228331438 -pn -5036326330.0 +fake 800.000000 54432.1608040676221759 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999816662373855615 -pn -4949492142.0 +fake 1200.000000 54442.2110552311124421 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.99964302668641098 -pn -4862659484.0 +fake 400.000000 54452.2613065788043287 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998872578856580476 -pn -4775828332.0 +fake 800.000000 54462.3115577385225695 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999944859944271937 -pn -4688998523.0 +fake 1200.000000 54472.3618090102764815 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999879509734557669 -pn -4602169981.0 +fake 400.000000 54482.4120603123242361 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999995082747862556 -pn -4515342546.0 +fake 800.000000 54492.4623115789709374 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999816373048808085 -pn -4428515929.0 +fake 1200.000000 54502.5125628608401852 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999774228300115434 -pn -4341689957.0 +fake 400.000000 54512.5628140679162617 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9997085408988969114 -pn -4254864383.0 +fake 800.000000 54522.6130652806756481 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9999465770119952785 -pn -4168038864.0 +fake 1200.000000 54532.6633165776520833 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999716794891414485 -pn -4081213200.0 +fake 400.000000 54542.7135678561986690 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999874148552329637 -pn -3994387138.0 +fake 800.000000 54552.7638191202778934 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999771394414748284 -pn -3907560367.0 +fake 1200.000000 54562.8140703119680556 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999811388630607196 -pn -3820732736.0 +fake 400.000000 54572.8643216100541089 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998608980168683345 -pn -3733904058.0 +fake 800.000000 54582.9145728791742246 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999906682211421074 -pn -3647074118.0 +fake 1200.000000 54592.9648241367796759 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999895041223835885 -pn -3560242862.0 +fake 400.000000 54603.0150753694817014 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998716615977875505 -pn -3473410218.0 +fake 800.000000 54613.0653266786111111 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000010231048680057 -pn -3386576091.0 +fake 1200.000000 54623.1155778719024421 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000004776792867532 -pn -3299740544.0 +fake 400.000000 54633.1658291538335648 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.99967172207249051 -pn -3212903628.0 +fake 800.000000 54643.2160804013887037 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000108122754010403 -pn -3126065367.0 +fake 1200.000000 54653.2663316774607060 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999980552841559337 -pn -3039225929.0 +fake 400.000000 54663.3165828598465394 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000068670936557959 -pn -2952385478.0 +fake 800.000000 54673.3668341711507523 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999679371789044481 -pn -2865544126.0 +fake 1200.000000 54683.4170853739764931 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999990687399299406 -pn -2778702126.0 +fake 400.000000 54693.4673367320890972 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9999720094319615475 -pn -2691859715.0 +fake 800.000000 54703.5175879366010185 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.00005208609530464 -pn -2605017061.0 +fake 1200.000000 54713.5678392190929050 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999987574852184185 -pn -2518174457.0 +fake 400.000000 54723.6180904375442015 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.99989236262914892 -pn -2431332170.0 +fake 800.000000 54733.6683416922952431 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9997956601147193045 -pn -2344490365.0 +fake 1200.000000 54743.7185929525367246 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999866982796133246 -pn -2257649329.0 +fake 400.000000 54753.7688441728829630 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999859586492998137 -pn -2170809292.0 +fake 800.000000 54763.8190955198084374 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998250224034890485 -pn -2083970360.0 +fake 1200.000000 54773.8693467226388774 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0000284739417204605 -pn -1997132752.0 +fake 400.000000 54783.9195980059268866 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0002369037720472405 -pn -1910296595.0 +fake 800.000000 54793.9698491886990046 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000065942886162499 -pn -1823461885.0 +fake 1200.000000 54804.0201005511857755 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999944891336582828 -pn -1736628712.0 +fake 400.000000 54814.0703517722123264 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999927984918873445 -pn -1649797063.0 +fake 800.000000 54824.1206030679627083 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999787109292992464 -pn -1562966790.0 +fake 1200.000000 54834.1708543070402546 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000085109166870542 -pn -1476137742.0 +fake 400.000000 54844.2211055407940278 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999892553417434964 -pn -1389309959.0 +fake 800.000000 54854.2713567972613657 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000030425752240226 -pn -1302483073.0 +fake 1200.000000 54864.3216080450358565 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000126721528663949 -pn -1215656918.0 +fake 400.000000 54874.3718593442969676 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999967224511925204 -pn -1128831243.0 +fake 800.000000 54884.4221106018054630 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999973596921750957 -pn -1042005722.0 +fake 1200.000000 54894.4723617696319560 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999859830658437138 -pn -955180146.0 +fake 400.000000 54904.5226130309290625 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999938186615571616 -pn -868354255.0 +fake 800.000000 54914.5728643355897106 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9998099117997843585 -pn -781527742.0 +fake 1200.000000 54924.6231155570151967 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999933195919728759 -pn -694700435.0 +fake 400.000000 54934.6733667954781017 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999879310638236722 -pn -607872144.0 +fake 800.000000 54944.7236180665839236 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.00000920829822565 -pn -521042643.0 +fake 1200.000000 54954.7738693193937037 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999937597913816637 -pn -434211854.0 +fake 400.000000 54964.8241205821242014 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999776063534331207 -pn -347379703.0 +fake 800.000000 54974.8743718712418403 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000173321152803697 -pn -260546076.0 +fake 1200.000000 54984.9246230913911458 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000124638396050736 -pn -173711016.0 +fake 400.000000 54994.9748743358192478 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000005764303560627 -pn -86874574.0 +fake 800.000000 55005.0251256172320255 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999930822646455541 -pn -36749.0 +fake 1200.000000 55015.0753769322097917 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999909343956333098 -pn 86802299.0 +fake 400.000000 55025.1256281353382986 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999997018102579935 -pn 173642409.0 +fake 800.000000 55035.1758794046910995 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999960011834497916 -pn 260483490.0 +fake 1200.000000 55045.2261306872286574 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000095822229545901 -pn 347325288.0 +fake 400.000000 55055.2763818888063657 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999909720239615588 -pn 434167570.0 +fake 800.000000 55065.3266331769281134 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000019166610654427 -pn 521010185.0 +fake 1200.000000 55075.3768844260846644 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999921504496390928 -pn 607852826.0 +fake 400.000000 55085.4271357309605325 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9999044840334824925 -pn 694695237.0 +fake 800.000000 55095.4773869083597221 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999911289996185105 -pn 781537250.0 +fake 1200.000000 55105.5276381787532987 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000066185921840471 -pn 868378565.0 +fake 400.000000 55115.5778894215778819 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999946260456715445 -pn 955218957.0 +fake 800.000000 55125.6281406752153588 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.99994596279539534 -pn 1042058303.0 +fake 1200.000000 55135.6783919489981597 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000072974981106686 -pn 1128896371.0 +fake 400.000000 55145.7286431653392593 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000047762124557193 -pn 1215733029.0 +fake 800.000000 55155.7788944329273264 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9999371452968319925 -pn 1302568256.0 +fake 1200.000000 55165.8291457748755208 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999967032042070764 -pn 1389401948.0 +fake 400.000000 55175.8793970056304398 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000127460559598744 -pn 1476234108.0 +fake 800.000000 55185.9296482781195716 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999949353445512128 -pn 1563064854.0 +fake 1200.000000 55195.9798995514579283 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000038469599295368 -pn 1649894233.0 +fake 400.000000 55206.0301507896968171 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000184488525271191 -pn 1736722388.0 +fake 800.000000 55216.0804020083592593 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000020397349509981 -pn 1823549567.0 +fake 1200.000000 55226.1306532466018171 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999832160375013049 -pn 1910375940.0 +fake 400.000000 55236.1809045240011805 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000220882169742772 -pn 1997201742.0 +fake 800.000000 55246.2311557515021991 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0000712212124736955 -pn 2084027296.0 +fake 1200.000000 55256.2814070796484143 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.99988297755272494 -pn 2170852822.0 +fake 400.000000 55266.3316582414325116 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000027401418097545 -pn 2257678566.0 +fake 800.000000 55276.3819095680059607 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000059497362002379 -pn 2344504852.0 +fake 1200.000000 55286.4321607618225463 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000242234020140315 -pn 2431331861.0 +fake 400.000000 55296.4824120603521874 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000334406591437854 -pn 2518159783.0 +fake 800.000000 55306.5326633376824653 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999886245295953048 -pn 2604988867.0 +fake 1200.000000 55316.5829146155144906 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000059986473515738 -pn 2691819199.0 +fake 400.000000 55326.6331658011917592 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000066715297739709 -pn 2778650861.0 +fake 800.000000 55336.6834170336900000 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000088876729311907 -pn 2865483994.0 +fake 1200.000000 55346.7336683820114005 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0001555801232671445 -pn 2952318559.0 +fake 400.000000 55356.7839196332011459 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999993792149525078 -pn 3039154520.0 +fake 800.000000 55366.8341708451521991 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0000175345260493365 -pn 3125991898.0 +fake 1200.000000 55376.8844221036844907 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999927402643881202 -pn 3212830533.0 +fake 400.000000 55386.9346733113793634 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000126990980531342 -pn 3299670286.0 +fake 800.000000 55396.9849246155898148 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.9999658517649233766 -pn 3386511074.0 +fake 1200.000000 55407.0351759348490509 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999947721079837201 -pn 3473352642.0 +fake 400.000000 55417.0854271475039583 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000089478169504019 -pn 3560194776.0 +fake 800.000000 55427.1356784344858333 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000100905303029587 -pn 3647037321.0 +fake 1200.000000 55437.1859296868953009 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0001193380051912945 -pn 3733879972.0 +fake 400.000000 55447.2361809303057870 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000020504533308162 -pn 3820722483.0 +fake 800.000000 55457.2864321273571759 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000068139775406002 -pn 3907564673.0 +fake 1200.000000 55467.3366834195775347 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000128256706893354 -pn 3994406245.0 +fake 400.000000 55477.3869346802707523 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000173487726014024 -pn 4081246971.0 +fake 800.000000 55487.4371859654377315 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999922519302911486 -pn 4168086708.0 +fake 1200.000000 55497.4874372176980671 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000077837822342276 -pn 4254925225.0 +fake 400.000000 55507.5376884554523611 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000207002326607559 -pn 4341762373.0 +fake 800.000000 55517.5879397130096992 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000239330535125347 -pn 4428598110.0 +fake 1200.000000 55527.6381908988099653 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000162814231407506 -pn 4515432327.0 +fake 400.000000 55537.6884421805568172 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000120149866405688 -pn 4602265002.0 +fake 800.000000 55547.7386934888747339 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000237508928458069 -pn 4689096237.0 +fake 1200.000000 55557.7889446945634028 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000141186570696267 -pn 4775926069.0 +fake 400.000000 55567.8391959928024653 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0002246170119403965 -pn 4862754615.0 +fake 800.000000 55577.8894472530937152 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000219027074768849 -pn 4949582119.0 +fake 1200.000000 55587.9396984947352199 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000268042086678121 -pn 5036408740.0 +fake 400.000000 55597.9899497214191552 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000040338675025169 -pn 5123234696.0 +fake 800.000000 55608.0402009677436111 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000263667756289297 -pn 5210060320.0 +fake 1200.000000 55618.0904522558979745 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000257447578972783 -pn 5296885822.0 +fake 400.000000 55628.1407035277368403 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000263034605704899 -pn 5383711450.0 +fake 800.000000 55638.1909548204844676 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000185262083079626 -pn 5470537541.0 +fake 1200.000000 55648.2412060724572569 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000217004341520001 -pn 5557364273.0 +fake 400.000000 55658.2914572935849074 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000134242018416416 -pn 5644191850.0 +fake 800.000000 55668.3417085787174190 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000200539617497755 -pn 5731020538.0 +fake 1200.000000 55678.3919597921981829 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000234690861819577 -pn 5817850422.0 +fake 400.000000 55688.4422110333989121 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000117595413135451 -pn 5904681611.0 +fake 800.000000 55698.4924623667830324 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000242095632078261 -pn 5991514255.0 +fake 1200.000000 55708.5427135788948612 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000156527178478722 -pn 6078348320.0 +fake 400.000000 55718.5929648814338888 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000298488753902373 -pn 6165183800.0 +fake 800.000000 55728.6432161318090741 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.00020249294367174 -pn 6252020717.0 +fake 1200.000000 55738.6934673512211690 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0004272154687253536 -pn 6338858925.0 +fake 400.000000 55748.7437186170883681 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000058013579587481 -pn 6425698306.0 +fake 800.000000 55758.7939698085431249 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000378942298283491 -pn 6512538772.0 +fake 1200.000000 55768.8442211370403588 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000301078007239742 -pn 6599380089.0 +fake 400.000000 55778.8944723171171297 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000222538928678034 -pn 6686222049.0 +fake 800.000000 55788.9447236245216435 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000197672546348252 -pn 6773064493.0 +fake 1200.000000 55798.9949748948978704 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000208636614789991 -pn 6859907132.0 +fake 400.000000 55809.0452261702578357 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000151785477645361 -pn 6946749714.0 +fake 800.000000 55819.0954773940804167 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000344632926116159 -pn 7033592054.0 +fake 1200.000000 55829.1457285995028009 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0003719383921038395 -pn 7120433862.0 +fake 400.000000 55839.1959798529702083 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000278484473008056 -pn 7207274894.0 +fake 800.000000 55849.2462311071718403 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000126392667082668 -pn 7294115004.0 +fake 1200.000000 55859.2964824622073033 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000133314170880341 -pn 7380953957.0 +fake 400.000000 55869.3467336315645370 1.000 gmrt -pp_dme 0.0001 -pp_dm 9.999932149006183948 -pn 7467791577.0 +fake 800.000000 55879.3969849065807871 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000301782830795302 -pn 7554627822.0 +fake 1200.000000 55889.4472361936712269 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000220044521693582 -pn 7641462564.0 +fake 400.000000 55899.4974873990938426 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0000454418055698845 -pn 7728295755.0 +fake 800.000000 55909.5477386928395254 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000062125494163085 -pn 7815127495.0 +fake 1200.000000 55919.5979899230785996 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000270080431560062 -pn 7901957794.0 +fake 400.000000 55929.6482412396903356 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0002914779061594 -pn 7988786754.0 +fake 800.000000 55939.6984924249845602 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.00028655900684835 -pn 8075614613.0 +fake 1200.000000 55949.7487436704994676 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0003792051184662885 -pn 8162441507.0 +fake 400.000000 55959.7989949440956827 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000325634328277777 -pn 8249267654.0 +fake 800.000000 55969.8492462052852083 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000326622727975257 -pn 8336093382.0 +fake 1200.000000 55979.8994974572198611 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000160297267337141 -pn 8422918888.0 +fake 400.000000 55989.9497487285055439 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.000299306219806401 -pn 8509744435.0 +fake 800.000000 56000.0000000023089352 1.000 gmrt -pp_dme 0.0001 -pp_dm 10.0003150571158367504 -pn 8596570357.0 From d215e1fe79fad0ea16665b2a3b53097c5fd0dad1 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 14:36:25 -0500 Subject: [PATCH 18/30] cleanup --- docs/examples/bayesian-wideband-example.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index c521381ac..ddc02a0a8 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -1,14 +1,13 @@ -from pint.models import get_model_and_toas -from pint.fitter import WidebandDownhillFitter -from pint.models.priors import Prior -from scipy.stats import uniform -from pint.bayesian import BayesianTiming -import emcee -import numpy as np -from pint.logging import setup as setup_log import corner +import emcee import matplotlib.pyplot as plt +import numpy as np + +from pint.bayesian import BayesianTiming from pint.config import examplefile +from pint.fitter import WidebandDownhillFitter +from pint.logging import setup as setup_log +from pint.models import get_model_and_toas setup_log(level="WARNING") From bfc30b3508e860229cd9ce3a26e13e05c6cf3ec7 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 14:53:59 -0500 Subject: [PATCH 19/30] comments --- docs/examples/bayesian-wideband-example.py | 34 ++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index ddc02a0a8..1d2706e53 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -1,3 +1,4 @@ +# %% import corner import emcee import matplotlib.pyplot as plt @@ -9,14 +10,27 @@ from pint.logging import setup as setup_log from pint.models import get_model_and_toas +# %% +# Turn off log messages. They can slow down the processing. setup_log(level="WARNING") +# %% +# This is a simulated dataset. m, t = get_model_and_toas(examplefile("test-wb-0.par"), examplefile("test-wb-0.tim")) +# %% +# Fit the model to the data to get the parameter uncertainties. ftr = WidebandDownhillFitter(t, m) ftr.fit_toas() m = ftr.model +# %% +# Now set the priors. +# I am cheating here by setting the priors around the maximum likelihood estimates. +# This is a bad idea for real datasets and can bias the estimates. I am doing this +# here just to make everything finish faster. In the real world, these priors should +# be informed by, e.g. previous (independent) timing solutions, pulsar search results, +# VLBI localization etc. Note that unbounded uniform priors don't work here. prior_info = {} for par in m.free_params: param = getattr(m, par) @@ -24,6 +38,10 @@ param_max = float(param.value + 10 * param.uncertainty_value) prior_info[par] = {"distr": "uniform", "pmin": param_min, "pmax": param_max} +# %% +# Set the EFAC and DMEFAC priors and unfreeze them. +# Don't do this before the fitting step. The fitter doesn't know +# how to deal with noise parameters. prior_info["EFAC1"] = {"distr": "normal", "mu": 1, "sigma": 0.1} prior_info["DMEFAC1"] = {"distr": "normal", "mu": 1, "sigma": 0.1} @@ -32,14 +50,20 @@ m.DMEFAC1.frozen = False m.DMEFAC1.uncertainty_value = 0.01 +# %% +# The likelihood function behaves better if `use_pulse_numbers==True`. bt = BayesianTiming(m, t, use_pulse_numbers=True, prior_info=prior_info) +# %% print("Number of parameters = ", bt.nparams) print("Likelihood method = ", bt.likelihood_method) +# %% nwalkers = 25 sampler = emcee.EnsembleSampler(nwalkers, bt.nparams, bt.lnposterior) +# %% +# Start the sampler close to the maximul likelihood estimate. maxlike_params = np.array([param.value for param in bt.params], dtype=float) maxlike_errors = [param.uncertainty_value for param in bt.params] start_points = ( @@ -47,6 +71,7 @@ + np.random.randn(nwalkers, bt.nparams) * maxlike_errors ) +# %% print("Running emcee...") chain_length = 1000 sampler.run_mcmc( @@ -55,8 +80,11 @@ progress=True, ) +# %% samples_emcee = sampler.get_chain(flat=True, discard=100) +# %% +# Plot the chains to make sure they have converged and the burn-in has been removed properly. for idx, param_chain in enumerate(samples_emcee.T): plt.subplot(bt.nparams, 1, idx + 1) plt.plot(param_chain) @@ -64,7 +92,13 @@ plt.autoscale() plt.show() +# %% fig = corner.corner( samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params ) plt.show() + +# %% +# It seems like EFAC1 and DMEFAC1 have been overestimated. This is not an issue with the +# Bayesian analysis, but with the simulations. PINT simulations sometimes inject more or less +# noise than expected From e8c5bb33e2fd90c7ab63a876933923b18ffd9a11 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 15:02:11 -0500 Subject: [PATCH 20/30] comments --- src/pint/bayesian.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/pint/bayesian.py b/src/pint/bayesian.py index 48e0c2d97..c83a3aeca 100644 --- a/src/pint/bayesian.py +++ b/src/pint/bayesian.py @@ -34,19 +34,32 @@ class BayesianTiming: 1. The `prior` attribute of each free parameter in the `model` object should be set to an instance of :class:`pint.models.priors.Prior`. - 2. The parameters of BayesianTiming.model will change for every likelihood function call. + 2. The parameters of `BayesianTiming.model` will change for every likelihood function call. These parameters in general will not be the best-fit values. Hence, it is NOT a good idea to save it as a par file. - 3. Only narow-band TOAs are supported at present. + 3. Both narrow-band and wide-band TOAs are supported. 4. Currently, only uniform and normal distributions are supported in prior_info. More general priors should be set directly in the TimingModel object before creating the - BayesianTiming object. Here is an example prior_info object: - - `prior_info = { "F0" : {"distr" : "normal", "mu" : 1, "sigma" : 0.00001}, "EFAC1" : {"distr" : "uniform", "pmin" : 0.5, "pmax" : 2.0} }` - - See examples/bayesian-example-NGC6440E.py for detailed example. + BayesianTiming object. Here is an example prior_info object:: + + ``` + prior_info = { + "F0" : { + "distr" : "normal", + "mu" : 1, + "sigma" : 0.00001 + }, + "EFAC1" : { + "distr" : "uniform", + "pmin" : 0.5, + "pmax" : 2.0 + } + } + ``` + + See `examples/bayesian-example-NGC6440E.py` and `examples/bayesian-wideband-example` for detailed examples. """ def __init__(self, model, toas, use_pulse_numbers=False, prior_info=None): From f66772e1a30b76419e6e07a496292c13415b5371 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 15:07:28 -0500 Subject: [PATCH 21/30] cleanup --- docs/examples/bayesian-wideband-example.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index 1d2706e53..85f949eea 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -97,8 +97,3 @@ samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params ) plt.show() - -# %% -# It seems like EFAC1 and DMEFAC1 have been overestimated. This is not an issue with the -# Bayesian analysis, but with the simulations. PINT simulations sometimes inject more or less -# noise than expected From 6863a5d6d501d4da51519a56861664f52d7b1adb Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 15:19:18 -0500 Subject: [PATCH 22/30] delete example script --- .../bayesian-example-J1614-2230-wb.py | 71 ------------------- 1 file changed, 71 deletions(-) delete mode 100644 docs/examples/bayesian-example-J1614-2230-wb.py diff --git a/docs/examples/bayesian-example-J1614-2230-wb.py b/docs/examples/bayesian-example-J1614-2230-wb.py deleted file mode 100644 index 35eb5492e..000000000 --- a/docs/examples/bayesian-example-J1614-2230-wb.py +++ /dev/null @@ -1,71 +0,0 @@ -from pint.models import get_model_and_toas -from pint.bayesian import BayesianTiming -from pint.models.priors import Prior -from pint.config import examplefile -from scipy.stats import uniform - -import numpy as np -import matplotlib.pyplot as plt -import nestle -import corner - -# Wideband par and tim files -parfile = examplefile("J1614-2230_NANOGrav_12yv3.wb.gls.par") -timfile = examplefile("J1614-2230_NANOGrav_12yv3.wb.tim") -model, toas = get_model_and_toas(parfile, timfile) - -# The par file has a lot of free parameters. -# Let us keep only a few of them free so that the sampling finishes in a reasonable time. -for param_label, param in model.get_params_dict().items(): - getattr(model, param_label).frozen = True - -model.ELAT.frozen = False -model.ELONG.frozen = False -for par in model.free_params: - param = getattr(model, par) - param_min = float(param.value - 3 * param.uncertainty_value) - param_span = float(6 * param.uncertainty_value) - param.prior = Prior(uniform(param_min, param_span)) - -model.DM.frozen = False -dm = float(model.DM.value) -model.DM.prior = Prior(uniform(dm - 0.05, dm + 0.05)) - -model.EFAC1.frozen = False -model.EFAC1.prior = Prior(uniform(0.1, 2)) - -model.DMEFAC1.frozen = False -model.DMEFAC1.prior = Prior(uniform(0.1, 5)) - -# Make sure the free parameters are what we expect. -print(model.free_params) - -# Create the BayesianTiming object -bt = BayesianTiming(model, toas) - -# Test out if the likelihood function works. -test_cube = 0.5 * np.ones(bt.nparams) -test_params = bt.prior_transform(test_cube) -test_lnl = bt.lnlikelihood(test_params) -print(test_lnl) - -# Do the sampling using nestle -# This takes about over 30 minutes to finish. -result_nestle = nestle.sample( - bt.lnlikelihood, - bt.prior_transform, - bt.nparams, - method="multi", - npoints=200, - dlogz=0.5, - callback=nestle.print_progress, -) - -# Visualize the posterior distribution -fig = corner.corner( - result_nestle.samples, - weights=result_nestle.weights, - labels=bt.param_labels, - range=[0.999] * bt.nparams, -) -plt.show() From f1c55cc63f636d1fa288751e30e93e1119855883 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 15:24:43 -0500 Subject: [PATCH 23/30] update conf.py --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index f7da720bc..f72e3f1cf 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -120,7 +120,7 @@ "examples/example_dmx_ranges.py", "examples/example_pulse_numbers.py", "examples/bayesian-example-NGC6440E.py", - "examples/bayesian-example-J1614-2230-wb.py", + "examples/bayesian-wideband-example.py", "conf.py", "_ext", ] From af0171904498b82c966de4d01c600d422fb898f4 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 18 Aug 2023 16:24:27 -0500 Subject: [PATCH 24/30] disable logging --- docs/examples/bayesian-example-NGC6440E.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py index daf571e8f..f19afd219 100644 --- a/docs/examples/bayesian-example-NGC6440E.py +++ b/docs/examples/bayesian-example-NGC6440E.py @@ -6,6 +6,7 @@ from pint.bayesian import BayesianTiming from pint.config import examplefile from pint.models.priors import Prior +from pint.logging import setup as setup_log from scipy.stats import uniform # %% @@ -16,6 +17,10 @@ import io import matplotlib.pyplot as plt +# %% +# Turn off log messages. They can slow down the processing. +setup_log(level="WARNING") + # %% # Read the par and tim files parfile = examplefile("NGC6440E.par.good") From 5e486ed9e28e35b7c8bfcb27356b8669f9b3bac4 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 21 Aug 2023 17:28:25 -0500 Subject: [PATCH 25/30] changelog --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 2a67299aa..f7bfd18da 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -19,6 +19,7 @@ the released changes. - `pint.print_info()` function for bug reporting - Added an autocorrelation function to check for chain convergence in `event_optimize` - Minor doc updates to explain default NHARMS and missing derivative functions +- Support for wideband data in `pint.bayesian` (no correlated noise). ### Fixed - Deleting JUMP1 from flag tables will not prevent fitting - Simulating TOAs from tim file when PLANET_SHAPIRO is true now works From 18dedef00aa9d337f495cfd727db491bf29740ad Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 21 Sep 2023 11:18:26 -0500 Subject: [PATCH 26/30] examples in rtd --- docs/conf.py | 4 +- docs/examples/bayesian-example-NGC6440E.py | 127 ++++++++++++--------- docs/examples/bayesian-wideband-example.py | 44 ++++--- 3 files changed, 100 insertions(+), 75 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 4eee9a879..5516fbf3e 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -119,8 +119,8 @@ "examples/compare_tempo2_B1855.py", "examples/example_dmx_ranges.py", "examples/example_pulse_numbers.py", - "examples/bayesian-example-NGC6440E.py", - "examples/bayesian-wideband-example.py", + # "examples/bayesian-example-NGC6440E.py", + # "examples/bayesian-wideband-example.py", "conf.py", "_ext", ] diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py index f19afd219..d371c7316 100644 --- a/docs/examples/bayesian-example-NGC6440E.py +++ b/docs/examples/bayesian-example-NGC6440E.py @@ -75,31 +75,41 @@ + np.random.randn(nwalkers, bt.nparams) * maxlike_errors ) +# %% +# ** IMPORTANT!!! ** +# This is used to exclude some of the following time-consuming steps from the readthedocs build. +# Set this to False while actually using this example. +rtd = True + # %% # Use longer chain_length for real runs. It is kept small here so that # the sampling finishes quickly (and because I know the burn in is short # because of the cheating priors above). -print("Running emcee...") -chain_length = 1000 -sampler.run_mcmc( - start_points, - chain_length, - progress=True, -) +if not rtd: + print("Running emcee...") -# %% -# Merge all the chains together after discarding the first 100 samples as 'burn-in'. -# The burn-in should be decided after looking at the chains in the real world. -samples_emcee = sampler.get_chain(flat=True, discard=100) + chain_length = 1000 + + sampler.run_mcmc( + start_points, + chain_length, + progress=True, + ) # %% -# Plot the MCMC chains to make sure that the burn-in has been removed properly. -# Otherwise, go back and discard more points. -for idx, param_chain in enumerate(samples_emcee.T): - plt.subplot(bt.nparams, 1, idx + 1) - plt.plot(param_chain, label=bt.param_labels[idx]) - plt.legend() -plt.show() +if not rtd: + # Merge all the chains together after discarding the first 100 samples as 'burn-in'. + # The burn-in should be decided after looking at the chains in the real world. + samples_emcee = sampler.get_chain(flat=True, discard=100) + + # %% + # Plot the MCMC chains to make sure that the burn-in has been removed properly. + # Otherwise, go back and discard more points. + for idx, param_chain in enumerate(samples_emcee.T): + plt.subplot(bt.nparams, 1, idx + 1) + plt.plot(param_chain, label=bt.param_labels[idx]) + plt.legend() + plt.show() # %% # Plot the posterior distribution. @@ -123,28 +133,30 @@ # `dlogz` is the target accuracy in the computed Bayesian evidence. # Increasing `npoints` or decreasing `dlogz` gives more accurate results, # but at the cost of time. -print("Running nestle...") -result_nestle_1 = nestle.sample( - bt.lnlikelihood, - bt.prior_transform, - bt.nparams, - method="multi", - npoints=150, - dlogz=0.5, - callback=nestle.print_progress, -) +if not rtd: + print("Running nestle...") + result_nestle_1 = nestle.sample( + bt.lnlikelihood, + bt.prior_transform, + bt.nparams, + method="multi", + npoints=150, + dlogz=0.5, + callback=nestle.print_progress, + ) # %% # Plot the posterior # The nested samples come with weights, which must be taken into account # while plotting. -fig = corner.corner( - result_nestle_1.samples, - weights=result_nestle_1.weights, - labels=bt.param_labels, - range=[0.999] * bt.nparams, -) -plt.show() +if not rtd: + fig = corner.corner( + result_nestle_1.samples, + weights=result_nestle_1.weights, + labels=bt.param_labels, + range=[0.999] * bt.nparams, + ) + plt.show() # %% [markdown] # Let us create a new model with an EFAC applied to all toas (all @@ -175,36 +187,41 @@ print(bt2.likelihood_method) # %% -result_nestle_2 = nestle.sample( - bt2.lnlikelihood, - bt2.prior_transform, - bt2.nparams, - method="multi", - npoints=150, - dlogz=0.5, - callback=nestle.print_progress, -) +if not rtd: + result_nestle_2 = nestle.sample( + bt2.lnlikelihood, + bt2.prior_transform, + bt2.nparams, + method="multi", + npoints=150, + dlogz=0.5, + callback=nestle.print_progress, + ) # %% # Plot the posterior. # The EFAC looks consistent with 1. -fig2 = corner.corner( - result_nestle_2.samples, - weights=result_nestle_2.weights, - labels=bt2.param_labels, - range=[0.999] * bt2.nparams, -) -plt.show() +if not rtd: + fig2 = corner.corner( + result_nestle_2.samples, + weights=result_nestle_2.weights, + labels=bt2.param_labels, + range=[0.999] * bt2.nparams, + ) + plt.show() # %% [markdown] # Now let us look at the evidences and compute the Bayes factor. # %% -print(f"Evidence without EFAC : {result_nestle_1.logz} +/- {result_nestle_1.logzerr}") -print(f"Evidence with EFAC : {result_nestle_2.logz} +/- {result_nestle_2.logzerr}") +if not rtd: + print( + f"Evidence without EFAC : {result_nestle_1.logz} +/- {result_nestle_1.logzerr}" + ) + print(f"Evidence with EFAC : {result_nestle_2.logz} +/- {result_nestle_2.logzerr}") -bf = np.exp(result_nestle_1.logz - result_nestle_2.logz) -print(f"Bayes factor : {bf} (in favor of no EFAC)") + bf = np.exp(result_nestle_1.logz - result_nestle_2.logz) + print(f"Bayes factor : {bf} (in favor of no EFAC)") # %% [markdown] # The Bayes factor tells us that the EFAC is unnecessary for this dataset. diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index 85f949eea..f576b9546 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -72,28 +72,36 @@ ) # %% -print("Running emcee...") -chain_length = 1000 -sampler.run_mcmc( - start_points, - chain_length, - progress=True, -) +# ** IMPORTANT!!! ** +# This is used to exclude the following time-consuming steps from the readthedocs build. +# Set this to False while actually using this example. +rtd = True # %% -samples_emcee = sampler.get_chain(flat=True, discard=100) +if not rtd: + print("Running emcee...") + chain_length = 1000 + sampler.run_mcmc( + start_points, + chain_length, + progress=True, + ) + + samples_emcee = sampler.get_chain(flat=True, discard=100) # %% # Plot the chains to make sure they have converged and the burn-in has been removed properly. -for idx, param_chain in enumerate(samples_emcee.T): - plt.subplot(bt.nparams, 1, idx + 1) - plt.plot(param_chain) - plt.ylabel(bt.param_labels[idx]) - plt.autoscale() -plt.show() +if not rtd: + for idx, param_chain in enumerate(samples_emcee.T): + plt.subplot(bt.nparams, 1, idx + 1) + plt.plot(param_chain) + plt.ylabel(bt.param_labels[idx]) + plt.autoscale() + plt.show() # %% -fig = corner.corner( - samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params -) -plt.show() +if not rtd: + fig = corner.corner( + samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params + ) + plt.show() From 474876c0cce37243adf55602c69597b9f568d1a4 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 21 Sep 2023 11:30:10 -0500 Subject: [PATCH 27/30] -- --- docs/examples/bayesian-example-NGC6440E.py | 15 +++++++++++++++ docs/examples/bayesian-wideband-example.py | 18 ++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py index d371c7316..f6d171608 100644 --- a/docs/examples/bayesian-example-NGC6440E.py +++ b/docs/examples/bayesian-example-NGC6440E.py @@ -1,3 +1,18 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.7 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + # %% [markdown] # # PINT Bayesian Interface Examples diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index f576b9546..785b28103 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -1,3 +1,21 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.4 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% [markdown] +# # PINT Bayesian Interface Example (Wideband) + # %% import corner import emcee From 928a296a16d18f0f3843240f546e45bb2fbf6ff7 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 21 Sep 2023 11:36:09 -0500 Subject: [PATCH 28/30] -- --- docs/examples/bayesian-example-NGC6440E.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py index f6d171608..0c45189ae 100644 --- a/docs/examples/bayesian-example-NGC6440E.py +++ b/docs/examples/bayesian-example-NGC6440E.py @@ -128,8 +128,9 @@ # %% # Plot the posterior distribution. -fig = corner.corner(samples_emcee, labels=bt.param_labels) -plt.show() +if not rtd: + fig = corner.corner(samples_emcee, labels=bt.param_labels) + plt.show() # %% [markdown] # ## Nested sampling with nestle From 8195855d0aadb8c7e08f0b3e3f637e4266999a82 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 21 Sep 2023 11:43:48 -0500 Subject: [PATCH 29/30] fix nb --- docs/examples/bayesian-example-NGC6440E.py | 1 + docs/examples/bayesian-wideband-example.py | 15 --------------- 2 files changed, 1 insertion(+), 15 deletions(-) diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py index 0c45189ae..7b1952736 100644 --- a/docs/examples/bayesian-example-NGC6440E.py +++ b/docs/examples/bayesian-example-NGC6440E.py @@ -118,6 +118,7 @@ samples_emcee = sampler.get_chain(flat=True, discard=100) # %% +if not rtd: # Plot the MCMC chains to make sure that the burn-in has been removed properly. # Otherwise, go back and discard more points. for idx, param_chain in enumerate(samples_emcee.T): diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py index 785b28103..6ab04bb32 100644 --- a/docs/examples/bayesian-wideband-example.py +++ b/docs/examples/bayesian-wideband-example.py @@ -1,18 +1,3 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.14.4 -# kernelspec: -# display_name: Python 3 (ipykernel) -# language: python -# name: python3 -# --- - # %% [markdown] # # PINT Bayesian Interface Example (Wideband) From e43f10fe8292f04476f4fddcdadfec61a4c1e769 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 21 Sep 2023 12:33:23 -0500 Subject: [PATCH 30/30] tutorials.rst --- docs/tutorials.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/tutorials.rst b/docs/tutorials.rst index fea519176..06f7237f5 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -67,6 +67,8 @@ are not included in the default build because they take too long, but you can do examples/check_phase_connection.ipynb examples/PINT_observatories.ipynb examples/solar_wind.ipynb + examples/bayesian-example-NGC6440E.py + examples/bayesian-wideband-example.py examples-rendered/paper_validation_example.ipynb .. _`Time a Pulsar`: examples/time_a_pulsar.html