AerSimulator.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _AER_SIMULATOR_H_
#define _AER_SIMULATOR_H_
 
#ifndef NO_QISKIT_AER
 
#ifdef INCLUDED_BY_FACTORY
 
#include "AerState.h"
 
namespace Simulators
{
// TODO: Maybe use the pimpl idiom
// https://en.cppreference.com/w/cpp/language/pimpl to hide the implementation
// for good but during development this should be good enough
namespace Private
{

class AerSimulator : public AerState
{
  public:
    AerSimulator() = default;
    // allow no copy or assignment
    AerSimulator(const AerSimulator &) = delete;
    AerSimulator &operator=(const AerSimulator &) = delete;
 
    // but allow moving
    AerSimulator(AerSimulator &&other) = default;
    AerSimulator &operator=(AerSimulator &&other) = default;

    void ApplyGenericOneQubitGate(Types::qubit_t qubit, const Eigen::Matrix2cd &gate) override
    {
        const AER::reg_t qubits = {qubit};
        AER::cmatrix_t gate_matrix(2, 2);
        gate_matrix(0, 0) = gate(0, 0);
        gate_matrix(0, 1) = gate(0, 1);
        gate_matrix(1, 0) = gate(1, 0);
        gate_matrix(1, 1) = gate(1, 1);
 
        state->apply_unitary(qubits, gate_matrix);
        NotifyObservers(qubits);
    }

    void ApplyGenericTwoQubitGate(Types::qubit_t qubit0, Types::qubit_t qubit1, const Eigen::Matrix4cd &gate) override
    {
        // Aer, like QCSim, maps the first entry to local matrix bit zero.
        const AER::reg_t qubits = {qubit0, qubit1};
        AER::cmatrix_t gate_matrix(4, 4);
        for (size_t i = 0; i < 4; ++i)
            for (size_t j = 0; j < 4; ++j)
                gate_matrix(i, j) = gate(i, j);
 
        state->apply_unitary(qubits, gate_matrix);
        NotifyObservers(qubits);
    }

    void ApplyP(Types::qubit_t qubit, double lambda) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("P gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "p";
            op.qubits = AER::reg_t{qubit};
            op.params = {lambda};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            const AER::cvector_t P = {{1, 0}, std::polar(1., lambda)};
            state->apply_diagonal_matrix(qubits, P);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyX(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
        state->apply_x(qubit);
        NotifyObservers(qubits);
    }

    void ApplyY(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
        state->apply_y(qubit);
        NotifyObservers(qubits);
    }

    void ApplyZ(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
        state->apply_z(qubit);
        NotifyObservers(qubits);
    }

    void ApplyH(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
        state->apply_h(qubit);
        NotifyObservers(qubits);
    }

    void ApplyS(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "s";
            op.qubits = AER::reg_t{qubit};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            const AER::cvector_t S = {{1, 0}, complex_t(0.0, 1)};
            state->apply_diagonal_matrix(qubits, S);
        }
 
        NotifyObservers(qubits);
    }

    void ApplySDG(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "sdg";
            op.qubits = AER::reg_t{qubit};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            const AER::cvector_t Sdg = {{1, 0}, complex_t(0.0, -1)};
            state->apply_diagonal_matrix(qubits, Sdg);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyT(Types::qubit_t qubit) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("T gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
        // state->apply_u(qubit, 0, 0, M_PI / 4.0);
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "t";
            op.qubits = AER::reg_t{qubit};
            state->buffer_op(std::move(op));
        }
        else
        {
            const AER::cvector_t T = {{1, 0}, std::polar<double>(1, M_PI / 4.)};
            state->apply_diagonal_matrix(qubits, T);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyTDG(Types::qubit_t qubit) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("TDG gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
        // state->apply_u(qubit, 0, 0, -M_PI / 4.0);
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "tdg";
            op.qubits = AER::reg_t{qubit};
            state->buffer_op(std::move(op));
        }
        else
        {
            const AER::cvector_t Tdg = {{1, 0}, std::polar<double>(1, -M_PI / 4.)};
            state->apply_diagonal_matrix(qubits, Tdg);
        }
 
        NotifyObservers(qubits);
    }

    void ApplySx(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "sx";
            op.qubits = AER::reg_t{qubit};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            // there is a difference in the global phase for this
            // changed to match qcsim behavior (and also the qiskit docs and
            // implementation - not the one exposed in 'contrib', though)
            // state->apply_mcrx({ qubit }, -M_PI / 4.0);
            state->apply_unitary(qubits, AER::Linalg::Matrix::SX);
        }
 
        NotifyObservers(qubits);
    }

    void ApplySxDAG(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "sxdg";
            op.qubits = AER::reg_t{qubit};
 
            state->buffer_op(std::move(op));
        }
        else
            state->apply_unitary(qubits, AER::Linalg::Matrix::SXDG);
 
