mbproj2 API reference

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

Models

Define hydrostatic and non-hydrostatic models.

ModelNullPot: non-hydrostatic model, equivalent to spectral fitting ModelHydro: hydrostatic model, parameterized by density profile

class mbproj2.model.Model(annuli, NH_1022pcm2=None)

Base class for different models.

The parameter to the models are provided by a dict mapping the parameter name to a Param object.

Parameters:

• annuli (Annuli) – annuli on sky
• NH_1022pcm2 (float) – absorbing column density in 10^22 cm^-2

computeMassProf(pars)

Compute mass profile.

Parameters:pars (dict[str, Param]) – input parameters Returns:g and potential arrays (CGS; both can be 0)

computeProfs(pars)

Compute profiles of physical parameters.

Parameters:pars (dict[str, Param]) – parameter values Returns:arrays of ne, T and Z

defPars()

Return type:dict[str,Param] Returns:default dict of parameters (names to Param objects).

prior(pars)

Compute prior, given parameters.

Returns log likelihood.

class mbproj2.model.ModelHydro(annuli, mass_cmpt, ne_cmpt, Z_cmpt, NH_1022pcm2=None)

This is a form of the model assuming hydrostatic equilibrium. The temperature is calculated from the density and pressure computed from the mass model.

The Pout_logergpcm3 parameter is the outer pressure in log10 erg/cm^3.

Parameters:

• annuli (Annuli) – annuli to analyse
• mass_cmpt (CmptMass) – dark matter mass component
• ne_cmpt (Cmpt) – density component
• Z_cmpt (Cmpt) – metallicity component
• NH_1022pcm2 (float) – absorbing column density in 10^22 cm^-2

computeMassProf(pars)

Compute g and potential given parameters.

computeProfs(pars)

Calculate profiles assuming hydrostatic equilibrium.

Returns:ne_pcm3, T_keV, Z_solar

class mbproj2.model.ModelHydroEntropy(annuli, mass_cmpt, K_cmpt, Z_cmpt, NH_1022pcm2=None, self_gravity=True)

This is a form of the model assuming hydrostatic equilibrium, but parameterising using the entropy (K), not density.

As the gas mass can’t be calculated directly, the routine iterates a number of times to get the density, updating the potential each time.

The Pout_logergpcm3 parameter is the outer pressure in log10 erg/cm^3.

Parameters:

• annuli (Annuli) – annuli to analyse
• mass_cmpt (CmptMass) – dark matter mass component
• K_cmpt (Cmpt) – density component
• Z_cmpt (Cmpt) – metallicity component
• NH_1022pcm2 (float) – absorbing column density in 10^22 cm^-2
• self_gravity (float) – include loop to iteratively calculate self-gravity of baryonic mass

computeMassProf(pars)

Compute g and potential given parameters.

computeProfs(pars)

Calculate profiles assuming hydrostatic equilibrium.

iterateComputeProfs(pars)

Iteratively compute output profiles.

This is iterative because if we want to include gas self gravity, this changes g.

class mbproj2.model.ModelNullPot(annuli, ne_cmpt, T_cmpt, Z_cmpt, NH_1022pcm2=None)

This is a form of the model without any gravitational potential parameters.

Density, temperature and abunance are all separately fit.

Parameters:

• annuli (Annuli) – annuli to analyse
• ne_cmpt (Cmpt) – density component
• T_cmpt (Cmpt) – temperature component
• Z_cmpt (Cmpt) – metallicity component
• NH_1022pcm2 (float) – absorbing column density in 10^22 cm^-2

mbproj2.model.computeGasAccn(annuli, ne_prof)

Compute acceleration due to gas mass for density profile given.

Model parameters

The parameters to the model.

class mbproj2.param.Param(val, minval=-1e+99, maxval=1e+99, frozen=False)

Model parameter with minimum and maximum value.

Parameters:

• val (float) – value of parameter
• minval (float) – minimum allowed value
• maxval (float) – maximum allowed value
• frozen (bool) – whether parameter is allowed to vary

