mgefit 5.0.15

MgeFit: Multi-Gaussian Expansion Fitting of Galactic Images

The MgeFit Package

MgeFit: Multi-Gaussian Expansion Fitting of Galactic Images

MgeFit is a Python implementation of the robust and efficient Multi-Gaussian Expansion (MGE) fitting algorithm for galactic images of Cappellari (2002).

The MGE parameterization is useful in the construction of realistic dynamical models of galaxies (see JAM modelling), for PSF deconvolution of images, for the correction and estimation of dust absorption effects, or for galaxy photometry.

Contents

• The MgeFit Package

  • Attribution

  • Installation

  • Usage Examples

• mge_fit_sectors

  • Purpose

  • Calling Sequence

  • Input Parameters

  • Optional Keywords

  • Output Parameters

• mge_fit_1d

  • Purpose

  • Calling Sequence

  • Input Parameters

  • Optional Keywords

  • Output Parameters

  • Example

• find_galaxy

  • Purpose

  • Calling Sequence

  • Input Parameters

  • Optional Keywords

  • Output Parameters

  • Example

• License

• Changelog

Attribution

If you use this software for your research, please cite Cappellari (2002). The BibTeX entry for the paper is:

@Article{Cappellari2002,
    author = {{Cappellari}, M.},
    title = {Efficient multi-Gaussian expansion of galaxies},
    journal = {MNRAS},
    eprint = {arXiv:astro-ph/0201430}
    year = {2002},
    volume = {333},
    pages = {400-410},
    doi = {10.1046/j.1365-8711.2002.05412.x}
}

Installation

install with:

pip install mgefit

Without writing access to the global site-packages directory, use:

pip install --user mgefit

To upgrade MgeFit to the latest version use:

pip install --upgrade mgefit

Usage Examples

To learn how to use the MgeFit package, copy, modify and run the example programs in the mgefit/examples directory. It can be found within the main MgeFit package installation folder inside site-packages. The detailed documentation is contained in the docstring of each file, or for the main procedure on PyPi.

---

mge_fit_sectors

Purpose

Approximates the surface brightness of a galaxy with a Multi-Gaussian Expansion (MGE) model, using the robust and automated fitting method of Cappellari (2002).

The measurements are taken along sectors with a previous call to the procedure sectors_photometry in the MgeFit package . All measurements within this program are in the instrumental units of pixels and counts. This routine fits MGE models with constant position angle and common center.

Calling Sequence

from mgefit.mge_fit_sectors import mge_fit_sectors

m = mge_fit_sectors(radius, angle, counts, eps,
         bulge_disk=False, fignum=1, linear=False, negative=False,
         ngauss=None, normpsf=1., outer_slope=4, plot=False, qbounds=None,
         quiet=False, rbounds=None, scale=1., sigmapsf=0., sol=0)

total_counts, sigma, q_obs = m.sol  # assign the solution to variables
print(m.sol.T)  # Print a table of best-fitting MGE parameters

Example programs are in the mgefit/examples directory. It can be found within the main MgeFit package installation folder inside site-packages.

A PDF manual readme_mge_fit_sectors.pdf is contained in the main installation directory mgefit of the package inside site-packages.

Input Parameters

radius: array_like with shape (n,)

Vector containing the radius of the surface brightness measurements, taken from the galaxy center. This is given in units of PIXELS (!) of the image. When fitting multiple images simultaneously this is in pixels units of the high-resolution image.

angle: array_like with shape (n,)

Vector containing the polar angle of the surface brightness measurements, measured from the galaxy major axis (in degrees).

counts: array_like with shape (n,)

Vector containing the surface brightness measurements in COUNTS/PIXEL (!) at the polar coordinates specified by the vectors radius and angle. When fitting multiple images simultaneously this is in counts/pixels units of the high-resolution image.

eps:

Crude estimate for the galaxy characteristic ellipticity eps = 1 - b/a = 1 - q'

Optional Keywords

bulge_disk:

Set bulge_disk=True to perform a non-parametric bulge/disk decomposition using MGE. When bulge_disk=True the Gaussians are divided into two sets, each with a unique axial ratio. The two sets are meant to describe and model the contribution of bulge and disks in spiral or lenticular galaxies, or nuclear disk in ellipticals.

