Merge pull request #1557 from abhisrkckl/phase-offset

Explicit phase offset as a fittable parameter

CHANGELOG.md

@@ -13,6 +13,10 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 - `validate_component_types` method for more rigorous validation of timing model components.
 - roundtrip test to make sure clock corrections are not written to tim files
 - `calc_phase_mean` and `calc_time_mean` methods in `Residuals` class to compute the residual mean.
+- `PhaseOffset` component (overall phase offset between physical and TZR toas)
+- `tzr` attribute in `TOAs` class to identify TZR TOAs
+- Documentation: Explanation for offsets
+- Example: `phase_offset_example.py`
 ### Fixed
 - fixed docstring for `add_param_from_top`
 - Gridded calculations now respect logger settings

docs/examples/fit_NGC6440E.py

@@ -61,7 +61,7 @@
 # ```python
 # # Define the timing model
 # m = get_model(parfile)
-# # Read in the TOAs, using the solar system epemeris and other things from the model
+# # Read in the TOAs, using the solar system ephemeris and other things from the model
 # t = pint.toa.get_TOAs(timfile, model=m)
 # ```
 

docs/examples/phase_offset_example.py

@@ -0,0 +1,139 @@
+#! /usr/bin/env python
+# ---
+# jupyter:
+#   jupytext:
+#     cell_metadata_filter: -all
+#     formats: ipynb,py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.14.4
+#   kernelspec:
+#     display_name: .env
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Demonstrate phase offset
+#
+# This notebook is primarily designed to operate as a plain `.py` script.
+# You should be able to run the `.py` script that occurs in the
+# `docs/examples/` directory in order to carry out a simple fit of a
+# timing model to some data. You should also be able to run the notebook
+# version as it is here (it may be necessary to `make notebooks` to
+# produce a `.ipynb` version using `jupytext`).
+
+# %%
+from pint.models import get_model_and_toas, PhaseOffset
+from pint.residuals import Residuals
+from pint.config import examplefile
+from pint.fitter import DownhillWLSFitter
+
+import matplotlib.pyplot as plt
+from astropy.visualization import quantity_support
+
+quantity_support()
+
+# %%
+# Read the TOAs and the model
+parfile = examplefile("J1028-5819-example.par")
+timfile = examplefile("J1028-5819-example.tim")
+model, toas = get_model_and_toas(parfile, timfile)
+
+# %%
+# Create a Residuals object
+res = Residuals(toas, model)
+
+# %%
+# By default, the residuals are mean-subtracted.
+resids1 = res.calc_phase_resids().to("")
+
+# We can disable mean subtraction by setting `subtract_mean` to False.
+resids2 = res.calc_phase_resids(subtract_mean=False).to("")
+
+# %%
+# Let us plot the residuals with and without mean subtraction.
+# In the bottom plot, there is clearly an offset between the two cases
+# although it is not so clear in the top plot.
+
+mjds = toas.get_mjds()
+errors = toas.get_errors() * model.F0.quantity
+
+plt.subplot(211)
+plt.errorbar(mjds, resids1, errors, ls="", marker="x", label="Mean subtracted")
+plt.errorbar(mjds, resids2, errors, ls="", marker="x", label="Not mean subtracted")
+plt.xlabel("MJD")
+plt.ylabel("Phase residuals")
+plt.axhline(0, ls="dotted", color="grey")
+plt.legend()
+
+plt.subplot(212)
+plt.plot(mjds, resids2 - resids1, ls="", marker="x")
+plt.xlabel("MJD")
+plt.ylabel("Phase residual difference")
+plt.show()
+
+# %%
+# This phase offset that gets subtracted implicitly can be computed
+# using the `calc_phase_mean` function. There is also a similar function
+# `calc_time_mean` for time offsets.
+
+implicit_offset = res.calc_phase_mean().to("")
+print("Implicit offset = ", implicit_offset)
+
+# %%
+# Now let us look at the design matrix.
+T, Tparams, Tunits = model.designmatrix(toas)
+print("Design matrix params :", Tparams)
+
+# The implicit offset is represented as "Offset".
+
+# %%
+# We can explicitly fit for this offset using the "PHOFF" parameter.
+# This is available in the "PhaseOffset" component
+
+po = PhaseOffset()
+model.add_component(po)
+assert hasattr(model, "PHOFF")
+
+# %%
+# Let us fit this now.
+
+model.PHOFF.frozen = False
+ftr = DownhillWLSFitter(toas, model)
+ftr.fit_toas()
+print(
+    f"PHOFF fit value = {ftr.model.PHOFF.value} +/- {ftr.model.PHOFF.uncertainty_value}"
+)
+
+# This is consistent with the implicit offset we got earlier.
+
+# %%
+# Let us plot the post-fit residuals.
+
+mjds = ftr.toas.get_mjds()
+errors = ftr.toas.get_errors() * model.F0.quantity
+resids = ftr.resids.calc_phase_resids().to("")
+
+plt.errorbar(mjds, resids, errors, ls="", marker="x", label="After fitting PHOFF")
+plt.xlabel("MJD")
+plt.ylabel("Phase residuals")
+plt.axhline(0, ls="dotted", color="grey")
+plt.legend()
+plt.show()
+
+# %%
+# Let us compute the phase residual mean again.
+phase_mean = ftr.resids.calc_phase_mean().to("")
+print("Phase residual mean = ", phase_mean)
+
+# i.e., we have successfully gotten rid of the offset by fitting PHOFF.
+
+# %%
+# Now let us look at the design matrix again.
+T, Tparams, Tunits = model.designmatrix(toas)
+print("Design matrix params :", Tparams)
+
+# The explicit offset "PHOFF" has replaced the implicit "Offset".

