Fix u.cycle bugs on astropy 4.0

requirements.txt

@@ -1,5 +1,5 @@
 numpy>=1.11.0
-astropy>=2.0, <4.0
+astropy>=2.0
 scipy>=0.18.1
 jplephem>=2.6
 matplotlib>=1.5.3

requirements_dev.txt

@@ -3,8 +3,8 @@ pip>=9.0.1
 coverage>=4.3.4
 Sphinx>=1.8.5; python_version < "3.0"
 Sphinx>=2.2; python_version >= "3.0"
-astropy>=2.0,<4.0; python_version < "3.0"
-astropy>=3.2,<4.0; python_version >="3.0"
+astropy>=2.0; python_version < "3.0"
+astropy>=3.2; python_version >="3.0"
 #astropy-helpers>=1.3
 sphinx-rtd-theme>=0.1.9
 coveralls>=1.1

setup.cfg

@@ -26,7 +26,7 @@ package_dir =
 include_package_data = True
 install_requires =
     astropy>=2.0; python_version < "3.0"
-    astropy>=3.2, <4.0; python_version >= "3.0"
+    astropy>=3.2; python_version >= "3.0"
     numpy>=1.11.0
     scipy>=0.18.1
     jplephem>=2.6

src/pint/mcmc_fitter.py

@@ -291,7 +291,8 @@ def get_event_phases(self):
         """
         phases = self.model.phase(self.toas)[1]
         # ensure all positive
-        return np.where(phases < 0.0 * u.cycle, phases + 1.0 * u.cycle, phases)
+        phases = phases.to(u.cycle).value
+        return np.where(phases < 0.0, phases + 1.0, phases)
 
     def lnposterior(self, theta):
         """
@@ -641,7 +642,8 @@ def get_event_phases(self, index=None):
             print("Showing all %d phases" % len(phases))
         else:
             phases = self.model.phase(self.toas_list[index])[1]
-        return np.where(phases < 0.0 * u.cycle, phases + 1.0 * u.cycle, phases)
+        phases = phases.to(u.cycle).value
+        return np.where(phases < 0.0, phases + 1.0, phases)
 
     def get_template_vals(self, phases, index):
         if self.templates[index] is None:

src/pint/models/ifunc.py

@@ -11,7 +11,7 @@ class IFunc(PhaseComponent):
     """This class implements tabulated delays.
 
     These mimic a tempo2 feature, which supports piecewise, linear, and sinc
-    interpolation.  The implementation here currently only supports the 
+    interpolation.  The implementation here currently only supports the
     first two formulae.
 
     For consistency with tempo2, although the IFuncs represent time series,
@@ -29,7 +29,7 @@ class IFunc(PhaseComponent):
     0 == piecewise (no interpolation)
     1 == sinc (not supported)
     2 == linear
-    
+
     NB that the trailing 0.0s are necessary for accurate tempo2 parsing.
     NB also that tempo2 has a static setting MAX_IFUNC whose default value
     is 1000.
@@ -130,5 +130,5 @@ def ifunc_phase(self, toas, delays):
         else:
             raise ValueError("Interpolation type %d not supported.".format(itype))
 
-        phase = ((times * u.s) * self.F0.quantity * 2 * np.pi).to(u.cycle)
+        phase = (times * u.s) * self.F0.quantity * u.cycle
         return phase

src/pint/models/solar_wind_dispersion.py

@@ -62,7 +62,7 @@ def solar_wind_delay(self, toas, acc_delay=None):
             rsa = tbl["obs_sun_pos"].quantity
             pos = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"].astype(np.float64))
             r = np.sqrt(np.sum(rsa * rsa, axis=1))
-            cos_theta = np.sum(rsa * pos, axis=1) / r
+            cos_theta = (np.sum(rsa * pos, axis=1) / r).to(u.Unit("")).value
             ret = (
                 const.au ** 2.0
                 * np.arccos(cos_theta)

src/pint/models/stand_alone_psr_binaries/DDK_model.py

@@ -5,13 +5,10 @@
 import astropy.units as u
 import numpy as np
 from astropy import log
-
 from pint import GMsun, Tsun, ls
 
 from .DD_model import DDmodel
 
-u.set_enabled_equivalencies(u.dimensionless_angles())
-
 
 class DDKmodel(DDmodel):
     """DDK model, a Kopeikin method corrected DD model.

src/pint/models/stand_alone_psr_binaries/ELL1H_model.py

