ExecutionCost.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef __EXECUTION_COST_H_
#define __EXECUTION_COST_H_
 
#include "../Circuit/Factory.h"
#include "../Simulators/Core/Factory.h"
#include "../Simulators/QCSim/QcsimPauliPropagator.h"
#include "../Utils/Interpolation/MultipleLinearRegression.h"
#include "../Utils/Interpolation/MultivariateHermiteInterpolation.h"
#include "../Utils/LogFile.h"
 
#include <algorithm>
#include <cstddef>
#include <iomanip>
 
namespace Estimators
{
 
// this is not for time estimation, it can be used for something very, very,
// VERY rough, to decide just if it would execute in a reasonable time or not it
// cannot be used to compare runtime for two simulators of the same circuit, but
// it sort of can be used to compare two circuits in the same simulator
// (especially if the same number of samples is used for both circuits) for
// something better we have the estimators which got the constants from
// benchmarking and rely on implementation details of the simulators, so they
// are more accurate, but also less general
 
// Except for pauli propagation, where the operation cost depends on the
// position in the circuit (if the circuit is non-clifford), using a neural
// network (or even a simpler regressor) to learn the relation between the
// estimated cost and the actual execution time (single threaded, multithreading
// rises some other issues) would be probably possible
class ExecutionCost
{
  public:
    struct CircuitInfo
    {
        size_t nrQubits = 0;
        size_t nrOneQubitOps = 0;
        size_t nrTwoQubitOps = 0;
        size_t nrThreeQubitOps = 0;
        size_t nrMiddleMeasurementOps = 0;
        size_t nrEndMeasurementOps = 0;
        size_t nrOneQubitOpsExecutedOnce = 0;
        size_t nrTwoQubitOpsExecutedOnce = 0;
        size_t nrThreeQubitOpsExecutedOnce = 0;
 
        CircuitInfo() = default;
        CircuitInfo(const CircuitInfo &other) = default;
 
        CircuitInfo &operator=(const CircuitInfo &other) = default;
 
        double getFieldValue(size_t index) const
        {
            switch (index)
            {
            case 0:
                return static_cast<double>(nrQubits);
            case 1:
                return static_cast<double>(nrOneQubitOps);
            case 2:
                return static_cast<double>(nrTwoQubitOps);
            case 3:
                return static_cast<double>(nrThreeQubitOps);
            case 4:
                return static_cast<double>(nrMiddleMeasurementOps);
            case 5:
                return static_cast<double>(nrEndMeasurementOps);
            case 6:
                return static_cast<double>(nrOneQubitOpsExecutedOnce);
            case 7:
                return static_cast<double>(nrTwoQubitOpsExecutedOnce);
            case 8:
                return static_cast<double>(nrThreeQubitOpsExecutedOnce);
            default:
                throw std::out_of_range("Invalid index for CircuitInfo field");
            }
        }
    };
 
    struct ExecutionInfo : public CircuitInfo
    {
        size_t nrSamples = 0;
        size_t nrQubitsSampled = 0;
        size_t maxBondDim = 0;
        size_t nrPauliOps = 0;
        double executionCost = 0;
        double runtime = 0;
 
        ExecutionInfo() = default;
        ExecutionInfo(const ExecutionInfo &other) = default;
 
        ExecutionInfo(const CircuitInfo &circuitInfo) : CircuitInfo(circuitInfo)
        {
        }
 
        ExecutionInfo &operator=(const ExecutionInfo &other) = default;
 
        ExecutionInfo &operator=(const CircuitInfo &circuitInfo)
        {
            CircuitInfo::operator=(circuitInfo);
            return *this;
        }
 
        double getFieldValue(size_t index) const
        {
            if (index < 9)
            {
                return CircuitInfo::getFieldValue(index);
            }
            else
            {
                switch (index)
                {
                case 9:
                    return static_cast<double>(nrSamples);
                case 10:
                    return static_cast<double>(nrQubitsSampled);
                case 11:
                    return static_cast<double>(maxBondDim);
                case 12:
                    return static_cast<double>(nrPauliOps);
                case 13:
                    return static_cast<double>(executionCost);
                case 14:
                    return static_cast<double>(runtime);
                default:
                    throw std::out_of_range("Invalid index for ExecutionInfo field");
                }
            }
        }
    };
 
