State.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _SIMULATOR_STATE_H_
#define _SIMULATOR_STATE_H_
 
#include <Eigen/Eigen>
#include <cmath>
#include <complex>
#include <cstdint>
#include <random>
#include <stdexcept>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
#ifndef NO_QISKIT_AER
#include "framework/linalg/vector.hpp"
#endif
 
#include "../../Utils/RandomStream.h"
#include "../Core/QuantumChannel.h"
#include "../TensorNetworks/TensorQueries.h"
#include "SimulatorObserver.h"
 
namespace Circuits
{
template <typename Time> class IOperation;
}
 
namespace Simulators
{

/*
template<class T> class avoid_init_allocator : public std::allocator<T>
{
public:
        using std::allocator<T>::allocator;
 
        template <class U, class... Args> void construct(U*, Args&&...) {}
};
*/

enum class SimulatorType : int
{
#ifndef NO_QISKIT_AER
    kQiskitAer, 
#endif
    kQCSim, 
#ifndef NO_QISKIT_AER
    kCompositeQiskitAer, 
#endif
    kCompositeQCSim, 
    kGpuSim,         
    kQuestSim,       
    kDistGpuSim,     
    kDistMpiGpuSim   
};
 
inline bool IsDistributedGpuSimulator(SimulatorType type)
{
    return type == SimulatorType::kDistGpuSim || type == SimulatorType::kDistMpiGpuSim;
}
 
inline bool IsGpuSimulator(SimulatorType type)
{
    return type == SimulatorType::kGpuSim || IsDistributedGpuSimulator(type);
}

enum class SimulationType : int
{
    kStatevector,           
    kMatrixProductState,    
    kStabilizer,            
    kTensorNetwork,         
    kPauliPropagator,       
    kExtendedStabilizer,    
    kPathIntegral,          
    kDensityMatrix,         
    kMatrixProductOperator, 
    kOther                  
};

class IState
{
  public:
    virtual ~IState() = default;

    virtual void SetSeed(uint64_t seed)
    {
        SeedAuxiliaryRng(seed);
        const std::string value = std::to_string(seed);
        Configure("seed", value.c_str());
    }

    double RandomUniform()
    {
        return auxUniform(auxRng);
    }
 
    static uint64_t DeriveSeed(uint64_t seed, uint64_t stream)
    {
        return Utils::RandomStream::Derive(seed, stream);
    }

    virtual void Initialize() = 0;

    virtual void InitializeState(size_t num_qubits, std::vector<std::complex<double>> &amplitudes) = 0;

#ifndef NO_QISKIT_AER
    virtual void InitializeState(size_t num_qubits, AER::Vector<std::complex<double>> &amplitudes) = 0;
#endif

    virtual void InitializeState(size_t num_qubits, Eigen::VectorXcd &amplitudes) = 0;

    virtual void InitializeToBasisState(size_t num_qubits, Types::qubit_t basisState)
    {
        (void)num_qubits;
        (void)basisState;
        throw std::runtime_error("This simulator does not support initialization to a computational "
                                 "basis state");
    }

    virtual void InitializeToBasisState(size_t num_qubits, const std::vector<bool> &basisState)
    {
        (void)num_qubits;
        (void)basisState;
        throw std::runtime_error("This simulator does not support initialization to a computational "
                                 "basis state");
    }

    virtual void InitializeToMixtureOfBasisStates(size_t num_qubits, const std::vector<std::pair<Types::qubit_t, double>> &mixture)
    {
        (void)num_qubits;
        (void)mixture;
        throw std::runtime_error("This simulator does not support initialization to a mixture of "
                                 "computational basis states");
    }

    virtual void InitializeToMixtureOfBasisStates(size_t num_qubits, const std::vector<std::pair<std::vector<bool>, double>> &mixture)
    {
        (void)num_qubits;
        (void)mixture;
        throw std::runtime_error("This simulator does not support initialization to a mixture of "
                                 "computational basis states");
    }

    virtual void Reset() = 0;

    virtual bool SupportsMPSSwapOptimization() const
    {
        return false;
    }
 
