diff --git a/pytket/conanfile.py b/pytket/conanfile.py index b64029ee6f..8156e6834a 100644 --- a/pytket/conanfile.py +++ b/pytket/conanfile.py @@ -32,7 +32,7 @@ def package(self): cmake.install() def requirements(self): - self.requires("tket/1.2.108@tket/stable") + self.requires("tket/1.2.109@tket/stable") self.requires("tklog/0.3.3@tket/stable") self.requires("tkrng/0.3.3@tket/stable") self.requires("tkassert/0.3.4@tket/stable") diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt index 5d534e4316..0185f4ece6 100644 --- a/tket/CMakeLists.txt +++ b/tket/CMakeLists.txt @@ -225,7 +225,6 @@ target_sources(tket src/ZX/MBQCRewrites.cpp src/ZX/ZXRWSequences.cpp src/Converters/ChoiMixTableauConverters.cpp - src/Converters/PauliGadget.cpp src/Converters/PauliGraphConverters.cpp src/Converters/Gauss.cpp src/Converters/PhasePoly.cpp @@ -379,7 +378,6 @@ target_sources(tket include/tket/ZX/ZXGenerator.hpp include/tket/Converters/Converters.hpp include/tket/Converters/Gauss.hpp - include/tket/Converters/PauliGadget.hpp include/tket/Converters/PhasePoly.hpp include/tket/Converters/UnitaryTableauBox.hpp include/tket/Placement/NeighbourPlacements.hpp diff --git a/tket/conanfile.py b/tket/conanfile.py index b7e268ffba..a564cddd92 100644 --- a/tket/conanfile.py +++ b/tket/conanfile.py @@ -23,7 +23,7 @@ class TketConan(ConanFile): name = "tket" - version = "1.2.108" + version = "1.2.109" package_type = "library" license = "Apache 2" homepage = "https://github.com/CQCL/tket" diff --git a/tket/include/tket/Circuit/CircUtils.hpp b/tket/include/tket/Circuit/CircUtils.hpp index ae04b26ada..281ce94eb7 100644 --- a/tket/include/tket/Circuit/CircUtils.hpp +++ b/tket/include/tket/Circuit/CircUtils.hpp @@ -111,13 +111,13 @@ std::pair decompose_2cx_DV(const Eigen::Matrix4cd& U); * Construct a phase gadget * * @param n_qubits number of qubits - * @param t phase parameter + * @param angle phase parameter * @param cx_config CX configuration * * @return phase gadget implementation wrapped in a ConjugationBox */ Circuit phase_gadget( - unsigned n_qubits, const Expr& t, + unsigned n_qubits, const Expr& angle, CXConfigType cx_config = CXConfigType::Snake); /** @@ -127,13 +127,29 @@ Circuit phase_gadget( * \f$ e^{-\frac12 i \pi t \sigma_0 \otimes \sigma_1 \otimes \cdots} \f$ * where \f$ \sigma_i \in \{I,X,Y,Z\} \f$ are the Pauli operators. * - * @param paulis Pauli operators - * @param t angle in half-turns + * @param pauli Pauli operators; coefficient gives rotation angle in half-turns * @param cx_config CX configuration * @return Pauli gadget implementation wrapped in a ConjugationBox */ Circuit pauli_gadget( - const std::vector& paulis, const Expr& t, + SpSymPauliTensor pauli, CXConfigType cx_config = CXConfigType::Snake); + +/** + * Construct a circuit realising a pair of Pauli gadgets with the fewest + * two-qubit gates. + * + * The returned circuit implements the unitary e^{-i pi angle1 paulis1 / 2} + * e^{-i pi angle0 paulis0 / 2}, i.e. a gadget of angle0 about paulis0 followed + * by a gadget of angle1 about paulis1. + * + * @param paulis0 Pauli operators for first gadget; coefficient gives rotation + * angle in half-turns + * @param paulis1 Pauli operators for second gadget; coefficient gives rotation + * angle in half-turns + * @param cx_config CX configuration + */ +Circuit pauli_gadget_pair( + SpSymPauliTensor paulis0, SpSymPauliTensor paulis1, CXConfigType cx_config = CXConfigType::Snake); /** diff --git a/tket/include/tket/Circuit/PauliExpBoxes.hpp b/tket/include/tket/Circuit/PauliExpBoxes.hpp index 325828f577..151bb12849 100644 --- a/tket/include/tket/Circuit/PauliExpBoxes.hpp +++ b/tket/include/tket/Circuit/PauliExpBoxes.hpp @@ -275,4 +275,46 @@ class TermSequenceBox : public Box { CXConfigType cx_configuration_; }; +/** + * Constructs a PauliExpBox for a single pauli gadget and appends it to a + * circuit. + * + * @param circ The circuit to append the box to + * @param pauli The pauli operator of the gadget; coefficient gives the rotation + * angle in half-turns + * @param cx_config The CX configuration to be used during synthesis + */ +void append_single_pauli_gadget_as_pauli_exp_box( + Circuit &circ, const SpSymPauliTensor &pauli, CXConfigType cx_config); + +/** + * Constructs a PauliExpPairBox for a pair of pauli gadgets and appends it to a + * circuit. The pauli gadgets may or may not commute, so the ordering matters. + * + * @param circ The circuit to append the box to + * @param pauli0 The pauli operator of the first gadget; coefficient gives the + * rotation angle in half-turns + * @param pauli1 The pauli operator of the second gadget; coefficient gives the + * rotation angle in half-turns + * @param cx_config The CX configuration to be used during synthesis + */ +void append_pauli_gadget_pair_as_box( + Circuit &circ, const SpSymPauliTensor &pauli0, + const SpSymPauliTensor &pauli1, CXConfigType cx_config); + +/** + * Constructs a PauliExpCommutingSetBox for a set of mutually commuting pauli + * gadgets and appends it to a circuit. As the pauli gadgets all commute, the + * ordering does not matter semantically, but may yield different synthesised + * circuits. + * + * @param circ The circuit to append the box to + * @param gadgets Description of the pauli gadgets; coefficients give the + * rotation angles in half-turns + * @param cx_config The CX configuration to be used during synthesis + */ +void append_commuting_pauli_gadget_set_as_box( + Circuit &circ, const std::list &gadgets, + CXConfigType cx_config); + } // namespace tket diff --git a/tket/include/tket/Clifford/ChoiMixTableau.hpp b/tket/include/tket/Clifford/ChoiMixTableau.hpp index ea2e0bfc51..8ecb7c8808 100644 --- a/tket/include/tket/Clifford/ChoiMixTableau.hpp +++ b/tket/include/tket/Clifford/ChoiMixTableau.hpp @@ -45,7 +45,8 @@ class ChoiMixTableau { * When mapped to a sparse readable representation, independent * SpPauliStabiliser objects are used for each segment, so we no longer expect * their individual phases to be +-1, instead only requiring this on their - * product. + * product. get_row() will automatically transpose the input segment term so + * it is presented as RxS s.t. SCR = C. * * Columns of the tableau are indexed by pair of Qubit id and a tag to * distinguish input vs output. Rows are not maintained in any particular @@ -93,6 +94,7 @@ class ChoiMixTableau { * Construct a tableau directly from its rows. * Each row is represented as a product of SpPauliStabilisers where the first * is over the input qubits and the second is over the outputs. + * A row RxS is a pair s.t. SCR = C */ explicit ChoiMixTableau(const std::list& rows); /** @@ -122,13 +124,23 @@ class ChoiMixTableau { * Get the number of boundaries representing outputs from the process. */ unsigned get_n_outputs() const; + /** + * Get all qubit names present in the input segment. + */ + qubit_vector_t input_qubits() const; + /** + * Get all qubit names present in the output segment. + */ + qubit_vector_t output_qubits() const; /** - * Read off a row as a Pauli string + * Read off a row as a Pauli string. + * Returns a pair of Pauli strings RxS such that SCR = C */ row_tensor_t get_row(unsigned i) const; /** - * Combine rows into a single row + * Combine rows into a single row. + * Returns a pair of Pauli strings RxS such that SCR = C */ row_tensor_t get_row_product(const std::vector& rows) const; @@ -142,7 +154,10 @@ class ChoiMixTableau { * outputs. */ void apply_S(const Qubit& qb, TableauSegment seg = TableauSegment::Output); + void apply_Z(const Qubit& qb, TableauSegment seg = TableauSegment::Output); void apply_V(const Qubit& qb, TableauSegment seg = TableauSegment::Output); + void apply_X(const Qubit& qb, TableauSegment seg = TableauSegment::Output); + void apply_H(const Qubit& qb, TableauSegment seg = TableauSegment::Output); void apply_CX( const Qubit& control, const Qubit& target, TableauSegment seg = TableauSegment::Output); diff --git a/tket/include/tket/Clifford/SymplecticTableau.hpp b/tket/include/tket/Clifford/SymplecticTableau.hpp index 961bbbef06..b353e3b5dd 100644 --- a/tket/include/tket/Clifford/SymplecticTableau.hpp +++ b/tket/include/tket/Clifford/SymplecticTableau.hpp @@ -20,12 +20,6 @@ namespace tket { -// Forward declare friend classes for converters -class ChoiMixTableau; -class UnitaryTableau; -class UnitaryRevTableau; -class Circuit; - /** * Boolean encoding of Pauli * = ==> I @@ -136,10 +130,13 @@ class SymplecticTableau { void row_mult(unsigned ra, unsigned rw, Complex coeff = 1.); /** - * Applies an S/V/CX gate to the given qubit(s) + * Applies a chosen gate to the given qubit(s) */ void apply_S(unsigned qb); + void apply_Z(unsigned qb); void apply_V(unsigned qb); + void apply_X(unsigned qb); + void apply_H(unsigned qb); void apply_CX(unsigned qc, unsigned qt); void apply_gate(OpType type, const std::vector &qbs); @@ -173,29 +170,19 @@ class SymplecticTableau { */ void gaussian_form(); - private: - /** - * Number of rows - */ - unsigned n_rows_; - - /** - * Number of qubits in each row - */ - unsigned n_qubits_; - /** * Tableau contents */ - MatrixXb xmat_; - MatrixXb zmat_; - VectorXb phase_; + MatrixXb xmat; + MatrixXb zmat; + VectorXb phase; /** * Complex conjugate of the state by conjugating rows */ SymplecticTableau conjugate() const; + private: /** * Helper methods for manipulating the tableau when applying gates */ @@ -206,18 +193,6 @@ class SymplecticTableau { void col_mult( const MatrixXb::ColXpr &a, const MatrixXb::ColXpr &b, bool flip, MatrixXb::ColXpr &w, VectorXb &pw); - - friend class UnitaryTableau; - friend class ChoiMixTableau; - friend Circuit unitary_tableau_to_circuit(const UnitaryTableau &tab); - friend std::pair cm_tableau_to_circuit( - const ChoiMixTableau &tab); - friend std::ostream &operator<<(std::ostream &os, const UnitaryTableau &tab); - friend std::ostream &operator<<( - std::ostream &os, const UnitaryRevTableau &tab); - - friend void to_json(nlohmann::json &j, const SymplecticTableau &tab); - friend void from_json(const nlohmann::json &j, SymplecticTableau &tab); }; JSON_DECL(SymplecticTableau) diff --git a/tket/include/tket/Clifford/UnitaryTableau.hpp b/tket/include/tket/Clifford/UnitaryTableau.hpp index 0c3b780a5c..1740eed478 100644 --- a/tket/include/tket/Clifford/UnitaryTableau.hpp +++ b/tket/include/tket/Clifford/UnitaryTableau.hpp @@ -101,8 +101,14 @@ class UnitaryTableau { */ void apply_S_at_front(const Qubit& qb); void apply_S_at_end(const Qubit& qb); + void apply_Z_at_front(const Qubit& qb); + void apply_Z_at_end(const Qubit& qb); void apply_V_at_front(const Qubit& qb); void apply_V_at_end(const Qubit& qb); + void apply_X_at_front(const Qubit& qb); + void apply_X_at_end(const Qubit& qb); + void apply_H_at_front(const Qubit& qb); + void apply_H_at_end(const Qubit& qb); void apply_CX_at_front(const Qubit& control, const Qubit& target); void apply_CX_at_end(const Qubit& control, const Qubit& target); void apply_gate_at_front(OpType type, const qubit_vector_t& qbs); @@ -244,8 +250,14 @@ class UnitaryRevTableau { */ void apply_S_at_front(const Qubit& qb); void apply_S_at_end(const Qubit& qb); + void apply_Z_at_front(const Qubit& qb); + void apply_Z_at_end(const Qubit& qb); void apply_V_at_front(const Qubit& qb); void apply_V_at_end(const Qubit& qb); + void apply_X_at_front(const Qubit& qb); + void apply_X_at_end(const Qubit& qb); + void apply_H_at_front(const Qubit& qb); + void apply_H_at_end(const Qubit& qb); void apply_CX_at_front(const Qubit& control, const Qubit& target); void apply_CX_at_end(const Qubit& control, const Qubit& target); void apply_gate_at_front(OpType type, const qubit_vector_t& qbs); diff --git a/tket/include/tket/Converters/Converters.hpp b/tket/include/tket/Converters/Converters.hpp index ff141a2432..a3de8c9146 100644 --- a/tket/include/tket/Converters/Converters.hpp +++ b/tket/include/tket/Converters/Converters.hpp @@ -47,13 +47,88 @@ ChoiMixTableau circuit_to_cm_tableau(const Circuit &circ); /** * Constructs a circuit producing the same effect as a ChoiMixTableau. - * Uses a naive synthesis method until we develop a good heuristic. * Since Circuit does not support distinct qubit addresses for inputs and * outputs, also returns a map from the output qubit IDs in the tableau to their - * corresponding outputs in the circuit + * corresponding outputs in the circuit. + * + * The circuit produced will be the (possibly non-unitary) channel whose + * stabilisers are exactly those of the tableau and no more, using + * initialisations, post-selections, discards, resets, and collapses to ensure + * this. It will automatically reuse qubits so no more qubits will be needed + * than max(tab.get_n_inputs(), tab.get_n_outputs()). + * + * Example 1: + * ZXI -> () + * YYZ -> () + * This becomes a diagonalisation circuit followed by post-selections. + * + * Example 2: + * Z -> ZZ + * X -> IY + * Z -> -XX + * Combining the first and last rows reveals an initialisation is required for I + * -> YY. Since there are two output qubits, at least one of them does not + * already exist in the input fragment so we can freely add an extra qubit on + * the input side, initialise it and apply a unitary mapping IZ -> YY. + * + * Example 3: + * ZX -> IZ + * II -> ZI + * We require an initialised qubit for the final row, but both input and output + * spaces only have q[0] and q[1], of which both inputs need to be open for the + * first row. We can obtain an initialised qubit by resetting a qubit after + * reducing the first row to only a single qubit. */ -std::pair cm_tableau_to_circuit( - const ChoiMixTableau &circ); +std::pair cm_tableau_to_exact_circuit( + const ChoiMixTableau &tab, CXConfigType cx_config = CXConfigType::Snake); + +/** + * We define a unitary extension of a ChoiMixTableau to be a unitary circuit + * whose stabilizer group contain all the rows of the ChoiMixTableau and + * possibly more. This is useful when we are treating the ChoiMixTableau as a + * means to encode a diagonalisation problem, since we are generally looking for + * a unitary as we may wish to apply the inverse afterwards (e.g. conjugating + * some rotations to implement a set of Pauli gadgets). + * + * Not every ChoiMixTableau can be extended to a unitary by just adding rows, + * e.g. if it requires any initialisation or post-selections. In this case, the + * unitary circuit is extended with additional input qubits which are assumed to + * be zero-initialised, and additional output qubits which are assumed to be + * post-selected. The synthesis guarantees that, if we take the unitary, + * initialise all designated inputs, and post-select on all designated outputs, + * every row from the original tableau is a stabiliser for the remaining + * projector. When not enough additional qubit names are provided, an error is + * thrown. + * + * + * Example 1: + * ZXI -> () + * YYZ -> () + * Since, in exact synthesis, at least two post-selections would be required, we + * pick two names from post_names. This is then a diagonalisation circuit which + * maps each row to an arbitrary diagonal string over post_names. + * + * Example 2: + * Z -> ZZ + * X -> IY + * Z -> -XX + * Combining the first and last rows reveals an initialisation is required for I + * -> YY. We extend the inputs with a qubit from init_names. The initialisation + * can manifest as either altering the first row to ZZ -> ZZ or the last row to + * ZZ -> -XX. + * + * Example 3: + * ZX -> IZ + * II -> ZI + * We require an initialised qubit for the final row, but both input and output + * spaces only have q[0] and q[1], of which both inputs need to be open for the + * first row. Unlike exact synthesis, we cannot reuse qubits, so the returned + * circuit will be over 3 qubits, extending with a name from init_names. + */ +std::pair cm_tableau_to_unitary_extension_circuit( + const ChoiMixTableau &tab, const std::vector &init_names = {}, + const std::vector &post_names = {}, + CXConfigType cx_config = CXConfigType::Snake); PauliGraph circuit_to_pauli_graph(const Circuit &circ); diff --git a/tket/include/tket/Converters/PauliGadget.hpp b/tket/include/tket/Converters/PauliGadget.hpp deleted file mode 100644 index ef09cba1a7..0000000000 --- a/tket/include/tket/Converters/PauliGadget.hpp +++ /dev/null @@ -1,83 +0,0 @@ -// Copyright 2019-2024 Cambridge Quantum 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. - -#pragma once - -#include - -#include "tket/Circuit/Circuit.hpp" -#include "tket/Utils/PauliTensor.hpp" - -namespace tket { - -class ImplicitPermutationNotAllowed : public std::logic_error { - public: - explicit ImplicitPermutationNotAllowed(const std::string& message) - : std::logic_error(message) {} -}; - -/** - * Append a Pauli gadget to the end of a given circuit. - * Automatically uses Snake CX configuration - * - * @param circ circuit to append to - * @param pauli Pauli operators and their respective qubits; coefficient gives - * rotation angle in half-turns - * @param cx_config which type of CX configuration to decompose into - */ -void append_single_pauli_gadget( - Circuit& circ, const SpSymPauliTensor& pauli, - CXConfigType cx_config = CXConfigType::Snake); - -/** - * Append a Pauli gadget to the end of a given circuit as a - * PauliExpBox. - * Automatically uses Snake CX configuration - * - * @param circ circuit to append to - * @param pauli Pauli operators and their respective qubits; coefficient gives - * rotation angle in half-turns - * @param cx_config which type of CX configuration to decompose into - */ -void append_single_pauli_gadget_as_pauli_exp_box( - Circuit& circ, const SpSymPauliTensor& pauli, CXConfigType cx_config); - -/** - * Append a pair of Pauli gadgets to the end of a given circuit. - * (shallow) Uses an adapted arrangement of CX that gives balanced trees - * over the matching qubits to improve depth. Better performance - * is not guaranteed as CXs may not align for cancellation and - * it can be harder to route. - * (!shallow) Uses the original method with naive arrangement of CXs. - * - * @param circ circuit to append to - * @param pauli0 first Pauli string; coefficient gives rotation angle in - * half-turns - * @param pauli1 second Pauli string; coefficient gives rotation angle in - * half-turns - * @param cx_config which type of CX configuration to decompose into - */ -void append_pauli_gadget_pair( - Circuit& circ, SpSymPauliTensor pauli0, SpSymPauliTensor pauli1, - CXConfigType cx_config = CXConfigType::Snake); - -void append_pauli_gadget_pair_as_box( - Circuit& circ, const SpSymPauliTensor& pauli0, - const SpSymPauliTensor& pauli1, CXConfigType cx_config); - -void append_commuting_pauli_gadget_set_as_box( - Circuit& circ, const std::list& gadgets, - CXConfigType cx_config); - -} // namespace tket diff --git a/tket/include/tket/Diagonalisation/Diagonalisation.hpp b/tket/include/tket/Diagonalisation/Diagonalisation.hpp index 3d43ccd260..93fabbe92c 100644 --- a/tket/include/tket/Diagonalisation/Diagonalisation.hpp +++ b/tket/include/tket/Diagonalisation/Diagonalisation.hpp @@ -60,13 +60,51 @@ void apply_conjugations( SpSymPauliTensor &qps, const Conjugations &conjugations); /** - * Given two qubits on which to conjugate a CX gate, try to conjugate with a - * XXPhase3 instead. If successful, undoes conjugations that must be undone and - * replaces it with XXPhase3 conjugation. Returns true if successful and false - * otherwise. + * Given a Pauli tensor P, produces a short Clifford circuit C which maps P to Z + * on a single qubit, i.e. Z_i C P = C. This can be viewed as the components + * required to synthesise a single Pauli gadget C^dag RZ(a)_i C = exp(-i pi a + * P/2) (up to global phase), or as a diagonalisation of a single Pauli string + * along with CXs to reduce it to a single qubit. Returns the circuit C and the + * qubit i where the Z ends up. */ -bool conjugate_with_xxphase3( - const Qubit &qb_a, const Qubit &qb_b, Conjugations &conjugations, - Circuit &cliff_circ); +std::pair reduce_pauli_to_z( + const SpPauliStabiliser &pauli, CXConfigType cx_config); + +/** + * Given a pair of (either commuting or anticommuting) Pauli tensors P0, P1, + * produces a short Clifford circuit C which maps P0 and P1 to strings which + * overlap on at most one qubit (which is also returned). If P0 and P1 + * anticommute, a mismatching qubit is always left. If they commute, contain at + * least one matching qubit and no mismatching qubits, then the final matching + * qubit is returned if allow_matching_final is true, otherwise P0 and P1 are + * reduced to non-overlapping strings. + */ +std::pair> reduce_overlap_of_paulis( + SpPauliStabiliser &pauli0, SpPauliStabiliser &pauli1, + CXConfigType cx_config, bool allow_matching_final = false); + +/** + * Given a pair of anticommuting Pauli tensors P0, P1, produces a short Clifford + * circuit C which maps P0 to Z and P1 to X on the same qubit, i.e. Z_i C P0 = C + * = X_i C P1. This can be viewed as the components required to synthesise a + * pair of noncommuting Pauli gadgets C^dag RX(b)_i RZ(a)_i C = exp(-i pi b + * P1/2) exp(-i pi a P0/2) (up to global phase). This is not strictly a + * diagonalisation because anticommuting strings cannot be simultaneously + * diagonalised. Returns the circuit C and the qubit i where the Z and X end up. + */ +std::pair reduce_anticommuting_paulis_to_z_x( + SpPauliStabiliser pauli0, SpPauliStabiliser pauli1, CXConfigType cx_config); + +/** + * Given a pair of commuting Pauli tensors P0, P1, produces a short Clifford + * circuit C which maps P0 and P1 to Z on different qubits, i.e. Z_i C P0 = C = + * Z_j C P1. This can be viewed as the components required to synthesise a pair + * of commuting Pauli gadgets C^dag RZ(b)_j RZ(a)_i C = exp(-i pi b P1/2) exp(-i + * pi a P0/2) (up to global phase), or as a mutual diagonalisation of two Pauli + * strings along with CXs to reduce them to independent, individual qubits. + * Returns the circuit C and the qubits i and j where the Zs end up. + */ +std::tuple reduce_commuting_paulis_to_zi_iz( + SpPauliStabiliser pauli0, SpPauliStabiliser pauli1, CXConfigType cx_config); } // namespace tket diff --git a/tket/src/Circuit/CircUtils.cpp b/tket/src/Circuit/CircUtils.cpp index 45c820b858..aa29bb9ee1 100644 --- a/tket/src/Circuit/CircUtils.cpp +++ b/tket/src/Circuit/CircUtils.cpp @@ -22,6 +22,7 @@ #include "tket/Circuit/CircPool.hpp" #include "tket/Circuit/Circuit.hpp" #include "tket/Circuit/ConjugationBox.hpp" +#include "tket/Diagonalisation/Diagonalisation.hpp" #include "tket/Gate/GatePtr.hpp" #include "tket/Gate/GateUnitaryMatrixImplementations.hpp" #include "tket/Gate/Rotation.hpp" @@ -269,61 +270,69 @@ std::pair decompose_2cx_DV(const Eigen::Matrix4cd &U) { } Circuit phase_gadget(unsigned n_qubits, const Expr &t, CXConfigType cx_config) { - // Handle n_qubits==0 as a special case, or the calculations below - // go badly wrong. - Circuit new_circ(n_qubits); - Circuit compute(n_qubits); - Circuit action(n_qubits); - Circuit uncompute(n_qubits); - if (n_qubits == 0) { - new_circ.add_phase(-t / 2); - return new_circ; + return pauli_gadget( + SpSymPauliTensor(DensePauliMap(n_qubits, Pauli::Z), t), cx_config); +} + +Circuit pauli_gadget(SpSymPauliTensor paulis, CXConfigType cx_config) { + Circuit circ; + for (const std::pair &qp : paulis.string) + circ.add_qubit(qp.first); + if (SpPauliString(paulis.string) == SpPauliString()) { + circ.add_phase(-paulis.coeff / 2); + return circ; + } + Circuit compute(circ); + Circuit action(circ); + qubit_vector_t qubits; + for (const std::pair &qp : paulis.string) { + switch (qp.second) { + case Pauli::I: + break; + case Pauli::X: + compute.add_op(OpType::H, {qp.first}); + qubits.push_back(qp.first); + break; + case Pauli::Y: + compute.add_op(OpType::V, {qp.first}); + qubits.push_back(qp.first); + break; + case Pauli::Z: + qubits.push_back(qp.first); + break; + } } + unsigned n_qubits = qubits.size(); switch (cx_config) { case CXConfigType::Snake: { for (unsigned i = n_qubits - 1; i != 0; --i) { unsigned j = i - 1; - compute.add_op(OpType::CX, {i, j}); - } - action.add_op(OpType::Rz, t, {0}); - for (unsigned i = 0; i != n_qubits - 1; ++i) { - unsigned j = i + 1; - uncompute.add_op(OpType::CX, {j, i}); + compute.add_op(OpType::CX, {qubits[i], qubits[j]}); } + action.add_op(OpType::Rz, paulis.coeff, {qubits[0]}); break; } case CXConfigType::Star: { for (unsigned i = n_qubits - 1; i != 0; --i) { - compute.add_op(OpType::CX, {i, 0}); - } - action.add_op(OpType::Rz, t, {0}); - for (unsigned i = 1; i != n_qubits; ++i) { - uncompute.add_op(OpType::CX, {i, 0}); + compute.add_op(OpType::CX, {qubits[i], qubits[0]}); } + action.add_op(OpType::Rz, paulis.coeff, {qubits[0]}); break; } case CXConfigType::Tree: { unsigned complete_layers = floor(log2(n_qubits)); unsigned dense_end = pow(2, complete_layers); for (unsigned i = 0; i < n_qubits - dense_end; i++) - compute.add_op( - OpType::CX, {dense_end + i, dense_end - 1 - i}); + compute.add_op( + OpType::CX, {qubits[dense_end + i], qubits[dense_end - 1 - i]}); for (unsigned step_size = 1; step_size < dense_end; step_size *= 2) { for (unsigned i = 0; i < dense_end; i += 2 * step_size) - compute.add_op(OpType::CX, {i + step_size, i}); - } - action.add_op(OpType::Rz, t, {0}); - for (unsigned step_size = dense_end / 2; step_size >= 1; step_size /= 2) { - for (unsigned i = 0; i < dense_end; i += 2 * step_size) - uncompute.add_op(OpType::CX, {i + step_size, i}); + compute.add_op(OpType::CX, {qubits[i + step_size], qubits[i]}); } - for (unsigned i = 0; i < n_qubits - dense_end; i++) - uncompute.add_op( - OpType::CX, {dense_end + i, dense_end - 1 - i}); + action.add_op(OpType::Rz, paulis.coeff, {qubits[0]}); break; } case CXConfigType::MultiQGate: { - std::vector> conjugations; int sign_correction = 1; for (int q = n_qubits - 1; q > 0; q -= 2) { if (q - 1 > 0) { @@ -331,89 +340,72 @@ Circuit phase_gadget(unsigned n_qubits, const Expr &t, CXConfigType cx_config) { // this is only equal to the CX decompositions above // up to phase, but phase differences are cancelled out by // its dagger XXPhase(-1/2) below. - compute.add_op(OpType::H, {i}); - compute.add_op(OpType::H, {j}); - compute.add_op(OpType::XXPhase3, 0.5, {i, j, 0}); + compute.add_op(OpType::H, {qubits[i]}); + compute.add_op(OpType::H, {qubits[j]}); + compute.add_op( + OpType::XXPhase3, 0.5, {qubits[i], qubits[j], qubits[0]}); sign_correction *= -1; - conjugations.push_back({i, j, 0}); } else { unsigned i = q; - compute.add_op(OpType::CX, {i, 0}); - conjugations.push_back({i, 0}); - } - } - action.add_op(OpType::Rz, sign_correction * t, {0}); - for (const auto &conj : conjugations) { - if (conj.size() == 2) { - uncompute.add_op(OpType::CX, conj); - } else { - TKET_ASSERT(conj.size() == 3); - uncompute.add_op(OpType::XXPhase3, -0.5, conj); - uncompute.add_op(OpType::H, {conj[0]}); - uncompute.add_op(OpType::H, {conj[1]}); + compute.add_op(OpType::CX, {qubits[i], qubits[0]}); } } + action.add_op( + OpType::Rz, sign_correction * paulis.coeff, {qubits[0]}); break; } } + // ConjugationBox components must be in the default register + compute.flatten_registers(); + action.flatten_registers(); ConjugationBox box( - std::make_shared(compute), std::make_shared(action), - std::make_shared(uncompute)); - new_circ.add_box(box, new_circ.all_qubits()); - return new_circ; + std::make_shared(compute), std::make_shared(action)); + circ.add_box(box, circ.all_qubits()); + return circ; } -Circuit pauli_gadget( - const std::vector &paulis, const Expr &t, CXConfigType cx_config) { - unsigned n = paulis.size(); - Circuit circ(n); - Circuit compute(n); - Circuit action(n); - Circuit uncompute(n); - std::vector qubits; - for (unsigned i = 0; i < n; i++) { - switch (paulis[i]) { - case Pauli::I: - break; - case Pauli::X: - compute.add_op(OpType::H, {i}); - qubits.push_back(i); - break; - case Pauli::Y: - compute.add_op(OpType::V, {i}); - qubits.push_back(i); - break; - case Pauli::Z: - qubits.push_back(i); - break; - } - } - if (qubits.empty()) { - circ.add_phase(-t / 2); +Circuit pauli_gadget_pair( + SpSymPauliTensor paulis0, SpSymPauliTensor paulis1, + CXConfigType cx_config) { + Circuit circ; + for (const std::pair &qp : paulis0.string) + circ.add_qubit(qp.first); + for (const std::pair &qp : paulis1.string) + circ.add_qubit(qp.first, false); + if (SpPauliString(paulis0.string) == SpPauliString{}) { + circ.append(pauli_gadget(paulis1, cx_config)); + circ.add_phase(-paulis0.coeff / 2); + return circ; + } else if (SpPauliString(paulis1.string) == SpPauliString{}) { + circ.append(pauli_gadget(paulis0, cx_config)); + circ.add_phase(-paulis1.coeff / 2); return circ; } - Vertex v = action.add_op(OpType::PhaseGadget, t, qubits); - Circuit cx_gadget = phase_gadget(action.n_in_edges(v), t, cx_config); - Subcircuit sub = {action.get_in_edges(v), action.get_all_out_edges(v), {v}}; - action.substitute(cx_gadget, sub, Circuit::VertexDeletion::Yes); - for (unsigned i = 0; i < n; i++) { - switch (paulis[i]) { - case Pauli::I: - break; - case Pauli::X: - uncompute.add_op(OpType::H, {i}); - break; - case Pauli::Y: - uncompute.add_op(OpType::Vdg, {i}); - break; - case Pauli::Z: - break; - } - } - ConjugationBox box( - std::make_shared(compute), std::make_shared(action), - std::make_shared(uncompute)); - circ.add_box(box, circ.all_qubits()); + paulis0.compress(); + paulis1.compress(); + Circuit u(circ); + Circuit v(circ); + + // Reduce the overlap down to at most 1 qubit, which may be matching or + // mismatching; allow both gadgets to build on that qubit + SpPauliStabiliser p0stab(paulis0.string); + SpPauliStabiliser p1stab(paulis1.string); + u.append(reduce_overlap_of_paulis(p0stab, p1stab, cx_config, true).first); + paulis0 = SpSymPauliTensor(p0stab) * SpSymPauliTensor({}, paulis0.coeff); + paulis1 = SpSymPauliTensor(p1stab) * SpSymPauliTensor({}, paulis1.coeff); + + /* + * Combine circuits to give final result + */ + v.append(pauli_gadget(paulis0)); + v.append(pauli_gadget(paulis1)); + // ConjugationBox components must be in the default register + qubit_vector_t all_qubits = u.all_qubits(); + u.flatten_registers(); + v.flatten_registers(); + ConjugationBox cjbox( + std::make_shared(u), std::make_shared(v)); + circ.add_box(cjbox, all_qubits); return circ; } diff --git a/tket/src/Circuit/PauliExpBoxes.cpp b/tket/src/Circuit/PauliExpBoxes.cpp index 9fcd33a06f..c219426c34 100644 --- a/tket/src/Circuit/PauliExpBoxes.cpp +++ b/tket/src/Circuit/PauliExpBoxes.cpp @@ -18,7 +18,6 @@ #include "tket/Circuit/CircUtils.hpp" #include "tket/Circuit/ConjugationBox.hpp" -#include "tket/Converters/PauliGadget.hpp" #include "tket/Converters/PhasePoly.hpp" #include "tket/Diagonalisation/Diagonalisation.hpp" #include "tket/Ops/OpJsonFactory.hpp" @@ -61,7 +60,11 @@ Op_ptr PauliExpBox::symbol_substitution( } void PauliExpBox::generate_circuit() const { - Circuit circ = pauli_gadget(paulis_.string, paulis_.coeff, cx_config_); + // paulis_ gets cast to a sparse form, so circuit from pauli_gadget will only + // contain qubits with {X, Y, Z}; appending it to a blank circuit containing + // all qubits makes the size of the circuit fixed + Circuit circ(paulis_.size()); + circ.append(pauli_gadget(paulis_, cx_config_)); circ_ = std::make_shared(circ); } @@ -149,8 +152,12 @@ Op_ptr PauliExpPairBox::symbol_substitution( } void PauliExpPairBox::generate_circuit() const { - Circuit circ = Circuit(paulis0_.size()); - append_pauli_gadget_pair(circ, paulis0_, paulis1_, cx_config_); + // paulis0_ and paulis1_ gets cast to a sparse form, so circuit from + // pauli_gadget_pair will only contain qubits with {X, Y, Z} on at least one; + // appending it to a blank circuit containing all qubits makes the size of the + // circuit fixed + Circuit circ(paulis0_.size()); + circ.append(pauli_gadget_pair(paulis0_, paulis1_, cx_config_)); circ_ = std::make_shared(circ); } @@ -306,7 +313,7 @@ void PauliExpCommutingSetBox::generate_circuit() const { Circuit phase_poly_circ(n_qubits); for (const SpSymPauliTensor &pgp : gadgets) { - append_single_pauli_gadget(phase_poly_circ, pgp); + phase_poly_circ.append(pauli_gadget(pgp, CXConfigType::Snake)); } phase_poly_circ.decompose_boxes_recursively(); PhasePolyBox ppbox(phase_poly_circ); @@ -638,4 +645,79 @@ void TermSequenceBox::generate_circuit() const { REGISTER_OPFACTORY(TermSequenceBox, TermSequenceBox) +void append_single_pauli_gadget_as_pauli_exp_box( + Circuit &circ, const SpSymPauliTensor &pauli, CXConfigType cx_config) { + std::vector string; + std::vector mapping; + for (const std::pair &term : pauli.string) { + string.push_back(term.second); + mapping.push_back(term.first); + } + PauliExpBox box(SymPauliTensor(string, pauli.coeff), cx_config); + circ.add_box(box, mapping); +} + +void append_pauli_gadget_pair_as_box( + Circuit &circ, const SpSymPauliTensor &pauli0, + const SpSymPauliTensor &pauli1, CXConfigType cx_config) { + std::vector mapping; + std::vector paulis0; + std::vector paulis1; + QubitPauliMap p1map = pauli1.string; + // add paulis for qubits in pauli0_string + for (const std::pair &term : pauli0.string) { + mapping.push_back(term.first); + paulis0.push_back(term.second); + auto found = p1map.find(term.first); + if (found == p1map.end()) { + paulis1.push_back(Pauli::I); + } else { + paulis1.push_back(found->second); + p1map.erase(found); + } + } + // add paulis for qubits in pauli1_string that weren't in pauli0_string + for (const std::pair &term : p1map) { + mapping.push_back(term.first); + paulis1.push_back(term.second); + paulis0.push_back(Pauli::I); // If pauli0_string contained qubit, would + // have been handled above + } + PauliExpPairBox box( + SymPauliTensor(paulis0, pauli0.coeff), + SymPauliTensor(paulis1, pauli1.coeff), cx_config); + circ.add_box(box, mapping); +} + +void append_commuting_pauli_gadget_set_as_box( + Circuit &circ, const std::list &gadgets, + CXConfigType cx_config) { + // Translate from QubitPauliTensors to vectors of Paulis of same length + // Preserves ordering of qubits + + std::set all_qubits; + for (const SpSymPauliTensor &gadget : gadgets) { + for (const std::pair &qubit_pauli : gadget.string) { + all_qubits.insert(qubit_pauli.first); + } + } + + std::vector mapping; + for (const auto &qubit : all_qubits) { + mapping.push_back(qubit); + } + + std::vector pauli_gadgets; + for (const SpSymPauliTensor &gadget : gadgets) { + SymPauliTensor &new_gadget = + pauli_gadgets.emplace_back(DensePauliMap{}, gadget.coeff); + for (const Qubit &qubit : mapping) { + new_gadget.string.push_back(gadget.get(qubit)); + } + } + + PauliExpCommutingSetBox box(pauli_gadgets, cx_config); + circ.add_box(box, mapping); +} + } // namespace tket diff --git a/tket/src/Clifford/ChoiMixTableau.cpp b/tket/src/Clifford/ChoiMixTableau.cpp index 77c8ef02a4..ae9698b9c2 100644 --- a/tket/src/Clifford/ChoiMixTableau.cpp +++ b/tket/src/Clifford/ChoiMixTableau.cpp @@ -102,12 +102,14 @@ ChoiMixTableau::ChoiMixTableau(const std::list& rows) MatrixXb zmat = MatrixXb::Zero(n_rows, n_qbs); VectorXb phase = VectorXb::Zero(n_rows); unsigned r = 0; + unsigned n_ys = 0; for (const row_tensor_t& row : rows) { for (const std::pair& qb : row.first.string) { unsigned c = col_index_.left.at(col_key_t{qb.first, TableauSegment::Input}); if (qb.second == Pauli::X || qb.second == Pauli::Y) xmat(r, c) = true; if (qb.second == Pauli::Z || qb.second == Pauli::Y) zmat(r, c) = true; + if (qb.second == Pauli::Y) ++n_ys; } for (const std::pair& qb : row.second.string) { unsigned c = @@ -115,7 +117,8 @@ ChoiMixTableau::ChoiMixTableau(const std::list& rows) if (qb.second == Pauli::X || qb.second == Pauli::Y) xmat(r, c) = true; if (qb.second == Pauli::Z || qb.second == Pauli::Y) zmat(r, c) = true; } - phase(r) = row.first.is_real_negative() ^ row.second.is_real_negative(); + phase(r) = row.first.is_real_negative() ^ row.second.is_real_negative() ^ + (n_ys % 2 == 1); ++r; } tab_ = SymplecticTableau(xmat, zmat, phase); @@ -125,22 +128,30 @@ unsigned ChoiMixTableau::get_n_rows() const { return tab_.get_n_rows(); } unsigned ChoiMixTableau::get_n_boundaries() const { return col_index_.size(); } -unsigned ChoiMixTableau::get_n_inputs() const { - unsigned n = 0; +unsigned ChoiMixTableau::get_n_inputs() const { return input_qubits().size(); } + +unsigned ChoiMixTableau::get_n_outputs() const { + return output_qubits().size(); +} + +qubit_vector_t ChoiMixTableau::input_qubits() const { + qubit_vector_t ins; BOOST_FOREACH ( tableau_col_index_t::left_const_reference entry, col_index_.left) { - if (entry.first.second == TableauSegment::Input) ++n; + if (entry.first.second == TableauSegment::Input) + ins.push_back(entry.first.first); } - return n; + return ins; } -unsigned ChoiMixTableau::get_n_outputs() const { - unsigned n = 0; +qubit_vector_t ChoiMixTableau::output_qubits() const { + qubit_vector_t outs; BOOST_FOREACH ( tableau_col_index_t::left_const_reference entry, col_index_.left) { - if (entry.first.second == TableauSegment::Output) ++n; + if (entry.first.second == TableauSegment::Output) + outs.push_back(entry.first.first); } - return n; + return outs; } ChoiMixTableau::row_tensor_t ChoiMixTableau::stab_to_row_tensor( @@ -173,17 +184,22 @@ PauliStabiliser ChoiMixTableau::row_tensor_to_stab( } ChoiMixTableau::row_tensor_t ChoiMixTableau::get_row(unsigned i) const { - return stab_to_row_tensor(tab_.get_pauli(i)); + ChoiMixTableau::row_tensor_t res = stab_to_row_tensor(tab_.get_pauli(i)); + res.first.transpose(); + res.second.coeff = (res.first.coeff + res.second.coeff) % 4; + res.first.coeff = 0; + return res; } ChoiMixTableau::row_tensor_t ChoiMixTableau::get_row_product( const std::vector& rows) const { row_tensor_t result = {{}, {}}; for (unsigned i : rows) { - row_tensor_t row_i = get_row(i); + row_tensor_t row_i = stab_to_row_tensor(tab_.get_pauli(i)); result.first = result.first * row_i.first; result.second = result.second * row_i.second; } + result.first.transpose(); result.second.coeff = (result.first.coeff + result.second.coeff) % 4; result.first.coeff = 0; return result; @@ -194,11 +210,26 @@ void ChoiMixTableau::apply_S(const Qubit& qb, TableauSegment seg) { tab_.apply_S(col); } +void ChoiMixTableau::apply_Z(const Qubit& qb, TableauSegment seg) { + unsigned col = col_index_.left.at(col_key_t{qb, seg}); + tab_.apply_Z(col); +} + void ChoiMixTableau::apply_V(const Qubit& qb, TableauSegment seg) { unsigned col = col_index_.left.at(col_key_t{qb, seg}); tab_.apply_V(col); } +void ChoiMixTableau::apply_X(const Qubit& qb, TableauSegment seg) { + unsigned col = col_index_.left.at(col_key_t{qb, seg}); + tab_.apply_X(col); +} + +void ChoiMixTableau::apply_H(const Qubit& qb, TableauSegment seg) { + unsigned col = col_index_.left.at(col_key_t{qb, seg}); + tab_.apply_H(col); +} + void ChoiMixTableau::apply_CX( const Qubit& control, const Qubit& target, TableauSegment seg) { unsigned uc = col_index_.left.at(col_key_t{control, seg}); @@ -210,20 +241,16 @@ void ChoiMixTableau::apply_gate( OpType type, const qubit_vector_t& qbs, TableauSegment seg) { switch (type) { case OpType::Z: { - apply_S(qbs.at(0), seg); - apply_S(qbs.at(0), seg); + apply_Z(qbs.at(0), seg); break; } case OpType::X: { - apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); + apply_X(qbs.at(0), seg); break; } case OpType::Y: { - apply_S(qbs.at(0), seg); - apply_S(qbs.at(0), seg); - apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); + apply_Z(qbs.at(0), seg); + apply_X(qbs.at(0), seg); break; } case OpType::S: { @@ -232,8 +259,7 @@ void ChoiMixTableau::apply_gate( } case OpType::Sdg: { apply_S(qbs.at(0), seg); - apply_S(qbs.at(0), seg); - apply_S(qbs.at(0), seg); + apply_Z(qbs.at(0), seg); break; } case OpType::SX: @@ -244,14 +270,11 @@ void ChoiMixTableau::apply_gate( case OpType::SXdg: case OpType::Vdg: { apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); + apply_X(qbs.at(0), seg); break; } case OpType::H: { - apply_S(qbs.at(0), seg); - apply_V(qbs.at(0), seg); - apply_S(qbs.at(0), seg); + apply_H(qbs.at(0), seg); break; } case OpType::CX: { @@ -263,25 +286,19 @@ void ChoiMixTableau::apply_gate( apply_S(qbs.at(1), seg); apply_CX(qbs.at(0), qbs.at(1), seg); apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); + apply_Z(qbs.at(1), seg); } else { apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); + apply_Z(qbs.at(1), seg); apply_CX(qbs.at(0), qbs.at(1), seg); apply_S(qbs.at(1), seg); } break; } case OpType::CZ: { - apply_S(qbs.at(1), seg); - apply_V(qbs.at(1), seg); - apply_S(qbs.at(1), seg); + apply_H(qbs.at(1), seg); apply_CX(qbs.at(0), qbs.at(1), seg); - apply_S(qbs.at(1), seg); - apply_V(qbs.at(1), seg); - apply_S(qbs.at(1), seg); + apply_H(qbs.at(1), seg); break; } case OpType::ZZMax: { @@ -292,16 +309,14 @@ void ChoiMixTableau::apply_gate( } case OpType::ECR: { if (seg == TableauSegment::Input) { - apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); + apply_X(qbs.at(0), seg); apply_S(qbs.at(0), seg); apply_V(qbs.at(1), seg); apply_CX(qbs.at(0), qbs.at(1), seg); } else { apply_CX(qbs.at(0), qbs.at(1), seg); apply_S(qbs.at(0), seg); - apply_V(qbs.at(0), seg); - apply_V(qbs.at(0), seg); + apply_X(qbs.at(0), seg); apply_V(qbs.at(1), seg); } break; @@ -312,8 +327,7 @@ void ChoiMixTableau::apply_gate( apply_CX(qbs.at(0), qbs.at(1), seg); apply_V(qbs.at(0), seg); apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); - apply_S(qbs.at(1), seg); + apply_Z(qbs.at(1), seg); apply_CX(qbs.at(0), qbs.at(1), seg); apply_V(qbs.at(0), seg); apply_V(qbs.at(1), seg); @@ -341,28 +355,25 @@ void ChoiMixTableau::apply_gate( // reinsert qubit initialised to maximally mixed state (no coherent // stabilizers) col_index_.insert({{qbs.at(0), TableauSegment::Input}, col}); - tab_.xmat_.conservativeResize(rows, col + 1); - tab_.xmat_.col(col) = MatrixXb::Zero(rows, 1); - tab_.zmat_.conservativeResize(rows, col + 1); - tab_.zmat_.col(col) = MatrixXb::Zero(rows, 1); - ++tab_.n_qubits_; + tab_.xmat.conservativeResize(rows, col + 1); + tab_.xmat.col(col) = MatrixXb::Zero(rows, 1); + tab_.zmat.conservativeResize(rows, col + 1); + tab_.zmat.col(col) = MatrixXb::Zero(rows, 1); } else { discard_qubit(qbs.at(0), TableauSegment::Output); unsigned col = get_n_boundaries(); unsigned rows = get_n_rows(); // reinsert qubit initialised to |0> (add a Z stabilizer) col_index_.insert({{qbs.at(0), TableauSegment::Output}, col}); - tab_.xmat_.conservativeResize(rows + 1, col + 1); - tab_.xmat_.col(col) = MatrixXb::Zero(rows + 1, 1); - tab_.xmat_.row(rows) = MatrixXb::Zero(1, col + 1); - tab_.zmat_.conservativeResize(rows + 1, col + 1); - tab_.zmat_.col(col) = MatrixXb::Zero(rows + 1, 1); - tab_.zmat_.row(rows) = MatrixXb::Zero(1, col + 1); - tab_.zmat_(rows, col) = true; - tab_.phase_.conservativeResize(rows + 1); - tab_.phase_(rows) = false; - ++tab_.n_rows_; - ++tab_.n_qubits_; + tab_.xmat.conservativeResize(rows + 1, col + 1); + tab_.xmat.col(col) = MatrixXb::Zero(rows + 1, 1); + tab_.xmat.row(rows) = MatrixXb::Zero(1, col + 1); + tab_.zmat.conservativeResize(rows + 1, col + 1); + tab_.zmat.col(col) = MatrixXb::Zero(rows + 1, 1); + tab_.zmat.row(rows) = MatrixXb::Zero(1, col + 1); + tab_.zmat(rows, col) = true; + tab_.phase.conservativeResize(rows + 1); + tab_.phase(rows) = false; } break; } @@ -399,10 +410,10 @@ void ChoiMixTableau::post_select(const Qubit& qb, TableauSegment seg) { unsigned n_cols = get_n_boundaries(); unsigned col = col_index_.left.at(col_key_t{qb, seg}); for (unsigned r = 0; r < n_rows; ++r) { - if (tab_.zmat_(r, col)) { + if (tab_.zmat(r, col)) { bool only_z = true; for (unsigned c = 0; c < n_cols; ++c) { - if ((tab_.xmat_(r, c) || tab_.zmat_(r, c)) && (c != col)) { + if ((tab_.xmat(r, c) || tab_.zmat(r, c)) && (c != col)) { only_z = false; break; } @@ -410,7 +421,7 @@ void ChoiMixTableau::post_select(const Qubit& qb, TableauSegment seg) { if (!only_z) break; // Not deterministic // From here, we know we are in a deterministic case // If deterministically fail, throw an exception - if (tab_.phase_(r)) + if (tab_.phase(r)) throw std::logic_error( "Post-selecting a tableau fails deterministically"); // Otherwise, we succeed and remove the stabilizer @@ -423,7 +434,7 @@ void ChoiMixTableau::post_select(const Qubit& qb, TableauSegment seg) { // Isolate a single row with an X (if one exists) std::optional x_row = std::nullopt; for (unsigned r = 0; r < n_rows; ++r) { - if (tab_.xmat_(r, col)) { + if (tab_.xmat(r, col)) { if (x_row) { // Already found another row with an X, so combine them tab_.row_mult(*x_row, r); @@ -446,7 +457,7 @@ void ChoiMixTableau::discard_qubit(const Qubit& qb, TableauSegment seg) { // Isolate a single row with an X (if one exists) std::optional x_row = std::nullopt; for (unsigned r = 0; r < get_n_rows(); ++r) { - if (tab_.xmat_(r, col)) { + if (tab_.xmat(r, col)) { if (x_row) { // Already found another row with an X, so combine them tab_.row_mult(*x_row, r); @@ -464,7 +475,7 @@ void ChoiMixTableau::discard_qubit(const Qubit& qb, TableauSegment seg) { // Isolate a single row with a Z (if one exists) std::optional z_row = std::nullopt; for (unsigned r = 0; r < get_n_rows(); ++r) { - if (tab_.zmat_(r, col)) { + if (tab_.zmat(r, col)) { if (z_row) { // Already found another row with a Z, so combine them tab_.row_mult(*z_row, r); @@ -487,7 +498,7 @@ void ChoiMixTableau::collapse_qubit(const Qubit& qb, TableauSegment seg) { // Isolate a single row with an X (if one exists) std::optional x_row = std::nullopt; for (unsigned r = 0; r < get_n_rows(); ++r) { - if (tab_.xmat_(r, col)) { + if (tab_.xmat(r, col)) { if (x_row) { // Already found another row with an X, so combine them tab_.row_mult(*x_row, r); @@ -512,14 +523,13 @@ void ChoiMixTableau::remove_row(unsigned row) { unsigned n_rows = get_n_rows(); unsigned n_cols = get_n_boundaries(); if (row < n_rows - 1) { - tab_.xmat_.row(row) = tab_.xmat_.row(n_rows - 1); - tab_.zmat_.row(row) = tab_.zmat_.row(n_rows - 1); - tab_.phase_(row) = tab_.phase_(n_rows - 1); + tab_.xmat.row(row) = tab_.xmat.row(n_rows - 1); + tab_.zmat.row(row) = tab_.zmat.row(n_rows - 1); + tab_.phase(row) = tab_.phase(n_rows - 1); } - tab_.xmat_.conservativeResize(n_rows - 1, n_cols); - tab_.zmat_.conservativeResize(n_rows - 1, n_cols); - tab_.phase_.conservativeResize(n_rows - 1); - --tab_.n_rows_; + tab_.xmat.conservativeResize(n_rows - 1, n_cols); + tab_.zmat.conservativeResize(n_rows - 1, n_cols); + tab_.phase.conservativeResize(n_rows - 1); } void ChoiMixTableau::remove_col(unsigned col) { @@ -530,11 +540,11 @@ void ChoiMixTableau::remove_col(unsigned col) { unsigned n_rows = get_n_rows(); unsigned n_cols = get_n_boundaries(); if (col < n_cols - 1) { - tab_.xmat_.col(col) = tab_.xmat_.col(n_cols - 1); - tab_.zmat_.col(col) = tab_.zmat_.col(n_cols - 1); + tab_.xmat.col(col) = tab_.xmat.col(n_cols - 1); + tab_.zmat.col(col) = tab_.zmat.col(n_cols - 1); } - tab_.xmat_.conservativeResize(n_rows, n_cols - 1); - tab_.zmat_.conservativeResize(n_rows, n_cols - 1); + tab_.xmat.conservativeResize(n_rows, n_cols - 1); + tab_.zmat.conservativeResize(n_rows, n_cols - 1); col_index_.right.erase(col); if (col < n_cols - 1) { tableau_col_index_t::right_iterator it = col_index_.right.find(n_cols - 1); @@ -542,7 +552,6 @@ void ChoiMixTableau::remove_col(unsigned col) { col_index_.right.erase(it); col_index_.insert({last, col}); } - --tab_.n_qubits_; } void ChoiMixTableau::canonical_column_order(TableauSegment first) { @@ -579,10 +588,10 @@ void ChoiMixTableau::canonical_column_order(TableauSegment first) { for (unsigned j = 0; j < i; ++j) { col_key_t key = new_index.right.at(j); unsigned c = col_index_.left.at(key); - xmat.col(j) = tab_.xmat_.col(c); - zmat.col(j) = tab_.zmat_.col(c); + xmat.col(j) = tab_.xmat.col(c); + zmat.col(j) = tab_.zmat.col(c); } - tab_ = SymplecticTableau(xmat, zmat, tab_.phase_); + tab_ = SymplecticTableau(xmat, zmat, tab_.phase); col_index_ = new_index; } @@ -620,12 +629,12 @@ ChoiMixTableau ChoiMixTableau::compose( } MatrixXb fullx(f_rows + s_rows, f_cols + s_cols), fullz(f_rows + s_rows, f_cols + s_cols); - fullx << first.tab_.xmat_, MatrixXb::Zero(f_rows, s_cols), - MatrixXb::Zero(s_rows, f_cols), second.tab_.xmat_; - fullz << first.tab_.zmat_, MatrixXb::Zero(f_rows, s_cols), - MatrixXb::Zero(s_rows, f_cols), second.tab_.zmat_; + fullx << first.tab_.xmat, MatrixXb::Zero(f_rows, s_cols), + MatrixXb::Zero(s_rows, f_cols), second.tab_.xmat; + fullz << first.tab_.zmat, MatrixXb::Zero(f_rows, s_cols), + MatrixXb::Zero(s_rows, f_cols), second.tab_.zmat; VectorXb fullph(f_rows + s_rows); - fullph << first.tab_.phase_, second.tab_.phase_; + fullph << first.tab_.phase, second.tab_.phase; ChoiMixTableau combined(fullx, fullz, fullph, 0); // For each connecting pair of qubits, compose via a Bell post-selection for (unsigned i = 0; i < f_cols; ++i) { diff --git a/tket/src/Clifford/SymplecticTableau.cpp b/tket/src/Clifford/SymplecticTableau.cpp index 5c2cd7513a..4c15b44d1f 100644 --- a/tket/src/Clifford/SymplecticTableau.cpp +++ b/tket/src/Clifford/SymplecticTableau.cpp @@ -62,56 +62,51 @@ const std::map, std::pair> SymplecticTableau::SymplecticTableau( const MatrixXb &xmat, const MatrixXb &zmat, const VectorXb &phase) - : n_rows_(xmat.rows()), - n_qubits_(xmat.cols()), - xmat_(xmat), - zmat_(zmat), - phase_(phase) { - if (zmat.rows() != n_rows_ || phase_.size() != n_rows_) + : xmat(xmat), zmat(zmat), phase(phase) { + if (zmat.rows() != xmat.rows() || phase.size() != xmat.rows()) throw std::invalid_argument( "Tableau must have the same number of rows in each component."); - if (zmat.cols() != n_qubits_) + if (zmat.cols() != xmat.cols()) throw std::invalid_argument( "Tableau must have the same number of columns in x and z components."); } SymplecticTableau::SymplecticTableau(const PauliStabiliserVec &rows) { - n_rows_ = rows.size(); - if (n_rows_ == 0) - n_qubits_ = 0; - else - n_qubits_ = rows[0].string.size(); - xmat_ = MatrixXb::Zero(n_rows_, n_qubits_); - zmat_ = MatrixXb::Zero(n_rows_, n_qubits_); - phase_ = VectorXb::Zero(n_rows_); - for (unsigned i = 0; i < n_rows_; ++i) { + unsigned n_rows = rows.size(); + unsigned n_qubits = 0; + if (n_rows != 0) n_qubits = rows[0].string.size(); + xmat = MatrixXb::Zero(n_rows, n_qubits); + zmat = MatrixXb::Zero(n_rows, n_qubits); + phase = VectorXb::Zero(n_rows); + for (unsigned i = 0; i < n_rows; ++i) { const PauliStabiliser &stab = rows[i]; - if (stab.string.size() != n_qubits_) + if (stab.string.size() != n_qubits) throw std::invalid_argument( "Tableau must have the same number of qubits in each row."); - for (unsigned q = 0; q < n_qubits_; ++q) { + for (unsigned q = 0; q < n_qubits; ++q) { const Pauli &p = stab.get(q); - xmat_(i, q) = (p == Pauli::X) || (p == Pauli::Y); - zmat_(i, q) = (p == Pauli::Z) || (p == Pauli::Y); + xmat(i, q) = (p == Pauli::X) || (p == Pauli::Y); + zmat(i, q) = (p == Pauli::Z) || (p == Pauli::Y); } - phase_(i) = stab.is_real_negative(); + phase(i) = stab.is_real_negative(); } } -unsigned SymplecticTableau::get_n_rows() const { return n_rows_; } -unsigned SymplecticTableau::get_n_qubits() const { return n_qubits_; } +unsigned SymplecticTableau::get_n_rows() const { return xmat.rows(); } +unsigned SymplecticTableau::get_n_qubits() const { return xmat.cols(); } PauliStabiliser SymplecticTableau::get_pauli(unsigned i) const { - std::vector str(n_qubits_); - for (unsigned q = 0; q < n_qubits_; ++q) { - str[q] = BoolPauli{xmat_(i, q), zmat_(i, q)}.to_pauli(); + unsigned n_qubits = get_n_qubits(); + std::vector str(n_qubits); + for (unsigned q = 0; q < n_qubits; ++q) { + str[q] = BoolPauli{xmat(i, q), zmat(i, q)}.to_pauli(); } - return PauliStabiliser(str, phase_(i) ? 2 : 0); + return PauliStabiliser(str, phase(i) ? 2 : 0); } std::ostream &operator<<(std::ostream &os, const SymplecticTableau &tab) { - for (unsigned i = 0; i < tab.n_rows_; ++i) { - os << tab.xmat_.row(i) << " " << tab.zmat_.row(i) << " " << tab.phase_(i) + for (unsigned i = 0; i < tab.get_n_rows(); ++i) { + os << tab.xmat.row(i) << " " << tab.zmat.row(i) << " " << tab.phase(i) << std::endl; } return os; @@ -120,40 +115,58 @@ std::ostream &operator<<(std::ostream &os, const SymplecticTableau &tab) { bool SymplecticTableau::operator==(const SymplecticTableau &other) const { // Need this to short-circuit before matrix checks as comparing matrices of // different sizes will throw an exception - return (this->n_rows_ == other.n_rows_) && - (this->n_qubits_ == other.n_qubits_) && (this->xmat_ == other.xmat_) && - (this->zmat_ == other.zmat_) && (this->phase_ == other.phase_); + return (this->get_n_rows() == other.get_n_rows()) && + (this->get_n_qubits() == other.get_n_qubits()) && + (this->xmat == other.xmat) && (this->zmat == other.zmat) && + (this->phase == other.phase); } void SymplecticTableau::row_mult(unsigned ra, unsigned rw, Complex coeff) { - MatrixXb::RowXpr xa = xmat_.row(ra); - MatrixXb::RowXpr za = zmat_.row(ra); - MatrixXb::RowXpr xw = xmat_.row(rw); - MatrixXb::RowXpr zw = zmat_.row(rw); - row_mult(xa, za, phase_(ra), xw, zw, phase_(rw), coeff, xw, zw, phase_(rw)); + MatrixXb::RowXpr xa = xmat.row(ra); + MatrixXb::RowXpr za = zmat.row(ra); + MatrixXb::RowXpr xw = xmat.row(rw); + MatrixXb::RowXpr zw = zmat.row(rw); + row_mult(xa, za, phase(ra), xw, zw, phase(rw), coeff, xw, zw, phase(rw)); } void SymplecticTableau::apply_S(unsigned qb) { - MatrixXb::ColXpr xcol = xmat_.col(qb); - MatrixXb::ColXpr zcol = zmat_.col(qb); - col_mult(xcol, zcol, false, zcol, phase_); + MatrixXb::ColXpr xcol = xmat.col(qb); + MatrixXb::ColXpr zcol = zmat.col(qb); + col_mult(xcol, zcol, false, zcol, phase); +} + +void SymplecticTableau::apply_Z(unsigned qb) { + for (unsigned i = 0; i < get_n_rows(); ++i) phase(i) = phase(i) ^ xmat(i, qb); } void SymplecticTableau::apply_V(unsigned qb) { - MatrixXb::ColXpr xcol = xmat_.col(qb); - MatrixXb::ColXpr zcol = zmat_.col(qb); - col_mult(zcol, xcol, true, xcol, phase_); + MatrixXb::ColXpr xcol = xmat.col(qb); + MatrixXb::ColXpr zcol = zmat.col(qb); + col_mult(zcol, xcol, true, xcol, phase); +} + +void SymplecticTableau::apply_X(unsigned qb) { + for (unsigned i = 0; i < get_n_rows(); ++i) phase(i) = phase(i) ^ zmat(i, qb); +} + +void SymplecticTableau::apply_H(unsigned qb) { + for (unsigned i = 0; i < get_n_rows(); ++i) { + phase(i) = phase(i) ^ (xmat(i, qb) && zmat(i, qb)); + bool temp = xmat(i, qb); + xmat(i, qb) = zmat(i, qb); + zmat(i, qb) = temp; + } } void SymplecticTableau::apply_CX(unsigned qc, unsigned qt) { if (qc == qt) throw std::logic_error( "Attempting to apply a CX with equal control and target in a tableau"); - for (unsigned i = 0; i < n_rows_; ++i) { - phase_(i) = phase_(i) ^ (xmat_(i, qc) && zmat_(i, qt) && - !(xmat_(i, qt) ^ zmat_(i, qc))); - xmat_(i, qt) = xmat_(i, qc) ^ xmat_(i, qt); - zmat_(i, qc) = zmat_(i, qc) ^ zmat_(i, qt); + for (unsigned i = 0; i < get_n_rows(); ++i) { + phase(i) = + phase(i) ^ (xmat(i, qc) && zmat(i, qt) && !(xmat(i, qt) ^ zmat(i, qc))); + xmat(i, qt) = xmat(i, qc) ^ xmat(i, qt); + zmat(i, qc) = zmat(i, qc) ^ zmat(i, qt); } } @@ -161,20 +174,16 @@ void SymplecticTableau::apply_gate( OpType type, const std::vector &qbs) { switch (type) { case OpType::Z: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); + apply_Z(qbs.at(0)); break; } case OpType::X: { - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); + apply_X(qbs.at(0)); break; } case OpType::Y: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); + apply_Z(qbs.at(0)); + apply_X(qbs.at(0)); break; } case OpType::S: { @@ -183,8 +192,7 @@ void SymplecticTableau::apply_gate( } case OpType::Sdg: { apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); + apply_Z(qbs.at(0)); break; } case OpType::V: @@ -195,14 +203,11 @@ void SymplecticTableau::apply_gate( case OpType::Vdg: case OpType::SXdg: { apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); + apply_X(qbs.at(0)); break; } case OpType::H: { - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); + apply_H(qbs.at(0)); break; } case OpType::CX: { @@ -211,63 +216,43 @@ void SymplecticTableau::apply_gate( } case OpType::CY: { apply_S(qbs.at(1)); - apply_S(qbs.at(1)); - apply_S(qbs.at(1)); + apply_Z(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); apply_S(qbs.at(1)); break; } case OpType::CZ: { - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); + apply_H(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); + apply_H(qbs.at(1)); break; } case OpType::ZZMax: { + apply_H(qbs.at(1)); apply_S(qbs.at(0)); - apply_S(qbs.at(1)); - apply_S(qbs.at(1)); - apply_S(qbs.at(1)); apply_V(qbs.at(1)); - apply_S(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); + apply_H(qbs.at(1)); break; } case OpType::ECR: { apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - apply_V(qbs.at(1)); - apply_V(qbs.at(1)); + apply_X(qbs.at(0)); apply_V(qbs.at(1)); + apply_X(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); break; } case OpType::ISWAPMax: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_S(qbs.at(1)); - apply_S(qbs.at(1)); apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(1)); - apply_V(qbs.at(1)); apply_V(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); apply_V(qbs.at(0)); apply_S(qbs.at(1)); - apply_V(qbs.at(1)); + apply_Z(qbs.at(1)); apply_CX(qbs.at(0), qbs.at(1)); apply_V(qbs.at(0)); + apply_V(qbs.at(1)); break; } case OpType::SWAP: { @@ -294,7 +279,8 @@ void SymplecticTableau::apply_gate( void SymplecticTableau::apply_pauli_gadget( const PauliStabiliser &pauli, unsigned half_pis) { - if (pauli.string.size() != n_qubits_) { + unsigned n_qubits = get_n_qubits(); + if (pauli.string.size() != n_qubits) { throw std::invalid_argument( "Cannot apply pauli gadget to SymplecticTableau; string and tableau " "have different numbers of qubits"); @@ -326,39 +312,40 @@ void SymplecticTableau::apply_pauli_gadget( // From here, half_pis == 1 or 3 // They act the same except for a phase flip on the product term - MatrixXb pauli_xrow = MatrixXb::Zero(1, n_qubits_); - MatrixXb pauli_zrow = MatrixXb::Zero(1, n_qubits_); - for (unsigned i = 0; i < n_qubits_; ++i) { + MatrixXb pauli_xrow = MatrixXb::Zero(1, n_qubits); + MatrixXb pauli_zrow = MatrixXb::Zero(1, n_qubits); + for (unsigned i = 0; i < n_qubits; ++i) { Pauli p = pauli.string.at(i); pauli_xrow(i) = (p == Pauli::X) || (p == Pauli::Y); pauli_zrow(i) = (p == Pauli::Z) || (p == Pauli::Y); } - bool phase = pauli.is_real_negative() ^ (half_pis == 3); + bool phase_flip = pauli.is_real_negative() ^ (half_pis == 3); - for (unsigned i = 0; i < n_rows_; ++i) { + for (unsigned i = 0; i < get_n_rows(); ++i) { bool anti = false; - MatrixXb::RowXpr xr = xmat_.row(i); - MatrixXb::RowXpr zr = zmat_.row(i); - for (unsigned q = 0; q < n_qubits_; ++q) { + MatrixXb::RowXpr xr = xmat.row(i); + MatrixXb::RowXpr zr = zmat.row(i); + for (unsigned q = 0; q < n_qubits; ++q) { anti ^= (xr(q) && pauli_zrow(q)); anti ^= (zr(q) && pauli_xrow(q)); } if (anti) { row_mult( - xr, zr, phase_(i), pauli_xrow.row(0), pauli_zrow.row(0), phase, i_, - xr, zr, phase_(i)); + xr, zr, phase(i), pauli_xrow.row(0), pauli_zrow.row(0), phase_flip, + i_, xr, zr, phase(i)); } } } MatrixXb SymplecticTableau::anticommuting_rows() const { - MatrixXb res = MatrixXb::Zero(n_rows_, n_rows_); - for (unsigned i = 0; i < n_rows_; ++i) { + unsigned n_rows = get_n_rows(); + MatrixXb res = MatrixXb::Zero(n_rows, n_rows); + for (unsigned i = 0; i < n_rows; ++i) { for (unsigned j = 0; j < i; ++j) { bool anti = false; - for (unsigned q = 0; q < n_qubits_; ++q) { - anti ^= (xmat_(i, q) && zmat_(j, q)); - anti ^= (xmat_(j, q) && zmat_(i, q)); + for (unsigned q = 0; q < get_n_qubits(); ++q) { + anti ^= (xmat(i, q) && zmat(j, q)); + anti ^= (xmat(j, q) && zmat(i, q)); } res(i, j) = anti; res(j, i) = anti; @@ -372,32 +359,33 @@ unsigned SymplecticTableau::rank() const { SymplecticTableau copy(*this); copy.gaussian_form(); unsigned empty_rows = 0; - for (unsigned i = 0; i < n_rows_; ++i) { - if (copy.xmat_.row(n_rows_ - 1 - i).isZero() && - copy.zmat_.row(n_rows_ - 1 - i).isZero()) + unsigned n_rows = get_n_rows(); + for (unsigned i = 0; i < n_rows; ++i) { + if (copy.xmat.row(n_rows - 1 - i).isZero() && + copy.zmat.row(n_rows - 1 - i).isZero()) ++empty_rows; else break; } - return n_rows_ - empty_rows; + return n_rows - empty_rows; } SymplecticTableau SymplecticTableau::conjugate() const { SymplecticTableau conj(*this); - for (unsigned i = 0; i < n_rows_; ++i) { + for (unsigned i = 0; i < get_n_rows(); ++i) { unsigned sum = 0; - for (unsigned j = 0; j < n_qubits_; ++j) { - if (xmat_(i, j) && zmat_(i, j)) ++sum; + for (unsigned j = 0; j < get_n_qubits(); ++j) { + if (xmat(i, j) && zmat(i, j)) ++sum; } - if (sum % 2 == 1) conj.phase_(i) ^= true; + if (sum % 2 == 1) conj.phase(i) ^= true; } return conj; } void SymplecticTableau::gaussian_form() { - MatrixXb fullmat = MatrixXb::Zero(n_rows_, 2 * n_qubits_); - fullmat(Eigen::all, Eigen::seq(0, Eigen::last, 2)) = xmat_; - fullmat(Eigen::all, Eigen::seq(1, Eigen::last, 2)) = zmat_; + MatrixXb fullmat = MatrixXb::Zero(get_n_rows(), 2 * get_n_qubits()); + fullmat(Eigen::all, Eigen::seq(0, Eigen::last, 2)) = xmat; + fullmat(Eigen::all, Eigen::seq(1, Eigen::last, 2)) = zmat; std::vector> row_ops = gaussian_elimination_row_ops(fullmat); for (const std::pair &op : row_ops) { @@ -411,7 +399,7 @@ void SymplecticTableau::row_mult( Complex phase, MatrixXb::RowXpr &xw, MatrixXb::RowXpr &zw, bool &pw) { if (pa) phase *= -1; if (pb) phase *= -1; - for (unsigned i = 0; i < n_qubits_; i++) { + for (unsigned i = 0; i < get_n_qubits(); i++) { std::pair res = BoolPauli::mult_lut.at({{xa(i), za(i)}, {xb(i), zb(i)}}); xw(i) = res.first.x; @@ -424,18 +412,18 @@ void SymplecticTableau::row_mult( void SymplecticTableau::col_mult( const MatrixXb::ColXpr &a, const MatrixXb::ColXpr &b, bool flip, MatrixXb::ColXpr &w, VectorXb &pw) { - for (unsigned i = 0; i < n_rows_; i++) { + for (unsigned i = 0; i < get_n_rows(); i++) { pw(i) = pw(i) ^ (a(i) && (b(i) ^ flip)); w(i) = a(i) ^ b(i); } } void to_json(nlohmann::json &j, const SymplecticTableau &tab) { - j["nrows"] = tab.n_rows_; - j["nqubits"] = tab.n_qubits_; - j["xmat"] = tab.xmat_; - j["zmat"] = tab.zmat_; - j["phase"] = tab.phase_; + j["nrows"] = tab.get_n_rows(); + j["nqubits"] = tab.get_n_qubits(); + j["xmat"] = tab.xmat; + j["zmat"] = tab.zmat; + j["phase"] = tab.phase; } void from_json(const nlohmann::json &j, SymplecticTableau &tab) { diff --git a/tket/src/Clifford/UnitaryTableau.cpp b/tket/src/Clifford/UnitaryTableau.cpp index 3de53c2777..e73a1d1203 100644 --- a/tket/src/Clifford/UnitaryTableau.cpp +++ b/tket/src/Clifford/UnitaryTableau.cpp @@ -155,6 +155,16 @@ void UnitaryTableau::apply_S_at_end(const Qubit& qb) { tab_.apply_S(uqb); } +void UnitaryTableau::apply_Z_at_front(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + tab_.phase(uqb) = !tab_.phase(uqb); +} + +void UnitaryTableau::apply_Z_at_end(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + tab_.apply_Z(uqb); +} + void UnitaryTableau::apply_V_at_front(const Qubit& qb) { unsigned uqb = qubits_.left.at(qb); tab_.row_mult(uqb, uqb + qubits_.size(), -i_); @@ -165,6 +175,31 @@ void UnitaryTableau::apply_V_at_end(const Qubit& qb) { tab_.apply_V(uqb); } +void UnitaryTableau::apply_X_at_front(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + tab_.phase(uqb + qubits_.size()) = !tab_.phase(uqb + qubits_.size()); +} + +void UnitaryTableau::apply_X_at_end(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + tab_.apply_X(uqb); +} + +void UnitaryTableau::apply_H_at_front(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + unsigned n_qubits = qubits_.size(); + bool temp = tab_.phase(uqb); + tab_.phase(uqb) = tab_.phase(uqb + n_qubits); + tab_.phase(uqb + n_qubits) = temp; + tab_.xmat.row(uqb).swap(tab_.xmat.row(uqb + n_qubits)); + tab_.zmat.row(uqb).swap(tab_.zmat.row(uqb + n_qubits)); +} + +void UnitaryTableau::apply_H_at_end(const Qubit& qb) { + unsigned uqb = qubits_.left.at(qb); + tab_.apply_H(uqb); +} + void UnitaryTableau::apply_CX_at_front( const Qubit& control, const Qubit& target) { unsigned uc = qubits_.left.at(control); @@ -184,20 +219,16 @@ void UnitaryTableau::apply_gate_at_front( OpType type, const qubit_vector_t& qbs) { switch (type) { case OpType::Z: { - apply_S_at_front(qbs.at(0)); - apply_S_at_front(qbs.at(0)); + apply_Z_at_front(qbs.at(0)); break; } case OpType::X: { - apply_V_at_front(qbs.at(0)); - apply_V_at_front(qbs.at(0)); + apply_X_at_front(qbs.at(0)); break; } case OpType::Y: { - apply_S_at_front(qbs.at(0)); - apply_S_at_front(qbs.at(0)); - apply_V_at_front(qbs.at(0)); - apply_V_at_front(qbs.at(0)); + apply_Z_at_front(qbs.at(0)); + apply_X_at_front(qbs.at(0)); break; } case OpType::S: { @@ -206,24 +237,22 @@ void UnitaryTableau::apply_gate_at_front( } case OpType::Sdg: { apply_S_at_front(qbs.at(0)); - apply_S_at_front(qbs.at(0)); - apply_S_at_front(qbs.at(0)); + apply_Z_at_front(qbs.at(0)); break; } - case OpType::V: { + case OpType::V: + case OpType::SX: { apply_V_at_front(qbs.at(0)); break; } - case OpType::Vdg: { - apply_V_at_front(qbs.at(0)); - apply_V_at_front(qbs.at(0)); + case OpType::Vdg: + case OpType::SXdg: { apply_V_at_front(qbs.at(0)); + apply_X_at_front(qbs.at(0)); break; } case OpType::H: { - apply_S_at_front(qbs.at(0)); - apply_V_at_front(qbs.at(0)); - apply_S_at_front(qbs.at(0)); + apply_H_at_front(qbs.at(0)); break; } case OpType::CX: { @@ -234,18 +263,41 @@ void UnitaryTableau::apply_gate_at_front( apply_S_at_front(qbs.at(1)); apply_CX_at_front(qbs.at(0), qbs.at(1)); apply_S_at_front(qbs.at(1)); - apply_S_at_front(qbs.at(1)); - apply_S_at_front(qbs.at(1)); + apply_Z_at_front(qbs.at(1)); break; } case OpType::CZ: { - apply_S_at_front(qbs.at(1)); + apply_H_at_front(qbs.at(1)); + apply_CX_at_front(qbs.at(0), qbs.at(1)); + apply_H_at_front(qbs.at(1)); + break; + } + case OpType::ZZMax: { + apply_H_at_front(qbs.at(1)); + apply_S_at_front(qbs.at(0)); apply_V_at_front(qbs.at(1)); - apply_S_at_front(qbs.at(1)); apply_CX_at_front(qbs.at(0), qbs.at(1)); - apply_S_at_front(qbs.at(1)); + apply_H_at_front(qbs.at(1)); + break; + } + case OpType::ECR: { + apply_CX_at_front(qbs.at(0), qbs.at(1)); + apply_X_at_front(qbs.at(0)); + apply_S_at_front(qbs.at(0)); apply_V_at_front(qbs.at(1)); + apply_X_at_front(qbs.at(1)); + break; + } + case OpType::ISWAPMax: { + apply_V_at_front(qbs.at(0)); + apply_V_at_front(qbs.at(1)); + apply_CX_at_front(qbs.at(0), qbs.at(1)); + apply_V_at_front(qbs.at(0)); apply_S_at_front(qbs.at(1)); + apply_Z_at_front(qbs.at(1)); + apply_CX_at_front(qbs.at(0), qbs.at(1)); + apply_V_at_front(qbs.at(0)); + apply_V_at_front(qbs.at(1)); break; } case OpType::SWAP: { @@ -383,8 +435,8 @@ UnitaryTableau UnitaryTableau::dagger() const { for (unsigned j = 0; j < nqb; ++j) { // Take effect of some input on some output and invert auto inv_cell = invert_cell_map().at( - {BoolPauli{tab_.xmat_(i, j), tab_.zmat_(i, j)}, - BoolPauli{tab_.xmat_(i + nqb, j), tab_.zmat_(i + nqb, j)}}); + {BoolPauli{tab_.xmat(i, j), tab_.zmat(i, j)}, + BoolPauli{tab_.xmat(i + nqb, j), tab_.zmat(i + nqb, j)}}); // Transpose tableau and fill in cell dxx(j, i) = inv_cell.first.x; dxz(j, i) = inv_cell.first.z; @@ -399,9 +451,9 @@ UnitaryTableau UnitaryTableau::dagger() const { // Correct phases for (unsigned i = 0; i < nqb; ++i) { SpPauliStabiliser xr = dag.get_xrow(qubits_.right.at(i)); - dag.tab_.phase_(i) = get_row_product(xr).is_real_negative(); + dag.tab_.phase(i) = get_row_product(xr).is_real_negative(); SpPauliStabiliser zr = dag.get_zrow(qubits_.right.at(i)); - dag.tab_.phase_(i + nqb) = get_row_product(zr).is_real_negative(); + dag.tab_.phase(i + nqb) = get_row_product(zr).is_real_negative(); } return dag; @@ -422,14 +474,14 @@ std::ostream& operator<<(std::ostream& os, const UnitaryTableau& tab) { unsigned nqs = tab.qubits_.size(); for (unsigned i = 0; i < nqs; ++i) { Qubit qi = tab.qubits_.right.at(i); - os << "X@" << qi.repr() << "\t->\t" << tab.tab_.xmat_.row(i) << " " - << tab.tab_.zmat_.row(i) << " " << tab.tab_.phase_(i) << std::endl; + os << "X@" << qi.repr() << "\t->\t" << tab.tab_.xmat.row(i) << " " + << tab.tab_.zmat.row(i) << " " << tab.tab_.phase(i) << std::endl; } os << "--" << std::endl; for (unsigned i = 0; i < nqs; ++i) { Qubit qi = tab.qubits_.right.at(i); - os << "Z@" << qi.repr() << "\t->\t" << tab.tab_.xmat_.row(i + nqs) << " " - << tab.tab_.zmat_.row(i + nqs) << " " << tab.tab_.phase_(i + nqs) + os << "Z@" << qi.repr() << "\t->\t" << tab.tab_.xmat.row(i + nqs) << " " + << tab.tab_.zmat.row(i + nqs) << " " << tab.tab_.phase(i + nqs) << std::endl; } return os; @@ -446,13 +498,13 @@ bool UnitaryTableau::operator==(const UnitaryTableau& other) const { for (unsigned j = 0; j < nq; ++j) { Qubit qj = qubits_.right.at(j); unsigned oj = other.qubits_.left.at(qj); - if (tab_.xmat_(i, j) != other.tab_.xmat_(oi, oj)) return false; - if (tab_.zmat_(i, j) != other.tab_.zmat_(oi, oj)) return false; - if (tab_.xmat_(i + nq, j) != other.tab_.xmat_(oi + nq, oj)) return false; - if (tab_.zmat_(i + nq, j) != other.tab_.zmat_(oi + nq, oj)) return false; + if (tab_.xmat(i, j) != other.tab_.xmat(oi, oj)) return false; + if (tab_.zmat(i, j) != other.tab_.zmat(oi, oj)) return false; + if (tab_.xmat(i + nq, j) != other.tab_.xmat(oi + nq, oj)) return false; + if (tab_.zmat(i + nq, j) != other.tab_.zmat(oi + nq, oj)) return false; } - if (tab_.phase_(i) != other.tab_.phase_(oi)) return false; - if (tab_.phase_(i + nq) != other.tab_.phase_(oi + nq)) return false; + if (tab_.phase(i) != other.tab_.phase(oi)) return false; + if (tab_.phase(i + nq) != other.tab_.phase(oi + nq)) return false; } return true; @@ -529,6 +581,14 @@ void UnitaryRevTableau::apply_S_at_end(const Qubit& qb) { tab_.apply_pauli_at_front(SpPauliStabiliser(qb, Pauli::Z), 3); } +void UnitaryRevTableau::apply_Z_at_front(const Qubit& qb) { + tab_.apply_Z_at_end(qb); +} + +void UnitaryRevTableau::apply_Z_at_end(const Qubit& qb) { + tab_.apply_Z_at_front(qb); +} + void UnitaryRevTableau::apply_V_at_front(const Qubit& qb) { tab_.apply_pauli_at_end(SpPauliStabiliser(qb, Pauli::X), 3); } @@ -537,6 +597,22 @@ void UnitaryRevTableau::apply_V_at_end(const Qubit& qb) { tab_.apply_pauli_at_front(SpPauliStabiliser(qb, Pauli::X), 3); } +void UnitaryRevTableau::apply_X_at_front(const Qubit& qb) { + tab_.apply_X_at_end(qb); +} + +void UnitaryRevTableau::apply_X_at_end(const Qubit& qb) { + tab_.apply_X_at_front(qb); +} + +void UnitaryRevTableau::apply_H_at_front(const Qubit& qb) { + tab_.apply_H_at_end(qb); +} + +void UnitaryRevTableau::apply_H_at_end(const Qubit& qb) { + tab_.apply_H_at_front(qb); +} + void UnitaryRevTableau::apply_CX_at_front( const Qubit& control, const Qubit& target) { tab_.apply_CX_at_end(control, target); @@ -549,14 +625,54 @@ void UnitaryRevTableau::apply_CX_at_end( void UnitaryRevTableau::apply_gate_at_front( OpType type, const qubit_vector_t& qbs) { - if (type != OpType::Phase) - tab_.apply_gate_at_end(get_op_ptr(type)->dagger()->get_type(), qbs); + // Handle types whose dagger is not an optype + switch (type) { + case OpType::ZZMax: { + tab_.apply_gate_at_end(OpType::ZZMax, qbs); + tab_.apply_gate_at_end(OpType::Z, {qbs.at(0)}); + tab_.apply_gate_at_end(OpType::Z, {qbs.at(1)}); + break; + } + case OpType::ISWAPMax: { + tab_.apply_gate_at_end(OpType::ISWAPMax, qbs); + tab_.apply_gate_at_end(OpType::Z, {qbs.at(0)}); + tab_.apply_gate_at_end(OpType::Z, {qbs.at(1)}); + break; + } + case OpType::Phase: { + break; + } + default: { + tab_.apply_gate_at_end(get_op_ptr(type)->dagger()->get_type(), qbs); + break; + } + } } void UnitaryRevTableau::apply_gate_at_end( OpType type, const qubit_vector_t& qbs) { - if (type != OpType::Phase) - tab_.apply_gate_at_front(get_op_ptr(type)->dagger()->get_type(), qbs); + // Handle types whose dagger is not an optype + switch (type) { + case OpType::ZZMax: { + tab_.apply_gate_at_front(OpType::ZZMax, qbs); + tab_.apply_gate_at_front(OpType::Z, {qbs.at(0)}); + tab_.apply_gate_at_front(OpType::Z, {qbs.at(1)}); + break; + } + case OpType::ISWAPMax: { + tab_.apply_gate_at_front(OpType::ISWAPMax, qbs); + tab_.apply_gate_at_front(OpType::Z, {qbs.at(0)}); + tab_.apply_gate_at_front(OpType::Z, {qbs.at(1)}); + break; + } + case OpType::Phase: { + break; + } + default: { + tab_.apply_gate_at_front(get_op_ptr(type)->dagger()->get_type(), qbs); + break; + } + } } void UnitaryRevTableau::apply_pauli_at_front( @@ -598,16 +714,16 @@ std::ostream& operator<<(std::ostream& os, const UnitaryRevTableau& tab) { unsigned nqs = tab.tab_.qubits_.size(); for (unsigned i = 0; i < nqs; ++i) { Qubit qi = tab.tab_.qubits_.right.at(i); - os << tab.tab_.tab_.xmat_.row(i) << " " << tab.tab_.tab_.zmat_.row(i) - << " " << tab.tab_.tab_.phase_(i) << "\t->\t" << "X@" << qi.repr() + os << tab.tab_.tab_.xmat.row(i) << " " << tab.tab_.tab_.zmat.row(i) + << " " << tab.tab_.tab_.phase(i) << "\t->\t" << "X@" << qi.repr() << std::endl; } os << "--" << std::endl; for (unsigned i = 0; i < nqs; ++i) { Qubit qi = tab.tab_.qubits_.right.at(i); - os << tab.tab_.tab_.xmat_.row(i + nqs) << " " - << tab.tab_.tab_.zmat_.row(i + nqs) << " " - << tab.tab_.tab_.phase_(i + nqs) << "\t->\t" << "Z@" << qi.repr() + os << tab.tab_.tab_.xmat.row(i + nqs) << " " + << tab.tab_.tab_.zmat.row(i + nqs) << " " + << tab.tab_.tab_.phase(i + nqs) << "\t->\t" << "Z@" << qi.repr() << std::endl; } return os; diff --git a/tket/src/Converters/ChoiMixTableauConverters.cpp b/tket/src/Converters/ChoiMixTableauConverters.cpp index 14bba5fcc2..624b0dc657 100644 --- a/tket/src/Converters/ChoiMixTableauConverters.cpp +++ b/tket/src/Converters/ChoiMixTableauConverters.cpp @@ -30,637 +30,732 @@ ChoiMixTableau circuit_to_cm_tableau(const Circuit& circ) { qubit_vector_t qbs = {args.begin(), args.end()}; tab.apply_gate(com.get_op_ptr()->get_type(), qbs); } + tab.rename_qubits( + circ.implicit_qubit_permutation(), + ChoiMixTableau::TableauSegment::Output); for (const Qubit& q : circ.discarded_qubits()) { tab.discard_qubit(q); } + tab.canonical_column_order(); + tab.gaussian_form(); return tab; } -std::pair cm_tableau_to_circuit(const ChoiMixTableau& t) { +struct ChoiMixBuilder { /** - * THE PLAN: - * We first identify and solve all post-selections by Gaussian elimination and - * diagonalisation of the input-only rows. This provides us with better - * guarantees about the form of the stabilizers generated by the remaining - * tableau - specifically that every possible stabilizer will involve at least - * some output portion, so out_qb will always be set in later sections of the - * synthesis. Diagonalisation of the input-only rows reduces them to just Z - * strings but not necessarily over a minimal number of qubits. This can be - * achieved by taking the Z matrix of these rows and performing row-wise - * Gaussian elimination. The leading 1s indicate which qubits we are isolating - * them onto and the rest gives the CXs required to reduce them to just the - * leading qubits. After all post-selections have been identified, we enter - * the main loop of handling all other inputs and reducing them one at a time - * to either an identity wire or Z-decoherence (OpType::Collapse) connected to - * an output, or to a discard. Let in_qb be this qubit to be solved. Pick a - * row containing X_in_qb (if one exists). Since it must contain some - * non-identity component in the output segment, we can pick one of these to - * be out_qb and apply unitary gates at the input and output to reduce the row - * to X_in_qb X_out_qb. We do the same with Z (if a row exists) to necessarily - * give Z_in_qb Z_out_qb by commutativity properties, or we pick out_qb here - * if no X row was found. If we had both an X row and a Z row, we have reduced - * it to an identity wire just fine. If we have only one of them, e.g. X_in_qb - * X_out_qb, we wish to make it the only row with X on either in_qb or out_qb - * to leave a decoherence channel. We note that any other row containing - * X_out_qb must have some other P_out2, since their combination cannot leave - * an empty output segment. So we can apply a unitary gate to eliminate the - * X_out_qb on the other row. By doing output-first gaussian elimination on - * the output sub-tableau ignoring the target row and X_out_qb column, for - * each such row with X_out_qb there is some other output column P_out2 for - * which it is the unique entry, so we can be sure that applying the unitary - * gate does not add X_out_qb onto other rows. Having this strategy available - * to make it the unique row with X_out_qb still allows us to make it the - * unique row with X_in_qb by row combinations. Once all rows with inputs have - * been eliminated, any unsolved inputs must have been discarded and the - * remaining tableau is an inverse diagonalisation tableau (only outputs). + * We will consider applying gates to either side of the tableau to reduce + * qubits down to one of a few simple states (identity, collapse, zero + * initialise, mixed initialised, post-selected, or discarded), allowing us to + * remove that qubit and continue until the tableau contains no qubits left. + * This gradually builds up a set of operations both before and after the + * working tableau. We should ask for the following combination of temporary + * states to compose to form a channel equivalent to that described by the + * input tableau: + * - in_circ: a unitary circuit (without implicit permutations) + * - post_selected: a set of post-selection actions to apply to the outputs of + * in_circ + * - discarded: a set of discard actions to apply to the outputs of in_circ + * - collapsed: a set of collapse actions to apply to the outputs of in_circ + * (i.e. decoherence in the Z basis) + * - tab: the remaining tableau still to be solved. Acting as an identity on + * any qubits not contained within tab + * - in_out_permutation: a permutation of the qubits, read as a map from the + * input qubit name to the output qubit it is sent to + * - zero_initialised: a set of initialisations of fresh output qubits + * (in_out_permutation may join these qubits onto input qubits that have been + * post-selected or discarded to reuse qubits) + * - mix_initialised: a set of initialisations of output qubits into + * maximally-mixed states (in_out_permutation may similarly join these onto + * input qubits no longer in use) + * - out_circ_tp: the transpose of a unitary circuit (without implicit + * permuations); we store this as the transpose so we can build it up in + * reverse */ + Circuit in_circ; + std::set post_selected; + std::set discarded; + std::set collapsed; + ChoiMixTableau tab; + boost::bimap in_out_permutation; + std::set zero_initialised; + std::set mix_initialised; + Circuit out_circ_tp; + + // The CXConfigType preferred when invoking diagonalisation techniques + CXConfigType cx_config; + + // Additional qubit names (distinct from qubits already on the respective + // segment of the tableau) than can be used for zero initialisations and + // post-selections when synthesising a unitary extension + qubit_vector_t unitary_init_names; + qubit_vector_t unitary_post_names; + + // Initialises the builder with a tableau + explicit ChoiMixBuilder(const ChoiMixTableau& tab, CXConfigType cx_config); + // For synthesis of a unitary extension, initialises the builder with a + // tableau and some additional qubit names which the resulting circuit may + // optionally use, representing zero-initialised or post-selected qubits. + // These are the qubits on which we can freely add Zs to the rows of the given + // tableau to guarantee a unitary extension; none of these names should appear + // in tab + explicit ChoiMixBuilder( + const ChoiMixTableau& tab, CXConfigType cx_config, + const qubit_vector_t& init_names, const qubit_vector_t& post_names); + + // Debug method: applies all staged operations back onto tab to provide the + // tableau that the synthesis result is currently aiming towards. In exact + // synthesis, this should remain invariant during synthesis. For synthesis of + // a unitary extension, this should at least span the rows of the original + // tableau up to additional Zs on spare input and output qubits. + ChoiMixTableau realised_tableau() const; - // Operate on a copy of the tableau to track the remainder as we extract gates - ChoiMixTableau tab(t); + /** + * STAGES OF SYNTHESIS + */ - // Canonicalise tableau and perform output-first gaussian elimination to - // isolate post-selected subspace in lower rows - tab.canonical_column_order(ChoiMixTableau::TableauSegment::Output); - tab.gaussian_form(); + // Match up pairs of generators that anti-commute in the input segment but + // commute with all others; such pairs of rows reduce to an identity wire + // between a pair of qubits; we solve this by pairwise Pauli reduction methods + void solve_id_subspace(); + // After removing the identity subspace, all remaining rows mutually commute + // within each tableau segment; diagonalise each segment individually + void diagonalise_segments(); + // Solve the post-selected subspace which has already been diagonalised + void solve_postselected_subspace(); + // Solve the zero-initialised subspace which has already been diagonalised + void solve_initialised_subspace(); + // All remaining rows are in the collapsed subspace (each row is the unique + // stabilizer passing through some Collapse gate); solve it + void solve_collapsed_subspace(); + + // Simplifies the tableau by removing qubits on which all rows have I; such + // qubits are either discarded inputs or mixed-initialised outputs + void remove_unused_qubits(); + // For synthesis of a unitary extension, match up qubits from + // post-selected/zero-initialised with unitary_post_names/unitary_init_names + // and add them to in_out_permutation + void assign_init_post_names(); + // Fill out in_out_permutation to map all qubits; this typically takes the + // form of a standard qubit reuse pattern (e.g. discard and reinitialise) + void assign_remaining_names(); + // Once tab has been completely reduced to no rows and no qubits, compose the + // staged operations to build the output circuit and return the renaming map + // from output names of original tableau to the qubits of the returned circuit + // they are mapped to + std::pair output_circuit(); + std::pair unitary_output_circuit(); +}; + +std::pair cm_tableau_to_exact_circuit( + const ChoiMixTableau& tab, CXConfigType cx_config) { + ChoiMixBuilder builder(tab, cx_config); + builder.remove_unused_qubits(); + builder.solve_id_subspace(); + builder.diagonalise_segments(); + builder.solve_postselected_subspace(); + builder.solve_initialised_subspace(); + builder.solve_collapsed_subspace(); + builder.remove_unused_qubits(); + builder.assign_remaining_names(); + return builder.output_circuit(); +} - // Set up circuits for extracting gates into (in_circ is set up after - // diagonalising the post-selected subspace) - qubit_vector_t input_qubits, output_qubits; +std::pair cm_tableau_to_unitary_extension_circuit( + const ChoiMixTableau& tab, const std::vector& init_names, + const std::vector& post_names, CXConfigType cx_config) { + ChoiMixBuilder builder(tab, cx_config, init_names, post_names); + builder.remove_unused_qubits(); + builder.solve_id_subspace(); + builder.diagonalise_segments(); + builder.solve_postselected_subspace(); + builder.solve_initialised_subspace(); + builder.solve_collapsed_subspace(); + builder.remove_unused_qubits(); + builder.assign_init_post_names(); + builder.assign_remaining_names(); + return builder.unitary_output_circuit(); +} + +ChoiMixBuilder::ChoiMixBuilder(const ChoiMixTableau& t, CXConfigType cx) + : ChoiMixBuilder(t, cx, {}, {}) {} + +ChoiMixBuilder::ChoiMixBuilder( + const ChoiMixTableau& t, CXConfigType cx, const qubit_vector_t& inits, + const qubit_vector_t& posts) + : in_circ(), + post_selected(), + discarded(), + collapsed(), + tab(t), + in_out_permutation(), + zero_initialised(), + mix_initialised(), + out_circ_tp(), + cx_config(cx), + unitary_init_names(inits), + unitary_post_names(posts) { for (unsigned i = 0; i < tab.get_n_boundaries(); ++i) { ChoiMixTableau::col_key_t key = tab.col_index_.right.at(i); if (key.second == ChoiMixTableau::TableauSegment::Input) - input_qubits.push_back(key.first); + in_circ.add_qubit(key.first); else - output_qubits.push_back(key.first); + out_circ_tp.add_qubit(key.first); } - Circuit out_circ_tp(output_qubits, {}); - unit_map_t join_permutation; - std::set post_selected; - std::map> in_x_row; - std::map> in_z_row; - boost::bimap matched_qubits; - - // Call diagonalisation methods to diagonalise post-selected subspace - std::list to_diag; - for (unsigned r = tab.get_n_rows(); r > 0;) { - --r; - ChoiMixTableau::row_tensor_t rten = tab.get_row(r); - if (rten.second.size() != 0) { - // Reached the rows with non-empty output segment - break; - } - // Else, we add the row to the subspace - to_diag.push_back(rten.first); + for (const Qubit& init_q : unitary_init_names) { + if (tab.col_index_.left.find(ChoiMixTableau::col_key_t{ + init_q, ChoiMixTableau::TableauSegment::Input}) != + tab.col_index_.left.end()) + throw std::logic_error( + "Free qubit name for initialisation conflicts with existing live " + "input of ChoiMixTableau"); } - unsigned post_selected_size = to_diag.size(); - std::set diag_ins{input_qubits.begin(), input_qubits.end()}; - Circuit in_circ = mutual_diagonalise(to_diag, diag_ins, CXConfigType::Tree); - // Extract the dagger of each gate in order from tab - for (const Command& com : in_circ) { - auto args = com.get_args(); - qubit_vector_t qbs = {args.begin(), args.end()}; - tab.apply_gate( - com.get_op_ptr()->dagger()->get_type(), qbs, - ChoiMixTableau::TableauSegment::Input); + for (const Qubit& post_q : unitary_post_names) { + if (tab.col_index_.left.find(ChoiMixTableau::col_key_t{ + post_q, ChoiMixTableau::TableauSegment::Output}) != + tab.col_index_.left.end()) + throw std::logic_error( + "Free qubit name for post-selection conflicts with existing live " + "output of ChoiMixTableau"); } +} - // Diagonalised rows should still be at the bottom of the tableau - reduce - // them to a minimal set of qubits for post-selection by first reducing to - // upper echelon form - std::vector> row_ops = - gaussian_elimination_row_ops( - tab.tab_.zmat_.bottomRows(post_selected_size)); - for (const std::pair& op : row_ops) { - tab.tab_.row_mult(op.first, op.second); - } - // Obtain CX instructions as column operations - std::vector> col_ops = - gaussian_elimination_col_ops( - tab.tab_.zmat_.bottomRows(post_selected_size)); - // These gates will also swap qubits to isolate the post-selections on the - // first few qubits - this is fine as can be cleaned up later with peephole - // optimisations - for (const std::pair& op : col_ops) { - tab.tab_.apply_CX(op.first, op.second); - ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(op.first); - ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(op.second); - in_circ.add_op(OpType::CX, {ctrl.first, trgt.first}); +// DEBUG METHOD: Ignore this for coverage checks +// GCOVR_EXCL_START +ChoiMixTableau ChoiMixBuilder::realised_tableau() const { + ChoiMixTableau in_tab = circuit_to_cm_tableau(in_circ); + for (const Qubit& q : post_selected) + in_tab.post_select(q, ChoiMixTableau::TableauSegment::Output); + for (const Qubit& q : discarded) + in_tab.discard_qubit(q, ChoiMixTableau::TableauSegment::Output); + for (const Qubit& q : collapsed) + in_tab.collapse_qubit(q, ChoiMixTableau::TableauSegment::Output); + ChoiMixTableau out_tab = circuit_to_cm_tableau(out_circ_tp.transpose()); + for (const Qubit& q : zero_initialised) + out_tab.post_select(q, ChoiMixTableau::TableauSegment::Input); + for (const Qubit& q : mix_initialised) + out_tab.discard_qubit(q, ChoiMixTableau::TableauSegment::Input); + qubit_map_t out_in_permutation{}; + using perm_entry = boost::bimap::left_const_reference; + BOOST_FOREACH (perm_entry entry, in_out_permutation.left) { + out_in_permutation.insert({entry.second, entry.first}); } + out_tab.rename_qubits( + out_in_permutation, ChoiMixTableau::TableauSegment::Output); + return ChoiMixTableau::compose(ChoiMixTableau::compose(in_tab, tab), out_tab); +} +// GCOVR_EXCL_STOP - // Post-select rows - // TODO Post-selection op is not yet available in tket - replace this once - // implemented; an implementation hint is provided below - if (post_selected_size != 0) - throw std::logic_error( - "Not yet implemented: post-selection required during ChoiMixTableau " - "synthesis"); - // for (unsigned r = 0; r < post_selected_size; ++r) { - // ChoiMixTableau::row_tensor_t row = tab.get_row(tab.get_n_rows() - 1); - // if (row.second.string.map.size() != 0 || row.first.string.map.size() != 1 - // || row.first.string.map.begin()->second != Pauli::Z) throw - // std::logic_error("Unexpected error during post-selection identification - // in ChoiMixTableau synthesis"); Qubit post_selected_qb = - // row.first.string.map.begin()->first; if (row.second.coeff == -1.) - // in_circ.add_op(OpType::X, {post_selected_qb}); - // tab.remove_row(tab.get_n_rows() - 1); - // post_selected.insert(post_selected_qb); - // throw std::logic_error("Not yet implemented: post-selection required - // during ChoiMixTableau synthesis"); - // } - +void ChoiMixBuilder::solve_id_subspace() { // Input-first gaussian elimination to solve input-sides of remaining rows tab.canonical_column_order(ChoiMixTableau::TableauSegment::Input); tab.gaussian_form(); - // Iterate through remaining inputs and reduce output portion to a single - // qubit - for (const Qubit& in_qb : input_qubits) { - // Skip post-selected qubits - if (post_selected.find(in_qb) != post_selected.end()) continue; - - unsigned col = tab.col_index_.left.at(ChoiMixTableau::col_key_t{ - in_qb, ChoiMixTableau::TableauSegment::Input}); - std::optional out_qb = std::nullopt; + std::set solved_rows; + std::set solved_ins, solved_outs; + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + if (solved_rows.find(r) != solved_rows.end()) continue; - // Find the row with X_in_qb (if one exists) - std::optional x_row = std::nullopt; - for (unsigned r = 0; r < tab.get_n_rows(); ++r) { - if (tab.tab_.xmat_(r, col)) { - x_row = r; - break; - } + // Look for a row which anticommutes with row r over the inputs + std::list xcols, zcols; + for (unsigned c = 0; c < tab.get_n_inputs(); ++c) { + if (tab.tab_.xmat(r, c)) xcols.push_back(c); + if (tab.tab_.zmat(r, c)) zcols.push_back(c); } - - if (x_row) { - // A possible optimisation could involve row multiplications to reduce the - // Hamming weight of the row before applying gates, but minimising this - // would solve minimum weight/distance of a binary linear code whose - // decision problem is NP-complete (shown by Vardy, "The intractability of - // computing the minimum distance of a code", 1997). Just settle on using - // the first row for now, reducing the input and output to a single qubit - ChoiMixTableau::row_tensor_t row_paulis = tab.get_row(*x_row); - for (const std::pair& p : row_paulis.second.string) { - if (matched_qubits.right.find(p.first) == matched_qubits.right.end()) { - out_qb = p.first; - matched_qubits.insert({in_qb, *out_qb}); - break; - } + for (unsigned r2 = r + 1; r2 < tab.get_n_rows(); ++r2) { + if (solved_rows.find(r2) != solved_rows.end()) continue; + bool anti = false; + for (const unsigned c : xcols) anti ^= tab.tab_.zmat(r2, c); + for (const unsigned c : zcols) anti ^= tab.tab_.xmat(r2, c); + if (!anti) continue; + + // Found a candidate pair of rows. Because of the Gaussian elimination, it + // is more likely that the first mismatching qubit is X for r and Z for + // r2, so favour reducing r2 to Z and r to X + ChoiMixTableau::row_tensor_t row_r = tab.get_row(r); + ChoiMixTableau::row_tensor_t row_r2 = tab.get_row(r2); + std::pair in_diag_circ = + reduce_anticommuting_paulis_to_z_x( + row_r2.first, row_r.first, cx_config); + in_circ.append(in_diag_circ.first); + for (const Command& com : in_diag_circ.first) { + auto args = com.get_args(); + qubit_vector_t qbs = {args.begin(), args.end()}; + tab.apply_gate( + com.get_op_ptr()->dagger()->get_type(), qbs, + ChoiMixTableau::TableauSegment::Input); } - // Reduce input string to just X_in_qb - if (row_paulis.first.get(in_qb) == Pauli::Y) { - // If it is a Y, extract an Sdg gate so the Pauli is exactly X - in_circ.add_op(OpType::Sdg, {in_qb}); - tab.apply_S(in_qb, ChoiMixTableau::TableauSegment::Input); - } - for (const std::pair& qbp : row_paulis.first.string) { - if (qbp.first == in_qb) continue; - // Extract an entangling gate to eliminate the qubit - switch (qbp.second) { - case Pauli::X: { - in_circ.add_op(OpType::CX, {in_qb, qbp.first}); - tab.apply_CX( - in_qb, qbp.first, ChoiMixTableau::TableauSegment::Input); - break; - } - case Pauli::Y: { - in_circ.add_op(OpType::CY, {in_qb, qbp.first}); - tab.apply_gate( - OpType::CY, {in_qb, qbp.first}, - ChoiMixTableau::TableauSegment::Input); - break; - } - case Pauli::Z: { - in_circ.add_op(OpType::CZ, {in_qb, qbp.first}); - tab.apply_gate( - OpType::CZ, {in_qb, qbp.first}, - ChoiMixTableau::TableauSegment::Input); - break; - } - default: { - break; - } - } + // Since the full rows must commute but they anticommute over the inputs, + // they must also anticommute over the outputs; we similarly reduce these + // down to Z and X + std::pair out_diag_circ_dag = + reduce_anticommuting_paulis_to_z_x( + row_r2.second, row_r.second, cx_config); + out_circ_tp.append(out_diag_circ_dag.first.dagger().transpose()); + for (const Command& com : out_diag_circ_dag.first) { + auto args = com.get_args(); + qubit_vector_t qbs = {args.begin(), args.end()}; + tab.apply_gate( + com.get_op_ptr()->get_type(), qbs, + ChoiMixTableau::TableauSegment::Output); } - // And then the same for X_out_qb - if (row_paulis.second.get(*out_qb) == Pauli::Y) { - // If it is a Y, extract an Sdg gate so the Pauli is exactly X - out_circ_tp.add_op(OpType::Sdg, {*out_qb}); - tab.apply_S(*out_qb, ChoiMixTableau::TableauSegment::Output); - } else if (row_paulis.second.get(*out_qb) == Pauli::Z) { - // If it is a Z, extract an Vdg and Sdg gate so the Pauli is exactly X - out_circ_tp.add_op(OpType::Vdg, {*out_qb}); - out_circ_tp.add_op(OpType::Sdg, {*out_qb}); - tab.apply_V(*out_qb, ChoiMixTableau::TableauSegment::Output); - tab.apply_S(*out_qb, ChoiMixTableau::TableauSegment::Output); + // Check that rows have been successfully reduced + row_r = tab.get_row(r); + row_r2 = tab.get_row(r2); + if (row_r.first.size() != 1 || + row_r.first.string.begin()->second != Pauli::X || + row_r.second.size() != 1 || + row_r.second.string.begin()->second != Pauli::X || + row_r2.first.size() != 1 || + row_r2.first.string.begin()->second != Pauli::Z || + row_r2.second.size() != 1 || + row_r2.second.string.begin()->second != Pauli::Z) + throw std::logic_error( + "Unexpected error during identity reduction in ChoiMixTableau " + "synthesis"); + // Solve phases + if (row_r.second.is_real_negative()) { + in_circ.add_op(OpType::Z, {in_diag_circ.second}); + tab.apply_gate( + OpType::Z, {in_diag_circ.second}, + ChoiMixTableau::TableauSegment::Input); } - for (const std::pair& qbp : - row_paulis.second.string) { - if (qbp.first == *out_qb) continue; - // Extract an entangling gate to eliminate the qubit - switch (qbp.second) { - case Pauli::X: { - out_circ_tp.add_op(OpType::CX, {*out_qb, qbp.first}); - tab.apply_CX( - *out_qb, qbp.first, ChoiMixTableau::TableauSegment::Output); - break; - } - case Pauli::Y: { - // CY does not have a transpose OpType defined so decompose - out_circ_tp.add_op(OpType::S, {qbp.first}); - out_circ_tp.add_op(OpType::CX, {*out_qb, qbp.first}); - out_circ_tp.add_op(OpType::Sdg, {qbp.first}); - tab.apply_gate( - OpType::CY, {*out_qb, qbp.first}, - ChoiMixTableau::TableauSegment::Output); - break; - } - case Pauli::Z: { - out_circ_tp.add_op(OpType::CZ, {*out_qb, qbp.first}); - tab.apply_gate( - OpType::CZ, {*out_qb, qbp.first}, - ChoiMixTableau::TableauSegment::Output); - break; - } - default: { - break; - } - } + if (row_r2.second.is_real_negative()) { + in_circ.add_op(OpType::X, {in_diag_circ.second}); + tab.apply_gate( + OpType::X, {in_diag_circ.second}, + ChoiMixTableau::TableauSegment::Input); } - } - - // Find the row with Z_in_qb (if one exists) - std::optional z_row = std::nullopt; - for (unsigned r = 0; r < tab.get_n_rows(); ++r) { - if (tab.tab_.zmat_(r, col)) { - z_row = r; - break; + // Connect in permutation + in_out_permutation.insert( + {in_diag_circ.second, out_diag_circ_dag.second}); + solved_rows.insert(r); + solved_rows.insert(r2); + + // Remove these solved qubits from other rows; by commutation of rows, a + // row contains Z@in_diag_circ.second iff it contains + // Z@out_diag_circ_dag.second and similarly for X + unsigned in_c = tab.col_index_.left.at(ChoiMixTableau::col_key_t{ + in_diag_circ.second, ChoiMixTableau::TableauSegment::Input}); + for (unsigned r3 = 0; r3 < tab.get_n_rows(); ++r3) { + if (r3 != r && tab.tab_.xmat(r3, in_c)) tab.tab_.row_mult(r, r3); + if (r3 != r2 && tab.tab_.zmat(r3, in_c)) tab.tab_.row_mult(r2, r3); } - } - if (z_row) { - // If both an X and Z row exist, then out_qb should have a - // value and the rows should have anticommuting Paulis on out_qb to - // preserve commutativity of rows - ChoiMixTableau::row_tensor_t row_paulis = tab.get_row(*z_row); - if (!x_row) { - for (const std::pair& p : - row_paulis.second.string) { - if (matched_qubits.right.find(p.first) == - matched_qubits.right.end()) { - out_qb = p.first; - matched_qubits.insert({in_qb, *out_qb}); - break; - } - } - } - - // Reduce input string to just Z_in_qb. - // No need to consider different paulis on in_qb: if we had X or Y, this - // row would have been identified as x_row instead of Z row, and if - // another row was already chosen as x_row then canonical gaussian form - // would imply all other rows do not contain X_in_qb or Y_in_qb - for (const std::pair& qbp : row_paulis.first.string) { - if (qbp.first == in_qb) continue; - // Extract an entangling gate to eliminate the qubit - switch (qbp.second) { - case Pauli::X: { - in_circ.add_op(OpType::H, {qbp.first}); - in_circ.add_op(OpType::CX, {qbp.first, in_qb}); - tab.apply_gate( - OpType::H, {qbp.first}, ChoiMixTableau::TableauSegment::Input); - tab.apply_CX( - qbp.first, in_qb, ChoiMixTableau::TableauSegment::Input); - break; - } - case Pauli::Y: { - in_circ.add_op(OpType::Vdg, {qbp.first}); - in_circ.add_op(OpType::CX, {qbp.first, in_qb}); - tab.apply_V(qbp.first, ChoiMixTableau::TableauSegment::Input); - tab.apply_CX( - qbp.first, in_qb, ChoiMixTableau::TableauSegment::Input); - break; - } - case Pauli::Z: { - in_circ.add_op(OpType::CX, {qbp.first, in_qb}); - tab.apply_CX( - qbp.first, in_qb, ChoiMixTableau::TableauSegment::Input); - break; - } - default: { - break; - } - } - } - - // And then reduce output string to just Z_out_qb - if (row_paulis.second.get(*out_qb) == Pauli::Y) { - // If it is a Y, extract a Vdg gate so the Pauli is exactly Z - out_circ_tp.add_op(OpType::Vdg, {*out_qb}); - tab.apply_V(*out_qb, ChoiMixTableau::TableauSegment::Output); - } else if (row_paulis.second.get(*out_qb) == Pauli::X) { - // If it is an X, extract an Sdg and Vdg gate so the Pauli is exactly Z - // We do not need to care about messing up the X row here since if we - // solved an X row then this row can't also have X by commutativity - out_circ_tp.add_op(OpType::Sdg, {*out_qb}); - out_circ_tp.add_op(OpType::Vdg, {*out_qb}); - tab.apply_S(*out_qb, ChoiMixTableau::TableauSegment::Output); - tab.apply_V(*out_qb, ChoiMixTableau::TableauSegment::Output); - } - for (const std::pair& qbp : - row_paulis.second.string) { - if (qbp.first == *out_qb) continue; - // Extract an entangling gate to eliminate the qubit - switch (qbp.second) { - case Pauli::X: { - out_circ_tp.add_op(OpType::H, {qbp.first}); - out_circ_tp.add_op(OpType::CX, {qbp.first, *out_qb}); - tab.apply_gate( - OpType::H, {qbp.first}, ChoiMixTableau::TableauSegment::Output); - tab.apply_CX( - qbp.first, *out_qb, ChoiMixTableau::TableauSegment::Output); - break; - } - case Pauli::Y: { - out_circ_tp.add_op(OpType::Vdg, {qbp.first}); - out_circ_tp.add_op(OpType::CX, {qbp.first, *out_qb}); - tab.apply_V(qbp.first, ChoiMixTableau::TableauSegment::Output); - tab.apply_CX( - qbp.first, *out_qb, ChoiMixTableau::TableauSegment::Output); - break; - } - case Pauli::Z: { - out_circ_tp.add_op(OpType::CX, {qbp.first, *out_qb}); - tab.apply_CX( - qbp.first, *out_qb, ChoiMixTableau::TableauSegment::Output); - break; - } - default: { - break; - } - } - } - } - - in_x_row.insert({in_qb, x_row}); - in_z_row.insert({in_qb, z_row}); - } - - // X and Z row are guaranteed to be the unique rows with X_in_qb and Z_in_qb. - // If the Z row exists, then all rows must commute with Z_in_qb Z_out_qb, so - // if any row has X_out_qb it must also have X_in_qb. Hence if both exist then - // the X row is also the unique row with X_out_qb, and by similar argument the - // Z row is the unique row with Z_out_qb. However, if only one row exists, - // this uniqueness may not hold but can be forced by applying some unitary - // gates to eliminate e.g. Z_out_qb from all other rows. We use a gaussian - // elimination subroutine to identify a combination of gates that won't add - // Z_out_qb onto other qubits. Since we have already removed them from all - // other rows containing input components, we only need to consider the - // output-only rows, which all exist at the bottom of the tableau. - unsigned out_stabs = 0; - while (out_stabs < tab.get_n_rows()) { - ChoiMixTableau::row_tensor_t rten = - tab.get_row(tab.get_n_rows() - 1 - out_stabs); - if (rten.first.size() != 0) { - // Reached the rows with non-empty input segment + solved_ins.insert({in_diag_circ.second}); + solved_outs.insert({out_diag_circ_dag.second}); break; } - ++out_stabs; - } - MatrixXb out_rows = MatrixXb::Zero(out_stabs, 2 * output_qubits.size()); - out_rows(Eigen::all, Eigen::seq(0, Eigen::last, 2)) = tab.tab_.xmat_.block( - tab.get_n_rows() - out_stabs, input_qubits.size(), out_stabs, - output_qubits.size()); - out_rows(Eigen::all, Eigen::seq(1, Eigen::last, 2)) = tab.tab_.zmat_.block( - tab.get_n_rows() - out_stabs, input_qubits.size(), out_stabs, - output_qubits.size()); - // Identify other qubits to apply gates to for removing the extra matched - // output terms - MatrixXb unmatched_outs = MatrixXb::Zero( - out_stabs, 2 * (output_qubits.size() - matched_qubits.size())); - std::vector> col_lookup; - for (unsigned out_col = 0; out_col < output_qubits.size(); ++out_col) { - unsigned tab_col = input_qubits.size() + out_col; - Qubit out_qb = tab.col_index_.right.at(tab_col).first; - if (matched_qubits.right.find(out_qb) != matched_qubits.right.end()) - continue; - unmatched_outs.col(col_lookup.size()) = out_rows.col(2 * out_col); - col_lookup.push_back({out_qb, Pauli::X}); - unmatched_outs.col(col_lookup.size()) = out_rows.col(2 * out_col + 1); - col_lookup.push_back({out_qb, Pauli::Z}); - } - row_ops = gaussian_elimination_row_ops(unmatched_outs); - for (const std::pair& op : row_ops) { - for (unsigned c = 0; c < unmatched_outs.cols(); ++c) { - unmatched_outs(op.second, c) = - unmatched_outs(op.second, c) ^ unmatched_outs(op.first, c); - } - tab.tab_.row_mult( - tab.get_n_rows() - out_stabs + op.first, - tab.get_n_rows() - out_stabs + op.second); - } - // Go through the output-only rows and remove terms from matched qubits - for (unsigned r = 0; r < out_stabs; ++r) { - unsigned leading_col = 0; - for (unsigned c = 0; c < unmatched_outs.cols(); ++c) { - if (unmatched_outs(r, c)) { - leading_col = c; - break; - } - } - std::pair alternate_qb = col_lookup.at(leading_col); - - if (alternate_qb.second != Pauli::X) { - // Make the alternate point of contact X so we only need one set of rule - // for eliminating Paulis - tab.apply_gate( - OpType::H, {alternate_qb.first}, - ChoiMixTableau::TableauSegment::Output); - out_circ_tp.add_op(OpType::H, {alternate_qb.first}); - } - - ChoiMixTableau::row_tensor_t row_paulis = - tab.get_row(tab.get_n_rows() - out_stabs + r); - - for (const std::pair& qbp : row_paulis.second.string) { - if (matched_qubits.right.find(qbp.first) == matched_qubits.right.end()) - continue; - // Alternate point is guaranteed to be unmatched, so always needs an - // entangling gate - switch (qbp.second) { - case Pauli::X: { - out_circ_tp.add_op( - OpType::CX, {alternate_qb.first, qbp.first}); - tab.apply_CX( - alternate_qb.first, qbp.first, - ChoiMixTableau::TableauSegment::Output); - break; - } - case Pauli::Z: { - out_circ_tp.add_op( - OpType::CZ, {alternate_qb.first, qbp.first}); - tab.apply_gate( - OpType::CZ, {alternate_qb.first, qbp.first}, - ChoiMixTableau::TableauSegment::Output); - break; - } - default: { - // Don't have to care about Y since any matched qubit has a row that - // is reduced to either X or Z and all other rows must commute with - // that - break; - } - } - } } - // Now that X_in_qb X_out_qb (or Zs) is the unique row for each of X_in_qb and - // X_out_qb, we can actually link up the qubit wires and remove the rows - for (const Qubit& in_qb : input_qubits) { - if (post_selected.find(in_qb) != post_selected.end()) continue; - - std::optional x_row = in_x_row.at(in_qb); - std::optional z_row = in_z_row.at(in_qb); - auto found = matched_qubits.left.find(in_qb); - std::optional out_qb = (found == matched_qubits.left.end()) - ? std::nullopt - : std::optional{found->second}; - // Handle phases and resolve qubit connections - if (x_row) { - if (z_row) { - // Hook up with an identity wire - if (tab.tab_.phase_(*z_row)) in_circ.add_op(OpType::X, {in_qb}); - if (tab.tab_.phase_(*x_row)) in_circ.add_op(OpType::Z, {in_qb}); - join_permutation.insert({*out_qb, in_qb}); - } else { - // Just an X row, so must be connected to out_qb via a decoherence - if (tab.tab_.phase_(*x_row)) in_circ.add_op(OpType::Z, {in_qb}); - in_circ.add_op(OpType::H, {in_qb}); - in_circ.add_op(OpType::Collapse, {in_qb}); - in_circ.add_op(OpType::H, {in_qb}); - join_permutation.insert({*out_qb, in_qb}); - } - } else { - if (z_row) { - // Just a Z row, so must be connected to out_qb via a decoherence - if (tab.tab_.phase_(*z_row)) in_circ.add_op(OpType::X, {in_qb}); - in_circ.add_op(OpType::Collapse, {in_qb}); - join_permutation.insert({*out_qb, in_qb}); - } else { - // No rows involving this input, so it is discarded - in_circ.qubit_discard(in_qb); + // Remove solved rows and qubits from tableau; since removing rows/columns + // replaces them with the row/column from the end, remove in reverse order + for (auto it = solved_rows.rbegin(); it != solved_rows.rend(); ++it) + tab.remove_row(*it); + for (auto it = solved_ins.rbegin(); it != solved_ins.rend(); ++it) + tab.discard_qubit(*it, ChoiMixTableau::TableauSegment::Input); + for (auto it = solved_outs.rbegin(); it != solved_outs.rend(); ++it) + tab.discard_qubit(*it, ChoiMixTableau::TableauSegment::Output); +} + +// Given a matrix that is already in upper echelon form, use the fact that the +// leading columns are already unique to give column operations that reduce it +// down to identity over the leading columns, eliminating extra swap gates to +// move to the first spaces +static std::vector> +leading_column_gaussian_col_ops(const MatrixXb& source) { + std::vector col_list; + std::set non_leads; + for (unsigned r = 0; r < source.rows(); ++r) { + bool leading_found = false; + for (unsigned c = 0; c < source.cols(); ++c) { + if (source(r, c)) { + if (leading_found) + non_leads.insert(c); + else { + leading_found = true; + col_list.push_back(c); + } } } } + for (const unsigned& c : non_leads) col_list.push_back(c); + MatrixXb reordered = MatrixXb::Zero(source.rows(), col_list.size()); + for (unsigned c = 0; c < col_list.size(); ++c) + reordered.col(c) = source.col(col_list.at(c)); + std::vector> reordered_ops = + gaussian_elimination_col_ops(reordered); + std::vector> res; + for (const std::pair& op : reordered_ops) + res.push_back({col_list.at(op.first), col_list.at(op.second)}); + return res; +} - // Remove rows with inputs from the tableau - tab.tab_ = SymplecticTableau( - tab.tab_.xmat_.bottomRows(out_stabs), - tab.tab_.zmat_.bottomRows(out_stabs), tab.tab_.phase_.tail(out_stabs)); - - // Can't use a template with multiple parameters within a macro since the - // comma will register as an argument delimiter for the macro - using match_entry = boost::bimap::left_const_reference; - BOOST_FOREACH (match_entry entry, matched_qubits.left) { - tab.discard_qubit(entry.first, ChoiMixTableau::TableauSegment::Input); - tab.discard_qubit(entry.second, ChoiMixTableau::TableauSegment::Output); - } +void ChoiMixBuilder::diagonalise_segments() { + // Canonicalise tableau tab.canonical_column_order(ChoiMixTableau::TableauSegment::Output); + tab.gaussian_form(); - // Only remaining rows must be completely over the outputs. Call - // diagonalisation methods to diagonalise coherent subspace - to_diag.clear(); + // Set up diagonalisation tasks + std::list to_diag_ins, to_diag_outs; for (unsigned r = 0; r < tab.get_n_rows(); ++r) { ChoiMixTableau::row_tensor_t rten = tab.get_row(r); - to_diag.push_back(rten.second); + if (!rten.first.string.empty()) to_diag_ins.push_back(rten.first); + if (!rten.second.string.empty()) to_diag_outs.push_back(rten.second); } - std::set diag_outs; - for (const Qubit& out : output_qubits) { - if (matched_qubits.right.find(out) == matched_qubits.right.end()) - diag_outs.insert(out); + qubit_vector_t input_qubits = tab.input_qubits(); + std::set diag_ins{input_qubits.begin(), input_qubits.end()}; + Circuit in_diag_circ = mutual_diagonalise(to_diag_ins, diag_ins, cx_config); + for (const Command& com : in_diag_circ) { + auto args = com.get_args(); + in_circ.add_op(com.get_op_ptr(), args); + qubit_vector_t qbs = {args.begin(), args.end()}; + tab.apply_gate( + com.get_op_ptr()->dagger()->get_type(), qbs, + ChoiMixTableau::TableauSegment::Input); } + qubit_vector_t output_qubits = tab.output_qubits(); + std::set diag_outs{output_qubits.begin(), output_qubits.end()}; Circuit out_diag_circ = - mutual_diagonalise(to_diag, diag_outs, CXConfigType::Tree); - // Extract the dagger of each gate in order from tab + mutual_diagonalise(to_diag_outs, diag_outs, cx_config); for (const Command& com : out_diag_circ) { auto args = com.get_args(); + out_circ_tp.add_op(com.get_op_ptr()->dagger()->transpose(), args); qubit_vector_t qbs = {args.begin(), args.end()}; - tab.apply_gate(com.get_op_ptr()->get_type(), qbs); - out_circ_tp.add_op(com.get_op_ptr()->transpose()->dagger(), qbs); + tab.apply_gate( + com.get_op_ptr()->get_type(), qbs, + ChoiMixTableau::TableauSegment::Output); } - // All rows are diagonalised so we can just focus on the Z matrix. Reduce them - // to a minimal set of qubits for initialisation by first reducing to upper - // echelon form - row_ops = gaussian_elimination_row_ops(tab.tab_.zmat_); - for (const std::pair& op : row_ops) { - tab.tab_.row_mult(op.first, op.second); + // All rows are diagonalised, so we can just focus on the Z matrix + if (tab.tab_.xmat != MatrixXb::Zero(tab.get_n_rows(), tab.get_n_boundaries())) + throw std::logic_error( + "Diagonalisation in ChoiMixTableau synthesis failed"); +} + +void ChoiMixBuilder::solve_postselected_subspace() { + // As column order is currently output first, gaussian form will reveal the + // post-selected space at the bottom of the tableau and the submatrix of those + // rows will already be in upper echelon form + tab.gaussian_form(); + // Reduce them to a minimal set of qubits using CX gates + unsigned n_postselected = 0; + for (; n_postselected < tab.get_n_rows(); ++n_postselected) { + if (!tab.get_row(tab.get_n_rows() - 1 - n_postselected) + .second.string.empty()) + break; } - // Obtain CX instructions as column operations - col_ops = gaussian_elimination_col_ops(tab.tab_.zmat_); + unsigned n_ins = tab.get_n_inputs(); + unsigned n_outs = tab.get_n_outputs(); + MatrixXb subtableau = tab.tab_.zmat.bottomRightCorner(n_postselected, n_ins); + std::vector> col_ops = + leading_column_gaussian_col_ops(subtableau); for (const std::pair& op : col_ops) { - tab.tab_.apply_CX(op.second, op.first); - ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(op.second); - ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(op.first); - out_circ_tp.add_op(OpType::CX, {ctrl.first, trgt.first}); + unsigned tab_ctrl_col = n_outs + op.second; + unsigned tab_trgt_col = n_outs + op.first; + tab.tab_.apply_CX(tab_ctrl_col, tab_trgt_col); + ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(tab_ctrl_col); + ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(tab_trgt_col); + in_circ.add_op(OpType::CX, {ctrl.first, trgt.first}); } + // Postselect rows + for (unsigned r = 0; r < n_postselected; ++r) { + unsigned final_row = tab.get_n_rows() - 1; + ChoiMixTableau::row_tensor_t row = tab.get_row(final_row); + if (row.second.size() != 0 || row.first.size() != 1 || + row.first.string.begin()->second != Pauli::Z) + throw std::logic_error( + "Unexpected error during post-selection identification in " + "ChoiMixTableau synthesis"); + Qubit post_selected_qb = row.first.string.begin()->first; + // Multiply other rows to remove Z_qb components + unsigned qb_col = tab.col_index_.left.at(ChoiMixTableau::col_key_t{ + post_selected_qb, ChoiMixTableau::TableauSegment::Input}); + for (unsigned s = 0; s < final_row; ++s) + if (tab.tab_.zmat(s, qb_col)) tab.tab_.row_mult(final_row, s); + // Post-select on correct phase + if (row.second.is_real_negative()) + in_circ.add_op(OpType::X, {post_selected_qb}); + tab.remove_row(final_row); + post_selected.insert(post_selected_qb); + tab.discard_qubit(post_selected_qb, ChoiMixTableau::TableauSegment::Input); + } +} - // Fix phases of zero_initialised qubits - std::set zero_initialised; - Circuit out_circ(output_qubits, {}); - for (unsigned r = 0; r < tab.get_n_rows(); ++r) { +void ChoiMixBuilder::solve_initialised_subspace() { + // Input-first gaussian elimination now sorts the remaining rows into the + // collapsed subspace followed by the zero-initialised subspace and the + // collapsed subspace rows are in upper echelon form over the inputs, giving + // unique leading columns and allowing us to solve them with CXs by + // column-wise gaussian elimination; same for zero-initialised rows over the + // outputs + tab.canonical_column_order(ChoiMixTableau::TableauSegment::Input); + tab.gaussian_form(); + + // Reduce the zero-initialised space to a minimal set of qubits using CX gates + unsigned n_collapsed = 0; + for (; n_collapsed < tab.get_n_rows(); ++n_collapsed) { + if (tab.get_row(n_collapsed).first.string.empty()) break; + } + unsigned n_ins = tab.get_n_inputs(); + unsigned n_outs = tab.get_n_outputs(); + MatrixXb subtableau = + tab.tab_.zmat.bottomRightCorner(tab.get_n_rows() - n_collapsed, n_outs); + std::vector> col_ops = + leading_column_gaussian_col_ops(subtableau); + for (const std::pair& op : col_ops) { + unsigned tab_ctrl_col = n_ins + op.second; + unsigned tab_trgt_col = n_ins + op.first; + tab.tab_.apply_CX(tab_ctrl_col, tab_trgt_col); + ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(tab_ctrl_col); + ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(tab_trgt_col); + out_circ_tp.add_op(OpType::CX, {ctrl.first, trgt.first}); + } + // Initialise rows + for (unsigned r = tab.get_n_rows(); r-- > n_collapsed;) { + // r always refers to the final row in the tableau ChoiMixTableau::row_tensor_t row = tab.get_row(r); if (row.first.size() != 0 || row.second.size() != 1 || row.second.string.begin()->second != Pauli::Z) throw std::logic_error( - "Unexpected error during zero initialisation in ChoiMixTableau " - "synthesis"); + "Unexpected error during initialisation identification in " + "ChoiMixTableau synthesis"); Qubit initialised_qb = row.second.string.begin()->first; - out_circ.qubit_create(initialised_qb); - if (row.second.is_real_negative()) { - out_circ.add_op(OpType::X, {initialised_qb}); - } + // Multiply other rows to remove Z_qb components + unsigned qb_col = tab.col_index_.left.at(ChoiMixTableau::col_key_t{ + initialised_qb, ChoiMixTableau::TableauSegment::Output}); + for (unsigned s = 0; s < r; ++s) + if (tab.tab_.zmat(s, qb_col)) tab.tab_.row_mult(r, s); + // Initialise with correct phase + if (row.second.is_real_negative()) + out_circ_tp.add_op(OpType::X, {initialised_qb}); + tab.remove_row(r); zero_initialised.insert(initialised_qb); + tab.discard_qubit(initialised_qb, ChoiMixTableau::TableauSegment::Output); } +} - // Remaining outputs that aren't zero initialised or matched need to be - // initialised in the maximally-mixed state. - // Also match up unmatched outputs to either unmatched inputs or reusable - // output names (ones that are already matched up to other input names), - // preferring the qubit of the same name - std::list reusable_names; - for (const Qubit& out_qb : output_qubits) { - if (matched_qubits.right.find(out_qb) != matched_qubits.right.end() && - matched_qubits.left.find(out_qb) == matched_qubits.left.end()) - reusable_names.push_back(out_qb); - } - for (const Qubit& out_qb : output_qubits) { - if (matched_qubits.right.find(out_qb) == matched_qubits.right.end()) { - if (zero_initialised.find(out_qb) == zero_initialised.end()) { - out_circ.qubit_create(out_qb); - out_circ.add_op(OpType::H, {out_qb}); - out_circ.add_op(OpType::Collapse, {out_qb}); +void ChoiMixBuilder::solve_collapsed_subspace() { + // Solving the initialised subspace will have preserved the upper echelon form + // of the collapsed subspace; reduce the inputs of the collapsed space to a + // minimal set of qubits using CX gates + unsigned n_ins = tab.get_n_inputs(); + unsigned n_outs = tab.get_n_outputs(); + MatrixXb subtableau = tab.tab_.zmat.topLeftCorner(tab.get_n_rows(), n_ins); + std::vector> col_ops = + leading_column_gaussian_col_ops(subtableau); + for (const std::pair& op : col_ops) { + tab.tab_.apply_CX(op.second, op.first); + ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(op.second); + ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(op.first); + in_circ.add_op(OpType::CX, {ctrl.first, trgt.first}); + } + // Since row multiplications will unsolve the inputs, we cannot get the output + // segment into upper echelon form for the same CX-saving trick; instead we + // accept just removing any qubits that are now unused after solving the + // initialised subspace + remove_unused_qubits(); + tab.canonical_column_order(ChoiMixTableau::TableauSegment::Input); + // Solve the output segment using CX gates + n_ins = tab.get_n_inputs(); + n_outs = tab.get_n_outputs(); + col_ops = gaussian_elimination_col_ops( + tab.tab_.zmat.topRightCorner(tab.get_n_rows(), n_outs)); + for (const std::pair& op : col_ops) { + tab.tab_.apply_CX(n_ins + op.second, n_ins + op.first); + ChoiMixTableau::col_key_t ctrl = tab.col_index_.right.at(n_ins + op.second); + ChoiMixTableau::col_key_t trgt = tab.col_index_.right.at(n_ins + op.first); + out_circ_tp.add_op(OpType::CX, {ctrl.first, trgt.first}); + } + // Connect up and remove rows and columns + for (unsigned r = tab.get_n_rows(); r-- > 0;) { + // r refers to the final row + // Check that row r has been successfully reduced + ChoiMixTableau::row_tensor_t row_r = tab.get_row(r); + if (row_r.first.size() != 1 || + row_r.first.string.begin()->second != Pauli::Z || + row_r.second.size() != 1 || + row_r.second.string.begin()->second != Pauli::Z) + throw std::logic_error( + "Unexpected error during collapsed subspace reduction in " + "ChoiMixTableau synthesis"); + Qubit in_q = row_r.first.string.begin()->first; + Qubit out_q = row_r.second.string.begin()->first; + // Solve phase + if (row_r.second.is_real_negative()) { + in_circ.add_op(OpType::X, {in_q}); + tab.apply_gate(OpType::X, {in_q}, ChoiMixTableau::TableauSegment::Input); + } + // Connect in permutation + in_out_permutation.insert({in_q, out_q}); + collapsed.insert(in_q); + tab.remove_row(r); + tab.discard_qubit(in_q, ChoiMixTableau::TableauSegment::Input); + tab.discard_qubit(out_q, ChoiMixTableau::TableauSegment::Output); + } +} + +void ChoiMixBuilder::remove_unused_qubits() { + // Since removing a column replaces it with the last column, remove in reverse + // order to examine each column exactly once + for (unsigned c = tab.get_n_boundaries(); c-- > 0;) { + bool used = false; + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + if (tab.tab_.zmat(r, c) || tab.tab_.xmat(r, c)) { + used = true; + break; } - if (matched_qubits.left.find(out_qb) == matched_qubits.left.end()) { - matched_qubits.insert({out_qb, out_qb}); - join_permutation.insert({out_qb, out_qb}); - } else { - // Since the matching is a bijection, matched_qubits.left - - // matched_qubits.right (set difference) is the same size as - // matched_qubits.right - matched_qubits.left, so there will be exactly - // the right number of reusable names to pull from here - Qubit name = reusable_names.front(); - reusable_names.pop_front(); - matched_qubits.insert({name, out_qb}); - join_permutation.insert({out_qb, name}); + } + if (used) continue; + ChoiMixTableau::col_key_t col = tab.col_index_.right.at(c); + if (col.second == ChoiMixTableau::TableauSegment::Input) + discarded.insert(col.first); + else + mix_initialised.insert(col.first); + tab.discard_qubit(col.first, col.second); + } +} + +void ChoiMixBuilder::assign_init_post_names() { + auto it = unitary_post_names.begin(); + for (const Qubit& ps : post_selected) { + if (it == unitary_post_names.end()) + throw std::logic_error( + "Not enough additional qubit names for unitary extension of " + "ChoiMixTableau to safely handle post-selected subspace"); + in_out_permutation.insert({ps, *it}); + ++it; + } + unitary_post_names = {it, unitary_post_names.end()}; + + it = unitary_init_names.begin(); + for (const Qubit& zi : zero_initialised) { + if (it == unitary_init_names.end()) + throw std::logic_error( + "Not enough additional qubit names for unitary extension of " + "ChoiMixTableau to safely handle initialised subspace"); + in_out_permutation.insert({*it, zi}); + ++it; + } + unitary_init_names = {it, unitary_init_names.end()}; +} + +void ChoiMixBuilder::assign_remaining_names() { + // Some post-selected or initialised qubits might have already been matched up + // for unitary synthesis, so we only need to match up the remainder + std::set unsolved_ins = discarded; + for (const Qubit& q : post_selected) { + if (in_out_permutation.left.find(q) == in_out_permutation.left.end()) + unsolved_ins.insert(q); + } + std::set unsolved_outs = mix_initialised; + for (const Qubit& q : zero_initialised) { + if (in_out_permutation.right.find(q) == in_out_permutation.right.end()) + unsolved_outs.insert(q); + } + // If there are more unsolved_ins than unsolved_outs, we want to pad out + // unsolved_outs with extra names that don't appear as output names of the + // original tableau; between unsolved_ins and the inputs already in + // in_out_permutation, there will be at least enough of these + if (unsolved_ins.size() > unsolved_outs.size()) { + for (const Qubit& q : unsolved_ins) { + if (in_out_permutation.right.find(q) == in_out_permutation.right.end()) { + unsolved_outs.insert(q); + if (unsolved_ins.size() == unsolved_outs.size()) break; + } + } + if (unsolved_ins.size() > unsolved_outs.size()) { + using perm_entry = boost::bimap::left_const_reference; + BOOST_FOREACH (perm_entry entry, in_out_permutation.left) { + if (in_out_permutation.right.find(entry.first) == + in_out_permutation.right.end()) { + unsolved_outs.insert(entry.first); + if (unsolved_ins.size() == unsolved_outs.size()) break; + } + } + } + } else if (unsolved_ins.size() < unsolved_outs.size()) { + for (const Qubit& q : unsolved_outs) { + if (in_out_permutation.left.find(q) == in_out_permutation.left.end()) { + unsolved_ins.insert(q); + if (unsolved_ins.size() == unsolved_outs.size()) break; + } + } + if (unsolved_ins.size() < unsolved_outs.size()) { + using perm_entry = boost::bimap::left_const_reference; + BOOST_FOREACH (perm_entry entry, in_out_permutation.left) { + if (in_out_permutation.left.find(entry.second) == + in_out_permutation.left.end()) { + unsolved_ins.insert(entry.second); + if (unsolved_ins.size() == unsolved_outs.size()) break; + } } } } + // Prefer to connect qubits with the same names + for (auto in_it = unsolved_ins.begin(); in_it != unsolved_ins.end();) { + auto temp_it = in_it++; + auto out_it = unsolved_outs.find(*temp_it); + if (out_it != unsolved_outs.end()) { + in_out_permutation.insert({*temp_it, *temp_it}); + unsolved_ins.erase(temp_it); + unsolved_outs.erase(out_it); + } + } + // Pair up remainders; by our earlier padding, they should have the exact same + // number of elements, so pair them up exactly + for (const Qubit& in : unsolved_ins) { + auto it = unsolved_outs.begin(); + in_out_permutation.insert({in, *it}); + unsolved_outs.erase(it); + } +} - // Initialise qubits with stabilizer rows and stitch subcircuits together +std::pair ChoiMixBuilder::output_circuit() { + if (tab.get_n_rows() != 0 || tab.get_n_boundaries() != 0) + throw std::logic_error( + "Unexpected error during ChoiMixTableau synthesis, reached the end " + "with a non-empty tableau remaining"); + if (!post_selected.empty()) { + throw std::logic_error( + "Not yet implemented: post-selection required during ChoiMixTableau " + "synthesis"); + } + for (const Qubit& q : discarded) in_circ.qubit_discard(q); + for (const Qubit& q : collapsed) in_circ.add_op(OpType::Collapse, {q}); + Circuit out_circ(out_circ_tp.all_qubits(), {}); + for (const Qubit& q : zero_initialised) out_circ.qubit_create(q); + for (const Qubit& q : mix_initialised) { + out_circ.qubit_create(q); + out_circ.add_op(OpType::H, {q}); + out_circ.add_op(OpType::Collapse, {q}); + } out_circ.append(out_circ_tp.transpose()); - in_circ.append_with_map(out_circ, join_permutation); - return {in_circ, join_permutation}; + qubit_map_t return_perm; + unit_map_t append_perm; + using perm_entry = boost::bimap::left_const_reference; + BOOST_FOREACH (perm_entry entry, in_out_permutation.left) { + return_perm.insert({entry.second, entry.first}); + append_perm.insert({entry.second, entry.first}); + } + in_circ.append_with_map(out_circ, append_perm); + return {in_circ, return_perm}; +} + +std::pair ChoiMixBuilder::unitary_output_circuit() { + if (tab.get_n_rows() != 0 || tab.get_n_boundaries() != 0) + throw std::logic_error( + "Unexpected error during ChoiMixTableau synthesis, reached the end " + "with a non-empty tableau remaining"); + qubit_map_t return_perm; + unit_map_t append_perm; + using perm_entry = boost::bimap::left_const_reference; + BOOST_FOREACH (perm_entry entry, in_out_permutation.left) { + return_perm.insert({entry.second, entry.first}); + append_perm.insert({entry.second, entry.first}); + } + in_circ.append_with_map(out_circ_tp.transpose(), append_perm); + return {in_circ, return_perm}; } } // namespace tket diff --git a/tket/src/Converters/PauliGadget.cpp b/tket/src/Converters/PauliGadget.cpp deleted file mode 100644 index 941f90ea9b..0000000000 --- a/tket/src/Converters/PauliGadget.cpp +++ /dev/null @@ -1,381 +0,0 @@ -// Copyright 2019-2024 Cambridge Quantum 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. - -#include "tket/Converters/PauliGadget.hpp" - -#include "tket/Circuit/CircUtils.hpp" -#include "tket/Circuit/ConjugationBox.hpp" -#include "tket/Circuit/PauliExpBoxes.hpp" - -namespace tket { - -void append_single_pauli_gadget( - Circuit &circ, const SpSymPauliTensor &pauli, CXConfigType cx_config) { - std::vector string; - unit_map_t mapping; - unsigned i = 0; - for (const std::pair &term : pauli.string) { - string.push_back(term.second); - mapping.insert({Qubit(q_default_reg(), i), term.first}); - i++; - } - Circuit gadget = pauli_gadget(string, pauli.coeff, cx_config); - circ.append_with_map(gadget, mapping); -} - -void append_single_pauli_gadget_as_pauli_exp_box( - Circuit &circ, const SpSymPauliTensor &pauli, CXConfigType cx_config) { - std::vector string; - std::vector mapping; - for (const std::pair &term : pauli.string) { - string.push_back(term.second); - mapping.push_back(term.first); - } - PauliExpBox box(SymPauliTensor(string, pauli.coeff), cx_config); - circ.add_box(box, mapping); -} - -void append_pauli_gadget_pair_as_box( - Circuit &circ, const SpSymPauliTensor &pauli0, - const SpSymPauliTensor &pauli1, CXConfigType cx_config) { - std::vector mapping; - std::vector paulis0; - std::vector paulis1; - QubitPauliMap p1map = pauli1.string; - // add paulis for qubits in pauli0_string - for (const std::pair &term : pauli0.string) { - mapping.push_back(term.first); - paulis0.push_back(term.second); - auto found = p1map.find(term.first); - if (found == p1map.end()) { - paulis1.push_back(Pauli::I); - } else { - paulis1.push_back(found->second); - p1map.erase(found); - } - } - // add paulis for qubits in pauli1_string that weren't in pauli0_string - for (const std::pair &term : p1map) { - mapping.push_back(term.first); - paulis1.push_back(term.second); - paulis0.push_back(Pauli::I); // If pauli0_string contained qubit, would - // have been handled above - } - PauliExpPairBox box( - SymPauliTensor(paulis0, pauli0.coeff), - SymPauliTensor(paulis1, pauli1.coeff), cx_config); - circ.add_box(box, mapping); -} - -void append_commuting_pauli_gadget_set_as_box( - Circuit &circ, const std::list &gadgets, - CXConfigType cx_config) { - // Translate SpSymPauliTensors to vectors of Paulis of same length - // Preserves ordering of qubits - - std::set all_qubits; - for (const SpSymPauliTensor &gadget : gadgets) { - for (const std::pair &qubit_pauli : gadget.string) { - all_qubits.insert(qubit_pauli.first); - } - } - - std::vector mapping; - for (const Qubit &qubit : all_qubits) { - mapping.push_back(qubit); - } - - std::vector pauli_gadgets; - for (const SpSymPauliTensor &gadget : gadgets) { - SymPauliTensor &new_gadget = - pauli_gadgets.emplace_back(DensePauliMap{}, gadget.coeff); - for (const Qubit &qubit : mapping) { - new_gadget.string.push_back(gadget.get(qubit)); - } - } - - PauliExpCommutingSetBox box(pauli_gadgets, cx_config); - circ.add_box(box, mapping); -} - -static void reduce_shared_qs_by_CX_snake( - Circuit &circ, std::set &match, SpSymPauliTensor &pauli0, - SpSymPauliTensor &pauli1) { - unsigned match_size = match.size(); - while (match_size > 1) { // We allow one match left over - auto it = --match.end(); - Qubit to_eliminate = *it; - match.erase(it); - Qubit helper = *match.rbegin(); - // extend CX snake - circ.add_op(OpType::CX, {to_eliminate, helper}); - pauli0.string.erase(to_eliminate); - pauli1.string.erase(to_eliminate); - match_size--; - } -} - -static void reduce_shared_qs_by_CX_star( - Circuit &circ, std::set &match, SpSymPauliTensor &pauli0, - SpSymPauliTensor &pauli1) { - std::set::iterator iter = match.begin(); - for (std::set::iterator next = match.begin(); match.size() > 1; - iter = next) { - ++next; - Qubit to_eliminate = *iter; - circ.add_op(OpType::CX, {to_eliminate, *match.rbegin()}); - pauli0.string.erase(to_eliminate); - pauli1.string.erase(to_eliminate); - match.erase(iter); - } -} - -static void reduce_shared_qs_by_CX_tree( - Circuit &circ, std::set &match, SpSymPauliTensor &pauli0, - SpSymPauliTensor &pauli1) { - while (match.size() > 1) { - std::set remaining; - std::set::iterator it = match.begin(); - while (it != match.end()) { - Qubit maintained = *it; - it++; - remaining.insert(maintained); - if (it != match.end()) { - Qubit to_eliminate = *it; - it++; - circ.add_op(OpType::CX, {to_eliminate, maintained}); - pauli0.string.erase(to_eliminate); - pauli1.string.erase(to_eliminate); - } - } - match = remaining; - } -} - -static void reduce_shared_qs_by_CX_multiqgate( - Circuit &circ, std::set &match, SpSymPauliTensor &pauli0, - SpSymPauliTensor &pauli1) { - if (match.size() <= 1) { - return; - } - // last qubit is target - Qubit target = *match.rbegin(); - while (match.size() > 1) { - std::set::iterator iter = match.begin(); - if (match.size() == 2) { - // use CX - Qubit to_eliminate = *iter; - match.erase(iter); - pauli0.string.erase(to_eliminate); - pauli1.string.erase(to_eliminate); - - circ.add_op(OpType::CX, {to_eliminate, target}); - } else { - // use XXPhase3 - Qubit to_eliminate1 = *iter; - match.erase(iter++); - pauli0.string.erase(to_eliminate1); - pauli1.string.erase(to_eliminate1); - - Qubit to_eliminate2 = *iter; - match.erase(iter); - pauli0.string.erase(to_eliminate2); - pauli1.string.erase(to_eliminate2); - - circ.add_op(OpType::H, {to_eliminate1}); - circ.add_op(OpType::H, {to_eliminate2}); - circ.add_op( - OpType::XXPhase3, 0.5, {to_eliminate1, to_eliminate2, target}); - circ.add_op(OpType::X, {target}); - } - } -} - -void append_pauli_gadget_pair( - Circuit &circ, SpSymPauliTensor pauli0, SpSymPauliTensor pauli1, - CXConfigType cx_config) { - /* - * Cowtan, Dilkes, Duncan, Simmons, Sivarajah: Phase Gadget Synthesis for - * Shallow Circuits, Lemma 4.9 - * Let s and t be Pauli strings; then there exists a Clifford unitary U such - * that - * P(a, s) . P(b, t) = U . P(a, s') . P(b, t') . U^\dagger - * where s' and t' are Pauli strings with intersection at most 1. - * - * Follows the procedure to reduce the intersection of the gadgets and then - * synthesises the remainder individually. - */ - pauli0.compress(); - pauli1.compress(); - - /* - * Step 1: Partition qubits into those just affected by pauli0 (just0) and - * pauli1 (just1), and those in both which either match or don't - */ - std::set just0 = pauli0.own_qubits(pauli1); - std::set just1 = pauli1.own_qubits(pauli0); - std::set match = pauli0.common_qubits(pauli1); - std::set mismatch = pauli0.conflicting_qubits(pauli1); - - /* - * Step 2: Build the unitary U that minimises the intersection of the gadgets. - */ - Circuit u; - for (const Qubit &qb : just0) u.add_qubit(qb); - for (const Qubit &qb : just1) u.add_qubit(qb); - for (const Qubit &qb : match) u.add_qubit(qb); - for (const Qubit &qb : mismatch) u.add_qubit(qb); - Circuit v(u); - - /* - * Step 2.i: Remove (almost) all matches by converting to Z basis and applying - * CXs - */ - for (const Qubit &qb : match) { - switch (pauli0.get(qb)) { - case Pauli::X: - u.add_op(OpType::H, {qb}); - pauli0.set(qb, Pauli::Z); - pauli1.set(qb, Pauli::Z); - break; - case Pauli::Y: - u.add_op(OpType::V, {qb}); - pauli0.set(qb, Pauli::Z); - pauli1.set(qb, Pauli::Z); - break; - default: - break; - } - } - switch (cx_config) { - case CXConfigType::Snake: { - reduce_shared_qs_by_CX_snake(u, match, pauli0, pauli1); - break; - } - case CXConfigType::Star: { - reduce_shared_qs_by_CX_star(u, match, pauli0, pauli1); - break; - } - case CXConfigType::Tree: { - reduce_shared_qs_by_CX_tree(u, match, pauli0, pauli1); - break; - } - case CXConfigType::MultiQGate: { - reduce_shared_qs_by_CX_multiqgate(u, match, pauli0, pauli1); - break; - } - default: - throw std::logic_error( - "Unknown CXConfigType received when decomposing gadget."); - } - /* - * Step 2.ii: Convert mismatches to Z in pauli0 and X in pauli1 - */ - for (const Qubit &qb : mismatch) { - switch (pauli0.get(qb)) { - case Pauli::X: { - switch (pauli1.get(qb)) { - case Pauli::Y: - u.add_op(OpType::Sdg, {qb}); - u.add_op(OpType::Vdg, {qb}); - break; - case Pauli::Z: - u.add_op(OpType::H, {qb}); - break; - default: - break; // Cannot hit this case - } - break; - } - case Pauli::Y: { - switch (pauli1.get(qb)) { - case Pauli::X: - u.add_op(OpType::V, {qb}); - break; - case Pauli::Z: - u.add_op(OpType::V, {qb}); - u.add_op(OpType::S, {qb}); - break; - default: - break; // Cannot hit this case - } - break; - } - default: { // Necessarily Z - if (pauli1.get(qb) == Pauli::Y) u.add_op(OpType::Sdg, {qb}); - // No need to act if already X - } - } - pauli0.set(qb, Pauli::Z); - pauli1.set(qb, Pauli::X); - } - - /* - * Step 2.iii: Remove the final matching qubit against a mismatch if one - * exists, otherwise allow both gadgets to build it - */ - if (!match.empty()) { - Qubit last_match = *match.begin(); - match.erase(last_match); - if (!mismatch.empty()) { - Qubit mismatch_used = - *mismatch.rbegin(); // Prefer to use the one that may be left over - // after reducing pairs - u.add_op(OpType::S, {mismatch_used}); - u.add_op(OpType::CX, {last_match, mismatch_used}); - u.add_op(OpType::Sdg, {mismatch_used}); - pauli0.string.erase(last_match); - pauli1.string.erase(last_match); - } else { - just0.insert(last_match); - just1.insert(last_match); - } - } - - /* - * Step 2.iv: Reduce pairs of mismatches to different qubits. - * Allow both gadgets to build a remaining qubit if it exists. - */ - std::set::iterator mis_it = mismatch.begin(); - while (mis_it != mismatch.end()) { - Qubit z_in_0 = *mis_it; - just0.insert(z_in_0); - mis_it++; - if (mis_it == mismatch.end()) { - just1.insert(z_in_0); - } else { - Qubit x_in_1 = *mis_it; - u.add_op(OpType::CX, {x_in_1, z_in_0}); - pauli0.string.erase(x_in_1); - pauli1.string.erase(z_in_0); - just1.insert(x_in_1); - mis_it++; - } - } - - /* - * Step 3: Combine circuits to give final result - */ - append_single_pauli_gadget(v, pauli0); - append_single_pauli_gadget(v, pauli1); - // ConjugationBox components must be in the default register - qubit_vector_t all_qubits = u.all_qubits(); - u.flatten_registers(); - v.flatten_registers(); - ConjugationBox cjbox( - std::make_shared(u), std::make_shared(v)); - circ.add_box(cjbox, all_qubits); -} - -} // namespace tket diff --git a/tket/src/Converters/PauliGraphConverters.cpp b/tket/src/Converters/PauliGraphConverters.cpp index e04d58ff55..46a3126433 100644 --- a/tket/src/Converters/PauliGraphConverters.cpp +++ b/tket/src/Converters/PauliGraphConverters.cpp @@ -15,7 +15,6 @@ #include "tket/Circuit/Boxes.hpp" #include "tket/Circuit/PauliExpBoxes.hpp" #include "tket/Converters/Converters.hpp" -#include "tket/Converters/PauliGadget.hpp" #include "tket/Converters/PhasePoly.hpp" #include "tket/Diagonalisation/Diagonalisation.hpp" #include "tket/Gate/Gate.hpp" diff --git a/tket/src/Converters/UnitaryTableauConverters.cpp b/tket/src/Converters/UnitaryTableauConverters.cpp index 22ad4dfbef..8631cfa1a3 100644 --- a/tket/src/Converters/UnitaryTableauConverters.cpp +++ b/tket/src/Converters/UnitaryTableauConverters.cpp @@ -44,7 +44,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * Step 1: Use Hadamards (in our case, Vs) to make C (z rows of xmat_) have * full rank */ - MatrixXb echelon = tabl.xmat_.block(size, 0, size, size); + MatrixXb echelon = tabl.xmat.block(size, 0, size, size); std::map leading_val_to_col; for (unsigned i = 0; i < size; i++) { for (unsigned j = 0; j < size; j++) { @@ -64,9 +64,8 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { continue; // Independent of previous cols c.add_op(OpType::V, {i}); tabl.apply_V(i); - tabl.apply_V(i); - tabl.apply_V(i); - echelon.col(i) = tabl.zmat_.block(size, i, size, 1); + tabl.apply_X(i); + echelon.col(i) = tabl.zmat.block(size, i, size, 1); for (unsigned j = 0; j < size; j++) { if (echelon(j, i)) { if (leading_val_to_col.find(j) == leading_val_to_col.end()) { @@ -89,7 +88,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * / A B \ * \ I D / */ - MatrixXb to_reduce = tabl.xmat_.block(size, 0, size, size); + MatrixXb to_reduce = tabl.xmat.block(size, 0, size, size); for (const std::pair& qbs : gaussian_elimination_col_ops(to_reduce)) { c.add_op(OpType::CX, {qbs.first, qbs.second}); @@ -103,13 +102,12 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * for some invertible M. */ std::pair zp_z_llt = - binary_LLT_decomposition(tabl.zmat_.block(size, 0, size, size)); + binary_LLT_decomposition(tabl.zmat.block(size, 0, size, size)); for (unsigned i = 0; i < size; i++) { if (zp_z_llt.second(i, i)) { c.add_op(OpType::S, {i}); tabl.apply_S(i); - tabl.apply_S(i); - tabl.apply_S(i); + tabl.apply_Z(i); } } @@ -140,8 +138,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { for (unsigned i = 0; i < size; i++) { c.add_op(OpType::S, {i}); tabl.apply_S(i); - tabl.apply_S(i); - tabl.apply_S(i); + tabl.apply_Z(i); } /* @@ -150,7 +147,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * \ I 0 / * By commutativity relations, IB^T = A0^T + I, therefore B = I. */ - to_reduce = tabl.xmat_.block(size, 0, size, size); + to_reduce = tabl.xmat.block(size, 0, size, size); for (const std::pair& qbs : gaussian_elimination_col_ops(to_reduce)) { c.add_op(OpType::CX, {qbs.first, qbs.second}); @@ -164,9 +161,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { */ for (unsigned i = 0; i < size; i++) { c.add_op(OpType::H, {i}); - tabl.apply_S(i); - tabl.apply_V(i); - tabl.apply_S(i); + tabl.apply_H(i); } /* @@ -175,13 +170,12 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * some invertible N. */ std::pair xp_z_llt = - binary_LLT_decomposition(tabl.zmat_.block(0, 0, size, size)); + binary_LLT_decomposition(tabl.zmat.block(0, 0, size, size)); for (unsigned i = 0; i < size; i++) { if (xp_z_llt.second(i, i)) { c.add_op(OpType::S, {i}); tabl.apply_S(i); - tabl.apply_S(i); - tabl.apply_S(i); + tabl.apply_Z(i); } } @@ -210,8 +204,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { for (unsigned i = 0; i < size; i++) { c.add_op(OpType::S, {i}); tabl.apply_S(i); - tabl.apply_S(i); - tabl.apply_S(i); + tabl.apply_Z(i); } /* @@ -220,7 +213,7 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * \ 0 I / */ for (const std::pair& qbs : - gaussian_elimination_col_ops(tabl.xmat_.block(0, 0, size, size))) { + gaussian_elimination_col_ops(tabl.xmat.block(0, 0, size, size))) { c.add_op(OpType::CX, {qbs.first, qbs.second}); tabl.apply_CX(qbs.first, qbs.second); } @@ -229,15 +222,13 @@ Circuit unitary_tableau_to_circuit(const UnitaryTableau& tab) { * DELAYED STEPS: Set all phases to 0 by applying Z or X gates */ for (unsigned i = 0; i < size; i++) { - if (tabl.phase_(i)) { + if (tabl.phase(i)) { c.add_op(OpType::Z, {i}); - tabl.apply_S(i); - tabl.apply_S(i); + tabl.apply_Z(i); } - if (tabl.phase_(i + size)) { + if (tabl.phase(i + size)) { c.add_op(OpType::X, {i}); - tabl.apply_V(i); - tabl.apply_V(i); + tabl.apply_X(i); } } diff --git a/tket/src/Diagonalisation/Diagonalisation.cpp b/tket/src/Diagonalisation/Diagonalisation.cpp index 7ab9d876db..6b433e10a3 100644 --- a/tket/src/Diagonalisation/Diagonalisation.cpp +++ b/tket/src/Diagonalisation/Diagonalisation.cpp @@ -342,4 +342,413 @@ void apply_conjugations( qps.coeff *= cast_coeff(stab.coeff); } +std::pair reduce_pauli_to_z( + const SpPauliStabiliser &pauli, CXConfigType cx_config) { + Circuit circ; + qubit_vector_t qubits; + for (const std::pair &qp : pauli.string) { + circ.add_qubit(qp.first); + if (qp.second != Pauli::I) qubits.push_back(qp.first); + switch (qp.second) { + case Pauli::X: { + circ.add_op(OpType::H, {qp.first}); + break; + } + case Pauli::Y: { + circ.add_op(OpType::V, {qp.first}); + break; + } + default: { + break; + } + } + } + unsigned n_qubits = qubits.size(); + if (n_qubits == 0) throw std::logic_error("Cannot reduce identity to Z"); + switch (cx_config) { + case CXConfigType::Snake: { + for (unsigned i = n_qubits - 1; i != 0; --i) { + circ.add_op(OpType::CX, {qubits.at(i), qubits.at(i - 1)}); + } + break; + } + case CXConfigType::Star: { + for (unsigned i = n_qubits - 1; i != 0; --i) { + circ.add_op(OpType::CX, {qubits.at(i), qubits.front()}); + } + break; + } + case CXConfigType::Tree: { + for (unsigned step_size = 1; step_size < n_qubits; step_size *= 2) { + for (unsigned i = 0; step_size + i < n_qubits; i += 2 * step_size) { + circ.add_op( + OpType::CX, {qubits.at(step_size + i), qubits.at(i)}); + } + } + break; + } + case CXConfigType::MultiQGate: { + bool flip_phase = false; + for (unsigned i = n_qubits - 1; i != 0; --i) { + if (i == 1) { + circ.add_op(OpType::CX, {qubits.at(i), qubits.front()}); + } else { + /** + * This is only equal to the CX decompositions above up to phase, + * but phase differences are cancelled out by its dagger + */ + circ.add_op(OpType::H, {qubits.at(i)}); + circ.add_op(OpType::H, {qubits.at(i - 1)}); + circ.add_op( + OpType::XXPhase3, 0.5, + {qubits.at(i), qubits.at(i - 1), qubits.front()}); + --i; + flip_phase = !flip_phase; + } + } + if (flip_phase) circ.add_op(OpType::X, {qubits.front()}); + break; + } + } + return {circ, qubits.front()}; +} + +static void reduce_shared_qs_by_CX_snake( + Circuit &circ, std::set &match, SpPauliStabiliser &pauli0, + SpPauliStabiliser &pauli1) { + unsigned match_size = match.size(); + while (match_size > 1) { // We allow one match left over + auto it = --match.end(); + Qubit to_eliminate = *it; + match.erase(it); + Qubit helper = *match.rbegin(); + // extend CX snake + circ.add_op(OpType::CX, {to_eliminate, helper}); + pauli0.string.erase(to_eliminate); + pauli1.string.erase(to_eliminate); + match_size--; + } +} + +static void reduce_shared_qs_by_CX_star( + Circuit &circ, std::set &match, SpPauliStabiliser &pauli0, + SpPauliStabiliser &pauli1) { + std::set::iterator iter = match.begin(); + for (std::set::iterator next = match.begin(); match.size() > 1; + iter = next) { + ++next; + Qubit to_eliminate = *iter; + circ.add_op(OpType::CX, {to_eliminate, *match.rbegin()}); + pauli0.string.erase(to_eliminate); + pauli1.string.erase(to_eliminate); + match.erase(iter); + } +} + +static void reduce_shared_qs_by_CX_tree( + Circuit &circ, std::set &match, SpPauliStabiliser &pauli0, + SpPauliStabiliser &pauli1) { + while (match.size() > 1) { + std::set remaining; + std::set::iterator it = match.begin(); + while (it != match.end()) { + Qubit maintained = *it; + it++; + remaining.insert(maintained); + if (it != match.end()) { + Qubit to_eliminate = *it; + it++; + circ.add_op(OpType::CX, {to_eliminate, maintained}); + pauli0.string.erase(to_eliminate); + pauli1.string.erase(to_eliminate); + } + } + match = remaining; + } +} + +static void reduce_shared_qs_by_CX_multiqgate( + Circuit &circ, std::set &match, SpPauliStabiliser &pauli0, + SpPauliStabiliser &pauli1) { + if (match.size() <= 1) { + return; + } + // last qubit is target + Qubit target = *match.rbegin(); + while (match.size() > 1) { + std::set::iterator iter = match.begin(); + if (match.size() == 2) { + // use CX + Qubit to_eliminate = *iter; + match.erase(iter); + pauli0.string.erase(to_eliminate); + pauli1.string.erase(to_eliminate); + + circ.add_op(OpType::CX, {to_eliminate, target}); + } else { + // use XXPhase3 + Qubit to_eliminate1 = *iter; + match.erase(iter++); + pauli0.string.erase(to_eliminate1); + pauli1.string.erase(to_eliminate1); + + Qubit to_eliminate2 = *iter; + match.erase(iter); + pauli0.string.erase(to_eliminate2); + pauli1.string.erase(to_eliminate2); + + circ.add_op(OpType::H, {to_eliminate1}); + circ.add_op(OpType::H, {to_eliminate2}); + circ.add_op( + OpType::XXPhase3, 0.5, {to_eliminate1, to_eliminate2, target}); + circ.add_op(OpType::X, {target}); + } + } +} + +std::pair> reduce_overlap_of_paulis( + SpPauliStabiliser &pauli0, SpPauliStabiliser &pauli1, + CXConfigType cx_config, bool allow_matching_final) { + /* + * Cowtan, Dilkes, Duncan, Simmons, Sivarajah: Phase Gadget Synthesis for + * Shallow Circuits, Lemma 4.9 + * Let s and t be Pauli strings; then there exists a Clifford unitary U such + * that + * P(a, s) . P(b, t) = U . P(a, s') . P(b, t') . U^\dagger + * where s' and t' are Pauli strings with intersection at most 1. + * + * Follows the procedure to reduce the intersection of the gadgets and then + * synthesises the remainder individually. + */ + + /* + * Step 1: Identify qubits in both pauli0 and pauli1 which either match or + * don't + */ + std::set match = pauli0.common_qubits(pauli1); + std::set mismatch = pauli0.conflicting_qubits(pauli1); + + /* + * Step 2: Build the unitary U that minimises the intersection of the gadgets. + */ + Circuit u; + for (const std::pair &qp : pauli0.string) + u.add_qubit(qp.first); + for (const std::pair &qp : pauli1.string) { + if (!u.contains_unit(qp.first)) u.add_qubit(qp.first); + } + + /* + * Step 2.i: Remove (almost) all matches by converting to Z basis and applying + * CXs + */ + for (const Qubit &qb : match) { + switch (pauli0.get(qb)) { + case Pauli::X: + u.add_op(OpType::H, {qb}); + pauli0.set(qb, Pauli::Z); + pauli1.set(qb, Pauli::Z); + break; + case Pauli::Y: + u.add_op(OpType::V, {qb}); + pauli0.set(qb, Pauli::Z); + pauli1.set(qb, Pauli::Z); + break; + default: + break; + } + } + switch (cx_config) { + case CXConfigType::Snake: { + reduce_shared_qs_by_CX_snake(u, match, pauli0, pauli1); + break; + } + case CXConfigType::Star: { + reduce_shared_qs_by_CX_star(u, match, pauli0, pauli1); + break; + } + case CXConfigType::Tree: { + reduce_shared_qs_by_CX_tree(u, match, pauli0, pauli1); + break; + } + case CXConfigType::MultiQGate: { + reduce_shared_qs_by_CX_multiqgate(u, match, pauli0, pauli1); + break; + } + default: + throw std::logic_error( + "Unknown CXConfigType received when decomposing gadget."); + } + /* + * Step 2.ii: Convert mismatches to Z in pauli0 and X in pauli1 + */ + for (const Qubit &qb : mismatch) { + switch (pauli0.get(qb)) { + case Pauli::X: { + switch (pauli1.get(qb)) { + case Pauli::Y: + u.add_op(OpType::Sdg, {qb}); + u.add_op(OpType::Vdg, {qb}); + break; + case Pauli::Z: + u.add_op(OpType::H, {qb}); + break; + default: + TKET_ASSERT(false); + break; // Cannot hit this case + } + break; + } + case Pauli::Y: { + switch (pauli1.get(qb)) { + case Pauli::X: + u.add_op(OpType::V, {qb}); + break; + case Pauli::Z: + u.add_op(OpType::V, {qb}); + u.add_op(OpType::S, {qb}); + break; + default: + TKET_ASSERT(false); + break; // Cannot hit this case + } + break; + } + default: { // Necessarily Z + if (pauli1.get(qb) == Pauli::Y) u.add_op(OpType::Sdg, {qb}); + // No need to act if already X + } + } + pauli0.set(qb, Pauli::Z); + pauli1.set(qb, Pauli::X); + } + + /* + * Step 2.iii: Remove the final matching qubit against a mismatch if one + * exists, otherwise remove into another qubit or allow final overlap to be + * matching + */ + std::optional last_overlap = std::nullopt; + if (!match.empty()) { + Qubit last_match = *match.begin(); + if (!mismatch.empty()) { + Qubit mismatch_used = + *mismatch.rbegin(); // Prefer to use the one that may be left over + // after reducing pairs + u.add_op(OpType::S, {mismatch_used}); + u.add_op(OpType::CX, {last_match, mismatch_used}); + u.add_op(OpType::Sdg, {mismatch_used}); + pauli0.string.erase(last_match); + pauli1.string.erase(last_match); + } else if (!allow_matching_final) { + std::optional> other; + for (const std::pair &qp : pauli0.string) { + if (qp.first != last_match && qp.second != Pauli::I) { + other = qp; + pauli0.string.erase(last_match); + break; + } + } + if (!other) { + for (const std::pair &qp : pauli1.string) { + if (qp.first != last_match && qp.second != Pauli::I) { + other = qp; + pauli1.string.erase(last_match); + break; + } + } + if (!other) + throw std::logic_error( + "Cannot reduce identical Paulis to different qubits"); + } + if (other->second == Pauli::X) { + u.add_op(OpType::H, {other->first}); + u.add_op(OpType::CX, {last_match, other->first}); + u.add_op(OpType::H, {other->first}); + } else { + u.add_op(OpType::CX, {last_match, other->first}); + } + } else { + last_overlap = last_match; + } + } + + /* + * Step 2.iv: Reduce pairs of mismatches to different qubits. + * Allow both gadgets to build a remaining qubit if it exists. + */ + std::set::iterator mis_it = mismatch.begin(); + while (mis_it != mismatch.end()) { + Qubit z_in_0 = *mis_it; + mis_it++; + if (mis_it != mismatch.end()) { + Qubit x_in_1 = *mis_it; + u.add_op(OpType::CX, {x_in_1, z_in_0}); + pauli0.string.erase(x_in_1); + pauli1.string.erase(z_in_0); + mis_it++; + } else { + last_overlap = z_in_0; + } + } + + return {u, last_overlap}; +} + +std::pair reduce_anticommuting_paulis_to_z_x( + SpPauliStabiliser pauli0, SpPauliStabiliser pauli1, + CXConfigType cx_config) { + std::pair> reduced_overlap = + reduce_overlap_of_paulis(pauli0, pauli1, cx_config); + Circuit &u = reduced_overlap.first; + if (!reduced_overlap.second) + throw std::logic_error("No overlap for anti-commuting paulis"); + Qubit &last_mismatch = *reduced_overlap.second; + + /** + * Reduce each remaining Pauli to the shared mismatching qubit. + * Since reduce_pauli_to_Z does not allow us to pick the final qubit, we + * reserve the mismatching qubit, call reduce_pauli_to_Z on the rest, and add + * a CX. + */ + pauli0.string.erase(last_mismatch); + pauli0.compress(); + if (!pauli0.string.empty()) { + std::pair diag0 = reduce_pauli_to_z(pauli0, cx_config); + u.append(diag0.first); + u.add_op(OpType::CX, {diag0.second, last_mismatch}); + } + pauli1.compress(); + pauli1.string.erase(last_mismatch); + if (!pauli1.string.empty()) { + std::pair diag1 = reduce_pauli_to_z(pauli1, cx_config); + u.append(diag1.first); + u.add_op(OpType::H, {last_mismatch}); + u.add_op(OpType::CX, {diag1.second, last_mismatch}); + u.add_op(OpType::H, {last_mismatch}); + } + + return {u, last_mismatch}; +} + +std::tuple reduce_commuting_paulis_to_zi_iz( + SpPauliStabiliser pauli0, SpPauliStabiliser pauli1, + CXConfigType cx_config) { + std::pair> reduced_overlap = + reduce_overlap_of_paulis(pauli0, pauli1, cx_config); + Circuit &u = reduced_overlap.first; + if (reduced_overlap.second) + throw std::logic_error("Overlap remaining for commuting paulis"); + + /** + * Reduce each remaining Pauli to a single qubit. + */ + std::pair diag0 = reduce_pauli_to_z(pauli0, cx_config); + u.append(diag0.first); + std::pair diag1 = reduce_pauli_to_z(pauli1, cx_config); + u.append(diag1.first); + + return {u, diag0.second, diag1.second}; +} + } // namespace tket diff --git a/tket/src/Transformations/PauliOptimisation.cpp b/tket/src/Transformations/PauliOptimisation.cpp index 928fa73c58..eac7732d04 100644 --- a/tket/src/Transformations/PauliOptimisation.cpp +++ b/tket/src/Transformations/PauliOptimisation.cpp @@ -14,8 +14,8 @@ #include "tket/Transformations/PauliOptimisation.hpp" +#include "tket/Circuit/CircUtils.hpp" #include "tket/Converters/Converters.hpp" -#include "tket/Converters/PauliGadget.hpp" #include "tket/OpType/OpType.hpp" #include "tket/OpType/OpTypeInfo.hpp" #include "tket/Ops/Op.hpp" @@ -168,14 +168,14 @@ Transform pairwise_pauli_gadgets(CXConfigType cx_config) { // Synthesise pairs of Pauli Gadgets unsigned g = 0; while (g + 1 < pauli_gadgets.size()) { - append_pauli_gadget_pair( - gadget_circ, pauli_gadgets[g], pauli_gadgets[g + 1], cx_config); + gadget_circ.append( + pauli_gadget_pair(pauli_gadgets[g], pauli_gadgets[g + 1], cx_config)); g += 2; } // As we synthesised Pauli gadgets 2 at a time, if there were an odd // number, we will have one left over, so add that one on its own if (g < pauli_gadgets.size()) { - append_single_pauli_gadget(gadget_circ, pauli_gadgets[g], cx_config); + gadget_circ.append(pauli_gadget(pauli_gadgets[g], cx_config)); } // Stitch gadget circuit and Clifford circuit together circ = gadget_circ >> clifford_circ; diff --git a/tket/test/CMakeLists.txt b/tket/test/CMakeLists.txt index 5b1ccb620d..65cf8824ef 100644 --- a/tket/test/CMakeLists.txt +++ b/tket/test/CMakeLists.txt @@ -124,6 +124,7 @@ add_executable(test-tket src/Circuit/test_DummyBox.cpp src/test_UnitaryTableau.cpp src/test_ChoiMixTableau.cpp + src/test_Diagonalisation.cpp src/test_PhasePolynomials.cpp src/test_PauliGraph.cpp src/test_Architectures.cpp diff --git a/tket/test/src/test_ChoiMixTableau.cpp b/tket/test/src/test_ChoiMixTableau.cpp index 49ad200d54..ac8ae94817 100644 --- a/tket/test/src/test_ChoiMixTableau.cpp +++ b/tket/test/src/test_ChoiMixTableau.cpp @@ -15,6 +15,7 @@ #include #include "testutil.hpp" +#include "tket/Circuit/Simulation/CircuitSimulator.hpp" #include "tket/Converters/Converters.hpp" namespace tket { @@ -60,6 +61,12 @@ static ChoiMixTableau get_tableau_with_gates_applied_at_front() { OpType::CX, {Qubit(0), Qubit(1)}, ChoiMixTableau::TableauSegment::Input); return tab; } +static qubit_map_t inv_perm(const qubit_map_t& perm) { + qubit_map_t inv; + for (const std::pair& qp : perm) + inv.insert({qp.second, qp.first}); + return inv; +} SCENARIO("Correct creation of ChoiMixTableau") { GIVEN( @@ -89,7 +96,7 @@ SCENARIO("Correct creation of ChoiMixTableau") { tab.get_row_product({0, 1}) == ChoiMixTableau::row_tensor_t{ SpPauliStabiliser(Qubit(0), Pauli::Y), - SpPauliStabiliser(Qubit(0), Pauli::Y, 2)}); + SpPauliStabiliser(Qubit(0), Pauli::Y)}); THEN("Serialize and deserialize") { nlohmann::json j_tab = tab; ChoiMixTableau tab2{{}}; @@ -174,12 +181,10 @@ SCENARIO("Correct creation of ChoiMixTableau") { SpPauliStabiliser(Qubit(0), Pauli::Z), SpPauliStabiliser(Qubit(0), Pauli::Z)}); // Affecting the input segment should give the same effect as for - // UnitaryRevTableau (since lhs is transposed, +Y is flipped to -Y, and - // phase is returned on rhs) + // UnitaryRevTableau REQUIRE( - tab.get_row(2) == - ChoiMixTableau::row_tensor_t{ - SpPauliStabiliser(Qubit(1), Pauli::Y), SpPauliStabiliser({}, 2)}); + tab.get_row(2) == ChoiMixTableau::row_tensor_t{ + SpPauliStabiliser(Qubit(1), Pauli::Y), {}}); // Affecting the output segment should give the same effect as for // UnitaryTableau REQUIRE( @@ -197,9 +202,8 @@ SCENARIO("Correct creation of ChoiMixTableau") { SpPauliStabiliser(Qubit(0), Pauli::Z), SpPauliStabiliser(Qubit(0), Pauli::Y, 2)}); REQUIRE( - tab.get_row(2) == - ChoiMixTableau::row_tensor_t{ - SpPauliStabiliser(Qubit(1), Pauli::Y), SpPauliStabiliser({}, 2)}); + tab.get_row(2) == ChoiMixTableau::row_tensor_t{ + SpPauliStabiliser(Qubit(1), Pauli::Y), {}}); REQUIRE( tab.get_row(3) == ChoiMixTableau::row_tensor_t{ {}, SpPauliStabiliser(Qubit(2), Pauli::Y, 2)}); @@ -215,9 +219,8 @@ SCENARIO("Correct creation of ChoiMixTableau") { SpPauliStabiliser(Qubit(0), Pauli::Z), SpPauliStabiliser(Qubit(0), Pauli::Z, 2)}); REQUIRE( - tab.get_row(2) == - ChoiMixTableau::row_tensor_t{ - SpPauliStabiliser(Qubit(1), Pauli::Y), SpPauliStabiliser({}, 2)}); + tab.get_row(2) == ChoiMixTableau::row_tensor_t{ + SpPauliStabiliser(Qubit(1), Pauli::Y), {}}); REQUIRE( tab.get_row(3) == ChoiMixTableau::row_tensor_t{ {}, SpPauliStabiliser(Qubit(2), Pauli::Y, 2)}); @@ -449,17 +452,21 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { GIVEN("An identity circuit") { Circuit circ(3); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; + Circuit res = cm_tableau_to_exact_circuit(tab).first; + REQUIRE(res == circ); + res = cm_tableau_to_unitary_extension_circuit(tab).first; REQUIRE(res == circ); } GIVEN("Just some Pauli gates for phase tests") { Circuit circ(4); circ.add_op(OpType::X, {1}); - circ.add_op(OpType::X, {2}); circ.add_op(OpType::Z, {2}); + circ.add_op(OpType::X, {2}); circ.add_op(OpType::Z, {3}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; + Circuit res = cm_tableau_to_exact_circuit(tab).first; + REQUIRE(res == circ); + res = cm_tableau_to_unitary_extension_circuit(tab).first; REQUIRE(res == circ); } GIVEN("Iterate through single-qubit Cliffords with all entanglements") { @@ -480,7 +487,9 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { if ((i / 9) % 3 == 1) circ.add_op(OpType::S, {0}); if ((i / 9) % 3 == 2) circ.add_op(OpType::Sdg, {0}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; + Circuit res = cm_tableau_to_exact_circuit(tab).first; + Circuit res_uni = cm_tableau_to_unitary_extension_circuit(tab).first; + REQUIRE(res == res_uni); ChoiMixTableau res_tab = circuit_to_cm_tableau(res); REQUIRE(res_tab == tab); REQUIRE(test_unitary_comparison(circ, res, true)); @@ -489,10 +498,15 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { GIVEN("A unitary circuit") { Circuit circ = get_test_circ(); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; - ChoiMixTableau res_tab = circuit_to_cm_tableau(res); + std::pair res = cm_tableau_to_exact_circuit(tab); + res.first.permute_boundary_output(inv_perm(res.second)); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit(tab); + res_uni.first.permute_boundary_output(inv_perm(res_uni.second)); + REQUIRE(res.first == res_uni.first); + ChoiMixTableau res_tab = circuit_to_cm_tableau(res.first); REQUIRE(res_tab == tab); - REQUIRE(test_unitary_comparison(circ, res, true)); + REQUIRE(test_unitary_comparison(circ, res.first, true)); } GIVEN("Check unitary equivalence by calculating matrix") { Circuit circ(4); @@ -501,21 +515,69 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { circ.add_op(OpType::ISWAPMax, {0, 3}); circ.add_op(OpType::SX, {1}); circ.add_op(OpType::SXdg, {2}); + circ.add_op(OpType::CY, {1, 3}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; - REQUIRE(test_unitary_comparison(circ, res, true)); + std::pair res = cm_tableau_to_exact_circuit(tab); + res.first.permute_boundary_output(inv_perm(res.second)); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit(tab); + res_uni.first.permute_boundary_output(inv_perm(res_uni.second)); + REQUIRE(res.first == res_uni.first); + REQUIRE(test_unitary_comparison(circ, res.first, true)); + THEN("Build the tableau manually for apply_gate coverage on inputs") { + ChoiMixTableau rev_tab(4); + rev_tab.apply_gate( + OpType::CY, {Qubit(1), Qubit(3)}, + ChoiMixTableau::TableauSegment::Input); + rev_tab.apply_gate( + OpType::SXdg, {Qubit(2)}, ChoiMixTableau::TableauSegment::Input); + rev_tab.apply_gate( + OpType::SX, {Qubit(1)}, ChoiMixTableau::TableauSegment::Input); + rev_tab.apply_gate( + OpType::ISWAPMax, {Qubit(0), Qubit(3)}, + ChoiMixTableau::TableauSegment::Input); + rev_tab.apply_gate( + OpType::ECR, {Qubit(2), Qubit(3)}, + ChoiMixTableau::TableauSegment::Input); + rev_tab.apply_gate( + OpType::ZZMax, {Qubit(0), Qubit(1)}, + ChoiMixTableau::TableauSegment::Input); + rev_tab.canonical_column_order(); + rev_tab.gaussian_form(); + REQUIRE(tab == rev_tab); + } } GIVEN("A Clifford state") { Circuit circ = get_test_circ(); circ.qubit_create_all(); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; + Circuit res = cm_tableau_to_exact_circuit(tab).first; + ChoiMixTableau res_tab = circuit_to_cm_tableau(res); + REQUIRE(res_tab == tab); + Circuit res_uni = + cm_tableau_to_unitary_extension_circuit(tab, circ.all_qubits()).first; + REQUIRE(test_statevector_comparison(res, res_uni, true)); + } + GIVEN("A partial Clifford state (tests mixed initialisations)") { + Circuit circ(3); + add_ops_list_one_to_circuit(circ); + circ.add_op(OpType::Collapse, {1}); + circ.qubit_create_all(); + ChoiMixTableau tab = circuit_to_cm_tableau(circ); + Circuit res = cm_tableau_to_exact_circuit(tab).first; + CHECK(res.created_qubits().size() == 3); + CHECK(res.discarded_qubits().size() == 0); + CHECK(res.count_gates(OpType::Collapse) == 1); ChoiMixTableau res_tab = circuit_to_cm_tableau(res); - tab.canonical_column_order(); - tab.gaussian_form(); - res_tab.canonical_column_order(); - res_tab.gaussian_form(); REQUIRE(res_tab == tab); + Circuit res_uni = + cm_tableau_to_unitary_extension_circuit(tab, circ.all_qubits()).first; + Eigen::VectorXcd res_sv = tket_sim::get_statevector(res_uni); + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + Eigen::MatrixXcd outmat = rrow.second.to_sparse_matrix(3); + CHECK((outmat * res_sv).isApprox(res_sv)); + } } GIVEN("A total diagonalisation circuit") { Circuit circ = get_test_circ(); @@ -523,9 +585,14 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { circ.add_op(OpType::Collapse, {i}); } ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; + Circuit res = cm_tableau_to_exact_circuit(tab).first; ChoiMixTableau res_tab = circuit_to_cm_tableau(res); REQUIRE(res_tab == tab); + // Test unitary synthesis by statevector of dagger + Circuit as_state = get_test_circ().dagger(); + Circuit res_uni_dag = + cm_tableau_to_unitary_extension_circuit(tab).first.dagger(); + REQUIRE(test_statevector_comparison(as_state, res_uni_dag, true)); } GIVEN("A partial diagonalisation circuit") { Circuit circ = get_test_circ(); @@ -534,18 +601,26 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { } circ.qubit_discard(Qubit(0)); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - std::pair res = cm_tableau_to_circuit(tab); + std::pair res = cm_tableau_to_exact_circuit(tab); ChoiMixTableau res_tab = circuit_to_cm_tableau(res.first); qubit_map_t perm; - for (const std::pair& p : res.second) { - perm.insert({Qubit(p.second), Qubit(p.first)}); + for (const std::pair& p : res.second) { + perm.insert({p.second, p.first}); } res_tab.rename_qubits(perm, ChoiMixTableau::TableauSegment::Output); - tab.canonical_column_order(); - tab.gaussian_form(); res_tab.canonical_column_order(); res_tab.gaussian_form(); REQUIRE(res_tab == tab); + Circuit res_uni_dag = + cm_tableau_to_unitary_extension_circuit(tab).first.dagger(); + Eigen::VectorXcd as_state = tket_sim::get_statevector(res_uni_dag); + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + CmplxSpMat rmat = rrow.first.to_sparse_matrix(3); + if (rrow.second.is_real_negative()) rmat *= -1.; + Eigen::MatrixXcd rmatd = rmat; + CHECK((rmat * as_state).isApprox(as_state)); + } } GIVEN("Another circuit for extra test coverage in row reductions") { Circuit circ(5); @@ -564,13 +639,24 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { circ.add_op(OpType::Collapse, {4}); circ.add_op(OpType::H, {4}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - Circuit res = cm_tableau_to_circuit(tab).first; - ChoiMixTableau res_tab = circuit_to_cm_tableau(res); - tab.canonical_column_order(); - tab.gaussian_form(); - res_tab.canonical_column_order(); - res_tab.gaussian_form(); + std::pair res = cm_tableau_to_exact_circuit(tab); + res.first.permute_boundary_output(inv_perm(res.second)); + ChoiMixTableau res_tab = circuit_to_cm_tableau(res.first); REQUIRE(res_tab == tab); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit(tab); + res_uni.first.permute_boundary_output(inv_perm(res_uni.second)); + res_tab = circuit_to_cm_tableau(res_uni.first); + res_tab.tab_.row_mult(0, 1); + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + CmplxSpMat inmat = rrow.first.to_sparse_matrix(5); + Eigen::MatrixXcd inmatd = inmat; + CmplxSpMat outmat = rrow.second.to_sparse_matrix(5); + Eigen::MatrixXcd outmatd = outmat; + CHECK((outmatd * res_u * inmatd).isApprox(res_u)); + } } GIVEN("An isometry") { Circuit circ(5); @@ -586,18 +672,36 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { circ.add_op(OpType::CX, {1, 2}); circ.add_op(OpType::CX, {1, 0}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - std::pair res = cm_tableau_to_circuit(tab); + std::pair res = cm_tableau_to_exact_circuit(tab); ChoiMixTableau res_tab = circuit_to_cm_tableau(res.first); qubit_map_t perm; - for (const std::pair& p : res.second) { - perm.insert({Qubit(p.second), Qubit(p.first)}); + for (const std::pair& p : res.second) { + perm.insert({p.second, p.first}); } res_tab.rename_qubits(perm, ChoiMixTableau::TableauSegment::Output); - tab.canonical_column_order(); - tab.gaussian_form(); res_tab.canonical_column_order(); res_tab.gaussian_form(); REQUIRE(res_tab == tab); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit( + tab, {Qubit(1), Qubit(2), Qubit(3)}); + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + Eigen::MatrixXcd init_proj = Eigen::MatrixXcd::Zero(32, 32); + init_proj.block(0, 0, 2, 2) = Eigen::MatrixXcd::Identity(2, 2); + init_proj.block(16, 16, 2, 2) = Eigen::MatrixXcd::Identity(2, 2); + Eigen::MatrixXcd res_iso = res_u * init_proj; + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + CmplxSpMat inmat = rrow.first.to_sparse_matrix(5); + Eigen::MatrixXcd inmatd = inmat; + QubitPauliMap outstr; + for (const std::pair& qp : rrow.second.string) + outstr.insert({res_uni.second.at(qp.first), qp.second}); + CmplxSpMat outmat = SpPauliString(outstr).to_sparse_matrix(5); + Eigen::MatrixXcd outmatd = outmat; + if (rrow.second.is_real_negative()) outmatd *= -1.; + CHECK((outmatd * res_iso * inmatd).isApprox(res_iso)); + } } GIVEN("Extra coverage for isometries") { Circuit circ(5); @@ -615,18 +719,211 @@ SCENARIO("Synthesis of circuits from ChoiMixTableaus") { circ.add_op(OpType::CX, {1, 2}); circ.add_op(OpType::CX, {1, 0}); ChoiMixTableau tab = circuit_to_cm_tableau(circ); - std::pair res = cm_tableau_to_circuit(tab); + std::pair res = cm_tableau_to_exact_circuit(tab); ChoiMixTableau res_tab = circuit_to_cm_tableau(res.first); qubit_map_t perm; - for (const std::pair& p : res.second) { - perm.insert({Qubit(p.second), Qubit(p.first)}); + for (const std::pair& p : res.second) { + perm.insert({p.second, p.first}); } res_tab.rename_qubits(perm, ChoiMixTableau::TableauSegment::Output); + res_tab.canonical_column_order(); + res_tab.gaussian_form(); + REQUIRE(res_tab == tab); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit( + tab, {Qubit(1), Qubit(2), Qubit(3)}); + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + Eigen::MatrixXcd init_proj = Eigen::MatrixXcd::Zero(32, 32); + init_proj.block(0, 0, 2, 2) = Eigen::MatrixXcd::Identity(2, 2); + init_proj.block(16, 16, 2, 2) = Eigen::MatrixXcd::Identity(2, 2); + Eigen::MatrixXcd res_iso = res_u * init_proj; + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + CmplxSpMat inmat = rrow.first.to_sparse_matrix(5); + Eigen::MatrixXcd inmatd = inmat; + QubitPauliMap outstr; + for (const std::pair& qp : rrow.second.string) + outstr.insert({res_uni.second.at(qp.first), qp.second}); + CmplxSpMat outmat = SpPauliString(outstr).to_sparse_matrix(5); + Eigen::MatrixXcd outmatd = outmat; + if (rrow.second.is_real_negative()) outmatd *= -1.; + CHECK((outmatd * res_iso * inmatd).isApprox(res_iso)); + } + } + GIVEN("Synthesising a tableau requiring post-selection") { + Circuit circ = get_test_circ(); + ChoiMixTableau tab = circuit_to_cm_tableau(circ); + tab.post_select(Qubit(0), ChoiMixTableau::TableauSegment::Output); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit(tab, {}, {Qubit(0)}); + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + // q[0] was removed from the tableau by postselection so need to infer + // position in res_uni.second from the other qubits + SpPauliString zzz({Pauli::Z, Pauli::Z, Pauli::Z}); + zzz.set(res_uni.second.at(Qubit(1)), Pauli::I); + zzz.set(res_uni.second.at(Qubit(2)), Pauli::I); + Eigen::MatrixXcd z0 = zzz.to_sparse_matrix(3); + Eigen::MatrixXcd res_proj = 0.5 * (res_u + (z0 * res_u)); + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + CmplxSpMat inmat = rrow.first.to_sparse_matrix(3); + Eigen::MatrixXcd inmatd = inmat; + QubitPauliMap outstr; + for (const std::pair& qp : rrow.second.string) + outstr.insert({res_uni.second.at(qp.first), qp.second}); + CmplxSpMat outmat = SpPauliString(outstr).to_sparse_matrix(3); + Eigen::MatrixXcd outmatd = outmat; + if (rrow.second.is_real_negative()) outmatd *= -1.; + CHECK((outmatd * res_proj * inmatd).isApprox(res_proj)); + } + } + GIVEN("Synthesising a tableau with all post-selections") { + Circuit circ = get_test_circ(); + ChoiMixTableau tab = circuit_to_cm_tableau(circ); + tab.post_select(Qubit(0), ChoiMixTableau::TableauSegment::Output); + tab.post_select(Qubit(1), ChoiMixTableau::TableauSegment::Output); + tab.post_select(Qubit(2), ChoiMixTableau::TableauSegment::Output); + Circuit res = cm_tableau_to_unitary_extension_circuit( + tab, {}, {Qubit(0), Qubit(1), Qubit(2)}) + .first.dagger(); + Eigen::VectorXcd res_sv = tket_sim::get_statevector(res); + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + Eigen::MatrixXcd inmat = rrow.first.to_sparse_matrix(3); + if (rrow.second.is_real_negative()) inmat *= -1.; + CHECK((inmat * res_sv).isApprox(res_sv)); + } + } + GIVEN("Initialisations, collapses, discards and post-selections") { + Circuit circ(5); + circ.qubit_create(Qubit(1)); + circ.qubit_create(Qubit(2)); + circ.add_op(OpType::H, {4}); + circ.add_op(OpType::Collapse, {4}); + circ.add_op(OpType::CX, {4, 1}); + circ.add_op(OpType::CX, {4, 2}); + circ.add_op(OpType::CX, {4, 3}); + circ.add_op(OpType::H, {4}); + circ.add_op(OpType::H, {1}); + circ.add_op(OpType::V, {2}); + circ.add_op(OpType::CX, {1, 2}); + circ.add_op(OpType::CX, {1, 0}); + circ.qubit_discard(Qubit(0)); + ChoiMixTableau tab = circuit_to_cm_tableau(circ); + tab.post_select(Qubit(3), ChoiMixTableau::TableauSegment::Output); tab.canonical_column_order(); tab.gaussian_form(); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit(tab, {Qubit(1)}, {Qubit(0)}); + // First rebuild tableau by initialising, post-selecting, etc. + ChoiMixTableau res_tab = circuit_to_cm_tableau(res_uni.first); + qubit_map_t perm; + for (const std::pair& p : res_uni.second) + perm.insert({p.second, p.first}); + res_tab.rename_qubits(perm, ChoiMixTableau::TableauSegment::Output); + // Post-select/initialise + res_tab.post_select(Qubit(1), ChoiMixTableau::TableauSegment::Input); + res_tab.post_select(Qubit(0), ChoiMixTableau::TableauSegment::Output); + // Collapsing q[4] in X basis as per circ + res_tab.apply_gate( + OpType::H, {Qubit(4)}, ChoiMixTableau::TableauSegment::Output); + res_tab.collapse_qubit(Qubit(4), ChoiMixTableau::TableauSegment::Output); + res_tab.apply_gate( + OpType::H, {Qubit(4)}, ChoiMixTableau::TableauSegment::Output); + // Discarding q[0] also removes Z row for q[0], so recreate this by + // XCollapse at input + res_tab.apply_gate( + OpType::H, {Qubit(0)}, ChoiMixTableau::TableauSegment::Input); + res_tab.collapse_qubit(Qubit(0), ChoiMixTableau::TableauSegment::Input); + res_tab.apply_gate( + OpType::H, {Qubit(0)}, ChoiMixTableau::TableauSegment::Input); res_tab.canonical_column_order(); res_tab.gaussian_form(); REQUIRE(res_tab == tab); + + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + qubit_vector_t res_qbs = res_uni.first.all_qubits(); + // q[1] has no input terms, so initialise it + SpPauliString z1({Qubit(1), Pauli::Z}); + Eigen::MatrixXcd z1u = z1.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (res_u * z1u)); + // q[0] has no output terms, so postselect it + SpPauliString z0({res_uni.second.at(Qubit(0)), Pauli::Z}); + Eigen::MatrixXcd z0u = z0.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (z0u * res_u)); + + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + Eigen::MatrixXcd inmat = rrow.first.to_sparse_matrix(res_qbs); + QubitPauliMap outstr; + for (const std::pair& qp : rrow.second.string) + outstr.insert({res_uni.second.at(qp.first), qp.second}); + Eigen::MatrixXcd outmat = SpPauliString(outstr).to_sparse_matrix(res_qbs); + if (rrow.second.is_real_negative()) outmat *= -1.; + CHECK((outmat * res_u * inmat).isApprox(res_u)); + } + } + GIVEN( + "A custom tableau with overlapping initialised and post-selected " + "qubits") { + std::list rows{ + {SpPauliStabiliser({Pauli::Z, Pauli::X, Pauli::I}), {}}, + {SpPauliStabiliser({Pauli::X, Pauli::Y, Pauli::Z}), {}}, + {{}, SpPauliStabiliser({Pauli::X, Pauli::X, Pauli::I})}, + {{}, SpPauliStabiliser({Pauli::I, Pauli::X, Pauli::X})}, + {SpPauliStabiliser({Pauli::I, Pauli::I, Pauli::Z}), + SpPauliStabiliser({Pauli::Z, Pauli::Z, Pauli::Z})}, + {SpPauliStabiliser({Pauli::Z, Pauli::I, Pauli::X}), + SpPauliStabiliser({Pauli::I, Pauli::I, Pauli::X})}, + }; + ChoiMixTableau tab(rows); + REQUIRE_THROWS(cm_tableau_to_unitary_extension_circuit(tab)); + std::pair res_uni = + cm_tableau_to_unitary_extension_circuit( + tab, {Qubit(3), Qubit(4)}, {Qubit(3), Qubit(4)}); + + ChoiMixTableau res_tab = circuit_to_cm_tableau(res_uni.first); + qubit_map_t perm; + for (const std::pair& p : res_uni.second) + perm.insert({p.second, p.first}); + res_tab.rename_qubits(perm, ChoiMixTableau::TableauSegment::Output); + res_tab.post_select(Qubit(3), ChoiMixTableau::TableauSegment::Input); + res_tab.post_select(Qubit(4), ChoiMixTableau::TableauSegment::Input); + res_tab.post_select(Qubit(3), ChoiMixTableau::TableauSegment::Output); + res_tab.post_select(Qubit(4), ChoiMixTableau::TableauSegment::Output); + res_tab.canonical_column_order(); + res_tab.gaussian_form(); + tab.canonical_column_order(); + tab.gaussian_form(); + REQUIRE(res_tab == tab); + + Eigen::MatrixXcd res_u = tket_sim::get_unitary(res_uni.first); + qubit_vector_t res_qbs = res_uni.first.all_qubits(); + // initialise q[3] and q[4] + SpPauliString z3i{Qubit(3), Pauli::Z}; + Eigen::MatrixXcd z3iu = z3i.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (res_u * z3iu)); + SpPauliString z4i{Qubit(4), Pauli::Z}; + Eigen::MatrixXcd z4iu = z4i.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (res_u * z4iu)); + // post-select q[3] and q[4] + SpPauliString z3o{res_uni.second.at(Qubit(3)), Pauli::Z}; + Eigen::MatrixXcd z3ou = z3o.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (z3ou * res_u)); + SpPauliString z4o{res_uni.second.at(Qubit(4)), Pauli::Z}; + Eigen::MatrixXcd z4ou = z4o.to_sparse_matrix(res_qbs); + res_u = 0.5 * (res_u + (z4ou * res_u)); + + for (unsigned r = 0; r < tab.get_n_rows(); ++r) { + ChoiMixTableau::row_tensor_t rrow = tab.get_row(r); + Eigen::MatrixXcd inmat = rrow.first.to_sparse_matrix(res_qbs); + QubitPauliMap outstr; + for (const std::pair& qp : rrow.second.string) + outstr.insert({res_uni.second.at(qp.first), qp.second}); + Eigen::MatrixXcd outmat = SpPauliString(outstr).to_sparse_matrix(res_qbs); + if (rrow.second.is_real_negative()) outmat *= -1.; + CHECK((outmat * res_u * inmat).isApprox(res_u)); + } } } diff --git a/tket/test/src/test_Diagonalisation.cpp b/tket/test/src/test_Diagonalisation.cpp new file mode 100644 index 0000000000..a2778dc07a --- /dev/null +++ b/tket/test/src/test_Diagonalisation.cpp @@ -0,0 +1,148 @@ +// Copyright 2019-2024 Cambridge Quantum 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. + +#include + +#include "tket/Circuit/Simulation/CircuitSimulator.hpp" +#include "tket/Diagonalisation/Diagonalisation.hpp" + +namespace tket { + +namespace test_Diagonalisation { + +SCENARIO("Matrix tests for reducing a Pauli to Z") { + std::list test_paulis{Pauli::I, Pauli::X, Pauli::Y, Pauli::Z}; + std::list test_configs = { + CXConfigType::Snake, CXConfigType::Tree, CXConfigType::Star, + CXConfigType::MultiQGate}; + for (const Pauli& p : test_paulis) { + // If p is Pauli::I, it is dropped from the sparse representation in the + // constructor, so need to add Qubit(3) to the circuit and sparse matrix + // explicitly + SpPauliStabiliser pt({Pauli::X, Pauli::Y, Pauli::Z, p}); + CmplxSpMat pt_u = pt.to_sparse_matrix(4); + Eigen::MatrixXcd pt_ud = pt_u; + for (const CXConfigType& config : test_configs) { + std::pair diag = reduce_pauli_to_z(pt, config); + if (p == Pauli::I) diag.first.add_qubit(Qubit(3)); + Eigen::MatrixXcd diag_u = tket_sim::get_unitary(diag.first); + CmplxSpMat z_u = SpPauliString(diag.second, Pauli::Z).to_sparse_matrix(4); + Eigen::MatrixXcd z_ud = z_u; + CHECK((z_ud * diag_u * pt_ud).isApprox(diag_u)); + } + } +} + +SCENARIO("Matrix tests for reducing two anticommuting Paulis to Z X") { + std::list non_trivials{Pauli::X, Pauli::Y, Pauli::Z}; + std::list test_configs = { + CXConfigType::Snake, CXConfigType::Tree, CXConfigType::Star, + CXConfigType::MultiQGate}; + // Loop through all commuting options for two qubits + for (const Pauli& p0 : non_trivials) { + for (const Pauli& p1 : non_trivials) { + SpPauliStabiliser p({Pauli::Z, p0, p1, Pauli::Z}); + CmplxSpMat p_u = p.to_sparse_matrix(); + Eigen::MatrixXcd p_ud = p_u; + for (const Pauli& q0 : non_trivials) { + for (const Pauli& q1 : non_trivials) { + SpPauliStabiliser q({Pauli::X, q0, q1, Pauli::Z}); + if (p.commutes_with(q)) continue; + CmplxSpMat q_u = q.to_sparse_matrix(); + Eigen::MatrixXcd q_ud = q_u; + for (const CXConfigType& config : test_configs) { + std::pair diag = + reduce_anticommuting_paulis_to_z_x(p, q, config); + Eigen::MatrixXcd diag_u = tket_sim::get_unitary(diag.first); + CmplxSpMat z_u = + SpPauliString(diag.second, Pauli::Z).to_sparse_matrix(4); + Eigen::MatrixXcd z_ud = z_u; + CmplxSpMat x_u = + SpPauliString(diag.second, Pauli::X).to_sparse_matrix(4); + Eigen::MatrixXcd x_ud = x_u; + CHECK((z_ud * diag_u * p_ud).isApprox(diag_u)); + CHECK((x_ud * diag_u * q_ud).isApprox(diag_u)); + } + } + } + } + } +} + +SCENARIO("Matrix tests for reducing two commuting Paulis to Z X") { + std::list paulis{Pauli::I, Pauli::X, Pauli::Y, Pauli::Z}; + std::list test_configs = { + CXConfigType::Snake, CXConfigType::Tree, CXConfigType::Star, + CXConfigType::MultiQGate}; + // Loop through all commuting options for two qubits + for (const Pauli& p0 : paulis) { + for (const Pauli& p1 : paulis) { + SpPauliStabiliser p({Pauli::Z, p0, p1, Pauli::Z}); + CmplxSpMat p_u = p.to_sparse_matrix(4); + Eigen::MatrixXcd p_ud = p_u; + for (const Pauli& q0 : paulis) { + for (const Pauli& q1 : paulis) { + SpPauliStabiliser q({Pauli::Z, q0, q1, Pauli::I}); + if (!p.commutes_with(q)) continue; + CmplxSpMat q_u = q.to_sparse_matrix(4); + Eigen::MatrixXcd q_ud = q_u; + for (const CXConfigType& config : test_configs) { + std::tuple diag = + reduce_commuting_paulis_to_zi_iz(p, q, config); + Circuit& circ = std::get<0>(diag); + // In cases with matching Pauli::Is, the circuit produced may not + // include all qubits + for (unsigned i = 0; i < 4; ++i) circ.add_qubit(Qubit(i), false); + Eigen::MatrixXcd diag_u = tket_sim::get_unitary(circ); + CmplxSpMat zi_u = + SpPauliString(std::get<1>(diag), Pauli::Z).to_sparse_matrix(4); + Eigen::MatrixXcd zi_ud = zi_u; + CmplxSpMat iz_u = + SpPauliString(std::get<2>(diag), Pauli::Z).to_sparse_matrix(4); + Eigen::MatrixXcd iz_ud = iz_u; + CHECK((zi_ud * diag_u * p_ud).isApprox(diag_u)); + CHECK((iz_ud * diag_u * q_ud).isApprox(diag_u)); + } + } + } + } + } +} + +SCENARIO("Reducing shared qubits to no matches") { + GIVEN("Strings with no mismatch and second completely contains the first") { + SpPauliStabiliser pauli0({{Qubit(0), Pauli::X}, {Qubit(1), Pauli::Y}}); + SpPauliStabiliser pauli1( + {{Qubit(0), Pauli::X}, {Qubit(1), Pauli::Y}, {Qubit(2), Pauli::Z}}); + SpPauliStabiliser p0_orig = pauli0; + SpPauliStabiliser p1_orig = pauli1; + std::pair> sol = + reduce_overlap_of_paulis(pauli0, pauli1, CXConfigType::Snake, false); + // Check there is no final overlapping qubit returned + CHECK(!sol.second); + CHECK(pauli0.common_qubits(pauli1).empty()); + CHECK(pauli0.conflicting_qubits(pauli1).empty()); + // Check the strings are updated correctly to match the unitary + Eigen::MatrixXcd diag_u = tket_sim::get_unitary(sol.first); + Eigen::MatrixXcd p0_u = pauli0.to_sparse_matrix(3); + Eigen::MatrixXcd p0_o_u = p0_orig.to_sparse_matrix(3); + CHECK((p0_u * diag_u * p0_o_u).isApprox(diag_u)); + Eigen::MatrixXcd p1_u = pauli1.to_sparse_matrix(3); + Eigen::MatrixXcd p1_o_u = p1_orig.to_sparse_matrix(3); + CHECK((p1_u * diag_u * p1_o_u).isApprox(diag_u)); + } +} + +} // namespace test_Diagonalisation +} // namespace tket \ No newline at end of file diff --git a/tket/test/src/test_PauliGraph.cpp b/tket/test/src/test_PauliGraph.cpp index 2895899f09..7ddd257586 100644 --- a/tket/test/src/test_PauliGraph.cpp +++ b/tket/test/src/test_PauliGraph.cpp @@ -18,10 +18,10 @@ #include "CircuitsForTesting.hpp" #include "testutil.hpp" #include "tket/Circuit/Boxes.hpp" +#include "tket/Circuit/CircUtils.hpp" #include "tket/Circuit/PauliExpBoxes.hpp" #include "tket/Circuit/Simulation/CircuitSimulator.hpp" #include "tket/Converters/Converters.hpp" -#include "tket/Converters/PauliGadget.hpp" #include "tket/Diagonalisation/Diagonalisation.hpp" #include "tket/Gate/SymTable.hpp" #include "tket/PauliGraph/ConjugatePauliFunctions.hpp" @@ -920,13 +920,13 @@ SCENARIO("Diagonalise a pair of gadgets") { Circuit correct; for (unsigned i = 0; i < 2; ++i) { - append_single_pauli_gadget(correct, gadgets.at(i)); + correct.append(pauli_gadget(gadgets.at(i))); } auto u_correct = tket_sim::get_unitary(correct); GIVEN("Snake configuration") { CXConfigType config = CXConfigType::Snake; - append_pauli_gadget_pair(circ, gadgets.at(0), gadgets.at(1), config); + circ.append(pauli_gadget_pair(gadgets.at(0), gadgets.at(1), config)); THEN("Unitary is correct") { auto u_res = tket_sim::get_unitary(circ); REQUIRE((u_correct - u_res).cwiseAbs().sum() < ERR_EPS); @@ -934,7 +934,7 @@ SCENARIO("Diagonalise a pair of gadgets") { } GIVEN("Star configuration") { CXConfigType config = CXConfigType::Star; - append_pauli_gadget_pair(circ, gadgets.at(0), gadgets.at(1), config); + circ.append(pauli_gadget_pair(gadgets.at(0), gadgets.at(1), config)); THEN("Unitary is correct") { auto u_res = tket_sim::get_unitary(circ); REQUIRE((u_correct - u_res).cwiseAbs().sum() < ERR_EPS); @@ -942,7 +942,7 @@ SCENARIO("Diagonalise a pair of gadgets") { } GIVEN("Tree configuration") { CXConfigType config = CXConfigType::Tree; - append_pauli_gadget_pair(circ, gadgets.at(0), gadgets.at(1), config); + circ.append(pauli_gadget_pair(gadgets.at(0), gadgets.at(1), config)); THEN("Unitary is correct") { auto u_res = tket_sim::get_unitary(circ); REQUIRE((u_correct - u_res).cwiseAbs().sum() < ERR_EPS); @@ -950,7 +950,7 @@ SCENARIO("Diagonalise a pair of gadgets") { } GIVEN("MultiQGate configuration") { CXConfigType config = CXConfigType::MultiQGate; - append_pauli_gadget_pair(circ, gadgets.at(0), gadgets.at(1), config); + circ.append(pauli_gadget_pair(gadgets.at(0), gadgets.at(1), config)); circ.decompose_boxes_recursively(); THEN("XXPhase3 were used") { REQUIRE(circ.count_gates(OpType::XXPhase3) == 2); diff --git a/tket/test/src/test_PhaseGadget.cpp b/tket/test/src/test_PhaseGadget.cpp index 9be60aa248..56f5182843 100644 --- a/tket/test/src/test_PhaseGadget.cpp +++ b/tket/test/src/test_PhaseGadget.cpp @@ -149,7 +149,8 @@ SCENARIO("Constructing Pauli gadgets") { 1.*i_, 0., 0., 0.; // clang-format on Eigen::Matrix4cd m = (+0.5 * PI * i_ * t * a).exp(); - Circuit circ = pauli_gadget({Pauli::X, Pauli::Y}, t); + Circuit circ = + pauli_gadget(SpSymPauliTensor(DensePauliMap{Pauli::X, Pauli::Y}, t)); const Eigen::Matrix4cd u = tket_sim::get_unitary(circ); Eigen::Matrix4cd v = m * u; REQUIRE((v - Eigen::Matrix4cd::Identity()).cwiseAbs().sum() < ERR_EPS); @@ -286,7 +287,7 @@ SCENARIO("Decompose phase gadgets") { for (unsigned i = 0; i < n_qubits; ++i) { nZ.push_back(Pauli::Z); } - Circuit pauli_gadget_circ = pauli_gadget(nZ, 0.2); + Circuit pauli_gadget_circ = pauli_gadget(SpSymPauliTensor(nZ, 0.2)); pauli_gadget_circ.decompose_boxes_recursively(); Circuit phase_gadget_circ = phase_gadget(n_qubits, 0.2, config); phase_gadget_circ.decompose_boxes_recursively(); diff --git a/tket/test/src/test_UnitaryTableau.cpp b/tket/test/src/test_UnitaryTableau.cpp index 67b31fab98..6693ff658d 100644 --- a/tket/test/src/test_UnitaryTableau.cpp +++ b/tket/test/src/test_UnitaryTableau.cpp @@ -16,6 +16,7 @@ #include #include "testutil.hpp" +#include "tket/Circuit/Simulation/CircuitSimulator.hpp" #include "tket/Clifford/UnitaryTableau.hpp" #include "tket/Converters/Converters.hpp" #include "tket/Converters/UnitaryTableauBox.hpp" @@ -246,11 +247,59 @@ SCENARIO("Correct creation of UnitaryTableau") { CHECK(tab0 == tab4); CHECK(tab0 == tab5); } + GIVEN("A single Z gate") { + UnitaryTableau tab0(3); + UnitaryTableau tab1(3); + UnitaryTableau tab2(3); + UnitaryTableau tab3(3); + UnitaryTableau tab4(3); + UnitaryTableau tab5(3); + tab0.apply_gate_at_front(OpType::Z, {Qubit(0)}); + tab1.apply_gate_at_end(OpType::Z, {Qubit(0)}); + tab2.apply_gate_at_front(OpType::S, {Qubit(0)}); + tab2.apply_gate_at_front(OpType::S, {Qubit(0)}); + tab3.apply_gate_at_end(OpType::Sdg, {Qubit(0)}); + tab3.apply_gate_at_end(OpType::Sdg, {Qubit(0)}); + tab4.apply_Z_at_front(Qubit(0)); + tab5.apply_Z_at_end(Qubit(0)); + CHECK(tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z)); + CHECK(tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X, 2)); + CHECK(tab0 == tab1); + CHECK(tab0 == tab2); + CHECK(tab0 == tab3); + CHECK(tab0 == tab4); + CHECK(tab0 == tab5); + } + GIVEN("A single X gate") { + UnitaryTableau tab0(3); + UnitaryTableau tab1(3); + UnitaryTableau tab2(3); + UnitaryTableau tab3(3); + UnitaryTableau tab4(3); + UnitaryTableau tab5(3); + tab0.apply_gate_at_front(OpType::X, {Qubit(0)}); + tab1.apply_gate_at_end(OpType::X, {Qubit(0)}); + tab2.apply_gate_at_front(OpType::V, {Qubit(0)}); + tab2.apply_gate_at_front(OpType::V, {Qubit(0)}); + tab3.apply_gate_at_end(OpType::Vdg, {Qubit(0)}); + tab3.apply_gate_at_end(OpType::Vdg, {Qubit(0)}); + tab4.apply_X_at_front(Qubit(0)); + tab5.apply_X_at_end(Qubit(0)); + CHECK(tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z, 2)); + CHECK(tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X)); + CHECK(tab0 == tab1); + CHECK(tab0 == tab2); + CHECK(tab0 == tab3); + CHECK(tab0 == tab4); + CHECK(tab0 == tab5); + } GIVEN("A single H gate") { UnitaryTableau tab0(3); UnitaryTableau tab1(3); UnitaryTableau tab2(3); UnitaryTableau tab3(3); + UnitaryTableau tab4(3); + UnitaryTableau tab5(3); tab0.apply_gate_at_front(OpType::H, {Qubit(0)}); tab1.apply_gate_at_end(OpType::H, {Qubit(0)}); tab2.apply_gate_at_front(OpType::S, {Qubit(0)}); @@ -259,11 +308,15 @@ SCENARIO("Correct creation of UnitaryTableau") { tab3.apply_gate_at_end(OpType::Vdg, {Qubit(0)}); tab3.apply_gate_at_end(OpType::Sdg, {Qubit(0)}); tab3.apply_gate_at_end(OpType::Vdg, {Qubit(0)}); + tab4.apply_H_at_front(Qubit(0)); + tab5.apply_H_at_end(Qubit(0)); CHECK(tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X)); CHECK(tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z)); CHECK(tab0 == tab1); CHECK(tab0 == tab2); CHECK(tab0 == tab3); + CHECK(tab0 == tab4); + CHECK(tab0 == tab5); } GIVEN("A single CX gate") { UnitaryTableau tab0(3); @@ -297,6 +350,21 @@ SCENARIO("Correct creation of UnitaryTableau") { UnitaryTableau tab = circuit_to_unitary_tableau(circ); UnitaryTableau rev_tab = get_tableau_with_gates_applied_at_front(); REQUIRE(tab == rev_tab); + Eigen::MatrixXcd circ_u = tket_sim::get_unitary(circ); + for (unsigned q = 0; q < 3; ++q) { + CmplxSpMat xq = SpPauliString(Qubit(q), Pauli::X).to_sparse_matrix(3); + Eigen::MatrixXcd xqd = xq; + PauliStabiliser xrow = tab.get_xrow(Qubit(q)); + CmplxSpMat xrowmat = xrow.to_sparse_matrix(3); + Eigen::MatrixXcd xrowmatd = xrowmat; + CHECK((xrowmatd * circ_u * xqd).isApprox(circ_u)); + CmplxSpMat zq = SpPauliString(Qubit(q), Pauli::Z).to_sparse_matrix(3); + Eigen::MatrixXcd zqd = zq; + PauliStabiliser zrow = tab.get_zrow(Qubit(q)); + CmplxSpMat zrowmat = zrow.to_sparse_matrix(3); + Eigen::MatrixXcd zrowmatd = zrowmat; + CHECK((zrowmatd * circ_u * zqd).isApprox(circ_u)); + } } GIVEN("A PI/2 rotation") { Circuit circ = get_test_circ(); @@ -359,6 +427,22 @@ SCENARIO("Synthesis of circuits from UnitaryTableau") { UnitaryTableau res_tab = circuit_to_unitary_tableau(res); REQUIRE(res_tab == tab); } + GIVEN("Additional gate coverage, check unitary") { + Circuit circ(4); + circ.add_op(OpType::ZZMax, {0, 1}); + circ.add_op(OpType::ECR, {2, 3}); + circ.add_op(OpType::ISWAPMax, {1, 2}); + circ.add_op(OpType::noop, {0}); + UnitaryTableau tab = circuit_to_unitary_tableau(circ); + UnitaryTableau rev_tab(4); + rev_tab.apply_gate_at_front(OpType::noop, {Qubit(0)}); + rev_tab.apply_gate_at_front(OpType::ISWAPMax, {Qubit(1), Qubit(2)}); + rev_tab.apply_gate_at_front(OpType::ECR, {Qubit(2), Qubit(3)}); + rev_tab.apply_gate_at_front(OpType::ZZMax, {Qubit(0), Qubit(1)}); + REQUIRE(tab == rev_tab); + Circuit res = unitary_tableau_to_circuit(tab); + REQUIRE(test_unitary_comparison(circ, res, true)); + } } SCENARIO("Correct creation of UnitaryRevTableau") { @@ -418,14 +502,52 @@ SCENARIO("Correct creation of UnitaryRevTableau") { CHECK(tab0 == tab4); CHECK(tab0 == tab5); } + GIVEN("A single Z gate") { + UnitaryRevTableau tab0(3); + UnitaryRevTableau tab1(3); + UnitaryRevTableau tab2(3); + UnitaryRevTableau tab3(3); + tab0.apply_gate_at_end(OpType::Z, {Qubit(0)}); + tab1.apply_gate_at_front(OpType::Z, {Qubit(0)}); + tab2.apply_Z_at_end(Qubit(0)); + tab3.apply_Z_at_front(Qubit(0)); + REQUIRE(tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z)); + REQUIRE( + tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X, 2)); + REQUIRE(tab0 == tab1); + REQUIRE(tab0 == tab2); + REQUIRE(tab0 == tab3); + } + GIVEN("A single X gate") { + UnitaryRevTableau tab0(3); + UnitaryRevTableau tab1(3); + UnitaryRevTableau tab2(3); + UnitaryRevTableau tab3(3); + tab0.apply_gate_at_end(OpType::X, {Qubit(0)}); + tab1.apply_gate_at_front(OpType::X, {Qubit(0)}); + tab2.apply_X_at_end(Qubit(0)); + tab3.apply_X_at_front(Qubit(0)); + REQUIRE( + tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z, 2)); + REQUIRE(tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X)); + REQUIRE(tab0 == tab1); + REQUIRE(tab0 == tab2); + REQUIRE(tab0 == tab3); + } GIVEN("A single H gate") { UnitaryRevTableau tab0(3); UnitaryRevTableau tab1(3); + UnitaryRevTableau tab2(3); + UnitaryRevTableau tab3(3); tab0.apply_gate_at_end(OpType::H, {Qubit(0)}); tab1.apply_gate_at_front(OpType::H, {Qubit(0)}); + tab2.apply_H_at_end(Qubit(0)); + tab3.apply_H_at_front(Qubit(0)); REQUIRE(tab0.get_zrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::X)); REQUIRE(tab0.get_xrow(Qubit(0)) == SpPauliStabiliser(Qubit(0), Pauli::Z)); REQUIRE(tab0 == tab1); + REQUIRE(tab0 == tab2); + REQUIRE(tab0 == tab3); } GIVEN("A single CX gate") { UnitaryRevTableau tab0(3); @@ -538,6 +660,29 @@ SCENARIO("Synthesis of circuits from UnitaryRevTableau") { REQUIRE(test_unitary_comparison(circ, res, true)); UnitaryRevTableau res_tab = circuit_to_unitary_rev_tableau(res); REQUIRE(res_tab == tab); + // Manually call apply_gate methods + UnitaryRevTableau tab_front(2); + tab_front.apply_gate_at_front(OpType::ISWAPMax, {Qubit(0), Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SXdg, {Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SX, {Qubit(0)}); + tab_front.apply_gate_at_front(OpType::ECR, {Qubit(0), Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SXdg, {Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SX, {Qubit(0)}); + tab_front.apply_gate_at_front(OpType::ZZMax, {Qubit(0), Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SXdg, {Qubit(1)}); + tab_front.apply_gate_at_front(OpType::SX, {Qubit(0)}); + REQUIRE(tab == tab_front); + UnitaryRevTableau tab_end(2); + tab_end.apply_gate_at_end(OpType::SX, {Qubit(0)}); + tab_end.apply_gate_at_end(OpType::SXdg, {Qubit(1)}); + tab_end.apply_gate_at_end(OpType::ZZMax, {Qubit(0), Qubit(1)}); + tab_end.apply_gate_at_end(OpType::SX, {Qubit(0)}); + tab_end.apply_gate_at_end(OpType::SXdg, {Qubit(1)}); + tab_end.apply_gate_at_end(OpType::ECR, {Qubit(0), Qubit(1)}); + tab_end.apply_gate_at_end(OpType::SX, {Qubit(0)}); + tab_end.apply_gate_at_end(OpType::SXdg, {Qubit(1)}); + tab_end.apply_gate_at_end(OpType::ISWAPMax, {Qubit(0), Qubit(1)}); + REQUIRE(tab == tab_end); } }