Composite.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _COMPOSITE_H_
#define _COMPOSITE_H_
 
#ifdef INCLUDED_BY_FACTORY
 
#include "Individual.h"
 
#include <vector>
 
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 : public ISimulator {
public:
  CompositeSimulator(SimulatorType type =
#ifdef NO_QISKIT_AER
                         Simulators::SimulatorType::kQCSim
#else
                         Simulators::SimulatorType::kQiskitAer
#endif
                     ) noexcept
      : type(type) {
    // 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 {
    qubitsMap.resize(nrQubits);
 
    // 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);
      sim->Initialize();
      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 {
    // don't allow chaning the method, it should stay statevector
    if (std::string("method") == key)
      return;
 
    for (auto &[id, simulator] : simulators)
      simulator->Configure(key, value);
  }

  std::string GetConfiguration(const char *key) const override {
    if (simulators.empty())
      return "";
 
    return simulators.begin()->second->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 {
    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;
  }

  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::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 {
    // 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;
  }

  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();
  }
 
  // 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() { 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 = Clone(); }

  void RestoreState() override {
    if (savedState) {
      const CompositeSimulator *savedStatePtr =
          static_cast<CompositeSimulator *>(savedState.get());
 
      type =
          savedStatePtr->type; 
      qubitsMap =
          savedStatePtr
              ->qubitsMap; 
 
      nrQubits =
          savedStatePtr->nrQubits; 
      nextId = savedStatePtr
                   ->nextId; 
      enableMultithreading =
          savedStatePtr
              ->enableMultithreading; 
 
      simulators.clear();
      for (auto &[id, simulator] : savedStatePtr->simulators) {
        auto isim = simulator->Clone();
        simulators[id] = std::unique_ptr<IndividualSimulator>(
            static_cast<IndividualSimulator *>(isim.release()));
      }
    }
  }

  Types::qubit_t MeasureNoCollapse() override {
    Types::qubit_t res = 0;
 
    for (auto &[id, simulator] : simulators) {
      const Types::qubit_t meas = simulator->MeasureNoCollapse();
      res |= meas;
    }
 
    return res;
  }

  std::unique_ptr<ISimulator> Clone() override {
    auto clone = std::make_unique<CompositeSimulator>(type);
 
    clone->type = type; 
    clone->qubitsMap =
        qubitsMap; 
 
    clone->nrQubits = nrQubits; 
    clone->nextId = nextId; 
    clone->enableMultithreading =
        enableMultithreading; 
 
    for (auto &[id, simulator] : simulators) {
      auto isim = simulator->Clone();
      clone->simulators[id] = std::unique_ptr<IndividualSimulator>(
          static_cast<IndividualSimulator *>(isim.release()));
    }
 
    if (savedState)
      clone->savedState = savedState->Clone();
 
    return clone;
  }
 
private:
  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;
 
    qubitsMap[qubit] = nextId; // the qubit will be in the new simulator
    simulators[nextId] = sim->Split(qubit, qubitOutcome, enableMultithreading);
 
    ++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; 
};
 
} // namespace Private
} // namespace Simulators
 
#endif
#endif