Classes and functions

beansp.beans.prior_func(theta_in)[source]

This function is the default prior and implements a simple box prior for all the parameters. Notably the systematic error on the burst times is limited to 10, which gives a maximum absolute error on the times of 10*BSTART_ERR = 100 min

Parameters:

theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, and (optionally) mass, radius & f_t

Returns:

prior probability

beansp.beans.prior_1808(theta_in)[source]

This function implements a simple box prior for all the parameters excluding mass and radius, which comes instead from Beans.mr_prior() and the metallicity, which comes instead from Beans.lnZprior()

This prior is explicitly intended for use with SAX J1808.4-3658, and should only be used with extreme caution in other cases!

Parameters:

theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, mass, radius, and (optionally) f_t

Returns:

prior probability

beansp.beans.lnZprior(z)[source]

This beta function for the metallicity prior is from Andy Casey and is an approximation of the metallicity of a mock galaxy at 2.5-4.5 kpc for the location of SAX J1808.4-3658. Assuming ZCNO = 0.01 is average value.

Parameters:

z – CNO metallicity (mass fraction)

Returns:

prior probability

beansp.beans.corr_goodwin19(burst, **kwargs)[source]

This is an example Settle correction function that applies the correction of Goodwin et al. (2019), multiplying the recurrence time and burst energy by 0.65. With this correction there is no need to modify the alpha value. The **kwargs can be used to incorporate other parameters into the correction; from settle we pass the base flux F, mdot M, H-fraction X, metallicity Z, and neutron star radius R and mass M (none used for this example)

Parameters:

burst – 3-element tuple output from settle, with alpha, trec [hr], and burst energy [1e39 erg], all values in the observer frame

Returns:

tuple with any necessary corrections performed

beansp.burstrain.generate_burst_train(bean, base, x_0, z, dist, xi_p, mass, radius, full_model=False, debug=False)[source]

This routine generates a simulated burst train based on the model input parameters, and the mdot history inferred from the persistent flux measurements (tobs, pflux, pfluxe)

Parameters:

  • bean – Beans object, from which the remaining parameters are drawn: bstart, pflux, pfluxe, tobs, numburstssim, ref_ind

  • base – base flux [MeV/nucleon]

  • x_0 – accreted H-fraction

  • z – accreted CNO metallicity

  • dist – source distance (kpc)

  • xi_p – anisotropy of persistent emission

  • mass – NS mass (M_sun)

  • radius – NS radius (km)

  • full_model – if set to True, include all the parameters in the dict that is returned

  • debug – set to True to show additional debugging information

Returns:

a dictionary with the following keys: ['time', 'mdot', 'alpha', 'e_b', 'iref'], and optionally (if full_model is set to True) also: ['base', 'z', 'x_0', 'dist', 'xi_p', 'mdot_max', 'mass', 'radius', 'forward', 'backward'] (all but the last 2 are copies of the input parameters).

The main outputs are the time, mdot, alpha, and e_b elements (numpy arrays) which are the predicted properties of the bursts and the accretion rate inferred for each interval.

We have one more element of the time array than the other arrays, because we can’t determine the properties for that burst, as we may not have enough data to back project. So the ith element of e_b, alpha etc. belongs with the (i+1)th element of time:

time  [ 0  1  2  3  4  5  6  7  ...  n ]
mdot     [ 0  1  2  3  4  5  6  ... n-1 ]
alpha    [ 0  1  2  3  4  5  6  ... n-1 ]
e_b      [ 0  1  2  3  4  5  6  ... n-1 ]
beansp.burstrain.burstensemble(bean, base, x_0, z, dist, xi_p, mass, radius, full_model=False)[source]

This routine generates as many burst predictions as there are burst measurements. Written initially by Luke Waterson, 2021

Parameters:

  • bean – Beans object, from which the remaining parameters are drawn: bstart, pflux, numburstsobs

  • base – base flux [MeV/nucleon]

  • z – accreted CNO metallicity

  • x_0 – accreted H-fraction

  • dist – source distance (kpc)

  • xi_p – anisotropy of persistent emission

  • mass – NS mass (M_sun)

  • radius – NS radius (km)

Returns:

a dictionary with the following keys: ['base', 'z', 'x_0', 'dist', 'xi_p', 'time', 'mdot_max', 'mdot', 'iref', 'alpha', 'e_b', 'mass', 'radius']

The main outputs (as for burstrain.generate_burst_train()) are the time, e_b, alpha, and mdot elements (numpy arrays) which are the predicted properties of the bursts and the accretion rate inferred for each interval.

