Program Listing for File MeasurementsBase.hpp

Return to documentation for file (pennylane_lightning/core/src/measurements/MeasurementsBase.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 <random>
#include <string>
#include <vector>

#include "Observables.hpp"

#include "CPUMemoryModel.hpp"

#include "Util.hpp"

namespace {
using namespace Pennylane::Observables;
using namespace Pennylane::Util;
} // namespace

namespace Pennylane::Measures {
template <class StateVectorT, class Derived> class MeasurementsBase {
  private:
    using PrecisionT = typename StateVectorT::PrecisionT;
    using ComplexT = typename StateVectorT::ComplexT;

  protected:
#ifdef _ENABLE_PLGPU
    StateVectorT &_statevector;
#else
    const StateVectorT &_statevector;
#endif
    std::mt19937 rng;

  public:
#ifdef _ENABLE_PLGPU
    explicit MeasurementsBase(StateVectorT &statevector)
        : _statevector{statevector} {};
#else
    explicit MeasurementsBase(const StateVectorT &statevector)
        : _statevector{statevector} {};
#endif

    void setSeed(const std::size_t seed) { rng.seed(seed); }

    void setRandomSeed() {
        std::random_device rd;
        setSeed(rd());
    }

    auto expval(const Observable<StateVectorT> &obs) -> PrecisionT {
        return static_cast<Derived *>(this)->expval(obs);
    }

    auto var(const Observable<StateVectorT> &obs) -> PrecisionT {
        return static_cast<Derived *>(this)->var(obs);
    }

    auto probs() -> std::vector<PrecisionT> {
        return static_cast<Derived *>(this)->probs();
    };

    auto probs(const std::vector<std::size_t> &wires)
        -> std::vector<PrecisionT> {
        return static_cast<Derived *>(this)->probs(wires);
    };

    auto generate_samples(size_t num_samples) -> std::vector<std::size_t> {
        return static_cast<Derived *>(this)->generate_samples(num_samples);
    };

    auto expval(const Observable<StateVectorT> &obs,
                const std::size_t &num_shots,
                const std::vector<std::size_t> &shot_range = {}) -> PrecisionT {
        PrecisionT result{0.0};

        if (obs.getObsName().find("SparseHamiltonian") != std::string::npos) {
            // SparseHamiltonian does not support shot measurement in pennylane.
            PL_ABORT("SparseHamiltonian observables do not support shot "
                     "measurement.");
        } 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++) {
                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);
            result =
                std::accumulate(obs_samples.begin(), obs_samples.end(), 0.0);
            result /= obs_samples.size();
        }
        return result;
    }

    auto measure_with_samples(const Observable<StateVectorT> &obs,
                              const std::size_t &num_shots,
                              const std::vector<std::size_t> &shot_range)
        -> std::vector<PrecisionT> {
        const std::size_t num_qubits = _statevector.getTotalNumQubits();
        std::vector<std::size_t> obs_wires;
        std::vector<std::vector<PrecisionT>> eigenValues;

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

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

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

        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::size_t idx = 0;
            std::size_t wire_idx = 0;
            for (auto &obs_wire : obs_wires) {
                idx += sub_samples[i * num_qubits + obs_wire]
                       << (obs_wires.size() - 1 - wire_idx);
                wire_idx++;
            }

            obs_samples[i] = eigenVals[idx];
        }
        return obs_samples;
    }

    auto var(const Observable<StateVectorT> &obs, const std::size_t &num_shots)
        -> PrecisionT {
        PrecisionT result{0.0};
        if (obs.getObsName().find("SparseHamiltonian") != std::string::npos) {
            // SparseHamiltonian does not support shot measurement in pennylane.
            PL_ABORT("SparseHamiltonian observables do not support shot "
                     "measurement.");
        } else if (obs.getObsName().find("Hamiltonian") != std::string::npos) {
            // Branch for Hamiltonian observables
            auto coeffs = obs.getCoeffs();
            auto obs_terms = obs.getObs();

            std::size_t obs_term_idx = 0;
            for (const auto &coeff : coeffs) {
                result +=
                    coeff * coeff * var(*obs_terms[obs_term_idx], num_shots);
                obs_term_idx++;
            }
        } else {
            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();
            result = mean_square - square_mean * square_mean;
        }
        return result;
    }

    auto probs(const Observable<StateVectorT> &obs, std::size_t num_shots = 0)
        -> std::vector<PrecisionT> {
        PL_ABORT_IF(
            obs.getObsName().find("Hamiltonian") != std::string::npos,
            "Hamiltonian and Sparse Hamiltonian do not support samples().");
        std::vector<std::size_t> obs_wires;
        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 > std::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 > std::size_t{0}) {
                return measure.probs(obs_wires, num_shots);
            }
            return measure.probs(obs_wires);
        }
    }

    auto probs(const std::vector<std::size_t> &wires, std::size_t num_shots)
        -> std::vector<PrecisionT> {
        auto counts_map = counts(num_shots);

        std::size_t num_wires = _statevector.getTotalNumQubits();

        std::vector<PrecisionT> prob_shots(size_t{1} << wires.size(), 0.0);

        for (auto &it : counts_map) {
            std::size_t bitVal = 0;
            for (size_t bit = 0; bit < wires.size(); bit++) {
                // Mapping the value of wires[bit]th bit to local [bit]th bit of
                // the output
                bitVal +=
                    ((it.first >> (num_wires - std::size_t{1} - wires[bit])) &
                     std::size_t{1})
                    << (wires.size() - std::size_t{1} - bit);
            }

            prob_shots[bitVal] +=
                it.second / static_cast<PrecisionT>(num_shots);
        }

        return prob_shots;
    }

    auto probs(size_t num_shots) -> std::vector<PrecisionT> {
        auto counts_map = counts(num_shots);

        std::size_t num_wires = _statevector.getTotalNumQubits();

        std::vector<PrecisionT> prob_shots(size_t{1} << num_wires, 0.0);

        for (auto &it : counts_map) {
            prob_shots[it.first] =
                it.second / static_cast<PrecisionT>(num_shots);
        }

        return prob_shots;
    }

    auto sample(const Observable<StateVectorT> &obs,
                const std::size_t &num_shots) -> std::vector<PrecisionT> {
        PL_ABORT_IF(
            obs.getObsName().find("Hamiltonian") != std::string::npos,
            "Hamiltonian and Sparse Hamiltonian do not support samples().");
        std::vector<std::size_t> obs_wires;
        std::vector<std::size_t> shot_range = {};

        return measure_with_samples(obs, num_shots, shot_range);
    }

    auto sample(const std::size_t &num_shots) -> std::vector<std::size_t> {
        Derived measure(_statevector);
        return measure.generate_samples(num_shots);
    }

    auto counts(const Observable<StateVectorT> &obs,
                const std::size_t &num_shots)
        -> std::unordered_map<PrecisionT, std::size_t> {
        std::unordered_map<PrecisionT, std::size_t> outcome_map;
        auto sample_data = sample(obs, num_shots);
        for (size_t i = 0; i < num_shots; i++) {
            auto key = sample_data[i];
            auto it = outcome_map.find(key);
            if (it != outcome_map.end()) {
                it->second += 1;
            } else {
                outcome_map[key] = 1;
            }
        }
        return outcome_map;
    }

    auto counts(const std::size_t &num_shots)
        -> std::unordered_map<size_t, std::size_t> {
        std::unordered_map<size_t, std::size_t> outcome_map;
        auto sample_data = sample(num_shots);

        std::size_t num_wires = _statevector.getTotalNumQubits();
        for (size_t i = 0; i < num_shots; i++) {
            std::size_t key = 0;
            for (size_t j = 0; j < num_wires; j++) {
                key += sample_data[i * num_wires + j] << (num_wires - 1 - j);
            }

            auto it = outcome_map.find(key);
            if (it != outcome_map.end()) {
                it->second += 1;
            } else {
                outcome_map[key] = 1;
            }
        }
        return outcome_map;
    }

  private:
    auto _sample_state(const Observable<StateVectorT> &obs,
                       const std::size_t &num_shots,
                       const std::vector<std::size_t> &shot_range,
                       std::vector<std::size_t> &obs_wires,
                       std::vector<std::vector<PrecisionT>> &eigenValues)
        -> std::vector<std::size_t> {
        const std::size_t num_qubits = _statevector.getTotalNumQubits();
        std::vector<std::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<std::size_t> sub_samples(shot_range.size() *
                                                 num_qubits);
            // Get a slice of samples based on the shot_range vector
            std::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