From 239c0040487f26fe73fdf19452ce11b8e8fbd845 Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Mon, 4 May 2020 10:02:18 -0400
Subject: [PATCH 1/8] Small doc updates

---
 README.rst                  | 2 ++
 docs/basic-installation.rst | 7 +++++++
 2 files changed, 9 insertions(+)

diff --git a/README.rst b/README.rst
index 0cf403681..6973ea20c 100644
--- a/README.rst
+++ b/README.rst
@@ -58,6 +58,8 @@ If you want access to the source code, example notebooks, and tests, you can ins
 cloning the source repository from GitHub, then install
 it, ensuring that all dependencies needed to run PINT are available::
 
+    $ git checkout https://github.com/nanograv/PINT.git
+    $ cd PINT
     $ pip install .
 
 Complete installation instructions are availble here_.
diff --git a/docs/basic-installation.rst b/docs/basic-installation.rst
index 770dcb6be..5fd94dda3 100644
--- a/docs/basic-installation.rst
+++ b/docs/basic-installation.rst
@@ -13,4 +13,11 @@ By default this will install in your system site-packages.  Depending on your sy
 to install it for just yourself (e.g. if you don't have permission to write in the system site-packages), or you may want to create a 
 virtualenv to work on PINT (using a virtualenv is highly recommended by the PINT developers).
 
+If you want access to the source code, example notebooks, and tests, you can install from source, by 
+cloning the source repository from GitHub, then install it, ensuring that all dependencies needed to run PINT are available::
+
+    $ git checkout https://github.com/nanograv/PINT.git
+    $ cd PINT
+    $ pip install .
+
 If this fails, or for more explicit installation and troubleshooting instructions see :ref:`Installation`.

From d10d157be897e2ee8b2a8b10b0bdb2bc1b092960 Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Mon, 4 May 2020 11:29:13 -0400
Subject: [PATCH 2/8] Update author list in docs

---
 AUTHORS.rst | 34 +++++++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/AUTHORS.rst b/AUTHORS.rst
index 87e19c1a6..8c549c176 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -6,24 +6,40 @@ Development of PINT is partially supported by the NANOGrav_ pulsar timing array
 
 .. _NANOGrav: http://nanograv.org/
 
+If you use PINT for your work, please cite the PINT paper and the ASCL entry_ .
+
+.. _entry: http://ascl.net/1902.007
+
 Development Leads
 -----------------
 
-In alphabetical order:
+As listed in the PINT paper:
 
-* Anne Archibald
-* Matteo Bachetti
-* Paul Demorest
-* Justin A. Ellis
 * Luo Jing
-* Rick Jenet
 * Scott Ransom
+* Paul Demorest
 * Paul Ray
-* Chris Sheehy
-* Michele Vallisneri
+* Anne Archibald
+* Matthew Kerr
+* Ross Jennings
+* Matteo Bachetti
 * Rutger van Haasteren
+* Jonathan Colen
+* Camryn Phillips
+* Kevin Stovall
+* Michael Lam
+* Rick Jenet
 
 Other contributors
 ------------------
 
-None yet. Why not be the first?
+* Bastian Beischer
+* Chris Deil
+* Julia Deneva
+* Justin A. Ellis
+* Maarten van Kerkwijk
+* Matt Pitkin
+* Ben Shapiro-Albert
+* Chris Sheehy
+* Renee Spiewak
+* Michele Vallisneri

From da09877e8d4fcd3cecb8482636275c666b612998 Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Tue, 5 May 2020 07:38:45 -0400
Subject: [PATCH 3/8] Add units on returned masses

---
 src/pint/utils.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/pint/utils.py b/src/pint/utils.py
index 1c9f50593..2d77df44f 100644
--- a/src/pint/utils.py
+++ b/src/pint/utils.py
@@ -1101,6 +1101,10 @@ def pulsar_mass(pb, x, mc, inc):
         Companion mass in u.solMass
     inc : Angle
         Inclination angle, in u.deg or u.rad
+
+    Returns
+    -------
+    mass : Quantity in u.solMass
     """
     massfunct = mass_funct(pb, x)
 
@@ -1108,7 +1112,7 @@ def pulsar_mass(pb, x, mc, inc):
     def localmf(mp, mc=mc, mf=massfunct, i=inc):
         return (mass_funct2(mp * u.solMass, mc, i) - mf).value
 
-    return zeros.bisect(localmf, 0.0, 1000.0)
+    return zeros.bisect(localmf, 0.0, 1000.0) * u.solMass
 
 
 def companion_mass(pb, x, inc=60.0 * u.deg, mpsr=1.4 * u.solMass):
@@ -1127,6 +1131,10 @@ def companion_mass(pb, x, inc=60.0 * u.deg, mpsr=1.4 * u.solMass):
         Inclination angle, in u.deg or u.rad. Default is 60 deg.
     mpsr : Quantity, optional
         Pulsar mass in u.solMass. Default is 1.4 Msun
+
+    Returns
+    -------
+    mass : Quantity in u.solMass
     """
     massfunct = mass_funct(pb, x)
 
@@ -1134,7 +1142,7 @@ def companion_mass(pb, x, inc=60.0 * u.deg, mpsr=1.4 * u.solMass):
     def localmf(mc, mp=mpsr, mf=massfunct, i=inc):
         return (mass_funct2(mp, mc * u.solMass, i) - mf).value
 
-    return zeros.bisect(localmf, 0.001, 1000.1)
+    return zeros.bisect(localmf, 0.001, 1000.1) * u.solMass
 
 
 def ELL1_check(A1, E, TRES, NTOA, outstring=True):

From 6fcbd62519200f7b50dfeed88c84de54d8ddab01 Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Tue, 5 May 2020 07:41:15 -0400
Subject: [PATCH 4/8] Improve print statement

---
 src/pint/fitter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pint/fitter.py b/src/pint/fitter.py
index cc31c8ce6..782245148 100644
--- a/src/pint/fitter.py
+++ b/src/pint/fitter.py
@@ -197,7 +197,7 @@ def get_summary(self, nodmx=False):
         s = "Fitted model using {} method with {} free parameters to {} TOAs\n".format(
             self.method, len(self.get_fitparams()), self.toas.ntoas
         )
-        s += "Prefit residuals {}, Postfit residuals {}\n".format(
+        s += "Prefit residuals Wrms = {}, Postfit residuals Wrms = {}\n".format(
             self.resids_init.rms_weighted(), self.resids.rms_weighted()
         )
         s += "Chisq = {:.3f} for {} d.o.f. for reduced Chisq of {:.3f}\n".format(

From acb21c83eb1e391801790fe370be17c281c6de62 Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Wed, 6 May 2020 09:26:07 -0400
Subject: [PATCH 5/8] Clean up unselect

---
 src/pint/toa.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/pint/toa.py b/src/pint/toa.py
index d4005ed8d..edc4b3e8e 100644
--- a/src/pint/toa.py
+++ b/src/pint/toa.py
@@ -971,11 +971,10 @@ def select(self, selectarray):
 
     def unselect(self):
         """Return to previous selected version of the TOA table (stored in stack)."""
-        if hasattr(self, "table_selects"):
-            # may raise an exception about an empty list
+        try:
             self.table = self.table_selects.pop()
-        else:
-            raise ValueError("No previous TOA table found.  No changes made.")
+        except (AttributeError, IndexError) as e:
+            log.error("No previous TOA table found.  No changes made.")
 
     def pickle(self, filename=None):
         """Write the TOAs to a .pickle file with optional filename."""

From 291c2dcaddebf0fadac95bcbbdd875756dbc900c Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Wed, 6 May 2020 09:27:05 -0400
Subject: [PATCH 6/8] Fix bugs in dmx_ranges and add other algorithm

---
 src/pint/utils.py | 169 +++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 143 insertions(+), 26 deletions(-)

diff --git a/src/pint/utils.py b/src/pint/utils.py
index 2d77df44f..37bffe309 100644
--- a/src/pint/utils.py
+++ b/src/pint/utils.py
@@ -487,38 +487,37 @@ def pmtot(model):
 
 
 class dmxrange:
+    """Internal class for building DMX ranges"""
+
     def __init__(self, lofreqs, hifreqs):
+        """lofreqs and hifreqs are lists of MJDs that are in the low or high band respectively"""
         self.los = lofreqs
         self.his = hifreqs
-        self.min = min(lofreqs + hifreqs)
-        self.max = max(lofreqs + hifreqs)
+        self.min = min(lofreqs + hifreqs) - 0.001 * u.d
+        self.max = max(lofreqs + hifreqs) + 0.001 * u.d
 
     def sum_print(self):
         print(
-            "{:8.2f}-{:8.2f} ({:.2f}): ".format(
-                self.min.value, self.max.value, self.max - self.min
-            ),
-            end="",
+            "{:8.2f}-{:8.2f} ({:8.2f}): NLO={:5d} NHI={:5d}".format(
+                self.min.value,
+                self.max.value,
+                self.max - self.min,
+                len(self.los),
+                len(self.his),
+            )
         )
-        print("lofs: ", end="")
-        [print(f.value, end=" ") for f in self.los]
-        print("hifs: ")
-        [print(f.value, end=" ") for f in self.his]
-        print("")
 
 
 def dmx_ranges(
     toas,
     divide_freq=1000.0 * u.MHz,
-    offset=0.1 * u.d,
+    offset=0.01 * u.d,
     max_diff=15.0 * u.d,
     verbose=False,
 ):
     """Compute initial DMX ranges for a set of TOAs
     
-    This is a rudimentary translation of $TEMPO/utils/dmx_ranges/DMX_ranges2.py,
-    which at least reads the MJDs and freqs from a TOAs object but that is it.
-    It should be updated to fill the resulting DMX ranges into a model object.
+    This is a rudimentary translation of $TEMPO/utils/dmx_ranges/DMX_ranges2.py
 
     Parameters
     ----------
@@ -650,7 +649,7 @@ def dmx_ranges(
             dmx_comp.DMX_0001.value = 0.0
             dmx_comp.DMX_0001.frozen = False
             dmx_comp.DMXR1_0001.value = (DMX.min - offset).value
-            dmx_comp.DMXR2_0001.value = (DMX.min + offset).value
+            dmx_comp.DMXR2_0001.value = (DMX.max + offset).value
 
         else:
             # Add the DMX parameters
@@ -674,7 +673,118 @@ def dmx_ranges(
             dmxr2_par = pint.models.parameter.prefixParameter(
                 parameter_type="mjd",
                 name="DMXR2_{:04d}".format(ii + 1),
-                value=(DMX.min + offset).value,
+                value=(DMX.max + offset).value,
+                units=u.d,
+            )
+            dmx_comp.add_param(dmxr2_par, setup=True)
+    # Validate component
+    dmx_comp.validate()
+
+    return mask, dmx_comp
+
+
+def dmx_ranges2(toas, divide_freq=1000.0 * u.MHz, binwidth=15.0 * u.d, verbose=False):
+    """Compute initial DMX ranges for a set of TOAs
+    
+    This is an alternative algorithm for computing DMX ranges
+
+    Parameters
+    ----------
+    divide_freq : Quantity, MHz
+        Requires TOAs above and below this freq for a good DMX range
+    offset : Quantity, days
+        The buffer to include around each DMX range. Warning, may cause bins to overlap?!?
+    max_diff : Quantity, days
+        Maximum duration of a DMX bin
+    verbose : bool
+        If True, print out verbose information about the DMX ranges including par file lines.
+
+    Returns
+    -------
+    mask : bool array
+        Array with True for all TOAs that got assigned to a DMX bin
+    component : TimingModel.Component object
+        A DMX Component class with the DMX ranges included
+    """
+    from pint.models.timing_model import Component
+    import pint.models.parameter
+
+    MJDs = toas.get_mjds()
+    freqs = toas.table["freq"]
+
+    DMXs = []
+
+    prevbinR2 = MJDs[0] - 0.001 * u.d
+    while True:
+        # Consider all TOAs with times after the last bin up through a total span of binwidth
+        # Get indexes that should be in this bin
+        # If there are no more MJDs to process, we are done.
+        if not np.any(MJDs > prevbinR2):
+            break
+        startMJD = MJDs[MJDs > prevbinR2][0]
+        binidx = np.logical_and(MJDs > prevbinR2, MJDs <= startMJD + binwidth)
+        if not np.any(binidx):
+            break
+        binMJDs = MJDs[binidx]
+        binfreqs = freqs[binidx]
+        loMJDs = binMJDs[binfreqs < divide_freq]
+        hiMJDs = binMJDs[binfreqs >= divide_freq]
+        # If we have freqs below and above the divide, this is a good bin
+        if np.any(binfreqs < divide_freq) and np.any(binfreqs > divide_freq):
+            DMXs.append(dmxrange(list(loMJDs), list(hiMJDs)))
+        else:
+            # These TOAs cannot be used
+            pass
+        prevbinR2 = binMJDs.max()
+
+    if verbose:
+        print(
+            "\n These are the good DMX ranges with number of TOAs above/below the dividing freq:"
+        )
+        for DMX in DMXs:
+            DMX.sum_print()
+
+    # Init mask to all False
+    mask = np.zeros_like(MJDs.value, dtype=np.bool)
+    # Mark TOAs as True if they are in any DMX bin
+    for DMX in DMXs:
+        mask[np.logical_and(MJDs >= DMX.min, MJDs <= DMX.max)] = True
+    log.info("{} out of {} TOAs are in a DMX bin".format(mask.sum(), len(mask)))
+    # Instantiate a DMX component
+    dmx_class = Component.component_types["DispersionDMX"]
+    dmx_comp = dmx_class()
+    # Add parameters
+    for ii, DMX in enumerate(DMXs):
+        if ii == 0:
+            # Already have DMX_0001 in component, so just set parameters
+            dmx_comp.DMX_0001.value = 0.0
+            dmx_comp.DMX_0001.frozen = False
+            dmx_comp.DMXR1_0001.value = DMX.min.value
+            dmx_comp.DMXR2_0001.value = DMX.max.value
+
+        else:
+            # Add the DMX parameters
+            dmx_par = pint.models.parameter.prefixParameter(
+                parameter_type="float",
+                name="DMX_{:04d}".format(ii + 1),
+                value=0.0,
+                units=u.pc / u.cm ** 3,
+                frozen=False,
+            )
+            dmx_comp.add_param(dmx_par, setup=True)
+
+            dmxr1_par = pint.models.parameter.prefixParameter(
+                parameter_type="mjd",
+                name="DMXR1_{:04d}".format(ii + 1),
+                value=DMX.min.value,
+                units=u.d,
+            )
+            dmx_comp.add_param(dmxr1_par, setup=True)
+
+            dmxr2_par = pint.models.parameter.prefixParameter(
+                parameter_type="mjd",
+                name="DMXR2_{:04d}".format(ii + 1),
+                value=DMX.max.value,
                 units=u.d,
             )
             dmx_comp.add_param(dmxr2_par, setup=True)
@@ -704,16 +814,23 @@ def dmxstats(fitter):
             mjds.value > getattr(model, "DMXR1_{:04d}".format(ii)).value,
             mjds.value < getattr(model, "DMXR2_{:04d}".format(ii)).value,
         )
-        mjds_in_bin = mjds[mmask]
-        freqs_in_bin = freqs[mmask]
-        span = (mjds_in_bin.max() - mjds_in_bin.min()).to(u.d)
-        # Warning: min() and max() seem to strip the units
-        freqspan = freqs_in_bin.max() - freqs_in_bin.min()
-        print(
-            "DMX_{:04d}: NTOAS={:5d}, MJDSpan={:14.4f}, FreqSpan={:8.3f}-{:8.3f}".format(
-                ii, mmask.sum(), span, freqs_in_bin.min(), freqs_in_bin.max()
+        if np.any(mmask):
+            mjds_in_bin = mjds[mmask]
+            freqs_in_bin = freqs[mmask]
+            span = (mjds_in_bin.max() - mjds_in_bin.min()).to(u.d)
+            # Warning: min() and max() seem to strip the units
+            freqspan = freqs_in_bin.max() - freqs_in_bin.min()
+            print(
+                "DMX_{:04d}: NTOAS={:5d}, MJDSpan={:14.4f}, FreqSpan={:8.3f}-{:8.3f}".format(
+                    ii, mmask.sum(), span, freqs_in_bin.min(), freqs_in_bin.max()
+                )
+            )
+        else:
+            print(
+                "DMX_{:04d}: NTOAS={:5d}, MJDSpan={:14.4f}, FreqSpan={:8.3f}-{:8.3f}".format(
+                    ii, mmask.sum(), 0 * u.d, 0 * u.MHz, 0 * u.MHz
+                )
             )
-        )
         ii += 1
 
     return

From b8954d662bff879fcf2c89a5e6bd4b198f71832f Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Mon, 11 May 2020 16:03:31 -0400
Subject: [PATCH 7/8] Correct POSEPOCH and DMEPOCH to be TDB timescale

---
 src/pint/models/absolute_phase.py   | 4 +++-
 src/pint/models/astrometry.py       | 6 +++++-
 src/pint/models/dispersion_model.py | 4 +++-
 src/pint/models/parameter.py        | 4 ++--
 4 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/src/pint/models/absolute_phase.py b/src/pint/models/absolute_phase.py
index 10ff2ce78..9c27eae53 100644
--- a/src/pint/models/absolute_phase.py
+++ b/src/pint/models/absolute_phase.py
@@ -26,7 +26,9 @@ class AbsPhase(PhaseComponent):
     def __init__(self):
         super(AbsPhase, self).__init__()
         self.add_param(
-            MJDParameter(name="TZRMJD", description="Epoch of the zero phase.")
+            MJDParameter(
+                name="TZRMJD", description="Epoch of the zero phase.", time_scale="utc"
+            )
         )
         self.add_param(
             strParameter(
diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py
index 8e11ce50f..e70d1e0a1 100644
--- a/src/pint/models/astrometry.py
+++ b/src/pint/models/astrometry.py
@@ -31,7 +31,11 @@ class Astrometry(DelayComponent):
     def __init__(self):
         super(Astrometry, self).__init__()
         self.add_param(
-            MJDParameter(name="POSEPOCH", description="Reference epoch for position")
+            MJDParameter(
+                name="POSEPOCH",
+                description="Reference epoch for position",
+                time_scale="tdb",
+            )
         )
 
         self.add_param(
diff --git a/src/pint/models/dispersion_model.py b/src/pint/models/dispersion_model.py
index 07807995c..b2aaf7963 100644
--- a/src/pint/models/dispersion_model.py
+++ b/src/pint/models/dispersion_model.py
@@ -84,7 +84,9 @@ def __init__(self):
             )
         )
         self.add_param(
-            MJDParameter(name="DMEPOCH", description="Epoch of DM measurement")
+            MJDParameter(
+                name="DMEPOCH", description="Epoch of DM measurement", time_scale="tdb"
+            )
         )
 
         self.dm_value_funcs += [self.base_dm]
diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py
index eca269b12..7500bb20e 100644
--- a/src/pint/models/parameter.py
+++ b/src/pint/models/parameter.py
@@ -753,7 +753,7 @@ class MJDParameter(Parameter):
     aliases : list, optional
         An optional list of strings specifying alternate names that can also
         be accepted for this parameter.
-    time_scale : str, optional, default 'utc'
+    time_scale : str, optional, default 'tdb'
         MJD parameter time scale.
 
     Example::
@@ -772,7 +772,7 @@ def __init__(
         frozen=True,
         continuous=True,
         aliases=None,
-        time_scale="utc",
+        time_scale="tdb",
         **kwargs
     ):
         self._time_scale = time_scale

From 8699bb6a46fb1a1bd792af4e0cf4bc560c29931a Mon Sep 17 00:00:00 2001
From: Paul Ray <paul.ray@nrl.navy.mil>
Date: Sun, 3 May 2020 07:12:07 -0400
Subject: [PATCH 8/8] Correct keyword name in notebook

---
 docs/examples/build_model_from_scratch.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/examples/build_model_from_scratch.md b/docs/examples/build_model_from_scratch.md
index f0ab1a30f..69d149a9a 100644
--- a/docs/examples/build_model_from_scratch.md
+++ b/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)