GCC Code Coverage Report

Directory: ./
File: Utils/MatrixAnalysis.cpp
Date: 2022-10-15 05:10:18
Warnings: 2 unchecked decisions!
Exec Total Coverage
Lines: 296 312 94.9%
Functions: 27 30 90.0%
Branches: 445 794 56.0%
Decisions: 86 94 91.5%
Line Branch Decision Exec Source
1 // Copyright 2019-2022 Cambridge Quantum Computing
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 #include "MatrixAnalysis.hpp"
16
17 #include <fstream>
18 #include <iostream>
19 #include <limits>
20 #include <map>
21 #include <optional>
22 #include <sstream>
23 #include <tkassert/Assert.hpp>
24 #include <utility>
25 #include <vector>
26
27 #include "Utils/EigenConfig.hpp"
28
29 namespace tket {
30
31 12159 bool is_unitary(const Eigen::MatrixXcd &U, double tol) {
32 12159 unsigned m = U.rows(), n = U.cols();
33
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12159 times.
1/2
✓ Decision 'true' taken 12159 times.
✗ Decision 'false' not taken.
 12159 if (m != n) return false;
34
3/6
✓ Branch 2 taken 12159 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12159 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12159 times.
✗ Branch 9 not taken.
 12159 return Eigen::MatrixXcd::Identity(n, n).isApprox(U.adjoint() * U, tol);
35 }
36
37 2 bool is_projector(const Eigen::MatrixXcd &P, double tol) {
38 2 unsigned m = P.rows(), n = P.cols();
39
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
1/2
✓ Decision 'true' taken 2 times.
✗ Decision 'false' not taken.
 2 if (m != n) return false;
40
6/12
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
 2 return P.isApprox(P * P, tol) && P.isApprox(P.adjoint());
41 }
42
43 86 std::pair<MatrixXb, MatrixXb> binary_LLT_decomposition(const MatrixXb &a) {
44 /*
45 * Aaronson-Gottesman: Improved Simulation of Stabilizer Circuits, Lemma 7
46 * For any symmetric, binary matrix A, there exists a diagonal matrix D
47 * and an invertible lower-triangular matrix L such that A + D = LL^T.
48 * Given A, returns the pair <L, D>
49 */
50 86 unsigned n = a.rows();
51
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 MatrixXb lo = MatrixXb::Identity(n, n);
52
2/2
✓ Branch 0 taken 278 times.
✓ Branch 1 taken 86 times.
2/2
✓ Decision 'true' taken 278 times.
✓ Decision 'false' taken 86 times.
 364 for (unsigned j = 0; j < n; j++) {
53
2/2
✓ Branch 0 taken 338 times.
✓ Branch 1 taken 278 times.
2/2
✓ Decision 'true' taken 338 times.
✓ Decision 'false' taken 278 times.
 616 for (unsigned i = j + 1; i < n; i++) {
54
1/2
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
 338 bool sum = a(i, j);
55
2/2
✓ Branch 0 taken 190 times.
✓ Branch 1 taken 338 times.
2/2
✓ Decision 'true' taken 190 times.
✓ Decision 'false' taken 338 times.
 528 for (unsigned k = 0; k < j; k++) {
56
5/8
✓ Branch 1 taken 190 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 180 times.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
 190 sum ^= lo(i, k) && lo(j, k);
57 }
58
1/2
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
 338 lo(i, j) = sum;
59 }
60 }
61
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 MatrixXb d = MatrixXb::Zero(n, n);
62
2/2
✓ Branch 0 taken 278 times.
✓ Branch 1 taken 86 times.
2/2
✓ Decision 'true' taken 278 times.
✓ Decision 'false' taken 86 times.
 364 for (unsigned i = 0; i < n; i++) {
63
1/2
✓ Branch 1 taken 278 times.
✗ Branch 2 not taken.
 278 bool sum = a(i, i);
64
2/2
✓ Branch 0 taken 954 times.
✓ Branch 1 taken 278 times.
2/2
✓ Decision 'true' taken 954 times.
✓ Decision 'false' taken 278 times.
 1232 for (unsigned k = 0; k < n; k++) {
65
1/2
✓ Branch 1 taken 954 times.
✗ Branch 2 not taken.
 954 sum ^= lo(i, k);
66 } // diagonal element of LL^T is just parity of row of L
67
1/2
✓ Branch 1 taken 278 times.
✗ Branch 2 not taken.
 278 d(i, i) = sum;
68 }
69
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 172 return {lo, d};
70 86 }
71
72 215 std::vector<std::pair<unsigned, unsigned>> gaussian_elimination_col_ops(
73 const MatrixXb &a, unsigned blocksize) {
74
2/4
✓ Branch 2 taken 215 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 215 times.
✗ Branch 6 not taken.
 215 return gaussian_elimination_row_ops(a.transpose(), blocksize);
75 }
76
77 struct MatrixXbBlockCompare {
78 22038 bool operator()(
79 const MatrixXb::BlockXpr &lhs, const MatrixXb::BlockXpr &rhs) const {
80 TKET_ASSERT(lhs.rows() == rhs.rows());
81 TKET_ASSERT(lhs.cols() == rhs.cols());
82
2/2
✓ Branch 1 taken 22038 times.
✓ Branch 2 taken 406 times.
2/2
✓ Decision 'true' taken 22038 times.
✓ Decision 'false' taken 406 times.
 22444 for (unsigned r = 0; r < lhs.rows(); ++r) {
83
2/2
✓ Branch 1 taken 49988 times.
✓ Branch 2 taken 406 times.
2/2
✓ Decision 'true' taken 49988 times.
✓ Decision 'false' taken 406 times.
 50394 for (unsigned c = 0; c < lhs.cols(); ++c) {
84
2/2
✓ Branch 2 taken 12596 times.
✓ Branch 3 taken 37392 times.
2/2
✓ Decision 'true' taken 37392 times.
✓ Decision 'false' taken 12596 times.
 49988 if (lhs(r, c) < rhs(r, c)) return true;
85
2/2
✓ Branch 2 taken 9036 times.
✓ Branch 3 taken 28356 times.
2/2
✓ Decision 'true' taken 406 times.
✓ Decision 'false' taken 36986 times.
 37392 if (lhs(r, c) > rhs(r, c)) return false;
86 }
87 }
88 406 return false;
89 }
90 };
91
92 /* see https://web.eecs.umich.edu/~imarkov/pubs/jour/qic08-cnot.pdf for a full
93 * explanation of this technique */
94 /* K. Patel, I. Markov, J. Hayes. Optimal Synthesis of Linear Reversible
95 Circuits. QIC 2008 */
96 532 std::vector<std::pair<unsigned, unsigned>> gaussian_elimination_row_ops(
97 const MatrixXb &a, unsigned blocksize) {
98 532 std::vector<std::pair<unsigned, unsigned>> ops;
99
1/2
✓ Branch 1 taken 532 times.
✗ Branch 2 not taken.
 532 MatrixXb m = a;
100 532 unsigned rows = m.rows();
101 532 unsigned cols = m.cols();
102 532 std::vector<unsigned> pcols; // we know which columns have non-zero entries
103 // (to save actually transposing the matrix)
104 532 unsigned pivot_row = 0;
105 532 unsigned ceiling =
106 532 (cols + blocksize - 1) / blocksize; // ceil(cols/blocksize)
107
108 // Get to upper echelon form
109
2/2
✓ Branch 0 taken 645 times.
✓ Branch 1 taken 532 times.
2/2
✓ Decision 'true' taken 645 times.
✓ Decision 'false' taken 532 times.
 1177 for (unsigned sec = 0; sec < ceiling; ++sec) {
110 // determine column range for this section
111 // if it hits cols don't go any higher
112 645 unsigned i0 = sec * blocksize;
113 645 unsigned i1 = std::min(cols, (sec + 1) * blocksize);
114
115 /* first, try to eliminate sub-rows (ie chunks), for greater speed than
116 * naively doing gaussian elim. */
117 645 std::map<MatrixXb::BlockXpr, unsigned, MatrixXbBlockCompare> chunks;
118
2/2
✓ Branch 0 taken 3053 times.
✓ Branch 1 taken 645 times.
2/2
✓ Decision 'true' taken 3053 times.
✓ Decision 'false' taken 645 times.
 3698 for (unsigned r = pivot_row; r < rows; ++r) {
119
1/2
✓ Branch 1 taken 3053 times.
✗ Branch 2 not taken.
 3053 MatrixXb::BlockXpr t = m.block(r, i0, 1, i1 - i0);
120
4/6
✓ Branch 1 taken 3053 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3053 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 413 times.
✓ Branch 7 taken 2640 times.
2/2
✓ Decision 'true' taken 2640 times.
✓ Decision 'false' taken 413 times.
 3053 if (t == MatrixXb::Zero(1, i1 - i0)) continue;
121 /* if first copy of pattern, save. If duplicate then remove by adding
122 * rows*/
123
1/2
✓ Branch 1 taken 2640 times.
✗ Branch 2 not taken.
 2640 auto chunk_it = chunks.find(t);
124
2/2
✓ Branch 2 taken 152 times.
✓ Branch 3 taken 2488 times.
2/2
✓ Decision 'true' taken 1937 times.
✓ Decision 'false' taken 703 times.
 2640 if (chunk_it != chunks.end()) {
125
2/2
✓ Branch 0 taken 1785 times.
✓ Branch 1 taken 152 times.
2/2
✓ Decision 'true' taken 1785 times.
✓ Decision 'false' taken 152 times.
 1937 for (unsigned k = 0; k < cols; ++k) {
126
2/4
✓ Branch 2 taken 1785 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1785 times.
✗ Branch 6 not taken.
 1785 m(r, k) ^= m(chunk_it->second, k);
127 }
128
1/2
✓ Branch 3 taken 152 times.
✗ Branch 4 not taken.
 152 ops.push_back({chunk_it->second, r});
129 } else {
130
1/2
✓ Branch 2 taken 2488 times.
✗ Branch 3 not taken.
 2488 chunks.insert({t, r});
131 }
132 }
133 /* do gaussian elim. on remaining entries */
134
2/2
✓ Branch 0 taken 2591 times.
✓ Branch 1 taken 645 times.
2/2
✓ Decision 'true' taken 2591 times.
✓ Decision 'false' taken 645 times.
 3236 for (unsigned col = i0; col < i1; ++col) {
135 // find first 1 element in column after pivot_row
136 2591 unsigned first_1 = pivot_row;
137
7/8
✓ Branch 0 taken 3022 times.
✓ Branch 1 taken 113 times.
✓ Branch 3 taken 3022 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 544 times.
✓ Branch 6 taken 2478 times.
✓ Branch 7 taken 544 times.
✓ Branch 8 taken 2591 times.
0/1
? Decision couldn't be analyzed.
 3135 while (first_1 < rows && !m(first_1, col)) {
138 544 ++first_1;
139 }
140
2/2
✓ Branch 0 taken 113 times.
✓ Branch 1 taken 2478 times.
2/2
✓ Decision 'true' taken 2478 times.
✓ Decision 'false' taken 113 times.
 2591 if (first_1 == rows) continue;
141
142 // pull back to pivot
143
2/2
✓ Branch 0 taken 116 times.
✓ Branch 1 taken 2362 times.
2/2
✓ Decision 'true' taken 767 times.
✓ Decision 'false' taken 1711 times.
 2478 if (first_1 != pivot_row) {
144
2/2
✓ Branch 0 taken 651 times.
✓ Branch 1 taken 116 times.
2/2
✓ Decision 'true' taken 651 times.
✓ Decision 'false' taken 116 times.
 767 for (unsigned k = 0; k < cols; ++k) {
145
2/4
✓ Branch 1 taken 651 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 651 times.
✗ Branch 5 not taken.
 651 m(pivot_row, k) ^= m(first_1, k);
146 }
147
1/2
✓ Branch 2 taken 116 times.
✗ Branch 3 not taken.
 116 ops.push_back({first_1, pivot_row});
148 }
149
150 // clear all entries below pivot
151
2/2
✓ Branch 1 taken 7116 times.
✓ Branch 2 taken 2478 times.
2/2
✓ Decision 'true' taken 7116 times.
✓ Decision 'false' taken 2478 times.
 9594 for (unsigned r = std::max(pivot_row + 1, first_1); r < rows; ++r) {
152
3/4
✓ Branch 1 taken 7116 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 772 times.
✓ Branch 4 taken 6344 times.
2/2
✓ Decision 'true' taken 5920 times.
✓ Decision 'false' taken 1196 times.
 7116 if (m(r, col)) {
153
2/2
✓ Branch 0 taken 5148 times.
✓ Branch 1 taken 772 times.
2/2
✓ Decision 'true' taken 5148 times.
✓ Decision 'false' taken 772 times.
 5920 for (unsigned k = 0; k < cols; ++k) {
154
2/4
✓ Branch 1 taken 5148 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5148 times.
✗ Branch 5 not taken.
 5148 m(r, k) ^= m(pivot_row, k);
155 }
156
1/2
✓ Branch 2 taken 772 times.
✗ Branch 3 not taken.
 772 ops.push_back({pivot_row, r});
157 }
158 }
159
160 // record that we pivoted for this column
161
1/2
✓ Branch 1 taken 2478 times.
✗ Branch 2 not taken.
 2478 pcols.push_back(col);
162 2478 ++pivot_row;
163 }
164 645 }
165
166 /* matrix is now in upper triangular form; matrix is reduced to diagonal */
167
168 532 --pivot_row;
169
170
2/2
✓ Branch 0 taken 645 times.
✓ Branch 1 taken 532 times.
2/2
✓ Decision 'true' taken 645 times.
✓ Decision 'false' taken 532 times.
 1177 for (unsigned sec = ceiling; sec-- > 0;) {
171 645 unsigned i0 = sec * blocksize;
172 645 unsigned i1 = std::min(cols, (sec + 1) * blocksize);
173
174 645 std::map<MatrixXb::BlockXpr, unsigned, MatrixXbBlockCompare> chunks;
175
2/2
✓ Branch 0 taken 3120 times.
✓ Branch 1 taken 645 times.
2/2
✓ Decision 'true' taken 3120 times.
✓ Decision 'false' taken 645 times.
 3765 for (unsigned r = pivot_row + 1; r-- > 0;) {
176
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 MatrixXb::BlockXpr t = m.block(r, i0, 1, i1 - i0);
177
4/6
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3120 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 557 times.
✓ Branch 7 taken 2563 times.
2/2
✓ Decision 'true' taken 2563 times.
✓ Decision 'false' taken 557 times.
 3120 if (t == MatrixXb::Zero(1, i1 - i0)) continue;
178
179
1/2
✓ Branch 1 taken 2563 times.
✗ Branch 2 not taken.
 2563 auto chunk_it = chunks.find(t);
180
2/2
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 2512 times.
2/2
✓ Decision 'true' taken 672 times.
✓ Decision 'false' taken 1891 times.
 2563 if (chunk_it != chunks.end()) {
181
2/2
✓ Branch 0 taken 621 times.
✓ Branch 1 taken 51 times.
2/2
✓ Decision 'true' taken 621 times.
✓ Decision 'false' taken 51 times.
 672 for (unsigned k = 0; k < cols; ++k) {
182
2/4
✓ Branch 2 taken 621 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 621 times.
✗ Branch 6 not taken.
 621 m(r, k) ^= m(chunk_it->second, k);
183 }
184
1/2
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
 51 ops.push_back({chunk_it->second, r});
185 } else {
186
1/2
✓ Branch 2 taken 2512 times.
✗ Branch 3 not taken.
 2512 chunks.insert({t, r});
187 }
188 }
189
7/8
✓ Branch 1 taken 2582 times.
✓ Branch 2 taken 541 times.
✓ Branch 4 taken 2478 times.
✓ Branch 5 taken 104 times.
✓ Branch 7 taken 2478 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2478 times.
✓ Branch 10 taken 645 times.
0/1
? Decision couldn't be analyzed.
 3123 while (!pcols.empty() && i0 <= pcols.back() && pcols.back() < i1) {
190 2478 unsigned pcol = pcols.back();
191 2478 pcols.pop_back();
192
2/2
✓ Branch 0 taken 6877 times.
✓ Branch 1 taken 2478 times.
2/2
✓ Decision 'true' taken 6877 times.
✓ Decision 'false' taken 2478 times.
 9355 for (unsigned r = 0; r < pivot_row; ++r) {
193
3/4
✓ Branch 1 taken 6877 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 632 times.
✓ Branch 4 taken 6245 times.
2/2
✓ Decision 'true' taken 4960 times.
✓ Decision 'false' taken 1917 times.
 6877 if (m(r, pcol)) {
194
2/2
✓ Branch 0 taken 4328 times.
✓ Branch 1 taken 632 times.
2/2
✓ Decision 'true' taken 4328 times.
✓ Decision 'false' taken 632 times.
 4960 for (unsigned k = 0; k < cols; ++k) {
195
2/4
✓ Branch 1 taken 4328 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4328 times.
✗ Branch 5 not taken.
 4328 m(r, k) ^= m(pivot_row, k);
196 }
197
1/2
✓ Branch 2 taken 632 times.
✗ Branch 3 not taken.
 632 ops.push_back({pivot_row, r});
198 }
199 }
200 2478 --pivot_row;
201 }
202 645 }
203
204 1064 return ops;
205 532 }
206
207 3 static Eigen::PermutationMatrix<Eigen::Dynamic> qubit_permutation(
208 unsigned n_qubits) {
209 3 Eigen::PermutationMatrix<Eigen::Dynamic> perm(1u << n_qubits);
210
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 3 times.
2/2
✓ Decision 'true' taken 20 times.
✓ Decision 'false' taken 3 times.
 23 for (unsigned i = 0; i < (1u << n_qubits); i++) {
211 20 unsigned rev = 0;
212 20 unsigned forwards = i;
213
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 20 times.
2/2
✓ Decision 'true' taken 56 times.
✓ Decision 'false' taken 20 times.
 76 for (unsigned q = 0; q < n_qubits; q++) {
214 56 rev <<= 1;
215 56 rev += forwards % 2;
216 56 forwards >>= 1;
217 }
218
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 perm.indices()[i] = rev;
219 }
220 3 return perm;
221 }
222
223 5868 Eigen::PermutationMatrix<Eigen::Dynamic> lift_perm(
224 const std::map<unsigned, unsigned> &p) {
225 5868 unsigned n = p.size();
226 5868 Eigen::PermutationMatrix<Eigen::Dynamic> pm(1u << n);
227
2/2
✓ Branch 0 taken 129346 times.
✓ Branch 1 taken 5868 times.
2/2
✓ Decision 'true' taken 129346 times.
✓ Decision 'false' taken 5868 times.
 135214 for (unsigned i = 0; i < (1u << n); ++i) {
228 129346 unsigned target = 0;
229 129346 unsigned mask = 1u << n;
230
2/2
✓ Branch 0 taken 907740 times.
✓ Branch 1 taken 129346 times.
2/2
✓ Decision 'true' taken 907740 times.
✓ Decision 'false' taken 129346 times.
 1037086 for (unsigned q = 0; q < n; ++q) {
231 907740 mask /= 2;
232
2/2
✓ Branch 0 taken 453870 times.
✓ Branch 1 taken 453870 times.
2/2
✓ Decision 'true' taken 453870 times.
✓ Decision 'false' taken 453870 times.
 907740 if (i & mask) {
233
1/2
✓ Branch 1 taken 453870 times.
✗ Branch 2 not taken.
 453870 target |= 1u << (n - 1 - p.at(q));
234 }
235 }
236
1/2
✓ Branch 2 taken 129346 times.
✗ Branch 3 not taken.
 129346 pm.indices()[i] = target;
237 }
238 5868 return pm;
239 }
240
241 5868 static Eigen::PermutationMatrix<Eigen::Dynamic> qubit_permutation(
242 const qubit_map_t &qmap) {
243 5868 std::map<Qubit, unsigned> q_indices;
244 5868 std::map<unsigned, unsigned> uq_map;
245 5868 unsigned qi = 0;
246
2/2
✓ Branch 4 taken 16188 times.
✓ Branch 5 taken 5868 times.
2/2
✓ Decision 'true' taken 16188 times.
✓ Decision 'false' taken 5868 times.
 22056 for (const std::pair<const Qubit, Qubit> &pair : qmap) {
247
1/2
✓ Branch 2 taken 16188 times.
✗ Branch 3 not taken.
 16188 q_indices.insert({pair.first, qi});
248 16188 ++qi;
249 }
250
2/2
✓ Branch 4 taken 16188 times.
✓ Branch 5 taken 5868 times.
2/2
✓ Decision 'true' taken 16188 times.
✓ Decision 'false' taken 5868 times.
 22056 for (const std::pair<const Qubit, Qubit> &pair : qmap) {
251
1/2
✓ Branch 1 taken 16188 times.
✗ Branch 2 not taken.
 16188 unsigned in = q_indices[pair.first];
252
1/2
✓ Branch 1 taken 16188 times.
✗ Branch 2 not taken.
 16188 unsigned out = q_indices[pair.second];
253
1/2
✓ Branch 2 taken 16188 times.
✗ Branch 3 not taken.
 16188 uq_map.insert({in, out});
254 }
255
1/2
✓ Branch 1 taken 5868 times.
✗ Branch 2 not taken.
 11736 return lift_perm(uq_map);
256 5868 }
257
258 static const Eigen::PermutationMatrix<4, 4> SWAP = qubit_permutation(2);
259
260 1 Eigen::Matrix4cd reverse_indexing(const Eigen::Matrix4cd &m) {
261
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 return SWAP * m * SWAP;
262 }
263
264 ✗ Matrix8cd reverse_indexing(const Matrix8cd &m) {
265 return static_cast<const Matrix8cd &>(
266 ✗ reverse_indexing(static_cast<const Eigen::MatrixXcd &>(m)));
267 }
268
269 1 Eigen::MatrixXcd reverse_indexing(const Eigen::MatrixXcd &m) {
270 1 unsigned dim = m.rows();
271
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 unsigned n_qubits = get_number_of_qubits(dim);
272
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::PermutationMatrix<Eigen::Dynamic> perm = qubit_permutation(n_qubits);
273
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 2 return perm * m * perm;
274 1 }
275
276 1 Eigen::VectorXcd reverse_indexing(const Eigen::VectorXcd &v) {
277 1 unsigned dim = v.size();
278
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 unsigned n_qubits = get_number_of_qubits(dim);
279
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::PermutationMatrix<Eigen::Dynamic> perm = qubit_permutation(n_qubits);
280
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 return perm * v;
281 1 }
282
283 5868 Eigen::MatrixXcd apply_qubit_permutation(
284 const Eigen::MatrixXcd &m, const qubit_map_t &perm) {
285
1/2
✓ Branch 1 taken 5868 times.
✗ Branch 2 not taken.
 5868 Eigen::PermutationMatrix<Eigen::Dynamic> perm_m = qubit_permutation(perm);
286
2/4
✓ Branch 1 taken 5868 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5868 times.
✗ Branch 5 not taken.
 11736 return perm_m * m;
287 5868 }
288
289 ✗ Eigen::VectorXcd apply_qubit_permutation(
290 const Eigen::VectorXcd &v, const qubit_map_t &perm) {
291 ✗ Eigen::PermutationMatrix<Eigen::Dynamic> perm_m = qubit_permutation(perm);
292 ✗ return perm_m * v;
293 }
294
295 866 double trace_fidelity(double a, double b, double c) {
296 866 constexpr double g = PI / 2;
297 866 a *= g;
298 866 b *= g;
299 866 c *= g;
300 866 double trace_sq = 16. * (pow(cos(a) * cos(b) * cos(c), 2) +
301 866 pow(sin(a) * sin(b) * sin(c), 2));
302 866 return (4. + trace_sq) / 20.;
303 }
304
305 9951 inline double mod(double d, double max) { return d - max * floor(d / max); }
306
307 std::tuple<Eigen::Matrix4cd, std::array<double, 3>, Eigen::Matrix4cd>
308 3317 get_information_content(const Eigen::Matrix4cd &X) {
309 using ExpGate = std::array<double, 3>;
310 using Mat4 = Eigen::Matrix4cd;
311 using Vec4 = Eigen::Vector4cd;
312
313
3/6
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 3317 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3317 times.
 3317 if (!is_unitary(X)) {
314 ✗ throw std::invalid_argument(
315 ✗ "Non-unitary matrix passed to get_information_content");
316 }
317
318
3/6
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
 3317 Eigen::Matrix2cd PauliX, PauliY, PauliZ;
319
4/8
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 3317 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 3317 times.
✗ Branch 15 not taken.
 3317 PauliX << 0, 1, 1, 0;
320
4/8
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 3317 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3317 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
 3317 PauliY << 0, -i_, i_, 0;
321
4/8
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 3317 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 3317 times.
✗ Branch 15 not taken.
 3317 PauliZ << 1, 0, 0, -1;
322
323 // change of basis for SU(2) x SU(2) -> SO(4)
324
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 Mat4 MagicM;
325
16/32
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 3317 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 3317 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 3317 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 3317 times.
✗ Branch 25 not taken.
✓ Branch 28 taken 3317 times.
✗ Branch 29 not taken.
✓ Branch 32 taken 3317 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 3317 times.
✗ Branch 36 not taken.
✓ Branch 39 taken 3317 times.
✗ Branch 40 not taken.
✓ Branch 43 taken 3317 times.
✗ Branch 44 not taken.
✓ Branch 47 taken 3317 times.
✗ Branch 48 not taken.
✓ Branch 51 taken 3317 times.
✗ Branch 52 not taken.
✓ Branch 55 taken 3317 times.
✗ Branch 56 not taken.
✓ Branch 59 taken 3317 times.
✗ Branch 60 not taken.
 3317 MagicM << 1, 0, 0, i_, 0, i_, 1, 0, 0, i_, -1, 0, 1, 0, 0, -i_;
326
1/2
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
 3317 MagicM /= sqrt(2.); // unitary
327
328 // change to magic basis and make sure U in SU(4)
329
2/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
 3317 const auto norm_X = pow(X.determinant(), 0.25);
330
5/10
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
 3317 const Mat4 Xprime = MagicM.adjoint() * X * MagicM / norm_X;
331
332 // Find a common eigendecomposition of Re( X'.adjoint * X' ) and Im(
333 // X'.adjoint * X' ) Use pseudorandom linear comb to avoid issues with
334 // multiplicities > 1
335
3/6
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
 3317 Mat4 X2 = Xprime.transpose() * Xprime;
336 // For Clifford matrix, SelfAdjointEigenSolver seems to have a higher chance
337 // to produce eigenvectors that eventually lead to non-clifford angles when
338 // there are rounding errors. In this case, the matrix X2 should consist of
339 // 0, 1, -1, i and -i entries only, so clamp to these values.
340
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 X2 = X2.unaryExpr([](Complex x) -> Complex {
341
2/2
✓ Branch 1 taken 22142 times.
✓ Branch 2 taken 30930 times.
2/2
✓ Decision 'true' taken 22142 times.
✓ Decision 'false' taken 30930 times.
 53072 if (std::abs(x) < 1e-14) {
342 22142 return 0.;
343
2/2
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 30525 times.
2/2
✓ Decision 'true' taken 405 times.
✓ Decision 'false' taken 30525 times.
 30930 } else if (std::abs(x - 1.) < 1e-14) {
344 405 return 1.;
345
2/2
✓ Branch 2 taken 399 times.
✓ Branch 3 taken 30126 times.
2/2
✓ Decision 'true' taken 399 times.
✓ Decision 'false' taken 30126 times.
 30525 } else if (std::abs(x + 1.) < 1e-14) {
346 399 return -1.;
347
2/2
✓ Branch 2 taken 3568 times.
✓ Branch 3 taken 26558 times.
2/2
✓ Decision 'true' taken 3568 times.
✓ Decision 'false' taken 26558 times.
 30126 } else if (std::abs(x - i_) < 1e-14) {
348 3568 return i_;
349
2/2
✓ Branch 2 taken 652 times.
✓ Branch 3 taken 25906 times.
2/2
✓ Decision 'true' taken 652 times.
✓ Decision 'false' taken 25906 times.
 26558 } else if (std::abs(x + i_) < 1e-14) {
350 652 return -i_;
351 } else {
352 25906 return x;
353 }
354
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 });
355
2/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
 3317 const Eigen::Matrix4d X2real = X2.real();
356
2/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
 3317 const Eigen::Matrix4d X2imag = X2.imag();
357
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 Mat4 eigv;
358
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 Vec4 eigs;
359 3317 double r = 0.;
360 while (true) {
361 3317 r += 1 / PI;
362
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3317 times.
1/2
✓ Decision 'true' taken 3317 times.
✗ Decision 'false' not taken.
 3317 if (r >= 1) r -= 1;
363 Eigen::SelfAdjointEigenSolver<Eigen::Matrix4d> ces(
364
4/8
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
 3317 r * X2real + (1 - r) * X2imag);
365
2/4
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3317 times.
✗ Branch 6 not taken.
 3317 eigv = ces.eigenvectors().cast<Complex>();
366
5/10
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
 3317 eigs = (eigv.transpose() * X2 * eigv).diagonal();
367
368
7/14
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3317 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3317 times.
✗ Branch 20 not taken.
 3317 if (std::abs((X2 - eigv * eigs.asDiagonal() * eigv.adjoint()).sum()) <
369 EPS) {
370 // success!
371 3317 break;
372 }
373 }
374
375 Eigen::Vector4d thetas =
376
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 16585 eigs.unaryExpr([](const Complex lambda) { return std::arg(lambda) / 2.; })
377
2/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
 3317 .real();
378 // make sure exp(i thetas) is in SU(4) ie Σ thetas = 0
379
3/6
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
 3317 thetas(3) = -thetas.head<3>().sum();
380
381
3/6
✓ Branch 2 taken 3317 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3317 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3317 times.
✗ Branch 9 not taken.
 3317 const Mat4 tmp = (-i_ * thetas).asDiagonal();
382
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 const auto eigs_sqrt_inv = tmp.exp();
383
384
2/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
 3317 Mat4 Q2 = eigv.transpose();
385
4/8
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
 3317 Mat4 Q1 = Xprime * Q2.transpose() * eigs_sqrt_inv;
386
387 // we now have Xprime = Q1 * exp(i thetas) * Q2
388 // with Q1,Q2 in SO(4), exp(i thetas) in SU(4)
389
390 // transform to Pauli basis exp(i k_i σ_ii)
391
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 Eigen::Matrix4d basis_change;
392
15/30
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3317 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3317 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 3317 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 3317 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 3317 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 3317 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 3317 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 3317 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 3317 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 3317 times.
✗ Branch 44 not taken.
 3317 basis_change << 1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1,
393
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 1; // k3 = 0 always
394 Eigen::Vector3d k =
395
5/10
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3317 times.
✗ Branch 14 not taken.
 3317 (-.5 * basis_change * thetas / PI).head<3>().unaryExpr([](double d) {
396 9951 return mod(d, 4);
397
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 });
398
399 // we need to ensure that det(Q) == 1 so that K1,K2 \in SU(2) x SU(2)
400
3/4
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2010 times.
✓ Branch 5 taken 1307 times.
 3317 if (Q2.determinant().real() < 0) {
401
4/8
✓ Branch 1 taken 2010 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2010 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2010 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2010 times.
✗ Branch 11 not taken.
 2010 Q2.row(3) = -Q2.row(3);
402
4/8
✓ Branch 1 taken 2010 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2010 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2010 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2010 times.
✗ Branch 11 not taken.
 2010 Q1 = Xprime * Q2.transpose() * eigs_sqrt_inv;
403 }
404
405 // these are our local operations (left and right)
406
4/8
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
 3317 Mat4 K1 = MagicM * Q1 * MagicM.adjoint();
407
4/8
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3317 times.
✗ Branch 11 not taken.
 3317 Mat4 K2 = MagicM * Q2 * MagicM.adjoint();
408
409 // fix phase
410
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 3317 K1 *= norm_X;
411
412 // Finally, we got our ks
413
3/6
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3317 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3317 times.
✗ Branch 8 not taken.
 3317 ExpGate A{k(0), k(1), k(2)};
414
415 // K1,K2 in SU(2)xSU(2), A = Exp(i aσ_XX + i bσ_YY + i cσ_ZZ)
416
1/2
✓ Branch 1 taken 3317 times.
✗ Branch 2 not taken.
 6634 return std::tuple<Mat4, ExpGate, Mat4>{K1, A, K2};
417 }
418
419 6404 std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd> kronecker_decomposition(
420 Eigen::Matrix4cd &U) {
421 using Mat4 = Eigen::Matrix4cd;
422
4/8
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
 6404 Mat4 Up = U / pow(U.determinant(), 0.25);
423
1/2
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
 6404 Mat4 U_mod;
424
14/28
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6404 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6404 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6404 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6404 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 6404 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 6404 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 6404 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 6404 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 6404 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 6404 times.
✗ Branch 41 not taken.
 12808 U_mod << Up(0, 0), Up(1, 0), Up(0, 1), Up(1, 1), Up(2, 0), Up(3, 0), Up(2, 1),
425
14/28
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6404 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6404 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6404 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6404 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 6404 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 6404 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 6404 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 6404 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 6404 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 6404 times.
✗ Branch 41 not taken.
 6404 Up(3, 1), Up(0, 2), Up(1, 2), Up(0, 3), Up(1, 3), Up(2, 2), Up(3, 2),
426
4/8
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
 6404 Up(2, 3), Up(3, 3);
427 Eigen::JacobiSVD<Mat4, Eigen::NoQRPreconditioner> svd(
428 U_mod, Eigen::DecompositionOptions::ComputeFullU |
429
1/2
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
 6404 Eigen::DecompositionOptions::ComputeFullV);
430
2/4
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
 6404 Eigen::Vector4cd sings = svd.singularValues();
431
2/4
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
 6404 Mat4 l = svd.matrixU();
432
3/6
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
 6404 Mat4 r = svd.matrixV().conjugate();
433
2/4
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
 6404 Eigen::Matrix2cd u0, u1;
434
8/16
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6404 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6404 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6404 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6404 times.
✗ Branch 23 not taken.
 6404 u0 << l(0, 0), l(2, 0), l(1, 0), l(3, 0);
435
8/16
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6404 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6404 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6404 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6404 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6404 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6404 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6404 times.
✗ Branch 23 not taken.
 6404 u1 << r(0, 0), r(2, 0), r(1, 0), r(3, 0);
436
2/4
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 6404 times.
✗ Branch 6 not taken.
 6404 u0 *= sqrt(sings[0]);
437
2/4
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 6404 times.
✗ Branch 6 not taken.
 6404 u1 *= sqrt(sings[0]);
438
1/2
✓ Branch 1 taken 6404 times.
✗ Branch 2 not taken.
 12808 return {u0, u1};
439 6404 }
440
441 600451 unsigned get_matrix_size(unsigned number_of_qubits) {
442 // Left-shifting by the number of bits stored by the int type,
443 // or more, is UNDEFINED by the C++ standard!
444 // It does NOT necessarily give zero!
445
2/2
✓ Branch 0 taken 599450 times.
✓ Branch 1 taken 1001 times.
 600451 if (number_of_qubits < std::numeric_limits<unsigned>::digits) {
446 599450 unsigned size = 1;
447 599450 size <<= number_of_qubits;
448 599450 return size;
449 }
450
1/2
✓ Branch 1 taken 1001 times.
✗ Branch 2 not taken.
 1001 std::stringstream ss;
451
3/6
✓ Branch 1 taken 1001 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1001 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1001 times.
✗ Branch 8 not taken.
 1001 ss << "get_matrix_size for " << number_of_qubits << " qubits; overflow!";
452
2/4
✓ Branch 2 taken 1001 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1001 times.
✗ Branch 6 not taken.
 1001 throw std::runtime_error(ss.str());
453 1001 }
454
455 204569 unsigned get_number_of_qubits(unsigned matrix_size) {
456 204569 unsigned n_qubits = (unsigned)log2(matrix_size);
457
3/4
✓ Branch 1 taken 204569 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 201347 times.
✓ Branch 4 taken 3222 times.
 204569 if (get_matrix_size(n_qubits) == matrix_size) {
458 201347 return n_qubits;
459 }
460
1/2
✓ Branch 1 taken 3222 times.
✗ Branch 2 not taken.
 3222 std::stringstream ss;
461
2/4
✓ Branch 1 taken 3222 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3222 times.
✗ Branch 5 not taken.
 3222 ss << "get_number_of_qubits: matrix size " << matrix_size
462
1/2
✓ Branch 1 taken 3222 times.
✗ Branch 2 not taken.
 3222 << " is not a power of two";
463
2/4
✓ Branch 2 taken 3222 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3222 times.
✗ Branch 6 not taken.
 3222 throw std::runtime_error(ss.str());
464 3222 }
465
466 193336 SparseMatrixXcd get_sparse_matrix(
467 const std::vector<TripletCd> &triplets, unsigned rows, unsigned cols) {
468 193336 SparseMatrixXcd matr(rows, cols);
469
1/2
✓ Branch 3 taken 193336 times.
✗ Branch 4 not taken.
 193336 matr.setFromTriplets(triplets.cbegin(), triplets.cend());
470 193336 return matr;
471 }
472
473 192750 SparseMatrixXcd get_sparse_square_matrix(
474 const std::vector<TripletCd> &triplets, unsigned rows) {
475 192750 return get_sparse_matrix(triplets, rows, rows);
476 }
477
478 ✗ std::vector<TripletCd> get_triplets(
479 const SparseMatrixXcd &matr, double abs_epsilon) {
480 ✗ std::vector<TripletCd> result;
481 ✗ for (unsigned kk = 0; kk < matr.outerSize(); ++kk) {
482 ✗ for (SparseMatrixXcd::InnerIterator it(matr, kk); it; ++it) {
483 ✗ if (std::abs(it.value()) > abs_epsilon) {
484 ✗ const unsigned row = it.row();
485 ✗ const unsigned col = it.col();
486 ✗ result.emplace_back(row, col, it.value());
487 }
488 }
489 }
490 ✗ return result;
491 }
492
493 191431 std::vector<TripletCd> get_triplets(
494 const Eigen::MatrixXcd &matr, double abs_epsilon) {
495 191431 std::vector<TripletCd> triplets;
496
2/2
✓ Branch 1 taken 550257 times.
✓ Branch 2 taken 191431 times.
 741688 for (unsigned jj = 0; jj < matr.cols(); ++jj) {
497
2/2
✓ Branch 1 taken 1781185 times.
✓ Branch 2 taken 550257 times.
 2331442 for (unsigned ii = 0; ii < matr.rows(); ++ii) {
498
3/4
✓ Branch 1 taken 1781185 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 702533 times.
✓ Branch 5 taken 1078652 times.
 1781185 if (!(std::abs(matr(ii, jj)) <= abs_epsilon)) {
499
2/4
✓ Branch 1 taken 702533 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 702533 times.
✗ Branch 5 not taken.
 702533 triplets.emplace_back(ii, jj, matr(ii, jj));
500 }
501 }
502 }
503 191431 return triplets;
504 }
505
506 261 bool in_weyl_chamber(const std::array<Expr, 3> &k) {
507 261 bool is_symbolic = true;
508 261 double last_val = .5;
509
2/2
✓ Branch 1 taken 768 times.
✓ Branch 2 taken 246 times.
 1014 for (unsigned i = 0; i < k.size(); ++i) {
510
1/2
✓ Branch 2 taken 768 times.
✗ Branch 3 not taken.
 768 std::optional<double> eval = eval_expr_mod(k[i], 4);
511
2/2
✓ Branch 1 taken 728 times.
✓ Branch 2 taken 40 times.
 768 if (eval) {
512 728 is_symbolic = false;
513
2/2
✓ Branch 1 taken 244 times.
✓ Branch 2 taken 484 times.
 728 if (i + 1 == k.size()) {
514 244 double abs_eval = std::min(*eval, -(*eval) + 4);
515
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 241 times.
 244 if (abs_eval - last_val > EPS) {
516 15 return false;
517 }
518 } else {
519
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 473 times.
 484 if (*eval - last_val > EPS) {
520 11 return false;
521 }
522 }
523 714 last_val = *eval;
524
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 39 times.
 40 } else if (!is_symbolic) {
525 1 return false;
526 }
527 }
528 246 return true;
529 }
530
531 3096 Eigen::Matrix2cd nth_root(const Eigen::Matrix2cd &u, unsigned n) {
532
4/6
✓ Branch 1 taken 3096 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3096 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3094 times.
 3096 if (u.isApprox(Eigen::Matrix2cd::Identity(), EPS)) {
533
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 return Eigen::Matrix2cd::Identity();
534
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3093 times.
 3094 } else if (n == 0) {
535
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 throw std::invalid_argument("Non-identity matrix does not have a 0th root");
536 }
537
1/2
✓ Branch 1 taken 3093 times.
✗ Branch 2 not taken.
 3093 Eigen::ComplexEigenSolver<Eigen::Matrix2cd> eigen_solver(u);
538
3/6
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3093 times.
✗ Branch 9 not taken.
 3093 return std::pow(eigen_solver.eigenvalues()[0], 1. / n) *
539
2/4
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
 6186 eigen_solver.eigenvectors().col(0) *
540
3/6
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3093 times.
✗ Branch 9 not taken.
 6186 eigen_solver.eigenvectors().col(0).adjoint() +
541
3/6
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3093 times.
✗ Branch 9 not taken.
 3093 std::pow(eigen_solver.eigenvalues()[1], 1. / n) *
542
2/4
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
 6186 eigen_solver.eigenvectors().col(1) *
543
3/6
✓ Branch 2 taken 3093 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3093 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3093 times.
✗ Branch 9 not taken.
 9279 eigen_solver.eigenvectors().col(1).adjoint();
544 }
545
546 } // namespace tket
547