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

Renaming functions for clarity of twoF and fix issue of F -> likelihood

Fixes #4 .
parent 6e0b510a
...@@ -576,19 +576,9 @@ class ComputeFstat(BaseSearchClass): ...@@ -576,19 +576,9 @@ class ComputeFstat(BaseSearchClass):
self.windowRange.tauBand = 0 self.windowRange.tauBand = 0
self.windowRange.dtau = 1 self.windowRange.dtau = 1
def compute_fullycoherent_det_stat_single_point( def get_fullycoherent_twoF(self, tstart, tend, F0, F1, F2, Alpha, Delta,
self, F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, asini=None, period=None, ecc=None, tp=None,
tp=None, argp=None): argp=None):
""" Compute the fully-coherent det. statistic at a single point """
return self.run_computefstatistic_single_point(
self.minStartTime, self.maxStartTime, F0, F1, F2, Alpha, Delta,
asini, period, ecc, tp, argp)
def run_computefstatistic_single_point(self, tstart, tend, F0, F1,
F2, Alpha, Delta, asini=None,
period=None, ecc=None, tp=None,
argp=None):
""" Returns twoF or ln(BSGL) fully-coherently at a single point """ """ Returns twoF or ln(BSGL) fully-coherently at a single point """
self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0]) self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0])
...@@ -684,7 +674,7 @@ class ComputeFstat(BaseSearchClass): ...@@ -684,7 +674,7 @@ class ComputeFstat(BaseSearchClass):
self.transient = True self.transient = True
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
for tau in taus: for tau in taus:
twoFs.append(self.run_computefstatistic_single_point( twoFs.append(self.get_fullycoherent_twoF(
tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2, tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2,
Alpha=Alpha, Delta=Delta, asini=asini, period=period, ecc=ecc, Alpha=Alpha, Delta=Delta, asini=asini, period=period, ecc=ecc,
tp=tp, argp=argp)) tp=tp, argp=argp))
...@@ -875,7 +865,7 @@ class SemiCoherentSearch(ComputeFstat): ...@@ -875,7 +865,7 @@ class SemiCoherentSearch(ComputeFstat):
'Semi-coherent end time {} after last SFT timestamp {}' 'Semi-coherent end time {} after last SFT timestamp {}'
.format(self.tboundaries[-1], self.SFT_timestamps[-1])) .format(self.tboundaries[-1], self.SFT_timestamps[-1]))
def run_semi_coherent_computefstatistic_single_point( def get_semicoherent_twoF(
self, F0, F1, F2, Alpha, Delta, asini=None, self, F0, F1, F2, Alpha, Delta, asini=None,
period=None, ecc=None, tp=None, argp=None, period=None, ecc=None, tp=None, argp=None,
record_segments=False): record_segments=False):
...@@ -992,7 +982,7 @@ class SemiCoherentGlitchSearch(ComputeFstat): ...@@ -992,7 +982,7 @@ class SemiCoherentGlitchSearch(ComputeFstat):
self.binary = False self.binary = False
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
def compute_nglitch_fstat(self, F0, F1, F2, Alpha, Delta, *args): def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
""" Returns the semi-coherent glitch summed twoF """ """ Returns the semi-coherent glitch summed twoF """
args = list(args) args = list(args)
...@@ -1013,7 +1003,7 @@ class SemiCoherentGlitchSearch(ComputeFstat): ...@@ -1013,7 +1003,7 @@ class SemiCoherentGlitchSearch(ComputeFstat):
for i, theta_i_at_tref in enumerate(thetas): for i, theta_i_at_tref in enumerate(thetas):
ts, te = tboundaries[i], tboundaries[i+1] ts, te = tboundaries[i], tboundaries[i+1]
twoFVal = self.run_computefstatistic_single_point( twoFVal = self.get_fullycoherent_twoF(
ts, te, theta_i_at_tref[1], theta_i_at_tref[2], ts, te, theta_i_at_tref[1], theta_i_at_tref[2],
theta_i_at_tref[3], Alpha, Delta) theta_i_at_tref[3], Alpha, Delta)
twoFSum += twoFVal twoFSum += twoFVal
...@@ -1039,14 +1029,14 @@ class SemiCoherentGlitchSearch(ComputeFstat): ...@@ -1039,14 +1029,14 @@ class SemiCoherentGlitchSearch(ComputeFstat):
theta_post_glitch = self._shift_coefficients( theta_post_glitch = self._shift_coefficients(
theta_post_glitch_at_glitch, tref - tglitch) theta_post_glitch_at_glitch, tref - tglitch)
twoFsegA = self.run_computefstatistic_single_point( twoFsegA = self.get_fullycoherent_twoF(
self.minStartTime, tglitch, theta[0], theta[1], theta[2], Alpha, self.minStartTime, tglitch, theta[0], theta[1], theta[2], Alpha,
Delta) Delta)
if tglitch == self.maxStartTime: if tglitch == self.maxStartTime:
return twoFsegA return twoFsegA
twoFsegB = self.run_computefstatistic_single_point( twoFsegB = self.get_fullycoherent_twoF(
tglitch, self.maxStartTime, theta_post_glitch[0], tglitch, self.maxStartTime, theta_post_glitch[0],
theta_post_glitch[1], theta_post_glitch[2], Alpha, theta_post_glitch[1], theta_post_glitch[2], Alpha,
Delta) Delta)
......
...@@ -60,7 +60,7 @@ class GridSearch(BaseSearchClass): ...@@ -60,7 +60,7 @@ class GridSearch(BaseSearchClass):
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources, injectSources=self.injectSources,
assumeSqrtSX=self.assumeSqrtSX) assumeSqrtSX=self.assumeSqrtSX)
self.search.get_det_stat = self.search.run_computefstatistic_single_point self.search.get_det_stat = self.search.get_fullycoherent_twoF
else: else:
self.search = SemiCoherentSearch( self.search = SemiCoherentSearch(
label=self.label, outdir=self.outdir, tref=self.tref, label=self.label, outdir=self.outdir, tref=self.tref,
...@@ -71,7 +71,7 @@ class GridSearch(BaseSearchClass): ...@@ -71,7 +71,7 @@ class GridSearch(BaseSearchClass):
injectSources=self.injectSources) injectSources=self.injectSources)
def cut_out_tstart_tend(*vals): def cut_out_tstart_tend(*vals):
return self.search.run_semi_coherent_computefstatistic_single_point(*vals[2:]) return self.search.get_semicoherent_twoF(*vals[2:])
self.search.get_det_stat = cut_out_tstart_tend self.search.get_det_stat = cut_out_tstart_tend
def get_array_from_tuple(self, x): def get_array_from_tuple(self, x):
...@@ -397,7 +397,7 @@ class FrequencySlidingWindow(GridSearch): ...@@ -397,7 +397,7 @@ class FrequencySlidingWindow(GridSearch):
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources) injectSources=self.injectSources)
self.search.get_det_stat = ( self.search.get_det_stat = (
self.search.run_computefstatistic_single_point) self.search.get_fullycoherent_twoF)
def get_input_data_array(self): def get_input_data_array(self):
arrays = [] arrays = []
......
...@@ -160,9 +160,9 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -160,9 +160,9 @@ class MCMCSearch(core.BaseSearchClass):
def logl(self, theta, search): def logl(self, theta, search):
for j, theta_i in enumerate(self.theta_idxs): for j, theta_i in enumerate(self.theta_idxs):
self.fixed_theta[theta_i] = theta[j] self.fixed_theta[theta_i] = theta[j]
FS = search.compute_fullycoherent_det_stat_single_point( twoF = search.get_fullycoherent_twoF(
*self.fixed_theta) self.minStartTime, self.maxStartTime, *self.fixed_theta)
return FS + self.likelihoodcoef return twoF/2.0 + self.likelihoodcoef
def _unpack_input_theta(self): def _unpack_input_theta(self):
full_theta_keys = ['F0', 'F1', 'F2', 'Alpha', 'Delta'] full_theta_keys = ['F0', 'F1', 'F2', 'Alpha', 'Delta']
...@@ -1244,11 +1244,11 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -1244,11 +1244,11 @@ class MCMCSearch(core.BaseSearchClass):
""" """
if any(np.isposinf(self.lnlikes)): if any(np.isposinf(self.lnlikes)):
logging.info('twoF values contain positive infinite values') logging.info('lnlike values contain positive infinite values')
if any(np.isneginf(self.lnlikes)): if any(np.isneginf(self.lnlikes)):
logging.info('twoF values contain negative infinite values') logging.info('lnlike values contain negative infinite values')
if any(np.isnan(self.lnlikes)): if any(np.isnan(self.lnlikes)):
logging.info('twoF values contain nan') logging.info('lnlike values contain nan')
idxs = np.isfinite(self.lnlikes) idxs = np.isfinite(self.lnlikes)
jmax = np.nanargmax(self.lnlikes[idxs]) jmax = np.nanargmax(self.lnlikes[idxs])
maxlogl = self.lnlikes[jmax] maxlogl = self.lnlikes[jmax]
...@@ -1262,7 +1262,7 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -1262,7 +1262,7 @@ class MCMCSearch(core.BaseSearchClass):
maxtwoF = self.logl(p, self.search) maxtwoF = self.logl(p, self.search)
self.search.BSGL = self.BSGL self.search.BSGL = self.BSGL
else: else:
maxtwoF = maxlogl - self.likelihoodcoef maxtwoF = (maxlogl - self.likelihoodcoef)*2
repeats = [] repeats = []
for i, k in enumerate(self.theta_keys): for i, k in enumerate(self.theta_keys):
...@@ -1633,8 +1633,8 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1633,8 +1633,8 @@ class MCMCGlitchSearch(MCMCSearch):
for j, theta_i in enumerate(self.theta_idxs): for j, theta_i in enumerate(self.theta_idxs):
self.fixed_theta[theta_i] = theta[j] self.fixed_theta[theta_i] = theta[j]
FS = search.compute_nglitch_fstat(*self.fixed_theta) twoF = search.get_semicoherent_nglitch_twoF(*self.fixed_theta)
return FS + self.likelihoodcoef return twoF/2.0 + self.likelihoodcoef
def _unpack_input_theta(self): def _unpack_input_theta(self):
glitch_keys = ['delta_F0', 'delta_F1', 'tglitch'] glitch_keys = ['delta_F0', 'delta_F1', 'tglitch']
...@@ -1839,9 +1839,9 @@ class MCMCSemiCoherentSearch(MCMCSearch): ...@@ -1839,9 +1839,9 @@ class MCMCSemiCoherentSearch(MCMCSearch):
def logl(self, theta, search): def logl(self, theta, search):
for j, theta_i in enumerate(self.theta_idxs): for j, theta_i in enumerate(self.theta_idxs):
self.fixed_theta[theta_i] = theta[j] self.fixed_theta[theta_i] = theta[j]
FS = search.run_semi_coherent_computefstatistic_single_point( twoF = search.get_semicoherent_twoF(
*self.fixed_theta) *self.fixed_theta)
return FS + self.likelihoodcoef return twoF/2.0 + self.likelihoodcoef
class MCMCFollowUpSearch(MCMCSemiCoherentSearch): class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
...@@ -2152,8 +2152,8 @@ class MCMCTransientSearch(MCMCSearch): ...@@ -2152,8 +2152,8 @@ class MCMCTransientSearch(MCMCSearch):
in_theta[1] = in_theta[0] + in_theta[1] in_theta[1] = in_theta[0] + in_theta[1]
if in_theta[1] > self.maxStartTime: if in_theta[1] > self.maxStartTime:
return -np.inf return -np.inf
FS = search.run_computefstatistic_single_point(*in_theta) twoF = search.get_fullycoherent_twoF(*in_theta)
return FS + self.likelihoodcoef return twoF/2.0 + self.likelihoodcoef
def _unpack_input_theta(self): def _unpack_input_theta(self):
full_theta_keys = ['transient_tstart', full_theta_keys = ['transient_tstart',
......
...@@ -97,13 +97,9 @@ class TestComputeFstat(Test): ...@@ -97,13 +97,9 @@ class TestComputeFstat(Test):
search = pyfstat.ComputeFstat( search = pyfstat.ComputeFstat(
tref=Writer.tref, tref=Writer.tref,
sftfilepattern='{}/*{}*sft'.format(Writer.outdir, Writer.label)) sftfilepattern='{}/*{}*sft'.format(Writer.outdir, Writer.label))
FS = search.run_computefstatistic_single_point(Writer.tstart, FS = search.get_fullycoherent_twoF(
Writer.tend, Writer.tstart, Writer.tend, Writer.F0, Writer.F1, Writer.F2,
Writer.F0, Writer.Alpha, Writer.Delta)
Writer.F1,
Writer.F2,
Writer.Alpha,
Writer.Delta)
print predicted_FS, FS print predicted_FS, FS
self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2) self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2)
...@@ -116,13 +112,9 @@ class TestComputeFstat(Test): ...@@ -116,13 +112,9 @@ class TestComputeFstat(Test):
search = pyfstat.ComputeFstat( search = pyfstat.ComputeFstat(
tref=Writer.tref, assumeSqrtSX=1, tref=Writer.tref, assumeSqrtSX=1,
sftfilepattern='{}/*{}*sft'.format(Writer.outdir, Writer.label)) sftfilepattern='{}/*{}*sft'.format(Writer.outdir, Writer.label))
FS = search.run_computefstatistic_single_point(Writer.tstart, FS = search.get_fullycoherent_twoF(
Writer.tend, Writer.tstart, Writer.tend, Writer.F0, Writer.F1, Writer.F2,
Writer.F0, Writer.Alpha, Writer.Delta)
Writer.F1,
Writer.F2,
Writer.Alpha,
Writer.Delta)
print predicted_FS, FS print predicted_FS, FS
self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2) self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2)
...@@ -137,13 +129,9 @@ class TestComputeFstat(Test): ...@@ -137,13 +129,9 @@ class TestComputeFstat(Test):
minCoverFreq=28, maxCoverFreq=32, minStartTime=Writer.tstart, minCoverFreq=28, maxCoverFreq=32, minStartTime=Writer.tstart,
maxStartTime=Writer.tstart+Writer.duration, maxStartTime=Writer.tstart+Writer.duration,
detectors=Writer.detectors) detectors=Writer.detectors)
FS = search.run_computefstatistic_single_point(Writer.tstart, FS = search.get_fullycoherent_twoF(
Writer.tend, Writer.tstart, Writer.tend, Writer.F0, Writer.F1, Writer.F2,
Writer.F0, Writer.Alpha, Writer.Delta)
Writer.F1,
Writer.F2,
Writer.Alpha,
Writer.Delta)
Writer.make_data() Writer.make_data()
predicted_FS = Writer.predict_fstat() predicted_FS = Writer.predict_fstat()
print predicted_FS, FS print predicted_FS, FS
...@@ -153,7 +141,7 @@ class TestComputeFstat(Test): ...@@ -153,7 +141,7 @@ class TestComputeFstat(Test):
class TestSemiCoherentGlitchSearch(Test): class TestSemiCoherentGlitchSearch(Test):
label = "TestSemiCoherentGlitchSearch" label = "TestSemiCoherentGlitchSearch"
def test_compute_nglitch_fstat(self): def test_get_semicoherent_nglitch_twoF(self):
duration = 10*86400 duration = 10*86400
dtglitch = .5*duration dtglitch = .5*duration
delta_F0 = 0 delta_F0 = 0
...@@ -171,10 +159,9 @@ class TestSemiCoherentGlitchSearch(Test): ...@@ -171,10 +159,9 @@ class TestSemiCoherentGlitchSearch(Test):
tref=Writer.tref, minStartTime=Writer.tstart, tref=Writer.tref, minStartTime=Writer.tstart,
maxStartTime=Writer.tend, nglitch=1) maxStartTime=Writer.tend, nglitch=1)
FS = search.compute_nglitch_fstat(Writer.F0, Writer.F1, Writer.F2, FS = search.get_semicoherent_nglitch_twoF(
Writer.Alpha, Writer.Delta, Writer.F0, Writer.F1, Writer.F2, Writer.Alpha, Writer.Delta,
Writer.delta_F0, Writer.delta_F1, Writer.delta_F0, Writer.delta_F1, search.minStartTime+dtglitch)
search.minStartTime+dtglitch)
# Compute the predicted semi-coherent glitch Fstat # Compute the predicted semi-coherent glitch Fstat
minStartTime = Writer.tstart minStartTime = Writer.tstart
......
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