IntelSimulator.cpp

Go to the documentation of this file.

//##############################################################################
//##############################################################################

#include "Simulator.hpp"
#include "GateWriter.hpp"
#include "qureg/qureg.hpp"
#include "util/tinymatrix.hpp"
#include <cstdlib>
#include <iostream>
#include <limits>

#ifdef ENABLE_MPI
    #include "mpi.h"
#endif

namespace QNLP{

class IntelSimulator : public SimulatorGeneral<IntelSimulator> {
    public:
    using TMDP = qhipster::TinyMatrix<ComplexDP, 2, 2, 32>;
    using QRDP = QubitRegister<ComplexDP>;
    using CST = const std::size_t;

    IntelSimulator(int numQubits, bool useFusion=false) : SimulatorGeneral<IntelSimulator>(), 
                                    numQubits(numQubits), 
                                    qubitRegister(QubitRegister<ComplexDP> (numQubits, "base", 0)),
                                    gates(5), uid( reinterpret_cast<std::size_t>(this) ){

        //Define Pauli X
        gates[0](0,0) = ComplexDP(0.,0.);       gates[0](0,1) = ComplexDP(1.,0.);
        gates[0](1,0) = ComplexDP(1.,0.);       gates[0](1,1) = ComplexDP(0.,0.);

        //Define Pauli Y
        gates[1](0,0) = ComplexDP(0.,0.);       gates[1](0,1) = -ComplexDP(0.,1.);
        gates[1](1,0) = ComplexDP(0.,1.);       gates[1](1,1) = ComplexDP(0.,0.);

        //Define Pauli Z
        gates[2](0,0) = ComplexDP(1.,0.);       gates[2](0,1) = ComplexDP(0.,0.);
        gates[2](1,0) = ComplexDP(0.,0.);       gates[2](1,1) = ComplexDP(-1.,0.);

        //Define I
        gates[3](0,0) = ComplexDP(1.,0.);       gates[3](0,1) = ComplexDP(0.,0.);
        gates[3](1,0) = ComplexDP(0.,0.);       gates[3](1,1) = ComplexDP(1.,0.);

        //Define Pauli H
        double coeff = (1./sqrt(2.));
        gates[4](0,0) = coeff*ComplexDP(1.,0.);   gates[4](0,1) = coeff*ComplexDP(1.,0.);
        gates[4](1,0) = coeff*ComplexDP(1.,0.);   gates[4](1,1) = -coeff*ComplexDP(1.,0.);

        //Ensure the cache maps are populated before use.
        this->initCaches();

        #ifdef ENABLE_MPI //If for some strange reason the MPI environement is not enable through the Base CRTP class, enable using Intel-QS
            int mpi_is_init;
            MPI_Initialized(&mpi_is_init);
            if (! mpi_is_init){ // Attempt init using Intel-QS MPI env
                int argc_tmp = 0;
                char** argv_tmp = new char*[argc_tmp];

                qhipster::mpi::Environment env(argc_tmp, argv_tmp); //we do not expect to pass any params here
                delete argv_tmp;
            }
            MPI_Comm_rank(MPI_COMM_WORLD, &rank);
        #endif

        /* Set up random number generator for randomly collapsing qubit to 0 or 1
         *
         * Note: a random number will be generated on rank 0 and then broadcasted
         * to all ranks so that each rank collapses the respective qubit to the
         * same value.
         */
        std::mt19937 mt_(rd()); 
        std::uniform_real_distribution<double> dist_(0.0,1.0);
        mt = mt_;
        dist = dist_;
        if(useFusion == true){
            qubitRegister.TurnOnFusion();
            std::cerr << "Warning: enabling fusion may cause inconsistent results." << std::endl;
        }
        gate_count_1qubit = 0;
        gate_count_2qubit = 0;
    }

    ~IntelSimulator(){ }

    // 1 qubit
    inline void applyGateU(const TMDP& U, CST qubitIndex, std::string label="U"){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.Apply1QubitGate(qubitIndex, U);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall(label, U.tostr(), qubitIndex);
        #endif
    }

    inline void applyGateI(std::size_t qubitIndex){
        #ifndef RESOURCE_ESTIMATE
        applyGateU(getGateI(), qubitIndex, "I");
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("I", getGateI().tostr(), qubitIndex);
        #endif
    }

    inline void applyGatePhaseShift(std::size_t qubit_idx, double angle){
        //Phase gate is identity with 1,1 index modulated by angle
        TMDP U(gates[3]);
        U(1, 1) = ComplexDP(cos(angle), sin(angle));

        #ifndef RESOURCE_ESTIMATE
        applyGateU(U, qubit_idx, "Phase:=" + std::to_string(angle));
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("PShift(theta=" + std::to_string(angle) + ")", U.tostr(), qubitIndex);
        #endif

    }