When this keyword is set one should increase ngauss. One should also either set qbounds=None or specify four elements in qbounds for the even/odd Gaussians.

fignum:

Optional number to assign to the plotting window.

linear:

Set this keyword to start with the fully linear algorithm and then optimize the fit with the nonlinear method [see Section 3.4 of Cappellari (2002) for details]. Nowadays using this keyword instroduces a small speed panelty but produces more robust fits and is always reccomended.

negative:

Set this keyword to allow for negative Gaussians in the fit. Use this only if everything else didn’t work or if there is clear evidence that negative Gaussians are actually needed. Negative Gaussians are needed e.g. when fitting a boxy bulge.

ngauss:

Maximum number of Gaussians allowed in the MGE fit. Typical values are in the range 10-20 when linear=False (default: ngauss=12) and 20**2-40**2 when linear=True (default: ngauss=30**2=900).

normpsf: array_like with shape (p,)

This is optional if only a scalar is given for sigmapsf, otherwise it must contain the normalization of each MGE component of the PSF, whose sigma is given by sigmapsf. The vector needs to have the same number of elements as sigmapsf and it must be normalized as normpsf.sum() = 1 (default: normpsf=1).

outer_slope:

This scalar forces the surface brightness profile of the MGE model to decrease at least as fast as R**(-outer_slope) at the largest measured radius (Default: outer_slope=2).

plot:

Set plot=True to plot the best-fitting MGE profiles.

qbounds:

This can be either a two or a four elements vector.

If qbounds has two elements, it gives the minimum and maximum axial ratio q allowed in the MGE fit.

If qbounds has four elements [[qMin1, qMax1], [qMin2, qMax2]], then the first two elements give the limits on q for the even Gaussians, while the last two elements give the limits on q for the odd Gaussians. In this way qbounds can be used to perform disk/bulges decompositions in a way similar to the bulge_disk keyword.

quiet:

Set quiet = True to suppress printed output.

rbounds:

Two elements vector giving the minimum and maximum sigma allowed in the MGE fit. This is in the same PIXELS units as radius.

scale:

The pixel scale in arcsec/pixels. This is only used for the scale on the plots. It has no influence on the output. (default: 1)

sigmapsf: array_like with shape (p,)

Scalar giving the sigma of the PSF, or vector with the sigma of an MGE model for the circular PSF. This is in the same PIXELS units as radius. When fitting multiple images simultaneously this is the PSF of the high-resolution image [see pg. 406 of Cappellari (2002) for details]. (Default: no convolution)

sol:

If this keyword has at least 6 elements in input, the sigma and q_obs will be used as starting point for the optimization. This is useful for testing but is never needed otherwise. The format is described in the Output Parameters below.

Output Parameters

Stored as attributes of the MgeFit class:

.sol: array_like with shape (3, ngauss)

Array with the best-fitting MGE solution:

.sol[0]: array_like with shape (ngauss)

total_counts of the Gaussian components;

.sol[1]: array_like with shape (ngauss)

sigma is the dispersion of the best-fitting Gaussians in pixels;

.sol[2]: array_like with shape (ngauss)

q_obs is the observed axial ratio of the best-fitting Gaussian components [q’ in the notation of Cappellari (2002)].

The relation below gives the Gaussians peak surface brightness surf:

total_counts, q_obs, sigma = m.sol  # Assign MGE solution to variables
surf = total_counts/(2*np.pi*q_obs*sigma**2)

absdev:

Mean absolute deviation between the fitted MGE and the data expressed as a fraction. Good fits to high S/N images can reach values of absdev < 0.02 = 2%.

---

mge_fit_1d

Purpose

Approximates the 1-dim radial profile of the density or surface-brightness of a galaxy with a Multi-Gaussian Expansion (MGE) model, using the robust and automated fitting method of Cappellari (2002).

Calling Sequence

p = mge_fit_1d(x, y, negative=False, ngauss=12, rbounds=None,
               inner_slope=2, outer_slope=3, quiet=False, fignum=1)

total_counts, sigma = p.sol  # Assign the solution to variables
print(m.sol.T)  # Print a table of best-fitting MGE parameters

See usage example at the end of this file.

Input Parameters

x: array_like with shape (n,)

Vector of the logarithmically sampled (Important!) abscissa for the profile that has to be fitted.

y: array_like with shape (n,)

Ordinate of the profile evaluated at the abscissa x.

Optional Keywords

ngauss:

Number of Gaussian on want to fit. Typical values are 10-20.

negative:

Set negative=True to allow for negative Gaussians in the fit. Use this only if there is clear evidence that negative Gaussians are actually needed e.g. to fit a radially increasing profile.

inner_slope:

This scalar forces the surface brightness profile of the MGE model to decrease as R**(-inner_slope) at the smallest measured radius (Default: inner_slope=2). The default is generally ok, but it can be changed e.g. to allow for steeper inner profiles.

outer_slope:

This scalar forces the surface brightness profile of the MGE model to decrease as R^(**outer_slope) at the largest measured radius (Default: outer_slope=3). The default is generally ok, but it can be changed e.g. to allow for a sharper outer truncation.

rbounds:

Two elements vector giving the minimum and maximum sigma allowed in the MGE fit. This is in the same units of x. With this keyword, inner_slope and outer_slope are ignored.

quiet:

Set quiet=True to suppress plots and printed output.

Output Parameters

Stored as attributes of the mge_fit_1d class:

.sol: array_like with shape (2, ngauss)

Contains an array with the best-fitting solution:

.sol[0]: array_like with shape (ngauss)

TotalCounts (area under the 1-dim curve) of the Gaussians.

.sol[1]: array_like with shape (ngauss)

sigma, is the dispersion of the Gaussians in units of x.

IMPORTANT: Given the TotalCounts definition:

• If the fit was done to the projected surface density profile (Msun/pc^2) the Gaussians peak surface densities (Msun/pc^2) are given by:

  surf = TotalCounts/(np.sqrt(2*np.pi)*Sigma)

• If the fit was done to the intrinsic density profile (Msun/pc^3) the Gaussians peak densities (Msun/pc^3) are still given by:

  dens = TotalCounts/(np.sqrt(2*np.pi)*Sigma)

• In either case, one can convert between densities in (Msun/pc^3) and surface densities in (Msun/pc^2) using the formula:

  dens = surf/(np.sqrt(2*np.pi)*Sigma)

  This is eq. (6) of the file readme_mge_fit_sectors.pdf contained in the main installation directory mgefit of the package inside site-packages.

Example

The sequence of commands below was used to generate the one-dimensional MGE fit of Fig.3 in Cappellari (2002)

import numpy as np
from mgefit.mge_fit_1d import mge_fit_1d

n = 300 # Number of sampled points
R = np.geomspace(0.01, 300, n)  # radii must be logarithmically sampled
rho = (1 + R)**(-4)
p = mge_fit_1d(R, rho, ngauss=16)

---

find_galaxy

Purpose

Find the largest region of connected pixels (after smoothing) lying above a given intensity level of the image. This is useful to automatically identify the location and orientation of a galaxy in an image, assuming it is the largest positive fluctuation. The conventions used by this routine are the same as for the rest of the MgeFit`` package.

This procedure uses the weighted first and second moments of the intensity distribution for the computation of the galaxy center and position angle.

Calling Sequence

from mgefit.find_galaxy import find_galaxy

f = find_galaxy(img, binning=5, fraction=0.1, level=None,
                nblob=1, plot=False, quiet=False)
print('Ellipticity: %.3f' % f.eps)

Input Parameters

Img: array_like with shape (nr, nc)

The galaxy images as a 2-dim array.

Optional Keywords

binning: float, optional

pixel scale for the image smoothing applied before selection.

fraction: float, optional

This corresponds (approximately) to the fraction [0 < fraction < 1] of the image pixels that one wants to consider to estimate galaxy parameters (default fraction = 0.1 = 10%)

level: float, optional

Level above which to select pixels to consider in the estimate of the galaxy parameters. This is an alternative to the use of the fraction keyword.

nblob: int, optional

If nblob=1 (default) the procedure selects the largest feature in the image, if nblob=2 the second largest is selected, and so on for increasing values of nblob. This is useful when the galaxy is not the largest feature in the image.

plot: boolean, optional

set plot=True to display an image in the current graphic window showing the pixels used in the computation of the moments.

quiet: boolean, optional

Set quiet = True to suppress printed output.

Output Parameters

Stored as attributes of the find_galaxy class:

.eps:

The galaxy “average” ellipticity Eps = 1 - b/a = 1 - q'.

