linop package

linop 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/linop/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=False, **kwargs):
        self.__nargin = nargin
        self.__nargout = nargout
        self.__symmetric = symmetric
        self.__shape = (nargout, nargin)
        self.__dtype = kwargs.get("dtype", np.float64)
        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):
        """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

    @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

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/linop/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/linop/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/linop/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")

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

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

IdentityOperator

Bases: LinearOperator

Class representing the identity operator of size nargin.

Source code in brahmap/linop/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
        )

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/linop/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))
        for j in range(m):
            ej = np.zeros(m)
            ej[j] = 1.0
            H[:, j] = self * ej
        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)

    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/linop/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/linop/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/linop/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)

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/linop/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/linop/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 = issubclass(np.dtype(matrix.dtype).type, np.complex)

        symmetric = (
            np.all(matrix == matrix.conj().T)
            if iscomplex
            else np.all(matrix == matrix.T)
        )

        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,
        )

ShapeError

Bases: Exception

Exception class for handling shape mismatch errors.

Exception raised when defining a linear operator of the wrong shape or multiplying a linear operator with a vector of the wrong shape.

Source code in brahmap/linop/linop.py

class ShapeError(Exception):

    """
    Exception class for handling shape mismatch errors.

    Exception raised when defining a linear operator of the wrong shape or
    multiplying a linear operator with a vector of the wrong shape.

    """

    def __init__(self, value):
        super(ShapeError, self).__init__()
        self.value = value

    def __str__(self):
        return repr(self.value)

ZeroOperator

Bases: LinearOperator

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

Source code in brahmap/linop/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
        )

PysparseLinearOperator(A)

Return a linear operator from a Pysparse sparse matrix.

.. deprecated:: 0.6 Use :func:aslinearoperator instead.

Source code in brahmap/linop/linop.py

def PysparseLinearOperator(A):
    """
    Return a linear operator from a Pysparse sparse matrix.

    .. deprecated:: 0.6
        Use :func:`aslinearoperator` instead.

    """

    nargout, nargin = A.shape
    try:
        symmetric = A.issym
    except Exception:
        symmetric = A.isSymmetric()

    def matvec(x):
        if x.shape != (nargin,):
            msg = "Input has shape " + str(x.shape)
            msg += " instead of (%d,)" % nargin
            raise ValueError(msg)
        if hasattr(A, "__mul__"):
            return A * x
        Ax = np.empty(nargout)
        A.matvec(x, Ax)
        return Ax

    def rmatvec(y):
        if y.shape != (nargout,):
            msg = "Input has shape " + str(y.shape)
            msg += " instead of (%d,)" % nargout
            raise ValueError(msg)
        if hasattr(A, "__rmul__"):
            return y * A
        ATy = np.empty(nargin)
        A.rmatvec(y, ATy)
        return ATy

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

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/linop/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/linop/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/linop/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")

linop_from_ndarray(A)

Return a linear operator from a Numpy ndarray.

.. deprecated:: 0.4 Use :class:MatrixLinearOperator or :func:aslinearoperator instead.

Source code in brahmap/linop/linop.py

def linop_from_ndarray(A):
    """
    Return a linear operator from a Numpy `ndarray`.

    .. deprecated:: 0.4
        Use :class:`MatrixLinearOperator` or :func:`aslinearoperator` instead.

    """
    return LinearOperator(
        A.shape[1],
        A.shape[0],
        lambda v: np.dot(A, v),
        rmatvec=lambda u: np.dot(A.T, u),
        symmetric=False,
        dtype=A.dtype,
    )

blkop module

BlockDiagonalLinearOperator

Bases: LinearOperator

A block diagonal linear operator.

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

Source code in brahmap/linop/blkop.py

class BlockDiagonalLinearOperator(LinearOperator):

    """
    A block diagonal linear operator.

    Each block must be a linear operator.
    The blocks may 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")

        symmetric = reduce(lambda x, y: x and y, [blk.symmetric for blk in blocks])

        self._blocks = blocks

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

        nargins = [blk.shape[-1] for blk in blocks]
        log.debug("nargins = " + repr(nargins))

        nargouts = [blk.shape[0] for blk in blocks]
        log.debug("nargouts = " + repr(nargouts))

        nargins = np.array(nargins)
        nargouts = np.array(nargouts)
        nargin = sum(nargins)
        nargout = sum(nargouts)

        # Create blocks of transpose operator.
        blocksT = [blk.T for blk in blocks]

        def blk_matvec(x, blks):
            nx = len(x)
            nargins = [blk.shape[-1] for blk in blocks]
            nargin = sum(nargins)
            nargouts = [blk.shape[0] for blk in blocks]
            nargout = sum(nargouts)
            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.empty(int(nargout), dtype=result_type)

            nblks = len(blks)

            row_start = col_start = 0
            for blk in range(nblks):
                row_end = row_start + int(nargouts[blk])
                yout = y[row_start:row_end]

                col_end = col_start + int(nargins[blk])
                xin = x[col_start:col_end]

                B = blks[blk]
                yout[:] = B * xin

                col_start = col_end
                row_start = row_end

            return y

        blk_dtypes = [blk.dtype for blk in blocks]
        op_dtype = np.result_type(*blk_dtypes[0:10])

        super(BlockDiagonalLinearOperator, 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 diagonal operator."""
        return self._blocks

    def __getitem__(self, idx):
        blks = self._blocks[idx]
        if isinstance(idx, slice):
            return BlockDiagonalLinearOperator(blks, symmetric=self.symmetric)
        return blks

    def __setitem__(self, idx, ops):
        if not isinstance(ops, BaseLinearOperator):
            if isinstance(ops, list) or isinstance(ops, tuple):
                for op in ops:
                    if not isinstance(op, BaseLinearOperator):
                        msg = "Block operators can only contain"
                        msg += " linear operators"
                        raise ValueError(msg)
        self._blocks[idx] = ops

blocks property

The list of blocks defining the block diagonal operator.

BlockDiagonalPreconditioner

Bases: BlockDiagonalLinearOperator

An alias for BlockDiagonalLinearOperator.

Holds an additional solve method equivalent to __mul__.

Source code in brahmap/linop/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/linop/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/linop/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/linop/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/linop/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/linop/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/linop/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
        )