    // 2 qubit
    inline void applyGateSqrtSwap(  std::size_t qubit_idx0, std::size_t qubit_idx1){    
        std::cerr << "NOT YET IMPLEMENTED" << std::endl; 
        std::abort();
    }

    // 3 qubit
    inline void applyGateCCX(std::size_t ctrl_qubit0, std::size_t ctrl_qubit1, std::size_t target_qubit){
        this->applyGateNCU(this->getGateX(), std::vector<std::size_t> {ctrl_qubit0, ctrl_qubit1}, target_qubit, "X");
    }

    inline void applyGateCSwap(std::size_t ctrl_qubit, std::size_t qubit_swap0, std::size_t qubit_swap1){
        //V = sqrt(X)
        TMDP V;
        V(0,0) = {0.5,  0.5};
        V(0,1) = {0.5, -0.5};
        V(1,0) = {0.5, -0.5};
        V(1,1) = {0.5,  0.5};
        
        TMDP V_dag;
        V_dag(0,0) = {0.5, -0.5};
        V_dag(0,1) = {0.5,  0.5};
        V_dag(1,0) = {0.5,  0.5};
        V_dag(1,1) = {0.5, -0.5};
        
        applyGateCX(qubit_swap1, qubit_swap0);
        applyGateCU(V, qubit_swap0, qubit_swap1, "X");
        applyGateCU(V, ctrl_qubit, qubit_swap1, "X");
        
        applyGateCX(ctrl_qubit, qubit_swap0);
        applyGateCU(V_dag, qubit_swap0, qubit_swap1, "X");
        applyGateCX(qubit_swap1, qubit_swap0);
        applyGateCX(ctrl_qubit, qubit_swap0);
    }

    //#################################################

    inline void applyGateX(CST qubitIndex){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyPauliX(qubitIndex);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("X", getGateX().tostr(), qubitIndex);
        #endif
    }

    inline void applyGateY(CST qubitIndex){ 
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyPauliY(qubitIndex);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("Y", getGateY().tostr(), qubitIndex);
        #endif
    }

    inline void applyGateZ(CST qubitIndex){ 
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyPauliZ(qubitIndex);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("Z", getGateZ().tostr(), qubitIndex);
        #endif
    }

    inline void applyGateH(CST qubitIndex){ 
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyHadamard(qubitIndex);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall("H", getGateH().tostr(), qubitIndex);
        #endif
    }

   inline void applyGateSqrtX(CST qubitIndex){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyPauliSqrtX(qubitIndex);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall(
            "\\sqrt[2]{X}", 
            matrixSqrt<decltype(getGateX())>(getGateX()).tostr(), 
            qubitIndex
        );
        #endif
    };

    inline void applyGateRotX(CST qubitIndex, double angle) {
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyRotationX(qubitIndex, angle);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall(
            "R_X(\\theta=" + std::to_string(angle) + ")", 
            getGateI().tostr(), 
            qubitIndex
        );
        #endif
    };

    inline void applyGateRotY(CST qubitIndex, double angle) {
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyRotationY(qubitIndex, angle);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall(
            "R_Y(\\theta=" + std::to_string(angle) + ")", 
            getGateI().tostr(), 
            qubitIndex
        );
        #endif
    };

    inline void applyGateRotZ(CST qubitIndex, double angle) {
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyRotationZ(qubitIndex, angle);
        #endif

        gate_count_1qubit++;

        #ifdef GATE_LOGGING
        writer.oneQubitGateCall(
            "R_Z(\\theta=" + std::to_string(angle) + ")", 
            getGateI().tostr(), 
            qubitIndex
        );
        #endif
    };

    inline TMDP getGateX(){ return gates[0]; }

    inline TMDP getGateY(){ return gates[1]; }

    inline TMDP getGateZ(){ return gates[2]; }

    inline TMDP getGateI(){ return gates[3]; }

    inline TMDP getGateH(){ return gates[4]; }

    inline void applyGateCU(const TMDP& U, CST control, CST target, std::string label="U"){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyControlled1QubitGate(control, target, U);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( label, U.tostr(), control, target );
        #endif
    }

    inline void applyGateCX(CST control, CST target){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCPauliX(control, target);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "X", getGateX().tostr(), control, target );
        #endif
    }

    inline void applyGateCY(CST control, CST target){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCPauliY(control, target);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "Y", getGateY().tostr(), control, target );
        #endif
    }

    inline void applyGateCZ(CST control, CST target){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCPauliZ(control, target);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "Z", getGateZ().tostr(), control, target );
        #endif
    }

    inline void applyGateCH(CST control, CST target){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCHadamard(control, target);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "H", getGateH().tostr(), control, target );
        #endif
    }

    inline void applyGateCPhaseShift(double angle, unsigned int control, unsigned int target){
        TMDP U(gates[3]);
        U(1, 1) = ComplexDP(cos(angle), sin(angle));

        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyControlled1QubitGate(control, target, U);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "CPhase", U.tostr(), control, target );
        #endif
    }

    inline void applyGateCRotX(CST control, CST target, const double theta){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCRotationX(control, target, theta);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "CR_X", getGateI().tostr(), control, target );
        #endif
    }

    inline void applyGateCRotY(CST control, CST target, double theta){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCRotationY(control, target, theta);
        #endif

        gate_count_2qubit++;

        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "CR_Y", getGateI().tostr(), control, target );
        #endif
    }

    inline void applyGateCRotZ(CST control, CST target, const double theta){
        #ifndef RESOURCE_ESTIMATE
        qubitRegister.ApplyCRotationZ(control, target, theta);
        #endif
        
        gate_count_2qubit++;
        
        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "CR_Z", getGateI().tostr(), control, target );
        #endif
    }

    inline void applyGateSwap(CST qubit_idx0, CST qubit_idx1){
        qubitRegister.ApplySwap(qubit_idx0, qubit_idx1);
        #ifdef GATE_LOGGING
        writer.twoQubitGateCall( "SWAP", getGateI().tostr(), qubit_idx0, qubit_idx1 );
        #endif
    }

    inline QubitRegister<ComplexDP>& getQubitRegister() { 
        return this->qubitRegister; 
    }

    inline const QubitRegister<ComplexDP>& getQubitRegister() const { 
        return this->qubitRegister; 
    };

    std::size_t getNumQubits() { 
        return numQubits; 
    }

    void initRegister(){
        this->qubitRegister.Initialize("base",0);
        this->initCaches();
        gate_count_1qubit = 0;
        gate_count_2qubit = 0;
    }

    inline void applyAmplitudeNorm(){
        this->qubitRegister.Normalize();
    }

    bool applyMeasurement(CST target, bool normalize=true){
        double rand;
        bool bit_val;

        #ifdef ENABLE_MPI
            if(rank == 0){
                rand = dist(mt);
            }
            MPI_Bcast(&rand, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Barrier(MPI_COMM_WORLD);
        #else
            rand = dist(mt);
        #endif
        collapseQubit(target, (bit_val = ( rand < getStateProbability(target) ) ) );
        if(normalize){
            applyAmplitudeNorm();
        }
        return bit_val;
    }

    void collapseToBasisZ(CST target, bool collapseValue){
        collapseQubit(target, collapseValue);
        applyAmplitudeNorm();
    }

    inline void PrintStates(std::string x, std::vector<std::size_t> qubits = {}){
        qubitRegister.Print(x,qubits);
    }

    #ifdef GATE_LOGGING

    GateWriter& getGateWriter(){
        return writer;
    } 
    #endif

    std::pair<std::size_t, std::size_t> getGateCounts(){
        std::cout << "######### Gate counts #########" << std::endl;
        std::cout << "1 qubit = " << gate_count_1qubit << std::endl;
        std::cout << "2 qubit = " << gate_count_2qubit << std::endl;
        std::cout << "total = " << gate_count_1qubit + gate_count_2qubit << std::endl;
        std::cout << "###############################" << std::endl;
        return std::make_pair(gate_count_1qubit, gate_count_2qubit);
    }

    inline complex<double> overlap( IntelSimulator &sim){
        if(sim.uid != this->uid){
            return qubitRegister.ComputeOverlap(sim.qubitRegister);
        }
        else{
            return std::numeric_limits<double>::quiet_NaN();
        }
    }

    private:
    //inline static std::size_t suid = 0; //works in C++17.
    const std::size_t uid;

    std::size_t numQubits = 0;
    QRDP qubitRegister;
    std::vector<TMDP> gates;
    #ifdef ENABLE_MPI
        int rank;
    #endif

    std::size_t gate_count_1qubit;
    std::size_t gate_count_2qubit;

    std::random_device rd; 
    std::mt19937 mt;
    std::uniform_real_distribution<double> dist;


    // Measurement methods
    inline void collapseQubit(CST target, bool collapseValue){
        qubitRegister.CollapseQubit(target, collapseValue);
    }

    inline double getStateProbability(CST target){
        return qubitRegister.GetProbability(target);
    }

};

};