    // Named SWAP gates exchange logical labels without moving tensor sites.
    // Internal routing swaps still require tensor contractions.
    virtual bool SwapGatesRelabelQubits() const
    {
        return false;
    }
 
    // True when external meeting-position lookahead can use an upcoming list.
    // Requires a live MPS/MPO, enabled lookahead and an authoritative qubit map.
    virtual bool IsRoutingLookaheadEnabled() const
    {
        return false;
    }

    virtual void SetInitialQubitsMap(const std::vector<long long int> &initialMap)
    {
    }

    virtual void SetUseOptimalMeetingPosition(bool /*enable*/)
    {
    }

    virtual void SetLookaheadDepth(int /*depth*/)
    {
    }
 
    // Network execution calls this after computing the initial layout and
    // resolving its automatic depth. Explicit SetLookaheadDepth calls retain
    // their original meaning, including zero for immediate-cost routing.
    virtual void SetAutomaticLookaheadDepth(int depth)
    {
        SetLookaheadDepth(depth);
    }

    virtual void SetLookaheadDepthWithHeuristic(int /*depth*/)
    {
    }

    virtual void SetUpcomingGates(const std::vector<std::shared_ptr<Circuits::IOperation<double>>> & /*gates*/)
    {
    }
 
    // Prepared interactions used by initial-layout optimization. An active
    // routing plan may retain its complete indexed window during execution;
    // this is not a remaining-operations/suffix API. Inactive mid-execution
    // views are empty. Empty means the caller should retain its circuit view.
    virtual const std::vector<std::shared_ptr<Circuits::IOperation<double>>> &GetUpcomingRoutingOperations() const
    {
        static const std::vector<std::shared_ptr<Circuits::IOperation<double>>> empty;
        return empty;
    }

    virtual long long int GetGatesCounter() const
    {
        return 0;
    }

    virtual void SetGatesCounter(long long int /*counter*/)
    {
    }

    virtual void IncrementGatesCounter()
    {
    }
 
    // the following four functions are also for MPS swaps optimizations, might be
    // removed in the future
 
    virtual double getGrowthFactorSwap() const
    {
        return 0.;
    }
 
    virtual double getGrowthFactorGate() const
    {
        return 0.;
    }
 
    virtual void setGrowthFactorSwap(double factor)
    {
    }
 
    virtual void setGrowthFactorGate(double factor)
    {
    }

    virtual void Configure(const char *key, const char *value) = 0;

    virtual std::string GetConfiguration(const char *key) const = 0;

    virtual size_t AllocateQubits(size_t num_qubits) = 0;

    virtual size_t GetNumberOfQubits() const = 0;

    virtual void Clear() = 0;

    virtual size_t Measure(const Types::qubits_vector &qubits) = 0;

    virtual std::vector<bool> MeasureMany(const Types::qubits_vector &qubits) = 0;

    virtual void ApplyReset(const Types::qubits_vector &qubits) = 0;

    virtual bool SupportsQuantumChannels() const
    {
        return false;
    }

    virtual void ApplyQuantumChannel(const Types::qubits_vector &targets, const QuantumChannel &channel)
    {
        (void)targets;
        (void)channel;
        throw std::runtime_error("This simulator does not support exact quantum-channel evolution");
    }

    virtual double ProbabilityBits(const std::vector<bool> &)
    {
        throw std::runtime_error("This simulator does not support ProbabilityBits");
    }
 
    virtual std::complex<double> DensityMatrixElementBits(const std::vector<bool> &, const std::vector<bool> &) const
    {
        throw std::runtime_error("This simulator does not support DensityMatrixElementBits");
    }
 
    virtual Eigen::MatrixXcd GetDensityMatrix(bool = true) const
    {
        throw std::runtime_error("This simulator does not support GetDensityMatrix");
    }
 
    virtual std::complex<double> ExpectationValueComplex(const std::string &, bool = true) const
    {
        throw std::runtime_error("This simulator does not support ExpectationValueComplex");
    }
 
