la_forge package

Submodules

la_forge.core module

class la_forge.core.Core(chaindir=None, corepath=None, burn=0.25, label=None, fancy_par_names=None, chain=None, params=None, pt_chains=False, skiprows=0)

Bases: object

An object that stores the parameters and chains from a bayesian analysis. Currently configured specifically for posteriors produced by PTMCMCSampler.

Parameters

chaindirstr

Directory with chains and file with parameter names. Currently supports chains as {‘chain_1.txt’,’chain.fit’} and parameters as {‘pars.txt’,’params.txt’,’pars.npy’} . If chains are stored in a FITS file it is assumed that the parameters are listed as the column names.

labelstr

Name of the core.

burnint, optional

Number of samples burned from beginning of chain. Used when calculating statistics and plotting histograms.

fancy_par_nameslist of str

List of strings provided as names to be used when plotting parameters. Must be the same length as the parameter list associated with the chains.

chainarray, optional

Array that contains samples from an MCMC chain that is samples x param in shape. Loaded from file if dir given as chaindir.

paramslist, optional

List of parameters that corresponds to the parameters in the chain. Is loaded automatically if in the chain directory given above.

corepathstr

Path to an already saved core. Assumed to be an hdf5 file made with la_forge.

pt_chainsbool

Whether to load all higher temperature chains from a parallel tempering (PT) analysis.

skiprowsint

Number of rows to skip while loading a chain text file. This effectively acts as a burn in, which can not be changed once the file is loaded (unless loade again). Useful when dealing with large chains and loading multiple times.

credint(param, onesided=False, interval=68)

Returns credible interval of parameter given.

Parameters

paramstr, list of str

onesidedbool, optional

Whether to calculate a one sided or two sided credible interval. The onesided option gives an upper limit.

interval: float, optional

Width of interval in percent. Default set to 68%.

get_map_dict()

Return a dictionary of the max a postori values for the parameters in the core. The keys are the appropriate parameter names.

get_map_param(param)

Returns maximum a posteri value of samples for the parameter given.

get_param(param, to_burn=True)

Returns array of samples for the parameter given.

param can either be a single list or list of strings.

get_param_credint(param, onesided=False, interval=68)

Returns credible interval of parameter given.

Parameters

paramstr, list of str

onesidedbool, optional

Whether to calculate a one sided or two sided credible interval. The onesided option gives an upper limit.

interval: float, optional

Width of interval in percent. Default set to 68%.

get_param_median(param)

Returns median of parameter given.

property map_idx

Maximum a posteri parameter values. From burned chain.

property map_params

Return all Maximum a posteri parameters.

median(param)

Returns median of parameter provided. Can be given as a string or list of strings.

save(filepath)

Save Core object as HDF5.

set_burn(burn)

Set number of samples to burn.

Parameters

burnint, float

An integer designating the number of samples to remove from beginning of chain. A float between 0 and 1 can also be used, which will estimate the integer value as a fraction of the full array length.

set_fancy_par_names(names_list)

Set fancy_par_names.

set_rn_freqs(freqs=None, Tspan=None, nfreqs=30, log=False, partimdir=None, psr=None, freq_path='fourier_components.txt')

Set gaussian process red noise frequency array.

Parameters

freqslist or array of floats

List or array of frequencies used for red noise gaussian process.

Tspanfloat, optional

Timespan of the data set. Used for calculating frequencies. Linear array is calculated as [1/Tspan, … ,nfreqs/Tspan].

nfreqsint, optional

Number of frequencies used for red noise gaussian process.

logbool, optional

Whether to use a log-linear space when calculating the frequency array.

partimdirstr, optional

Directory with pulsar data (assumed the same for tim and par files.) Calls the utils.get_Tspan() method which loads an enterprise.Pulsar(psr,partimdir) and extracts the timespan.

psrstr, optional

Puslar name, used when get the time span by loading enterprise.Pulsar() as in the documentation of partimdir above. It is assumed that there is only one par and tim file in the directory with this pulsar name in the file name.

freq_pathstr, optional

Path to a txt file containing the rednoise frequencies to be used.

Returns

Array of red noise frequencies.

class la_forge.core.HyperModelCore(label=None, param_dict=None, chaindir=None, burn=0.25, corepath=None, fancy_par_names=None, chain=None, params=None, pt_chains=False, skiprows=0)

Bases: Core

A class to make cores for the chains made by the enterprise_extensions HyperModel framework.

model_core(nmodel)

Return a core that only contains the parameters and samples from a single HyperModel model.

