# pyfstat package¶

These pages document the full API and set of classes provided by PyFstat.

See here for installation instructions and other general information.

## 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 utils.get_ephemeris_files() are used.

Parameters:
• earth_ephem (str) – Paths of the two files containing positions of Earth and Sun, respectively at evenly spaced times, as passed to CreateFstatInput

• sun_ephem (str) – Paths of the two files containing positions of Earth and Sun, respectively at evenly spaced times, as passed to CreateFstatInput

pprint_init_params_dict()[source]

Pretty-print a parameters dictionary for output file headers.

Returns:

pretty_init_parameters – A list of lines to be printed, including opening/closing “{” and “}”, consistent indentation, as well as end-of-line commas, but no comment markers at start of lines.

Return type:

list

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:

Return type:

list

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]

Base search class providing an interface to lalpulsar.ComputeFstat.

In most cases, users should be using one of the higher-level search classes from the grid_based_searches or mcmc_based_searches modules instead.

See the lalpulsar documentation at https://lscsoft.docs.ligo.org/lalsuite/lalpulsar/group___compute_fstat__h.html and R. Prix, The F-statistic and its implementation in ComputeFstatistic_v2 ( https://dcc.ligo.org/T0900149/public ) for details of the lalpulsar module and the meaning of various technical concepts as embodied by some of the class’s parameters.

Normally this will read in existing data through the sftfilepattern argument, but if that option is None and the necessary alternative arguments are used, it can also generate simulated data (including noise and/or signals) on the fly.

NOTE that the detection statistics that can be computed from an instance of this class depend on the BSGL, BtSG and transientWindowType arguments given at initialisation. See get_fullycoherent_detstat() and get_transient_detstats() for details. To change what you want to compute, you may need to initialise a new instance with different options.

NOTE for GPU users (tCWFstatMapVersion=”pycuda”): This class tries to conveniently deal with GPU context management behind the scenes. A known problematic case is if you try to instantiate it twice from the same session/script. If you then get some messages like RuntimeError: make_default_context() and invalid device context, that is because the GPU is still blocked from the first instance when you try to initiate the second. To avoid this problem, use context management:

with pyfstat.ComputeFstat(
[...],
tCWFstatMapVersion="pycuda",
) as search:
search.get_fullycoherent_detstat([...])


or manually call the search.finalizer_() method where needed.

Parameters:
• tref (int) – GPS seconds of the reference time.

• sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.

• minStartTime (int) – Only use SFTs with timestamps starting from within this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime].

• maxStartTime (int) – Only use SFTs with timestamps starting from within this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime].

• Tsft (int) – SFT duration in seconds. Only required if sftfilepattern=None and hence simulted data is generated on the fly.

• binary (bool) – If true, search over binary parameters.

• singleFstats (bool) – If true, also compute the single-detector twoF values.

• BSGL (bool) –

If true, compute the log10BSGL statistic rather than the twoF value. For details, see Keitel et al (PRD 89, 064023, 2014): https://arxiv.org/abs/1311.5738 Note this automatically sets singleFstats=True as well. Tuning parameters are currently hardcoded:

• Fstar0=15 for coherent searches.

• A p-value of 1e-6 and correspondingly recalculated Fstar0 for semicoherent searches.

• Uniform per-detector prior line-vs-Gaussian odds.

• BtSG (bool) – If true and transientWindowType is not None, compute the transient $$\ln\mathcal{B}_{\mathrm{tS}/\mathrm{G}}$$ statistic from Prix, Giampanis & Messenger (PRD 84, 023007, 2011) (tCWFstatMap marginalised over uniform t0, tau priors). rather than the maxTwoF value.

• transientWindowType (str) – If rect or exp, allow for the Fstat to be computed over a transient range. (none instead of None explicitly calls the transient-window function, but with the full range, for debugging.) (If not None, will also force atoms regardless of computeAtoms option.)

