Quantum++  v0.8
C++11 quantum computing library
functions.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 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 FUNCTIONS_H_
28 #define FUNCTIONS_H_
29 
30 // Collection of generic quantum computing functions
31 
32 namespace qpp
33 {
34 
35 // Eigen function wrappers
43 template<typename Derived>
45  const Eigen::MatrixBase<Derived>& A)
46 {
48 
49  // check zero-size
51  throw Exception("qpp::transpose()", Exception::Type::ZERO_SIZE);
52 
53  return rA.transpose();
54 }
55 
63 template<typename Derived>
65  const Eigen::MatrixBase<Derived>& A)
66 {
68 
69  // check zero-size
71  throw Exception("qpp::conjugate()", Exception::Type::ZERO_SIZE);
72 
73  return rA.conjugate();
74 }
75 
83 template<typename Derived>
84 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A)
85 {
87 
88  // check zero-size
90  throw Exception("qpp::adjoint()", Exception::Type::ZERO_SIZE);
91 
92  return rA.adjoint();
93 }
94 
102 template<typename Derived>
103 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A)
104 {
105  const dyn_mat<typename Derived::Scalar>& rA = A;
106 
107  // check zero-size
109  throw Exception("qpp::inverse()", Exception::Type::ZERO_SIZE);
110 
111  return rA.inverse();
112 }
113 
120 template<typename Derived>
121 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A)
122 {
123  const dyn_mat<typename Derived::Scalar>& rA = A;
124 
125  // check zero-size
127  throw Exception("qpp::trace()", Exception::Type::ZERO_SIZE);
128 
129  return rA.trace();
130 }
131 
139 template<typename Derived>
140 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A)
141 {
142  const dyn_mat<typename Derived::Scalar>& rA = A;
143 
144  // check zero-size
146  throw Exception("qpp::det()", Exception::Type::ZERO_SIZE);
147 
148  return rA.determinant();
149 }
150 
160 template<typename Derived>
161 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A)
162 {
163  const dyn_mat<typename Derived::Scalar>& rA = A;
164 
165  // check zero-size
167  throw Exception("qpp::logdet()", Exception::Type::ZERO_SIZE);
168 
169  // check square matrix
171  throw Exception("qpp::logdet()",
173 
174  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
176  lu.matrixLU().template triangularView<Eigen::Upper>();
177  typename Derived::Scalar result = std::log(U(0, 0));
178 
179  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
180  result += std::log(U(i, i));
181 
182  return result;
183 
184 }
185 
193 template<typename Derived>
194 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A)
195 {
196  const dyn_mat<typename Derived::Scalar>& rA = A;
197 
198  // check zero-size
200  throw Exception("qpp::sum()", Exception::Type::ZERO_SIZE);
201 
202  return rA.sum();
203 }
204 
212 template<typename Derived>
213 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A)
214 {
215  const dyn_mat<typename Derived::Scalar>& rA = A;
216 
217  // check zero-size
219  throw Exception("qpp::prod()", Exception::Type::ZERO_SIZE);
220 
221  return rA.prod();
222 }
223 
230 template<typename Derived>
231 double norm(const Eigen::MatrixBase<Derived>& A)
232 {
233  const dyn_mat<typename Derived::Scalar>& rA = A;
234 
235  // check zero-size
237  throw Exception("qpp::norm()", Exception::Type::ZERO_SIZE);
238 
239  // convert matrix to complex then return its norm
240  return (rA.template cast<cplx>()).norm();
241 }
242 
251 template<typename Derived>
252 std::pair<dyn_col_vect<cplx>, cmat> eig(const Eigen::MatrixBase<Derived>& A)
253 {
254  const dyn_mat<typename Derived::Scalar>& rA = A;
255 
256  // check zero-size
258  throw Exception("qpp::eig()", Exception::Type::ZERO_SIZE);
259 
260  // check square matrix
262  throw Exception("qpp::eig()", Exception::Type::MATRIX_NOT_SQUARE);
263 
264  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
265 
266  return std::make_pair(es.eigenvalues(), es.eigenvectors());
267 }
268 
276 template<typename Derived>
277 dyn_col_vect<cplx> evals(const Eigen::MatrixBase<Derived>& A)
278 {
279  const dyn_mat<typename Derived::Scalar>& rA = A;
280 
281  // check zero-size
283  throw Exception("qpp::evals()", Exception::Type::ZERO_SIZE);
284 
285  // check square matrix
287  throw Exception("qpp::evals()", Exception::Type::MATRIX_NOT_SQUARE);
288 
289  return eig(rA).first;
290 }
291 
299 template<typename Derived>
300 cmat evects(const Eigen::MatrixBase<Derived>& A)
301 {
302  const dyn_mat<typename Derived::Scalar>& rA = A;
303 
304  // check zero-size
306  throw Exception("qpp::evects()", Exception::Type::ZERO_SIZE);
307 
308  // check square matrix
310  throw Exception("qpp::evects()", Exception::Type::MATRIX_NOT_SQUARE);
311 
312  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
313 
314  return eig(rA).second;
315 }
316 
325 template<typename Derived>
326 std::pair<dyn_col_vect<double>, cmat> heig(
327  const Eigen::MatrixBase<Derived>& A)
328 {
329  const dyn_mat<typename Derived::Scalar>& rA = A;
330 
331  // check zero-size
333  throw Exception("qpp::heig()", Exception::Type::ZERO_SIZE);
334 
335  // check square matrix
337  throw Exception("qpp::heig()", Exception::Type::MATRIX_NOT_SQUARE);
338 
339  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
340 
341  return std::make_pair(es.eigenvalues(), es.eigenvectors());
342 }
343 
351 template<typename Derived>
352 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A)
353 {
354  const dyn_mat<typename Derived::Scalar>& rA = A;
355 
356  // check zero-size
358  throw Exception("qpp::hevals()", Exception::Type::ZERO_SIZE);
359 
360  // check square matrix
362  throw Exception("qpp::hevals()", Exception::Type::MATRIX_NOT_SQUARE);
363 
364  return heig(rA).first;
365 }
366 
374 template<typename Derived>
375 cmat hevects(const Eigen::MatrixBase<Derived>& A)
376 {
377  const dyn_mat<typename Derived::Scalar>& rA = A;
378 
379  // check zero-size
381  throw Exception("qpp::hevects()", Exception::Type::ZERO_SIZE);
382 
383  // check square matrix
385  throw Exception("qpp::hevects()",
387 
388  return heig(rA).second;
389 }
390 
400 template<typename Derived>
401 std::tuple<cmat, dyn_col_vect<double>, cmat>
402 svd(const Eigen::MatrixBase<Derived>& A)
403 {
404  const dyn_mat<typename Derived::Scalar>& rA = A;
405 
406  // check zero-size
408  throw Exception("qpp::svd()", Exception::Type::ZERO_SIZE);
409 
410  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
411  sv(rA,
412  Eigen::DecompositionOptions::ComputeFullU |
413  Eigen::DecompositionOptions::ComputeFullV);
414 
415  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
416 }
417 
425 template<typename Derived>
426 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A)
427 {
428  const dyn_mat<typename Derived::Scalar>& rA = A;
429 
430  // check zero-size
432  throw Exception("qpp::svals()", Exception::Type::ZERO_SIZE);
433 
434  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
435 
436  return sv.singularValues();
437 }
438 
446 template<typename Derived>
447 cmat svdU(const Eigen::MatrixBase<Derived>& A)
448 {
449  const dyn_mat<typename Derived::Scalar>& rA = A;
450 
451  // check zero-size
453  throw Exception("qpp::svdU()", Exception::Type::ZERO_SIZE);
454 
455  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
456  sv(rA, Eigen::DecompositionOptions::ComputeFullU);
457 
458  return sv.matrixU();
459 }
460 
468 template<typename Derived>
469 cmat svdV(const Eigen::MatrixBase<Derived>& A)
470 {
471  const dyn_mat<typename Derived::Scalar>& rA = A;
472 
473  // check zero-size
475  throw Exception("qpp::svdV()", Exception::Type::ZERO_SIZE);
476 
477  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
478  sv(rA, Eigen::DecompositionOptions::ComputeFullV);
479 
480  return sv.matrixV();
481 }
482 
483 // Matrix functional calculus
484 
492 template<typename Derived>
493 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (* f)(const cplx&))
494 {
495  const dyn_mat<typename Derived::Scalar>& rA = A;
496 
497  // check zero-size
499  throw Exception("qpp::funm()", Exception::Type::ZERO_SIZE);
500 
501  // check square matrix
503  throw Exception("qpp::funm()", Exception::Type::MATRIX_NOT_SQUARE);
504 
505  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
506  cmat evects = es.eigenvectors();
507  cmat evals = es.eigenvalues();
508  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
509  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
510 
511  cmat evalsdiag = evals.asDiagonal();
512 
513  return evects * evalsdiag * evects.inverse();
514 }
515 
522 template<typename Derived>
523 cmat sqrtm(const Eigen::MatrixBase<Derived>& A)
524 {
525  const dyn_mat<typename Derived::Scalar>& rA = A;
526 
527  // check zero-size
529  throw Exception("qpp::sqrtm()", Exception::Type::ZERO_SIZE);
530 
531  // check square matrix
533  throw Exception("qpp::sqrtm()", Exception::Type::MATRIX_NOT_SQUARE);
534 
535  return funm(rA, &std::sqrt);
536 }
537 
544 template<typename Derived>
545 cmat absm(const Eigen::MatrixBase<Derived>& A)
546 {
547  const dyn_mat<typename Derived::Scalar>& rA = A;
548 
549  // check zero-size
551  throw Exception("qpp::absm()", Exception::Type::ZERO_SIZE);
552 
553  // check square matrix
555  throw Exception("qpp::absm()", Exception::Type::MATRIX_NOT_SQUARE);
556 
557  return sqrtm(adjoint(rA) * rA);
558 }
559 
566 template<typename Derived>
567 cmat expm(const Eigen::MatrixBase<Derived>& A)
568 {
569  const dyn_mat<typename Derived::Scalar>& rA = A;
570 
571  // check zero-size
573  throw Exception("qpp::expm()", Exception::Type::ZERO_SIZE);
574 
575  // check square matrix
577  throw Exception("qpp::expm()", Exception::Type::MATRIX_NOT_SQUARE);
578 
579  return funm(rA, &std::exp);
580 }
581 
588 template<typename Derived>
589 cmat logm(const Eigen::MatrixBase<Derived>& A)
590 {
591  const dyn_mat<typename Derived::Scalar>& rA = A;
592 
593  // check zero-size
595  throw Exception("qpp::logm()", Exception::Type::ZERO_SIZE);
596 
597  // check square matrix
599  throw Exception("qpp::logm()", Exception::Type::MATRIX_NOT_SQUARE);
600 
601  return funm(rA, &std::log);
602 }
603 
610 template<typename Derived>
611 cmat sinm(const Eigen::MatrixBase<Derived>& A)
612 {
613  const dyn_mat<typename Derived::Scalar>& rA = A;
614 
615  // check zero-size
617  throw Exception("qpp::sinm()", Exception::Type::ZERO_SIZE);
618 
619  // check square matrix
621  throw Exception("qpp::sinm()", Exception::Type::MATRIX_NOT_SQUARE);
622 
623  return funm(rA, &std::sin);
624 }
625 
632 template<typename Derived>
633 cmat cosm(const Eigen::MatrixBase<Derived>& A)
634 {
635  const dyn_mat<typename Derived::Scalar>& rA = A;
636 
637  // check zero-size
639  throw Exception("qpp::cosm()", Exception::Type::ZERO_SIZE);
640 
641  // check square matrix
643  throw Exception("qpp::cosm()", Exception::Type::MATRIX_NOT_SQUARE);
644 
645  return funm(rA, &std::cos);
646 }
647 
659 template<typename Derived>
660 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z)
661 {
662  const dyn_mat<typename Derived::Scalar>& rA = A;
663 
664  // check zero-size
666  throw Exception("qpp::spectralpowm()", Exception::Type::ZERO_SIZE);
667 
668  // check square matrix
670  throw Exception("qpp::spectralpowm()",
672 
673  // Define A^0 = Id, for z IDENTICALLY zero
674  if (real(z) == 0 && imag(z) == 0)
675  return cmat::Identity(rA.rows(), rA.rows());
676 
677  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
678  cmat evects = es.eigenvectors();
679  cmat evals = es.eigenvalues();
680  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
681  evals(i) = std::pow(evals(i), z);
682 
683  cmat evalsdiag = evals.asDiagonal();
684 
685  return evects * evalsdiag * evects.inverse();
686 }
687 
700 template<typename Derived>
701 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
702  idx n)
703 {
705 
706  // check zero-size
708  throw Exception("qpp::powm()", Exception::Type::ZERO_SIZE);
709 
710  // check square matrix
712  throw Exception("qpp::powm()", Exception::Type::MATRIX_NOT_SQUARE);
713 
716  cA.rows(), cA.rows());
717 
718  // fast matrix power
719  for (; n > 0; n /= 2)
720  {
721  if (n % 2)
722  result = (result * cA).eval();
723  cA = (cA * cA).eval();
724  }
725 
726  return result;
727 }
728 
737 template<typename Derived>
738 double schatten(const Eigen::MatrixBase<Derived>& A, double p)
739 {
740  const dyn_mat<typename Derived::Scalar>& rA = A;
741 
742  // check zero-size
744  throw Exception("qpp::schatten()", Exception::Type::ZERO_SIZE);
745  if (p < 1)
746  throw Exception("qpp::schatten()", Exception::Type::OUT_OF_RANGE);
747 
748  if (p == infty) // infinity norm (largest singular value)
749  return svals(rA)(0);
750 
751  return std::pow(trace(powm(absm(rA), p)).real(), 1. / p);
752 }
753 
754 // other functions
755 
764 template<typename OutputScalar, typename Derived>
765 dyn_mat<OutputScalar> cwise(const Eigen::MatrixBase<Derived>& A,
766  OutputScalar (* f)(const typename Derived::Scalar&))
767 {
768  const dyn_mat<typename Derived::Scalar>& rA = A;
769 
770  // check zero-size
772  throw Exception("qpp::cwise()", Exception::Type::ZERO_SIZE);
773 
774  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
775 
776 #pragma omp parallel for collapse(2)
777  // column major order for speed
778  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
779  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
780  result(i, j) = (*f)(rA(i, j));
781 
782  return result;
783 }
784 
785 // Kronecker product of multiple matrices, preserve return type
786 // variadic template
797 template<typename T>
799 {
800  return head;
801 }
802 
813 template<typename T, typename ... Args>
814 dyn_mat<typename T::Scalar> kron(const T& head, const Args& ... tail)
815 {
816  return internal::_kron2(head, kron(tail...));
817 }
818 
828 template<typename Derived>
829 dyn_mat<typename Derived::Scalar> kron(const std::vector <Derived>& As)
830 {
831  if (As.size() == 0)
832  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
833 
834  for (auto&& it : As)
836  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
837 
838  dyn_mat<typename Derived::Scalar> result = As[0];
839  for (idx i = 1; i < As.size(); ++i)
840  {
841  result = kron(result, As[i]);
842  }
843 
844  return result;
845 }
846 
847 // Kronecker product of a list of matrices, preserve return type
848 // deduce the template parameters from initializer_list
859 template<typename Derived>
861  const std::initializer_list<Derived>& As)
862 {
863  return kron(std::vector<Derived>(As));
864 }
865 
875 template<typename Derived>
876 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
877  idx n)
878 {
879  const dyn_mat<typename Derived::Scalar>& rA = A;
880 
881  // check zero-size
883  throw Exception("qpp::kronpow()", Exception::Type::ZERO_SIZE);
884 
885  // check out of range
886  if (n == 0)
887  throw Exception("qpp::kronpow()", Exception::Type::OUT_OF_RANGE);
888 
889  std::vector <dyn_mat<typename Derived::Scalar>> As(n, rA);
890 
891  return kron(As);
892 }
893 
894 // Direct sum of multiple matrices, preserve return type
895 // variadic template
906 template<typename T>
908 {
909  return head;
910 }
911 
922 template<typename T, typename ... Args>
923 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args& ... tail)
924 {
925  return internal::_dirsum2(head, dirsum(tail...));
926 }
927 
937 template<typename Derived>
938 dyn_mat<typename Derived::Scalar> dirsum(const std::vector <Derived>& As)
939 {
940  if (As.size() == 0)
941  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
942 
943  for (auto&& it : As)
945  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
946 
947  idx total_rows = 0, total_cols = 0;
948  for (idx i = 0; i < As.size(); ++i)
949  {
950  total_rows += static_cast<idx>(As[i].rows());
951  total_cols += static_cast<idx>(As[i].cols());
952  }
954  dyn_mat < typename Derived::Scalar > ::Zero(total_rows, total_cols);
955 
956  idx cur_row = 0, cur_col = 0;
957  for (idx i = 0; i < As.size(); ++i)
958  {
959  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
960  cur_row += static_cast<idx>(As[i].rows());
961  cur_col += static_cast<idx>(As[i].cols());
962  }
963 
964  return result;
965 }
966 
967 // Direct sum of a list of matrices, preserve return type
968 // deduce the template parameters from initializer_list
979 template<typename Derived>
981  const std::initializer_list<Derived>& As)
982 {
983  return dirsum(std::vector<Derived>(As));
984 }
985 
995 template<typename Derived>
997  const Eigen::MatrixBase<Derived>& A, idx n)
998 {
999  const dyn_mat<typename Derived::Scalar>& rA = A;
1000 
1001  // check zero-size
1003  throw Exception("qpp::dirsumpow()", Exception::Type::ZERO_SIZE);
1004 
1005  // check out of range
1006  if (n == 0)
1007  throw Exception("qpp::dirsumpow()", Exception::Type::OUT_OF_RANGE);
1008 
1009  std::vector <dyn_mat<typename Derived::Scalar>> As(n, rA);
1010 
1011  return dirsum(As);
1012 }
1013 
1025 template<typename Derived>
1026 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1027  idx rows, idx cols)
1028 {
1029  const dyn_mat<typename Derived::Scalar>& rA = A;
1030 
1031  idx Arows = static_cast<idx>(rA.rows());
1032  idx Acols = static_cast<idx>(rA.cols());
1033 
1034  // check zero-size
1036  throw Exception("qpp::reshape()", Exception::Type::ZERO_SIZE);
1037 
1038  if (Arows * Acols != rows * cols)
1039  throw Exception("qpp::reshape()",
1041 
1042  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1043  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1044 }
1045 
1058 template<typename Derived1, typename Derived2>
1059 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1060  const Eigen::MatrixBase<Derived2>& B)
1061 {
1062  const dyn_mat<typename Derived1::Scalar>& rA = A;
1063  const dyn_mat<typename Derived2::Scalar>& rB = B;
1064 
1065  // EXCEPTION CHECKS
1066 
1067  // check types
1068  if (!std::is_same<typename Derived1::Scalar,
1069  typename Derived2::Scalar>::value)
1070  throw Exception("qpp::comm()", Exception::Type::TYPE_MISMATCH);
1071 
1072  // check zero-size
1075  throw Exception("qpp::comm()", Exception::Type::ZERO_SIZE);
1076 
1077  // check square matrices
1079  throw Exception("qpp::comm()", Exception::Type::MATRIX_NOT_SQUARE);
1080 
1081  // check equal dimensions
1082  if (rA.rows() != rB.rows())
1083  throw Exception("qpp::comm()", Exception::Type::DIMS_NOT_EQUAL);
1084 
1085  return rA * rB - rB * rA;
1086 }
1087 
1100 template<typename Derived1, typename Derived2>
1102  const Eigen::MatrixBase<Derived1>& A,
1103  const Eigen::MatrixBase<Derived2>& B)
1104 {
1105  const dyn_mat<typename Derived1::Scalar>& rA = A;
1106  const dyn_mat<typename Derived2::Scalar>& rB = B;
1107 
1108  // EXCEPTION CHECKS
1109 
1110  // check types
1111  if (!std::is_same<typename Derived1::Scalar,
1112  typename Derived2::Scalar>::value)
1113  throw Exception("qpp::anticomm()", Exception::Type::TYPE_MISMATCH);
1114 
1115  // check zero-size
1118  throw Exception("qpp::anticomm()", Exception::Type::ZERO_SIZE);
1119 
1120  // check square matrices
1122  throw Exception("qpp::anticomm()",
1124 
1125  // check equal dimensions
1126  if (rA.rows() != rB.rows())
1127  throw Exception("qpp::anticomm()", Exception::Type::DIMS_NOT_EQUAL);
1128 
1129  return rA * rB + rB * rA;
1130 }
1131 
1142 template<typename Derived>
1143 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& V)
1144 {
1145  const dyn_mat<typename Derived::Scalar>& rV = V;
1146 
1147  // check zero-size
1149  throw Exception("qpp::prj()", Exception::Type::ZERO_SIZE);
1150 
1151  // check column vector
1152  if (!internal::_check_cvector(rV))
1153  throw Exception("qpp::prj()", Exception::Type::MATRIX_NOT_CVECTOR);
1154 
1155  double normV = norm(rV);
1156  if (normV > eps)
1157  return rV * adjoint(rV) / (normV * normV);
1158  else
1160  ::Zero(rV.rows(), rV.rows());
1161 
1162 }
1163 
1171 template<typename Derived>
1172 dyn_mat<typename Derived::Scalar> grams(const std::vector <Derived>& Vs)
1173 {
1174  // check empty list
1176  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1177 
1178  for (auto&& it : Vs)
1180  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1181 
1182  // check that Vs[0] is a column vector
1183  if (!internal::_check_cvector(Vs[0]))
1184  throw Exception("qpp::grams()", Exception::Type::MATRIX_NOT_CVECTOR);
1185 
1186  // now check that all the rest match the size of the first vector
1187  for (auto&& it : Vs)
1188  if (it.rows() != Vs[0].rows() || it.cols() != 1)
1189  throw Exception("qpp::grams()", Exception::Type::DIMS_NOT_EQUAL);
1190 
1192  dyn_mat < typename Derived::Scalar > ::Identity(Vs[0].rows(),
1193  Vs[0].rows());
1194 
1196  dyn_mat < typename Derived::Scalar > ::Zero(Vs[0].rows(), 1);
1197 
1198  std::vector <dyn_mat<typename Derived::Scalar>> outvecs;
1199  // find the first non-zero vector in the list
1200  idx pos = 0;
1201  for (pos = 0; pos < Vs.size(); ++pos)
1202  {
1203  if (norm(Vs[pos]) > eps) // add it as the first element
1204  {
1205  outvecs.push_back(Vs[pos]);
1206  break;
1207  }
1208  }
1209 
1210  // start the process
1211  for (idx i = pos + 1; i < Vs.size(); ++i)
1212  {
1213  cut -= prj(outvecs[i - 1 - pos]);
1214  vi = cut * Vs[i];
1215  outvecs.push_back(vi);
1216  }
1217 
1218  dyn_mat<typename Derived::Scalar> result(Vs[0].rows(), outvecs.size());
1219 
1220  idx cnt = 0;
1221  for (auto&& it : outvecs)
1222  {
1223  double normV = norm(it);
1224  if (normV > eps) // we add only the non-zero vectors
1225  {
1226  result.col(cnt) = it / normV;
1227  cnt++;
1228  }
1229  }
1230 
1231  return result.block(0, 0, Vs[0].rows(), cnt);
1232 }
1233 
1234 // deduce the template parameters from initializer_list
1242 template<typename Derived>
1244  const std::initializer_list<Derived>& Vs)
1245 {
1246  return grams(std::vector<Derived>(Vs));
1247 }
1248 
1256 template<typename Derived>
1257 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A)
1258 {
1259  const dyn_mat<typename Derived::Scalar>& rA = A;
1260 
1262  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1263 
1264  std::vector <dyn_mat<typename Derived::Scalar>> input;
1265 
1266  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1267  input.push_back(rA.col(i));
1268 
1269  return grams<dyn_mat<typename Derived::Scalar>>(input);
1270 }
1271 
1282 inline std::vector <idx> n2multiidx(idx n, const std::vector <idx>& dims)
1283 {
1284  if (!internal::_check_dims(dims))
1285  throw Exception("qpp::n2multiidx()", Exception::Type::DIMS_INVALID);
1286 
1287  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1288  static_cast<idx>(1), std::multiplies<idx>()))
1289  throw Exception("qpp::n2multiidx()", Exception::Type::OUT_OF_RANGE);
1290 
1291  // double the size for matrices reshaped as vectors
1292  idx result[2 * maxn];
1293  internal::_n2multiidx(n, dims.size(), dims.data(), result);
1294 
1295  return std::vector<idx>(result, result + dims.size());
1296 }
1297 
1308 inline idx multiidx2n(const std::vector <idx>& midx,
1309  const std::vector <idx>& dims)
1310 {
1311  if (!internal::_check_dims(dims))
1312  throw Exception("qpp::multiidx2n()", Exception::Type::DIMS_INVALID);
1313 
1314  for (idx i = 0; i < dims.size(); ++i)
1315  if (midx[i] >= dims[i])
1316  throw Exception("qpp::multiidx2n()",
1318 
1319  return internal::_multiidx2n(midx.data(), dims.size(), dims.data());
1320 }
1321 
1334 inline ket mket(const std::vector <idx>& mask,
1335  const std::vector <idx>& dims)
1336 {
1337  idx n = mask.size();
1338 
1339  idx D = std::accumulate(std::begin(dims), std::end(dims),
1340  static_cast<idx>(1), std::multiplies<idx>());
1341 
1342  // check zero size
1343  if (n == 0)
1344  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1345  // check valid dims
1346  if (!internal::_check_dims(dims))
1347  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1348  // check mask and dims have the same size
1349  if (mask.size() != dims.size())
1350  throw Exception("qpp::mket()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1351  // check mask is a valid vector
1352  for (idx i = 0; i < n; ++i)
1353  if (mask[i] >= dims[i])
1354  throw Exception("qpp::mket()",
1356 
1357  ket result = ket::Zero(D);
1358  idx pos = multiidx2n(mask, dims);
1359  result(pos) = 1;
1360 
1361  return result;
1362 }
1363 
1376 inline ket mket(const std::vector <idx>& mask, idx d = 2)
1377 {
1378  idx n = mask.size();
1379  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1380 
1381  // check zero size
1382  if (n == 0)
1383  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1384  // check valid dims
1385  if (d == 0)
1386  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1387  // check mask is a valid vector
1388  for (idx i = 0; i < n; ++i)
1389  if (mask[i] >= d)
1390  throw Exception("qpp::mket()",
1392 
1393  ket result = ket::Zero(D);
1394  std::vector <idx> dims(n, d);
1395  idx pos = multiidx2n(mask, dims);
1396  result(pos) = 1;
1397 
1398  return result;
1399 }
1400 
1415 inline cmat mprj(const std::vector <idx>& mask,
1416  const std::vector <idx>& dims)
1417 {
1418  idx n = mask.size();
1419 
1420  idx D = std::accumulate(std::begin(dims), std::end(dims),
1421  static_cast<idx>(1), std::multiplies<idx>());
1422 
1423  // check zero size
1424  if (n == 0)
1425  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1426  // check valid dims
1427  if (!internal::_check_dims(dims))
1428  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1429  // check mask and dims have the same size
1430  if (mask.size() != dims.size())
1431  throw Exception("qpp::mprj()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1432  // check mask is a valid vector
1433  for (idx i = 0; i < n; ++i)
1434  if (mask[i] >= dims[i])
1435  throw Exception("qpp::mprj()",
1437 
1438  cmat result = cmat::Zero(D, D);
1439  idx pos = multiidx2n(mask, dims);
1440  result(pos, pos) = 1;
1441 
1442  return result;
1443 }
1444 
1459 inline cmat mprj(const std::vector <idx>& mask, idx d = 2)
1460 {
1461  idx n = mask.size();
1462  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1463 
1464  // check zero size
1465  if (n == 0)
1466  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1467  // check valid dims
1468  if (d == 0)
1469  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1470  // check mask is a valid vector
1471  for (idx i = 0; i < n; ++i)
1472  if (mask[i] >= d)
1473  throw Exception("qpp::mprj()",
1475 
1476  cmat result = cmat::Zero(D, D);
1477  std::vector <idx> dims(n, d);
1478  idx pos = multiidx2n(mask, dims);
1479  result(pos, pos) = 1;
1480 
1481  return result;
1482 }
1483 
1491 template<typename InputIterator>
1492 std::vector <double> abssq(InputIterator first, InputIterator last)
1493 {
1494  std::vector <double> weights(std::distance(first, last));
1495  std::transform(first, last, std::begin(weights),
1496  [](cplx z) -> double
1497  {
1498  return std::norm(z);
1499  });
1500 
1501  return weights;
1502 }
1503 
1510 template<typename Derived>
1511 std::vector <double> abssq(const Eigen::MatrixBase<Derived>& V)
1512 {
1513  const dyn_mat<typename Derived::Scalar>& rV = V;
1514 
1515  // check zero-size
1517  throw Exception("qpp::abssq()", Exception::Type::ZERO_SIZE);
1518 
1519  // check column vector
1520  if (!internal::_check_cvector(rV))
1521  throw Exception("qpp::abssq()", Exception::Type::MATRIX_NOT_CVECTOR);
1522 
1523  return abssq(rV.data(), rV.data() + rV.rows());
1524 }
1525 
1534 template<typename InputIterator>
1535 typename std::iterator_traits<InputIterator>::value_type
1536 sum(InputIterator first, InputIterator last)
1537 {
1538  using value_type =
1539  typename std::iterator_traits<InputIterator>::value_type;
1540 
1541  return std::accumulate(first, last, static_cast<value_type>(0));
1542 }
1543 
1551 template<typename Container>
1552 typename Container::value_type
1553 sum(const Container& c)
1554 {
1555  using value_type = typename Container::value_type;
1556 
1557  return std::accumulate(std::begin(c), std::end(c),
1558  static_cast<value_type>(0));
1559 }
1560 
1569 template<typename InputIterator>
1570 typename std::iterator_traits<InputIterator>::value_type
1571 prod(InputIterator first, InputIterator last)
1572 {
1573  using value_type =
1574  typename std::iterator_traits<InputIterator>::value_type;
1575 
1576  return std::accumulate(first, last, static_cast<value_type>(1),
1577  std::multiplies<value_type>());
1578 }
1579 
1580 
1588 template<typename Container>
1589 typename Container::value_type
1590 prod(const Container& c)
1591 {
1592  using value_type = typename Container::value_type;
1593 
1594  return std::accumulate(std::begin(c), std::end(c),
1595  static_cast<value_type>(1),
1596  std::multiplies<value_type>());
1597 }
1598 
1611 template<typename Derived>
1613  const Eigen::MatrixBase<Derived>& A)
1614 {
1615  const dyn_mat<typename Derived::Scalar>& rA = A;
1616 
1617  // EXCEPTION CHECKS
1618  // check zero-size
1620  throw Exception("qpp::rho2pure()", Exception::Type::ZERO_SIZE);
1621  // check square matrix
1622  if (!internal::_check_square_mat(rA))
1623  throw Exception("qpp::rho2pure()", Exception::Type::MATRIX_NOT_SQUARE);
1624  // END EXPCEPTION CHECKS
1625 
1626  dyn_col_vect<double> tmp_evals = hevals(rA);
1627  cmat tmp_evects = hevects(rA);
1630  // find the non-zero eigenvector
1631  // there is only one, assuming the state is pure
1632  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1633  {
1634  if (std::abs(tmp_evals(k)) > eps)
1635  {
1636  result = tmp_evects.col(k);
1637  break;
1638  }
1639  }
1640 
1641  return result;
1642 }
1643 
1652 template<typename T>
1653 std::vector <T> complement(std::vector <T> subsys, idx N)
1654 {
1655  if (N < subsys.size())
1656  throw Exception("qpp::complement()", Exception::Type::OUT_OF_RANGE);
1657 
1658  std::vector <T> all(N);
1659  std::vector <T> subsys_bar(N - subsys.size());
1660 
1661  std::iota(std::begin(all), std::end(all), 0);
1662  std::sort(std::begin(subsys), std::end(subsys));
1663  std::set_difference(std::begin(all), std::end(all),
1664  std::begin(subsys), std::end(subsys),
1665  std::begin(subsys_bar));
1666 
1667  return subsys_bar;
1668 }
1669 
1680 template<typename Derived>
1681 std::vector <double> rho2bloch(
1682  const Eigen::MatrixBase<Derived>& A)
1683 {
1684  const dyn_mat<typename Derived::Scalar>& rA = A;
1685 
1686  // check qubit matrix
1688  throw Exception("qpp::rho2bloch()", Exception::Type::NOT_QUBIT_MATRIX);
1689 
1690  std::vector <double> result(3);
1691  cmat X(2, 2), Y(2, 2), Z(2, 2);
1692  X << 0, 1, 1, 0;
1693  Y << 0, -1_i, 1_i, 0;
1694  Z << 1, 0, 0, -1;
1695  result[0] = std::real(trace(rA * X));
1696  result[1] = std::real(trace(rA * Y));
1697  result[2] = std::real(trace(rA * Z));
1698 
1699  return result;
1700 }
1701 
1710 inline cmat bloch2rho(const std::vector <double>& r)
1711 {
1712  // check 3-dimensional vector
1713  if (r.size() != 3)
1714  throw Exception("qpp::bloch2rho", "r is not a 3-dimensional vector!");
1715 
1716  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1717  X << 0, 1, 1, 0;
1718  Y << 0, -1_i, 1_i, 0;
1719  Z << 1, 0, 0, -1;
1720  Id2 << 1, 0, 0, 1;
1721 
1722  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1723 }
1724 
1725 } /* namespace qpp */
1726 
1727 #endif /* FUNCTIONS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:326
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:161
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:82
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolut value.
Definition: functions.h:545
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:265
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:567
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:213
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:75
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:907
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:84
std::vector< double > rho2bloch(const Eigen::MatrixBase< Derived > &A)
Computes the 3-dimensional real Bloch vector corresponding to the qubit density matrix A...
Definition: functions.h:1681
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:57
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1101
Quantum++ main namespace.
Definition: codes.h:30
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:252
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:231
cmat spectralpowm(const Eigen::MatrixBase< Derived > &A, const cplx z)
Matrix power.
Definition: functions.h:660
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:352
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:765
dyn_col_vect< typename Derived::Scalar > rho2pure(const Eigen::MatrixBase< Derived > &A)
Finds the pure state representation of a matrix proportional to a projector onto a pure state...
Definition: functions.h:1612
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:996
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1653
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:300
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:493
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:738
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:447
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1710
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1282
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1308
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:426
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:798
bool _check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:213
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:52
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:54
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1415
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:140
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:44
dyn_mat< typename Derived::Scalar > kronpow(const Eigen::MatrixBase< Derived > &A, idx n)
Kronecker power.
Definition: functions.h:876
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Matrix power.
Definition: functions.h:701
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &V)
Projector.
Definition: functions.h:1143
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:277
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:375
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:64
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:96
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:305
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:121
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:469
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:633
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:96
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:611
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of a range of complex numbers.
Definition: functions.h:1492
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1059
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:67
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1334
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:41
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:103
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:589
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:523
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:119
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:194
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1026
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:402
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &Vs)
Gram-Schmidt orthogonalization.
Definition: functions.h:1172