Toa-sim crashing

[tsmith2]

I am attempting to use toa-sim in order to take existing NANOGrav tim and par files to construct a ‘synthetic’ PTA with intrinsic red and white noise parameters of my choice (no GWB yet). I have attempted to read in each tim and par file sequentially and then generate the synthetic files. However, after a few pulsars the code balks and prints (on the terminal) either a message that reads ‘Diagonal element of covariance matrix = 0’ when saving the par files. I can then restart the kernel and get par files for all of the pulsars, but his is quite cumbersome, and I am worried that the files I am saving are possibly incorrect.

My first question is: since I am using the NANOGrav par files to generate the synthetic TOAs, can I use the NANOGrav par files when analyzing the synthetic TOAs? If this is the case, the error when saving the par files from within TOA-sim is moot.

If I do need to use the par files saved from within TOA-sim, is there a way to avoid this error from occurring?

[vallis]

Can you share your code? I haven't seen such an error; I think it comes from the tempo2 layer, which is harder to diagnose and fit.

Usually you'd use the new par file because the parameters are re-fit after modifying residuals. But the original par file may not be too far, and it may be OK to refit when you load the synthetic data.

[tsmith2]

I see-- that is helpful to know! Here is the snippet of code which reads in the NG par and tim files and then generates the synthetic TOAs. Since the code crashes every few pulsars I resorted to running this pulsar by pulsar by supplying a 'psrstring' name for each pulsar by hand. The code crashes with the error I listed above when saving the par file. 'LT' refers to libstempo.toasim.

parfiles = sorted(glob.glob(os.path.join(os.path.join(datadir,'par'), psrstring + '*.par')))
timfiles = sorted(glob.glob(os.path.join(os.path.join(datadir,'tim'), psrstring + '*.tim')))

parfiles = [x for x in parfiles if 'J1713+0747_NANOGrav_12yv3.gls.par' not in x]

psr = T.tempopulsar(parfile = parfiles[0],timfile = timfiles[0])
psrname = parfiles[0].split('/')[-1].split('_')[0]

LT.make_ideal(psr)
LT.add_efac(psr,efac=1.0)
LT.add_rednoise(psr,A_RN,gamma_RN,components=30)
LT.add_jitter(psr, 1e-7)
LT.add_equad(psr, 1e-7)
psr.savepar('./sim_PTA_RN_dist/par/{}-simulate.par'.format(psrname))
psr.savetim('./sim_PTA_RN_dist/tim/{}-simulate.tim'.format(psrname))
T.purgetim('./sim_PTA_RN_dist/tim/{}-simulate.tim'.format(psrname))

[jellis18]

@tsmith2, One thought is that you care seeing covariance errors because you are adding 100ns or EQUAD and jitter (which is a lot for some pulsars). When tempo tries to fit the model it probably can't resolve some binary parameters so you get covariance matrix problems.

This is just a guess. Does it fail randomly or for certain pulsars?

[tsmith2]

@jellis18 Yes! I appears to fail randomly.

If I want to simulate realistic TOAs, can you suggest what white noise parameters would be appropriate? I took 100ns by looking at the mean value for EQUAD and ECORR from the distribution of white noise parameters for the NG 12.5yr pulsars.

[jellis18]

@tsmith2 if they fail randomly then I could be wrong.

Its been a long time since I've actually worked with pulsar data but I believe the best way to simulate realistic data is to use the actual efac, equad and ecorr values from the par file and use the flag and flagid arguments in toasim; however, I don't quite remember how that works.

@vallis do you know? I know we used to do this but I can't find any examples in any github pages....

Either, way once this is figured out it would be good to probably add this as an example in the repo here.

[tsmith2]

@jellis18 @vallis I think I see. So, for J1713, the par file gives the white noise parameters for each backend:

