def _defgate(self, program, gate_name, gate_matrix):
        """Defines a gate named gate_name with matrix gate_matrix in program. In addition, updates
         self.defined_gates to track what has been defined.

        :param Program program: Pyquil Program to add the defgate and gate to.
        :param str gate_name: The name of the gate to add to program.
        :param numpy.ndarray gate_matrix: The array corresponding to the gate to define.
        :return: the modified Program.
        :retype: Program
        """
        new_program = pq.Program()
        new_program += program
        if gate_name not in self.defined_gates:
            new_program.defgate(gate_name, gate_matrix)
            self.defined_gates.add(gate_name)
        return new_program 

def state_tomography_programs(state_prep, qubits=None,
                              rotation_generator=tomography.default_rotations):
    """
    Yield tomographic sequences that prepare a state with Quil program `state_prep` and then append
    tomographic rotations on the specified `qubits`. If `qubits is None`, it assumes all qubits in
    the program should be tomographically rotated.

    :param Program state_prep: The program to prepare the state to be tomographed.
    :param list|NoneType qubits: A list of Qubits or Numbers, to perform the tomography on. If
    `None`, performs it on all in state_prep.
    :param generator rotation_generator: A generator that yields tomography rotations to perform.
    :return: Program for state tomography.
    :rtype: Program
    """
    if qubits is None:
        qubits = state_prep.get_qubits()
    for tomography_program in rotation_generator(*qubits):
        state_tomography_program = Program(Pragma("PRESERVE_BLOCK"))
        state_tomography_program.inst(state_prep)
        state_tomography_program.inst(tomography_program)
        state_tomography_program.inst(Pragma("END_PRESERVE_BLOCK"))
        yield state_tomography_program 

def do_state_tomography(preparation_program, nsamples, cxn, qubits=None, use_run=False):
    """
    Method to perform both a QPU and QVM state tomography, and use the latter as
    as reference to calculate the fidelity of the former.

    :param Program preparation_program: Program to execute.
    :param int nsamples: Number of samples to take for the program.
    :param QVMConnection|QPUConnection cxn: Connection on which to run the program.
    :param list qubits: List of qubits for the program.
    to use in the tomography analysis.
    :param bool use_run: If ``True``, use append measurements on all qubits and use ``cxn.run``
        instead of ``cxn.run_and_measure``.
    :return: The state tomogram.
    :rtype: StateTomography
    """
    return tomography._do_tomography(preparation_program, nsamples, cxn, qubits,
                                     tomography.MAX_QUBITS_STATE_TOMO,
                                     StateTomography, state_tomography_programs,
                                     DEFAULT_STATE_TOMO_SETTINGS, use_run=use_run) 

def basis_state_preps(*qubits):
    """
    Generate a sequence of programs that prepares the measurement
    basis states of some set of qubits in the order such that the qubit
    with highest index is iterated over the most quickly:
    E.g., for ``qubits=(0, 1)``, it returns the circuits::

        I_0 I_1
        I_0 X_1
        X_0 I_1
        X_0 X_1

    :param list qubits: Each qubit to include in the basis state preparation.
    :return: Yields programs for each basis state preparation.
    :rtype: Program
    """
    for prep in cartesian_product([I, X], repeat=len(qubits)):
        basis_prep = Program(Pragma("PRESERVE_BLOCK"))
        for gate, qubit in zip(prep, qubits):
            basis_prep.inst(gate(qubit))
        basis_prep.inst(Pragma("END_PRESERVE_BLOCK"))
        yield basis_prep 

def density(self, pyquil_program):
        """
        Run program and compute the density matrix

        Loads and checks program if all gates are within the stabilizer set.
        Then executes program and returns the final density matrix for the
        stabilizer

        :param Program pyquil_program: a pyquil Program containing only
                                       CNOT-H-S-MEASUREMENT operations
        :return:
        """
        self.load_program(pyquil_program)
        self.tableau = self._n_qubit_tableau(self.num_qubits)
        self.kernel()
        stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau())
        pauli_ops = reduce(lambda x, y: x * y, [0.5 * (sI(0) + b) for b in stabilizers])
        return tensor_up(pauli_ops, self.num_qubits) 

