GpuState.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _GPUSTATE_H_
#define _GPUSTATE_H_
 
#ifdef __linux__
 
#ifdef INCLUDED_BY_FACTORY
 
#include "MPSDummySimulator.h"
 
#include <cstdint>
#include <iomanip>
#include <limits>
#include <sstream>
#include <utility>
 
#include "Configuration.h"
 
namespace Simulators {
// TODO: Maybe use the pimpl idiom
// https://en.cppreference.com/w/cpp/language/pimpl to hide the implementation
// for good but during development this should be good enough
namespace Private {

class GpuState : public ISimulator {
 public:
  void Initialize() override {
    if (nrQubits) {
      auto initializationLock = GpuLibrary::GetInstance()->LockInitialization();
      const int gpuDevice = configuration.IsSet("gpu_device")
          ? Configuration::ParseGpuDevice(configuration.GetConfiguration("gpu_device"))
          : SimulatorsFactory::ResolveGpuDevice();
      configuration.SetConfiguration("gpu_device", std::to_string(gpuDevice));
      if (!SimulatorsFactory::GetGpuLibrary(gpuDevice))
        throw std::runtime_error("GpuState::Initialize: Unable to initialize GPU device " +
                                 std::to_string(gpuDevice));
      if (simulationType == SimulationType::kStatevector) {
        state = SimulatorsFactory::CreateGpuLibStateVectorSim(gpuDevice);
        if (state) {
          // ensure the config settings are applied, they need to be applied
          // after the simulator is created
          for (const auto& [key, value] : configuration.GetConfigMap())
            if (key != "method") Configure(key.c_str(), value.c_str());
 
          const bool res = state->Create(nrQubits);
          if (!res)
            throw std::runtime_error(
                "GpuState::Initialize: Failed to create "
                "and initialize the statevector state.");
        } else
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the statevector state.");
      } else if (simulationType == SimulationType::kDensityMatrix) {
        densityMatrix = SimulatorsFactory::CreateGpuDensityMatrix(gpuDevice);
        if (!densityMatrix)
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the density matrix state.");
        // Precision must be selected before native density-matrix storage is
        // allocated.
        for (const auto& [key, value] : configuration.GetConfigMap())
          if (key != "method") Configure(key.c_str(), value.c_str());
        if (!densityMatrix->Create(nrQubits))
          throw std::runtime_error(
              "GpuState::Initialize: Failed to initialize the density matrix state.");
      } else if (simulationType == SimulationType::kMatrixProductOperator) {
        mpo = SimulatorsFactory::CreateGpuMPO(gpuDevice);
        if (!mpo)
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the matrix product "
              "operator state.");
        mpo->SetCallbackContext((void*)this);
        curMaxBondDim = 1;
        mpo->SetBondDimensionsCallback(&GpuState::BondDimCallback);
 
        // Precision and truncation controls must be selected before native
        // storage is allocated.
        for (const auto& [key, value] : configuration.GetConfigMap())
          if (key != "method") Configure(key.c_str(), value.c_str());
        if (!mpo->Create(nrQubits))
          throw std::runtime_error(
              "GpuState::Initialize: Failed to initialize the matrix product "
              "operator state.");
        // default is true
        if (!useOptimalMeetingPosition)
          mpo->SetUseOptimalMeetingPosition(false);
      } else if (simulationType == SimulationType::kMatrixProductState) {
        mps = SimulatorsFactory::CreateGpuLibMPSSim(gpuDevice);
        if (mps) {
          mps->SetCallbackContext((void*)this);
          curMaxBondDim = 1;
          mps->SetBondDimensionsCallback(&GpuState::BondDimCallback);
 
          // ensure the config settings are applied, they need to be applied
          // after the simulator is created but before the state is created
          for (const auto& [key, value] : configuration.GetConfigMap())
            if (key != "method") Configure(key.c_str(), value.c_str());
 
          const bool res = mps->Create(nrQubits);
          if (!res)
            throw std::runtime_error(
                "GpuState::Initialize: Failed to create "
                "and initialize the MPS state.");
        } else
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the MPS state.");
        // default is true
        if (!useOptimalMeetingPosition)
          mps->SetUseOptimalMeetingPosition(false);
      } else if (simulationType == SimulationType::kTensorNetwork) {
        tn = SimulatorsFactory::CreateGpuLibTensorNetSim(gpuDevice);
        if (tn) {
          // ensure the config settings are applied, they need to be applied
          // after the simulator is created but before the state is created
          for (const auto& [key, value] : configuration.GetConfigMap())
            if (key != "method") Configure(key.c_str(), value.c_str());
 
          const bool res = tn->Create(nrQubits);
          if (!res)
            throw std::runtime_error(
                "GpuState::Initialize: Failed to create "
                "and initialize the tensor network state.");
        } else
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the tensor network "
              "state.");
      } else if (simulationType == SimulationType::kPauliPropagator) {
        pp = SimulatorsFactory::CreateGpuPauliPropagatorSimulatorUnique(gpuDevice);
        if (pp) {
          // ensure the config settings are applied, they need to be applied
          // after the simulator is created but before the state is created
          for (const auto& [key, value] : configuration.GetConfigMap())
            if (key != "method") Configure(key.c_str(), value.c_str());
 
          const bool res = pp->CreateSimulator(nrQubits);
          if (!res)
            throw std::runtime_error(
                "GpuState::Initialize: Failed to create "
                "and initialize the Pauli propagator state.");
 
          pp->SetWillUseSampling(true);  // TODO: check setting
          if (!pp->AllocateMemory(0.9))
            throw std::runtime_error(
                "GpuState::Initialize: Failed to allocate memory for the "
                "Pauli propagator state.");
        } else
          throw std::runtime_error(
              "GpuState::Initialize: Failed to create the Pauli propagator "
              "state.");
      } else
        throw std::runtime_error(
            "GpuState::Initialize: Invalid simulation "
            "type for initializing the state.");
      if (GetGpuDevice() != gpuDevice)
        throw std::runtime_error("GpuState::Initialize: GPU plugin did not confirm the requested device; update the GPU library");
    }
  }

  void InitializeState(size_t num_qubits,
                       std::vector<std::complex<double>> &amplitudes) override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    if (simulationType != SimulationType::kStatevector &&
        simulationType != SimulationType::kDensityMatrix &&
        simulationType != SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    const bool created =
        simulationType == SimulationType::kDensityMatrix
            ? densityMatrix->CreateWithState(
                  nrQubits,
                  reinterpret_cast<const double *>(amplitudes.data()))
            : simulationType == SimulationType::kMatrixProductOperator
                  ? mpo->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()))
                  : state->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()));
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeState: Failed to initialize the state.");
  }

#ifndef NO_QISKIT_AER
  void InitializeState(size_t num_qubits,
                       AER::Vector<std::complex<double>> &amplitudes) override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    if (simulationType != SimulationType::kStatevector &&
        simulationType != SimulationType::kDensityMatrix &&
        simulationType != SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    const bool created =
        simulationType == SimulationType::kDensityMatrix
            ? densityMatrix->CreateWithState(
                  nrQubits,
                  reinterpret_cast<const double *>(amplitudes.data()))
            : simulationType == SimulationType::kMatrixProductOperator
                  ? mpo->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()))
                  : state->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()));
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeState: Failed to initialize the state.");
  }
#endif

  void InitializeState(size_t num_qubits,
                       Eigen::VectorXcd &amplitudes) override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    if (simulationType != SimulationType::kStatevector &&
        simulationType != SimulationType::kDensityMatrix &&
        simulationType != SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    const bool created =
        simulationType == SimulationType::kDensityMatrix
            ? densityMatrix->CreateWithState(
                  nrQubits,
                  reinterpret_cast<const double *>(amplitudes.data()))
            : simulationType == SimulationType::kMatrixProductOperator
                  ? mpo->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()))
                  : state->CreateWithState(
                        nrQubits,
                        reinterpret_cast<const double *>(amplitudes.data()));
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeState: Failed to initialize the state.");
  }

  void InitializeToBasisState(size_t num_qubits,
                              Types::qubit_t basisState) override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    bool created = true;
    if (simulationType == SimulationType::kDensityMatrix)
      created = densityMatrix->CreateWithBasisState(
          nrQubits, static_cast<unsigned long long>(basisState));
    else if (simulationType == SimulationType::kMatrixProductOperator)
      created = mpo->CreateWithBasisState(
          nrQubits, static_cast<unsigned long long>(basisState));
    else if (simulationType == SimulationType::kMatrixProductState)
      created = mps->CreateWithBasisState(
          nrQubits, static_cast<unsigned long long>(basisState));
    else
      for (size_t q = 0; q < num_qubits; ++q)
        if ((basisState >> q) & 1ULL) ApplyX(static_cast<Types::qubit_t>(q));
 
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeToBasisState: Failed to initialize the "
          "state.");
  }

  void InitializeToBasisState(size_t num_qubits,
                              const std::vector<bool> &basisState) override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    bool created = true;
    if (simulationType == SimulationType::kMatrixProductOperator ||
        simulationType == SimulationType::kMatrixProductState) {
      std::vector<unsigned char> stateBits(num_qubits, 0);
      for (size_t q = 0; q < num_qubits && q < basisState.size(); ++q)
        stateBits[q] = basisState[q] ? 1 : 0;
      created = simulationType == SimulationType::kMatrixProductOperator
                    ? mpo->CreateWithBasisStateBits(nrQubits, stateBits)
                    : mps->CreateWithBasisStateBits(nrQubits, stateBits);
    } else
      for (size_t q = 0; q < num_qubits && q < basisState.size(); ++q)
        if (basisState[q]) ApplyX(static_cast<Types::qubit_t>(q));
 
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeToBasisState: Failed to initialize the "
          "state.");
  }

  void InitializeToMixtureOfBasisStates(
      size_t num_qubits,
      const std::vector<std::pair<Types::qubit_t, double>> &mixture)
      override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    if (simulationType != SimulationType::kDensityMatrix &&
        simulationType != SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::InitializeToMixtureOfBasisStates: Invalid simulation "
          "type for initializing to a mixture of basis states.");
 
    std::vector<std::pair<unsigned long long, double>> converted;
    converted.reserve(mixture.size());
    for (const auto &[basisState, weight] : mixture)
      converted.emplace_back(static_cast<unsigned long long>(basisState),
                             weight);
 
    const bool created =
        simulationType == SimulationType::kDensityMatrix
            ? densityMatrix->CreateWithMixtureOfBasisStates(nrQubits,
                                                            converted)
            : mpo->CreateWithMixtureOfBasisStates(nrQubits, converted);
 
    if (!created)
      throw std::runtime_error(
          "GpuState::InitializeToMixtureOfBasisStates: Failed to initialize "
          "the state.");
  }

  void InitializeToMixtureOfBasisStates(
      size_t num_qubits,
      const std::vector<std::pair<std::vector<bool>, double>> &mixture)
      override {
    if (num_qubits == 0) return;
    Clear();
    nrQubits = num_qubits;
    Initialize();
 
    if (simulationType != SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::InitializeToMixtureOfBasisStates: Invalid simulation "
          "type for initializing to a mixture of basis states.");
 
    std::vector<double> weights;
    weights.reserve(mixture.size());
    std::vector<unsigned char> stateBitsFlat;
    stateBitsFlat.reserve(mixture.size() * num_qubits);
    for (const auto &[basisState, weight] : mixture) {
      weights.push_back(weight);
      for (size_t q = 0; q < num_qubits; ++q)
        stateBitsFlat.push_back(
            (q < basisState.size() && basisState[q]) ? 1 : 0);
    }
 
    if (!mpo->CreateWithMixtureOfBasisStatesBits(nrQubits, stateBitsFlat,
                                                 weights))
      throw std::runtime_error(
          "GpuState::InitializeToMixtureOfBasisStates: Failed to initialize "
          "the state.");
  }

  void Reset() override {
    if (state)
      state->Reset();
    else if (densityMatrix)
      densityMatrix->Reset();
    else if (mpo) {
      mpo->Reset();
      curMaxBondDim = 1;
    } else if (mps) {
      mps->Reset();
      curMaxBondDim = 1;
    } else if (tn)
      tn->Reset();
    else if (pp)
      pp->ClearOperators();
 
    upcomingGateIndex = 0;
  }

  bool SupportsMPSSwapOptimization() const override { return true; }

  void SetInitialQubitsMap(
      const std::vector<long long int> &initialMap) override {
    if (mps || mpo) {
      if (mps) mps->SetInitialQubitsMap(initialMap);
      else mpo->SetInitialQubitsMap(initialMap);
      if (!dummySim || dummySim->getNrQubits() != initialMap.size()) {
        dummySim =
            std::make_unique<Simulators::MPSDummySimulator>(initialMap.size());
        dummySim->SetMaxBondDimension(
            configuration.GetConfigurationAsInt(MaxBondDimensionConfigKey()));
      }
      dummySim->setGrowthFactorGate(growthFactorGate);
      dummySim->setGrowthFactorSwap(growthFactorSwap);
      dummySim->SetInitialQubitsMap(initialMap);
    }
  }
 
  void SetUseOptimalMeetingPosition(bool enable) override {
    useOptimalMeetingPosition = enable;
    if (mps || mpo) {
      if (mps) mps->SetUseOptimalMeetingPosition(enable);
      else mpo->SetUseOptimalMeetingPosition(enable);
 
      if (enable) {
        // Register an observer that advances the gate index
        ClearObservers();  // for now we only have this observer, so this should
                           // be fine
        gateCounterObserver =
            std::make_shared<GateCounterObserver>(upcomingGateIndex);
        RegisterObserver(gateCounterObserver);
 
        // Set up a meeting position callback that uses MPSDummySimulator
        // for lookahead evaluation with actual bond dimensions
        // the callback is called only for two qubits gates and only if
        // executing them would require a swap
        if (mps)
          mps->SetMeetingPositionCallback(&GpuState::FindBestMeetingPosition);
        else
          mpo->SetMeetingPositionCallback(&GpuState::FindBestMeetingPosition);
      }
    }
  }
 
  void SetLookaheadDepth(int depth) override {
    lookaheadDepth = depth;
    if (depth > 0 && !useOptimalMeetingPosition) {
      if (mps) mps->SetUseOptimalMeetingPosition(true);
      else if (mpo) mpo->SetUseOptimalMeetingPosition(true);
    }
  }
 
  void SetLookaheadDepthWithHeuristic(int depth) override {
    lookaheadDepthWithHeuristic = depth;
    if (lookaheadDepth < depth) SetLookaheadDepth(depth);
  }
 
  void SetUpcomingGates(
      const std::vector<std::shared_ptr<Circuits::IOperation<double>>> &gates)
      override {
    upcomingGates = gates;
    upcomingGateIndex = 0;
 
    if (!mps && !mpo) return;
 
    // Register an observer that advances the gate index
    ClearObservers();  // for now we only have this observer, so this should be
                       // fine
    gateCounterObserver =
        std::make_shared<GateCounterObserver>(upcomingGateIndex);
    RegisterObserver(gateCounterObserver);
 
    // Set up a meeting position callback that uses MPSDummySimulator
    // for lookahead evaluation with actual bond dimensions
    // the callback is called only for two qubits gates and only if executing
    // them would require a swap
    if (mps)
      mps->SetMeetingPositionCallback(&GpuState::FindBestMeetingPosition);
    else
      mpo->SetMeetingPositionCallback(&GpuState::FindBestMeetingPosition);
  }

  long long int GetGatesCounter() const override { return upcomingGateIndex; }

  void SetGatesCounter(long long int counter) override {
    upcomingGateIndex = counter;
  }

  void IncrementGatesCounter() override { ++upcomingGateIndex; }
 
  double getGrowthFactorSwap() const override { return growthFactorSwap; }
  double getGrowthFactorGate() const override { return growthFactorGate; }
 
  void setGrowthFactorSwap(double factor) override {
    growthFactorSwap = factor;
    if (dummySim) dummySim->setGrowthFactorSwap(factor);
  }
 
  void setGrowthFactorGate(double factor) override {
    growthFactorGate = factor;
    if (dummySim) dummySim->setGrowthFactorGate(factor);
  }

  void Configure(const char *key, const char *value) override {
    if (!key || !value) return;
    if (std::string("gpu_device") == key) {
      const int device = Configuration::ParseGpuDevice(value);
      if ((state || densityMatrix || mpo || mps || tn || pp) &&
          device != Configuration::ParseGpuDevice(configuration.GetConfiguration(key)))
        throw std::invalid_argument("gpu_device cannot change after initialization; clear the simulator first");
      configuration.SetConfiguration(key, std::to_string(device));
      return;
    }
    const auto svdGroup = Configuration::GpuSvdSettingGroup(key);
    if (!svdGroup.empty()) {
      const bool enabled = Configuration::ParseGpuSvdFlag(value);
      const bool gesvd = svdGroup == key;
      const char algorithm = std::string(key).back();
      const auto apply = [algorithm, enabled](auto& backend) {
        if (!backend) return true; // Applied after native object creation.
        if (algorithm == 'j') return backend->SetGesvdJ(enabled);
        if (algorithm == 'p') return backend->SetGesvdP(enabled);
        return backend->SetGesvdR(enabled);
      };
      const auto applyGesvd = [enabled](auto& backend) {
        // GESVD has no native flag. A false entry is an inactive selector,
        // so replaying it must not clear a subsequently selected algorithm.
        if (!enabled || !backend) return true;
        return backend->SetGesvdJ(false) && backend->SetGesvdP(false) &&
               backend->SetGesvdR(false);
      };
      bool applied = true;
      if (svdGroup == "matrix_product_state_use_gesvd")
        applied = gesvd ? applyGesvd(mps) : apply(mps);
      else if (svdGroup == "matrix_product_operator_use_gesvd")
        applied = gesvd ? applyGesvd(mpo) : apply(mpo);
      else if (svdGroup == "tensor_network_use_gesvd")
        applied = gesvd ? applyGesvd(tn) : apply(tn);
      if (!applied)
        throw std::runtime_error(std::string("GPU library cannot apply ") + key +
                                 "; an updated GPU plugin may be required");
      configuration.SetConfiguration(key, value);
      return;
    }
    if (std::string("method") == key) {
      if (std::string("statevector") == value)
        simulationType = SimulationType::kStatevector;
      else if (std::string("matrix_product_state") == value)
        simulationType = SimulationType::kMatrixProductState;
      else if (std::string("density_matrix") == value)
        simulationType = SimulationType::kDensityMatrix;
      else if (std::string("matrix_product_operator") == value)
        simulationType = SimulationType::kMatrixProductOperator;
      else if (std::string("tensor_network") == value)
        simulationType = SimulationType::kTensorNetwork;
      else if (std::string("pauli_propagator") == value)
        simulationType = SimulationType::kPauliPropagator;
 
      // Match the GPU library's default cap.
      if ((simulationType == SimulationType::kMatrixProductState ||
           simulationType == SimulationType::kMatrixProductOperator) &&
          !configuration.IsSet("matrix_product_state_max_bond_dimension"))
        configuration.SetConfiguration(
            "matrix_product_state_max_bond_dimension", "128");
      if (simulationType == SimulationType::kMatrixProductOperator) {
        const auto bondDimension =
            configuration.GetConfiguration(MaxBondDimensionConfigKey());
        configuration.SetConfiguration(
            "matrix_product_state_max_bond_dimension", bondDimension);
        configuration.SetConfiguration(
            "matrix_product_operator_max_bond_dimension", bondDimension);
      }
    }
 
    if (std::string("use_double_precision") == key && densityMatrix &&
        densityMatrix->IsCreated())
      throw std::runtime_error(
          "GpuState::Configure: Density-matrix precision must be configured "
          "before initialization.");
 
    if (std::string("use_double_precision") == key && mpo && mpo->IsCreated())
      throw std::runtime_error(
          "GpuState::Configure: Matrix-product-operator precision must be "
          "configured before initialization.");
 
    if (!configuration.WasApplied(key, value))
        configuration.SetConfiguration(key, value);
 
    if (std::string("seed") == key) {
      const uint64_t seed = std::stoull(value);
      nextSeedStream = 0;
      if (state) state->SetSeed(seed);
      if (densityMatrix) densityMatrix->SetSeed(seed);
      if (mpo) mpo->SetSeed(seed);
      if (mps) mps->SetSeed(seed);
      if (tn) tn->SetSeed(seed);
      if (pp) pp->SetSeed(seed);
      return;
    }
 
    if (std::string("matrix_product_state_truncation_threshold") == key ||
        std::string("matrix_product_operator_truncation_threshold") == key) {
      // SetCutoff() sets the numeric threshold value. How that number is interpreted --
      // relative to the largest singular value at each bond/split, or as a cumulative
      // discarded-weight budget -- is a separate, independently configurable setting; see
      // matrix_product_state_truncation_mode / matrix_product_operator_truncation_mode
      // below. All three GPU backends default to discarded-weight (matching Qiskit Aer's
      // and ITensor's convention) unless relative_max is explicitly requested -- see
      // TruncationMode in maestro-gpu-simulators' lib/truncationmode.hpp and its use in
      // mpsimpl.cu/mpo.cu/tensornet.cu.
      const double singularValueThreshold = std::stod(value);
      if (singularValueThreshold > 0.) {
        if (mps) mps->SetCutoff(singularValueThreshold);
        if (tn) tn->SetCutoff(singularValueThreshold);
        if (mpo) mpo->SetCutoff(singularValueThreshold);
      }
    } else if (std::string("matrix_product_state_truncation_mode") == key ||
               std::string("matrix_product_operator_truncation_mode") == key) {
      // "relative_max" -> TruncationMode::RelativeToMax (0), "discarded_weight" ->
      // TruncationMode::DiscardedWeight (1, the default -- see lib/truncationmode.hpp in
      // maestro-gpu-simulators).
      int truncationMode = -1;
      if (std::string("relative_max") == value)
        truncationMode = 0;
      else if (std::string("discarded_weight") == value)
        truncationMode = 1;
      if (truncationMode >= 0) {
        if (mps) mps->SetTruncationMode(truncationMode);
        if (tn) tn->SetTruncationMode(truncationMode);
        if (mpo) mpo->SetTruncationMode(truncationMode);
      }
    } else if (std::string("matrix_product_state_max_bond_dimension") == key ||
               std::string("matrix_product_operator_max_bond_dimension") ==
                   key) {
      const long long int chi = std::stoi(value);
      if (simulationType == SimulationType::kMatrixProductOperator) {
        // Both names address the same MPO setting; the latest write wins.
        const std::string bondDimension(value);
        for (const char* alias : {"matrix_product_state_max_bond_dimension",
                                  "matrix_product_operator_max_bond_dimension"})
          if (!configuration.WasApplied(alias, bondDimension))
            configuration.SetConfiguration(alias, bondDimension);
      }
      if (chi > 0) {
        if (mps) mps->SetMaxExtent(chi);
        if (tn) tn->SetMaxExtent(chi);
        if (mpo) mpo->SetMaxExtent(chi);
        if (dummySim) dummySim->SetMaxBondDimension(chi);
      }
 
    } else if (std::string("matrix_product_operator_kraus_completeness_check") == key) {
      int mode = -1;
      if (std::string(value) == "ignore") mode = 0;
      else if (std::string(value) == "warn") mode = 1;
      else if (std::string(value) == "strict") mode = 2;
      if (mpo && mode >= 0 && !mpo->SetKrausCompletenessCheck(mode))
        throw std::runtime_error("Invalid GPU MPO Kraus completeness mode");
    } else if (std::string("use_double_precision") == key) {
      const bool useDoublePrecision =
          (std::string("1") == value || std::string("true") == value);
      if (mps) mps->SetDataType(useDoublePrecision);
      if (tn) tn->SetDataType(useDoublePrecision);
      if (state) state->SetDataType(useDoublePrecision);
      if (densityMatrix) densityMatrix->SetDataType(useDoublePrecision);
      if (mpo) mpo->SetDataType(useDoublePrecision);
    }
    
    if (pp) {
      if (std::string("pauli_propagator_coefficient_threshold") == key) {
        const double coefficientThreshold = std::stod(value);
        pp->SetCoefficientTruncationCutoff(coefficientThreshold);
      } else if (std::string("pauli_propagator_pauli_weight_threshold") ==
                 key) {
        const double pauliWeightThreshold = std::stod(value);
        pp->SetWeightTruncationCutoff(pauliWeightThreshold);
      } else if (std::string("pauli_propagator_steps_between_trims") == key) {
        const int stepsBetweenTrims = std::stoi(value);
        pp->SetNumGatesBetweenTruncations(stepsBetweenTrims);
      } else if (std::string("pauli_propagator_num_gates_between_deduplications") ==
                 key) {
        const int numGatesBetweenDeduplications = std::stoi(value);
        pp->SetNumGatesBetweenDeduplications(numGatesBetweenDeduplications);
      }
    }
  }

  std::string GetConfiguration(const char *key) const override {
    if (!key) return {};
    const auto svdGroup = Configuration::GpuSvdSettingGroup(key);
    if (!svdGroup.empty()) {
      if (svdGroup == key) {
        const auto readGesvd = [](const auto& backend) {
          return !backend->GetGesvdJ() && !backend->GetGesvdP() &&
                 !backend->GetGesvdR();
        };
        if (svdGroup == "matrix_product_state_use_gesvd" && mps)
          return readGesvd(mps) ? "true" : "false";
        if (svdGroup == "matrix_product_operator_use_gesvd" && mpo)
          return readGesvd(mpo) ? "true" : "false";
        if (svdGroup == "tensor_network_use_gesvd" && tn)
          return readGesvd(tn) ? "true" : "false";
      }
      const char algorithm = std::string(key).back();
      const auto read = [algorithm](const auto& backend) {
        if (algorithm == 'j') return backend->GetGesvdJ();
        if (algorithm == 'p') return backend->GetGesvdP();
        return backend->GetGesvdR();
      };
      if (svdGroup == "matrix_product_state_use_gesvd" && mps)
        return read(mps) ? "true" : "false";
      if (svdGroup == "matrix_product_operator_use_gesvd" && mpo)
        return read(mpo) ? "true" : "false";
      if (svdGroup == "tensor_network_use_gesvd" && tn)
        return read(tn) ? "true" : "false";
    }
    if (std::string("method") == key) {
      switch (simulationType) {
        case SimulationType::kStatevector:
          return "statevector";
        case SimulationType::kMatrixProductState:
          return "matrix_product_state";
        case SimulationType::kDensityMatrix:
          return "density_matrix";
        case SimulationType::kMatrixProductOperator:
          return "matrix_product_operator";
        case SimulationType::kTensorNetwork:
          return "tensor_network";
        case SimulationType::kPauliPropagator:
          return "pauli_propagator";
        default:
          return "other";
      }
    }
 
    return configuration.GetConfiguration(key);
  }

  size_t AllocateQubits(size_t num_qubits) override {
    if ((simulationType == SimulationType::kStatevector && state) ||
        (simulationType == SimulationType::kDensityMatrix && densityMatrix) ||
        (simulationType == SimulationType::kMatrixProductOperator && mpo) ||
        (simulationType == SimulationType::kMatrixProductState && mps) ||
        (simulationType == SimulationType::kPauliPropagator && pp))
      return 0;
 
    const size_t oldNrQubits = nrQubits;
    nrQubits += num_qubits;
 
    return oldNrQubits;
  }

  size_t GetNumberOfQubits() const override { return nrQubits; }

  void Clear() override {
    state = nullptr;
    densityMatrix = nullptr;
    mpo = nullptr;
    mps = nullptr;
    tn = nullptr;
    pp = nullptr;
    nrQubits = 0;
    dummySim = nullptr;
    upcomingGateIndex = 0;
    upcomingGates.clear();
  }

  size_t Measure(const Types::qubits_vector &qubits) override {
    // TODO: this is inefficient, maybe implement it better in gpu sim
    // for now it has the possibility of measuring a qubits interval, but not a
    // list of qubits
    if (qubits.size() > sizeof(size_t) * 8)
      std::cerr
          << "Warning: The number of qubits to measure is larger than the "
             "number of bits in the size_t type, the outcome will be undefined"
          << std::endl;
 
    size_t res = 0;
    size_t mask = 1ULL;
 
    DontNotify();
    if (simulationType == SimulationType::kStatevector) {
      // TODO: measure all qubits in one shot?
      for (size_t qubit : qubits) {
        if (state->MeasureQubitCollapse(static_cast<int>(qubit))) res |= mask;
        mask <<= 1;
      }
    } else if (simulationType == SimulationType::kDensityMatrix) {
      for (size_t qubit : qubits) {
        if (densityMatrix->Measure(static_cast<unsigned int>(qubit))) res |= mask;
        mask <<= 1;
      }
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      for (size_t qubit : qubits) {
        if (mpo->Measure(static_cast<unsigned int>(qubit))) res |= mask;
        mask <<= 1;
      }
    } else if (simulationType == SimulationType::kMatrixProductState) {
      // TODO: measure all qubits in one shot?
      for (size_t qubit : qubits) {
        if (mps->Measure(static_cast<unsigned int>(qubit))) res |= mask;
        mask <<= 1;
      }
    } else if (simulationType == SimulationType::kTensorNetwork) {
      // TODO: measure all qubits in one shot?
      for (size_t qubit : qubits) {
        if (tn->Measure(static_cast<unsigned int>(qubit))) res |= mask;
        mask <<= 1;
      }
    } else if (simulationType == SimulationType::kPauliPropagator) {
      // TODO: measure all qubits in one shot?
      for (size_t qubit : qubits) {
        if (pp->MeasureQubit(static_cast<int>(qubit))) res |= mask;
        mask <<= 1;
      }
    }
 
    Notify();
    NotifyObservers(qubits);
 
    return res;
  }

  std::vector<bool> MeasureMany(const Types::qubits_vector &qubits) override {
    std::vector<bool> res(qubits.size(), false);
 
    DontNotify();
    if (simulationType == SimulationType::kStatevector) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = state->MeasureQubitCollapse(static_cast<int>(qubits[i]));
    } else if (simulationType == SimulationType::kDensityMatrix) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = densityMatrix->Measure(qubits[i]);
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = mpo->Measure(static_cast<unsigned int>(qubits[i]));
    } else if (simulationType == SimulationType::kMatrixProductState) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = mps->Measure(static_cast<unsigned int>(qubits[i]));
    } else if (simulationType == SimulationType::kTensorNetwork) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = tn->Measure(static_cast<unsigned int>(qubits[i]));
    } else if (simulationType == SimulationType::kPauliPropagator) {
      for (size_t i = 0; i < qubits.size(); ++i)
        res[i] = pp->MeasureQubit(static_cast<int>(qubits[i]));
    }
    Notify();
    NotifyObservers(qubits);
 
    return res;
  }

  void ApplyReset(const Types::qubits_vector &qubits) override {
    DontNotify();
    if (simulationType == SimulationType::kStatevector) {
      for (size_t qubit : qubits)
        if (state->MeasureQubitCollapse(static_cast<int>(qubit)))
          state->ApplyX(static_cast<int>(qubit));
    } else if (simulationType == SimulationType::kDensityMatrix) {
      for (size_t qubit : qubits) densityMatrix->ApplyReset(qubit);
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      for (size_t qubit : qubits)
        mpo->ApplyReset(static_cast<int>(qubit));
    } else if (simulationType == SimulationType::kMatrixProductState) {
      for (size_t qubit : qubits)
        if (mps->Measure(static_cast<unsigned int>(qubit)))
          mps->ApplyX(static_cast<unsigned int>(qubit));
    } else if (simulationType == SimulationType::kTensorNetwork) {
      for (size_t qubit : qubits)
        if (tn->Measure(static_cast<unsigned int>(qubit)))
          tn->ApplyX(static_cast<unsigned int>(qubit));
    } else if (simulationType == SimulationType::kPauliPropagator) {
      for (size_t qubit : qubits)
        if (pp->MeasureQubit(static_cast<int>(qubit)))
          pp->ApplyX(static_cast<int>(qubit));
    }
 
    Notify();
    NotifyObservers(qubits);
  }
 
  bool SupportsQuantumChannels() const override {
    return simulationType == SimulationType::kDensityMatrix ||
           simulationType == SimulationType::kMatrixProductOperator;
  }
 
  void ApplyQuantumChannel(const Types::qubits_vector& targets,
                           const QuantumChannel& channel) override {
    if (!densityMatrix && !mpo)
      throw std::runtime_error(
          "GPU quantum channels require an initialized density matrix or "
          "matrix product operator");
    if (targets.size() != channel.GetNumberOfQubits() || targets.empty() ||
        targets.size() > 2)
      throw std::invalid_argument(
          "GPU density matrices and matrix product operators support one- "
          "and two-qubit local channels");
    std::vector<int> gpuTargets;
    gpuTargets.reserve(targets.size());
    for (auto target : targets) {
      if (target >= nrQubits ||
          std::find(gpuTargets.begin(), gpuTargets.end(), target) !=
              gpuTargets.end())
        throw std::invalid_argument("Invalid GPU quantum-channel target");
      gpuTargets.push_back(static_cast<int>(target));
    }
    const auto& kraus = channel.GetKrausOperators();
    std::vector<double> interleaved;
    interleaved.reserve(kraus.size() * kraus.front().size() * 2);
    for (const auto& op : kraus)
      for (Eigen::Index i = 0; i < op.size(); ++i) {
        interleaved.push_back(op.data()[i].real());
        interleaved.push_back(op.data()[i].imag());
      }
    const bool applied =
        densityMatrix
            ? densityMatrix->ApplyKraus(gpuTargets, kraus.size(),
                                        interleaved.data())
            : mpo->ApplyKraus(gpuTargets, kraus.size(), interleaved.data());
    if (!applied)
      throw std::runtime_error(
          "GPU density-matrix/matrix-product-operator channel application "
          "failed");
    NotifyObservers(targets);
  }
 
  std::complex<double> DensityMatrixTrace() const override {
    if (densityMatrix) return densityMatrix->Trace();
    if (mpo) return mpo->Trace();
    throw std::runtime_error("GPU mixed-state diagnostics require density_matrix or matrix_product_operator");
  }
  double DensityMatrixPurity() const override {
    if (densityMatrix) return densityMatrix->Purity();
    if (mpo) return mpo->Purity();
    throw std::runtime_error("GPU mixed-state diagnostics require density_matrix or matrix_product_operator");
  }
  std::complex<double> DensityMatrixTraceOfSquare() const override {
    if (mpo) return mpo->TraceOfSquare();
    if (densityMatrix) {
      const double tr = densityMatrix->Trace();
      return densityMatrix->Purity() * tr * tr;
    }
    throw std::runtime_error("GPU mixed-state diagnostics require density_matrix or matrix_product_operator");
  }
  std::complex<double> DensityMatrixOverlap(const IState &other) const override {
    const auto *rhs = dynamic_cast<const GpuState *>(&other);
    if (!rhs) throw std::invalid_argument("Density-matrix overlap requires matching GPU backends");
    if (densityMatrix && rhs->densityMatrix)
      return densityMatrix->HilbertSchmidtOverlap(*rhs->densityMatrix);
    if (mpo && rhs->mpo) return mpo->HilbertSchmidtOverlap(*rhs->mpo);
    throw std::invalid_argument("Density-matrix overlap requires two density matrices or two MPOs");
  }
  double DensityMatrixHermiticityResidual() const override {
    if (mpo) return mpo->HermiticityResidual();
    if (densityMatrix) return densityMatrix->IsHermitian() ? 0. : std::numeric_limits<double>::infinity();
    throw std::runtime_error("GPU mixed-state diagnostics require density_matrix or matrix_product_operator");
  }
  bool IsDensityMatrixHermitian(double eps = 1e-10) const override {
    if (densityMatrix) return densityMatrix->IsHermitian(eps);
    if (mpo) return mpo->IsHermitian(eps);
    throw std::runtime_error("GPU mixed-state diagnostics require density_matrix or matrix_product_operator");
  }
  Eigen::MatrixXcd PartialTrace(const Types::qubits_vector &qubits) const override {
    std::vector<int> keep(qubits.begin(), qubits.end());
    const auto values = densityMatrix ? densityMatrix->PartialTrace(keep) :
                        mpo ? mpo->PartialTrace(keep) : throw std::runtime_error("GPU partial trace requires a mixed-state backend");
    const Eigen::Index dim = static_cast<Eigen::Index>(size_t{1} << keep.size());
    return Eigen::Map<const Eigen::MatrixXcd>(values.data(), dim, dim);
  }
  double FidelityWithStatevector(const Eigen::VectorXcd &psi) const override {
    std::vector<double> raw(2 * static_cast<size_t>(psi.size()));
    for (Eigen::Index i = 0; i < psi.size(); ++i) { raw[2*i] = psi[i].real(); raw[2*i+1] = psi[i].imag(); }
    if (densityMatrix) return densityMatrix->FidelityWithStatevector(raw.data());
    if (mpo) return mpo->FidelityWithStatevector(raw.data());
    throw std::runtime_error("GPU mixed-state fidelity requires density_matrix or matrix_product_operator");
  }
  void RestoreDensityMatrixTrace() override {
    if (!mpo) throw std::runtime_error("Trace restoration is only available for GPU MPO");
    mpo->RestoreTrace();
  }
  void HermitizeDensityMatrix() override {
    if (!mpo) throw std::runtime_error("Hermitization is only available for GPU MPO");
    mpo->Hermitize();
  }
  void TrimMatrixProductOperator() override {
    if (!mpo) throw std::runtime_error("GPU MPO is not initialized");
    mpo->Trim();
  }
  void ReCanonicalizeMatrixProductOperator() override {
    if (!mpo) throw std::runtime_error("GPU MPO is not initialized");
    mpo->ReCanonicalize();
  }

  double Probability(Types::qubit_t outcome) override {
    if (simulationType == SimulationType::kStatevector)
      return state->BasisStateProbability(outcome);
    else if (simulationType == SimulationType::kDensityMatrix)
      return densityMatrix->Probability(outcome);
    else if (simulationType == SimulationType::kMatrixProductOperator)
      return mpo->Probability(outcome);
    else if (simulationType == SimulationType::kMatrixProductState ||
             simulationType == SimulationType::kTensorNetwork) {
      const auto ampl = Amplitude(outcome);
      return std::norm(ampl);
    } else if (simulationType == SimulationType::kPauliPropagator) {
      return pp->Probability(outcome);
    }
 
    return 0.0;
  }

  std::complex<double> Amplitude(Types::qubit_t outcome) override {
    double real = 0.0;
    double imag = 0.0;
 
    if (simulationType == SimulationType::kStatevector)
      state->Amplitude(outcome, &real, &imag);
    else if (simulationType == SimulationType::kDensityMatrix)
      throw std::runtime_error(
          "GpuState::Amplitude: Amplitudes are not defined for density matrices.");
    else if (simulationType == SimulationType::kMatrixProductOperator)
      throw std::runtime_error(
          "GpuState::Amplitude: Amplitudes are not defined for matrix "
          "product operators.");
    else if (simulationType == SimulationType::kMatrixProductState ||
             simulationType == SimulationType::kTensorNetwork) {
      std::vector<long int> fixedValues(nrQubits);
      for (size_t i = 0; i < nrQubits; ++i)
        fixedValues[i] = (outcome & (1ULL << i)) ? 1 : 0;
      if (simulationType == SimulationType::kMatrixProductState)
        mps->Amplitude(nrQubits, fixedValues.data(), &real, &imag);
      else if (simulationType == SimulationType::kTensorNetwork)
        tn->Amplitude(nrQubits, fixedValues.data(), &real, &imag);
    } else if (simulationType == SimulationType::kPauliPropagator) {
      // Pauli propagator does not support amplitude calculation
      throw std::runtime_error(
          "GpuState::Amplitude: Invalid simulation type for amplitude "
          "calculation.");
    }
 
    return std::complex<double>(real, imag);
  }

  std::complex<double> ProjectOnZero() override {
    if (simulationType == SimulationType::kMatrixProductState)
      return mps->ProjectOnZero();
 
    return Amplitude(0);
  }

  std::vector<double> AllProbabilities() override {
    if (nrQubits == 0) return {};
    const size_t numStates = 1ULL << nrQubits;
    std::vector<double> result(numStates);
 
    if (simulationType == SimulationType::kStatevector)
      state->AllProbabilities(result.data());
    else if (simulationType == SimulationType::kDensityMatrix)
      densityMatrix->AllProbabilities(result.data());
    else if (simulationType == SimulationType::kMatrixProductOperator)
      mpo->AllProbabilities(result.data());
    else if (simulationType == SimulationType::kMatrixProductState ||
             simulationType == SimulationType::kTensorNetwork) {
      // this is very slow, it should be used only for tests!
      for (Types::qubit_t i = 0; i < (Types::qubit_t)numStates; ++i) {
        const auto val = Amplitude(i);
        result[i] = std::norm(std::complex<double>(val.real(), val.imag()));
      }
    } else if (simulationType == SimulationType::kPauliPropagator) {
      // this is very slow, it should be used only for tests!
      for (Types::qubit_t i = 0; i < (Types::qubit_t)numStates; ++i) {
        result[i] = pp->Probability(i);
      }
    }
 
    return result;
  }

  std::vector<double> Probabilities(
      const Types::qubits_vector &qubits) override {
    std::vector<double> result(qubits.size());
 
    if (simulationType == SimulationType::kStatevector) {
      for (size_t i = 0; i < qubits.size(); ++i)
        result[i] = state->BasisStateProbability(qubits[i]);
    } else if (simulationType == SimulationType::kDensityMatrix) {
      for (size_t i = 0; i < qubits.size(); ++i)
        result[i] = densityMatrix->Probability(qubits[i]);
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      for (size_t i = 0; i < qubits.size(); ++i)
        result[i] = mpo->Probability(qubits[i]);
    } else if (simulationType == SimulationType::kMatrixProductState ||
               simulationType == SimulationType::kTensorNetwork) {
      for (size_t i = 0; i < qubits.size(); ++i) {
        const auto ampl = Amplitude(qubits[i]);
        result[i] = std::norm(ampl);
      }
    } else if (simulationType == SimulationType::kPauliPropagator) {
      for (size_t i = 0; i < qubits.size(); ++i)
        result[i] = pp->Probability(qubits[i]);
    }
 
    return result;
  }

  std::unordered_map<Types::qubit_t, Types::qubit_t> SampleCounts(
      const Types::qubits_vector &qubits, size_t shots = 1000) override {
    if (qubits.empty() || shots == 0) return {};
 
    if (qubits.size() > sizeof(Types::qubit_t) * 8)
      std::cerr
          << "Warning: The number of qubits to measure is larger than the "
             "number of bits in the Types::qubit_t type, the outcome will be "
             "undefined"
          << std::endl;
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> result;
 
    DontNotify();
 
    if (simulationType == SimulationType::kStatevector) {
      std::vector<long int> samples(shots);
      state->SampleAll(shots, samples.data());
 
      for (auto outcome : samples) {
        // qubits might not be in order, translate the outcome to the correct
        // order
        Types::qubit_t translatedOutcome = 0;
        Types::qubit_t mask = 1ULL;
        for (size_t i = 0; i < qubits.size(); ++i) {
          if (outcome & (1ULL << qubits[i])) translatedOutcome |= mask;
          mask <<= 1;
        }
        ++result[translatedOutcome];
      }
    } else if (simulationType == SimulationType::kDensityMatrix) {
      std::vector<long int> samples(shots);
      if (!densityMatrix->SampleAll(shots, samples.data())) {
        Notify();
        throw std::runtime_error(
            "GpuState::SampleCounts: Density-matrix sampling failed.");
      }
      for (auto outcome : samples) {
        Types::qubit_t translatedOutcome = 0;
        for (size_t i = 0; i < qubits.size(); ++i)
          if (outcome & (1ULL << qubits[i])) translatedOutcome |= 1ULL << i;
        ++result[translatedOutcome];
      }
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      std::vector<long int> samples(shots);
      if (!mpo->SampleAll(shots, samples.data())) {
        Notify();
        throw std::runtime_error(
            "GpuState::SampleCounts: Matrix-product-operator sampling "
            "failed.");
      }
      for (auto outcome : samples) {
        Types::qubit_t translatedOutcome = 0;
        for (size_t i = 0; i < qubits.size(); ++i)
          if (outcome & (1ULL << qubits[i])) translatedOutcome |= 1ULL << i;
        ++result[translatedOutcome];
      }
    } else if (simulationType == SimulationType::kMatrixProductState) {
      std::unordered_map<std::vector<bool>, int64_t> *map =
          mps->GetMapForSample();
 
      std::vector<unsigned int> qubitsIndices(qubits.begin(), qubits.end());
 
      mps->Sample(shots, qubitsIndices.size(), qubitsIndices.data(), map);
      const auto positions = SampleBitPositions(qubits);
 
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) {
        Types::qubit_t outcome = 0;
        Types::qubit_t mask = 1ULL;
        for (Types::qubit_t q = 0; q < qubits.size(); ++q) {
          if (meas[positions[q]]) outcome |= mask;
          mask <<= 1;
        }
 
        result[outcome] += cnt;
      }
 
      mps->FreeMapForSample(map);
    } else if (simulationType == SimulationType::kTensorNetwork) {
      std::unordered_map<std::vector<bool>, int64_t> *map =
          tn->GetMapForSample();
      std::vector<unsigned int> qubitsIndices(qubits.begin(), qubits.end());
      tn->Sample(shots, qubitsIndices.size(), qubitsIndices.data(), map);
      const auto positions = SampleBitPositions(qubits);
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) {
        Types::qubit_t outcome = 0;
        Types::qubit_t mask = 1ULL;
        for (Types::qubit_t q = 0; q < qubits.size(); ++q) {
          if (meas[positions[q]]) outcome |= mask;
          mask <<= 1;
        }
        result[outcome] += cnt;
      }
      tn->FreeMapForSample(map);
    } else if (simulationType == SimulationType::kPauliPropagator) {
      std::vector<int> qb(qubits.begin(), qubits.end());
      for (size_t shot = 0; shot < shots; ++shot) {
        size_t meas = 0;
        auto res = pp->SampleQubits(qb);
        for (size_t i = 0; i < qubits.size(); ++i) {
          if (res[i]) meas |= (1ULL << i);
        }
        ++result[meas];
      }
    }
 
    Notify();
    NotifyObservers(qubits);
 
    return result;
  }

  std::unordered_map<std::vector<bool>, Types::qubit_t> SampleCountsMany(
      const Types::qubits_vector &qubits, size_t shots = 1000) override {
    if (qubits.empty() || shots == 0) return {};
 
    std::unordered_map<std::vector<bool>, Types::qubit_t> result;
 
    DontNotify();
 
    if (simulationType == SimulationType::kStatevector) {
      std::vector<long int> samples(shots);
      state->SampleAll(shots, samples.data());
 
      std::vector<bool> outcomeVec(qubits.size());
      for (auto outcome : samples) {
        for (size_t i = 0; i < qubits.size(); ++i)
          outcomeVec[i] = ((outcome >> qubits[i]) & 1) == 1;
        ++result[outcomeVec];
      }
    } else if (simulationType == SimulationType::kDensityMatrix) {
      std::vector<long int> samples(shots);
      if (!densityMatrix->SampleAll(shots, samples.data())) {
        Notify();
        throw std::runtime_error(
            "GpuState::SampleCountsMany: Density-matrix sampling failed.");
      }
      std::vector<bool> outcomeVec(qubits.size());
      for (auto outcome : samples) {
        for (size_t i = 0; i < qubits.size(); ++i)
          outcomeVec[i] = ((outcome >> qubits[i]) & 1) != 0;
        ++result[outcomeVec];
      }
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      std::vector<long int> samples(shots);
      if (!mpo->SampleAll(shots, samples.data())) {
        Notify();
        throw std::runtime_error(
            "GpuState::SampleCountsMany: Matrix-product-operator sampling "
            "failed.");
      }
      std::vector<bool> outcomeVec(qubits.size());
      for (auto outcome : samples) {
        for (size_t i = 0; i < qubits.size(); ++i)
          outcomeVec[i] = ((outcome >> qubits[i]) & 1) != 0;
        ++result[outcomeVec];
      }
    } else if (simulationType == SimulationType::kMatrixProductState) {
      std::unordered_map<std::vector<bool>, int64_t> *map =
          mps->GetMapForSample();
 
      std::vector<unsigned int> qubitsIndices(qubits.begin(), qubits.end());
      mps->Sample(shots, qubitsIndices.size(), qubitsIndices.data(), map);
      const auto positions = SampleBitPositions(qubits);
 
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) {
        std::vector<bool> ordered(qubits.size());
        for (size_t q = 0; q < qubits.size(); ++q) ordered[q] = meas[positions[q]];
        result[ordered] += cnt;
      }
 
      mps->FreeMapForSample(map);
    } else if (simulationType == SimulationType::kTensorNetwork) {
      std::unordered_map<std::vector<bool>, int64_t> *map =
          tn->GetMapForSample();
      std::vector<unsigned int> qubitsIndices(qubits.begin(), qubits.end());
      tn->Sample(shots, qubitsIndices.size(), qubitsIndices.data(), map);
      const auto positions = SampleBitPositions(qubits);
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) {
        std::vector<bool> ordered(qubits.size());
        for (size_t q = 0; q < qubits.size(); ++q) ordered[q] = meas[positions[q]];
        result[ordered] += cnt;
      }
      tn->FreeMapForSample(map);
    } else if (simulationType == SimulationType::kPauliPropagator) {
      std::vector<int> qb(qubits.begin(), qubits.end());
      for (size_t shot = 0; shot < shots; ++shot) {
        const auto res = pp->SampleQubits(qb);
        ++result[res];
      }
    }
 
    Notify();
    NotifyObservers(qubits);
 
    return result;
  }

  double ExpectationValue(const std::string &pauliString) override {
    double result = 0.0;
 
    if (simulationType == SimulationType::kStatevector)
      result = state->ExpectationValue(pauliString);
    else if (simulationType == SimulationType::kDensityMatrix)
      result = densityMatrix->ExpectationValue(pauliString);
    else if (simulationType == SimulationType::kMatrixProductOperator)
      result = mpo->ExpectationValue(pauliString);
    else if (simulationType == SimulationType::kMatrixProductState)
      result = mps->ExpectationValue(pauliString);
    else if (simulationType == SimulationType::kTensorNetwork)
      result = tn->ExpectationValue(pauliString);
    else if (simulationType == SimulationType::kPauliPropagator)
      result = pp->ExpectationValue(pauliString);
    else
      throw std::runtime_error(
          "GpuState::ExpectationValue: Invalid simulation type for expectation "
          "value calculation.");
 
    return result;
  }

  int GetGpuDevice() const override {
    if (state) return state->GetGpuDevice();
    if (densityMatrix) return densityMatrix->GetGpuDevice();
    if (mpo) return mpo->GetGpuDevice();
    if (mps) return mps->GetGpuDevice();
    if (tn) return tn->GetGpuDevice();
    if (pp) return pp->GetGpuDevice();
    return -1;
  }
 
  SimulatorType GetType() const override { return SimulatorType::kGpuSim; }

  SimulationType GetSimulationType() const override { return simulationType; }

  void Flush() override {}

  void SaveStateToInternalDestructive() override {
    if (simulationType == SimulationType::kStatevector)
      state->SaveStateDestructive();
    else if (simulationType == SimulationType::kPauliPropagator)
      return;
    else
      throw std::runtime_error(
          "GpuState::SaveStateToInternalDestructive: Invalid simulation type "
          "for saving the state destructively.");
  }

  void RestoreInternalDestructiveSavedState() override {
    if (simulationType == SimulationType::kStatevector)
      state->RestoreStateFreeSaved();
    else if (simulationType == SimulationType::kPauliPropagator)
      return;
    else
      throw std::runtime_error(
          "GpuState::RestoreInternalDestructiveSavedState: Invalid simulation "
          "type for restoring the state destructively.");
  }

  void SaveState() override {
    if (simulationType == SimulationType::kStatevector)
      state->SaveState();
    else if (simulationType == SimulationType::kDensityMatrix)
      densityMatrix->SaveState();
    else if (simulationType == SimulationType::kMatrixProductOperator)
      mpo->SaveState();
    else if (simulationType == SimulationType::kMatrixProductState)
      mps->SaveState();
    else if (simulationType == SimulationType::kTensorNetwork)
      tn->SaveState();
    else if (simulationType == SimulationType::kPauliPropagator)
      pp->SaveState();
  }

  void RestoreState() override {
    if (simulationType == SimulationType::kStatevector)
      state->RestoreStateNoFreeSaved();
    else if (simulationType == SimulationType::kDensityMatrix)
      densityMatrix->RestoreState();
    else if (simulationType == SimulationType::kMatrixProductOperator)
      mpo->RestoreState();
    else if (simulationType == SimulationType::kMatrixProductState)
      mps->RestoreState();
    else if (simulationType == SimulationType::kTensorNetwork)
      tn->RestoreState();
    else if (simulationType == SimulationType::kPauliPropagator)
      pp->RestoreState();
  }

  std::complex<double> AmplitudeRaw(Types::qubit_t outcome) override {
    return Amplitude(outcome);
  }

  void SetMultithreading(bool multithreading = true) override {
    // don't do anything here, the multithreading is always enabled
  }

  bool GetMultithreading() const override { return true; }

  bool IsQcsim() const override { return false; }

  Types::qubit_t MeasureNoCollapse() override {
    if (simulationType == SimulationType::kStatevector)
      return state->MeasureAllQubitsNoCollapse();
    else if (simulationType == SimulationType::kDensityMatrix) {
      std::vector<long int> samples(1);
      if (!densityMatrix->SampleAll(1, samples.data()))
        throw std::runtime_error(
            "GpuState::MeasureNoCollapse: Density-matrix sampling failed.");
      return static_cast<Types::qubit_t>(samples.front());
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      std::vector<long int> samples(1);
      if (!mpo->SampleAll(1, samples.data()))
        throw std::runtime_error(
            "GpuState::MeasureNoCollapse: Matrix-product-operator sampling "
            "failed.");
      return static_cast<Types::qubit_t>(samples.front());
    } else if (simulationType == SimulationType::kMatrixProductState ||
             simulationType == SimulationType::kTensorNetwork ||
             simulationType == SimulationType::kPauliPropagator) {
      if (nrQubits > sizeof(Types::qubit_t) * 8)
        std::cerr
            << "Warning: The number of qubits to measure is larger than the "
               "number of bits in the Types::qubit_t type, the outcome will be "
               "undefined"
            << std::endl;
 
      Types::qubits_vector fixedValues(nrQubits);
      std::iota(fixedValues.begin(), fixedValues.end(), 0);
      const auto res = SampleCounts(fixedValues, 1);
      if (res.empty()) return 0;
      return res.begin()
          ->first;  // return the first outcome, as it is the only one
    }
 
    throw std::runtime_error(
        "GpuState::MeasureNoCollapse: Invalid simulation type for measuring "
        "all the qubits without collapsing the state.");
 
    return 0;
  }

  std::vector<bool> MeasureNoCollapseMany() override {
    if (simulationType == SimulationType::kStatevector) {
      const auto meas = state->MeasureAllQubitsNoCollapse();
      std::vector<bool> result(nrQubits, false);
      for (size_t i = 0; i < nrQubits; ++i) result[i] = ((meas >> i) & 1) == 1;
      return result;
    } else if (simulationType == SimulationType::kDensityMatrix) {
      std::vector<long int> samples(1);
      if (!densityMatrix->SampleAll(1, samples.data()))
        throw std::runtime_error(
            "GpuState::MeasureNoCollapseMany: Density-matrix sampling failed.");
      std::vector<bool> result(nrQubits, false);
      for (size_t i = 0; i < nrQubits; ++i)
        result[i] = ((samples.front() >> i) & 1) != 0;
      return result;
    } else if (simulationType == SimulationType::kMatrixProductOperator) {
      std::vector<long int> samples(1);
      if (!mpo->SampleAll(1, samples.data()))
        throw std::runtime_error(
            "GpuState::MeasureNoCollapseMany: Matrix-product-operator "
            "sampling failed.");
      std::vector<bool> result(nrQubits, false);
      for (size_t i = 0; i < nrQubits; ++i)
        result[i] = ((samples.front() >> i) & 1) != 0;
      return result;
    } else if (simulationType == SimulationType::kMatrixProductState ||
               simulationType == SimulationType::kTensorNetwork ||
               simulationType == SimulationType::kPauliPropagator) {
      Types::qubits_vector fixedValues(nrQubits);
      std::iota(fixedValues.begin(), fixedValues.end(), 0);
      const auto res = SampleCountsMany(fixedValues, 1);
      if (res.empty()) return std::vector<bool>(nrQubits, false);
      return res.begin()
          ->first;  // return the first outcome, as it is the only one
    }
 
    throw std::runtime_error(
        "GpuState::MeasureNoCollapseMany: Invalid simulation type for "
        "measuring "
        "all the qubits without collapsing the state.");
 
    return std::vector<bool>(nrQubits, false);
  }

  size_t GetCurrentMaxBondDimension() const override { return curMaxBondDim; }
 
  
  const Configuration& GetConfiguration() const { return configuration; }
 
    const std::unordered_map<std::string, std::string>& GetConfigMap()
      const override {
    return configuration.GetConfigMap();
  }
 
 protected:
  // MPSSample/TNSample consume a set and return bits in ascending logical
  // qubit order. Maestro's sampling interface preserves the caller's order.
  static std::vector<size_t> SampleBitPositions(const Types::qubits_vector& qubits) {
    auto sorted = qubits;
    std::sort(sorted.begin(), sorted.end());
    sorted.erase(std::unique(sorted.begin(), sorted.end()), sorted.end());
    std::vector<size_t> positions;
    positions.reserve(qubits.size());
    for (const auto q : qubits)
      positions.push_back(std::lower_bound(sorted.begin(), sorted.end(), q) - sorted.begin());
    return positions;
  }
 
  // MPO accepts both names; Configure keeps their values synchronized.
  // The fallback also handles configuration before the method is selected.
  const char* MaxBondDimensionConfigKey() const {
    return simulationType == SimulationType::kMatrixProductOperator &&
                   configuration.IsSet(
                       "matrix_product_operator_max_bond_dimension")
               ? "matrix_product_operator_max_bond_dimension"
               : "matrix_product_state_max_bond_dimension";
  }
 
  static int64_t FindBestMeetingPosition(void* thisPtr, const int64_t* bondDims) {
    GpuState* self = static_cast<GpuState*>(thisPtr);
 
    return self->FindBestMeetingPositionFunc(bondDims);
  };
 
  int64_t FindBestMeetingPositionFunc(const int64_t* bondDims)
  {
    const size_t nQ = GetNumberOfQubits();
 
    if (lookaheadDepth <= 0 || lookaheadDepth == std::numeric_limits<int>::max()) 
       return -1;
 
    if (!dummySim || dummySim->getNrQubits() != nQ) {
      dummySim = std::make_unique<Simulators::MPSDummySimulator>(nQ);
      dummySim->SetMaxBondDimension(
          configuration.GetConfigurationAsInt(MaxBondDimensionConfigKey()));
      dummySim->setGrowthFactorGate(growthFactorGate);
      dummySim->setGrowthFactorSwap(growthFactorSwap);
    }
 
    dummySim->setTotalSwappingCost(0);
 
    // Convert actual bond dims to doubles
    std::vector<double> bondDimsD(bondDims, bondDims + nrQubits - 1);
    dummySim->SetCurrentBondDimensions(bondDimsD);
 
      // display bond dimensions for debugging
#ifdef LOG_CALLBACK_INFO
    std::cerr << "Bond dimensions before swapping and applying the gate:";
    for (size_t i = 0; i < nrQubits - 1; ++i) {
      std::cerr << bondDims[i] << " ";
    }
    std::cerr << std::endl;
#endif
 
    if (upcomingGates.size() <= static_cast<size_t>(upcomingGateIndex)) {
      return -1;  // will fallback
    }
 
    const auto &op = upcomingGates[upcomingGateIndex];
    const auto qbits = op->AffectedQubits();
 
    if (qbits.size() != 2) {
      std::cerr << "Error: Meeting position callback called for a gate "
                   "that does not have exactly 2 qubits."
                << std::endl;
 
      return -1;  // will fallback
    }
 
#ifdef LOG_CALLBACK_INFO
    const auto& qmap = dummySim->getQubitsMap();
 
    std::cerr << "Applying 2-qubit gate on physical qubits " << qmap[qbits[0]]
              << " and " << qmap[qbits[1]] << std::endl;
 
    std::cerr << "Finding best meeting position for upcoming gates starting at index "
              << upcomingGateIndex << " with lookahead depth " << lookaheadDepth
              << " and heuristic depth " << lookaheadDepthWithHeuristic
              << std::endl;
 
    std::cerr << "Affected qubits: ";
    for (const auto& q : qbits) std::cerr << q << " ";
    std::cerr << std::endl;
#endif
 
    double bestCost = std::numeric_limits<double>::infinity();
    int64_t res = dummySim->FindBestMeetingPosition(
        upcomingGates, upcomingGateIndex, lookaheadDepth,
        lookaheadDepthWithHeuristic, 0, bestCost);
 
#ifdef LOG_CALLBACK_INFO
    std::cerr << "Swapping the two qubits on position: " << res << " and "
              << (res + 1) << std::endl;
#endif
 
    dummySim->SwapQubitsToPosition(qbits[0], qbits[1], res);
    dummySim->ApplyGate(op);
 
    // display the expected bond dimensions after applying the gate for
    // debugging
 
#ifdef LOG_CALLBACK_INFO
    const auto& expectedBondDims = dummySim->getCurrentBondDimensions();
    std::cerr << "Expected bond dimensions after swapping and applying "
                 "the gate: ";
    for (size_t i = 0; i < expectedBondDims.size(); ++i) {
      std::cerr << expectedBondDims[i] << " ";
    }
    std::cerr << std::endl;
 
    std::cerr << "Best meeting position: " << res
              << " with estimated cost: " << bestCost << std::endl;
#endif
 
 
    return res;
  }
 
  static void BondDimCallback(void* thisPtr, const int64_t* bondDims) {
    GpuState* self = static_cast<GpuState*>(thisPtr);
    
    return self->BondDimCallbackFunc(bondDims);
  }
 
  void BondDimCallbackFunc(const int64_t* bondDims)
  {
    if (bondDims) {
      const size_t nQ = GetNumberOfQubits();
      for (int i = 0; i < static_cast<int>(nQ) - 1; ++i)
        if (static_cast<size_t>(bondDims[i]) > curMaxBondDim) curMaxBondDim = static_cast<size_t>(bondDims[i]);
    }
  }
 
 
  SimulationType simulationType =
      SimulationType::kStatevector; 
  uint64_t nextSeedStream = 0;
 
  std::unique_ptr<GpuLibStateVectorSim>
      state;                         
  std::unique_ptr<GpuDensityMatrix>
      densityMatrix;                 
  std::unique_ptr<GpuMPO>
      mpo; 
  std::unique_ptr<GpuLibMPSSim> mps; 
  std::unique_ptr<GpuLibTNSim> tn;   
  std::unique_ptr<GpuPauliPropagator>
      pp; 
 
  size_t nrQubits = 0; 
 
  int lookaheadDepth = 0;
  int lookaheadDepthWithHeuristic = 0;
  bool useOptimalMeetingPosition = true;
  std::vector<std::shared_ptr<Circuits::IOperation<>>> upcomingGates;
  long long int upcomingGateIndex = 0;
  double growthFactorSwap = 1.;
  double growthFactorGate = 0.65; 
 
  std::unique_ptr<Simulators::MPSDummySimulator> dummySim;
 
  // Observer that counts applied gates to track position in upcomingGates
  class GateCounterObserver : public ISimulatorObserver {
   public:
    GateCounterObserver(long long int &indexRef) : index(indexRef) {}
    void Update(const Types::qubits_vector &) override { ++index; }
 
   private:
    long long int &index;
  };
 
  std::shared_ptr<GateCounterObserver> gateCounterObserver;
  size_t curMaxBondDim = 0;
 
  Configuration configuration; 
};
 
}  // namespace Private
}  // namespace Simulators
 
#endif
#endif
#endif