Decompose a 2-qubit unitary composed of two 1-qubit local gates.

    Uses the "Nearest Kronecker Product" algorithm. Will give erratic
    results if the gate is not the direct product of two 1-qubit gates.
979520
def kronecker_decomposition(gate: Gate) -> Circuit:
    
    # An alternative approach would be to take partial traces, but
    # this approach appears to be more robust.

    if gate.qubit_nb != 2:
        raise ValueError('Expected 2-qubit gate')

    U = asarray(gate.asoperator())
    rank = 2**gate.qubit_nb
    U /= np.linalg.det(U) ** (1/rank)

    R = np.stack([U[0:2, 0:2].reshape(4),
                  U[0:2, 2:4].reshape(4),
                  U[2:4, 0:2].reshape(4),
                  U[2:4, 2:4].reshape(4)])
    u, s, vh = np.linalg.svd(R)
    v = vh.transpose()
    A = (np.sqrt(s[0]) * u[:, 0]).reshape(2, 2)
    B = (np.sqrt(s[0]) * v[:, 0]).reshape(2, 2)

    q0, q1 = gate.qubits
    g0 = Gate(A, qubits=[q0])
    g1 = Gate(B, qubits=[q1])

    if not gates_close(gate, Circuit([g0, g1]).asgate()):
        raise ValueError("Gate cannot be decomposed into two 1-qubit gates")

    circ = Circuit()
    circ += zyz_decomposition(g0)
    circ += zyz_decomposition(g1)

    assert gates_close(gate, circ.asgate())  # Sanity check

    return circ
1243227
def do_gate_matrix(self, matrix: np.ndarray,
                       qubits: Sequence[int]) -> 'AbstractQuantumSimulator':
        
        # e.g. 2-qubit matrix is 4x4; turns into (2,2,2,2) tensor.
        tensor = np.reshape(matrix, (2,) * len(qubits) * 2)

        # Note to developers: you can use either einsum- or tensordot- based functions.
        # tensordot seems a little faster, but feel free to experiment.
        # self.wf = targeted_einsum(gate=gate_matrix, wf=self.wf, wf_target_inds=qubits)
        self.wf = targeted_tensordot(gate=tensor, wf=self.wf, wf_target_inds=qubits)
        return self
491081
private boolean _decompose() {
        double h[] = QH.data;

        for( int k = 0; k < N-2; k++ ) { // k = column
            u[k*2] = 0;
            u[k*2+1] = 0;
            double max = QrHelperFunctions_ZDRM.extractColumnAndMax(QH,k+1,N,k,u,0);

            if( max > 0 ) {
                // -------- set up the reflector Q_k

                double gamma = QrHelperFunctions_ZDRM.computeTauGammaAndDivide(k+1,N,u,max,tau);
                gammas[k] = gamma;

                // divide u by u_0
                double real_u_0 = u[(k+1)*2]   + tau.real;
                double imag_u_0 = u[(k+1)*2+1] + tau.imaginary;
                QrHelperFunctions_ZDRM.divideElements(k + 2, N, u, 0, real_u_0,imag_u_0 );

                // write the reflector into the lower left column of the matrix
                for (int i = k+2; i < N; i++) {
                    h[(i*N+k)*2]   = u[i*2];
                    h[(i*N+k)*2+1] = u[i*2+1];
                }

                u[(k+1)*2]   = 1;
                u[(k+1)*2+1] = 0;

                // ---------- multiply on the left by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultR(QH, u,0, gamma, k + 1, k + 1, N, b);

                // ---------- multiply on the right by Q_k
                QrHelperFunctions_ZDRM.rank1UpdateMultL(QH, u,0, gamma, 0, k + 1, N);

                // since the first element in the householder vector is known to be 1
                // store the full upper hessenberg
                h[((k+1)*N+k)*2]   = -tau.real*max;
                h[((k+1)*N+k)*2+1] = -tau.imaginary*max;

            } else {
                gammas[k] = 0;
            }

        }

        return true;
    }
876006
public void decompose( DMatrixRMaj A ) {
        init(A);

        for( int k = 1; k < N; k++ ) {
            similarTransform(k);
//            System.out.println("k=="+k);
//            QT.print();
        }
    }
446316
def partial_trace(self, qubits: Qubits) -> 'QubitVector':
        
        N = self.qubit_nb
        R = self.rank
        if R == 1:
            raise ValueError('Cannot take trace of vector')

        new_qubits: List[Qubit] = list(self.qubits)
        for q in qubits:
            new_qubits.remove(q)

        if not new_qubits:
            raise ValueError('Cannot remove all qubits with partial_trace.')

        indices = [self.qubits.index(qubit) for qubit in qubits]
        subscripts = list(EINSUM_SUBSCRIPTS)[0:N*R]
        for idx in indices:
            for r in range(1, R):
                subscripts[r * N + idx] = subscripts[idx]
        subscript_str = ''.join(subscripts)

        # Only numpy's einsum works with repeated subscripts
        tensor = self.asarray()
        tensor = np.einsum(subscript_str, tensor)

        return QubitVector(tensor, new_qubits)
272641
def zyz_decomposition(gate: Gate) -> Circuit:
    
    if gate.qubit_nb != 1:
        raise ValueError('Expected 1-qubit gate')

    q, = gate.qubits

    U = asarray(gate.asoperator())
    U /= np.linalg.det(U) ** (1/2)    # SU(2)

    if abs(U[0, 0]) > abs(U[1, 0]):
        theta1 = 2 * np.arccos(min(abs(U[0, 0]), 1))
    else:
        theta1 = 2 * np.arcsin(min(abs(U[1, 0]), 1))

    cos_halftheta1 = np.cos(theta1/2)
    if not np.isclose(cos_halftheta1, 0.0):
        phase = U[1, 1] / cos_halftheta1
        theta0_plus_theta2 = 2 * np.arctan2(np.imag(phase), np.real(phase))
    else:
        theta0_plus_theta2 = 0.0

    sin_halftheta1 = np.sin(theta1/2)
    if not np.isclose(sin_halftheta1, 0.0):
        phase = U[1, 0] / sin_halftheta1
        theta0_sub_theta2 = 2 * np.arctan2(np.imag(phase), np.real(phase))
    else:
        theta0_sub_theta2 = 0.0

    theta0 = (theta0_plus_theta2 + theta0_sub_theta2) / 2
    theta2 = (theta0_plus_theta2 - theta0_sub_theta2) / 2

    t0 = theta0/np.pi
    t1 = theta1/np.pi
    t2 = theta2/np.pi

    circ1 = Circuit()
    circ1 += TZ(t2, q)
    circ1 += TY(t1, q)
    circ1 += TZ(t0, q)

    return circ1
1008123
def do_gate(self, gate: Gate):
        
        unitary = lifted_gate(gate=gate, n_qubits=self.n_qubits)
        self.wf = unitary.dot(self.wf)
        return self
1302773
@Override
    public boolean decompose( ZMatrixRMaj A )
    {
        if( A.numRows != A.numCols )
            throw new IllegalArgumentException("A must be square.");
        if( A.numRows <= 0 )
            return false;

        QH = A;

        N = A.numCols;

        if( b.length < N*2 ) {
            b = new double[ N*2 ];
            gammas = new double[ N ];
            u = new double[ N*2 ];
        }
        return _decompose();
    }
24892
def _get_gate_tensor_and_qubits(gate: Gate):
    
    if len(gate.params) > 0:
        matrix = QUANTUM_GATES[gate.name](*gate.params)
    else:
        matrix = QUANTUM_GATES[gate.name]

    qubit_inds = [q.index for q in gate.qubits]

    # e.g. 2-qubit matrix is 4x4; turns into (2,2,2,2) tensor.
    tensor = np.reshape(matrix, (2,) * len(qubit_inds) * 2)

    return tensor, qubit_inds
321087
def join_states(*states: State) -> State:
    """"""
    vectors = [ket.vec for ket in states]
    vec = reduce(outer_product, vectors)
    return State(vec.tensor, vec.qubits)
78542
def _disjoint_qubits(op1: ops.Operation, op2: ops.Operation) -> bool:
    """"""
    return not set(op1.qubits) & set(op2.qubits)
377574
def quaternion_product(quat1, quat2):
    """"""
    return np.array([
        quat1[0]*quat2[0] - np.dot(quat1[1:], quat2[1:]),
        quat1[0]*quat2[1] + quat2[0]*quat1[1] + quat1[2]*quat2[3] - quat1[3]*quat2[2],
        quat1[0]*quat2[2] + quat2[0]*quat1[2] + quat1[3]*quat2[1] - quat1[1]*quat2[3],
        quat1[0]*quat2[3] + quat2[0]*quat1[3] + quat1[1]*quat2[2] - quat1[2]*quat2[1]
    ], float)
801583
def xor_nonlocal_decompose(qubits: Iterable[raw_types.Qid],
                           onto_qubit: raw_types.Qid
                          ) -> Iterable[raw_types.Operation]:
    """"""
    for qubit in qubits:
        if qubit != onto_qubit:
            yield common_gates.CNOT(qubit, onto_qubit)
670662
def factoring_qaoa(n_step, num, minimizer=None, sampler=None, verbose=True):
    
    def get_nbit(n):
        m = 1
        while 2**m < n:
            m += 1
        return m

    n1_bits = get_nbit(int(num**0.5)) - 1
    n2_bits = get_nbit(int(num**0.5))

    def mk_expr(offset, n):
        expr = pauli.Expr.from_number(1)
        for i in range(n):
            expr = expr + 2**(i + 1) * q(i + offset)
        return expr

    def bitseparator(bits):
        assert len(bits) == n1_bits + n2_bits
        p = 1
        m = 1
        for b in bits[:n1_bits]:
            if b:
                p += 2**m
            m += 1
        q = 1
        m = 1
        for b in bits[n1_bits:]:
            if b:
                q += 2**m
            m += 1
        return p, q

    hamiltonian = (num - mk_expr(0, n1_bits) * mk_expr(n1_bits, n2_bits))**2
    return vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, n_step), minimizer, sampler), bitseparator
734044
def pauli_product(*elements: Pauli) -> Pauli:
    """"""
    result_terms = []

    for terms in product(*elements):
        coeff = reduce(mul, [term[1] for term in terms])
        ops = (term[0] for term in terms)
        out = []
        key = itemgetter(0)
        for qubit, qops in groupby(heapq.merge(*ops, key=key), key=key):
            res = next(qops)[1]  # Operator: X Y Z
            for op in qops:
                pair = res + op[1]
                res, rescoeff = PAULI_PROD[pair]
                coeff *= rescoeff
            if res != 'I':
                out.append((qubit, res))

        p = Pauli(((tuple(out), coeff),))
        result_terms.append(p)

    return pauli_sum(*result_terms)
476951
def converted_gate_set(circuit: circuits.Circuit,
                       no_clifford_gates: bool = False,
                       atol: float = 1e-8,
                       ) -> circuits.Circuit:
    
    conv_circuit = circuits.Circuit(circuit)
    optimizers.ConvertToCzAndSingleGates().optimize_circuit(conv_circuit)
    optimizers.MergeSingleQubitGates().optimize_circuit(conv_circuit)
    ConvertToPauliStringPhasors(ignore_failures=True,
                                keep_clifford=not no_clifford_gates,
                                atol=atol,
                                ).optimize_circuit(conv_circuit)
    optimizers.DropEmptyMoments().optimize_circuit(conv_circuit)
    return conv_circuit
1048946
def decoherence_noise_with_asymmetric_ro(gates: Sequence[Gate], p00=0.975, p11=0.911):
    
    noise_model = _decoherence_noise_model(gates)
    aprobs = np.array([[p00, 1 - p00],
                       [1 - p11, p11]])
    aprobs = {q: aprobs for q in noise_model.assignment_probs.keys()}
    return NoiseModel(noise_model.gates, aprobs)