ImmediateComposite.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _IMMEDIATE_COMPOSITE_H_
#define _IMMEDIATE_COMPOSITE_H_
 
#ifdef INCLUDED_BY_FACTORY
 
#include "Individual.h"
 
#include <vector>
 
#include "../Core/Configuration.h"
 
namespace Simulators
{
 
// TODO: Maybe use the pimpl idiom
// https://en.cppreference.com/w/cpp/language/pimpl to hide the implementation
// for good but during development this should be good enough
namespace Private
{

class ImmediateCompositeSimulator : public ISimulator
{
  public:
    ImmediateCompositeSimulator(SimulatorType type =
#ifdef NO_QISKIT_AER
                                    Simulators::SimulatorType::kQCSim
#else
                                    Simulators::SimulatorType::kQiskitAer
#endif
                                )
        : type(type)
    {
        SeedAuxiliaryRng(samplingSeed);
        // just in case somebody tries to use a composite simulator composed of
        // composite simulators or a gpu simulator - this one would work, but until
        // we implement some supporting functionality with cuda, especially for
        // composite probably it shouldn't be used splitting and merging is slow and
        // it should be done in the videocard memory instead of using the cpu, if
        // the gpu simulators are used otherwise the transfer from videocard memory
        // and host memory back and forth is going to slow down the simulation quite
        // a bit
        if (
#ifndef NO_QISKIT_AER
            type != SimulatorType::kQiskitAer &&
#endif
            type != SimulatorType::kQCSim)
            type =
#ifdef NO_QISKIT_AER
                Simulators::SimulatorType::kQCSim;
#else
                Simulators::SimulatorType::kQiskitAer;
#endif
    }

    void Initialize() override
    {
        simulators.clear();
        qubitsMap.resize(nrQubits);
        const uint64_t initializationSeed = nextInitialization++ == 0 ? samplingSeed : DeriveSeed(samplingSeed, 0x494E4954ULL + nextInitialization);
 
        // start with one qubit simulators:
        for (size_t q = 0; q < nrQubits; ++q)
        {
            qubitsMap[q] = q;
            auto sim = std::make_unique<IndividualSimulator>(type);
            sim->AllocateQubits(1);
 
            // this is for tests when we default to matrix product state simulator
#if 0
                    sim->Configure("method", "statevector");
#endif
 
            sim->SetMultithreading(enableMultithreading);
            if (config.IsSet("precision"))
                sim->Configure("precision", config.GetConfiguration("precision").c_str());
            sim->Initialize();
            for (const auto &[key, value] : config.GetConfigMap())
                if (key != "seed" && key != "method")
                    sim->Configure(key.c_str(), value.c_str());
            sim->SetSeed(DeriveSeed(initializationSeed, q));
            sim->GetQubitsMap()[q] = 0;
            simulators[q] = std::move(sim);
        }
        nextId = nrQubits;
    }

    void Reset() override
    {
        Initialize();
    }

    void InitializeState(size_t num_qubits, std::vector<std::complex<double>> &amplitudes) override
    {
        throw std::runtime_error("CompositeSimulator::InitializeState not supported");
    }

    /*
    void InitializeState(size_t num_qubits, std::vector<std::complex<double>,
    avoid_init_allocator<std::complex<double>>>& amplitudes) override
    {
            throw std::runtime_error("CompositeSimulator::InitializeState not
    supported");
    }
    */

#ifndef NO_QISKIT_AER
    void InitializeState(size_t num_qubits, AER::Vector<std::complex<double>> &amplitudes) override
    {
        throw std::runtime_error("CompositeSimulator::InitializeState not supported");
    }
#endif

    void InitializeState(size_t num_qubits, Eigen::VectorXcd &amplitudes) override
    {
        throw std::runtime_error("CompositeSimulator::InitializeState not supported");
    }

    void Configure(const char *key, const char *value) override
    {
        Utils::Sampling::ValidateSetting(key, value);
        // don't allow chaning the method, it should stay statevector
        if (std::string("method") == key)
            return;
 
        if (std::string(key) == "precision" || std::string(key) == "use_double_precision")
        {
            Configuration::ParsePrecision(key, value);
            // Fixed once the children exist; a reapplied precision is ignored until
            // Clear, so later children cannot end up with a different one.
            if (!simulators.empty())
                return;
            config.SetConfiguration(key, value);
            return;
        }
 
        config.SetConfiguration(key, value);
 
        if (std::string("seed") == key)
        {
            const uint64_t seed = std::stoull(value);
            samplingSeed = seed;
            nextCloneStream = nextRestoreStream = nextInitialization = nextSplitStream = 0;
            SeedAuxiliaryRng(seed);
            for (auto &[id, simulator] : simulators)
                simulator->SetSeed(DeriveSeed(seed, id));
        }
        else
        {
            for (auto &[id, simulator] : simulators)
                simulator->Configure(key, value);
        }
    }

    std::string GetConfiguration(const char *key) const override
    {
        if (std::string(key) == "sampling_seed")
            return std::to_string(samplingSeed);
        if (std::string(key) == "sampling_policy" && type == SimulatorType::kQCSim && !config.IsSet(key))
            return "reproducible_v1";
        if (std::string(key) == "precision" || std::string(key) == "use_double_precision")
        {
            if (!simulators.empty())
                return simulators.begin()->second->GetConfiguration(key);
            const bool useDouble = type == SimulatorType::kQCSim || !config.IsSet("precision") || config.GetConfiguration("precision") == "double";
            return std::string(key) == "precision" ? (useDouble ? "double" : "single") : (useDouble ? "1" : "0");
        }
        if (simulators.empty())
            return "";
 
        return config.GetConfiguration(key);
    }

    size_t AllocateQubits(size_t num_qubits) override
    {
        if (!simulators.empty())
            return 0;
 
        const size_t oldNrQubits = nrQubits;
        nrQubits += num_qubits;
        return oldNrQubits;
    }

    size_t GetNumberOfQubits() const override
    {
        return nrQubits;
    }

    void Clear() override
    {
        simulators.clear();
        qubitsMap.clear();
        nrQubits = 0;
    }

    void SaveStateToInternalDestructive() override
    {
        for (auto &[id, simulator] : simulators)
            simulator->SaveStateToInternalDestructive();
    }

    void RestoreInternalDestructiveSavedState() override
    {
        for (auto &[id, simulator] : simulators)
            simulator->RestoreInternalDestructiveSavedState();
    }

    std::complex<double> AmplitudeRaw(Types::qubit_t outcome) override
    {
        std::complex<double> res = 1.0;
 
        for (auto &[id, simulator] : simulators)
            res *= simulator->AmplitudeRaw(simulator->ConvertOutcomeFromGlobal(outcome));
 
        return res;
    }

    size_t Measure(const Types::qubits_vector &qubits) override
    {
        if (qubits.size() > sizeof(size_t) * 8)
            std::cerr << "Warning: The number of qubits to measure is larger than the "
                         "number of bits in the size_t type, the outcome will be undefined"
                      << std::endl;
 
        size_t res = 0;
        size_t mask = 1ULL;
 
        DontNotify();
        for (Types::qubit_t qubit : qubits)
        {
            const bool outcome = GetSimulator(qubit)->Measure({qubit}) != 0;
            if (outcome)
                res |= mask;
            mask <<= 1;
 
            Split(qubit, outcome);
        }
        Notify();
 
        NotifyObservers(qubits);
 
        return res;
    }

    std::vector<bool> MeasureMany(const Types::qubits_vector &qubits) override
    {
        std::vector<bool> res(qubits.size(), false);
 
        DontNotify();
        for (size_t q = 0; q < qubits.size(); ++q)
        {
            Types::qubit_t qubit = qubits[q];
            const bool outcome = GetSimulator(qubit)->Measure({qubit}) != 0;
            if (outcome)
                res[q] = true;
 
            Split(qubit, outcome);
        }
        Notify();
 
        NotifyObservers(qubits);
 
        return res;
    }

    void ApplyReset(const Types::qubits_vector &qubits) override
    {
        DontNotify();
        for (Types::qubit_t qubit : qubits)
        {
            GetSimulator(qubit)->ApplyReset({qubit});
            Split(qubit, false);
        }
        Notify();
 
        NotifyObservers(qubits);
    }

    double Probability(Types::qubit_t outcome) override
    {
        return std::norm(Amplitude(outcome));
    }

    std::complex<double> Amplitude(Types::qubit_t outcome) override
    {
        std::complex<double> res = 1.0;
 
        for (auto &[id, simulator] : simulators)
            res *= simulator->Amplitude(outcome);
 
        return res;
    }

    std::complex<double> ProjectOnZero() override
    {
        std::complex<double> res = 1.0;
 
        for (auto &[id, simulator] : simulators)
            res *= simulator->ProjectOnZero();
 
        return res;
    }

    std::vector<double> AllProbabilities() override
    {
        const size_t nrBasisStates = 1ULL << nrQubits;
        std::vector<double> result;
 
        for (size_t i = 0; i < nrBasisStates; ++i)
            result.emplace_back(Probability(i));
 
        return result;
    }

    std::vector<double> Probabilities(const Types::qubits_vector &qubits) override
    {
        std::vector<double> result;
 
        for (size_t i = 0; i < qubits.size(); ++i)
            result.emplace_back(Probability(qubits[i]));
 
        return result;
    }

    std::unordered_map<Types::qubit_t, Types::qubit_t> SampleCounts(const Types::qubits_vector &qubits, size_t shots = 1000) override
    {
        if (type == SimulatorType::kQCSim && config.GetConfiguration("sampling_policy") != "legacy")
            return SamplePrepared<false>(qubits, shots);
        if (qubits.empty() || !shots)
            return {};
        if (GetNumberOfQubits() > sizeof(Types::qubit_t) * 8)
            std::cerr << "Warning: The number of qubits to measure is larger than the "
                         "number of bits in the Types::qubit_t type, the outcome will be "
                         "undefined"
                      << std::endl;
        // TODO: improve it as for the qcsim statevector simulator case!
        std::unordered_map<Types::qubit_t, Types::qubit_t> result;
        DontNotify();
 
        SaveStateToInternalDestructive();
 
        if (shots > 1)
        {
            InitializeAlias();
 
            for (size_t shot = 0; shot < shots; ++shot)
            {
                size_t measRaw = 0;
                for (auto &[id, simulator] : simulators)
                    measRaw |= simulator->SampleFromAlias();
 
                size_t meas = 0;
                size_t mask = 1ULL;
                for (auto q : qubits)
                {
                    const size_t qubitMask = 1ULL << q;
                    if ((measRaw & qubitMask) != 0)
                        meas |= mask;
                    mask <<= 1ULL;
                }
 
                ++result[meas];
            }
 
            ClearAlias();
        }
        else
            for (size_t shot = 0; shot < shots; ++shot)
            {
                const auto measRaw = MeasureNoCollapse();
 
                size_t meas = 0;
                size_t mask = 1ULL;
                for (auto q : qubits)
                {
                    const size_t qubitMask = 1ULL << q;
                    if ((measRaw & qubitMask) != 0)
                        meas |= mask;
                    mask <<= 1ULL;
                }
 
                ++result[meas];
            }
 
        RestoreInternalDestructiveSavedState();
 
        Notify();
        NotifyObservers(qubits);
 
        return result;
    }

    std::unordered_map<std::vector<bool>, Types::qubit_t> SampleCountsMany(const Types::qubits_vector &qubits, size_t shots = 1000) override
    {
        if (type == SimulatorType::kQCSim && config.GetConfiguration("sampling_policy") != "legacy")
            return SamplePrepared<true>(qubits, shots);
        if (qubits.empty() || !shots)
            return {};
        std::unordered_map<std::vector<bool>, Types::qubit_t> result;
        DontNotify();
 
        SaveStateToInternalDestructive();
        std::vector<bool> meas(qubits.size());
        if (shots > 1)
        {
            InitializeAlias();
 
            for (size_t shot = 0; shot < shots; ++shot)
            {
                size_t measRaw = 0;
                for (auto &[id, simulator] : simulators)
                    measRaw |= simulator->SampleFromAlias();
 
                for (size_t i = 0; i < qubits.size(); ++i)
                    meas[i] = ((measRaw >> qubits[i]) & 1) == 1;
 
                ++result[meas];
            }
 
            ClearAlias();
        }
        else
        {
            for (size_t shot = 0; shot < shots; ++shot)
            {
                const auto measRaw = MeasureNoCollapseMany();
 
                for (size_t i = 0; i < qubits.size(); ++i)
                    meas[i] = measRaw[qubits[i]];
 
                ++result[meas];
            }
        }
 
        RestoreInternalDestructiveSavedState();
 
        Notify();
        NotifyObservers(qubits);
 
        return result;
    }

    double ExpectationValue(const std::string &pauliString) override
    {
        std::unordered_map<size_t, std::string> pauliStrings;
 
        for (size_t q = 0; q < pauliString.size(); ++q)
        {
            const char op = toupper(pauliString[q]);
            if (op == 'I')
                continue;
 
            const size_t simId = qubitsMap[q];
            const size_t localQubit = simulators[simId]->GetQubitsMap().at(q);
 
            if (pauliStrings[simId].size() <= localQubit)
                pauliStrings[simId].resize(localQubit + 1, 'I');
 
            pauliStrings[simId][localQubit] = op;
        }
 
        double result = 1.0;
        for (auto &[id, localPauliString] : pauliStrings)
            result *= simulators[id]->ExpectationValue(localPauliString);
 
        return result;
    }

    SimulatorType GetType() const override
    {
#ifdef NO_QISKIT_AER
        return SimulatorType::kCompositeQCSim;
#else
        if (type != SimulatorType::kQiskitAer)
            return SimulatorType::kCompositeQCSim;
 
        return SimulatorType::kCompositeQiskitAer;
#endif
    }

    SimulationType GetSimulationType() const override
    {
        return SimulationType::kStatevector;
    }

    void Flush() override
    {
        for (auto &[id, simulator] : simulators)
            simulator->Flush();
    }

    void ApplyGenericOneQubitGate(Types::qubit_t qubit, const Eigen::Matrix2cd &gate) override
    {
        GetSimulator(qubit)->ApplyGenericOneQubitGate(qubit, gate);
        NotifyObservers({qubit});
    }

    void ApplyGenericTwoQubitGate(Types::qubit_t qubit0, Types::qubit_t qubit1, const Eigen::Matrix4cd &gate) override
    {
        JoinIfNeeded(qubit0, qubit1);
        GetSimulator(qubit0)->ApplyGenericTwoQubitGate(qubit0, qubit1, gate);
        NotifyObservers({qubit0, qubit1});
    }
 
    void ApplyGenericThreeQubitGate(Types::qubit_t q0, Types::qubit_t q1, Types::qubit_t q2, const Matrix8cd &gate) override
    {
        JoinIfNeeded(q0, q1);
        JoinIfNeeded(q0, q2);
        GetSimulator(q0)->ApplyGenericThreeQubitGate(q0, q1, q2, gate);
        NotifyObservers({q0, q1, q2});
    }
 
    // YES, all one qubit gates are that easy:

    void ApplyP(Types::qubit_t qubit, double lambda) override
    {
        GetSimulator(qubit)->ApplyP(qubit, lambda);
        NotifyObservers({qubit});
    }

    void ApplyX(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyX(qubit);
        NotifyObservers({qubit});
    }

    void ApplyY(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyY(qubit);
        NotifyObservers({qubit});
    }

    void ApplyZ(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyZ(qubit);
        NotifyObservers({qubit});
    }

    void ApplyH(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyH(qubit);
        NotifyObservers({qubit});
    }

    void ApplyS(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyS(qubit);
        NotifyObservers({qubit});
    }

    void ApplySDG(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplySDG(qubit);
        NotifyObservers({qubit});
    }

    void ApplyT(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyT(qubit);
        NotifyObservers({qubit});
    }

    void ApplyTDG(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyTDG(qubit);
        NotifyObservers({qubit});
    }

    void ApplySx(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplySx(qubit);
        NotifyObservers({qubit});
    }

    void ApplySxDAG(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplySxDAG(qubit);
        NotifyObservers({qubit});
    }

    void ApplyK(Types::qubit_t qubit) override
    {
        GetSimulator(qubit)->ApplyK(qubit);
        NotifyObservers({qubit});
    }

    void ApplyRx(Types::qubit_t qubit, double theta) override
    {
        GetSimulator(qubit)->ApplyRx(qubit, theta);
        NotifyObservers({qubit});
    }

    void ApplyRy(Types::qubit_t qubit, double theta) override
    {
        GetSimulator(qubit)->ApplyRy(qubit, theta);
        NotifyObservers({qubit});
    }

    void ApplyRz(Types::qubit_t qubit, double theta) override
    {
        GetSimulator(qubit)->ApplyRz(qubit, theta);
        NotifyObservers({qubit});
    }

    void ApplyU(Types::qubit_t qubit, double theta, double phi, double lambda, double gamma) override
    {
        GetSimulator(qubit)->ApplyU(qubit, theta, phi, lambda, gamma);
        NotifyObservers({qubit});
    }
 
    // the gates that operate on more than one qubit need joining of the
    // simulators, if the qubits are in different simulators

    void ApplyCX(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCX(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCY(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCY(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCZ(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCZ(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCP(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double lambda) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCP(ctrl_qubit, tgt_qubit, lambda);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCRx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCRx(ctrl_qubit, tgt_qubit, theta);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCRy(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCRy(ctrl_qubit, tgt_qubit, theta);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCRz(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCRz(ctrl_qubit, tgt_qubit, theta);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCH(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCH(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCSx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCSx(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplyCSxDAG(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCSxDAG(ctrl_qubit, tgt_qubit);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }

    void ApplySwap(Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        JoinIfNeeded(qubit0, qubit1);
        GetSimulator(qubit0)->ApplySwap(qubit0, qubit1);
        NotifyObservers({qubit1, qubit0});
    }

    void ApplyCCX(Types::qubit_t qubit0, Types::qubit_t qubit1, Types::qubit_t qubit2) override
    {
        // TODO: See if it's worth optimizing to joing all three simulators at once,
        // if needed (that is, there is a different one for each qubit)
        JoinIfNeeded(qubit0, qubit1);
        JoinIfNeeded(qubit0, qubit2);
        GetSimulator(qubit0)->ApplyCCX(qubit0, qubit1, qubit2);
        NotifyObservers({qubit2, qubit1, qubit0});
    }

    void ApplyCSwap(Types::qubit_t ctrl_qubit, Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        // TODO: See if it's worth optimizing to joing all three simulators at once,
        // if needed (that is, there is a different one for each qubit)
        JoinIfNeeded(ctrl_qubit, qubit0);
        JoinIfNeeded(ctrl_qubit, qubit1);
        GetSimulator(ctrl_qubit)->ApplyCSwap(ctrl_qubit, qubit0, qubit1);
        NotifyObservers({qubit1, qubit0, ctrl_qubit});
    }

    void ApplyCU(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta, double phi, double lambda, double gamma) override
    {
        JoinIfNeeded(ctrl_qubit, tgt_qubit);
        GetSimulator(ctrl_qubit)->ApplyCU(ctrl_qubit, tgt_qubit, theta, phi, lambda, gamma);
        NotifyObservers({tgt_qubit, ctrl_qubit});
    }
 
    void ApplyNop() override
    {
        GetSimulator(0)->ApplyNop();
    }

    void SetMultithreading(bool multithreading = true) override
    {
        enableMultithreading = multithreading;
        for (auto &[id, simulator] : simulators)
            simulator->SetMultithreading(multithreading);
    }

    bool GetMultithreading() const override
    {
        return enableMultithreading;
    }

    bool IsQcsim() const override
    {
        return type == Simulators::SimulatorType::kQCSim;
    }

    void SaveState() override
    {
        savedState = CloneStateForExecution(DeriveSeed(samplingSeed, 0x53415645ULL), false);
    }

    void RestoreState() override
    {
        if (savedState)
        {
            const ImmediateCompositeSimulator *savedStatePtr = static_cast<ImmediateCompositeSimulator *>(savedState.get());
 
            type = savedStatePtr->type;           
            qubitsMap = savedStatePtr->qubitsMap; 
 
            nrQubits = savedStatePtr->nrQubits;                         
            nextId = savedStatePtr->nextId;                             
            enableMultithreading = savedStatePtr->enableMultithreading; 
 
            simulators.clear();
            const uint64_t restoreSeed = DeriveSeed(samplingSeed, 0x52455354ULL + nextRestoreStream++);
            for (auto &[id, simulator] : savedStatePtr->simulators)
            {
                auto isim = simulator->CloneForExecution(DeriveSeed(restoreSeed, id));
                simulators[id] = std::unique_ptr<IndividualSimulator>(static_cast<IndividualSimulator *>(isim.release()));
            }
        }
    }

    Types::qubit_t MeasureNoCollapse() override
    {
        Types::qubit_t res = 0;
 
        if (GetNumberOfQubits() > sizeof(Types::qubit_t) * 8)
            std::cerr << "Warning: The number of qubits to measure is larger than the "
                         "number of bits in the Types::qubit_t type, the outcome will be "
                         "undefined"
                      << std::endl;
 
        for (auto &[id, simulator] : simulators)
        {
            const Types::qubit_t meas = simulator->MeasureNoCollapse();
            res |= meas;
        }
 
        return res;
    }

    std::vector<bool> MeasureNoCollapseMany() override
    {
        std::vector<bool> res(nrQubits, false);
        for (auto &[id, simulator] : simulators)
        {
            const std::vector<bool> meas = simulator->MeasureNoCollapseMany();
            for (size_t i = 0; i < meas.size(); ++i)
                if (meas[i])
                    res[i] = true;
        }
        return res;
    }

    std::unique_ptr<ISimulator> Clone() override
    {
        if (nextCloneStream == std::numeric_limits<uint64_t>::max())
            throw std::overflow_error("Composite clone stream exhausted");
        return CloneForExecution(DeriveSeed(samplingSeed, nextCloneStream++));
    }
 
    std::unique_ptr<ISimulator> CloneForExecution(uint64_t seed) override
    {
        return CloneStateForExecution(seed, true);
    }
 
  private:
    std::unique_ptr<ISimulator> CloneStateForExecution(uint64_t seed, bool includeSnapshot)
    {
        auto clone = std::make_unique<ImmediateCompositeSimulator>(type);
 
        clone->type = type;           
        clone->qubitsMap = qubitsMap; 
 
        clone->nrQubits = nrQubits;                         
        clone->nextId = nextId;                             
        clone->enableMultithreading = enableMultithreading; 
 
        clone->config = config; 
 
        for (auto &[id, simulator] : simulators)
        {
            auto isim = simulator->CloneForExecution(DeriveSeed(seed, id));
            clone->simulators[id] = std::unique_ptr<IndividualSimulator>(static_cast<IndividualSimulator *>(isim.release()));
        }
 
        if (includeSnapshot && savedState)
            clone->savedState = static_cast<ImmediateCompositeSimulator *>(savedState.get())->CloneStateForExecution(DeriveSeed(seed, 0x53415645ULL), false);
 
        clone->SetSeed(seed);
        return clone;
    }
 
  public:
    const Configuration &GetConfiguration() const
    {
        return config;
    }
 
    const std::unordered_map<std::string, std::string> &GetConfigMap() const override
    {
        return config.GetConfigMap();
    }
 
  private:
    template <bool Many> Utils::Sampling::Counts<Many> SamplePrepared(const Types::qubits_vector &qubits, size_t shots)
    {
        using namespace Utils::Sampling;
        using Plan = QCSimStatevectorSimulator::SamplingPlan;
        Counts<Many> result;
        if (qubits.empty() || !shots)
            return result;
        ValidateQubits(qubits, nrQubits, Many);
        if (simulators.size() == 1)
        {
            auto &component = *simulators.begin()->second;
            if constexpr (Many)
                result = component.SampleCountsMany(qubits, shots);
            else
                result = component.SampleCounts(qubits, shots);
            NotifyObservers(qubits);
            return result;
        }
 
        // Duplicate selections share a joint bit; output positions are expanded
        // only when inserting counts. No global qubit ID is used as a shift.
        std::unordered_map<size_t, size_t> uniqueBits;
        std::vector<size_t> outputBits;
        for (const auto qubit : qubits)
            outputBits.push_back(uniqueBits.emplace(qubit, uniqueBits.size()).first->second);
        constexpr size_t wordBits = std::numeric_limits<size_t>::digits;
        const size_t words = (uniqueBits.size() + wordBits - 1) / wordBits;
 
        struct Component
        {
            size_t id;
            IndividualSimulator *simulator;
            std::vector<std::pair<size_t, size_t>> bits; // local bit -> joint bit
            std::unique_ptr<Plan> plan;
            std::vector<size_t> projection;
            size_t minimumBytes = 0;
        };
 
        std::vector<Component> components;
        components.reserve(simulators.size());
        for (auto &[id, simulator] : simulators)
        {
            Component component{id, simulator.get()};
            for (auto [qubit, local] : simulator->GetQubitsMap())
                if (const auto found = uniqueBits.find(qubit); found != uniqueBits.end())
                    component.bits.emplace_back(local, found->second);
            components.push_back(std::move(component));
        }
        // Stable budget assignment; neither hash iteration order nor worker
        // count is allowed to choose a component's sampling algorithm.
        std::sort(components.begin(), components.end(), [](const auto &a, const auto &b) { return a.id < b.id; });
        auto project = [&](const Component &component, size_t row, size_t *joint) {
            if (!component.projection.empty())
                joint[0] |= component.projection[row];
            else
                for (const auto [local, output] : component.bits)
                    joint[output / wordBits] |= ((row >> local) & size_t{1}) << (output % wordBits);
        };
        std::vector<bool> key;
        if constexpr (Many)
            key.resize(qubits.size());
        auto insert = [&](const size_t *joint, size_t count) {
            if constexpr (Many)
            {
                for (size_t bit = 0; bit < outputBits.size(); ++bit)
                    key[bit] = ((joint[outputBits[bit] / wordBits] >> (outputBits[bit] % wordBits)) & 1) != 0;
                result[key] += count;
            }
            else
            {
                Types::qubit_t packed = 0;
                for (size_t bit = 0; bit < outputBits.size(); ++bit)
                    packed |= Types::qubit_t((joint[outputBits[bit] / wordBits] >> (outputBits[bit] % wordBits)) & 1) << bit;
                result[packed] += count;
            }
        };
        if (shots == 1)
        {
            std::vector<size_t> joint(words, 0);
            for (auto &component : components)
                if (component.bits.empty())
                    component.simulator->DiscardLocalMeasurementDraw();
                else
                    project(component, component.simulator->SampleLocalOnce(), joint.data());
            insert(joint.data(), 1);
            NotifyObservers(qubits);
            return result;
        }
 
        const Options options = components.front().simulator->SamplingOptions();
        const size_t sharedReserve = options.memoryBytes / 4;
        size_t remaining = options.memoryBytes - sharedReserve, reserved = 0;
        const size_t selected = std::count_if(components.begin(), components.end(), [](const auto &c) { return !c.bits.empty(); });
        const size_t componentBudget = remaining / std::max<size_t>(1, selected);
        for (auto &component : components)
            if (!component.bits.empty())
            {
                component.minimumBytes = Plan::MinimumMemoryBytes(component.simulator->SamplingRows(), shots, componentBudget);
                if (component.minimumBytes > remaining - reserved)
                    throw std::length_error("Composite probability indexes exceed the sampling scratch budget");
                reserved += component.minimumBytes;
            }
        for (auto &component : components)
            if (!component.bits.empty())
            {
                reserved -= component.minimumBytes;
                auto localOptions = options;
                localOptions.memoryBytes = remaining - reserved;
                component.plan = std::make_unique<Plan>(component.simulator->PrepareSampling(shots, localOptions));
                remaining -= localOptions.memoryBytes - component.plan->ScratchAvailable();
            }
        // Small local projection tables avoid per-shot bit permutation work.
        // Their memory comes from the same shared budget as sampling tables.
        for (auto &component : components)
            if (words == 1 && component.bits.size() > 1 && component.plan && component.plan->Size() <= 65536 &&
                shots >= std::max<size_t>(256, component.plan->Size()) && component.plan->Size() <= remaining / sizeof(size_t))
            {
                std::vector<size_t> table(component.plan->Size(), 0);
                for (size_t row = 0; row < table.size(); ++row)
                    project(component, row, &table[row]);
                remaining -= table.size() * sizeof(size_t);
                component.projection = std::move(table);
            }
        remaining += sharedReserve;
        const size_t categories = uniqueBits.size() <= 20 ? size_t{1} << uniqueBits.size() : 0;
        const bool histogram = categories && shots >= std::max<size_t>(1024, categories / 4) && categories <= remaining / (2 * sizeof(size_t));
        std::vector<size_t> bins(histogram ? categories : 0, 0);
        remaining -= bins.size() * sizeof(size_t);
        // Uniforms, sampled rows, sorted-CDF scratch, and joint output words.
        if (words > (std::numeric_limits<size_t>::max() - 32) / sizeof(size_t))
            throw std::length_error("Composite sampling output is too wide");
        const size_t bytesPerShot = 32 + words * sizeof(size_t);
        const size_t batch = std::min({shots, options.batchSize, remaining / bytesPerShot});
        if (!batch)
            throw std::length_error("Composite sampling batch exceeds the scratch budget");
        std::vector<size_t> joint(batch * words), sampled;
        std::vector<double> uniforms;
        // Unselected components still consume their original one word per shot.
        // Omitting their probability preparation cannot restart later streams.
        for (auto &component : components)
            if (component.bits.empty())
                component.simulator->SamplingRng().discard(shots);
        for (size_t begin = 0; begin < shots;)
        {
            const size_t count = std::min(batch, shots - begin);
            std::fill(joint.begin(), joint.begin() + count * words, 0);
            for (auto &component : components)
            {
                if (!component.plan)
                    continue;
                const auto &plan = *component.plan;
                auto &rng = component.simulator->SamplingRng();
                const int threads = plan.Categories() <= 64 && plan.Method() == Family::Alias ? 1 : Threads(options.multithreading, count);
                if (threads == 1 && (plan.Method() == Family::Alias || plan.Size() <= 65536 || count < 128))
                {
                    for (size_t shot = 0; shot < count; ++shot)
                        project(component, plan.Sample(Utils::RandomStream::Uniform(rng)), &joint[shot * words]);
                }
                else
                {
                    uniforms.resize(count);
                    for (double &value : uniforms)
                        value = Utils::RandomStream::Uniform(rng);
                    plan.Draw(uniforms, sampled, threads);
#ifdef _OPENMP
#pragma omp parallel for if (threads > 1) num_threads(threads) schedule(static)
#endif
                    for (int64_t shot = 0; shot < static_cast<int64_t>(count); ++shot)
                        project(component, sampled[shot], &joint[size_t(shot) * words]);
                }
            }
            for (size_t shot = 0; shot < count; ++shot)
                if (histogram)
                    ++bins[joint[shot]];
                else
                    insert(&joint[shot * words], 1);
            begin += count;
        }
        for (size_t row = 0; row < bins.size(); ++row)
            if (bins[row])
                insert(&row, bins[row]);
        NotifyObservers(qubits);
        return result;
    }
 
    void InitializeAlias()
    {
        for (auto &[id, simulator] : simulators)
            simulator->InitializeAlias();
    }
 
    void ClearAlias()
    {
        for (auto &[id, simulator] : simulators)
            simulator->ClearAlias();
    }

    inline std::unique_ptr<IndividualSimulator> &GetSimulator(size_t qubit)
    {
        assert(qubitsMap.size() > qubit);
        assert(simulators.find(qubitsMap[qubit]) != simulators.end());
 
        // assume it was called with a valid qubit
        return simulators[qubitsMap[qubit]];
    }
 
    // if executing a gate on two or three qubits, use this to see if it can be
    // executed directly if it returns true, join the two simulators together and
    // execute the gate (for three qubits, do that once again for the third qubit,
    // before execution)
    inline bool JoinNeeded(size_t qubit1, size_t qubit2)
    {
        // TODO: Not really needed for all gates that act on multiple qubits, is
        // this worth pursuing? the obvious example is the identity gate, which does
        // nothing, so a join is not needed
        return qubitsMap[qubit1] != qubitsMap[qubit2];
    }

    inline void JoinIfNeeded(size_t qubit1, size_t qubit2)
    {
        if (JoinNeeded(qubit1, qubit2))
        {
            const size_t simId = qubitsMap[qubit1];
            const size_t eraseSimId = qubitsMap[qubit2];
            auto &sim1 = GetSimulator(qubit1);
            const auto &sim2 = GetSimulator(qubit2);
 
            sim1->Join(simId, sim2, qubitsMap, enableMultithreading);
 
            simulators.erase(eraseSimId);
        }
    }

    inline void Split(size_t qubit, bool qubitOutcome = false)
    {
        auto &sim = GetSimulator(qubit);   // get the simulator for the qubit, this one will be split
        if (sim->GetNumberOfQubits() == 1) // no need to split it, it's already for a single qubit
            return;
 
        if (nextSplitStream == std::numeric_limits<uint64_t>::max())
            throw std::overflow_error("Composite split stream exhausted");
        // A detached component gets its own stream. Restores and register
        // rebuilds must not rewind this ordinal and reuse a child's seed.
        const uint64_t splitSeed = DeriveSeed(DeriveSeed(samplingSeed, 0x53504C4954ULL), nextSplitStream++);
        qubitsMap[qubit] = nextId; // the qubit will be in the new simulator
        simulators[nextId] = sim->Split(qubit, qubitOutcome, enableMultithreading, splitSeed);
 
        ++nextId;
    }
 
    SimulatorType type;                                                          
    std::vector<size_t> qubitsMap;                                               
    std::unordered_map<size_t, std::unique_ptr<IndividualSimulator>> simulators; 
    size_t nrQubits = 0;                                                         
    size_t nextId = 0;                                                           
    bool enableMultithreading = true;                                            
 
    std::unique_ptr<ISimulator> savedState; 
 
    Configuration config; 
    uint64_t samplingSeed = Utils::RandomStream::FreshSeed();
    uint64_t nextCloneStream = 0, nextRestoreStream = 0, nextInitialization = 0, nextSplitStream = 0;
};
 
} // namespace Private
} // namespace Simulators
 
#endif
#endif