This Jupyter notebook can be downloaded from paper_validation_example.ipynb, or viewed as a python script at paper_validation_example.py.

Validation Example for PINT paper

A comparison between PINT result and Tempo/Tempo2 result. This example is presented in the PINT paper to validate that PINT is able to process the PSR J1600-3053 NANOGrav 11-year data set, which uses DD binary model, and get a comparable result with TEMPO/TEMPO2. For more discussion see: https://arxiv.org/abs/2012.00074

This comparison includes the following steps: * PINT run on PSR J1600-3053 NANOGrav 11-year data set * Create PINT pre-fit residuals. * Re-fit data using PINT. * Create PINT post-fit residuals. * Obtain the PINT post-fit parameter values. * TEMPO run on PSR J1600-3053 NANOGrav 11-year data set * Create TEMPO pre-fit residuals. * Re-fit data using TEMPO. * Create TEMPO post-fit residuals. * Obtain the TEMPO post-fit parameter values. * Compare the pre-fit and post-fit residuals between PINT and TEMPO. * Compare the pre-fit and post-fit parameters between PINT and TEMPO.
* TEMPO2 run on PSR J1600-3053 NANOGrav 11-year data set * Create TEMPO2 pre-fit residuals. * Re-fit data using TEMPO2. * Create TEMPO2 post-fit residuals. * Obtain the TEMPO2 post-fit parameter values. * Compare the pre-fit and post-fit residuals between PINT and TEMPO2. * Compare the pre-fit and post-fit parameters between PINT and TEMPO2.

Method of this comparison:

The method we used in this comparison starts from the published data set with a converged timing model. The data in this comparison are produced by TEMPO. 0. The PINT pre-fit residuals and post-fit residuals should be close and no visible difference. If the post-fit and pre-fit residuals are very different, it means PINT has a bug or some different definitions of parameters(e.g., PINT DDK model has different definition of KOM parameter).
1. Compare the pre-fit residuals between PINT-TEMPO/TEMPO2, to see if there are outstanding discrepancies(residual difference > 10 ns, since we are aim to have our machine precision in the 1 ns level). These discrepancies are likely caused by different code implements and we try to verify the reasons(some of them are bugs in PINT or in TEMPO/TEMPO2, but some of them are just implementation choices). 2. Fit the data and compare the chi2 value and post-fit residuals between PINT and TEMPO/TEMPO2 to exam if the fitter returns a reasonable result(i.e., chi2 and post-fit rms do not deviate from TEMPO/TEMPO2’s value too much, etc.). If the fit gives a “non-sense” result, that is worth to do more investigation. One of the useful indicators is the discrepancy of the post-fit residuals between PINT and TEMPO/TEMPO2. 3. Compare the post-fit parameters between PINT and TEMPO/TEMPO2. Since the initial parameter are from the same data, there should be not difference. The PINT post-fit parameters should change within the original uncertainties(TEMPO uncertainties), and the PINT uncertainty vs TEMPO/TEMPO2 uncertainty should be close.

Requirement: * Data set: PRS J1600-3053 NANOGrav 11-year data * One copy of PRS J1600-3053 NANOGrav 11-year data is included in the PINT source code docs/examples/J1600-3053_NANOGrav_11yv1.gls.par and docs/examples/J1600-3053_NANOGrav_11yv1.tim, which is the default data path in this notebook. Note, this requires the user to download the PINT source code from github. * The offical NANOGrav 11-year data can be downloaded at: https://data.nanograv.org/. The data path should be changed to the data location. * PINT version: 0.8.0 or higher. * TEMPO and its python utils tempo_utils. * TEMPO version for current comparison: 13.101 (2020-11-04 c5fbddf) * TEMPO2 and its python utils tempo2_utils. * TEMPO2 version for current comparison: 2019.01.1 * TEMPO_utils and TEMPO2_utils are packaged together and can be downloaded from https://github.com/demorest/tempo_utils. * TEMPO2 general2 plugins.

[1]:

import pint
import sys
from pint import toa
from pint import models
from pint.fitter import GLSFitter
import os
import matplotlib.pyplot as plt
import astropy.units as u
import tempo2_utils as t2u
import tempo_utils
import tempo2_utils
import numpy as np
from astropy.table import Table
from astropy.io import ascii
import subprocess
import tempfile
from pint import ls
import astropy.constants as ct
from pint.solar_system_ephemerides import objPosVel_wrt_SSB
from astropy.time import Time

