SimpleDisconnectedNetwork.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _SIMPLE_NETWORK_H_
#define _SIMPLE_NETWORK_H_
 
#include "QubitRegister.h"
#include "SimpleController.h"
#include "SimpleHost.h"
#include "../Estimators/SimulatorsEstimatorInterface.h"
#include "NetworkJob.h"
 
#include "../Simulators/MPSDummySimulator.h"
 
#include "Configuration.h"
 
namespace Network {

template <typename Time = Types::time_type,
          class Controller = SimpleController<Time>>
class SimpleDisconnectedNetwork : public INetwork<Time> {
 public:
  using BaseClass = INetwork<Time>; 
  using ExecuteResults =
      typename BaseClass::ExecuteResults; 

  SimpleDisconnectedNetwork(const std::vector<Types::qubit_t> &qubits = {},
                            const std::vector<size_t> &cbits = {}) {
    configuration.SetConfiguration("use_double_precision", "0");
 
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQCSim,
         Simulators::SimulationType::kStatevector});
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQCSim,
         Simulators::SimulationType::kStabilizer});
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQCSim,
         Simulators::SimulationType::kMatrixProductState});
    // this needs more work until it becomes useful - it's typically way too
    // slow to be timed
    // simulatorsForOptimizations.insert({
    // SimulatorType::kQCSim,
    // SimulationType::kTensorNetwork });
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kCompositeQCSim,
         Simulators::SimulationType::kStatevector});
 
#ifndef NO_QISKIT_AER
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQiskitAer,
         Simulators::SimulationType::kStatevector});
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQiskitAer,
         Simulators::SimulationType::kStabilizer});
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kQiskitAer,
         Simulators::SimulationType::kMatrixProductState});
    simulatorsForOptimizations.insert(
        {Simulators::SimulatorType::kCompositeQiskitAer,
         Simulators::SimulationType::kStatevector});
#endif
 
#ifdef __linux__
    { // GPU candidates are resolved lazily when simulation is requested.
      simulatorsForOptimizations.insert(
          {Simulators::SimulatorType::kGpuSim,
           Simulators::SimulationType::kStatevector});
      simulatorsForOptimizations.insert(
          {Simulators::SimulatorType::kGpuSim,
           Simulators::SimulationType::kMatrixProductState});
      simulatorsForOptimizations.insert(
          {Simulators::SimulatorType::kGpuSim,
           Simulators::SimulationType::kTensorNetwork});
    }
