mangadap.util.covariance module
Defines a class used to store and interface with covariance matrices.
Todo
Allow for calculation of the inverse of the covariance matrix.
Instead of 3D covariance cubes being an array of sparse objects, make the whole thing a sparse array.
The usage examples need to be updated!
Usage examples
You can calculate the covariance matrix for a given wavelength
channel in a mangadap.drpfits.DRPFits
object:
# Access the DRP RSS file
from mangadap.drpfits import DRPFits
drpf = DRPFits(7495, 12703, 'RSS', read=True)
# Calculate a single covariance matrix
C = drpf.covariance_matrix(2281)
# Show the result in an image
C.show()
# Access specific elements
print(C[0,0])
# Convert to a 'dense' array
dense_C = C.toarray()
# Get the triplets of the non-zero elements
i, j, v = C.find()
# Write it to disk (overwrite existing file)
C.write('test_covariance.fits', overwrite=True)
The covariance matrix is stored in “coordinate” format in a fits binary table. Since the covariance matrix is symmetric by definition, only those non-zero elements in the upper triangle (\(C_{ij}\) where \(i\leq j\)) are saved (in memory or on disk). You can read an existing covariance matrix fits file:
from mangadap.util.covariance import Covariance
C = Covariance(ifile='test_covariance.fits')
You can calculate a set of covariance matrices or the full covariance cube:
# Access the DRP RSS file
from mangadap.drpfits import DRPFits
drpf = DRPFits(7495, 12703, 'RSS', read=True)
# Calculate the full cube
CC = drpf.covariance_cube() # BEWARE: This may require a LOT of memory
# Calculate fewer but many covariance matrices
channels = [ 0, 1000, 2000, 3000, 4000 ]
C = drpf.covariance_cube(channels=channels)
# Access individual elements
print(C[0,0,0])
print(C[0,0,1000])
# Write the full cube, or set of channels
CC.write('full_covariance_cube.fits') # BEWARE: This will be a BIG file
C.write('covariance_channel_set.fits')
Although you can access the data in the covariance matrix as
explained above, this is generally inefficient because the
Covariance.__getitem__()
function currently cannot handle
slices. If you need to perform operations with the covariance matrix,
you’re better off working with the Covariance.cov
attribute
directly. To do this, you won’t be able to use the aliasing of the
channel indices for 3D covariance matrices.
The Covariance
class also allows you to toggle between
accessing the matrix as a true covariance matrix, or by splitting the
matrix into its variance and correlation components. For example,:
# Get the covariance matrix for a single wavelength channel
from mangadap.drpfits import DRPFits
drpf = DRPFits(7495, 12703, 'RSS', read=True)
C = drpf.covariance_matrix(2281)
# Show the covariance matrix before and after changing it to a
# correlation matrix
C.show()
C.to_correlation()
C.show()
print(C.is_correlation)
C.revert_correlation()
Covariance matrices that have been converted to correlation matrices
can be written and read in without issue. See
Covariance.write()
and Covariance.read()
. For example:
# Get the covariance matrix for a single wavelength channel
import numpy
from mangadap.util.covariance import Covariance
from mangadap.drpfits import DRPFits
drpf = DRPFits(7495, 12703, 'RSS', read=True)
channels = [ 0, 1000, 2000, 3000, 4000 ]
Cov = drpf.covariance_cube(channels=channels)
Cov.to_correlation()
Cov.show(channel=2000)
Cov.write('correlation_matrix.fits', overwrite=True)
Cov.revert_correlation()
Corr = Covariance(ifile='correlation_matrix.fits')
Corr.revert_correlation()
assert not (numpy.abs(numpy.sum(Cov.toarray(channel=2000)
- Corr.toarray(channel=2000))) > 0.0)
Copyright © 2019, SDSS-IV/MaNGA Pipeline Group
- class mangadap.util.covariance.Covariance(inp, input_indx=None, impose_triu=False, correlation=False, raw_shape=None)[source]
Bases:
object
A general utility for storing, manipulating, and file I/O of sparse covariance matrices.
Todo
Can create empty object. Does this make sense?
Warning
Although possible to abstract further, use of
raw_shape
can only be 2D in the current implementation.- Parameters:
inp (scipy.sparse.csr_matrix, numpy.ndarray) – Covariance matrix to store. Input must be covariance data, not correlation data. Data type can be either a single scipy.sparse.csr_matrix object or a 1-D array of them. Assumes all sparse matrices in the input ndarray have the same size. If None, the covariance object is instantiated empty.
input_indx (numpy.ndarray, optional) – If inp is an array of scipy.sparse.csr_matrix objects, this is an integer array specifying a pseudo-set of indices to use instead of the direct index in the array. I.e., if inp is an array with 5 elements, one can provide a 5-element array for
input_indx
that are used instead as the designated index for the covariance matrices. See:__getitem__()
,_grab_true_index()
.impose_triu (
bool
, optional) – Flag to force the scipy.sparse.csr_matrix object to only be the upper triangle of the covariance matrix. The covariance matrix is symmetric such that C_ij = C_ji, so it’s not necessary to keep both values. This flag will force a call to scipy.sparse.triu when setting the covariance matrix. Otherwise, the input matrix is assumed to only have the upper triangle of numbers.correlation (
bool
, optional) – Convert the input to a correlation matrix. The input must be the covariance matrix.raw_shape (
tuple
, optional) – The covariance data is for a higher dimensional array with this shape. For example, if the covariance data is for a 2D image with shape(nx,ny)
– the shape of the covariance array is be(nx*ny, nx*ny)
– provideraw_shape=(nx,ny)
. This is primarily used for reading and writing; see alsotranspose_raw_shape()
. If None, any higher dimensionality is ignored. When theCovariance
object contains multiple covariance arrays, this raw shape should be applicable to all of them.
- Raises:
TypeError – Raised if the input array is not one-dimensional or the input covariance matrix are not scipy.sparse.csr_matrix objects.
ValueError – Raised if the
input_indx
either does not have the correct size or is not required due to the covariance object being a single matrix.
- cov
The covariance matrix stored in sparse format.
- shape
Shape of the full array.
- Type:
tuple
- raw_shape
The covariance data is for a higher dimensional array with this shape. For example, if the covariance data is for a 2D image, this would be
(nx,ny)
and the shape of the covariance array would be(nx*ny, nx*ny)
.- Type:
tuple
- dim
The number of dimensions in the covariance matrix. This can either be 2 (a single covariance matrix) or 3 (multiple covariance matrices).
- Type:
int
- nnz
The number of non-zero covariance matrix elements.
- Type:
int
- input_indx
If
cov
is an array of scipy.sparse.csr_matrix objects, this is an integer array specifying a pseudo-set of indices to use instead of the direct index in the array. I.e., ifcov
is an array with 5 elements, one can provide a 5-element array forinput_indx
that are used instead as the designated index for the covariance matrices. See:__getitem__()
,_grab_true_index()
.- Type:
- inv
The inverse of the covariance matrix. This is not currently calculated!
- Type:
- var
Array with the variance provided by the diagonal of the/each covariance matrix. This is only populated if necessary, either by being requested (
variance()
) or if needed to convert between covariance and correlation matrices.- Type:
numpy.ndarray
- is_correlation
Flag that the covariance matrix has been saved as a variance vector and a correlation matrix.
- Type:
bool
- _grab_true_index(inp)[source]
In the case of a 3D array, return the true array-element index given the pseudo index.
Todo
Search operation is inefficient. Apparently a better option is a feature request for numpy 2.0
- Parameters:
inp (
int
) – Requested index to convert to the real index in the stored array of covariance matrices.- Returns:
The index in the stored array of the requested covariance matrix.
- Return type:
int
- _impose_upper_triangle()[source]
Force
cov
to only contain non-zero elements in its upper triangle.
- apply_new_variance(var)[source]
Using the same correlation coefficients, return a new
Covariance
object with the provided variance.- Parameters:
var (numpy.ndarray) – Variance vector. Must have a length that matches the shape of this
Covariance
instance.- Returns:
A covariance matrix with the same shape and correlation coefficients and this object, but with different variance.
- Return type:
- Raises:
ValueError – Raised if the length of the variance vector is incorrect.
- bin_to_spaxel_covariance(bin_indx)[source]
Propagate the covariance matrix data for the stacked spectra into the full cube.
This (self) should be the covariance/correlation matrix for the binned spectra.
Todo
Does this still need to be tested?
Generalize the nomenclature.
- Parameters:
bin_indx (numpy.ndarray) – The integer vector with the bin associated with each spectrum in the DRP cube. This is the flattened BINID array.
- Returns:
Covariance/Correlation matrix for the spaxelized binned data.
- Return type:
- coordinate_data(reshape=False)[source]
Construct data arrays with the non-zero covariance components in coordinate format.
This procedure is primarily used when constructing the data arrays to write to a fits file.
Regardless of whether or not the current internal data is a covariance matrix or a correlation matrix, the data is always returned as a correlation matrix.
Matching the class convention, the returned data only includes the upper triangle.
If the covariance matrix is 2-dimensional, four columns are output:
\(i,j\): The row and column indices, respectively, of the covariance matrix.
\(rho_{ij}\): The correlation coefficient between pixels \(i\) and \(j\).
\(V_i\): The variance in each pixel; i.e., the value of \(C_{ii} \forall i\). If reshape is True, this will be output as a two-dimenional array.
If the covariance matrix is 3-dimensional, one additional column is output:
\(k\): The indices of the channels associated with each covariance matrix. This is either just the index of the covariance matrix or the provided pseudo-indices of each channel.
If using the reshape option, the \(i,j\) indices are converted to two columns each providing the indices in the associated reshaped array with coordinates \(c_1,c_2\).
- Parameters:
reshape (
bool
, optional) – Reshape the output in the case when \(i,j\) are actually pixels in a two-dimensional image. The shape of the image is expected to be square, such that the shape of the covariance matrix is \(N_x\times N_y\), where \(N_x = N_y\). Each of theI
andJ
output columns are then split into two columns according to associate coordinates, such that there are 8 output columns.- Returns:
Four to seven numpy.ndarray objects depending on if the data has been reshaped into an image. In detail:
If 2D and not reshaping: The four returned objects contain the indices along the first and second axes, the correlation coefficient, and the variance vector.
If 2D and reshaping: The six returned objects are the unraveled indices in the 2D source array (two indices each for the two axes of the covariance matrix, ordered as the first and second indices for the first covariance index then for the second covariance index), the correlation coefficient, and the 2D variance array.
If 3D and not reshaping: The five returned objects contain the indices along the first and second axes, the index of the relevant correlation matrix (see
input_index
), the correlation coefficient, and the variance vector.If 3D and reshaping: The seven returned objects are the unraveled indices in the 2D source array (two indices each for the two axes of the covariance matrix, ordered as the first and second indices for the first covariance index then for the second covariance index), the index of the relevant correlation matrix (see
input_index
), the correlation coefficient, and the 2D variance array.
- Return type:
tuple
- find(channel=None)[source]
Find the non-zero values in the full covariance matrix (not just the upper triangle).
This is a simple wrapper for
full()
and scipy.sparse.find.- Parameters:
channel (
int
, optional) – The pseudo-index of the covariance matrix to plot. Required if the covariance object is 3D.- Returns:
A tuple of arrays
i
,j
, andc
. The arraysi
andj
contain the index coordinates of the non-zero values, andc
contains the values themselves.- Return type:
tuple
- classmethod from_array(covar, cov_tol=None, rho_tol=None, raw_shape=None)[source]
Define a covariance object using a dense array.
- Parameters:
covar (array-like) – Array with the covariance data. The shape of the array must be square. Input can be any object that can be converted to a dense array using the object method
toarray
or usingnumpy.atleast_2d
.cov_tol (
float
, optional) – Any covariance value less than this is assumed to be equivalent to (and set to) 0.rho_tol (
float
, optional) – Any correlation coefficient less than this is assumed to be equivalent to (and set to) 0.raw_shape (
tuple
, optional) – The covariance data is for a higher dimensional array with this shape. For example, if the covariance data is for a 2D image with shape(nx,ny)
– the shape of the covariance array is be(nx*ny, nx*ny)
– provideraw_shape=(nx,ny)
. This is primarily used for reading and writing; see alsotranspose_raw_shape()
. If None, any higher dimensionality is ignored. When theCovariance
object contains multiple covariance arrays, this raw shape should be applicable to all of them.
- Returns:
The covariance matrix built using the provided array.
- Return type:
- Raises:
ValueError – Raised if
covar
could not be converted to a dense array.
- classmethod from_fits(source, ivar_ext='IVAR', transpose_ivar=False, covar_ext='CORREL', impose_triu=False, correlation=False, quiet=False)[source]
Read covariance data from a fits file.
This read operation matches the data saved to a fits file using
write()
. The class can read covariance data written by other programs as long as they have a commensurate format. See the description of thewrite()
method.If the extension names and column names are correct,
Covariance
can read fits files that were not produced explicitly by this method. This is useful for MaNGA DAP products that include the covariance data as one extension among others. The methodsoutput_hdus()
andcoordinate_data()
are provided to produce the data that can be placed in any fits file.The method determines if the output data were reshaped by checking the number of columns in the binary table.
When the covariance data are for a higher dimensional array, the memory order of the higher dimensional array (particularly whether it’s constructed row- or column-major) is important to ensuring the loaded data is correct. This is why we have chosen the specific format for the binary table used to store the covariance data. Regardless of whether the data was written by a row-major or column-major language, the format should be such that this class can properly read and recover the covariance matrix. However, in some cases, you still may need to transpose the inverse variance data; see
transpose_ivar
.- Parameters:
source (
str
, astropy.io.fits.HDUList) – Initialize the object using an astropy.io.fits.HDUList object or path to a fits file.ivar_ext (
str
, optional) – If reading the data fromsource
, this is the name of the extension with the inverse variance data. Default is'IVAR'
. If None, the variance is taken as unity.transpose_ivar (
bool
, optional) – Flag to transpose the inverse variance data before rescaling the correlation matrix. Should only be necessary in some cases when the covariance data was written by a method other thanwrite()
.covar_ext (
str
, optional) – If reading the data fromsource
, this is the name of the extension with covariance data. Default is'CORREL'
.impose_triu (
bool
, optional) – Flag to force the scipy.sparse.csr_matrix object to only be the upper triangle of the covariance matrix. The covariance matrix is symmetric such that \(C_{ij} = C_{ji}\), so it’s not necessary to keep both values. This flag will force a call to scipy.sparse.triu when setting the covariance matrix. Otherwise, the input matrix is assumed to only have the upper triangle of numbers.correlation (
bool
, optional) – Return the matrix as a correlation matrix. Default (False) is to use the data (always saved in correlation format; seewrite()
) to construct the covariance matrix.quiet (
bool
, optional) – Suppress terminal output.
- classmethod from_matrix_multiplication(T, Sigma)[source]
Construct the covariance matrix that results from a matrix multiplication.
The matrix multiplication should be of the form:
\[{\mathbf T} \times {\mathbf X} = {\mathbf Y}\]where \({\mathbf T}\) is a transfer matrix of size \(N_y\times N_x\), \({\mathbf X}\) is a vector of size \(N_x\), and \({\mathbf Y}\) is the vector of length \({N_y}\) that results from the multiplication.
The covariance matrix is then
\[{\mathbf C} = {\mathbf T} \times {\mathbf \Sigma} \times {\mathbf T}^{\rm T},\]where \({\mathbf \Sigma}\) is the covariance matrix for the elements of \({\mathbf X}\). If Sigma is provided as a vector of length \(N_x\), it is assumed that the elements of \({\mathbf X}\) are independent and the provided vector gives the variance in each element; i.e., the provided data represent the diagonal of \({\mathbf \Sigma}\).
- Parameters:
T (scipy.sparse.csr_matrix, numpy.ndarray) – Transfer matrix. See above.
Sigma (scipy.sparse.csr_matrix, numpy.ndarray) – Covariance matrix. See above.
- classmethod from_samples(samples, cov_tol=None, rho_tol=None)[source]
Define a covariance object using descrete samples.
The covariance is generated using numpy.cov for a set of discretely sampled data for an \(N\)-dimensional parameter space.
- Parameters:
samples (numpy.ndarray) – Array with samples drawn from an \(N\)-dimensional parameter space. The shape of the input array must be \(N_{\rm par}\times N_{\rm samples}\).
cov_tol (
float
, optional) – Any covariance value less than this is assumed to be equivalent to (and set to) 0.rho_tol (
float
, optional) – Any correlation coefficient less than this is assumed to be equivalent to (and set to) 0.
- Returns:
An \(N_{\rm par}\times N_{\rm par}\) covariance matrix built using the provided samples.
- Return type:
- Raises:
ValueError – Raised if input array is not 2D or if the number of samples (length of the second axis) is less than 2.
- classmethod from_variance(variance, correlation=False)[source]
Construct a diagonal covariance matrix using the provided variance.
- Parameters:
variance (numpy.ndarray) – The variance vector.
correlation (
bool
, optional) – Upon instantiation, convert theCovariance
object to a correlation matrix.
- full(channel=None)[source]
Return a scipy.sparse.csr_matrix object with both its upper and lower triangle filled, ensuring that they are symmetric.
This method is essentially equivalent to
toarray()
except that it returns a sparse array and can only be used for one channel at a time.- Parameters:
channel (
int
, optional) – The pseudo-index of the covariance matrix to return.- Returns:
The sparse matrix with both the upper and lower triangles filled (with symmetric information).
- Return type:
- Raises:
ValueError – Raised if the object is 3D and channel is not provided.
- output_hdus(reshape=False, hdr=None)[source]
Construct the output HDUs and header that contain the covariance data.
- Parameters:
reshape (
bool
, optional) – Reshape the output in the case when \(i,j\) are actually pixels in a two-dimensional image. The shape of the image is expected to be square, such that the shape of the covariance matrix is \(N_x\times N_y\), where \(N_x = N_y\). Each of theI
andJ
output columns are then split into two columns according to associate coordinates, such that there are 8 output columns.hdr (astropy.io.fits.Header, optional) – astropy.io.fits.Header instance to which to add covariance keywords. If None, a new astropy.io.fits.Header instance is returned.
- Returns:
Returns three objects:
The header for the primary HDU.
An astropy.io.fits.ImageHDU object with the variance vector.
An astropy.io.fits.BinTableHDU object with the correlation coefficients.
- Return type:
tuple
- revert_correlation()[source]
Revert the object from a correlation matrix back to a full covariance matrix.
This function does nothing if the correlation flag has not been flipped. The variances must have already been calculated!
- show(channel=None, zoom=None, ofile=None, log10=False)[source]
Show a covariance/correlation matrix data.
This converts the (selected) covariance matrix to a filled array and plots the array using pyplot.imshow. If an output file is provided, the image is redirected to the designated output file; otherwise, the image is plotted to the screen.
- Parameters:
channel (
int
, optional) – The pseudo-index of the covariance matrix to plot. Required if the covariance object is 3D.zoom (
float
, optional) – Factor by which to zoom in on the center of the image by removing the other regions of the array. E.g. zoom=2 will show only the central quarter of the covariance matrix.ofile (
str
, optional) – If provided, the array is output to this file instead of being plotted to the screen.
- spaxel_to_bin_covariance(bin_indx)[source]
Opposite of
bin_to_spaxel_covariance()
: Revert the covariance matrix to the covariance between the unique binned spectra.This (self) should be the covariance/correlation matrix for the binned spectra redistributed to the size of the original spaxels map.
Todo
Does this still need to be tested?
Warning
This does NOT propagate the covariance between the spaxels into the covariance in the binned data. That operation is done by, e.g.,
mangadap.proc.spectralstack.SpectralStack._stack_with_covariance()
. This only performs the inverse operation ofbin_to_spaxel_covariance()
.- Parameters:
bin_indx (numpy.ndarray) – The integer vector with the bin associated with each spectrum in the DRP cube. This is the flattened BINID array.
- Returns:
Covariance/Correlation matrix for the stacked spectra.
- Return type:
- static square_shape(size)[source]
Determine the shape of a 2D square array resulting from reshaping a vector.
- Parameters:
size (
int
) – Size of the vector.- Returns:
Length of one axis of the square array.
- Return type:
int
- to_correlation()[source]
Convert the covariance matrix into a correlation matrix by dividing each element by the variances.
If the matrix is a correlation matrix already (see
is_correlation
), no operations are performed. Otherwise, the variance vectors are computed, if necessary, and used to normalize the covariance values.A
Covariance
object can be reverted from a correlation matrix usingrevert_correlation()
.
- toarray(channel=None)[source]
Convert the sparse covariance matrix to a dense array, filled with zeros where appropriate.
- Parameters:
channel (
int
, optional) – The pseudo-index of the covariance matrix to plot. If None and the object is 3D, the full 3D matrix is returned.- Returns:
Dense array with the full covariance matrix.
- Return type:
- transpose_raw_shape()[source]
Reorder the covariance array to account for a transpose in the ravel ordering of the source array.
For a higher dimensional source array (see
raw_shape
), the indices in the covariance matrix relevant to source array can be found using numpy.ravel_multi_index (or numpy.unravel_index for vice versa). However, if the source array is transposed, these operations will be invalid.This method constructs a new Covariance object from the existing data but with the indices rearranged for the transposed source array.
Warning
If
raw_shape
it is not defined, a warning is issued and this method simply returns a copy of the current covariance data.
- variance(copy=True)[source]
Return the variance vector(s) of the covariance matrix.
- Parameters:
copy (
bool
, optional) – Return a copy instead of a reference.
- write(ofile, reshape=False, hdr=None, overwrite=False)[source]
Write the covariance object to a fits file.
Objects written using this function can be reinstantiated using
from_fits()
.The covariance matrix (matrices) are stored in “coordinate” format using fits binary tables; see scipy.sparse.coo_matrix. The matrix is always stored as a correlation matrix, even if the object is currently in the state holding the covariance data.
Independent of the dimensionality of the covariance matrix, the written file has a
PRIMARY
extension with the keywordCOVSHAPE
that specifies the original dimensions of the covariance matrix; seeshape
.The correlation data are written to the
CORREL
extension. The number of columns in this extension depends on the provided keywords; seecoordinate_data()
. The column names are:INDXI
,INDXJ
,INDXK
: indices in the covariance matrix. The last index is provided only if the object is 3D.INDXI
andINDXJ
are separated into two columns if the output is reshaped; these columns areINDXI_C1
,INDXI_C2
,INDXJ_C1
,INDXJ_C2
.RHOIJ
: The non-zero correlation coefficients located the specified coordinates
The inverse of the variance along the diagonal of the covariance matrix is output in an ImageHDU in extension
IVAR
.For 3D matrices, if pseudo-indices have been provided, these are used in the
INDXK
column; however, theCOVSHAPE
header keyword only gives the shape of the unique indices of the covariance channels.- Parameters:
ofile (
str
) – File name for the output.reshape (
bool
, optional) – Reshape the output in the case when \(i,j\) are actually pixels in a two-dimensional image. The shape of the image is expected to be square, such that the shape of the covariance matrix is \(N_x\times N_y\), where \(N_x = N_y\). Each of theI
andJ
output columns are then split into two columns according to associate coordinates.hdr (astropy.io.fits.Header, optional) – A header object to include in the PRIMARY extension. The SHAPE keyword will be added/overwritten.
overwrite (
bool
, optional) – Overwrite any existing file.
- Raises:
FileExistsError – Raised if the output file already exists and overwrite is False.
TypeError – Raised if the input
hdr
does not have the correct type.