diff --git a/examples/glitch_robust_search.py b/examples/glitch_robust_search.py new file mode 100644 index 0000000000000000000000000000000000000000..86e16b7736587f834b4549443c1824993bc0fbf5 --- /dev/null +++ b/examples/glitch_robust_search.py @@ -0,0 +1,54 @@ +import numpy as np +import pyfstat + +outdir = 'data' +label = 'glitch_robust_search' + +# Properties of the GW data +tstart = 1000000000 +Tspan = 60 * 86400 + +# Fixed properties of the signal +F0s = 30 +F1s = -1e-8 +F2s = 0 +Alpha = np.radians(83.6292) +Delta = np.radians(22.0144) + +tref = tstart + .5 * Tspan + +sftfilepath = 'data/*glitching_signal*sft' + +F0_width = np.sqrt(3)/(np.pi*Tspan) +F1_width = np.sqrt(45/4.)/(np.pi*Tspan**2) +DeltaF0 = 50 * F0_width +DeltaF1 = 50 * F1_width + +theta_prior = {'F0': {'type': 'unif', + 'lower': F0s-DeltaF0, + 'upper': F0s+DeltaF0}, + 'F1': {'type': 'unif', + 'lower': F1s-DeltaF1, + 'upper': F1s+DeltaF1}, + 'F2': F2s, + 'delta_F0': {'type': 'unif', + 'lower': 0, + 'upper': 1e-5}, + 'delta_F1': {'type': 'unif', + 'lower': -1e-11, + 'upper': 1e-11}, + 'tglitch': {'type': 'unif', + 'lower': tstart+0.1*Tspan, + 'upper': tstart+0.9*Tspan}, + 'Alpha': Alpha, + 'Delta': Delta, + } + +search = pyfstat.MCMCGlitchSearch( + label=label, outdir=outdir, sftfilepath=sftfilepath, + theta_prior=theta_prior, nglitch=1, tref=tref, nsteps=[500, 500], + ntemps=3, log10temperature_min=-0.5, minStartTime=tstart, + maxStartTime=tstart+Tspan) +search.run() +search.plot_corner(label_offset=0.8, add_prior=True) +search.print_summary() diff --git a/examples/glitch_robust_search_make_simulated_data.py b/examples/glitch_robust_search_make_simulated_data.py new file mode 100644 index 0000000000000000000000000000000000000000..106eddaea2821d34c433490f53f69597261c64e6 --- /dev/null +++ b/examples/glitch_robust_search_make_simulated_data.py @@ -0,0 +1,40 @@ +import numpy as np +import pyfstat + +outdir = 'data' +label = 'simulated_glitching_signal' + +# Properties of the GW data +tstart = 1000000000 +Tspan = 60 * 86400 + +tref = tstart + .5 * Tspan + +# Fixed properties of the signal +F0s = 30 +F1s = -1e-8 +F2s = 0 +Alpha = np.radians(83.6292) +Delta = np.radians(22.0144) +h0 = 1e-25 +sqrtSX = 1e-24 +psi = -0.1 +phi = 0 +cosi = 0.5 + +# Glitch properties +dtglitch = 0.45 * Tspan # time (in secs) after minStartTime +dF0 = 5e-6 +dF1 = 1e-12 + + +detectors = 'H1' + +glitch_data = pyfstat.Writer( + label=label, outdir=outdir, tref=tref, tstart=tstart, + F0=F0s, F1=F1s, F2=F2s, duration=Tspan, Alpha=Alpha, + Delta=Delta, sqrtSX=sqrtSX, dtglitch=dtglitch, + h0=h0, cosi=cosi, phi=phi, psi=psi, + delta_F0=dF0, delta_F1=dF1, add_noise=True) + +glitch_data.make_data() diff --git a/examples/transient_search_using_MCMC.py b/examples/transient_search_using_MCMC.py index c5c54c3c2fc318318e7b11c7b1c14025b6ded97f..12125f0c237cbbd3f6bc6e58539a741f737aab22 100644 --- a/examples/transient_search_using_MCMC.py +++ b/examples/transient_search_using_MCMC.py @@ -1,64 +1,20 @@ -import pyfstat -import numpy as np -import matplotlib.pyplot as plt +#!/usr/bin/env python -plt.style.use('thesis') +import pyfstat F0 = 30.0 F1 = -1e-10 F2 = 0 -Alpha = 5e-3 -Delta = 6e-2 - -tstart = 1000000000 -duration = 100*86400 -data_tstart = tstart - duration -data_tend = data_tstart + 3*duration -tref = .5*(data_tstart+data_tend) - -h0 = 1e-23 -sqrtSX = 1e-22 - -transient = pyfstat.Writer( - label='transient', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1, - F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, - minStartTime=data_tstart, maxStartTime=data_tend) -transient.make_data() -print transient.predict_fstat() - +Alpha = 0.5 +Delta = 1 +minStartTime = 1000000000 +maxStartTime = minStartTime + 200*86400 +Tspan = maxStartTime - minStartTime +tref = minStartTime DeltaF0 = 6e-7 DeltaF1 = 1e-13 -VF0 = (np.pi * duration * DeltaF0)**2 / 3.0 -VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45. -print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1) - -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 - } - -ntemps = 3 -log10temperature_min = -1 -nwalkers = 100 -nsteps = [750, 250] - -mcmc = pyfstat.MCMCSearch( - label='transient_search_initial_stage', outdir='data', - sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref, - minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps, - nwalkers=nwalkers, ntemps=ntemps, - log10temperature_min=log10temperature_min) -mcmc.run() -mcmc.plot_cumulative_max() -mcmc.print_summary() theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2., @@ -70,22 +26,24 @@ theta_prior = {'F0': {'type': 'unif', 'Alpha': Alpha, 'Delta': Delta, 'transient_tstart': {'type': 'unif', - 'lower': data_tstart, - 'upper': data_tend}, + 'lower': minStartTime, + 'upper': maxStartTime}, 'transient_duration': {'type': 'halfnorm', - 'loc': 0, - 'scale': 0.5*duration} + 'loc': 0.001*Tspan, + 'scale': 0.5*Tspan} } -nwalkers = 500 -nsteps = [200, 200] +ntemps = 2 +log10temperature_min = -1 +nwalkers = 100 +nsteps = [100, 100] mcmc = pyfstat.MCMCTransientSearch( label='transient_search', outdir='data', - sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref, - minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps, - nwalkers=nwalkers, ntemps=ntemps, + sftfilepath='data/*simulated_transient_signal*sft', + theta_prior=theta_prior, tref=tref, minStartTime=minStartTime, + maxStartTime=maxStartTime, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps, log10temperature_min=log10temperature_min) mcmc.run() -mcmc.plot_corner(add_prior=True) +mcmc.plot_corner(label_offset=0.7) mcmc.print_summary() diff --git a/examples/transient_search_using_MCMC_make_simulated_data.py b/examples/transient_search_using_MCMC_make_simulated_data.py new file mode 100644 index 0000000000000000000000000000000000000000..4fa3dd1d1157446b5de54dfb9b7f3b132889fdc1 --- /dev/null +++ b/examples/transient_search_using_MCMC_make_simulated_data.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +import pyfstat + +F0 = 30.0 +F1 = -1e-10 +F2 = 0 +Alpha = 0.5 +Delta = 1 + +minStartTime = 1000000000 +maxStartTime = minStartTime + 200*86400 + +transient_tstart = minStartTime + 50*86400 +transient_duration = 100*86400 +tref = minStartTime + +h0 = 1e-23 +sqrtSX = 1e-22 + +transient = pyfstat.Writer( + label='simulated_transient_signal', outdir='data', tref=tref, + tstart=transient_tstart, F0=F0, F1=F1, F2=F2, duration=transient_duration, + Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, minStartTime=minStartTime, + maxStartTime=maxStartTime) +transient.make_data()