From 82864ef2875ca6fcf53fc680b8578d72c15f5f0c Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Sun, 10 Dec 2023 16:03:28 -0600 Subject: [PATCH 1/6] moved get_derived_params --- CHANGELOG-unreleased.md | 1 + src/pint/fitter.py | 212 ++------------------------------ src/pint/models/timing_model.py | 194 ++++++++++++++++++++++++++++- 3 files changed, 203 insertions(+), 204 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 4a43942b2..e89255740 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,6 +9,7 @@ the released changes. ## Unreleased ### Changed +- Moved `get_derived_params` to `timing_model` ### Added - Added numdifftools to setup.cfg to match requirements.txt - Documentation: Added `convert_parfile` to list of command-line tools in RTD diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 287f7f467..6dae2191c 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -464,211 +464,17 @@ def get_summary(self, nodmx=False): ufloat(par.value, par.uncertainty.value), par.units, ) - s += "\n" + self.get_derived_params() + s += "\n" + self.model.get_derived_params() return s - def get_derived_params(self): - """Return a string with various derived parameters from the fitted model""" - - import uncertainties.umath as um - from uncertainties import ufloat - - # Now print some useful derived parameters - s = "Derived Parameters:\n" - if hasattr(self.model, "F0"): - F0 = self.model.F0.quantity - if not self.model.F0.frozen: - p, perr = pint.derived_quantities.pferrs(F0, self.model.F0.uncertainty) - s += f"Period = {p.to(u.s)} +/- {perr.to(u.s)}\n" - else: - s += f"Period = {(1.0 / F0).to(u.s)}\n" - if hasattr(self.model, "F1"): - F1 = self.model.F1.quantity - if not any([self.model.F1.frozen, self.model.F0.frozen]): - p, perr, pd, pderr = pint.derived_quantities.pferrs( - F0, self.model.F0.uncertainty, F1, self.model.F1.uncertainty - ) - s += f"Pdot = {pd.to(u.dimensionless_unscaled)} +/- {pderr.to(u.dimensionless_unscaled)}\n" - if F1.value < 0.0: # spinning-down - brakingindex = 3 - s += f"Characteristic age = {pint.derived_quantities.pulsar_age(F0, F1, n=brakingindex):.4g} (braking index = {brakingindex})\n" - s += f"Surface magnetic field = {pint.derived_quantities.pulsar_B(F0, F1):.3g}\n" - s += f"Magnetic field at light cylinder = {pint.derived_quantities.pulsar_B_lightcyl(F0, F1):.4g}\n" - I_NS = I = 1.0e45 * u.g * u.cm**2 - s += f"Spindown Edot = {pint.derived_quantities.pulsar_edot(F0, F1, I=I_NS):.4g} (I={I_NS})\n" - else: - s += "Not computing Age, B, or Edot since F1 > 0.0\n" - - if hasattr(self.model, "PX") and not self.model.PX.frozen: - s += "\n" - px = self.model.PX.as_ufloat(u.arcsec) - s += f"Parallax distance = {1.0/px:.3uP} pc\n" - - # Now binary system derived parameters - if self.model.is_binary: - for x in self.model.components: - if x.startswith("Binary"): - binary = x - - s += f"\nBinary model {binary}\n" - - btx = False - if ( - hasattr(self.model, "FB0") - and self.model.FB0.quantity is not None - and self.model.FB0.value != 0.0 - ): - btx = True - FB0 = self.model.FB0.quantity - if not self.model.FB0.frozen: - pb, pberr = pint.derived_quantities.pferrs( - FB0, self.model.FB0.uncertainty - ) - s += f"Orbital Period (PB) = {pb.to(u.d)} +/- {pberr.to(u.d)}\n" - else: - s += f"Orbital Period (PB) = {(1.0 / FB0).to(u.d)}\n" - - if ( - hasattr(self.model, "FB1") - and self.model.FB1.quantity is not None - and self.model.FB1.value != 0.0 - ): - FB1 = self.model.FB1.quantity - if not any([self.model.FB1.frozen, self.model.FB0.frozen]): - pb, pberr, pbd, pbderr = pint.derived_quantities.pferrs( - FB0, self.model.FB0.uncertainty, FB1, self.model.FB1.uncertainty - ) - s += f"Orbital Pdot (PBDOT) = {pbd.to(u.dimensionless_unscaled)} +/- {pbderr.to(u.dimensionless_unscaled)}\n" - - ell1 = False - if binary.startswith("BinaryELL1"): - ell1 = True - eps1 = self.model.EPS1.as_ufloat() - eps2 = self.model.EPS2.as_ufloat() - tasc = ufloat( - # This is a time in MJD - self.model.TASC.quantity.mjd, - self.model.TASC.uncertainty.to(u.d).value, - ) - if hasattr(self.model, "PB") and self.model.PB.value is not None: - pb = self.model.PB.as_ufloat(u.d) - elif hasattr(self.model, "FB0") and self.model.FB0.value is not None: - pb = 1 / self.model.FB0.as_ufloat(1 / u.d) - s += "Conversion from ELL1 parameters:\n" - ecc = um.sqrt(eps1**2 + eps2**2) - s += "ECC = {:P}\n".format(ecc) - om = um.atan2(eps1, eps2) * 180.0 / np.pi - if om < 0.0: - om += 360.0 - s += f"OM = {om:P} deg\n" - t0 = tasc + pb * om / 360.0 - s += f"T0 = {t0:SP}\n" - - a1 = ( - self.model.A1.quantity - if self.model.A1.quantity is not None - else 0 * pint.ls - ) - if self.is_wideband: - s += pint.utils.ELL1_check( - a1, - ecc.nominal_value * u.s / u.s, - self.resids.toa.rms_weighted(), - self.toas.ntoas, - outstring=True, - ) - else: - s += pint.utils.ELL1_check( - a1, - ecc.nominal_value * u.s / u.s, - self.resids.rms_weighted(), - self.toas.ntoas, - outstring=True, - ) - s += "\n" - if hasattr(self.model, "FB0") and self.model.FB0.value is not None: - pb, pberr = pint.derived_quantities.pferrs( - self.model.FB0.quantity, self.model.FB0.uncertainty - ) - # Masses and inclination - pb = pb.to(u.d) if btx else self.model.PB.quantity - pberr = pberr.to(u.d) if btx else self.model.PB.uncertainty - if not self.model.A1.frozen: - pbs = ufloat( - pb.to(u.s).value, - pberr.to(u.s).value, - ) - a1 = self.model.A1.as_ufloat(pint.ls) - # This is the mass function, done explicitly so that we get - # uncertainty propagation automatically. - # TODO: derived quantities funcs should take uncertainties - fm = 4.0 * np.pi**2 * a1**3 / (4.925490947e-6 * pbs**2) - s += f"Mass function = {fm:SP} Msun\n" - mcmed = pint.derived_quantities.companion_mass( - pb, - self.model.A1.quantity, - i=60.0 * u.deg, - mp=1.4 * u.solMass, - ) - mcmin = pint.derived_quantities.companion_mass( - pb, - self.model.A1.quantity, - i=90.0 * u.deg, - mp=1.4 * u.solMass, - ) - s += f"Min / Median Companion mass (assuming Mpsr = 1.4 Msun) = {mcmin.value:.4f} / {mcmed.value:.4f} Msun\n" - - if ( - hasattr(self.model, "OMDOT") - and self.model.OMDOT.quantity is not None - and self.model.OMDOT.value != 0.0 - ): - omdot = self.model.OMDOT.quantity - omdot_err = self.model.OMDOT.uncertainty - ecc = ( - ecc.n * u.dimensionless_unscaled - if ell1 - else self.model.ECC.quantity - ) - Mtot = pint.derived_quantities.omdot_to_mtot(omdot, pb, ecc) - # Assume that the uncertainty on OMDOT dominates the Mtot uncertainty - # This is probably a good assumption until we can get the uncertainties module - # to work with quantities. - Mtot_hi = pint.derived_quantities.omdot_to_mtot( - omdot + omdot_err, - pb, - ecc, - ) - Mtot_lo = pint.derived_quantities.omdot_to_mtot( - omdot - omdot_err, - pb, - ecc, - ) - Mtot_err = max(abs(Mtot_hi - Mtot), abs(Mtot - Mtot_lo)) - mt = ufloat(Mtot.value, Mtot_err.value) - s += f"Total mass, assuming GR, from OMDOT is {mt:SP} Msun\n" - - if ( - hasattr(self.model, "SINI") - and self.model.SINI.quantity is not None - and (self.model.SINI.value >= 0.0 and self.model.SINI.value < 1.0) - ): - with contextlib.suppress(TypeError, ValueError): - # Put this in a try in case SINI is UNSET or an illegal value - if not self.model.SINI.frozen: - si = self.model.SINI.as_ufloat() - s += f"From SINI in model:\n" - s += f" cos(i) = {um.sqrt(1 - si**2):SP}\n" - s += f" i = {um.asin(si) * 180.0 / np.pi:SP} deg\n" - - psrmass = pint.derived_quantities.pulsar_mass( - pb, - self.model.A1.quantity, - self.model.M2.quantity, - np.arcsin(self.model.SINI.quantity), - ) - s += f"Pulsar mass (Shapiro Delay) = {psrmass}" - return s + def get_derived_params(self, returndict=False): + return self.model.get_derived_params( + rms=self.resids.toa.rms_weighted() + if self.is_wideband + else self.resids.rms_weighted(), + ntoas=self.toas.ntoas, + returndict=returndict, + ) def print_summary(self): """Write a summary of the TOAs to stdout.""" diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 6dd7506bc..7a2764a3e 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -31,13 +31,14 @@ import abc import copy import inspect +import contextlib from collections import OrderedDict, defaultdict from functools import wraps from warnings import warn from uncertainties import ufloat import astropy.time as time -import astropy.units as u +from astropy import units as u, constants as c import numpy as np from astropy.utils.decorators import lazyproperty import astropy.coordinates as coords @@ -2947,6 +2948,197 @@ def as_ICRS(self, epoch=None): return new_model + def get_derived_params(self, rms=None, ntoas=None, returndict=False): + """Return a string with various derived parameters from the fitted model""" + + import uncertainties.umath as um + from uncertainties import ufloat + + outdict = {} + + # Now print some useful derived parameters + s = "Derived Parameters:\n" + if hasattr(self, "F0"): + F0 = self.F0.as_ufloat() + p = 1 / F0 + s += f"Period = {p:P} s\n" + outdict["P (s)"] = p + if hasattr(self, "F1"): + F1 = self.F1.as_ufloat() + pdot = -F1 / F0**2 + outdict["Pdot (s/s)"] = pdot + s += f"Pdot = {pdot:P}\n" + if self.F1.value < 0.0: # spinning-down + brakingindex = 3 + s += f"Characteristic age = {pint.derived_quantities.pulsar_age(self.F0.quantity, self.F1.quantity, n=brakingindex):.4g} (braking index = {brakingindex})\n" + s += f"Surface magnetic field = {pint.derived_quantities.pulsar_B(self.F0.quantity, self.F1.quantity):.3g}\n" + s += f"Magnetic field at light cylinder = {pint.derived_quantities.pulsar_B_lightcyl(self.F0.quantity, self.F1.quantity):.4g}\n" + I_NS = I = 1.0e45 * u.g * u.cm**2 + s += f"Spindown Edot = {pint.derived_quantities.pulsar_edot(self.F0.quantity, self.F1.quantity, I=I_NS):.4g} (I={I_NS})\n" + outdict["age"] = pint.derived_quantities.pulsar_age( + self.F0.quantity, self.F1.quantity, n=brakingindex + ) + outdict["B"] = pint.derived_quantities.pulsar_B( + self.F0.quantity, self.F1.quantity + ) + outdict["Blc"] = pint.derived_quantities.pulsar_B_lightcyl( + self.F0.quantity, self.F1.quantity + ) + outdict["Edot"] = pint.derived_quantities.pulsar_B_lightcyl( + self.F0.quantity, self.F1.quantity + ) + else: + s += "Not computing Age, B, or Edot since F1 > 0.0\n" + + if hasattr(self, "PX") and not self.PX.frozen: + s += "\n" + px = self.PX.as_ufloat(u.arcsec) + s += f"Parallax distance = {1.0/px:.3uP} pc\n" + outdict["Dist (pc):P"] = 1.0 / px + # Now binary system derived parameters + if self.is_binary: + for x in self.components: + if x.startswith("Binary"): + binary = x + + s += f"\nBinary model {binary}\n" + outdict["Binary"] = binary + + btx = False + if ( + hasattr(self, "FB0") + and self.FB0.quantity is not None + and self.FB0.value != 0.0 + ): + btx = True + pb = 1 / self.FB0.as_ufloat(1 / u.d) + s += f"Orbital Period (PB) = {pb:P} (d)\n" + else: + pb = self.PB.as_ufloat(u.d) + outdict["PB (d)"] = pb + + pbdot = None + if ( + hasattr(self, "FB1") + and self.FB1.quantity is not None + and self.FB1.value != 0.0 + ): + pbdot = -self.FB1.as_ufloat(u.Hz / u.s) / self.FB0.as_ufloat(u.Hz) ** 2 + elif ( + hasattr(self, "PBDOT") + and self.PBDOT.quantity is not None + and self.PBDOT.value != 0 + ): + pbdot = self.PBDOT.as_ufloat(u.s / u.s) + + if pbdot is not None: + s += f"Orbital Pdot (PBDOT) = {pbdot:P} (s/s)\n" + outdict["PBDOT (s/s)"] = pbdot + + ell1 = False + if binary.startswith("BinaryELL1"): + ell1 = True + eps1 = self.EPS1.as_ufloat() + eps2 = self.EPS2.as_ufloat() + tasc = ufloat( + # This is a time in MJD + self.TASC.quantity.mjd, + self.TASC.uncertainty.to(u.d).value + if self.TASC.uncertainty is not None + else 0, + ) + s += "Conversion from ELL1 parameters:\n" + ecc = um.sqrt(eps1**2 + eps2**2) + s += "ECC = {:P}\n".format(ecc) + outdict["ECC"] = ecc + om = um.atan2(eps1, eps2) * 180.0 / np.pi + if om < 0.0: + om += 360.0 + s += f"OM = {om:P} deg\n" + outdict["OM (deg)"] = om + t0 = tasc + pb * om / 360.0 + s += f"T0 = {t0:SP}\n" + outdict["T0"] = t0 + + a1 = self.A1.quantity if self.A1.quantity is not None else 0 * pint.ls + if rms is not None and ntoas is not None: + s += pint.utils.ELL1_check( + a1, + ecc.nominal_value * u.s / u.s, + rms, + ntoas, + outstring=True, + ) + s += "\n" + # Masses and inclination + if not self.A1.frozen: + a1 = self.A1.as_ufloat(pint.ls) + # This is the mass function, done explicitly so that we get + # uncertainty propagation automatically. + # TODO: derived quantities funcs should take uncertainties + fm = 4.0 * np.pi**2 * a1**3 / (4.925490947e-6 * pb**2) + s += f"Mass function = {fm:SP} Msun\n" + outdict["Mass Function (Msun)"] = fm + mcmed = pint.derived_quantities.companion_mass( + pb.n * u.d, + self.A1.quantity, + i=60.0 * u.deg, + mp=1.4 * u.solMass, + ) + mcmin = pint.derived_quantities.companion_mass( + pb.n * u.d, + self.A1.quantity, + i=90.0 * u.deg, + mp=1.4 * u.solMass, + ) + s += f"Min / Median Companion mass (assuming Mpsr = 1.4 Msun) = {mcmin.value:.4f} / {mcmed.value:.4f} Msun\n" + outdict["Mc,med (Msun)"] = mcmed.value + outdict["Mc,min (Msun)"] = mcmin.value + + if ( + hasattr(self, "OMDOT") + and self.OMDOT.quantity is not None + and self.OMDOT.value != 0.0 + ): + omdot = self.OMDOT.as_ufloat(u.rad / u.s) + e = ecc if ell1 else self.ECC.as_ufloat() + mt = ( + ( + omdot + / ( + 3 + * (c.G * u.Msun / c.c**3).to_value(u.s) ** (2.0 / 3) + * ((pb * 86400 / 2 / np.pi)) ** (-5.0 / 3) + * (1 - e**2) ** -1 + ) + ) + ) ** (3.0 / 2) + s += f"Total mass, assuming GR, from OMDOT is {mt:SP} Msun\n" + outdict["Mtot (Msun)"] = mt + + if ( + hasattr(self, "SINI") + and self.SINI.quantity is not None + and (self.SINI.value >= 0.0 and self.SINI.value < 1.0) + ): + with contextlib.suppress(TypeError, ValueError): + # Put this in a try in case SINI is UNSET or an illegal value + if not self.SINI.frozen: + si = self.SINI.as_ufloat() + s += f"From SINI in model:\n" + s += f" cos(i) = {um.sqrt(1 - si**2):SP}\n" + s += f" i = {um.asin(si) * 180.0 / np.pi:SP} deg\n" + + psrmass = pint.derived_quantities.pulsar_mass( + pb.n * u.d, + self.A1.quantity, + self.M2.quantity, + np.arcsin(self.SINI.quantity), + ) + s += f"Pulsar mass (Shapiro Delay) = {psrmass}" + outdict["Mp (Msun)"] = psrmass + return s if not returndict else s, outdict + class ModelMeta(abc.ABCMeta): """Ensure timing model registration. From c2feff440ac2fc9a5a8ce6eb2bcbc86a7f26425a Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Sun, 10 Dec 2023 16:23:14 -0600 Subject: [PATCH 2/6] added docstrings --- src/pint/fitter.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 6dae2191c..b550225e9 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -468,6 +468,23 @@ def get_summary(self, nodmx=False): return s def get_derived_params(self, returndict=False): + """Return a string with various derived parameters from the fitted model + + Parameters + ---------- + returndict : bool, optional + Whether to only return the string of results or also a dictionary + + Returns + ------- + results : str + parameters : dict, optional + + See Also + -------- + :func:`pint.models.timing_model.TimingModel.get_derived_params` + """ + return self.model.get_derived_params( rms=self.resids.toa.rms_weighted() if self.is_wideband From 58588ab9f945760930c5518d56bc801ec899f2a1 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Sun, 10 Dec 2023 17:35:48 -0600 Subject: [PATCH 3/6] fix bug --- src/pint/models/timing_model.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 7a2764a3e..79ab9b813 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -2949,7 +2949,22 @@ def as_ICRS(self, epoch=None): return new_model def get_derived_params(self, rms=None, ntoas=None, returndict=False): - """Return a string with various derived parameters from the fitted model""" + """Return a string with various derived parameters from the fitted model + + Parameters + ---------- + rms : astropy.units.Quantity, optional + RMS of fit for checking ELL1 validity + ntoas : int, optional + Number of TOAs for checking ELL1 validity + returndict : bool, optional + Whether to only return the string of results or also a dictionary + + Returns + ------- + results : str + parameters : dict, optional + """ import uncertainties.umath as um from uncertainties import ufloat @@ -3137,7 +3152,9 @@ def get_derived_params(self, rms=None, ntoas=None, returndict=False): ) s += f"Pulsar mass (Shapiro Delay) = {psrmass}" outdict["Mp (Msun)"] = psrmass - return s if not returndict else s, outdict + if not returndict: + return s + return s, outdict class ModelMeta(abc.ABCMeta): From b1a8d81e125df9aa22f17098dca35ade68ccc9bb Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Sun, 10 Dec 2023 19:07:49 -0600 Subject: [PATCH 4/6] fixed units --- src/pint/models/timing_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 79ab9b813..870a46652 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -3091,7 +3091,7 @@ def get_derived_params(self, rms=None, ntoas=None, returndict=False): # This is the mass function, done explicitly so that we get # uncertainty propagation automatically. # TODO: derived quantities funcs should take uncertainties - fm = 4.0 * np.pi**2 * a1**3 / (4.925490947e-6 * pb**2) + fm = 4.0 * np.pi**2 * a1**3 / (4.925490947e-6 * (pb * 86400) ** 2) s += f"Mass function = {fm:SP} Msun\n" outdict["Mass Function (Msun)"] = fm mcmed = pint.derived_quantities.companion_mass( From 67e1d21e4ddfb259deb5005c8cb5188b06418d5f Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 12 Dec 2023 18:35:43 -0600 Subject: [PATCH 5/6] tests --- tests/test_derivedparams.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 tests/test_derivedparams.py diff --git a/tests/test_derivedparams.py b/tests/test_derivedparams.py new file mode 100644 index 000000000..3e7ae8680 --- /dev/null +++ b/tests/test_derivedparams.py @@ -0,0 +1,15 @@ +import glob +import pytest +from pint.models import get_model +from pinttestdata import datadir + + +@pytest.mark.parametrize("parname", list(glob.glob(str(datadir) + "/*.par"))) +def test_derived_params(parname): + try: + m = get_model(parname) + readin = True + except: + readin = False + if readin: + out = m.get_derived_params() From 21d1b7ad1112ad07413a5a0327fc965f7e3de014 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 13 Dec 2023 09:15:04 -0600 Subject: [PATCH 6/6] fixed format of output --- src/pint/models/timing_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 870a46652..79876de27 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -3009,7 +3009,7 @@ def get_derived_params(self, rms=None, ntoas=None, returndict=False): s += "\n" px = self.PX.as_ufloat(u.arcsec) s += f"Parallax distance = {1.0/px:.3uP} pc\n" - outdict["Dist (pc):P"] = 1.0 / px + outdict["Dist (pc)"] = 1.0 / px # Now binary system derived parameters if self.is_binary: for x in self.components: