PCA

class turbustat.statistics.PCA(cube, n_eigs=None, distance=None)

Bases: turbustat.statistics.base_statistic.BaseStatisticMixIn

Implementation of Principal Component Analysis (Heyer & Brunt, 2002)

Parameters:

cube : numpy.ndarray or astropy.io.fits.PrimaryHDU or astropy.io.fits.ImageHDU or SpectralCube

Data cube.

n_eigs : int

Deprecated. Input using compute_pca or run.

distance : Quantity, optional

Distance to object in physical units. The output spatial widths will be converted to the units given here.

Examples

>>> from turbustat.statistics import PCA
>>> from spectral_cube import SpectralCube
>>> import astropy.units as u
>>> cube = SpectralCube.read("adv.fits") # doctest: +SKIP
>>> pca = PCA(cube, distance=250 * u.pc) # doctest: +SKIP
>>> pca.run(verbose=True) # doctest: +SKIP

Attributes Summary

data
distance
eigvals	All eigenvalues.
eigvecs	All eigenvectors.
gamma	Slope of the size-linewidth relation with correction from Brunt & Heyer 2002
gamma_error_range	One-sigma error bounds on gamma.
header
index	Power-law index.
index_error_range	One-sigma error bounds on the index.
n_eigs
need_header_flag
no_data_flag
total_variance	Total variance of all eigenvalues.
var_proportion	Proportion of variance described by the first n_eigs eigenvalues.

Methods Summary

autocorr_images(self[, n_eigs])	Create the autocorrelation of the eigenimages.
autocorr_spec(self[, n_eigs])	Create the autocorrelation spectra of the eigenvectors.
compute_pca(self[, mean_sub, n_eigs, …])	Create the covariance matrix and its eigenvalues.
eigimages(self[, n_eigs])	Create eigenimages up to the n_eigs.
find_spatial_widths(self[, method, …])	Derive the spatial widths using the autocorrelation of the eigenimages.
find_spectral_widths(self[, method])	Derive the spectral widths using the autocorrelation of the eigenvectors.
fit_plaw(self[, xlow, xhigh, fit_method, …])	Fit the size-linewidth relation.
input_data_header(self, data, header[, …])	Check if the header is given separately from the data type.
intercept(self[, unit])	Intercept from the fits, converted out of the log value.
intercept_error_range(self[, unit])	One-sigma error bounds on the intercept.
load_beam(self[, beam])	Try loading the beam from the header or a given object.
load_results(pickle_file)	Load in a saved pickle file.
model(self, x)	Model with the fit parameters from fit_plaw
noise_ACF(self[, n_eigs])	Create the noise autocorrelation function based off of the eigenvalues beyond PCA.n_eigs.
plot_fit(self[, save_name, show_cov_bar, …])	Plot the covariance matrix, bar plot of eigenvalues, and the fitted size-line width relation.
run(self[, show_progress, verbose, …])	Run the decomposition and fitting in one step.
save_results(self, output_name[, keep_data])	Save the results of the SCF to avoid re-computing.
sonic_length(self[, T_k, mu, use_gamma, unit])	Estimate of the sonic length based on a given temperature.
spatial_width(self[, unit])	Spatial widths for the first n_eigs components.
spatial_width_error(self[, unit])	The 1-sigma error bounds on the spatial widths for the first n_eigs components.
spectral_width(self[, unit])	Spectral widths for the first n_eigs components.
spectral_width_error(self[, unit])	The error bounds on the spectral widths for the first n_eigs components.

Attributes Documentation

data

distance

eigvals

All eigenvalues.

eigvecs

All eigenvectors.

gamma

Slope of the size-linewidth relation with correction from Brunt & Heyer 2002

gamma_error_range

One-sigma error bounds on gamma.

header

index

Power-law index.

index_error_range

One-sigma error bounds on the index.

n_eigs

need_header_flag = True

no_data_flag = False

total_variance

Total variance of all eigenvalues.