Redefine the Tempo2_util function for larger number of observations

[3]:

_nobs = 30000
def newpar2(parfile,timfile):
    """
    Run tempo2, return new parfile (as list of lines).  input parfile
    can be either lines or a filename.
    """
    orig_dir = os.getcwd()
    try:
        temp_dir = tempfile.mkdtemp(prefix="tempo2")
        try:
            lines = open(parfile,'r').readlines()
        except:
            lines = parfile
        open("%s/pulsar.par" % temp_dir, 'w').writelines(lines)
        timpath = os.path.abspath(timfile)
        os.chdir(temp_dir)
        cmd = "tempo2 -nobs %d -newpar -f pulsar.par %s -norescale" % (_nobs, timpath)
        os.system(cmd + " > /dev/null")
        outparlines = open('new.par').readlines()
    finally:
        os.chdir(orig_dir)
    os.system("rm -rf %s" % temp_dir)
    for l in outparlines:
        if l.startswith('TRES'): rms = float(l.split()[1])
        elif l.startswith('CHI2R'): (foo, chi2r, ndof) = l.split()
    return float(chi2r)*float(ndof), int(ndof), rms, outparlines

Set up date file path for PSR J1600-3053.

  • Note

    • This path only works when PINT is installed from source code, which has docs and example directories.

[4]:

psr = "J1600-3053"
par_file = os.path.join('../examples', psr + "_NANOGrav_11yv1.gls.par")
tim_file = os.path.join('../examples', psr + "_NANOGrav_11yv1.tim")

PINT run

Load TOAs to PINT

[5]:

t = toa.get_TOAs(tim_file, ephem="DE436", bipm_version="BIPM2015")

INFO: Applying clock corrections (include_gps = True, include_bipm = True) [pint.toa]
INFO: Observatory gbt, loading clock file
        /home/luo/.local/lib/python3.6/site-packages/pint/datafiles/time.dat [pint.observatory.topo_obs]
INFO: Applying observatory clock corrections. [pint.observatory.topo_obs]
INFO: Applying GPS to UTC clock correction (~few nanoseconds) [pint.observatory.topo_obs]
INFO: Observatory gbt, loading GPS clock file
        /home/luo/.local/lib/python3.6/site-packages/pint/datafiles/gps2utc.clk [pint.observatory.topo_obs]
INFO: Applying TT(TAI) to TT(BIPM2015) clock correction (~27 us) [pint.observatory.topo_obs]
INFO: Observatory gbt, loading BIPM clock file
        /home/luo/.local/lib/python3.6/site-packages/pint/datafiles/tai2tt_bipm2015.clk [pint.observatory.topo_obs]
INFO: Computing TDB columns. [pint.toa]
INFO: Using EPHEM = DE436 for TDB calculation. [pint.toa]
INFO: Computing PosVels of observatories and Earth, using DE436 [pint.toa]
INFO: Set solar system ephemeris to link:
        https://data.nanograv.org/static/data/ephem/de436.bsp [pint.solar_system_ephemerides]

[6]:

print("There are {} TOAs in the dataset.".format(t.ntoas))

There are 12433 TOAs in the dataset.

Load timing model from .par file

Since PINT only uses the IAU 2000 version of precession-nutation model but NANOGrav 11-year data uses old precession-nutation model, You will see a UserWarning: PINT only supports 'T2CMETHOD IAU2000B'.

[7]:

m = models.get_model(par_file)

INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]
INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]

/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")

Make the General Least Square fitter

[8]:

f = GLSFitter(model=m, toas=t)

Fit TOAs for 9 iterations.

The expected chi2 value should be close to TEMPO and TEMPO2, but not the same.

[9]:

chi2 = f.fit_toas(9)
print("Postfit Chi2: ", chi2)
print("Degree of freedom: ", f.resids.dof)

/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")
/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")

Postfit Chi2:  12368.094375740437515
Degree of freedom:  12307

The weighted RMS value for pre-fit and post-fit residuals

[10]:

print("Pre-fit residual weighted RMS:", f.resids_init.rms_weighted())
print("Post-fit residual weighted RMS:", f.resids.rms_weighted())

Pre-fit residual weighted RMS: 0.944170684867224 us
Post-fit residual weighted RMS: 0.9441138383219785 us

Plot the pre-fit and post-fit residuals

[11]:

pint_prefit = f.resids_init.time_resids.to_value(u.us)
pint_postfit = f.resids.time_resids.to_value(u.us)