class la_forge.core.TimingCore(chaindir=None, burn=0.25, label=None, fancy_par_names=None, chain=None, params=None, pt_chains=False, tm_pars_path=None)

Bases: Core

A class for cores that use the enterprise_extensions timing framework. The Cores for timing objects need special attention because they are sampled in a standard format, rather than using the real parameter ranges. These Cores allow for automatic handling of the parameters.

get_param(param, to_burn=True, tm_convert=True)

Returns array of samples for the parameter given. Will convert timing parameters to physical units based on TimingCore.tm_pars_orig entries. Will also accept shortened timing model parameter names, like PX.

param can either be a single list or list of strings.

mass_function(PB, A1)

Computes Keplerian mass function, given projected size and orbital period.

Parameters

PBfloat

Orbital period [days]

A1float

Projected semimajor axis [lt-s]

Returns

mass function

Mass function [solar mass]

mass_pulsar()

Computes the companion mass from the Keplerian mass function. This function uses a Newton-Raphson method since the equation is transcendental.

la_forge.core.load_Core(filepath)

la_forge.diagnostics module

la_forge.gp module

class la_forge.gp.Signal_Reconstruction(psrs, pta, chain=None, burn=None, p_list='all', core=None)

Bases: object

Class for building Gaussian process realizations from enterprise models.

reconstruct_signal(gp_type='achrom_rn', det_signal=False, mlv=False, idx=None, condition=False, eps=1e-16)

Parameters

gp_typestr, {‘achrom_rn’,’gw’,’DM’,’none’,’all’,timing parameters}

Type of gaussian process signal to be reconstructed. In addition any GP in psr.fitpars or Signal_Reconstruction.gp_types may be called.

[‘achrom_rn’,’red_noise’] : Return the achromatic red noise.

[‘DM’] : Return the timing-model parts of dispersion model.

[timing parameters] : Any of the timing parameters from the linear timing model. A list is available as psr.fitpars.

[‘timing’] : Return the entire timing model.

[‘gw’] : Gravitational wave signal. Works with common process in full PTAs.

[‘none’] : Returns no Gaussian processes. Meant to be used for returning only a deterministic signal.

[‘all’] : Returns all Gaussian processes.

det_signalbool

Whether to include the deterministic signals in the reconstruction.

mlvbool

Whether to use the maximum likelihood value for the reconstruction.

idxint, optional

Index of the chain array to use.

Returns

wavearray

A reconstruction of a single gaussian process signal realization.

sample_params(index)

sample_posterior(samp_idx, array_params=['alphas', 'rho', 'nE'])

la_forge.rednoise module

la_forge.rednoise.determine_if_limit(vals, threshold=0.1, minval=-10, lower_q=0.3)

Function to determine if an array or list of values is sufficiently

separate from the minimum value.

Parameters

valsarray or list

threshold: float

Threshold above minval for determining whether to count as twosided interval.

minval: float

Minimum possible value for posterior.

lower_q: float

Percentile value to evaluate lower bound.

la_forge.rednoise.get_Tspan(pulsar, filepath=None, fourier_components=None, datadir=None)

Function for getting timespan of a set of pulsar dataself.

Parameters

pulsarstr

filepathstr

Filepath to a txt file with pulsar name and timespan in two columns. If supplied this file is used to return the timespan.

fourier_componentslist or array

Frequencies used in gaussian process modeling. If given 1/numpy.amin(fourier_components) is retruned as timespan.

datadirstr

Directory with pulsar data (assumed the same for tim and par files.) Calls the utils.get_Tspan() method which loads an enterprise.Pulsar() and extracts the timespan.

la_forge.rednoise.get_rn_freqs(core)

Get red noise frequency array from a core, with error message if noise array has not been included.

la_forge.rednoise.plot_free_spec(core, axis, parname_root, prior_min=None, ci=95, violin=False, Color='k', Fillstyle='full', plot_ul=False, verbose=True, Tspan=None)

Plots red noise free spectral parameters in units of residual time. Determines whether the posteriors should be considered as a fit a parameter or as upper limits of the given parameter and plots accordingly.

Parameters

corelist

la_forge.core.Core() object which contains the posteriors for the relevant red noise parameters to be plotted.

axismatplotlib.pyplot.Axis

Matplotlib.pyplot axis object to append various red noise parameter plots.

parname_rootstr

Name of red noise free spectral coefficient parameters.

prior_minfloat, None, ‘bayes’

Minimum value for uniform or log-uniform prior used in search over free spectral coefficients. If ‘bayes’ is used then a gorilla_bf calculation is done to determine if confidence interval should be plotted.

verbosebool, optional

Colorstr, optional

