Commit f7311ed6 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Merge branch 'develop-GA' of gitlab.aei.uni-hannover.de:GregAshton/PyFstat into develop-GA

parents adb2e3d2 f93e27b1
test_app: test_app:
script: script:
- . /home/greg/lalsuite-install/etc/lalapps-user-env.sh - . /home/greg/lalsuite-install/etc/lalapps-user-env.sh
- /home/greg/anaconda2/bin/python tests.py TestWriter - /home/greg/anaconda2/bin/python tests.py Writer
- /home/greg/anaconda2/bin/python tests.py TestBaseSearchClass - /home/greg/anaconda2/bin/python tests.py par
- /home/greg/anaconda2/bin/python tests.py TestComputeFstat - /home/greg/anaconda2/bin/python tests.py BaseSearchClass
- /home/greg/anaconda2/bin/python tests.py TestSemiCoherentGlitchSearch - /home/greg/anaconda2/bin/python tests.py ComputeFstat
- /home/greg/anaconda2/bin/python tests.py SemiCoherentSearch
- /home/greg/anaconda2/bin/python tests.py SemiCoherentGlitchSearch
- /home/greg/anaconda2/bin/python tests.py MCMCSearch
- /home/greg/anaconda2/bin/python tests.py GridSearch
...@@ -47,8 +47,8 @@ theta_prior = {'F0': {'type': 'unif', ...@@ -47,8 +47,8 @@ theta_prior = {'F0': {'type': 'unif',
'Delta': Delta 'Delta': Delta
} }
ntemps = 1 ntemps = 2
log10beta_min = -1 log10beta_min = -0.5
nwalkers = 100 nwalkers = 100
nsteps = [300, 300] nsteps = [300, 300]
...@@ -57,6 +57,9 @@ mcmc = pyfstat.MCMCSearch( ...@@ -57,6 +57,9 @@ mcmc = pyfstat.MCMCSearch(
sftfilepattern='{}/*{}*sft'.format(outdir, label), theta_prior=theta_prior, sftfilepattern='{}/*{}*sft'.format(outdir, label), theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
nwalkers=nwalkers, ntemps=ntemps, log10beta_min=log10beta_min) nwalkers=nwalkers, ntemps=ntemps, log10beta_min=log10beta_min)
mcmc.run(subtractions=[F0, F1]) mcmc.transform_dictionary = dict(
F0=dict(subtractor=F0, symbol='$f-f^\mathrm{s}$'),
F1=dict(subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$'))
mcmc.run()
mcmc.plot_corner(add_prior=True) mcmc.plot_corner(add_prior=True)
mcmc.print_summary() mcmc.print_summary()
import numpy as np
from pyfstat import MCMCSearch
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},
'F1': {'type': 'unif', 'lower': F1*(1+1e-3), 'upper': F1*(1-1e-3)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
}
ntemps = 2
log10beta_min = -0.01
nwalkers = 100
nsteps = [500, 500]
mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',
sftfilepattern='data/*_glitch*.sft',
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
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
sftfilepattern = '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, sftfilepattern=sftfilepattern,
theta_prior=theta_prior, nglitch=1, tref=tref, nsteps=[500, 500],
ntemps=3, log10beta_min=-0.5, minStartTime=tstart,
maxStartTime=tstart+Tspan)
search.run()
search.plot_corner(label_offset=0.8, add_prior=True)
search.print_summary()
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()
from pyfstat import Writer, GlitchWriter from pyfstat import Writer, GlitchWriter
import numpy as np import numpy as np
outdir = 'data'
# First, we generate data with a reasonably strong smooth signal # First, we generate data with a reasonably strong smooth signal
# Define parameters of the Crab pulsar as an example # Define parameters of the Crab pulsar as an example
...@@ -9,19 +10,19 @@ F1 = -1e-10 ...@@ -9,19 +10,19 @@ F1 = -1e-10
F2 = 0 F2 = 0
Alpha = np.radians(83.6292) Alpha = np.radians(83.6292)
Delta = np.radians(22.0144) Delta = np.radians(22.0144)
tref = 362750407.0
# Signal strength # Signal strength
h0 = 1e-23 h0 = 5e-24
# Properties of the GW data # Properties of the GW data
sqrtSX = 1e-22 sqrtSX = 1e-22
tstart = 1000000000 tstart = 1000000000
duration = 100*86400 duration = 100*86400
tend = tstart+duration tend = tstart+duration
tref = tstart + 0.5*duration
data = Writer( data = Writer(
label='basic', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1, label='0_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX) F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
data.make_data() data.make_data()
...@@ -31,11 +32,11 @@ print 'Predicted twoF value: {}\n'.format(twoF) ...@@ -31,11 +32,11 @@ print 'Predicted twoF value: {}\n'.format(twoF)
# Next, taking the same signal parameters, we include a glitch half way through # Next, taking the same signal parameters, we include a glitch half way through
dtglitch = duration/2.0 dtglitch = duration/2.0
delta_F0 = 4e-5 delta_F0 = 5e-6
delta_F1 = 0 delta_F1 = 0
glitch_data = GlitchWriter( glitch_data = GlitchWriter(
label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1, label='1_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1) dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
glitch_data.make_data() glitch_data.make_data()
...@@ -49,7 +50,7 @@ delta_F1 = [0, 0] ...@@ -49,7 +50,7 @@ delta_F1 = [0, 0]
delta_F2 = [0, 0] delta_F2 = [0, 0]
two_glitch_data = GlitchWriter( two_glitch_data = GlitchWriter(
label='twoglitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1, label='2_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_phi=delta_phi, delta_F0=delta_F0, dtglitch=dtglitch, delta_phi=delta_phi, delta_F0=delta_F0,
delta_F1=delta_F1, delta_F2=delta_F2) delta_F1=delta_F1, delta_F2=delta_F2)
......
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'delta_F0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-5*F0},
'delta_F1': 0,
'tglitch': {'type': 'unif',
'lower': tstart+0.1*duration,
'upper': tstart+0.9*duration},
}
ntemps = 4
log10beta_min = -1
nwalkers = 100
nsteps = [5000, 1000, 1000]
mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_glitch_search_using_MCMC', 'data',
sftfilepattern='data/*_glitch*sft', theta_prior=theta_prior, tref=tref,
tstart=tstart, tend=tend, nsteps=nsteps, nwalkers=nwalkers,
scatter_val=1e-10, nglitch=1, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 5e-3
Delta = 6e-2
tref = 362750407.0
tstart = 1000000000
duration = 100*86400
tend = tstart + duration
theta_prior = {'F0': {'type': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)},
'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'delta_F0_0': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0},
'delta_F0_1': {'type': 'halfnorm', 'loc': 0,
'scale': 1e-7*F0},
'delta_F1_0': 0,
'delta_F1_1': 0,
'tglitch_0': {'type': 'unif',
'lower': tstart+0.01*duration,
'upper': tstart+0.5*duration},
'tglitch_1': {'type': 'unif',
'lower': tstart+0.5*duration,
'upper': tstart+0.99*duration},
}
nwalkers = 100
nsteps = [1000, 1000, 5000]
mcmc = pyfstat.MCMCGlitchSearch(
'semi_coherent_twoglitch_search', 'data', sftfilepattern='data/*twoglitch*sft',
theta_prior=theta_prior, tref=tref, tstart=tstart,
tend=tend, nsteps=nsteps, nwalkers=nwalkers, scatter_val=1e-10, nglitch=2)
mcmc.run()
mcmc.plot_corner(add_prior=True, tglitch_ratio=True, figsize=(10, 10))
mcmc.print_summary()
import numpy as np
import matplotlib.pyplot as plt
import pyfstat
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
plt.style.use('paper')
label = 'semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch'
Nstar = 100
F0_width = np.sqrt(Nstar)*np.sqrt(12)/(np.pi*duration)
F1_width = np.sqrt(Nstar)*np.sqrt(180)/(np.pi*duration**2)
theta_prior = {
'F0': {'type': 'unif',
'lower': F0-F0_width/2.,
'upper': F0+F0_width/2.},
'F1': {'type': 'unif',
'lower': F1-F1_width/2.,
'upper': F1+F1_width/2.},
'F2': F2,
'delta_F0': {'type': 'unif',
'lower': 0,
'upper': 1e-5},
'delta_F1': 0,
'tglitch': {'type': 'unif',
'lower': tstart+0.1*duration,
'upper': tstart+0.9*duration},
'Alpha': Alpha,
'Delta': Delta,
}
ntemps = 2
log10beta_min = -0.5
nwalkers = 100
nsteps = [500, 500]
mcmc = pyfstat.MCMCGlitchSearch(
label=label, sftfilepattern='data/*1_glitch*sft', theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tstart+duration,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min, nglitch=1)
mcmc.transform_dictionary['F0'] = dict(
subtractor=F0, symbol='$f-f^\mathrm{s}$')
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$')
mcmc.run()
mcmc.plot_corner()
mcmc.print_summary()
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
try:
from gridcorner import gridcorner
except ImportError:
raise ImportError(
"Python module 'gridcorner' not found, please install from "
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner")
label = 'semicoherent_glitch_robust_directed_grid_search_on_1_glitch'
plt.style.use('paper')
Nstar = 150
F0_width = np.sqrt(Nstar)*np.sqrt(12)/(np.pi*duration)
F1_width = np.sqrt(Nstar)*np.sqrt(180)/(np.pi*duration**2)
N = 61
F0s = [F0-F0_width/2., F0+F0_width/2., F0_width/N]
F1s = [F1-F1_width/2., F1+F1_width/2., F1_width/N]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
max_delta_F0 = 1e-5
tglitchs = [tstart+0.1*duration, tstart+0.9*duration, 0.8*float(duration)/N]
delta_F0s = [0, max_delta_F0, max_delta_F0/N]
delta_F1s = [0]
search = pyfstat.GridGlitchSearch(
label, outdir, 'data/*1_glitch*sft', F0s=F0s, F1s=F1s, F2s=F2s,
Alphas=Alphas, Deltas=Deltas, tref=tref, minStartTime=tstart,
maxStartTime=tstart+duration, tglitchs=tglitchs, delta_F0s=delta_F0s,
delta_F1s=delta_F1s)
search.run()
F0_vals = np.unique(search.data[:, 0]) - F0
F1_vals = np.unique(search.data[:, 1]) - F1
delta_F0s_vals = np.unique(search.data[:, 5])
tglitch_vals = np.unique(search.data[:, 7])
tglitch_vals_days = (tglitch_vals-tstart) / 86400.
twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals),
len(delta_F0s_vals), len(tglitch_vals)))
xyz = [F0_vals, F1_vals, delta_F0s_vals, tglitch_vals_days]
labels = ['$f - f_0$\n[Hz]', '$\dot{f} - \dot{f}_0$\n[Hz/s]',
'$\delta f$\n[Hz]', '$t^g_0$\n[days]', '$\widehat{2\mathcal{F}}$']
fig, axes = gridcorner(
twoF, xyz, projection='log_mean', whspace=0.1, factor=1.2, labels=labels)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label),
bbox_inches='tight')
import numpy as np
import matplotlib.pyplot as plt
import pyfstat
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
plt.style.use('paper')
label = 'standard_directed_MCMC_search_on_1_glitch'
Nstar = 10000
F0_width = np.sqrt(Nstar)*np.sqrt(12)/(np.pi*duration)
F1_width = np.sqrt(Nstar)*np.sqrt(180)/(np.pi*duration**2)
theta_prior = {
'F0': {'type': 'unif',
'lower': F0-F0_width/2.,
'upper': F0+F0_width/2.},
'F1': {'type': 'unif',
'lower': F1-F1_width/2.,
'upper': F1+F1_width/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
}
ntemps = 2
log10beta_min = -0.5
nwalkers = 100
nsteps = [500, 2000]
mcmc = pyfstat.MCMCSearch(
label=label, sftfilepattern='data/*1_glitch*sft', theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tstart+duration,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.transform_dictionary['F0'] = dict(
subtractor=F0, symbol='$f-f^\mathrm{s}$')
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$')
mcmc.run()
mcmc.plot_corner()
mcmc.print_summary()
import pyfstat import pyfstat
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from projection_matrix import projection_matrix
plt.style.use('paper') try:
from gridcorner import gridcorner
except ImportError:
raise ImportError(
"Python module 'gridcorner' not found, please install from "
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner")
F0 = 30.0 F0 = 30.0
F1 = 1e-10 F1 = 1e-10
...@@ -55,6 +59,6 @@ twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals), len(F2_vals))) ...@@ -55,6 +59,6 @@ twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals), len(F2_vals)))
xyz = [F0_vals, F1_vals, F2_vals] xyz = [F0_vals, F1_vals, F2_vals]
labels = ['$f - f_0$', '$\dot{f} - \dot{f}_0$', '$\ddot{f} - \ddot{f}_0$', labels = ['$f - f_0$', '$\dot{f} - \dot{f}_0$', '$\ddot{f} - \ddot{f}_0$',
'$\widetilde{2\mathcal{F}}$'] '$\widetilde{2\mathcal{F}}$']
fig, axes = projection_matrix(twoF, xyz, projection='log_mean', labels=labels, fig, axes = gridcorner(
whspace=0.1, factor=1.8) twoF, xyz, projection='log_mean', labels=labels, whspace=0.1, factor=1.8)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label)) fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label))
...@@ -2,8 +2,6 @@ import pyfstat ...@@ -2,8 +2,6 @@ import pyfstat
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.style.use('paper')
F0 = 30.0 F0 = 30.0
F1 = 0 F1 = 0
F2 = 0 F2 = 0
......
...@@ -56,7 +56,7 @@ mcmc = pyfstat.MCMCSearch( ...@@ -56,7 +56,7 @@ mcmc = pyfstat.MCMCSearch(
sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref, sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
ntemps=ntemps, log10beta_min=log10beta_min) ntemps=ntemps, log10beta_min=log10beta_min)
mcmc.run(context='paper', subtractions=[30, -1e-10]) mcmc.run()
mcmc.plot_corner(add_prior=True) mcmc.plot_corner(add_prior=True)
mcmc.print_summary() mcmc.print_summary()
......
...@@ -2,8 +2,6 @@ import pyfstat ...@@ -2,8 +2,6 @@ import pyfstat
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.style.use('./paper-style.mplstyle')
F0 = 30.0 F0 = 30.0