mbproj2d API reference

Detailed below are the various classes which make up mbproj2d. Note that all the submodules are imported into the main mbproj2d module, so they don’t need to be imported separately.

Models

class mbproj2d.model.BackModelBase(name, pars, images, expmap=None)

Base background model. Does nothing.

Parameters:

• name – name of model

• pars – Pars() object

• images – list of Image() objects

• expmap – exposure map name to lookup in Image

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

class mbproj2d.model.BackModelFlat(name, pars, images, log=False, normarea=False, defval=0.0, expmap=None)

Flat background model.

Parameters:

• name – name of model

• pars – dict of parameters

• images – list of data.Image objects

• log (bool) – apply log scaling to value of background

• normarea (bool) – normalise background to per sq arcsec

• defval – default parameter

• expmap – name or index of exposure map to use (if any)

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

class mbproj2d.model.BackModelVigNoVig(name, pars, images, defval=0.0, expmap='expmap', expmap_novig='expmap_novig', rescale_novig=True)

A background model with vignetted and non-vignetted components.

The X_vf holds the vignetted fraction. This is an unbounded parameter which is transformed to the vignetting fraction via the sigmoid function. By default it has a Gaussian prior about 0.

All parameters are log units per square arcsec

Note that the unvignetted exposure map is multiplied by median of the ratio of the vignetted to unvignetted if rescale_novig is True. This helps to make a seamless transition between the two as a function of X_vf. This reduces the degeneracies between the X_vf and normalisation parameters.

Parameters:

• name – name of model

• pars – dict of parameters

• images – list of data.Image objects

• defval – default parameter

• expmap – name or index of vignetted exposure map to use

• expmap_novig – name or index of unvignetted exposure map to use

• rescale_vig – rescale non-vignetted map by ratio of vig to non-vig

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

estimatePars(pars, images, binup=8)

Take images and estimate initial background values.

exception mbproj2d.model.ModelEvaluationException

class mbproj2d.model.SrcModelBase(name, pars, images, cx=0.0, cy=0.0)

Base class for source models.

Parameters:

• name – name of model

• pars – dict of parameters

• images – list of data.Image objects

• cx – initial centre x coordinate

• cy – initial centre y coordinate

class mbproj2d.model.TotalModel(pars, images, src_models=None, src_expmap=None, back_models=None)

Combined model for data.

Parameters:

• pars – Pars object (currently unused here)

• images – list of Image objects

• src_models – list of source models

• src_expmap – name of exposure map to use for sources

• back_models – list of background models

compute(pars, apply_psf=True, apply_expmap=True, apply_src=True, apply_back=True)

Return a list of image arrays for the models. :param pars: Pars() object with parameters :param apply_psf: whether to convolve with PSF :param apply_expmap: whether to multiply by exposure map :param apply_src: apply source models :param apply_back: apply background models

compute_separate(pars, apply_psf=True, apply_expmap=True)

Compute model for each component separately and combined. :param pars: Pars() object with parameters :param apply_psf: whether to convolve with PSF :param apply_expmap: whether to multiply by exposure map

prior(pars)

Given parameters, compute prior.

This includes the contribution from parameters, the object models and background models.

writeToFits(filename, pars, mask=True, apply_psf=True, apply_expmap=True, apply_src=True, apply_back=True, trim_original=False)

Write model as a series of HDUs in a FITS file.

Parameters:

• filename – FITS filename

• pars – pars dictionary

• mask – mask out masked pixels

• apply_psf – convolve with PSF

• apply_expmap – multiply by exposure map

• apply_src – include source models

• apply_back – include background models

• trim_original – write to original size

Cluster modelling

class mbproj2d.cluster.ClusterBase(name, pars, images, cosmo=None, NH_1022pcm2=0.0, maxradius_kpc=3000.0, cx=0.0, cy=0.0)

Base cluster model class.

This does quite a lot of work for the child classes, as it takes the computed ne,T,Z profiles and makes the images.

compute(pars, imgarrs)

Add cluster model to images.

computeMassProfile(pars, radii)

Return g and Phi profiles for cluster, if available.

computeProfiles(pars, radii)

Compute plasma profiles for use in physical profile computation

Parameters:

• pars – Pars object with parameters

• radii – Radii object with radii to compute for

Returns (ne_prof, T_prof, Z_prof)

class mbproj2d.cluster.ClusterHydro(name, pars, images, cosmo=None, NH_1022pcm2=0.0, ne_prof=None, mass_prof=None, Z_prof=None, gas_has_mass=True, maxradius_kpc=3000.0, cx=0.0, cy=0.0)

Hydrostatic model for cluster, given density, mass model and metallicity profile.

computeMassProfile(pars, radii)

Compute g and potential given parameters.

computeProfiles(pars, radii)

Compute plasma profiles assuming hydrostatic equilibrium

Parameters:

• pars – Pars object with parameters

• radii – Radii object with radii to compute for

Returns (ne_prof, T_prof, Z_prof)

class mbproj2d.cluster.ClusterHydroPressureMass(name, pars, images, cosmo=None, NH_1022pcm2=0.0, p_prof=None, mass_prof=None, Z_prof=None, maxradius_kpc=3000.0, cx=0.0, cy=0.0)

Hydrostatic model for cluster, given pressure and mass model and metallicity profile.

computeMassProfile(pars, radii)

Compute g and potential given parameters.

computeProfiles(pars, radii)

Compute plasma profiles assuming hydrostatic equilibrium

Parameters:

• pars – Pars object with parameters

• radii – Radii object with radii to compute for

Returns (ne_prof, T_prof, Z_prof)

class mbproj2d.cluster.ClusterNonHydro(name, pars, images, cosmo=None, NH_1022pcm2=0.0, ne_prof=None, T_prof=None, Z_prof=None, maxradius_kpc=3000.0, cx=0.0, cy=0.0)

Model for a cluster, given density, temperature and metallicity profiles.

computeProfiles(pars, radii)

Compute plasma profiles.

Parameters:

• pars – Pars object with parameters

• radii – Radii object with radii to compute for

Returns (ne_prof, T_prof, Z_prof)

class mbproj2d.cluster.EmissionMeasureCluster(name, pars, images, cosmo=None, NH_1022pcm2=0.0, emiss_prof=None, T_prof=None, Z_prof=None, maxradius_kpc=3000.0, cx=0.0, cy=0.0)

Less physical cluster model, where we parameterize using a projected emission measure profile.

Here T_prof and Z_prof are also projected quantities emiss_prof has units of cm^-5

compute(pars, imgarrs)

Compute cluster images.

This completely overrides the base class.

mbproj2d.cluster.computeGasAccn(radii, ne_prof)

Compute acceleration due to gas mass for density profile given.

Soft X-ray background modelling

class mbproj2d.sxbkg.BackModelImage(name, pars, images, backimgarrs, expmap=None)

Background model based on images.

Parameters:

• name – name of model

• pars – dict of parameters

• images – list of data.Image objects

• backimgarrs – images to use as background for each Image

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

class mbproj2d.sxbkg.BackModelXspecModel(name, pars, images, xcms, scale0, expmap=None, usearf=True)

Background model based on an xspec model(s) rate and an exposure map.

Parameters:

• name – model name

• images – list of Image objects

• xcms – list/tuple of xcm filenames

• scale0 – scaling value to apply to rates computed from xcm model

• expmap – exposure map to use (if any)

• usearf – use ARF in computation of rates

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

class mbproj2d.sxbkg.Cxb(name, pars, images, expmap=None, NH_1022pcm2=0, gamma=1.41, norm_scale=2.7e-10)

Model for calculating unresolved cosmic X-ray background

gamma = 1.41 (De Luca & Molendi 2004) norm_scale=11.6 photon/cm2/s/sr

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

class mbproj2d.sxbkg.Sxbkg(name, pars, images, expmap=None, NH_1022pcm2=0, T_unabsorbed_keV=0.1, T_absorbed_keV=0.28, scale_unabsorbed=1.6e-09, scale_absorbed=6.81e-11)

Model for calculating soft X-ray background

Units are scaled according to some reasonable values (norms taken from UDS field)

compute(pars, imgarrs, apply_expmap=True)

Compute the background model given parameters and input images. If apply_expmap is True, the expmap should be applied to the output.

Point source modelling

class mbproj2d.point.PointBase(name, pars, images, NH_1022pcm2=0.0, cx=0.0, cy=0.0)

class mbproj2d.point.PointPowerlaw(name, pars, images, NH_1022pcm2=0.0, cx=0.0, cy=0.0, gamma=1.7)

Powerlaw model at a particular position.

gamma is a fitting parameter (default fixed)

mbproj2d.point.addPointToImage(imgarr, cy, cx, rate)

Add point to image array, splitting across pixels as appropriate.

Model parameters

class mbproj2d.par.Par(val, prior=None, frozen=False, xform=None, linked=None, minval=-inf, maxval=inf, soft=False)

Parameter for model.

Parameters:

• val (float) – parameter value

• prior – prior object or None for flat prior

• frozen – whether to leave parameter frozen

• xform – function to transform value for model or ‘exp’ for an exp(x) scaling

• linked – another Par object to link this parameter to another

• minval (float) – minimum value for default flat prior

• maxval (float) – maximum value for default flat prior

• soft – use a soft flat prior instead of a sharp one

calcPrior()

Calculate prior.

isFree()

Is the parameter free?

property v

Value for using in model, after transformation or linking, if any.

class mbproj2d.par.Pars

Parameters for a model.

This is based around a dictionary class. Each parameter has a name.

bounds()

Return lower,upper bounds for free parameters.

calcPrior()

Return total prior of parameters.

copy()

Return a deep copy of self.

freeKeys()

Return sorted list of keys of parameters which are free

freeVals()

Return list of values for parameters which are free in sorted key order.

load(filename, skip=False)

Load parameters from file.

Parameters:

• filename – filename to load from

• skip – if set, then we continue if file not found

match(pattern, use_re=False)

Returns a dictionary of parameters whose names match a pattern.

Parameters:

• pattern – glob-style parameter match, e.g. “ne_*” or “abc_???_alpha” (default), a regular expression string (if use_re)

• use_re – if set, treat pattern as a regular expression

Returns {‘name’: par, …}

matchFreeze(pattern, use_re=False)

Freeze parameters which match the name given.

Parameters:

• pattern – glob-style parameter match, e.g. “ne_*” or “abc_???_alpha”

• use_re – if set, treat pattern as a regular expression

matchSet(pattern, val, use_re=False)

Set values for parameters which match the name given.

Parameters:

• pattern – glob-style parameter match, e.g. “ne_*” or “abc_???_alpha”.

• val – constant (to set to same value) or iterable (to set to sequence)

matchThaw(pattern, use_re=False)

Thaw parameters which match the name given.

Parameters:

• pattern – glob-style parameter match, e.g. “ne_*” or “abc_???_alpha”.

• use_re – if set, treat pattern as a regular expression

numFree()

Return number of free parameters

save(filename)

Saves the parameters as a Python pickle.

Note: upgrading the source code of mbproj2d or your prior may prevent the saved file from being loadable again. Take care before relying on this for long term storage.

Parameters:

filename – output filename

setFree(vals)

Given a list of values, set those which are free.

Note: number of free parameters should be number of vals

