diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index f5cfe8075b..d40aad0c43 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -2,6 +2,9 @@ ### New features since last release +* Add shots support for expectation value calculation for given observables (`NamedObs`, `TensorProd` and `Hamiltonian`) based on Pauli words, `Identity` and `Hadamard` in the C++ layer by adding `measure_with_samples` to the measurement interface. All Lightning backends support this support feature. +[(#556)](https://github.com/PennyLaneAI/pennylane-lightning/pull/556) + * `qml.QubitUnitary` operators can be included in a circuit differentiated with the adjoint method. Lightning handles circuits with arbitrary non-differentiable `qml.QubitUnitary` operators. 1,2-qubit `qml.QubitUnitary` operators with differentiable parameters can be differentiated using decomposition. [(#540)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/540) diff --git a/mpitests/test_expval.py b/mpitests/test_expval.py index ad76da1aa5..9db1e76b31 100644 --- a/mpitests/test_expval.py +++ b/mpitests/test_expval.py @@ -31,8 +31,6 @@ class TestExpval: def test_identity_expectation(self, theta, phi, tol): """Test that identity expectation value (i.e. the trace) is 1""" dev = qml.device(device_name, mpi=True, wires=3) - if device_name == "lightning.gpu" and dev.R_DTYPE == np.float32: - pytest.skip("Skipped FP32 tests for expval in lightning.gpu") O1 = qml.Identity(wires=[0]) O2 = qml.Identity(wires=[1]) @@ -49,9 +47,6 @@ def test_pauliz_expectation(self, theta, phi, tol): """Test that PauliZ expectation value is correct""" dev = qml.device(device_name, mpi=True, wires=3) - if device_name == "lightning.gpu" and dev.R_DTYPE == np.float32: - pytest.skip("Skipped FP32 tests for expval in lightning.gpu") - O1 = qml.PauliZ(wires=[0]) O2 = qml.PauliZ(wires=[1]) @@ -67,9 +62,6 @@ def test_paulix_expectation(self, theta, phi, tol): """Test that PauliX expectation value is correct""" dev = qml.device(device_name, mpi=True, wires=3) - if device_name == "lightning.gpu" and dev.R_DTYPE == np.float32: - pytest.skip("Skipped FP32 tests for expval in lightning.gpu") - O1 = qml.PauliX(wires=[0]) O2 = qml.PauliX(wires=[1]) @@ -89,9 +81,6 @@ def test_pauliy_expectation(self, theta, phi, tol): """Test that PauliY expectation value is correct""" dev = qml.device(device_name, mpi=True, wires=3) - if device_name == "lightning.gpu" and dev.R_DTYPE == np.float32: - pytest.skip("Skipped FP32 tests for expval in lightning.gpu") - O1 = qml.PauliY(wires=[0]) O2 = qml.PauliY(wires=[1]) @@ -130,8 +119,6 @@ def test_hermitian_expectation(self, n_wires, theta, phi, tol): n_qubits = 7 dev_def = qml.device("default.qubit", wires=n_qubits) dev = qml.device(device_name, mpi=True, wires=n_qubits) - if device_name == "lightning.gpu" and dev.R_DTYPE == np.float32: - pytest.skip("Skipped FP32 tests for expval in lightning.gpu") comm = MPI.COMM_WORLD m = 2**n_wires diff --git a/pennylane_lightning/core/_version.py b/pennylane_lightning/core/_version.py index 9daab1b916..226ada7853 100644 --- a/pennylane_lightning/core/_version.py +++ b/pennylane_lightning/core/_version.py @@ -16,4 +16,4 @@ Version number (major.minor.patch[-label]) """ -__version__ = "0.34.0-dev7" +__version__ = "0.34.0-dev8" diff --git a/pennylane_lightning/core/src/measurements/MeasurementsBase.hpp b/pennylane_lightning/core/src/measurements/MeasurementsBase.hpp index 9dbf6f9958..879377b3d7 100644 --- a/pennylane_lightning/core/src/measurements/MeasurementsBase.hpp +++ b/pennylane_lightning/core/src/measurements/MeasurementsBase.hpp @@ -17,10 +17,13 @@ */ #pragma once +#include #include #include "Observables.hpp" +#include "CPUMemoryModel.hpp" + /// @cond DEV namespace { using namespace Pennylane::Observables; @@ -111,6 +114,190 @@ template class MeasurementsBase { auto generate_samples(size_t num_samples) -> std::vector { return static_cast(this)->generate_samples(num_samples); }; -}; -} // namespace Pennylane::Measures \ No newline at end of file + /** + * @brief Calculate the expectation value for a general Observable. + * + * @param obs Observable. + * @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. + * + * @return Expectation value with respect to the given observable. + */ + auto expval(const Observable &obs, const size_t &num_shots, + const std::vector &shot_range = {}) -> PrecisionT { + PrecisionT result{0.0}; + + if (obs.getObsName().find("SparseHamiltonian") != std::string::npos) { + // SparseHamiltonian does not support samples in pennylane. + PL_ABORT("For SparseHamiltonian Observables, expval calculation is " + "not supported by shots"); + } else if (obs.getObsName().find("Hermitian") != std::string::npos) { + // TODO support. This support requires an additional method to solve + // eigenpair and unitary matrices, and the results of eigenpair and + // unitary matrices data need to be added to the Hermitian class and + // public methods are need to access eigen values. Note the + // assumption that eigen values are -1 and 1 in the + // `measurement_with_sample` method should be updated as well. + PL_ABORT("For Hermitian Observables, expval calculation is not " + "supported by shots"); + } else if (obs.getObsName().find("Hamiltonian") != std::string::npos) { + auto coeffs = obs.getCoeffs(); + 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(); + } + } else { + auto obs_samples = measure_with_samples(obs, num_shots, shot_range); + result = + std::accumulate(obs_samples.begin(), obs_samples.end(), 0.0); + result /= obs_samples.size(); + } + return result; + } + + /** + * @brief Calculate the expectation value for a general Observable. + * + * @param obs Observable. + * @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 &obs, + const size_t &num_shots, + const std::vector &shot_range, + size_t term_idx = 0) { + const size_t num_qubits = _statevector.getTotalNumQubits(); + std::vector obs_wires; + std::vector identity_wires; + + auto sub_samples = _sample_state(obs, num_shots, shot_range, obs_wires, + identity_wires, term_idx); + + size_t num_samples = shot_range.empty() ? num_shots : shot_range.size(); + + std::vector 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++; + } + } + } + + for (size_t i = 0; i < num_samples; i++) { + std::vector local_sample(obs_wires.size()); + size_t idx = 0; + for (auto &obs_wire : obs_wires) { + local_sample[idx] = sub_samples[i * num_qubits + obs_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( + 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; + } + } + return obs_samples; + } + + 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 &obs, + std::vector &obs_wires, + std::vector &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 + * + * @param obs The observable to sample + * @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 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 std::vector samples in std::vector + */ + auto _sample_state(const Observable &obs, + const size_t &num_shots, + const std::vector &shot_range, + std::vector &obs_wires, + std::vector &identity_wires, + const size_t &term_idx = 0) { + 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); + + if (!shot_range.empty()) { + std::vector sub_samples(shot_range.size() * num_qubits); + // Get a slice of samples based on the shot_range vector + size_t shot_idx = 0; + for (const auto &i : shot_range) { + for (size_t j = i * num_qubits; j < (i + 1) * num_qubits; j++) { + // TODO some extra work to make it cache-friendly + sub_samples[shot_idx * num_qubits + j - i * num_qubits] = + samples[j]; + } + shot_idx++; + } + return sub_samples; + } + return samples; + } +}; +} // namespace Pennylane::Measures diff --git a/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp b/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp index 7626ea835c..2f6d8adce8 100644 --- a/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp +++ b/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp @@ -193,6 +193,123 @@ TEST_CASE("Expval - NamedObs", "[MeasurementsBase][Observables]") { } } +template void testNamedObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + + // Defining the State Vector that will be measured. + auto statevector_data = createNonTrivialState(); + StateVectorT statevector(statevector_data.data(), + statevector_data.size()); + + // Initializing the measures class. + // This object attaches to the statevector allowing several measures. + Measurements Measurer(statevector); + + std::vector> wires_list = {{0}, {1}, {2}}; + std::vector obs_name = {"PauliX", "PauliY", "PauliZ", + "Hadamard", "Identity"}; + // Expected results calculated with Pennylane default.qubit: + std::vector> exp_values_ref = { + {0.49272486, 0.42073549, 0.28232124}, + {-0.64421768, -0.47942553, -0.29552020}, + {0.58498357, 0.77015115, 0.91266780}, + {0.7620549436, 0.8420840225, 0.8449848566}, + {1.0, 1.0, 1.0}}; + for (size_t ind_obs = 0; ind_obs < obs_name.size(); ind_obs++) { + DYNAMIC_SECTION(obs_name[ind_obs] + << " - Varying wires" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + for (size_t ind_wires = 0; ind_wires < wires_list.size(); + ind_wires++) { + NamedObs obs(obs_name[ind_obs], + wires_list[ind_wires]); + PrecisionT expected = exp_values_ref[ind_obs][ind_wires]; + PrecisionT result = + Measurer.expval(obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + } + } + + for (size_t ind_obs = 0; ind_obs < obs_name.size(); ind_obs++) { + DYNAMIC_SECTION(obs_name[ind_obs] + << " - Varying wires-with shots_range" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + for (size_t ind_wires = 0; ind_wires < wires_list.size(); + ind_wires++) { + NamedObs obs(obs_name[ind_obs], + wires_list[ind_wires]); + PrecisionT expected = exp_values_ref[ind_obs][ind_wires]; + PrecisionT result = + Measurer.expval(obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + } + } + testNamedObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot- NamedObs", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testNamedObsExpvalShot(); + } +} + +template void testHermitianObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + using MatrixT = std::vector; + + // Defining the State Vector that will be measured. + auto statevector_data = createNonTrivialState(); + StateVectorT statevector(statevector_data.data(), + statevector_data.size()); + + // Initializing the measures class. + // This object attaches to the statevector allowing several measures. + Measurements Measurer(statevector); + + const PrecisionT theta = M_PI / 2; + const PrecisionT real_term = std::cos(theta); + const PrecisionT imag_term = std::sin(theta); + + DYNAMIC_SECTION("Failed for Hermitian" + << StateVectorToName::name) { + MatrixT Hermitian_matrix{real_term, ComplexT{0, imag_term}, + ComplexT{0, -imag_term}, real_term}; + + HermitianObs obs(Hermitian_matrix, {0}); + size_t num_shots = 1000; + std::vector shots_range = {}; + REQUIRE_THROWS_WITH( + Measurer.expval(obs, num_shots, shots_range), + Catch::Matchers::Contains( + "expval calculation is not supported by shots")); + REQUIRE(obs.getCoeffs().size() == 0); + } + + testHermitianObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot - HermitianObs ", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testHermitianObsExpvalShot(); + } +} + template void testHermitianObsExpval() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; @@ -270,6 +387,95 @@ TEST_CASE("Expval - HermitianObs", "[MeasurementsBase][Observables]") { } } +template void testTensorProdObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = StateVectorT::ComplexT; + + // Defining the State Vector that will be measured. + std::vector 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}}; + StateVectorT statevector(statevector_data.data(), + statevector_data.size()); + + // Initializing the measures class. + // This object attaches to the statevector allowing several measures. + Measurements Measurer(statevector); + + DYNAMIC_SECTION(" - Without shots_range" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = TensorProdObs::create({X0, Z1}); + auto expected = PrecisionT(-0.36); + auto result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" - With Identity but no shots_range" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = TensorProdObs::create({X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With shots_range" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = TensorProdObs::create({X0, Z1}); + auto expected = PrecisionT(-0.36); + auto result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With Identity and shots_range" + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = TensorProdObs::create({X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + testTensorProdObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot- TensorProdObs", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testTensorProdObsExpvalShot(); + } +} + template void testNamedObsVar() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; @@ -456,4 +662,101 @@ TEST_CASE("Samples", "[MeasurementsBase]") { if constexpr (BACKEND_FOUND) { testSamples(); } +} + +template void testHamiltonianObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + + // Defining the State Vector that will be measured. + std::vector 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}}; + StateVectorT statevector(statevector_data.data(), + statevector_data.size()); + + // Initializing the measures class. + // This object attaches to the statevector allowing several measures. + Measurements Measurer(statevector); + + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + + auto ob = Hamiltonian::create({0.3, 0.5}, {X0, Z1}); + + DYNAMIC_SECTION("Without shots_range " + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + + auto res = Measurer.expval(*ob, num_shots, shots_range); + auto expected = PrecisionT(-0.086); + REQUIRE(expected == Approx(res).margin(5e-2)); + } + + DYNAMIC_SECTION("With shots_range " + << StateVectorToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + + auto res = Measurer.expval(*ob, num_shots, shots_range); + auto expected = PrecisionT(-0.086); + REQUIRE(expected == Approx(res).margin(5e-2)); + } + + testHamiltonianObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot - HamiltonianObs ", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testHamiltonianObsExpvalShot(); + } +} + +template void testSparseHObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using ComplexT = typename StateVectorT::ComplexT; + + // Defining the State Vector that will be measured. + auto statevector_data = createNonTrivialState(); + StateVectorT statevector(statevector_data.data(), + statevector_data.size()); + + // Initializing the measures class. + // This object attaches to the statevector allowing several measures. + Measurements Measurer(statevector); + + auto sparseH = SparseHamiltonian::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, {0, 1, 2}); + + DYNAMIC_SECTION("Failed for SparseH " + << StateVectorToName::name) { + size_t num_shots = 1000; + std::vector shots_range = {}; + REQUIRE_THROWS_WITH( + Measurer.expval(*sparseH, num_shots, shots_range), + Catch::Matchers::Contains( + "expval calculation is not supported by shots")); + } + + testSparseHObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot - SparseHObs ", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testSparseHObsExpvalShot(); + } } \ No newline at end of file diff --git a/pennylane_lightning/core/src/measurements/tests/mpi/Test_MeasurementsBaseMPI.cpp b/pennylane_lightning/core/src/measurements/tests/mpi/Test_MeasurementsBaseMPI.cpp index b459f762a5..7996c029c8 100644 --- a/pennylane_lightning/core/src/measurements/tests/mpi/Test_MeasurementsBaseMPI.cpp +++ b/pennylane_lightning/core/src/measurements/tests/mpi/Test_MeasurementsBaseMPI.cpp @@ -187,6 +187,82 @@ TEST_CASE("Expval - NamedObs", "[MeasurementsBase][Observables]") { } } +template void testNamedObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + + // Defining the State Vector that will be measured. + auto statevector_data = + createNonTrivialState>(); + + size_t num_qubits = 3; + + MPIManager mpi_manager(MPI_COMM_WORLD); + REQUIRE(mpi_manager.getSize() == 2); + + size_t mpi_buffersize = 1; + + size_t nGlobalIndexBits = + std::bit_width(static_cast(mpi_manager.getSize())) - 1; + size_t nLocalIndexBits = num_qubits - nGlobalIndexBits; + + int nDevices = 0; + cudaGetDeviceCount(&nDevices); + REQUIRE(nDevices >= 2); + int deviceId = mpi_manager.getRank() % nDevices; + cudaSetDevice(deviceId); + DevTag dt_local(deviceId, 0); + mpi_manager.Barrier(); + + auto sv_data_local = mpi_manager.scatter(statevector_data, 0); + + StateVectorT sv(mpi_manager, dt_local, mpi_buffersize, nGlobalIndexBits, + nLocalIndexBits); + sv.CopyHostDataToGpu(sv_data_local.data(), sv_data_local.size(), false); + mpi_manager.Barrier(); + // Initializing the measurements class. + // This object attaches to the statevector allowing several measures. + MeasurementsMPI Measurer(sv); + + std::vector> wires_list = {{0}, {1}, {2}}; + std::vector obs_name = {"PauliX", "PauliY", "PauliZ", + "Hadamard"}; + // Expected results calculated with Pennylane default.qubit: + std::vector> exp_values_ref = { + {0.49272486, 0.42073549, 0.28232124}, + {-0.64421768, -0.47942553, -0.29552020}, + {0.58498357, 0.77015115, 0.91266780}, + {0.7620549436, 0.8420840225, 0.8449848566}}; + + size_t num_shots = 10000; + std::vector shots_range = {}; + + for (size_t ind_obs = 0; ind_obs < obs_name.size(); ind_obs++) { + DYNAMIC_SECTION(obs_name[ind_obs] + << " - Varying wires" + << StateVectorMPIToName::name) { + for (size_t ind_wires = 0; ind_wires < wires_list.size(); + ind_wires++) { + NamedObsMPI obs(obs_name[ind_obs], + wires_list[ind_wires]); + PrecisionT expected = exp_values_ref[ind_obs][ind_wires]; + PrecisionT result = + Measurer.expval(obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + } + } + testNamedObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot- NamedObs", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testNamedObsExpvalShot(); + } +} + template void testHermitianObsExpval() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; @@ -289,6 +365,234 @@ TEST_CASE("Expval - HermitianObs", "[MeasurementsBase][Observables]") { } } +template void testTensorProdObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using ComplexT = typename StateVectorT::ComplexT; + using PrecisionT = typename StateVectorT::PrecisionT; + + // Defining the State Vector that will be measured. + std::vector 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}}; + + size_t num_qubits = 3; + + MPIManager mpi_manager(MPI_COMM_WORLD); + REQUIRE(mpi_manager.getSize() == 2); + + size_t mpi_buffersize = 1; + + size_t nGlobalIndexBits = + std::bit_width(static_cast(mpi_manager.getSize())) - 1; + size_t nLocalIndexBits = num_qubits - nGlobalIndexBits; + + int nDevices = 0; + cudaGetDeviceCount(&nDevices); + REQUIRE(nDevices >= 2); + int deviceId = mpi_manager.getRank() % nDevices; + cudaSetDevice(deviceId); + DevTag dt_local(deviceId, 0); + mpi_manager.Barrier(); + + auto sv_data_local = mpi_manager.scatter(statevector_data, 0); + + StateVectorT sv(mpi_manager, dt_local, mpi_buffersize, nGlobalIndexBits, + nLocalIndexBits); + sv.CopyHostDataToGpu(sv_data_local.data(), sv_data_local.size(), false); + mpi_manager.Barrier(); + // Initializing the measurements class. + // This object attaches to the statevector allowing several measures. + MeasurementsMPI Measurer(sv); + + DYNAMIC_SECTION(" - Without shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = TensorProdObsMPI::create({X0, Z1}); + auto expected = PrecisionT(-0.36); + auto result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" - With Identity but no shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = TensorProdObsMPI::create({X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = TensorProdObsMPI::create({X0, Z1}); + auto expected = PrecisionT(-0.36); + auto result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With Identity and shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = TensorProdObsMPI::create({X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + testTensorProdObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot- TensorProdObs", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testTensorProdObsExpvalShot(); + } +} + +template void testHamiltonianObsExpvalShot() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using ComplexT = typename StateVectorT::ComplexT; + using PrecisionT = typename StateVectorT::PrecisionT; + + // Defining the State Vector that will be measured. + std::vector 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}}; + + size_t num_qubits = 3; + + MPIManager mpi_manager(MPI_COMM_WORLD); + REQUIRE(mpi_manager.getSize() == 2); + + size_t mpi_buffersize = 1; + + size_t nGlobalIndexBits = + std::bit_width(static_cast(mpi_manager.getSize())) - 1; + size_t nLocalIndexBits = num_qubits - nGlobalIndexBits; + + int nDevices = 0; + cudaGetDeviceCount(&nDevices); + REQUIRE(nDevices >= 2); + int deviceId = mpi_manager.getRank() % nDevices; + cudaSetDevice(deviceId); + DevTag dt_local(deviceId, 0); + mpi_manager.Barrier(); + + auto sv_data_local = mpi_manager.scatter(statevector_data, 0); + + StateVectorT sv(mpi_manager, dt_local, mpi_buffersize, nGlobalIndexBits, + nLocalIndexBits); + sv.CopyHostDataToGpu(sv_data_local.data(), sv_data_local.size(), false); + mpi_manager.Barrier(); + // Initializing the measurements class. + // This object attaches to the statevector allowing several measures. + MeasurementsMPI Measurer(sv); + + DYNAMIC_SECTION(" - Without shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = + HamiltonianMPI::create({0.3, 0.5}, {X0, Z1}); + PrecisionT expected = PrecisionT(-0.086); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" - With Identity but no shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range = {}; + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = + HamiltonianMPI::create({0.3, 0.5}, {X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto Z1 = std::make_shared>( + "PauliZ", std::vector{1}); + auto obs = + HamiltonianMPI::create({0.3, 0.5}, {X0, Z1}); + PrecisionT expected = PrecisionT(-0.086); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + DYNAMIC_SECTION(" With Identity and shots_range" + << StateVectorMPIToName::name) { + size_t num_shots = 10000; + std::vector shots_range; + for (size_t i = 0; i < num_shots; i += 2) { + shots_range.push_back(i); + } + auto X0 = std::make_shared>( + "PauliX", std::vector{0}); + auto I1 = std::make_shared>( + "Identity", std::vector{1}); + auto obs = + HamiltonianMPI::create({0.3, 0.5}, {X0, I1}); + PrecisionT expected = Measurer.expval(*obs); + PrecisionT result = Measurer.expval(*obs, num_shots, shots_range); + REQUIRE(expected == Approx(result).margin(5e-2)); + } + + testHamiltonianObsExpvalShot(); + } +} + +TEST_CASE("Expval Shot- HamiltonianObs", "[MeasurementsBase][Observables]") { + if constexpr (BACKEND_FOUND) { + testHamiltonianObsExpvalShot(); + } +} + template void testNamedObsVar() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; diff --git a/pennylane_lightning/core/src/observables/Observables.hpp b/pennylane_lightning/core/src/observables/Observables.hpp index 183a9fcd82..45d8aedf4a 100644 --- a/pennylane_lightning/core/src/observables/Observables.hpp +++ b/pennylane_lightning/core/src/observables/Observables.hpp @@ -61,6 +61,22 @@ template class Observable { */ virtual void applyInPlace(StateVectorT &sv) const = 0; + /** + * @brief Apply unitaries of an observable to the given statevector in + * place. + * + * @param sv Reference to StateVector object. + * @param identity_wires Reference to a std::vector object which stores + * wires of Identity gates in the observable. + * @param ob_wires Reference to a std::vector object which stores wires of + * the observable. + * @param term_idx Index of a Hamiltonian term. + */ + virtual void applyInPlaceShots(StateVectorT &sv, + std::vector &identity_wires, + std::vector &ob_wires, + size_t term_idx = 0) const = 0; + /** * @brief Get the name of the observable */ @@ -71,6 +87,13 @@ template class Observable { */ [[nodiscard]] virtual auto getWires() const -> std::vector = 0; + /** + * @brief Get the coefficients of a Hamiltonian observable. + */ + [[nodiscard]] virtual auto getCoeffs() const -> std::vector { + return {}; + }; + /** * @brief Test whether this object is equal to another object */ @@ -140,6 +163,32 @@ class NamedObsBase : public Observable { void applyInPlace(StateVectorT &sv) const override { sv.applyOperation(obs_name_, wires_, false, params_); } + + void + applyInPlaceShots(StateVectorT &sv, std::vector &identity_wire, + std::vector &ob_wires, + [[maybe_unused]] size_t term_idx = 0) const override { + ob_wires.clear(); + identity_wire.clear(); + ob_wires.push_back(wires_[0]); + + if (obs_name_ == "PauliX") { + sv.applyOperation("Hadamard", wires_, false); + } else if (obs_name_ == "PauliY") { + sv.applyOperations({"PauliZ", "S", "Hadamard"}, + {wires_, wires_, wires_}, {false, false, false}); + } else if (obs_name_ == "Hadamard") { + const PrecisionT theta = -M_PI / 4.0; + sv.applyOperation("RY", wires_, false, {theta}); + } else if (obs_name_ == "PauliZ") { + } else if (obs_name_ == "Identity") { + identity_wire.push_back(wires_[0]); + } else { + PL_ABORT("Provided NamedObs does not supported for shots " + "calculation. Supported NamedObs are PauliX, PauliY, " + "PauliZ, Identity and Hadamard."); + } + } }; /** @@ -192,6 +241,15 @@ class HermitianObsBase : public Observable { void applyInPlace(StateVectorT &sv) const override { sv.applyMatrix(matrix_, wires_); } + + void + applyInPlaceShots([[maybe_unused]] StateVectorT &sv, + [[maybe_unused]] std::vector &identity_wire, + [[maybe_unused]] std::vector &ob_wires, + [[maybe_unused]] size_t term_idx = 0) const override { + PL_ABORT("Hermitian observables do not support applyInPlaceShots " + "method."); + } }; /** @@ -224,6 +282,7 @@ class TensorProdObsBase : public Observable { } public: + using PrecisionT = typename StateVectorT::PrecisionT; /** * @brief Create a tensor product of observables * @@ -301,6 +360,23 @@ class TensorProdObsBase : public Observable { } } + void + applyInPlaceShots(StateVectorT &sv, std::vector &identity_wires, + std::vector &ob_wires, + [[maybe_unused]] size_t term_idx = 0) const override { + identity_wires.clear(); + ob_wires.clear(); + for (const auto &ob : obs_) { + std::vector identity_wire; + std::vector ob_wire; + ob->applyInPlaceShots(sv, identity_wire, ob_wire); + if (!identity_wire.empty()) { + identity_wires.push_back(identity_wire[0]); + } + ob_wires.push_back(ob_wire[0]); + } + } + [[nodiscard]] auto getObsName() const -> std::string override { using Util::operator<<; std::ostringstream obs_stream; @@ -386,6 +462,15 @@ class HamiltonianBase : public Observable { "defined at the backend level."); } + void + applyInPlaceShots([[maybe_unused]] StateVectorT &sv, + [[maybe_unused]] std::vector &identity_wires, + [[maybe_unused]] std::vector &ob_wires, + [[maybe_unused]] size_t term_idx = 0) const override { + PL_ABORT("For Hamiltonian Observables, the applyInPlace method must be " + "defined at the backend level."); + } + [[nodiscard]] auto getWires() const -> std::vector override { std::unordered_set wires; @@ -412,6 +497,13 @@ class HamiltonianBase : public Observable { ss << "]}"; return ss.str(); } + + /** + * @brief Get the coefficients of the observable. + */ + [[nodiscard]] auto getCoeffs() const -> std::vector override { + return coeffs_; + }; }; /** @@ -496,6 +588,15 @@ class SparseHamiltonianBase : public Observable { "defined at the backend level."); } + void + applyInPlaceShots([[maybe_unused]] StateVectorT &sv, + [[maybe_unused]] std::vector &identity_wire, + [[maybe_unused]] std::vector &ob_wires, + [[maybe_unused]] size_t term_idx = 0) const override { + PL_ABORT("SparseHamiltonian observables do not the applyInPlaceShots " + "method."); + } + [[nodiscard]] auto getObsName() const -> std::string override { using Pennylane::Util::operator<<; std::ostringstream ss; diff --git a/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp b/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp index d0bb9799b1..c928fba3d4 100644 --- a/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp +++ b/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp @@ -21,6 +21,7 @@ #include #include +#include #include /** * @file @@ -82,6 +83,7 @@ template struct StateVectorToName {}; template void testNamedObsBase() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; using NamedObsT = NamedObsBase; DYNAMIC_SECTION("Name of the Observable must be correct - " @@ -119,6 +121,24 @@ template void testNamedObsBase() { REQUIRE(ob1 != ob3); } + DYNAMIC_SECTION("Unsupported NamedObs for applyInPlaceShots") { + std::mt19937_64 re{1337}; + const size_t num_qubits = 3; + auto init_state = + createRandomStateVectorData(re, num_qubits); + + StateVectorT state_vector(init_state.data(), init_state.size()); + auto obs = NamedObsT("RY", {0}, {0.4}); + + std::vector identity_wire; + std::vector ob_wires; + + REQUIRE_THROWS_WITH( + obs.applyInPlaceShots(state_vector, identity_wire, ob_wires), + Catch::Matchers::Contains( + "Provided NamedObs does not supported for shots")); + } + testNamedObsBase(); } } @@ -133,6 +153,7 @@ template void testHermitianObsBase() { if constexpr (!std::is_same_v) { using StateVectorT = typename TypeList::Type; using ComplexT = typename StateVectorT::ComplexT; + using PrecisionT = typename StateVectorT::PrecisionT; using HermitianObsT = HermitianObsBase; DYNAMIC_SECTION("HermitianObs only accepts correct arguments - " @@ -183,6 +204,26 @@ template void testHermitianObsBase() { REQUIRE(ob2 != ob3); } + DYNAMIC_SECTION("Failed for HermitianObs for applyInPlaceShots - " + << StateVectorToName::name) { + std::mt19937_64 re{1337}; + const size_t num_qubits = 3; + auto init_state = + createRandomStateVectorData(re, num_qubits); + + StateVectorT state_vector(init_state.data(), init_state.size()); + auto obs = + HermitianObsT{std::vector{1.0, 0.0, -1.0, 0.0}, {0}}; + + std::vector identity_wire; + std::vector ob_wires; + + REQUIRE_THROWS_WITH( + obs.applyInPlaceShots(state_vector, identity_wire, ob_wires), + Catch::Matchers::Contains("Hermitian observables do not " + "support applyInPlaceShots method.")); + } + testHermitianObsBase(); } } @@ -453,6 +494,22 @@ template void testHamiltonianBase() { REQUIRE_THROWS_AS(ham->applyInPlace(state_vector), LightningException); } + + DYNAMIC_SECTION("applyInPlaceShots must fail - " + << StateVectorToName::name) { + auto ham = + HamiltonianT::create({PrecisionT{1.0}, h, h}, {zz, x1, x2}); + auto st_data = createZeroState(2); + + StateVectorT state_vector(st_data.data(), st_data.size()); + + std::vector identity_wires; + std::vector ob_wires; + + REQUIRE_THROWS_AS(ham->applyInPlaceShots( + state_vector, identity_wires, ob_wires), + LightningException); + } } testHamiltonianBase(); } @@ -535,6 +592,23 @@ template void testSparseHamiltonianBase() { LightningException); } + DYNAMIC_SECTION("SparseHamiltonianBase - applyInPlaceShots must fail - " + << StateVectorToName::name) { + auto init_state = + createRandomStateVectorData(re, num_qubits); + + StateVectorT state_vector(init_state.data(), init_state.size()); + + std::vector identity_wire; + std::vector ob_wires; + + REQUIRE_THROWS_WITH( + sparseH->applyInPlaceShots(state_vector, identity_wire, + ob_wires), + Catch::Matchers::Contains("SparseHamiltonian observables do " + "not the applyInPlaceShots method.")); + } + testSparseHamiltonianBase(); } } diff --git a/pennylane_lightning/core/src/simulators/base/StateVectorBase.hpp b/pennylane_lightning/core/src/simulators/base/StateVectorBase.hpp index 727447f8e6..0baa89d0ae 100644 --- a/pennylane_lightning/core/src/simulators/base/StateVectorBase.hpp +++ b/pennylane_lightning/core/src/simulators/base/StateVectorBase.hpp @@ -66,6 +66,15 @@ template class StateVectorBase { return num_qubits_; } + /** + * @brief Get the total number of qubits of the simulated system. + * + * @return std::size_t + */ + [[nodiscard]] auto getTotalNumQubits() const -> size_t { + return num_qubits_; + } + /** * @brief Get the size of the statevector * diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaMPI.hpp index ca2c564e63..268d7332ae 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaMPI.hpp @@ -38,6 +38,8 @@ #include "cuGates_host.hpp" #include "cuda_helpers.hpp" +#include "CPUMemoryModel.hpp" + #include "LinearAlg.hpp" /// @cond DEV @@ -97,6 +99,7 @@ class StateVectorCudaMPI final StateVectorCudaMPI>::CFP_t; using PrecisionT = Precision; using ComplexT = std::complex; + using MemoryStorageT = Pennylane::Util::MemoryStorageLocation::Undefined; StateVectorCudaMPI() = delete; @@ -204,6 +207,8 @@ class StateVectorCudaMPI final handle_.get(), mpi_manager_, 0, BaseType::getData(), numLocalQubits_, localStream_.get())), gate_cache_(true, other.getDataBuffer().getDevTag()) { + PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); + mpi_manager_.Barrier(); BaseType::CopyGpuDataToGpuIn(other.getData(), other.getLength(), false); PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); mpi_manager_.Barrier(); diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaManaged.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaManaged.hpp index faddef3710..7f1fac2652 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaManaged.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/StateVectorCudaManaged.hpp @@ -34,6 +34,8 @@ #include "cuGates_host.hpp" #include "cuda_helpers.hpp" +#include "CPUMemoryModel.hpp" + #include "cuError.hpp" #include "LinearAlg.hpp" @@ -81,6 +83,7 @@ class StateVectorCudaManaged using CFP_t = typename StateVectorCudaBase>::CFP_t; + using MemoryStorageT = Pennylane::Util::MemoryStorageLocation::Undefined; StateVectorCudaManaged() = delete; StateVectorCudaManaged(size_t num_qubits) diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp index 8ca6eacc69..78ab5b53e8 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp @@ -350,6 +350,21 @@ class Measurements final return static_cast(expect); } + /** + * @brief Expectation value for a Observable with shots + * + * @param obs Observable. + * @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. + * @return Floating point expected value of the observable. + */ + + auto expval(const Observable &obs, const size_t &num_shots, + const std::vector &shot_range) -> PrecisionT { + return BaseType::expval(obs, num_shots, shot_range); + } + /** * @brief Expected value of an observable. * diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp index ff101654df..c4bac79435 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp @@ -87,6 +87,7 @@ class MeasurementsMPI final } else { data_type_ = CUDA_C_32F; } + mpi_manager_.Barrier(); }; /** @@ -291,6 +292,7 @@ class MeasurementsMPI final precumulative = 0; } PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); + mpi_manager_.Barrier(); // Ensure the 'custatevecSamplerApplySubSVOffset' function can be called // successfully without reducing accuracy. @@ -320,6 +322,7 @@ class MeasurementsMPI final if (mpi_manager_.getRank() == 0) { preshotOffset = 0; } + mpi_manager_.Barrier(); int nSubShots = shotOffset - preshotOffset; if (nSubShots > 0) { @@ -341,6 +344,9 @@ class MeasurementsMPI final PL_CUDA_IS_SUCCESS(cudaFree(extraWorkspace)); } + PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); + mpi_manager_.Barrier(); + mpi_manager_.Allreduce(localBitStrings, globalBitStrings, "sum"); @@ -350,6 +356,7 @@ class MeasurementsMPI final (globalBitStrings[i] >> j) & 1U; } } + mpi_manager_.Barrier(); return samples; } @@ -486,6 +493,24 @@ class MeasurementsMPI final return static_cast(expect); } + /** + * @brief Expectation value for a Observable with shots + * + * @param obs Observable. + * @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. + * @return Floating point expected value of the observable. + */ + + auto expval(const Observable &obs, const size_t &num_shots, + const std::vector &shot_range) -> PrecisionT { + mpi_manager_.Barrier(); + PrecisionT result = BaseType::expval(obs, num_shots, shot_range); + mpi_manager_.Barrier(); + return result; + } + /** * @brief Expected value of an observable. * diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp index 6d0d0a94e7..b5a6c3deb0 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp @@ -207,6 +207,16 @@ class Hamiltonian final : public HamiltonianBase { } sv.updateData(std::move(buffer)); } + + // to work with + void applyInPlaceShots(StateVectorT &sv, + std::vector &identity_wires, + std::vector &ob_wires, + size_t term_idx) const override { + ob_wires.clear(); + this->obs_[term_idx]->applyInPlaceShots(sv, identity_wires, ob_wires, + term_idx); + } }; /** diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp index 955365915a..9595783421 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp @@ -214,6 +214,16 @@ class HamiltonianMPI final : public HamiltonianBase { PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); mpi_manager.Barrier(); } + + // to work with + void applyInPlaceShots(StateVectorT &sv, + std::vector &identity_wires, + std::vector &ob_wires, + size_t term_idx) const override { + ob_wires.clear(); + this->obs_[term_idx]->applyInPlaceShots(sv, identity_wires, ob_wires, + term_idx); + } }; /** diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp index 02064e811f..f33d3a8730 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp @@ -35,6 +35,8 @@ #include "StateVectorBase.hpp" #include "Util.hpp" +#include "CPUMemoryModel.hpp" + /// @cond DEV namespace { using Pennylane::Gates::GateOperation; @@ -90,6 +92,7 @@ class StateVectorKokkos final Kokkos::View>; using TeamPolicy = Kokkos::TeamPolicy<>; + using MemoryStorageT = Pennylane::Util::MemoryStorageLocation::Undefined; StateVectorKokkos() = delete; StateVectorKokkos(size_t num_qubits, diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp index c9414921f6..e4212d5b15 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp @@ -275,6 +275,20 @@ class Measurements final return expected_value_list; } + /** + * @brief Expectation value for a Observable with shots + * + * @param obs Observable. + * @param num_shots Number of shots. + * @param shots_range Vector of shot number to measurement. + * @return Floating point expected value of the observable. + */ + + auto expval(const Observable &obs, const size_t &num_shots, + const std::vector &shot_range) -> PrecisionT { + return BaseType::expval(obs, num_shots, shot_range); + } + /** * @brief Expected value of a Sparse Hamiltonian. * diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp index c3fae6b3ea..44357eb9e5 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp @@ -198,6 +198,16 @@ class Hamiltonian final : public HamiltonianBase { } sv.updateData(buffer); } + + // to work with + void applyInPlaceShots(StateVectorT &sv, + std::vector &identity_wires, + std::vector &ob_wires, + size_t term_idx) const override { + ob_wires.clear(); + this->obs_[term_idx]->applyInPlaceShots(sv, identity_wires, ob_wires, + term_idx); + } }; /** diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp index 534e00d112..e1d3ef551d 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp @@ -268,6 +268,20 @@ class Measurements final return result; } + /** + * @brief Expectation value for a Observable with shots + * + * @param obs Observable + * @param num_shots Number of shots. + * @param shots_range Vector of shot number to measurement + * @return Floating point expected value of the observable. + */ + + auto expval(const Observable &obs, const size_t &num_shots, + const std::vector &shot_range) -> PrecisionT { + return BaseType::expval(obs, num_shots, shot_range); + } + /** * @brief Variance value for a general Observable * diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp index b659969037..ce1aa1cca5 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp @@ -358,6 +358,16 @@ class Hamiltonian final : public HamiltonianBase { StateVectorT, Pennylane::Util::use_openmp>::run(this->coeffs_, this->obs_, sv); } + + // to work with + void applyInPlaceShots(StateVectorT &sv, + std::vector &identity_wires, + std::vector &ob_wires, + size_t term_idx) const override { + ob_wires.clear(); + this->obs_[term_idx]->applyInPlaceShots(sv, identity_wires, ob_wires, + term_idx); + } }; /**