Circuit.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _CIRCUIT_H_
#define _CIRCUIT_H_
 
#define _USE_MATH_DEFINES
#include <math.h>
#include <set>
 
#include "Conditional.h"
#include "Delay.h"
#include "Operations.h"
#include "QuantumChannelOperation.h"
#include "QuantumGates.h"
#include "Reset.h"
#include <vector>
 
namespace Circuits {

template <typename Time = Types::time_type>
class Circuit : public IOperation<Time> {
 public:
  using ExecuteResults =
      std::unordered_map<std::vector<bool>,
                         size_t>; 
  using BitMapping =
      std::unordered_map<Types::qubit_t,
                         Types::qubit_t>; 
 
  using Operation = IOperation<Time>;              
  using OperationPtr = std::shared_ptr<Operation>; 
  using OperationsVector =
      std::vector<OperationPtr>; 
 
  using value_type = typename OperationsVector::value_type;
  using allocator_type = typename OperationsVector::allocator_type;
  using pointer = typename OperationsVector::pointer;
  using const_pointer = typename OperationsVector::const_pointer;
  using reference = typename OperationsVector::reference;
  using const_reference = typename OperationsVector::const_reference;
  using size_type = typename OperationsVector::size_type;
  using difference_type = typename OperationsVector::difference_type;
 
  using iterator = typename OperationsVector::iterator;
  using const_iterator = typename OperationsVector::const_iterator;
  using reverse_iterator = typename OperationsVector::reverse_iterator;
  using const_reverse_iterator =
      typename OperationsVector::const_reverse_iterator;

  Circuit(const OperationsVector &ops = {}) : Operation(), operations(ops) {}

  void Execute(const std::shared_ptr<Simulators::ISimulator>& sim,
               OperationState& state) const override {
    ExecuteBD(sim, state);
  }

  void ExecuteBD(const std::shared_ptr<Simulators::ISimulator> &sim,
               OperationState &state, size_t* curMaxBondDim = nullptr) const {
    state.Reset();
    if (!sim) return;
 
    for (const auto& op : operations) {
      op->Execute(sim, state);
      if (curMaxBondDim) {
        const auto bondDim = sim->GetCurrentMaxBondDimension();
        if (bondDim > *curMaxBondDim) *curMaxBondDim = bondDim;
      }
    }
    // sim->Flush();
  }

  OperationType GetType() const override { return OperationType::kComposite; }

  void AddOperation(const OperationPtr &op) { operations.push_back(op); }

  void Delay(Types::qubit_t qubit, Time duration) {
    operations.push_back(std::make_shared<Circuits::Delay<Time>>(qubit, duration));
  }

  void ReplaceOperation(size_t index, const OperationPtr &op) {
    if (index >= operations.size()) return;
    operations[index] = op;
  }

  void SetOperations(const OperationsVector &ops) { operations = ops; }

  void AddOperations(const OperationsVector &ops) {
    operations.insert(operations.end(), ops.begin(), ops.end());
  }

  void AddCircuit(const std::shared_ptr<Circuit<Time>> &circuit) {
    AddOperations(circuit->GetOperations());
  }

  const OperationsVector &GetOperations() const { return operations; }

  void Clear() { operations.clear(); }

  OperationPtr Clone() const override {
    OperationsVector newops;
 
    for (auto &op : operations) newops.emplace_back(op->Clone());
 
    return std::make_shared<Circuit<Time>>(newops);
  }

  OperationPtr CloneFlyweight() const {
    OperationsVector newops;
 
    for (auto &op : operations) newops.push_back(op);
 
    return std::make_shared<Circuit<Time>>(newops);
  }

  OperationPtr Remap(const BitMapping &qubitsMap,
                     const BitMapping &bitsMap = {}) const override {
    OperationsVector newops;
 
    for (const auto &op : operations)
      newops.emplace_back(op->Remap(qubitsMap, bitsMap));
 
    return std::make_shared<Circuit<Time>>(newops);
  }

  std::shared_ptr<Circuit<Time>> RemapToContinuous(BitMapping &newQubitsMap,
                                                   BitMapping &reverseBitsMap,
                                                   size_t &nrQubits,
                                                   size_t &nrCbits) const {
    OperationsVector newops;
 
    BitMapping newBitsMap;
 
    nrQubits = 0;
    nrCbits = 0;
 
    for (const auto &op : operations) {
      const auto affectedBits = op->AffectedBits();
      const auto affectedQubits = op->AffectedQubits();
 
      for (const auto qubit : affectedQubits) {
        const auto it = newQubitsMap.find(qubit);
        if (it == newQubitsMap.end()) {
          newQubitsMap[qubit] = nrQubits;
          ++nrQubits;
        }
      }
 
      for (const auto bit : affectedBits) {
        const auto it = newBitsMap.find(bit);
        if (it == newBitsMap.end()) {
          newBitsMap[bit] = nrCbits;
          reverseBitsMap[nrCbits] = bit;
          ++nrCbits;
        }
      }
 
      newops.emplace_back(op->Remap(newQubitsMap, newBitsMap));
    }
 
    return std::make_shared<Circuit<Time>>(newops);
  }

  static ExecuteResults RemapResultsBack(const ExecuteResults &results,
                                         const BitMapping &bitsMap = {},
                                         bool ignoreNotMapped = false,
                                         size_t sz = 0) {
    ExecuteResults newResults;
 
    if (!ignoreNotMapped && sz == 0) {
      for (const auto &[from, to] : bitsMap)
        if (to > sz) sz = to;
 
      ++sz;
    }
 
    for (const auto &res : results) {
      Circuits::OperationState mappedState(res.first);
 
      mappedState.Remap(bitsMap, ignoreNotMapped,
                        ignoreNotMapped ? bitsMap.size() : sz);
      newResults[mappedState.GetAllBits()] += res.second;
    }
 
    return newResults;
  }

  static void AccumulateResults(ExecuteResults &results,
                                const ExecuteResults &newResults) {
    for (const auto &res : newResults) results[res.first] += res.second;
  }

  static void AccumulateResultsWithRemapBack(ExecuteResults &results,
                                             const ExecuteResults &newResults,
                                             const BitMapping &bitsMap = {},
                                             bool ignoreNotMapped = true,
                                             size_t sz = 0) {
    if (!ignoreNotMapped && sz == 0) {
      for (const auto &[from, to] : bitsMap)
        if (to > sz) sz = to;
 
      ++sz;
    }
 
    for (const auto &res : newResults) {
      Circuits::OperationState mappedState(res.first);
 
      mappedState.Remap(bitsMap, ignoreNotMapped,
                        ignoreNotMapped ? bitsMap.size() : sz);
      results[mappedState.GetAllBits()] += res.second;
      /*
      const auto it = bitsMap.find(res.first);
      if (it != bitsMap.end())
              results[it->second] += res.second;
      */
    }
  }
 
  // TODO: This converts all swap, cswap and ccnot gates
  // it's not really needed to convert them all,
  // only those that are not applied locally, that is, on a single host
  // the local ones can remain as they are
  // so use network topology to decide which ones to convert
  // that's for later, when we have the network part implemented
  // and also the code that splits the circuit

  void ConvertForDistribution() {
    // this will make the circuit better for distributed computing
    // TODO: if composite operations will be implemented, those need to be
    // optimized as well
 
    ReplaceThreeQubitAndSwapGates();
  }

  void ConvertForCutting() { ReplaceThreeQubitAndSwapGates(true); }
 
  /*
   * @brief Splits the measurements to measurements on individual qubits and
   * tries to order them as needed by the following clasically conditional
   * gates.
   *
   * Splits the measurements to measurements on individual qubits and tries to
   * order them as needed by the following clasically conditional gates. This is
   * needed for netqasm, which requires the measurements to be in the right
   * order, if the following clasically conditional gates need sending to
   * another host. This wouldn't be needed if they are local, but we don't know
   * that at this point.
   *
   * So in short, all measurements on more than one qubit are converted to one
   * qubit measurements, and then all measurements that are grouped together
   * (not separated by some other operation) are ordered in the same order as
   * the conditions on the following clasically conditional gates.
   */
  void EnsureProperOrderForMeasurements() {
    // TODO: Maybe this should be moved at the netqasm level circuit conversion,
    // since it's not needed for all kinds of distributed computing currently
    // there is a virtual function in controller that does nothing in the base
    // class, but for netqasm it calls this method
    OperationsVector newops;
 
    for (size_t i = 0; i < operations.size(); ++i) {
      const auto op = operations[i];
      if (op->GetType() != OperationType::kMeasurement) {
        newops.emplace_back(op);
        continue;
      }
 
      // ok, if it's a measurement, look ahead, accumulate all measurements and
      // then add them in the right order
      std::unordered_set<size_t> bits;
      std::unordered_map<size_t, Types::qubit_t> measQubits;
      std::unordered_map<size_t, Time> measDelays;
 
      auto affectedBits = op->AffectedBits();
      auto affectedQubits = op->AffectedQubits();
 
      for (size_t q = 0; q < affectedQubits.size(); ++q) {
        bits.insert(affectedBits[q]);
        measQubits[affectedBits[q]] = affectedQubits[q];
        measDelays[affectedBits[q]] = op->GetDelay();
      }
 
      size_t j = i + 1;
      for (; j < operations.size(); ++j) {
        const auto op2 = operations[j];
 
        if (op2->GetType() != OperationType::kMeasurement) break;
 
        affectedQubits = op2->AffectedQubits();
 
        const auto meas =
            std::static_pointer_cast<MeasurementOperation<Time>>(op2);
        affectedBits = meas->GetBitsIndices();
        for (size_t q = 0; q < affectedBits.size(); ++q) {
          bits.insert(affectedBits[q]);
          measQubits[affectedBits[q]] = affectedQubits[q];
          measDelays[affectedBits[q]] = op2->GetDelay();
        }
      }
 
      i = j - 1;
 
      // the right order is the one following in the classically controlled
      // gates
      for (; j < operations.size(); ++j) {
        const auto op2 = operations[j];
        if (op2->GetType() == OperationType::kConditionalGate ||
            op2->GetType() == OperationType::kConditionalMeasurement ||
            op2->GetType() == OperationType::kConditionalRandomGen) {
          auto condop =
              std::static_pointer_cast<IConditionalOperation<Time>>(op2);
          const auto condbits = condop->AffectedBits();
          for (const auto bit : condbits)
            if (bits.find(bit) != bits.end()) {
              newops.emplace_back(std::make_shared<MeasurementOperation<Time>>(
                  std::vector{std::make_pair(measQubits[bit], bit)},
                  measDelays[bit]));
              bits.erase(bit);
            }
        }
        if (bits.empty()) break;
      }
 
      // now add the measurements that were left in any order
      for (auto bit : bits)
        newops.emplace_back(std::make_shared<MeasurementOperation<Time>>(
            std::vector{std::make_pair(measQubits[bit], bit)},
            measDelays[bit]));
    }
 
    operations.swap(newops);
  }

