Quantum++  v0.6
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>> measure(
45  const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks)
46 {
48 
49  // EXCEPTION CHECKS
50  // check zero-size
52  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
53 
54  // check the Kraus operators
55  if (Ks.size() == 0)
56  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
57  if (!internal::_check_square_mat(Ks[0]))
58  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
59  if (Ks[0].rows() != rA.rows())
60  throw Exception("qpp::measure()",
62  for (auto&& it : Ks)
63  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
64  throw Exception("qpp::measure()", Exception::Type::DIMS_NOT_EQUAL);
65  // END EXCEPTION CHECKS
66 
67  // probabilities
68  std::vector<double> prob(Ks.size());
69  // resulting states
70  std::vector<cmat> outstates(Ks.size());
71 
72  //************ density matrix ************//
73  if (internal::_check_square_mat(rA)) // square matrix
74  {
75  for (idx i = 0; i < Ks.size(); ++i)
76  {
77  outstates[i] = cmat::Zero(rA.rows(), rA.rows());
78  cmat tmp = Ks[i] * rA * adjoint(Ks[i]); // un-normalized;
79  prob[i] = std::abs(trace(tmp)); // probability
80  if (prob[i] > eps)
81  outstates[i] = tmp / prob[i]; // normalized
82  }
83  }
84  //************ ket ************//
85  else if (internal::_check_cvector(rA)) // column vector
86  {
87  for (idx i = 0; i < Ks.size(); ++i)
88  {
89  outstates[i] = ket::Zero(rA.rows());
90  ket tmp = Ks[i] * rA; // un-normalized;
91  // probability
92  prob[i] = std::pow(norm(tmp), 2);
93  if (prob[i] > eps)
94  outstates[i] = tmp / std::sqrt(prob[i]); // normalized
95  }
96  }
97  else
98  throw Exception("qpp::measure()",
100 
101  // sample from the probability distribution
102  std::discrete_distribution<idx> dd(std::begin(prob),
103  std::end(prob));
104  idx result = dd(RandomDevices::get_instance()._rng);
105 
106  return std::make_tuple(result, prob, outstates);
107 }
108 
109 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
110 // http://stackoverflow.com
111 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
121 template<typename Derived>
122 std::tuple<idx, std::vector<double>, std::vector<cmat>> measure(
123  const Eigen::MatrixBase<Derived>& A,
124  const std::initializer_list<cmat>& Ks)
125 {
126  return measure(A, std::vector<cmat>(Ks));
127 }
128 
139 template<typename Derived>
140 std::tuple<idx, std::vector<double>, std::vector<cmat>> measure(
141  const Eigen::MatrixBase<Derived>& A, const cmat& U)
142 {
143  const dyn_mat<typename Derived::Scalar>& rA = A;
144 
145  // EXCEPTION CHECKS
146  // check zero-size
148  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
149 
150  // check the unitary basis matrix U
152  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
154  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
155  if (U.rows() != rA.rows())
156  throw Exception("qpp::measure()",
158  // END EXCEPTION CHECKS
159 
160  std::vector<cmat> Ks(U.rows());
161  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
162  Ks[i] = U.col(i) * adjoint(U.col(i));
163 
164  return measure(rA, Ks);
165 }
166 
167 // partial measurements
185 template<typename Derived>
186 std::tuple<idx, std::vector<double>, std::vector<cmat>>
188  const Eigen::MatrixBase<Derived>& A,
189  const std::vector<cmat>& Ks,
190  const std::vector<idx>& subsys,
191  const std::vector<idx>& dims)
192 {
193  const cmat& rA = A;
194 
195  // EXCEPTION CHECKS
196 
197  // check zero-size
199  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
200 
201  // check that dimension is valid
202  if (!internal::_check_dims(dims))
203  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
204 
205  // check that dims match rho matrix
206  if (!internal::_check_dims_match_mat(dims, rA))
207  throw Exception("qpp::measure()",
209 
210  // check subsys is valid w.r.t. dims
211  if (!internal::_check_subsys_match_dims(subsys, dims))
212  throw Exception("qpp::measure()",
214 
215  std::vector<idx> subsys_dims(subsys.size());
216  for (idx i = 0; i < subsys.size(); ++i)
217  subsys_dims[i] = dims[subsys[i]];
218 
219  idx D = prod(dims.begin(), dims.end());
220  idx Dsubsys = prod(subsys_dims.begin(), subsys_dims.end());
221  idx Dbar = D / Dsubsys;
222 
223  // check the Kraus operators
224  if (Ks.size() == 0)
225  throw Exception("qpp::measure()",
227  if (!internal::_check_square_mat(Ks[0]))
228  throw Exception("qpp::measure()",
230  if (Dsubsys != static_cast<idx>(Ks[0].rows()))
231  throw Exception("qpp::measure()",
233  for (auto&& it : Ks)
234  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
235  throw Exception("qpp::measure()",
237  // END EXCEPTION CHECKS
238 
239  // probabilities
240  std::vector<double> prob(Ks.size());
241  // resulting states
242  std::vector<cmat> outstates(Ks.size());
243 
244  //************ density matrix ************//
245  if (internal::_check_square_mat(rA)) // square matrix
246  {
247  for (idx i = 0; i < Ks.size(); ++i)
248  {
249  outstates[i] = cmat::Zero(Dbar, Dbar);
250  cmat tmp = apply(rA, Ks[i], subsys, dims);
251  tmp = ptrace(tmp, subsys, dims);
252  prob[i] = std::abs(trace(tmp)); // probability
253  if (prob[i] > eps)
254  {
255  // normalized output state
256  // corresponding to measurement result i
257  outstates[i] = tmp / prob[i];
258  }
259  }
260  }
261  //************ ket ************//
262  else if (internal::_check_cvector(rA)) // column vector
263  {
264  for (idx i = 0; i < Ks.size(); ++i)
265  {
266  outstates[i] = cmat::Zero(Dbar, Dbar);
267  ket tmp = apply(rA, Ks[i], subsys, dims);
268  prob[i] = std::pow(norm(tmp), 2);
269  if (prob[i] > eps)
270  {
271  // normalized output state
272  // corresponding to measurement result i
273  tmp /= std::sqrt(prob[i]);
274  outstates[i] = ptrace(tmp, subsys, dims);
275  }
276  }
277  }
278  else
279  throw Exception("qpp::measure()",
281 
282  // sample from the probability distribution
283  std::discrete_distribution<idx> dd(std::begin(prob),
284  std::end(prob));
285  idx result = dd(RandomDevices::get_instance()._rng);
286 
287  return std::make_tuple(result, prob, outstates);
288 }
289 
290 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
291 // http://stackoverflow.com
292 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
310 template<typename Derived>
311 std::tuple<idx, std::vector<double>, std::vector<cmat>>
313  const Eigen::MatrixBase<Derived>& A,
314  const std::initializer_list<cmat>& Ks,
315  const std::vector<idx>& subsys,
316  const std::vector<idx>& dims)
317 {
318  return measure(A, std::vector<cmat>(Ks), subsys, dims);
319 }
320 
338 template<typename Derived>
339 std::tuple<idx, std::vector<double>, std::vector<cmat>>
341  const Eigen::MatrixBase<Derived>& A,
342  const std::vector<cmat>& Ks,
343  const std::vector<idx>& subsys,
344  const idx d = 2)
345 {
346  const cmat& rA = A;
347 
348  // check zero size
350  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
351 
352  idx n =
353  static_cast<idx>(std::llround(std::log2(rA.rows()) /
354  std::log2(d)));
355  std::vector<idx> dims(n, d); // local dimensions vector
356 
357  return measure(rA, Ks, subsys, dims);
358 }
359 
360 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
361 // http://stackoverflow.com
362 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
380 template<typename Derived>
381 std::tuple<idx, std::vector<double>, std::vector<cmat>>
383  const Eigen::MatrixBase<Derived>& A,
384  const std::initializer_list<cmat>& Ks,
385  const std::vector<idx>& subsys,
386  const idx d = 2)
387 {
388  return measure(A, std::vector<cmat>(Ks), subsys, d);
389 }
390 
408 template<typename Derived>
409 std::tuple<idx, std::vector<double>, std::vector<cmat>>
411  const Eigen::MatrixBase<Derived>& A,
412  const cmat& U,
413  const std::vector<idx>& subsys,
414  const std::vector<idx>& dims)
415 {
416  const cmat& rA = A;
417 
418  // EXCEPTION CHECKS
419 
420  // check zero-size
422  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
423 
424  // check that dimension is valid
425  if (!internal::_check_dims(dims))
426  throw Exception("qpp::measure()", Exception::Type::DIMS_INVALID);
427 
428  // check that dims match rho matrix
429  if (!internal::_check_dims_match_mat(dims, rA))
430  throw Exception("qpp::measure()",
432 
433  // check subsys is valid w.r.t. dims
434  if (!internal::_check_subsys_match_dims(subsys, dims))
435  throw Exception("qpp::measure()",
437 
438  std::vector<idx> subsys_dims(subsys.size());
439  for (idx i = 0; i < subsys.size(); ++i)
440  subsys_dims[i] = dims[subsys[i]];
441 
442  idx Dsubsys = prod(subsys_dims.begin(), subsys_dims.end());
443 
444  // check the unitary basis matrix U
446  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
448  throw Exception("qpp::measure()", Exception::Type::MATRIX_NOT_SQUARE);
449  if (Dsubsys != static_cast<idx>(U.rows()))
450  throw Exception("qpp::measure()",
452  // END EXCEPTION CHECKS
453 
454  std::vector<cmat> Ks(U.rows());
455  for (idx i = 0; i < static_cast<idx>(U.rows()); i++)
456  Ks[i] = U.col(i) * adjoint(U.col(i));
457 
458  return measure(rA, Ks, subsys, dims);
459 }
460 
478 template<typename Derived>
479 std::tuple<idx, std::vector<double>, std::vector<cmat>>
481  const Eigen::MatrixBase<Derived>& A,
482  const cmat& U,
483  const std::vector<idx>& subsys,
484  const idx d = 2)
485 {
486  const cmat& rA = A;
487 
488  // check zero size
490  throw Exception("qpp::measure()", Exception::Type::ZERO_SIZE);
491 
492  idx n =
493  static_cast<idx>(std::llround(std::log2(rA.rows()) /
494  std::log2(d)));
495  std::vector<idx> dims(n, d); // local dimensions vector
496 
497  return measure(rA, U, subsys, dims);
498 }
499 
514 template<typename Derived>
515 std::tuple<std::vector<idx>, double, cmat>
517  const Eigen::MatrixBase<Derived>& A,
518  std::vector<idx> subsys,
519  std::vector<idx> dims)
520 {
522 
523  // EXCEPTION CHECKS
524 
525  // check zero-size
527  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
528 
529  // check that dimension is valid
530  if (!internal::_check_dims(dims))
531  throw Exception("qpp::measure_seq()", Exception::Type::DIMS_INVALID);
532 
533  // check that dims match rho matrix
534  if (!internal::_check_dims_match_mat(dims, cA))
535  throw Exception("qpp::measure_seq()",
537 
538  // check subsys is valid w.r.t. dims
539  if (!internal::_check_subsys_match_dims(subsys, dims))
540  throw Exception("qpp::measure_seq()",
542 
543  // END EXCEPTION CHECKS
544 
545  std::vector<idx> result;
546  double prob = 1;
547 
548  // sort subsys in decreasing order,
549  // the order of measurements does not matter
550  std::sort(std::begin(subsys), std::end(subsys), std::greater<idx>{});
551 
552  //************ density matrix or column vector ************//
554  {
555  while (subsys.size() > 0)
556  {
557  auto tmp = measure(cA, Gates::get_instance().Id(dims[subsys[0]]),
558  {subsys[0]}, dims);
559  result.push_back(std::get<0>(tmp));
560  prob *= std::get<1>(tmp)[std::get<0>(tmp)];
561  cA = std::get<2>(tmp)[std::get<0>(tmp)];
562 
563  // remove the subsystem
564  dims.erase(std::next(dims.begin(), subsys[0]));
565  subsys.erase(subsys.begin());
566  }
567  // order result in increasing order with respect to subsys
568  std::reverse(std::begin(result), std::end(result));
569  }
570  else
571  throw Exception("qpp::measure_seq()",
573 
574  return std::make_tuple(result, prob, cA);
575 }
576 
591 template<typename Derived>
592 std::tuple<std::vector<idx>, double, cmat>
594  const Eigen::MatrixBase<Derived>& A,
595  std::vector<idx> subsys, idx d = 2)
596 {
597  const cmat& rA = A;
598 
599  // check zero size
601  throw Exception("qpp::measure_seq()", Exception::Type::ZERO_SIZE);
602 
603  idx n =
604  static_cast<idx>(std::llround(std::log2(rA.rows()) /
605  std::log2(d)));
606  std::vector<idx> dims(n, d); // local dimensions vector
607 
608  return measure_seq(rA, subsys, dims);
609 }
610 
611 } /* namespace qpp */
612 
613 #endif /* INSTRUMENTS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:110
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:189
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
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:73
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:46
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: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:516
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:44
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:119
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:56
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:125