write(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print out parameters.

class mbproj2d.par.PriorBase

Base class for all Priors.

bounds()

Return upper and lower bounds.

makeValidValue(val)

Transform value to pass to model.

This can, for example, ensure that a model parameter is valid, even if the bounds are soft

paramFromUnit(unit)

Compute a parameter value to an input 0…1.

class mbproj2d.par.PriorBoundedGaussian(mu, sigma, minval=None, maxval=None)

Gaussian prior Bounded

Parameters:

• mu – Gaussian centre

• sigma – Gaussian width

• minval – minimum allowed value

• maxval – maximum allowed value

bounds()

Return upper and lower bounds.

paramFromUnit(unit)

Compute a parameter value to an input 0…1.

class mbproj2d.par.PriorBoundedGaussianSoft(mu, sigma, minval=None, maxval=None, width=0.01)

Gaussian prior Bounded with soft limits

Parameters:

• mu – Gaussian centre

• sigma – Gaussian width

• minval – minimum allowed value

• maxval – maximum allowed value

makeValidValue(val)

Transform value to pass to model.

This can, for example, ensure that a model parameter is valid, even if the bounds are soft

paramFromUnit(unit)

Compute a parameter value to an input 0…1.

class mbproj2d.par.PriorFlat(minval, maxval)

Flat prior.

Parameters:

• minval – minimum allowed value

• minval – maximum allowed value

bounds()

Return upper and lower bounds.

paramFromUnit(unit)

Compute a parameter value to an input 0…1.

class mbproj2d.par.PriorFlatSoft(minval, maxval, width=0.01)

Flat prior, where the likelihood is increase sharply increased at the bounds.

The prior is exponential in width beyond the edges

Parameters:

• minval – soft minimum value

• maxval – soft maximum value

• width – width of transition beyond minval/maxval

makeValidValue(val)

Transform value to pass to model.

This can, for example, ensure that a model parameter is valid, even if the bounds are soft

class mbproj2d.par.PriorGaussian(mu, sigma)

Gaussian prior

Parameters:

• mu – Gaussian centre

• sigma – Gaussian width

paramFromUnit(unit)

Compute a parameter value to an input 0…1.

Profiles

class mbproj2d.profile.ProfileBeta(name, pars)

Beta model.

Parameterised by: logn0: log_e of n0 beta: beta value logrc: log_e of rc_kpc.

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileBinned(name, pars, rbin_edges_kpc, defval=0.0, log=False)

Profile made up of constant values between particular radial edges.

Parameters:

• rbin_edges_kpc – array of bin edges, kpc

• defval – default value

• log – where to apply exp to output.

compute(radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileFlat(name, pars, defval=0.0, log=False, minval=-inf, maxval=inf)

Constant value profile.

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileGNFWDensity(name, pars)

Density model following the formalism of Eq.A1 in Nagai+07, which is more flexible than S&A+07

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileGNFWPressure(name, pars)

Pressure profile with Nagai+07 formalism, and is used in UPP (Arnaud+10)

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileInterpol(name, pars, rcent_kpc, defval=0.0, log=False, extrapolate=False, mode='linear')

Create interpolated profile between fixed values

Parameters:

• rcent_kpc – where to interpolate between in kpc (sorted)

• defval – initial parameter values

• log – log profile

• extrapolate – extrapolate profile beyond endpoints (for linear)

• mode – interpolation type (‘linear’, ‘cubic’, ‘quadratic’…)

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileMcDonaldT(name, pars)

Temperature model from McDonald+14, equation 1

Log values are are log_e

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.ProfileSum(name, pars, subprofiles)

Add one or more profiles to make a total profile.

compute(pars, radii)

Compute profile at centres of bins, given edges.

prior(pars)

Return any prior associated with this profile.

class mbproj2d.profile.ProfileVikhDensity(name, pars, mode='double')

Density model from Vikhlinin+06, Eqn 3.

Modes: ‘double’: all components ‘single’: only first component ‘betacore’: only first two terms of 1st cmpt (beta, with powerlaw core)

Densities and radii are are log base 10

compute(pars, radii)

Compute profile at centres of bins, given edges.

class mbproj2d.profile.Radii(rshell_kpc, num)

Define an equally-spaced set of shells/annuli and project between the two. :param rshell_kpc: width of shells/annuli in kpc :param num: number of shells/annuli

project(emissivity_pkpc3)

Project from emissivity profile to surface brightness (per kpc2).

Mass profiles

Define mass profiles.

class mbproj2d.profile_mass.CmptMassKing(name, cosmo, pars)

King potential.

This is the modified Hubble potential, where rho = rho0 / (1+(r/r0)**2)**1.5

r0 = sqrt(9*sigma**2/(4*Pi*G*rho0))

We define rho in terms of r0 and sigma

Model parameters are king_sigma_logkmps (sigma in log10 km/s) and king_rcore_logkpc (rcore in log10 kpc).

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassBase(name, cosmo, pars)

Base mass profile class.

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassGNFW(name, cosmo, pars)

Generalised NFW.

This is an NFW with a free inner slope (alpha).

rho(r) = rho0 / ( (r/rs)**alpha * (1+r/rs)**(3-alpha) )

For details see Schmidt & Allen (2007) http://adsabs.harvard.edu/doi/10.1111/j.1365-2966.2007.11928.x

Model parameters are gnfw_logconc (log10 concentration), gnfw_r200_logMpc (log10 r200 in Mpc) and gnfw_alpha (alpha parameter; 1 is standard NFW).

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassNFW(name, cosmo, pars)

NFW profile. Useful detals here: http://nedwww.ipac.caltech.edu/level5/Sept04/Brainerd/Brainerd5.html and Lisa Voigt’s thesis

Model parameters are

X_logconc (log_e concentration) and X_r200_logMpc (log_e r200 in Mpc)

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassNFWfnM(name, cosmo, pars, overdensity=500)

NFW profile which is a function of mass at some overdensity

Useful detals here: http://nedwww.ipac.caltech.edu/level5/Sept04/Brainerd/Brainerd5.html and Lisa Voigt’s thesis

This version parameterizes M (for some given overdensity) instead of R200

Model parameters are

X_logconc (log_e concentration) and X_M_logMsun (log_e Msun at overdensity)

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassPoint(name, cosmo, pars)

Point mass.

Model parameter is pt_M_logMsun, which is point mass in log solar masses.

compute(pars, radii)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

class mbproj2d.profile_mass.ProfileMassSum(name, cosmo, pars, cmpts)

Add together different profiles.

compute(pars)

Get gravitational field and potential profiles

Parameters:

• pars – Pars() parameters objects()

• radii – Radii() object

Returns: g, Phi

Data

class mbproj2d.data.Image(img_id, imagearr, emin_keV=0.5, emax_keV=2.0, rmf='image.rmf', arf='image.arf', pixsize_as=1.0, expmaps=None, mask=None, psf=None, origin=(0, 0), wcs=None, optimal_size=True, pad=0)

Image class details images to fit.

Parameters:

• img_id – unique id for image (str or int)

• imagearr – numpy image array for image

• emin_keV (float) – minimum energy

• emax_keV (float) – maximum energy

• rmf – response matrix file

• arf – ancillary response matrix file

• pixsize_as – size of pixels in arcsec

• expmaps – dict of numpy exposure map arrays (different components can use different exposure maps, if needed)

• mask – numpy mask array (None means no mask)

• psf – PSF object

• origin – position in pixels (y,x) coordinates are measured relative to (should be same position in all images)

• wcs – optional WCS stored with this image

• optimal_size – expand images to be optimal size for PSF convolution

• pad – pad image sizes by size given to ensure convolution doesn’t wrap

binUp(factor)

Return binned copy of image.

Parameters:

factor – integer bin factor

normalizeExposurePosn(ra, dec, rad_arcsec, exposure, expmaps=None)

Normalize exposure maps to have the mean exposure over the region.

The data are scaled so that the pixels within the circle (ra, dec) within the radius have the exposure given.

Requires a WCS to set on the input image. The mask is used to remove pixels when averaging.

Parameters:

• ra – ra (deg)

• dec – dec (deg)

• rad_arcsec – radius (arcsec)

• exposure – exposure (s)

• expmaps – which exposure maps to process: list or by default, all

class mbproj2d.data.PSF(img, pixsize_as=1.0, origin=None)

PSF modelling class.

Parameters:

• img – 2D PSF image

• pixsize_as – size of pixels in arcsec

• origin – (y, x) origin of PSF centre

applyTo(inimg, minval=1e-10)

Convolve image with PSF.

matchImage(imgshape, img_pixsize_as)

Adjust PSF to different pixel size and image shape.

Cosmology

Routines for calculating distances from cosmology.

TODO: Replace with astopy’s versions?

class mbproj2d.cosmo.Cosmology(z, H0=70.0, q0=0.5, WM=0.3, WV=0.7)

Cosmology calculation object.

Parameters:

• H0 (float) – Hubble’s constant (km/s/Mpc)

• q0 (float) – Deceleration parameter

• WM (float) – Matter density

• WV (float) – Vacuum density

property D_A

Get angular diameter distance in Mpc.

property D_L

Get luminosity distance in Mpc.

property as_kpc

Get number of kpc per arcsec.

property rho_c

Get critical density of universe (cgs).

Fitting

class mbproj2d.fit.Fit(images, totalmodel, pars)

For fitting a total model to the data.

like()

Get Likelihood object for current fit.

printHeader()

Show line with list of thawed/non-linked params.

printPars(like)

Print thawed/non-linked parameters to output on single line.

run(verbose=True, sigdelta=0.01, maxloops=10, maxfititer=None, methods=None, fitoptions=None)

Parameters:

• verbose – show fit progress

• sigdelta – what is a significant change in fit statistic

• maxfititer – maximum number of iterations to pass to fit routine

• maxloops – fit a maximum number of times with improvement before exiting

• methods – iterate through these fitting methods to look for improvement

• fitoptions – options to pass to minimizer

Default methods are

(‘LN_SBPLX’, ‘LN_BOBYQA’): if nlopt is installed (‘Nelder-Mead’, ‘Powell’): if nlopt is not installed

Returns (fit_like, success) where success indicates less than maximum number of iterations were done.

save(filename)

Saves the current fit state as a Python pickle.

This includes the input data, parameters and model.

Note: upgrading the source code of mbproj2d or your custom model may prevent the saved file from being loadable again. Take care before relying on this for long term storage.

Parameters:

filename – output filename

class mbproj2d.fit.Likelihood(images, model, pars)

Calculate a likelihood for the model and data.

Members: - prior - prior value - images - list of likelihoods for each input image - total - combined likelihood

Markov Chain Monte Carlo

class mbproj2d.mcmc.MCMC(fit, nwalkers=50, processes=1, initspread=0.01, moves=None, verbose=True, sampler='emcee')

Handles MCMC analysis of Fit. Please note this is deprecated for new code - please see MCMCSampler/MCMCSamplerEmcee and MCMCStore/MCMCStoreHDF5.

Parameters:

• fit (Fit) – Fit object to use for mcmc

• walkers (int) – number of walkers to use

• processes (int) – number of simultaneous processes to compute likelihoods

• initspread (float) – random Gaussian width added to create initial parameters

• moves – moves parameter to emcee.EnsembleSampler

• verbose – whether to write progress

• sampler – ‘emcee’ or ‘zeus’

burnIn(length, autorefit=True, minfrac=0.2, minimprove=0.01)

Burn in, restarting fit and burn if necessary.

Parameters:

• autorefit (bool) – refit position if new minimum is found during burn in

• minfrac (float) – minimum fraction of burn in to do if new minimum found

• minimprove (float) – minimum improvement in fit statistic to do a new fit

run(length, progress=True)

Run chain.

Parameters:

length (int) – length of chain

save(outfilename, thin=1, discard=0)

Save chain to HDF5 file.

Parameters:

• outfilename (str) – output hdf5 filename

• thin (int) – save every N samples from chain

• discard (int) – discard first N samples

verbprint(*args, **argsv)

Output if verbose.

Physical quantities

class mbproj2d.phys.Phys(pars, model, rmin_kpc=0.5, rmax_kpc=2000, rsteps=256, linbin_kpc=0.5, binning='log', average='midpt', fluxrange_keV=((0.5, 2.0),), luminrange_keV=((0.5, 2.0),), rate_rmf=None, rate_arf=None, rate_bands=((0.3, 2.3), (0.5, 2.0)))

Calculate physical profiles for a cluster model.

Parameters:

• pars – parameters to apply

• model – Cluster model component (not TotalModel!)

• rmin_kpc – minimum radius to compute profiles for

• rmax_kpc – maximum radius to compute profiles for

• rsteps – number of steps to compute profiles for

• binning – should be ‘log’ or ‘lin’ for bin size

• average – should be ‘midpt’, ‘volume’, ‘mean’ for how to convert from input to output bins

• linbin_kpc – internal linear bin size to use before rebinning

• fluxrange_keV – list of energy ranges to compute (unabsorbed) fluxes between

• luminrange_keV – list of rest frame energy ranges to compute luminosities between

• rate_rmf – RMF file to calculate count rates (if required)

• rate_arf – ARF file to calculate count rates (if required)

• rate_bands – Energy ranges to calculate rates within

calc(ne_pcm3, T_keV, Z_solar, g_cmps2, phi_ergpg)

Given input profiles, calculate output profiles.

These are assumed to have the output radii spacing

calcPhysChainStats(physchains, confint=68.269)

Take physical chains and compute profiles with 1 sigma ranges.

Parameters:

physchains – physical chains computed from computePhysChains

calc_cuml_midpt(vals)

For calculating cumulative quantities around midpoint, so result is independent of binning.

chainFileToStatsFile(chainfname, statfname, burn=0, thin=10, randsamples=None, confint=68.269, mode='hdf5')

Simple wrapper to convert chain file to phys stats file.

Parameters:

• chainfname – input chain file name

• statfname – output statistics file name

• burn – how many iterations to remove off input

• thin – discard every N entries

• randsamples – select N random samples after burn (ignores thin)

• confint – confidence interval to output

• mode – write file mode (only hdf5 supported so far)

computePhysChains(chain)

Compute set of chains for each physical quantity.

Parameters:

chain – list/array of free parameter values (can be loaded using loadChainFromFile)

computePhysPars(pars)

Compute physical quanities for an individual set of parameters. Returned as a dict.

Parameters:

pars – Pars() object with parameters to compute physical parameters for.

loadChainFromFile(chainfname, burn=0, thin=10, randsamples=None)

Get list of parameter values from chain.

Parameters:

• chainfname – input chain HDF5 file

• burn – how many iterations to remove off input

• thin – discard every N entries

• randsamples – randomly sample N samples if set

loadPhysChains(filename, burn=0, thin=1)

Load a Phys chain from the HDF5 file given.

Parameters:

• filename (str) – input hdf5 filename

• burn (int) – discard first N samples

• thin (int) – load very N samples

Returns a dict.

savePhysChains(filename, physchain, thin=1, discard=0)

Save the Phys chains (from computePhysChains) in an HDF5 file.

Parameters:

• filename (str) – output hdf5 filename

• thin (int) – save every N samples from chain

• discard (int) – discard first N samples

writeStatsToFile(fname, stats, mode='hdf5')

Write computed statistics to output filename.

Each physical variable is written with a column, with the 1-sigma +- error bars.

Parameters:

mode – write file mode (only hdf5 supported so far)

mbproj2d.phys.fracMassHalf(snum, edges)

Fraction of mass of shell which is in the inner and outer halves of (r1+r2)/2.

mbproj2d.phys.make2dtuples(v)

Convert (a,b) to ((a,b),), but keep ((a,b),…).

Surface brightness

Convenience routines

mbproj2d.convenience.PSFLoad(filename, hdu)

Load a PSF image from filename given index/name of HDU.

PSF should have CRPIX1/2 as origin and CDELT1 as pixel size (arcsec)

Parameters:

• filename – FITS filename containing PSF

• hdu – HDU name/index for PSF

mbproj2d.convenience.imageLoad(img_fname, exp_fname, rmf='rmf.fits', arf='arf.fits', emin_keV=0.5, emax_keV=2.0, exp_novig_fname=None, origin=None, pix_origin=None, mask_fname=None, psf=None, band_name=None, pad=0)

A helper to construct an Image from input fits images

Parameters:

• img_fname – input fits image filename

• exp_fname – input exposure fits image filename

• rmf – response filename

• arf – arf filename

• emin_keV – minimum energy (keV)

• emax_keV – maximum energy (keV)

• exp_novig_fname – unvignetted exposure filename (optional)

• origin – origin to measure coordinates from (ra_deg, dec_deg) or SkyCoord. By default uses FITS reference coordinate.

• pix_origin – override above with origin (y,x) in pixels (optional)

• mask_fname – name of filename to get mask image (uses exposure>0 by default)

• band_name – name for band (default “X.XX_Y.YY” with emin/emax)

• pad – number of pixels to pad image size with zeros at edge (to avoid convolution wrapping)

Exposure maps loaded in the Band object are given names ‘expmap’ and optionally ‘expmap_novig’

Utilities

Collection of useful functions.

class mbproj2d.utils.AtomicWriteFile(filename)

Write to a file, renaming to final name when finished.

class mbproj2d.utils.CacheOnKey(key)

Cacheing class which writes to a pickled file in a subdirectory.

Items are indexed via md5 hash.

The first two characters are used to create a subdirectory within cachedir to reduce the number of files per directory.

exists()

Return whether cached value exists.

read()

Read cached value, or returns KeyError if not found.

write(val)

Update or write the cached value.

class mbproj2d.utils.ConvPSFHelper(psf)

Helper class for doing convolution of image with PSF.

Parameters:

psf – input PSF array (2D image)

doConv(inimg, outimg)

Do convolution with PSF.

Parameters:

• inimg – where to get input

• outimg – where to place output

updatePSF(psf)

Set convolution code to use new PSF.

class mbproj2d.utils.WithLock(filename)

Hacky lockfile class.

mbproj2d.utils.binImage(img, factor, mean=False)

Bin up image by factor.

Parameters:

• img – input image array

• factor – integer bin factor for both dimensions

• mean – If True, calculate means, otherwise calculate sums

mbproj2d.utils.calcChi2(model, data, error)

Calculate chi2 between model and data.

mbproj2d.utils.calcMedianErrors(results)

Take a set of repeated results, and calculate the median and errors (from perecentiles).

mbproj2d.utils.calcNeSqdToNormPerKpc3(cosmo)

Compute factor to convert a ne-squared to a norm per kpc3.

mbproj2d.utils.cashLogLikelihood(data, model)

Calculate log likelihood of Cash statistic.

mbproj2d.utils.diffCube(a, b)

Calculate a**3-b**3.

mbproj2d.utils.diffQuart(a, b)

Calculate a**4-b**4.

mbproj2d.utils.diffSqr(a, b)

Calculate a**2-b**2.

mbproj2d.utils.loadChainFromFile(chainfname, pars, burn=0, thin=10, randsamples=None)

Get list of parameter values from chain. Walkers are collapsed into one dimension. Output dimensions are (numsamples, numpars).

Parameters:

• chainfname – input chain HDF5 file

• pars – Pars() object to check against chain parameters

• burn – how many iterations to remove off input

• thin – discard every N entries

• randsamples – randomly sample N samples if set (ignores thin)

mbproj2d.utils.projectionVolume(R1, R2, y1, y2)

Return the projected volume of a shell of radius R1->R2 onto an annulus on the sky of y1->y2.

this is the integral: Int(y=y1,y2) Int(x=sqrt(R1^2-y^2),sqrt(R2^2-y^2)) 2*pi*y dx dy = Int(y=y1,y2) 2*pi*y*( sqrt(R2^2-y^2) - sqrt(R1^2-y^2) ) dy

This is half the total volume (front only)

mbproj2d.utils.projectionVolumeMatrix(radii)

Calculate volumes (front and back) using a matrix calculation.

Dot matrix with emissivity array to compute projected surface brightnesses.

Output looks like this: >>> utils.projectionVolumeMatrix(N.arange(5)) array([[ 4.1887902 , 7.55593906, 6.57110358, 6.4200197 ],

[ 0. , 21.76559237, 26.1838121 , 21.27257712], [ 0. , 0. , 46.83209821, 49.71516053], [ 0. , 0. , 0. , 77.57748023]])

mbproj2d.utils.run_irfft2(x)

Call numpy or pyfftw for irfft2.

mbproj2d.utils.run_rfft2(x)

Call numpy or pyfftw for rfft2.

mbproj2d.utils.symmetriseErrors(data)

Take numpy-format data,+,- and convert to data,+-.

mbproj2d.utils.uprint(*args, file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Unbuffered print.

XSpec helper routines

Module to interrogate xspec to get count rates and luminosities given model parameters.

class mbproj2d.xspechelper.XSpecContext

Manager to handle creating xspec instance and cleaning up.

class mbproj2d.xspechelper.XSpecHelper

A helper to get count rates for temperature and densities.

changeResponse(rmf, arf, minenergy_keV, maxenergy_keV)

Create a fake spectrum using the response and use energy range given.

dummyResponse()

Make a wide-energy band dummy response.

getFlux(emin_keV, emax_keV)

Get flux from current model in energy range given.

getRate()

Get rate from current model in current noticed band.

loadXCM(xcmfile)

Load an xcm file into xspec.

readResult()

Return result from xspec.

setAbund(abund)

Set model abundance.

setApec(NH_1022pcm2, T_keV, Z_solar, redshift, norm)

Set an absorbed apec model.

setPlaw(NH_1022pcm2, gamma, norm)

Set an absorbed powerlaw model.

throwAwayOutput()

Ignore output from program until no more data available.

mbproj2d.xspechelper.deleteFile(filename)

Delete file, ignoring errors.

Parallel routines

Module to run a forked parallel process for evalulating functions, one at a time (ForkParallel) or using a queue of input data (ForkQueue).

class mbproj2d.forkparallel.ForkQueuePool(func, ninstances, initfunc=None, argc=None, argv=None)

Execute function in multiple forked processes.

childLoop()

Wait for commands on the socket and execute.

execute(argslist)

Execute the list of items on the queue.

This version cheats by just splitting the input up into equal-sized chunks.

Maybe the chunksize should be smaller, if some chunks would take much less time

finish()

Close child forks and close sockets.

map(func, args)

Allow Pool to be used like multiprocessing.Pool

Note func is ignored here.

mbproj2d.forkparallel.recvItem(sock)

Receive pickled item from socket.

mbproj2d.forkparallel.recvLen(sock, length)

Receive exactly length bytes from socket.

mbproj2d.forkparallel.sendItem(sock, item)

Pickle and send item to socket using size + pickled protocol.