Color of the free spectral coefficient markers.

Fillstylestr, optional

Fillstyle for the free spectral coefficient markers.

Tspanfloat, optional

Timespan of the data set. Used for converting amplitudes to residual time. Calculated from lowest red noise frequency if not provided.

la_forge.rednoise.plot_powerlaw(core, axis, amp_par, gam_par, verbose=True, Color='k', Linestyle='-', n_realizations=0, Tspan=None, to_resid=True)

Plots a power law line from the given parmeters in units of residual time.

Parameters

corelist

la_forge.core.Core() object which contains the posteriors for the relevant red noise parameters to be plotted.

axismatplotlib.pyplot.Axis

Matplotlib.pyplot axis object to append various red noise parameter plots.

amp_parstr

Name of red noise powerlaw amplitude parameter.

gam_parstr

Name of red noise powerlaw spectral index parameter (gamma).

verbosebool, optional

n_realizationsint, optional

Number of realizations to plot.

Colorlist, optional

Color to make the plot.

Tspanfloat, optional

Timespan of the data set. Used for converting amplitudes to residual time. Calculated from lowest red noise frequency if not provided.

la_forge.rednoise.plot_rednoise_spectrum(pulsar, cores, show_figure=True, rn_types=None, plot_2d_hist=True, verbose=True, Tspan=None, title_suffix='', freq_yr=1, plotpath=None, cmap='gist_rainbow', n_plaw_realizations=0, n_tproc_realizations=1000, n_bplaw_realizations=100, Colors=None, bins=30, labels=None, legend=True, legend_loc=None, leg_alpha=1.0, Bbox_anchor=(0.5, -0.25, 1.0, 0.2), freq_xtra=None, free_spec_min=None, free_spec_ci=95, free_spec_violin=False, free_spec_ul=False, ncol=None, plot_density=None, plot_contours=None, add_2d_scatter=None, bplaw_kwargs={}, return_plot=False, excess_noise=False, levels=(0.39346934, 0.86466472, 0.988891))

Function to plot various red noise parameters in the same figure.

Parameters

pulsarstr

coreslist

List of la_forge.core.Core() objects which contain the posteriors for the relevant red noise parameters to be plotted.

Tspanfloat, optional

Timespan of the data set. Used for converting amplitudes to residual time. Calculated from lowest red noise frequency if not provided.

show_figurebool

rn_typeslist {‘’,’_dm_gp’,’_chrom_gp’,’_red_noise’}

List of strings to choose which type of red noise parameters are used in each of the plots.

plot_2d_histbool, optional

Whether to include two dimensional histogram of powerlaw red noise parameters.

verbosebool, optional

title_suffixstr, optional

Added to title of red noise plot as: ‘Red Noise Spectrum: ‘ + pulsar + ‘ ‘ + title_suffix

freq_yrint , optional

Number of 1/year harmonics to include in plot.

plotpathstr, optional

Path and file name to which plot will be saved.

cmapstr, optional

Color map from which to cycle plot colrs, if not given in Colors kwarg.

n_plaw_realizationsint, optional

Number of powerlaw realizations to plot.

n_tproc_realizationsint, optional

Number of T-process realizations to plot.

Colorslist, optional

List of colors to cycle through in plots.

labelslist, optional

Labels of various plots, for legend.

legend_loctuple or str, optional

Legend location with respect to Bbox_anchor.

leg_alphafloat, optional

Opacity of legend background.

Bbox_anchortuple, optional

This is the bbox_to_anchor value for the legend.

la_forge.rednoise.plot_tprocess(core, axis, alpha_parname_root, amp_par, gam_par, Color='k', n_realizations=100, Tspan=None)

Plots a power law line from the given parmeters in units of residual time.

Parameters

corelist

la_forge.core.Core() object which contains the posteriors for the relevant red noise parameters to be plotted.

axismatplotlib.pyplot.Axis

Matplotlib.pyplot axis object to append various red noise parameter plots.

alpha_parname_rootstr

Root of the t-process coefficient names, i.e. for J1713+0747_red_noise_alphas_0 give: ‘J1713+0747_red_noise_alphas’.

amp_parstr

Name of red noise powerlaw amplitude parameter.

gam_parstr

Name of red noise powerlaw spectral index parameter (gamma).

n_realizationsint, optional

Number of realizations to plot.

Colorlist, optional

Color to make the plot.

Tspanfloat, optional

Timespan of the data set. Used for converting amplitudes to residual time. Calculated from lowest red noise frequency if not provided.

la_forge.slices module

