32 #ifndef INSTRUMENTS_H_ 33 #define INSTRUMENTS_H_ 46 template <
typename Derived>
47 dyn_col_vect<typename Derived::Scalar>
48 ip(
const Eigen::MatrixBase<Derived>& phi,
const Eigen::MatrixBase<Derived>& psi,
49 const std::vector<idx>& subsys,
const std::vector<idx>& dims) {
84 std::vector<idx> subsys_dims(subsys.size());
85 for (
idx i = 0; i < subsys.size(); ++i)
86 subsys_dims[i] = dims[subsys[i]];
91 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
93 idx D =
static_cast<idx>(rpsi.rows());
94 idx Dsubsys_bar = D / Dsubsys;
97 idx n_subsys = subsys.size();
98 idx n_subsys_bar = n - n_subsys;
106 std::vector<idx> subsys_bar =
complement(subsys, n);
107 std::copy(std::begin(subsys_bar), std::end(subsys_bar),
108 std::begin(Csubsys_bar));
110 for (
idx i = 0; i < n; ++i) {
113 for (
idx i = 0; i < n_subsys; ++i) {
114 Csubsys[i] = subsys[i];
115 Cdimssubsys[i] = dims[subsys[i]];
117 for (
idx i = 0; i < n_subsys_bar; ++i) {
118 Cdimssubsys_bar[i] = dims[subsys_bar[i]];
121 auto worker = [&](
idx b) noexcept->typename Derived::Scalar {
130 for (
idx k = 0; k < n_subsys_bar; ++k) {
131 Cmidxrow[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
134 typename Derived::Scalar result = 0;
135 for (
idx a = 0; a < Dsubsys; ++a) {
139 for (
idx k = 0; k < n_subsys; ++k) {
140 Cmidxrow[Csubsys[k]] = Cmidxrowsubsys[k];
145 result += std::conj(rphi(a)) * rpsi(i);
151 dyn_col_vect<typename Derived::Scalar> result(Dsubsys_bar);
153 #pragma omp parallel for 154 #endif // WITH_OPENMP_ 155 for (
idx m = 0; m < Dsubsys_bar; ++m)
156 result(m) = worker(m);
171 template <
typename Derived>
172 dyn_col_vect<typename Derived::Scalar>
173 ip(
const Eigen::MatrixBase<Derived>& phi,
const Eigen::MatrixBase<Derived>& psi,
174 const std::vector<idx>& subsys,
idx d = 2) {
189 std::vector<idx> dims(n, d);
190 return ip(phi, psi, subsys, dims);
204 template <
typename Derived>
205 std::tuple<idx, std::vector<double>, std::vector<cmat>>
206 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks) {
220 if (Ks[0].rows() != rA.rows())
222 for (
auto&& elem : Ks)
223 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
228 std::vector<double> prob(Ks.size());
230 std::vector<cmat> outstates(Ks.size());
235 for (
idx i = 0; i < Ks.size(); ++i) {
236 outstates[i] = cmat::Zero(rA.rows(), rA.rows());
238 prob[i] = std::abs(
trace(tmp));
240 outstates[i] = tmp / prob[i];
246 for (
idx i = 0; i < Ks.size(); ++i) {
247 outstates[i] = ket::Zero(rA.rows());
248 ket tmp = Ks[i] * rA;
250 prob[i] = std::pow(
norm(tmp), 2);
252 outstates[i] = tmp / std::sqrt(prob[i]);
258 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
260 #ifdef NO_THREAD_LOCAL_ 261 RandomDevices::get_instance().get_prng();
263 RandomDevices::get_thread_local_instance().get_prng();
265 idx result = dd(gen);
267 return std::make_tuple(result, prob, outstates);
283 template <
typename Derived>
284 std::tuple<idx, std::vector<double>, std::vector<cmat>>
286 const std::initializer_list<cmat>& Ks) {
287 return measure(A, std::vector<cmat>(Ks));
300 template <
typename Derived>
301 std::tuple<idx, std::vector<double>, std::vector<cmat>>
316 if (U.rows() != rA.rows())
320 std::vector<cmat> Ks(U.rows());
321 for (
idx i = 0; i < static_cast<idx>(U.rows()); ++i)
322 Ks[i] = U.col(i) *
adjoint(U.col(i));
344 template <
typename Derived>
345 std::tuple<idx, std::vector<double>, std::vector<cmat>>
346 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
347 const std::vector<idx>& target,
const std::vector<idx>& dims) {
348 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
364 std::vector<idx> subsys_dims(target.size());
365 for (
idx i = 0; i < target.size(); ++i)
366 subsys_dims[i] = dims[target[i]];
368 idx D =
prod(std::begin(dims), std::end(dims));
369 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
370 idx Dsubsys_bar = D / Dsubsys;
377 if (Dsubsys != static_cast<idx>(Ks[0].rows()))
379 for (
auto&& elem : Ks)
380 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
385 std::vector<double> prob(Ks.size());
387 std::vector<cmat> outstates(Ks.size(),
388 cmat::Zero(Dsubsys_bar, Dsubsys_bar));
396 for (
idx i = 0; i < Ks.size(); ++i) {
397 cmat tmp =
apply(rA, Ks[i], target, dims);
398 tmp =
ptrace(tmp, target, dims);
399 prob[i] = std::abs(
trace(tmp));
403 outstates[i] = tmp / prob[i];
413 for (
idx i = 0; i < Ks.size(); ++i) {
414 ket tmp =
apply(rA, Ks[i], target, dims);
415 prob[i] = std::pow(
norm(tmp), 2);
419 tmp /= std::sqrt(prob[i]);
420 outstates[i] =
ptrace(tmp, target, dims);
427 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
429 #ifdef NO_THREAD_LOCAL_ 430 RandomDevices::get_instance().get_prng();
432 RandomDevices::get_thread_local_instance().get_prng();
434 idx result = dd(gen);
436 return std::make_tuple(result, prob, outstates);
458 template <
typename Derived>
459 std::tuple<idx, std::vector<double>, std::vector<cmat>>
461 const std::initializer_list<cmat>& Ks,
const std::vector<idx>& target,
462 const std::vector<idx>& dims) {
463 return measure(A, std::vector<cmat>(Ks), target, dims);
482 template <
typename Derived>
483 std::tuple<idx, std::vector<double>, std::vector<cmat>>
484 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
485 const std::vector<idx>& target,
idx d = 2) {
486 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
500 std::vector<idx> dims(n, d);
502 return measure(rA, Ks, target, dims);
524 template <
typename Derived>
525 std::tuple<idx, std::vector<double>, std::vector<cmat>>
527 const std::initializer_list<cmat>& Ks,
const std::vector<idx>& target,
529 return measure(A, std::vector<cmat>(Ks), target, d);
550 template <
typename Derived>
551 std::tuple<idx, std::vector<double>, std::vector<cmat>>
553 const std::vector<idx>& target,
const std::vector<idx>& dims) {
554 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
570 std::vector<idx> subsys_dims(target.size());
571 for (
idx i = 0; i < target.size(); ++i)
572 subsys_dims[i] = dims[target[i]];
574 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
579 if (Dsubsys != static_cast<idx>(V.rows()))
584 idx M =
static_cast<idx>(V.cols());
588 const ket& rpsi = A.derived();
593 std::vector<double> prob(M);
594 std::vector<cmat> outstates(M);
597 #pragma omp parallel for 598 #endif // WITH_OPENMP_ 599 for (
idx i = 0; i < M; ++i)
601 ip(static_cast<const ket&>(V.col(i)), rpsi, target, dims);
603 for (
idx i = 0; i < M; ++i) {
604 double tmp =
norm(outstates[i]);
614 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
616 #ifdef NO_THREAD_LOCAL_ 617 RandomDevices::get_instance().get_prng();
619 RandomDevices::get_thread_local_instance().get_prng();
621 idx result = dd(gen);
623 return std::make_tuple(result, prob, outstates);
631 std::vector<cmat> Ks(M);
632 for (
idx i = 0; i < M; ++i)
633 Ks[i] = V.col(i) *
adjoint(V.col(i));
635 return measure(rA, Ks, target, dims);
659 template <
typename Derived>
660 std::tuple<idx, std::vector<double>, std::vector<cmat>>
662 const std::vector<idx>& target,
idx d = 2) {
663 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
677 std::vector<idx> dims(n, d);
679 return measure(rA, V, target, dims);
695 template <
typename Derived>
696 std::tuple<std::vector<idx>, double,
cmat>
697 measure_seq(
const Eigen::MatrixBase<Derived>& A, std::vector<idx> target,
698 std::vector<idx> dims) {
731 std::vector<idx> result;
736 std::sort(std::begin(target), std::end(target), std::greater<idx>{});
739 while (target.size() > 0) {
740 auto tmp =
measure(rA, Gates::get_instance().Id(dims[target[0]]),
742 result.push_back(std::get<0>(tmp));
743 prob *= std::get<1>(tmp)[std::get<0>(tmp)];
744 rA = std::get<2>(tmp)[std::get<0>(tmp)];
747 dims.erase(std::next(std::begin(dims), target[0]));
748 target.erase(std::begin(target));
751 std::reverse(std::begin(result), std::end(result));
753 return std::make_tuple(result, prob, rA);
769 template <
typename Derived>
770 std::tuple<std::vector<idx>, double,
cmat>
771 measure_seq(
const Eigen::MatrixBase<Derived>& A, std::vector<idx> target,
773 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
787 std::vector<idx> dims(n, d);
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimensions not equal exception.
Definition: exception.h:284
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:67
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:219
std::tuple< idx, std::vector< double >, std::vector< cmat > > measure(const Eigen::MatrixBase< Derived > &A, const cmat &V, const std::vector< idx > &target, idx d=2)
Measures the part target of the multi-partite state vector or density matrix A in the orthonormal bas...
Definition: instruments.h:661
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
dyn_col_vect< typename Derived::Scalar > ip(const Eigen::MatrixBase< Derived > &phi, const Eigen::MatrixBase< Derived > &psi, const std::vector< idx > &subsys, idx d=2)
Generalized inner product.
Definition: instruments.h:173
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:228
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
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
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Matrix is not a column vector exception.
Definition: exception.h:164
Quantum++ main namespace.
Definition: circuits.h:35
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:248
Matrix is not square nor column vector exception.
Definition: exception.h:209
Invalid dimension(s) exception.
Definition: exception.h:269
idx get_num_subsys(idx D, idx d)
Definition: util.h:370
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1187
dyn_mat< typename Derived1::Scalar > apply(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &target, const std::vector< idx > &dims)
Applies the gate A to the part target of the multi-partite state vector or density matrix state...
Definition: operations.h:444
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:89
Dimension(s) mismatch column vector size exception.
Definition: exception.h:316
std::tuple< std::vector< idx >, double, cmat > measure_seq(const Eigen::MatrixBase< Derived > &A, std::vector< idx > target, idx d=2)
Sequentially measures the part target of the multi-partite state vector or density matrix A in the co...
Definition: instruments.h:771
dyn_col_vect< typename Derived::Scalar > ip(const Eigen::MatrixBase< Derived > &phi, const Eigen::MatrixBase< Derived > &psi, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Generalized inner product.
Definition: instruments.h:48
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:130
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
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
std::size_t idx
Non-negative integer index, make sure you use an unsigned type.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:149
std::vector< idx > complement(std::vector< idx > subsys, idx n)
Constructs the complement of a subsystem vector.
Definition: functions.h:1780
Object has zero size exception.
Definition: exception.h:134