class mbproj2.param.ParamBase(val, frozen=False)

Base class for parameters.

Parameters:

• val (float) – value of parameter
• frozen (bool) – whether parameter is allowed to vary

prior()

Log prior on parameter.

class mbproj2.param.ParamGaussian(val, prior_mu, prior_sigma, frozen=False)

Parameter with a Gaussian/normal prior.

Parameters:

• val (float) – value of parameter
• prior_mu (float) – centre of prior
• prior_sigma (float) – width of prior
• frozen (bool) – whether parameter is allowed to vary

Components

Components which make up a profile. Each component has a set of parameters (of type ParamXXX).

class mbproj2.cmpt.Cmpt(name, annuli)

Parametrize a profile.

Parameters:

• name – prepended to each model parameter name
• annuli (Annuli) – annuli to analyse with component

computeProf(pars)

Compute profile given model parameters.

Parameters:pars (dict[str, Param]) – parameter values Returns:output profile

defPars()

Return type:dict[str,Param] Returns:default dict of parameters (names to Param objects).

prior(pars)

Given parameters, compute prior.

Parameters:pars (dict[str, Param]) – parameter values Returns:log prior

class mbproj2.cmpt.CmptBeta(name, annuli)

Beta model.

Model parameters are XX_n0 (log base 10), XX_beta and XX_rc (log10 kpc)

Parameters:

• name – prepended to each model parameter name
• annuli (Annuli) – annuli to analyse with component

class mbproj2.cmpt.CmptBinWidthIncr(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, nradbins=5, log=False)

Component where bins have increasing widths. Not sure this is very useful.

Adds _dw_* parameters which are delta-widths

class mbproj2.cmpt.CmptBinned(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, binning=1, interpolate=False, log=False)

A profile made of bins with a parameter for every N bin.

name: name (used as start of parameter names) annuli: Annuli object defval: default value minval: minimum value maxval: maximum value binning: factor to bin annuli (how many bins per parameter) interpolate: interpolate values in intermediate bins log: use 10**values to convert to physical quantity

class mbproj2.cmpt.CmptBinnedJumpPrior(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, binning=1, interpolate=False, log=False, priorjump=0.0)

A binned component using a prior that the values shouldn’t jump by more than the factor given.

Parameters:

• name – name (used as start of parameter names)
• annuli – Annuli object
• defval – default value
• minval – minimum value
• maxval – maximum value
• binning – factor to bin annuli (how many bins per parameter)
• interpolate – interpolate values in intermediate bins
• log – use 10**values to convert to physical quantity
• priorjump – fractional difference allowed to jump between bins, implemented as a prior

class mbproj2.cmpt.CmptBinnedMoveRad(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, nradbins=5, log=False)

Binned data with movable radii.

Parameters:

• name – name (used as start of parameter names)
• annuli (Annuli) – Annuli object
• defval – default value
• minval – minimum value
• maxval – maximum value
• nradbins – number of control points (“bins”) to use
• log – use 10**values to convert to physical quantity

Model parameters: ‘XX_YYY’ and ‘XX_r_YYY’ where YYY goes from 000...999, based on the number of radial bins.

class mbproj2.cmpt.CmptDoubleBeta(name, annuli)

Double beta model.

Model parameters are XX_n0_N (log base 10), XX_beta_N and XX_rc_N (log10 kpc), where N is 1 and 2

Parameters:

• name – prepended to each model parameter name
• annuli (Annuli) – annuli to analyse with component

class mbproj2.cmpt.CmptFlat(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, log=False)

A flat profile.

Parameters:

• name – name (used as parameter name)
• annuli – Annuli object
• defval – default value
• minval – minimum value
• maxval – maximum value
• log – use 10**value to convert to physical quantity

class mbproj2.cmpt.CmptIncr(name, annuli, defval=0.0, minval=-5, maxval=5)

An increasing-inward log component.

class mbproj2.cmpt.CmptIncrMoveRad(name, annuli, defval=0.0, defouter=0.0, minval=-5.0, maxval=5.0, nradbins=5, log=False)

