Processing Script for J0645+5158¶
Import statements and load notebook utilities¶
In [1]:
!git checkout master
!git branch
In [2]:
import os
import sys
sys.path.append("/home/jovyan/work/shared/nanograv_timing_2017/util/")
from nanograv_utils import *
sys.path.append("/home/jovyan/work/shared/")
from pypulse.tim import Tim
%matplotlib inline
# Load the notebook utilities
nanoutils = NANOGravUtils()
tempoutils = TempoUtils()
pintutils = PINTUtils()
enterpriseutils = EnterpriseUtils()
User-defined constants here¶
In [3]:
FLAGS = "-G" # Use "-G" for Generalized least squares fit
NITER = 50000 #10000 # How many iterations to run PAL2 for
Set pulsar name, create working par and tim files, working directories¶
In [4]:
EXCISEPARFILE = "J0645+5158.excise.par"
EXCISETIMFILE = "J0645+5158.excise.tim"
PARFILE = "J0645+5158.working.par"
TIMFILE = "J0645+5158.working.tim"
# Copy initial par/tim files into working files
call("cp J0645+5158.dmx.par %s"%EXCISEPARFILE)
call("cp J0645+5158.dmx.tim %s"%EXCISETIMFILE)
Initial tempo run and residuals plot¶
In [5]:
tempoutils.tempo(EXCISEPARFILE,EXCISETIMFILE,c=True)
tempoutils.plot_resids()
TOA excision¶
In [6]:
# TOA excision from outlier analysis
outlierfilenames = np.loadtxt("outlier/J0645+5158.dmx-outliers.txt",unpack=True,usecols=(0,),dtype=np.str,ndmin=1)
outlierchans = np.loadtxt("outlier/J0645+5158.dmx-outliers.txt",unpack=True,usecols=(6,),dtype=np.int,ndmin=1)
outliersubints = np.loadtxt("outlier/J0645+5158.dmx-outliers.txt",unpack=True,usecols=(8,),dtype=np.int,ndmin=1)
outlierpouts = np.loadtxt("outlier/J0645+5158.dmx-outliers.txt",unpack=True,usecols=(12,),dtype=np.float,ndmin=1)
tim = Tim(EXCISETIMFILE)
for toa in tim.toas:
filename,chan,subint = toa.getFilename(),int(toa.get('chan')),int(toa.get('subint'))
for i,outlierfilename in enumerate(outlierfilenames):
if outlierfilename == filename and outlierchans[i] == chan and outliersubints[i] == subint and outlierpouts[i] >= 0.10:
toa.comment(cut="outlier10")
break
tim.save(EXCISETIMFILE)
nanoutils.plot_multipage_pdf("outlier/J0645+5158.dmx-residuals.pdf")#,height=400)
Clean up the par file¶
In [7]:
# Ensure correct DMX binning, create PARFILE and TIMFILE
! ../../util/./dmx_fixer.py J0645+5158.excise.par J0645+5158.excise.tim excise
# Center the EPOCHs
! $TEMPO/util/center_epoch/center_epoch.py J0645+5158.excise.par J0645+5158.excise.tim J0645+5158.excise.par
In [8]:
# Noise parameter addition from outlier analysis
with open("outlier/J0645+5158.dmx-noise.txt",'r') as FILE:
outlierlines = FILE.readlines()
with open(EXCISEPARFILE,'r') as FILE:
lines = FILE.readlines()
output = ""
for line in lines:
if "EFAC" in line or "EQUAD" in line or "ECORR" in line or "RNAMP" in line or "RNIDX" in line: #remove these parameters
continue
output += line
rnoutput = ""
for line in outlierlines: #add these parameters
if "efac" in line:
output += line.replace("efac-","T2EFAC -f ")
elif "equad" in line:
value = 10**(float(line.split()[-1].strip()))
string = line.replace("equad-","T2EQUAD -f ")
output += "%s %0.5f\n"%(" ".join(string.split()[:-1]),value)
elif "jitter" in line:
value = 10**(float(line.split()[-1]))
string = line.replace("jitter-","ECORR -f ")
output += "%s %0.5f\n"%(" ".join(string.split()[:-1]),value)
elif "RN-Amplitude" in line:
value = float(line.split()[1])
RNAMP = tempoutils.convert_to_RNAMP(value)
RNAMP = str("%0.4e"%RNAMP).replace('e','D')
rnoutput += "RNAMP %s\n"%RNAMP
elif "RN-spectral-index" in line:
gamma = float(line.split()[1])
rnoutput += "RNIDX %0.4f\n"%(-gamma)
if not (gamma < 1.0):
output += rnoutput
with open(EXCISEPARFILE,'w') as FILE:
FILE.write(output)
In [9]:
# Copy excise files to working files
call("cp %s %s"%(EXCISEPARFILE,PARFILE))
call("cp %s %s"%(EXCISETIMFILE,TIMFILE))
Run finalize_timing.py¶
In [10]:
! ../../util/./finalize_timing.py -G J0645+5158.working.par J0645+5158.working.tim
In [11]:
nanoutils.plot_multipage_pdf("J0645+5158.summary.pdf")
Setup noise model directory¶
In [12]:
NOISEDIR = "noisemodel/"
CHAINSDIR = "noisemodel/chains/J0645+5158/"
call("mkdir %s"%NOISEDIR)
call("mkdir noise_output")
#call("mkdir %s"%CHAINSDIR)
call("cp %s %s%s"%(TIMFILE,NOISEDIR,TIMFILE))
call("cp %s %s%s"%(PARFILE,NOISEDIR,PARFILE))
Edit par files to remove noise model terms¶
In [13]:
with open("%s%s"%(NOISEDIR,PARFILE),'r') as FILE:
lines = FILE.readlines()
output = ""
for line in lines:
write = True
for param in ["T2EFAC","T2EQUAD","ECORR","RNAMP","RNIDX"]:
if param in line:
write = False
continue
if write:
output += line
with open("%s%s"%(NOISEDIR,PARFILE),'w') as FILE:
FILE.write(output)
Run noise analysis using ENTERPRISE (currently only supporting tempo)¶
In [14]:
enterpriseutils.setup_model(NOISEDIR+PARFILE,NOISEDIR+TIMFILE)
In [15]:
enterpriseutils.setup_sampler(CHAINSDIR)
In [16]:
enterpriseutils.run_sampler(N=NITER)
In [17]:
enterpriseutils.process_noise('%s/chain_1.txt'%CHAINSDIR,"noise_output/t1_noise.txt") #need to add in T2 noise!
In [18]:
enterpriseutils.plot_corner('%s/chain_1.txt'%CHAINSDIR)
Copy the noise parameters into the par file¶
In [19]:
# No need to comment out low-spectral-index red noise as this is taken care of in EnterpriseUtils.process_noise() now
with open("noise_output/t1_noise.txt",'r') as FILE:
noiselines = FILE.readlines()
with open("%s%s"%(NOISEDIR,PARFILE),'r') as FILE: #read in the noisemodel version
lines = FILE.readlines()
output = ""
for line in lines:
if "JUMP" in line:
for noiseline in noiselines:
output += noiseline
output += line
with open("%s"%PARFILE,'w') as FILE: #write out to the main version
FILE.write(output)
Re-run tempo with noise model¶
In [20]:
tempoutils.tempo(PARFILE,TIMFILE,G=True,c=True)
call("mv J0645+5158.par %s"%PARFILE)
Re-run finalize_timing.py¶
In [21]:
! ../../util/./finalize_timing.py -G J0645+5158.working.par J0645+5158.working.tim
In [22]:
nanoutils.plot_multipage_pdf("J0645+5158.summary.pdf")
Push to git¶
In [23]:
!git add J0645+5158.working.par J0645+5158.working.tim process.ipynb
!git commit -m "J0645+5158 run" J0645+5158.working.par J0645+5158.working.tim process.ipynb
!git pull
!git push