Program Listing for File MeasurementsLQubit.hpp¶
↰ Return to documentation for file (pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp
)
// Copyright 2018-2023 Xanadu Quantum Technologies Inc.
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#pragma once
#include <algorithm>
#include <complex>
#include <cstdio>
#include <random>
#include <stack>
#include <type_traits>
#include <unordered_map>
#include <vector>
#include "LinearAlgebra.hpp"
#include "MeasurementsBase.hpp"
#include "NDPermuter.hpp"
#include "Observables.hpp"
#include "SparseLinAlg.hpp"
#include "StateVectorLQubitManaged.hpp"
#include "StateVectorLQubitRaw.hpp"
#include "TransitionKernels.hpp"
#include "Util.hpp" //transpose_state_tensor, sorting_indices
namespace {
using namespace Pennylane::Measures;
using namespace Pennylane::Observables;
using Pennylane::LightningQubit::StateVectorLQubitManaged;
using Pennylane::LightningQubit::Util::innerProdC;
namespace PUtil = Pennylane::Util;
} // namespace
namespace Pennylane::LightningQubit::Measures {
template <class StateVectorT>
class Measurements final
: public MeasurementsBase<StateVectorT, Measurements<StateVectorT>> {
private:
using PrecisionT = typename StateVectorT::PrecisionT;
using ComplexT = typename StateVectorT::ComplexT;
using BaseType = MeasurementsBase<StateVectorT, Measurements<StateVectorT>>;
public:
explicit Measurements(const StateVectorT &statevector)
: BaseType{statevector} {};
std::vector<PrecisionT> probs() {
const ComplexT *arr_data = this->_statevector.getData();
std::vector<PrecisionT> basis_probs(this->_statevector.getLength(), 0);
std::transform(
arr_data, arr_data + this->_statevector.getLength(),
basis_probs.begin(),
[](const ComplexT &z) -> PrecisionT { return std::norm(z); });
return basis_probs;
};
std::vector<PrecisionT>
probs(const std::vector<std::size_t> &wires,
[[maybe_unused]] const std::vector<std::size_t> &device_wires = {}) {
// Determining index that would sort the vector.
// This information is needed later.
const auto sorted_ind_wires = Pennylane::Util::sorting_indices(wires);
// Sorting wires.
std::vector<std::size_t> sorted_wires(wires.size());
for (size_t pos = 0; pos < wires.size(); pos++) {
sorted_wires[pos] = wires[sorted_ind_wires[pos]];
}
// Determining probabilities for the sorted wires.
const ComplexT *arr_data = this->_statevector.getData();
std::size_t num_qubits = this->_statevector.getNumQubits();
const std::vector<std::size_t> all_indices =
Gates::generateBitPatterns(sorted_wires, num_qubits);
const std::vector<std::size_t> all_offsets = Gates::generateBitPatterns(
Gates::getIndicesAfterExclusion(sorted_wires, num_qubits),
num_qubits);
std::vector<PrecisionT> probabilities(all_indices.size(), 0);
std::size_t ind_probs = 0;
for (auto index : all_indices) {
for (auto offset : all_offsets) {
probabilities[ind_probs] += std::norm(arr_data[index + offset]);
}
ind_probs++;
}
// Permute the data according to the required wire ordering
if (wires != sorted_wires) {
static constexpr std::size_t CACHE_SIZE = 8;
PUtil::Permuter<PUtil::DefaultPermuter<CACHE_SIZE>> p{};
std::vector<std::size_t> shape(wires.size(), 2);
std::vector<std::string> wire_labels_old(sorted_wires.size(), "");
std::vector<std::string> wire_labels_new(wires.size(), "");
std::transform(sorted_wires.begin(), sorted_wires.end(),
wire_labels_old.begin(), [](std::size_t index) {
return std::to_string(index);
});
std::transform(
wires.begin(), wires.end(), wire_labels_new.begin(),
[](std::size_t index) { return std::to_string(index); });
auto probs_sorted = probabilities;
p.Transpose(probabilities, shape, probs_sorted, wire_labels_old,
wire_labels_new);
return probs_sorted;
}
return probabilities;
}
std::vector<PrecisionT> probs(const Observable<StateVectorT> &obs,
std::size_t num_shots = 0) {
return BaseType::probs(obs, num_shots);
}
std::vector<PrecisionT> probs(size_t num_shots) {
return BaseType::probs(num_shots);
}
std::vector<PrecisionT> probs(const std::vector<std::size_t> &wires,
std::size_t num_shots) {
return BaseType::probs(wires, num_shots);
}
PrecisionT expval(const std::vector<ComplexT> &matrix,
const std::vector<std::size_t> &wires) {
// Copying the original state vector, for the application of the
// observable operator.
StateVectorLQubitManaged<PrecisionT> operator_statevector(
this->_statevector);
operator_statevector.applyMatrix(matrix, wires);
ComplexT expected_value = innerProdC(this->_statevector.getData(),
operator_statevector.getData(),
this->_statevector.getLength());
return std::real(expected_value);
};
PrecisionT expval(const std::string &operation,
const std::vector<std::size_t> &wires) {
// Copying the original state vector, for the application of the
// observable operator.
StateVectorLQubitManaged<PrecisionT> operator_statevector(
this->_statevector);
operator_statevector.applyOperation(operation, wires);
ComplexT expected_value = innerProdC(this->_statevector.getData(),
operator_statevector.getData(),
this->_statevector.getLength());
return std::real(expected_value);
};
template <class index_type>
PrecisionT expval(const index_type *row_map_ptr,
const index_type row_map_size,
const index_type *entries_ptr, const ComplexT *values_ptr,
const index_type numNNZ) {
PL_ABORT_IF(
(this->_statevector.getLength() != (size_t(row_map_size) - 1)),
"Statevector and Hamiltonian have incompatible sizes.");
auto operator_vector = Util::apply_Sparse_Matrix(
this->_statevector.getData(),
static_cast<index_type>(this->_statevector.getLength()),
row_map_ptr, row_map_size, entries_ptr, values_ptr, numNNZ);
ComplexT expected_value =
innerProdC(this->_statevector.getData(), operator_vector.data(),
this->_statevector.getLength());
return std::real(expected_value);
};
template <typename op_type>
std::vector<PrecisionT>
expval(const std::vector<op_type> &operations_list,
const std::vector<std::vector<std::size_t>> &wires_list) {
PL_ABORT_IF(
(operations_list.size() != wires_list.size()),
"The lengths of the list of operations and wires do not match.");
std::vector<PrecisionT> expected_value_list;
for (size_t index = 0; index < operations_list.size(); index++) {
expected_value_list.emplace_back(
expval(operations_list[index], wires_list[index]));
}
return expected_value_list;
}
auto expval(const Observable<StateVectorT> &ob) -> PrecisionT {
PrecisionT result{};
if constexpr (std::is_same_v<typename StateVectorT::MemoryStorageT,
MemoryStorageLocation::Internal>) {
StateVectorT sv(this->_statevector);
result = calculateObsExpval(sv, ob, this->_statevector);
} else if constexpr (std::is_same_v<
typename StateVectorT::MemoryStorageT,
MemoryStorageLocation::External>) {
std::vector<ComplexT> data_storage(
this->_statevector.getData(),
this->_statevector.getData() + this->_statevector.getLength());
StateVectorT sv(data_storage.data(), data_storage.size());
result = calculateObsExpval(sv, ob, this->_statevector);
} else {
PL_ABORT("Undefined memory storage location for StateVectorT.");
}
return result;
}
auto expval(const Observable<StateVectorT> &obs,
const std::size_t &num_shots,
const std::vector<std::size_t> &shot_range) -> PrecisionT {
return BaseType::expval(obs, num_shots, shot_range);
}
auto var(const Observable<StateVectorT> &obs, const std::size_t &num_shots)
-> PrecisionT {
return BaseType::var(obs, num_shots);
}
auto var(const Observable<StateVectorT> &ob) -> PrecisionT {
PrecisionT result{};
if constexpr (std::is_same_v<typename StateVectorT::MemoryStorageT,
MemoryStorageLocation::Internal>) {
StateVectorT sv(this->_statevector);
result = calculateObsVar(sv, ob, this->_statevector);
} else if constexpr (std::is_same_v<
typename StateVectorT::MemoryStorageT,
MemoryStorageLocation::External>) {
std::vector<ComplexT> data_storage(
this->_statevector.getData(),
this->_statevector.getData() + this->_statevector.getLength());
StateVectorT sv(data_storage.data(), data_storage.size());
result = calculateObsVar(sv, ob, this->_statevector);
} else {
PL_ABORT("Undefined memory storage location for StateVectorT.");
}
return result;
}
PrecisionT var(const std::string &operation,
const std::vector<std::size_t> &wires) {
// Copying the original state vector, for the application of the
// observable operator.
StateVectorLQubitManaged<PrecisionT> operator_statevector(
this->_statevector);
operator_statevector.applyOperation(operation, wires);
const std::complex<PrecisionT> *opsv_data =
operator_statevector.getData();
std::size_t orgsv_len = this->_statevector.getLength();
PrecisionT mean_square =
std::real(innerProdC(opsv_data, opsv_data, orgsv_len));
PrecisionT squared_mean = std::real(
innerProdC(this->_statevector.getData(), opsv_data, orgsv_len));
squared_mean = static_cast<PrecisionT>(std::pow(squared_mean, 2));
return (mean_square - squared_mean);
};
PrecisionT var(const std::vector<ComplexT> &matrix,
const std::vector<std::size_t> &wires) {
// Copying the original state vector, for the application of the
// observable operator.
StateVectorLQubitManaged<PrecisionT> operator_statevector(
this->_statevector);
operator_statevector.applyMatrix(matrix, wires);
const std::complex<PrecisionT> *opsv_data =
operator_statevector.getData();
std::size_t orgsv_len = this->_statevector.getLength();
PrecisionT mean_square =
std::real(innerProdC(opsv_data, opsv_data, orgsv_len));
PrecisionT squared_mean = std::real(
innerProdC(this->_statevector.getData(), opsv_data, orgsv_len));
squared_mean = static_cast<PrecisionT>(std::pow(squared_mean, 2));
return (mean_square - squared_mean);
};
template <typename op_type>
std::vector<PrecisionT>
var(const std::vector<op_type> &operations_list,
const std::vector<std::vector<std::size_t>> &wires_list) {
PL_ABORT_IF(
(operations_list.size() != wires_list.size()),
"The lengths of the list of operations and wires do not match.");
std::vector<PrecisionT> expected_value_list;
for (size_t index = 0; index < operations_list.size(); index++) {
expected_value_list.emplace_back(
var(operations_list[index], wires_list[index]));
}
return expected_value_list;
};
std::vector<std::size_t>
generate_samples_metropolis(const std::string &kernelname,
std::size_t num_burnin,
std::size_t num_samples) {
std::size_t num_qubits = this->_statevector.getNumQubits();
std::uniform_real_distribution<PrecisionT> distrib(0.0, 1.0);
std::vector<std::size_t> samples(num_samples * num_qubits, 0);
std::unordered_map<size_t, std::size_t> cache;
this->setRandomSeed();
TransitionKernelType transition_kernel = TransitionKernelType::Local;
if (kernelname == "NonZeroRandom") {
transition_kernel = TransitionKernelType::NonZeroRandom;
}
auto tk =
kernel_factory(transition_kernel, this->_statevector.getData(),
this->_statevector.getNumQubits());
std::size_t idx = 0;
// Burn In
for (size_t i = 0; i < num_burnin; i++) {
idx = metropolis_step(this->_statevector, tk, this->rng, distrib,
idx); // Burn-in.
}
// Sample
for (size_t i = 0; i < num_samples; i++) {
idx = metropolis_step(this->_statevector, tk, this->rng, distrib,
idx);
if (cache.contains(idx)) {
std::size_t cache_id = cache[idx];
auto it_temp = samples.begin() + cache_id * num_qubits;
std::copy(it_temp, it_temp + num_qubits,
samples.begin() + i * num_qubits);
}
// If not cached, compute
else {
for (size_t j = 0; j < num_qubits; j++) {
samples[i * num_qubits + (num_qubits - 1 - j)] =
(idx >> j) & 1U;
}
cache[idx] = i;
}
}
return samples;
}
template <class index_type>
PrecisionT var(const index_type *row_map_ptr, const index_type row_map_size,
const index_type *entries_ptr, const ComplexT *values_ptr,
const index_type numNNZ) {
PL_ABORT_IF(
(this->_statevector.getLength() != (size_t(row_map_size) - 1)),
"Statevector and Hamiltonian have incompatible sizes.");
auto operator_vector = Util::apply_Sparse_Matrix(
this->_statevector.getData(),
static_cast<index_type>(this->_statevector.getLength()),
row_map_ptr, row_map_size, entries_ptr, values_ptr, numNNZ);
const PrecisionT mean_square =
std::real(innerProdC(operator_vector.data(), operator_vector.data(),
operator_vector.size()));
const auto squared_mean = static_cast<PrecisionT>(
std::pow(std::real(innerProdC(operator_vector.data(),
this->_statevector.getData(),
operator_vector.size())),
2));
return (mean_square - squared_mean);
};
std::vector<std::size_t> generate_samples(size_t num_samples) {
const std::size_t num_qubits = this->_statevector.getNumQubits();
auto &&probabilities = probs();
std::vector<std::size_t> samples(num_samples * num_qubits, 0);
std::uniform_real_distribution<PrecisionT> distribution(0.0, 1.0);
std::unordered_map<size_t, std::size_t> cache;
this->setRandomSeed();
const std::size_t N = probabilities.size();
std::vector<double> bucket(N);
std::vector<std::size_t> bucket_partner(N);
std::stack<std::size_t> overfull_bucket_ids;
std::stack<std::size_t> underfull_bucket_ids;
for (size_t i = 0; i < N; i++) {
bucket[i] = N * probabilities[i];
bucket_partner[i] = i;
if (bucket[i] > 1.0) {
overfull_bucket_ids.push(i);
}
if (bucket[i] < 1.0) {
underfull_bucket_ids.push(i);
}
}
// Run alias algorithm
while (!underfull_bucket_ids.empty() && !overfull_bucket_ids.empty()) {
// get an overfull bucket
std::size_t i = overfull_bucket_ids.top();
// get an underfull bucket
std::size_t j = underfull_bucket_ids.top();
underfull_bucket_ids.pop();
// underfull bucket is partned with an overfull bucket
bucket_partner[j] = i;
bucket[i] = bucket[i] + bucket[j] - 1;
// if overfull bucket is now underfull
// put in underfull stack
if (bucket[i] < 1) {
overfull_bucket_ids.pop();
underfull_bucket_ids.push(i);
}
// if overfull bucket is full -> remove
else if (bucket[i] == 1.0) {
overfull_bucket_ids.pop();
}
}
// Pick samples
for (size_t i = 0; i < num_samples; i++) {
PrecisionT pct = distribution(this->rng) * N;
auto idx = static_cast<std::size_t>(pct);
if (pct - idx > bucket[idx]) {
idx = bucket_partner[idx];
}
// If cached, retrieve sample from cache
if (cache.contains(idx)) {
std::size_t cache_id = cache[idx];
auto it_temp = samples.begin() + cache_id * num_qubits;
std::copy(it_temp, it_temp + num_qubits,
samples.begin() + i * num_qubits);
}
// If not cached, compute
else {
for (size_t j = 0; j < num_qubits; j++) {
samples[i * num_qubits + (num_qubits - 1 - j)] =
(idx >> j) & 1U;
}
cache[idx] = i;
}
}
return samples;
}
private:
auto inline calculateObsExpval(StateVectorT &bra,
const Observable<StateVectorT> &obs,
const StateVectorT &ket) -> PrecisionT {
obs.applyInPlace(bra);
return std::real(
innerProdC(bra.getData(), ket.getData(), ket.getLength()));
}
auto inline calculateObsVar(StateVectorT &bra,
const Observable<StateVectorT> &obs,
const StateVectorT &ket) -> PrecisionT {
obs.applyInPlace(bra);
PrecisionT mean_square = std::real(
innerProdC(bra.getData(), bra.getData(), bra.getLength()));
auto squared_mean = static_cast<PrecisionT>(
std::pow(std::real(innerProdC(bra.getData(), ket.getData(),
ket.getLength())),
2));
return (mean_square - squared_mean);
}
std::size_t
metropolis_step(const StateVectorT &sv,
const std::unique_ptr<TransitionKernel<PrecisionT>> &tk,
std::mt19937 &gen,
std::uniform_real_distribution<PrecisionT> &distrib,
std::size_t init_idx) {
auto init_plog = std::log(
(sv.getData()[init_idx] * std::conj(sv.getData()[init_idx]))
.real());
auto init_qratio = tk->operator()(init_idx);
// transition kernel outputs these two
auto &trans_idx = init_qratio.first;
auto &trans_qratio = init_qratio.second;
auto trans_plog = std::log(
(sv.getData()[trans_idx] * std::conj(sv.getData()[trans_idx]))
.real());
auto alph = std::min<PrecisionT>(
1., trans_qratio * std::exp(trans_plog - init_plog));
auto ran = distrib(gen);
if (ran < alph) {
return trans_idx;
}
return init_idx;
}
}; // class Measurements
} // namespace Pennylane::LightningQubit::Measures
api/program_listing_file_pennylane_lightning_core_src_simulators_lightning_qubit_measurements_MeasurementsLQubit.hpp
Download Python script
Download Notebook
View on GitHub