Quantum++  v1.0-rc3
A modern C++11 quantum computing library
functions.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 FUNCTIONS_H_
33 #define FUNCTIONS_H_
34 
35 namespace qpp {
36 // Eigen function wrappers
44 template <typename Derived>
45 dyn_mat<typename Derived::Scalar>
46 transpose(const Eigen::MatrixBase<Derived>& A) {
47  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
48 
49  // check zero-size
51  throw exception::ZeroSize("qpp::transpose()");
52 
53  return rA.transpose();
54 }
55 
63 template <typename Derived>
64 dyn_mat<typename Derived::Scalar>
65 conjugate(const Eigen::MatrixBase<Derived>& A) {
66  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
67 
68  // check zero-size
70  throw exception::ZeroSize("qpp::conjugate()");
71 
72  return rA.conjugate();
73 }
74 
82 template <typename Derived>
83 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A) {
84  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
85 
86  // EXCEPTION CHECKS
87 
88  // check zero-size
90  throw exception::ZeroSize("qpp::adjoint()");
91  // END EXCEPTION CHECKS
92 
93  return rA.adjoint();
94 }
95 
103 template <typename Derived>
104 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A) {
105  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
106 
107  // EXCEPTION CHECKS
108 
109  // check zero-size
111  throw exception::ZeroSize("qpp::inverse()");
112  // END EXCEPTION CHECKS
113 
114  return rA.inverse();
115 }
116 
123 template <typename Derived>
124 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A) {
125  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
126 
127  // EXCEPTION CHECKS
128 
129  // check zero-size
131  throw exception::ZeroSize("qpp::trace()");
132  // END EXCEPTION CHECKS
133 
134  return rA.trace();
135 }
136 
144 template <typename Derived>
145 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A) {
146  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
147 
148  // EXCEPTION CHECKS
149 
150  // check zero-size
152  throw exception::ZeroSize("qpp::det()");
153  // END EXCEPTION CHECKS
154 
155  return rA.determinant();
156 }
157 
167 template <typename Derived>
168 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A) {
169  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
170 
171  // EXCEPTION CHECKS
172 
173  // check zero-size
175  throw exception::ZeroSize("qpp::logdet()");
176 
177  // check square matrix
179  throw exception::MatrixNotSquare("qpp::logdet()");
180  // END EXCEPTION CHECKS
181 
182  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
183  dyn_mat<typename Derived::Scalar> U =
184  lu.matrixLU().template triangularView<Eigen::Upper>();
185  typename Derived::Scalar result = std::log(U(0, 0));
186 
187  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
188  result += std::log(U(i, i));
189 
190  return result;
191 }
192 
200 template <typename Derived>
201 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A) {
202  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
203 
204  // EXCEPTION CHECKS
205 
206  // check zero-size
208  throw exception::ZeroSize("qpp::sum()");
209  // END EXCEPTION CHECKS
210 
211  return rA.sum();
212 }
213 
221 template <typename Derived>
222 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A) {
223  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
224 
225  // EXCEPTION CHECKS
226 
227  // check zero-size
229  throw exception::ZeroSize("qpp::prod()");
230  // END EXCEPTION CHECKS
231 
232  return rA.prod();
233 }
234 
241 template <typename Derived>
242 double norm(const Eigen::MatrixBase<Derived>& A) {
243  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
244 
245  // EXCEPTION CHECKS
246 
247  // check zero-size
249  throw exception::ZeroSize("qpp::norm()");
250  // END EXCEPTION CHECKS
251 
252  // convert matrix to complex then return its norm
253  return (rA.template cast<cplx>()).norm();
254 }
255 
264 template <typename Derived>
265 std::pair<dyn_col_vect<cplx>, cmat>
266 
267 eig(const Eigen::MatrixBase<Derived>& A) {
268  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
269 
270  // EXCEPTION CHECKS
271 
272  // check zero-size
274  throw exception::ZeroSize("qpp::eig()");
275 
276  // check square matrix
278  throw exception::MatrixNotSquare("qpp::eig()");
279  // END EXCEPTION CHECKS
280 
281  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
282 
283  return std::make_pair(es.eigenvalues(), es.eigenvectors());
284 }
285 
293 template <typename Derived>
294 dyn_col_vect<cplx> evals(const Eigen::MatrixBase<Derived>& A) {
295  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
296 
297  // EXCEPTION CHECKS
298 
299  // check zero-size
301  throw exception::ZeroSize("qpp::evals()");
302 
303  // check square matrix
305  throw exception::MatrixNotSquare("qpp::evals()");
306  // END EXCEPTION CHECKS
307 
308  return eig(rA).first;
309 }
310 
318 template <typename Derived>
319 cmat evects(const Eigen::MatrixBase<Derived>& A) {
320  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
321 
322  // EXCEPTION CHECKS
323 
324  // check zero-size
326  throw exception::ZeroSize("qpp::evects()");
327 
328  // check square matrix
330  throw exception::MatrixNotSquare("qpp::evects()");
331  // END EXCEPTION CHECKS
332 
333  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
334 
335  return eig(rA).second;
336 }
337 
346 template <typename Derived>
347 std::pair<dyn_col_vect<double>, cmat>
348 
349 heig(const Eigen::MatrixBase<Derived>& A) {
350  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
351 
352  // EXCEPTION CHECKS
353 
354  // check zero-size
356  throw exception::ZeroSize("qpp::heig()");
357 
358  // check square matrix
360  throw exception::MatrixNotSquare("qpp::heig()");
361  // END EXCEPTION CHECKS
362 
363  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
364 
365  return std::make_pair(es.eigenvalues(), es.eigenvectors());
366 }
367 
375 template <typename Derived>
376 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A) {
377  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
378 
379  // EXCEPTION CHECKS
380 
381  // check zero-size
383  throw exception::ZeroSize("qpp::hevals()");
384 
385  // check square matrix
387  throw exception::MatrixNotSquare("qpp::hevals()");
388  // END EXCEPTION CHECKS
389 
390  return heig(rA).first;
391 }
392 
400 template <typename Derived>
401 cmat hevects(const Eigen::MatrixBase<Derived>& A) {
402  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
403 
404  // EXCEPTION CHECKS
405 
406  // check zero-size
408  throw exception::ZeroSize("qpp::hevects()");
409 
410  // check square matrix
412  throw exception::MatrixNotSquare("qpp::hevects()");
413  // END EXCEPTION CHECKS
414 
415  return heig(rA).second;
416 }
417 
427 template <typename Derived>
428 std::tuple<cmat, dyn_col_vect<double>, cmat>
429 
430 svd(const Eigen::MatrixBase<Derived>& A) {
431  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
432 
433  // EXCEPTION CHECKS
434 
435  // check zero-size
437  throw exception::ZeroSize("qpp::svd()");
438  // END EXCEPTION CHECKS
439 
440  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
441  rA, Eigen::DecompositionOptions::ComputeFullU |
442  Eigen::DecompositionOptions::ComputeFullV);
443 
444  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
445 }
446 
454 template <typename Derived>
455 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A) {
456  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
457 
458  // EXCEPTION CHECKS
459 
460  // check zero-size
462  throw exception::ZeroSize("qpp::svals()");
463  // END EXCEPTION CHECKS
464 
465  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
466 
467  return sv.singularValues();
468 }
469 
477 template <typename Derived>
478 cmat svdU(const Eigen::MatrixBase<Derived>& A) {
479  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
480 
481  // EXCEPTION CHECKS
482 
483  // check zero-size
485  throw exception::ZeroSize("qpp::svdU()");
486  // END EXCEPTION CHECKS
487 
488  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
489  rA, Eigen::DecompositionOptions::ComputeFullU);
490 
491  return sv.matrixU();
492 }
493 
501 template <typename Derived>
502 cmat svdV(const Eigen::MatrixBase<Derived>& A) {
503  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
504 
505  // EXCEPTION CHECKS
506 
507  // check zero-size
509  throw exception::ZeroSize("qpp::svdV()");
510  // END EXCEPTION CHECKS
511 
512  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
513  rA, Eigen::DecompositionOptions::ComputeFullV);
514 
515  return sv.matrixV();
516 }
517 
518 // Matrix functional calculus
519 
527 template <typename Derived>
528 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (*f)(const cplx&)) {
529  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
530 
531  // EXCEPTION CHECKS
532 
533  // check zero-size
535  throw exception::ZeroSize("qpp::funm()");
536 
537  // check square matrix
539  throw exception::MatrixNotSquare("qpp::funm()");
540  // END EXCEPTION CHECKS
541 
542  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
543  cmat evects = es.eigenvectors();
544  cmat evals = es.eigenvalues();
545  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
546  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
547 
548  cmat evalsdiag = evals.asDiagonal();
549 
550  return evects * evalsdiag * evects.inverse();
551 }
552 
559 template <typename Derived>
560 cmat sqrtm(const Eigen::MatrixBase<Derived>& A) {
561  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
562 
563  // EXCEPTION CHECKS
564 
565  // check zero-size
567  throw exception::ZeroSize("qpp::sqrtm()");
568 
569  // check square matrix
571  throw exception::MatrixNotSquare("qpp::sqrtm()");
572  // END EXCEPTION CHECKS
573 
574  return funm(rA, &std::sqrt);
575 }
576 
583 template <typename Derived>
584 cmat absm(const Eigen::MatrixBase<Derived>& A) {
585  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
586 
587  // EXCEPTION CHECKS
588 
589  // check zero-size
591  throw exception::ZeroSize("qpp::absm()");
592 
593  // check square matrix
595  throw exception::MatrixNotSquare("qpp::absm()");
596  // END EXCEPTION CHECKS
597 
598  return sqrtm(adjoint(rA) * rA);
599 }
600 
607 template <typename Derived>
608 cmat expm(const Eigen::MatrixBase<Derived>& A) {
609  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
610 
611  // EXCEPTION CHECKS
612 
613  // check zero-size
615  throw exception::ZeroSize("qpp::expm()");
616 
617  // check square matrix
619  throw exception::MatrixNotSquare("qpp::expm()");
620  // END EXCEPTION CHECKS
621 
622  return funm(rA, &std::exp);
623 }
624 
631 template <typename Derived>
632 cmat logm(const Eigen::MatrixBase<Derived>& A) {
633  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
634 
635  // EXCEPTION CHECKS
636 
637  // check zero-size
639  throw exception::ZeroSize("qpp::logm()");
640 
641  // check square matrix
643  throw exception::MatrixNotSquare("qpp::logm()");
644  // END EXCEPTION CHECKS
645 
646  return funm(rA, &std::log);
647 }
648 
655 template <typename Derived>
656 cmat sinm(const Eigen::MatrixBase<Derived>& A) {
657  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
658 
659  // EXCEPTION CHECKS
660 
661  // check zero-size
663  throw exception::ZeroSize("qpp::sinm()");
664 
665  // check square matrix
667  throw exception::MatrixNotSquare("qpp::sinm()");
668  // END EXCEPTION CHECKS
669 
670  return funm(rA, &std::sin);
671 }
672 
679 template <typename Derived>
680 cmat cosm(const Eigen::MatrixBase<Derived>& A) {
681  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
682 
683  // EXCEPTION CHECKS
684 
685  // check zero-size
687  throw exception::ZeroSize("qpp::cosm()");
688 
689  // check square matrix
691  throw exception::MatrixNotSquare("qpp::cosm()");
692  // END EXCEPTION CHECKS
693 
694  return funm(rA, &std::cos);
695 }
696 
708 template <typename Derived>
709 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z) {
710  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
711 
712  // EXCEPTION CHECKS
713 
714  // check zero-size
716  throw exception::ZeroSize("qpp::spectralpowm()");
717 
718  // check square matrix
720  throw exception::MatrixNotSquare("qpp::spectralpowm()");
721  // END EXCEPTION CHECKS
722 
723  // Define A^0 = Id, for z IDENTICALLY zero
724  if (real(z) == 0 && imag(z) == 0)
725  return cmat::Identity(rA.rows(), rA.rows());
726 
727  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
728  cmat evects = es.eigenvectors();
729  cmat evals = es.eigenvalues();
730  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
731  evals(i) = std::pow(evals(i), z);
732 
733  cmat evalsdiag = evals.asDiagonal();
734 
735  return evects * evalsdiag * evects.inverse();
736 }
737 
750 template <typename Derived>
751 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
752  idx n) {
753  // EXCEPTION CHECKS
754 
755  // check zero-size
757  throw exception::ZeroSize("qpp::powm()");
758 
759  // check square matrix
761  throw exception::MatrixNotSquare("qpp::powm()");
762  // END EXCEPTION CHECKS
763 
764  // if n = 1, return the matrix unchanged
765  if (n == 1)
766  return A;
767 
768  dyn_mat<typename Derived::Scalar> result =
769  dyn_mat<typename Derived::Scalar>::Identity(A.rows(), A.rows());
770 
771  // if n = 0, return the identity (as just prepared in result)
772  if (n == 0)
773  return result;
774 
775  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
776 
777  // fast matrix power
778  for (; n > 0; n /= 2) {
779  if (n % 2)
780  result = (result * cA).eval();
781  cA = (cA * cA).eval();
782  }
783 
784  return result;
785 }
786 
795 template <typename Derived>
796 double schatten(const Eigen::MatrixBase<Derived>& A, double p) {
797  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
798 
799  // EXCEPTION CHECKS
800 
801  // check zero-size
803  throw exception::ZeroSize("qpp::schatten()");
804  if (p < 1)
805  throw exception::OutOfRange("qpp::schatten()");
806  // END EXCEPTION CHECKS
807 
808  if (p == infty) // infinity norm (largest singular value)
809  return svals(rA)(0);
810 
811  const dyn_col_vect<double> sv = svals(rA);
812  double result = 0;
813  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
814  result += std::pow(sv[i], p);
815 
816  return std::pow(result, 1. / p);
817 }
818 
819 // other functions
820 
829 template <typename OutputScalar, typename Derived>
830 dyn_mat<OutputScalar>
831 cwise(const Eigen::MatrixBase<Derived>& A,
832  OutputScalar (*f)(const typename Derived::Scalar&)) {
833  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
834 
835  // EXCEPTION CHECKS
836 
837  // check zero-size
839  throw exception::ZeroSize("qpp::cwise()");
840  // END EXCEPTION CHECKS
841 
842  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
843 
844 #ifdef WITH_OPENMP_
845 #pragma omp parallel for collapse(2)
846 #endif // WITH_OPENMP_
847  // column major order for speed
848  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
849  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
850  result(i, j) = (*f)(rA(i, j));
851 
852  return result;
853 }
854 
855 // Kronecker product of multiple matrices, preserve return type
856 // variadic template
867 template <typename T>
868 dyn_mat<typename T::Scalar> kron(const T& head) {
869  return head;
870 }
871 
882 template <typename T, typename... Args>
883 dyn_mat<typename T::Scalar> kron(const T& head, const Args&... tail) {
884  return internal::kron2(head, kron(tail...));
885 }
886 
896 template <typename Derived>
897 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As) {
898  // EXCEPTION CHECKS
899 
900  if (As.size() == 0)
901  throw exception::ZeroSize("qpp::kron()");
902 
903  for (auto&& it : As)
905  throw exception::ZeroSize("qpp::kron()");
906  // END EXCEPTION CHECKS
907 
908  dyn_mat<typename Derived::Scalar> result = As[0].derived();
909  for (idx i = 1; i < As.size(); ++i) {
910  result = kron(result, As[i]);
911  }
912 
913  return result;
914 }
915 
916 // Kronecker product of a list of matrices, preserve return type
917 // deduce the template parameters from initializer_list
928 template <typename Derived>
929 dyn_mat<typename Derived::Scalar>
930 kron(const std::initializer_list<Derived>& As) {
931  return kron(std::vector<Derived>(As));
932 }
933 
943 template <typename Derived>
944 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
945  idx n) {
946  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
947 
948  // EXCEPTION CHECKS
949 
950  // check zero-size
952  throw exception::ZeroSize("qpp::kronpow()");
953 
954  // check out of range
955  if (n == 0)
956  throw exception::OutOfRange("qpp::kronpow()");
957  // END EXCEPTION CHECKS
958 
959  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
960 
961  return kron(As);
962 }
963 
964 // Direct sum of multiple matrices, preserve return type
965 // variadic template
976 template <typename T>
977 dyn_mat<typename T::Scalar> dirsum(const T& head) {
978  return head;
979 }
980 
991 template <typename T, typename... Args>
992 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args&... tail) {
993  return internal::dirsum2(head, dirsum(tail...));
994 }
995 
1005 template <typename Derived>
1006 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As) {
1007  // EXCEPTION CHECKS
1008 
1009  if (As.size() == 0)
1010  throw exception::ZeroSize("qpp::dirsum()");
1011 
1012  for (auto&& it : As)
1014  throw exception::ZeroSize("qpp::dirsum()");
1015  // END EXCEPTION CHECKS
1016 
1017  idx total_rows = 0, total_cols = 0;
1018  for (idx i = 0; i < As.size(); ++i) {
1019  total_rows += static_cast<idx>(As[i].rows());
1020  total_cols += static_cast<idx>(As[i].cols());
1021  }
1022  dyn_mat<typename Derived::Scalar> result =
1023  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1024 
1025  idx cur_row = 0, cur_col = 0;
1026  for (idx i = 0; i < As.size(); ++i) {
1027  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1028  cur_row += static_cast<idx>(As[i].rows());
1029  cur_col += static_cast<idx>(As[i].cols());
1030  }
1031 
1032  return result;
1033 }
1034 
1035 // Direct sum of a list of matrices, preserve return type
1036 // deduce the template parameters from initializer_list
1047 template <typename Derived>
1048 dyn_mat<typename Derived::Scalar>
1049 dirsum(const std::initializer_list<Derived>& As) {
1050  return dirsum(std::vector<Derived>(As));
1051 }
1052 
1062 template <typename Derived>
1063 dyn_mat<typename Derived::Scalar> dirsumpow(const Eigen::MatrixBase<Derived>& A,
1064  idx n) {
1065  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1066 
1067  // EXCEPTION CHECKS
1068 
1069  // check zero-size
1071  throw exception::ZeroSize("qpp::dirsumpow()");
1072 
1073  // check out of range
1074  if (n == 0)
1075  throw exception::OutOfRange("qpp::dirsumpow()");
1076  // END EXCEPTION CHECKS
1077 
1078  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1079 
1080  return dirsum(As);
1081 }
1082 
1094 template <typename Derived>
1095 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1096  idx rows, idx cols) {
1097  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1098 
1099  idx Arows = static_cast<idx>(rA.rows());
1100  idx Acols = static_cast<idx>(rA.cols());
1101 
1102  // EXCEPTION CHECKS
1103 
1104  // check zero-size
1106  throw exception::ZeroSize("qpp::reshape()");
1107 
1108  if (Arows * Acols != rows * cols)
1109  throw exception::DimsMismatchMatrix("qpp::reshape()");
1110  // END EXCEPTION CHECKS
1111 
1112  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1113  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1114 }
1115 
1128 template <typename Derived1, typename Derived2>
1129 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1130  const Eigen::MatrixBase<Derived2>& B) {
1131  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1132  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1133 
1134  // EXCEPTION CHECKS
1135 
1136  // check types
1137  if (!std::is_same<typename Derived1::Scalar,
1138  typename Derived2::Scalar>::value)
1139  throw exception::TypeMismatch("qpp::comm()");
1140 
1141  // check zero-size
1143  throw exception::ZeroSize("qpp::comm()");
1144 
1145  // check square matrices
1147  throw exception::MatrixNotSquare("qpp::comm()");
1148 
1149  // check equal dimensions
1150  if (rA.rows() != rB.rows())
1151  throw exception::DimsNotEqual("qpp::comm()");
1152  // END EXCEPTION CHECKS
1153 
1154  return rA * rB - rB * rA;
1155 }
1156 
1169 template <typename Derived1, typename Derived2>
1170 dyn_mat<typename Derived1::Scalar>
1171 anticomm(const Eigen::MatrixBase<Derived1>& A,
1172  const Eigen::MatrixBase<Derived2>& B) {
1173  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1174  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1175 
1176  // EXCEPTION CHECKS
1177 
1178  // check types
1179  if (!std::is_same<typename Derived1::Scalar,
1180  typename Derived2::Scalar>::value)
1181  throw exception::TypeMismatch("qpp::anticomm()");
1182 
1183  // check zero-size
1185  throw exception::ZeroSize("qpp::anticomm()");
1186 
1187  // check square matrices
1189  throw exception::MatrixNotSquare("qpp::anticomm()");
1190 
1191  // check equal dimensions
1192  if (rA.rows() != rB.rows())
1193  throw exception::DimsNotEqual("qpp::anticomm()");
1194  // END EXCEPTION CHECKS
1195 
1196  return rA * rB + rB * rA;
1197 }
1198 
1209 template <typename Derived>
1210 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A) {
1211  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1212 
1213  // EXCEPTION CHECKS
1214 
1215  // check zero-size
1217  throw exception::ZeroSize("qpp::prj()");
1218 
1219  // check column vector
1220  if (!internal::check_cvector(rA))
1221  throw exception::MatrixNotCvector("qpp::prj()");
1222  // END EXCEPTION CHECKS
1223 
1224  double normA = norm(rA);
1225  if (normA > eps)
1226  return rA * adjoint(rA) / (normA * normA);
1227  else
1228  return dyn_mat<typename Derived::Scalar>::Zero(rA.rows(), rA.rows());
1229 }
1230 
1238 template <typename Derived>
1239 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As) {
1240  // EXCEPTION CHECKS
1241 
1242  // check empty list
1244  throw exception::ZeroSize("qpp::grams()");
1245 
1246  for (auto&& it : As)
1248  throw exception::ZeroSize("qpp::grams()");
1249 
1250  // check that As[0] is a column vector
1251  if (!internal::check_cvector(As[0]))
1252  throw exception::MatrixNotCvector("qpp::grams()");
1253 
1254  // now check that all the rest match the size of the first vector
1255  for (auto&& it : As)
1256  if (it.rows() != As[0].rows() || it.cols() != 1)
1257  throw exception::DimsNotEqual("qpp::grams()");
1258  // END EXCEPTION CHECKS
1259 
1260  dyn_mat<typename Derived::Scalar> cut =
1261  dyn_mat<typename Derived::Scalar>::Identity(As[0].rows(), As[0].rows());
1262 
1263  dyn_mat<typename Derived::Scalar> vi =
1264  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1265 
1266  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1267  // find the first non-zero vector in the list
1268  idx pos = 0;
1269  for (pos = 0; pos < As.size(); ++pos) {
1270  if (norm(As[pos]) > eps) // add it as the first element
1271  {
1272  outvecs.push_back(As[pos]);
1273  break;
1274  }
1275  }
1276 
1277  // start the process
1278  for (idx i = pos + 1; i < As.size(); ++i) {
1279  cut -= prj(outvecs[i - 1 - pos]);
1280  vi = cut * As[i];
1281  outvecs.push_back(vi);
1282  }
1283 
1284  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1285 
1286  idx cnt = 0;
1287  for (auto&& it : outvecs) {
1288  double normA = norm(it);
1289  if (normA > eps) // we add only the non-zero vectors
1290  {
1291  result.col(cnt) = it / normA;
1292  cnt++;
1293  }
1294  }
1295 
1296  return result.block(0, 0, As[0].rows(), cnt);
1297 }
1298 
1299 // deduce the template parameters from initializer_list
1307 template <typename Derived>
1308 dyn_mat<typename Derived::Scalar>
1309 grams(const std::initializer_list<Derived>& As) {
1310  return grams(std::vector<Derived>(As));
1311 }
1312 
1320 template <typename Derived>
1321 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A) {
1322  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1323 
1324  // EXCEPTION CHECKS
1325 
1327  throw exception::ZeroSize("qpp::grams()");
1328  // END EXCEPTION CHECKS
1329 
1330  std::vector<dyn_mat<typename Derived::Scalar>> input;
1331 
1332  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1333  input.push_back(rA.col(i));
1334 
1335  return grams<dyn_mat<typename Derived::Scalar>>(input);
1336 }
1337 
1348 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims) {
1349  // EXCEPTION CHECKS
1350 
1351  if (!internal::check_dims(dims))
1352  throw exception::DimsInvalid("qpp::n2multiidx()");
1353 
1354  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1355  static_cast<idx>(1), std::multiplies<idx>()))
1356  throw exception::OutOfRange("qpp::n2multiidx()");
1357  // END EXCEPTION CHECKS
1358 
1359  // double the size for matrices reshaped as vectors
1360  idx result[2 * maxn];
1361  internal::n2multiidx(n, dims.size(), dims.data(), result);
1362 
1363  return std::vector<idx>(result, result + dims.size());
1364 }
1365 
1376 inline idx multiidx2n(const std::vector<idx>& midx,
1377  const std::vector<idx>& dims) {
1378  // EXCEPTION CHECKS
1379 
1380  if (!internal::check_dims(dims))
1381  throw exception::DimsInvalid("qpp::multiidx2n()");
1382 
1383  for (idx i = 0; i < dims.size(); ++i)
1384  if (midx[i] >= dims[i])
1385  throw exception::OutOfRange("qpp::multiidx2n()");
1386  // END EXCEPTION CHECKS
1387 
1388  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1389 }
1390 
1403 inline ket mket(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1404  idx N = mask.size();
1405 
1406  idx D = std::accumulate(std::begin(dims), std::end(dims),
1407  static_cast<idx>(1), std::multiplies<idx>());
1408 
1409  // EXCEPTION CHECKS
1410 
1411  // check zero size
1412  if (N == 0)
1413  throw exception::ZeroSize("qpp::mket()");
1414  // check valid dims
1415  if (!internal::check_dims(dims))
1416  throw exception::DimsInvalid("qpp::mket()");
1417  // check mask and dims have the same size
1418  if (mask.size() != dims.size())
1419  throw exception::SubsysMismatchDims("qpp::mket()");
1420  // check mask is a valid vector
1421  for (idx i = 0; i < N; ++i)
1422  if (mask[i] >= dims[i])
1423  throw exception::SubsysMismatchDims("qpp::mket()");
1424  // END EXCEPTION CHECKS
1425 
1426  ket result = ket::Zero(D);
1427  idx pos = multiidx2n(mask, dims);
1428  result(pos) = 1;
1429 
1430  return result;
1431 }
1432 
1445 inline ket mket(const std::vector<idx>& mask, idx d = 2) {
1446  idx N = mask.size();
1447  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1448 
1449  // EXCEPTION CHECKS
1450 
1451  // check zero size
1452  if (N == 0)
1453  throw exception::ZeroSize("qpp::mket()");
1454 
1455  // check valid dims
1456  if (d == 0)
1457  throw exception::DimsInvalid("qpp::mket()");
1458 
1459  // check mask is a valid vector
1460  for (idx i = 0; i < N; ++i)
1461  if (mask[i] >= d)
1462  throw exception::SubsysMismatchDims("qpp::mket()");
1463  // END EXCEPTION CHECKS
1464 
1465  ket result = ket::Zero(D);
1466  std::vector<idx> dims(N, d);
1467  idx pos = multiidx2n(mask, dims);
1468  result(pos) = 1;
1469 
1470  return result;
1471 }
1472 
1473 // Idea taken from http://techblog.altplus.co.jp/entry/2017/11/08/130921
1483 template <char... Bits>
1484 ket operator"" _q() {
1485  constexpr idx n = sizeof...(Bits);
1486  constexpr char bits[n + 1] = {Bits..., '\0'};
1487  qpp::ket q = qpp::ket::Zero(std::pow(2, n));
1488  idx pos = 0;
1489 
1490  // EXCEPTION CHECKS
1491 
1492  // check valid multi-partite qubit state
1493  for (idx i = 0; i < n; ++i) {
1494  if (bits[i] != '0' && bits[i] != '1')
1495  throw exception::OutOfRange(R"xxx("qpp::operator "" _q())xxx");
1496  }
1497  // END EXCEPTION CHECKS
1498 
1499  pos = std::stoi(bits, nullptr, 2);
1500  q(pos) = 1;
1501 
1502  return q;
1503 }
1504 
1519 inline cmat mprj(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1520  idx N = mask.size();
1521 
1522  idx D = std::accumulate(std::begin(dims), std::end(dims),
1523  static_cast<idx>(1), std::multiplies<idx>());
1524 
1525  // EXCEPTION CHECKS
1526 
1527  // check zero size
1528  if (N == 0)
1529  throw exception::ZeroSize("qpp::mprj()");
1530  // check valid dims
1531  if (!internal::check_dims(dims))
1532  throw exception::DimsInvalid("qpp::mprj()");
1533  // check mask and dims have the same size
1534  if (mask.size() != dims.size())
1535  throw exception::SubsysMismatchDims("qpp::mprj()");
1536  // check mask is a valid vector
1537  for (idx i = 0; i < N; ++i)
1538  if (mask[i] >= dims[i])
1539  throw exception::SubsysMismatchDims("qpp::mprj()");
1540  // END EXCEPTION CHECKS
1541 
1542  cmat result = cmat::Zero(D, D);
1543  idx pos = multiidx2n(mask, dims);
1544  result(pos, pos) = 1;
1545 
1546  return result;
1547 }
1548 
1563 inline cmat mprj(const std::vector<idx>& mask, idx d = 2) {
1564  idx N = mask.size();
1565  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1566 
1567  // EXCEPTION CHECKS
1568 
1569  // check zero size
1570  if (N == 0)
1571  throw exception::ZeroSize("qpp::mprj()");
1572 
1573  // check valid dims
1574  if (d == 0)
1575  throw exception::DimsInvalid("qpp::mprj()");
1576 
1577  // check mask is a valid vector
1578  for (idx i = 0; i < N; ++i)
1579  if (mask[i] >= d)
1580  throw exception::SubsysMismatchDims("qpp::mprj()");
1581  // END EXCEPTION CHECKS
1582 
1583  cmat result = cmat::Zero(D, D);
1584  std::vector<idx> dims(N, d);
1585  idx pos = multiidx2n(mask, dims);
1586  result(pos, pos) = 1;
1587 
1588  return result;
1589 }
1590 
1599 template <typename InputIterator>
1600 std::vector<double> abssq(InputIterator first, InputIterator last) {
1601  std::vector<double> weights(std::distance(first, last));
1602  std::transform(first, last, std::begin(weights),
1603  [](cplx z) -> double { return std::norm(z); });
1604 
1605  return weights;
1606 }
1607 
1614 template <typename Container>
1615 std::vector<double>
1616 abssq(const Container& c,
1617  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1618 // we need the std::enable_if to SFINAE out Eigen expressions
1619 // that will otherwise match, instead of matching
1620 // the overload below:
1621 
1622 // template<typename Derived>
1623 // abssq(const Eigen::MatrixBase<Derived>& A)
1624 {
1625  return abssq(std::begin(c), std::end(c));
1626 }
1627 
1634 template <typename Derived>
1635 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A) {
1636  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1637 
1638  // EXCEPTION CHECKS
1639 
1640  // check zero-size
1642  throw exception::ZeroSize("qpp::abssq()");
1643  // END EXCEPTION CHECKS
1644 
1645  return abssq(rA.data(), rA.data() + rA.size());
1646 }
1647 
1656 template <typename InputIterator>
1657 typename std::iterator_traits<InputIterator>::value_type
1658 sum(InputIterator first, InputIterator last) {
1659  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1660 
1661  return std::accumulate(first, last, static_cast<value_type>(0));
1662 }
1663 
1671 template <typename Container>
1672 typename Container::value_type
1673 sum(const Container& c,
1674  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1675  return sum(std::begin(c), std::end(c));
1676 }
1677 
1686 template <typename InputIterator>
1687 typename std::iterator_traits<InputIterator>::value_type
1688 prod(InputIterator first, InputIterator last) {
1689  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1690 
1691  return std::accumulate(first, last, static_cast<value_type>(1),
1692  std::multiplies<value_type>());
1693 }
1694 
1702 template <typename Container>
1703 typename Container::value_type
1704 prod(const Container& c,
1705  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1706  return prod(std::begin(c), std::end(c));
1707 }
1708 
1721 template <typename Derived>
1722 dyn_col_vect<typename Derived::Scalar>
1723 rho2pure(const Eigen::MatrixBase<Derived>& A) {
1724  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1725 
1726  // EXCEPTION CHECKS
1727  // check zero-size
1729  throw exception::ZeroSize("qpp::rho2pure()");
1730  // check square matrix
1731  if (!internal::check_square_mat(rA))
1732  throw exception::MatrixNotSquare("qpp::rho2pure()");
1733  // END EXCEPTION CHECKS
1734 
1735  dyn_col_vect<double> tmp_evals = hevals(rA);
1736  cmat tmp_evects = hevects(rA);
1737  dyn_col_vect<typename Derived::Scalar> result =
1738  dyn_col_vect<typename Derived::Scalar>::Zero(rA.rows());
1739  // find the non-zero eigenvector
1740  // there is only one, assuming the state is pure
1741  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k) {
1742  if (std::abs(tmp_evals(k)) > eps) {
1743  result = tmp_evects.col(k);
1744  break;
1745  }
1746  }
1747 
1748  return result;
1749 }
1750 
1759 template <typename T>
1760 std::vector<T> complement(std::vector<T> subsys, idx N) {
1761  // EXCEPTION CHECKS
1762 
1763  if (N < subsys.size())
1764  throw exception::OutOfRange("qpp::complement()");
1765  // END EXCEPTION CHECKS
1766 
1767  std::vector<T> all(N);
1768  std::vector<T> subsys_bar(N - subsys.size());
1769 
1770  std::iota(std::begin(all), std::end(all), 0);
1771  std::sort(std::begin(subsys), std::end(subsys));
1772  std::set_difference(std::begin(all), std::end(all), std::begin(subsys),
1773  std::end(subsys), std::begin(subsys_bar));
1774 
1775  return subsys_bar;
1776 }
1777 
1788 template <typename Derived>
1789 std::vector<double> rho2bloch(const Eigen::MatrixBase<Derived>& A) {
1790  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1791 
1792  // EXCEPTION CHECKS
1793 
1794  // check qubit matrix
1796  throw exception::NotQubitMatrix("qpp::rho2bloch()");
1797  // END EXCEPTION CHECKS
1798 
1799  std::vector<double> result(3);
1800  cmat X(2, 2), Y(2, 2), Z(2, 2);
1801  X << 0, 1, 1, 0;
1802  Y << 0, -1_i, 1_i, 0;
1803  Z << 1, 0, 0, -1;
1804  result[0] = std::real(trace(rA * X));
1805  result[1] = std::real(trace(rA * Y));
1806  result[2] = std::real(trace(rA * Z));
1807 
1808  return result;
1809 }
1810 
1819 inline cmat bloch2rho(const std::vector<double>& r) {
1820  // EXCEPTION CHECKS
1821 
1822  // check 3-dimensional vector
1823  if (r.size() != 3)
1824  throw exception::CustomException("qpp::bloch2rho",
1825  "r is not a 3-dimensional vector!");
1826  // END EXCEPTION CHECKS
1827 
1828  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1829  X << 0, 1, 1, 0;
1830  Y << 0, -1_i, 1_i, 0;
1831  Z << 1, 0, 0, -1;
1832  Id2 << 1, 0, 0, 1;
1833 
1834  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1835 }
1836 
1837 } /* namespace qpp */
1838 
1839 #endif /* FUNCTIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:73
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
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:66
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
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: codes.h:35
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:87
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
std::size_t idx
Definition: experimental.h:43
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:317
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:236