#endif
 
    if (!qubits.empty()) CreateNetwork(qubits, cbits);
  }

  void CreateNetwork(const std::vector<Types::qubit_t> &qubits,
                     const std::vector<size_t> &cbits) {
    size_t qubitsOffset = 0;
    size_t cbitsOffset = 0;
 
    for (size_t i = 0; i < qubits.size(); ++i) {
      const size_t numQubits = qubits[i];
      const size_t numBits = (i < cbits.size() ? cbits[i] : 0);
      hosts.emplace_back(std::make_shared<SimpleHost<Time>>(
          i, qubitsOffset, numQubits, cbitsOffset, numBits));
      qubitsOffset += numQubits;
      cbitsOffset += numBits;
    }
 
    for (size_t i = 0; i < hosts.size(); ++i) {
      std::static_pointer_cast<SimpleHost<Time>>(hosts[i])->SetEntangledQubitId(
          qubitsOffset);
      std::static_pointer_cast<SimpleHost<Time>>(hosts[i])
          ->SetEntangledQubitMeasurementBit(cbitsOffset);
      ++qubitsOffset;
      ++cbitsOffset;
    }
 
    controller = std::make_shared<Controller>();
  }

  void Execute(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    recreateIfNeeded = false;
 
    const auto res = RepeatedExecute(circuit, 1);
 
    recreateIfNeeded = recreate;
 
    // put the results in the state
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    if (recreate &&
        (!simulator ||
         (simulator && (simType != simulator->GetType() ||
                        method != simulator->GetSimulationType() ||
                        simulator->GetNumberOfQubits() != numQubits))))
      CreateSimulator(simType, method);
  }

  void ExecuteOnHost(const std::shared_ptr<Circuits::Circuit<Time>> &circuit,
                     size_t hostId) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    recreateIfNeeded = false;
 
    const auto res = RepeatedExecuteOnHost(circuit, hostId, 1);
 
    recreateIfNeeded = recreate;
 
    // put the results in the state
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    if (recreate &&
        (!simulator ||
         (simulator && (simType != simulator->GetType() ||
                        method != simulator->GetSimulationType() ||
                        simulator->GetNumberOfQubits() != numQubits))))
      CreateSimulator(simType, method);
  }

  std::vector<double> ExecuteExpectations(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit,
      const std::vector<std::string> &paulis) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    recreateIfNeeded = false;
 
    pauliStrings = &paulis;
    const auto res = RepeatedExecute(circuit, 1);
    pauliStrings = nullptr;
 
    recreateIfNeeded = recreate;
 
    // put the results in the state
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    std::vector<double> expectations(paulis.size());
    if (simulator) {
      // translate the pauli strings to the mapped order of qubits
      const size_t numOps = simulator->GetNumberOfQubits();
 
      auto optimiser = controller->GetOptimiser();
      if (optimiser) {
        // convert the classical state results back to the expected order
        const auto &qubitsMap = optimiser->GetQubitsMap();
 
        for (size_t i = 0; i < paulis.size(); ++i) {
          std::string translated(numOps, 'I');
 
          for (size_t j = 0; j < paulis[i].size(); ++j) {
            const auto pos = qubitsMap.find(j);
            if (pos != qubitsMap.end())
              translated[pos->second] = paulis[i][j];
            else
              translated[j] = paulis[i][j];
          }
 
          expectations[i] = simulator->ExpectationValue(translated);
        }
      } else {
        for (size_t i = 0; i < paulis.size(); ++i)
          expectations[i] = simulator->ExpectationValue(paulis[i]);
      }
    }
 
    if (recreate && (!simulator || simType != simulator->GetType() ||
                     method != simulator->GetSimulationType() ||
                     simulator->GetNumberOfQubits() != numQubits))
      CreateSimulator(simType, method);
 
    return expectations;
  }

  std::vector<double> ExecuteOnHostExpectations(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t hostId,
      const std::vector<std::string> &paulis) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    // RAII: restore recreateIfNeeded and clear pauliStrings on any exit path.
    struct ScopedRestore {
      bool &flag, saved;
      const std::vector<std::string> **ps;
      ScopedRestore(bool &f, const std::vector<std::string> **p)
          : flag(f), saved(f), ps(p) {
        flag = false;
      }
      ~ScopedRestore() {
        flag = saved;
        *ps = nullptr;
      }
    } restoreGuard(recreateIfNeeded, &pauliStrings);
 
    pauliStrings = &paulis;
    const auto res = RepeatedExecuteOnHost(circuit, hostId, 1);
 
    // put the results in the state
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    // for (const auto& m : qubitsMapOnHost)
    //  std::cout << "Mapping qubit " << m.first << " to " << m.second <<
    // std::endl;
 
    const size_t offsetBase = qubitsMapOnHost.size();
 
    std::vector<double> expectations(paulis.size(), 1.);
    if (simulator) {
      // translate the pauli strings to the mapped order of qubits
      const size_t numOps = simulator->GetNumberOfQubits();
 
      // convert the pauli strings to the actual qubits order
      for (size_t i = 0; i < paulis.size(); ++i) {
        std::string translated(std::max(numOps, paulis[i].size()), 'I');
 
        size_t offset = offsetBase;
 
        for (size_t j = 0; j < paulis[i].size(); ++j) {
          auto pos = qubitsMapOnHost.find(j);
          if (pos != qubitsMapOnHost.end())
            translated[pos->second] = paulis[i][j];
          else {
            translated[offset] = paulis[i][j];
            ++offset;
          }
        }
 
        // std::cout << "Translated pauli string: " << translated << std::endl;
 
        expectations[i] = simulator->ExpectationValue(translated);
      }
    } else {
      throw std::runtime_error(
          "ExecuteOnHostExpectations: no simulator available after execution.");
    }
 
    if (recreate && (!simulator || simType != simulator->GetType() ||
                     method != simulator->GetSimulationType() ||
                     simulator->GetNumberOfQubits() != numQubits))
      CreateSimulator(simType, method);
 
    return expectations;
  }

  std::vector<std::complex<double>> ExecuteOnHostAmplitudes(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit,
      size_t hostId) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    // RAII: restore recreateIfNeeded on any exit path (including exceptions).
    struct ScopedRestoreFlag {
      bool &flag, saved;
      ScopedRestoreFlag(bool &f) : flag(f), saved(f) { flag = false; }
      ~ScopedRestoreFlag() { flag = saved; }
    } restoreGuard(recreateIfNeeded);
 
    const auto res = RepeatedExecuteOnHost(circuit, hostId, 1);
 
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    if (!simulator)
      throw std::runtime_error(
          "ExecuteOnHostAmplitudes: no simulator available after execution.");
 
    std::vector<std::complex<double>> amplitudes;
    const size_t n = simulator->GetNumberOfQubits();
    const size_t dim = 1ULL << n;
    amplitudes.resize(dim);
    for (size_t state = 0; state < dim; ++state)
      amplitudes[state] = simulator->Amplitude(state);
 
    // Remap amplitudes back to the original qubit ordering if qubits were
    // remapped during execution on the host.
    if (!qubitsMapOnHost.empty()) {
      const size_t offsetBase = qubitsMapOnHost.size();
 
      // Build reverse mapping: simulator qubit position -> original qubit
      // position.
      std::vector<size_t> simToOrig(n);
      size_t offset = offsetBase;
 
      for (size_t qbit = 0; qbit < n; ++qbit) {
        auto pos = qubitsMapOnHost.find(qbit);
        if (pos != qubitsMapOnHost.end())
          simToOrig[pos->second] = pos->first;
        else
          simToOrig[qbit] = offset++;
      }
 
      std::vector<std::complex<double>> remapped(dim);
 
      for (size_t sim_state = 0; sim_state < dim; ++sim_state) {
        size_t orig_state = 0;
        for (size_t qbit = 0; qbit < n; ++qbit) {
          if (sim_state & (1ULL << qbit))
            orig_state |= (1ULL << simToOrig[qbit]);
        }
        if (orig_state < dim) remapped[orig_state] = amplitudes[sim_state];
      }
      amplitudes.swap(remapped);
    }
 
    if (recreate && (!simulator || simType != simulator->GetType() ||
                     method != simulator->GetSimulationType() ||
                     simulator->GetNumberOfQubits() != numQubits))
      CreateSimulator(simType, method);
 
    return amplitudes;
  }

  std::complex<double> ExecuteOnHostProjectOnZero(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit,
      size_t hostId) override {
    const auto recreate = recreateIfNeeded;
 
    auto simType = Simulators::SimulatorType::kQCSim;
    auto method = Simulators::SimulationType::kMatrixProductState;
    size_t numQubits = 2;
    if (simulator) {
      simType = simulator->GetType();
      method = simulator->GetSimulationType();
      numQubits = simulator->GetNumberOfQubits();
    }
 
    // RAII: restore recreateIfNeeded on any exit path (including exceptions).
    struct ScopedRestoreFlag {
      bool &flag, saved;
      ScopedRestoreFlag(bool &f) : flag(f), saved(f) { flag = false; }
      ~ScopedRestoreFlag() { flag = saved; }
    } restoreGuard(recreateIfNeeded);
 
    const auto res = RepeatedExecuteOnHost(circuit, hostId, 1);
 
    if (!res.empty()) {
      const auto &first = *res.begin();
      GetState().SetResultsInOrder(first.first);
    }
 
    if (!simulator)
      throw std::runtime_error(
          "ExecuteOnHostProjectOnZero: no simulator available after "
          "execution.");
 
    const std::complex<double> result = simulator->ProjectOnZero();
 
    if (recreate && (!simulator || simType != simulator->GetType() ||
                     method != simulator->GetSimulationType() ||
                     simulator->GetNumberOfQubits() != numQubits))
      CreateSimulator(simType, method);
 
    return result;
  }

  ExecuteResults RepeatedExecute(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit,
      size_t shots = 1000) override {
    if (!controller || !circuit) return {};
 
    distCirc = controller->DistributeCircuit(BaseClass::getptr(), circuit);
    if (!distCirc) return {};
 
#ifdef _DEBUG
    for (auto q : distCirc->AffectedQubits()) {
      if (q >= GetNumQubits()) {
        std::cout
            << "This is a distributed circuit, using entanglement or cutting"
            << std::endl;
        break;
      }
    }
#endif
 
    if (!simulator) return {};
 
    auto simType = simulator->GetType();
    if (distCirc->HasOpsAfterMeasurements() &&
        (
#ifndef NO_QISKIT_AER
            simType == Simulators::SimulatorType::kCompositeQiskitAer ||
#endif
            simType == Simulators::SimulatorType::kCompositeQCSim))
      distCirc->MoveMeasurementsAndResets();
 
    auto method = simulator->GetSimulationType();
 
    const auto saveSimType = simType;
    const auto saveMethod = method;
 
    if (GetOptimizeSimulator() && distCirc->IsClifford() &&
        method != Simulators::SimulationType::kStabilizer
    // this is for the gpu simulator, as it doesn't support stabilizer
#ifdef __linux__
        && !Simulators::IsGpuSimulator(simType)
#endif
    ) {
      method = Simulators::SimulationType::kStabilizer;
 
      if (simType == Simulators::SimulatorType::kCompositeQCSim)
        simType = Simulators::SimulatorType::kQCSim;
#ifndef NO_QISKIT_AER
      else if (simType == Simulators::SimulatorType::kCompositeQiskitAer)
        simType = Simulators::SimulatorType::kQiskitAer;
#endif
    }
 
    ExecuteResults res;
    const size_t nrQubits = GetNumQubits() + GetNumNetworkEntangledQubits();
    const size_t nrCbitsResults = GetNumClassicalBits();
 
    CaptureSimulatorConfiguration();
 
    // do that only if the optimization for simulator is on and the estimator is
    // available, ortherwise an 'optimal' simulator won't be created
    if (optimizeSimulator && simulatorsEstimator &&
        simulatorsEstimator->IsInitialized()) {
      simulator->Clear();
      GetState().Clear();
    }
 
    curMaxBondDim = 0;
 
    std::vector<bool> executed;
    auto optSim =
        ChooseBestSimulator(distCirc, shots, nrQubits, nrQubits, nrCbitsResults,
                            simType, method, executed);
 
    lastSimulatorType = simType;
    lastMethod = method;
    lastGpuDevice = -1;
 
    size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
    if (Simulators::IsGpuSimulator(simType))
      nrThreads = 1;
    else
#endif
        if (((method == Simulators::SimulationType::kStatevector || method == Simulators::SimulationType::kPathIntegral) &&
             !distCirc->HasOpsAfterMeasurements()) ||
            simType == Simulators::SimulatorType::kQuestSim)
      nrThreads = 1;
 
    nrThreads = std::min(nrThreads, std::max<size_t>(shots, 1ULL));
 
    std::mutex resultsMutex;
 
    auto dcirc = distCirc;
 
    if (nrThreads > 1) {
      // since it's going to execute on multiple threads, free the memory from
      // the network's simulator and state, it's going to use other ones,
      // created in the threads if optimization already exists, it will be
      // cloned in the threads, otherwise a new one will be created in the
      // threads
      if (!optimizeSimulator || !simulatorsEstimator ||
          !simulatorsEstimator
               ->IsInitialized())  // otherwise it was already cleared
      {
        simulator->Clear();
        GetState().Clear();
      }
 
      const size_t cntPerThread = std::max<size_t>(shots / nrThreads, 1ULL);
 
      threadsPool.Resize(nrThreads);
      threadsPool.SetFinishLimit(shots);
 
      uint64_t jobStream = 0;
      while (shots > 0) {
        const size_t curCnt = std::min(cntPerThread, shots);
 
        shots -= curCnt;
 
        auto job = std::make_shared<ExecuteJob<Time>>(
            dcirc, res, curCnt, nrQubits, nrQubits, nrCbitsResults, simType,
            method, resultsMutex);
        job->optimiseMultipleShotsExecution = GetOptimizeSimulator();
 
        job->network = BaseClass::getptr();
        job->curMaxBondDim = &curMaxBondDim;
 
        job->config = ExecutionConfiguration(simType, nrQubits);
        if (configuration.IsSet("seed")) {
          const uint64_t childSeed = Simulators::IState::DeriveSeed(
              std::stoull(configuration.GetConfiguration("seed")), jobStream++);
          job->config.SetConfiguration("seed", std::to_string(childSeed));
        }
 
        if (optSim) {
          job->optSim = optSim->Clone();
          job->executedGates = executed;
          if (job->config.IsSet("seed"))
            job->optSim->SetSeed(
                std::stoull(job->config.GetConfiguration("seed")));
        }
 
        threadsPool.AddRunJob(std::move(job));
      }
 
      threadsPool.WaitForFinish();
      threadsPool.Stop();
    } else {
      const size_t curCnt = shots;
 
      auto job = std::make_shared<ExecuteJob<Time>>(
          dcirc, res, curCnt, nrQubits, nrQubits, nrCbitsResults, simType,
          method, resultsMutex);
      job->optimiseMultipleShotsExecution = GetOptimizeSimulator();
 
      job->network = BaseClass::getptr();
      job->curMaxBondDim = &curMaxBondDim;
 
      job->config = ExecutionConfiguration(simType, nrQubits);
 
      if (optSim) {
        optSim->SetMultithreading(true);
        job->optSim = optSim;
        job->executedGates = executed;
      } else {
        if (simulator && method == saveMethod && simType == saveSimType) {
          // use the already created simulator
          optSim = simulator;
          job->optSim = optSim;
          OptimizeMPSInitialQubitsMap(optSim, dcirc,
                                      optSim->GetNumberOfQubits());
          job->executedGates.resize(dcirc->size(),
                                    false);  // no gates executed yet
          simulator = nullptr;
        }
      }
 
      job->DoWorkNoLock();
      lastGpuDevice = job->optSim ? job->optSim->GetGpuDevice() : -1;
      if (!recreateIfNeeded) simulator = job->optSim;
    }
 
    if (recreateIfNeeded) CreateSimulator(saveSimType, saveMethod);
 
    ConvertBackResults(res);
 
    return res;
  }

  ExecuteResults RepeatedExecuteOnHost(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t hostId,
      size_t shots = 1000) override {
    if (!circuit || hostId >= GetNumHosts()) return {};
 
    size_t nrQubits = 0;
    size_t nrCbits = 0;
 
    const bool distributed = simulator &&
        Simulators::IsDistributedGpuSimulator(simulator->GetType());
    auto mappingCircuit = circuit;
    if (!distributed && GetController()->GetOptimizeCircuit()) {
      mappingCircuit = std::static_pointer_cast<Circuits::Circuit<Time>>(circuit->Clone());
      mappingCircuit->Optimize();
    }
    const auto reverseQubitsMap = MapCircuitOnHost(
        mappingCircuit, hostId, nrQubits, nrCbits, true);
    // Resolve indexing before optimization can remove a wire that disambiguates
    // local from global numbering (for example a cancelling pair on qubit 0).
    if (distributed && distCirc && GetController()->GetOptimizeCircuit())
      distCirc->Optimize();
    if (nrCbits == 0) nrCbits = nrQubits;
 
    if (!simulator || !distCirc) return {};
 
    auto simType = simulator->GetType();
 
    CaptureSimulatorConfiguration();
 
    if (distCirc->HasOpsAfterMeasurements() &&
        (
#ifndef NO_QISKIT_AER
            simType == Simulators::SimulatorType::kCompositeQiskitAer ||
#endif
            simType == Simulators::SimulatorType::kCompositeQCSim))
      distCirc->MoveMeasurementsAndResets();
 
    auto method = simulator->GetSimulationType();
    const auto saveSimType = simType;
    const auto saveMethod = method;
 
    if (GetOptimizeSimulator() && distCirc->IsClifford() &&
        method != Simulators::SimulationType::kStabilizer
    // this is for the gpu simulator, as it doesn't support stabilizer
#ifdef __linux__
        && !Simulators::IsGpuSimulator(simType)
#endif
    ) {
      method = Simulators::SimulationType::kStabilizer;
 
      if (simType == Simulators::SimulatorType::kCompositeQCSim)
        simType = Simulators::SimulatorType::kQCSim;
#ifndef NO_QISKIT_AER
      else if (simType == Simulators::SimulatorType::kCompositeQiskitAer)
        simType = Simulators::SimulatorType::kQiskitAer;
#endif
    }
 
    ExecuteResults res;
 
    // since it's going to execute on multiple threads, free the memory from the
    // network's simulator and state, it's going to use other ones, created in
    // the threads
    simulator->Clear();
    GetState().Clear();
 
    curMaxBondDim = 0;
 
    std::vector<bool> executed;
    auto optSim = ChooseBestSimulator(distCirc, shots, nrQubits, nrCbits,
                                      nrCbits, simType, method, executed);
 
    lastSimulatorType = simType;
    lastMethod = method;
    lastGpuDevice = -1;
 
    size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
    if (Simulators::IsGpuSimulator(simType))
      nrThreads = 1;
    else
#endif
        if (((method == Simulators::SimulationType::kStatevector ||
              method == Simulators::SimulationType::kPathIntegral) &&
             !distCirc->HasOpsAfterMeasurements()) ||
            simType == Simulators::SimulatorType::kQuestSim)
      nrThreads = 1;
 
    nrThreads = std::min(nrThreads, std::max<size_t>(shots, 1ULL));
 
    // WARNING: be sure to not put this above ChooseBestSimulator, as that one
    // can change the shots value!
 
    std::mutex resultsMutex;
 
    const auto dcirc = distCirc;
 
    if (nrThreads > 1) {
      // this rounds up, rounding down is better
      // const size_t cntPerThread = static_cast<size_t>((shots - 1) / nrThreads
      // + 1);
      const size_t cntPerThread = std::max<size_t>(shots / nrThreads, 1ULL);
 
      threadsPool.Resize(nrThreads);
      threadsPool.SetFinishLimit(shots);
 
      uint64_t jobStream = 0;
      while (shots > 0) {
        const size_t curCnt = std::min(cntPerThread, shots);
        shots -= curCnt;
 
        auto job = std::make_shared<ExecuteJob<Time>>(
            dcirc, res, curCnt, nrQubits, nrCbits, nrCbits, simType, method,
            resultsMutex);
        job->optimiseMultipleShotsExecution = GetOptimizeSimulator();
 
        job->network = BaseClass::getptr();
        job->curMaxBondDim = &curMaxBondDim;
 
        job->config = ExecutionConfiguration(simType, nrQubits);
        if (configuration.IsSet("seed")) {
          const uint64_t childSeed = Simulators::IState::DeriveSeed(
              std::stoull(configuration.GetConfiguration("seed")), jobStream++);
          job->config.SetConfiguration("seed", std::to_string(childSeed));
        }
 
        if (optSim) {
          job->optSim = optSim->Clone();
          job->executedGates = executed;
          if (job->config.IsSet("seed"))
            job->optSim->SetSeed(
                std::stoull(job->config.GetConfiguration("seed")));
        }
 
        threadsPool.AddRunJob(std::move(job));
      }
 
      threadsPool.WaitForFinish();
      threadsPool.Stop();
    } else {
      const size_t curCnt = shots;
 
      auto job = std::make_shared<ExecuteJob<Time>>(
          dcirc, res, curCnt, nrQubits, nrCbits, nrCbits, simType, method,
          resultsMutex);
      job->optimiseMultipleShotsExecution = GetOptimizeSimulator();
 
      job->network = BaseClass::getptr();
      job->curMaxBondDim = &curMaxBondDim;
 
      job->config = ExecutionConfiguration(simType, nrQubits);
 
      if (optSim) {
        optSim->SetMultithreading(true);
        job->optSim = optSim;
        job->executedGates = executed;
      }
 
      job->DoWorkNoLock();
      lastGpuDevice = job->optSim ? job->optSim->GetGpuDevice() : -1;
      if (!recreateIfNeeded) simulator = job->optSim;
    }
 
    if (recreateIfNeeded) CreateSimulator(saveSimType, saveMethod);
 
    if (!reverseQubitsMap.empty()) ConvertBackResults(res, reverseQubitsMap);
 
    return res;
  }

  size_t GetNumberOfGatesDistributedOrCut(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit) const override {
    if (!circuit) return 0;
 
    size_t distgates = 0;
 
    for (const auto &op : circuit->GetOperations())
      if (!IsLocalOperation(op)) ++distgates;
 
    return distgates;
  }

  std::vector<ExecuteResults> ExecuteScheduled(
      const std::vector<Schedulers::ExecuteCircuit<Time>> &circuits) override {
    // create a default one if not set
    if (!GetScheduler()) {
      CreateScheduler();
 
      if (!GetScheduler()) return {};
    }
 
    return GetScheduler()->ExecuteScheduled(circuits);
  }

  void CreateSimulator(
      Simulators::SimulatorType simType = Simulators::SimulatorType::kQCSim,
      Simulators::SimulationType simExecType =
          Simulators::SimulationType::kMatrixProductState,
      size_t nrQubits = 0) override {
    classicalState.Clear();
    classicalState.AllocateBits(GetNumClassicalBits() +
                                GetNumNetworkEntangledQubits());
 
    simulator =
        Simulators::SimulatorsFactory::CreateSimulator(simType, simExecType);
 
    if (simulator) {
      const size_t allocationQubits = nrQubits == 0
          ? GetNumQubits() + GetNumNetworkEntangledQubits() : nrQubits;
      ExecutionConfiguration(simType, allocationQubits)
          .ApplyConfigurationToSimulator(simulator);
 
      simulator->AllocateQubits(allocationQubits);
      simulator->Initialize();
 
      // Pin the resolved default as well as explicit selections. Cloning or
      // recreating this network must not follow later process-default changes.
      if (simType == Simulators::SimulatorType::kGpuSim)
        configuration.SetConfiguration(
            "gpu_device", std::to_string(simulator->GetGpuDevice()));
 
      if (Simulators::IsDistributedGpuSimulator(simType) && nrQubits == 0) {
        resolvedDistributedDevices[simType] =
            simulator->GetConfiguration("distributed_shard_devices");
      }
      simulator->setGrowthFactorGate(growthFactorGate);
      simulator->setGrowthFactorSwap(growthFactorSwap);
      simulator->SetLookaheadDepth(lookaheadDepth);
      simulator->SetLookaheadDepthWithHeuristic(lookaheadDepthWithHeuristic);
    }
  }

  void Configure(const char *key, const char *value) override {
    if (!key || !value) return;
 
    if (std::string("distributed_host_qubit_indexing") == key) {
      const std::string mode(value);
      if (mode != "auto" && mode != "local" && mode != "global")
        throw std::invalid_argument(
            "distributed_host_qubit_indexing must be auto, local or global");
      distributedHostQubitIndexing = mode;
      return;
    }
 
    if (std::string("distributed_devices") == key ||
        std::string("gpu_device") == key)
      resolvedDistributedDevices.clear();
 
    if (std::string("max_simulators") == key)
      maxSimulators = std::stoull(value);
 
    if (std::string("gpu_device") == key) {
      Simulators::Configuration::ParseGpuDevice(value);
      if (simulator) simulator->Configure(key, value);
      configuration.SetConfiguration(key, value);
      return;
    }
    if (!Simulators::Configuration::GpuSvdSettingGroup(key).empty()) {
      if (simulator) simulator->Configure(key, value);
      configuration.SetConfiguration(key, value);
      return;
    }
    configuration.SetConfiguration(key, value);
 
    if (simulator) simulator->Configure(key, value);
  }

  std::shared_ptr<Simulators::ISimulator> GetSimulator() const override {
    return simulator;
  }

  Circuits::OperationState &GetState() override { return classicalState; }

  void CreateScheduler(
      SchedulerType schType =
          SchedulerType::kNoEntanglementQubitsParallel) override {
    if (!controller) return;
 
    controller->CreateScheduler(BaseClass::getptr(), schType);
  }

  std::shared_ptr<Schedulers::IScheduler<Time>> GetScheduler() const override {
    if (!controller) return nullptr;
 
    return controller->GetScheduler();
  }

  const std::shared_ptr<IHost<Time>> GetHost(size_t hostId) const override {
    if (hostId >= hosts.size()) return nullptr;
 
    return hosts[hostId];
  }

  const std::shared_ptr<IController<Time>> GetController() const override {
    return controller;
  }

  size_t GetNumHosts() const override { return hosts.size(); }

  size_t GetNumQubits() const override {
    size_t res = 0;
 
    for (const auto &host : hosts) res += host->GetNumQubits();
 
    return res;
  }

  size_t GetNumQubitsForHost(size_t hostId) const override {
    if (hostId >= hosts.size()) return 0;
 
    return hosts[hostId]->GetNumQubits();
  }

  size_t GetNumNetworkEntangledQubits() const override {
    size_t res = 0;
 
    for (const auto &host : hosts) res += host->GetNumNetworkEntangledQubits();
 
    return res;
  }

  size_t GetNumNetworkEntangledQubitsForHost(size_t hostId) const override {
    if (hostId >= hosts.size()) return 0;
 
    return hosts[hostId]->GetNumNetworkEntangledQubits();
  }

  size_t GetNumClassicalBits() const override {
    size_t res = 0;
 
    for (const auto &host : hosts) res += host->GetNumClassicalBits();
 
    return res;
  }

  size_t GetNumClassicalBitsForHost(size_t hostId) const override {
    if (hostId >= hosts.size()) return 0;
 
    return hosts[hostId]->GetNumClassicalBits();
  }

  std::vector<std::shared_ptr<IHost<Time>>> &GetHosts() { return hosts; }

  void SetController(const std::shared_ptr<IController<Time>> &cntrl) {
    controller = cntrl;
  }

  bool SendPacket(size_t fromHostId, size_t toHostId,
                  const std::vector<uint8_t> &packet) override {
    return false;
  }

  NetworkType GetType() const override {
    return NetworkType::kSimpleDisconnectedNetwork;
  }

  bool IsLocalOperation(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    const auto qubits = op->AffectedQubits();
 
    if (qubits.empty()) return true;
 
    size_t firstQubit = qubits[0];
 
    for (size_t q = 1; q < qubits.size(); ++q)
      if (!AreQubitsOnSameHost(firstQubit, qubits[q])) return false;
 
    return true;
  }

  bool IsDistributedOperation(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    const auto qubits = op->AffectedQubits();
 
    if (qubits.empty()) return false;
 
    // grab the first qubit that is on a host (skip over network entangled
    // qubits)
    size_t firstQubit = qubits[0];
    size_t q = 1;
    for (; IsNetworkEntangledQubit(firstQubit) && q < qubits.size(); ++q)
      firstQubit = qubits[q];
 
    // check to see one of the other qubits is on a different host, but ignore
    // the network entangled qubits
    for (; q < qubits.size(); ++q)
      if (!IsNetworkEntangledQubit(qubits[q]) &&
          !AreQubitsOnSameHost(firstQubit, qubits[q]))
        return true;
 
    return false;
  }

  bool OperatesWithNetworkEntangledQubit(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    const auto qubits = op->AffectedQubits();
 
    if (qubits.empty()) return false;
 
    for (size_t q = 0; q < qubits.size(); ++q)
      if (IsNetworkEntangledQubit(qubits[q])) return true;
 
    return false;
  }

  bool IsEntanglingGate(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    if (op->GetType() != Circuits::OperationType::kGate) return false;
    const auto qubits = op->AffectedQubits();
    if (qubits.size() != 2) return false;
 
    return IsNetworkEntangledQubit(qubits[0]) &&
           IsNetworkEntangledQubit(qubits[1]);
  }

  bool ExpectsClassicalBitFromOtherHost(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    if (!op->IsConditional()) return false;
 
    const auto qubits = op->AffectedQubits();
 
    const std::shared_ptr<Circuits::IConditionalOperation<Time>> condOp =
        std::static_pointer_cast<Circuits::IConditionalOperation<Time>>(op);
    const auto &classicalBits = condOp->GetCondition()->GetBitsIndices();
 
    if (qubits.empty() && classicalBits.empty())
      throw std::runtime_error(
          "No classical bits specified!");  // this would be odd!
 
    // consider it on the host where it has the first qubit (or bit, if there
    // are no qubits)
 
    const size_t hostId = GetHostIdForAnyQubit(qubits[0]);
 
    // now check the classical bits
    for (const auto bit : classicalBits)
      if (hostId != GetHostIdForClassicalBit(bit)) return true;
 
    return false;
  }

  size_t GetHostIdForClassicalControl(
      const std::shared_ptr<Circuits::IOperation<Time>> &op) const override {
    if (!op->IsConditional())
      throw std::runtime_error("Operation is not conditional!");
 
    std::shared_ptr<Circuits::IConditionalOperation<Time>> condOp =
        std::static_pointer_cast<Circuits::IConditionalOperation<Time>>(op);
    const auto classicalBits = condOp->AffectedBits();
 
    if (classicalBits.empty())
      throw std::runtime_error("No classical bits specified!");
 
    return GetHostIdForClassicalBit(classicalBits[0]);
  }

  bool AreQubitsOnSameHost(size_t qubitId1, size_t qubitId2) const override {
    for (const auto &host : hosts) {
      const bool present1 = host->IsQubitOnHost(qubitId1);
      const bool present2 = host->IsQubitOnHost(qubitId2);
 
      if (present1 && present2)
        return true;
      else if (present1 || present2)
        return false;
    }
 
    return false;
  }

  bool AreClassicalBitsOnSameHost(size_t bitId1, size_t bitId2) const override {
    for (const auto &host : hosts) {
      const bool present1 = host->IsClassicalBitOnHost(bitId1);
      const bool present2 = host->IsClassicalBitOnHost(bitId2);
 
      if (present1 && present2)
        return true;
      else if (present1 || present2)
        return false;
    }
 
    return false;
  }

  bool AreQubitAndClassicalBitOnSameHost(size_t qubitId,
                                         size_t bitId) const override {
    for (const auto &host : hosts) {
      const bool present1 = host->IsQubitOnHost(qubitId);
      const bool present2 = host->IsClassicalBitOnHost(bitId);
 
      if (present1 && present2)
        return true;
      else if (present1 || present2)
        return false;
    }
 
    return false;
  }

  size_t GetHostIdForQubit(size_t qubitId) const override {
    for (const auto &host : hosts)
      if (host->IsQubitOnHost(qubitId)) return host->GetId();
 
    return std::numeric_limits<size_t>::max();
  }

  size_t GetHostIdForEntangledQubit(size_t qubitId) const override {
    for (const auto &host : hosts)
      if (host->IsEntangledQubitOnHost(qubitId)) return host->GetId();
 
    return std::numeric_limits<size_t>::max();
  }

  size_t GetHostIdForAnyQubit(size_t qubitId) const override {
    if (IsNetworkEntangledQubit(qubitId))
      return GetHostIdForEntangledQubit(qubitId);
 
    return GetHostIdForQubit(qubitId);
  }

  size_t GetHostIdForClassicalBit(size_t classicalBitId) const override {
    for (const auto &host : hosts)
      if (host->IsClassicalBitOnHost(classicalBitId)) return host->GetId();
 
    return std::numeric_limits<size_t>::max();
  }

  std::vector<size_t> GetQubitsIds(size_t hostId) const override {
    if (hostId >= hosts.size()) return std::vector<size_t>();
 
    return hosts[hostId]->GetQubitsIds();
  }

  std::vector<size_t> GetNetworkEntangledQubitsIds(
      size_t hostId) const override {
    if (hostId >= hosts.size()) return std::vector<size_t>();
 
    return hosts[hostId]->GetNetworkEntangledQubitsIds();
  }

  std::vector<size_t> GetClassicalBitsIds(size_t hostId) const override {
    if (hostId >= hosts.size()) return std::vector<size_t>();
 
    return hosts[hostId]->GetClassicalBitsIds();
  }

  std::vector<size_t> GetEntangledQubitMeasurementBitIds(
      size_t hostId) const override {
    if (hostId >= hosts.size()) return std::vector<size_t>();
 
    return hosts[hostId]->GetEntangledQubitMeasurementBitIds();
  }

  bool IsNetworkEntangledQubit(size_t qubitId) const override {
    return qubitId >= GetNumQubits();
  }

  bool IsEntanglementQubitBusy(size_t qubitId) const override { return false; }

  bool AreEntanglementQubitsBusy(size_t qubitId1,
                                 size_t qubitId2) const override {
    return false;
  }

  void MarkEntangledQubitsBusy(size_t qubitId1, size_t qubitId2) override {
    throw std::runtime_error(
        "Entanglement between hosts is not supported in the simple network");
  }

  void MarkEntangledQubitFree(size_t qubitId) override {
    throw std::runtime_error(
        "Entanglement between hosts is not supported in the simple network");
  }

  void ClearEntanglements() override {
    throw std::runtime_error(
        "Entanglement between hosts is not supported in the simple network");
  }

  std::shared_ptr<Circuits::Circuit<Time>> GetDistributedCircuit()
      const override {
    return distCirc;
  }

  int GetLastGpuDevice() const override { return lastGpuDevice; }
 
  Simulators::SimulatorType GetLastSimulatorType() const override {
    return lastSimulatorType;
  }

  Simulators::SimulationType GetLastSimulationType() const override {
    return lastMethod;
  }

  size_t GetMaxSimulators() const override { return maxSimulators; }

  void SetMaxSimulators(size_t val) override {
    if (val < 1)
      val = 1;
    else if (val > (size_t)QC::QubitRegisterCalculator<>::GetNumberOfThreads())
      val = (size_t)QC::QubitRegisterCalculator<>::GetNumberOfThreads();
 
    maxSimulators = val;
  }

  void SetOptimizeSimulator(bool optimize = true) override {
    optimizeSimulator = optimize;
  }

  bool GetOptimizeSimulator() const override { return optimizeSimulator; }

  const typename BaseClass::SimulatorsSet &GetSimulatorsSet() const override {
    return simulatorsForOptimizations;
  }

  void AddOptimizationSimulator(Simulators::SimulatorType type,
                                Simulators::SimulationType kind) override {
    simulatorsForOptimizations.insert({type, kind});
  }

  void RemoveOptimizationSimulator(Simulators::SimulatorType type,
                                   Simulators::SimulationType kind) override {
    simulatorsForOptimizations.erase({type, kind});
  }

  void RemoveAllOptimizationSimulatorsAndAdd(
      Simulators::SimulatorType type,
      Simulators::SimulationType kind) override {
    simulatorsForOptimizations.clear();
    simulatorsForOptimizations.insert({type, kind});
  }

  bool OptimizationSimulatorExists(
      Simulators::SimulatorType type,
      Simulators::SimulationType kind) const override {
    if (simulatorsForOptimizations.empty()) return true;
 
    return simulatorsForOptimizations.find({type, kind}) !=
           simulatorsForOptimizations.end();
  }

  std::shared_ptr<INetwork<Time>> Clone() const override {
    const size_t numHosts = GetNumHosts();
 
    std::vector<Types::qubit_t> qubits(numHosts);
    std::vector<size_t> cbits(numHosts);
 
    for (size_t h = 0; h < numHosts; ++h) {
      qubits[h] = GetNumQubitsForHost(h);
      cbits[h] = GetNumClassicalBitsForHost(h);
    }
 
    const auto cloned =
        std::make_shared<SimpleDisconnectedNetwork<Time, Controller>>(qubits,
                                                                      cbits);
    
    cloned->configuration = configuration;
    cloned->distributedHostQubitIndexing = distributedHostQubitIndexing;
    cloned->resolvedDistributedDevices = resolvedDistributedDevices;
 
    cloned->maxSimulators = maxSimulators;
 
    cloned->optimizeSimulator = optimizeSimulator;
    cloned->simulatorsForOptimizations = simulatorsForOptimizations;
 
    cloned->SetMPSOptimizeSwaps(GetMPSOptimizeSwaps());
 
    cloned->SetMPSOptimizationBondDimensionThreshold(GetMPSOptimizationBondDimensionThreshold());
    cloned->SetMPSOptimizationQubitsNumberThreshold(GetMPSOptimizationQubitsNumberThreshold());
 
    cloned->SetLookaheadDepth(GetLookaheadDepth());
    cloned->SetLookaheadDepthWithHeuristic(GetLookaheadDepthWithHeuristic());
 
    cloned->setGrowthFactorGate(getGrowthFactorGate());
    cloned->setGrowthFactorSwap(getGrowthFactorSwap());
 
    if (GetSimulator())
      cloned->CreateSimulator(GetSimulator()->GetType(),
                              GetSimulator()->GetSimulationType());
 
    return cloned;
  }
 
  std::shared_ptr<Simulators::ISimulator> ChooseBestSimulator(
      std::shared_ptr<Circuits::Circuit<Time>> &dcirc, size_t &counts,
      size_t nrQubits, size_t nrCbits, size_t nrResultCbits,
      Simulators::SimulatorType &simType, Simulators::SimulationType &method,
      std::vector<bool> &executed, bool multithreading = false,
      bool dontRunCircuitStart = false) override {
    // Distribution is an explicit execution choice. Timing-based backend
    // selection must not replace it or diverge between MPI ranks.
    if (Simulators::IsDistributedGpuSimulator(simType)) return nullptr;
    if (!optimizeSimulator) return nullptr;
 
    if ((!simulatorsEstimator || !simulatorsEstimator->IsInitialized()) &&
        simulatorsForOptimizations.size() != 1)
      return nullptr;
 
    // when multithreading is set to true it means it needs a multithreaded
    // simulator
 
    std::vector<
        std::pair<Simulators::SimulatorType, Simulators::SimulationType>>
        simulatorTypes;
 
    const bool checkTensorNetwork =
        method == Simulators::SimulationType::kTensorNetwork;
 
    // the others are to be picked between statevector, composite, tensor
    // networks and mps, for now at least for tensor networks in the future it's
    // worth checking different contractors!!!!
    //
    // clifford was decided at higher level
    if (method == Simulators::SimulationType::kStabilizer) {
      // compare qcsim with qiskit aer if qiskit aer is available, let the best
      // one win
      if (OptimizationSimulatorExists(Simulators::SimulatorType::kQCSim,
                                      Simulators::SimulationType::kStabilizer))
        simulatorTypes.emplace_back(Simulators::SimulatorType::kQCSim,
                                    Simulators::SimulationType::kStabilizer);
 
#ifndef NO_QISKIT_AER
      // if the number of shots is too small, probably it's not worth it, it's
      // going to be better to just execute them multithreading
      if (OptimizationSimulatorExists(Simulators::SimulatorType::kQiskitAer,
                                      Simulators::SimulationType::kStabilizer))
        simulatorTypes.emplace_back(Simulators::SimulatorType::kQiskitAer,
                                    Simulators::SimulationType::kStabilizer);
#endif
    }
 
    if (OptimizationSimulatorExists(Simulators::SimulatorType::kQCSim,
                                    Simulators::SimulationType::kStatevector))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQCSim,
                                  Simulators::SimulationType::kStatevector);
 
    if (OptimizationSimulatorExists(Simulators::SimulatorType::kCompositeQCSim,
                                    Simulators::SimulationType::kStatevector))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kCompositeQCSim,
                                  Simulators::SimulationType::kStatevector);
 
    if (checkTensorNetwork &&
        OptimizationSimulatorExists(Simulators::SimulatorType::kQCSim,
                                    Simulators::SimulationType::kTensorNetwork))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQCSim,
                                  Simulators::SimulationType::kTensorNetwork);
 
    const long long int maxBondDim = configuration.GetConfigurationAsInt(
        "matrix_product_state_max_bond_dimension");
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kQCSim,
            Simulators::SimulationType::kMatrixProductState) &&
        (nrQubits <= 4 || maxBondDim > 0))
      simulatorTypes.emplace_back(
          Simulators::SimulatorType::kQCSim,
          Simulators::SimulationType::kMatrixProductState);
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kQCSim,
            Simulators::SimulationType::kPauliPropagator))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQCSim,
                                  Simulators::SimulationType::kPauliPropagator);
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kQCSim,
            Simulators::SimulationType::kPathIntegral))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQCSim,
                                  Simulators::SimulationType::kPathIntegral);
 