docs/explanation.rst

@@ -16,6 +16,8 @@ you will find PINT sufficient for all your needs!
 .. _TEMPO: http://tempo.sourceforge.net/
 .. _TEMPO2: https://www.atnf.csiro.au/research/pulsar/tempo2/
 
+.. sectnum::
+
 Time
 ----
 
@@ -215,12 +217,49 @@ The total DM and dispersion slope predicted by a given timing model (:class:`pin
 for a given set of TOAs (:class:`pint.toa.TOAs`) can be computed using :func:`pint.models.timing_model.TimingModel.total_dm`
 and :func:`pint.models.timing_model.TimingModel.total_dispersion_slope` methods respectively.
 
+Offsets in pulsar timing
+------------------------
+Offsets arise in pulsar timing models for a variety of reasons. The different types of
+offsets are listed below:
+
+Overall phase offset
+''''''''''''''''''''
+The pulse phase corresponding to the TOAs are usually computed in reference to an arbitrary 
+fiducial TOA known as the TZR TOA (see :class:`pint.models.absolute_phase.AbsPhase`). Since the 
+choice of the TZR TOA is arbitrary, there can be an overall phase offset between the TZR TOA and 
+the measured TOAs. There are three ways to account for this offset: (1) subtract the weighted mean
+from the timing residuals, (2) make the TZR TOA (given by the `TZRMJD` parameter) fittable, or 
+(3) introduce a fittable phase offset parameter between measured TOAs and the TZR TOA.
+Traditionally, pulsar timing packages have opted to implicitly subtract the residual mean, and this
+is the default behavior of `PINT`. Option (2) is hard to implement because the TZR TOA may be 
+specified at any observatory, and computing the TZR phase requires the application of the clock 
+corrections. The explicit phase offset (option 3) can be invoked by adding the `PHOFF` parameter, 
+(implemented in :class:`pint.models.phase_offset.PhaseOffset`). If the explicit offset `PHOFF`
+is given, the implicit residual mean subtraction behavior will be disabled.
+
+System-dependent delays
+'''''''''''''''''''''''
+It is very common to have TOAs for the same pulsar obtained using different observatories, 
+telescope receivers, backend systems, and data processing pipelines, especially in long-running 
+campaigns. Delays can arise between the TOAs measured using such different systems due to, among
+other reasons, instrumental delays, differences in algorithms used for RFI mitigation, folding, TOA 
+measurement etc., and the choice of different template profiles used for TOA measurement. Such 
+offsets are usually modeled using phase jumps (the `JUMP` parameter, see :class:`pint.models.jump.PhaseJump`) 
+between TOAs generated from different systems.
+
+System-dependent DM offsets
+'''''''''''''''''''''''''''
+Similar to system-dependent delays, offsets can arise between wideband DM values measured using 
+different systems due to the choice of template portraits with different fiducial DMs. This is 
+usually modeled using DM jumps (the `DMJUMP` parameter, see :class:`pint.models.dispersion_model.DispersionJump`).
+
 Observatories
 -------------
 
-PINT comes with a number of defined observatories.  Those on the surface of the Earth are :class:`~pint.observatory.topo_obs.TopoObs` 
-instances.  It can also pull in observatories from ``astropy``, 
-and you can define your own.  Observatories are generally referenced when reading TOA files, but can also be accessed directly::
+PINT comes with a number of defined observatories.  Those on the surface of the Earth 
+are :class:`~pint.observatory.topo_obs.TopoObs` instances.  It can also pull in observatories 
+from ``astropy``, and you can define your own.  Observatories are generally referenced when 
+reading TOA files, but can also be accessed directly::
 
   import pint.observatory
   gbt = pint.observatory.get_observatory("gbt")  
@@ -246,12 +285,14 @@ The observatory data are stored in JSON format.  A simple example is::
         "origin": "The Robert C. Byrd Green Bank Telescope.\nThis data was obtained by Joe Swiggum from Ryan Lynch in 2021 September.\n"
     }
 
-The observatory is defined by its name (``gbt``) and its position.  This can be given as geocentric coordinates in the 
-International_Terrestrial_Reference_System_ (ITRF) through the ``itrf_xyz`` triple (units as ``m``), or geodetic coordinates 
-(WGS84_ assumed) through ``lat``, ``lon``, ``alt`` 
-(units are ``deg`` and ``m``).  Conversion is done through Astropy_EarthLocation_.
+The observatory is defined by its name (``gbt``) and its position.  This can be given as 
+geocentric coordinates in the International_Terrestrial_Reference_System_ (ITRF) through 
+the ``itrf_xyz`` triple (units as ``m``), or geodetic coordinates (WGS84_ assumed) through 
+``lat``, ``lon``, ``alt`` (units are ``deg`` and ``m``).  Conversion is done through 
+Astropy_EarthLocation_.
 
-Other attributes are optional.  Here we have also specified the ``tempo_code`` and ``itoa_code``, and a human-readable ``origin`` string.
+Other attributes are optional.  Here we have also specified the ``tempo_code`` and 
+``itoa_code``, and a human-readable ``origin`` string.
 
 A more complex/complete example is::
 
@@ -282,19 +323,22 @@ A more complex/complete example is::
         ]
     }
 
-Here we have included additional explicit ``aliases``, specified the clock format via ``clock_fmt``, and specified that the last entry in the 
-clock file is bogus (``bogus_last_correction``).  There are two clock files included in ``clock_file``:
+Here we have included additional explicit ``aliases``, specified the clock format via 
+``clock_fmt``, and specified that the last entry in the clock file is bogus (``bogus_last_correction``).  
+There are two clock files included in ``clock_file``:
 
 * ``jbroach2jb.clk`` (where we also specify that it is ``valid_beyond_ends``)
 * ``jb2gps.clk``
 
-These are combined to reference this particular telescope/instrument combination.  For the full set of options, see :class:`~pint.observatory.topo_obs.TopoObs`.
+These are combined to reference this particular telescope/instrument combination.  
+For the full set of options, see :class:`~pint.observatory.topo_obs.TopoObs`.
 
 
 Adding New Observatories
 ''''''''''''''''''''''''
 
-In addition to modifying ``pint.config.runtimefile("observatories.json")``, there are other ways to add new observatories.  
+In addition to modifying ``pint.config.runtimefile("observatories.json")``, there are other 
+ways to add new observatories.  
 **Make sure you define any new observatory before you load any TOAs.**
 
 1. You can define them pythonically:
@@ -304,7 +348,8 @@ In addition to modifying ``pint.config.runtimefile("observatories.json")``, ther
     import astropy.coordinates
     newobs = pint.observatory.topo_obs.TopoObs("newobs", location=astropy.coordinates.EarthLocation.of_site("keck"), origin="another way to get Keck")
 
-This can be done by specifying the ITRF coordinates, (``lat``, ``lon``, ``alt``), or a :class:`~astropy.coordinates.EarthLocation` instance.
+This can be done by specifying the ITRF coordinates, (``lat``, ``lon``, ``alt``), or a 
+:class:`~astropy.coordinates.EarthLocation` instance.
 
 2. You can include them just for the duration of your python session:
 ::
@@ -327,22 +372,26 @@ This can be done by specifying the ITRF coordinates, (``lat``, ``lon``, ``alt``)
         }"""
     load_observatories(io.StringIO(fakeGBT), overwrite=True)
 
-Note that since we are overwriting an existing observatory (rather than defining a completely new one) we specify ``overwrite=True``.  
+Note that since we are overwriting an existing observatory (rather than defining a 
+completely new one) we specify ``overwrite=True``.  
 
-3. You can define them in a different file on disk.  If you took the JSON above and put it into a file ``/home/user/anothergbt.json``, 
+3. You can define them in a different file on disk.  If you took the JSON above and 
+put it into a file ``/home/user/anothergbt.json``, 
 you could then do::
 
     export $PINT_OBS_OVERRIDE=/home/user/anothergbt.json
 
-(or the equivalent in your shell of choice) before you start any PINT scripts.  By default this will overwrite any existing definitions.
+(or the equivalent in your shell of choice) before you start any PINT scripts.  
+By default this will overwrite any existing definitions.
 
 4. You can rely on ``astropy``.  For instance:
 ::
 
     import pint.observatory
     keck = pint.observatory.Observatory.get("keck")
 
-will find Keck.  :func:`astropy.coordinates.EarthLocation.get_site_names` will return a list of potential observatories.
+will find Keck.  :func:`astropy.coordinates.EarthLocation.get_site_names` will return a list 
+of potential observatories.
 
 .. _International_Terrestrial_Reference_System: https://en.wikipedia.org/wiki/International_Terrestrial_Reference_System_and_Frame
 .. _WGS84: https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84

src/pint/data/examples/J1028-5819-example.par

@@ -0,0 +1,25 @@
+# Created: 2023-05-17T14:02:21.738913
+# PINT_version: 0.9.5+145.ga1930cdf.dirty
+# User: Abhimanyu Susobhanan (abhimanyu)
+# Host: abhimanyu-VirtualBox
+# OS: Linux-5.19.0-41-generic-x86_64-with-glibc2.35
+# Format: pint
+# Converted from tempo2 example file "example1.par"
+PSRJ                           J1028-5819
+UNITS                                 TDB
+DILATEFREQ                              N
+DMDATA                                  N
+NTOA                                    0
+CHI2                                  0.0
+RAJ                     10:28:28.00000000 1 0.00000000000000000000
+DECJ                   -58:19:05.20000000 1 0.00000000000000000000
+PMRA                                  0.0
+PMDEC                                 0.0
+PX                                    0.0
+POSEPOCH           54561.9998229616586976
+F0                  10.940532469635118635 1 0.0
+F1              -1.9300000598500644258e-12 1 0.0
+PEPOCH             54561.9998229616586976
+PLANET_SHAPIRO                          N
+DM                  96.525001496639228994
+DMEPOCH            54561.9998229616586976