PathIntegralSimulator.h

Go to the documentation of this file.

 
#pragma once
 
#ifndef _PATH_INTEGRAL_SIMULATOR_H_
#define _PATH_INTEGRAL_SIMULATOR_H_
 
#include <optional>
#include <type_traits>
#include <utility>
 
#include "../../Circuit/Circuit.h"
#include "PathIntegral.h"
 
namespace Simulators
{
 
class PathIntegralSimulator
{
  public:
    using AmplitudeMap = std::remove_reference_t<decltype(std::declval<QC::PathIntegral::PathIntegralSimulator &>().GetAmplitudes())>;
 
    void SetSeed(uint64_t seed)
    {
        simulator.SetSeed(seed);
    }
 
    void SetTrimValue(double val)
    {
        simulator.SetTrimValue(val);
    }
 
    double GetTrimValue() const
    {
        return simulator.GetTrimValue();
    }
 
    void SetMaxDoublingsForBackwardPaths(size_t doublings)
    {
        simulator.SetMaxDoublingsForBackwardPaths(doublings);
    }
 
    size_t GetMaxDoublingsForBackwardPaths() const
    {
        return simulator.GetMaxDoublingsForBackwardPaths();
    }
 
    void Reset()
    {
        simulator.Reset();
    }
 
    bool SetCircuit(const std::shared_ptr<Circuits::Circuit<>> &circuit)
    {
        const auto convertedCircuit = ConvertCircuit(circuit);
        if (convertedCircuit.empty())
            return false;
        simulator.SetCircuit(convertedCircuit);
        return true;
    }
 
    std::complex<double> AmplitudeFromZero(const std::vector<bool> &endState)
    {
        return simulator.Propagate(endState);
    }
 
    std::complex<double> Amplitude(const std::vector<bool> &startState, const std::vector<bool> &endState)
    {
        return simulator.Propagate(startState, endState);
    }
 