#ifndef NO_QISKIT_AER
    // tensor networks are out of the picture for now for qiskit aer, since they
    // are available with cuda library, and work only on linux (obviously when
    // compiled properly and if there is the right hw an driver installed)
 
    if (OptimizationSimulatorExists(Simulators::SimulatorType::kQiskitAer,
                                    Simulators::SimulationType::kStatevector))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQiskitAer,
                                  Simulators::SimulationType::kStatevector);
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kCompositeQiskitAer,
            Simulators::SimulationType::kStatevector))
      simulatorTypes.emplace_back(
          Simulators::SimulatorType::kCompositeQiskitAer,
          Simulators::SimulationType::kStatevector);
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kQiskitAer,
            Simulators::SimulationType::kMatrixProductState) &&
        (nrQubits <= 4 || maxBondDim > 0))
      simulatorTypes.emplace_back(
          Simulators::SimulatorType::kQiskitAer,
          Simulators::SimulationType::kMatrixProductState);
#endif
 
#ifdef __linux__
    const int gpuDevice = configuration.IsSet("gpu_device")
        ? Simulators::Configuration::ParseGpuDevice(configuration.GetConfiguration("gpu_device"))
        : Simulators::SimulatorsFactory::ResolveGpuDevice();
    Simulators::SimulatorsFactory::ScopedGpuDevice gpuDeviceScope(gpuDevice);
    const bool hasGpuCandidate = std::any_of(
        simulatorsForOptimizations.begin(), simulatorsForOptimizations.end(),
        [](const auto& candidate) { return candidate.first == Simulators::SimulatorType::kGpuSim; });
    if (configuration.IsSet("gpu_device") && hasGpuCandidate &&
        !Simulators::SimulatorsFactory::IsGpuLibraryAvailable(gpuDevice))
      throw std::runtime_error("Unable to initialize requested GPU device " + std::to_string(gpuDevice));
    if (hasGpuCandidate && Simulators::SimulatorsFactory::IsGpuLibraryAvailable(gpuDevice)) {
      if (OptimizationSimulatorExists(Simulators::SimulatorType::kGpuSim,
                                      Simulators::SimulationType::kStatevector))
        simulatorTypes.emplace_back(Simulators::SimulatorType::kGpuSim,
                                    Simulators::SimulationType::kStatevector);
      if (OptimizationSimulatorExists(
              Simulators::SimulatorType::kGpuSim,
              Simulators::SimulationType::kMatrixProductState))
        simulatorTypes.emplace_back(
            Simulators::SimulatorType::kGpuSim,
            Simulators::SimulationType::kMatrixProductState);
      if (OptimizationSimulatorExists(
              Simulators::SimulatorType::kGpuSim,
              Simulators::SimulationType::kTensorNetwork))
        simulatorTypes.emplace_back(Simulators::SimulatorType::kGpuSim,
                                    Simulators::SimulationType::kTensorNetwork);
      if (OptimizationSimulatorExists(
              Simulators::SimulatorType::kGpuSim,
              Simulators::SimulationType::kPauliPropagator))
        simulatorTypes.emplace_back(
            Simulators::SimulatorType::kGpuSim,
            Simulators::SimulationType::kPauliPropagator);
    }
