<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def matrix_representation(op):
 """Return a matrix representation of a linear operator. Parameters op : `Operator` The linear operator of which one wants a matrix representation. If the domain or range is a `ProductSpace`, it must be a power-space. Returns ------- matrix : `numpy.ndarray` The matrix representation of the operator. The shape will be ``op.domain.shape + op.range.shape`` and the dtype is the promoted (greatest) dtype of the domain and range. Examples -------- Approximate a matrix on its own: array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) It also works with `ProductSpace`'s and higher dimensional `TensorSpace`'s. In this case, the returned "matrix" will also be higher dimensional: True Since the "matrix" is now higher dimensional, we need to use e.g. `numpy.tensordot` if we want to compute with the matrix representation: ProductSpace(uniform_discr([ 0., 0.], [ 2., 2.], (2, 2)), 2).element([ <BLANKLINE> [[ 2. , 2. ], [-2.75, -6.75]], <BLANKLINE> [[ 4. , -4.75], [ 4. , -6.75]] ]) array([[[ 2. , 2. ], [-2.75, -6.75]], <BLANKLINE> [[ 4. , -4.75], [ 4. , -6.75]]]) Notes The algorithm works by letting the operator act on all unit vectors, and stacking the output as a matrix. """


    if not op.is_linear:
        raise ValueError('the operator is not linear')

    if not (isinstance(op.domain, TensorSpace) or
            (isinstance(op.domain, ProductSpace) and
             op.domain.is_power_space and
             all(isinstance(spc, TensorSpace) for spc in op.domain))):
        raise TypeError('operator domain {!r} is neither `TensorSpace` '
                        'nor `ProductSpace` with only equal `TensorSpace` '
                        'components'.format(op.domain))

    if not (isinstance(op.range, TensorSpace) or
            (isinstance(op.range, ProductSpace) and
             op.range.is_power_space and
             all(isinstance(spc, TensorSpace) for spc in op.range))):
        raise TypeError('operator range {!r} is neither `TensorSpace` '
                        'nor `ProductSpace` with only equal `TensorSpace` '
                        'components'.format(op.range))

    # Generate the matrix
    dtype = np.promote_types(op.domain.dtype, op.range.dtype)
    matrix = np.zeros(op.range.shape + op.domain.shape, dtype=dtype)
    tmp_ran = op.range.element()  # Store for reuse in loop
    tmp_dom = op.domain.zero()  # Store for reuse in loop

    for j in nd_iterator(op.domain.shape):
        tmp_dom[j] = 1.0

        op(tmp_dom, out=tmp_ran)
        matrix[(Ellipsis,) + j] = tmp_ran.asarray()

        tmp_dom[j] = 0.0

    return matrix