Timing Part 2: Refining a known timing solution¶

Most of the time in PTA science, we're not starting from scratch -- someone else has already built a good timing solution for the pulsars we're analyzing.

However, every time we add more data, we need to re-fit the timing model. This is partly to make sure nothing has happened that makes our model less accurate than we expected, but it's also partly because the more data we add, the more detailed our timing model can become.

Remember, when we first discovered a pulsar, we only had a few pieces of information, but we could add more just by having a couple hundred TOAs. PTA pulsars can have thousands or tens of thousands of TOAs -- this allows us to make our datasets more and more robust.

In the first section of the tutorial, we'll load a small set of TOAs with a slightly out-of-date timing model and update the model. In the second part of the tutorial, we'll load a real IPTA dataset and update the timing model.

After you complete a guided version of Tutorial 2, you can use Tutorial 2.5 (a clean version with fewer comments) to check solutions for the additional pulsars you will have in your PTA. There are par, tim, and config files included for 12 pulsars, but you'll probably want to pick one pulsar per group member rather than timing all of the pulsars.

Loading the relevant sofware¶

In [ ]:

# This cell imports the software we'll use in our analysis
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from io import StringIO
import pint.fitter
from pint.models import get_model
from pint.toa import get_TOAs
from pint.residuals import Residuals
from pint.simulation import make_fake_toas_uniform
import pint.logging
# Advanced User Note: if you would like to see the DEBUG and INFO statements associated with PINT functions, 
# comment out the following line. These are not necessary for learners, but can be useful for experts.
pint.logging.setup("ERROR")

Timing a known pulsar (Simple version)¶

When we already have a pretty good timing solution, we won't build our par file from scratch like we did in the first tutorial. Instead, we'll use PINT tools to load a par and tim file. (Along the way, we are essentially going to re-build the function do_timing that we had in Tutorial 1).

We'll begin by setting up variables that represent the par and tim file names. Our files live in folders called par and tim respectively and are named NGC6440E.par and NGC6440E.tim.

Question: \ Why is it valuable to use variables for the par and tim filenames rather than typing them into the commands for loading models and data directly?

In [ ]:

parfile = 
timfile = 

Now use PINT to load the timing model and the TOAs with get_model and get_TOAs respectively.

Note that we want to specify an ephemeris in get_TOAs so that t = get_TOAs(your-tim, ephem='DE440'). We don't need any additional options with get_model.

Activity:\ Load the TOAs and model found in the tim and par file referenced above.

Each of these objects has properties we can access. Within the TOAs, we can find the modified Julian Days that correspond to the observations with toas.get_mjds(). We can also find the error in each TOA using toa.get_errors(). For better units, formatting, we more commonly use the format toas.get_errors().toas(u.us).value. This will save the values in microseconds

Activity \ In the cell below, find the MJDs and the errors for the TOAs imported above.

Questions:

What is the first MJD in the tim file? What is the last MJD?
What is the average TOA error value?

We have tons of info in our model too -- perhaps the easiest way to access it is just to call the model in our notebook. This will show us all the parameters, what kind of parameter they are, the value for the parameter, if they're fit, and the units of each parameter

Activity\ Show our model below (hint -- don't do anything too sophisticated, just type the model's variable name and run the cell).

When we get a par file for a PTA pulsar, we know that at some time in the past, this model was a good fit to the data (we can often figure out roughly when by looking at the START and END parameters in the par file). But we also know that in any model fitting problem, if we add more data points, the best fit can change. So, our first task is to see how well the old model fits this data.

To plot any model in PINT, we have to set up a fitter object. It is this fitter object that allows us to view residuals as well as to do the actual fitting process. PINT fitters are found in the PINT class pint.fitter.

There are several types of fitter implemented in PINT, most of which are variations on a weighted least squares fit. In PTA analyses, the best option is generally the DownhillGLSFitter. This is a generalized least squares fitter (which unlike the very similar DownhillWLSFitter is compatible with correlated noise parameters like ECORR). The "Downhill" part refers to the degree to which the fitter requires convergence.

The syntax for instantiating a fitter is my_fitter = pint.fitter.DownhillGLSFitter(toas,model).

If you're not sure what fitter to use, you can also use pint.fitter.Fitter.auto(toas, model.

Activity \ Using the TOAs and model you've already loaded, set up a Downhill generalized least squares fitter in the space below.

