Base Operators

linop sub-module

BaseLinearOperator

Bases: object

Base class defining the common interface shared by all linear operators.

A linear operator is a linear mapping x -> A(x) such that the size of the input vector x is nargin and the size of the output is nargout. It can be visualized as a matrix of shape (nargout, nargin). Its type is any valid Numpy dtype. By default, it has dtype numpy.float but this can be changed to, e.g., numpy.complex via the dtype keyword argument and attribute.

A logger may be attached to the linear operator via the logger keyword argument.

Source code in brahmap/base/linop.py

class BaseLinearOperator(object):
    """
    Base class defining the common interface shared by all linear operators.

    A linear operator is a linear mapping x -> A(x) such that the size of the
    input vector x is `nargin` and the size of the output is `nargout`. It can
    be visualized as a matrix of shape (`nargout`, `nargin`). Its type is any
    valid Numpy `dtype`. By default, it has `dtype` `numpy.float` but this can
    be changed to, e.g., `numpy.complex` via the `dtype` keyword argument and
    attribute.

    A logger may be attached to the linear operator via the `logger` keyword
    argument.

    """

    def __init__(
        self, nargin, nargout, symmetric: bool = False, dtype=np.float64, **kwargs
    ):
        self.__nargin = nargin
        self.__nargout = nargout
        self.__symmetric = symmetric
        self.__shape = (nargout, nargin)
        self.dtype = dtype
        self._nMatvec = 0

        # Log activity.
        self.logger = kwargs.get("logger", null_log)
        self.logger.info("New linear operator with shape " + str(self.shape))
        return

    @property
    def nargin(self):
        """The size of an input vector."""
        return self.__nargin

    @property
    def nargout(self):
        """The size of an output vector."""
        return self.__nargout

    @property
    def symmetric(self) -> bool:
        """Indicate whether the operator is symmetric or not."""
        return self.__symmetric

    @property
    def shape(self):
        """The shape of the operator."""
        return self.__shape

    @property
    def dtype(self):
        """The data type of the operator."""
        return self.__dtype

    @dtype.setter
    def dtype(self, dtype):
        self.__dtype = dtype

    @property
    def nMatvec(self):
        """The number of products with vectors computed so far."""
        return self._nMatvec

    def reset_counters(self):
        """Reset operator/vector product counter to zero."""
        self._nMatvec = 0

    def __call__(self, *args, **kwargs):
        # An alias for __mul__.
        return self.__mul__(*args, **kwargs)

    def __mul__(self, x):
        raise NotImplementedError("Please subclass to implement __mul__.")

    def __repr__(self):
        if self.symmetric:
            s = "Symmetric"
        else:
            s = "Unsymmetric"
        s += " <" + self.__class__.__name__ + ">"
        s += " of type %s" % self.dtype
        s += " with shape (%d,%d)" % (self.nargout, self.nargin)
        return s

    def dot(self, x):
        """Numpy-like dot() method."""
        return self.__mul__(x)

dtype property writable

The data type of the operator.

nMatvec property

The number of products with vectors computed so far.

nargin property

The size of an input vector.

nargout property

The size of an output vector.

shape property

The shape of the operator.

symmetric property

Indicate whether the operator is symmetric or not.

dot(x)

Numpy-like dot() method.

Source code in brahmap/base/linop.py

def dot(self, x):
    """Numpy-like dot() method."""
    return self.__mul__(x)

reset_counters()

Reset operator/vector product counter to zero.

Source code in brahmap/base/linop.py

def reset_counters(self):
    """Reset operator/vector product counter to zero."""
    self._nMatvec = 0

DiagonalOperator

Bases: LinearOperator

Class representing a diagonal operator.

A diagonal linear operator defined by its diagonal diag (a Numpy array.) The type must be specified in the diag argument, e.g., np.ones(5, dtype=np.complex) or np.ones(5).astype(np.complex).

Source code in brahmap/base/linop.py

