From beaf4c59137554942259a40d9c1ff3aa26bd76c4 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Fri, 8 Nov 2024 14:42:47 +0000 Subject: [PATCH 01/19] Copy files over to a new branch just for this feature --- tket/CMakeLists.txt | 3 + tket/include/tket/Clifford/APState.hpp | 146 ++++ tket/include/tket/Converters/Converters.hpp | 31 + tket/src/Clifford/APState.cpp | 861 ++++++++++++++++++++ tket/src/Converters/APStateConverters.cpp | 176 ++++ tket/test/CMakeLists.txt | 1 + tket/test/src/test_APState.cpp | 347 ++++++++ 7 files changed, 1565 insertions(+) create mode 100644 tket/include/tket/Clifford/APState.hpp create mode 100644 tket/src/Clifford/APState.cpp create mode 100644 tket/src/Converters/APStateConverters.cpp create mode 100644 tket/test/src/test_APState.cpp diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt index 41808ece11..e7cfc3a681 100644 --- a/tket/CMakeLists.txt +++ b/tket/CMakeLists.txt @@ -168,9 +168,11 @@ target_sources(tket src/Circuit/SubcircuitFinder.cpp src/Circuit/ThreeQubitConversion.cpp src/Circuit/ToffoliBox.cpp + src/Clifford/APState.cpp src/Clifford/ChoiMixTableau.cpp src/Clifford/SymplecticTableau.cpp src/Clifford/UnitaryTableau.cpp + src/Converters/APStateConverters.cpp src/Converters/ChoiMixTableauConverters.cpp src/Converters/Gauss.cpp src/Converters/PauliGraphConverters.cpp @@ -322,6 +324,7 @@ target_sources(tket include/tket/Circuit/StatePreparation.hpp include/tket/Circuit/ThreeQubitConversion.hpp include/tket/Circuit/ToffoliBox.hpp + include/tket/Clifford/APState.hpp include/tket/Clifford/ChoiMixTableau.hpp include/tket/Clifford/SymplecticTableau.hpp include/tket/Clifford/UnitaryTableau.hpp diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp new file mode 100644 index 0000000000..402106b11c --- /dev/null +++ b/tket/include/tket/Clifford/APState.hpp @@ -0,0 +1,146 @@ +// 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 "tket/OpType/OpType.hpp" +#include "tket/Utils/MatrixAnalysis.hpp" + +namespace tket { + +/** + * APState class for Clifford states with global phase tracking + * + * The "affine with phases" form of a Clifford state from ZX calculus (see + * Kissinger & van de Wetering, "Picturing Quantum Software") represents n-qubit + * Clifford states uniquely with the following data: + * - A binary (k,n) matrix A in reduced row echelon form. The leading columns of + * each row are referred to as "leading qubits", and the others as "free + * qubits". We apply H gates to all free qubits, then for each row of A we apply + * a CX targeted on the leading qubit and control by each other entry in the + * row. + * - A binary (k) vector B describing X gates applied to the leading qubits. + * - A Clifford phase polynomial Phi over the free qubits. Wlog this can be in + * the form of a symmetric, zero-diagonal, binary (n-k,n-k) matrix E describing + * CZ gates and a (n-k) vector P of integers mod 4 describing S gates. + * + * This gives a canonical statevector: + * \sum_{x, Ax=B} i^{\Phi(x)} \ket{x} + * + * This canonical statevector fixes a reference phase from which we can track + * global phase with an additional parameter. + * + * When applying gates, qubits may switch between being leading and free + * (including those that aren't involved in the gate due to the need to reduce + * to the canonical form, e.g. reduced row echelon form for A). The updates are + * easiest to prove in the ZX calculus for the gateset (CZ, S i.e. pi/2 green + * phase, SX i.e. pi/2 red phase). + * + * To save on resizing the matrices and vectors, we will keep them at full size + * and just assert correctness that, e.g. every entry in A is either on the + * diagonal (to indicate the qubit is leading) or in the row of a leading qubit + * and a column of a later following qubit. + */ +class APState { + public: + /** + * An upper triangular matrix whose diagonal entries indicate leading qubits + * and off-diagonal entries represent a CX gate from the column qubit + * (necessarily a free qubit) to the row qubit. + * + * If a diagonal entry is 0 (it is a free qubit), then every entry in the row + * is 0. If a diagonal entry is 1 (it is a leading qubit), all other entries + * in the column are 0. + */ + MatrixXb A; + + /** + * A vector indicating X gates on leading qubits. + * + * Must be 0 on every free qubit. + */ + VectorXb B; + + /** + * A symmetric, zero diagonal matrix whose entries indicate CZs between free + * qubits. + * + * Every diagonal entry must be 0. + * The column and row for each leading qubit must be 0. + */ + MatrixXb E; + + /** + * A vector indicating S^{P(i)} on qubit i which must be free. + * + * Must be 0 on every leading qubit. + */ + Eigen::VectorXi P; + + /** + * The global phase term (in half-turns). + */ + Expr phase; + + /** + * Constructs a state in AP form from given data. + */ + APState( + const MatrixXb& A_, const VectorXb& B_, const MatrixXb& E_, + const Eigen::VectorXi& P_, const Expr& phase_); + + /** + * Constructs the state \ket{0}^{\otimes n_qubits} in AP form. + */ + APState(unsigned n_qubits); + + /** + * Constructs the state in AP form from a given statevector. + */ + APState(const Eigen::VectorXcd& sv); + + /** + * Verifies the internal correctness of the data structure. Throws an + * exception if the data structure is invalid. + */ + void verify(); + + /** + * Calculates the statevector of the state. + */ + Eigen::VectorXcd to_statevector(); + + /** + * Applies a CZ gate to the state. O(n^2) wrt number of qubits in the state. + */ + void apply_CZ(unsigned ctrl, unsigned trgt); + + /** + * Applies an S gate to the state. O(n^2) wrt number of qubits in the state. + */ + void apply_S(unsigned q); + + /** + * Applies a V gate to the state. O(n^2) wrt number of qubits in the state. + */ + void apply_V(unsigned q); + + /** + * Applies an unparameterised Clifford gate to the state on the chosen qubits. + * O(n^2) wrt number of qubits in the state. + */ + void apply_gate(OpType type, const std::vector& qbs); +}; + +} // namespace tket \ No newline at end of file diff --git a/tket/include/tket/Converters/Converters.hpp b/tket/include/tket/Converters/Converters.hpp index a3de8c9146..0a9f9809cc 100644 --- a/tket/include/tket/Converters/Converters.hpp +++ b/tket/include/tket/Converters/Converters.hpp @@ -15,6 +15,7 @@ #pragma once #include "tket/Circuit/Circuit.hpp" +#include "tket/Clifford/APState.hpp" #include "tket/Clifford/ChoiMixTableau.hpp" #include "tket/Clifford/UnitaryTableau.hpp" #include "tket/PauliGraph/PauliGraph.hpp" @@ -45,6 +46,20 @@ Circuit unitary_rev_tableau_to_circuit(const UnitaryRevTableau &tab); */ ChoiMixTableau circuit_to_cm_tableau(const Circuit &circ); +/** + * Construct an APState for a given circuit. + * Assumes all input qubits are initialised in the |0> state. + * Will throw an exception if it contains non-Clifford gates. + */ +APState circuit_to_apstate(const Circuit &circ); + +/** + * Construct a circuit producing the state described by the APState. + * Uses the standard circuit form implied by the ZX description of reduced + * AP-form (layered as X-H-CX-CZ-S). + */ +Circuit apstate_to_circuit(const APState &ap); + /** * Constructs a circuit producing the same effect as a ChoiMixTableau. * Since Circuit does not support distinct qubit addresses for inputs and @@ -130,6 +145,22 @@ std::pair cm_tableau_to_unitary_extension_circuit( const std::vector &post_names = {}, CXConfigType cx_config = CXConfigType::Snake); +/** + * Convert a SymplecticTableau describing a stabiliser state to reduced AP-form + * (i.e. an APState object). Since tableau methods do not track global phase, we + * set it to 0. Throws an exception if the tableau does not describe a pure + * state (i.e. it must have n commuting rows for n qubits). + */ +APState tableau_to_apstate(SymplecticTableau tab); + +/** + * Convert an APState describing a stabiliser state into a tableau form, + * discarding the global phase information. This allows us to reuse tableau + * synthesis methods to synthesise APState objects (the result would need + * re-simulating as an APState to compare and deduce the required global phase). + */ +SymplecticTableau apstate_to_tableau(const APState &ap); + PauliGraph circuit_to_pauli_graph(const Circuit &circ); /** diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp new file mode 100644 index 0000000000..a21e4b446b --- /dev/null +++ b/tket/src/Clifford/APState.cpp @@ -0,0 +1,861 @@ +// 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/Clifford/APState.hpp" + +#include + +#include "tket/OpType/OpTypeInfo.hpp" + +namespace tket { + +APState::APState( + const MatrixXb& A_, const VectorXb& B_, const MatrixXb& E_, + const Eigen::VectorXi& P_, const Expr& phase_) + : A(A_), B(B_), E(E_), P(P_), phase(phase_) { + verify(); +} + +APState::APState(unsigned n_qubits) + : A(MatrixXb::Identity(n_qubits, n_qubits)), + B(VectorXb::Zero(n_qubits)), + E(MatrixXb::Zero(n_qubits, n_qubits)), + P(Eigen::VectorXi(n_qubits)), + phase(0) {} + +unsigned clifford_phase(const Complex& c) { + unsigned res = 0; + unsigned n_possibles = 0; + if (c.real() > EPS) { + res = 0; + ++n_possibles; + } + if (c.imag() > EPS) { + res = 1; + ++n_possibles; + } + if (c.real() < -EPS) { + res = 2; + ++n_possibles; + } + if (c.imag() < -EPS) { + res = 3; + ++n_possibles; + } + TKET_ASSERT(n_possibles == 1); + return res; +} + +APState::APState(const Eigen::VectorXcd& sv) { + unsigned n_qbs = 0; + while (sv.size() > (1 << n_qbs)) ++n_qbs; + TKET_ASSERT(sv.size() == (1 << n_qbs)); + + // Find non-zero entries as a vector space and offset + unsigned z0 = 0; + for (unsigned x = 0; x < sv.size(); ++x) { + if (sv(x) != 0.) { + z0 = x; + break; + } + } + std::vector offsets; + unsigned n_non_zero = 0; + for (unsigned x = 1; x < sv.size(); ++x) { + if (sv(z0 ^ x) != 0.) { + ++n_non_zero; + if (n_non_zero == (1 << offsets.size())) offsets.push_back(x); + } + } + + // Find A as the dual space + MatrixXb offset_mat(offsets.size(), n_qbs); + for (unsigned r = 0; r < offsets.size(); ++r) { + unsigned off = offsets.at(r); + for (unsigned c = 0; c < n_qbs; ++c) { + // Binary encoding of offsets in reverse order to guarantee free qubits + // are the later ones, meaning we produce A in row echelon form + offset_mat(r, c) = (off & (1 << c)) != 0; + } + } + std::vector> row_ops = + gaussian_elimination_row_ops(offset_mat); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + offset_mat(op.second, j) ^= offset_mat(op.first, j); + } + } + std::map mat_leaders; + A = MatrixXb::Zero(n_qbs, n_qbs); + for (unsigned c = 0; c < n_qbs; ++c) { + bool free_qubit = false; + for (unsigned r = 0; r < offsets.size(); ++r) { + if (offset_mat(r, c)) { + auto inserted = mat_leaders.insert({r, c}); + if (inserted.second) { + free_qubit = true; + break; + } else { + // Reverse bit orderings back to normal + A(n_qbs - 1 - c, n_qbs - 1 - inserted.first->second) = true; + } + } + } + A(n_qbs - 1 - c, n_qbs - 1 - c) = !free_qubit; + } + + VectorXb z0_vec(n_qbs); + for (unsigned i = 0; i < n_qbs; ++i) { + z0_vec(i) = ((z0 & (1 << (n_qbs - 1 - i))) != 0); + } + B = A * z0_vec; + + E = MatrixXb::Zero(n_qbs, n_qbs); + P = Eigen::VectorXi::Zero(n_qbs); + unsigned neutral_z0 = z0; // Index with 0 for all free qubits + std::map offset_for_free; + for (const std::pair& row_qfree : mat_leaders) { + unsigned offset = 0; + for (unsigned i = 0; i < n_qbs; ++i) + if (offset_mat(row_qfree.first, i)) offset += (1 << i); + unsigned qfree = n_qbs - 1 - row_qfree.second; + offset_for_free.insert({qfree, offset}); + if ((neutral_z0 & (1 << row_qfree.second)) != 0) + neutral_z0 = neutral_z0 ^ offset; + } + for (const std::pair& qfree_offset : + offset_for_free) { + unsigned qfree = qfree_offset.first; + unsigned offset = qfree_offset.second; + auto complex_phase = sv(neutral_z0 ^ offset) / sv(neutral_z0); + P(qfree) = clifford_phase(complex_phase); + for (const std::pair& qfree_offset2 : + offset_for_free) { + if (qfree_offset == qfree_offset2) break; // Only go up to solved phases + unsigned qfree2 = qfree_offset2.first; + unsigned offset2 = qfree_offset2.second; + auto complex_phase_cz = + sv(neutral_z0 ^ offset ^ offset2) / sv(neutral_z0); + E(qfree, qfree2) = E(qfree2, qfree) = + ((clifford_phase(complex_phase_cz) - P(qfree) - P(qfree2)) % 4) == 2; + } + } + + phase = Expr(std::arg(sv(neutral_z0)) / PI); +} + +void APState::verify() { + // Check dimensions all agree + unsigned n_qubits = A.rows(); + TKET_ASSERT(A.cols() == n_qubits); + TKET_ASSERT(B.size() == n_qubits); + TKET_ASSERT(E.rows() == n_qubits); + TKET_ASSERT(E.cols() == n_qubits); + TKET_ASSERT(P.size() == n_qubits); + + for (unsigned r = 0; r < n_qubits; ++r) { + for (unsigned c = 0; c < r; ++c) { + // Check A is upper triangular + TKET_ASSERT(!A(r, c)); + // Check E is symmetric + TKET_ASSERT(E(r, c) == E(c, r)); + } + if (A(r, r)) { + // Check leaders are eliminated in A (row echelon) + for (unsigned r2 = 0; r2 < r; ++r2) { + TKET_ASSERT(!A(r2, r)); + } + } else { + // Check A only has entries in rows with a leader + for (unsigned c = r + 1; c < n_qubits; ++c) { + TKET_ASSERT(!A(r, c)); + } + } + } + + for (unsigned q = 0; q < n_qubits; ++q) { + if (A(q, q)) { + // Check P is zero on leaders + TKET_ASSERT(P(q) % 4 == 0); + // Check E is zero on leaders (simplified by symmetry) + for (unsigned q2 = 0; q2 < n_qubits; ++q2) { + TKET_ASSERT(!E(q, q2)); + } + } else { + // Check B is zero on free qubits + TKET_ASSERT(!B(q)); + // Check diagonals of E are still zero + TKET_ASSERT(!E(q, q)); + } + } +} + +VectorXb z2_mult(const MatrixXb& M, const VectorXb& v) { + VectorXb res = VectorXb::Zero(M.rows()); + for (unsigned i = 0; i < M.cols(); ++i) { + if (v(i)) { + for (unsigned j = 0; j < M.rows(); ++j) res(j) ^= M(j, i); + } + } + return res; +} + +Eigen::VectorXcd APState::to_statevector() { + unsigned n_qubits = A.cols(); + Eigen::VectorXcd sv = Eigen::VectorXcd::Zero(1 << n_qubits); + unsigned n_terms = 0; + Complex g_phase = std::exp(i_ * PI * eval_expr(phase).value()); + for (unsigned x = 0; x < 1 << n_qubits; ++x) { + VectorXb x_binary = VectorXb::Zero(n_qubits); + for (unsigned i = 0; i < n_qubits; ++i) { + unsigned mask = 1 << (n_qubits - 1 - i); + x_binary(i) = ((x & mask) != 0); + } + + if (z2_mult(A, x_binary) == B) { + ++n_terms; + unsigned i_phases = 0; + for (unsigned q = 0; q < n_qubits; ++q) { + if (x_binary(q)) { + i_phases += P(q); + for (unsigned q2 = q; q2 < n_qubits; ++q2) { + if (E(q, q2) && x_binary(q2)) i_phases += 2; + } + } + } + switch (i_phases % 4) { + case 0: { + sv(x) = g_phase; + break; + } + case 1: { + sv(x) = i_ * g_phase; + break; + } + case 2: { + sv(x) = -g_phase; + break; + } + case 3: { + sv(x) = -i_ * g_phase; + break; + } + default: { + TKET_ASSERT(false); + } + } + } + } + return pow(n_terms, -0.5) * sv; +} + +void APState::apply_CZ(unsigned ctrl, unsigned trgt) { + if (ctrl == trgt) + throw std::logic_error( + "APState::apply_CZ: cannot apply CZ when control and target qubits " + "match"); + unsigned n_qubits = A.cols(); + if (ctrl >= n_qubits || trgt >= n_qubits) + throw std::logic_error("APState::apply_CZ: qubit indices out of range"); + + if (A(ctrl, ctrl)) { + if (A(trgt, trgt)) { + // ctrl and trgt are both leading + + std::set just_ctrl, just_trgt, both; + // Add phases and determine connectivity + for (unsigned q = 0; q < n_qubits; ++q) { + if (A(ctrl, q)) { + if (B(trgt)) P(q) += 2; + if (A(trgt, q)) { + if (!B(ctrl)) P(q) += 2; + both.insert(q); + } else { + just_ctrl.insert(q); + } + } else if (A(trgt, q)) { + if (B(ctrl)) P(q) += 2; + just_trgt.insert(q); + } + } + // ctrl and trgt were included in the previous loop, so reset them + just_ctrl.erase(ctrl); + just_trgt.erase(trgt); + P(ctrl) = 0; + P(trgt) = 0; + // Pivoting-style update between those connected to ctrl, trgt, and both + for (const unsigned q1 : just_ctrl) { + for (const unsigned q2 : just_trgt) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + for (const unsigned q2 : both) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + } + for (const unsigned q1 : just_trgt) { + for (const unsigned q2 : both) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + } + + if (B(ctrl) && B(trgt)) phase += 1; + } else { + // ctrl is leading, trgt is free + + // Flip connectivity between trgt and every q on ctrl + for (unsigned q = ctrl + 1; q < n_qubits; ++q) { + if (A(ctrl, q)) { + E(trgt, q) ^= true; + E(q, trgt) ^= true; + // Note that if q = trgt, this does nothing, preserving the zero + // diagonal of E + } + } + // Add phase to trgt + if (A(ctrl, trgt) ^ B(ctrl)) P(trgt) += 2; + + // No global phase change + } + } else { + if (A(trgt, trgt)) { + // trgt is leading, ctrl is free + + // CZ is symmetric, so reuse the previous case + apply_CZ(trgt, ctrl); + } else { + // ctrl and trgt are both free + + // Just add the CZ gate to E + E(ctrl, trgt) ^= true; + E(trgt, ctrl) ^= true; + + // No global phase change + } + } +} + +void APState::apply_S(unsigned q) { + unsigned n_qubits = A.cols(); + if (q >= n_qubits) + throw std::logic_error("APState::apply_S: qubit index out of range"); + + if (A(q, q)) { + // q is leading + + // Local complementation update to E + for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { + if (A(q, q1)) { + for (unsigned q2 = q + 1; q2 < n_qubits; ++q2) { + E(q1, q2) ^= ((q1 != q2) && A(q, q2)); + } + } + } + + // Update global phase + if (B(q)) phase += .5; + + // Update local phases + unsigned local_phase_change = B(q) ? 3 : 1; + for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { + if (A(q, q1)) P(q1) += local_phase_change; + } + } else { + // q is free + + // Add to local phase + P(q) += 1; + + // No global phase change + } +} + +void APState::apply_V(unsigned q) { + unsigned n_qubits = A.cols(); + if (q >= n_qubits) + throw std::logic_error("APState::apply_V: qubit index out of range"); + + if (A(q, q)) { + // q is leading, but will become free + + // Local complementation for E + for (unsigned q1 = q; q1 < n_qubits; ++q1) { + if (A(q, q1)) { + for (unsigned q2 = q1 + 1; q2 < n_qubits; ++q2) { + E(q1, q2) ^= A(q, q2); + E(q2, q1) ^= A(q, q2); + } + + // Local phases + P(q1) += B(q) ? 1 : 3; + } + } + + // Global phase update + if (B(q)) phase += -.5; + + // q is no longer a leader + for (unsigned q1 = q; q1 < n_qubits; ++q1) { + A(q, q1) = false; + } + B(q) = false; + } else { + // q is free + + std::optional last_connected_leader; + for (unsigned q1 = q; q1 > 0;) { + --q1; + if (A(q1, q)) { + last_connected_leader = q1; + break; + } + } + + if (last_connected_leader.has_value()) { + // q is free and connected to a leader + + // For each other leader connected to q, add LCL's row in A + std::list other_leaders; + for (unsigned q1 = 0; q1 < *last_connected_leader; ++q1) { + if (A(q1, q)) { + other_leaders.push_back(q1); + for (unsigned q2 = *last_connected_leader; q2 < n_qubits; ++q2) { + A(q1, q2) ^= A(*last_connected_leader, q2); + } + } + } + + // Sort the connected frees into just q, just LCL, and both + std::list just_q, just_lcl, both; + for (unsigned q1 = 0; q1 < n_qubits; ++q1) { + if (q1 == q || q1 == *last_connected_leader) continue; + if (E(q, q1)) { + if (A(*last_connected_leader, q1)) + both.push_back(q1); + else + just_q.push_back(q1); + } else if (A(*last_connected_leader, q1)) + just_lcl.push_back(q1); + } + + // Split case by phase on q + // In each case: + // - Perform appropriate complementations between connected frees + // - Set LCL's connectivity in E appropriately + // - Modify q's connectivity in E + // - Add all local and global phases by a case matrix + if (P(q) % 2 == 0) { + // Handle complementations between neighbours + for (const unsigned q1 : just_lcl) { + // Local complement within group + for (const unsigned q2 : just_lcl) { + E(q1, q2) ^= true; + } + E(q1, q1) ^= true; + // Complement between groups + for (const unsigned q2 : just_q) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + // Add connections with q + E(q, q1) ^= true; + E(q1, q) ^= true; + } + for (const unsigned q1 : both) { + // Local complement within group + for (const unsigned q2 : both) { + E(q1, q2) ^= true; + } + E(q1, q1) ^= true; + // Complement between groups + for (const unsigned q2 : just_q) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + } + // Remove connections with q + for (const unsigned q1 : just_q) { + E(q, q1) ^= true; + E(q1, q) ^= true; + } + + // Move LCL from A to E + E(q, *last_connected_leader) = true; + E(*last_connected_leader, q) = true; + for (const unsigned q1 : just_q) { + E(q1, *last_connected_leader) = true; + E(*last_connected_leader, q1) = true; + } + for (const unsigned q1 : just_lcl) { + E(q1, *last_connected_leader) = true; + E(*last_connected_leader, q1) = true; + A(*last_connected_leader, q1) = false; + } + for (const unsigned q1 : both) A(*last_connected_leader, q1) = false; + A(*last_connected_leader, q) = false; + A(*last_connected_leader, *last_connected_leader) = false; + + // Phases + bool a = (P(q) % 4 == 2); + bool b = B(*last_connected_leader); + P(q) = b ? 1 : 3; + P(*last_connected_leader) = (a ^ b) ? 1 : 3; + B(*last_connected_leader) = false; + for (const unsigned q1 : other_leaders) B(q1) ^= b; + for (const unsigned q1 : both) P(q1) += a ? 3 : 1; + for (const unsigned q1 : just_lcl) P(q1) += (b ^ a) ? 1 : 3; + for (const unsigned q1 : just_q) P(q1) += b ? 2 : 0; + if (b) phase += a ? .5 : -.5; + } else { + // Handle complementations between neighbours + for (const unsigned q1 : just_lcl) { + // Complement between groups + for (const unsigned q2 : just_q) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + for (const unsigned q2 : both) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + // Add connections with q + E(q, q1) = true; + E(q1, q) = true; + } + for (const unsigned q1 : just_q) { + // Complement between groups + for (const unsigned q2 : both) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + // Remove connections with q + E(q, q1) = false; + E(q1, q) = false; + } + + // Move LCL from A to E + E(q, *last_connected_leader) = true; + E(*last_connected_leader, q) = true; + for (const unsigned q1 : just_q) { + E(q1, *last_connected_leader) = true; + E(*last_connected_leader, q1) = true; + } + for (const unsigned q1 : both) { + E(q1, *last_connected_leader) = true; + E(*last_connected_leader, q1) = true; + A(*last_connected_leader, q1) = false; + } + for (const unsigned q1 : just_lcl) + A(*last_connected_leader, q1) = false; + A(*last_connected_leader, q) = false; + A(*last_connected_leader, *last_connected_leader) = false; + + // Phases + bool a = (P(q) % 4 == 3); + bool b = B(*last_connected_leader); + P(q) = b ? 1 : 3; + P(*last_connected_leader) = a ? 2 : 0; + B(*last_connected_leader) = false; + for (const unsigned q1 : other_leaders) B(q1) ^= b; + for (const unsigned q1 : just_q) P(q1) += b ? 2 : 0; + for (const unsigned q1 : just_lcl) P(q1) += a ? 2 : 0; + for (const unsigned q1 : both) P(q1) += (a ^ b) ? 0 : 2; + if (a && b) phase += 1.; + } + } else if (P(q) % 2 == 0) { + // q has phase 0/pi and no connected leader + + // Local complementation update to E + std::list connected; + for (unsigned q1 = 0; q1 < n_qubits; ++q1) { + if (E(q, q1)) { + for (const unsigned q2 : connected) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + connected.push_back(q1); + + // Add local phase + P(q1) += (P(q) % 4 == 0) ? 1 : 3; + } + } + + // Local phase on q remains the same + + // Global phase + phase += (P(q) % 4 == 0) ? -0.25 : 0.25; + } else { + // q has phase +-pi/2 and no connected leader + + std::optional first_connected_free = std::nullopt; + for (unsigned q1 = 0; q1 < q; ++q1) { + if (E(q, q1)) { + first_connected_free = q1; + break; + } + } + + if (first_connected_free.has_value()) { + // q has phase +-pi/2, no connected leader, but a previous free which + // will become the new leader + + // Make FCF a leader + for (unsigned q1 = *first_connected_free; q1 < n_qubits; ++q1) { + if (E(q, q1)) { + A(*first_connected_free, q1) = true; + // q ends up disconnected in E + E(q, q1) = false; + E(q1, q) = false; + } + } + A(*first_connected_free, q) = true; + B(*first_connected_free) = (P(q) % 4 == 3); + + // Set phase of q to -pi/2 + P(q) = 3; + + // No global phase change + + // Make A row reduced + for (unsigned q1 = 0; q1 < *first_connected_free; ++q1) { + if (A(q1, *first_connected_free)) { + for (unsigned q2 = *first_connected_free; q2 < n_qubits; ++q2) { + A(q1, q2) ^= A(*first_connected_free, q2); + } + B(q1) ^= B(*first_connected_free); + } + } + + // Need to apply rules for S and CZ gates on FCF to remove + unsigned s_gates = P(*first_connected_free); + P(*first_connected_free) = 0; + std::list cz_targets; + E(*first_connected_free, q) = false; + E(q, *first_connected_free) = false; + for (unsigned q1 = 0; q1 < n_qubits; ++q1) { + if (E(*first_connected_free, q1)) { + cz_targets.push_back(q1); + E(*first_connected_free, q1) = false; + E(q1, *first_connected_free) = false; + } + } + for (unsigned i = 0; i < s_gates; ++i) apply_S(*first_connected_free); + for (const unsigned trgt : cz_targets) + apply_CZ(*first_connected_free, trgt); + } else { + // q has phase +-pi/2, no connected leader, and no previous free, so q + // will become the new leader + A(q, q) = true; + + std::list connected; + for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { + if (E(q, q1)) { + // Local complementation + for (const unsigned q2 : connected) { + E(q1, q2) ^= true; + E(q2, q1) ^= true; + } + connected.push_back(q1); + + // Add to A + A(q, q1) = true; + + // Reset E rows + E(q, q1) = false; + E(q1, q) = false; + } + } + + // Update local and global phases + if (P(q) % 4 == 1) { + for (const unsigned q1 : connected) P(q1) += 3; + } else { + for (const unsigned q1 : connected) P(q1) += 1; + B(q) = true; + phase += -.5; + } + P(q) = 0; + } + } + } +} + +void APState::apply_gate(OpType type, const std::vector& qbs) { + switch (type) { + case OpType::Z: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::X: { + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + phase += .5; + break; + } + case OpType::Y: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + phase += 1.; + break; + } + case OpType::S: { + apply_S(qbs.at(0)); + break; + } + case OpType::Sdg: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::V: { + apply_V(qbs.at(0)); + break; + } + case OpType::SX: { + apply_V(qbs.at(0)); + phase += .25; + break; + } + case OpType::Vdg: { + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + phase += 1.; + break; + } + case OpType::SXdg: { + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + phase += .75; + break; + } + case OpType::H: { + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::CX: { + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + break; + } + case OpType::CY: { + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(1)); + apply_V(qbs.at(1)); + apply_V(qbs.at(1)); + phase += 1.; + break; + } + case OpType::CZ: { + apply_CZ(qbs.at(0), qbs.at(1)); + break; + } + case OpType::ZZMax: { + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + phase -= .25; + break; + } + case OpType::ECR: { + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + phase += .25; + break; + } + case OpType::ISWAPMax: { + apply_V(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(0)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(0)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + phase += 1.; + break; + } + case OpType::SWAP: { + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::BRIDGE: { + apply_S(qbs.at(2)); + apply_V(qbs.at(2)); + apply_S(qbs.at(2)); + apply_CZ(qbs.at(0), qbs.at(2)); + apply_S(qbs.at(2)); + apply_V(qbs.at(2)); + apply_S(qbs.at(2)); + break; + } + case OpType::noop: { + break; + } + case OpType::Phase: { + throw std::logic_error( + "OpType::Phase cannot be applied via APState::apply_gate"); + } + default: { + throw BadOpType( + "Cannot be applied to a APState: not a Clifford gate", type); + } + } +} + +} // namespace tket \ No newline at end of file diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp new file mode 100644 index 0000000000..a56e96da79 --- /dev/null +++ b/tket/src/Converters/APStateConverters.cpp @@ -0,0 +1,176 @@ +// 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/Converters.hpp" + +namespace tket { + +APState circuit_to_apstate(const Circuit& circ) { + APState aps(circ.n_qubits()); + std::map qb_ordering; + for (const Qubit& q : circ.all_qubits()) + qb_ordering.insert({q, qb_ordering.size()}); + for (const Command& com : circ) { + auto args = com.get_args(); + std::vector qbs; + for (const UnitID& q : args) qbs.push_back(qb_ordering.at(q)); + aps.apply_gate(com.get_op_ptr()->get_type(), qbs); + } + return aps; +} + +Circuit apstate_to_circuit(const APState& ap) { + unsigned n_qbs = ap.A.rows(); + Circuit circ(n_qbs); + circ.qubit_create_all(); + for (unsigned q = 0; q < n_qbs; ++q) { + if (!ap.A(q, q)) + circ.add_op(OpType::H, {q}); + else if (ap.B(q)) + circ.add_op(OpType::X, {q}); + } + for (unsigned trgt = 0; trgt < n_qbs; ++trgt) { + if (ap.A(trgt, trgt)) { + for (unsigned ctrl = trgt + 1; ctrl < n_qbs; ++ctrl) + if (ap.A(trgt, ctrl)) circ.add_op(OpType::CX, {ctrl, trgt}); + } + } + for (unsigned q1 = 0; q1 < n_qbs; ++q1) { + if (ap.A(q1, q1)) continue; + for (unsigned q2 = q1 + 1; q2 < n_qbs; ++q2) { + if (ap.E(q1, q2)) circ.add_op(OpType::CZ, {q1, q2}); + } + switch (ap.P(q1) % 4) { + case 1: { + circ.add_op(OpType::S, {q1}); + break; + } + case 2: { + circ.add_op(OpType::Z, {q1}); + break; + } + case 3: { + circ.add_op(OpType::Sdg, {q1}); + break; + } + default: { + break; + } + } + } + circ.add_phase(ap.phase); + return circ; +} + +APState tableau_to_apstate(SymplecticTableau tab) { + unsigned n_qbs = tab.get_n_qubits(); + if (tab.get_n_rows() != n_qbs) + throw std::logic_error( + "tableau_to_apstate requires a tableau with n commuting rows for n " + "qubits"); + MatrixXb fullmat = MatrixXb::Zero(n_qbs, 2 * n_qbs); + /** + * Gaussian elimination by the x matrix first ensures the bottom rows are only + * Zs, i.e. describing rows of A. Reversing the columns of the x matrix + * guarantees that each row has an X on at most one free qubit, simplifying + * the code for finding E and P + */ + for (unsigned c = 0; c < n_qbs; ++c) { + fullmat.col(c) = tab.xmat.col(n_qbs - 1 - c); + } + fullmat.block(0, n_qbs, n_qbs, n_qbs) = tab.zmat; + std::vector> row_ops = + gaussian_elimination_row_ops(fullmat); + for (const std::pair& op : row_ops) + tab.row_mult(op.first, op.second); + + MatrixXb A = MatrixXb::Zero(n_qbs, n_qbs); + VectorXb B = VectorXb::Zero(n_qbs); + MatrixXb E = MatrixXb::Zero(n_qbs, n_qbs); + Eigen::VectorXi P = Eigen::VectorXi::Zero(n_qbs); + + unsigned n_leading = 0; + for (unsigned r = n_qbs; r > 0;) { + --r; + bool is_a_row = true; + for (unsigned c = 0; c < n_qbs; ++c) { + if (tab.xmat(r, c)) { + is_a_row = false; + break; + } + } + if (!is_a_row) break; + unsigned first = 0; + for (unsigned c = 0; c < n_qbs; ++c) { + if (tab.zmat(r, c)) { + first = c; + break; + } + } + A.row(first) = tab.zmat.row(r); + B(first) = tab.phase(r); + ++n_leading; + } + /** + * Each free qubit is after all leaders connected to it, by reduced + * row-echelon of A. Therefore Gaussian elimination of the x matrix in reverse + * order gives rows corresponding to the columns of A of free qubits (plus the + * free qubit itself). + * + * Then the corresponding row of the z matrix is exactly the free qubit's row + * of E, from which we take the diagonal and phase vector to determine P. + */ + for (unsigned r = 0; r < n_qbs - n_leading; ++r) { + unsigned free = 0; + for (unsigned c = n_qbs; c > 0;) { + --c; + if (tab.xmat(r, c)) { + free = c; + break; + } + } + E.row(free) = tab.zmat.row(r); + if (tab.zmat(r, free)) P(free) += 1; + if (tab.phase(r)) P(free) += 2; + E(free, free) = false; + } + + return APState(A, B, E, P, 0); +} + +SymplecticTableau apstate_to_tableau(const APState& ap) { + unsigned n_qbs = ap.A.rows(); + MatrixXb xmat = MatrixXb::Zero(n_qbs, n_qbs); + MatrixXb zmat = MatrixXb::Zero(n_qbs, n_qbs); + VectorXb phase = VectorXb::Zero(n_qbs); + + for (unsigned q = 0; q < n_qbs; ++q) { + if (ap.A(q, q)) { + // Leading qubit + zmat.row(q) = ap.A.row(q); + phase(q) = ap.B(q); + } else { + // Free qubit + xmat.row(q) = ap.A.col(q); + xmat(q, q) = true; + zmat.row(q) = ap.E.row(q); + zmat(q, q) = (ap.P(q) % 2 == 1); + phase(q) = (ap.P(q) % 4 > 1); + } + } + + return SymplecticTableau(xmat, zmat, phase); +} + +} // namespace tket \ No newline at end of file diff --git a/tket/test/CMakeLists.txt b/tket/test/CMakeLists.txt index 4e020c4d16..fe7011d6c5 100644 --- a/tket/test/CMakeLists.txt +++ b/tket/test/CMakeLists.txt @@ -102,6 +102,7 @@ add_executable(test-tket src/Simulation/test_CircuitSimulator.cpp src/Simulation/test_PauliExpBoxUnitaryCalculator.cpp src/test_AASRoute.cpp + src/test_APState.cpp src/test_ArchitectureAwareSynthesis.cpp src/test_Architectures.cpp src/test_Assertion.cpp diff --git a/tket/test/src/test_APState.cpp b/tket/test/src/test_APState.cpp new file mode 100644 index 0000000000..12320e8a90 --- /dev/null +++ b/tket/test/src/test_APState.cpp @@ -0,0 +1,347 @@ +// 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 "Simulation/ComparisonFunctions.hpp" +#include "testutil.hpp" +#include "tket/Circuit/Simulation/CircuitSimulator.hpp" +#include "tket/Clifford/APState.hpp" +#include "tket/Converters/Converters.hpp" + +namespace tket { +namespace test_APState { + +void test_apply_gate( + const MatrixXb& A, const VectorXb& B, const MatrixXb& E, + const Eigen::VectorXi& P, OpType ot, const std::vector& args) { + APState ap(A, B, E, P, 0); + ap.verify(); + auto sv_before = ap.to_statevector(); + + Circuit circ(A.rows()); + circ.add_op(ot, args); + auto gate_u = tket_sim::get_unitary(circ); + + ap.apply_gate(ot, args); + + ap.verify(); + auto sv_after = ap.to_statevector(); + + CHECK(tket_sim::compare_statevectors_or_unitaries( + gate_u * sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL)); +} + +SCENARIO("CZ cases") { + GIVEN("CZ on free qubits") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 3) = true; + E(2, 3) = E(3, 2) = true; + test_apply_gate(A, B, E, P, OpType::CZ, {1, 2}); + } + GIVEN("CZ on one leading qubit and connected free") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = A(1, 4) = true; + B(1) = (b == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 3}); + } + } + GIVEN("CZ on one leading qubit and unconnected free") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = A(1, 4) = true; + B(1) = (b == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 4}); + } + } + GIVEN("CZ on leading qubits") { + for (unsigned b1 = 0; b1 < 2; ++b1) { + for (unsigned b2 = 0; b2 < 2; ++b2) { + MatrixXb A = MatrixXb::Zero(8, 8); + VectorXb B = VectorXb::Zero(8); + MatrixXb E = MatrixXb::Zero(8, 8); + Eigen::VectorXi P = Eigen::VectorXi::Zero(8); + A(0, 0) = A(0, 2) = A(0, 3) = A(0, 4) = A(0, 5) = true; + A(1, 1) = A(1, 4) = A(1, 5) = A(1, 6) = A(1, 7) = true; + B(0) = (b1 == 1); + B(1) = (b2 == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 1}); + } + } + } +} + +SCENARIO("S cases") { + GIVEN("S on free qubit") { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + A(0, 0) = A(0, 1) = A(0, 2) = true; + E(1, 2) = E(2, 1) = true; + test_apply_gate(A, B, E, P, OpType::S, {2}); + } + GIVEN("S on leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 2) = true; + B(0) = (b == 1); + E(1, 3) = E(3, 1) = true; + test_apply_gate(A, B, E, P, OpType::S, {0}); + } + } + GIVEN("S on disconnected leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(1, 1); + VectorXb B = VectorXb::Zero(1); + MatrixXb E = MatrixXb::Zero(1, 1); + Eigen::VectorXi P = Eigen::VectorXi::Zero(1); + A(0, 0) = true; + B(0) = (b == 1); + test_apply_gate(A, B, E, P, OpType::S, {0}); + } + } +} + +SCENARIO("V cases") { + GIVEN("V on leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = true; + B(0) = (b == 1); + test_apply_gate(A, B, E, P, OpType::V, {0}); + } + } + GIVEN("V on free qubit with some leading") { + for (unsigned b = 0; b < 2; ++b) { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(9, 9); + VectorXb B = VectorXb::Zero(9); + MatrixXb E = MatrixXb::Zero(9, 9); + Eigen::VectorXi P = Eigen::VectorXi::Zero(9); + A(0, 0) = A(0, 2) = A(0, 4) = A(0, 5) = A(0, 7) = true; + A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = A(1, 5) = A(1, 6) = true; + B(1) = (b == 1); + E(4, 5) = E(5, 4) = true; + E(4, 6) = E(6, 4) = true; + E(4, 7) = E(7, 4) = true; + E(4, 8) = E(8, 4) = true; + P(4) = p; + test_apply_gate(A, B, E, P, OpType::V, {4}); + } + } + } + GIVEN("V on free qubit with some earlier connected free") { + for (unsigned p1 = 0; p1 < 4; ++p1) { + for (unsigned p2 = 0; p2 < 4; ++p2) { + MatrixXb A = MatrixXb::Zero(9, 9); + VectorXb B = VectorXb::Zero(9); + MatrixXb E = MatrixXb::Zero(9, 9); + Eigen::VectorXi P = Eigen::VectorXi::Zero(9); + A(0, 0) = true; + A(1, 1) = A(1, 4) = A(1, 6) = true; + A(2, 2) = A(2, 4) = A(2, 7) = true; + A(3, 3) = A(3, 4) = A(3, 8) = A(3, 8) = true; + E(4, 5) = E(5, 4) = true; + E(4, 7) = E(7, 4) = true; + E(4, 8) = E(8, 4) = true; + E(5, 6) = E(6, 5) = true; + E(5, 7) = E(7, 5) = true; + E(5, 8) = E(8, 5) = true; + P(4) = p1; + P(5) = p2; + test_apply_gate(A, B, E, P, OpType::V, {5}); + } + } + } + GIVEN("V on free qubit with no earlier connected free") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = true; + E(1, 2) = E(2, 1) = true; + E(1, 3) = E(3, 1) = true; + P(1) = p; + test_apply_gate(A, B, E, P, OpType::V, {1}); + } + } + GIVEN("V on disconnected free qubit") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(1, 1); + VectorXb B = VectorXb::Zero(1); + MatrixXb E = MatrixXb::Zero(1, 1); + Eigen::VectorXi P = Eigen::VectorXi::Zero(1, 1); + P(0) = p; + test_apply_gate(A, B, E, P, OpType::V, {0}); + } + } +} + +SCENARIO("Gate encodings") { + std::list>> test_gates = { + {OpType::Z, {0}}, {OpType::X, {0}}, + {OpType::Y, {0}}, {OpType::S, {0}}, + {OpType::Sdg, {0}}, {OpType::V, {0}}, + {OpType::Vdg, {0}}, {OpType::SX, {0}}, + {OpType::SXdg, {0}}, {OpType::H, {0}}, + {OpType::CX, {0, 1}}, {OpType::CY, {0, 1}}, + {OpType::CZ, {0, 1}}, {OpType::ZZMax, {0, 1}}, + {OpType::ECR, {0, 1}}, {OpType::ISWAPMax, {0, 1}}, + {OpType::SWAP, {0, 1}}, {OpType::BRIDGE, {0, 1, 2}}, + {OpType::noop, {0}}, + }; + GIVEN("Check Z actions") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Identity(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + test_apply_gate(A, B, E, P, com.first, com.second); + } + } + GIVEN("Check X actions") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + test_apply_gate(A, B, E, P, com.first, com.second); + } + } +} + +SCENARIO("Loading from a statevector") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 2) = true; + B(0) = true; + E(2, 3) = E(3, 2) = true; + P(2) = 1; + P(3) = 2; + APState ap(A, B, E, P, 0); + auto sv = ap.to_statevector(); + APState reconstructed(sv); + auto sv2 = reconstructed.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv, sv2, tket_sim::MatrixEquivalence::EQUAL)); +} + +SCENARIO("Converting from/to a circuit") { + GIVEN("A circuit in the standard AP form") { + Circuit circ(4); + circ.qubit_create_all(); + circ.add_op(OpType::H, {2}); + circ.add_op(OpType::H, {3}); + circ.add_op(OpType::CX, {2, 0}); + circ.add_op(OpType::CX, {2, 1}); + circ.add_op(OpType::CX, {3, 1}); + circ.add_op(OpType::S, {2}); + circ.add_op(OpType::Z, {3}); + APState ap = circuit_to_apstate(circ); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + Circuit reconstructed = apstate_to_circuit(ap); + CHECK(circ == reconstructed); + } + GIVEN("A generic circuit") { + Circuit circ(4); + circ.qubit_create_all(); + circ.add_op(OpType::V, {0}); + circ.add_op(OpType::CX, {0, 1}); + circ.add_op(OpType::CY, {1, 3}); + circ.add_op(OpType::H, {3}); + circ.add_op(OpType::ZZMax, {2, 3}); + APState ap = circuit_to_apstate(circ); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + Circuit reconstructed = apstate_to_circuit(ap); + auto sv_rec = tket_sim::get_statevector(reconstructed); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_rec, tket_sim::MatrixEquivalence::EQUAL)); + } +} + +SCENARIO("Converting from/to a tableau") { + GIVEN("Check up to global phase using circuit") { + Circuit circ(8); + circ.qubit_create_all(); + circ.add_op(OpType::X, {1}); + circ.add_op(OpType::X, {5}); + circ.add_op(OpType::H, {2}); + circ.add_op(OpType::H, {4}); + circ.add_op(OpType::H, {6}); + circ.add_op(OpType::H, {7}); + circ.add_op(OpType::CX, {2, 1}); + circ.add_op(OpType::CX, {4, 0}); + circ.add_op(OpType::CX, {4, 3}); + circ.add_op(OpType::CX, {6, 0}); + circ.add_op(OpType::CX, {6, 1}); + circ.add_op(OpType::CX, {7, 5}); + circ.add_op(OpType::CZ, {2, 6}); + circ.add_op(OpType::CZ, {4, 6}); + circ.add_op(OpType::CZ, {4, 7}); + circ.add_op(OpType::CZ, {6, 7}); + circ.add_op(OpType::S, {2}); + circ.add_op(OpType::Sdg, {4}); + circ.add_op(OpType::Z, {7}); + ChoiMixTableau cmt = circuit_to_cm_tableau(circ); + APState ap = tableau_to_apstate(cmt.tab_); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + SymplecticTableau tab2 = apstate_to_tableau(ap); + ChoiMixTableau cmt2(tab2.xmat, tab2.zmat, tab2.phase); + auto [circ2, perm] = cm_tableau_to_exact_circuit(cmt2); + qubit_map_t inv; + for (const std::pair& qp : perm) + inv.insert({qp.second, qp.first}); + circ2.permute_boundary_output(inv); + auto sv_circ2 = tket_sim::get_statevector(circ2); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_circ2, + tket_sim::MatrixEquivalence::EQUAL_UP_TO_GLOBAL_PHASE)); + } +} + +} // namespace test_APState +} // namespace tket \ No newline at end of file From 8dfb97a319c4eb37f17dabbf471f8f57c7c6e86e Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Wed, 13 Nov 2024 09:12:34 +0000 Subject: [PATCH 02/19] Attempt signed shift fix --- tket/src/Clifford/APState.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp index a21e4b446b..2dfa489611 100644 --- a/tket/src/Clifford/APState.cpp +++ b/tket/src/Clifford/APState.cpp @@ -75,7 +75,7 @@ APState::APState(const Eigen::VectorXcd& sv) { for (unsigned x = 1; x < sv.size(); ++x) { if (sv(z0 ^ x) != 0.) { ++n_non_zero; - if (n_non_zero == (1 << offsets.size())) offsets.push_back(x); + if (n_non_zero == (1u << offsets.size())) offsets.push_back(x); } } @@ -216,7 +216,7 @@ Eigen::VectorXcd APState::to_statevector() { Eigen::VectorXcd sv = Eigen::VectorXcd::Zero(1 << n_qubits); unsigned n_terms = 0; Complex g_phase = std::exp(i_ * PI * eval_expr(phase).value()); - for (unsigned x = 0; x < 1 << n_qubits; ++x) { + for (unsigned x = 0; x < (1u << n_qubits); ++x) { VectorXb x_binary = VectorXb::Zero(n_qubits); for (unsigned i = 0; i < n_qubits; ++i) { unsigned mask = 1 << (n_qubits - 1 - i); From e51bd2e09f2bdf90060a0b93efad18cd9a168ae1 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Wed, 13 Nov 2024 09:12:44 +0000 Subject: [PATCH 03/19] Bump tket version --- tket/conanfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tket/conanfile.py b/tket/conanfile.py index 147e674cb2..31633eef95 100644 --- a/tket/conanfile.py +++ b/tket/conanfile.py @@ -23,7 +23,7 @@ class TketConan(ConanFile): name = "tket" - version = "1.3.39" + version = "1.3.40" package_type = "library" license = "Apache 2" homepage = "https://github.com/CQCL/tket" From dc60ac30a8230546bfd7661c445b8fde9b51dad8 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Wed, 13 Nov 2024 09:21:56 +0000 Subject: [PATCH 04/19] Forgot to change the version number in pytket, oops! --- pytket/conanfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytket/conanfile.py b/pytket/conanfile.py index e8b79f648b..9b76b1cb1a 100644 --- a/pytket/conanfile.py +++ b/pytket/conanfile.py @@ -38,7 +38,7 @@ def requirements(self): self.requires("pybind11_json/0.2.14") self.requires("symengine/0.13.0") self.requires("tkassert/0.3.4@tket/stable") - self.requires("tket/1.3.46@tket/stable") + self.requires("tket/1.3.47@tket/stable") self.requires("tklog/0.3.3@tket/stable") self.requires("tkrng/0.3.3@tket/stable") self.requires("tktokenswap/0.3.9@tket/stable") From f1881a743614c3a5495ff57a02f480d5ba46f5d7 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 14 Nov 2024 09:12:22 +0000 Subject: [PATCH 05/19] Doxygen and type warnings fix --- tket/include/tket/Clifford/APState.hpp | 2 +- tket/src/Converters/APStateConverters.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp index 402106b11c..4d9577c4e8 100644 --- a/tket/include/tket/Clifford/APState.hpp +++ b/tket/include/tket/Clifford/APState.hpp @@ -36,7 +36,7 @@ namespace tket { * CZ gates and a (n-k) vector P of integers mod 4 describing S gates. * * This gives a canonical statevector: - * \sum_{x, Ax=B} i^{\Phi(x)} \ket{x} + * \sum_{x, Ax=B} i^{\Phi(x)} |{x}> * * This canonical statevector fixes a reference phase from which we can track * global phase with an additional parameter. diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp index a56e96da79..d73dcbc53d 100644 --- a/tket/src/Converters/APStateConverters.cpp +++ b/tket/src/Converters/APStateConverters.cpp @@ -20,7 +20,7 @@ APState circuit_to_apstate(const Circuit& circ) { APState aps(circ.n_qubits()); std::map qb_ordering; for (const Qubit& q : circ.all_qubits()) - qb_ordering.insert({q, qb_ordering.size()}); + qb_ordering.insert({q, (unsigned)qb_ordering.size()}); for (const Command& com : circ) { auto args = com.get_args(); std::vector qbs; From 8d7dc861c3c61bb1bfdf6ea4da41a413914d5623 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 14 Nov 2024 09:40:24 +0000 Subject: [PATCH 06/19] Missed a second doxygen error --- tket/include/tket/Clifford/APState.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp index 4d9577c4e8..63b90ecef2 100644 --- a/tket/include/tket/Clifford/APState.hpp +++ b/tket/include/tket/Clifford/APState.hpp @@ -36,7 +36,7 @@ namespace tket { * CZ gates and a (n-k) vector P of integers mod 4 describing S gates. * * This gives a canonical statevector: - * \sum_{x, Ax=B} i^{\Phi(x)} |{x}> + * \sum_{x, Ax=B} i^{\Phi(x)} |x> * * This canonical statevector fixes a reference phase from which we can track * global phase with an additional parameter. @@ -101,7 +101,7 @@ class APState { const Eigen::VectorXi& P_, const Expr& phase_); /** - * Constructs the state \ket{0}^{\otimes n_qubits} in AP form. + * Constructs the state |0>^{\otimes n_qubits} in AP form. */ APState(unsigned n_qubits); From 9d486910a0deb53d95965584b572e782abee1dfd Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 14 Nov 2024 12:25:49 +0000 Subject: [PATCH 07/19] Wow, doxygen really hates maths, eh? --- tket/include/tket/Clifford/APState.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp index 63b90ecef2..b1eb808371 100644 --- a/tket/include/tket/Clifford/APState.hpp +++ b/tket/include/tket/Clifford/APState.hpp @@ -101,7 +101,7 @@ class APState { const Eigen::VectorXi& P_, const Expr& phase_); /** - * Constructs the state |0>^{\otimes n_qubits} in AP form. + * Constructs the state |0>^{x n_qubits} in AP form. */ APState(unsigned n_qubits); From 4f9b65854d084b8851e260933a230bb3d0177bc2 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 14:48:44 +0000 Subject: [PATCH 08/19] Add new version as new files --- tket/CMakeLists.txt | 7 +- tket/include/tket/Clifford/APState2.hpp | 369 +++++ tket/include/tket/Converters/Converters.hpp | 77 +- tket/src/Clifford/APState2.cpp | 1470 +++++++++++++++++++ tket/src/Clifford/ChoiMixTableau.cpp | 2 +- tket/src/Converters/APStateConverters2.cpp | 433 ++++++ tket/test/CMakeLists.txt | 2 +- tket/test/src/test_APState2.cpp | 514 +++++++ 8 files changed, 2849 insertions(+), 25 deletions(-) create mode 100644 tket/include/tket/Clifford/APState2.hpp create mode 100644 tket/src/Clifford/APState2.cpp create mode 100644 tket/src/Converters/APStateConverters2.cpp create mode 100644 tket/test/src/test_APState2.cpp diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt index e7cfc3a681..4656592384 100644 --- a/tket/CMakeLists.txt +++ b/tket/CMakeLists.txt @@ -168,11 +168,11 @@ target_sources(tket src/Circuit/SubcircuitFinder.cpp src/Circuit/ThreeQubitConversion.cpp src/Circuit/ToffoliBox.cpp - src/Clifford/APState.cpp + src/Clifford/APState2.cpp src/Clifford/ChoiMixTableau.cpp src/Clifford/SymplecticTableau.cpp src/Clifford/UnitaryTableau.cpp - src/Converters/APStateConverters.cpp + src/Converters/APStateConverters2.cpp src/Converters/ChoiMixTableauConverters.cpp src/Converters/Gauss.cpp src/Converters/PauliGraphConverters.cpp @@ -324,7 +324,8 @@ target_sources(tket include/tket/Circuit/StatePreparation.hpp include/tket/Circuit/ThreeQubitConversion.hpp include/tket/Circuit/ToffoliBox.hpp - include/tket/Clifford/APState.hpp + # include/tket/Clifford/APState.hpp + include/tket/Clifford/APState2.hpp include/tket/Clifford/ChoiMixTableau.hpp include/tket/Clifford/SymplecticTableau.hpp include/tket/Clifford/UnitaryTableau.hpp diff --git a/tket/include/tket/Clifford/APState2.hpp b/tket/include/tket/Clifford/APState2.hpp new file mode 100644 index 0000000000..3a3460b9bb --- /dev/null +++ b/tket/include/tket/Clifford/APState2.hpp @@ -0,0 +1,369 @@ +// 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 "tket/OpType/OpType.hpp" +#include "tket/Utils/MatrixAnalysis.hpp" + +namespace tket { + +/** + * APState class, giving a unique form for a (possibly mixed) Clifford state, + * through which we can track global phase when appropriate. + * + * The "affine with phases" form of a Clifford state from ZX calculus (see + * Kissinger & van de Wetering, "Picturing Quantum Software") represents n-qubit + * Clifford states uniquely with the following data: + * - A binary (n,n) matrix A. + * - A binary (n) vector B. + * - A symmetric, zero-diagonal, binary (n,n) matrix E. + * - A n vector P of integers mod 4 describing S gates. + * + * This gives a canonical statevector (up to a normalisation scalar): + * \sum_{x, Ax=B} i^{\Phi(x)} |x> + * + * This canonical statevector fixes a reference phase from which we can track + * global phase with an additional parameter. + * + * We can generalise this to mixed qubit states by adding a binary (n,n) matrix + * C, with which we now have the canonical density matrix (up to a normalisation + * scalar): + * \sum_{x_1, x_2; Ax_1 = B = Ax_2, Cx_1 = Cx_2} + * i^{x_1^T E x_1 + P^T x_1} |x_1>^{x n_qubits} in AP form. + */ + APState(unsigned n_qubits); + + /** + * Constructs the state in AP form from a given statevector. + */ + APState(const Eigen::VectorXcd& sv); + + /** + * Constructs the state in AP form from a given density matrix. + */ + APState(const Eigen::MatrixXcd& dm); + + /** + * Other required constructors + */ + APState(const APState& other) = default; + APState(APState&& other) = default; + APState& operator=(const APState& other) = default; + APState& operator=(APState& other) = default; + + /** + * Verifies the internal correctness of the data structure. Throws an + * exception if the data structure is invalid. + */ + void verify() const; + + /** + * Calculates the statevector of the state. + * + * Throws an exception if C is non-zero (i.e. it represents a mixed state). + */ + Eigen::VectorXcd to_statevector() const; + + /** + * Calculates the density matrix of the state. + */ + Eigen::MatrixXcd to_density_matrix() const; + + /** + * Applies a CZ gate to the state. O(1). + */ + void apply_CZ(unsigned ctrl, unsigned trgt); + + /** + * Applies an S gate to the state. O(1). + */ + void apply_S(unsigned q); + + /** + * Applies a V gate to the state. O(n^2) wrt number of qubits in the state. + */ + void apply_V(unsigned q); + + /** + * Applies an X gate to the state, faster than applying V twice. O(n) wrt + * number of qubits in the state. + */ + void apply_X(unsigned q); + + /** + * Applies an unparameterised Clifford gate to the state on the chosen qubits. + * O(n^2) wrt number of qubits in the state. + */ + void apply_gate(OpType type, const std::vector& qbs); + + /** + * Initialises a new qubit in the |0> state. Returns the index of the new + * qubit. O(n) wrt number of qubits in the state. + */ + unsigned init_qubit(); + + /** + * Post-selects the chosen qubit to <0|, removing it from the state. Moves the + * final qubit into the place of q, returning its old index. O(n^3) wrt + * number of qubits in the state. + */ + unsigned post_select(unsigned q); + + /** + * Collapses the given qubit in the Z basis. O(n^3) wrt number of qubits in + * the state. + */ + void collapse_qubit(unsigned q); + + /** + * Discards the given qubit, removing it from the state. Moves the final qubit + * into the place of q, returning its old index. O(n^3) wrt number of qubits + * in the state. + */ + unsigned discard_qubit(unsigned q); + + /** + * Manipulates the state into the following normal form: + * - A is in reduced row-echelon form. We refer to the leading columns of each + * row of A as "leading qubits". + * - C is in reduced row-echelon form, and furthermore is zero in any column + * of a leading qubit. We refer to the leading columns of each row of C as + * "mixed qubits". + * - Each entry of E is either between a mixed qubit and a free qubit, or + * between two free qubits. + * - For each leading or mixed qubit, the index of P is zero. + * + * Takes time O(n^4) wrt number of qubits in the state. + */ + void normal_form(); + + bool operator==(const APState& other) const; +}; + +/** + * A wrapper for APState to provide the same interface as ChoiMixTableau wrt + * qubit indexing and distinction between input vs output. + * + * When applying gates, the methods automatically transpose anything occurring + * over the input subspace according to the CJ-isomorphism. + */ +class ChoiAPState { + public: + enum class TableauSegment { Input, Output }; + typedef std::pair col_key_t; + typedef boost::bimap tableau_col_index_t; + + /** + * The internal APState. + */ + APState ap_; + + /** + * Map between column indices and the corresponding qubit ID and type. + */ + tableau_col_index_t col_index_; + + /** + * Construct the identity unitary over n qubits (given default qubit names). + */ + explicit ChoiAPState(unsigned n); + /** + * Construct the identity unitary over specific qubits. + */ + explicit ChoiAPState(const qubit_vector_t& qbs); + /** + * Construct the APState from the underlying binary matrices. + * Qubits are given default names and mapped such that the first columns are + * inputs and the last columns are outputs. + */ + explicit ChoiAPState( + const MatrixXb& A, const VectorXb& b, const MatrixXb& C, + const MatrixXb& E, const Eigen::VectorXi& P, const Expr& phase = 0, + unsigned n_ins = 0); + /** + * Other required constructors + */ + ChoiAPState(const ChoiAPState& other) = default; + ChoiAPState(ChoiAPState&& other) = default; + ChoiAPState& operator=(const ChoiAPState& other) = default; + ChoiAPState& operator=(ChoiAPState& other) = default; + + /** + * Get the total number of qubits/boundaries (sum of inputs and outputs) of + * the process. + */ + unsigned get_n_boundaries() const; + + /** + * Get the number of boundaries representing inputs to the process. + */ + unsigned get_n_inputs() const; + + /** + * Get the number of boundaries representing outputs to 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; + + /** + * Transforms the state according to consuming a Clifford gate at either end + * of the circuit. Applying a gate to the inputs really consumes its transpose + * via the CJ isomorphism. This method handles the transposition internally by + * reversing the order of basic self-transpose gates. For multi-qubit gates, + * the qubits applied must either be all inputs or all outputs. + */ + void apply_gate( + OpType type, const qubit_vector_t& qbs, + TableauSegment seg = TableauSegment::Output); + + /** + * Initialises a new qubit in the |0> state (or <0| for inputs). O(n) wrt + * number of qubits in the state. + */ + void init_qubit(const Qubit& qb, TableauSegment seg = TableauSegment::Output); + + /** + * Post-selects the chosen qubit to <0| (or applies |0> on inputs), removing + * it from the state. This DOES NOT CHECK whether or not the post-selection + * would be successful (e.g. we can post-select <0| on a |1> state). O(n^3) + * wrt number of qubits in the state. + */ + void post_select( + const Qubit& qb, TableauSegment seg = TableauSegment::Output); + + /** + * Discards the given qubit, removing it from the state. O(n^3) wrt number of + * qubits in the state. + */ + void discard_qubit( + const Qubit& qb, TableauSegment seg = TableauSegment::Output); + + /** + * Collapses the given qubit in the Z basis. O(n^3) wrt number of qubits in + * the state. + */ + void collapse_qubit( + const Qubit& qb, TableauSegment seg = TableauSegment::Output); + + /** + * Permutes columns into a canonical order: chosen segment first, subordered + * in ILO. + */ + void canonical_column_order(TableauSegment first = TableauSegment::Input); + + /** + * Reduces the underlying APState to its normal form wrt the current ordering + * of qubits in its columns. + */ + void normal_form(); + + /** + * Renames qubits. + */ + void rename_qubits(const qubit_map_t& qmap, TableauSegment seg); + + bool operator==(const ChoiAPState& other) const; +}; + +JSON_DECL(APState) +JSON_DECL(ChoiAPState::TableauSegment) +JSON_DECL(ChoiAPState) + +} // namespace tket \ No newline at end of file diff --git a/tket/include/tket/Converters/Converters.hpp b/tket/include/tket/Converters/Converters.hpp index 0a9f9809cc..5061204781 100644 --- a/tket/include/tket/Converters/Converters.hpp +++ b/tket/include/tket/Converters/Converters.hpp @@ -15,7 +15,7 @@ #pragma once #include "tket/Circuit/Circuit.hpp" -#include "tket/Clifford/APState.hpp" +#include "tket/Clifford/APState2.hpp" #include "tket/Clifford/ChoiMixTableau.hpp" #include "tket/Clifford/UnitaryTableau.hpp" #include "tket/PauliGraph/PauliGraph.hpp" @@ -46,20 +46,6 @@ Circuit unitary_rev_tableau_to_circuit(const UnitaryRevTableau &tab); */ ChoiMixTableau circuit_to_cm_tableau(const Circuit &circ); -/** - * Construct an APState for a given circuit. - * Assumes all input qubits are initialised in the |0> state. - * Will throw an exception if it contains non-Clifford gates. - */ -APState circuit_to_apstate(const Circuit &circ); - -/** - * Construct a circuit producing the state described by the APState. - * Uses the standard circuit form implied by the ZX description of reduced - * AP-form (layered as X-H-CX-CZ-S). - */ -Circuit apstate_to_circuit(const APState &ap); - /** * Constructs a circuit producing the same effect as a ChoiMixTableau. * Since Circuit does not support distinct qubit addresses for inputs and @@ -146,10 +132,46 @@ std::pair cm_tableau_to_unitary_extension_circuit( CXConfigType cx_config = CXConfigType::Snake); /** - * Convert a SymplecticTableau describing a stabiliser state to reduced AP-form - * (i.e. an APState object). Since tableau methods do not track global phase, we - * set it to 0. Throws an exception if the tableau does not describe a pure - * state (i.e. it must have n commuting rows for n qubits). + * Construct an APState for a given circuit. + * Assumes all input qubits are initialised in the |0> state. + * Will throw an exception if it contains non-Clifford gates. + */ +APState circuit_to_apstate(const Circuit &circ); + +/** + * Construct a circuit producing the state described by the APState. + * Uses the standard circuit form implied by the ZX description of reduced + * AP-form (layered as X-H-CX-(Collapse-CX)*-CZ-S). + */ +Circuit apstate_to_circuit(const APState &ap); + +/** + * Construct a ChoiAPState for a given circuit. + * Will incorporate qubit initialisations and discarding into the circuit. + * Will throw an exception if it contains non-Clifford gates. + */ +ChoiAPState circuit_to_choi_apstate(const Circuit &circ); + +/** + * Constructs a circuit producing the same effect as a ChoiAPState. + * Uses the same synthesis method as cm_tableau_to_exact_circuit. + */ +std::pair choi_apstate_to_exact_circuit( + ChoiAPState ap, CXConfigType cx_config = CXConfigType::Snake); + +/** + * Constructs a unitary circuit whose stabilizer group contains all the + * stabilizers of the ChoiAPState and possibly more. See + * cm_tableau_to_unitary_extension_circuit for more details. + */ +std::pair choi_apstate_to_unitary_extension_circuit( + ChoiAPState ap, const std::vector &init_names = {}, + const std::vector &post_names = {}, + CXConfigType cx_config = CXConfigType::Snake); + +/** + * Convert a SymplecticTableau describing a stabiliser state to AP-form. Since + * tableau methods do not track global phase, we set it to 0. */ APState tableau_to_apstate(SymplecticTableau tab); @@ -159,7 +181,22 @@ APState tableau_to_apstate(SymplecticTableau tab); * synthesis methods to synthesise APState objects (the result would need * re-simulating as an APState to compare and deduce the required global phase). */ -SymplecticTableau apstate_to_tableau(const APState &ap); +SymplecticTableau apstate_to_tableau(APState ap); + +/** + * Convert a ChoiMixTableau describing a stabiliser state to AP-form. Since + * tableau methods do not track global phase, we set it to 0. + */ +ChoiAPState cm_tableau_to_choi_apstate(const ChoiMixTableau &tab); + +/** + * Convert a ChoiAPState describing a stabiliser state into a tableau form, + * discarding the global phase information. This allows us to reuse tableau + * synthesis methods to synthesise ChoiAPState objects (the result would need + * re-simulating as a ChoiAPState to compare and deduce the required global + * phase). + */ +ChoiMixTableau choi_apstate_to_cm_tableau(const ChoiAPState &ap); PauliGraph circuit_to_pauli_graph(const Circuit &circ); diff --git a/tket/src/Clifford/APState2.cpp b/tket/src/Clifford/APState2.cpp new file mode 100644 index 0000000000..6fae7249ba --- /dev/null +++ b/tket/src/Clifford/APState2.cpp @@ -0,0 +1,1470 @@ +// 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/Clifford/APState2.hpp" + +#include +#include + +#include "tket/OpType/OpTypeInfo.hpp" + +namespace tket { + +/*************************** APState IMPLEMENTATION ***************************/ + +APState::APState( + const MatrixXb& A_, const VectorXb& B_, const MatrixXb& C_, + const MatrixXb& E_, const Eigen::VectorXi& P_, const Expr& phase_) + : A(A_), B(B_), C(C_), E(E_), P(P_), phase(phase_) { + verify(); +} + +APState::APState(unsigned n_qubits) + : A(MatrixXb::Identity(n_qubits, n_qubits)), + B(VectorXb::Zero(n_qubits)), + C(MatrixXb::Zero(n_qubits, n_qubits)), + E(MatrixXb::Zero(n_qubits, n_qubits)), + P(Eigen::VectorXi::Zero(n_qubits)), + phase(0) {} + +unsigned clifford_phase(const Complex& c) { + unsigned res = 0; + unsigned n_possibles = 0; + if (c.real() > EPS) { + res = 0; + ++n_possibles; + } + if (c.imag() > EPS) { + res = 1; + ++n_possibles; + } + if (c.real() < -EPS) { + res = 2; + ++n_possibles; + } + if (c.imag() < -EPS) { + res = 3; + ++n_possibles; + } + TKET_ASSERT(n_possibles == 1); + return res; +} + +APState::APState(const Eigen::VectorXcd& sv) { + unsigned n_qbs = 0; + while (sv.size() > (1 << n_qbs)) ++n_qbs; + TKET_ASSERT(sv.size() == (1 << n_qbs)); + + // Find non-zero entries as a vector space and offset + unsigned z0 = 0; + for (unsigned x = 0; x < sv.size(); ++x) { + if (sv(x) != 0.) { + z0 = x; + break; + } + } + std::vector offsets; + unsigned n_non_zero = 0; + for (unsigned x = 1; x < sv.size(); ++x) { + if (sv(z0 ^ x) != 0.) { + ++n_non_zero; + if (n_non_zero == (1u << offsets.size())) offsets.push_back(x); + } + } + + // Find A as the dual space + MatrixXb offset_mat(offsets.size(), n_qbs); + for (unsigned r = 0; r < offsets.size(); ++r) { + unsigned off = offsets.at(r); + for (unsigned c = 0; c < n_qbs; ++c) { + // Binary encoding of offsets in reverse order to guarantee free qubits + // are the later ones, meaning we produce A in row echelon form + offset_mat(r, c) = (off & (1 << c)) != 0; + } + } + std::vector> row_ops = + gaussian_elimination_row_ops(offset_mat); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + offset_mat(op.second, j) ^= offset_mat(op.first, j); + } + } + std::map mat_leaders; + A = MatrixXb::Zero(n_qbs, n_qbs); + for (unsigned c = 0; c < n_qbs; ++c) { + bool free_qubit = false; + for (unsigned r = 0; r < offsets.size(); ++r) { + if (offset_mat(r, c)) { + auto inserted = mat_leaders.insert({r, c}); + if (inserted.second) { + free_qubit = true; + break; + } else { + // Reverse bit orderings back to normal + A(n_qbs - 1 - c, n_qbs - 1 - inserted.first->second) = true; + } + } + } + A(n_qbs - 1 - c, n_qbs - 1 - c) = !free_qubit; + } + + VectorXb z0_vec(n_qbs); + for (unsigned i = 0; i < n_qbs; ++i) { + z0_vec(i) = ((z0 & (1 << (n_qbs - 1 - i))) != 0); + } + B = A * z0_vec; + + E = MatrixXb::Zero(n_qbs, n_qbs); + P = Eigen::VectorXi::Zero(n_qbs); + unsigned neutral_z0 = z0; // Index with 0 for all free qubits + std::map offset_for_free; + for (const std::pair& row_qfree : mat_leaders) { + unsigned offset = 0; + for (unsigned i = 0; i < n_qbs; ++i) + if (offset_mat(row_qfree.first, i)) offset += (1 << i); + unsigned qfree = n_qbs - 1 - row_qfree.second; + offset_for_free.insert({qfree, offset}); + if ((neutral_z0 & (1 << row_qfree.second)) != 0) + neutral_z0 = neutral_z0 ^ offset; + } + for (const std::pair& qfree_offset : + offset_for_free) { + unsigned qfree = qfree_offset.first; + unsigned offset = qfree_offset.second; + auto complex_phase = sv(neutral_z0 ^ offset) / sv(neutral_z0); + P(qfree) = clifford_phase(complex_phase); + for (const std::pair& qfree_offset2 : + offset_for_free) { + if (qfree_offset == qfree_offset2) break; // Only go up to solved phases + unsigned qfree2 = qfree_offset2.first; + unsigned offset2 = qfree_offset2.second; + auto complex_phase_cz = + sv(neutral_z0 ^ offset ^ offset2) / sv(neutral_z0); + E(qfree, qfree2) = E(qfree2, qfree) = + ((unsigned)(clifford_phase(complex_phase_cz) - P(qfree) - P(qfree2)) % + 4) == 2; + } + } + + phase = Expr(std::arg(sv(neutral_z0)) / PI); +} + +APState::APState(const Eigen::MatrixXcd& dm) { + Eigen::VectorXcd dm_as_vec = dm.reshaped(); + APState pure_double(dm_as_vec); + unsigned n_qbs = pure_double.A.cols() / 2; + A = MatrixXb::Zero(n_qbs, n_qbs); + B = VectorXb::Zero(n_qbs); + C = MatrixXb::Zero(n_qbs, n_qbs); + + MatrixXb full_mat(2 * n_qbs, 2 * n_qbs + 1); + full_mat.block(0, 0, 2 * n_qbs, 2 * n_qbs) = pure_double.A; + full_mat.col(2 * n_qbs) = pure_double.B; + + // This full matrix is generated by AI, IA, and CC. Gaussian elimination puts + // zeros at the bottom, followed by IA. + std::vector> row_ops = + gaussian_elimination_row_ops(full_mat); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < 2 * n_qbs + 1; ++j) { + full_mat(op.second, j) ^= full_mat(op.first, j); + } + } + + unsigned first_zero = 0; + for (unsigned r = n_qbs; r > 0;) { + --r; + bool is_zero = true; + // By symmetry of AI, IA, and CC, the bottom non-empty row must have some + // component on the right side + for (unsigned c = n_qbs; c < 2 * n_qbs; ++c) { + if (full_mat(r, c)) { + is_zero = false; + break; + } + } + if (!is_zero) { + first_zero = r + 1; + break; + } + } + unsigned first_right = 0; + for (unsigned r = first_zero; r > 0;) { + --r; + bool is_right_only = true; + for (unsigned c = 0; c < n_qbs; ++c) { + if (full_mat(r, c)) { + is_right_only = false; + break; + } + } + if (!is_right_only) { + first_right = r + 1; + break; + } + } + A.block(0, 0, first_zero - first_right, n_qbs) = + full_mat.block(first_right, n_qbs, first_zero - first_right, n_qbs); + B.head(first_zero - first_right) = + full_mat.block(first_right, 2 * n_qbs, first_zero - first_right, n_qbs); + + // Flip the column order and reduce the remaining rows to obtain CC above AI; + // get the row combinations from the reordered matrix but apply them to the + // matrix with the correct column ordering + MatrixXb remaining_rows(first_right, 2 * n_qbs); + for (unsigned r = 0; r < first_right; ++r) { + for (unsigned c = 0; c < 2 * n_qbs; ++c) { + remaining_rows(r, c) = full_mat(r, 2 * n_qbs - 1 - c); + } + } + row_ops = gaussian_elimination_row_ops(remaining_rows); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < 2 * n_qbs; ++j) { + full_mat(op.second, j) ^= full_mat(op.first, j); + } + } + + unsigned first_left = 0; + for (unsigned r = first_right; r > 0;) { + --r; + bool is_left_only = true; + for (unsigned c = n_qbs; c < 2 * n_qbs; ++c) { + if (full_mat(r, c)) { + is_left_only = false; + break; + } + } + if (!is_left_only) { + first_left = r + 1; + break; + } + } + C.block(0, 0, first_left, n_qbs) = full_mat.block(0, 0, first_left, n_qbs); + + /** + * Recall that pure_double.A is generated by AI, IA, and CC. The statevector + * constructor builds this in normal form, so it is in reduced row-echelon + * form. The first block of qubits may have leaders for leaders or mixed + * qubits in the density matrix, but the second block only has them for true + * leaders in the density matrix. This means pure_double.E and pure_double.P + * are not just the direct sum of two copies of the true E and P, but the + * bottom segments contain them exactly. + */ + E = pure_double.E.block(n_qbs, n_qbs, n_qbs, n_qbs); + P = pure_double.P.tail(n_qbs); + + phase = 0; +} + +void APState::verify() const { + // Check dimensions all agree + unsigned n_qubits = A.rows(); + TKET_ASSERT(A.cols() == n_qubits); + TKET_ASSERT(B.size() == n_qubits); + TKET_ASSERT(C.rows() == n_qubits); + TKET_ASSERT(C.cols() == n_qubits); + TKET_ASSERT(E.rows() == n_qubits); + TKET_ASSERT(E.cols() == n_qubits); + TKET_ASSERT(P.size() == n_qubits); + + for (unsigned r = 0; r < n_qubits; ++r) { + for (unsigned c = 0; c < r; ++c) { + // Check E is symmetric + TKET_ASSERT(E(r, c) == E(c, r)); + } + // Check diagonals of E are zero + TKET_ASSERT(!E(r, r)); + } +} + +VectorXb z2_mult(const MatrixXb& M, const VectorXb& v) { + VectorXb res = VectorXb::Zero(M.rows()); + for (unsigned i = 0; i < M.cols(); ++i) { + if (v(i)) { + for (unsigned j = 0; j < M.rows(); ++j) res(j) ^= M(j, i); + } + } + return res; +} + +Eigen::VectorXcd APState::to_statevector() const { + unsigned n_qubits = A.cols(); + Eigen::VectorXcd sv = Eigen::VectorXcd::Zero(1 << n_qubits); + unsigned n_terms = 0; + Complex g_phase = std::exp(i_ * PI * eval_expr(phase).value()); + for (unsigned x = 0; x < (1u << n_qubits); ++x) { + VectorXb x_binary = VectorXb::Zero(n_qubits); + for (unsigned i = 0; i < n_qubits; ++i) { + unsigned mask = 1 << (n_qubits - 1 - i); + x_binary(i) = ((x & mask) != 0); + } + + if (z2_mult(A, x_binary) == B) { + ++n_terms; + unsigned i_phases = 0; + for (unsigned q = 0; q < n_qubits; ++q) { + if (x_binary(q)) { + i_phases += P(q); + for (unsigned q2 = q; q2 < n_qubits; ++q2) { + if (E(q, q2) && x_binary(q2)) i_phases += 2; + } + } + } + switch (i_phases % 4) { + case 0: { + sv(x) = g_phase; + break; + } + case 1: { + sv(x) = i_ * g_phase; + break; + } + case 2: { + sv(x) = -g_phase; + break; + } + case 3: { + sv(x) = -i_ * g_phase; + break; + } + default: { + TKET_ASSERT(false); + } + } + } + } + return pow(n_terms, -0.5) * sv; +} + +Eigen::MatrixXcd APState::to_density_matrix() const { + unsigned n_qubits = A.cols(); + Eigen::MatrixXcd dm = Eigen::MatrixXcd::Zero(1 << n_qubits, 1 << n_qubits); + for (unsigned x = 0; x < (1u << n_qubits); ++x) { + VectorXb x_binary = VectorXb::Zero(n_qubits); + for (unsigned i = 0; i < n_qubits; ++i) { + unsigned mask = 1 << (n_qubits - 1 - i); + x_binary(i) = ((x & mask) != 0); + } + + if (z2_mult(A, x_binary) == B) { + MatrixXb Cx = z2_mult(C, x_binary); + for (unsigned y = 0; y < (1u << n_qubits); ++y) { + VectorXb y_binary = VectorXb::Zero(n_qubits); + for (unsigned i = 0; i < n_qubits; ++i) { + unsigned mask = 1 << (n_qubits - 1 - i); + y_binary(i) = ((y & mask) != 0); + } + + if ((z2_mult(A, y_binary) == B) && (z2_mult(C, y_binary) == Cx)) { + unsigned i_phases = 0; + for (unsigned q = 0; q < n_qubits; ++q) { + if (x_binary(q)) { + i_phases += P(q); + for (unsigned q2 = q; q2 < n_qubits; ++q2) { + if (E(q, q2) && x_binary(q2)) i_phases += 2; + } + } + if (y_binary(q)) { + i_phases -= P(q); + for (unsigned q2 = q; q2 < n_qubits; ++q2) { + if (E(q, q2) && y_binary(q2)) i_phases += 2; + } + } + } + switch (i_phases % 4) { + case 0: { + dm(x, y) = 1; + break; + } + case 1: { + dm(x, y) = i_; + break; + } + case 2: { + dm(x, y) = -1; + break; + } + case 3: { + dm(x, y) = -i_; + break; + } + default: { + TKET_ASSERT(false); + } + } + } + } + } + } + return dm / dm.trace(); +} + +void APState::apply_CZ(unsigned ctrl, unsigned trgt) { + E(ctrl, trgt) ^= true; + E(trgt, ctrl) ^= true; +} + +void APState::apply_S(unsigned q) { P(q) += 1; } + +void APState::apply_V(unsigned q) { + unsigned n_qbs = A.cols(); + std::list A_rows; + std::list C_rows; + for (unsigned r = 0; r < n_qbs; ++r) { + if (A(r, q)) A_rows.push_back(r); + if (C(r, q)) C_rows.push_back(r); + } + if ((unsigned)P(q) % 2 == 0) { + bool a = ((unsigned)P(q) % 4) == 2; + if (A_rows.empty()) { + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (E(q, q2)) { + // Update local phase on neighbours + P(q2) += a ? 3 : 1; + // Local complementation between neighbours + for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) { + if (E(q, q3)) { + E(q2, q3) ^= true; + E(q3, q2) ^= true; + } + } + // Add connections to all C_rows + for (unsigned r : C_rows) { + C(r, q2) ^= true; + } + } + } + // Global phase + phase += a ? .25 : -.25; + } else { + // Choose one of the red spiders to remove + unsigned r = A_rows.back(); + A_rows.pop_back(); + // Stratify neighbourhoods of r (in A) and q (in E) + std::list just_r; + std::list just_q; + std::list both; + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (q2 == q) continue; + if (A(r, q2)) { + if (E(q, q2)) + both.push_back(q2); + else + just_r.push_back(q2); + } else if (E(q, q2)) + just_q.push_back(q2); + } + // Update A and B + for (unsigned ar : A_rows) { + A(ar, q) = false; + for (unsigned q2 : just_r) A(ar, q2) ^= true; + for (unsigned q2 : both) A(ar, q2) ^= true; + B(ar) ^= B(r); + } + // Update C + for (unsigned cr : C_rows) { + C(cr, q) = false; + for (unsigned q2 : just_r) C(cr, q2) ^= true; + for (unsigned q2 : both) C(cr, q2) ^= true; + } + // Update E and P + for (unsigned rn : just_r) { + // Complementation within just_r + for (unsigned rn2 : just_r) { + E(rn, rn2) ^= true; + } + // Reset diagonal + E(rn, rn) = false; + // Complementation between just_r and just_q + for (unsigned qn : just_q) { + E(rn, qn) ^= true; + E(qn, rn) ^= true; + } + // Connect to q + E(rn, q) = true; + E(q, rn) = true; + // Local phases + P(rn) += (a ^ B(r)) ? 1 : 3; + } + for (unsigned bn : both) { + // Complementation within both + for (unsigned bn2 : both) { + E(bn, bn2) ^= true; + } + // Reset diagonal + E(bn, bn) = false; + // Complementation between both and just_q + for (unsigned qn : just_q) { + E(bn, qn) ^= true; + E(qn, bn) ^= true; + } + // Local phases + P(bn) += a ? 3 : 1; + } + for (unsigned qn : just_q) { + // Remove connection to q + E(q, qn) = false; + E(qn, q) = false; + // Local phases + P(qn) += B(r) ? 2 : 0; + } + // Local phase on q + P(q) = B(r) ? 1 : 3; + // Global phase + if (B(r)) phase += a ? .5 : 1.5; + // Remove row r from A and B + for (unsigned i = 0; i < n_qbs; ++i) { + A(r, i) = false; + } + B(r) = false; + } + } else { + bool a = ((unsigned)P(q) % 4) == 3; + if (!A_rows.empty()) { + // Choose one of the red spiders to remove + unsigned r = A_rows.back(); + A_rows.pop_back(); + // Stratify neighbourhoods of r (in A) and q (in E) + std::list just_r; + std::list just_q; + std::list both; + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (q2 == q) continue; + if (A(r, q2)) { + if (E(q, q2)) + both.push_back(q2); + else + just_r.push_back(q2); + } else if (E(q, q2)) + just_q.push_back(q2); + } + // Update A and B + for (unsigned ar : A_rows) { + A(ar, q) = false; + for (unsigned q2 : just_r) A(ar, q2) ^= true; + for (unsigned q2 : both) A(ar, q2) ^= true; + B(ar) ^= B(r); + } + // Update C + for (unsigned cr : C_rows) { + C(cr, q) = false; + for (unsigned q2 : just_r) C(cr, q2) ^= true; + for (unsigned q2 : both) C(cr, q2) ^= true; + } + // Update E and P + for (unsigned rn : just_r) { + // Complementation between just_r and just_q + for (unsigned qn : just_q) { + E(rn, qn) ^= true; + E(qn, rn) ^= true; + } + // Complementation between just_r and both + for (unsigned bn : both) { + E(rn, bn) ^= true; + E(bn, rn) ^= true; + } + // Connect to q + E(rn, q) = true; + E(q, rn) = true; + // Local phases + P(rn) += a ? 2 : 0; + } + for (unsigned bn : both) { + // Complementation between both and just_q + for (unsigned qn : just_q) { + E(bn, qn) ^= true; + E(qn, bn) ^= true; + } + // Local phases + P(bn) += (a ^ B(r)) ? 0 : 2; + } + for (unsigned qn : just_q) { + // Remove connection to q + E(qn, q) = false; + E(q, qn) = false; + // Local phases + P(qn) += B(r) ? 2 : 0; + } + // Local phase on q + P(q) = B(r) ? 1 : 3; + // Global phase + if (a && B(r)) phase += 1; + // Remove row r from A and B + for (unsigned i = 0; i < n_qbs; ++i) { + A(r, i) = false; + } + B(r) = false; + } else if (!C_rows.empty()) { + // Choose one of the discards to pivot around + unsigned d = C_rows.back(); + C_rows.pop_back(); + // Stratify neighbourhoods of d (in C) and q (in E) + std::list just_d; + std::list just_q; + std::list both; + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (q2 == q) continue; + if (C(d, q2)) { + if (E(q, q2)) + both.push_back(q2); + else + just_d.push_back(q2); + } else if (E(q, q2)) + just_q.push_back(q2); + } + // Update C + for (unsigned cr : C_rows) { + C(cr, q) = false; + for (unsigned q2 : just_d) C(cr, q2) ^= true; + for (unsigned q2 : both) C(cr, q2) ^= true; + } + // Update E and P + for (unsigned dn : just_d) { + // Complementation between just_d and just_q + for (unsigned qn : just_q) { + E(dn, qn) ^= true; + E(qn, dn) ^= true; + } + // Complementation between just_d and both + for (unsigned bn : both) { + E(dn, bn) ^= true; + E(bn, dn) ^= true; + } + // Connect to q + E(dn, q) = true; + E(q, dn) = true; + // Remove connection with d + C(d, dn) = false; + // Local phases + P(dn) += a ? 2 : 0; + } + for (unsigned bn : both) { + // Complementation between both and just_q + for (unsigned qn : just_q) { + E(bn, qn) ^= true; + E(qn, bn) ^= true; + } + // Local phases + P(bn) += a ? 0 : 2; + } + for (unsigned qn : just_q) { + // Connect to d + C(d, qn) = true; + // Remove connection to q + E(qn, q) = false; + E(q, qn) = false; + // No local phase change + } + // Local phase on q + P(q) += 3; + // No global phase change + } else { + VectorXb new_row = VectorXb::Zero(n_qbs); + new_row(q) = true; + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (E(q, q2)) { + // Connect to new red spider + new_row(q2) = true; + // Local phase on neighbours + P(q2) += a ? 1 : 3; + // Local complementation between neighbours + for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) { + if (E(q, q3)) { + E(q2, q3) ^= true; + E(q3, q2) ^= true; + } + } + // Remove connection with q + E(q, q2) = false; + E(q2, q) = false; + } + } + // Reset local phase on q + P(q) = 0; + // Global phase + phase += a ? 1.5 : 0; + // Add new_row to A and B + MatrixXb combined_mat = MatrixXb::Zero(n_qbs + 1, n_qbs + 1); + combined_mat.block(0, 0, n_qbs, n_qbs) = A; + combined_mat.block(0, n_qbs, n_qbs, 1) = B; + combined_mat.block(n_qbs, 0, 1, n_qbs) = new_row.transpose(); + combined_mat(n_qbs, n_qbs) = a; + std::vector> row_ops = + gaussian_elimination_row_ops(combined_mat); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs + 1; ++j) { + combined_mat(op.second, j) ^= combined_mat(op.first, j); + } + } + A = combined_mat.block(0, 0, n_qbs, n_qbs); + B = combined_mat.block(0, n_qbs, n_qbs, 1); + } + } +} + +void APState::apply_X(unsigned q) { + // Push through the local phase + phase += P(q) * .5; + P(q) = -P(q); + // Pushing through the CZs adds Zs on neighbours + for (unsigned q2 = 0; q2 < E.cols(); ++q2) { + if (E(q, q2)) P(q2) += 2; + } + // Pushing through adjacency matrix adds Xs onto connected reds + for (unsigned r = 0; r < A.rows(); ++r) { + if (A(r, q)) B(r) ^= true; + } +} + +void APState::apply_gate(OpType type, const std::vector& qbs) { + switch (type) { + case OpType::Z: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::X: { + apply_X(qbs.at(0)); + break; + } + case OpType::Y: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + apply_X(qbs.at(0)); + phase += .5; + break; + } + case OpType::S: { + apply_S(qbs.at(0)); + break; + } + case OpType::Sdg: { + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::V: { + apply_V(qbs.at(0)); + break; + } + case OpType::SX: { + apply_V(qbs.at(0)); + phase += .25; + break; + } + case OpType::Vdg: { + apply_V(qbs.at(0)); + apply_X(qbs.at(0)); + phase += .5; + break; + } + case OpType::SXdg: { + apply_V(qbs.at(0)); + apply_X(qbs.at(0)); + phase += .25; + break; + } + case OpType::H: { + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::CX: { + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + break; + } + case OpType::CY: { + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(1)); + apply_X(qbs.at(1)); + phase += .5; + break; + } + case OpType::CZ: { + apply_CZ(qbs.at(0), qbs.at(1)); + break; + } + case OpType::ZZMax: { + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + phase -= .25; + break; + } + case OpType::ECR: { + apply_S(qbs.at(0)); + apply_X(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + phase += -.25; + break; + } + case OpType::ISWAPMax: { + apply_V(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(0)); + apply_V(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_V(qbs.at(0)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + phase += 1.; + break; + } + case OpType::SWAP: { + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + apply_S(qbs.at(1)); + apply_V(qbs.at(1)); + apply_S(qbs.at(1)); + apply_CZ(qbs.at(0), qbs.at(1)); + apply_S(qbs.at(0)); + apply_V(qbs.at(0)); + apply_S(qbs.at(0)); + break; + } + case OpType::BRIDGE: { + apply_S(qbs.at(2)); + apply_V(qbs.at(2)); + apply_S(qbs.at(2)); + apply_CZ(qbs.at(0), qbs.at(2)); + apply_S(qbs.at(2)); + apply_V(qbs.at(2)); + apply_S(qbs.at(2)); + break; + } + case OpType::noop: { + break; + } + case OpType::Reset: { + unsigned q = qbs.at(0); + discard_qubit(q); + unsigned new_qb = init_qubit(); + // Swap qubits to preserve the ordering + if (q != new_qb) { + unsigned n_qbs = A.cols(); + for (unsigned r = 0; r < n_qbs; ++r) { + bool temp = A(r, q); + A(r, q) = A(r, new_qb); + A(r, new_qb) = temp; + } + for (unsigned r = 0; r < n_qbs; ++r) { + bool temp = C(r, q); + C(r, q) = C(r, new_qb); + C(r, new_qb) = temp; + } + for (unsigned r = 0; r < n_qbs; ++r) { + bool temp = E(r, q); + E(r, q) = E(r, new_qb); + E(r, new_qb) = temp; + } + for (unsigned c = 0; c < n_qbs; ++c) { + bool temp = E(q, c); + E(q, c) = E(new_qb, c); + E(new_qb, c) = temp; + } + int temp = P(q); + P(q) = P(new_qb); + P(new_qb) = temp; + } + break; + } + case OpType::Collapse: { + collapse_qubit(qbs.at(0)); + break; + } + case OpType::Phase: { + throw std::logic_error( + "OpType::Phase cannot be applied via APState::apply_gate"); + } + default: { + throw BadOpType( + "Cannot be applied to a APState: not a Clifford gate", type); + } + } +} + +unsigned APState::init_qubit() { + unsigned n_qbs = A.cols(); + A.conservativeResize(n_qbs + 1, n_qbs + 1); + A.col(n_qbs) = VectorXb::Zero(n_qbs + 1); + A.row(n_qbs) = VectorXb::Zero(n_qbs + 1); + A(n_qbs, n_qbs) = true; + B.conservativeResize(n_qbs + 1); + B(n_qbs) = false; + C.conservativeResize(n_qbs + 1, n_qbs + 1); + C.col(n_qbs) = VectorXb::Zero(n_qbs + 1); + C.row(n_qbs) = VectorXb::Zero(n_qbs + 1); + E.conservativeResize(n_qbs + 1, n_qbs + 1); + E.col(n_qbs) = VectorXb::Zero(n_qbs + 1); + E.row(n_qbs) = VectorXb::Zero(n_qbs + 1); + P.conservativeResize(n_qbs + 1); + P(n_qbs) = 0; + return n_qbs; +} + +unsigned APState::post_select(unsigned q) { + unsigned n_qbs = A.cols(); + if (q >= n_qbs) + throw std::invalid_argument("APState: post-selecting a non-existent qubit"); + MatrixXb AB = MatrixXb::Zero(n_qbs, n_qbs); + AB.block(0, 0, n_qbs, n_qbs) = A; + AB.col(q) = A.col(n_qbs - 1); + AB.col(n_qbs - 1) = B; + std::vector> row_ops = + gaussian_elimination_row_ops(AB); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + AB(op.second, j) ^= AB(op.first, j); + } + } + A.resize(n_qbs - 1, n_qbs - 1); + A = AB.block(0, 0, n_qbs - 1, n_qbs - 1); + B.resize(n_qbs - 1); + B = AB.block(0, n_qbs - 1, n_qbs - 1, 1); + C.col(q) = C.col(n_qbs - 1); + row_ops = gaussian_elimination_row_ops(C); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + C(op.second, j) ^= C(op.first, j); + } + } + C.conservativeResize(n_qbs - 1, n_qbs - 1); + E.col(q) = E.col(n_qbs - 1); + E.row(q) = E.row(n_qbs - 1); + E.conservativeResize(n_qbs - 1, n_qbs - 1); + P(q) = P(n_qbs - 1); + P.conservativeResize(n_qbs - 1); + return n_qbs - 1; +} + +void APState::collapse_qubit(unsigned q) { + unsigned n_qbs = A.cols(); + MatrixXb C_ext = MatrixXb::Zero(n_qbs + 1, n_qbs); + C_ext.block(0, 0, n_qbs, n_qbs) = C; + C_ext(n_qbs, q) = 1; + std::vector> row_ops = + gaussian_elimination_row_ops(C_ext); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + C_ext(op.second, j) ^= C_ext(op.first, j); + } + } + C = C_ext.block(0, 0, n_qbs, n_qbs); +} + +unsigned APState::discard_qubit(unsigned q) { + collapse_qubit(q); + apply_V(q); + collapse_qubit(q); + return post_select(q); +} + +void APState::normal_form() { + unsigned n_qbs = A.cols(); + // Get A into reduced row-echelon form + std::vector> row_ops = + gaussian_elimination_row_ops(A); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + A(op.second, j) ^= A(op.first, j); + } + B(op.second) ^= B(op.first); + } + // Identify leading qubits + std::map leader_to_row; + for (unsigned r = 0; r < n_qbs; ++r) { + bool leader_found = false; + for (unsigned c = 0; c < n_qbs; ++c) { + if (A(r, c)) { + leader_found = true; + leader_to_row.insert({c, r}); + break; + } + } + if (!leader_found) break; + } + // Remove leaders from C + for (unsigned r = 0; r < n_qbs; ++r) { + for (const std::pair& lr : leader_to_row) { + if (C(r, lr.first)) { + for (unsigned c = 0; c < n_qbs; ++c) { + C(r, c) ^= A(lr.second, c); + } + } + } + } + // Remove leaders from E + for (const std::pair& lr : leader_to_row) { + for (unsigned q2 = lr.first; q2 < n_qbs; ++q2) { + if (E(lr.first, q2)) { + E(lr.first, q2) = false; + E(q2, lr.first) = false; + for (unsigned q3 = lr.first + 1; q3 < n_qbs; ++q3) { + if (A(lr.second, q3)) { + E(q2, q3) ^= true; + E(q3, q2) ^= true; + } + } + P(q2) += (B(lr.second) ^ A(lr.second, q2)) ? 2 : 0; + } + } + } + // Remove leaders from P + for (const std::pair& lr : leader_to_row) { + if ((unsigned)P(lr.first) % 2 == 1) { + // Complementation within neighbours under A (excluding the leader) + for (unsigned q2 = lr.first + 1; q2 < n_qbs; ++q2) { + if (A(lr.second, q2)) { + for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) { + if (A(lr.second, q3)) { + E(q2, q3) ^= true; + E(q3, q2) ^= true; + } + } + // Local phases of neighbours + P(q2) += (P(lr.first) % 4) + (B(lr.second) ? 2 : 0); + } + } + } else if ((unsigned)P(lr.first) % 4 == 2) { + // Local phases of neighbours + for (unsigned q2 = lr.first + 1; q2 < n_qbs; ++q2) { + if (A(lr.second, q2)) P(q2) += 2; + } + } + // Global phase + if (B(lr.second)) phase += P(lr.first) * .5; + // Reset P + P(lr.first) = 0; + } + // Get C into reduced row-echelon form + row_ops = gaussian_elimination_row_ops(C); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + C(op.second, j) ^= C(op.first, j); + } + } + // Identify mixed qubits + std::map mixed_to_row; + for (unsigned r = 0; r < n_qbs; ++r) { + bool mixed_found = false; + for (unsigned c = 0; c < n_qbs; ++c) { + if (C(r, c)) { + mixed_found = true; + mixed_to_row.insert({c, r}); + break; + } + } + if (!mixed_found) break; + } + // Remove E connections between mixed qubits + for (auto mr1 = mixed_to_row.begin(); mr1 != mixed_to_row.end(); ++mr1) { + for (auto mr2 = mixed_to_row.begin(); mr2 != mr1; ++mr2) { + if (E(mr1->first, mr2->first)) { + // Local complementation along the sum of their rows in C + // mr2->first < mr1->first, and they are the first entries in their rows + for (unsigned i = mr2->first; i < n_qbs; ++i) { + if (C(mr1->second, i) ^ C(mr2->second, i)) { + for (unsigned j = i + 1; j < n_qbs; ++j) { + if (C(mr1->second, j) ^ C(mr2->second, j)) { + E(i, j) ^= true; + E(j, i) ^= true; + } + } + // Add local phases around neighbourhood + P(i) += 1; + } + } + } + } + } + // Remove mixed qubits from P + for (const std::pair& mr : mixed_to_row) { + unsigned Pm = (unsigned)P(mr.first) % 4; + if (Pm % 2 == 1) { + // Complementation within row of C (including the mixed qubit) + for (unsigned q2 = mr.first; q2 < n_qbs; ++q2) { + if (C(mr.second, q2)) { + for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) { + if (C(mr.second, q3)) { + E(q2, q3) ^= true; + E(q3, q2) ^= true; + } + } + // Local phases (we saved to Pm earlier so we can reset P(mr.first) + // here) + P(q2) += -Pm; + } + } + } else if (Pm == 2) { + // Local phases within row of C + for (unsigned q2 = mr.first; q2 < n_qbs; ++q2) { + if (C(mr.second, q2)) { + P(q2) += 2; + } + } + } + // No global phase change + } +} + +bool APState::operator==(const APState& other) const { + for (unsigned i = 0; i < P.size(); ++i) { + if (((unsigned)P(i) % 4) != ((unsigned)other.P(i) % 4)) return false; + } + return (A == other.A) && (B == other.B) && (C == other.C) && (E == other.E) && + (phase == other.phase); +} + +void to_json(nlohmann::json& j, const APState& aps) { + j["nqubits"] = aps.A.cols(); + j["A"] = aps.A; + j["B"] = aps.B; + j["C"] = aps.C; + j["E"] = aps.E; + j["P"] = aps.P; + j["phase"] = aps.phase; +} + +void from_json(const nlohmann::json& j, APState& aps) { + unsigned n_qbs = j.at("nqubits").get(); + MatrixXb A(n_qbs, n_qbs); + VectorXb B(n_qbs); + MatrixXb C(n_qbs, n_qbs); + MatrixXb E(n_qbs, n_qbs); + Eigen::VectorXi P(n_qbs); + from_json(j.at("A"), A); + from_json(j.at("B"), B); + from_json(j.at("C"), C); + from_json(j.at("E"), E); + from_json(j.at("P"), P); + Expr phase = j.at("phase").get(); + aps = APState(A, B, C, E, P, phase); +} + +/************************* ChoiAPState IMPLEMENTATION *************************/ + +static APState id_aps(unsigned n) { + MatrixXb A(2 * n, 2 * n); + A << MatrixXb::Identity(n, n), MatrixXb::Identity(n, n), + MatrixXb::Zero(n, 2 * n); + return APState( + A, VectorXb::Zero(2 * n), MatrixXb::Zero(2 * n, 2 * n), + MatrixXb::Zero(2 * n, 2 * n), Eigen::VectorXi::Zero(2 * n), 0); +} + +ChoiAPState::ChoiAPState(unsigned n) : ap_(id_aps(n)), col_index_() { + for (unsigned i = 0; i < n; ++i) { + col_index_.insert({{Qubit(i), TableauSegment::Input}, i}); + col_index_.insert({{Qubit(i), TableauSegment::Output}, n + i}); + } +} + +ChoiAPState::ChoiAPState(const qubit_vector_t& qbs) + : ap_(id_aps(qbs.size())), col_index_() { + unsigned n = qbs.size(); + unsigned i = 0; + for (const Qubit& qb : qbs) { + col_index_.insert({{qb, TableauSegment::Input}, i}); + col_index_.insert({{qb, TableauSegment::Output}, n + i}); + ++i; + } +} + +ChoiAPState::ChoiAPState( + const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E, + const Eigen::VectorXi& P, const Expr& phase, unsigned n_ins) + : ap_(A, B, C, E, P, phase), col_index_() { + unsigned n_qbs = A.cols(); + if (n_ins > n_qbs) + throw std::invalid_argument( + "Number of inputs of a ChoiAPState cannot be larger than the number of " + "qubits"); + for (unsigned i = 0; i < n_ins; ++i) { + col_index_.insert({{Qubit(i), TableauSegment::Input}, i}); + } + for (unsigned i = 0; i < n_qbs - n_ins; ++i) { + col_index_.insert({{Qubit(i), TableauSegment::Output}, n_ins + i}); + } +} + +unsigned ChoiAPState::get_n_boundaries() const { return ap_.A.cols(); } + +unsigned ChoiAPState::get_n_inputs() const { return input_qubits().size(); } + +unsigned ChoiAPState::get_n_outputs() const { return output_qubits().size(); } + +qubit_vector_t ChoiAPState::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) + ins.push_back(entry.first.first); + } + return ins; +} + +qubit_vector_t ChoiAPState::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) + outs.push_back(entry.first.first); + } + return outs; +} + +void ChoiAPState::apply_gate( + OpType type, const qubit_vector_t& qbs, TableauSegment seg) { + std::vector u_qbs; + for (const Qubit& q : qbs) { + u_qbs.push_back(col_index_.left.at(col_key_t{q, seg})); + } + if (seg == TableauSegment::Output) { + ap_.apply_gate(type, u_qbs); + } + if (seg == TableauSegment::Input) { + switch (type) { + case OpType::Z: + case OpType::X: + case OpType::S: + case OpType::Sdg: + case OpType::SX: + case OpType::V: + case OpType::SXdg: + case OpType::Vdg: + case OpType::H: + case OpType::CX: + case OpType::CZ: + case OpType::ZZMax: + case OpType::ISWAPMax: + case OpType::SWAP: + case OpType::BRIDGE: + case OpType::noop: { + ap_.apply_gate(type, u_qbs); + break; + } + case OpType::Y: { + ap_.apply_gate(OpType::Y, u_qbs); + ap_.phase += 1.; + break; + } + case OpType::CY: { + ap_.apply_V(u_qbs.at(1)); + ap_.apply_X(u_qbs.at(1)); + ap_.apply_CZ(u_qbs.at(0), u_qbs.at(1)); + ap_.apply_V(u_qbs.at(1)); + break; + } + case OpType::ECR: { + ap_.apply_S(u_qbs.at(1)); + ap_.apply_V(u_qbs.at(1)); + ap_.apply_S(u_qbs.at(1)); + ap_.apply_CZ(u_qbs.at(0), u_qbs.at(1)); + ap_.apply_V(u_qbs.at(1)); + ap_.apply_S(u_qbs.at(1)); + ap_.apply_X(u_qbs.at(0)); + ap_.apply_S(u_qbs.at(0)); + ap_.phase += .25; + break; + } + case OpType::Reset: { + unsigned q = u_qbs.at(0); + unsigned removed_q = ap_.post_select(q); + tableau_col_index_t::right_iterator it = + col_index_.right.find(removed_q); + col_key_t removed_key = it->second; + col_index_.right.erase(it); + it = col_index_.right.find(u_qbs.at(0)); + col_index_.right.erase(it); + col_index_.insert({removed_key, q}); + unsigned new_q = ap_.init_qubit(); + ap_.apply_V(new_q); + ap_.collapse_qubit(new_q); + col_index_.insert({{qbs.at(0), seg}, new_q}); + break; + } + case OpType::Phase: { + throw std::logic_error( + "OpType::Phase cannot be applied via ChoiAPState::apply_gate"); + } + default: { + throw BadOpType( + "Cannot be applied to a ChoiAPState: not a Clifford gate", type); + } + } + } +} + +void ChoiAPState::init_qubit(const Qubit& qb, TableauSegment seg) { + unsigned u_qb = ap_.init_qubit(); + col_index_.insert({{qb, seg}, u_qb}); +} + +void ChoiAPState::post_select(const Qubit& qb, TableauSegment seg) { + tableau_col_index_t::left_iterator lit = + col_index_.left.find(col_key_t{qb, seg}); + unsigned u_qb = lit->second; + unsigned removed_u_qb = ap_.post_select(u_qb); + col_index_.left.erase(lit); + tableau_col_index_t::right_iterator rit = col_index_.right.find(removed_u_qb); + col_key_t removed_key = rit->second; + col_index_.right.erase(rit); + col_index_.insert({removed_key, u_qb}); +} + +void ChoiAPState::discard_qubit(const Qubit& qb, TableauSegment seg) { + collapse_qubit(qb, seg); + apply_gate(OpType::V, {qb}, seg); + collapse_qubit(qb, seg); + post_select(qb, seg); +} + +void ChoiAPState::collapse_qubit(const Qubit& qb, TableauSegment seg) { + unsigned u_qb = col_index_.left.at(col_key_t{qb, seg}); + ap_.collapse_qubit(u_qb); +} + +void ChoiAPState::canonical_column_order(TableauSegment first) { + std::set ins; + std::set outs; + BOOST_FOREACH ( + tableau_col_index_t::left_const_reference entry, col_index_.left) { + if (entry.first.second == TableauSegment::Input) + ins.insert(entry.first.first); + else + outs.insert(entry.first.first); + } + tableau_col_index_t new_index; + unsigned i = 0; + if (first == TableauSegment::Input) { + for (const Qubit& q : ins) { + new_index.insert({{q, TableauSegment::Input}, i}); + ++i; + } + } + for (const Qubit& q : outs) { + new_index.insert({{q, TableauSegment::Output}, i}); + ++i; + } + if (first == TableauSegment::Output) { + for (const Qubit& q : ins) { + new_index.insert({{q, TableauSegment::Input}, i}); + ++i; + } + } + MatrixXb A = MatrixXb::Zero(i, i); + MatrixXb C = MatrixXb::Zero(i, i); + // E requires reordering both columns and rows + // Reorder columns to get Etemp, then reorder rows for E + MatrixXb Etemp = MatrixXb::Zero(i, i); + MatrixXb E = MatrixXb::Zero(i, i); + Eigen::VectorXi P = Eigen::VectorXi(i); + for (unsigned j = 0; j < i; ++j) { + col_key_t key = new_index.right.at(j); + unsigned c = col_index_.left.at(key); + A.col(j) = ap_.A.col(c); + C.col(j) = ap_.C.col(c); + Etemp.col(j) = ap_.E.col(c); + P(j) = ap_.P(c); + } + for (unsigned j = 0; j < i; ++j) { + col_key_t key = new_index.right.at(j); + unsigned c = col_index_.left.at(key); + E.row(j) = Etemp.row(c); + } + ap_ = APState(A, ap_.B, C, E, P, ap_.phase); + col_index_ = new_index; +} + +void ChoiAPState::normal_form() { ap_.normal_form(); } + +void ChoiAPState::rename_qubits(const qubit_map_t& qmap, TableauSegment seg) { + tableau_col_index_t new_index; + BOOST_FOREACH ( + tableau_col_index_t::left_const_reference entry, col_index_.left) { + auto found = qmap.find(entry.first.first); + if (entry.first.second == seg && found != qmap.end()) + new_index.insert({{found->second, seg}, entry.second}); + else + new_index.insert({entry.first, entry.second}); + } + col_index_ = new_index; +} + +bool ChoiAPState::operator==(const ChoiAPState& other) const { + return (col_index_ == other.col_index_) && (ap_ == other.ap_); +} + +void to_json(nlohmann::json& j, const ChoiAPState::TableauSegment& seg) { + j = (seg == ChoiAPState::TableauSegment::Input) ? "In" : "Out"; +} + +void from_json(const nlohmann::json& j, ChoiAPState::TableauSegment& seg) { + const std::string str_seg = j.get(); + seg = (str_seg == "In") ? ChoiAPState::TableauSegment::Input + : ChoiAPState::TableauSegment::Output; +} + +void to_json(nlohmann::json& j, const ChoiAPState& cap) { + j["aps"] = cap.ap_; + std::vector qbs; + for (unsigned i = 0; i < cap.get_n_boundaries(); ++i) { + qbs.push_back(cap.col_index_.right.at(i)); + } + j["qubits"] = qbs; +} + +void from_json(const nlohmann::json& j, ChoiAPState& cap) { + j.at("aps").get_to(cap.ap_); + std::vector qbs = + j.at("qubits").get>(); + if ((unsigned)qbs.size() != (unsigned)cap.ap_.A.cols()) + throw std::invalid_argument( + "Number of qubits in json ChoiAPState does not match APState " + "size."); + cap.col_index_.clear(); + for (unsigned i = 0; i < qbs.size(); ++i) { + cap.col_index_.insert({qbs.at(i), i}); + } +} + +} // namespace tket \ No newline at end of file diff --git a/tket/src/Clifford/ChoiMixTableau.cpp b/tket/src/Clifford/ChoiMixTableau.cpp index ef0b08410e..9a1b1b6271 100644 --- a/tket/src/Clifford/ChoiMixTableau.cpp +++ b/tket/src/Clifford/ChoiMixTableau.cpp @@ -71,7 +71,7 @@ ChoiMixTableau::ChoiMixTableau( col_index_.insert({{Qubit(i), TableauSegment::Input}, i}); } for (unsigned i = 0; i < n_bounds - n_ins; ++i) { - col_index_.insert({{Qubit(i), TableauSegment::Output}, i}); + col_index_.insert({{Qubit(i), TableauSegment::Output}, n_ins + i}); } } diff --git a/tket/src/Converters/APStateConverters2.cpp b/tket/src/Converters/APStateConverters2.cpp new file mode 100644 index 0000000000..54cc01c2d7 --- /dev/null +++ b/tket/src/Converters/APStateConverters2.cpp @@ -0,0 +1,433 @@ +// 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/Converters.hpp" + +namespace tket { + +APState circuit_to_apstate(const Circuit& circ) { + APState aps(circ.n_qubits()); + std::map qb_ordering; + for (const Qubit& q : circ.all_qubits()) + qb_ordering.insert({q, (unsigned)qb_ordering.size()}); + for (const Command& com : circ) { + auto args = com.get_args(); + std::vector qbs; + for (const UnitID& q : args) qbs.push_back(qb_ordering.at(q)); + aps.apply_gate(com.get_op_ptr()->get_type(), qbs); + } + return aps; +} + +Circuit apstate_to_circuit(const APState& ap) { + unsigned n_qbs = ap.A.cols(); + Circuit circ(n_qbs); + circ.qubit_create_all(); + + MatrixXb ABgauss(n_qbs, n_qbs + 1); + ABgauss.block(0, 0, n_qbs, n_qbs) = ap.A; + ABgauss.col(n_qbs) = ap.B; + std::vector> row_ops = + gaussian_elimination_row_ops(ABgauss); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs + 1; ++j) { + ABgauss(op.second, j) ^= ABgauss(op.first, j); + } + } + std::map leader_to_row; + for (unsigned r = 0; r < n_qbs; ++r) { + bool empty_row = true; + for (unsigned c = 0; c < n_qbs; ++c) { + if (ABgauss(r, c)) { + leader_to_row.insert({c, r}); + empty_row = false; + break; + } + } + if (empty_row) break; + } + + for (unsigned q = 0; q < n_qbs; ++q) { + auto lr = leader_to_row.find(q); + if (lr == leader_to_row.end()) + circ.add_op(OpType::H, {q}); + else if (ABgauss(lr->second, n_qbs)) + circ.add_op(OpType::X, {q}); + } + for (const std::pair& lr : leader_to_row) { + for (unsigned ctrl = lr.first + 1; ctrl < n_qbs; ++ctrl) { + if (ABgauss(lr.second, ctrl)) + circ.add_op(OpType::CX, {ctrl, lr.first}); + } + } + for (unsigned d = 0; d < n_qbs; ++d) { + std::optional first = std::nullopt; + for (unsigned c = 0; c < n_qbs; ++c) { + if (ap.C(d, c)) { + if (first) { + circ.add_op(OpType::CX, {c, *first}); + } else { + first = c; + } + } + } + if (first) { + circ.add_op(OpType::Collapse, {*first}); + for (unsigned c = n_qbs; c > *first + 1;) { + --c; + if (ap.C(d, c)) circ.add_op(OpType::CX, {c, *first}); + } + } + } + for (unsigned q1 = 0; q1 < n_qbs; ++q1) { + for (unsigned q2 = q1 + 1; q2 < n_qbs; ++q2) { + if (ap.E(q1, q2)) circ.add_op(OpType::CZ, {q1, q2}); + } + switch ((unsigned)ap.P(q1) % 4) { + case 1: { + circ.add_op(OpType::S, {q1}); + break; + } + case 2: { + circ.add_op(OpType::Z, {q1}); + break; + } + case 3: { + circ.add_op(OpType::Sdg, {q1}); + break; + } + default: { + break; + } + } + } + circ.add_phase(ap.phase); + return circ; +} + +ChoiAPState circuit_to_choi_apstate(const Circuit& circ) { + ChoiAPState ap(circ.all_qubits()); + for (const Qubit& q : circ.created_qubits()) { + ap.post_select(q, ChoiAPState::TableauSegment::Input); + } + for (const Command& com : circ) { + auto args = com.get_args(); + qubit_vector_t qbs = {args.begin(), args.end()}; + ap.apply_gate(com.get_op_ptr()->get_type(), qbs); + } + ap.rename_qubits( + circ.implicit_qubit_permutation(), ChoiAPState::TableauSegment::Output); + for (const Qubit& q : circ.discarded_qubits()) { + ap.discard_qubit(q); + } + ap.canonical_column_order(); + ap.normal_form(); + return ap; +} + +std::pair choi_apstate_to_exact_circuit( + ChoiAPState ap, CXConfigType cx_config) { + ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap); + std::pair res = + cm_tableau_to_exact_circuit(tab, cx_config); + ChoiAPState reconstructed = circuit_to_choi_apstate(res.first); + reconstructed.rename_qubits(res.second, ChoiAPState::TableauSegment::Output); + ap.canonical_column_order(); + ap.normal_form(); + reconstructed.canonical_column_order(); + reconstructed.normal_form(); + res.first.add_phase(ap.ap_.phase - reconstructed.ap_.phase); + return res; +} + +std::pair choi_apstate_to_unitary_extension_circuit( + ChoiAPState ap, const std::vector& init_names, + const std::vector& post_names, CXConfigType cx_config) { + ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap); + std::pair res = cm_tableau_to_unitary_extension_circuit( + tab, init_names, post_names, cx_config); + ChoiAPState reconstructed = circuit_to_choi_apstate(res.first); + reconstructed.rename_qubits(res.second, ChoiAPState::TableauSegment::Output); + ap.canonical_column_order(); + ap.normal_form(); + reconstructed.canonical_column_order(); + reconstructed.normal_form(); + res.first.add_phase(ap.ap_.phase - reconstructed.ap_.phase); + return res; +} + +APState tableau_to_apstate(SymplecticTableau tab) { + unsigned n_qbs = tab.get_n_qubits(); + unsigned n_rows = tab.get_n_rows(); + if (n_rows > n_qbs) + throw std::logic_error( + "tableau_to_apstate requires a tableau with up to n commuting rows for " + "n qubits"); + MatrixXb fullmat = MatrixXb::Zero(n_rows, 2 * n_qbs); + /** + * Gaussian elimination by the x matrix first ensures the bottom rows are only + * Zs, i.e. describing rows of A. Reversing the columns of the x matrix + * guarantees that each row has an X on at most one free qubit, simplifying + * the code for finding E and P + */ + for (unsigned c = 0; c < n_qbs; ++c) { + fullmat.col(c) = tab.xmat.col(n_qbs - 1 - c); + } + fullmat.block(0, n_qbs, n_rows, n_qbs) = tab.zmat; + std::vector> row_ops = + gaussian_elimination_row_ops(fullmat); + for (const std::pair& op : row_ops) + tab.row_mult(op.first, op.second); + + MatrixXb A = MatrixXb::Zero(n_qbs, n_qbs); + VectorXb B = VectorXb::Zero(n_qbs); + MatrixXb C = MatrixXb::Zero(n_qbs, n_qbs); + MatrixXb E = MatrixXb::Zero(n_qbs, n_qbs); + Eigen::VectorXi P = Eigen::VectorXi::Zero(n_qbs); + + unsigned n_leading = 0; + for (unsigned r = n_rows; r > 0;) { + --r; + bool is_a_row = true; + for (unsigned c = 0; c < n_qbs; ++c) { + if (tab.xmat(r, c)) { + is_a_row = false; + break; + } + } + if (!is_a_row) break; + ++n_leading; + } + A.block(0, 0, n_leading, n_qbs) = + tab.zmat.block(n_rows - n_leading, 0, n_leading, n_qbs); + B.head(n_leading) = tab.phase.tail(n_leading); + std::map leader_to_row; + for (unsigned r = 0; r < n_leading; ++r) { + for (unsigned c = r; c < n_qbs; ++c) { + if (A(r, c)) { + leader_to_row.insert({c, r}); + break; + } + } + } + + /** + * Each free qubit q is after all leaders connected to it, by reduced + * row-echelon of A. Therefore Gaussian elimination of the x matrix in reverse + * order gives rows corresponding to the columns of A of free qubits (plus the + * free qubit itself). + * + * Then the corresponding row q of the z matrix is the free qubit's row of E, + * plus the rows of any E for any mixed qubit connected to q in C, plus an + * extra flip at q if P(q) is odd. We can first look at the column for each + * mixed qubit to identify their rows of E and subtract from this matrix to + * leave the interactions between free qubits and their local phases. + */ + + // Start by identifying the free and mixed qubits + std::map free_to_row; + for (unsigned r = 0; r < n_qbs - n_leading; ++r) { + for (unsigned c = n_qbs - r; c > 0;) { + --c; + if (tab.xmat(r, c)) { + free_to_row.insert({c, r}); + break; + } + } + } + std::map mixed_to_row; + for (unsigned q = 0; q < n_qbs; ++q) { + if (leader_to_row.find(q) == leader_to_row.end() && + free_to_row.find(q) == free_to_row.end()) { + unsigned r = mixed_to_row.size(); + mixed_to_row.insert({q, r}); + C(r, q) = true; + } + } + // Identify C and the mixed rows of E by looking at what mixed qubits appear + // in the rows of each free qubit + for (const std::pair& fr : free_to_row) { + for (const std::pair& mr : mixed_to_row) { + if (tab.xmat(fr.second, mr.first)) C(mr.second, fr.first) = true; + if (tab.zmat(fr.second, mr.first)) { + E(mr.first, fr.first) = true; + E(fr.first, mr.first) = true; + } + } + } + // Identify connections in E between free qubits + for (auto fr1 = free_to_row.begin(); fr1 != free_to_row.end(); ++fr1) { + for (auto fr2 = free_to_row.begin(); fr2 != fr1; ++fr2) { + unsigned n_shared_mixed = 0; + for (const std::pair& mr : mixed_to_row) { + if (C(mr.second, fr1->first) && C(mr.second, fr2->first)) + ++n_shared_mixed; + } + if (tab.zmat(fr1->second, fr2->first) ^ (n_shared_mixed % 2 == 1)) { + E(fr1->first, fr2->first) = true; + E(fr2->first, fr1->first) = true; + } + } + } + // Identify P + for (const std::pair& fr : free_to_row) { + unsigned n_mixed_in_c_and_e = 0; + for (const std::pair& mr : mixed_to_row) { + if (C(mr.second, fr.first) && E(mr.first, fr.first)) ++n_mixed_in_c_and_e; + } + if (tab.zmat(fr.second, fr.first) ^ (n_mixed_in_c_and_e % 2 == 1)) + P(fr.first) += 1; + if (tab.phase(fr.second) ^ (n_mixed_in_c_and_e % 2 == 1)) P(fr.first) += 2; + } + + return APState(A, B, C, E, P, 0); +} + +SymplecticTableau apstate_to_tableau(APState ap) { + unsigned n_qbs = ap.A.cols(); + // Want A and C in reduced row-echelon form to identify leaders and mixed + // qubits, but don't need the rest in normal form + std::vector> row_ops = + gaussian_elimination_row_ops(ap.A); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + ap.A(op.second, j) ^= ap.A(op.first, j); + } + ap.B(op.second) ^= ap.B(op.first); + } + std::map leader_to_row; + for (unsigned r = 0; r < n_qbs; ++r) { + bool leader_found = false; + for (unsigned c = 0; c < n_qbs; ++c) { + if (ap.A(r, c)) { + leader_found = true; + leader_to_row.insert({c, r}); + break; + } + } + if (!leader_found) break; + } + for (unsigned r = 0; r < n_qbs; ++r) { + for (const std::pair& lr : leader_to_row) { + if (ap.C(r, lr.first)) { + for (unsigned c = 0; c < n_qbs; ++c) { + ap.C(r, c) ^= ap.A(lr.second, c); + } + } + } + } + row_ops = gaussian_elimination_row_ops(ap.C); + for (const std::pair& op : row_ops) { + for (unsigned j = 0; j < n_qbs; ++j) { + ap.C(op.second, j) ^= ap.C(op.first, j); + } + } + std::map mixed_to_row; + for (unsigned r = 0; r < n_qbs; ++r) { + bool mixed_found = false; + for (unsigned c = 0; c < n_qbs; ++c) { + if (ap.C(r, c)) { + mixed_found = true; + mixed_to_row.insert({c, r}); + break; + } + } + if (!mixed_found) break; + } + unsigned n_rows = n_qbs - mixed_to_row.size(); + MatrixXb xmat = MatrixXb::Zero(n_rows, n_qbs); + MatrixXb zmat = MatrixXb::Zero(n_rows, n_qbs); + VectorXb phase = VectorXb::Zero(n_rows); + + // One stabiliser per leader, with Z on that qubit and Z on every neighbour in + // A + unsigned n_leading = leader_to_row.size(); + zmat.block(0, 0, n_leading, n_qbs) = ap.A.block(0, 0, n_leading, n_qbs); + phase.head(n_leading) = ap.B.head(n_leading); + + // One stabiliser per free, with X on that qubit, every neighbour in C, and + // the odd neighbourhood of this set in A; pushing this through the phase + // polynomial adds Zs + unsigned r = n_leading; + for (unsigned q = 0; q < n_qbs; ++q) { + if (leader_to_row.find(q) == leader_to_row.end() && + mixed_to_row.find(q) == mixed_to_row.end()) { + // Calculate the Xs + xmat(r, q) = true; + for (const std::pair& mr : mixed_to_row) { + xmat(r, mr.first) = ap.C(mr.second, q); + } + for (const std::pair& lr : leader_to_row) { + unsigned n_mixed_in_between = 0; + for (const std::pair& mr : mixed_to_row) { + if (ap.A(lr.second, mr.first) && ap.C(mr.second, q)) + ++n_mixed_in_between; + } + xmat(r, lr.first) = ap.A(lr.second, q) ^ (n_mixed_in_between % 2 == 1); + } + // Push through the phase polynomial to calculate the Zs + for (unsigned q2 = 0; q2 < n_qbs; ++q2) { + if (xmat(r, q2)) { + // Pushing an X through each CZ creates a Z + for (unsigned q3 = 0; q3 < n_qbs; ++q3) { + if (ap.E(q2, q3)) { + zmat(r, q3) ^= true; + // If both qubits of a CZ have Xs, we will need to reorder the X + // and Z on one to match the other + if (q2 < q3 && xmat(r, q3)) phase(r) ^= true; + } + } + // Pushing an X through the local phase + zmat(r, q2) ^= ((unsigned)ap.P(q2) % 2 == 1); + phase(r) ^= ((unsigned)ap.P(q2) % 4 > 1); + } + } + ++r; + } + } + + return SymplecticTableau(xmat, zmat, phase); +} + +ChoiAPState cm_tableau_to_choi_apstate(const ChoiMixTableau& tab) { + ChoiAPState ap(0); + ap.ap_ = tableau_to_apstate(tab.tab_); + BOOST_FOREACH ( + ChoiMixTableau::tableau_col_index_t::left_const_reference entry, + tab.col_index_.left) { + ChoiAPState::TableauSegment seg = + (entry.first.second == ChoiMixTableau::TableauSegment::Input) + ? ChoiAPState::TableauSegment::Input + : ChoiAPState::TableauSegment::Output; + ap.col_index_.insert({{entry.first.first, seg}, entry.second}); + } + return ap; +} + +ChoiMixTableau choi_apstate_to_cm_tableau(const ChoiAPState& ap) { + ChoiMixTableau tab(0); + tab.tab_ = apstate_to_tableau(ap.ap_); + BOOST_FOREACH ( + ChoiAPState::tableau_col_index_t::left_const_reference entry, + ap.col_index_.left) { + ChoiMixTableau::TableauSegment seg = + (entry.first.second == ChoiAPState::TableauSegment::Input) + ? ChoiMixTableau::TableauSegment::Input + : ChoiMixTableau::TableauSegment::Output; + tab.col_index_.insert({{entry.first.first, seg}, entry.second}); + } + return tab; +} + +} // namespace tket \ No newline at end of file diff --git a/tket/test/CMakeLists.txt b/tket/test/CMakeLists.txt index fe7011d6c5..5e2d38c21c 100644 --- a/tket/test/CMakeLists.txt +++ b/tket/test/CMakeLists.txt @@ -102,7 +102,7 @@ add_executable(test-tket src/Simulation/test_CircuitSimulator.cpp src/Simulation/test_PauliExpBoxUnitaryCalculator.cpp src/test_AASRoute.cpp - src/test_APState.cpp + src/test_APState2.cpp src/test_ArchitectureAwareSynthesis.cpp src/test_Architectures.cpp src/test_Assertion.cpp diff --git a/tket/test/src/test_APState2.cpp b/tket/test/src/test_APState2.cpp new file mode 100644 index 0000000000..dcaeceeb01 --- /dev/null +++ b/tket/test/src/test_APState2.cpp @@ -0,0 +1,514 @@ +// 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 "Simulation/ComparisonFunctions.hpp" +#include "testutil.hpp" +#include "tket/Circuit/Simulation/CircuitSimulator.hpp" +#include "tket/Clifford/APState2.hpp" +#include "tket/Converters/Converters.hpp" + +namespace tket { +namespace test_APState2 { + +void test_apply_gate( + const MatrixXb& A, const VectorXb& B, const MatrixXb& E, + const Eigen::VectorXi& P, OpType ot, const std::vector& args) { + MatrixXb C = MatrixXb::Zero(A.cols(), A.cols()); + APState ap(A, B, C, E, P, 0); + ap.verify(); + auto sv_before = ap.to_statevector(); + + Circuit circ(A.cols()); + circ.add_op(ot, args); + auto gate_u = tket_sim::get_unitary(circ); + + ap.apply_gate(ot, args); + + ap.verify(); + auto sv_after = ap.to_statevector(); + + CHECK(tket_sim::compare_statevectors_or_unitaries( + gate_u * sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL)); +} + +void test_apply_gate_dm( + const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E, + const Eigen::VectorXi& P, OpType ot, const std::vector& args) { + APState ap(A, B, C, E, P, 0); + ap.verify(); + auto dm_before = ap.to_density_matrix(); + + Circuit circ(A.cols()); + circ.add_op(ot, args); + auto gate_u = tket_sim::get_unitary(circ); + + ap.apply_gate(ot, args); + + ap.verify(); + auto dm_after = ap.to_density_matrix(); + + CHECK(dm_after.isApprox(gate_u * dm_before * gate_u.adjoint(), EPS)); +} + +SCENARIO("CZ cases") { + GIVEN("CZ on free qubits") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 3) = true; + E(2, 3) = E(3, 2) = true; + test_apply_gate(A, B, E, P, OpType::CZ, {1, 2}); + } + GIVEN("CZ on one leading qubit and connected free") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = A(1, 4) = true; + B(1) = (b == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 3}); + } + } + GIVEN("CZ on one leading qubit and unconnected free") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = A(1, 4) = true; + B(1) = (b == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 4}); + } + } + GIVEN("CZ on leading qubits") { + for (unsigned b1 = 0; b1 < 2; ++b1) { + for (unsigned b2 = 0; b2 < 2; ++b2) { + MatrixXb A = MatrixXb::Zero(8, 8); + VectorXb B = VectorXb::Zero(8); + MatrixXb E = MatrixXb::Zero(8, 8); + Eigen::VectorXi P = Eigen::VectorXi::Zero(8); + A(0, 0) = A(0, 2) = A(0, 3) = A(0, 4) = A(0, 5) = true; + A(1, 1) = A(1, 4) = A(1, 5) = A(1, 6) = A(1, 7) = true; + B(0) = (b1 == 1); + B(1) = (b2 == 1); + test_apply_gate(A, B, E, P, OpType::CZ, {0, 1}); + } + } + } + GIVEN("CZ on mixed state") { + for (unsigned b1 = 0; b1 < 2; ++b1) { + for (unsigned b2 = 0; b2 < 2; ++b2) { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb C = MatrixXb::Zero(5, 5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 1) = A(0, 3) = A(0, 4) = true; + A(1, 0) = A(1, 2) = A(1, 3) = true; + B(0) = (b1 == 1); + B(1) = (b1 == 1); + C(0, 0) = C(0, 2) = true; + E(0, 3) = E(3, 0) = true; + E(1, 3) = E(3, 1) = true; + P(0) = 3; + P(1) = 1; + test_apply_gate_dm(A, B, C, E, P, OpType::CZ, {0, 1}); + } + } + } +} + +SCENARIO("S cases") { + GIVEN("S on free qubit") { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + A(0, 0) = A(0, 1) = A(0, 2) = true; + E(1, 2) = E(2, 1) = true; + test_apply_gate(A, B, E, P, OpType::S, {2}); + } + GIVEN("S on leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 2) = true; + B(0) = (b == 1); + E(1, 3) = E(3, 1) = true; + test_apply_gate(A, B, E, P, OpType::S, {0}); + } + } + GIVEN("S on disconnected leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(1, 1); + VectorXb B = VectorXb::Zero(1); + MatrixXb E = MatrixXb::Zero(1, 1); + Eigen::VectorXi P = Eigen::VectorXi::Zero(1); + A(0, 0) = true; + B(0) = (b == 1); + test_apply_gate(A, B, E, P, OpType::S, {0}); + } + } + GIVEN("S on a mixed state") { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb C = MatrixXb::Zero(3, 3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + A(0, 0) = A(0, 1) = A(0, 2) = true; + C(0, 1) = C(0, 2) = true; + E(1, 2) = E(2, 1) = true; + test_apply_gate_dm(A, B, C, E, P, OpType::S, {2}); + } +} + +SCENARIO("V cases") { + GIVEN("V on leading qubit") { + for (unsigned b = 0; b < 2; ++b) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 3) = true; + B(0) = (b == 1); + test_apply_gate(A, B, E, P, OpType::V, {0}); + } + } + GIVEN("V on free qubit with some leading") { + for (unsigned b = 0; b < 2; ++b) { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(9, 9); + VectorXb B = VectorXb::Zero(9); + MatrixXb E = MatrixXb::Zero(9, 9); + Eigen::VectorXi P = Eigen::VectorXi::Zero(9); + A(0, 0) = A(0, 2) = A(0, 4) = A(0, 5) = A(0, 7) = true; + A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = A(1, 5) = A(1, 6) = true; + B(1) = (b == 1); + E(4, 5) = E(5, 4) = true; + E(4, 6) = E(6, 4) = true; + E(4, 7) = E(7, 4) = true; + E(4, 8) = E(8, 4) = true; + P(4) = p; + test_apply_gate(A, B, E, P, OpType::V, {4}); + } + } + } + GIVEN("V on free qubit with some earlier connected free") { + for (unsigned p1 = 0; p1 < 4; ++p1) { + for (unsigned p2 = 0; p2 < 4; ++p2) { + MatrixXb A = MatrixXb::Zero(9, 9); + VectorXb B = VectorXb::Zero(9); + MatrixXb E = MatrixXb::Zero(9, 9); + Eigen::VectorXi P = Eigen::VectorXi::Zero(9); + A(0, 0) = true; + A(1, 1) = A(1, 4) = A(1, 6) = true; + A(2, 2) = A(2, 4) = A(2, 7) = true; + A(3, 3) = A(3, 4) = A(3, 8) = A(3, 8) = true; + E(4, 5) = E(5, 4) = true; + E(4, 7) = E(7, 4) = true; + E(4, 8) = E(8, 4) = true; + E(5, 6) = E(6, 5) = true; + E(5, 7) = E(7, 5) = true; + E(5, 8) = E(8, 5) = true; + P(4) = p1; + P(5) = p2; + test_apply_gate(A, B, E, P, OpType::V, {5}); + } + } + } + GIVEN("V on free qubit with no earlier connected free") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = true; + E(1, 2) = E(2, 1) = true; + E(1, 3) = E(3, 1) = true; + P(1) = p; + test_apply_gate(A, B, E, P, OpType::V, {1}); + } + } + GIVEN("V on disconnected free qubit") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(1, 1); + VectorXb B = VectorXb::Zero(1); + MatrixXb E = MatrixXb::Zero(1, 1); + Eigen::VectorXi P = Eigen::VectorXi::Zero(1, 1); + P(0) = p; + test_apply_gate(A, B, E, P, OpType::V, {0}); + } + } + GIVEN("V on a qubit involved in A and C") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(7, 7); + VectorXb B = VectorXb::Zero(7); + MatrixXb C = MatrixXb::Zero(7, 7); + MatrixXb E = MatrixXb::Zero(7, 7); + Eigen::VectorXi P = Eigen::VectorXi::Zero(7); + A(0, 0) = A(0, 2) = true; + A(1, 0) = A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = true; + C(0, 0) = C(0, 2) = C(0, 3) = true; + C(1, 0) = C(1, 1) = C(1, 2) = C(1, 5) = C(1, 6) = true; + E(0, 3) = E(3, 0) = true; + E(0, 4) = E(4, 0) = true; + E(0, 5) = E(5, 0) = true; + E(0, 6) = E(6, 0) = true; + P(0) = p; + test_apply_gate_dm(A, B, C, E, P, OpType::V, {0}); + } + } +} + +SCENARIO("Gate encodings") { + std::list>> test_gates = { + {OpType::Z, {0}}, {OpType::X, {0}}, + {OpType::Y, {0}}, {OpType::S, {0}}, + {OpType::Sdg, {0}}, {OpType::V, {0}}, + {OpType::Vdg, {0}}, {OpType::SX, {0}}, + {OpType::SXdg, {0}}, {OpType::H, {0}}, + {OpType::CX, {0, 1}}, {OpType::CY, {0, 1}}, + {OpType::CZ, {0, 1}}, {OpType::ZZMax, {0, 1}}, + {OpType::ECR, {0, 1}}, {OpType::ISWAPMax, {0, 1}}, + {OpType::SWAP, {0, 1}}, {OpType::BRIDGE, {0, 1, 2}}, + {OpType::noop, {0}}, + }; + GIVEN("Check Z actions") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Identity(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + test_apply_gate(A, B, E, P, com.first, com.second); + } + } + GIVEN("Check X actions") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + test_apply_gate(A, B, E, P, com.first, com.second); + } + } +} + +SCENARIO("Loading from a statevector") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 2) = true; + B(0) = true; + E(2, 3) = E(3, 2) = true; + P(2) = 1; + P(3) = 2; + APState ap(A, B, C, E, P, 0); + auto sv = ap.to_statevector(); + APState reconstructed(sv); + auto sv2 = reconstructed.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv, sv2, tket_sim::MatrixEquivalence::EQUAL)); +} + +SCENARIO("Loading from a density matrix") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 2) = A(0, 3) = true; + C(0, 1) = C(0, 2) = true; + C(1, 0) = true; + B(0) = true; + E(2, 3) = E(3, 2) = true; + P(0) = 3; + P(1) = 2; + P(2) = 1; + P(3) = 2; + APState ap(A, B, C, E, P, 0); + auto dm = ap.to_density_matrix(); + APState reconstructed(dm); + auto dm2 = reconstructed.to_density_matrix(); + CHECK(dm2.isApprox(dm, EPS)); + // This state is mixed, so only check equality of normal forms up to phase + ap.normal_form(); + reconstructed.normal_form(); + reconstructed.phase = ap.phase; + CHECK(ap == reconstructed); + THEN("Test serialisation") { + nlohmann::json j_ap = ap; + APState ap2(0); + j_ap.get_to(ap2); + REQUIRE(ap == ap2); + } +} + +SCENARIO("Converting from/to a circuit") { + GIVEN("A pure circuit in the standard AP form") { + Circuit circ(4); + circ.qubit_create_all(); + circ.add_op(OpType::H, {2}); + circ.add_op(OpType::H, {3}); + circ.add_op(OpType::CX, {2, 0}); + circ.add_op(OpType::CX, {2, 1}); + circ.add_op(OpType::CX, {3, 1}); + circ.add_op(OpType::S, {2}); + circ.add_op(OpType::Z, {3}); + APState ap = circuit_to_apstate(circ); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + ap.normal_form(); + sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + Circuit reconstructed = apstate_to_circuit(ap); + CHECK(circ == reconstructed); + } + GIVEN("A generic pure circuit") { + Circuit circ(4); + circ.qubit_create_all(); + circ.add_op(OpType::V, {0}); + circ.add_op(OpType::CX, {0, 1}); + circ.add_op(OpType::CY, {1, 3}); + circ.add_op(OpType::H, {3}); + circ.add_op(OpType::ZZMax, {2, 3}); + APState ap = circuit_to_apstate(circ); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + Circuit reconstructed = apstate_to_circuit(ap); + auto sv_rec = tket_sim::get_statevector(reconstructed); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_rec, tket_sim::MatrixEquivalence::EQUAL)); + } + 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)); + ChoiAPState ap = circuit_to_choi_apstate(circ); + ap.post_select(Qubit(3), ChoiAPState::TableauSegment::Output); + ap.canonical_column_order(); + ap.normal_form(); + // Define correct form from a hand calculation + MatrixXb A = MatrixXb::Zero(6, 6); + VectorXb B = VectorXb::Zero(6); + MatrixXb C = MatrixXb::Zero(6, 6); + C(0, 0) = C(0, 3) = true; + C(1, 1) = true; + MatrixXb E = MatrixXb::Zero(6, 6); + E(1, 2) = E(2, 1) = true; + E(1, 4) = E(4, 1) = true; + E(1, 5) = E(5, 1) = true; + E(3, 4) = E(4, 3) = true; + Eigen::VectorXi P = Eigen::VectorXi::Zero(6); + P(3) = P(4) = 3; + // Ignore phase by setting them to match + APState correct(A, B, C, E, P, ap.ap_.phase); + CHECK(ap.ap_ == correct); + THEN("Serialize and deserialize") { + nlohmann::json j_ap = ap; + ChoiAPState ap2({}); + j_ap.get_to(ap2); + CHECK(ap == ap2); + } + THEN("Check conversion to/from a tableau") { + ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap); + ChoiMixTableau tab2 = circuit_to_cm_tableau(circ); + tab2.post_select(Qubit(3), ChoiMixTableau::TableauSegment::Output); + tab.canonical_column_order(); + tab.gaussian_form(); + tab2.canonical_column_order(); + tab2.gaussian_form(); + REQUIRE(tab == tab2); + ChoiAPState ap2 = cm_tableau_to_choi_apstate(tab); + ap2.canonical_column_order(); + ap2.normal_form(); + // Converting to a tableau drops phase, so ignore this in equivalence + // check + ap2.ap_.phase = ap.ap_.phase; + REQUIRE(ap == ap2); + } + } +} + +SCENARIO("Converting from/to a tableau") { + GIVEN("Check pure state up to global phase using circuit") { + Circuit circ(8); + circ.qubit_create_all(); + circ.add_op(OpType::X, {1}); + circ.add_op(OpType::X, {5}); + circ.add_op(OpType::H, {2}); + circ.add_op(OpType::H, {4}); + circ.add_op(OpType::H, {6}); + circ.add_op(OpType::H, {7}); + circ.add_op(OpType::CX, {2, 1}); + circ.add_op(OpType::CX, {4, 0}); + circ.add_op(OpType::CX, {4, 3}); + circ.add_op(OpType::CX, {6, 0}); + circ.add_op(OpType::CX, {6, 1}); + circ.add_op(OpType::CX, {7, 5}); + circ.add_op(OpType::CZ, {2, 6}); + circ.add_op(OpType::CZ, {4, 6}); + circ.add_op(OpType::CZ, {4, 7}); + circ.add_op(OpType::CZ, {6, 7}); + circ.add_op(OpType::S, {2}); + circ.add_op(OpType::Sdg, {4}); + circ.add_op(OpType::Z, {7}); + ChoiMixTableau cmt = circuit_to_cm_tableau(circ); + ChoiAPState ap = cm_tableau_to_choi_apstate(cmt); + auto sv_circ = tket_sim::get_statevector(circ); + auto sv_ap = ap.ap_.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); + ChoiMixTableau cmt2 = choi_apstate_to_cm_tableau(ap); + auto [circ2, perm] = cm_tableau_to_exact_circuit(cmt2); + qubit_map_t inv; + for (const std::pair& qp : perm) + inv.insert({qp.second, qp.first}); + circ2.permute_boundary_output(inv); + auto sv_circ2 = tket_sim::get_statevector(circ2); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_circ, sv_circ2, + tket_sim::MatrixEquivalence::EQUAL_UP_TO_GLOBAL_PHASE)); + } +} + +} // namespace test_APState2 +} // namespace tket \ No newline at end of file From 7d48368c9d97ea30804c6bcfc497db5d31891967 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 14:49:42 +0000 Subject: [PATCH 09/19] Remove old version --- tket/include/tket/Clifford/APState.hpp | 146 ----- tket/src/Clifford/APState.cpp | 861 ------------------------- tket/test/src/test_APState.cpp | 347 ---------- 3 files changed, 1354 deletions(-) delete mode 100644 tket/include/tket/Clifford/APState.hpp delete mode 100644 tket/src/Clifford/APState.cpp delete mode 100644 tket/test/src/test_APState.cpp diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp deleted file mode 100644 index b1eb808371..0000000000 --- a/tket/include/tket/Clifford/APState.hpp +++ /dev/null @@ -1,146 +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 "tket/OpType/OpType.hpp" -#include "tket/Utils/MatrixAnalysis.hpp" - -namespace tket { - -/** - * APState class for Clifford states with global phase tracking - * - * The "affine with phases" form of a Clifford state from ZX calculus (see - * Kissinger & van de Wetering, "Picturing Quantum Software") represents n-qubit - * Clifford states uniquely with the following data: - * - A binary (k,n) matrix A in reduced row echelon form. The leading columns of - * each row are referred to as "leading qubits", and the others as "free - * qubits". We apply H gates to all free qubits, then for each row of A we apply - * a CX targeted on the leading qubit and control by each other entry in the - * row. - * - A binary (k) vector B describing X gates applied to the leading qubits. - * - A Clifford phase polynomial Phi over the free qubits. Wlog this can be in - * the form of a symmetric, zero-diagonal, binary (n-k,n-k) matrix E describing - * CZ gates and a (n-k) vector P of integers mod 4 describing S gates. - * - * This gives a canonical statevector: - * \sum_{x, Ax=B} i^{\Phi(x)} |x> - * - * This canonical statevector fixes a reference phase from which we can track - * global phase with an additional parameter. - * - * When applying gates, qubits may switch between being leading and free - * (including those that aren't involved in the gate due to the need to reduce - * to the canonical form, e.g. reduced row echelon form for A). The updates are - * easiest to prove in the ZX calculus for the gateset (CZ, S i.e. pi/2 green - * phase, SX i.e. pi/2 red phase). - * - * To save on resizing the matrices and vectors, we will keep them at full size - * and just assert correctness that, e.g. every entry in A is either on the - * diagonal (to indicate the qubit is leading) or in the row of a leading qubit - * and a column of a later following qubit. - */ -class APState { - public: - /** - * An upper triangular matrix whose diagonal entries indicate leading qubits - * and off-diagonal entries represent a CX gate from the column qubit - * (necessarily a free qubit) to the row qubit. - * - * If a diagonal entry is 0 (it is a free qubit), then every entry in the row - * is 0. If a diagonal entry is 1 (it is a leading qubit), all other entries - * in the column are 0. - */ - MatrixXb A; - - /** - * A vector indicating X gates on leading qubits. - * - * Must be 0 on every free qubit. - */ - VectorXb B; - - /** - * A symmetric, zero diagonal matrix whose entries indicate CZs between free - * qubits. - * - * Every diagonal entry must be 0. - * The column and row for each leading qubit must be 0. - */ - MatrixXb E; - - /** - * A vector indicating S^{P(i)} on qubit i which must be free. - * - * Must be 0 on every leading qubit. - */ - Eigen::VectorXi P; - - /** - * The global phase term (in half-turns). - */ - Expr phase; - - /** - * Constructs a state in AP form from given data. - */ - APState( - const MatrixXb& A_, const VectorXb& B_, const MatrixXb& E_, - const Eigen::VectorXi& P_, const Expr& phase_); - - /** - * Constructs the state |0>^{x n_qubits} in AP form. - */ - APState(unsigned n_qubits); - - /** - * Constructs the state in AP form from a given statevector. - */ - APState(const Eigen::VectorXcd& sv); - - /** - * Verifies the internal correctness of the data structure. Throws an - * exception if the data structure is invalid. - */ - void verify(); - - /** - * Calculates the statevector of the state. - */ - Eigen::VectorXcd to_statevector(); - - /** - * Applies a CZ gate to the state. O(n^2) wrt number of qubits in the state. - */ - void apply_CZ(unsigned ctrl, unsigned trgt); - - /** - * Applies an S gate to the state. O(n^2) wrt number of qubits in the state. - */ - void apply_S(unsigned q); - - /** - * Applies a V gate to the state. O(n^2) wrt number of qubits in the state. - */ - void apply_V(unsigned q); - - /** - * Applies an unparameterised Clifford gate to the state on the chosen qubits. - * O(n^2) wrt number of qubits in the state. - */ - void apply_gate(OpType type, const std::vector& qbs); -}; - -} // namespace tket \ No newline at end of file diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp deleted file mode 100644 index 2dfa489611..0000000000 --- a/tket/src/Clifford/APState.cpp +++ /dev/null @@ -1,861 +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/Clifford/APState.hpp" - -#include - -#include "tket/OpType/OpTypeInfo.hpp" - -namespace tket { - -APState::APState( - const MatrixXb& A_, const VectorXb& B_, const MatrixXb& E_, - const Eigen::VectorXi& P_, const Expr& phase_) - : A(A_), B(B_), E(E_), P(P_), phase(phase_) { - verify(); -} - -APState::APState(unsigned n_qubits) - : A(MatrixXb::Identity(n_qubits, n_qubits)), - B(VectorXb::Zero(n_qubits)), - E(MatrixXb::Zero(n_qubits, n_qubits)), - P(Eigen::VectorXi(n_qubits)), - phase(0) {} - -unsigned clifford_phase(const Complex& c) { - unsigned res = 0; - unsigned n_possibles = 0; - if (c.real() > EPS) { - res = 0; - ++n_possibles; - } - if (c.imag() > EPS) { - res = 1; - ++n_possibles; - } - if (c.real() < -EPS) { - res = 2; - ++n_possibles; - } - if (c.imag() < -EPS) { - res = 3; - ++n_possibles; - } - TKET_ASSERT(n_possibles == 1); - return res; -} - -APState::APState(const Eigen::VectorXcd& sv) { - unsigned n_qbs = 0; - while (sv.size() > (1 << n_qbs)) ++n_qbs; - TKET_ASSERT(sv.size() == (1 << n_qbs)); - - // Find non-zero entries as a vector space and offset - unsigned z0 = 0; - for (unsigned x = 0; x < sv.size(); ++x) { - if (sv(x) != 0.) { - z0 = x; - break; - } - } - std::vector offsets; - unsigned n_non_zero = 0; - for (unsigned x = 1; x < sv.size(); ++x) { - if (sv(z0 ^ x) != 0.) { - ++n_non_zero; - if (n_non_zero == (1u << offsets.size())) offsets.push_back(x); - } - } - - // Find A as the dual space - MatrixXb offset_mat(offsets.size(), n_qbs); - for (unsigned r = 0; r < offsets.size(); ++r) { - unsigned off = offsets.at(r); - for (unsigned c = 0; c < n_qbs; ++c) { - // Binary encoding of offsets in reverse order to guarantee free qubits - // are the later ones, meaning we produce A in row echelon form - offset_mat(r, c) = (off & (1 << c)) != 0; - } - } - std::vector> row_ops = - gaussian_elimination_row_ops(offset_mat); - for (const std::pair& op : row_ops) { - for (unsigned j = 0; j < n_qbs; ++j) { - offset_mat(op.second, j) ^= offset_mat(op.first, j); - } - } - std::map mat_leaders; - A = MatrixXb::Zero(n_qbs, n_qbs); - for (unsigned c = 0; c < n_qbs; ++c) { - bool free_qubit = false; - for (unsigned r = 0; r < offsets.size(); ++r) { - if (offset_mat(r, c)) { - auto inserted = mat_leaders.insert({r, c}); - if (inserted.second) { - free_qubit = true; - break; - } else { - // Reverse bit orderings back to normal - A(n_qbs - 1 - c, n_qbs - 1 - inserted.first->second) = true; - } - } - } - A(n_qbs - 1 - c, n_qbs - 1 - c) = !free_qubit; - } - - VectorXb z0_vec(n_qbs); - for (unsigned i = 0; i < n_qbs; ++i) { - z0_vec(i) = ((z0 & (1 << (n_qbs - 1 - i))) != 0); - } - B = A * z0_vec; - - E = MatrixXb::Zero(n_qbs, n_qbs); - P = Eigen::VectorXi::Zero(n_qbs); - unsigned neutral_z0 = z0; // Index with 0 for all free qubits - std::map offset_for_free; - for (const std::pair& row_qfree : mat_leaders) { - unsigned offset = 0; - for (unsigned i = 0; i < n_qbs; ++i) - if (offset_mat(row_qfree.first, i)) offset += (1 << i); - unsigned qfree = n_qbs - 1 - row_qfree.second; - offset_for_free.insert({qfree, offset}); - if ((neutral_z0 & (1 << row_qfree.second)) != 0) - neutral_z0 = neutral_z0 ^ offset; - } - for (const std::pair& qfree_offset : - offset_for_free) { - unsigned qfree = qfree_offset.first; - unsigned offset = qfree_offset.second; - auto complex_phase = sv(neutral_z0 ^ offset) / sv(neutral_z0); - P(qfree) = clifford_phase(complex_phase); - for (const std::pair& qfree_offset2 : - offset_for_free) { - if (qfree_offset == qfree_offset2) break; // Only go up to solved phases - unsigned qfree2 = qfree_offset2.first; - unsigned offset2 = qfree_offset2.second; - auto complex_phase_cz = - sv(neutral_z0 ^ offset ^ offset2) / sv(neutral_z0); - E(qfree, qfree2) = E(qfree2, qfree) = - ((clifford_phase(complex_phase_cz) - P(qfree) - P(qfree2)) % 4) == 2; - } - } - - phase = Expr(std::arg(sv(neutral_z0)) / PI); -} - -void APState::verify() { - // Check dimensions all agree - unsigned n_qubits = A.rows(); - TKET_ASSERT(A.cols() == n_qubits); - TKET_ASSERT(B.size() == n_qubits); - TKET_ASSERT(E.rows() == n_qubits); - TKET_ASSERT(E.cols() == n_qubits); - TKET_ASSERT(P.size() == n_qubits); - - for (unsigned r = 0; r < n_qubits; ++r) { - for (unsigned c = 0; c < r; ++c) { - // Check A is upper triangular - TKET_ASSERT(!A(r, c)); - // Check E is symmetric - TKET_ASSERT(E(r, c) == E(c, r)); - } - if (A(r, r)) { - // Check leaders are eliminated in A (row echelon) - for (unsigned r2 = 0; r2 < r; ++r2) { - TKET_ASSERT(!A(r2, r)); - } - } else { - // Check A only has entries in rows with a leader - for (unsigned c = r + 1; c < n_qubits; ++c) { - TKET_ASSERT(!A(r, c)); - } - } - } - - for (unsigned q = 0; q < n_qubits; ++q) { - if (A(q, q)) { - // Check P is zero on leaders - TKET_ASSERT(P(q) % 4 == 0); - // Check E is zero on leaders (simplified by symmetry) - for (unsigned q2 = 0; q2 < n_qubits; ++q2) { - TKET_ASSERT(!E(q, q2)); - } - } else { - // Check B is zero on free qubits - TKET_ASSERT(!B(q)); - // Check diagonals of E are still zero - TKET_ASSERT(!E(q, q)); - } - } -} - -VectorXb z2_mult(const MatrixXb& M, const VectorXb& v) { - VectorXb res = VectorXb::Zero(M.rows()); - for (unsigned i = 0; i < M.cols(); ++i) { - if (v(i)) { - for (unsigned j = 0; j < M.rows(); ++j) res(j) ^= M(j, i); - } - } - return res; -} - -Eigen::VectorXcd APState::to_statevector() { - unsigned n_qubits = A.cols(); - Eigen::VectorXcd sv = Eigen::VectorXcd::Zero(1 << n_qubits); - unsigned n_terms = 0; - Complex g_phase = std::exp(i_ * PI * eval_expr(phase).value()); - for (unsigned x = 0; x < (1u << n_qubits); ++x) { - VectorXb x_binary = VectorXb::Zero(n_qubits); - for (unsigned i = 0; i < n_qubits; ++i) { - unsigned mask = 1 << (n_qubits - 1 - i); - x_binary(i) = ((x & mask) != 0); - } - - if (z2_mult(A, x_binary) == B) { - ++n_terms; - unsigned i_phases = 0; - for (unsigned q = 0; q < n_qubits; ++q) { - if (x_binary(q)) { - i_phases += P(q); - for (unsigned q2 = q; q2 < n_qubits; ++q2) { - if (E(q, q2) && x_binary(q2)) i_phases += 2; - } - } - } - switch (i_phases % 4) { - case 0: { - sv(x) = g_phase; - break; - } - case 1: { - sv(x) = i_ * g_phase; - break; - } - case 2: { - sv(x) = -g_phase; - break; - } - case 3: { - sv(x) = -i_ * g_phase; - break; - } - default: { - TKET_ASSERT(false); - } - } - } - } - return pow(n_terms, -0.5) * sv; -} - -void APState::apply_CZ(unsigned ctrl, unsigned trgt) { - if (ctrl == trgt) - throw std::logic_error( - "APState::apply_CZ: cannot apply CZ when control and target qubits " - "match"); - unsigned n_qubits = A.cols(); - if (ctrl >= n_qubits || trgt >= n_qubits) - throw std::logic_error("APState::apply_CZ: qubit indices out of range"); - - if (A(ctrl, ctrl)) { - if (A(trgt, trgt)) { - // ctrl and trgt are both leading - - std::set just_ctrl, just_trgt, both; - // Add phases and determine connectivity - for (unsigned q = 0; q < n_qubits; ++q) { - if (A(ctrl, q)) { - if (B(trgt)) P(q) += 2; - if (A(trgt, q)) { - if (!B(ctrl)) P(q) += 2; - both.insert(q); - } else { - just_ctrl.insert(q); - } - } else if (A(trgt, q)) { - if (B(ctrl)) P(q) += 2; - just_trgt.insert(q); - } - } - // ctrl and trgt were included in the previous loop, so reset them - just_ctrl.erase(ctrl); - just_trgt.erase(trgt); - P(ctrl) = 0; - P(trgt) = 0; - // Pivoting-style update between those connected to ctrl, trgt, and both - for (const unsigned q1 : just_ctrl) { - for (const unsigned q2 : just_trgt) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - for (const unsigned q2 : both) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - } - for (const unsigned q1 : just_trgt) { - for (const unsigned q2 : both) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - } - - if (B(ctrl) && B(trgt)) phase += 1; - } else { - // ctrl is leading, trgt is free - - // Flip connectivity between trgt and every q on ctrl - for (unsigned q = ctrl + 1; q < n_qubits; ++q) { - if (A(ctrl, q)) { - E(trgt, q) ^= true; - E(q, trgt) ^= true; - // Note that if q = trgt, this does nothing, preserving the zero - // diagonal of E - } - } - // Add phase to trgt - if (A(ctrl, trgt) ^ B(ctrl)) P(trgt) += 2; - - // No global phase change - } - } else { - if (A(trgt, trgt)) { - // trgt is leading, ctrl is free - - // CZ is symmetric, so reuse the previous case - apply_CZ(trgt, ctrl); - } else { - // ctrl and trgt are both free - - // Just add the CZ gate to E - E(ctrl, trgt) ^= true; - E(trgt, ctrl) ^= true; - - // No global phase change - } - } -} - -void APState::apply_S(unsigned q) { - unsigned n_qubits = A.cols(); - if (q >= n_qubits) - throw std::logic_error("APState::apply_S: qubit index out of range"); - - if (A(q, q)) { - // q is leading - - // Local complementation update to E - for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { - if (A(q, q1)) { - for (unsigned q2 = q + 1; q2 < n_qubits; ++q2) { - E(q1, q2) ^= ((q1 != q2) && A(q, q2)); - } - } - } - - // Update global phase - if (B(q)) phase += .5; - - // Update local phases - unsigned local_phase_change = B(q) ? 3 : 1; - for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { - if (A(q, q1)) P(q1) += local_phase_change; - } - } else { - // q is free - - // Add to local phase - P(q) += 1; - - // No global phase change - } -} - -void APState::apply_V(unsigned q) { - unsigned n_qubits = A.cols(); - if (q >= n_qubits) - throw std::logic_error("APState::apply_V: qubit index out of range"); - - if (A(q, q)) { - // q is leading, but will become free - - // Local complementation for E - for (unsigned q1 = q; q1 < n_qubits; ++q1) { - if (A(q, q1)) { - for (unsigned q2 = q1 + 1; q2 < n_qubits; ++q2) { - E(q1, q2) ^= A(q, q2); - E(q2, q1) ^= A(q, q2); - } - - // Local phases - P(q1) += B(q) ? 1 : 3; - } - } - - // Global phase update - if (B(q)) phase += -.5; - - // q is no longer a leader - for (unsigned q1 = q; q1 < n_qubits; ++q1) { - A(q, q1) = false; - } - B(q) = false; - } else { - // q is free - - std::optional last_connected_leader; - for (unsigned q1 = q; q1 > 0;) { - --q1; - if (A(q1, q)) { - last_connected_leader = q1; - break; - } - } - - if (last_connected_leader.has_value()) { - // q is free and connected to a leader - - // For each other leader connected to q, add LCL's row in A - std::list other_leaders; - for (unsigned q1 = 0; q1 < *last_connected_leader; ++q1) { - if (A(q1, q)) { - other_leaders.push_back(q1); - for (unsigned q2 = *last_connected_leader; q2 < n_qubits; ++q2) { - A(q1, q2) ^= A(*last_connected_leader, q2); - } - } - } - - // Sort the connected frees into just q, just LCL, and both - std::list just_q, just_lcl, both; - for (unsigned q1 = 0; q1 < n_qubits; ++q1) { - if (q1 == q || q1 == *last_connected_leader) continue; - if (E(q, q1)) { - if (A(*last_connected_leader, q1)) - both.push_back(q1); - else - just_q.push_back(q1); - } else if (A(*last_connected_leader, q1)) - just_lcl.push_back(q1); - } - - // Split case by phase on q - // In each case: - // - Perform appropriate complementations between connected frees - // - Set LCL's connectivity in E appropriately - // - Modify q's connectivity in E - // - Add all local and global phases by a case matrix - if (P(q) % 2 == 0) { - // Handle complementations between neighbours - for (const unsigned q1 : just_lcl) { - // Local complement within group - for (const unsigned q2 : just_lcl) { - E(q1, q2) ^= true; - } - E(q1, q1) ^= true; - // Complement between groups - for (const unsigned q2 : just_q) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - // Add connections with q - E(q, q1) ^= true; - E(q1, q) ^= true; - } - for (const unsigned q1 : both) { - // Local complement within group - for (const unsigned q2 : both) { - E(q1, q2) ^= true; - } - E(q1, q1) ^= true; - // Complement between groups - for (const unsigned q2 : just_q) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - } - // Remove connections with q - for (const unsigned q1 : just_q) { - E(q, q1) ^= true; - E(q1, q) ^= true; - } - - // Move LCL from A to E - E(q, *last_connected_leader) = true; - E(*last_connected_leader, q) = true; - for (const unsigned q1 : just_q) { - E(q1, *last_connected_leader) = true; - E(*last_connected_leader, q1) = true; - } - for (const unsigned q1 : just_lcl) { - E(q1, *last_connected_leader) = true; - E(*last_connected_leader, q1) = true; - A(*last_connected_leader, q1) = false; - } - for (const unsigned q1 : both) A(*last_connected_leader, q1) = false; - A(*last_connected_leader, q) = false; - A(*last_connected_leader, *last_connected_leader) = false; - - // Phases - bool a = (P(q) % 4 == 2); - bool b = B(*last_connected_leader); - P(q) = b ? 1 : 3; - P(*last_connected_leader) = (a ^ b) ? 1 : 3; - B(*last_connected_leader) = false; - for (const unsigned q1 : other_leaders) B(q1) ^= b; - for (const unsigned q1 : both) P(q1) += a ? 3 : 1; - for (const unsigned q1 : just_lcl) P(q1) += (b ^ a) ? 1 : 3; - for (const unsigned q1 : just_q) P(q1) += b ? 2 : 0; - if (b) phase += a ? .5 : -.5; - } else { - // Handle complementations between neighbours - for (const unsigned q1 : just_lcl) { - // Complement between groups - for (const unsigned q2 : just_q) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - for (const unsigned q2 : both) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - // Add connections with q - E(q, q1) = true; - E(q1, q) = true; - } - for (const unsigned q1 : just_q) { - // Complement between groups - for (const unsigned q2 : both) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - // Remove connections with q - E(q, q1) = false; - E(q1, q) = false; - } - - // Move LCL from A to E - E(q, *last_connected_leader) = true; - E(*last_connected_leader, q) = true; - for (const unsigned q1 : just_q) { - E(q1, *last_connected_leader) = true; - E(*last_connected_leader, q1) = true; - } - for (const unsigned q1 : both) { - E(q1, *last_connected_leader) = true; - E(*last_connected_leader, q1) = true; - A(*last_connected_leader, q1) = false; - } - for (const unsigned q1 : just_lcl) - A(*last_connected_leader, q1) = false; - A(*last_connected_leader, q) = false; - A(*last_connected_leader, *last_connected_leader) = false; - - // Phases - bool a = (P(q) % 4 == 3); - bool b = B(*last_connected_leader); - P(q) = b ? 1 : 3; - P(*last_connected_leader) = a ? 2 : 0; - B(*last_connected_leader) = false; - for (const unsigned q1 : other_leaders) B(q1) ^= b; - for (const unsigned q1 : just_q) P(q1) += b ? 2 : 0; - for (const unsigned q1 : just_lcl) P(q1) += a ? 2 : 0; - for (const unsigned q1 : both) P(q1) += (a ^ b) ? 0 : 2; - if (a && b) phase += 1.; - } - } else if (P(q) % 2 == 0) { - // q has phase 0/pi and no connected leader - - // Local complementation update to E - std::list connected; - for (unsigned q1 = 0; q1 < n_qubits; ++q1) { - if (E(q, q1)) { - for (const unsigned q2 : connected) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - connected.push_back(q1); - - // Add local phase - P(q1) += (P(q) % 4 == 0) ? 1 : 3; - } - } - - // Local phase on q remains the same - - // Global phase - phase += (P(q) % 4 == 0) ? -0.25 : 0.25; - } else { - // q has phase +-pi/2 and no connected leader - - std::optional first_connected_free = std::nullopt; - for (unsigned q1 = 0; q1 < q; ++q1) { - if (E(q, q1)) { - first_connected_free = q1; - break; - } - } - - if (first_connected_free.has_value()) { - // q has phase +-pi/2, no connected leader, but a previous free which - // will become the new leader - - // Make FCF a leader - for (unsigned q1 = *first_connected_free; q1 < n_qubits; ++q1) { - if (E(q, q1)) { - A(*first_connected_free, q1) = true; - // q ends up disconnected in E - E(q, q1) = false; - E(q1, q) = false; - } - } - A(*first_connected_free, q) = true; - B(*first_connected_free) = (P(q) % 4 == 3); - - // Set phase of q to -pi/2 - P(q) = 3; - - // No global phase change - - // Make A row reduced - for (unsigned q1 = 0; q1 < *first_connected_free; ++q1) { - if (A(q1, *first_connected_free)) { - for (unsigned q2 = *first_connected_free; q2 < n_qubits; ++q2) { - A(q1, q2) ^= A(*first_connected_free, q2); - } - B(q1) ^= B(*first_connected_free); - } - } - - // Need to apply rules for S and CZ gates on FCF to remove - unsigned s_gates = P(*first_connected_free); - P(*first_connected_free) = 0; - std::list cz_targets; - E(*first_connected_free, q) = false; - E(q, *first_connected_free) = false; - for (unsigned q1 = 0; q1 < n_qubits; ++q1) { - if (E(*first_connected_free, q1)) { - cz_targets.push_back(q1); - E(*first_connected_free, q1) = false; - E(q1, *first_connected_free) = false; - } - } - for (unsigned i = 0; i < s_gates; ++i) apply_S(*first_connected_free); - for (const unsigned trgt : cz_targets) - apply_CZ(*first_connected_free, trgt); - } else { - // q has phase +-pi/2, no connected leader, and no previous free, so q - // will become the new leader - A(q, q) = true; - - std::list connected; - for (unsigned q1 = q + 1; q1 < n_qubits; ++q1) { - if (E(q, q1)) { - // Local complementation - for (const unsigned q2 : connected) { - E(q1, q2) ^= true; - E(q2, q1) ^= true; - } - connected.push_back(q1); - - // Add to A - A(q, q1) = true; - - // Reset E rows - E(q, q1) = false; - E(q1, q) = false; - } - } - - // Update local and global phases - if (P(q) % 4 == 1) { - for (const unsigned q1 : connected) P(q1) += 3; - } else { - for (const unsigned q1 : connected) P(q1) += 1; - B(q) = true; - phase += -.5; - } - P(q) = 0; - } - } - } -} - -void APState::apply_gate(OpType type, const std::vector& qbs) { - switch (type) { - case OpType::Z: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - break; - } - case OpType::X: { - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - phase += .5; - break; - } - case OpType::Y: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - phase += 1.; - break; - } - case OpType::S: { - apply_S(qbs.at(0)); - break; - } - case OpType::Sdg: { - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - apply_S(qbs.at(0)); - break; - } - case OpType::V: { - apply_V(qbs.at(0)); - break; - } - case OpType::SX: { - apply_V(qbs.at(0)); - phase += .25; - break; - } - case OpType::Vdg: { - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - phase += 1.; - break; - } - case OpType::SXdg: { - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - phase += .75; - break; - } - case OpType::H: { - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - break; - } - case OpType::CX: { - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - break; - } - case OpType::CY: { - apply_V(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_V(qbs.at(1)); - apply_V(qbs.at(1)); - apply_V(qbs.at(1)); - phase += 1.; - break; - } - case OpType::CZ: { - apply_CZ(qbs.at(0), qbs.at(1)); - break; - } - case OpType::ZZMax: { - apply_S(qbs.at(0)); - apply_S(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - phase -= .25; - break; - } - case OpType::ECR: { - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - phase += .25; - break; - } - case OpType::ISWAPMax: { - apply_V(qbs.at(0)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_V(qbs.at(0)); - apply_V(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_V(qbs.at(0)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - phase += 1.; - break; - } - case OpType::SWAP: { - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - apply_S(qbs.at(1)); - apply_V(qbs.at(1)); - apply_S(qbs.at(1)); - apply_CZ(qbs.at(0), qbs.at(1)); - apply_S(qbs.at(0)); - apply_V(qbs.at(0)); - apply_S(qbs.at(0)); - break; - } - case OpType::BRIDGE: { - apply_S(qbs.at(2)); - apply_V(qbs.at(2)); - apply_S(qbs.at(2)); - apply_CZ(qbs.at(0), qbs.at(2)); - apply_S(qbs.at(2)); - apply_V(qbs.at(2)); - apply_S(qbs.at(2)); - break; - } - case OpType::noop: { - break; - } - case OpType::Phase: { - throw std::logic_error( - "OpType::Phase cannot be applied via APState::apply_gate"); - } - default: { - throw BadOpType( - "Cannot be applied to a APState: not a Clifford gate", type); - } - } -} - -} // namespace tket \ No newline at end of file diff --git a/tket/test/src/test_APState.cpp b/tket/test/src/test_APState.cpp deleted file mode 100644 index 12320e8a90..0000000000 --- a/tket/test/src/test_APState.cpp +++ /dev/null @@ -1,347 +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 - -#include "Simulation/ComparisonFunctions.hpp" -#include "testutil.hpp" -#include "tket/Circuit/Simulation/CircuitSimulator.hpp" -#include "tket/Clifford/APState.hpp" -#include "tket/Converters/Converters.hpp" - -namespace tket { -namespace test_APState { - -void test_apply_gate( - const MatrixXb& A, const VectorXb& B, const MatrixXb& E, - const Eigen::VectorXi& P, OpType ot, const std::vector& args) { - APState ap(A, B, E, P, 0); - ap.verify(); - auto sv_before = ap.to_statevector(); - - Circuit circ(A.rows()); - circ.add_op(ot, args); - auto gate_u = tket_sim::get_unitary(circ); - - ap.apply_gate(ot, args); - - ap.verify(); - auto sv_after = ap.to_statevector(); - - CHECK(tket_sim::compare_statevectors_or_unitaries( - gate_u * sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL)); -} - -SCENARIO("CZ cases") { - GIVEN("CZ on free qubits") { - MatrixXb A = MatrixXb::Zero(4, 4); - VectorXb B = VectorXb::Zero(4); - MatrixXb E = MatrixXb::Zero(4, 4); - Eigen::VectorXi P = Eigen::VectorXi::Zero(4); - A(0, 0) = A(0, 1) = A(0, 3) = true; - E(2, 3) = E(3, 2) = true; - test_apply_gate(A, B, E, P, OpType::CZ, {1, 2}); - } - GIVEN("CZ on one leading qubit and connected free") { - for (unsigned b = 0; b < 2; ++b) { - MatrixXb A = MatrixXb::Zero(5, 5); - VectorXb B = VectorXb::Zero(5); - MatrixXb E = MatrixXb::Zero(5, 5); - Eigen::VectorXi P = Eigen::VectorXi::Zero(5); - A(0, 0) = A(0, 2) = A(0, 3) = true; - A(1, 1) = A(1, 3) = A(1, 4) = true; - B(1) = (b == 1); - test_apply_gate(A, B, E, P, OpType::CZ, {0, 3}); - } - } - GIVEN("CZ on one leading qubit and unconnected free") { - for (unsigned b = 0; b < 2; ++b) { - MatrixXb A = MatrixXb::Zero(5, 5); - VectorXb B = VectorXb::Zero(5); - MatrixXb E = MatrixXb::Zero(5, 5); - Eigen::VectorXi P = Eigen::VectorXi::Zero(5); - A(0, 0) = A(0, 2) = A(0, 3) = true; - A(1, 1) = A(1, 3) = A(1, 4) = true; - B(1) = (b == 1); - test_apply_gate(A, B, E, P, OpType::CZ, {0, 4}); - } - } - GIVEN("CZ on leading qubits") { - for (unsigned b1 = 0; b1 < 2; ++b1) { - for (unsigned b2 = 0; b2 < 2; ++b2) { - MatrixXb A = MatrixXb::Zero(8, 8); - VectorXb B = VectorXb::Zero(8); - MatrixXb E = MatrixXb::Zero(8, 8); - Eigen::VectorXi P = Eigen::VectorXi::Zero(8); - A(0, 0) = A(0, 2) = A(0, 3) = A(0, 4) = A(0, 5) = true; - A(1, 1) = A(1, 4) = A(1, 5) = A(1, 6) = A(1, 7) = true; - B(0) = (b1 == 1); - B(1) = (b2 == 1); - test_apply_gate(A, B, E, P, OpType::CZ, {0, 1}); - } - } - } -} - -SCENARIO("S cases") { - GIVEN("S on free qubit") { - MatrixXb A = MatrixXb::Zero(3, 3); - VectorXb B = VectorXb::Zero(3); - MatrixXb E = MatrixXb::Zero(3, 3); - Eigen::VectorXi P = Eigen::VectorXi::Zero(3); - A(0, 0) = A(0, 1) = A(0, 2) = true; - E(1, 2) = E(2, 1) = true; - test_apply_gate(A, B, E, P, OpType::S, {2}); - } - GIVEN("S on leading qubit") { - for (unsigned b = 0; b < 2; ++b) { - MatrixXb A = MatrixXb::Zero(4, 4); - VectorXb B = VectorXb::Zero(4); - MatrixXb E = MatrixXb::Zero(4, 4); - Eigen::VectorXi P = Eigen::VectorXi::Zero(4); - A(0, 0) = A(0, 1) = A(0, 2) = true; - B(0) = (b == 1); - E(1, 3) = E(3, 1) = true; - test_apply_gate(A, B, E, P, OpType::S, {0}); - } - } - GIVEN("S on disconnected leading qubit") { - for (unsigned b = 0; b < 2; ++b) { - MatrixXb A = MatrixXb::Zero(1, 1); - VectorXb B = VectorXb::Zero(1); - MatrixXb E = MatrixXb::Zero(1, 1); - Eigen::VectorXi P = Eigen::VectorXi::Zero(1); - A(0, 0) = true; - B(0) = (b == 1); - test_apply_gate(A, B, E, P, OpType::S, {0}); - } - } -} - -SCENARIO("V cases") { - GIVEN("V on leading qubit") { - for (unsigned b = 0; b < 2; ++b) { - MatrixXb A = MatrixXb::Zero(4, 4); - VectorXb B = VectorXb::Zero(4); - MatrixXb E = MatrixXb::Zero(4, 4); - Eigen::VectorXi P = Eigen::VectorXi::Zero(4); - A(0, 0) = A(0, 2) = A(0, 3) = true; - A(1, 1) = A(1, 3) = true; - B(0) = (b == 1); - test_apply_gate(A, B, E, P, OpType::V, {0}); - } - } - GIVEN("V on free qubit with some leading") { - for (unsigned b = 0; b < 2; ++b) { - for (unsigned p = 0; p < 4; ++p) { - MatrixXb A = MatrixXb::Zero(9, 9); - VectorXb B = VectorXb::Zero(9); - MatrixXb E = MatrixXb::Zero(9, 9); - Eigen::VectorXi P = Eigen::VectorXi::Zero(9); - A(0, 0) = A(0, 2) = A(0, 4) = A(0, 5) = A(0, 7) = true; - A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = A(1, 5) = A(1, 6) = true; - B(1) = (b == 1); - E(4, 5) = E(5, 4) = true; - E(4, 6) = E(6, 4) = true; - E(4, 7) = E(7, 4) = true; - E(4, 8) = E(8, 4) = true; - P(4) = p; - test_apply_gate(A, B, E, P, OpType::V, {4}); - } - } - } - GIVEN("V on free qubit with some earlier connected free") { - for (unsigned p1 = 0; p1 < 4; ++p1) { - for (unsigned p2 = 0; p2 < 4; ++p2) { - MatrixXb A = MatrixXb::Zero(9, 9); - VectorXb B = VectorXb::Zero(9); - MatrixXb E = MatrixXb::Zero(9, 9); - Eigen::VectorXi P = Eigen::VectorXi::Zero(9); - A(0, 0) = true; - A(1, 1) = A(1, 4) = A(1, 6) = true; - A(2, 2) = A(2, 4) = A(2, 7) = true; - A(3, 3) = A(3, 4) = A(3, 8) = A(3, 8) = true; - E(4, 5) = E(5, 4) = true; - E(4, 7) = E(7, 4) = true; - E(4, 8) = E(8, 4) = true; - E(5, 6) = E(6, 5) = true; - E(5, 7) = E(7, 5) = true; - E(5, 8) = E(8, 5) = true; - P(4) = p1; - P(5) = p2; - test_apply_gate(A, B, E, P, OpType::V, {5}); - } - } - } - GIVEN("V on free qubit with no earlier connected free") { - for (unsigned p = 0; p < 4; ++p) { - MatrixXb A = MatrixXb::Zero(4, 4); - VectorXb B = VectorXb::Zero(4); - MatrixXb E = MatrixXb::Zero(4, 4); - Eigen::VectorXi P = Eigen::VectorXi::Zero(4); - A(0, 0) = A(0, 2) = true; - E(1, 2) = E(2, 1) = true; - E(1, 3) = E(3, 1) = true; - P(1) = p; - test_apply_gate(A, B, E, P, OpType::V, {1}); - } - } - GIVEN("V on disconnected free qubit") { - for (unsigned p = 0; p < 4; ++p) { - MatrixXb A = MatrixXb::Zero(1, 1); - VectorXb B = VectorXb::Zero(1); - MatrixXb E = MatrixXb::Zero(1, 1); - Eigen::VectorXi P = Eigen::VectorXi::Zero(1, 1); - P(0) = p; - test_apply_gate(A, B, E, P, OpType::V, {0}); - } - } -} - -SCENARIO("Gate encodings") { - std::list>> test_gates = { - {OpType::Z, {0}}, {OpType::X, {0}}, - {OpType::Y, {0}}, {OpType::S, {0}}, - {OpType::Sdg, {0}}, {OpType::V, {0}}, - {OpType::Vdg, {0}}, {OpType::SX, {0}}, - {OpType::SXdg, {0}}, {OpType::H, {0}}, - {OpType::CX, {0, 1}}, {OpType::CY, {0, 1}}, - {OpType::CZ, {0, 1}}, {OpType::ZZMax, {0, 1}}, - {OpType::ECR, {0, 1}}, {OpType::ISWAPMax, {0, 1}}, - {OpType::SWAP, {0, 1}}, {OpType::BRIDGE, {0, 1, 2}}, - {OpType::noop, {0}}, - }; - GIVEN("Check Z actions") { - for (const auto& com : test_gates) { - MatrixXb A = MatrixXb::Identity(3, 3); - VectorXb B = VectorXb::Zero(3); - MatrixXb E = MatrixXb::Zero(3, 3); - Eigen::VectorXi P = Eigen::VectorXi::Zero(3); - test_apply_gate(A, B, E, P, com.first, com.second); - } - } - GIVEN("Check X actions") { - for (const auto& com : test_gates) { - MatrixXb A = MatrixXb::Zero(3, 3); - VectorXb B = VectorXb::Zero(3); - MatrixXb E = MatrixXb::Zero(3, 3); - Eigen::VectorXi P = Eigen::VectorXi::Zero(3); - test_apply_gate(A, B, E, P, com.first, com.second); - } - } -} - -SCENARIO("Loading from a statevector") { - MatrixXb A = MatrixXb::Zero(4, 4); - VectorXb B = VectorXb::Zero(4); - MatrixXb E = MatrixXb::Zero(4, 4); - Eigen::VectorXi P = Eigen::VectorXi::Zero(4); - A(0, 0) = A(0, 2) = A(0, 3) = true; - A(1, 1) = A(1, 2) = true; - B(0) = true; - E(2, 3) = E(3, 2) = true; - P(2) = 1; - P(3) = 2; - APState ap(A, B, E, P, 0); - auto sv = ap.to_statevector(); - APState reconstructed(sv); - auto sv2 = reconstructed.to_statevector(); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv, sv2, tket_sim::MatrixEquivalence::EQUAL)); -} - -SCENARIO("Converting from/to a circuit") { - GIVEN("A circuit in the standard AP form") { - Circuit circ(4); - circ.qubit_create_all(); - circ.add_op(OpType::H, {2}); - circ.add_op(OpType::H, {3}); - circ.add_op(OpType::CX, {2, 0}); - circ.add_op(OpType::CX, {2, 1}); - circ.add_op(OpType::CX, {3, 1}); - circ.add_op(OpType::S, {2}); - circ.add_op(OpType::Z, {3}); - APState ap = circuit_to_apstate(circ); - auto sv_circ = tket_sim::get_statevector(circ); - auto sv_ap = ap.to_statevector(); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); - Circuit reconstructed = apstate_to_circuit(ap); - CHECK(circ == reconstructed); - } - GIVEN("A generic circuit") { - Circuit circ(4); - circ.qubit_create_all(); - circ.add_op(OpType::V, {0}); - circ.add_op(OpType::CX, {0, 1}); - circ.add_op(OpType::CY, {1, 3}); - circ.add_op(OpType::H, {3}); - circ.add_op(OpType::ZZMax, {2, 3}); - APState ap = circuit_to_apstate(circ); - auto sv_circ = tket_sim::get_statevector(circ); - auto sv_ap = ap.to_statevector(); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); - Circuit reconstructed = apstate_to_circuit(ap); - auto sv_rec = tket_sim::get_statevector(reconstructed); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv_circ, sv_rec, tket_sim::MatrixEquivalence::EQUAL)); - } -} - -SCENARIO("Converting from/to a tableau") { - GIVEN("Check up to global phase using circuit") { - Circuit circ(8); - circ.qubit_create_all(); - circ.add_op(OpType::X, {1}); - circ.add_op(OpType::X, {5}); - circ.add_op(OpType::H, {2}); - circ.add_op(OpType::H, {4}); - circ.add_op(OpType::H, {6}); - circ.add_op(OpType::H, {7}); - circ.add_op(OpType::CX, {2, 1}); - circ.add_op(OpType::CX, {4, 0}); - circ.add_op(OpType::CX, {4, 3}); - circ.add_op(OpType::CX, {6, 0}); - circ.add_op(OpType::CX, {6, 1}); - circ.add_op(OpType::CX, {7, 5}); - circ.add_op(OpType::CZ, {2, 6}); - circ.add_op(OpType::CZ, {4, 6}); - circ.add_op(OpType::CZ, {4, 7}); - circ.add_op(OpType::CZ, {6, 7}); - circ.add_op(OpType::S, {2}); - circ.add_op(OpType::Sdg, {4}); - circ.add_op(OpType::Z, {7}); - ChoiMixTableau cmt = circuit_to_cm_tableau(circ); - APState ap = tableau_to_apstate(cmt.tab_); - auto sv_circ = tket_sim::get_statevector(circ); - auto sv_ap = ap.to_statevector(); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL)); - SymplecticTableau tab2 = apstate_to_tableau(ap); - ChoiMixTableau cmt2(tab2.xmat, tab2.zmat, tab2.phase); - auto [circ2, perm] = cm_tableau_to_exact_circuit(cmt2); - qubit_map_t inv; - for (const std::pair& qp : perm) - inv.insert({qp.second, qp.first}); - circ2.permute_boundary_output(inv); - auto sv_circ2 = tket_sim::get_statevector(circ2); - CHECK(tket_sim::compare_statevectors_or_unitaries( - sv_circ, sv_circ2, - tket_sim::MatrixEquivalence::EQUAL_UP_TO_GLOBAL_PHASE)); - } -} - -} // namespace test_APState -} // namespace tket \ No newline at end of file From c8bf6b3c9c53824cc35098163b660afc0b0578be Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 14:52:02 +0000 Subject: [PATCH 10/19] Missed one --- tket/src/Converters/APStateConverters.cpp | 176 ---------------------- 1 file changed, 176 deletions(-) delete mode 100644 tket/src/Converters/APStateConverters.cpp diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp deleted file mode 100644 index d73dcbc53d..0000000000 --- a/tket/src/Converters/APStateConverters.cpp +++ /dev/null @@ -1,176 +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/Converters.hpp" - -namespace tket { - -APState circuit_to_apstate(const Circuit& circ) { - APState aps(circ.n_qubits()); - std::map qb_ordering; - for (const Qubit& q : circ.all_qubits()) - qb_ordering.insert({q, (unsigned)qb_ordering.size()}); - for (const Command& com : circ) { - auto args = com.get_args(); - std::vector qbs; - for (const UnitID& q : args) qbs.push_back(qb_ordering.at(q)); - aps.apply_gate(com.get_op_ptr()->get_type(), qbs); - } - return aps; -} - -Circuit apstate_to_circuit(const APState& ap) { - unsigned n_qbs = ap.A.rows(); - Circuit circ(n_qbs); - circ.qubit_create_all(); - for (unsigned q = 0; q < n_qbs; ++q) { - if (!ap.A(q, q)) - circ.add_op(OpType::H, {q}); - else if (ap.B(q)) - circ.add_op(OpType::X, {q}); - } - for (unsigned trgt = 0; trgt < n_qbs; ++trgt) { - if (ap.A(trgt, trgt)) { - for (unsigned ctrl = trgt + 1; ctrl < n_qbs; ++ctrl) - if (ap.A(trgt, ctrl)) circ.add_op(OpType::CX, {ctrl, trgt}); - } - } - for (unsigned q1 = 0; q1 < n_qbs; ++q1) { - if (ap.A(q1, q1)) continue; - for (unsigned q2 = q1 + 1; q2 < n_qbs; ++q2) { - if (ap.E(q1, q2)) circ.add_op(OpType::CZ, {q1, q2}); - } - switch (ap.P(q1) % 4) { - case 1: { - circ.add_op(OpType::S, {q1}); - break; - } - case 2: { - circ.add_op(OpType::Z, {q1}); - break; - } - case 3: { - circ.add_op(OpType::Sdg, {q1}); - break; - } - default: { - break; - } - } - } - circ.add_phase(ap.phase); - return circ; -} - -APState tableau_to_apstate(SymplecticTableau tab) { - unsigned n_qbs = tab.get_n_qubits(); - if (tab.get_n_rows() != n_qbs) - throw std::logic_error( - "tableau_to_apstate requires a tableau with n commuting rows for n " - "qubits"); - MatrixXb fullmat = MatrixXb::Zero(n_qbs, 2 * n_qbs); - /** - * Gaussian elimination by the x matrix first ensures the bottom rows are only - * Zs, i.e. describing rows of A. Reversing the columns of the x matrix - * guarantees that each row has an X on at most one free qubit, simplifying - * the code for finding E and P - */ - for (unsigned c = 0; c < n_qbs; ++c) { - fullmat.col(c) = tab.xmat.col(n_qbs - 1 - c); - } - fullmat.block(0, n_qbs, n_qbs, n_qbs) = tab.zmat; - std::vector> row_ops = - gaussian_elimination_row_ops(fullmat); - for (const std::pair& op : row_ops) - tab.row_mult(op.first, op.second); - - MatrixXb A = MatrixXb::Zero(n_qbs, n_qbs); - VectorXb B = VectorXb::Zero(n_qbs); - MatrixXb E = MatrixXb::Zero(n_qbs, n_qbs); - Eigen::VectorXi P = Eigen::VectorXi::Zero(n_qbs); - - unsigned n_leading = 0; - for (unsigned r = n_qbs; r > 0;) { - --r; - bool is_a_row = true; - for (unsigned c = 0; c < n_qbs; ++c) { - if (tab.xmat(r, c)) { - is_a_row = false; - break; - } - } - if (!is_a_row) break; - unsigned first = 0; - for (unsigned c = 0; c < n_qbs; ++c) { - if (tab.zmat(r, c)) { - first = c; - break; - } - } - A.row(first) = tab.zmat.row(r); - B(first) = tab.phase(r); - ++n_leading; - } - /** - * Each free qubit is after all leaders connected to it, by reduced - * row-echelon of A. Therefore Gaussian elimination of the x matrix in reverse - * order gives rows corresponding to the columns of A of free qubits (plus the - * free qubit itself). - * - * Then the corresponding row of the z matrix is exactly the free qubit's row - * of E, from which we take the diagonal and phase vector to determine P. - */ - for (unsigned r = 0; r < n_qbs - n_leading; ++r) { - unsigned free = 0; - for (unsigned c = n_qbs; c > 0;) { - --c; - if (tab.xmat(r, c)) { - free = c; - break; - } - } - E.row(free) = tab.zmat.row(r); - if (tab.zmat(r, free)) P(free) += 1; - if (tab.phase(r)) P(free) += 2; - E(free, free) = false; - } - - return APState(A, B, E, P, 0); -} - -SymplecticTableau apstate_to_tableau(const APState& ap) { - unsigned n_qbs = ap.A.rows(); - MatrixXb xmat = MatrixXb::Zero(n_qbs, n_qbs); - MatrixXb zmat = MatrixXb::Zero(n_qbs, n_qbs); - VectorXb phase = VectorXb::Zero(n_qbs); - - for (unsigned q = 0; q < n_qbs; ++q) { - if (ap.A(q, q)) { - // Leading qubit - zmat.row(q) = ap.A.row(q); - phase(q) = ap.B(q); - } else { - // Free qubit - xmat.row(q) = ap.A.col(q); - xmat(q, q) = true; - zmat.row(q) = ap.E.row(q); - zmat(q, q) = (ap.P(q) % 2 == 1); - phase(q) = (ap.P(q) % 4 > 1); - } - } - - return SymplecticTableau(xmat, zmat, phase); -} - -} // namespace tket \ No newline at end of file From 37f8e2d9c9f2e95785c4d47268fe95470c03cd4b Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 14:52:53 +0000 Subject: [PATCH 11/19] Rename files --- tket/include/tket/Clifford/{APState2.hpp => APState.hpp} | 0 tket/src/Clifford/{APState2.cpp => APState.cpp} | 0 .../Converters/{APStateConverters2.cpp => APStateConverters.cpp} | 0 tket/test/src/{test_APState2.cpp => test_APState.cpp} | 0 4 files changed, 0 insertions(+), 0 deletions(-) rename tket/include/tket/Clifford/{APState2.hpp => APState.hpp} (100%) rename tket/src/Clifford/{APState2.cpp => APState.cpp} (100%) rename tket/src/Converters/{APStateConverters2.cpp => APStateConverters.cpp} (100%) rename tket/test/src/{test_APState2.cpp => test_APState.cpp} (100%) diff --git a/tket/include/tket/Clifford/APState2.hpp b/tket/include/tket/Clifford/APState.hpp similarity index 100% rename from tket/include/tket/Clifford/APState2.hpp rename to tket/include/tket/Clifford/APState.hpp diff --git a/tket/src/Clifford/APState2.cpp b/tket/src/Clifford/APState.cpp similarity index 100% rename from tket/src/Clifford/APState2.cpp rename to tket/src/Clifford/APState.cpp diff --git a/tket/src/Converters/APStateConverters2.cpp b/tket/src/Converters/APStateConverters.cpp similarity index 100% rename from tket/src/Converters/APStateConverters2.cpp rename to tket/src/Converters/APStateConverters.cpp diff --git a/tket/test/src/test_APState2.cpp b/tket/test/src/test_APState.cpp similarity index 100% rename from tket/test/src/test_APState2.cpp rename to tket/test/src/test_APState.cpp From da6cf2d64fe5beb0016a32797a3dbc4441dee160 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 14:59:56 +0000 Subject: [PATCH 12/19] Change file references back --- tket/CMakeLists.txt | 5 ++--- tket/include/tket/Converters/Converters.hpp | 2 +- tket/src/Clifford/APState.cpp | 2 +- tket/test/CMakeLists.txt | 2 +- tket/test/src/test_APState.cpp | 6 +++--- 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt index 4656592384..650a1d61c6 100644 --- a/tket/CMakeLists.txt +++ b/tket/CMakeLists.txt @@ -168,7 +168,7 @@ target_sources(tket src/Circuit/SubcircuitFinder.cpp src/Circuit/ThreeQubitConversion.cpp src/Circuit/ToffoliBox.cpp - src/Clifford/APState2.cpp + src/Clifford/APState.cpp src/Clifford/ChoiMixTableau.cpp src/Clifford/SymplecticTableau.cpp src/Clifford/UnitaryTableau.cpp @@ -324,8 +324,7 @@ target_sources(tket include/tket/Circuit/StatePreparation.hpp include/tket/Circuit/ThreeQubitConversion.hpp include/tket/Circuit/ToffoliBox.hpp - # include/tket/Clifford/APState.hpp - include/tket/Clifford/APState2.hpp + include/tket/Clifford/APState.hpp include/tket/Clifford/ChoiMixTableau.hpp include/tket/Clifford/SymplecticTableau.hpp include/tket/Clifford/UnitaryTableau.hpp diff --git a/tket/include/tket/Converters/Converters.hpp b/tket/include/tket/Converters/Converters.hpp index 5061204781..2ccdd3cc31 100644 --- a/tket/include/tket/Converters/Converters.hpp +++ b/tket/include/tket/Converters/Converters.hpp @@ -15,7 +15,7 @@ #pragma once #include "tket/Circuit/Circuit.hpp" -#include "tket/Clifford/APState2.hpp" +#include "tket/Clifford/APState.hpp" #include "tket/Clifford/ChoiMixTableau.hpp" #include "tket/Clifford/UnitaryTableau.hpp" #include "tket/PauliGraph/PauliGraph.hpp" diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp index 6fae7249ba..de582405a4 100644 --- a/tket/src/Clifford/APState.cpp +++ b/tket/src/Clifford/APState.cpp @@ -12,7 +12,7 @@ // See the License for the specific language governing permissions and // limitations under the License. -#include "tket/Clifford/APState2.hpp" +#include "tket/Clifford/APState.hpp" #include #include diff --git a/tket/test/CMakeLists.txt b/tket/test/CMakeLists.txt index 5e2d38c21c..fe7011d6c5 100644 --- a/tket/test/CMakeLists.txt +++ b/tket/test/CMakeLists.txt @@ -102,7 +102,7 @@ add_executable(test-tket src/Simulation/test_CircuitSimulator.cpp src/Simulation/test_PauliExpBoxUnitaryCalculator.cpp src/test_AASRoute.cpp - src/test_APState2.cpp + src/test_APState.cpp src/test_ArchitectureAwareSynthesis.cpp src/test_Architectures.cpp src/test_Assertion.cpp diff --git a/tket/test/src/test_APState.cpp b/tket/test/src/test_APState.cpp index dcaeceeb01..6d94602994 100644 --- a/tket/test/src/test_APState.cpp +++ b/tket/test/src/test_APState.cpp @@ -17,11 +17,11 @@ #include "Simulation/ComparisonFunctions.hpp" #include "testutil.hpp" #include "tket/Circuit/Simulation/CircuitSimulator.hpp" -#include "tket/Clifford/APState2.hpp" +#include "tket/Clifford/APState.hpp" #include "tket/Converters/Converters.hpp" namespace tket { -namespace test_APState2 { +namespace test_APState { void test_apply_gate( const MatrixXb& A, const VectorXb& B, const MatrixXb& E, @@ -510,5 +510,5 @@ SCENARIO("Converting from/to a tableau") { } } -} // namespace test_APState2 +} // namespace test_APState } // namespace tket \ No newline at end of file From e83c0a675997efe01e6855b75d656f0e32a1a310 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 15:00:11 +0000 Subject: [PATCH 13/19] The last commit missed this one as well --- tket/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt index 650a1d61c6..e7cfc3a681 100644 --- a/tket/CMakeLists.txt +++ b/tket/CMakeLists.txt @@ -172,7 +172,7 @@ target_sources(tket src/Clifford/ChoiMixTableau.cpp src/Clifford/SymplecticTableau.cpp src/Clifford/UnitaryTableau.cpp - src/Converters/APStateConverters2.cpp + src/Converters/APStateConverters.cpp src/Converters/ChoiMixTableauConverters.cpp src/Converters/Gauss.cpp src/Converters/PauliGraphConverters.cpp From 3d2ab89f5301cc6f715cda4a2a0d1e538be8b0f3 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 15:01:17 +0000 Subject: [PATCH 14/19] Bump tket version --- pytket/conanfile.py | 2 +- tket/conanfile.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pytket/conanfile.py b/pytket/conanfile.py index 0bc83653b6..c50c6c9909 100644 --- a/pytket/conanfile.py +++ b/pytket/conanfile.py @@ -38,7 +38,7 @@ def requirements(self): self.requires("pybind11_json/0.2.14") self.requires("symengine/0.13.0") self.requires("tkassert/0.3.4@tket/stable") - self.requires("tket/1.3.52@tket/stable") + self.requires("tket/1.3.53@tket/stable") self.requires("tklog/0.3.3@tket/stable") self.requires("tkrng/0.3.3@tket/stable") self.requires("tktokenswap/0.3.9@tket/stable") diff --git a/tket/conanfile.py b/tket/conanfile.py index 60aceeff08..2602d319f1 100644 --- a/tket/conanfile.py +++ b/tket/conanfile.py @@ -23,7 +23,7 @@ class TketConan(ConanFile): name = "tket" - version = "1.3.52" + version = "1.3.53" package_type = "library" license = "Apache 2" homepage = "https://github.com/CQCL/tket" From fb097934a6e6019ea1694625d6ab1e6c6ae6fa8c Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 15:46:39 +0000 Subject: [PATCH 15/19] Doxygen math blocks --- tket/include/tket/Clifford/APState.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp index 3a3460b9bb..3e49111205 100644 --- a/tket/include/tket/Clifford/APState.hpp +++ b/tket/include/tket/Clifford/APState.hpp @@ -32,7 +32,7 @@ namespace tket { * - A n vector P of integers mod 4 describing S gates. * * This gives a canonical statevector (up to a normalisation scalar): - * \sum_{x, Ax=B} i^{\Phi(x)} |x> + * \f[\sum_{x, Ax=B} i^{\Phi(x)} |x>\f] * * This canonical statevector fixes a reference phase from which we can track * global phase with an additional parameter. @@ -40,8 +40,8 @@ namespace tket { * We can generalise this to mixed qubit states by adding a binary (n,n) matrix * C, with which we now have the canonical density matrix (up to a normalisation * scalar): - * \sum_{x_1, x_2; Ax_1 = B = Ax_2, Cx_1 = Cx_2} - * i^{x_1^T E x_1 + P^T x_1} |x_1> Date: Thu, 28 Nov 2024 15:46:56 +0000 Subject: [PATCH 16/19] Wrong size of matrix block --- tket/src/Clifford/APState.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp index de582405a4..1878e78b4e 100644 --- a/tket/src/Clifford/APState.cpp +++ b/tket/src/Clifford/APState.cpp @@ -217,7 +217,7 @@ APState::APState(const Eigen::MatrixXcd& dm) { A.block(0, 0, first_zero - first_right, n_qbs) = full_mat.block(first_right, n_qbs, first_zero - first_right, n_qbs); B.head(first_zero - first_right) = - full_mat.block(first_right, 2 * n_qbs, first_zero - first_right, n_qbs); + full_mat.block(first_right, 2 * n_qbs, first_zero - first_right, 1); // Flip the column order and reduce the remaining rows to obtain CC above AI; // get the row combinations from the reordered matrix but apply them to the From b3b0b6a986e4e912309c7a1c3b9b75fe7b101534 Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 28 Nov 2024 17:12:00 +0000 Subject: [PATCH 17/19] Fixed another bounds issue --- tket/src/Converters/APStateConverters.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp index 54cc01c2d7..36c1cf7db6 100644 --- a/tket/src/Converters/APStateConverters.cpp +++ b/tket/src/Converters/APStateConverters.cpp @@ -237,7 +237,7 @@ APState tableau_to_apstate(SymplecticTableau tab) { // Start by identifying the free and mixed qubits std::map free_to_row; - for (unsigned r = 0; r < n_qbs - n_leading; ++r) { + for (unsigned r = 0; r < n_rows - n_leading; ++r) { for (unsigned c = n_qbs - r; c > 0;) { --c; if (tab.xmat(r, c)) { From 981515536fc2f9813ca1099ded0b55c89c5f1b5a Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Wed, 4 Dec 2024 17:28:00 +0000 Subject: [PATCH 18/19] Increase coverage and fix some extra bugs --- tket/src/Clifford/APState.cpp | 11 +- tket/src/Converters/APStateConverters.cpp | 4 + tket/test/src/test_APState.cpp | 330 ++++++++++++++++++++++ 3 files changed, 341 insertions(+), 4 deletions(-) diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp index 1878e78b4e..c9d96d241b 100644 --- a/tket/src/Clifford/APState.cpp +++ b/tket/src/Clifford/APState.cpp @@ -667,7 +667,7 @@ void APState::apply_V(unsigned q) { // No local phase change } // Local phase on q - P(q) += 3; + P(q) = 3; // No global phase change } else { VectorXb new_row = VectorXb::Zero(n_qbs); @@ -1346,9 +1346,12 @@ void ChoiAPState::post_select(const Qubit& qb, TableauSegment seg) { unsigned removed_u_qb = ap_.post_select(u_qb); col_index_.left.erase(lit); tableau_col_index_t::right_iterator rit = col_index_.right.find(removed_u_qb); - col_key_t removed_key = rit->second; - col_index_.right.erase(rit); - col_index_.insert({removed_key, u_qb}); + if (rit != col_index_.right.end()) { + // Ignore if post-selected last column so none others need changing + col_key_t removed_key = rit->second; + col_index_.right.erase(rit); + col_index_.insert({removed_key, u_qb}); + } } void ChoiAPState::discard_qubit(const Qubit& qb, TableauSegment seg) { diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp index 36c1cf7db6..d8070ff35c 100644 --- a/tket/src/Converters/APStateConverters.cpp +++ b/tket/src/Converters/APStateConverters.cpp @@ -158,6 +158,10 @@ std::pair choi_apstate_to_unitary_extension_circuit( std::pair res = cm_tableau_to_unitary_extension_circuit( tab, init_names, post_names, cx_config); ChoiAPState reconstructed = circuit_to_choi_apstate(res.first); + for (const Qubit& q : init_names) + reconstructed.post_select(q, ChoiAPState::TableauSegment::Input); + for (const Qubit& q : post_names) + reconstructed.post_select(q, ChoiAPState::TableauSegment::Output); reconstructed.rename_qubits(res.second, ChoiAPState::TableauSegment::Output); ap.canonical_column_order(); ap.normal_form(); diff --git a/tket/test/src/test_APState.cpp b/tket/test/src/test_APState.cpp index 6d94602994..064c3830df 100644 --- a/tket/test/src/test_APState.cpp +++ b/tket/test/src/test_APState.cpp @@ -63,6 +63,124 @@ void test_apply_gate_dm( CHECK(dm_after.isApprox(gate_u * dm_before * gate_u.adjoint(), EPS)); } +void test_apply_gate_dm_input( + const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E, + const Eigen::VectorXi& P, OpType ot, const qubit_vector_t& args) { + ChoiAPState ap(A, B, C, E, P, 0, A.cols()); + ap.ap_.verify(); + auto dm_before = ap.ap_.to_density_matrix(); + + Circuit circ(A.cols()); + circ.add_op(ot, args); + auto gate_u = tket_sim::get_unitary(circ); + + ap.apply_gate(ot, args, ChoiAPState::TableauSegment::Input); + + ap.ap_.verify(); + auto dm_after = ap.ap_.to_density_matrix(); + + CHECK(dm_after.isApprox( + gate_u.transpose() * dm_before * gate_u.conjugate(), EPS)); +} + +SCENARIO("Normal form") { + GIVEN("Make A reduced row-echelon form (pure state)") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 2) = A(0, 3) = true; + A(1, 0) = A(1, 1) = A(1, 2) = true; + A(2, 0) = A(2, 1) = A(2, 2) = true; + A(3, 0) = A(3, 1) = true; + APState ap(A, B, C, E, P, 0); + auto sv_before = ap.to_statevector(); + ap.normal_form(); + auto sv_after = ap.to_statevector(); + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL)); + MatrixXb corrA = MatrixXb::Zero(4, 4); + corrA(0, 0) = corrA(0, 1) = true; + corrA(1, 2) = true; + corrA(2, 3) = true; + CHECK(ap.A == corrA); + } + GIVEN("Make A and C reduced row-echelon form") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 2) = A(0, 3) = true; + A(1, 0) = A(1, 1) = A(1, 2) = true; + C(0, 0) = C(0, 1) = C(0, 2) = true; + C(1, 0) = true; + C(3, 2) = true; + APState ap(A, B, C, E, P, 0); + auto dm_before = ap.to_density_matrix(); + ap.normal_form(); + auto dm_after = ap.to_density_matrix(); + CHECK(dm_after.isApprox(dm_before, EPS)); + MatrixXb corrA = MatrixXb::Zero(4, 4); + MatrixXb corrC = MatrixXb::Zero(4, 4); + corrA(0, 0) = corrA(0, 1) = corrA(0, 3) = true; + corrA(1, 2) = corrA(1, 3) = true; + corrC(0, 1) = true; + corrC(1, 3) = true; + CHECK(ap.A == corrA); + CHECK(ap.C == corrC); + } + GIVEN("Removing leaders from E and P (pure state)") { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb C = MatrixXb::Zero(5, 5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + A(0, 0) = A(0, 2) = A(0, 3) = true; + A(1, 1) = A(1, 2) = true; + E(0, 1) = E(1, 0) = true; + E(0, 3) = E(3, 0) = true; + E(0, 4) = E(4, 0) = true; + for (unsigned b = 0; b < 2; ++b) { + for (unsigned p = 0; p < 4; ++p) { + B(0) = (b == 1); + P(0) = p; + APState ap(A, B, C, E, P, 0); + auto sv_before = ap.to_statevector(); + ap.normal_form(); + auto sv_after = ap.to_statevector(); + // Just check using statevector; too much changes in each case to nicely + // test the matrices + CHECK(tket_sim::compare_statevectors_or_unitaries( + sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL)); + } + } + } + GIVEN("Removing mixed qubits from E and P") { + MatrixXb A = MatrixXb::Zero(5, 5); + VectorXb B = VectorXb::Zero(5); + MatrixXb C = MatrixXb::Zero(5, 5); + MatrixXb E = MatrixXb::Zero(5, 5); + Eigen::VectorXi P = Eigen::VectorXi::Zero(5); + C(0, 0) = C(0, 2) = C(0, 3) = true; + C(1, 1) = true; + E(0, 1) = E(1, 0) = true; + E(0, 2) = E(2, 0) = true; + E(0, 4) = E(4, 0) = true; + for (unsigned p = 0; p < 4; ++p) { + P(0) = p; + APState ap(A, B, C, E, P, 0); + auto dm_before = ap.to_density_matrix(); + ap.normal_form(); + auto dm_after = ap.to_density_matrix(); + // Just check using statevector; too much changes in each case to nicely + // test the matrices + CHECK(dm_after.isApprox(dm_before, EPS)); + } + } +} + SCENARIO("CZ cases") { GIVEN("CZ on free qubits") { MatrixXb A = MatrixXb::Zero(4, 4); @@ -278,6 +396,165 @@ SCENARIO("V cases") { test_apply_gate_dm(A, B, C, E, P, OpType::V, {0}); } } + GIVEN("V on a mixed state with zero A") { + for (unsigned p = 0; p < 4; ++p) { + MatrixXb A = MatrixXb::Zero(7, 7); + VectorXb B = VectorXb::Zero(7); + MatrixXb C = MatrixXb::Zero(7, 7); + MatrixXb E = MatrixXb::Zero(7, 7); + Eigen::VectorXi P = Eigen::VectorXi::Zero(7); + C(0, 0) = C(0, 2) = true; + C(1, 0) = C(1, 1) = C(1, 2) = C(1, 3) = C(1, 4) = true; + C(2, 0) = C(2, 2) = C(0, 3) = true; + C(3, 0) = C(3, 1) = C(3, 2) = C(3, 5) = C(3, 6) = true; + E(0, 3) = E(3, 0) = true; + E(0, 4) = E(4, 0) = true; + E(0, 5) = E(5, 0) = true; + E(0, 6) = E(6, 0) = true; + P(0) = p; + test_apply_gate_dm(A, B, C, E, P, OpType::V, {0}); + } + } +} + +SCENARIO("Qubit Reset") { + GIVEN("Reset a qubit with a local state") { + // Qubit 0 is in the |0> state + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb C = MatrixXb::Zero(3, 3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + A(0, 0) = true; + A(1, 1) = A(1, 2) = true; + B(1) = true; + P(2) = 1; + APState correct(A, B, C, E, P, 0); + // correct is already in normal form + for (unsigned s = 0; s < 6; ++s) { + // s = 0,1,2,3: XY basis states + // s = 4: |0> + // s = 5: |1> + A(0, 0) = (s < 4) ? false : true; + P(0) = (s < 4) ? s : 0; + B(0) = (s == 5); + APState ap(A, B, C, E, P, 0); + ap.apply_gate(OpType::Reset, {0}); + // Check up to global phase + ap.normal_form(); + ap.phase = correct.phase; + CHECK(ap == correct); + } + } + GIVEN("Reset one side of a Bell state") { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb C = MatrixXb::Zero(3, 3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + A(0, 0) = A(0, 1) = true; + APState ap(A, B, C, E, P, 0); + ap.apply_gate(OpType::Reset, {0}); + ap.normal_form(); + // Qubit 0 ends in |0> + A(0, 1) = false; + // Qubit 1 ends in maximally-mixed state + C(0, 1) = true; + APState correct(A, B, C, E, P, 0); + CHECK(ap == correct); + } + GIVEN("Reset on a normal form state") { + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 3) = true; + B(0) = true; + C(0, 1) = C(0, 2) = true; + E(1, 3) = E(3, 1) = true; + E(2, 3) = E(3, 2) = true; + P(2) = 1; + P(3) = 2; + APState ap(A, B, C, E, P, 0); + WHEN("Apply to Qubit 0") { + ap.apply_gate(OpType::Reset, {0}); + ap.normal_form(); + // Qubit 0 in state |0> + A(0, 1) = A(0, 3) = false; + B(0) = false; + // A row becomes a C row (combine with other row for gaussian form) + C(1, 2) = C(1, 3) = true; + // More gaussian steps + C(0, 2) = false; + C(0, 3) = true; + // LC about C(1, -) to remove P(2) + E(2, 3) = E(3, 2) = false; + P(2) = 0; + P(3) = 1; + APState correct(A, B, C, E, P, {0}); + // Check equality up to global phase + ap.phase = correct.phase; + CHECK(ap == correct); + } + WHEN("Apply to Qubit 1") { + ap.apply_gate(OpType::Reset, {1}); + ap.normal_form(); + // Correct form verified by hand + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 1) = true; + C(0, 0) = C(0, 3) = true; + C(1, 2) = true; + E(0, 3) = E(3, 0) = true; + E(2, 3) = E(3, 2) = true; + P(3) = 2; + APState correct(A, B, C, E, P, {0}); + // Check equality up to global phase + ap.phase = correct.phase; + CHECK(ap == correct); + } + WHEN("Apply to Qubit 2") { + ap.apply_gate(OpType::Reset, {2}); + ap.normal_form(); + // Correct form verified by hand + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 0) = A(0, 1) = A(0, 3) = true; + A(1, 2) = true; + B(0) = true; + C(0, 1) = true; + C(1, 3) = true; + APState correct(A, B, C, E, P, {0}); + // Check equality up to global phase + ap.phase = correct.phase; + CHECK(ap == correct); + } + WHEN("Apply to Qubit 3") { + ap.apply_gate(OpType::Reset, {3}); + ap.normal_form(); + // Correct form verified by hand + MatrixXb A = MatrixXb::Zero(4, 4); + VectorXb B = VectorXb::Zero(4); + MatrixXb C = MatrixXb::Zero(4, 4); + MatrixXb E = MatrixXb::Zero(4, 4); + Eigen::VectorXi P = Eigen::VectorXi::Zero(4); + A(0, 3) = true; + C(0, 0) = C(0, 2) = true; + C(1, 1) = C(1, 2) = true; + P(2) = 1; + APState correct(A, B, C, E, P, {0}); + // Check equality up to global phase + ap.phase = correct.phase; + CHECK(ap == correct); + } + } } SCENARIO("Gate encodings") { @@ -311,6 +588,30 @@ SCENARIO("Gate encodings") { test_apply_gate(A, B, E, P, com.first, com.second); } } + GIVEN("Check Z actions on inputs of ChoiAPState") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Identity(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb C = MatrixXb::Zero(3, 3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + qubit_vector_t qbs; + for (unsigned q : com.second) qbs.push_back(Qubit(q)); + test_apply_gate_dm_input(A, B, C, E, P, com.first, qbs); + } + } + GIVEN("Check X actions on inputs of ChoiAPState") { + for (const auto& com : test_gates) { + MatrixXb A = MatrixXb::Zero(3, 3); + VectorXb B = VectorXb::Zero(3); + MatrixXb C = MatrixXb::Zero(3, 3); + MatrixXb E = MatrixXb::Zero(3, 3); + Eigen::VectorXi P = Eigen::VectorXi::Zero(3); + qubit_vector_t qbs; + for (unsigned q : com.second) qbs.push_back(Qubit(q)); + test_apply_gate_dm_input(A, B, C, E, P, com.first, qbs); + } + } } SCENARIO("Loading from a statevector") { @@ -442,6 +743,35 @@ SCENARIO("Converting from/to a circuit") { // Ignore phase by setting them to match APState correct(A, B, C, E, P, ap.ap_.phase); CHECK(ap.ap_ == correct); + std::pair res_uni = + choi_apstate_to_unitary_extension_circuit(ap, {Qubit(1)}, {Qubit(0)}); + // Rebuild state by initialising, post-selecting, etc. + ChoiAPState res_ap = circuit_to_choi_apstate(res_uni.first); + qubit_map_t perm; + for (const std::pair& p : res_uni.second) + perm.insert({p.second, p.first}); + res_ap.rename_qubits(perm, ChoiAPState::TableauSegment::Output); + // Post-select/initialise + res_ap.post_select(Qubit(1), ChoiAPState::TableauSegment::Input); + res_ap.post_select(Qubit(0), ChoiAPState::TableauSegment::Output); + // Collapsing q[4] in X basis as per circ + res_ap.apply_gate( + OpType::H, {Qubit(4)}, ChoiAPState::TableauSegment::Output); + res_ap.collapse_qubit(Qubit(4), ChoiAPState::TableauSegment::Output); + res_ap.apply_gate( + OpType::H, {Qubit(4)}, ChoiAPState::TableauSegment::Output); + // Discarding q[0] also removes Z row for q[0], so recreate this by + // XCollapse at input + res_ap.apply_gate( + OpType::H, {Qubit(0)}, ChoiAPState::TableauSegment::Input); + res_ap.collapse_qubit(Qubit(0), ChoiAPState::TableauSegment::Input); + res_ap.apply_gate( + OpType::H, {Qubit(0)}, ChoiAPState::TableauSegment::Input); + res_ap.canonical_column_order(); + res_ap.normal_form(); + // Mixed state, so only guaranteed up to phase + res_ap.ap_.phase = ap.ap_.phase; + CHECK(res_ap == ap); THEN("Serialize and deserialize") { nlohmann::json j_ap = ap; ChoiAPState ap2({}); From d00087adaa8743c851a0d05416b281bf3016b36a Mon Sep 17 00:00:00 2001 From: Will Simmons Date: Thu, 5 Dec 2024 09:13:13 +0000 Subject: [PATCH 19/19] Bump tket version --- pytket/conanfile.py | 2 +- tket/conanfile.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pytket/conanfile.py b/pytket/conanfile.py index c50c6c9909..bc553d5788 100644 --- a/pytket/conanfile.py +++ b/pytket/conanfile.py @@ -38,7 +38,7 @@ def requirements(self): self.requires("pybind11_json/0.2.14") self.requires("symengine/0.13.0") self.requires("tkassert/0.3.4@tket/stable") - self.requires("tket/1.3.53@tket/stable") + self.requires("tket/1.3.54@tket/stable") self.requires("tklog/0.3.3@tket/stable") self.requires("tkrng/0.3.3@tket/stable") self.requires("tktokenswap/0.3.9@tket/stable") diff --git a/tket/conanfile.py b/tket/conanfile.py index 2602d319f1..dcc80f7bb7 100644 --- a/tket/conanfile.py +++ b/tket/conanfile.py @@ -23,7 +23,7 @@ class TketConan(ConanFile): name = "tket" - version = "1.3.53" + version = "1.3.54" package_type = "library" license = "Apache 2" homepage = "https://github.com/CQCL/tket"