From 7cb7b01e26f46441bf6bdf99b68dc53ec7640331 Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Wed, 21 Sep 2016 11:23:59 +0200
Subject: [PATCH] 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.
---
pyfstat.py | 320 -------------------------------------------------
tests/tests.py | 87 --------------
2 files changed, 407 deletions(-)
diff --git a/pyfstat.py b/pyfstat.py
index 048ef55..130846b 100755
--- a/pyfstat.py
+++ b/pyfstat.py
@@ -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 """
diff --git a/tests/tests.py b/tests/tests.py
index 7f90a7d..852eb92 100644
--- a/tests/tests.py
+++ b/tests/tests.py
@@ -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
--
GitLab