class DiagonalOperator(LinearOperator):
    """
    Class representing a diagonal operator.

    A diagonal linear operator defined by its diagonal `diag` (a Numpy array.)
    The type must be specified in the `diag` argument, e.g.,
    `np.ones(5, dtype=np.complex)` or `np.ones(5).astype(np.complex)`.

    """

    def __init__(self, diag, **kwargs):
        if "symmetric" in kwargs:
            kwargs.pop("symmetric")
        if "matvec" in kwargs:
            kwargs.pop("matvec")
        if "dtype" in kwargs:
            kwargs.pop("dtype")

        self.diag = np.asarray(diag)
        if self.diag.ndim != 1:
            msg = "diag array must be 1-d"
            raise ValueError(msg)

        super(DiagonalOperator, self).__init__(
            self.diag.shape[0],
            self.diag.shape[0],
            symmetric=True,
            matvec=lambda x: self.diag * x,
            dtype=self.diag.dtype,
            **kwargs,
        )

IdentityOperator

Bases: LinearOperator

Class representing the identity operator of size nargin.

Source code in brahmap/base/linop.py

class IdentityOperator(LinearOperator):
    """Class representing the identity operator of size `nargin`."""

    def __init__(self, nargin, **kwargs):
        if "symmetric" in kwargs:
            kwargs.pop("symmetric")
        if "matvec" in kwargs:
            kwargs.pop("matvec")

        super(IdentityOperator, self).__init__(
            nargin, nargin, symmetric=True, matvec=lambda x: x, **kwargs
        )

InverseLO

Bases: LinearOperator

Construct the inverse operator of a matrix :math:A, as a linear operator.

Parameters

• A : {linear operator} the linear operator of the linear system to invert;
• method : {function } the method to compute A^-1 (see below);
• P : {linear operator } (optional) the preconditioner for the computation of the inverse operator.

Source code in brahmap/base/linop.py

class InverseLO(LinearOperator):
    r"""
    Construct the inverse operator of a matrix :math:`A`, as a linear operator.

    **Parameters**

    - ``A`` : {linear operator}
        the linear operator of the linear system to invert;
    - ``method`` : {function }
        the method to compute ``A^-1`` (see below);
    - ``P`` : {linear operator } (optional)
        the preconditioner for the computation of the inverse operator.

    """

    def mult(self, x):
        r"""
        It returns  :math:`y=A^{-1}x` by solving the linear system :math:`Ay=x`
        with a certain :mod:`scipy` routine (e.g. :func:`scipy.sparse.linalg.cg`)
        defined above as ``method``.
        """

        y, info = self.method(self.A, x, M=self.preconditioner)
        self.isconverged(info)
        return y

    def isconverged(self, info):
        r"""
        It returns a Boolean value  depending on the
        exit status of the solver.

        **Parameters**

        - ``info`` : {int}
            output of the solver method (usually :func:`scipy.sparse.cg`).
        """
        self.__converged = info
        if info == 0:
            return True
        else:
            return False

    def __init__(self, A, method=None, preconditioner=None):
        super(InverseLO, self).__init__(
            nargin=A.shape[0], nargout=A.shape[1], matvec=self.mult, symmetric=True
        )
        self.A = A
        self.__method = method
        self.__preconditioner = preconditioner
        self.__converged = None

    @property
    def method(self):
        r"""
        The method to compute the inverse of A. \
        It can be any :mod:`scipy.sparse.linalg` solver, namely :func:`scipy.sparse.linalg.cg`,
        :func:`scipy.sparse.linalg.bicg`, etc.

        """
        return self.__method

    @property
    def converged(self):
        r"""
        provides convergence information:

        - 0 : successful exit;
        - >0 : convergence to tolerance not achieved, number of iterations;
        - <0 : illegal input or breakdown.

        """
        return self.__converged

    @property
    def preconditioner(self):
        """
        Preconditioner for the solver.
        """
        return self.__preconditioner

converged property

provides convergence information:

• 0 : successful exit;

• 0 : convergence to tolerance not achieved, number of iterations;

• <0 : illegal input or breakdown.

method property

