Add support for n-controlled gates

.github/workflows/cpp_build_with_cmake.yml

@@ -26,4 +26,4 @@ jobs:
       run: cmake --build build
 
     - name: unit test
-      run: ./build/bin/utest
+      run: ./build/bin/utest --gtest_death_test_style=threadsafe

examples/CMakeLists.txt

@@ -18,6 +18,10 @@ target_link_libraries(noisy_circuit_test.exe iqs)
 add_executable(test_of_custom_gates.exe test_of_custom_gates.cpp)
 target_link_libraries(test_of_custom_gates.exe iqs)
 
+add_executable(test_of_ncu_gates.exe test_of_ncu_gates.cpp)
+target_link_libraries(test_of_ncu_gates.exe iqs)
+
+
 ################################################################################
 
 set_target_properties( benchgates.exe
@@ -26,6 +30,7 @@ set_target_properties( benchgates.exe
                        heisenberg_dynamics.exe
                        noisy_circuit_test.exe
                        test_of_custom_gates.exe
+                       test_of_ncu_gates.exe
     PROPERTIES
     RUNTIME_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/bin"
 )

examples/test_of_ncu_gates.cpp

@@ -0,0 +1,105 @@
+/**
+ * @file test_of_ncu_gates.cpp
+ * @author Lee J. O'Riordan ([email protected])
+ * @author Myles Doyle ([email protected])
+ * @brief Application to show functionality of ApplyNCU gate call.
+ * @version 0.2
+ * @date 2020-06-12
+ */
+
+#include "../include/qureg.hpp"
+
+/////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////
+
+int main(int argc, char **argv)
+{
+    unsigned myrank=0, nprocs=1;
+    qhipster::mpi::Environment env(argc, argv);
+    myrank = env.GetStateRank();
+    nprocs = qhipster::mpi::Environment::GetStateSize();
+    if (env.IsUsefulRank() == false) return 0;
+    int num_threads = 1;
+#ifdef _OPENMP
+#pragma omp parallel
+    {
+        num_threads = omp_get_num_threads();
+    }
+#endif
+
+    // number of qubits in compute register.
+    std::size_t num_qubits_compute = 4;
+    if(argc != 2){
+        fprintf(stderr, "usage: %s <num_qubits_compute> \n", argv[0]);
+        exit(1);
+    }
+    else{
+        num_qubits_compute = atoi(argv[1]);
+    }
+
+    // pauliX gate that will be applied in NCU
+    ComplexDP zero = {0.,0.};
+    ComplexDP one = {1.,0.};
+    TM2x2<ComplexDP> pauliX{{zero,one,one,zero}};
+
+    // Setup vector to store compute and auxiliary quantum register indices.
+    // |compute reg>|auxiliary reg>
+    // Set number of auxiliary qubits to equal number of controlqubits
+    // if there are more than 4 control qubits, else auxiliary register
+    // will be empty as it will not be used in the NCU routine for
+    // optimisations.
+    std::size_t num_qubits_control = num_qubits_compute - 1;
+    std::size_t num_qubits_auxiliary = (num_qubits_compute - 2) * (num_qubits_control > 4);
+    std::vector<std::size_t> reg_compute(num_qubits_compute);
+    std::vector<std::size_t> reg_auxiliary(num_qubits_auxiliary);
+
+    // Set qubit indices of registers
+    {
+        std::size_t qubit_index = 0;
+        for(std::size_t i = 0; i < num_qubits_compute; i++){
+            reg_compute[i] = qubit_index;
+            qubit_index++;
+        }
+        for(std::size_t i = 0; i < num_qubits_auxiliary; i++){
+            reg_auxiliary[i] = qubit_index;
+            qubit_index++;
+        }
+    }
+
+    // Set qubit indices for qubits acting as control
+    std::vector<std::size_t> control_ids(num_qubits_control);
+
+    // Set vector containing indices of the qubits acting as
+    // control for the NCU gate.
+    for(std::size_t i = 0; i < num_qubits_control; i++){
+        control_ids[i] = reg_compute[i];
+    }
+
+    // Set index of target qubit
+    std::size_t target_id = num_qubits_compute - 1;
+
+    QubitRegister<ComplexDP> psi(num_qubits_compute + num_qubits_auxiliary);
+    psi.Initialize("base", 0);
+
+    {
+        psi.EnableStatistics();
+
+        // Apply a Hadamard gate to first num_qubits_compute-1
+        // qubits in the compute register.
+        for(std::size_t qubit_id = 0; qubit_id < num_qubits_compute-1; qubit_id++){
+            psi.ApplyHadamard(reg_compute[qubit_id]);
+        }
+
+        psi.Print("Before NCU");
+        psi.GetStatistics();
+
+        // Apply NCU
+        psi.ApplyNCU(pauliX, control_ids, reg_auxiliary, target_id);
+
+        // Observe only state with the first num_qubits_compute-1 
+        // qubits in the compute register set to 1 executes PauliX
+        // on the target qubit.
+        psi.Print("After NCU");
+    }
+}

include/GateCache.hpp

@@ -0,0 +1,147 @@
+/**
+ * @file Gates.hpp
+ * @author Lee J. O'Riordan ([email protected])
+ * @author Myles Doyle ([email protected])
+ * @brief Gates wrapper/storage class. Adapted from QNLP.
+ * @version 0.2
+ * @date 2020-06-01
+ */
+
+#ifndef GATECACHE
+#define GATECACHE
+
+#include <unordered_map>
+
+#include <complex>
+#include <vector>
+#include <qureg.hpp>
+
+#include <mat_ops.hpp>
+
+#include "qureg.hpp"
+
+/**
+ * @brief Class to cache intermediate matrix values used within other parts of the computation. 
+ * Heavily depended upon by NCU to store sqrt matrix values following Barenco et al. (1995) decomposition.
+ * 
+ * @tparam SimulatorType The simulator type with SimulatorGeneral as base class
+ */
+template <class Type>
+class GateCache {
+    private:
+    //using GateType = TM2x2<Type>;
+    std::size_t cache_depth;
+
+    public:
+    GateCache() : cache_depth(0) { };
+
+    GateCache(std::size_t default_depth) : cache_depth(default_depth) {
+        initCache(cache_depth);
+    }
+
+    ~GateCache(){ clearCache(); }
+
+    //Take the 2x2 matrix type from the template SimulatorType
+    using GateType = TM2x2<Type>;
+
+    //Maintain a map for each gate label (X,Y,Z, etc.), and use vectors to store sqrt (indexed by 1/2^(i), and pairing matrix and adjoint)
+    std::unordered_map<std::string, std::vector< std::pair<GateType, GateType> > > gateCacheMap;
+
+    void clearCache(){
+        gateCacheMap.clear();
+        cache_depth = 0;
+    }
+
+    /**
+     * @brief Initialise the gate cache with PauliX,Y,Z and H up to a given sqrt depth
+     * 
+     * @param sqrt_depth The depth to which calculate sqrt matrices and their respective adjoints
+     */
+    void initCache(const std::size_t sqrt_depth){
+        // If we do not have a sufficient circuit depth, clear and rebuild up to given depth.
+
+        if(cache_depth < sqrt_depth ){
+            gateCacheMap.clear();
+            cache_depth = 0;
+        }
+
+        if (gateCacheMap.empty()){
+            gateCacheMap["X"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_x(), adjointMatrix( construct_pauli_x() ) ) };
+            gateCacheMap["Y"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_y(), adjointMatrix( construct_pauli_y() ) ) };
+            gateCacheMap["Z"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_z(), adjointMatrix( construct_pauli_z() ) ) };
+            gateCacheMap["H"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_hadamard(), adjointMatrix( construct_hadamard() ) ) };
+
+            for( std::size_t depth = 1; depth <= sqrt_depth; depth++ ){
+                for( auto& kv : gateCacheMap ){
+                    kv.second.reserve(sqrt_depth + 1);
+                    auto m = matrixSqrt<GateType>(kv.second[depth-1].first);
+                    kv.second.emplace(kv.second.begin() + depth, std::make_pair( m, adjointMatrix( m ) ) );
+                }
+            }
+            cache_depth = sqrt_depth;
+        }
+    }
+
+    /**
+     * @brief Adds new gate to the cache up to a given sqrt depth
+     * 
+     * @param gateLabel Label of gate to index into map
+     * @param gate Gate matrix
+     * @param max_depth Depth of calculations for sqrt and associate adjoints
+     */
+    void addToCache(const std::string gateLabel, const GateType& gate, std::size_t max_depth){
+        if(max_depth <= cache_depth && gateCacheMap.find(gateLabel) != gateCacheMap.end() ){
+            return;
+        }
+        else if(max_depth > cache_depth){
+            initCache(max_depth);
+        }
+
+        std::vector< std::pair<GateType, GateType> > v;
+
+        v.reserve(max_depth + 1);
+        v.push_back(std::make_pair( gate, adjointMatrix( gate ) ) );
+        for( std::size_t depth = 1; depth <= max_depth; depth++ ){
+            auto m = matrixSqrt<GateType>( v[depth-1].first );
+            v.emplace(v.begin() + depth, std::make_pair( m, adjointMatrix( m ) ) );
+        }
+        gateCacheMap.emplace(std::make_pair(gateLabel, v) );
+    }
+
+    constexpr GateType construct_pauli_x(){
+        GateType px;
+        px(0, 0) = Type(0., 0.);
+        px(0, 1) = Type(1., 0.);
+        px(1, 0) = Type(1., 0.);
+        px(1, 1) = Type(0., 0.);
+        return px;
+    }
+    constexpr GateType construct_pauli_y(){
+        GateType py;
+        py(0, 0) = Type(0., 0.);
+        py(0, 1) = Type(0., -1.);
+        py(1, 0) = Type(0., 1.);
+        py(1, 1) = Type(0., 0.);
+        return py;
+    }
+    constexpr GateType construct_pauli_z(){
+        GateType pz;
+        pz(0, 0) = Type(1., 0.);
+        pz(0, 1) = Type(0., 0.);
+        pz(1, 0) = Type(0., 0.);
+        pz(1, 1) = Type(-1., 0.);
+        return pz;
+    }
+    constexpr GateType construct_hadamard(){
+        GateType h;
+        auto f = Type(1. / std::sqrt(2.), 0.0);
+        h(0, 0) = h(0, 1) = h(1, 0) = f;
+        h(1, 1) = -f;
+        return h;
+    }
+};
+
+template class GateCache<ComplexDP>;
+template class GateCache<ComplexSP>;
+
+#endif

include/mat_ops.hpp

@@ -0,0 +1,65 @@
+
+/**
+ * @file mat_ops.hpp
+ * @author Lee J. O'Riordan ([email protected])
+ * @brief Templated methods to manipulate small matrices. Adapted from QNLP.
+ * @version 0.2
+ * @date 2020-06-01
+ */
+
+#ifndef MAT_OPS
+#define MAT_OPS
+
+#include <complex>
+#include <vector>
+#include <iostream>
+
+/**
+ * @brief Calculates the unitary matrix square root (U == VV, where V is returned)
+ * 
+ * @tparam Type ComplexDP or ComplexSP
+ * @param U Unitary matrix to be rooted
+ * @return Matrix V such that VV == U
+ */
+template <class Mat2x2Type>            
+const 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;
+}
+
+/**
+ * @brief Function to calculate the adjoint of an input matrix
+ * 
+ * @tparam Type ComplexDP or ComplexSP
+ * @param U Unitary matrix to be adjointed
+ * @return qhipster::TinyMatrix<Type, 2, 2, 32> U^{\dagger}
+ */
+template <class Mat2x2Type>            
+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;
+}
+#endif