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 Nsubsys = subsys.size();
98 idx Nsubsys_bar = N - Nsubsys;
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 < Nsubsys; ++i) {
114 Csubsys[i] = subsys[i];
115 Cdimssubsys[i] = dims[subsys[i]];
117 for (
idx i = 0; i < Nsubsys_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 < Nsubsys_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 < Nsubsys; ++k) {
140 Cmidxrow[Csubsys[k]] = Cmidxrowsubsys[k];
145 result += std::conj(rphi(a)) * rpsi(i);
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>
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);
203 template <
typename Derived>
204 std::tuple<idx, std::vector<double>, std::vector<cmat>>
205 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks) {
219 if (Ks[0].rows() != rA.rows())
222 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
227 std::vector<double> prob(Ks.size());
229 std::vector<cmat> outstates(Ks.size());
234 for (
idx i = 0; i < Ks.size(); ++i) {
235 outstates[i] = cmat::Zero(rA.rows(), rA.rows());
237 prob[i] = std::abs(
trace(tmp));
239 outstates[i] = tmp / prob[i];
245 for (
idx i = 0; i < Ks.size(); ++i) {
246 outstates[i] = ket::Zero(rA.rows());
247 ket tmp = Ks[i] * rA;
249 prob[i] = std::pow(
norm(tmp), 2);
251 outstates[i] = tmp / std::sqrt(prob[i]);
257 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
260 return std::make_tuple(result, prob, outstates);
275 template <
typename Derived>
276 std::tuple<idx, std::vector<double>, std::vector<cmat>>
278 const std::initializer_list<cmat>& Ks) {
279 return measure(A, std::vector<cmat>(Ks));
292 template <
typename Derived>
293 std::tuple<idx, std::vector<double>, std::vector<cmat>>
308 if (U.rows() != rA.rows())
312 std::vector<cmat> Ks(U.rows());
313 for (
idx i = 0; i < static_cast<idx>(U.rows()); ++i)
314 Ks[i] = U.col(i) *
adjoint(U.col(i));
337 template <
typename Derived>
338 std::tuple<idx, std::vector<double>, std::vector<cmat>>
339 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
340 const std::vector<idx>& subsys,
const std::vector<idx>& dims) {
341 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
361 std::vector<idx> subsys_dims(subsys.size());
362 for (
idx i = 0; i < subsys.size(); ++i)
363 subsys_dims[i] = dims[subsys[i]];
365 idx D =
prod(std::begin(dims), std::end(dims));
366 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
367 idx Dsubsys_bar = D / Dsubsys;
374 if (Dsubsys != static_cast<idx>(Ks[0].rows()))
377 if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
382 std::vector<double> prob(Ks.size());
384 std::vector<cmat> outstates(Ks.size(),
385 cmat::Zero(Dsubsys_bar, Dsubsys_bar));
390 for (
idx i = 0; i < Ks.size(); ++i) {
391 cmat tmp =
apply(rA, Ks[i], subsys, dims);
392 tmp =
ptrace(tmp, subsys, dims);
393 prob[i] = std::abs(
trace(tmp));
397 outstates[i] = tmp / prob[i];
404 for (
idx i = 0; i < Ks.size(); ++i) {
405 ket tmp =
apply(rA, Ks[i], subsys, dims);
406 prob[i] = std::pow(
norm(tmp), 2);
410 tmp /= std::sqrt(prob[i]);
411 outstates[i] =
ptrace(tmp, subsys, dims);
418 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
421 return std::make_tuple(result, prob, outstates);
444 template <
typename Derived>
445 std::tuple<idx, std::vector<double>, std::vector<cmat>>
447 const std::initializer_list<cmat>& Ks,
const std::vector<idx>& subsys,
448 const std::vector<idx>& dims) {
449 return measure(A, std::vector<cmat>(Ks), subsys, dims);
469 template <
typename Derived>
470 std::tuple<idx, std::vector<double>, std::vector<cmat>>
471 measure(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
472 const std::vector<idx>& subsys,
idx d = 2) {
473 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
487 std::vector<idx> dims(N, d);
489 return measure(rA, Ks, subsys, dims);
512 template <
typename Derived>
513 std::tuple<idx, std::vector<double>, std::vector<cmat>>
515 const std::initializer_list<cmat>& Ks,
const std::vector<idx>& subsys,
517 return measure(A, std::vector<cmat>(Ks), subsys, d);
538 template <
typename Derived>
539 std::tuple<idx, std::vector<double>, std::vector<cmat>>
541 const std::vector<idx>& subsys,
const std::vector<idx>& dims) {
542 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
558 std::vector<idx> subsys_dims(subsys.size());
559 for (
idx i = 0; i < subsys.size(); ++i)
560 subsys_dims[i] = dims[subsys[i]];
562 idx Dsubsys =
prod(std::begin(subsys_dims), std::end(subsys_dims));
567 if (Dsubsys != static_cast<idx>(V.rows()))
572 idx M =
static_cast<idx>(V.cols());
576 const ket& rpsi = A.derived();
581 std::vector<double> prob(M);
582 std::vector<cmat> outstates(M);
585 #pragma omp parallel for 586 #endif // WITH_OPENMP_ 587 for (
idx i = 0; i < M; ++i)
589 ip(static_cast<const ket&>(V.col(i)), rpsi, subsys, dims);
591 for (
idx i = 0; i < M; ++i) {
592 double tmp =
norm(outstates[i]);
602 std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
605 return std::make_tuple(result, prob, outstates);
613 std::vector<cmat> Ks(M);
614 for (
idx i = 0; i < M; ++i)
615 Ks[i] = V.col(i) *
adjoint(V.col(i));
617 return measure(rA, Ks, subsys, dims);
641 template <
typename Derived>
642 std::tuple<idx, std::vector<double>, std::vector<cmat>>
644 const std::vector<idx>& subsys,
idx d = 2) {
645 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
659 std::vector<idx> dims(N, d);
661 return measure(rA, V, subsys, dims);
678 template <
typename Derived>
679 std::tuple<std::vector<idx>, double,
cmat>
680 measure_seq(
const Eigen::MatrixBase<Derived>& A, std::vector<idx> subsys,
681 std::vector<idx> dims) {
715 std::vector<idx> result;
720 std::sort(std::begin(subsys), std::end(subsys), std::greater<idx>{});
723 while (subsys.size() > 0) {
726 result.push_back(std::get<0>(tmp));
727 prob *= std::get<1>(tmp)[std::get<0>(tmp)];
728 cA = std::get<2>(tmp)[std::get<0>(tmp)];
731 dims.erase(std::next(std::begin(dims), subsys[0]));
732 subsys.erase(std::begin(subsys));
735 std::reverse(std::begin(result), std::end(result));
737 return std::make_tuple(result, prob, cA);
754 template <
typename Derived>
755 std::tuple<std::vector<idx>, double,
cmat>
756 measure_seq(
const Eigen::MatrixBase<Derived>& A, std::vector<idx> subsys,
758 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
772 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:75
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:222
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:68
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
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:440
Matrix is not a column vector exception.
Definition: exception.h:164
Quantum++ main namespace.
Definition: codes.h:35
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:242
idx get_num_subsys(idx sz, idx d)
Definition: util.h:366
Matrix is not square nor column vector exception.
Definition: exception.h:209
Invalid dimension(s) exception.
Definition: exception.h:269
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1733
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
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:680
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:83
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:205
Dimension(s) mismatch column vector size exception.
Definition: exception.h:316
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
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:1176
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:124
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
static RandomDevices & get_instance() noexcept(std::is_nothrow_constructible< RandomDevices >::value)
Definition: singleton.h:92
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.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:149
Object has zero size exception.
Definition: exception.h:134