#endif
 
    if (OptimizationSimulatorExists(Simulators::SimulatorType::kQuestSim,
                                    Simulators::SimulationType::kStatevector))
      simulatorTypes.emplace_back(Simulators::SimulatorType::kQuestSim,
                                  Simulators::SimulationType::kStatevector);
 
 
    // Honor a singleton optimization set (e.g. density_matrix) even if it is
    // not in the hardcoded candidate list. Skip backends that cannot actually
    // be constructed, such as GPU MPS when the GPU library is missing: the
    // caller then keeps the network simulator instead of recording 0 shots.
    if (simulatorTypes.empty() && simulatorsForOptimizations.size() == 1) {
      const auto candidate = *simulatorsForOptimizations.begin();
      if (candidate.first == Simulators::SimulatorType::kGpuSim &&
          !Simulators::SimulatorsFactory::IsGpuLibraryAvailable()) return nullptr;
      if (Simulators::SimulatorsFactory::CreateSimulator(candidate.first,
                                                         candidate.second))
        simulatorTypes.push_back(candidate);
    }
 
    if (simulatorTypes.empty())
      return nullptr;
    else if (simulatorTypes.size() == 1) {
      const auto candidateType = simulatorTypes[0].first;
      const auto candidateMethod = simulatorTypes[0].second;
 
      std::shared_ptr<Simulators::ISimulator> sim =
          Simulators::SimulatorsFactory::CreateSimulator(candidateType,
                                                         candidateMethod);
      if (sim) {
        simType = candidateType;
        method = candidateMethod;
        configuration.ApplyConfigurationToSimulator(sim);
        
        if (method == Simulators::SimulationType::kMatrixProductState) {
          sim->AllocateQubits(nrQubits);
          sim->Initialize();
 
          sim->setGrowthFactorGate(growthFactorGate);
          sim->setGrowthFactorSwap(growthFactorSwap);
          sim->SetLookaheadDepth(lookaheadDepth);
          sim->SetLookaheadDepthWithHeuristic(lookaheadDepthWithHeuristic);
 
          OptimizeMPSInitialQubitsMap(sim, dcirc, nrQubits);
        } else {
          sim->AllocateQubits(nrQubits);
          sim->Initialize();
        }
 
        if (!dontRunCircuitStart) {
          sim->SetMultithreading(true);
          Estimators::SimulatorsEstimatorInterface<
              Time>::ExecuteUpToMeasurements(dcirc, nrQubits, nrCbits,
                                             nrResultCbits, sim, executed, &curMaxBondDim);
        }
        sim->SetMultithreading(multithreading || GetMaxSimulators() == 1);
 
        return sim;
      }
 
      // Singleton backend cannot be constructed (typical for GPU on CI).
      // Leave simType/method unchanged so RepeatedExecuteOnHost keeps the
      // network simulator instead of recording empty counts.
      return nullptr;
    }
 
    const double singularValueThreshold =
        configuration.GetConfigurationAsDouble(
            "matrix_product_state_truncation_threshold");
 
    const std::string truncationMode = configuration.GetConfiguration(
        "matrix_product_state_truncation_mode");
 
    const std::string mpsSample = configuration.GetConfiguration(
        "mps_sample_measure_algorithm");
 
    std::shared_ptr<Simulators::ISimulator> sim =
        simulatorsEstimator->ChooseBestSimulator(
            simulatorTypes, dcirc, counts, nrQubits, nrCbits, nrResultCbits,
            simType, method, executed, maxBondDim, singularValueThreshold,
            truncationMode, mpsSample, GetMaxSimulators(), pauliStrings, multithreading);
 
    if (sim) {
      sim->AllocateQubits(nrQubits);
      sim->Initialize();
 
      sim->setGrowthFactorGate(growthFactorGate);
      sim->setGrowthFactorSwap(growthFactorSwap);
      sim->SetLookaheadDepth(lookaheadDepth);
      sim->SetLookaheadDepthWithHeuristic(lookaheadDepthWithHeuristic);
 
      OptimizeMPSInitialQubitsMap(sim, dcirc, nrQubits);
 
      if (!dontRunCircuitStart) {
        sim->SetMultithreading(true);
        Estimators::SimulatorsEstimatorInterface<Time>::ExecuteUpToMeasurements(
            dcirc, nrQubits, nrCbits, nrResultCbits, sim, executed, &curMaxBondDim);
      }
      sim->SetMultithreading(multithreading || GetMaxSimulators() == 1);
    }
 
    return sim;
  }
 
  void SetInitialQubitsMapOptimization(bool optimize = true) override {
    optimizeInitialQubitsMap = optimize;
  }
 
  bool GetInitialQubitsMapOptimization() const override {
    return optimizeInitialQubitsMap;
  }
 
  void SetMPSOptimizeSwaps(bool optimize = true) override {
    mpsOptimizeSwaps = optimize;
 
    if (simulator) {
      simulator->SetLookaheadDepth(0);
      simulator->SetLookaheadDepthWithHeuristic(0);
    }
  }
 
  bool GetMPSOptimizeSwaps() const override { return mpsOptimizeSwaps; }
 
  void SetMPSOptimizationBondDimensionThreshold(size_t threshold) override {
    mpsOptimizationBondDimensionThreshold = threshold;
 
    if (simulator &&
        std::stoull(simulator->GetConfiguration(
            "matrix_product_state_max_bond_dimension")) < threshold) {
      simulator->SetLookaheadDepth(0);
      simulator->SetLookaheadDepthWithHeuristic(0);
    }
  }
 
  size_t GetMPSOptimizationBondDimensionThreshold() const override {
    return mpsOptimizationBondDimensionThreshold;
  }
 
  void SetMPSOptimizationQubitsNumberThreshold(size_t threshold) override {
    mpsOptimizationQubitsNumberThreshold = threshold;
 
    if (GetNumQubits() < threshold && simulator) {
      simulator->SetLookaheadDepth(0);
      simulator->SetLookaheadDepthWithHeuristic(0);
    }
  }
 
  size_t GetMPSOptimizationQubitsNumberThreshold() const override {
    return mpsOptimizationQubitsNumberThreshold;
  }
 
  void SetLookaheadDepth(int depth) override {
    if (depth < 0) depth = std::numeric_limits<int>::max();
 
    lookaheadDepth = depth;
 
    if (simulator && lookaheadDepth != std::numeric_limits<int>::max()) {
      simulator->SetLookaheadDepth(0);
      simulator->SetLookaheadDepthWithHeuristic(0);
    }
  }
 
  int GetLookaheadDepth() const override { return lookaheadDepth; }
 
  void SetLookaheadDepthWithHeuristic(int depth) override {
    if (depth < 0) depth = std::numeric_limits<int>::max();
 
    if (depth > lookaheadDepth) depth = lookaheadDepth;
 
    lookaheadDepthWithHeuristic = depth;
 
    if (simulator && lookaheadDepthWithHeuristic != std::numeric_limits<int>::max())
      simulator->SetLookaheadDepthWithHeuristic(depth);
  }
 
  int GetLookaheadDepthWithHeuristic() const override {
    return lookaheadDepthWithHeuristic;
  }
 
  double getGrowthFactorSwap() const override { return growthFactorSwap; }
  double getGrowthFactorGate() const override { return growthFactorGate; }
 
  void setGrowthFactorSwap(double factor) override {
    growthFactorSwap = factor;
 
    if (simulator) simulator->setGrowthFactorSwap(factor);
  }
 
  void setGrowthFactorGate(double factor) override {
    growthFactorGate = factor;
 
    if (simulator) simulator->setGrowthFactorGate(factor);
  }

  size_t GetCurrentMaxBondDimension() const override { return curMaxBondDim; }
 
 protected:
  // Resolved placement is separate from user settings. Importing a smaller
  // simulator must not turn its automatic shard subset into an explicit choice.
  void CaptureSimulatorConfiguration() {
    const auto requested = configuration.GetConfiguration("distributed_devices");
    configuration.ApplyConfigurationFromSimulator(simulator);
    if (simulator && Simulators::IsDistributedGpuSimulator(simulator->GetType()))
      configuration.SetConfiguration("distributed_devices", requested);
  }
 
  Configuration<Time> ExecutionConfiguration(Simulators::SimulatorType type,
                                             size_t qubits) const {
    auto result = configuration;
    const auto found = resolvedDistributedDevices.find(type);
    if (!Simulators::IsDistributedGpuSimulator(type) ||
        !configuration.GetConfiguration("distributed_devices").empty() ||
        found == resolvedDistributedDevices.end())
      return result;
    auto devices = found->second;
    // MPI shard count belongs to the communicator. Explicit global-qubit
    // settings also constrain the shard count and must be validated unchanged.
    if (type == Simulators::SimulatorType::kDistGpuSim &&
        !configuration.IsSet("distributed_global_qubits")) {
      size_t count = 1;
      for (char c : devices) if (c == ',') ++count;
      size_t bits = 0;
      for (size_t n = count; n > 1; n >>= 1) ++bits;
      while (count > 1 && bits >= qubits) {
        count /= 2;
        --bits;
      }
      size_t end = 0;
      for (size_t i = 0; i < count; ++i) {
        end = devices.find(',', end);
        if (end == std::string::npos) break;
        if (i + 1 < count) ++end;
      }
      devices = devices.substr(0, end);
    }
    result.SetConfiguration("distributed_devices", devices);
    return result;
  }
 
  void OptimizeMPSInitialQubitsMap(
      std::shared_ptr<Simulators::ISimulator> &sim,
      std::shared_ptr<Circuits::Circuit<Time>> &dcirc, size_t nrQubits) const {
    if (sim->GetSimulationType() ==
            Simulators::SimulationType::kMatrixProductState &&
        (optimizeInitialQubitsMap || mpsOptimizeSwaps) &&
        sim->SupportsMPSSwapOptimization()) {
      if (mpsOptimizationQubitsNumberThreshold <= nrQubits) {
        const auto maxBondDimValue =
            configuration.GetConfigurationAsInt("matrix_product_state_max_bond_dimension");
 
        if (maxBondDimValue <= 0 ||
            static_cast<int>(mpsOptimizationBondDimensionThreshold) <= maxBondDimValue) {
          // need to be sure the circuit is correctly converted
          dcirc->ConvertForCutting();  // convert the three qubit gates
          auto layers = dcirc->ToMultipleQubitsLayersNoClone();
 
          Simulators::MPSDummySimulator dummySim(nrQubits);
          dummySim.setGrowthFactorGate(growthFactorGate);
          dummySim.setGrowthFactorSwap(growthFactorSwap);
 
          if (maxBondDimValue > 0)
            dummySim.SetMaxBondDimension(maxBondDimValue);
 
          if (optimizeInitialQubitsMap) {
            const auto optimalMap = dummySim.ComputeOptimalQubitsMap(layers);
            sim->SetInitialQubitsMap(optimalMap);
          }
 
          auto optCirc = Circuits::Circuit<Time>::LayersToCircuit(layers);
          dcirc->SetOperations(optCirc->GetOperations());
 
          if (mpsOptimizeSwaps) {
            // TODO: come up with something better!
            int lookaheadDepthLocal = lookaheadDepth;
 
            if (lookaheadDepthLocal == std::numeric_limits<int>::max()) {
              double avgTwoQubitGatesPerLayer = 0.0;
              for (const auto &layer : layers) {
                int twoQubitGates = 0;
                for (const auto &op : layer->GetOperations()) {
                  if (op->AffectedQubits().size() >= 2) {
                    ++twoQubitGates;
                  }
                }
                avgTwoQubitGatesPerLayer += twoQubitGates;
              }
              avgTwoQubitGatesPerLayer /= layers.size();
 
              int lookaheadVal = static_cast<int>(4. * avgTwoQubitGatesPerLayer);
              if (lookaheadVal > 15) lookaheadVal = 15;
 
              lookaheadDepthLocal =
                  layers.size() < 8 || nrQubits <= 10 ? 0
                  : layers.size() < 15 ? static_cast<int>(lookaheadVal)
                  : layers.size() < 25 ? static_cast<int>(1.5 * lookaheadVal)
                                       : 2 * lookaheadVal;
            }
 
            int lookaheadHeuristicDepthLocal = lookaheadDepthWithHeuristic;
 
            if (lookaheadHeuristicDepthLocal == std::numeric_limits<int>::max())
              lookaheadHeuristicDepthLocal =
                  layers.size() < 10 || nrQubits <= 10 ? 0
                  : layers.size() < 20 ? lookaheadDepthLocal - 1
                                       : lookaheadDepthLocal - 2;
            
            if (lookaheadHeuristicDepthLocal < 0)
              lookaheadHeuristicDepthLocal = 0;
 
            sim->setGrowthFactorGate(growthFactorGate);
            sim->setGrowthFactorSwap(growthFactorSwap);
            sim->SetUseOptimalMeetingPosition(true);
            sim->SetLookaheadDepth(lookaheadDepthLocal);
            sim->SetLookaheadDepthWithHeuristic(lookaheadHeuristicDepthLocal);
            sim->SetUpcomingGates(dcirc->GetOperations());
          }
        }
      }
    }
  }

  void ConvertBackState() {
    auto optimiser = controller->GetOptimiser();
    if (optimiser) {
      // convert the classical state results back to the expected order
      const auto &qubitsMap = optimiser->GetReverseQubitsMap();
 
      ConvertBackState(qubitsMap);
    }
  }

  void ConvertBackState(
      const std::unordered_map<Types::qubit_t, Types::qubit_t> &qubitsMap) {
    // might not be the one stored in the network, might exist in the DES
    Circuits::OperationState &theClassicalState = GetState();
 
    theClassicalState.Remap(qubitsMap);
  }

  void ConvertBackResults(ExecuteResults &res) {
    auto optimiser = controller->GetOptimiser();
    if (optimiser) {
      // convert the classical state results back to the expected order
      const auto &qubitsMap = optimiser->GetReverseQubitsMap();
 
      ConvertBackResults(res, qubitsMap);
    }
  }

  void ConvertBackResults(
      ExecuteResults &res,
      const std::unordered_map<Types::qubit_t, Types::qubit_t> &bitsMap) const {
    ExecuteResults translatedRes;
 
    size_t numClassicalBits = 0;
    for (const auto &[q, b] : bitsMap)
      if (b >= numClassicalBits) numClassicalBits = b + 1;
 
    numClassicalBits = std::max(numClassicalBits, GetNumClassicalBits());
 
    for (const auto &r : res) {
      Circuits::OperationState translatedState(r.first);
 
      translatedState.Remap(bitsMap, false, numClassicalBits);
      translatedRes[translatedState.GetAllBits()] = r.second;
    }
 
    res.swap(translatedRes);
  }

  std::unordered_map<Types::qubit_t, Types::qubit_t> MapCircuitOnHost(
      const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t hostId,
      size_t &nrQubits, size_t &nrCbits, bool useSeparateSimForHosts = false) {
    qubitsMapOnHost.clear();
    nrQubits = 0;
    nrCbits = 0;
    if (!circuit) return {};
 
    const auto host =
        std::static_pointer_cast<SimpleHost<Time>>(GetHost(hostId));
    const size_t hostNrQubits = host->GetNumQubits();
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> reverseQubitsMap;
 
    if (!useSeparateSimForHosts) {
      size_t mxq = 0;
      size_t mnq = std::numeric_limits<size_t>::max();
      size_t mxb = 0;
      size_t mnb = std::numeric_limits<size_t>::max();
 
      for (const auto &op : circuit->GetOperations()) {
        const auto qbits = op->AffectedQubits();
        for (auto q : qbits) {
          if (q > mxq) mxq = q;
          if (q < mnq) mnq = q;
        }
        const auto cbits = op->AffectedBits();
        for (auto b : cbits) {
          if (b > mxb) mxb = b;
          if (b < mnb) mnb = b;
        }
      }
 
      if (mnq > mxq) mnq = 0;
      if (mnb > mxb) mnb = 0;
 
      nrQubits = mxq - mnq + 1;
      nrCbits = mxb - mnb + 1;
      if (nrCbits < nrQubits) nrCbits = nrQubits;
 
      const size_t startQubit = host->GetStartQubitId();
 
      if (mnq < startQubit || mxq >= startQubit + hostNrQubits) {
        if (nrQubits >
            hostNrQubits +
                1)  // the host has an additional 'special' qubit for the
                    // entanglement or other operations (like those for cutting)
          throw std::runtime_error("Circuit does not fit on the host!");
 
        for (size_t i = 0; i < nrCbits; ++i) {
          const size_t mapFrom = mnq + i;
          const size_t mapTo = startQubit + i;
 
          qubitsMapOnHost[mapFrom] = mapTo;
          reverseQubitsMap[mapTo] = mapFrom;
        }
 
        distCirc = std::static_pointer_cast<Circuits::Circuit<Time>>(
            circuit->Remap(qubitsMapOnHost, qubitsMapOnHost));
      }
 
      return reverseQubitsMap;
    }
 
    distCirc = circuit->RemapToContinuous(qubitsMapOnHost, reverseQubitsMap,
                                          nrQubits, nrCbits);
 
    if (simulator && Simulators::IsDistributedGpuSimulator(simulator->GetType())) {
      // Distribution settings refer to register qubits. Keep their numbering
      // and idle wires when the network creates a smaller per-host simulator.
      // Classical results retain RemapToContinuous's independent mapping.
      if (!hostNrQubits)
        throw std::runtime_error("Circuit does not fit on a host with no qubits!");
      const size_t start = host->GetStartQubitId();
      bool fitsLocal = true, fitsGlobal = true;
      for (const auto& entry : qubitsMapOnHost) {
        const auto q = entry.first;
        fitsLocal = fitsLocal && q < hostNrQubits;
        fitsGlobal = fitsGlobal && q >= start && q - start < hostNrQubits;
      }
      // Auto follows the host API's already-mapped precedence. Sparse local
      // circuits in the overlap must explicitly select local indexing.
      const bool global = distributedHostQubitIndexing == "global" ||
          (distributedHostQubitIndexing == "auto" && fitsGlobal);
      if (!(global ? fitsGlobal : fitsLocal))
        throw std::runtime_error("Circuit does not fit on the host!");
      const size_t offset = global ? start : 0;
      if (pauliStrings)
        for (const auto& pauli : *pauliStrings)
          if (pauli.size() > hostNrQubits)
            throw std::invalid_argument(
                "Host Pauli strings use local indices and must fit the host");
      std::unordered_map<Types::qubit_t, Types::qubit_t> restoreQubits;
      for (const auto& [original, compact] : qubitsMapOnHost) {
        if (original < offset || original - offset >= hostNrQubits)
          throw std::runtime_error("Circuit does not fit on the host!");
        restoreQubits[compact] = original - offset;
      }
      distCirc = std::static_pointer_cast<Circuits::Circuit<Time>>(
          distCirc->Remap(restoreQubits));
      qubitsMapOnHost.clear();
      for (size_t q = 0; q < hostNrQubits; ++q) qubitsMapOnHost[q] = q;
      nrQubits = hostNrQubits;
    } else if (pauliStrings) {
      // Observables can mention idle qubits absent from the circuit. Allocate
      // them in |0> and include them in the expectation remapping.
      size_t width = 0;
      for (const auto& pauli : *pauliStrings) width = std::max(width, pauli.size());
      for (size_t q = 0; q < width; ++q)
        if (!qubitsMapOnHost.count(q)) qubitsMapOnHost[q] = nrQubits++;
    }
 
    assert(nrQubits == qubitsMapOnHost.size());
 
    if (nrQubits == 0) nrQubits = 1;
 
    if (nrQubits >
        hostNrQubits +
            1)  // the host has an additional 'special' qubit for the
                // entanglement or other operations (like those for cutting)
      throw std::runtime_error("Circuit does not fit on the host!");
 
    return reverseQubitsMap;
  }
 
  bool optimizeSimulator = true; 
  typename BaseClass::SimulatorsSet simulatorsForOptimizations;
 
  int lastGpuDevice = -1;
  Simulators::SimulatorType lastSimulatorType =
      Simulators::SimulatorType::kQCSim; 
  Simulators::SimulationType lastMethod =
      Simulators::SimulationType::kStatevector; 
 
  Configuration<Time> configuration;
  std::string distributedHostQubitIndexing = "auto";
  std::unordered_map<Simulators::SimulatorType, std::string>
      resolvedDistributedDevices;
 
  size_t maxSimulators = QC::QubitRegisterCalculator<>::
      GetNumberOfThreads(); 
 
  Circuits::OperationState
      classicalState; 
  std::shared_ptr<Simulators::ISimulator>
      simulator; 
 
  std::shared_ptr<Circuits::Circuit<Time>>
      distCirc; 
 
  std::shared_ptr<IController<Time>>
      controller; 
  // TODO: depending on the network topology, we will have adiacency lists, etc.
  // or simply a vector of hosts for a totally connected network (or where the
  // communication details do not matter so much)
  std::vector<std::shared_ptr<IHost<Time>>>
      hosts; 
 
  std::unique_ptr<Estimators::SimulatorsEstimatorInterface<Time>>
      simulatorsEstimator; 
 
 private:
  Utils::ThreadsPool<ExecuteJob<Time>>
      threadsPool; 
  bool recreateIfNeeded =
      true; 
  std::unordered_map<Types::qubit_t, Types::qubit_t>
      qubitsMapOnHost; 
  const std::vector<std::string> *pauliStrings =
      nullptr; 
 
  bool optimizeInitialQubitsMap = true; 
  bool mpsOptimizeSwaps = true; 
  size_t mpsOptimizationBondDimensionThreshold =
      32; 
  size_t mpsOptimizationQubitsNumberThreshold =
      12; 
 
  int lookaheadDepth =
      std::numeric_limits<int>::max(); 
  int lookaheadDepthWithHeuristic = std::numeric_limits<int>::max(); 
 
  double growthFactorSwap = 1.;
  double growthFactorGate = 0.7;
  size_t curMaxBondDim = 0;
};
 
}  // namespace Network
 
#endif  // !_SIMPLE_NETWORK_H_