var_proportion

Proportion of variance described by the first n_eigs eigenvalues.

Methods Documentation

autocorr_images(self, n_eigs=None)

Create the autocorrelation of the eigenimages.

Parameters:

n_eigs : None or int

The number of autocorrelation images to create. When n_eigs is negative, the last -n_eig autocorrelation images are created. If None is given, the number in n_eigs will be returned.

Returns:

acors : np.ndarray

3D array, where the first dimension if the number of eigenvalues.

autocorr_spec(self, n_eigs=None)

Create the autocorrelation spectra of the eigenvectors.

Parameters:

n_eigs : None or int

The number of autocorrelation vectors to create. When n_eigs is negative, the last -n_eig autocorrelation vectors are created. If None is given, the number in n_eigs will be returned.

Returns:

acors : np.ndarray

2D array, where the first dimension if the number of eigenvalues.

compute_pca(self, mean_sub=False, n_eigs='auto', min_eigval=None, eigen_cut_method='value', show_progress=True)

Create the covariance matrix and its eigenvalues.

If mean_sub is disabled, the first eigenvalue is dominated by the mean of the data, not the variance.

Parameters:

mean_sub : bool, optional

When enabled, subtracts the means of the channels before calculating the covariance. By default, this is disabled to match the Heyer & Brunt method.

n_eigs : {int, ‘auto’}, optional

Number of eigenvalues to compute. The default setting is ‘auto’, which requires min_eigval to be set. Otherwise, the number of eigenvalues used can be set using an int. Setting to -1 will use all of the eigenvalues.

min_eigval : float, optional

The cut-off value to determine the number of important eigenvalues. When eigen_cut_method is proportional, min_eigval is the total proportion of variance described up to the Nth eigenvalue. When eigen_cut_method is value, min_eigval is the minimum variance described by that eigenvalue.

eigen_cut_method : {‘proportion’, ‘value’}, optional

Set whether min_eigval is the proportion of variance determined up to the Nth eigenvalue (proportion) or the minimum value of variance (value).

show_progress : bool, optional

Show a progress bar during the creation of the covariance matrix.

eigimages(self, n_eigs=None)

Create eigenimages up to the n_eigs.

Parameters:

n_eigs : None or int

The number of eigenimages to create. When n_eigs is negative, the last -n_eig eigenimages are created. If None is given, the number in n_eigs will be returned.

Returns:

eigimgs : ndarray

3D array, where the first dimension if the number of eigenvalues.

find_spatial_widths(self, method='contour', brunt_beamcorrect=True, beam_fwhm=None, distance=None, diagnosticplots=False, **fit_kwargs)

Derive the spatial widths using the autocorrelation of the eigenimages.

Parameters:

methods : {‘contour’, ‘fit’, ‘interpolate’, ‘xinterpolate’}, optional

Spatial fitting method to use. The default method is ‘contour’ (fits an ellipse to the 1/e contour about the peak; this is the method used by the Heyer & Brunt works). See WidthEstimate2D for a description of all methods.

brunt_beamcorrect : bool, optional

Apply the beam correction described in Chris Brunt’s thesis. A beam will be searched for in the given header (looking for “BMAJ”). If this fails, the value must be given in beam_fwhm with angular units.

beam_fwhm : None of Quantity, optional

When beam correction is enabled, the FWHM beam size can be given here.

distance : Quantity, optional

Distance to object in physical units. The output spatial widths will be converted to the units given here.

diagnosticplots : bool, optional

Plot the first 9 autocorrelation images with the contour fits. Only implemented for method='contour'.

fit_kwargs : dict, optional

Used when method is ‘contour’. Passed to turbustat.statistics.stats_utils.EllipseModel.estimate_stderrs.

find_spectral_widths(self, method='walk-down')

Derive the spectral widths using the autocorrelation of the eigenvectors.

Parameters:

method : {“walk-down”, “fit”, “interpolate”}, optional

