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 (looking in ~/.pyfstat or $LALPULSAR_DATADIR)

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.

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.

  • 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 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.

  • 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.

  • 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().

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 (twoF or log10BSGL) fully-coherently at a single point.

These are also stored to self.twoF and self.log10BSGL respectively. As the basic statistic of this class, self.twoF is always computed. If self.BSGL, additionally the single-detector 2F-stat values are saved in self.twoFX.

If transient parameters are enabled (self.transientWindowType is set), the full transient-F-stat map will also be 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_detstat(), 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_detstat(), i.e. only used if self.transientWindowType is set.

Returns

stat – A single value of the detection statistic (twoF or log10BSGL) at the input parameter values. Also stored as self.twoF or self.log10BSGL.

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_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.

The full transient-F-stat map will also be 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_detstat(), 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_detstat(), i.e. only used if self.transientWindowType is set.

Returns

maxTwoF – A single value of the detection statistic (twoF or log10BSGL) 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 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, 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.

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.BSGL, additionally the single-detector 2F-stat values are saved in self.twoFX.

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.BSGL=True) 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

compute_glitch_fstat_single(F0, F1, F2, Alpha, Delta, delta_F0, delta_F1, tglitch)[source]

Returns the semi-coherent glitch summed twoF for nglitch=1.

NOTE: OBSOLETE, used only for testing.

class pyfstat.core.DeprecatedClass(*args, **kwargs)[source]

Bases: object

Outdated classes are marked for future removal by inheriting from this.

class pyfstat.core.DefunctClass(*args, **kwargs)[source]

Bases: object

Removed classes are retained for a while but marked by inheriting from this.

last_supported_version = None
pr_welcome = True

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}$', 'log10BSGL': '$\\log_{10}\\mathcal{B}_{\\mathrm{SGL}}$', '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:

  1. Output file with matching name found in outdir.

  2. Output file is not older than SFT files matching sftfilepattern.

  3. Parameters string in file header matches current search setup.

  4. 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

  1. none (standard, persistent CW signal)

  2. rect (rectangular: constant amplitude within [t0,t0+tau], zero outside)

  3. 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.

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:

  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

  2. 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.max_slice(D, axis)[source]

Return the slice along the given axis

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.gridcorner.projection_2D(ax, x, y, D, xidx, yidx, projection, lines=None, **kwargs)[source]
pyfstat.gridcorner.projection_1D(ax, x, D, xidx, projection, showDvals=True, lines=None, **kwargs)[source]

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 and next in the $LALPULSAR_DATADIR environment variable.

If neither source provides the necessary files, a warning is emitted and the user can still run PyFstat searches, but must then include ephemeris options manually on each class instantiation.

The ‘DE405’ ephemerides version provided with lalpulsar is expected.

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_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)[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.

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).

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.make_sfts module

PyFstat tools to generate and manipulate data in the form of SFTs.

exception pyfstat.make_sfts.KeyboardInterruptError[source]

Bases: Exception

Custom exception to allow overriding interrupts.

class pyfstat.make_sfts.InjectionParametersGenerator(priors=None, seed=None)[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 – Argument to be fed to numpy.random.default_rng, with all of its accepted types.

set_priors(new_priors)[source]

Set priors to draw parameter space points from.

Parameters

new_priors (dict) – The new set of priors to update the object with.

set_seed(seed)[source]

Set the random seed for subsequent draws.

Parameters

seed – Argument to be fed to numpy.random.default_rng, with all of its accepted types.

draw()[source]

Draw a single multi-dimensional parameter space point from the given priors.

Returns

injection_parameters – Dictionary with parameter names as keys and their numeric values.

Return type

dict

class pyfstat.make_sfts.AllSkyInjectionParametersGenerator(priors=None, seed=None)[source]

Bases: pyfstat.make_sfts.InjectionParametersGenerator

Like InjectionParametersGenerator, but with hardcoded all-sky priors.

This ensures uniform coverage of the 2D celestial sphere: uniform distribution in Alpha and sine distribution for Delta.

It assumes 1) PyFstat notation and 2) equatorial coordinates.

Alpha and Delta are given ‘restricted’ status to stop the user from changing them as long as using this special class.

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 – Argument to be fed to numpy.random.default_rng, with all of its accepted types.

set_priors(new_priors)[source]

Set priors to draw parameter space points from.

Parameters