Once we create a fitter, we are able to assess how well our model fits our data. We haven't done a fit yet, but this is a little quirk of PINT: that all the information about combining the model and data can be found in the fitter object. It is also possible to evaluate the some of this information (like the residuals between a model and data) without a fitter object, but frankly, it can be kind of a waste of time if we know we eventually want to fit a new model to the data!

To access the residuals without setting up a fitter, we would use the Residuals class, e.g. resids = Residuals(toas, model). Alternatively, any fitter object creates a residuals object, so we can access the residuals via the my_fitter.resids.

Either way, residuals can be expressed either in terms of time (seconds) through time_resids or in terms of phase (turns) through phase_resids.

A sample call, for a set of TOAs called toas and a fitter called my_fitter, presented in units of microseconds (u.us) would look like

resids = my_fitter.resids.time_resids.toas(u.us).value

Activity\ Generate pre-fit residuals for your fittter.

Find the average pre-fit residual value.

In [ ]:

rs_prefit =

Next, we'd like to make a plot. We want the plot to have Modified Julian Days on the x-axis, residuals on the y-axis, and we want the data points to have y-error bars. Because we want error bars, we'll use plt.errorbar(x,y,yerr=errors) instead of plt.plot(x,y).

Activity\ Make a residuals plot. Hint: You can look at Tutorial 1 for help on syntax.

Assess the fit of your model to your data. \ Question\ How does this look? Do you think this is a good fit? Why or why not?

Now, let's actually run the fit. This is pretty straightforward -- we just type my_fitter.fit_toas(). We can also set a higher number of iterations for the fiter with the option maxiter

Activity \ Run the fit!

Now, save the post-fit residuals and plot again to see how our fitter is doing.

Question\ Does this look like a better fit? Why or why not?

We need to understand how the parameters changed, and we'd really like a more definite answer for "How good was the fit?" than "Hmmm, looks nice?" Fortunately, PINT has a built-in option for this.

The fitter has a method print_summary, which we can access via my_fitter.print_summary(). This will print the number of TOAs present in the data, the number of parameters fit, the pre-fit weighted RMS residual & the post-fit weighted RMS residuals, the reduced chi-squared value for the fit, and a comparison of the initial and final values.

Activity\ In the cell below, print the summary of your fit.

Question

Based on the results of print_summary is your current solution a good fit to this data? How do you know?
What parameters changed in the fitting process? By how much?

If we are content with our analysis, it's time to save our new model. This will allow us to pass along our solution to another PTA member or to come back to these data later and not have to duplicate effort.

The command in PINT to save a new par file is generally model.write_parfile('filename.par'). However, we want to be sure we write out the new model we found with fitting, not our original input model. Therefore, we instead use the model inside the fitter -- accessed via my_fitter.model.

Activity\ In the space below, save your new model as a par file.

Doing this analysis "for real" -- using PINT Pal¶

If you actually get involved in either NANOGrav PINT-based timing or IPTA data combination, you'll use a supplemental software package called PINT-Pal. PINT-Pal includes tools that make it easier to handle a large set of data, to make plots, to make nice summary documents and so on.

One of the biggest changes in a PINT-Pal framework vs. a pure PINT framework is the addition of the "config" file. The "config" file is a text file in YAML (yet another markup language) format which saves information about the desired par file, tim file(s), fitting parameters, conventions for analysis, bad data points, and even noise analysis.

In addition to inputing file paths, we also use the config file to tell PINT what parameters to fit. All parameters we want fit should be put in the line free-params as well as having a 1 next to the quantity in the par file. This seems a little clunky at first, but it has the huge advantage of providing an "at a glance" list of fit parameters.

The PINT functionality is the same when using PINT-Pal; we've just hidden some of the details so that a big/difficult problem is easier to get our head around.

This portion of the tutorial is just a quick re-do of the above tutorial with PINT-Pal functions, but if you plan to get involved in data combination, make sure to check out the separate dr3_student_workshop tutorial available here: https://github.com/gooddc/dr3_student_workshop.

In [ ]:

import pint_pal as pp
import pint_pal.lite_utils as lu
import pint_pal.plot_utils as pu
from pint_pal.timingconfiguration import TimingConfiguration
from astropy import log
from pint.fitter import ConvergenceFailure
import pint.fitter
import os
import sys
import copy
from astropy.visualization import quantity_support
quantity_support()

