From 86c81408a803c29b7a03d6d2960b65f66bd8ea99 Mon Sep 17 00:00:00 2001 From: Deven Bhakta Date: Thu, 17 Oct 2024 13:17:06 -0400 Subject: [PATCH 1/4] Phases calculation using the Designmatrix feature --- src/pint/plot_utils.py | 6 +++-- src/pint/scripts/event_optimize.py | 35 ++++++++++++++++++++++++++++-- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/src/pint/plot_utils.py b/src/pint/plot_utils.py index 45eba3977..917412b71 100644 --- a/src/pint/plot_utils.py +++ b/src/pint/plot_utils.py @@ -279,11 +279,13 @@ def plot_priors( for i in range(len(keys[:-1])): values[i] = values[i][burnin:].flatten() x_range.append(np.linspace(values[i].min(), values[i].max(), num=bins)) - priors.append(getattr(model, keys[i]).prior.pdf(x_range[i])) + priors.append(getattr(model, keys[i]).prior.logpdf(x_range[i])) a, x = np.histogram(values[i], bins=bins, density=True) counts.append(a) - fig, axs = plt.subplots(len(keys), figsize=(8, 11), constrained_layout=True) + fig, axs = plt.subplots( + len(keys), figsize=(8, len(keys) * 1.5), constrained_layout=True + ) for i, p in enumerate(keys): if i != len(keys[:-1]): diff --git a/src/pint/scripts/event_optimize.py b/src/pint/scripts/event_optimize.py index e4c9dd4ff..1b66b2229 100755 --- a/src/pint/scripts/event_optimize.py +++ b/src/pint/scripts/event_optimize.py @@ -414,6 +414,16 @@ def __init__( self.model, phs, phserr ) self.n_fit_params = len(self.fitvals) + self.M, _, _ = self.model.designmatrix(self.toas) + self.M = self.M.transpose() * -self.model.F0.value + self.phases = self.get_event_phases() + self.calc_phase = False + + def calc_phase_matrix(self, theta): + d_phs = np.zeros(len(self.toas)) + for i in range(len(theta) - 1): + d_phs += self.M[i + 1] * (self.fitvals[i] - theta[i]) + return (self.phases - d_phs) % 1 def get_event_phases(self): """ @@ -446,7 +456,10 @@ def lnposterior(self, theta): return -np.inf, -np.inf, -np.inf # Call PINT to compute the phases - phases = self.get_event_phases() + if self.calc_phase: + phases = self.calc_phase_matrix(theta) + else: + phases = self.get_event_phases() lnlikelihood = profile_likelihood( theta[-1], self.xtemp, phases, self.template, self.weights ) @@ -686,6 +699,13 @@ def main(argv=None): action="store_true", dest="noautocorr", ) + parser.add_argument( + "--calc_phase", + help="Calculates the phase at each MCMC step using the designmatrix", + default=False, + action="store_true", + dest="calc_phase", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -862,6 +882,9 @@ def main(argv=None): # This way, one walker should always be in a good position pos[0] = ftr.fitvals + # How phase will be calculated at each step (either with the designmatrix or ) + ftr.calc_phase = True if args.calc_phase else False + import emcee # Setting up a backend to save the chains into an h5 file @@ -925,7 +948,7 @@ def chains_to_dict(names, sampler): def plot_chains(chain_dict, file=False): npts = len(chain_dict) - fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, 9)) + fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, npts * 1.5)) for ii, name in enumerate(chain_dict.keys()): axes[ii].plot(chain_dict[name], color="k", alpha=0.3) axes[ii].set_ylabel(name) @@ -950,6 +973,7 @@ def plot_chains(chain_dict, file=False): lnprior_samps = blobs["lnprior"] lnlikelihood_samps = blobs["lnlikelihood"] lnpost_samps = lnprior_samps + lnlikelihood_samps + maxpost = lnpost_samps[:][burnin:].max() ind = np.unravel_index( np.argmax(lnpost_samps[:][burnin:]), lnpost_samps[:][burnin:].shape ) @@ -1000,8 +1024,15 @@ def plot_chains(chain_dict, file=False): ] ftr.set_param_uncertainties(dict(zip(ftr.fitkeys[:-1], errors[:-1]))) + # Calculating the AIC and BIC + n_params = len(ftr.model.free_params) + AIC = 2 * (n_params - maxpost) + BIC = n_params * np.log(len(ts)) - 2 * maxpost + ftr.model.NTOA.value = ts.ntoas f = open(filename + "_post.par", "w") f.write(ftr.model.as_parfile()) + f.write(f"\n#The AIC is {AIC}") + f.write(f"\n#The BIC is {BIC}") f.close() # Print the best MCMC values and ranges From eb43e6c0e5a8a232060bc706c39b01b2f297142b Mon Sep 17 00:00:00 2001 From: Deven Bhakta Date: Thu, 24 Oct 2024 11:24:42 -0400 Subject: [PATCH 2/4] Renamed calc_phase option to linearize_model option within event_optimze --- CHANGELOG-unreleased.md | 2 ++ src/pint/scripts/event_optimize.py | 12 ++++++------ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 2899b99f8..991f5caad 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,6 +9,8 @@ the released changes. ## Unreleased ### Changed +- Updated the `plot_priors` function in `plot_utils.py` and `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters. ### Added +- Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix. ### Fixed ### Removed diff --git a/src/pint/scripts/event_optimize.py b/src/pint/scripts/event_optimize.py index 1b66b2229..071860a8d 100755 --- a/src/pint/scripts/event_optimize.py +++ b/src/pint/scripts/event_optimize.py @@ -417,7 +417,7 @@ def __init__( self.M, _, _ = self.model.designmatrix(self.toas) self.M = self.M.transpose() * -self.model.F0.value self.phases = self.get_event_phases() - self.calc_phase = False + self.linearize_model = False def calc_phase_matrix(self, theta): d_phs = np.zeros(len(self.toas)) @@ -456,7 +456,7 @@ def lnposterior(self, theta): return -np.inf, -np.inf, -np.inf # Call PINT to compute the phases - if self.calc_phase: + if self.linearize_model: phases = self.calc_phase_matrix(theta) else: phases = self.get_event_phases() @@ -700,11 +700,11 @@ def main(argv=None): dest="noautocorr", ) parser.add_argument( - "--calc_phase", + "--linearize_model", help="Calculates the phase at each MCMC step using the designmatrix", default=False, action="store_true", - dest="calc_phase", + dest="linearize_model", ) args = parser.parse_args(argv) @@ -882,8 +882,8 @@ def main(argv=None): # This way, one walker should always be in a good position pos[0] = ftr.fitvals - # How phase will be calculated at each step (either with the designmatrix or ) - ftr.calc_phase = True if args.calc_phase else False + # How phase will be calculated at each step (either with the designmatrix orthe exact phase calculation) + ftr.linearize_model = args.linearize_model import emcee From 827bcbbabe0b4977feecf8b1d14646f746ac2762 Mon Sep 17 00:00:00 2001 From: Deven Bhakta Date: Thu, 24 Oct 2024 11:26:55 -0400 Subject: [PATCH 3/4] Changelog entry update --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 991f5caad..366c272e0 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -12,5 +12,6 @@ the released changes. - Updated the `plot_priors` function in `plot_utils.py` and `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters. ### Added - Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix. +- Added AIC and BIC calculation to be written in the post fit parfile from `event_optimize` ### Fixed ### Removed From 85b27a22fab1584099361b74f57828058fa44be8 Mon Sep 17 00:00:00 2001 From: Deven Bhakta Date: Thu, 14 Nov 2024 10:39:34 -0500 Subject: [PATCH 4/4] Revert plot_priors change --- CHANGELOG-unreleased.md | 2 +- src/pint/plot_utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index f5892a3f9..09561bfca 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,7 +9,7 @@ the released changes. ## Unreleased ### Changed -- Updated the `plot_priors` function in `plot_utils.py` and `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters. +- Updated the `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters. ### Added - Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix. - Added AIC and BIC calculation to be written in the post fit parfile from `event_optimize` diff --git a/src/pint/plot_utils.py b/src/pint/plot_utils.py index 917412b71..289655f6d 100644 --- a/src/pint/plot_utils.py +++ b/src/pint/plot_utils.py @@ -279,7 +279,7 @@ def plot_priors( for i in range(len(keys[:-1])): values[i] = values[i][burnin:].flatten() x_range.append(np.linspace(values[i].min(), values[i].max(), num=bins)) - priors.append(getattr(model, keys[i]).prior.logpdf(x_range[i])) + priors.append(getattr(model, keys[i]).prior.pdf(x_range[i])) a, x = np.histogram(values[i], bins=bins, density=True) counts.append(a)