XACC
instruments.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2020 Vlad Gheorghiu.
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef INSTRUMENTS_H_
33 #define INSTRUMENTS_H_
34 
35 namespace qpp {
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) {
50  const dyn_col_vect<typename Derived::Scalar>& rphi = phi.derived();
51  const dyn_col_vect<typename Derived::Scalar>& rpsi = psi.derived();
52 
53  // EXCEPTION CHECKS
54 
55  // check zero-size
56  if (!internal::check_nonzero_size(rphi))
57  throw exception::ZeroSize("qpp::ip()");
58 
59  // check zero-size
60  if (!internal::check_nonzero_size(rpsi))
61  throw exception::ZeroSize("qpp::ip()");
62 
63  // check column vector
64  if (!internal::check_cvector(rphi))
65  throw exception::MatrixNotCvector("qpp::ip()");
66 
67  // check column vector
68  if (!internal::check_cvector(rpsi))
69  throw exception::MatrixNotCvector("qpp::ip()");
70 
71  // check that dims is a valid dimension vector
72  if (!internal::check_dims(dims))
73  throw exception::DimsInvalid("qpp::ip()");
74 
75  // check that subsys are valid w.r.t. dims
76  if (!internal::check_subsys_match_dims(subsys, dims))
77  throw exception::SubsysMismatchDims("qpp::ip()");
78 
79  // check that dims match psi column vector
80  if (!internal::check_dims_match_cvect(dims, rpsi))
81  throw exception::DimsMismatchCvector("qpp::ip()");
82 
83  // check that subsys match pho column vector
84  std::vector<idx> subsys_dims(subsys.size());
85  for (idx i = 0; i < subsys.size(); ++i)
86  subsys_dims[i] = dims[subsys[i]];
87  if (!internal::check_dims_match_cvect(subsys_dims, rphi))
88  throw exception::DimsMismatchCvector("qpp::ip()");
89  // END EXCEPTION CHECKS
90 
91  idx Dsubsys = prod(std::begin(subsys_dims), std::end(subsys_dims));
92 
93  idx D = static_cast<idx>(rpsi.rows());
94  idx Dsubsys_bar = D / Dsubsys;
95 
96  idx n = dims.size();
97  idx n_subsys = subsys.size();
98  idx n_subsys_bar = n - n_subsys;
99 
100  idx Cdims[maxn];
101  idx Csubsys[maxn];
102  idx Cdimssubsys[maxn];
103  idx Csubsys_bar[maxn];
104  idx Cdimssubsys_bar[maxn];
105 
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));
109 
110  for (idx i = 0; i < n; ++i) {
111  Cdims[i] = dims[i];
112  }
113  for (idx i = 0; i < n_subsys; ++i) {
114  Csubsys[i] = subsys[i];
115  Cdimssubsys[i] = dims[subsys[i]];
116  }
117  for (idx i = 0; i < n_subsys_bar; ++i) {
118  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
119  }
120 
121  auto worker = [&](idx b) noexcept->typename Derived::Scalar {
122  idx Cmidxrow[maxn];
123  idx Cmidxrowsubsys[maxn];
124  idx Cmidxcolsubsys_bar[maxn];
125 
126  /* get the col multi-indexes of the complement */
127  internal::n2multiidx(b, n_subsys_bar, Cdimssubsys_bar,
128  Cmidxcolsubsys_bar);
129  /* write it in the global row multi-index */
130  for (idx k = 0; k < n_subsys_bar; ++k) {
131  Cmidxrow[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
132  }
133 
134  typename Derived::Scalar result = 0;
135  for (idx a = 0; a < Dsubsys; ++a) {
136  /* get the row multi-indexes of the subsys */
137  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxrowsubsys);
138  /* write it in the global row multi-index */
139  for (idx k = 0; k < n_subsys; ++k) {
140  Cmidxrow[Csubsys[k]] = Cmidxrowsubsys[k];
141  }
142  // compute the row index
143  idx i = internal::multiidx2n(Cmidxrow, n, Cdims);
144 
145  result += std::conj(rphi(a)) * rpsi(i);
146  }
147 
148  return result;
149  }; /* end worker */
150 
151  dyn_col_vect<typename Derived::Scalar> result(Dsubsys_bar);
152 #ifdef WITH_OPENMP_
153 #pragma omp parallel for
154 #endif // WITH_OPENMP_
155  for (idx m = 0; m < Dsubsys_bar; ++m)
156  result(m) = worker(m);
157 
158  return result;
159 }
160 
171 template <typename Derived>
172 dyn_col_vect<typename Derived::Scalar>
173 ip(const Eigen::MatrixBase<Derived>& phi, const Eigen::MatrixBase<Derived>& psi,
174  const std::vector<idx>& subsys, idx d = 2) {
175  const dyn_col_vect<typename Derived::Scalar>& rphi = phi.derived();
176  const dyn_col_vect<typename Derived::Scalar>& rpsi = psi.derived();
177 
178  // EXCEPTION CHECKS
179 
180  if (!internal::check_nonzero_size(rpsi))
181  throw exception::ZeroSize("qpp::ip()");
182 
183  // check valid dims
184  if (d < 2)
185  throw exception::DimsInvalid("qpp::ip()");
186  // END EXCEPTION CHECKS
187 
188  idx n = internal::get_num_subsys(static_cast<idx>(rpsi.rows()), d);
189  std::vector<idx> dims(n, d); // local dimensions vector
190 
191  return ip(phi, psi, subsys, dims);
192 }
193 
194 // full measurements
204 template <typename Derived>
205 std::tuple<idx, std::vector<double>, std::vector<cmat>>
206 measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
207  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
208 
209  // EXCEPTION CHECKS
210 
211  // check zero-size
212  if (!internal::check_nonzero_size(rA))
213  throw exception::ZeroSize("qpp::measure()");
214 
215  // check the Kraus operators
216  if (Ks.empty())
217  throw exception::ZeroSize("qpp::measure()");
218  if (!internal::check_square_mat(Ks[0]))
219  throw exception::MatrixNotSquare("qpp::measure()");
220  if (Ks[0].rows() != rA.rows())
221  throw exception::DimsMismatchMatrix("qpp::measure()");
222  for (auto&& elem : Ks)
223  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
224  throw exception::DimsNotEqual("qpp::measure()");
225  // END EXCEPTION CHECKS
226 
227  // probabilities
228  std::vector<double> prob(Ks.size());
229  // resulting states
230  std::vector<cmat> outstates(Ks.size());
231 
232  //************ density matrix ************//
233  if (internal::check_square_mat(rA)) // square matrix
234  {
235  for (idx i = 0; i < Ks.size(); ++i) {
236  outstates[i] = cmat::Zero(rA.rows(), rA.rows());
237  cmat tmp = Ks[i] * rA * adjoint(Ks[i]); // un-normalized;
238  prob[i] = std::abs(trace(tmp)); // probability
239  if (prob[i] > 0)
240  outstates[i] = tmp / prob[i]; // normalized
241  }
242  }
243  //************ ket ************//
244  else if (internal::check_cvector(rA)) // column vector
245  {
246  for (idx i = 0; i < Ks.size(); ++i) {
247  outstates[i] = ket::Zero(rA.rows());
248  ket tmp = Ks[i] * rA; // un-normalized;
249  // probability
250  prob[i] = std::pow(norm(tmp), 2);
251  if (prob[i] > 0)
252  outstates[i] = tmp / std::sqrt(prob[i]); // normalized
253  }
254  } else
255  throw exception::MatrixNotSquareNorCvector("qpp::measure()");
256 
257  // sample from the probability distribution
258  std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
259  auto& gen =
260 #ifdef NO_THREAD_LOCAL_
261  RandomDevices::get_instance().get_prng();
262 #else
263  RandomDevices::get_thread_local_instance().get_prng();
264 #endif
265  idx result = dd(gen);
266 
267  return std::make_tuple(result, prob, outstates);
268 }
269 
270 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
271 // http://stackoverflow.com
272 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
282 template <typename Derived>
283 std::tuple<idx, std::vector<double>, std::vector<cmat>>
284 measure(const Eigen::MatrixBase<Derived>& A,
285  const std::initializer_list<cmat>& Ks) {
286  return measure(A, std::vector<cmat>(Ks));
287 }
288 
301 template <typename Derived>
302 std::tuple<idx, std::vector<double>, std::vector<cmat>>
303 measure(const Eigen::MatrixBase<Derived>& A, const cmat& U) {
304  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
305 
306  // EXCEPTION CHECKS
307 
308  // check zero-size
309  if (!internal::check_nonzero_size(rA))
310  throw exception::ZeroSize("qpp::measure()");
311 
312  // check the unitary basis matrix U
313  if (!internal::check_nonzero_size(U))
314  throw exception::ZeroSize("qpp::measure()");
315  if (!internal::check_square_mat(U))
316  throw exception::MatrixNotSquare("qpp::measure()");
317  if (U.rows() != rA.rows())
318  throw exception::DimsMismatchMatrix("qpp::measure()");
319  // END EXCEPTION CHECKS
320 
321  std::vector<cmat> Ks(U.rows());
322  for (idx i = 0; i < static_cast<idx>(U.rows()); ++i)
323  Ks[i] = U.col(i) * adjoint(U.col(i));
324 
325  return measure(rA, Ks);
326 }
327 
328 // partial measurements
346 template <typename Derived>
347 std::tuple<idx, std::vector<double>, std::vector<cmat>>
348 measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
349  const std::vector<idx>& target, const std::vector<idx>& dims,
350  bool destructive = true) {
351  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
352 
353  // EXCEPTION CHECKS
354 
355  // check zero-size
356  if (!internal::check_nonzero_size(rA))
357  throw exception::ZeroSize("qpp::measure()");
358 
359  // check that dimension is valid
360  if (!internal::check_dims(dims))
361  throw exception::DimsInvalid("qpp::measure()");
362 
363  // check that target is valid w.r.t. dims
364  if (!internal::check_subsys_match_dims(target, dims))
365  throw exception::SubsysMismatchDims("qpp::measure()");
366 
367  // check valid state and matching dimensions
368  if (internal::check_cvector(rA)) {
369  if (!internal::check_dims_match_cvect(dims, rA))
370  throw exception::DimsMismatchCvector("qpp::measure()");
371  } else if (internal::check_square_mat(rA)) {
372  if (!internal::check_dims_match_mat(dims, rA))
373  throw exception::DimsMismatchMatrix("qpp::measure()");
374  } else
375  throw exception::MatrixNotSquareNorCvector("qpp::measure()");
376 
377  std::vector<idx> subsys_dims(target.size());
378  for (idx i = 0; i < target.size(); ++i)
379  subsys_dims[i] = dims[target[i]];
380 
381  idx D = prod(std::begin(dims), std::end(dims));
382  idx Dsubsys = prod(std::begin(subsys_dims), std::end(subsys_dims));
383  idx Dsubsys_bar = D / Dsubsys;
384 
385  // check the Kraus operators
386  if (Ks.empty())
387  throw exception::ZeroSize("qpp::measure()");
388  if (!internal::check_square_mat(Ks[0]))
389  throw exception::MatrixNotSquare("qpp::measure()");
390  if (Dsubsys != static_cast<idx>(Ks[0].rows()))
391  throw exception::DimsMismatchMatrix("qpp::measure()");
392  for (auto&& elem : Ks)
393  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
394  throw exception::DimsNotEqual("qpp::measure()");
395  // END EXCEPTION CHECKS
396 
397  // probabilities
398  std::vector<double> prob(Ks.size());
399  // resulting states
400  std::vector<cmat> outstates;
401 
402  if (destructive)
403  outstates.resize(Ks.size(), cmat::Zero(Dsubsys_bar, Dsubsys_bar));
404  else
405  outstates.resize(Ks.size(), cmat::Zero(D, D));
406 
407  //************ ket ************//
408  if (internal::check_cvector(rA)) // column vector
409  {
410  for (idx i = 0; i < Ks.size(); ++i) {
411  ket tmp = apply(rA, Ks[i], target, dims);
412  prob[i] = std::pow(norm(tmp), 2);
413  if (prob[i] > 0) {
414  // normalized output state
415  // corresponding to measurement result i
416  tmp /= std::sqrt(prob[i]);
417  if (destructive)
418  outstates[i] = ptrace(tmp, target, dims);
419  else
420  outstates[i] = tmp;
421  }
422  }
423  }
424  //************ density matrix ************//
425  else // square matrix
426  {
427  for (idx i = 0; i < Ks.size(); ++i) {
428  cmat tmp = apply(rA, Ks[i], target, dims);
429  if (destructive)
430  tmp = ptrace(tmp, target, dims);
431  prob[i] = std::abs(trace(tmp)); // probability
432  if (prob[i] > 0) {
433  // normalized output state
434  // corresponding to measurement result i
435  outstates[i] = tmp / prob[i];
436  }
437  }
438  }
439  // sample from the probability distribution
440  std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
441  auto& gen =
442 #ifdef NO_THREAD_LOCAL_
443  RandomDevices::get_instance().get_prng();
444 #else
445  RandomDevices::get_thread_local_instance().get_prng();
446 #endif
447  idx result = dd(gen);
448 
449  return std::make_tuple(result, prob, outstates);
450 }
451 
452 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
453 // http://stackoverflow.com
454 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
472 template <typename Derived>
473 std::tuple<idx, std::vector<double>, std::vector<cmat>>
474 measure(const Eigen::MatrixBase<Derived>& A,
475  const std::initializer_list<cmat>& Ks, const std::vector<idx>& target,
476  const std::vector<idx>& dims, bool destructive = true) {
477  return measure(A, std::vector<cmat>(Ks), target, dims, destructive);
478 }
479 
497 template <typename Derived>
498 std::tuple<idx, std::vector<double>, std::vector<cmat>>
499 measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
500  const std::vector<idx>& target, idx d = 2, bool destructive = true) {
501  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
502 
503  // EXCEPTION CHECKS
504 
505  // check zero size
506  if (!internal::check_nonzero_size(rA))
507  throw exception::ZeroSize("qpp::measure()");
508 
509  // check valid dims
510  if (d < 2)
511  throw exception::DimsInvalid("qpp::measure()");
512  // END EXCEPTION CHECKS
513 
514  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
515  std::vector<idx> dims(n, d); // local dimensions vector
516 
517  return measure(rA, Ks, target, dims, destructive);
518 }
519 
520 // std::initializer_list overload, avoids ambiguity for 2-element lists, see
521 // http://stackoverflow.com
522 // /questions/26750039/ambiguity-when-using-initializer-list-as-parameter
540 template <typename Derived>
541 std::tuple<idx, std::vector<double>, std::vector<cmat>>
542 measure(const Eigen::MatrixBase<Derived>& A,
543  const std::initializer_list<cmat>& Ks, const std::vector<idx>& target,
544  idx d = 2, bool destructive = true) {
545  return measure(A, std::vector<cmat>(Ks), target, d, destructive);
546 }
547 
571 template <typename Derived>
572 std::tuple<idx, std::vector<double>, std::vector<cmat>>
573 measure(const Eigen::MatrixBase<Derived>& A, const cmat& V,
574  const std::vector<idx>& target, const std::vector<idx>& dims,
575  bool destructive = true) {
576  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
577 
578  // EXCEPTION CHECKS
579 
580  // check zero-size
581  if (!internal::check_nonzero_size(rA))
582  throw exception::ZeroSize("qpp::measure()");
583 
584  // check that dimension is valid
585  if (!internal::check_dims(dims))
586  throw exception::DimsInvalid("qpp::measure()");
587 
588  // check that target is valid w.r.t. dims
589  if (!internal::check_subsys_match_dims(target, dims))
590  throw exception::SubsysMismatchDims("qpp::measure()");
591 
592  // check valid state and matching dimensions
593  if (internal::check_cvector(rA)) {
594  if (!internal::check_dims_match_cvect(dims, rA))
595  throw exception::DimsMismatchCvector("qpp::measure()");
596  } else if (internal::check_square_mat(rA)) {
597  if (!internal::check_dims_match_mat(dims, rA))
598  throw exception::DimsMismatchMatrix("qpp::measure()");
599  } else
600  throw exception::MatrixNotSquareNorCvector("qpp::measure()");
601 
602  std::vector<idx> subsys_dims(target.size());
603  for (idx i = 0; i < target.size(); ++i)
604  subsys_dims[i] = dims[target[i]];
605 
606  idx Dsubsys = prod(std::begin(subsys_dims), std::end(subsys_dims));
607 
608  // check the matrix V
609  if (!internal::check_nonzero_size(V))
610  throw exception::ZeroSize("qpp::measure()");
611  if (Dsubsys != static_cast<idx>(V.rows()))
612  throw exception::DimsMismatchMatrix("qpp::measure()");
613  // END EXCEPTION CHECKS
614 
615  // number of basis elements or number of rank-1 projectors
616  idx M = static_cast<idx>(V.cols());
617 
618  //************ ket ************//
619  if (internal::check_cvector(rA)) {
620  const ket& rpsi = A.derived();
621  std::vector<double> prob(M); // probabilities
622  std::vector<cmat> outstates(M); // resulting states
623 
624 #ifdef WITH_OPENMP_
625 #pragma omp parallel for
626 #endif // WITH_OPENMP_
627  for (idx i = 0; i < M; ++i) {
628  if (destructive)
629  outstates[i] =
630  ip(static_cast<const ket&>(V.col(i)), rpsi, target, dims);
631  else
632  outstates[i] = apply(rpsi, prj(V.col(i)), target, dims);
633  }
634 
635  for (idx i = 0; i < M; ++i) {
636  double tmp = norm(outstates[i]);
637  prob[i] = tmp * tmp;
638  if (prob[i] > 0) {
639  // normalized output state
640  // corresponding to measurement result m
641  outstates[i] /= tmp;
642  }
643  }
644 
645  // sample from the probability distribution
646  std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
647  auto& gen =
648 #ifdef NO_THREAD_LOCAL_
649  RandomDevices::get_instance().get_prng();
650 #else
651  RandomDevices::get_thread_local_instance().get_prng();
652 #endif
653  idx result = dd(gen);
654 
655  return std::make_tuple(result, prob, outstates);
656  }
657  //************ density matrix ************//
658  else {
659  std::vector<cmat> Ks(M);
660  for (idx i = 0; i < M; ++i)
661  Ks[i] = V.col(i) * adjoint(V.col(i));
662 
663  return measure(rA, Ks, target, dims, destructive);
664  }
665 }
666 
690 template <typename Derived>
691 std::tuple<idx, std::vector<double>, std::vector<cmat>>
692 measure(const Eigen::MatrixBase<Derived>& A, const cmat& V,
693  const std::vector<idx>& target, idx d = 2, bool destructive = true) {
694  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
695 
696  // EXCEPTION CHECKS
697 
698  // check zero size
699  if (!internal::check_nonzero_size(rA))
700  throw exception::ZeroSize("qpp::measure()");
701 
702  // check valid dims
703  if (d < 2)
704  throw exception::DimsInvalid("qpp::measure()");
705  // END EXCEPTION CHECKS
706 
707  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
708  std::vector<idx> dims(n, d); // local dimensions vector
709 
710  return measure(rA, V, target, dims, destructive);
711 }
712 
730 template <typename Derived>
731 std::tuple<std::vector<idx>, double, cmat>
732 measure_seq(const Eigen::MatrixBase<Derived>& A, std::vector<idx> target,
733  std::vector<idx> dims, bool destructive = true) {
734  // typename std::remove_const<
735  // typename Eigen::MatrixBase<Derived>::EvalReturnType
736  // >::type rA = A.derived();
737 
738  dyn_mat<typename Derived::Scalar> rA = A.derived();
739 
740  // EXCEPTION CHECKS
741 
742  // check zero-size
743  if (!internal::check_nonzero_size(rA))
744  throw exception::ZeroSize("qpp::measure_seq()");
745 
746  // check that dimension is valid
747  if (!internal::check_dims(dims))
748  throw exception::DimsInvalid("qpp::measure_seq()");
749 
750  // check valid state and matching dimensions
751  if (internal::check_cvector(rA)) {
752  if (!internal::check_dims_match_cvect(dims, rA))
753  throw exception::DimsMismatchCvector("qpp::measure_seq()");
754  } else if (internal::check_square_mat(rA)) {
755  if (!internal::check_dims_match_mat(dims, rA))
756  throw exception::DimsMismatchMatrix("qpp::measure_seq()");
757  } else
758  throw exception::MatrixNotSquareNorCvector("qpp::measure_seq()");
759 
760  // check that target is valid w.r.t. dims
761  if (!internal::check_subsys_match_dims(target, dims))
762  throw exception::SubsysMismatchDims("qpp::measure_seq()");
763  // END EXCEPTION CHECKS
764 
765  std::vector<idx> result;
766  double prob = 1;
767 
768  // sort target in decreasing order,
769  // the order of measurements does not matter
770  std::sort(std::begin(target), std::end(target), std::greater<idx>{});
771 
772  //************ density matrix or column vector ************//
773  while (target.size() > 0) {
774  auto tmp = measure(rA, Gates::get_instance().Id(dims[target[0]]),
775  {target[0]}, dims, destructive);
776  result.emplace_back(std::get<0>(tmp));
777  prob *= std::get<1>(tmp)[std::get<0>(tmp)];
778  rA = std::get<2>(tmp)[std::get<0>(tmp)];
779 
780  if (destructive) {
781  // remove the subsystem
782  dims.erase(std::next(std::begin(dims), target[0]));
783  }
784  target.erase(std::begin(target));
785  }
786  // order result in increasing order with respect to target
787  std::reverse(std::begin(result), std::end(result));
788 
789  return std::make_tuple(result, prob, rA);
790 }
791 
809 template <typename Derived>
810 std::tuple<std::vector<idx>, double, cmat>
811 measure_seq(const Eigen::MatrixBase<Derived>& A, std::vector<idx> target,
812  idx d = 2, bool destructive = true) {
813  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
814 
815  // EXCEPTION CHECKS
816 
817  // check zero size
818  if (!internal::check_nonzero_size(rA))
819  throw exception::ZeroSize("qpp::measure_seq()");
820 
821  // check valid dims
822  if (d < 2)
823  throw exception::DimsInvalid("qpp::measure_seq()");
824  // END EXCEPTION CHECKS
825 
826  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
827  std::vector<idx> dims(n, d); // local dimensions vector
828 
829  return measure_seq(rA, target, dims, destructive);
830 }
831 
843 template <typename Derived>
844 dyn_mat<typename Derived::Scalar> reset(const Eigen::MatrixBase<Derived>& A,
845  std::vector<idx> target,
846  std::vector<idx> dims) {
847  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
848 
849  // EXCEPTION CHECKS
850 
851  // check zero-size
852  if (!internal::check_nonzero_size(rA))
853  throw exception::ZeroSize("qpp::reset()");
854 
855  // check that dimension is valid
856  if (!internal::check_dims(dims))
857  throw exception::DimsInvalid("qpp::reset()");
858 
859  // check valid state and matching dimensions
860  if (internal::check_cvector(rA)) {
861  if (!internal::check_dims_match_cvect(dims, rA))
862  throw exception::DimsMismatchCvector("qpp::reset()");
863  } else if (internal::check_square_mat(rA)) {
864  if (!internal::check_dims_match_mat(dims, rA))
865  throw exception::DimsMismatchMatrix("qpp::reset()");
866  } else
867  throw exception::MatrixNotSquareNorCvector("qpp::reset()");
868 
869  // check that target is valid w.r.t. dims
870  if (!internal::check_subsys_match_dims(target, dims))
871  throw exception::SubsysMismatchDims("qpp::reset()");
872  // END EXCEPTION CHECKS
873 
874  dyn_mat<typename Derived::Scalar> result;
875  std::vector<idx> resZ;
876 
877  std::tie(resZ, std::ignore, result) = measure_seq(rA, target, dims, false);
878  for (idx i = 0; i < target.size(); ++i) {
879  cmat correction =
880  powm(Gates::get_instance().Xd(dims[i]), dims[i] - resZ[i]);
881  result = apply(result, correction, {target[i]}, dims);
882  }
883 
884  return result;
885 }
886 
898 template <typename Derived>
899 dyn_mat<typename Derived::Scalar> reset(const Eigen::MatrixBase<Derived>& A,
900  std::vector<idx> target, idx d = 2) {
901  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
902 
903  // EXCEPTION CHECKS
904 
905  // check zero size
906  if (!internal::check_nonzero_size(rA))
907  throw exception::ZeroSize("qpp::reset()");
908 
909  // check valid dims
910  if (d < 2)
911  throw exception::DimsInvalid("qpp::reset()");
912  // END EXCEPTION CHECKS
913 
914  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
915  std::vector<idx> dims(n, d); // local dimensions vector
916 
917  return reset(rA, target, dims);
918 }
919 
930 template <typename Derived>
931 dyn_mat<typename Derived::Scalar> discard(const Eigen::MatrixBase<Derived>& A,
932  std::vector<idx> target,
933  std::vector<idx> dims) {
934  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
935 
936  // EXCEPTION CHECKS
937 
938  // check zero-size
939  if (!internal::check_nonzero_size(rA))
940  throw exception::ZeroSize("qpp::discard()");
941 
942  // check that dimension is valid
943  if (!internal::check_dims(dims))
944  throw exception::DimsInvalid("qpp::discard()");
945 
946  // check valid state and matching dimensions
947  if (internal::check_cvector(rA)) {
948  if (!internal::check_dims_match_cvect(dims, rA))
949  throw exception::DimsMismatchCvector("qpp::discard()");
950  } else if (internal::check_square_mat(rA)) {
951  if (!internal::check_dims_match_mat(dims, rA))
952  throw exception::DimsMismatchMatrix("qpp::discard()");
953  } else
954  throw exception::MatrixNotSquareNorCvector("qpp::discard()");
955 
956  // check that target is valid w.r.t. dims
957  if (!internal::check_subsys_match_dims(target, dims))
958  throw exception::SubsysMismatchDims("qpp::discard()");
959  // END EXCEPTION CHECKS
960 
961  dyn_mat<typename Derived::Scalar> result;
962  std::tie(std::ignore, std::ignore, result) = measure_seq(rA, target, dims);
963 
964  return result;
965 }
966 
977 template <typename Derived>
978 dyn_mat<typename Derived::Scalar> discard(const Eigen::MatrixBase<Derived>& A,
979  std::vector<idx> target, idx d = 2) {
980  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
981 
982  // EXCEPTION CHECKS
983 
984  // check zero size
985  if (!internal::check_nonzero_size(rA))
986  throw exception::ZeroSize("qpp::discard()");
987 
988  // check valid dims
989  if (d < 2)
990  throw exception::DimsInvalid("qpp::discard()");
991  // END EXCEPTION CHECKS
992 
993  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
994  std::vector<idx> dims(n, d); // local dimensions vector
995 
996  return discard(rA, target, dims);
997 }
998 
999 } /* namespace qpp */
1000 
1001 #endif /* INSTRUMENTS_H_ */
qpp
Quantum++ main namespace.
Definition: circuits.h:35
Generated by   doxygen 1.8.11