XACC
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 - 2020 Vlad Gheorghiu.
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.empty())
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 vector has no duplicates
210 inline bool check_no_duplicates(std::vector<idx> v) {
211  std::sort(std::begin(v), std::end(v));
212  if (std::unique(std::begin(v), std::end(v)) != std::end(v))
213  return false;
214  else
215  return true;
216 }
217 
218 // check that subsys is valid with respect to valid dims
219 inline bool check_subsys_match_dims(const std::vector<idx>& subsys,
220  const std::vector<idx>& dims) {
221  // subsys can be empty
222 
223  // check valid number of subsystems
224  if (subsys.size() > dims.size())
225  return false;
226 
227  // check no duplicates
228  if (!check_no_duplicates(subsys))
229  return false;
230 
231  // check range of subsystems
232  return std::find_if(std::begin(subsys), std::end(subsys),
233  [dims](idx i) -> bool {
234  return i + 1 > dims.size();
235  }) == std::end(subsys);
236 }
237 
238 // check matrix is 2 x 2
239 template <typename Derived>
240 bool check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept {
241  return A.rows() == 2 && A.cols() == 2;
242 }
243 
244 // check column vector is 2 x 1
245 template <typename Derived>
246 bool check_qubit_cvector(const Eigen::MatrixBase<Derived>& A) noexcept {
247  return A.rows() == 2 && A.cols() == 1;
248 }
249 
250 // check row vector is 1 x 2
251 template <typename Derived>
252 bool check_qubit_rvector(const Eigen::MatrixBase<Derived>& A) noexcept {
253  return A.rows() == 1 && A.cols() == 2;
254 }
255 
256 // check row vector is 1 x 2 or 2 x 1
257 template <typename Derived>
258 bool check_qubit_vector(const Eigen::MatrixBase<Derived>& A) noexcept {
259  return (A.rows() == 1 && A.cols() == 2) || (A.rows() == 2 && A.cols() == 1);
260 }
261 
262 // check valid permutation
263 inline bool check_perm(const std::vector<idx>& perm) {
264  if (perm.empty())
265  return false;
266 
267  std::vector<idx> ordered(perm.size());
268  std::iota(std::begin(ordered), std::end(ordered), 0);
269 
270  return std::is_permutation(std::begin(ordered), std::end(ordered),
271  std::begin(perm));
272 }
273 
274 // Kronecker product of 2 matrices, preserve return type
275 // internal function for the variadic template function wrapper kron()
276 template <typename Derived1, typename Derived2>
277 dyn_mat<typename Derived1::Scalar> kron2(const Eigen::MatrixBase<Derived1>& A,
278  const Eigen::MatrixBase<Derived2>& B) {
279  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
280  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
281 
282  // EXCEPTION CHECKS
283 
284  // check types
285  if (!std::is_same<typename Derived1::Scalar,
286  typename Derived2::Scalar>::value)
287  throw exception::TypeMismatch("qpp::kron()");
288 
289  // check zero-size
290  if (!internal::check_nonzero_size(rA))
291  throw exception::ZeroSize("qpp::kron()");
292 
293  // check zero-size
294  if (!internal::check_nonzero_size(rB))
295  throw exception::ZeroSize("qpp::kron()");
296  // END EXCEPTION CHECKS
297 
298  idx Acols = static_cast<idx>(rA.cols());
299  idx Arows = static_cast<idx>(rA.rows());
300  idx Bcols = static_cast<idx>(rB.cols());
301  idx Brows = static_cast<idx>(rB.rows());
302 
303  dyn_mat<typename Derived1::Scalar> result;
304  result.resize(Arows * Brows, Acols * Bcols);
305 
306 #ifdef WITH_OPENMP_
307 #pragma omp parallel for collapse(2)
308 #endif // WITH_OPENMP_
309  // column major order for speed
310  for (idx j = 0; j < Acols; ++j)
311  for (idx i = 0; i < Arows; ++i)
312  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
313 
314  return result;
315 }
316 
317 // Direct sum of 2 matrices, preserve return type
318 // internal function for the variadic template function wrapper dirsum()
319 template <typename Derived1, typename Derived2>
320 dyn_mat<typename Derived1::Scalar>
321 dirsum2(const Eigen::MatrixBase<Derived1>& A,
322  const Eigen::MatrixBase<Derived2>& B) {
323  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
324  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
325 
326  // EXCEPTION CHECKS
327 
328  // check types
329  if (!std::is_same<typename Derived1::Scalar,
330  typename Derived2::Scalar>::value)
331  throw exception::TypeMismatch("qpp::dirsum()");
332 
333  // check zero-size
334  if (!internal::check_nonzero_size(rA))
335  throw exception::ZeroSize("qpp::dirsum()");
336 
337  // check zero-size
338  if (!internal::check_nonzero_size(rB))
339  throw exception::ZeroSize("qpp::dirsum()");
340  // END EXCEPTION CHECKS
341 
342  idx Acols = static_cast<idx>(rA.cols());
343  idx Arows = static_cast<idx>(rA.rows());
344  idx Bcols = static_cast<idx>(rB.cols());
345  idx Brows = static_cast<idx>(rB.rows());
346 
347  dyn_mat<typename Derived1::Scalar> result =
348  dyn_mat<typename Derived1::Scalar>::Zero(Arows + Brows, Acols + Bcols);
349 
350  result.block(0, 0, Arows, Acols) = rA;
351  result.block(Arows, Acols, Brows, Bcols) = rB;
352 
353  return result;
354 }
355 
356 // may be useful, extracts variadic template argument pack into a std::vector
357 template <typename T>
358 // ends the recursion
359 void variadic_vector_emplace(std::vector<T>&) {}
360 
361 // may be useful, extracts variadic template argument pack into a std::vector
362 template <typename T, typename First, typename... Args>
363 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&&... args) {
364  v.emplace_back(std::forward<First>(first));
365  variadic_vector_emplace(v, std::forward<Args>(args)...);
366 }
367 
368 // returns the number of subsystems (each subsystem assumed of the same
369 // dimension d) from an object (ket/bra/density matrix) of size D
370 inline idx get_num_subsys(idx D, idx d) {
371 // error checks only in DEBUG version
372 #ifndef NDEBUG
373  assert(D > 0);
374  assert(d > 1);
375 #endif
376  return static_cast<idx>(std::llround(std::log2(D) / std::log2(d)));
377 }
378 
379 // returns the dimension of a subsystem (each subsystem assumed of the same
380 // dimension d) from an object (ket/bra/density matrix) of size sz consisting
381 // of N subsystems
382 inline idx get_dim_subsys(idx sz, idx N) {
383 // error checks only in DEBUG version
384 #ifndef NDEBUG
385  assert(N > 0);
386  assert(sz > 0);
387 #endif
388  if (N == 2)
389  return static_cast<idx>(std::llround(std::sqrt(sz)));
390 
391  return static_cast<idx>(std::llround(std::pow(sz, 1. / N)));
392 }
393 
394 // chops a floating point or complex number to zero
395 template <typename T,
396  typename std::enable_if<std::numeric_limits<T>::is_iec559 ||
397  is_complex<T>::value>::type* = nullptr>
398 T abs_chop(const T& x, double chop = qpp::chop) {
399  if (std::abs(x) < chop)
400  return 0;
401 
402  return x;
403 }
404 
405 // returns it unchanged otherwise
406 template <typename T,
407  typename std::enable_if<!(std::numeric_limits<T>::is_iec559 ||
408  is_complex<T>::value)>::type* = nullptr>
409 T abs_chop(const T& x, double QPP_UNUSED_ chop = qpp::chop) {
410  return x;
411 }
412 
413 // implementation details for pretty formatting
414 struct Display_Impl_ {
415  template <typename T>
416  // T must support rows(), cols(), operator()(idx, idx) const
417  std::ostream& display_impl_(const T& A, std::ostream& os,
418  double chop = qpp::chop) const {
419  std::ostringstream ostr;
420  ostr.copyfmt(os); // copy os' state
421 
422  std::vector<std::string> vstr;
423  std::string strA;
424 
425  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
426  for (idx j = 0; j < static_cast<idx>(A.cols()); ++j) {
427  strA.clear(); // clear the temporary string
428  ostr.clear();
429  ostr.str(std::string{}); // clear the ostringstream
430 
431  // convert to complex
432  double re = static_cast<cplx>(A(i, j)).real();
433  double im = static_cast<cplx>(A(i, j)).imag();
434 
435  if (std::abs(re) < chop && std::abs(im) < chop) {
436  ostr << "0 "; // otherwise segfault on destruction
437  // if using only vstr.emplace_back("0 ");
438  // bug in MATLAB libmx
439  vstr.emplace_back(ostr.str());
440  } else if (std::abs(re) < chop) {
441  ostr << im;
442  vstr.emplace_back(ostr.str() + "i");
443  } else if (std::abs(im) < chop) {
444  ostr << re;
445  vstr.emplace_back(ostr.str() + " ");
446  } else {
447  ostr << re;
448  strA = ostr.str();
449 
450  strA += (im > 0 ? " + " : " - ");
451  ostr.clear();
452  ostr.str(std::string()); // clear
453  ostr << std::abs(im);
454  strA += ostr.str();
455  strA += "i";
456  vstr.emplace_back(strA);
457  }
458  }
459  }
460 
461  // determine the maximum lenght of the entries in each column
462  std::vector<idx> maxlengthcols(A.cols(), 0);
463 
464  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
465  for (idx j = 0; j < static_cast<idx>(A.cols()); ++j)
466  if (vstr[i * A.cols() + j].size() > maxlengthcols[j])
467  maxlengthcols[j] = vstr[i * A.cols() + j].size();
468 
469  // finally display it!
470  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
471  os << std::setw(static_cast<int>(maxlengthcols[0])) << std::right
472  << vstr[i * A.cols()]; // display first column
473  // then the rest
474  for (idx j = 1; j < static_cast<idx>(A.cols()); ++j)
475  os << std::setw(static_cast<int>(maxlengthcols[j] + 2))
476  << std::right << vstr[i * A.cols() + j];
477 
478  if (i < static_cast<idx>(A.rows()) - 1)
479  os << std::endl;
480  }
481 
482  return os;
483  }
484 };
485 
486 } /* namespace internal */
487 } /* namespace qpp */
488 
489 #endif /* INTERNAL_UTIL_H_ */
Quantum++ main namespace.
Definition: circuits.h:35