plt.figure(figsize=(8,5), dpi=150)
plt.subplot(2, 1, 1)
plt.errorbar(t.get_mjds().to_value(u.day), f.resids_init.time_resids.to_value(u.us),
             yerr=t.get_errors().to_value(u.us), fmt='x')

plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (us)')
plt.title('PINT pre-fit residuals for PSR J1600-3053 NANOGrav 11-year data')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.errorbar(t.get_mjds().to_value(u.day), f.resids.time_resids.to_value(u.us),
             yerr=t.get_errors().to_value(u.us), fmt='x')
plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (us)')
plt.title('PINT post-fit residuals for PSR J1600-3053 NANOGrav 11-year data')
plt.grid(True)
plt.tight_layout()
plt.savefig("J1600_PINT")
../_images/examples-rendered_paper_validation_example_20_0.png

TEMPO run

Use tempo_utils to analysis the same data set.

[12]:

tempo_toa = tempo_utils.read_toa_file(tim_file)
tempo_chi2, ndof, rms_t, tempo_par = tempo_utils.run_tempo(tempo_toa ,par_file, get_output_par=True,
                                                   gls=True)

[13]:

print("TEMPO postfit chi2: ", tempo_chi2)
print("TEMPO postfit weighted rms: ", rms_t)

TEMPO postfit chi2:  12368.46
TEMPO postfit weighted rms:  0.944

Write the TEMPO postfit model to a new .par file, for comparison later

[14]:

# Write out the post fit tempo parfile.
tempo_parfile = open(psr + '_tempo.par', 'w')
for line in tempo_par:
    tempo_parfile.write(line)
tempo_parfile.close()

Get the TEMPO residuals

[15]:

tempo_prefit = tempo_toa.get_prefit()
tempo_postfit = tempo_toa.get_resids()
mjds = tempo_toa.get_mjd()
freqs = tempo_toa.get_freq()
errs = tempo_toa.get_resid_err()

Plot the PINT - TEMPO residual difference.

[16]:

tp_diff_pre = (pint_prefit - tempo_prefit) * u.us
tp_diff_post = (pint_postfit - tempo_postfit) * u.us

[17]:

plt.figure(figsize=(8,5), dpi=150)
plt.subplot(2, 1, 1)
plt.plot(mjds, (tp_diff_pre - tp_diff_pre.mean()).to_value(u.ns), '+')
plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (ns)')
plt.title('PSR J1600-3053 prefit residual differences between PINT and TEMPO')
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(mjds, (tp_diff_post - tp_diff_post.mean()).to_value(u.ns), '+')
plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (ns)')
plt.title('PSR J1600-3053 postfit residual differences between PINT and TEMPO')
plt.grid(True)
plt.tight_layout()
plt.savefig("J1600_PINT_tempo.eps")
../_images/examples-rendered_paper_validation_example_30_0.png

The PINT-TEMPO pre-fit residual discrepancy is due to the different precession-nutation model in the two packages. * TEMPO: IAU6501976 precession and IAU 1980 nutation. * PINT: IAU2000B precession-nutation.

Compare the parameter between TEMPO and PINT

  • Reported quantities

    • TEMPO value

    • TEMPO uncertainty

    • Parameter units

    • TEMPO parameter value - PINT parameter value

    • TEMPO/PINT parameter absolute difference divided by TEMPO uncertainty

    • PINT uncertainty divided by TEMPO uncertainty, if TEMPO provides the uncertainty value

[18]:

# Create the parameter compare table
tv = []
tu = []
tv_pv = []
tv_pv_tc = []
tc_pc = []
units = []
names = []
no_t_unc = []
tempo_new_model = models.get_model(psr + '_tempo.par')
for param in tempo_new_model.params:
    t_par = getattr(tempo_new_model, param)
    pint_par = getattr(f.model, param)
    tempoq = t_par.quantity
    pintq = pint_par.quantity
    try:
        diffq =  tempoq - pintq
        if t_par.uncertainty_value != 0.0:
            diff_tcq = np.abs(diffq) / t_par.uncertainty
            uvsu = pint_par.uncertainty / t_par.uncertainty
            no_t_unc.append(False)
        else:
            diff_tcq = np.abs(diffq) / pint_par.uncertainty
            uvsu = t_par.uncertainty
            no_t_unc.append(True)
    except TypeError:
        continue
    uvsu = pint_par.uncertainty / t_par.uncertainty
    tv.append(tempoq.value)
    tu.append(t_par.uncertainty.value)
    tv_pv.append(diffq.value)
    tv_pv_tc.append(diff_tcq.value)
    tc_pc.append(uvsu)
    units.append(t_par.units)
    names.append(param)