class beansp.Beans(prior=<function prior_func>, lnlike=<function lnlike_sys>, corr=None, config_file=None, run_id='test', obsname=None, burstname=None, gtiname=None, continuous=True, maxgap=2, interp='linear', smooth=0.02, model=<function settle>, theta=(0.58, 0.013, 0.4, 3.5, 1.0, 1.0, 1.5, 11.8), sampler='emcee', fluen=True, alpha=False, pflux=True, numburstssim=3, bc=2.21, ref_ind=1, threads=4, test_model=True, restart=False, **kwargs)[source]

The main object class that includes the basic functionality required for beans. The code will read in burst (and observation) data and attempt to simulate bursts to match the observed burst properties. There are two principle modes; the original function generates a “train” of individual bursts, observed (for example) during a transient outburst, as for the original application to the 2002 outburst of SAX J1808.4-3658, observed with RXTE/PCA. The alternative is to match to a set of non-contiguous bursts (“ensemble” mode)

__init__(prior=<function prior_func>, lnlike=<function lnlike_sys>, corr=None, config_file=None, run_id='test', obsname=None, burstname=None, gtiname=None, continuous=True, maxgap=2, interp='linear', smooth=0.02, model=<function settle>, theta=(0.58, 0.013, 0.4, 3.5, 1.0, 1.0, 1.5, 11.8), sampler='emcee', fluen=True, alpha=False, pflux=True, numburstssim=3, bc=2.21, ref_ind=1, threads=4, test_model=True, restart=False, **kwargs)[source]

Initialise a Beans object. In addition to the parameters below, the user can specify additional parameter relevant to the various samplers: nwalkers, nsteps, stretch_a for emcee; and nlive, dlogz for dynesty. If these are not specified default values will be applied (when Beans.do_run() is called).

Parameters:

  • prior – prior function to use

  • corr – correction function for bursts, or None

  • config_file – file to read in configuration parameters (in which case the keyword params below are ignored)

  • run_id – string identifier for the run, used to label all the result files, and where you want output to be saved

  • obsname – name of the file including the flux history, from which the mdot is estimated to generate to generate the burst train (set obsname=None for a non-contiguous, or “ensemble” mode run)

  • burstname – name of the burst data file, listing the bursts

  • gtiname – name of the GTI file, set to None to turn off checking

  • continuous – for burst “train” modes, set to True to generate a continuous train (i.e. using generate_burst_train rather than the newer punkt_train)

  • interp – interpolation mode for the flux; possible values are ‘linear’, or ‘spline’

  • smooth – smoothing factor for spline interpolation

  • model – burst model to use, one of settle or grid_interp

  • theta – initial centroid values for walker model parameters, with X, Z, Q_b, d, xi_b, xi_p, and (optionally) mass, radius, & f_t

  • fluen – set to True (default) to include the fluences in the data for comparison, or False to omit

  • alpha – set to True to include the alphas in the data for comparison, or False (default) to omit

  • pflux – set to True (default) to include the peak fluxes in the data for comparison, if they have been supplied

  • numburstssim – number of bursts to simulate, for the “train” mode, both earlier and later than the reference burst; i.e. set to half the total number of bursts you want to simulate. Don’t forget to account for missed bursts!

  • bc – bolometric correction to adopt for the flux history (set to 1.0 if the fluxes are already bolometric):

  • ref_ind – rank of “index” burst, against which the other burst times are relative to. For the “train” mode, should be around the middle of the predicted burst train. This burst will not be simulated but will be used as a reference to predict the other bursts.

  • threads – number of threads for emcee to use (e.g. number of cores your computer has). Set to None to use all available

  • test_model – flag to test the model during the setup process

  • restart – set to True to continue a previously interrupted run

Result:

Beans object including all the required data

burst_table(show=True, nsep=None, predicted=False, key=None)[source]

This method creates and, optionally, displays a table from the model and observed bursts.

The table is returned as an astropy Table object that can be saved as an MRT file, like so: tab.write('test.dat',format='mrt') although you will have to then edit the file afterwards to complete the metadata

Parameters:

  • show – set to True to display each row (in LaTeX format)

  • nsep – number of inferred bursts between each pair, for use when model data is unavailable and the automatic code doesn’t return the desired values

  • predicted – set to True to show predicted bursts also (not yet implemented)

  • key – for the case where we have more than one set of predictions (e.g. for different numbers of predicted bursts), the key will identify which to return

Returns:

astropy table of burst properties

do_analysis(options=['autocor', 'posteriors'], part=None, truths=False, burnin=2000, ilab=None, title=None, savefig=False)[source]

This method is for running standard analysis and displaying the results. Nothing is returned, but various options for analysis are available, which will (optionally) populate some of the beansp.Beans attributes, including

reader - sampler object read in from HDF file
sampler - chain data
nsteps_completed - number of steps completed
samples - flattened samples, omitting first samples_burnin
last - last walker positions
probs - likelihoods (total, prior, and contributions from each parameter) for the last walker positions
cc - ChainConsumer object with samples and derived cosi, gravity, redshift
cc_parameters - plot labels for cc object
model_pred - dictionary with model realisations read in from the “blobs”

By default the method will also create several files, labeled by the run_id; drawn from

{}_autocorrelationtimes.pdf (via plot_autocorr)
{}_posteriors.pdf
{}_parameterconstraints_pred.txt

TODO: need to reorganise a bit, and add more options

Parameters:

  • options – array of strings corresponding to various analysis options, listed in the analyses dict below

  • part – string or array “partition” dividing the set of samples into two or more separate groups for analysis

  • truths – parameter vector to overplot on (one of the) corner plots (TODO: need to check if >1 corner plot are selected +truths, which will likely result in an error due to incompatible number of parameters)

  • burnin – number of steps to discard when plotting the posteriors. If -ve, discard all but that number

  • title – set to a string to add a title to the plot, or False to omit; if set to None will print some generic information

  • ilab – set to a list of (integer) inclination values to label in the anisotropy model plot (‘fig8’)

  • savefig – set to True to save figures to .pdf files, False to skip

Returns:

none

do_run(sampler=None, plot=False, analyse=True, burnin=2000, **kwargs)[source]

This routine runs the chain for as many steps as is specified in the init call. Emcee parameters are defined in runemcee module. Have previously used multiprocessing to speed things up, but that’s not currently active

Parameters:

  • sampler – define the sampler to use; default is emcee

  • plot – set to True to do the plot of the initial guess. This seems redundant since it’s also plotted at the __init__ stage

  • analyse – set to True to call do_analysis automatically once the chains finish

  • burnin – number of steps to ignore at the start of the run, passed to do_analysis

Returns:

lnprob(theta_in, x, y, yerr)[source]

The full log-probability function incorporating the priors (via the Beans.lnprior attribute), and model likelihood (via Beans.lnlike()), that is passed to runemcee when creating the sampler (in the Beans.do_run() method).

Calls self.lnprior and self.lnlike; these are pointers to functions that can be modified, and the theta_in is not unpacked here, so you can add your own combinations to customise the code

Parameters:

  • theta_in – parameter vector, with X, Z, Q_b, d, xi_b, xi_p, mass, radius, and (optionally) f_t

  • x – the “independent” variable, passed to Beans.lnlike()

  • y – the “dependent” variable (i.e. measurements), passed to Beans.lnlike()

  • yerr – erorr estimates on y

Returns:

total (prior+model) likelihood, prior likelihood, model array (from Beans.lnlike())

plot(show_model=True, model=None, mdot=True, imatch=None, title=None, savefig=False)[source]

Display a plot of the data and model results, for a burst train calculated with burstrain.generate_burst_train(); or burst “ensemble” data calculated with burstrain.burstensemble(). Adapted from the example at https://matplotlib.org/gallery/api/two_scales.html

Parameters:

  • show_model – set to False to skip the model generation step, in which case only the observed data will be plotted

  • model – array of packed model prediction, OR dict giving full model results

  • mdot – flag to show mdot rather than flux (only possible if you’re passing the full model)

  • imatch – can pass a “matching” array for train models, to use instead of the automatically-calculated version

  • title – add a title, if required

  • savefig – set to True to save the figure to PDF

prune(key=None, nwalkers=None, scale=0.0, add_param=None, savefile=None)[source]

This method will “prune” the walkers to keep only one set of solutions. It is necessary that Beans.do_analysis() has already been called, with the comparison option, to identify the model classes.

The method will return the new set of positions, and optionally save them to a pickle file

Parameters:

  • key – model class to keep, usually labeled by the number of predicted bursts

  • nwalkers – number of samples to generate, if not the current number of walkers

  • scale – in case it’s necessary to distribute the walker positions around the seed values, this parameter sets the Gaussian scale

  • add_param – number of additional parameters to add in the saved array

  • savefile – name of pickle file to save to

Returns:

array of walker positions, dimensions (nwalkers, ndim)

sim_data(file=None, ensemble=False)[source]

This method generates a table of simulated data based on the current theta vector, that can be used for tests of the MCMC method

Errors are taken as the mean of the observational errors; bursts matching observed bursts (according to the matching algorithm) are flagged as True in the 6th column

In lieu of an output file, you can just paste the results into a text file and use as the input (specify as burstname)

Parameters:

  • file – file to save the results to

  • ensemble – set to True to output “ensemble” data, i.e. including persistent fluxes, and with recurrence times rather than burst times

Returns: