Reference/API

Fitting

mcfine.fitting.chi_square(observed, model, observed_error=None)

Calculate standard chi-square.

If errors are provided, then this is the sum of \(({\rm obs}-{\rm model})^2/{\rm error}^2\). Else the sum of \(({\rm obs}-{\rm model})^2/{\rm model}^2\).

Parameters:

• observed (np.ndarray) – Observed values.

• model (np.ndarray) – Model values.

• observed_error (float or np.ndarray) – The error in the observed values. Defaults to None.

Returns:

The chi-square value.

Return type:

float

mcfine.fitting.gaussian(x, amp, centre, width)

Evaluate a Gaussian on a 1D grid.

Calculates a Gaussian, using \(f(x) = A \exp[-0.5(x - \mu)^2/\sigma^2]\).

Parameters:

• x (np.ndarray) – Grid to calculate Gaussian on.

• amp (float or np.ndarray) – Height(s) of curve peak(s), \(A\).

• centre (float or np.ndarray) – Peak centre(s), \(\mu\).

• width (float or np.ndarray) – Standard deviation(s), \(\sigma\).

Returns:

Gaussian model array

Return type:

np.ndarray

mcfine.fitting.get_nearest_value(data, value)

Get the nearest value in a dataset

Parameters:

• data – array of values to hunt through

• value – value to find nearest in data to

Returns:

Nearest value

mcfine.fitting.get_nearest_values(dataset, keys, values)

Function to parallelise get_nearest_value

Parameters:

• dataset (dict) – Dataset to hunt through

• keys (list) – List of keys to search with

• values (list) – List of values to find in the dataset

Returns

list of the nearest values

mcfine.fitting.get_radex_multiple_components(theta, vel, v_lines, qn_ul)

Wrapper around RADEX to get out the multiple components

mcfine.fitting.hyperfine_structure_lte(t_ex, tau, v_centre, line_width, strength_lines, v_lines, vel, return_hyperfine_components=False, log_tau=True)

Create hyperfine intensity profile.

Takes line strengths and relative velocity centres, along with excitation temperature and optical depth to produce a hyperfine intensity profile.

Parameters:

• t_ex (float) – Excitation temperature (K).

• tau (float) – Total optical depth of the line.

• v_centre (float) – Central velocity of the strongest component (km/s).

• strength_lines (np.ndarray) – Array of relative line strengths

• v_lines (np.ndarray) – Array of relative velocity shifts for the lines (km/s)

• vel (np.ndarray) – Velocity array (km/s)

• line_width (float) – Width of components (assumed to be the same for each hyperfine component; km/s).

• return_hyperfine_components (bool) – Return the intensity for each hyperfine component. Defaults to False.

• log_tau (bool) – If True, will assume tau is in a logarithmic scale. Else is linear. Defaults to True.

Returns:

Intensity for each individual hyperfine component (if return_hyperfine_components is True), and the total

intensity for all components

mcfine.fitting.hyperfine_structure_radex(t_ex, tau, v_centre, line_width, v_lines, vel, return_hyperfine_components=False)

Create hyperfine intensity profile using RADEX.

Takes line strengths and relative velocity centres calculated with RADEX, and produces a spectrum.

Parameters:

• t_ex (float) – Excitation temperature (K).

• tau (float) – Total optical depth of the line.

• v_centre (float) – Central velocity of the strongest component (km/s).

• line_width (float) – Width of components (assumed to be the same for each hyperfine component; km/s).

• v_lines (float) – Velocity of the various components (km/s).

• vel (np.ndarray) – Velocity array (km/s)

• return_hyperfine_components (bool) – Return the intensity for each hyperfine component. Defaults to False.

Returns:

Intensity for each individual hyperfine component (if return_hyperfine_components is True), and the total

intensity for all components

mcfine.fitting.initial_lmfit(params, intensity, intensity_err, vel, strength_lines, v_lines, props, n_comp=1, fit_type='lte', log_tau=True)

Get initial guess for MC walkers from lmfit

This is the residual calculation for LMFIT. Just to get our initial parameters to instantiate the MC

Parameters:

• params – List of lmfit parameters

• intensity – Measured intensity

• intensity_err – Measured error on intensity

• vel – Velocity axis

• strength_lines – Relative line strength

• v_lines – Relative line velocity

• props – List of properties for the line profiles

• n_comp – Number of components to fit. Defaults to 1

• fit_type – Type of fit to do. Either ‘lte’ or ‘radex’. Defaults to ‘lte’

• log_tau – Whether to fit tau in log-space (default) or linear space

Returns:

Residual (chisq) value from the fit