class la_forge.slices.SlicesCore(label=None, slicedirs=None, pars2pull=None, params=None, corepath=None, fancy_par_names=None, verbose=True, burn=0.25, parfile='pars.txt')

Bases: Core

A class to make a la_forge.core.Core object that contains a subset of parameters from different chains. Currently this supports a list of strings for multiple columns of a given txt file or a single string.

Parameters

la_forge.slices.get_col(col, filename)

la_forge.slices.get_idx(par, filename)

la_forge.slices.store_chains(filepaths, idxs, verbose=True)

la_forge.utils module

la_forge.utils.bayes_fac(samples, ntol=200, logAmin=-18, logAmax=-12, nsamples=100, smallest_dA=0.01, largest_dA=0.1)

Computes the Savage Dickey Bayes Factor and uncertainty. Based on code in enterprise_extensions. Expanded to include more options for when there are very few samples at lower amplitudes.

Parameters

• samples – MCMC samples of GWB (or common red noise) amplitude

• ntol – Tolerance on number of samples in bin

• logAmin – Minimum log amplitude being considered.

• logAmax – Maximum log amplitude being considered.

Returns

(bayes factor, 1-sigma bayes factor uncertainty)

la_forge.utils.compute_rho(log10_A, gamma, f, T)

Converts from power to residual RMS.

la_forge.utils.epoch_ave_resid(psr, correction=None, dt=10)

Epoch averaged residuals organized by receiver.

Parameters

psrenterprise.pulsar.Pulsar

correctionarray, optional

Numpy array which gives a correction to the residuals. Used for adding various Gaussian process realizations or timing model perturbations.

dtfloat

Coarse graining time [sec]. Sets filter size for TOAs.

Returns

fe_residsdict of arrays

Dictionary where each entry is an array of epoch averaged TOAS, residuals and TOA errors. Keys are the various receivers.

fe_maskdict of arrays

Dictionary where each entry is an array that asks as a mask for the receiver used as a key.

la_forge.utils.figsize(scale)

la_forge.utils.getMax2d(samples1, samples2, weights=None, smooth=True, bins=[40, 40], x_range=None, y_range=None, logx=False, logy=False, logz=False)

Function to return the maximum likelihood values by interpolating over a two dimensional histogram made of two sets of samples.

Parameters

samples1, samples2array or list

Arrays or lists from which to find two dimensional maximum likelihood values.

weightsarray of floats

Weights to use in histogram.

binslist of ints

List of 2 integers which dictates number of bins for samples1 and samples2.

x_rangetuple, optional

Range of samples1

y_rangetuple, optional

Range of samples2

logxbool, optional

A value of True use log10 scale for samples1.

logybool, optional

A value of True use log10 scale for samples2.

logzbool, optional

A value of True indicates that the z axis is in log10.

la_forge.utils.get_Tspan(pulsar, datadir)

Returns timespan of a pulsars dataset by loading the pulsar as an enterprise.Pulsar() object.

Parameters

pulsarstr

datadirstr

Directory where par and tim files are found.

la_forge.utils.get_params_2d_mlv(core, par1, par2)

Convenience function for finding two dimensional maximum likelihood value for any two parameters.

la_forge.utils.get_rn_noise_params_2d_mlv(core, pulsar)

Convenience function to find 2d rednoise maximum likelihood values.

la_forge.utils.powerlaw(freqs, log10_A=-16, gamma=5)

la_forge.utils.quantize_fast(toas, residuals, toaerrs, dt=0.1)

Function to quantize and average TOAs by observation epoch. Used especially for NANOGrav multiband data.

Based on [3].

Parameters

toasarray

residualsarray

toaerrsarray

dtfloat

Coarse graining time [sec].

la_forge.utils.rn_power(amp, gamma=None, freqs=None, T=None, sum_freqs=True)

Calculate the power in a red noise signal assuming the P=A^2(f/f_yr)^-gamma form.

la_forge.utils.set_publication_params(param_dict=None, scale=0.5)

la_forge.utils.weighted_quantile(values, quantiles, sample_weight=None, values_sorted=False, old_style=False)

Very close to numpy.percentile, but supports weights. NOTE: quantiles should be in [0, 1]! [From Max Ghenis via Stack Overflow: https://stackoverflow.com/questions/21844024/weighted-percentile-using-numpy]

Parameters

valuesnumpy.array

The data.

quantilesarray-like

Many quantiles needed.

sample_weightarray-like

Samples weights. The same length as array.

values_sortedbool

If True, then will avoid sorting of initial array.

old_style: bool

If True, will correct output to be consistent with numpy.percentile.

Returns

computed quantilesnumpy.array

Module contents

Top-level package for La Forge.