    // Same normalization and validation as the scalar complex query. Backends
    // may share contractions; other backends retain their scalar behavior.
    virtual std::vector<std::complex<double>> ExpectationValuesComplex(const std::vector<std::string> &paulis, bool normalized = true) const
    {
        std::vector<std::complex<double>> values;
        values.reserve(paulis.size());
        for (const auto &pauli : paulis)
            values.push_back(ExpectationValueComplex(pauli, normalized));
        return values;
    }
 
    virtual void ApplyOperator(const Types::qubits_vector &, const Eigen::MatrixXcd &, bool = false)
    {
        throw std::runtime_error("This simulator does not support ApplyOperator");
    }
 
    virtual void MoveAtBeginningOfChain(const Types::qubits_vector &)
    {
        throw std::runtime_error("This simulator does not support MoveAtBeginningOfChain");
    }

    virtual std::complex<double> DensityMatrixTrace() const
    {
        throw std::runtime_error("This simulator does not expose a density-matrix trace");
    }
 
    virtual double DensityMatrixPurity() const
    {
        throw std::runtime_error("This simulator does not expose density-matrix purity");
    }
 
    virtual std::complex<double> DensityMatrixTraceOfSquare() const
    {
        throw std::runtime_error("This simulator does not expose Tr(rho^2)");
    }
 
    virtual std::complex<double> DensityMatrixOverlap(const IState &) const
    {
        throw std::runtime_error("This simulator does not support density-matrix overlap");
    }
 
    virtual double DensityMatrixHermiticityResidual() const
    {
        throw std::runtime_error("This simulator does not expose a Hermiticity residual");
    }
 
    virtual bool IsDensityMatrixHermitian(double = 1e-10) const
    {
        throw std::runtime_error("This simulator does not expose a Hermiticity test");
    }
 
    virtual Eigen::MatrixXcd PartialTrace(const Types::qubits_vector &) const
    {
        throw std::runtime_error("This simulator does not support partial trace");
    }
 
    virtual double FidelityWithStatevector(const Eigen::VectorXcd &) const
    {
        throw std::runtime_error("This simulator does not support mixed-state fidelity");
    }
 
    virtual void RestoreDensityMatrixTrace()
    {
        throw std::runtime_error("This simulator cannot restore density-matrix trace");
    }
 
    virtual void HermitizeDensityMatrix()
    {
        throw std::runtime_error("This simulator cannot hermitize its density matrix");
    }

    virtual void Trim()
    {
        throw std::runtime_error("This simulator is not a matrix-product state or operator");
    }

    virtual void ReCanonicalize()
    {
        throw std::runtime_error("This simulator is not a matrix-product state or operator");
    }

    void ApplyKrausChannel(const Types::qubits_vector &targets, const QuantumChannel::KrausOperators &krausOperators)
    {
        ApplyQuantumChannel(targets, QuantumChannel(krausOperators));
    }

    void ApplyChannel(const Types::qubits_vector &targets, const QuantumChannel::KrausOperators &krausOperators)
    {
        ApplyKrausChannel(targets, krausOperators);
    }

    void ApplyPauliChannel(const Types::qubits_vector &targets, const std::vector<double> &probabilities)
    {
        ApplyQuantumChannel(targets, QuantumChannel::Pauli(probabilities));
    }