.pa:

Standard astronomical PA measured counter-clockwise from the image Y axis (assumed to coincide with North). Note: f.pa = 270 - f.theta.

.theta:

Position angle measured clock-wise from the image X axis.

.majoraxis:

Maximum distance of any of the the selected pixels from (xmed, ymed). For a standar galaxy surface brightness this corresponds to the major axis of the selected isophote.

.xpeak:

First index (row) in Img, of the pixel containing the galaxy center. To be precise this coordinate represents the first index of the brightest pixels within a 40x40 pixels box centered on the galaxy average center.

IMPORTANT: The brightest galaxy pixel is the element Img[xpeak, ypeak]. In the plot produced by find_galaxy, the brightest pixel has coordinates (Ypeak, Xpeak), with the two coordinates swapped!

.ypeak:

Second index (column) in Img, of the brightest galaxy pixel.

.xmed:

X coordinate of luminosity weighted galaxy center.

.ymed:

Y coordinate of luminosity weighted galaxy center.

Example

The command below locates the position and orientation of a galaxy in the image Img`` and produces a plot showing the obtained results and the region of the image used in the computation:

f = find_galaxy(img, plot=True)
print(f'Ellipticity: {f.eps:.3f}')

The command below only uses 2% of the image pixels to estimate the intensity weighted moments and show the results:

f = find_galaxy(img, fraction=0.02, plot=True)
print(f'Coordinates of the peak intensity: {f.xpeak} {f.xpeak}')

---

License

Other/Proprietary License

Copyright (c) 1999-2023 Michele Cappellari

This software is provided as is without any warranty whatsoever. Permission to use, for non-commercial purposes is granted. Permission to modify for personal or internal use is granted, provided this copyright and disclaimer are included in all copies of the software. All other rights are reserved. In particular, redistribution of the code is not allowed.

Changelog

V5.0.15: MC, Oxford, 31 March 2023

• mge_print_contours, mge_print_contours_twist: New keyword minlevel.

• mge_print_contours, mge_print_contours_twist: Analytic integral of the central pixel.

• mge_fit_examples: Updated to use minlevel keyword.

• Included documentation of more functions on PyPi.

V5.0.14: MC, Oxford, 24 June 2021

• Formatted documentation as reStructuredText.

V5.0.13: MC, Oxford, 1 October 2020

• With negative=True use Numpy linalg.lstsq instead of Scipy for a better default in the criterion for rank deficiency. Fixed ignoring negative weights in output with negative=True.

V5.0.12: MC, Oxford, 1 October 2018

• Fixed clock DeprecationWarning in Python 3.7. Use SciPy 1.1 maxiter keyword in nnls.

V5.0.11: MC, Oxford, 12 May 2018

• Dropped legacy Python 2.7 support.

V5.0.10: MC, Oxford, 17 April 2018

• Changed imports for mgefit package.

V5.0.9: MC, Oxford, 21 November 2017

• changed sigmapsf and normpsf keywords to lowercase.

V5.0.8: MC, Oxford, 25 May 2017

• _fitfunc() does not return unused status any more, for consistency with the corresponding change to cap_mpfit.

V5.0.7: MC, Oxford, 14 February 2017

• Make plot() callable after the program terminates. Included fignum keyword and removed obsolete debug keyword. Use line colors from current color cycle.

V5.0.6: MC, Oxford, 24 January 2017

• Improved labelling for Matplotlib 2.0.

V5.0.5: MC, Oxford, 18 June 2015

• Fixed plotting issue when combining profiles from multiple images. Thanks to Arianna Picotti (MPIA) for the bug report with examples. Only plot profiles for the best-fitting MGE.

V5.0.4: MC, Atlantic Ocean, 6 June 2015

• Fully broadcast _fitfunc.

V5.0.3: MC, Atlantic Ocean, 28 March 2015

• Make sure qbounds is a Numpy array. Include absdev in the class attributes. Nicely format printed solution.

V5.0.2: MC, Oxford, 24 September 2014

• Improved plotting.

V5.0.1: MC, Oxford, 25 May 2014

• Support both Python 2.7 and Python 3.

V5.0.0: MC, Aspen Airport, 8 February 2014

• Translated from IDL into Python.

V4.1.3: MC, Oxford, 23 January 2013

