How to use the qutip.qeye function in qutip

"You may pass a non-trivial `full_space`")
    if mapping is not None:
        if expr in mapping:
            ret = mapping[expr]
            if isinstance(ret, qutip.Qobj):
                return ret
            else:
                assert callable(ret)
                return ret(expr)
    if expr is IdentityOperator:
        local_spaces = full_space.local_factors
        if len(local_spaces) == 0:
            raise ValueError("full_space %s does not have local factors"
                             % full_space)
        else:
            return qutip.tensor(*[qutip.qeye(s.dimension)
                                  for s in local_spaces])
    elif expr is ZeroOperator:
        return qutip.tensor(
            *[qutip.Qobj(csr_matrix((s.dimension, s.dimension)))
              for s in full_space.local_factors]
        )
    elif isinstance(expr, LocalOperator):
        return _convert_local_operator_to_qutip(expr, full_space, mapping)
    elif (isinstance(expr, Operator) and isinstance(expr, Operation)):
        return _convert_operator_operation_to_qutip(expr, full_space, mapping)
    elif isinstance(expr, OperatorTrace):
        raise NotImplementedError('Cannot convert OperatorTrace to '
                                  'qutip')
        # actually, this is perfectly doable in principle, but requires a bit
        # of work
    elif isinstance(expr, State):

ss_conf = stats.sections.get('config')
            if ss_conf is None:
                ss_conf = stats.add_section('config')

        c, nu = self._calc_matsubara_params()

        if renorm:
            norm_plus, norm_minus = self._calc_renorm_factors()
            if stats:
                stats.add_message('options', 'renormalisation', ss_conf)
        # Dimensions et by system
        N_temp = 1
        for i in H_sys.dims[0]:
            N_temp *= i
        sup_dim = N_temp**2
        unit_sys = qeye(N_temp)


        # Use shorthands (mainly as in referenced PRL)
        lam0 = self.coup_strength
        gam = self.cut_freq
        N_c = self.N_cut
        N_m = self.N_exp
        Q = coup_op # Q as shorthand for coupling operator
        beta = 1.0/(self.boltzmann*self.temperature)

        # Ntot is the total number of ancillary elements in the hierarchy
        # Ntot = factorial(N_c + N_m) / (factorial(N_c)*factorial(N_m))
        # Turns out to be the same as nstates from state_number_enumerate
        N_he, he2idx, idx2he = enr_state_dictionaries([N_c + 1]*N_m , N_c)

        unit_helems = fast_identity(N_he)

def _convert_ket_to_qutip(expr, full_space, mapping):
    n = full_space.dimension
    if full_space != expr.space:
        all_spaces = full_space.local_factors
        own_space_index = all_spaces.index(expr.space)
        factors = (
            [qutip.qeye(s.dimension) for s in all_spaces[:own_space_index]] +
            [convert_to_qutip(expr, expr.space, mapping=mapping), ] +
            [qutip.qeye(s.dimension) for s in all_spaces[own_space_index + 1:]]
        )
        return qutip.tensor(*factors)
    if isinstance(expr, BasisKet):
        return qutip.basis(n, expr.index)
    elif isinstance(expr, CoherentStateKet):
        return qutip.coherent(n, complex(expr.operands[1]))
    elif isinstance(expr, KetPlus):
        return sum((convert_to_qutip(op, full_space, mapping=mapping)
                    for op in expr.operands), 0)
    elif isinstance(expr, TensorKet):
        if any(len(op.space) > 1 for op in expr.operands):
            se = expr.expand()
            if se == expr:
                raise ValueError("Cannot represent as QuTiP "
                                 "object: {!s}".format(expr))
            return convert_to_qutip(se, full_space, mapping=mapping)

.format(expr))
            return convert_to_qutip(se, full_space, mapping=mapping)
        all_spaces = full_space.local_factors
        by_space = []
        ck = 0
        for ls in all_spaces:
            # group factors by associated local space
            ls_ops = [convert_to_qutip(o, o.space, mapping=mapping)
                      for o in expr.operands if o.space == ls]
            if len(ls_ops):
                # compute factor associated with local space
                by_space.append(reduce(lambda a, b: a * b, ls_ops, 1))
                ck += len(ls_ops)
            else:
                # if trivial action, take identity matrix
                by_space.append(qutip.qeye(ls.dimension))
        assert ck == len(expr.operands)
        # combine local factors in tensor product
        return qutip.tensor(*by_space)
    elif isinstance(expr, Adjoint):
        return convert_to_qutip(qutip.dag(expr.operands[0]), full_space,
                                mapping=mapping)
    elif isinstance(expr, ScalarTimesOperator):
        try:
            coeff = complex(expr.coeff)
        except TypeError:
            raise TypeError("Scalar coefficient '%s' is not numerical" %
                            expr.coeff)
        return coeff * convert_to_qutip(expr.term, full_space=full_space,
                                        mapping=mapping)
    elif isinstance(expr, PseudoInverse):
        mo = convert_to_qutip(expr.operand, full_space=full_space,

if options is None:
        options = Options()

    dot_energy, dot_state = Hsys.eigenstates(sparse=sparse)
    deltaE = dot_energy[1] - dot_energy[0]
    if (w_th < deltaE/2):
        warnings.warn("Given w_th might not provide accurate results")
    gamma = deltaE / (2 * np.pi * wc)
    wa = 2 * np.pi * gamma * wc  # reaction coordinate frequency
    g = np.sqrt(np.pi * wa * alpha / 2.0)  # reaction coordinate coupling
    nb = (1 / (np.exp(wa/w_th) - 1))

    # Reaction coordinate hamiltonian/operators

    dimensions = dims(Q)
    a = tensor(destroy(N), qeye(dimensions[1]))
    unit = tensor(qeye(N), qeye(dimensions[1]))
    Nmax = N * dimensions[1][0]
    Q_exp = tensor(qeye(N), Q)
    Hsys_exp = tensor(qeye(N), Hsys)
    e_ops_exp = [tensor(qeye(N), kk) for kk in e_ops]

    na = a.dag() * a
    xa = a.dag() + a

    # decoupled Hamiltonian
    H0 = wa * a.dag() * a + Hsys_exp
    # interaction
    H1 = (g * (a.dag() + a) * Q_exp)
    H = H0 + H1
    L = 0
    PsipreEta = 0

def to_qutip(self, full_space):
        return qutip.tensor(*[qutip.qeye(HilbertSpace((s,)).dimension) for s in full_space])

def _generator(k, H, L1, L2, S=None, c_ops_markov=None):
    """
    Create a Liouvillian for a cascaded chain of k system copies
    """
    id = qt.qeye(H.dims[0][0])
    Id = qt.sprepost(id, id)
    if S is None:
        S = np.identity(len(L1))
    # create Lindbladian
    L = qt.Qobj()
    E0 = Id
    # first system
    L += qt.liouvillian(None, [_localop(c, 1, k) for c in L2])
    for l in range(1, k):
        # Identiy superoperator
        E0 = qt.composite(E0, Id)
        # Bare Hamiltonian
        Hl = _localop(H, l, k)
        L += qt.liouvillian(Hl, [])
        # Markovian Decay channels
        if c_ops_markov is not None:

else:
                inter_hop = Qobj(inter_hop, dims=dim_ih, type='oper')
            self._H_inter_list = [inter_hop]
            self._H_inter = inter_hop
            is_real = is_real and np.isreal(inter_hop).all()

        elif inter_hop is None:      # inter_hop is the default None)
            # So, we set self._H_inter_list from cell_num_site and
            # cell_site_dof
            if self._length_of_unit_cell == 1:
                inter_hop = Qobj([[-1]], type='oper')
            else:
                bNm = basis(cell_num_site, cell_num_site-1)
                bN0 = basis(cell_num_site, 0)
                siteT = -bNm * bN0.dag()
                inter_hop = tensor(Qobj(siteT), qeye(self.cell_site_dof))
            dih = inter_hop.dims[0]
            if all(x == 1 for x in dih):
                dih = [1]
            else:
                while 1 in dih:
                    dih.remove(1)
            self._H_inter_list = [Qobj(inter_hop, dims=[dih, dih],
                                       type='oper')]
            self._H_inter = Qobj(inter_hop, dims=[dih, dih], type='oper')
        else:
            raise Exception("inter_hop is required to be a Qobj or a \
                            list of Qobjs.")

        self.positions_of_sites = [(i/self.cell_num_site) for i in
                                   range(self.cell_num_site)]
        self._inter_vec_list = [[1] for i in range(len(self._H_inter_list))]

def _localop(op, l, k):
    """
    Create a local operator on the l'th system by tensoring
    with identity operators on all the other k-1 systems
    """
    if l < 1 or l > k:
        raise IndexError('index l out of range')
    h = op
    I = qt.qeye(op.shape[0])
    I.dims = op.dims
    for i in range(1, l):
        h = qt.tensor(I, h)
    for i in range(l+1, k+1):
        h = qt.tensor(h, I)
    return h

def _convert_local_operator_to_qutip(expr, full_space, mapping):
    """Convert a LocalOperator instance to qutip"""
    n = full_space.dimension
    if full_space != expr.space:
        all_spaces = full_space.local_factors
        own_space_index = all_spaces.index(expr.space)
        return qutip.tensor(
            *([qutip.qeye(s.dimension)
               for s in all_spaces[:own_space_index]] +
              [convert_to_qutip(expr, expr.space, mapping=mapping), ] +
              [qutip.qeye(s.dimension)
               for s in all_spaces[own_space_index + 1:]])
        )
    if isinstance(expr, Create):
        return qutip.create(n)
    elif isinstance(expr, Jz):
        return qutip.jmat((expr.space.dimension-1)/2., "z")
    elif isinstance(expr, Jplus):
        return qutip.jmat((expr.space.dimension-1)/2., "+")
    elif isinstance(expr, Jminus):
        return qutip.jmat((expr.space.dimension-1)/2., "-")
    elif isinstance(expr, Destroy):
        return qutip.destroy(n)
    elif isinstance(expr, Phase):
        arg = complex(expr.operands[1]) * arange(n)
        d = np_cos(arg) + 1j * np_sin(arg)
        return qutip.Qobj(np_diag(d))