%matplotlib inline

# Set logging level (PINT uses loguru)
log.setLevel("INFO") # Set desired verbosity of log statements (DEBUG/INFO/WARNING/ERROR)
pint.logging.setup(level="ERROR", usecolors=True)

Load data with PINT-Pal¶

Several sample config files are provided in the configs folder. Pick whichever one you like. Open that file and take a look at what it contains.

Then, fill in that config file name below and run the cell.

This cell's primary function is get_model_and_toas, which will (as it says) combine get_model and get_TOAs from PINT. It also has some additional options, like using a pickle file to load a previously used set of model & TOAs. (This cell can take a while for large datasets).

In [ ]:

config = "configs/[yourconfigname].yaml"  # fill in actual path
par_directory = None   # default location
tim_directory = None   # default location
tc = TimingConfiguration(config, par_directory=par_directory, tim_directory=tim_directory)

# To combine TOAs, assumption is that cuts have already been applied properly
mo,to = tc.get_model_and_toas(apply_initial_cuts=False,usepickle=False, )

Now let's do some data housekeeping tasks.

In [ ]:

# Apply manual cuts. This line allows the "bad TOAs" to be flagged out in the config file.
tc.manual_cuts(to)

# Computing pulse numbers ensures param changes in the model will not break phase connection
to.compute_pulse_numbers(mo)

# Set non-binary epochs to the center of the data span
lu.center_epochs(mo,to)

# Summarize TOAs present
summary = to.get_summary()

Question: Based on the summary above, what telescope do you think was used to collect these TOAs? How many TOAs are there?

Build a fitter and plot¶

Next, we'll construct a fitter object with our model and TOAs (loaded from the par and tim files) and make our first plots. Though we've made a fitter object using the PINT Pal tool construct_fitter, we haven't yet executed a fit -- that's why our plots are labeled as "pre-fit."

These plotting functions plot_residuals_time and plot_residuals_orb are really helpful PINT-Pal tools. The first plots residuals as a function of time, and has built in options to apply different colors and labels based on information in the tim file (like what telescope collected the data or what PTA the data are from). The second does much the same but as a function of binary orbital phase (if the pulsar is in a binary). If the pulsar is not in a binary, this cell creates only one plot.

It's not default, but you can also add save=True to these calls to automatically save the plots.

In [ ]:

# Define the fitter object and plot pre-fit residuals
fo = tc.construct_fitter(to,mo)
pu.plot_residuals_time(fo, restype='prefit', legend=True)
if mo.is_binary:
    pu.plot_residuals_orb(fo, restype='prefit', legend=True)

Question: What if any interesting features do you see in this residual plot? How does it compare to the plot you generated without using PINT-pal?

Fit a new model and evaluate¶

Now, we actually want to fit a new model to these TOAs.

We'll first create a list of free parameters, then we'll go ahead and fit the model, using the fitting tools in PINT. One of the options available in our config is to set the type of fitter used.

Question: Open the config file. What fitter are we using in the following cell?

In [ ]:

# Set free params based on list in the config file (want to update JUMP handling differently soon)
fo.model.free_params = tc.get_free_params(fo)

# Do the fit
try:
    fo.fit_toas(maxiter=tc.get_niter())
    fo.model.CHI2.value = fo.resids.chi2
except ConvergenceFailure:
    run_Ftest = False
    log.warning('Failed to converge; moving on with best result, but should address before final version.')

Our fitter should have converged okay. This step is fairly fast for our small data set, but it can take a long time for a larger dataset!

Now let's examine how our fitter did. First, plots. Notice that now, they're labeled as "postfit".

In [ ]:

# Plot post-fit residuals, print summary of results, write prenoise solution
pu.plot_residuals_time(fo, restype='postfit', legend=True)
if mo.is_binary:
    pu.plot_residuals_orb(fo, restype='postfit', legend=True)

Now we can print a quick summary of our results. This is a great way to examine the difference between our pre-fit and post-fit parameters at a glance. Because all we did was add a bit more data, we shouldn't see a major difference in any of the parameters. If we do, that's a red flag.

The PINT-Pal function check_convergence makes sure that our fitter has indeed converged. It will also warn us if we're not fitting parameters we should be fitting.