compare_table = Table((names, tv, tu, units, tv_pv, tv_pv_tc, tc_pc, no_t_unc), names = ('name', 'Tempo Value', 'Tempo uncertainty', 'units',
                                                                                      'Tempo_V-PINT_V',
                                                                                      'Tempo_PINT_diff/unct',
                                                                                      'PINT_unct/Tempo_unct',
                                                                                      'no_t_unc'))
compare_table.sort('Tempo_PINT_diff/unct')
compare_table = compare_table[::-1]
compare_table.write('parameter_compare.t.html', format='html', overwrite=True)

INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]
INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]

/home/luo/.local/lib/python3.6/site-packages/pint/models/timing_model.py:304: UserWarning: PINT only supports 'T2CMETHOD IAU2000B'
  warn("PINT only supports 'T2CMETHOD IAU2000B'")

If one wants the Latex output please use the cell below.

[20]:

#ascii.write(compare_table, sys.stdout, Writer = ascii.Latex,
#            latexdict = {'tabletype': 'table*'})

Check out the maximum DMX difference

[21]:

max_dmx = 0
max_dmx_index = 0
for ii, row in enumerate(compare_table):
    if row['name'].startswith('DMX_'):
        if row['Tempo_PINT_diff/unct'] > max_dmx:
            max_dmx = row['Tempo_PINT_diff/unct']
            max_dmx_index = ii

dmx_max = compare_table[max_dmx_index]['name']

compare_table[max_dmx_index]

[21]:
Row index=6
nameTempo ValueTempo uncertaintyunitsTempo_V-PINT_VTempo_PINT_diff/unctPINT_unct/Tempo_unctno_t_unc
str8str48float128objectfloat128float128float128bool
DMX_00100.000669275610.00020051850499999999489pc / cm3-5.0888210248638612865e-060.0253783112180287866170.9999978605214362437False

Output the table in the paper

[22]:

paper_params = ['F0', 'F1', 'FD1', 'FD2', 'JUMP1', 'PX',
                'ELONG', 'ELAT', 'PMELONG', 'PMELAT', 'PB',
                'A1', 'A1DOT', 'ECC', 'T0', 'OM', 'OMDOT', 'M2',
                'SINI', dmx_max]
# Get the table index of the parameters above
paper_param_index = []
for pp in paper_params:
    # We assume the parameter name are unique in the table
    idx = np.where(compare_table['name'] == pp)[0][0]
    paper_param_index.append(idx)
paper_param_index = np.array(paper_param_index)
compare_table[paper_param_index]

[22]:
Table length=20
nameTempo ValueTempo uncertaintyunitsTempo_V-PINT_VTempo_PINT_diff/unctPINT_unct/Tempo_unctno_t_unc
str8str48float128objectfloat128float128float128bool
F0277.93771124297461485.186e-13Hz-1.471045507628332416e-140.0283657058933346011571.0000737060332136081False
F1-7.338737472765e-164.619148184227e-21Hz / s6.3620434143467849095e-230.0137731961838142522691.0001125817340140925False
FD13.98314325e-051.6566479199999999207e-06s-2.5460168459788480762e-090.00153684848496888118960.9999972210804803918False
FD2-1.47296057e-051.1922595999999999884e-06s1.3701818969060749554e-090.00114923117155531803560.99999858894148907495False
JUMP1-8.789e-061.2999999999999999941e-07s-4.6499257614788852208e-100.00357686597036837314671.0037094619649047367False
PX0.5040.07349999999999999589mas-0.00207140404968053637450.0281823680228644422860.99982583641940803165False
ELONG244.3476778440795.9573e-09deg-5.923646156924533557e-100.099435082284332391470.99997665908491340815False
ELAT-10.07183902536513.36103e-08deg-3.1908024911500576515e-090.094935257678451476231.0000723294528560605False
PMELONG0.46260.010399999999999999523mas / yr0.000711869555339084136850.068448995705681164871.0031591658256824307False
PMELAT-7.15550.058200000000000001732mas / yr-0.000504848952748382373560.0086743806314155038480.9992173822396303029False
PB14.348465725503022.1222661e-06d-3.4570101856666590745e-080.0162892400046660457631.0000725690593949082False
A18.8016531228.1100000000000004906e-07ls1.4913034362962207524e-080.0183884517422468647670.9844174919313783967False
A1DOT-4e-156.260000000000000155e-16ls / s8.912568233346704436e-180.0142373294462407422310.99986821703571171494False
ECC0.00017372948.9000000000000002855e-09-2.385671795204248602e-100.0268053010697106575131.0022775158955208319False
T055878.26189804510000000.0005167676d-1.0512794592798524462e-050.0203433701973547189520.9909423736001105048False
OM181.849568165780.01296546975deg-0.000263761398727566098720.0203433738856678215390.9909564544661490358False
OMDOT0.00522090.0013554deg / yr-2.2108721249702442383e-050.0163115842184612973171.0000993017977168461False
M20.2718940.089418999999999998485solMass-0.00164113817268091155550.0183533496536632222130.97866231479642251667False
SINI0.9062850.033993000000000002380.000543558683816058874070.0159903122353443026550.9838880245552520387False
DMX_00100.000669275610.00020051850499999999489pc / cm3-5.0888210248638612865e-060.0253783112180287866170.9999978605214362437False

