def projected(self, range_basis, source_basis, product=None, name=None):
        assert source_basis is None or source_basis in self.source
        assert range_basis is None or range_basis in self.range
        assert product is None or product.source == product.range == self.range

        if len(self.interpolation_dofs) == 0:
            return ZeroOperator(self.source, self.range, self.name).projected(range_basis, source_basis, product, name)
        elif not hasattr(self, 'restricted_operator') or source_basis is None:
            return super().projected(range_basis, source_basis, product, name)
        else:
            name = name or self.name + '_projected'

            if range_basis is not None:
                if product is None:
                    projected_collateral_basis = NumpyVectorSpace.make_array(self.collateral_basis.dot(range_basis),
                                                                             self.range.id)
                else:
                    projected_collateral_basis = NumpyVectorSpace.make_array(product.apply2(self.collateral_basis,
                                                                                            range_basis),
                                                                             self.range.id)
            else:
                projected_collateral_basis = self.collateral_basis

            return ProjectedEmpiciralInterpolatedOperator(self.restricted_operator, self.interpolation_matrix,
                                                          NumpyVectorSpace.make_array(source_basis.components(self.source_dofs)),
                                                          projected_collateral_basis, self.triangular,
                                                          self.source.id, None, name)

    def projected_to_subbasis(self, dim_range=None, dim_source=None, dim_collateral=None, name=None):
        assert dim_source is None or dim_source <= self.source.dim
        assert dim_range is None or dim_range <= self.range.dim
        assert dim_collateral is None or dim_collateral <= self.restricted_operator.range.dim
        if not isinstance(self.projected_collateral_basis.space, NumpyVectorSpace):
            raise NotImplementedError
        name = name or '{}_projected_to_subbasis'.format(self.name)

        interpolation_matrix = self.interpolation_matrix[:dim_collateral, :dim_collateral]

        if dim_collateral is not None:
            restricted_operator, source_dofs = self.restricted_operator.restricted(np.arange(dim_collateral))
        else:
            restricted_operator = self.restricted_operator

        old_pcb = self.projected_collateral_basis
        projected_collateral_basis = NumpyVectorSpace.make_array(old_pcb.data[:dim_collateral, :dim_range],
                                                                 old_pcb.space.id)

        old_sbd = self.source_basis_dofs
        source_basis_dofs = NumpyVectorSpace.make_array(old_sbd.data[:dim_source]) if dim_collateral is None \
            else NumpyVectorSpace.make_array(old_sbd.data[:dim_source, source_dofs])

        return ProjectedEmpiciralInterpolatedOperator(restricted_operator, interpolation_matrix,
                                                      source_basis_dofs, projected_collateral_basis, self.triangular,
                                                      self.source.id, solver_options=self.solver_options, name=name)

    def visualize(self, U, codim=2, **kwargs):
        """Visualize scalar data associated to the grid as a patch plot.

        Parameters
        ----------
        U
            |NumPy array| of the data to visualize. If `U.dim == 2 and len(U) > 1`, the
            data is visualized as a time series of plots. Alternatively, a tuple of
            |Numpy arrays| can be provided, in which case a subplot is created for
            each entry of the tuple. The lengths of all arrays have to agree.
        codim
            The codimension of the entities the data in `U` is attached to (either 0 or 2).
        kwargs
            See :func:`~pymor.gui.qt.visualize_patch`
        """
        from pymor.gui.qt import visualize_patch
        from pymor.vectorarrays.interfaces import VectorArrayInterface
        from pymor.vectorarrays.numpy import NumpyVectorSpace, NumpyVectorArray
        if isinstance(U, (np.ndarray, VectorArrayInterface)):
            U = (U,)
        assert all(isinstance(u, (np.ndarray, VectorArrayInterface)) for u in U)
        U = tuple(NumpyVectorSpace.make_array(u) if isinstance(u, np.ndarray) else
                  u if isinstance(u, NumpyVectorArray) else
                  NumpyVectorSpace.make_array(u.data)
                  for u in U)
        bounding_box = kwargs.pop('bounding_box', self.domain)
        visualize_patch(self, U, codim=codim, bounding_box=bounding_box, **kwargs)

def test_project_array():
    np.random.seed(123)
    U = NumpyVectorSpace.from_numpy(np.random.random((2, 10)))
    basis = NumpyVectorSpace.from_numpy(np.random.random((3, 10)))
    U_p = project_array(U, basis, orthonormal=False)
    onb = gram_schmidt(basis)
    U_p2 = project_array(U, onb, orthonormal=True)
    assert np.all(relative_error(U_p, U_p2) < 1e-10)

def test_identity_lincomb():
    space = NumpyVectorSpace(10)
    identity = IdentityOperator(space)
    ones = space.ones()
    idid = (identity + identity)
    assert almost_equal(ones * 2, idid.apply(ones))
    assert almost_equal(ones * 2, idid.apply_adjoint(ones))
    assert almost_equal(ones * 0.5, idid.apply_inverse(ones))
    assert almost_equal(ones * 0.5, idid.apply_inverse_adjoint(ones))

def test_ext(extension_alg):
    size = 5
    ident = np.identity(size)
    current = ident[0]
    for i in range(1, size):
        c = NumpyVectorSpace.from_data(current)
        n, _ = extension_alg(c, NumpyVectorSpace.from_data(ident[i]))
        assert np.allclose(n.data, ident[0:i+1])
        current = ident[0:i+1]

def test_axpy():
    x = NumpyVectorSpace.from_numpy(np.array([1.]))
    y = NumpyVectorSpace.from_numpy(np.array([1.]))
    y.axpy(1 + 1j, x)
    assert y.to_numpy()[0, 0] == 2 + 1j

    x = NumpyVectorSpace.from_numpy(np.array([1 + 1j]))
    y = NumpyVectorSpace.from_numpy(np.array([1.]))
    y.axpy(-1, x)
    assert y.to_numpy()[0, 0] == -1j

def test_project_array_with_product():
    np.random.seed(123)
    U = NumpyVectorSpace.from_numpy(np.random.random((1, 10)))
    basis = NumpyVectorSpace.from_numpy(np.random.random((3, 10)))
    product = np.random.random((10, 10))
    product = NumpyMatrixOperator(product.T.dot(product))
    U_p = project_array(U, basis, product=product, orthonormal=False)
    onb = gram_schmidt(basis, product=product)
    U_p2 = project_array(U, onb, product=product, orthonormal=True)
    assert np.all(relative_error(U_p, U_p2, product) < 1e-10)

 def __init__(self, array, adjoint=False, space_id=None, name=None):
     self._array = array.copy()
     if adjoint:
         self.source = array.space
         self.range = NumpyVectorSpace(len(array), space_id)
     else:
         self.source = NumpyVectorSpace(len(array), space_id)
         self.range = array.space
     self.adjoint = adjoint
     self.space_id = space_id
     self.name = name

 def __init__(self, array, transposed=False, space_id=None, name=None):
     self._array = array.copy()
     if transposed:
         self.source = array.space
         self.range = NumpyVectorSpace(len(array), space_id)
     else:
         self.source = NumpyVectorSpace(len(array), space_id)
         self.range = array.space
     self.transposed = transposed
     self.space_id = space_id
     self.name = name

