Individual.h

Go to the documentation of this file.

 
#ifndef _INDIVIDUAL_H
#define _INDIVIDUAL_H
 
#ifdef INCLUDED_BY_FACTORY
 
#include "../../Utils/Sampling/Alias.h"
#include "../Aer/AerSimulator.h"
#include "../Core/Factory.h"
#include "../QCSim/QCSimStatevectorSimulator.h"
#include <unordered_map>
 
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;

class IndividualSimulator : public ISimulator
{
    friend class ImmediateCompositeSimulator;
 
  public:
    IndividualSimulator(SimulatorType type =
#ifdef NO_QISKIT_AER
                            Simulators::SimulatorType::kQCSim
#else
                            Simulators::SimulatorType::kQiskitAer
#endif
                        ) noexcept
        : simulator(SimulatorsFactory::CreateImmediateSimulatorUnique(type))
    {
    }

    void Reset() override
    {
        simulator->Reset();
    }

    inline void Join(size_t simId, const std::unique_ptr<IndividualSimulator> &other, std::vector<size_t> &qubitsMapToSim, bool enableMultithreading)
    {
        // 1. grab the state of both simulators (in the first phase, using
        // Amplitude, but maybe something faster can be done)
        // 2. join the states by a tensor product
        const size_t nrQubits1 = GetNumberOfQubits();
        const size_t nrBasisStates1 = 1ULL << nrQubits1;
        const size_t nrQubits2 = other->GetNumberOfQubits();
        const size_t nrBasisStates2 = 1ULL << nrQubits2;
 
        const size_t newNrQubits = nrQubits1 + nrQubits2;
        const size_t nrBasisStates = 1ULL << newNrQubits;
 
        SaveStateToInternalDestructive();
        other->SaveStateToInternalDestructive();
 
        if (GetType() == SimulatorType::kQCSim)
        {
            if (enableMultithreading && nrBasisStates > OmpLimitJoin)
                JoinOmpQcsim(nrQubits1, nrBasisStates1, nrBasisStates2, newNrQubits, nrBasisStates, other, enableMultithreading);
            else
            {
                Eigen::VectorXcd newAmplitudes;
                newAmplitudes.resize(nrBasisStates);
 
                for (size_t state2 = 0; state2 < nrBasisStates2; ++state2)
                {
                    const auto ampl2 = other->AmplitudeRaw(state2);
                    const size_t state2Mask = state2 << nrQubits1;
                    for (size_t state1 = 0; state1 < nrBasisStates1; ++state1)
                        newAmplitudes[state2Mask | state1] = AmplitudeRaw(state1) * ampl2;
                }
 
                // 3. set the state of the current simulator to the joined state
                // the original qubits of this simulator get mapped as they are
                // the other ones get shifted to the left by the number of qubits of
                // this simulator so transfer mapping keeping this in mind
 
                // simulator->SetMultithreading(enableMultithreading);
                static_cast<QCSimStatevectorSimulator *>(simulator.get())
                    ->ReplaceState(newNrQubits,
                                   newAmplitudes); // this will end up by swapping the data from
                                                   // newAmplitudes to the simulator, no allocation
                                                   // and copying is done
            }
        }
#ifndef NO_QISKIT_AER
        else
        {
            if (enableMultithreading && nrBasisStates > OmpLimitJoin)
                JoinOmpAer(nrQubits1, nrBasisStates1, nrBasisStates2, newNrQubits, nrBasisStates, other, enableMultithreading);
            else
            {
                AER::Vector<std::complex<double>> newAmplitudes(nrBasisStates, false); // the false here avoids data
                                                                                       // initialization, it will be set anyway
 
                for (size_t state2 = 0; state2 < nrBasisStates2; ++state2)
                {
                    const auto ampl2 = other->AmplitudeRaw(state2);
                    const size_t state2Mask = state2 << nrQubits1;
                    for (size_t state1 = 0; state1 < nrBasisStates1; ++state1)
                        newAmplitudes[state2Mask | state1] = AmplitudeRaw(state1) * ampl2;
                }
 
                // 3. set the state of the current simulator to the joined state
                // the original qubits of this simulator get mapped as they are
                // the other ones get shifted to the left by the number of qubits of
                // this simulator so transfer mapping keeping this in mind
                // simulator->SetMultithreading(enableMultithreading);
                simulator->InitializeState(newNrQubits,
                                           newAmplitudes); // this will move the data from newAmplitudes to
                                                           // the simulator, no allocation and copying is done
            }
        }
#endif
 
        for (auto [origq, mapq] : other->GetQubitsMap())
        {
            qubitsMap[origq] = mapq + nrQubits1;
            qubitsMapToSim[origq] = simId;
        }
    }

