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"
 
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 = {}) {
    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__
        && simType != Simulators::SimulatorType::kGpuSim
#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();
 
    maxBondDim =
        simulator->GetConfiguration("matrix_product_state_max_bond_dimension");
    singularValueThreshold = simulator->GetConfiguration(
        "matrix_product_state_truncation_threshold");
    mpsSample = simulator->GetConfiguration("mps_sample_measure_algorithm");
 
    // 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();
    }
 
    std::vector<bool> executed;
    auto optSim =
        ChooseBestSimulator(distCirc, shots, nrQubits, nrQubits, nrCbitsResults,
                            simType, method, executed);
 
    lastSimulatorType = simType;
    lastMethod = method;
 
    size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
    if (simType == Simulators::SimulatorType::kGpuSim)
      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);
 
      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->maxBondDim = maxBondDim;
        job->mpsSample = mpsSample;
        job->singularValueThreshold = singularValueThreshold;
 
        job->network = BaseClass::getptr();
 
        if (optSim) {
          job->optSim = optSim->Clone();
          job->executedGates = executed;
        }
 
        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->maxBondDim = maxBondDim;
      job->mpsSample = mpsSample;
      job->singularValueThreshold = singularValueThreshold;
 
      job->network = BaseClass::getptr();
 
      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();
      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;
 
    std::shared_ptr<Circuits::Circuit<Time>> optCircuit;
    if (GetController()->GetOptimizeCircuit()) {
      optCircuit =
          std::static_pointer_cast<Circuits::Circuit<Time>>(circuit->Clone());
      optCircuit->Optimize();
    }
    const auto reverseQubitsMap = MapCircuitOnHost(
        GetController()->GetOptimizeCircuit() ? optCircuit : circuit, hostId,
        nrQubits, nrCbits, true);
    if (nrCbits == 0) nrCbits = nrQubits;
 
    if (!simulator || !distCirc) return {};
 
    auto simType = simulator->GetType();
 
    maxBondDim =
        simulator->GetConfiguration("matrix_product_state_max_bond_dimension");
    singularValueThreshold = simulator->GetConfiguration(
        "matrix_product_state_truncation_threshold");
    mpsSample = simulator->GetConfiguration("mps_sample_measure_algorithm");
 
    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__
        && simType != Simulators::SimulatorType::kGpuSim
#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();
 
    std::vector<bool> executed;
    auto optSim = ChooseBestSimulator(distCirc, shots, nrQubits, nrCbits,
                                      nrCbits, simType, method, executed);
 
    lastSimulatorType = simType;
    lastMethod = method;
 
    size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
    if (simType == Simulators::SimulatorType::kGpuSim)
      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);
 
      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->maxBondDim = maxBondDim;
        job->mpsSample = mpsSample;
        job->singularValueThreshold = singularValueThreshold;
 
        job->network = BaseClass::getptr();
 
        if (optSim) {
          job->optSim = optSim->Clone();
          job->executedGates = executed;
        }
 
        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->maxBondDim = maxBondDim;
      job->mpsSample = mpsSample;
      job->singularValueThreshold = singularValueThreshold;
 
      job->network = BaseClass::getptr();
 
      if (optSim) {
        optSim->SetMultithreading(true);
        job->optSim = optSim;
        job->executedGates = executed;
      }
 
      job->DoWorkNoLock();
      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) {
      if (!maxBondDim.empty())
        simulator->Configure("matrix_product_state_max_bond_dimension",
                             maxBondDim.c_str());
      if (!singularValueThreshold.empty())
        simulator->Configure("matrix_product_state_truncation_threshold",
                             singularValueThreshold.c_str());
      if (!mpsSample.empty())
        simulator->Configure("mps_sample_measure_algorithm", mpsSample.c_str());
      if (useDoublePrecision) simulator->Configure("use_double_precision", "1");
 
      simulator->AllocateQubits(
          nrQubits == 0 ? GetNumQubits() + GetNumNetworkEntangledQubits()
                        : nrQubits);
      simulator->Initialize();
 
      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("matrix_product_state_max_bond_dimension") == key)
      maxBondDim = value;
    else if (std::string("matrix_product_state_truncation_threshold") == key)
      singularValueThreshold = value;
    else if (std::string("mps_sample_measure_algorithm") == key)
      mpsSample = value;
    else if (std::string("use_double_precision") == key)
      useDoublePrecision =
          (std::string("1") == value || std::string("true") == value);
    else if (std::string("max_simulators") == key)
      maxSimulators = std::stoull(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;
  }

  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->maxBondDim = maxBondDim;
    cloned->singularValueThreshold = singularValueThreshold;
    cloned->mpsSample = mpsSample;
 
    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) const override {
    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);
 
    if (OptimizationSimulatorExists(
            Simulators::SimulatorType::kQCSim,
            Simulators::SimulationType::kMatrixProductState) &&
        (nrQubits <= 4 || !maxBondDim.empty()))
      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.empty()))
      simulatorTypes.emplace_back(
          Simulators::SimulatorType::kQiskitAer,
          Simulators::SimulationType::kMatrixProductState);
#endif
 
#ifdef __linux__
    if (Simulators::SimulatorsFactory::IsGpuLibraryAvailable()) {
      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);
 
    if (simulatorTypes.empty())
      return nullptr;
    else if (simulatorTypes.size() == 1) {
      simType = simulatorTypes[0].first;
      method = simulatorTypes[0].second;
 
      std::shared_ptr<Simulators::ISimulator> sim =
          Simulators::SimulatorsFactory::CreateSimulator(simType, method);
      if (sim) {
        if (method == Simulators::SimulationType::kMatrixProductState) {
          if (!maxBondDim.empty())
            sim->Configure("matrix_product_state_max_bond_dimension",
                           maxBondDim.c_str());
          if (!singularValueThreshold.empty())
            sim->Configure("matrix_product_state_truncation_threshold",
                           singularValueThreshold.c_str());
          if (!mpsSample.empty())
            sim->Configure("mps_sample_measure_algorithm", mpsSample.c_str());
 
          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);
        }
        sim->SetMultithreading(multithreading || GetMaxSimulators() == 1);
 
        return sim;
      }
    }
 
    std::shared_ptr<Simulators::ISimulator> sim =
        simulatorsEstimator->ChooseBestSimulator(
            simulatorTypes, dcirc, counts, nrQubits, nrCbits, nrResultCbits,
            simType, method, executed, maxBondDim, singularValueThreshold,
            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);
      }
      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);
  }
 
 protected:
  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 =
            maxBondDim.empty() ? 0 : std::stoi(maxBondDim);
 
        if (maxBondDim.empty() ||
            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 (!maxBondDim.empty())
            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);
 
    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;
 
  Simulators::SimulatorType lastSimulatorType =
      Simulators::SimulatorType::kQCSim; 
  Simulators::SimulationType lastMethod =
      Simulators::SimulationType::kStatevector; 
 
  std::string maxBondDim;
  std::string singularValueThreshold;
  std::string mpsSample;
  bool useDoublePrecision = false;
 
  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;
};
 
}  // namespace Network
 
#endif  // !_SIMPLE_NETWORK_H_