Repository: rigetticomputing/pyquil
Pull request #354 (2018-03-07 22:44:29)

Revamp qubit placeholder handling
On the basis that automated assignment of QubitPlaceholder to integer was never actually used by anyone, I've revamped the way that works.

Now you can use QubitPlaceholder's somewhat effectively:
 - Programs must be defined with all placeholder (or all integers)
 - You can manually control the mapping of placeholders to integers which is useful for running on real devices. In particular: it's useful for potentially running the same program on different devices with different qubit indices
 - PauliSums now work with QubitPlaceholders


This PR makes a lot of things illegal. Rather than trying to do very general things badly, you can now do very specific things that work in a sensible way and I've left many places in the code for making dealing with placeholders more general

Still to do: 
 - tests
 - ~~think about the workflow where you might want to call synthesize multiple times with different mappings~~ done

mpharrigan: This video was posted yesterday https://www.youtube.com/watch?v=nc2q4OOK6K8 and is both fascinating and relevant :) 
willzeng: I haven't looked at this in full detail yet so I apologize if I'm missing something basic. It seems to me that it is a rare use case that you will actually only want to assign one `QubitPlaceholder` at a time. In the current implementation you could of course do `qubit_register = [QubitPlaceholder for qq in range(10)]` or something. However, it seems like a register allocation is actually the conceptual primitive where a single qubit is the special case.

Would something like this be a bit cleaner:
```
    p = Program()
    qreg = QubitRegister()
    p.inst(H(qreg[0]))
    p.inst(CNOT(qreg[1], qreg[2)) 
    p.inst(H(qreg[3]))
    p.inst(X(qreg[3)) # this X is on the same one as the preceding H
    p = address_qubits(p, [0,3,4,0,2])
```
It means the user doesn't have to name every qubit with a variable name (saves like 5 lines from the testing example that I took this from). You could of course still do this by using a variable for the index, e.g. `qxxx = 3` and `p += H(qreg[qxxx])`.

Under the hood this could be implemented by overriding the `__getitem__` method on the `QubitRegister` to either produce a new `QubitPlaceholder` or pass something that points to an existing one.

Thoughts @mpharrigan @stevenheidel ?
stevenheidel: I like the idea of a QubitRegister, basically every other quantum package has something like this so unless we know something they don't this seems like a good idea.
mpharrigan: That sounds good, but I think `QubitRegister()` should just return a list of `QubitPlaceholder` objects instead of introducing a bespoke class with exactly the same semantics as a list
mpharrigan: In particular, I don't think we should be creating qubit placeholders on the fly implicitly unless there's actually user demand for that kind of trickery. For now, we should require explicit allocation of a number of qubits (and if you need more, you allocate some more)
mpharrigan: This is ready for review @stevenheidel et al. You can now request a register of qubits

```python
q = QubitPlaceholder.register(8)
H = sum(sZ(q[i]) * sZ(q[j]) for i, j in itertools.product(range(8), repeat=2))
```
jotterbach: Pinged you already offline, but want to bring this to the top of your inbox again @stevenheidel 
stevenheidel: @mpharrigan - overall this looks ready to be merged. Sorry I missed your ping 2 weeks ago, must have archived the email. If you want to push this through let's just meet in person to address the last few things and get it in. (plus do a 1.9 release)

Patch:
diff --git a/docs/Makefile b/docs/Makefile
index b549f72..69fe55e 100644
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -4,7 +4,6 @@
 # You can set these variables from the command line.
 SPHINXOPTS    =
 SPHINXBUILD   = sphinx-build
-SPHINXPROJ    = redwood
 SOURCEDIR     = source
 BUILDDIR      = build
 
diff --git a/docs/README.md b/docs/README.md
index 481c041..d3116c0 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -10,4 +10,8 @@ To build the docs, run:
 
 `make docs`
 
-from the `pyquil` root folder.
+from the `pyquil` root folder. Or run:
+
+`make html`
+
+from this folder.
diff --git a/docs/source/changes.rst b/docs/source/changes.rst
index 444e7da..b5103ef 100644
--- a/docs/source/changes.rst
+++ b/docs/source/changes.rst
@@ -13,3 +13,23 @@ v1.9
    backwards compatibility. However, this function should not be used for comparing term-type like
    it has been used previously. Use :py:func:`PauliTerm.operations_as_set()` instead. In the future,
    ``sort_ops`` will default to False and will eventually be removed.
+ - :py:func:`Program.alloc()` has been deprecated. Please instantiate :py:class:`QubitPlaceholder()`
+   directly or request a "register" (list) of ``n`` placeholders by using the class constructor
+   :py:func:`QubitPlaceholder.register(n)`.
+ - Programs must contain either (1) all instantiated qubits with integer indexes or (2) all
+   placeholder qubits of type :py:class:`QubitPlaceholder`. We have found that most users use
+   (1) but (2) will become useful with larger and more diverse devices.
+ - Programs that contain qubit placeholders must be **explicitly addressed** prior to execution.
+   Previously, qubits would be assigned "under the hood" to integers 0...N. Now, you must use
+   :py:func:`address_qubits` which returns a new program with all qubits indexed depending
+   on the ``qubit_mapping`` argument. The original program is unaffected and can be "readdressed"
+   multiple times.
+ - In light of the above, ``shift_quantum_gates`` has been removed. Users who relied on this
+   functionality should use :py:class:`QubitPlaceholder` and :py:func:`address_qubits` to
+   achieve the same result. Users should also double-check data resulting from use of this function
+   as there were several edge cases which would cause the shift to be applied incorrectly resulting
+   in badly-addressed qubits.
+ - :py:class:`PauliTerm` can now accept :py:class:`QubitPlaceholder` in addition to integers.
+ - :py:class:`QubitPlaceholder` is no longer a subclass of :py:class:`Qubit`.
+   :py:class:`LabelPlaceholder` is no longer a subclass of :py:class:`Label`.
+ - :py:class:`QuilAtom` subclasses' hash functions have changed.
diff --git a/pyquil/api/qpu.py b/pyquil/api/qpu.py
index f88820d..0bcbbfe 100644
--- a/pyquil/api/qpu.py
+++ b/pyquil/api/qpu.py
@@ -14,9 +14,11 @@
 #    limitations under the License.
 ##############################################################################
 import warnings
+import time
 
 from six import integer_types
 
+from pyquil.api import errors
 from pyquil.api.job import Job
 from pyquil.device import Device
 from pyquil.gates import MEASURE
@@ -190,7 +192,14 @@ with the former, the device.
             classical_addresses = get_classical_addresses_from_program(quil_program)
 
         payload = self._run_payload(quil_program, classical_addresses, trials, needs_compilation=needs_compilation, isa=isa)
-        response = post_json(self.session, self.async_endpoint + "/job", self._wrap_program(payload))
+        response = None
+        while response is None:
+            try:
+                response = post_json(self.session, self.async_endpoint + "/job", self._wrap_program(payload))
+            except errors.DeviceRetuningError:
+                print("QPU is retuning. Will try to reconnect in 10 seconds...")
+                time.sleep(10)
+
         return get_job_id(response)
 
     def _run_payload(self, quil_program, classical_addresses, trials, needs_compilation, isa):
diff --git a/pyquil/parameters.py b/pyquil/parameters.py
index e660024..679258f 100644
--- a/pyquil/parameters.py
+++ b/pyquil/parameters.py
@@ -31,9 +31,11 @@ def format_parameter(element):
             out += repr(r)
 
         if i == 1:
-            out += 'i'
+            assert np.isclose(r, 0, atol=1e-14)
+            out = 'i'
         elif i == -1:
-            out += '-i'
+            assert np.isclose(r, 0, atol=1e-14)
+            out = '-i'
         elif i < 0:
             out += repr(i) + 'i'
         else:
@@ -136,6 +138,7 @@ class Parameter(QuilAtom, Expression):
     """
     Parameters in Quil are represented as a label like '%x' for the parameter named 'x'.
     """
+
     def __init__(self, name):
         self.name = name
 
@@ -145,6 +148,15 @@ class Parameter(QuilAtom, Expression):
     def _substitute(self, d):
         return d.get(self, self)
 
+    def __str__(self):
+        return '%' + self.name
+
+    def __hash__(self):
+        return hash(self.name)
+
+    def __eq__(self, other):
+        return isinstance(other, Parameter) and other.name == self.name
+
 
 class Function(Expression):
     """
diff --git a/pyquil/paulis.py b/pyquil/paulis.py
index 17e61c2..2fc0143 100644
--- a/pyquil/paulis.py
+++ b/pyquil/paulis.py
@@ -21,6 +21,9 @@ from __future__ import division
 from itertools import product
 import numpy as np
 import copy
+
+from pyquil.quilatom import QubitPlaceholder
+
 from .quil import Program
 from .gates import H, RZ, RX, CNOT, X, PHASE, STANDARD_GATES
 from numbers import Number
@@ -57,6 +60,11 @@ can't use np.isclose() for hashing terms though.
 """
 
 
+def _valid_qubit(index):
+    return ((isinstance(index, integer_types) and index >= 0) or
+            isinstance(index, QubitPlaceholder))
+
+
 class PauliTerm(object):
     """A term is a product of Pauli operators operating on different qubits.
     """
