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>
 
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) {
      if (simulationType == SimulationType::kStatevector) {
        state = SimulatorsFactory::CreateGpuLibStateVectorSim();
        if (state) {
          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::kMatrixProductState) {
        mps = SimulatorsFactory::CreateGpuLibMPSSim();
        if (mps) {
          if (useDoublePrecision) mps->SetDataType(true);
          if (limitEntanglement && singularValueThreshold > 0.)
            mps->SetCutoff(singularValueThreshold);
          if (limitSize && chi > 0) mps->SetMaxExtent(chi);
          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();
        if (tn) {
          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();
        if (pp) {
          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.");
    }
  }

  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)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    if (nrQubits) {
      if (simulationType == SimulationType::kStatevector) {
        state = SimulatorsFactory::CreateGpuLibStateVectorSim();
        if (state)
          state->CreateWithState(
              nrQubits, reinterpret_cast<const double *>(amplitudes.data()));
      }
    }
  }

#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)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    if (nrQubits) {
      if (simulationType == SimulationType::kStatevector) {
        state = SimulatorsFactory::CreateGpuLibStateVectorSim();
        if (state)
          state->CreateWithState(
              nrQubits, reinterpret_cast<const double *>(amplitudes.data()));
      }
    }
  }
#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)
      throw std::runtime_error(
          "GpuState::InitializeState: Invalid simulation "
          "type for initializing the state.");
 
    if (nrQubits) {
      if (simulationType == SimulationType::kStatevector) {
        state = SimulatorsFactory::CreateGpuLibStateVectorSim();
        if (state)
          state->CreateWithState(
              nrQubits, reinterpret_cast<const double *>(amplitudes.data()));
      }
    }
  }

  void Reset() override {
    if (state)
      state->Reset();
    else if (mps)
      mps->Reset();
    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) {
      mps->SetInitialQubitsMap(initialMap);
      if (!dummySim || dummySim->getNrQubits() != initialMap.size()) {
        dummySim =
            std::make_unique<Simulators::MPSDummySimulator>(initialMap.size());
        dummySim->SetMaxBondDimension(
            limitSize ? static_cast<long long int>(chi) : 0);
      }
      dummySim->setGrowthFactorGate(growthFactorGate);
      dummySim->setGrowthFactorSwap(growthFactorSwap);
      dummySim->SetInitialQubitsMap(initialMap);
    }
  }
 
  void SetUseOptimalMeetingPosition(bool enable) override {
    useOptimalMeetingPosition = enable;
    if (mps) mps->SetUseOptimalMeetingPosition(enable);
  }
 
  void SetLookaheadDepth(int depth) override {
    lookaheadDepth = depth;
    if (mps && depth > 0 && !useOptimalMeetingPosition)
      mps->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 || lookaheadDepth <= 0) 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
    mps->SetCallbackContext((void*)this);
    mps->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 (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("tensor_network") == value)
        simulationType = SimulationType::kTensorNetwork;
      else if (std::string("pauli_propagator") == value)
        simulationType = SimulationType::kPauliPropagator;
    } else if (std::string("matrix_product_state_truncation_threshold") ==
               key) {
      singularValueThreshold = std::stod(value);
      if (singularValueThreshold > 0.) {
        limitEntanglement = true;
        if (mps) mps->SetCutoff(singularValueThreshold);
        if (tn) tn->SetCutoff(singularValueThreshold);
      } else
        limitEntanglement = false;
    } else if (std::string("matrix_product_state_max_bond_dimension") == key) {
      chi = std::stoi(value);
      if (chi > 0) {
        limitSize = true;
        if (mps) mps->SetMaxExtent(chi);
        if (tn) tn->SetMaxExtent(chi);
      } else
        limitSize = false;
    } else if (std::string("use_double_precision") == key) {
      useDoublePrecision =
          (std::string("1") == value || std::string("true") == value);
    }
    // TODO: add pauli propagator configuration options
  }

  std::string GetConfiguration(const char *key) const override {
    if (std::string("method") == key) {
      switch (simulationType) {
        case SimulationType::kStatevector:
          return "statevector";
        case SimulationType::kMatrixProductState:
          return "matrix_product_state";
        case SimulationType::kTensorNetwork:
          return "tensor_network";
        case SimulationType::kPauliPropagator:
          return "pauli_propagator";
        default:
          return "other";
      }
    } else if (std::string("matrix_product_state_truncation_threshold") ==
               key) {
      if (limitEntanglement && singularValueThreshold > 0.) {
        std::ostringstream oss;
        oss << std::setprecision(std::numeric_limits<double>::max_digits10)
            << singularValueThreshold;
        return oss.str();
      }
    } else if (std::string("matrix_product_state_max_bond_dimension") == key) {
      if (limitSize && limitSize > 0) return std::to_string(chi);
    }
    // TODO: add pauli propagator configuration options
 
    return "";
  }

  size_t AllocateQubits(size_t num_qubits) override {
    if ((simulationType == SimulationType::kStatevector && state) ||
        (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;
    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::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::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::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);
  }

  double Probability(Types::qubit_t outcome) override {
    if (simulationType == SimulationType::kStatevector)
      return state->BasisStateProbability(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::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::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::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::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);
 
      // 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[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);
      // 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[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::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);
 
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) result[meas] += 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);
      // put the results in the result map
      for (const auto &[meas, cnt] : *map) result[meas] += 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::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;
  }

  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
      throw std::runtime_error(
          "GpuState::SaveStateToInternalDestructive: Invalid simulation type "
          "for saving the state destructively.");
  }

  void RestoreInternalDestructiveSavedState() override {
    if (simulationType == SimulationType::kStatevector)
      state->RestoreStateFreeSaved();
    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::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::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::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(
        "QCSimState::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::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(
        "QCSimState::MeasureNoCollapseMany: Invalid simulation type for "
        "measuring "
        "all the qubits without collapsing the state.");
 
    return std::vector<bool>(nrQubits, false);
  }
 
 protected:
  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 (!dummySim || dummySim->getNrQubits() != nQ) {
      dummySim = std::make_unique<Simulators::MPSDummySimulator>(nQ);
      dummySim->SetMaxBondDimension(
          limitSize ? static_cast<long long int>(chi) : 0);
      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);
 
    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
    }
 
    double bestCost = std::numeric_limits<double>::infinity();
    int64_t res = dummySim->FindBestMeetingPosition(
        upcomingGates, upcomingGateIndex, lookaheadDepth,
        lookaheadDepthWithHeuristic, 0, bestCost);
 
    dummySim->SwapQubitsToPosition(qbits[0], qbits[1], res);
    dummySim->ApplyGate(op);
 
    return res;
  }
 
  SimulationType simulationType =
      SimulationType::kStatevector; 
 
  std::unique_ptr<GpuLibStateVectorSim>
      state;                         
  std::unique_ptr<GpuLibMPSSim> mps; 
  std::unique_ptr<GpuLibTNSim> tn;   
  std::unique_ptr<GpuPauliPropagator>
      pp; 
 
  size_t nrQubits = 0; 
  bool limitSize = false;
  bool limitEntanglement = false;
  Eigen::Index chi = 10;               // if limitSize is true
  double singularValueThreshold = 0.;  // if limitEntanglement is true
  bool useDoublePrecision = false;     // if true, GPU MPS uses float64
 
  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;
};
 
}  // namespace Private
}  // namespace Simulators
 
#endif
#endif
#endif