brahmap.core.PointingLO

Bases: LinearOperator

A linear operator representing the pointing matrix (projection operator) \(P\).

This class encapsulates the highly sparse projection/de-projection operations (scatter/gather) required to map time samples onto sky pixels or vice versa, according to the given pointing information and map-making solver configuration. It excludes bad pointing samples and pathological pixels while performing projection/de-projection operations. The shape of the operator is [nsamples, new_npix*ncomponents] where ncomponents depends on the number of components being projected on the sky pixels. For instance, for IQU map-making, ncomponents = 3.

Parameters:

Name	Type	Description	Default
processed_samples	ProcessTimeSamples | SharedMemProcessTimeSamples	The pre-processed time samples object containing pointing and map-making metadata	required
solver_type	SolverType | None	The map-making solver configuration to use. If None, it falls back to the solver_type of processed_samples, by default None	None
return_copy	bool	If True, the transposed operator (rmatvec) returns a copy of the shared memory buffer. If False, it returns the shared memory buffer directly. This argument is ignored if processed_samples is not a SharedMemProcessTimeSamples object, by default True	True

Attributes:

Name	Type	Description
solver_type	SolverType	The current map-making solver configuration
return_copy	bool	Whether the transposed operator returns a copy of the shared memory buffer

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.

Source code in brahmap/core/linearoperators.py

class PointingLO(LinearOperator):
    """A linear operator representing the pointing matrix (projection operator) $P$.

    This class encapsulates the highly sparse projection/de-projection
    operations (scatter/gather) required to map time samples onto sky
    pixels or vice versa, according to the given pointing information and
    map-making solver configuration. It excludes bad pointing samples and
    pathological pixels while performing projection/de-projection
    operations. The shape of the operator is `[nsamples, new_npix*ncomponents]`
    where `ncomponents` depends on the number of components being
    projected on the sky pixels. For instance, for `IQU` map-making,
    `ncomponents = 3`.

    Parameters
    ----------
    processed_samples : ProcessTimeSamples | SharedMemProcessTimeSamples
        The pre-processed time samples object containing pointing and
        map-making metadata
    solver_type : SolverType | None, optional
        The map-making solver configuration to use. If `None`, it falls
        back to the `solver_type` of `processed_samples`, by default `None`
    return_copy : bool, optional
        If `True`, the transposed operator (`rmatvec`) returns a copy of the
        shared memory buffer. If `False`, it returns the shared memory buffer
        directly. This argument is ignored if `processed_samples` is not a
        [`SharedMemProcessTimeSamples`][brahmap.core.SharedMemProcessTimeSamples] object, by default `True`

    Attributes
    ----------
    solver_type : SolverType
        The current map-making solver configuration
    return_copy : bool
        Whether the transposed operator returns a copy of the shared memory
        buffer
    """

    def __init__(
        self,
        processed_samples: ProcessTimeSamples | SharedMemProcessTimeSamples,
        solver_type: None | SolverType = None,
        return_copy: bool = True,
    ) -> None:
        ### Some of the functionalities of this class are implemented with C++
        ### extensions. A corresponding full Python implementation is provided in
        ### `tests/py_PointingLO.py` for reference.

        if solver_type is None:
            self.__solver_type = processed_samples.solver_type
        else:
            if int(processed_samples.solver_type) < int(solver_type):
                raise ValueError(
                    "`solver_type` must be lower than or equal to the "
                    "`solver_type` of `processed_samples` object"
                )
            self.__solver_type = solver_type

        self.__return_copy = return_copy

        self.new_npix = processed_samples.new_npix
        self.ncols = processed_samples.new_npix * self.solver_type
        self.nrows = processed_samples.nsamples

        self.pointings = processed_samples.pointings
        self.pointings_flag = processed_samples.pointings_flag

        if self.solver_type > 1:
            self.sin2phi = processed_samples.sin2phi
            self.cos2phi = processed_samples.cos2phi

        self._is_shmem = hasattr(processed_samples, "shared_mem_manager")

        if self._is_shmem:
            assert isinstance(processed_samples, SharedMemProcessTimeSamples)
            self.__shared_mem_mgr = processed_samples.shared_mem_manager
            mgr = self.__shared_mem_mgr

            # Allocated node-level shared memory arrays for transposed product
            self._node_prod, self._win_node_prod = mgr.alloc_shared_node(
                self.ncols,
                processed_samples.dtype_float,
            )

            self._grp_prod, self._win_grp_prod = mgr.alloc_shared_comm(
                self.ncols,
                processed_samples.dtype_float,
                comm=mgr.tree_grp_comm,
                comm_root=0,
            )

            rmatvec_I = self._rmult_I_shmem
            rmatvec_QU = self._rmult_QU_shmem
            rmatvec_IQU = self._rmult_IQU_shmem
        else:
            rmatvec_I = self._rmult_I
            rmatvec_QU = self._rmult_QU
            rmatvec_IQU = self._rmult_IQU

        if self.solver_type == 1:
            super().__init__(
                nargin=self.ncols,
                nargout=self.nrows,
                symmetric=False,
                matvec=self._mult_I,
                rmatvec=rmatvec_I,
                dtype=processed_samples.dtype_float,
            )
        elif self.solver_type == 2:
            super().__init__(
                nargin=self.ncols,
                nargout=self.nrows,
                symmetric=False,
                matvec=self._mult_QU,
                rmatvec=rmatvec_QU,
                dtype=processed_samples.dtype_float,
            )
        else:
            super().__init__(
                nargin=self.ncols,
                nargout=self.nrows,
                symmetric=False,
                matvec=self._mult_IQU,
                rmatvec=rmatvec_IQU,
                dtype=processed_samples.dtype_float,
            )

    def _mult_I(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the matrix-vector product $Pv$ for temperature-only ($I$)
        map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector $v$ of size `new_npix`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector of size `nsamples`
        """

        prod = np.zeros(self.nrows, dtype=self.dtype)

        PointingLO_tools.PLO_mult_I(
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            vec=vec,
            prod=prod,
        )

        return prod

    def _rmult_I(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        temperature-only ($I$) map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector of size `new_npix`
        """

        prod = np.zeros(self.ncols, dtype=self.dtype)

        PointingLO_tools.PLO_rmult_I(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            vec=vec,
            prod=prod,
            comm=MPI_UTILS.comm,
        )

        return prod

    def _rmult_I_shmem(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        temperature-only ($I$) map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector of size `new_npix`
        """

        if self.__shared_mem_mgr.tree_grp_rank == 0:
            self._grp_prod[:] = 0
        if self.__shared_mem_mgr.node_rank == 0:
            self._node_prod[:] = 0

        self._win_grp_prod.Fence(0)
        self._win_node_prod.Fence(0)

        PointingLO_tools.shmem_PLO_rmult_I(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            vec=vec,
            grp_prod=self._grp_prod,
            win_grp_prod=self._win_grp_prod,
            node_prod=self._node_prod,
            win_node_prod=self._win_node_prod,
            node_root=self.__shared_mem_mgr.node_root,
            tree_grp_comm=self.__shared_mem_mgr.tree_grp_comm,
            tree_grp_root_comm=self.__shared_mem_mgr.tree_grp_root_comm,
            node_comm=self.__shared_mem_mgr.node_comm,
            node_root_comm=self.__shared_mem_mgr.node_root_comm,
        )

        return self._node_prod.copy() if self.return_copy else self._node_prod

    def _mult_QU(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the matrix-vector product $Pv$ for linear
        polarization ($QU$) map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector $v$ of size `2*new_npix`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector of size `nsamples`
        """

        prod = np.zeros(self.nrows, dtype=self.dtype)

        PointingLO_tools.PLO_mult_QU(
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            prod=prod,
        )

        return prod

    def _rmult_QU(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        linear polarization ($QU$) map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting array of size `2*new_npix`
        """

        prod = np.zeros(self.ncols, dtype=self.dtype)

        PointingLO_tools.PLO_rmult_QU(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            prod=prod,
            comm=MPI_UTILS.comm,
        )

        return prod

    def _rmult_QU_shmem(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        linear polarization ($QU$) map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting array of size `2*new_npix`
        """

        if self.__shared_mem_mgr.tree_grp_rank == 0:
            self._grp_prod[:] = 0
        if self.__shared_mem_mgr.node_rank == 0:
            self._node_prod[:] = 0

        self._win_grp_prod.Fence(0)
        self._win_node_prod.Fence(0)

        PointingLO_tools.shmem_PLO_rmult_QU(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            grp_prod=self._grp_prod,
            win_grp_prod=self._win_grp_prod,
            node_prod=self._node_prod,
            win_node_prod=self._win_node_prod,
            node_root=self.__shared_mem_mgr.node_root,
            tree_grp_comm=self.__shared_mem_mgr.tree_grp_comm,
            tree_grp_root_comm=self.__shared_mem_mgr.tree_grp_root_comm,
            node_comm=self.__shared_mem_mgr.node_comm,
            node_root_comm=self.__shared_mem_mgr.node_root_comm,
        )

        return self._node_prod.copy() if self.return_copy else self._node_prod

    def _mult_IQU(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the matrix-vector product $Pv$ for temperature and
        linear polarization map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `3*new_npix`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector of size `nsamples`
        """

        prod = np.zeros(self.nrows, dtype=self.dtype)

        PointingLO_tools.PLO_mult_IQU(
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            prod=prod,
        )

        return prod

    def _rmult_IQU(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        temperature and linear polarization map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting array of size `3*new_npix`
        """

        prod = np.zeros(self.ncols, dtype=self.dtype)

        PointingLO_tools.PLO_rmult_IQU(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            prod=prod,
            comm=MPI_UTILS.comm,
        )

        return prod

    def _rmult_IQU_shmem(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the transposed matrix-vector product $P^T v$ for
        temperature and linear polarization map-making.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector of size `nsamples`

        Returns
        -------
        npt.NDArray[np.number]
            The resulting array of size `3*new_npix`
        """

        if self.__shared_mem_mgr.tree_grp_rank == 0:
            self._grp_prod[:] = 0
        if self.__shared_mem_mgr.node_rank == 0:
            self._node_prod[:] = 0

        self._win_grp_prod.Fence(0)
        self._win_node_prod.Fence(0)

        PointingLO_tools.shmem_PLO_rmult_IQU(
            new_npix=self.new_npix,
            nsamples=self.nrows,
            pointings=self.pointings,
            pointings_flag=self.pointings_flag,
            sin2phi=self.sin2phi,
            cos2phi=self.cos2phi,
            vec=vec,
            grp_prod=self._grp_prod,
            win_grp_prod=self._win_grp_prod,
            node_prod=self._node_prod,
            win_node_prod=self._win_node_prod,
            node_root=self.__shared_mem_mgr.node_root,
            tree_grp_comm=self.__shared_mem_mgr.tree_grp_comm,
            tree_grp_root_comm=self.__shared_mem_mgr.tree_grp_root_comm,
            node_comm=self.__shared_mem_mgr.node_comm,
            node_root_comm=self.__shared_mem_mgr.node_root_comm,
        )

        return self._node_prod.copy() if self.return_copy else self._node_prod

    @property
    def solver_type(self) -> SolverType:
        """The current map-making solver configuration.

        Returns
        -------
        SolverType
            The map-making solver type
        """
        return self.__solver_type

    @property
    def return_copy(self) -> bool:
        """Whether the transposed operator returns a copy of the shared memory
        buffer.

        Returns
        -------
        bool
            `True` if a copy is returned, `False` otherwise.
        """
        return self.__return_copy

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

solver_type: SolverType property

The current map-making solver configuration.

Returns:

Type	Description
SolverType	The map-making solver type

return_copy: bool property

Whether the transposed operator returns a copy of the shared memory buffer.

Returns:

Type	Description
bool	True if a copy is returned, False otherwise.

Methods:

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