Quantum++  v1.0
A modern C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2018 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef INTERNAL_UTIL_H_
33 #define INTERNAL_UTIL_H_
34 
35 namespace qpp {
40 namespace internal {
41 // integer index to multi-index, use C-style array for speed
42 // standard lexicographical order, e.g. 00, 01, 10, 11
43 inline void n2multiidx(idx n, idx numdims, const idx* const dims,
44  idx* result) noexcept {
45 // error checks only in DEBUG version
46 #ifndef NDEBUG
47  if (numdims > 0) // numdims equal zero is a no-op
48  {
49  idx D = 1;
50  for (idx i = 0; i < numdims; ++i)
51  D *= dims[i];
52  assert(n < D);
53  }
54 #endif
55  // no error checks in release version to improve speed
56  for (idx i = 0; i < numdims; ++i) {
57  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
58  n /= (dims[numdims - i - 1]);
59  }
60 }
61 
62 // silence g++4.9 bogus warning -Warray-bounds and -Wmaybe-uninitialized
63 // in qpp::internal::multiidx2n()
64 #if (__GNUC__ && !__clang__)
65 #pragma GCC diagnostic push
66 #pragma GCC diagnostic ignored "-Warray-bounds"
67 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
68 #endif
69 // multi-index to integer index, use C-style array for speed,
70 // standard lexicographical order, e.g. 00->0, 01->1, 10->2, 11->3
71 inline idx multiidx2n(const idx* const midx, idx numdims,
72  const idx* const dims) noexcept {
73 // error checks only in DEBUG version
74 #ifndef NDEBUG
75  assert(numdims > 0);
76 #endif
77 
78  // no error checks in release version to improve speed
79 
80  // Static allocation for speed!
81  // double the size for matrices reshaped as vectors
82  idx part_prod[2 * maxn];
83 
84  idx result = 0;
85  part_prod[numdims - 1] = 1;
86  for (idx i = 1; i < numdims; ++i) {
87  part_prod[numdims - i - 1] = part_prod[numdims - i] * dims[numdims - i];
88  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
89  }
90 
91  return result + midx[numdims - 1];
92 }
93 #if (__GNUC__ && !__clang__)
94 #pragma GCC diagnostic pop
95 #endif
96 
97 // check square matrix
98 template <typename Derived>
99 bool check_square_mat(const Eigen::MatrixBase<Derived>& A) {
100  return A.rows() == A.cols();
101 }
102 
103 // check whether input is a vector or not
104 template <typename Derived>
105 bool check_vector(const Eigen::MatrixBase<Derived>& A) {
106  return A.rows() == 1 || A.cols() == 1;
107 }
108 
109 // check whether input is a row vector or not
110 template <typename Derived>
111 bool check_rvector(const Eigen::MatrixBase<Derived>& A) {
112  return A.rows() == 1;
113 }
114 
115 // check whether input is a column vector or not
116 template <typename Derived>
117 bool check_cvector(const Eigen::MatrixBase<Derived>& A) {
118  return A.cols() == 1;
119 }
120 
121 // check non-zero size of object that supports size() function
122 template <typename T>
123 bool check_nonzero_size(const T& x) noexcept {
124  return x.size() != 0;
125 }
126 
127 // check that all sizes match
128 template <typename T1, typename T2>
129 bool check_matching_sizes(const T1& lhs, const T2& rhs) noexcept {
130  return lhs.size() == rhs.size();
131 }
132 
133 // check that dims is a valid dimension vector
134 inline bool check_dims(const std::vector<idx>& dims) {
135  if (dims.size() == 0)
136  return false;
137 
138  return std::find_if(std::begin(dims), std::end(dims),
139  [dims](idx i) -> bool {
140  if (i == 0)
141  return true;
142  else
143  return false;
144  }) == std::end(dims);
145 }
146 
147 // check that valid dims match the dimensions
148 // of valid (non-zero sized) square matrix
149 template <typename Derived>
150 bool check_dims_match_mat(const std::vector<idx>& dims,
151  const Eigen::MatrixBase<Derived>& A) {
152 // error checks only in DEBUG version
153 #ifndef NDEBUG
154  assert(dims.size() > 0);
155  assert(A.rows() == A.cols());
156 #endif
157  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
158  static_cast<idx>(1), std::multiplies<idx>());
159 
160  return proddim == static_cast<idx>(A.rows());
161 }
162 
163 // check that valid dims match the dimensions of valid column vector
164 template <typename Derived>
165 bool check_dims_match_cvect(const std::vector<idx>& dims,
166  const Eigen::MatrixBase<Derived>& A) {
167 // error checks only in DEBUG version
168 #ifndef NDEBUG
169  assert(dims.size() > 0);
170  assert(A.rows() > 0);
171  assert(A.cols() == 1);
172 #endif
173  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
174  static_cast<idx>(1), std::multiplies<idx>());
175 
176  return proddim == static_cast<idx>(A.rows());
177 }
178 
179 // check that valid dims match the dimensions of valid row vector
180 template <typename Derived>
181 bool check_dims_match_rvect(const std::vector<idx>& dims,
182  const Eigen::MatrixBase<Derived>& A) {
183 // error checks only in DEBUG version
184 #ifndef NDEBUG
185  assert(dims.size() > 0);
186  assert(A.cols() > 0);
187  assert(A.rows() == 1);
188 #endif
189  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
190  static_cast<idx>(1), std::multiplies<idx>());
191  ;
192 
193  return proddim == static_cast<idx>(A.cols());
194 }
195 
196 // check that all elements in valid dims equal to dim
197 inline bool check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept {
198 // error checks only in DEBUG version
199 #ifndef NDEBUG
200  assert(dims.size() > 0);
201 #endif
202  for (idx i : dims)
203  if (i != dim)
204  return false;
205 
206  return true;
207 }
208 
209 // check that subsys is valid with respect to valid dims
210 inline bool check_subsys_match_dims(const std::vector<idx>& subsys,
211  const std::vector<idx>& dims) {
212  // subsys can be empty
213 
214  // check valid number of subsystems
215  if (subsys.size() > dims.size())
216  return false;
217 
218  // sort the subsystems
219  std::vector<idx> subsyssort = subsys;
220  std::sort(std::begin(subsyssort), std::end(subsyssort));
221 
222  // check duplicates
223  if (std::unique(std::begin(subsyssort), std::end(subsyssort)) !=
224  std::end(subsyssort))
225  return false;
226 
227  // check range of subsystems
228  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
229  [dims](idx i) -> bool {
230  return i > dims.size() - 1;
231  }) == std::end(subsyssort);
232 }
233 
234 // check matrix is 2 x 2
235 template <typename Derived>
236 bool check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept {
237  return A.rows() == 2 && A.cols() == 2;
238 }
239 
240 // check column vector is 2 x 1
241 template <typename Derived>
242 bool check_qubit_cvector(const Eigen::MatrixBase<Derived>& A) noexcept {
243  return A.rows() == 2 && A.cols() == 1;
244 }
245 
246 // check row vector is 1 x 2
247 template <typename Derived>
248 bool check_qubit_rvector(const Eigen::MatrixBase<Derived>& A) noexcept {
249  return A.rows() == 1 && A.cols() == 2;
250 }
251 
252 // check row vector is 1 x 2 or 2 x 1
253 template <typename Derived>
254 bool check_qubit_vector(const Eigen::MatrixBase<Derived>& A) noexcept {
255  return (A.rows() == 1 && A.cols() == 2) || (A.rows() == 2 && A.cols() == 1);
256 }
257 
258 // check valid permutation
259 inline bool check_perm(const std::vector<idx>& perm) {
260  if (perm.size() == 0)
261  return false;
262 
263  std::vector<idx> ordered(perm.size());
264  std::iota(std::begin(ordered), std::end(ordered), 0);
265 
266  return std::is_permutation(std::begin(ordered), std::end(ordered),
267  std::begin(perm));
268 }
269 
270 // Kronecker product of 2 matrices, preserve return type
271 // internal function for the variadic template function wrapper kron()
272 template <typename Derived1, typename Derived2>
273 dyn_mat<typename Derived1::Scalar> kron2(const Eigen::MatrixBase<Derived1>& A,
274  const Eigen::MatrixBase<Derived2>& B) {
275  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
276  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
277 
278  // EXCEPTION CHECKS
279 
280  // check types
281  if (!std::is_same<typename Derived1::Scalar,
282  typename Derived2::Scalar>::value)
283  throw exception::TypeMismatch("qpp::kron()");
284 
285  // check zero-size
287  throw exception::ZeroSize("qpp::kron()");
288 
289  // check zero-size
291  throw exception::ZeroSize("qpp::kron()");
292  // END EXCEPTION CHECKS
293 
294  idx Acols = static_cast<idx>(rA.cols());
295  idx Arows = static_cast<idx>(rA.rows());
296  idx Bcols = static_cast<idx>(rB.cols());
297  idx Brows = static_cast<idx>(rB.rows());
298 
300  result.resize(Arows * Brows, Acols * Bcols);
301 
302 #ifdef WITH_OPENMP_
303 #pragma omp parallel for collapse(2)
304 #endif // WITH_OPENMP_
305  // column major order for speed
306  for (idx j = 0; j < Acols; ++j)
307  for (idx i = 0; i < Arows; ++i)
308  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
309 
310  return result;
311 }
312 
313 // Direct sum of 2 matrices, preserve return type
314 // internal function for the variadic template function wrapper dirsum()
315 template <typename Derived1, typename Derived2>
317 dirsum2(const Eigen::MatrixBase<Derived1>& A,
318  const Eigen::MatrixBase<Derived2>& B) {
319  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
320  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
321 
322  // EXCEPTION CHECKS
323 
324  // check types
325  if (!std::is_same<typename Derived1::Scalar,
326  typename Derived2::Scalar>::value)
327  throw exception::TypeMismatch("qpp::dirsum()");
328 
329  // check zero-size
331  throw exception::ZeroSize("qpp::dirsum()");
332 
333  // check zero-size
335  throw exception::ZeroSize("qpp::dirsum()");
336  // END EXCEPTION CHECKS
337 
338  idx Acols = static_cast<idx>(rA.cols());
339  idx Arows = static_cast<idx>(rA.rows());
340  idx Bcols = static_cast<idx>(rB.cols());
341  idx Brows = static_cast<idx>(rB.rows());
342 
344  dyn_mat<typename Derived1::Scalar>::Zero(Arows + Brows, Acols + Bcols);
345 
346  result.block(0, 0, Arows, Acols) = rA;
347  result.block(Arows, Acols, Brows, Bcols) = rB;
348 
349  return result;
350 }
351 
352 // may be useful, extracts variadic template argument pack into a std::vector
353 template <typename T>
354 // ends the recursion
355 void variadic_vector_emplace(std::vector<T>&) {}
356 
357 // may be useful, extracts variadic template argument pack into a std::vector
358 template <typename T, typename First, typename... Args>
359 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&&... args) {
360  v.emplace_back(std::forward<First>(first));
361  variadic_vector_emplace(v, std::forward<Args>(args)...);
362 }
363 
364 // returns the number of subsystems (each subsystem assumed of the same
365 // dimension d) from an object (ket/bra/density matrix) of size sz
366 inline idx get_num_subsys(idx sz, idx d) {
367 // error checks only in DEBUG version
368 #ifndef NDEBUG
369  assert(sz > 0);
370  assert(d > 1);
371 #endif
372  return static_cast<idx>(std::llround(std::log2(sz) / std::log2(d)));
373 }
374 
375 // returns the dimension of a subsystem (each subsystem assumed of the same
376 // dimension d) from an object (ket/bra/density matrix) of size sz consisting
377 // of N subsystems
378 inline idx get_dim_subsys(idx sz, idx N) {
379 // error checks only in DEBUG version
380 #ifndef NDEBUG
381  assert(N > 0);
382  assert(sz > 0);
383 #endif
384  if (N == 2)
385  return static_cast<idx>(std::llround(std::sqrt(sz)));
386 
387  return static_cast<idx>(std::llround(std::pow(sz, 1. / N)));
388 }
389 
390 // implementation details for pretty formatting
392  template <typename T>
393  // T must support rows(), cols(), operator()(idx, idx) const
394  std::ostream& display_impl_(const T& A, std::ostream& os,
395  double chop = qpp::chop) const {
396  std::ostringstream ostr;
397  ostr.copyfmt(os); // copy os' state
398 
399  std::vector<std::string> vstr;
400  std::string strA;
401 
402  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
403  for (idx j = 0; j < static_cast<idx>(A.cols()); ++j) {
404  strA.clear(); // clear the temporary string
405  ostr.clear();
406  ostr.str(std::string{}); // clear the ostringstream
407 
408  // convert to complex
409  double re = static_cast<cplx>(A(i, j)).real();
410  double im = static_cast<cplx>(A(i, j)).imag();
411 
412  if (std::abs(re) < chop && std::abs(im) < chop) {
413  ostr << "0 "; // otherwise segfault on destruction
414  // if using only vstr.push_back("0 ");
415  // bug in MATLAB libmx
416  vstr.push_back(ostr.str());
417  } else if (std::abs(re) < chop) {
418  ostr << im;
419  vstr.push_back(ostr.str() + "i");
420  } else if (std::abs(im) < chop) {
421  ostr << re;
422  vstr.push_back(ostr.str() + " ");
423  } else {
424  ostr << re;
425  strA = ostr.str();
426 
427  strA += (im > 0 ? " + " : " - ");
428  ostr.clear();
429  ostr.str(std::string()); // clear
430  ostr << std::abs(im);
431  strA += ostr.str();
432  strA += "i";
433  vstr.push_back(strA);
434  }
435  }
436  }
437 
438  // determine the maximum lenght of the entries in each column
439  std::vector<idx> maxlengthcols(A.cols(), 0);
440 
441  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
442  for (idx j = 0; j < static_cast<idx>(A.cols()); ++j)
443  if (vstr[i * A.cols() + j].size() > maxlengthcols[j])
444  maxlengthcols[j] = vstr[i * A.cols() + j].size();
445 
446  // finally display it!
447  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
448  os << std::setw(static_cast<int>(maxlengthcols[0])) << std::right
449  << vstr[i * A.cols()]; // display first column
450  // then the rest
451  for (idx j = 1; j < static_cast<idx>(A.cols()); ++j)
452  os << std::setw(static_cast<int>(maxlengthcols[j] + 2))
453  << std::right << vstr[i * A.cols() + j];
454 
455  if (i < static_cast<idx>(A.rows()) - 1)
456  os << std::endl;
457  }
458 
459  return os;
460  }
461 };
462 
463 } /* namespace internal */
464 } /* namespace qpp */
465 
466 #endif /* INTERNAL_UTIL_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
constexpr double chop
Used in qpp::disp() for setting to zero numbers that have their absolute value smaller than qpp::chop...
Definition: constants.h:60
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:75
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:273
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
std::ostream & display_impl_(const T &A, std::ostream &os, double chop=qpp::chop) const
Definition: util.h:394
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
Quantum++ main namespace.
Definition: codes.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:259
idx get_num_subsys(idx sz, idx d)
Definition: util.h:366
bool check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:111
bool check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:197
bool check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
Definition: util.h:391
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
bool check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:181
bool check_matching_sizes(const T1 &lhs, const T2 &rhs) noexcept
Definition: util.h:129
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:355
Type mismatch exception.
Definition: exception.h:530
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:378
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:317
bool check_qubit_rvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:248
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
bool check_qubit_vector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:254
std::size_t idx
Non-negative integer index.
Definition: types.h:39
bool check_qubit_cvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:242
Object has zero size exception.
Definition: exception.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:236