TEMPO2 run

Before TEMPO2 run, the .par file has to be modified for a more accurate TEMPO2 vs PINT comparison. We save the modified .par file in a file named “[PSR name]_tempo2.par”. In this case, “J1600-3053_tempo2.par”

  • Modified parameters in the .par file:

    • ECL IERS2010 —-> ECL IERS 2003 (In this version of TEMPO2, the IERS 2003 Obliquity angle is hardcoded in its code. To match TEMPO2’s default value, we change the ECL to IERS 2003 in the .par file )

    • T2CMETHOD TEMPO —-> # T2CMETHOD TEMPO (TEMPO2 supports both IAU 2000 precession-nutation model and old TEMPO-style model. To make TEMPO2 ues its default precession and nutation model, IAU 2000, this line in the .par file has to be commmented out.)

  • Note, this modified .par file is provided in the docs/examples directory. If PINT is not installed from source code, one have to modify the .par file from the NANOGrav 11-year data.

[23]:

tempo2_par = os.path.join('../examples',"J1600-3053_tempo2.par")

PINT refit using the modified tempo2-style parfile

[24]:

m_t2 = models.get_model(tempo2_par)

INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]
INFO: Parameter A1DOT's value will be scaled by 1e-12 [pint.models.parameter]

[25]:

f_t2 = GLSFitter(toas=t, model=m_t2)
f_t2.fit_toas()

[25]:

12368.092265853045861

Tempo2 fit

[26]:

tempo2_chi2, ndof, rms_t2, tempo2_new_par = newpar2(tempo2_par, tim_file)
print("TEMPO2 chi2: ", tempo2_chi2)
print("TEMPO2 rms: ", rms_t2)

TEMPO2 chi2:  12265.156200000001
TEMPO2 rms:  0.944

Get TEMPO2 residuals, toa value, observing frequencies, and data error

[27]:

tempo2_result = t2u.general2(tempo2_par, tim_file, ['sat', 'pre', 'post', 'freq', 'err'])
# TEMPO2's residual unit is second
tp2_diff_pre = f_t2.resids_init.time_resids - tempo2_result['pre'] * u.s
tp2_diff_post = f_t2.resids.time_resids - tempo2_result['post'] * u.s

Plot the TEMPO2 - PINT residual difference

[28]:

plt.figure(figsize=(8,5), dpi=150)
plt.subplot(2, 1, 1)
plt.plot(mjds, (tp2_diff_pre - tp2_diff_pre.mean()).to_value(u.ns), '+')
plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (ns)')
plt.title('PSR J1600-3053 prefit residual differences between PINT and TEMPO2')
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(mjds, (tp2_diff_post - tp2_diff_post.mean()).to_value(u.ns), '+')
plt.xlabel('MJD (day)')
plt.ylabel('Time Residuals (ns)')
plt.title('PSR J1600-3053 postfit residual differences between PINT and TEMPO2')
plt.grid(True)
plt.tight_layout()
plt.savefig("J1600_PINT_tempo2")
../_images/examples-rendered_paper_validation_example_52_0.png

In this comparison, PINT and TEMPO2’s results, both pre-fit and post-fit, agree with each other within the level of 5 ns.

Write out the TEMPO2 postfit parameter to a new file

  • Note, since the ECL parameter is hard coded in tempo2, we will have to add it manually

[29]:

# Write out the post fit tempo parfile.
tempo2_parfile = open(psr + '_new_tempo2.2.par', 'w')
for line in tempo2_new_par:
    tempo2_parfile.write(line)
tempo2_parfile.write("ECL IERS2003")
tempo2_parfile.close()

Compare the parameter between TEMPO2 and PINT

  • Reported quantities

    • TEMPO2 value

    • TEMPO2 uncertainty

    • Parameter units

    • TEMPO2 parameter value - PINT parameter value

    • TEMPO2/PINT parameter absolute difference divided by TEMPO2 uncertainty

    • PINT uncertainty divided by TEMPO2 uncertainty, if TEMPO2 provides the uncertainty value

[30]:

# Create the parameter compare table
tv = []
t2_unc = []
tv_pv = []
tv_pv_tc = []
tc_pc = []
units = []
names = []
no_t2_unc = []
tempo2_new_model = models.get_model(psr + '_new_tempo2.2.par')
for param in tempo2_new_model.params:
    t2_par = getattr(tempo2_new_model, param)
    pint2_par = getattr(f_t2.model, param)
    tempo2q = t2_par.quantity
    pint2q = pint2_par.quantity
    try:
        diff2q =  tempo2q - pint2q
        if t2_par.uncertainty_value != 0.0:
            diff_tcq = np.abs(diff2q) / t2_par.uncertainty
            uvsu = pint2_par.uncertainty / t2_par.uncertainty
            no_t2_unc.append(False)
        else:
            diff_tcq = np.abs(diff2q) / pint2_par.uncertainty
            uvsu = t2_par.uncertainty
            no_t2_unc.append(True)
    except TypeError:
        continue
    uvsu = pint2_par.uncertainty / t2_par.uncertainty
    tv.append(tempo2q.value)
    t2_unc.append(t2_par.uncertainty.value)
    tv_pv.append(diff2q.value)
    tv_pv_tc.append(diff_tcq.value)
    tc_pc.append(uvsu)
    units.append(t2_par.units)
    names.append(param)

compare_table2 = Table((names, tv, t2_unc,units, tv_pv, tv_pv_tc, tc_pc, no_t2_unc), names = ('name', 'Tempo2 Value', 'T2 unc','units',
                                                                                      'Tempo2_V-PINT_V',
                                                                                      'Tempo2_PINT_diff/unct',
                                                                                      'PINT_unct/Tempo2_unct',
                                                                                      'no_t_unc'))
compare_table2.sort('Tempo2_PINT_diff/unct')
compare_table2 = compare_table2[::-1]
compare_table2.write('parameter_compare.t2.html', format='html', overwrite=True)

WARNING: EPHVER 5 does nothing in PINT [pint.models.timing_model]
WARNING: Unrecognized parfile line 'NE_SW          0' [pint.models.timing_model]
WARNING: Unrecognized parfile line 'NE_SW2 0.000' [pint.models.timing_model]
/home/luo/.local/lib/python3.6/site-packages/astropy/units/quantity.py:477: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

Print the parameter difference in a table.

The table is sorted by relative difference in descending order.

[31]:

compare_table2