The method to compute the inverse of A. \ It can be any :mod:scipy.sparse.linalg solver, namely :func:scipy.sparse.linalg.cg, :func:scipy.sparse.linalg.bicg, etc.

preconditioner property

Preconditioner for the solver.

isconverged(info)

It returns a Boolean value depending on the exit status of the solver.

Parameters

• info : {int} output of the solver method (usually :func:scipy.sparse.cg).

Source code in brahmap/base/linop.py

def isconverged(self, info):
    r"""
    It returns a Boolean value  depending on the
    exit status of the solver.

    **Parameters**

    - ``info`` : {int}
        output of the solver method (usually :func:`scipy.sparse.cg`).
    """
    self.__converged = info
    if info == 0:
        return True
    else:
        return False

mult(x)

It returns :math:y=A^{-1}x by solving the linear system :math:Ay=x with a certain :mod:scipy routine (e.g. :func:scipy.sparse.linalg.cg) defined above as method.

Source code in brahmap/base/linop.py

def mult(self, x):
    r"""
    It returns  :math:`y=A^{-1}x` by solving the linear system :math:`Ay=x`
    with a certain :mod:`scipy` routine (e.g. :func:`scipy.sparse.linalg.cg`)
    defined above as ``method``.
    """

    y, info = self.method(self.A, x, M=self.preconditioner)
    self.isconverged(info)
    return y

LinearOperator

Bases: BaseLinearOperator

Generic linear operator class.

A linear operator constructed from a matvec and (possibly) a rmatvec function. If symmetric is True, rmatvec is ignored. All other keyword arguments are passed directly to the superclass.

Source code in brahmap/base/linop.py