    static double EstimateExecutionCost(Simulators::SimulationType method, size_t nrQubits, const std::shared_ptr<Circuits::Circuit<>> &circuit,
                                        size_t maxBondDim)
    {
        if (method == Simulators::SimulationType::kPauliPropagator)
            return Simulators::QcsimPauliPropagator::GetCost(circuit);
        else if (method == Simulators::SimulationType::kStatevector)
        {
            const double opOrder = exp2(nrQubits);
 
            double cost = 0;
            for (const auto &op : *circuit)
            {
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement)
                    cost += opOrder * affectedQubits.size();
                if (op->GetType() == Circuits::OperationType::kReset)
                    cost += 2. * opOrder * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += opOrder;
                    else if (affectedQubits.size() == 2)
                        cost += 4 * opOrder;
                    else if (affectedQubits.size() == 3)
                        cost += 16 * opOrder;
                }
            }
            return cost;
        }
        else if (method == Simulators::SimulationType::kMatrixProductState)
        {
            const double twoQubitOpOrder = pow(maxBondDim, 3);
            const double oneQubitOpOrder = pow(maxBondDim, 2);
 
            double cost = 0;
            for (const auto &op : *circuit)
            {
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    // it's not that simple at all
                    // a qubit measurement works by applying a projector on the qubit
                    // tensor, which would cost as a one quibit gate but then we need to
                    // propagate the effect of the measurement along the chain, left and
                    // right, which is like applying a two qubit gate (SVD is the costlier
                    // operation there) on all the qubits that are entangled with the
                    // measured one, and in the worst case this can be all the other
                    // qubits if more than one qubits are measured at the same time, then
                    // we can have some optimizations, but let's just say that it's like
                    // measuring them one by one, so the cost is multiplied by the number
                    // of measured qubits
                    cost += twoQubitOpOrder * affectedQubits.size() * nrQubits / 2.;
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += oneQubitOpOrder;
                    else if (affectedQubits.size() == 2)
                        // I wish it were simple, but applying a gate involves swapping the
                        // qubits next to each other, then applying the gate
                        cost += twoQubitOpOrder * nrQubits / 3.;
                    else if (affectedQubits.size() == 3)
                        // qiskit aer has three qubit ops, qcsim has not, they are
                        // decomposed into one and two qubit gates
                        cost += twoQubitOpOrder * nrQubits * 2; // very rough estimation
                }
            }
            return cost;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            const double measOrder = pow(nrQubits, 2);
            const double opOrder = static_cast<double>(nrQubits);
 
            double cost = 0;
            for (const auto &op : *circuit)
            {
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    cost += measOrder * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                    cost += opOrder * affectedQubits.size();
            }
            return cost;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            double cost = 0;
            double doublingCost = 1;
 
