Quantum++  v0.8.2
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 
46 measure(const Eigen::MatrixBase<Derived>& A, const std::vector <cmat>& Ks)
47 {
49 
50  // EXCEPTION CHECKS
51  // check zero-size
53  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
54 
55  // check the Kraus operators
56  if (Ks.size() == 0)
57  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
58  if (!internal::_check_square_mat(Ks[0]))
59  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
60  if (Ks[0].rows() != rA.rows())
61  throw Exception("qpp::measure()",
63  for (auto&& it : Ks)
64  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
65  throw Exception("qpp::measure()", Exception::Type::DIMS_NOT_EQUAL);
66  // END EXCEPTION CHECKS
67 
68  // probabilities
69  std::vector <double> prob(Ks.size());
70  // resulting states
71  std::vector <cmat> outstates(Ks.size());
72 
73  //************ density matrix ************//
74  if (internal::_check_square_mat(rA)) // square matrix
75  {
76  for (idx i = 0; i < Ks.size(); ++i)
77  {
78  outstates[i] = cmat::Zero(rA.rows(), rA.rows());
79  cmat tmp = Ks[i] * rA * adjoint(Ks[i]); // un-normalized;
80  prob[i] = std::abs(trace(tmp)); // probability
81  if (prob[i] > eps)
82  outstates[i] = tmp / prob[i]; // normalized
83  }
84  }
85  //************ ket ************//
86  else if (internal::_check_cvector(rA)) // column vector
87  {
88  for (idx i = 0; i < Ks.size(); ++i)
89  {
90  outstates[i] = ket::Zero(rA.rows());
91  ket tmp = Ks[i] * rA; // un-normalized;
92  // probability
93  prob[i] = std::pow(norm(tmp), 2);
94  if (prob[i] > eps)
95  outstates[i] = tmp / std::sqrt(prob[i]); // normalized
96  }
97  }
98  else
99  throw Exception("qpp::measure()",
101 
102  // sample from the probability distribution
103  std::discrete_distribution<idx> dd(std::begin(prob),
104  std::end(prob));
105  idx result = dd(RandomDevices::get_instance()._rng);
106 
107  return std::make_tuple(result, prob, outstates);
108 }
109 
110 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
111 // http://stackoverflow.com
112 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
122 template<typename Derived>
123 std::tuple<idx, std::vector < double>, std::vector<cmat>>
124 
125 measure(const Eigen::MatrixBase<Derived>& A,
126  const std::initializer_list<cmat>& Ks)
127 {
128  return measure(A, std::vector<cmat>(Ks));
129 }
130 
141 template<typename Derived>
142 std::tuple<idx, std::vector < double>, std::vector<cmat>>
143 
144 measure(const Eigen::MatrixBase<Derived>& A, const cmat& U)
145 {
146  const dyn_mat<typename Derived::Scalar>& rA = A;
147 
148  // EXCEPTION CHECKS
149  // check zero-size
151  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
152 
153  // check the unitary basis matrix U
155  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
157  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
158  if (U.rows() != rA.rows())
159  throw Exception("qpp::measure()",
161  // END EXCEPTION CHECKS
162 
163  std::vector <cmat> Ks(U.rows());
164  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
165  Ks[i] = U.col(i) * adjoint(U.col(i));
166 
167  return measure(rA, Ks);
168 }
169 
170 // partial measurements
188 template<typename Derived>
189 std::tuple<idx, std::vector < double>, std::vector<cmat>>
190 
191 measure(const Eigen::MatrixBase<Derived>& A,
192  const std::vector <cmat>& Ks,
193  const std::vector <idx>& subsys,
194  const std::vector <idx>& dims)
195 {
196  const cmat& rA = A;
197 
198  // EXCEPTION CHECKS
199 
200  // check zero-size
202  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
203 
204  // check that dimension is valid
205  if (!internal::_check_dims(dims))
206  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
207 
208  // check that dims match rho matrix
209  if (!internal::_check_dims_match_mat(dims, rA))
210  throw Exception("qpp::measure()",
212 
213  // check subsys is valid w.r.t. dims
214  if (!internal::_check_subsys_match_dims(subsys, dims))
215  throw Exception("qpp::measure()",
217 
218  std::vector <idx> subsys_dims(subsys.size());
219  for (idx i = 0; i < subsys.size(); ++i)
220  subsys_dims[i] = dims[subsys[i]];
221 
222  idx D = prod(std::begin(dims), std::end(dims));
223  idx Dsubsys = prod(std::begin(subsys_dims), std::end(subsys_dims));
224  idx Dbar = D / Dsubsys;
225 
226  // check the Kraus operators
227  if (Ks.size() == 0)
228  throw Exception("qpp::measure()",
230  if (!internal::_check_square_mat(Ks[0]))
231  throw Exception("qpp::measure()",
233  if (Dsubsys != static_cast<idx>(Ks[0].rows()))
234  throw Exception("qpp::measure()",
236  for (auto&& it : Ks)
237  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
238  throw Exception("qpp::measure()",
240  // END EXCEPTION CHECKS
241 
242  // probabilities
243  std::vector <double> prob(Ks.size());
244  // resulting states
245  std::vector <cmat> outstates(Ks.size());
246 
247  //************ density matrix ************//
248  if (internal::_check_square_mat(rA)) // square matrix
249  {
250  for (idx i = 0; i < Ks.size(); ++i)
251  {
252  outstates[i] = cmat::Zero(Dbar, Dbar);
253  cmat tmp = apply(rA, Ks[i], subsys, dims);
254  tmp = ptrace(tmp, subsys, dims);
255  prob[i] = std::abs(trace(tmp)); // probability
256  if (prob[i] > eps)
257  {
258  // normalized output state
259  // corresponding to measurement result i
260  outstates[i] = tmp / prob[i];
261  }
262  }
263  }
264  //************ ket ************//
265  else if (internal::_check_cvector(rA)) // column vector
266  {
267  for (idx i = 0; i < Ks.size(); ++i)
268  {
269  outstates[i] = cmat::Zero(Dbar, Dbar);
270  ket tmp = apply(rA, Ks[i], subsys, dims);
271  prob[i] = std::pow(norm(tmp), 2);
272  if (prob[i] > eps)
273  {
274  // normalized output state
275  // corresponding to measurement result i
276  tmp /= std::sqrt(prob[i]);
277  outstates[i] = ptrace(tmp, subsys, dims);
278  }
279  }
280  }
281  else
282  throw Exception("qpp::measure()",
284 
285  // sample from the probability distribution
286  std::discrete_distribution<idx> dd(std::begin(prob),
287  std::end(prob));
288  idx result = dd(RandomDevices::get_instance()._rng);
289 
290  return std::make_tuple(result, prob, outstates);
291 }
292 
293 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
294 // http://stackoverflow.com
295 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
313 template<typename Derived>
314 std::tuple<idx, std::vector < double>, std::vector<cmat>>
315 
316 measure(const Eigen::MatrixBase<Derived>& A,
317  const std::initializer_list<cmat>& Ks,
318  const std::vector <idx>& subsys,
319  const std::vector <idx>& dims)
320 {
321  return measure(A, std::vector<cmat>(Ks), subsys, dims);
322 }
323 
341 template<typename Derived>
342 std::tuple<idx, std::vector < double>, std::vector<cmat>>
343 
344 measure(const Eigen::MatrixBase<Derived>& A,
345  const std::vector <cmat>& Ks,
346  const std::vector <idx>& subsys,
347  const idx d = 2)
348 {
349  const cmat& rA = A;
350 
351  // check zero size
353  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
354 
355  idx n =
356  static_cast<idx>(std::llround(std::log2(rA.rows()) /
357  std::log2(d)));
358  std::vector <idx> dims(n, d); // local dimensions vector
359 
360  return measure(rA, Ks, subsys, dims);
361 }
362 
363 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
364 // http://stackoverflow.com
365 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
383 template<typename Derived>
384 std::tuple<idx, std::vector < double>, std::vector<cmat>>
385 
386 measure(const Eigen::MatrixBase<Derived>& A,
387  const std::initializer_list<cmat>& Ks,
388  const std::vector <idx>& subsys,
389  const idx d = 2)
390 {
391  return measure(A, std::vector<cmat>(Ks), subsys, d);
392 }
393 
411 template<typename Derived>
412 std::tuple<idx, std::vector < double>, std::vector<cmat>>
413 
414 measure(const Eigen::MatrixBase<Derived>& A,
415  const cmat& U,
416  const std::vector <idx>& subsys,
417  const std::vector <idx>& dims)
418 {
419  const cmat& rA = A;
420 
421  // EXCEPTION CHECKS
422 
423  // check zero-size
425  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
426 
427  // check that dimension is valid
428  if (!internal::_check_dims(dims))
429  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
430 
431  // check that dims match rho matrix
432  if (!internal::_check_dims_match_mat(dims, rA))
433  throw Exception("qpp::measure()",
435 
436  // check subsys is valid w.r.t. dims
437  if (!internal::_check_subsys_match_dims(subsys, dims))
438  throw Exception("qpp::measure()",
440 
441  std::vector <idx> subsys_dims(subsys.size());
442  for (idx i = 0; i < subsys.size(); ++i)
443  subsys_dims[i] = dims[subsys[i]];
444 
445  idx Dsubsys = prod(std::begin(subsys_dims), std::end(subsys_dims));
446 
447  // check the unitary basis matrix U
449  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
451  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
452  if (Dsubsys != static_cast<idx>(U.rows()))
453  throw Exception("qpp::measure()",
455  // END EXCEPTION CHECKS
456 
457  std::vector <cmat> Ks(U.rows());
458  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
459  Ks[i] = U.col(i) * adjoint(U.col(i));
460 
461  return measure(rA, Ks, subsys, dims);
462 }
463 
481 template<typename Derived>
482 std::tuple<idx, std::vector < double>, std::vector<cmat>>
483 
484 measure(const Eigen::MatrixBase<Derived>& A,
485  const cmat& U,
486  const std::vector <idx>& subsys,
487  const idx d = 2)
488 {
489  const cmat& rA = A;
490 
491  // check zero size
493  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
494 
495  idx n =
496  static_cast<idx>(std::llround(std::log2(rA.rows()) /
497  std::log2(d)));
498  std::vector <idx> dims(n, d); // local dimensions vector
499 
500  return measure(rA, U, subsys, dims);
501 }
502 
517 template<typename Derived>
518 std::tuple<std::vector < idx>, double, cmat>
519 
520 measure_seq(const Eigen::MatrixBase<Derived>& A,
521  std::vector <idx> subsys,
522  std::vector <idx> dims)
523 {
525 
526  // EXCEPTION CHECKS
527 
528  // check zero-size
530  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
531 
532  // check that dimension is valid
533  if (!internal::_check_dims(dims))
534  throw Exception("qpp::measure_seq()", Exception::Type::DIMS_INVALID);
535 
536  // check that dims match rho matrix
537  if (!internal::_check_dims_match_mat(dims, cA))
538  throw Exception("qpp::measure_seq()",
540 
541  // check subsys is valid w.r.t. dims
542  if (!internal::_check_subsys_match_dims(subsys, dims))
543  throw Exception("qpp::measure_seq()",
545 
546  // END EXCEPTION CHECKS
547 
548  std::vector <idx> result;
549  double prob = 1;
550 
551  // sort subsys in decreasing order,
552  // the order of measurements does not matter
553  std::sort(std::begin(subsys), std::end(subsys), std::greater<idx>{});
554 
555  //************ density matrix or column vector ************//
557  {
558  while (subsys.size() > 0)
559  {
560  auto tmp = measure(cA, Gates::get_instance().Id(dims[subsys[0]]),
561  {subsys[0]}, dims);
562  result.push_back(std::get<0>(tmp));
563  prob *= std::get<1>(tmp)[std::get<0>(tmp)];
564  cA = std::get<2>(tmp)[std::get<0>(tmp)];
565 
566  // remove the subsystem
567  dims.erase(std::next(std::begin(dims), subsys[0]));
568  subsys.erase(std::begin(subsys));
569  }
570  // order result in increasing order with respect to subsys
571  std::reverse(std::begin(result), std::end(result));
572  }
573  else
574  throw Exception("qpp::measure_seq()",
576 
577  return std::make_tuple(result, prob, cA);
578 }
579 
594 template<typename Derived>
595 std::tuple<std::vector < idx>, double, cmat>
596 
597 measure_seq(const Eigen::MatrixBase<Derived>& A,
598  std::vector <idx> subsys, idx d = 2)
599 {
600  const cmat& rA = A;
601 
602  // check zero size
604  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
605 
606  idx n =
607  static_cast<idx>(std::llround(std::log2(rA.rows()) /
608  std::log2(d)));
609  std::vector <idx> dims(n, d); // local dimensions vector
610 
611  return measure_seq(rA, subsys, dims);
612 }
613 
614 } /* namespace qpp */
615 
616 #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:190
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:142
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:73
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:83
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:520
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: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:66
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:126