class LinearOperator(BaseLinearOperator):
    """
    Generic linear operator class.

    A linear operator constructed from a `matvec` and (possibly) a
    `rmatvec` function. If `symmetric` is `True`, `rmatvec` is
    ignored. All other keyword arguments are passed directly to the superclass.

    """

    def __init__(self, nargin, nargout, matvec, rmatvec=None, **kwargs):
        super(LinearOperator, self).__init__(nargin, nargout, **kwargs)
        adjoint_of = kwargs.get("adjoint_of", None) or kwargs.get("transpose_of", None)
        rmatvec = rmatvec or kwargs.get("matvec_transp", None)

        self.__matvec = matvec

        if self.symmetric:
            self.__H = self
        else:
            if adjoint_of is None:
                if rmatvec is not None:
                    # Create 'pointer' to transpose operator.
                    self.__H = LinearOperator(
                        nargout,
                        nargin,
                        matvec=rmatvec,
                        rmatvec=matvec,
                        adjoint_of=self,
                        **kwargs,
                    )
                else:
                    self.__H = None
            else:
                # Use operator supplied as transpose operator.
                if isinstance(adjoint_of, BaseLinearOperator):
                    self.__H = adjoint_of
                else:
                    msg = "kwarg adjoint_of / transpose_of must be of type LinearOperator."
                    msg += " Got " + str(adjoint_of.__class__)
                    raise ValueError(msg)

    @property
    def T(self):
        """The transpose operator.

        .. note:: this is an alias to the adjoint operator

        """
        return self.__H

    @property
    def H(self):
        """The adjoint operator."""
        return self.__H

    def matvec(self, x):
        """
        Matrix-vector multiplication.

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

        """
        x = np.asanyarray(x)
        M, N = self.shape

        # check input data consistency
        N = int(N)
        try:
            x = x.reshape(N)
        except ValueError:
            msg = "input array size incompatible with operator dimensions"
            raise ValueError(msg)

        y = self.__matvec(x)

        # check output data consistency
        M = int(M)
        try:
            y = y.reshape(M)
        except ValueError:
            msg = "output array size incompatible with operator dimensions"
            raise ValueError(msg)

        return y

    def to_array(self):
        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

    def __mul_scalar(self, x):
        """Product between a linear operator and a scalar."""
        result_type = np.result_type(self.dtype, type(x))

        if x != 0:

            def matvec(y):
                return x * (self(y))

            def rmatvec(y):
                return x * (self.H(y))

            return LinearOperator(
                self.nargin,
                self.nargout,
                symmetric=self.symmetric,
                matvec=matvec,
                rmatvec=rmatvec,
                dtype=result_type,
            )
        else:
            return ZeroOperator(self.nargin, self.nargout, dtype=result_type)

    def __mul_linop(self, op):
        """Product between two linear operators."""
        if self.nargin != op.nargout:
            raise ShapeError("Cannot multiply operators together")

        def matvec(x):
            return self(op(x))

        def rmatvec(x):
            return op.T(self.H(x))

        result_type = np.result_type(self.dtype, op.dtype)

        return LinearOperator(
            op.nargin,
            self.nargout,
            symmetric=False,  # Generally.
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __mul_vector(self, x):
        """Product between a linear operator and a vector."""
        self._nMatvec += 1
        result_type = np.result_type(self.dtype, x.dtype)
        return self.matvec(x).astype(result_type, copy=False)

    def __mul__(self, x):
        if np.isscalar(x):
            return self.__mul_scalar(x)
        elif isinstance(x, BaseLinearOperator):
            return self.__mul_linop(x)
        elif isinstance(x, np.ndarray):
            return self.__mul_vector(x)
        else:
            raise ValueError("Cannot multiply")

    def __rmul__(self, x):
        if np.isscalar(x):
            return self.__mul__(x)
        raise ValueError("Cannot multiply")

    def __add__(self, other):
        if not isinstance(other, BaseLinearOperator):
            raise ValueError("Cannot add")
        if self.shape != other.shape:
            raise ShapeError("Cannot add")

        def matvec(x):
            return self(x) + other(x)

        def rmatvec(x):
            return self.H(x) + other.T(x)

        result_type = np.result_type(self.dtype, other.dtype)

        return LinearOperator(
            self.nargin,
            self.nargout,
            symmetric=self.symmetric and other.symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __neg__(self):
        return self * (-1)

    def __sub__(self, other):
        if not isinstance(other, BaseLinearOperator):
            raise ValueError("Cannot add")
        if self.shape != other.shape:
            raise ShapeError("Cannot add")

        def matvec(x):
            return self(x) - other(x)

        def rmatvec(x):
            return self.H(x) - other.T(x)

        result_type = np.result_type(self.dtype, other.dtype)

        return LinearOperator(
            self.nargin,
            self.nargout,
            symmetric=self.symmetric and other.symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __truediv__(self, other):
        if np.isscalar(other):
            return self * (1 / other)
        else:
            raise ValueError("Cannot divide")

    def __pow__(self, other):
        if not isinstance(other, int):
            raise ValueError("Can only raise to integer power")
        if other < 0:
            raise ValueError("Can only raise to nonnegative power")
        if self.nargin != self.nargout:
            raise ShapeError("Can only raise square operators to a power")
        if other == 0:
            return IdentityOperator(self.nargin)
        if other == 1:
            return self
        return self * self ** (other - 1)

H property

The adjoint operator.

T property

The transpose operator.

.. note:: this is an alias to the adjoint operator

__mul_linop(op)

Product between two linear operators.

Source code in brahmap/base/linop.py

def __mul_linop(self, op):
    """Product between two linear operators."""
    if self.nargin != op.nargout:
        raise ShapeError("Cannot multiply operators together")

    def matvec(x):
        return self(op(x))

    def rmatvec(x):
        return op.T(self.H(x))

    result_type = np.result_type(self.dtype, op.dtype)

    return LinearOperator(
        op.nargin,
        self.nargout,
        symmetric=False,  # Generally.
        matvec=matvec,
        rmatvec=rmatvec,
        dtype=result_type,
    )

__mul_scalar(x)

Product between a linear operator and a scalar.

Source code in brahmap/base/linop.py

def __mul_scalar(self, x):
    """Product between a linear operator and a scalar."""
    result_type = np.result_type(self.dtype, type(x))

    if x != 0:

        def matvec(y):
            return x * (self(y))

        def rmatvec(y):
            return x * (self.H(y))

        return LinearOperator(
            self.nargin,
            self.nargout,
            symmetric=self.symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )
    else:
        return ZeroOperator(self.nargin, self.nargout, dtype=result_type)

__mul_vector(x)

Product between a linear operator and a vector.

Source code in brahmap/base/linop.py

def __mul_vector(self, x):
    """Product between a linear operator and a vector."""
    self._nMatvec += 1
    result_type = np.result_type(self.dtype, x.dtype)
    return self.matvec(x).astype(result_type, copy=False)

matvec(x)

Matrix-vector multiplication.

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

Source code in brahmap/base/linop.py

def matvec(self, x):
    """
    Matrix-vector multiplication.

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

    """
    x = np.asanyarray(x)
    M, N = self.shape

    # check input data consistency
    N = int(N)
    try:
        x = x.reshape(N)
    except ValueError:
        msg = "input array size incompatible with operator dimensions"
        raise ValueError(msg)

    y = self.__matvec(x)

    # check output data consistency
    M = int(M)
    try:
        y = y.reshape(M)
    except ValueError:
        msg = "output array size incompatible with operator dimensions"
        raise ValueError(msg)

    return y

MatrixLinearOperator

Bases: LinearOperator

Class representing a matrix operator.

A linear operator wrapping the multiplication with a matrix and its transpose (real) or conjugate transpose (complex). The operator's dtype is the same as the specified matrix argument.

.. versionadded:: 0.3

Source code in brahmap/base/linop.py

class MatrixLinearOperator(LinearOperator):
    """
    Class representing a matrix operator.

    A linear operator wrapping the multiplication with a matrix and its
    transpose (real) or conjugate transpose (complex). The operator's dtype
    is the same as the specified `matrix` argument.

    .. versionadded:: 0.3

    """

    def __init__(self, matrix, **kwargs):
        if "symmetric" in kwargs:
            kwargs.pop("symmetric")
        if "matvec" in kwargs:
            kwargs.pop("matvec")
        if "dtype" in kwargs:
            kwargs.pop("dtype")

        if not hasattr(matrix, "shape"):
            matrix = np.asanyarray(matrix)

        if matrix.ndim != 2:
            msg = "matrix must be 2-d (shape can be [M, N], [M, 1] or [1, N])"
            raise ValueError(msg)

        matvec = matrix.dot
        iscomplex = np.iscomplexobj(matrix)

        if matrix.shape[0] == matrix.shape[1]:
            symmetric = np.all(matrix == matrix.conj().T)
        else:
            symmetric = False

        if not symmetric:
            rmatvec = matrix.conj().T.dot if iscomplex else matrix.T.dot
        else:
            rmatvec = None

        super(MatrixLinearOperator, self).__init__(
            matrix.shape[1],
            matrix.shape[0],
            symmetric=symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=matrix.dtype,
            **kwargs,
        )

ZeroOperator

Bases: LinearOperator

Class representing the zero operator of shape nargout-by-nargin.

Source code in brahmap/base/linop.py

class ZeroOperator(LinearOperator):
    """Class representing the zero operator of shape `nargout`-by-`nargin`."""

    def __init__(self, nargin, nargout, **kwargs):
        if "matvec" in kwargs:
            kwargs.pop("matvec")
        if "rmatvec" in kwargs:
            kwargs.pop("rmatvec")

        def matvec(x):
            if x.shape != (nargin,):
                msg = "Input has shape " + str(x.shape)
                msg += " instead of (%d,)" % self.nargin
                raise ValueError(msg)
            return np.zeros(nargout)

        def rmatvec(x):
            if x.shape != (nargout,):
                msg = "Input has shape " + str(x.shape)
                msg += " instead of (%d,)" % self.nargout
                raise ValueError(msg)
            return np.zeros(nargin)

        super(ZeroOperator, self).__init__(
            nargin, nargout, matvec=matvec, rmatvec=rmatvec, **kwargs
        )

ReducedLinearOperator(op, row_indices, col_indices)

Implement reduction of a linear operator (non symmetrical).

Reduce a linear operator by limiting its input to col_indices and its output to row_indices.

Source code in brahmap/base/linop.py

def ReducedLinearOperator(op, row_indices, col_indices):
    """
    Implement reduction of a linear operator (non symmetrical).

    Reduce a linear operator by limiting its input to `col_indices` and its
    output to `row_indices`.

    """

    nargin, nargout = len(col_indices), len(row_indices)
    m, n = op.shape  # Shape of non-reduced operator.

    def matvec(x):
        z = np.zeros(n, dtype=x.dtype)
        z[col_indices] = x[:]
        y = op * z
        return y[row_indices]

    def rmatvec(x):
        z = np.zeros(m, dtype=x.dtype)
        z[row_indices] = x[:]
        y = op.H * z
        return y[col_indices]

    return LinearOperator(
        nargin, nargout, matvec=matvec, symmetric=False, rmatvec=rmatvec
    )

SymmetricallyReducedLinearOperator(op, indices)

Implement reduction of a linear operator (symmetrical).

Reduce a linear operator symmetrically by reducing boths its input and output to indices.

Source code in brahmap/base/linop.py

def SymmetricallyReducedLinearOperator(op, indices):
    """
    Implement reduction of a linear operator (symmetrical).

    Reduce a linear operator symmetrically by reducing boths its input and
    output to `indices`.

    """

    nargin = len(indices)
    m, n = op.shape  # Shape of non-reduced operator.

    def matvec(x):
        z = np.zeros(n, dtype=x.dtype)
        z[indices] = x[:]
        y = op * z
        return y[indices]

    def rmatvec(x):
        z = np.zeros(m, dtype=x.dtype)
        z[indices] = x[:]
        y = op * z
        return y[indices]

    return LinearOperator(
        nargin, nargin, matvec=matvec, symmetric=op.symmetric, rmatvec=rmatvec
    )

aslinearoperator(A)

Return A as a LinearOperator.

'A' may be any of the following types: - linop.LinearOperator - scipy.LinearOperator - ndarray - matrix - sparse matrix (e.g. csr_matrix, lil_matrix, etc.) - any object with .shape and .matvec attributes

See the :class:LinearOperator documentation for additonal information.

.. versionadded:: 0.4

Source code in brahmap/base/linop.py

def aslinearoperator(A):
    """Return A as a LinearOperator.

    'A' may be any of the following types:
    - linop.LinearOperator
    - scipy.LinearOperator
    - ndarray
    - matrix
    - sparse matrix (e.g. csr_matrix, lil_matrix, etc.)
    - any object with .shape and .matvec attributes

    See the :class:`LinearOperator` documentation for additonal information.

    .. versionadded:: 0.4

    """
    if isinstance(A, LinearOperator):
        return A

    try:
        import numpy as np

        if isinstance(A, np.ndarray) or isinstance(A, np.matrix):
            return MatrixLinearOperator(A)
    except ImportError:
        pass

    try:
        import scipy.sparse as ssp

        if ssp.isspmatrix(A):
            return MatrixLinearOperator(A)
    except ImportError:
        pass

    if hasattr(A, "shape"):
        nargout, nargin = A.shape
        matvec = None
        rmatvec = None
        dtype = None
        symmetric = False
        if hasattr(A, "matvec"):
            matvec = A.matvec
            if hasattr(A, "rmatvec"):
                rmatvec = A.rmatvec
            elif hasattr(A, "matvec_transp"):
                rmatvec = A.matvec_transp
            if hasattr(A, "dtype"):
                dtype = A.dtype
            if hasattr(A, "symmetric"):
                symmetric = A.symmetric
        elif hasattr(A, "__mul__"):

            def matvec(x):
                return A * x

            if hasattr(A, "__rmul__"):

                def rmatvec(x):
                    return x * A

            if hasattr(A, "dtype"):
                dtype = A.dtype
            try:
                symmetric = A.isSymmetric()
            except Exception:
                symmetric = False
        return LinearOperator(
            nargin,
            nargout,
            symmetric=symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=dtype,
        )
    else:
        raise TypeError("unsupported object type")

blkop sub-module

BlockDiagonalPreconditioner

Bases: BlockDiagonalLinearOperator

An alias for BlockDiagonalLinearOperator.

Holds an additional solve method equivalent to __mul__.

Source code in brahmap/base/blkop.py

class BlockDiagonalPreconditioner(BlockDiagonalLinearOperator):
    """
    An alias for ``BlockDiagonalLinearOperator``.

    Holds an additional ``solve`` method equivalent to ``__mul__``.

    """

    def solve(self, x):
        """An alias to __call__."""
        return self.__call__(x)

solve(x)

An alias to call.

Source code in brahmap/base/blkop.py

def solve(self, x):
    """An alias to __call__."""
    return self.__call__(x)

BlockHorizontalLinearOperator

Bases: BlockLinearOperator

A block horizontal linear operator.

Each block must be a linear operator. The blocks must be specified as one list, e.g., [A, B, C].

Source code in brahmap/base/blkop.py

class BlockHorizontalLinearOperator(BlockLinearOperator):
    """
    A block horizontal linear operator.

    Each block must be a linear operator.
    The blocks must be specified as one list, e.g., `[A, B, C]`.

    """

    def __init__(self, blocks, **kwargs):
        try:
            for block in blocks:
                __ = block.shape
        except (TypeError, AttributeError):
            raise ValueError("blocks should be a flattened list of operators")

        blocks = [[blk for blk in blocks]]

        super(BlockHorizontalLinearOperator, self).__init__(
            blocks=blocks, symmetric=False, **kwargs
        )

BlockLinearOperator

Bases: LinearOperator

A linear operator defined by blocks. Each block must be a linear operator.

blocks should be a list of lists describing the blocks row-wise. If there is only one block row, it should be specified as [[b1, b2, ..., bn]], not as [b1, b2, ..., bn].

If the overall linear operator is symmetric, only its upper triangle need be specified, e.g., [[A,B,C], [D,E], [F]], and the blocks on the diagonal must be square and symmetric.

Source code in brahmap/base/blkop.py

class BlockLinearOperator(LinearOperator):
    """
    A linear operator defined by blocks. Each block must be a linear operator.

    `blocks` should be a list of lists describing the blocks row-wise.
    If there is only one block row, it should be specified as
    `[[b1, b2, ..., bn]]`, not as `[b1, b2, ..., bn]`.

    If the overall linear operator is symmetric, only its upper triangle
    need be specified, e.g., `[[A,B,C], [D,E], [F]]`, and the blocks on the
    diagonal must be square and symmetric.

    """

    def __init__(self, blocks, symmetric=False, **kwargs):
        # If building a symmetric operator, fill in the blanks.
        # They're just references to existing objects.
        try:
            for block_row in blocks:
                for block_col in block_row:
                    __ = block_col.shape
        except (TypeError, AttributeError):
            raise ValueError("blocks should be a nested list of operators")

        if symmetric:
            nrow = len(blocks)
            ncol = len(blocks[0])
            if nrow != ncol:
                raise ShapeError("Inconsistent shape.")

            for block_row in blocks:
                if not block_row[0].symmetric:
                    raise ValueError("Blocks on diagonal must be symmetric.")

            self._blocks = blocks[:]
            for i in range(1, nrow):
                for j in range(i - 1, -1, -1):
                    self._blocks[i].insert(0, self._blocks[j][i].T)

        else:
            self._blocks = blocks

        log = kwargs.get("logger", null_log)
        log.debug("Building new BlockLinearOperator")

        nargins = [[blk.shape[-1] for blk in row] for row in self._blocks]
        log.debug("nargins = " + repr(nargins))
        nargins_by_row = [nargin[0] for nargin in nargins]
        if min(nargins_by_row) != max(nargins_by_row):
            raise ShapeError("Inconsistent block shapes")

        nargouts = [[blk.shape[0] for blk in row] for row in self._blocks]
        log.debug("nargouts = " + repr(nargouts))
        for row in nargouts:
            if min(row) != max(row):
                raise ShapeError("Inconsistent block shapes")

        nargin = sum(nargins[0])
        nargout = sum([out[0] for out in nargouts])

        # Create blocks of transpose operator.
        blocksT = list(map(lambda *row: [blk.T for blk in row], *self._blocks))

        def blk_matvec(x, blks):
            nargins = [[blk.shape[-1] for blk in blkrow] for blkrow in blks]
            nargouts = [[blk.shape[0] for blk in blkrow] for blkrow in blks]
            nargin = sum(nargins[0])
            nargout = sum([out[0] for out in nargouts])
            nx = len(x)
            self.logger.debug("Multiplying with a vector of size %d" % nx)
            self.logger.debug("nargin=%d, nargout=%d" % (nargin, nargout))
            if nx != nargin:
                raise ShapeError("Multiplying with vector of wrong shape.")

            result_type = np.result_type(self.dtype, x.dtype)
            y = np.zeros(nargout, dtype=result_type)

            nblk_row = len(blks)
            nblk_col = len(blks[0])

            row_start = col_start = 0
            for row in range(nblk_row):
                row_end = row_start + nargouts[row][0]
                yout = y[row_start:row_end]
                for col in range(nblk_col):
                    col_end = col_start + nargins[0][col]
                    xin = x[col_start:col_end]
                    B = blks[row][col]
                    yout[:] += B * xin
                    col_start = col_end
                row_start = row_end
                col_start = 0

            return y

        flat_blocks = list(itertools.chain(*blocks))
        blk_dtypes = [blk.dtype for blk in flat_blocks]
        op_dtype = np.result_type(*blk_dtypes)

        super(BlockLinearOperator, self).__init__(
            nargin,
            nargout,
            symmetric=symmetric,
            matvec=lambda x: blk_matvec(x, self._blocks),
            rmatvec=lambda x: blk_matvec(x, blocksT),
            dtype=op_dtype,
            **kwargs,
        )

        self.H._blocks = blocksT

    @property
    def blocks(self):
        """The list of blocks defining the block operator."""
        return self._blocks

    def __getitem__(self, indices):
        blks = np.matrix(self._blocks, dtype=object)[indices]
        # If indexing narrowed it down to a single block, return it.
        if isinstance(blks, BaseLinearOperator):
            return blks
        # Otherwise, we have a matrix of blocks.
        return BlockLinearOperator(blks.tolist(), symmetric=False)

    def __contains__(self, op):
        flat_blocks = list(itertools.chain(*self.blocks))
        return op in flat_blocks

    def __iter__(self):
        for block in self._blocks:
            yield block

blocks property

The list of blocks defining the block operator.

BlockPreconditioner

Bases: BlockLinearOperator

An alias for BlockLinearOperator.

Holds an additional solve method equivalent to __mul__.

Source code in brahmap/base/blkop.py

class BlockPreconditioner(BlockLinearOperator):
    """An alias for ``BlockLinearOperator``.

    Holds an additional ``solve`` method equivalent to ``__mul__``.

    """

    def solve(self, x):
        """An alias to __call__."""
        return self.__call__(x)

solve(x)

An alias to call.

Source code in brahmap/base/blkop.py

def solve(self, x):
    """An alias to __call__."""
    return self.__call__(x)

BlockVerticalLinearOperator

Bases: BlockLinearOperator

A block vertical linear operator.

Each block must be a linear operator. The blocks must be specified as one list, e.g., [A, B, C].

Source code in brahmap/base/blkop.py

class BlockVerticalLinearOperator(BlockLinearOperator):
    """
    A block vertical linear operator.

    Each block must be a linear operator.
    The blocks must be specified as one list, e.g., `[A, B, C]`.

    """

    def __init__(self, blocks, **kwargs):
        try:
            for block in blocks:
                __ = block.shape
        except (TypeError, AttributeError):
            raise ValueError("blocks should be a flattened list of operators")

        blocks = [[blk] for blk in blocks]

        super(BlockVerticalLinearOperator, self).__init__(
            blocks=blocks, symmetric=False, **kwargs
        )