GCC Code Coverage Report

Directory: ./
File: Simulation/PauliExpBoxUnitaryCalculator.cpp
Date: 2022-10-15 05:10:18
Exec Total Coverage
Lines: 76 78 97.4%
Functions: 9 9 100.0%
Branches: 40 70 57.1%
Decisions: 15 16 93.8%
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 "PauliExpBoxUnitaryCalculator.hpp"
16
17 #include <algorithm>
18 #include <tkassert/Assert.hpp>
19
20 #include "Circuit/Boxes.hpp"
21 #include "Utils/PauliStrings.hpp"
22
23 namespace tket {
24 namespace tket_sim {
25 namespace internal {
26
27 // A nonzero entry in the matrix containing only 1,0,-1 values,
28 // which we will build up for the tensor product
29 // (taking out a factor of i from Y).
30 typedef Eigen::Triplet<int> Entry;
31
32 1 static std::map<Pauli, std::array<Entry, 2>> get_pauli_map() {
33 1 std::map<Pauli, std::array<Entry, 2>> result;
34 {
35
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto& i_matr = result[Pauli::I];
36 1 i_matr[0] = Entry(0, 0, 1);
37 1 i_matr[1] = Entry(1, 1, 1);
38 }
39 {
40
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto& x_matr = result[Pauli::X];
41 1 x_matr[0] = Entry(0, 1, 1);
42 1 x_matr[1] = Entry(1, 0, 1);
43 }
44 {
45 // Remove a factor of i.
46
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto& y_matr = result[Pauli::Y];
47 1 y_matr[0] = Entry(0, 1, -1);
48 1 y_matr[1] = Entry(1, 0, 1);
49 }
50 {
51
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto& z_matr = result[Pauli::Z];
52 1 z_matr[0] = Entry(0, 0, 1);
53 1 z_matr[1] = Entry(1, 1, -1);
54 }
55 1 return result;
56 }
57
58 // The tensor product of the left and right entries.
59 45660 static Entry get_combined_entry(const Entry& left, const Entry& right) {
60 91320 return Entry(
61 45660 2 * left.row() + right.row(), 2 * left.col() + right.col(),
62 91320 left.value() * right.value());
63 }
64
65 namespace {
66 // To save repeated reallocation, build the data up inside here.
67 struct PauliExpBoxUnitaryCalculator {
68 // The two nonzero entries of the 2x2 matrix.
69 const std::map<Pauli, std::array<Entry, 2>> pauli_map;
70
71 // The tensor product matrix, all factors of i removed.
72 // Although QubitPauliString and QubitPauliTensor could probably
73 // be made to do the work for us, they are more complicated
74 // than we need.
75 std::vector<Entry> sparse_matrix;
76
77 // The number of Y which occurred.
78 unsigned power_of_i;
79
80 // A work vector to avoid repeated reallocation.
81 std::vector<TripletCd> triplets;
82
83 // We need to remember which diagonal entries were filled,
84 // to construct the final triplets.
85 std::vector<unsigned> set_diagonals;
86
87 // In the tensor product of the given entry in "sparse_matrix"
88 // with the new 2x2 Pauli matrix on the right,
89 // add the new entries to "sparse_matrix"
90 // (including changing the existing entry).
91 // I.e., each entry in "sparse_matrix" is a 1 or -1,
92 // which adds two more +/-1 values to the list of nonzero entries
93 // (whilst being removed itself),
94 // since we take out a factor of i from Y.
95 void add_entries(unsigned sparse_matrix_index, Pauli pauli);
96
97 // Call this to begin entering a new Pauli string "manually".
98 void clear();
99
100 // Add a single Pauli to the RIGHT of the current tensor product.
101 void append(Pauli pauli);
102
103 void fill_triplets(double phase);
104
105 PauliExpBoxUnitaryCalculator();
106
107 static PauliExpBoxUnitaryCalculator& get();
108 };
109
110 1845 PauliExpBoxUnitaryCalculator& PauliExpBoxUnitaryCalculator::get() {
111
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1844 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 1845 static PauliExpBoxUnitaryCalculator calculator;
112 1845 return calculator;
113 }
114
115 1 PauliExpBoxUnitaryCalculator::PauliExpBoxUnitaryCalculator()
116 1 : pauli_map(get_pauli_map()) {}
117
118 1845 void PauliExpBoxUnitaryCalculator::clear() {
119 1845 sparse_matrix.resize(1);
120 1845 sparse_matrix[0] = Entry(0, 0, 1);
121 1845 power_of_i = 0;
122 1845 triplets.clear();
123 1845 }
124
125 22830 void PauliExpBoxUnitaryCalculator::add_entries(
126 unsigned sparse_matrix_index, Pauli pauli) {
127 TKET_ASSERT(sparse_matrix_index < sparse_matrix.size());
128 22830 const auto& single_pauli = pauli_map.at(pauli);
129
1/2
✓ Branch 1 taken 22830 times.
✗ Branch 2 not taken.
 22830 sparse_matrix.push_back(
130 22830 get_combined_entry(sparse_matrix[sparse_matrix_index], single_pauli[0]));
131 // Overwrite existing entry.
132 // The reference is valid until after the new entry returns!
133 22830 auto& existing_entry = sparse_matrix[sparse_matrix_index];
134 22830 existing_entry = get_combined_entry(existing_entry, single_pauli[1]);
135 22830 }
136
137 6668 void PauliExpBoxUnitaryCalculator::append(Pauli pauli) {
138 6668 const auto current_size = sparse_matrix.size();
139
2/2
✓ Branch 0 taken 1662 times.
✓ Branch 1 taken 5006 times.
2/2
✓ Decision 'true' taken 1662 times.
✓ Decision 'false' taken 5006 times.
 6668 if (pauli == Pauli::Y) {
140 1662 ++power_of_i;
141 }
142
2/2
✓ Branch 0 taken 22830 times.
✓ Branch 1 taken 6668 times.
2/2
✓ Decision 'true' taken 22830 times.
✓ Decision 'false' taken 6668 times.
 29498 for (unsigned ii = 0; ii < current_size; ++ii) {
143 22830 add_entries(ii, pauli);
144 }
145 6668 }
146
147 1845 void PauliExpBoxUnitaryCalculator::fill_triplets(double phase) {
148
1/2
✓ Branch 2 taken 1845 times.
✗ Branch 3 not taken.
 1845 triplets.reserve(sparse_matrix.size());
149 1845 triplets.clear();
150
151 // If M^2 = I then exp(itM) = cos(t)I + i.sin(t)M.
152 1845 const double angle = -0.5 * PI * phase;
153 1845 const double cc = std::cos(angle);
154 1845 const double ss = std::sin(angle);
155
156 // TODO: efficiency: (i) check if cos(t) or sin(t) is nearly zero;
157 // (ii): fill the I coefficients FIRST, so we don't have to bother with
158 // "set_diagonals" (as we can overwrite diagonal entries directly).
159 1845 const std::complex<double> identity_coefficient(cc);
160
161 // We took out factors of i from M (coming from Pauli Y),
162 // so have to put them back in.
163 1845 power_of_i = power_of_i % 4;
164 const auto matrix_coefficient =
165 1845 std::polar(1.0, 0.5 * (power_of_i + 1) * PI) * ss;
166
167 // When we form the sparse representation of the matrix exp(itM),
168 // we must check which diagonal entries were missed out,
169 // coming from the cos(t)I term.
170 1845 set_diagonals.clear();
171
172
2/2
✓ Branch 5 taken 24675 times.
✓ Branch 6 taken 1845 times.
2/2
✓ Decision 'true' taken 24675 times.
✓ Decision 'false' taken 1845 times.
 26520 for (const auto& entry : sparse_matrix) {
173 // In Utils\PauliStrings.hpp, there is
174 // QubitPauliTensor operator*(Complex a, const QubitPauliTensor &qpt),
175 // which messes up ordinary std statements like std::complex * int.
176 // Not really a problem, but can be unexpected.
177 std::complex<double> value =
178
1/2
✓ Branch 3 taken 24675 times.
✗ Branch 4 not taken.
 24675 matrix_coefficient * std::complex<double>(entry.value());
179
180
2/2
✓ Branch 2 taken 1779 times.
✓ Branch 3 taken 22896 times.
2/2
✓ Decision 'true' taken 1779 times.
✓ Decision 'false' taken 22896 times.
 24675 if (entry.row() == entry.col()) {
181
1/2
✓ Branch 2 taken 1779 times.
✗ Branch 3 not taken.
 1779 set_diagonals.push_back(entry.row());
182 1779 value += identity_coefficient;
183 }
184
1/2
✓ Branch 3 taken 24675 times.
✗ Branch 4 not taken.
 24675 triplets.emplace_back(entry.row(), entry.col(), value);
185 }
186 // Complexity: let there be k diagonals set, out of N, and you want
187 // the other N-k (remembering that N = 2^n can be large).
188 // Here, we storing the diagonals we DON'T want in a std::vector, sort them,
189 // then select the OTHER elements. This takes time
190 // O(k.log k) + O(N.log k) time (with pretty small O constants,
191 // since sorting a std::vector is fast).
192 //
193 // If instead we use std::set to store all diagonals,
194 // remove each as it is set, and finally iterate through the
195 // REMAINING elements, it would be
196 // O(N.log N) + O((N-k).log N) + O(k.log k) time.
197 // This is probably a bit worse.
198 //
199 // If we really wanted to speed it up, we'd use bit operations.
200
1/2
✓ Branch 3 taken 1845 times.
✗ Branch 4 not taken.
 1845 std::sort(set_diagonals.begin(), set_diagonals.end());
201
2/2
✓ Branch 1 taken 24675 times.
✓ Branch 2 taken 1845 times.
2/2
✓ Decision 'true' taken 24675 times.
✓ Decision 'false' taken 1845 times.
 26520 for (unsigned ii = 0; ii < sparse_matrix.size(); ++ii) {
202
3/4
✓ Branch 3 taken 24675 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 22896 times.
✓ Branch 6 taken 1779 times.
2/2
✓ Decision 'true' taken 22896 times.
✓ Decision 'false' taken 1779 times.
 24675 if (!std::binary_search(set_diagonals.cbegin(), set_diagonals.cend(), ii)) {
203
1/2
✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
 22896 triplets.emplace_back(ii, ii, identity_coefficient);
204 }
205 }
206 1845 }
207 } // namespace
208
209 1845 std::vector<TripletCd> get_triplets(const PauliExpBox& box) {
210
2/4
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1845 times.
✗ Branch 5 not taken.
 1845 const auto phase_optional = eval_expr(box.get_phase());
211
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1845 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 1845 times.
 1845 if (!phase_optional) {
212 ✗ throw SymbolsNotSupported(
213 "PauliExpBoxUnitaryCalculator called "
214 ✗ "with symbolic phase parameter");
215 }
216
1/2
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
 1845 const auto pauli_string = box.get_paulis();
217
1/2
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
 1845 auto& calculator = PauliExpBoxUnitaryCalculator::get();
218
1/2
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
 1845 calculator.clear();
219
2/2
✓ Branch 5 taken 6668 times.
✓ Branch 6 taken 1845 times.
2/2
✓ Decision 'true' taken 6668 times.
✓ Decision 'false' taken 1845 times.
 8513 for (auto pauli : pauli_string) {
220
1/2
✓ Branch 1 taken 6668 times.
✗ Branch 2 not taken.
 6668 calculator.append(pauli);
221 }
222
2/4
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1845 times.
✗ Branch 5 not taken.
 1845 calculator.fill_triplets(phase_optional.value());
223
1/2
✓ Branch 1 taken 1845 times.
✗ Branch 2 not taken.
 3690 return calculator.triplets;
224 1845 }
225
226 } // namespace internal
227 } // namespace tket_sim
228 } // namespace tket
229