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 Dsubsysbar = D / Dsubsys;
102 idx nsubsys = subsys.size();
103 idx nsubsysbar = 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(Csubsysbar));
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 < nsubsysbar; ++i )
126 Cdimssubsysbar[i] = dims[subsys_bar[i]];
129 auto worker = [=](
idx b) noexcept
130 ->
typename Derived::Scalar
138 Cdimssubsysbar, Cmidxcolsubsysbar);
140 for (
idx k = 0; k < nsubsysbar; ++k )
142 Cmidxrow[Csubsysbar[k]] = Cmidxcolsubsysbar[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
169 for (
idx m = 0; m < Dsubsysbar; ++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 static_cast<idx>(std::llround(std::log2(rpsi.rows()) /
208 std::vector<idx> dims(n, d);
209 return ip(phi, psi, subsys, dims);
222 template<
typename Derived>
223 std::tuple<idx, std::vector<double>, std::vector<cmat>>
224 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks)
235 if ( Ks.size() == 0 )
239 if ( Ks[0].rows() != rA.rows())
242 for (
auto&& it : Ks )
243 if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
248 std::vector<double> prob(Ks.size());
250 std::vector<cmat> outstates(Ks.size());
255 for (
idx i = 0; i < Ks.size(); ++i )
257 outstates[i] = cmat::Zero(rA.rows(), rA.rows());
259 prob[i] = std::abs(
trace(tmp));
261 outstates[i] = tmp / prob[i];
267 for (
idx i = 0; i < Ks.size(); ++i )
269 outstates[i] = ket::Zero(rA.rows());
270 ket tmp = Ks[i] * rA;
272 prob[i] = std::pow(
norm(tmp), 2);
274 outstates[i] = tmp / std::sqrt(prob[i]);
282 std::discrete_distribution<idx> dd(std::begin(prob),
286 return std::make_tuple(result, prob, outstates);
301 template<
typename Derived>
302 std::tuple<idx, std::vector<double>, std::vector<cmat>>
304 const std::initializer_list<cmat>& Ks)
306 return measure(A, std::vector<cmat>(Ks));
319 template<
typename Derived>
320 std::tuple<idx, std::vector<double>, std::vector<cmat>>
336 if ( U.rows() != rA.rows())
341 std::vector<cmat> Ks(U.rows());
342 for (
idx i = 0; i < static_cast<idx>(U.rows()); ++i )
343 Ks[i] = U.col(i) *
adjoint(U.col(i));
366 template<
typename Derived>
367 std::tuple<idx, std::vector<double>, std::vector<cmat>>
369 const std::vector<cmat>& Ks,
370 const std::vector<idx>& subsys,
371 const std::vector<idx>& dims)
373 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
396 std::vector<idx> subsys_dims(subsys.size());
397 for (
idx i = 0; i < subsys.size(); ++i )
398 subsys_dims[i] = dims[subsys[i]];
400 idx D =
prod(std::begin(dims), std::end(dims));
401 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
402 idx Dsubsysbar = D / Dsubsys;
405 if ( Ks.size() == 0 )
411 if ( Dsubsys != static_cast<idx>(Ks[0].rows()))
414 for (
auto&& it : Ks )
415 if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
421 std::vector<double> prob(Ks.size());
423 std::vector<cmat> outstates(Ks.size(), cmat::Zero(Dsubsysbar, Dsubsysbar));
428 for (
idx i = 0; i < Ks.size(); ++i )
430 cmat tmp =
apply(rA, Ks[i], subsys, dims);
431 tmp =
ptrace(tmp, subsys, dims);
432 prob[i] = std::abs(
trace(tmp));
437 outstates[i] = tmp / prob[i];
444 for (
idx i = 0; i < Ks.size(); ++i )
446 ket tmp =
apply(rA, Ks[i], subsys, dims);
447 prob[i] = std::pow(
norm(tmp), 2);
452 tmp /= std::sqrt(prob[i]);
453 outstates[i] =
ptrace(tmp, subsys, dims);
462 std::discrete_distribution<idx> dd(std::begin(prob),
466 return std::make_tuple(result, prob, outstates);
489 template<
typename Derived>
490 std::tuple<idx, std::vector<double>, std::vector<cmat>>
492 const std::initializer_list<cmat>& Ks,
493 const std::vector<idx>& subsys,
494 const std::vector<idx>& dims)
496 return measure(A, std::vector<cmat>(Ks), subsys, dims);
516 template<
typename Derived>
517 std::tuple<idx, std::vector<double>, std::vector<cmat>>
519 const std::vector<cmat>& Ks,
520 const std::vector<idx>& subsys,
523 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
538 static_cast<idx>(std::llround(std::log2(rA.rows()) /
540 std::vector<idx> dims(n, d);
542 return measure(rA, Ks, subsys, dims);
565 template<
typename Derived>
566 std::tuple<idx, std::vector<double>, std::vector<cmat>>
568 const std::initializer_list<cmat>& Ks,
569 const std::vector<idx>& subsys,
572 return measure(A, std::vector<cmat>(Ks), subsys, d);
593 template<
typename Derived>
594 std::tuple<idx, std::vector<double>, std::vector<cmat>>
598 const std::vector<idx>& subsys,
599 const std::vector<idx>& dims)
601 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
619 std::vector<idx> subsys_dims(subsys.size());
620 for (
idx i = 0; i < subsys.size(); ++i )
621 subsys_dims[i] = dims[subsys[i]];
623 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
630 if ( Dsubsys != static_cast<idx>(V.rows()))
636 idx M =
static_cast<idx>(V.cols());
641 const ket& rpsi = A.derived();
647 std::vector<double> prob(M);
648 std::vector<cmat> outstates(M);
651 #pragma omp parallel for
653 for (
idx i = 0; i < M; ++i )
654 outstates[i] =
ip(static_cast<const ket&>(V.col(i)),
657 for (
idx i = 0; i < M; ++i )
659 double tmp =
norm(outstates[i]);
670 std::discrete_distribution<idx> dd(std::begin(prob),
674 return std::make_tuple(result, prob, outstates);
684 std::vector<cmat> Ks(M);
685 for (
idx i = 0; i < M; ++i )
686 Ks[i] = V.col(i) *
adjoint(V.col(i));
688 return measure(rA, Ks, subsys, dims);
713 template<
typename Derived>
714 std::tuple<idx, std::vector<double>, std::vector<cmat>>
718 const std::vector<idx>& subsys,
721 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
736 static_cast<idx>(std::llround(std::log2(rA.rows()) /
738 std::vector<idx> dims(n, d);
740 return measure(rA, V, subsys, dims);
757 template<
typename Derived>
758 std::tuple<std::vector<idx>, double,
cmat>
760 std::vector<idx> subsys,
761 std::vector<idx> dims)
763 typename Eigen::MatrixBase<Derived>::EvalReturnType cA = A;
786 std::vector<idx> result;
791 std::sort(std::begin(subsys), std::end(subsys), std::greater<idx>{});
796 while (subsys.size() > 0)
802 result.push_back(std::get<0>(tmp));
803 prob *= std::get<1>(tmp)[std::get<0>(tmp)];
804 cA = std::get<2>(tmp)[std::get<0>(tmp)];
807 dims.erase(std::next(std::begin(dims), subsys[0]));
808 subsys.erase(std::begin(subsys));
811 std::reverse(std::begin(result), std::end(result));
817 return std::make_tuple(result, prob, cA);
834 template<
typename Derived>
835 std::tuple<std::vector<idx>, double,
cmat>
837 std::vector<idx> subsys,
840 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
855 static_cast<idx>(std::llround(std::log2(rA.rows()) /
857 std::vector<idx> dims(n, d);
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
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_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
bool _check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:152
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:83
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:56
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:492
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:252
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1818
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:759
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:224
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:61
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
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
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:112
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:1218
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:95
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:66
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:48
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:125