Quantum++  v1.0-rc2
A modern C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef INTERNAL_UTIL_H_
28 #define INTERNAL_UTIL_H_
29 
30 namespace qpp
31 {
36 namespace internal
37 {
38 // integer index to multi-index, use C-style array for speed
39 // standard lexicographical order, e.g. 00, 01, 10, 11
40 inline void n2multiidx(idx n, idx numdims, const idx* const dims, idx* result)
41 noexcept
42 {
43  // error checks only in DEBUG version
44 #ifndef NDEBUG
45  if (numdims > 0) // numdims equal zero is a no-op
46  {
47  idx D = 1;
48  for (idx i = 0; i < numdims; ++i)
49  D *= dims[i];
50  assert(n < D);
51  }
52 #endif
53  // no error checks in release version to improve speed
54  for (idx i = 0; i < numdims; ++i)
55  {
56  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
57  n /= (dims[numdims - i - 1]);
58  }
59 }
60 
61 // silence g++4.9 bogus warning -Warray-bounds and -Wmaybe-uninitialized
62 // in qpp::internal::multiidx2n()
63 #if (__GNUC__ && !__clang__)
64 #pragma GCC diagnostic push
65 #pragma GCC diagnostic ignored "-Warray-bounds"
66 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
67 #endif
68 // multi-index to integer index, use C-style array for speed,
69 // standard lexicographical order, e.g. 00->0, 01->1, 10->2, 11->3
70 inline idx multiidx2n(const idx* const midx, idx numdims, const idx* const dims)
71 noexcept
72 {
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  {
88  part_prod[numdims - i - 1] =
89  part_prod[numdims - i] * dims[numdims - i];
90  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
91  }
92 
93  return result + midx[numdims - 1];
94 }
95 #if (__GNUC__ && !__clang__)
96 #pragma GCC diagnostic pop
97 #endif
98 
99 // check square matrix
100 template<typename Derived>
101 bool check_square_mat(const Eigen::MatrixBase<Derived>& A)
102 {
103  return A.rows() == A.cols();
104 }
105 
106 // check whether input is a vector or not
107 template<typename Derived>
108 bool check_vector(const Eigen::MatrixBase<Derived>& A)
109 {
110  return A.rows() == 1 || A.cols() == 1;
111 }
112 
113 // check whether input is a row vector or not
114 template<typename Derived>
115 bool check_rvector(const Eigen::MatrixBase<Derived>& A)
116 {
117  return A.rows() == 1;
118 }
119 
120 // check whether input is a column vector or not
121 template<typename Derived>
122 bool check_cvector(const Eigen::MatrixBase<Derived>& A)
123 {
124  return A.cols() == 1;
125 }
126 
127 // check non-zero size of object that supports size() function
128 template<typename T>
129 bool check_nonzero_size(const T& x) noexcept
130 {
131  return x.size() != 0;
132 }
133 
134 // check that all sizes match
135 template<typename T1, typename T2>
136 bool check_matching_sizes(const T1& lhs, const T2& rhs) noexcept
137 {
138  return lhs.size() == rhs.size();
139 }
140 
141 // check that dims is a valid dimension vector
142 inline bool check_dims(const std::vector<idx>& dims)
143 {
144  if (dims.size() == 0)
145  return false;
146 
147  return std::find_if(std::begin(dims), std::end(dims),
148  [dims](idx i) -> bool
149  {
150  if (i == 0) return true;
151  else return false;
152  }) == std::end(dims);
153 }
154 
155 // check that valid dims match the dimensions
156 // of valid (non-zero sized) square matrix
157 template<typename Derived>
158 bool check_dims_match_mat(const std::vector<idx>& dims,
159  const Eigen::MatrixBase<Derived>& A)
160 {
161 // error checks only in DEBUG version
162 #ifndef NDEBUG
163  assert(dims.size() > 0);
164  assert(A.rows() == A.cols());
165 #endif
166  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
167  static_cast<idx>(1), std::multiplies<idx>());
168 
169  return proddim == static_cast<idx>(A.rows());
170 }
171 
172 // check that valid dims match the dimensions of valid column vector
173 template<typename Derived>
174 bool check_dims_match_cvect(const std::vector<idx>& dims,
175  const Eigen::MatrixBase<Derived>& A)
176 {
177 // error checks only in DEBUG version
178 #ifndef NDEBUG
179  assert(dims.size() > 0);
180  assert(A.rows() > 0);
181  assert(A.cols() == 1);
182 #endif
183  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
184  static_cast<idx>(1), std::multiplies<idx>());
185 
186  return proddim == static_cast<idx>(A.rows());
187 }
188 
189 // check that valid dims match the dimensions of valid row vector
190 template<typename Derived>
191 bool check_dims_match_rvect(const std::vector<idx>& dims,
192  const Eigen::MatrixBase<Derived>& A)
193 {
194 // error checks only in DEBUG version
195 #ifndef NDEBUG
196  assert(dims.size() > 0);
197  assert(A.cols() > 0);
198  assert(A.rows() == 1);
199 #endif
200  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
201  static_cast<idx>(1), std::multiplies<idx>());;
202 
203  return proddim == static_cast<idx>(A.cols());
204 }
205 
206 // check that all elements in valid dims equal to dim
207 inline bool check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept
208 {
209 // error checks only in DEBUG version
210 #ifndef NDEBUG
211  assert(dims.size() > 0);
212 #endif
213  for (idx i : dims)
214  if (i != dim)
215  return false;
216 
217  return true;
218 }
219 
220 // check that subsys is valid with respect to valid dims
221 inline bool check_subsys_match_dims(const std::vector<idx>& subsys,
222  const std::vector<idx>& dims)
223 {
224  // subsys can be empty
225 
226  // check valid number of subsystems
227  if (subsys.size() > dims.size())
228  return false;
229 
230  // sort the subsystems
231  std::vector<idx> subsyssort = subsys;
232  std::sort(std::begin(subsyssort), std::end(subsyssort));
233 
234  // check duplicates
235  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
236  != std::end(subsyssort))
237  return false;
238 
239  // check range of subsystems
240  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
241  [dims](idx i) -> bool
242  {
243  return i > dims.size() - 1;
244  }) == std::end(subsyssort);
245 }
246 
247 // check matrix is 2 x 2
248 template<typename Derived>
249 bool check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
250 {
251  return A.rows() == 2 && A.cols() == 2;
252 }
253 
254 // check column vector is 2 x 1
255 template<typename Derived>
256 bool check_qubit_cvector(const Eigen::MatrixBase<Derived>& A) noexcept
257 {
258  return A.rows() == 2 && A.cols() == 1;
259 }
260 
261 // check row vector is 1 x 2
262 template<typename Derived>
263 bool check_qubit_rvector(const Eigen::MatrixBase<Derived>& A) noexcept
264 {
265  return A.rows() == 1 && A.cols() == 2;
266 }
267 
268 // check row vector is 1 x 2 or 2 x 1
269 template<typename Derived>
270 bool check_qubit_vector(const Eigen::MatrixBase<Derived>& A) noexcept
271 {
272  return (A.rows() == 1 && A.cols() == 2) ||
273  (A.rows() == 2 && A.cols() == 1);
274 }
275 
276 
277 // check valid permutation
278 inline bool check_perm(const std::vector<idx>& perm)
279 {
280  if (perm.size() == 0)
281  return false;
282 
283  std::vector<idx> ordered(perm.size());
284  std::iota(std::begin(ordered), std::end(ordered), 0);
285 
286  return std::is_permutation(std::begin(ordered), std::end(ordered),
287  std::begin(perm));
288 }
289 
290 // Kronecker product of 2 matrices, preserve return type
291 // internal function for the variadic template function wrapper kron()
292 template<typename Derived1, typename Derived2>
293 dyn_mat<typename Derived1::Scalar> kron2(const Eigen::MatrixBase<Derived1>& A,
294  const Eigen::MatrixBase<Derived2>& B)
295 {
296  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
297  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
298 
299  // EXCEPTION CHECKS
300 
301  // check types
302  if (!std::is_same<typename Derived1::Scalar,
303  typename Derived2::Scalar>::value)
304  throw exception::TypeMismatch("qpp::kron()");
305 
306  // check zero-size
308  throw exception::ZeroSize("qpp::kron()");
309 
310  // check zero-size
312  throw exception::ZeroSize("qpp::kron()");
313  // END EXCEPTION CHECKS
314 
315  idx Acols = static_cast<idx>(rA.cols());
316  idx Arows = static_cast<idx>(rA.rows());
317  idx Bcols = static_cast<idx>(rB.cols());
318  idx Brows = static_cast<idx>(rB.rows());
319 
321  result.resize(Arows * Brows, Acols * Bcols);
322 
323 #ifdef WITH_OPENMP_
324 #pragma omp parallel for collapse(2)
325 #endif // WITH_OPENMP_
326  for (idx j = 0; j < Acols; ++j) // column major order for speed
327  for (idx i = 0; i < Arows; ++i)
328  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
329 
330  return result;
331 }
332 
333 // Direct sum of 2 matrices, preserve return type
334 // internal function for the variadic template function wrapper dirsum()
335 template<typename Derived1, typename Derived2>
337  const Eigen::MatrixBase<Derived1>& A,
338  const Eigen::MatrixBase<Derived2>& B)
339 {
340  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
341  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
342 
343  // EXCEPTION CHECKS
344 
345  // check types
346  if (!std::is_same<typename Derived1::Scalar,
347  typename Derived2::Scalar>::value)
348  throw exception::TypeMismatch("qpp::dirsum()");
349 
350  // check zero-size
352  throw exception::ZeroSize("qpp::dirsum()");
353 
354  // check zero-size
356  throw exception::ZeroSize("qpp::dirsum()");
357  // END EXCEPTION CHECKS
358 
359  idx Acols = static_cast<idx>(rA.cols());
360  idx Arows = static_cast<idx>(rA.rows());
361  idx Bcols = static_cast<idx>(rB.cols());
362  idx Brows = static_cast<idx>(rB.rows());
363 
366  Acols + Bcols);
367 
368  result.block(0, 0, Arows, Acols) = rA;
369  result.block(Arows, Acols, Brows, Bcols) = rB;
370 
371  return result;
372 }
373 
374 // may be useful, extracts variadic template argument pack into a std::vector
375 template<typename T>
376 // ends the recursion
377 void variadic_vector_emplace(std::vector<T>&)
378 {
379 }
380 
381 // may be useful, extracts variadic template argument pack into a std::vector
382 template<typename T, typename First, typename ... Args>
383 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&& ... args)
384 {
385  v.emplace_back(std::forward<First>(first));
386  variadic_vector_emplace(v, std::forward<Args>(args)...);
387 }
388 
389 // returns the number of subsystems (each subsystem assumed of the same
390 // dimension d) from an object (ket/bra/density matrix) of size sz
391 inline idx get_num_subsys(idx sz, idx d)
392 {
393 // error checks only in DEBUG version
394 #ifndef NDEBUG
395  assert(sz > 0);
396  assert(d > 1);
397 #endif
398  return static_cast<idx>(std::llround(std::log2(sz) / std::log2(d)));
399 }
400 
401 // returns the dimension of a subsystem (each subsystem assumed of the same
402 // dimension d) from an object (ket/bra/density matrix) of size sz consisting
403 // of N subsystems
404 inline idx get_dim_subsys(idx sz, idx N)
405 {
406 // error checks only in DEBUG version
407 #ifndef NDEBUG
408  assert(N > 0);
409  assert(sz > 0);
410 #endif
411  if (N == 2)
412  return static_cast<idx>(std::llround(std::sqrt(sz)));
413 
414  return static_cast<idx>(std::llround(std::pow(sz, 1. / N)));
415 }
416 
417 // implementation details for pretty formatting
419 {
420  template<typename T>
421  // T must support rows(), cols(), operator()(idx, idx) const
422  std::ostream& display_impl_(const T& A,
423  std::ostream& os,
424  double chop = qpp::chop) const
425  {
426  std::ostringstream ostr;
427  ostr.copyfmt(os); // copy os' state
428 
429  std::vector<std::string> vstr;
430  std::string strA;
431 
432  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
433  {
434  for (idx j = 0;
435  j < static_cast<idx>(A.cols()); ++j)
436  {
437  strA.clear(); // clear the temporary string
438  ostr.clear();
439  ostr.str(std::string {}); // clear the ostringstream
440 
441  // convert to complex
442  double re = static_cast<cplx>(A(i, j)).real();
443  double im = static_cast<cplx>(A(i, j)).imag();
444 
445  if (std::abs(re) < chop && std::abs(im) < chop)
446  {
447  ostr << "0 "; // otherwise segfault on destruction
448  // if using only vstr.push_back("0 ");
449  // bug in MATLAB libmx
450  vstr.push_back(ostr.str());
451  } else if (std::abs(re) < chop)
452  {
453  ostr << im;
454  vstr.push_back(ostr.str() + "i");
455  } else if (std::abs(im) < chop)
456  {
457  ostr << re;
458  vstr.push_back(ostr.str() + " ");
459  } else
460  {
461  ostr << re;
462  strA = ostr.str();
463 
464  strA += (im > 0 ? " + " : " - ");
465  ostr.clear();
466  ostr.str(std::string()); // clear
467  ostr << std::abs(im);
468  strA += ostr.str();
469  strA += "i";
470  vstr.push_back(strA);
471  }
472  }
473  }
474 
475  // determine the maximum lenght of the entries in each column
476  std::vector<idx> maxlengthcols(A.cols(), 0);
477 
478  for (idx i = 0; i < static_cast<idx>(A.rows());
479  ++i)
480  for (idx j = 0;
481  j < static_cast<idx>(A.cols()); ++j)
482  if (vstr[i * A.cols() + j].size() > maxlengthcols[j])
483  maxlengthcols[j] = vstr[i * A.cols() + j].size();
484 
485  // finally display it!
486  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
487  {
488  os << std::setw(static_cast<int>(maxlengthcols[0])) << std::right
489  << vstr[i * A.cols()]; // display first column
490  // then the rest
491  for (idx j = 1;
492  j < static_cast<idx>(A.cols()); ++j)
493  os << std::setw(static_cast<int>(maxlengthcols[j] + 2))
494  << std::right << vstr[i * A.cols() + j];
495 
496  if (i < static_cast<idx>(A.rows()) - 1)
497  os << std::endl;
498  }
499 
500  return os;
501  }
502 };
503 
504 } /* namespace internal */
505 } /* namespace qpp */
506 
507 #endif /* INTERNAL_UTIL_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
constexpr double chop
Used in qpp::disp() for setting to zero numbers that have their absolute value smaller than qpp::chop...
Definition: constants.h:56
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:71
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:221
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:293
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:122
std::ostream & display_impl_(const T &A, std::ostream &os, double chop=qpp::chop) const
Definition: util.h:422
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:77
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:158
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:70
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:40
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:278
idx get_num_subsys(idx sz, idx d)
Definition: util.h:391
bool check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:115
bool check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:207
bool check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:108
Definition: util.h:418
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:174
bool check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:191
bool check_matching_sizes(const T1 &lhs, const T2 &rhs) noexcept
Definition: util.h:136
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:45
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:377
Type mismatch exception.
Definition: exception.h:585
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:404
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:336
bool check_qubit_rvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:263
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:142
bool check_qubit_vector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:270
std::size_t idx
Non-negative integer index.
Definition: types.h:35
bool check_qubit_cvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:256
Object has zero size exception.
Definition: exception.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:249