""" Searches using grid-based methods """
from __future__ import division, absolute_import, print_function
import os
import logging
import itertools
from collections import OrderedDict
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.misc import logsumexp
import pyfstat.helper_functions as helper_functions
from pyfstat.core import (BaseSearchClass, ComputeFstat,
SemiCoherentGlitchSearch, SemiCoherentSearch, tqdm,
args, read_par)
import lalpulsar
import lal
class GridSearch(BaseSearchClass):
""" Gridded search using ComputeFstat """
tex_labels = {'F0': '$f$', 'F1': '$\dot{f}$', 'F2': '$\ddot{f}$',
'Alpha': r'$\alpha$', 'Delta': r'$\delta$'}
tex_labels0 = {'F0': '$-f_0$', 'F1': '$-\dot{f}_0$', 'F2': '$-\ddot{f}_0$',
'Alpha': r'$-\alpha_0$', 'Delta': r'$-\delta_0$'}
@helper_functions.initializer
def __init__(self, label, outdir, sftfilepattern, F0s, F1s, F2s, Alphas,
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):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
F0s, F1s, F2s, delta_F0s, delta_F1s, tglitchs, Alphas, Deltas: tuple
Length 3 tuple describing the grid for each parameter, e.g
[F0min, F0max, dF0], for a fixed value simply give [F0]. Unless
input_arrays == True, then these are the values to search at.
tref, minStartTime, maxStartTime: int
GPS seconds of the reference time, start time and end time
input_arrays: bool
if true, use the F0s, F1s, etc as is
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
self.set_out_file()
self.keys = ['_', '_', 'F0', 'F1', 'F2', 'Alpha', 'Delta']
self.search_keys = [x+'s' for x in self.keys[2:]]
for k in self.search_keys:
setattr(self, k, np.atleast_1d(getattr(self, k)))
def inititate_search_object(self):
logging.info('Setting up search object')
if self.nsegs == 1:
self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=False,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources,
assumeSqrtSX=self.assumeSqrtSX)
self.search.get_det_stat = self.search.get_fullycoherent_twoF
else:
self.search = SemiCoherentSearch(
label=self.label, outdir=self.outdir, tref=self.tref,
nsegs=self.nsegs, sftfilepattern=self.sftfilepattern,
BSGL=self.BSGL, minStartTime=self.minStartTime,
maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, detectors=self.detectors,
injectSources=self.injectSources)
def cut_out_tstart_tend(*vals):
return self.search.get_semicoherent_twoF(*vals[2:])
self.search.get_det_stat = cut_out_tstart_tend
def get_array_from_tuple(self, x):
if len(x) == 1:
return np.array(x)
elif len(x) == 3 and self.input_arrays is False:
return np.arange(x[0], x[1], x[2])
else:
logging.info('Using tuple as is')
return np.array(x)
def get_input_data_array(self):
logging.info("Generating input data array")
coord_arrays = []
for tup in ([self.minStartTime], [self.maxStartTime], self.F0s,
self.F1s, self.F2s, self.Alphas, self.Deltas):
coord_arrays.append(self.get_array_from_tuple(tup))
input_data = []
for vals in itertools.product(*coord_arrays):
input_data.append(vals)
self.input_data = np.array(input_data)
self.coord_arrays = coord_arrays
def check_old_data_is_okay_to_use(self):
if args.clean:
return False
if os.path.isfile(self.out_file) is False:
logging.info('No old data found, continuing with grid search')
return False
if self.sftfilepattern is not None:
oldest_sft = min([os.path.getmtime(f) for f in
self._get_list_of_matching_sfts()])
if os.path.getmtime(self.out_file) < oldest_sft:
logging.info('Search output data outdates sft files,'
+ ' continuing with grid search')
return False
data = np.atleast_2d(np.genfromtxt(self.out_file, delimiter=' '))
if np.all(data[:, 0:-1] == self.input_data):
logging.info(
'Old data found with matching input, no search performed')
return data
else:
logging.info(
'Old data found, input differs, continuing with grid search')
return False
return False
def run(self, return_data=False):
self.get_input_data_array()
old_data = self.check_old_data_is_okay_to_use()
if old_data is not False:
self.data = old_data
return
if hasattr(self, 'search') is False:
self.inititate_search_object()
data = []
for vals in tqdm(self.input_data):
FS = self.search.get_det_stat(*vals)
data.append(list(vals) + [FS])
data = np.array(data, dtype=np.float)
if return_data:
return data
else:
logging.info('Saving data to {}'.format(self.out_file))
np.savetxt(self.out_file, data, delimiter=' ')
self.data = data
def convert_F0_to_mismatch(self, F0, F0hat, Tseg):
DeltaF0 = F0[1] - F0[0]
m_spacing = (np.pi*Tseg*DeltaF0)**2 / 12.
N = len(F0)
return np.arange(-N*m_spacing/2., N*m_spacing/2., m_spacing)
def convert_F1_to_mismatch(self, F1, F1hat, Tseg):
DeltaF1 = F1[1] - F1[0]
m_spacing = (np.pi*Tseg**2*DeltaF1)**2 / 720.
N = len(F1)
return np.arange(-N*m_spacing/2., N*m_spacing/2., m_spacing)
def add_mismatch_to_ax(self, ax, x, y, xkey, ykey, xhat, yhat, Tseg):
axX = ax.twiny()
axX.zorder = -10
axY = ax.twinx()
axY.zorder = -10
if xkey == 'F0':
m = self.convert_F0_to_mismatch(x, xhat, Tseg)
axX.set_xlim(m[0], m[-1])
if ykey == 'F1':
m = self.convert_F1_to_mismatch(y, yhat, Tseg)
axY.set_ylim(m[0], m[-1])
def plot_1D(self, xkey, ax=None, x0=None, xrescale=1, savefig=True,
xlabel=None, ylabel='$\widetilde{2\mathcal{F}}$'):
if ax is None:
fig, ax = plt.subplots()
xidx = self.keys.index(xkey)
x = np.unique(self.data[:, xidx])
if x0:
x = x - x0
x = x * xrescale
z = self.data[:, -1]
ax.plot(x, z)
if x0:
ax.set_xlabel(self.tex_labels[xkey]+self.tex_labels0[xkey])
else:
ax.set_xlabel(self.tex_labels[xkey])
if xlabel:
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
if savefig:
fig.tight_layout()
fig.savefig('{}/{}_1D.png'.format(self.outdir, self.label))
else:
return fig, ax
def plot_2D(self, xkey, ykey, ax=None, save=True, vmin=None, vmax=None,
add_mismatch=None, xN=None, yN=None, flat_keys=[],
rel_flat_idxs=[], flatten_method=np.max, title=None,
predicted_twoF=None, cm=None, cbarkwargs={}, x0=None, y0=None,
colorbar=False):
""" Plots a 2D grid of 2F values
Parameters
----------
add_mismatch: tuple (xhat, yhat, Tseg)
If not None, add a secondary axis with the metric mismatch from the
point xhat, yhat with duration Tseg
flatten_method: np.max
Function to use in flattening the flat_keys
"""
if ax is None:
fig, ax = plt.subplots()
xidx = self.keys.index(xkey)
yidx = self.keys.index(ykey)
flat_idxs = [self.keys.index(k) for k in flat_keys]
x = np.unique(self.data[:, xidx])
if x0:
x = x-x0
y = np.unique(self.data[:, yidx])
if y0:
y = y-y0
flat_vals = [np.unique(self.data[:, j]) for j in flat_idxs]
z = self.data[:, -1]
Y, X = np.meshgrid(y, x)
shape = [len(x), len(y)] + [len(v) for v in flat_vals]
Z = z.reshape(shape)
if len(rel_flat_idxs) > 0:
Z = flatten_method(Z, axis=tuple(rel_flat_idxs))
if predicted_twoF:
Z = (predicted_twoF - Z) / (predicted_twoF + 4)
if cm is None:
cm = plt.cm.viridis_r
else:
if cm is None:
cm = plt.cm.viridis
pax = ax.pcolormesh(X, Y, Z, cmap=cm, vmin=vmin, vmax=vmax)
if colorbar:
cb = plt.colorbar(pax, ax=ax, **cbarkwargs)
cb.set_label('$2\mathcal{F}$')
if add_mismatch:
self.add_mismatch_to_ax(ax, x, y, xkey, ykey, *add_mismatch)
ax.set_xlim(x[0], x[-1])
ax.set_ylim(y[0], y[-1])
if x0:
ax.set_xlabel(self.tex_labels[xkey]+self.tex_labels0[xkey])
else:
ax.set_xlabel(self.tex_labels[xkey])
if y0:
ax.set_ylabel(self.tex_labels[ykey]+self.tex_labels0[ykey])
else:
ax.set_ylabel(self.tex_labels[ykey])
if title:
ax.set_title(title)
if xN:
ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(xN))
if yN:
ax.yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(yN))
if save:
fig.tight_layout()
fig.savefig('{}/{}_2D.png'.format(self.outdir, self.label))
else:
return ax
def get_max_twoF(self):
twoF = self.data[:, -1]
idx = np.argmax(twoF)
v = self.data[idx, :]
d = OrderedDict(minStartTime=v[0], maxStartTime=v[1], F0=v[2], F1=v[3],
F2=v[4], Alpha=v[5], Delta=v[6], twoF=v[7])
return d
def print_max_twoF(self):
d = self.get_max_twoF()
print('Max twoF values for {}:'.format(self.label))
for k, v in d.iteritems():
print(' {}={}'.format(k, v))
def set_out_file(self, extra_label=None):
if self.detectors:
dets = self.detectors.replace(',', '')
else:
dets = 'NA'
if extra_label:
self.out_file = '{}/{}_{}_{}_{}.txt'.format(
self.outdir, self.label, dets, type(self).__name__,
extra_label)
else:
self.out_file = '{}/{}_{}_{}.txt'.format(
self.outdir, self.label, dets, type(self).__name__)
class SliceGridSearch(GridSearch):
""" Slice gridded search using ComputeFstat """
@helper_functions.initializer
def __init__(self, label, outdir, sftfilepattern, F0s, F1s, F2s, Alphas,
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, Lambda0=None):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
F0s, F1s, F2s, delta_F0s, delta_F1s, tglitchs, Alphas, Deltas: tuple
Length 3 tuple describing the grid for each parameter, e.g
[F0min, F0max, dF0], for a fixed value simply give [F0]. Unless
input_arrays == True, then these are the values to search at.
tref, minStartTime, maxStartTime: int
GPS seconds of the reference time, start time and end time
input_arrays: bool
if true, use the F0s, F1s, etc as is
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
self.set_out_file()
self.keys = ['_', '_', 'F0', 'F1', 'F2', 'Alpha', 'Delta']
self.ndim = 0
self.thetas = [F0s, F1s, Alphas, Deltas]
self.ndim = 4
self.search_keys = ['F0', 'F1', 'Alpha', 'Delta']
self.Lambda0 = np.array(Lambda0)
if len(self.Lambda0) != len(self.search_keys):
raise ValueError(
'Lambda0 must be of length {}'.format(len(self.search_keys)))
def run(self, factor=2, max_n_ticks=4, whspace=0.07, save=True,
**kwargs):
lbdim = 0.5 * factor # size of left/bottom margin
trdim = 0.4 * factor # size of top/right margin
plotdim = factor * self.ndim + factor * (self.ndim - 1.) * whspace
dim = lbdim + plotdim + trdim
fig, axes = plt.subplots(self.ndim, self.ndim, figsize=(dim, dim))
# Format the figure.
lb = lbdim / dim
tr = (lbdim + plotdim) / dim
fig.subplots_adjust(left=lb, bottom=lb, right=tr, top=tr,
wspace=whspace, hspace=whspace)
search = GridSearch(
self.label, self.outdir, self.sftfilepattern,
F0s=self.Lambda0[0], F1s=self.Lambda0[1], F2s=self.F2s[0],
Alphas=self.Lambda0[2], Deltas=self.Lambda0[3], tref=self.tref,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime)
for i, ikey in enumerate(self.search_keys):
setattr(search, ikey+'s', self.thetas[i])
search.label = '{}_{}'.format(self.label, ikey)
search.set_out_file()
search.run()
axes[i, i] = search.plot_1D(ikey, ax=axes[i, i], savefig=False,
x0=self.Lambda0[i]
)
setattr(search, ikey+'s', [self.Lambda0[i]])
axes[i, i].yaxis.tick_right()
axes[i, i].yaxis.set_label_position("right")
axes[i, i].set_xlabel('')
for j, jkey in enumerate(self.search_keys):
ax = axes[i, j]
if j > i:
ax.set_frame_on(False)
ax.set_xticks([])
ax.set_yticks([])
continue
ax.get_shared_x_axes().join(axes[self.ndim-1, j], ax)
if i < self.ndim - 1:
ax.set_xticklabels([])
if j < i:
ax.get_shared_y_axes().join(axes[i, i-1], ax)
if j > 0:
ax.set_yticklabels([])
if j == i:
continue
ax.xaxis.set_major_locator(
matplotlib.ticker.MaxNLocator(max_n_ticks, prune="upper"))
ax.yaxis.set_major_locator(
matplotlib.ticker.MaxNLocator(max_n_ticks, prune="upper"))
setattr(search, ikey+'s', self.thetas[i])
setattr(search, jkey+'s', self.thetas[j])
search.label = '{}_{}'.format(self.label, ikey+jkey)
search.set_out_file()
search.run()
ax = search.plot_2D(jkey, ikey, ax=ax, save=False,
y0=self.Lambda0[i], x0=self.Lambda0[j],
**kwargs)
setattr(search, ikey+'s', [self.Lambda0[i]])
setattr(search, jkey+'s', [self.Lambda0[j]])
ax.grid(lw=0.2, ls='--', zorder=10)
ax.set_xlabel('')
ax.set_ylabel('')
for i, ikey in enumerate(self.search_keys):
axes[-1, i].set_xlabel(
self.tex_labels[ikey]+self.tex_labels0[ikey])
if i > 0:
axes[i, 0].set_ylabel(
self.tex_labels[ikey]+self.tex_labels0[ikey])
axes[i, i].set_ylabel("$2\mathcal{F}$")
if save:
fig.savefig(
'{}/{}_slice_projection.png'.format(self.outdir, self.label))
else:
return fig, axes
class GridUniformPriorSearch():
@helper_functions.initializer
def __init__(self, theta_prior, NF0, NF1, label, outdir, sftfilepattern,
tref, minStartTime, maxStartTime, minCoverFreq=None,
maxCoverFreq=None, BSGL=False, detectors=None, nsegs=1,
SSBprec=None, injectSources=None):
dF0 = (theta_prior['F0']['upper'] - theta_prior['F0']['lower'])/NF0
dF1 = (theta_prior['F1']['upper'] - theta_prior['F1']['lower'])/NF1
F0s = [theta_prior['F0']['lower'], theta_prior['F0']['upper'], dF0]
F1s = [theta_prior['F1']['lower'], theta_prior['F1']['upper'], dF1]
self.search = GridSearch(
label, outdir, sftfilepattern, F0s=F0s, F1s=F1s, tref=tref,
Alphas=[theta_prior['Alpha']], Deltas=[theta_prior['Delta']],
minStartTime=minStartTime, maxStartTime=maxStartTime, BSGL=BSGL,
detectors=detectors, minCoverFreq=minCoverFreq,
injectSources=injectSources, maxCoverFreq=maxCoverFreq,
nsegs=nsegs, SSBprec=SSBprec)
def run(self):
self.search.run()
def get_2D_plot(self, **kwargs):
return self.search.plot_2D('F0', 'F1', **kwargs)
class GridGlitchSearch(GridSearch):
""" Grid search using the SemiCoherentGlitchSearch """
@helper_functions.initializer
def __init__(self, label, outdir, sftfilepattern=None, F0s=[0],
F1s=[0], F2s=[0], delta_F0s=[0], delta_F1s=[0], tglitchs=None,
Alphas=[0], Deltas=[0], tref=None, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
write_after=1000):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
F0s, F1s, F2s, delta_F0s, delta_F1s, tglitchs, Alphas, Deltas: tuple
Length 3 tuple describing the grid for each parameter, e.g
[F0min, F0max, dF0], for a fixed value simply give [F0].
tref, minStartTime, maxStartTime: int
GPS seconds of the reference time, start time and end time
For all other parameters, see pyfstat.ComputeFStat.
"""
if tglitchs is None:
self.tglitchs = [self.maxStartTime]
self.search = SemiCoherentGlitchSearch(
label=label, outdir=outdir, sftfilepattern=self.sftfilepattern,
tref=tref, minStartTime=minStartTime, maxStartTime=maxStartTime,
minCoverFreq=minCoverFreq, maxCoverFreq=maxCoverFreq,
BSGL=self.BSGL)
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
self.set_out_file()
self.keys = ['F0', 'F1', 'F2', 'Alpha', 'Delta', 'delta_F0',
'delta_F1', 'tglitch']
def get_input_data_array(self):
arrays = []
for tup in (self.F0s, self.F1s, self.F2s, self.Alphas, self.Deltas,
self.delta_F0s, self.delta_F1s, self.tglitchs):
arrays.append(self.get_array_from_tuple(tup))
input_data = []
for vals in itertools.product(*arrays):
input_data.append(vals)
self.arrays = arrays
self.input_data = np.array(input_data)
class FrequencySlidingWindow(GridSearch):
""" A sliding-window search over the Frequency """
@helper_functions.initializer
def __init__(self, label, outdir, sftfilepattern, F0s, F1, F2,
Alpha, Delta, tref, minStartTime=None,
maxStartTime=None, window_size=10*86400, window_delta=86400,
BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detectors=None, SSBprec=None, injectSources=None):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
F0s: array
Frequency range
F1, F2, Alpha, Delta: float
Fixed values to compute twoF(F) over
tref, minStartTime, maxStartTime: int
GPS seconds of the reference time, start time and end time
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
self.set_out_file()
self.nsegs = 1
self.F1s = [F1]
self.F2s = [F2]
self.Alphas = [Alpha]
self.Deltas = [Delta]
self.input_arrays = False
def inititate_search_object(self):
logging.info('Setting up search object')
self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=True,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources)
self.search.get_det_stat = (
self.search.get_fullycoherent_twoF)
def get_input_data_array(self):
arrays = []
tstarts = [self.minStartTime]
while tstarts[-1] + self.window_size < self.maxStartTime:
tstarts.append(tstarts[-1]+self.window_delta)
arrays = [tstarts]
for tup in (self.F0s, self.F1s, self.F2s,
self.Alphas, self.Deltas):
arrays.append(self.get_array_from_tuple(tup))
input_data = []
for vals in itertools.product(*arrays):
input_data.append(vals)
input_data = np.array(input_data)
input_data = np.insert(
input_data, 1, input_data[:, 0] + self.window_size, axis=1)
self.arrays = arrays
self.input_data = np.array(input_data)
def plot_sliding_window(self, F0=None, ax=None, savefig=True,
colorbar=True, timestamps=False):
data = self.data
if ax is None:
ax = plt.subplot()
tstarts = np.unique(data[:, 0])
tends = np.unique(data[:, 1])
frequencies = np.unique(data[:, 2])
twoF = data[:, -1]
tmids = (tstarts + tends) / 2.0
dts = (tmids - self.minStartTime) / 86400.
if F0:
frequencies = frequencies - F0
ax.set_ylabel('Frequency - $f_0$ [Hz] \n $f_0={:0.2f}$'.format(F0))
else:
ax.set_ylabel('Frequency [Hz]')
twoF = twoF.reshape((len(tmids), len(frequencies)))
Y, X = np.meshgrid(frequencies, dts)
pax = ax.pcolormesh(X, Y, twoF)
if colorbar:
cb = plt.colorbar(pax, ax=ax)
cb.set_label('$2\mathcal{F}$')
ax.set_xlabel(
r'Mid-point (days after $t_\mathrm{{start}}$={})'.format(
self.minStartTime))
ax.set_title(
'Sliding window length = {} days in increments of {} days'
.format(self.window_size/86400, self.window_delta/86400),
)
if timestamps:
axT = ax.twiny()
axT.set_xlim(tmids[0]*1e-9, tmids[-1]*1e-9)
axT.set_xlabel('Mid-point timestamp [GPS $10^{9}$ s]')
ax.set_title(ax.get_title(), y=1.18)
if savefig:
plt.tight_layout()
plt.savefig(
'{}/{}_sliding_window.png'.format(self.outdir, self.label))
else:
return ax
class EarthTest(GridSearch):
""" """
tex_labels = {'deltaRadius': '$\Delta R$ [m]',
'phaseOffset': 'phase-offset [rad]',
'deltaPspin': '$\Delta P_\mathrm{spin}$ [s]'}
@helper_functions.initializer
def __init__(self, label, outdir, sftfilepattern, deltaRadius,
phaseOffset, deltaPspin, F0, F1, F2, Alpha,
Delta, tref=None, minStartTime=None, maxStartTime=None,
BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detectors=None, injectSources=None,
assumeSqrtSX=None):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
F0, F1, F2, Alpha, Delta: float
tref, minStartTime, maxStartTime: int
GPS seconds of the reference time, start time and end time
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
self.nsegs = 1
self.F0s = [F0]
self.F1s = [F1]
self.F2s = [F2]
self.Alphas = [Alpha]
self.Deltas = [Delta]
self.duration = maxStartTime - minStartTime
self.deltaRadius = np.atleast_1d(deltaRadius)
self.phaseOffset = np.atleast_1d(phaseOffset)
self.phaseOffset = self.phaseOffset + 1e-12 # Hack to stop cached data being used
self.deltaPspin = np.atleast_1d(deltaPspin)
self.set_out_file()
self.SSBprec = lalpulsar.SSBPREC_RELATIVISTIC
self.keys = ['deltaRadius', 'phaseOffset', 'deltaPspin']
self.prior_widths = [
np.max(self.deltaRadius)-np.min(self.deltaRadius),
np.max(self.phaseOffset)-np.min(self.phaseOffset),
np.max(self.deltaPspin)-np.min(self.deltaPspin)]
if hasattr(self, 'search') is False:
self.inititate_search_object()
def get_input_data_array(self):
logging.info("Generating input data array")
coord_arrays = [self.deltaRadius, self.phaseOffset, self.deltaPspin]
input_data = []
for vals in itertools.product(*coord_arrays):
input_data.append(vals)
self.input_data = np.array(input_data)
self.coord_arrays = coord_arrays
def run_special(self):
vals = [self.minStartTime, self.maxStartTime, self.F0, self.F1,
self.F2, self.Alpha, self.Delta]
self.special_data = {'zero': [0, 0, 0]}
for key, (dR, dphi, dP) in self.special_data.iteritems():
rescaleRadius = (1 + dR / lal.REARTH_SI)
rescalePeriod = (1 + dP / lal.DAYSID_SI)
lalpulsar.BarycenterModifyEarthRotation(
rescaleRadius, dphi, rescalePeriod, self.tref)
FS = self.search.get_det_stat(*vals)
self.special_data[key] = list([dR, dphi, dP]) + [FS]
def run(self):
self.run_special()
self.get_input_data_array()
old_data = self.check_old_data_is_okay_to_use()
if old_data is not False:
self.data = old_data
return
data = []
vals = [self.minStartTime, self.maxStartTime, self.F0, self.F1,
self.F2, self.Alpha, self.Delta]
for (dR, dphi, dP) in tqdm(self.input_data):
rescaleRadius = (1 + dR / lal.REARTH_SI)
rescalePeriod = (1 + dP / lal.DAYSID_SI)
lalpulsar.BarycenterModifyEarthRotation(
rescaleRadius, dphi, rescalePeriod, self.tref)
FS = self.search.get_det_stat(*vals)
data.append(list([dR, dphi, dP]) + [FS])
data = np.array(data, dtype=np.float)
logging.info('Saving data to {}'.format(self.out_file))
np.savetxt(self.out_file, data, delimiter=' ')
self.data = data
def marginalised_bayes_factor(self, prior_widths=None):
if prior_widths is None:
prior_widths = self.prior_widths
ndims = self.data.shape[1] - 1
params = np.array([np.unique(self.data[:, j]) for j in range(ndims)])
twoF = self.data[:, -1].reshape(tuple([len(p) for p in params]))
F = twoF / 2.0
for i, x in enumerate(params[::-1]):
if len(x) > 1:
dx = x[1] - x[0]
F = logsumexp(F, axis=-1)+np.log(dx)-np.log(prior_widths[-1-i])
else:
F = np.squeeze(F, axis=-1)
marginalised_F = np.atleast_1d(F)[0]
F_at_zero = self.special_data['zero'][-1]/2.0
max_idx = np.argmax(self.data[:, -1])
max_F = self.data[max_idx, -1]/2.0
max_F_params = self.data[max_idx, :-1]
logging.info('F at zero = {:.1f}, marginalised_F = {:.1f},'
' max_F = {:.1f} ({})'.format(
F_at_zero, marginalised_F, max_F, max_F_params))
return F_at_zero - marginalised_F, (F_at_zero - max_F) / F_at_zero
def plot_corner(self, prior_widths=None, fig=None, axes=None,
projection='log_mean'):
Bsa, FmaxMismatch = self.marginalised_bayes_factor(prior_widths)
data = self.data[:, -1].reshape(
(len(self.deltaRadius), len(self.phaseOffset),
len(self.deltaPspin)))
xyz = [self.deltaRadius/lal.REARTH_SI, self.phaseOffset/(np.pi),
self.deltaPspin/60.]
labels = [r'$\frac{\Delta R}{R_\mathrm{Earth}}$',
r'$\frac{\Delta \phi}{\pi}$',
r'$\Delta P_\mathrm{spin}$ [min]',
r'$2\mathcal{F}$']
from projection_matrix import projection_matrix
fig, axes = projection_matrix(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},'
r' $\frac{{\mathcal{{F}}_0-\mathcal{{F}}_\mathrm{{max}}}}'
r'{{\mathcal{{F}}_0}}={:.1e}$'
.format(self.duration/86400, self.F0, Bsa, FmaxMismatch), y=0.99,
size=14)
fig.savefig('{}/{}_projection_matrix.png'.format(
self.outdir, self.label))
def plot(self, key, prior_widths=None):
Bsa, FmaxMismatch = self.marginalised_bayes_factor(prior_widths)
rescales_defaults = {'deltaRadius': 1/lal.REARTH_SI,
'phaseOffset': 1/np.pi,
'deltaPspin': 1}
labels = {'deltaRadius': r'$\frac{\Delta R}{R_\mathrm{Earth}}$',
'phaseOffset': r'$\frac{\Delta \phi}{\pi}$',
'deltaPspin': r'$\Delta P_\mathrm{spin}$ [s]'
}
fig, ax = self.plot_1D(key, xrescale=rescales_defaults[key],
xlabel=labels[key], savefig=False)
ax.set_title(
'T={} days, $f$={} Hz, $\log\mathcal{{B}}_{{S/A}}$={:.1f}'
.format(self.duration/86400, self.F0, Bsa))
fig.tight_layout()
fig.savefig('{}/{}_1D.png'.format(self.outdir, self.label))
class DMoff_NO_SPIN(GridSearch):
""" DMoff test using SSBPREC_NO_SPIN """
@helper_functions.initializer
def __init__(self, par, label, outdir, sftfilepattern, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
detectors=None, injectSources=None, assumeSqrtSX=None):
"""
Parameters
----------
par: dict, str
Either a par dictionary (containing 'F0', 'F1', 'Alpha', 'Delta'
and 'tref') or a path to a .par file to read in the F0, F1 etc
label, outdir: str
A label and directory to read/write data from/to
sftfilepattern: str
Pattern to match SFTs using wildcards (*?) and ranges [0-9];
mutiple patterns can be given separated by colons.
minStartTime, maxStartTime: int
GPS seconds of the start time and end time
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
if os.path.isdir(outdir) is False:
os.mkdir(outdir)
if type(par) == dict:
self.par = par
elif type(par) == str and os.path.isfile(par):
self.par = read_par(filename=par)
else:
raise ValueError('The .par file does not exist')
self.nsegs = 1
self.BSGL = False
self.tref = self.par['tref']
self.F1s = [self.par.get('F1', 0)]
self.F2s = [self.par.get('F2', 0)]
self.Alphas = [self.par['Alpha']]
self.Deltas = [self.par['Delta']]
self.Re = 6.371e6
self.c = 2.998e8
a0 = self.Re/self.c # *np.cos(self.par['Delta'])
self.m0 = np.max([4, int(np.ceil(2*np.pi*self.par['F0']*a0))])
logging.info(
'Setting up DMoff_NO_SPIN search with m0 = {}'.format(self.m0))
def get_results(self):
""" Compute the three summed detection statistics
Returns
-------
m0, twoF_SUM, twoFstar_SUM_SIDEREAL, twoFstar_SUM_TERRESTRIAL
"""
self.SSBprec = lalpulsar.SSBPREC_RELATIVISTIC
self.set_out_file('SSBPREC_RELATIVISTIC')
self.F0s = [self.par['F0']+j/lal.DAYSID_SI for j in range(-4, 5)]
self.run()
twoF_SUM = np.sum(self.data[:, -1])
self.SSBprec = lalpulsar.SSBPREC_NO_SPIN
self.set_out_file('SSBPREC_NO_SPIN')
self.F0s = [self.par['F0']+j/lal.DAYSID_SI
for j in range(-self.m0, self.m0+1)]
self.run()
twoFstar_SUM = np.sum(self.data[:, -1])
self.set_out_file('SSBPREC_NO_SPIN_TERRESTRIAL')
self.F0s = [self.par['F0']+j/lal.DAYJUL_SI
for j in range(-self.m0, self.m0+1)]
self.run()
twoFstar_SUM_terrestrial = np.sum(self.data[:, -1])
return self.m0, twoF_SUM, twoFstar_SUM, twoFstar_SUM_terrestrial