mcfine.fitting.ln_like(theta, intensity, intensity_err, vel, strength_lines, v_lines, props, n_comp=1, fit_type='lte', log_tau=True)

Calculate the (negative) log-likelihood for the model

This is the main meat of the MCMC fitting, given the parameters we calculate the likelihood of the model and return that

Parameters:

• theta – Parameters for the fit

• intensity – Observed intensity

• intensity_err – Intensity uncertainty

• vel – Observed velocity

• strength_lines – Relative line strength

• v_lines – Relative line velocities

• props – List of properties being fit

• n_comp – Number of components to fit. Defaults to 1

• fit_type – Either ‘lte’ (default) or ‘radex’

• log_tau – Whether to fit tau in log-space (default) or linear space

Returns:

Negative ln-likelihood for the model.

mcfine.fitting.ln_prior(theta, vel, props, bounds, n_comp=1)

Apply flat priors to the emcee fitting

This also ensures any velocity components are strictly increasing, to avoid degeneracies there

Parameters:

• theta – Parameters for the fit

• vel – Observed velocity

• props – List of properties being fit

• bounds – Bounds on parameters

• n_comp – Number of components to fit. Defaults to 1

Returns:

0 if within bounds, -infinity otherwise

mcfine.fitting.ln_prob(theta, intensity, intensity_err, vel, strength_lines, v_lines, props, bounds, n_comp=1, fit_type='lte')

Calculate the ln-probability for emcee

This combines the prior (generally flat) with the ln-likelihood to return to emcee

Parameters:

• theta – Parameters for the fit

• intensity – Observed intensity

• intensity_err – Intensity uncertainty

• vel – Observed velocity

• strength_lines – Relative line strength

• v_lines – Relative line velocities

• props – List of properties being fit

• bounds – Bounds on parameters

• n_comp – Number of components to fit. Defaults to 1

• fit_type – Either ‘lte’ (default) or ‘radex’

Returns:

Combined probability (likelihood+prior) for the model.

Return type:

float

mcfine.fitting.multiple_components(theta, vel, strength_lines, v_lines, props, n_comp, fit_type='lte', log_tau=True)

Sum intensities for multiple lines.

Takes n_comp distinct lines, and calculates the total intensity of their various hyperfine lines.

Parameters:

• theta – [t_ex, tau, vel, vel_width] for each component. Should have a length of 4*`n_comp`

• vel – velocity array

• strength_lines (np.ndarray) – Array for relative line strengths

• props – List of properties for the line profiles

• n_comp (int) – Number of distinct components to calculate intensities for

• fit_type (str) – Should be one of ALLOWED_FIT_TYPES. Defaults to ‘lte’

• log_tau (bool) – Whether to fit tau in log or linear space. Defaults to True (log space)

Returns:

The sum of the intensities for all the distinct components.

Return type:

np.ndarray

mcfine.fitting.radex_grid_interp(theta, qn_ul, labels=None)

Interpolate generated RADEX grid to get useful values out without weird edge effects

Parameters:

• theta (list) – property values

• qn_ul (list) – Names for the transitions

• labels (list) – Labels for proeprties

Returns:

tau and t_ex

mcfine.fitting.residual(observed, model, observed_error=None)

Calculate standard residual.

If errors are provided, then this is the sum of \(({\rm obs}-{\rm model})/{\rm error}\). Else the sum of \(({\rm obs}-{\rm model})/{\rm model}\).

Parameters:

• observed (np.ndarray) – Observed values.

• model (np.ndarray) – Model values.

• observed_error (float or np.ndarray) – The error in the observed values. Defaults to None.

Returns:

The residual value.

Return type:

float

mcfine.fitting.round_nearest(x, a)

Round x to the nearest a

Parameters:

• x – value to round

• a – value to round to

HyperfineFitter

class mcfine.fitting.HyperfineFitter(data=None, vel=None, error=None, mask=None, config_file=None, local_file=None)

create_fit_cube(fit_dict_filename=None, n_comp_filename=None, cube_filename=None)

Create upper/lower errors for spectral fit plots

Parameters:

• fit_dict_filename (str) – Name for the filename of fitted parameter dictionary. Defaults to None, which will pull from config.toml

• n_comp_filename (str) – Name for the filename of component number map. Defaults to None, which will pull from config.toml

• cube_filename (str) – Name for the filename of output cube. Defaults to None, which will pull from config.toml

delta_bic_looper(data, error, save=False, overwrite=True, progress=False)

Increase spectral complexity until we hit diminishing returns

This will iteratively build up the spectral model one component at a time, before looping backwards to remove the weaker components which might just be fits to the noise

Parameters:

• data – Observed data

• error – Observed error

• save – Whether to save out or not, defaults to False

• overwrite – Whether to overwrite existing fits. Defaults to False

• progress – Whether to display progress bars or not. Defaults to False

Returns:

fitted number of components, likelihood, and the emcee sampler

emcee_wrapper(data, error, pos, n_dims, n_comp=1, progress=False)

Light wrapper around emcee

The runs the emcee part, with some custom moves and passes the arguments neatly to functions

Parameters:

• data – Observed data

• error – Observed uncertainty

• pos – Position for the walkers

• n_dims – Number of dimensions for the problem

• n_comp – Number of components to fit. Defaults to None

• progress – Whether or not to display progress bar. Defaults to False

encourage_spatial_coherence(input_dir='fit', output_dir='fit_coherence', fit_dict_filename=None, n_comp_filename=None, likelihood_filename=None, reverse_direction=False)

Loop over fits to encourage spatial coherence

This will loop over RA/Dec to potentially replace fits with neighbouring fits. This helps to encourage spatial coherence and can also help in the case of catastrophic misfits

Parameters:

• input_dir – directory for the input fits. Defaults to ‘fit’

• output_dir – directory for the output of the coherence checking. Defaults to ‘fit_coherence’

• fit_dict_filename – Filename structure for the fit dictionary. Default to None, which will choose some generic name

• n_comp_filename – Filename for the n_comp map. Defaults to None, which will choose some generic name

• likelihood_filename – Filename for the likelihood map. Defaults to None, which will choose some generic name

• reverse_direction – Whether to reverse how we step through x/y in the coherence encouragement. Defaults to False, but likely you should run one case forward, then another backward

generate_radex_grid(output_file=None)

Pre-generate RADEX grid, given various bounds and steps

This will be the grid we interpolate over for speed later. As such, the parameters should be set to cover the parameter range you care about. Parameters are stored in config.toml

Parameters:

output_file (str) – Output filename. Defaults to None, which will pull from config.toml

get_fits_from_samples(samples, vel, n_draws=100, n_comp=1)

Get a number of fit lines from an MCMC run

Parameters:

• samples – emcee output

• vel – Velocity grid to evaluate the fit on

• n_draws (int) – Number of draws to pull from samples

• n_comp (int) – Number of components in the fit

Returns:

array of best fit line

make_parameter_maps(fit_dict_filename=None, n_comp_filename=None, maps_filename=None)

Make maps of fitted parameters

Parameters:

• fit_dict_filename (str) – Name for the filename of fitted parameter dictionary. Defaults to None, which will pull from config.toml

• n_comp_filename (str) – Name for the filename of component number map. Defaults to None, which will pull from config.toml

• maps_filename (str) – Name for the filename of output maps. Defaults to None, which will pull from config.toml

multicomponent_fitter(fit_dict_filename=None, n_comp_filename=None, likelihood_filename=None)

Run the multicomponent fitter

This runs everything, essentially, and is the part that you should call. It’ll do LMFIT to get guesses, emcee to properly sample parameter space and then iteratively add components before removing them.

Parameters:

• fit_dict_filename – Filename to save the fitted emcee walkers to. Defaults to None, which will not save anything

• n_comp_filename – Filename to save out fitted number of components. Defaults to None, which will not save anything

• likelihood_filename – Filename to save out likelihood for the best fit. Defaults to None, which will not save anything

parallel_fit_samples(ij, fit_dict_filename=None, n_comp=None)

Pull fit percentiles from a single pixel

parallel_fitting(ij, fit_dict_filename='fit_dict', save=True, overwrite=False)

Parallel function for MCMC fitting.

Wraps up the MCMC fitting to pass off to multiple cores. Because of overheads, it’s easier to farm out multiple fits to multiple cores, rather than run the MCMC with multiple threads.

Parameters:

• ij (tuple) – tuple containing (i, j) coordinates of the pixel to fit.

• fit_dict_filename (str) – Base filename for MCMC ft pickle. Will append coordinates on afterwards. Defaults to fit_dict.

• save (bool) – Save out files? Defaults to True.

• overwrite (bool) – Overwrite existing files? Defaults to False.

Returns:

Number of fitted components and the best-fit likelihood.

parallel_map_making(ij, n_comp=None, fit_dict_filename='fit_dict', n_samples=500)

Pull parameters for map out of a single pixel

run_mcmc(data, error, n_comp=1, save=True, overwrite=False, progress=False, fit_dict_filename='fit_dict.pkl', p0_fit=None)

Run emcee to get a fit out

Parameters:

• data – Observed data

• error – Observed uncertainty

• n_comp – Number of components to fit. Defaults to 1

• save – Whether to save the sampler out. Defaults to True

• overwrite – Whether to overwrite existing fit. Defaults to False

• progress – Whether to display progress bar. Defaults to False

• fit_dict_filename – Name for the sampler. Defaults to fit_dict.pkl

• p0_fit – Initial guess for the fit. Defaults to None, which will use some basic parameters that are likely to be suboptimal

Plotting

HyperfinePlotter

class mcfine.plotting.HyperfinePlotter(data=None, vel=None, error=None, mask=None, config_file=None, local_file=None)

corner(flat_samples, plot_name, n_comp=1)

Make a corner plot

fit(flat_samples, data, error, n_comp=1, plot_name='fit', show_individual_components=True, show_hyperfine_components=False, n_points=1000, n_draws=100, figsize=(10, 4), x_label='Velocity (km s$^{-1}$)', y_label='Intensity (K)')

Plot a fit spectrum

parallel_corner(ij, plot_name='corner', fit_dict_filename=None, n_comp=None)

Wrapper to parallelise corner plotting

parallel_plot_fit(ij, fit_dict_filename=None, n_comp=None, plot_name='fit', show_individual_components=True, show_hyperfine_components=False, n_points=1000, n_draws=100, figsize=(10, 4), x_label='Velocity (km s$^{-1}$)', y_label='Intensity (K)')

Wrapper to parallelise fit spectrum plotting

parallel_step(ij, plot_name='step', fit_dict_filename='fit_dict', n_comp=None)

Wrapper to parallelise step plotting

plot_corner(fit_dict_filename=None, n_comp_filename=None, plot_name=None, grid=None)

Create corner plot for either a spectrum or a (subset of) a cube

Parameters:

• fit_dict_filename (str) – Name for the file containing the fit parameters. Defaults to None, which will pull from config.toml

• n_comp_filename (str) – Name for the file containing the n_comp map. Only used for cubes. Defaults to None, which will pull from config.toml

• plot_name (str) – Output plot name. Defaults to None, which will pull from config.toml

• grid (np.ndarray) – If fitting a cube, can pass a ndarray of 1 (plot) and 0 (do not plot) to avoid making loads of plots. Defaults to None, which will then use the fitting mask and plot everything

plot_fit(fit_dict_filename=None, n_comp_filename=None, plot_name=None, grid=None)

Create fit spectrum plot for either a single spectrum or a (subset of) a cube

Parameters:

• fit_dict_filename (str) – Name for the file containing the fit parameters. Defaults to None, which will pull from config.toml

• n_comp_filename (str) – Name for the file containing the n_comp map. Only used for cubes. Defaults to None, which will pull from config.toml

• plot_name (str) – Output plot name. Defaults to None, which will pull from config.toml

• grid (np.ndarray) – If fitting a cube, can pass a ndarray of 1 (plot) and 0 (do not plot) to avoid making loads of plots. Defaults to None, which will then use the fitting mask and plot everything

plot_step(fit_dict_filename=None, n_comp_filename=None, plot_name=None, grid=None)

Create step plot for either a spectrum or a (subset of) a cube

Parameters:

• fit_dict_filename (str) – Name for the file containing the fit parameters. Defaults to None, which will pull from config.toml

• n_comp_filename (str) – Name for the file containing the n_comp map. Only used for cubes. Defaults to None, which will pull from config.toml

• plot_name (str) – Output plot name. Defaults to None, which will pull from config.toml

• grid (np.ndarray) – If fitting a cube, can pass a ndarray of 1 (plot) and 0 (do not plot) to avoid making loads of plots. Defaults to None, which will then use the fitting mask and plot everything

step(samples, plot_name='step_plot', n_comp=1)

Make a step plot

Utilities

mcfine.utils.check_overwrite(val_dict, key)

Check if overwrite is in a dict.

Parameters:

• val_dict (dict) – Dictionary of overwrites.

• key (str) – Key to check.

Returns:

Bool of whether to overwrite or not.

mcfine.utils.get_dict_val(val_dict, default_val_dict, table, key, logger=None)

Get value from a dictionary, falling to default if not present.

Parameters:

• val_dict (dict) – Provided dictionary of values

• default_val_dict (dict) – Dictionary of fallback values

• table (str) – Table key to search for

• key (str) – Dictionary key to search for

• logger – Optional logger instance

Returns:

Dictionary value

mcfine.utils.load_fit_dict(file_name)

Load a pickled fit dict

Parameters:

file_name (str) – file name to load

Returns:

Fit dictionary

mcfine.utils.print_config_params()

Print out all config parameters, and default values

mcfine.utils.save_fit_dict(fit_dict, file_name)

Save a fit dictionary

Parameters:

• fit_dict – fit dictionary to save

• file_name (str) – file to save to