A profile with control points, using interpolation to find the values in between.

The radii of the control points are model parameters XX_r_YYY (in log kpc), where YYY is the index of the annulus 000...999. The y values (XX_YYY) are gradients in log (cm^-3 log kpc^-1).

This model forces the density profile to increase inwards.

class mbproj2.cmpt.CmptInterpolMoveRad(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, nradbins=5, log=False, intbeyond=False)

A profile with control points, using interpolation to find the values in between.

The radii of the control points are parameters (XX_r_999 in log kpc)

Parameters:

• name – used as start of parameter names
• annuli (Annuli) – Annuli object
• defval – default value
• minval – minimum value
• maxval – maximum value
• nradbins – number of control points (“bins”) to use
• interpolate – interpolate values in intermediate bins
• log – use 10**values to convert to physical quantity
• intbeyond – powerlaw interpolate inside and outside radii (assumes constant values if False)

class mbproj2.cmpt.CmptMoveRadBase(name, annuli, defval=0.0, minval=-1e+99, maxval=1e+99, nradbins=5, log=False)

Base class for components with bins which can move.

Parameters:

• name – name (used as start of parameter names)
• annuli (Annuli) – Annuli object
• defval – default value
• minval – minimum value
• maxval – maximum value
• nradbins – number of control points (“bins”) to use
• log – use 10**values to convert to physical quantity

Model parameters: ‘XX_YYY’ and ‘XX_r_YYY’ where YYY goes from 000...999, based on the number of radial bins.

class mbproj2.cmpt.CmptVikhDensity(name, annuli, mode='double')

Density model from Vikhlinin+06, Eqn 3.

Use double mode for 2nd component or single otherwise.

Densities and radii are are log base 10

mbproj2.cmpt.betaprof(rin_cm, rout_cm, n0, beta, rc)

Return beta function density profile

Calculates average density in each shell.

Mass components

CmptMass objects define dark matter potentials.

class mbproj2.cmpt_mass.CmptMass(name, annuli, suffix='')

Base mass component.

Parameters:

• name – prefix for parameters
• annuli (Annuli) – annuli to examine
• suffix – suffix to add to name if set

computeProf(pars)

Compute g_cmps2 and potential_ergpg profiles.

class mbproj2.cmpt_mass.CmptMassArb(annuli, nradbins, suffix=None)

Parametrise mass density using interpolation.

Model parameters are arb_rho_YYY (mass density in log10 g cm^-3) and arb_r_YYY (radii in log10 kpc) for annuli YYY.

Parameters:

• annuli (Annuli) – annuli used
• suffix – suffix to append to name arb in parameters

class mbproj2.cmpt_mass.CmptMassGNFW(annuli, suffix=None)

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

Parameters:

• annuli (Annuli) – Annuli object
• suffix – suffix to append to name gnfw in parameters

class mbproj2.cmpt_mass.CmptMassKing(annuli, suffix=None)

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

Parameters:

• annuli (Annuli) – Annuli object
• suffix – suffix to append to name king in parameters

class mbproj2.cmpt_mass.CmptMassMulti(name, annuli, cmpts, suffix=None)

Multi-component mass profile made up CmptMass objects.

Parameters:

• name – name of component
• annuli (Annuli) – annuli to use
• cmpts (list[CmptMass]) – components to sum

class mbproj2.cmpt_mass.CmptMassNFW(annuli, suffix=None)

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

Model parameters are nfw_logconc (log10 concentration) and nfw_r200_logMpc (log10 r200 in Mpc)

Parameters:

• annuli (Annuli) – Annuli object
• suffix – suffix to append to name nfw in parameters

class mbproj2.cmpt_mass.CmptMassPoint(annuli, suffix=None)

Point mass.

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

annuli: Annuli object suffix: suffix to append to name pt in parameters

Data

Contains classes to represent data to be fit.

class mbproj2.data.Annuli(edges_arcmin, cosmology)

Geometric information about the annuli on the sky.

Parameters:

• edges_arcmin – edges of annuli in arcmin (if N annuli, should be N+1 edges)
• cosmology (Cosmology) – Cosmology object

update(edges_arcmin, cosmology)

Change the annuli. Useful for recalculating models with new grid.

Parameters:

• edges_arcmin – edges of annuli in arcmin
• cosmology (Cosmology) – Cosmology object

class mbproj2.data.Band(emin_keV, emax_keV, cts, rmf, arf, exposures, backrates=None, areascales=None, psfmatrix=None)

Count profile in a band.

Parameters:

• emin_keV – minimum energy of band in keV
• emax_keV – maximum energy of band in keV
• cts – numpy array of counts in each annulus
• rmf – response matrix filename
• arf – ancillary response matrix filename
• exposures – numpy array of exposures in each annulus

optionally: :param backrates: numpy array of rates of cts/s/arcmin^2 in each annulus :param areascales: numpy array of scaling factors to convert from geometric area in annulus to real area (including pixels) :param psfmatrix: matrix to convolve to account for PSF, usually calculated using functions in psfconvolve submodule

calcProjProfile(annuli, ne_prof, T_prof, Z_prof, NH_1022pcm2, backscale=1.0)

Predict profile given cluster profiles.

class mbproj2.data.Data(bands, annuli)

Dataset class.

bands: list of Band objects annuli: Annuli object

mbproj2.data.expandlist(x, length)

If x is a list, check it has the length length. Otherwise, expand item to be a list with length given.

mbproj2.data.loadAnnuli(filename, cosmology, centrecol=0, hwcol=1)

Helper to load annuli from data file.

Data file is in text with whitespace-separated columns.

Parameters:

• cosmology (Cosmology) – cosmology to use
• centrecol (int) – index of column giving centre (arcmin) of annulus
• hwcol (int) – index of column giving half-width (arcmin) of annulus

mbproj2.data.loadBand(filename, emin_keV, emax_keV, rmf, arf, radiuscol=0, hwcol=1, ctcol=2, areacol=3, expcol=4)

Load a band using standard data format.

Cosmology

Routines for calculating distances from cosmology.

TODO: Replace with astopy’s versions?

class mbproj2.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

D_A

Get angular diameter distance in Mpc.

D_L

Get luminosity distance in Mpc.

kpc_per_arcsec

Get number of kpc per arcsec.

Fitting

A fit to the model.

class mbproj2.fit.Fit(pars, model, data)

Class to help fitting model, by keeping track of thawed parameters.

Parameters:

• pars (dict[str,ParamBase]) – parameters for model
• model (Model) – Model to fit
• data (Data) – Data to fit

The parameters pars are for the model, but a parameter called backscale can be included, which controls the scaling of the background

calcProfiles()

Predict model profiles for each band.

doFitting(silent=False, maxiter=10)

Optimize parameters to increase likelihood. Uses scipy’s Nelder-Mead and Powell optimizers, repeating if a new minimum is found.

Parameters:

• silent (bool) – print output during fitting
• maxiter (int) – maximum number of iterations

Returns:

log likelihood

getLikelihood(vals=None)

Get likelihood for parameters given.

Also include are the priors from the various components

likeFromProfs(predprofs)

Given predicted profiles, calculate log likelihood (excluding prior).

Parameters:predprofs (list[numpy.array]) – input profiles

plotProfiles()

Plot surface brightness profiles of model and data if Veusz is installed.

refreshThawed()

Call this after making changes to which parameters are thawed.

save(filename)

Save fit in Python pickle format.

thawedParVals()

Return list of numeric values of thawed parameters.

updateThawed(vals)

Update values of parameter ParamBase objects which are thawed. :param list[float] vals: numerical values of parameters

mbproj2.fit.genericPopulationMinimizer(function, gennewparams, popnum=1000, keepfrac=0.8, maxiter=1000, sigmabreak=0.001)

Minimize using a set of populations.

Markov Chain Monte Carlo

For doing MCMC given Fit object.

class mbproj2.mcmc.MCMC(fit, walkers=100, processes=1, initspread=0.01)

