From 0bbdaf3548446ecb82893b9ec9a08a05aa746e77 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 10:47:03 -0500 Subject: [PATCH 01/19] to_value --- src/pint/residuals.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/pint/residuals.py b/src/pint/residuals.py index 736c5798b..a7876a4e9 100644 --- a/src/pint/residuals.py +++ b/src/pint/residuals.py @@ -577,15 +577,11 @@ def calc_chi2(self): # This is more correct way, but it is the slowest. # return (((self.time_resids / self.toas.get_errors()).decompose()**2.0).sum()).value - # This method is faster then the method above but not the most correct way - # return ((self.time_resids.to(u.s) / self.toas.get_errors().to(u.s)).value**2.0).sum() - # This the fastest way, but highly depend on the assumption of time_resids and # error units. Ensure only a pure number is returned. - try: - return ((self.time_resids / toa_errors.to(u.s)) ** 2.0).sum().value - except ValueError: - return ((self.time_resids / toa_errors.to(u.s)) ** 2.0).sum() + return ( + (self.time_resids.to_value(u.s) / toa_errors.to_value(u.s)) ** 2.0 + ).sum() def ecorr_average(self, use_noise_model=True): """Uses the ECORR noise model time-binning to compute "epoch-averaged" residuals. @@ -763,10 +759,12 @@ def calc_chi2(self): data_errors = self.get_data_error() if (data_errors == 0.0).any(): return np.inf - try: - return ((self.resids / data_errors) ** 2.0).sum().decompose().value - except ValueError: - return ((self.resids / data_errors) ** 2.0).sum().decompose() + + return ( + ((self.resids / data_errors) ** 2.0) + .sum() + .to_value(u.dimensionless_unscaled) + ) def rms_weighted(self): """Compute weighted RMS of the residuals in time.""" From da8df8dfea6d56f056f6320df82ae0c04c71e5ca Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 13:03:04 -0500 Subject: [PATCH 02/19] optimize astrometry (icrs) derivatives --- src/pint/models/astrometry.py | 99 ++++++++++++++++++++++++----------- 1 file changed, 69 insertions(+), 30 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 1d737849c..2eeb6e9ef 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -145,30 +145,68 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): ) return delay * u.second - def get_d_delay_quantities(self, toas): + def get_d_delay_quantities(self, toas, include="all"): """Calculate values needed for many d_delay_d_param functions""" # TODO: Should delay not have units of u.second? - delay = self._parent.delay(toas) + # @abhisrkckl: Commenting out unused computation + # This was slowing things significantly. + # delay = self._parent.delay(toas) + + if include == "all": + include = [ + "epoch", + "ssb_obs_r", + "ssb_obs_z", + "ssb_obs_xy", + "ssb_obs_x", + "ssb_obs_y", + "in_psr_obs", + "earth_dec", + "earth_ra", + ] # TODO: tbl['tdbld'].quantity should have units of u.day # NOTE: Do we need to include the delay here? tbl = toas.table + rd = {"epoch": tbl["tdbld"].quantity * u.day} + # Distance from SSB to observatory, and from SSB to psr ssb_obs = tbl["ssb_obs_pos"].quantity - ssb_psr = self.ssb_to_psb_xyz_ICRS(epoch=np.array(rd["epoch"])) + + ssb_obs_value = ssb_obs.value + ssb_obs_unit = ssb_obs.unit # Cartesian coordinates, and derived quantities - rd["ssb_obs_r"] = np.sqrt(np.sum(ssb_obs**2, axis=1)) - rd["ssb_obs_z"] = ssb_obs[:, 2] - rd["ssb_obs_xy"] = np.sqrt(ssb_obs[:, 0] ** 2 + ssb_obs[:, 1] ** 2) - rd["ssb_obs_x"] = ssb_obs[:, 0] - rd["ssb_obs_y"] = ssb_obs[:, 1] - rd["in_psr_obs"] = np.sum(ssb_obs * ssb_psr, axis=1) + if "ssb_obs_r" in include: + rd["ssb_obs_r"] = ( + np.sqrt(np.sum(ssb_obs_value**2, axis=1)) << ssb_obs_unit + ) + + if "ssb_obs_z" in include: + rd["ssb_obs_z"] = ssb_obs[:, 2] + + if "ssb_obs_xy" in include or "earth_dec" in include: + rd["ssb_obs_xy"] = ( + np.sqrt(ssb_obs_value[:, 0] ** 2 + ssb_obs_value[:, 1] ** 2) + << ssb_obs_unit + ) + + if "earth_dec" in include: + rd["earth_dec"] = np.arctan2(ssb_obs[:, 2], rd["ssb_obs_xy"]) + + if "ssb_obs_x" in include: + rd["ssb_obs_x"] = ssb_obs[:, 0] + + if "ssb_obs_y" in include: + rd["ssb_obs_y"] = ssb_obs[:, 1] + + if "in_psr_obs" in include: + ssb_psr = self.ssb_to_psb_xyz_ICRS(epoch=np.array(rd["epoch"])) + rd["in_psr_obs"] = np.sum(ssb_obs * ssb_psr, axis=1) - # Earth right ascension and declination - rd["earth_dec"] = np.arctan2(rd["ssb_obs_z"], rd["ssb_obs_xy"]) - rd["earth_ra"] = np.arctan2(rd["ssb_obs_y"], rd["ssb_obs_x"]) + if "earth_ra" in include: + rd["earth_ra"] = np.arctan2(ssb_obs[:, 1], ssb_obs[:, 0]) return rd @@ -200,17 +238,16 @@ def d_delay_astrometry_d_PX(self, toas, param="", acc_delay=None): t_d = 0.5 * px_r * delta'/ c, and delta = delta' * px_r / (1 AU) """ - rd = self.get_d_delay_quantities(toas) + rd = self.get_d_delay_quantities(toas, include=["ssb_obs_r", "in_psr_obs"]) - px_r = np.sqrt(rd["ssb_obs_r"] ** 2 - rd["in_psr_obs"] ** 2) - dd_dpx = 0.5 * (px_r**2 / (u.AU * const.c)) * (u.mas / u.radian) + px_r_2 = rd["ssb_obs_r"] ** 2 - rd["in_psr_obs"] ** 2 + dd_dpx = 0.5 * px_r_2 / (u.AU * const.c * u.radian) - # We want to return sec / mas - return dd_dpx.decompose(u.si.bases) / u.mas + return dd_dpx.to(u.s / u.mas) def d_delay_astrometry_d_POSEPOCH(self, toas, param="", acc_delay=None): """Calculate the derivative wrt POSEPOCH""" - pass + raise NotImplementedError def change_posepoch(self, new_epoch): """Change POSEPOCH to a new value and update the position accordingly. @@ -405,7 +442,7 @@ def d_delay_astrometry_d_RAJ(self, toas, param="", acc_delay=None): delay = r*[cos(de)*cos(dp)*cos(ae-aa)+sin(de)*sin(dp)]/c """ - rd = self.get_d_delay_quantities(toas) + rd = self.get_d_delay_quantities(toas, ["earth_dec", "earth_ra", "ssb_obs_r"]) psr_ra = self.RAJ.quantity psr_dec = self.DECJ.quantity @@ -415,14 +452,14 @@ def d_delay_astrometry_d_RAJ(self, toas, param="", acc_delay=None): ) dd_draj = rd["ssb_obs_r"] * geom / (const.c * u.radian) - return dd_draj.decompose(u.si.bases) + return dd_draj.to(u.s / self.RAJ.units) def d_delay_astrometry_d_DECJ(self, toas, param="", acc_delay=None): """Calculate the derivative wrt DECJ Definitions as in d_delay_d_RAJ """ - rd = self.get_d_delay_quantities(toas) + rd = self.get_d_delay_quantities(toas, ["earth_dec", "earth_ra", "ssb_obs_r"]) psr_ra = self.RAJ.quantity psr_dec = self.DECJ.quantity @@ -432,7 +469,7 @@ def d_delay_astrometry_d_DECJ(self, toas, param="", acc_delay=None): ) - np.sin(rd["earth_dec"]) * np.cos(psr_dec) dd_ddecj = rd["ssb_obs_r"] * geom / (const.c * u.radian) - return dd_ddecj.decompose(u.si.bases) + return dd_ddecj.to(u.s / self.DECJ.units) def d_delay_astrometry_d_PMRA(self, toas, param="", acc_delay=None): """Calculate the derivative wrt PMRA @@ -440,18 +477,19 @@ def d_delay_astrometry_d_PMRA(self, toas, param="", acc_delay=None): Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for the pulsar RA """ - rd = self.get_d_delay_quantities(toas) + rd = self.get_d_delay_quantities( + toas, ["epoch", "earth_dec", "earth_ra", "ssb_obs_r"] + ) psr_ra = self.RAJ.quantity te = rd["epoch"] - self.POSEPOCH.quantity.tdb.mjd_long * u.day geom = np.cos(rd["earth_dec"]) * np.sin(psr_ra - rd["earth_ra"]) - deriv = rd["ssb_obs_r"] * geom * te / (const.c * u.radian) - dd_dpmra = deriv * u.mas / u.year + dd_dpmra = rd["ssb_obs_r"] * geom * te / (const.c * u.radian) # We want to return sec / (mas / yr) - return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year) + return dd_dpmra.to(u.s / self.PMRA.units) def d_delay_astrometry_d_PMDEC(self, toas, param="", acc_delay=None): """Calculate the derivative wrt PMDEC @@ -459,7 +497,9 @@ def d_delay_astrometry_d_PMDEC(self, toas, param="", acc_delay=None): Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for the pulsar DEC """ - rd = self.get_d_delay_quantities(toas) + rd = self.get_d_delay_quantities( + toas, ["epoch", "earth_dec", "earth_ra", "ssb_obs_r"] + ) psr_ra = self.RAJ.quantity psr_dec = self.DECJ.quantity @@ -469,11 +509,10 @@ def d_delay_astrometry_d_PMDEC(self, toas, param="", acc_delay=None): psr_ra - rd["earth_ra"] ) - np.cos(psr_dec) * np.sin(rd["earth_dec"]) - deriv = rd["ssb_obs_r"] * geom * te / (const.c * u.radian) - dd_dpmdec = deriv * u.mas / u.year + dd_dpmdec = rd["ssb_obs_r"] * geom * te / (const.c * u.radian) # We want to return sec / (mas / yr) - return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year) + return dd_dpmdec.to(u.s / self.PMDEC.units) def change_posepoch(self, new_epoch): """Change POSEPOCH to a new value and update the position accordingly. From 1e84f06196fb8d81ce453d8a8ddb67b957e2ed8e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 13:51:28 -0500 Subject: [PATCH 03/19] comment --- src/pint/models/astrometry.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 2eeb6e9ef..ce85d312b 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -63,6 +63,16 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): If epochs (MJD) are given, proper motion is included in the calculation. """ # TODO: would it be better for this to return a 6-vector (pos, vel)? + + # @abhisrkckl: The following operations are stupid slow for some reason. + # This is somewhat uunderstandable when the proper motion is non-zero, and + # we'll have to do this for every epoch. But this is slow even when there is + # no proper motion. The underlying operations should just be some trig functions + # and simple 3x3 matrix operations; so I don't understand why this has to be so slow. + # Anyway, I am not going down this rabbit hole right now. This rant is here to remind + # people that there is a performance bottleneck here. Same problem is there with a similar + # function in AstrometryEcliptic. + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): @@ -173,7 +183,6 @@ def get_d_delay_quantities(self, toas, include="all"): # Distance from SSB to observatory, and from SSB to psr ssb_obs = tbl["ssb_obs_pos"].quantity - ssb_obs_value = ssb_obs.value ssb_obs_unit = ssb_obs.unit @@ -382,6 +391,7 @@ def get_psr_coords(self, epoch=None): obstime=self.POSEPOCH.quantity, frame=coords.ICRS, ) + newepoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) From 4188453ee1733e8ec01d475bf4c3d74ff3ecda83 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 13:58:04 -0500 Subject: [PATCH 04/19] update --- src/pint/residuals.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/pint/residuals.py b/src/pint/residuals.py index a7876a4e9..52c75d2ea 100644 --- a/src/pint/residuals.py +++ b/src/pint/residuals.py @@ -547,7 +547,7 @@ def _calc_gls_chi2(self): return chi2 - def calc_chi2(self): + def calc_chi2(self, update=False): """Return the weighted chi-squared for the model and toas. If the errors on the TOAs are independent this is a straightforward @@ -566,6 +566,9 @@ def calc_chi2(self): a minimizer for example - may need to be checked to confirm that they correctly return infinity. """ + if update: + self.update() + if self.model.has_correlated_errors: return self._calc_gls_chi2() else: @@ -755,7 +758,10 @@ def calc_resids(self): resids -= resids.mean() return resids - def calc_chi2(self): + def calc_chi2(self, update=False): + if update: + self.update() + data_errors = self.get_data_error() if (data_errors == 0.0).any(): return np.inf @@ -952,7 +958,7 @@ def chi2(self): assert self._chi2 is not None return self._chi2 - def calc_chi2(self, full_cov=False): + def calc_chi2(self, full_cov=False, update=False): """Return the weighted chi-squared for the model and toas. If the errors on the TOAs are independent this is a straightforward @@ -971,6 +977,10 @@ def calc_chi2(self, full_cov=False): a minimizer for example - may need to be checked to confirm that they correctly return infinity. """ + if update: + self.toa.update() + self.dm.update() + log.debug("Using wideband GLS fitter to compute residual chi2") # Use GLS but don't actually fit from pint.fitter import WidebandTOAFitter From 1820ba680901ee7a0f59fb8dbdf6d663cbedbb99 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 15:17:05 -0500 Subject: [PATCH 05/19] cached _ssb_to_psb_xyz_ICRS --- src/pint/models/astrometry.py | 58 ++++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 4 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index ce85d312b..958dfcd11 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -8,10 +8,10 @@ import astropy.units as u import numpy as np from astropy.time import Time - +from erfa import ErfaWarning from loguru import logger as log +from memoization import cached -from erfa import ErfaWarning from pint import ls from pint.models.parameter import ( AngleParameter, @@ -23,7 +23,6 @@ from pint.pulsar_ecliptic import OBL, PulsarEcliptic from pint.utils import add_dummy_distance, remove_dummy_distance - astropy_version = sys.modules["astropy"].__version__ mas_yr = u.mas / u.yr @@ -72,8 +71,18 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # Anyway, I am not going down this rabbit hole right now. This rant is here to remind # people that there is a performance bottleneck here. Same problem is there with a similar # function in AstrometryEcliptic. + # As a stopgap measure, I am moving this computation to a cached pure function below. - return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + + return _ssb_to_psb_xyz_ICRS( + self.RAJ.quantity, + self.DECJ.quantity, + self.PMRA.quantity, + self.PMDEC.quantity, + self.POSEPOCH.quantity, + np.array(epoch), + ) def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates. @@ -1124,3 +1133,44 @@ def as_ICRS(self, epoch=None): m_eq.PMDEC.frozen = self.PMELAT.frozen return m_eq + + +@cached +def _ssb_to_psb_xyz_ICRS(ra, dec, pmra, pmdec, posepoch, epoch): + """A cached version of `AstrometryEquatorial.ssb_to_psb_xyz_ICRS`. + This is used to avoid repeated calls since this function is expensive. + """ + if epoch is None or (pmra.value == 0.0 and pmdec.value == 0.0): + coords_as_ICRS = coords.SkyCoord( + ra=ra, + dec=dec, + pm_ra_cosdec=pmra, + pm_dec=pmdec, + obstime=posepoch, + frame=coords.ICRS, + ) + else: + newepoch = ( + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") + ) + position_now = add_dummy_distance( + coords.SkyCoord( + ra=ra, + dec=dec, + pm_ra_cosdec=pmra, + pm_dec=pmdec, + obstime=posepoch, + frame=coords.ICRS, + ) + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ErfaWarning) + # for the most part the dummy distance should remove any potential erfa warnings + # but for some very large proper motions that does not quite work + # so we catch the warnings + position_then = position_now.apply_space_motion(new_obstime=newepoch) + position_then = remove_dummy_distance(position_then) + + coords_as_ICRS = position_then + + return coords_as_ICRS.cartesian.xyz.transpose() From 3c4fe07683d596a1bc5b8beaa2d343ad7b7e2322 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 15:36:26 -0500 Subject: [PATCH 06/19] _ssb_to_psb_xyz_ICRS_eq --- src/pint/models/astrometry.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 958dfcd11..43a08734a 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -73,16 +73,17 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # function in AstrometryEcliptic. # As a stopgap measure, I am moving this computation to a cached pure function below. - # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - - return _ssb_to_psb_xyz_ICRS( - self.RAJ.quantity, - self.DECJ.quantity, - self.PMRA.quantity, - self.PMDEC.quantity, - self.POSEPOCH.quantity, - np.array(epoch), - ) + if isinstance(self, AstrometryEcliptic): + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + else: + return _ssb_to_psb_xyz_ICRS_eq( + self.RAJ.quantity, + self.DECJ.quantity, + self.PMRA.quantity, + self.PMDEC.quantity, + self.POSEPOCH.quantity, + np.array(epoch), + ) def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates. @@ -1136,7 +1137,7 @@ def as_ICRS(self, epoch=None): @cached -def _ssb_to_psb_xyz_ICRS(ra, dec, pmra, pmdec, posepoch, epoch): +def _ssb_to_psb_xyz_ICRS_eq(ra, dec, pmra, pmdec, posepoch, epoch): """A cached version of `AstrometryEquatorial.ssb_to_psb_xyz_ICRS`. This is used to avoid repeated calls since this function is expensive. """ From 2807721711e1ab806847ae5e97241e08bbeae01a Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 18:19:31 -0500 Subject: [PATCH 07/19] ssb_to_psb_xyz_ICRS --- src/pint/models/astrometry.py | 118 +++++++++++++++++++--------------- 1 file changed, 65 insertions(+), 53 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 43a08734a..7750e269f 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -10,7 +10,6 @@ from astropy.time import Time from erfa import ErfaWarning from loguru import logger as log -from memoization import cached from pint import ls from pint.models.parameter import ( @@ -73,17 +72,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # function in AstrometryEcliptic. # As a stopgap measure, I am moving this computation to a cached pure function below. - if isinstance(self, AstrometryEcliptic): - return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - else: - return _ssb_to_psb_xyz_ICRS_eq( - self.RAJ.quantity, - self.DECJ.quantity, - self.PMRA.quantity, - self.PMDEC.quantity, - self.POSEPOCH.quantity, - np.array(epoch), - ) + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates. @@ -439,6 +428,35 @@ def coords_as_GAL(self, epoch=None): pos_icrs = self.get_psr_coords(epoch=epoch) return pos_icrs.transform_to(coords.Galactic) + def ssb_to_psb_xyz_ICRS(self, epoch=None): + """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. + + If epochs (MJD) are given, proper motion is included in the calculation. + """ + + # There is a the version of this function in `Astrometry` that gives the same result. + # This override is significantly faster. + + ra0, dec0 = self.RAJ.quantity, self.DECJ.quantity + pmra, pmdec = self.PMRA.quantity, self.PMDEC.quantity + posepoch = self.POSEPOCH.quantity + + if epoch is None or (pmra == 0 and pmdec == 0): + ra, dec = ra0, dec0 + else: + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="mjd") + ) + dt = epoch - posepoch + ra = ra0 + pmra * dt / np.cos(dec0) + dec = dec0 + pmdec * dt + + return u.Quantity( + [np.cos(dec) * np.cos(ra), np.cos(dec) * np.sin(ra), np.sin(dec)] + ).T + def get_params_as_ICRS(self): return { "RAJ": self.RAJ.quantity, @@ -813,6 +831,41 @@ def coords_as_ECL(self, epoch=None, ecl=None): pos_ecl = pos_ecl.transform_to(PulsarEcliptic(ecl=ecl)) return pos_ecl + def ssb_to_psb_xyz_ICRS(self, epoch=None): + """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. + + If epochs (MJD) are given, proper motion is included in the calculation. + """ + + # There is a the version of this function in `Astrometry` that gives the same result. + # This override is significantly faster. + + elong0, elat0 = self.ELONG.quantity, self.ELAT.quantity + pmelong, pmelat = self.PMELONG.quantity, self.PMELAT.quantity + posepoch = self.POSEPOCH.quantity + obl = OBL[self.ECL.value] + + if epoch is None or (pmelong == 0 and pmelat == 0): + elong, elat = elong0, elat0 + else: + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="mjd") + ) + dt = epoch - posepoch + elong = elong0 + pmelong * dt / np.cos(elat0) + elat = elat0 + pmelat * dt + + xhat_ecl = u.Quantity( + [np.cos(elat) * np.cos(elong), np.cos(elat) * np.sin(elong), np.sin(elat)] + ) + R = np.array( + [[1, 0, 0], [0, np.cos(obl), -np.sin(obl)], [0, np.sin(obl), np.cos(obl)]] + ) + + return (R @ xhat_ecl).T + def get_d_delay_quantities_ecliptical(self, toas): """Calculate values needed for many d_delay_d_param functions.""" # TODO: Move all these calculations in a separate class for elegance @@ -1134,44 +1187,3 @@ def as_ICRS(self, epoch=None): m_eq.PMDEC.frozen = self.PMELAT.frozen return m_eq - - -@cached -def _ssb_to_psb_xyz_ICRS_eq(ra, dec, pmra, pmdec, posepoch, epoch): - """A cached version of `AstrometryEquatorial.ssb_to_psb_xyz_ICRS`. - This is used to avoid repeated calls since this function is expensive. - """ - if epoch is None or (pmra.value == 0.0 and pmdec.value == 0.0): - coords_as_ICRS = coords.SkyCoord( - ra=ra, - dec=dec, - pm_ra_cosdec=pmra, - pm_dec=pmdec, - obstime=posepoch, - frame=coords.ICRS, - ) - else: - newepoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") - ) - position_now = add_dummy_distance( - coords.SkyCoord( - ra=ra, - dec=dec, - pm_ra_cosdec=pmra, - pm_dec=pmdec, - obstime=posepoch, - frame=coords.ICRS, - ) - ) - with warnings.catch_warnings(): - warnings.simplefilter("ignore", ErfaWarning) - # for the most part the dummy distance should remove any potential erfa warnings - # but for some very large proper motions that does not quite work - # so we catch the warnings - position_then = position_now.apply_space_motion(new_obstime=newepoch) - position_then = remove_dummy_distance(position_then) - - coords_as_ICRS = position_then - - return coords_as_ICRS.cartesian.xyz.transpose() From bda77a10cb544ef39f5f3731b6e0224cd2ecfc12 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 18:26:36 -0500 Subject: [PATCH 08/19] gitignore --- tests/datafile/.gitignore | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 tests/datafile/.gitignore diff --git a/tests/datafile/.gitignore b/tests/datafile/.gitignore new file mode 100644 index 000000000..242803bc9 --- /dev/null +++ b/tests/datafile/.gitignore @@ -0,0 +1,4 @@ +_test_pintempo.out +par*_*.par +fake_toas.tim +*.converted.par \ No newline at end of file From f4fa8cef8f6a8120eb55d4d58ddf13105e80d30b Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 18:38:54 -0500 Subject: [PATCH 09/19] test_fast_ssb_to_psb_xyz_ICRS --- tests/test_astrometry.py | 46 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/test_astrometry.py b/tests/test_astrometry.py index 152cd7c27..58ca7061d 100644 --- a/tests/test_astrometry.py +++ b/tests/test_astrometry.py @@ -7,6 +7,7 @@ from astropy.coordinates import Latitude, Longitude from pint.models import get_model +from pint.models.astrometry import AstrometryEcliptic, AstrometryEquatorial from pinttestdata import datadir @@ -81,3 +82,48 @@ def test_pm_acquires_posepoch_equatorial(): def test_pm_one_set_other_not(): m = get_model(StringIO("\n".join([par_basic_equatorial, "PMRA 7"]))) assert m.POSEPOCH.quantity == m.PEPOCH.quantity + + +def test_fast_ssb_to_psb_xyz_ICRS(): + m = get_model(StringIO(par_basic_ecliptic)) + ts = np.linspace(55000, 55500, 10) + + assert np.allclose( + m.ssb_to_psb_xyz_ICRS(), + super( + AstrometryEcliptic, m.components["AstrometryEcliptic"] + ).ssb_to_psb_xyz_ICRS(), + ) + + m.POSEPOCH.quantity = m.PEPOCH.quantity + m.PMELAT.value = 1e-3 + m.PMELONG.value = 1e-3 + + assert np.allclose( + m.ssb_to_psb_xyz_ICRS(epoch=ts), + super( + AstrometryEcliptic, m.components["AstrometryEcliptic"] + ).ssb_to_psb_xyz_ICRS(epoch=ts), + ) + + # ========= + + m = get_model(StringIO(par_basic_equatorial)) + + assert np.allclose( + m.ssb_to_psb_xyz_ICRS(), + super( + AstrometryEquatorial, m.components["AstrometryEquatorial"] + ).ssb_to_psb_xyz_ICRS(), + ) + + m.POSEPOCH.quantity = m.PEPOCH.quantity + m.PMRA.value = 1e-3 + m.PMDEC.value = 1e-3 + + assert np.allclose( + m.ssb_to_psb_xyz_ICRS(epoch=ts), + super( + AstrometryEquatorial, m.components["AstrometryEquatorial"] + ).ssb_to_psb_xyz_ICRS(epoch=ts), + ) From bf8d86fc85c2de1ec9a1b0f475e727930367da8e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 18:45:05 -0500 Subject: [PATCH 10/19] CHANGELOG --- CHANGELOG-unreleased.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 1bebb820a..1f716e837 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -12,12 +12,15 @@ the released changes. - `WAVE` parameters can be added to a `Wave` model with `add_wave_component()` in `wave.py` - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). +- Optionally avoid certain computations in `Astrometry.get_d_delay_quantities()` method. ### Added - Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters - `Parameter.as_latex` method for latex representation of a parameter. - `pint.output.publish` module and `pintpublish` script for generating publication (LaTeX) output. - Added radial velocity methods for binary models - Added `DMWaveX` model (Fourier representation of DM noise) +- Faster implementations of `ssb_to_psb_xyz_ICRS()` method in `AstrometryEquatorial` and `AstrometryEcliptic` without using `astropy.coordinates`. Original version in `Astrometry` is retained for comparison. +- `update` option in `Residuals.calc_chi2` and `WidebandTOAResiduals.calc_chi2` ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. From d3c7be944c25d80139f2cc1c7bad675a11bfcc57 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sat, 23 Sep 2023 19:39:40 -0500 Subject: [PATCH 11/19] solar_system_shapiro --- src/pint/models/solar_system_shapiro.py | 58 ++++++++++++------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/src/pint/models/solar_system_shapiro.py b/src/pint/models/solar_system_shapiro.py index 1cd850065..860abbf1b 100644 --- a/src/pint/models/solar_system_shapiro.py +++ b/src/pint/models/solar_system_shapiro.py @@ -2,7 +2,7 @@ import astropy.constants as const import astropy.units as u -import numpy +import numpy as np from loguru import logger as log from pint import ( @@ -69,8 +69,8 @@ def ss_obj_shapiro_delay(obj_pos, psr_dir, T_obj): T_obj : mass of object in seconds (GM/c^3) """ # TODO: numpy.sum currently loses units in some cases... - r = (numpy.sqrt(numpy.sum(obj_pos**2, axis=1))) * obj_pos.unit - rcostheta = numpy.sum(obj_pos * psr_dir, axis=1) + r = (np.sqrt(np.sum(obj_pos**2, axis=1))) * obj_pos.unit + rcostheta = np.sum(obj_pos * psr_dir, axis=1) # This is the 2nd to last term from Eqn 4.6 in Backer & # Hellings, ARAA, 1986 with gamma = 1 (as defined by GR). We # have the opposite sign of the cos(theta) term, since our @@ -79,7 +79,7 @@ def ss_obj_shapiro_delay(obj_pos, psr_dir, T_obj): # pulsar (as described after Eqn 4.3 in the paper). # See also https://en.wikipedia.org/wiki/Shapiro_time_delay # where \Delta t = \frac{2GM}{c^3}\log(1-\vec{R}\cdot\vec{x}) - return -2.0 * T_obj * numpy.log((r - rcostheta) / const.au).value + return -2.0 * T_obj * np.log((r - rcostheta) / const.au).value def solar_system_shapiro_delay(self, toas, acc_delay=None): """ @@ -96,28 +96,28 @@ def solar_system_shapiro_delay(self, toas, acc_delay=None): """ # Start out with 0 delay with units of seconds tbl = toas.table - delay = numpy.zeros(len(tbl)) - for key, grp in toas.get_obs_groups(): - if key.lower() == "barycenter": - log.debug("Skipping Shapiro delay for Barycentric TOAs") - continue - psr_dir = self._parent.ssb_to_psb_xyz_ICRS( - epoch=tbl[grp]["tdbld"].astype(numpy.float64) - ) - delay[grp] += self.ss_obj_shapiro_delay( - tbl[grp]["obs_sun_pos"], psr_dir, self._ss_mass_sec["sun"] - ) - try: - if self.PLANET_SHAPIRO.value: - for pl in ("jupiter", "saturn", "venus", "uranus", "neptune"): - delay[grp] += self.ss_obj_shapiro_delay( - tbl[grp][f"obs_{pl}_pos"], - psr_dir, - self._ss_mass_sec[pl], - ) - except KeyError as e: - raise KeyError( - "Planet positions not found when trying to compute Solar System Shapiro delay. " - "Make sure that you include `planets=True` in your `get_TOAs()` call, or use `get_model_and_toas()`." - ) from e - return delay * u.second + delay = np.zeros(len(tbl)) + + non_bary_mask = toas.get_obss() != "barycenter" + + tbl_mask = tbl if np.all(non_bary_mask) else tbl[non_bary_mask] + + psr_dir = self._parent.ssb_to_psb_xyz_ICRS(epoch=tbl_mask["tdbld"]) + delay[non_bary_mask] += self.ss_obj_shapiro_delay( + tbl_mask["obs_sun_pos"], psr_dir, self._ss_mass_sec["sun"] + ) + try: + if self.PLANET_SHAPIRO.value: + for pl in ("jupiter", "saturn", "venus", "uranus", "neptune"): + delay[non_bary_mask] += self.ss_obj_shapiro_delay( + tbl_mask[f"obs_{pl}_pos"], + psr_dir, + self._ss_mass_sec[pl], + ) + except KeyError as e: + raise KeyError( + "Planet positions not found when trying to compute Solar System Shapiro delay. " + "Make sure that you include `planets=True` in your `get_TOAs()` call, or use `get_model_and_toas()`." + ) from e + + return u.Quantity(delay, u.s) From c22a006c649cf17bc67cdab4f15e6c22e48bbd58 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sun, 24 Sep 2023 16:17:31 -0500 Subject: [PATCH 12/19] ssb_to_psb_xyz_ICRS --- src/pint/models/astrometry.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 7750e269f..b4e81bb5d 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -453,9 +453,11 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): ra = ra0 + pmra * dt / np.cos(dec0) dec = dec0 + pmdec * dt - return u.Quantity( - [np.cos(dec) * np.cos(ra), np.cos(dec) * np.sin(ra), np.sin(dec)] - ).T + ra = ra.to_value(u.rad) + dec = dec.to_value(u.rad) + + cos_dec = np.cos(dec) + return u.Quantity([cos_dec * np.cos(ra), cos_dec * np.sin(ra), np.sin(dec)]).T def get_params_as_ICRS(self): return { From b46f613151766fc37dae9a022f7f576e45e68261 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Sun, 24 Sep 2023 17:19:15 -0500 Subject: [PATCH 13/19] doc --- CHANGELOG-unreleased.md | 1 + src/pint/models/astrometry.py | 9 +++++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 1f716e837..56f9dc5a5 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -13,6 +13,7 @@ the released changes. - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). - Optionally avoid certain computations in `Astrometry.get_d_delay_quantities()` method. +- Removed an unnecessary for loop in `solar_system_shapiro_delay` ### Added - Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters - `Parameter.as_latex` method for latex representation of a parameter. diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index b4e81bb5d..b80b2cdc4 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -63,14 +63,15 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # TODO: would it be better for this to return a 6-vector (pos, vel)? # @abhisrkckl: The following operations are stupid slow for some reason. - # This is somewhat uunderstandable when the proper motion is non-zero, and + # This is somewhat understandable when the proper motion is non-zero, and # we'll have to do this for every epoch. But this is slow even when there is # no proper motion. The underlying operations should just be some trig functions # and simple 3x3 matrix operations; so I don't understand why this has to be so slow. # Anyway, I am not going down this rabbit hole right now. This rant is here to remind - # people that there is a performance bottleneck here. Same problem is there with a similar - # function in AstrometryEcliptic. - # As a stopgap measure, I am moving this computation to a cached pure function below. + # people that there is a performance bottleneck here. + # + # As a stopgap measure, I am overriding this function in the subclasses without + # using `astropy.coordinates`. return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() From 02af78040139bc345f3fdbbb57d765d17563e5a1 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 25 Sep 2023 09:19:42 -0500 Subject: [PATCH 14/19] docs --- src/pint/models/astrometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index b80b2cdc4..f178777ac 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -62,11 +62,11 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): """ # TODO: would it be better for this to return a 6-vector (pos, vel)? - # @abhisrkckl: The following operations are stupid slow for some reason. + # @abhisrkckl: The following operations are very slow for some reason. # This is somewhat understandable when the proper motion is non-zero, and # we'll have to do this for every epoch. But this is slow even when there is # no proper motion. The underlying operations should just be some trig functions - # and simple 3x3 matrix operations; so I don't understand why this has to be so slow. + # and a 3x3 matrix multiplication; so I don't understand why this has to be so slow. # Anyway, I am not going down this rabbit hole right now. This rant is here to remind # people that there is a performance bottleneck here. # From b5068254ea5b69f795d01427353479239edae8d0 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 25 Sep 2023 09:35:13 -0500 Subject: [PATCH 15/19] sun_angle --- src/pint/models/astrometry.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index f178777ac..672a8bd45 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -115,15 +115,17 @@ def sun_angle(self, toas, heliocenter=True, also_distance=False): """ tbl = toas.table - if heliocenter: - osv = tbl["obs_sun_pos"].quantity.copy() - else: - osv = -tbl["ssb_obs_pos"].quantity.copy() + osv = ( + tbl["obs_sun_pos"].quantity if heliocenter else -tbl["ssb_obs_pos"].quantity + ) + psr_vec = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"]) r = (osv**2).sum(axis=1) ** 0.5 - osv /= r[:, None] - cos = (osv * psr_vec).sum(axis=1) - return (np.arccos(cos), r) if also_distance else np.arccos(cos) + + osvn = osv / r[:, None] + + cos_sa = (osvn * psr_vec).sum(axis=1) + return (np.arccos(cos_sa), r) if also_distance else np.arccos(cos_sa) def barycentric_radio_freq(self, toas): raise NotImplementedError From 007df395076d682922e08a627671a9962db946cf Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 25 Sep 2023 10:26:52 -0500 Subject: [PATCH 16/19] to_value --- src/pint/fitter.py | 22 +++++++++++----------- src/pint/pint_matrix.py | 4 +++- src/pint/plot_utils.py | 4 ++-- src/pint/residuals.py | 4 ++-- src/pint/toa.py | 16 ++++++++-------- src/pint/utils.py | 2 +- 6 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 61bcd9532..ef58758d4 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -531,7 +531,7 @@ def get_derived_params(self): tasc = ufloat( # This is a time in MJD self.model.TASC.quantity.mjd, - self.model.TASC.uncertainty.to(u.d).value, + self.model.TASC.uncertainty.to_value(u.d), ) if hasattr(self.model, "PB") and self.model.PB.value is not None: pb = self.model.PB.as_ufloat(u.d) @@ -578,8 +578,8 @@ def get_derived_params(self): 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, + pb.to_value(u.s), + pberr.to_value(u.s), ) a1 = self.model.A1.as_ufloat(pint.ls) # This is the mass function, done explicitly so that we get @@ -1289,8 +1289,8 @@ def step(self): toas=self.fitter.toas, incfrozen=False, incoffset=True ) # Get residuals and TOA uncertainties in seconds - Nvec = self.model.scaled_toa_uncertainty(self.fitter.toas).to(u.s).value - scaled_resids = self.resids.time_resids.to(u.s).value / Nvec + Nvec = self.model.scaled_toa_uncertainty(self.fitter.toas).to_value(u.s) + scaled_resids = self.resids.time_resids.to_value(u.s) / Nvec # "Whiten" design matrix and residuals by dividing by uncertainties M = M / Nvec.reshape((-1, 1)) @@ -1442,7 +1442,7 @@ def step(self): covariance_matrix_labels = [covariance_matrix_labels] * 2 self.parameter_covariance_matrix_labels = covariance_matrix_labels - residuals = self.resids.time_resids.to(u.s).value + residuals = self.resids.time_resids.to_value(u.s) # get any noise design matrices and weight vectors if not self.full_cov: @@ -1470,7 +1470,7 @@ def step(self): phiinv /= norm**2 # Why are we scaling residuals by the *square* of the uncertainty? Nvec = ( - self.model.scaled_toa_uncertainty(self.fitter.toas).to(u.s).value ** 2 + self.model.scaled_toa_uncertainty(self.fitter.toas).to_value(u.s) ** 2 ) cinv = 1 / Nvec mtcm = np.dot(M.T, cinv[:, None] * M) @@ -1989,8 +1989,8 @@ def fit_toas(self, maxiter=1, threshold=None, debug=False): M, params, units = self.get_designmatrix() # Get residuals and TOA uncertainties in seconds self.update_resids() - residuals = self.resids.time_resids.to(u.s).value - Nvec = self.model.scaled_toa_uncertainty(self.toas).to(u.s).value + residuals = self.resids.time_resids.to_value(u.s) + Nvec = self.model.scaled_toa_uncertainty(self.toas).to_value(u.s) # "Whiten" design matrix and residuals by dividing by uncertainties M = M / Nvec.reshape((-1, 1)) @@ -2161,7 +2161,7 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False, debug=False): if i == 0: # Why is this here? self.update_resids() - residuals = self.resids.time_resids.to(u.s).value + residuals = self.resids.time_resids.to_value(u.s) # get any noise design matrices and weight vectors if not full_cov: @@ -2188,7 +2188,7 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False, debug=False): else: phiinv /= norm**2 - Nvec = self.model.scaled_toa_uncertainty(self.toas).to(u.s).value ** 2 + Nvec = self.model.scaled_toa_uncertainty(self.toas).to_value(u.s) ** 2 cinv = 1 / Nvec mtcm = np.dot(M.T, cinv[:, None] * M) mtcm += np.diag(phiinv) diff --git a/src/pint/pint_matrix.py b/src/pint/pint_matrix.py index 1900249ff..ad6aac88b 100644 --- a/src/pint/pint_matrix.py +++ b/src/pint/pint_matrix.py @@ -478,7 +478,9 @@ def __call__(self, data, model, derivative_params, offset=True, offset_padding=1 else: param_unit = getattr(model, param).units # Since this is the phase derivative, we know the quantity unit. - q = deriv_func(data, delay, param).to(u.Unit("") / param_unit) + q = deriv_func(data, delay, param).to( + u.dimensionless_unscaled / param_unit + ) # NOTE Here we have negative sign here. Since in pulsar timing # the residuals are calculated as (Phase - int(Phase)), which is different diff --git a/src/pint/plot_utils.py b/src/pint/plot_utils.py index e95dae99e..0d5b96426 100644 --- a/src/pint/plot_utils.py +++ b/src/pint/plot_utils.py @@ -45,7 +45,7 @@ def phaseogram( else: mjds = mjds_in - years = (mjds.to(u.d).value - 51544.0) / 365.25 + 2000.0 + years = (mjds.to_value(u.d) - 51544.0) / 365.25 + 2000.0 phss = phases + rotate phss[phss >= 1.0] -= 1.0 plt.figure(figsize=(width, 8)) @@ -133,7 +133,7 @@ def phaseogram_binned( else: mjds = mjds_in - years = (mjds.to(u.d).value - 51544.0) / 365.25 + 2000.0 + years = (mjds.to_value(u.d) - 51544.0) / 365.25 + 2000.0 phss = phases + rotate phss[phss >= 1.0] -= 1.0 plt.figure(figsize=(width, 8)) diff --git a/src/pint/residuals.py b/src/pint/residuals.py index 52c75d2ea..2977e9f19 100644 --- a/src/pint/residuals.py +++ b/src/pint/residuals.py @@ -506,7 +506,7 @@ def _calc_gls_chi2(self): If the system is singular, it uses singular value decomposition instead.""" self.update() - residuals = self.time_resids.to(u.s).value + residuals = self.time_resids.to_value(u.s) M = self.model.noise_model_designmatrix(self.toas) phi = self.model.noise_model_basis_weight(self.toas) @@ -520,7 +520,7 @@ def _calc_gls_chi2(self): M, norm = normalize_designmatrix(M, None) phiinv /= norm**2 - Nvec = self.model.scaled_toa_uncertainty(self.toas).to(u.s).value ** 2 + Nvec = self.model.scaled_toa_uncertainty(self.toas).to_value(u.s) ** 2 cinv = 1 / Nvec mtcm = np.dot(M.T, cinv[:, None] * M) mtcm += np.diag(phiinv) diff --git a/src/pint/toa.py b/src/pint/toa.py index 682ca60ce..2731a649e 100644 --- a/src/pint/toa.py +++ b/src/pint/toa.py @@ -632,7 +632,7 @@ def format_toa_line( freq = 0.0 * u.MHz flagstring = "" if dm != 0.0 * pint.dmu: - flagstring += "-dm {:.5f}".format(dm.to(pint.dmu).value) + flagstring += "-dm {:.5f}".format(dm.to_value(pint.dmu)) # Here I need to append any actual flags for flag in flags.keys(): v = flags[flag] @@ -651,9 +651,9 @@ def format_toa_line( obscode = obs.name out = "%s %f %s %.3f %s %s\n" % ( name, - freq.to(u.MHz).value, + freq.to_value(u.MHz), toa_str, - toaerr.to(u.us).value, + toaerr.to_value(u.us), alias_translation.get(obscode, obscode), flagstring, ) @@ -677,17 +677,17 @@ def format_toa_line( if dm != 0.0 * pint.dmu: out = obs.tempo_code + " %13s%9.3f%20s%9.2f %9.4f\n" % ( name, - freq.to(u.MHz).value, + freq.to_value(u.MHz), toa_str, - toaerr.to(u.us).value, - dm.to(pint.dmu).value, + toaerr.to_value(u.us), + dm.to_value(pint.dmu), ) else: out = obs.tempo_code + " %13s%9.3f%20s%9.2f\n" % ( name, - freq.to(u.MHz).value, + freq.to_value(u.MHz), toa_str, - toaerr.to(u.us).value, + toaerr.to_value(u.us), ) else: raise ValueError(f"Unknown TOA format ({format})") diff --git a/src/pint/utils.py b/src/pint/utils.py index 785ae36be..a1c505648 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -2476,7 +2476,7 @@ def divide_times(t, t0, offset=0.5): """ dt = t - t0 - values = (dt.to(u.yr).value + offset) // 1 + values = (dt.to_value(u.yr) + offset) // 1 indices = np.digitize(values, np.unique(values), right=True) return indices From 2186fbb751301ac5f5deef401846bf41dc7bbb90 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 25 Sep 2023 10:45:42 -0500 Subject: [PATCH 17/19] to_value --- src/pint/models/astrometry.py | 6 +++--- src/pint/models/binary_ell1.py | 2 +- src/pint/models/dispersion_model.py | 2 +- src/pint/models/fdjump.py | 4 ++-- src/pint/models/frequency_dependent.py | 2 +- src/pint/models/ifunc.py | 2 +- src/pint/models/noise_model.py | 16 ++++++++-------- src/pint/models/parameter.py | 14 +++++++------- src/pint/models/pulsar_binary.py | 2 +- .../stand_alone_psr_binaries/DDK_model.py | 18 +++++++++++------- .../stand_alone_psr_binaries/DD_model.py | 4 ++-- .../stand_alone_psr_binaries/ELL1H_model.py | 4 ++-- .../stand_alone_psr_binaries/ELL1_model.py | 12 ++++++------ .../stand_alone_psr_binaries/binary_generic.py | 8 ++++---- .../stand_alone_psr_binaries/binary_orbits.py | 10 +++++----- src/pint/models/timing_model.py | 2 +- src/pint/models/troposphere_delay.py | 4 ++-- 17 files changed, 58 insertions(+), 54 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 672a8bd45..df9d5fb78 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -144,7 +144,7 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): if np.any(c): L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"][c].astype(np.float64)) re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) - delay[c] = -re_dot_L.to(ls).value + delay[c] = -re_dot_L.to_value(ls) if self.PX.value != 0.0: L = (1.0 / self.PX.value) * u.kpc # TODO: np.sum currently loses units in some cases... @@ -153,8 +153,8 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): * tbl["ssb_obs_pos"].unit ** 2 ) delay[c] += ( - (0.5 * (re_sqr / L) * (1.0 - re_dot_L**2 / re_sqr)).to(ls).value - ) + 0.5 * (re_sqr / L) * (1.0 - re_dot_L**2 / re_sqr) + ).to_value(ls) return delay * u.second def get_d_delay_quantities(self, toas, include="all"): diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index 8e3c93b4e..29044049d 100644 --- a/src/pint/models/binary_ell1.py +++ b/src/pint/models/binary_ell1.py @@ -262,7 +262,7 @@ def change_binary_epoch(self, new_epoch): tasc_ld = self.TASC.quantity.tdb.mjd_long dt = (new_epoch.tdb.mjd_long - tasc_ld) * u.day d_orbits = dt / PB - PBDOT * dt**2 / (2.0 * PB**2) - n_orbits = np.round(d_orbits.to(u.Unit(""))) + n_orbits = np.round(d_orbits.to(u.dimensionless_unscaled)) if n_orbits == 0: return dt_integer_orbits = PB * n_orbits + PB * PBDOT * n_orbits**2 / 2.0 diff --git a/src/pint/models/dispersion_model.py b/src/pint/models/dispersion_model.py index 6ef81258e..f7eea8bc2 100644 --- a/src/pint/models/dispersion_model.py +++ b/src/pint/models/dispersion_model.py @@ -266,7 +266,7 @@ def d_dm_d_DMs( else: DMEPOCH = self.DMEPOCH.value dt = (toas["tdbld"] - DMEPOCH) * u.day - dt_value = (dt.to(u.yr)).value + dt_value = dt.to_value(u.yr) return taylor_horner(dt_value, dm_terms) * (self.DM.units / par.units) def change_dmepoch(self, new_epoch): diff --git a/src/pint/models/fdjump.py b/src/pint/models/fdjump.py index 9cec3354e..752cceda4 100644 --- a/src/pint/models/fdjump.py +++ b/src/pint/models/fdjump.py @@ -139,9 +139,9 @@ def get_freq_y(self, toas): freq = tbl["freq"] y = ( - np.log(freq.to(u.GHz).value) + np.log(freq.to_value(u.GHz)) if self.FDJUMPLOG.value - else freq.to(u.GHz).value + else freq.to_value(u.GHz) ) non_finite = np.invert(np.isfinite(y)) y[non_finite] = 0.0 diff --git a/src/pint/models/frequency_dependent.py b/src/pint/models/frequency_dependent.py index 3094a8a9a..3b55be380 100644 --- a/src/pint/models/frequency_dependent.py +++ b/src/pint/models/frequency_dependent.py @@ -85,7 +85,7 @@ def FD_delay(self, toas, acc_delay=None): def FD_delay_frequency(self, freq): """Compute the FD delay at an array of frequencies.""" FD_mapping = self.get_prefix_mapping_component("FD") - log_freq = np.log(freq.to(u.GHz).value) + log_freq = np.log(freq.to_value(u.GHz)) non_finite = np.invert(np.isfinite(log_freq)) log_freq[non_finite] = 0.0 FD_coeff = [ diff --git a/src/pint/models/ifunc.py b/src/pint/models/ifunc.py index 91a3b79cb..9ff0ce441 100644 --- a/src/pint/models/ifunc.py +++ b/src/pint/models/ifunc.py @@ -110,7 +110,7 @@ def ifunc_phase(self, toas, delays): # the MJDs(x) and offsets (y) of the interpolation points x, y = np.asarray([t.quantity for t in terms]).T # form barycentered times - ts = toas.table["tdbld"] - delays.to(u.day).value + ts = toas.table["tdbld"] - delays.to_value(u.day) times = np.zeros(len(ts)) # Determine what type of interpolation we are doing. diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 49036fd8a..0654e4149 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -160,7 +160,7 @@ def scale_toa_sigma(self, toas): return sigma_scaled def sigma_scaled_cov_matrix(self, toas): - scaled_sigma = self.scale_toa_sigma(toas).to(u.s).value ** 2 + scaled_sigma = self.scale_toa_sigma(toas).to_value(u.s) ** 2 return np.diag(scaled_sigma) @@ -337,7 +337,7 @@ def get_noise_basis(self, toas): A quantization matrix maps TOAs to observing epochs. """ tbl = toas.table - t = (tbl["tdbld"].quantity * u.day).to(u.s).value + t = (tbl["tdbld"].quantity * u.day).to_value(u.s) ecorrs = self.get_ecorrs() umats = [] for ec in ecorrs: @@ -363,7 +363,7 @@ def get_noise_weights(self, toas, nweights=None): """ ecorrs = self.get_ecorrs() if nweights is None: - ts = (toas.table["tdbld"].quantity * u.day).to(u.s).value + ts = (toas.table["tdbld"].quantity * u.day).to_value(u.s) nweights = [ get_ecorr_nweights(ts[ec.select_toa_mask(toas)]) for ec in ecorrs ] @@ -371,7 +371,7 @@ def get_noise_weights(self, toas, nweights=None): weights = np.zeros(nc) nctot = 0 for ec, nn in zip(ecorrs, nweights): - weights[nctot : nn + nctot] = ec.quantity.to(u.s).value ** 2 + weights[nctot : nn + nctot] = ec.quantity.to_value(u.s) ** 2 nctot += nn return weights @@ -462,7 +462,7 @@ def get_noise_basis(self, toas): See the documentation for pl_dm_basis_weight_pair function for details.""" tbl = toas.table - t = (tbl["tdbld"].quantity * u.day).to(u.s).value + t = (tbl["tdbld"].quantity * u.day).to_value(u.s) freqs = self._parent.barycentric_radio_freq(toas).to(u.MHz) fref = 1400 * u.MHz D = (fref.value / freqs.value) ** 2 @@ -476,7 +476,7 @@ def get_noise_weights(self, toas): See the documentation for pl_dm_basis_weight_pair for details.""" tbl = toas.table - t = (tbl["tdbld"].quantity * u.day).to(u.s).value + t = (tbl["tdbld"].quantity * u.day).to_value(u.s) amp, gam, nf = self.get_pl_vals() Ffreqs = get_rednoise_freqs(t, nf) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] @@ -592,7 +592,7 @@ def get_noise_basis(self, toas): See the documentation for pl_rn_basis_weight_pair function for details.""" tbl = toas.table - t = (tbl["tdbld"].quantity * u.day).to(u.s).value + t = (tbl["tdbld"].quantity * u.day).to_value(u.s) nf = self.get_pl_vals()[2] return create_fourier_design_matrix(t, nf) @@ -602,7 +602,7 @@ def get_noise_weights(self, toas): See the documentation for pl_rn_basis_weight_pair for details.""" tbl = toas.table - t = (tbl["tdbld"].quantity * u.day).to(u.s).value + t = (tbl["tdbld"].quantity * u.day).to_value(u.s) amp, gam, nf = self.get_pl_vals() Ffreqs = get_rednoise_freqs(t, nf) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index 5bf0e9d48..8af5aa1d7 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -409,7 +409,7 @@ def _print_uncertainty(self, uncertainty): This converts the :class:`~astropy.units.Quantity` provided to the appropriate units, extracts the value, and converts that to a string. """ - return str(uncertainty.to(self.units).value) + return str(uncertainty.to_value(self.units)) def __repr__(self): out = "{0:16s}{1:20s}".format(f"{self.__class__.__name__}(", self.name) @@ -783,7 +783,7 @@ def _set_uncertainty(self, val): def str_quantity(self, quan): """Quantity as a string (for floating-point values).""" - v = quan.to(self.units).value + v = quan.to_value(self.units) if self._long_double and not isinstance(v, np.longdouble): raise ValueError( "Parameter is supposed to contain long double values but contains a float" @@ -798,9 +798,9 @@ def _get_value(self, quan): return quan elif isinstance(quan, list): # for pairParamters - return [x.to(self.units).value for x in quan] + return [x.to_value(self.units) for x in quan] else: - return quan.to(self.units).value + return quan.to_value(self.units) def as_ufloat(self, units=None): """Return the parameter as a :class:`uncertainties.ufloat` @@ -1606,7 +1606,7 @@ def str_quantity(self, quan): return self.param_comp.str_quantity(quan) def _print_uncertainty(self, uncertainty): - return str(uncertainty.to(self.units).value) + return str(uncertainty.to_value(self.units)) def name_matches(self, name): return self.param_comp.name_matches(name) @@ -2239,8 +2239,8 @@ def str_quantity(self, quan): if len(quan) != 2: raise ValueError(f"Don't know how to print this as a pair: {quan}") - v0 = quan[0].to(self.units).value - v1 = quan[1].to(self.units).value + v0 = quan[0].to_value(self.units) + v1 = quan[1].to_value(self.units) if self._long_double: if not isinstance(v0, np.longdouble): raise TypeError( diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index b27dc923c..725abcf21 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -508,7 +508,7 @@ def change_binary_epoch(self, new_epoch): t0_ld = self.T0.quantity.tdb.mjd_long dt = (new_epoch.tdb.mjd_long - t0_ld) * u.day d_orbits = dt / PB - PBDOT * dt**2 / (2.0 * PB**2) - n_orbits = np.round(d_orbits.to(u.Unit(""))) + n_orbits = np.round(d_orbits.to(u.dimensionless_unscaled)) if n_orbits == 0: return dt_integer_orbits = PB * n_orbits + PB * PBDOT * n_orbits**2 / 2.0 diff --git a/src/pint/models/stand_alone_psr_binaries/DDK_model.py b/src/pint/models/stand_alone_psr_binaries/DDK_model.py index 709033aae..4c94f8119 100644 --- a/src/pint/models/stand_alone_psr_binaries/DDK_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DDK_model.py @@ -175,31 +175,35 @@ def kin(self): def d_SINI_d_KIN(self): # with u.set_enabled_equivalencies(u.dimensionless_angles()): - return np.cos(self.kin()).to(u.Unit("") / self.KIN.unit) + return np.cos(self.kin()).to(u.dimensionless_unscaled / self.KIN.unit) def d_SINI_d_KOM(self): if not self.K96: - return np.cos(self.kin()) * u.Unit("") / self.KOM.unit + return np.cos(self.kin()) * u.dimensionless_unscaled / self.KOM.unit d_si_d_kom = ( (-self.PMLONG_DDK * self.cos_KOM - self.PMLAT_DDK * self.sin_KOM) * self.tt0 * np.cos(self.kin()) ) # with u.set_enabled_equivalencies(u.dimensionless_angles()): - return d_si_d_kom.to(u.Unit("") / self.KOM.unit) + return d_si_d_kom.to(u.dimensionless_unscaled / self.KOM.unit) def d_SINI_d_T0(self): if not self.K96: - return np.ones(len(self.tt0)) * u.Unit("") / self.T0.unit + return np.ones(len(self.tt0)) * u.dimensionless_unscaled / self.T0.unit d_si_d_kom = -(-self.PMLONG_DDK * self.sin_KOM + self.PMLAT_DDK * self.cos_KOM) - return d_si_d_kom.to(u.Unit("") / self.T0.unit) + return d_si_d_kom.to(u.dimensionless_unscaled / self.T0.unit) def d_SINI_d_par(self, par): par_obj = getattr(self, par) try: ko_func = getattr(self, f"d_SINI_d_{par}") except Exception: - ko_func = lambda: np.zeros(len(self.tt0)) * u.Unit("") / par_obj.unit + ko_func = ( + lambda: np.zeros(len(self.tt0)) + * u.dimensionless_unscaled + / par_obj.unit + ) return ko_func() def d_kin_proper_motion_d_KOM(self): @@ -394,7 +398,7 @@ def delta_sini_parallax(self): / PX_kpc * (self.delta_I0() * self.sin_KOM - self.delta_J0() * self.cos_KOM) ) - return delta_sini.to("") + return delta_sini.to(u.dimensionless_unscaled) def delta_a1_parallax(self): """ diff --git a/src/pint/models/stand_alone_psr_binaries/DD_model.py b/src/pint/models/stand_alone_psr_binaries/DD_model.py index e9b82bdc3..e15b7b085 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -149,7 +149,7 @@ def d_omega_d_OMDOT(self): n = 2*pi/PB dOmega/dOMDOT = PB/2*pi*nu """ - PB = (self.pb()).to("second") + PB = (self.pb()).to(u.s) nu = self.nu() return PB / (2 * np.pi * u.rad) * nu @@ -568,7 +568,7 @@ def nhat(self): n = 2*pi/PB # should here be M() """ cosE = np.cos(self.E()) - return 2.0 * np.pi / self.pb().to("second") / (1 - self.ecc() * cosE) + return 2.0 * np.pi / self.pb().to(u.s) / (1 - self.ecc() * cosE) def d_nhat_d_par(self, par): """Derivative. diff --git a/src/pint/models/stand_alone_psr_binaries/ELL1H_model.py b/src/pint/models/stand_alone_psr_binaries/ELL1H_model.py index a976baef1..c2478fbc0 100644 --- a/src/pint/models/stand_alone_psr_binaries/ELL1H_model.py +++ b/src/pint/models/stand_alone_psr_binaries/ELL1H_model.py @@ -46,8 +46,8 @@ def __init__(self): { "H3": 0.0 * u.second, "H4": 0.0 * u.second, - "STIGMA": 0.0 * u.Unit(""), - "NHARMS": 3 * u.Unit(""), + "STIGMA": 0.0 * u.dimensionless_unscaled, + "NHARMS": 3 * u.dimensionless_unscaled, } ) self.binary_params = list(self.param_default_value.keys()) diff --git a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py index b3fd8133d..1191a424b 100644 --- a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py +++ b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py @@ -45,7 +45,7 @@ def ttasc(self): t = self.t if not hasattr(self.t, "unit") or self.t.unit is None: t = self.t * u.day - return (t - self.TASC.value * u.day).to("second") + return (t - self.TASC.value * u.day).to(u.s) def a1(self): """ELL1 model a1 calculation. @@ -103,14 +103,14 @@ def Phi(self): return self.M() def orbits_ELL1(self): - PB = (self.pb()).to("second") + PB = (self.pb()).to(u.s) PBDOT = self.pbdot() ttasc = self.ttasc() return (ttasc / PB - 0.5 * PBDOT * (ttasc / PB) ** 2).decompose() def d_Phi_d_TASC(self): """dPhi/dTASC""" - PB = self.pb().to("second") + PB = self.pb().to(u.s) PBDOT = self.pbdot() ttasc = self.ttasc() return (PBDOT * ttasc / PB - 1.0) * 2 * np.pi * u.rad / PB @@ -157,7 +157,7 @@ def delayI(self): Dre = self.delayR() Drep = self.Drep() Drepp = self.Drepp() - PB = self.pb().to("second") + PB = self.pb().to(u.s) nhat = 2 * np.pi / self.pb() return ( Dre @@ -182,7 +182,7 @@ def d_delayI_d_par(self, par): Dre = self.delayR() Drep = self.Drep() Drepp = self.Drepp() - PB = self.pb().to("second") + PB = self.pb().to(u.s) nhat = 2 * np.pi / self.pb() d_delayI_d_Dre = ( @@ -214,7 +214,7 @@ def ELL1_ecc(self): def ELL1_T0(self): return self.TASC + self.pb() / (2 * np.pi) * ( np.arctan(self.eps1() / self.eps2()) - ).to(u.Unit(""), equivalencies=u.dimensionless_angles()) + ).to(u.dimensionless_unscaled, equivalencies=u.dimensionless_angles()) ############################### def d_delayR_da1(self): diff --git a/src/pint/models/stand_alone_psr_binaries/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index c62f174a3..5cad0a48f 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -376,7 +376,7 @@ def get_tt0(self, barycentricTOA): T0 = self.T0 if not hasattr(barycentricTOA, "unit") or barycentricTOA.unit is None: barycentricTOA = barycentricTOA * u.day - return (barycentricTOA - T0).to("second") + return (barycentricTOA - T0).to(u.s) def ecc(self): """Calculate eccentricity with EDOT""" @@ -473,7 +473,7 @@ def d_M_d_par(self, par): par_obj = getattr(self, par) result = self.orbits_cls.d_orbits_d_par(par) with u.set_enabled_equivalencies(u.dimensionless_angles()): - result = result.to(u.Unit("") / par_obj.unit) + result = result.to(u.dimensionless_unscaled / par_obj.unit) return result ############################################### @@ -629,7 +629,7 @@ def d_nu_d_par(self, par): return np.zeros(len(self.tt0)) * nu.unit / par_obj.unit def omega(self): - PB = self.pb().to("second") + PB = self.pb().to(u.s) OMDOT = self.OMDOT OM = self.OM return OM + OMDOT * self.tt0 @@ -682,7 +682,7 @@ def pbprime(self): return self.pb() - self.pbdot() * self.tt0 def P(self): - return self.P0 + self.P1 * (self.t - self.PEPOCH).to("second") + return self.P0 + self.P1 * (self.t - self.PEPOCH).to(u.s) def t0(self): return self.t - self.PEPOCH diff --git a/src/pint/models/stand_alone_psr_binaries/binary_orbits.py b/src/pint/models/stand_alone_psr_binaries/binary_orbits.py index 3a1d60bbb..68704321c 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_orbits.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_orbits.py @@ -97,7 +97,7 @@ def __init__(self, parent, orbit_params=["PB", "PBDOT", "XPBDOT", "T0"]): def orbits(self): """Orbital phase (number of orbits since T0).""" - PB = self.PB.to("second") + PB = self.PB.to(u.s) PBDOT = self.PBDOT XPBDOT = self.XPBDOT return ( @@ -114,14 +114,14 @@ def pbdot_orbit(self): def d_orbits_d_T0(self): """The derivatve of orbits with respect to T0.""" - PB = self.PB.to("second") + PB = self.PB.to(u.s) PBDOT = self.PBDOT XPBDOT = self.XPBDOT return ((PBDOT - XPBDOT) * self.tt0 / PB - 1.0) * 2 * np.pi * u.rad / PB def d_orbits_d_PB(self): """dM/dPB this could be a generic function""" - PB = self.PB.to("second") + PB = self.PB.to(u.s) PBDOT = self.PBDOT XPBDOT = self.XPBDOT return ( @@ -133,12 +133,12 @@ def d_orbits_d_PB(self): def d_orbits_d_PBDOT(self): """dM/dPBDOT this could be a generic function""" - PB = self.PB.to("second") + PB = self.PB.to(u.s) return -np.pi * u.rad * self.tt0**2 / PB**2 def d_orbits_d_XPBDOT(self): """dM/dPBDOT this could be a generic function""" - PB = self.PB.to("second") + PB = self.PB.to(u.s) return -np.pi * u.rad * self.tt0**2 / PB**2 def d_pbprime_d_PB(self): diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index b6abe6270..e1d1ba0d3 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -1470,7 +1470,7 @@ def toa_covariance_matrix(self, toas): """ result = np.zeros((len(toas), len(toas))) if "ScaleToaError" not in self.components: - result += np.diag(toas.table["error"].quantity.to(u.s).value ** 2) + result += np.diag(toas.table["error"].quantity.to_value(u.s) ** 2) for nf in self.covariance_matrix_funcs: result += nf(toas) diff --git a/src/pint/models/troposphere_delay.py b/src/pint/models/troposphere_delay.py index 911f2ad12..9b482ef26 100644 --- a/src/pint/models/troposphere_delay.py +++ b/src/pint/models/troposphere_delay.py @@ -236,7 +236,7 @@ def pressure_from_altitude(self, H): gph = self.EARTH_R * H / (self.EARTH_R + H) # geopotential height if gph > 11 * u.km: log.warning("Pressure approximation invalid for elevations above 11 km") - T = 288.15 - 0.0065 * H.to(u.m).value # temperature lapse + T = 288.15 - 0.0065 * H.to_value(u.m) # temperature lapse return 101.325 * (288.15 / T) ** -5.25575 * u.kPa def zenith_delay(self, lat, H): @@ -314,7 +314,7 @@ def mapping_function(self, alt, lat, H, mjd): # now add in the mapping correction based on height fcorrection = self._herring_map(alt, self.A_HT, self.B_HT, self.C_HT) - return baseMap + (1 / np.sin(alt) - fcorrection) * H.to(u.km).value + return baseMap + (1 / np.sin(alt) - fcorrection) * H.to_value(u.km) def wet_map(self, alt, lat): """This is very similar to the normal mapping function except it uses different From 7737661d4d3e62af1ec18d33b96a1f6ef7aeb1fc Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 25 Sep 2023 10:53:33 -0500 Subject: [PATCH 18/19] changelog --- CHANGELOG-unreleased.md | 1 + src/pint/observatory/clock_file.py | 2 +- src/pint/observatory/satellite_obs.py | 2 +- src/pint/observatory/topo_obs.py | 2 +- src/pint/output/publish.py | 6 ++++-- 5 files changed, 8 insertions(+), 5 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 56f9dc5a5..ace043a77 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -30,4 +30,5 @@ the released changes. - Fixed an incorrect docstring in `pbprime()` functions. - Fix ICRS -> ECL conversion when parameter uncertainties are not set. - `get_TOAs` raises an exception upon finding mixed narrowband and wideband TOAs in a tim file. `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag. +- Optimized `astropy.units` operations ### Removed diff --git a/src/pint/observatory/clock_file.py b/src/pint/observatory/clock_file.py index 2d3fe0667..73b184ade 100644 --- a/src/pint/observatory/clock_file.py +++ b/src/pint/observatory/clock_file.py @@ -178,7 +178,7 @@ def evaluate(self, t, limits="warn"): raise ClockCorrectionOutOfRange(msg) # Can't pass Times directly to np.interp. This should be OK: - return np.interp(t.mjd, self.time.mjd, self.clock.to(u.us).value) * u.us + return np.interp(t.mjd, self.time.mjd, self.clock.to_value(u.us)) * u.us def last_correction_mjd(self): """Last MJD for which corrections are available.""" diff --git a/src/pint/observatory/satellite_obs.py b/src/pint/observatory/satellite_obs.py index 2021c0b48..a38ef2858 100644 --- a/src/pint/observatory/satellite_obs.py +++ b/src/pint/observatory/satellite_obs.py @@ -70,7 +70,7 @@ def load_Fermi_FT2(ft2_filename): dt = mjds_TT[1] - mjds_TT[0] log.info(f"FT2 spacing is {str(dt.to(u.s))}") # Use "spacing" argument for gradient to handle nonuniform entries - tt = mjds_TT.to(u.s).value + tt = mjds_TT.to_value(u.s) Vx = np.gradient(X.value, tt) * u.m / u.s Vy = np.gradient(Y.value, tt) * u.m / u.s Vz = np.gradient(Z.value, tt) * u.m / u.s diff --git a/src/pint/observatory/topo_obs.py b/src/pint/observatory/topo_obs.py index b40e5ccc5..0af8507d7 100644 --- a/src/pint/observatory/topo_obs.py +++ b/src/pint/observatory/topo_obs.py @@ -391,7 +391,7 @@ def _get_TDB_ephem(self, t, ephem): topo_tdb_tt = geo_tdb_tt - topo_time_corr return Time( t.tt.jd1 - JD_MJD, - t.tt.jd2 - topo_tdb_tt.to(u.day).value, + t.tt.jd2 - topo_tdb_tt.to_value(u.day), format="pulsar_mjd", scale="tdb", location=self.earth_location_itrf(), diff --git a/src/pint/output/publish.py b/src/pint/output/publish.py index 2cb82b0ac..a092f04ff 100644 --- a/src/pint/output/publish.py +++ b/src/pint/output/publish.py @@ -17,6 +17,8 @@ ) from pint.toa import TOAs from pint.residuals import Residuals, WidebandTOAResiduals +from pint import dmu +from astropy import units as u from io import StringIO import numpy as np @@ -242,11 +244,11 @@ def publish( f"Fitting method \\dotfill & {fit_method} \\\\ \n" ) tex.write( - f"RMS TOA residuals ($\\mu s$) \\dotfill & {toares.calc_time_resids().to('us').value.std():.2f} \\\\ \n" + f"RMS TOA residuals ($\\mu s$) \\dotfill & {toares.calc_time_resids().to_value(u.microsecond).std():.2f} \\\\ \n" ) if toas.is_wideband(): tex.write( - f"RMS DM residuals (pc / cm3) \\dotfill & {dmres.calc_resids().to('pc/cm^3').value.std():.2f} \\\\ \n" + f"RMS DM residuals (pc / cm3) \\dotfill & {dmres.calc_resids().to_value(dmu).std():.2f} \\\\ \n" ) tex.write( f"$\\chi^2$ \\dotfill & {res.chi2:.2f} \\\\ \n" From 72354c340e589554bfc87321cc723c1cde04ad05 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 28 Sep 2023 16:02:48 -0500 Subject: [PATCH 19/19] = --- src/pint/derived_quantities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py index 104b73d8d..d8d9cdf3b 100644 --- a/src/pint/derived_quantities.py +++ b/src/pint/derived_quantities.py @@ -279,7 +279,7 @@ def pulsar_B_lightcyl(f: u.Hz, fdot: u.Hz / u.s): 2.9e8 * u.G * p.to_value(u.s) ** (-5.0 / 2.0) - * np.sqrt(pd.to(u.dimensionless_unscaled).value) + * np.sqrt(pd.to_value(u.dimensionless_unscaled)) )