From 98a7011833406d117aa7db55334d022cccdcc611 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 20 May 2026 11:02:55 -0500 Subject: [PATCH 01/11] testing --- src/pint/__init__.py | 2 +- src/pint/fitter.py | 45 ++++---- src/pint/models/astrometry.py | 102 +++++++++--------- src/pint/models/binary_bt.py | 2 +- src/pint/models/dispersion_model.py | 2 +- src/pint/models/jump.py | 16 +-- src/pint/models/parameter.py | 11 +- src/pint/models/spindown.py | 4 +- .../binary_generic.py | 4 +- src/pint/models/timing_model.py | 4 +- src/pint/pulsar_mjd.py | 86 ++++++++++----- src/pint/toa.py | 4 +- src/pint/utils.py | 55 +++++----- tests/test_binconvert.py | 28 +++-- tests/test_change_epoch.py | 3 +- tests/test_downhill_fitter.py | 2 - tests/test_precision.py | 11 +- 17 files changed, 219 insertions(+), 162 deletions(-) diff --git a/src/pint/__init__.py b/src/pint/__init__.py index 8ab94f9279..bb053dd591 100644 --- a/src/pint/__init__.py +++ b/src/pint/__init__.py @@ -68,7 +68,7 @@ # define equivalency for astropy units light_second_equivalency = [(ls, si.second, lambda x: x, lambda x: x)] # hourangle_second unit -hourangle_second = u.def_unit("hourangle_second", u.hourangle / np.longdouble(3600.0)) +hourangle_second = u.def_unit("hourangle_second", u.hourangle / np.float64(3600.0)) # Following are from here: # http://ssd.jpl.nasa.gov/?constants (grabbed on 30 Dec 2013) diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 3df684fa96..a85461cca5 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -1269,10 +1269,10 @@ def take_step(self, step, lambda_=1) -> "WLSState": def parameter_covariance_matrix(self) -> CovarianceMatrix: # make sure we compute the SVD self.step - # Sigma = np.dot(Vt.T / s, U.T) + # Sigma = np.matmul(Vt.T / s, U.T) # The post-fit parameter covariance matrix # Sigma = V s^-2 V^T - Sigma = np.dot(self.Vt.T / (self.s**2), self.Vt) + Sigma = np.matmul(self.Vt.T / (self.s**2), self.Vt) return CovarianceMatrix( (Sigma / self.fac).T / self.fac, self.parameter_covariance_matrix_labels ) @@ -1589,8 +1589,8 @@ def mtcm_mtcy_mtcmplain(self): ).matrix cf = scipy.linalg.cho_factor(cov) cm = scipy.linalg.cho_solve(cf, self.M) - mtcm = np.dot(self.M.T, cm) - mtcy = np.dot(cm.T, residuals) + mtcm = np.matmul(self.M.T, cm) + mtcy = np.matmul(cm.T, residuals) mtcmplain = mtcm else: Nvec = ( @@ -1607,10 +1607,10 @@ def mtcm_mtcy_mtcmplain(self): ** 2 ) cinv = 1 / Nvec - mtcm = np.dot(self.M.T, cinv[:, None] * self.M) + mtcm = np.matmul(self.M.T, cinv[:, None] * self.M) mtcmplain = mtcm mtcm += np.diag(self.phiinv) - mtcy = np.dot(self.M.T, cinv * residuals) + mtcy = np.matmul(self.M.T, cinv * residuals) return mtcm, mtcy, mtcmplain @cached_property @@ -1629,7 +1629,7 @@ def mtcmplain(self): def U_s_Vt_xhat(self): U, s, Vt = scipy.linalg.svd(self.mtcm, full_matrices=False) s = apply_Sdiag_threshold(s, Vt, self.threshold, self.params) - xhat = np.dot(Vt.T, np.dot(U.T, self.mtcy) / s) + xhat = np.matmul(Vt.T, np.matmul(U.T, self.mtcy) / s) return U, s, Vt, xhat @cached_property @@ -1661,7 +1661,7 @@ def take_step(self, step, lambda_=1) -> "WidebandState": @cached_property def parameter_covariance_matrix(self) -> CovarianceMatrix: # make sure we compute the SVD - xvar = np.dot(self.Vt.T / self.s, self.Vt) + xvar = np.matmul(self.Vt.T / self.s, self.Vt) # is this the best place to do this? covariance_matrix_labels = { param: (i, i + 1, unit) @@ -2337,16 +2337,16 @@ def fit_toas( cov = self.get_noise_covariancematrix().matrix cf = scipy.linalg.cho_factor(cov) cm = scipy.linalg.cho_solve(cf, M) - mtcm = np.dot(M.T, cm) - mtcy = np.dot(cm.T, residuals) + mtcm = np.matmul(M.T, cm) + mtcy = np.matmul(cm.T, residuals) # mtcm, mtcy = get_gls_mtcm_mtcy_fullcov(cov, M, residuals) else: phiinv /= norm**2 Nvec = self.scaled_all_sigma() ** 2 cinv = 1 / Nvec - mtcm = np.dot(M.T, cinv[:, None] * M) + mtcm = np.matmul(M.T, cinv[:, None] * M) mtcm += np.diag(phiinv) - mtcy = np.dot(M.T, cinv * residuals) + mtcy = np.matmul(M.T, cinv * residuals) # mtcm, mtcy = get_gls_mtcm_mtcy(phiinv, Nvec, M, residuals) if threshold <= 0: @@ -2357,12 +2357,12 @@ def fit_toas( else: xvar, xhat = _solve_svd(mtcm, mtcy, threshold, params) - newres = residuals - np.dot(M, xhat) + newres = residuals - np.matmul(M, xhat) # compute linearized chisq if full_cov: - chi2 = np.dot(newres, scipy.linalg.cho_solve(cf, newres)) + chi2 = np.matmul(newres, scipy.linalg.cho_solve(cf, newres)) else: - chi2 = np.dot(newres, cinv * newres) + np.dot(xhat, phiinv * xhat) + chi2 = np.matmul(newres, cinv * newres) + np.matmul(xhat, phiinv * xhat) # compute absolute estimates, normalized errors, covariance matrix dpars = xhat / norm @@ -2460,7 +2460,7 @@ def fit_toas( s = apply_Sdiag_threshold(s, Vt, threshold, current_state.params) - dx = np.dot(Vt.T, np.dot(U.T, b) / s) + dx = np.matmul(Vt.T, np.matmul(U.T, b) / s) step = dx / current_state.norm @@ -2676,6 +2676,7 @@ def fit_wls_svd( # M2 = U S V^T # Both U and V^T are orthogonal matrices. + print(np.any(np.isnan(M2))) U, Sdiag, VT = scipy.linalg.svd(M2, full_matrices=False) # Deal with degeneracies by replacing very small singular @@ -2704,8 +2705,8 @@ def get_gls_mtcm_mtcy_fullcov( """ cf = scipy.linalg.cho_factor(cov) cm = scipy.linalg.cho_solve(cf, M) - mtcm = np.dot(M.T, cm) - mtcy = np.dot(cm.T, residuals) + mtcm = np.matmul(M.T, cm) + mtcy = np.matmul(cm.T, residuals) return mtcm, mtcy @@ -2720,9 +2721,9 @@ def get_gls_mtcm_mtcy( correlated noise basis, and residuals y (`residuals`). """ cinv = 1 / Nvec - mtcm = np.dot(M.T, cinv[:, None] * M) + mtcm = np.matmul(M.T, cinv[:, None] * M) mtcm += np.diag(phiinv) - mtcy = np.dot(M.T, cinv * residuals) + mtcy = np.matmul(M.T, cinv * residuals) return mtcm, mtcy @@ -2751,8 +2752,8 @@ def _solve_svd( corresponding parameters.""" U, s, Vt = scipy.linalg.svd(mtcm, full_matrices=False) s = apply_Sdiag_threshold(s, Vt, threshold, params) - xvar = np.dot(Vt.T / s, Vt) - xhat = np.dot(Vt.T, np.dot(U.T, mtcy) / s) + xvar = np.matmul(Vt.T / s, Vt) + xhat = np.matmul(Vt.T, np.matmul(U.T, mtcy) / s) return xvar, xhat diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index c08a6fd494..b0f8b03aca 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -382,10 +382,10 @@ def get_psr_coords(self, epoch: Optional[time_like] = None) -> coords.SkyCoord: """ if epoch is None or (self.PMRA.value == 0.0 and self.PMDEC.value == 0.0): return coords.SkyCoord( - ra=self.RAJ.quantity, - dec=self.DECJ.quantity, - pm_ra_cosdec=self.PMRA.quantity, - pm_dec=self.PMDEC.quantity, + ra=self.RAJ.quantity.astype(np.float64), + dec=self.DECJ.quantity.astype(np.float64), + pm_ra_cosdec=self.PMRA.quantity.astype(np.float64), + pm_dec=self.PMDEC.quantity.astype(np.float64), obstime=self.POSEPOCH.quantity, frame=coords.ICRS, ) @@ -511,17 +511,19 @@ def ssb_to_psb_xyz_ICRS(self, epoch: Optional[time_like] = None) -> u.Quantity: warnings.simplefilter("ignore", ErfaWarning) # note that starpm wants mu_alpha not mu_alpha * cos(delta) starpmout = pmsafe( - self.RAJ.quantity.to_value(u.radian), - self.DECJ.quantity.to_value(u.radian), - self.PMRA.quantity.to_value(u.radian / u.yr) - / np.cos(self.DECJ.quantity).value, - self.PMDEC.quantity.to_value(u.radian / u.yr), - self.PX.quantity.to_value(u.arcsec), + np.float64(self.RAJ.quantity.to_value(u.radian)), + np.float64(self.DECJ.quantity.to_value(u.radian)), + np.float64( + self.PMRA.quantity.to_value(u.radian / u.yr) + / np.cos(self.DECJ.quantity).value + ), + np.float64(self.PMDEC.quantity.to_value(u.radian / u.yr)), + np.float64(self.PX.quantity.to_value(u.arcsec)), 0.0, - self.POSEPOCH.quantity.jd1, - self.POSEPOCH.quantity.jd2, + np.float64(self.POSEPOCH.quantity.jd1), + np.float64(self.POSEPOCH.quantity.jd2), jd1, - jd2, + np.float64(jd2), ) # ra,dec now in radians ra, dec = starpmout[0], starpmout[1] @@ -702,19 +704,19 @@ def as_ECL( # put it in here as pm_ra_cosdec since astropy complains otherwise dt = 1 * u.yr c = coords.SkyCoord( - ra=self.RAJ.quantity, - dec=self.DECJ.quantity, + ra=self.RAJ.quantity.astype(np.float64), + dec=self.DECJ.quantity.astype(np.float64), obstime=self.POSEPOCH.quantity, pm_ra_cosdec=( self.RAJ.uncertainty * np.cos(self.DECJ.quantity) / dt if self.RAJ.uncertainty is not None else 0 * self.RAJ.units / dt - ), + ).astype(np.float64), pm_dec=( self.DECJ.uncertainty / dt if self.DECJ.uncertainty is not None else 0 * self.DECJ.units / dt - ), + ).astype(np.float64), frame=coords.ICRS, ) c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl)) @@ -723,19 +725,19 @@ def as_ECL( # use fake proper motions to convert uncertainties on proper motion # assume that the PM_RA _does_ include cos(DEC) c = coords.SkyCoord( - ra=self.RAJ.quantity, - dec=self.DECJ.quantity, + ra=self.RAJ.quantity.astype(np.float64), + dec=self.DECJ.quantity.astype(np.float64), obstime=self.POSEPOCH.quantity, pm_ra_cosdec=( self.PMRA.uncertainty if self.PMRA.uncertainty is not None else 0 * self.PMRA.units - ), + ).astype(np.float64), pm_dec=( self.PMDEC.uncertainty if self.PMDEC.uncertainty is not None else 0 * self.PMDEC.units - ), + ).astype(np.float64), frame=coords.ICRS, ) c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl)) @@ -875,10 +877,10 @@ def get_psr_coords(self, epoch: Optional[time_like] = None) -> coords.SkyCoord: # Compute only once return coords.SkyCoord( obliquity=obliquity, - lon=self.ELONG.quantity, - lat=self.ELAT.quantity, - pm_lon_coslat=self.PMELONG.quantity, - pm_lat=self.PMELAT.quantity, + lon=self.ELONG.quantity.astype(np.float64), + lat=self.ELAT.quantity.astype(np.float64), + pm_lon_coslat=self.PMELONG.quantity.astype(np.float64), + pm_lat=self.PMELAT.quantity.astype(np.float64), obstime=self.POSEPOCH.quantity, frame=PulsarEcliptic, ) @@ -1010,16 +1012,16 @@ def ssb_to_psb_xyz_ECL( warnings.simplefilter("ignore", ErfaWarning) # note that pmsafe wants mu_lon not mu_lon * cos(lat) starpmout = pmsafe( - lon, - lat, - pm_lon, - pm_lat, - self.PX.quantity.to_value(u.arcsec), + np.float64(lon), + np.float64(lat), + np.float64(pm_lon), + np.float64(pm_lat), + np.float64(self.PX.quantity.to_value(u.arcsec)), 0.0, - self.POSEPOCH.quantity.jd1, - self.POSEPOCH.quantity.jd2, - jd1, - jd2, + np.float64(self.POSEPOCH.quantity.jd1), + np.float64(self.POSEPOCH.quantity.jd2), + np.float64(jd1), + np.float64(jd2), ) # lon,lat now in radians lon, lat = starpmout[0], starpmout[1] @@ -1242,20 +1244,20 @@ def as_ECL( # put it in here as pm_ra_cosdec since astropy complains otherwise dt = 1 * u.yr c = coords.SkyCoord( - lon=self.ELONG.quantity, - lat=self.ELAT.quantity, + lon=self.ELONG.quantity.astype(np.float64), + lat=self.ELAT.quantity.astype(np.float64), obliquity=OBL[self.ECL.value], obstime=self.POSEPOCH.quantity, pm_lon_coslat=( self.ELONG.uncertainty * np.cos(self.ELAT.quantity) / dt if self.ELONG.uncertainty is not None else 0 * self.ELONG.units / dt - ), + ).astype(np.float64), pm_lat=( self.ELAT.uncertainty / dt if self.ELAT.uncertainty is not None else 0 * self.ELAT.units / dt - ), + ).astype(np.float64), frame=PulsarEcliptic, ) c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl)) @@ -1264,20 +1266,20 @@ def as_ECL( # use fake proper motions to convert uncertainties on proper motion # assume that PMELONG uncertainty includes cos(DEC) c = coords.SkyCoord( - lon=self.ELONG.quantity, - lat=self.ELAT.quantity, + lon=self.ELONG.quantity.astype(np.float64), + lat=self.ELAT.quantity.astype(np.float64), obliquity=OBL[self.ECL.value], obstime=self.POSEPOCH.quantity, pm_lon_coslat=( self.PMELONG.uncertainty if self.PMELONG.uncertainty is not None else 0 * self.PMELONG.units - ), + ).astype(np.float64), pm_lat=( self.PMELAT.uncertainty if self.PMELAT.uncertainty is not None else 0 * self.PMELAT.units - ), + ).astype(np.float64), frame=PulsarEcliptic, ) c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl)) @@ -1321,20 +1323,20 @@ def as_ICRS(self, epoch: Optional[time_like] = None) -> "AstrometryEquatorial": # put it in as pm_lon_coslat since astropy complains otherwise dt = 1 * u.yr c = coords.SkyCoord( - lon=self.ELONG.quantity, - lat=self.ELAT.quantity, + lon=self.ELONG.quantity.astype(np.float64), + lat=self.ELAT.quantity.astype(np.float64), obliquity=OBL[self.ECL.value], obstime=self.POSEPOCH.quantity, pm_lon_coslat=( self.ELONG.uncertainty * np.cos(self.ELAT.quantity) / dt if self.ELONG.uncertainty is not None else 0 * self.ELONG.units / dt - ), + ).astype(np.float64), pm_lat=( self.ELAT.uncertainty / dt if self.ELAT.uncertainty is not None else 0 * self.ELAT.units / dt - ), + ).astype(np.float64), frame=PulsarEcliptic, ) c_ICRS = c.transform_to(coords.ICRS) @@ -1344,20 +1346,20 @@ def as_ICRS(self, epoch: Optional[time_like] = None) -> "AstrometryEquatorial": # use fake proper motions to convert uncertainties on proper motion # assume that PMELONG uncertainty includes cos(DEC) c = coords.SkyCoord( - lon=self.ELONG.quantity, - lat=self.ELAT.quantity, + lon=self.ELONG.quantity.astype(np.float64), + lat=self.ELAT.quantity.astype(np.float64), obliquity=OBL[self.ECL.value], obstime=self.POSEPOCH.quantity, pm_lon_coslat=( self.PMELONG.uncertainty if self.PMELONG.uncertainty is not None else 0 * self.PMELONG.units - ), + ).astype(np.float64), pm_lat=( self.PMELAT.uncertainty if self.PMELAT.uncertainty is not None else 0 * self.PMELAT.units - ), + ).astype(np.float64), frame=PulsarEcliptic, ) c_ICRS = c.transform_to(coords.ICRS) diff --git a/src/pint/models/binary_bt.py b/src/pint/models/binary_bt.py index e055bd3993..c4a5abf64d 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -271,7 +271,7 @@ def add_piecewise_param(self, piece_index, **kwargs): param_unit = u.d if isinstance(paramx, u.quantity.Quantity): paramx = paramx.value - elif isinstance(paramx, np.float128): + elif isinstance(paramx, np.longdouble): paramx = paramx elif isinstance(paramx, Time): paramx = paramx.mjd diff --git a/src/pint/models/dispersion_model.py b/src/pint/models/dispersion_model.py index 39508d515a..e1b538fd6b 100644 --- a/src/pint/models/dispersion_model.py +++ b/src/pint/models/dispersion_model.py @@ -228,7 +228,7 @@ def base_dm(self, toas): dt_value = np.zeros(len(toas), dtype=np.longdouble) dm_terms_value = [d.value for d in dm_terms] dm = taylor_horner(dt_value, dm_terms_value) - return dm * self.DM.units + return dm.astype(np.float64) * self.DM.units def constant_dispersion_delay(self, toas, acc_delay=None): """This is a wrapper function for interacting with the TimingModel class""" diff --git a/src/pint/models/jump.py b/src/pint/models/jump.py index 8469aa4a7a..24c3ea4053 100644 --- a/src/pint/models/jump.py +++ b/src/pint/models/jump.py @@ -1,7 +1,7 @@ -"""Phase jumps. """ +"""Phase jumps.""" import astropy.units as u -import numpy +import numpy as np from loguru import logger as log from pint.models.parameter import maskParameter @@ -50,7 +50,7 @@ def jump_delay(self, toas, acc_delay=None): in the unit of seconds. """ tbl = toas.table - jdelay = numpy.zeros(len(tbl)) + jdelay = np.zeros(len(tbl)) for jump in self.jumps: jump_par = getattr(self, jump) mask = jump_par.select_toa_mask(toas) @@ -61,7 +61,7 @@ def jump_delay(self, toas, acc_delay=None): def d_delay_d_jump(self, toas, jump_param, acc_delay=None): tbl = toas.table - d_delay_d_j = numpy.zeros(len(tbl)) + d_delay_d_j = np.zeros(len(tbl)) jpar = getattr(self, jump_param) mask = jpar.select_toa_mask(toas) d_delay_d_j[mask] = -1.0 @@ -122,7 +122,7 @@ def jump_phase(self, toas, delay): """ tbl = toas.table # base this on the first available jump (doesn't have to be JUMP1) - jphase = numpy.zeros(len(tbl)) * ( + jphase = np.zeros(len(tbl)) * ( getattr(self, self.get_params_of_type("maskParameter")[0]).units * self._parent.F0.units ) @@ -131,13 +131,15 @@ def jump_phase(self, toas, delay): mask = jump_par.select_toa_mask(toas) # NOTE: Currently parfile jump value has opposite sign with our # phase calculation. - jphase[mask] += jump_par.quantity * self._parent.F0.quantity + jphase[mask] += jump_par.quantity * self._parent.F0.quantity.astype( + np.float64 + ) return jphase def d_phase_d_jump(self, toas, jump_param, delay): tbl = toas.table jpar = getattr(self, jump_param) - d_phase_d_j = numpy.zeros(len(tbl)) + d_phase_d_j = np.zeros(len(tbl)) mask = jpar.select_toa_mask(toas) d_phase_d_j[mask] = self._parent.F0.value return (d_phase_d_j * self._parent.F0.units).to(1 / u.second) diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index d560fae436..18a78b48b6 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -1186,7 +1186,9 @@ def _set_quantity(self, val): mjd float mjd string (in pulsar_mjd format) """ - if isinstance(val, numbers.Number): + if isinstance(val, numbers.Number) or ( + isinstance(val, np.ndarray) and val.ndim == 0 and not hasattr(val, "unit") + ): val = np.longdouble(val) result = time_from_longdouble(val, self.time_scale) elif isinstance(val, (str, bytes)): @@ -1349,7 +1351,12 @@ def _set_quantity(self, val): 2. float 3. number string """ - if isinstance(val, numbers.Number): + + # do not know why this is needed, but for some values + # that are QuadPrecision they are stored as ndim=0 arrays + if isinstance(val, numbers.Number) or ( + isinstance(val, np.ndarray) and val.ndim == 0 and not hasattr(val, "unit") + ): result = Angle(np.longdouble(val) * self.units) elif isinstance(val, str): # FIXME: what if the user included a unit suffix? diff --git a/src/pint/models/spindown.py b/src/pint/models/spindown.py index 76408d228c..79538efa95 100644 --- a/src/pint/models/spindown.py +++ b/src/pint/models/spindown.py @@ -170,9 +170,9 @@ def change_pepoch(self, new_epoch, toas=None, delay=None): first pulsar frame toa. """ if isinstance(new_epoch, Time): - new_epoch = Time(new_epoch, scale="tdb", precision=9) + new_epoch = Time(new_epoch, scale="tdb") else: - new_epoch = Time(new_epoch, scale="tdb", format="mjd", precision=9) + new_epoch = Time(new_epoch, scale="tdb", format="mjd") if self.PEPOCH.value is None: if toas is None or delay is None: 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 571a98182b..b7ff2ba2f1 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -350,14 +350,14 @@ def compute_eccentric_anomaly(self, eccentricity, mean_anomaly): """ if hasattr(eccentricity, "unit"): # FIXME: isn't this an error? - e = np.longdouble(eccentricity).value + e = np.longdouble(eccentricity) else: e = eccentricity if any(e < 0) or any(e >= 1): raise ValueError("Eccentricity should be in the range of [0,1).") if hasattr(mean_anomaly, "unit"): - ma = np.longdouble(mean_anomaly).value + ma = np.longdouble(mean_anomaly) else: ma = mean_anomaly k = lambda E: E - e * np.sin(E) - ma # Kepler Equation diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index a5bdd5cd61..ac48f95b0c 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -1600,7 +1600,7 @@ def delay( # Do NOT cycle through delay_funcs - cycle through components until cutoff for dc in self.DelayComponent_list[:idx]: for df in dc.delay_funcs_component: - delay += df(toas, delay) + delay += df(toas, delay).astype(np.float64) return delay def phase(self, toas: TOAs, abs_phase: Optional[bool] = None) -> Phase: @@ -2152,7 +2152,7 @@ def d_delay_d_param( ) -> u.Quantity: """Return the derivative of delay with respect to the parameter.""" par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) + result = np.longdouble(np.zeros(toas.ntoas)) << (u.s / par.units) delay_derivs = self.delay_deriv_funcs if param not in list(delay_derivs.keys()): raise AttributeError( diff --git a/src/pint/pulsar_mjd.py b/src/pint/pulsar_mjd.py index 6b3f274d7f..6747e90689 100644 --- a/src/pint/pulsar_mjd.py +++ b/src/pint/pulsar_mjd.py @@ -35,6 +35,7 @@ import numpy as np from astropy.time import Time from astropy.time.formats import TimeFormat +import pint try: maketrans = str.maketrans @@ -43,22 +44,6 @@ from string import maketrans -# This check is implemented in pint.utils, but we want to avoid circular imports -if np.finfo(np.longdouble).eps > 2e-19: - import warnings - - def readable_warning(message, category, filename, lineno, line=None): - return "%s: %s\n" % (category.__name__, message) - - warnings.formatwarning = readable_warning - - msg = ( - "This platform does not support extended precision " - "floating-point, and PINT will run at reduced precision." - ) - warnings.warn(msg, RuntimeWarning) - - __all__ = [ "Time", "PulsarMJD", @@ -82,6 +67,45 @@ def readable_warning(message, category, filename, lineno, line=None): "jds_to_mjds_pulsar", ] +use_npq = False +eps = np.finfo(np.longdouble).eps +longdoubledtype = np.longdouble +# This check is implemented in pint.utils, but we want to avoid circular imports +if eps > 2e-19: + import warnings + + def readable_warning(message, category, filename, lineno, line=None): + return "%s: %s\n" % (category.__name__, message) + + warnings.formatwarning = readable_warning + ok = False + try: + import numpy_quaddtype as npq + + ok = npq.epsilon < 2e-19 + except ImportError: + pass + + if not ok: + msg = ( + "This platform does not support extended precision " + "floating-point, and PINT will run at reduced precision." + ) + warnings.warn(msg, RuntimeWarning) + else: + use_npq = True + +if use_npq: + from numpy_quaddtype import QuadPrecDType, QuadPrecision + + # needed to register the dtype for string representations + # see https://github.com/nanograv/PINT/issues/1944#issuecomment-3774471273 + np._core.sctypeDict[str(QuadPrecDType())] = QuadPrecDType() + + np.longdouble = QuadPrecision + longdoubledtype = QuadPrecDType + eps = npq.epsilon + class PulsarMJD(TimeFormat): """Change handling of days with leap seconds. @@ -118,7 +142,7 @@ class MJDLong(TimeFormat): name = "mjd_long" def _check_val_type(self, val1, val2): - if val1.dtype != np.longdouble: + if val1.dtype != longdoubledtype: raise ValueError( "mjd_long requires a long double number but got {!r} of type {}".format( val1, val1.dtype @@ -126,7 +150,7 @@ def _check_val_type(self, val1, val2): ) if val2 is None: val2 = 0 - elif val2.dtype != np.longdouble: + elif val2.dtype != longdoubledtype: raise ValueError( "mjd_long requires a long double number but got {!r} of type {}".format( val2, val2.dtype @@ -263,10 +287,11 @@ def time_from_mjd_string(s, scale="utc", format="pulsar_mjd"): def time_from_longdouble(t, scale="utc", format="pulsar_mjd"): - t = np.longdouble(t) - i = float(np.floor(t)) - f = float(t - i) - return astropy.time.Time(val=i, val2=f, format=format, scale=scale) + if isinstance(t, str): + t = np.longdouble(t) + i = np.floor(t) + f = t - i + return astropy.time.Time(val=float(i), val2=float(f), format=format, scale=scale) def time_to_mjd_string(t): @@ -280,7 +305,7 @@ def time_to_mjd_string(t): return t.mjd_string -longdouble_mjd_eps = (70000 * u.day * np.finfo(np.longdouble).eps).to(u.ns) +longdouble_mjd_eps = (70000 * u.day * eps).to(u.ns) def time_to_longdouble(t): @@ -327,7 +352,14 @@ def data2longdouble(data): np.longdouble """ - return str2longdouble(data) if type(data) is str else np.longdouble(data) + if type(data) is str: + return str2longdouble(data) + elif isinstance(data, (np.ndarray, np.number)): + return data.astype(longdoubledtype) + elif isinstance(data, list): + return data2longdouble(np.array(data)) + else: + return np.longdouble(data) def quantity2longdouble_withunit(data): @@ -459,15 +491,15 @@ def _str_to_mjds(s): num, expon = ss.split("e") expon = int(expon) if expon < 0: - imjd, fmjd = 0, np.longdouble(ss) + imjd, fmjd = 0, (ss).astype(np.longdouble) else: mjd_s = num.split(".") # If input was given as an integer, add floating "0" if len(mjd_s) == 1: mjd_s.append("0") imjd_s, fmjd_s = mjd_s - imjd = np.longdouble(int(imjd_s)) - fmjd = np.longdouble(f"0.{fmjd_s}") + imjd = (int(imjd_s)).astype(np.longdouble) + fmjd = (f"0.{fmjd_s}").astype(np.longdouble) if ss.startswith("-"): fmjd = -fmjd imjd *= 10**expon diff --git a/src/pint/toa.py b/src/pint/toa.py index 69af944849..eab1b895c8 100644 --- a/src/pint/toa.py +++ b/src/pint/toa.py @@ -1978,7 +1978,7 @@ def phase_columns_from_flags(self): self.table["delta_pulse_number"] += dphs # Then, add pulse_number as a table column if possible - pns = np.array(self.get_flag_value("pn", np.nan, float)[0]) + pns = np.array(self.get_flag_value("pn", np.nan, float)[0]).astype(np.float64) if np.all(np.isnan(pns)): raise ValueError("No pulse numbers found") self.table["pulse_number"] = pns @@ -2005,7 +2005,7 @@ def compute_pulse_numbers(self, model: "TimingModel") -> None: # paulr: I think pulse numbers should be computed with abs_phase=True! delta_pulse_numbers = Phase(self.table["delta_pulse_number"]) phases = model.phase(self, abs_phase=True) + delta_pulse_numbers - self.table["pulse_number"] = phases.int + self.table["pulse_number"] = phases.int.astype(np.float64) self.table["pulse_number"].unit = u.dimensionless_unscaled def remove_pulse_numbers(self): diff --git a/src/pint/utils.py b/src/pint/utils.py index 8b117245b1..7d1471a9cc 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -163,7 +163,16 @@ def check_longdouble_precision() -> bool: Returns True if long doubles have enough precision to use PINT for sub-microsecond timing on this machine. """ - return np.finfo(np.longdouble).eps < 2e-19 + return True + if np.finfo(np.longdouble).eps < 2e-19: + return True + else: + try: + import numpy_quaddtype as npq + + return npq.epsilon < 2e-19 + except ImportError: + return False def require_longdouble_precision() -> None: @@ -2209,43 +2218,31 @@ def add_dummy_distance( return c if isinstance(c.frame, coords.builtin_frames.icrs.ICRS): - return ( - coords.SkyCoord( - ra=c.ra, - dec=c.dec, - pm_ra_cosdec=c.pm_ra_cosdec, - pm_dec=c.pm_dec, - obstime=c.obstime, - distance=distance, - frame=coords.ICRS, - ) - if hasattr(c, "pm_ra_cosdec") - else coords.SkyCoord( - ra=c.ra, - dec=c.dec, - pm_ra_cosdec=c.pm_ra, - pm_dec=c.pm_dec, - obstime=c.obstime, - distance=distance, - frame=coords.ICRS, - ) + return coords.SkyCoord( + ra=(c.ra).astype(np.float64), + dec=(c.dec).astype(np.float64), + pm_ra_cosdec=(c.pm_ra_cosdec).astype(np.float64), + pm_dec=(c.pm_dec).astype(np.float64), + obstime=c.obstime, + distance=distance, + frame=coords.ICRS, ) elif isinstance(c.frame, coords.builtin_frames.galactic.Galactic): return coords.SkyCoord( - l=c.l, - b=c.b, - pm_l_cosb=c.pm_l_cosb, - pm_b=c.pm_b, + l=(c.l).astype(np.float64), + b=(c.b).astype(np.float64), + pm_l_cosb=c.pm_l_cosb.astype(np.float64), + pm_b=c.pm_b.astype(np.float64), obstime=c.obstime, distance=distance, frame=coords.Galactic, ) elif isinstance(c.frame, pint.pulsar_ecliptic.PulsarEcliptic): return coords.SkyCoord( - lon=c.lon, - lat=c.lat, - pm_lon_coslat=c.pm_lon_coslat, - pm_lat=c.pm_lat, + lon=(c.lon).astype(np.float64), + lat=(c.lat).astype(np.float64), + pm_lon_coslat=c.pm_lon_coslat.astype(np.float64), + pm_lat=c.pm_lat.astype(np.float64), obstime=c.obstime, distance=distance, obliquity=c.obliquity, diff --git a/tests/test_binconvert.py b/tests/test_binconvert.py index 25c5ba87f5..4c4b6a183e 100644 --- a/tests/test_binconvert.py +++ b/tests/test_binconvert.py @@ -193,9 +193,17 @@ def test_ELL1_roundtrip(output): rtol=0.2, ), f"{p} uncertainty: {getattr(m, p).uncertainty_value} does not match {getattr(mback, p).uncertainty_value}" else: - assert ( - getattr(m, p).value == getattr(mback, p).value - ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" + # this used to be an equivalence test, but it seems like with QuadPrecision + # we need to allow for a little slop + if isinstance(getattr(m, p).value, np.longdouble): + assert np.isclose( + getattr(m, p).value, getattr(mback, p).value, atol=1e-15 + ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" + print(p, getattr(m, p).value - getattr(mback, p).value) + else: + assert ( + getattr(m, p).value == getattr(mback, p).value + ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" @pytest.mark.parametrize("output", ["ELL1", "ELL1H", "ELL1k", "DD", "BT", "DDS", "DDK"]) @@ -229,9 +237,17 @@ def test_ELL1_roundtripFB0(output): rtol=0.2, ), f"{p} uncertainty: {getattr(m, p).uncertainty_value} does not match {getattr(mback, p).uncertainty_value}" else: - assert ( - getattr(m, p).value == getattr(mback, p).value - ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" + # this used to be an equivalence test, but it seems like with QuadPrecision + # we need to allow for a little slop + if isinstance(getattr(m, p).value, np.longdouble): + assert np.isclose( + getattr(m, p).value, getattr(mback, p).value, atol=1e-15 + ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" + print(p, getattr(m, p).value - getattr(mback, p).value) + else: + assert ( + getattr(m, p).value == getattr(mback, p).value + ), f"{p}: {getattr(m, p).value} does not match {getattr(mback, p).value}" @pytest.mark.parametrize( diff --git a/tests/test_change_epoch.py b/tests/test_change_epoch.py index 41431d2bda..2d799cd21f 100644 --- a/tests/test_change_epoch.py +++ b/tests/test_change_epoch.py @@ -33,7 +33,8 @@ def test_change_pepoch(model): epoch_diff = (t0.mjd_long - model.PEPOCH.quantity.mjd_long) * u.day F0_at_t0 = model.F0.quantity + model.F1.quantity * epoch_diff.to(u.s) model.change_pepoch(t0) - assert model.F0.quantity == F0_at_t0 + # previous test was for equality. Now get slight differences at high precision + assert np.isclose(model.F0.quantity, F0_at_t0, atol=1e-15) def test_change_posepoch(model): diff --git a/tests/test_downhill_fitter.py b/tests/test_downhill_fitter.py index 9574daf34a..05a67fb927 100644 --- a/tests/test_downhill_fitter.py +++ b/tests/test_downhill_fitter.py @@ -48,7 +48,6 @@ def model_eccentric_toas_ecorr(): model_eccentric = get_model( io.StringIO("\n".join([par_eccentric, "ECORR tel @ 2"])) ) - toas = merge_TOAs( [ make_fake_toas_uniform( @@ -70,7 +69,6 @@ def model_eccentric_toas_wb(): model_eccentric = get_model( io.StringIO("\n".join([par_eccentric, "ECORR tel @ 2"])) ) - toas = merge_TOAs( [ make_fake_toas_uniform( diff --git a/tests/test_precision.py b/tests/test_precision.py index 9f6b46544a..c9c1281d41 100644 --- a/tests/test_precision.py +++ b/tests/test_precision.py @@ -39,10 +39,11 @@ time_to_longdouble, time_to_mjd_string, two_sum, + eps, ) from pint.utils import PosVel -time_eps = (np.finfo(float).eps * u.day).to(u.ns) +time_eps = (eps * u.day).to(u.ns) leap_sec_mjds = { @@ -187,7 +188,7 @@ def test_str2longdouble_rejects_bytes(s): ("1", 1), (1.0, 1.0), (np.float32(1.0), 1.0), - (np.longdouble(1.0), 1.0), + (np.float64(1.0).astype(np.longdouble), 1.0), (np.array([1]), np.array([1.0])), (np.array([[1]]), np.array([[1.0]])), ], @@ -205,7 +206,7 @@ def test_data2longdouble_accepts_types(d, ld): [1], [1, 2, 3], [1.5, 2], - np.ones(5, dtype=np.longdouble) + np.finfo(np.longdouble).eps, + np.ones(5, dtype=np.longdouble) + eps, np.random.randn(2, 3, 4), ], ) @@ -477,9 +478,9 @@ def test_pulsar_mjd_never_differs_too_much_from_mjd_utc(i_f): def test_posvel_respects_longdouble(): pos = np.ones(3, dtype=np.longdouble) - pos[0] += np.finfo(np.longdouble).eps + pos[0] += eps vel = np.ones(3, dtype=np.longdouble) - vel[1] += np.finfo(np.longdouble).eps + vel[1] += eps pv = PosVel(pos, vel) assert_array_equal(pv.pos, pos) assert_array_equal(pv.vel, vel) From 7d4a11c400c877b6da198c4b68f638073fe04ed3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 26 May 2026 10:29:53 -0500 Subject: [PATCH 02/11] fix time conversion --- src/pint/pulsar_mjd.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/pint/pulsar_mjd.py b/src/pint/pulsar_mjd.py index 6747e90689..026699cce3 100644 --- a/src/pint/pulsar_mjd.py +++ b/src/pint/pulsar_mjd.py @@ -358,6 +358,8 @@ def data2longdouble(data): return data.astype(longdoubledtype) elif isinstance(data, list): return data2longdouble(np.array(data)) + elif isinstance(data, astropy.time.Time): + return astropy.time.Time(data2longdouble(data.mjd), format="mjd") else: return np.longdouble(data) @@ -459,7 +461,7 @@ def mjds_to_jds_pulsar(mjd1, mjd2): # then back to astropy/ERFA-convention jd1,jd2 using the # ERFA dtf2d() routine which handles leap seconds. v1, v2 = day_frac(mjd1, mjd2) - (y, mo, d, f) = erfa.jd2cal(erfa.DJM0 + v1, v2) + y, mo, d, f = erfa.jd2cal(erfa.DJM0 + v1, v2) # Fractional day to HMS. Uses 86400-second day always. # Seems like there should be a ERFA routine for this.. # There is: erfa.d2tf. Unfortunately it takes a "number of @@ -536,7 +538,7 @@ def str_to_mjds(s): def _mjds_to_str(mjd1, mjd2): - (imjd, fmjd) = day_frac(mjd1, mjd2) + imjd, fmjd = day_frac(mjd1, mjd2) imjd = int(imjd) assert np.fabs(fmjd) < 1.0 while fmjd < 0.0: From 61db44b0a959f1cda2bb3c4baec1b4f7dacb42e6 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 26 May 2026 10:33:15 -0500 Subject: [PATCH 03/11] fix par matching --- tests/test_parameters.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_parameters.py b/tests/test_parameters.py index d144a5a8db..17cbca3716 100644 --- a/tests/test_parameters.py +++ b/tests/test_parameters.py @@ -359,8 +359,8 @@ def test_start_finish_in_par(self): assert hasattr(m1, "FINISH") assert isinstance(m1.FINISH, MJDParameter) - assert m1.START.value == start_preval - assert m1.FINISH.value == finish_preval + assert np.isclose(m1.START.value, start_preval) + assert np.isclose(m1.FINISH.value, finish_preval) assert m1.START.frozen == True assert m1.FINISH.frozen == True From 957bf51883126ec4f993aee6d95b4ed28f3d307e Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 26 May 2026 10:41:51 -0500 Subject: [PATCH 04/11] avoid broken np.average --- tests/test_phase_offset.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_phase_offset.py b/tests/test_phase_offset.py index 61d46bd5fd..3454d1b5a5 100644 --- a/tests/test_phase_offset.py +++ b/tests/test_phase_offset.py @@ -46,8 +46,11 @@ def test_phase_offset(): mean1 = res.calc_phase_mean().value # The new mean should be equal to (0.2 - -0.1). assert np.isclose(mean1 - mean0, 0.3, atol=1e-3) + print(res.calc_phase_resids(), res.calc_phase_resids().__class__) + print(1 / t.get_errors() ** 2, t.get_errors().__class__) assert np.isclose( - np.average(res.calc_phase_resids(), weights=1 / t.get_errors() ** 2), + np.sum(res.calc_phase_resids() / t.get_errors().value ** 2) + / np.sum(1 / t.get_errors().value ** 2), 0.3, atol=1e-3, ) From eac8a187ffa7cf304834be2e47dc8eb013e9e332 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 26 May 2026 12:47:34 -0500 Subject: [PATCH 05/11] float128 -> longdouble --- src/pint/simulation.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/pint/simulation.py b/src/pint/simulation.py index 97de3d3e16..422c6baa84 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -1,5 +1,4 @@ -"""Functions related to simulating TOAs and models -""" +"""Functions related to simulating TOAs and models""" from collections import OrderedDict from copy import deepcopy @@ -587,7 +586,7 @@ def calculate_random_models( Nmjd = len(toas) phases_i = np.zeros((Nmodels, Nmjd)) phases_f = np.zeros((Nmodels, Nmjd)) - freqs = np.zeros((Nmodels, Nmjd), dtype=np.float128) * u.Hz + freqs = np.zeros((Nmodels, Nmjd), dtype=np.longdouble) * u.Hz cov_matrix = fitter.parameter_covariance_matrix # this is a list of the parameter names in the order they appear in the covariance matrix From 7393407df8d12f382dcdc9f7362a8389606ffca2 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 26 May 2026 13:03:34 -0500 Subject: [PATCH 06/11] fixed typing --- src/pint/models/troposphere_delay.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/pint/models/troposphere_delay.py b/src/pint/models/troposphere_delay.py index 3a068709ed..ffa877b886 100644 --- a/src/pint/models/troposphere_delay.py +++ b/src/pint/models/troposphere_delay.py @@ -134,13 +134,13 @@ def _get_target_skycoord(self): """return the sky coordinates for the target, either from equatorial or ecliptic coordinates""" try: radec = SkyCoord( - self._parent.RAJ.value * self._parent.RAJ.units, - self._parent.DECJ.value * self._parent.DECJ.units, + np.float64(self._parent.RAJ.value) * self._parent.RAJ.units, + np.float64(self._parent.DECJ.value) * self._parent.DECJ.units, ) # just do this once instead of adjusting over time except AttributeError: radec = SkyCoord( - self._parent.ELONG.value * self._parent.ELONG.units, - self._parent.ELAT.value * self._parent.ELAT.units, + np.float64(self._parent.ELONG.value) * self._parent.ELONG.units, + np.float64(self._parent.ELAT.value) * self._parent.ELAT.units, frame="barycentricmeanecliptic", ) return radec @@ -156,7 +156,6 @@ def troposphere_delay(self, toas, acc_delay=None): # if not correcting for troposphere, return the default zero delay if self.CORRECT_TROPOSPHERE.value: radec = self._get_target_skycoord() - # the only python for loop is to iterate through the unique observatory locations # all other math is computed through numpy for key, grp in toas.get_obs_groups(): @@ -229,7 +228,7 @@ def delay_model(self, alt, lat, H, mjd): # modify the delay if any of the altitudes are invalid if not np.all(altIsValid): - delay *= altIsValid # this will make the invalid delays zero + delay *= np.float64(altIsValid) # this will make the invalid delays zero return delay def pressure_from_altitude(self, H): From 27bc64e71d98f1eda25f567ee0c5108d0aae27d4 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 27 May 2026 09:53:00 -0500 Subject: [PATCH 07/11] monkeypatched astropy.time.core._make_array --- src/pint/pulsar_mjd.py | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/pint/pulsar_mjd.py b/src/pint/pulsar_mjd.py index 026699cce3..86bf1e79e0 100644 --- a/src/pint/pulsar_mjd.py +++ b/src/pint/pulsar_mjd.py @@ -35,6 +35,7 @@ import numpy as np from astropy.time import Time from astropy.time.formats import TimeFormat +import astropy.time.core import pint try: @@ -103,9 +104,47 @@ def readable_warning(message, category, filename, lineno, line=None): np._core.sctypeDict[str(QuadPrecDType())] = QuadPrecDType() np.longdouble = QuadPrecision - longdoubledtype = QuadPrecDType + longdoubledtype = QuadPrecDType(backend="sleef") eps = npq.epsilon + def _make_array(val, copy=None): + """ + Take ``val`` and convert/reshape to an array. If ``copy`` is `True` + then copy input values. + + Returns + ------- + val : ndarray + Array version of ``val``. + """ + if isinstance(val, (tuple, list)) and len(val) > 0 and isinstance(val[0], Time): + dtype = object + else: + dtype = None + + val = np.array(val, copy=copy, subok=True, dtype=dtype) + + # Allow only float64, string or object arrays as input + # (object is for datetime, maybe add more specific test later?) + # This also ensures the right byteorder for float64 (closes #2942). + # also allow for QuadPrecision now + if ( + val.dtype.kind == "f" + and val.dtype.itemsize >= np.dtype(np.float64).itemsize + ): + pass + elif val.dtype.kind in "OSUMaV": + pass + elif val.dtype == longdoubledtype: + pass + else: + val = np.asanyarray(val, dtype=np.float64) + + return val + + # monkeypatch this into astropy + astropy.time.core._make_array = _make_array + class PulsarMJD(TimeFormat): """Change handling of days with leap seconds. From d8a82177c3ae5755993ad9d59dc2f31a30dae8ec Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 27 May 2026 12:12:13 -0500 Subject: [PATCH 08/11] a few more fixes --- src/pint/models/astrometry.py | 4 +++- src/pint/pulsar_mjd.py | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index c1356f04a3..a4090d0293 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -955,7 +955,9 @@ def get_psr_coords(self, epoch: Optional[time_like] = None) -> coords.SkyCoord: ) # Compute for each time because there is proper motion newepoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") + epoch + if isinstance(epoch, Time) + else Time(epoch.astype(np.float64), scale="tdb", format="mjd") ) position_now = add_dummy_distance(self.get_psr_coords()) with warnings.catch_warnings(): diff --git a/src/pint/pulsar_mjd.py b/src/pint/pulsar_mjd.py index 86bf1e79e0..fa38a8a7fc 100644 --- a/src/pint/pulsar_mjd.py +++ b/src/pint/pulsar_mjd.py @@ -532,15 +532,15 @@ def _str_to_mjds(s): num, expon = ss.split("e") expon = int(expon) if expon < 0: - imjd, fmjd = 0, (ss).astype(np.longdouble) + imjd, fmjd = 0, np.longdouble(ss) else: mjd_s = num.split(".") # If input was given as an integer, add floating "0" if len(mjd_s) == 1: mjd_s.append("0") imjd_s, fmjd_s = mjd_s - imjd = (int(imjd_s)).astype(np.longdouble) - fmjd = (f"0.{fmjd_s}").astype(np.longdouble) + imjd = int(imjd_s) + fmjd = np.longdouble(f"0.{fmjd_s}") if ss.startswith("-"): fmjd = -fmjd imjd *= 10**expon From c8ff5210458f5e2cf7d3f2d0b3694bfe416e1de5 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 27 May 2026 12:39:58 -0500 Subject: [PATCH 09/11] added explicit casting --- src/pint/simulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/simulation.py b/src/pint/simulation.py index 422c6baa84..ea586b7f1d 100644 --- a/src/pint/simulation.py +++ b/src/pint/simulation.py @@ -66,7 +66,7 @@ def zero_residuals( maxresid = np.abs(resids).max() if abs(resids).max() < tolerance: break - ts.adjust_TOAs(-time.TimeDelta(resids)) + ts.adjust_TOAs(-time.TimeDelta(np.float64(resids))) else: raise ValueError( f"Unable to make fake residuals - left over errors are {abs(resids).max()}" @@ -314,7 +314,7 @@ def make_fake_toas_uniform( clk_version = get_fake_toa_clock_versions(model, include_bipm=include_bipm) ts = pint.toa.get_TOAs_array( - times, + np.float64(times), obs=obs, scale=get_observatory(obs).timescale, freqs=freq_array, From 29a7c2dfcbe65ed643fd56ebe83764e7a8d7a0a6 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 27 May 2026 12:42:25 -0500 Subject: [PATCH 10/11] fixed test --- tests/test_wideband_dm_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_wideband_dm_data.py b/tests/test_wideband_dm_data.py index 8083948fba..351bdaf523 100644 --- a/tests/test_wideband_dm_data.py +++ b/tests/test_wideband_dm_data.py @@ -1,4 +1,4 @@ -""" Various of tests on the wideband DM data""" +"""Various of tests on the wideband DM data""" import io import os @@ -67,7 +67,7 @@ def wb_toas_all(wb_model): r = Residuals(toas, wb_model) if np.all(r.time_resids < 1 * u.ns): break - toas.adjust_TOAs(TimeDelta(-r.time_resids)) + toas.adjust_TOAs(TimeDelta(-np.float64(r.time_resids))) else: raise ValueError return toas From 631054246031f5423ba20ed5deaee8f44fb9c9cb Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 27 May 2026 12:50:41 -0500 Subject: [PATCH 11/11] fix precision --- src/pint/fermi_toas.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/fermi_toas.py b/src/pint/fermi_toas.py index 54debd9a20..94436a6c4e 100644 --- a/src/pint/fermi_toas.py +++ b/src/pint/fermi_toas.py @@ -303,12 +303,12 @@ def get_Fermi_TOAs( t = Time( val=mjds[:, 0], val2=mjds[:, 1], - format="mjd", + format="mjd_long", scale=scale, location=location, ) else: - t = Time(mjds, format="mjd", scale=scale, location=location) + t = Time(mjds, format="mjd_long", scale=scale, location=location) if weightcolumn is None: flags = [{"energy": str(e)} for e in energies.to_value(u.MeV)] else: