Remove u.cycle units from phases

CHANGELOG.md

@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project, at least loosely, adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Changed
+- Changed units of Phase to be u.dimensionless_unscaled instead of u.cycle, which was confusing
 
 ## [0.6.3] - 2020-05-04
 ### Added

docs/examples/build_model_from_scratch.md

@@ -271,7 +271,7 @@ The prefix type of parameters have to use `prefixParameter` class from `pint.mod
 ```python
 # Add prefix parameters
 dmx_0003 = p.prefixParameter(
-    parameter_type="float", name="DMX_0003", value=None, unit=u.pc / u.cm ** 3
+    parameter_type="float", name="DMX_0003", value=None, units=u.pc / u.cm ** 3
 )
 
 tm.components["DispersionDMX"].add_param(dmx_0003, setup=True)
@@ -290,10 +290,10 @@ However only adding DMX_0003 is not enough, since one DMX parameter also need a
 
 ```python
 dmxr1_0003 = p.prefixParameter(
-    parameter_type="mjd", name="DMXR1_0003", value=None, unit=u.day
+    parameter_type="mjd", name="DMXR1_0003", value=None, units=u.day
 )  # DMXR1 is a type of MJD parameter internally.
 dmxr2_0003 = p.prefixParameter(
-    parameter_type="mjd", name="DMXR2_0003", value=None, unit=u.day
+    parameter_type="mjd", name="DMXR2_0003", value=None, units=u.day
 )  # DMXR1 is a type of MJD parameter internally.
 
 tm.components["DispersionDMX"].add_param(dmxr1_0003, setup=True)

src/pint/__init__.py

@@ -17,7 +17,6 @@
     "ls",
     "dmu",
     "light_second_equivalency",
-    "dimensionless_cycles",
     "hourangle_second",
     "pulsar_mjd",
     "GMsun",
@@ -50,7 +49,6 @@
 
 # define equivalency for astropy units
 light_second_equivalency = [(ls, si.second, lambda x: x, lambda x: x)]
-dimensionless_cycles = [(u.cycle, None)]
 # hourangle_second unit
 hourangle_second = u.def_unit("hourangle_second", u.hourangle / np.longdouble(3600.0))
 
@@ -84,7 +82,7 @@
     "Tsun": Tsun,
     "GMsun": GMsun,
     "MJD": u.day,
-    "pulse phase": u.cycle,
+    "pulse phase": u.dimensionless_unscaled,
     "hourangle_second": hourangle_second,
 }
 

src/pint/mcmc_fitter.py

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

src/pint/models/glitch.py

@@ -5,8 +5,8 @@
 
 import astropy.units as u
 import numpy as np
+from astropy import log
 
-from pint import dimensionless_cycles
 from pint.models.parameter import prefixParameter
 from pint.models.timing_model import MissingParameter, PhaseComponent
 from pint.utils import split_prefixed_name
@@ -155,43 +155,39 @@ def glitch_phase(self, toas, delay):
         returns an array of phases in long double
         """
         tbl = toas.table
-        phs = np.zeros_like(tbl, dtype=np.longdouble) * u.cycle
+        phs = u.Quantity(np.zeros_like(tbl, dtype=np.longdouble))
         glepnames = [x for x in self.params if x.startswith("GLEP_")]
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            for glepnm in glepnames:
-                glep = getattr(self, glepnm)
-                eph = glep.value
-                idx = glep.index
-                dphs = getattr(self, "GLPH_%d" % idx).quantity
-                dF0 = getattr(self, "GLF0_%d" % idx).quantity
-                dF1 = getattr(self, "GLF1_%d" % idx).quantity
-                dF2 = getattr(self, "GLF2_%d" % idx).quantity
-                dt = (tbl["tdbld"] - eph) * u.day - delay
-                dt = dt.to(u.second)
-                affected = dt > 0.0  # TOAs affected by glitch
-                # decay term
-                dF0D = getattr(self, "GLF0D_%d" % idx).quantity
-                if dF0D != 0.0:
-                    tau = getattr(self, "GLTD_%d" % idx).quantity
-                    decayterm = (
-                        dF0D
-                        * tau
-                        * (1.0 - np.exp(-(dt[affected] / tau).to(u.Unit(""))))
-                    )
-                else:
-                    decayterm = 0.0
+        for glepnm in glepnames:
+            glep = getattr(self, glepnm)
+            eph = glep.value
+            idx = glep.index
+            dphs = getattr(self, "GLPH_%d" % idx).quantity
+            dF0 = getattr(self, "GLF0_%d" % idx).quantity
+            dF1 = getattr(self, "GLF1_%d" % idx).quantity
+            dF2 = getattr(self, "GLF2_%d" % idx).quantity
+            dt = (tbl["tdbld"] - eph) * u.day - delay
+            dt = dt.to(u.second)
+            affected = dt > 0.0  # TOAs affected by glitch
+            # decay term
+            dF0D = getattr(self, "GLF0D_%d" % idx).quantity
+            if dF0D != 0.0:
+                tau = getattr(self, "GLTD_%d" % idx).quantity
+                decayterm = dF0D * tau * (1.0 - np.exp(-(dt[affected] / tau)))
+            else:
+                decayterm = 0.0 * u.Unit("")
 
-                phs[affected] += (
-                    dphs
-                    + dt[affected]
-                    * (
-                        dF0
-                        + 0.5 * dt[affected] * dF1
-                        + 1.0 / 6.0 * dt[affected] * dt[affected] * dF2
-                    )
-                    + decayterm
+            log.info("{} {} ".format(dphs, dphs.unit))
+            phs[affected] += (
+                dphs
+                + dt[affected]
+                * (
+                    dF0
+                    + 0.5 * dt[affected] * dF1
+                    + 1.0 / 6.0 * dt[affected] * dt[affected] * dF2
                 )
-            return phs.to(u.cycle)
+                + decayterm
+            )
+        return phs
 
     def d_phase_d_GLPH(self, toas, param, delay):
         """Calculate the derivative wrt GLPH"""
@@ -206,8 +202,8 @@ def d_phase_d_GLPH(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLPH = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLPH.units
-        dpdGLPH[affected] += 1.0 * u.cycle / par_GLPH.units
+        dpdGLPH = np.zeros(len(tbl), dtype=np.longdouble) / par_GLPH.units
+        dpdGLPH[affected] += 1.0 / par_GLPH.units
         return dpdGLPH
 
     def d_phase_d_GLF0(self, toas, param, delay):
@@ -223,9 +219,8 @@ def d_phase_d_GLF0(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF0 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF0.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF0[affected] = dt[affected]
+        dpdGLF0 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0.units
+        dpdGLF0[affected] = dt[affected]
         return dpdGLF0
 
     def d_phase_d_GLF1(self, toas, param, delay):
@@ -241,9 +236,8 @@ def d_phase_d_GLF1(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF1 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF1.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF1[affected] += np.longdouble(0.5) * dt[affected] * dt[affected]
+        dpdGLF1 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF1.units
+        dpdGLF1[affected] += np.longdouble(0.5) * dt[affected] * dt[affected]
         return dpdGLF1
 
     def d_phase_d_GLF2(self, toas, param, delay):
@@ -259,11 +253,10 @@ def d_phase_d_GLF2(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF2 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF2.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF2[affected] += (
-                np.longdouble(1.0) / 6.0 * dt[affected] * dt[affected] * dt[affected]
-            )
+        dpdGLF2 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF2.units
+        dpdGLF2[affected] += (
+            np.longdouble(1.0) / 6.0 * dt[affected] * dt[affected] * dt[affected]
+        )
         return dpdGLF2
 
     def d_phase_d_GLF0D(self, toas, param, delay):
@@ -281,11 +274,8 @@ def d_phase_d_GLF0D(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF0D = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF0D.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF0D[affected] += tau * (
-                np.longdouble(1.0) - np.exp(-dt[affected] / tau)
-            )
+        dpdGLF0D = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0D.units
+        dpdGLF0D[affected] += tau * (np.longdouble(1.0) - np.exp(-dt[affected] / tau))
         return dpdGLF0D
 
     def d_phase_d_GLTD(self, toas, param, delay):
@@ -300,17 +290,14 @@ def d_phase_d_GLTD(self, toas, param, delay):
         eph = np.longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLTD = getattr(self, param)
         if par_GLTD.value == 0.0:
-            return np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLTD.units
+            return np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units
         glf0d = getattr(self, "GLF0D_" + ids).quantity
         tau = par_GLTD.quantity
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLTD = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLTD.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLTD[affected] += glf0d * (
-                np.longdouble(1.0) - np.exp(-dt[affected] / tau)
-            ) + glf0d * tau * (-np.exp(-dt[affected] / tau)) * dt[affected] / (
-                tau * tau
-            )
+        dpdGLTD = np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units
+        dpdGLTD[affected] += glf0d * (
+            np.longdouble(1.0) - np.exp(-dt[affected] / tau)
+        ) + glf0d * tau * (-np.exp(-dt[affected] / tau)) * dt[affected] / (tau * tau)
         return dpdGLTD

src/pint/models/ifunc.py

@@ -132,5 +132,5 @@ def ifunc_phase(self, toas, delays):
         else:
             raise ValueError("Interpolation type %d not supported.".format(itype))
 
-        phase = (times * u.s) * self.F0.quantity * u.cycle
+        phase = ((times * u.s) * self.F0.quantity).to(u.dimensionless_unscaled)
         return phase

src/pint/models/jump.py

@@ -6,7 +6,6 @@
 import astropy.units as u
 import numpy
 
-from pint import dimensionless_cycles
 from pint.models.parameter import maskParameter
 from pint.models.timing_model import DelayComponent, MissingParameter, PhaseComponent
 
@@ -120,8 +119,7 @@ def d_phase_d_jump(self, toas, jump_param, delay):
         d_phase_d_j = numpy.zeros(len(tbl))
         mask = jpar.select_toa_mask(toas)
         d_phase_d_j[mask] = self.F0.value
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            return (d_phase_d_j * self.F0.units).to(u.cycle / u.second)
+        return (d_phase_d_j * self.F0.units).to(1 / u.second)
 
     def print_par(self):
         result = ""

src/pint/models/spindown.py

@@ -7,7 +7,6 @@
 import numpy
 
 import pint.toa as toa
-from pint import dimensionless_cycles
 from pint.models.parameter import MJDParameter, floatParameter, prefixParameter
 from pint.models.timing_model import MissingParameter, PhaseComponent
 from pint.pulsar_mjd import Time
@@ -127,10 +126,9 @@ def spindown_phase(self, toas, delay):
         """
         dt = self.get_dt(toas, delay)
         # Add the [0.0] because that is the constant phase term
-        fterms = [0.0 * u.cycle] + self.get_spin_terms()
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            phs = taylor_horner(dt.to(u.second), fterms)
-            return phs.to(u.cycle)
+        fterms = [0.0 * u.dimensionless_unscaled] + self.get_spin_terms()
+        phs = taylor_horner(dt.to(u.second), fterms)
+        return phs.to(u.dimensionless_unscaled)
 
     def change_pepoch(self, new_epoch, toas=None, delay=None):
         """Move PEPOCH to a new time and change the related paramters.
@@ -197,13 +195,11 @@ def d_phase_d_F(self, toas, param, delay):
         fterms = [ft * numpy.longdouble(0.0) / unit for ft in fterms]
         fterms[order] += numpy.longdouble(1.0)
         dt = self.get_dt(toas, delay)
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            d_pphs_d_f = taylor_horner(dt.to(u.second), fterms)
-            return d_pphs_d_f.to(u.cycle / unit)
+        d_pphs_d_f = taylor_horner(dt.to(u.second), fterms)
+        return d_pphs_d_f.to(1 / unit)
 
     def d_spindown_phase_d_delay(self, toas, delay):
         dt = self.get_dt(toas, delay)
         fterms = [0.0] + self.get_spin_terms()
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            d_pphs_d_delay = taylor_horner_deriv(dt.to(u.second), fterms)
-            return -d_pphs_d_delay.to(u.cycle / u.second)
+        d_pphs_d_delay = taylor_horner_deriv(dt.to(u.second), fterms)
+        return -d_pphs_d_delay.to(1 / u.second)

src/pint/models/timing_model.py

@@ -16,7 +16,6 @@
 from astropy import log
 
 import pint
-from pint import dimensionless_cycles
 from pint.models.parameter import strParameter, maskParameter
 from pint.phase import Phase
 from pint.utils import PrefixError, interesting_lines, lines_of, split_prefixed_name
@@ -906,8 +905,7 @@ def d_phase_d_toa(self, toas, sample_step=None):
         dp = sample_phase[1] - sample_phase[0]
         d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step)
         del copy_toas
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            return d_phase_d_toa.to(u.Hz)
+        return d_phase_d_toa.to(u.Hz)
 
     def d_phase_d_tpulsar(self, toas):
         """Return the derivative of phase wrt time at the pulsar.
@@ -922,7 +920,7 @@ def d_phase_d_param(self, toas, delay, param):
         # Is it safe to assume that any param affecting delay only affects
         # phase indirectly (and vice-versa)??
         par = getattr(self, param)
-        result = np.longdouble(np.zeros(toas.ntoas)) * u.cycle / par.units
+        result = np.longdouble(np.zeros(toas.ntoas)) / par.units
         param_phase_derivs = []
         phase_derivs = self.phase_deriv_funcs
         delay_derivs = self.delay_deriv_funcs
@@ -940,7 +938,7 @@ def d_phase_d_param(self, toas, delay, param):
             #                         d_delay_d_param
 
             d_delay_d_p = self.d_delay_d_param(toas, param)
-            dpdd_result = np.longdouble(np.zeros(toas.ntoas)) * u.cycle / u.second
+            dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second
             for dpddf in self.d_phase_d_delay_funcs:
                 dpdd_result += dpddf(toas, delay)
             result = dpdd_result * d_delay_d_p
@@ -978,8 +976,12 @@ def d_phase_d_param_num(self, toas, param, step=1e-2):
             h = ori_value * step
         parv = [par.value - h, par.value + h]
 
-        phase_i = np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.cycle
-        phase_f = np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.cycle
+        phase_i = (
+            np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled
+        )
+        phase_f = (
+            np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled
+        )
         for ii, val in enumerate(parv):
             par.value = val
             ph = self.phase(toas)

src/pint/models/wave.py

@@ -111,5 +111,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 * u.cycle
+        phase = ((times) * self.F0.quantity).to(u.dimensionless_unscaled)
         return phase