In [ ]:

    

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext Cython
%load_ext autoreload
%autoreload 2

from __future__ import division
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['savefig.dpi'] = 1.5 * matplotlib.rcParams['savefig.dpi']

import numpy as np
import sys,os,glob,h5py
import cPickle as pickle

import libstempo as T2
import libstempo.plot as LP

import NX01_psr
import NX01_datafile

try:
    from IPython.core.display import clear_output
    have_ipython = True
except ImportError:
    have_ipython = False

Msol = 1.98855*10.0**30.0


    

In [ ]:

    

parpath = '/Users/staylor/Research/IPTA_MDC/IPTA_Challenge1_open/Challenge_Data/Dataset1/'
timpath = '/Users/staylor/Research/IPTA_MDC/IPTA_Challenge1_open/Challenge_Data/Dataset1/'


    

In [ ]:

    

parfiles = sorted(glob.glob(parpath+'*.par'))
timfiles = sorted(glob.glob(timpath+'*.tim'))


    

In [ ]:

    

parfiles = [par for par in parfiles if 'red' not in par]


    

In [ ]:

    

len(parfiles), len(timfiles)


    

In [ ]:

    

#######################################
# PASSING THROUGH TEMPO2 VIA libstempo
#######################################

t2psr = []
for ii in range(len(parfiles)):
    
    t2psr.append( T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii], ) )
    t2psr[ii].fit(iters=10)
    
    if np.any(np.isfinite(t2psr[ii].residuals())==False)==True:
        t2psr[ii] = T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii])
                 
    clear_output()
    print '\r', '{0} of {1}'.format(ii+1,len(parfiles))
    sys.stdout.flush()


    

In [ ]:

    

# Check out some plots if you want
ind = 0
print t2psr[ind].name
plt.errorbar(t2psr[ind].toas(),t2psr[ind].residuals(),1e-6*t2psr[ind].toaerrs,fmt='.',alpha=0.2)
plt.show()


    

In [ ]:

    

psr = [NX01_psr.PsrObj(p) for p in t2psr]
[p.grab_all_vars(jitterbin=1.0,makeGmat=False,
                 fastDesign=True,planetssb=True) for p in psr] # according to the 9 year paper,
                                                               # the jitterbin used is 1s


    

In [ ]:

    

for ii in range(len(psr)):
    psr[ii].parfile = parfiles[ii] 
    psr[ii].timfile = timfiles[ii]


    

In [ ]:

    

dirname = os.getcwd() + '/ipta_mdc_hdf5files/'
if not os.path.exists(dirname):
    os.makedirs(dirname)


    

In [ ]:

    

for ii,p in enumerate(psr):
    df = NX01_datafile.DataFile(dirname + '/' + p.name + '.hdf5')
    df.addTempoPulsar(p)
    
    clear_output()
    print '\r', '{0} of {1}'.format(ii+1,len(psr))
    sys.stdout.flush()


    

In [ ]:

    

fil = open(dirname + '/psrList_iptamdc.txt','w')
print >>fil, "NAME", "HDF5-PATH", "PARFILE-PATH", "TIMFILE-PATH"
print >>fil, "#############################################"
for p in psr:
    print >>fil, p.name, dirname+'/'+p.name+'.hdf5'
fil.close()


    

In [ ]: