<jupyter_start><jupyter_text>DSA simulations<jupyter_code>% matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext line_profiler
%load_ext autoreload
%autoreload 2

from __future__ import division

import numpy as np
import glob, os
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['figure.dpi'] = 2.5 * matplotlib.rcParams['figure.dpi']

import astropy
from astropy.time import Time

import enterprise
from enterprise.pulsar import Pulsar

import enterprise_extensions
from enterprise_extensions import models, model_utils

import libstempo as T2, libstempo.toasim as LT, libstempo.plot as LP
from ephem import Ecliptic, Equatorial

datadir = '../partim_no_noise/'
def figsize(scale):
    fig_width_pt = 513.17 #469.755                  # Get this from LaTeX using \the\textwidth
    inches_per_pt = 1.0/72.27                       # Convert pt to inch
    golden_mean = (np.sqrt(5.0)-1.0)/2.0            # Aesthetic ratio (you could change this)
    fig_width = fig_width_pt*inches_per_pt*scale    # width in inches
    fig_height = fig_width*golden_mean              # height in inches
    fig_size = [fig_width,fig_height]
    return fig_size

#plt.rcParams.update(plt.rcParamsDefault)
params = {'backend': 'pdf',
        'axes.labelsize': 10,
        'lines.markersize': 4,
        'font.size': 10,
        'xtick.major.size':6,
        'xtick.minor.size':3,  
        'ytick.major.size':6,
        'ytick.minor.size':3, 
        'xtick.major.width':0.5,
        'ytick.major.width':0.5,
        'xtick.minor.width':0.5,
        'ytick.minor.width':0.5,
        'lines.markeredgewidth':1,
        'axes.linewidth':1.2,
        'legend.fontsize': 7,
        'xtick.labelsize': 10,
        'ytick.labelsize': 10,
        'savefig.dpi':200,
        'path.simplify':True,
        'font.family': 'serif',
        'font.serif':'Times',
        'text.latex.preamble': [r'\usepackage{amsmath}',r'\usepackage{amsbsy}',
                                r'\DeclareMathAlphabet{\mathcal}{OMS}{cmsy}{m}{n}'],
        'text.usetex':True,
        'figure.figsize': figsize(0.5)}
plt.rcParams.update(params)<jupyter_output><empty_output><jupyter_text>Useful functions<jupyter_code>import pandas as pd
# Get (equatorial) position of pulsars from name
def pos_from_name(name):
    pos = []
    for p in name:
        if '+' in p:
            tmp = p.split('+')
            raj = '{0}:{1}:00.00'.format(tmp[0][1:3],tmp[0][3:])
            decj = '+{0}:{1}:00.00'.format(tmp[1][:2],tmp[1][2:])
        else:
            tmp = p.split('-')
            raj = '{0}:{1}:00.00'.format(tmp[0][1:3],tmp[0][3:])
            decj = '-{0}:{1}:00.00'.format(tmp[1][:2],tmp[1][2:])

        eq = Equatorial(raj, decj)
        if 'B' in p: epoch = '1950'
        elif 'J' in p: epoch= '2000'
        ec = Ecliptic(eq, epoch=str(epoch))
        pos.append([float(eq.ra), float(eq.dec)])
    
    return np.array(pos)

# Get name of pulsar from (equatorial) position
def name_from_pos(pos):
    name = []
    for p in pos:
        eq = Equatorial(p[0], p[1])
        ra = str(eq.ra) 
        dec = str(eq.dec)
        #
        if float(ra.split(':')[0]) < 10:
            ra = '0' + ''.join(ra.split(':')[:2])
        else:
            ra = ''.join(ra.split(':')[:2])
        #
        if np.abs(float(dec.split(':')[0])) < 10:
            dec = ':'.join(['-0'+dec.split(':')[0][1:],dec.split(':')[1]])
        else:
            dec = ':'.join(dec.split(':')[:2])
        if float(dec.split(':')[0]) < 0:
            dec = ''.join(dec.split(':')[:2])
        elif float(dec.split(':')[0]) > 0:
            dec = '+' + ''.join(dec.split(':')[:2])
        else:
            dec = ''.join(dec.split(':')[:2])
        #
        name.append('J' + ra + dec)
    
    return np.array(name)
# Draw new random positions based on a sample
def invtran_sample(pos, size):
    x, y = pos[:,0], pos[:,1]                     
    hist, xbin, ybin = np.histogram2d(x, y, bins=(10, 10))
    xbinc = xbin[:-1] + np.diff(xbin)/2.0
    ybinc = ybin[:-1] + np.diff(ybin)/2.0
    cdf = np.cumsum(hist.ravel())
    cdf = cdf / cdf[-1]

    values = np.random.rand(size)
    value_bins = np.searchsorted(cdf, values)
    x_idx, y_idx = np.unravel_index(value_bins,
                                    (len(xbinc),
                                     len(ybinc)))
    delta_x = np.diff(xbin)[0]
    delta_y = np.diff(ybin)[0]
    if size == 1:
        random_from_cdf = [xbinc[x_idx][0] + np.random.uniform(-delta_x/2.0, delta_x/2.0),
                           ybinc[y_idx][0] + np.random.uniform(-delta_y/2.0, delta_y/2.0)]
    else:
        random_from_cdf = np.column_stack((xbinc[x_idx] + np.random.uniform(-delta_x/2.0, delta_x/2.0),
                                           ybinc[y_idx] + np.random.uniform(-delta_y/2.0, delta_y/2.0)))
    
    return random_from_cdf
def year2mjd(year):
    # rounds to nearest year
    return float(Time("{}-01-01T00:00:00".format(str(int(np.rint(year)))), 
                      format='isot').mjd)

def mjd2year(mjd):
    return float(Time(mjd, format='mjd').decimalyear)<jupyter_output><empty_output><jupyter_text>Process data<jupyter_code>data = pd.read_csv('../data/v2/RMSonlyvsTime2018-DSA2000-LoseAO.csv',header=0,skip_blank_lines=True,)
# Correcting for mistake in spreadsheet
for ii,name in data.iterrows():
    print name
    #if name.Observatory == 'GBT' and name.PSR == 'NEWPSR' and name.Epoch2 >= 2025:
    #    data.iloc[ii,6] = name.RMS3
    #    data.iloc[ii,7] = np.nan
    #    data.iloc[ii,8] = np.nan
# names of all real gbt and ao pulsars
real_pulsars_gbt = np.array([name.PSR for ii,name in data.iterrows() if 'NEWPSR' not in name.PSR and name.Observatory=='GBT'])
real_pulsars_ao = np.array([name.PSR for ii,name in data.iterrows() if 'NEWPSR' not in name.PSR and name.Observatory=='AO'])

# names of all fake gbt and ao pulsars (these are all NEWPSR)
fake_pulsars_gbt = np.array([name.PSR for ii,name in data.iterrows() if 'NEWPSR' in name.PSR and name.Observatory=='GBT'])
fake_pulsars_ao = np.array([name.PSR for ii,name in data.iterrows() if 'NEWPSR' in name.PSR and name.Observatory=='AO'])
# Get equatorial positions (ra, dec...in radians) from names
pos_pulsars_gbt = pos_from_name(real_pulsars_gbt)
pos_pulsars_ao = pos_from_name(real_pulsars_ao)
# Generate new fake pulsar positions
fakepos_pulsars_gbt = np.array([invtran_sample(pos_pulsars_gbt, size=1) for ii in range(len(fake_pulsars_gbt))])
fakepos_pulsars_ao = np.array([invtran_sample(pos_pulsars_ao, size=1) for ii in range(len(fake_pulsars_ao))])

# Generate new fake pulsar names
fake_pulsars_gbt = name_from_pos(fakepos_pulsars_gbt)
fake_pulsars_ao = name_from_pos(fakepos_pulsars_ao)
# Make copy of data frame and replace NEWPSR with new names

data_copy = data.copy(deep=True)
ct_gbt = 0
ct_ao = 0
for ii,name in data_copy.iterrows():
    
    # New GBT pulsars
    if name.Observatory == 'GBT' and name.PSR=='NEWPSR':
        #tmp_pos = invtran_sample(pos_pulsars_gbt, size=1)
        #tmp_name = name_from_pos(tmp_pos)
        #data_copy.iloc[ii,0] = tmp_name[0]
        data_copy.iloc[ii,0] = fake_pulsars_gbt[ct_gbt]
        ct_gbt += 1
        
    # New AO pulsars
    if name.Observatory == 'AO' and name.PSR=='NEWPSR':
        #tmp_pos = invtran_sample(pos_pulsars_ao, size=1)
        #tmp_name = name_from_pos(tmp_pos)
        #data_copy.iloc[ii,0] = tmp_name[0]
        data_copy.iloc[ii,0] = fake_pulsars_ao[ct_ao]
        ct_ao += 1
# Create new data columns for RAJ and DECJ

data_copy = data_copy.assign(RAJ=pd.Series(np.concatenate([pos_pulsars_gbt[:,0],
                                                           fakepos_pulsars_gbt[:,0],
                                                           pos_pulsars_ao[:,0],
                                                           fakepos_pulsars_ao[:,0]])).values)

data_copy = data_copy.assign(DECJ=pd.Series(np.concatenate([pos_pulsars_gbt[:,1],
                                                           fakepos_pulsars_gbt[:,1],
                                                           pos_pulsars_ao[:,1],
                                                           fakepos_pulsars_ao[:,1]])).values)<jupyter_output><empty_output><jupyter_text>Dealing with Burning Dumpster sim<jupyter_code>data = pd.read_csv('../data/v2/RMSonlyvsTime2018-DSA2000-BurningDumpster.csv',header=0,skip_blank_lines=True,)
data_copy = data.copy(deep=True)
final_data = pd.read_csv('../data/v2/status_quo.csv')
data_copy.PSR = final_data.PSR
data_copy = data_copy.assign(RAJ=final_data.RAJ)
data_copy = data_copy.assign(DECJ=final_data.DECJ)
# Write to csv
data_copy.to_csv('../data/v2/burning_dumpster.csv')<jupyter_output><empty_output><jupyter_text>Creating par files<jupyter_code>sim_type = 'status_quo' # 'status_quo', 'lose_ao', or 'burning_dumpster'
# Read from csv
data_copy = pd.read_csv('../data/v2/{}.csv'.format(sim_type))
for ii,name in data_copy.iterrows():
    
    # Get coordinates for par file
    tmp = Equatorial(data_copy.iloc[ii].RAJ, data_copy.iloc[ii].DECJ, epoch='2000')
    
    # Set PEPOCHs to be 5 years in
    rms = np.array([name.RMS1, name.RMS2, name.RMS3, 
                    name.RMS4, name.RMS5])
    rms[rms=='gap'] = np.inf
    rms = np.array(rms,dtype=float)
    epoch = np.array([name.Epoch1, name.Epoch2, name.Epoch3, 
                      name.Epoch4, name.Epoch5])
    start = epoch[np.where(~np.isnan(rms))[0][0]]
    
    with open('../data/template.par', 'r') as fil:
        pardata = fil.read()
        
    with open('../data/v2/par/' + name.PSR + '.par', 'w') as filnew:
        for line in pardata.split('\n'):
            if 'PSR' in line:
                print >>filnew, '\t\t'.join([line.split()[0], name.PSR])
            elif 'RAJ' in line:
                print >>filnew, '\t\t'.join([line.split()[0], str(tmp.ra), 
                                         line.split()[2], line.split()[3]])
            elif 'DECJ' in line:
                print >>filnew, '\t\t'.join([line.split()[0], str(tmp.dec), 
                                         line.split()[2], line.split()[3]])
            elif 'PEPOCH' in line or 'POSEPOCH' in line or 'DMEPOCH' in line:
                print >>filnew, '\t\t'.join([line.split()[0], str(int(year2mjd(start+5)))])
            else:
                print >>filnew, line
    print name.PSR, start<jupyter_output>/Users/staylor/miniconda2/lib/python2.7/site-packages/ipykernel_launcher.py:9: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if __name__ == '__main__':<jupyter_text>Creating tim files<jupyter_code>sim_type = 'burning_dumpster' # 'status_quo', 'lose_ao', or 'burning_dumptser'
# Read from csv
data_copy = pd.read_csv('../data/v2/{}.csv'.format(sim_type))
dsa_sims = []
start_data = []
for jj,name in data_copy.iterrows():
    
    psrname = name.PSR
    
    rms = np.array([name.RMS1, name.RMS2, name.RMS3, 
                    name.RMS4, name.RMS5])
    rms[rms=='gap'] = np.inf
    rms = np.array(rms,dtype=float)
    epoch = np.array([name.Epoch1, name.Epoch2, name.Epoch3, 
                      name.Epoch4, name.Epoch5])
    
    ### Start and End year
    start_yr = epoch[np.where(~np.isnan(rms))[0][0]]
    start_yr_mjd = year2mjd(start_yr)
    #
    end_yr = 2045
    end_yr_mjd = year2mjd(end_yr)
    
    ### Spacing and obstimes
    spacing = 365.25 / 20.0 # days between observations
    #
    obstimes = np.arange(start_yr_mjd, end_yr_mjd, spacing)
    # removing data gaps
    for kk,rmss in enumerate(rms):
        if np.isinf(rmss):
            mask = np.logical_and(obstimes >= year2mjd(epoch[kk]), 
                                  obstimes <= year2mjd(epoch[kk+1]))
            obstimes = obstimes[~mask]
    
    ### Segmenting obstimes based on hardware/telescope switches
    stops = list(epoch[np.where(~np.isnan(rms))[0]]) + [end_yr]
    stops = [year2mjd(yr) for yr in stops]
    errors = list(rms[np.where(~np.isnan(rms))[0]])
    
    ### Masking sections of data based on these stops
    masks = []
    for kk,stop in enumerate(stops):
        if kk < len(stops)-1:
            masks.append(np.logical_and(obstimes >= stops[kk], 
                                        obstimes <= stops[kk+1]))
    
    ### Applying RMS errors
    toa_errs = np.ones_like(obstimes)
    for kk,mask in enumerate(masks):
        toa_errs[mask] *= float(errors[kk])
    
    ### Make fake dataset
    par = '../data/v2/par/' + psrname + '.par'
    dsa_sims.append(LT.fakepulsar(parfile=par, obstimes=obstimes, 
                                  toaerr=toa_errs, 
                                  observatory=name.Observatory.lower()))
    # white noise
    LT.add_efac(dsa_sims[jj])
    # save .tim
    dsa_sims[jj].savetim('../data/v2/tim_{}/'.format(sim_type) 
                         + dsa_sims[jj].name + '.tim')
    ###
    
    start_data.append([psrname, start_yr, start_yr_mjd])
    print psrname, par, start_yr_mjd, end_yr_mjd, len(stops), len(masks), len(errors)
start_data = np.array(start_data)
start_data[start_data[:,1].argsort()]

fil = open('sims_psr_startdata_{}.txt'.format(sim_type),'w')
for line in start_data[start_data[:,1].argsort()]:
    print >>fil, line[0], line[1], line[2]
fil.close()<jupyter_output><empty_output><jupyter_text>Read In And Check Pulsars<jupyter_code>import enterprise
from enterprise.pulsar import Pulsar
from enterprise.signals import parameter
from enterprise.signals import white_signals
from enterprise.signals import gp_signals
from enterprise.signals import signal_base

import enterprise_extensions
from enterprise_extensions import models, model_utils
psr1 = Pulsar('../data/v2/par/J2234+0944.par', '../data/v2/tim_status_quo/J2234+0944.tim', ephem='DE436')
psr2 = Pulsar('../data/v2/par/J2234+0944.par', '../data/v2/tim_lose_ao/J2234+0944.tim', ephem='DE436')
psr3 = Pulsar('../data/v2/par/J2234+0944.par', '../data/v2/tim_burning_dumpster/J2234+0944.tim', ephem='DE436')
plt.errorbar([mjd2year(p) for p in psr1.toas/86400.0], 
             psr1.residuals/1e-6, psr1.toaerrs/1e-6, 
             alpha=0.3, fmt='.')

#plt.errorbar([mjd2year(p) for p in psr2.toas/86400.0], 
#             psr2.residuals/1e-6, psr2.toaerrs/1e-6, 
#             alpha=0.3, fmt='.')

plt.errorbar([mjd2year(p) for p in psr3.toas/86400.0], 
             psr3.residuals/1e-6, psr3.toaerrs/1e-6, 
             alpha=0.3, fmt='.')

plt.xlabel(r'Year')
plt.ylabel(r'Residuals [$\mu$s]')
plt.title(psr1.name)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15242811941119824, 'fraction_numerical': 0.028889806025588115, 'mean_word_length': 2.6041150223103617, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '15811', 'license': 'MIT', 'n_tokens_mistral': 5628, 'n_tokens_neox': 5399, 'n_words': 955, 'path': 'code/nano_dsa_sims_v2.ipynb', 'repo_name': 'LBJ-Wade/pta_simulations', 'size': '14538'}