  size_t GetMaxQubitIndex() const {
    size_t mx = 0;
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      for (auto q : qbits)
        if (q > mx) mx = q;
    }
 
    return mx;
  }

  size_t GetMinQubitIndex() const {
    size_t mn = std::numeric_limits<size_t>::max();
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      for (auto q : qbits)
        if (q < mn) mn = q;
    }
 
    return mn;
  }

  size_t GetMaxCbitIndex() const {
    size_t mx = 0;
    for (const auto &op : operations) {
      const auto cbits = op->AffectedBits();
      for (auto q : cbits)
        if (q > mx) mx = q;
    }
 
    return mx;
  }

  size_t GetMinCbitIndex() const {
    size_t mn = std::numeric_limits<size_t>::max();
    for (const auto &op : operations) {
      const auto cbits = op->AffectedBits();
      for (auto q : cbits)
        if (q < mn) mn = q;
    }
 
    return mn;
  }

  std::set<size_t> GetQubits() const {
    std::set<size_t> qubits;
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      qubits.insert(qbits.begin(), qbits.end());
    }
 
    return qubits;
  }

  std::set<size_t> GetBits() const {
    std::set<size_t> cbits;
    for (const auto &op : operations) {
      const auto bits = op->AffectedBits();
      cbits.insert(bits.begin(), bits.end());
    }
 
    return cbits;
  }

  Types::qubits_vector AffectedQubits() const override {
    auto qubits = GetQubits();
 
    Types::qubits_vector qubitsVec;
    qubitsVec.reserve(qubits.size());
 
    for (auto q : qubits) qubitsVec.emplace_back(q);
 
    return qubitsVec;
  }

  std::vector<size_t> AffectedBits() const override {
    auto bits = GetBits();
 
    std::vector<size_t> bitsVec;
    bitsVec.reserve(bits.size());
 
    for (auto b : bits) bitsVec.emplace_back(b);
 
    return bitsVec;
  }

  bool NeedsEntanglementForDistribution() const override {
    for (const auto &op : operations)
      if (op->NeedsEntanglementForDistribution()) return true;
 
    return false;
  }

  bool CanAffectQuantumState() const override {
    for (const auto &op : operations)
      if (op->CanAffectQuantumState()) return true;
 
    return false;
  }

  std::unordered_map<size_t, OperationPtr> GetLastOperationsOnQubits() const {
    std::unordered_map<size_t, OperationPtr> lastOps;
 
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      for (auto q : qbits) lastOps[q] = op;
    }
 
    return lastOps;
  }

  std::unordered_map<size_t, OperationPtr> GetFirstOperationsOnQubits() const {
    std::unordered_map<size_t, OperationPtr> firstOps;
 
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      for (auto q : qbits) {
        if (firstOps.find(q) == firstOps.end()) firstOps[q] = op;
      }
    }
 
    return firstOps;
  }

  void AddResetsIfNeeded(Time delay = 0) {
    const auto GetLastOps = GetLastOperationsOnQubits();
 
    for (const auto &[q, op] : GetLastOps)
      if (op->GetType() !=
          OperationType::kReset)  // don't add it if there is already a reset
                                  // operation on the qubit
        operations.emplace_back(
            std::make_shared<Reset<Time>>(Types::qubits_vector{q}, delay));
  }

  void AddResetsAtBeginningIfNeeded(Time delay = 0) {
    const auto GetFirstOps = GetFirstOperationsOnQubits();
 
    for (const auto &[q, op] : GetFirstOps)
      if (op->GetType() !=
          OperationType::kReset)  // don't add it if there is already a reset
                                  // operation on the qubit
        operations.insert(
            operations.begin(),
            std::make_shared<Reset<Time>>(Types::qubits_vector{q}, delay));
  }

  void Optimize(bool optimizeRotationGates = true) {
    // Some ideas, from simple to more complex:
    //
    // IMPORTANT: Focus on the gates that are added for distributed computing,
    // either for the one with entanglement or the one with cutting the reason
    // is that maybe the provided circuit is not that bad, but due of the
    // supplementary gates added, duplicates (for example) will occur the most
    // important one qubit ones added are hadamard then X, S, Sdag and Z
    //
    // 1. First, one qubit gates can be combined into a single gate or even no
    // gate a) Straightforward for those that are their own inverse (that is,
    // hermitian/involutory): Hadamard and Pauli gates for example, if one finds
    // two of them in sequence, they can be removed other ones are the one that
    // are followed by their 'dag' (inverse, since the gates are unitary) in the
    // circuit, those can be removed, too
 
    // several resets can be also changed into a single one, also repeated
    // measurements of the same qubit, with result in the same cbit can be
    // replaced by a single measurement
 
    // b) other ones can be combined into a single gate, for example phase shift
    // gates or rotation gates two phase gates (not the general phase shift we
    // have, but the one with 1, i on the diagonal) can be combined into a Z
    // gate, for example, or two sqrtNot gates can be combined into a single X
    // gate combinations of phase gates and hadamard can be replaced by pauli
    // gates in some cases, and so on even the U gate could be used to join
    // together several one qubit gates
 
    // c) even more complex... three or more one qubit gates could be combined
    // into a single gate... an example is HXH = Z other examples SXS^t = Y,
    // SZS^t = Z
 
    // 2. Two qubit gates can be optimized, too
    // for example two CNOT gates in sequence can be removed if the control
    // qubit is the same, the same goes for two CZ or CY or SWAP gates (this
    // goes for the three qubit gates, CCX, CCY, CCZ, CSWAP, too)
 
    // 3. Some gates commute, the reorder can give some opportunities for more
    // optimization
 
    // 4. Groups of gates
    // the possibilities are endless, but might be easier to focus first on
    // Clifford gates for example a X sandwiched between CNOTs can be replaced
    // with two X on each qubit if the original X is on the control qubit or
    // with an X on the target qubit if the original X is on the target qubit a
    // similar thing happens if Z is sandwiched between CNOTs, but this time Z x
    // Z appears if the original Z is on the target qubit and if original Z is
    // on the control qubit, then the CNOTs dissapear and the Z remains on the
    // control qubit
 
    // three CNOT gates with the one in the middle turned upside down compare
    // with the other two can be replaced by a single swap gate
 
    // first stage, take out duplicates of H, X, Y, Z
 
    bool changed;
 
    do {
      changed = false;
 
      std::vector<std::shared_ptr<IOperation<Time>>> newops;
      newops.reserve(operations.size());
 
      for (int i = 0; i < static_cast<int>(operations.size()); ++i) {
        const std::shared_ptr<IOperation<Time>> &op = operations[i];
 
        const auto type = op->GetType();
        if (type == OperationType::kNoOp)
          continue;
        else if (type == OperationType::kGate) {
          std::shared_ptr<IQuantumGate<Time>> gate =
              std::static_pointer_cast<IQuantumGate<Time>>(op);
          const auto qubits = gate->AffectedQubits();
 
          if (qubits.size() == 1) {
            // TODO: HXH = Z, SXS^t = Y, SZS^t = Z ????
 
            auto gateType = gate->GetGateType();
            bool replace = false;
 
            // if it's one of the interesting gates, look ahead to see if it's
            // followed by the same gate on the same qubit if yes, replace the
            // next one with a nop and skip the current one (or replace the pair
            // with a single gate, depending on the type) set changed to true if
            // something was changed
            switch (gateType) {
              case QuantumGateType::kRxGateType:
                [[fallthrough]];
              case QuantumGateType::kRyGateType:
                [[fallthrough]];
              case QuantumGateType::kRzGateType:
                if (!optimizeRotationGates) {
                  newops.push_back(op);
                  break;
                }
                [[fallthrough]];
              case QuantumGateType::kPhaseGateType:
                replace = true;
                [[fallthrough]];
              // those above will be replaced the pair with a single gate, all
              // the following are the ones that get removed
              case QuantumGateType::kHadamardGateType:
                [[fallthrough]];
              case QuantumGateType::kXGateType:
                [[fallthrough]];
              case QuantumGateType::kYGateType:
                [[fallthrough]];
              case QuantumGateType::kZGateType:
                [[fallthrough]];
              case QuantumGateType::kKGateType:
                [[fallthrough]];
              case QuantumGateType::kSGateType:
                [[fallthrough]];
              case QuantumGateType::kSdgGateType:
                [[fallthrough]];
              case QuantumGateType::kTGateType:
                [[fallthrough]];
              case QuantumGateType::kTdgGateType:
                [[fallthrough]];
              case QuantumGateType::kSxGateType:
                [[fallthrough]];
              case QuantumGateType::kSxDagGateType: {
                bool found = false;
 
                if (gateType == QuantumGateType::kSGateType)
                  gateType = QuantumGateType::kSdgGateType;
                else if (gateType == QuantumGateType::kSdgGateType)
                  gateType = QuantumGateType::kSGateType;
                else if (gateType == QuantumGateType::kTGateType)
                  gateType = QuantumGateType::kTdgGateType;
                else if (gateType == QuantumGateType::kTdgGateType)
                  gateType = QuantumGateType::kTGateType;
                else if (gateType == QuantumGateType::kSxGateType)
                  gateType = QuantumGateType::kSxDagGateType;
                else if (gateType == QuantumGateType::kSxDagGateType)
                  gateType = QuantumGateType::kSxGateType;
 
                for (size_t j = i + 1; j < operations.size(); ++j) {
                  auto &nextOp = operations[j];
                  if (!nextOp->CanAffectQuantumState()) continue;
 
                  const auto nextQubits = nextOp->AffectedQubits();
                  bool hasQubit = false;
 
                  for (auto q : nextQubits)
                    if (q == qubits[0]) {
                      hasQubit = true;
                      break;
                    }
 
                  if (!hasQubit)
                    continue;  // an op that does not touch the current qubit
                               // can be skipped
                  else if (nextQubits.size() != 1)
                    break;  // if it touches the current qubit and it's
                            // something else than a single qubit gate, stop
 
                  const auto nextType = nextOp->GetType();
                  if (nextType != OperationType::kGate)
                    break;  // could be a classically conditioned gate, stop
 
                  const auto &nextGate =
                      std::static_pointer_cast<SingleQubitGate<Time>>(nextOp);
                  if (nextGate->GetGateType() == gateType) {
                    if (replace) {
                      const auto params1 = gate->GetParams();
                      const auto params2 = nextGate->GetParams();
 
                      const double param = params1[0] + params2[0];
                      const auto delay =
                          gate->GetDelay() + nextGate->GetDelay();
 
                      if (gateType == QuantumGateType::kPhaseGateType)
                        newops.push_back(std::make_shared<PhaseGate<Time>>(
                            qubits[0], param, delay));
                      else if (gateType == QuantumGateType::kRxGateType)
                        newops.push_back(std::make_shared<RxGate<Time>>(
                            qubits[0], param, delay));
                      else if (gateType == QuantumGateType::kRyGateType)
                        newops.push_back(std::make_shared<RyGate<Time>>(
                            qubits[0], param, delay));
                      else
                        newops.push_back(std::make_shared<RzGate<Time>>(
                            qubits[0], param, delay));
                    }
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if ((gateType == QuantumGateType::kSGateType &&
                              nextGate->GetGateType() ==
                                  QuantumGateType::kSdgGateType) ||
                             (gateType == QuantumGateType::kSdgGateType &&
                              nextGate->GetGateType() ==
                                  QuantumGateType::kSGateType)) {
                    // if expecting an S gate (or a Sdg gate) and found the
                    // original one instead, replace the pair with a Z gate (S *
                    // S = Z, Sdag * Sdag = Z)
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    newops.push_back(
                        std::make_shared<ZGate<Time>>(qubits[0], delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if ((gateType == QuantumGateType::kSxGateType &&
                              nextGate->GetGateType() ==
                                  QuantumGateType::kSxDagGateType) ||
                             (gateType == QuantumGateType::kSxDagGateType &&
                              nextGate->GetGateType() ==
                                  QuantumGateType::kSxGateType)) {
                    // if expecting an S gate (or a Sdg gate) and found the
                    // original one instead, replace the pair with a X gate (Sx
                    // * Sx = X, SXdag * SXdag = X)
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    newops.push_back(
                        std::make_shared<XGate<Time>>(qubits[0], delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if (gateType == QuantumGateType::kTGateType &&
                             nextGate->GetGateType() ==
                                 QuantumGateType::kTdgGateType) {
                    // if expecting a T gate and found the Tdgate instead,
                    // replace the pair with a Sdag gate (Tdg * Tdg = Sdag)
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    newops.push_back(
                        std::make_shared<SdgGate<Time>>(qubits[0], delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if (gateType == QuantumGateType::kTdgGateType &&
                             nextGate->GetGateType() ==
                                 QuantumGateType::kTGateType) {
                    // if expecting a Tdg gate and found the T gate instead,
                    // replace the pair with a S gate (T * T = S)
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    newops.push_back(
                        std::make_shared<SGate<Time>>(qubits[0], delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if (gateType == QuantumGateType::kPhaseGateType &&
                             (nextGate->GetGateType() ==
                                  QuantumGateType::kSGateType ||
                              nextGate->GetGateType() ==
                                  QuantumGateType::kSdgGateType ||
                              nextGate->GetGateType() ==
                                  QuantumGateType::kTGateType ||
                              nextGate->GetGateType() ==
                                  QuantumGateType::kTdgGateType)) {
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    double param2;
                    if (nextGate->GetGateType() == QuantumGateType::kSGateType)
                      param2 = 0.5 * M_PI;
                    else if (nextGate->GetGateType() ==
                             QuantumGateType::kSdgGateType)
                      param2 = -0.5 * M_PI;
                    else if (nextGate->GetGateType() ==
                             QuantumGateType::kTGateType)
                      param2 = 0.25 * M_PI;
                    else
                      param2 = -0.25 * M_PI;
 
                    const auto param = gate->GetParams()[0] + param2;
                    newops.push_back(std::make_shared<PhaseGate<Time>>(
                        qubits[0], param, delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else if (nextGate->GetGateType() ==
                                 QuantumGateType::kPhaseGateType &&
                             (gateType == QuantumGateType::kSGateType ||
                              gateType == QuantumGateType::kSdgGateType ||
                              gateType == QuantumGateType::kTGateType ||
                              gateType == QuantumGateType::kTdgGateType)) {
                    const auto delay = gate->GetDelay() + nextGate->GetDelay();
                    double param1;
                    if (gateType == QuantumGateType::kSGateType)
                      param1 = -0.5 * M_PI;
                    else if (gateType == QuantumGateType::kSdgGateType)
                      param1 = 0.5 * M_PI;
                    else if (gateType == QuantumGateType::kTGateType)
                      param1 = -0.25 * M_PI;
                    else
                      param1 = 0.25 * M_PI;
 
                    const auto param = nextGate->GetParams()[0] + param1;
                    newops.push_back(std::make_shared<PhaseGate<Time>>(
                        qubits[0], param, delay));
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else
                    break;  // not the expected gate, acting on same qubit, bail
                            // out
                }
 
                if (!found) newops.push_back(op);
              } break;
              default:
                // if no, just add it
                newops.push_back(op);
                break;
            }
          } else if (qubits.size() == 2) {
            auto gateType = gate->GetGateType();
            bool replace = false;
 
            // if it's one of the interesting gates, look ahead to see if it's
            // followed by the same gate on the same qubit if yes, replace the
            // next one with a nop and skip the current one (or replace the pair
            // with a single gate, depending on the type) set changed to true if
            // something was changed
            switch (gateType) {
              case QuantumGateType::kCRxGateType:
                [[fallthrough]];
              case QuantumGateType::kCRyGateType:
                [[fallthrough]];
              case QuantumGateType::kCRzGateType:
                if (!optimizeRotationGates) {
                  newops.push_back(op);
                  break;
                }
                [[fallthrough]];
              case QuantumGateType::kCPGateType:
                replace = true;
                [[fallthrough]];
                // those above will be replaced the pair with a single gate, all
                // the following are the ones that get removed
              case QuantumGateType::kCXGateType:
                [[fallthrough]];
              case QuantumGateType::kCYGateType:
                [[fallthrough]];
              case QuantumGateType::kCZGateType:
                [[fallthrough]];
              case QuantumGateType::kCHGateType:
                [[fallthrough]];
              case QuantumGateType::kSwapGateType:
                [[fallthrough]];
              case QuantumGateType::kCSxGateType:
                [[fallthrough]];
              case QuantumGateType::kCSxDagGateType: {
                bool found = false;
 
                if (gateType == QuantumGateType::kCSxGateType)
                  gateType = QuantumGateType::kCSxDagGateType;
                else if (gateType == QuantumGateType::kCSxDagGateType)
                  gateType = QuantumGateType::kCSxGateType;
 
                // looking forward for the next operation that acts on the same
                // qubits
                for (size_t j = i + 1; j < operations.size(); ++j) {
                  auto &nextOp = operations[j];
                  if (!nextOp->CanAffectQuantumState()) continue;
 
                  const auto nextQubits = nextOp->AffectedQubits();
 
                  bool hasQubit = false;
 
                  for (auto q : nextQubits)
                    if (q == qubits[0] || q == qubits[1]) {
                      hasQubit = true;
                      break;
                    }
 
                  if (!hasQubit)
                    continue;  // an op that does not touch the current qubit
                               // can be skipped
                  else if (nextQubits.size() != 2)
                    break;  // if it touches a current qubit and it's something
                            // else than a two qubits gate, stop
                  // if it's not the same qubits, bail out
                  else if (gateType == QuantumGateType::kSwapGateType &&
                           !((qubits[0] == nextQubits[0] &&
                              qubits[1] == nextQubits[1]) ||
                             (qubits[0] == nextQubits[1] &&
                              qubits[1] == nextQubits[0])))
                    break;
                  else if (!(qubits[0] == nextQubits[0] &&
                             qubits[1] == nextQubits[1]))
                    break;
 
                  const auto nextType = nextOp->GetType();
                  if (nextType != OperationType::kGate)
                    break;  // could be a classically conditioned gate, stop
 
                  const auto &nextGate =
                      std::static_pointer_cast<TwoQubitsGate<Time>>(nextOp);
                  if (nextGate->GetGateType() == gateType) {
                    if (replace) {
                      const auto params1 = gate->GetParams();
                      const auto params2 = nextGate->GetParams();
                      const double param = params1[0] + params2[0];
                      const auto delay =
                          gate->GetDelay() + nextGate->GetDelay();
 
                      if (gateType == QuantumGateType::kCPGateType)
                        newops.push_back(std::make_shared<CPGate<Time>>(
                            qubits[0], qubits[1], param, delay));
                      else if (gateType == QuantumGateType::kCRxGateType)
                        newops.push_back(std::make_shared<CRxGate<Time>>(
                            qubits[0], qubits[1], param, delay));
                      else if (gateType == QuantumGateType::kCRyGateType)
                        newops.push_back(std::make_shared<CRyGate<Time>>(
                            qubits[0], qubits[1], param, delay));
                      else
                        newops.push_back(std::make_shared<CRzGate<Time>>(
                            qubits[0], qubits[1], param, delay));
                    }
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;  // continue merging gates, we found one
                                     // that was merged/removed
                    found = true;    // don't put op in the new operations, we
                                     // handled it
                    break;
                  } else
                    break;  // not the expected gate, acting on same qubits,
                            // bail out
                }           // end for of looking forward
 
                if (!found) newops.push_back(op);
              } break;
              default:
                // if no, just add it
                newops.push_back(op);
                break;
            }
          } else if (qubits.size() == 3) {
            auto gateType = gate->GetGateType();
 
            // if it's one of the interesting gates, look ahead to see if it's
            // followed by the same gate on the same qubit if yes, replace the
            // next one with a nop and skip the current one (or replace the pair
            // with a single gate, depending on the type) set changed to true if
            // something was changed
            switch (gateType) {
              case QuantumGateType::kCSwapGateType:
                [[fallthrough]];
              case QuantumGateType::kCCXGateType: {
                bool found = false;
 
                for (size_t j = i + 1; j < operations.size(); ++j) {
                  auto &nextOp = operations[j];
                  if (!nextOp->CanAffectQuantumState()) continue;
 
                  const auto nextQubits = nextOp->AffectedQubits();
 
                  bool hasQubit = false;
 
                  for (auto q : nextQubits)
                    if (q == qubits[0] || q == qubits[1] || q == qubits[2]) {
                      hasQubit = true;
                      break;
                    }
 
                  if (!hasQubit)
                    continue;  // an op that does not touch the current qubit
                               // can be skipped
                  else if (nextQubits.size() != 3)
                    break;  // if it touches a current qubit and it's something
                            // else than a three qubits gate, stop
                  // if it's not the same qubits, bail out
                  else if (gateType == QuantumGateType::kCSwapGateType &&
                           (qubits[0] != nextQubits[0] ||
                            !((qubits[1] == nextQubits[1] &&
                               qubits[2] == nextQubits[2]) ||
                              (qubits[1] == nextQubits[2] &&
                               qubits[2] == nextQubits[1]))))
                    break;
                  else if (gateType == QuantumGateType::kCCXGateType &&
                           (qubits[2] != nextQubits[2] ||
                            !(qubits[1] == nextQubits[1] &&
                              qubits[2] == nextQubits[2]) ||
                            !(qubits[1] == nextQubits[2] &&
                              qubits[2] == nextQubits[1])))
                    break;
 
                  const auto nextType = nextOp->GetType();
                  if (nextType != OperationType::kGate)
                    break;  // could be a classically conditioned gate, stop
 
                  const auto &nextGate =
                      std::static_pointer_cast<ThreeQubitsGate<Time>>(nextOp);
                  if (nextGate->GetGateType() == gateType) {
                    nextOp = std::make_shared<NoOperation<Time>>();
                    changed = true;
                    found = true;
                    break;
                  } else
                    break;  // not the expected gate, acting on same qubits,
                            // bail out
                }
 
                if (!found) newops.push_back(op);
              } break;
              default:
                // if no, just add it
                newops.push_back(op);
                break;
            }
          } else
            newops.push_back(op);
        } else
          newops.push_back(op);
      }  // end for on circuit operations
 
      operations.swap(newops);
    } while (changed);
  }

  void MoveMeasurementsAndResets() {
    OperationsVector newops;
    newops.reserve(operations.size());
 
    size_t qubitsNo = std::max(GetMaxQubitIndex(), GetMaxCbitIndex()) + 1;
 
    std::unordered_map<Types::qubit_t, std::vector<OperationPtr>> qubitOps;
 
    std::vector<OperationPtr> lastOps(qubitsNo);
 
    std::unordered_map<OperationPtr, std::unordered_set<OperationPtr>>
        dependenciesMap;
 
    for (const auto &op : operations) {
      std::unordered_set<OperationPtr> dependencies;
 
      const auto cbits = op->AffectedBits();
      for (auto c : cbits) {
        const auto lastOp = lastOps[c];
        if (lastOp) dependencies.insert(lastOp);
      }
 
      const auto qubits = op->AffectedQubits();
      for (auto q : qubits) {
        qubitOps[q].push_back(op);
 
        const auto lastOp = lastOps[q];
        if (lastOp) dependencies.insert(lastOp);
 
        lastOps[q] = op;
      }
 
      for (auto c : cbits) lastOps[c] = op;
 
      dependenciesMap[op] = dependencies;
    }
    lastOps.clear();
 
    std::vector<Types::qubit_t> indices(qubitsNo, 0);
 
    while (!dependenciesMap.empty()) {
      OperationPtr nextOp;
 
      // try to locate a 'next' gate for a qubit that is either a measurement or
      // a reset
      for (size_t q = 0; q < qubitsNo; ++q) {
        if (qubitOps.find(q) ==
            qubitOps.end())  // no operation left on this qubit
          continue;
 
        // grab the current operation for this qubit
        const auto &ops = qubitOps[q];
        const auto &op = ops[indices[q]];
 
        // consider only measurements and resets
        if (op->GetType() == OperationType::kMeasurement ||
            op->GetType() == OperationType::kReset) {
          bool hasDependencies = false;
 
          for (const auto &opd : dependenciesMap[op])
            if (dependenciesMap.find(opd) != dependenciesMap.end()) {
              hasDependencies = true;
              break;
            }
 
          if (!hasDependencies) {
            nextOp = op;
            break;
          }
        }
      }
 
      if (nextOp) {
        dependenciesMap.erase(nextOp);
 
        const auto qubits = nextOp->AffectedQubits();
        for (auto q : qubits) {
          ++indices[q];
          if (indices[q] >= qubitOps[q].size()) qubitOps.erase(q);
        }
 
        newops.emplace_back(std::move(nextOp));
        continue;
      }
 
      // if there is no measurement or reset, add the next gate
      for (Types::qubit_t q = 0; q < qubitsNo; ++q) {
        if (qubitOps.find(q) ==
            qubitOps.end())  // no operation left on this qubit
          continue;
 
        // grab the current operation for this qubit
        const auto &ops = qubitOps[q];
        const auto &op = ops[indices[q]];
 
        bool hasDependencies = false;
 
        for (const auto &opd : dependenciesMap[op])
          if (dependenciesMap.find(opd) != dependenciesMap.end()) {
            hasDependencies = true;
            break;
          }
 
        if (!hasDependencies) {
          nextOp = op;
          break;
        }
      }
 
      if (nextOp) {
        dependenciesMap.erase(nextOp);
 
        const auto qubits = nextOp->AffectedQubits();
        for (auto q : qubits) {
          ++indices[q];
          if (indices[q] >= qubitOps[q].size()) qubitOps.erase(q);
        }
 
        newops.emplace_back(std::move(nextOp));
      }
    }
 
    assert(newops.size() == operations.size());
 
    operations.swap(newops);
  }

  std::pair<std::vector<size_t>, std::vector<Time>> GetDepth() const {
    size_t maxDepth;
    Time maxTime;
 
    size_t qubitsNo = GetMaxQubitIndex() + 1;
    std::vector<Time> qubitTimes(qubitsNo, 0);
    std::vector<size_t> qubitDepths(qubitsNo, 0);
 
    std::unordered_map<size_t, size_t> fromQubits;
 
    for (const auto &op : operations) {
      const auto qbits = op->AffectedQubits();
      const auto delay = op->GetDelay();
 
      maxTime = 0;
      maxDepth = 0;
      for (auto q : qbits) {
        qubitTimes[q] += delay;
        ++qubitDepths[q];
        if (qubitTimes[q] > maxTime) maxTime = qubitTimes[q];
        if (qubitDepths[q] > maxDepth) maxDepth = qubitDepths[q];
      }
 
      const auto t = op->GetType();
      std::vector<size_t> condbits;
 
      // TODO: deal with 'random gen' operations, those do not affect qubits
      // directly, but they do affect the classical bits and can be used in
      // conditional operations
 
      if (t == OperationType::kConditionalGate ||
          t == OperationType::kConditionalMeasurement) {
        condbits = op->AffectedBits();
 
        for (auto bit : condbits) {
          if (fromQubits.find(bit) != fromQubits.end()) bit = fromQubits[bit];
 
          bool found = false;
          for (auto q : qbits) {
            if (q == bit) {
              found = true;
              break;
            }
          }
          if (found || bit >= qubitsNo) continue;
 
          qubitTimes[bit] += delay;
          ++qubitDepths[bit];
          if (qubitTimes[bit] > maxTime) maxTime = qubitTimes[bit];
          if (qubitDepths[bit] > maxDepth) maxDepth = qubitDepths[bit];
        }
 
        if (t == OperationType::kConditionalMeasurement) {
          const auto condMeas =
              std::static_pointer_cast<ConditionalMeasurement<Time>>(op);
          const auto meas = condMeas->GetOperation();
          const auto measQubits = meas->AffectedQubits();
          const auto measBits = meas->AffectedBits();
 
          for (size_t i = 0; i < measQubits.size(); ++i) {
            if (i < measBits.size())
              fromQubits[measBits[i]] = measQubits[i];
            else
              fromQubits[measQubits[i]] = measQubits[i];
          }
        }
      } else if (t == OperationType::kMeasurement) {
        condbits = op->AffectedBits();
 
        for (size_t i = 0; i < qbits.size(); ++i) {
          if (i < condbits.size())
            fromQubits[condbits[i]] = qbits[i];
          else
            fromQubits[qbits[i]] = qbits[i];
        }
      }
 
      for (auto q : qbits) {
        qubitTimes[q] = maxTime;
        qubitDepths[q] = maxDepth;
      }
 
      for (auto bit : condbits) {
        qubitTimes[bit] = maxTime;
        qubitDepths[bit] = maxDepth;
      }
    }
 
    return std::make_pair(qubitDepths, qubitTimes);
  }

  std::pair<size_t, Time> GetMaxDepth() const {
    auto [qubitDepths, qubitTimes] = GetDepth();
 
    Time maxTime = 0;
    size_t maxDepth = 0;
    for (size_t qubit = 0; qubit < qubitDepths.size(); ++qubit) {
      if (qubitTimes[qubit] > maxTime) maxTime = qubitTimes[qubit];
      if (qubitDepths[qubit] > maxDepth) maxDepth = qubitDepths[qubit];
    }
 
    return std::make_pair(maxDepth, maxTime);
  }

  size_t GetNumberOfOperations() const { return operations.size(); }

  OperationPtr GetOperation(size_t pos) const {
    if (pos >= operations.size()) return nullptr;
 
    return operations[pos];
  }

  std::shared_ptr<Circuit<Time>> GetCircuitCut(Types::qubit_t startQubit,
                                               Types::qubit_t endQubit) const {
    OperationsVector newops;
    newops.reserve(operations.size());
 
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
      bool containsOutsideQubits = false;
      bool containsInsideQubits = false;
      for (const auto q : qubits) {
        if (q < startQubit || q > endQubit) {
          containsOutsideQubits = true;
          if (containsInsideQubits) break;
        } else {
          containsInsideQubits = true;
          if (containsOutsideQubits) break;
        }
      }
 
      if (containsInsideQubits) {
        if (containsOutsideQubits)
          throw std::runtime_error(
              "Cannot cut the circuit with the specified interval");
        newops.emplace_back(op->Clone());
      }
    }
 
    return std::make_shared<Circuit<Time>>(newops);
  }

  bool HasOpsAfterMeasurements() const {
    std::unordered_set<Types::qubit_t> measuredQubits;
    std::unordered_set<Types::qubit_t> affectedQubits;
    std::unordered_set<Types::qubit_t> resetQubits;
 
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
 
      if (op->GetType() == OperationType::kMeasurement) {
        for (const auto qbit : qubits)
          if (resetQubits.find(qbit) !=
              resetQubits.end())  // there is a reset on this qubit already and
                                  // it's not at the beginning of the circuit
            return true;
 
        /*
        const auto bits = op->AffectedBits();
        if (bits.size() != qubits.size())
                return true;
 
        for (size_t b = 0; b < bits.size(); ++b)
                if (bits[b] != qubits[b])
                        return true;
        */
        measuredQubits.insert(qubits.begin(), qubits.end());
      } else if (op->GetType() == OperationType::kConditionalGate ||
                 op->GetType() == OperationType::kConditionalMeasurement ||
                 op->GetType() == OperationType::kRandomGen ||
                 op->GetType() == OperationType::kConditionalRandomGen)
        return true;
      else if (op->GetType() == OperationType::kReset) {
        // resets in the middle of the circuit are treated as measurements
        for (const auto qbit : qubits) {
          // if there is already a gate applied on the qubit but no measurement
          // yet, it's considered in the middle
          // if there is no gate applied,
          // then it's the first operation on the qubit
          if (affectedQubits.find(qbit) != affectedQubits.end() ||
              measuredQubits.find(qbit) != measuredQubits.end())
            resetQubits.insert(qbit);
 
          affectedQubits.insert(qbit);
        }
      } else {
        for (const auto qbit : qubits) {
          if (measuredQubits.find(qbit) !=
              measuredQubits
                  .end())  // there is a measurement on this qubit already
            return true;
 
          if (resetQubits.find(qbit) !=
              resetQubits.end())  // there is a reset on this qubit already and
                                  // it's not at the beginning of the circuit
            return true;
 
          affectedQubits.insert(qbit);
        }
      }
    }
 
    return false;
  }

  std::vector<bool> ExecuteNonMeasurements(
      const std::shared_ptr<Simulators::ISimulator> &sim,
      OperationState &state, size_t* curMaxBondDim = nullptr) const {
    std::vector<bool> executedOps;
    executedOps.reserve(operations.size());
 
    std::unordered_set<Types::qubit_t> measuredQubits;
    std::unordered_set<Types::qubit_t> affectedQubits;
 
    bool executionStopped = false;
 
    for (size_t i = 0; i < operations.size(); ++i) {
      auto &op = operations[i];
      const auto qubits = op->AffectedQubits();
 
      bool executed = false;
 
      if (op->GetType() == OperationType::kMeasurement ||
          op->GetType() == OperationType::kConditionalMeasurement ||
          op->GetType() == OperationType::kRandomGen ||
          op->GetType() == OperationType::kConditionalRandomGen)
        measuredQubits.insert(qubits.begin(), qubits.end());
      else if (op->GetType() == OperationType::kReset) {
        // if it's the first op on qubit(s), execute it, otherwise treat it as a
        // measurement
        executed = true;
        for (auto qubit : qubits)
          if (affectedQubits.find(qubit) != affectedQubits.end()) {
            executed = false;
            break;
          }
 
        if (executed) {
          if (sim) {
            op->Execute(sim, state);
            if (curMaxBondDim) {
              const auto bondDim = sim->GetCurrentMaxBondDimension();
              if (bondDim > *curMaxBondDim) *curMaxBondDim = bondDim;
            }
          }
        } else
          measuredQubits.insert(qubits.begin(), qubits.end());
      } else  // regular gate or conditional gate
      {
        const auto bits = op->AffectedBits();
 
        // a measurement on a qubit prevents execution of any following gate
        // than affects the same qubit also a gate that's not executed and it
        // would affect certain qubits will prevent the execution of any
        // following gate that affects those qubits
 
        bool canExecute = op->GetType() == OperationType::kGate ||
                          op->GetType() == OperationType::kQuantumChannel ||
                          op->GetType() == OperationType::kDelay;
 
        if (canExecute)  // a conditional gate cannot be executed, it needs
                         // something executed at each shot, either a
                         // measurement or a random number generated
        {
          for (auto bit : bits)
            if (measuredQubits.find(bit) != measuredQubits.end()) {
              canExecute = false;
              break;
            }
 
          for (auto qubit : qubits)
            if (measuredQubits.find(qubit) != measuredQubits.end()) {
              canExecute = false;
              break;
            }
        }
 
        if (canExecute) {
          if (sim) {
            op->Execute(sim, state);
            if (curMaxBondDim) {
              const auto bondDim = sim->GetCurrentMaxBondDimension();
              if (bondDim > *curMaxBondDim) *curMaxBondDim = bondDim;
            }
          }
          executed = true;
        } else {
          // this is a 'trick', if it cannot execute, then neither can any
          // following gate that affects any of the already involved qubits
          measuredQubits.insert(bits.begin(), bits.end());
          measuredQubits.insert(qubits.begin(), qubits.end());
        }
      }
 
      affectedQubits.insert(qubits.begin(), qubits.end());
 
      if (!executed) {
        executionStopped = true;
        if (sim && sim->GetSimulationType() ==
                       Simulators::SimulationType::kMatrixProductState && sim->SupportsMPSSwapOptimization() && 
            op->GetType() != OperationType::kRandomGen &&
            op->GetType() != OperationType::kConditionalRandomGen && op->GetType() != OperationType::kNoOp)
          sim->SetGatesCounter(sim->GetGatesCounter() + 1);
      }
      if (executionStopped) executedOps.emplace_back(executed);
    }
 
    // if (sim) sim->Flush();
 
    return executedOps;
  }

  void ExecuteMeasurements(const std::shared_ptr<Simulators::ISimulator> &sim,
                           OperationState &state,
                           const std::vector<bool> &executedOps, size_t* curMaxBondDim = nullptr) const {
    state.Reset();
    if (!sim) return;
 
    // if (executedOps.empty() && !operations.empty()) throw
    // std::runtime_error("The executed operations vector is empty");
 
    const size_t dif = operations.size() - executedOps.size();
 
    for (size_t i = dif; i < operations.size(); ++i)
      if (!executedOps[i - dif]) {
        operations[i]->Execute(sim, state);
        if (curMaxBondDim) {
          const auto bondDim = sim->GetCurrentMaxBondDimension();
          if (bondDim > *curMaxBondDim) *curMaxBondDim = bondDim;
        }
      }
 
    // sim->Flush();
  }
 

  std::shared_ptr<Circuits::Circuit<Time>> RemoveExecutedOperations(
      std::vector<bool> &executedOps) const {
    if (executedOps.empty()) {
      executedOps.resize(size(), false);
      return std::make_shared<Circuit<Time>>(operations);
    }
 
    OperationsVector newops;
    newops.reserve(operations.size());
 
    const size_t dif = operations.size() - executedOps.size();
    for (size_t i = dif; i < operations.size(); ++i)
      if (!executedOps[i - dif]) newops.emplace_back(operations[i]);
 
    std::vector<bool> newExecutedOps(newops.size(), false);
    executedOps.swap(newExecutedOps);
 
    return std::make_shared<Circuit<Time>>(newops);
  }
 
 
 
  // used internally to optimize measurements in the case of having measurements
  // only at the end of the circuit
  std::shared_ptr<MeasurementOperation<Time>> GetLastMeasurements(
      const std::vector<bool> &executedOps, bool sort = true) const {
    const size_t dif = operations.size() - executedOps.size();
    std::vector<std::pair<Types::qubit_t, size_t>> measurements;
    measurements.reserve(dif);
 
    for (size_t i = dif; i < operations.size(); ++i)
      if (!executedOps[i - dif] &&
          operations[i]->GetType() == OperationType::kMeasurement) {
        auto measOp =
            std::static_pointer_cast<MeasurementOperation<Time>>(operations[i]);
        const auto &qubits = measOp->GetQubits();
        const auto &bits = measOp->GetBitsIndices();
 
        for (size_t j = 0; j < qubits.size(); ++j)
          measurements.emplace_back(qubits[j], bits[j]);
      }
 
    // qiskit aer expects sometimes to have them in sorted order, so...
    if (sort)
      std::sort(
          measurements.begin(), measurements.end(),
          [](const auto &p1, const auto &p2) { return p1.first < p2.first; });
 
    return std::make_shared<MeasurementOperation<Time>>(measurements);
  }

  bool HasConditionalOperations() const {
    for (const auto &op : operations)
      if (op->GetType() == OperationType::kConditionalGate ||
          op->GetType() == OperationType::kConditionalMeasurement ||
          op->GetType() == OperationType::kConditionalRandomGen)
        return true;
 
    return false;
  }

  bool ActsOnlyOnAdjacentQubits() const {
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
      if (qubits.size() <= 1) continue;
 
      if (qubits.size() == 2) {
        if (std::abs(qubits[0] - qubits[1]) != 1) return false;
      } else {
        Types::qubit_t minQubit = qubits[0];
        Types::qubit_t maxQubit = qubits[0];
 
        for (size_t i = 1; i < qubits.size(); ++i) {
          if (qubits[i] < minQubit)
            minQubit = qubits[i];
          else if (qubits[i] > maxQubit)
            maxQubit = qubits[i];
        }
 
        if (maxQubit - minQubit >= qubits.size()) return false;
      }
    }
 
    return true;
  }

  bool IsForest() const {
    std::unordered_map<Types::qubit_t, size_t> qubits;
    std::unordered_map<Types::qubit_t, Types::qubits_vector> lastQubits;
 
    for (const auto &op : operations) {
      const auto q = op->AffectedQubits();
      // one qubit gates or other operations that do not affect qubits do not
      // change anything
      if (q.size() <= 1) continue;
 
      bool allInTheLastQubits = true;
 
      for (const auto qubit : q) {
        if (lastQubits.find(qubit) == lastQubits.end()) {
          allInTheLastQubits = false;
          break;
        } else {
          const auto &lastQ = lastQubits[qubit];
 
          for (const auto q1 : q)
            if (std::find(lastQ.cbegin(), lastQ.cend(), q1) == lastQ.cend()) {
              allInTheLastQubits = false;
              break;
            }
 
          if (!allInTheLastQubits) break;
        }
      }
 
      if (allInTheLastQubits) continue;
 
      for (const auto qubit : q) {
        if (qubits[qubit] > 1)  // if the qubit is affected again...
          return false;
 
        ++qubits[qubit];
 
        lastQubits[qubit] = q;
      }
    }
 
    return true;
  }

  bool IsClifford() const override {
    for (const auto &op : operations)
      if (!op->IsClifford()) return false;
 
    return true;
  }

  double CliffordPercentage() const {
    size_t cliffordOps = 0;
    for (const auto &op : operations)
      if (op->IsClifford()) ++cliffordOps;
 
    return static_cast<double>(cliffordOps) / operations.size();
  }

  std::unordered_set<Types::qubit_t> GetCliffordQubits() const {
    std::unordered_set<Types::qubit_t> cliffordQubits;
    std::unordered_set<Types::qubit_t> nonCliffordQubits;
 
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
      if (op->IsClifford()) {
        for (const auto q : qubits) cliffordQubits.insert(q);
      } else {
        for (const auto q : qubits) nonCliffordQubits.insert(q);
      }
    }
 
    for (const auto q : nonCliffordQubits) cliffordQubits.erase(q);
 
    return cliffordQubits;
  }

  std::vector<std::shared_ptr<Circuit<Time>>> SplitCircuit() const {
    std::vector<std::shared_ptr<Circuit<Time>>> circuits;
 
    // find how many disjoint circuits we have in this circuit
 
    std::unordered_map<Types::qubit_t, std::unordered_set<Types::qubit_t>>
        circuitsMap;
    auto allQubits = GetQubits();
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitCircuitMap;
 
    // start with a bunch of disjoint sets of qubits, containing each a single
    // qubit
 
    for (auto qubit : allQubits) {
      circuitsMap[qubit] = std::unordered_set<Types::qubit_t>{qubit};
      qubitCircuitMap[qubit] = qubit;
    }
 
    // then the gates will join them together into circuits
 
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
 
      if (qubits.empty()) continue;
 
      auto qubitIt = qubits.cbegin();
      auto firstQubit = *qubitIt;
      // where is the first qubit in the disjoint sets?
      auto firstQubitCircuit = qubitCircuitMap[firstQubit];
 
      ++qubitIt;
 
      for (; qubitIt != qubits.cend(); ++qubitIt) {
        auto qubit = *qubitIt;
 
        // where is the qubit in the disjoint sets?
        auto qubitCircuit = qubitCircuitMap[qubit];
 
        // join the circuits / qubits sets
 
        if (firstQubitCircuit != qubitCircuit) {
          // join the circuits
          circuitsMap[firstQubitCircuit].insert(
              circuitsMap[qubitCircuit].begin(),
              circuitsMap[qubitCircuit].end());
 
          // update the qubit to circuit map
          for (auto q : circuitsMap[qubitCircuit])
            qubitCircuitMap[q] = firstQubitCircuit;
 
          // remove the joined circuit
          circuitsMap.erase(qubitCircuit);
        }
      }
    }
 
    size_t circSize = 1ULL;
    // static cast added to prevent compiler error on MacOS
    circSize = std::max(circSize, static_cast<size_t>(circuitsMap.size()));
    circuits.resize(circSize);
 
    for (size_t i = 0; i < circuits.size(); ++i)
      circuits[i] = std::make_shared<Circuit<Time>>();
 
    std::unordered_map<Types::qubit_t, size_t> qubitsSetsToCircuit;
 
    size_t circuitNo = 0;
    for (const auto &[id, qubitSet] : circuitsMap) {
      qubitsSetsToCircuit[id] = circuitNo;
 
      ++circuitNo;
    }
 
    // now fill them up with the operations
 
    for (const auto &op : operations) {
      const auto qubits = op->AffectedQubits();
 
      if (qubits.empty()) {
        circuits[0]->AddOperation(op->Clone());
        continue;
      }
 
      const auto circ = qubitsSetsToCircuit[qubitCircuitMap[*qubits.cbegin()]];
 
      circuits[circ]->AddOperation(op->Clone());
    }
 
    return circuits;
  }

  std::vector<std::shared_ptr<Circuits::Circuit<Time>>> ToLayers() const {
    std::vector<std::shared_ptr<Circuits::Circuit<Time>>> layers;
    layers.emplace_back(std::make_shared<Circuits::Circuit<Time>>());
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
    std::unordered_map<size_t, size_t> classicalBitLayer;
 
    for (const auto &op : GetOperations()) {
      // check the instruction, see if a new layer is needed
 
      // only qubits matter here, the others can be classicaly sent if they are
      // needed, even if they are shared... but that should be handled
      // somewhere, when not done implicitely (as for the 'simple network' case)
      if (op->CanAffectQuantumState()) {
        const auto qubits = op->AffectedQubits();
        size_t maxLevel = 0;
 
        for (Types::qubit_t qbit : qubits) {
          ++qubitsUsed[qbit];
          maxLevel = std::max(maxLevel, static_cast<size_t>(qubitsUsed[qbit]));
 
          if (layers.size() < qubitsUsed[qbit]) {
            auto circ = std::make_shared<Circuits::Circuit<Time>>();
            layers.push_back(std::move(circ));
          }
        }
 
        if (op->IsConditional()) {
          const auto bits = op->AffectedBits();
          for (const auto bit : bits)
            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
        }
 
        // now set all the qubits in the instruction to the max level
        for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
 
        // ensure enough layers exist for the computed level
        while (layers.size() <= layerIdx)
          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
 
        layers[layerIdx]->AddOperation(op->Clone());
 
        const auto writtenBits = op->AffectedBits();
        if (!writtenBits.empty() && !op->IsConditional()) {
          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
          for (const auto bit : writtenBits)
            classicalBitLayer[bit] =
                std::max(classicalBitLayer[bit], writtenLevel);
        }
      } else
        // add the instruction to the last layer
        layers.back()->AddOperation(op->Clone());
    }
 
    return layers;
  }

  std::vector<std::shared_ptr<Circuits::Circuit<Time>>> ToLayersNoClone()
      const {
    std::vector<std::shared_ptr<Circuits::Circuit<Time>>> layers;
    layers.emplace_back(std::make_shared<Circuits::Circuit<Time>>());
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
    std::unordered_map<size_t, size_t> classicalBitLayer;
 
    for (const auto &op : GetOperations()) {
      // check the instruction, see if a new layer is needed
 
      // only qubits matter here, the others can be classicaly sent if they are
      // needed, even if they are shared... but that should be handled
      // somewhere, when not done implicitely (as for the 'simple network' case)
      if (op->CanAffectQuantumState()) {
        const auto qubits = op->AffectedQubits();
        size_t maxLevel = 0;
 
        for (Types::qubit_t qbit : qubits) {
          ++qubitsUsed[qbit];
          maxLevel = std::max(maxLevel, static_cast<size_t>(qubitsUsed[qbit]));
 
          if (layers.size() < qubitsUsed[qbit]) {
            auto circ = std::make_shared<Circuits::Circuit<Time>>();
            layers.push_back(std::move(circ));
          }
        }
 
        if (op->IsConditional()) {
          const auto bits = op->AffectedBits();
          for (const auto bit : bits)
            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
        }
 
        // now set all the qubits in the instruction to the max level
        for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
 
        // ensure enough layers exist for the computed level
        while (layers.size() <= layerIdx)
          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
 
        layers[layerIdx]->AddOperation(op);
 
        const auto writtenBits = op->AffectedBits();
        if (!writtenBits.empty() && !op->IsConditional()) {
          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
          for (const auto bit : writtenBits)
            classicalBitLayer[bit] =
                std::max(classicalBitLayer[bit], writtenLevel);
        }
      } else
        // add the instruction to the last layer
        layers.back()->AddOperation(op);
    }
 
    return layers;
  }

  std::vector<std::shared_ptr<Circuits::Circuit<Time>>> ToMultipleQubitsLayers()
      const {
    std::vector<std::shared_ptr<Circuits::Circuit<Time>>> layers;
    layers.emplace_back(std::make_shared<Circuits::Circuit<Time>>());
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
    std::unordered_map<size_t, size_t> classicalBitLayer;
 
    for (const auto &op : GetOperations()) {
      // check the instruction, see if a new layer is needed
 
      // only qubits matter here, the others can be classicaly sent if they are
      // needed, even if they are shared... but that should be handled
      // somewhere, when not done implicitely (as for the 'simple network' case)
      if (op->CanAffectQuantumState()) {
        const auto qubits = op->AffectedQubits();
        size_t maxLevel = 0;
 
        for (Types::qubit_t qbit : qubits) {
          if (qubits.size() > 1) ++qubitsUsed[qbit];
          maxLevel = std::max(maxLevel, static_cast<size_t>(qubitsUsed[qbit]));
 
          if (layers.size() < qubitsUsed[qbit]) {
            auto circ = std::make_shared<Circuits::Circuit<Time>>();
            layers.push_back(std::move(circ));
          }
        }
 
        if (op->IsConditional()) {
          const auto bits = op->AffectedBits();
          for (const auto bit : bits)
            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
        }
 
        // now set all the qubits in the instruction to the max level
        for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
 
        // ensure enough layers exist for the computed level
        while (layers.size() <= layerIdx)
          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
 
        layers[layerIdx]->AddOperation(op->Clone());
 
        const auto writtenBits = op->AffectedBits();
        if (!writtenBits.empty() && !op->IsConditional()) {
          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
          for (const auto bit : writtenBits)
            classicalBitLayer[bit] =
                std::max(classicalBitLayer[bit], writtenLevel);
        }
      } else
        // add the instruction to the last layer
        layers.back()->AddOperation(op->Clone());
    }
 
    return layers;
  }

  std::vector<std::shared_ptr<Circuits::Circuit<Time>>>
  ToMultipleQubitsLayersNoClone() const {
    std::vector<std::shared_ptr<Circuits::Circuit<Time>>> layers;
    layers.emplace_back(std::make_shared<Circuits::Circuit<Time>>());
 
    std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
 
    // Track which layer (1-based, like qubitsUsed) each classical bit was last
    // written to, so that conditional operations reading those bits are placed
    // in the same layer or later.
    std::unordered_map<size_t, size_t> classicalBitLayer;
 
    for (const auto &op : GetOperations()) {
      // check the instruction, see if a new layer is needed
 
      // only qubits matter here, the others can be classicaly sent if they are
      // needed, even if they are shared... but that should be handled
      // somewhere, when not done implicitely (as for the 'simple network' case)
      if (op->CanAffectQuantumState()) {
        const auto qubits = op->AffectedQubits();
        size_t maxLevel = 0;
 
        for (Types::qubit_t qbit : qubits) {
          if (qubits.size() > 1) ++qubitsUsed[qbit];
          maxLevel = std::max(maxLevel, static_cast<size_t>(qubitsUsed[qbit]));
 
          if (layers.size() < qubitsUsed[qbit]) {
            auto circ = std::make_shared<Circuits::Circuit<Time>>();
            layers.push_back(std::move(circ));
          }
        }
 
        // For conditional operations, ensure this op is placed no earlier than
        // the layer where the classical bits it depends on were written.
        if (op->IsConditional()) {
          const auto bits = op->AffectedBits();
          for (const auto bit : bits)
            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
        }
 
        // now set all the qubits in the instruction to the max level
        for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
 
        // ensure enough layers exist for the computed level
        while (layers.size() <= layerIdx)
          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
 
        layers[layerIdx]->AddOperation(op);
 
        // Record the layer for classical bits written by measurements / resets
        // so that later conditional ops respect the dependency.
        const auto writtenBits = op->AffectedBits();
        if (!writtenBits.empty() && !op->IsConditional()) {
          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
          for (const auto bit : writtenBits)
            classicalBitLayer[bit] =
                std::max(classicalBitLayer[bit], writtenLevel);
        }
      } else
        // add the instruction to the last layer
        layers.back()->AddOperation(op);
    }
 
    return layers;
  }

  static std::shared_ptr<Circuits::Circuit<Time>> LayersToCircuit(
      const std::vector<std::shared_ptr<Circuits::Circuit<Time>>> &layers) {
    auto circuit{std::make_shared<Circuits::Circuit<Time>>()};
 
    for (const auto &layer : layers)
      circuit->AddOperations(layer->GetOperations());
 
    return circuit;
  }

  bool IsBranching() const override {
    for (const auto &op : GetOperations())
      if (op->IsBranching()) return true;
 
    return false;
  }

  iterator begin() noexcept { return operations.begin(); }

  iterator end() noexcept { return operations.end(); }
 

  const_iterator begin() const noexcept { return operations.begin(); }

  const_iterator end() const noexcept { return operations.end(); }

  const_iterator cbegin() const noexcept { return operations.cbegin(); }

  const_iterator cend() const noexcept { return operations.cend(); }

  reverse_iterator rbegin() noexcept { return operations.rbegin(); }

  reverse_iterator rend() noexcept { return operations.rend(); }

  const_reverse_iterator crbegin() const noexcept {
    return operations.crbegin();
  }

  const_reverse_iterator crend() const noexcept { return operations.crend(); }

  auto size() const { return operations.size(); }

  auto empty() const { return operations.empty(); }

  auto &operator[](size_t pos) { return operations[pos]; }

  const auto &operator[](size_t pos) const { return operations[pos]; }

  void resize(size_t size) {
    if (size < operations.size()) operations.resize(size);
  }
 
 private:
  void ReplaceThreeQubitAndSwapGates(bool onlyThreeQubits = false) {
    // just replace all three qubit gates(just ccnot and cswap will exist in the
    // first phase) with several gates on less qubits also replace swap gates
    // with three cnots (in the first phase) the controlled ones must be
    // replaced as well
 
    // TODO: if composite operations will be implemented, those need to be
    // optimized as well
 
    std::vector<std::shared_ptr<IOperation<Time>>> newops;
    newops.reserve(operations.size());
 
    for (std::shared_ptr<IOperation<Time>> op : operations) {
      if (op->GetType() == OperationType::kGate) {
        std::shared_ptr<IQuantumGate<Time>> gate =
            std::static_pointer_cast<IQuantumGate<Time>>(op);
 
        if (NeedsConversion(gate, onlyThreeQubits)) {
          std::vector<std::shared_ptr<IGateOperation<Time>>> newgates =
              ConvertGate(gate, onlyThreeQubits);
          newops.insert(newops.end(), newgates.begin(), newgates.end());
        } else
          newops.push_back(op);
      } else if (op->GetType() == OperationType::kConditionalGate) {
        std::shared_ptr<ConditionalGate<Time>> condgate =
            std::static_pointer_cast<ConditionalGate<Time>>(op);
        std::shared_ptr<IQuantumGate<Time>> gate =
            std::static_pointer_cast<IQuantumGate<Time>>(
                condgate->GetOperation());
 
        if (NeedsConversion(gate, onlyThreeQubits)) {
          std::vector<std::shared_ptr<IGateOperation<Time>>> newgates =
              ConvertGate(gate, onlyThreeQubits);
          std::shared_ptr<ICondition> cond = condgate->GetCondition();
 
          for (auto gate : newgates)
            newops.push_back(
                std::make_shared<ConditionalGate<Time>>(gate, cond));
        } else
          newops.push_back(op);
      } else
        newops.push_back(op);
    }
 
    operations.swap(newops);
  }

  static bool NeedsConversion(const std::shared_ptr<IQuantumGate<Time>> &gate,
                              bool onlyThreeQubits = false) {
    const bool hasThreeQubits = gate->GetNumQubits() == 3;
    if (onlyThreeQubits) return hasThreeQubits;
 
    return hasThreeQubits ||
           gate->GetGateType() == QuantumGateType::kSwapGateType;
  }

  static std::vector<std::shared_ptr<IGateOperation<Time>>> ConvertGate(
      std::shared_ptr<IQuantumGate<Time>> &gate, bool onlyThreeQubits = false) {
    // TODO: if delays are used, how to transfer delays from the converted gate
    // to the resulting gates?
    std::vector<std::shared_ptr<IGateOperation<Time>>> newops;
 
    if (gate->GetNumQubits() == 3) {
      // must be converted no matter what
      if (gate->GetGateType() == QuantumGateType::kCCXGateType) {
        const size_t q1 = gate->GetQubit(0);  // control 1
        const size_t q2 = gate->GetQubit(1);  // control 2
        const size_t q3 = gate->GetQubit(2);  // target
 
        // Sleator-Weinfurter decomposition
        newops.push_back(std::make_shared<CSxGate<Time>>(q2, q3));
        newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
        newops.push_back(std::make_shared<CSxDagGate<Time>>(q2, q3));
        newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
        newops.push_back(std::make_shared<CSxGate<Time>>(q1, q3));
      } else if (gate->GetGateType() == QuantumGateType::kCSwapGateType) {
        const size_t q1 = gate->GetQubit(0);  // control 1
        const size_t q2 = gate->GetQubit(1);  // control 2
        const size_t q3 = gate->GetQubit(2);  // target
 
        // TODO: find a better decomposition
        // this one I've got with the qiskit transpiler
        newops.push_back(std::make_shared<CXGate<Time>>(q3, q2));
 
        newops.push_back(std::make_shared<CSxGate<Time>>(q2, q3));
        newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
        newops.push_back(std::make_shared<PhaseGate<Time>>(q3, M_PI));
 
        newops.push_back(std::make_shared<PhaseGate<Time>>(q2, -M_PI_2));
 
        newops.push_back(std::make_shared<CSxGate<Time>>(q2, q3));
        newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
        newops.push_back(std::make_shared<PhaseGate<Time>>(q3, M_PI));
 
        newops.push_back(std::make_shared<CSxGate<Time>>(q1, q3));
 
        newops.push_back(std::make_shared<CXGate<Time>>(q3, q2));
      } else
        newops.push_back(gate);
    } else if (!onlyThreeQubits &&
               gate->GetGateType() == QuantumGateType::kSwapGateType) {
      // must be converted no matter what
      const size_t q1 = gate->GetQubit(0);
      const size_t q2 = gate->GetQubit(1);
 
      // for now replace it with three cnots, but maybe later make it
      // configurable there are other possibilities, for example three cy gates
      newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
      newops.push_back(std::make_shared<CXGate<Time>>(q2, q1));
      newops.push_back(std::make_shared<CXGate<Time>>(q1, q2));
    } else
      newops.push_back(gate);
 
    return newops;
  }
 
  OperationsVector operations; 
};

template <typename Time = Types::time_type>
class ComparableCircuit : public Circuit<Time> {
 public:
  using BaseClass = Circuit<Time>;                 
  using Operation = IOperation<Time>;              
  using OperationPtr = std::shared_ptr<Operation>; 
  using OperationsVector =
      std::vector<OperationPtr>; 

  ComparableCircuit(const OperationsVector &ops = {}) : BaseClass(ops) {}

  ComparableCircuit(const BaseClass &circ) : BaseClass(circ.GetOperations()) {}

  ComparableCircuit &operator=(const BaseClass &circ) {
    BaseClass::SetOperations(circ.GetOperations());
 
    return *this;
  }

  bool operator==(const BaseClass &rhs) const {
    if (BaseClass::GetOperations().size() != rhs.GetOperations().size())
      return false;
 
    for (size_t i = 0; i < BaseClass::GetOperations().size(); ++i) {
      if (BaseClass::GetOperations()[i]->GetType() !=
          rhs.GetOperations()[i]->GetType())
        return false;
 
      switch (BaseClass::GetOperations()[i]->GetType()) {
        case OperationType::kGate:
          if (std::static_pointer_cast<IQuantumGate<Time>>(
                  BaseClass::GetOperations()[i])
                      ->GetGateType() !=
                  std::static_pointer_cast<IQuantumGate<Time>>(
                      rhs.GetOperations()[i])
                      ->GetGateType() ||
              BaseClass::GetOperations()[i]->AffectedBits() !=
                  rhs.GetOperations()[i]->AffectedBits())
            return false;
          if (approximateParamsCheck) {
            const auto params1 = std::static_pointer_cast<IQuantumGate<Time>>(
                                     BaseClass::GetOperations()[i])
                                     ->GetParams();
            const auto params2 = std::static_pointer_cast<IQuantumGate<Time>>(
                                     rhs.GetOperations()[i])
                                     ->GetParams();
            if (params1.size() != params2.size()) return false;
 
            for (size_t j = 0; j < params1.size(); ++j)
              if (std::abs(params1[j] - params2[j]) > paramsEpsilon)
                return false;
          } else if (std::static_pointer_cast<IQuantumGate<Time>>(
                         BaseClass::GetOperations()[i])
                         ->GetParams() !=
                     std::static_pointer_cast<IQuantumGate<Time>>(
                         rhs.GetOperations()[i])
                         ->GetParams())
            return false;
          break;
        case OperationType::kMeasurement:
          if (BaseClass::GetOperations()[i]->AffectedQubits() !=
                  rhs.GetOperations()[i]->AffectedQubits() ||
              BaseClass::GetOperations()[i]->AffectedBits() !=
                  rhs.GetOperations()[i]->AffectedBits())
            return false;
          break;
        case OperationType::kRandomGen:
          if (BaseClass::GetOperations()[i]->AffectedBits() !=
              rhs.GetOperations()[i]->AffectedBits())
            return false;
          break;
        case OperationType::kConditionalGate:
        case OperationType::kConditionalMeasurement:
        case OperationType::kConditionalRandomGen: {
          // first, check the conditions
          const auto leftCondition =
              std::static_pointer_cast<IConditionalOperation<Time>>(
                  BaseClass::GetOperations()[i])
                  ->GetCondition();
          const auto rightCondition =
              std::static_pointer_cast<IConditionalOperation<Time>>(
                  rhs.GetOperations()[i])
                  ->GetCondition();
          if (leftCondition->GetBitsIndices() !=
              rightCondition->GetBitsIndices())
            return false;
 
          const auto leftEqCondition =
              std::static_pointer_cast<EqualCondition>(leftCondition);
          const auto rightEqCondition =
              std::static_pointer_cast<EqualCondition>(rightCondition);
          if (!leftEqCondition || !rightEqCondition) return false;
 
          if (leftEqCondition->GetAllBits() != rightEqCondition->GetAllBits())
            return false;
 
          // now check the operations
          const auto leftOp =
              std::static_pointer_cast<IConditionalOperation<Time>>(
                  BaseClass::GetOperations()[i])
                  ->GetOperation();
          const auto rightOp =
              std::static_pointer_cast<IConditionalOperation<Time>>(
                  rhs.GetOperations()[i])
                  ->GetOperation();
 
          ComparableCircuit<Time> leftCircuit;
          BaseClass rightCircuit;
          leftCircuit.SetApproximateParamsCheck(approximateParamsCheck);
          leftCircuit.AddOperation(leftOp);
          rightCircuit.AddOperation(rightOp);
 
          if (leftCircuit != rightCircuit) return false;
        } break;
        case OperationType::kReset:
          if (BaseClass::GetOperations()[i]->AffectedQubits() !=
                  rhs.GetOperations()[i]->AffectedQubits() ||
              std::static_pointer_cast<Reset<Time>>(
                  BaseClass::GetOperations()[i])
                      ->GetResetTargets() !=
                  std::static_pointer_cast<Reset<Time>>(rhs.GetOperations()[i])
                      ->GetResetTargets())
            return false;
          break;
        case OperationType::kQuantumChannel: {
          const auto left =
              std::static_pointer_cast<QuantumChannelOperation<Time>>(
                  BaseClass::GetOperations()[i]);
          const auto right =
              std::static_pointer_cast<QuantumChannelOperation<Time>>(
                  rhs.GetOperations()[i]);
          if (left->AffectedQubits() != right->AffectedQubits() ||
              !left->GetChannel().IsApprox(
                  right->GetChannel(),
                  approximateParamsCheck ? paramsEpsilon : 0.0))
            return false;
        } break;
        case OperationType::kDelay: {
          const auto left =
              std::static_pointer_cast<Delay<Time>>(
                  BaseClass::GetOperations()[i]);
          const auto right =
              std::static_pointer_cast<Delay<Time>>(
                  rhs.GetOperations()[i]);
          if (left->GetQubit() != right->GetQubit()) return false;
          if (approximateParamsCheck) {
            if (std::abs(left->GetDuration() - right->GetDuration()) > paramsEpsilon)
              return false;
          } else if (left->GetDuration() != right->GetDuration()) {
            return false;
          }
        } break;
        case OperationType::kNoOp:
          break;
        default:
          return false;
      }
 
      if (BaseClass::GetOperations()[i]->GetDelay() !=
          rhs.GetOperations()[i]->GetDelay())
        return false;
    }
 
    return true;
  }

  bool operator!=(const BaseClass &rhs) const { return !(*this == rhs); }

  void SetApproximateParamsCheck(bool check) { approximateParamsCheck = check; }

  bool GetApproximateParamsCheck() const { return approximateParamsCheck; }

  void SetParamsEpsilon(double eps) { paramsEpsilon = eps; }

  double GetParamsEpsilon() const { return paramsEpsilon; }
 
 private:
  bool approximateParamsCheck =
      false; 
  double paramsEpsilon = 1e-8; 
};
 
}  // namespace Circuits
 
#endif  // !_CIRCUIT_H_