new_priors (dict) – The new set of priors to update the object with.

set_seed(seed)[source]

Set the random seed for subsequent draws.

Parameters

seed – Argument to be fed to numpy.random.default_rng, with all of its accepted types.

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.

  • duration (int) – Duration (in GPS seconds) of the total data set.

  • 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.

  • 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.

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.

tend()[source]

Simple method to return tstart+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, 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 you want to use noiseSFTs but auto-estimate a minimal band, call helper_functions.get_covering_band() yourself and pass the results to Writer as fmin, Band.

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.

make_data(verbose=False)[source]

A convenience wrapper to generate a cff file and then SFTs.

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 and cosi will be ignored. cosi is used to rescale h0 with the usual sqrt(cosi**4 + 6*cosi**2 + 1) relation.

Parameters
  • label (string) – A human-readable label to be used in naming the output files.

  • tstart (int) – Starting GPS epoch of the data set.

  • duration (int) – Duration (in GPS seconds) of the total data set.

  • 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.

  • 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.

mfd = 'lalapps_Makefakedata_v4'

The executable (older version that supports the –lineFeature option).

required_mfd_line_parameters = ['Freq', 'phi0', 'h0', 'cosi', 'transientWindowType', 'transientStartTime', 'transientTau']

Other signal parameters will be removed before passing to MFDv4.

calculate_fmin_Band()[source]

Set fmin and Band for the output SFTs to cover.

This method adapts Writer.calculate_fmin_Band to work with MakeFakeData_v4, as MakeFakeData_v4 requires fmin and Band to be given.

The overriding prevents (fmin, Band) from not being set when self.noiseSFTs evaluates to true.

See Writer.calculate_fmin_Band for further details.

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.

make_cff(verbose=False)[source]

Generates a .cff file including signal injection parameters, including a glitch.

This will be saved to self.config_file_name.

Parameters

verbose (boolean) – If true, increase logging verbosity.

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

concatenate_sft_files()[source]

Merges the individual SFT files via lalapps_splitSFTs executable.

pre_compute_evolution()[source]

Precomputes evolution parameters for the artifact.

This computes midtimes, frequencies, phases and amplitudes over the list of SFT timestamps.

make_ith_sft(i)[source]

Call MFDv4 to create a single SFT with evolved artifact parameters.

make_data()[source]

Create a full multi-SFT data set.

This loops over SFTs and generate them serially or in parallel, then contatenates the results together at the end.

run_makefakedata_v4(mid_time, lineFreq, linePhi, h0, tmp_outdir)[source]

Generate SFT data using the MFDv4 code with the –lineFeature option.

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.

get_h0(t)[source]

Evaluate the artifact amplitude at a given time.

NOTE: Here it’s actually changing over time!

Parameters

t (float) – Time stamp to evaluate at.

Returns

h0 – Amplitude at time t.

Return type

float

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.

  • 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.

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.

Keep in mind that one has to explicitely set savefig=True to output the figure!

See pyfstat.core.ComputeFstat.plot_twoF_cumulative for a comprehensive list of accepted arguments and their default values.

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

export_samples_to_disk()[source]

Export MCMC samples into a text file using numpy.savetxt.

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”.

generate_loudest()[source]

Use lalapps_ComputeFstatistic_v2 to produce a .loudest file

write_prior_table()[source]

Generate a .tex file of the prior

print_summary()[source]

Prints a summary of the max twoF found to the terminal

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

write_evidence_file_from_dict(EvidenceDict, evidence_file_name)[source]

Write evidence dict to a file

Parameters
  • EvidenceDict (dict) – Dictionary to dump into a file.

  • evidence_file_name (str) – File name to dump dict into.

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=False)[source]

Override MCMCSearch.plot_cumulative_max implementation to dea with the split at gltiches.

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

read_setup_input_file(run_setup_input_file)[source]
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.

  • 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.

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.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)[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

The class can be initialized with the following:

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.

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

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().

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(wantCuda=False, 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:

  1. 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
  • wantCuda (bool) – Only if this is True and it was possible to import pycuda, a CUDA device context is created and returned.

  • 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)[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.

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)[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.

Returns

FstatMap – The computed results, see the class definition for details.

Return type

pyTransientFstatMap

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)[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.

Returns

FstatMap – The computed results, see the class definition for details.

Return type

pyTransientFstatMap

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

Module contents