    inline std::unique_ptr<IndividualSimulator> Split(size_t qubit, bool qubitOutcome, bool enableMultithreading, uint64_t splitSeed)
    {
        const size_t oldNrQubits = GetNumberOfQubits();
        const size_t newNrQubits = oldNrQubits - 1;
        const size_t nrBasisStates = 1ULL << newNrQubits;
        const size_t localQubit = qubitsMap[qubit];
 
        // the new simulator split from this one
        // a one qubit one, containing the qubit that was measured or reset
        // TODO: this can be optimized a little bit by directly initializing with
        // the proper amplitudes, but I'm not sure if it's worth it
        auto newSimulator = std::make_unique<IndividualSimulator>(GetType());
        newSimulator->AllocateQubits(1);
        newSimulator->GetQubitsMap()[qubit] = 0; // the qubit is mapped to the only local qubit in the new simulator,
                                                 // which is 0
        newSimulator->SetMultithreading(enableMultithreading);
 
        for (const auto &[key, value] : GetConfigMap())
            if (key != "seed")
                newSimulator->Configure(key.c_str(), value.c_str());
        newSimulator->SetSeed(splitSeed);
 
        newSimulator->Initialize();
        if (qubitOutcome)
        {
            newSimulator->ApplyX(qubit);
            // newSimulator->Flush();
        }
 
        qubitsMap.erase(qubit); // the qubit is removed from the current simulator
 
        SaveStateToInternalDestructive();
 
        if (GetType() == SimulatorType::kQCSim)
        {
            /*
            if (nrBasisStates > OmpLimitSplit) // parallelization for assignment and
            some bit manipulations, I must do some benchmarks to see if it's worth it
            and find where the limit is SplitOmpQcsim(localQubit, newNrQubits,
            nrBasisStates, qubitOutcome); else
            */
            {
                // now the adjusted current simulator, without the removed qubit
                Eigen::VectorXcd newAmplitudes;
                newAmplitudes.resize(nrBasisStates);
 
                // compute the new amplitudes
 
                const size_t localQubitMask = 1ULL << localQubit;
                const size_t maskLow = localQubitMask - 1ULL;
                const size_t maskHigh = ~maskLow;
                const size_t qubitMask = qubitOutcome ? localQubitMask : 0ULL;
 
                for (size_t state = 0; state < nrBasisStates; ++state)
                {
                    const size_t stateLow = state & maskLow;
                    const size_t stateHigh = (state & maskHigh) << 1ULL;
 
                    newAmplitudes[state] = AmplitudeRaw(stateLow | stateHigh | qubitMask);
                }
 
                // simulator->SetMultithreading(enableMultithreading);
                static_cast<QCSimStatevectorSimulator *>(simulator.get())
                    ->ReplaceState(newNrQubits,
                                   newAmplitudes); // this will end up by swapping the data from
                                                   // newAmplitudes to the simulator, no allocation
                                                   // and copying is done
            }
        }
#ifndef NO_QISKIT_AER
        else
        {
            /*
            if (nrBasisStates > OmpLimitSplit) // parallelization for assignment and
            some bit manipulations, I must do some benchmarks to see if it's worth it
            and find where the limit is SplitOmpAer(localQubit, newNrQubits,
            nrBasisStates, qubitOutcome); else
            */
            {
                // now the adjusted current simulator, without the removed qubit
                AER::Vector<std::complex<double>> newAmplitudes(nrBasisStates, false); // the false here avoids data
                                                                                       // initialization, it will be set anyway
 
                // compute the new amplitudes
                const size_t localQubitMask = 1ULL << localQubit;
                const size_t maskLow = localQubitMask - 1ULL;
                const size_t maskHigh = ~maskLow;
                const size_t qubitMask = qubitOutcome ? localQubitMask : 0ULL;
 
                for (size_t state = 0; state < nrBasisStates; ++state)
                {
                    const size_t stateLow = state & maskLow;
                    const size_t stateHigh = (state & maskHigh) << 1ULL;
 
                    newAmplitudes[state] = AmplitudeRaw(stateLow | stateHigh | qubitMask);
                }
 
                // simulator->SetMultithreading(enableMultithreading);
                simulator->InitializeState(newNrQubits,
                                           newAmplitudes); // this will move the data from newAmplitudes to
                                                           // the simulator, no allocation and copying is done
            }
        }
#endif
 
        // now adjust the local qubits map
        for (auto &mapped : qubitsMap)
            if (mapped.second > localQubit)
                --mapped.second;
 
        return newSimulator;
    }

    std::complex<double> AmplitudeRaw(Types::qubit_t outcome) override
    {
        return simulator->AmplitudeRaw(outcome);
    }

    void SaveStateToInternalDestructive() override
    {
        simulator->SaveStateToInternalDestructive();
    }

    void RestoreInternalDestructiveSavedState() override
    {
        simulator->RestoreInternalDestructiveSavedState();
    }

    inline Types::qubits_vector ConvertQubits(const Types::qubits_vector &qubits)
    {
        Types::qubits_vector converted;
        converted.reserve(qubits.size());
 
        for (auto qubit : qubits)
            if (HasQubit(qubit))
                converted.emplace_back(qubitsMap[qubit]);
 
        return converted;
    }

    inline bool HasQubit(Types::qubit_t qubit) const
    {
        return qubitsMap.find(qubit) != qubitsMap.end();
    }

    inline Types::qubit_t ConvertOutcomeFromLocal(Types::qubit_t outcome) const
    {
        Types::qubit_t res = 0;
 
        for (auto [origQubit, localQubit] : qubitsMap)
            if (outcome & (1ULL << localQubit))
                res |= (1ULL << origQubit);
 
        return res;
    }

    inline Types::qubit_t ConvertOutcomeFromGlobal(Types::qubit_t outcome) const
    {
        Types::qubit_t res = 0;
 
        for (auto [origQubit, localQubit] : qubitsMap)
            if (outcome & (1ULL << origQubit))
                res |= (1ULL << localQubit);
 
        return res;
    }

    inline std::vector<bool> ConvertOutcomeFromLocal(const std::vector<bool> &outcome) const
    {
        std::vector<bool> res;
 
        size_t maxQubit = 0;
        for (auto [origQubit, localQubit] : qubitsMap)
            if (origQubit > maxQubit)
                maxQubit = origQubit;
        res.resize(maxQubit + 1, false);
 
        for (auto [origQubit, localQubit] : qubitsMap)
            if (outcome[localQubit])
                res[origQubit] = true;
 
        return res;
    }

    inline std::vector<bool> ConvertOutcomeFromGlobal(const std::vector<bool> &outcome) const
    {
        const size_t nrQubits = simulator->GetNumberOfQubits();
        std::vector<bool> res(nrQubits, false);
 
        for (auto [origQubit, localQubit] : qubitsMap)
            if (outcome[origQubit])
                res[localQubit] = true;
 
        return res;
    }

    void SaveState() override
    {
        if (!simulator)
            return;
        const size_t nrBasisStates = 1ULL << simulator->GetNumberOfQubits();
        savedState.resize(nrBasisStates);
 
        for (Types::qubit_t state = 0; state < nrBasisStates; ++state)
            savedState[state] = simulator->Amplitude(state);
    }

    void ClearSavedState()
    {
        savedState.clear();
    }

    void RestoreState() override
    {
        if (!simulator || savedState.empty())
            return;
        const size_t nrQubits = simulator->GetNumberOfQubits();
 
        if (GetType() == SimulatorType::kQCSim)
        {
            Eigen::VectorXcd amplitudes(Eigen::Map<const Eigen::VectorXcd>(savedState.data(), savedState.size()));
            static_cast<QCSimStatevectorSimulator *>(simulator.get())->ReplaceState(nrQubits, amplitudes);
        }
        else
        {
            simulator->Clear();
            simulator->InitializeState(nrQubits, savedState);
        }
        ClearSavedState();
    }

