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

Clean up of the earth/sun ephemeris handling 

- Provide a function in the BaseSearchClass to read in the ephemeris
  files
- Some clean up of docstrings
parent df586f0a
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
pyfstat/core.py
+ 36
− 46
View file @ 03a4bddd
...@@ -5,7 +5,6 @@ import os ...@@ -5,7 +5,6 @@ import os
import logging import logging
import copy import copy
import glob import glob
import numpy as np import numpy as np
import scipy.special import scipy.special
...@@ -27,15 +26,14 @@ else: ...@@ -27,15 +26,14 @@ else:
helper_functions.set_up_matplotlib_defaults() helper_functions.set_up_matplotlib_defaults()
args, tqdm = helper_functions.set_up_command_line_arguments() args, tqdm = helper_functions.set_up_command_line_arguments()
earth_ephem, sun_ephem = helper_functions.set_up_ephemeris_configuration()
detector_colors = {'h1': 'C0', 'l1': 'C1'} detector_colors = {'h1': 'C0', 'l1': 'C1'}
class Bunch(object): class Bunch(object):
""" Turns dictionary into object with attribute-style access """ Turns dictionary into object with attribute-style access
Parameter Parameters
--------- ----------
dict dict
Input dictionary Input dictionary
...@@ -178,18 +176,7 @@ def predict_fstat(h0, cosi, psi, Alpha, Delta, Freq, sftfilepattern, ...@@ -178,18 +176,7 @@ def predict_fstat(h0, cosi, psi, Alpha, Delta, Freq, sftfilepattern,
class BaseSearchClass(object): class BaseSearchClass(object):
""" The base search class providing parent methods to other searches """ The base search class providing parent methods to other searches """
Attributes
----------
earth_ephem_default, sun_ephem_default : str
Paths to the earth and sun ephemeris files to use if not specified
elsewhere
"""
earth_ephem_default = earth_ephem
sun_ephem_default = sun_ephem
def _add_log_file(self): def _add_log_file(self):
""" Log output to a file, requires class to have outdir and label """ """ Log output to a file, requires class to have outdir and label """
...@@ -306,19 +293,37 @@ class BaseSearchClass(object): ...@@ -306,19 +293,37 @@ class BaseSearchClass(object):
raise IOError('No sfts found matching {}'.format( raise IOError('No sfts found matching {}'.format(
self.sftfilepattern)) self.sftfilepattern))
def set_ephemeris_files(self, earth_ephem=None, sun_ephem=None):
""" Set the ephemeris files to use for the Earth and Sun
class ComputeFstat(object): Parameters
""" Base class providing interface to `lalpulsar.ComputeFstat` """ ----------
earth_ephem, sun_ephem: str
Paths of the two files containing positions of Earth and Sun,
respectively at evenly spaced times, as passed to CreateFstatInput
earth_ephem_default = earth_ephem Note: If not manually set, default values in ~/.pyfstat are used
sun_ephem_default = sun_ephem
"""
earth_ephem_default, sun_ephem_default = (
helper_functions.get_ephemeris_files())
if earth_ephem is None:
self.earth_ephem = earth_ephem_default
if sun_ephem is None:
self.sun_ephem = sun_ephem_default
class ComputeFstat(BaseSearchClass):
""" Base class providing interface to `lalpulsar.ComputeFstat` """
@helper_functions.initializer @helper_functions.initializer
def __init__(self, tref, sftfilepattern=None, minStartTime=None, def __init__(self, tref, sftfilepattern=None, minStartTime=None,
maxStartTime=None, binary=False, transient=True, BSGL=False, maxStartTime=None, binary=False, transient=True, BSGL=False,
detectors=None, minCoverFreq=None, maxCoverFreq=None, detectors=None, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, injectSources=None, injectSources=None, injectSqrtSX=None, assumeSqrtSX=None,
injectSqrtSX=None, assumeSqrtSX=None, SSBprec=None): SSBprec=None):
""" """
Parameters Parameters
---------- ----------
...@@ -343,10 +348,6 @@ class ComputeFstat(object): ...@@ -343,10 +348,6 @@ class ComputeFstat(object):
The min and max cover frequency passed to CreateFstatInput, if The min and max cover frequency passed to CreateFstatInput, if
either is None the range of frequencies in the SFT less 1Hz is either is None the range of frequencies in the SFT less 1Hz is
used. used.
earth_ephem, sun_ephem : str
Paths of the two files containing positions of Earth and Sun,
respectively at evenly spaced times, as passed to CreateFstatInput.
If None defaults will be used.
injectSources : dict or str injectSources : dict or str
Either a dictionary of the values to inject, or a string pointing Either a dictionary of the values to inject, or a string pointing
to the .cff file to inject to the .cff file to inject
...@@ -362,11 +363,7 @@ class ComputeFstat(object): ...@@ -362,11 +363,7 @@ class ComputeFstat(object):
""" """
if earth_ephem is None: self.set_ephemeris_files()
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
def _get_SFTCatalog(self): def _get_SFTCatalog(self):
...@@ -827,15 +824,15 @@ class ComputeFstat(object): ...@@ -827,15 +824,15 @@ class ComputeFstat(object):
return ax return ax
class SemiCoherentSearch(BaseSearchClass, ComputeFstat): class SemiCoherentSearch(ComputeFstat):
""" A semi-coherent search """ """ A semi-coherent search """
@helper_functions.initializer @helper_functions.initializer
def __init__(self, label, outdir, tref, nsegs=None, sftfilepattern=None, def __init__(self, label, outdir, tref, nsegs=None, sftfilepattern=None,
binary=False, BSGL=False, minStartTime=None, binary=False, BSGL=False, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
detectors=None, earth_ephem=None, sun_ephem=None, detectors=None, injectSources=None, assumeSqrtSX=None,
injectSources=None, assumeSqrtSX=None, SSBprec=None): SSBprec=None):
""" """
Parameters Parameters
---------- ----------
...@@ -853,10 +850,7 @@ class SemiCoherentSearch(BaseSearchClass, ComputeFstat): ...@@ -853,10 +850,7 @@ class SemiCoherentSearch(BaseSearchClass, ComputeFstat):
""" """
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
if self.earth_ephem is None: self.set_ephemeris_files()
self.earth_ephem = self.earth_ephem_default
if self.sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
self.transient = True self.transient = True
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
self.init_semicoherent_parameters() self.init_semicoherent_parameters()
...@@ -957,7 +951,7 @@ class SemiCoherentSearch(BaseSearchClass, ComputeFstat): ...@@ -957,7 +951,7 @@ class SemiCoherentSearch(BaseSearchClass, ComputeFstat):
return d_detStat return d_detStat
class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): class SemiCoherentGlitchSearch(ComputeFstat):
""" A semi-coherent glitch search """ A semi-coherent glitch search
This implements a basic `semi-coherent glitch F-stat in which the data This implements a basic `semi-coherent glitch F-stat in which the data
...@@ -970,8 +964,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): ...@@ -970,8 +964,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
def __init__(self, label, outdir, tref, minStartTime, maxStartTime, def __init__(self, label, outdir, tref, minStartTime, maxStartTime,
nglitch=0, sftfilepattern=None, theta0_idx=0, BSGL=False, nglitch=0, sftfilepattern=None, theta0_idx=0, BSGL=False,
minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None, minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None,
detectors=None, earth_ephem=None, sun_ephem=None, detectors=None, SSBprec=None, injectSources=None):
SSBprec=None, injectSources=None):
""" """
Parameters Parameters
---------- ----------
...@@ -994,10 +987,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): ...@@ -994,10 +987,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
""" """
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
if self.earth_ephem is None: self.set_ephemeris_files()
self.earth_ephem = self.earth_ephem_default
if self.sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
self.transient = True self.transient = True
self.binary = False self.binary = False
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
......
pyfstat/grid_based_searches.py
+ 6
− 32
View file @ 03a4bddd
...@@ -13,7 +13,7 @@ import matplotlib.pyplot as plt ...@@ -13,7 +13,7 @@ import matplotlib.pyplot as plt
import pyfstat.helper_functions as helper_functions import pyfstat.helper_functions as helper_functions
from pyfstat.core import (BaseSearchClass, ComputeFstat, from pyfstat.core import (BaseSearchClass, ComputeFstat,
SemiCoherentGlitchSearch, SemiCoherentSearch, tqdm, SemiCoherentGlitchSearch, SemiCoherentSearch, tqdm,
args, earth_ephem, sun_ephem, read_par) args, read_par)
class GridSearch(BaseSearchClass): class GridSearch(BaseSearchClass):
...@@ -22,9 +22,8 @@ class GridSearch(BaseSearchClass): ...@@ -22,9 +22,8 @@ class GridSearch(BaseSearchClass):
def __init__(self, label, outdir, sftfilepattern, F0s=[0], F1s=[0], F2s=[0], def __init__(self, label, outdir, sftfilepattern, F0s=[0], F1s=[0], F2s=[0],
Alphas=[0], Deltas=[0], tref=None, minStartTime=None, Alphas=[0], Deltas=[0], tref=None, minStartTime=None,
maxStartTime=None, nsegs=1, BSGL=False, minCoverFreq=None, maxStartTime=None, nsegs=1, BSGL=False, minCoverFreq=None,
maxCoverFreq=None, earth_ephem=None, sun_ephem=None, maxCoverFreq=None, detectors=None, SSBprec=None,
detectors=None, SSBprec=None, injectSources=None, injectSources=None, input_arrays=False, assumeSqrtSX=None):
input_arrays=False, assumeSqrtSX=None):
""" """
Parameters Parameters
---------- ----------
...@@ -45,11 +44,6 @@ class GridSearch(BaseSearchClass): ...@@ -45,11 +44,6 @@ class GridSearch(BaseSearchClass):
For all other parameters, see `pyfstat.ComputeFStat` for details For all other parameters, see `pyfstat.ComputeFStat` for details
""" """
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
self.set_out_file() self.set_out_file()
...@@ -61,7 +55,6 @@ class GridSearch(BaseSearchClass): ...@@ -61,7 +55,6 @@ class GridSearch(BaseSearchClass):
self.search = ComputeFstat( self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
detectors=self.detectors, transient=False, detectors=self.detectors, transient=False,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
...@@ -75,7 +68,6 @@ class GridSearch(BaseSearchClass): ...@@ -75,7 +68,6 @@ class GridSearch(BaseSearchClass):
BSGL=self.BSGL, minStartTime=self.minStartTime, BSGL=self.BSGL, minStartTime=self.minStartTime,
maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq, maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, detectors=self.detectors, maxCoverFreq=self.maxCoverFreq, detectors=self.detectors,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
injectSources=self.injectSources) injectSources=self.injectSources)
def cut_out_tstart_tend(*vals): def cut_out_tstart_tend(*vals):
...@@ -313,7 +305,7 @@ class GridGlitchSearch(GridSearch): ...@@ -313,7 +305,7 @@ class GridGlitchSearch(GridSearch):
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, minStartTime=None, Alphas=[0], Deltas=[0], tref=None, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
write_after=1000, earth_ephem=None, sun_ephem=None): write_after=1000):
""" """
Parameters Parameters
...@@ -333,16 +325,11 @@ class GridGlitchSearch(GridSearch): ...@@ -333,16 +325,11 @@ class GridGlitchSearch(GridSearch):
""" """
if tglitchs is None: if tglitchs is None:
self.tglitchs = [self.maxStartTime] self.tglitchs = [self.maxStartTime]
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
self.search = SemiCoherentGlitchSearch( self.search = SemiCoherentGlitchSearch(
label=label, outdir=outdir, sftfilepattern=self.sftfilepattern, label=label, outdir=outdir, sftfilepattern=self.sftfilepattern,
tref=tref, minStartTime=minStartTime, maxStartTime=maxStartTime, tref=tref, minStartTime=minStartTime, maxStartTime=maxStartTime,
minCoverFreq=minCoverFreq, maxCoverFreq=maxCoverFreq, minCoverFreq=minCoverFreq, maxCoverFreq=maxCoverFreq,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
BSGL=self.BSGL) BSGL=self.BSGL)
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
...@@ -372,8 +359,7 @@ class FrequencySlidingWindow(GridSearch): ...@@ -372,8 +359,7 @@ class FrequencySlidingWindow(GridSearch):
Alpha, Delta, tref, minStartTime=None, Alpha, Delta, tref, minStartTime=None,
maxStartTime=None, window_size=10*86400, window_delta=86400, maxStartTime=None, window_size=10*86400, window_delta=86400,
BSGL=False, minCoverFreq=None, maxCoverFreq=None, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, detectors=None, detectors=None, SSBprec=None, injectSources=None):
SSBprec=None, injectSources=None):
""" """
Parameters Parameters
---------- ----------
...@@ -392,11 +378,6 @@ class FrequencySlidingWindow(GridSearch): ...@@ -392,11 +378,6 @@ class FrequencySlidingWindow(GridSearch):
For all other parameters, see `pyfstat.ComputeFStat` for details For all other parameters, see `pyfstat.ComputeFStat` for details
""" """
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
self.set_out_file() self.set_out_file()
...@@ -411,7 +392,6 @@ class FrequencySlidingWindow(GridSearch): ...@@ -411,7 +392,6 @@ class FrequencySlidingWindow(GridSearch):
self.search = ComputeFstat( self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
detectors=self.detectors, transient=True, detectors=self.detectors, transient=True,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
...@@ -487,8 +467,7 @@ class DMoff_NO_SPIN(GridSearch): ...@@ -487,8 +467,7 @@ class DMoff_NO_SPIN(GridSearch):
@helper_functions.initializer @helper_functions.initializer
def __init__(self, par, label, outdir, sftfilepattern, minStartTime=None, def __init__(self, par, label, outdir, sftfilepattern, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, detectors=None, detectors=None, injectSources=None, assumeSqrtSX=None):
injectSources=None, assumeSqrtSX=None):
""" """
Parameters Parameters
---------- ----------
...@@ -506,11 +485,6 @@ class DMoff_NO_SPIN(GridSearch): ...@@ -506,11 +485,6 @@ class DMoff_NO_SPIN(GridSearch):
For all other parameters, see `pyfstat.ComputeFStat` for details For all other parameters, see `pyfstat.ComputeFStat` for details
""" """
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
......
pyfstat/helper_functions.py
+ 1
− 1
View file @ 03a4bddd
...@@ -80,7 +80,7 @@ def set_up_command_line_arguments(): ...@@ -80,7 +80,7 @@ def set_up_command_line_arguments():
return args, tqdm return args, tqdm
def set_up_ephemeris_configuration(): def get_ephemeris_files():
""" Returns the earth_ephem and sun_ephem """ """ Returns the earth_ephem and sun_ephem """
config_file = os.path.expanduser('~')+'/.pyfstat.conf' config_file = os.path.expanduser('~')+'/.pyfstat.conf'
if os.path.isfile(config_file): if os.path.isfile(config_file):
......
pyfstat/make_sfts.py
+ 3
− 12
View file @ 03a4bddd
...@@ -13,8 +13,6 @@ import lalpulsar ...@@ -13,8 +13,6 @@ import lalpulsar
from pyfstat.core import BaseSearchClass, tqdm, args from pyfstat.core import BaseSearchClass, tqdm, args
import pyfstat.helper_functions as helper_functions import pyfstat.helper_functions as helper_functions
earth_ephem, sun_ephem = helper_functions.set_up_ephemeris_configuration()
class KeyboardInterruptError(Exception): class KeyboardInterruptError(Exception):
pass pass
...@@ -49,6 +47,7 @@ class Writer(BaseSearchClass): ...@@ -49,6 +47,7 @@ class Writer(BaseSearchClass):
see `lalapps_Makefakedata_v5 --help` for help with the other paramaters see `lalapps_Makefakedata_v5 --help` for help with the other paramaters
""" """
self.set_ephemeris_files()
self.tstart = int(tstart) self.tstart = int(tstart)
self.duration = int(duration) self.duration = int(duration)
...@@ -372,15 +371,11 @@ class GlitchWriter(Writer): ...@@ -372,15 +371,11 @@ class GlitchWriter(Writer):
class FrequencyModulatedArtifactWriter(Writer): class FrequencyModulatedArtifactWriter(Writer):
""" Instance object for generating SFTs containing artifacts """ """ Instance object for generating SFTs containing artifacts """
earth_ephem_default = earth_ephem
sun_ephem_default = sun_ephem
@helper_functions.initializer @helper_functions.initializer
def __init__(self, label, outdir=".", tstart=700000000, def __init__(self, label, outdir=".", tstart=700000000,
data_duration=86400, F0=30, F1=0, tref=None, h0=10, Tsft=1800, data_duration=86400, F0=30, F1=0, tref=None, h0=10, Tsft=1800,
sqrtSX=1, Band=4, Pmod=lal.DAYSID_SI, Pmod_phi=0, Pmod_amp=1, sqrtSX=1, Band=4, Pmod=lal.DAYSID_SI, Pmod_phi=0, Pmod_amp=1,
Alpha=None, Delta=None, IFO='H1', earth_ephem=None, Alpha=None, Delta=None, IFO='H1'):
sun_ephem=None):
""" """
Parameters Parameters
---------- ----------
...@@ -399,6 +394,7 @@ class FrequencyModulatedArtifactWriter(Writer): ...@@ -399,6 +394,7 @@ class FrequencyModulatedArtifactWriter(Writer):
see `lalapps_Makefakedata_v4 --help` for help with the other paramaters see `lalapps_Makefakedata_v4 --help` for help with the other paramaters
""" """
self.set_ephemeris_files()
self.tstart = int(tstart) self.tstart = int(tstart)
self.data_duration = int(data_duration) self.data_duration = int(data_duration)
...@@ -414,11 +410,6 @@ class FrequencyModulatedArtifactWriter(Writer): ...@@ -414,11 +410,6 @@ class FrequencyModulatedArtifactWriter(Writer):
self.cosi = 0 self.cosi = 0
self.Fmax = F0 self.Fmax = F0
if self.earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if self.sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if Alpha is not None and Delta is not None: if Alpha is not None and Delta is not None:
self.n = np.array([np.cos(Alpha)*np.cos(Delta), self.n = np.array([np.cos(Alpha)*np.cos(Delta),
np.sin(Alpha)*np.cos(Delta), np.sin(Alpha)*np.cos(Delta),
......
pyfstat/mcmc_based_searches.py
+ 20
− 49
View file @ 03a4bddd
...@@ -16,8 +16,8 @@ import corner ...@@ -16,8 +16,8 @@ import corner
import dill as pickle import dill as pickle
import pyfstat.core as core import pyfstat.core as core
from pyfstat.core import tqdm, args, earth_ephem, sun_ephem, read_par from pyfstat.core import tqdm, args, read_par
from pyfstat.optimal_setup_functions import get_Nstar_estimate, get_optimal_setup import pyfstat.optimal_setup_functions as optimal_setup_functions
import pyfstat.helper_functions as helper_functions import pyfstat.helper_functions as helper_functions
...@@ -33,17 +33,17 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -33,17 +33,17 @@ class MCMCSearch(core.BaseSearchClass):
asini='', period='s', ecc='', tp='', argp='') asini='', period='s', ecc='', tp='', argp='')
rescale_dictionary = {} rescale_dictionary = {}
@helper_functions.initializer @helper_functions.initializer
def __init__(self, label, outdir, theta_prior, tref, minStartTime, def __init__(self, label, outdir, theta_prior, tref, minStartTime,
maxStartTime, sftfilepattern=None, nsteps=[100, 100], maxStartTime, sftfilepattern=None, nsteps=[100, 100],
nwalkers=100, ntemps=1, log10temperature_min=-5, nwalkers=100, ntemps=1, log10temperature_min=-5,
theta_initial=None, scatter_val=1e-10, rhohatmax=1000, theta_initial=None, scatter_val=1e-10, rhohatmax=1000,
binary=False, BSGL=False, minCoverFreq=None, SSBprec=None, binary=False, BSGL=False, minCoverFreq=None, SSBprec=None,
maxCoverFreq=None, detectors=None, earth_ephem=None, maxCoverFreq=None, detectors=None,
sun_ephem=None, injectSources=None, assumeSqrtSX=None): injectSources=None, assumeSqrtSX=None):
""" """
Parameters Parameters
----------
label, outdir: str label, outdir: str
A label and directory to read/write data from/to A label and directory to read/write data from/to
sftfilepattern: str sftfilepattern: str
...@@ -83,10 +83,6 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -83,10 +83,6 @@ class MCMCSearch(core.BaseSearchClass):
minCoverFreq, maxCoverFreq: float minCoverFreq, maxCoverFreq: float
Minimum and maximum instantaneous frequency which will be covered Minimum and maximum instantaneous frequency which will be covered
over the SFT time span as passed to CreateFstatInput over the SFT time span as passed to CreateFstatInput
earth_ephem, sun_ephem: str
Paths of the two files containing positions of Earth and Sun,
respectively at evenly spaced times, as passed to CreateFstatInput
If None defaults defined in BaseSearchClass will be used
""" """
...@@ -108,11 +104,6 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -108,11 +104,6 @@ class MCMCSearch(core.BaseSearchClass):
else: else:
self.betas = None self.betas = None
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if args.clean and os.path.isfile(self.pickle_path): if args.clean and os.path.isfile(self.pickle_path):
os.rename(self.pickle_path, self.pickle_path+".old") os.rename(self.pickle_path, self.pickle_path+".old")
...@@ -137,7 +128,6 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -137,7 +128,6 @@ class MCMCSearch(core.BaseSearchClass):
self.search = core.ComputeFstat( self.search = core.ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
detectors=self.detectors, BSGL=self.BSGL, transient=False, detectors=self.detectors, BSGL=self.BSGL, transient=False,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
binary=self.binary, injectSources=self.injectSources, binary=self.binary, injectSources=self.injectSources,
...@@ -1516,7 +1506,7 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1516,7 +1506,7 @@ class MCMCGlitchSearch(MCMCSearch):
theta_initial=None, scatter_val=1e-10, rhohatmax=1000, theta_initial=None, scatter_val=1e-10, rhohatmax=1000,
dtglitchmin=1*86400, theta0_idx=0, detectors=None, dtglitchmin=1*86400, theta0_idx=0, detectors=None,
BSGL=False, minCoverFreq=None, maxCoverFreq=None, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, injectSources=None): injectSources=None):
""" """
Parameters Parameters
---------- ----------
...@@ -1570,10 +1560,6 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1570,10 +1560,6 @@ class MCMCGlitchSearch(MCMCSearch):
minCoverFreq, maxCoverFreq: float minCoverFreq, maxCoverFreq: float
Minimum and maximum instantaneous frequency which will be covered Minimum and maximum instantaneous frequency which will be covered
over the SFT time span as passed to CreateFstatInput over the SFT time span as passed to CreateFstatInput
earth_ephem, sun_ephem: str
Paths of the two files containing positions of Earth and Sun,
respectively at evenly spaced times, as passed to CreateFstatInput
If None defaults defined in BaseSearchClass will be used
""" """
...@@ -1590,11 +1576,6 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1590,11 +1576,6 @@ class MCMCGlitchSearch(MCMCSearch):
self.betas = np.logspace(0, self.log10temperature_min, self.ntemps) self.betas = np.logspace(0, self.log10temperature_min, self.ntemps)
else: else:
self.betas = None self.betas = None
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if args.clean and os.path.isfile(self.pickle_path): if args.clean and os.path.isfile(self.pickle_path):
os.rename(self.pickle_path, self.pickle_path+".old") os.rename(self.pickle_path, self.pickle_path+".old")
...@@ -1611,8 +1592,7 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1611,8 +1592,7 @@ class MCMCGlitchSearch(MCMCSearch):
label=self.label, outdir=self.outdir, sftfilepattern=self.sftfilepattern, label=self.label, outdir=self.outdir, sftfilepattern=self.sftfilepattern,
tref=self.tref, minStartTime=self.minStartTime, tref=self.tref, minStartTime=self.minStartTime,
maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq, maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem, maxCoverFreq=self.maxCoverFreq, detectors=self.detectors, BSGL=self.BSGL,
sun_ephem=self.sun_ephem, detectors=self.detectors, BSGL=self.BSGL,
nglitch=self.nglitch, theta0_idx=self.theta0_idx, nglitch=self.nglitch, theta0_idx=self.theta0_idx,
injectSources=self.injectSources) injectSources=self.injectSources)
...@@ -1778,8 +1758,7 @@ class MCMCSemiCoherentSearch(MCMCSearch): ...@@ -1778,8 +1758,7 @@ class MCMCSemiCoherentSearch(MCMCSearch):
theta_initial=None, scatter_val=1e-10, rhohatmax=1000, theta_initial=None, scatter_val=1e-10, rhohatmax=1000,
detectors=None, BSGL=False, minStartTime=None, detectors=None, BSGL=False, minStartTime=None,
maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, injectSources=None, injectSources=None, assumeSqrtSX=None):
assumeSqrtSX=None):
""" """
""" """
...@@ -1797,11 +1776,6 @@ class MCMCSemiCoherentSearch(MCMCSearch): ...@@ -1797,11 +1776,6 @@ class MCMCSemiCoherentSearch(MCMCSearch):
self.betas = np.logspace(0, self.log10temperature_min, self.ntemps) self.betas = np.logspace(0, self.log10temperature_min, self.ntemps)
else: else:
self.betas = None self.betas = None
if earth_ephem is None:
self.earth_ephem = self.earth_ephem_default
if sun_ephem is None:
self.sun_ephem = self.sun_ephem_default
if args.clean and os.path.isfile(self.pickle_path): if args.clean and os.path.isfile(self.pickle_path):
os.rename(self.pickle_path, self.pickle_path+".old") os.rename(self.pickle_path, self.pickle_path+".old")
...@@ -1827,11 +1801,10 @@ class MCMCSemiCoherentSearch(MCMCSearch): ...@@ -1827,11 +1801,10 @@ class MCMCSemiCoherentSearch(MCMCSearch):
logging.info('Setting up search object') logging.info('Setting up search object')
self.search = core.SemiCoherentSearch( self.search = core.SemiCoherentSearch(
label=self.label, outdir=self.outdir, tref=self.tref, label=self.label, outdir=self.outdir, tref=self.tref,
nsegs=self.nsegs, sftfilepattern=self.sftfilepattern, binary=self.binary, nsegs=self.nsegs, sftfilepattern=self.sftfilepattern,
BSGL=self.BSGL, minStartTime=self.minStartTime, binary=self.binary, BSGL=self.BSGL, minStartTime=self.minStartTime,
maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq, maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, detectors=self.detectors, maxCoverFreq=self.maxCoverFreq, detectors=self.detectors,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
injectSources=self.injectSources, assumeSqrtSX=self.assumeSqrtSX) injectSources=self.injectSources, assumeSqrtSX=self.assumeSqrtSX)
def logp(self, theta_vals, theta_prior, theta_keys, search): def logp(self, theta_vals, theta_prior, theta_keys, search):
...@@ -1888,7 +1861,7 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch): ...@@ -1888,7 +1861,7 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
return True return True
else: else:
logging.info( logging.info(
'Old setup does not match one of NstarMax, Nsegs0 or prior') "Old setup doesn't match one of NstarMax, Nsegs0 or prior")
except KeyError as e: except KeyError as e:
logging.info( logging.info(
'Error found when comparing with old setup: {}'.format(e)) 'Error found when comparing with old setup: {}'.format(e))
...@@ -1926,13 +1899,13 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch): ...@@ -1926,13 +1899,13 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
generate_setup = True generate_setup = True
if generate_setup: if generate_setup:
nsegs_vals, Nstar_vals = get_optimal_setup( nsegs_vals, Nstar_vals = (
optimal_setup_functions.get_optimal_setup(
NstarMax, Nsegs0, self.tref, self.minStartTime, NstarMax, Nsegs0, self.tref, self.minStartTime,
self.maxStartTime, self.theta_prior, self.maxStartTime, self.theta_prior,
self.search.detector_names, self.earth_ephem, self.search.detector_names))
self.sun_ephem) self.write_setup_input_file(run_setup_input_file, NstarMax,
self.write_setup_input_file(run_setup_input_file, NstarMax, Nsegs0, Nsegs0, nsegs_vals, Nstar_vals,
nsegs_vals, Nstar_vals,
self.theta_prior) self.theta_prior)
run_setup = [((self.nsteps[0], 0), nsegs, False) run_setup = [((self.nsteps[0], 0), nsegs, False)
...@@ -1954,10 +1927,9 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch): ...@@ -1954,10 +1927,9 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
if args.no_template_counting: if args.no_template_counting:
Nstar_vals.append([1, 1, 1]) Nstar_vals.append([1, 1, 1])
else: else:
Nstar = get_Nstar_estimate( Nstar = optimal_setup_functions.get_Nstar_estimate(
rs[1], self.tref, self.minStartTime, self.maxStartTime, rs[1], self.tref, self.minStartTime, self.maxStartTime,
self.theta_prior, self.search.detector_names, self.theta_prior, self.search.detector_names)
self.earth_ephem, self.sun_ephem)
Nstar_vals.append(Nstar) Nstar_vals.append(Nstar)
if log_table: if log_table:
...@@ -2137,7 +2109,6 @@ class MCMCTransientSearch(MCMCSearch): ...@@ -2137,7 +2109,6 @@ class MCMCTransientSearch(MCMCSearch):
self.search = core.ComputeFstat( self.search = core.ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
detectors=self.detectors, transient=True, detectors=self.detectors, transient=True,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, binary=self.binary, BSGL=self.BSGL, binary=self.binary,
......
pyfstat/optimal_setup_functions.py
+ 9
− 13
View file @ 03a4bddd
...@@ -15,7 +15,7 @@ import pyfstat.helper_functions as helper_functions ...@@ -15,7 +15,7 @@ import pyfstat.helper_functions as helper_functions
def get_optimal_setup( def get_optimal_setup(
NstarMax, Nsegs0, tref, minStartTime, maxStartTime, prior, NstarMax, Nsegs0, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem): detector_names):
""" Using an optimisation step, calculate the optimal setup ladder """ Using an optimisation step, calculate the optimal setup ladder
Parameters Parameters
...@@ -29,7 +29,6 @@ def get_optimal_setup( ...@@ -29,7 +29,6 @@ def get_optimal_setup(
Prior dictionary, each item must either be a fixed scalar value, or Prior dictionary, each item must either be a fixed scalar value, or
a uniform prior. a uniform prior.
detector_names : list of str detector_names : list of str
earth_ephem, sun_ephem : str
Returns Returns
------- -------
...@@ -43,7 +42,7 @@ def get_optimal_setup( ...@@ -43,7 +42,7 @@ def get_optimal_setup(
Nstar_0 = get_Nstar_estimate( Nstar_0 = get_Nstar_estimate(
Nsegs0, tref, minStartTime, maxStartTime, prior, Nsegs0, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem) detector_names)
logging.info( logging.info(
'Stage {}, nsegs={}, Nstar={}'.format(0, Nsegs0, int(Nstar_0))) 'Stage {}, nsegs={}, Nstar={}'.format(0, Nsegs0, int(Nstar_0)))
...@@ -55,7 +54,7 @@ def get_optimal_setup( ...@@ -55,7 +54,7 @@ def get_optimal_setup(
while nsegs_i > 1: while nsegs_i > 1:
nsegs_i, Nstar_i = _get_nsegs_ip1( nsegs_i, Nstar_i = _get_nsegs_ip1(
nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior, nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem) detector_names)
nsegs_vals.append(nsegs_i) nsegs_vals.append(nsegs_i)
Nstar_vals.append(Nstar_i) Nstar_vals.append(Nstar_i)
i += 1 i += 1
...@@ -66,13 +65,13 @@ def get_optimal_setup( ...@@ -66,13 +65,13 @@ def get_optimal_setup(
def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior, def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem): detector_names):
""" Calculate Nsegs_{i+1} given Nsegs_{i} """ """ Calculate Nsegs_{i+1} given Nsegs_{i} """
log10NstarMax = np.log10(NstarMax) log10NstarMax = np.log10(NstarMax)
log10Nstari = np.log10(get_Nstar_estimate( log10Nstari = np.log10(get_Nstar_estimate(
nsegs_i, tref, minStartTime, maxStartTime, prior, nsegs_i, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem)) detector_names))
def f(nsegs_ip1): def f(nsegs_ip1):
if nsegs_ip1[0] > nsegs_i: if nsegs_ip1[0] > nsegs_i:
...@@ -83,8 +82,7 @@ def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior, ...@@ -83,8 +82,7 @@ def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior,
if nsegs_ip1 == 0: if nsegs_ip1 == 0:
nsegs_ip1 = 1 nsegs_ip1 = 1
Nstarip1 = get_Nstar_estimate( Nstarip1 = get_Nstar_estimate(
nsegs_ip1, tref, minStartTime, maxStartTime, prior, nsegs_ip1, tref, minStartTime, maxStartTime, prior, detector_names)
detector_names, earth_ephem, sun_ephem)
if Nstarip1 is None: if Nstarip1 is None:
return 1e6 return 1e6
else: else:
...@@ -99,7 +97,7 @@ def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior, ...@@ -99,7 +97,7 @@ def _get_nsegs_ip1(nsegs_i, NstarMax, tref, minStartTime, maxStartTime, prior,
if res.success: if res.success:
return nsegs_ip1, get_Nstar_estimate( return nsegs_ip1, get_Nstar_estimate(
nsegs_ip1, tref, minStartTime, maxStartTime, prior, nsegs_ip1, tref, minStartTime, maxStartTime, prior,
detector_names, earth_ephem, sun_ephem) detector_names)
else: else:
raise ValueError('Optimisation unsuccesful') raise ValueError('Optimisation unsuccesful')
...@@ -160,8 +158,7 @@ def _extract_data_from_prior(prior): ...@@ -160,8 +158,7 @@ def _extract_data_from_prior(prior):
def get_Nstar_estimate( def get_Nstar_estimate(
nsegs, tref, minStartTime, maxStartTime, prior, nsegs, tref, minStartTime, maxStartTime, prior, detector_names):
detector_names, earth_ephem, sun_ephem):
""" Returns N* estimated from the super-sky metric """ Returns N* estimated from the super-sky metric
Parameters Parameters
...@@ -176,8 +173,6 @@ def get_Nstar_estimate( ...@@ -176,8 +173,6 @@ def get_Nstar_estimate(
The prior dictionary The prior dictionary
detector_names : array detector_names : array
Array of detectors to average over Array of detectors to average over
earth_ephem, sun_ephem : str
Paths to the ephemeris files
Returns Returns
------- -------
...@@ -187,6 +182,7 @@ def get_Nstar_estimate( ...@@ -187,6 +182,7 @@ def get_Nstar_estimate(
thickness is unity. thickness is unity.
""" """
earth_ephem, sun_ephem = helper_functions.get_ephemeris_files()
in_phys, spindowns, sky, fiducial_freq = _extract_data_from_prior(prior) in_phys, spindowns, sky, fiducial_freq = _extract_data_from_prior(prior)
out_rssky = np.zeros(in_phys.shape) out_rssky = np.zeros(in_phys.shape)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint