Quantum++  v1.0.0-beta3
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 // silence g++4.9 bogus warning -Warray-bounds and -Wmaybe-uninitialized
31 // in qpp::util::_multiidx2n()
32 #if (__GNUC__ && !__clang__)
33 #pragma GCC diagnostic ignored "-Warray-bounds"
34 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
35 #endif
36 
37 namespace qpp
38 {
43 namespace internal
44 {
45 
46 // integer index to multi-index, use C-style array for speed
47 // standard lexicographical order, e.g. 00, 01, 10, 11
48 inline void n2multiidx(idx n, idx numdims, const idx* const dims, idx* result)
49 noexcept
50 {
51  // error checks only in DEBUG version
52 #ifndef NDEBUG
53  if(numdims > 0) // numdims equal zero is a no-op
54  {
55  idx D = 1;
56  for (idx i = 0; i < numdims; ++i)
57  D *= dims[i];
58  assert(n < D);
59  }
60 #endif
61  // no error checks in release version to improve speed
62  for (idx i = 0; i < numdims; ++i)
63  {
64  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
65  n /= (dims[numdims - i - 1]);
66  }
67 }
68 
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, const idx* const dims)
72 noexcept
73 {
74  // error checks only in DEBUG version
75 #ifndef NDEBUG
76  assert(numdims > 0);
77 #endif
78 
79  // no error checks in release version to improve speed
80 
81  // Static allocation for speed!
82  // double the size for matrices reshaped as vectors
83  idx part_prod[2 * maxn];
84 
85  idx result = 0;
86  part_prod[numdims - 1] = 1;
87  for (idx i = 1; i < numdims; ++i)
88  {
89  part_prod[numdims - i - 1] =
90  part_prod[numdims - i] * dims[numdims - i];
91  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
92  }
93 
94  return result + midx[numdims - 1];
95 }
96 
97 // check square matrix
98 template<typename Derived>
99 bool check_square_mat(const Eigen::MatrixBase<Derived>& A)
100 {
101  return A.rows() == A.cols();
102 }
103 
104 // check whether input is a vector or not
105 template<typename Derived>
106 bool check_vector(const Eigen::MatrixBase<Derived>& A)
107 {
108  return A.rows() == 1 || A.cols() == 1;
109 }
110 
111 // check whether input is a row vector or not
112 template<typename Derived>
113 bool check_rvector(const Eigen::MatrixBase<Derived>& A)
114 {
115  return A.rows() == 1;
116 }
117 
118 // check whether input is a column vector or not
119 template<typename Derived>
120 bool check_cvector(const Eigen::MatrixBase<Derived>& A)
121 {
122  return A.cols() == 1;
123 }
124 
125 // check non-zero size of object that supports size() function
126 template<typename T>
127 bool check_nonzero_size(const T& x) noexcept
128 {
129  return x.size() != 0;
130 }
131 
132 // check that all sizes match
133 template<typename T1, typename T2>
134 bool check_matching_sizes(const T1& lhs, const T2& rhs) noexcept
135 {
136  return lhs.size() == rhs.size();
137 }
138 
139 // check that dims is a valid dimension vector
140 inline bool check_dims(const std::vector<idx>& dims)
141 {
142  if (dims.size() == 0)
143  return false;
144 
145  return std::find_if(std::begin(dims), std::end(dims),
146  [dims](idx i) -> bool
147  {
148  if (i == 0) return true;
149  else return false;
150  }) == std::end(dims);
151 }
152 
153 // check that valid dims match the dimensions
154 // of valid (non-zero sized) square matrix
155 template<typename Derived>
156 bool check_dims_match_mat(const std::vector<idx>& dims,
157  const Eigen::MatrixBase<Derived>& A)
158 {
159  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
160  static_cast<idx>(1), std::multiplies<idx>());
161 
162  return proddim == static_cast<idx>(A.rows());
163 }
164 
165 // check that valid dims match the dimensions of valid column vector
166 template<typename Derived>
167 bool check_dims_match_cvect(const std::vector<idx>& dims,
168  const Eigen::MatrixBase<Derived>& A)
169 {
170  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
171  static_cast<idx>(1), std::multiplies<idx>());
172 
173  return proddim == static_cast<idx>(A.rows());
174 }
175 
176 // check that valid dims match the dimensions of valid row vector
177 template<typename Derived>
178 bool check_dims_match_rvect(const std::vector<idx>& dims,
179  const Eigen::MatrixBase<Derived>& A)
180 {
181  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
182  static_cast<idx>(1), std::multiplies<idx>());;
183 
184  return proddim == static_cast<idx>(A.cols());
185 }
186 
187 // check that all elements in valid dims equal to dim
188 inline bool check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept
189 {
190  for (idx i : dims)
191  if (i != dim)
192  return false;
193 
194  return true;
195 }
196 
197 // check that subsys is valid with respect to valid dims
198 inline bool check_subsys_match_dims(const std::vector<idx>& subsys,
199  const std::vector<idx>& dims)
200 {
201  // // check non-zero sized subsystems
202  // if (subsys.size() == 0)
203  // return false;
204 
205  // check valid number of subsystems
206  if (subsys.size() > dims.size())
207  return false;
208 
209  // sort the subsystems
210  std::vector<idx> subsyssort = subsys;
211  std::sort(std::begin(subsyssort), std::end(subsyssort));
212 
213  // check duplicates
214  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
215  != std::end(subsyssort))
216  return false;
217 
218  // check range of subsystems
219  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
220  [dims](idx i) -> bool
221  {
222  return i > dims.size() - 1;
223  }) == std::end(subsyssort);
224 }
225 
226 // check matrix is 2 x 2
227 template<typename Derived>
228 bool check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
229 {
230  return A.rows() == 2 && A.cols() == 2;
231 }
232 
233 // check column vector is 2 x 1
234 template<typename Derived>
235 bool check_qubit_cvector(const Eigen::MatrixBase<Derived>& A) noexcept
236 {
237  return A.rows() == 2 && A.cols() == 1;
238 }
239 
240 // check row vector is 1 x 2
241 template<typename Derived>
242 bool check_qubit_rvector(const Eigen::MatrixBase<Derived>& A) noexcept
243 {
244  return A.rows() == 1 && A.cols() == 2;
245 }
246 
247 // check row vector is 1 x 2 or 2 x 1
248 template<typename Derived>
249 bool check_qubit_vector(const Eigen::MatrixBase<Derived>& A) noexcept
250 {
251  return (A.rows() == 1 && A.cols() == 2) ||
252  (A.rows() == 2 && A.cols() == 1);
253 }
254 
255 
256 // check valid permutation
257 inline bool check_perm(const std::vector<idx>& perm)
258 {
259  if (perm.size() == 0)
260  return false;
261 
262  std::vector<idx> ordered(perm.size());
263  std::iota(std::begin(ordered), std::end(ordered), 0);
264 
265  return std::is_permutation(std::begin(ordered), std::end(ordered),
266  std::begin(perm));
267 }
268 
269 // Kronecker product of 2 matrices, preserve return type
270 // internal function for the variadic template function wrapper kron()
271 template<typename Derived1, typename Derived2>
272 dyn_mat<typename Derived1::Scalar> kron2(const Eigen::MatrixBase<Derived1>& A,
273  const Eigen::MatrixBase<Derived2>& B)
274 {
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("qpp::kron()", Exception::Type::TYPE_MISMATCH);
284 
285  // check zero-size
287  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
288 
289  // check zero-size
291  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
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  for (idx j = 0; j < Acols; ++j) // column major order for speed
306  for (idx i = 0; i < Arows; ++i)
307  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
308 
309  return result;
310 }
311 
312 // Direct sum of 2 matrices, preserve return type
313 // internal function for the variadic template function wrapper dirsum()
314 template<typename Derived1, typename Derived2>
316  const Eigen::MatrixBase<Derived1>& A,
317  const Eigen::MatrixBase<Derived2>& B)
318 {
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("qpp::dirsum()", Exception::Type::TYPE_MISMATCH);
328 
329  // check zero-size
331  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
332 
333  // check zero-size
335  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
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 
345  Acols + Bcols);
346 
347  result.block(0, 0, Arows, Acols) = rA;
348  result.block(Arows, Acols, Brows, Bcols) = rB;
349 
350  return result;
351 }
352 
353 // may be useful, extracts variadic template argument pack into a std::vector
354 template<typename T>
355 // ends the recursion
356 void variadic_vector_emplace(std::vector<T>&)
357 {
358 }
359 
360 // may be useful, extracts variadic template argument pack into a std::vector
361 template<typename T, typename First, typename ... Args>
362 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&& ... args)
363 {
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 sz
370 inline idx get_num_subsys(idx sz, idx d)
371 {
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 {
380  if (N == 2)
381  return static_cast<idx>(std::llround(std::sqrt(sz)));
382 
383  return static_cast<idx>(std::llround(std::pow(sz, 1. / N)));
384 }
385 
386 // implementation details for pretty formatting
388 {
389  template<typename T>
390  // T must support rows(), cols(), operator()(idx, idx) const
391  std::ostream& display_impl_(const T& A,
392  std::ostream& os,
393  double chop = qpp::chop) const
394  {
395  std::ostringstream ostr;
396  ostr.copyfmt(os); // copy os' state
397 
398  std::vector<std::string> vstr;
399  std::string strA;
400 
401  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
402  {
403  for (idx j = 0;
404  j < static_cast<idx>(A.cols()); ++j)
405  {
406  strA.clear(); // clear the temporary string
407  ostr.clear();
408  ostr.str(std::string {}); // clear the ostringstream
409 
410  // convert to complex
411  double re = static_cast<cplx>(A(i, j)).real();
412  double im = static_cast<cplx>(A(i, j)).imag();
413 
414  if (std::abs(re) < chop && std::abs(im) < chop)
415  {
416  ostr << "0 "; // otherwise segfault on destruction
417  // if using only vstr.push_back("0 ");
418  // bug in MATLAB libmx
419  vstr.push_back(ostr.str());
420  } else if (std::abs(re) < chop)
421  {
422  ostr << im;
423  vstr.push_back(ostr.str() + "i");
424  } else if (std::abs(im) < chop)
425  {
426  ostr << re;
427  vstr.push_back(ostr.str() + " ");
428  } else
429  {
430  ostr << re;
431  strA = ostr.str();
432 
433  strA += (im > 0 ? " + " : " - ");
434  ostr.clear();
435  ostr.str(std::string()); // clear
436  ostr << std::abs(im);
437  strA += ostr.str();
438  strA += "i";
439  vstr.push_back(strA);
440  }
441  }
442  }
443 
444  // determine the maximum lenght of the entries in each column
445  std::vector<idx> maxlengthcols(A.cols(), 0);
446 
447  for (idx i = 0; i < static_cast<idx>(A.rows());
448  ++i)
449  for (idx j = 0;
450  j < static_cast<idx>(A.cols()); ++j)
451  if (vstr[i * A.cols() + j].size() > maxlengthcols[j])
452  maxlengthcols[j] = vstr[i * A.cols() + j].size();
453 
454  // finally display it!
455  for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
456  {
457  os << std::setw(static_cast<int>(maxlengthcols[0])) << std::right
458  << vstr[i * A.cols()]; // display first column
459  // then the rest
460  for (idx j = 1;
461  j < static_cast<idx>(A.cols()); ++j)
462  os << std::setw(static_cast<int>(maxlengthcols[j] + 2))
463  << std::right << vstr[i * A.cols() + j];
464 
465  if (i < static_cast<idx>(A.rows()) - 1)
466  os << std::endl;
467  }
468 
469  return os;
470  }
471 };
472 
473 } /* namespace internal */
474 } /* namespace qpp */
475 
476 #endif /* INTERNAL_UTIL_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:127
constexpr double chop
Used in qpp::disp() for setting to zero numbers that have their absolute value smaller than qpp::chop...
Definition: constants.h:59
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:74
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:198
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:272
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:120
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:78
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:156
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
std::ostream & display_impl_(const T &A, std::ostream &os, double chop=qpp::chop) const
Definition: util.h:391
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:48
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:257
idx get_num_subsys(idx sz, idx d)
Definition: util.h:370
bool check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:113
bool check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:188
bool check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:106
Definition: util.h:387
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:167
bool check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:178
bool check_matching_sizes(const T1 &lhs, const T2 &rhs) noexcept
Definition: util.h:134
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:46
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:356
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:315
bool check_qubit_rvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:242
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:140
bool check_qubit_vector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:249
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool check_qubit_cvector(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:235
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:228