Simulator.hpp

Go to the documentation of this file.

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

#ifndef QNLP_SIMULATOR_H
#define QNLP_SIMULATOR_H
#include <cstddef>
#include <utility> //std::declval
#include <vector>
#include <iostream>

// Include all additional modules to be used within simulator
#include "GateWriter.hpp"
#include "ncu.hpp"
#include "oracle.hpp"
#include "diffusion.hpp"
#include "qft.hpp"
#include "arithmetic.hpp"
#include "bin_into_superpos.hpp"
#include "hamming.hpp"
#include "bit_group.hpp"

#if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
//pybind11 fails to compile with icpc using cpp17 support, so fallback to 14 if necessary
#include <experimental/any>
#else
#include <any>
#endif

#ifdef VIRTUAL_INTERFACE
#include "ISimulator.hpp"
#endif

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

namespace QNLP{

    /*
     * @brief Returns the boolean True if the bit at index 'bit' in 'byte' is set to 1, else returns False
     * 
     * @param byte An integer whose bit at index 'bit' is being tested if it is set
     * @param bit The index into the integer 'byte' being tested if it is set
     */
    #define IS_SET(byte,bit) (((byte) & (1UL << (bit))) >> (bit))

    template <class DerivedType>
    #ifdef VIRTUAL_INTERFACE
    class SimulatorGeneral : virtual public ISimulator {
    #else
    class SimulatorGeneral {
    #endif
    
    private:
    SimulatorGeneral(){ 
    //If we are building MPI support, ensure that it is init'd here before subclasses.
    #ifdef ENABLE_MPI
        int mpi_is_init;
        MPI_Initialized(&mpi_is_init);
        if (! mpi_is_init){
            int argc_tmp = 0;
            char** argv_tmp = new char*[argc_tmp];
            MPI_Init(&argc_tmp, &argv_tmp);
            delete argv_tmp;
        }
    #endif

        sim_ncu =  NCU<DerivedType>(static_cast<DerivedType&>(*this));
    }; 
    
    friend DerivedType;
    
    protected:
    #ifdef GATE_LOGGING
    GateWriter writer;
    #endif


    #if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
    std::experimental::any sim_ncu;
    #else
    std::any sim_ncu;
    #endif

    public:
        //using Mat2x2Type = decltype(std::declval<DerivedType>().getGateX());
        virtual ~SimulatorGeneral(){ }
        //##############################################################################
        //                           Single qubit gates
        //##############################################################################

        void applyGateX(std::size_t qubit_idx){ 
            static_cast<DerivedType&>(*this).applyGateX(qubit_idx);
        };
        void applyGateY(std::size_t qubit_idx){
            static_cast<DerivedType&>(*this).applyGateY(qubit_idx);
        }
        void applyGateZ(std::size_t qubit_idx){
            static_cast<DerivedType&>(*this).applyGateZ(qubit_idx);
        }
        void applyGateI(std::size_t qubit_idx){
            static_cast<DerivedType&>(*this).applyGateI(qubit_idx); 
        }
        void applyGateH(std::size_t qubit_idx){
            static_cast<DerivedType&>(*this).applyGateH(qubit_idx);
        }
        void applyGateSqrtX(std::size_t qubit_idx){
            static_cast<DerivedType&>(*this).applyGateSqrtX(qubit_idx);
        }
        void applyGateRotX(std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateRotX(qubit_idx, angle_rad);
        }
        void applyGateRotY(std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateRotY(qubit_idx, angle_rad);
        }
        void applyGateRotZ(std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateRotZ(qubit_idx, angle_rad);
        }

        template<class Mat2x2Type>
        void applyGateU(const Mat2x2Type &U, std::size_t qubit_idx){
            static_cast<DerivedType*>(this)->applyGateU(U, qubit_idx);
        }

        decltype(auto) getGateX() {
            return static_cast<DerivedType&>(*this).getGateX();
        }

        decltype(auto) getGateY(){
            return static_cast<DerivedType&>(*this).getGateY();
        }
        decltype(auto) getGateZ(){
            return static_cast<DerivedType&>(*this).getGateZ();
        }
        decltype(auto) getGateI(){
            return static_cast<DerivedType*>(this)->getGateI();
        }
        decltype(auto) getGateH(){
            return static_cast<DerivedType*>(this)->getGateH();
        }

        //##############################################################################
        //                                  2 qubit gates
        //##############################################################################

        template<class Mat2x2Type>
        void applyGateCU(const Mat2x2Type &U, const std::size_t control, const std::size_t target, std::string label="CU"){
            static_cast<DerivedType*>(this)->applyGateCU(U, control, target, label);
        }

        void applyGateCX(const std::size_t control, const std::size_t target){
            static_cast<DerivedType*>(this)->applyGateCX(control, target);
        }

        void applyGateCY(const std::size_t control, const std::size_t target){
            static_cast<DerivedType*>(this)->applyGateCY(control, target);
        }

        void applyGateCZ(const std::size_t control, const std::size_t target){
            static_cast<DerivedType*>(this)->applyGateCZ(control, target);
        }

        void applyGateCH(const std::size_t control, const std::size_t target){
            static_cast<DerivedType*>(this)->applyGateCH(control, target);
        }

        void applyGateSwap(std::size_t qubit_idx0, std::size_t qubit_idx1){
            static_cast<DerivedType*>(this)->applyGateSwap(qubit_idx0,qubit_idx1);
        }

        void applyGateSqrtSwap(std::size_t qubit_idx0, std::size_t qubit_idx1){
            static_cast<DerivedType*>(this)->applyGateSqrtSwap(qubit_idx0,qubit_idx1);
        }
        
        void applyGatePhaseShift(double angle, std::size_t qubit_idx){
            static_cast<DerivedType*>(this)->applyGatePhaseShift(angle, qubit_idx);
        }

        void applyGateCPhaseShift(double angle, std::size_t control, std::size_t target){
            static_cast<DerivedType*>(this)->applyGateCPhaseShift(angle, control, target);
        }

        void applyGateCCX(std::size_t ctrl_qubit0, std::size_t ctrl_qubit1, std::size_t target_qubit){
            static_cast<DerivedType*>(this)->applyGateCCX(ctrl_qubit0, ctrl_qubit1, target_qubit);
        }

        void applyGateCSwap(std::size_t ctrl_qubit, std::size_t qubit_swap0, std::size_t qubit_swap1){
            assert( static_cast<DerivedType*>(this)->getNumQubits() > 2 );
            //The requested qubit indices must be available to use within the register range
            assert( (ctrl_qubit < getNumQubits() ) && (qubit_swap0 < getNumQubits() ) && (qubit_swap1 < getNumQubits()) );
            //The qubits must be different from one another
            assert( ctrl_qubit != qubit_swap0 && ctrl_qubit != qubit_swap1 && qubit_swap0 != qubit_swap1 );
            static_cast<DerivedType*>(this)->applyGateCSwap(ctrl_qubit, qubit_swap0, qubit_swap1);
        }

        void applyGateCRotX(std::size_t ctrl_qubit, std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateCRotX(ctrl_qubit, qubit_idx, angle_rad);
        }

        void applyGateCRotY(std::size_t ctrl_qubit, std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateCRotY(ctrl_qubit, qubit_idx, angle_rad);
        }

        void applyGateCRotZ(std::size_t ctrl_qubit, std::size_t qubit_idx, double angle_rad){
            static_cast<DerivedType&>(*this).applyGateCRotZ(ctrl_qubit, qubit_idx, angle_rad);
        }

        decltype(auto) getQubitRegister(){ 
           return static_cast<DerivedType*>(this)->getQubitRegister();
        }

        std::size_t getNumQubits(){
            return static_cast<DerivedType*>(this)->getNumQubits();    
        }

        void applyQFT(std::size_t minIdx, std::size_t maxIdx){
            QFT<decltype(static_cast<DerivedType&>(*this))>::applyQFT(static_cast<DerivedType&>(*this), minIdx, maxIdx);
        }

        void applyIQFT(std::size_t minIdx, std::size_t maxIdx){
            QFT<decltype(static_cast<DerivedType&>(*this))>::applyIQFT(static_cast<DerivedType&>(*this), minIdx, maxIdx);
        }


        void sumReg(std::size_t r0_minIdx, std::size_t r0_maxIdx, std::size_t r1_minIdx, std::size_t r1_maxIdx){
            Arithmetic<decltype(static_cast<DerivedType&>(*this))>::sum_reg(static_cast<DerivedType&>(*this), r0_minIdx, r0_maxIdx, r1_minIdx, r1_maxIdx);
        }

        void subReg(std::size_t r0_minIdx, std::size_t r0_maxIdx, std::size_t r1_minIdx, std::size_t r1_maxIdx){
            Arithmetic<decltype(static_cast<DerivedType&>(*this))>::sub_reg(static_cast<DerivedType&>(*this), r0_minIdx, r0_maxIdx, r1_minIdx, r1_maxIdx);
        }