For running Markov Chain Monte Carlo.

Parameters:

• fit (Fit) – Fit object to use for mcmc
• walkers (int) – number of emcee walkers to use
• processes (int) – number of simultaneous processes to compute likelihoods
• initspread (float) – random Gaussian width added to create initial parameters

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)

Run chain.

Parameters:length (int) – length of chain

save(outfilename, thin=1)

Save chain to HDF5 file.

Parameters:

• outfilename (str) – output hdf5 filename
• thin (int) – save every N samples from chain

Physical quantities

Compute physical quantities from Model and parameters, or the MCMC chains.

mbproj2.phys.computePhysChains(chainfilename, model, pars, burn=0, thin=10, randsample=False)

Compute set of chains for each physical quantity.

Parameters:

• chainfilename – input chain filename
• model (Model) – model to use
• pars (dict[str, ParamBase]) – parameters to apply
• burn – skip initial N items in chain
• thin – skip every N iterations in chain
• randsample – randomly sample from chain at thin interval

Returns:

tuple with dict of name with profiles, radial bins in arcmin, with widths, radial bins in kpc, with widths

mbproj2.phys.fracMassHalf(snum, annuli)

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

mbproj2.phys.physFromProfs(model, pars)

Given model and parameters, calculate physical quantities.

Parameters:

• model (Model) – model to use
• pars (dict[str, ParamBase]) – parameters to apply

mbproj2.phys.replayChainPhys(chainfilename, model, pars, burn=0, thin=10, confint=68.269, randsample=False)

Replay chain, computing median physical quantity profiles.

Parameters:

• chainfilename – input physical chain filename
• model (Model) – input model
• pars (dict[str, ParamBase]) – parameters used in model
• confint – total confidence interval (percentage)
• burn – skip initial N items in chain
• thin – skip every N iterations in chain
• randsample – randomly sample chain when thinning

Returns:

medians and confidence interval percentiles

mbproj2.phys.savePhysChain(infilename, outfilename, model, pars, burn=0, thin=10, randsample=False)

Convert parameter chain to physical chain, written to HDF5.

Parameters:

• infilename – input HDF5 chain filename
• outfilename – output HDF5 chain filename
• model (Model) – input model
• pars (dict[str, ParamBase]) – parameters used in model
• burn – throw away initial N parameters
• thin – throw away every N parameters
• randsample – sample values randomly from the chain when thinning

mbproj2.phys.savePhysProfilesHDF5(outfilename, profiles)

Given median profiles from replayChainPhys, save output profiles to hdf5.

mbproj2.phys.savePhysProfilesText(outfilename, profiles)

Given median profiles from replayChainPhys, save output profiles to text.

Helper routines

Various useful functions for the user.

mbproj2.helpers.autoRadialBins(annuli, data, minsn, minbins=2, maxbins=100)

Take radial count profiles and choose bins using number of projected counts.

mbproj2.helpers.estimateDensityProfile(inmodel, data, modelpars)

Do a simplified fit to get the initial density profile.

This assumes an isothermal cluster with constant 0.3 solar metallicity and does not use hydrostatic equilibrium.

mbproj2.helpers.fitBeta(annuli, data, NH_1022pcm2, Z_solar, T_keV, silent=True)

Fit beta density model with isothermal cluster.

Return ne_cmpt, T_cmpt, Z_cmpt, pars

mbproj2.helpers.initialNeCmptBinnedFromBeta(annuli, data, NH_1022pcm2=0.01, Z_solar=0.3, T_keV=3.0)

Return ne component and initial parameters.

mbproj2.helpers.initialNeCmptInterpolMoveRadFromBeta(annuli, data, mode, NH_1022pcm2=0.01, Z_solar=0.3, T_keV=3.0, nradbins=10, minsn=30, minbins=2, maxbins=100, silent=True)

Create a density profile with the possibility to move interpolated bin radii based on an isothermal beta model initial fit.

returns: density component, default density parameters

mode should be: ‘lognbins’, ‘minsn’