Individual.h

Go to the documentation of this file.

 
#ifndef _INDIVIDUAL_H
#define _INDIVIDUAL_H
 
#ifdef INCLUDED_BY_FACTORY
 
#include "../Utils/Alias.h"
#include "Factory.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 CompositeSimulator;
 
class IndividualSimulator : public ISimulator {
  friend class CompositeSimulator;
 
 public:
  IndividualSimulator(SimulatorType type =
#ifdef NO_QISKIT_AER
                          Simulators::SimulatorType::kQCSim
#else
                          Simulators::SimulatorType::kQiskitAer
#endif
                      ) noexcept
      : simulator(SimulatorsFactory::CreateSimulatorUnique(
            type, Simulators::SimulationType::kStatevector)) {
  }
 
  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);
        simulator->InitializeState(
            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) {
    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);
    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);
        simulator->InitializeState(
            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.reserve(nrBasisStates);
 
    for (Types::qubit_t state = 0; state < nrBasisStates; ++state)
      savedState.emplace_back(simulator->Amplitude(state));
  }
 
  void ClearSavedState() { savedState.clear(); }
 
  void RestoreState() override {
    if (!simulator) return;
    const size_t nrQubits = simulator->GetNumberOfQubits();
 
    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 {
    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 {
    std::unordered_map<Types::qubit_t, Types::qubit_t> res;
 
    const auto sc = simulator->SampleCounts(ConvertQubits(qubits), shots);
 
    for (auto [outcome, count] : sc)
      res[ConvertOutcomeFromLocal(outcome)] += count;
 
    return res;
  }
 
  std::unordered_map<std::vector<bool>, Types::qubit_t> SampleCountsMany(
      const Types::qubits_vector &qubits, size_t shots = 1000) override {
    std::unordered_map<std::vector<bool>, Types::qubit_t> res;
 
    const auto sc = simulator->SampleCountsMany(ConvertQubits(qubits), shots);
 
    for (auto [outcome, count] : sc)
      res[ConvertOutcomeFromLocal(outcome)] = count;
 
    return res;
  }
 
  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 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;
  }
 
  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)
      QCSimSimulator *qcsim = dynamic_cast<QCSimSimulator *>(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);
  }
 
 private:
  void InitializeAlias() {
    // TODO: implement it!
    if (GetType() == SimulatorType::kQCSim) {
      // qcsim - convert 'simulator' to qcsim simulator and access 'state' (from
      // there the statevector is accessible)
      QCSimSimulator *qcsim = dynamic_cast<QCSimSimulator *>(simulator.get());
 
      alias = std::unique_ptr<Utils::Alias>(
          new 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);
    simulator->InitializeState(
        newNrQubits, newAmplitudes);  // this will end up by swapping 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);
          }
 
          simulator->InitializeState(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