    inline std::unordered_map<Types::qubit_t, Types::qubit_t> &GetQubitsMap()
    {
        return qubitsMap;
    }

    inline const std::unordered_map<Types::qubit_t, Types::qubit_t> &GetQubitsMap() const
    {
        return qubitsMap;
    }

    void Initialize() override
    {
        simulator->Initialize();
    }

    void InitializeState(size_t num_qubits, std::vector<std::complex<double>> &amplitudes) override
    {
        simulator->InitializeState(num_qubits, amplitudes);
    }

    /*
    void InitializeState(size_t num_qubits, std::vector<std::complex<double>,
    avoid_init_allocator<std::complex<double>>>& amplitudes) override
    {
            simulator->InitializeState(num_qubits, amplitudes);
    }
    */

#ifndef NO_QISKIT_AER
    void InitializeState(size_t num_qubits, AER::Vector<std::complex<double>> &amplitudes) override
    {
        simulator->InitializeState(num_qubits, amplitudes);
    }
#endif

    void InitializeState(size_t num_qubits, Eigen::VectorXcd &amplitudes) override
    {
        simulator->InitializeState(num_qubits, amplitudes);
    }

    void Configure(const char *key, const char *value) override
    {
        if (std::string("seed") == key)
        {
            const uint64_t seed = std::stoull(value);
            SeedAuxiliaryRng(seed);
            simulator->SetSeed(seed);
        }
        else
        {
            simulator->Configure(key, value);
        }
    }

    std::string GetConfiguration(const char *key) const override
    {
        if (!simulator)
            return "";
 
        return simulator->GetConfiguration(key);
    }

    size_t AllocateQubits(size_t num_qubits) override
    {
        return simulator->AllocateQubits(num_qubits);
    }

    size_t GetNumberOfQubits() const override
    {
        return simulator->GetNumberOfQubits();
    }

    void Clear() override
    {
        simulator->Clear();
    }

    size_t Measure(const Types::qubits_vector &qubits) override
    {
        return ConvertOutcomeFromLocal(simulator->Measure(ConvertQubits(qubits)));
    }

    std::vector<bool> MeasureMany(const Types::qubits_vector &qubits) override
    {
        return ConvertOutcomeFromLocal(simulator->MeasureMany(ConvertQubits(qubits)));
    }

    void ApplyReset(const Types::qubits_vector &qubits) override
    {
        simulator->ApplyReset(ConvertQubits(qubits));
    }

    double Probability(Types::qubit_t outcome) override
    {
        return simulator->Probability(ConvertOutcomeFromGlobal(outcome));
    }

    std::complex<double> Amplitude(Types::qubit_t outcome) override
    {
        return simulator->Amplitude(ConvertOutcomeFromGlobal(outcome));
    }

    std::complex<double> ProjectOnZero() override
    {
        return Amplitude(0);
    }

    std::vector<double> AllProbabilities() override
    {
        return simulator->AllProbabilities();
    }

    std::vector<double> Probabilities(const Types::qubits_vector &qubits) override
    {
        return simulator->Probabilities(ConvertQubits(qubits));
    }

    std::unordered_map<Types::qubit_t, Types::qubit_t> SampleCounts(const Types::qubits_vector &qubits, size_t shots = 1000) override
    {
        // Results are indexed by positions in the requested list, not by local
        // register IDs. Convert the input IDs only; preserve the returned order.
        return simulator->SampleCounts(ConvertQubits(qubits), shots);
    }

    std::unordered_map<std::vector<bool>, Types::qubit_t> SampleCountsMany(const Types::qubits_vector &qubits, size_t shots = 1000) override
    {
        return simulator->SampleCountsMany(ConvertQubits(qubits), shots);
    }

    double ExpectationValue(const std::string &pauliString) override
    {
        return simulator->ExpectationValue(pauliString);
    }

    SimulatorType GetType() const override
    {
        return simulator->GetType();
    }

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

    void Flush() override
    {
        simulator->Flush();
    }
 
    // WARNING: for all the following functions, the call is supposed to be made
    // after the proper joining (or splitting)!

