In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In [2]:
from __future__ import print_function
import sys, math, numpy as N, matplotlib.pyplot as P
Load the libstempo Python extension. It requires a source installation of tempo2, as well as current Python and compiler, and the numpy and Cython packages.
(Both Python 2.7 and 3.4 are supported; this means that in Python 2.7 all returned strings will be unicode strings, while in Python 3 all function arguments should be default unicode strings rather than bytes. This should work transparently, although there are limitations to what characters can be passed to tempo2; you should probably restrain yourself to ASCII.
In [3]:
from libstempo.libstempo import *
In [4]:
import libstempo
In [5]:
libstempo.__path__
Out[5]:
In [6]:
import libstempo as T
T.data = T.__path__[0] + '/data/' # example files
In [7]:
print("Python version :",sys.version.split()[0])
print("libstempo version:",T.__version__)
print("Tempo2 version :",T.libstempo.tempo2version())
We load a single-pulsar object. Doing this will automatically run the tempo2 fit routine once.
In [8]:
psr = T.tempopulsar(parfile = T.data + '/J1909-3744_NANOGrav_dfg+12.par',
timfile = T.data + '/J1909-3744_NANOGrav_dfg+12.tim')
Let's start simple: what is the name of this pulsar? (You can change it, by the way.)
In [9]:
psr.name
Out[9]:
Next, let's look at observations: there are psr.nobs of them; we can get numpy arrays of the site TOAs [in MJDs] with psr.stoas, of the TOA measurement errors [in microseconds] with psr.toaerrs, and of the measurement frequencies with psr.freqs. These arrays are views of the tempo2 data, so you can write to them (but you cannot currently change the number of observations).
In [10]:
psr.nobs
Out[10]:
In [11]:
psr.stoas
Out[11]:
In [12]:
psr.toaerrs
Out[12]:
In [13]:
psr.freqs
Out[13]:
By contrast, barycentric TOAs and frequencies are computed on the basis of current pulsar parameters, so you get them by calling psr methods (with parentheses), and you get a copy of the current values. Writing to it has no effect on the tempo2 data.
In [14]:
psr.toas()
Out[14]:
In [15]:
psr.ssbfreqs()
Out[15]:
Residuals (in seconds) are returned by residuals(). The method takes a few options... I'll let its docstring help describe them. libstempo is fully documented in this way (try help(T.tempopulsar)).
In [16]:
help(psr.residuals)
In [17]:
psr.residuals()
Out[17]:
We can plot TOAs vs. residuals, but we should first sort the arrays; otherwise the array follow the order in the tim file, which may not be chronological.
In [18]:
# get sorted array of indices
i = N.argsort(psr.toas())
# use numpy fancy indexing to order residuals
P.errorbar(psr.toas()[i],psr.residuals()[i],yerr=1e-6*psr.toaerrs[i],fmt='.',alpha=0.2);
We can also see what flags have been set on the observations, and what their values are. The latter returns a numpy vector of strings. Flags are not currently writable.
In [19]:
psr.flags()
Out[19]:
In [20]:
psr.flagvals('chanid')
Out[20]:
In fact, there's a commodity routine in libstempo.plot to plot residuals, taking flags into account.
In [21]:
import libstempo.plot as LP
LP.plotres(psr,group='pta',alpha=0.2)
Timing-model parameters can be accessed by using psr as a Python dictionary. Each parameter is a special object with properties val, err (as well as fit, which is true is the parameter is currently being fitted, and set, which is true if the parameter was assigned a value).
In [22]:
psr['RAJ'].val, psr['RAJ'].err, psr['RAJ'].fit, psr['RAJ'].set
Out[22]:
The names of all fitted parameters, of all set parameters, and of all parameters are returned by psr.pars(which='fit'). We show only the first few.
In [23]:
fitpars = psr.pars() # defaults to fitted parameters
setpars = psr.pars(which='set')
allpars = psr.pars(which='all')
print(len(fitpars),len(setpars),len(allpars))
print(fitpars[:10])
The number of fitting parameters is psr.ndim.
In [24]:
psr.ndim
Out[24]:
Changing the parameter values results in different residuals.
In [25]:
# look +/- 3 sigmas around the current value
x0, dx = psr['RAJ'].val, psr['RAJ'].err
xs = x0 + dx * N.linspace(-3,3,20)
res = []
for x in xs:
psr['RAJ'].val = x
res.append(psr.rms()/1e-6)
psr['RAJ'].val = x0 # restore the original value
P.plot(xs,res)
Out[25]:
We can also call a least-squares fitting routine, which will fit around the current parameter values, replacing them with their new best values. Individual parameters can be included or excluded in the fitting by setting their 'fit' field. (Note: as of version 2.3.0, libstempo provides its own fit, although it does call tempo2 to compute the design matrix.)
In [26]:
psr['DM'].fit
Out[26]:
In [27]:
psr['DM'].fit = True
print(psr['DM'].val)
In [35]:
ret = psr.fit()
In [37]:
print(psr['DM'].val,psr['DM'].err)
The fit returns a tuple consisting of best-fit vector, standard errors, covariance matrix, and linearized chisq. Note that these vectors and matrix are (ndim+1)- or (ndim+1)x(ndim+1)-dimensional, with the first row/column corresponding to a constant phase offset referenced to the first TOA (even if that point is not used).
The exact chisq can be recomputed by psr.chisq() (which evaluates N.sum(psr.residuals()**2 / (1e-12 * psr.toaerrs**2))).
The pulsar parameters can be read in bulk by calling psr.vals(which='fit'), which will default to fitted parameters, but can also be given 'all', 'set', or even a list of parameter names.
In [38]:
fitvals = psr.vals()
print(fitvals)
In [39]:
psr.vals(which=['RAJ','DECJ','PMRA'])
Out[39]:
To set parameter values in bulk, you give a first argument to vals. Or call it with a dictionary.
In [40]:
psr.vals([5.1,-0.6],which=['RAJ','DECJ','PMRA'])
psr.vals({'PMRA': -9.5})
print(psr.vals(which=['RAJ','DECJ','PMRA']))
# restore original values
psr.vals(fitvals)
Be careful about loss of precision; tempopar.val is a numpy longdouble, so you should be careful about assigning it a regular Python double. By contrast, doing arithmetics with numpy longdoubles will preserve their nature and precision.
You can access errors in a similar way with psr.errs(...).
It's also possible to obtain the design matrix computed at the current parameter values, which has shape psr.nobs * (len(psr.pars) + 1), since a constant offset is always included among the fitting parameters.
In [41]:
d = psr.designmatrix()
These, for instance, are the derivatives with respect to RAJ and DECJ, evaluated at the TOAs.
In [42]:
# we need the sorted-index array compute above
P.plot(psr.toas()[i]/365.25,d[i,1],'-x'); P.hold(True)
P.plot(psr.toas()[i]/365.25,d[i,2],'-x')
Out[42]:
It's easy to save the current timing-model to a new par file. Omitting the argument will overwrite the original parfile.
In [43]:
psr.savepar('./foo.par')
In [44]:
!head foo.par
Same for writing tim files.
In [45]:
psr.savetim('./foo.tim')
In [46]:
!head foo.tim
With libstempo, it's easy to replicate some of the "toasim" plugin functionality. By subtracting the residuals from the site TOAs (psr.stoas, vs. the barycentered psr.toas) and refitting, we can create a "perfect" timing solution. (Note that 1 ns is roughly tempo2's claimed accuracy.)
In [47]:
print(math.sqrt(N.mean(psr.residuals()**2)) / 1e-6)
In [48]:
psr.stoas[:] -= psr.residuals() / 86400.0
ret = psr.fit(iters = 4)
In [49]:
print(math.sqrt(N.mean(psr.residuals()**2)) / 1e-6)
Then we can add, e.g., homoskedastic white measurement noise at 100 ns (remember the tempo units: days for TOAs, us for errors, s for residuals).
In [50]:
psr.stoas[:] += 0.1e-6 * N.random.randn(psr.nobs) / 86400.0
psr.toaerrs[:] = 0.1
ret = psr.fit()
In [51]:
i = N.argsort(psr.toas())
P.errorbar(psr.toas()[i],psr.residuals()[i],yerr=1e-6*psr.toaerrs[i],fmt='.')
Out[51]: