Python operators.qeye函数代码示例

def _jc_liouvillian(N):
    from qutip.tensor import tensor
    from qutip.operators import destroy, qeye
    from qutip.superoperator import liouvillian
    wc = 1.0  * 2 * np.pi  # cavity frequency
    wa = 1.0  * 2 * np.pi  # atom frequency
    g  = 0.05 * 2 * np.pi  # coupling strength
    kappa = 0.005          # cavity dissipation rate
    gamma = 0.05           # atom dissipation rate
    n_th_a = 1           # temperature in frequency units
    use_rwa = 0
    # operators
    a  = tensor(destroy(N), qeye(2))
    sm = tensor(qeye(N), destroy(2))
    # Hamiltonian
    if use_rwa:
        H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
    else:
        H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * (sm + sm.dag())
    c_op_list = []

    rate = kappa * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * a)

    rate = kappa * n_th_a
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * a.dag())

    rate = gamma
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sm)

    return liouvillian(H, c_op_list)

def test_sp_bandwidth():
    "Sparse: Bandwidth"
    # Bandwidth test 1
    A = create(25)+destroy(25)+qeye(25)
    band = sp_bandwidth(A.data)
    assert_equal(band[0], 3)
    assert_equal(band[1] == band[2] == 1, 1)
    # Bandwidth test 2
    A = np.array([[1, 0, 0, 0, 1, 0, 0, 0],
                  [0, 1, 1, 0, 0, 1, 0, 1],
                  [0, 1, 1, 0, 1, 0, 0, 0],
                  [0, 0, 0, 1, 0, 0, 1, 0],
                  [1, 0, 1, 0, 1, 0, 0, 0],
                  [0, 1, 0, 0, 0, 1, 0, 1],
                  [0, 0, 0, 1, 0, 0, 1, 0],
                  [0, 1, 0, 0, 0, 1, 0, 1]], dtype=np.int32)
    A = sp.csr_matrix(A)
    out1 = sp_bandwidth(A)
    assert_equal(out1[0], 13)
    assert_equal(out1[1] == out1[2] == 6, 1)
    # Bandwidth test 3
    perm = reverse_cuthill_mckee(A)
    B = sp_permute(A, perm, perm)
    out2 = sp_bandwidth(B)
    assert_equal(out2[0], 5)
    assert_equal(out2[1] == out2[2] == 2, 1)
    # Asymmetric bandwidth test
    A = destroy(25)+qeye(25)
    out1 = sp_bandwidth(A.data)
    assert_equal(out1[0], 2)
    assert_equal(out1[1], 0)
    assert_equal(out1[2], 1)

def _subsystem_apply_reference(state, channel, mask):
    if isket(state):
        state = ket2dm(state)

    if isoper(channel):
        full_oper = tensor([channel if mask[j]
                            else qeye(state.dims[0][j])
                            for j in range(len(state.dims[0]))])
        return full_oper * state * full_oper.dag()
    else:
        # Go to Choi, then Kraus
        # chan_mat = array(channel.data.todense())
        choi_matrix = super_to_choi(channel)
        vals, vecs = eig(choi_matrix.full())
        vecs = list(map(array, zip(*vecs)))
        kraus_list = [sqrt(vals[j]) * vec2mat(vecs[j])
                      for j in range(len(vals))]
        # Kraus operators to be padded with identities:
        k_qubit_kraus_list = product(kraus_list, repeat=sum(mask))
        rho_out = Qobj(inpt=zeros(state.shape), dims=state.dims)
        for operator_iter in k_qubit_kraus_list:
            operator_iter = iter(operator_iter)
            op_iter_list = [next(operator_iter).conj().T if mask[j]
                            else qeye(state.dims[0][j])
                            for j in range(len(state.dims[0]))]
            full_oper = tensor(list(map(Qobj, op_iter_list)))
            rho_out = rho_out + full_oper * state * full_oper.dag()
        return Qobj(rho_out)

def test_unitarity_known():
    """
    Metrics: Unitarity for known cases.
    """
    def case(q_oper, known_unitarity):
        assert_almost_equal(unitarity(q_oper), known_unitarity)

    yield case, to_super(sigmax()), 1.0
    yield case, sum(map(
        to_super, [qeye(2), sigmax(), sigmay(), sigmaz()]
    )) / 4, 0.0
    yield case, sum(map(
        to_super, [qeye(2), sigmax()]
    )) / 2, 1 / 3.0

def test_rand_unitary_haar_unitarity():
    """
    Random Qobjs: Tests that unitaries are actually unitary.
    """
    U = rand_unitary_haar(5)
    I = qeye(5)

    assert_(U * U.dag() == I)

 def test_stinespring_dims(self):
     """
     Stinespring: Check that dims of channels are preserved.
     """
     # FIXME: not the most general test, since this assumes a map
     #        from square matrices to square matrices on the same space.
     chan = super_tensor(to_super(sigmax()), to_super(qeye(3)))
     A, B = to_stinespring(chan)
     assert_equal(A.dims, [[2, 3, 1], [2, 3]])
     assert_equal(B.dims, [[2, 3, 1], [2, 3]])

 def _prop_identity(self, U, tol=1e-6):
     """
     Returns True if and only if U is proportional to the
     identity.
     """
     if U[0, 0] != 0:
         norm_U = U / U[0, 0]
         return (qeye(U.dims[0]) - norm_U).norm() <= tol
     else:
         return False

def _powers(op, N):
    """
    Generator that yields powers of an operator `op`,
    through to `N`.
    """
    acc = qeye(op.dims[0])
    yield acc

    for _ in range(N - 1):
        acc *= op
        yield acc

def _opto_liouvillian(N):
    from qutip.tensor import tensor
    from qutip.operators import destroy, qeye
    from qutip.superoperator import liouvillian
    Nc = 5                      # Number of cavity states
    Nm = N                     # Number of mech states
    kappa = 0.3                 # Cavity damping rate
    E = 0.1                     # Driving Amplitude         
    g0 = 2.4*kappa              # Coupling strength
    Qm = 1e4                    # Mech quality factor
    gamma = 1/Qm                # Mech damping rate
    n_th = 1                    # Mech bath temperature
    delta = -0.43               # Detuning
    a = tensor(destroy(Nc), qeye(Nm))
    b = tensor(qeye(Nc), destroy(Nm))
    num_b = b.dag()*b
    num_a = a.dag()*a
    H = -delta*(num_a)+num_b+g0*(b.dag()+b)*num_a+E*(a.dag()+a)
    cc = np.sqrt(kappa)*a
    cm = np.sqrt(gamma*(1.0 + n_th))*b
    cp = np.sqrt(gamma*n_th)*b.dag()
    c_ops = [cc,cm,cp]

    return liouvillian(H, c_ops)

def _spin_hamiltonian(N):
    from qutip.tensor import tensor
    from qutip.operators import qeye, sigmax, sigmay, sigmaz
    # array of spin energy splittings and coupling strengths. here we use
    # uniform parameters, but in general we don't have too
    h  = 1.0 * 2 * np.pi * np.ones(N) 
    Jz = 0.1 * 2 * np.pi * np.ones(N)
    Jx = 0.1 * 2 * np.pi * np.ones(N)
    Jy = 0.1 * 2 * np.pi * np.ones(N)
    # dephasing rate
    gamma = 0.01 * np.ones(N)

    si = qeye(2)
    sx = sigmax()
    sy = sigmay()
    sz = sigmaz()

    sx_list = []
    sy_list = []
    sz_list = []

    for n in range(N):
        op_list = []
        for m in range(N):
            op_list.append(si)

        op_list[n] = sx
        sx_list.append(tensor(op_list))

        op_list[n] = sy
        sy_list.append(tensor(op_list))

        op_list[n] = sz
        sz_list.append(tensor(op_list))

    # construct the hamiltonian
    H = 0

    # energy splitting terms
    for n in range(N):
        H += - 0.5 * h[n] * sz_list[n]

    # interaction terms
    for n in range(N-1):
        H += - 0.5 * Jx[n] * sx_list[n] * sx_list[n+1]
        H += - 0.5 * Jy[n] * sy_list[n] * sy_list[n+1]
        H += - 0.5 * Jz[n] * sz_list[n] * sz_list[n+1]
    return H

def test_QobjUnitaryOper():
    "Qobj unitarity"
    # Check some standard operators
    Sx = sigmax()
    Sy = sigmay()
    assert_unitarity(qeye(4), True, "qeye(4) should be unitary.")
    assert_unitarity(Sx, True, "sigmax() should be unitary.")
    assert_unitarity(Sy, True, "sigmax() should be unitary.")
    assert_unitarity(sigmam(), False, "sigmam() should NOT be unitary.")
    assert_unitarity(destroy(10), False, "destroy(10) should NOT be unitary.")
    # Check multiplcation of unitary is unitary
    assert_unitarity(Sx*Sy, True, "sigmax()*sigmay() should be unitary.")
    # Check some other operations clear unitarity
    assert_unitarity(Sx+Sy, False, "sigmax()+sigmay() should NOT be unitary.")
    assert_unitarity(4*Sx, False, "4*sigmax() should NOT be unitary.")
    assert_unitarity(Sx*4, False, "sigmax()*4 should NOT be unitary.")
    assert_unitarity(4+Sx, False, "4+sigmax() should NOT be unitary.")
    assert_unitarity(Sx+4, False, "sigmax()+4 should NOT be unitary.")

def test_fidelity_known_cases():
    """
    Metrics: Checks fidelity against known cases.
    """
    ket0 = basis(2, 0)
    ket1 = basis(2, 1)
    ketp = (ket0 + ket1).unit()
    # A state that almost overlaps with |+> should serve as
    # a nice test case, especially since we can analytically
    # calculate its overlap with |+>.
    ketpy = (ket0 + np.exp(1j * np.pi / 4) * ket1).unit()

    mms = qeye(2).unit()

    assert_almost_equal(fidelity(ket0, ketp), 1 / np.sqrt(2))
    assert_almost_equal(fidelity(ket0, ket1), 0)
    assert_almost_equal(fidelity(ket0, mms),  1 / np.sqrt(2))
    assert_almost_equal(fidelity(ketp, ketpy),
        np.sqrt(
            (1 / 8) + (1 / 2 + 1 / (2 * np.sqrt(2))) ** 2
        )
    )

def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """

    L = H if issuper(H) else liouvillian(H, c_ops)

    tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
    N = np.prod(L.dims[0][0])

    A = L.full()
    b = spre(b_op).full()
    a = spre(a_op).full()

    tr_vec = np.transpose(mat2vec(tr_mat.full()))

    rho_ss = steadystate(L)
    rho = np.transpose(mat2vec(rho_ss.full()))

    I = np.identity(N * N)
    P = np.kron(np.transpose(rho), tr_vec)
    Q = I - P

    spectrum = np.zeros(len(wlist))

    for idx, w in enumerate(wlist):
        if use_pinv:
            MMR = np.linalg.pinv(-1.0j * w * I + A)
        else:
            MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))

        s = np.dot(tr_vec,
                   np.dot(a, np.dot(MMR, np.dot(b, np.transpose(rho)))))
        spectrum[idx] = -2 * np.real(s[0, 0])

    return spectrum

    def test_chi_known(self):
        """
        Superoperator: Chi-matrix for known cases is correct.
        """
        def case(S, chi_expected, silent=True):
            chi_actual = to_chi(S)
            chiq = Qobj(chi_expected, dims=[[[2], [2]], [[2], [2]]], superrep='chi')
            if not silent:
                print(chi_actual)
                print(chi_expected)
            assert_almost_equal((chi_actual - chiq).norm('tr'), 0)

        yield case, sigmax(), [
            [0, 0, 0, 0],
            [0, 4, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, to_super(sigmax()), [
            [0, 0, 0, 0],
            [0, 4, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, qeye(2), [
            [4, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, (-1j * sigmax() * pi / 4).expm(), [
            [2, 2j, 0, 0],
            [-2j, 2, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]

def qudit_swap(dim):
    # We should likely generalize this and include it in qip.gates.
    W = qeye([dim, dim])
    return tensor_swap(W, (0, 1))