def get_compiled_prog(theta):
    return Program([
        RZgate(-pi/2, 0),
        RXgate(-pi/2, 0),
        RZgate(-pi/2, 1),
        RXgate( pi/2, 1),
        CZgate(1, 0),
        RZgate(-pi/2, 1),
        RXgate(-pi/2, 1),
        RZgate(theta, 1),
        RXgate( pi/2, 1),
        CZgate(1, 0),
        RXgate( pi/2, 0),
        RZgate( pi/2, 0),
        RZgate(-pi/2, 1),
        RXgate( pi/2, 1),
        RZgate(-pi/2, 1),
    ]) 

def test_exp_circuit(qvm):
    true_wf = np.array([0.54030231-0.84147098j,
                        0.00000000+0.j,
                        0.00000000+0.j,
                        0.00000000+0.j,
                        0.00000000+0.j,
                        0.00000000+0.j,
                        0.00000000+0.j,
                        0.00000000+0.j])

    create2kill1 = PauliTerm("X", 1, -0.25)*PauliTerm("Y", 2)
    create2kill1 += PauliTerm("Y", 1, 0.25)*PauliTerm("Y", 2)
    create2kill1 += PauliTerm("Y", 1, 0.25)*PauliTerm("X", 2)
    create2kill1 += PauliTerm("X", 1, 0.25)*PauliTerm("X", 2)
    create2kill1 += PauliTerm("I", 0, 1.0)
    prog = Program()
    for term in create2kill1.terms:
        single_exp_prog = exponentiate(term)
        prog += single_exp_prog

    wf, _ = qvm.wavefunction(prog)
    wf = np.reshape(wf.amplitudes, -1)
    assert np.allclose(wf.dot(np.conj(wf).T), true_wf.dot(np.conj(true_wf).T)) 

def test_errors(qvm):
    # NOP unsupported
    prog = Program(NOP)
    with pytest.raises(TypeError):
        qvm.wavefunction(prog)
    with pytest.raises(TypeError):
        qvm.run(prog)
    with pytest.raises(TypeError):
        qvm.run_and_measure(prog)

    # WAIT unsupported
    prog = Program(WAIT)
    with pytest.raises(TypeError):
        qvm.wavefunction(prog)
    with pytest.raises(TypeError):
        qvm.run(prog)
    with pytest.raises(TypeError):
        qvm.run_and_measure(prog) 

def tests_against_cloud(qvm):

    # simple program
    p = Program(H(0))
    cloud_results = [[0], [0], [0], [0], [0], [0], [0], [0], [0], [0]]
    local_results = qvm.run(p, classical_addresses=[0], trials=10)
    assert len(cloud_results) == len(local_results)
    
    cloud_wf = dc.HADAMARD_WF
    local_wf, _ = qvm.wavefunction(p)
    assert np.allclose(cloud_wf, local_wf.amplitudes)

    # complex program
    p = Program(dc.QFT_8_INSTRUCTIONS)
    cloud_wf = dc.QFT_8_WF_PROBS
    local_wf, _ = qvm.wavefunction(p)
    assert np.allclose(cloud_wf, local_wf.amplitudes) 

def test_random_gates(qvm_unitary):
    p = Program().inst([H(0), H(1), H(0)])
    test_unitary = qvm_unitary.unitary(p)
    actual_unitary = np.kron(gate_matrix['H'], np.eye(2 ** 1))
    assert np.allclose(test_unitary, actual_unitary)

    p = Program().inst([H(0), X(1), Y(2), Z(3)])
    test_unitary = qvm_unitary.unitary(p)
    actual_unitary = np.kron(gate_matrix['Z'],
                             np.kron(gate_matrix['Y'],
                                     np.kron(gate_matrix['X'],
                                             gate_matrix['H'])))
    assert np.allclose(test_unitary, actual_unitary)

    p = Program().inst([X(2), CNOT(2, 1), CNOT(1, 0)])
    test_unitary = qvm_unitary.unitary(p)
    # gates are multiplied in 'backwards' order
    actual_unitary = np.kron(np.eye(2 ** 1), gate_matrix['CNOT']).dot(
                     np.kron(gate_matrix['CNOT'], np.eye(2 ** 1))).dot(
                     np.kron(gate_matrix['X'], np.eye(2 ** 2)))
    assert np.allclose(test_unitary, actual_unitary) 

def test_qaoa_unitary(qvm_unitary):
    wf_true = [0.00167784 + 1.00210180e-05*1j, 0.50000000 - 4.99997185e-01*1j,
               0.50000000 - 4.99997185e-01*1j, 0.00167784 + 1.00210180e-05*1j]

    prog = Program()
    prog.inst([RY(np.pi/2, 0), RX(np.pi, 0),
               RY(np.pi/2, 1), RX(np.pi, 1),
               CNOT(0, 1), RX(-np.pi/2, 1), RY(4.71572463191, 1),
               RX(np.pi/2, 1), CNOT(0, 1),
               RX(-2*2.74973750579, 0), RX(-2*2.74973750579, 1)])

    test_unitary = qvm_unitary.unitary(prog)
    wf_test = np.zeros(4)
    wf_test[0] = 1.0
    wf_test = test_unitary.dot(wf_test)
    assert np.allclose(wf_test, wf_true) 

def test_tensor_gates_single_qubit():
    prog = Program().inst([Hgate(0)])
    test_unitary = tensor_gates(gate_matrix, {}, prog.instructions[0], 1).toarray()
    true_unitary = gate_matrix['H']
    assert np.allclose(test_unitary, true_unitary)

    prog = Program().inst([Hgate(0)])
    test_unitary = tensor_gates(gate_matrix, {}, prog.instructions[0], 5).toarray()
    true_unitary = np.kron(np.eye(2**4), gate_matrix['H'])
    assert np.allclose(test_unitary, true_unitary)

    prog = Program().inst([RXgate(0.2, 3)])
    test_unitary = tensor_gates(gate_matrix, {}, prog.instructions[0], 5).toarray()
    true_unitary = np.kron(np.eye(2**1), np.kron(gate_matrix['RX'](0.2),  np.eye(2**3)))
    assert np.allclose(test_unitary, true_unitary)

    prog = Program().inst([RXgate(0.5, 4)])
    test_unitary = tensor_gates(gate_matrix, {}, prog.instructions[0], 5).toarray()
    true_unitary = np.kron(np.eye(2**0), np.kron(gate_matrix['RX'](0.5),  np.eye(2**4)))
    assert np.allclose(test_unitary, true_unitary) 

def test_simulation_phase():
    """
    Test if the Phase gate is applied correctly to the tableau

    S|0> = |0>

    S|+> = |R>
    """
    prog = Program().inst(S(0))
    qvmstab = QVM_Stabilizer(num_qubits=1)
    qvmstab._apply_phase(prog.instructions[0])
    true_stab = np.array([[1, 1, 0],
                          [0, 1, 0]])
    assert np.allclose(true_stab, qvmstab.tableau)

    prog = Program().inst([H(0), S(0)])
    qvmstab = QVM_Stabilizer(num_qubits=1)
    qvmstab._apply_hadamard(prog.instructions[0])
    qvmstab._apply_phase(prog.instructions[1])
    true_stab = np.array([[0, 1, 0],
                          [1, 1, 0]])
    assert np.allclose(true_stab, qvmstab.tableau) 

def run(self, pyquil_program, classical_addresses=None, trials=1):
        """
        Run a pyQuil program multiple times, accumulating the values deposited
        in a list of classical addresses.

        This uses numpy's inverse sampling method to calculate bit string
        outcomes

        :param Program pyquil_program: A pyQuil program.
        :param list classical_addresses: A list of classical addresses.
                                         This is ignored but kept to have
                                         similar input as Forest QVM.
        :param int trials: Number of shots to collect.

        :return: A list of lists of bits. Each sublist corresponds to the
                 values in `classical_addresses`.
        :rtype: list
        """
        results = []
        for trial in range(trials):
            _ = self.density(pyquil_program)
            results.append(self.classical_memory[classical_addresses])
        return list(results) 

def run(self, pyquil_program, classical_addresses=None, trials=1):
        """
        Run a pyQuil program multiple times, accumulating the values deposited
        in a list of classical addresses.

        :param Program pyquil_program: A pyQuil program.
        :param list classical_addresses: A list of classical addresses.
        :param int trials: Number of shots to collect.

        :return: A list of lists of bits. Each sublist corresponds to the
                 values in `classical_addresses`.
        :rtype: list
        """
        results = []
        for trial in range(trials):
            _, classical_vals = self.wavefunction(pyquil_program,
                                                  classical_addresses=classical_addresses)
            results.append(classical_vals)

        return results 

def apply(self, operations, **kwargs):
        """Run the QVM"""
        # pylint: disable=attribute-defined-outside-init
        if self.parametric_compilation and "pyqvm" not in self.qc.name:
            self.apply_parametric_program(operations, **kwargs)
        else:
            super().apply(operations, **kwargs)

        prag = Program(Pragma("INITIAL_REWIRING", ['"PARTIAL"']))

        if self.active_reset:
            prag += RESET()

        self.prog = prag + self.prog

        qubits = sorted(self.wiring.values())
        ro = self.prog.declare("ro", "BIT", len(qubits))
        for i, q in enumerate(qubits):
            self.prog.inst(MEASURE(q, ro[i]))

        self.prog.wrap_in_numshots_loop(self.shots) 

def build(self):
        """Builds this controlled gate.

        :return: The controlled gate, defined by this object.
        :rtype: Program
        """
        self.defined_gates = set(STANDARD_GATE_NAMES)
        prog = self._recursive_builder(self.operation,
                                       self.gate_name,
                                       self.control_qubits,
                                       self.target_qubit)
        return prog 

def sample_bad_readout(program, num_samples, assignment_probs, cxn):
    """
    Generate `n` samples of measuring all outcomes of a Quil `program`
    assuming the assignment probabilities `assignment_probs` by simulating the
    wave function on a qvm QVMConnection `cxn`

    :param pyquil.quil.Program program: The program.
    :param int num_samples: The number of samples
    :param numpy.ndarray assignment_probs: A matrix of assignment probabilities
    :param QVMConnection cxn: the QVM connection.
    :return: The resulting sampled outcomes from assignment_probs applied to cxn, one dimensional.
    :rtype: numpy.ndarray
    """
    wf = cxn.wavefunction(program)
    return sample_outcomes(assignment_probs.dot(abs(wf.amplitudes.ravel())**2), num_samples) 

def process_tomography_programs(process, qubits=None,
                                pre_rotation_generator=tomography.default_rotations,
                                post_rotation_generator=tomography.default_rotations):
    """
    Generator that yields tomographic sequences that wrap a process encoded by a QUIL program `proc`
    in tomographic rotations on the specified `qubits`.

    If `qubits is None`, it assumes all qubits in the program should be
    tomographically rotated.

    :param Program process: A Quil program
    :param list|NoneType qubits: The specific qubits for which to generate the tomographic sequences
    :param pre_rotation_generator: A generator that yields tomographic pre-rotations to perform.
    :param post_rotation_generator: A generator that yields tomographic post-rotations to perform.
    :return: Program for process tomography.
    :rtype: Program
    """
    if qubits is None:
        qubits = process.get_qubits()
    for tomographic_pre_rotation in pre_rotation_generator(*qubits):
        for tomography_post_rotation in post_rotation_generator(*qubits):
            process_tomography_program = Program(Pragma("PRESERVE_BLOCK"))
            process_tomography_program.inst(tomographic_pre_rotation)
            process_tomography_program.inst(process)
            process_tomography_program.inst(tomography_post_rotation)
            process_tomography_program.inst(Pragma("END_PRESERVE_BLOCK"))

            yield process_tomography_program 

def test_1_qubit_control():
    prog = Program()
    qubit = Qubit(0)
    control_qubit = Qubit(1)
    prog += (ControlledProgramBuilder()
             .with_controls([control_qubit])
             .with_target(qubit)
             .with_operation(SIGMA_Z)
             .with_gate_name(SIGMA_Z_NAME).build())
    # This should be one "CZ" instruction, from control_qubit to qubit.
    assert len(prog) == 1
    instruction = prog.instructions[0]
    assert instruction.name == (ControlledProgramBuilder()
                                .format_gate_name("C", SIGMA_Z_NAME))
    assert instruction.qubits == [control_qubit, qubit] 

def test_2_qubit_control():
    """Test that ControlledProgramBuilder builds the program correctly all the way through."""
    prog = Program()
    qubit = Qubit(0)
    control_qubit_one = Qubit(1)
    control_qubit_two = Qubit(2)
    prog += (ControlledProgramBuilder()
             .with_controls([control_qubit_one, control_qubit_two])
             .with_target(qubit)
             .with_operation(SIGMA_Z)
             .with_gate_name(SIGMA_Z_NAME).build())
    # This should be one "CZ" instruction, from control_qubit to qubit.
    assert len(prog) == 5
    # Run tests
    double_control_test(prog.instructions, qubit, control_qubit_one, control_qubit_two) 

def test_sample_bad_readout():
    np.random.seed(234)
    assignment_probs = .3*np.random.rand(4, 4) + np.eye(4)
    assignment_probs /= assignment_probs.sum(axis=0)[np.newaxis, :]
    cxn = Mock()
    cxn.wavefunction.return_value.amplitudes = 0.5j * np.ones(4)
    with patch("grove.tomography.utils.sample_outcomes") as so:
        ut.sample_bad_readout(Program(X(0), X(1), X(0), X(1)), 10000, assignment_probs, cxn)
        assert np.allclose(so.call_args[0][0], 0.25 * assignment_probs.sum(axis=1)) 

def test_append_measure_register():
    q0 = QubitPlaceholder()
    p = Program(H(q0), RX(np.pi/2, 0))
    p = append_measure_register(p)
    assert str(p[-1]) == "MEASURE 0 ro[1]" 

def run(self, pyquil_program, classical_addresses=None, trials=1):
        """
        Run pyquil program and return the results

        Loads and checks program if all gates are within the stabilizer set.
        Then executes program.  If measurements are requested then the measured
        results are returned

        :param Program pyquil_program: a pyquil Program containing only
                                       CNOT-H-S-MEASUREMENT operations
        :param classical_addresses: classical addresses to return
        :param trials: number of times to repeat the execution of the program
        :return: list of lists of classical memory after each run
        """
        self.load_program(pyquil_program)

        if classical_addresses is None:
            classical_addresses = get_classical_addresses_from_program(pyquil_program)

        results = []
        for _ in range(trials):
            # set up stabilizers
            self.tableau = self._n_qubit_tableau(self.num_qubits)
            self.kernel()
            results.append(list(map(int, self.classical_memory[classical_addresses])))
            # results.append([int(b) for b in self.classical_memory[classical_addresses]])
            assert results[-1] == [int(b) for b in self.classical_memory[classical_addresses]]

            # reset qvm
            self.memory_reset()
            self.program_counter = 0

        return results 

def memory_reset(self):
        if self.program is None:
            raise TypeError("Program must be loaded to call reset")

        # setup quantum and classical memory
        q_max, c_max = self.identify_bits()
        if c_max <= 512:  # allocate at least 512 cbits (as floor)
            c_max = 512
        self.num_qubits = q_max
        self.classical_memory = np.zeros(c_max).astype(bool) 

def expectation(self, pyquil_program, operator_programs=[Program()]):
        """
        Calculate the expectation value given a state prepared.

        :param pyquil_program: (pyquil.Program) object containing only protoQuil
                               instructions.
        :param operator_programs: (optional, list) of PauliTerms. Default is
                                  Identity operator.

        :return: expectation value of the operators.
        :rtype: float
        """
        # TODO
        raise NotImplementedError() 

def test_generate_arbitrary_states(qvm):
    for k in list(dc.ARBITRARY_STATE_GEN_INSTRUCTIONS.keys()):

        v = np.asarray(dc.ARBITRARY_STATE_GEN_WF[k])
        norm = np.sqrt(np.sum(np.multiply(np.conj(v), v)))

        p = Program(dc.ARBITRARY_STATE_GEN_INSTRUCTIONS[k])
        wf, _ = qvm.wavefunction(p)

        # check actual part of wavefunction
        assert np.allclose(v.reshape(-1), wf.amplitudes[:len(v)] * norm)

        # check remaining zeros part of wavefunction
        assert np.allclose(np.zeros((wf.amplitudes.shape[0] - len(v), 1)),
                           wf.amplitudes[len(v):] * norm) 

def test_qaoa_density():
    wf_true = [0.00167784 + 1.00210180e-05*1j, 0.50000000 - 4.99997185e-01*1j,
               0.50000000 - 4.99997185e-01*1j, 0.00167784 + 1.00210180e-05*1j]
    wf_true = np.reshape(np.array(wf_true), (4, 1))
    rho_true = np.dot(wf_true, np.conj(wf_true).T)
    prog = Program()
    prog.inst([RYgate(np.pi/2, 0), RXgate(np.pi, 0),
               RYgate(np.pi/2, 1), RXgate(np.pi, 1),
               CNOTgate(0, 1), RXgate(-np.pi/2, 1), RYgate(4.71572463191, 1),
               RXgate(np.pi/2, 1), CNOTgate(0, 1),
               RXgate(-2*2.74973750579, 0), RXgate(-2*2.74973750579, 1)])

    qvm = QVMConnection(type_trans='density')
    wf_rho = qvm.density(prog)
    assert np.isclose(rho_true, wf_rho.todense()).all() 

def test_kraus_through_qvm_t2():
    noise = NoiseModel(T1=INFINITY, ro_fidelity=1.0)
    qvm = QVMConnection(type_trans='density', noise_model=noise)
    rho = random_1q_density()
    identity_program = Program().inst(Igate(0))
    qvm._density = sps.csc_matrix(rho)
    qvm.num_qubits = 1
    qvm.load_program(identity_program)
    qvm.kernel()
    p = 1 - np.exp(-noise.gate_time_1q/noise.T2)
    final_density = np.array([[rho[0, 0], (1 - p) * rho[0, 1]],
                              [(1 - p) * rho[1, 0],  rho[1, 1]]])
    assert np.allclose(final_density, qvm._density.todense()) 

def test_kraus_through_qvm_t1t2():
    noise = NoiseModel(ro_fidelity=1.0)
    qvm = QVMConnection(type_trans='density', noise_model=noise)
    rho = random_1q_density()
    identity_program = Program().inst(Igate(0))
    qvm._density = sps.csc_matrix(rho)
    qvm.num_qubits = 1
    qvm.load_program(identity_program)
    qvm.kernel()
    p = 1 - np.exp(-noise.gate_time_1q/noise.T2)
    p1 = 1 - np.exp(-noise.gate_time_1q/noise.T1)
    final_density = np.array([[rho[0, 0] + p1 * rho[1, 1], np.sqrt(1 - p1) * (1 - p) * rho[0, 1]],
                              [np.sqrt(1 - p1) * (1 - p) * rho[1, 0],  (1 - p1) * rho[1, 1]]])
    assert np.allclose(final_density, qvm._density.todense())