WIP: Incorporate quaddtype

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)

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:

src/pint/fitter.py

@@ -1256,10 +1256,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
         )
@@ -2224,16 +2224,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:
@@ -2244,12 +2244,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
@@ -2366,7 +2366,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 / norm
 
@@ -2582,6 +2582,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
@@ -2610,8 +2611,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
 
 
@@ -2661,8 +2662,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
 
 

src/pint/models/astrometry.py

@@ -447,10 +447,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,
             )
@@ -576,17 +576,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]
@@ -771,19 +773,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))
@@ -792,19 +794,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))
@@ -944,16 +946,18 @@ 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,
             )
             # 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():
@@ -1079,16 +1083,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]
@@ -1326,20 +1330,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))
@@ -1348,20 +1352,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))
@@ -1405,20 +1409,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)
@@ -1428,20 +1432,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)

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

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"""