@@ -414,11 +414,13 @@ def d_delayS_d_par(self, par):
         d_Phi_d_par = self.prtl_der("Phi", par)
         d_STIGMA_d_par = self.prtl_der("STIGMA", par)
 
-        return (
-            d_delayS_d_H3 * d_H3_d_par
-            + d_delayS_d_Phi * d_Phi_d_par
-            + d_delayS_d_STIGMA * d_STIGMA_d_par
-        )
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            d_delayS_d_par = (
+                d_delayS_d_H3 * d_H3_d_par
+                + d_delayS_d_Phi * d_Phi_d_par
+                + d_delayS_d_STIGMA * d_STIGMA_d_par
+            )
+        return d_delayS_d_par
 
     def d_ELL1Hdelay_d_par(self, par):
         return self.d_delayI_d_par(par) + self.d_delayS_d_par(par)

src/pint/models/stand_alone_psr_binaries/ELL1_model.py

@@ -166,18 +166,20 @@ def d_Dre_d_par(self, par):
         d_Dre_d_eps1 = a1 / c.c * (-0.5 * np.cos(2 * Phi))
         d_Dre_d_eps2 = a1 / c.c * (0.5 * np.sin(2 * Phi))
 
-        return (
-            d_a1_d_par
-            / c.c
-            * (
-                np.sin(Phi)
-                - 0.5 * eps1 * np.cos(2 * Phi)
-                + 0.5 * eps2 * np.sin(2 * Phi)
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            d_Dre_d_par = (
+                d_a1_d_par
+                / c.c
+                * (
+                    np.sin(Phi)
+                    - 0.5 * eps1 * np.cos(2 * Phi)
+                    + 0.5 * eps2 * np.sin(2 * Phi)
+                )
+                + d_Dre_d_Phi * d_Phi_d_par
+                + d_Dre_d_eps1 * self.prtl_der("eps1", par)
+                + d_Dre_d_eps2 * self.prtl_der("eps2", par)
             )
-            + d_Dre_d_Phi * d_Phi_d_par
-            + d_Dre_d_eps1 * self.prtl_der("eps1", par)
-            + d_Dre_d_eps2 * self.prtl_der("eps2", par)
-        )
+        return d_Dre_d_par
 
     def Drep(self):
         """ dDre/dPhi
@@ -215,14 +217,16 @@ def d_Drep_d_par(self, par):
         d_Drep_d_eps1 = a1 / c.c * np.sin(2.0 * Phi)
         d_Drep_d_eps2 = a1 / c.c * np.cos(2.0 * Phi)
 
-        return (
-            d_a1_d_par
-            / c.c
-            * (np.cos(Phi) + eps1 * np.sin(2.0 * Phi) + eps2 * np.cos(2.0 * Phi))
-            + d_Drep_d_Phi * d_Phi_d_par
-            + d_Drep_d_eps1 * self.prtl_der("eps1", par)
-            + d_Drep_d_eps2 * self.prtl_der("eps2", par)
-        )
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            d_Drep_d_par = (
+                d_a1_d_par
+                / c.c
+                * (np.cos(Phi) + eps1 * np.sin(2.0 * Phi) + eps2 * np.cos(2.0 * Phi))
+                + d_Drep_d_Phi * d_Phi_d_par
+                + d_Drep_d_eps1 * self.prtl_der("eps1", par)
+                + d_Drep_d_eps2 * self.prtl_der("eps2", par)
+            )
+        return d_Drep_d_par
 
     def Drepp(self):
         a1 = self.a1()
@@ -265,17 +269,19 @@ def d_Drepp_d_par(self, par):
         d_Drepp_d_eps1 = a1 / c.c * 2.0 * np.cos(2.0 * Phi)
         d_Drepp_d_eps2 = -a1 / c.c * 2.0 * np.sin(2.0 * Phi)
 
-        return (
-            d_a1_d_par
-            / c.c
-            * (
-                -np.sin(Phi)
-                + 2.0 * (eps1 * np.cos(2.0 * Phi) - eps2 * np.sin(2.0 * Phi))
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            d_Drepp_d_par = (
+                d_a1_d_par
+                / c.c
+                * (
+                    -np.sin(Phi)
+                    + 2.0 * (eps1 * np.cos(2.0 * Phi) - eps2 * np.sin(2.0 * Phi))
+                )
+                + d_Drepp_d_Phi * d_Phi_d_par
+                + d_Drepp_d_eps1 * self.prtl_der("eps1", par)
+                + d_Drepp_d_eps2 * self.prtl_der("eps2", par)
             )
-            + d_Drepp_d_Phi * d_Phi_d_par
-            + d_Drepp_d_eps1 * self.prtl_der("eps1", par)
-            + d_Drepp_d_eps2 * self.prtl_der("eps2", par)
-        )
+        return d_Drepp_d_par
 
     def delayR(self):
         """ELL1 Roemer delay in proper time. Ch. Lange,1 F. Camilo, 2001 eq. A6 """
@@ -407,11 +413,13 @@ def d_delayS_d_par(self, par):
         d_SINI_d_par = self.prtl_der("SINI", par)
         d_Phi_d_par = self.prtl_der("Phi", par)
 
-        return (
-            d_delayS_d_TM2 * d_TM2_d_par
-            + d_delayS_d_SINI * d_SINI_d_par
-            + d_delayS_d_Phi * d_Phi_d_par
-        )
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            d_delayS_d_par = (
+                d_delayS_d_TM2 * d_TM2_d_par
+                + d_delayS_d_SINI * d_SINI_d_par
+                + d_delayS_d_Phi * d_Phi_d_par
+            )
+        return d_delayS_d_par
 
     def ELL1delay(self):
         # TODO need add aberration delay

src/pint/models/stand_alone_psr_binaries/binary_generic.py

@@ -483,7 +483,9 @@ def d_M_d_par(self, par):
 
         par_obj = getattr(self, par)
         result = self.orbits_cls.d_orbits_d_par(par)
-        return result.to(u.Unit("") / par_obj.unit)
+        with u.set_enabled_equivalencies(u.dimensionless_angles()):
+            result = result.to(u.Unit("") / par_obj.unit)
+        return result
 
     ###############################################
 

src/pint/models/wave.py

@@ -108,5 +108,5 @@ def wave_phase(self, toas, delays):
             times += wave_a * np.sin(wave_phase)
             times += wave_b * np.cos(wave_phase)
 
-        phase = ((times) * self.F0.quantity * 2 * np.pi).to(u.cycle)
+        phase = (times) * self.F0.quantity * u.cycle
         return phase

src/pint/scripts/event_optimize.py

@@ -279,7 +279,8 @@ def get_event_phases(self):
         """
         phss = self.model.phase(self.toas)[1]
         # ensure all postive
-        return np.where(phss < 0.0 * u.cycle, phss + 1.0 * u.cycle, phss)
+        phss = phss.to(u.cycle).value
+        return np.where(phss < 0.0, phss + 1.0, phss)
 
     def lnprior(self, theta):
         """

src/pint/scripts/fermiphase.py

@@ -117,7 +117,8 @@ def main(argv=None):
     # Compute model phase for each TOA
     iphss, phss = modelin.phase(ts, abs_phase=True)
     # ensure all postive
-    phases = np.where(phss < 0.0 * u.cycle, phss + 1.0 * u.cycle, phss)
+    phss = phss.to(u.cycle).value
+    phases = np.where(phss < 0.0, phss + 1.0, phss)
     mjds = ts.get_mjds()
     weights = np.array([w["weight"] for w in ts.table["flags"]])
     h = float(hmw(phases, weights))

src/pint/scripts/photonphase.py

@@ -199,8 +199,9 @@ def main(argv=None):
     # Compute model phase for each TOA
     iphss, phss = modelin.phase(ts, abs_phase=True)
     # ensure all postive
-    negmask = phss < 0.0 * u.cycle
-    phases = np.where(negmask, phss + 1.0 * u.cycle, phss)
+    phss = phss.to(u.cycle).value
+    negmask = phss < 0.0
+    phases = np.where(negmask, phss + 1.0, phss)
     h = float(hm(phases))
     print("Htest : {0:.2f} ({1:.2f} sigma)".format(h, h2sig(h)))
     if args.plot:

tests/test_precision.py

@@ -483,6 +483,7 @@ def test_posvel_respects_longdouble():
 @example(i_f=(43875, -1.000000000000002))
 @example(i_f=(48803, 1.0769154457079824e-09))
 @example(i_f=(48803, 1.0000079160299438e-06))
+@settings(deadline=1000)
 @pytest.mark.parametrize("format", ["pulsar_mjd", "mjd"])
 def test_time_from_mjd_string_accuracy_vs_longdouble(format, i_f):
     i, f = i_f