27 #ifndef INSTRUMENTS_H_
28 #define INSTRUMENTS_H_
43 template<
typename Derived>
45 const Eigen::MatrixBase<Derived>& phi,
46 const Eigen::MatrixBase<Derived>& psi,
47 const std::vector<idx>& subsys,
48 const std::vector<idx>& dims)
88 std::vector<idx> subsys_dims(subsys.size());
89 for (
idx i = 0; i < subsys.size(); ++i)
90 subsys_dims[i] = dims[subsys[i]];
96 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
98 idx D =
static_cast<idx>(rpsi.rows());
99 idx Dsubsys_bar = D / Dsubsys;
102 idx Nsubsys = subsys.size();
103 idx Nsubsys_bar = N - Nsubsys;
111 std::vector<idx> subsys_bar =
complement(subsys, N);
112 std::copy(std::begin(subsys_bar), std::end(subsys_bar),
113 std::begin(Csubsys_bar));
115 for (
idx i = 0; i < N; ++i)
119 for (
idx i = 0; i < Nsubsys; ++i)
121 Csubsys[i] = subsys[i];
122 Cdimssubsys[i] = dims[subsys[i]];
124 for (
idx i = 0; i < Nsubsys_bar; ++i)
126 Cdimssubsys_bar[i] = dims[subsys_bar[i]];
129 auto worker = [=](
idx b) noexcept
130 ->
typename Derived::Scalar
138 Cdimssubsys_bar, Cmidxcolsubsys_bar);
140 for (
idx k = 0; k < Nsubsys_bar; ++k)
142 Cmidxrow[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
145 typename Derived::Scalar result = 0;
146 for (
idx a = 0; a < Dsubsys; ++a)
150 Cdimssubsys, Cmidxrowsubsys);
152 for (
idx k = 0; k < Nsubsys; ++k)
154 Cmidxrow[Csubsys[k]] = Cmidxrowsubsys[k];
159 result += std::conj(rphi(a)) * rpsi(i);
167 #pragma omp parallel for
168 #endif // WITH_OPENMP_
169 for (
idx m = 0; m < Dsubsys_bar; ++m)
170 result(m) = worker(m);
185 template<
typename Derived>
187 const Eigen::MatrixBase<Derived>& phi,
188 const Eigen::MatrixBase<Derived>& psi,
189 const std::vector<idx>& subsys,
206 std::vector<idx> dims(N, d);
207 return ip(phi, psi, subsys, dims);
220 template<
typename Derived>
221 std::tuple<idx, std::vector<double>, std::vector<cmat>>
222 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks)
237 if (Ks[0].rows() != rA.rows())
241 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
246 std::vector<double> prob(Ks.size());
248 std::vector<cmat> outstates(Ks.size());
253 for (
idx i = 0; i < Ks.size(); ++i)
255 outstates[i] = cmat::Zero(rA.rows(), rA.rows());
257 prob[i] = std::abs(
trace(tmp));
259 outstates[i] = tmp / prob[i];
265 for (
idx i = 0; i < Ks.size(); ++i)
267 outstates[i] = ket::Zero(rA.rows());
268 ket tmp = Ks[i] * rA;
270 prob[i] = std::pow(
norm(tmp), 2);
272 outstates[i] = tmp / std::sqrt(prob[i]);
279 std::discrete_distribution<idx> dd(std::begin(prob),
283 return std::make_tuple(result, prob, outstates);
298 template<
typename Derived>
299 std::tuple<idx, std::vector<double>, std::vector<cmat>>
301 const std::initializer_list<cmat>& Ks)
303 return measure(A, std::vector<cmat>(Ks));
316 template<
typename Derived>
317 std::tuple<idx, std::vector<double>, std::vector<cmat>>
333 if (U.rows() != rA.rows())
338 std::vector<cmat> Ks(U.rows());
339 for (
idx i = 0; i < static_cast<idx>(U.rows()); ++i)
340 Ks[i] = U.col(i) *
adjoint(U.col(i));
363 template<
typename Derived>
364 std::tuple<idx, std::vector<double>, std::vector<cmat>>
366 const std::vector<cmat>& Ks,
367 const std::vector<idx>& subsys,
368 const std::vector<idx>& dims)
370 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
393 std::vector<idx> subsys_dims(subsys.size());
394 for (
idx i = 0; i < subsys.size(); ++i)
395 subsys_dims[i] = dims[subsys[i]];
397 idx D =
prod(std::begin(dims), std::end(dims));
398 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
399 idx Dsubsys_bar = D / Dsubsys;
408 if (Dsubsys != static_cast<idx>(Ks[0].rows()))
412 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
418 std::vector<double> prob(Ks.size());
420 std::vector<cmat> outstates(Ks.size(), cmat::Zero(Dsubsys_bar, Dsubsys_bar));
425 for (
idx i = 0; i < Ks.size(); ++i)
427 cmat tmp =
apply(rA, Ks[i], subsys, dims);
428 tmp =
ptrace(tmp, subsys, dims);
429 prob[i] = std::abs(
trace(tmp));
434 outstates[i] = tmp / prob[i];
441 for (
idx i = 0; i < Ks.size(); ++i)
443 ket tmp =
apply(rA, Ks[i], subsys, dims);
444 prob[i] = std::pow(
norm(tmp), 2);
449 tmp /= std::sqrt(prob[i]);
450 outstates[i] =
ptrace(tmp, subsys, dims);
458 std::discrete_distribution<idx> dd(std::begin(prob),
462 return std::make_tuple(result, prob, outstates);
485 template<
typename Derived>
486 std::tuple<idx, std::vector<double>, std::vector<cmat>>
488 const std::initializer_list<cmat>& Ks,
489 const std::vector<idx>& subsys,
490 const std::vector<idx>& dims)
492 return measure(A, std::vector<cmat>(Ks), subsys, dims);
512 template<
typename Derived>
513 std::tuple<idx, std::vector<double>, std::vector<cmat>>
515 const std::vector<cmat>& Ks,
516 const std::vector<idx>& subsys,
519 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
534 std::vector<idx> dims(N, d);
536 return measure(rA, Ks, subsys, dims);
559 template<
typename Derived>
560 std::tuple<idx, std::vector<double>, std::vector<cmat>>
562 const std::initializer_list<cmat>& Ks,
563 const std::vector<idx>& subsys,
566 return measure(A, std::vector<cmat>(Ks), subsys, d);
587 template<
typename Derived>
588 std::tuple<idx, std::vector<double>, std::vector<cmat>>
591 const std::vector<idx>& subsys,
592 const std::vector<idx>& dims)
594 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
612 std::vector<idx> subsys_dims(subsys.size());
613 for (
idx i = 0; i < subsys.size(); ++i)
614 subsys_dims[i] = dims[subsys[i]];
616 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
623 if (Dsubsys != static_cast<idx>(V.rows()))
629 idx M =
static_cast<idx>(V.cols());
634 const ket& rpsi = A.derived();
640 std::vector<double> prob(M);
641 std::vector<cmat> outstates(M);
644 #pragma omp parallel for
645 #endif // WITH_OPENMP_
646 for (
idx i = 0; i < M; ++i)
647 outstates[i] =
ip(static_cast<const ket&>(V.col(i)),
650 for (
idx i = 0; i < M; ++i)
652 double tmp =
norm(outstates[i]);
663 std::discrete_distribution<idx> dd(std::begin(prob),
667 return std::make_tuple(result, prob, outstates);
677 std::vector<cmat> Ks(M);
678 for (
idx i = 0; i < M; ++i)
679 Ks[i] = V.col(i) *
adjoint(V.col(i));
681 return measure(rA, Ks, subsys, dims);
706 template<
typename Derived>
707 std::tuple<idx, std::vector<double>, std::vector<cmat>>
710 const std::vector<idx>& subsys,
713 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
728 std::vector<idx> dims(N, d);
730 return measure(rA, V, subsys, dims);
747 template<
typename Derived>
748 std::tuple<std::vector<idx>, double,
cmat>
750 std::vector<idx> subsys,
751 std::vector<idx> dims)
780 std::vector<idx> result;
785 std::sort(std::begin(subsys), std::end(subsys), std::greater<idx> {});
790 while (subsys.size() > 0)
796 result.push_back(std::get<0>(tmp));
797 prob *= std::get<1>(tmp)[std::get<0>(tmp)];
798 cA = std::get<2>(tmp)[std::get<0>(tmp)];
801 dims.erase(std::next(std::begin(dims), subsys[0]));
802 subsys.erase(std::begin(subsys));
805 std::reverse(std::begin(result), std::end(result));
810 return std::make_tuple(result, prob, cA);
827 template<
typename Derived>
828 std::tuple<std::vector<idx>, double,
cmat>
830 std::vector<idx> subsys,
833 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
848 std::vector<idx> dims(N, d);
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:112
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:74
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:231
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:67
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:141
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:51
dyn_mat< typename Derived1::Scalar > apply(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Applies the gate A to the part subsys of the multi-partite state vector or density matrix state...
Definition: operations.h:485
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:252
idx get_num_subsys(idx sz, idx d)
Definition: util.h:355
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:152
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1814
std::tuple< std::vector< idx >, double, cmat > measure_seq(const Eigen::MatrixBase< Derived > &A, std::vector< idx > subsys, std::vector< idx > dims)
Sequentially measures the part subsys of the multi-partite state vector or density matrix A in the co...
Definition: instruments.h:749
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
std::tuple< idx, std::vector< double >, std::vector< cmat > > measure(const Eigen::MatrixBase< Derived > &A, const std::vector< cmat > &Ks)
Measures the state A using the set of Kraus operators Ks.
Definition: instruments.h:222
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
void n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:48
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:44
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1209
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:127
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:90
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:90
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:125
std::size_t idx
Non-negative integer index.
Definition: types.h:36
idx multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:61