Skip to content

Commit

Permalink
assert not all probs are zero when measuring
Browse files Browse the repository at this point in the history
  • Loading branch information
@vsoftco
vsoftco committed Jan 10, 2021
1 parent dfb30c0 commit edab322
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 24 deletions.
1 change: 1 addition & 0 deletions include/classes/noise.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,7 @@ class NoiseBase {
"qpp::NoiseBase::compute_state_()");

// now do the actual noise generation
assert(probs_ != decltype(probs_)(probs_.size(), 0)); // not all zeros
std::discrete_distribution<idx> dd{std::begin(probs_),
std::end(probs_)};
auto& gen =
Expand Down
51 changes: 28 additions & 23 deletions include/instruments.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
// END EXCEPTION CHECKS

// probabilities
std::vector<double> prob(Ks.size());
std::vector<double> probs(Ks.size());
// resulting states
std::vector<cmat> outstates(Ks.size());

Expand All @@ -238,9 +238,10 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
for (idx i = 0; i < Ks.size(); ++i) {
outstates[i] = cmat::Zero(Dout, Dout);
cmat tmp = Ks[i] * rA * adjoint(Ks[i]); // un-normalized;
prob[i] = std::abs(trace(tmp)); // probability
if (prob[i] > 0)
outstates[i] = tmp / prob[i]; // normalized
probs[i] = std::abs(trace(tmp)); // probability
if (probs[i] > 0) {
outstates[i] = tmp / probs[i]; // normalized
}
}
}
//************ ket ************//
Expand All @@ -250,15 +251,17 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
outstates[i] = ket::Zero(Dout);
ket tmp = Ks[i] * rA; // un-normalized;
// probability
prob[i] = std::pow(norm(tmp), 2);
if (prob[i] > 0)
outstates[i] = tmp / std::sqrt(prob[i]); // normalized
probs[i] = std::pow(norm(tmp), 2);
if (probs[i] > 0) {
outstates[i] = tmp / std::sqrt(probs[i]); // normalized
}
}
} else
throw exception::MatrixNotSquareNorCvector("qpp::measure()");

// sample from the probability distribution
std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
assert(probs != decltype(probs)(probs.size(), 0)); // not all zeros
std::discrete_distribution<idx> dd(std::begin(probs), std::end(probs));
auto& gen =
#ifdef NO_THREAD_LOCAL_
RandomDevices::get_instance().get_prng();
Expand All @@ -267,7 +270,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
#endif
idx result = dd(gen);

return std::make_tuple(result, prob, outstates);
return std::make_tuple(result, probs, outstates);
}

// std::initializer_list overload, avoids ambiguity for 2-element lists, see
Expand Down Expand Up @@ -398,7 +401,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
// END EXCEPTION CHECKS

// probabilities
std::vector<double> prob(Ks.size());
std::vector<double> probs(Ks.size());
// resulting states
std::vector<cmat> outstates;

Expand All @@ -412,11 +415,11 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
{
for (idx i = 0; i < Ks.size(); ++i) {
ket tmp = apply(rA, Ks[i], target, dims);
prob[i] = std::pow(norm(tmp), 2);
if (prob[i] > 0) {
probs[i] = std::pow(norm(tmp), 2);
if (probs[i] > 0) {
// normalized output state
// corresponding to measurement result i
tmp /= std::sqrt(prob[i]);
tmp /= std::sqrt(probs[i]);
if (destructive)
outstates[i] = ptrace(tmp, target, dims);
else
Expand All @@ -431,16 +434,17 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
cmat tmp = apply(rA, Ks[i], target, dims);
if (destructive)
tmp = ptrace(tmp, target, dims);
prob[i] = std::abs(trace(tmp)); // probability
if (prob[i] > 0) {
probs[i] = std::abs(trace(tmp)); // probability
if (probs[i] > 0) {
// normalized output state
// corresponding to measurement result i
outstates[i] = tmp / prob[i];
outstates[i] = tmp / probs[i];
}
}
}
// sample from the probability distribution
std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
assert(probs != decltype(probs)(probs.size(), 0)); // not all zeros
std::discrete_distribution<idx> dd(std::begin(probs), std::end(probs));
auto& gen =
#ifdef NO_THREAD_LOCAL_
RandomDevices::get_instance().get_prng();
Expand All @@ -449,7 +453,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
#endif
idx result = dd(gen);

return std::make_tuple(result, prob, outstates);
return std::make_tuple(result, probs, outstates);
}

// std::initializer_list overload, avoids ambiguity for 2-element lists, see
Expand Down Expand Up @@ -621,7 +625,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const cmat& V,
//************ ket ************//
if (internal::check_cvector(rA)) {
const ket& rpsi = A.derived();
std::vector<double> prob(M); // probabilities
std::vector<double> probs(M); // probabilities
std::vector<cmat> outstates(M); // resulting states

#ifdef HAS_OPENMP
Expand All @@ -638,16 +642,17 @@ measure(const Eigen::MatrixBase<Derived>& A, const cmat& V,

for (idx i = 0; i < M; ++i) {
double tmp = norm(outstates[i]);
prob[i] = tmp * tmp;
if (prob[i] > 0) {
probs[i] = tmp * tmp;
if (probs[i] > 0) {
// normalized output state
// corresponding to measurement result m
outstates[i] /= tmp;
}
}

// sample from the probability distribution
std::discrete_distribution<idx> dd(std::begin(prob), std::end(prob));
assert(probs != decltype(probs)(probs.size(), 0)); // not all zeros
std::discrete_distribution<idx> dd(std::begin(probs), std::end(probs));
auto& gen =
#ifdef NO_THREAD_LOCAL_
RandomDevices::get_instance().get_prng();
Expand All @@ -656,7 +661,7 @@ measure(const Eigen::MatrixBase<Derived>& A, const cmat& V,
#endif
idx result = dd(gen);

return std::make_tuple(result, prob, outstates);
return std::make_tuple(result, probs, outstates);
}
//************ density matrix ************//
else {
Expand Down
2 changes: 1 addition & 1 deletion include/internal/util.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -399,7 +399,7 @@ T abs_chop(const T& x, double chop = qpp::chop) {
template <typename T,
typename std::enable_if<!(std::numeric_limits<T>::is_iec559 ||
is_complex<T>::value)>::type* = nullptr>
T abs_chop(const T& x, double QPP_UNUSED_ chop = qpp::chop) {
T abs_chop(const T& x, QPP_UNUSED_ double chop = qpp::chop) {
return x;
}

Expand Down

0 comments on commit edab322

Please sign in to comment.