    std::optional<QC::Gates::AppliedGate<>> ConvertGate(const std::shared_ptr<Circuits::IQuantumGate<>> &gate)
    {
        if (!gate)
            return std::nullopt;
 
        switch (gate->GetGateType())
        {
            // single qubit gates
        case Circuits::QuantumGateType::kPhaseGateType: {
            const auto g = std::static_pointer_cast<Circuits::PhaseGate<>>(gate);
            pgate.SetPhaseShift(g->GetLambda());
            return QC::Gates::AppliedGate<>(pgate.getRawOperatorMatrix(), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kXGateType:
            return QC::Gates::AppliedGate<>(xgate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kYGateType:
            return QC::Gates::AppliedGate<>(ygate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kZGateType:
            return QC::Gates::AppliedGate<>(zgate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kHadamardGateType:
            return QC::Gates::AppliedGate<>(h.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kSGateType:
            return QC::Gates::AppliedGate<>(sgate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kSdgGateType:
            return QC::Gates::AppliedGate<>(sdggate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kTGateType:
            return QC::Gates::AppliedGate<>(tgate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kTdgGateType:
            return QC::Gates::AppliedGate<>(tdggate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kSxGateType:
            return QC::Gates::AppliedGate<>(sxgate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kSxDagGateType:
            return QC::Gates::AppliedGate<>(sxdaggate.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kKGateType:
            return QC::Gates::AppliedGate<>(k.getRawOperatorMatrix(), gate->GetQubit(0));
        case Circuits::QuantumGateType::kRxGateType: {
            const auto g = std::static_pointer_cast<Circuits::RxGate<>>(gate);
            rxgate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(rxgate.getRawOperatorMatrix(), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kRyGateType: {
            const auto g = std::static_pointer_cast<Circuits::RyGate<>>(gate);
            rygate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(rygate.getRawOperatorMatrix(), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kRzGateType: {
            const auto g = std::static_pointer_cast<Circuits::RzGate<>>(gate);
            rzgate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(rzgate.getRawOperatorMatrix(), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kUGateType: {
            const auto g = std::static_pointer_cast<Circuits::UGate<>>(gate);
            ugate.SetParams(g->GetTheta(), g->GetPhi(), g->GetLambda(), g->GetGamma());
            return QC::Gates::AppliedGate<>(ugate.getRawOperatorMatrix(), gate->GetQubit(0));
        }
 
        // two qubit gates
        case Circuits::QuantumGateType::kSwapGateType:
            return QC::Gates::AppliedGate<>(swapgate.getRawOperatorMatrix(), gate->GetQubit(0), gate->GetQubit(1));
        case Circuits::QuantumGateType::kCXGateType:
            return QC::Gates::AppliedGate<>(cxgate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCYGateType:
            return QC::Gates::AppliedGate<>(cygate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCZGateType:
            return QC::Gates::AppliedGate<>(czgate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCPGateType: {
            const auto g = std::static_pointer_cast<Circuits::CPGate<>>(gate);
            cpgate.SetPhaseShift(g->GetLambda());
            return QC::Gates::AppliedGate<>(cpgate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kCRxGateType: {
            const auto g = std::static_pointer_cast<Circuits::CRxGate<>>(gate);
            crxgate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(crxgate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kCRyGateType: {
            const auto g = std::static_pointer_cast<Circuits::CRyGate<>>(gate);
            crygate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(crygate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kCRzGateType: {
            const auto g = std::static_pointer_cast<Circuits::CRzGate<>>(gate);
            crzgate.SetTheta(g->GetTheta());
            return QC::Gates::AppliedGate<>(crzgate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        }
        case Circuits::QuantumGateType::kCHGateType:
            return QC::Gates::AppliedGate<>(ch.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCSxGateType:
            return QC::Gates::AppliedGate<>(csx.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCSxDagGateType:
            return QC::Gates::AppliedGate<>(csxdag.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCUGateType: {
            const auto g = std::static_pointer_cast<Circuits::CUGate<>>(gate);
            cugate.SetParams(g->GetTheta(), g->GetPhi(), g->GetLambda(), g->GetGamma());
            return QC::Gates::AppliedGate<>(cugate.getRawOperatorMatrix(), gate->GetQubit(1), gate->GetQubit(0));
        }
 
        // three qubit gates
        case Circuits::QuantumGateType::kCCXGateType:
            return QC::Gates::AppliedGate<>(ccxgate.getRawOperatorMatrix(), gate->GetQubit(2), gate->GetQubit(1), gate->GetQubit(0));
        case Circuits::QuantumGateType::kCSwapGateType:
            return QC::Gates::AppliedGate<>(cswapgate.getRawOperatorMatrix(), gate->GetQubit(2), gate->GetQubit(1), gate->GetQubit(0));
 
        default:
            return std::nullopt;
        }
    }
 
    std::vector<QC::Gates::AppliedGate<>> ConvertCircuit(const std::shared_ptr<Circuits::Circuit<>> &circuit)
    {
        std::vector<QC::Gates::AppliedGate<>> result;
        if (!circuit)
            return result;
 
        const auto &operations = circuit->GetOperations();
        result.reserve(operations.size());
 
        for (const auto &op : operations)
        {
            if (!op || op->GetType() != Circuits::OperationType::kGate)
                continue;
 
            const auto &gate = std::static_pointer_cast<Circuits::IQuantumGate<>>(op);
            if (auto applied = ConvertGate(gate))
                result.emplace_back(std::move(*applied));
        }
 
        return result;
    }
 
    void PropagateStep(const QC::Gates::AppliedGate<> &gate, AmplitudeMap &currentAmplitudes)
    {
        simulator.PropagateStep(gate, currentAmplitudes);
    }
 
    static size_t GetBranchingForQcsimCircuit(const std::vector<QC::Gates::AppliedGate<>> &circuit)
    {
        size_t doublings = 0;
        for (const auto &gate : circuit)
            if (gate.isBranching())
                ++doublings;
 
        return doublings;
    }
 
    size_t GetBranchingForMaestroCircuit(const std::shared_ptr<Circuits::Circuit<>> &circuit)
    {
        /*
        const auto circ = ConvertCircuit(circuit);
        return GetBranchingForQcsimCircuit(circ);
        */
 
        return GetBranchingForMaestroCircuitFromPos(circuit, 0);
    }
 
    size_t GetBranchingForMaestroCircuitFromPos(const std::shared_ptr<Circuits::Circuit<>> &circuit, size_t pos)
    {
        if (!circuit)
            return 0;
 
        size_t doublings = 0;
        const auto &operations = circuit->GetOperations();
        for (size_t i = pos; i < operations.size(); ++i)
            if (operations[i]->IsBranching())
                ++doublings;
 
        return doublings;
    }
 
    double QubitProbability(size_t qubit, bool value = true)
    {
        return simulator.QubitProbability(qubit, value);
    }
 
    std::complex<double> AmplitudeForOutcome(size_t outcome)
    {
        const auto &amplitudes = simulator.GetAmplitudes();
        if (amplitudes.empty())
            return 0.0;
 
        QC::PathIntegral::FastVectorBool state(amplitudes.begin()->first.size());
        for (size_t q = 0; q < state.size(); ++q)
            state.set(q, (outcome >> q) & 1);
        if (auto it = amplitudes.find(state); it != amplitudes.end())
            return it->second;
        return 0.0;
    }
 
    double Probability(size_t outcome)
    {
        return std::norm(AmplitudeForOutcome(outcome));
    }
 
    QC::PathIntegral::FastVectorBool MeasureNoCollapse()
    {
        return simulator.MeasureNoCollapse();
    }
 
    bool MeasureQubit(size_t qubit)
    {
        return simulator.MeasureQubit(qubit);
    }
 
    std::unique_ptr<PathIntegralSimulator> Clone() const
    {
        auto clone = std::make_unique<PathIntegralSimulator>();
        clone->simulator = simulator.Clone();
        return clone;
    }
 
    AmplitudeMap &Amplitudes()
    {
        return simulator.GetAmplitudes();
    }
 
    void SetStartZeroState(size_t numQubits)
    {
        auto &amplitudes = simulator.GetAmplitudes();
 
        amplitudes.clear();
        QC::PathIntegral::FastVectorBool zeroState(numQubits);
        amplitudes[zeroState] = 1.0;
    }
 
    double ExpectationValue(const std::string &pauliStringOrig)
    {
        // save current amplitudes
        const auto &currentAmplitudes = simulator.GetAmplitudes();
        auto amplitudes = currentAmplitudes;
 
        // now apply the pauli string as gates
        for (size_t i = 0; i < pauliStringOrig.size(); ++i)
        {
            const char c = pauliStringOrig[i];
            switch (c)
            {
            case 'X':
                [[fallthrough]];
            case 'x':
                PropagateStep(QC::Gates::AppliedGate<>(xgate.getRawOperatorMatrix(), i), amplitudes);
                break;
            case 'Y':
                [[fallthrough]];
            case 'y':
                PropagateStep(QC::Gates::AppliedGate<>(ygate.getRawOperatorMatrix(), i), amplitudes);
                break;
            case 'Z':
                [[fallthrough]];
            case 'z':
                PropagateStep(QC::Gates::AppliedGate<>(zgate.getRawOperatorMatrix(), i), amplitudes);
                break;
            default:
                break;
            }
        }
 
        // use the new amplitudes to calculate the expectation value
        double expectation = 0.0;
 
        for (const auto &[state, amplitude] : amplitudes)
        {
            if (currentAmplitudes.find(state) == currentAmplitudes.end())
                continue;
 
            const auto &oldAmplitude = currentAmplitudes.at(state);
            expectation += std::real(std::conj(oldAmplitude) * amplitude);
        }
 
        return expectation;
    }
 
    void ApplyGate(const QC::Gates::AppliedGate<> &gate)
    {
        simulator.PropagateStep(gate, simulator.GetAmplitudes());
    }
 
    void SaveState()
    {
        simulator.SaveAmplitudes();
    }
 
    void RestoreState()
    {
        simulator.RestoreAmplitudes();
    }
 
  private:
    QC::PathIntegral::PathIntegralSimulator simulator;
 
    // These gates are never modified, and their matrix accessors return const
    // references. Sharing them avoids allocation on construction and cloning.
    inline static const QC::Gates::PauliXGate<> xgate;
    inline static const QC::Gates::PauliYGate<> ygate;
    inline static const QC::Gates::PauliZGate<> zgate;
    inline static const QC::Gates::HadamardGate<> h;
    inline static const QC::Gates::SGate<> sgate;
    inline static const QC::Gates::SDGGate<> sdggate;
    inline static const QC::Gates::TGate<> tgate;
    inline static const QC::Gates::TDGGate<> tdggate;
    inline static const QC::Gates::SquareRootNOTGate<> sxgate;
    inline static const QC::Gates::SquareRootNOTDagGate<> sxdaggate;
    inline static const QC::Gates::HyGate<> k;
    inline static const QC::Gates::CNOTGate<> cxgate;
    inline static const QC::Gates::ControlledYGate<> cygate;
    inline static const QC::Gates::ControlledZGate<> czgate;
    inline static const QC::Gates::ControlledHadamardGate<> ch;
    inline static const QC::Gates::ControlledSquareRootNOTGate<> csx;
    inline static const QC::Gates::ControlledSquareRootNOTDagGate<> csxdag;
    inline static const QC::Gates::SwapGate<> swapgate;
    inline static const QC::Gates::ToffoliGate<> ccxgate;
    inline static const QC::Gates::FredkinGate<> cswapgate;
 
    // Setters update these matrices before each application. Keep them local
    // to each simulator so independent instances and clones cannot interfere.
    QC::Gates::PhaseShiftGate<> pgate;
    QC::Gates::RxGate<> rxgate;
    QC::Gates::RyGate<> rygate;
    QC::Gates::RzGate<> rzgate;
    QC::Gates::UGate<> ugate;
    QC::Gates::ControlledPhaseShiftGate<> cpgate;
    QC::Gates::ControlledRxGate<> crxgate;
    QC::Gates::ControlledRyGate<> crygate;
    QC::Gates::ControlledRzGate<> crzgate;
    QC::Gates::ControlledUGate<> cugate;
};
 
} // namespace Simulators
 
#endif // _PATH_INTEGRAL_SIMULATOR_H_