            for (const auto &op : *circuit)
            {
                if (op->IsBranching())
                    doublingCost *= 2;
 
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    cost += doublingCost * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += doublingCost;
                    else if (affectedQubits.size() == 2)
                        cost += doublingCost * 2;
                    else if (affectedQubits.size() == 3)
                        cost += doublingCost * 4;
                }
            }
 
            return cost;
        }
 
        // for tensor network is hard to guess, it depends on contraction path
 
        return std::numeric_limits<double>::infinity();
    }
 
    static double EstimateSamplingCost(Simulators::SimulationType method, size_t nrQubits, size_t nrQubitsSampled, size_t samples,
                                       const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t maxBondDim)
    {
        // sampling works very differenty for such circuits
        Circuits::OperationState dummyState;
        const bool hasMeasurementsInTheMiddle = circuit->HasOpsAfterMeasurements();
        const std::vector<bool> executedOps = circuit->ExecuteNonMeasurements(nullptr, dummyState);
 
        const size_t dif = circuit->size() - executedOps.size();
 
        if (method == Simulators::SimulationType::kPauliPropagator)
            return Simulators::QcsimPauliPropagator::GetSamplingCost(circuit, nrQubitsSampled, samples);
 
        if (method == Simulators::SimulationType::kStatevector)
        {
            const double opOrder = exp2(nrQubits);
 
            double cost = 0;
            for (size_t i = 0; i < circuit->size(); ++i)
            {
                if (i >= dif && !executedOps[i - dif])
                    continue;
 
                const auto &op = (*circuit)[i];
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement)
                    cost += opOrder * affectedQubits.size();
                if (op->GetType() == Circuits::OperationType::kReset)
                    cost += 2. * opOrder * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += opOrder;
                    else if (affectedQubits.size() == 2)
                        cost += 4 * opOrder;
                    else if (affectedQubits.size() == 3)
                        cost += 16 * opOrder;
                }
            }
 
            // the sampling cost depends on the simulator
            // qiskit aer constucts an index gowing over the statevector, qcsim does
            // something similar to have the samples then O(nrQubits), but in maestro
            // for qcsim (and quest) this is replaced by alias sampling with O(1) for
            // each sample
 
            if (hasMeasurementsInTheMiddle)
            {
                double samplingCost = 0;
                for (size_t i = dif; i < circuit->size(); ++i)
                {
                    if (executedOps[i - dif])
                        continue;
 
                    const auto &op = (*circuit)[i];
                    const auto affectedQubits = op->AffectedQubits();
                    if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement)
                        samplingCost += opOrder * affectedQubits.size();
                    if (op->GetType() == Circuits::OperationType::kReset)
                        samplingCost += 2. * opOrder * affectedQubits.size();
                    else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                    {
                        if (affectedQubits.size() == 1)
                            samplingCost += opOrder;
                        else if (affectedQubits.size() == 2)
                            samplingCost += 4 * opOrder;
                        else if (affectedQubits.size() == 3)
                            samplingCost += 16 * opOrder;
                    }
                }
 
                samplingCost += opOrder * nrQubitsSampled;
 
                return cost + samplingCost * samples;
            }
 
            // some dummy cost, it's not going to fit all anyways
            return cost + 30 * opOrder + samples * nrQubits;
        }
        else if (method == Simulators::SimulationType::kMatrixProductState)
        {
            const double oneQubitOpOrder = static_cast<double>(maxBondDim * maxBondDim);
            const double twoQubitOpOrder = static_cast<double>(oneQubitOpOrder * maxBondDim);
 
            double cost = 0;
            for (size_t i = 0; i < circuit->size(); ++i)
            {
                if (i >= dif && !executedOps[i - dif])
                    continue;
 
                const auto &op = (*circuit)[i];
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    // it's not that simple at all
                    // a qubit measurement works by applying a projector on the qubit
                    // tensor, which would cost as a one quibit gate but then we need to
                    // propagate the effect of the measurement along the chain, left and
                    // right, which is like applying a two qubit gate (SVD is the costlier
                    // operation there) on all the qubits that are entangled with the
                    // measured one, and in the worst case this can be all the other
                    // qubits if more than one qubits are measured at the same time, then
                    // we can have some optimizations, but let's just say that it's like
                    // measuring them one by one, so the cost is multiplied by the number
                    // of measured qubits
                    cost += twoQubitOpOrder * affectedQubits.size() * nrQubits / 2.;
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += oneQubitOpOrder;
                    else if (affectedQubits.size() == 2)
                        // I wish it were simple, but applying a gate involves swapping the
                        // qubits next to each other, then applying the gate
                        cost += twoQubitOpOrder * nrQubits / 3.;
                    else if (affectedQubits.size() == 3)
                        // qiskit aer has three qubit ops, qcsim has not, they are
                        // decomposed into one and two qubit gates
                        cost += twoQubitOpOrder * nrQubits * 2; // very rough estimation
                }
            }
 
            if (hasMeasurementsInTheMiddle)
            {
                double samplingCost = 0;
                for (size_t i = dif; i < circuit->size(); ++i)
                {
                    if (executedOps[i - dif])
                        continue;
 
                    const auto &op = (*circuit)[i];
                    const auto affectedQubits = op->AffectedQubits();
                    if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                        op->GetType() == Circuits::OperationType::kReset)
                        samplingCost += twoQubitOpOrder * affectedQubits.size() * nrQubits / 2.;
                    else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                    {
                        if (affectedQubits.size() == 1)
                            samplingCost += oneQubitOpOrder;
                        else if (affectedQubits.size() == 2)
                            // I wish it were simple, but applying a gate involves swapping
                            // the qubits next to each other, then applying the gate
                            samplingCost += twoQubitOpOrder * nrQubits / 3.;
                        else if (affectedQubits.size() == 3)
                            // qiskit aer has three qubit ops, qcsim has not, they are
                            // decomposed into one and two qubit gates
                            samplingCost += twoQubitOpOrder * nrQubits * 2; // very rough estimation
                    }
                }
 
                samplingCost += twoQubitOpOrder * nrQubits * nrQubits;
 
                return cost + samplingCost * samples;
            }
 
            // sampling can be done here (and for other simulator types as well)
            // either by saving the state, measuring then restoring the state or
            // simply going along the chain, computing probabilities, throwing biased
            // coins and doing matrix multiplications (qcsim does this if all qubits
            // are sampled, otherwise the measurement method is used)
 
            // some dummy cost, it's not going to fit all anyways
            return cost + samples * twoQubitOpOrder * nrQubits * nrQubits;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            const double measOrder = pow(nrQubits, 2);
            const double opOrder = static_cast<double>(nrQubits);
 
            double cost = 0;
            for (size_t i = 0; i < circuit->size(); ++i)
            {
                if (i >= dif && !executedOps[i - dif])
                    continue;
 
                const auto &op = (*circuit)[i];
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    cost += measOrder * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                    cost += opOrder * affectedQubits.size();
            }
 
            if (hasMeasurementsInTheMiddle)
            {
                double samplingCost = 0;
                for (size_t i = dif; i < circuit->size(); ++i)
                {
                    if (executedOps[i - dif])
                        continue;
 
                    const auto &op = (*circuit)[i];
                    const auto affectedQubits = op->AffectedQubits();
                    if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                        op->GetType() == Circuits::OperationType::kReset)
                        samplingCost += measOrder * affectedQubits.size();
                    else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                        samplingCost += opOrder * affectedQubits.size();
                }
 
                samplingCost += measOrder * nrQubitsSampled;
 
                return cost + samplingCost * samples;
            }
 
            // TODO: to be changed to handle the all qubits sampling case, also the
            // sampling for qcsim which is efficient for sampling less than all
            // qubits!
 
            // sampling is done with saving state, measuring and restoring state
            // the overhead for saving / restoring is not here (it's of the same order
            // as measuring), but measurements are not all O(n^2) anyways
            return cost + samples * measOrder * nrQubitsSampled;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            double cost = 0;
            double doublingCost = 1;
            for (size_t i = 0; i < circuit->size(); ++i)
            {
                if (i >= dif && !executedOps[i - dif])
                    continue;
                const auto &op = (*circuit)[i];
                if (op->IsBranching())
                    doublingCost *= 2;
 
                const auto affectedQubits = op->AffectedQubits();
                if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                    op->GetType() == Circuits::OperationType::kReset)
                    cost += doublingCost * affectedQubits.size();
                else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                {
                    if (affectedQubits.size() == 1)
                        cost += doublingCost;
                    else if (affectedQubits.size() == 2)
                        cost += doublingCost * 2;
                    else if (affectedQubits.size() == 3)
                        cost += doublingCost * 4;
                }
            }
 
            if (hasMeasurementsInTheMiddle)
            {
                double samplingCost = 0;
                for (size_t i = dif; i < circuit->size(); ++i)
                {
                    if (executedOps[i - dif])
                        continue;
 
                    const auto &op = (*circuit)[i];
                    if (op->IsBranching())
                        doublingCost *= 2;
 
                    const auto affectedQubits = op->AffectedQubits();
 
                    if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement ||
                        op->GetType() == Circuits::OperationType::kReset)
                        samplingCost += doublingCost * affectedQubits.size();
                    else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate)
                    {
                        if (affectedQubits.size() == 1)
                            samplingCost += doublingCost;
                        else if (affectedQubits.size() == 2)
                            samplingCost += doublingCost * 2;
                        else if (affectedQubits.size() == 3)
                            samplingCost += doublingCost * 4;
                    }
                }
 
                samplingCost += doublingCost * nrQubitsSampled;
 
                return cost + samplingCost * samples;
            }
 
            // some dummy cost, it's not going to fit all anyways
            return cost + 30 * doublingCost + samples * nrQubits;
        }
 
        // for tensor network is hard to guess, it depends on contraction path
 
        return std::numeric_limits<double>::infinity();
    }
 
    static double EstimatePauliExpectationCost(const std::string &pauliString, Simulators::SimulationType method, size_t nrQubits,
                                               const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t maxBondDim)
    {
        // a pauli string propagated
        if (method == Simulators::SimulationType::kPauliPropagator)
            return Simulators::QcsimPauliPropagator::GetCost(circuit);
 
        double cost = EstimateExecutionCost(method, nrQubits, circuit, maxBondDim);
 
        size_t pauliCnt = 0;
        for (char c : pauliString)
        {
            if (c == 'X' || c == 'x' || c == 'Y' || c == 'y' || c == 'Z' || c == 'z')
                ++pauliCnt;
        }
 
        if (method == Simulators::SimulationType::kStatevector)
        {
            const double opOrder = exp2(nrQubits);
 
            // each pauli matrix is applied to the statevector, then the product with
            // the original is computed
            cost += opOrder * (pauliCnt + 1);
 
            return cost;
        }
        else if (method == Simulators::SimulationType::kMatrixProductState)
        {
            const double twoQubitOpOrder = pow(maxBondDim, 3);
            const double oneQubitOpOrder = pow(maxBondDim, 2);
 
            // each pauli is applied to the chain then the resulted chain is
            // contracted with the original one
            cost += oneQubitOpOrder * pauliCnt + nrQubits * twoQubitOpOrder;
 
            return cost;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            cost += nrQubits * nrQubits;
            return cost;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            double doublingCost = 1;
            for (const auto &op : *circuit)
            {
                if (op->IsBranching())
                    doublingCost *= 2;
            }
            cost += 4 * doublingCost * pauliCnt;
 
            return cost;
        }
 
        return std::numeric_limits<double>::infinity();
    }
 
    static std::shared_ptr<Circuits::Circuit<>> GenerateRandomCircuit(size_t nrQubits, size_t depth, double measureInsideProbability = 0.,
                                                                      size_t nrMeasAtEnd = 0, bool isClifford = false, size_t nrNonCliffordGatesLimit = 0,
                                                                      size_t nrBranchingGatesLimit = 0)
    {
        auto circuit = std::make_shared<Circuits::Circuit<>>();
        std::random_device rdev;
        std::mt19937 rng(rdev());
        std::uniform_real_distribution<double> dist(0.0, 1.0);
        std::uniform_real_distribution<double> paramDist(-2 * M_PI, 2 * M_PI);
        std::uniform_int_distribution<Types::qubit_t> qubitDist(0, nrQubits - 1);
        std::uniform_int_distribution<int> gateDist(0, static_cast<int>(Circuits::QuantumGateType::kCCXGateType));
 
        std::uniform_int_distribution<int> gateDistOneQubit(0, static_cast<int>(Circuits::QuantumGateType::kUGateType));
 
        std::vector<Types::qubit_t> qubits(nrQubits);
        std::iota(qubits.begin(), qubits.end(), 0);
 
        size_t nrNonCliffordGates = 0;
        if (nrNonCliffordGatesLimit == 0 && !isClifford)
            nrNonCliffordGatesLimit = depth;
 
        size_t nrBranchingGates = 0;
        if (nrBranchingGatesLimit == 0)
            nrBranchingGatesLimit = depth;
 
        for (size_t i = 0; i < depth; ++i)
        {
            if (dist(rng) < measureInsideProbability)
            {
                Types::qubit_t q = qubitDist(rng);
                std::vector<std::pair<Types::qubit_t, size_t>> qs = {{q, q}};
                circuit->AddOperation(std::make_shared<Circuits::MeasurementOperation<>>(qs));
                continue;
            }
 
            std::shuffle(qubits.begin(), qubits.end(), rng);
            const auto q1 = qubits[0];
            const auto q2 = qubits[1];
            const auto q3 = qubits[2];
 
            auto gateType = static_cast<Circuits::QuantumGateType>(nrNonCliffordGatesLimit < depth ? gateDistOneQubit(rng) // avoid three qubit gates, they are
                                                                                                                           // non-clifford and they cost a lot
                                                                                                   : gateDist(rng));
            auto param1 = paramDist(rng);
            auto param2 = paramDist(rng);
            auto param3 = paramDist(rng);
            auto param4 = paramDist(rng);
 
            auto theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
            if (isClifford)
            {
                while (!theGate->IsClifford())
                {
                    gateType = static_cast<Circuits::QuantumGateType>(gateDist(rng));
                    theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
                }
            }
 
            if (!theGate->IsClifford())
                ++nrNonCliffordGates;
            if (theGate->IsBranching())
                ++nrBranchingGates;
 
            if (nrNonCliffordGates > nrNonCliffordGatesLimit)
            {
                // replace the non clifford gate with a clifford one
                gateType = static_cast<Circuits::QuantumGateType>(gateDistOneQubit(rng));
                theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
                while (!theGate->IsClifford())
                {
                    gateType = static_cast<Circuits::QuantumGateType>(gateDistOneQubit(rng));
                    theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
                }
                --nrNonCliffordGates;
            }
 
            if (nrBranchingGates > nrBranchingGatesLimit)
            {
                // replace the branching gate with a non branching one
                gateType = static_cast<Circuits::QuantumGateType>(gateDist(rng));
 
                theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
                while (theGate->IsBranching())
                {
                    gateType = static_cast<Circuits::QuantumGateType>(gateDist(rng));
                    theGate = Circuits::CircuitFactory<>::CreateGate(gateType, q1, q2, q3, param1, param2, param3, param4);
                }
                --nrBranchingGates;
            }
 
            circuit->AddOperation(theGate);
        }
 
        std::shuffle(circuit->begin(), circuit->end(), rng);
 
        if (nrMeasAtEnd > 0)
        {
            if (nrMeasAtEnd > nrQubits)
                nrMeasAtEnd = nrQubits;
            std::shuffle(qubits.begin(), qubits.end(), rng);
 
            for (size_t i = 0; i < nrMeasAtEnd; ++i)
            {
                std::vector<std::pair<Types::qubit_t, size_t>> qs = {{qubits[i], i}};
                circuit->AddOperation(std::make_shared<Circuits::MeasurementOperation<>>(qs));
            }
        }
 
        return circuit;
    }
 
    static double MeasureExecutionTime(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrQubits,
                                       const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t nrReps, size_t maxBondDim)
    {
        auto sim = GetSimulator(simType, method, nrQubits, maxBondDim);
 
        Circuits::OperationState dummyState(nrQubits);
        auto start = std::chrono::high_resolution_clock::now();
        for (size_t i = 0; i < nrReps; ++i)
            circuit->Execute(sim, dummyState);
        auto end = std::chrono::high_resolution_clock::now();
        return std::chrono::duration<double>(end - start).count();
    }
 
    static double MeasureSamplingTime(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrQubits,
                                      const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t nrQubitsSampled, size_t nrSamples, size_t nrReps,
                                      size_t maxBondDim)
    {
        auto sim = GetSimulator(simType, method, nrQubits, maxBondDim);
        Circuits::OperationState dummyState(nrQubits);
 
        if (nrQubitsSampled > nrQubits)
            nrQubitsSampled = nrQubits;
        Types::qubits_vector qubitsSampled(nrQubitsSampled);
        std::iota(qubitsSampled.begin(), qubitsSampled.end(), 0);
 
        const bool hasMeasurementsInTheMiddle = circuit->HasOpsAfterMeasurements();
 
        auto start = std::chrono::high_resolution_clock::now();
 
        if (hasMeasurementsInTheMiddle)
        {
            for (size_t i = 0; i < nrReps; ++i)
            {
                const auto executedOps = circuit->ExecuteNonMeasurements(sim, dummyState); // execute the circuit up to measurements
 
                // now sample
                for (size_t sample = 0; sample < nrSamples; ++sample)
                {
                    circuit->ExecuteMeasurements(sim, dummyState, executedOps); // execute the measurements
                }
            }
        }
        else
        {
            for (size_t i = 0; i < nrReps; ++i)
            {
                circuit->Execute(sim, dummyState); // execute the circuit first
                sim->SampleCountsMany(qubitsSampled, nrSamples);
            }
        }
        auto end = std::chrono::high_resolution_clock::now();
 
        return std::chrono::duration<double>(end - start).count();
    }
 
    static double MeasurePauliExpectationTime(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrQubits,
                                              const std::shared_ptr<Circuits::Circuit<>> &circuit, const std::string &pauliString, size_t nrReps,
                                              size_t maxBondDim)
    {
        auto sim = GetSimulator(simType, method, nrQubits, maxBondDim);
        Circuits::OperationState dummyState(nrQubits);
        auto start = std::chrono::high_resolution_clock::now();
        for (size_t i = 0; i < nrReps; ++i)
        {
            circuit->Execute(sim, dummyState); // execute the circuit first
            sim->ExpectationValue(pauliString);
        }
        auto end = std::chrono::high_resolution_clock::now();
        return std::chrono::duration<double>(end - start).count();
    }
 
    static std::shared_ptr<Simulators::ISimulator> GetSimulator(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrQubits,
                                                                size_t maxBondDim)
    {
        auto sim = Simulators::SimulatorsFactory::CreateSimulator(simType, method);
        if (method == Simulators::SimulationType::kMatrixProductState)
        {
            const auto strVal = std::to_string(maxBondDim);
            sim->Configure("matrix_product_state_max_bond_dimension", strVal.c_str());
        }
        sim->AllocateQubits(nrQubits);
        sim->SetMultithreading(false);
        sim->Initialize();
        return sim;
    }
 
    static CircuitInfo GetCircuitInfo(const std::shared_ptr<Circuits::Circuit<>> &circuit)
    {
        CircuitInfo info;
 
        info.nrQubits = circuit->GetMaxQubitIndex() + 1;
        Circuits::OperationState dummyState(info.nrQubits);
        const bool hasMeasurementsInTheMiddle = circuit->HasOpsAfterMeasurements();
        const std::vector<bool> executedOps = circuit->ExecuteNonMeasurements(nullptr, dummyState);
 
        const size_t dif = circuit->size() - executedOps.size();
 
        size_t i = 0;
        for (const auto &op : *circuit)
        {
            const auto affectedQubits = op->AffectedQubits();
            if (op->GetType() == Circuits::OperationType::kMeasurement || op->GetType() == Circuits::OperationType::kConditionalMeasurement)
            {
                if (hasMeasurementsInTheMiddle)
                    ++info.nrMiddleMeasurementOps;
                else
                    ++info.nrEndMeasurementOps;
            }
            else if (op->GetType() == Circuits::OperationType::kGate || op->GetType() == Circuits::OperationType::kConditionalGate ||
                     op->GetType() == Circuits::OperationType::kQuantumChannel)
            {
                if (affectedQubits.size() == 1)
                {
                    ++info.nrOneQubitOps;
                    if (i < dif || executedOps[i - dif])
                        ++info.nrOneQubitOpsExecutedOnce;
                }
                else if (affectedQubits.size() == 2)
                {
                    ++info.nrTwoQubitOps;
                    if (i < dif || executedOps[i - dif])
                        ++info.nrTwoQubitOpsExecutedOnce;
                }
                else if (affectedQubits.size() == 3)
                {
                    ++info.nrThreeQubitOps;
                    if (i < dif || executedOps[i - dif])
                        ++info.nrThreeQubitOpsExecutedOnce;
                }
            }
            ++i;
        }
 
        return info;
    }
 
    static std::vector<ExecutionInfo> ReadLog(const std::string &logFilePath)
    {
        std::vector<ExecutionInfo> executionInfos;
 
        std::ifstream logFile(logFilePath);
        if (!logFile.is_open())
        {
            std::cerr << "Failed to open log file: " << logFilePath << std::endl;
            return executionInfos;
        }
 
        std::string line;
 
        while (std::getline(logFile, line))
        {
            std::stringstream ss(line);
            std::string value;
            ExecutionInfo info;
            std::getline(ss, value, ',');
            info.nrQubits = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrOneQubitOps = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrTwoQubitOps = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrThreeQubitOps = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrMiddleMeasurementOps = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrEndMeasurementOps = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrOneQubitOpsExecutedOnce = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrTwoQubitOpsExecutedOnce = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrThreeQubitOpsExecutedOnce = std::stoul(value);
 
            std::getline(ss, value, ',');
            info.nrSamples = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrQubitsSampled = std::stoul(value);
            std::getline(ss, value, ',');
            info.maxBondDim = std::stoul(value);
            std::getline(ss, value, ',');
            info.nrPauliOps = std::stoul(value);
 
            std::getline(ss, value, ',');
            info.executionCost = std::stod(value);
            std::getline(ss, value, ',');
            info.runtime = std::stod(value);
 
            if (info.nrSamples < 1)
                info.nrSamples = 1;
 
            executionInfos.push_back(std::move(info));
        }
 
        std::random_device rdev;
        std::mt19937 rng(rdev());
        std::shuffle(executionInfos.begin(), executionInfos.end(), rng);
 
        return executionInfos;
    }
 
    static void BenchmarkAndLogExecution(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrReps, size_t nrMinQubits,
                                         size_t nrMaxQubits, size_t stepQubits, size_t depthMin, size_t depthMax, size_t stepDepth,
                                         double measureInsideProbability, size_t nrMeasAtEndMin, size_t nrMeasAtEndMax, size_t stepMeasAtEnd,
                                         size_t nrRandomCircuitsPerConfig, const std::string &logFilePath, size_t startBondDim = 16, size_t endBondDim = 16)
    {
        bool isClifford = (method == Simulators::SimulationType::kStabilizer);
        int nrNonCliffordGates = static_cast<int>(depthMax); // a limit
        int nrBranchingGates = static_cast<int>(depthMax);
 
        // but if it's pauli...
        if (method == Simulators::SimulationType::kPauliPropagator)
        {
            nrNonCliffordGates = 1;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            nrNonCliffordGates = 0;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            nrBranchingGates = 8;
        }
 
        std::cout << "Benchmarking execution for simType: " << static_cast<int>(simType) << ", method: " << static_cast<int>(method) << std::endl;
 
        Utils::LogFile log(logFilePath);
 
        for (size_t nrQubits = nrMinQubits; nrQubits <= nrMaxQubits; nrQubits += stepQubits)
        {
            std::cout << "  Qubits: " << nrQubits << std::endl;
            for (size_t depth = depthMin; depth <= depthMax; depth += stepDepth)
            {
                std::cout << "    Depth: " << depth << std::endl;
                for (size_t nrMeasAtEnd = nrMeasAtEndMin; nrMeasAtEnd <= std::min(nrMeasAtEndMax, nrQubits); nrMeasAtEnd += stepMeasAtEnd)
                {
                    std::cout << "      Measurements at end: " << nrMeasAtEnd << std::endl;
                    for (size_t i = 0; i < nrRandomCircuitsPerConfig; ++i)
                    {
                        std::cout << "        Random circuit: " << i + 1 << "/" << nrRandomCircuitsPerConfig << std::endl;
                        if (method == Simulators::SimulationType::kPauliPropagator)
                            isClifford = !isClifford;
 
                        const auto circuit =
                            GenerateRandomCircuit(nrQubits, depth, measureInsideProbability, nrMeasAtEnd, isClifford, nrNonCliffordGates, nrBranchingGates);
 
                        for (size_t maxBondDim = startBondDim; maxBondDim <= endBondDim; maxBondDim *= 2)
                        {
                            BenchmarkAndLogExecution(simType, method, circuit, nrReps, maxBondDim, log);
                        }
                    }
                }
            }
        }
    }
 
    static std::string GeneratePauliString(size_t nrQubits)
    {
        std::string pauli;
        pauli.resize(nrQubits);
        std::random_device rd;
        std::mt19937 g(rd());
        std::uniform_int_distribution<int> dist(0, 3);
 
        for (size_t i = 0; i < nrQubits; ++i)
        {
            const int v = dist(g);
            switch (v)
            {
            case 0:
                pauli[i] = 'I';
                break;
            case 1:
                pauli[i] = 'X';
                break;
            case 2:
                pauli[i] = 'Y';
                break;
            case 3:
                pauli[i] = 'Z';
                break;
            }
        }
 
        return pauli;
    }
 
    static void BenchmarkAndLogPauliExpectation(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrReps, size_t nrMinQubits,
                                                size_t nrMaxQubits, size_t stepQubits, size_t depthMin, size_t depthMax, size_t stepDepth,
                                                size_t nrRandomCircuitsPerConfig, const std::string &logFilePath, size_t startBondDim = 16,
                                                size_t endBondDim = 16)
    {
        bool isClifford = (method == Simulators::SimulationType::kStabilizer);
        int nrNonCliffordGates = static_cast<int>(depthMax); // a limit
        int nrBranchingGates = static_cast<int>(depthMax);
 
        // but if it's pauli...
        if (method == Simulators::SimulationType::kPauliPropagator)
        {
            nrNonCliffordGates = 1;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            nrNonCliffordGates = 0;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            nrBranchingGates = 8;
        }
        Utils::LogFile log(logFilePath);
 
        std::cout << "Benchmarking Pauli expectation for simType: " << static_cast<int>(simType) << ", method: " << static_cast<int>(method) << std::endl;
 
        for (size_t nrQubits = nrMinQubits; nrQubits <= nrMaxQubits; nrQubits += stepQubits)
        {
            std::cout << "  Qubits: " << nrQubits << std::endl;
            for (size_t depth = depthMin; depth <= depthMax; depth += stepDepth)
            {
                std::cout << "    Depth: " << depth << std::endl;
                for (size_t i = 0; i < nrRandomCircuitsPerConfig; ++i)
                {
                    if (method == Simulators::SimulationType::kPauliPropagator)
                        isClifford = !isClifford;
 
                    const auto circuit = GenerateRandomCircuit(nrQubits, depth, 0., 0, isClifford, nrNonCliffordGates, nrBranchingGates);
                    const std::string pauliString = GeneratePauliString(nrQubits);
 
                    std::cout << "      Random circuit: " << i + 1 << "/" << nrRandomCircuitsPerConfig << " Pauli string: " << pauliString << std::endl;
                    for (size_t maxBondDim = startBondDim; maxBondDim <= endBondDim; maxBondDim *= 2)
                    {
                        std::cout << "        Max bond dimension: " << maxBondDim << std::endl;
                        BenchmarkAndLogPauliExpectation(simType, method, circuit, pauliString, nrReps, maxBondDim, log);
                    }
                }
            }
        }
    }
 
    static void BenchmarkAndLogSampling(Simulators::SimulatorType simType, Simulators::SimulationType method, size_t nrReps, size_t nrMinQubits,
                                        size_t nrMaxQubits, size_t stepQubits, size_t depthMin, size_t depthMax, size_t stepDepth, size_t nrMeasAtEndMin,
                                        size_t nrMeasAtEndMax, size_t stepMeasAtEnd, size_t nrSamplesMin, size_t nrSamplesMax, size_t multiplierSamples,
                                        size_t nrRandomCircuitsPerConfig, const std::string &logFilePath, size_t startBondDim = 16, size_t endBondDim = 16)
    {
        bool isClifford = (method == Simulators::SimulationType::kStabilizer);
        int nrNonCliffordGates = static_cast<int>(depthMax); // a limit
        int nrBranchingGates = static_cast<int>(depthMax);
        // but if it's pauli...
        if (method == Simulators::SimulationType::kPauliPropagator)
        {
            nrNonCliffordGates = 1;
        }
        else if (method == Simulators::SimulationType::kStabilizer)
        {
            nrNonCliffordGates = 0;
        }
        else if (method == Simulators::SimulationType::kPathIntegral)
        {
            nrBranchingGates = 8;
        }
        Utils::LogFile log(logFilePath);
 
        std::cout << "Benchmarking sampling for simType: " << static_cast<int>(simType) << ", method: " << static_cast<int>(method) << std::endl;
 
        for (size_t nrQubits = nrMinQubits; nrQubits <= nrMaxQubits; nrQubits += stepQubits)
        {
            std::cout << "  Qubits: " << nrQubits << std::endl;
            for (size_t depth = depthMin; depth <= depthMax; depth += stepDepth)
            {
                std::cout << "    Depth: " << depth << std::endl;
                for (size_t nrSamples = nrSamplesMin; nrSamples <= nrSamplesMax; nrSamples *= multiplierSamples)
                {
                    std::cout << "      Samples: " << nrSamples << std::endl;
                    for (size_t i = 0; i < nrRandomCircuitsPerConfig; ++i)
                    {
                        if (method == Simulators::SimulationType::kPauliPropagator)
                            isClifford = !isClifford;
                        const auto circuit = GenerateRandomCircuit(nrQubits, depth, 0., 0, isClifford, nrNonCliffordGates, nrBranchingGates);
                        for (size_t nrQubitsSampled = nrQubits; nrQubitsSampled >= 1; nrQubitsSampled /= 2)
                        {
                            std::cout << "        Random circuit: " << i + 1 << "/" << nrRandomCircuitsPerConfig << " Qubits sampled: " << nrQubitsSampled
                                      << std::endl;
                            for (size_t maxBondDim = startBondDim; maxBondDim <= endBondDim; maxBondDim *= 2)
                            {
                                BenchmarkAndLogSampling(simType, method, circuit, nrQubitsSampled, nrSamples, nrReps, maxBondDim, log);
                            }
                        }
                    }
                }
            }
        }
    }
 
    static void BenchmarkAndLogExecution(Simulators::SimulatorType simType, Simulators::SimulationType method,
                                         const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t nrReps, size_t maxBondDim, Utils::LogFile &log)
    {
        const auto info = GetCircuitInfo(circuit);
        const double estimatedCost = EstimateExecutionCost(method, info.nrQubits, circuit, maxBondDim);
        const double executionTime = MeasureExecutionTime(simType, method, info.nrQubits, circuit, nrReps, maxBondDim) / nrReps;
 
        std::stringstream ss;
 
        ss << info.nrQubits << "," << info.nrOneQubitOps << "," << info.nrTwoQubitOps << "," << info.nrThreeQubitOps << "," << info.nrMiddleMeasurementOps
           << "," << info.nrEndMeasurementOps << "," << info.nrOneQubitOpsExecutedOnce << "," << info.nrTwoQubitOpsExecutedOnce << ","
           << info.nrThreeQubitOpsExecutedOnce << ","
           << "0,0," << maxBondDim << ","
           << "0," << estimatedCost << "," << executionTime;
 
        log.Log(ss.str());
    }
 
    static void BenchmarkAndLogSampling(Simulators::SimulatorType simType, Simulators::SimulationType method,
                                        const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t nrQubitsSampled, size_t nrSamples, size_t nrReps,
                                        size_t maxBondDim, Utils::LogFile &log)
    {
        const auto info = GetCircuitInfo(circuit);
        const double estimatedCost = EstimateSamplingCost(method, info.nrQubits, nrQubitsSampled, nrSamples, circuit, maxBondDim);
        const double samplingTime = MeasureSamplingTime(simType, method, info.nrQubits, circuit, nrQubitsSampled, nrSamples, nrReps, maxBondDim) / nrReps;
        std::stringstream ss;
        ss << info.nrQubits << "," << info.nrOneQubitOps << "," << info.nrTwoQubitOps << "," << info.nrThreeQubitOps << "," << info.nrMiddleMeasurementOps
           << "," << info.nrEndMeasurementOps << "," << info.nrOneQubitOpsExecutedOnce << "," << info.nrTwoQubitOpsExecutedOnce << ","
           << info.nrThreeQubitOpsExecutedOnce << "," << nrSamples << "," << nrQubitsSampled << "," << maxBondDim << ","
           << "0," << estimatedCost << "," << samplingTime;
        log.Log(ss.str());
    }
 
    static void BenchmarkAndLogPauliExpectation(Simulators::SimulatorType simType, Simulators::SimulationType method,
                                                const std::shared_ptr<Circuits::Circuit<>> &circuit, const std::string &pauliString, size_t nrReps,
                                                size_t maxBondDim, Utils::LogFile &log)
    {
        auto info = GetCircuitInfo(circuit);
        info.nrQubits = std::max(info.nrQubits, pauliString.size());
        const double estimatedCost = EstimatePauliExpectationCost(pauliString, method, info.nrQubits, circuit, maxBondDim);
        const double expectationTime = MeasurePauliExpectationTime(simType, method, info.nrQubits, circuit, pauliString, nrReps, maxBondDim) / nrReps;
 
        size_t cntPauli = 0;
        for (char c : pauliString)
        {
            if (c == 'X' || c == 'x' || c == 'Y' || c == 'y' || c == 'Z' || c == 'z')
                ++cntPauli;
        }
 
        std::stringstream ss;
        ss << info.nrQubits << "," << info.nrOneQubitOps << "," << info.nrTwoQubitOps << "," << info.nrThreeQubitOps << "," << info.nrMiddleMeasurementOps
           << "," << info.nrEndMeasurementOps << "," << info.nrOneQubitOpsExecutedOnce << "," << info.nrTwoQubitOpsExecutedOnce << ","
           << info.nrThreeQubitOpsExecutedOnce << ","
           << "0,0," << maxBondDim << "," << cntPauli << "," << estimatedCost << "," << expectationTime;
        log.Log(ss.str());
    }
 
    static std::shared_ptr<Utils::MultipleLinearRegression> GetRegressor(const std::string &logFilePath, const std::vector<size_t> &featureIndices)
    {
        const auto executionInfos = ReadLog(logFilePath);
 
        std::vector<std::vector<double>> features;
        features.reserve(executionInfos.size());
 
        std::vector<double> targetValues;
        targetValues.reserve(executionInfos.size());
 
        for (const auto &info : executionInfos)
        {
            std::vector<double> featureVector(featureIndices.size());
            for (size_t i = 0; i < featureVector.size(); ++i)
                featureVector[i] = info.getFieldValue(featureIndices[i]);
            features.push_back(std::move(featureVector));
 
            targetValues.push_back(info.runtime);
        }
 
        auto regressor = std::make_shared<Utils::MultipleLinearRegression>();
        regressor->SetSamples(features, targetValues);
 
        return regressor;
    }
};
 
} // namespace Estimators
 
#endif