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"
 
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) {
    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 = RepeatedExecuteOnHost(circuit, hostId, 1);
    pauliStrings = nullptr;
 
    recreateIfNeeded = recreate;
 
    // 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);
      }
    }
 
    if (recreate && (!simulator || simType != simulator->GetType() ||
                     method != simulator->GetSimulationType() ||
                     simulator->GetNumberOfQubits() != numQubits))
      CreateSimulator(simType, method);
 
    return expectations;
  }

  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);
 
#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 || distCirc->empty())
      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) {
      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;
#endif
    if (method == Simulators::SimulationType::kStatevector &&
        !distCirc->HasOpsAfterMeasurements())
      nrThreads = 1;
 
    nrThreads = std::min(nrThreads, std::max<size_t>(shots, 1ULL));
 
    std::mutex resultsMutex;
 
    const 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) // 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;
 
        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;
 
      if (optSim) {
        optSim->SetMultithreading(true);
        job->optSim = optSim;
        job->executedGates = executed;
      } else {
        if (simulator && method == saveMethod && simType == saveSimType) {
          optSim = simulator;
          job->optSim = optSim;
          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->empty())
      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 &&
            !distCirc->HasOpsAfterMeasurements())
      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;
 
        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;
 
      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());
 
      simulator->AllocateQubits(
          nrQubits == 0 ? GetNumQubits() + GetNumNetworkEntangledQubits()
                        : nrQubits);
      simulator->Initialize();
    }
  }

  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("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;
 
    if (GetSimulator())
      cloned->CreateSimulator(GetSimulator()->GetType(),
                              GetSimulator()->GetSimulationType());
 
    return cloned;
  }
 
  std::shared_ptr<Simulators::ISimulator>
  ChooseBestSimulator(const 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 && 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);
 
#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);
    }
#endif
 
    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->SetMultithreading(true);
        if (!dontRunCircuitStart)
          Estimators::SimulatorsEstimatorInterface<
              Time>::ExecuteUpToMeasurements(dcirc, nrQubits, nrCbits,
                                             nrResultCbits, sim, executed,
                                             multithreading);
 
        return sim;
      }
    }
 
    return simulatorsEstimator->ChooseBestSimulator(
        simulatorTypes, dcirc, counts, nrQubits, nrCbits, nrResultCbits,
        simType, method, executed, maxBondDim, singularValueThreshold,
        mpsSample, GetMaxSimulators(), pauliStrings, multithreading,
        dontRunCircuitStart);
  }
 
protected:
  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;
 
  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; 
};
 
} // namespace Network
 
#endif // !_SIMPLE_NETWORK_H_