pyfstat package
These pages document the full API and set of classes provided by PyFstat.
See here for installation instructions and other general information.
Submodules
pyfstat.core module
The core tools used in pyfstat
- class pyfstat.core.BaseSearchClass(*args, **kwargs)[source]
Bases:
object
The base class providing parent methods to other PyFstat classes.
This does not actually have any ‘search’ functionality, which needs to be added by child classes along with full initialization and any other custom methods.
- set_ephemeris_files(earth_ephem=None, sun_ephem=None)[source]
Set the ephemeris files to use for the Earth and Sun.
NOTE: If not given explicit arguments, default values from helper_functions.get_ephemeris_files() are used.
- Parameters
earth_ephem (str) – Paths of the two files containing positions of Earth and Sun, respectively at evenly spaced times, as passed to CreateFstatInput
sun_ephem (str) – Paths of the two files containing positions of Earth and Sun, respectively at evenly spaced times, as passed to CreateFstatInput
- pprint_init_params_dict()[source]
Pretty-print a parameters dictionary for output file headers.
- Returns
pretty_init_parameters – A list of lines to be printed, including opening/closing “{” and “}”, consistent indentation, as well as end-of-line commas, but no comment markers at start of lines.
- Return type
list
- get_output_file_header()[source]
Constructs a meta-information header for text output files.
This will include PyFstat and LALSuite versioning, information about when/where/how the code was run, and input parameters of the instantiated class.
- Returns
header – A list of formatted header lines.
- Return type
list
- read_par(filename=None, label=None, outdir=None, suffix='par', raise_error=True)[source]
Read a key=val file and return a dictionary.
- Parameters
filename (str or None) – Filename (path) containing rows of key=val data to read in.
label (str or None) – If filename is None, form the file to read as outdir/label.suffix.
outdir (str or None) – If filename is None, form the file to read as outdir/label.suffix.
suffix (str or None) – If filename is None, form the file to read as outdir/label.suffix.
raise_error (bool) – If True, raise an error for lines which are not comments, but cannot be read.
- Returns
params_dict – A dictionary of the parsed key=val pairs.
- Return type
dict
- static translate_keys_to_lal(dictionary)[source]
Convert input keys into lalpulsar convention.
In PyFstat’s convention, input keys (search parameter names) are F0, F1, F2, …, while lalpulsar functions prefer to use Freq, f1dot, f2dot, ….
Since lalpulsar keys are only used internally to call lalpulsar routines, this function is provided so the keys can be translated on the fly.
- Parameters
dictionary (dict) – Dictionary to translate. A copy will be made (and returned) before translation takes place.
- Returns
translated_dict – Copy of “dictionary” with new keys according to lalpulsar convention.
- Return type
dict
- class pyfstat.core.ComputeFstat(*args, **kwargs)[source]
Bases:
pyfstat.core.BaseSearchClass
Base search class providing an interface to lalpulsar.ComputeFstat.
In most cases, users should be using one of the higher-level search classes from the grid_based_searches or mcmc_based_searches modules instead.
See the lalpulsar documentation at https://lscsoft.docs.ligo.org/lalsuite/lalpulsar/group___compute_fstat__h.html and R. Prix, The F-statistic and its implementation in ComputeFstatistic_v2 ( https://dcc.ligo.org/T0900149/public ) for details of the lalpulsar module and the meaning of various technical concepts as embodied by some of the class’s parameters.
Normally this will read in existing data through the sftfilepattern argument, but if that option is None and the necessary alternative arguments are used, it can also generate simulated data (including noise and/or signals) on the fly.
NOTE that the detection statistics that can be computed from an instance of this class depend on the BSGL, BtSG and transientWindowType arguments given at initialisation. See get_fullycoherent_detstat() and get_transient_detstats() for details. To change what you want to compute, you may need to initialise a new instance with different options.
NOTE for GPU users (tCWFstatMapVersion=”pycuda”): This class tries to conveniently deal with GPU context management behind the scenes. A known problematic case is if you try to instantiate it twice from the same session/script. If you then get some messages like RuntimeError: make_default_context() and invalid device context, that is because the GPU is still blocked from the first instance when you try to initiate the second. To avoid this problem, use context management:
with pyfstat.ComputeFstat( [...], tCWFstatMapVersion="pycuda", ) as search: search.get_fullycoherent_detstat([...])
or manually call the search.finalizer_() method where needed.
- Parameters
tref (int) – GPS seconds of the reference time.
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.
minStartTime (int) – Only use SFTs with timestamps starting from within this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime].
maxStartTime (int) – Only use SFTs with timestamps starting from within this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime].
Tsft (int) – SFT duration in seconds. Only required if sftfilepattern=None and hence simulted data is generated on the fly.
binary (bool) – If true, search over binary parameters.
singleFstats (bool) – If true, also compute the single-detector twoF values.
BSGL (bool) –
If true, compute the log10BSGL statistic rather than the twoF value. For details, see Keitel et al (PRD 89, 064023, 2014): https://arxiv.org/abs/1311.5738 Note this automatically sets singleFstats=True as well. Tuning parameters are currently hardcoded:
Fstar0=15 for coherent searches.
A p-value of 1e-6 and correspondingly recalculated Fstar0 for semicoherent searches.
Uniform per-detector prior line-vs-Gaussian odds.
BtSG (bool) – If true and transientWindowType is not None, compute the transient \(\ln\mathcal{B}_{\mathrm{tS}/\mathrm{G}}\) statistic from Prix, Giampanis & Messenger (PRD 84, 023007, 2011) (tCWFstatMap marginalised over uniform t0, tau priors). rather than the maxTwoF value.
transientWindowType (str) – If rect or exp, allow for the Fstat to be computed over a transient range. (none instead of None explicitly calls the transient-window function, but with the full range, for debugging.) (If not None, will also force atoms regardless of computeAtoms option.)
t0Band (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave Fstat with t0=minStartTime, tau=maxStartTime-minStartTime.
tauBand (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave Fstat with t0=minStartTime, tau=maxStartTime-minStartTime.
tauMin (int) – Minimum transient duration to cover, defaults to 2*Tsft.
dt0 (int) – Grid resolution in transient start-time, defaults to Tsft.
dtau (int) – Grid resolution in transient duration, defaults to Tsft.
detectors (str) – Two-character references to the detectors for which to use data. Specify None for no constraint. For multiple detectors, separate by commas.
minCoverFreq (float) – The min and max cover frequency passed to lalpulsar.CreateFstatInput. For negative values, these will be used as offsets from the min/max frequency contained in the sftfilepattern. If either is None, the search_ranges argument is used to estimate them. If the automatic estimation fails and you do not have a good idea what to set these two options to, setting both to -0.5 will reproduce the default behaviour of PyFstat <=1.4 and may be a reasonably safe fallback in many cases.
maxCoverFreq (float) – The min and max cover frequency passed to lalpulsar.CreateFstatInput. For negative values, these will be used as offsets from the min/max frequency contained in the sftfilepattern. If either is None, the search_ranges argument is used to estimate them. If the automatic estimation fails and you do not have a good idea what to set these two options to, setting both to -0.5 will reproduce the default behaviour of PyFstat <=1.4 and may be a reasonably safe fallback in many cases.
search_ranges (dict) – Dictionary of ranges in all search parameters, only used to estimate frequency band passed to lalpulsar.CreateFstatInput, if minCoverFreq, maxCoverFreq are not specified (==`None`). For actually running searches, grids/points will have to be passed separately to the .run() method. The entry for each parameter must be a list of length 1, 2 or 3: [single_value], [min,max] or [min,max,step].
injectSources (dict or str) – Either a dictionary of the signal parameters to inject, or a string pointing to a .cff file defining a signal.
injectSqrtSX (float or list or str) – Single-sided PSD values for generating fake Gaussian noise on the fly. Single float or str value: use same for all IFOs. List or comma-separated string: must match len(detectors) and/or the data in sftfilepattern. Detectors will be paired to list elements following alphabetical order.
randSeed (int or None) – random seed for on-the-fly noise generation using injectSqrtSX. Setting this to 0 or None is equivalent; both will randomise the seed, following the behaviour of XLALAddGaussianNoise(), while any number not equal to 0 will produce a reproducible noise realisation.
assumeSqrtSX (float or list or str) – Don’t estimate noise-floors but assume this (stationary) single-sided PSD. Single float or str value: use same for all IFOs. List or comma-separated string: must match len(detectors) and/or the data in sftfilepattern. Detectors will be paired to list elements following alphabetical order. If working with signal-only data, please set assumeSqrtSX=1 .
SSBprec (int) – Flag to set the Solar System Barycentring (SSB) calculation in lalpulsar: 0=Newtonian, 1=relativistic, 2=relativistic optimised, 3=DMoff, 4=NO_SPIN
RngMedWindow (int) – Running-Median window size for F-statistic noise normalization (number of SFT bins).
tCWFstatMapVersion (str) – Choose between implementations of the transient F-statistic funcionality: standard lal implementation, pycuda for GPU version, and some others only for devel/debug.
cudaDeviceName (str) – GPU name to be matched against drv.Device output, only for tCWFstatMapVersion=pycuda.
computeAtoms (bool) – Request calculation of ‘F-statistic atoms’ regardless of transientWindowType.
earth_ephem (str) – Earth ephemeris file path. If None, will check standard sources as per helper_functions.get_ephemeris_files().
sun_ephem (str) – Sun ephemeris file path. If None, will check standard sources as per helper_functions.get_ephemeris_files().
allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength]
- init_computefstatistic()[source]
Initialization step for the F-stastic computation internals.
This sets up the special input and output structures the lalpulsar module needs, the ephemerides, optional on-the-fly signal injections, and extra options for multi-detector consistency checks and transient searches.
All inputs are taken from the pre-initialized object, so this function does not have additional arguments of its own.
- estimate_min_max_CoverFreq()[source]
Extract spanned spin-range at reference -time from the template bank.
To use this method, self.search_ranges must be a dictionary of lists per search parameter which can be either [single_value], [min,max] or [min,max,step].
- get_fullycoherent_detstat(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, tstart=None, tend=None)[source]
Computes the detection statistic(s) fully-coherently at a single point.
Currently supported statistics:
twoF (CW)
log10BSGL (CW or transient)
maxTwoF (transient)
lnBtSG (transient)
All computed statistics are stored as attributes, but only one statistic is returned.
As the basic statistic of this class, twoF is always computed and stored as self.twoF as well, and it is the default return value.
If self.singleFstats, additionally the single-detector 2F-stat values are stored in self.twoFX.
If self.BSGL, the log10BSGL statistic for CWs is additionally stored, and it is returned instead of twoF.
If transient parameters are enabled (self.transientWindowType is set), maxTwoF will always be computed and stored, and returned by default. Depending on the self.BSGL and self.BtSG options, either log10BSGL (a transient version of it, superseding the CW version) or lnBtSG will also be computed, stored, and returned instead of maxTwoF. The full transient-F-stat map is also computed here, but stored in self.FstatMap, not returned.
- Parameters
F0 (float) – Parameters at which to compute the statistic.
F1 (float) – Parameters at which to compute the statistic.
F2 (float) – Parameters at which to compute the statistic.
Alpha (float) – Parameters at which to compute the statistic.
Delta (float) – Parameters at which to compute the statistic.
asini (float, optional) – Optional: Binary parameters at which to compute the statistic.
period (float, optional) – Optional: Binary parameters at which to compute the statistic.
ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.
tp (float, optional) – Optional: Binary parameters at which to compute the statistic.
argp (float, optional) – Optional: Binary parameters at which to compute the statistic.
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. This is only passed on to self.get_transient_detstats(), i.e. only used if self.transientWindowType is set.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. This is only passed on to self.get_transient_detstats(), i.e. only used if self.transientWindowType is set.
- Returns
stat – A single value of the main detection statistic at the input parameter values.
- Return type
float
- get_fullycoherent_twoF(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None)[source]
Computes the fully-coherent 2F statistic at a single point.
NOTE: This always uses the full data set as defined when initialising the search object. If you want to restrict the range of data used for a single 2F computation, you need to set a self.transientWindowType and then call self.get_fullycoherent_detstat() with tstart and tend options instead of this funcion.
- Parameters
F0 (float) – Parameters at which to compute the statistic.
F1 (float) – Parameters at which to compute the statistic.
F2 (float) – Parameters at which to compute the statistic.
Alpha (float) – Parameters at which to compute the statistic.
Delta (float) – Parameters at which to compute the statistic.
asini (float, optional) – Optional: Binary parameters at which to compute the statistic.
period (float, optional) – Optional: Binary parameters at which to compute the statistic.
ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.
tp (float, optional) – Optional: Binary parameters at which to compute the statistic.
argp (float, optional) – Optional: Binary parameters at which to compute the statistic.
- Returns
twoF – A single value of the fully-coherent 2F statistic at the input parameter values. Also stored as self.twoF.
- Return type
float
- get_fullycoherent_single_IFO_twoFs()[source]
Computes single-detector F-stats at a single point.
This requires self.get_fullycoherent_twoF() to be run first.
- Returns
twoFX – A list of the single-detector detection statistics twoF. Also stored as self.twoFX.
- Return type
list
- get_fullycoherent_log10BSGL()[source]
Computes the line-robust statistic log10BSGL at a single point.
This requires self.get_fullycoherent_twoF() and self.get_fullycoherent_single_IFO_twoFs() to be run first.
- Returns
log10BSGL – A single value of the detection statistic log10BSGL at the input parameter values. Also stored as self.log10BSGL.
- Return type
float
- get_transient_detstats(tstart=None, tend=None)[source]
Computes one or more transient detection statistics at a single point.
This requires self.get_fullycoherent_twoF() to be run first.
All computed statistics will be stored as attributes of self, but only one (twoF, log10BSGL or lnBtSG) will be the return value.
The full transient-F-stat map will also be computed here, but stored in self.FstatMap, not returned.
- Parameters
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.
- Returns
detstat – A single value of the main chosen detection statistic (maxTwoF, log10BSGL or lnBtSG) at the input parameter values.
- Return type
float
- get_transient_maxTwoFstat(tstart=None, tend=None)[source]
Computes the transient maxTwoF statistic at a single point.
This requires self.get_fullycoherent_twoF() to be run first, and is itself now only a backwards compatibility / convenience wrapper around the more general get_transient_detstats().
The full transient-F-stat map will also be computed here, but stored in self.FstatMap, not returned.
- Parameters
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.
- Returns
maxTwoF – A single value of the detection statistic at the input parameter values. Also stored as self.maxTwoF.
- Return type
float
- get_transient_log10BSGL()[source]
Computes a transient detection statistic log10BSGL at a single point.
This should normally be called through get_transient_detstats(), but if called stand-alone, it requires self.get_transient_maxTwoFstat() to be run first.
The single-detector 2F-stat values used for that computation (at the index of maxTwoF) are saved in self.twoFXatMaxTwoF, not returned.
- Returns
log10BSGL – A single value of the detection statistic log10BSGL at the input parameter values. Also stored as self.log10BSGL.
- Return type
float
- calculate_twoF_cumulative(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, tstart=None, tend=None, transient_tstart=None, transient_duration=None, num_segments=1000)[source]
Calculate the cumulative twoF over subsets of the observation span.
This means that we consider sub-“segments” of the [tstart,tend] interval, each starting at the overall tstart and with increasing durations, and compute the 2F for each of these, which for a true CW signal should increase roughly with duration towards the full value.
- Parameters
F0 (float) – Parameters at which to compute the cumulative twoF.
F1 (float) – Parameters at which to compute the cumulative twoF.
F2 (float) – Parameters at which to compute the cumulative twoF.
Alpha (float) – Parameters at which to compute the cumulative twoF.
Delta (float) – Parameters at which to compute the cumulative twoF.
asini (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.
period (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.
ecc (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.
tp (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.
argp (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime;. If outside those: auto-truncated.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime;. If outside those: auto-truncated.
num_segments (int) – Number of segments to split [tstart,tend] into.
transient_tstart (float or None) – These are not actually used by this function, but just included so a parameters dict can be safely passed.
transient_duration (float or None) – These are not actually used by this function, but just included so a parameters dict can be safely passed.
- Returns
cumulative_durations (ndarray of shape (num_segments,)) – Offsets of each segment’s tend from the overall tstart.
twoFs (ndarray of shape (num_segments,)) – Values of twoF computed over [[tstart,tstart+duration] for duration in cumulative_durations].
- predict_twoF_cumulative(F0, Alpha, Delta, h0, cosi, psi, tstart=None, tend=None, num_segments=10, **predict_fstat_kwargs)[source]
Calculate expected 2F, with uncertainty, over subsets of the observation span.
This yields the expected behaviour that calculate_twoF_cumulative() can be compared against: 2F for CW signals increases with duration as we take longer and longer subsets of the total observation span.
- Parameters
F0 (float) – Parameters at which to compute the cumulative predicted twoF.
Alpha (float) – Parameters at which to compute the cumulative predicted twoF.
Delta (float) – Parameters at which to compute the cumulative predicted twoF.
h0 (float) – Parameters at which to compute the cumulative predicted twoF.
cosi (float) – Parameters at which to compute the cumulative predicted twoF.
psi (float) – Parameters at which to compute the cumulative predicted twoF.
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.
num_segments (int) – Number of segments to split [tstart,tend] into.
predict_fstat_kwargs – Other kwargs to be passed to helper_functions.predict_fstat().
- Returns
tstart (int) – GPS start time of the observation span.
cumulative_durations (ndarray of shape (num_segments,)) – Offsets of each segment’s tend from the overall tstart.
pfs (ndarray of size (num_segments,)) – Predicted 2F for each segment.
pfs_sigma (ndarray of size (num_segments,)) – Standard deviations of predicted 2F.
- plot_twoF_cumulative(CFS_input, PFS_input=None, tstart=None, tend=None, num_segments_CFS=1000, num_segments_PFS=10, custom_ax_kwargs=None, savefig=False, label=None, outdir=None, **PFS_kwargs)[source]
Plot how 2F accumulates over time.
This compares the accumulation on the actual data set (‘CFS’, from self.calculate_twoF_cumulative()) against (optionally) the average expectation (‘PFS’, from self.predict_twoF_cumulative()).
- Parameters
CFS_input (dict) – Input arguments for self.calculate_twoF_cumulative() (besides [tstart, tend, num_segments]).
PFS_input (dict) – Input arguments for self.predict_twoF_cumulative() (besides [tstart, tend, num_segments]). If None: do not calculate predicted 2F.
tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.
tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.
num_segments_(CFS|PFS) (int) – Number of time segments to (compute|predict) twoF.
custom_ax_kwargs (dict) – Optional axis formatting options.
savefig (bool) – If true, save the figure in outdir. If false, return an axis object without saving to disk.
label (str) – Output filename will be constructed by appending _twoFcumulative.png to this label. (Ignored unless savefig=true.)
outdir (str) – Output folder (ignored unless savefig=true).
PFS_kwargs (dict) – Other kwargs to be passed to self.predict_twoF_cumulative().
- Returns
ax – The axes object containing the plot.
- Return type
matplotlib.axes._subplots_AxesSubplot, optional
- write_atoms_to_file(fnamebase='')[source]
Save F-statistic atoms (time-dependent quantities) for a given parameter-space point.
- Parameters
fnamebase (str) – Basis for output filename, full name will be {fnamebase}_Fstatatoms_{dopplerName}.dat where dopplerName is a canonical lalpulsar formatting of the ‘Doppler’ parameter space point (frequency-evolution parameters).
- class pyfstat.core.SemiCoherentSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.ComputeFstat
A simple semi-coherent search class.
This will split the data set into multiple segments, run a coherent F-stat search over each, and produce a final semi-coherent detection statistic as the sum over segments.
This does not include any concept of refinement between the two steps, as some grid-based semi-coherent search algorithms do; both the per-segment coherent F-statistics and the incoherent sum are done at the same parameter space point.
The implementation is based on a simple trick using the transient F-stat map functionality: basic F-stat atoms are computed only once over the full data set, then the transient code with rectangular ‘windows’ is used to compute the per-segment F-stats, and these are summed to get the semi-coherent result.
Only parameters with a special meaning for SemiCoherentSearch itself are explicitly documented here. For all other parameters inherited from pyfstat.ComputeFStat see the documentation of that class.
- Parameters
label (str) – A label and directory to read/write data from/to.
outdir (str) – A label and directory to read/write data from/to.
tref (int) – GPS seconds of the reference time.
nsegs (int) – The (fixed) number of segments to split the data set into.
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; multiple patterns can be given separated by colons.
minStartTime (int) – Only use SFTs with timestamps starting from this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime]. Also used to set up segment boundaries, i.e. maxStartTime-minStartTime will be divided by nsegs to obtain the per-segment coherence time Tcoh.
maxStartTime (int) – Only use SFTs with timestamps starting from this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime]. Also used to set up segment boundaries, i.e. maxStartTime-minStartTime will be divided by nsegs to obtain the per-segment coherence time Tcoh.
- init_semicoherent_parameters()[source]
Set up a list of equal-length segments and the corresponding transient windows.
For a requested number of segments self.nsegs, self.tboundaries will have self.nsegs+1 entries covering [self.minStartTime,self.maxStartTime] and self.Tcoh will be the total duration divided by self.nsegs.
Each segment is required to be at least two SFTs long.f
- get_semicoherent_det_stat(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, record_segments=False)[source]
Computes the detection statistic (twoF or log10BSGL) semi-coherently at a single point.
As the basic statistic of this class, self.twoF is always computed. If self.singleFstats, additionally the single-detector 2F-stat values are saved in self.twoFX and (optionally) self.twoFX_per_segment.
- Parameters
F0 (float) – Parameters at which to compute the statistic.
F1 (float) – Parameters at which to compute the statistic.
F2 (float) – Parameters at which to compute the statistic.
Alpha (float) – Parameters at which to compute the statistic.
Delta (float) – Parameters at which to compute the statistic.
asini (float, optional) – Optional: Binary parameters at which to compute the statistic.
period (float, optional) – Optional: Binary parameters at which to compute the statistic.
ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.
tp (float, optional) – Optional: Binary parameters at which to compute the statistic.
argp (float, optional) – Optional: Binary parameters at which to compute the statistic.
record_segments (boolean) – If True, store the per-segment F-stat values as self.twoF_per_segment and (if self.singleFstats) the per-detector per-segment F-stats as self.twoFX_per_segment.
- Returns
stat – A single value of the detection statistic (semi-coherent twoF or log10BSGL) at the input parameter values. Also stored as self.twoF or self.log10BSGL.
- Return type
float
- get_semicoherent_twoF(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, record_segments=False)[source]
Computes the semi-coherent twoF statistic at a single point.
- Parameters
F0 (float) – Parameters at which to compute the statistic.
F1 (float) – Parameters at which to compute the statistic.
F2 (float) – Parameters at which to compute the statistic.
Alpha (float) – Parameters at which to compute the statistic.
Delta (float) – Parameters at which to compute the statistic.
asini (float, optional) – Optional: Binary parameters at which to compute the statistic.
period (float, optional) – Optional: Binary parameters at which to compute the statistic.
ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.
tp (float, optional) – Optional: Binary parameters at which to compute the statistic.
argp (float, optional) – Optional: Binary parameters at which to compute the statistic.
record_segments (boolean) – If True, store the per-segment F-stat values as self.twoF_per_segment.
- Returns
twoF – A single value of the semi-coherent twoF statistic at the input parameter values. Also stored as self.twoF.
- Return type
float
- get_semicoherent_single_IFO_twoFs(record_segments=False)[source]
Computes the semi-coherent single-detector F-statss at a single point.
This requires self.get_semicoherent_twoF() to be run first.
- Parameters
record_segments (boolean) – If True, store the per-detector per-segment F-stat values as self.twoFX_per_segment.
- Returns
twoFX – A list of the single-detector detection statistics twoF. Also stored as self.twoFX.
- Return type
list
- get_semicoherent_log10BSGL()[source]
Computes the semi-coherent log10BSGL statistic at a single point.
This requires self.get_semicoherent_twoF() and self.get_semicoherent_single_IFO_twoFs() to be run first.
- Returns
log10BSGL – A single value of the semi-coherent log10BSGL statistic at the input parameter values. Also stored as self.log10BSGL.
- Return type
float
- class pyfstat.core.SearchForSignalWithJumps(*args, **kwargs)[source]
Bases:
pyfstat.core.BaseSearchClass
Internal helper class with some useful methods for glitches or timing noise.
Users should never need to interact with this class, just with the derived search classes.
- class pyfstat.core.SemiCoherentGlitchSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.SearchForSignalWithJumps
,pyfstat.core.ComputeFstat
A semi-coherent search for CW signals from sources with timing glitches.
This implements a basic semi-coherent F-stat search in which the data is divided into segments either side of the proposed glitch epochs and the fully-coherent F-stat in each segment is summed to give the semi-coherent F-stat.
Only parameters with a special meaning for SemiCoherentGlitchSearch itself are explicitly documented here. For all other parameters inherited from pyfstat.ComputeFStat see the documentation of that class.
- Parameters
label (str) – A label and directory to read/write data from/to.
outdir (str) – A label and directory to read/write data from/to.
tref (int) – GPS seconds of the reference time, and start and end of the data.
minStartTime (int) – GPS seconds of the reference time, and start and end of the data.
maxStartTime (int) – GPS seconds of the reference time, and start and end of the data.
nglitch (int) – The (fixed) number of glitches. This is also allowed to be zero, but occasionally this causes issues, in which case please use the basic ComputeFstat class instead.
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; multiple patterns can be given separated by colons.
theta0_idx (int) – Index (zero-based) of which segment the theta (searched parameters) refer to. This is useful if providing a tight prior on theta to allow the signal to jump to theta (and not just from).
- get_semicoherent_nglitch_twoF(F0, F1, F2, Alpha, Delta, *args)[source]
Returns the semi-coherent glitch summed twoF.
- Parameters
F0 (float) – Parameters at which to compute the statistic.
F1 (float) – Parameters at which to compute the statistic.
F2 (float) – Parameters at which to compute the statistic.
Alpha (float) – Parameters at which to compute the statistic.
Delta (float) – Parameters at which to compute the statistic.
args (dict) – Additional arguments for the glitch parameters; see the source code for full details.
- Returns
twoFSum – A single value of the semi-coherent summed detection statistic at the input parameter values.
- Return type
float
pyfstat.grid_based_searches module
PyFstat search classes using grid-based methods.
- class pyfstat.grid_based_searches.GridSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.BaseSearchClass
A search evaluating the F-statistic over a regular grid in parameter space.
This implements a simple ‘square box’ grid with fixed spacing and ranges in each dimension, i.e. for each parameter there’s a simple 1D list of grid points and the total grid is just the Cartesian product of these.
For N parameter space dimensions and a total of M points in the product grid, the basic output is a (N+1,M)-dimensional array with the detection statistic (twoF or log10BSGL) appended.
NOTE: if a large number of grid points are used, checks against cached data may be slow as the array is loaded into memory. To avoid this, run with the clean option which uses a generator instead.
Most parameters are the same as for the core.ComputeFstat class, only the additional ones are documented here:
- Parameters
label (str) – Output filenames will be constructed using this label.
outdir (str) – Output directory.
F0s (tuple) – A length 3 tuple describing the grid for each parameter, e.g [F0min, F0max, dF0]. Alternatively, for a fixed value simply give [F0]. Unless input_arrays=True, then these are the exact arrays to search over.
F1s (tuple) – A length 3 tuple describing the grid for each parameter, e.g [F0min, F0max, dF0]. Alternatively, for a fixed value simply give [F0]. Unless input_arrays=True, then these are the exact arrays to search over.
F2s (tuple) – A length 3 tuple describing the grid for each parameter, e.g [F0min, F0max, dF0]. Alternatively, for a fixed value simply give [F0]. Unless input_arrays=True, then these are the exact arrays to search over.
Alphas (tuple) – A length 3 tuple describing the grid for each parameter, e.g [F0min, F0max, dF0]. Alternatively, for a fixed value simply give [F0]. Unless input_arrays=True, then these are the exact arrays to search over.
Deltas (tuple) – A length 3 tuple describing the grid for each parameter, e.g [F0min, F0max, dF0]. Alternatively, for a fixed value simply give [F0]. Unless input_arrays=True, then these are the exact arrays to search over.
nsegs (int) – Number of segments to split the data set into. If nsegs=1, the basic ComputeFstat class is used. If nsegs>1, the SemiCoherentSearch class is used.
input_arrays (bool) – If true, use the F0s, F1s, etc as arrays just as they are given (do not interpret as 3-tuples of [min,max,step]).
- tex_labels = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'lnBtSG': '$\\ln\\mathcal{B}_{\\mathrm{tS/G}}$', 'log10BSGL': '$\\log_{10}\\mathcal{B}_{\\mathrm{S/GL}}$', 'maxTwoF': '$\\max\\widetilde{2\\mathcal{F}}$', 'twoF': '$\\widetilde{2\\mathcal{F}}$'}
Formatted labels used for plot annotations.
- tex_labels0 = {'Alpha': '$-\\alpha_0$', 'Delta': '$-\\delta_0$', 'F0': '$-f_0$', 'F1': '$-\\dot{f}_0$', 'F2': '$-\\ddot{f}_0$'}
Formatted labels used for annotating central values in plots.
- fmt_detstat = '%.9g'
Standard output precision for detection statistics.
- check_old_data_is_okay_to_use()[source]
Check if an existing output file matches this search and reuse the results.
Results will be loaded from old output file, and no new search run, if all of the following checks pass:
Output file with matching name found in outdir.
Output file is not older than SFT files matching sftfilepattern.
Parameters string in file header matches current search setup.
Data in old file can be loaded successfully, its input parts (i.e. minus the detection statistic columns) matches in dimension with current grid, and the values in those input columns match with the current grid.
Through helper_functions.read_txt_file_with_header(), the existing file is read in with np.genfromtxt().
- run(return_data=False)[source]
Execute the actual search over the full grid.
This iterates over all points in the multi-dimensional product grid and the end result is either returned as a numpy array or saved to disk.
- Parameters
return_data (boolean) – If true, the final inputs+outputs data set is returned as a numpy array. If false, it is saved to disk and nothing is returned.
- Returns
data – The final inputs+outputs data set. Only if return_data=true.
- Return type
np.ndarray
- save_array_to_disk()[source]
Save the results array to a txt file.
This includes a header with version and parameters information.
It should be flexible enough to be reused by child classes, as long as the _get_savetxt_fmt_dict() method is suitably overridden to account for any additional parameters.
- plot_1D(xkey, ax=None, x0=None, xrescale=1, savefig=True, xlabel=None, ylabel=None, agg_chunksize=None)[source]
Make a plot of the detection statistic over a single grid dimension.
- Parameters
xkey (str) – The name of the search parameter to plot against.
ax (matplotlib.axes._subplots_AxesSubplot or None) – An optional pre-existing axes set to draw into.
x0 (float or None) – Plot x values relative to this central value.
xrescale (float) – Rescale all x values by this factor.
savefig (bool) – If true, save the figure in self.outdir. If false, return an axis object without saving to disk.
xlabel (str or None) – Override default text label for the x-axis.
ylabel (str or None) – Override default text label for the y-axis.
agg_chunksize (int or None) – Set this to some high value to work around matplotlib ‘Exceeded cell block limit’ errors.
- Returns
ax – The axes object containing the plot, only if savefig=false.
- Return type
matplotlib.axes._subplots_AxesSubplot, optional
- plot_2D(xkey, ykey, ax=None, savefig=True, vmin=None, vmax=None, add_mismatch=None, xN=None, yN=None, flat_keys=[], rel_flat_idxs=[], flatten_method=<function amax>, title=None, predicted_twoF=None, cm=None, cbarkwargs={}, x0=None, y0=None, colorbar=False, xrescale=1, yrescale=1, xlabel=None, ylabel=None, zlabel=None)[source]
Plots the detection statistic over a 2D grid.
FIXME: this will currently fail if the search went over >2 dimensions.
- Parameters
xkey (str) – The name of the first search parameter to plot against.
ykey (str) – The name of the second search parameter to plot against.
ax (matplotlib.axes._subplots_AxesSubplot or None) – An optional pre-existing axes set to draw into.
savefig (bool) – If true, save the figure in self.outdir. If false, return an axis object without saving to disk.
vmin (float or None) – Cutoffs for rescaling the colormap.
vmax (float or None) – Cutoffs for rescaling the colormap.
add_mismatch (tuple or None) – If given a tuple (xhat, yhat, Tseg), add a secondary axis with the metric mismatch from the point (xhat, yhat) with duration Tseg.
xN (int or None) – Number of tick label intervals.
yN (int or None) – Number of tick label intervals.
flat_keys (list) – Keys to be used in flattening higher-dimensional arrays.
rel_flat_idxs (list) – Indices to be used in flattening higher-dimensional arrays.
flatten_method (numpy function) – Function to use in flattening the flat_keys, default: np.max.
title (str or None) – Optional plot title text.
predicted_twoF (float or None) – Expected/predicted value of twoF, used to rescale the z-axis.
cm (matplotlib.colors.ListedColormap or None) – Override standard (viridis) colormap.
cbarkwargs (dict) – Additional arguments for colorbar formatting.
x0 (float) – Plot x values relative to this central value.
y0 (float) – Plot y values relative to this central value.
xrescale (float) – Rescale all x values by this factor.
yrescale (float) – Rescale all y values by this factor.
xlabel (str) – Override default text label for the x-axis.
ylabel (str) – Override default text label for the y-axis.
zlabel (str) – Override default text label for the z-axis.
- Returns
ax – The axes object containing the plot, only if savefig=false.
- Return type
matplotlib.axes._subplots_AxesSubplot, optional
- get_max_det_stat()[source]
Get the maximum detection statistic over the grid.
This requires the run() method to have been called before.
- Returns
d – Dictionary containing parameters and detection statistic at the maximum.
- Return type
dict
- get_max_twoF()[source]
Get the maximum twoF over the grid.
This requires the run() method to have been called before.
- Returns
d – Dictionary containing parameters and twoF value at the maximum.
- Return type
dict
- print_max_twoF()[source]
Get and print the maximum twoF point over the grid.
This prints out the full dictionary from get_max_twoF(), i.e. the maximum value and its corresponding parameters.
- set_out_file(extra_label=None)[source]
Set (or reset) the name of the main output file.
File will always be stored in self.outdir and the base of the name be determined from self.label and other parts of the search setup, but this method allows to attach an extra_label bit if desired.
- Parameters
extra_label (str) – Additional text bit to be attached at the end of the filename (but before the extension).
- class pyfstat.grid_based_searches.TransientGridSearch(*args, **kwargs)[source]
Bases:
pyfstat.grid_based_searches.GridSearch
A search for transient CW-like signals using the F-statistic.
This is based on the transient signal model and transient-F-stat algorithm from Prix, Giampanis & Messenger (PRD 84, 023007, 2011): https://arxiv.org/abs/1104.1704
The frequency evolution parameters are searched over in a grid just like in the normal GridSearch, then at each point the time-dependent ‘atoms’ are used to evaluate partial sums of the F-statistic over a 2D array in transient start times t0 and duration parameters tau.
The signal templates are modulated by a ‘transient window function’ which can be
none (standard, persistent CW signal)
rect (rectangular: constant amplitude within [t0,t0+tau], zero outside)
exp (exponential decay over [t0,t0+3*tau], zero outside)
This class currently only supports fully-coherent searches (nsegs=1 is hardcoded).
Also see Keitel & Ashton (CQG 35, 205003, 2018): https://arxiv.org/abs/1805.05652 for a detailed discussion of the GPU implementation.
NOTE for GPU users (tCWFstatMapVersion=”pycuda”): The underlying ComputeFstat class tries to conveniently deal with GPU context management behind the scenes. A known problematic case is if you try to instantiate it twice from the same session/script. If you then get some messages like RuntimeError: make_default_context() and invalid device context, that is because the GPU is still blocked from the first instance when you try to initiate the second. To avoid this problem, use context management:
with pyfstat.TransientGridSearch( [...], tCWFstatMapVersion="pycuda", ) as search: search.search.run()
or manually call the search.search.finalizer_() method where needed.
Most parameters are the same as for GridSearch and the core.ComputeFstat class, only the additional ones are documented here:
- Parameters
transientWindowType (str) – If rect or exp, allow for the Fstat to be computed over a transient range. (none instead of None explicitly calls the transient-window function, but with the full range, for debugging.)
t0Band (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave F-stat with t0=minStartTime, tau=maxStartTime-minStartTime.
tauBand (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave F-stat with t0=minStartTime, tau=maxStartTime-minStartTime.
tauMin (int) – Minimum transient duration to cover, defaults to 2*Tsft.
dt0 (int) – Grid resolution in transient start-time, defaults to Tsft.
dtau (int) – Grid resolution in transient duration, defaults to Tsft.
outputTransientFstatMap (bool) – If true, write additional output files for (t0,tau) F-stat maps. (One file for each grid point!)
outputAtoms (bool) – If true, write additional output files for the F-stat atoms. (One file for each grid point!)
tCWFstatMapVersion (str) – Choose between implementations of the transient F-statistic funcionality: standard lal implementation, pycuda for GPU version, and some others only for devel/debug.
cudaDeviceName (str) – GPU name to be matched against drv.Device output, only for tCWFstatMapVersion=pycuda.
- run(return_data=False)[source]
Execute the actual search over the full grid.
This iterates over all points in the multi-dimensional product grid and the end result is either returned as a numpy array or saved to disk.
If the outputTransientFstatMap or outputAtoms options have been set when initiating the search, additional files are written for each frequency-evolution parameter-space point (‘Doppler’ point).
- Parameters
return_data (boolean) – If true, the final inputs+outputs data set is returned as a numpy array. If false, it is saved to disk and nothing is returned.
- Returns
data – The final inputs+outputs data set. Only if return_data=true.
- Return type
np.ndarray
- get_transient_fstat_map_filename(param_point)[source]
Filename convention for given grid point: freq_alpha_delta_f1dot_f2dot
- Parameters
param_point (tuple, dict, list, np.void or np.ndarray) – A multi-dimensional parameter point. If not a type with named fields (e.g. a plain tuple or list), the order must match that of self.output_keys.
- Returns
f – The constructed filename.
- Return type
str
- class pyfstat.grid_based_searches.SliceGridSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
- class pyfstat.grid_based_searches.GridUniformPriorSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
- class pyfstat.grid_based_searches.GridGlitchSearch(*args, **kwargs)[source]
Bases:
pyfstat.grid_based_searches.GridSearch
A grid search using the SemiCoherentGlitchSearch class.
This implements a basic semi-coherent F-stat search in which the data is divided into segments either side of the proposed glitch epochs and the fully-coherent F-stat in each segment is summed to give the semi-coherent F-stat.
This class currently only works for a single glitch in the observing time.
Most parameters are the same as for GridSearch and the core.SemiCoherentGlitchSearch class, only the additional ones are documented here:
- Parameters
delta_F0s (tuple) – A length 3 tuple describing the grid of frequency jumps, or just [delta_F0] for a fixed value.
delta_F1s (tuple) – A length 3 tuple describing the grid of spindown parameter jumps, or just [delta_F1] for a fixed value.
tglitchs (tuple) – A length 3 tuple describing the grid of glitch epochs, or just [tglitch] for a fixed value. These are relative time offsets, referenced to zero at minStartTime.
- class pyfstat.grid_based_searches.SlidingWindow(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
- class pyfstat.grid_based_searches.FrequencySlidingWindow(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
- class pyfstat.grid_based_searches.EarthTest(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
- class pyfstat.grid_based_searches.DMoff_NO_SPIN(*args, **kwargs)[source]
Bases:
pyfstat.core.DefunctClass
- last_supported_version = '1.9.0'
pyfstat.gridcorner module
A corner plotting tool for an array (grid) of dependent values.
Given an N-dimensional set of data (i.e. some function evaluated over a grid of coordinates), plot all possible 1D and 2D projections in the style of a ‘corner’ plot.
This code has been copied from Gregory Ashton’s repository https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner and it uses both the central idea and some specific code from Daniel Foreman-Mackey’s https://github.com/dfm/corner.py re-used under the following licence requirements:
Copyright (c) 2013-2020 Daniel Foreman-Mackey
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
The views and conclusions contained in the software and documentation are those of the authors and should not be interpreted as representing official policies, either expressed or implied, of the FreeBSD Project.
- pyfstat.gridcorner.log_mean(loga, axis)[source]
Calculate the log(<a>) mean
Given N logged value log, calculate the log_mean log(<loga>)=log(sum(np.exp(loga))) - log(N). Useful for marginalizing over logged likelihoods for example.
- Parameters
loga (array_like) – Input_array.
axies (None or int or type of ints, optional) – Axis or axes over which the sum is taken. By default axis is None, and all elements are summed.
- Returns
log_mean – The logged average value (shape loga.shape)
- Return type
ndarry
- pyfstat.gridcorner.idx_array_slice(D, axis, slice_idx)[source]
Return the slice along the given axis
- pyfstat.gridcorner.gridcorner(D, xyz, labels=None, projection='max_slice', max_n_ticks=4, factor=2, whspace=0.05, showDvals=True, lines=None, label_offset=0.4, **kwargs)[source]
Generate a grid corner plot
- Parameters
D (array_like) – N-dimensional data to plot, D.shape should be (n1, n2,…, nN), where N, is the number of grid points along dimension i.
xyz (list) – List of 1-dimensional arrays of coordinates. xyz[i] should have length N (see help for D).
labels (list) – N+1 length list of labels; the first N correspond to the coordinates labels, the final label is for the dependent (D) variable.
projection (str or func) – If a string, one of {“log_mean”, “max_slice”} to use inbuilt functions to calculate either the logged mean or maximum slice projection. Else a function to use for projection, must take an `axis argument. Default is gridcorner.max_slice(), to project out a slice along the maximum.
max_n_ticks (int) – Number of ticks for x and y axis of the pcolormesh plots.
factor (float) – Controls the size of one window.
showDvals (bool) – If true (default) show the D values on the right-hand-side of the 1D plots and add a label.
lines (array_like) – N-dimensional list of values to delineate.
- Returns
The figure and NxN set of axes
- Return type
fig, axes
pyfstat.helper_functions module
A collection of helpful functions to facilitate ease-of-use of PyFstat.
Most of these are used internally by other parts of the package and are of interest mostly only for developers, but others can also be helpful for end users.
- pyfstat.helper_functions.set_up_optional_tqdm()[source]
Provides local replacement for tqdm if it cannot be imported.
- pyfstat.helper_functions.set_up_matplotlib_defaults()[source]
Sets some defaults for matplotlib plotting.
- pyfstat.helper_functions.set_up_command_line_arguments()[source]
Parse global commandline arguments.
- pyfstat.helper_functions.get_ephemeris_files()[source]
Set the ephemeris files to use for the Earth and Sun.
This looks first for a configuration file ~/.pyfstat.conf giving individual earth/sun file paths like this:
` earth_ephem = '/my/path/earth00-40-DE405.dat.gz' sun_ephem = '/my/path/sun00-40-DE405.dat.gz' `
If such a file is not found or does not conform to that format, then we rely on lal’s recently improved ability to find proper default fallback paths for the [earth/sun]00-40-DE405 ephemerides with both pip- and conda-installed packages,
Alternatively, ephemeris options can be set manually on each class instantiation.
NOTE that the $LALPULSAR_DATADIR environment variable is no longer supported!
- Returns
earth_ephem, sun_ephem – Paths of the two files containing positions of Earth and Sun.
- Return type
str
- pyfstat.helper_functions.round_to_n(x, n)[source]
Simple rounding function for getting a fixed number of digits.
- Parameters
x (float) – The number to round.
n (int) – The number of digits to round to (before plus after the decimal separator).
- Returns
rounded – The rounded number.
- Return type
float
- pyfstat.helper_functions.texify_float(x, d=2)[source]
Format float numbers nicely for LaTeX output, including rounding.
Numbers with absolute values between 0.01 and 100 will be returned in plain float format, while smaller or larger numbers will be returned in powers-of-ten notation.
- Parameters
x (float) – The number to round and format.
n (int) – The number of digits to round to (before plus after the decimal separator).
- Returns
formatted – The formatted string.
- Return type
str
- pyfstat.helper_functions.initializer(func)[source]
Decorator to automatically assign the parameters of a class instantiation to self.
- pyfstat.helper_functions.get_peak_values(frequencies, twoF, threshold_2F, F0=None, F0range=None)[source]
Find a set of local peaks of twoF values over a 1D frequency list.
- Parameters
frequencies (np.ndarray) – 1D array of frequency values.
twoF (np.ndarray) – Corresponding 1D array of detection statistic values.
threshold_2F (float) – Only report peaks above this threshold.
F0 (float or None) – Only report peaks within [F0-F0range,F0+F0range].
F0range (float or None) – Only report peaks within [F0-F0range,F0+F0range].
- Returns
F0maxs (np.ndarray) – 1D array of peak frequencies.
twoFmaxs (np.ndarray) – 1D array of peak twoF values.
freq_errs (np.ndarray) – 1D array of peak frequency estimation error, taken as the spacing between neighbouring input frequencies.
- pyfstat.helper_functions.get_comb_values(F0, frequencies, twoF, period, N=4)[source]
Check which points of a [frequencies,twoF] set correspond to a certain comb of frequencies.
- Parameters
F0 (float) – Base frequency for the comb.
frequencies (np.ndarray) – 1D array of frequency values.
twoF (np.ndarray) – Corresponding 1D array of detection statistic values.
period (str) – Modulation type of the comb, either sidereal or terrestrial.
N (int) – Number of comb harmonics.
- Returns
comb_frequencies (np.ndarray) – 1D array of relative frequency offsets for the comb harmonics.
comb_twoFs (np.ndarray) – 1D array of twoF values at the closest-matching frequencies.
freq_errs (np.ndarray) – 1D array of frequency match error, taken as the spacing between neighbouring input frequencies.
- pyfstat.helper_functions.run_commandline(cl, log_level=20, raise_error=True, return_output=True)[source]
Run a string cmd as a subprocess, check for errors and return output.
- Parameters
cl (str) – Command to run
log_level (int) – Sets the logging level for some of this function’s messages. See https://docs.python.org/library/logging.html#logging-levels Default is ‘20’ (INFO). FIXME: Not used for all messages.
raise_error (bool) – If True, raise an error if the subprocess fails. If False, continue and just return 0.
return_output (bool) – If True, return the captured output of the subprocess (stdout and stderr). If False, return nothing on successful execution.
- Returns
out – The captured output of the subprocess (stdout and stderr) if return_output=True. 0 on failed execution if raise_error=False.
- Return type
str or int, optional
- pyfstat.helper_functions.convert_array_to_gsl_matrix(array)[source]
Convert a numpy array to a LAL-wrapped GSL matrix.
- Parameters
array (np.ndarray) – The array to convert. array.shape must have 2 dimensions.
- Returns
gsl_matrix – The LAL-wrapped GSL matrix object.
- Return type
lal.gsl_matrix
- pyfstat.helper_functions.get_sft_as_arrays(sftfilepattern, fMin=None, fMax=None, constraints=None)[source]
- Parameters
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; multiple patterns can be given separated by colons.
fMin (float or None) – Restrict frequency range to [fMin, fMax]. If None, retreive the full frequency range.
fMax (float or None) – Restrict frequency range to [fMin, fMax]. If None, retreive the full frequency range.
constraints (lalpulsar.SFTConstraints() or None) – Constrains to be fed into XLALSFTdataFind to specify detector, GPS time range or timestamps to be retrieved.
- Returns
freqs (np.ndarray) – The frequency bins in each SFT. These will be the same for each SFT, so only a single 1D array is returned.
times (dict of np.ndarray) – The SFT start times as a dictionary of 1D arrays, one for each detector. Keys correspond to the official detector names as returned by XLALlalpulsar.ListIFOsInCatalog.
data (dict of np.ndarray) – A dictionary of 2D arrays of the complex Fourier amplitudes of the SFT data for each detector in each frequency bin at each timestamp. Keys correspond to the official detector names as returned by XLALlalpulsar.ListIFOsInCatalog.
- pyfstat.helper_functions.get_sft_array(sftfilepattern, F0=None, dF0=None)[source]
Return the raw data (absolute values) from a set of SFTs.
FIXME: currently only returns data for first detector.
- Parameters
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; multiple patterns can be given separated by colons.
F0 (float or None) – Restrict frequency range to [F0-dF0,F0+dF0].
dF0 (float or None) – Restrict frequency range to [F0-dF0,F0+dF0].
- Returns
times (np.ndarray) – The SFT starttimes as a 1D array.
freqs (np.ndarray) – The frequency bins in each SFT. These will be the same for each SFT, so only a single 1D array is returned.
data (np.ndarray) – A 2D array of the absolute values of the SFT data in each frequency bin at each timestamp.
- pyfstat.helper_functions.get_covering_band(tref, tstart, tend, F0, F1, F2, F0band=0.0, F1band=0.0, F2band=0.0, maxOrbitAsini=0.0, minOrbitPeriod=0.0, maxOrbitEcc=0.0)[source]
Get the covering band for CW signals for given time and parameter ranges.
This uses the lalpulsar function XLALCWSignalCoveringBand(), accounting for the spin evolution of the signals within the given [F0,F1,F2] ranges, the maximum possible Dopper modulation due to detector motion (i.e. for the worst-case sky locations), and for worst-case binary orbital motion.
- Parameters
tref (int) – Reference time (in GPS seconds) for the signal parameters.
tstart (int) – Start time (in GPS seconds) for the signal evolution to consider.
tend (int) – End time (in GPS seconds) for the signal evolution to consider.
F0 (float) – Minimum frequency and spin-down of signals to be covered.
F1 (float) – Minimum frequency and spin-down of signals to be covered.
F1 – Minimum frequency and spin-down of signals to be covered.
F0band (float) – Ranges of frequency and spin-down of signals to be covered.
F1band (float) – Ranges of frequency and spin-down of signals to be covered.
F1band – Ranges of frequency and spin-down of signals to be covered.
maxOrbitAsini (float) – Largest orbital projected semi-major axis to be covered.
minOrbitPeriod (float) – Shortest orbital period to be covered.
maxOrbitEcc (float) – Highest orbital eccentricity to be covered.
- Returns
minCoverFreq, maxCoverFreq – Estimates of the minimum and maximum frequencies of the signals from the given parameter ranges over the [tstart,tend] duration.
- Return type
float
- pyfstat.helper_functions.match_commandlines(cl1, cl2, be_strict_about_full_executable_path=False)[source]
Check if two commandline strings match element-by-element, regardless of order.
- Parameters
cl1 (str) – Commandline strings of executable –key1=val1 –key2=val2 etc format.
cl2 (str) – Commandline strings of executable –key1=val1 –key2=val2 etc format.
be_strict_about_full_executable_path (bool) – If False (default), only checks the basename of the executable. If True, requires its full path to match.
- Returns
match – Whether the executable and all key=val pairs of the two strings matched.
- Return type
bool
- pyfstat.helper_functions.get_version_string()[source]
Get the canonical version string of the currently running PyFstat instance.
- Returns
version – The version string.
- Return type
str
- pyfstat.helper_functions.get_doppler_params_output_format(keys)[source]
Set a canonical output precision for frequency evolution parameters.
This uses the same format (%.16g) as the write_FstatCandidate_to_fp() function of lalapps_ComputeFstatistic_v2.
This assigns that format to each parameter name in keys which matches a hardcoded list of known standard ‘Doppler’ parameters, and ignores any others.
- Parameters
keys (dict) – The parameter keys for which to select formats.
- Returns
fmt – A dictionary assigning the default format to each parameter key from the hardcoded list of standard ‘Doppler’ parameters.
- Return type
dict
- pyfstat.helper_functions.read_txt_file_with_header(f, names=True, comments='#')[source]
Wrapper to np.genfromtxt with smarter handling of variable-length commented headers.
The header is identified as an uninterrupted block of lines from the beginning of the file, each starting with the given comments character.
After identifying a header of length Nhead, this function then tells np.genfromtxt() to skip Nhead-1 lines (to allow for reading field names from the last commented line before the actual data starts).
- Parameters
f (str) – Name of the file to read.
names (bool) – Passed on to np.genfromtxt(): If True, the field names are read from the last header line.
comments (str) – The character used to indicate the start of a comment. Also passed on to np.genfromtxt().
- Returns
data – The data array read from the file after skipping the header.
- Return type
np.ndarray
- pyfstat.helper_functions.get_lalapps_commandline_from_SFTDescriptor(descriptor)[source]
Extract a lalapps commandline from the ‘comment’ entry of a SFT descriptor.
Most SFT creation tools save their commandline into that entry, so we can extract it and reuse it to reproduce that data.
- Parameters
descriptor (SFTDescriptor) – Element of a lalpulsar.SFTCatalog structure.
- Returns
cmd – A lalapps commandline string, or an empty string if ‘lalapps’ not found in comment.
- Return type
str
- pyfstat.helper_functions.read_parameters_dict_lines_from_file_header(outfile, comments='#', strip_spaces=True)[source]
Load a list of pretty-printed parameters dictionary lines from a commented file header.
Returns a list of lines from a commented file header that match the pretty-printed parameters dictionary format as generated by BaseSearchClass.get_output_file_header(). The opening/closing bracket lines ({,`}`) are not included. Newline characters at the end of each line are stripped.
- Parameters
outfile (str) – Name of a PyFstat-produced output file.
comments (str) – Comment character used to start header lines.
strip_spaces (bool) – Whether to strip leading/trailing spaces.
- Returns
dict_lines – A list of unparsed pprinted dictionary entries.
- Return type
list
- pyfstat.helper_functions.get_parameters_dict_from_file_header(outfile, comments='#', eval_values=False)[source]
Load a parameters dict from a commented file header.
Returns a parameters dictionary, as generated by BaseSearchClass.get_output_file_header(), from an output file header. Always returns a proper python dictionary, but the values will be unparsed strings if not requested otherwise.
- Parameters
outfile (str) – Name of a PyFstat-produced output file.
comments (str) – Comment character used to start header lines.
eval_values (bool) – If False, return dictionary values as unparsed strings. If True, evaluate each of them. DANGER! Only do this if you trust the source of the file!
- Returns
params_dict – A dictionary of parameters (with values either as unparsed strings, or evaluated).
- Return type
dictionary
- pyfstat.helper_functions.read_par(filename=None, label=None, outdir=None, suffix='par', comments=['%', '#'], raise_error=False)[source]
Read in a .par or .loudest file, returns a dictionary of the key=val pairs.
Notes
This can also be used to read in .loudest files produced by lalapps_ComputeFstatistic_v2, or any file which has rows of key=val data (in which the val can be understood using eval(val)).
- Parameters
filename (str) – Filename (path) containing rows of key=val data to read in.
label (str, optional) – If filename is None, form the file to read as outdir/label.suffix.
outdir (str, optional) – If filename is None, form the file to read as outdir/label.suffix.
suffix (str, optional) – If filename is None, form the file to read as outdir/label.suffix.
comments (str or list of strings, optional) – Characters denoting that a row is a comment.
raise_error (bool, optional) – If True, raise an error for lines which are not comments, but cannot be read.
- Returns
d – The key=val pairs as a dictionary.
- Return type
dict
- pyfstat.helper_functions.get_dictionary_from_lines(lines, comments, raise_error)[source]
Return a dictionary of key=val pairs for each line in a list.
- Parameters
lines (list of strings) – The list of lines to parse.
comments (str or list of strings) – Characters denoting that a row is a comment.
raise_error (bool) – If True, raise an error for lines which are not comments, but cannot be read. Note that CFSv2 “loudest” files contain complex numbers which fill raise an error unless this is set to False.
- Returns
d – The key=val pairs as a dictionary.
- Return type
dict
- pyfstat.helper_functions.get_predict_fstat_parameters_from_dict(signal_parameters, transientWindowType=None)[source]
Extract a subset of parameters as needed for predicting F-stats.
Given a dictionary with arbitrary signal parameters, this extracts only those ones required by helper_functions.predict_fstat(): Freq, Alpha, Delta, h0, cosi, psi. Also preserves transient parameters, if included in the input dict.
- Parameters
signal_parameters (dict) – Dictionary containing at least those signal parameters required by helper_functions.predict_fstat. This dictionary’s keys must follow the PyFstat convention (e.g. F0 instead of Freq).
transientWindowType (str) – Transient window type to store in the output dict. Currently required because the typical input dicts produced by various PyFstat functions tend not to store this property. If there is a key with this name already, its value will be overwritten.
- Returns
predict_fstat_params – The dictionary of selected parameters.
- Return type
dict
- pyfstat.helper_functions.predict_fstat(h0=None, cosi=None, psi=None, Alpha=None, Delta=None, F0=None, sftfilepattern=None, timestampsFiles=None, minStartTime=None, duration=None, IFOs=None, assumeSqrtSX=None, tempory_filename='fs.tmp', earth_ephem=None, sun_ephem=None, transientWindowType='none', transientStartTime=None, transientTau=None)[source]
Wrapper to lalapps_PredictFstat for predicting expected F-stat values.
- Parameters
h0 (float) – Signal parameters, see lalapps_PredictFstat –help for more info.
cosi (float) – Signal parameters, see lalapps_PredictFstat –help for more info.
psi (float) – Signal parameters, see lalapps_PredictFstat –help for more info.
Alpha (float) – Signal parameters, see lalapps_PredictFstat –help for more info.
Delta (float) – Signal parameters, see lalapps_PredictFstat –help for more info.
F0 (float or None) – Signal frequency. Only needed for noise floor estimation when given sftfilepattern but assumeSqrtSX=None. The actual F-stat prediction is frequency-independent.
sftfilepattern (str or None) – Pattern matching the SFT files to use for inferring detectors, timestamps and/or estimating the noise floor.
timestampsFiles (str or None) – Comma-separated list of per-detector files containing timestamps to use. Only used if sftfilepattern=None.
minStartTime (int or None) – If sftfilepattern given: used as optional constraints. If timestampsFiles given: ignored. If neither given: used as the interval for prediction.
duration (int or None) – If sftfilepattern given: used as optional constraints. If timestampsFiles given: ignored. If neither given: used as the interval for prediction.
IFOs (str or None) – Comma-separated list of detectors. Required if sftfilepattern=None, ignored otherwise.
assumeSqrtSX (float or str) – Assume stationary per-detector noise-floor instead of estimating from SFTs. Single float or str value: use same for all IFOs. Comma-separated string: must match len(IFOs) and/or the data in sftfilepattern. Detectors will be paired to list elements following alphabetical order. Required if sftfilepattern=None, optional otherwise..
tempory_filename (str) – Temporary file used for lalapps_PredictFstat output, will be deleted at the end.
earth_ephem (str or None) – Ephemerides files, defaults will be used if None.
sun_ephem (str or None) – Ephemerides files, defaults will be used if None.
transientWindowType (str) – Optional parameter for transient signals, see lalapps_PredictFstat –help. Default of none means a classical Continuous Wave signal.
transientStartTime (int or None) – Optional parameters for transient signals, see lalapps_PredictFstat –help.
transientTau (int or None) – Optional parameters for transient signals, see lalapps_PredictFstat –help.
- Returns
twoF_expected, twoF_sigma – The expectation and standard deviation of 2F.
- Return type
float
- pyfstat.helper_functions.parse_list_of_numbers(val)[source]
Convert a number, list of numbers or comma-separated str into a list of numbers.
This is useful e.g. for sqrtSX inputs where the user can be expected to try either type of input.
- Parameters
val (float, list or str) – The input to be parsed.
- Returns
out – The parsed list.
- Return type
list
- pyfstat.helper_functions.generate_loudest_file(max_params, tref, outdir, label, sftfilepattern, minStartTime=None, maxStartTime=None, transientWindowType=None, earth_ephem=None, sun_ephem=None)[source]
Use lalapps_ComputeFstatistic_v2 to produce a .loudest file.
- Parameters
max_params (dict) – Dictionary of a single parameter-space point. This needs to already have been translated to lal conventions and must NOT include detection statistic entries!
tref (int) – Reference time of the max_params.
outdir (str) – Directory to place the .loudest file in.
label (str) – Search name bit to be used in the output filename.
sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.
minStartTime (int or None) – GPS seconds of the start time and end time; default: use al available data.
maxStartTime (int or None) – GPS seconds of the start time and end time; default: use al available data.
transientWindowType (str or None) – optional: transient window type, needs to go with t0 and tau parameters inside max_params.
earth_ephem (str or None) – optional: user-set Earth ephemeris file
sun_ephem (str or None) – optional: user-set Sun ephemeris file
- Returns
loudest_file – The filename of the CFSv2 output file.
- Return type
str
- pyfstat.helper_functions.gps_to_datestr_utc(gps)[source]
Convert an integer count of GPS seconds to a UTC date-time string.
This uses the locale’s default string formatting as per datetime.strftime(). It is intended just for informing the user and may not be as reliable in all situations as lal[apps]_tconvert. If you want to do any postprocessing of the date-time string, for safety you should probably call that commandline tool.
- Parameters
gps (int) – Integer seconds since GPS seconds.
- Returns
dtstr – A string representation of date-time in UTC and locale format.
- Return type
str
pyfstat.injection_parameters module
Generate injection parameters drawn from different prior populations
- class pyfstat.injection_parameters.InjectionParametersGenerator(*, seed: Union[None, int] = None, rng: numpy.random._generator.Generator = NOTHING, priors: dict = NOTHING)[source]
Bases:
object
Draw injection parameter samples from priors and return in dictionary format.
- Parameters
priors (dict) –
Each key refers to one of the signal’s parameters (following the PyFstat convention). Priors can be given as values in three formats (by order of evaluation):
1. Callable without required arguments: {“ParameterA”: np.random.uniform}.
2. Dict containing numpy.random distribution as key and kwargs in a dict as value: {“ParameterA”: {“uniform”: {“low”: 0, “high”:1}}}.
3. Constant value to be returned as is: {“ParameterA”: 1.0}.
seed (None, int) – Argument to be fed to numpy.random.default_rng, with all of its accepted types.
_rng (np.random.Generator) – Alternatively, this class accepts an already-initialized np.Generator, in which case the seed argument will be ignored.
Method generated by attrs for class InjectionParametersGenerator.
- seed: Union[None, int]
- priors: dict
- class pyfstat.injection_parameters.AllSkyInjectionParametersGenerator(*, seed: Union[None, int] = None, rng: numpy.random._generator.Generator = NOTHING, priors: dict = NOTHING)[source]
Bases:
pyfstat.injection_parameters.InjectionParametersGenerator
Draw injection parameter samples from priors and return in dictionary format. This class works in exactly the same way as InjectionParametersGenerator, but including by default two extra keys, Alpha and Delta (sky position’s right ascension and declination in radians), which are sample isotropically across the celesetial sphere.
Alpha’s distribution is Uniform(0, 2 pi), and sin(Delta)’s distribution is Uniform(-1, 1).
The only reason this class exists is because, using the notation we specified in the base class, there is no way to generate arcsine distributed numbers using a specific seed, as numpy does not have such a number generator and hence has to be constructed by applying a function to a uniform number.
Method generated by attrs for class InjectionParametersGenerator.
- seed: Union[None, int]
- priors: dict
pyfstat.make_sfts module
PyFstat tools to generate and manipulate data in the form of SFTs.
- class pyfstat.make_sfts.Writer(*args, **kwargs)[source]
Bases:
pyfstat.core.BaseSearchClass
The main class for generating data in the form of SFTs.
Short Fourier Transforms (SFTs) are a standard data format used in LALSuite, containing the Fourier transform of strain data over a duration Tsft.
SFT data can be generated from scratch, including Gaussian noise and/or simulated CW signals or transient signals. Existing SFTs (real data or previously simulated) can also be reused through the noiseSFTs option, allowing to ‘inject’ additional signals into them.
This class currently relies on the lalapps_Makefakedata_v5 executable which will be run in a subprocess. See lalapps_Makefakedata_v5 –help for more detailed help with some of the parameters.
- Parameters
label (string) – A human-readable label to be used in naming the output files.
tstart (int) – Starting GPS epoch of the data set. NOTE: mutually exclusive with timestamps.
duration (int) – Duration (in GPS seconds) of the total data set. NOTE: mutually exclusive with timestamps.
tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.
F0 (float or None) – Frequency of a signal to inject. Also used (if Band is not None) as center of frequency band. Also needed when noise-only (h0=None or h0==0) but no noiseSFTs given, in which case it is also used as center of frequency band.
F1 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
F2 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Alpha (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Delta (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
h0 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
cosi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
psi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
phi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.
outdir (str) – The directory where files are written to. Default: current working directory.
sqrtSX (float or list or str or None) – Single-sided PSD values for generating fake Gaussian noise. Single float or str value: use same for all detectors. List or comma-separated string: must match len(detectors). Detectors will be paired to list elements following alphabetical order.
noiseSFTs (str or None) – Existing SFT files on top of which signals will be injected. If not None, additional constraints can be applied using the arguments tstart and duration. NOTE: mutually exclusive with timestamps.
SFTWindowType (str or None) – LAL name of the windowing function to apply to the data.
SFTWindowBeta (float) – Optional parameter for some windowing functions.
Band (float or None) – If float, and F0 is also not None, then output SFTs cover [F0-Band/2,F0+Band/2]. If None and noiseSFTs given, use their bandwidth. If None and no noiseSFTs given, a minimal covering band for a perfectly-matched single-template ComputeFstat analysis is estimated.
detectors (str or None) – Comma-separated list of detectors to generate data for.
earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
sun_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
transientWindowType (str) – If none, a fully persistent CW signal is simulated. If rect or exp, a transient signal with the corresponding amplitude evolution is simulated.
transientStartTime (int or None) – Start time for a transient signal.
transientTau (int or None) – Duration (rect case) or decay time (exp case) of a transient signal.
randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.
timestamps (str or dict) – Dictionary of timestamps (each key must refer to a detector), list of timestamps (detectors should be set), or comma-separated list of per-detector timestamps files (simple text files, lines with <seconds> <nanoseconds>, comments should use %, each time stamp gives the start time of one SFT). WARNING: In that last case, order must match that of detectors! NOTE: mutually exclusive with [tstart,`duration`] and with noiseSFTs.
- mfd = 'lalapps_Makefakedata_v5'
The executable; can be overridden by child classes.
- signal_parameter_labels = ['tref', 'F0', 'F1', 'F2', 'Alpha', 'Delta', 'h0', 'cosi', 'psi', 'phi', 'transientWindowType', 'transientStartTime', 'transientTau']
Default convention of labels for the various signal parameters.
- gps_time_and_string_formats_as_LAL = {'refTime': ':10.9f', 'transientStartTime': ':10.0f', 'transientTau': ':10.0f', 'transientWindowType': ':s'}
Dictionary to ensure proper format handling for some special parameters.
GPS times should NOT be parsed using scientific notation. LAL routines would silently parse them wrongly.
- required_signal_parameters = ['F0', 'Alpha', 'Delta', 'cosi']
List of parameters required for a successful execution of Makefakedata_v5. The rest of available parameters are not required as they have default values silently given by Makefakedata_v5
- property tend
` Defined as self.start + self.duration.
If stored as an attribute, there would be the risk of it going out of sync with the other two values.
- calculate_fmin_Band()[source]
Set fmin and Band for the output SFTs to cover.
Either uses the user-provided Band and puts F0 in the middle; does nothing to later reuse the full bandwidth of noiseSFTs (only if using MFDv5); or if F0!=None, noiseSFTs=None and Band=None it estimates a minimal band for just the injected signal: F-stat covering band plus extra bins for demod default parameters. This way a perfectly matched single-template ComputeFstat analysis should run through perfectly on the returned SFTs. For any wider-band or mismatched search, one needs to set Band manually. If using MFDv4, at least F0 is required even if noiseSFTs!=None.
- make_cff(verbose=False)[source]
Generates a .cff file including signal injection parameters.
This will be saved to self.config_file_name.
- Parameters
verbose (boolean) – If true, increase logging verbosity.
- check_cached_data_okay_to_use(cl_mfd)[source]
Check if SFT files already exist that can be re-used.
This does not check the actual data contents of the SFTs, but only the following criteria:
filename
if injecting a signal, that the .cff file is older than the SFTs (but its contents should have been checked separately)
that the commandline stored in the (first) SFT header matches
- Parameters
cl_mfd (str) – The commandline we’d execute if not finding matching files.
- run_makefakedata()[source]
Generate the SFT data calling lalapps_Makefakedata_v5.
This first builds the full commandline, then calls check_cached_data_okay_to_use() to see if equivalent data files already exist, and else runs the actual generation code.
- predict_fstat(assumeSqrtSX=None)[source]
Predict the expected F-statistic value for the injection parameters.
Through helper_functions.predict_fstat(), this wraps the lalapps_PredictFstat executable.
- Parameters
assumeSqrtSX (float, str or None) – If None, PSD is estimated from self.sftfilepath. Else, assume this stationary per-detector noise-floor instead. Single float or str value: use same for all IFOs. Comma-separated string: must match len(self.detectors) and the data in self.sftfilepath. Detectors will be paired to list elements following alphabetical order.
- class pyfstat.make_sfts.BinaryModulatedWriter(*args, **kwargs)[source]
Bases:
pyfstat.make_sfts.Writer
Special Writer variant for simulating a CW signal for a source in a binary system.
Most parameters are the same as for the basic Writer class, only the additional ones are documented here:
- Parameters
tp – binary orbit parameters
argp – binary orbit parameters
asini – binary orbit parameters
ecc – binary orbit parameters
period – binary orbit parameters
- class pyfstat.make_sfts.LineWriter(*args, **kwargs)[source]
Bases:
pyfstat.make_sfts.Writer
Inject a simulated line-like detector artifact into SFT data.
A (transient) line is defined as a constant amplitude and constant excess power artifact in the data.
In practice, it corresponds to a CW without Doppler or antenna-patern-induced amplitude modulation.
NOTE: This functionality is implemented via lalapps_MakeFakeData_v4’s lineFeature option. This version of MFD only supports one interferometer at a time.
NOTE: All signal parameters except for h0, Freq, phi0 and transient parameters will be ignored.
- Parameters
label (string) – A human-readable label to be used in naming the output files.
tstart (int) – Starting GPS epoch of the data set. NOTE: mutually exclusive with timestamps.
duration (int) – Duration (in GPS seconds) of the total data set. NOTE: mutually exclusive with timestamps.
tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.
F0 (float or None) – Frequency of a signal to inject. Also used (if Band is not None) as center of frequency band. Also needed when noise-only (h0=None or h0==0) but no noiseSFTs given, in which case it is also used as center of frequency band.
F1 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
F2 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Alpha (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Delta (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
h0 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
cosi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
psi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
phi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.
Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.
outdir (str) – The directory where files are written to. Default: current working directory.
sqrtSX (float or list or str or None) – Single-sided PSD values for generating fake Gaussian noise. Single float or str value: use same for all detectors. List or comma-separated string: must match len(detectors). Detectors will be paired to list elements following alphabetical order.
noiseSFTs (str or None) – Existing SFT files on top of which signals will be injected. If not None, additional constraints can be applied using the arguments tstart and duration. NOTE: mutually exclusive with timestamps.
SFTWindowType (str or None) – LAL name of the windowing function to apply to the data.
SFTWindowBeta (float) – Optional parameter for some windowing functions.
Band (float or None) – If float, and F0 is also not None, then output SFTs cover [F0-Band/2,F0+Band/2]. If None and noiseSFTs given, use their bandwidth. If None and no noiseSFTs given, a minimal covering band for a perfectly-matched single-template ComputeFstat analysis is estimated.
detectors (str or None) – Comma-separated list of detectors to generate data for.
earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
sun_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
transientWindowType (str) – If none, a fully persistent CW signal is simulated. If rect or exp, a transient signal with the corresponding amplitude evolution is simulated.
transientStartTime (int or None) – Start time for a transient signal.
transientTau (int or None) – Duration (rect case) or decay time (exp case) of a transient signal.
randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.
timestamps (str or dict) – Dictionary of timestamps (each key must refer to a detector), list of timestamps (detectors should be set), or comma-separated list of per-detector timestamps files (simple text files, lines with <seconds> <nanoseconds>, comments should use %, each time stamp gives the start time of one SFT). WARNING: In that last case, order must match that of detectors! NOTE: mutually exclusive with [tstart,`duration`] and with noiseSFTs.
- mfd = 'lalapps_Makefakedata_v4'
The executable (older version that supports the –lineFeature option).
- required_signal_parameters = ['F0', 'phi', 'h0']
Required parameters for Makefakedata_v4 to success. Any other parameter is silently given a default value by Makefakedata_v4.
- signal_parameters_labels = ['F0', 'phi', 'h0', 'transientWindowType', 'transientStartTime', 'transientTau']
Other signal parameters will be removed before passing to Makefakedata_v4.
- class pyfstat.make_sfts.GlitchWriter(*args, **kwargs)[source]
Bases:
pyfstat.core.SearchForSignalWithJumps
,pyfstat.make_sfts.Writer
Special Writer variant for simulating a CW signal containing a timing glitch.
Most parameters are the same as for the basic Writer class, only the additional ones are documented here:
- Parameters
dtglitch (float or None) – Time (in GPS seconds) of the glitch after tstart. To create data without a glitch, set dtglitch=None.
delta_phi (float) – Instantaneous glitch magnitudes in rad, Hz, and Hz/s respectively.
delta_F0 (float) – Instantaneous glitch magnitudes in rad, Hz, and Hz/s respectively.
delta_F1 (float) – Instantaneous glitch magnitudes in rad, Hz, and Hz/s respectively.
- class pyfstat.make_sfts.FrequencyModulatedArtifactWriter(*args, **kwargs)[source]
Bases:
pyfstat.make_sfts.Writer
Specialized Writer variant to generate SFTs containing simulated instrumental artifacts.
Contrary to the main Writer class, this calls the older lalapps_Makefakedata_v4 executable which supports the special –lineFeature option. See lalapps_Makefakedata_v4 –help for more detailed help with some of the parameters.
- Parameters
label (string) – A human-readable label to be used in naming the output files.
outdir (str) – The directory where files are written to. Default: current working directory.
tstart (int) – Starting GPS epoch of the data set.
duration (int) – Duration (in GPS seconds) of the total data set.
F0 (float) – Frequency of the artifact.
F1 (float) – Frequency drift of the artifact.
tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.
h0 (float) – Amplitude of the artifact.
Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.
sqrtSX (float) – Background detector noise level.
Band (float) – Output SFTs cover [F0-Band/2,F0+Band/2].
Pmod (float) – Modulation period of the artifact.
Pmod_phi (float) – Additional parameters for modulation of the artifact.
Pmod_amp (float) – Additional parameters for modulation of the artifact.
Alpha (float or None) – If not none: add an orbital modulation to the artifact corresponding to a signal from that sky position, in radians.
Delta (float or None) – If not none: add an orbital modulation to the artifact corresponding to a signal from that sky position, in radians.
detectors (str or None) – Comma-separated list of detectors to generate data for.
earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
sun_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.
- get_frequency(t)[source]
Evolve the artifact frequency in time.
This includes a drift term and optionally, if Alpha and Delta are not None, a simulated orbital modulation.
- Parameters
t (float) – Time stamp to evaluate the frequency at.
- Returns
f – Frequency at time t.
- Return type
float
- get_h0(t)[source]
Evaluate the artifact amplitude at a given time.
NOTE: Here it’s actually implemented as a constant!
- Parameters
t (float) – Time stamp to evaluate at.
- Returns
h0 – Amplitude at time t.
- Return type
float
- pre_compute_evolution()[source]
Precomputes evolution parameters for the artifact.
This computes midtimes, frequencies, phases and amplitudes over the list of SFT timestamps.
- class pyfstat.make_sfts.FrequencyAmplitudeModulatedArtifactWriter(*args, **kwargs)[source]
Bases:
pyfstat.make_sfts.FrequencyModulatedArtifactWriter
A variant of FrequencyModulatedArtifactWriter with evolving amplitude.
- Parameters
label (string) – A human-readable label to be used in naming the output files.
outdir (str) – The directory where files are written to. Default: current working directory.
tstart (int) – Starting GPS epoch of the data set.
duration (int) – Duration (in GPS seconds) of the total data set.
F0 (float) – Frequency of the artifact.
F1 (float) – Frequency drift of the artifact.
tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.
h0 (float) – Amplitude of the artifact.
Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.
sqrtSX (float) – Background detector noise level.
Band (float) – Output SFTs cover [F0-Band/2,F0+Band/2].
Pmod (float) – Modulation period of the artifact.
Pmod_phi (float) – Additional parameters for modulation of the artifact.
Pmod_amp (float) – Additional parameters for modulation of the artifact.
Alpha (float or None) – If not none: add an orbital modulation to the artifact corresponding to a signal from that sky position, in radians.
Delta (float or None) – If not none: add an orbital modulation to the artifact corresponding to a signal from that sky position, in radians.
detectors (str or None) – Comma-separated list of detectors to generate data for.
earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
sun_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per helper_functions.get_ephemeris_files().
randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.
- class pyfstat.make_sfts.InjectionParametersGenerator(*args, **kwargs)[source]
Bases:
pyfstat.injection_parameters.InjectionParametersGenerator
Method generated by attrs for class InjectionParametersGenerator.
- seed: Union[None, int]
- priors: dict
- class pyfstat.make_sfts.AllSkyInjectionParametersGenerator(*args, **kwargs)[source]
Bases:
pyfstat.injection_parameters.AllSkyInjectionParametersGenerator
Method generated by attrs for class InjectionParametersGenerator.
- seed: Union[None, int]
- priors: dict
pyfstat.mcmc_based_searches module
PyFstat search & follow-up classes using MCMC-based methods
The general approach is described in Ashton & Prix (PRD 97, 103020, 2018): https://arxiv.org/abs/1802.05450 and we use the ptemcee sampler described in Vousden et al. (MNRAS 455, 1919-1937, 2016): https://arxiv.org/abs/1501.05823 and based on Foreman-Mackey et al. (PASP 125, 306, 2013): https://arxiv.org/abs/1202.3665
Defining the prior
The MCMC based searches (i.e. pyfstat.MCMC*) require a prior specification for each model parameter, implemented via a python dictionary. This is best explained through a simple example, here is the prior for a directed search with a uniform prior on the frequency and a normal prior on the frequency derivative:
theta_prior = {'F0': {'type': 'unif',
'lower': 29.9,
'upper': 30.1},
'F1': {'type': 'norm',
'loc': 0,
'scale': 1e-10},
'F2': 0,
'Alpha': 2.3,
'Delta': 1.8
}
For the sky positions Alpha
and Delta
, we give the fixed values (i.e. they are considered known by
the MCMC simulation), the same is true for F2
, the second derivative of the frequency which we fix at 0
.
Meanwhile, for the frequency F0
and first frequency derivative F1
we give a dictionary specifying their
prior distribution. This dictionary must contain three arguments: the type
(in this case either unif
or
norm
) which specifies the type of distribution, then two shape arguments. The shape parameters will depend
on the type
of distribution, but here we use lower
and upper
, required for the unif
prior while
loc
and scale
are required for the norm
prior.
Currently, two other types of prior are implemented: halfnorm
, neghalfnorm
(both of which require loc
and scale
shape parameters). Further priors can be added by modifying pyfstat.MCMCSearch._generic_lnprior
.
- class pyfstat.mcmc_based_searches.MCMCSearch(*args, **kwargs)[source]
Bases:
pyfstat.core.BaseSearchClass
MCMC search using ComputeFstat.
Evaluates the coherent F-statistic across a parameter space region corresponding to an isolated/binary-modulated CW signal.
- Parameters
theta_prior (dict) – Dictionary of priors and fixed values for the search parameters. For each parameters (key of the dict), if it is to be held fixed the value should be the constant float, if it is be searched, the value should be a dictionary of the prior.
tref (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
minStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
maxStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
label (str) – A label and output directory (optional, default is data) to name files
outdir (str) – A label and output directory (optional, default is data) to name files
sftfilepattern (str, optional) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.
detectors (str, optional) – Two character reference to the detectors to use, specify None for no contraint and comma separated strings for multiple references.
nsteps (list (2,), optional) – Number of burn-in and production steps to take, [nburn, nprod]. See pyfstat.MCMCSearch.setup_initialisation() for details on adding initialisation steps.
nwalkers (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.
ntemps (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.
log10beta_min (float < 0, optional) – The log_10(beta) value. If given, the set of betas passed to PTSampler are generated from np.logspace(0, log10beta_min, ntemps) (given in descending order to ptemcee).
theta_initial (dict, array, optional) – A dictionary of distribution about which to distribute the initial walkers about.
rhohatmax (float, optional) – Upper bound for the SNR scale parameter (required to normalise the Bayes factor) - this needs to be carefully set when using the evidence.
binary (bool, optional) – If true, search over binary orbital parameters.
BSGL (bool, optional) – If true, use the BSGL statistic.
BtSG (bool, optional) – If true, use the transient lnBtSG statistic. (Only for transient searches.)
SSBPrec (int, optional) – SSBPrec (SSB precision) to use when calling ComputeFstat. See core.ComputeFstat.
RngMedWindow (int, optional) – Running-Median window size (number of bins) for ComputeFstat. See core.ComputeFstat.
minCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.
maxCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.
injectSources (dict, optional) – If given, inject these properties into the SFT files before running the search. See core.ComputeFstat.
assumeSqrtSX (float or list or str) – Don’t estimate noise-floors, but assume (stationary) per-IFO sqrt{SX}. See core.ComputeFstat.
transientWindowType (str) – If ‘rect’ or ‘exp’, compute atoms so that a transient (t0,tau) map can later be computed. (‘none’ instead of None explicitly calls the transient-window function, but with the full range, for debugging). See core.ComputeFstat. Currently only supported for nsegs=1.
tCWFstatMapVersion (str) – Choose between standard ‘lal’ implementation, ‘pycuda’ for gpu, and some others for devel/debug.
allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength].
- symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'argp': 'argp', 'asini': 'asini', 'ecc': 'ecc', 'period': 'P', 'tp': 'tp'}
Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots
- unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'argp': '', 'asini': '', 'ecc': '', 'period': 's', 'tp': 's'}
Key, val pairs of the parameters (i.e. F0, F1), and the units (i.e. Hz)
- transform_dictionary = {}
Key, val pairs of the parameters (i.e. F0, F1), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.
- setup_initialisation(nburn0, scatter_val=1e-10)[source]
Add an initialisation step to the MCMC run
If called prior to run(), adds an intial step in which the MCMC simulation is run for nburn0 steps. After this, the MCMC simulation continues in the usual manner (i.e. for nburn and nprod steps), but the walkers are reset scattered around the maximum likelihood position of the initialisation step.
- Parameters
nburn0 (int) – Number of initialisation steps to take.
scatter_val (float) – Relative number to scatter walkers around the maximum likelihood position after the initialisation step. If the maximum likelihood point is located at p, the new walkers are randomly drawn from a multivariate gaussian distribution centered at p with standard deviation diag(scatter_val * p).
- run(proposal_scale_factor=2, save_pickle=True, export_samples=True, save_loudest=True, plot_walkers=True, walker_plot_args=None, window=50)[source]
Run the MCMC simulatation
- Parameters
proposal_scale_factor (float) – The proposal scale factor a > 1 used by the sampler. See Goodman & Weare (Comm App Math Comp Sci, Vol 5, No. 1, 2010): 10.2140/camcos.2010.5.65. The bigger the value, the wider the range to draw proposals from. If the acceptance fraction is too low, you can raise it by decreasing the a parameter; and if it is too high, you can reduce it by increasing the a parameter. See Foreman-Mackay et al. (PASP 125 306, 2013): https://arxiv.org/abs/1202.3665.
save_pickle (bool) – If true, save a pickle file of the full sampler state.
export_samples (bool) – If true, save ASCII samples file to disk. See MCMCSearch.export_samples_to_disk.
save_loudest (bool) – If true, save a CFSv2 .loudest file to disk. See MCMCSearch.generate_loudest.
plot_walkers (bool) – If true, save trace plots of the walkers.
walker_plot_args – Dictionary passed as kwargs to _plot_walkers to control the plotting. Histogram of sampled detection statistic values can be retrieved setting “plot_det_stat” to True. Parameters corresponding to an injected signal can be passed through “injection_parameters” as a dictionary containing the parameters of said signal. All parameters being searched for must be present, otherwise this option is ignored. If both “fig” and “axes” entries are set, the plot is not saved to disk directly, but (fig, axes) are returned.
window (int) – The minimum number of autocorrelation times needed to trust the result when estimating the autocorrelation time (see ptemcee.Sampler.get_autocorr_time for further details.
- plot_corner(figsize=(10, 10), add_prior=False, nstds=None, label_offset=0.4, dpi=300, rc_context={}, tglitch_ratio=False, fig_and_axes=None, save_fig=True, **kwargs)[source]
Generate a corner plot of the posterior
Using the corner package (https://pypi.python.org/pypi/corner/), generate estimates of the posterior from the production samples.
- Parameters
figsize (tuple (7, 7)) – Figure size in inches (passed to plt.subplots)
add_prior (bool, str) – If true, plot the prior as a red line. If ‘full’ then for uniform priors plot the full extent of the prior.
nstds (float) – The number of standard deviations to plot centered on the median. Standard deviation is computed from the samples using numpy.std.
label_offset (float) – Offset the labels from the plot: useful to prevent overlapping the tick labels with the axis labels. This option is passed to ax.[x|y]axis.set_label_coords.
dpi (int) – Passed to plt.savefig.
rc_context (dict) – Dictionary of rc values to set while generating the figure (see matplotlib rc for more details).
tglitch_ratio (bool) – If true, and tglitch is a parameter, plot posteriors as the fractional time at which the glitch occurs instead of the actual time.
fig_and_axes (tuple) – (fig, axes) tuple to plot on. The axes must be of the right shape, namely (ndim, ndim)
save_fig (bool) – If true, save the figure, else return the fig, axes.
**kwargs – Passed to corner.corner. Use “truths” to plot the true parameters of a signal.
- Returns
The matplotlib figure and axes, only returned if save_fig = False.
- Return type
fig, axes
- plot_chainconsumer(save_fig=True, label_offset=0.25, dpi=300, **kwargs)[source]
Generate a corner plot of the posterior using the chaniconsumer package.
chainconsumer is an optional dependency of PyFstat. See https://samreay.github.io/ChainConsumer/.
Parameters are akin to the ones described in MCMCSearch.plot_corner. Only the differing parameters are explicitly described.
- Parameters
**kwargs – Passed to chainconsumer.plotter.plot. Use “truths” to plot the true parameters of a signal.
- plot_prior_posterior(normal_stds=2, injection_parameters=None, fig_and_axes=None, save_fig=True)[source]
Plot the prior and posterior probability distributions in the same figure
- Parameters
normal_stds (int) – Bounds of priors in terms of their standard deviation. Only used if norm, halfnorm, neghalfnorm or lognorm priors are given, otherwise ignored.
injection_parameters (dict) – Dictionary containing the parameters of a signal. All parameters being searched must be present as dictionary keys, otherwise this option is ignored.
fig_and_axes (tuple) – (fig, axes) tuple to plot on.
save_fig (bool) – If true, save the figure, else return the fig, axes.
- Returns
(fig, ax) – If save_fig evaluates to False, return figure and axes.
- Return type
(matplotlib.pyplot.figure, matplotlib.pyplot.axes)
- plot_cumulative_max(**kwargs)[source]
Plot the cumulative twoF for the maximum posterior estimate.
This method accepts the same arguments as pyfstat.core.ComputeFstat.plot_twoF_cumulative, except for CFS_input, which is taken from the loudest candidate; and label and outdir, which are taken from the instance of this class.
For example, one can pass signal arguments to predic_twoF_cumulative through PFS_kwargs, or set the number of segments using num_segments_(CFS|PFS). The same applies for other options such as tstart, tend or savefig. Every single of these arguments will be passed to pyfstat.core.ComputeFstat.plot_twoF_cumulative as they are, using their default argument otherwise.
See pyfstat.core.ComputeFstat.plot_twoF_cumulative for a comprehensive list of accepted arguments and their default values.
Unlike the core function, here savefig=True is the default, for consistency with other MCMC plotting functions.
- get_saved_data_dictionary()[source]
Read the data saved in self.pickel_path and return it as a dictionary.
- Returns
d – Dictionary containing the data saved in the pickle self.pickle_path.
- Return type
dict
- get_max_twoF()[source]
Get the max. likelihood (loudest) sample and the compute its corresponding detection statistic.
The employed detection statistic depends on self.BSGL (i.e. 2F if self.BSGL evaluates to False, log10BSGL otherwise).
- Returns
d (dict) – Parameters of the loudest sample.
maxtwoF (float) – Detection statistic (2F or log10BSGL) corresponding to the loudest sample.
- get_summary_stats()[source]
Returns a dict of point estimates for all production samples.
Point estimates are computed on the MCMC samples using numpy.mean, numpy.std and numpy.quantiles with q=[0.005, 0.05, 0.25, 0.5, 0.75, 0.95, 0.995].
- Returns
d – Dictionary containing point estimates corresponding to [“mean”, “std”, “lower99”, “lower90”, “lower50”, “median”, “upper50”, “upper90”, “upper99”].
- Return type
dict
- check_if_samples_are_railing(threshold=0.01)[source]
Returns a boolean estimate of if the samples are railing
- Parameters
threshold (float [0, 1]) – Fraction of the uniform prior to test (at upper and lower bound)
- Returns
return_flag – IF true, the samples are railing
- Return type
bool
- write_par(method='median')[source]
Writes a .par of the best-fit params with an estimated std
- Parameters
method (str) – How to select the best-fit params. Available methods: “median”, “mean”, “twoFmax”.
- get_p_value(delta_F0=0, time_trials=0)[source]
Gets the p-value for the maximum twoFhat value assuming Gaussian noise
- Parameters
delta_F0 (float) – Frequency variation due to a glitch.
time_trials (int, optional) – Number of trials in each glitch + 1.
- compute_evidence(make_plots=False, write_to_file=None)[source]
Computes the evidence/marginal likelihood for the model.
- Parameters
make_plots (bool) – Plot the results and save them to os.path.join(self.outdir, self.label + “_beta_lnl.png”)
write_to_file (str) – If given, dump evidence and uncertainty estimation to the specified path.
- Returns
log10evidence (float) – Estimation of the log10 evidence.
log10evidence_err (float) – Log10 uncertainty of the evidence estimation.
- static read_evidence_file_to_dict(evidence_file_name='Evidences.txt')[source]
Read evidence file and put it into an OrderedDict
An evidence file contains paris (log10evidence, log10evidence_err) for each considered model. These pairs are prepended by the self.label variable.
- Parameters
evidence_file_name (str) – Filename to read.
- Returns
EvidenceDict – Dictionary with the contents of evidence_file_name
- Return type
dict
- class pyfstat.mcmc_based_searches.MCMCGlitchSearch(*args, **kwargs)[source]
Bases:
pyfstat.mcmc_based_searches.MCMCSearch
MCMC search using the SemiCoherentGlitchSearch
See parent MCMCSearch for a list of all additional parameters, here we list only the additional init parameters of this class.
- Parameters
nglitch (int) – The number of glitches to allow
dtglitchmin (int) – The minimum duration (in seconds) of a segment between two glitches or a glitch and the start/end of the data
theta0_idx – Index (zero-based) of which segment the theta refers to - useful if providing a tight prior on theta to allow the signal to jump too theta (and not just from)
int – Index (zero-based) of which segment the theta refers to - useful if providing a tight prior on theta to allow the signal to jump too theta (and not just from)
- symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'delta_F0': '$\\delta f$', 'delta_F1': '$\\delta \\dot{f}$', 'tglitch': '$t_\\mathrm{glitch}$'}
Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots
- glitch_symbol_dictionary = {'delta_F0': '$\\delta f$', 'delta_F1': '$\\delta \\dot{f}$', 'tglitch': '$t_\\mathrm{glitch}$'}
Key, val pairs of glitch parameters (dF0, dF1, tglitch), to LaTeX math symbols for plots. This dictionary included within self.symbol_dictionary.
- unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'delta_F0': 'Hz', 'delta_F1': 'Hz/s', 'tglitch': 's'}
Key, val pairs of the parameters (F0, F1, …, including glitch parameters), and the units (Hz, Hz/s, …).
- transform_dictionary = {'tglitch': {'label': '$t^{g}_0$ \\n [d]', 'multiplier': 1.1574074074074073e-05, 'subtractor': 'minStartTime', 'unit': 'day'}}
Key, val pairs of the parameters (F0, F1, …), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.
- plot_cumulative_max(savefig=True)[source]
Override MCMCSearch.plot_cumulative_max implementation to deal with the split at glitches.
- Parameters
savefig (boolean) – included for consistency with core plot_twoF_cumulative() function. If true, save the figure in outdir. If false, return an axis object.
- class pyfstat.mcmc_based_searches.MCMCSemiCoherentSearch(*args, **kwargs)[source]
Bases:
pyfstat.mcmc_based_searches.MCMCSearch
MCMC search for a signal using the semicoherent ComputeFstat.
Evaluates the semicoherent F-statistic acros a parameter space region corresponding to an isolated/binary-modulated CW signal.
See MCMCSearch for a list of additional parameters, here we list only the additional init parameters of this class.
- Parameters
nsegs (int) – The number of segments into which the input datastream will be devided. Coherence time is computed internally as (maxStartTime - minStarTime) / nsegs.
- class pyfstat.mcmc_based_searches.MCMCFollowUpSearch(*args, **kwargs)[source]
Bases:
pyfstat.mcmc_based_searches.MCMCSemiCoherentSearch
Hierarchical follow-up procedure
Executes MCMC runs with increasing coherence times in order to follow up a parameter space region. The main idea is to use an MCMC run to identify an interesting parameter space region to then zoom-in said region using a finer “effective resolution” by increasing the coherence time. See Ashton & Prix (PRD 97, 103020, 2018): https://arxiv.org/abs/1802.05450
See MCMCSemiCoherentSearch for a list of additional parameters, here we list only the additional init parameters of this class.
- Parameters
nsegs (int) – The number of segments into which the input datastream will be devided. Coherence time is computed internally as (maxStartTime - minStarTime) / nsegs.
- run(run_setup=None, proposal_scale_factor=2, NstarMax=10, Nsegs0=None, save_pickle=True, export_samples=True, save_loudest=True, plot_walkers=True, walker_plot_args=None, log_table=True, gen_tex_table=True, window=50)[source]
Run the follow-up with the given run_setup.
See MCMCSearch.run’s docstring for a description of the remainder arguments.
- Parameters
run_setup – See MCMCFollowUpSearch.init_run_setup.
log_table – See MCMCFollowUpSearch.init_run_setup.
gen_tex_table – See MCMCFollowUpSearch.init_run_setup.
NstarMax – See pyfstat.optimal_setup_functions.get_optimal_setup.
Nsegs0 – See pyfstat.optimal_setup_functions.get_optimal_setup.
- init_run_setup(run_setup=None, NstarMax=1000, Nsegs0=None, log_table=True, gen_tex_table=True)[source]
Initialize the setup of the follow-up run computing the required quantities fro, NstarMax and Nsegs0.
- Parameters
NstarMax (int) – Required parameters to create a new follow-up setup. See pyfstat.optimal_setup_functions.get_optimal_setup.
Nsegs0 (int) – Required parameters to create a new follow-up setup. See pyfstat.optimal_setup_functions.get_optimal_setup.
run_setup (optional) – If None, a new setup will be created from NstarMax and Nsegs0. Use MCMCFollowUpSearch.read_setup_input_file to read a previous setup file.
log_table (bool) – Log follow-up setup using logging.info as a table.
gen_tex_table (bool) – Dump follow-up setup into a text file as a tex table. File is constructed as os.path.join(self.outdir, self.label + “_run_setup.tex”).
- Returns
run_setup – List containing the setup of the follow-up run.
- Return type
list
- class pyfstat.mcmc_based_searches.MCMCTransientSearch(*args, **kwargs)[source]
Bases:
pyfstat.mcmc_based_searches.MCMCSearch
MCMC search for a transient signal using ComputeFstat
- Parameters
theta_prior (dict) – Dictionary of priors and fixed values for the search parameters. For each parameters (key of the dict), if it is to be held fixed the value should be the constant float, if it is be searched, the value should be a dictionary of the prior.
tref (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
minStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
maxStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.
label (str) – A label and output directory (optional, default is data) to name files
outdir (str) – A label and output directory (optional, default is data) to name files
sftfilepattern (str, optional) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.
detectors (str, optional) – Two character reference to the detectors to use, specify None for no contraint and comma separated strings for multiple references.
nsteps (list (2,), optional) – Number of burn-in and production steps to take, [nburn, nprod]. See pyfstat.MCMCSearch.setup_initialisation() for details on adding initialisation steps.
nwalkers (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.
ntemps (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.
log10beta_min (float < 0, optional) – The log_10(beta) value. If given, the set of betas passed to PTSampler are generated from np.logspace(0, log10beta_min, ntemps) (given in descending order to ptemcee).
theta_initial (dict, array, optional) – A dictionary of distribution about which to distribute the initial walkers about.
rhohatmax (float, optional) – Upper bound for the SNR scale parameter (required to normalise the Bayes factor) - this needs to be carefully set when using the evidence.
binary (bool, optional) – If true, search over binary orbital parameters.
BSGL (bool, optional) – If true, use the BSGL statistic.
BtSG (bool, optional) – If true, use the transient lnBtSG statistic. (Only for transient searches.)
SSBPrec (int, optional) – SSBPrec (SSB precision) to use when calling ComputeFstat. See core.ComputeFstat.
RngMedWindow (int, optional) – Running-Median window size (number of bins) for ComputeFstat. See core.ComputeFstat.
minCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.
maxCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.
injectSources (dict, optional) – If given, inject these properties into the SFT files before running the search. See core.ComputeFstat.
assumeSqrtSX (float or list or str) – Don’t estimate noise-floors, but assume (stationary) per-IFO sqrt{SX}. See core.ComputeFstat.
transientWindowType (str) – If ‘rect’ or ‘exp’, compute atoms so that a transient (t0,tau) map can later be computed. (‘none’ instead of None explicitly calls the transient-window function, but with the full range, for debugging). See core.ComputeFstat. Currently only supported for nsegs=1.
tCWFstatMapVersion (str) – Choose between standard ‘lal’ implementation, ‘pycuda’ for gpu, and some others for devel/debug.
allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength].
- symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'transient_duration': '$\\Delta T$', 'transient_tstart': '$t_\\mathrm{start}$'}
Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots
- unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'transient_duration': 's', 'transient_tstart': 's'}
Key, val pairs of the parameters (F0, F1, …, including glitch parameters), and the units (Hz, Hz/s, …).
- transform_dictionary = {'transient_duration': {'multiplier': 1.1574074074074073e-05, 'symbol': 'Transient duration', 'unit': 'day'}, 'transient_tstart': {'label': 'Transient start-time \n days after minStartTime', 'multiplier': 1.1574074074074073e-05, 'subtractor': 'minStartTime', 'unit': 'day'}}
Key, val pairs of the parameters (F0, F1, …), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.
pyfstat.optimal_setup_functions module
Provides functions to aid in calculating the optimal setup for zoom follow up
- pyfstat.optimal_setup_functions.get_optimal_setup(NstarMax, Nsegs0, tref, minStartTime, maxStartTime, prior, detector_names)[source]
Using an optimisation step, calculate the optimal setup ladder
The details of the methods are described in Sec Va of arXiv:1802.05450. Here we provide implementation details. All equation numbers refer to arXiv:1802.05450.
- Parameters
NstarMax (float) – The ratio of the size at the old coherence time to the new coherence time for each step, see Eq. (31). Larger values allow a more rapid “zoom” of the search space at the cost of convergence. Smaller values are more conservative at the cost of additional computing time. The exact choice should be optimized for the problem in hand, but values of 100-1000 are typically okay.
Nsegs0 (int) – The number of segments for the initial step of the ladder
tref (int) – GPS times of the reference, start, and end time.
minStartTime (int) – GPS times of the reference, start, and end time.
maxStartTime (int) – GPS times of the reference, start, and end time.
prior (dict) – Prior dictionary, each item must either be a fixed scalar value, or a uniform prior.
detector_names (list or str) – Names of the detectors to use
- Returns
nsegs, Nstar – Ladder of segment numbers and the corresponding Nstar
- Return type
list
- pyfstat.optimal_setup_functions.get_Nstar_estimate(nsegs, tref, minStartTime, maxStartTime, prior, detector_names)[source]
Returns N* estimated from the super-sky metric
- Parameters
nsegs (int) – Number of semi-coherent segments
tref (int) – Reference time in GPS seconds
minStartTime (int) – Minimum and maximum SFT timestamps
maxStartTime (int) – Minimum and maximum SFT timestamps
prior (dict) – The prior dictionary
detector_names (array) – Array of detectors to average over
- Returns
Nstar – The estimated approximate number of templates to cover the prior parameter space at a mismatch of unity, assuming the normalised thickness is unity.
- Return type
int
pyfstat.snr module
- class pyfstat.snr.SignalToNoiseRatio(*, detector_states: lalpulsar.MultiDetectorStateSeries, noise_weights: Optional[lalpulsar.MultiNoiseWeights] = None, assumeSqrtSX: Optional[float] = None)[source]
Bases:
object
Compute the optimal SNR of a CW signal as expected in Gaussian noise.
The definition of SNR (shortcut for “optimal signal-to-noise ratio”) is taken from Eq. (76) of https://dcc.ligo.org/T0900149-v6/public and is such that \(\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2\), where \(\langle 2\mathcal{F}\rangle\) represents the expected value over noise realizations of twice the F-statistic of a template perfectly matched to an existing signal in the data.
Computing \(\textrm{SNR}^2\) requires two quantities:
- The antenna pattern matrix \(\mathcal{M}\), which depends on the sky position \(\vec{n}\)and polarization angle \(\psi\) and encodes the effect of the detector’s antenna pattern responseover the course of the observing run.
- The JKS amplitude parameters \((\mathcal{A}^0, \mathcal{A}^1, \mathcal{A}^2, \mathcal{A}^3)\)[JKS1998] which are functions of the CW’s amplitude parameters \((h_0,\cos\iota, \psi, \phi_0)\) or,alternatively, \((A_{+}, A_{\times}, \psi, \phi_0)\).
- Parameters
detector_states (lalpulsar.MultiDetectorStateSeries) – MultiDetectorStateSeries as produced by DetectorStates. Provides the required information to compute the antenna pattern contribution.
noise_weights (Union[lalpulsar.MultiNoiseWeights, None]) – Optional, incompatible with assumeSqrtSX. Can be computed from SFTs using SignalToNoiseRatio.from_sfts. Noise weights to account for a varying noise floor or unequal noise floors in different detectors.
assumeSqrtSX (float) – Optional, incompatible with noise_weights. Single-sided amplitude spectral density (ASD) of the detector noise. This value is used for all detectors, meaning it’s not currently possible to manually specify different noise floors without creating SFT files. (To be improved in the future; developer note: will require SWIG constructor for MultiNoiseWeights.)
Method generated by attrs for class SignalToNoiseRatio.
- detector_states: lalpulsar.MultiDetectorStateSeries
- noise_weights: Optional[lalpulsar.MultiNoiseWeights]
- assumeSqrtSX: float
- classmethod from_sfts(F0, sftfilepath, time_offset=None, running_median_window=101, sft_constraint=None)[source]
Alternative constructor to retrieve detector states and noise weights from SFT files. This method is based on
DetectorStates.multi_detector_states_from_sfts()
. This is currently the other way in which varying / different noise floors can be used when computing SNRs.- Parameters
F0 (float) – Central frequency [Hz] to retrieve from the SFT files to compute noise weights.
sftfilepath (str) – Path to SFT files in a format compatible with XLALSFTdataFind.
time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).
running_median_window (int) – Window used to compute the running-median noise floor estimation. Default value is consistent with that used in lalapps_PredictFstat.
sft_constraint (lalpulsar.SFTConstraint) – Optional argument to specify further constraints in XLALSFTdataFind.
- compute_snr2(Alpha, Delta, psi, phi, h0=None, cosi=None, aPlus=None, aCross=None)[source]
Compute the \(\textrm{SNR}^2\) of a CW signal using XLALComputeOptimalSNR2FromMmunu. Parameters correspond to the standard ones used to describe a CW (see e.g. Eqs. (16), (26), (30) of https://dcc.ligo.org/T0900149-v6/public ).
Mind that this function returns squared SNR (Eq. (76) of https://dcc.ligo.org/T0900149-v6/public ), which can be directly related to the expected F-statistic as \(\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2\).
- Parameters
Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.
Delta (float) – Declination (equatorial latitude) of the signal in radians.
psi (float) – Polarization angle.
h0 (float) – Nominal GW amplitude. Must be given together with cosi and conflicts with aPlus and aCross.
cosi (float) – Cosine of the source inclination w.r.t. line of sight. Must be given together with h0 and conflicts with aPlus and aCross.
aPlus (float) – Plus polarization amplitude. Must be given with aCross and conflicts with h0 and cosi.
aCross (float) – Cross polarization amplitude. Must be given with aPlus and conflicts with h0 and cosi.
- Returns
SNR^2 – Squared signal-to-noise ratio of a CW signal consistent with the specified parameters in the specified detector network.
- Return type
float
- compute_h0_from_snr2(Alpha, Delta, psi, phi, cosi, snr2)[source]
Convert the \(\textrm{SNR}^2\) of a CW signal to a corresponding amplitude \(h_0\) given the source orientation. Parameters correspond to the standard ones used to describe a CW (see e.g. Eqs. (16), (26), (30) of https://dcc.ligo.org/T0900149-v6/public ).
This function returns “inverts” Eq. (77) of https://dcc.ligo.org/T0900149-v6/public by computing the overall prefactor on \(h_0\) using self.compute_snr2(h0=1, …).
- Parameters
Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.
Delta (float) – Declination (equatorial latitude) of the signal in radians.
psi (float) – Polarization angle.
cosi (float) – Cosine of the source inclination w.r.t. line of sight. Must be given together with h0 and conflicts with aPlus and aCross.
snr2 (float) – Squared signal-to-noise ratio of a CW signal in the specified detector network.
- Returns
h0 – Nominal GW amplitude.
- Return type
float
- compute_twoF(*args, **kwargs)[source]
Compute the expected \(2\mathcal{F}\) value of a CW signal from the result of compute_snr2.
\[\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2\]\[\sigma_{2\mathcal{F}} = \sqrt{8 + 4 \textrm{SNR}^2}\]Input parameters are passed untouched to self.compute_snr2. See corresponding docstring for a list of valid parameters.
- Returns
expected_2F – Expected value of a non-central chi-squared distribution with four degrees of freedom and non-centrality parameter given by SNR^2.
stdev_2F – Standard deviation of a non-central chi-squared distribution with four degrees of freedom and non-centrality parameter given by SNR^2.
- compute_Mmunu(Alpha, Delta)[source]
Compute Mmunu matrix at a specific sky position using the detector states (and possible noise weights) given at initialization time. If no noise weigths were given, unit weights are assumed.
- Parameters
Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.
Delta (float) – Declination (equatorial latitude) of the signal in radians.
- Returns
Mmunu – Mmunu matrix encoding the response of the given detector network to a CW at the specified sky position.
- Return type
lalpulsar.AntennaPatternMatrix
- class pyfstat.snr.DetectorStates[source]
Bases:
object
Python interface to XLALGetMultiDetectorStates and XLALGetMultiDetectorStatesFromMultiSFTs.
- get_multi_detector_states(timestamps, Tsft, detectors=None, time_offset=None)[source]
- Parameters
timestamps (array-like or dict) – GPS timestamps at which detector states will be retrieved. If array, use the same set of timestamps for all detectors, which must be explicitly given by the user via detectors. If dictionary, each key should correspond to a valid detector name to be parsed by XLALParseMultiLALDetector and the associated value should be an array-like set of GPS timestamps for each individual detector.
Tsft (float) – Timespan covered by each timestamp. It does not need to coincide with the separation between consecutive timestamps.
detectors (list[str] or comma-separated string) – List of detectors to be parsed using XLALParseMultiLALDetector. Conflicts with dictionary of timestamps, required otherwise.
time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).
- Returns
multi_detector_states – Resulting multi-detector states produced by XLALGetMultiDetectorStates
- Return type
lalpulsar.MultiDetectorStateSeries
- get_multi_detector_states_from_sfts(sftfilepath, central_frequency, time_offset=None, frequency_wing_bins=1, sft_constraint=None, return_sfts=False)[source]
- Parameters
sftfilepath (str) – Path to SFT files in a format compatible with XLALSFTdataFind.
central_frequency (float) – Frequency [Hz] around which SFT data will be retrieved. This option is only relevant if further information is to be retrieved from the SFTs (i.e. return_sfts=True).
time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).
frequency_wing_bins (int) – Frequency bins around the central frequency to retrieve from SFT data. Bin size is determined using the SFT baseline time as obtained from the catalog. This option is only relevant if further information is to be retrieved from the SFTs (i.e. return_sfts=True).
sft_constraint (lalpulsar.SFTConstraint) – Optional argument to specify further constraints in XLALSFTdataFind.
return_sfts (bool) – If True, also return the loaded SFTs. This is useful to compute further quantities such as noise weights.
- Returns
multi_detector_states (lalpulsar.MultiDetectorStateSeries) – Resulting multi-detector states produced by XLALGetMultiDetectorStatesFromMultiSFTs
multi_sfts (lalpulsar.MultiSFTVector) – Only if return_sfts is True. MultiSFTVector produced by XLALLoadMultiSFTs along the specified frequency band.
pyfstat.tcw_fstat_map_funcs module
Additional helper functions dealing with transient-CW F(t0,tau) maps.
See Prix, Giampanis & Messenger (PRD 84, 023007, 2011): https://arxiv.org/abs/1104.1704 for the algorithm in general and Keitel & Ashton (CQG 35, 205003, 2018): https://arxiv.org/abs/1805.05652 for a detailed discussion of the GPU implementation.
- class pyfstat.tcw_fstat_map_funcs.pyTransientFstatMap(N_t0Range=None, N_tauRange=None, transientFstatMap_t=None, from_file=None)[source]
Bases:
object
Simplified object class for a F(t0,tau) F-stat map.
This is based on LALSuite’s transientFstatMap_t type, replacing the gsl matrix with a numpy array.
Here, t0 is a transient start time, tau is a transient duration parameter, and F(t0,tau) is the F-statistic (not 2F)! evaluated for a signal with those parameters (and an implicit window function, which is not stored inside this object).
The ‘map’ covers a range of different (t0,tau) pairs.
- F_mn
2D array of F values (not 2F!), with m an index over start-times t0, and n an index over duration parameters tau, in steps of dt0 in [t0, t0+t0Band], and dtau in [tau, tau+tauBand].
- Type
np.ndarray
- maxF
Maximum of F (not 2F!) over the array.
- Type
float
- t0_ML
Maximum likelihood estimate of the transient start time t0.
- Type
float
- tau_ML
Maximum likelihood estimate of the transient duration tau.
- Type
float
- lnBtSG
Natural log of the marginalised transient Bayes factor. NOTE: This is always initialised as nan, and you have to call get_lnBtSG() to get its actual value.
- Type
float
The class can be initialized with either a pair of (N_t0Range,N_tauRange), from a lalpulsar object, or reading from a file.
- Parameters
N_t0Range (int) – Number of t0 values covered.
N_tauRange (int) – Number of tau values covered.
transientFstatMap_t (lalpulsar.transientFstatMap_t) – pre-allocated matrix from lalpulsar to initialise from.
from_file (str or None) – Text file, compatible with lalpulsar.write_transientFstatMap_to_fp() format, to load and initialise from.
- read_from_file(file)[source]
Read F_mn map from a text file and set all other fields.
Apart from optional header lines (# comments), the format has to be consistent with lalpulsar.write_transientFstatMap_to_fp() and the write_F_mn_to_file() method of this class itself: with the columns [t0[s], tau[s], 2F]. NOTE that the file is expected to provide 2F, so the values will be halved to obtain F for storage in this class.
- Parameters
file (str) – Name of the file to load from.
- get_maxF_idx()[source]
Gets the 2D-unravelled index pair of the maximum of the F_mn map
- Returns
idx – The m,n indices of the map entry with maximal F value.
- Return type
tuple
- get_lnBtSG()[source]
Compute (natural log of the) transient-CW Bayes-factor B_tSG = P(x|HyptS)/P(x|HypG).
Here HypG = Gaussian noise hypothesis, HyptS = transient-CW signal hypothesis.
B_tSG is marginalized over start-time and timescale of transient CW signal, using given type and parameters of transient window range.
This is a python port of the lalpulsar.ComputeTransientBstat implementation, replacing for loops by numpy operations.
- write_F_mn_to_file(tCWfile, windowRange, header=[])[source]
Format a 2D transient-F-stat matrix over (t0,tau) and write as a text file.
Apart from the optional extra header lines, the format is consistent with lalpulsar.write_transientFstatMap_to_fp(), with the columns [t0[s], tau[s], 2F]. NOTE that the output is 2F, not F like stored in this class itself!
- Parameters
tCWfile (str) – Name of the file to write to.
windowRange (lalpulsar.transientWindowRange_t) – A lalpulsar structure containing the transient parameters.
header (list) – A list of additional header lines to print at the start of the file.
- pyfstat.tcw_fstat_map_funcs.fstatmap_versions = {'lal': <function <lambda>>, 'pycuda': <function <lambda>>}
Dictionary of the actual callable transient F-stat map functions this module supports.
Actual runtime availability depends on the corresponding external modules being available.
- pyfstat.tcw_fstat_map_funcs.init_transient_fstat_map_features(feature='lal', cudaDeviceName=None)[source]
Initialization of available modules (or ‘features’) for computing transient F-stat maps.
Currently, two implementations are supported and checked for through the _optional_import() method:
lal: requires both lal and lalpulsar packages to be importable.
2. pycuda: requires the pycuda package to be importable along with its modules driver, gpuarray, tools and compiler.
- Parameters
feature (str) – Set the transient F-stat map implementation.
cudaDeviceName (str or None) – Request a CUDA device with this name. Partial matches are allowed.
- Returns
features (dict) – A dictionary of available method names, to match fstatmap_versions. Each key’s value is set to True only if all required modules are importable on this system.
gpu_context (pycuda.driver.Context or None) – A CUDA device context object, if assigned.
- pyfstat.tcw_fstat_map_funcs.call_compute_transient_fstat_map(version, features, multiFstatAtoms=None, windowRange=None, BtSG=False)[source]
Call a version of the ComputeTransientFstatMap function.
This checks if the requested version is available, and if so, executes the computation of a transient F-statistic map over the windowRange.
- Parameters
version (str) – Name of the method to call (currently supported: ‘lal’ or ‘pycuda’).
features (dict) – Dictionary of available features, as obtained from init_transient_fstat_map_features().
multiFstatAtoms (lalpulsar.MultiFstatAtomVector or None) – The time-dependent F-stat atoms previously computed by ComputeFstat.
windowRange (lalpulsar.transientWindowRange_t or None) – The structure defining the transient parameters.
BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using the appropriate implementation for each feature, and store it in FstatMap.lnBtSG.
- Returns
FstatMap (pyTransientFstatMap or lalpulsar.transientFstatMap_t) – The output of the called function, including the evaluated transient F-statistic map over the windowRange.
timingFstatMap (float) – Execution time of the called function.
- pyfstat.tcw_fstat_map_funcs.lalpulsar_compute_transient_fstat_map(multiFstatAtoms, windowRange, BtSG=False)[source]
Wrapper for the standard lalpulsar function for computing a transient F-statistic map.
See https://lscsoft.docs.ligo.org/lalsuite/lalpulsar/_transient_c_w__utils_8h.html for the wrapped function.
- Parameters
multiFstatAtoms (lalpulsar.MultiFstatAtomVector) – The time-dependent F-stat atoms previously computed by ComputeFstat.
windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.
BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using the corresponding lalpulsar function, and store it in FstatMap.lnBtSG.
- Returns
FstatMap – The computed results, see the class definition for details.
- Return type
- pyfstat.tcw_fstat_map_funcs.reshape_FstatAtomsVector(atomsVector)[source]
Make a dictionary of ndarrays out of an F-stat atoms ‘vector’ structure.
- Parameters
atomsVector (lalpulsar.FstatAtomVector) – The atoms in a ‘vector’-like structure: iterating over timestamps as the higher hierarchical level, with a set of ‘atoms’ quantities defined at each timestamp.
- Returns
atomsDict – A dictionary with an entry for each quantity, which then is a 1D ndarray over timestamps for that one quantity.
- Return type
dict
- pyfstat.tcw_fstat_map_funcs.pycuda_compute_transient_fstat_map(multiFstatAtoms, windowRange, BtSG=False)[source]
GPU version of computing a transient F-statistic map.
This is based on XLALComputeTransientFstatMap from LALSuite, (C) 2009 Reinhard Prix, licensed under GPL.
The ‘map’ consists of F-statistics evaluated over a range of different (t0,tau) pairs (transient start-times and duration parameters).
This is a high-level wrapper function; the actual CUDA computations are performed in one of the functions pycuda_compute_transient_fstat_map_rect() or pycuda_compute_transient_fstat_map_exp(), depending on the window functon defined in windowRange.
- Parameters
multiFstatAtoms (lalpulsar.MultiFstatAtomVector) – The time-dependent F-stat atoms previously computed by ComputeFstat.
windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.
BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using a CPU python port of the corresponding lalpulsar function, and store it in FstatMap.lnBtSG.
- Returns
FstatMap – The computed results, see the class definition for details.
- Return type
- pyfstat.tcw_fstat_map_funcs.pycuda_compute_transient_fstat_map_rect(atomsInputMatrix, windowRange, tCWparams)[source]
GPU computation of the transient F-stat map for rectangular windows.
As discussed in Keitel & Ashton (CQG 35, 205003, 2018): https://arxiv.org/abs/1805.05652 this version only does GPU parallization for the outer loop, keeping the partial sums of the inner loop local to each individual kernel using the ‘memory trick’.
- Parameters
atomsInputMatrix (np.ndarray) – A 2D array of stacked named columns containing the F-stat atoms.
windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.
tCWparams (dict) – A dictionary of miscellaneous parameters.
- Returns
F_mn – A 2D array of the computed transient F-stat map over the [t0,tau] range.
- Return type
np.ndarray
- pyfstat.tcw_fstat_map_funcs.pycuda_compute_transient_fstat_map_exp(atomsInputMatrix, windowRange, tCWparams)[source]
GPU computation of the transient F-stat map for exponential windows.
As discussed in Keitel & Ashton (CQG 35, 205003, 2018): https://arxiv.org/abs/1805.05652 this version does full GPU parallization of both the inner and outer loop.
- Parameters
atomsInputMatrix (np.ndarray) – A 2D array of stacked named columns containing the F-stat atoms.
windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.
tCWparams (dict) – A dictionary of miscellaneous parameters.
- Returns
F_mn – A 2D array of the computed transient F-stat map over the [t0,tau] range.
- Return type
np.ndarray