Adds ability to return multi-qubit expectation values to all devi… (#13)

* added ability to return multi-qubit expectation values to all devices

* cleaned up by adding separate probability method

* typo

* added multimode hermitian expval test to pyqvm

* added suggested changes

* undo conftest change

* update to work with new pennylane

* correct _expval_queue -> -obs_queue

* remove deprecated pytest message kwarg

* swap out U2 matrix in tests

* fix two incorrect exception messages in tests

* fix typo postive->positive

pennylane_forest/_version.py

@@ -2,4 +2,4 @@
    Version number (major.minor.patch[-label])
 """
 
-__version__ = '0.2.0'
+__version__ = '0.3.0'

pennylane_forest/device.py

@@ -164,7 +164,7 @@ class ForestDevice(Device):
     Args:
         wires (int): the number of modes to initialize the device in
         shots (int): Number of circuit evaluations/random samples used
-            to estimate expectation values of expectations.
+            to estimate expectation values of observables.
             For simulator devices, 0 means the exact EV is returned.
 
     Keyword args:
@@ -234,7 +234,7 @@ def apply(self, operation, wires, par):
             self.active_wires = set(range(self.num_wires))
 
     @abc.abstractmethod
-    def pre_expval(self): #pragma no cover
+    def pre_measure(self): #pragma no cover
         """Run the QVM or QPU"""
         raise NotImplementedError
 

pennylane_forest/numpy_wavefunction.py

@@ -19,13 +19,15 @@
 Code details
 ~~~~~~~~~~~~
 """
+import itertools
+
 import numpy as np
 
 from pyquil.pyqvm import PyQVM
 from pyquil.numpy_simulator import NumpyWavefunctionSimulator
 
 from .device import ForestDevice
-from .wavefunction import expectation_map, spectral_decomposition_qubit
+from .wavefunction import observable_map, spectral_decomposition_qubit
 from ._version import __version__
 
 
@@ -35,12 +37,12 @@ class NumpyWavefunctionDevice(ForestDevice):
     Args:
         wires (int): the number of qubits to initialize the device in
         shots (int): Number of circuit evaluations/random samples used
-            to estimate expectation values of expectations.
+            to estimate expectation values of observables.
     """
     name = 'pyQVM NumpyWavefunction Simulator Device'
     short_name = 'forest.numpy_wavefunction'
 
-    expectations = {'PauliX', 'PauliY', 'PauliZ', 'Hadamard', 'Hermitian', 'Identity'}
+    observables = {'PauliX', 'PauliY', 'PauliZ', 'Hadamard', 'Hermitian', 'Identity'}
 
     def __init__(self, wires, *, shots=0, **kwargs):
         super().__init__(wires, shots, **kwargs)
@@ -51,20 +53,19 @@ def pre_apply(self):
         self.reset()
         self.qc.wf_simulator.reset()
 
-    def pre_expval(self):
+    def pre_measure(self):
         # TODO: currently, the PyQVM considers qubit 0 as the leftmost bit and therefore
         # returns amplitudes in the opposite of the Rigetti Lisp QVM (which considers qubit
         # 0 as the rightmost bit). This may change in the future, so in the future this
-        # might need to get udpated to be similar to the pre_expval function of
+        # might need to get udpated to be similar to the pre_measure function of
         #pennylane_forest/wavefunction.py
         self.state = self.qc.execute(self.prog).wf_simulator.wf.flatten()
 
-    def expval(self, expectation, wires, par):
-        # measurement/expectation value <psi|A|psi>
-        if expectation == 'Hermitian':
+    def expval(self, observable, wires, par):
+        if observable == 'Hermitian':
             A = par[0]
         else:
-            A = expectation_map[expectation]
+            A = observable_map[observable]
 
         if self.shots == 0:
             # exact expectation value
@@ -90,25 +91,50 @@ def ev(self, A, wires):
           float: expectation value :math:`\left\langle{A}\right\rangle = \left\langle{\psi}\mid A\mid{\psi}\right\rangle`
         """
         # Expand the Hermitian observable over the entire subsystem
-        A = self.expand_one(A, wires)
+        A = self.expand(A, wires)
         return np.vdot(self.state, A @ self.state).real
 
-    def expand_one(self, U, wires):
-        r"""Expand a one-qubit operator into a full system operator.
+    def expand(self, U, wires):
+        r"""Expand a multi-qubit operator into a full system operator.
 
         Args:
-          U (array): :math:`2\times 2` matrix
-          wires (Sequence[int]): target subsystem
+            U (array): :math:`2^n \times 2^n` matrix where n = len(wires).
+            wires (Sequence[int]): Target subsystems (order matters! the
+                left-most Hilbert space is at index 0).
 
         Returns:
-          array: :math:`2^n\times 2^n` matrix
+            array: :math:`2^N\times 2^N` matrix. The full system operator.
         """
-        if U.shape != (2, 2):
-            raise ValueError('2x2 matrix required.')
-        if len(wires) != 1:
-            raise ValueError('One target subsystem required.')
-        wires = wires[0]
-        before = 2**wires
-        after = 2**(self.num_wires-wires-1)
-        U = np.kron(np.kron(np.eye(before), U), np.eye(after))
-        return U
+        if self.num_wires == 1:
+            # total number of wires is 1, simply return the matrix
+            return U
+
+        N = self.num_wires
+        wires = np.asarray(wires)
+
+        if np.any(wires < 0) or np.any(wires >= N) or len(set(wires)) != len(wires):
+            raise ValueError("Invalid target subsystems provided in 'wires' argument.")
+
+        if U.shape != (2**len(wires), 2**len(wires)):
+            raise ValueError("Matrix parameter must be of size (2**len(wires), 2**len(wires))")
+
+        # generate N qubit basis states via the cartesian product
+        tuples = np.array(list(itertools.product([0, 1], repeat=N)))
+
+        # wires not acted on by the operator
+        inactive_wires = list(set(range(N))-set(wires))
+
+        # expand U to act on the entire system
+        U = np.kron(U, np.identity(2**len(inactive_wires)))
+
+        # move active wires to beginning of the list of wires
+        rearranged_wires = np.array(list(wires)+inactive_wires)
+
+        # convert to computational basis
+        # i.e., converting the list of basis state bit strings into
+        # a list of decimal numbers that correspond to the computational
+        # basis state. For example, [0, 1, 0, 1, 1] = 2^3+2^1+2^0 = 11.
+        perm = np.ravel_multi_index(tuples[:, rearranged_wires].T, [2]*N)
+
+        # permute U to take into account rearranged wires
+        return U[:, perm][perm]

pennylane_forest/qpu.py

@@ -33,7 +33,7 @@ class QPUDevice(QVMDevice):
     Args:
         device (str): the name of the device to initialise.
         shots (int): number of circuit evaluations/random samples used
-            to estimate expectation values of expectations.
+            to estimate expectation values of observables.
         active_reset (bool): whether to actively reset qubits instead of waiting for
             for qubits to decay to the ground state naturally.
             Setting this to ``True`` results in a significantly faster expectation value

pennylane_forest/qvm.py

@@ -50,7 +50,7 @@ class QVMDevice(ForestDevice):
             * Graph topology representing the device architecture.
 
         shots (int): number of circuit evaluations/random samples used
-            to estimate expectation values of expectations.
+            to estimate expectation values of observables.
         noisy (bool): set to ``True`` to add noise models to your QVM.
 
     Keyword args:
@@ -66,7 +66,7 @@ class QVMDevice(ForestDevice):
     """
     name = 'Forest QVM Device'
     short_name = 'forest.qvm'
-    expectations = {'PauliX', 'PauliY', 'PauliZ', 'Identity', 'Hadamard', 'Hermitian'}
+    observables = {'PauliX', 'PauliY', 'PauliZ', 'Identity', 'Hadamard', 'Hermitian'}
 
     def __init__(self, device, *, shots=1024, noisy=False, **kwargs):
         self._eigs = {}
@@ -106,25 +106,25 @@ def __init__(self, device, *, shots=1024, noisy=False, **kwargs):
 
         self.active_reset = False
 
-    def pre_expval(self):
+    def pre_measure(self):
         """Run the QVM"""
         # pylint: disable=attribute-defined-outside-init
-        for e in self.expval_queue:
-            wire = [e.wires[0]]
+        for e in self.obs_queue:
+            wires = e.wires
 
             if e.name == 'PauliX':
                 # X = H.Z.H
-                self.apply('Hadamard', wire, [])
+                self.apply('Hadamard', wires, [])
 
             elif e.name == 'PauliY':
                 # Y = (HS^)^.Z.(HS^) and S^=SZ
-                self.apply('PauliZ', wire, [])
-                self.apply('S', wire, [])
-                self.apply('Hadamard', wire, [])
+                self.apply('PauliZ', wires, [])
+                self.apply('S', wires, [])
+                self.apply('Hadamard', wires, [])
 
             elif e.name == 'Hadamard':
                 # H = Ry(-pi/4)^.Z.Ry(-pi/4)
-                self.apply('RY', wire, [-np.pi/4])
+                self.apply('RY', wires, [-np.pi/4])
 
             elif e.name == 'Hermitian':
                 # For arbitrary Hermitian matrix H, let U be the unitary matrix
@@ -143,7 +143,7 @@ def pre_expval(self):
                     self._eigs[Hkey] = {'eigval': w, 'eigvec': U}
 
                 # Perform a change of basis before measuring by applying U^ to the circuit
-                self.apply('QubitUnitary', wire, [U.conj().T])
+                self.apply('QubitUnitary', wires, [U.conj().T])
 
         prag = Program(Pragma('INITIAL_REWIRING', ['"PARTIAL"']))
 
@@ -169,22 +169,75 @@ def pre_expval(self):
         for i, q in enumerate(qubits):
             self.state[q] = bitstring_array[:, i]
 
-    def expval(self, expectation, wires, par):
-        evZ = np.mean(1-2*self.state[wires[0]])
+    def expval(self, observable, wires, par):
+        if len(wires) == 1:
+            # 1 qubit observable
+            evZ = np.mean(1-2*self.state[wires[0]])
 
-        # for single qubit state probabilities |psi|^2 = (p0, p1),
-        # we know that p0+p1=1 and that <Z>=p0-p1
-        p0 = (1+evZ)/2
-        p1 = (1-evZ)/2
+            # for single qubit state probabilities |psi|^2 = (p0, p1),
+            # we know that p0+p1=1 and that <Z>=p0-p1
+            p0 = (1+evZ)/2
+            p1 = (1-evZ)/2
 
-        if expectation == 'Identity':
-            # <I> = \sum_i p_i
-            return p0+p1
+            if observable == 'Identity':
+                # <I> = \sum_i p_i
+                return p0 + p1
 
-        if expectation == 'Hermitian':
-            # <H> = \sum_i w_i p_i
+            if observable == 'Hermitian':
+                # <H> = \sum_i w_i p_i
+                Hkey = tuple(par[0].flatten().tolist())
+                w = self._eigs[Hkey]['eigval']
+                return w[0]*p0 + w[1]*p1
+
+            return evZ
+
+        # Multi-qubit observable
+        # ----------------------
+        # Currently, we only support qml.expval.Hermitian(A, wires),
+        # where A is a 2^N x 2^N matrix acting on N wires.
+        #
+        # Eventually, we will also support tensor products of Pauli
+        # matrices in the PennyLane UI.
+
+        probs = self.probabilities(wires)
+
+        if observable == 'Hermitian':
             Hkey = tuple(par[0].flatten().tolist())
             w = self._eigs[Hkey]['eigval']
-            return w[0]*p0 + w[1]*p1
-
-        return evZ
+            # <A> = \sum_i w_i p_i
+            return w @ probs
+
+    def probabilities(self, wires):
+        """Returns the (marginal) probabilities of the quantum state.
+
+        Args:
+            wires (Sequence[int]): sequence of wires to return
+                marginal probabilities for. Wires not provided
+                are traced out of the system.
+
+        Returns:
+            array: array of shape ``[2**len(wires)]`` containing
+            the probabilities of each computational basis state
+        """
+
+        # create an array of size [2^len(wires), 2] to store
+        # the resulting probability of each computational basis state
+        probs = np.zeros([2**len(wires), 2])
+        probs[:, 0] = np.arange(2**len(wires))
+
+        # extract the measured samples
+        res = np.array([self.state[w] for w in wires]).T
+        for i in res:
+            # for each sample, calculate which
+            # computational basis state it corresponds to
+            cb = np.sum(2**np.arange(len(wires)-1, -1, -1)*i)
+            # add a tally for this computational basis state
+            # to our array of basis probabilities
+            probs[cb, 1] += 1
+
+        # sort the probabilities by the first column,
+        # and divide by the number of shots
+        probs = probs[probs[:, 0].argsort()] / self.shots
+        probs = probs[:, 1]
+
+        return probs