27 #ifndef INSTRUMENTS_H_ 28 #define INSTRUMENTS_H_ 42 template<
typename Derived>
44 const Eigen::MatrixBase<Derived>& phi,
45 const Eigen::MatrixBase<Derived>& psi,
46 const std::vector<idx>& subsys,
47 const std::vector<idx>& dims)
83 std::vector<idx> subsys_dims(subsys.size());
84 for (
idx i = 0; i < subsys.size(); ++i)
85 subsys_dims[i] = dims[subsys[i]];
90 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
92 idx D =
static_cast<idx>(rpsi.rows());
93 idx Dsubsys_bar = D / Dsubsys;
96 idx Nsubsys = subsys.size();
97 idx Nsubsys_bar = N - Nsubsys;
105 std::vector<idx> subsys_bar =
complement(subsys, N);
106 std::copy(std::begin(subsys_bar), std::end(subsys_bar),
107 std::begin(Csubsys_bar));
109 for (
idx i = 0; i < N; ++i)
113 for (
idx i = 0; i < Nsubsys; ++i)
115 Csubsys[i] = subsys[i];
116 Cdimssubsys[i] = dims[subsys[i]];
118 for (
idx i = 0; i < Nsubsys_bar; ++i)
120 Cdimssubsys_bar[i] = dims[subsys_bar[i]];
123 auto worker = [&](
idx b) noexcept
124 ->
typename Derived::Scalar
132 Cdimssubsys_bar, Cmidxcolsubsys_bar);
134 for (
idx k = 0; k < Nsubsys_bar; ++k)
136 Cmidxrow[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
139 typename Derived::Scalar result = 0;
140 for (
idx a = 0; a < Dsubsys; ++a)
144 Cdimssubsys, Cmidxrowsubsys);
146 for (
idx k = 0; k < Nsubsys; ++k)
148 Cmidxrow[Csubsys[k]] = Cmidxrowsubsys[k];
153 result += std::conj(rphi(a)) * rpsi(i);
161 #pragma omp parallel for 162 #endif // WITH_OPENMP_ 163 for (
idx m = 0; m < Dsubsys_bar; ++m)
164 result(m) = worker(m);
179 template<
typename Derived>
181 const Eigen::MatrixBase<Derived>& phi,
182 const Eigen::MatrixBase<Derived>& psi,
183 const std::vector<idx>& subsys,
200 std::vector<idx> dims(N, d);
201 return ip(phi, psi, subsys, dims);
214 template<
typename Derived>
215 std::tuple<idx, std::vector<double>, std::vector<cmat>>
216 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks)
231 if (Ks[0].rows() != rA.rows())
234 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
239 std::vector<double> prob(Ks.size());
241 std::vector<cmat> outstates(Ks.size());
246 for (
idx i = 0; i < Ks.size(); ++i)
248 outstates[i] = cmat::Zero(rA.rows(), rA.rows());
250 prob[i] = std::abs(
trace(tmp));
252 outstates[i] = tmp / prob[i];
258 for (
idx i = 0; i < Ks.size(); ++i)
260 outstates[i] = ket::Zero(rA.rows());
261 ket tmp = Ks[i] * rA;
263 prob[i] = std::pow(
norm(tmp), 2);
265 outstates[i] = tmp / std::sqrt(prob[i]);
271 std::discrete_distribution<idx> dd(std::begin(prob),
275 return std::make_tuple(result, prob, outstates);
290 template<
typename Derived>
291 std::tuple<idx, std::vector<double>, std::vector<cmat>>
293 const std::initializer_list<cmat>& Ks)
295 return measure(A, std::vector<cmat>(Ks));
308 template<
typename Derived>
309 std::tuple<idx, std::vector<double>, std::vector<cmat>>
325 if (U.rows() != rA.rows())
329 std::vector<cmat> Ks(U.rows());
330 for (
idx i = 0; i < static_cast<idx>(U.rows()); ++i)
331 Ks[i] = U.col(i) *
adjoint(U.col(i));
354 template<
typename Derived>
355 std::tuple<idx, std::vector<double>, std::vector<cmat>>
357 const std::vector<cmat>& Ks,
358 const std::vector<idx>& subsys,
359 const std::vector<idx>& dims)
361 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
381 std::vector<idx> subsys_dims(subsys.size());
382 for (
idx i = 0; i < subsys.size(); ++i)
383 subsys_dims[i] = dims[subsys[i]];
385 idx D =
prod(std::begin(dims), std::end(dims));
386 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
387 idx Dsubsys_bar = D / Dsubsys;
394 if (Dsubsys != static_cast<idx>(Ks[0].rows()))
397 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
402 std::vector<double> prob(Ks.size());
404 std::vector<cmat> outstates(Ks.size(),
405 cmat::Zero(Dsubsys_bar, Dsubsys_bar));
410 for (
idx i = 0; i < Ks.size(); ++i)
412 cmat tmp =
apply(rA, Ks[i], subsys, dims);
413 tmp =
ptrace(tmp, subsys, dims);
414 prob[i] = std::abs(
trace(tmp));
419 outstates[i] = tmp / prob[i];
426 for (
idx i = 0; i < Ks.size(); ++i)
428 ket tmp =
apply(rA, Ks[i], subsys, dims);
429 prob[i] = std::pow(
norm(tmp), 2);
434 tmp /= std::sqrt(prob[i]);
435 outstates[i] =
ptrace(tmp, subsys, dims);
442 std::discrete_distribution<idx> dd(std::begin(prob),
446 return std::make_tuple(result, prob, outstates);
469 template<
typename Derived>
470 std::tuple<idx, std::vector<double>, std::vector<cmat>>
472 const std::initializer_list<cmat>& Ks,
473 const std::vector<idx>& subsys,
474 const std::vector<idx>& dims)
476 return measure(A, std::vector<cmat>(Ks), subsys, dims);
496 template<
typename Derived>
497 std::tuple<idx, std::vector<double>, std::vector<cmat>>
499 const std::vector<cmat>& Ks,
500 const std::vector<idx>& subsys,
503 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
518 std::vector<idx> dims(N, d);
520 return measure(rA, Ks, subsys, dims);
543 template<
typename Derived>
544 std::tuple<idx, std::vector<double>, std::vector<cmat>>
546 const std::initializer_list<cmat>& Ks,
547 const std::vector<idx>& subsys,
550 return measure(A, std::vector<cmat>(Ks), subsys, d);
571 template<
typename Derived>
572 std::tuple<idx, std::vector<double>, std::vector<cmat>>
575 const std::vector<idx>& subsys,
576 const std::vector<idx>& dims)
578 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
595 std::vector<idx> subsys_dims(subsys.size());
596 for (
idx i = 0; i < subsys.size(); ++i)
597 subsys_dims[i] = dims[subsys[i]];
599 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
604 if (Dsubsys != static_cast<idx>(V.rows()))
609 idx M =
static_cast<idx>(V.cols());
614 const ket& rpsi = A.derived();
619 std::vector<double> prob(M);
620 std::vector<cmat> outstates(M);
623 #pragma omp parallel for 624 #endif // WITH_OPENMP_ 625 for (
idx i = 0; i < M; ++i)
626 outstates[i] =
ip(static_cast<const ket&>(V.col(i)),
629 for (
idx i = 0; i < M; ++i)
631 double tmp =
norm(outstates[i]);
642 std::discrete_distribution<idx> dd(std::begin(prob),
646 return std::make_tuple(result, prob, outstates);
655 std::vector<cmat> Ks(M);
656 for (
idx i = 0; i < M; ++i)
657 Ks[i] = V.col(i) *
adjoint(V.col(i));
659 return measure(rA, Ks, subsys, dims);
683 template<
typename Derived>
684 std::tuple<idx, std::vector<double>, std::vector<cmat>>
687 const std::vector<idx>& subsys,
690 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
705 std::vector<idx> dims(N, d);
707 return measure(rA, V, subsys, dims);
724 template<
typename Derived>
725 std::tuple<std::vector<idx>, double,
cmat>
727 std::vector<idx> subsys,
728 std::vector<idx> dims)
766 std::vector<idx> result;
771 std::sort(std::begin(subsys), std::end(subsys), std::greater<idx> {});
774 while (subsys.size() > 0)
780 result.push_back(std::get<0>(tmp));
781 prob *= std::get<1>(tmp)[std::get<0>(tmp)];
782 cA = std::get<2>(tmp)[std::get<0>(tmp)];
785 dims.erase(std::next(std::begin(dims), subsys[0]));
786 subsys.erase(std::begin(subsys));
789 std::reverse(std::begin(result), std::end(result));
791 return std::make_tuple(result, prob, cA);
808 template<
typename Derived>
809 std::tuple<std::vector<idx>, double,
cmat>
811 std::vector<idx> subsys,
814 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
829 std::vector<idx> dims(N, d);
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
Dimensions not equal exception.
Definition: exception.h:304
Dimension(s) mismatch matrix size exception.
Definition: exception.h:322
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:71
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:221
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:122
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:227
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:64
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:77
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:158
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:70
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:40
Subsystems mismatch dimensions exception.
Definition: exception.h:395
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:50
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:473
Matrix is not a column vector exception.
Definition: exception.h:168
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:248
idx get_num_subsys(idx sz, idx d)
Definition: util.h:391
Matrix is not square nor column vector exception.
Definition: exception.h:219
Invalid dimension(s) exception.
Definition: exception.h:287
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1811
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:174
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:726
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:81
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:216
Dimension(s) mismatch column vector size exception.
Definition: exception.h:340
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:43
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:1241
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:89
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:89
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:142
std::size_t idx
Non-negative integer index.
Definition: types.h:35
Matrix is not square exception.
Definition: exception.h:151
Object has zero size exception.
Definition: exception.h:134