How To

Quick solutions to common tasks (taken from #pint and elsewhere)

How to upgrade PINT:

With pip:

pip install -U pint-pulsar

With conda:

conda update pint-pulsar

How to go to a specific version of PINT:

With pip:

pip install -U pint-pulsar==0.8.4

or similar.

With conda:

conda install pint-pulsar=0.8.4

or similar.

Find data files for testing/tutorials

The data files (par and tim) associated with the tutorials and other examples can be located via pint.config.examplefile() (available via the pint.config module):

import pint.config
fullfilename = pint.config.examplefile(filename)

For example, the file NGC6440E.par from the Time a Pulsar notebook can be found via:

import pint
fullfilename = pint.config.examplefile("NGC6440E.par")

Load a par file

from pint.models import get_model
m = get_model(parfile)

Load a tim file

from pint.toa import get_TOAs
t = get_TOAs(timfile)

Load a tim and par file together

from pint.models import get_model_and_toas
m, t = get_model_and_toas(parfile, timfile)

Get the red noise basis functions and the corresponding coefficients out of a PINT fitter object

Select TOAs

You can index by column name into the TOAs object, so you can do toas["observatory"] or whatever the column is called; and that's an array, so you can do toas["observatory"]=="arecibo" to get a Boolean array; and you can index with boolean arrays, so you can do toas[toas["observatory"]=="arecibo"] to get a new TOAs object referencing a subset.

Avoid "KeyError: 'obs_jupiter_pos' error when trying to grab residuals?"

You need to have the TOAs object compute the positions of the planets and add them to the table.

ts.compute_posvels(ephem,planets=True)

This should be done automatically if you load your TOAs with the get_TOAs() or get_model_and_TOAs()

Convert from ELAT/ELONG <-> RA/DEC if I have a timing model

model.as_ICRS(epoch=epoch)

which will give it to you as a model with AstrometryEquatorial components at the requested epoch. Similarly,

model.as_ECL(epoch=epoch)

does the same for AstrometryEcliptic (with an optional specification of the obliquity).

Add a jump programmatically

PINT can handle jumps in the model outside a par file. An example is:

import numpy as np
from astropy import units as u, constants as c
from pint.models import get_model, get_model_and_toas, parameter
from pint import fitter
from pint.models import PhaseJump
import pint.config

m, t = get_model_and_toas(pint.config.examplefile("NGC6440E.par"),
                          pint.config.examplefile("NGC6440E.tim"))

# fit the nominal model
f = fitter.WLSFitter(toas=t, model=m)
f.fit_toas()

# group TOAs: find clusters with gaps of <2h
clusters = t.get_clusters(add_column=True)

# put in the pulse numbers based on the previous fit
t.compute_pulse_numbers(f.model)
# just for a test, add an offset to a set of TOAs
t['delta_pulse_number'][clusters==3]+=3

# now fit without a jump
fnojump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fnojump.fit_toas()


# add the Jump Component to the model
m.add_component(PhaseJump(), validate=False)

# now add the actual jump
# it can be keyed on any parameter that maskParameter will accept
# here we will use a range of MJDs
par = parameter.maskParameter(
    "JUMP",
    key="mjd",
    value=0.0,
    key_value=[t[clusters==3].get_mjds().min().value,
               t[clusters==3].get_mjds().max().value],
    units=u.s,
    frozen=False,
    )
m.components['PhaseJump'].add_param(par, setup=True)

# you can also do it indirectly through the flags as:
# m.components["PhaseJump"].add_jump_and_flags(t.table["flags"][clusters == 3])

# and fit with a jump
fjump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fjump.fit_toas()

print(f"Original chi^2 = {f.resids.calc_chi2():.2f} for {f.resids.dof} DOF")
print(f"After adding 3 rotations to some TOAs, chi^2 = {fnojump.resids.calc_chi2():.2f} for {fnojump.resids.dof} DOF")
print(f"Then after adding a jump to those TOAs, chi^2 = {fjump.resids.calc_chi2():.2f} for {fjump.resids.dof} DOF")
print(f"Best-fit value of the jump is {fjump.model.JUMP1.quantity} +/- {fjump.model.JUMP1.uncertainty} ({(fjump.model.JUMP1.quantity*fjump.model.F0.quantity).decompose():.3f} +/- {(fjump.model.JUMP1.uncertainty*fjump.model.F0.quantity).decompose():.3f} rotations)")

which returns:

Original chi^2 = 59.57 for 56 DOF
After adding 3 rotations to some TOAs, chi^2 = 19136746.30 for 56 DOF
Then after adding a jump to those TOAs, chi^2 = 56.60 for 55 DOF
Best-fit value of the jump is -0.048772786677935796 s +/- 1.114921182802775e-05 s (-2.999 +/- 0.001 rotations)

showing that the offset we applied has been absorbed by the jump (plus a little extra, so chi^2 has actually improved).

See maskParamter documentation on the ways to select the TOAs.

new pulsar timing tools

below are several tools i use for new pulsar timing, including cleaning TOAs, adding wraps, and fitting parameters

various packages I use as well as some little convenience functions

import numpy as np
import matplotlib.pyplot as plt
import pint.residuals as res
import copy
from pint.models import BinaryELL1, BinaryDD, PhaseJump, parameter, get_model
from pint.simulation import make_fake_toas_uniform as mft
from astropy import units as u, constants as c
def dot(l1,l2):
    return np.array([v1 and v2 for v1,v2 in zip(l1,l2)])
def inv(l):
    return np.array([not i for i in l])

zapping TOAs on given MJDs, before/after some MJD, or within a window of days

def mask_toas(toas,before=None,after=None,on=None,window=None):
    cnd=np.array([True for t in toas.get_mjds()])
    if before is not None:
        cnd = dot(cnd,toas.get_mjds().value > before)
    if after is not None:
        cnd = dot(cnd,toas.get_mjds().value < after)
    if on is not None:
        on=np.array(on)
        for i,m in enumerate(on):
            m=m*u.day
            if type(m) is int:
                cnd = dot(cnd,inv(np.abs(toas.get_mjds()-m).astype(int) == np.abs((toas.get_mjds()-m)).min().astype(int)))
            else:
                cnd = dot(cnd,inv(np.abs(toas.get_mjds()-m) == np.abs((toas.get_mjds()-m)).min()))
    if window is not None:
        if len(window)!=2:
            raise ValueError("window must be a 2 element list/array")
        window = window*u.day
        lower = window[0]
        upper = window[1]
        cnd = dot(cnd,toas.get_mjds() < lower)+dot(cnd,toas.get_mjds() > upper)
    print(f'{sum(cnd)}/{len(cnd)} TOAs selected')
    return toas[cnd]

add integer phase wraps between TOAs

def add_wraps(toas,mjd,sign,nwrap=1):
    cnd = toas.table['mjd'] > Time(mjd,scale='utc',format='mjd')
    if sign == '-':
        toas.table['pulse_number'][cnd] -= nwrap
    elif sign == '+':
        toas.table['pulse_number'][cnd] += nwrap
    else:
        raise TypeError('sign must be "+" or "-"')

plot residuals in phase

def plot_fit(toas,model,track_mode="use_pulse_numbers",title=None,xlim=None,ylim=None):
    rs=res.Residuals(toas,model,track_mode=track_mode)
    fig, ax = plt.subplots(figsize=(12,10))
    if xlim is not None:
        ax.set_xlim(xlim)
    if ylim is not None:
        ax.set_ylim(ylim)
    ax.errorbar(rs.toas.get_mjds().value,rs.calc_phase_resids(), \
                 yerr=(rs.toas.get_errors()*model.F0.quantity).decompose().value,fmt='x')
    ax.tick_params(labelsize=15)

    if title is None:
        ax.set_title('%s Residuals, %s toas' %(model.PSR.value,len(toas.get_mjds())),fontsize=18)
    else:
        ax.set_title(title,fontsize=18)
    ax.set_xlabel('MJD',fontsize=15)
    ax.set_ylabel(f'Residuals [phase, P0 = {((1/model.F0.quantity).to(u.ms)).value:2.0f} ms]',fontsize=15)
    ax.grid()
    return fig, ax

model uncertainties in pulsar phase over a given range of time

def calculate_phase_uncertainties(model, MJDmin, MJDmax, Nmodels=100, params = 'all', error=1*u.us):
    mjds = np.arange(MJDmin,MJDmax)
    Nmjd = len(mjds)
    phases_i = np.zeros((Nmodels,Nmjd))
    phases_f = np.zeros((Nmodels, Nmjd))
    tnew = mft(MJDmin,MJDmax,Nmjd,model=model, error=error)
    pnew = {}
    if params == 'all':
        params = model.free_params
    for p in params:
        pnew[p] = getattr(model,p).quantity + np.random.normal(size=Nmodels) * getattr(model,p).uncertainty
    for imodel in range(Nmodels):
        m2 = copy.deepcopy(model)
        for p in params:
            getattr(m2,p).quantity=pnew[p][imodel]
        phase = m2.phase(tnew, abs_phase=True)
        phases_i[imodel] = phase.int
        phases_f[imodel] = phase.frac
    phases = phases_i+ phases_f
    phases0 = model.phase(tnew, abs_phase = True)
    dphase = phases - (phases0.int + phases0.frac)
    return tnew, dphase

plot uncertainties in phase (uses calculate_phase_uncertainties())

def plot_phase_unc(model,start,end,params='all'):
    if params == 'all':
        print("calculating phase uncertainty due to all parameters")
        plab = 'All params'
        t, dp = calculate_phase_uncertainties(model, start, end)
    else:
        if type(params) is list:
            print("calculating phase uncertainty due to params "+str(params))
            plab = str(params)
            t, dp = calculate_phase_uncertainties(model, start, end, params = params)
        else:
            raise TypeError('"params" should be either list or "all"')

    plt.gcf().set_size_inches(12,10)
    plt.plot(t.get_mjds(),dp.std(axis=0),'.',label=plab)
    dt = t.get_mjds() - model.PEPOCH.value*u.d
    plt.plot(t.get_mjds(), np.sqrt((model.F0.uncertainty * dt)**2 + (0.5*model.F1.uncertainty*dt**2)**2).decompose(),label='Analytic')
    plt.xlabel('MJD')
    plt.ylabel('Phase Uncertainty (cycles)')
    plt.legend()

workflow

When I'm working on new pulsar solutions, my work flow usually follows something like this:

• load pulsar model/TOAs via get_TOAs() and get_model() model = get_model('parfile.par') toas = get_toas('toas.tim')

• make a copy of the TOAs that I will edit (for easy resets) newtoas = toas

• use plot_fit(newtoas), identify bad toas/missing wraps and mask those data, i.e.

  newtoas = mask_toas(newtoas,before=59000,on=59054.321) newtoas.compute_pulse_numbers(model) add_wraps(newtoas,59166,'-') plot_fit(newtoas,model)

• once i have the TOAs i want, i fit the data and look at the resids and see how the model changed:

  f=WLSFitter(newtoas,model,track_mode='use_pulse_numbers') f.model.free_params=['F0'] # XX TOAs, [lowMJD,highMJD] f.fit_toas() plot_fit(newtoas,f.model) f.print_summary() f.model.compare(model,verbosity='check')

• when the model appears to be a good fit, i update the model and add new observations

  model=f.model newtoas=mask_toas(toas,before=58000,on=59054.321) plot_fit(newtoas,model)

• iterate these last two steps until the fit breaks/i run out of data.

some other, less frequently used snippets

Plot frequency vs residuals

rs=res.Residuals(newtoas,f.model)
fig,ax = plt.subplots(figsize=(12,10))
ax.tick_params(labelsize=15)
ax.set_ylabel('Frequency [MHz]',fontsize=18)
ax.set_xlabel('Phase residuals',fontsize=18)

y = newtoas.get_freqs().to('MHz').value
x = rs.calc_phase_resids()

ax.errorbar(x,y,xerr=newtoas.get_errors().to('s').value*f.model.F0.value,elinewidth=2,lw=0,marker='+')

Plot residuals as a function of orbital phase

x = f.model.orbital_phase(newtoas.get_mjds()).value
rs=res.Residuals(newtoas,f.model)
y = rs.calc_phase_resids()

fig, ax = plt.subplots(figsize=(12,10))
ax.tick_params(labelsize=15)
ax.set_xlabel('Orbital Phase',fontsize=18)
ax.set_ylabel('Phase Residuals',fontsize=18)
ax.grid()
for mjd in np.unique(newtoas.get_mjds().astype(int)):
    cnd = dot(newtoas.get_mjds().astype(int) == mjd,newtoas.get_errors().astype(int) <= 125*u.us)
    ax.errorbar(x[cnd],y[cnd],yerr=(newtoas.get_errors().to('s')*f.model.F0.quantity).value[cnd],elinewidth=2,lw=0,marker='+',label=mjd.value)


ax.legend(fontsize=15)

removing/adding ELL1 binary component

if 'BinaryELL1' in model.components:
    model.remove_component('BinaryELL1')

cmp = BinaryELL1()
cmp.PB.value = 10
cmp.EPS1.value = 1e-5
cmp.EPS2.value = 1e-5 
cmp.TASC.value = 59200
cmp.A1.value = 10

model.add_component(cmp,setup=True)

adding DD binary model

cmp = BinaryDD()
cmp.PB.value = 10
cmp.ECC.value = 0.8
cmp.T0.value = 59251.
cmp.OM.value = 269.
cmp.A1.value = 136.9

model.add_component(cmp,setup=True)

adding spindown component

n = 2
model.components['Spindown'].add_param(
    parameter.prefixParameter(
        name='F'+str(n),
        value=0.0,
        units=u.Hz/u.s**n,
        frozen=False,
        parameter_type="float",
        longdouble=True
    ),
    setup=True
)
model.components['Spindown'].validate()