    void ApplyGenericOneQubitGate(Types::qubit_t qubit, const Eigen::Matrix2cd &gate) override
    {
        simulator->ApplyGenericOneQubitGate(qubitsMap[qubit], gate);
    }

    void ApplyGenericTwoQubitGate(Types::qubit_t qubit0, Types::qubit_t qubit1, const Eigen::Matrix4cd &gate) override
    {
        simulator->ApplyGenericTwoQubitGate(qubitsMap[qubit0], qubitsMap[qubit1], gate);
    }
 
    void ApplyGenericThreeQubitGate(Types::qubit_t q0, Types::qubit_t q1, Types::qubit_t q2, const Matrix8cd &gate) override
    {
        simulator->ApplyGenericThreeQubitGate(qubitsMap.at(q0), qubitsMap.at(q1), qubitsMap.at(q2), gate);
    }

    void ApplyP(Types::qubit_t qubit, double lambda) override
    {
        simulator->ApplyP(qubitsMap[qubit], lambda);
    }

    void ApplyX(Types::qubit_t qubit) override
    {
        simulator->ApplyX(qubitsMap[qubit]);
    }

    void ApplyY(Types::qubit_t qubit) override
    {
        simulator->ApplyY(qubitsMap[qubit]);
    }

    void ApplyZ(Types::qubit_t qubit) override
    {
        simulator->ApplyZ(qubitsMap[qubit]);
    }

    void ApplyH(Types::qubit_t qubit) override
    {
        simulator->ApplyH(qubitsMap[qubit]);
    }

    void ApplyS(Types::qubit_t qubit) override
    {
        simulator->ApplyS(qubitsMap[qubit]);
    }

    void ApplySDG(Types::qubit_t qubit) override
    {
        simulator->ApplySDG(qubitsMap[qubit]);
    }

    void ApplyT(Types::qubit_t qubit) override
    {
        simulator->ApplyT(qubitsMap[qubit]);
    }

    void ApplyTDG(Types::qubit_t qubit) override
    {
        simulator->ApplyTDG(qubitsMap[qubit]);
    }

    void ApplySx(Types::qubit_t qubit) override
    {
        simulator->ApplySx(qubitsMap[qubit]);
    }

    void ApplySxDAG(Types::qubit_t qubit) override
    {
        simulator->ApplySxDAG(qubitsMap[qubit]);
    }

    void ApplyK(Types::qubit_t qubit) override
    {
        simulator->ApplyK(qubitsMap[qubit]);
    }

    void ApplyRx(Types::qubit_t qubit, double theta) override
    {
        simulator->ApplyRx(qubitsMap[qubit], theta);
    }

    void ApplyRy(Types::qubit_t qubit, double theta) override
    {
        simulator->ApplyRy(qubitsMap[qubit], theta);
    }

    void ApplyRz(Types::qubit_t qubit, double theta) override
    {
        simulator->ApplyRz(qubitsMap[qubit], theta);
    }

    void ApplyU(Types::qubit_t qubit, double theta, double phi, double lambda, double gamma) override
    {
        simulator->ApplyU(qubitsMap[qubit], theta, phi, lambda, gamma);
    }

