import pyfstat import numpy as np import os import sys import time ID = sys.argv[1] outdir = sys.argv[2] label = 'run_{}'.format(ID) data_label = '{}_data'.format(label) results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID) # Properties of the GW data sqrtSX = 1e-23 tstart = 1000000000 Tspan = 100*86400 tend = tstart + Tspan # Fixed properties of the signal F0_center = 30 F1_center = -1e-10 F2 = 0 Alpha = np.radians(83.6292) Delta = np.radians(22.0144) tref = .5*(tstart+tend) VF0 = VF1 = 200 dF0 = np.sqrt(3)/(np.pi*Tspan) dF1 = np.sqrt(45/4.)/(np.pi*Tspan**2) DeltaF0 = VF0 * dF0 DeltaF1 = VF1 * dF1 nsteps = 25 run_setup = [((nsteps, 0), 20, False), ((nsteps, 0), 7, False), ((nsteps, 0), 2, False), ((nsteps, nsteps), 1, False)] h0 = 0 F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0 F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1 psi = np.random.uniform(-np.pi/4, np.pi/4) phi = np.random.uniform(0, 2*np.pi) cosi = np.random.uniform(-1, 1) # Next, taking the same signal parameters, we include a glitch half way through dtglitch = Tspan/2.0 delta_F0 = 0.25*DeltaF0 delta_F1 = -0.1*DeltaF1 glitch_data = pyfstat.Writer( label=data_label, outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1) glitch_data.make_data() startTime = time.time() theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2., 'upper': F0+DeltaF0/2}, 'F1': {'type': 'unif', 'lower': F1-DeltaF1/2., 'upper': F1+DeltaF1/2}, 'F2': F2, 'Alpha': Alpha, 'Delta': Delta, 'tglitch': {'type': 'unif', 'lower': tstart+0.1*Tspan, 'upper': tend-0.1*Tspan}, 'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0}, 'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1}, } ntemps = 2 log10temperature_min = -0.1 nwalkers = 100 nsteps = [500, 500] glitch_mcmc = pyfstat.MCMCGlitchSearch( label=label, outdir=outdir, sftfilepath='{}/*{}*sft'.format(outdir, data_label), theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, log10temperature_min=log10temperature_min) glitch_mcmc.run(run_setup=run_setup, create_plots=False, log_table=False, gen_tex_table=False) glitch_mcmc.print_summary() d, maxtwoF = glitch_mcmc.get_max_twoF() dF0 = F0 - d['F0'] dF1 = F1 - d['F1'] tglitch = d['tglitch'] R = (tglitch - tstart) / Tspan delta_F0 = d['delta_F0'] delta_F1 = d['delta_F1'] runTime = time.time() - startTime with open(results_file_name, 'a') as f: f.write('{:1.8e} {:1.8e} {} {:1.8e} {:1.8e} {:1.8e} {}\n' .format(dF0, dF1, R, delta_F0, delta_F1, maxtwoF, runTime)) os.system('rm {}/*{}*'.format(outdir, label))