Circuit.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _CIRCUIT_H_
#define _CIRCUIT_H_
 
#define _USE_MATH_DEFINES
#include <algorithm>
#include <math.h>
#include <numeric>
#include <set>
 
#include "Operations/Conditional.h"
#include "Operations/Delay.h"
#include "Operations/Operations.h"
#include "Operations/QuantumChannelOperation.h"
#include "Operations/QuantumGates.h"
#include "Operations/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)
            ExecuteOperation(op, sim, state, curMaxBondDim);
        if (sim->IsGateFusionEnabled())
            sim->FlushPendingGates();
        if (curMaxBondDim)
            *curMaxBondDim = std::max(*curMaxBondDim, sim->GetExecutedMaxBondDimension());
    }

    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);
    }

    bool HasUnseededRandomOperations() const
    {
        for (const auto &op : operations)
        {
            if (op->GetType() == OperationType::kRandomGen && !std::static_pointer_cast<Random<Time>>(op)->HasExplicitSeed())
                return true;
            if (op->GetType() == OperationType::kConditionalRandomGen &&
                !std::static_pointer_cast<Random<Time>>(std::static_pointer_cast<ConditionalRandomGen<Time>>(op)->GetOperation())->HasExplicitSeed())
                return true;
            if (op->GetType() == OperationType::kComposite && std::static_pointer_cast<Circuit<Time>>(op)->HasUnseededRandomOperations())
                return true;
        }
        return false;
    }

    std::shared_ptr<Circuit<Time>> CloneForExecution(uint64_t stream, uint64_t defaultSeed = 0) const
    {
        OperationsVector newops;
        newops.reserve(operations.size());
        for (size_t i = 0; i < operations.size(); ++i)
        {
            const auto &op = operations[i];
            switch (op->GetType())
            {
            case OperationType::kRandomGen:
                newops.push_back(std::static_pointer_cast<Random<Time>>(op)->CloneForExecution(Simulators::IState::DeriveSeed(stream, i), defaultSeed));
                break;
            case OperationType::kConditionalRandomGen: {
                auto copy = std::static_pointer_cast<ConditionalRandomGen<Time>>(op->Clone());
                copy->SetOperation(
                    std::static_pointer_cast<Random<Time>>(copy->GetOperation())->CloneForExecution(Simulators::IState::DeriveSeed(stream, i), defaultSeed));
                newops.push_back(std::move(copy));
                break;
            }
            case OperationType::kComposite:
                newops.push_back(std::static_pointer_cast<Circuit<Time>>(op)->CloneForExecution(Simulators::IState::DeriveSeed(stream, i), defaultSeed));
                break;
            default:
                newops.push_back(op);
                break;
            }
        }
        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);
    }

    std::shared_ptr<Circuit<Time>> RemoveEndMeasurements() const
    {
        OperationsVector retained;
        retained.reserve(operations.size());
 
        std::unordered_set<Types::qubit_t> activeQubits;
        std::unordered_set<size_t> neededBits;
 
        for (auto it = operations.rbegin(); it != operations.rend(); ++it)
        {
            const auto &op = *it;
            if (!op)
                continue;
 
            const auto type = op->GetType();
            if (type == OperationType::kMeasurement || type == OperationType::kConditionalMeasurement)
            {
                std::shared_ptr<MeasurementOperation<Time>> meas;
                std::shared_ptr<ConditionalMeasurement<Time>> condMeas;
 
                if (type == OperationType::kMeasurement)
                {
                    meas = std::static_pointer_cast<MeasurementOperation<Time>>(op);
                }
                else
                {
                    condMeas = std::static_pointer_cast<ConditionalMeasurement<Time>>(op);
                    meas = std::static_pointer_cast<MeasurementOperation<Time>>(condMeas->GetOperation());
                }
 
                if (!meas)
                    continue;
 
                const auto &qubits = meas->GetQubits();
                const auto &bits = meas->GetBitsIndices();
                const size_t count = std::min(qubits.size(), bits.size());
                if (count == 0)
                    continue;
 
                std::vector<bool> keepFlags(count, false);
                size_t keptCount = 0;
 
                for (size_t idx = count; idx-- > 0;)
                {
                    const auto q = qubits[idx];
                    const auto b = bits[idx];
                    if (activeQubits.count(q) > 0 || neededBits.count(b) > 0)
                    {
                        keepFlags[idx] = true;
                        ++keptCount;
                        if (type == OperationType::kMeasurement)
                            neededBits.erase(b);
                    }
                }
 
                if (keptCount == 0)
                    continue;
 
                if (type == OperationType::kConditionalMeasurement)
                {
                    const auto condBits = condMeas->AffectedBits();
                    neededBits.insert(condBits.begin(), condBits.end());
                }
 
                if (keptCount == count)
                {
                    retained.push_back(op);
                }
                else
                {
                    std::vector<std::pair<Types::qubit_t, size_t>> keptPairs;
                    keptPairs.reserve(keptCount);
                    for (size_t i = 0; i < count; ++i)
                    {
                        if (keepFlags[i])
                            keptPairs.emplace_back(qubits[i], bits[i]);
                    }
 
                    auto filteredMeas = std::make_shared<MeasurementOperation<Time>>(keptPairs, meas->GetDelay());
                    if (meas->HasReadout())
                    {
                        const auto &rates = meas->GetReadout();
                        std::vector<ReadoutRates> keptRates;
                        keptRates.reserve(keptCount);
                        for (size_t i = 0; i < count; ++i)
                        {
                            if (keepFlags[i])
                                keptRates.push_back(i < rates.size() ? rates[i] : ReadoutRates{});
                        }
                        filteredMeas->SetReadout(std::move(keptRates));
                    }
 
                    if (type == OperationType::kConditionalMeasurement)
                    {
                        retained.push_back(std::make_shared<ConditionalMeasurement<Time>>(filteredMeas, condMeas->GetCondition(), condMeas->GetDelay()));
                    }
                    else
                    {
                        retained.push_back(std::move(filteredMeas));
                    }
                }
            }
            else
            {
                retained.push_back(op);
 
                switch (type)
                {
                case OperationType::kGate:
                case OperationType::kReset:
                case OperationType::kQuantumChannel: {
                    const auto qs = op->AffectedQubits();
                    activeQubits.insert(qs.begin(), qs.end());
                    break;
                }
                case OperationType::kConditionalGate: {
                    const auto qs = op->AffectedQubits();
                    activeQubits.insert(qs.begin(), qs.end());
                    const auto bs = op->AffectedBits();
                    neededBits.insert(bs.begin(), bs.end());
                    break;
                }
                case OperationType::kRandomGen: {
                    const auto bs = op->AffectedBits();
                    for (const auto b : bs)
                        neededBits.erase(b);
                    break;
                }
                case OperationType::kConditionalRandomGen: {
                    const auto bs = op->AffectedBits();
                    neededBits.insert(bs.begin(), bs.end());
                    break;
                }
                case OperationType::kComposite: {
                    const auto qs = op->AffectedQubits();
                    activeQubits.insert(qs.begin(), qs.end());
                    const auto bs = GetReferencedBits(op);
                    neededBits.insert(bs.begin(), bs.end());
                    break;
                }
                case OperationType::kDelay:
                case OperationType::kNoOp:
                default:
                    break;
                }
            }
        }
 
        std::reverse(retained.begin(), retained.end());
        return std::make_shared<Circuit<Time>>(retained);
    }

    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;
 
        // Number qubits in ascending order, not in order of first use: chain
        // simulators (MPS/MPO) take the numbering as their site order, and a
        // reordered gate sequence would otherwise scramble it.
        std::vector<Types::qubit_t> usedQubits;
        for (const auto &op : operations)
            for (const auto qubit : op->AffectedQubits())
                usedQubits.push_back(qubit);
        std::sort(usedQubits.begin(), usedQubits.end());
        usedQubits.erase(std::unique(usedQubits.begin(), usedQubits.end()), usedQubits.end());
        for (const auto qubit : usedQubits)
            if (newQubitsMap.find(qubit) == newQubitsMap.end())
            {
                newQubitsMap[qubit] = nrQubits;
                ++nrQubits;
            }
 
        for (const auto &op : operations)
        {
            const auto affectedBits = GetReferencedBits(op);
 
            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;
            // Only for bits whose source measurement carried readout rates.
            std::unordered_map<size_t, ReadoutRates> measReadout;
 
            auto affectedBits = op->AffectedBits();
            auto affectedQubits = op->AffectedQubits();
 
            const auto firstMeas = std::static_pointer_cast<MeasurementOperation<Time>>(op);
            const auto &firstRates = firstMeas->GetReadout();
 
            for (size_t q = 0; q < affectedQubits.size(); ++q)
            {
                bits.insert(affectedBits[q]);
                measQubits[affectedBits[q]] = affectedQubits[q];
                measDelays[affectedBits[q]] = op->GetDelay();
                if (q < firstRates.size())
                    measReadout[affectedBits[q]] = firstRates[q];
            }
 
            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();
                const auto &rates = meas->GetReadout();
                for (size_t q = 0; q < affectedBits.size(); ++q)
                {
                    bits.insert(affectedBits[q]);
                    measQubits[affectedBits[q]] = affectedQubits[q];
                    measDelays[affectedBits[q]] = op2->GetDelay();
                    if (q < rates.size())
                        measReadout[affectedBits[q]] = rates[q];
                    else
                        measReadout.erase(affectedBits[q]);
                }
            }
 
            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())
                        {
                            auto rebuilt = std::make_shared<MeasurementOperation<Time>>(std::vector{std::make_pair(measQubits[bit], bit)}, measDelays[bit]);
                            const auto rateit = measReadout.find(bit);
                            if (rateit != measReadout.end())
                                rebuilt->SetReadout({rateit->second});
                            newops.emplace_back(rebuilt);
                            bits.erase(bit);
                        }
                }
                if (bits.empty())
                    break;
            }
 
            // now add the measurements that were left in any order
            for (auto bit : bits)
            {
                auto rebuilt = std::make_shared<MeasurementOperation<Time>>(std::vector{std::make_pair(measQubits[bit], bit)}, measDelays[bit]);
                const auto rateit = measReadout.find(bit);
                if (rateit != measReadout.end())
                    rebuilt->SetReadout({rateit->second});
                newops.emplace_back(rebuilt);
            }
        }
 
        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 = GetReferencedBits(op);
            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 = GetReferencedBits(op);
            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 = GetReferencedBits(op);
            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()
    {
        // The dependency planner does not inspect nested controls or measurements.
        // In particular, a classical-only composite has no qubit queue to drain.
        if (HasCompositeOperations())
            return;
        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 = GetReferencedBits(op);
            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 HasCompositeOperations() const
    {
        return std::any_of(operations.begin(), operations.end(), [](const auto &op) { return op->GetType() == OperationType::kComposite; });
    }

    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)
        {
            // Prefix execution and terminal-measurement extraction treat composites
            // as opaque operations. Run them per shot, including gate-only
            // composites, until those optimizations can inspect their contents
            // consistently.
            if (op->GetType() == OperationType::kComposite)
                return true;
            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);
                }
            }
        }
 
        // A trailing mid-circuit reset still has to run. ExecuteNonMeasurements
        // defers it like a measurement, and the terminal-measurement path would
        // drop it (e.g. sampled T1 noise before an expectation value).
        return !resetQubits.empty();
    }

    std::vector<bool> ExecuteNonMeasurements(const std::shared_ptr<Simulators::ISimulator> &sim, OperationState &state, size_t *curMaxBondDim = nullptr) const
    {
        // With gate fusion, the look-ahead routing of an MPS/MPO follows the
        // simulator's upcoming operations in order, and stepping over an operation
        // that is deferred to the shots breaks that plan for the rest of this
        // prefix. Plan the order that actually runs instead: the operations
        // executed here first, then the deferred ones.
        bool planExecutionOrder = sim && sim->IsGateFusionEnabled() && sim->IsRoutingLookaheadEnabled() && sim->GetGatesCounter() == 0;
        std::vector<bool> mask;
        if (planExecutionOrder)
        {
            OperationState dryState(state);
            mask = ExecuteNonMeasurements(nullptr, dryState);
            planExecutionOrder = std::find(mask.begin(), mask.end(), true) != mask.end() && std::find(mask.begin(), mask.end(), false) != mask.end();
        }
        if (planExecutionOrder)
        {
            const size_t dif = operations.size() - mask.size();
            OperationsVector ordered(operations.begin(), operations.begin() + dif);
            ordered.reserve(operations.size());
            for (size_t i = dif; i < operations.size(); ++i)
                if (mask[i - dif])
                    ordered.emplace_back(operations[i]);
            for (size_t i = dif; i < operations.size(); ++i)
                if (!mask[i - dif])
                    ordered.emplace_back(operations[i]);
            sim->SetUpcomingGates(ordered);
        }
 
        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)
                    {
                        ExecuteOperation(op, sim, state, curMaxBondDim);
                    }
                }
                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)
                    {
                        ExecuteOperation(op, sim, state, curMaxBondDim);
                    }
                    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->IsGateFusionEnabled())
                {
                    // with the execution order planned, deferred operations are at the
                    // end of the upcoming operations and are not stepped over here
                    if (!planExecutionOrder)
                        sim->IncrementGatesCounter();
                }
                else if (sim &&
                         (sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductState ||
                          sim->GetSimulationType() == Simulators::SimulationType::kMatrixProductOperator) &&
                         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)
        {
            if (sim->IsGateFusionEnabled())
                sim->FlushPendingGates();
            if (curMaxBondDim)
                *curMaxBondDim = std::max(*curMaxBondDim, sim->GetExecutedMaxBondDimension());
        }
 
        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])
                ExecuteOperation(operations[i], sim, state, curMaxBondDim);
 
        if (sim->IsGateFusionEnabled())
            sim->FlushPendingGates();
        if (curMaxBondDim)
            *curMaxBondDim = std::max(*curMaxBondDim, sim->GetExecutedMaxBondDimension());
    }

    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;
        std::vector<ReadoutRates> readout;
        std::unordered_map<size_t, size_t> lastWrite;
        bool anyReadout = false;
        measurements.reserve(dif);
        readout.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();
                const auto &rates = measOp->GetReadout();
 
                for (size_t j = 0; j < qubits.size(); ++j)
                {
                    lastWrite[bits[j]] = measurements.size();
                    measurements.emplace_back(qubits[j], bits[j]);
                    readout.push_back(j < rates.size() ? rates[j] : ReadoutRates{});
                }
 
                anyReadout = anyReadout || measOp->HasReadout();
            }
 
        // Only the final write to each classical bit is observable in a terminal
        // measurement batch. Drop overwritten writes before sorting for a backend,
        // otherwise sorting can change which qubit (and readout rate) wins.
        std::vector<size_t> order;
        order.reserve(lastWrite.size());
        for (size_t i = 0; i < measurements.size(); ++i)
            if (lastWrite.at(measurements[i].second) == i)
                order.push_back(i);
 
        // Keep each surviving readout rate attached to its measurement pair.
        if (sort)
            std::sort(order.begin(), order.end(), [&](size_t a, size_t b) { return measurements[a].first < measurements[b].first; });
 
        if (sort || order.size() != measurements.size())
        {
            std::vector<std::pair<Types::qubit_t, size_t>> sortedMeas;
            std::vector<ReadoutRates> sortedReadout;
            sortedMeas.reserve(order.size());
            sortedReadout.reserve(order.size());
            for (const auto idx : order)
            {
                sortedMeas.push_back(measurements[idx]);
                sortedReadout.push_back(readout[idx]);
            }
            measurements.swap(sortedMeas);
            readout.swap(sortedReadout);
        }
 
        auto merged = std::make_shared<MeasurementOperation<Time>>(measurements);
        if (anyReadout)
            merged->SetReadout(std::move(readout));
 
        return merged;
    }

    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;
        // The last layer reading each bit: a later write (a measurement) must not
        // move ahead of a read of the old value or of an earlier write.
        std::unordered_map<size_t, size_t> classicalBitReadLayer;
 
        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]);
                }
                else
                {
                    // the bits are written: keep the order of the writes and reads
                    for (const auto bit : op->AffectedBits())
                        maxLevel = std::max(maxLevel, std::max(classicalBitLayer[bit], classicalBitReadLayer[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 usedBits = op->AffectedBits();
                const size_t usedLevel = maxLevel > 0 ? maxLevel : 1;
                auto &bitLayers = op->IsConditional() ? classicalBitReadLayer : classicalBitLayer;
                for (const auto bit : usedBits)
                    bitLayers[bit] = std::max(bitLayers[bit], usedLevel);
            }
            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;
        // The last layer reading each bit: a later write (a measurement) must not
        // move ahead of a read of the old value or of an earlier write.
        std::unordered_map<size_t, size_t> classicalBitReadLayer;
 
        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]);
                }
                else
                {
                    // the bits are written: keep the order of the writes and reads
                    for (const auto bit : op->AffectedBits())
                        maxLevel = std::max(maxLevel, std::max(classicalBitLayer[bit], classicalBitReadLayer[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 usedBits = op->AffectedBits();
                const size_t usedLevel = maxLevel > 0 ? maxLevel : 1;
                auto &bitLayers = op->IsConditional() ? classicalBitReadLayer : classicalBitLayer;
                for (const auto bit : usedBits)
                    bitLayers[bit] = std::max(bitLayers[bit], usedLevel);
            }
            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;
        // The last layer reading each bit: a later write (a measurement) must not
        // move ahead of a read of the old value or of an earlier write.
        std::unordered_map<size_t, size_t> classicalBitReadLayer;
 
        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]);
                }
                else
                {
                    // the bits are written: keep the order of the writes and reads
                    for (const auto bit : op->AffectedBits())
                        maxLevel = std::max(maxLevel, std::max(classicalBitLayer[bit], classicalBitReadLayer[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 usedBits = op->AffectedBits();
                const size_t usedLevel = maxLevel > 0 ? maxLevel : 1;
                auto &bitLayers = op->IsConditional() ? classicalBitReadLayer : classicalBitLayer;
                for (const auto bit : usedBits)
                    bitLayers[bit] = std::max(bitLayers[bit], usedLevel);
            }
            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;
        // Likewise the last layer reading each bit: a later write (a measurement)
        // must not move ahead of a read of the old value or of an earlier write.
        std::unordered_map<size_t, size_t> classicalBitReadLayer;
 
        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]);
                }
                else
                {
                    // the bits are written: keep the order of the writes and reads
                    for (const auto bit : op->AffectedBits())
                        maxLevel = std::max(maxLevel, std::max(classicalBitLayer[bit], classicalBitReadLayer[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, and the layer
                // of the bits read by conditional ops so that later writes do.
                const auto usedBits = op->AffectedBits();
                const size_t usedLevel = maxLevel > 0 ? maxLevel : 1;
                auto &bitLayers = op->IsConditional() ? classicalBitReadLayer : classicalBitLayer;
                for (const auto bit : usedBits)
                    bitLayers[bit] = std::max(bitLayers[bit], usedLevel);
            }
            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:
    // A conditional's AffectedBits() is deliberately predicate-only. Register
    // sizing, mapping and dependencies also need the wrapped operation's outputs.
    // Nested circuits collect their referenced bits through GetBits().
    static std::vector<size_t> GetReferencedBits(const OperationPtr &op)
    {
        auto bits = op->AffectedBits();
        if (op->IsConditional())
        {
            const auto conditional = std::static_pointer_cast<IConditionalOperation<Time>>(op);
            if (const auto operation = conditional->GetOperation())
            {
                const auto outputs = GetReferencedBits(operation);
                bits.insert(bits.end(), outputs.begin(), outputs.end());
            }
        }
        return bits;
    }
 
    static void ExecuteOperation(const OperationPtr &op, const std::shared_ptr<Simulators::ISimulator> &sim, OperationState &state, size_t *curMaxBondDim)
    {
        const bool boundary = sim->IsGateFusionEnabled() && op->GetType() != OperationType::kGate;
        const auto sourcePosition = sim->GetGatesCounter();
        if (boundary)
            sim->FlushPendingGates();
        if (op->GetType() == OperationType::kComposite)
        {
            // Reset classical bits once at the start of the shot, not when entering
            // a nested circuit: its conditions and writes share the enclosing state.
            const auto circuit = std::static_pointer_cast<Circuit<Time>>(op);
            for (const auto &nested : circuit->GetOperations())
                ExecuteOperation(nested, sim, state, curMaxBondDim);
            if (boundary)
            {
                sim->FlushPendingGates();
                sim->SetGatesCounter(sourcePosition + 1);
            }
            return;
        }
        op->Execute(sim, state);
        // A circuit boundary is one source operation, even if it invokes no
        // simulator primitive (delay/classical op) or several (reset/conditional).
        if (boundary)
        {
            sim->FlushPendingGates();
            sim->SetGatesCounter(sourcePosition + 1);
        }
        if (curMaxBondDim)
        {
            const auto bondDim = sim->GetExecutedMaxBondDimension();
            if (bondDim > *curMaxBondDim)
                *curMaxBondDim = bondDim;
        }
    }

    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<Circuits::Delay<Time>>(BaseClass::GetOperations()[i]);
                const auto right = std::static_pointer_cast<Circuits::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_