ncu.hpp

Go to the documentation of this file.

#ifndef QNLP_NCU
#define QNLP_NCU

#include <unordered_map>
#include <complex>
#include <cassert>
#include <utility>
#include <vector>
#include <iostream>

#include "GateCache.hpp"
#include "mat_ops.hpp"

namespace QNLP{
    template <class SimulatorType>
    class NCU{
        private:
        std::unordered_set<std::string> default_gates {"X", "Y", "Z", "I", "H"};
        struct OptParamsCX {
            std::size_t n;
            std::size_t m;
            std::size_t l;
            std::size_t num_ops;
        };

        std::unordered_map<std::size_t, std::size_t> op_call_counts_CX ;
        std::unordered_map<std::size_t, OptParamsCX> opt_op_call_params_CX ;

        //Goal: replace above hashmaps, as well as sqrtMAtrices etc with the GateCache object, which holds them all.
        GateCache<SimulatorType> gate_cache;

        protected:
            //Take the 2x2 matrix type from the template SimulatorType
            using Mat2x2Type = decltype(std::declval<SimulatorType>().getGateX());

            static std::size_t num_gate_ops;

        public:
            NCU() {
                //Optimised building block routines for nCX
                if(op_call_counts_CX.size() == 0){
                    op_call_counts_CX[3] = 13;
                    op_call_counts_CX[5] = 60;
                }
            };

            NCU(SimulatorType& qSim) : NCU(){
                gate_cache = GateCache<SimulatorType>(qSim);
            }

            ~NCU(){
                clearMaps();
                gate_cache.clearCache();
            };

            void initialiseMaps( SimulatorType& qSim,  std::size_t num_ctrl_lines){
                gate_cache.initCache(qSim, num_ctrl_lines);
            }

            void addToMaps( SimulatorType& qSim, std::string U_label, const Mat2x2Type& U, std::size_t num_ctrl_lines){
                gate_cache.addToCache(qSim, U_label, U, num_ctrl_lines);
            }

            GateCache<SimulatorType>& getGateCache(){
                return gate_cache;
            }

            void clearMaps(){
                gate_cache.clearCache();
            }


            void applyNQubitControl(SimulatorType& qSim, 
                    const std::vector<std::size_t> ctrlIndices,
                    const std::vector<std::size_t> auxIndices,
                    const unsigned int qTarget,
                    const std::string gateLabel,
                    const Mat2x2Type& U,
                    const std::size_t depth
            ){
                //No safety checks; be aware of what is physically possible (qTarget not in control_indices)
                int local_depth = depth + 1;

                //Determine the range over which the qubits exist; consider as a count of the control ops, hence +1 since extremeties included
                std::size_t cOps = ctrlIndices.size();

                //This can be generalised to find a more optimal set of splitting parameters. Using default example breakdown given in Barenco et al ('95) for now.
/*                if( (cOps >= 6) && (auxIndices.size() >= 1) && (gateLabel == "X") && (depth == 0) ){
                    auto params = find_optimal_params( ctrlIndices.size() );

                    //Need at least 1 aux qubit to perform decomposition
                    assert( auxIndices.size() > 0);
                    
                    //Gate step 1 setup
                    assert( ctrlIndices.size() >= params.m );
                    std::vector<std::size_t> subMCtrlIndicesNCX(ctrlIndices.begin(), ctrlIndices.begin() + params.m);
                    std::vector<std::size_t> subMAuxIndicesNCX(ctrlIndices.begin() + params.m,ctrlIndices.end());
                    subMAuxIndicesNCX.push_back(qTarget);

                    //Gate step 2 setup
                    std::vector<std::size_t> subLCtrlIndicesNCX(ctrlIndices.begin() + params.m, ctrlIndices.end());
                    subLCtrlIndicesNCX.push_back(auxIndices[0]);
                    std::vector<std::size_t> subLAuxIndicesNCX(ctrlIndices.begin(), ctrlIndices.begin() + params.m);

                    applyNQubitControl(qSim, subMCtrlIndicesNCX, subMAuxIndicesNCX, auxIndices[0], "X", qSim.getGateX(), 0 );
                    applyNQubitControl(qSim, subLCtrlIndicesNCX, subLAuxIndicesNCX, qTarget, "X", qSim.getGateX(), 0 );
                    applyNQubitControl(qSim, subMCtrlIndicesNCX, subMAuxIndicesNCX, auxIndices[0], "X", qSim.getGateX(), 0 );
                    applyNQubitControl(qSim, subLCtrlIndicesNCX, subLAuxIndicesNCX, qTarget, "X", qSim.getGateX(), 0 );

                }
*/
                if( (cOps >= 5) && ( auxIndices.size() >= cOps-2 ) && (gateLabel == "X") && (depth == 0) ){ //161 -> 60 2-qubit gate calls
                    qSim.applyGateCCX( ctrlIndices.back(), *(auxIndices.begin() + ctrlIndices.size() - 3), qTarget);

                    for (std::size_t i = ctrlIndices.size()-2; i >= 2; i--){
                        qSim.applyGateCCX( *(ctrlIndices.begin()+i), *(auxIndices.begin() + (i-2)), *(auxIndices.begin() + (i-1)));
                    }
                    qSim.applyGateCCX( *(ctrlIndices.begin()), *(ctrlIndices.begin()+1), *(auxIndices.begin()) );
                    
                    for (std::size_t i = 2; i <= ctrlIndices.size()-2; i++){
                        qSim.applyGateCCX( *(ctrlIndices.begin()+i), *(auxIndices.begin()+(i-2)), *(auxIndices.begin()+(i-1)));
                    }
                    qSim.applyGateCCX( ctrlIndices.back(), *(auxIndices.begin() + ctrlIndices.size() - 3), qTarget);

                    for (std::size_t i = ctrlIndices.size()-2; i >= 2; i--){
                        qSim.applyGateCCX( *(ctrlIndices.begin()+i), *(auxIndices.begin() + (i-2)), *(auxIndices.begin() + (i-1)));
                    }
                    qSim.applyGateCCX( *(ctrlIndices.begin()), *(ctrlIndices.begin()+1), *(auxIndices.begin()) );
                    for (std::size_t i = 2; i <= ctrlIndices.size()-2; i++){
                        qSim.applyGateCCX( *(ctrlIndices.begin()+i), *(auxIndices.begin()+(i-2)), *(auxIndices.begin()+(i-1)));
                    }
                }

                else if(cOps == 3){ //Optimisation for replacing 17 with 13 2-qubit gate calls
                    //Apply the 13 2-qubit gate calls
                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].first, ctrlIndices[0], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[0], ctrlIndices[1]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].second, ctrlIndices[1], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[0], ctrlIndices[1]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].first, ctrlIndices[1], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[1], ctrlIndices[2]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].second, ctrlIndices[2], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[0], ctrlIndices[2]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].first, ctrlIndices[2], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[1], ctrlIndices[2]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].second, ctrlIndices[2], qTarget, gateLabel );
                    qSim.applyGateCX( ctrlIndices[0], ctrlIndices[2]);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth+1].first, ctrlIndices[2], qTarget, gateLabel );
                }

                else if (cOps >= 2 && cOps !=3){
                    std::vector<std::size_t> subCtrlIndices(ctrlIndices.begin(), ctrlIndices.end()-1);

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth].first, ctrlIndices.back(), qTarget, gateLabel );

                    applyNQubitControl(qSim, subCtrlIndices, auxIndices, ctrlIndices.back(), "X", qSim.getGateX(), 0 );

                    qSim.applyGateCU( gate_cache.gateCacheMap[gateLabel][local_depth].second, ctrlIndices.back(), qTarget, gateLabel );

                    applyNQubitControl(qSim, subCtrlIndices, auxIndices, ctrlIndices.back(), "X", qSim.getGateX(), 0 );

                    applyNQubitControl(qSim, subCtrlIndices, auxIndices, qTarget, gateLabel, gate_cache.gateCacheMap[gateLabel][local_depth+1].first, local_depth );
                }

                //If the number of control qubits is less than 2, assume we have decomposed sufficiently
                else{
                    qSim.applyGateCU(gate_cache.gateCacheMap[gateLabel][depth].first, ctrlIndices[0], qTarget, gateLabel); //The first decomposed matrix value is used here
                }
            }

        std::size_t get_op_calls(std::size_t num_ctrl_lines){
            std::size_t num_ops = 0;
            auto it = op_call_counts_CX.find(num_ctrl_lines);
            if ( it != op_call_counts_CX.end() ){
                std::cout << "get_op_calls CTRL=" << num_ctrl_lines << "  OPS=" << it->second << std::endl;

                return it->second;
            }
            else if (num_ctrl_lines >= 2){
                num_ops += (2 + 3*get_op_calls(num_ctrl_lines-1));
            }
            else{
                num_ops += 1;
            }
            op_call_counts_CX[num_ctrl_lines] = num_ops;
            std::cout << "get_op_calls CTRL=" << num_ctrl_lines << "  OPS=" << num_ops << std::endl;

            return num_ops;
        }

        inline std::size_t get_ops_for_params(std::size_t l, std::size_t m){
            return 2*(get_op_calls(m) + get_op_calls(l));
        }

        OptParamsCX find_optimal_params(std::size_t num_ctrl_lines){

            //If entry exists, return it
            auto it = opt_op_call_params_CX.find(num_ctrl_lines);
            if ( it != opt_op_call_params_CX.end() ){
                std::cout << "FOUND OPTIMAL PARAMS SEARCH FOR CTRL=" << num_ctrl_lines << "  M=" << it->second.m << " L=" << it->second.l << " OPS=" << it->second.num_ops << std::endl;

                return it->second;
            }
            
            //Determine number of comparisons to make.
            std::size_t num_comp = static_cast<std::size_t>(std::floor((float)(num_ctrl_lines)/2.0)) +1;

            OptParamsCX optimal_params;
            std::size_t tmp_num_ops;
            
            for(std::size_t _m = 2; _m <= num_comp; ++_m){
                std::size_t _l = (num_ctrl_lines - _m + 1);
                tmp_num_ops = get_ops_for_params(_l, _m);
                if (_m == 2){
                    optimal_params = OptParamsCX{num_ctrl_lines, _m, _l, tmp_num_ops};
                }
                else if ( tmp_num_ops < optimal_params.num_ops){
                    optimal_params = OptParamsCX{num_ctrl_lines, _m, _l, tmp_num_ops};
                }
                std::cout << "OPTIMAL PARAMS SEARCH FOR CTRL=" << num_ctrl_lines << "  M=" << _m << " L=" << _l << " OPS=" << tmp_num_ops << std::endl;
            }
            opt_op_call_params_CX[num_ctrl_lines] = optimal_params;
            op_call_counts_CX[num_ctrl_lines] = optimal_params.num_ops;
            std::cout << "OPTIMAL PARAMS FOR CTRL=" << num_ctrl_lines << "  M=" << optimal_params.m << " L=" << optimal_params.l << " OPS=" << optimal_params.num_ops << std::endl;
            return optimal_params;
        }
    };

};
#endif