brahmap.lbsim.LBSim_InvNoiseCovLO_Circulant

Bases: BlockDiagInvNoiseCovLO

A block-diagonal linear operator where each diagonal block represents the inverse of a circulant noise covariance matrix \(N^{-1}\).

Parameters:

Name	Type	Description	Default
obs	Observation | List[Observation]	An instance of the Observation class or a list of the same	required
input	dict | ArrayLike	The input array or data defining the operator. It can be a dictionary that maps the detector name to the corresponding input array OR a single array that is used for all the detectors	required
input_type	Literal['covariance', 'power_spectrum']	Specifies whether the input is a covariance array or a power spectrum array, by default "power_spectrum"	'power_spectrum'
dtype	DTypeFloat	The data type of the operator, by default np.float64	float64

Methods:

Name	Description
reset_counters	Resets matrix-vector product counter to zero.
dot	Numpy-like dot() method.
matvec	Matrix-vector multiplication method.
to_array	Returns the dense form of the linear operator as a 2D NumPy array.
get_inverse	Returns the inverse block-diagonal covariance operator.

Attributes:

Name	Type	Description
dtype	DTypeLike	The data type of the operator.
nargin	int	Size of the input vector \(x\), i.e. the number of columns of the operator
nargout	int	Size of the output vector \(A(x)\), i.e. the number of rows of the operator
symmetric	bool	Indicates whether the operator is symmetric or not
shape	tuple[int, int]	A tuple (nargout, nargin) representing the shape of the operator
nMatvec	int	The number of matrix-vector multiplications computed so far
T	LinearOperator	The transpose operator
H	LinearOperator	The adjoint operator
block_list	list[LinearOperator]	A list of linear operators representing the individual diagonal blocks.
num_blocks	int	The total number of diagonal blocks in the operator.
row_size	NDArray[integer]	Array containing the number of rows for each block.
col_size	NDArray[integer]	Array containing the number of columns for each block.
size	int	Array containing the number of rows/columns for each block
diag	NDArray[number]	Array containing the diagonal of the operator

Source code in brahmap/lbsim/lbsim_noise_operators.py

class LBSim_InvNoiseCovLO_Circulant(BlockDiagInvNoiseCovLO):
    """A block-diagonal linear operator where each diagonal block
    represents the inverse of a circulant noise covariance matrix $N^{-1}$.

    Parameters
    ----------
    obs : lbs.Observation | List[lbs.Observation]
        An instance of the `Observation` class or a list of the same
    input : dict | npt.ArrayLike
        The input array or data defining the operator. It can be a
        dictionary that maps the detector name to the corresponding
        input array OR a single array that is used for all the detectors
    input_type : Literal["covariance", "power_spectrum"], optional
        Specifies whether the `input` is a covariance array or a power
        spectrum array, by default `"power_spectrum"`
    dtype : DTypeFloat, optional
        The data type of the operator, by default `np.float64`
    """

    def __init__(
        self,
        obs: lbs.Observation | List[lbs.Observation],
        input: dict | npt.ArrayLike,
        input_type: Literal["covariance", "power_spectrum"] = "power_spectrum",
        dtype: DTypeFloat = np.float64,
    ) -> None:
        if isinstance(obs, lbs.Observation):
            obs_list = [obs]
        else:
            obs_list = obs

        block_size = []
        block_input: Any = None

        if isinstance(input, dict):
            block_input = []

            for obs in obs_list:
                # if input is a dict
                for det_idx in range(obs.n_detectors):
                    block_size.append(obs.n_samples)

                    resized_input = self.__resize_input(
                        new_size=obs.n_samples,
                        input=input[obs.name[det_idx]],
                        input_type=input_type,
                        dtype=dtype,
                    )

                    block_input.append(resized_input)

        elif isinstance(input, (np.ndarray, list)):
            block_input_dict: dict = {}

            for obs in obs_list:
                for det_idx in range(obs.n_detectors):
                    # if input is an array or a list, it will be taken as same for all
                    # the detectors available in the observation
                    block_size.append(obs.n_samples)

                    if obs.n_samples not in block_input_dict:
                        resized_input = self.__resize_input(
                            new_size=obs.n_samples,
                            input=input,
                            input_type=input_type,
                            dtype=dtype,
                        )

                        block_input_dict[obs.n_samples] = resized_input
            block_input = block_input_dict
        else:
            raise ValueError(
                "The input must be an array or a list or a dictionary that maps detector "
                "names to their covariance/power spectrum"
            )

        super(LBSim_InvNoiseCovLO_Circulant, self).__init__(
            InvNoiseCovLO_Circulant,
            block_size=block_size,
            block_input=block_input,
            input_type=input_type,
            dtype=dtype,
        )

    def __resize_input(
        self, new_size, input, input_type, dtype
    ) -> npt.NDArray[np.number]:
        if input_type == "covariance":
            # if the size of the returned array is smaller than new_size, it
            # will be captured by the InvNoiseCovLO_Circulant class
            # automatically

            # Slicing the input array here is probably not the best choice as
            # it breaks the symmetry of the covariance and renders it
            # non-circulant. Same goes for slicing the covariance computed
            # through power spectrum. The best solution would be to create per
            # observation, per detector operators independently and supply them
            # to `BlockDiagInvNoiseCovLO`
            return input[:new_size]
        elif input_type == "power_spectrum":
            input_size = len(input)
            if input_size > new_size:
                new_input = scipy.fft.ifft(
                    input,
                    workers=MPI_UTILS.nthreads_per_process,
                )[:new_size]  # new covariance
                new_input = scipy.fft.fft(  # type: ignore
                    new_input,
                    workers=MPI_UTILS.nthreads_per_process,
                ).real.astype(
                    dtype=dtype,
                    copy=False,
                )  # new ps
                return new_input
            else:
                # If input size is equal to expected size, it will be fine.
                # If it is smaller, InvNoiseCovLO_Circulant class will
                # throw an error automatically
                return input
        return np.asarray(input)

Attributes

dtype: npt.DTypeLike property writable

The data type of the operator.

Returns:

Type	Description
DTypeLike	The NumPy data type of the operator

nargin: int property

Size of the input vector \(x\), i.e. the number of columns of the operator

Returns:

Type	Description
int	The number of input columns

nargout: int property

Size of the output vector \(A(x)\), i.e. the number of rows of the operator

Returns:

Type	Description
int	The number of output rows

symmetric: bool property

Indicates whether the operator is symmetric or not

Returns:

Type	Description
bool	True if symmetric, False otherwise

shape: Tuple[int, int] property

A tuple (nargout, nargin) representing the shape of the operator

Returns:

Type	Description
tuple[int, int]	A tuple (nrows, ncols)

nMatvec: int property

The number of matrix-vector multiplications computed so far

Returns:

Type	Description
int	The number of matrix-vector multiplications performed

T: LinearOperator property

The transpose operator

Returns:

Type	Description
LinearOperator	The transpose of this linear operator

H: LinearOperator property

The adjoint operator

Returns:

Type	Description
LinearOperator	The Hermitian adjoint of this linear operator

block_list: List[LinearOperator] property

A list of linear operators representing the individual diagonal blocks.

Returns:

Type	Description
list[LinearOperator]	A list of linear operators representing the individual diagonal blocks

num_blocks: int property

The total number of diagonal blocks in the operator.

Returns:

Type	Description
int	The total number of diagonal blocks in the operator

row_size: npt.NDArray[np.integer] property

Array containing the number of rows for each block.

Returns:

Type	Description
NDArray[integer]	Array containing the number of rows for each block

col_size: npt.NDArray[np.integer] property

Array containing the number of columns for each block.

Returns:

Type	Description
NDArray[integer]	Array containing the number of columns for each block

size: int property

Array containing the number of rows/columns for each block

Returns:

Type	Description
int	Array containing the number of rows/columns for each block

diag: npt.NDArray[np.number] property

Array containing the diagonal of the operator

Returns:

Type	Description
NDArray[number]	Array containing the diagonal of the operator

Functions

reset_counters() -> None

Resets matrix-vector product counter to zero.

Source code in brahmap/base/linop.py

def reset_counters(self) -> None:
    """Resets matrix-vector product counter to zero."""
    self._nMatvec = 0

dot(x) -> npt.NDArray[np.number]

Numpy-like dot() method.

Parameters:

Name	Type	Description	Default
x	Any	The input vector or object to multiply with.	required

Returns:

Type	Description
NDArray[number]	The result of the dot product.

