Add support for samples (#18)

* added sampling support

* added tests

* suggested changes

* added wires options

* added wires options

* added wires options

* added wires options

* added wires options

* undo

pennylane_forest/_version.py

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

pennylane_forest/numpy_wavefunction.py

@@ -26,12 +26,11 @@
 from pyquil.pyqvm import PyQVM
 from pyquil.numpy_simulator import NumpyWavefunctionSimulator
 
-from .device import ForestDevice
-from .wavefunction import observable_map, spectral_decomposition_qubit
+from .wavefunction import WavefunctionDevice
 from ._version import __version__
 
 
-class NumpyWavefunctionDevice(ForestDevice):
+class NumpyWavefunctionDevice(WavefunctionDevice):
     r"""NumpyWavefunction simulator device for PennyLane.
 
     Args:
@@ -45,7 +44,7 @@ class NumpyWavefunctionDevice(ForestDevice):
     observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity"}
 
     def __init__(self, wires, *, shots=0, **kwargs):
-        super().__init__(wires, shots, **kwargs)
+        super(WavefunctionDevice, self).__init__(wires, shots, **kwargs)
         self.qc = PyQVM(n_qubits=wires, quantum_simulator_type=NumpyWavefunctionSimulator)
         self.state = None
 
@@ -60,45 +59,3 @@ def pre_measure(self):
         # 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, observable, wires, par):
-        if observable == "Hermitian":
-            A = par[0]
-        else:
-            A = observable_map[observable]
-
-        if self.shots == 0:
-            # exact expectation value
-            ev = self.ev(A, wires)
-        else:
-            # estimate the ev
-            # sample Bernoulli distribution n_eval times / binomial distribution once
-            a, P = spectral_decomposition_qubit(A)
-            p0 = self.ev(P[0], wires)  # probability of measuring a[0]
-            n0 = np.random.binomial(self.shots, p0)
-            ev = (n0 * a[0] + (self.shots - n0) * a[1]) / self.shots
-
-        return ev
-
-    def var(self, observable, wires, par):
-        if observable == "Hermitian":
-            A = par[0]
-        else:
-            A = observable_map[observable]
-
-        var = self.ev(A @ A, wires) - self.ev(A, wires) ** 2
-        return var
-
-    def ev(self, A, wires):
-        r"""Evaluates a one-qubit expectation in the current state.
-
-        Args:
-          A (array): :math:`2\times 2` Hermitian matrix corresponding to the expectation
-          wires (Sequence[int]): target subsystem
-
-        Returns:
-          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
-        As = self.mat_vec_product(A, self.state, wires)
-        return np.vdot(self.state, As).real

pennylane_forest/qvm.py

@@ -19,7 +19,9 @@
 Code details
 ~~~~~~~~~~~~
 """
+import itertools
 import re
+
 import numpy as np
 
 import networkx as nx
@@ -174,42 +176,37 @@ def pre_measure(self):
             self.state[q] = bitstring_array[:, i]
 
     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
+        return np.mean(self.sample(observable, wires, par))
 
-            if observable == "Identity":
-                # <I> = \sum_i p_i
-                return p0 + p1
+    def var(self, observable, wires, par):
+        return np.var(self.sample(observable, wires, par))
 
-            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
+    def sample(self, observable, wires, par, n=None):
+        if n is None:
+            n = self.shots
 
-            return evZ
+        if n == 0:
+            raise ValueError("Calling sample with n = 0 is not possible.")
 
-        # 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.
+        if n < 0 or not isinstance(n, int):
+            raise ValueError("The number of samples must be a positive integer.")
 
-        probs = self.probabilities(wires)
+        if observable == "Identity":
+            return np.ones([n])
 
         if observable == "Hermitian":
             Hkey = tuple(par[0].flatten().tolist())
-            w = self._eigs[Hkey]["eigval"]
-            # <A> = \sum_i w_i p_i
-            return w @ probs
+            eigvals = self._eigs[Hkey]["eigval"]
+            res = np.array([self.state[i] for i in wires]).T
+
+            samples = np.zeros([n])
+
+            for w, b in zip(eigvals, itertools.product([0, 1], repeat=len(wires))):
+                samples = np.where(np.all(res == b, axis=1), w, samples)
+
+            return samples
+
+        return 1 - 2 * self.state[wires[0]]
 
     def probabilities(self, wires):
         """Returns the (marginal) probabilities of the quantum state.
@@ -223,7 +220,6 @@ def probabilities(self, wires):
             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])
@@ -245,41 +241,3 @@ def probabilities(self, wires):
         probs = probs[:, 1]
 
         return probs
-
-    def var(self, observable, wires, par):
-        if len(wires) == 1:
-            # 1 qubit observable
-            if observable == "Identity":
-                return 0
-
-            varZ = np.var(1 - 2 * self.state[wires[0]])
-
-            if observable in {"PauliX", "PauliY", "PauliZ", "Hadamard"}:
-                return varZ
-
-            # for single qubit state probabilities |psi|^2 = (p0, p1),
-            # we know that p0+p1=1 and that <Z>=p0-p1
-            evZ = np.mean(1 - 2 * self.state[wires[0]])
-            probs = np.array([(1 + evZ) / 2, (1 - evZ) / 2])
-
-            if observable == "Hermitian":
-                # <H> = \sum_i w_i p_i
-                Hkey = tuple(par[0].flatten().tolist())
-                w = self._eigs[Hkey]["eigval"]
-                return (w ** 2) @ probs - (w @ probs) ** 2
-
-        # 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"]
-            # <A> = \sum_i w_i p_i
-            return (w ** 2) @ probs - (w @ probs) ** 2

pennylane_forest/wavefunction.py

@@ -29,6 +29,7 @@
 import itertools
 
 import numpy as np
+from numpy.linalg import eigh
 
 from pyquil.api import WavefunctionSimulator
 
@@ -46,24 +47,6 @@
 observable_map = {"PauliX": X, "PauliY": Y, "PauliZ": Z, "Identity": I, "Hadamard": H}
 
 
-def spectral_decomposition_qubit(A):
-    r"""Spectral decomposition of a :math:`2\times 2` Hermitian matrix.
-
-    Args:
-        A (array): :math:`2\times 2` Hermitian matrix
-
-    Returns:
-        (vector[float], list[array[complex]]): (a, P): eigenvalues and hermitian projectors
-        such that :math:`A = \sum_k a_k P_k`.
-    """
-    d, v = np.linalg.eigh(A)
-    P = []
-    for k in range(2):
-        temp = v[:, k]
-        P.append(np.outer(temp, temp.conj()))
-    return d, P
-
-
 class WavefunctionDevice(ForestDevice):
     r"""Wavefunction simulator device for PennyLane.
 
@@ -136,11 +119,7 @@ def expval(self, observable, wires, par):
             ev = self.ev(A, wires)
         else:
             # estimate the ev
-            # sample Bernoulli distribution n_eval times / binomial distribution once
-            a, P = spectral_decomposition_qubit(A)
-            p0 = self.ev(P[0], wires)  # probability of measuring a[0]
-            n0 = np.random.binomial(self.shots, p0)
-            ev = (n0 * a[0] + (self.shots - n0) * a[1]) / self.shots
+            ev = np.mean(self.sample(observable, wires, par, self.shots))
 
         return ev
 
@@ -150,14 +129,48 @@ def var(self, observable, wires, par):
         else:
             A = observable_map[observable]
 
-        var = self.ev(A @ A, wires) - self.ev(A, wires) ** 2
+        if self.shots == 0:
+            # exact variance value
+            var = self.ev(A @ A, wires) - self.ev(A, wires)**2
+        else:
+            # estimate the variance
+            var = np.var(self.sample(observable, wires, par, self.shots))
+
         return var
 
+    def sample(self, observable, wires, par, n=None):
+        if n is None:
+            n = self.shots
+
+        if n == 0:
+            raise ValueError("Calling sample with n = 0 is not possible.")
+
+        if n < 0 or not isinstance(n, int):
+            raise ValueError("The number of samples must be a positive integer.")
+
+        if observable == "Hermitian":
+            A = par[0]
+        else:
+            A = observable_map[observable]
+
+        # Calculate the probability p of observing
+        # eigenvalue a for the observable A.
+        a, P = self.spectral_decomposition(A)
+
+        p = np.zeros(a.shape)
+
+        for idx, Pi in enumerate(P):
+            p[idx] = self.ev(Pi, wires)
+
+        # randomly sample from the eigenvalues of A
+        # according to the computed probabilities
+        return np.random.choice(a, n, p=p)
+
     def ev(self, A, wires):
-        r"""Evaluates a one-qubit expectation in the current state.
+        r"""Evaluates an expectation in the current state.
 
         Args:
-          A (array): :math:`2\times 2` Hermitian matrix corresponding to the expectation
+          A (array): :math:`2^N\times 2^N` Hermitian matrix corresponding to the expectation
           wires (Sequence[int]): target subsystem
 
         Returns:
@@ -166,3 +179,21 @@ def ev(self, A, wires):
         # Expand the Hermitian observable over the entire subsystem
         As = self.mat_vec_product(A, self.state, wires)
         return np.vdot(self.state, As).real
+
+    @staticmethod
+    def spectral_decomposition(A):
+        r"""Spectral decomposition of a Hermitian matrix.
+
+        Args:
+            A (array): Hermitian matrix
+
+        Returns:
+            (vector[float], list[array[complex]]): (a, P): eigenvalues and Hermitian projectors
+                such that :math:`A = \sum_k a_k P_k`.
+        """
+        d, v = eigh(A)
+        P = []
+        for k in range(d.shape[0]):
+            temp = v[:, k]
+            P.append(np.outer(temp, temp.conj()))
+        return d, P