Merge branch 'master' into chromatic

AUTHORS.rst

@@ -20,7 +20,7 @@ Active developers are indicated by (*). Authors of the PINT paper are indicated
 * Anne Archibald (#*)
 * Matteo Bachetti (#)
 * Bastian Beischer 
-* Deven Bhakta 
+* Deven Bhakta (*)
 * Chloe Champagne (#)
 * Jonathan Colen (#)
 * Thankful Cromartie 

CHANGELOG-unreleased.md

@@ -9,18 +9,7 @@ the released changes.
 
 ## Unreleased
 ### Changed
-- Third-order Roemer delay terms to ELL1 model
-- Made the addition of a TZR TOA (`AbsPhase`) in the `TimingModel` explicit in `Residuals` class.
-- Updated `CONTRIBUTING.rst` with the latest information.
-- Made `TimingModel.params` and `TimingModel.ordered_params` identical.
 ### Added
-- Third-order Roemer delay terms to ELL1 model
-- Options to add a TZR TOA (`AbsPhase`) during the creation of a `TimingModel` using `ModelBuilder.__call__`, `get_model`, and `get_model_and_toas`
-- `pint.print_info()` function for bug reporting
 ### Fixed
-- Deleting JUMP1 from flag tables will not prevent fitting
-- Simulating TOAs from tim file when PLANET_SHAPIRO is true now works
-- Docstrings for `get_toas()` and `get_model_and_toas()`
-- Set `DelayComponent_list` and `NoiseComponent_list` to empty list if such components are absent
-- Fix invalid access of `PLANET_SHAPIRO` in models without `Astrometry`
+- Fixed RTD by specifying theme explicitly.
 ### Removed

CHANGELOG.md

@@ -7,13 +7,36 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 This file contains the released changes to the codebase. See CHANGELOG-unreleased.md for
 the unreleased changes. This file should only be changed while tagging a new version.
 
+## [0.9.7] 2023-08-24
+### Changed
+- Third-order Roemer delay terms to ELL1 model
+- Made the addition of a TZR TOA (`AbsPhase`) in the `TimingModel` explicit in `Residuals` class.
+- Updated `CONTRIBUTING.rst` with the latest information.
+- Made `TimingModel.params` and `TimingModel.ordered_params` identical. Deprecated `TimingModel.ordered_params`.
+### Added
+- Third-order Roemer delay terms to ELL1 model
+- Options to add a TZR TOA (`AbsPhase`) during the creation of a `TimingModel` using `ModelBuilder.__call__`, `get_model`, and `get_model_and_toas`
+- `pint.print_info()` function for bug reporting
+- Added an autocorrelation function to check for chain convergence in `event_optimize`
+- A hacky implementation of system-dependent FD parameters (FDJUMP)
+- Minor doc updates to explain default NHARMS and missing derivative functions
+### Fixed
+- Deleting JUMP1 from flag tables will not prevent fitting
+- Simulating TOAs from tim file when PLANET_SHAPIRO is true now works
+- Docstrings for `get_toas()` and `get_model_and_toas()`
+- Set `DelayComponent_list` and `NoiseComponent_list` to empty list if such components are absent
+- Fix invalid access of `PLANET_SHAPIRO` in models without `Astrometry`
+### Removed
+
+
 ## [0.9.6] 2023-06-22
 ### Changed
 - Applied `sourcery` refactors to the entire codebase
 - Changed threshold for `test_model_derivatives` test to avoid CI failures
 - Unreleased CHANGELOG entries should now be entered in `CHANGELOG-unreleased.md` instead of `CHANGELOG.md`. Updated documentation accordingly.
 - Changed tests to remove `unittest` and use pure pytest format
 - Changed deprecated `sampler.chain` usage
+- Download data automatically in the profiling script `high_level_benchmark.py` instead of silently giving wrong results.
 ### Added
 - `SpindownBase` as the abstract base class for `Spindown` and `PeriodSpindown` in the `How_to_build_a_timing_model_component.py` example.
 - `SolarWindDispersionBase` as the abstract base class for solar wind dispersion components.

docs/conf.py

@@ -182,8 +182,10 @@ def setup(app):
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 on_rtd = os.environ.get("READTHEDOCS", None) == "True"
-if not on_rtd:  # only import and set the theme if we're building docs locally
-    html_theme = "sphinx_rtd_theme"
+# if not on_rtd:  # only import and set the theme if we're building docs locally
+#     html_theme = "sphinx_rtd_theme"
+# it seems that we need to specify it anyway
+html_theme = "sphinx_rtd_theme"
 
 # Theme options are theme-specific and customize the look and feel of a
 # theme further.  For a list of options available for each theme, see the

profiling/.gitignore

@@ -0,0 +1,2 @@
+J0740+6620.cfr+19.tim
+bench_*_summary

profiling/high_level_benchmark.py

@@ -16,6 +16,7 @@
 import sys
 import os
 import platform
+import urllib.request
 from prfparser import parse_file
 
 
@@ -67,11 +68,19 @@ def get_results(script, outfile):
     parser = argparse.ArgumentParser(
         description="High-level summary of python file timing."
     )
+
+    if not os.path.isfile("J0740+6620.cfr+19.tim"):
+        print("Downloading data file J0740+6620.cfr+19.tim ...")
+        urllib.request.urlretrieve(
+            "https://data.nanograv.org/static/data/J0740+6620.cfr+19.tim",
+            "J0740+6620.cfr+19.tim",
+        )
+
+    script1 = "bench_load_TOAs.py"
     script2 = "bench_chisq_grid.py"
     script3 = "bench_chisq_grid_WLSFitter.py"
     script4 = "bench_MCMC.py"
 
-    script1 = "bench_load_TOAs.py"
     # time scripts
     output1 = bench_file(script1)
     output2 = bench_file(script2)

requirements_dev.txt

@@ -3,10 +3,10 @@ pip>=9.0.1
 setuptools>=41.0
 coverage>=4.3.4
 traitlets
-Sphinx>=2.2,!=5.1.0
+Sphinx==6.2.1
 astropy>=4.0,!=4.0.1,!=4.0.1.post1
 #astropy-helpers>=1.3
-sphinx-rtd-theme>=1.1.1
+sphinx-rtd-theme==1.2.2
 coveralls>=1.1
 wheel>=0.29.0
 pytest>=4.3

src/pint/binaryconvert.py

@@ -560,6 +560,10 @@ def convert_binary(model, output, NHARMS=3, useSTIGMA=False, KOM=0 * u.deg):
     Returns
     -------
     outmodel : pint.models.timing_model.TimingModel
+
+    Notes
+    -----
+    Default value in `pint` for `NHARMS` is 7, while in `tempo2` it is 4.
     """
     # Do initial checks
     if output not in binary_types:

src/pint/fitter.py

@@ -363,7 +363,7 @@ def get_summary(self, nodmx=False):
 
         # to handle all parameter names, determine the longest length for the first column
         longestName = 0  # optionally specify the minimum length here instead of 0
-        for pn in self.model.params_ordered:
+        for pn in self.model.params:
             if nodmx and pn.startswith("DMX"):
                 continue
             if len(pn) > longestName:
@@ -378,7 +378,7 @@ def get_summary(self, nodmx=False):
         s += ("{:<" + spacingName + "s} {:>20s} {:>28s} {}\n").format(
             "=" * longestName, "=" * 20, "=" * 28, "=" * 5
         )
-        for pn in self.model.params_ordered:
+        for pn in self.model.params:
             if nodmx and pn.startswith("DMX"):
                 continue
             prefitpar = getattr(self.model_init, pn)

src/pint/models/__init__.py

@@ -37,6 +37,7 @@
 from pint.models.solar_system_shapiro import SolarSystemShapiro
 from pint.models.solar_wind_dispersion import SolarWindDispersion, SolarWindDispersionX
 from pint.models.spindown import Spindown
+from pint.models.fdjump import FDJump
 
 # Import the main timing model classes
 from pint.models.timing_model import DEFAULT_ORDER, TimingModel

src/pint/models/binary_ell1.py

@@ -312,10 +312,12 @@ class BinaryELL1H(BinaryELL1):
     .. paramtable::
         :class: pint.models.binary_ell1.BinaryELL1H
 
-    Note
-    ----
+    Notes
+    -----
     Only the Medium-inclination case model is implemented.
 
+    Default value in `pint` for `NHARMS` is 7, while in `tempo2` it is 4.
+
     References
     ----------
     - Freire & Wex (2010), MNRAS, 409 (1), 199-212 [1]_

src/pint/models/fdjump.py

@@ -0,0 +1,205 @@
+"""System and frequency dependent delays to model profile evolution."""
+
+import re
+from warnings import warn
+
+import astropy.units as u
+import numpy as np
+
+from pint.models.parameter import boolParameter, maskParameter
+from pint.models.timing_model import DelayComponent
+
+fdjump_max_index = 20
+
+
+class FDJump(DelayComponent):
+    """A timing model for system-dependent frequency evolution of pulsar
+    profiles.
+
+    This model expresses the delay as a polynomial function of the
+    observing frequency/logarithm of observing frequency in the SSB frame.
+    This is intended to compensate for the delays introduced by frequency-dependent
+    profile structure when a different profiles are used for different systems.
+
+    The default behavior is to have FDJUMPs as polynomials of the observing
+    frequency (rather than log-frequency). This is different from the convention
+    used for global FD parameters. This choice is made to be compatible with tempo2.
+    This is controlled using the FDJUMPLOG parameter. "FDJUMPLOG Y" may not be
+    tempo2-compatible.
+
+    Note
+    ----
+    FDJUMPs have two indices: the polynomial/FD/prefix index and the system/mask
+    index. i.e., they have properties of both maskParameters such as JUMPs and
+    prefixParameters such as FDs. There is currently no elegant way in PINT to implement
+    such parameters due to the way parameter indexing is implemented; there is no way to
+    distinguish between mask and prefix indices.
+
+    Hence, they are implemented here as maskParameters as a stopgap measure.
+    This means that there must be an upper limit for the FD indices. This is controlled
+    using the `pint.models.fdjump.fdjump_max_index` variable, and is 20 by default.
+    Note that this is strictly a limitation of the implementation and not a property
+    of FDJUMPs themselves.
+
+    FDJUMPs appear in tempo2-format par files as "FDJUMPp", where p is the FD index.
+    The mask index is not explicitly mentioned in par files similar to JUMPs.
+    PINT understands both "FDJUMPp" and "FDpJUMP" as the same parameter in par files,
+    but the internal representation is always "FDpJUMPq", where q is the mask index.
+
+    PINT understands 'q' as the mask parameter just fine, but the identification of 'p'
+    as the prefix parameter is done in a hacky way.
+
+    This implementation may be overhauled in the future.
+
+    Parameters supported:
+
+    .. paramtable::
+        :class: pint.models.fdjump.FDJump
+    """
+
+    register = True
+    category = "fdjump"
+
+    def __init__(self):
+        super().__init__()
+
+        # Matches "FDpJUMPq" where p and q are integers.
+        self.param_regex = re.compile("^FD(\\d+)JUMP(\\d+)")
+
+        self.add_param(
+            boolParameter(
+                name="FDJUMPLOG",
+                value=False,
+                description="Whether to use log-frequency (Y) or linear-frequency (N) for computing FDJUMPs.",
+            )
+        )
+        for j in range(1, fdjump_max_index + 1):
+            self.add_param(
+                maskParameter(
+                    name=f"FD{j}JUMP",
+                    units="second",
+                    description=f"System-dependent FD parameter of polynomial index {j}",
+                )
+            )
+
+        self.delay_funcs_component += [self.fdjump_delay]
+
+    def setup(self):
+        super().setup()
+
+        self.fdjumps = [
+            mask_par
+            for mask_par in self.get_params_of_type("maskParameter")
+            if self.param_regex.match(mask_par)
+        ]
+
+        for fdj in self.fdjumps:
+            # prevents duplicates from being added to phase_deriv_funcs
+            if fdj in self.deriv_funcs.keys():
+                del self.deriv_funcs[fdj]
+            self.register_deriv_funcs(self.d_delay_d_FDJUMP, fdj)
+
+    def get_fd_index(self, par):
+        """Extract the FD index from an FDJUMP parameter name. In a parameter name
+        "FDpJUMPq", p is the FD/prefix index and q is the mask index.
+
+        Parameters
+        ----------
+        par: Parameter name (str)
+
+        Returns
+        -------
+        FD index (int)
+        """
+        if m := self.param_regex.match(par):
+            return int(m.groups()[0])
+        else:
+            raise ValueError(
+                f"The given parameter {par} does not correspond to an FDJUMP."
+            )
+
+    def get_freq_y(self, toas):
+        """Get frequency or log-frequency in GHz based on the FDJUMPLOG value.
+        Returns (freq/1_GHz) if FDJUMPLOG==N and log(freq/1_GHz) if FDJUMPLOG==Y.
+        Any non-finite values are replaced by zero.
+
+        Parameters
+        ----------
+        toas: pint.toa.TOAs
+
+        Returns
+        -------
+        (freq/1_GHz) or log(freq/1_GHz) depending on the value of FDJUMPLOG (float).
+        """
+        tbl = toas.table
+        try:
+            freq = self._parent.barycentric_radio_freq(toas)
+        except AttributeError:
+            warn("Using topocentric frequency for frequency dependent delay!")
+            freq = tbl["freq"]
+
+        y = (
+            np.log(freq.to(u.GHz).value)
+            if self.FDJUMPLOG.value
+            else freq.to(u.GHz).value
+        )
+        non_finite = np.invert(np.isfinite(y))
+        y[non_finite] = 0.0
+
+        return y
+
+    def fdjump_delay(self, toas, acc_delay=None):
+        """Calculate frequency dependent delay.
+
+        If FDJUMPLOG is Y, use the following expression (similar to global FD parameters):
+
+            FDJUMP_delay = sum_i(c_i * (log(obs_freq/1GHz))^i)
+
+        If FDJUMPLOG is N, use the following expression (same as in tempo2, default):
+
+            FDJUMP_delay = sum_i(c_i * (obs_freq/1GHz)^i)
+        """
+        y = self.get_freq_y(toas)
+
+        delay = np.zeros_like(y)
+        for fdjump in self.fdjumps:
+            fdj = getattr(self, fdjump)
+            if fdj.quantity is not None:
+                mask = fdj.select_toa_mask(toas)
+                ymask = y[mask]
+                fdidx = self.get_fd_index(fdjump)
+                fdcoeff = fdj.value
+                delay[mask] += fdcoeff * ymask**fdidx
+
+        return delay * u.s
+
+    def d_delay_d_FDJUMP(self, toas, param, acc_delay=None):
+        """Derivative of delay w.r.t. FDJUMP parameters."""
+        assert (
+            bool(self.param_regex.match(param))
+            and hasattr(self, param)
+            and getattr(self, param).quantity is not None
+        ), f"{param} is not present in the FDJUMP model."
+
+        y = self.get_freq_y(toas)
+        mask = getattr(self, param).select_toa_mask(toas)
+        ymask = y[mask]
+        fdidx = self.get_fd_index(param)
+
+        delay_derivative = np.zeros_like(y)
+        delay_derivative[mask] = ymask**fdidx
+
+        return delay_derivative * u.dimensionless_unscaled
+
+    def print_par(self, format="pint"):
+        par = super().print_par(format)
+
+        if format != "tempo2":
+            return par
+
+        for fdjump in self.fdjumps:
+            if getattr(self, fdjump).quantity is not None:
+                j = self.get_fd_index(fdjump)
+                par = par.replace(f"FD{j}JUMP", f"FDJUMP{j}")
+
+        return par