Source code in brahmap/base/linop.py

def dot(self, x) -> npt.NDArray[np.number]:
    """Numpy-like dot() method.

    Parameters
    ----------
    x : Any
        The input vector or object to multiply with.
    Returns
    -------
    npt.NDArray[np.number]
        The result of the dot product.
    """
    return self.__mul__(x)

matvec(x) -> npt.NDArray[np.number]

Matrix-vector multiplication method.

The matvec method encapsulates the matvec routine specified at construct time, to ensure the consistency of the input and output arrays with the operator's shape.

Parameters:

Name	Type	Description	Default
x	NDArray[number]	The input vector \(x\) to be multiplied by the operator	required

Returns:

Type	Description
NDArray[number]	The result of the matrix-vector multiplication \(A(x)\)

Source code in brahmap/base/linop.py

def matvec(self, x) -> npt.NDArray[np.number]:
    """
    Matrix-vector multiplication method.

    The `matvec` method encapsulates the `matvec`
    routine specified at construct time, to ensure the
    consistency of the input and output arrays with the
    operator's shape.

    Parameters
    ----------
    x : npt.NDArray[np.number]
        The input vector $x$ to be multiplied by the operator

    Returns
    -------
    npt.NDArray[np.number]
        The result of the matrix-vector multiplication $A(x)$
    """
    x = np.asanyarray(x, dtype=self.dtype)
    M, N = self.shape

    # check input data consistency
    N = int(N)
    try:
        x = x.reshape(N)
    except ValueError:
        msg = (
            f"The size of the input array is incompatible with the "
            f"dimensions required by the operator\n"
            f"size of the input array: {x.size}\n"
            f"shape of the operator: {self.shape}"
        )
        msg = f"{self.__class__.__name__}: " + msg
        raise ValueError(msg)

    y = self.__matvec(x)

    # check output data consistency
    M = int(M)
    try:
        y = y.reshape(M)
    except ValueError:
        msg = (
            f"The size of the output array is incompatible with the "
            f"dimensions required by the operator\n"
            f"size of the output array: {y.size}\n"
            f"shape of the operator: {self.shape}"
        )
        msg = f"{self.__class__.__name__}: " + msg
        raise ValueError(msg)

    return y

to_array() -> npt.NDArray[np.number]

Returns the dense form of the linear operator as a 2D NumPy array.

Warning

This method first allocates a NumPy array of shape self.shape and data-type self.dtype, and then fills them with numbers. As such, for a large linear operator, it can occupy an enormous amount of memory and crash your system. Don't use it unless you understand the risk!

Returns:

Type	Description
NDArray[number]	The dense 2D array representation of the linear operator

Source code in brahmap/base/linop.py

def to_array(self) -> npt.NDArray[np.number]:
    """Returns the dense form of the linear operator as a 2D NumPy array.

    !!! Warning

        This method first allocates a NumPy array of shape `self.shape`
        and data-type `self.dtype`, and then fills them with numbers. As
        such, for a large linear operator, it can occupy an enormous
        amount of memory and crash your system. Don't use it unless you
        understand the risk!

    Returns
    -------
    npt.NDArray[np.number]
        The dense 2D array representation of the linear operator
    """
    n, m = self.shape
    H = np.empty((n, m), dtype=self.dtype)
    ej = np.zeros(m, dtype=self.dtype)
    for j in range(m):
        ej[j] = 1.0
        H[:, j] = self * ej
        ej[j] = 0.0
    return H

get_inverse() -> BaseBlockDiagInvNoiseCovLinearOperator

Returns the inverse block-diagonal covariance operator.

Returns:

Type	Description
BaseBlockDiagInvNoiseCovLinearOperator	The inverse block-diagonal covariance operator

Source code in brahmap/base/noise_ops.py

def get_inverse(self) -> "BaseBlockDiagInvNoiseCovLinearOperator":
    """Returns the inverse block-diagonal covariance operator.

    Returns
    -------
    BaseBlockDiagInvNoiseCovLinearOperator
        The inverse block-diagonal covariance operator
    """
    inverse_list = [
        cast(NoiseCovLinearOperator, block).get_inverse()
        for block in self.block_list
    ]
    return BaseBlockDiagInvNoiseCovLinearOperator(
        block_list=inverse_list,
    )