    void ApplyPauliChannel(Types::qubit_t qubit, double px, double py, double pz)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::Pauli(px, py, pz));
    }
 
    void ApplyBitFlipNoise(Types::qubit_t qubit, double probability)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::BitFlip(probability));
    }
 
    void ApplyBitPhaseFlipNoise(Types::qubit_t qubit, double probability)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::BitPhaseFlip(probability));
    }
 
    void ApplyPhaseFlipNoise(Types::qubit_t qubit, double probability)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::PhaseFlip(probability));
    }

    void ApplyDephasingNoise(Types::qubit_t qubit, double probability)
    {
        ApplyPhaseFlipNoise(qubit, probability);
    }

    void ApplyDepolarizingNoise(Types::qubit_t qubit, double errorProbability)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::Depolarizing(errorProbability));
    }

    void ApplyDepolarizingMixingNoise(Types::qubit_t qubit, double mixingProbability)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::DepolarizingMixing(mixingProbability));
    }
 
    void ApplyAmplitudeDamping(Types::qubit_t qubit, double gamma)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::AmplitudeDamping(gamma));
    }

    void ApplyT1Relaxation(Types::qubit_t qubit, double gamma)
    {
        ApplyAmplitudeDamping(qubit, gamma);
    }

    void ApplyT1RelaxationFromTime(Types::qubit_t qubit, double duration, double t1)
    {
        if (!std::isfinite(duration) || duration < 0.0)
            throw std::invalid_argument("T1-relaxation duration must be finite and nonnegative");
        if (std::isnan(t1) || t1 <= 0.0)
            throw std::invalid_argument("T1 must be positive (infinity is allowed)");
        const double gamma = std::isinf(t1) ? 0.0 : -std::expm1(-duration / t1);
        ApplyAmplitudeDamping(qubit, gamma);
    }

    void ApplyPhaseDamping(Types::qubit_t qubit, double gamma)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::PhaseDamping(gamma));
    }

    void ApplyPhaseDampingFromTime(Types::qubit_t qubit, double duration, double tPhi)
    {
        if (!std::isfinite(duration) || duration < 0.0)
            throw std::invalid_argument("Phase-damping duration must be finite and nonnegative");
        if (std::isnan(tPhi) || tPhi <= 0.0)
            throw std::invalid_argument("T_phi must be positive (infinity is allowed)");
        const double gamma = std::isinf(tPhi) ? 0.0 : -std::expm1(-2.0 * duration / tPhi);
        ApplyPhaseDamping(qubit, gamma);
    }
 
    void ApplyGeneralizedAmplitudeDamping(Types::qubit_t qubit, double gamma, double excitedStatePopulation)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::GeneralizedAmplitudeDamping(gamma, excitedStatePopulation));
    }
 
    void ApplyThermalRelaxation(Types::qubit_t qubit, double duration, double t1, double t2, double excitedStatePopulation = 0.0)
    {
        ApplyQuantumChannel({qubit}, QuantumChannel::ThermalRelaxation(duration, t1, t2, excitedStatePopulation));
    }
 
    void ApplyCorrelatedPhaseFlipNoise(Types::qubit_t qubit0, Types::qubit_t qubit1, double probability)
    {
        ApplyQuantumChannel({qubit0, qubit1}, QuantumChannel::CorrelatedPhaseFlip(probability));
    }
 
    void ApplyCorrelatedPhaseFlipNoise(Types::qubit_t qubit0, Types::qubit_t qubit1, double probability, double correlation)
    {
        ApplyQuantumChannel({qubit0, qubit1}, QuantumChannel::CorrelatedPhaseFlip(probability, correlation));
    }
 
    void ApplyTwoQubitDepolarizingNoise(Types::qubit_t qubit0, Types::qubit_t qubit1, double errorProbability)
    {
        ApplyQuantumChannel({qubit0, qubit1}, QuantumChannel::TwoQubitDepolarizing(errorProbability));
    }
 
    void ApplyTwoQubitDepolarizingMixingNoise(Types::qubit_t qubit0, Types::qubit_t qubit1, double mixingProbability)
    {
        ApplyQuantumChannel({qubit0, qubit1}, QuantumChannel::TwoQubitDepolarizingMixing(mixingProbability));
    }

    virtual double Probability(Types::qubit_t outcome) = 0;

    virtual std::complex<double> Amplitude(Types::qubit_t outcome) = 0;
 
    // Full pure-state readout, indexed by logical qubit (q0 is the low bit).
    // Backends with a bulk reader override the scalar fallback.
    virtual std::vector<std::complex<double>> GetStateVector()
    {
        if (GetSimulationType() == SimulationType::kDensityMatrix || GetSimulationType() == SimulationType::kMatrixProductOperator)
            throw std::runtime_error("Mixed states have no unique statevector");
        if (GetNumberOfQubits() == 0)
            return {};
        std::vector<std::complex<double>> result(TensorQueries::DenseSize(GetNumberOfQubits()));
        for (size_t i = 0; i < result.size(); ++i)
            result[i] = Amplitude(i);
        return result;
    }
 
    // <psi|O_(k-1)...O_0|psi>. Each matrix is 2x2; repeated targets are allowed.
    // This is a read-only query, including for non-Hermitian/non-unitary
    // operators.
    virtual std::complex<double> ExpectationValueOperators(const Types::qubits_vector &, const std::vector<Eigen::MatrixXcd> &)
    {
        throw std::runtime_error("This simulator does not support ExpectationValueOperators");
    }

    virtual std::complex<double> ProjectOnZero() = 0;

    virtual std::vector<double> AllProbabilities() = 0;

    virtual std::vector<double> Probabilities(const Types::qubits_vector &qubits) = 0;

    virtual std::unordered_map<Types::qubit_t, Types::qubit_t> SampleCounts(const Types::qubits_vector &qubits, size_t shots = 1000) = 0;

    virtual std::unordered_map<std::vector<bool>, Types::qubit_t> SampleCountsMany(const Types::qubits_vector &qubits, size_t shots = 1000) = 0;

    virtual double ExpectationValue(const std::string &pauliString) = 0;
 
    // Backends may share contractions across this read-only batch. The default
    // retains the existing behavior for every other simulator.
    virtual std::vector<double> ExpectationValues(const std::vector<std::string> &paulis)
    {
        std::vector<double> values;
        values.reserve(paulis.size());
        for (const auto &pauli : paulis)
            values.push_back(ExpectationValue(pauli));
        return values;
    }

    void RegisterObserver(const std::shared_ptr<ISimulatorObserver> &observer)
    {
        observers.insert(observer);
    }

    void UnregisterObserver(const std::shared_ptr<ISimulatorObserver> &observer)
    {
        observers.erase(observer);
    }

    void ClearObservers()
    {
        observers.clear();
    }

    virtual SimulatorType GetType() const = 0;

    virtual SimulationType GetSimulationType() const = 0;

    virtual void Flush() = 0;

    virtual void SaveStateToInternalDestructive() = 0;

    virtual void RestoreInternalDestructiveSavedState() = 0;

    virtual void SaveState() = 0;

    virtual void RestoreState() = 0;

    virtual std::complex<double> AmplitudeRaw(Types::qubit_t outcome) = 0;

    virtual void SetMultithreading(bool multithreading = true) = 0;

    virtual bool GetMultithreading() const = 0;

    virtual bool IsQcsim() const = 0;

    virtual Types::qubit_t MeasureNoCollapse() = 0;

    virtual std::vector<bool> MeasureNoCollapseMany() = 0;

    virtual size_t GetCurrentMaxBondDimension() const
    {
        return 0;
    }
 
    // Internal telemetry must not synchronize a pending gate after every call.
    virtual size_t GetExecutedMaxBondDimension() const
    {
        return GetCurrentMaxBondDimension();
    }
 
    virtual const std::unordered_map<std::string, std::string> &GetConfigMap() const = 0;
 
  protected:
    void SeedAuxiliaryRng(uint64_t seed)
    {
        auxRng.seed(DeriveSeed(seed, kAuxRngStream));
        auxUniform.reset();
    }

    void DontNotify()
    {
        notifyObservers = false;
    }

    void Notify()
    {
        notifyObservers = true;
    }

    void NotifyObservers(const Types::qubits_vector &affectedQubits)
    {
        if (!notifyObservers)
            return;
 
        for (auto &observer : observers)
        {
            observer->Update(affectedQubits);
        }
    }
 
  protected:
    bool NotificationsEnabled() const noexcept
    {
        return notifyObservers;
    }
 
  private:
    static constexpr uint64_t kAuxRngStream = 0x52454144ULL; // 'READ'
 
    std::unordered_set<std::shared_ptr<ISimulatorObserver>> observers; 
    bool notifyObservers = true;                                       
    std::mt19937_64 auxRng{std::random_device{}()};                    
    std::uniform_real_distribution<double> auxUniform{0.0, 1.0};
};
 
} // namespace Simulators
 
#endif // !_SIMULATOR_STATE_H_