Commit 4153cc00 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Various minor improvements:

1) Adds BSGL labels to cumulative plot (if used)
2) Switches on transient and re-init before computing cumulative
3) Adds cumulative plot to MCMC search
4) Adds cumulative plot for MCMCGlitchSearch
5) Adds support for neghalfnorm in plot prior posterior
parent 9cb7c9f0
...@@ -266,7 +266,6 @@ class ComputeFstat(object): ...@@ -266,7 +266,6 @@ class ComputeFstat(object):
'Loaded {} data files from detectors {} spanning {} to {}'.format( 'Loaded {} data files from detectors {} spanning {} to {}'.format(
len(SFT_timestamps), names, int(SFT_timestamps[0]), len(SFT_timestamps), names, int(SFT_timestamps[0]),
int(SFT_timestamps[-1]))) int(SFT_timestamps[-1])))
self.SFT_timestamps = SFT_timestamps
logging.info('Initialising ephems') logging.info('Initialising ephems')
ephems = lalpulsar.InitBarycenter(self.earth_ephem, self.sun_ephem) ephems = lalpulsar.InitBarycenter(self.earth_ephem, self.sun_ephem)
...@@ -418,6 +417,9 @@ class ComputeFstat(object): ...@@ -418,6 +417,9 @@ class ComputeFstat(object):
duration = tend - tstart duration = tend - tstart
taus = np.linspace(minfraction*duration, duration, Npoints) taus = np.linspace(minfraction*duration, duration, Npoints)
twoFs = [] twoFs = []
if self.transient is False:
self.transient = True
self.init_computefstatistic_single_point()
for tau in taus: for tau in taus:
twoFs.append(self.run_computefstatistic_single_point( twoFs.append(self.run_computefstatistic_single_point(
tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2, tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2,
...@@ -434,11 +436,15 @@ class ComputeFstat(object): ...@@ -434,11 +436,15 @@ class ComputeFstat(object):
ax.plot(taus/86400., twoFs, label=label, color=c) ax.plot(taus/86400., twoFs, label=label, color=c)
ax.set_xlabel(r'Days from $t_{{\rm start}}={:.0f}$'.format( ax.set_xlabel(r'Days from $t_{{\rm start}}={:.0f}$'.format(
kwargs['tstart'])) kwargs['tstart']))
ax.set_ylabel(r'$\widetilde{2\mathcal{F}}_{\rm cumulative}$') if self.BSGL:
ax.set_ylabel(r'$\log_{10}(\mathrm{BSGL})_{\rm cumulative}$')
else:
ax.set_ylabel(r'$\widetilde{2\mathcal{F}}_{\rm cumulative}$')
ax.set_xlim(0, taus[-1]/86400) ax.set_xlim(0, taus[-1]/86400)
ax.set_title(title) ax.set_title(title)
if savefig: if savefig:
plt.savefig('{}/{}_twoFcumulative.png'.format(outdir, label)) plt.savefig('{}/{}_twoFcumulative.png'.format(outdir, label))
return taus, twoFs
else: else:
return ax return ax
...@@ -923,6 +929,11 @@ class MCMCSearch(BaseSearchClass): ...@@ -923,6 +929,11 @@ class MCMCSearch(BaseSearchClass):
upper = prior_dict['loc'] + normal_stds * prior_dict['scale'] upper = prior_dict['loc'] + normal_stds * prior_dict['scale']
x = np.linspace(lower, upper, N) x = np.linspace(lower, upper, N)
prior = [prior_func(xi) for xi in x] prior = [prior_func(xi) for xi in x]
elif prior_dict['type'] == 'neghalfnorm':
upper = prior_dict['loc']
lower = prior_dict['loc'] - normal_stds * prior_dict['scale']
x = np.linspace(lower, upper, N)
prior = [prior_func(xi) for xi in x]
else: else:
raise ValueError('Not implemented for prior type {}'.format( raise ValueError('Not implemented for prior type {}'.format(
prior_dict['type'])) prior_dict['type']))
...@@ -946,6 +957,16 @@ class MCMCSearch(BaseSearchClass): ...@@ -946,6 +957,16 @@ class MCMCSearch(BaseSearchClass):
fig.savefig('{}/{}_prior_posterior.png'.format( fig.savefig('{}/{}_prior_posterior.png'.format(
self.outdir, self.label)) self.outdir, self.label))
def plot_cumulative_max(self):
d, maxtwoF = self.get_max_twoF()
for key, val in self.theta_prior.iteritems():
if key not in d:
d[key] = val
self.search.plot_twoF_cumulative(
self.label, self.outdir, F0=d['F0'], F1=d['F1'], F2=d['F2'],
Alpha=d['Alpha'], Delta=d['Delta'], tstart=self.tstart,
tend=self.tend)
def generic_lnprior(self, **kwargs): def generic_lnprior(self, **kwargs):
""" Return a lambda function of the pdf """ Return a lambda function of the pdf
...@@ -1538,6 +1559,46 @@ _ sftfilepath: str ...@@ -1538,6 +1559,46 @@ _ sftfilepath: str
axis=2) axis=2)
return p0 return p0
def plot_cumulative_max(self):
fig, ax = plt.subplots()
d, maxtwoF = self.get_max_twoF()
for key, val in self.theta_prior.iteritems():
if key not in d:
d[key] = val
delta_F0s = [d['delta_F0_{}'.format(i)] for i in range(self.nglitch)]
delta_F0s.insert(self.theta0_idx, 0)
delta_F0s = np.array(delta_F0s)
delta_F0s[:self.theta0_idx] *= -1
tglitches = [d['tglitch_{}'.format(i)] for i in range(self.nglitch)]
tbounderies = [self.tstart] + tglitches + [self.tend]
for j in range(self.nglitch+1):
ts = tbounderies[j]
te = tbounderies[j+1]
if (te - ts)/86400 < 5:
logging.info('Period too short to perform cumulative search')
continue
if j < self.theta0_idx:
summed_deltaF0 = np.sum(delta_F0s[j:self.theta0_idx])
F0_j = d['F0'] - summed_deltaF0
taus, twoFs = self.search.calculate_twoF_cumulative(
F0_j, F1=d['F1'], F2=d['F2'], Alpha=d['Alpha'],
Delta=d['Delta'], tstart=ts, tend=te)
elif j >= self.theta0_idx:
summed_deltaF0 = np.sum(delta_F0s[self.theta0_idx:j+1])
F0_j = d['F0'] + summed_deltaF0
taus, twoFs = self.search.calculate_twoF_cumulative(
F0_j, F1=d['F1'], F2=d['F2'], Alpha=d['Alpha'],
Delta=d['Delta'], tstart=ts, tend=te)
ax.plot(ts+taus, twoFs)
ax.set_xlabel('GPS time')
fig.savefig('{}/{}_twoFcumulative.png'.format(self.outdir, self.label))
class GridSearch(BaseSearchClass): class GridSearch(BaseSearchClass):
""" Gridded search using ComputeFstat """ """ Gridded search using ComputeFstat """
......
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