Skip to content

Commit

Permalink
Refactor obs base and safer LQRaw call for measurement shots (#570)
Browse files Browse the repository at this point in the history 
* initial commit

* Auto update version

* add changelog

* fix typo

* Trigger MGPU CI

* quick fix namespace

* fix unit tests

* make format

* add more unit tests for better codecov

* update docstring

* fix typo

---------

Co-authored-by: Dev version update bot <github-actions[bot]@users.noreply.github.com>
  • Loading branch information
@multiphaseCFD @github-actions
multiphaseCFD and github-actions[bot] committed Dec 6, 2023
1 parent fc84399 commit a8f5f9f
Show file tree
Hide file tree
Showing 12 changed files with 247 additions and 206 deletions.
5 changes: 5 additions & 0 deletions .github/CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
# Release 0.34.0-dev

### New features since last release
* Shot-noise related methods now accommodate observable objects with arbitrary eigenvalues. Add a Kronecker product method for two diagonal matrices.
[(#570)](https://github.com/PennyLaneAI/pennylane-lightning/pull/570)

* Add shot-noise support for probs in the C++ layer. Probabilities are calculated from generated samples. All Lightning backends support this feature. Please note that target wires should be sorted in ascending manner.
[(#568)](https://github.com/PennyLaneAI/pennylane-lightning/pull/568)
Expand All @@ -21,6 +23,9 @@

### Improvements

* Refactor shot-noise related methods of MeasurementsBase class in the C++ layer and eigenvalues are not limited to `1` and `-1`. Add `getObs()` method to Observables class. Refactor `applyInPlaceShots` to allow users to get eigenvalues of Observables object. Deprecated `_preprocess_state` method in `MeasurementsBase` class for safer use of the `LightningQubitRaw` backend.
[(#570)](https://github.com/PennyLaneAI/pennylane-lightning/pull/570)

* Modify `setup.py` to use backend-specific build directory (`f"build_{backend}"`) to accelerate rebuilding backends in alternance.
[(#540)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/540)

Expand Down
2 changes: 1 addition & 1 deletion pennylane_lightning/core/_version.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,4 +16,4 @@
Version number (major.minor.patch[-label])
"""

__version__ = "0.34.0-dev12"
__version__ = "0.34.0-dev13"
185 changes: 76 additions & 109 deletions pennylane_lightning/core/src/measurements/MeasurementsBase.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,12 @@

#include "CPUMemoryModel.hpp"

#include "Util.hpp"

/// @cond DEV
namespace {
using namespace Pennylane::Observables;
using namespace Pennylane::Util;
} // namespace
/// @endcond

Expand Down Expand Up @@ -144,15 +147,11 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
"supported by shots");
} else if (obs.getObsName().find("Hamiltonian") != std::string::npos) {
auto coeffs = obs.getCoeffs();
auto obsTerms = obs.getObs();
for (size_t obs_term_idx = 0; obs_term_idx < coeffs.size();
obs_term_idx++) {
auto obs_samples = measure_with_samples(
obs, num_shots, shot_range, obs_term_idx);
PrecisionT result_per_term = std::accumulate(
obs_samples.begin(), obs_samples.end(), 0.0);

result +=
coeffs[obs_term_idx] * result_per_term / obs_samples.size();
result += coeffs[obs_term_idx] * expval(*obsTerms[obs_term_idx],
num_shots, shot_range);
}
} else {
auto obs_samples = measure_with_samples(obs, num_shots, shot_range);
Expand All @@ -170,60 +169,39 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
* @param num_shots Number of shots used to generate samples
* @param shot_range The range of samples to use. All samples are used
* by default.
* @param term_idx Index of a Hamiltonian term
*
* @return Expectation value with respect to the given observable.
*/
auto measure_with_samples(const Observable<StateVectorT> &obs,
const size_t &num_shots,
const std::vector<size_t> &shot_range,
size_t term_idx = 0) {
const std::vector<size_t> &shot_range) {
const size_t num_qubits = _statevector.getTotalNumQubits();
std::vector<size_t> obs_wires;
std::vector<size_t> identity_wires;
std::vector<std::vector<PrecisionT>> eigenValues;

auto sub_samples = _sample_state(obs, num_shots, shot_range, obs_wires,
identity_wires, term_idx);
auto sub_samples =
_sample_state(obs, num_shots, shot_range, obs_wires, eigenValues);

size_t num_samples = shot_range.empty() ? num_shots : shot_range.size();

std::vector<PrecisionT> obs_samples(num_samples, 0);

size_t num_identity_obs = identity_wires.size();
if (!identity_wires.empty()) {
size_t identity_obs_idx = 0;
for (size_t i = 0; i < obs_wires.size(); i++) {
if (identity_wires[identity_obs_idx] == obs_wires[i]) {
std::swap(obs_wires[identity_obs_idx], obs_wires[i]);
identity_obs_idx++;
}
}
std::vector<PrecisionT> eigenVals = eigenValues[0];

for (size_t i = 1; i < eigenValues.size(); i++) {
eigenVals = kronProd(eigenVals, eigenValues[i]);
}

for (size_t i = 0; i < num_samples; i++) {
std::vector<size_t> local_sample(obs_wires.size());
size_t idx = 0;
size_t wire_idx = 0;
for (auto &obs_wire : obs_wires) {
local_sample[idx] = sub_samples[i * num_qubits + obs_wire];
idx++;
idx += sub_samples[i * num_qubits + obs_wire]
<< (obs_wires.size() - 1 - wire_idx);
wire_idx++;
}

if (num_identity_obs != obs_wires.size()) {
// eigen values are `1` and `-1` for PauliX, PauliY, PauliZ,
// Hadamard gates the eigen value for a eigen vector |00001> is
// -1 since sum of the value at each bit position is odd
size_t bitSum = static_cast<size_t>(
std::accumulate(local_sample.begin() + num_identity_obs,
local_sample.end(), 0));
if ((bitSum & size_t{1}) == 1) {
obs_samples[i] = -1;
} else {
obs_samples[i] = 1;
}
} else {
// eigen value for Identity gate is `1`
obs_samples[i] = 1;
}
obs_samples[i] = eigenVals[idx];
}
return obs_samples;
}
Expand All @@ -238,35 +216,27 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
*/
auto var(const Observable<StateVectorT> &obs, const size_t &num_shots) {
if (obs.getObsName().find("Hamiltonian") == std::string::npos) {
// Branch for NamedObs and TensorProd observables
auto square_mean = expval(obs, num_shots, {});
PrecisionT result =
1 - square_mean *
square_mean; //`1` used here is because Eigenvalues for
// Paulis, Hadamard and Identity are {-1,
// 1}. Need to change based on eigen values
// when add Hermitian support.
auto obs_samples = measure_with_samples(obs, num_shots, {});
auto square_mean =
std::accumulate(obs_samples.begin(), obs_samples.end(), 0.0) /
obs_samples.size();
auto mean_square =
std::accumulate(obs_samples.begin(), obs_samples.end(), 0.0,
[](PrecisionT acc, PrecisionT element) {
return acc + element * element;
}) /
obs_samples.size();
PrecisionT result = mean_square - square_mean * square_mean;
return result;
}
// Branch for Hamiltonian observables
auto coeffs = obs.getCoeffs();
auto obs_terms = obs.getObs();

PrecisionT result{0.0};
size_t obs_term_idx = 0;
for (const auto &coeff : coeffs) {
std::vector<size_t> shot_range = {};
auto obs_samples =
measure_with_samples(obs, num_shots, shot_range, obs_term_idx);
PrecisionT expval_per_term =
std::accumulate(obs_samples.begin(), obs_samples.end(), 0.0);
auto term_mean = expval_per_term / obs_samples.size();

result +=
coeff * coeff *
(1 - term_mean *
term_mean); //`1` used here is because Eigenvalues for
// Paulis, Hadamard and Identity are {-1,
// 1}. Need to change based on eigen values
// when add Hermitian support.
result += coeff * coeff * var(*obs_terms[obs_term_idx], num_shots);
obs_term_idx++;
}
return result;
Expand All @@ -287,13 +257,30 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
obs.getObsName().find("Hamiltonian") != std::string::npos,
"Hamiltonian and Sparse Hamiltonian do not support samples().");
std::vector<size_t> obs_wires;
std::vector<size_t> identity_wires;
auto sv = _preprocess_state(obs, obs_wires, identity_wires);
Derived measure(sv);
if (num_shots != size_t{0}) {
return measure.probs(obs_wires, num_shots);
std::vector<std::vector<PrecisionT>> eigenvalues;
if constexpr (std::is_same_v<
typename StateVectorT::MemoryStorageT,
Pennylane::Util::MemoryStorageLocation::External>) {
std::vector<ComplexT> data_storage(
this->_statevector.getData(),
this->_statevector.getData() + this->_statevector.getLength());
StateVectorT sv(data_storage.data(), data_storage.size());
sv.updateData(data_storage.data(), data_storage.size());
obs.applyInPlaceShots(sv, eigenvalues, obs_wires);
Derived measure(sv);
if (num_shots > size_t{0}) {
return measure.probs(obs_wires, num_shots);
}
return measure.probs(obs_wires);
} else {
StateVectorT sv(_statevector);
obs.applyInPlaceShots(sv, eigenvalues, obs_wires);
Derived measure(sv);
if (num_shots > size_t{0}) {
return measure.probs(obs_wires, num_shots);
}
return measure.probs(obs_wires);
}
return measure.probs(obs_wires);
}

/**
Expand Down Expand Up @@ -366,11 +353,9 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
obs.getObsName().find("Hamiltonian") != std::string::npos,
"Hamiltonian and Sparse Hamiltonian do not support samples().");
std::vector<size_t> obs_wires;
std::vector<size_t> identity_wires;
std::vector<size_t> shot_range = {};
size_t term_idx = 0;

return measure_with_samples(obs, num_shots, shot_range, term_idx);
return measure_with_samples(obs, num_shots, shot_range);
}

/**
Expand Down Expand Up @@ -442,35 +427,6 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
}

private:
/**
* @brief Return preprocess state with a observable
*
* @param obs The observable to sample
* @param obs_wires Observable wires.
* @param identity_wires Wires of Identity gates
* @param term_idx Index of a Hamiltonian term. For other observables, its
* value is 0, which is set as default.
*
* @return a StateVectorT object
*/
auto _preprocess_state(const Observable<StateVectorT> &obs,
std::vector<size_t> &obs_wires,
std::vector<size_t> &identity_wires,
const size_t &term_idx = 0) {
if constexpr (std::is_same_v<
typename StateVectorT::MemoryStorageT,
Pennylane::Util::MemoryStorageLocation::External>) {
StateVectorT sv(_statevector.getData(), _statevector.getLength());
sv.updateData(_statevector.getData(), _statevector.getLength());
obs.applyInPlaceShots(sv, identity_wires, obs_wires, term_idx);
return sv;
} else {
StateVectorT sv(_statevector);
obs.applyInPlaceShots(sv, identity_wires, obs_wires, term_idx);
return sv;
}
}

/**
* @brief Return samples of a observable
*
Expand All @@ -479,22 +435,33 @@ template <class StateVectorT, class Derived> class MeasurementsBase {
* @param shot_range The range of samples to use. All samples are used by
* default.
* @param obs_wires Observable wires.
* @param identity_wires Wires of Identity gates
* @param term_idx Index of a Hamiltonian term. For other observables, its
* value is 0, which is set as default.
* @param eigenValues eigenvalues of the observable.
*
* @return std::vector<size_t> samples in std::vector
*/
auto _sample_state(const Observable<StateVectorT> &obs,
const size_t &num_shots,
const std::vector<size_t> &shot_range,
std::vector<size_t> &obs_wires,
std::vector<size_t> &identity_wires,
const size_t &term_idx = 0) {
std::vector<std::vector<PrecisionT>> &eigenValues) {
const size_t num_qubits = _statevector.getTotalNumQubits();
auto sv = _preprocess_state(obs, obs_wires, identity_wires, term_idx);
Derived measure(sv);
auto samples = measure.generate_samples(num_shots);
std::vector<size_t> samples;
if constexpr (std::is_same_v<
typename StateVectorT::MemoryStorageT,
Pennylane::Util::MemoryStorageLocation::External>) {
std::vector<ComplexT> data_storage(
this->_statevector.getData(),
this->_statevector.getData() + this->_statevector.getLength());
StateVectorT sv(data_storage.data(), data_storage.size());
obs.applyInPlaceShots(sv, eigenValues, obs_wires);
Derived measure(sv);
samples = measure.generate_samples(num_shots);
} else {
StateVectorT sv(_statevector);
obs.applyInPlaceShots(sv, eigenValues, obs_wires);
Derived measure(sv);
samples = measure.generate_samples(num_shots);
}

if (!shot_range.empty()) {
std::vector<size_t> sub_samples(shot_range.size() * num_qubits);
Expand Down
40 changes: 38 additions & 2 deletions pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,10 +217,12 @@ template <typename TypeList> void testProbabilitiesObs() {

// Defining the Statevector that will be measured.
auto statevector_data = createNonTrivialState<StateVectorT>();
auto sv_data = createNonTrivialState<StateVectorT>();

StateVectorT statevector(statevector_data.data(),
statevector_data.size());

StateVectorT sv(statevector_data.data(), statevector_data.size());
StateVectorT sv(sv_data.data(), sv_data.size());

DYNAMIC_SECTION("Test PauliX"
<< StateVectorToName<StateVectorT>::name) {
Expand Down Expand Up @@ -370,10 +372,11 @@ template <typename TypeList> void testProbabilitiesObsShots() {

// Defining the Statevector that will be measured.
auto statevector_data = createNonTrivialState<StateVectorT>();
auto sv_data = createNonTrivialState<StateVectorT>();
StateVectorT statevector(statevector_data.data(),
statevector_data.size());

StateVectorT sv(statevector_data.data(), statevector_data.size());
StateVectorT sv(sv_data.data(), sv_data.size());

DYNAMIC_SECTION("Test TensorProd XYZ"
<< StateVectorToName<StateVectorT>::name) {
Expand Down Expand Up @@ -1192,9 +1195,16 @@ template <typename TypeList> void testHamiltonianObsExpvalShot() {
std::vector<ComplexT> statevector_data{
{0.0, 0.0}, {0.0, 0.1}, {0.1, 0.1}, {0.1, 0.2},
{0.2, 0.2}, {0.3, 0.3}, {0.3, 0.4}, {0.4, 0.5}};

std::vector<ComplexT> sv_data{{0.0, 0.0}, {0.0, 0.1}, {0.1, 0.1},
{0.1, 0.2}, {0.2, 0.2}, {0.3, 0.3},
{0.3, 0.4}, {0.4, 0.5}};

StateVectorT statevector(statevector_data.data(),
statevector_data.size());

StateVectorT sv(sv_data.data(), sv_data.size());

// Initializing the measures class.
// This object attaches to the statevector allowing several measures.
Measurements<StateVectorT> Measurer(statevector);
Expand Down Expand Up @@ -1229,6 +1239,32 @@ template <typename TypeList> void testHamiltonianObsExpvalShot() {
REQUIRE(expected == Approx(res).margin(5e-2));
}

DYNAMIC_SECTION("TensorProd with shots_range "
<< StateVectorToName<StateVectorT>::name) {
auto X0 = std::make_shared<NamedObs<StateVectorT>>(
"PauliX", std::vector<size_t>{0});
auto Z1 = std::make_shared<NamedObs<StateVectorT>>(
"PauliZ", std::vector<size_t>{1});
auto obs0 = TensorProdObs<StateVectorT>::create({X0, Z1});

auto Y0 = std::make_shared<NamedObs<StateVectorT>>(
"PauliY", std::vector<size_t>{0});
auto H1 = std::make_shared<NamedObs<StateVectorT>>(
"Hadamard", std::vector<size_t>{1});
auto obs1 = TensorProdObs<StateVectorT>::create({Y0, H1});

auto obs =
Hamiltonian<StateVectorT>::create({0.1, 0.3}, {obs0, obs1});

Measurements<StateVectorT> Measurer_analytic(sv);
auto expected = Measurer_analytic.expval(*obs);

size_t num_shots = 2000;
auto res = Measurer.expval(*obs, num_shots, {});

REQUIRE(expected == Approx(res).margin(5e-2));
}

testHamiltonianObsExpvalShot<typename TypeList::Next>();
}
}
Expand Down
Loading

0 comments on commit a8f5f9f

Please sign in to comment.