[31]:
Table length=125
nameTempo2 ValueT2 uncunitsTempo2_V-PINT_VTempo2_PINT_diff/unctPINT_unct/Tempo2_unctno_t_unc
str8str48float128objectfloat128float128float128bool
ECC0.000173729661575211688.922286680669999241e-094.158524756102052744e-110.00466082844559504490961.0000400945849889922False
DMX_00980.00133946131224894170.00019579968831114546654pc / cm3-5.203720844662671624e-070.0026576757550264501830.999999270898346615False
DMX_0070-0.000237479639065179730.00019767137320477682749pc / cm3-4.67328571424551701e-070.00236416919581180140151.0000006149827960211False
DMX_00970.00139283306619874460.00019620100461426303326pc / cm3-4.3986846819175744183e-070.0022419277060102306420.9999998585691036723False
DMX_0055-0.00053077049044036210.00019675128861832102923pc / cm3-3.9826563099699431592e-070.00202420850096485076261.0000000155754873443False
DMX_0063-0.000484105710728255740.00019894769104906708185pc / cm3-3.8863084052005031355e-070.0019534322739347656681.0000001764728303488False
DMX_00790.000189767952940002160.00019490725481464179483pc / cm3-3.6899448284041791446e-070.0018931798264324964811.000000186999052243False
DMX_00100.000674033569559790.00020051850482404336064pc / cm3-3.7725252445390284467e-070.00188138508605450180750.99999987947730495375False
F1-7.3387383041227678664e-164.619148404392432094e-21Hz / s-8.1916189706268323405e-240.00177340458748570907890.9999982063295332385False
DMX_00860.000295253466908306440.0001961188165133768578pc / cm3-3.4520664969226941277e-070.0017601913769896913431.0000003907045813545False
........................
DMX_00920.00132072951385398940.00019585216459019454951pc / cm3-3.4341162511537445812e-080.000175342266874576870011.0000000634867269866False
DMX_0058-0.00053775814687447930.00019927530538964258904pc / cm3-3.1767895307833332597e-080.000159417120178125611311.0000003635542755731False
DMX_00170.00017898677292651840.00021197504981685315979pc / cm32.854267636918103937e-080.0001346511129203270150.9999995235765041235False
DMX_00890.00074956144462958460.00021586616414944812654pc / cm32.4131246270523223907e-080.0001117879977420487084541.0000000332176062212False
DMX_0027-0.000182880825351814140.00019391445756469536201pc / cm3-1.7893991355739631913e-089.227775783437751106e-051.000000166662465162False
DMX_0040-0.00052424493853935320.00020212647115737782458pc / cm3-1.6804078603323106822e-088.313645663085500665e-050.99999999749641554914False
DMX_0032-6.265675469663684e-050.00019561483985536690729pc / cm31.5367171419820358814e-087.855831097059144273e-051.0000001097328992117False
DMX_00010.00164843721682323250.00022434462780433157077pc / cm3-1.4668063582322365956e-086.538183564224023551e-051.0000004801524928766False
DMX_00838.544780315309648e-060.00020486177918444288125pc / cm3-9.68929340697572577e-094.729673561143965028e-051.0000001554708373153False
DMX_0044-0.00033900236624910280.00021062295971768858391pc / cm3-4.0461326309814554802e-091.9210311337399994e-051.000000188731715367False

If one wants to get the latex version, please use the line below.

[32]:

#ascii.write(compare_table2, sys.stdout, Writer = ascii.Latex,
#            latexdict = {'tabletype': 'table*'})

Check out the maximum DMX difference

[33]:

max_dmx = 0
max_dmx_index = 0
for ii, row in enumerate(compare_table2):
    if row['name'].startswith('DMX_'):
        if row['Tempo2_PINT_diff/unct'] > max_dmx:
            max_dmx = row['Tempo2_PINT_diff/unct']
            max_dmx_index = ii

dmx_max2 = compare_table2[max_dmx_index]['name']

compare_table2[max_dmx_index]

[33]:
Row index=1
nameTempo2 ValueT2 uncunitsTempo2_V-PINT_VTempo2_PINT_diff/unctPINT_unct/Tempo2_unctno_t_unc
str8str48float128objectfloat128float128float128bool
DMX_00980.00133946131224894170.00019579968831114546654pc / cm3-5.203720844662671624e-070.0026576757550264501830.999999270898346615False

Output the table in the paper

[34]:

paper_params = ['F0', 'F1', 'FD1', 'FD2', 'JUMP1', 'PX',
                'ELONG', 'ELAT', 'PMELONG', 'PMELAT', 'PB',
                'A1', 'A1DOT', 'ECC', 'T0', 'OM', 'OMDOT', 'M2',
                'SINI', dmx_max]
# Get the table index of the parameters above
paper_param_index = []
for pp in paper_params:
    # We assume the parameter name are unique in the table
    idx = np.where(compare_table2['name'] == pp)[0][0]
    paper_param_index.append(idx)
paper_param_index = np.array(paper_param_index)
compare_table2[paper_param_index]