Spectral fitting method to use. The default method is ‘walk-down’ (starting at the peak, continue until reaching 1/e of the peak; this is the method used by the Heyer & Brunt works). See WidthEstimate1D for a description of all methods.

fit_plaw(self, xlow=None, xhigh=None, fit_method='odr', verbose=False, **kwargs)

Fit the size-linewidth relation. This is done through Orthogonal Distance Regression (via scipy), or through MCMC (requires installing emcee).

Parameters:

xlow : Quantity, optional

Lower spatial scale limit to consider in the fit.

xhigh : Quantity, optional

Upper spatial scale value limit to consider in the fit.

fit_method : {‘odr’, ‘bayes’}, optional

Set the type of fitting to perform. Options are ‘odr’ (orthogonal distance regression) or ‘bayes’ (MCMC). Note that ‘bayes’ requires the emcee package to be installed.

verbose : bool, optional

Prints out additional information about the fitting and plots the solution.

**kwargs

Passed to bayes_linear when fit_method is bayes.

input_data_header(self, data, header, need_copy=False)

Check if the header is given separately from the data type.

intercept(self, unit=Unit("pix"))

Intercept from the fits, converted out of the log value.

Parameters:

unit : Unit, optional

The spectral unit to convert the widths to. Can be in pixels, or a spectral unit equivalent to the unit specified in the header.

intercept_error_range(self, unit=Unit("pix"))

One-sigma error bounds on the intercept.

Parameters:

unit : Unit, optional

The spectral unit to convert the widths to. Can be in pixels, or a spectral unit equivalent to the unit specified in the header.

load_beam(self, beam=None)

Try loading the beam from the header or a given object.

Parameters:

beam : Beam, optional

The beam.

static load_results(pickle_file)

Load in a saved pickle file.

Parameters:

pickle_file : str

Name of filename to load in.

Returns:

self : Save statistic class

Statistic instance with saved results.

Examples

Load saved results. >>> stat = Statistic.load_results(“stat_saved.pkl”) # doctest: +SKIP

model(self, x)

Model with the fit parameters from fit_plaw

noise_ACF(self, n_eigs=-10)

Create the noise autocorrelation function based off of the eigenvalues beyond PCA.n_eigs. By default the final 10 eigenvectors whose eigenvalues are above the machine precision limit of the data cube’s dtype are used.

Parameters:

n_eigs : int, optional

The number of eigenvalues to use for estimating the noise ACF. The default is to use the last 10 eigenvectors.

plot_fit(self, save_name=None, show_cov_bar=True, show_sl_fit=True, n_eigs=None, color='r', fit_color='k', symbol='o', cov_cmap='viridis', spatial_unit=Unit("pix"), spectral_unit=Unit("pix"), show_residual=True)

Plot the covariance matrix, bar plot of eigenvalues, and the fitted size-line width relation.

Parameters:

save_name : str, optional

Save name for the figure. Enables saving the plot.

show_cov_bar : bool, optional

Show the covariance matrix and eigenvalue variance bar plot.

show_sl_fit : bool, optional

Show the size-line width relation, if fit.

n_eigs : int, optional

Number of eigenvalues to show in the bar plot. Defaults to the automatically-set value (PCA.n_eigs).

color : {str, RGB tuple}, optional

Color to use in the plots. Defaults to red.

fit_color : {str, RBG tuple}, optional

Colour to show the fit line in. Defaults to color when None is given.

symbol : str, optional

Marker shape to plot the data.

cov_cmap : {str, matplotlib colormap}, optional

Colormap to show the covariance matrix in.

show_residual : bool, optional

Plot the fit residuals.

run(self, show_progress=True, verbose=False, save_name=None, mean_sub=False, decomp_only=False, n_eigs='auto', min_eigval=None, eigen_cut_method='value', spatial_method='contour', spectral_method='walk-down', xlow=None, xhigh=None, fit_method='odr', beam_fwhm=None, brunt_beamcorrect=True, spatial_output_unit=Unit("pix"), spectral_output_unit=Unit("pix"))

