Add Wigner Function Functionality and Documentation.

README.md

@@ -145,7 +145,7 @@ Visit the [IBM Q experience community](https://quantumexperience.ng.bluemix.net/
 
 ## Authors (alphabetical)
 
-Jim Challenger, Andrew Cross, Ismael Faro, Jay Gambetta, Juan Gomez, Paco Martin, Antonio Mezzacapo, Jesus Perez, and John Smolin, Erick Winston, Chris Wood.
+Jim Challenger, Andrew Cross, Vincent Dwyer, Mark Everitt, Ismael Faro, Jay Gambetta, Juan Gomez, Paco Martin, Antonio Mezzacapo, Jesus Perez, Russell Rundle, Todd Tilma, John Smolin, Erick Winston, Chris Wood
 
 In future releases, anyone who contributes with code to this project is welcome to include their name here.
 

qiskit/tools/qcvv/tomography.py

@@ -711,3 +711,99 @@ def fit_tomography_data(data, method=None, options=None):
         return rho
     else:
         print('error: method unknown reconstruction method "%s"' % method)
+
+
+###############################################################
+# Wigner function tomography
+###############################################################
+
+def build_wigner_circuits(Q_program, name, phis, thetas, qubits,
+                          qreg, creg, silent=False):
+
+    """Create the circuits to rotate to points in phase space
+
+    Args:
+        Q_program (QuantumProgram): A quantum program to store the circuits.
+        name (string): The name of the base circuit to be appended.
+        phis (np.matrix[[complex]]):
+        thetas (np.matrix[[complex]]):
+        qubits (list[int]): a list of the qubit indexes of qreg to be measured.
+        qreg (QuantumRegister): the quantum register containing qubits to be
+                                measured.
+        creg (ClassicalRegister): the classical register containing bits to
+                                    store measurement outcomes.
+        silent (bool, optional): hide verbose output.
+
+    Returns: A list of names of the added wigner function circuits.
+
+    """
+
+    orig = Q_program.get_circuit(name)
+    labels = []
+    points = len(phis[0])
+
+
+    for point in range(points):
+        label = '_wigner_phase_point'
+        label += str(point)
+        circuit = Q_program.create_circuit(label, [qreg], [creg])
+        c_index = 0
+
+        for qubit in range(len(qubits)):
+            circuit.u3(thetas[qubit][point], 0,
+                       phis[qubit][point],qreg[qubits[qubit]])
+            circuit.measure(qreg[qubits[qubit]],creg[qubits[qubit]])
+
+        Q_program.add_circuit(name+label, orig+circuit)
+        labels.append(name+label)
+
+    if not silent:
+        print('>> created Wigner function circuits for "%s"' % name)
+    return labels
+
+
+def wigner_data(Q_result, name, meas_qubits, labels, shots=None):
+
+    """Get the value of the Wigner function from measurement results.
+
+    Args:
+        Q_result (Result): Results from execution of a state tomography
+                            circuits on a backend.
+        name (string): The name of the base state preparation circuit.
+        meas_qubits (list[int]): a list of the qubit indexes measured.
+        labels : a list of names of the circuits
+
+    Returns: The values of the Wigner function at measured points in
+            phase space
+
+    """
+
+    num = len(meas_qubits)
+
+    dim = 2**num
+    P = [0.5+0.5*np.sqrt(3),0.5-0.5*np.sqrt(3)]
+    parity = 1
+
+    for i in range(num):
+        parity = np.kron(parity,P)
+
+    W = [0]*len(labels)
+    wpt = 0
+    counts = [marginal_counts(Q_result.get_counts(circ), meas_qubits)
+              for circ in labels]
+    for entry in counts:
+        x =[0]*dim
+
+        for i in range(dim):
+            if bin(i)[2:].zfill(num) in entry:
+                x[i] = float(entry[bin(i)[2:].zfill(num)])
+
+            if shots is None:
+                shots = np.sum(x)
+
+                for i in range(dim):
+                W[wpt] = W[wpt]+(x[i]/shots)*parity[i]
+
+            wpt += 1
+
+    return W

qiskit/tools/visualization.py

@@ -23,7 +23,9 @@
 from scipy import linalg as la
 from collections import Counter
 import matplotlib.pyplot as plt
+from matplotlib import cm
 from mpl_toolkits.mplot3d import proj3d
+from mpl_toolkits.mplot3d import Axes3D
 from matplotlib.patches import FancyArrowPatch
 from qiskit.tools.qi.pauli import pauli_group, pauli_singles
 
@@ -447,5 +449,241 @@ def plot_state(rho, method='city'):
                                np.dot(x.to_matrix(), rho))),
                                pauli_singles(i, num)))
             plot_bloch_vector(bloch_state, "qubit " + str(i))
+    elif method == "wigner":
+        plot_wigner_function(rho)
+    else:
+        print("No method given")
+
+###############################################################
+# Plotting Wigner functions
+###############################################################
+
+def plot_wigner_function(state, res=100):
+
+    """Plot the equal angle slice spin Wigner function of an arbitrary
+    quantum state
+
+    Args:
+        state (np.matrix[[complex]]): - Matrix of 2**n x 2**n complex
+                                        numbers
+                                    - State Vector of 2**n x 1 complex
+                                        numbers
+    res (int) : number of theta and phi values in meshgrid
+                on sphere (creates a res x res grid of points)
+    Returns:
+        none: plot is shown with matplotlib on the screen
+    References:
+        [1] T. Tilma, M. J. Everitt, J. H. Samson, W. J. Munro,
+        and K. Nemoto, Phys. Rev. Lett. 117, 180401 (2016).
+        [2] R. P. Rundle, P. W. Mills, T. Tilma, J. H. Samson, and
+        M. J. Everitt, Phys. Rev. A 96, 022117 (2017).
+
+    """"
+    state = np.array(state)
+    if state.ndim == 1:
+        state = np.outer(state,state) # turns state vector to a density matrix
+    state = np.matrix(state)
+    num = int(np.log2(len(state))) # number of qubits
+    phi_vals = np.linspace(0,np.pi,num=res,
+                           dtype = np.complex_)
+    theta_vals = np.linspace(0,0.5*np.pi,num=res,
+                             dtype = np.complex_) # phi and theta values for WF
+    W = np.empty([res,res])
+    harr = np.sqrt(3)
+    Delta_su2 = np.zeros((2,2),dtype = np.complex_)
+    #create the spin Wigner function
+
+    for theta in range(res):
+        costheta  = harr*np.cos(2*theta_vals[theta])
+        sintheta  = harr*np.sin(2*theta_vals[theta])
+
+        for phi in range(res):
+            Delta_su2[0,0] =  0.5*(1+costheta)
+            Delta_su2[0,1] = -0.5*(np.exp(2j*phi_vals[phi])*sintheta)
+            Delta_su2[1,0] = -0.5*(np.exp(-2j*phi_vals[phi])*sintheta)
+            Delta_su2[1,1] =  0.5*(1-costheta)
+            kernel = 1
+            for i in range(num):
+                kernel = np.kron(kernel,Delta_su2) # creates phase point kernel              
+
+            W[phi,theta] = np.real(np.trace(state*kernel)) # The Wigner function
+
+    # Plot a sphere (x,y,z) with Wigner function facecolor data stored in Wc
+    fig = plt.figure(figsize=(11,9))
+    ax = fig.gca(projection = '3d')
+    Wmax = np.amax(W)
+    #color data for plotting
+    Wc = cm.seismic_r((W+Wmax)/(2*Wmax)) # color data for sphere
+    Wc2 = cm.seismic_r((W[0:res, int(res/2):res]+Wmax)/(2*Wmax)) # bottom
+    Wc3 = cm.seismic_r((W[int(res/4):int(3*res/4), 0:res]+Wmax)/(2*Wmax)) #side
+    Wc4 = cm.seismic_r((W[int(res/2):res, 0:res]+Wmax)/(2*Wmax)) #  back
+
+    u = np.linspace(0, 2 * np.pi, res)
+    v = np.linspace(0, np.pi, res)
+    x = np.outer(np.cos(u), np.sin(v))
+    y = np.outer(np.sin(u), np.sin(v))
+    z = np.outer(np.ones(np.size(u)), np.cos(v)) # creates a sphere mesh
+
+    ax.plot_surface(x,y,z, facecolors=Wc, 
+                    vmin=-Wmax, vmax=Wmax,
+                    rcount=res, ccount=res,
+                    linewidth=0, zorder=0.5,
+                    antialiased=False) # plots Wigner Bloch sphere
+
+    ax.plot_surface(x[0:res, int(res/2):res],
+                    y[0:res, int(res/2):res],
+                    -1.5*np.ones((res,int(res/2))),
+                    facecolors=Wc2,
+                    vmin=-Wmax, vmax=Wmax,
+                    rcount=res/2, ccount=res/2,
+                    linewidth=0, zorder=0.5,
+                    antialiased=False) # plots bottom reflection
+
+    ax.plot_surface(-1.5*np.ones((int(res/2), res)),
+                    y[int(res/4):int(3*res/4), 0:res],
+                    z[int(res/4):int(3*res/4), 0:res],
+                    facecolors=Wc3,
+                    vmin=-Wmax, vmax=Wmax,
+                    rcount=res/2, ccount=res/2,
+                    linewidth=0, zorder=0.5,
+                    antialiased=False) # plots side reflection
+
+    ax.plot_surface(x[int(res/2):res, 0:res],
+                    1.5*np.ones((int(res/2), res)),
+                    z[int(res/2):res, 0:res],
+                    facecolors=Wc4,
+                    vmin=-Wmax, vmax=Wmax,
+                    rcount=res/2, ccount=res/2,
+                    linewidth=0, zorder=0.5,
+                    antialiased=False) # plots back reflection
+
+    ax.w_xaxis.set_pane_color((0.4, 0.4, 0.4, 1.0))
+    ax.w_yaxis.set_pane_color((0.4, 0.4, 0.4, 1.0))
+    ax.w_zaxis.set_pane_color((0.4, 0.4, 0.4, 1.0))
+    ax.set_xticks([], [])
+    ax.set_yticks([], [])
+    ax.set_zticks([], [])
+    ax.grid(False)
+    ax.xaxis.pane.set_edgecolor('black')
+    ax.yaxis.pane.set_edgecolor('black')
+    ax.zaxis.pane.set_edgecolor('black')
+    ax.set_xlim(-1.5, 1.5)
+    ax.set_ylim(-1.5, 1.5)
+    ax.set_zlim(-1.5, 1.5)
+    m = cm.ScalarMappable(cmap=cm.seismic_r)
+    m.set_array([-Wmax, Wmax])
+    plt.colorbar(m, shrink=0.5, aspect=10)
+
+    plt.show()
+
+def plot_wigner_curve(wigner_data, xaxis=None):
+    """Plots a curve for points in phase space of the spin Wigner
+        function
+
+    Args:
+        wigner_data(np.array): an array of points to plot as a 2d curve
+        xaxis (np.array):  the range of the x axis
+
+    Returns:
+        none: plot is shown with matplotlib to the screen
+
+    """
+    if xaxis is None:
+        xaxis = np.linspace(0,len(wigner_data)-1,num=len(wigner_data))
+
+    plt.plot(xaxis,wigner_data)
+    plt.show()
+
+def plot_wigner_plaquette(wigner_data, maxWigner='local'):
+    """Plots plaquette of wigner function data, the plaquette will
+    consist of cicles each colored to match the value of the Wigner
+    function at the given point in phase space.
+
+    Args:
+        wigner_data (matrix): array of Wigner function data where the
+                            rows are plotted along the x axis and the
+                            columns are plotted along the y axis
+        maxWigner:  - 'local' puts the maximum value to maximum of the points
+                    - 'unit' sets maximum to 1
+                    - float for a custom maximum.
+
+    Returns:
+        none: plot is shown with matplotlib to the screen
+
+    """
+    wigner_data = np.matrix(wigner_data)
+    dim = wigner_data.shape
+
+    if maxWigner == 'local':
+        Wmax = np.amax(wigner_data)
+    elif maxWigner == 'unit':
+        Wmax = 1
+    else:
+        Wmax = maxWigner #For a float input
+
+    cmap = matplotlib.cm.get_cmap('seismic_r')
+
+    xax = dim[1]-0.5
+    yax = dim[0]-0.5
+    norm = np.amax(dim)
+
+    fig = plt.figure(figsize=((xax+0.5)*6/norm,(yax+0.5)*6/norm))
+    ax = fig.gca()
+
+    for x in range(int(dim[1])):
+        for y in range(int(dim[0])):
+
+            circle = plt.Circle((x,y),0.49,
+                                color =
+                                cmap((wigner_data[y,x]+Wmax)/(2*Wmax)))
+            ax.add_artist(circle)
+
+    ax.set_xlim(-1, xax+0.5)
+    ax.set_ylim(-1, yax+0.5)
+    ax.set_xticks([], [])
+    ax.set_yticks([], [])
+    m = cm.ScalarMappable(cmap=cm.seismic_r)
+    m.set_array([-Wmax, Wmax])
+    plt.colorbar(m, shrink=0.5, aspect=10)
+    plt.show()
+
+def plot_wigner_data(wigner_data, phis=None,
+                     thetas=None, method=None,
+                     text_out=None):
+    """Plots Wigner results in appropriate format
+
+    Args:
+        wigner_data: Output returned from the wigner_data function
+        phis: Values of phi
+        thetas: Values of theta
+        method: how the data is to be plotted,
+            methods are:
+                point: a single point in phase space
+                curve: a two dimensional curve
+                plaquette: points plotted as circles
+
+    Returns:
+        none: plot is shown with matplotlib to the screen
+    """
+
+    if method is None:
+        wigDim = len(np.shape(wigner_data))
+        if wigDim == 1:
+            if np.shape(wigner_data) == 1:
+                method = 'point'
+            else:
+                method = 'curve'
+        elif wigDim ==2:
+            method ='plaquette'
+
+    if method == 'curve':
+        plot_wigner_curve(wigner_data, xaxis=phis)
+    elif method == 'plaquette':
+        plot_wigner_plaquette(wigner_data)
+    elif method == 'state':
+        wigner_function(wigner_data, text_out)
+    elif method == 'point':
+        plot_wigner_plaquette(wigner_data)
+        print('point in phase space is '+str(wigner_data))
     else:
         print("No method given")