@@ -70,7 +78,7 @@ class PauliTerm(object):
         :param float coefficient: The coefficient multiplying the operator, e.g. 1.5 * Z_1
         """
         assert op in PAULI_OPS
-        assert isinstance(index, integer_types) and index >= 0
+        assert _valid_qubit(index)
 
         self._ops = OrderedDict()
         if op != "I":
@@ -308,7 +316,7 @@ class PauliTerm(object):
         assert all([op[0] in PAULI_OPS for op in terms_list])
 
         indices = [op[1] for op in terms_list]
-        assert all([isinstance(index, integer_types) and index >= 0 for index in indices])
+        assert all(_valid_qubit(index) for index in indices)
 
         # this is because from_list doesn't call simplify in order to be more efficient.
         if len(set(indices)) != len(indices):
diff --git a/pyquil/quil.py b/pyquil/quil.py
index a9c055e..ff26d21 100644
--- a/pyquil/quil.py
+++ b/pyquil/quil.py
@@ -16,10 +16,11 @@
 """
 Module for creating and defining Quil programs.
 """
+import itertools
+import types
 import warnings
-from itertools import count
+from collections import OrderedDict
 from math import pi
-import types
 
 import numpy as np
 from six import string_types
@@ -28,21 +29,23 @@ from pyquil._parser.PyQuilListener import run_parser
 from pyquil.noise import _check_kraus_ops, _create_kraus_pragmas
 from pyquil.parameters import format_parameter
 from pyquil.quilatom import LabelPlaceholder, QubitPlaceholder, unpack_qubit
-from .gates import MEASURE, STANDARD_GATES, H
-from .quilbase import (DefGate, Gate, Measurement, Pragma, AbstractInstruction, Qubit,
-                       Jump, Label, JumpConditional, JumpTarget, JumpUnless, JumpWhen, Addr)
+from pyquil.gates import MEASURE, STANDARD_GATES, H
+from pyquil.quilbase import (DefGate, Gate, Measurement, Pragma, AbstractInstruction, Qubit,
+                             Jump, Label, JumpConditional, JumpTarget, JumpUnless, JumpWhen, Addr)
 
 
 class Program(object):
     def __init__(self, *instructions):
         self._defined_gates = []
-        # Implementation note: the key difference between the private _instructions and the public instructions
-        # property below is that the private _instructions list may contain placeholder values
+        # Implementation note: the key difference between the private _instructions and
+        # the public instructions property below is that the private _instructions list
+        # may contain placeholder labels.
         self._instructions = []
 
-        # Performance optimization: as stated above _instructions may contain placeholder values so the program must
-        # first be synthesized. _synthesized_instructions is simply a cache on the result of the _synthesize() method.
-        # It is marked as None whenever new instructions are added.
+        # Performance optimization: as stated above _instructions may contain placeholder
+        # labels so the program must first be have its labels instantiated.
+        # _synthesized_instructions is simply a cache on the result of the _synthesize()
+        # method.  It is marked as None whenever new instructions are added.
         self._synthesized_instructions = None
 
         self.inst(*instructions)
@@ -60,7 +63,7 @@ class Program(object):
         Fill in any placeholders and return a list of quil AbstractInstructions.
         """
         if self._synthesized_instructions is None:
-            self._synthesized_instructions = self._synthesize()
+            self._synthesize()
 
         return self._synthesized_instructions
 
@@ -125,11 +128,11 @@ class Program(object):
                     self.inst(instr)
 
             # Implementation note: these two base cases are the only ones which modify the program
-
             elif isinstance(instruction, DefGate):
                 defined_gate_names = [gate.name for gate in self._defined_gates]
                 if instruction.name in defined_gate_names:
-                    warnings.warn("Gate {} has already been defined in this program".format(instruction.name))
+                    warnings.warn("Gate {} has already been defined in this program"
+                                  .format(instruction.name))
 
                 self._defined_gates.append(instruction)
             elif isinstance(instruction, AbstractInstruction):
@@ -353,24 +356,33 @@ class Program(object):
         :return: A qubit.
         :rtype: Qubit
         """
+        warnings.warn("`alloc` is deprecated and will be removed in a future version of pyQuil. "
+                      "Please create a `QubitPlaceholder` directly", DeprecationWarning)
         return QubitPlaceholder()
 
-    def out(self):
+    def _out(self, allow_placeholders):
         """
         Converts the Quil program to a readable string.
 
-        :return: String form of a program
-        :rtype: string
+        :param allow_placeholders: Whether to complain if the program contains placeholders.
         """
-        s = ""
-        for dg in self._defined_gates:
-            s += dg.out()
-            s += "\n"
-        for instr in self.instructions:
-            s += instr.out() + "\n"
-        return s
+        return '\n'.join(itertools.chain(
+            (dg.out() for dg in self._defined_gates),
+            (instr.out(allow_placeholders=allow_placeholders) for instr in self.instructions),
+            [''],
+        ))
+
+    def out(self):
+        """
+        Serializes the Quil program to a string suitable for submitting to the QVM or QPU.
+        """
+        return '\n'.join(itertools.chain(
+            (dg.out() for dg in self._defined_gates),
+            (instr.out() for instr in self.instructions),
+            [''],
+        ))
 
-    def get_qubits(self):
+    def get_qubits(self, indices=True):
         """
         Returns all of the qubit indices used in this program, including gate applications and
         allocated qubits. e.g.
@@ -384,15 +396,15 @@ class Program(object):
             >>> len(p.get_qubits())
             2
 
+        :param indices: Return qubit indices as integers intead of the
+            wrapping :py:class:`Qubit` object
         :return: A set of all the qubit indices used in this program
         :rtype: set
         """
         qubits = set()
         for instr in self.instructions:
-            if isinstance(instr, Gate):
-                qubits |= {q.index for q in instr.qubits}
-            elif isinstance(instr, Measurement):
-                qubits.add(instr.qubit.index)
+            if isinstance(instr, (Gate, Measurement)):
+                qubits |= instr.get_qubits(indices=indices)
         return qubits
 
     def is_protoquil(self):
@@ -460,94 +472,19 @@ class Program(object):
 
     def _synthesize(self):
         """
-        Takes a program which may contain placeholders and assigns them all defined values.
+        Assigns all placeholder labels to actual values.
 
-        For qubit placeholders:
-        1. We look through the program to find all the known indexes of qubits and add them to a set
-        2. We create a mapping from undefined qubits to their newly assigned index
-        3. For every qubit placeholder in the program, if it's not already been assigned then look through the set of
-            known indexes and find the lowest available one
+        Changed in 1.9: Either all qubits must be defined or all undefined. If qubits are
+        undefined, this method will not help you. You must explicitly call `address_qubits`
+        which will return a new Program.
 
-        For label placeholders:
-        1. Start a counter at 1
-        2. For every label placeholder in the program, replace it with a defined label using the counter and increment
-            the counter
+        Changed in 1.9: This function now returns ``self`` and updates
+        ``self._synthesized_instructions``.
 
-        :return: List of AbstractInstructions with all placeholders removed
+        :return: This object with the ``_synthesized_instructions`` member set.
         """
-        used_indexes = set()
-        for instr in self._instructions:
-            if isinstance(instr, Gate):
-                for q in instr.qubits:
-                    if not isinstance(q, QubitPlaceholder):
-                        used_indexes.add(q.index)
-            elif isinstance(instr, Measurement):
-                if not isinstance(instr.qubit, QubitPlaceholder):
-                    used_indexes.add(instr.qubit.index)
-
-        def find_available_index():
-            # Just do a linear search.
-            for i in count(start=0, step=1):
-                if i not in used_indexes:
-                    return i
-
-        qubit_mapping = dict()
-
-        def remap_qubit(qubit):
-            if not isinstance(qubit, QubitPlaceholder):
-                return qubit
-            if id(qubit) in qubit_mapping:
-                return qubit_mapping[id(qubit)]
-            else:
-                available_index = find_available_index()
-                used_indexes.add(available_index)
-                remapped_qubit = Qubit(available_index)
-                qubit_mapping[id(qubit)] = remapped_qubit
-                return remapped_qubit
-
-        label_mapping = dict()
-        label_counter = 1
-
-        def remap_label(placeholder):
-            if id(placeholder) in label_mapping:
-                return label_mapping[id(placeholder)]
-            else:
-                label = Label(placeholder.prefix + str(label_counter))
-                label_mapping[id(placeholder)] = label
-                return label
-
-        result = []
-        for instr in self._instructions:
-            # Remap qubits on Gate and Measurement instructions
-            if isinstance(instr, Gate):
-                remapped_qubits = [remap_qubit(q) for q in instr.qubits]
-                result.append(Gate(instr.name, instr.params, remapped_qubits))
-            elif isinstance(instr, Measurement):
-                result.append(Measurement(remap_qubit(instr.qubit), instr.classical_reg))
-
-            # Remap any label placeholders on jump or target instructions
-            elif isinstance(instr, Jump) and isinstance(instr.target, LabelPlaceholder):
-                result.append(Jump(remap_label(instr.target)))
-                label_counter += 1
-            elif isinstance(instr, JumpTarget) and isinstance(instr.label, LabelPlaceholder):
-                result.append(JumpTarget(remap_label(instr.label)))
-                label_counter += 1
-            elif isinstance(instr, JumpConditional) and isinstance(instr.target, LabelPlaceholder):
-                new_label = remap_label(instr.target)
-                if isinstance(instr, JumpWhen):
-                    result.append(JumpWhen(new_label, instr.condition))
-                elif isinstance(instr, JumpUnless):
-                    result.append(JumpUnless(new_label, instr.condition))
-                else:
-                    raise TypeError("Encountered a JumpConditional that wasn't JumpWhen or JumpUnless: {} {}"
-                                    .format(type(instr), instr))
-                label_counter += 1
-
-            # Otherwise simply add it to the result
-            else:
-                result.append(instr)
-
-        return result
+        self._synthesized_instructions = instantiate_labels(self._instructions)
+        return self
 
     def __add__(self, other):
         """
@@ -589,7 +526,171 @@ class Program(object):
         return len(self._instructions)
 
     def __str__(self):
-        return self.out()
+        """
+        A string representation of the Quil program for inspection.
+
+        This may not be suitable for submission to a QPU or QVM for example if
+        your program contains unaddressed QubitPlaceholders
+        """
+        return '\n'.join(itertools.chain(
+            (str(dg) for dg in self._defined_gates),
+            (str(instr) for instr in self.instructions),
+            [''],
+        ))
+
+
+def _what_type_of_qubit_does_it_use(program):
+    """Helper function to peruse through a program's qubits.
+
+    This function will also enforce the condition that a Program uses either all placeholders
+    or all instantiated qubits to avoid accidentally mixing the two. This function will warn
+    if your program doesn't use any qubits.
+
+    :return: tuple of (whether the program uses placeholder qubits, whether the program uses
+        real qubits, a list of qubits ordered by their first appearance in the program)
+    """
+    has_placeholders = False
+    has_real_qubits = False
+
+    # We probably want to index qubits in the order they are encountered in the program
+    # so an ordered set would be nice. Python doesn't *have* an ordered set. Use the keys
+    # of an ordered dictionary instead
+    qubits = OrderedDict()
+
+    for instr in program:
+        if isinstance(instr, Gate):
+            for q in instr.qubits:
+                qubits[q] = 1
+                if isinstance(q, QubitPlaceholder):
+                    has_placeholders = True
+                elif isinstance(q, Qubit):
+                    has_real_qubits = True
+                else:
+                    raise ValueError("Unknown qubit type {}".format(q))
+        elif isinstance(instr, Measurement):
+            qubits[instr.qubit] = 1
+            if isinstance(instr.qubit, QubitPlaceholder):
+                has_placeholders = True
+            elif isinstance(instr.qubit, Qubit):
+                has_real_qubits = True
+            else:
+                raise ValueError("Unknown qubit type {}".format(instr.qubit))
+
+    if not (has_placeholders or has_real_qubits):
+        warnings.warn("Your program doesn't use any qubits")
+
+    if has_placeholders and has_real_qubits:
+        raise ValueError("Your program mixes instantiated qubits with placeholders")
+
+    return has_placeholders, has_real_qubits, list(qubits.keys())
+
+
+def get_default_qubit_mapping(program):
+    """
+    Takes a program which contains qubit placeholders and provides a mapping to the integers
+    0 through N-1.
+
+    The output of this function is suitable for input to :py:func:`address_qubits`.
+
+    :param program: A program containing qubit placeholders
+    :return: A dictionary mapping qubit placeholder to an addressed qubit from 0 through N-1.
+    """
+    fake_qubits, real_qubits, qubits = _what_type_of_qubit_does_it_use(program)
+    if real_qubits:
+        warnings.warn("This program contains integer qubits, "
+                      "so getting a mapping doesn't make sense.")
+        return {q: q for q in qubits}
+    return {qp: Qubit(i) for i, qp in enumerate(qubits)}
+
+
+def address_qubits(program, qubit_mapping=None):
+    """
+    Takes a program which contains placeholders and assigns them all defined values.
+
+    Either all qubits must be defined or all undefined. If qubits are
+    undefined, you may provide a qubit mapping to specify how placeholders get mapped
+    to actual qubits. If a mapping is not provided, integers 0 through N are used.
+
+    This function will also instantiate any label placeholders.
+
+    :param program: The program.
+    :param qubit_mapping: A dictionary-like object that maps from :py:class:`QubitPlaceholder`
+        to :py:class:`Qubit` or ``int`` (but not both).
+    :return: A new Program with all qubit and label placeholders assigned to real qubits and labels.
+    """
+    fake_qubits, real_qubits, qubits = _what_type_of_qubit_does_it_use(program)
+    if real_qubits:
+        if qubit_mapping is not None:
+            warnings.warn("A qubit mapping was provided but the program does not "
+                          "contain any placeholders to map!")
+        return program
+
+    if qubit_mapping is None:
+        qubit_mapping = {qp: Qubit(i) for i, qp in enumerate(qubits)}
+    else:
+        if all(isinstance(v, Qubit) for v in qubit_mapping.values()):
+            pass  # we good
+        elif all(isinstance(v, int) for v in qubit_mapping.values()):
+            qubit_mapping = {k: Qubit(v) for k, v in qubit_mapping.items()}
+        else:
+            raise ValueError("Qubit mapping must map to type Qubit or int (but not both)")
+
+    result = []
+    for instr in program:
+        # Remap qubits on Gate and Measurement instructions
+        if isinstance(instr, Gate):
+            remapped_qubits = [qubit_mapping[q] for q in instr.qubits]
+            result.append(Gate(instr.name, instr.params, remapped_qubits))
+        elif isinstance(instr, Measurement):
+            result.append(Measurement(qubit_mapping[instr.qubit], instr.classical_reg))
+
+        # Otherwise simply add it to the result
+        else:
+            result.append(instr)
+
+    return Program(result)
+
+
+def _get_label(placeholder, label_mapping, label_i):
+    """Helper function to either get the appropriate label for a given placeholder or generate
+    a new label and update the mapping.
+
+    See :py:func:`instantiate_labels` for usage.
+    """
+    if placeholder in label_mapping:
+        return label_mapping[placeholder], label_mapping, label_i
+
+    new_target = Label("{}{}".format(placeholder.prefix, label_i))
+    label_i += 1
+    label_mapping[placeholder] = new_target
+    return new_target, label_mapping, label_i
+
+
+def instantiate_labels(instructions):
+    """
+    Takes an iterable of instructions which may contain label placeholders and assigns
+    them all defined values.
+
+    :return: list of instructions with all label placeholders assigned to real labels.
+    """
+    label_i = 1
+    result = []
+    label_mapping = dict()
+    for instr in instructions:
+        if isinstance(instr, Jump) and isinstance(instr.target, LabelPlaceholder):
+            new_target, label_mapping, label_i = _get_label(instr.target, label_mapping, label_i)
+            result.append(Jump(new_target))
+        elif isinstance(instr, JumpConditional) and isinstance(instr.target, LabelPlaceholder):
+            new_target, label_mapping, label_i = _get_label(instr.target, label_mapping, label_i)
+            cls = instr.__class__  # Make the correct subclass
+            result.append(cls(new_target, instr.condition))
+        elif isinstance(instr, JumpTarget) and isinstance(instr.label, LabelPlaceholder):
+            new_label, label_mapping, label_i = _get_label(instr.label, label_mapping, label_i)
+            result.append(JumpTarget(new_label))
+        else:
+            result.append(instr)
+
+    return result
 
 
 def merge_programs(prog_list):
@@ -603,26 +704,6 @@ def merge_programs(prog_list):
     return sum(prog_list, Program())
 
 
-def shift_quantum_gates(program, shift_offset):
-    """
-    Shifts a quantum gates in a quil program so that all qubit indices change by a certain offset
-    :param program: a pyquil Program
-    :param shift_offset: integer
-    :return: pyquil Program with shifted qubit indices
-    """
-    if not isinstance(shift_offset, int):
-        raise ValueError("shift_offset must be an integer")
-    if not isinstance(program, Program):
-        raise ValueError("Program must be a pyquil Program instance")
-    shifted_program = Program()
-    shifted_program.inst(program)
-    for instruct in shifted_program:
-        if isinstance(instruct, Gate):
-            for qubit in instruct.qubits:
-                qubit.index += shift_offset
-    return shifted_program
-
-
 def get_classical_addresses_from_program(program):
     """
     Returns a sorted list of classical addresses found in the MEASURE instructions in the program.
diff --git a/pyquil/quilatom.py b/pyquil/quilatom.py
index d496f3e..13f3094 100644
--- a/pyquil/quilatom.py
+++ b/pyquil/quilatom.py
@@ -5,20 +5,21 @@ class QuilAtom(object):
     """
     Abstract class for atomic elements of Quil.
     """
+
     def out(self):
-        pass
+        raise NotImplementedError()
 
     def __str__(self):
-        return self.out()
+        raise NotImplementedError()
 
     def __eq__(self, other):
-        return isinstance(other, self.__class__) and self.out() == other.out()
+        raise NotImplementedError()
 
     def __ne__(self, other):
         return not self.__eq__(other)
 
     def __hash__(self):
-        return hash(str(self))
+        raise NotImplementedError()
 
 
 class Qubit(QuilAtom):
@@ -27,6 +28,7 @@ class Qubit(QuilAtom):
 
     :param int index: Index of the qubit.
     """
+
     def __init__(self, index):
         if not (isinstance(index, integer_types) and index >= 0):
             raise TypeError("Addr index must be a non-negative int")
@@ -35,17 +37,22 @@ class Qubit(QuilAtom):
     def out(self):
         return str(self.index)
 
+    def __str__(self):
+        return str(self.index)
+
     def __repr__(self):
         return "<Qubit {0}>".format(self.index)
 
+    def __hash__(self):
+        return hash(self.index)
+
+    def __eq__(self, other):
+        return isinstance(other, Qubit) and other.index == self.index
 
-class QubitPlaceholder(Qubit):
-    def __init__(self):
-        pass
 
-    @property
-    def index(self):
-        raise RuntimeError("Qubit has not been assigned an index")
+class QubitPlaceholder(QuilAtom):
+    def out(self):
+        raise RuntimeError("Qubit {} has not been assigned an index".format(self))
 
     def __str__(self):
         return repr(self)
@@ -53,6 +60,31 @@ class QubitPlaceholder(Qubit):
     def __repr__(self):
         return "<QubitPlaceholder {}>".format(id(self))
 
+    def __hash__(self):
+        return hash(id(self))
+
+    def __eq__(self, other):
+        return isinstance(other, QubitPlaceholder) and id(other) == id(self)
+
+    @classmethod
+    def register(cls, n):
+        """Return a 'register' of ``n`` QubitPlaceholders.
+
+        >>> qs = QubitPlaceholder.register(8) # a qubyte
+        >>> prog = Program(H(q) for q in qs)
+        >>> address_qubits(prog).out()
+        H 0
+        H 1
+        ...
+        >>>
+
+        The returned register is a Python list of QubitPlaceholder objects, so all
+        normal list semantics apply.
+
+        :param n: The number of qubits in the register
+        """
+        return [cls() for _ in range(n)]
+
 
 def unpack_qubit(qubit):
     """
@@ -65,6 +97,8 @@ def unpack_qubit(qubit):
         return Qubit(qubit)
     elif isinstance(qubit, Qubit):
         return qubit
+    elif isinstance(qubit, QubitPlaceholder):
+        return qubit
     else:
         raise TypeError("qubit should be an int or Qubit instance")
 
@@ -101,11 +135,20 @@ class Addr(QuilAtom):
         self.address = value
 
     def out(self):
-        return "[{0}]".format(self.address)
+        return "[{}]".format(self.address)
+
+    def __str__(self):
+        return "[{}]".format(self.address)
 
     def __repr__(self):
         return "<Addr {0}>".format(self.address)
 
+    def __eq__(self, other):
+        return isinstance(other, Addr) and other.address == self.address
+
+    def __hash__(self):
+        return hash(self.address)
+
 
 class Label(QuilAtom):
     """
@@ -118,18 +161,26 @@ class Label(QuilAtom):
         self.name = label_name
 
     def out(self):
-        return "@" + str(self.name)
+        return "@{name}".format(name=self.name)
+
+    def __str__(self):
+        return "@{name}".format(name=self.name)
 
     def __repr__(self):
         return "<Label {0}>".format(repr(self.name))
 
+    def __eq__(self, other):
+        return isinstance(other, Label) and other.name == self.name
 
-class LabelPlaceholder(Label):
+    def __hash__(self):
+        return hash(self.name)
+
+
+class LabelPlaceholder(QuilAtom):
     def __init__(self, prefix="L"):
         self.prefix = prefix
 
-    @property
-    def name(self):
+    def out(self):
         raise RuntimeError("Label has not been assigned a name")
 
     def __str__(self):
@@ -137,3 +188,9 @@ class LabelPlaceholder(Label):
 
     def __repr__(self):
         return "<LabelPlaceholder {} {}>".format(self.prefix, id(self))
+
+    def __eq__(self, other):
+        return isinstance(other, LabelPlaceholder) and id(other) == id(self)
+
+    def __hash__(self):
+        return hash(id(self))
diff --git a/pyquil/quilbase.py b/pyquil/quilbase.py
index e9c6bd5..64a2f2c 100644
--- a/pyquil/quilbase.py
+++ b/pyquil/quilbase.py
@@ -21,7 +21,7 @@ import numpy as np
 from six import integer_types, string_types
 
 from pyquil.parameters import Expression, _contained_parameters, format_parameter
-from pyquil.quilatom import Qubit, Addr, Label, unpack_qubit
+from pyquil.quilatom import Qubit, Addr, Label, unpack_qubit, QubitPlaceholder, LabelPlaceholder
 
 
 class AbstractInstruction(object):
@@ -51,6 +51,12 @@ RESERVED_WORDS = ['DEFGATE', 'DEFCIRCUIT', 'MEASURE',
                   'FALSE', 'TRUE', 'NOT', 'AND', 'OR', 'MOVE', 'EXCHANGE']
 
 
+def _extract_qubit_index(qubit, index=True):
+    if (not index) or isinstance(qubit, QubitPlaceholder):
+        return qubit
+    return qubit.index
+
+
 class Gate(AbstractInstruction):
     """
     This is the pyQuil object for a quantum gate instruction.
@@ -69,19 +75,22 @@ class Gate(AbstractInstruction):
         if not isinstance(qubits, list) or not qubits:
             raise TypeError("Gate arguments must be a non-empty list")
         for qubit in qubits:
-            if not isinstance(qubit, Qubit):
+            if not isinstance(qubit, (Qubit, QubitPlaceholder)):
                 raise TypeError("Gate arguments must all be Qubits")
 
         self.name = name
         self.params = params
         self.qubits = qubits
 
+    def get_qubits(self, indices=True):
+        return {_extract_qubit_index(q, indices) for q in self.qubits}
+
     def out(self):
         def format_params(params):
             return "(" + ",".join(map(format_parameter, params)) + ")"
 
         def format_qubits(qubits):
-            return " ".join([str(qubit) for qubit in qubits])
+            return " ".join([qubit.out() for qubit in qubits])
 
         if self.params:
             return "{}{} {}".format(self.name, format_params(self.params), format_qubits(self.qubits))
@@ -89,7 +98,19 @@ class Gate(AbstractInstruction):
             return "{} {}".format(self.name, format_qubits(self.qubits))
 
     def __repr__(self):
-        return "<Gate " + self.out() + ">"
+        return "<Gate " + str(self) + ">"
+
+    def __str__(self):
+        def format_params(params):
+            return "(" + ",".join(map(format_parameter, params)) + ")"
+
+        def format_qubits(qubits):
+            return " ".join([str(qubit) for qubit in qubits])
+
+        if self.params:
+            return "{}{} {}".format(self.name, format_params(self.params), format_qubits(self.qubits))
+        else:
+            return "{} {}".format(self.name, format_qubits(self.qubits))
 
 
 class Measurement(AbstractInstruction):
@@ -98,7 +119,7 @@ class Measurement(AbstractInstruction):
     """
 
     def __init__(self, qubit, classical_reg=None):
-        if not isinstance(qubit, Qubit):
+        if not isinstance(qubit, (Qubit, QubitPlaceholder)):
             raise TypeError("qubit should be a Qubit")
         if classical_reg and not isinstance(classical_reg, Addr):
             raise TypeError("classical_reg should be None or an Addr instance")
@@ -112,6 +133,9 @@ class Measurement(AbstractInstruction):
         else:
             return "MEASURE {}".format(self.qubit)
 
+    def get_qubits(self, indices=True):
+        return {_extract_qubit_index(self.qubit, indices)}
+
 
 class DefGate(AbstractInstruction):
     """
@@ -222,7 +246,7 @@ class JumpTarget(AbstractInstruction):
     """
 
     def __init__(self, label):
-        if not isinstance(label, Label):
+        if not isinstance(label, (Label, LabelPlaceholder)):
             raise TypeError("label must be a Label")
         self.label = label
 
@@ -239,7 +263,7 @@ class JumpConditional(AbstractInstruction):
     """
 
     def __init__(self, target, condition):
-        if not isinstance(target, Label):
+        if not isinstance(target, (Label, LabelPlaceholder)):
             raise TypeError("target should be a Label")
         if not isinstance(condition, Addr):
             raise TypeError("condition should be an Addr")
@@ -366,7 +390,7 @@ class Jump(AbstractInstruction):
     """
 
     def __init__(self, target):
-        if not isinstance(target, Label):
+        if not isinstance(target, (Label, LabelPlaceholder)):
             raise TypeError("target should be a Label")
         self.target = target
 


Test patch:
diff --git a/pyquil/tests/test_parameters.py b/pyquil/tests/test_parameters.py
index c9298d2..7921057 100644
--- a/pyquil/tests/test_parameters.py
+++ b/pyquil/tests/test_parameters.py
@@ -14,6 +14,8 @@ def test_format_parameter():
         (1j, 'i'),
         (0 + 1j, 'i'),
         (-1j, '-i'),
+        (1e-15 + 1j, 'i'),
+        (1e-15 - 1j, '-i')
     ]
 
     for test_case in test_cases:
diff --git a/pyquil/tests/test_paulis_with_placeholders.py b/pyquil/tests/test_paulis_with_placeholders.py
new file mode 100644
index 0000000..d3b8661
--- /dev/null
+++ b/pyquil/tests/test_paulis_with_placeholders.py
@@ -0,0 +1,610 @@
+##############################################################################
+# Copyright 2016-2018 Rigetti Computing
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+##############################################################################
+
+import math
+import warnings
+from functools import reduce
+from itertools import product
+from operator import mul
+
+import numpy as np
+import pytest
+import re
+from six.moves import range
+
+from pyquil.gates import RX, RZ, CNOT, H, X, PHASE
+from pyquil.paulis import PauliTerm, PauliSum, exponential_map, exponentiate_commuting_pauli_sum, \
+    ID, UnequalLengthWarning, exponentiate, trotterize, is_zero, check_commutation, commuting_sets, \
+    term_with_coeff, sI, sX, sY, sZ, ZERO, is_identity
+from pyquil.quil import Program, address_qubits, get_default_qubit_mapping
+from pyquil.quilatom import QubitPlaceholder, Qubit
+
+
+def test_simplify_terms():
+    q = QubitPlaceholder.register(1)
+    term = PauliTerm('Z', q[0]) * -1.0 * PauliTerm('Z', q[0])
+    assert term.id() == ''
+    assert term.coefficient == -1.0
+
+    term = PauliTerm('Z', q[0]) + PauliTerm('Z', q[0], 1.0)
+    assert str(term).startswith('(2+0j)*Z<QubitPlaceholder ')
+
+
+def test_get_qubits():
+    q = QubitPlaceholder.register(2)
+    term = PauliTerm('Z', q[0]) * PauliTerm('X', q[1])
+    assert term.get_qubits() == q
+
+    q10 = QubitPlaceholder()
+    sum_term = PauliTerm('X', q[0], 0.5) + 0.5j * PauliTerm('Y', q10) * PauliTerm('Y', q[0], 0.5j)
+    assert sum_term.get_qubits() == [q[0], q10]
+
+
+def test_simplify_term_single():
+    q0, q1, q2 = QubitPlaceholder.register(3)
+    term = (PauliTerm('Z', q0) * PauliTerm('I', q1) *
+            PauliTerm('X', q2, 0.5j) * PauliTerm('Z', q0, 1.0))
+    assert term.id() == 'X{}'.format(q2)
+    assert term.coefficient == 0.5j
+
+
+def test_simplify_term_xz():
+    q0 = QubitPlaceholder()
+    term1 = (-0.5 * PauliTerm('X', q0)) * (-1.0 * PauliTerm('Z', q0))
+    term2 = -0.5 * PauliTerm('X', q0) * (-1.0) * PauliTerm('Z', q0)
+    term3 = 0.5 * PauliTerm('X', q0) * PauliTerm('Z', q0)
+    for term in [term1, term2, term3]:
+        assert term.id() == 'Y{}'.format(q0)
+        assert term.coefficient == -0.5j
+
+
+def test_simplify_term_multindex():
+    q0, q2 = QubitPlaceholder.register(2)
+    term = (PauliTerm('X', q0, coefficient=-0.5) * PauliTerm('Z', q0, coefficient=-1.0) *
+            PauliTerm('X', q2, 0.5))
+    assert term.id(sort_ops=False) == 'Y{q0}X{q2}'.format(q0=q0, q2=q2)
+    assert term.coefficient == -0.25j
+
+
+def test_simplify_sum_terms():
+    q0 = QubitPlaceholder()
+    sum_term = PauliSum([PauliTerm('X', q0, 0.5), PauliTerm('Z', q0, 0.5j)])
+    str_sum_term = str(sum_term + sum_term)
+    assert (str_sum_term == '(1+0j)*X{q0} + 1j*Z{q0}'.format(q0=q0) or
+            str_sum_term == '1j*Z{q0} + (1+0j)*X{q0}'.format(q0=q0))
+    sum_term = PauliSum([PauliTerm('X', q0, 0.5), PauliTerm('X', q0, 0.5)])
+    assert str(sum_term.simplify()) == '(1+0j)*X{q0}'.format(q0=q0)
+
+    # test the simplify on multiplication
+    sum_term = PauliSum([PauliTerm('X', q0, 0.5), PauliTerm('X', q0, 0.5)])
+    assert str(sum_term * sum_term) == '(1+0j)*I'
+
+
+def test_copy():
+    q0, q1 = QubitPlaceholder.register(2)
+    term = PauliTerm('X', q0, 0.5) * PauliTerm('X', q1, 0.5)
+    new_term = term.copy()
+
+    q2 = QubitPlaceholder()
+    term = term * PauliTerm('X', q2, 0.5)
+    new_term = new_term * PauliTerm('X', q2, 0.5)
+
+    assert term == new_term  # value equality
+    assert term is not new_term  # ref inequality
+    assert term._ops is not new_term._ops
+
+    term = PauliTerm('X', q0, 0.5) * PauliTerm('X', q1, 0.5)
+    new_term = term * PauliTerm('X', q2, 0.5)
+    assert term != new_term
+    assert term is not new_term
+    assert term._ops is not new_term._ops
+
+
+def test_len():
+    q0, q1 = QubitPlaceholder.register(2)
+    term = PauliTerm("Z", q0, 1.0) * PauliTerm("Z", q1, 1.0)
+    assert len(term) == 2
+
+
+def test_sum_len():
+    q0, q1 = QubitPlaceholder.register(2)
+    pauli_sum = PauliTerm("Z", q0, 1.0) + PauliTerm("Z", q1, 1.0)
+    assert len(pauli_sum) == 2
+
+
+def test_enumerate():
+    q0, q1, q5 = QubitPlaceholder.register(3)
+    term = PauliTerm("Z", q0, 1.0) * PauliTerm("Z", q1, 1.0) * PauliTerm("X", q5, 5)
+    position_op_pairs = [(q0, "Z"), (q1, "Z"), (q5, "X")]
+    for key, val in term:
+        assert (key, val) in position_op_pairs
+
+
+def test_getitem():
+    q = QubitPlaceholder.register(6)
+    term = PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0) * PauliTerm("X", q[5], 5)
+    assert term[q[0]] == "Z"
+    assert term[q[1]] == "Z"
+    assert term[q[2]] == "I"
+    assert term[q[3]] == "I"
+    assert term[q[4]] == "I"
+    assert term[q[5]] == "X"
+    assert len(term) == 3
+
+
+def test_ids():
+    q = QubitPlaceholder.register(6)
+    term_1 = PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0) * PauliTerm("X", q[5], 5)
+    term_2 = PauliTerm("X", q[5], 5) * PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    with pytest.raises(TypeError):
+        # Not sortable
+        t = term_1.id() == term_2.id()
+
+
+def test_ids_no_sort():
+    q = QubitPlaceholder.register(6)
+    term_1 = PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0) * PauliTerm("X", q[5], 5)
+    term_2 = PauliTerm("X", q[5], 5) * PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    assert re.match('Z.+Z.+X.+', term_1.id(sort_ops=False))
+    assert re.match('X.+Z.+Z.+', term_2.id(sort_ops=False))
+
+
+def test_operations_as_set():
+    q = QubitPlaceholder.register(6)
+    term_1 = PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0) * PauliTerm("X", q[5], 5)
+    term_2 = PauliTerm("X", q[5], 5) * PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    assert term_1.operations_as_set() == term_2.operations_as_set()
+
+
+def test_pauli_sum():
+    q = QubitPlaceholder.register(8)
+    q_plus = 0.5 * PauliTerm('X', q[0]) + 0.5j * PauliTerm('Y', q[0])
+    the_sum = q_plus * PauliSum([PauliTerm('X', q[0])])
+    term_strings = [str(x) for x in the_sum.terms]
+    assert '(0.5+0j)*I' in term_strings
+    assert len(term_strings) == 2
+    assert len(the_sum.terms) == 2
+
+    the_sum = q_plus * PauliTerm('X', q[0])
+    term_strings = [str(x) for x in the_sum.terms]
+    assert '(0.5+0j)*I' in term_strings
+    assert len(term_strings) == 2
+    assert len(the_sum.terms) == 2
+
+    the_sum = PauliTerm('X', q[0]) * q_plus
+    term_strings = [str(x) for x in the_sum.terms]
+    assert '(0.5+0j)*I' in term_strings
+    assert len(term_strings) == 2
+    assert len(the_sum.terms) == 2
+
+    with pytest.raises(ValueError):
+        _ = PauliSum(sI(q[0]))
+    with pytest.raises(ValueError):
+        _ = PauliSum([1, 1, 1, 1])
+    with pytest.raises(ValueError):
+        _ = the_sum * []
+
+
+def test_pauliterm_sub():
+    q = QubitPlaceholder.register(8)
+    assert str(sX(q[1]) - 2.0) == str(sX(q[1]) + -2.0)
+    assert str(1.4 - sZ(q[1])) == str(1.4 + -1.0 * sZ(q[1]))
+
+
+def test_ps_sub():
+    q0 = QubitPlaceholder()
+    term = 3 * ID()
+    b = term - 1.0
+    assert str(b) == "(2+0j)*I"
+    assert str(b - 1.0) == "(1+0j)*I"
+    assert str(1.0 - b) == "(-1+0j)*I"
+    b = 1.0 - sX(q0)
+    assert re.match(r"\(1\+0j\)\*I \+ \(-1\+0j\)\*X<.+>", str(b))
+    b = sX(q0) - 1.0
+    assert re.match(r"\(1\+0j\)\*X<.+> \+ \(-1\+0j\)\*I", str(b))
+
+
+def test_exponentiate_1():
+    # test rotation of single qubit
+    q = QubitPlaceholder.register(8)
+    generator = PauliTerm("Z", q[0], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst(RZ(2.0)(q[0]))
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_2():
+    # testing general 2-circuit
+    q = QubitPlaceholder.register(8)
+    generator = PauliTerm("Z", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst(CNOT(q[0], q[1])).inst(RZ(2.0)(q[1])).inst(CNOT(q[0], q[1]))
+
+    mapping = get_default_qubit_mapping(prog)
+    assert address_qubits(prog, mapping) == address_qubits(result_prog, mapping)
+
+
+def test_exponentiate_bp0_ZX():
+    q = QubitPlaceholder.register(8)
+    # testing change of basis position 0
+    generator = PauliTerm("X", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    param_prog = exponential_map(generator)
+    prog = param_prog(1)
+    result_prog = Program().inst(
+        [H(q[0]), CNOT(q[0], q[1]), RZ(2.0)(q[1]), CNOT(q[0], q[1]), H(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_bp1_XZ():
+    # testing change of basis position 1
+    q = QubitPlaceholder.register(8)
+    generator = PauliTerm("Z", q[0], 1.0) * PauliTerm("X", q[1], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst(
+        [H(q[1]), CNOT(q[0], q[1]), RZ(2.0)(q[1]), CNOT(q[0], q[1]), H(q[1])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_bp0_ZY():
+    # testing change of basis position 0
+    q = QubitPlaceholder.register(8)
+    generator = PauliTerm("Y", q[0], 1.0) * PauliTerm("Z", q[1], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([RX(math.pi / 2.0)(q[0]), CNOT(q[0], q[1]), RZ(2.0)(q[1]),
+                                  CNOT(q[0], q[1]), RX(-math.pi / 2)(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_bp1_YZ():
+    q = QubitPlaceholder.register(8)
+    # testing change of basis position 1
+    generator = PauliTerm("Z", q[0], 1.0) * PauliTerm("Y", q[1], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([RX(math.pi / 2.0)(q[1]), CNOT(q[0], q[1]),
+                                  RZ(2.0)(q[1]), CNOT(q[0], q[1]), RX(-math.pi / 2.0)(q[1])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_3cob():
+    # testing circuit for 3-terms with change of basis
+    q = QubitPlaceholder.register(8)
+    generator = PauliTerm("Z", q[0], 1.0) * PauliTerm("Y", q[1], 1.0) * PauliTerm("X", q[2], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([RX(math.pi / 2.0)(q[1]), H(q[2]), CNOT(q[0], q[1]),
+                                  CNOT(q[1], q[2]), RZ(2.0)(q[2]), CNOT(q[1], q[2]),
+                                  CNOT(q[0], q[1]), RX(-math.pi / 2.0)(q[1]), H(q[2])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_3ns():
+    # testing circuit for 3-terms non-sequential
+    q = QubitPlaceholder.register(8)
+    generator = (PauliTerm("Y", q[0], 1.0) *
+                 PauliTerm("I", q[1], 1.0) *
+                 PauliTerm("Y", q[2], 1.0) *
+                 PauliTerm("Y", q[3], 1.0))
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([RX(math.pi / 2.0)(q[0]), RX(math.pi / 2.0)(q[2]),
+                                  RX(math.pi / 2.0)(q[3]), CNOT(q[0], q[2]),
+                                  CNOT(q[2], q[3]), RZ(2.0)(q[3]), CNOT(q[2], q[3]),
+                                  CNOT(q[0], q[2]), RX(-math.pi / 2.0)(q[0]),
+                                  RX(-math.pi / 2.0)(q[2]), RX(-math.pi / 2.0)(q[3])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_commuting_pauli_sum():
+    q = QubitPlaceholder.register(8)
+    pauli_sum = PauliSum([PauliTerm('Z', q[0], 0.5), PauliTerm('Z', q[1], 0.5)])
+    prog = Program().inst(RZ(1.)(q[0])).inst(RZ(1.)(q[1]))
+    result_prog = exponentiate_commuting_pauli_sum(pauli_sum)(1.)
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_prog():
+    q = QubitPlaceholder.register(8)
+    ham = PauliTerm("Z", q[0])
+    result_prog = Program(RZ(2.0, q[0]))
+    prog = exponentiate(ham)
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_exponentiate_identity():
+    q = QubitPlaceholder.register(11)
+    mapping = {qp: Qubit(i) for i, qp in enumerate(q)}
+    generator = PauliTerm("I", q[1], 0.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([X(q[0]), PHASE(-0.0)(q[0]), X(q[0]), PHASE(-0.0)(q[0])])
+    assert address_qubits(prog, mapping) == address_qubits(result_prog, mapping)
+
+    generator = PauliTerm("I", q[1], 1.0)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([X(q[0]), PHASE(-1.0)(q[0]), X(q[0]), PHASE(-1.0)(q[0])])
+    assert address_qubits(prog, mapping) == address_qubits(result_prog, mapping)
+
+    generator = PauliTerm("I", q[10], 0.08)
+    para_prog = exponential_map(generator)
+    prog = para_prog(1)
+    result_prog = Program().inst([X(q[0]), PHASE(-0.08)(q[0]), X(q[0]), PHASE(-0.08)(q[0])])
+    assert address_qubits(prog, mapping) == address_qubits(result_prog, mapping)
+
+
+def test_trotterize():
+    q = QubitPlaceholder.register(8)
+    term_one = PauliTerm("X", q[0], 1.0)
+    term_two = PauliTerm("Z", q[0], 1.0)
+
+    with pytest.raises(ValueError):
+        trotterize(term_one, term_two, trotter_order=0)
+    with pytest.raises(ValueError):
+        trotterize(term_one, term_two, trotter_order=5)
+
+    prog = trotterize(term_one, term_one)
+    result_prog = Program().inst([H(q[0]), RZ(2.0)(q[0]), H(q[0]), H(q[0]),
+                                  RZ(2.0)(q[0]), H(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+    # trotter_order 1 steps 1
+    prog = trotterize(term_one, term_two, trotter_steps=1)
+    result_prog = Program().inst([H(q[0]), RZ(2.0)(q[0]), H(q[0]), RZ(2.0)(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+    # trotter_order 1 steps 2
+    prog = trotterize(term_one, term_two, trotter_steps=2)
+    result_prog = Program().inst([H(q[0]), RZ(1.0)(q[0]), H(q[0]), RZ(1.0)(q[0]),
+                                  H(q[0]), RZ(1.0)(q[0]), H(q[0]), RZ(1.0)(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+    # trotter_order 2 steps 1
+    prog = trotterize(term_one, term_two, trotter_order=2)
+    result_prog = Program().inst([H(q[0]), RZ(1.0)(q[0]), H(q[0]), RZ(2.0)(q[0]),
+                                  H(q[0]), RZ(1.0)(q[0]), H(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+    # trotter_order 2 steps 2
+    prog = trotterize(term_one, term_two, trotter_order=2, trotter_steps=2)
+    result_prog = Program().inst([H(q[0]), RZ(0.5)(q[0]), H(q[0]), RZ(1.0)(q[0]),
+                                  H(q[0]), RZ(0.5)(q[0]), H(q[0]),
+                                  H(q[0]), RZ(0.5)(q[0]), H(q[0]), RZ(1.0)(q[0]),
+                                  H(q[0]), RZ(0.5)(q[0]), H(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+    # trotter_order 3 steps 1
+    prog = trotterize(term_one, term_two, trotter_order=3, trotter_steps=1)
+    result_prog = Program().inst([H(q[0]), RZ(14.0 / 24)(q[0]), H(q[0]), RZ(4.0 / 3.0)(q[0]),
+                                  H(q[0]), RZ(1.5)(q[0]), H(q[0]), RZ(-4.0 / 3.0)(q[0]),
+                                  H(q[0]), RZ(-2.0 / 24)(q[0]), H(q[0]), RZ(2.0)(q[0])])
+    assert address_qubits(prog) == address_qubits(result_prog)
+
+
+def test_is_zeron():
+    q = QubitPlaceholder.register(8)
+    with pytest.raises(TypeError):
+        is_zero(1)
+
+    p_term = PauliTerm("X", q[0])
+    ps_term = p_term + PauliTerm("Z", q[1])
+
+    assert not is_zero(p_term)
+    assert is_zero(p_term + -1 * p_term)
+    assert not is_zero(ps_term)
+
+
+def test_check_commutation():
+    q = QubitPlaceholder.register(8)
+    term1 = PauliTerm("X", q[0]) * PauliTerm("X", q[1])
+    term2 = PauliTerm("Y", q[0]) * PauliTerm("Y", q[1])
+    term3 = PauliTerm("Y", q[0]) * PauliTerm("Z", q[2])
+    # assert check_commutation(PauliSum([term1]), term2)
+    assert check_commutation([term2], term3)
+    assert check_commutation([term2], term3)
+    assert not check_commutation([term1], term3)
+
+
+def _commutator(t1, t2):
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore',
+                                message=r"The term .+ will be combined with .+, "
+                                        r"but they have different orders of operations.*",
+                                category=UserWarning)
+        return t1 * t2 + -1 * t2 * t1
+
+
+def test_check_commutation_rigorous():
+    # more rigorous test.  Get all operators in Pauli group
+    p_n_group = ("I", "X", "Y", "Z")
+    pauli_list = list(product(p_n_group, repeat=3))
+    pauli_ops = [list(zip(x, QubitPlaceholder.register(3))) for x in pauli_list]
+    pauli_ops_pq = [reduce(mul, (PauliTerm(*x) for x in op)) for op in pauli_ops]
+
+    non_commuting_pairs = []
+    commuting_pairs = []
+    for x in range(len(pauli_ops_pq)):
+        for y in range(x, len(pauli_ops_pq)):
+
+            tmp_op = _commutator(pauli_ops_pq[x], pauli_ops_pq[y])
+            assert len(tmp_op.terms) == 1
+            if is_identity(tmp_op.terms[0]):
+                commuting_pairs.append((pauli_ops_pq[x], pauli_ops_pq[y]))
+            else:
+                non_commuting_pairs.append((pauli_ops_pq[x], pauli_ops_pq[y]))
+
+    # now that we have our sets let's check against our code.
+    for t1, t2 in non_commuting_pairs:
+        assert not check_commutation([t1], t2)
+
+    for t1, t2 in commuting_pairs:
+        assert check_commutation([t1], t2)
+
+
+def test_commuting_sets():
+    q = QubitPlaceholder.register(8)
+    term1 = PauliTerm("X", q[0]) * PauliTerm("X", q[1])
+    term2 = PauliTerm("Y", q[0]) * PauliTerm("Y", q[1])
+    term3 = PauliTerm("Y", q[0]) * PauliTerm("Z", q[2])
+    pauli_sum = term1 + term2 + term3
+    commuting_sets(pauli_sum)
+
+
+def test_paulisum_iteration():
+    q = QubitPlaceholder.register(8)
+    term_list = [sX(q[2]), sZ(q[4])]
+    pauli_sum = sum(term_list)
+    for ii, term in enumerate(pauli_sum):
+        assert term_list[ii] == term
+
+
+def test_paulisum_indexing():
+    q = QubitPlaceholder.register(8)
+    pauli_sum = 0.5 * sX(q[0]) + 0.1 * sZ(q[1])
+    assert pauli_sum[0] == 0.5 * sX(q[0])
+    for ii, term in enumerate(pauli_sum.terms):
+        assert pauli_sum[ii] == term
+
+
+def test_term_powers():
+    for qubit in QubitPlaceholder.register(2):
+        pauli_terms = [sI(qubit), sX(qubit), sY(qubit), sZ(qubit)]
+        for pauli_term in pauli_terms:
+            assert pauli_term ** 0 == sI(qubit)
+            assert pauli_term ** 1 == pauli_term
+            assert pauli_term ** 2 == sI(qubit)
+            assert pauli_term ** 3 == pauli_term
+    with pytest.raises(ValueError):
+        pauli_terms[0] ** -1
+
+
+def test_term_large_powers():
+    # Test to make sure large powers can be computed
+    q = QubitPlaceholder.register(2)
+    (PauliTerm('X', q[0], 2) * PauliTerm('Y', q[0], 2)) ** 400
+
+
+def test_sum_power():
+    q = QubitPlaceholder.register(8)
+    pauli_sum = (sY(q[0]) - sX(q[0])) * (1.0 / np.sqrt(2))
+    assert pauli_sum ** 2 == PauliSum([sI(q[0])])
+    with pytest.raises(ValueError):
+        _ = pauli_sum ** -1
+    pauli_sum = sI(q[0]) + sI(q[1])
+    assert pauli_sum ** 0 == sI(q[0])
+    # Test to make sure large powers can be computed
+    pauli_sum ** 400
+
+
+def test_term_equality():
+    q0, q10 = QubitPlaceholder.register(2)
+    with pytest.raises(TypeError):
+        sI(q0) != 0
+    assert sI(q0) == sI(q0)
+    assert PauliTerm('X', q10, 1 + 1.j) == PauliTerm('X', q10, 1 + 1.j)
+    assert PauliTerm('X', q10, 1 + 1.j) + PauliTerm('X', q10, 1 + 1.j) != PauliTerm('X', q10,
+                                                                                    1 + 1.j)
+    assert PauliTerm('X', q10, 1 + 1.j) != PauliTerm('X', q10, 1 + 1.j) + PauliTerm('X', q10,
+                                                                                    1 + 1.j)
+
+
+def test_term_with_coeff():
+    q0 = QubitPlaceholder()
+    assert PauliTerm('X', q0, 1.j) == term_with_coeff(sX(q0), 1.j)
+    assert PauliTerm('X', q0, -1.0) == term_with_coeff(sX(q0), -1)
+    with pytest.raises(ValueError):
+        term_with_coeff(sI(q0), None)
+
+
+def test_sum_equality():
+    q0, q1 = QubitPlaceholder.register(2)
+    pauli_sum = sY(q0) - sX(q0)
+    assert pauli_sum != 2 * pauli_sum
+    with pytest.warns(UnequalLengthWarning):
+        assert pauli_sum != pauli_sum + sZ(q0)
+    with pytest.warns(UnequalLengthWarning):
+        assert pauli_sum + sZ(q0) != pauli_sum
+    assert pauli_sum != sY(q1) - sX(q1)
+    assert pauli_sum == -1.0 * sX(q0) + sY(q0)
+    assert pauli_sum == pauli_sum * 1.0
+    with pytest.raises(TypeError):
+        assert pauli_sum != 0
+
+
+def test_zero_term():
+    qubit_id = QubitPlaceholder()
+    coefficient = 10
+    ps = sI(qubit_id) + sX(qubit_id)
+    assert coefficient * ZERO() == ZERO()
+    assert ZERO() * coefficient == ZERO()
+    assert ZERO() * ID() == ZERO()
+    assert ZERO() + ID() == ID()
+    assert ZERO() + ps == ps
+    assert ps + ZERO() == ps
+
+
+def test_from_list():
+    q = QubitPlaceholder.register(8)
+    terms_list = [("X", q[0]), ("Y", q[1]), ("Z", q[5])]
+    term = reduce(lambda x, y: x * y, [PauliTerm(*x) for x in terms_list])
+
+    pterm = PauliTerm.from_list(terms_list)
+    assert pterm == term
+
+    with pytest.raises(ValueError):
+        # terms are not on disjoint qubits
+        pterm = PauliTerm.from_list([("X", q[0]), ("Y", q[0])])
+
+
+def test_ordered():
+    q = QubitPlaceholder.register(8)
+    mapping = {x: i for i, x in enumerate(q)}
+    term = sZ(q[3]) * sZ(q[2]) * sZ(q[1])
+    prog = address_qubits(exponential_map(term)(0.5), mapping)
+    assert prog.out() == "CNOT 3 2\n" \
+                         "CNOT 2 1\n" \
+                         "RZ(1.0) 1\n" \
+                         "CNOT 2 1\n" \
+                         "CNOT 3 2\n"
+
+
+def test_simplify():
+    q = QubitPlaceholder.register(8)
+    t1 = sZ(q[0]) * sZ(q[1])
+    t2 = sZ(q[0]) * sZ(q[1])
+    assert (t1 + t2) == 2 * sZ(q[0]) * sZ(q[1])
+
+
+def test_dont_simplify():
+    q = QubitPlaceholder.register(8)
+    t1 = sZ(q[0]) * sZ(q[1])
+    t2 = sZ(q[2]) * sZ(q[3])
+    assert (t1 + t2) != 2 * sZ(q[0]) * sZ(q[1])
+
+
+def test_simplify_warning():
+    q = QubitPlaceholder.register(8)
+    t1 = sZ(q[0]) * sZ(q[1])
+    t2 = sZ(q[1]) * sZ(q[0])
+    with pytest.warns(UserWarning) as e:
+        tsum = t1 + t2
+
+    assert tsum == 2 * sZ(q[0]) * sZ(q[1])
+    assert 'will be combined with' in str(e[0].message)
diff --git a/pyquil/tests/test_quil.py b/pyquil/tests/test_quil.py
index 4edddaf..c610b48 100644
--- a/pyquil/tests/test_quil.py
+++ b/pyquil/tests/test_quil.py
@@ -14,7 +14,7 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 ##############################################################################
-
+import re
 from math import pi
 
 import numpy as np
@@ -24,8 +24,10 @@ from pyquil.gates import I, X, Y, Z, H, T, S, RX, RY, RZ, CNOT, CCNOT, PHASE, CP
     CPHASE10, CPHASE, SWAP, CSWAP, ISWAP, PSWAP, MEASURE, HALT, WAIT, NOP, RESET, \
     TRUE, FALSE, NOT, AND, OR, MOVE, EXCHANGE
 from pyquil.parameters import Parameter, quil_sin, quil_cos
-from pyquil.quil import Program, merge_programs, shift_quantum_gates, \
+from pyquil.paulis import exponential_map, sZ
+from pyquil.quil import Program, merge_programs, address_qubits, \
     get_classical_addresses_from_program
+from pyquil.quilatom import QubitPlaceholder
 from pyquil.quilbase import DefGate, Gate, Addr, Qubit, JumpWhen
 
 
@@ -134,7 +136,9 @@ def test_iteration():
     program = Program(gate_list)
     for ii, instruction in enumerate(program):
         assert instruction == gate_list[ii]
-    gate_generator = (gate_list[ii] for ii in range(3))   # https://github.com/rigetticomputing/pyquil/issues/265
+
+    # https://github.com/rigetticomputing/pyquil/issues/265
+    gate_generator = (gate_list[ii] for ii in range(3))
     program = Program(gate_generator)
     for ii, instruction in enumerate(program):
         assert instruction == gate_list[ii]
@@ -238,8 +242,8 @@ def test_dagger():
                        H(0), CNOT(0, 1), CCNOT(0, 1, 2),
                        SWAP(0, 1), CSWAP(0, 1, 2))
     assert p.dagger().out() == 'CSWAP 0 1 2\nSWAP 0 1\n' \
-        'CCNOT 0 1 2\nCNOT 0 1\nH 0\n' \
-        'Z 0\nY 0\nX 0\nI 0\n'
+                               'CCNOT 0 1 2\nCNOT 0 1\nH 0\n' \
+                               'Z 0\nY 0\nX 0\nI 0\n'
 
     # these gates require negating a parameter
     p = Program().inst(PHASE(pi, 0), RX(pi, 0), RY(pi, 0),
@@ -435,16 +439,98 @@ def test_alloc():
 
     p.inst(X(q3))  # X 4
 
+    with pytest.raises(RuntimeError) as e:
+        _ = p.out()
+    assert e.match(r'Qubit <QubitPlaceholder \d+> has not been assigned an index')
+
+
+def test_alloc_2():
+    p = Program()
+
+    p.inst(H(0))  # H 0
+
+    q1 = p.alloc()  # q1 = 1
+    q2 = p.alloc()  # q2 = 3
+
+    p.inst(CNOT(q1, q2))  # CNOT 1 3
+
+    p.inst(H(2))
+
+    q3 = p.alloc()  # q3 = 4
+
+    p.inst(X(q3))  # X 4
+    with pytest.raises(ValueError) as e:
+        _ = address_qubits(p, {
+            q1: 1,
+            q2: 3,
+            q3: 4,
+        })
+
+    assert e.match('Your program mixes instantiated qubits with placeholders')
+
+
+def test_alloc_new():
+    p = Program()
+
+    q0 = QubitPlaceholder()
+    p.inst(H(q0))  # H 0
+
+    q1 = QubitPlaceholder()
+    q2 = QubitPlaceholder()
+
+    p.inst(CNOT(q1, q2))  # CNOT 1 3
+
+    qxxx = QubitPlaceholder()
+    p.inst(H(qxxx))
+
+    q3 = QubitPlaceholder()
+
+    p.inst(X(q3))  # X 4
+    p = address_qubits(p, {
+        q1: 1,
+        q2: 3,
+        q3: 4,
+        q0: 0,
+        qxxx: 2,
+    })
+
     assert p.out() == "H 0\n" \
                       "CNOT 1 3\n" \
                       "H 2\n" \
                       "X 4\n"
 
 
+def test_multiaddress():
+    p = Program()
+    q0, q1 = [QubitPlaceholder() for _ in range(2)]
+    p += exponential_map(sZ(q0) * sZ(q1))(0.5)
+
+    map1 = {q0: 0, q1: 1}
+    map2 = {q0: 9, q1: 10}
+
+    p1 = address_qubits(p, map1)
+
+    with pytest.raises(RuntimeError):
+        _ = p.out()  # make sure the original isn't affected
+
+    assert p1.out() == "CNOT 0 1\n" \
+                       "RZ(1.0) 1\n" \
+                       "CNOT 0 1\n"
+
+    p2 = address_qubits(p, map2)
+    assert p1.out() == "CNOT 0 1\n" \
+                       "RZ(1.0) 1\n" \
+                       "CNOT 0 1\n"
+    assert p2.out() == "CNOT 9 10\n" \
+                       "RZ(1.0) 10\n" \
+                       "CNOT 9 10\n"
+
+
 def test_multiple_instantiate():
     p = Program()
     q = p.alloc()
     p.inst(H(q))
+    p = address_qubits(p)
     assert p.out() == 'H 0\n'
     assert p.out() == 'H 0\n'
 
@@ -456,6 +542,7 @@ def test_reuse_alloc():
     p.inst(H(q1))
     p.inst(H(q2))
     p.inst(CNOT(q1, q2))
+    p = address_qubits(p)
     assert p.out() == 'H 0\nH 1\nCNOT 0 1\n'
 
 
@@ -465,18 +552,15 @@ def test_prog_merge():
     assert merge_programs([prog_0, prog_1]).out() == (prog_0 + prog_1).out()
 
 
-def test_quantum_gate_shift():
-    prog = Program(X(0), CNOT(0, 4), MEASURE(5, [5]))
-    assert shift_quantum_gates(prog, 5) == Program(X(5), CNOT(5, 9), MEASURE(5, [5]))
-
-
 def test_get_qubits():
     pq = Program(X(0), CNOT(0, 4), MEASURE(5, [5]))
     assert pq.get_qubits() == {0, 4, 5}
 
+    q = [QubitPlaceholder() for _ in range(6)]
+    pq = Program(X(q[0]), CNOT(q[0], q[4]), MEASURE(q[5], [5]))
     qq = pq.alloc()
-    pq.inst(Y(2), X(qq))
-    assert pq.get_qubits() == {0, 1, 2, 4, 5}  # this synthesizes the allocation
+    pq.inst(Y(q[2]), X(qq))
+    assert address_qubits(pq).get_qubits() == {0, 1, 2, 3, 4}
 
     qubit_index = 1
     p = Program(("H", qubit_index))
@@ -484,7 +568,20 @@ def test_get_qubits():
     q1 = p.alloc()
     q2 = p.alloc()
     p.inst(("CNOT", q1, q2))
-    assert p.get_qubits() == {qubit_index, 0, 2}
+    with pytest.raises(ValueError) as e:
+        _ = address_qubits(p).get_qubits()
+    assert e.match('Your program mixes instantiated qubits with placeholders')
+
+
+def test_get_qubit_placeholders():
+    qs = QubitPlaceholder.register(8)
+    pq = Program(X(qs[0]), CNOT(qs[0], qs[4]), MEASURE(qs[5], [5]))
+    assert pq.get_qubits() == {qs[i] for i in [0, 4, 5]}
+
+
+def test_get_qubits_not_as_indices():
+    pq = Program(X(0), CNOT(0, 4), MEASURE(5, [5]))
+    assert pq.get_qubits(indices=False) == {Qubit(i) for i in [0, 4, 5]}
 
 
 def test_eq():
@@ -492,6 +589,7 @@ def test_eq():
     q1 = p1.alloc()
     q2 = p1.alloc()
     p1.inst([H(q1), CNOT(q1, q2)])
+    p1 = address_qubits(p1)
 
     p2 = Program()
     p2.inst([H(0), CNOT(0, 1)])
@@ -630,7 +728,7 @@ def test_allocating_qubits_on_multiple_programs():
     qubit1 = q.alloc()
     q.inst(X(qubit1))
 
-    assert (p + q).out() == "X 0\nX 1\n"
+    assert address_qubits(p + q).out() == "X 0\nX 1\n"
 
 
 # https://github.com/rigetticomputing/pyquil/issues/163
@@ -652,7 +750,7 @@ def test_nesting_a_program_inside_itself():
 def test_inline_alloc():
     p = Program()
     p += H(p.alloc())
-    assert p.out() == "H 0\n"
+    assert address_qubits(p).out() == "H 0\n"
 
 
 # https://github.com/rigetticomputing/pyquil/issues/138
@@ -662,6 +760,19 @@ def test_defgate_integer_input():
     assert dg.out() == "DEFGATE TEST:\n    1, 0\n    0, 1\n"
 
 
+def test_out_vs_str():
+    qs = QubitPlaceholder.register(6)
+    pq = Program(X(qs[0]), CNOT(qs[0], qs[4]), MEASURE(qs[5], [5]))
+
+    with pytest.raises(RuntimeError) as e:
+        pq.out()
+    assert e.match(r'Qubit <.*> has not been assigned an index')
+
+    string_version = str(pq)
+    should_be_re = (r'X <.*>\nCNOT <.*> <.*>\nMEASURE <.*> \[5\]\n')
+    assert re.fullmatch(should_be_re, string_version, flags=re.MULTILINE)
+
+
 def test_get_classical_addresses_from_program():
     p = Program([H(i) for i in range(4)])
     assert get_classical_addresses_from_program(p) == []