    void ApplyCX(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCX(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplyCY(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCY(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplyCZ(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCZ(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplyCP(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double lambda) override
    {
        simulator->ApplyCP(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit], lambda);
    }

    void ApplyCRx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        simulator->ApplyCRx(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit], theta);
    }

    void ApplyCRy(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        simulator->ApplyCRy(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit], theta);
    }

    void ApplyCRz(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        simulator->ApplyCRz(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit], theta);
    }

    void ApplyCH(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCH(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplyCSx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCSx(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplyCSxDAG(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        simulator->ApplyCSxDAG(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit]);
    }

    void ApplySwap(Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        simulator->ApplySwap(qubitsMap[qubit0], qubitsMap[qubit1]);
    }

    void ApplyCCX(Types::qubit_t qubit0, Types::qubit_t qubit1, Types::qubit_t qubit2) override
    {
        simulator->ApplyCCX(qubitsMap[qubit0], qubitsMap[qubit1], qubitsMap[qubit2]);
    }

    void ApplyCSwap(Types::qubit_t ctrl_qubit, Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        simulator->ApplyCSwap(qubitsMap[ctrl_qubit], qubitsMap[qubit0], qubitsMap[qubit1]);
    }

    void ApplyCU(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta, double phi, double lambda, double gamma) override
    {
        simulator->ApplyCU(qubitsMap[ctrl_qubit], qubitsMap[tgt_qubit], theta, phi, lambda, gamma);
    }
 
    void ApplyNop() override
    {
        simulator->ApplyNop();
    }

    void SetMultithreading(bool multithreading = true) override
    {
        if (simulator)
            simulator->SetMultithreading(multithreading);
 
        processor_count = multithreading ? QC::QubitRegister<>::GetNumberOfThreads() : 1;
    }

    bool GetMultithreading() const override
    {
        if (simulator)
            return simulator->GetMultithreading();
 
        return false;
    }

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

    Types::qubit_t MeasureNoCollapse() override
    {
        return ConvertOutcomeFromLocal(simulator->MeasureNoCollapse());
    }

    std::vector<bool> MeasureNoCollapseMany() override
    {
        auto res = simulator->MeasureNoCollapseMany();
        return ConvertOutcomeFromLocal(res);
    }

    std::unique_ptr<ISimulator> Clone() override
    {
        auto cloned = std::make_unique<IndividualSimulator>();
 
        cloned->qubitsMap = qubitsMap;          
        cloned->savedState = savedState;        
        cloned->simulator = simulator->Clone(); 
 
        return cloned;
    }
 
    std::unique_ptr<ISimulator> CloneForExecution(uint64_t seed) override
    {
        auto cloned = std::make_unique<IndividualSimulator>();
        cloned->qubitsMap = qubitsMap;
        cloned->savedState = savedState;
        cloned->simulator = simulator->CloneForExecution(seed);
        return cloned;
    }
 
    Types::qubit_t SampleFromAlias()
    {
        if (!alias || !simulator)
            return 0;
 
        double prob = 0.0;
        if (GetType() == SimulatorType::kQCSim)
        {
            // qcsim - convert 'simulator' to qcsim simulator and access 'state' (from
            // there the statevector is accessible)
            QCSimStatevectorSimulator *qcsim = dynamic_cast<QCSimStatevectorSimulator *>(simulator.get());
            prob = 1. - qcsim->uniformZeroOne(qcsim->rng);
        }
#ifndef NO_QISKIT_AER
        else
        {
            // qiskit aer - convert 'simulator' to qiskit aer simulator and access
            // 'savedAmplitudes' (assumes destructive saving of the state)
#ifndef NO_QISKIT_AER
            AerSimulator *aer = dynamic_cast<AerSimulator *>(simulator.get());
            prob = 1 - aer->uniformZeroOne(aer->rng);
#else
            return 0;
#endif
        }
#endif
 
        const size_t measRaw = alias->Sample(prob);
 
        return ConvertOutcomeFromLocal(measRaw);
    }
 
    const std::unordered_map<std::string, std::string> &GetConfigMap() const override
    {
        if (!simulator)
        {
            static const std::unordered_map<std::string, std::string> emptyMap;
            return emptyMap;
        }
 
        return simulator->GetConfigMap();
    }
 
  private:
    // These helpers operate on component-local rows. The composite projects
    // directly into requested output positions, including global IDs above 63.
    QCSimStatevectorSimulator &SamplingStatevector() const
    {
        return *static_cast<QCSimStatevectorSimulator *>(simulator.get());
    }
 
    size_t SamplingRows() const
    {
        return size_t{1} << SamplingStatevector().samplingSupportQubits;
    }
 
    Utils::Sampling::Options SamplingOptions() const
    {
        return SamplingStatevector().SamplingOptions();
    }
 
    std::mt19937_64 &SamplingRng()
    {
        return SamplingStatevector().rng;
    }
 
    QCSimStatevectorSimulator::SamplingPlan PrepareSampling(size_t shots, Utils::Sampling::Options options) const
    {
        return {SamplingRows(), shots, {&SamplingStatevector().state->getRegisterStorage()}, options};
    }
 
    size_t SampleLocalOnce()
    {
        auto &backend = *SamplingStatevector().state;
 
        struct RestoreThreading
        {
            QC::QubitRegister<> &backend;
            bool enabled;
 
            ~RestoreThreading()
            {
                backend.SetMultithreading(enabled);
            }
        } restore{backend, simulator->GetMultithreading()};
 
        backend.SetMultithreading(false);
        return backend.MeasureNoCollapse();
    }
 
    void DiscardLocalMeasurementDraw()
    {
        SamplingStatevector().state->DiscardMeasurementDraw();
    }
 
    void InitializeAlias()
    {
        if (GetType() == SimulatorType::kQCSim)
        {
            // qcsim - convert 'simulator' to qcsim simulator and access 'state' (from
            // there the statevector is accessible)
            QCSimStatevectorSimulator *qcsim = dynamic_cast<QCSimStatevectorSimulator *>(simulator.get());
            alias = std::make_unique<Utils::Alias>(qcsim->state->getRegisterStorage());
        }
#ifndef NO_QISKIT_AER
        else
        {
            // qiskit aer - convert 'simulator' to qiskit aer simulator and access
            // 'savedAmplitudes' (assumes destructive saving of the state)
#ifndef NO_QISKIT_AER
            AerSimulator *aer = dynamic_cast<AerSimulator *>(simulator.get());
            if (aer)
            {
                alias = std::unique_ptr<Utils::Alias>(new Utils::Alias(aer->savedAmplitudes));
            }
#else
            // If we are here, it means it's not QCSim, but Qiskit is disabled.
            throw std::runtime_error("Qiskit Aer is disabled in this build.");
#endif
        }
#endif
    }
 
    void ClearAlias()
    {
        alias = nullptr;
    }

#ifndef NO_QISKIT_AER
    inline void JoinOmpAer(size_t nrQubits1, size_t nrBasisStates1, size_t nrBasisStates2, size_t newNrQubits, size_t nrBasisStates,
                           const std::unique_ptr<IndividualSimulator> &other, bool enableMultithreading)
    {
        AER::Vector<std::complex<double>> newAmplitudes(nrBasisStates, false); // the false here avoids data initialization, it
                                                                               // will be set anyway
 
        /*
        const size_t state1Mask = nrBasisStates1 - 1ULL;
 
    #pragma omp parallel for num_threads(processor_count) schedule(static,
    OmpLimitJoin / divSchedule) for (long long int state = 0; state <
    static_cast<long long int>(nrBasisStates); ++state) newAmplitudes[state] =
    AmplitudeRaw(state & state1Mask) * other->AmplitudeRaw(state >> nrQubits1);
        */
 
        // TODO: check if this is better
#pragma omp parallel for num_threads(processor_count)
        for (long long int state2 = 0; state2 < static_cast<long long int>(nrBasisStates2); ++state2)
        {
            const auto ampl2 = other->AmplitudeRaw(state2);
            const size_t state2Mask = state2 << nrQubits1;
            for (size_t state1 = 0; state1 < nrBasisStates1; ++state1)
                newAmplitudes[state2Mask | state1] = AmplitudeRaw(state1) * ampl2;
        }
 
        // 3. set the state of the current simulator to the joined state
        // the original qubits of this simulator get mapped as they are
        // the other ones get shifted to the left by the number of qubits of this
        // simulator so transfer mapping keeping this in mind
        // simulator->SetMultithreading(enableMultithreading);
        simulator->InitializeState(newNrQubits,
                                   newAmplitudes); // this will move the data from newAmplitudes to the
                                                   // simulator, no allocation and copying is done
    }
#endif

    /*
    inline void SplitOmpAer(size_t localQubit, size_t newNrQubits, size_t
  nrBasisStates, bool qubitOutcome = false)
    {
            // now the adjusted current simulator, without the removed qubit
            AER::Vector<std::complex<double>> newAmplitudes(nrBasisStates, false);
  // the false here avoids data initialization, it will be set anyway
 
            // compute the new amplitudes
 
            const size_t localQubitMask = 1ULL << localQubit;
            const size_t maskLow = localQubitMask - 1ULL;
            const size_t maskHigh = ~maskLow;
            const size_t qubitMask = qubitOutcome ? localQubitMask : 0ULL;
 
  #pragma omp parallel for num_threads(processor_count) schedule(static,
  OmpLimitSplit / divSchedule) for (long long int state = 0; state <
  static_cast<long long int>(nrBasisStates); ++state)
            {
                    const size_t stateLow = state & maskLow;
                    const size_t stateHigh = (state & maskHigh) << 1ULL;
 
                    newAmplitudes[state] = AmplitudeRaw(stateLow | stateHigh |
  qubitMask);
            }
 
            simulator->InitializeState(newNrQubits, newAmplitudes); // this will
  move the data from newAmplitudes to the simulator, no allocation and copying is
  done
    }
    */

    inline void JoinOmpQcsim(size_t nrQubits1, size_t nrBasisStates1, size_t nrBasisStates2, size_t newNrQubits, size_t nrBasisStates,
                             const std::unique_ptr<IndividualSimulator> &other, bool enableMultithreading)
    {
        Eigen::VectorXcd newAmplitudes;
        newAmplitudes.resize(nrBasisStates);
 
        /*
        const size_t state1Mask = nrBasisStates1 - 1ULL;
 
    #pragma omp parallel for num_threads(processor_count) schedule(static,
    OmpLimitJoin / divSchedule) for (long long int state = 0; state <
    static_cast<long long int>(nrBasisStates); ++state) newAmplitudes[state] =
    AmplitudeRaw(state & state1Mask) * other->AmplitudeRaw(state >> nrQubits1);
        */
 
        // TODO: check if this is better
#pragma omp parallel for num_threads(processor_count)
        for (long long int state2 = 0; state2 < static_cast<long long int>(nrBasisStates2); ++state2)
        {
            const auto ampl2 = other->AmplitudeRaw(state2);
            const size_t state2Mask = state2 << nrQubits1;
            for (size_t state1 = 0; state1 < nrBasisStates1; ++state1)
                newAmplitudes[state2Mask | state1] = AmplitudeRaw(state1) * ampl2;
        }
 
        // 3. set the state of the current simulator to the joined state
        // the original qubits of this simulator get mapped as they are
        // the other ones get shifted to the left by the number of qubits of this
        // simulator so transfer mapping keeping this in mind
        // simulator->SetMultithreading(enableMultithreading);
        static_cast<QCSimStatevectorSimulator *>(simulator.get())->ReplaceState(newNrQubits, newAmplitudes); // swaps the data
                                                                                                             // from newAmplitudes to the simulator, no
                                                                                                             // allocation and copying is done
    }

    /*
    inline void SplitOmpQcsim(size_t localQubit, size_t newNrQubits, size_t
  nrBasisStates, bool qubitOutcome = false)
    {
            // now the adjusted current simulator, without the removed qubit
            Eigen::VectorXcd newAmplitudes;
            newAmplitudes.resize(nrBasisStates);
 
            // compute the new amplitudes
            const size_t localQubitMask = 1ULL << localQubit;
            const size_t maskLow = localQubitMask - 1ULL;
            const size_t maskHigh = ~maskLow;
            const size_t qubitMask = qubitOutcome ? localQubitMask : 0ULL;
 
  #pragma omp parallel for num_threads(processor_count) schedule(static,
  OmpLimitSplit / divSchedule) for (long long int state = 0; state <
  static_cast<long long int>(nrBasisStates); ++state)
            {
                    const size_t stateLow = state & maskLow;
                    const size_t stateHigh = (state & maskHigh) << 1ULL;
 
                    newAmplitudes[state] = AmplitudeRaw(stateLow | stateHigh |
  qubitMask);
            }
 
            static_cast<QCSimStatevectorSimulator *>(simulator.get())->ReplaceState(newNrQubits, newAmplitudes); // this will
  end up by swapping the data from newAmplitudes to the simulator, no allocation
  and copying is done
    }
    */
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsMap; 
    std::unique_ptr<ISimulator> simulator;                        
    std::vector<std::complex<double>> savedState;                 
 
    std::unique_ptr<Utils::Alias> alias; 
 
    int processor_count = QC::QubitRegister<>::GetNumberOfThreads(); 
 
    // constexpr static int divSchedule = 4;
    constexpr static size_t OmpLimitJoin = 4096 * 2;
    // constexpr static size_t OmpLimitSplit = OmpLimitJoin * 16;
};
 
} // namespace Private
} // namespace Simulators
 
#endif
#endif