[34]:
Table length=20
nameTempo2 ValueT2 uncunitsTempo2_V-PINT_VTempo2_PINT_diff/unctPINT_unct/Tempo2_unctno_t_unc
str8str48float128objectfloat128float128float128bool
F0277.937711242974627885.1859268946902080184e-13Hz-6.6613381477509392425e-160.00128450290237823877811.0000081311054239701False
F1-7.3387383041227678664e-164.619148404392432094e-21Hz / s-8.1916189706268323405e-240.00177340458748570907890.9999982063295332385False
FD13.983282287426775e-051.6566478062738200598e-06s-1.6361212107284211874e-090.00098760955981853026061.00000000325577032False
FD2-1.4729805752137882e-051.1922596055992699934e-06s1.4357622807706935195e-090.00120423628715411431061.0000000135934370427False
JUMP1-8.7887456483184e-060.0s-4.9036265484266971897e-110.00037580263061884339052infTrue
PX0.50612420123220640.07348886965486496614mas1.8776131796016670705e-050.000255496266090326118631.0000000214259234799False
ELONG244.347677842553825.95727548431e-09deg9.1233687271596863866e-120.00153146664967708780971.0000013787890671494False
ELAT-10.0718390470430653.361025894297e-08deg-1.44861900253090425394e-110.00043100501099647159270.99999279185393830885False
PMELONG0.46190960156254910.010433361011620021289mas / yr7.4203184182164427796e-060.00071121074119376856761.0000025748109051538False
PMELAT-7.1551456742758220.058156247552489513664mas / yr-7.1705189014892312116e-050.00123297482269938514160.99999266738091230344False
PB14.3484657546613667862.12226632065849e-06d-1.9243153815198810186e-090.00090672662652574659070.99999757655418149095False
A18.801653122864638.114047416773300209e-07ls-8.196980871844061767e-100.00101022097244580079220.99998413921662954174False
A1DOT-4.008979189463729e-156.2586911221949290846e-16ls / s-1.0335518582171720488e-180.001651386588726085591.0000003602893390298False
ECC0.000173729661575211688.922286680669999241e-094.158524756102052744e-110.00466082844559504490961.0000400945849889922False
T055878.26189947384950700.00051676746764245482d4.8831633153723075225e-070.000944944026304478853571.0000114310932686899False
OM181.849604015494514780.01296564244572522874deg1.225778713571934464e-050.000945405303827481071871.0000120115751260204False
OMDOT0.00523955285176455407780.00135543635075636363deg / yr-1.228657154944533207e-060.000906466138567786814870.9999975920019996213False
M20.27176338143833560.08941866471282471085solMass0.000104263806940840808580.00116601838414488148851.0000265499521452384False
SINI0.90642005682258460.03399283139781983376-4.2776490711826653524e-050.00125839740182896817761.0000279809291119371False
DMX_00100.000674033569559790.00020051850482404336064pc / cm3-3.7725252445390284467e-070.00188138508605450180750.99999987947730495375False

The residual difference between PINT and TEMPO2 is at the level of ~1ns

  • We believe the discrepancy is mainly from the solar system geometric delay.

  • We will use the tempo2 postfit parameters, which are wrote out to J1600-3053_new_tempo2.2.par

[35]:

tempo2_result2 = t2u.general2('J1600-3053_new_tempo2.2.par', tim_file, ['sat', 'pre', 'post', 'freq', 'err'])
m_t22 = models.get_model('J1600-3053_new_tempo2.2.par')
f_t22 = GLSFitter(toas=t, model=m_t22)
f_t22.fit_toas()
tp2_diff_pre2 = f_t22.resids_init.time_resids - tempo2_result2['pre'] * u.s
tp2_diff_post2 = f_t22.resids.time_resids - tempo2_result2['post'] * u.s

WARNING: EPHVER 5 does nothing in PINT [pint.models.timing_model]
WARNING: Unrecognized parfile line 'NE_SW          0' [pint.models.timing_model]
WARNING: Unrecognized parfile line 'NE_SW2 0.000' [pint.models.timing_model]

[36]:

PINT_solar = m_t22.solar_system_geometric_delay(t)
tempo2_solar = t2u.general2('J1600-3053_new_tempo2.2.par', tim_file, ['roemer'])

[37]:

diff_solar = PINT_solar + tempo2_solar['roemer'] * u.s
plt.figure(figsize=(8,2.5), dpi=150)
plt.plot(mjds, (tp2_diff_post2 - tp2_diff_post2.mean()).to_value(u.ns), '+')
plt.plot(mjds, (diff_solar - diff_solar.mean()).to_value(u.ns, equivalencies=[(ls, u.s)]), 'x')

plt.xlabel('MJD (day)')
plt.ylabel('Discrepancies (ns)')
#plt.title('PSR J1600-3053 postfit residual differences between PINT and TEMPO2')
plt.grid(True)
plt.legend(['Postfit Residual Differences', 'Solar System Geometric Delay Difference'],
           loc='upper center', bbox_to_anchor=(0.5, -0.3), shadow=True, ncol=2)
plt.tight_layout()
plt.savefig("solar_geo")
../_images/examples-rendered_paper_validation_example_69_0.png