Skip to content
Snippets Groups Projects
Commit 889b0ab3 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Various fixes to the pyfstat module

1) Fixes grid search and adds example
2) Automates min/max Cover Freq to simplify input
parent 7b14a868
No related branches found
No related tags found
No related merge requests found
from pyfstat import GridGlitchSearch, Writer
import numpy as np
phi = 0
F0 = 2.4343
F1 = -1e-10
F2 = 0
delta_phi = 0
delta_F0 = 1e-7
delta_F1 = 0
duration = 100*86400
dtglitch = 100*43200
Alpha = 5e-3
Delta = 6e-2
tstart = 700000000
tend = tstart+duration
tglitch = tstart + dtglitch
tref = tglitch
glitch_data = Writer(label='glitch', outdir='data', delta_phi=delta_phi,
delta_F0=delta_F0, tref=tref, tstart=tstart,
delta_F1=delta_F1, phi=phi, F0=F0, F1=F1, F2=F2,
duration=duration, dtglitch=dtglitch, Alpha=Alpha,
Delta=Delta)
glitch_data.make_data()
F0s = [F0]
F1s = [F1]
F2s = [F2]
m = 1e-3
dF = np.sqrt(12 * m) / (np.pi * duration)
delta_F0s = [-1e-6*F0, 1e-6*F0, dF]
delta_F1s = [delta_F1]
dT = duration / 10.
tglitchs = [tstart+dT, tend-dT, duration/100.]
Alphas = [Alpha]
Deltas = [Delta]
grid = GridGlitchSearch('grid', 'data', glitch_data.label, glitch_data.outdir,
F0s, F1s, F2s, delta_F0s, delta_F1s, tglitchs, Alphas,
Deltas, tref, tstart, tend)
grid.run()
grid.plot_2D('delta_F0', 'tglitch')
print grid.get_max_twoF()
...@@ -169,7 +169,9 @@ class ComputeFstat(object): ...@@ -169,7 +169,9 @@ class ComputeFstat(object):
self.sft_filepath = self.sftdir+'/*_'+self.sftlabel+"*sft" self.sft_filepath = self.sftdir+'/*_'+self.sftlabel+"*sft"
SFTCatalog = lalpulsar.SFTdataFind(self.sft_filepath, constraints) SFTCatalog = lalpulsar.SFTdataFind(self.sft_filepath, constraints)
names = list(set([d.header.name for d in SFTCatalog.data])) names = list(set([d.header.name for d in SFTCatalog.data]))
logging.info('Loaded data from detectors {}'.format(names)) logging.info(
'Loaded data from detectors {} matching pattern {}'.format(
names, self.sft_filepath))
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)
...@@ -185,6 +187,9 @@ class ComputeFstat(object): ...@@ -185,6 +187,9 @@ class ComputeFstat(object):
fB = fA + (numBins-1)*SFTCatalog.data[0].header.deltaF fB = fA + (numBins-1)*SFTCatalog.data[0].header.deltaF
self.minCoverFreq = fA + 0.5 self.minCoverFreq = fA + 0.5
self.maxCoverFreq = fB - 0.5 self.maxCoverFreq = fB - 0.5
logging.info('Min/max cover freqs not provided, using '
'{} and {}, est. from SFTs'.format(
self.minCoverFreq, self.maxCoverFreq))
self.FstatInput = lalpulsar.CreateFstatInput(SFTCatalog, self.FstatInput = lalpulsar.CreateFstatInput(SFTCatalog,
self.minCoverFreq, self.minCoverFreq,
...@@ -347,9 +352,9 @@ class MCMCGlitchSearch(BaseSearchClass): ...@@ -347,9 +352,9 @@ class MCMCGlitchSearch(BaseSearchClass):
@initializer @initializer
def __init__(self, label, outdir, sftlabel, sftdir, theta, tref, tstart, def __init__(self, label, outdir, sftlabel, sftdir, theta, tref, tstart,
tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1, tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1,
nglitch=0, minCoverFreq=29, maxCoverFreq=31, scatter_val=1e-4, nglitch=0, minCoverFreq=None, maxCoverFreq=None,
betas=None, detector=None, dtglitchmin=20*86400, scatter_val=1e-4, betas=None, detector=None,
earth_ephem=None, sun_ephem=None): dtglitchmin=20*86400, earth_ephem=None, sun_ephem=None):
""" """
Parameters Parameters
label, outdir: str label, outdir: str
...@@ -904,7 +909,7 @@ class GridGlitchSearch(BaseSearchClass): ...@@ -904,7 +909,7 @@ class GridGlitchSearch(BaseSearchClass):
def __init__(self, label, outdir, sftlabel=None, sftdir=None, F0s=[0], def __init__(self, label, outdir, sftlabel=None, sftdir=None, F0s=[0],
F1s=[0], F2s=[0], delta_F0s=[0], delta_F1s=[0], tglitchs=None, F1s=[0], F2s=[0], delta_F0s=[0], delta_F1s=[0], tglitchs=None,
Alphas=[0], Deltas=[0], tref=None, tstart=None, tend=None, Alphas=[0], Deltas=[0], tref=None, tstart=None, tend=None,
minCoverFreq=29, maxCoverFreq=31, write_after=1000, minCoverFreq=None, maxCoverFreq=None, write_after=1000,
earth_ephem=None, sun_ephem=None): earth_ephem=None, sun_ephem=None):
""" """
Parameters Parameters
...@@ -946,8 +951,8 @@ class GridGlitchSearch(BaseSearchClass): ...@@ -946,8 +951,8 @@ class GridGlitchSearch(BaseSearchClass):
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
self.out_file = '{}/{}_gridFS.txt'.format(self.outdir, self.label) self.out_file = '{}/{}_gridFS.txt'.format(self.outdir, self.label)
self.keys = ['F0', 'F1', 'F2', 'delta_F0', 'delta_F1', 'tglitch', self.keys = ['F0', 'F1', 'F2', 'Alpha', 'Delta', 'delta_F0',
'Alpha', 'Delta'] 'delta_F1', 'tglitch']
def get_array_from_tuple(self, x): def get_array_from_tuple(self, x):
if len(x) == 1: if len(x) == 1:
...@@ -995,7 +1000,7 @@ class GridGlitchSearch(BaseSearchClass): ...@@ -995,7 +1000,7 @@ class GridGlitchSearch(BaseSearchClass):
counter = 0 counter = 0
data = [] data = []
for vals in self.input_data: for vals in self.input_data:
FS = self.search.compute_glitch_fstat(*vals) FS = self.search.compute_glitch_fstat_single(*vals)
data.append(list(vals) + [FS]) data.append(list(vals) + [FS])
if counter > self.write_after: if counter > self.write_after:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment