Skip to content

Commit

Permalink
Merge pull request #1692 from dlakaplan/derivedpars
Browse files Browse the repository at this point in the history 
Moved `get_derived_params` to `timing_model`
  • Loading branch information
@abhisrkckl
abhisrkckl authored Dec 18, 2023
2 parents 47026ba + d6634b7 commit 2d8e13d
Show file tree
Hide file tree
Showing 4 changed files with 248 additions and 204 deletions.
1 change: 1 addition & 0 deletions CHANGELOG-unreleased.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
225 changes: 22 additions & 203 deletions src/pint/fitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,215 +469,34 @@ 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"""
def get_derived_params(self, returndict=False):
"""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"
Parameters
----------
returndict : bool, optional
Whether to only return the string of results or also a dictionary
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"
Returns
-------
results : str
parameters : dict, optional
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
if self.model.OMDOT.uncertainty is not None
else 0 * self.model.OMDOT.quantity.unit
)
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"
See Also
--------
:func:`pint.models.timing_model.TimingModel.get_derived_params`
"""

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
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."""
Expand Down
Loading

0 comments on commit 2d8e13d

Please sign in to comment.