Run the decomposition and fitting in one step.

Parameters:

show_progress : bool, optional

Show a progress bar during the creation of the covariance matrix. Enabled by default.

verbose : bool, optional

Enables plotting of the results.

save_name : str,optional

Save the figure when a file name is given.

mean_sub : bool, optional

See compute_pca

decomp_only : bool, optional

Only run the PCA decomposition, not the entire procedure to derive the size-linewidth relation. This should be enabled when using PCA_Distance.

n_eigs : {“auto”, int}, optional

See compute_pca

min_eigval : float, optional

See compute_pca

eigen_cut_method : {‘proportion’, ‘value’}, optional

See compute_pca

spatial_method : str, optional

See fit_spatial_widths.

spectral_method : str, optional

See fit_spectral_widths.

xlow : Quantity, optional

See fit_plaw.

xhigh : Quantity, optional

See fit_plaw.

fit_method : str, optional

See fit_plaw.

beam_fwhm : None of Quantity, optional

See fit_spatial_widths.

brunt_beamcorrect : bool, optional

See fit_spatial_widths.

spatial_output_unit : astropy.units.Unit, optional

Pixel, anglular, or physical unit to convert the spatial sizes to when plotting. Defaults to pixels. Physical unit conversion requires a distance to be given.

spectral_output_unit : astropy.units.Unit, optional

Pixel or spectral unit to convert spectral sizes to when plotting. Defaults to pixels. The spectral unit MUST match the spectral unit defined in the data cube.

save_results(self, output_name, keep_data=False)

Save the results of the SCF to avoid re-computing. The pickled file will not include the data cube by default.

Parameters:

output_name : str

Name of the outputted pickle file.

keep_data : bool, optional

Save the data cube in the pickle file when enabled.

sonic_length(self, T_k=<Quantity 10. K>, mu=1.36, use_gamma=True, unit=Unit("pix"))

Estimate of the sonic length based on a given temperature. Uses the intercept from the fit.

Based on sonic.pro used in the Heyer & Brunt PCA implementation.

Because error from the MCMC fit need not be symmetric, the MCMC samples are needed to provide the correct CIs for the sonic length.

Parameters:

T_k : Quantity, optional

Temperature given in units convertible to Kelvin.

mu : float, optional

Factor to multiply by m_H to account for He and metals.

use_gamma : bool, optional

Toggle whether to use gamma or index. See link given in gamma.

Returns:

lambd : Quantity

Value of the sonic length. If distance was provided, this will be in the units given in the distance. Otherwise, the result will be in pixel units.

lambd_error_range : Quantity

The 1-sigma bounds on the sonic length. The units will match lambd.

unit : Unit, optional

The spatial unit to convert the widths to. Can be in pixels, an angular unit, or (if distance is given) a physical unit.

spatial_width(self, unit=Unit("pix"))

Spatial widths for the first n_eigs components.

Parameters:

unit : Unit, optional

The spatial unit to convert the widths to. Can be in pixels, an angular unit, or (if distance is given) a physical unit.

spatial_width_error(self, unit=Unit("pix"))

The 1-sigma error bounds on the spatial widths for the first n_eigs components.

Parameters:

unit : Unit, optional

The spatial unit to convert the widths to. Can be in pixels, an angular unit, or (if distance is given) a physical unit.

spectral_width(self, unit=Unit("pix"))

Spectral widths for the first n_eigs components.

Parameters:

unit : Unit, optional

The spectral unit to convert the widths to. Can be in pixels, or a spectral unit equivalent to the unit specified in the header.

spectral_width_error(self, unit=Unit("pix"))

The error bounds on the spectral widths for the first n_eigs components.

Parameters:

unit : Unit, optional

The spectral unit to convert the widths to. Can be in pixels, or a spectral unit equivalent to the unit specified in the header.