Commit 856dd6c2 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Adds binary search parameters to basic search and MCMC search

Note: this also adds a check for transient, if not called, then the
whole data is searched. This does not however mean only transient, but
also glitch like searches which use the transient tools.

Additionally fixes the plot_walkers command to allow it to plot walkers
when ndim =1 and to switch off the use of offsets
parent 8435f54d
...@@ -150,7 +150,8 @@ class ComputeFstat(object): ...@@ -150,7 +150,8 @@ class ComputeFstat(object):
@initializer @initializer
def __init__(self, tref, sftlabel=None, sftdir=None, def __init__(self, tref, sftlabel=None, sftdir=None,
minCoverFreq=None, maxCoverFreq=None, minCoverFreq=None, maxCoverFreq=None,
detector=None, earth_ephem=None, sun_ephem=None): detector=None, earth_ephem=None, sun_ephem=None,
binary=False, transient=True):
if earth_ephem is None: if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default self.earth_ephem = self.earth_ephem_default
...@@ -178,7 +179,11 @@ class ComputeFstat(object): ...@@ -178,7 +179,11 @@ class ComputeFstat(object):
logging.info('Initialising FstatInput') logging.info('Initialising FstatInput')
dFreq = 0 dFreq = 0
self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET if self.transient:
self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET
else:
self.whatToCompute = lalpulsar.FSTATQ_2F
FstatOptionalArgs = lalpulsar.FstatOptionalArgsDefaults FstatOptionalArgs = lalpulsar.FstatOptionalArgsDefaults
if self.minCoverFreq is None or self.maxCoverFreq is None: if self.minCoverFreq is None or self.maxCoverFreq is None:
...@@ -210,20 +215,29 @@ class ComputeFstat(object): ...@@ -210,20 +215,29 @@ class ComputeFstat(object):
logging.info('Initialising FstatResults') logging.info('Initialising FstatResults')
self.FstatResults = lalpulsar.FstatResults() self.FstatResults = lalpulsar.FstatResults()
self.windowRange = lalpulsar.transientWindowRange_t() if self.transient:
self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR self.windowRange = lalpulsar.transientWindowRange_t()
self.windowRange.t0Band = 0 self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR
self.windowRange.dt0 = 1 self.windowRange.t0Band = 0
self.windowRange.tauBand = 0 self.windowRange.dt0 = 1
self.windowRange.dtau = 1 self.windowRange.tauBand = 0
self.windowRange.dtau = 1
def run_computefstatistic_single_point(self, tstart, tend, F0, F1, def run_computefstatistic_single_point(self, tstart, tend, F0, F1,
F2, Alpha, Delta): F2, Alpha, Delta, asini=None,
period=None, ecc=None, tp=None,
argp=None):
""" Compute the F-stat fully-coherently at a single point """ """ Compute the F-stat 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])
self.PulsarDopplerParams.Alpha = Alpha self.PulsarDopplerParams.Alpha = Alpha
self.PulsarDopplerParams.Delta = Delta self.PulsarDopplerParams.Delta = Delta
if self.binary:
self.PulsarDopplerParams.asini = asini
self.PulsarDopplerParams.period = period
self.PulsarDopplerParams.ecc = ecc
self.PulsarDopplerParams.tp = tp
self.PulsarDopplerParams.argp = argp
lalpulsar.ComputeFstat(self.FstatResults, lalpulsar.ComputeFstat(self.FstatResults,
self.FstatInput, self.FstatInput,
...@@ -232,6 +246,9 @@ class ComputeFstat(object): ...@@ -232,6 +246,9 @@ class ComputeFstat(object):
self.whatToCompute self.whatToCompute
) )
if self.transient is False:
return self.FstatResults.twoF[0]
self.windowRange.t0 = int(tstart) # TYPE UINT4 self.windowRange.t0 = int(tstart) # TYPE UINT4
self.windowRange.tau = int(tend - tstart) # TYPE UINT4 self.windowRange.tau = int(tend - tstart) # TYPE UINT4
FS = lalpulsar.ComputeTransientFstatMap( FS = lalpulsar.ComputeTransientFstatMap(
...@@ -351,8 +368,8 @@ class MCMCSearch(BaseSearchClass): ...@@ -351,8 +368,8 @@ class MCMCSearch(BaseSearchClass):
@initializer @initializer
def __init__(self, label, outdir, sftlabel, sftdir, theta_prior, tref, def __init__(self, label, outdir, sftlabel, sftdir, theta_prior, tref,
tstart, tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1, tstart, tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1,
theta_initial=None, minCoverFreq=None, theta_initial=None, minCoverFreq=None, maxCoverFreq=None,
maxCoverFreq=None, scatter_val=1e-4, betas=None, scatter_val=1e-4, betas=None, binary=False,
detector=None, earth_ephem=None, sun_ephem=None): detector=None, earth_ephem=None, sun_ephem=None):
""" """
Parameters Parameters
...@@ -385,6 +402,8 @@ class MCMCSearch(BaseSearchClass): ...@@ -385,6 +402,8 @@ class MCMCSearch(BaseSearchClass):
Paths of the two files containing positions of Earth and Sun, Paths of the two files containing positions of Earth and Sun,
respectively at evenly spaced times, as passed to CreateFstatInput respectively at evenly spaced times, as passed to CreateFstatInput
If None defaults defined in BaseSearchClass will be used If None defaults defined in BaseSearchClass will be used
binary: Bool
If true, search over binary parameters
""" """
...@@ -415,7 +434,7 @@ class MCMCSearch(BaseSearchClass): ...@@ -415,7 +434,7 @@ class MCMCSearch(BaseSearchClass):
tref=self.tref, sftlabel=self.sftlabel, tref=self.tref, sftlabel=self.sftlabel,
sftdir=self.sftdir, minCoverFreq=self.minCoverFreq, sftdir=self.sftdir, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem, maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem,
sun_ephem=self.sun_ephem, detector=self.detector) sun_ephem=self.sun_ephem, detector=self.detector, transient=False)
def logp(self, theta_vals, theta_prior, theta_keys, search): def logp(self, theta_vals, theta_prior, theta_keys, search):
H = [self.generic_lnprior(**theta_prior[key])(p) for p, key in H = [self.generic_lnprior(**theta_prior[key])(p) for p, key in
...@@ -431,10 +450,17 @@ class MCMCSearch(BaseSearchClass): ...@@ -431,10 +450,17 @@ class MCMCSearch(BaseSearchClass):
def unpack_input_theta(self): def unpack_input_theta(self):
full_theta_keys = ['tstart', 'tend', 'F0', 'F1', 'F2', 'Alpha', full_theta_keys = ['tstart', 'tend', 'F0', 'F1', 'F2', 'Alpha',
'Delta'] 'Delta']
if self.binary:
full_theta_keys += [
'asini', 'period', 'ecc', 'tp', 'argp']
full_theta_keys_copy = copy.copy(full_theta_keys) full_theta_keys_copy = copy.copy(full_theta_keys)
full_theta_symbols = ['_', '_', '$f$', '$\dot{f}$', '$\ddot{f}$', full_theta_symbols = ['_', '_', '$f$', '$\dot{f}$', '$\ddot{f}$',
r'$\alpha$', r'$\delta$'] r'$\alpha$', r'$\delta$']
if self.binary:
full_theta_symbols += [
'asini', 'period', 'period', 'ecc', 'tp', 'argp']
self.theta_keys = [] self.theta_keys = []
fixed_theta_dict = {} fixed_theta_dict = {}
for key, val in self.theta_prior.iteritems(): for key, val in self.theta_prior.iteritems():
...@@ -685,12 +711,18 @@ class MCMCSearch(BaseSearchClass): ...@@ -685,12 +711,18 @@ class MCMCSearch(BaseSearchClass):
if ndim > 1: if ndim > 1:
for i in range(ndim): for i in range(ndim):
axes[i].plot(chain[:, start:stop, i].T, color="k", cs = chain[:, start:stop, i].T
alpha=alpha) axes[i].plot(cs, color="k", alpha=alpha)
if symbols: if symbols:
axes[i].set_ylabel(symbols[i]) axes[i].set_ylabel(symbols[i])
if draw_vline is not None: if draw_vline is not None:
axes[i].axvline(draw_vline, lw=2, ls="--") axes[i].axvline(draw_vline, lw=2, ls="--")
axes[i].ticklabel_format(useOffset=False, axis='y')
else:
cs = chain[:, start:stop, 0].T
axes.plot(cs, color='k', alpha=alpha)
axes.ticklabel_format(useOffset=False, axis='y')
return fig, axes return fig, axes
...@@ -744,7 +776,7 @@ class MCMCSearch(BaseSearchClass): ...@@ -744,7 +776,7 @@ class MCMCSearch(BaseSearchClass):
logging.warning("The sampler has produced nan's") logging.warning("The sampler has produced nan's")
p = pF[np.nanargmax(lnp)] p = pF[np.nanargmax(lnp)]
p0 = self._generate_scattered_p0(p) p0 = self.generate_scattered_p0(p)
return p0 return p0
......
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