        template<class Mat2x2Type>
        void applyGateNCU(const Mat2x2Type& U, const std::vector<std::size_t>& ctrlIndices, std::size_t target, std::string label){
            #if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
            std::experimental::any_cast<NCU<DerivedType>&>(sim_ncu).applyNQubitControl(static_cast<DerivedType&>(*this), ctrlIndices, {}, target, label, U, 0);
            #else
            std::any_cast<NCU<DerivedType>&>(sim_ncu).applyNQubitControl(static_cast<DerivedType&>(*this), ctrlIndices, {}, target, label, U, 0);
            #endif
        }

        template<class Mat2x2Type>
        void applyGateNCU(const Mat2x2Type& U, const std::vector<std::size_t>& ctrlIndices, const std::vector<std::size_t>& auxIndices, std::size_t target, std::string label){
            #if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
            std::experimental::any_cast<NCU<DerivedType>&>(sim_ncu).applyNQubitControl(static_cast<DerivedType&>(*this), ctrlIndices, auxIndices, target, label, U, 0);
            #else
            std::any_cast<NCU<DerivedType>&>(sim_ncu).applyNQubitControl(static_cast<DerivedType&>(*this), ctrlIndices, auxIndices, target, label, U, 0);
            #endif
        }

        template<class Mat2x2Type>
        void applyOracleU(std::size_t bit_pattern, const std::vector<std::size_t>& ctrlIndices, std::size_t target, const Mat2x2Type& U , std::string gateLabel){
            Oracle<DerivedType>::bitStringNCU(static_cast<DerivedType&>(*this), bit_pattern, ctrlIndices, target, U, gateLabel);
        }

        template<class Mat2x2Type>
        void applyOracleU(std::size_t bit_pattern, const std::vector<std::size_t>& ctrlIndices, const std::vector<std::size_t>& auxIndices, std::size_t target, const Mat2x2Type& U , std::string gateLabel){
            Oracle<DerivedType>::bitStringNCU(static_cast<DerivedType&>(*this), bit_pattern, ctrlIndices, auxIndices, target, U, gateLabel);
        }

        void applyOraclePhase(std::size_t bit_pattern, const std::vector<std::size_t>& ctrlIndices, std::size_t target){
            //applyOracleU<decltype(static_cast<DerivedType*>(this)->getGateZ())>(bit_pattern, ctrlIndices, target, static_cast<DerivedType*>(this)->getGateZ());
            Oracle<DerivedType> oracle;
            oracle.bitStringPhaseOracle(static_cast<DerivedType&>(*this), bit_pattern, ctrlIndices, target );
        }

        void applyDiffusion(const std::vector<std::size_t>& ctrlIndices, std::size_t target){
            Diffusion<DerivedType> diffusion;
            diffusion.applyOpDiffusion(static_cast<DerivedType&>(*this), ctrlIndices, target);
        }

        void encodeToRegister(std::size_t target_pattern, 
                const std::vector<std::size_t> target_register, 
                std::size_t len_bin_pattern){

            for(std::size_t i = 0; i < len_bin_pattern; i++){
                if(IS_SET(target_pattern,i)){
                    applyGateX(target_register[i]);
                }
            }
        }

        void encodeBinToSuperpos_unique(const std::vector<std::size_t>& reg_memory,
                const std::vector<std::size_t>& reg_auxiliary,
                const std::vector<std::size_t>& bin_patterns,
                const std::size_t len_bin_pattern){

            EncodeBinIntoSuperpos<DerivedType> encoder(bin_patterns.size(), len_bin_pattern);
            encoder.encodeBinInToSuperpos_unique(static_cast<DerivedType&>(*this), reg_memory, reg_auxiliary, bin_patterns);
        }

        void applyHammingDistanceRotY(std::size_t test_pattern, 
                const std::vector<std::size_t> reg_mem, 
                const std::vector<std::size_t> reg_auxiliary,  
                std::size_t len_bin_pattern){

            assert(len_bin_pattern < reg_auxiliary.size()-1);

            // Encode test pattern to auxiliary register
            encodeToRegister(test_pattern, reg_auxiliary, len_bin_pattern);

            HammingDistance<DerivedType> hamming_operator(len_bin_pattern);
            hamming_operator.computeHammingDistanceRotY(static_cast<DerivedType&>(*this), reg_mem, reg_auxiliary, len_bin_pattern);

            // Un-encode test pattern from auxiliary register
            encodeToRegister(test_pattern, reg_auxiliary, len_bin_pattern);
        }


        void applyHammingDistanceOverwrite(std::size_t test_pattern, 
                const std::vector<std::size_t> reg_mem, 
                const std::vector<std::size_t> reg_auxiliary,  
                std::size_t len_bin_pattern){

            assert(reg_mem.size() < reg_auxiliary.size()-1);

            // Encode test pattern to auxiliary register
            encodeToRegister(test_pattern, reg_auxiliary, len_bin_pattern);

            HammingDistance<DerivedType>::computeHammingDistanceOverwriteAux(static_cast<DerivedType&>(*this), reg_mem, reg_auxiliary);
        }

        bool applyMeasurement(std::size_t target, bool normalize=true){
            return static_cast<DerivedType*>(this)->applyMeasurement(target, normalize);
        }

        std::size_t applyMeasurementToRegister(std::vector<std::size_t> target_qubits, bool normalize=true){
            // Store current state of training register in it's integer format
            std::size_t val = 0; 
            for(int j = target_qubits.size() - 1; j > -1; j--){
                val |= (static_cast<DerivedType*>(this)->applyMeasurement(target_qubits[j], normalize) << j);
            } 
            return val;
        }

        void groupQubits(const std::vector<std::size_t> reg_auxiliary, bool lsb=true){
            assert( reg_auxiliary.size() >= 2);

            const std::vector<std::size_t> reg_ctrl ( reg_auxiliary.end()-2, reg_auxiliary.end() );
            const std::vector<std::size_t> sub_reg (reg_auxiliary.begin(), reg_auxiliary.end()-2 ) ;

            BitGroup<DerivedType>::bit_group(static_cast<DerivedType&>(*this), sub_reg, reg_ctrl, lsb);
        }

        void collapseToBasisZ(std::size_t target, bool collapseValue){
            static_cast<DerivedType*>(this)->collapseToBasisZ(target, collapseValue);
        }
                
        void initRegister(){
            static_cast<DerivedType&>(*this).initRegister();
        }

        void initCaches(){
            #if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
            std::experimental::any_cast<NCU<DerivedType>&>(sim_ncu).initialiseMaps(static_cast<DerivedType&>(*this), 16);
            #else
            std::any_cast<NCU<DerivedType>&>(sim_ncu).initialiseMaps(static_cast<DerivedType&>(*this), 16);
            #endif
        }

        template<class Mat2x2Type>
        void addUToCache(std::string gateLabel, const Mat2x2Type& U){
            #if defined(__INTEL_COMPILER) || defined(__INTEL_LLVM_COMPILER)
            std::experimental::any_cast<NCU<DerivedType>&>(sim_ncu).getGateCache().addToCache(static_cast<DerivedType&>(*this), gateLabel, U, 16);
            #else
            std::any_cast<NCU<DerivedType>&>(sim_ncu).getGateCache().addToCache(static_cast<DerivedType&>(*this), gateLabel, U, 16);
            #endif
        }

        void PrintStates(std::string x, std::vector<std::size_t> qubits = {}){
            static_cast<DerivedType*>(this)->PrintStates(x, qubits);
        }

        template<class Mat2x2Type>
        Mat2x2Type matrixSqrt(const Mat2x2Type& U){
            Mat2x2Type V(U);
            std::complex<double> delta = U(0,0)*U(1,1) - U(0,1)*U(1,0);
            std::complex<double> tau = U(0,0) + U(1,1);
            std::complex<double> s = sqrt(delta);
            std::complex<double> t = sqrt(tau + 2.0*s);

            //must be a way to vectorise these; TinyMatrix have a scale/shift option?
            V(0,0) += s;
            V(1,1) += s;
            std::complex<double> scale_factor(1.,0.);
            scale_factor/=t;
            V(0,0) *= scale_factor; //(std::complex<double>(1.,0.)/t);
            V(0,1) *= scale_factor; //(1/t);
            V(1,0) *= scale_factor; //(1/t);
            V(1,1) *= scale_factor; //(1/t);

            return V;
        }

        template<class Mat2x2Type>
        static Mat2x2Type adjointMatrix(const Mat2x2Type& U){
            Mat2x2Type Uadjoint(U);
            std::complex<double> tmp;
            tmp = Uadjoint(0,1);
            Uadjoint(0,1) = Uadjoint(1,0);
            Uadjoint(1,0) = tmp;
            Uadjoint(0,0) = std::conj(Uadjoint(0,0));
            Uadjoint(0,1) = std::conj(Uadjoint(0,1));
            Uadjoint(1,0) = std::conj(Uadjoint(1,0));
            Uadjoint(1,1) = std::conj(Uadjoint(1,1));
            return Uadjoint;
        }

        void InvertRegister(const unsigned int minIdx, const unsigned int maxIdx){
            unsigned int range2 = ((maxIdx - minIdx)%2 == 1) ? (maxIdx - minIdx)/2 +1 : (maxIdx - minIdx)/2;
            for(unsigned int idx = 0; idx < range2; idx++){
                applyGateSwap(minIdx+idx, maxIdx-idx);
            }
        }

        #ifdef GATE_LOGGING

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