• Explained optional usage of SOL in input. Removed stop when MPFIT reports over/underflow.

V4.1.2: MC, Oxford, 24 April 2012

• Small change to the treatment of the innermost unresolved Gaussians.

V4.1.1: MC, Oxford, 12 November 2010

• Added keyword /QUIET.

V4.1.0: MC, Oxford, 22 April 2010

• Allow QBOUNDS to have four elements, to perform bulge/disk decompositions similarly to the /BULGE_DISK option.

V4.0.1: MC, Oxford, 6 June 2009

• Added output keyword ABSDEV. Fixed display not being updated while iterating under Windows.

V4.0.0: MC, Windhoek, 5 October 2008

• Added /BULGE_DISK keyword to perform non-parametric bulge/disks decompositions using MGE. Updated MPFIT to version v1.52 2008/05/04, to fix a bug with the required parinfo.tied mechanism. In the new version of MPFIT, which I again renamed MGE_MPFIT, I implemented my previous important modification to improve convergence with MGE_FIT_SECTORS.

V3.9.5: MC, Oxford, 24 September 2008

• Force Gaussians smaller than the PSF, which have a degenerate axial ratio, to have the same axial ratio as the mean of the first two well determined Gaussians.

V3.9.4: MC, Oxford, 16 May 2008

• Use more robust la_least_squares (IDL 5.6) instead of SVDC with /NEGATIVE keyword.

V3.9.3: MC, Leiden, 18 October 2005

• Changed axes labels in plots.

V3.9.2: MC, Leiden, 11 October 2005

• Print iterations of the longer part at the end, not of the short “Gaussian cleaning” part.

V3.9.1: MC, Leiden, 1 May 2005

• Replaced LOGRANGE keyword in example with the new MAGRANGE.

V3.9.0: MC, Leiden, 23 October 2004

• Allow forcing the outer slope of the surface brightness profile of the MGE model to decrease at least as R**-n at the largest measured radius (cfr. version 3.8). Clean the solution at the end of the nonlinear fit as already done in the /LINEAR implementation. It’s almost always redundant, but quick.

V3.8.1: MC, Vicenza, 23 August 2004

• Make sure this routine uses the Nov/2000 version of Craig Markwardt MPFIT which was renamed MGE_MPFIT to prevent potential conflicts with more recent versions of the same routine.

V3.8.0: MC, Leiden, 8 May 2004

• Force the surface brightness of the MGE model to decrease at least as R**-2 at the largest measured radius.

V3.7.6: MC, Leiden, 20 March 2004

• Use updated calling sequence for BVLS.

V3.7.5: MC, Leiden 23 July 2003

• Corrected small bug introduced in V3.73. Thanks to Arend Sluis.

V3.7.4: MC, Leiden, 9 May 2003

• Use N_ELEMENTS instead of KEYWORD_SET to test non-logical keywords.

V3.7.3: MC, Leiden, 7 March 2003

• Force the input parameters to the given bounds if they fall outside the required range before starting the fit. After feedback from Remco van den Bosch.

V3.7.2: MC, Leiden, 13 October 2002

• Added ERRMSG keyword to MPFIT call.

V3.7.1: MC, Leiden 20 May 2002

• Added compilation options.

V3.7.0: MC, Leiden, 23 February 2002

• Added explicit stepsize (STEP) of numerical derivative in parinfo structure, after suggestion by Craig B. Markwardt.

V3.6.0: MC, Leiden, 23 October 2001

• Modified implementation of /NEGATIVE keyword.

V3.5.0: MC, Leiden, 8 October 2001

• Updated documentation.

V3.4.0: MC, Leiden, 20 September 2001

• Added /FASTNORM keyword

V3.3.0: MC, Leiden, 26 July 2001

• Added MGE PSF convolution, central pixel integration and changed program input parameters to make it independent from SECTORS_PHOTOMETRY

V3.2.0: MC, Leiden, 8 July 2001

• Graphical changes: always show about 7 sectors on the screen, and print plots with shared axes.

V3.1.0: MC, Leiden, 27 April 2001

• More robust definition of err in FITFUNC_MGE_SECTORS.

V3.0.0: MC, Padova, July 2000

• Significant changes.

V2.0.0: MC, Leiden, January 2000

• Major revisions.

V1.0.0: Padova, February 1999

• First implementation by Michele Cappellari.