• t0Band (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave Fstat with t0=minStartTime, tau=maxStartTime-minStartTime.

• tauBand (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave Fstat with t0=minStartTime, tau=maxStartTime-minStartTime.

• tauMin (int) – Minimum transient duration to cover, defaults to 2*Tsft.

• dt0 (int) – Grid resolution in transient start-time, defaults to Tsft.

• dtau (int) – Grid resolution in transient duration, defaults to Tsft.

• detectors (str) – Two-character references to the detectors for which to use data. Specify None for no constraint. For multiple detectors, separate by commas.

• minCoverFreq (float) – The min and max cover frequency passed to lalpulsar.CreateFstatInput. For negative values, these will be used as offsets from the min/max frequency contained in the sftfilepattern. If either is None, the search_ranges argument is used to estimate them. If the automatic estimation fails and you do not have a good idea what to set these two options to, setting both to -0.5 will reproduce the default behaviour of PyFstat <=1.4 and may be a reasonably safe fallback in many cases.

• maxCoverFreq (float) – The min and max cover frequency passed to lalpulsar.CreateFstatInput. For negative values, these will be used as offsets from the min/max frequency contained in the sftfilepattern. If either is None, the search_ranges argument is used to estimate them. If the automatic estimation fails and you do not have a good idea what to set these two options to, setting both to -0.5 will reproduce the default behaviour of PyFstat <=1.4 and may be a reasonably safe fallback in many cases.

• search_ranges (dict) – Dictionary of ranges in all search parameters, only used to estimate frequency band passed to lalpulsar.CreateFstatInput, if minCoverFreq, maxCoverFreq are not specified (==None). For actually running searches, grids/points will have to be passed separately to the .run() method. The entry for each parameter must be a list of length 1, 2 or 3: [single_value], [min,max] or [min,max,step].

• injectSources (dict or str) – Either a dictionary of the signal parameters to inject, or a string pointing to a .cff file defining a signal.

• injectSqrtSX (float or list or str) – Single-sided PSD values for generating fake Gaussian noise on the fly. Single float or str value: use same for all IFOs. List or comma-separated string: must match len(detectors) and/or the data in sftfilepattern. Detectors will be paired to list elements following alphabetical order.

• randSeed (int or None) – random seed for on-the-fly noise generation using injectSqrtSX. Setting this to 0 or None is equivalent; both will randomise the seed, following the behaviour of XLALAddGaussianNoise(), while any number not equal to 0 will produce a reproducible noise realisation.

• assumeSqrtSX (float or list or str) – Don’t estimate noise-floors but assume this (stationary) single-sided PSD. Single float or str value: use same for all IFOs. List or comma-separated string: must match len(detectors) and/or the data in sftfilepattern. Detectors will be paired to list elements following alphabetical order. If working with signal-only data, please set assumeSqrtSX=1 .

• SSBprec (int) – Flag to set the Solar System Barycentring (SSB) calculation in lalpulsar: 0=Newtonian, 1=relativistic, 2=relativistic optimised, 3=DMoff, 4=NO_SPIN

• RngMedWindow (int) – Running-Median window size for F-statistic noise normalization (number of SFT bins).

• tCWFstatMapVersion (str) – Choose between implementations of the transient F-statistic funcionality: standard lal implementation, pycuda for GPU version, and some others only for devel/debug.

• cudaDeviceName (str) – GPU name to be matched against drv.Device output, only for tCWFstatMapVersion=pycuda.

• computeAtoms (bool) – Request calculation of ‘F-statistic atoms’ regardless of transientWindowType.

• earth_ephem (str) – Earth ephemeris file path. If None, will check standard sources as per utils.get_ephemeris_files().

• sun_ephem (str) – Sun ephemeris file path. If None, will check standard sources as per utils.get_ephemeris_files().

• allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength]

init_computefstatistic()[source]

Initialization step for the F-stastic computation internals.

This sets up the special input and output structures the lalpulsar module needs, the ephemerides, optional on-the-fly signal injections, and extra options for multi-detector consistency checks and transient searches.

All inputs are taken from the pre-initialized object, so this function does not have additional arguments of its own.

estimate_min_max_CoverFreq()[source]

Extract spanned spin-range at reference -time from the template bank.

To use this method, self.search_ranges must be a dictionary of lists per search parameter which can be either [single_value], [min,max] or [min,max,step].

get_fullycoherent_detstat(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, tstart=None, tend=None)[source]

Computes the detection statistic(s) fully-coherently at a single point.

Currently supported statistics:

• twoF (CW)

• log10BSGL (CW or transient)

• maxTwoF (transient)

• lnBtSG (transient)

All computed statistics are stored as attributes, but only one statistic is returned.

As the basic statistic of this class, twoF is always computed and stored as self.twoF as well, and it is the default return value.

If self.singleFstats, additionally the single-detector 2F-stat values are stored in self.twoFX.

If self.BSGL, the log10BSGL statistic for CWs is additionally stored, and it is returned instead of twoF.

If transient parameters are enabled (self.transientWindowType is set), maxTwoF will always be computed and stored, and returned by default. Depending on the self.BSGL and self.BtSG options, either log10BSGL (a transient version of it, superseding the CW version) or lnBtSG will also be computed, stored, and returned instead of maxTwoF. The full transient-F-stat map is also computed here, but stored in self.FstatMap, not returned.

Parameters:
• F0 (float) – Parameters at which to compute the statistic.

• F1 (float) – Parameters at which to compute the statistic.

• F2 (float) – Parameters at which to compute the statistic.

• Alpha (float) – Parameters at which to compute the statistic.

• Delta (float) – Parameters at which to compute the statistic.

• asini (float, optional) – Optional: Binary parameters at which to compute the statistic.

• period (float, optional) – Optional: Binary parameters at which to compute the statistic.

• ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.

• tp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• argp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. This is only passed on to self.get_transient_detstats(), i.e. only used if self.transientWindowType is set.

• tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. This is only passed on to self.get_transient_detstats(), i.e. only used if self.transientWindowType is set.

Returns:

stat – A single value of the main detection statistic at the input parameter values.

Return type:

float

get_fullycoherent_twoF(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None)[source]

Computes the fully-coherent 2F statistic at a single point.

NOTE: This always uses the full data set as defined when initialising the search object. If you want to restrict the range of data used for a single 2F computation, you need to set a self.transientWindowType and then call self.get_fullycoherent_detstat() with tstart and tend options instead of this funcion.

Parameters:
• F0 (float) – Parameters at which to compute the statistic.

• F1 (float) – Parameters at which to compute the statistic.

• F2 (float) – Parameters at which to compute the statistic.

• Alpha (float) – Parameters at which to compute the statistic.

• Delta (float) – Parameters at which to compute the statistic.

• asini (float, optional) – Optional: Binary parameters at which to compute the statistic.

• period (float, optional) – Optional: Binary parameters at which to compute the statistic.

• ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.

• tp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• argp (float, optional) – Optional: Binary parameters at which to compute the statistic.

Returns:

twoF – A single value of the fully-coherent 2F statistic at the input parameter values. Also stored as self.twoF.

Return type:

float

get_fullycoherent_single_IFO_twoFs()[source]

Computes single-detector F-stats at a single point.

This requires self.get_fullycoherent_twoF() to be run first.

Returns:

twoFX – A list of the single-detector detection statistics twoF. Also stored as self.twoFX.

Return type:

list

get_fullycoherent_log10BSGL()[source]

Computes the line-robust statistic log10BSGL at a single point.

This requires self.get_fullycoherent_twoF() and self.get_fullycoherent_single_IFO_twoFs() to be run first.

Returns:

log10BSGL – A single value of the detection statistic log10BSGL at the input parameter values. Also stored as self.log10BSGL.

Return type:

float

get_transient_detstats(tstart=None, tend=None)[source]

Computes one or more transient detection statistics at a single point.

This requires self.get_fullycoherent_twoF() to be run first.

All computed statistics will be stored as attributes of self, but only one (twoF, log10BSGL or lnBtSG) will be the return value.

The full transient-F-stat map will also be computed here, but stored in self.FstatMap, not returned.

Parameters:
• tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.

• tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.

Returns:

detstat – A single value of the main chosen detection statistic (maxTwoF, log10BSGL or lnBtSG) at the input parameter values.

Return type:

float

get_transient_maxTwoFstat(tstart=None, tend=None)[source]

Computes the transient maxTwoF statistic at a single point.

This requires self.get_fullycoherent_twoF() to be run first, and is itself now only a backwards compatibility / convenience wrapper around the more general get_transient_detstats().

The full transient-F-stat map will also be computed here, but stored in self.FstatMap, not returned.

Parameters:
• tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.

• tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime.

Returns:

maxTwoF – A single value of the detection statistic at the input parameter values. Also stored as self.maxTwoF.

Return type:

float

get_transient_log10BSGL()[source]

Computes a transient detection statistic log10BSGL at a single point.

This should normally be called through get_transient_detstats(), but if called stand-alone, it requires self.get_transient_maxTwoFstat() to be run first.

The single-detector 2F-stat values used for that computation (at the index of maxTwoF) are saved in self.twoFXatMaxTwoF, not returned.

Returns:

log10BSGL – A single value of the detection statistic log10BSGL at the input parameter values. Also stored as self.log10BSGL.

Return type:

float

calculate_twoF_cumulative(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, tstart=None, tend=None, transient_tstart=None, transient_duration=None, num_segments=1000)[source]

Calculate the cumulative twoF over subsets of the observation span.

This means that we consider sub-“segments” of the [tstart,tend] interval, each starting at the overall tstart and with increasing durations, and compute the 2F for each of these, which for a true CW signal should increase roughly with duration towards the full value.

Parameters:
• F0 (float) – Parameters at which to compute the cumulative twoF.

• F1 (float) – Parameters at which to compute the cumulative twoF.

• F2 (float) – Parameters at which to compute the cumulative twoF.

• Alpha (float) – Parameters at which to compute the cumulative twoF.

• Delta (float) – Parameters at which to compute the cumulative twoF.

• asini (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.

• period (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.

• ecc (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.

• tp (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.

• argp (float, optional) – Optional: Binary parameters at which to compute the cumulative 2F.

• tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime;. If outside those: auto-truncated.

• tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime;. If outside those: auto-truncated.

• num_segments (int) – Number of segments to split [tstart,tend] into.

• transient_tstart (float or None) – These are not actually used by this function, but just included so a parameters dict can be safely passed.

• transient_duration (float or None) – These are not actually used by this function, but just included so a parameters dict can be safely passed.

Returns:

• cumulative_durations (ndarray of shape (num_segments,)) – Offsets of each segment’s tend from the overall tstart.

• twoFs (ndarray of shape (num_segments,)) – Values of twoF computed over [[tstart,tstart+duration] for duration in cumulative_durations].

predict_twoF_cumulative(F0, Alpha, Delta, h0, cosi, psi, tstart=None, tend=None, num_segments=10, **predict_fstat_kwargs)[source]

Calculate expected 2F, with uncertainty, over subsets of the observation span.

This yields the expected behaviour that calculate_twoF_cumulative() can be compared against: 2F for CW signals increases with duration as we take longer and longer subsets of the total observation span.

Parameters:
• F0 (float) – Parameters at which to compute the cumulative predicted twoF.

• Alpha (float) – Parameters at which to compute the cumulative predicted twoF.

• Delta (float) – Parameters at which to compute the cumulative predicted twoF.

• h0 (float) – Parameters at which to compute the cumulative predicted twoF.

• cosi (float) – Parameters at which to compute the cumulative predicted twoF.

• psi (float) – Parameters at which to compute the cumulative predicted twoF.

• tstart (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.

• tend (int or None) – GPS times to restrict the range of data used. If None: falls back to self.minStartTime and self.maxStartTime. If outside those: auto-truncated.

• num_segments (int) – Number of segments to split [tstart,tend] into.

• predict_fstat_kwargs – Other kwargs to be passed to utils.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: ComputeFstat

A simple semi-coherent search class.

This will split the data set into multiple segments, run a coherent F-stat search over each, and produce a final semi-coherent detection statistic as the sum over segments.

This does not include any concept of refinement between the two steps, as some grid-based semi-coherent search algorithms do; both the per-segment coherent F-statistics and the incoherent sum are done at the same parameter space point.

The implementation is based on a simple trick using the transient F-stat map functionality: basic F-stat atoms are computed only once over the full data set, then the transient code with rectangular ‘windows’ is used to compute the per-segment F-stats, and these are summed to get the semi-coherent result.

Only parameters with a special meaning for SemiCoherentSearch itself are explicitly documented here. For all other parameters inherited from pyfstat.ComputeFStat see the documentation of that class.

Parameters:
• label (str) – A label and directory to read/write data from/to.

• outdir (str) – A label and directory to read/write data from/to.

• tref (int) – GPS seconds of the reference time.

• nsegs (int) – The (fixed) number of segments to split the data set into.

• sftfilepattern (str) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; multiple patterns can be given separated by colons.

• minStartTime (int) – Only use SFTs with timestamps starting from this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime]. Also used to set up segment boundaries, i.e. maxStartTime-minStartTime will be divided by nsegs to obtain the per-segment coherence time Tcoh.

• maxStartTime (int) – Only use SFTs with timestamps starting from this range, following the XLALCWGPSinRange convention: half-open intervals [minStartTime,maxStartTime]. Also used to set up segment boundaries, i.e. maxStartTime-minStartTime will be divided by nsegs to obtain the per-segment coherence time Tcoh.

init_semicoherent_parameters()[source]

Set up a list of equal-length segments and the corresponding transient windows.

For a requested number of segments self.nsegs, self.tboundaries will have self.nsegs+1 entries covering [self.minStartTime,self.maxStartTime] and self.Tcoh will be the total duration divided by self.nsegs.

Each segment is required to be at least two SFTs long.f

get_semicoherent_det_stat(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, record_segments=False)[source]

Computes the detection statistic (twoF or log10BSGL) semi-coherently at a single point.

As the basic statistic of this class, self.twoF is always computed. If self.singleFstats, additionally the single-detector 2F-stat values are saved in self.twoFX and (optionally) self.twoFX_per_segment.

Parameters:
• F0 (float) – Parameters at which to compute the statistic.

• F1 (float) – Parameters at which to compute the statistic.

• F2 (float) – Parameters at which to compute the statistic.

• Alpha (float) – Parameters at which to compute the statistic.

• Delta (float) – Parameters at which to compute the statistic.

• asini (float, optional) – Optional: Binary parameters at which to compute the statistic.

• period (float, optional) – Optional: Binary parameters at which to compute the statistic.

• ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.

• tp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• argp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• record_segments (boolean) – If True, store the per-segment F-stat values as self.twoF_per_segment and (if self.singleFstats) the per-detector per-segment F-stats as self.twoFX_per_segment.

Returns:

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

Return type:

float

get_semicoherent_twoF(F0, F1, F2, Alpha, Delta, asini=None, period=None, ecc=None, tp=None, argp=None, record_segments=False)[source]

Computes the semi-coherent twoF statistic at a single point.

Parameters:
• F0 (float) – Parameters at which to compute the statistic.

• F1 (float) – Parameters at which to compute the statistic.

• F2 (float) – Parameters at which to compute the statistic.

• Alpha (float) – Parameters at which to compute the statistic.

• Delta (float) – Parameters at which to compute the statistic.

• asini (float, optional) – Optional: Binary parameters at which to compute the statistic.

• period (float, optional) – Optional: Binary parameters at which to compute the statistic.

• ecc (float, optional) – Optional: Binary parameters at which to compute the statistic.

• tp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• argp (float, optional) – Optional: Binary parameters at which to compute the statistic.

• record_segments (boolean) – If True, store the per-segment F-stat values as self.twoF_per_segment.

Returns:

twoF – A single value of the semi-coherent twoF statistic at the input parameter values. Also stored as self.twoF.

Return type:

float

get_semicoherent_single_IFO_twoFs(record_segments=False)[source]

Computes the semi-coherent single-detector F-statss at a single point.

This requires self.get_semicoherent_twoF() to be run first.

Parameters:

record_segments (boolean) – If True, store the per-detector per-segment F-stat values as self.twoFX_per_segment.

Returns:

twoFX – A list of the single-detector detection statistics twoF. Also stored as self.twoFX.

Return type:

list

get_semicoherent_log10BSGL()[source]

Computes the semi-coherent log10BSGL statistic at a single point.

This requires self.get_semicoherent_twoF() and self.get_semicoherent_single_IFO_twoFs() to be run first.

Returns:

log10BSGL – A single value of the semi-coherent log10BSGL statistic at the input parameter values. Also stored as self.log10BSGL.

Return type:

float

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

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]

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]

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

• clean (bool) – If true, ignore existing data and overwrite. Otherwise, re-use existing data if no inconsistencies are found.

tex_labels = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'lnBtSG': '$\\ln\\mathcal{B}_{\\mathrm{tS/G}}$', 'log10BSGL': '$\\log_{10}\\mathcal{B}_{\\mathrm{S/GL}}$', 'maxTwoF': '$\\max\\widetilde{2\\mathcal{F}}$', 'twoF': '$\\widetilde{2\\mathcal{F}}$'}

Formatted labels used for plot annotations.

tex_labels0 = {'Alpha': '$-\\alpha_0$', 'Delta': '$-\\delta_0$', 'F0': '$-f_0$', 'F1': '$-\\dot{f}_0$', 'F2': '$-\\ddot{f}_0$'}

Formatted labels used for annotating central values in plots.

fmt_detstat = '%.9g'

Standard output precision for detection statistics.

check_old_data_is_okay_to_use()[source]

Check if an existing output file matches this search and reuse the results.

Results will be loaded from old output file, and no new search run, if all of the following checks pass:

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.

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.

generate_loudest()[source]

Use ComputeFstatistic_v2 executable to produce a .loudest file

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

NOTE for GPU users (tCWFstatMapVersion=”pycuda”): The underlying ComputeFstat class tries to conveniently deal with GPU context management behind the scenes. A known problematic case is if you try to instantiate it twice from the same session/script. If you then get some messages like RuntimeError: make_default_context() and invalid device context, that is because the GPU is still blocked from the first instance when you try to initiate the second. To avoid this problem, use context management:

with pyfstat.TransientGridSearch(
[...],
tCWFstatMapVersion="pycuda",
) as search:
search.search.run()


or manually call the search.search.finalizer_() method where needed.

Most parameters are the same as for GridSearch and the core.ComputeFstat class, only the additional ones are documented here:

Parameters:
• transientWindowType (str) – If rect or exp, allow for the Fstat to be computed over a transient range. (none instead of None explicitly calls the transient-window function, but with the full range, for debugging.)

• t0Band (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave F-stat with t0=minStartTime, tau=maxStartTime-minStartTime.

• tauBand (int) – Search ranges for transient start-time t0 and duration tau. If >0, search t0 in (minStartTime,minStartTime+t0Band) and tau in (tauMin,2*Tsft+tauBand). If =0, only compute the continuous-wave F-stat with t0=minStartTime, tau=maxStartTime-minStartTime.

• tauMin (int) – Minimum transient duration to cover, defaults to 2*Tsft.

• dt0 (int) – Grid resolution in transient start-time, defaults to Tsft.

• dtau (int) – Grid resolution in transient duration, defaults to Tsft.

• outputTransientFstatMap (bool) – If true, write additional output files for (t0,tau) F-stat maps. (One file for each grid point!)

• outputAtoms (bool) – If true, write additional output files for the F-stat atoms. (One file for each grid point!)

• tCWFstatMapVersion (str) – Choose between implementations of the transient F-statistic funcionality: standard lal implementation, pycuda for GPU version, and some others only for devel/debug.

• cudaDeviceName (str) – GPU name to be matched against drv.Device output, only for tCWFstatMapVersion=pycuda.

run(return_data=False)[source]

Execute the actual search over the full grid.

This iterates over all points in the multi-dimensional product grid and the end result is either returned as a numpy array or saved to disk.

If the outputTransientFstatMap or outputAtoms options have been set when initiating the search, additional files are written for each frequency-evolution parameter-space point (‘Doppler’ point).

Parameters:

return_data (boolean) – If true, the final inputs+outputs data set is returned as a numpy array. If false, it is saved to disk and nothing is returned.

Returns:

data – The final inputs+outputs data set. Only if return_data=true.

Return type:

np.ndarray

get_transient_fstat_map_filename(param_point)[source]

Filename convention for given grid point: freq_alpha_delta_f1dot_f2dot

Parameters:

param_point (tuple, dict, list, np.void or np.ndarray) – A multi-dimensional parameter point. If not a type with named fields (e.g. a plain tuple or list), the order must match that of self.output_keys.

Returns:

f – The constructed filename.

Return type:

str

class pyfstat.grid_based_searches.SliceGridSearch(*args, **kwargs)[source]

Bases: DefunctClass

last_supported_version = '1.9.0'
class pyfstat.grid_based_searches.GridUniformPriorSearch(*args, **kwargs)[source]

Bases: DefunctClass

last_supported_version = '1.9.0'
class pyfstat.grid_based_searches.GridGlitchSearch(*args, **kwargs)[source]

Bases: 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: DefunctClass

last_supported_version = '1.9.0'
class pyfstat.grid_based_searches.FrequencySlidingWindow(*args, **kwargs)[source]

Bases: DefunctClass

last_supported_version = '1.9.0'
class pyfstat.grid_based_searches.EarthTest(*args, **kwargs)[source]

Bases: DefunctClass

last_supported_version = '1.9.0'
class pyfstat.grid_based_searches.DMoff_NO_SPIN(*args, **kwargs)[source]

Bases: 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:

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.injection_parameters module¶

Generate injection parameters drawn from different prior populations

class pyfstat.injection_parameters.InjectionParametersGenerator(*, seed=None, rng=NOTHING, priors=NOTHING)[source]

Bases: object

Draw injection parameter samples from priors and return in dictionary format.

Parameters:
• priors (dict) –

Each key refers to one of the signal’s parameters (following the PyFstat convention). Priors can be given as values in three formats (by order of evaluation):

1. Callable without required arguments: {“ParameterA”: np.random.uniform}.

2. Dict containing numpy.random distribution as key and kwargs in a dict as value: {“ParameterA”: {“uniform”: {“low”: 0, “high”:1}}}.

3. Constant value to be returned as is: {“ParameterA”: 1.0}.

• seed (None, int) – Argument to be fed to numpy.random.default_rng, with all of its accepted types.

• _rng (np.random.Generator) – Alternatively, this class accepts an already-initialized np.Generator, in which case the seed argument will be ignored.

Method generated by attrs for class InjectionParametersGenerator.

seed: Union[None, int]
priors: dict
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.injection_parameters.AllSkyInjectionParametersGenerator(*, seed=None, rng=NOTHING, priors=NOTHING)[source]

Draw injection parameter samples from priors and return in dictionary format. This class works in exactly the same way as InjectionParametersGenerator, but including by default two extra keys, Alpha and Delta (sky position’s right ascension and declination in radians), which are sample isotropically across the celesetial sphere.

Alpha’s distribution is Uniform(0, 2 pi), and sin(Delta)’s distribution is Uniform(-1, 1).

The only reason this class exists is because, using the notation we specified in the base class, there is no way to generate arcsine distributed numbers using a specific seed, as numpy does not have such a number generator and hence has to be constructed by applying a function to a uniform number.

Method generated by attrs for class InjectionParametersGenerator.

seed: Union[None, int]
priors: dict

## pyfstat.logging module¶

PyFstat’s logging implementation.

PyFstat main logger is called pyfstat and can be accessed via:

import logging
logger = logging.getLoger('pyfstat')


Basics of logging: For all our purposes, there are logger objects and handler objects. Loggers are the ones in charge of logging, hence you call them to emit a logging message with a specific logging level (e.g. logger.info); handlers are in charge of redirecting that message to a specific place (e.g. a file or your terminal, which is usually referred to as a stream).

The default behaviour upon importing pyfstat is to attach a logging.StreamHandler to the pyfstat logger, printing out to sys.stdout. This is only done if the root logger has no handlers attached yet; if it does have at least one handler already, then we inherit those and do not add any extra handlers by default. If, for any reason, logging cannot access sys.stdout at import time, the exception is reported via print and no handlers are attached (i.e. the logger won’t print to sys.stdout).

The user can modify the logger’s behaviour at run-time using set_up_logger. This function attaches extra logging.StreamHandler and logging.FileHandler handlers to the logger, allowing to redirect logging messages to a different stream or a specific output file specified using the outdir, label variables (with the same format as in the rest of the package).

Finally, logging can be disabled, or the level changed, at run-time by manually configuring the pyfstat logger. For example, the following block of code will suppress logging messages below WARNING::

import logging
logging.getLoger('pyfstat').setLevel(logging.WARNING)

pyfstat.logging.set_up_logger(outdir=None, label='pyfstat', log_level='INFO', streams=(<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>, ), append=True)[source]

Add file and stream handlers to the pyfstat logger.

Handler names generated from streams and outdir, label must be unique and no duplicated handler will be attached by this function.

Parameters:
• outdir (Optional[str]) – Path to outdir directory. If None, no file handler will be added.

• label (Optional[str]) – Label for the file output handler, i.e. the log file will be called label.log. Required, in conjunction with outdir, to add a file handler. Ignored otherwise.

• log_level (str) – Level of logging. This level is imposed on the logger itself and every single handler attached to it.

• streams (Optional[Iterable[TextIOWrapper]]) – Stream to which logging messages will be passed using a StreamHandler object. By default, log to sys.stdout. Other common streams include e.g. sys.stderr.

• append (bool) – If False, removes all handlers from the pyfstat logger before adding new ones. This removal is not propagated to handlers on the root logger.

Returns:

Configured instance of the logging.Logger class.

Return type:

obj

## pyfstat.make_sfts module¶

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

class pyfstat.make_sfts.Writer(*args, **kwargs)[source]

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 Makefakedata_v5 executable which will be run in a subprocess. See lalpulsar_Makefakedata_v5 –help for more detailed help with some of the parameters.

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

• tstart (int) – Starting GPS epoch of the data set. NOTE: mutually exclusive with timestamps.

• duration (int) – Duration (in GPS seconds) of the total data set. NOTE: mutually exclusive with timestamps.

• tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.

• F0 (float or None) – Frequency of a signal to inject. Also used (if Band is not None) as center of frequency band. Also needed when noise-only (h0=None or h0==0) but no noiseSFTs given, in which case it is also used as center of frequency band.

• F1 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• F2 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Alpha (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Delta (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• h0 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• cosi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• psi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• phi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.

• outdir (str) – The directory where files are written to. Default: current working directory.

• sqrtSX (float or list or str or None) – Single-sided PSD values for generating fake Gaussian noise. Single float or str value: use same for all detectors. List or comma-separated string: must match len(detectors). Detectors will be paired to list elements following alphabetical order.

• noiseSFTs (str or None) – Existing SFT files on top of which signals will be injected. If not None, additional constraints can be applied using the arguments tstart and duration. NOTE: mutually exclusive with timestamps.

• SFTWindowType (str or None) – LAL name of the windowing function to apply to the data.

• SFTWindowBeta (float) – Optional parameter for some windowing functions.

• Band (float or None) – If float, and F0 is also not None, then output SFTs cover [F0-Band/2,F0+Band/2]. If None and noiseSFTs given, use their bandwidth. If None and no noiseSFTs given, a minimal covering band for a perfectly-matched single-template ComputeFstat analysis is estimated.

• detectors (str or None) – Comma-separated list of detectors to generate data for.

• earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per utils.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 utils.get_ephemeris_files().

• transientWindowType (str) – If none, a fully persistent CW signal is simulated. If rect or exp, a transient signal with the corresponding amplitude evolution is simulated.

• transientStartTime (int or None) – Start time for a transient signal.

• transientTau (int or None) – Duration (rect case) or decay time (exp case) of a transient signal.

• randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.

• timestamps (str or dict) – Dictionary of timestamps (each key must refer to a detector), list of timestamps (detectors should be set), or comma-separated list of per-detector timestamps files (simple text files, lines with <seconds> <nanoseconds>, comments should use %, each time stamp gives the start time of one SFT). WARNING: In that last case, order must match that of detectors! NOTE: mutually exclusive with [tstart,duration] and with noiseSFTs.

mfd = 'lalpulsar_Makefakedata_v5'

The executable; can be overridden by child classes.

signal_parameter_labels = ['tref', 'F0', 'F1', 'F2', 'Alpha', 'Delta', 'h0', 'cosi', 'psi', 'phi', 'transientWindowType', 'transientStartTime', 'transientTau']

Default convention of labels for the various signal parameters.

gps_time_and_string_formats_as_LAL = {'refTime': ':10.9f', 'transientStartTime': ':10.0f', 'transientTau': ':10.0f', 'transientWindowType': ':s'}

Dictionary to ensure proper format handling for some special parameters.

GPS times should NOT be parsed using scientific notation. LAL routines would silently parse them wrongly.

required_signal_parameters = ['F0', 'Alpha', 'Delta', 'cosi']

List of parameters required for a successful execution of Makefakedata_v5. The rest of available parameters are not required as they have default values silently given by Makefakedata_v5

property tend

Defined as self.start + self.duration.

If stored as an attribute, there would be the risk of it going out of sync with the other two values.

calculate_fmin_Band()[source]

Set fmin and Band for the output SFTs to cover.

Either uses the user-provided Band and puts F0 in the middle; does nothing to later reuse the full bandwidth of noiseSFTs (only if using MFDv5); or if F0!=None, noiseSFTs=None and Band=None it estimates a minimal band for just the injected signal: F-stat covering band plus extra bins for demod default parameters. This way a perfectly matched single-template ComputeFstat analysis should run through perfectly on the returned SFTs. For any wider-band or mismatched search, one needs to set Band manually. If using MFDv4, at least F0 is required even if noiseSFTs!=None.

make_cff(verbose=False)[source]

Generates a .cff file including signal injection parameters.

This will be saved to self.config_file_name.

Parameters:

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

check_cached_data_okay_to_use(cl_mfd)[source]

Check if SFT files already exist that can be re-used.

This does not check the actual data contents of the SFTs, but only the following criteria:

• filename

• if injecting a signal, that the .cff file is older than the SFTs (but its contents should have been checked separately)

• that the commandline stored in the (first) SFT header matches

Parameters:

cl_mfd (str) – The commandline we’d execute if not finding matching files.

make_data(verbose=False)[source]

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

run_makefakedata()[source]

Generate the SFT data calling Makefakedata_v5 executable.

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 utils.predict_fstat(), this wraps the 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: 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: 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 Makefakedata_v4’s lineFeature option. This version of MFD only supports one interferometer at a time.

NOTE: All signal parameters except for h0, Freq, phi0 and transient parameters will be ignored.

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

• tstart (int) – Starting GPS epoch of the data set. NOTE: mutually exclusive with timestamps.

• duration (int) – Duration (in GPS seconds) of the total data set. NOTE: mutually exclusive with timestamps.

• tref (float or None) – Reference time for simulated signals. Default is None, which sets the reference time to tstart.

• F0 (float or None) – Frequency of a signal to inject. Also used (if Band is not None) as center of frequency band. Also needed when noise-only (h0=None or h0==0) but no noiseSFTs given, in which case it is also used as center of frequency band.

• F1 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• F2 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Alpha (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Delta (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• h0 (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• cosi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• psi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• phi (float or None) – Additional frequency evolution and amplitude parameters for a signal. If h0=None or h0=0, these are all ignored. If h0>0, then at least [Alpha,Delta,cosi] need to be set explicitly.

• Tsft (int) – The SFT duration in seconds. Will be ignored if noiseSFTs are given.

• outdir (str) – The directory where files are written to. Default: current working directory.

• sqrtSX (float or list or str or None) – Single-sided PSD values for generating fake Gaussian noise. Single float or str value: use same for all detectors. List or comma-separated string: must match len(detectors). Detectors will be paired to list elements following alphabetical order.

• noiseSFTs (str or None) – Existing SFT files on top of which signals will be injected. If not None, additional constraints can be applied using the arguments tstart and duration. NOTE: mutually exclusive with timestamps.

• SFTWindowType (str or None) – LAL name of the windowing function to apply to the data.

• SFTWindowBeta (float) – Optional parameter for some windowing functions.

• Band (float or None) – If float, and F0 is also not None, then output SFTs cover [F0-Band/2,F0+Band/2]. If None and noiseSFTs given, use their bandwidth. If None and no noiseSFTs given, a minimal covering band for a perfectly-matched single-template ComputeFstat analysis is estimated.

• detectors (str or None) – Comma-separated list of detectors to generate data for.

• earth_ephem (str or None) – Paths of the two files containing positions of Earth and Sun. If None, will check standard sources as per utils.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 utils.get_ephemeris_files().

• transientWindowType (str) – If none, a fully persistent CW signal is simulated. If rect or exp, a transient signal with the corresponding amplitude evolution is simulated.

• transientStartTime (int or None) – Start time for a transient signal.

• transientTau (int or None) – Duration (rect case) or decay time (exp case) of a transient signal.

• randSeed (int or None) – Optionally fix the random seed of Gaussian noise generation for reproducibility.

• timestamps (str or dict) – Dictionary of timestamps (each key must refer to a detector), list of timestamps (detectors should be set), or comma-separated list of per-detector timestamps files (simple text files, lines with <seconds> <nanoseconds>, comments should use %, each time stamp gives the start time of one SFT). WARNING: In that last case, order must match that of detectors! NOTE: mutually exclusive with [tstart,duration] and with noiseSFTs.

mfd = 'lalpulsar_Makefakedata_v4'

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

required_signal_parameters = ['F0', 'phi', 'h0']

Required parameters for Makefakedata_v4 to success. Any other parameter is silently given a default value by Makefakedata_v4.

signal_parameters_labels = ['F0', 'phi', 'h0', 'transientWindowType', 'transientStartTime', 'transientTau']

Other signal parameters will be removed before passing to Makefakedata_v4.

class pyfstat.make_sfts.GlitchWriter(*args, **kwargs)[source]

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: Writer

Specialized Writer variant to generate SFTs containing simulated instrumental artifacts.

Contrary to the main Writer class, this calls the older Makefakedata_v4 executable which supports the special –lineFeature option. See lalpulsar_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 utils.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 utils.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 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.

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.

Parameters:

num_processes (int) – Number threads to use when running in parallel. Verbatim implementation of the former args.N.

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]

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

class pyfstat.make_sfts.InjectionParametersGenerator(*, seed=None, rng=NOTHING, priors=NOTHING)[source]

Method generated by attrs for class InjectionParametersGenerator.

seed: Union[None, int]
priors: dict
class pyfstat.make_sfts.AllSkyInjectionParametersGenerator(*, seed=None, rng=NOTHING, priors=NOTHING)[source]

Method generated by attrs for class InjectionParametersGenerator.

seed: Union[None, int]
priors: dict

## pyfstat.mcmc_based_searches module¶

PyFstat search & follow-up classes using MCMC-based methods

The general approach is described in Ashton & Prix (PRD 97, 103020, 2018): https://arxiv.org/abs/1802.05450 and we use the ptemcee sampler described in Vousden et al. (MNRAS 455, 1919-1937, 2016): https://arxiv.org/abs/1501.05823 and based on Foreman-Mackey et al. (PASP 125, 306, 2013): https://arxiv.org/abs/1202.3665

### Defining the prior¶

The MCMC based searches (i.e. pyfstat.MCMC*) require a prior specification for each model parameter, implemented via a python dictionary. This is best explained through a simple example, here is the prior for a directed search with a uniform prior on the frequency and a normal prior on the frequency derivative:

theta_prior = {'F0': {'type': 'unif',
'lower': 29.9,
'upper': 30.1},
'F1': {'type': 'norm',
'loc': 0,
'scale': 1e-10},
'F2': 0,
'Alpha': 2.3,
'Delta': 1.8
}


For the sky positions Alpha and Delta, we give the fixed values (i.e. they are considered known by the MCMC simulation), the same is true for F2, the second derivative of the frequency which we fix at 0. Meanwhile, for the frequency F0 and first frequency derivative F1 we give a dictionary specifying their prior distribution. This dictionary must contain three arguments: the type (in this case either unif or norm) which specifies the type of distribution, then two shape arguments. The shape parameters will depend on the type of distribution, but here we use lower and upper, required for the unif prior while loc and scale are required for the norm prior.

Currently, two other types of prior are implemented: halfnorm, neghalfnorm (both of which require loc and scale shape parameters). Further priors can be added by modifying pyfstat.MCMCSearch._generic_lnprior.

class pyfstat.mcmc_based_searches.MCMCSearch(*args, **kwargs)[source]

MCMC search using ComputeFstat.

Evaluates the coherent F-statistic across a parameter space region corresponding to an isolated/binary-modulated CW signal.

Parameters:
• theta_prior (dict) – Dictionary of priors and fixed values for the search parameters. For each parameters (key of the dict), if it is to be held fixed the value should be the constant float, if it is be searched, the value should be a dictionary of the prior.

• tref (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• minStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• maxStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• label (str) – A label and output directory (optional, default is data) to name files

• outdir (str) – A label and output directory (optional, default is data) to name files

• sftfilepattern (str, optional) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.

• detectors (str, optional) – Two character reference to the detectors to use, specify None for no contraint and comma separated strings for multiple references.

• nsteps (list (2,), optional) – Number of burn-in and production steps to take, [nburn, nprod]. See pyfstat.MCMCSearch.setup_initialisation() for details on adding initialisation steps.

• nwalkers (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.

• ntemps (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.

• log10beta_min (float < 0, optional) – The log_10(beta) value. If given, the set of betas passed to PTSampler are generated from np.logspace(0, log10beta_min, ntemps) (given in descending order to ptemcee).

• theta_initial (dict, array, optional) – A dictionary of distribution about which to distribute the initial walkers about.

• rhohatmax (float, optional) – Upper bound for the SNR scale parameter (required to normalise the Bayes factor) - this needs to be carefully set when using the evidence.

• binary (bool, optional) – If true, search over binary orbital parameters.

• BSGL (bool, optional) – If true, use the BSGL statistic.

• BtSG (bool, optional) – If true, use the transient lnBtSG statistic. (Only for transient searches.)

• SSBPrec (int, optional) – SSBPrec (SSB precision) to use when calling ComputeFstat. See core.ComputeFstat.

• RngMedWindow (int, optional) – Running-Median window size (number of bins) for ComputeFstat. See core.ComputeFstat.

• minCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.

• maxCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.

• injectSources (dict, optional) – If given, inject these properties into the SFT files before running the search. See core.ComputeFstat.

• assumeSqrtSX (float or list or str) – Don’t estimate noise-floors, but assume (stationary) per-IFO sqrt{SX}. See core.ComputeFstat.

• transientWindowType (str) – If ‘rect’ or ‘exp’, compute atoms so that a transient (t0,tau) map can later be computed. (‘none’ instead of None explicitly calls the transient-window function, but with the full range, for debugging). See core.ComputeFstat. Currently only supported for nsegs=1.

• tCWFstatMapVersion (str) – Choose between standard ‘lal’ implementation, ‘pycuda’ for gpu, and some others for devel/debug.

• allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength].

• clean (bool) – If true, ignore existing data and overwrite. Otherwise, re-use existing data if no inconsistencies are found.

symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'argp': 'argp', 'asini': 'asini', 'ecc': 'ecc', 'period': 'P', 'tp': 'tp'}

Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots

unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'argp': '', 'asini': '', 'ecc': '', 'period': 's', 'tp': 's'}

Key, val pairs of the parameters (i.e. F0, F1), and the units (i.e. Hz)

transform_dictionary = {}

Key, val pairs of the parameters (i.e. F0, F1), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.

setup_initialisation(nburn0, scatter_val=1e-10)[source]

Add an initialisation step to the MCMC run

If called prior to run(), adds an intial step in which the MCMC simulation is run for nburn0 steps. After this, the MCMC simulation continues in the usual manner (i.e. for nburn and nprod steps), but the walkers are reset scattered around the maximum likelihood position of the initialisation step.

Parameters:
• nburn0 (int) – Number of initialisation steps to take.

• scatter_val (float) – Relative number to scatter walkers around the maximum likelihood position after the initialisation step. If the maximum likelihood point is located at p, the new walkers are randomly drawn from a multivariate gaussian distribution centered at p with standard deviation diag(scatter_val * p).

run(proposal_scale_factor=2, save_pickle=True, export_samples=True, save_loudest=True, plot_walkers=True, walker_plot_args=None, window=50)[source]

Run the MCMC simulatation

Parameters:
• proposal_scale_factor (float) – The proposal scale factor a > 1 used by the sampler. See Goodman & Weare (Comm App Math Comp Sci, Vol 5, No. 1, 2010): 10.2140/camcos.2010.5.65. The bigger the value, the wider the range to draw proposals from. If the acceptance fraction is too low, you can raise it by decreasing the a parameter; and if it is too high, you can reduce it by increasing the a parameter. See Foreman-Mackay et al. (PASP 125 306, 2013): https://arxiv.org/abs/1202.3665.

• save_pickle (bool) – If true, save a pickle file of the full sampler state.

• export_samples (bool) – If true, save ASCII samples file to disk. See MCMCSearch.export_samples_to_disk.

• save_loudest (bool) – If true, save a CFSv2 .loudest file to disk. See MCMCSearch.generate_loudest.

• plot_walkers (bool) – If true, save trace plots of the walkers.

• walker_plot_args – Dictionary passed as kwargs to _plot_walkers to control the plotting. Histogram of sampled detection statistic values can be retrieved setting “plot_det_stat” to True. Parameters corresponding to an injected signal can be passed through “injection_parameters” as a dictionary containing the parameters of said signal. All parameters being searched for must be present, otherwise this option is ignored. If both “fig” and “axes” entries are set, the plot is not saved to disk directly, but (fig, axes) are returned.

• window (int) – The minimum number of autocorrelation times needed to trust the result when estimating the autocorrelation time (see ptemcee.Sampler.get_autocorr_time for further details.

plot_corner(figsize=(10, 10), add_prior=False, nstds=None, label_offset=0.4, dpi=300, rc_context={}, tglitch_ratio=False, fig_and_axes=None, save_fig=True, **kwargs)[source]

Generate a corner plot of the posterior

Using the corner package (https://pypi.python.org/pypi/corner/), generate estimates of the posterior from the production samples.

Parameters:
• figsize (tuple (7, 7)) – Figure size in inches (passed to plt.subplots)

• add_prior (bool, str) – If true, plot the prior as a red line. If ‘full’ then for uniform priors plot the full extent of the prior.

• nstds (float) – The number of standard deviations to plot centered on the median. Standard deviation is computed from the samples using numpy.std.

• label_offset (float) – Offset the labels from the plot: useful to prevent overlapping the tick labels with the axis labels. This option is passed to ax.[x|y]axis.set_label_coords.

• dpi (int) – Passed to plt.savefig.

• rc_context (dict) – Dictionary of rc values to set while generating the figure (see matplotlib rc for more details).

• tglitch_ratio (bool) – If true, and tglitch is a parameter, plot posteriors as the fractional time at which the glitch occurs instead of the actual time.

• fig_and_axes (tuple) – (fig, axes) tuple to plot on. The axes must be of the right shape, namely (ndim, ndim)

• save_fig (bool) – If true, save the figure, else return the fig, axes.

• **kwargs – Passed to corner.corner. Use “truths” to plot the true parameters of a signal.

Returns:

The matplotlib figure and axes, only returned if save_fig = False.

Return type:

fig, axes

plot_chainconsumer(save_fig=True, label_offset=0.25, dpi=300, **kwargs)[source]

Generate a corner plot of the posterior using the chaniconsumer package.

chainconsumer is an optional dependency of PyFstat. See https://samreay.github.io/ChainConsumer/.

Parameters are akin to the ones described in MCMCSearch.plot_corner. Only the differing parameters are explicitly described.

Parameters:

**kwargs – Passed to chainconsumer.plotter.plot. Use “truths” to plot the true parameters of a signal.

plot_prior_posterior(normal_stds=2, injection_parameters=None, fig_and_axes=None, save_fig=True)[source]

Plot the prior and posterior probability distributions in the same figure

Parameters:
• normal_stds (int) – Bounds of priors in terms of their standard deviation. Only used if norm, halfnorm, neghalfnorm or lognorm priors are given, otherwise ignored.

• injection_parameters (dict) – Dictionary containing the parameters of a signal. All parameters being searched must be present as dictionary keys, otherwise this option is ignored.

• fig_and_axes (tuple) – (fig, axes) tuple to plot on.

• save_fig (bool) – If true, save the figure, else return the fig, axes.

Returns:

(fig, ax) – If save_fig evaluates to False, return figure and axes.

Return type:

(matplotlib.pyplot.figure, matplotlib.pyplot.axes)

plot_cumulative_max(**kwargs)[source]

Plot the cumulative twoF for the maximum posterior estimate.

This method accepts the same arguments as pyfstat.core.ComputeFstat.plot_twoF_cumulative, except for CFS_input, which is taken from the loudest candidate; and label and outdir, which are taken from the instance of this class.

For example, one can pass signal arguments to predic_twoF_cumulative through PFS_kwargs, or set the number of segments using num_segments_(CFS|PFS). The same applies for other options such as tstart, tend or savefig. Every single of these arguments will be passed to pyfstat.core.ComputeFstat.plot_twoF_cumulative as they are, using their default argument otherwise.

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

Unlike the core function, here savefig=True is the default, for consistency with other MCMC plotting functions.

get_saved_data_dictionary()[source]

Read the data saved in self.pickel_path and return it as a dictionary.

Returns:

d – Dictionary containing the data saved in the pickle self.pickle_path.

Return type:

dict

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 ComputeFstatistic_v2 executable 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.

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: MCMCSearch

MCMC search using the SemiCoherentGlitchSearch

See parent MCMCSearch for a list of all additional parameters, here we list only the additional init parameters of this class.

Parameters:
• nglitch (int) – The number of glitches to allow

• dtglitchmin (int) – The minimum duration (in seconds) of a segment between two glitches or a glitch and the start/end of the data

• theta0_idx – Index (zero-based) of which segment the theta refers to - useful if providing a tight prior on theta to allow the signal to jump too theta (and not just from)

• int – Index (zero-based) of which segment the theta refers to - useful if providing a tight prior on theta to allow the signal to jump too theta (and not just from)

symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'delta_F0': '$\\delta f$', 'delta_F1': '$\\delta \\dot{f}$', 'tglitch': '$t_\\mathrm{glitch}$'}

Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots

glitch_symbol_dictionary = {'delta_F0': '$\\delta f$', 'delta_F1': '$\\delta \\dot{f}$', 'tglitch': '$t_\\mathrm{glitch}$'}

Key, val pairs of glitch parameters (dF0, dF1, tglitch), to LaTeX math symbols for plots. This dictionary included within self.symbol_dictionary.

unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'delta_F0': 'Hz', 'delta_F1': 'Hz/s', 'tglitch': 's'}

Key, val pairs of the parameters (F0, F1, …, including glitch parameters), and the units (Hz, Hz/s, …).

transform_dictionary = {'tglitch': {'label': '$t^{g}_0$ \\n [d]', 'multiplier': 1.1574074074074073e-05, 'subtractor': 'minStartTime', 'unit': 'day'}}

Key, val pairs of the parameters (F0, F1, …), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.

plot_cumulative_max(savefig=True)[source]

Override MCMCSearch.plot_cumulative_max implementation to deal with the split at glitches.

Parameters:

savefig (boolean) – included for consistency with core plot_twoF_cumulative() function. If true, save the figure in outdir. If false, return an axis object.

class pyfstat.mcmc_based_searches.MCMCSemiCoherentSearch(*args, **kwargs)[source]

Bases: 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]

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, setup_only=False, no_template_counting=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 logger.info as a table.

• gen_tex_table (bool) – Dump follow-up setup into a text file as a tex table. File is constructed as os.path.join(self.outdir, self.label + “_run_setup.tex”).

Returns:

run_setup – List containing the setup of the follow-up run.

Return type:

list

class pyfstat.mcmc_based_searches.MCMCTransientSearch(*args, **kwargs)[source]

Bases: MCMCSearch

MCMC search for a transient signal using ComputeFstat

Parameters:
• theta_prior (dict) – Dictionary of priors and fixed values for the search parameters. For each parameters (key of the dict), if it is to be held fixed the value should be the constant float, if it is be searched, the value should be a dictionary of the prior.

• tref (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• minStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• maxStartTime (int) – GPS seconds of the reference time, start time and end time. While tref is requirede, minStartTime and maxStartTime default to None in which case all available data is used.

• label (str) – A label and output directory (optional, default is data) to name files

• outdir (str) – A label and output directory (optional, default is data) to name files

• sftfilepattern (str, optional) – Pattern to match SFTs using wildcards (*?) and ranges [0-9]; mutiple patterns can be given separated by colons.

• detectors (str, optional) – Two character reference to the detectors to use, specify None for no contraint and comma separated strings for multiple references.

• nsteps (list (2,), optional) – Number of burn-in and production steps to take, [nburn, nprod]. See pyfstat.MCMCSearch.setup_initialisation() for details on adding initialisation steps.

• nwalkers (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.

• ntemps (int, optional) – The number of walkers and temperates to use in the parallel tempered PTSampler.

• log10beta_min (float < 0, optional) – The log_10(beta) value. If given, the set of betas passed to PTSampler are generated from np.logspace(0, log10beta_min, ntemps) (given in descending order to ptemcee).

• theta_initial (dict, array, optional) – A dictionary of distribution about which to distribute the initial walkers about.

• rhohatmax (float, optional) – Upper bound for the SNR scale parameter (required to normalise the Bayes factor) - this needs to be carefully set when using the evidence.

• binary (bool, optional) – If true, search over binary orbital parameters.

• BSGL (bool, optional) – If true, use the BSGL statistic.

• BtSG (bool, optional) – If true, use the transient lnBtSG statistic. (Only for transient searches.)

• SSBPrec (int, optional) – SSBPrec (SSB precision) to use when calling ComputeFstat. See core.ComputeFstat.

• RngMedWindow (int, optional) – Running-Median window size (number of bins) for ComputeFstat. See core.ComputeFstat.

• minCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.

• maxCoverFreq (float, optional) – Minimum and maximum instantaneous frequency which will be covered over the SFT time span as passed to CreateFstatInput. See core.ComputeFstat.

• injectSources (dict, optional) – If given, inject these properties into the SFT files before running the search. See core.ComputeFstat.

• assumeSqrtSX (float or list or str) – Don’t estimate noise-floors, but assume (stationary) per-IFO sqrt{SX}. See core.ComputeFstat.

• transientWindowType (str) – If ‘rect’ or ‘exp’, compute atoms so that a transient (t0,tau) map can later be computed. (‘none’ instead of None explicitly calls the transient-window function, but with the full range, for debugging). See core.ComputeFstat. Currently only supported for nsegs=1.

• tCWFstatMapVersion (str) – Choose between standard ‘lal’ implementation, ‘pycuda’ for gpu, and some others for devel/debug.

• allowedMismatchFromSFTLength (float) – Maximum allowed mismatch from SFTs being too long [Default: what’s hardcoded in XLALFstatMaximumSFTLength].

• clean (bool) – If true, ignore existing data and overwrite. Otherwise, re-use existing data if no inconsistencies are found.

symbol_dictionary = {'Alpha': '$\\alpha$', 'Delta': '$\\delta$', 'F0': '$f$', 'F1': '$\\dot{f}$', 'F2': '$\\ddot{f}$', 'transient_duration': '$\\Delta T$', 'transient_tstart': '$t_\\mathrm{start}$'}

Key, val pairs of the parameters (F0, F1, …), to LaTeX math symbols for plots

unit_dictionary = {'Alpha': 'rad', 'Delta': 'rad', 'F0': 'Hz', 'F1': 'Hz/s', 'F2': 'Hz/s$^2$', 'transient_duration': 's', 'transient_tstart': 's'}

Key, val pairs of the parameters (F0, F1, …, including glitch parameters), and the units (Hz, Hz/s, …).

transform_dictionary = {'transient_duration': {'multiplier': 1.1574074074074073e-05, 'symbol': 'Transient duration', 'unit': 'day'}, 'transient_tstart': {'label': 'Transient start-time \n days after minStartTime', 'multiplier': 1.1574074074074073e-05, 'subtractor': 'minStartTime', 'unit': 'day'}}

Key, val pairs of the parameters (F0, F1, …), where the key is itself a dictionary which can item multiplier, subtractor, or unit by which to transform by and update the units.

## pyfstat.optimal_setup_functions module¶

Provides functions to aid in calculating the optimal setup for zoom follow up

pyfstat.optimal_setup_functions.get_optimal_setup(NstarMax, Nsegs0, tref, minStartTime, maxStartTime, prior, detector_names)[source]

Using an optimisation step, calculate the optimal setup ladder

The details of the methods are described in Sec Va of arXiv:1802.05450. Here we provide implementation details. All equation numbers refer to arXiv:1802.05450.

Parameters:
• NstarMax (float) – The ratio of the size at the old coherence time to the new coherence time for each step, see Eq. (31). Larger values allow a more rapid “zoom” of the search space at the cost of convergence. Smaller values are more conservative at the cost of additional computing time. The exact choice should be optimized for the problem in hand, but values of 100-1000 are typically okay.

• Nsegs0 (int) – The number of segments for the initial step of the ladder

• tref (int) – GPS times of the reference, start, and end time.

• minStartTime (int) – GPS times of the reference, start, and end time.

• maxStartTime (int) – GPS times of the reference, start, and end time.

• prior (dict) – Prior dictionary, each item must either be a fixed scalar value, or a uniform prior.

• detector_names (list or str) – Names of the detectors to use

Returns:

nsegs, Nstar – Ladder of segment numbers and the corresponding Nstar

Return type:

list

pyfstat.optimal_setup_functions.get_Nstar_estimate(nsegs, tref, minStartTime, maxStartTime, prior, detector_names)[source]

Returns N* estimated from the super-sky metric

Parameters:
• nsegs (int) – Number of semi-coherent segments

• tref (int) – Reference time in GPS seconds

• minStartTime (int) – Minimum and maximum SFT timestamps

• maxStartTime (int) – Minimum and maximum SFT timestamps

• prior (dict) – The prior dictionary

• detector_names (array) – Array of detectors to average over

Returns:

Nstar – The estimated approximate number of templates to cover the prior parameter space at a mismatch of unity, assuming the normalised thickness is unity.

Return type:

int

## pyfstat.snr module¶

class pyfstat.snr.SignalToNoiseRatio(*, detector_states, noise_weights=None, assumeSqrtSX=None)[source]

Bases: object

Compute the optimal SNR of a CW signal as expected in Gaussian noise.

The definition of SNR (shortcut for “optimal signal-to-noise ratio”) is taken from Eq. (76) of https://dcc.ligo.org/T0900149-v6/public and is such that $$\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2$$, where $$\langle 2\mathcal{F}\rangle$$ represents the expected value over noise realizations of twice the F-statistic of a template perfectly matched to an existing signal in the data.

Computing $$\textrm{SNR}^2$$ requires two quantities:

• The antenna pattern matrix $$\mathcal{M}$$, which depends on the sky position $$\vec{n}$$
and polarization angle $$\psi$$ and encodes the effect of the detector’s antenna pattern response
over the course of the observing run.
• The JKS amplitude parameters $$(\mathcal{A}^0, \mathcal{A}^1, \mathcal{A}^2, \mathcal{A}^3)$$
[JKS1998] which are functions of the CW’s amplitude parameters $$(h_0,\cos\iota, \psi, \phi_0)$$ or,
alternatively, $$(A_{+}, A_{\times}, \psi, \phi_0)$$.
Parameters:
• detector_states (lalpulsar.MultiDetectorStateSeries) – MultiDetectorStateSeries as produced by DetectorStates. Provides the required information to compute the antenna pattern contribution.

• noise_weights (Union[lalpulsar.MultiNoiseWeights, None]) – Optional, incompatible with assumeSqrtSX. Can be computed from SFTs using SignalToNoiseRatio.from_sfts. Noise weights to account for a varying noise floor or unequal noise floors in different detectors.

• assumeSqrtSX (float) – Optional, incompatible with noise_weights. Single-sided amplitude spectral density (ASD) of the detector noise. This value is used for all detectors, meaning it’s not currently possible to manually specify different noise floors without creating SFT files. (To be improved in the future; developer note: will require SWIG constructor for MultiNoiseWeights.)

Method generated by attrs for class SignalToNoiseRatio.

detector_states: MultiDetectorStateSeries
noise_weights: Optional[MultiNoiseWeights]
assumeSqrtSX: float
classmethod from_sfts(F0, sftfilepath, time_offset=None, running_median_window=101, sft_constraint=None)[source]

Alternative constructor to retrieve detector states and noise weights from SFT files. This method is based on DetectorStates.multi_detector_states_from_sfts(). This is currently the other way in which varying / different noise floors can be used when computing SNRs.

Parameters:
• F0 (float) – Central frequency [Hz] to retrieve from the SFT files to compute noise weights.

• sftfilepath (str) – Path to SFT files in a format compatible with XLALSFTdataFind.

• time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).

• running_median_window (int) – Window used to compute the running-median noise floor estimation. Default value is consistent with that used in PredictFstat executable.

• sft_constraint (lalpulsar.SFTConstraint) – Optional argument to specify further constraints in XLALSFTdataFind.

compute_snr2(Alpha, Delta, psi, phi, h0=None, cosi=None, aPlus=None, aCross=None)[source]

Compute the $$\textrm{SNR}^2$$ of a CW signal using XLALComputeOptimalSNR2FromMmunu. Parameters correspond to the standard ones used to describe a CW (see e.g. Eqs. (16), (26), (30) of https://dcc.ligo.org/T0900149-v6/public ).

Mind that this function returns squared SNR (Eq. (76) of https://dcc.ligo.org/T0900149-v6/public ), which can be directly related to the expected F-statistic as $$\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2$$.

Parameters:
• Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.

• Delta (float) – Declination (equatorial latitude) of the signal in radians.

• psi (float) – Polarization angle.

• h0 (float) – Nominal GW amplitude. Must be given together with cosi and conflicts with aPlus and aCross.

• cosi (float) – Cosine of the source inclination w.r.t. line of sight. Must be given together with h0 and conflicts with aPlus and aCross.

• aPlus (float) – Plus polarization amplitude. Must be given with aCross and conflicts with h0 and cosi.

• aCross (float) – Cross polarization amplitude. Must be given with aPlus and conflicts with h0 and cosi.

Returns:

SNR^2 – Squared signal-to-noise ratio of a CW signal consistent with the specified parameters in the specified detector network.

Return type:

float

compute_h0_from_snr2(Alpha, Delta, psi, phi, cosi, snr2)[source]

Convert the $$\textrm{SNR}^2$$ of a CW signal to a corresponding amplitude $$h_0$$ given the source orientation. Parameters correspond to the standard ones used to describe a CW (see e.g. Eqs. (16), (26), (30) of https://dcc.ligo.org/T0900149-v6/public ).

This function returns “inverts” Eq. (77) of https://dcc.ligo.org/T0900149-v6/public by computing the overall prefactor on $$h_0$$ using self.compute_snr2(h0=1, …).

Parameters:
• Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.

• Delta (float) – Declination (equatorial latitude) of the signal in radians.

• psi (float) – Polarization angle.

• cosi (float) – Cosine of the source inclination w.r.t. line of sight. Must be given together with h0 and conflicts with aPlus and aCross.

• snr2 (float) – Squared signal-to-noise ratio of a CW signal in the specified detector network.

Returns:

h0 – Nominal GW amplitude.

Return type:

float

compute_twoF(*args, **kwargs)[source]

Compute the expected $$2\mathcal{F}$$ value of a CW signal from the result of compute_snr2.

$\langle 2\mathcal{F}\rangle = 4 + \textrm{SNR}^2$
$\sigma_{2\mathcal{F}} = \sqrt{8 + 4 \textrm{SNR}^2}$

Input parameters are passed untouched to self.compute_snr2. See corresponding docstring for a list of valid parameters.

Returns:

• expected_2F – Expected value of a non-central chi-squared distribution with four degrees of freedom and non-centrality parameter given by SNR^2.

• stdev_2F – Standard deviation of a non-central chi-squared distribution with four degrees of freedom and non-centrality parameter given by SNR^2.

compute_Mmunu(Alpha, Delta)[source]

Compute Mmunu matrix at a specific sky position using the detector states (and possible noise weights) given at initialization time. If no noise weigths were given, unit weights are assumed.

Parameters:
• Alpha (float) – Right ascension (equatorial longitude) of the signal in radians.

• Delta (float) – Declination (equatorial latitude) of the signal in radians.

Returns:

Mmunu – Mmunu matrix encoding the response of the given detector network to a CW at the specified sky position.

Return type:

lalpulsar.AntennaPatternMatrix

class pyfstat.snr.DetectorStates[source]

Bases: object

Python interface to XLALGetMultiDetectorStates and XLALGetMultiDetectorStatesFromMultiSFTs.

get_multi_detector_states(timestamps, Tsft, detectors=None, time_offset=None)[source]
Parameters:
• timestamps (array-like or dict) – GPS timestamps at which detector states will be retrieved. If array, use the same set of timestamps for all detectors, which must be explicitly given by the user via detectors. If dictionary, each key should correspond to a valid detector name to be parsed by XLALParseMultiLALDetector and the associated value should be an array-like set of GPS timestamps for each individual detector.

• Tsft (float) – Timespan covered by each timestamp. It does not need to coincide with the separation between consecutive timestamps.

• detectors (list[str] or comma-separated string) – List of detectors to be parsed using XLALParseMultiLALDetector. Conflicts with dictionary of timestamps, required otherwise.

• time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).

Returns:

multi_detector_states – Resulting multi-detector states produced by XLALGetMultiDetectorStates

Return type:

lalpulsar.MultiDetectorStateSeries

get_multi_detector_states_from_sfts(sftfilepath, central_frequency, time_offset=None, frequency_wing_bins=1, sft_constraint=None, return_sfts=False)[source]
Parameters:
• sftfilepath (str) – Path to SFT files in a format compatible with XLALSFTdataFind.

• central_frequency (float) – Frequency [Hz] around which SFT data will be retrieved. This option is only relevant if further information is to be retrieved from the SFTs (i.e. return_sfts=True).

• time_offset (float) – Timestamp offset to retrieve detector states. Defaults to LALSuite’s default of using the central time of an STF (SFT’s timestamp + Tsft/2).

• frequency_wing_bins (int) – Frequency bins around the central frequency to retrieve from SFT data. Bin size is determined using the SFT baseline time as obtained from the catalog. This option is only relevant if further information is to be retrieved from the SFTs (i.e. return_sfts=True).

• sft_constraint (lalpulsar.SFTConstraint) – Optional argument to specify further constraints in XLALSFTdataFind.

• return_sfts (bool) – If True, also return the loaded SFTs. This is useful to compute further quantities such as noise weights.

Returns:

• multi_detector_states (lalpulsar.MultiDetectorStateSeries) – Resulting multi-detector states produced by XLALGetMultiDetectorStatesFromMultiSFTs

• multi_sfts (lalpulsar.MultiSFTVector) – Only if return_sfts is True. MultiSFTVector produced by XLALLoadMultiSFTs along the specified frequency band.

## pyfstat.tcw_fstat_map_funcs module¶

Additional helper functions dealing with transient-CW F(t0,tau) maps.

See Prix, Giampanis & Messenger (PRD 84, 023007, 2011): https://arxiv.org/abs/1104.1704 for the algorithm in general and Keitel & Ashton (CQG 35, 205003, 2018): https://arxiv.org/abs/1805.05652 for a detailed discussion of the GPU implementation.

class pyfstat.tcw_fstat_map_funcs.pyTransientFstatMap(N_t0Range=None, N_tauRange=None, transientFstatMap_t=None, from_file=None)[source]

Bases: object

Simplified object class for a F(t0,tau) F-stat map.

This is based on LALSuite’s transientFstatMap_t type, replacing the gsl matrix with a numpy array.

Here, t0 is a transient start time, tau is a transient duration parameter, and F(t0,tau) is the F-statistic (not 2F)! evaluated for a signal with those parameters (and an implicit window function, which is not stored inside this object).

The ‘map’ covers a range of different (t0,tau) pairs.

F_mn

2D array of F values (not 2F!), with m an index over start-times t0, and n an index over duration parameters tau, in steps of dt0 in [t0, t0+t0Band], and dtau in [tau, tau+tauBand].

Type:

np.ndarray

maxF

Maximum of F (not 2F!) over the array.

Type:

float

t0_ML

Maximum likelihood estimate of the transient start time t0.

Type:

float

tau_ML

Maximum likelihood estimate of the transient duration tau.

Type:

float

lnBtSG

Natural log of the marginalised transient Bayes factor. NOTE: This is always initialised as nan, and you have to call get_lnBtSG() to get its actual value.

Type:

float

The class can be initialized with either a pair of (N_t0Range,N_tauRange), from a lalpulsar object, or reading from a file.

Parameters:
• N_t0Range (int) – Number of t0 values covered.

• N_tauRange (int) – Number of tau values covered.

• transientFstatMap_t (lalpulsar.transientFstatMap_t) – pre-allocated matrix from lalpulsar to initialise from.

• from_file (str or None) – Text file, compatible with lalpulsar.write_transientFstatMap_to_fp() format, to load and initialise from.

Read F_mn map from a text file and set all other fields.

Apart from optional header lines (# comments), the format has to be consistent with lalpulsar.write_transientFstatMap_to_fp() and the write_F_mn_to_file() method of this class itself: with the columns [t0[s], tau[s], 2F]. NOTE that the file is expected to provide 2F, so the values will be halved to obtain F for storage in this class.

Parameters:

file (str) – Name of the file to load from.

get_maxF_idx()[source]

Gets the 2D-unravelled index pair of the maximum of the F_mn map

Returns:

idx – The m,n indices of the map entry with maximal F value.

Return type:

tuple

get_lnBtSG()[source]

Compute (natural log of the) transient-CW Bayes-factor B_tSG = P(x|HyptS)/P(x|HypG).

Here HypG = Gaussian noise hypothesis, HyptS = transient-CW signal hypothesis.

B_tSG is marginalized over start-time and timescale of transient CW signal, using given type and parameters of transient window range.

This is a python port of the lalpulsar.ComputeTransientBstat implementation, replacing for loops by numpy operations.

Format a 2D transient-F-stat matrix over (t0,tau) and write as a text file.

Apart from the optional extra header lines, the format is consistent with lalpulsar.write_transientFstatMap_to_fp(), with the columns [t0[s], tau[s], 2F]. NOTE that the output is 2F, not F like stored in this class itself!

Parameters:
• tCWfile (str) – Name of the file to write to.

• windowRange (lalpulsar.transientWindowRange_t) – A lalpulsar structure containing the transient parameters.

• header (list) – A list of additional header lines to print at the start of the file.

pyfstat.tcw_fstat_map_funcs.fstatmap_versions = {'lal': <function <lambda>>, 'pycuda': <function <lambda>>}

Dictionary of the actual callable transient F-stat map functions this module supports.

Actual runtime availability depends on the corresponding external modules being available.

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:
• feature (str) – Set the transient F-stat map implementation.

• cudaDeviceName (str or None) – Request a CUDA device with this name. Partial matches are allowed.

Returns:

• features (dict) – A dictionary of available method names, to match fstatmap_versions. Each key’s value is set to True only if all required modules are importable on this system.

• gpu_context (pycuda.driver.Context or None) – A CUDA device context object, if assigned.

pyfstat.tcw_fstat_map_funcs.call_compute_transient_fstat_map(version, features, multiFstatAtoms=None, windowRange=None, BtSG=False)[source]

Call a version of the ComputeTransientFstatMap function.

This checks if the requested version is available, and if so, executes the computation of a transient F-statistic map over the windowRange.

Parameters:
• version (str) – Name of the method to call (currently supported: ‘lal’ or ‘pycuda’).

• features (dict) – Dictionary of available features, as obtained from init_transient_fstat_map_features().

• multiFstatAtoms (lalpulsar.MultiFstatAtomVector or None) – The time-dependent F-stat atoms previously computed by ComputeFstat.

• windowRange (lalpulsar.transientWindowRange_t or None) – The structure defining the transient parameters.

• BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using the appropriate implementation for each feature, and store it in FstatMap.lnBtSG.

Returns:

• FstatMap (pyTransientFstatMap or lalpulsar.transientFstatMap_t) – The output of the called function, including the evaluated transient F-statistic map over the windowRange.

• timingFstatMap (float) – Execution time of the called function.

pyfstat.tcw_fstat_map_funcs.lalpulsar_compute_transient_fstat_map(multiFstatAtoms, windowRange, BtSG=False)[source]

Wrapper for the standard lalpulsar function for computing a transient F-statistic map.

See https://lscsoft.docs.ligo.org/lalsuite/lalpulsar/_transient_c_w__utils_8h.html for the wrapped function.

Parameters:
• multiFstatAtoms (lalpulsar.MultiFstatAtomVector) – The time-dependent F-stat atoms previously computed by ComputeFstat.

• windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.

• BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using the corresponding lalpulsar function, and store it in FstatMap.lnBtSG.

Returns:

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

Return type:

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

GPU version of computing a transient F-statistic map.

This is based on XLALComputeTransientFstatMap from LALSuite, (C) 2009 Reinhard Prix, licensed under GPL.

The ‘map’ consists of F-statistics evaluated over a range of different (t0,tau) pairs (transient start-times and duration parameters).

This is a high-level wrapper function; the actual CUDA computations are performed in one of the functions pycuda_compute_transient_fstat_map_rect() or pycuda_compute_transient_fstat_map_exp(), depending on the window functon defined in windowRange.

Parameters:
• multiFstatAtoms (lalpulsar.MultiFstatAtomVector) – The time-dependent F-stat atoms previously computed by ComputeFstat.

• windowRange (lalpulsar.transientWindowRange_t) – The structure defining the transient parameters.

• BtSG (boolean) – If true, also compute the lnBtSG transient Bayes factor statistic, using a CPU python port of the corresponding lalpulsar function, and store it in FstatMap.lnBtSG.

Returns:

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

Return type:

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