Expose QubitPauliTensor for docs and tableau methods (#986)

* Expose QubitPauliTensor and test tableau methods

* Remove QubitPauliMap from docstring

* Added to changelog

pytket/binders/pauli.cpp

@@ -50,7 +50,9 @@ PYBIND11_MODULE(pauli, m) {
           py::arg("qubits"), py::arg("paulis"))
       .def(
           py::init<QubitPauliMap>(),
-          "Construct a QubitPauliString from a QubitPauliMap.", py::arg("map"))
+          "Construct a QubitPauliString from a dictionary mapping "
+          ":py:class:`Qubit` to :py:class:`Pauli`.",
+          py::arg("map"))
       .def(
           "__hash__",
           [](const QubitPauliString &qps) { return hash_value(qps); })
@@ -211,7 +213,7 @@ PYBIND11_MODULE(pauli, m) {
       m, "PauliStabiliser",
       "A string of Pauli letters from the alphabet {I, X, Y, Z} "
       "with a +/- 1 coefficient.")
-      .def(py::init<>(), "Constructs an empty QubitPauliString.")
+      .def(py::init<>(), "Constructs an empty PauliStabiliser.")
       .def(
           py::init([](const std::vector<Pauli> &string, const int &coeff) {
             if (coeff == 1) {
@@ -237,6 +239,185 @@ PYBIND11_MODULE(pauli, m) {
           "The list of Pauli terms")
       .def("__eq__", &PauliStabiliser::operator==)
       .def("__ne__", &PauliStabiliser::operator!=);
+
+  py::class_<QubitPauliTensor>(
+      m, "QubitPauliTensor",
+      "A tensor formed by Pauli terms, consisting of a sparse map from "
+      ":py:class:`Qubit` to :py:class:`Pauli` (implemented as a "
+      ":py:class:`QubitPauliString`) and a complex coefficient.")
+      .def(
+          py::init<Complex>(),
+          "Constructs an empty QubitPauliTensor, representing the identity.",
+          py::arg("coeff") = 1.)
+      .def(
+          py::init<Qubit, Pauli, Complex>(),
+          "Constructs a QubitPauliTensor with a single Pauli term.",
+          py::arg("qubit"), py::arg("pauli"), py::arg("coeff") = 1.)
+      .def(
+          py::init([](const std::list<Qubit> &qubits,
+                      const std::list<Pauli> &paulis, const Complex &coeff) {
+            return QubitPauliTensor(QubitPauliString(qubits, paulis), coeff);
+          }),
+          "Constructs a QubitPauliTensor from two matching lists of "
+          "Qubits and Paulis.",
+          py::arg("qubits"), py::arg("paulis"), py::arg("coeff") = 1.)
+      .def(
+          py::init<QubitPauliMap, Complex>(),
+          "Construct a QubitPauliTensor from a dictionary mapping "
+          ":py:class:`Qubit` to :py:class:`Pauli`.",
+          py::arg("map"), py::arg("coeff") = 1.)
+      .def(
+          py::init<QubitPauliString, Complex>(),
+          "Construct a QubitPauliTensor from a QubitPauliString.",
+          py::arg("string"), py::arg("coeff") = 1.)
+      .def(
+          "__hash__",
+          [](const QubitPauliTensor &qps) { return hash_value(qps); })
+      .def("__repr__", &QubitPauliTensor::to_str)
+      .def("__eq__", &QubitPauliTensor::operator==)
+      .def("__ne__", &QubitPauliTensor::operator!=)
+      .def("__lt__", &QubitPauliTensor::operator<)
+      .def(
+          "__getitem__", [](const QubitPauliTensor &qpt,
+                            const Qubit &q) { return qpt.string.get(q); })
+      .def(
+          "__setitem__", [](QubitPauliTensor &qpt, const Qubit &q,
+                            Pauli p) { return qpt.string.set(q, p); })
+      .def(py::self * py::self)
+      .def(Complex() * py::self)
+      .def_readwrite(
+          "string", &QubitPauliTensor::string,
+          "The QubitPauliTensor's underlying :py:class:`QubitPauliString`")
+      .def_readwrite(
+          "coeff", &QubitPauliTensor::coeff,
+          "The global coefficient of the tensor")
+      .def(
+          "compress", &QubitPauliTensor::compress,
+          "Removes I terms to compress the sparse representation.")
+      .def(
+          "commutes_with", &QubitPauliTensor::commutes_with,
+          ":return: True if the two tensors commute, else False",
+          py::arg("other"))
+      .def(
+          "to_sparse_matrix",
+          [](const QubitPauliTensor &qpt) {
+            return (CmplxSpMat)(qpt.coeff * qpt.string.to_sparse_matrix());
+          },
+          "Represents the sparse string as a dense string (without "
+          "padding for extra qubits) and generates the matrix for the "
+          "tensor. Uses the ILO-BE convention, so ``Qubit(\"a\", 0)`` "
+          "is more significant that ``Qubit(\"a\", 1)`` and "
+          "``Qubit(\"b\")`` for indexing into the matrix."
+          "\n\n:return: a sparse matrix corresponding to the tensor")
+      .def(
+          "to_sparse_matrix",
+          [](const QubitPauliTensor &qpt, unsigned n_qubits) {
+            return (CmplxSpMat)(qpt.coeff *
+                                qpt.string.to_sparse_matrix(n_qubits));
+          },
+          "Represents the sparse string as a dense string over "
+          "`n_qubits` qubits (sequentially indexed from 0 in the "
+          "default register) and generates the matrix for the tensor. "
+          "Uses the ILO-BE convention, so ``Qubit(0)`` is the most "
+          "significant bit for indexing into the matrix."
+          "\n\n:param n_qubits: the number of qubits in the full "
+          "operator"
+          "\n:return: a sparse matrix corresponding to the operator",
+          py::arg("n_qubits"))
+      .def(
+          "to_sparse_matrix",
+          [](const QubitPauliTensor &qpt, const qubit_vector_t &qubits) {
+            return (CmplxSpMat)(qpt.coeff *
+                                qpt.string.to_sparse_matrix(qubits));
+          },
+          "Represents the sparse string as a dense string and generates "
+          "the matrix for the tensor. Orders qubits according to "
+          "`qubits` (padding with identities if they are not in the "
+          "sparse string), so ``qubits[0]`` is the most significant bit "
+          "for indexing into the matrix."
+          "\n\n:param qubits: the ordered list of qubits in the full "
+          "operator"
+          "\n:return: a sparse matrix corresponding to the operator",
+          py::arg("qubits"))
+      .def(
+          "dot_state",
+          [](const QubitPauliTensor &qpt, const Eigen::VectorXcd &state) {
+            return qpt.coeff * qpt.string.dot_state(state);
+          },
+          "Performs the dot product of the state with the pauli tensor. "
+          "Maps the qubits of the statevector with sequentially-indexed "
+          "qubits in the default register, with ``Qubit(0)`` being the "
+          "most significant qubit."
+          "\n\n:param state: statevector for qubits ``Qubit(0)`` to "
+          "``Qubit(n-1)``"
+          "\n:return: dot product of operator with state",
+          py::arg("state"))
+      .def(
+          "dot_state",
+          [](const QubitPauliTensor &qpt, const Eigen::VectorXcd &state,
+             const qubit_vector_t &qubits) {
+            return qpt.coeff * qpt.string.dot_state(state, qubits);
+          },
+          "Performs the dot product of the state with the pauli tensor. "
+          "Maps the qubits of the statevector according to the ordered "
+          "list `qubits`, with ``qubits[0]`` being the most significant "
+          "qubit."
+          "\n\n:param state: statevector"
+          "\n:param qubits: order of qubits in `state` from most to "
+          "least significant"
+          "\n:return: dot product of operator with state",
+          py::arg("state"), py::arg("qubits"))
+      .def(
+          "state_expectation",
+          [](const QubitPauliTensor &qpt, const Eigen::VectorXcd &state) {
+            return qpt.coeff * qpt.string.state_expectation(state);
+          },
+          "Calculates the expectation value of the state with the pauli "
+          "operator. Maps the qubits of the statevector with "
+          "sequentially-indexed qubits in the default register, with "
+          "``Qubit(0)`` being the most significant qubit."
+          "\n\n:param state: statevector for qubits ``Qubit(0)`` to "
+          "``Qubit(n-1)``"
+          "\n:return: expectation value with respect to state",
+          py::arg("state"))
+      .def(
+          "state_expectation",
+          [](const QubitPauliTensor &qpt, const Eigen::VectorXcd &state,
+             const qubit_vector_t &qubits) {
+            return qpt.coeff * qpt.string.state_expectation(state, qubits);
+          },
+          "Calculates the expectation value of the state with the pauli "
+          "operator. Maps the qubits of the statevector according to the "
+          "ordered list `qubits`, with ``qubits[0]`` being the most "
+          "significant qubit."
+          "\n\n:param state: statevector"
+          "\n:param qubits: order of qubits in `state` from most to "
+          "least significant"
+          "\n:return: expectation value with respect to state",
+          py::arg("state"), py::arg("qubits"))
+
+      .def(py::pickle(
+          [](const QubitPauliTensor &qpt) {
+            std::list<Qubit> qubits;
+            std::list<Pauli> paulis;
+            for (const std::pair<const Qubit, Pauli> &qp_pair :
+                 qpt.string.map) {
+              qubits.push_back(qp_pair.first);
+              paulis.push_back(qp_pair.second);
+            }
+            return py::make_tuple(qubits, paulis, qpt.coeff);
+          },
+          [](const py::tuple &t) {
+            if (t.size() != 3)
+              throw std::runtime_error(
+                  "Invalid state: tuple size: " + std::to_string(t.size()));
+            return QubitPauliTensor(
+                QubitPauliString(
+                    t[0].cast<std::list<Qubit>>(),
+                    t[1].cast<std::list<Pauli>>()),
+                t[2].cast<Complex>());
+          }));
+  ;
 }
 
 }  // namespace tket

pytket/binders/tableau.cpp

@@ -51,6 +51,14 @@ PYBIND11_MODULE(tableau, m) {
           "\n:param zph: The phases of the Z rows.",
           py::arg("xx"), py::arg("xz"), py::arg("xph"), py::arg("zx"),
           py::arg("zz"), py::arg("zph"))
+      .def(
+          py::init<>([](const Circuit& circ) {
+            return circuit_to_unitary_tableau(circ);
+          }),
+          "Constructs a :py:class:`UnitaryTableau` from a unitary "
+          ":py:class:`Circuit`. Throws an exception if the input contains "
+          "non-unitary operations."
+          "\n\n:param circ: The unitary circuit to convert to a tableau.")
       .def(
           "__repr__",
           [](const UnitaryTableau& tab) {
@@ -94,7 +102,14 @@ PYBIND11_MODULE(tableau, m) {
           "\n\n:param type: The :py:class:`OpType` of the gate to add. Must be "
           "an unparameterised Clifford gate type."
           "\n:param qbs: The qubits to apply the gate to. Length must match "
-          "the arity of the given gate type.");
+          "the arity of the given gate type.")
+      .def(
+          "to_circuit", &unitary_tableau_to_circuit,
+          "Synthesises a unitary :py:class:`Circuit` realising the same "
+          "unitary as the tableau. Uses the method from Aaronson & Gottesman: "
+          "\"Improved Simulation of Stabilizer Circuits\", Theorem 8. This is "
+          "not optimised for gate count, so is not recommended for "
+          "performance-sensitive usage.");
   py::class_<UnitaryTableauBox, std::shared_ptr<UnitaryTableauBox>, Op>(
       m, "UnitaryTableauBox",
       "A Clifford unitary specified by its actions on Paulis.")

pytket/docs/changelog.rst

@@ -18,6 +18,9 @@ Minor new features:
   be overridden using the ``always_squash_symbols`` parameter to
   ``SquashCustom``.
 * Add ``control_state`` argument to ``QControlBox``.
+* Add ``QubitPauliTensor`` (combining ``QubitPauliString`` with a complex
+  coefficient) to python binding. This is incorporated into ``UnitaryTableau`` 
+  row inspection for phase tracking.
 
 Fixes:
 

pytket/tests/tableau_test.py

@@ -14,6 +14,7 @@
 
 import pytest  # type: ignore
 from pytket.circuit import Circuit, OpType, Qubit  # type: ignore
+from pytket.pauli import Pauli, QubitPauliTensor  # type: ignore
 from pytket.tableau import UnitaryTableau, UnitaryTableauBox  # type: ignore
 from pytket.utils.results import compare_unitaries
 import numpy as np
@@ -81,3 +82,22 @@ def test_tableau_box_from_matrix() -> None:
     circ.X(2)
     circ.Sdg(2)
     assert compare_unitaries(circ.get_unitary(), np.eye(8, dtype=complex))
+
+
+def test_tableau_rows() -> None:
+    circ = Circuit(3)
+    circ.H(0)
+    circ.CX(0, 1)
+    circ.V(1)
+    circ.CZ(2, 1)
+    circ.Vdg(2)
+    tab = UnitaryTableau(circ)
+    assert tab.get_zrow(Qubit(0)) == QubitPauliTensor(
+        [Qubit(0), Qubit(1), Qubit(2)], [Pauli.X, Pauli.X, Pauli.Y], 1.0
+    )
+    assert tab.get_xrow(Qubit(1)) == QubitPauliTensor(
+        {Qubit(1): Pauli.X, Qubit(2): Pauli.Y}, 1.0
+    )
+    assert tab.get_row_product(
+        QubitPauliTensor(Qubit(0), Pauli.Z) * QubitPauliTensor(Qubit(1), Pauli.X, -1.0)
+    ) == QubitPauliTensor(Qubit(0), Pauli.X, -1.0)