        NotifyObservers(qubits);
    }

    void ApplyK(Types::qubit_t qubit) override
    {
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyZ(qubit);
            ApplyS(qubit);
            ApplyH(qubit);
            ApplyS(qubit);
        }
        else
        {
            static const cmatrix_t K =
                AER::Utils::make_matrix<complex_t>({{{1. / std::sqrt(2.), 0}, {0, -1. / std::sqrt(2.)}}, {{0, 1. / std::sqrt(2.)}, {-1. / std::sqrt(2.), 0}}});
 
            state->apply_unitary(qubits, K);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyRx(Types::qubit_t qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("Rx gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyH(qubit);
            ApplyRz(qubit, theta);
            ApplyH(qubit);
        }
        else
        {
            const cmatrix_t rx = AER::Linalg::Matrix::rx(theta);
            state->apply_unitary(qubits, rx);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyRy(Types::qubit_t qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("Ry gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplySDG(qubit);
            ApplyRx(qubit, theta);
            ApplyS(qubit);
        }
        else
        {
            const cmatrix_t ry = AER::Linalg::Matrix::ry(theta);
 
            state->apply_unitary(qubits, ry);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyRz(Types::qubit_t qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("Rz gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "rz";
            op.qubits = AER::reg_t{qubit};
            op.params = {theta};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            const cmatrix_t rz = AER::Linalg::Matrix::rz(theta);
 
            state->apply_unitary(qubits, rz);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyU(Types::qubit_t qubit, double theta, double phi, double lambda, double gamma) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("U gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyRz(qubit, lambda);
            ApplyRy(qubit, theta);
            ApplyRz(qubit, phi);
        }
        else
        {
            const cmatrix_t u = AER::Linalg::Matrix::u4(theta, phi, lambda, gamma);
 
            state->apply_unitary(qubits, u);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCX(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
        state->apply_cx(qubits);
        NotifyObservers(qubits);
    }

    void ApplyCY(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplySDG(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyS(tgt_qubit);
        }
        else
        {
            state->apply_cy(qubits);
        }
        NotifyObservers(qubits);
    }

    void ApplyCZ(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
        state->apply_cz(qubits);
        NotifyObservers(qubits);
    }

    void ApplyCP(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double lambda) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CP gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            const double halfAngle = lambda * 0.5;
            ApplyP(ctrl_qubit, halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyP(tgt_qubit, -halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyP(tgt_qubit, halfAngle);
        }
        else
        {
            const cmatrix_t CP = AER::Linalg::Matrix::cphase(lambda);
            // state->apply_mcphase(qubits, lambda);
            state->apply_unitary(qubits, CP);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCRx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CRx gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            const double halfAngle = theta * 0.5;
 
            ApplyH(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyRz(tgt_qubit, -halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyRz(tgt_qubit, halfAngle);
            ApplyH(tgt_qubit);
        }
        else
        {
            cmatrix_t mat(4, 4);
            mat(0, 0) = 1;
            mat(2, 2) = 1;
 
            const double t2 = theta * 0.5;
 
            const complex_t i(0., 1.);
            mat(1, 1) = std::cos(t2);
            mat(1, 3) = -i * std::sin(t2);
            mat(3, 1) = mat(1, 3);
            mat(3, 3) = mat(1, 1);
 
            // state->apply_mcrx(qubits, theta);
            state->apply_unitary(qubits, mat);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCRy(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CRy gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            const double halfAngle = theta * 0.5;
 
            ApplyRy(tgt_qubit, halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyRy(tgt_qubit, -halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
        }
        else
        {
            cmatrix_t mat(4, 4);
            mat(0, 0) = 1;
            mat(2, 2) = 1;
 
            const double t2 = theta * 0.5;
 
            mat(1, 1) = std::complex<double>(cos(t2), 0);
            mat(1, 3) = std::complex<double>(-sin(t2), 0);
            mat(3, 1) = std::complex<double>(sin(t2), 0);
            mat(3, 3) = mat(1, 1);
 
            // state->apply_mcry(qubits, theta);
            state->apply_unitary(qubits, mat);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCRz(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CRz gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            const double halfAngle = theta * 0.5;
            ApplyRz(tgt_qubit, halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyRz(tgt_qubit, -halfAngle);
            ApplyCX(ctrl_qubit, tgt_qubit);
        }
        else
        {
            const double t2 = theta * 0.5;
            AER::cvector_t v = {{1, 0}, std::polar(1., -t2), {1, 0}, std::polar(1., t2)};
 
            // state->apply_mcrz(qubits, theta);
            state->apply_diagonal_matrix(qubits, v);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCH(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CH gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyH(tgt_qubit);
            ApplySDG(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyH(tgt_qubit);
            ApplyT(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyT(tgt_qubit);
            ApplyH(tgt_qubit);
            ApplyS(tgt_qubit);
            ApplyX(tgt_qubit);
            ApplyS(ctrl_qubit);
        }
        else
        {
            const cmatrix_t CU = AER::Linalg::Matrix::cu(M_PI_2, 0, M_PI, 0);
            state->apply_unitary(qubits, CU);
            // state->apply_cu(qubits, M_PI_2, 0, M_PI, 0);
        }
        NotifyObservers(qubits);
    }

    void ApplyCSx(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CSx gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyH(tgt_qubit);
            ApplyT(ctrl_qubit);
            ApplyT(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyTDG(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyH(tgt_qubit);
        }
        else
        {
            static const cmatrix_t CSX = AER::Utils::make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0.5, 0.5}, {0, 0}, {0.5, -0.5}},
                                                                             {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
                                                                             {{0, 0}, {0.5, -0.5}, {0, 0}, {0.5, 0.5}}});
 
            state->apply_unitary(qubits, CSX);
        }
        NotifyObservers(qubits);
    }

    void ApplyCSxDAG(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CSxDAG gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            ApplyH(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyT(tgt_qubit);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyTDG(ctrl_qubit);
            ApplyTDG(tgt_qubit);
            ApplyH(tgt_qubit);
        }
        else
        {
            static const cmatrix_t CSXDG = AER::Utils::make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                               {{0, 0}, {0.5, -0.5}, {0, 0}, {0.5, 0.5}},
                                                                               {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
                                                                               {{0, 0}, {0.5, 0.5}, {0, 0}, {0.5, -0.5}}});
 
            state->apply_unitary(qubits, CSXDG);
        }
 
        NotifyObservers(qubits);
    }

    void ApplySwap(Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        const Types::qubits_vector qubits = {qubit0, qubit1};
 
        if (GetSimulationType() == SimulationType::kStabilizer || GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "swap";
            op.qubits = AER::reg_t{qubit0, qubit1};
 
            state->buffer_op(std::move(op));
        }
        else
            // state->apply_mcswap(qubits);
            state->apply_unitary(qubits, AER::Linalg::Matrix::SWAP);
 
        NotifyObservers(qubits);
    }

    void ApplyCCX(Types::qubit_t qubit0, Types::qubit_t qubit1, Types::qubit_t qubit2) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CCX gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {qubit0, qubit1, qubit2};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            AER::Operations::Op op;
            op.type = AER::Operations::OpType::gate;
            op.name = "ccx";
            op.qubits = AER::reg_t{qubit0, qubit1, qubit2};
 
            state->buffer_op(std::move(op));
        }
        else
        {
            static const cmatrix_t mat = AER::Utils::make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}});
 
            // state->apply_mcx(qubits);
            state->apply_unitary(qubits, mat);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCSwap(Types::qubit_t ctrl_qubit, Types::qubit_t qubit0, Types::qubit_t qubit1) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CSwap gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, qubit0, qubit1};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            const int q1 = static_cast<int>(ctrl_qubit); // control
            const int q2 = static_cast<int>(qubit0);
            const int q3 = static_cast<int>(qubit1);
 
            ApplyCX(q3, q2);
            ApplyCSx(q2, q3);
            ApplyCX(q1, q2);
 
            ApplyP(q3, M_PI);
            ApplyP(q2, -M_PI_2);
 
            ApplyCSx(q2, q3);
            ApplyCX(q1, q2);
 
            ApplyP(q3, M_PI);
            ApplyCSx(q1, q3);
            ApplyCX(q3, q2);
        }
        else
        {
            static const cmatrix_t mat = AER::Utils::make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}, {0, 0}},
                                                                             {{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {1, 0}}});
 
            // state->apply_mcswap(qubits);
            state->apply_unitary(qubits, mat);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyCU(Types::qubit_t ctrl_qubit, Types::qubit_t tgt_qubit, double theta, double phi, double lambda, double gamma) override
    {
        if (GetSimulationType() == SimulationType::kStabilizer)
            throw std::runtime_error("CU gate not supported in stabilizer simulation");
 
        const Types::qubits_vector qubits = {ctrl_qubit, tgt_qubit};
 
        if (GetSimulationType() == SimulationType::kExtendedStabilizer)
        {
            if (gamma != 0.0)
                ApplyP(ctrl_qubit, gamma);
 
            const double lambdaPlusPhiHalf = 0.5 * (lambda + phi);
            const double halfTheta = 0.5 * theta;
            ApplyP(tgt_qubit, 0.5 * (lambda - phi));
            ApplyP(ctrl_qubit, lambdaPlusPhiHalf);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyU(tgt_qubit, -halfTheta, 0, -lambdaPlusPhiHalf, 0);
            ApplyCX(ctrl_qubit, tgt_qubit);
            ApplyU(tgt_qubit, halfTheta, phi, 0, 0);
        }
        else
        {
            const cmatrix_t CU = AER::Linalg::Matrix::cu(theta, phi, lambda, gamma);
            // state->apply_mcu(qubits, theta, phi, lambda, gamma);
            state->apply_unitary(qubits, CU);
        }
 
        NotifyObservers(qubits);
    }

    void ApplyNop() override
    {
        AER::Operations::Op op;
        op.type = AER::Operations::OpType::barrier; // there is also 'nop' but that one
                                                    // doesn't work with stabilizer
        state->buffer_op(std::move(op));
    }

    std::unique_ptr<ISimulator> Clone() override
    {
        auto sim = std::make_unique<AerSimulator>();
 
        sim->configuration = configuration; // copy the configuration
        for (const auto &[key, value] : configuration.GetConfigMap())
            sim->Configure(key.c_str(), value.c_str());
 
        // the tricky part is cloning the state
        if (simulationType == SimulationType::kMatrixProductState)
            sim->Configure("method", "matrix_product_state");
        else if (simulationType == SimulationType::kStabilizer)
            sim->Configure("method", "stabilizer");
        else if (simulationType == SimulationType::kTensorNetwork)
            sim->Configure("method", "tensor_network");
        else if (simulationType == SimulationType::kExtendedStabilizer)
            sim->Configure("method", "extended_stabilizer");
        else if (simulationType == SimulationType::kDensityMatrix)
            sim->Configure("method", "density_matrix");
        else
            sim->Configure("method", "statevector");
 
        if (simulationType == SimulationType::kMatrixProductState)
        {
            std::string key = "matrix_product_state_max_bond_dimension";
            if (configuration.IsSet(key))
                sim->Configure(key.c_str(), configuration.GetConfiguration(key).c_str());
 
            key = "matrix_product_state_truncation_threshold";
            if (configuration.IsSet(key))
                sim->Configure(key.c_str(), configuration.GetConfiguration(key).c_str());
 
            key = "mps_sample_measure_algorithm";
            if (configuration.IsSet(key))
                sim->Configure(key.c_str(), configuration.GetConfiguration(key).c_str());
        }
 
        sim->SetMultithreading(enableMultithreading); 
 
        AER::Vector<complex_t> localSavedAmplitudes = savedAmplitudes; 
        AER::Data localSavedState = savedState;                        
        AER::cmatrix_t localSavedDensityMatrix = savedDensityMatrix;   
        auto localSavedExtendedStabilizerState = savedExtendedStabilizerState;
 
        // now the tricky part
        if (state && state->is_initialized())
        {
            SaveState(); // this also restores the state if it's destroyed in the
                         // saving process!
            // now the current state is saved in savedAmplitudes or savedState, so use
            // it to initialize the cloned state, like this
            sim->savedAmplitudes = std::move(savedAmplitudes);
            sim->savedState = std::move(savedState);
            sim->savedDensityMatrix = std::move(savedDensityMatrix);
            sim->savedExtendedStabilizerState = savedExtendedStabilizerState;
 
            sim->RestoreState(); // now the state is loaded in the cloned simulator
 
            // those saved previously can be an older state, so put them in the clone,
            // too
            sim->savedAmplitudes = localSavedAmplitudes;
            sim->savedState = localSavedState;
            sim->savedDensityMatrix = localSavedDensityMatrix;
            sim->savedExtendedStabilizerState = localSavedExtendedStabilizerState;
 
            // those saved previously can be an older state, so put them back
            savedAmplitudes = std::move(localSavedAmplitudes);
            savedState = std::move(localSavedState);
            savedDensityMatrix = std::move(localSavedDensityMatrix);
            savedExtendedStabilizerState = std::move(localSavedExtendedStabilizerState);
        }
        else
        {
            // std::cout << "Restored from non initialized state" << std::endl;
            sim->savedAmplitudes = std::move(localSavedAmplitudes);
            sim->savedState = std::move(localSavedState);
            sim->savedDensityMatrix = std::move(localSavedDensityMatrix);
            sim->savedExtendedStabilizerState = std::move(localSavedExtendedStabilizerState);
 
            sim->RestoreState(); // this might not be necessary, but sometimes, not
                                 // very often, an exception is thrown about the
                                 // state not being initialized, this might prevent
                                 // that
        }
 
        if (configuration.IsSet("seed"))
            sim->SetSeed(DeriveSeed(std::stoull(configuration.GetConfiguration("seed")), nextSeedStream++));
 
        return sim;
    }
};
 
} // namespace Private
} // namespace Simulators
 
#endif
 
#endif
 
#endif