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

Merge branch 'master' into develop-GA

parents d4b286d4 4870f7b0
......@@ -47,8 +47,8 @@ theta_prior = {'F0': {'type': 'unif',
'Delta': Delta
}
ntemps = 1
log10beta_min = -1
ntemps = 2
log10beta_min = -0.5
nwalkers = 100
nsteps = [300, 300]
......@@ -57,6 +57,9 @@ mcmc = pyfstat.MCMCSearch(
sftfilepattern='{}/*{}*sft'.format(outdir, label), theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
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.print_summary()
import pyfstat
import numpy as np
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
F1 = 1e-10
......@@ -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]
labels = ['$f - f_0$', '$\dot{f} - \dot{f}_0$', '$\ddot{f} - \ddot{f}_0$',
'$\widetilde{2\mathcal{F}}$']
fig, axes = projection_matrix(twoF, xyz, projection='log_mean', labels=labels,
whspace=0.1, factor=1.8)
fig, axes = gridcorner.gridcorner(
twoF, xyz, projection='log_mean', labels=labels, whspace=0.1, factor=1.8)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label))
......@@ -2,8 +2,6 @@ import pyfstat
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('paper')
F0 = 30.0
F1 = 0
F2 = 0
......
......@@ -56,7 +56,7 @@ mcmc = pyfstat.MCMCSearch(
sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
ntemps=ntemps, log10beta_min=log10beta_min)
mcmc.run(context='paper', subtractions=[30, -1e-10])
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
......
......@@ -2,8 +2,6 @@ import pyfstat
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('./paper-style.mplstyle')
F0 = 30.0
F1 = -1e-10
F2 = 0
......@@ -63,8 +61,8 @@ NstarMax = 1000
Nsegs0 = 100
fig, axes = plt.subplots(nrows=2, figsize=(3.4, 3.5))
fig, axes = mcmc.run(
NstarMax=NstarMax, Nsegs0=Nsegs0, subtractions=[F0, F1], labelpad=0.01,
plot_det_stat=False, return_fig=True, context='paper', fig=fig,
NstarMax=NstarMax, Nsegs0=Nsegs0, labelpad=0.01,
plot_det_stat=False, return_fig=True, fig=fig,
axes=axes)
for ax in axes:
ax.grid()
......
......@@ -58,6 +58,9 @@ mcmc = pyfstat.MCMCSemiCoherentSearch(
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
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.print_summary()
......@@ -25,9 +25,7 @@ theta_prior = {'F0': {'type': 'unif',
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'transient_tstart': {'type': 'unif',
'lower': minStartTime,
'upper': maxStartTime},
'transient_tstart': minStartTime,
'transient_duration': {'type': 'halfnorm',
'loc': 0.001*Tspan,
'scale': 0.5*Tspan}
......@@ -39,11 +37,12 @@ nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCTransientSearch(
label='transient_search', outdir='data',
sftfilepattern='data/*simulated_transient_signal*sft',
label='transient_search', outdir='data_l',
sftfilepattern='data_l/*simulated_transient_signal*sft',
theta_prior=theta_prior, tref=tref, minStartTime=minStartTime,
maxStartTime=maxStartTime, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
log10beta_min=log10beta_min,
transientWindowType='rect')
mcmc.run()
mcmc.plot_corner(label_offset=0.7)
mcmc.print_summary()
......@@ -19,8 +19,8 @@ h0 = 1e-23
sqrtSX = 1e-22
transient = pyfstat.Writer(
label='simulated_transient_signal', outdir='data', tref=tref,
label='simulated_transient_signal', outdir='data_l', 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)
maxStartTime=maxStartTime, transientWindowType='rect')
transient.make_data()
#!/usr/bin/env python
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 0.5
Delta = 1
minStartTime = 1000000000
maxStartTime = minStartTime + 2*86400
Tspan = maxStartTime - minStartTime
tref = minStartTime
Tsft = 1800
DeltaF0 = 1e-2
DeltaF1 = 1e-9
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,
'transient_tstart': {'type': 'unif',
'lower': minStartTime,
'upper': maxStartTime-2*Tsft},
'transient_duration': {'type': 'unif',
'lower': 2*Tsft,
'upper': Tspan-2*Tsft}
}
ntemps = 2
log10beta_min = -1
nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCTransientSearch(
label='transient_search', outdir='data_s',
sftfilepattern='data_s/*simulated_transient_signal*sft',
theta_prior=theta_prior, tref=tref, minStartTime=minStartTime,
maxStartTime=maxStartTime, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min,
transientWindowType='rect')
mcmc.run()
mcmc.plot_corner(label_offset=0.7)
mcmc.print_summary()
#!/usr/bin/env python
import pyfstat
import os
import numpy as np
import matplotlib.pyplot as plt
datadir = 'data_s'
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 0.5
Delta = 1
minStartTime = 1000000000
maxStartTime = minStartTime + 2*86400
Tspan = maxStartTime - minStartTime
tref = minStartTime
Tsft = 1800
m = 0.001
dF0 = np.sqrt(12*m)/(np.pi*Tspan)
DeltaF0 = 100*dF0
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
F1s = [F1]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
print('Standard CW search:')
search1 = pyfstat.GridSearch(
label='CW', outdir=datadir,
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
BSGL=False)
search1.run()
search1.print_max_twoF()
search1.plot_1D(xkey='F0',
xlabel='freq [Hz]', ylabel='$2\mathcal{F}$')
print('with t0,tau bands:')
search2 = pyfstat.GridSearch(
label='tCW', outdir=datadir,
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
transientWindowType='rect', t0Band=Tspan-2*Tsft, tauBand=Tspan,
BSGL=False)
search2.run()
search2.print_max_twoF()
search2.plot_1D(xkey='F0',
xlabel='freq [Hz]', ylabel='$2\mathcal{F}$')
#!/usr/bin/env python
import pyfstat
F0 = 30.0
F1 = -1e-10
F2 = 0
Alpha = 0.5
Delta = 1
minStartTime = 1000000000
maxStartTime = minStartTime + 2*86400
transient_tstart = minStartTime + 0.5*86400
transient_duration = 1*86400
tref = minStartTime
h0 = 1e-23
sqrtSX = 1e-22
detectors = 'H1,L1'
Tsft = 1800
transient = pyfstat.Writer(
label='simulated_transient_signal', outdir='data_s',
tref=tref, tstart=transient_tstart, duration=transient_duration,
F0=F0, F1=F1, F2=F2, Alpha=Alpha, Delta=Delta, h0=h0,
detectors=detectors,sqrtSX=sqrtSX,
minStartTime=minStartTime, maxStartTime=maxStartTime,
transientWindowType='rect', Tsft=Tsft)
transient.make_data()
......@@ -59,6 +59,6 @@ mcmc = pyfstat.MCMCSearch(
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.setup_initialisation(100, scatter_val=1e-10)
mcmc.run(subtractions=[F0, F1])
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
......@@ -330,7 +330,8 @@ class ComputeFstat(BaseSearchClass):
@helper_functions.initializer
def __init__(self, tref, sftfilepattern=None, minStartTime=None,
maxStartTime=None, binary=False, transient=True, BSGL=False,
maxStartTime=None, binary=False, BSGL=False,
transientWindowType=None, t0Band=None, tauBand=None,
detectors=None, minCoverFreq=None, maxCoverFreq=None,
injectSources=None, injectSqrtSX=None, assumeSqrtSX=None,
SSBprec=None):
......@@ -347,10 +348,18 @@ class ComputeFstat(BaseSearchClass):
this epoch
binary : bool
If true, search of binary parameters.
transient : bool
If true, allow for the Fstat to be computed over a transient range.
BSGL : bool
If true, compute the BSGL rather than the twoF value.
transientWindowType: str
If 'rect' or 'exp',
allow for the Fstat to be computed over a transient range.
('none' instead of None explicitly calls the transient-window
function, but with the full range, for debugging)
t0Band, tauBand: int
if >0, search t0 in (minStartTime,minStartTime+t0Band)
and tau in (2*Tsft,2*Tsft+tauBand).
if =0, only compute CW Fstat with t0=minStartTime,
tau=maxStartTime-minStartTime.
detectors : str
Two character reference to the data to use, specify None for no
contraint. If multiple-separate by comma.
......@@ -477,7 +486,7 @@ class ComputeFstat(BaseSearchClass):
logging.info('Initialising FstatInput')
dFreq = 0
if self.transient:
if self.transientWindowType:
self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET
else:
self.whatToCompute = lalpulsar.FSTATQ_2F
......@@ -593,14 +602,41 @@ class ComputeFstat(BaseSearchClass):
self.whatToCompute = (self.whatToCompute +
lalpulsar.FSTATQ_2F_PER_DET)
if self.transient:
if self.transientWindowType:
logging.info('Initialising transient parameters')
self.windowRange = lalpulsar.transientWindowRange_t()
self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR
self.windowRange.t0Band = 0
self.windowRange.dt0 = 1
self.windowRange.tauBand = 0
self.windowRange.dtau = 1
transientWindowTypes = {'none': lalpulsar.TRANSIENT_NONE,
'rect': lalpulsar.TRANSIENT_RECTANGULAR,
'exp': lalpulsar.TRANSIENT_EXPONENTIAL}
if self.transientWindowType in transientWindowTypes:
self.windowRange.type = transientWindowTypes[self.transientWindowType]
else:
raise ValueError(
'Unknown window-type ({}) passed as input, [{}] allows.'
.format(self.transientWindowType,
', '.join(transientWindowTypes)))
self.Tsft = int(1.0/SFTCatalog.data[0].header.deltaF)
if self.t0Band is None:
self.windowRange.t0Band = 0
self.windowRange.dt0 = 1
else:
if not isinstance(self.t0Band, int):
logging.warn('Casting non-integer t0Band={} to int...'
.format(self.t0Band))
self.t0Band = int(self.t0Band)
self.windowRange.t0Band = self.t0Band
self.windowRange.dt0 = self.Tsft
if self.tauBand is None:
self.windowRange.tauBand = 0
self.windowRange.dtau = 1
else:
if not isinstance(self.tauBand, int):
logging.warn('Casting non-integer tauBand={} to int...'
.format(self.tauBand))
self.tauBand = int(self.tauBand)
self.windowRange.tauBand = self.tauBand
self.windowRange.dtau = self.Tsft
def get_fullycoherent_twoF(self, tstart, tend, F0, F1, F2, Alpha, Delta,
asini=None, period=None, ecc=None, tp=None,
......@@ -624,7 +660,7 @@ class ComputeFstat(BaseSearchClass):
self.whatToCompute
)
if self.transient is False:
if not self.transientWindowType:
if self.BSGL is False:
return self.FstatResults.twoF[0]
......@@ -636,13 +672,20 @@ class ComputeFstat(BaseSearchClass):
return log10_BSGL/np.log10(np.exp(1))
self.windowRange.t0 = int(tstart) # TYPE UINT4
self.windowRange.tau = int(tend - tstart) # TYPE UINT4
if self.windowRange.tauBand == 0:
# true single-template search also in transient params:
# actual (t0,tau) window was set with tstart, tend before
self.windowRange.tau = int(tend - tstart) # TYPE UINT4
else:
# grid search: start at minimum tau required for nondegenerate
# F-stat computation
self.windowRange.tau = int(2*self.Tsft)
FS = lalpulsar.ComputeTransientFstatMap(
self.FstatResults.multiFatoms[0], self.windowRange, False)
twoF = 2*np.max(FS.F_mn.data)
if self.BSGL is False:
twoF = 2*FS.F_mn.data[0][0]
if np.isnan(twoF):
return 0
else:
......@@ -657,10 +700,17 @@ class ComputeFstat(BaseSearchClass):
FS1 = lalpulsar.ComputeTransientFstatMap(
FstatResults_single.multiFatoms[0], self.windowRange, False)
self.twoFX[0] = 2*FS0.F_mn.data[0][0]
self.twoFX[1] = 2*FS1.F_mn.data[0][0]
# for now, use the Doppler parameter with
# multi-detector F maximised over t0,tau
# to return BSGL
# FIXME: should we instead compute BSGL over the whole F_mn
# and return the maximum of that?
idx_maxTwoF = np.argmax(FS.F_mn.data)
self.twoFX[0] = 2*FS0.F_mn.data[idx_maxTwoF]
self.twoFX[1] = 2*FS1.F_mn.data[idx_maxTwoF]
log10_BSGL = lalpulsar.ComputeBSGL(
2*FS.F_mn.data[0][0], self.twoFX, self.BSGLSetup)
twoF, self.twoFX, self.BSGLSetup)
return log10_BSGL/np.log10(np.exp(1))
......@@ -696,8 +746,9 @@ class ComputeFstat(BaseSearchClass):
max_tau = SFTmaxStartTime - tstart
taus = np.linspace(min_tau, max_tau, npoints)
twoFs = []
if self.transient is False:
self.transient = True
if not self.transientWindowType:
# still call the transient-Fstat-map function, but using the full range
self.transientWindowType = 'none'
self.init_computefstatistic_single_point()
for tau in taus:
twoFs.append(self.get_fullycoherent_twoF(
......@@ -868,7 +919,9 @@ class SemiCoherentSearch(ComputeFstat):
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
self.set_ephemeris_files()
self.transient = True
self.transientWindowType = 'rect'
self.t0Band = None
self.tauBand = None
self.init_computefstatistic_single_point()
self.init_semicoherent_parameters()
......@@ -876,7 +929,7 @@ class SemiCoherentSearch(ComputeFstat):
logging.info(('Initialising semicoherent parameters from {} to {} in'
' {} segments').format(
self.minStartTime, self.maxStartTime, self.nsegs))
self.transient = True
self.transientWindowType = 'rect'
self.whatToCompute = lalpulsar.FSTATQ_2F+lalpulsar.FSTATQ_ATOMS_PER_DET
self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime,
self.nsegs+1)
......@@ -915,7 +968,7 @@ class SemiCoherentSearch(ComputeFstat):
self.whatToCompute
)
#if self.transient is False:
#if not self.transientWindowType:
# if self.BSGL is False:
# return self.FstatResults.twoF[0]
# twoF = np.float(self.FstatResults.twoF[0])
......@@ -1005,8 +1058,10 @@ class SemiCoherentGlitchSearch(ComputeFstat):
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
self.set_ephemeris_files()
self.transient = True
self.binary = False
self.transientWindowType = 'rect'
self.t0Band = None
self.tauBand = None
self.binary = False
self.init_computefstatistic_single_point()
def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
......
......@@ -31,7 +31,8 @@ class GridSearch(BaseSearchClass):
Deltas, tref=None, minStartTime=None, maxStartTime=None,
nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detectors=None, SSBprec=None, injectSources=None,
input_arrays=False, assumeSqrtSX=None):
input_arrays=False, assumeSqrtSX=None,
transientWindowType=None, t0Band=None, tauBand=None):
"""
Parameters
----------
......@@ -48,6 +49,16 @@ class GridSearch(BaseSearchClass):
GPS seconds of the reference time, start time and end time
input_arrays: bool
if true, use the F0s, F1s, etc as is
transientWindowType: str
If 'rect' or 'exp', compute atoms so that a transient (t0,tau) map
can later be computed. ('none' instead of None explicitly calls
the transient-window function, but with the full range, for
debugging). Currently only supported for nsegs=1.
t0Band, tauBand: int
if >0, search t0 in (minStartTime,minStartTime+t0Band)
and tau in (2*Tsft,2*Tsft+tauBand).
if =0, only compute CW Fstat with t0=minStartTime,
tau=maxStartTime-minStartTime.
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
......@@ -66,7 +77,9 @@ class GridSearch(BaseSearchClass):
self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=False,
detectors=self.detectors,
transientWindowType=self.transientWindowType,
t0Band=self.t0Band, tauBand=self.tauBand,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources,
......@@ -573,7 +586,7 @@ class FrequencySlidingWindow(GridSearch):
self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=True,
detectors=self.detectors, transientWindowType=self.transientWindowType,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources)
......@@ -779,10 +792,15 @@ class EarthTest(GridSearch):
r'$\Delta P_\mathrm{spin}$ [min]',
r'$2\mathcal{F}$']
from projection_matrix import projection_matrix
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")
fig, axes = projection_matrix(data, xyz, projection=projection,
factor=1.6, labels=labels)
fig, axes = gridcorner(data, xyz, projection=projection, factor=1.6,
labels=labels)
axes[-1][-1].axvline((lal.DAYJUL_SI - lal.DAYSID_SI)/60.0, color='C3')
plt.suptitle(
'T={:.1f} days, $f$={:.2f} Hz, $\log\mathcal{{B}}_{{S/A}}$={:.1f},'
......
......@@ -25,7 +25,8 @@ class Writer(BaseSearchClass):
tref=None, F0=30, F1=1e-10, F2=0, Alpha=5e-3,
Delta=6e-2, h0=0.1, cosi=0.0, psi=0.0, phi=0, Tsft=1800,
outdir=".", sqrtSX=1, Band=4, detectors='H1',
minStartTime=None, maxStartTime=None, add_noise=True):
minStartTime=None, maxStartTime=None, add_noise=True,
transientWindowType='none'):
"""
Parameters
----------
......@@ -91,9 +92,9 @@ class Writer(BaseSearchClass):
self.make_cff()
self.run_makefakedata()
def get_single_config_line(self, i, Alpha, Delta, h0, cosi, psi, phi, F0,
F1, F2, tref, tstart, duration_days):
template = (
def get_base_template(self, i, Alpha, Delta, h0, cosi, psi, phi, F0,
F1, F2, tref):
return (
"""[TS{}]
Alpha = {:1.18e}
Delta = {:1.18e}
......@@ -104,12 +105,35 @@ phi0 = {:1.18e}
Freq = {:1.18e}
f1dot = {:1.18e}
f2dot = {:1.18e}
refTime = {:10.6f}
transientWindowType=rect
transientStartTime={:10.3f}
transientTauDays={:1.3f}\n""")
refTime = {:10.6f}""")
def get_single_config_line_cw(
self, i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref):
template = (self.get_base_template(
i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref) + """\n""")
return template.format(
i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref)
def get_single_config_line_tcw(
self, i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref,
window, tstart, duration_days):
template = (self.get_base_template(
i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, F2, tref) + """
transientWindowType = {:s}
transientStartTime = {:10.3f}
transientTauDays = {:1.3f}\n""")
return template.format(i, Alpha, Delta, h0, cosi, psi, phi, F0, F1,
F2, tref, tstart, duration_days)