def test_complex():
    np.random.seed(0)
    I = np.eye(5)
    A = np.random.randn(5, 5)
    B = np.random.randn(5, 5)
    C = np.random.randn(3, 5)

    Iop = NumpyMatrixOperator(I)
    Aop = NumpyMatrixOperator(A)
    Bop = NumpyMatrixOperator(B)
    Cva = NumpyVectorSpace.from_numpy(C)

    # lincombs
    assert not np.iscomplexobj((Iop * 1 + Bop * 1).assemble().matrix)
    assert not np.iscomplexobj((Aop * 1 + Bop * 1).assemble().matrix)
    assert np.iscomplexobj((Aop * (1+0j) + Bop * (1+0j)).assemble().matrix)
    assert np.iscomplexobj((Aop * 1j + Bop * 1).assemble().matrix)
    assert np.iscomplexobj((Bop * 1 + Aop * 1j).assemble().matrix)

    # apply_inverse
    assert not np.iscomplexobj(Aop.apply_inverse(Cva).to_numpy())
    assert np.iscomplexobj((Aop * 1j).apply_inverse(Cva).to_numpy())
    assert np.iscomplexobj((Aop * 1 + Bop * 1j).assemble().apply_inverse(Cva).to_numpy())
    assert np.iscomplexobj(Aop.apply_inverse(Cva * 1j).to_numpy())

    # append
    for rsrv in (0, 10):
        for o_ind in (slice(None), [0]):
            va = NumpyVectorSpace(5).empty(reserve=rsrv)
            va.append(Cva)
            D = np.random.randn(1, 5) + 1j * np.random.randn(1, 5)
            Dva = NumpyVectorSpace.from_numpy(D)

            assert not np.iscomplexobj(va.to_numpy())
            assert np.iscomplexobj(Dva.to_numpy())
            va.append(Dva[o_ind])
            assert np.iscomplexobj(va.to_numpy())

    # scal
    assert not np.iscomplexobj(Cva.to_numpy())
    assert np.iscomplexobj((Cva * 1j).to_numpy())
    assert np.iscomplexobj((Cva * (1 + 0j)).to_numpy())

    # axpy
    assert not np.iscomplexobj(Cva.to_numpy())
    Cva[0].axpy(1, Dva)
    assert np.iscomplexobj(Cva.to_numpy())

    Cva = NumpyVectorSpace.from_numpy(C)
    assert not np.iscomplexobj(Cva.to_numpy())
    Cva[0].axpy(1j, Dva)
    assert np.iscomplexobj(Cva.to_numpy())

 def __init__(self, mapping, adjoint_mapping=None, dim_source=1, dim_range=1, linear=False, parameter_type=None,
              source_id=None, range_id=None, solver_options=None, name=None):
     self.source = NumpyVectorSpace(dim_source, source_id)
     self.range = NumpyVectorSpace(dim_range, range_id)
     self.solver_options = solver_options
     self.name = name
     self._mapping = mapping
     self._adjoint_mapping = adjoint_mapping
     self.linear = linear
     if parameter_type is not None:
         self.build_parameter_type(parameter_type)
     self.source_id = source_id  # needed for with_
     self.range_id = range_id

 def save(self):
     if not config.HAVE_PYVTK:
         msg = QMessageBox(QMessageBox.Critical, 'Error', 'VTK output disabled. Pleas install pyvtk.')
         msg.exec_()
         return
     filename = QFileDialog.getSaveFileName(self, 'Save as vtk file')[0]
     base_name = filename.split('.vtu')[0].split('.vtk')[0].split('.pvd')[0]
     if base_name:
         if len(self.U) == 1:
             write_vtk(self.grid, NumpyVectorSpace.make_array(self.U[0]), base_name, codim=self.codim)
         else:
             for i, u in enumerate(self.U):
                 write_vtk(self.grid, NumpyVectorSpace.make_array(u), '{}-{}'.format(base_name, i),
                           codim=self.codim)

 def restricted(self, dofs):
     assert all(0 <= c < self.range.dim for c in dofs)
     if not self.transposed:
         restricted_value = NumpyVectorSpace.make_array(self._array.components(dofs))
         return VectorArrayOperator(restricted_value, False), np.arange(self.source.dim, dtype=np.int32)
     else:
         raise NotImplementedError

    def with_cb_dim(self, dim):
        assert dim <= self.restricted_operator.range.dim

        interpolation_matrix = self.interpolation_matrix[:dim, :dim]

        restricted_operator, source_dofs = self.restricted_operator.restricted(np.arange(dim))

        old_pcb = self.projected_collateral_basis
        projected_collateral_basis = NumpyVectorSpace.make_array(old_pcb.to_numpy()[:dim, :])

        old_sbd = self.source_basis_dofs
        source_basis_dofs = NumpyVectorSpace.make_array(old_sbd.to_numpy()[:, source_dofs])

        return ProjectedEmpiciralInterpolatedOperator(restricted_operator, interpolation_matrix,
                                                      source_basis_dofs, projected_collateral_basis, self.triangular,
                                                      solver_options=self.solver_options, name=self.name)

def test_to_matrix():
    np.random.seed(0)
    A = np.random.randn(2, 2)
    B = np.random.randn(3, 3)
    C = np.random.randn(3, 3)

    X = np.bmat([[np.eye(2) + A, np.zeros((2, 3))], [np.zeros((3, 2)), B.dot(C.T)]])

    C = sps.csc_matrix(C)

    Aop = NumpyMatrixOperator(A)
    Bop = NumpyMatrixOperator(B)
    Cop = NumpyMatrixOperator(C)

    Xop = BlockDiagonalOperator([LincombOperator([IdentityOperator(NumpyVectorSpace(2)), Aop],
                                                 [1, 1]), Concatenation(Bop, AdjointOperator(Cop))])

    assert np.allclose(X, to_matrix(Xop))
    assert np.allclose(X, to_matrix(Xop, format='csr').toarray())

    np.random.seed(0)
    V = np.random.randn(10, 2)
    Vva = NumpyVectorSpace.make_array(V.T)
    Vop = VectorArrayOperator(Vva)
    assert np.allclose(V, to_matrix(Vop))
    Vop = VectorArrayOperator(Vva, transposed=True)
    assert np.allclose(V, to_matrix(Vop).T)

    def action_ProjectedEmpiciralInterpolatedOperator(self, op):
        if not isinstance(op.projected_collateral_basis.space, NumpyVectorSpace):
            raise NotImplementedError

        restricted_operator = op.restricted_operator

        old_pcb = op.projected_collateral_basis
        projected_collateral_basis = NumpyVectorSpace.make_array(old_pcb.to_numpy()[:, :self.dim_range],
                                                                 old_pcb.space.id)

        old_sbd = op.source_basis_dofs
        source_basis_dofs = NumpyVectorSpace.make_array(old_sbd.to_numpy()[:self.dim_source])

        return ProjectedEmpiciralInterpolatedOperator(restricted_operator, op.interpolation_matrix,
                                                      source_basis_dofs, projected_collateral_basis, op.triangular,
                                                      op.source.id, solver_options=op.solver_options, name=op.name)

class ComponentProjection(OperatorBase):
    """|Operator| representing the projection of a |VectorArray| on some of its components.

    Parameters
    ----------
    components
        List or 1D |NumPy array| of the indices of the vector
        :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.components` that ar
        to be extracted by the operator.
    source
        Source |VectorSpace| of the operator.
    name
        Name of the operator.
    """

    linear = True

    def __init__(self, components, source, name=None):
        assert all(0 <= c < source.dim for c in components)
        self.components = np.array(components, dtype=np.int32)
        self.range = NumpyVectorSpace(len(components))
        self.source = source
        self.name = name

    def apply(self, U, mu=None):
        assert U in self.source
        return self.range.make_array(U.components(self.components))

    def restricted(self, dofs):
        assert all(0 <= c < self.range.dim for c in dofs)
        source_dofs = self.components[dofs]
        return IdentityOperator(NumpyVectorSpace(len(source_dofs))), source_dofs

class RestrictedDuneSpaceOperator(OperatorBase):

    linear = False

    def __init__(self, impl, range_dofs):
        self._impl = impl
        self._range_dofs = range_dofs
        self.source = NumpyVectorSpace(impl.dimSource)
        self.range = NumpyVectorSpace(len(range_dofs))
        self.name = 'DuneBurgersSpaceOperator_restricted'

        self._source_vec = dune_module.Vector(impl.dimSource, 0.)
        self._range_vec = dune_module.Vector(impl.dimRange, 0.)
        self._source_array = np.frombuffer(self._source_vec.buffer())
        self._range_array = np.frombuffer(self._range_vec.buffer())
        self.build_parameter_type({'exponent': tuple()}, local_global=True)

    def apply(self, U, ind=None, mu=None):
        assert U in self.source
        mu = self.parse_parameter(mu)
        exponent = float(mu['exponent'])

        U = U.data if ind is None else \
            U.data[ind] if hasattr(ind, '__len__') else \
            U.data[ind:ind + 1]
        R = self.range.zeros(len(U))
        R_array = R.data

        for i, u in enumerate(U):
            self._source_array[:] = u
            self._range_array[:] = 0
            self._impl.apply(self._source_vec, self._range_vec, exponent, 1.)
            R_array[i] = self._range_array[:][self._range_dofs]

        return R

    def action_EmpiricalInterpolatedOperator(self, op):
        range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
        if len(op.interpolation_dofs) == 0:
            return self.apply(ZeroOperator(op.range, op.source, op.name))
        elif not hasattr(op, 'restricted_operator') or source_basis is None:
            raise RuleNotMatchingError('Has no restricted operator or source_basis is None')
        else:
            if range_basis is not None:
                projected_collateral_basis = NumpyVectorSpace.make_array(op.collateral_basis.inner(range_basis,
                                                                                                   product),
                                                                         op.range.id)
            else:
                projected_collateral_basis = op.collateral_basis

            return ProjectedEmpiciralInterpolatedOperator(op.restricted_operator, op.interpolation_matrix,
                                                          NumpyVectorSpace.make_array(source_basis.dofs(op.source_dofs)),
                                                          projected_collateral_basis, op.triangular,
                                                          op.source.id, None, op.name)