Quantum++  v0.8
C++11 quantum computing library
instruments.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef INSTRUMENTS_H_
28 #define INSTRUMENTS_H_
29 
30 // measurements
31 namespace qpp
32 {
33 // full measurements
43 template<typename Derived>
44 std::tuple<idx, std::vector < double>, std::vector<cmat>>
45 
47  const Eigen::MatrixBase<Derived>& A, const std::vector <cmat>& Ks)
48 {
50 
51  // EXCEPTION CHECKS
52  // check zero-size
54  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
55 
56  // check the Kraus operators
57  if (Ks.size() == 0)
58  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
59  if (!internal::_check_square_mat(Ks[0]))
60  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
61  if (Ks[0].rows() != rA.rows())
62  throw Exception("qpp::measure()",
64  for (auto&& it : Ks)
65  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
66  throw Exception("qpp::measure()", Exception::Type::DIMS_NOT_EQUAL);
67  // END EXCEPTION CHECKS
68 
69  // probabilities
70  std::vector <double> prob(Ks.size());
71  // resulting states
72  std::vector <cmat> outstates(Ks.size());
73 
74  //************ density matrix ************//
75  if (internal::_check_square_mat(rA)) // square matrix
76  {
77  for (idx i = 0; i < Ks.size(); ++i)
78  {
79  outstates[i] = cmat::Zero(rA.rows(), rA.rows());
80  cmat tmp = Ks[i] * rA * adjoint(Ks[i]); // un-normalized;
81  prob[i] = std::abs(trace(tmp)); // probability
82  if (prob[i] > eps)
83  outstates[i] = tmp / prob[i]; // normalized
84  }
85  }
86  //************ ket ************//
87  else if (internal::_check_cvector(rA)) // column vector
88  {
89  for (idx i = 0; i < Ks.size(); ++i)
90  {
91  outstates[i] = ket::Zero(rA.rows());
92  ket tmp = Ks[i] * rA; // un-normalized;
93  // probability
94  prob[i] = std::pow(norm(tmp), 2);
95  if (prob[i] > eps)
96  outstates[i] = tmp / std::sqrt(prob[i]); // normalized
97  }
98  }
99  else
100  throw Exception("qpp::measure()",
102 
103  // sample from the probability distribution
104  std::discrete_distribution<idx> dd(std::begin(prob),
105  std::end(prob));
106  idx result = dd(RandomDevices::get_instance()._rng);
107 
108  return std::make_tuple(result, prob, outstates);
109 }
110 
111 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
112 // http://stackoverflow.com
113 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
123 template<typename Derived>
124 std::tuple<idx, std::vector < double>, std::vector<cmat>>
125 
127  const Eigen::MatrixBase<Derived>& A,
128  const std::initializer_list<cmat>& Ks)
129 {
130  return measure(A, std::vector<cmat>(Ks));
131 }
132 
143 template<typename Derived>
144 std::tuple<idx, std::vector < double>, std::vector<cmat>>
145 
147  const Eigen::MatrixBase<Derived>& A, const cmat& U)
148 {
149  const dyn_mat<typename Derived::Scalar>& rA = A;
150 
151  // EXCEPTION CHECKS
152  // check zero-size
154  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
155 
156  // check the unitary basis matrix U
158  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
160  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
161  if (U.rows() != rA.rows())
162  throw Exception("qpp::measure()",
164  // END EXCEPTION CHECKS
165 
166  std::vector <cmat> Ks(U.rows());
167  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
168  Ks[i] = U.col(i) * adjoint(U.col(i));
169 
170  return measure(rA, Ks);
171 }
172 
173 // partial measurements
191 template<typename Derived>
192 std::tuple<idx, std::vector < double>, std::vector<cmat>>
193 
195  const Eigen::MatrixBase<Derived>& A,
196  const std::vector <cmat>& Ks,
197  const std::vector <idx>& subsys,
198  const std::vector <idx>& dims)
199 {
200  const cmat& rA = A;
201 
202  // EXCEPTION CHECKS
203 
204  // check zero-size
206  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
207 
208  // check that dimension is valid
209  if (!internal::_check_dims(dims))
210  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
211 
212  // check that dims match rho matrix
213  if (!internal::_check_dims_match_mat(dims, rA))
214  throw Exception("qpp::measure()",
216 
217  // check subsys is valid w.r.t. dims
218  if (!internal::_check_subsys_match_dims(subsys, dims))
219  throw Exception("qpp::measure()",
221 
222  std::vector <idx> subsys_dims(subsys.size());
223  for (idx i = 0; i < subsys.size(); ++i)
224  subsys_dims[i] = dims[subsys[i]];
225 
226  idx D = prod(dims.begin(), dims.end());
227  idx Dsubsys = prod(subsys_dims.begin(), subsys_dims.end());
228  idx Dbar = D / Dsubsys;
229 
230  // check the Kraus operators
231  if (Ks.size() == 0)
232  throw Exception("qpp::measure()",
234  if (!internal::_check_square_mat(Ks[0]))
235  throw Exception("qpp::measure()",
237  if (Dsubsys != static_cast<idx>(Ks[0].rows()))
238  throw Exception("qpp::measure()",
240  for (auto&& it : Ks)
241  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
242  throw Exception("qpp::measure()",
244  // END EXCEPTION CHECKS
245 
246  // probabilities
247  std::vector <double> prob(Ks.size());
248  // resulting states
249  std::vector <cmat> outstates(Ks.size());
250 
251  //************ density matrix ************//
252  if (internal::_check_square_mat(rA)) // square matrix
253  {
254  for (idx i = 0; i < Ks.size(); ++i)
255  {
256  outstates[i] = cmat::Zero(Dbar, Dbar);
257  cmat tmp = apply(rA, Ks[i], subsys, dims);
258  tmp = ptrace(tmp, subsys, dims);
259  prob[i] = std::abs(trace(tmp)); // probability
260  if (prob[i] > eps)
261  {
262  // normalized output state
263  // corresponding to measurement result i
264  outstates[i] = tmp / prob[i];
265  }
266  }
267  }
268  //************ ket ************//
269  else if (internal::_check_cvector(rA)) // column vector
270  {
271  for (idx i = 0; i < Ks.size(); ++i)
272  {
273  outstates[i] = cmat::Zero(Dbar, Dbar);
274  ket tmp = apply(rA, Ks[i], subsys, dims);
275  prob[i] = std::pow(norm(tmp), 2);
276  if (prob[i] > eps)
277  {
278  // normalized output state
279  // corresponding to measurement result i
280  tmp /= std::sqrt(prob[i]);
281  outstates[i] = ptrace(tmp, subsys, dims);
282  }
283  }
284  }
285  else
286  throw Exception("qpp::measure()",
288 
289  // sample from the probability distribution
290  std::discrete_distribution<idx> dd(std::begin(prob),
291  std::end(prob));
292  idx result = dd(RandomDevices::get_instance()._rng);
293 
294  return std::make_tuple(result, prob, outstates);
295 }
296 
297 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
298 // http://stackoverflow.com
299 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
317 template<typename Derived>
318 std::tuple<idx, std::vector < double>, std::vector<cmat>>
319 
321  const Eigen::MatrixBase<Derived>& A,
322  const std::initializer_list<cmat>& Ks,
323  const std::vector <idx>& subsys,
324  const std::vector <idx>& dims)
325 {
326  return measure(A, std::vector<cmat>(Ks), subsys, dims);
327 }
328 
346 template<typename Derived>
347 std::tuple<idx, std::vector < double>, std::vector<cmat>>
348 
350  const Eigen::MatrixBase<Derived>& A,
351  const std::vector <cmat>& Ks,
352  const std::vector <idx>& subsys,
353  const idx d = 2)
354 {
355  const cmat& rA = A;
356 
357  // check zero size
359  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
360 
361  idx n =
362  static_cast<idx>(std::llround(std::log2(rA.rows()) /
363  std::log2(d)));
364  std::vector <idx> dims(n, d); // local dimensions vector
365 
366  return measure(rA, Ks, subsys, dims);
367 }
368 
369 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
370 // http://stackoverflow.com
371 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
389 template<typename Derived>
390 std::tuple<idx, std::vector < double>, std::vector<cmat>>
391 
393  const Eigen::MatrixBase<Derived>& A,
394  const std::initializer_list<cmat>& Ks,
395  const std::vector <idx>& subsys,
396  const idx d = 2)
397 {
398  return measure(A, std::vector<cmat>(Ks), subsys, d);
399 }
400 
418 template<typename Derived>
419 std::tuple<idx, std::vector < double>, std::vector<cmat>>
420 
422  const Eigen::MatrixBase<Derived>& A,
423  const cmat& U,
424  const std::vector <idx>& subsys,
425  const std::vector <idx>& dims)
426 {
427  const cmat& rA = A;
428 
429  // EXCEPTION CHECKS
430 
431  // check zero-size
433  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
434 
435  // check that dimension is valid
436  if (!internal::_check_dims(dims))
437  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
438 
439  // check that dims match rho matrix
440  if (!internal::_check_dims_match_mat(dims, rA))
441  throw Exception("qpp::measure()",
443 
444  // check subsys is valid w.r.t. dims
445  if (!internal::_check_subsys_match_dims(subsys, dims))
446  throw Exception("qpp::measure()",
448 
449  std::vector <idx> subsys_dims(subsys.size());
450  for (idx i = 0; i < subsys.size(); ++i)
451  subsys_dims[i] = dims[subsys[i]];
452 
453  idx Dsubsys = prod(subsys_dims.begin(), subsys_dims.end());
454 
455  // check the unitary basis matrix U
457  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
459  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
460  if (Dsubsys != static_cast<idx>(U.rows()))
461  throw Exception("qpp::measure()",
463  // END EXCEPTION CHECKS
464 
465  std::vector <cmat> Ks(U.rows());
466  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
467  Ks[i] = U.col(i) * adjoint(U.col(i));
468 
469  return measure(rA, Ks, subsys, dims);
470 }
471 
489 template<typename Derived>
490 std::tuple<idx, std::vector < double>, std::vector<cmat>>
491 
493  const Eigen::MatrixBase<Derived>& A,
494  const cmat& U,
495  const std::vector <idx>& subsys,
496  const idx d = 2)
497 {
498  const cmat& rA = A;
499 
500  // check zero size
502  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
503 
504  idx n =
505  static_cast<idx>(std::llround(std::log2(rA.rows()) /
506  std::log2(d)));
507  std::vector <idx> dims(n, d); // local dimensions vector
508 
509  return measure(rA, U, subsys, dims);
510 }
511 
526 template<typename Derived>
527 std::tuple<std::vector < idx>, double, cmat>
528 
530  const Eigen::MatrixBase<Derived>& A,
531  std::vector <idx> subsys,
532  std::vector <idx> dims)
533 {
535 
536  // EXCEPTION CHECKS
537 
538  // check zero-size
540  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
541 
542  // check that dimension is valid
543  if (!internal::_check_dims(dims))
544  throw Exception("qpp::measure_seq()", Exception::Type::DIMS_INVALID);
545 
546  // check that dims match rho matrix
547  if (!internal::_check_dims_match_mat(dims, cA))
548  throw Exception("qpp::measure_seq()",
550 
551  // check subsys is valid w.r.t. dims
552  if (!internal::_check_subsys_match_dims(subsys, dims))
553  throw Exception("qpp::measure_seq()",
555 
556  // END EXCEPTION CHECKS
557 
558  std::vector <idx> result;
559  double prob = 1;
560 
561  // sort subsys in decreasing order,
562  // the order of measurements does not matter
563  std::sort(std::begin(subsys), std::end(subsys), std::greater<idx>{});
564 
565  //************ density matrix or column vector ************//
567  {
568  while (subsys.size() > 0)
569  {
570  auto tmp = measure(cA, Gates::get_instance().Id(dims[subsys[0]]),
571  {subsys[0]}, dims);
572  result.push_back(std::get<0>(tmp));
573  prob *= std::get<1>(tmp)[std::get<0>(tmp)];
574  cA = std::get<2>(tmp)[std::get<0>(tmp)];
575 
576  // remove the subsystem
577  dims.erase(std::next(dims.begin(), subsys[0]));
578  subsys.erase(subsys.begin());
579  }
580  // order result in increasing order with respect to subsys
581  std::reverse(std::begin(result), std::end(result));
582  }
583  else
584  throw Exception("qpp::measure_seq()",
586 
587  return std::make_tuple(result, prob, cA);
588 }
589 
604 template<typename Derived>
605 std::tuple<std::vector < idx>, double, cmat>
606 
608  const Eigen::MatrixBase<Derived>& A,
609  std::vector <idx> subsys, idx d = 2)
610 {
611  const cmat& rA = A;
612 
613  // check zero size
615  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
616 
617  idx n =
618  static_cast<idx>(std::llround(std::log2(rA.rows()) /
619  std::log2(d)));
620  std::vector <idx> dims(n, d); // local dimensions vector
621 
622  return measure_seq(rA, subsys, dims);
623 }
624 
625 } /* namespace qpp */
626 
627 #endif /* INSTRUMENTS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
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:135
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:213
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:75
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:84
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:529
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:57
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:412
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:231
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
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:1079
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:121
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:67
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:46
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:119