Commit 7cb7b01e authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Remove old methods

Removes the fully-coherent search which was just a wrapper to CFS_v2 and
all of the `slow' methods which used command line arguments rather than
the python interface.
parent 53a140d5
......@@ -141,272 +141,6 @@ class BaseSearchClass(object):
return thetas
class FullyCoherentNarrowBandSearch(BaseSearchClass):
""" Search over a narrow band of F0, F1, and F2 """
@initializer
def __init__(self, label, outdir, sftlabel=None, sftdir=None,
tglitch=None, tref=None, tstart=None, Alpha=None, Delta=None,
duration=None, Writer=None):
"""
Parameters
----------
label, outdir: str
A label and directory to read/write data from/to
sftlabel, sftdir: str
A label and directory in which to find the relevant sft file. If
None use label and outdir
"""
if self.sftlabel is None:
self.sftlabel = self.label
if self.sftdir is None:
self.sftdir = self.outdir
self.tend = self.tstart + self.duration
self.calculate_best_fit_F0_and_F1(Writer)
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
def calculate_best_fit_F0_and_F1(self, writer):
R = (writer.tglitch - writer.tstart) / float(writer.duration)
self.F0_min = (writer.F0 + writer.delta_F0*(1-R)**2*(2*R+1) +
3*writer.delta_phi*R*(1-R)/np.pi/writer.duration +
# .5*writer.duration*writer.F1*(1-R) +
.5*writer.duration*writer.delta_F1*(1-R)**3*(1+R))
self.F1_min = (writer.F1 + writer.delta_F1*(1-R)**3*(6*R**2+3*R+1) +
30*writer.delta_F0*R**2*(1-R)**2/writer.duration +
30*writer.delta_phi*R*(1-R)*(2*R-1)/np.pi/writer.duration**2
)
def get_grid_setup(self, m, n, searchF2=False):
""" Calc. the grid parameters of bands given the metric-mismatch
Parameters
----------
m: float in [0, 1]
The mismatch spacing between adjacent grid points
n: int
Number of grid points to search
"""
DeltaF0 = np.sqrt(12 * m) / (np.pi * self.duration)
DeltaF1 = np.sqrt(720 * m) / (np.pi * self.duration**2.0)
DeltaF2 = np.sqrt(100800 * m) / (np.pi * self.duration**3.0)
# Calculate the width of bands given n
F0Band = n * DeltaF0
F1Band = n * DeltaF1
F2Band = n * DeltaF2
# Search takes the lowest frequency in the band
F0_bottom = self.F0_min - .5 * F0Band
F1_bottom = self.F1_min - .5 * F1Band
if searchF2:
F2_bottom = self.F2_min-.5*self.F2Band # Not yet implemented
else:
F2_bottom = 0 # Broken functionality
F2Band = 0
Messg = ["Automated search for {}:".format(self.label),
"Grid parameters : m={}, n={}".format(m, n),
"Reference time: {}".format(self.tref),
"Analytic best-fit values : {}, {}".format(
self.F0_min, self.F1_min),
"F0Band : {} -- {}".format(F0_bottom,
F0_bottom + F0Band),
"F1Band : {} -- {}".format(F1_bottom,
F1_bottom + F1Band),
"F2Band : {} -- {}".format(F2_bottom,
F2_bottom + F2Band),
]
logging.info("\n ".join(Messg))
return (F0_bottom, DeltaF0, F0Band,
F1_bottom, DeltaF1, F1Band,
F2_bottom, DeltaF2, F2Band)
def run_computefstatistic_slow(self, m, n, search_F2=False):
""" Compute the f statistic fully-coherently over a grid """
(F0_bottom, DeltaF0, F0Band,
F1_bottom, DeltaF1, F1Band,
F2_bottom, DeltaF2, F2Band) = self.get_grid_setup(m, n)
c_l = []
c_l.append("lalapps_ComputeFstatistic_v2")
c_l.append("--Freq={}".format(F0_bottom))
c_l.append("--dFreq={}".format(DeltaF0))
c_l.append("--FreqBand={}".format(F0Band))
c_l.append("--f1dot={}".format(F1_bottom))
c_l.append("--df1dot={}".format(DeltaF1))
c_l.append("--f1dotBand={}".format(F1Band))
if search_F2:
c_l.append("--f2dot={}".format(F2_bottom))
c_l.append("--df2dot={}".format(DeltaF2))
c_l.append("--f2dotBand={}".format(F2Band))
else:
c_l.append("--f2dot={}".format(F2_bottom))
c_l.append("--DataFiles='{}'".format(
self.outdir+"/*SFT_"+self.label+"*sft"))
c_l.append("--refTime={:10.6f}".format(self.tref))
c_l.append("--outputFstat='{}'".format(self.fs_file_name))
c_l.append("--Alpha={}".format(self.Alpha))
c_l.append("--Delta={}".format(self.Delta))
c_l.append("--minStartTime={}".format(self.tstart))
c_l.append("--maxStartTime={}".format(self.tend))
logging.info("Executing: " + " ".join(c_l) + "\n")
os.system(" ".join(c_l))
self.read_in_fstat()
def run_computefstatistic(self, dFreq=0, numFreqBins=1):
""" Compute the f statistic fully-coherently over a grid """
constraints = lalpulsar.SFTConstraints()
FstatOptionalArgs = lalpulsar.FstatOptionalArgsDefaults
minCoverFreq = 29
maxCoverFreq = 31
SFTCatalog = lalpulsar.SFTdataFind(
self.sftdir+'/*_'+self.sftlabel+"*sft", constraints)
ephemerides = lalpulsar.InitBarycenter(
'~/lalsuite-install/share/lalpulsar/earth00-19-DE421.dat.gz',
'~/lalsuite-install/share/lalpulsar/sun00-19-DE421.dat.gz')
whatToCompute = lalpulsar.FSTATQ_2F
FstatInput = lalpulsar.CreateFstatInput(SFTCatalog,
minCoverFreq,
maxCoverFreq,
dFreq,
ephemerides,
FstatOptionalArgs
)
PulsarDopplerParams = lalpulsar.PulsarDopplerParams()
PulsarDopplerParams.refTime = self.tref
PulsarDopplerParams.Alpha = self.Alpha
PulsarDopplerParams.Delta = self.Delta
PulsarDopplerParams.fkdot = np.array([self.F0_min-dFreq*numFreqBins/2.,
self.F1_min, 0, 0, 0, 0, 0])
FstatResults = lalpulsar.FstatResults()
lalpulsar.ComputeFstat(FstatResults,
FstatInput,
PulsarDopplerParams,
numFreqBins,
whatToCompute
)
self.search_F0 = (np.linspace(0, dFreq * numFreqBins, numFreqBins)
+ FstatResults.doppler.fkdot[0])
self.search_F1 = np.zeros(numFreqBins) + self.F1_min
self.search_F2 = np.zeros(numFreqBins) + 0
self.search_FS = FstatResults.twoF
def read_in_fstat(self):
"""
Read in data from *_FS.dat file as produced by ComputeFStatistic_v2
"""
data = np.genfromtxt(self.fs_file_name, comments="%")
# If none of the components are varying:
if data.ndim == 1:
self.search_F0 = data[0]
self.search_F1 = data[3]
self.search_F2 = data[4]
self.search_FS = data[6]
return
search_F0 = data[:, 0]
search_F1 = data[:, 3]
search_F2 = data[:, 4]
search_FS = data[:, 6]
NF0 = len(np.unique(search_F0))
NF1 = len(np.unique(search_F1))
NF2 = len(np.unique(search_F2))
shape = (NF2, NF1, NF0)
self.data_shape = shape
self.search_F0 = np.squeeze(np.reshape(search_F0,
newshape=shape).transpose())
self.search_F1 = np.squeeze(np.reshape(search_F1,
newshape=shape).transpose())
self.search_F2 = np.squeeze(np.reshape(search_F2,
newshape=shape).transpose())
self.search_FS = np.squeeze(np.reshape(search_FS,
newshape=shape).transpose())
def get_FS_max(self):
""" Returns the maximum FS and the corresponding F0, F1, and F2 """
if np.shape(self.search_FS) == ():
return self.search_F0, self.search_F1, self.search_F2, self.search_FS
else:
max_idx = np.unravel_index(self.search_FS.argmax(), self.search_FS.shape)
return (self.search_F0[max_idx], self.search_F1[max_idx],
self.search_F2[max_idx], self.search_FS[max_idx])
def plot_output(self, output_type='FS', perfectlymatched_FS=None,
fig=None, ax=None, savefig=False, title=None):
"""
Plot the output of the *_FS.dat file as a contour plot
Parameters
----------
output_type: str
one of 'FS', 'rho', or 'mismatch'
perfectlymatched_FS: float
the 2F of a perfectly matched signal against which to
compute the mismatch
"""
resF0 = self.search_F0 - self.F0_min
resF1 = self.search_F1 - self.F1_min
if output_type == 'FS':
Z = self.search_FS
zlabel = r'$2\bar{\mathcal{F}}$'
elif output_type == 'rho':
Z = self.search_FS - 4
zlabel = r'\rho^{2}'
elif output_type == 'mismatch':
rho2 = self.search_FS - 4
perfectlymatched_rho2 = perfectlymatched_FS - 4
if perfectlymatched_FS:
Z = 1 - (rho2) / (perfectlymatched_rho2)
else:
raise ValueError('Plotting the mismatch requires a value for'
' the parameter perfectlymatched_rho2')
zlabel = 'mismatch'
if ax is None:
fig, ax = plt.subplots()
plt.rcParams['font.size'] = 12
pax = ax.pcolormesh(resF0, resF1, Z, cmap=plt.cm.viridis,
vmin=0, vmax=1)
fig.colorbar(pax, label=zlabel, ax=ax)
ax.set_xlim(np.min(resF0), np.max(resF0))
ax.set_ylim(np.min(resF1), np.max(resF1))
ax.set_xlabel(r'$f_0 - f_\textrm{min}$')
ax.set_ylabel(r'$\dot{f}_0 - \dot{f}_\textrm{min}$')
ax.set_title(self.label)
plt.tight_layout()
if savefig:
fig.savefig('output_{}.png'.format(self.label))
class SemiCoherentGlitchSearch(BaseSearchClass):
""" A semi-coherent glitch search
......@@ -415,8 +149,6 @@ class SemiCoherentGlitchSearch(BaseSearchClass):
fully-coherent F-stat in each segment is averaged to give the semi-coherent
F-stat
"""
# Copy methods
read_in_fstat = FullyCoherentNarrowBandSearch.__dict__['read_in_fstat']
@initializer
def __init__(self, label, outdir, tref, tstart, tend, sftlabel=None,
......@@ -591,58 +323,6 @@ class SemiCoherentGlitchSearch(BaseSearchClass):
useFReg)
return 2*FS.F_mn.data[0][0]
def compute_glitch_fstat_slow(self, F0, F1, F2, Alpha, Delta, delta_F0,
delta_F1, tglitch):
""" Compute the semi-coherent F-stat """
theta = [F0, F1, F2]
delta_theta = [delta_F0, delta_F1, 0]
tref = self.tref
theta_at_glitch = self.shift_coefficients(theta, tglitch - tref)
theta_post_glitch_at_glitch = theta_at_glitch + delta_theta
theta_post_glitch = self.shift_coefficients(
theta_post_glitch_at_glitch, tref - tglitch)
FsegA = self.run_computefstatistic_single_point_slow(
tref, self.tstart, tglitch, theta[0], theta[1], theta[2], Alpha,
Delta)
FsegB = self.run_computefstatistic_single_point_slow(
tref, tglitch, self.tend, theta_post_glitch[0],
theta_post_glitch[1], theta_post_glitch[2], Alpha,
Delta)
return (FsegA + FsegB) / 2.
def run_computefstatistic_single_point_slow(self, tref, tstart, tend, F0,
F1, F2, Alpha, Delta):
""" Compute the f statistic fully-coherently at a single point """
c_l = []
c_l.append("lalapps_ComputeFstatistic_v2")
c_l.append("--Freq={}".format(F0))
c_l.append("--f1dot={}".format(F1))
c_l.append("--f2dot={}".format(F2))
c_l.append("--DataFiles='{}'".format(
self.outdir+"/*SFT_"+self.label+"*sft"))
c_l.append("--refTime={:10.6f}".format(tref))
c_l.append("--outputFstat='{}'".format(self.fs_file_name))
c_l.append("--Alpha={}".format(Alpha))
c_l.append("--Delta={}".format(Delta))
c_l.append("--minStartTime={}".format(tstart))
c_l.append("--maxStartTime={}".format(tend))
logging.info("Executing: " + " ".join(c_l) + "\n")
os.system(" ".join(c_l))
self.read_in_fstat()
return self.search_FS
class MCMCGlitchSearch(BaseSearchClass):
""" MCMC search using the SemiCoherentGlitchSearch """
......
......@@ -75,41 +75,6 @@ class TestBaseSearchClass(unittest.TestCase):
rtol=1e-9, atol=1e-9))
class TestFullyCoherentNarrowBandSearch(unittest.TestCase):
label = "Test"
outdir = 'TestData'
def test_compute_fstat(self):
Writer = glitch_tools.Writer(self.label, outdir=self.outdir)
Writer.make_data()
search = glitch_searches.FullyCoherentNarrowBandSearch(
self.label, self.outdir, tref=Writer.tref, Alpha=Writer.Alpha,
Delta=Writer.Delta, duration=Writer.duration, tstart=Writer.tstart,
Writer=Writer)
search.run_computefstatistic_slow(m=1e-3, n=0)
_, _, _, FS_max_slow = search.get_FS_max()
search.run_computefstatistic(dFreq=0, numFreqBins=1)
_, _, _, FS_max = search.get_FS_max()
self.assertTrue(
np.abs(FS_max-FS_max_slow)/FS_max_slow < 0.1)
def test_compute_fstat_against_predict_fstat(self):
Writer = glitch_tools.Writer(self.label, outdir=self.outdir)
Writer.make_data()
Writer.run_makefakedata()
predicted_FS = Writer.predict_fstat()
search = glitch_searches.FullyCoherentNarrowBandSearch(
self.label, self.outdir, tref=Writer.tref, Alpha=Writer.Alpha,
Delta=Writer.Delta, duration=Writer.duration, tstart=Writer.tstart,
Writer=Writer)
search.run_computefstatistic(dFreq=0, numFreqBins=1)
_, _, _, FS_max = search.get_FS_max()
self.assertTrue(np.abs(predicted_FS-FS_max)/FS_max < 0.5)
class TestSemiCoherentGlitchSearch(unittest.TestCase):
label = "Test"
outdir = 'TestData'
......@@ -133,58 +98,6 @@ class TestSemiCoherentGlitchSearch(unittest.TestCase):
print predicted_FS, FS
self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.1)
def test_run_computefstatistic_single_point_slow(self):
Writer = pyfstat.Writer(self.label, outdir=self.outdir)
Writer.make_data()
predicted_FS = Writer.predict_fstat()
search = pyfstat.SemiCoherentGlitchSearch(
label=Writer.label, outdir=Writer.outdir, tref=Writer.tref,
tstart=Writer.tstart, tend=Writer.tend)
FS = search.run_computefstatistic_single_point_slow(search.tref,
search.tstart,
search.tend,
Writer.F0,
Writer.F1,
Writer.F2,
Writer.Alpha,
Writer.Delta)
self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.1)
def test_compute_glitch_fstat_slow(self):
duration = 100*86400
dtglitch = 100*43200
delta_F0 = 0
Writer = pyfstat.Writer(self.label, outdir=self.outdir,
duration=duration, dtglitch=dtglitch,
delta_F0=delta_F0)
Writer.make_data()
search = pyfstat.SemiCoherentGlitchSearch(
label=Writer.label, outdir=Writer.outdir, tref=Writer.tref,
tstart=Writer.tstart, tend=Writer.tend)
FS = search.compute_glitch_fstat_slow(Writer.F0, Writer.F1, Writer.F2,
Writer.Alpha, Writer.Delta,
Writer.delta_F0, Writer.delta_F1,
Writer.tglitch)
# Compute the predicted semi-coherent glitch Fstat
tstart = Writer.tstart
tend = Writer.tend
Writer.tend = tstart + dtglitch
FSA = Writer.predict_fstat()
Writer.tstart = tstart + dtglitch
Writer.tend = tend
FSB = Writer.predict_fstat()
predicted_FS = .5*(FSA + FSB)
print(predicted_FS, FS)
self.assertTrue(np.abs((FS - predicted_FS))/predicted_FS < 0.1)
def test_compute_nglitch_fstat(self):
duration = 100*86400
dtglitch = 100*43200
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment