Merge pull request #118 from rigetticomputing/feature/tomography-squa…

…shed

Feature/tomography

LICENSE

@@ -199,4 +199,35 @@
    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.
+   limitations under the License.
+----------------------------------------------------------------------------
+state_histogram in the tomography.utils module has adapted from the qutip library,
+ whose licence is below.
+
+ Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without 
+ modification, are permitted provided that the following conditions are met:
+
+ 1. Redistributions of source code must retain the above copyright notice, 
+ this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright notice, 
+ this list of conditions and the following disclaimer in the documentation 
+ and/or other materials provided with the distribution.
+
+ 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names of 
+ its contributors may be used to endorse or promote products derived from this 
+ software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 
+ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+ WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+ DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE 
+ FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
+ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
+ SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
+ OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
+ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

docs/_static/mathjax_macros.js

@@ -1,9 +1,14 @@
 MathJax.Hub.Config({
   TeX: {
     Macros: {
+      sket: ["\\left|\\left. #1 \\right\\rangle\\!\\right\\rangle",1],
+      sbra: ["\\left\\langle\\!\\left\\langle #1 \\right.\\right|",1],
+      sbraket: ["\\left\\langle\\!\\left\\langle #1 | #2 \\right\\rangle\\!\\right\\rangle",2],
       ket: ["\\left| #1 \\right\\rangle",1],
       bra: ["\\left\\langle #1 \\right|",1],
-      braket: ["\\left\\langle #1 | #2 \\right\\rangle",2]
+      braket: ["\\left\\langle #1 | #2 \\right\\rangle",2],
+      vec: ["\\text{vec}\\left(#1\\right)",1],
+      tr: ["\\text{Tr}\\left(#1\\right)",1]
     }
   }
 });

docs/index.rst

@@ -25,6 +25,7 @@ of which has its own self-contained documentation.
    qaoa
    qft
    phaseestimation
+   tomography
    grover
    bernstein_vazirani
    simon

grove/tests/tomography/test_operator_utils.py

@@ -0,0 +1,110 @@
+##############################################################################
+# Copyright 2017-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 matplotlib
+matplotlib.use('Agg')
+
+import numpy as np
+from scipy.sparse import csr_matrix
+import qutip as qt
+
+import grove.tomography.utils as ut
+from grove.tomography.operator_utils import to_realimag, FROBENIUS, is_projector, EPS, choi_matrix
+
+
+def test_operator_basis():
+    assert ut.PAULI_BASIS.all_hermitian()
+    assert ut.PAULI_BASIS.is_orthonormal()
+    assert is_projector(ut.GS)
+
+    two_qubit_pauli = ut.PAULI_BASIS.product(ut.PAULI_BASIS)
+    assert two_qubit_pauli.all_hermitian()
+    assert two_qubit_pauli.is_orthonormal()
+
+    sp = ut.PAULI_BASIS.super_basis()
+    assert sp.all_hermitian()
+    assert sp.is_orthonormal()
+
+    squared_pauli_basis = ut.PAULI_BASIS ** 2
+    for (l1, o1), (l2, o2) in zip(two_qubit_pauli, squared_pauli_basis):
+        assert l1 == l2
+        assert (o1 - o2).norm(FROBENIUS) < EPS
+
+    assert np.allclose(ut.PAULI_BASIS.basis_transform.T.toarray() * np.sqrt(2),
+                       np.array([[1., 0, 0, 1], [0, 1, 1, 0], [0, 1j, -1j, 0], [1, 0, 0, -1]]))
+
+    sX = qt.to_super(ut.QX)
+    tmX = ut.PAULI_BASIS.transfer_matrix(sX).toarray()
+    assert np.allclose(tmX, np.diag([1, 1, -1, -1]))
+    assert (sX - ut.PAULI_BASIS.super_from_tm(tmX)).norm(FROBENIUS) < EPS
+
+    pb3 = ut.PAULI_BASIS**3
+    assert pb3.dim == 4**3
+    assert pb3 == ut.n_qubit_pauli_basis(3)
+
+    assert ut.PAULI_BASIS**1 == ut.PAULI_BASIS
+
+    assert np.allclose(ut.PAULI_BASIS.project_op(ut.GS).toarray().ravel(),
+                       np.array([1, 0, 0, 1])/np.sqrt(2))
+
+    assert str(ut.PAULI_BASIS) == "<span[I,X,Y,Z]>"
+
+
+def test_super_operator_tools():
+    X, Y, Z, I = ut.QX, ut.QY, ut.QZ, ut.QI
+    bs = (I, X, Y, Z)
+
+    Xs = qt.sprepost(X, X)
+    # verify that Y+XYX==0 ( or XYX==-Y)
+    assert (Y + Xs(Y)).norm(FROBENIUS) < EPS
+
+    ptmX = np.array([[(bj * Xs(bk)).tr().real / 2 for bk in bs] for bj in bs])
+    assert np.allclose(ptmX, ut.PAULI_BASIS.transfer_matrix(Xs).toarray())
+
+    xchoi = qt.super_to_choi(Xs)
+    my_xchoi = choi_matrix(ptmX, ut.PAULI_BASIS)
+    assert (my_xchoi - xchoi).norm(FROBENIUS) < EPS
+
+    ys = qt.sprepost(Y, Y)
+    ptm_y = np.array([[(bj * ys(bk)).tr().real / 2 for bk in bs] for bj in bs])
+    assert np.allclose(ptm_y, ut.PAULI_BASIS.transfer_matrix(ys).toarray())
+
+    ychoi = qt.super_to_choi(ys)
+    my_ychoi = choi_matrix(ptm_y, ut.PAULI_BASIS)
+    assert (my_ychoi - ychoi).norm(FROBENIUS) < EPS
+
+    y2 = (-.25j * np.pi * Y).expm()
+    y2s = qt.sprepost(y2, y2.dag())
+    ptm_y2 = np.array([[(bj * y2s(bk)).tr().real / 2 for bk in bs] for bj in bs])
+    assert np.allclose(ptm_y2, ut.PAULI_BASIS.transfer_matrix(y2s).toarray())
+
+    y2choi = qt.super_to_choi(y2s)
+    my_y2choi = choi_matrix(ptm_y2, ut.PAULI_BASIS)
+    assert (my_y2choi - y2choi).norm(FROBENIUS) < EPS
+
+
+def test_to_realimag():
+    op = ut.QX + ut.QY
+    res = to_realimag(op)
+    assert isinstance(res, csr_matrix)
+    rd = res.toarray()
+    assert np.allclose(rd[:2, :2], [[0, 1], [1, 0]])
+    assert np.allclose(rd[:2, 2:], [[0, -1], [1, 0]])
+    assert np.allclose(rd[2:, :2], [[0, 1], [-1, 0]])
+    assert np.allclose(rd[2:, 2:], [[0, 1], [1, 0]])
+
+    res2 = to_realimag(op.data)
+    assert np.allclose(rd, res2.toarray())

grove/tests/tomography/test_process_tomography.py

@@ -0,0 +1,132 @@
+##############################################################################
+# Copyright 2017-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 matplotlib
+matplotlib.use('Agg')
+import pytest
+import os
+import numpy as np
+from mock import patch
+from mock import MagicMock, Mock
+import qutip as qt
+import json
+
+from pyquil.api import QVMConnection
+
+from grove.tomography.tomography import (MAX_QUBITS_PROCESS_TOMO,
+                                         default_channel_ops)
+from grove.tomography.process_tomography import (DEFAULT_PROCESS_TOMO_SETTINGS,
+                                                 process_tomography_programs,
+                                                 do_process_tomography, ProcessTomography,
+                                                 COMPLETELY_POSITIVE)
+from grove.tomography.process_tomography import (TRACE_PRESERVING)
+from grove.tomography.utils import (POVM_PI_BASIS, make_histogram,
+                                    sample_bad_readout, basis_state_preps,
+                                    estimate_assignment_probs, BAD_2Q_READOUT, SEED,
+                                    EPS, CNOT_PROGRAM)
+from grove.tomography.operator_utils import make_diagonal_povm
+
+
+SHOTS_PATH = os.path.join(os.path.dirname(__file__), 'process_shots.json')
+RESULTS_PATH = os.path.join(os.path.dirname(__file__), 'process_results.json')
+sample_bad_readout = MagicMock(sample_bad_readout)
+sample_bad_readout.side_effect = [np.array(shots) for shots in json.load(open(SHOTS_PATH, 'r'))]
+
+cxn = MagicMock(QVMConnection)
+cxn.run_and_measure.side_effect = json.load(open(RESULTS_PATH, 'r'))
+
+
+def test_process_tomography():
+    num_qubits = len(CNOT_PROGRAM.get_qubits())
+    dimension = 2 ** num_qubits
+
+    tomo_seq = list(process_tomography_programs(CNOT_PROGRAM))
+    nsamples = 3000
+
+    np.random.seed(SEED)
+    # We need more samples on the readout to ensure convergence.
+    state_prep_hists = [make_histogram(sample_bad_readout(p, 2 * nsamples, BAD_2Q_READOUT, cxn),
+                                       dimension) for p in basis_state_preps(*range(num_qubits))]
+    assignment_probs = estimate_assignment_probs(state_prep_hists)
+
+    histograms = np.zeros((len(tomo_seq), dimension))
+
+    for jj, p in enumerate(tomo_seq):
+        histograms[jj] = make_histogram(sample_bad_readout(p, nsamples, BAD_2Q_READOUT, cxn),
+                                        dimension)
+
+    channel_ops = list(default_channel_ops(num_qubits))
+    histograms = histograms.reshape((len(channel_ops), len(channel_ops), dimension))
+
+    povm = make_diagonal_povm(POVM_PI_BASIS ** num_qubits, assignment_probs)
+    cnot_ideal = qt.cnot()
+    for settings in [
+        DEFAULT_PROCESS_TOMO_SETTINGS,
+        DEFAULT_PROCESS_TOMO_SETTINGS._replace(constraints={TRACE_PRESERVING}),
+        DEFAULT_PROCESS_TOMO_SETTINGS._replace(constraints={TRACE_PRESERVING, COMPLETELY_POSITIVE}),
+    ]:
+
+        process_tomo = ProcessTomography.estimate_from_ssr(histograms, povm, channel_ops,
+                                                           channel_ops,
+                                                           settings)
+
+        assert abs(1 - process_tomo.avg_gate_fidelity(cnot_ideal)) < EPS
+
+        transfer_matrix = process_tomo.pauli_basis.transfer_matrix(qt.to_super(cnot_ideal))
+        assert abs(1 - process_tomo.avg_gate_fidelity(transfer_matrix)) < EPS
+        chi_rep = process_tomo.to_chi().data.toarray()
+        # When comparing to the identity, the chi representation is quadratically larger than the
+        # Hilbert space representation, so we take a square root.
+        probabilty_scale = np.sqrt(chi_rep.shape[0])
+        super_op_from_chi = np.zeros(process_tomo.pauli_basis.ops[0].shape, dtype=np.complex128)
+        for i, si in enumerate(process_tomo.pauli_basis.ops):
+            for j, sj in enumerate(process_tomo.pauli_basis.ops):
+                contribution = chi_rep[i][j] * si.data.toarray().conj().T.dot(sj.data.toarray())
+                super_op_from_chi += contribution / probabilty_scale
+        assert np.isclose(np.eye(process_tomo.pauli_basis.ops[0].shape[0]), super_op_from_chi,
+                          atol=EPS).all()
+        choi_rep = process_tomo.to_choi()
+
+        # Choi matrix should be a valid density matrix, scaled by the dimension of the system.
+        assert np.isclose(np.trace(choi_rep.data.toarray()) / probabilty_scale, 1, atol=EPS)
+
+        super_op = process_tomo.to_super()
+        # The map should be trace preserving.
+        assert np.isclose(np.sum(super_op[0]), 1, atol=EPS)
+
+    assert abs(1 - process_tomo.avg_gate_fidelity(qt.to_super(cnot_ideal))) < EPS
+
+    with patch("grove.tomography.utils.plot_pauli_transfer_matrix"), \
+         patch("grove.tomography.process_tomography.plt") as mplt:
+        mplt.subplots.return_value = Mock(), Mock()
+        process_tomo.plot()
+
+
+def test_do_process_tomography():
+    nsamples = 3000
+    qubits = list(range(MAX_QUBITS_PROCESS_TOMO + 1))
+    # Test with too many qubits.
+    with pytest.raises(ValueError):
+        _ = do_process_tomography(CNOT_PROGRAM, nsamples,
+                                  cxn, qubits)
+    process_tomo, assignment_probs, histograms = do_process_tomography(CNOT_PROGRAM, nsamples, cxn)
+    cnot_ideal = qt.cnot()
+    assert abs(1 - process_tomo.avg_gate_fidelity(cnot_ideal)) < EPS
+    for histogram_collection in histograms:
+        for histogram in histogram_collection:
+            assert np.sum(histogram) == nsamples
+    num_qubits = len(CNOT_PROGRAM.get_qubits())
+    assert np.isclose(assignment_probs, np.eye(2 ** num_qubits), atol=EPS).all()