efac_L-wide_ASP noisepar(val=1.039, flag='f', flagval='L-wide_ASP')
efac_L-wide_PUPPI noisepar(val=0.983, flag='f', flagval='L-wide_PUPPI')
efac_Rcvr1_2_GASP noisepar(val=1.068, flag='f', flagval='Rcvr1_2_GASP')
efac_Rcvr1_2_GUPPI noisepar(val=1.015, flag='f', flagval='Rcvr1_2_GUPPI')
efac_Rcvr_800_GASP noisepar(val=1.087, flag='f', flagval='Rcvr_800_GASP')
efac_Rcvr_800_GUPPI noisepar(val=1.043, flag='f', flagval='Rcvr_800_GUPPI')
efac_S-wide_ASP noisepar(val=1.133, flag='f', flagval='S-wide_ASP')
efac_S-wide_PUPPI noisepar(val=1.097, flag='f', flagval='S-wide_PUPPI')
equad_L-wide_ASP noisepar(val=0.028323139999999997, flag='f', flagval='L-wide_ASP')
equad_L-wide_PUPPI noisepar(val=0.00419741, flag='f', flagval='L-wide_PUPPI')
equad_Rcvr1_2_GASP noisepar(val=0.04172676, flag='f', flagval='Rcvr1_2_GASP')
equad_Rcvr1_2_GUPPI noisepar(val=0.0036844499999999997, flag='f', flagval='Rcvr1_2_GUPPI')
equad_Rcvr_800_GASP noisepar(val=0.13945122999999998, flag='f', flagval='Rcvr_800_GASP')
equad_Rcvr_800_GUPPI noisepar(val=0.04874981999999999, flag='f', flagval='Rcvr_800_GUPPI')
equad_S-wide_ASP noisepar(val=0.00398816, flag='f', flagval='S-wide_ASP')
equad_S-wide_PUPPI noisepar(val=0.00663685, flag='f', flagval='S-wide_PUPPI')
ecorr_L-wide_ASP noisepar(val=0.08671, flag='f', flagval='L-wide_ASP')
ecorr_L-wide_PUPPI noisepar(val=0.07685, flag='f', flagval='L-wide_PUPPI')
ecorr_Rcvr1_2_GASP noisepar(val=0.05595, flag='f', flagval='Rcvr1_2_GASP')
ecorr_Rcvr1_2_GUPPI noisepar(val=0.07734, flag='f', flagval='Rcvr1_2_GUPPI')
ecorr_Rcvr_800_GASP noisepar(val=0.00872, flag='f', flagval='Rcvr_800_GASP')
ecorr_Rcvr_800_GUPPI noisepar(val=0.19378, flag='f', flagval='Rcvr_800_GUPPI')
ecorr_S-wide_ASP noisepar(val=0.14538, flag='f', flagval='S-wide_ASP')
ecorr_S-wide_PUPPI noisepar(val=0.04972, flag='f', flagval='S-wide_PUPPI')
log10_ared -13.9914
gamma 1.29117
nred 100

So a realistic white noise realization would take the equad, ecorr, and efac for each of the backends, and that can be done using flagid and flags when generating the noise realization.

In this case, we would feed an array for the white noise parameter (one value for each backend)-- would I then set the flagid = 'f' and flags would be another array corresponding to what is called the 'flagval' in the above output?

[tsmith2]

After working on this a bit more, I decided to see if I could just save the par files of unmodified pulsars and found that the code still (randomly) crashes:

 parfiles = sorted(glob.glob(os.path.join(os.path.join(datadir,'par'), '*.par')))
 timfiles = sorted(glob.glob(os.path.join(os.path.join(datadir,'tim'), '*.tim')))

 parfiles = [x for x in parfiles if 'J1713+0747_NANOGrav_12yv3.gls.par' not in x]
 
 for p, t in zip(parfiles, timfiles):
     psrname = p.split('/')[-1].split('_')[0]
     psr = T.tempopulsar(parfile = p,timfile = t)
 
     psr.savepar('./temp/par/{}-simulate.par'.format(psrname))

This code still, randomly, crashes with the error 'Diagonal element of covariance matrix = 0' printed to the screen.

This makes me wonder if this is due to a memory leak in tempo2.

[sophiehourihane]

@tsmith2 Did you ever fix this error? I am running into it as well.