SimpleDisconnectedNetwork.h

Go to the documentation of this file.

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

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

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

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

    void Execute(const std::shared_ptr<Circuits::Circuit<Time>> &circuit) override
    {
        const auto recreate = recreateIfNeeded;
 
        auto simType = Simulators::SimulatorType::kQCSim;
        auto method = Simulators::SimulationType::kMatrixProductState;
        size_t numQubits = 2;
        if (simulator)
        {
            simType = simulator->GetType();
            method = simulator->GetSimulationType();
            numQubits = simulator->GetNumberOfQubits();
        }
 
        // Execute continues from the caller's prepared quantum state. Repeated
        // shots and expectation queries still start from zero, even for one shot.
        const auto reset = resetBeforeExecution;
        recreateIfNeeded = false;
        resetBeforeExecution = false;
 
        ExecuteResults res;
        try
        {
            res = RepeatedExecute(circuit, 1);
        }
        catch (...)
        {
            recreateIfNeeded = recreate;
            resetBeforeExecution = reset;
            throw;
        }
 
        recreateIfNeeded = recreate;
        resetBeforeExecution = reset;
 
        // 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))))
            RecreateSimulatorKeepingResults(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))))
            RecreateSimulatorKeepingResults(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;
 
        const auto cleanCircuit = circuit ? circuit->RemoveEndMeasurements() : circuit;
 
        pauliStrings = &paulis;
        const auto res = RepeatedExecute(cleanCircuit, 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();
                std::vector<std::string> translatedPaulis;
                translatedPaulis.reserve(paulis.size());
 
                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];
                    }
 
                    translatedPaulis.push_back(std::move(translated));
                }
                expectations = simulator->ExpectationValues(translatedPaulis);
            }
            else
            {
                expectations = simulator->ExpectationValues(paulis);
            }
        }
 
        if (recreate &&
            (!simulator || simType != simulator->GetType() || method != simulator->GetSimulationType() || simulator->GetNumberOfQubits() != numQubits))
            RecreateSimulatorKeepingResults(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);
 
        const auto cleanCircuit = circuit ? circuit->RemoveEndMeasurements() : circuit;
 
        pauliStrings = &paulis;
        const auto res = RepeatedExecuteOnHost(cleanCircuit, 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
            std::vector<std::string> translatedPaulis;
            translatedPaulis.reserve(paulis.size());
            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;
 
                translatedPaulis.push_back(std::move(translated));
            }
            expectations = simulator->ExpectationValues(translatedPaulis);
        }
        else
        {
            throw std::runtime_error("ExecuteOnHostExpectations: no simulator available after execution.");
        }
 
        if (recreate &&
            (!simulator || simType != simulator->GetType() || method != simulator->GetSimulationType() || simulator->GetNumberOfQubits() != numQubits))
            RecreateSimulatorKeepingResults(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.");
 
        const size_t n = simulator->GetNumberOfQubits();
        auto amplitudes = simulator->GetStateVector();
        const size_t dim = amplitudes.size();
 
        // Restore ascending original-qubit order within the compact register.
        // Original indices may be sparse, so use their ranks as output positions.
        if (!qubitsMapOnHost.empty())
        {
            // Build reverse mapping: simulator qubit position -> rank of its
            // original qubit among the mapped ones.
            std::vector<std::pair<Types::qubit_t, Types::qubit_t>> byOriginal(qubitsMapOnHost.begin(), qubitsMapOnHost.end());
            std::sort(byOriginal.begin(), byOriginal.end());
 
            std::vector<size_t> simToOrig(n);
            for (size_t pos = 0; pos < byOriginal.size(); ++pos)
                simToOrig[byOriginal[pos].second] = pos;
 
            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))
            RecreateSimulatorKeepingResults(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))
            RecreateSimulatorKeepingResults(simType, method);
 
        return result;
    }

    ExecuteResults RepeatedExecute(const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t shots = 1000) override
    {
        if (!controller || !circuit)
            return {};
 
        distCirc = controller->DistributeCircuit(BaseClass::getptr(), circuit);
        if (!distCirc)
            return {};
 
#ifdef _DEBUG
        for (auto q : distCirc->AffectedQubits())
        {
            if (q >= GetNumQubits())
            {
                std::cout << "This is a distributed circuit, using entanglement or cutting" << std::endl;
                break;
            }
        }
#endif
 
        if (!simulator)
            return {};
 
        auto simType = simulator->GetType();
        if (distCirc->HasOpsAfterMeasurements() && (
#ifndef NO_QISKIT_AER
                                                       simType == Simulators::SimulatorType::kCompositeQiskitAer ||
#endif
                                                       simType == Simulators::SimulatorType::kCompositeQCSim))
            distCirc->MoveMeasurementsAndResets();
 
        auto method = simulator->GetSimulationType();
 
        const auto saveSimType = simType;
        const auto saveMethod = method;
 
        if (GetOptimizeSimulator() && distCirc->IsClifford() && method != Simulators::SimulationType::kStabilizer
        // this is for the gpu simulator, as it doesn't support stabilizer
#ifdef __linux__
            && !Simulators::IsGpuSimulator(simType)
#endif
        )
        {
            method = Simulators::SimulationType::kStabilizer;
 
            if (simType == Simulators::SimulatorType::kCompositeQCSim)
                simType = Simulators::SimulatorType::kQCSim;
#ifndef NO_QISKIT_AER
            else if (simType == Simulators::SimulatorType::kCompositeQiskitAer)
                simType = Simulators::SimulatorType::kQiskitAer;
#endif
        }
 
        ExecuteResults res;
        const size_t nrQubits = GetNumQubits() + GetNumNetworkEntangledQubits();
        const size_t nrCbitsResults = GetNumClassicalBits();
        const size_t nrCbits = nrCbitsResults + GetNumNetworkEntangledQubits();
 
        CaptureSimulatorConfiguration();
 
        // do that only if the optimization for simulator is on and the estimator is
        // available, ortherwise an 'optimal' simulator won't be created
        if (optimizeSimulator && simulatorsEstimator && simulatorsEstimator->IsInitialized())
        {
            simulator->Clear();
            GetState().Clear();
        }
 
        curMaxBondDim = 0;
 
        std::vector<bool> executed;
        auto optSim = ChooseBestSimulator(distCirc, shots, nrQubits, nrCbits, nrCbitsResults, simType, method, executed);
 
        lastSimulatorType = simType;
        lastMethod = method;
        lastGpuDevice = -1;
        lastFusionWidth = 0;
        lastFusionEnabled = false;
 
        const bool reproducibleSampling = UsesReproducibleSampling(simType, method);
        size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
        if (Simulators::IsGpuSimulator(simType))
            nrThreads = 1;
        else
#endif
            if (((method == Simulators::SimulationType::kStatevector || method == Simulators::SimulationType::kPathIntegral ||
                  method == Simulators::SimulationType::kDensityMatrix) &&
                 !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 (reproducibleSampling && dcirc->HasOpsAfterMeasurements() &&
            (shots > 1 || (shots == 1 && configuration.GetConfiguration("trajectory_policy") == "shot_v1")))
        {
            ExecuteLogicalShotBlocks(dcirc, res, shots, nrThreads, nrQubits, nrCbits, nrCbitsResults, simType, method, optSim, executed, resultsMutex);
        }
        else if (nrThreads > 1)
        {
            // Clones are owned by outer shot workers; avoid constructing inner pools.
            if (optSim && method == Simulators::SimulationType::kPauliPropagator && optSim->GetType() == Simulators::SimulatorType::kQCSim)
                optSim->SetMultithreading(false);
            // since it's going to execute on multiple threads, free the memory from
            // the network's simulator and state, it's going to use other ones,
            // created in the threads if optimization already exists, it will be
            // cloned in the threads, otherwise a new one will be created in the
            // threads
            if (!optimizeSimulator || !simulatorsEstimator || !simulatorsEstimator->IsInitialized()) // otherwise it was already cleared
            {
                simulator->Clear();
                GetState().Clear();
            }
 
            const size_t cntPerThread = std::max<size_t>(shots / nrThreads, 1ULL);
 
            threadsPool.Resize(nrThreads);
            threadsPool.SetFinishLimit(shots);
 
            uint64_t jobStream = 0;
            while (shots > 0)
            {
                const size_t curCnt = std::min(cntPerThread, shots);
 
                shots -= curCnt;
 
                auto job = std::make_shared<ExecuteJob<Time>>(dcirc, res, curCnt, nrQubits, nrCbits, nrCbitsResults, simType, method, resultsMutex);
                job->optimiseMultipleShotsExecution = curCnt > 1 || GetOptimizeSimulator();
 
                job->network = BaseClass::getptr();
                job->onSimulatorReady = [this](const Simulators::ISimulator &actual) {
                    lastFusionWidth.store(actual.GetGateFusionMaxQubits());
                    lastFusionEnabled.store(actual.IsGateFusionEnabled());
                };
                job->curMaxBondDim = &curMaxBondDim;
 
                job->config = ExecutionConfiguration(simType, nrQubits);
                job->randomStream = jobStream++;
                if (configuration.IsSet("seed"))
                {
                    const uint64_t childSeed = Simulators::IState::DeriveSeed(std::stoull(configuration.GetConfiguration("seed")), job->randomStream);
                    job->config.SetConfiguration("seed", std::to_string(childSeed));
                }
 
                if (optSim)
                {
                    job->optSim = optSim->Clone();
                    job->executedGates = executed;
                    if (job->config.IsSet("seed"))
                        job->optSim->SetSeed(std::stoull(job->config.GetConfiguration("seed")));
                }
 
                threadsPool.AddRunJob(std::move(job));
            }
 
            // The workers stay for the next execution: creating them costs
            // more than a small batch of shots.
            threadsPool.WaitForFinish();
        }
        else
        {
            const size_t curCnt = shots;
 
            auto job = std::make_shared<ExecuteJob<Time>>(dcirc, res, curCnt, nrQubits, nrCbits, nrCbitsResults, simType, method, resultsMutex);
            job->optimiseMultipleShotsExecution = curCnt > 1 || GetOptimizeSimulator();
 
            job->network = BaseClass::getptr();
            job->onSimulatorReady = [this](const Simulators::ISimulator &actual) {
                lastFusionWidth.store(actual.GetGateFusionMaxQubits());
                lastFusionEnabled.store(actual.IsGateFusionEnabled());
            };
            job->curMaxBondDim = &curMaxBondDim;
 
            job->config = ExecutionConfiguration(simType, nrQubits);
 
            if (optSim)
            {
                optSim->SetMultithreading(true);
                job->optSim = optSim;
                job->executedGates = executed;
            }
            else
            {
                if (simulator && method == saveMethod && simType == saveSimType)
                {
                    // use the already created simulator
                    optSim = simulator;
                    job->optSim = optSim;
                    if (optSim->GetNumberOfQubits() == nrQubits)
                    {
                        if (resetBeforeExecution)
                        {
                            optSim->Reset();
                            optSim->SetGatesCounter(0);
                        }
                        job->config.ApplyConfigurationToSimulator(optSim);
                        // An initial mapping is only valid for a freshly reset register.
                        if (resetBeforeExecution)
                            OptimizeMPSInitialQubitsMap(optSim, dcirc, nrQubits);
                    }
                    job->executedGates.resize(dcirc->size(),
                                              false); // no gates executed yet
                    simulator = nullptr;
                }
            }
 
            job->DoWorkNoLock();
            lastGpuDevice = job->optSim ? job->optSim->GetGpuDevice() : -1;
            if (!recreateIfNeeded)
                simulator = job->optSim;
        }
 
        if (recreateIfNeeded)
            CreateSimulator(saveSimType, saveMethod);
 
        ConvertBackResults(res);
 
        return res;
    }

    ExecuteResults RepeatedExecuteOnHost(const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t hostId, size_t shots = 1000) override
    {
        if (!circuit || hostId >= GetNumHosts())
            return {};
 
        size_t nrQubits = 0;
        size_t nrCbits = 0;
 
        const bool distributed = simulator && Simulators::IsDistributedGpuSimulator(simulator->GetType());
        auto mappingCircuit = circuit;
        if (!distributed && GetController()->GetOptimizeCircuit())
        {
            mappingCircuit = std::static_pointer_cast<Circuits::Circuit<Time>>(circuit->Clone());
            mappingCircuit->Optimize();
        }
        const auto reverseQubitsMap = MapCircuitOnHost(mappingCircuit, hostId, nrQubits, nrCbits, true);
        // Resolve indexing before optimization can remove a wire that disambiguates
        // local from global numbering (for example a cancelling pair on qubit 0).
        if (distributed && distCirc && GetController()->GetOptimizeCircuit())
            distCirc->Optimize();
        if (nrCbits == 0)
            nrCbits = nrQubits;
 
        if (!simulator || !distCirc)
            return {};
 
        auto simType = simulator->GetType();
 
        CaptureSimulatorConfiguration();
 
        if (distCirc->HasOpsAfterMeasurements() && (
#ifndef NO_QISKIT_AER
                                                       simType == Simulators::SimulatorType::kCompositeQiskitAer ||
#endif
                                                       simType == Simulators::SimulatorType::kCompositeQCSim))
            distCirc->MoveMeasurementsAndResets();
 
        auto method = simulator->GetSimulationType();
        const auto saveSimType = simType;
        const auto saveMethod = method;
 
        if (GetOptimizeSimulator() && distCirc->IsClifford() && method != Simulators::SimulationType::kStabilizer
        // this is for the gpu simulator, as it doesn't support stabilizer
#ifdef __linux__
            && !Simulators::IsGpuSimulator(simType)
#endif
        )
        {
            method = Simulators::SimulationType::kStabilizer;
 
            if (simType == Simulators::SimulatorType::kCompositeQCSim)
                simType = Simulators::SimulatorType::kQCSim;
#ifndef NO_QISKIT_AER
            else if (simType == Simulators::SimulatorType::kCompositeQiskitAer)
                simType = Simulators::SimulatorType::kQiskitAer;
#endif
        }
 
        ExecuteResults res;
 
        // Backend selection can allocate a replacement. Fixed single-worker jobs
        // retain the existing simulator; release it below if workers are needed.
        if (GetOptimizeSimulator() && !distributed)
            simulator->Clear();
        GetState().Clear();
 
        curMaxBondDim = 0;
 
        std::vector<bool> executed;
        auto optSim = ChooseBestSimulator(distCirc, shots, nrQubits, nrCbits, nrCbits, simType, method, executed);
 
        lastSimulatorType = simType;
        lastMethod = method;
        lastGpuDevice = -1;
        lastFusionWidth = 0;
        lastFusionEnabled = false;
 
        const bool reproducibleSampling = UsesReproducibleSampling(simType, method);
        size_t nrThreads = GetMaxSimulators();
 
#ifdef __linux__
        if (Simulators::IsGpuSimulator(simType))
            nrThreads = 1;
        else
#endif
            if (((method == Simulators::SimulationType::kStatevector || method == Simulators::SimulationType::kPathIntegral ||
                  method == Simulators::SimulationType::kDensityMatrix) &&
                 !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;
 
        auto dcirc = distCirc;
 
        if (reproducibleSampling && dcirc->HasOpsAfterMeasurements() &&
            (shots > 1 || (shots == 1 && configuration.GetConfiguration("trajectory_policy") == "shot_v1")))
        {
            ExecuteLogicalShotBlocks(dcirc, res, shots, nrThreads, nrQubits, nrCbits, nrCbits, simType, method, optSim, executed, resultsMutex);
        }
        else if (nrThreads > 1)
        {
            // Clones are owned by outer shot workers; avoid constructing inner pools.
            if (optSim && method == Simulators::SimulationType::kPauliPropagator && optSim->GetType() == Simulators::SimulatorType::kQCSim)
                optSim->SetMultithreading(false);
            simulator->Clear();
            // this rounds up, rounding down is better
            // const size_t cntPerThread = static_cast<size_t>((shots - 1) / nrThreads
            // + 1);
            const size_t cntPerThread = std::max<size_t>(shots / nrThreads, 1ULL);
 
            threadsPool.Resize(nrThreads);
            threadsPool.SetFinishLimit(shots);
 
            uint64_t jobStream = 0;
            while (shots > 0)
            {
                const size_t curCnt = std::min(cntPerThread, shots);
                shots -= curCnt;
 
                auto job = std::make_shared<ExecuteJob<Time>>(dcirc, res, curCnt, nrQubits, nrCbits, nrCbits, simType, method, resultsMutex);
                job->optimiseMultipleShotsExecution = curCnt > 1 || GetOptimizeSimulator();
 
                job->network = BaseClass::getptr();
                job->onSimulatorReady = [this](const Simulators::ISimulator &actual) {
                    lastFusionWidth.store(actual.GetGateFusionMaxQubits());
                    lastFusionEnabled.store(actual.IsGateFusionEnabled());
                };
                job->curMaxBondDim = &curMaxBondDim;
 
                job->config = ExecutionConfiguration(simType, nrQubits);
                job->randomStream = jobStream++;
                if (configuration.IsSet("seed"))
                {
                    const uint64_t childSeed = Simulators::IState::DeriveSeed(std::stoull(configuration.GetConfiguration("seed")), job->randomStream);
                    job->config.SetConfiguration("seed", std::to_string(childSeed));
                }
 
                if (optSim)
                {
                    job->optSim = optSim->Clone();
                    job->executedGates = executed;
                    if (job->config.IsSet("seed"))
                        job->optSim->SetSeed(std::stoull(job->config.GetConfiguration("seed")));
                }
 
                threadsPool.AddRunJob(std::move(job));
            }
 
            // The workers stay for the next execution: creating them costs
            // more than a small batch of shots.
            threadsPool.WaitForFinish();
        }
        else
        {
            const size_t curCnt = shots;
 
            auto job = std::make_shared<ExecuteJob<Time>>(dcirc, res, curCnt, nrQubits, nrCbits, nrCbits, simType, method, resultsMutex);
            job->optimiseMultipleShotsExecution = curCnt > 1 || GetOptimizeSimulator();
 
            job->network = BaseClass::getptr();
            job->onSimulatorReady = [this](const Simulators::ISimulator &actual) {
                lastFusionWidth.store(actual.GetGateFusionMaxQubits());
                lastFusionEnabled.store(actual.IsGateFusionEnabled());
            };
            job->curMaxBondDim = &curMaxBondDim;
 
            job->config = ExecutionConfiguration(simType, nrQubits);
 
            if (optSim)
            {
                optSim->SetMultithreading(true);
                job->optSim = optSim;
                job->executedGates = executed;
            }
            else if (simulator && method == saveMethod && simType == saveSimType)
            {
                optSim = simulator;
                job->optSim = optSim;
                if (optSim->GetNumberOfQubits() == nrQubits)
                {
                    // Host execution always starts from zero, even when the previous
                    // invocation retained its state for an expectation/state query.
                    optSim->Reset();
                    optSim->SetGatesCounter(0);
                    job->config.ApplyConfigurationToSimulator(optSim);
                    OptimizeMPSInitialQubitsMap(optSim, dcirc, nrQubits);
                }
                // DoWorkNoLock reallocates/reinitializes when the host width changed.
                job->executedGates.resize(dcirc->size(), false);
            }
 
            job->DoWorkNoLock();
            lastGpuDevice = job->optSim ? job->optSim->GetGpuDevice() : -1;
            if (!recreateIfNeeded)
                simulator = job->optSim;
        }
 
        if (recreateIfNeeded)
            CreateSimulator(saveSimType, saveMethod);
 
        if (!reverseQubitsMap.empty())
            ConvertBackResults(res, reverseQubitsMap);
 
        return res;
    }

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

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

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

    void Configure(const char *key, const char *value) override
    {
        if (!key || !value)
            return;
 
        if (std::string("trajectory_policy") == key)
        {
            const std::string policy(value);
            if (policy != "block_v1" && policy != "shot_v1")
                throw std::invalid_argument("trajectory_policy must be block_v1 or shot_v1");
        }
        if (std::string("trajectory_max_memory_mb") == key)
        {
            size_t mb = 0;
            const std::string_view text(value);
            const auto parsed = std::from_chars(text.data(), text.data() + text.size(), mb);
            if (parsed.ec != std::errc{} || parsed.ptr != text.data() + text.size() || !mb || mb > (std::numeric_limits<size_t>::max() >> 20))
                throw std::invalid_argument("trajectory_max_memory_mb must be a positive representable MiB count");
        }
 
        if (std::string("gate_fusion") == key)
        {
            const std::string setting(value);
            if (setting != "true" && setting != "false" && setting != "0" && setting != "1" && setting != "auto")
                throw std::invalid_argument("gate_fusion must be true, false, 1, 0 or auto");
            configuration.SetConfiguration(key, value);
            if (simulator && (simulator->GetType() == Simulators::SimulatorType::kQCSim || simulator->GetType() == Simulators::SimulatorType::kCompositeQCSim ||
                              Simulators::IsGpuSimulator(simulator->GetType())))
                simulator->Configure(key, value);
            return;
        }
        if (std::string("distributed_host_qubit_indexing") == key)
        {
            const std::string mode(value);
            if (mode != "auto" && mode != "local" && mode != "global")
                throw std::invalid_argument("distributed_host_qubit_indexing must be auto, local or global");
            distributedHostQubitIndexing = mode;
            return;
        }
 
        if (std::string("distributed_devices") == key || std::string("gpu_device") == key)
            resolvedDistributedDevices.clear();
 
        if (std::string("max_simulators") == key)
            maxSimulators = std::stoull(value);
 
        if (std::string("gpu_device") == key)
        {
            Simulators::Configuration::ParseGpuDevice(value);
            if (simulator)
                simulator->Configure(key, value);
            configuration.SetConfiguration(key, value);
            return;
        }
        if (!Simulators::Configuration::GpuSvdSettingGroup(key).empty())
        {
            if (simulator)
                simulator->Configure(key, value);
            configuration.SetConfiguration(key, value);
            return;
        }
        configuration.SetConfiguration(key, value);
 
        if (simulator)
            simulator->Configure(key, value);
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    int GetLastGpuDevice() const override
    {
        return lastGpuDevice;
    }
 
    unsigned GetLastGateFusionMaxQubits() const override
    {
        return lastFusionWidth.load();
    }
 
    bool WasGateFusionEnabled() const override
    {
        return lastFusionEnabled.load();
    }
 
    bool GetGateFusion() const override
    {
        if (simulator)
        {
            const auto &options = simulator->GetConfigMap();
            const auto found = options.find("gate_fusion");
            if (found != options.end())
                return found->second != "false" && found->second != "0";
        }
        if (!configuration.IsSet("gate_fusion"))
            return true;
        const auto value = configuration.GetConfiguration("gate_fusion");
        return value != "false" && value != "0";
    }
 
    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 SetTrajectoryMinWorkerMicros(double micros)
    {
        trajectoryMinWorkerMicros = std::max(0., micros);
    }
 
    double GetTrajectoryMinWorkerMicros() const
    {
        return trajectoryMinWorkerMicros;
    }

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

    bool GetOptimizeSimulator() const override
    {
        return optimizeSimulator;
    }

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

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

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

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

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

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

    size_t GetCurrentMaxBondDimension() const override
    {
        return curMaxBondDim;
    }
 
  protected:
    // Resolved placement is separate from user settings. Importing a smaller
    // simulator must not turn its automatic shard subset into an explicit choice.
    void CaptureSimulatorConfiguration()
    {
        const auto requested = configuration.GetConfiguration("distributed_devices");
        configuration.ApplyConfigurationFromSimulator(simulator);
        if (simulator && Simulators::IsDistributedGpuSimulator(simulator->GetType()))
            configuration.SetConfiguration("distributed_devices", requested);
    }
 
    // Recreating the network simulator after a single execution must not drop
    // the classical results that execution stored in the network state.
    void RecreateSimulatorKeepingResults(Simulators::SimulatorType simType, Simulators::SimulationType method)
    {
        const auto results = GetState().GetAllBits();
        CreateSimulator(simType, method);
        GetState().SetResultsInOrder(results);
    }
 
    bool UsesReproducibleSampling(Simulators::SimulatorType type, Simulators::SimulationType method) const
    {
        return type == Simulators::SimulatorType::kQCSim && configuration.GetConfiguration("sampling_policy") != "legacy" &&
               (method == Simulators::SimulationType::kStatevector || method == Simulators::SimulationType::kDensityMatrix ||
                method == Simulators::SimulationType::kPathIntegral);
    }
 
    void ExecuteLogicalShotBlocks(const std::shared_ptr<Circuits::Circuit<Time>> &circuit, ExecuteResults &results, size_t shots, size_t workers, size_t qubits,
                                  size_t cbits, size_t resultBits, Simulators::SimulatorType type, Simulators::SimulationType method,
                                  const std::shared_ptr<Simulators::ISimulator> &prepared, const std::vector<bool> &executed, std::mutex &resultsMutex)
    {
        // Per-shot streams are the default. block_v1 retains the original fixed
        // logical ranges for replay. Terminal batches use one simulator instead.
        const bool independentShots = configuration.GetConfiguration("trajectory_policy") == "shot_v1";
        workers = std::max<size_t>(1, workers);
#ifdef _OPENMP
        // Outer jobs count against the same CPU budget as their internal teams.
        int threadBudget = omp_in_parallel() ? 1 : std::max(1, omp_get_max_threads());
        workers = std::min(workers, static_cast<size_t>(threadBudget));
#endif
        if (method == Simulators::SimulationType::kStatevector || method == Simulators::SimulationType::kDensityMatrix)
        {
            const size_t bitsPerQubit = method == Simulators::SimulationType::kDensityMatrix ? 2 : 1;
            if (qubits <= (std::numeric_limits<size_t>::digits - 5) / bitsPerQubit)
            {
                const size_t elements = size_t{1} << (bitsPerQubit * qubits);
                const size_t bytes = elements * sizeof(std::complex<double>);
                const size_t mb =
                    configuration.IsSet("trajectory_max_memory_mb") ? std::stoull(configuration.GetConfiguration("trajectory_max_memory_mb")) : 1024;
                const size_t slots = (mb << 20) / bytes;
                // Prepared executions retain a source and one shared checkpoint.
                // Without preparation each worker owns its own checkpoint too.
                const size_t memoryWorkers = prepared ? (slots > 2 ? slots - 2 : 1) : slots / 2;
                workers = std::min(workers, std::max<size_t>(1, memoryWorkers));
#ifdef _OPENMP
                // Large independent states contend for cache and memory bandwidth.
                // Reserve a small internal team instead of filling every CPU with
                // a separate state, and use no more threads in all than QCSim's
                // team for one such register (half the processors by default):
                // memory-bound work gains nothing from the second hardware
                // thread of a core, and with the default OpenMP limit filling
                // every logical CPU made these runs 1.3-1.5x slower. Logical RNG
                // streams do not depend on this choice.
                if (elements >= (size_t{1} << 20))
                {
                    threadBudget = std::min(threadBudget, std::max(1, QC::QubitRegisterCalculator<>::GetParallelTeamSize(elements)));
                    workers = std::min(workers, static_cast<size_t>(std::max(1, threadBudget / 4)));
                }
#endif
            }
            else
                workers = 1;
        }
        // shot_v1 streams do not depend on the worker count, so tiny states can
        // skip workers whose dispatch would cost more than their shots.
        if (independentShots && trajectoryMinWorkerMicros > 0)
            workers = std::min(workers, TrajectoryWorkersForWork(circuit, shots, qubits, method, executed));
        // shot_v1 can change physical job boundaries without changing any RNG
        // stream. Keep batches bounded while exposing short requests to workers.
        const size_t blockSize = independentShots ? std::min<size_t>(256, std::max<size_t>(1, shots / workers + (shots % workers != 0))) : 256;
        const uint64_t seed = configuration.IsSet("seed") ? std::stoull(configuration.GetConfiguration("seed"))
                                                          : Simulators::GenerateRandomSeed(type, configuration.GetConfigMap());
        const auto cloneMutex = std::make_shared<std::mutex>();
        // Execute() continues the live register when no optimizer replaced it.
        // A one-shot trajectory must not clear that state or invalidate handles
        // held by callers preparing and inspecting the register around Execute().
        const bool continueExecution = !resetBeforeExecution && shots == 1 && !prepared && simulator && simulator->GetType() == type &&
                                       simulator->GetSimulationType() == method && simulator->GetNumberOfQubits() == qubits;
        if (prepared)
            prepared->SetMultithreading(false);
        if (simulator && simulator != prepared && !continueExecution)
            simulator->Clear();
        workers = std::max<size_t>(1, std::min(workers, shots / blockSize + (shots % blockSize != 0)));
        if (workers > 1)
        {
            threadsPool.Resize(workers);
            threadsPool.SetFinishLimit(shots);
        }
        // The caller's OpenMP limit is a budget shared by concurrent shot jobs.
        // std::thread workers do not inherit the caller's OpenMP ICV, so pass
        // the budget explicitly and apply it only for the lifetime of each job.
        int internalThreads = 1;
#ifdef _OPENMP
        if (!omp_in_parallel() && method != Simulators::SimulationType::kPathIntegral)
            internalThreads = static_cast<int>(std::max<size_t>(1, static_cast<size_t>(threadBudget) / workers));
#endif
#if defined(_OPENMP) && defined(_MSC_VER) && !defined(__clang__) && !defined(_OPENMP_LLVM_RUNTIME)
        // Legacy MSVC OpenMP keeps this setting process-wide. Set it once before
        // starting jobs; per-worker save/restore races with other active jobs.
        struct RestoreThreadLimit
        {
            const int previous = omp_get_max_threads();
 
            ~RestoreThreadLimit()
            {
                omp_set_num_threads(previous);
            }
        } restoreThreadLimit;
 
        omp_set_num_threads(internalThreads);
#endif
        std::shared_ptr<ExecuteJob<Time>> last;
        for (size_t begin = 0, block = 0; begin < shots; ++block)
        {
            const size_t count = std::min(blockSize, shots - begin);
            auto job = std::make_shared<ExecuteJob<Time>>(circuit, results, count, qubits, cbits, resultBits, type, method, resultsMutex);
            // Even a one-shot tail must remove the same deterministic prefix:
            // classical instruction positions are part of shot_v1 stream IDs.
            job->optimiseMultipleShotsExecution = independentShots || count > 1 || GetOptimizeSimulator();
            // Reproducible measurement retains serial scan/reduction order;
            // gates and disjoint collapse updates may use the remaining budget.
#if defined(_OPENMP) && defined(_MSC_VER) && !defined(__clang__) && !defined(_OPENMP_LLVM_RUNTIME)
            job->internalThreads = 0; // The dispatcher owns the process-wide limit.
#else
            job->internalThreads = internalThreads;
#endif
            job->allowInternalMultithreading = internalThreads > 1;
            job->network = BaseClass::getptr();
            job->onSimulatorReady = [this](const Simulators::ISimulator &actual) {
                lastFusionWidth.store(actual.GetGateFusionMaxQubits());
                lastFusionEnabled.store(actual.IsGateFusionEnabled());
            };
            job->curMaxBondDim = &curMaxBondDim;
            job->config = ExecutionConfiguration(type, qubits);
            job->config.SetConfiguration("seed", std::to_string(Simulators::IState::DeriveSeed(seed, block)));
            job->randomStream = block;
            job->independentShotStreams = independentShots;
            job->trajectorySeed = Simulators::IState::DeriveSeed(seed, 0x53484f545f5631ULL);
            job->firstShot = begin;
            job->cloneSource = prepared;
            job->cloneMutex = cloneMutex;
            job->executedGates = executed;
            if (continueExecution)
            {
                job->optSim = simulator;
                job->config.ApplyConfigurationToSimulator(job->optSim);
                job->executedGates.assign(circuit->size(), false);
            }
            last = job;
            if (workers > 1)
                threadsPool.AddRunJob(std::move(job));
            else
                job->DoWork();
            begin += count;
        }
        if (workers > 1)
            threadsPool.WaitForFinish();
        if (!recreateIfNeeded && last)
            simulator = last->optSim;
    }
 
    // Rough serial cost of the shots of a statevector or density-matrix
    // trajectory batch (fixed per-shot cost plus the operations replayed after
    // the checkpoint over every amplitude), measured on QCSim; one worker per
    // trajectoryMinWorkerMicros of it.
    size_t TrajectoryWorkersForWork(const std::shared_ptr<Circuits::Circuit<Time>> &circuit, size_t shots, size_t qubits, Simulators::SimulationType method,
                                    const std::vector<bool> &executed) const
    {
        if (method != Simulators::SimulationType::kStatevector && method != Simulators::SimulationType::kDensityMatrix)
            return std::numeric_limits<size_t>::max();
        const size_t bitsPerQubit = method == Simulators::SimulationType::kDensityMatrix ? 2 : 1;
        if (qubits > (std::numeric_limits<size_t>::digits - 5) / bitsPerQubit)
            return std::numeric_limits<size_t>::max();
        const double elements = static_cast<double>(size_t{1} << (bitsPerQubit * qubits));
        const size_t done = static_cast<size_t>(std::count(executed.begin(), executed.end(), true));
        const double ops = static_cast<double>(circuit->size() > done ? circuit->size() - done : 1);
        const double micros = static_cast<double>(shots) * (6. + 0.4e-3 * ops * elements);
        return std::max<size_t>(1, static_cast<size_t>(micros / trajectoryMinWorkerMicros));
    }
 
    Configuration<Time> ExecutionConfiguration(Simulators::SimulatorType type, size_t qubits) const
    {
        auto result = configuration;
        if (type == Simulators::SimulatorType::kQCSim)
            result.SetConfiguration("reproducible_trajectory", distCirc && distCirc->HasOpsAfterMeasurements() ? "true" : "false");
        const auto found = resolvedDistributedDevices.find(type);
        if (!Simulators::IsDistributedGpuSimulator(type) || !configuration.GetConfiguration("distributed_devices").empty() ||
            found == resolvedDistributedDevices.end())
            return result;
        auto devices = found->second;
        // MPI shard count belongs to the communicator. Explicit global-qubit
        // settings also constrain the shard count and must be validated unchanged.
        if (type == Simulators::SimulatorType::kDistGpuSim && !configuration.IsSet("distributed_global_qubits"))
        {
            size_t count = 1;
            for (char c : devices)
                if (c == ',')
                    ++count;
            size_t bits = 0;
            for (size_t n = count; n > 1; n >>= 1)
                ++bits;
            while (count > 1 && bits >= qubits)
            {
                count /= 2;
                --bits;
            }
            size_t end = 0;
            for (size_t i = 0; i < count; ++i)
            {
                end = devices.find(',', end);
                if (end == std::string::npos)
                    break;
                if (i + 1 < count)
                    ++end;
            }
            devices = devices.substr(0, end);
        }
        result.SetConfiguration("distributed_devices", devices);
        return result;
    }
 
    void OptimizeMPSInitialQubitsMap(std::shared_ptr<Simulators::ISimulator> &sim, std::shared_ptr<Circuits::Circuit<Time>> &dcirc, size_t nrQubits) const
    {
        if ((sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductState ||
             sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator) &&
            (optimizeInitialQubitsMap || mpsOptimizeSwaps) && sim->SupportsMPSSwapOptimization() && !dcirc->HasCompositeOperations())
        {
            // an MPO site carries a two-qubit (operator) physical index, so its
            // bonds grow like those of an MPS twice as long
            const size_t routedQubits = sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator ? 2 * nrQubits : nrQubits;
            if (mpsOptimizationQubitsNumberThreshold <= routedQubits)
            {
                const auto maxBondDimValue =
                    configuration.GetConfigurationAsInt(sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator &&
                                                                configuration.IsSet("matrix_product_operator_max_bond_dimension")
                                                            ? "matrix_product_operator_max_bond_dimension"
                                                            : "matrix_product_state_max_bond_dimension");
 
                if (maxBondDimValue <= 0 || static_cast<int>(mpsOptimizationBondDimensionThreshold) <= maxBondDimValue)
                {
                    // need to be sure the circuit is correctly converted
                    dcirc->ConvertForCutting(); // convert the three qubit gates
                    auto layers = dcirc->ToMultipleQubitsLayersNoClone();
                    // Finalize source ordering before preparing fusion. Routing layers
                    // are only a cost model and must never replace executable operations.
                    auto ordered = Circuits::Circuit<Time>::LayersToCircuit(layers);
                    dcirc->SetOperations(ordered->GetOperations());
                    sim->SetUpcomingGates(dcirc->GetOperations());
                    if (sim->IsGateFusionEnabled())
                    {
                        Circuits::Circuit<Time> routingCircuit;
                        for (const auto &op : sim->GetUpcomingRoutingOperations())
                            if (op->GetType() == Circuits::OperationType::kGate)
                                routingCircuit.AddOperation(op);
                        layers = routingCircuit.ToMultipleQubitsLayersNoClone();
                    }
 
                    Simulators::MPSDummySimulator dummySim(nrQubits, sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator);
                    dummySim.SetSwapGatesRelabelQubits(sim->SwapGatesRelabelQubits());
                    dummySim.setGrowthFactorGate(growthFactorGate);
                    dummySim.setGrowthFactorSwap(growthFactorSwap);
 
                    if (maxBondDimValue > 0)
                    {
                        const bool virtualBonds = Simulators::RoutingConfig::VirtualBondsEnabled(configuration);
                        dummySim.SetMaxBondDimension(
                            virtualBonds ? Simulators::MPSDummySimulator::EffectiveMapMaxBondDimension(
                                               nrQubits, sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator, maxBondDimValue)
                                         : maxBondDimValue);
                        dummySim.SetRoutingCostModel(virtualBonds ? maxBondDimValue : 0, Simulators::RoutingConfig::TruncationWeight(configuration));
                    }
 
                    if (optimizeInitialQubitsMap)
                    {
                        const auto optimalMap = dummySim.ComputeOptimalQubitsMap(layers);
                        sim->SetInitialQubitsMap(optimalMap);
                    }
 
                    if (mpsOptimizeSwaps)
                    {
                        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;
                            }
                            if (!layers.empty())
                                avgTwoQubitGatesPerLayer /= layers.size();
 
                            int lookaheadVal = static_cast<int>(4. * avgTwoQubitGatesPerLayer);
                            if (lookaheadVal > 15)
                                lookaheadVal = 15;
 
                            lookaheadDepthLocal = layers.size() < 8 || routedQubits <= 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() < 8 || routedQubits <= 10 ? 0 : lookaheadDepthLocal - 1;
 
                        if (lookaheadHeuristicDepthLocal < 0)
                            lookaheadHeuristicDepthLocal = 0;
 
                        sim->setGrowthFactorGate(growthFactorGate);
                        sim->setGrowthFactorSwap(growthFactorSwap);
                        sim->SetUseOptimalMeetingPosition(true);
                        if (lookaheadDepth == std::numeric_limits<int>::max())
                            sim->SetAutomaticLookaheadDepth(lookaheadDepthLocal);
                        else
                            sim->SetLookaheadDepth(lookaheadDepthLocal);
                        sim->SetLookaheadDepthWithHeuristic(lookaheadHeuristicDepthLocal);
                        // The CPU window installed before layout planning has not
                        // changed. Its routing setters already refresh the plan.
                        if (sim->GetType() != Simulators::SimulatorType::kQCSim)
                            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)
        {
            // A circuit already using this host's qubit IDs still replaces the
            // previous circuit (and must be available on the first invocation).
            distCirc = circuit;
            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;
                }
            }
            for (auto b : circuit->GetBits())
            {
                if (b > mxb)
                    mxb = b;
                if (b < mnb)
                    mnb = b;
            }
 
            if (mnq > mxq)
                mnq = 0;
            if (mnb > mxb)
                mnb = 0;
 
            nrQubits = mxq - mnq + 1;
            nrCbits = mxb - mnb + 1;
            if (nrCbits < nrQubits)
                nrCbits = nrQubits;
 
            const size_t startQubit = host->GetStartQubitId();
 
            if (mnq < startQubit || mxq >= startQubit + hostNrQubits)
            {
                if (nrQubits > hostNrQubits + 1) // the host has an additional 'special' qubit for the
                                                 // entanglement or other operations (like those for cutting)
                    throw std::runtime_error("Circuit does not fit on the host!");
 
                for (size_t i = 0; i < nrCbits; ++i)
                {
                    const size_t mapFrom = mnq + i;
                    const size_t mapTo = startQubit + i;
 
                    qubitsMapOnHost[mapFrom] = mapTo;
                    reverseQubitsMap[mapTo] = mapFrom;
                }
 
                distCirc = std::static_pointer_cast<Circuits::Circuit<Time>>(circuit->Remap(qubitsMapOnHost, qubitsMapOnHost));
            }
 
            return reverseQubitsMap;
        }
 
        distCirc = circuit->RemapToContinuous(qubitsMapOnHost, reverseQubitsMap, nrQubits, nrCbits);
 
        if (simulator && Simulators::IsDistributedGpuSimulator(simulator->GetType()))
        {
            // Distribution settings refer to register qubits. Keep their numbering
            // and idle wires when the network creates a smaller per-host simulator.
            // Classical results retain RemapToContinuous's independent mapping.
            if (!hostNrQubits)
                throw std::runtime_error("Circuit does not fit on a host with no qubits!");
            const size_t start = host->GetStartQubitId();
            bool fitsLocal = true, fitsGlobal = true;
            for (const auto &entry : qubitsMapOnHost)
            {
                const auto q = entry.first;
                fitsLocal = fitsLocal && q < hostNrQubits;
                fitsGlobal = fitsGlobal && q >= start && q - start < hostNrQubits;
            }
            // Auto follows the host API's already-mapped precedence. Sparse local
            // circuits in the overlap must explicitly select local indexing.
            const bool global = distributedHostQubitIndexing == "global" || (distributedHostQubitIndexing == "auto" && fitsGlobal);
            if (!(global ? fitsGlobal : fitsLocal))
                throw std::runtime_error("Circuit does not fit on the host!");
            const size_t offset = global ? start : 0;
            if (pauliStrings)
                for (const auto &pauli : *pauliStrings)
                    if (pauli.size() > hostNrQubits)
                        throw std::invalid_argument("Host Pauli strings use local indices and must fit the host");
            std::unordered_map<Types::qubit_t, Types::qubit_t> restoreQubits;
            for (const auto &[original, compact] : qubitsMapOnHost)
            {
                if (original < offset || original - offset >= hostNrQubits)
                    throw std::runtime_error("Circuit does not fit on the host!");
                restoreQubits[compact] = original - offset;
            }
            distCirc = std::static_pointer_cast<Circuits::Circuit<Time>>(distCirc->Remap(restoreQubits));
            qubitsMapOnHost.clear();
            for (size_t q = 0; q < hostNrQubits; ++q)
                qubitsMapOnHost[q] = q;
            nrQubits = hostNrQubits;
        }
        else if (pauliStrings)
        {
            // Observables can mention idle qubits absent from the circuit. Allocate
            // them in |0> and include them in the expectation remapping.
            size_t width = 0;
            for (const auto &pauli : *pauliStrings)
                width = std::max(width, pauli.size());
            for (size_t q = 0; q < width; ++q)
                if (!qubitsMapOnHost.count(q))
                    qubitsMapOnHost[q] = nrQubits++;
        }
 
        assert(nrQubits == qubitsMapOnHost.size());
 
        if (nrQubits == 0)
            nrQubits = 1;
 
        if (nrQubits > hostNrQubits + 1) // the host has an additional 'special' qubit for the
                                         // entanglement or other operations (like those for cutting)
            throw std::runtime_error("Circuit does not fit on the host!");
 
        return reverseQubitsMap;
    }
 
    bool optimizeSimulator = true; 
    typename BaseClass::SimulatorsSet simulatorsForOptimizations;
 
    int lastGpuDevice = -1;
    std::atomic<unsigned> lastFusionWidth{0};
    std::atomic<bool> lastFusionEnabled{false};
    Simulators::SimulatorType lastSimulatorType = Simulators::SimulatorType::kQCSim;  
    Simulators::SimulationType lastMethod = Simulators::SimulationType::kStatevector; 
 
    Configuration<Time> configuration;
    std::string distributedHostQubitIndexing = "auto";
    std::unordered_map<Simulators::SimulatorType, std::string> resolvedDistributedDevices;
 
    size_t maxSimulators = QC::QubitRegisterCalculator<>::GetNumberOfThreads(); 
    double trajectoryMinWorkerMicros = 100.; 
 
    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; 
 
  protected:
    // Only Execute() preserves a prepared input state. Derived shot loops must
    // honor this independently of retaining the final simulator for queries.
    bool resetBeforeExecution = true;
 
    // Derived execution loops must retain the final state for state queries.
    bool recreateIfNeeded = true; 
 
  private:
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsMapOnHost; 
    const std::vector<std::string> *pauliStrings = nullptr;             
 
    bool optimizeInitialQubitsMap = true;              
    bool mpsOptimizeSwaps = true;                      
    size_t mpsOptimizationBondDimensionThreshold = 32; 
    size_t mpsOptimizationQubitsNumberThreshold = 12;  
 
    int lookaheadDepth = std::numeric_limits<int>::max();              
    int lookaheadDepthWithHeuristic = std::numeric_limits<int>::max(); 
 
    double growthFactorSwap = 1.;
    double growthFactorGate = 0.7;
    size_t curMaxBondDim = 0;
};
 
} // namespace Network
 
#endif // !_SIMPLE_NETWORK_H_