diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index c647ca1..0956ad1 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -10,7 +10,7 @@ on: pull_request_review: types: [submitted, edited] workflow_dispatch: - release: + jobs: diff --git a/.gitignore b/.gitignore index 252f33e..ab99f8f 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,6 @@ cmake *.egg-info test .vscode -thirdparty \ No newline at end of file +thirdparty +.pyd +MANIFEST.in diff --git a/CMakeLists.txt b/CMakeLists.txt index cf6b53c..73da6f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,10 +2,11 @@ cmake_minimum_required(VERSION 3.14...3.22) -project(qfvm LANGUAGES CXX C) +project(qfvm LANGUAGES CXX C) set (CMAKE_BUILD_TYPE Release) - +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CUDA_ARCHITECTURES 70;75;80;90) if(SKBUILD) execute_process( @@ -42,9 +43,6 @@ ExternalProject_Add(Eigen3 PREFIX ${EIGEN3_ROOT} GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git GIT_TAG 3.3.9 - # CONFIGURE_COMMAND cd ${EIGEN3_ROOT}/src/Eigen3 && cmake -B build -DCMAKE_INSTALL_PREFIX=${EIGEN3_ROOT} - # BUILD_COMMAND "" - # INSTALL_COMMAND cd ${EIGEN3_ROOT}/src/Eigen3 && cmake --build build --target install CONFIGURE_COMMAND "" BUILD_COMMAND "" @@ -72,6 +70,7 @@ if(CMAKE_HOST_SYSTEM_PROCESSOR STREQUAL "x86_64" OR CMAKE_HOST_SYSTEM_PROCESSOR endif() endif() +list (APPEND PRJ_INCLUDE_DIRS src/qfvm) pybind11_add_module(${PROJECT_NAME} MODULE src/${PROJECT_NAME}/${PROJECT_NAME}.cpp) add_dependencies(${PROJECT_NAME} Eigen3) #must add dependence for ninja target_compile_options(${PROJECT_NAME} PUBLIC ${PRJ_COMPILE_OPTIONS}) @@ -80,4 +79,49 @@ target_link_libraries(${PROJECT_NAME} PUBLIC ${PRJ_LIBRARIES}) set_target_properties(${PROJECT_NAME} PROPERTIES SUFFIX ${PYTHON_MODULE_EXTENSION}) target_compile_definitions(${PROJECT_NAME} PRIVATE VERSION_INFO=${PROJECT_VERSION}) -install(TARGETS ${PROJECT_NAME} DESTINATION .) +#GPU version +if (USE_GPU) + add_compile_definitions(_USE_GPU) + enable_language(CUDA) + set_source_files_properties(src/${PROJECT_NAME}/${PROJECT_NAME}.cpp PROPERTIES LANGUAGE CUDA) + target_link_libraries(${PROJECT_NAME} PUBLIC cudart) + target_compile_options(${PROJECT_NAME} PUBLIC $<$:--extended-lambda> ) + target_include_directories(${PROJECT_NAME} PUBLIC src/qfvm_gpu) + target_include_directories(${PROJECT_NAME} PUBLIC src/qfvm_gpu/cuda_utils) + target_include_directories(${PROJECT_NAME} PUBLIC ${CUDA_INCLUDE_DIRS}) + message("cuda_include" ${CUDA_INCLUDE_DIRS}) + if (USE_CUQUANTUM) + add_compile_definitions(_USE_CUQUANTUM) + function(set_with_fallback VARIABLE FALLBACK) + if (NOT DEFINED ${VARIABLE} OR ${VARIABLE} STREQUAL "") + set(${VARIABLE} $ENV{${VARIABLE}} CACHE INTERNAL ${VARIABLE}) + if (${VARIABLE} STREQUAL "") + if (NOT ${FALLBACK} STREQUAL "") + set(${VARIABLE} $ENV{${FALLBACK}} CACHE INTERNAL ${VARIABLE}) + endif () + endif () + endif () + endfunction() + + set_with_fallback(CUSTATEVEC_ROOT CUQUANTUM_ROOT) + + if (CUSTATEVEC_ROOT STREQUAL "") + message(FATAL_ERROR "Please set the environment variables CUSTATEVEC_ROOT or CUQUANTUM_ROOT to the path of the cuQuantum installation.") + endif () + + message(STATUS "Using CUSTATEVEC_ROOT = ${CUSTATEVEC_ROOT}") + + set(CMAKE_CUDA_FLAGS_ARCH_SM70 "-gencode arch=compute_70,code=sm_70") + set(CMAKE_CUDA_FLAGS_ARCH_SM75 "-gencode arch=compute_75,code=sm_75") + set(CMAKE_CUDA_FLAGS_ARCH_SM80 "-gencode arch=compute_80,code=sm_80 -gencode arch=compute_80,code=compute_80") + set(CMAKE_CUDA_FLAGS_ARCH_SM90 "-gencode arch=compute_90,code=sm_90 -gencode arch=compute_90,code=compute_90") + set(CMAKE_CUDA_FLAGS_ARCH "${CMAKE_CUDA_FLAGS_ARCH_SM70} ${CMAKE_CUDA_FLAGS_ARCH_SM75} ${CMAKE_CUDA_FLAGS_ARCH_SM80} ${CMAKE_CUDA_FLAGS_ARCH_SM90}") + set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} ${CMAKE_CUDA_FLAGS_ARCH}") + + target_include_directories(${PROJECT_NAME} PUBLIC ${CUDA_INCLUDE_DIRS} ${CUSTATEVEC_ROOT}/include) + target_link_directories(${PROJECT_NAME} PUBLIC ${CUSTATEVEC_ROOT}/lib ${CUSTATEVEC_ROOT}/lib64) + target_link_libraries(${PROJECT_NAME} PUBLIC -lcustatevec_static -lcublas ) + endif() +endif() + +install(TARGETS ${PROJECT_NAME} DESTINATION .) \ No newline at end of file diff --git a/README.md b/README.md index 30d92e6..caecdc0 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,18 @@ pip install -r requirements.txt python setup.py install ``` +## GPU support +To install PyQuafu with GPU-based circuit simulator, you need build from the source and make sure that [CUDA Toolkit](https://developer.nvidia.com/cuda-downloads) is installed. You can run + +``` +python setup.py install -DUSE_GPU=ON +``` +to install the GPU version. If you further have [cuQuantum](https://developer.nvidia.com/cuquantum-sdk) installed, you can install PyQuafu with cuQuantum support. +``` +python setup.py install -DUSE_GPU=ON -DUSE_CUQUANTUM=ON +``` + + ## Document Please see the website [docs](https://scq-cloud.github.io/). diff --git a/requirements.txt b/requirements.txt index 1b12b08..1419459 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,4 +6,6 @@ scipy>=1.8.1 setuptools>=58.0.4 sparse>=0.13.0 scikit-build>=0.16.1 -pybind11>=2.10.3 \ No newline at end of file +pybind11>=2.10.3 +ply~=3.11 +Pillow~=10.0.0 \ No newline at end of file diff --git a/setup.py b/setup.py index 9dc1c98..deecdf1 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ setup( name="pyquafu", - version="0.2.11", + version="0.3.0", author="ssli", author_email="ssli@iphy.ac.cn", url="https://github.com/ScQ-Cloud/pyquafu", @@ -45,4 +45,4 @@ zip_safe=False, setup_cfg=True, license="Apache-2.0 License" -) \ No newline at end of file +) diff --git a/src/qfvm/circuit.hpp b/src/qfvm/circuit.hpp index 214da35..4d576ee 100644 --- a/src/qfvm/circuit.hpp +++ b/src/qfvm/circuit.hpp @@ -89,7 +89,7 @@ class Circuit{ private: uint qubit_num_; vector gates_{0}; - + uint max_targe_num_; public: Circuit(); explicit Circuit(uint qubit_num); @@ -99,6 +99,7 @@ class Circuit{ void add_gate(QuantumOperator &gate); void compress_gates(); uint qubit_num() const { return qubit_num_; } + uint max_targe_num() const {return max_targe_num_;} vectorgates() const { return gates_; } }; @@ -121,16 +122,21 @@ void Circuit::add_gate(QuantumOperator &gate){ Circuit::Circuit(vector &gates) : - gates_(gates){ + gates_(gates), + max_targe_num_(0){ qubit_num_ = 0; for (auto gate : gates){ for (pos_t pos : gate.positions()){ + if (gate.targe_num() > max_targe_num_) + max_targe_num_ = gate.targe_num(); if (pos+1 > qubit_num_){ qubit_num_ = pos+1; } } } } Circuit::Circuit(py::object const&pycircuit) +: +max_targe_num_(0) { auto pygates = pycircuit.attr("gates"); auto used_qubits = pycircuit.attr("used_qubits").cast>(); @@ -139,6 +145,8 @@ Circuit::Circuit(py::object const&pycircuit) py::object pygate = py::reinterpret_borrow(pygate_h); QuantumOperator gate = from_pyops(pygate); if (gate){ + if (gate.targe_num() > max_targe_num_) + max_targe_num_ = gate.targe_num(); gates_.push_back(std::move(gate)); } } diff --git a/src/qfvm/qfvm.cpp b/src/qfvm/qfvm.cpp index 741cfea..e33245a 100644 --- a/src/qfvm/qfvm.cpp +++ b/src/qfvm/qfvm.cpp @@ -2,42 +2,111 @@ #include #include "simulator.hpp" +#ifdef _USE_GPU +#include +#endif + +#ifdef _USE_CUQUANTUM +#include +#endif + namespace py = pybind11; template -py::array_t to_numpy(std::vector &&src) { - vector* src_ptr = new std::vector(std::move(src)); - auto capsule = py::capsule(src_ptr, [](void* p) { delete reinterpret_cast*>(p); }); - return py::array_t( - src_ptr->size(), // shape of array - src_ptr->data(), // c-style contiguous strides for vector - capsule // numpy array references this parent - ); -} +py::array_t to_numpy(const std::tuple &src) { + auto src_ptr = std::get<0>(src); + auto src_size = std::get<1>(src); + auto capsule = py::capsule(src_ptr, [](void* p) { + delete [] reinterpret_cast(p); + }); + return py::array_t( + src_size, + src_ptr, + capsule + ); +} py::object execute(string qasm){ - return to_numpy(simulate(qasm).move_data()); + return to_numpy(simulate(qasm).move_data_to_python()); } -py::object simulate_circuit(py::object const&pycircuit, vector> const&inputstate){ - auto circuit = Circuit(pycircuit); - if (inputstate.size() == 0){ +py::object simulate_circuit(py::object const&pycircuit, py::array_t> &np_inputstate){ + auto circuit = Circuit(pycircuit); + py::buffer_info buf = np_inputstate.request(); + auto* data_ptr = reinterpret_cast*>(buf.ptr); + size_t data_size = buf.size; + + if (data_size == 0){ StateVector state; simulate(circuit, state); - return to_numpy(state.move_data()); + return to_numpy(state.move_data_to_python()); } else{ - StateVector state{inputstate}; + StateVector state(data_ptr, buf.size); simulate(circuit, state); - return to_numpy(state.move_data()); + state.move_data_to_python(); + return np_inputstate; + } +} + +#ifdef _USE_GPU +py::object simulate_circuit_gpu(py::object const&pycircuit, py::array_t> &np_inputstate){ + auto circuit = Circuit(pycircuit); + py::buffer_info buf = np_inputstate.request(); + auto* data_ptr = reinterpret_cast*>(buf.ptr); + size_t data_size = buf.size; + + + if (data_size == 0){ + StateVector state; + simulate_gpu(circuit, state); + return to_numpy(state.move_data_to_python()); + } + else{ + StateVector state(data_ptr, buf.size); + simulate_gpu(circuit, state); + state.move_data_to_python(); + return np_inputstate; + } +} +#endif + +#ifdef _USE_CUQUANTUM +py::object simulate_circuit_custate(py::object const&pycircuit, py::array_t> &np_inputstate){ + auto circuit = Circuit(pycircuit); + py::buffer_info buf = np_inputstate.request(); + auto* data_ptr = reinterpret_cast*>(buf.ptr); + size_t data_size = buf.size; + + + if (data_size == 0){ + StateVector state; + simulate_custate(circuit, state); + return to_numpy(state.move_data_to_python()); + } + else{ + StateVector state(data_ptr, buf.size); + simulate_custate(circuit, state); + state.move_data_to_python(); + return np_inputstate; } } +#endif + PYBIND11_MODULE(qfvm, m) { m.doc() = "Qfvm simulator"; m.def("execute", &execute, "Simulate with qasm"); m.def("simulate_circuit", &simulate_circuit, "Simulate with circuit", py::arg("circuit"), py::arg("inputstate")= py::array_t>(0)); + + #ifdef _USE_GPU + m.def("simulate_circuit_gpu", &simulate_circuit_gpu, "Simulate with circuit", py::arg("circuit"), py::arg("inputstate")= py::array_t>(0)); + #endif + + #ifdef _USE_CUQUANTUM + m.def("simulate_circuit_custate", &simulate_circuit_custate, "Simulate with circuit", py::arg("circuit"), py::arg("inputstate")= py::array_t>(0)); + #endif } diff --git a/src/qfvm/simulator.hpp b/src/qfvm/simulator.hpp index 49e19d7..568e895 100644 --- a/src/qfvm/simulator.hpp +++ b/src/qfvm/simulator.hpp @@ -66,11 +66,11 @@ void simulate(Circuit const& circuit, StateVector & state){ state.apply_rz(op.positions()[1], op.paras()[0]); state.apply_cnot(op.positions()[0], op.positions()[1]); break; - + //Other general gate - default: - { + default: + { if (op.targe_num() == 1){ auto mat_temp = op.mat(); complex *mat = mat_temp.data(); @@ -176,7 +176,7 @@ void simulate(string qasm, StateVector & state){ StateVector simulate(string qasm){ StateVectorstate; simulate(qasm, state); - return std::move(state); + return std::move(state); } diff --git a/src/qfvm/statevector.hpp b/src/qfvm/statevector.hpp index 977659c..a8a7f83 100644 --- a/src/qfvm/statevector.hpp +++ b/src/qfvm/statevector.hpp @@ -20,13 +20,14 @@ class StateVector{ private: uint num_; size_t size_; - vector> data_; + std::unique_ptr[]> data_; public: //construct function StateVector(); explicit StateVector(uint num); - explicit StateVector(vector> const&data); + explicit StateVector(complex *data, size_t data_size); + //Named gate function void apply_x(pos_t pos); @@ -48,7 +49,7 @@ class StateVector{ void apply_cry(pos_t control, pos_t targe, real_t theta); void apply_ccx(pos_t control1, pos_t control2, pos_t targe); void apply_swap(pos_t q1, pos_t q2); - + //General implementation //One-target gate, ctrl_num equal 2 represent multi-controlled gate template @@ -65,19 +66,25 @@ class StateVector{ complex operator[] (size_t j) const ; void set_num(uint num); - vector> move_data(){ return std::move(data_); } - void print_state(); + void print_state(); + std::tuple*, size_t> move_data_to_python() { + auto data_ptr = data_.release(); + return std::make_tuple(std::move(data_ptr), size_); + } + + complex* data(){ return data_.get(); } + size_t size(){ return size_; } + uint num(){ return num_; } }; //////// constructors /////// template -StateVector::StateVector(uint num) -: num_(num), -size_(std::pow(2, num)), -data_(size_) -{ +StateVector::StateVector(uint num) +: num_(num), +size_(1ULL<[]>(size_); data_[0] = complex(1., 0); }; @@ -85,14 +92,16 @@ template StateVector::StateVector() : StateVector(0){ } template -StateVector::StateVector(vector> const&data) +StateVector::StateVector(complex *data, size_t data_size) : -size_(data.size()), -data_(data) -{ +data_(data), +size_(data_size) +{ num_ = static_cast(std::log2(size_)); } + + //// useful functions ///// template std::complex StateVector::operator[] (size_t j) const{ @@ -101,15 +110,25 @@ std::complex StateVector::operator[] (size_t j) const{ template void StateVector::set_num(uint num){ - num_ = num; - size_ = std::pow(2, num_); - data_.resize(size_); - } + if (num_ > 0) { + // Initialized from statevector, + // should not resize + return; + } + num_ = num; + + if (size_ != 1ULL << num) { + data_.reset(); + size_ = 1ULL << num; + data_ = std::make_unique[]>(size_); + data_[0] = complex(1, 0); + } +} template void StateVector::print_state(){ - for (auto i : data_){ - std::cout << i << std::endl; + for (auto i=0;i::apply_x(pos_t pos){ const size_t offset = 1<>1; if (pos == 0){ //single step -#ifdef USE_SIMD +#ifdef USE_SIMD #pragma omp parallel for for(omp_i j = 0;j < size_;j+=2){ - double* ptr = (double*)(data_.data() + j); + double* ptr = (double*)(data_.get() + j); __m256d data = _mm256_loadu_pd(ptr); - data = _mm256_permute4x64_pd(data, 78); + data = _mm256_permute4x64_pd(data, 78); _mm256_storeu_pd(ptr, data); } #else @@ -141,8 +160,8 @@ void StateVector::apply_x(pos_t pos){ #pragma omp parallel for for(omp_i j = 0;j < rsize;j += 2){ size_t i = (j&(offset-1)) | (j>>pos<::apply_y(pos_t pos){ __m256d minus_half = _mm256_set_pd(1, -1, -1, 1); #pragma omp parallel for for(omp_i j = 0;j < size_;j+=2){ - double* ptr = (double*)(data_.data() + j); + double* ptr = (double*)(data_.get() + j); __m256d data = _mm256_loadu_pd(ptr); data = _mm256_permute4x64_pd(data, 27); data = _mm256_mul_pd(data, minus_half); @@ -179,7 +198,7 @@ void StateVector::apply_y(pos_t pos){ #else #pragma omp parallel for for(omp_i j = 0;j < size_;j+=2){ - complex temp = data_[j]; + complex temp = data_[j]; data_[j] = -im*data_[j+1]; data_[j+1] = im*temp; } @@ -189,16 +208,16 @@ void StateVector::apply_y(pos_t pos){ #ifdef USE_SIMD __m256d minus_even = _mm256_set_pd(1, -1, 1, -1); __m256d minus_odd = _mm256_set_pd(-1, 1, -1, 1); - + #pragma omp parallel for for(omp_i j = 0;j < rsize;j += 2){ size_t i = (j&(offset-1)) | (j>>pos<::apply_y(pos_t pos){ for(omp_i j = 0;j < rsize;j += 2){ size_t i = (j&(offset-1)) | (j>>pos< temp = data_[i]; + complex temp = data_[i]; data_[i] = -im*data_[i+offset]; data_[i+offset] = im*temp; - complex temp1 = data_[i1]; + complex temp1 = data_[i1]; data_[i1] = -im*data_[i1+offset]; data_[i1+offset] = im*temp1; } @@ -237,7 +256,7 @@ void StateVector::apply_z(pos_t pos){ #pragma omp parallel for for(omp_i j = 0;j < rsize;j += 2){ size_t i = (j&(offset-1)) | (j>>pos<::apply_cp(pos_t control, pos_t targe, real_t phase){ template void StateVector::apply_crx(pos_t control, pos_t targe, real_t theta){ complex mat[4] = {std::cos(theta/2), -imag_I*std::sin(theta/2), -imag_I*std::sin(theta/2), std::cos(theta/2)}; - + apply_one_targe_gate_general<1>(vector{control, targe}, mat); } template void StateVector::apply_cry(pos_t control, pos_t targe, real_t theta){ complex mat[4] = {std::cos(theta/2), -std::sin(theta/2),std::sin(theta/2), std::cos(theta/2)}; - + apply_one_targe_gate_real<1>(vector{control, targe}, mat); } @@ -361,7 +380,7 @@ void StateVector::apply_ccx(pos_t control1, pos_t control2, pos_t targe) template template void StateVector::apply_one_targe_gate_general(vector const& posv, complex *mat) -{ +{ std::function getind_func_near; std::function getind_func; size_t rsize; @@ -378,14 +397,14 @@ void StateVector::apply_one_targe_gate_general(vector const& posv getind_func_near = [&](size_t j)-> size_t { return 2*j; }; - + getind_func = [&](size_t j)-> size_t { return (j&(offset-1)) | (j>>targe<::apply_one_targe_gate_general(vector const& posv getind_func_near = getind_func; - + } else if(ctrl_num == 2){ has_control = true; @@ -426,14 +445,14 @@ void StateVector::apply_one_targe_gate_general(vector const& posv } return i; }; - getind_func_near = getind_func; + getind_func_near = getind_func; } - + const complex mat00 = mat[0]; const complex mat01 = mat[1]; const complex mat10 = mat[2]; const complex mat11 = mat[3]; - if (targe == 0){ + if (targe == 0){ #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func_near(j); @@ -462,11 +481,11 @@ void StateVector::apply_one_targe_gate_general(vector const& posv __m256d m_11re = _mm256_set_pd(mat[3].real(), mat[3].real(), mat[3].real(), mat[3].real()); __m256d m_11im = _mm256_set_pd(mat[3].imag(), -mat[3].imag(), mat[3].imag(), -mat[3].imag()); #pragma omp parallel for - for(omp_i j = 0;j < rsize; j+= 2){ + for(omp_i j = 0;j < rsize; j+= 2){ size_t i = getind_func(j); - - double* p0 = (double*)(data_.data()+i); - double* p1 = (double*)(data_.data()+i+offset); + + double* p0 = (double*)(data_.get()+i); + double* p1 = (double*)(data_.get()+i+offset); //load data __m256d data0 = _mm256_loadu_pd(p0); //lre_0, lim_0, rre_0, rim_0 __m256d data1 = _mm256_loadu_pd(p1); //lre_1, lim_1, rre_1, rim_1 @@ -508,7 +527,7 @@ void StateVector::apply_one_targe_gate_general(vector const& posv } #endif } -} +} template @@ -532,7 +551,7 @@ void StateVector::apply_one_targe_gate_x(vector const& posv) getind_func_near = [&](size_t j)-> size_t { return 2*j; }; - + getind_func = [&](size_t j)-> size_t { return (j&(offset-1)) | (j>>targe<::apply_one_targe_gate_x(vector const& posv) } return i; }; - getind_func_near = getind_func; + getind_func_near = getind_func; } - + if (targe == 0){ -#ifdef USE_SIMD +#ifdef USE_SIMD #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func_near(j); - double* ptr = (double*)(data_.data() + i); + double* ptr = (double*)(data_.get() + i); __m256d data = _mm256_loadu_pd(ptr); - data = _mm256_permute4x64_pd(data, 78); + data = _mm256_permute4x64_pd(data, 78); _mm256_storeu_pd(ptr, data); } #else @@ -594,7 +613,7 @@ void StateVector::apply_one_targe_gate_x(vector const& posv) size_t i = getind_func(j); std::swap(data_[i], data_[i+1]); } -#endif +#endif }else if (has_control && control == 0){ //single step #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ @@ -605,10 +624,10 @@ void StateVector::apply_one_targe_gate_x(vector const& posv) }else{//unroll to 2 #ifdef USE_SIMD #pragma omp parallel for - for(omp_i j = 0;j < rsize; j+= 2){ + for(omp_i j = 0;j < rsize; j+= 2){ size_t i = getind_func(j); - double* ptr0 = (double*)(data_.data() + i); - double* ptr1 = (double*)(data_.data() + i + offset); + double* ptr0 = (double*)(data_.get() + i); + double* ptr1 = (double*)(data_.get() + i + offset); __m256d data0 = _mm256_loadu_pd(ptr0); __m256d data1 = _mm256_loadu_pd(ptr1); _mm256_storeu_pd(ptr1, data0); @@ -624,7 +643,7 @@ void StateVector::apply_one_targe_gate_x(vector const& posv) } #endif } -} +} template template @@ -646,14 +665,14 @@ void StateVector::apply_one_targe_gate_real(vector const& posv, c getind_func_near = [&](size_t j)-> size_t { return 2*j; }; - + getind_func = [&](size_t j)-> size_t { return (j&(offset-1)) | (j>>targe<::apply_one_targe_gate_real(vector const& posv, c } return i; }; - getind_func_near = getind_func; + getind_func_near = getind_func; } - + const double mat00 = mat[0].real(); const double mat01 = mat[1].real(); const double mat10 = mat[2].real(); const double mat11 = mat[3].real(); - if (targe == 0){ + if (targe == 0){ #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func_near(j); @@ -706,9 +725,9 @@ void StateVector::apply_one_targe_gate_real(vector const& posv, c data_[i] = mat00*data_[i] + mat01*data_[i+1]; data_[i+1] = mat10*temp + mat11*data_[i+1]; } - + }else if (has_control && control == 0){ //single step - + #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func(j); @@ -723,11 +742,11 @@ void StateVector::apply_one_targe_gate_real(vector const& posv, c __m256d m_10re = _mm256_set_pd(mat[2].real(), mat[2].real(), mat[2].real(), mat[2].real()); __m256d m_11re = _mm256_set_pd(mat[3].real(), mat[3].real(), mat[3].real(), mat[3].real()); #pragma omp parallel for - for(omp_i j = 0;j < rsize; j+= 2){ + for(omp_i j = 0;j < rsize; j+= 2){ size_t i = getind_func(j); - - double* p0 = (double*)(data_.data()+i); - double* p1 = (double*)(data_.data()+i+offset); + + double* p0 = (double*)(data_.get()+i); + double* p1 = (double*)(data_.get()+i+offset); //load data __m256d data0 = _mm256_loadu_pd(p0); //lre_0, lim_0, rre_0, rim_0 __m256d data1 = _mm256_loadu_pd(p1); //lre_1, lim_1, rre_1, rim_1 @@ -761,7 +780,7 @@ void StateVector::apply_one_targe_gate_real(vector const& posv, c } #endif } -} +} template @@ -784,14 +803,14 @@ void StateVector::apply_one_targe_gate_diag(vector const& posv, c getind_func_near = [&](size_t j)-> size_t { return 2*j; }; - + getind_func = [&](size_t j)-> size_t { return (j&(offset-1)) | (j>>targe<::apply_one_targe_gate_diag(vector const& posv, c }; getind_func_near = getind_func; - + } else if(ctrl_num == 2){ has_control = true; @@ -831,19 +850,19 @@ void StateVector::apply_one_targe_gate_diag(vector const& posv, c } return i; }; - getind_func_near = getind_func; + getind_func_near = getind_func; } - - if (targe == 0){ + + if (targe == 0){ #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func_near(j); data_[i] *= mat[0]; data_[i+1] *= mat[1]; } - + }else if (has_control && control == 0){ //single step - + #pragma omp parallel for for(omp_i j = 0;j < rsize;j++){ size_t i = getind_func(j); @@ -859,11 +878,11 @@ void StateVector::apply_one_targe_gate_diag(vector const& posv, c __m256d m_11re = _mm256_set_pd(mat[1].real(), mat[1].real(), mat[1].real(), mat[1].real()); __m256d m_11im = _mm256_set_pd(mat[1].imag(), -mat[1].imag(), mat[1].imag(), -mat[1].imag()); #pragma omp parallel for - for(omp_i j = 0;j < rsize; j+= 2){ + for(omp_i j = 0;j < rsize; j+= 2){ size_t i = getind_func(j); - - double* p0 = (double*)(data_.data()+i); - double* p1 = (double*)(data_.data()+i+offset); + + double* p0 = (double*)(data_.get()+i); + double* p1 = (double*)(data_.get()+i+offset); //load data __m256d data0 = _mm256_loadu_pd(p0); //lre_0, lim_0, rre_0, rim_0 @@ -896,31 +915,30 @@ void StateVector::apply_one_targe_gate_diag(vector const& posv, c } #endif } -} +} template void StateVector::apply_multi_targe_gate_general(vector const& posv, uint control_num, RowMatrixXcd const& mat) { auto posv_sorted = posv; - auto targ_sorted = vector(posv.begin()+control_num, posv.end()); + auto targs = vector(posv.begin()+control_num, posv.end()); sort(posv_sorted.begin(), posv_sorted.end()); - sort(targ_sorted.begin(), targ_sorted.end()); size_t rsize = size_ >> posv.size(); - uint targe_num = targ_sorted.size(); + uint targe_num = targs.size(); size_t matsize= 1<< targe_num; std::vector targ_mask(matsize); //create target mask for (size_t m = 0; m < matsize;m++){ for (size_t j = 0; j < targe_num; j++){ if ((m>>j)&1){ - auto mask_pos = targ_sorted[j]; + auto mask_pos = targs[j]; targ_mask[m] |= 1ll< #include #include +#include #include #ifdef _MSC_VER @@ -29,4 +30,4 @@ using std::string; using RowMatrixXcd = Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; const complex imag_I = complex(0, 1.); -const double PI = 3.14159265358979323846; \ No newline at end of file +const double PI = 3.14159265358979323846; diff --git a/src/qfvm_gpu/apply_gate_custate.cuh b/src/qfvm_gpu/apply_gate_custate.cuh new file mode 100644 index 0000000..2d9c2e7 --- /dev/null +++ b/src/qfvm_gpu/apply_gate_custate.cuh @@ -0,0 +1,46 @@ + +#pragma once +#include +#include + + +void apply_gate_custate(cuDoubleComplex *psi_d, QuantumOperator &op, int n) +{ + + //get information form op + auto pos = op.positions(); + const int nTargets = op.targe_num(); + const int nControls = op.control_num(); + const int adjoint = 0; + + vector targets{pos.begin()+nControls, pos.end()}; + vector controls{pos.begin(), pos.begin()+nControls}; + + auto mat_temp = op.mat(); + cuDoubleComplex *mat = + reinterpret_cast(mat_temp.data()); + + // custatevec handle initialization + custatevecHandle_t handle; + custatevecCreate(&handle) ; + void* extraWorkspace = nullptr; + size_t extraWorkspaceSizeInBytes = 0; + + // check the size of external workspace + custatevecApplyMatrixGetWorkspaceSize( + handle, CUDA_C_64F, n, mat, CUDA_C_64F, CUSTATEVEC_MATRIX_LAYOUT_ROW, + adjoint, nTargets, nControls, CUSTATEVEC_COMPUTE_64F, &extraWorkspaceSizeInBytes) ; + + // allocate external workspace if necessary + if (extraWorkspaceSizeInBytes > 0) + cudaMalloc(&extraWorkspace, extraWorkspaceSizeInBytes); + + custatevecApplyMatrix( + handle, psi_d, CUDA_C_64F, n, mat, CUDA_C_64F, + CUSTATEVEC_MATRIX_LAYOUT_ROW, adjoint, targets.data(), nTargets, controls.data(), nullptr, + nControls, CUSTATEVEC_COMPUTE_64F, extraWorkspace, extraWorkspaceSizeInBytes); + + // destroy handle + custatevecDestroy(handle); +} + \ No newline at end of file diff --git a/src/qfvm_gpu/apply_gate_gpu.cuh b/src/qfvm_gpu/apply_gate_gpu.cuh new file mode 100644 index 0000000..e2d7be0 --- /dev/null +++ b/src/qfvm_gpu/apply_gate_gpu.cuh @@ -0,0 +1,411 @@ +#pragma once +#include +#include +#include +#include +#include + +struct targeIndex +{ + size_t ind0; + size_t ind1; +}; + + +__constant__ uint posv_d[50]; +__constant__ uint posv_sorted_d[50]; +__constant__ cuDoubleComplex mat_d_const[32*32]; //If target qubit < 5, use const memory; +__constant__ uint mat_mask_d_const[32]; + + +//-------------Single-target gate------------------------------- +template +__global__ void apply_one_targe_gate_kernel(cuDoubleComplex *psi_d, Func get_index, int rsize){ + cuDoubleComplex mat00 = mat_d_const[0]; + cuDoubleComplex mat01 = mat_d_const[1]; + cuDoubleComplex mat10 = mat_d_const[2]; + cuDoubleComplex mat11 = mat_d_const[3]; + + unsigned int gridSize = blockDim.x*gridDim.x; + for (int j = blockDim.x * blockIdx.x + threadIdx.x; j < rsize;j+= gridSize){ + targeIndex ind = get_index(j); + cuDoubleComplex temp = psi_d[ind.ind0]; + psi_d[ind.ind0] = cuCadd(cuCmul(mat00, psi_d[ind.ind0]), cuCmul(mat01, psi_d[ind.ind1])); + psi_d[ind.ind1] = cuCadd(cuCmul(mat10, temp), cuCmul(mat11, psi_d[ind.ind1])); + } +} + +template +void apply_one_targe_gate_gpu(cuDoubleComplex *psi_d, QuantumOperator &op, size_t size){ + //copy mat to device + auto mat_temp = op.mat(); + cuDoubleComplex *mat = + reinterpret_cast(mat_temp.data()); + checkCudaErrors(cudaMemcpyToSymbol(mat_d_const, mat, 4*sizeof(cuDoubleComplex))); + size_t rsize; + size_t offset; + size_t targe; + size_t control; + size_t setbit; + size_t poffset; + if (ctrl_num == 0){ + targe = op.positions()[0]; + offset = 1ll<>1; + auto getind_func = [offset, targe] __device__ (size_t j)-> targeIndex { + size_t ind0 = (j&(offset-1)) | (j>>targe<>>(psi_d, getind_func, rsize); + } + else if(ctrl_num == 1){ + control = op.positions()[0]; + targe = op.positions()[1]; + offset = 1ll<targe) { + control--; + } + poffset=1ll<>2; + auto getind_func = [control, targe, poffset, offset, setbit] __device__ (size_t j) -> targeIndex { + size_t ind0 = (j>>control<<(control+1))|(j&(poffset-1)); + ind0 = (ind0 >> targe << (targe+1))|(ind0 &(offset-1))|setbit; + size_t ind1 = ind0 + offset; + return {ind0, ind1}; + }; + + + size_t blockdim = rsize <= 1024 ? rsize : 1024; + size_t griddim = rsize / blockdim; + + apply_one_targe_gate_kernel<<>>(psi_d, getind_func, rsize); + } + else if(ctrl_num == 2){ + targe = op.positions().back(); + offset = 1ll<>psize; + + vector posv_sorted = op.positions(); + std::sort(posv_sorted.begin(), posv_sorted.end()); + //Copy pos to device + checkCudaErrors(cudaMemcpyToSymbol(posv_d, op.positions().data(), psize*sizeof(uint))); + checkCudaErrors(cudaMemcpyToSymbol(posv_sorted_d, posv_sorted.data(), psize*sizeof(uint))); + + auto getind_func = [offset, psize] __device__ (size_t j)-> targeIndex{ + size_t ind0 = j; + for (pos_t k=0;k < psize;k++) + { + pos_t _pos = posv_sorted_d[k]; + ind0 = (ind0&((1ll<<_pos)-1)) | (ind0>>_pos<<_pos<<1); + } + for (pos_t k=0;k < psize-1;k++){ + ind0 |= 1ll<>>(psi_d, getind_func, rsize); + } +} + +template +__global__ void apply_2to4_targe_gate_kernel(cuDoubleComplex *psi_d, uint ctrlnum, int psize){ + constexpr uint matlen = 1<>_pos<<_pos<<1); + } + // Set control + for (size_t k=0; k < ctrlnum;k++){ + i |= 1ll< +void apply_2to4_targe_gate_gpu_const(cuDoubleComplex *psi_d, QuantumOperator &op, size_t size){ + // uint targe_num = op.targe_num(); + uint matlen = 1<(pos.begin()+op.control_num(), pos.end()); + vector targ_mask(matlen); + //create target mask + for (size_t m = 0; m < matlen;m++){ + for (size_t j = 0; j < targe_num; j++){ + if ((m>>j)&1){ + auto mask_pos = targs[j]; + targ_mask[m] |= 1ll< posv_sorted = op.positions(); + uint psize = pos.size(); + std::sort(posv_sorted.begin(),posv_sorted.end()); + //Copy pos to device + checkCudaErrors(cudaMemcpyToSymbol(posv_d, pos.data(), psize*sizeof(uint))); + checkCudaErrors(cudaMemcpyToSymbol(posv_sorted_d, posv_sorted.data(), psize*sizeof(uint))); + + //copy mat to const memory + auto mat_temp = op.mat(); + cuDoubleComplex *mat = reinterpret_cast(mat_temp.data()); + + checkCudaErrors(cudaMemcpyToSymbol(mat_d_const, mat, matlen*matlen*sizeof(cuDoubleComplex))); + checkCudaErrors(cudaMemcpyToSymbol(mat_mask_d_const, targ_mask.data(), matlen*sizeof(uint))); + size_t rsize = size>>psize; + + uint max_thread_num = targe_num < 4 ? 1024 : 512; + size_t blockdim = rsize <= max_thread_num ? rsize : max_thread_num; + size_t griddim = rsize / blockdim; + apply_2to4_targe_gate_kernel<<>>(psi_d, op.control_num(), psize); +} + + + +// ------------Large target number gate--------------- + +__global__ void apply_5_targe_gate_kernel_const(cuDoubleComplex *psi_d, uint ctrlnum, int psize, size_t size){ + uint rsize = size>>psize; + uint targnum = psize-ctrlnum; + uint matlen = (1<>_pos<<_pos<<1); + } + // Set control + for (size_t k=0; k < ctrlnum;k++){ + i |= 1ll<>1;offset > 0;offset >>=1){ + v.x += __shfl_down_sync(0xFFFFFFFF, v.x, offset); + v.y += __shfl_down_sync(0xFFFFFFFF, v.y, offset); + } + __syncthreads(); + if (!idx) psi_d[ i | mat_mask_d_const[idy]] = v; +} + + +void apply_5_targe_gate_gpu_const(cuDoubleComplex *psi_d, QuantumOperator &op, size_t size){ + uint targe_num = op.targe_num(); + uint matlen = 1<(pos.begin()+op.control_num(), pos.end()); + vector targ_mask(matlen); + //create target mask + for (size_t m = 0; m < matlen;m++){ + for (size_t j = 0; j < targe_num; j++){ + if ((m>>j)&1){ + auto mask_pos = targs[j]; + targ_mask[m] |= 1ll< posv_sorted = op.positions(); + uint psize = pos.size(); + std::sort(posv_sorted.begin(),posv_sorted.end()); + //Copy pos to device + checkCudaErrors(cudaMemcpyToSymbol(posv_d, pos.data(), psize*sizeof(uint))); + checkCudaErrors(cudaMemcpyToSymbol(posv_sorted_d, posv_sorted.data(), psize*sizeof(uint))); + + //copy mat to const memory + auto mat_temp = op.mat(); + cuDoubleComplex *mat = reinterpret_cast(mat_temp.data()); + + checkCudaErrors(cudaMemcpyToSymbol(mat_d_const, mat, matlen*matlen*sizeof(cuDoubleComplex))); + checkCudaErrors(cudaMemcpyToSymbol(mat_mask_d_const, targ_mask.data(), matlen*sizeof(uint))); + size_t rsize = size>>psize; + uint thread_num = matlen > 32 ? 32 : matlen; + dim3 blockdim = dim3(thread_num, thread_num); + apply_5_targe_gate_kernel_const<<>>(psi_d, op.control_num(), psize, size); +} + + +//For target number 6-10 +__global__ void apply_multi_targe_gate_kernel_shared(cuDoubleComplex *psi_d, uint ctrlnum, cuDoubleComplex *mat_d, uint *mat_mask_d, int psize, size_t size){ + + uint rsize = size>>psize; + uint targnum = psize-ctrlnum; + uint matlen = (1<>_pos<<_pos<<1); + } + // Set control + for (size_t k=0; k < ctrlnum;k++){ + i |= 1ll<>1;offset > 0;offset >>=1){ + v.x += __shfl_down_sync(0xFFFFFFFF, v.x, offset); + v.y += __shfl_down_sync(0xFFFFFFFF, v.y, offset); + } + __syncthreads(); + if (!idx) local_sum[y] = cuCadd(local_sum[y], v); + } + } + + for (int y = idy; y < matlen;y+=blockDim.y){ + if (!idx) psi_d[ i | mat_mask_d[y]] = local_sum[y]; + } +} + + +void apply_multi_targe_gate_gpu_shared(cuDoubleComplex *psi_d, QuantumOperator &op, cuDoubleComplex *mat_d, uint *mat_mask_d, size_t size){ + uint targe_num = op.targe_num(); + uint matlen = 1<(pos.begin()+op.control_num(), pos.end()); + vector targ_mask(matlen); + //create target mask + for (size_t m = 0; m < matlen;m++){ + for (size_t j = 0; j < targe_num; j++){ + if ((m>>j)&1){ + auto mask_pos = targs[j]; + targ_mask[m] |= 1ll< posv_sorted = pos; + std::sort(posv_sorted.begin(),posv_sorted.end()); + //Copy pos to device + checkCudaErrors(cudaMemcpyToSymbol(posv_d, pos.data(), psize*sizeof(uint))); + checkCudaErrors(cudaMemcpyToSymbol(posv_sorted_d, posv_sorted.data(), psize*sizeof(uint))); + + + //copy mat to global memory + auto mat_temp = op.mat(); + cuDoubleComplex *mat = reinterpret_cast(mat_temp.data()); + checkCudaErrors(cudaMemcpy(mat_d, mat, matlen*matlen*sizeof(cuDoubleComplex), cudaMemcpyHostToDevice)); + checkCudaErrors(cudaMemcpy(mat_mask_d, targ_mask.data(), matlen*sizeof(uint), cudaMemcpyHostToDevice)); + + size_t rsize = size>>psize; + uint thread_num = matlen > 32 ? 32 : matlen; + dim3 blockdim = dim3(thread_num, thread_num); + + apply_multi_targe_gate_kernel_shared<<>>(psi_d, op.control_num(), mat_d, mat_mask_d, psize, size); +} + +//For target number > 10 +__global__ void apply_multi_targe_gate_kernel_global(cuDoubleComplex *psi_d, cuDoubleComplex *psi_d_copy, uint ctrlnum, cuDoubleComplex *mat_d, uint *mat_mask_d, int psize, size_t size){ + uint rsize = size>>psize; + uint targnum = psize-ctrlnum; + uint matlen = (1<>_pos<<_pos<<1); + } + // Set control + for (size_t k=0; k < ctrlnum;k++){ + i |= 1ll<>1;offset > 0;offset >>=1){ + v.x += __shfl_down_sync(0xFFFFFFFF, v.x, offset); + v.y += __shfl_down_sync(0xFFFFFFFF, v.y, offset); + } + __syncthreads(); + if (!idx) v_sum = cuCadd(v_sum, v); + } + if (!idx) psi_d[ i | mat_mask_d[y]] = v_sum; + } +} + +void apply_multi_targe_gate_gpu_global(cuDoubleComplex *psi_d, cuDoubleComplex *psi_d_copy, QuantumOperator &op, cuDoubleComplex *mat_d, uint *mat_mask_d, size_t size){ + uint targe_num = op.targe_num(); + uint matlen = 1<(pos.begin()+op.control_num(), pos.end()); + vector targ_mask(matlen); + //create target mask + for (size_t m = 0; m < matlen;m++){ + for (size_t j = 0; j < targe_num; j++){ + if ((m>>j)&1){ + auto mask_pos = targs[j]; + targ_mask[m] |= 1ll< posv_sorted = pos; + std::sort(posv_sorted.begin(),posv_sorted.end()); + //Copy pos to device + checkCudaErrors(cudaMemcpyToSymbol(posv_d, pos.data(), psize*sizeof(uint))); + checkCudaErrors(cudaMemcpyToSymbol(posv_sorted_d, posv_sorted.data(), psize*sizeof(uint))); + + + //copy mat to global memory + auto mat_temp = op.mat(); + cuDoubleComplex *mat = reinterpret_cast(mat_temp.data()); + checkCudaErrors(cudaMemcpy(mat_d, mat, matlen*matlen*sizeof(cuDoubleComplex), cudaMemcpyHostToDevice)); + checkCudaErrors(cudaMemcpy(mat_mask_d, targ_mask.data(), matlen*sizeof(uint), cudaMemcpyHostToDevice)); + + size_t rsize = size>>psize; + uint thread_num = matlen > 32 ? 32 : matlen; + dim3 blockdim = dim3(thread_num, thread_num); + + apply_multi_targe_gate_kernel_global<<>>(psi_d, psi_d_copy, op.control_num(), mat_d, mat_mask_d, psize, size); + checkCudaErrors(cudaMemcpy(psi_d_copy, psi_d, size*sizeof(cuDoubleComplex), cudaMemcpyDeviceToDevice)); +} \ No newline at end of file diff --git a/src/qfvm_gpu/cuda_simulator.cuh b/src/qfvm_gpu/cuda_simulator.cuh new file mode 100644 index 0000000..99ee5e4 --- /dev/null +++ b/src/qfvm_gpu/cuda_simulator.cuh @@ -0,0 +1,91 @@ +#pragma once +#include "cuda_statevector.cuh" +#include +#include +#include +#include +#include "apply_gate_gpu.cuh" + +void simulate_gpu(Circuit & circuit, CudaStateVector & psi_d){ + size_t size = psi_d.size(); + //initialize mat + cuDoubleComplex *mat_d; + uint *mat_mask_d; + CudaStateVector psi_d_copy{}; + if (circuit.max_targe_num() > 5) + { + uint max_matlen = 1< 10){ + psi_d_copy = psi_d; + } + } + + //apply_gate + for (auto gate : circuit.gates()){ + uint targnum = gate.targe_num(); + uint ctrlnum = gate.control_num(); + + if (targnum == 1){ + if (ctrlnum == 0){ + apply_one_targe_gate_gpu<0>(psi_d.data(), gate, size); + } + else if (ctrlnum == 1){ + apply_one_targe_gate_gpu<1>(psi_d.data(), gate, size); + } + else{ + apply_one_targe_gate_gpu<2>(psi_d.data(), gate, size); + } + } + else if(targnum > 1){ + if (targnum == 2){ + apply_2to4_targe_gate_gpu_const<2>(psi_d.data(), gate, size); + } + else if(targnum == 3){ + apply_2to4_targe_gate_gpu_const<3>(psi_d.data(), gate, size); + } + else if(targnum == 4){ + apply_2to4_targe_gate_gpu_const<4>(psi_d.data(), gate, size); + } + else if (targnum == 5){ + apply_5_targe_gate_gpu_const(psi_d.data(), gate, size); + }else if (targnum > 5 && targnum <= 10){ + apply_multi_targe_gate_gpu_shared(psi_d.data(), gate, mat_d, mat_mask_d, size); + } + else{ + apply_multi_targe_gate_gpu_global(psi_d.data(), psi_d_copy.data(), gate, mat_d, mat_mask_d, size); + } + } + else{ + throw "Invalid target number"; + } + } + + //free source + if (circuit.max_targe_num() > 5){ + checkCudaErrors(cudaFree(mat_d)); + checkCudaErrors(cudaFree(mat_mask_d)); + } +} + +void simulate_gpu(Circuit & circuit, StateVector & state){ + //initialize psi + state.set_num(circuit.qubit_num()); + size_t size = state.size(); + CudaStateVector psi_d(state); + + simulate_gpu(circuit, psi_d); + cudaDeviceSynchronize(); + + //copy back + complex* psi = reinterpret_cast*>(psi_d.data()); + checkCudaErrors(cudaMemcpy(state.data(), psi, size*sizeof(complex), cudaMemcpyDeviceToHost)); + psi=nullptr; +} + +StateVector simulate_gpu(Circuit & circuit){ + StateVector state(circuit.qubit_num()); + simulate_gpu(circuit, state); + return std::move(state); +} \ No newline at end of file diff --git a/src/qfvm_gpu/cuda_statevector.cuh b/src/qfvm_gpu/cuda_statevector.cuh new file mode 100644 index 0000000..19e8cb4 --- /dev/null +++ b/src/qfvm_gpu/cuda_statevector.cuh @@ -0,0 +1,58 @@ + +#pragma once +#include +#include +#include +#include + +class CudaStateVector +{ + protected: + uint num_; + size_t size_; + cuDoubleComplex *data_; + + public: + //construct function + CudaStateVector(){checkCudaErrors(cudaMalloc(&data_, 0));} + CudaStateVector(CudaStateVector const &other); + + explicit CudaStateVector(StateVector &sv); + ~CudaStateVector() { + checkCudaErrors(cudaFree(data_)); + } + + CudaStateVector &operator=(CudaStateVector const &other) + { + num_ = other.num(); + size_ = other.size(); + checkCudaErrors(cudaMalloc(&data_, size_*sizeof(cuDoubleComplex))); + checkCudaErrors(cudaMemcpy(data_, other.data(), size_*sizeof(cuDoubleComplex), cudaMemcpyDeviceToDevice)); + return *this; + } + + cuDoubleComplex* data() const { return data_;} + size_t size() const { return size_; } + uint num() const { return num_; } +}; + + +CudaStateVector::CudaStateVector(StateVector &sv) +: +num_(sv.num()), +size_(sv.size()) +{ + cuDoubleComplex *psi_h = reinterpret_cast(sv.data()); + checkCudaErrors(cudaMalloc(&data_, size_*sizeof(cuDoubleComplex))); + checkCudaErrors(cudaMemcpy(data_, psi_h, size_*sizeof(cuDoubleComplex), cudaMemcpyHostToDevice)); +} + +CudaStateVector::CudaStateVector(CudaStateVector const &other) +: +num_(other.num()), +size_(other.size()) +{ + checkCudaErrors(cudaMalloc(&data_, size_*sizeof(cuDoubleComplex))); + checkCudaErrors(cudaMemcpy(data_, other.data(), size_*sizeof(cuDoubleComplex), cudaMemcpyDeviceToDevice)); +} + diff --git a/src/qfvm_gpu/cuda_utils/CudaTexture.h b/src/qfvm_gpu/cuda_utils/CudaTexture.h new file mode 100644 index 0000000..0661d0d --- /dev/null +++ b/src/qfvm_gpu/cuda_utils/CudaTexture.h @@ -0,0 +1,39 @@ +#pragma once + + +#include "helper_cuda.h" +#include + + +struct CudaTexture { + cudaTextureObject_t tex; + + CudaTexture(CudaTexture const &) = delete; + CudaTexture(CudaTexture &&) = default; + CudaTexture &operator=(CudaTexture const &) = delete; + CudaTexture &operator=(CudaTexture &&) = default; + + template + CudaTexture(T *dataDev, int width, int height) { + cudaTextureObject_t tex; + cudaResourceDesc resDesc; + memset(&resDesc, 0, sizeof(resDesc)); + resDesc.resType = cudaResourceTypePitch2D; + resDesc.res.pitch2D.devPtr = dataDev; + resDesc.res.pitch2D.width = width; + resDesc.res.pitch2D.height = height; + resDesc.res.pitch2D.desc = cudaCreateChannelDesc(); + resDesc.res.pitch2D.pitchInBytes = width * sizeof(T); + cudaTextureDesc texDesc; + memset(&texDesc, 0, sizeof(texDesc)); + checkCudaErrors(cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL)); + } + + ~CudaTexture() { + checkCudaErrors(cudaDestroyTextureObject(tex)); + } + + constexpr operator cudaTextureObject_t() const { + return tex; + } +}; diff --git a/src/qfvm_gpu/cuda_utils/helper_cuda.h b/src/qfvm_gpu/cuda_utils/helper_cuda.h new file mode 100644 index 0000000..516dd57 --- /dev/null +++ b/src/qfvm_gpu/cuda_utils/helper_cuda.h @@ -0,0 +1,953 @@ +/** + * Copyright 1993-2017 NVIDIA Corporation. All rights reserved. + * + * Please refer to the NVIDIA end user license agreement (EULA) associated + * with this source code for terms and conditions that govern your use of + * this software. Any use, reproduction, disclosure, or distribution of + * this software and related documentation outside the terms of the EULA + * is strictly prohibited. + * + */ + +//////////////////////////////////////////////////////////////////////////////// +// These are CUDA Helper functions for initialization and error checking + +#ifndef COMMON_HELPER_CUDA_H_ +#define COMMON_HELPER_CUDA_H_ + +#pragma once + +#include +#include +#include +#include + +#include "helper_string.h" + +#ifndef EXIT_WAIVED +#define EXIT_WAIVED 2 +#endif + +// Note, it is required that your SDK sample to include the proper header +// files, please refer the CUDA examples for examples of the needed CUDA +// headers, which may change depending on which CUDA functions are used. + +// CUDA Runtime error messages +#ifdef __DRIVER_TYPES_H__ +static const char *_cudaGetErrorEnum(cudaError_t error) { + return cudaGetErrorName(error); +} +#endif + +#ifdef CUDA_DRIVER_API +// CUDA Driver API errors +static const char *_cudaGetErrorEnum(CUresult error) { + static char unknown[] = ""; + const char *ret = NULL; + cuGetErrorName(error, &ret); + return ret ? ret : unknown; +} +#endif + +#ifdef CUBLAS_API_H_ +// cuBLAS API errors +static const char *_cudaGetErrorEnum(cublasStatus_t error) { + switch (error) { + case CUBLAS_STATUS_SUCCESS: + return "CUBLAS_STATUS_SUCCESS"; + + case CUBLAS_STATUS_NOT_INITIALIZED: + return "CUBLAS_STATUS_NOT_INITIALIZED"; + + case CUBLAS_STATUS_ALLOC_FAILED: + return "CUBLAS_STATUS_ALLOC_FAILED"; + + case CUBLAS_STATUS_INVALID_VALUE: + return "CUBLAS_STATUS_INVALID_VALUE"; + + case CUBLAS_STATUS_ARCH_MISMATCH: + return "CUBLAS_STATUS_ARCH_MISMATCH"; + + case CUBLAS_STATUS_MAPPING_ERROR: + return "CUBLAS_STATUS_MAPPING_ERROR"; + + case CUBLAS_STATUS_EXECUTION_FAILED: + return "CUBLAS_STATUS_EXECUTION_FAILED"; + + case CUBLAS_STATUS_INTERNAL_ERROR: + return "CUBLAS_STATUS_INTERNAL_ERROR"; + + case CUBLAS_STATUS_NOT_SUPPORTED: + return "CUBLAS_STATUS_NOT_SUPPORTED"; + + case CUBLAS_STATUS_LICENSE_ERROR: + return "CUBLAS_STATUS_LICENSE_ERROR"; + } + + return ""; +} +#endif + +#ifdef _CUFFT_H_ +// cuFFT API errors +static const char *_cudaGetErrorEnum(cufftResult error) { + switch (error) { + case CUFFT_SUCCESS: + return "CUFFT_SUCCESS"; + + case CUFFT_INVALID_PLAN: + return "CUFFT_INVALID_PLAN"; + + case CUFFT_ALLOC_FAILED: + return "CUFFT_ALLOC_FAILED"; + + case CUFFT_INVALID_TYPE: + return "CUFFT_INVALID_TYPE"; + + case CUFFT_INVALID_VALUE: + return "CUFFT_INVALID_VALUE"; + + case CUFFT_INTERNAL_ERROR: + return "CUFFT_INTERNAL_ERROR"; + + case CUFFT_EXEC_FAILED: + return "CUFFT_EXEC_FAILED"; + + case CUFFT_SETUP_FAILED: + return "CUFFT_SETUP_FAILED"; + + case CUFFT_INVALID_SIZE: + return "CUFFT_INVALID_SIZE"; + + case CUFFT_UNALIGNED_DATA: + return "CUFFT_UNALIGNED_DATA"; + + case CUFFT_INCOMPLETE_PARAMETER_LIST: + return "CUFFT_INCOMPLETE_PARAMETER_LIST"; + + case CUFFT_INVALID_DEVICE: + return "CUFFT_INVALID_DEVICE"; + + case CUFFT_PARSE_ERROR: + return "CUFFT_PARSE_ERROR"; + + case CUFFT_NO_WORKSPACE: + return "CUFFT_NO_WORKSPACE"; + + case CUFFT_NOT_IMPLEMENTED: + return "CUFFT_NOT_IMPLEMENTED"; + + case CUFFT_LICENSE_ERROR: + return "CUFFT_LICENSE_ERROR"; + + case CUFFT_NOT_SUPPORTED: + return "CUFFT_NOT_SUPPORTED"; + } + + return ""; +} +#endif + +#ifdef CUSPARSEAPI +// cuSPARSE API errors +static const char *_cudaGetErrorEnum(cusparseStatus_t error) { + switch (error) { + case CUSPARSE_STATUS_SUCCESS: + return "CUSPARSE_STATUS_SUCCESS"; + + case CUSPARSE_STATUS_NOT_INITIALIZED: + return "CUSPARSE_STATUS_NOT_INITIALIZED"; + + case CUSPARSE_STATUS_ALLOC_FAILED: + return "CUSPARSE_STATUS_ALLOC_FAILED"; + + case CUSPARSE_STATUS_INVALID_VALUE: + return "CUSPARSE_STATUS_INVALID_VALUE"; + + case CUSPARSE_STATUS_ARCH_MISMATCH: + return "CUSPARSE_STATUS_ARCH_MISMATCH"; + + case CUSPARSE_STATUS_MAPPING_ERROR: + return "CUSPARSE_STATUS_MAPPING_ERROR"; + + case CUSPARSE_STATUS_EXECUTION_FAILED: + return "CUSPARSE_STATUS_EXECUTION_FAILED"; + + case CUSPARSE_STATUS_INTERNAL_ERROR: + return "CUSPARSE_STATUS_INTERNAL_ERROR"; + + case CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED: + return "CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED"; + } + + return ""; +} +#endif + +#ifdef CUSOLVER_COMMON_H_ +// cuSOLVER API errors +static const char *_cudaGetErrorEnum(cusolverStatus_t error) { + switch (error) { + case CUSOLVER_STATUS_SUCCESS: + return "CUSOLVER_STATUS_SUCCESS"; + case CUSOLVER_STATUS_NOT_INITIALIZED: + return "CUSOLVER_STATUS_NOT_INITIALIZED"; + case CUSOLVER_STATUS_ALLOC_FAILED: + return "CUSOLVER_STATUS_ALLOC_FAILED"; + case CUSOLVER_STATUS_INVALID_VALUE: + return "CUSOLVER_STATUS_INVALID_VALUE"; + case CUSOLVER_STATUS_ARCH_MISMATCH: + return "CUSOLVER_STATUS_ARCH_MISMATCH"; + case CUSOLVER_STATUS_MAPPING_ERROR: + return "CUSOLVER_STATUS_MAPPING_ERROR"; + case CUSOLVER_STATUS_EXECUTION_FAILED: + return "CUSOLVER_STATUS_EXECUTION_FAILED"; + case CUSOLVER_STATUS_INTERNAL_ERROR: + return "CUSOLVER_STATUS_INTERNAL_ERROR"; + case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED: + return "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED"; + case CUSOLVER_STATUS_NOT_SUPPORTED: + return "CUSOLVER_STATUS_NOT_SUPPORTED "; + case CUSOLVER_STATUS_ZERO_PIVOT: + return "CUSOLVER_STATUS_ZERO_PIVOT"; + case CUSOLVER_STATUS_INVALID_LICENSE: + return "CUSOLVER_STATUS_INVALID_LICENSE"; + } + + return ""; +} +#endif + +#ifdef CURAND_H_ +// cuRAND API errors +static const char *_cudaGetErrorEnum(curandStatus_t error) { + switch (error) { + case CURAND_STATUS_SUCCESS: + return "CURAND_STATUS_SUCCESS"; + + case CURAND_STATUS_VERSION_MISMATCH: + return "CURAND_STATUS_VERSION_MISMATCH"; + + case CURAND_STATUS_NOT_INITIALIZED: + return "CURAND_STATUS_NOT_INITIALIZED"; + + case CURAND_STATUS_ALLOCATION_FAILED: + return "CURAND_STATUS_ALLOCATION_FAILED"; + + case CURAND_STATUS_TYPE_ERROR: + return "CURAND_STATUS_TYPE_ERROR"; + + case CURAND_STATUS_OUT_OF_RANGE: + return "CURAND_STATUS_OUT_OF_RANGE"; + + case CURAND_STATUS_LENGTH_NOT_MULTIPLE: + return "CURAND_STATUS_LENGTH_NOT_MULTIPLE"; + + case CURAND_STATUS_DOUBLE_PRECISION_REQUIRED: + return "CURAND_STATUS_DOUBLE_PRECISION_REQUIRED"; + + case CURAND_STATUS_LAUNCH_FAILURE: + return "CURAND_STATUS_LAUNCH_FAILURE"; + + case CURAND_STATUS_PREEXISTING_FAILURE: + return "CURAND_STATUS_PREEXISTING_FAILURE"; + + case CURAND_STATUS_INITIALIZATION_FAILED: + return "CURAND_STATUS_INITIALIZATION_FAILED"; + + case CURAND_STATUS_ARCH_MISMATCH: + return "CURAND_STATUS_ARCH_MISMATCH"; + + case CURAND_STATUS_INTERNAL_ERROR: + return "CURAND_STATUS_INTERNAL_ERROR"; + } + + return ""; +} +#endif + +#ifdef NVJPEGAPI +// nvJPEG API errors +static const char *_cudaGetErrorEnum(nvjpegStatus_t error) { + switch (error) { + case NVJPEG_STATUS_SUCCESS: + return "NVJPEG_STATUS_SUCCESS"; + + case NVJPEG_STATUS_NOT_INITIALIZED: + return "NVJPEG_STATUS_NOT_INITIALIZED"; + + case NVJPEG_STATUS_INVALID_PARAMETER: + return "NVJPEG_STATUS_INVALID_PARAMETER"; + + case NVJPEG_STATUS_BAD_JPEG: + return "NVJPEG_STATUS_BAD_JPEG"; + + case NVJPEG_STATUS_JPEG_NOT_SUPPORTED: + return "NVJPEG_STATUS_JPEG_NOT_SUPPORTED"; + + case NVJPEG_STATUS_ALLOCATOR_FAILURE: + return "NVJPEG_STATUS_ALLOCATOR_FAILURE"; + + case NVJPEG_STATUS_EXECUTION_FAILED: + return "NVJPEG_STATUS_EXECUTION_FAILED"; + + case NVJPEG_STATUS_ARCH_MISMATCH: + return "NVJPEG_STATUS_ARCH_MISMATCH"; + + case NVJPEG_STATUS_INTERNAL_ERROR: + return "NVJPEG_STATUS_INTERNAL_ERROR"; + } + + return ""; +} +#endif + +#ifdef NV_NPPIDEFS_H +// NPP API errors +static const char *_cudaGetErrorEnum(NppStatus error) { + switch (error) { + case NPP_NOT_SUPPORTED_MODE_ERROR: + return "NPP_NOT_SUPPORTED_MODE_ERROR"; + + case NPP_ROUND_MODE_NOT_SUPPORTED_ERROR: + return "NPP_ROUND_MODE_NOT_SUPPORTED_ERROR"; + + case NPP_RESIZE_NO_OPERATION_ERROR: + return "NPP_RESIZE_NO_OPERATION_ERROR"; + + case NPP_NOT_SUFFICIENT_COMPUTE_CAPABILITY: + return "NPP_NOT_SUFFICIENT_COMPUTE_CAPABILITY"; + +#if ((NPP_VERSION_MAJOR << 12) + (NPP_VERSION_MINOR << 4)) <= 0x5000 + + case NPP_BAD_ARG_ERROR: + return "NPP_BAD_ARGUMENT_ERROR"; + + case NPP_COEFF_ERROR: + return "NPP_COEFFICIENT_ERROR"; + + case NPP_RECT_ERROR: + return "NPP_RECTANGLE_ERROR"; + + case NPP_QUAD_ERROR: + return "NPP_QUADRANGLE_ERROR"; + + case NPP_MEM_ALLOC_ERR: + return "NPP_MEMORY_ALLOCATION_ERROR"; + + case NPP_HISTO_NUMBER_OF_LEVELS_ERROR: + return "NPP_HISTOGRAM_NUMBER_OF_LEVELS_ERROR"; + + case NPP_INVALID_INPUT: + return "NPP_INVALID_INPUT"; + + case NPP_POINTER_ERROR: + return "NPP_POINTER_ERROR"; + + case NPP_WARNING: + return "NPP_WARNING"; + + case NPP_ODD_ROI_WARNING: + return "NPP_ODD_ROI_WARNING"; +#else + + // These are for CUDA 5.5 or higher + case NPP_BAD_ARGUMENT_ERROR: + return "NPP_BAD_ARGUMENT_ERROR"; + + case NPP_COEFFICIENT_ERROR: + return "NPP_COEFFICIENT_ERROR"; + + case NPP_RECTANGLE_ERROR: + return "NPP_RECTANGLE_ERROR"; + + case NPP_QUADRANGLE_ERROR: + return "NPP_QUADRANGLE_ERROR"; + + case NPP_MEMORY_ALLOCATION_ERR: + return "NPP_MEMORY_ALLOCATION_ERROR"; + + case NPP_HISTOGRAM_NUMBER_OF_LEVELS_ERROR: + return "NPP_HISTOGRAM_NUMBER_OF_LEVELS_ERROR"; + + case NPP_INVALID_HOST_POINTER_ERROR: + return "NPP_INVALID_HOST_POINTER_ERROR"; + + case NPP_INVALID_DEVICE_POINTER_ERROR: + return "NPP_INVALID_DEVICE_POINTER_ERROR"; +#endif + + case NPP_LUT_NUMBER_OF_LEVELS_ERROR: + return "NPP_LUT_NUMBER_OF_LEVELS_ERROR"; + + case NPP_TEXTURE_BIND_ERROR: + return "NPP_TEXTURE_BIND_ERROR"; + + case NPP_WRONG_INTERSECTION_ROI_ERROR: + return "NPP_WRONG_INTERSECTION_ROI_ERROR"; + + case NPP_NOT_EVEN_STEP_ERROR: + return "NPP_NOT_EVEN_STEP_ERROR"; + + case NPP_INTERPOLATION_ERROR: + return "NPP_INTERPOLATION_ERROR"; + + case NPP_RESIZE_FACTOR_ERROR: + return "NPP_RESIZE_FACTOR_ERROR"; + + case NPP_HAAR_CLASSIFIER_PIXEL_MATCH_ERROR: + return "NPP_HAAR_CLASSIFIER_PIXEL_MATCH_ERROR"; + +#if ((NPP_VERSION_MAJOR << 12) + (NPP_VERSION_MINOR << 4)) <= 0x5000 + + case NPP_MEMFREE_ERR: + return "NPP_MEMFREE_ERR"; + + case NPP_MEMSET_ERR: + return "NPP_MEMSET_ERR"; + + case NPP_MEMCPY_ERR: + return "NPP_MEMCPY_ERROR"; + + case NPP_MIRROR_FLIP_ERR: + return "NPP_MIRROR_FLIP_ERR"; +#else + + case NPP_MEMFREE_ERROR: + return "NPP_MEMFREE_ERROR"; + + case NPP_MEMSET_ERROR: + return "NPP_MEMSET_ERROR"; + + case NPP_MEMCPY_ERROR: + return "NPP_MEMCPY_ERROR"; + + case NPP_MIRROR_FLIP_ERROR: + return "NPP_MIRROR_FLIP_ERROR"; +#endif + + case NPP_ALIGNMENT_ERROR: + return "NPP_ALIGNMENT_ERROR"; + + case NPP_STEP_ERROR: + return "NPP_STEP_ERROR"; + + case NPP_SIZE_ERROR: + return "NPP_SIZE_ERROR"; + + case NPP_NULL_POINTER_ERROR: + return "NPP_NULL_POINTER_ERROR"; + + case NPP_CUDA_KERNEL_EXECUTION_ERROR: + return "NPP_CUDA_KERNEL_EXECUTION_ERROR"; + + case NPP_NOT_IMPLEMENTED_ERROR: + return "NPP_NOT_IMPLEMENTED_ERROR"; + + case NPP_ERROR: + return "NPP_ERROR"; + + case NPP_SUCCESS: + return "NPP_SUCCESS"; + + case NPP_WRONG_INTERSECTION_QUAD_WARNING: + return "NPP_WRONG_INTERSECTION_QUAD_WARNING"; + + case NPP_MISALIGNED_DST_ROI_WARNING: + return "NPP_MISALIGNED_DST_ROI_WARNING"; + + case NPP_AFFINE_QUAD_INCORRECT_WARNING: + return "NPP_AFFINE_QUAD_INCORRECT_WARNING"; + + case NPP_DOUBLE_SIZE_WARNING: + return "NPP_DOUBLE_SIZE_WARNING"; + + case NPP_WRONG_INTERSECTION_ROI_WARNING: + return "NPP_WRONG_INTERSECTION_ROI_WARNING"; + +#if ((NPP_VERSION_MAJOR << 12) + (NPP_VERSION_MINOR << 4)) >= 0x6000 + /* These are 6.0 or higher */ + case NPP_LUT_PALETTE_BITSIZE_ERROR: + return "NPP_LUT_PALETTE_BITSIZE_ERROR"; + + case NPP_ZC_MODE_NOT_SUPPORTED_ERROR: + return "NPP_ZC_MODE_NOT_SUPPORTED_ERROR"; + + case NPP_QUALITY_INDEX_ERROR: + return "NPP_QUALITY_INDEX_ERROR"; + + case NPP_CHANNEL_ORDER_ERROR: + return "NPP_CHANNEL_ORDER_ERROR"; + + case NPP_ZERO_MASK_VALUE_ERROR: + return "NPP_ZERO_MASK_VALUE_ERROR"; + + case NPP_NUMBER_OF_CHANNELS_ERROR: + return "NPP_NUMBER_OF_CHANNELS_ERROR"; + + case NPP_COI_ERROR: + return "NPP_COI_ERROR"; + + case NPP_DIVISOR_ERROR: + return "NPP_DIVISOR_ERROR"; + + case NPP_CHANNEL_ERROR: + return "NPP_CHANNEL_ERROR"; + + case NPP_STRIDE_ERROR: + return "NPP_STRIDE_ERROR"; + + case NPP_ANCHOR_ERROR: + return "NPP_ANCHOR_ERROR"; + + case NPP_MASK_SIZE_ERROR: + return "NPP_MASK_SIZE_ERROR"; + + case NPP_MOMENT_00_ZERO_ERROR: + return "NPP_MOMENT_00_ZERO_ERROR"; + + case NPP_THRESHOLD_NEGATIVE_LEVEL_ERROR: + return "NPP_THRESHOLD_NEGATIVE_LEVEL_ERROR"; + + case NPP_THRESHOLD_ERROR: + return "NPP_THRESHOLD_ERROR"; + + case NPP_CONTEXT_MATCH_ERROR: + return "NPP_CONTEXT_MATCH_ERROR"; + + case NPP_FFT_FLAG_ERROR: + return "NPP_FFT_FLAG_ERROR"; + + case NPP_FFT_ORDER_ERROR: + return "NPP_FFT_ORDER_ERROR"; + + case NPP_SCALE_RANGE_ERROR: + return "NPP_SCALE_RANGE_ERROR"; + + case NPP_DATA_TYPE_ERROR: + return "NPP_DATA_TYPE_ERROR"; + + case NPP_OUT_OFF_RANGE_ERROR: + return "NPP_OUT_OFF_RANGE_ERROR"; + + case NPP_DIVIDE_BY_ZERO_ERROR: + return "NPP_DIVIDE_BY_ZERO_ERROR"; + + case NPP_RANGE_ERROR: + return "NPP_RANGE_ERROR"; + + case NPP_NO_MEMORY_ERROR: + return "NPP_NO_MEMORY_ERROR"; + + case NPP_ERROR_RESERVED: + return "NPP_ERROR_RESERVED"; + + case NPP_NO_OPERATION_WARNING: + return "NPP_NO_OPERATION_WARNING"; + + case NPP_DIVIDE_BY_ZERO_WARNING: + return "NPP_DIVIDE_BY_ZERO_WARNING"; +#endif + +#if ((NPP_VERSION_MAJOR << 12) + (NPP_VERSION_MINOR << 4)) >= 0x7000 + /* These are 7.0 or higher */ + case NPP_OVERFLOW_ERROR: + return "NPP_OVERFLOW_ERROR"; + + case NPP_CORRUPTED_DATA_ERROR: + return "NPP_CORRUPTED_DATA_ERROR"; +#endif + } + + return ""; +} +#endif + +template +void check(T result, char const *const func, const char *const file, + int const line) { + if (result) { + fprintf(stderr, "CUDA error at %s:%d code=%d(%s) \"%s\" \n", file, line, + static_cast(result), _cudaGetErrorEnum(result), func); + exit(EXIT_FAILURE); + } +} + +#ifdef __DRIVER_TYPES_H__ +// This will output the proper CUDA error strings in the event +// that a CUDA host call returns an error +#define checkCudaErrors(val) check((val), #val, __FILE__, __LINE__) + +// This will output the proper error string when calling cudaGetLastError +#define getLastCudaError(msg) __getLastCudaError(msg, __FILE__, __LINE__) + +inline void __getLastCudaError(const char *errorMessage, const char *file, + const int line) { + cudaError_t err = cudaGetLastError(); + + if (cudaSuccess != err) { + fprintf(stderr, + "%s(%i) : getLastCudaError() CUDA error :" + " %s : (%d) %s.\n", + file, line, errorMessage, static_cast(err), + cudaGetErrorString(err)); + exit(EXIT_FAILURE); + } +} + +// This will only print the proper error string when calling cudaGetLastError +// but not exit program incase error detected. +#define printLastCudaError(msg) __printLastCudaError(msg, __FILE__, __LINE__) + +inline void __printLastCudaError(const char *errorMessage, const char *file, + const int line) { + cudaError_t err = cudaGetLastError(); + + if (cudaSuccess != err) { + fprintf(stderr, + "%s(%i) : getLastCudaError() CUDA error :" + " %s : (%d) %s.\n", + file, line, errorMessage, static_cast(err), + cudaGetErrorString(err)); + } +} +#endif + +#ifndef MAX +#define MAX(a, b) (a > b ? a : b) +#endif + +// Float To Int conversion +inline int ftoi(float value) { + return (value >= 0 ? static_cast(value + 0.5) + : static_cast(value - 0.5)); +} + +// Beginning of GPU Architecture definitions +inline int _ConvertSMVer2Cores(int major, int minor) { + // Defines for GPU Architecture types (using the SM version to determine + // the # of cores per SM + typedef struct { + int SM; // 0xMm (hexidecimal notation), M = SM Major version, + // and m = SM minor version + int Cores; + } sSMtoCores; + + sSMtoCores nGpuArchCoresPerSM[] = { + {0x30, 192}, + {0x32, 192}, + {0x35, 192}, + {0x37, 192}, + {0x50, 128}, + {0x52, 128}, + {0x53, 128}, + {0x60, 64}, + {0x61, 128}, + {0x62, 128}, + {0x70, 64}, + {0x72, 64}, + {0x75, 64}, + {0x80, 64}, + {0x86, 128}, + {0x87, 128}, + {-1, -1}}; + + int index = 0; + + while (nGpuArchCoresPerSM[index].SM != -1) { + if (nGpuArchCoresPerSM[index].SM == ((major << 4) + minor)) { + return nGpuArchCoresPerSM[index].Cores; + } + + index++; + } + + // If we don't find the values, we default use the previous one + // to run properly + printf( + "MapSMtoCores for SM %d.%d is undefined." + " Default to use %d Cores/SM\n", + major, minor, nGpuArchCoresPerSM[index - 1].Cores); + return nGpuArchCoresPerSM[index - 1].Cores; +} + +inline const char* _ConvertSMVer2ArchName(int major, int minor) { + // Defines for GPU Architecture types (using the SM version to determine + // the GPU Arch name) + typedef struct { + int SM; // 0xMm (hexidecimal notation), M = SM Major version, + // and m = SM minor version + const char* name; + } sSMtoArchName; + + sSMtoArchName nGpuArchNameSM[] = { + {0x30, "Kepler"}, + {0x32, "Kepler"}, + {0x35, "Kepler"}, + {0x37, "Kepler"}, + {0x50, "Maxwell"}, + {0x52, "Maxwell"}, + {0x53, "Maxwell"}, + {0x60, "Pascal"}, + {0x61, "Pascal"}, + {0x62, "Pascal"}, + {0x70, "Volta"}, + {0x72, "Xavier"}, + {0x75, "Turing"}, + {0x80, "Ampere"}, + {0x86, "Ampere"}, + {0x87, "Ampere"}, + {-1, "Graphics Device"}}; + + int index = 0; + + while (nGpuArchNameSM[index].SM != -1) { + if (nGpuArchNameSM[index].SM == ((major << 4) + minor)) { + return nGpuArchNameSM[index].name; + } + + index++; + } + + // If we don't find the values, we default use the previous one + // to run properly + printf( + "MapSMtoArchName for SM %d.%d is undefined." + " Default to use %s\n", + major, minor, nGpuArchNameSM[index - 1].name); + return nGpuArchNameSM[index - 1].name; +} + // end of GPU Architecture definitions + +#ifdef __CUDA_RUNTIME_H__ +// General GPU Device CUDA Initialization +inline int gpuDeviceInit(int devID) { + int device_count; + checkCudaErrors(cudaGetDeviceCount(&device_count)); + + if (device_count == 0) { + fprintf(stderr, + "gpuDeviceInit() CUDA error: " + "no devices supporting CUDA.\n"); + exit(EXIT_FAILURE); + } + + if (devID < 0) { + devID = 0; + } + + if (devID > device_count - 1) { + fprintf(stderr, "\n"); + fprintf(stderr, ">> %d CUDA capable GPU device(s) detected. <<\n", + device_count); + fprintf(stderr, + ">> gpuDeviceInit (-device=%d) is not a valid" + " GPU device. <<\n", + devID); + fprintf(stderr, "\n"); + return -devID; + } + + int computeMode = -1, major = 0, minor = 0; + checkCudaErrors(cudaDeviceGetAttribute(&computeMode, cudaDevAttrComputeMode, devID)); + checkCudaErrors(cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, devID)); + checkCudaErrors(cudaDeviceGetAttribute(&minor, cudaDevAttrComputeCapabilityMinor, devID)); + if (computeMode == cudaComputeModeProhibited) { + fprintf(stderr, + "Error: device is running in , no threads can use cudaSetDevice().\n"); + return -1; + } + + if (major < 1) { + fprintf(stderr, "gpuDeviceInit(): GPU device does not support CUDA.\n"); + exit(EXIT_FAILURE); + } + + checkCudaErrors(cudaSetDevice(devID)); + printf("gpuDeviceInit() CUDA Device [%d]: \"%s\n", devID, _ConvertSMVer2ArchName(major, minor)); + + return devID; +} + +// This function returns the best GPU (with maximum GFLOPS) +inline int gpuGetMaxGflopsDeviceId() { + int current_device = 0, sm_per_multiproc = 0; + int max_perf_device = 0; + int device_count = 0; + int devices_prohibited = 0; + + uint64_t max_compute_perf = 0; + checkCudaErrors(cudaGetDeviceCount(&device_count)); + + if (device_count == 0) { + fprintf(stderr, + "gpuGetMaxGflopsDeviceId() CUDA error:" + " no devices supporting CUDA.\n"); + exit(EXIT_FAILURE); + } + + // Find the best CUDA capable GPU device + current_device = 0; + + while (current_device < device_count) { + int computeMode = -1, major = 0, minor = 0; + checkCudaErrors(cudaDeviceGetAttribute(&computeMode, cudaDevAttrComputeMode, current_device)); + checkCudaErrors(cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, current_device)); + checkCudaErrors(cudaDeviceGetAttribute(&minor, cudaDevAttrComputeCapabilityMinor, current_device)); + + // If this GPU is not running on Compute Mode prohibited, + // then we can add it to the list + if (computeMode != cudaComputeModeProhibited) { + if (major == 9999 && minor == 9999) { + sm_per_multiproc = 1; + } else { + sm_per_multiproc = + _ConvertSMVer2Cores(major, minor); + } + int multiProcessorCount = 0, clockRate = 0; + checkCudaErrors(cudaDeviceGetAttribute(&multiProcessorCount, cudaDevAttrMultiProcessorCount, current_device)); + cudaError_t result = cudaDeviceGetAttribute(&clockRate, cudaDevAttrClockRate, current_device); + if (result != cudaSuccess) { + // If cudaDevAttrClockRate attribute is not supported we + // set clockRate as 1, to consider GPU with most SMs and CUDA Cores. + if(result == cudaErrorInvalidValue) { + clockRate = 1; + } + else { + fprintf(stderr, "CUDA error at %s:%d code=%d(%s) \n", __FILE__, __LINE__, + static_cast(result), _cudaGetErrorEnum(result)); + exit(EXIT_FAILURE); + } + } + uint64_t compute_perf = (uint64_t)multiProcessorCount * sm_per_multiproc * clockRate; + + if (compute_perf > max_compute_perf) { + max_compute_perf = compute_perf; + max_perf_device = current_device; + } + } else { + devices_prohibited++; + } + + ++current_device; + } + + if (devices_prohibited == device_count) { + fprintf(stderr, + "gpuGetMaxGflopsDeviceId() CUDA error:" + " all devices have compute mode prohibited.\n"); + exit(EXIT_FAILURE); + } + + return max_perf_device; +} + +// Initialization code to find the best CUDA Device +inline int findCudaDevice(int argc, const char **argv) { + int devID = 0; + + // If the command-line has a device number specified, use it + if (checkCmdLineFlag(argc, argv, "device")) { + devID = getCmdLineArgumentInt(argc, argv, "device="); + + if (devID < 0) { + printf("Invalid command line parameter\n "); + exit(EXIT_FAILURE); + } else { + devID = gpuDeviceInit(devID); + + if (devID < 0) { + printf("exiting...\n"); + exit(EXIT_FAILURE); + } + } + } else { + // Otherwise pick the device with highest Gflops/s + devID = gpuGetMaxGflopsDeviceId(); + checkCudaErrors(cudaSetDevice(devID)); + int major = 0, minor = 0; + checkCudaErrors(cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, devID)); + checkCudaErrors(cudaDeviceGetAttribute(&minor, cudaDevAttrComputeCapabilityMinor, devID)); + printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", + devID, _ConvertSMVer2ArchName(major, minor), major, minor); + + } + + return devID; +} + +inline int findIntegratedGPU() { + int current_device = 0; + int device_count = 0; + int devices_prohibited = 0; + + checkCudaErrors(cudaGetDeviceCount(&device_count)); + + if (device_count == 0) { + fprintf(stderr, "CUDA error: no devices supporting CUDA.\n"); + exit(EXIT_FAILURE); + } + + // Find the integrated GPU which is compute capable + while (current_device < device_count) { + int computeMode = -1, integrated = -1; + checkCudaErrors(cudaDeviceGetAttribute(&computeMode, cudaDevAttrComputeMode, current_device)); + checkCudaErrors(cudaDeviceGetAttribute(&integrated, cudaDevAttrIntegrated, current_device)); + // If GPU is integrated and is not running on Compute Mode prohibited, + // then cuda can map to GLES resource + if (integrated && (computeMode != cudaComputeModeProhibited)) { + checkCudaErrors(cudaSetDevice(current_device)); + + int major = 0, minor = 0; + checkCudaErrors(cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, current_device)); + checkCudaErrors(cudaDeviceGetAttribute(&minor, cudaDevAttrComputeCapabilityMinor, current_device)); + printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", + current_device, _ConvertSMVer2ArchName(major, minor), major, minor); + + return current_device; + } else { + devices_prohibited++; + } + + current_device++; + } + + if (devices_prohibited == device_count) { + fprintf(stderr, + "CUDA error:" + " No GLES-CUDA Interop capable GPU found.\n"); + exit(EXIT_FAILURE); + } + + return -1; +} + +// General check for CUDA GPU SM Capabilities +inline bool checkCudaCapabilities(int major_version, int minor_version) { + int dev; + int major = 0, minor = 0; + + checkCudaErrors(cudaGetDevice(&dev)); + checkCudaErrors(cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, dev)); + checkCudaErrors(cudaDeviceGetAttribute(&minor, cudaDevAttrComputeCapabilityMinor, dev)); + + if ((major > major_version) || + (major == major_version && + minor >= minor_version)) { + printf(" Device %d: <%16s >, Compute SM %d.%d detected\n", dev, + _ConvertSMVer2ArchName(major, minor), major, minor); + return true; + } else { + printf( + " No GPU device was found that can support " + "CUDA compute capability %d.%d.\n", + major_version, minor_version); + return false; + } +} +#endif + + // end of CUDA Helper Functions + +#endif // COMMON_HELPER_CUDA_H_ diff --git a/src/qfvm_gpu/cuda_utils/helper_custatevec.hpp b/src/qfvm_gpu/cuda_utils/helper_custatevec.hpp new file mode 100644 index 0000000..59fdbc3 --- /dev/null +++ b/src/qfvm_gpu/cuda_utils/helper_custatevec.hpp @@ -0,0 +1,31 @@ +/* + * Copyright (c) 2021-2022, NVIDIA CORPORATION & AFFILIATES. + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#define HANDLE_ERROR(x) \ +{ const auto err = x; \ + if (err != CUSTATEVEC_STATUS_SUCCESS ) { \ + printf("Error: %s in line %d\n", \ + custatevecGetErrorString(err), __LINE__); return err; } \ +}; + +#define HANDLE_CUDA_ERROR(x) \ +{ const auto err = x; \ + if (err != cudaSuccess ) { \ + printf("Error: %s in line %d\n", \ + cudaGetErrorString(err), __LINE__); return err; } \ +}; + +bool almost_equal(cuDoubleComplex x, cuDoubleComplex y) { + const double eps = 1.0e-5; + const cuDoubleComplex diff = cuCsub(x, y); + return (cuCabs(diff) < eps); +} + +bool almost_equal(double x, double y) { + const double eps = 1.0e-5; + const double diff = x - y; + return (abs(diff) < eps); +} diff --git a/src/qfvm_gpu/cuda_utils/helper_string.h b/src/qfvm_gpu/cuda_utils/helper_string.h new file mode 100644 index 0000000..77864b8 --- /dev/null +++ b/src/qfvm_gpu/cuda_utils/helper_string.h @@ -0,0 +1,683 @@ +/** + * Copyright 1993-2013 NVIDIA Corporation. All rights reserved. + * + * Please refer to the NVIDIA end user license agreement (EULA) associated + * with this source code for terms and conditions that govern your use of + * this software. Any use, reproduction, disclosure, or distribution of + * this software and related documentation outside the terms of the EULA + * is strictly prohibited. + * + */ + +// These are helper functions for the SDK samples (string parsing, timers, etc) +#ifndef COMMON_HELPER_STRING_H_ +#define COMMON_HELPER_STRING_H_ + +#include +#include +#include +#include + +#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) +#ifndef _CRT_SECURE_NO_DEPRECATE +#define _CRT_SECURE_NO_DEPRECATE +#endif +#ifndef STRCASECMP +#define STRCASECMP _stricmp +#endif +#ifndef STRNCASECMP +#define STRNCASECMP _strnicmp +#endif +#ifndef STRCPY +#define STRCPY(sFilePath, nLength, sPath) strcpy_s(sFilePath, nLength, sPath) +#endif + +#ifndef FOPEN +#define FOPEN(fHandle, filename, mode) fopen_s(&fHandle, filename, mode) +#endif +#ifndef FOPEN_FAIL +#define FOPEN_FAIL(result) (result != 0) +#endif +#ifndef SSCANF +#define SSCANF sscanf_s +#endif +#ifndef SPRINTF +#define SPRINTF sprintf_s +#endif +#else // Linux Includes +#include +#include + +#ifndef STRCASECMP +#define STRCASECMP strcasecmp +#endif +#ifndef STRNCASECMP +#define STRNCASECMP strncasecmp +#endif +#ifndef STRCPY +#define STRCPY(sFilePath, nLength, sPath) strcpy(sFilePath, sPath) +#endif + +#ifndef FOPEN +#define FOPEN(fHandle, filename, mode) (fHandle = fopen(filename, mode)) +#endif +#ifndef FOPEN_FAIL +#define FOPEN_FAIL(result) (result == NULL) +#endif +#ifndef SSCANF +#define SSCANF sscanf +#endif +#ifndef SPRINTF +#define SPRINTF sprintf +#endif +#endif + +#ifndef EXIT_WAIVED +#define EXIT_WAIVED 2 +#endif + +// CUDA Utility Helper Functions +inline int stringRemoveDelimiter(char delimiter, const char *string) { + int string_start = 0; + + while (string[string_start] == delimiter) { + string_start++; + } + + if (string_start >= static_cast(strlen(string) - 1)) { + return 0; + } + + return string_start; +} + +inline int getFileExtension(char *filename, char **extension) { + int string_length = static_cast(strlen(filename)); + + while (filename[string_length--] != '.') { + if (string_length == 0) break; + } + + if (string_length > 0) string_length += 2; + + if (string_length == 0) + *extension = NULL; + else + *extension = &filename[string_length]; + + return string_length; +} + +inline bool checkCmdLineFlag(const int argc, const char **argv, + const char *string_ref) { + bool bFound = false; + + if (argc >= 1) { + for (int i = 1; i < argc; i++) { + int string_start = stringRemoveDelimiter('-', argv[i]); + const char *string_argv = &argv[i][string_start]; + + const char *equal_pos = strchr(string_argv, '='); + int argv_length = static_cast( + equal_pos == 0 ? strlen(string_argv) : equal_pos - string_argv); + + int length = static_cast(strlen(string_ref)); + + if (length == argv_length && + !STRNCASECMP(string_argv, string_ref, length)) { + bFound = true; + continue; + } + } + } + + return bFound; +} + +// This function wraps the CUDA Driver API into a template function +template +inline bool getCmdLineArgumentValue(const int argc, const char **argv, + const char *string_ref, T *value) { + bool bFound = false; + + if (argc >= 1) { + for (int i = 1; i < argc; i++) { + int string_start = stringRemoveDelimiter('-', argv[i]); + const char *string_argv = &argv[i][string_start]; + int length = static_cast(strlen(string_ref)); + + if (!STRNCASECMP(string_argv, string_ref, length)) { + if (length + 1 <= static_cast(strlen(string_argv))) { + int auto_inc = (string_argv[length] == '=') ? 1 : 0; + *value = (T)atoi(&string_argv[length + auto_inc]); + } + + bFound = true; + i = argc; + } + } + } + + return bFound; +} + +inline int getCmdLineArgumentInt(const int argc, const char **argv, + const char *string_ref) { + bool bFound = false; + int value = -1; + + if (argc >= 1) { + for (int i = 1; i < argc; i++) { + int string_start = stringRemoveDelimiter('-', argv[i]); + const char *string_argv = &argv[i][string_start]; + int length = static_cast(strlen(string_ref)); + + if (!STRNCASECMP(string_argv, string_ref, length)) { + if (length + 1 <= static_cast(strlen(string_argv))) { + int auto_inc = (string_argv[length] == '=') ? 1 : 0; + value = atoi(&string_argv[length + auto_inc]); + } else { + value = 0; + } + + bFound = true; + continue; + } + } + } + + if (bFound) { + return value; + } else { + return 0; + } +} + +inline float getCmdLineArgumentFloat(const int argc, const char **argv, + const char *string_ref) { + bool bFound = false; + float value = -1; + + if (argc >= 1) { + for (int i = 1; i < argc; i++) { + int string_start = stringRemoveDelimiter('-', argv[i]); + const char *string_argv = &argv[i][string_start]; + int length = static_cast(strlen(string_ref)); + + if (!STRNCASECMP(string_argv, string_ref, length)) { + if (length + 1 <= static_cast(strlen(string_argv))) { + int auto_inc = (string_argv[length] == '=') ? 1 : 0; + value = static_cast(atof(&string_argv[length + auto_inc])); + } else { + value = 0.f; + } + + bFound = true; + continue; + } + } + } + + if (bFound) { + return value; + } else { + return 0; + } +} + +inline bool getCmdLineArgumentString(const int argc, const char **argv, + const char *string_ref, + char **string_retval) { + bool bFound = false; + + if (argc >= 1) { + for (int i = 1; i < argc; i++) { + int string_start = stringRemoveDelimiter('-', argv[i]); + char *string_argv = const_cast(&argv[i][string_start]); + int length = static_cast(strlen(string_ref)); + + if (!STRNCASECMP(string_argv, string_ref, length)) { + *string_retval = &string_argv[length + 1]; + bFound = true; + continue; + } + } + } + + if (!bFound) { + *string_retval = NULL; + } + + return bFound; +} + +////////////////////////////////////////////////////////////////////////////// +//! Find the path for a file assuming that +//! files are found in the searchPath. +//! +//! @return the path if succeeded, otherwise 0 +//! @param filename name of the file +//! @param executable_path optional absolute path of the executable +////////////////////////////////////////////////////////////////////////////// +inline char *sdkFindFilePath(const char *filename, + const char *executable_path) { + // defines a variable that is replaced with the name of the + // executable + + // Typical relative search paths to locate needed companion files (e.g. sample + // input data, or JIT source files) The origin for the relative search may be + // the .exe file, a .bat file launching an .exe, a browser .exe launching the + // .exe or .bat, etc + const char *searchPath[] = { + "./", // same dir + "./_data_files/", + "./common/", // "/common/" subdir + "./common/data/", // "/common/data/" subdir + "./data/", // "/data/" subdir + "./src/", // "/src/" subdir + "./src//data/", // "/src//data/" subdir + "./inc/", // "/inc/" subdir + "./0_Simple/", // "/0_Simple/" subdir + "./1_Utilities/", // "/1_Utilities/" subdir + "./2_Graphics/", // "/2_Graphics/" subdir + "./3_Imaging/", // "/3_Imaging/" subdir + "./4_Finance/", // "/4_Finance/" subdir + "./5_Simulations/", // "/5_Simulations/" subdir + "./6_Advanced/", // "/6_Advanced/" subdir + "./7_CUDALibraries/", // "/7_CUDALibraries/" subdir + "./8_Android/", // "/8_Android/" subdir + "./samples/", // "/samples/" subdir + + "./0_Simple//data/", // "/0_Simple//data/" + // subdir + "./1_Utilities//data/", // "/1_Utilities//data/" + // subdir + "./2_Graphics//data/", // "/2_Graphics//data/" + // subdir + "./3_Imaging//data/", // "/3_Imaging//data/" + // subdir + "./4_Finance//data/", // "/4_Finance//data/" + // subdir + "./5_Simulations//data/", // "/5_Simulations//data/" + // subdir + "./6_Advanced//data/", // "/6_Advanced//data/" + // subdir + "./7_CUDALibraries//", // "/7_CUDALibraries//" + // subdir + "./7_CUDALibraries//data/", // "/7_CUDALibraries//data/" + // subdir + + "../", // up 1 in tree + "../common/", // up 1 in tree, "/common/" subdir + "../common/data/", // up 1 in tree, "/common/data/" subdir + "../data/", // up 1 in tree, "/data/" subdir + "../src/", // up 1 in tree, "/src/" subdir + "../inc/", // up 1 in tree, "/inc/" subdir + + "../0_Simple//data/", // up 1 in tree, + // "/0_Simple//" + // subdir + "../1_Utilities//data/", // up 1 in tree, + // "/1_Utilities//" + // subdir + "../2_Graphics//data/", // up 1 in tree, + // "/2_Graphics//" + // subdir + "../3_Imaging//data/", // up 1 in tree, + // "/3_Imaging//" + // subdir + "../4_Finance//data/", // up 1 in tree, + // "/4_Finance//" + // subdir + "../5_Simulations//data/", // up 1 in tree, + // "/5_Simulations//" + // subdir + "../6_Advanced//data/", // up 1 in tree, + // "/6_Advanced//" + // subdir + "../7_CUDALibraries//data/", // up 1 in tree, + // "/7_CUDALibraries//" + // subdir + "../8_Android//data/", // up 1 in tree, + // "/8_Android//" + // subdir + "../samples//data/", // up 1 in tree, + // "/samples//" + // subdir + "../../", // up 2 in tree + "../../common/", // up 2 in tree, "/common/" subdir + "../../common/data/", // up 2 in tree, "/common/data/" subdir + "../../data/", // up 2 in tree, "/data/" subdir + "../../src/", // up 2 in tree, "/src/" subdir + "../../inc/", // up 2 in tree, "/inc/" subdir + "../../sandbox//data/", // up 2 in tree, + // "/sandbox//" + // subdir + "../../0_Simple//data/", // up 2 in tree, + // "/0_Simple//" + // subdir + "../../1_Utilities//data/", // up 2 in tree, + // "/1_Utilities//" + // subdir + "../../2_Graphics//data/", // up 2 in tree, + // "/2_Graphics//" + // subdir + "../../3_Imaging//data/", // up 2 in tree, + // "/3_Imaging//" + // subdir + "../../4_Finance//data/", // up 2 in tree, + // "/4_Finance//" + // subdir + "../../5_Simulations//data/", // up 2 in tree, + // "/5_Simulations//" + // subdir + "../../6_Advanced//data/", // up 2 in tree, + // "/6_Advanced//" + // subdir + "../../7_CUDALibraries//data/", // up 2 in tree, + // "/7_CUDALibraries//" + // subdir + "../../8_Android//data/", // up 2 in tree, + // "/8_Android//" + // subdir + "../../samples//data/", // up 2 in tree, + // "/samples//" + // subdir + "../../../", // up 3 in tree + "../../../src//", // up 3 in tree, + // "/src//" subdir + "../../../src//data/", // up 3 in tree, + // "/src//data/" + // subdir + "../../../src//src/", // up 3 in tree, + // "/src//src/" + // subdir + "../../../src//inc/", // up 3 in tree, + // "/src//inc/" + // subdir + "../../../sandbox//", // up 3 in tree, + // "/sandbox//" + // subdir + "../../../sandbox//data/", // up 3 in tree, + // "/sandbox//data/" + // subdir + "../../../sandbox//src/", // up 3 in tree, + // "/sandbox//src/" + // subdir + "../../../sandbox//inc/", // up 3 in tree, + // "/sandbox//inc/" + // subdir + "../../../0_Simple//data/", // up 3 in tree, + // "/0_Simple//" + // subdir + "../../../1_Utilities//data/", // up 3 in tree, + // "/1_Utilities//" + // subdir + "../../../2_Graphics//data/", // up 3 in tree, + // "/2_Graphics//" + // subdir + "../../../3_Imaging//data/", // up 3 in tree, + // "/3_Imaging//" + // subdir + "../../../4_Finance//data/", // up 3 in tree, + // "/4_Finance//" + // subdir + "../../../5_Simulations//data/", // up 3 in tree, + // "/5_Simulations//" + // subdir + "../../../6_Advanced//data/", // up 3 in tree, + // "/6_Advanced//" + // subdir + "../../../7_CUDALibraries//data/", // up 3 in tree, + // "/7_CUDALibraries//" + // subdir + "../../../8_Android//data/", // up 3 in tree, + // "/8_Android//" + // subdir + "../../../0_Simple//", // up 3 in tree, + // "/0_Simple//" + // subdir + "../../../1_Utilities//", // up 3 in tree, + // "/1_Utilities//" + // subdir + "../../../2_Graphics//", // up 3 in tree, + // "/2_Graphics//" + // subdir + "../../../3_Imaging//", // up 3 in tree, + // "/3_Imaging//" + // subdir + "../../../4_Finance//", // up 3 in tree, + // "/4_Finance//" + // subdir + "../../../5_Simulations//", // up 3 in tree, + // "/5_Simulations//" + // subdir + "../../../6_Advanced//", // up 3 in tree, + // "/6_Advanced//" + // subdir + "../../../7_CUDALibraries//", // up 3 in tree, + // "/7_CUDALibraries//" + // subdir + "../../../8_Android//", // up 3 in tree, + // "/8_Android//" + // subdir + "../../../samples//data/", // up 3 in tree, + // "/samples//" + // subdir + "../../../common/", // up 3 in tree, "../../../common/" subdir + "../../../common/data/", // up 3 in tree, "../../../common/data/" subdir + "../../../data/", // up 3 in tree, "../../../data/" subdir + "../../../../", // up 4 in tree + "../../../../src//", // up 4 in tree, + // "/src//" subdir + "../../../../src//data/", // up 4 in tree, + // "/src//data/" + // subdir + "../../../../src//src/", // up 4 in tree, + // "/src//src/" + // subdir + "../../../../src//inc/", // up 4 in tree, + // "/src//inc/" + // subdir + "../../../../sandbox//", // up 4 in tree, + // "/sandbox//" + // subdir + "../../../../sandbox//data/", // up 4 in tree, + // "/sandbox//data/" + // subdir + "../../../../sandbox//src/", // up 4 in tree, + // "/sandbox//src/" + // subdir + "../../../../sandbox//inc/", // up 4 in tree, + // "/sandbox//inc/" + // subdir + "../../../../0_Simple//data/", // up 4 in tree, + // "/0_Simple//" + // subdir + "../../../../1_Utilities//data/", // up 4 in tree, + // "/1_Utilities//" + // subdir + "../../../../2_Graphics//data/", // up 4 in tree, + // "/2_Graphics//" + // subdir + "../../../../3_Imaging//data/", // up 4 in tree, + // "/3_Imaging//" + // subdir + "../../../../4_Finance//data/", // up 4 in tree, + // "/4_Finance//" + // subdir + "../../../../5_Simulations//data/", // up 4 in tree, + // "/5_Simulations//" + // subdir + "../../../../6_Advanced//data/", // up 4 in tree, + // "/6_Advanced//" + // subdir + "../../../../7_CUDALibraries//data/", // up 4 in tree, + // "/7_CUDALibraries//" + // subdir + "../../../../8_Android//data/", // up 4 in tree, + // "/8_Android//" + // subdir + "../../../../0_Simple//", // up 4 in tree, + // "/0_Simple//" + // subdir + "../../../../1_Utilities//", // up 4 in tree, + // "/1_Utilities//" + // subdir + "../../../../2_Graphics//", // up 4 in tree, + // "/2_Graphics//" + // subdir + "../../../../3_Imaging//", // up 4 in tree, + // "/3_Imaging//" + // subdir + "../../../../4_Finance//", // up 4 in tree, + // "/4_Finance//" + // subdir + "../../../../5_Simulations//", // up 4 in tree, + // "/5_Simulations//" + // subdir + "../../../../6_Advanced//", // up 4 in tree, + // "/6_Advanced//" + // subdir + "../../../../7_CUDALibraries//", // up 4 in tree, + // "/7_CUDALibraries//" + // subdir + "../../../../8_Android//", // up 4 in tree, + // "/8_Android//" + // subdir + "../../../../samples//data/", // up 4 in tree, + // "/samples//" + // subdir + "../../../../common/", // up 4 in tree, "../../../common/" subdir + "../../../../common/data/", // up 4 in tree, "../../../common/data/" + // subdir + "../../../../data/", // up 4 in tree, "../../../data/" subdir + "../../../../../", // up 5 in tree + "../../../../../src//", // up 5 in tree, + // "/src//" + // subdir + "../../../../../src//data/", // up 5 in tree, + // "/src//data/" + // subdir + "../../../../../src//src/", // up 5 in tree, + // "/src//src/" + // subdir + "../../../../../src//inc/", // up 5 in tree, + // "/src//inc/" + // subdir + "../../../../../sandbox//", // up 5 in tree, + // "/sandbox//" + // subdir + "../../../../../sandbox//data/", // up 5 in tree, + // "/sandbox//data/" + // subdir + "../../../../../sandbox//src/", // up 5 in tree, + // "/sandbox//src/" + // subdir + "../../../../../sandbox//inc/", // up 5 in tree, + // "/sandbox//inc/" + // subdir + "../../../../../0_Simple//data/", // up 5 in tree, + // "/0_Simple//" + // subdir + "../../../../../1_Utilities//data/", // up 5 in tree, + // "/1_Utilities//" + // subdir + "../../../../../2_Graphics//data/", // up 5 in tree, + // "/2_Graphics//" + // subdir + "../../../../../3_Imaging//data/", // up 5 in tree, + // "/3_Imaging//" + // subdir + "../../../../../4_Finance//data/", // up 5 in tree, + // "/4_Finance//" + // subdir + "../../../../../5_Simulations//data/", // up 5 in tree, + // "/5_Simulations//" + // subdir + "../../../../../6_Advanced//data/", // up 5 in tree, + // "/6_Advanced//" + // subdir + "../../../../../7_CUDALibraries//data/", // up 5 in + // tree, + // "/7_CUDALibraries//" + // subdir + "../../../../../8_Android//data/", // up 5 in tree, + // "/8_Android//" + // subdir + "../../../../../samples//data/", // up 5 in tree, + // "/samples//" + // subdir + "../../../../../common/", // up 5 in tree, "../../../common/" subdir + "../../../../../common/data/", // up 5 in tree, "../../../common/data/" + // subdir + }; + + // Extract the executable name + std::string executable_name; + + if (executable_path != 0) { + executable_name = std::string(executable_path); + +#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) + // Windows path delimiter + size_t delimiter_pos = executable_name.find_last_of('\\'); + executable_name.erase(0, delimiter_pos + 1); + + if (executable_name.rfind(".exe") != std::string::npos) { + // we strip .exe, only if the .exe is found + executable_name.resize(executable_name.size() - 4); + } + +#else + // Linux & OSX path delimiter + size_t delimiter_pos = executable_name.find_last_of('/'); + executable_name.erase(0, delimiter_pos + 1); +#endif + } + + // Loop over all search paths and return the first hit + for (unsigned int i = 0; i < sizeof(searchPath) / sizeof(char *); ++i) { + std::string path(searchPath[i]); + size_t executable_name_pos = path.find(""); + + // If there is executable_name variable in the searchPath + // replace it with the value + if (executable_name_pos != std::string::npos) { + if (executable_path != 0) { + path.replace(executable_name_pos, strlen(""), + executable_name); + } else { + // Skip this path entry if no executable argument is given + continue; + } + } + +#ifdef _DEBUG + printf("sdkFindFilePath <%s> in %s\n", filename, path.c_str()); +#endif + + // Test if the file exists + path.append(filename); + FILE *fp; + FOPEN(fp, path.c_str(), "rb"); + + if (fp != NULL) { + fclose(fp); + // File found + // returning an allocated array here for backwards compatibility reasons + char *file_path = reinterpret_cast(malloc(path.length() + 1)); + STRCPY(file_path, path.length() + 1, path.c_str()); + return file_path; + } + + if (fp) { + fclose(fp); + } + } + + // File not found + return 0; +} + +#endif // COMMON_HELPER_STRING_H_ diff --git a/src/qfvm_gpu/cuda_utils/ticktock.h b/src/qfvm_gpu/cuda_utils/ticktock.h new file mode 100644 index 0000000..1adb8a9 --- /dev/null +++ b/src/qfvm_gpu/cuda_utils/ticktock.h @@ -0,0 +1,9 @@ +#pragma once + +//#include +//#define TICK(x) auto bench_##x = std::chrono::steady_clock::now(); +//#define TOCK(x) std::cout << #x ": " << std::chrono::duration_cast>(std::chrono::steady_clock::now() - bench_##x).count() << "s" << std::endl; + +#include +#define TICK(x) auto bench_##x = tbb::tick_count::now(); +#define TOCK(x) std::cout << #x ": " << (tbb::tick_count::now() - bench_##x).seconds() << "s" << std::endl; diff --git a/src/qfvm_gpu/custate_simu.cuh b/src/qfvm_gpu/custate_simu.cuh new file mode 100644 index 0000000..6456c7e --- /dev/null +++ b/src/qfvm_gpu/custate_simu.cuh @@ -0,0 +1,35 @@ +#pragma once +#include "cuda_statevector.cuh" +#include +#include +#include +#include "apply_gate_custate.cuh" + +void simulate_custate(Circuit & circuit, CudaStateVector & psi_d){ + size_t size = psi_d.size(); + int n = psi_d.num(); + for (auto gate : circuit.gates()){ + apply_gate_custate(psi_d.data(), gate, n); + } +} + +void simulate_custate(Circuit & circuit, StateVector & state){ + //initialize psi + state.set_num(circuit.qubit_num()); + size_t size = state.size(); + CudaStateVector psi_d(state); + + simulate_custate(circuit, psi_d); + cudaDeviceSynchronize(); + + //copy back + complex* psi = reinterpret_cast*>(psi_d.data()); + checkCudaErrors(cudaMemcpy(state.data(), psi, size*sizeof(complex), cudaMemcpyDeviceToHost)); + psi=nullptr; +} + +StateVector simulate_custate(Circuit & circuit){ + StateVector state(circuit.qubit_num()); + simulate_custate(circuit, state); + return std::move(state); +} \ No newline at end of file diff --git a/src/quafu/__init__.py b/src/quafu/__init__.py index b13f90a..7153ea2 100644 --- a/src/quafu/__init__.py +++ b/src/quafu/__init__.py @@ -1,10 +1,11 @@ from .circuits.quantum_circuit import QuantumCircuit -from .results.results import ExecResult, SimuResult +from .results.results import ExecResult, SimuResult from .tasks.tasks import Task from .users.userapi import User from .simulators.simulator import simulate __all__ = ["QuantumCircuit", "ExecResult", "Task", "User", "SimuResult", "simulate"] + def get_version(): - return "0.2.11" \ No newline at end of file + return "0.3.0" diff --git a/src/quafu/backends/backends.py b/src/quafu/backends/backends.py index cd30609..6587eaa 100644 --- a/src/quafu/backends/backends.py +++ b/src/quafu/backends/backends.py @@ -2,23 +2,28 @@ import json import re import networkx as nx -from quafu.users.userapi import load_account import numpy as np import matplotlib.pyplot as plt +from quafu.users.userapi import User + class Backend(object): - def __init__(self, backend_info : dict): + def __init__(self, backend_info: dict): self.name = backend_info['system_name'] self._valid_gates = backend_info['valid_gates'] self.qubit_num = backend_info['qubits'] self.system_id = backend_info['system_id'] + self.status = backend_info['status'] + self.qv = backend_info["QV"] + # self.task_in_queue = backend_info["task_in_queue"] - def get_chip_info(self): - api_token, url = load_account() + def get_chip_info(self, user=User()): + api_token = user.apitoken + url = user._url data = {"system_name": self.name.lower()} - headers={"api_token": api_token} - chip_info = requests.post(url = url + "qbackend/scq_get_chip_info/", data=data, - headers=headers) + headers = {"api_token": api_token} + chip_info = requests.post(url=url + user.chip_api, data=data, + headers=headers) chip_info = json.loads(chip_info.text) json_topo_struct = chip_info["topological_structure"] qubits_list = [] @@ -66,17 +71,17 @@ def get_chip_info(self): elarge = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] >= 0.9] esmall = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] < 0.9] - elarge_labels = {(u, v) : "%.3f" %d["weight"] for (u, v, d) in G.edges(data=True) if d["weight"] >= 0.9} - esmall_labels = {(u, v) : "%.3f" %d["weight"] for (u, v, d) in G.edges(data=True) if d["weight"] < 0.9} + elarge_labels = {(u, v): "%.3f" % d["weight"] for (u, v, d) in G.edges(data=True) if d["weight"] >= 0.9} + esmall_labels = {(u, v): "%.3f" % d["weight"] for (u, v, d) in G.edges(data=True) if d["weight"] < 0.9} - pos = nx.spring_layout(G, seed=1) + pos = nx.spring_layout(G, seed=1) fig, ax = plt.subplots() nx.draw_networkx_nodes(G, pos, node_size=400, ax=ax) nx.draw_networkx_edges(G, pos, edgelist=elarge, width=2, ax=ax) nx.draw_networkx_edges( G, pos, edgelist=esmall, width=2, alpha=0.5, style="dashed" - , ax=ax) + , ax=ax) nx.draw_networkx_labels(G, pos, font_size=14, font_family="sans-serif", ax=ax) edge_labels = nx.get_edge_attributes(G, "weight") @@ -85,33 +90,7 @@ def get_chip_info(self): fig.set_figwidth(14) fig.set_figheight(14) fig.tight_layout() - return {"mapping" : int_to_qubit, "topology_diagram": fig, "full_info": chip_info} - + return {"mapping": int_to_qubit, "topology_diagram": fig, "full_info": chip_info} def get_valid_gates(self): return self._valid_gates - - - - - -# class ScQ_P10(Backend): -# def __init__(self): -# super().__init__("ScQ-P10") -# self.valid_gates = ["cx", "cz", "rx", "ry", "rz", "x", "y", "z", "h", "sx", "sy", "id", "delay", "barrier", "cy", "cnot", "swap"] - -# class ScQ_P20(Backend): -# def __init__(self): -# super().__init__("ScQ-P20") -# self.valid_gates = ["cx", "cz", "rx", "ry", "rz", "x", "y", "z", "h", "sx", "sy", "id", "delay", "barrier", "cy", "cnot", "swap"] - -# class ScQ_P50(Backend): -# def __init__(self): -# super().__init__("ScQ-P50") -# self.valid_gates = ["cx", "cz", "rx", "ry", "rz", "x", "y", "z", "h"] - -# class ScQ_S41(Backend): -# def __init__(self): -# super().__init__("ScQ-S41") -# self.valid_gates = [ "rx", "ry", "rz", "x", "y", "z", "h", "sx", "sy", "id", "delay", "barrier", "xy"] - \ No newline at end of file diff --git a/src/quafu/benchmark/__init__.py b/src/quafu/benchmark/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/quafu/benchmark/adder.py b/src/quafu/benchmark/adder.py new file mode 100644 index 0000000..9e4ad37 --- /dev/null +++ b/src/quafu/benchmark/adder.py @@ -0,0 +1,126 @@ +from quafu.circuits.quantum_circuit import QuantumCircuit +import matplotlib.pyplot as plt + +""" +// quantum ripple-carry adder from Cuccaro et al, quant-ph/0410184 +OPENQASM 2.0; +include "qelib1.inc"; +""" + +n = 10 +qc = QuantumCircuit(n) + + +def majority(a, b, c): + """ + gate majority a,b,c + { + cx c,b; + cx c,a; + ccx a,b,c; + } + """ + qc.cnot(c, b) + qc.cnot(a, b) + qc.mcx([a, b], c) + + +def unmaj(a, b, c): + """ + gate unmaj a,b,c + { + ccx a,b,c; + cx c,a; + cx a,b; + } + """ + qc.mcx([a, b], c) + qc.cnot(c, a) + qc.cnot(a, b) + + +def qreg(_i, name): + """ + qreg cin[1]; + qreg a[4]; + qreg b[4]; + qreg cout[1]; + """ + if name == 'cin': + return _i + elif name == 'a': + return _i + 1 + elif name == 'b': + return _i + 5 + elif name == 'cout': + return _i + 9 + else: + raise ValueError('Unknown qreg name: {}'.format(name)) + + +def creg(_i, name): + """ + creg ans[5]; + """ + if name == 'ans': + return _i + else: + raise ValueError('Unknown creg name: {}'.format(name)) + + +""" +// set input states +x a[0]; // a = 0001 +x b; // b = 1111 +""" +qc.x(qreg(0, 'a')) +for i in range(4): + qc.x(qreg(i, 'b')) + +""" +// add a to b, storing result in b +majority cin[0],b[0],a[0]; +majority a[0],b[1],a[1]; +majority a[1],b[2],a[2]; +majority a[2],b[3],a[3]; +cx a[3],cout[0]; +unmaj a[2],b[3],a[3]; +unmaj a[1],b[2],a[2]; +unmaj a[0],b[1],a[1]; +unmaj cin[0],b[0],a[0]; +""" +majority(qreg(0, 'cin'), qreg(0, 'b'), qreg(0, 'a')) +majority(qreg(0, 'a'), qreg(1, 'b'), qreg(1, 'a')) +for i in range(1, 4): + majority(qreg(i-1, 'a'), qreg(i, 'b'), qreg(i, 'a')) +qc.cnot(qreg(3, 'a'), qreg(0, 'cout')) +unmaj(qreg(2, 'a'), qreg(3, 'b'), qreg(3, 'a')) +unmaj(qreg(1, 'a'), qreg(2, 'b'), qreg(2, 'a')) +unmaj(qreg(0, 'a'), qreg(1, 'b'), qreg(1, 'a')) + +""" +measure b[0] -> ans[0]; +measure b[1] -> ans[1]; +measure b[2] -> ans[2]; +measure b[3] -> ans[3]; +measure cout[0] -> ans[4]; +""" +measure_pos = [qreg(i, 'b') for i in range(4)] + [qreg(0, 'cout')] +measure_cbits = [creg(i, 'ans') for i in range(5)] +qc.measure(measure_pos, cbits=measure_cbits) +# qc.draw_circuit() +# print(qc.to_openqasm()) + +init_labels = dict.fromkeys(range(n)) +init_labels[0] = 'cin' +for i in range(4): + init_labels[i+1] = f'a_{i}' +for i in range(5): + init_labels[i+5] = f'b_{i}' + +end_labels = {i+5: f'ans_{i}' for i in range(5)} +qc.plot_circuit(title='Quantum ripple-carry adder', + init_labels=init_labels, + end_labels=end_labels, + ) +plt.show() diff --git a/src/quafu/benchmark/deutsch_jozsa.py b/src/quafu/benchmark/deutsch_jozsa.py new file mode 100644 index 0000000..ad327f3 --- /dev/null +++ b/src/quafu/benchmark/deutsch_jozsa.py @@ -0,0 +1,75 @@ +import random +import matplotlib.pyplot as plt +import numpy as np + +from quafu.circuits.quantum_circuit import QuantumCircuit +from quafu.visualisation.circuitPlot import CircuitPlotManager + + +def get_const_oracle(qc: QuantumCircuit): + n = qc.num - 1 + output = np.random.randint(2) + if output == 1: + qc.x(n) + qc.name = 'Constant Oracle' + return qc + + +def get_balanced_oracle(qc: QuantumCircuit): + n = qc.num - 1 + b_str = ''.join([random.choice('01') for _ in range(n)]) + + # Place X-qu_gate + for qubit in range(len(b_str)): + if b_str[qubit] == '1': + qc.x(qubit) + + # Use barrier as divider + qc.barrier() + + # Controlled-NOT qu_gate + for qubit in range(n): + qc.cnot(qubit, n) + + qc.barrier() + + # Place X-qu_gate + for qubit in range(len(b_str)): + if b_str[qubit] == '1': + qc.x(qubit) + + qc.name = 'Balanced Oracle' + return qc + + +def deutsch_jozsa(n: int, case: str): + circuit = QuantumCircuit(n + 1) # number of q-bit and c-bit + # Initialization + for qubit in range(n): + circuit.h(qubit) + circuit.x(n) + circuit.h(n) + + # Add oracle + ################################################# + if case == 'balanced': + get_balanced_oracle(circuit) + elif case == 'constant': + get_const_oracle(circuit) + else: + raise ValueError('undefined case: ' + case) + ################################################# + + # Repeat H-qu_gate + circuit.barrier() + for qubit in range(n): + circuit.h(qubit) + circuit.barrier() + + # Measure + circuit.measure(list(range(n)), list(range(n))) + return circuit + + +dj_qc = deutsch_jozsa(n=4, case='constant') +dj_qc.plot_circuit(title='Deutsch-Josza Circuit', show=True) diff --git a/src/quafu/benchmark/unitary_test.py b/src/quafu/benchmark/unitary_test.py new file mode 100644 index 0000000..f008c5d --- /dev/null +++ b/src/quafu/benchmark/unitary_test.py @@ -0,0 +1,11 @@ +from scipy.stats import unitary_group +from quafu import QuantumCircuit + +nqubit = 5 +qubits = list(range(nqubit)) +U0 = unitary_group.rvs(2 ** nqubit) + +# Using QSD to decompose the unitary +qc = QuantumCircuit(nqubit) +qc.unitary(U0, qubits) +qc.draw_circuit() diff --git a/src/quafu/benchmark/variational_n4.py b/src/quafu/benchmark/variational_n4.py new file mode 100644 index 0000000..f6213a6 --- /dev/null +++ b/src/quafu/benchmark/variational_n4.py @@ -0,0 +1,91 @@ +from quafu import QuantumCircuit + +n = 4 +qc = QuantumCircuit(n) + +qasm = """ +// Generated from Cirq v0.8.0 + +OPENQASM 2.0; +include "qelib1.inc"; + +// Qubits: [0, 1, 2, 3] +qreg q[4]; +creg c[4]; + +x q[0]; +x q[1]; + +// Gate: PhasedISWAP**0.9951774602384953 +rz(pi*0.25) q[1]; +rz(pi*-0.25) q[2]; +cx q[1],q[2]; +h q[1]; +cx q[2],q[1]; +rz(pi*0.4975887301) q[1]; +cx q[2],q[1]; +rz(pi*-0.4975887301) q[1]; +h q[1]; +cx q[1],q[2]; +rz(pi*-0.25) q[1]; +rz(pi*0.25) q[2]; + +rz(0) q[2]; + +// Gate: PhasedISWAP**-0.5024296754026449 +rz(pi*0.25) q[0]; +rz(pi*-0.25) q[1]; +cx q[0],q[1]; +h q[0]; +cx q[1],q[0]; +rz(pi*-0.2512148377) q[0]; +cx q[1],q[0]; +rz(pi*0.2512148377) q[0]; +h q[0]; +cx q[0],q[1]; +rz(pi*-0.25) q[0]; +rz(pi*0.25) q[1]; + +rz(0) q[1]; + +// Gate: PhasedISWAP**-0.49760685888033646 +rz(pi*0.25) q[2]; +rz(pi*-0.25) q[3]; +cx q[2],q[3]; +h q[2]; +cx q[3],q[2]; +rz(pi*-0.2488034294) q[2]; +cx q[3],q[2]; +rz(pi*0.2488034294) q[2]; +h q[2]; +cx q[2],q[3]; +rz(pi*-0.25) q[2]; +rz(pi*0.25) q[3]; + +rz(0) q[3]; + +// Gate: PhasedISWAP**0.004822678143889672 +rz(pi*0.25) q[1]; +rz(pi*-0.25) q[2]; +cx q[1],q[2]; +h q[1]; +cx q[2],q[1]; +rz(pi*0.0024113391) q[1]; +cx q[2],q[1]; +rz(pi*-0.0024113391) q[1]; +h q[1]; +cx q[1],q[2]; +rz(pi*-0.25) q[1]; +rz(pi*0.25) q[2]; + +rz(0) q[2]; + +measure q[0] -> c[0]; +measure q[1] -> c[1]; +measure q[2] -> c[2]; +measure q[3] -> c[3]; +""" + +qc.from_openqasm(qasm) +qc.draw_circuit() +qc.plot_circuit(show=True, title='Variational n4') diff --git a/src/quafu/circuits/para_circuit.py b/src/quafu/circuits/para_circuit.py index 8a31b96..04a85f1 100644 --- a/src/quafu/circuits/para_circuit.py +++ b/src/quafu/circuits/para_circuit.py @@ -1,8 +1,7 @@ - - from .quantum_circuit import QuantumCircuit from ..transpiler.Qcovercompiler import QcoverCompiler + class QAOACircuit(QuantumCircuit): def __init__(self, logical_qubits, physical_qubits, nodes, edges, params, p, gate="CNOT"): num = logical_qubits @@ -20,12 +19,9 @@ def compile_to_IOP(self): QASM from qcover directly """ qaoa_compiler = QcoverCompiler() - self.qasm = qaoa_compiler.graph_to_qasm(self.logical_qubits, self.physical_qubits, self.nodes, self.edges, self.paras, self.p, gate=self.gate) + self.qasm = qaoa_compiler.graph_to_qasm(self.logical_qubits, self.physical_qubits, self.nodes, self.edges, + self.paras, self.p, gate=self.gate) def upate_paras(self, paras): self.paras = paras pass - - - - diff --git a/src/quafu/circuits/quantum_circuit.py b/src/quafu/circuits/quantum_circuit.py index 827587f..db4c1f5 100644 --- a/src/quafu/circuits/quantum_circuit.py +++ b/src/quafu/circuits/quantum_circuit.py @@ -1,22 +1,28 @@ -from typing import Iterable +from typing import List import numpy as np -from ..elements.quantum_element import Barrier, Delay, ControlledGate, MultiQubitGate, ParaMultiQubitGate, QuantumGate, SingleQubitGate, XYResonance -from ..elements.element_gates import * + +import quafu.elements.element_gates.clifford +import quafu.elements.element_gates.pauli +from quafu.elements.quantum_element.pulses.quantum_pulse import QuantumPulse +from ..elements.quantum_element import Barrier, Delay, MultiQubitGate, QuantumGate, ControlledGate, \ + SingleQubitGate, XYResonance +import quafu.elements.element_gates as qeg from ..exceptions import CircuitError + class QuantumCircuit(object): def __init__(self, num: int): """ Initialize a QuantumCircuit object - - Args: + + Args: num (int): Total qubit number used """ self.num = num self.gates = [] self.openqasm = "" self.circuit = [] - self.measures = dict(zip(range(num), range(num))) + self.measures = {} self._used_qubits = [] @property @@ -24,14 +30,17 @@ def used_qubits(self) -> List: self.layered_circuit() return self._used_qubits - def add_gate(self, gate : QuantumGate): + def add_gate(self, gate: QuantumGate): + pos = np.array(gate.pos) + if np.any(pos >= self.num): + raise CircuitError(f"Gate position out of range: {gate.pos}") self.gates.append(gate) def layered_circuit(self) -> np.ndarray: """ Make layered circuit from the gate sequence self.gates. - Returns: + Returns: A layered list with left justed circuit. """ num = self.num @@ -39,7 +48,7 @@ def layered_circuit(self) -> np.ndarray: gateQlist = [[] for i in range(num)] used_qubits = [] for gate in gatelist: - if isinstance(gate, SingleQubitGate) or isinstance(gate, Delay): + if isinstance(gate, SingleQubitGate) or isinstance(gate, Delay) or isinstance(gate, QuantumPulse): gateQlist[gate.pos].append(gate) if gate.pos not in used_qubits: used_qubits.append(gate.pos) @@ -86,7 +95,7 @@ def layered_circuit(self) -> np.ndarray: self._used_qubits = list(used_qubits) return self.circuit - def draw_circuit(self, width : int=4, return_str : bool=False): + def draw_circuit(self, width: int = 4, return_str: bool = False): """ Draw layered circuit using ASCII, print in terminal. @@ -107,7 +116,7 @@ def draw_circuit(self, width : int=4, return_str : bool=False): maxlen = 1 + width for i in range(num): gate = layergates[i] - if isinstance(gate, SingleQubitGate) or isinstance(gate, Delay): + if isinstance(gate, SingleQubitGate) or isinstance(gate, Delay) or (isinstance(gate, QuantumPulse)): printlist[i * 2, l] = gate.symbol maxlen = max(maxlen, len(gate.symbol) + width) @@ -117,10 +126,10 @@ def draw_circuit(self, width : int=4, return_str : bool=False): printlist[2 * q1 + 1:2 * q2, l] = "|" printlist[q1 * 2, l] = "#" printlist[q2 * 2, l] = "#" - if isinstance(gate, ControlledGate): #Controlled-Multiqubit gate + if isinstance(gate, ControlledGate): # Controlled-Multiqubit gate for ctrl in gate.ctrls: printlist[reduce_map[ctrl] * 2, l] = "*" - + if gate.targ_name == "SWAP": printlist[reduce_map[gate.targs[0]] * 2, l] = "x" printlist[reduce_map[gate.targs[1]] * 2, l] = "x" @@ -134,8 +143,8 @@ def draw_circuit(self, width : int=4, return_str : bool=False): else: printlist[tq1 + tq2, l] = gate.symbol maxlen = max(maxlen, len(gate.symbol) + width) - - else: #Multiqubit gate + + else: # Multiqubit gate if gate.name == "SWAP": printlist[q1 * 2, l] = "x" printlist[q2 * 2, l] = "x" @@ -150,7 +159,6 @@ def draw_circuit(self, width : int=4, return_str : bool=False): q2 = reduce_map[max(pos)] printlist[2 * q1:2 * q2 + 1, l] = "||" maxlen = max(maxlen, len("||")) - printlist[-1, l] = maxlen @@ -166,13 +174,18 @@ def draw_circuit(self, width : int=4, return_str : bool=False): circuitstr.append("".ljust(6) + "".join( [printlist[j, l].center(int(printlist[-1, l]), " ") for l in range(depth)])) circuitstr = "\n".join(circuitstr) - + if return_str: return circuitstr else: print(circuitstr) - def from_openqasm(self, openqasm : str): + def plot_circuit(self, *args, **kwargs): + from quafu.visualisation.circuitPlot import CircuitPlotManager + cmp = CircuitPlotManager(self) + return cmp(*args, **kwargs) + + def from_openqasm(self, openqasm: str): """ Initialize the circuit from openqasm text. Args: @@ -310,50 +323,14 @@ def to_openqasm(self) -> str: """ Convert the circuit to openqasm text. - Returns: + Returns: openqasm text. """ qasm = "OPENQASM 2.0;\ninclude \"qelib1.inc\";\n" qasm += "qreg q[%d];\n" % self.num qasm += "creg meas[%d];\n" % len(self.measures) for gate in self.gates: - if isinstance(gate, SingleQubitGate): - if isinstance(gate, ParaSingleQubitGate): - if isinstance(gate.paras, Iterable): - qasm += "%s(" %gate.name.lower() + ",".join(["%s" %para for para in gate.paras]) + ") q[%d];\n" % (gate.pos) - else: - qasm += "%s(%s) " %(gate.name.lower(), gate.paras) + "q[%d];\n" % (gate.pos) - else: - if gate.name == "SY": - qasm += "ry(pi/2) q[%d];\n" %(gate.pos) - elif gate.name == "W": - qasm += "rz(-pi/4) q[%d];\nrx(pi) q[%d];\nrz(pi/4) q[%d];\n" %(gate.pos, gate.pos, gate.pos) - elif gate.name == "SW": - qasm += "rz(-pi/4) q[%d];\nrx(pi/2) q[%d];\nrz(pi/4) q[%d];\n" %(gate.pos, gate.pos, gate.pos) - else: - qasm += "%s q[%d];\n" % (gate.name.lower(), gate.pos) - - elif isinstance(gate, Delay): - qasm += "delay(%d%s) q[%d];\n" % (gate.duration, gate.unit, gate.pos) - elif isinstance(gate, XYResonance): - qasm += "xy(%d%s) " %(gate.duration, gate.unit) + ",".join(["q[%d]" % p for p in range(min(gate.pos), max(gate.pos)+1)]) + ";\n" - - elif isinstance(gate, Barrier) or isinstance(gate, MultiQubitGate): - if isinstance(gate, ParaMultiQubitGate) or (isinstance(gate, ControlledGate) and bool(gate.paras)): - if isinstance(gate.paras, Iterable): - qasm += "%s(" %gate.name.lower() + ",".join(["%s" %para for para in gate.paras]) + ") " + ",".join(["q[%d]" % p for p in gate.pos]) + ";\n" - else: - qasm += "%s(%s) " %(gate.name.lower(), gate.paras) + ",".join(["q[%d]" % p for p in gate.pos]) + ";\n" - - else: - if gate.name == "CS": - qasm += "cp(pi/2) " + "q[%d],q[%d];\n" % (gate.pos[0], gate.pos[1]) - elif gate.name == "CT": - qasm += "cp(pi/4) " + "q[%d],q[%d];\n" % (gate.pos[0], gate.pos[1]) - elif gate.name == "barrier": - qasm += "barrier " + ",".join(["q[%d]" % p for p in range(min(gate.pos), max(gate.pos)+1)]) + ";\n" - else: - qasm += "%s " %(gate.name.lower()) + ",".join(["q[%d]" % p for p in gate.pos]) + ";\n" + qasm += gate.to_qasm() + ";\n" for key in self.measures: qasm += "measure q[%d] -> meas[%d];\n" % (key, self.measures[key]) @@ -361,15 +338,15 @@ def to_openqasm(self) -> str: self.openqasm = qasm return qasm - def id(self, pos: int) -> "QuantumCircuit": """ Identity gate. Args: pos (int): qubit the gate act. - """ - self.gates.append(IdGate(pos)) + """ + gate = qeg.IdGate(pos) + self.add_gate(gate) return self def h(self, pos: int) -> "QuantumCircuit": @@ -379,7 +356,8 @@ def h(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(HGate(pos)) + gate = quafu.elements.element_gates.clifford.HGate(pos) + self.add_gate(gate) return self def x(self, pos: int) -> "QuantumCircuit": @@ -389,7 +367,8 @@ def x(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(XGate(pos)) + gate = qeg.XGate(pos) + self.add_gate(gate) return self def y(self, pos: int) -> "QuantumCircuit": @@ -399,7 +378,8 @@ def y(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(YGate(pos)) + gate = qeg.YGate(pos) + self.add_gate(gate) return self def z(self, pos: int) -> "QuantumCircuit": @@ -409,7 +389,7 @@ def z(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(ZGate(pos)) + self.add_gate(qeg.ZGate(pos)) return self def t(self, pos: int) -> "QuantumCircuit": @@ -419,9 +399,9 @@ def t(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(TGate(pos)) + self.add_gate(quafu.elements.element_gates.clifford.TGate(pos)) return self - + def tdg(self, pos: int) -> "QuantumCircuit": """ Tdg gate. (Inverse of T gate) @@ -429,7 +409,8 @@ def tdg(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(TdgGate(pos)) + self.add_gate(quafu.elements.element_gates.clifford.TdgGate(pos)) + return self def s(self, pos: int) -> "QuantumCircuit": """ @@ -438,7 +419,7 @@ def s(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(SGate(pos)) + self.add_gate(quafu.elements.element_gates.clifford.SGate(pos)) return self def sdg(self, pos: int) -> "QuantumCircuit": @@ -448,7 +429,7 @@ def sdg(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(SdgGate(pos)) + self.add_gate(quafu.elements.element_gates.clifford.SdgGate(pos)) return self def sx(self, pos: int) -> "QuantumCircuit": @@ -458,7 +439,18 @@ def sx(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(SXGate(pos)) + self.add_gate(quafu.elements.element_gates.pauli.SXGate(pos)) + return self + + def sxdg(self, pos: int) -> "QuantumCircuit": + """ + Inverse of √X gate. + + Args: + pos (int): qubit the gate act. + """ + gate = quafu.elements.element_gates.pauli.SXdgGate(pos) + self.add_gate(gate) return self def sy(self, pos: int) -> "QuantumCircuit": @@ -468,7 +460,18 @@ def sy(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(SYGate(pos)) + self.add_gate(quafu.elements.element_gates.pauli.SYGate(pos)) + return self + + def sydg(self, pos: int) -> "QuantumCircuit": + """ + Inverse of √Y gate. + + Args: + pos (int): qubit the gate act. + """ + gate = quafu.elements.element_gates.pauli.SYdgGate(pos) + self.add_gate(gate) return self def w(self, pos: int) -> "QuantumCircuit": @@ -478,9 +481,9 @@ def w(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(WGate(pos)) + self.add_gate(qeg.WGate(pos)) return self - + def sw(self, pos: int) -> "QuantumCircuit": """ √W gate. @@ -488,9 +491,9 @@ def sw(self, pos: int) -> "QuantumCircuit": Args: pos (int): qubit the gate act. """ - self.gates.append(SWGate(pos)) + self.add_gate(qeg.SWGate(pos)) return self - + def rx(self, pos: int, para: float) -> "QuantumCircuit": """ Single qubit rotation Rx gate. @@ -499,52 +502,56 @@ def rx(self, pos: int, para: float) -> "QuantumCircuit": pos (int): qubit the gate act. para (float): rotation angle """ - self.gates.append(RXGate(pos, para)) + self.add_gate(qeg.RXGate(pos, para)) return self def ry(self, pos: int, para: float) -> "QuantumCircuit": """ Single qubit rotation Ry gate. - + Args: pos (int): qubit the gate act. para (float): rotation angle """ - self.gates.append(RYGate(pos, para)) + self.add_gate(qeg.RYGate(pos, para)) return self def rz(self, pos: int, para: float) -> "QuantumCircuit": """ Single qubit rotation Rz gate. - + Args: pos (int): qubit the gate act. para (float): rotation angle """ - self.gates.append(RZGate(pos, para)) + self.add_gate(qeg.RZGate(pos, para)) return self def p(self, pos: int, para: float) -> "QuantumCircuit": """ Phase gate - + Args: pos (int): qubit the gate act. para (float): rotation angle """ - self.gates.append(PhaseGate(pos, para)) + self.add_gate(qeg.PhaseGate(pos, para)) + return self def cnot(self, ctrl: int, tar: int) -> "QuantumCircuit": """ CNOT gate. - + Args: ctrl (int): control qubit. tar (int): target qubit. """ - self.gates.append(CXGate(ctrl, tar)) + self.add_gate(qeg.CXGate(ctrl, tar)) return self + def cx(self, ctrl: int, tar: int) -> "QuantumCircuit": + return self.cnot(ctrl=ctrl, tar=tar) + def cy(self, ctrl: int, tar: int) -> "QuantumCircuit": """ Control-Y gate. @@ -553,18 +560,18 @@ def cy(self, ctrl: int, tar: int) -> "QuantumCircuit": ctrl (int): control qubit. tar (int): target qubit. """ - self.gates.append(CYGate(ctrl, tar)) + self.add_gate(qeg.CYGate(ctrl, tar)) return self def cz(self, ctrl: int, tar: int) -> "QuantumCircuit": """ Control-Z gate. - + Args: ctrl (int): control qubit. tar (int): target qubit. """ - self.gates.append(CZGate(ctrl, tar)) + self.add_gate(qeg.CZGate(ctrl, tar)) return self def cs(self, ctrl: int, tar: int) -> "QuantumCircuit": @@ -574,7 +581,7 @@ def cs(self, ctrl: int, tar: int) -> "QuantumCircuit": ctrl (int): control qubit. tar (int): target qubit. """ - self.gates.append(CSGate(ctrl, tar)) + self.add_gate(qeg.CSGate(ctrl, tar)) return self def ct(self, ctrl: int, tar: int) -> "QuantumCircuit": @@ -584,8 +591,8 @@ def ct(self, ctrl: int, tar: int) -> "QuantumCircuit": ctrl (int): control qubit. tar (int): target qubit. """ - - self.gates.append(CTGate(ctrl, tar)) + + self.add_gate(qeg.CTGate(ctrl, tar)) return self def cp(self, ctrl: int, tar: int, para) -> "QuantumCircuit": @@ -596,19 +603,29 @@ def cp(self, ctrl: int, tar: int, para) -> "QuantumCircuit": ctrl (int): control qubit. tar (int): target qubit. """ - self.gates.append(CPGate(ctrl, tar, para)) + self.add_gate(qeg.CPGate(ctrl, tar, para)) return self - def swap(self, q1: int, q2: int) -> "QuantumCircuit": """ SWAP gate - + Args: q1 (int): qubit the gate act. q2 (int): qubit the gate act. """ - self.gates.append(SwapGate(q1, q2)) + self.add_gate(qeg.SwapGate(q1, q2)) + return self + + def iswap(self, q1: int, q2: int) -> "QuantumCircuit": + """ + iSWAP gate + + Args: + q1 (int): qubit the gate act. + q2 (int): qubit the gate act. + """ + self.add_gate(qeg.ISwapGate(q1, q2)) return self def toffoli(self, ctrl1: int, ctrl2: int, targ: int) -> "QuantumCircuit": @@ -620,29 +637,31 @@ def toffoli(self, ctrl1: int, ctrl2: int, targ: int) -> "QuantumCircuit": ctrl2 (int): control qubit targ (int): target qubit """ - self.gates.append(ToffoliGate(ctrl1, ctrl2, targ)) + self.add_gate(qeg.ToffoliGate(ctrl1, ctrl2, targ)) return self - - def fredkin(self, ctrl: int, targ1:int , targ2: int) -> "QuantumCircuit": + + def fredkin(self, ctrl: int, targ1: int, targ2: int) -> "QuantumCircuit": """ Fredkin gate - + Args: ctrl (int): control qubit targ1 (int): target qubit targ2 (int): target qubit """ - self.gates.append(FredkinGate(ctrl, targ1, targ2)) + self.add_gate(qeg.FredkinGate(ctrl, targ1, targ2)) return self - def barrier(self, qlist: List[int]) -> "QuantumCircuit": + def barrier(self, qlist: List[int] = None) -> "QuantumCircuit": """ Add barrier for qubits in qlist. - + Args: qlist (list[int]): A list contain the qubit need add barrier. When qlist contain at least two qubit, the barrier will be added from minimum qubit to maximum qubit. For example: barrier([0, 2]) create barrier for qubits 0, 1, 2. To create discrete barrier, using barrier([0]), barrier([2]). """ - self.gates.append(Barrier(qlist)) + if qlist is None: + qlist = list(range(self.num)) + self.add_gate(Barrier(qlist)) return self def delay(self, pos, duration, unit="ns") -> "QuantumCircuit": @@ -652,12 +671,12 @@ def delay(self, pos, duration, unit="ns") -> "QuantumCircuit": Args: pos (int): qubit need delay. duration (int): duration of qubit delay, which represents integer times of unit. - unit (str): time unit for the duration. Can be "ns" and "us". + unit (str): time unit for the duration. Can be "ns" and "us". """ - self.gates.append(Delay(pos, duration, unit=unit)) + self.add_gate(Delay(pos, duration, unit=unit)) return self - def xy(self, qs: int, qe: int, duration: int, unit: str="ns") -> "QuantumCircuit": + def xy(self, qs: int, qe: int, duration: int, unit: str = "ns") -> "QuantumCircuit": """ XY resonance time evolution for quantum simulator Args: @@ -667,7 +686,7 @@ def xy(self, qs: int, qe: int, duration: int, unit: str="ns") -> "QuantumCircuit unit: time unit of duration. """ - self.gates.append(XYResonance(qs, qe, duration, unit=unit)) + self.add_gate(XYResonance(qs, qe, duration, unit=unit)) return self def rxx(self, q1: int, q2: int, theta): @@ -679,8 +698,8 @@ def rxx(self, q1: int, q2: int, theta): theta: rotation angle. """ - self.gates.append(RXXGate(q1, q2, theta)) - + self.add_gate(qeg.RXXGate(q1, q2, theta)) + def ryy(self, q1: int, q2: int, theta): """ Rotation about 2-qubit YY axis. @@ -690,7 +709,7 @@ def ryy(self, q1: int, q2: int, theta): theta: rotation angle. """ - self.gates.append(RYYGate(q1, q2, theta)) + self.add_gate(qeg.RYYGate(q1, q2, theta)) def rzz(self, q1: int, q2: int, theta): """ @@ -701,8 +720,8 @@ def rzz(self, q1: int, q2: int, theta): theta: rotation angle. """ - self.gates.append(RZZGate(q1, q2, theta)) - + self.add_gate(qeg.RZZGate(q1, q2, theta)) + def mcx(self, ctrls: List[int], targ: int): """ Multi-controlled X gate. @@ -710,8 +729,8 @@ def mcx(self, ctrls: List[int], targ: int): ctrls: A list of control qubits. targ: Target qubits. """ - self.gates.append(MCXGate(ctrls, targ)) - + self.add_gate(qeg.MCXGate(ctrls, targ)) + def mcy(self, ctrls: List[int], targ: int): """ Multi-controlled Y gate. @@ -719,8 +738,8 @@ def mcy(self, ctrls: List[int], targ: int): ctrls: A list of control qubits. targ: Target qubits. """ - self.gates.append(MCYGate(ctrls, targ)) - + self.add_gate(qeg.MCYGate(ctrls, targ)) + def mcz(self, ctrls: List[int], targ: int): """ Multi-controlled Z gate. @@ -728,25 +747,48 @@ def mcz(self, ctrls: List[int], targ: int): ctrls: A list of control qubits. targ: Target qubits. """ - self.gates.append(MCZGate(ctrls, targ)) - + self.add_gate(qeg.MCZGate(ctrls, targ)) + + def unitary(self, matrix: np.ndarray, pos: List[int]): + """ + Apply unitary to circuit on specified qubits. + + Args: + matrix (np.ndarray): unitary matrix. + pos (list[int]): qubits the gate act on. + """ + compiler = qeg.UnitaryDecomposer(array=matrix, qubits=pos) + compiler.apply_to_qc(self) - def measure(self, pos: List[int], cbits: List[int] = []) -> None: + def measure(self, pos: List[int] = None, cbits: List[int] = None) -> None: """ Measurement setting for experiment device. - + Args: pos: Qubits need measure. cbits: Classical bits keeping the measure results. """ - - self.measures = dict(zip(pos, range(len(pos)))) + if pos is None: + pos = list(range(self.num)) if cbits: - if len(cbits) == len(self.measures): - self.measures = dict(zip(pos, cbits)) - else: + if not len(cbits) == len(pos): raise CircuitError("Number of measured bits should equal to the number of classical bits") + else: + cbits = pos + newly_measures = dict(zip(pos, cbits)) + self.measures = {**self.measures, **newly_measures} + if not len(self.measures.values()) == len(set(self.measures.values())): + raise ValueError("Measured bits not uniquely assigned.") - + def add_pulse(self, + pulse: QuantumPulse, + pos: int = None) -> "QuantumCircuit": + """ + Add quantum gate from pulse. + """ + if pos is not None: + pulse.set_pos(pos) + self.add_gate(pulse) + return self diff --git a/src/quafu/dagcircuits/__init__.py b/src/quafu/dagcircuits/__init__.py new file mode 100755 index 0000000..8b13789 --- /dev/null +++ b/src/quafu/dagcircuits/__init__.py @@ -0,0 +1 @@ + diff --git a/src/quafu/dagcircuits/circuit_dag.py b/src/quafu/dagcircuits/circuit_dag.py new file mode 100644 index 0000000..172472c --- /dev/null +++ b/src/quafu/dagcircuits/circuit_dag.py @@ -0,0 +1,473 @@ +import numpy as np +from quafu import QuantumCircuit + +from quafu.elements.element_gates import * +from quafu.elements.quantum_element import Barrier, Delay, Measure, XYResonance +from quafu.pulses.quantum_pulse import GaussianPulse, RectPulse, FlattopPulse + +import networkx as nx +from typing import Dict, Any, List, Union +import copy + +from quafu.dagcircuits.instruction_node import InstructionNode # instruction_node.py in the same folder as circuit_dag.py now +from dag_circuit import DAGCircuit # dag_circuit.py in the same folder as circuit_dag.py now + +# import pygraphviz as pgv +from networkx.drawing.nx_pydot import write_dot +from IPython.display import Image, SVG + + + +# transform a gate in quantumcircuit of quafu(not include measure_gate), +# into a node in the graph, with specific label. +def gate_to_node(input_gate,specific_label: str): + ''' + transform a gate in quantumcircuit of quafu(not include measure_gate), + into a node in the graph, with specific label. + + Args: + inputgate: a gate in quantumcircuit of quafu(not include measure_gate) + label: the label of the node in the graph + + Returns: + node: a node in the graph, with specific label. A GateWrapper object + + ''' + + import copy + gate = copy.deepcopy(input_gate) # avoid modifying the original gate + if not isinstance(gate.pos, list): # if gate.pos is not a list, make it a list + gate.pos = [gate.pos] + + # use getattr check 'paras' and other attributes if exist. if the attr doesn't exist,return None + gate.paras = getattr(gate, 'paras', None) or None + gate.duration = getattr(gate, 'duration', None) or None + gate.unit = getattr(gate, 'unit', None) or None + gate.channel = getattr(gate, 'channel', None) or None + gate.time_func = getattr(gate, 'time_func', None) or None + + if gate.paras and not isinstance(gate.paras, list): # if paras is True and not a list, make it a list + gate.paras = [gate.paras] + + # hashable_gate = InstructionNode(gate.name, gate.pos, gate.paras,gate.matrix,gate.duration,gate.unit, label=i) + hashable_gate = InstructionNode(gate.name, gate.pos, gate.paras,gate.duration, gate.unit,gate.channel,gate.time_func, label=specific_label) + return hashable_gate + +# Building a DAG Graph using DAGCircuit from a QuantumCircuit +def circuit_to_dag(circuit, measure_flag = True): + ''' + Building a DAG Graph using DAGCircui from a QuantumCircuit + + Args: + circuit: a QuantumCircuit object + measure_flag: whether to add measure_gate node to the dag graph + + Returns: + g: a DAGCircuit object + + example: + .. jupyter-execute:: + + from circuit_dag import circuit_to_dag, dag_to_circuit, draw_dag + from quafu import QuantumCircuit + + # Create a quantum circuit as an example that you can modify as needed + circuit = QuantumCircuit(2) + circuit.h(0) + circuit.cnot(0, 1) + + # Build the dag graph + dep_graph = circuit_to_dag(circuit) # dag graph + ''' + + # Starting Label Index + i = 0 + + # A dictionary to store the last use of any qubit + qubit_last_use = {} + + # g = nx.MultiDiGraph() # two nodes can have multiple edges + # g = nx.DiGraph() # two nodes can only have one edge + g = DAGCircuit() # two nodes can only have one edge + + # Add the start node + # g.add_node(-1,{"color": "green"}) + g.add_nodes_from([(-1, {"color": "green"})]) + + # deepcopy the circuit to avoid modifying the original circuit + # gates = copy.deepcopy(circuit.gates) # need to import copy + # change to: gate = copy.deepcopy(input_gate) in gate_to_node() + + for gate in circuit.gates: + # transform gate to node + hashable_gate = gate_to_node(gate,specific_label=i) + i += 1 + + g.add_node(hashable_gate,color="blue") + + # Add edges based on qubit_last_use; update last use + for qubit in hashable_gate.pos: + if qubit in qubit_last_use: + g.add_edge(qubit_last_use[qubit], hashable_gate,label=f'q{qubit}') + else: + g.add_edge(-1, hashable_gate,label=f'q{qubit}',color="green") + + qubit_last_use[qubit] = hashable_gate + + if measure_flag: + # Add measure_gate node + qm = Any + qm.name = "measure" + qm.paras, qm.duration, qm.unit = [None,None,None] + qm.channel, qm.time_func = [None,None] + qm.pos = copy.deepcopy(circuit.measures) # circuit.measures is a dict + measure_gate = InstructionNode(qm.name, qm.pos, qm.paras, qm.duration, qm.unit, qm.channel, qm.time_func, label="m") + g.add_node(measure_gate,color="blue") + # Add edges from qubit_last_use[qubit] to measure_gate + for qubit in measure_gate.pos: + if qubit in qubit_last_use: + g.add_edge(qubit_last_use[qubit], measure_gate,label=f'q{qubit}') + else: + g.add_edge(-1, measure_gate,label=f'q{qubit}',color="green") + + qubit_last_use[qubit] = measure_gate + + # Add the end node + # g.add_node(float('inf'),{"color": "red"}) + g.add_nodes_from([(float('inf'), {"color": "red"})]) + + for qubit in qubit_last_use: + g.add_edge(qubit_last_use[qubit], float('inf'),label=f'q{qubit}',color="red") + + # update qubits_used, cbits_used, num_instruction_nodes + g.update_qubits_used() + g.update_cbits_used() + g.update_num_instruction_nodes() + + return g + + +# transform gate in dag nodes to gate in circuit which can be added to circuit +gate_classes = { + "x": XGate, + "y": YGate, + "z": ZGate, + "h": HGate, + "s": SGate, + "sdg": SdgGate, + "t": TGate, + "tdg": TdgGate, + "rx": RXGate, + "ry": RYGate, + "rz": RZGate, + "id": IdGate, + "sx": SXGate, + "sy": SYGate, + "w": WGate, + "sw": SWGate, + "p": PhaseGate, + "delay": Delay, + "barrier": Barrier, + "cx": CXGate, + "cnot": CXGate, + "cp": CPGate, + "swap": SwapGate, + "rxx": RXXGate, + "ryy": RYYGate, + "rzz": RZZGate, + "cy": CYGate, + "cz": CZGate, + "cs": CSGate, + "ct": CTGate, + "xy": XYResonance, + "ccx": ToffoliGate, + "toffoli": ToffoliGate, + "cswap": FredkinGate, + "fredkin": FredkinGate, + "mcx": MCXGate, + "mcy": MCYGate, + "mcz": MCZGate, + "gaussian": GaussianPulse, + "rect" : RectPulse, + "flattop" : FlattopPulse, + "measure" : Measure +} + +def node_to_gate(gate_in_dag): + """ + transform gate in dag graph, to gate in circuit which can be added to circuit + + Args: + gate_in_dag: a node in dag graph , gate_in_dag is a GateWrapper object. + in instruction_node, gate_in_dag.name is uppercase, gate_in_dag.pos is a list or a dict + gate_transform support gate with one qubit or more qubits, not measures! + and you should exculde nodes [-1 ,float('inf') , measure_gate] in dag graph + + Returns: + gate: gate which can be added to circuit in quafu + + example: + import networkx as nx + from quafu import QuantumCircuit + qcircuit = QuantumCircuit(n) + + for gate in nx.topological_sort(dep_graph): + + if gate not in [-1, float('inf')]: + # measure gate to do + if gate.name == "measure": + qcircuit.measures = gate.pos + + else: + # use gate_transform to transform gate in dag graph to gate in circuit + qcircuit.gates.append(node_to_gate(gate)) + return qcircuit + + """ + + gate_name = gate_in_dag.name.lower() + gate_class = gate_classes.get(gate_name) + + if not gate_class: + raise ValueError("gate is not supported") + + if gate_name == "barrier": + return gate_class(gate_in_dag.pos) + + # 从gate_in_dag获取参数列表 + args = gate_in_dag.pos + if gate_in_dag.paras: + args += gate_in_dag.paras + + # 处理 gate.duration 和 gate.unit + if gate_name in ["delay", "xy"]: + args.append(gate_in_dag.duration) + args.append(gate_in_dag.unit) + + # 处理多量子比特门 + if gate_name in ["mcx", "mcy", "mcz"]: + control_qubits = gate_in_dag.pos[:-1] + target_qubit = gate_in_dag.pos[-1] + return gate_class(control_qubits, target_qubit) + + # pulse node + if gate_name in ["gaussian", "rect", "flattop"]: + duration = gate_in_dag.duration + unit = gate_in_dag.unit + paras = gate_in_dag.paras + pos = gate_in_dag.pos[0] + channel = gate_in_dag.channel + # time_func = gate_in_dag.time_func + + return gate_class(pos, *paras, duration, unit, channel) + + # measure node + if gate_name == "measure": + return gate_class(gate_in_dag.pos) + + return gate_class(*args) + + + +# From DAG with Hashable Gates to quafu Gates added to circuit +def dag_to_circuit(dep_graph, n: int): + ''' + From DAG with Hashable Gates to quafu Gates added to circuit + + Args: + dep_graph (DAG): DAG with Hashable Gates + n (int): number of qubits + + Returns: + qcircuit (QuantumCircuit): quafu QuantumCircuit + + example: + .. jupyter-execute:: + + from circuit_dag import circuit_to_dag, dag_to_circuit, draw_dag + from quafu import QuantumCircuit + + # Create a quantum circuit as an example that you can modify as needed + circuit = QuantumCircuit(2) + circuit.h(0) + circuit.cnot(0, 1) + + # Build the dag graph + dep_graph = circuit_to_dag(circuit) # dag graph + + # use dag_to_circuit to transform dag graph to a new circuit + reconstructed_circuit = dag_to_circuit(dep_graph, circuit.num) + + + ''' + + qcircuit = QuantumCircuit(n) + + for gate in nx.topological_sort(dep_graph): + + if gate not in [-1, float('inf')]: + # measure gate to do + if gate.name == "measure": + qcircuit.measures = gate.pos + + else: + # use gate_transform to transform gate in dag graph to gate in circuit + qcircuit.gates.append(node_to_gate(gate)) + return qcircuit + +# Helper function to visualize the DAG,check the example in the docstring +def draw_dag(dep_g, output_format="png"): + ''' + Helper function to visualize the DAG + + Args: + dep_g (DAG): DAG with Hashable Gates + output_format (str): output format, "png" or "svg" + + Returns: + img (Image or SVG): show the image of DAG, which is Image(filename="dag.png") or SVG(filename="dag.svg") + + example: + .. jupyter-execute:: + ex1: + # directly draw PNG picture + draw_dag(dep_g, output_format="png") # save a png picture "dag.png" and show it in jupyter notebook + + # directly draw SVG picture + draw_dag(dep_g, output_format="svg") # save a svg picture "dag.svg" and show it in jupyter notebook + + ex2: + # generate PNG picture + img_png = draw_dag(dep_g, output_format="png") + + # generate SVG picture + img_svg = draw_dag(dep_g, output_format="svg") + + # show PNG picture + img_png + + # show SVG picture + img_svg + + + ''' + import pygraphviz + write_dot(dep_g, "dag.dot") + G = pygraphviz.AGraph("dag.dot") + G.layout(prog="dot") + + if output_format == "png": + G.draw("dag.png") + return Image(filename="dag.png") + elif output_format == "svg": + G.draw("dag.svg") + return SVG(filename="dag.svg") + else: + raise ValueError("Unsupported output format: choose either 'png' or 'svg'") + + +def nodelist_to_dag(op_nodes: List[Any]) -> DAGCircuit: + # Starting Label Index + i = 0 + + # A dictionary to store the last use of any qubit + qubit_last_use = {} + + # g = nx.MultiDiGraph() # two nodes can have multiple edges + # g = nx.DiGraph() # two nodes can only have one edge + g = DAGCircuit() + + # Add the start node + # g.add_node(-1,{"color": "green"}) + g.add_nodes_from([(-1, {"color": "green"})]) + + # deepcopy the circuit to avoid modifying the original circuit + # gates = copy.deepcopy(circuit.gates) # need to import copy + # change to: gate = copy.deepcopy(input_gate) in gate_to_node() + + for op_node in op_nodes: + # transform gate to node + hashable_gate = copy.deepcopy(op_node) + g.add_node(hashable_gate,color="blue") + + # Add edges based on qubit_last_use; update last use + for qubit in hashable_gate.pos: + if qubit in qubit_last_use: + g.add_edge(qubit_last_use[qubit], hashable_gate,label=f'q{qubit}') + else: + g.add_edge(-1, hashable_gate,label=f'q{qubit}',color="green") + + qubit_last_use[qubit] = hashable_gate + + + # Add the end node + # g.add_node(float('inf'),{"color": "red"}) + g.add_nodes_from([(float('inf'), {"color": "red"})]) + + for qubit in qubit_last_use: + g.add_edge(qubit_last_use[qubit], float('inf'),label=f'q{qubit}',color="red") + + # update the qubits_used, cbits_used, num_instruction_nodes + g.qubits_used = g.update_qubits_used() + g.cbits_used = g.update_cbits_used() + g.num_instruction_nodes = g.update_num_instruction_nodes() + + return g + +# nodes_qubit_mapping_dict +def nodelist_qubit_mapping_dict(nodes_list): + ''' + Args: + nodes_list: a list of nodes + Returns: + nodes_qubit_mapping_dict: a dict about keys are the qubits used by the nodes and values are the new qubits + ''' + nodes_list_qubits_used = set() + for node in nodes_list: + if hasattr(node, 'pos') and node.pos is not None: + nodes_list_qubits_used = nodes_list_qubits_used | set(node.pos) + + mapping_pos = list(range(len(nodes_list_qubits_used))) + # mapping, get a dict + nodes_qubit_mapping_dict = dict(zip(sorted(list(nodes_list_qubits_used)), mapping_pos)) + nodes_qubit_mapping_dict + + return nodes_qubit_mapping_dict + +def nodelist_qubit_mapping_dict_reverse(nodes_list): + ''' + Args: + nodes_list: a list of nodes + Returns: + nodes_qubit_mapping_dict_reverse: a dict about keys are the new qubits and values are the qubits used by the nodes + ''' + nodes_qubit_mapping_dict = nodelist_qubit_mapping_dict(nodes_list) + # reverse mapping, get a dict + nodes_qubit_mapping_dict_reverse = {value: key for key, value in nodes_qubit_mapping_dict.items()} + + return nodes_qubit_mapping_dict_reverse + +# a function to map nodes_list +def nodes_list_mapping(nodes_list, nodes_qubit_mapping_dict): + ''' + Args: + nodes_list: the nodes list of instruction nodes + nodes_qubit_mapping_dict: the dict of the mapping qubits + + return: + nodes_list_mapping: the nodes_list after mapping qubits + ''' + nodes_qubit_mapping_dict + nodes_list_mapping = [] + for node in nodes_list: + node_new = copy.deepcopy(node) + if hasattr(node, 'pos') and node.pos is not None: + if isinstance(node.pos, list): + node_new.pos = [nodes_qubit_mapping_dict[qubit] for qubit in node.pos] + elif isinstance(node.pos, dict): + node_new.pos = {} + # the values of the dict are void, so we need to copy the values from the original dict + for qubit in node.pos: + node_new.pos[nodes_qubit_mapping_dict[qubit]] = copy.deepcopy(node.pos[qubit]) + nodes_list_mapping.append(node_new) + return nodes_list_mapping \ No newline at end of file diff --git a/src/quafu/dagcircuits/dag_circuit.py b/src/quafu/dagcircuits/dag_circuit.py new file mode 100644 index 0000000..58941ce --- /dev/null +++ b/src/quafu/dagcircuits/dag_circuit.py @@ -0,0 +1,286 @@ +import networkx as nx +from typing import Dict, Any, List + +from instruction_node import InstructionNode + +from networkx.classes.multidigraph import MultiDiGraph + +class DAGCircuit(MultiDiGraph): + def __init__(self,qubits_used=None, cbits_used=None, incoming_graph_data=None, **attr): + super().__init__(incoming_graph_data, **attr) + + if qubits_used is None: + self.qubits_used = set() + elif isinstance(qubits_used, set): + self.qubits_used = qubits_used + else: + raise ValueError('qubits_used should be a set or None') + + if cbits_used is None: + self.cbits_used = set() + elif isinstance(cbits_used, set): + self.cbits_used = cbits_used + else: + raise ValueError('cbits_used should be a set or None') + + # num of instruction nodes + self.num_instruction_nodes = 0 + + + # add new methods or override existing methods here. + + def update_qubits_used(self): + ''' + qubits_used is a set of qubits used in DAGCircuit + based on node -1's edges' labels, the qubits is the integer part of the label + + return: + qubits_used: set of qubits used in DAGCircuit + ''' + if -1 not in self.nodes: + raise ValueError('-1 should be in DAGCircuit, please add it first') + self.qubits_used = set([int(edge[2]['label'][1:]) for edge in self.out_edges(-1, data=True)]) + return self.qubits_used + + def update_cbits_used(self): + ''' + cbits_used is a set of cbits used in DAGCircuit + calculated by measurenode's cbits + return: + cbits_used: set of cbits used in DAGCircuit + ''' + for node in self.nodes: + # if node.has a attribute 'name' and node.name == 'measure' + if hasattr(node, 'name') and node.name == 'measure': + self.cbits_used = set(node.pos.values()) + return self.cbits_used + + def update_num_instruction_nodes(self): + ''' + num_instruction_nodes is the number of instruction nodes in DAGCircuit + ''' + if -1 not in self.nodes: + raise ValueError('-1 should be in DAGCircuit, please add it first') + if float('inf') not in self.nodes: + raise ValueError('float("inf") should be in DAGCircuit, please add it first') + self.num_instruction_nodes = len(self.nodes) - 2 + + for node in self.nodes: + if hasattr(node, 'name') and node.name == 'measure': + self.num_instruction_nodes -= 1 + return self.num_instruction_nodes + + + + def nodes_dict(self): + ''' + nodes_dict is a dictionary of nodes with the node label as key and the node as value. + without -1 and float('inf') + ''' + nodes_dict = {} + for node in nx.topological_sort(self): + if node != -1 and node != float('inf'): + nodes_dict[node.label] = node + return nodes_dict + + + def nodes_list(self): + ''' + nodes_list is a list of nodes without -1 and float('inf') + ''' + nodes_list = [] + for node in nx.topological_sort(self): + if node != -1 and node != float('inf'): + nodes_list.append(node) + return nodes_list + + def node_qubits_predecessors(self, node:InstructionNode): + ''' + node_qubits_predecessors is a dict of {qubits -> predecessors }of node + Args: + node in DAGCircuit, node should not be -1 + Returns: + node_qubits_predecessors: dict of {qubits -> predecessors }of node + ''' + # for edge in self.in_edges(node, data=True): + # print(edge[0], edge[1], edge[2]) + + if node not in self.nodes: + raise ValueError('node should be in DAGCircuit') + if node in [-1]: + raise ValueError('-1 has no predecessors') + + predecessor_nodes = [edge[0] for edge in self.in_edges(node, data=True)] + qubits_labels = [int(edge[2]['label'][1:]) for edge in self.in_edges(node, data=True)] + node_qubits_predecessors = dict(zip(qubits_labels, predecessor_nodes)) + return node_qubits_predecessors + + def node_qubits_successors(self, node:InstructionNode): + ''' + node_qubits_successors is a dict of {qubits -> successors }of node + Args: + node in DAGCircuit, node should not be float('inf') + Returns: + node_qubits_successors: dict of {qubits -> successors }of node + + + ''' + if node not in self.nodes: + raise ValueError('node should be in DAGCircuit') + if node in [float('inf')]: + raise ValueError('float("inf") has no successors') + successor_nodes = [edge[1] for edge in self.out_edges(node, data=True)] + qubits_labels = [int(edge[2]['label'][1:]) for edge in self.out_edges(node, data=True)] + node_qubits_successors = dict(zip(qubits_labels, successor_nodes)) + return node_qubits_successors + + def node_qubits_inedges(self, node:InstructionNode): + ''' + node_qubits_inedges is a dict of {qubits -> inedges }of node + Args: + node in DAGCircuit, node should not be -1 + Returns: + node_qubits_inedges: dict of {qubits -> inedges }of node + ''' + if node not in self.nodes: + raise ValueError('node should be in DAGCircuit') + if node in [-1]: + raise ValueError('-1 has no predecessors') + + inedges = [edge for edge in self.in_edges(node)] + qubits_labels = [int(edge[2]['label'][1:]) for edge in self.in_edges(node, data=True)] + node_qubits_inedges = dict(zip(qubits_labels, inedges)) + return node_qubits_inedges + + def node_qubits_outedges(self, node:InstructionNode): + ''' + node_qubits_outedges is a dict of {qubits -> outedges }of node + Args: + node in DAGCircuit, node should not be float('inf') + Returns: + node_qubits_outedges: dict of {qubits -> outedges }of node + ''' + if node not in self.nodes: + raise ValueError('node should be in DAGCircuit') + if node in [float('inf')]: + raise ValueError('float("inf") has no successors') + outedges = [edge for edge in self.out_edges(node)] + qubits_labels = [int(edge[2]['label'][1:]) for edge in self.out_edges(node, data=True)] + node_qubits_outedges = dict(zip(qubits_labels, outedges)) + return node_qubits_outedges + + def remove_instruction_node(self, gate:InstructionNode): + ''' + remove a gate from DAGCircuit, and all edges connected to it. + add new edges about qubits of removed gate between all predecessors and successors of removed gate. + Args: + gate: InstructionNode, gate should be in DAGCircuit, gate should not be -1 or float('inf') + ''' + + if gate not in self.nodes: + raise ValueError('gate should be in DAGCircuit') + if gate in [-1, float('inf')]: + raise ValueError('gate should not be -1 or float("inf")') + + qubits_predecessors = self.node_qubits_predecessors(gate) + qubits_successors = self.node_qubits_successors(gate) + for qubit in gate.pos: + if qubits_predecessors[qubit] != -1 and qubits_successors[qubit] != float('inf'): + self.add_edge(qubits_predecessors[qubit], qubits_successors[qubit], label=f'q{qubit}') + elif qubits_predecessors[qubit] == -1 and qubits_successors[qubit] != float('inf'): + self.add_edge(qubits_predecessors[qubit], qubits_successors[qubit], label=f'q{qubit}',color='green') + else: + self.add_edge(qubits_predecessors[qubit], qubits_successors[qubit], label=f'q{qubit}',color='red') + + self.remove_node(gate) + + # update qubits + self.qubits_used = self.update_qubits_used() + + + def merge_dag(self, other_dag): + ''' + merge other_dag into self + Args: + other_dag: DAGCircuit + Returns: + self: DAGCircuit + ''' + if not isinstance(other_dag, DAGCircuit): + raise ValueError('other_dag should be a DAGCircuit') + if other_dag == None: + return self + if self == None: + return other_dag + + # for the same qubits (intersection), + # remove the outgoing edges from the final node of the original DAG and the incoming edges from the initial node of the other DAG, + # then connect the corresponding tail and head nodes by adding edges + other_dag_qubits_used = other_dag.update_qubits_used() + self_qubits_used = self.update_qubits_used() + + insect_qubits = self_qubits_used & other_dag_qubits_used + end_edges_labels_1 = self.node_qubits_inedges(float('inf')) + start_edges_labels_2 = other_dag.node_qubits_outedges(-1) + + if len(insect_qubits) != 0: + for insect_qubit in insect_qubits: + self.remove_edges_from([end_edges_labels_1[insect_qubit]]) + other_dag.remove_edges_from([start_edges_labels_2[insect_qubit]]) + self.add_edge(end_edges_labels_1[insect_qubit][0], start_edges_labels_2[insect_qubit][1], label=f'q{insect_qubit}') + + # add other_dag's nodes and edges into self + # ！if we add edges, we don't need to add nodes again + self.add_edges_from(other_dag.edges(data=True)) + + # remove the edges between -1 and float('inf') + self.remove_edges_from([edge for edge in self.edges(data=True) if edge[0] == -1 and edge[1] == float('inf')]) + + # update qubits + self.qubits_used = self.update_qubits_used() + + def add_instruction_node(self, gate:InstructionNode,predecessors_dict:Dict[int,InstructionNode],successors_dict:Dict[int,InstructionNode]): + ''' + add a gate into DAGCircuit, and all edges connected to it. + add new edges about qubits of new gate between all predecessors and successors of new gate. + Args: + gate: InstructionNode, gate should not be -1 or float('inf') + predecessors_dict: dict of {qubits -> predecessors }of gate + successors_dict: dict of {qubits -> successors }of gate + ''' + if gate in [-1, float('inf')]: + raise ValueError('gate should not be -1 or float("inf")') + + #remove the edges between the predessors,successors about the qubits used by the added node + qubits_pre_out_edges = [] + qubits_suc_in_edges = [] + for qubit in gate.pos: + pre_out_edges = self.node_qubits_outedges(predecessors_dict[qubit]) + qubits_pre_out_edges.append(pre_out_edges[qubit]) + + suc_in_edges = self.node_qubits_inedges(successors_dict[qubit]) + qubits_suc_in_edges.append(suc_in_edges[qubit]) + + self.remove_edges_from(qubits_pre_out_edges) + self.remove_edges_from(qubits_suc_in_edges) + + # add the new node and edges + for qubit in gate.pos: + if predecessors_dict[qubit] == -1: + self.add_edge(predecessors_dict[qubit], gate, label=f'q{qubit}',color='green') + else: + self.add_edge(predecessors_dict[qubit], gate, label=f'q{qubit}') + if successors_dict[qubit] == float('inf'): + self.add_edge(gate, successors_dict[qubit], label=f'q{qubit}',color='red') + else: + self.add_edge(gate, successors_dict[qubit], label=f'q{qubit}') + + # update qubits + self.qubits_used = self.update_qubits_used() + + + def is_dag(self): + ''' + is_dag is a bool value to check if DAGCircuit is a DAG + ''' + return nx.is_directed_acyclic_graph(self) \ No newline at end of file diff --git a/src/quafu/dagcircuits/instruction_node.py b/src/quafu/dagcircuits/instruction_node.py new file mode 100644 index 0000000..3c55778 --- /dev/null +++ b/src/quafu/dagcircuits/instruction_node.py @@ -0,0 +1,40 @@ +from typing import Dict, Any, List, Union +import dataclasses + +@dataclasses.dataclass +class InstructionNode: + name:Any # gate.name + pos:Union[List[Any], Dict[Any,Any]] # gate.pos | Dict[Any,Any] for measure + paras:List[Any] # gate.paras + # matrix:List[Any] # for gate in [QuantumGate] + duration: Union[float,int] # for gate in [Delay,XYResonance,QuantumPulse] in quafu + unit:str # for gate in [Delay,XYResonance] in quafu + channel:str # for gate in [QuantumPulse] in quafu + time_func: Any # for gate in [QuantumPulse] in quafu + label:str # used for specifying the instruction node + + def __hash__(self): + return hash((type(self.name), tuple(self.pos) ,self.label)) + + def __str__(self): + if self.name == 'measure': + args = ','.join(str(q) for q in self.pos.keys()) + args += f'=>{",".join(str(c) for c in self.pos.values())}' + else: + args = ','.join(str(q) for q in self.pos) + + if self.paras == None: + return f'{self.label}{{{self.name}({args})}}' # no paras + else: + # if self.paras not a list, then make it a list of str of .3f float + if not isinstance(self.paras, list): + formatted_paras = [f'{self.paras:.3f}'] + else: + formatted_paras = [f'{p:.3f}' for p in self.paras] + + formatted_paras_str = ','.join(formatted_paras) + + return f'{self.label}{{{self.name}({args})}}({formatted_paras_str})' + + def __repr__(self): + return str(self) \ No newline at end of file diff --git a/src/quafu/elements/element_gates.py b/src/quafu/elements/element_gates.py deleted file mode 100644 index bef02d1..0000000 --- a/src/quafu/elements/element_gates.py +++ /dev/null @@ -1,247 +0,0 @@ -# This is the file for concrete element quantum gates -from .quantum_element import ControlledGate, FixedSingleQubitGate, ParaMultiQubitGate, ParaSingleQubitGate,\ - ControlledGate, FixedMultiQubitGate -import numpy as np -from typing import List -from scipy.linalg import sqrtm - - -def _u2matrix(_phi=0., _lambda=0.): - "OpenQASM 3.0 specification" - return np.array([[1., np.exp(-1.j * _lambda)], - [np.exp(1.j * _phi), np.exp((_phi + _lambda) * 1.j)]], dtype=complex) - - -def _u3matrix(_theta=0., _phi=0., _lambda=0.): - "OpenQASM 3.0 specification" - return np.array([[np.cos(0.5 * _theta), -np.exp(_lambda * 1.j) * np.sin(0.5 * _theta)], - [np.exp(_phi * 1.j) * np.sin(0.5 * _theta), - np.exp((_phi + _lambda) * 1.j) * np.cos(0.5 * _theta)]], dtype=complex) - - -def _rxmatrix(theta): - return np.array([[np.cos(0.5 * theta), -1.j * np.sin(0.5 * theta)], - [-1.j * np.sin(0.5 * theta), np.cos(0.5 * theta)]], dtype=complex) - - -def _rymatrix(theta): - return np.array([[np.cos(0.5 * theta), - np.sin(0.5 * theta)], - [np.sin(0.5 * theta), np.cos(0.5 * theta)]], dtype=complex) - - -def _rzmatrix(theta): - return np.array([[np.exp(-0.5j * theta), 0.], - [0., np.exp(0.5j * theta)]], dtype=complex) - -def _pmatrix(labda): - return np.array([[1, 0], - [0, np.exp(1j*labda)]] ,dtype=complex) - -def _rxxmatrix(theta): - """Unitary evolution of XX interaction""" - return np.array([[np.cos(theta/2), 0, 0, -1j*np.sin(theta/2)], - [0, np.cos(theta/2), -1j*np.sin(theta/2), 0], - [0, -1j*np.sin(theta/2), np.cos(theta/2), 0], - [-1j*np.sin(theta/2), 0, 0, np.cos(theta/2)] - ]) - -def _ryymatrix(theta): - """ Unitary evolution of YY interaction""" - return np.array([[np.cos(theta/2), 0, 0, 1j*np.sin(theta/2)], - [0, np.cos(theta/2), -1j*np.sin(theta/2), 0], - [0, -1j*np.sin(theta/2), np.cos(theta/2), 0], - [1j*np.sin(theta/2), 0, 0, np.cos(theta/2)] - ]) - -def _rzzmatrix(theta): - return np.array([[np.exp(-1j*theta/2), 0, 0, 0], - [0, np.exp(1j*theta/2), 0, 0], - [0, 0, np.exp(1j*theta/2), 0], - [0, 0, 0, np.exp(-1j*theta/2)] - ]) - -class IdGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("Id", pos, matrix = np.array([[1., 0.], - [0., 1.]], dtype=complex)) - -class HGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("H", pos, matrix=1 / np.sqrt(2) * np.array([[1., 1.], - [1., -1.]], dtype=complex)) - - -class XGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("X", pos, matrix=np.array([[0., 1.], - [1., 0.]], dtype=complex)) - - -class YGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("Y", pos, matrix=np.array([[0., -1.j], - [1.j, 0.]], dtype=complex)) - - -class ZGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("Z", pos, matrix=np.array([[1., 0.], - [0., -1.]], dtype=complex)) - -class SGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("S", pos, matrix=np.array([[1., 0.], - [0., 1.j]], dtype=complex)) - -class SdgGate(FixedSingleQubitGate): - def __init__(sell, pos: int): - super().__init__("Sdg", pos, matrix = np.array([[1., 0.], - [0., -1.j]], dtype=complex)) - -class TGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("T", pos, matrix=np.array([[1., 0.], - [0., np.exp(1.j*np.pi/4)]], dtype=complex)) - -class TdgGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("Tdg", pos, matrix=np.array([[1., 0.], - [0, np.exp(-1.j*np.pi/4)]] ,dtype=complex)) - -class WGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("W", pos, matrix=np.zeros((2, 2), dtype=complex)) - self.matrix = (XGate(0).matrix + YGate(0).matrix)/np.sqrt(2) - -class SXGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("SX", pos, matrix=np.zeros((2, 2), dtype=complex)) - self.matrix = sqrtm(XGate(0).matrix) - self.symbol = "√X" - -class SYGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("SY", pos, matrix=np.zeros((2, 2), dtype=complex)) - self.matrix = sqrtm(YGate(0).matrix) - self.symbol = "√Y" - -class SWGate(FixedSingleQubitGate): - def __init__(self, pos: int): - super().__init__("SW", pos, matrix=np.zeros((2, 2), dtype=complex)) - # sqrtm(WGate(0).matrix) - self.matrix = np.array([[0.5+0.5j, -np.sqrt(0.5)*1j], - [np.sqrt(0.5), 0.5+ 0.5j]], dtype=complex) - self.symbol = "√W" - -class RXGate(ParaSingleQubitGate): - def __init__(self, pos: int, paras): - super().__init__("RX", pos, paras, matrix=_rxmatrix) - - -class RYGate(ParaSingleQubitGate): - def __init__(self, pos: int, paras): - super().__init__("RY", pos, paras, matrix=_rymatrix) - - -class RZGate(ParaSingleQubitGate): - def __init__(self, pos: int, paras): - super().__init__("RZ", pos, paras, matrix=_rzmatrix) - - -class PhaseGate(ParaSingleQubitGate): - def __init__(self, pos: int, paras): - super().__init__("P", pos, paras, matrix=_pmatrix) - - -class iSwapGate(FixedMultiQubitGate): - def __init__(self, q1:int, q2:int): - super().__init__("iSWAP", [q1, q2], matrix=np.array([[1., 0., 0., 0.], - [0., 0., 1.j, 0.], - [0., 1.j, 0., 0.], - [0., 0., 0., 1.]], dtype=complex)) - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - - -class SwapGate(FixedMultiQubitGate): - def __init__(self, q1:int, q2:int): - super().__init__("SWAP", [q1, q2], matrix=np.array([[1., 0., 0., 0.], - [0., 0., 1., 0.], - [0., 1., 0., 0.], - [0., 0., 0., 1.]], dtype=complex)) - self.symbol = "x" - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - -class CXGate(ControlledGate): - def __init__(self, ctrl:int, targ:int): - super().__init__("CX", "X", [ctrl], [targ], None, matrix=XGate(0).matrix) - self.symbol = "+" - - -class CYGate(ControlledGate): - def __init__(self, ctrl:int, targ:int): - super().__init__("CY", "Y", [ctrl], [targ], None, matrix=YGate(0).matrix) - - -class CZGate(ControlledGate): - def __init__(self, ctrl:int, targ:int): - super().__init__("CZ", "Z", [ctrl], [targ], None, matrix=ZGate(0).matrix) - -class CSGate(ControlledGate): - def __init__(self, ctrl:int, targ:int): - super().__init__("CS", "S", [ctrl], [targ], None, matrix=SGate(0).matrix) - - -class CTGate(ControlledGate): - def __init__(self, ctrl:int, targ:int): - super().__init__("CT", "T", [ctrl], [targ], None, matrix=TGate(0).matrix) - -class CPGate(ControlledGate): - def __init__(self, ctrl:int, targ:int, paras): - super().__init__("CP", "P", [ctrl], [targ], paras, matrix=PhaseGate(0, paras).matrix) - -class ToffoliGate(ControlledGate): - def __init__(self, ctrl1:int, ctrl2:int, targ:int): - super().__init__("CCX", "X", [ctrl1, ctrl2], [targ], None, matrix=XGate(0).matrix) - -class FredkinGate(ControlledGate): - def __init__(self, ctrl:int, targ1:int, targ2:int): - super().__init__("CSWAP", "SWAP", [ctrl], [targ1, targ2], None, matrix=SwapGate(0, 1).matrix) - - -class RXXGate(ParaMultiQubitGate): - def __init__(self, q1:int, q2:int, theta): - super().__init__("RXX", [q1, q2], theta, matrix=_rxxmatrix(theta)) - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - -class RYYGate(ParaMultiQubitGate): - def __init__(self, q1:int, q2:int, theta): - super().__init__("RYY", [q1, q2], theta, matrix=_ryymatrix(theta)) - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - -class RZZGate(ParaMultiQubitGate): - def __init__(self, q1:int, q2:int, theta): - super().__init__("RZZ", [q1, q2], theta, matrix=_rzzmatrix(theta)) - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - -class MCXGate(ControlledGate): - def __init__(self, ctrls, targ:int): - super().__init__("MCX", "X", ctrls, [targ], None, matrix=XGate(0).matrix) - -class MCYGate(ControlledGate): - def __init__(self, ctrls, targ:int): - super().__init__("MCY", "Y", ctrls, [targ], None, matrix=YGate(0).matrix) - -class MCZGate(ControlledGate): - def __init__(self, ctrls, targ:int): - super().__init__("MCZ", "Z", ctrls, [targ], None, matrix=ZGate(0).matrix) - diff --git a/src/quafu/elements/element_gates/__init__.py b/src/quafu/elements/element_gates/__init__.py new file mode 100644 index 0000000..434d0e6 --- /dev/null +++ b/src/quafu/elements/element_gates/__init__.py @@ -0,0 +1,20 @@ +from .pauli import XGate, YGate, ZGate, IdGate, WGate, SWGate, SXGate, SXdgGate, SYGate, SYdgGate +from .clifford import HGate, SGate, SdgGate, TGate, TdgGate +from .phase import PhaseGate +from .rotation import RXGate, RYGate, RZGate, RXXGate, RYYGate, RZZGate +from .swap import SwapGate, ISwapGate +from .c11 import CXGate, CYGate, CZGate, CSGate, CTGate, CPGate +from .c21 import ToffoliGate +from .c12 import FredkinGate +from .cm1 import MCXGate, MCYGate, MCZGate, ControlledU +from .unitary import UnitaryDecomposer + +__all__ = ['XGate', 'YGate', 'ZGate', 'IdGate', 'WGate', 'SWGate', + 'PhaseGate', + 'RXGate', 'RYGate', 'RZGate', 'RXXGate', 'RYYGate', 'RZZGate', + 'SwapGate', 'ISwapGate', + 'CXGate', 'CYGate', 'CZGate', 'CSGate', 'CTGate', 'CPGate', + 'ToffoliGate', + 'FredkinGate', + 'MCXGate', 'MCYGate', 'MCZGate', 'ControlledU', + 'UnitaryDecomposer'] diff --git a/src/quafu/elements/element_gates/c11.py b/src/quafu/elements/element_gates/c11.py new file mode 100644 index 0000000..78c9b62 --- /dev/null +++ b/src/quafu/elements/element_gates/c11.py @@ -0,0 +1,64 @@ +from ..quantum_element import ControlledGate +from abc import ABC +from .matrices import XMatrix, YMatrix, ZMatrix, SMatrix, TMatrix, pmatrix + + +class _C11Gate(ControlledGate, ABC): + ct_dims = (1, 1, 2) + + +class CXGate(_C11Gate): + name = "CX" + + def __init__(self, ctrl: int, targ: int): + super().__init__("X", [ctrl], [targ], None, tar_matrix=XMatrix) + self.symbol = "+" + + +class CYGate(_C11Gate): + name = "CY" + + def __init__(self, ctrl: int, targ: int): + super().__init__("Y", [ctrl], [targ], None, tar_matrix=YMatrix) + + +class CZGate(_C11Gate): + name = "CZ" + + def __init__(self, ctrl: int, targ: int): + super().__init__("Z", [ctrl], [targ], None, tar_matrix=ZMatrix) + + +class CSGate(_C11Gate): + name = "CS" + + def __init__(self, ctrl: int, targ: int): + super().__init__("S", [ctrl], [targ], None, tar_matrix=SMatrix) + + def to_qasm(self): + return "cp(pi/2) " + "q[%d],q[%d]" % (self.pos[0], self.pos[1]) + + +class CTGate(_C11Gate): + name = "CT" + + def __init__(self, ctrl: int, targ: int): + super().__init__("T", [ctrl], [targ], None, tar_matrix=TMatrix) + + def to_qasm(self): + return "cp(pi/4) " + "q[%d],q[%d]" % (self.pos[0], self.pos[1]) + + +class CPGate(_C11Gate): + name = "CP" + + def __init__(self, ctrl: int, targ: int, paras): + super().__init__("P", [ctrl], [targ], paras, tar_matrix=pmatrix(paras)) + + +ControlledGate.register_gate(CXGate) +ControlledGate.register_gate(CYGate) +ControlledGate.register_gate(CZGate) +ControlledGate.register_gate(CSGate) +ControlledGate.register_gate(CTGate) +ControlledGate.register_gate(CPGate) diff --git a/src/quafu/elements/element_gates/c12.py b/src/quafu/elements/element_gates/c12.py new file mode 100644 index 0000000..c51e9e6 --- /dev/null +++ b/src/quafu/elements/element_gates/c12.py @@ -0,0 +1,12 @@ +from ..quantum_element import ControlledGate +from .matrices import SwapMatrix + + +class FredkinGate(ControlledGate): + name = "CSWAP" + + def __init__(self, ctrl: int, targ1: int, targ2: int): + super().__init__("SWAP", [ctrl], [targ1, targ2], None, tar_matrix=SwapMatrix) + + +ControlledGate.register_gate(FredkinGate) diff --git a/src/quafu/elements/element_gates/c21.py b/src/quafu/elements/element_gates/c21.py new file mode 100644 index 0000000..5419cbd --- /dev/null +++ b/src/quafu/elements/element_gates/c21.py @@ -0,0 +1,13 @@ +from ..quantum_element import ControlledGate +from .matrices import XMatrix + + +class ToffoliGate(ControlledGate): + name = "CCX" + + def __init__(self, ctrl1: int, ctrl2: int, targ: int): + super().__init__("X", [ctrl1, ctrl2], [targ], None, tar_matrix=XMatrix) + + +ControlledGate.register_gate(ToffoliGate) + diff --git a/src/quafu/elements/element_gates/clifford.py b/src/quafu/elements/element_gates/clifford.py new file mode 100644 index 0000000..ea10936 --- /dev/null +++ b/src/quafu/elements/element_gates/clifford.py @@ -0,0 +1,53 @@ +import numpy as np + +from quafu.elements.element_gates.matrices import HMatrix +from quafu.elements.quantum_element import FixedSingleQubitGate + + +class HGate(FixedSingleQubitGate): + name = "H" + matrix = HMatrix + + def __init__(self, pos: int): + super().__init__(pos) + + +class SGate(FixedSingleQubitGate): + name = "S" + matrix = np.array([[1., 0.], + [0., 1.j]], dtype=complex) + + def __init__(self, pos: int): + super().__init__(pos) + + +class SdgGate(FixedSingleQubitGate): + name = "Sdg" + matrix = SGate.matrix.conj().T + + def __init__(self, pos: int): + super().__init__(pos) + + +class TGate(FixedSingleQubitGate): + name = "T" + matrix = np.array([[1., 0.], + [0., np.exp(1.j * np.pi / 4)]], dtype=complex) + + def __init__(self, pos: int): + super().__init__(pos) + + +class TdgGate(FixedSingleQubitGate): + name = "Tdg" + matrix = TGate.matrix.conj().T + + def __init__(self, pos: int): + super().__init__(pos) + + +FixedSingleQubitGate.register_gate(HGate) +FixedSingleQubitGate.register_gate(SGate) +FixedSingleQubitGate.register_gate(SdgGate) +FixedSingleQubitGate.register_gate(TGate) +FixedSingleQubitGate.register_gate(TdgGate) diff --git a/src/quafu/elements/element_gates/cm1.py b/src/quafu/elements/element_gates/cm1.py new file mode 100644 index 0000000..231d24b --- /dev/null +++ b/src/quafu/elements/element_gates/cm1.py @@ -0,0 +1,47 @@ +from typing import List, Union + +from .matrices import XMatrix, YMatrix, ZMatrix +from ..quantum_element import ControlledGate, SingleQubitGate, MultiQubitGate + + +class MCXGate(ControlledGate): + name = "MCX" + + def __init__(self, ctrls, targ: int): + super().__init__("X", ctrls, [targ], None, tar_matrix=XMatrix) + + +class MCYGate(ControlledGate): + name = "MCY" + + def __init__(self, ctrls, targ: int): + super().__init__("Y", ctrls, [targ], None, tar_matrix=YMatrix) + + +class MCZGate(ControlledGate): + name = "MCZ" + + def __init__(self, ctrls, targ: int): + super().__init__("Z", ctrls, [targ], None, tar_matrix=ZMatrix) + + +class ControlledU(ControlledGate): + """ Controlled gate class, where the matrix act non-trivially on target qubits""" + name = 'CU' + + def __init__(self, ctrls: List[int], u: Union[SingleQubitGate, MultiQubitGate]): + self.targ_gate = u + targs = u.pos + if isinstance(targs, int): + targs = [targs] + + super().__init__(u.name, ctrls, targs, u.paras, tar_matrix=self.targ_gate.get_targ_matrix()) + + def get_targ_matrix(self, reverse_order=False): + return self.targ_gate.get_targ_matrix(reverse_order) + + +ControlledGate.register_gate(MCXGate) +ControlledGate.register_gate(MCYGate) +ControlledGate.register_gate(MCZGate) +ControlledGate.register_gate(ControlledU) diff --git a/src/quafu/elements/element_gates/matrices/__init__.py b/src/quafu/elements/element_gates/matrices/__init__.py new file mode 100644 index 0000000..f69bd9f --- /dev/null +++ b/src/quafu/elements/element_gates/matrices/__init__.py @@ -0,0 +1,2 @@ +from .mat_lib import * +from .mat_utils import stack_matrices, is_zero, is_hermitian diff --git a/src/quafu/elements/element_gates/matrices/mat_lib.py b/src/quafu/elements/element_gates/matrices/mat_lib.py new file mode 100644 index 0000000..b40ebec --- /dev/null +++ b/src/quafu/elements/element_gates/matrices/mat_lib.py @@ -0,0 +1,128 @@ +import numpy as np + +IdMatrix = np.eye(2, dtype=complex) +XMatrix = np.array([[0., 1.], [1., 0.]], dtype=complex) +YMatrix = np.array([[0., -1.j], [1.j, 0.]], dtype=complex) +ZMatrix = np.array([[1., 0.], [0., -1.]], dtype=complex) +SMatrix = np.array([[1., 0.], [0., 1.j]], dtype=complex) +SXMatrix = np.array([[1., 1.j], [1.j, 1.]], dtype=complex) / np.sqrt(2) +SYMatrix = np.array([[1., -1.], [1., 1.]], dtype=complex) / np.sqrt(2) +TMatrix = np.array([[1., 0.], [0., np.exp(1.j * np.pi / 4)]], dtype=complex) +WMatrix = (XMatrix + YMatrix) / np.sqrt(2) +SWMatrix = np.array([[0.5 + 0.5j, -np.sqrt(0.5) * 1j], + [np.sqrt(0.5), 0.5 + 0.5j]], dtype=complex) +HMatrix = (XMatrix + ZMatrix) / np.sqrt(2) +SwapMatrix = np.array([[1., 0., 0., 0.], + [0., 0., 1., 0.], + [0., 1., 0., 0.], + [0., 0., 0., 1.]], dtype=complex) +ISwapMatrix = np.array([[1., 0., 0., 0.], + [0., 0., 1.j, 0.], + [0., 1.j, 0., 0.], + [0., 0., 0., 1.]], dtype=complex) +CXMatrix = np.array([[1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 0., 1.], + [0., 0., 1., 0.]], dtype=complex) +CYMatrix = np.array([[1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 0., -1.j], + [0., 0., 1.j, 0.]], dtype=complex) +CZMatrix = np.array([[1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 1., 0.], + [0., 0., 0., -1.]], dtype=complex) +ToffoliMatrix = np.array([[1., 0., 0., 0., 0., 0., 0., 0.], + [0., 1., 0., 0., 0., 0., 0., 0.], + [0., 0., 1., 0., 0., 0., 0., 0.], + [0., 0., 0., 1., 0., 0., 0., 0.], + [0., 0., 0., 0., 1., 0., 0., 0.], + [0., 0., 0., 0., 0., 1., 0., 0.], + [0., 0., 0., 0., 0., 0., 0., 1.], + [0., 0., 0., 0., 0., 0., 1., 0.]], dtype=complex) +FredkinMatrix = np.array([[1., 0., 0., 0., 0., 0., 0., 0.], + [0., 1., 0., 0., 0., 0., 0., 0.], + [0., 0., 1., 0., 0., 0., 0., 0.], + [0., 0., 0., 1., 0., 0., 0., 0.], + [0., 0., 0., 0., 1., 0., 0., 0.], + [0., 0., 0., 0., 0., 0., 1., 0.], + [0., 0., 0., 0., 0., 1., 0., 0.], + [0., 0., 0., 0., 0., 0., 0., 1.]], dtype=complex) + + +def u2matrix(_phi=0., _lambda=0.): + """OpenQASM 3.0 specification""" + return np.array([[1., np.exp(-1.j * _lambda)], + [np.exp(1.j * _phi), np.exp((_phi + _lambda) * 1.j)]], dtype=complex) + + +def u3matrix(_theta=0., _phi=0., _lambda=0.): + """OpenQASM 3.0 specification""" + return np.array([[np.cos(0.5 * _theta), -np.exp(_lambda * 1.j) * np.sin(0.5 * _theta)], + [np.exp(_phi * 1.j) * np.sin(0.5 * _theta), + np.exp((_phi + _lambda) * 1.j) * np.cos(0.5 * _theta)]], dtype=complex) + + +def rx_mat(theta): + return np.array([[np.cos(0.5 * theta), -1.j * np.sin(0.5 * theta)], + [-1.j * np.sin(0.5 * theta), np.cos(0.5 * theta)]], dtype=complex) + + +def ry_mat(theta): + return np.array([[np.cos(0.5 * theta), - np.sin(0.5 * theta)], + [np.sin(0.5 * theta), np.cos(0.5 * theta)]], dtype=complex) + + +def rz_mat(theta): + return np.array([[np.exp(-0.5j * theta), 0.], + [0., np.exp(0.5j * theta)]], dtype=complex) + + +def pmatrix(labda): + return np.array([[1, 0], + [0, np.exp(1j * labda)]], dtype=complex) + + +def rxx_mat(theta): + """Unitary evolution of XX interaction""" + return np.array([[np.cos(theta / 2), 0, 0, -1j * np.sin(theta / 2)], + [0, np.cos(theta / 2), -1j * np.sin(theta / 2), 0], + [0, -1j * np.sin(theta / 2), np.cos(theta / 2), 0], + [-1j * np.sin(theta / 2), 0, 0, np.cos(theta / 2)] + ]) + + +def ryy_mat(theta): + """ Unitary evolution of YY interaction""" + c = np.cos(theta / 2) + s = 1j * np.sin(theta / 2) + return np.array([[c, 0, 0, +s], + [0, c, -s, 0], + [0, -s, c, 0], + [+s, 0, 0, c] + ]) + + +def rzz_mat(theta): + return np.array([[np.exp(-1j * theta / 2), 0, 0, 0], + [0, np.exp(1j * theta / 2), 0, 0], + [0, 0, np.exp(1j * theta / 2), 0], + [0, 0, 0, np.exp(-1j * theta / 2)] + ]) + +# def su2_matrix(gamma: float, beta: float, delta: float): +# """ +# SU = Rz(beta)Ry(gamma)Rz(delta). +# +# Symbol convention is the same as in the textbook +# of Chuang and Nielsen. +# """ +# s, c = np.sin(gamma / 2), np.cos(gamma / 2) +# alpha1, alpha2 = (delta + beta) / 2, (delta - beta) / 2 +# su2_mat = np.array([[np.exp(-1.j * alpha1) * c, -np.exp(-1.j * alpha2) * s], +# [np.exp(1.j * alpha2) * s, np.exp(1.j * alpha1) * c]]) +# return su2_mat +# +# +# def u2_matrix(alpha: float, gamma: float, beta: float, delta: float): +# return np.exp(1.j * alpha) * su2_matrix(gamma, beta, delta) diff --git a/src/quafu/elements/element_gates/matrices/mat_utils.py b/src/quafu/elements/element_gates/matrices/mat_utils.py new file mode 100644 index 0000000..479799c --- /dev/null +++ b/src/quafu/elements/element_gates/matrices/mat_utils.py @@ -0,0 +1,121 @@ +import cmath + +import numpy as np +from numpy import ndarray + +from .mat_lib import IdMatrix + + +def split_matrix(matrix: ndarray): + """ + Evenly split a matrix into 4 sub-matrices. + """ + top, bottom = np.vsplit(matrix, 2) + t_left, t_right = np.hsplit(top, 2) + b_left, b_right = np.hsplit(bottom, 2) + return t_left, t_right, b_left, b_right + + +def stack_matrices(t_left, t_right, b_left, b_right): + """ + Stack 4 sub-matrices into a matrix. + """ + top = np.hstack((t_left, t_right)) + bottom = np.hstack((b_left, b_right)) + mat = np.vstack((top, bottom)) + return mat + + +def multi_kron(op1, op2, ind1, ind2, nspin): + tmp = 1 + for i in range(nspin): + if i == ind1: + tmp = np.kron(tmp, op1) + elif i == ind2: + tmp = np.kron(tmp, op2) + else: + tmp = np.kron(tmp, IdMatrix) + return tmp + + +def general_kron(op, ind, nqubit): + tmp = 1 + for i in range(nqubit): + if i == ind: + tmp = np.kron(tmp, op) + else: + tmp = np.kron(tmp, IdMatrix) + return tmp + + +####################################################### +def is_zero(a): + return not np.any(np.absolute(a) > 1e-8) + + +def is_approx(a, b, thres=1e-6): + # TODO: seems there are some very small elements that cannot be compared correctly + # if not np.allclose(a, b, rtol=thres, atol=thres): + # print(np.sum(a-b)) + return np.allclose(a, b, rtol=thres, atol=thres) + + +def is_unitary(matrix): + mat_dg = np.conjugate(matrix).T + id_mat = np.eye(matrix.shape[0]) + return is_approx(mat_dg @ matrix, id_mat) and is_approx(matrix @ mat_dg, id_mat) + + +def is_hermitian(matrix): + tmp = np.conjuate(matrix).T + return is_approx(tmp, matrix) + + +def is_diagonal(matrix: ndarray): + diag = np.diag(matrix) + diag_mat = np.diag(diag) + return is_approx(matrix, diag_mat) + + +def is_kron_with_id2(matrix): + """ + Check if the matrix is a Kronecker product of a matrix and identity matrix. + """ + nsize = matrix.shape[0] + + a_cond = is_zero(matrix[0:nsize:2, 1:nsize:2]) + b_cond = is_zero(matrix[1:nsize:2, 0:nsize:2]) + c_cond = is_approx(matrix[0, :-1], matrix[1, 1:]) + d_cond = is_approx(matrix[-2, :-1], matrix[-1, 1:]) + + return a_cond and b_cond and c_cond and d_cond + + +####################################################### +def get_global_phase(unitary): + """ Get the global phase of arbitrary unitary, and get the special unitary. + + Args: + unitary (np.array): arbitrary unitary + Returns: + global_phase: the global phase of arbitrary unitary + special_unitary (np.array) + """ + coefficient = np.linalg.det(unitary) ** (-0.5) + global_phase = -cmath.phase(coefficient) + special_unitary = coefficient * unitary + return global_phase, special_unitary + + +def matrix_distance_squared(unitary1, unitary2): + """ Used to compare the distance of two matrices. The global phase is ignored. + + Args: + unitary1 (np.array): A unitary matrix in the form of a numpy ndarray. + unitary2 (np.array): Another unitary matrix. + + Returns: + Float : A single value between 0 and 1 indicating how closely unitary1 and unitary2 match. + A value close to 0 indicates that unitary1 and unitary2 are the same unitary. + """ + return np.abs(1 - np.abs(np.sum(np.multiply(unitary1, np.conj(unitary2)))) / unitary1.shape[0]) diff --git a/src/quafu/elements/element_gates/pauli.py b/src/quafu/elements/element_gates/pauli.py new file mode 100644 index 0000000..9f513e6 --- /dev/null +++ b/src/quafu/elements/element_gates/pauli.py @@ -0,0 +1,113 @@ +import numpy as np + +from .matrices import XMatrix, YMatrix, ZMatrix, WMatrix, SWMatrix +from ..quantum_element import FixedSingleQubitGate + + +class IdGate(FixedSingleQubitGate): + name = "Id" + matrix = XMatrix + + def __init__(self, pos: int): + super().__init__(pos) + + +class XGate(FixedSingleQubitGate): + name = "X" + matrix = XMatrix + + def __init__(self, pos: int): + super().__init__(pos) + + +class YGate(FixedSingleQubitGate): + name = "Y" + matrix = YMatrix + + def __init__(self, pos: int): + super().__init__(pos) + + +class ZGate(FixedSingleQubitGate): + name = "Z" + matrix = ZMatrix + + def __init__(self, pos: int): + super().__init__(pos) + + +class WGate(FixedSingleQubitGate): + name = "W" + matrix = WMatrix + + def __init__(self, pos: int): + super().__init__(pos) + self.symbol = "W" + + def to_qasm(self): + return "rz(-pi/4) q[%d];\nrx(pi) q[%d];\nrz(pi/4) q[%d]" % (self.pos, self.pos, self.pos) + + +class SWGate(FixedSingleQubitGate): + name = "SW" + matrix = SWMatrix + + def __init__(self, pos: int): + super().__init__(pos) + self.symbol = "√W" + + def to_qasm(self): + return "rz(-pi/4) q[%d];\nrx(pi/2) q[%d];\nrz(pi/4) q[%d]" % (self.pos, self.pos, self.pos) + + +class SXGate(FixedSingleQubitGate): + name = "SX" + matrix = np.array([[0.5 + 0.5j, 0.5 - 0.5j], + [0.5 - 0.5j, 0.5 + 0.5j]], dtype=complex) + + def __init__(self, pos: int): + super().__init__(pos) + + +class SXdgGate(FixedSingleQubitGate): + name = "SXdg" + matrix = SXGate.matrix.conj().T + + def __init__(self, pos: int): + super().__init__(pos) + self.symbol = "√X" + + +class SYGate(FixedSingleQubitGate): + name = "SY" + matrix = np.array([[0.5 + 0.5j, -0.5 - 0.5j], + [0.5 + 0.5j, 0.5 + 0.5j]], dtype=complex) + + def __init__(self, pos: int): + super().__init__(pos) + self.symbol = "√Y" + + +class SYdgGate(FixedSingleQubitGate): + name = "SYdg" + matrix = SYGate.matrix.conj().T + + def __init__(self, pos: int): + super().__init__(pos) + self.symbol = "√Y†" + + def to_qasm(self): + # TODO: this seems incorrect + return "ry(pi/2) q[%d]" % self.pos + + +FixedSingleQubitGate.register_gate(IdGate) +FixedSingleQubitGate.register_gate(XGate) +FixedSingleQubitGate.register_gate(YGate) +FixedSingleQubitGate.register_gate(ZGate) +FixedSingleQubitGate.register_gate(WGate) +FixedSingleQubitGate.register_gate(SWGate) +FixedSingleQubitGate.register_gate(SXGate) +FixedSingleQubitGate.register_gate(SXdgGate) +FixedSingleQubitGate.register_gate(SYGate) +FixedSingleQubitGate.register_gate(SYdgGate) diff --git a/src/quafu/elements/element_gates/phase.py b/src/quafu/elements/element_gates/phase.py new file mode 100644 index 0000000..f40d119 --- /dev/null +++ b/src/quafu/elements/element_gates/phase.py @@ -0,0 +1,16 @@ +from ..quantum_element import ParaSingleQubitGate +from .matrices import pmatrix + + +class PhaseGate(ParaSingleQubitGate): + name = "P" + + def __init__(self, pos: int, paras: float = 0.): + super().__init__(pos, paras=paras) + + @property + def matrix(self): + return pmatrix(self.paras) + + +ParaSingleQubitGate.register_gate(PhaseGate) diff --git a/src/quafu/elements/element_gates/rotation.py b/src/quafu/elements/element_gates/rotation.py new file mode 100644 index 0000000..69ee7f7 --- /dev/null +++ b/src/quafu/elements/element_gates/rotation.py @@ -0,0 +1,85 @@ +from ..quantum_element import ParaMultiQubitGate, ParaSingleQubitGate +from .matrices import rx_mat, ry_mat, rz_mat, rxx_mat, ryy_mat, rzz_mat + + +class RXGate(ParaSingleQubitGate): + name = "RX" + + def __init__(self, pos: int, paras: float = 0.): + super().__init__(pos, paras) + + @property + def matrix(self): + return rx_mat(self.paras) + + +class RYGate(ParaSingleQubitGate): + name = "RY" + + def __init__(self, pos: int, paras: float = 0.): + super().__init__(pos, paras) + + @property + def matrix(self): + return ry_mat(self.paras) + + +class RZGate(ParaSingleQubitGate): + name = "RZ" + + def __init__(self, pos: int, paras: float = 0.): + super().__init__(pos, paras) + + @property + def matrix(self): + return rz_mat(self.paras) + + +class RXXGate(ParaMultiQubitGate): + name = "RXX" + + def __init__(self, q1: int, q2: int, theta: float = 0.): + super().__init__([q1, q2], paras=theta) + + @property + def matrix(self): + return rxx_mat(self.paras) + + def get_targ_matrix(self, reverse_order=False): + return self.matrix + + +class RYYGate(ParaMultiQubitGate): + name = "RYY" + + def __init__(self, q1: int, q2: int, theta: float = 0.): + super().__init__([q1, q2], paras=theta) + + @property + def matrix(self): + return ryy_mat(self.paras) + + def get_targ_matrix(self, reverse_order=False): + return self.matrix + + +class RZZGate(ParaMultiQubitGate): + name = "RZZ" + + def __init__(self, q1: int, q2: int, theta: float = 0.): + super().__init__([q1, q2], paras=theta) + + @property + def matrix(self): + return rzz_mat(self.paras) + + def get_targ_matrix(self, reverse_order=False): + return self.matrix + + +ParaSingleQubitGate.register_gate(RXGate) +ParaSingleQubitGate.register_gate(RYGate) +ParaSingleQubitGate.register_gate(RZGate) +ParaMultiQubitGate.register_gate(RXXGate) +ParaMultiQubitGate.register_gate(RYYGate) +ParaMultiQubitGate.register_gate(RZZGate) diff --git a/src/quafu/elements/element_gates/swap.py b/src/quafu/elements/element_gates/swap.py new file mode 100644 index 0000000..1e7ebcd --- /dev/null +++ b/src/quafu/elements/element_gates/swap.py @@ -0,0 +1,29 @@ +from ..quantum_element import FixedMultiQubitGate +from .matrices import ISwapMatrix, SwapMatrix + + +class ISwapGate(FixedMultiQubitGate): + name = "iSWAP" + matrix = ISwapMatrix + + def __init__(self, q1: int, q2: int): + super().__init__([q1, q2]) + + def get_targ_matrix(self, reverse_order=False): + return self.matrix + + +class SwapGate(FixedMultiQubitGate): + name = "SWAP" + matrix = SwapMatrix + + def __init__(self, q1: int, q2: int): + super().__init__([q1, q2]) + self.symbol = "x" + + def get_targ_matrix(self, reverse_order=False): + return self.matrix + + +FixedMultiQubitGate.register_gate(ISwapGate) +FixedMultiQubitGate.register_gate(SwapGate) diff --git a/src/quafu/elements/element_gates/unitary/__init__.py b/src/quafu/elements/element_gates/unitary/__init__.py new file mode 100644 index 0000000..22d04e9 --- /dev/null +++ b/src/quafu/elements/element_gates/unitary/__init__.py @@ -0,0 +1,2 @@ +from .decomposer import UnitaryDecomposer + diff --git a/src/quafu/elements/element_gates/unitary/decomposer.py b/src/quafu/elements/element_gates/unitary/decomposer.py new file mode 100644 index 0000000..ddbb42b --- /dev/null +++ b/src/quafu/elements/element_gates/unitary/decomposer.py @@ -0,0 +1,358 @@ +import cmath +import math + +import numpy as np +import scipy +from numpy import ndarray + +from ..matrices import rz_mat, ry_mat, CXMatrix +from ..matrices import mat_utils as mu + + +class UnitaryDecomposer(object): + def __init__(self, array: ndarray, qubits, verbose: bool = False): + self.array = array + self.qubit_num = len(qubits) + self._check_unitary(array, self.qubit_num) + + self.qubits = qubits + self.verbose = verbose + + self.Mk_table = genMk_table(self.qubit_num) # initialize the general M^k lookup table + self.gate_list = [] + + def __call__(self, *args, **kwargs): + _matrix = self.array + self._decompose_matrix(_matrix, self.qubits) + if self.verbose: + print("Decomposition done.") + + @staticmethod + def _check_unitary(array: ndarray, qubits): + assert len(array.shape) == 2 + assert mu.is_unitary(array) + assert len(array) == 2 ** qubits + + def _decompose_matrix(self, _matrix, qubits): + qubit_num = len(qubits) + self._check_unitary(_matrix, qubit_num) + if qubit_num == 1: + # self.gates_list.append((_matrix, qubits, 'U')) + # ZYZ decomposition for single-qubit gate + gamma, beta, alpha, global_phase = zyz_decomposition(_matrix) + self.gate_list.append((rz_mat(gamma), qubits, 'RZ', gamma)) + self.gate_list.append((ry_mat(beta), qubits, 'RY', beta)) + self.gate_list.append((rz_mat(alpha), qubits, 'RZ', alpha)) + return None + + # if num_qubit == 2: + # self.circuit.append((_matrix, qubits)) + # print("Two qubit gate do not need to decompose.") + # zyz_decomp(matrix) + U00, U01, U10, U11 = mu.split_matrix(_matrix) + + # if bottomLeftCorner(n)==0 and topRightCorner(n)==0: + if mu.is_zero(U01) and mu.is_zero(U10): + if mu.is_approx(U00, U11): + if self.verbose: + print('Optimization: Unitaries are equal, ' + 'skip one step in the recursion for unitaries of size') + self._decompose_matrix(U00, qubits[1:]) + else: + if self.verbose: + print("Optimization: q2 is zero, only demultiplexing will be performed.") + V, D, W = demultiplexing(U00, U11) + self._decompose_matrix(W, qubits[1:]) + self.multi_controlled_z(D, qubits[1:], qubits[0]) + self._decompose_matrix(V, qubits[1:]) + # Check if the kronecker product of a bigger matrix and the identity matrix. + # By checking if the first row is equal to the second row one over, and if the last two rows are equal + # Which means the last qubit is not affected by this gate + elif mu.is_kron_with_id2(_matrix): + # print("The last qubits not affect.") + nsize = len(_matrix) + self._decompose_matrix(_matrix[0:nsize:2, 0:nsize:2], qubits[:-1]) + else: + # print("CSD decomposition.") + L0, L1, R0, R1, c, s = fat_csd(_matrix) + V, D, W = demultiplexing(R0, R1) + + self._decompose_matrix(W, qubits[1:]) + self.multi_controlled_z(D, qubits[1:], qubits[0]) + self._decompose_matrix(V, qubits[1:]) + self.multi_controlled_y(s, qubits[1:], qubits[0]) + + V, D, W = demultiplexing(L0, L1) + self._decompose_matrix(W, qubits[1:]) + self.multi_controlled_z(D, qubits[1:], qubits[0]) + self._decompose_matrix(V, qubits[1:]) + + def multi_controlled_z(self, D, qubits, target_qubit): + print(D.shape[0]) + assert len(qubits) == int(math.log(D.shape[0], 2)) + num_qubit = len(qubits) + # print("The size of D matrix {}".format(len(qubits))) + # print(qubits) + + # alphas = -2*1j*np.log(np.diag(D)) + alphas = 2 * 1j * np.log(np.diag(D)) + Mk = self.Mk_table[num_qubit - 1] + thetas = np.linalg.solve(Mk, alphas) + # print(thetas) + assert len(thetas) == 2 ** num_qubit + + index = get_multi_control_index(num_qubit) + assert len(index) == len(thetas) + + for i in range(len(index)): + control_qubit = qubits[index[i]] + self.gate_list.append((rz_mat(thetas[i]), [target_qubit], 'RZ', thetas[i])) + self.gate_list.append((CXMatrix, [control_qubit, target_qubit], 'CX')) + + def multi_controlled_y(self, ss, qubits, target_qubit): + assert len(qubits) == int(math.log(ss.shape[0], 2)) + num_qubit = len(qubits) + + alphas = -2 * np.arcsin(np.diag(ss)) + Mk = self.Mk_table[num_qubit - 1] + thetas = np.linalg.solve(Mk, alphas) + # print(thetas) + assert len(thetas) == 2 ** num_qubit + + index = get_multi_control_index(num_qubit) + assert len(index) == len(thetas) + + for i in range(len(index)): + control_qubit = qubits[index[i]] + self.gate_list.append((ry_mat(thetas[i]), [target_qubit], 'RY', thetas[i])) + self.gate_list.append((CXMatrix, [control_qubit, target_qubit], 'CX')) + + def apply_to_qc(self, qc): + if len(self.gate_list) == 0: + self() + + for g in self.gate_list: + if g[2] == 'CX': + qc.cnot(g[1][0], g[1][1]) + elif g[2] == 'RY': + qc.ry(g[1][0], g[3].real) + elif g[2] == 'RZ': + qc.rz(g[1][0], g[3].real) + elif g[2] == 'U': + gamma, beta, alpha, global_phase = zyz_decomposition(g[0]) + qc.rz(g[1][0], gamma) + qc.ry(g[1][0], beta) + qc.rz(g[1][0], alpha) + else: + raise Exception("Unknown gate type or incorrect str: {}".format(g[2])) + + return qc + + +def zyz_decomposition(unitary): + """ ZYZ decomposition of arbitrary single-qubit gate (unitary). + SU = Rz(gamma) * Ry(beta) * Rz(alpha) + + Args: + unitary (np.array): arbitrary unitary + Returns: + global_phase: the global phase of arbitrary unitary + special_unitary (np.array) + """ + if unitary.shape[0] == 2: + global_phase, special_unitary = mu.get_global_phase(unitary) + beta = 2 * math.atan2(abs(special_unitary[1, 0]), abs(special_unitary[0, 0])) + t1 = cmath.phase(special_unitary[1, 1]) + t2 = cmath.phase(special_unitary[1, 0]) + alpha = t1 + t2 + gamma = t1 - t2 + else: + raise Exception("ZYZ decomposition only applies to single-qubit gate.") + return gamma, beta, alpha, global_phase + + +# # # # # # # Cosine-Sine Decomposition # # # # # # # # # # # # +def _thin_csd(q1, q2): + p = q1.shape[0] + print("the size of q1/q2: {}".format(p)) + + u1, c, v1 = np.linalg.svd(q1) + v1d = np.conjugate(v1).T + c = np.flip(c) + + cm = np.zeros((p, p), dtype=complex) + np.fill_diagonal(cm, c) + + u1 = np.fliplr(u1) + v1d = np.fliplr(v1d) + + q2 = q2 @ v1d + + # find the biggest index of c[k] <= 1/np.sqrt(2) + k = 0 + for i in range(1, p): + if c[i] <= 1 / np.sqrt(2): + k = i + + k = k + 1 + print("the k size: {}".format(k)) + + u2, _ = np.linalg.qr(q2[:, 0:k], mode='complete') + # u2, _= np.linalg.qr(q2, mode='complete') + print("the size of u2: {}".format(u2.shape)) + # print("the u2 matrix: {}".format(u2)) + s = np.conjugate(u2).T @ q2 + print("the size of s: {}".format(s.shape)) + # print("the s matrix: {}".format(np.real(s))) + + if k < p: + r2 = s[k:p, k:p] + print("the size of rs: {}".format(r2.shape)) + ut, ss, vt = np.linalg.svd(r2) + vtd = np.conjugate(vt).T + s[k:p, k:p] = np.diag(ss) + cm[:, k:p] = cm[:, k:p] @ vtd + u2[:, k:p] = u2[:, k:p] @ ut + v1d[:, k:p] = v1d[:, k:p] @ vtd + + w = cm[k:p, k:p] + z, r = np.linalg.qr(w, mode='complete') + cm[k:p, k:p] = r + u1[:, k:p] = u1[:, k:p] @ z + + for i in range(p): + if np.real(cm[i, i]) < 0: + cm[i, i] = -cm[i, i] + u1[:, i] = -u1[:, i] + if np.real(s[i, i]) < 0: + s[i, i] = -s[i, i] + u2[:, i] = -u2[:, i] + + return u1, u2, v1d, cm, s + + +def fat_csd(matrix): + """ + U = [U00, U01] = [u1 ][c s][v1 ] + [U10, U11] = [ u2][-s c][ v2] + """ + # print(matrix) + U00, U01, U10, U11 = mu.split_matrix(matrix) + + L0, L1, R0, cc, ss = _thin_csd(U00, U10) + R0 = np.conjugate(R0).T + ss = -ss + + # get the v2 + R1 = np.zeros_like(R0) + p = R1.shape[0] + for j in range(p): + if np.abs(ss[j, j]) > np.abs(cc[j, j]): + L0d = np.conjugate(L0).T + tmp = L0d @ U01 + R1[j, :] = tmp[j, :] / ss[j, j] + else: + L1d = np.conjugate(L1).T + tmp = L1d @ U11 + R1[j, :] = tmp[j, :] / cc[j, j] + + assert mu.is_approx(L0 @ cc @ R0, U00) + assert mu.is_approx(-L1 @ ss @ R0, U10) + assert mu.is_approx(L0 @ ss @ R1, U01) + assert mu.is_approx(L1 @ cc @ R1, U11) + + zeros_m = np.zeros_like(L0) + L = mu.stack_matrices(L0, zeros_m, zeros_m, L1) + D = mu.stack_matrices(cc, ss, -ss, cc) + R = mu.stack_matrices(R0, zeros_m, zeros_m, R1) + assert mu.is_approx(matrix, L @ D @ R) + + return L0, L1, R0, R1, cc, ss # L0, L1 is unitary + + +# # # # # # # # # # # # # # # # # # # # # # # # # # +def demultiplexing(u1, u2): + """ + U = [U1 0] = [V 0][D 0][W 0] + [0 U2] = [0 V][0 D*][0 W] + """ + assert mu.is_unitary(u1) + assert mu.is_unitary(u2) + + u2_dg = np.conjugate(u2).T + + [d2, v] = scipy.linalg.schur(u1 @ u2_dg) + assert mu.is_diagonal(d2) + assert mu.is_approx(v @ d2 @ np.conjugate(v).T, u1 @ u2_dg) + + d_tmp = np.sqrt(np.diag(d2)) + d = np.diag(d_tmp) + + assert mu.is_approx(d @ d, d2) + v_dg = np.conjugate(v).T + + w = d @ v_dg @ u2 + + assert mu.is_approx(u1, v @ d @ w) + assert mu.is_approx(u2, v @ np.conjugate(d).T @ w) + + zm = np.zeros_like(w) + vv = mu.stack_matrices(v, zm, zm, v) + dd = mu.stack_matrices(d, zm, zm, np.conjugate(d).T) + ww = mu.stack_matrices(w, zm, zm, w) + uu = mu.stack_matrices(u1, zm, zm, u2) + + assert mu.is_approx(vv @ dd @ ww, uu) + assert mu.is_unitary(v) + return v, d, w + + +def _graycode(n): + for i in range(1 << n): + gray = i ^ (i >> 1) + yield "{0:0{1}b}".format(gray, n) + + +def get_multi_control_index(n): + gray_codes = list(_graycode(n)) + size = 2 ** n + + index_list = [] + for i in range(size): + str1 = gray_codes[i] + str2 = gray_codes[(i + 1) % size] + + tmp = [k for k in range(len(str1)) if str1[k] != str2[k]] + assert len(tmp) == 1 + index_list.append(tmp[0]) + + return index_list + + +def genMk_table(nqubits): + """ + TODO: add docstring + """ + import re + + def bin2gray(num): + return num ^ int((num >> 1)) + + def genMk(k): + + Mk = np.zeros((2 ** k, 2 ** k)) + for i in range(2 ** k): + for j in range(2 ** k): + p = i & bin2gray(j) + strbin = "{0:b}".format(p) + tmp = [m.start() for m in re.finditer('1', strbin)] + Mk[i, j] = (-1) ** len(tmp) + + return Mk + + genMk_lookuptable = [] + for n in range(1, nqubits + 1): + tmp = genMk(n) + genMk_lookuptable.append(tmp) + + return genMk_lookuptable diff --git a/src/quafu/elements/quantum_element.py b/src/quafu/elements/quantum_element.py deleted file mode 100644 index 7dc3c4e..0000000 --- a/src/quafu/elements/quantum_element.py +++ /dev/null @@ -1,281 +0,0 @@ -# This is the file for abstract quantum gates class -from typing import Union, Callable, List, Tuple, Iterable, Any, Optional -import numpy as np -from functools import reduce -import copy - -class Barrier(object): - def __init__(self, pos): - self.name = "barrier" - self.__pos = pos - self.symbol = "||" - - @property - def pos(self): - return self.__pos - - @pos.setter - def pos(self, pos): - self.__pos = pos - - def __repr__(self): - return f"{self.__class__.__name__}" - -class Delay(object): - def __init__(self, pos : int, duration : int, unit="ns"): - self.name = "delay" - if isinstance(duration, int): - self.duration = duration - else: - raise TypeError("duration must be int") - self.unit=unit - self.pos=pos - self.symbol = "Delay(%d%s)" %(duration, unit) - - def __repr__(self): - return f"{self.__class__.__name__}" - -class XYResonance(object): - def __init__(self, qs : int, qe : int, duration : int, unit="ns"): - self.name = "XY" - if isinstance(duration, int): - self.duration = duration - else: - raise TypeError("duration must be int") - self.unit=unit - self.pos=list(range(qs, qe+1)) - self.symbol = "XY(%d%s)" %(duration, unit) - -class QuantumGate(object): - def __init__(self, name: str, pos: Union[int, List[int]], paras: Union[None,float, List[float]], matrix): - self.name = name - self.pos = pos - self.paras = paras - self.matrix = matrix - - if paras: - if isinstance(paras, Iterable): - self.symbol = "%s(" %self.name + ",".join(["%.3f" %para for para in self.paras]) + ")" - else: - self.symbol = "%s(%.3f)" % (self.name, paras) - else: - self.symbol = "%s" %self.name - - @property - def name(self): - return self.__name - - @name.setter - def name(self, _name): - self.__name = _name - - @property - def pos(self): - return self.__pos - - @pos.setter - def pos(self, _pos): - self.__pos = _pos - - @property - def paras(self): - return self.__paras - - @paras.setter - def paras(self, _paras): - self.__paras = _paras - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix): - self._matrix = matrix - - def __str__(self): - properties_names = ['pos', 'paras', 'matrix'] - properties_values = [getattr(self, x) for x in properties_names] - return "%s:\n%s" % (self.__class__.__name__, '\n'.join( - [f"{x} = {repr(properties_values[i])}" for i, x in enumerate(properties_names)])) - - def __repr__(self): - return f"{self.__class__.__name__}" - - -class SingleQubitGate(QuantumGate): - def __init__(self, name: str, pos: int, paras, matrix): - super().__init__(name, pos, paras=paras, matrix=matrix) - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix): - if isinstance(matrix, (np.ndarray, List)): - if np.shape(matrix) == (2, 2): - self._matrix = np.asarray(matrix, dtype=complex) - else: - raise ValueError(f'`{self.__class__.__name__}.matrix.shape` must be (2, 2)') - elif isinstance(matrix, type(None)): - self._matrix = matrix - else: - raise TypeError("Unsupported `matrix` type") - - def get_targ_matrix(self, reverse_order=False): - return self.matrix - -class FixedSingleQubitGate(SingleQubitGate): - def __init__(self, name, pos, matrix): - super().__init__(name, pos, paras=None, matrix=matrix) - - -class ParaSingleQubitGate(SingleQubitGate): - def __init__(self, name, pos, paras, matrix): - super().__init__(name, pos, paras=paras, matrix=matrix) - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix): - if isinstance(matrix, Callable): - self._matrix = matrix(self.paras) - elif isinstance(matrix, (np.ndarray, List)): - if np.shape(matrix) == (2, 2): - self._matrix = np.asarray(matrix, dtype=complex) - else: - raise ValueError(f'`{self.__class__.__name__}.matrix.shape` must be (2, 2)') - elif isinstance(matrix, type(None)): - self._matrix = matrix - else: - raise TypeError("Unsupported `matrix` type") - - -class MultiQubitGate(QuantumGate): - def __init__(self, name: str, pos: List[int], paras, matrix): - super().__init__(name, pos, paras=paras, matrix=matrix) - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix): - if isinstance(matrix, np.ndarray): - self._matrix = matrix - elif isinstance(matrix, List): - self._matrix = np.array(matrix, dtype=complex) - else: - raise TypeError("Unsupported `matrix` type") - - self.reorder_matrix() - - def reorder_matrix(self): - """Reorder the input sorted matrix to the pos order """ - qnum = len(self.pos) - dim = 2**qnum - inds = np.argsort(self.pos) - inds = np.concatenate([inds, inds+qnum]) - tensorm = self._matrix.reshape([2]*2*qnum) - self._matrix = np.transpose(tensorm, inds).reshape([dim, dim]) - - def get_targ_matrix(self, reverse_order=False): - targ_matrix = self._matrix - if reverse_order and (len(self.pos) > 1): - qnum = len(self.pos) - order = np.array(range(len(self.pos))[::-1]) - order = np.concatenate([order, order+qnum]) - dim = 2**qnum - tensorm = self._matrix.reshape([2]*2*qnum) - targ_matrix = np.transpose(tensorm, order).reshape([dim, dim]) - - return targ_matrix - - -class FixedMultiQubitGate(MultiQubitGate): - def __init__(self, name: str, pos: List[int], matrix): - super().__init__(name, pos, paras=None, matrix=matrix) - - -class ParaMultiQubitGate(MultiQubitGate): - def __init__(self, name, pos, paras, matrix): - super().__init__(name, pos, paras, matrix=matrix) - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix): - if isinstance(matrix, Callable): - self._matrix = matrix(self.paras) - self.reorder_matrix() - elif isinstance(matrix, (np.ndarray, List)): - self._matrix = matrix - self.reorder_matrix() - else: - raise TypeError("Unsupported `matrix` type") - - -class ControlledGate(MultiQubitGate): - """ Controlled gate class, where the matrix act non-trivallly on target qubits""" - def __init__(self, name, targe_name, ctrls: List[int], targs: List[int], paras, matrix): - self.ctrls = ctrls - self.targs = targs - self.targ_name = targe_name - super().__init__(name, ctrls+targs, paras, matrix) - self._targ_matrix = matrix - - if paras: - if isinstance(paras, Iterable): - self.symbol = "%s(" %self.targ_name + ",".join(["%.3f" %para for para in self.paras]) + ")" - else: - self.symbol = "%s(%.3f)" % (self.targ_name, paras) - else: - self.symbol = "%s" %self.targ_name - - - @property - def matrix(self): - return self._matrix - - @matrix.setter - def matrix(self, matrix : Union[np.ndarray, Callable]): - targ_dim = 2**(len(self.targs)) - qnum = len(self.pos) - dim = 2**(qnum) - if isinstance(matrix, Callable): - matrix = matrix(self.paras) - - if matrix.shape[0] != targ_dim: - raise ValueError("Dimension dismatch") - else: - self._matrix = np.zeros((dim , dim), dtype=complex) - control_dim = 2**len(self.pos) - targ_dim - for i in range(control_dim): - self._matrix[i, i] = 1. - - self._matrix[control_dim:, control_dim:] = matrix - self.reorder_matrix() - self._targ_matrix = matrix - - def get_targ_matrix(self, reverse_order=False): - return self._targ_matrix - -class ControlledU(ControlledGate): - def __init__(self, name, ctrls: List[int], U: Union[SingleQubitGate, MultiQubitGate]): - self.targ_gate = U - targs = U.pos - if isinstance(targs, int): - targs = [targs] - - super().__init__(name, U.name, ctrls, targs, U.paras, matrix=self.targ_gate.matrix) - - def get_targ_matrix(self, reverse_order=False): - return self.targ_gate.get_targ_matrix(reverse_order) - - - diff --git a/src/quafu/elements/quantum_element/__init__.py b/src/quafu/elements/quantum_element/__init__.py new file mode 100644 index 0000000..7dd3936 --- /dev/null +++ b/src/quafu/elements/quantum_element/__init__.py @@ -0,0 +1,4 @@ +from .quantum_element import Barrier, Delay, XYResonance, Measure +from .quantum_gate import QuantumGate, SingleQubitGate, MultiQubitGate, \ + ParaSingleQubitGate, ParaMultiQubitGate, ControlledGate, FixedSingleQubitGate, FixedMultiQubitGate +from .instruction import Instruction diff --git a/src/quafu/elements/quantum_element/instruction.py b/src/quafu/elements/quantum_element/instruction.py new file mode 100644 index 0000000..1ccb6db --- /dev/null +++ b/src/quafu/elements/quantum_element/instruction.py @@ -0,0 +1,75 @@ +# (C) Copyright 2023 Beijing Academy of Quantum Information Sciences +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from abc import ABC, abstractmethod +from typing import Union, List + + +PosType = Union[int, List[int]] +ParaType = dict[str, Union[float, int]] + + +class Instruction(ABC): + ins_classes = {} + + @property + @abstractmethod + def name(self) -> str: + raise NotImplementedError('name is not implemented for %s' % self.__class__.__name__ + + ', this should never happen.') + + @name.setter + def name(self, _): + import warnings + warnings.warn("Invalid assignment, names of standard instructions are not alterable.") + + @classmethod + def register_ins(cls, subclass, name: str = None): + assert issubclass(subclass, cls) + + if name is None: + name = subclass.name + if name in cls.ins_classes: + raise ValueError(f"Name {name} already exists.") + cls.ins_classes[name] = subclass + + def __init__(self, pos: PosType, paras: ParaType = None, *args, **kwargs): + self.pos = pos + self.paras = paras + + # @_ins_id.setter + # def sd_name(self, name: str): + # if self.sd_name is None: + # self.sd_name = name + # else: + # import warnings + # warnings.warn(message='Invalid assignment, names of standard ' + # 'instructions are not alterable.') + + # def to_dag_node(self): + # name = self.get_ins_id() + # label = self.__repr__() + # + # pos = self.pos + # paras = self.paras + # paras = {} if paras is None else paras + # duration = paras.get('duration', None) + # unit = paras.get('unit', None) + # channel = paras.get('channel', None) + # time_func = paras.get('time_func', None) + # + # return InstructionNode(name, pos, paras, duration, unit, channel, time_func, label) + + + diff --git a/src/quafu/elements/quantum_element/pulses/__init__.py b/src/quafu/elements/quantum_element/pulses/__init__.py new file mode 100644 index 0000000..34be485 --- /dev/null +++ b/src/quafu/elements/quantum_element/pulses/__init__.py @@ -0,0 +1,11 @@ +# !/usr/bin/env python +# -*- coding:utf-8 -*- + +# @File: __init__.py.py +# @Version: 1.0 +# @Author: Yun-Hao Shi +# @Email: yhshi@iphy.ac.cn +# @Reminders: |0> has been in the past, |1> is still in the future +# @License: Copyright (c) 2023, IOP, CAS + +from .quantum_pulse import QuantumPulse diff --git a/src/quafu/elements/quantum_element/pulses/quantum_pulse.py b/src/quafu/elements/quantum_element/pulses/quantum_pulse.py new file mode 100644 index 0000000..1ed07c0 --- /dev/null +++ b/src/quafu/elements/quantum_element/pulses/quantum_pulse.py @@ -0,0 +1,236 @@ +from abc import ABC, abstractmethod +from copy import deepcopy +from typing import Union, Optional + +import matplotlib.pyplot as plt +import numpy as np + +from quafu.elements.quantum_element.instruction import Instruction, PosType + +TimeType = Union[np.ndarray, float, int] + + +class QuantumPulse(Instruction, ABC): + pulse_classes = {} + + def __init__(self, + pos: PosType, + paras: list, + duration: Union[float, int], + unit: str, + channel: str, + ): + """ + Quantum Pulse for generating a quantum gate. + + Args: + pos (int): Qubit position. + paras (list): Parameters of the pulse. + duration (float, int): Pulse duration. + unit (str): Duration unit. + """ + super().__init__() + self.pos = pos + self.paras = paras + self.duration = duration + self.unit = unit + if channel in ["XY", "Z"]: + self.channel = channel + else: + raise ValueError("channel must be 'XY' or 'Z'") + + @property + def symbol(self): + return "%s(%d%s, %s)" % (self.name, self.duration, self.unit, self.channel) + + @abstractmethod + def time_func(self, t: Union[np.ndarray, float, int], **kwargs): + """ + Return the pulse data. + + Args: + t (np.ndarray, float, int): Time list. + kwargs (dict): Keyword arguments for the pulse. + """ + pass + + @classmethod + def register_pulse(cls, subclass, name: str = None): + assert issubclass(subclass, cls) + + if name is None: + name = subclass.name + if name in cls.pulse_classes: + raise ValueError(f"Name {name} already exists.") + cls.pulse_classes[name] = subclass + Instruction.register_ins(subclass, name) + + def __repr__(self): + return self.__str__() + + def __str__(self): + symbol = "%s(%d%s" % (self.name, self.duration, self.unit) + for para in self.paras: + symbol += ", %s" % para + symbol += ", %s" % self.channel + symbol += ")" + return symbol + + def __call__(self, + t: TimeType, + shift: Union[float, int] = 0., + offset: Union[float, int] = 0., + args: dict = None + ): + """ + Return pulse data. + + Args: + t (np.ndarray, float, int): Time list. + shift (float, int): Time shift. + offset (float, int): Pulse amplitude offset. + """ + window = np.logical_and(0 <= t, t <= self.duration) + if args is None: + return window * self.time_func(t - shift) + else: + return window * self.time_func(t - shift, **args) + + def __copy__(self): + """ Return a deepcopy of the pulse """ + return deepcopy(self) + + def to_qasm(self): + return self.__str__() + " q[%d]" % self.pos + + def plot(self, + t: Optional[np.ndarray] = None, + shift: Union[float, int] = 0., + offset: Union[float, int] = 0., + plot_real: bool = True, + plot_imag: bool = True, + fig=None, + ax=None, + **plot_kws): + """ + Plot the pulse waveform. + + Args: + t (np.ndarray): Time list of the plot. + shift (float, int): Time shift of the pulse. + offset (float, int): Offset of the pulse. + plot_real (bool): Plot real of the pulse. + plot_imag (bool): Plot imag of the pulse. + fig (Figure): Figure of the plot. + ax (Axes): Axes of the plot. + plot_kws (dict, optional): Plot kwargs of `ax.plot`. + """ + if t is None: + t = np.linspace(0, self.duration, 101) + if ax is None: + fig, ax = plt.subplots(1, 1, num=fig) + pulse_data = self(t, shift=shift, offset=offset) + if plot_real: + ax.plot(np.real(pulse_data), label='real', **plot_kws) + if plot_imag: + ax.plot(np.imag(pulse_data), label='imag', **plot_kws) + ax.set_xlabel("Time (%s)" % self.unit) + ax.set_ylabel("Pulse Amp (a.u.)") + ax.legend() + plt.show() + + def set_pos(self, pos: int): + """ Set qubit position """ + self.pos = pos + return self + + def set_unit(self, unit="ns"): + """ Set duration unit """ + self.unit = unit + return self + + +class RectPulse(QuantumPulse): + name = "rect" + + def __init__(self, pos, amp, duration, unit, channel): + self.amp = amp + + super().__init__(pos, [amp], duration, unit, channel) + + def time_func(self, t: Union[np.ndarray, float, int], **kwargs): + """ rect_time_func """ + amp_ = kwargs["amp"] + return amp_ * np.ones(np.array(t).shape) + + def __call__(self, + t: TimeType, + shift: Union[float, int] = 0, + offset: Union[float, int] = 0, + *args, + **kwargs): + args = {"amp": self.amp} + return super().__call__(t, shift, offset, args) + + +class FlattopPulse(QuantumPulse): + name = "flattop" + + def __init__(self, pos, amp, fwhm, duration, unit, channel): + self.amp = amp + self.fwhm = fwhm + + super().__init__(pos, [amp, fwhm], duration, unit, channel) + + def time_func(self, t, **kws): + """ flattop_time_func """ + from scipy.special import erf + + amp_, fwhm_ = kws["amp"], kws["fwhm"] + sigma_ = fwhm_ / (2 * np.sqrt(np.log(2))) + return amp_ * (erf((self.duration - t) / sigma_) + erf(t / sigma_) - 1.) + + def __call__(self, + t: TimeType, + shift: Union[float, int] = 0, + offset: Union[float, int] = 0, + *args, + **kwargs): + args = {"amp": self.amp, "fwhm": self.fwhm} + return super().__call__(t, shift, offset, args) + + +class GaussianPulse(QuantumPulse): + name = "gaussian" + + def __init__(self, pos, amp, fwhm, phase, duration, unit, channel): + self.amp = amp + if fwhm is None: + self.fwhm = 0.5 * duration + else: + self.fwhm = fwhm + self.phase = phase + + super().__init__(pos, [amp, fwhm, phase], duration, unit, channel) + + def time_func(self, t, **kws): + """ gaussian_time_func """ + amp_, fwhm_, phase_ = kws["amp"], kws["fwhm"], kws["phase"] + # start: t = 0, center: t = 0.5 * duration, end: t = duration + sigma_ = fwhm_ / np.sqrt(8 * np.log(2)) # fwhm to std. deviation + return amp_ * np.exp( + -(t - 0.5 * self.duration) ** 2 / (2 * sigma_ ** 2) + 1j * phase_) + + def __call__(self, + t: TimeType, + shift: Union[float, int] = 0, + offset: Union[float, int] = 0, + *args, + **kwargs): + args = {"amp": self.amp, "fwhm": self.fwhm, "phase": self.phase} + return super().__call__(t, shift, offset, args) + + +QuantumPulse.register_pulse(RectPulse) +QuantumPulse.register_pulse(FlattopPulse) +QuantumPulse.register_pulse(GaussianPulse) diff --git a/src/quafu/elements/quantum_element/quantum_element.py b/src/quafu/elements/quantum_element/quantum_element.py new file mode 100644 index 0000000..f182fbe --- /dev/null +++ b/src/quafu/elements/quantum_element/quantum_element.py @@ -0,0 +1,93 @@ +# (C) Copyright 2023 Beijing Academy of Quantum Information Sciences +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from quafu.elements.quantum_element.instruction import Instruction + +""" +Base classes for ALL kinds of possible instructions on superconducting +quantum circuits. +""" + + +class Barrier(Instruction): + name = "barrier" + + def __init__(self, pos): + super().__init__(pos) + self.symbol = "||" + + @property + def pos(self): + return self.__pos + + @pos.setter + def pos(self, pos): + self.__pos = pos + + def __repr__(self): + return f"{self.__class__.__name__}" + + def to_qasm(self): + return "barrier " + ",".join(["q[%d]" % p for p in range(min(self.pos), max(self.pos) + 1)]) + + +class Delay(Instruction): + name = "delay" + + def __init__(self, pos: int, duration: int, unit="ns"): + if isinstance(duration, int): + self.duration = duration + else: + raise TypeError("duration must be int") + super().__init__(pos) + self.unit = unit + self.symbol = "Delay(%d%s)" % (duration, unit) + + def __repr__(self): + return f"{self.__class__.__name__}" + + def to_qasm(self): + return "delay(%d%s) q[%d]" % (self.duration, self.unit, self.pos) + + +class XYResonance(Instruction): + name = "XY" + + def __init__(self, qs: int, qe: int, duration: int, unit="ns"): + if isinstance(duration, int): + self.duration = duration + else: + raise TypeError("duration must be int") + super().__init__(list(range(qs, qe + 1))) + self.unit = unit + self.symbol = "XY(%d%s)" % (duration, unit) + + def to_qasm(self): + return "xy(%d%s) " % (self.duration, self.unit) + ",".join( + ["q[%d]" % p for p in range(min(self.pos), max(self.pos) + 1)]) + + +class Measure(Instruction): + name = "measure" + + def __init__(self, bitmap: dict): + super().__init__(list(bitmap.keys())) + self.qbits = bitmap.keys() + self.cbits = bitmap.values() + + +Instruction.register_ins(Barrier) +Instruction.register_ins(Delay) +Instruction.register_ins(XYResonance) +Instruction.register_ins(Measure) diff --git a/src/quafu/elements/quantum_element/quantum_gate.py b/src/quafu/elements/quantum_element/quantum_gate.py new file mode 100644 index 0000000..5b60236 --- /dev/null +++ b/src/quafu/elements/quantum_element/quantum_gate.py @@ -0,0 +1,183 @@ +from abc import ABC, abstractmethod +from typing import List, Union, Iterable + +import numpy as np + +from quafu.elements.quantum_element.instruction import Instruction, PosType + + +def reorder_matrix(matrix: np.ndarray, pos: List): + """Reorder the input sorted matrix to the pos order """ + qnum = len(pos) + dim = 2 ** qnum + inds = np.argsort(pos) + inds = np.concatenate([inds, inds + qnum]) + tensorm = matrix.reshape([2] * 2 * qnum) + return np.transpose(tensorm, inds).reshape([dim, dim]) + + +class QuantumGate(Instruction): + gate_classes = {} + + def __init__(self, + pos: PosType, + paras: Union[float, List[float]] = None, + ): + super().__init__(pos, paras) + + if paras is not None: + if isinstance(paras, Iterable): + self.symbol = "%s(" % self.name + ",".join(["%.3f" % para for para in self.paras]) + ")" + else: + self.symbol = "%s(%.3f)" % (self.name, paras) + else: + self.symbol = "%s" % self.name + + @property + @abstractmethod + def matrix(self): + raise NotImplementedError("Matrix is not implemented for %s" % self.__class__.__name__ + + ", this should never happen.") + + @classmethod + def register_gate(cls, subclass, name: str = None): + assert issubclass(subclass, cls) + + if name is None: + name = subclass.name + if name in cls.gate_classes: + raise ValueError(f"Name {name} already exists.") + cls.gate_classes[name] = subclass + Instruction.register_ins(subclass, name) + + def __str__(self): + properties_names = ['pos', 'paras', 'matrix'] + properties_values = [getattr(self, x) for x in properties_names] + return "%s:\n%s" % (self.__class__.__name__, '\n'.join( + [f"{x} = {repr(properties_values[i])}" for i, x in enumerate(properties_names)])) + + def __repr__(self): + return f"{self.__class__.__name__}" + + def to_qasm(self): + qstr = "%s" % self.name.lower() + + if self.paras: + if isinstance(self.paras, Iterable): + qstr += "(" + ",".join(["%s" % para for para in self.paras]) + ")" + else: + qstr += "(%s)" % self.paras + qstr += " " + if isinstance(self.pos, Iterable): + qstr += ",".join(["q[%d]" % p for p in self.pos]) + else: + qstr += "q[%d]" % self.pos + + return qstr + + +class SingleQubitGate(QuantumGate, ABC): + def __init__(self, pos: int, paras: float = None): + super().__init__(pos, paras=paras) + + def get_targ_matrix(self): + return self.matrix + + +class MultiQubitGate(QuantumGate, ABC): + def __init__(self, pos: List[int], paras: float = None): + super().__init__(pos, paras) + + def get_targ_matrix(self, reverse_order=False): + targ_matrix = self.matrix + + if reverse_order and (len(self.pos) > 1): + qnum = len(self.pos) + dim = 2 ** qnum + order = np.array(range(len(self.pos))[::-1]) + order = np.concatenate([order, order + qnum]) + tensorm = targ_matrix.reshape([2] * 2 * qnum) + targ_matrix = np.transpose(tensorm, order).reshape([dim, dim]) + return targ_matrix + + +class ParaSingleQubitGate(SingleQubitGate, ABC): + def __init__(self, pos, paras: float): + if paras is None: + raise ValueError("`paras` can not be None for ParaSingleQubitGate") + elif isinstance(paras, int): + paras = float(paras) + print(paras) + + if not isinstance(paras, float): + raise TypeError(f"`paras` must be float or int for ParaSingleQubitGate, instead of {type(paras)}") + super().__init__(pos, paras=paras) + + +class FixedSingleQubitGate(SingleQubitGate, ABC): + def __init__(self, pos): + super().__init__(pos=pos, paras=None) + + +class FixedMultiQubitGate(MultiQubitGate, ABC): + def __init__(self, pos: List[int]): + super().__init__(pos=pos, paras=None) + + +class ParaMultiQubitGate(MultiQubitGate, ABC): + def __init__(self, pos, paras): + if paras is None: + raise ValueError("`paras` can not be None for ParaMultiQubitGate") + super().__init__(pos, paras) + + +class ControlledGate(MultiQubitGate, ABC): + """ Controlled gate class, where the matrix act non-trivaly on target qubits""" + + def __init__(self, targe_name, ctrls: List[int], targs: List[int], paras, tar_matrix): + super().__init__(ctrls + targs, paras) + self.ctrls = ctrls + self.targs = targs + self.targ_name = targe_name + self._targ_matrix = tar_matrix + + # set matrix + # TODO: change matrix according to control-type 0/1 + c_n, t_n, n = self.ct_nums + targ_dim = 2 ** t_n + dim = 2 ** n + ctrl_dim = dim - targ_dim + self._matrix = np.eye(dim, dtype=complex) + self._matrix[ctrl_dim:, ctrl_dim:] = tar_matrix + self._matrix = reorder_matrix(self._matrix, self.pos) + + if paras is not None: + if isinstance(paras, Iterable): + self.symbol = "%s(" % self.targ_name + ",".join(["%.3f" % para for para in self.paras]) + ")" + else: + self.symbol = "%s(%.3f)" % (self.targ_name, paras) + else: + self.symbol = "%s" % self.targ_name + + @property + def matrix(self): + # TODO: update matrix when paras of controlled-gate changed + return self._matrix + + @property + def ct_nums(self): + targ_num = len(self.targs) + ctrl_num = len(self.ctrls) + num = targ_num + ctrl_num + return ctrl_num, targ_num, num + + def get_targ_matrix(self, reverse_order=False): + targ_matrix = self._targ_matrix + if reverse_order and (len(self.targs) > 1): + qnum = len(self.targs) + order = np.array(range(len(self.targs))[::-1]) + order = np.concatenate([order, order + qnum]) + dim = 2 ** qnum + tensorm = targ_matrix.reshape([2] * 2 * qnum) + targ_matrix = np.transpose(tensorm, order).reshape([dim, dim]) + return targ_matrix diff --git a/src/quafu/exceptions/__init__.py b/src/quafu/exceptions/__init__.py new file mode 100644 index 0000000..0490b03 --- /dev/null +++ b/src/quafu/exceptions/__init__.py @@ -0,0 +1,2 @@ +from .quafu_error import * + diff --git a/src/quafu/exceptions/circuit_error.py b/src/quafu/exceptions/circuit_error.py new file mode 100644 index 0000000..8dda42f --- /dev/null +++ b/src/quafu/exceptions/circuit_error.py @@ -0,0 +1,19 @@ +from .quafu_error import CircuitError + + +# TODO +class IndexOutOfRangeError(CircuitError): + """ """ + pass + + +class InvalidParaError(CircuitError): + """ """ + pass + + +class UnsupportedYet(CircuitError): + """ Un-supported instructions occurs""" + pass + + diff --git a/src/quafu/exceptions.py b/src/quafu/exceptions/quafu_error.py similarity index 80% rename from src/quafu/exceptions.py rename to src/quafu/exceptions/quafu_error.py index 939989c..221a7f8 100644 --- a/src/quafu/exceptions.py +++ b/src/quafu/exceptions/quafu_error.py @@ -1,27 +1,31 @@ - -""" -Exceptions for errors raised while building circuit. -""" - -class QuafuError(Exception): - """Base class for errors raised by Quafu.""" - - def __init__(self, *message): - """Set the error message.""" - super().__init__(" ".join(message)) - self.message = " ".join(message) - - def __str__(self): - """Return the message.""" - return repr(self.message) - -class CircuitError(QuafuError): - """Exceptions for errors raised while building circuit.""" - pass - -class ServerError(QuafuError): - pass - -class CompileError(QuafuError): - pass - +""" +Exceptions for errors raised while building circuit. +""" + + +class QuafuError(Exception): + """Base class for errors raised by Quafu.""" + + def __init__(self, *message): + """Set the error message.""" + super().__init__(" ".join(message)) + self.message = " ".join(message) + + def __str__(self): + """Return the message.""" + return repr(self.message) + + +class CircuitError(QuafuError): + """Exceptions for errors raised while building circuit.""" + pass + + +class ServerError(QuafuError): + """ Exceptions for errors raised while connecting to server.""" + pass + + +class CompileError(QuafuError): + """ Exceptions for errors raised while compiling circuit. """ + pass diff --git a/src/quafu/qfasm/__init__.py b/src/quafu/qfasm/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/quafu/qfasm/qfasm_convertor.py b/src/quafu/qfasm/qfasm_convertor.py new file mode 100644 index 0000000..9ccebce --- /dev/null +++ b/src/quafu/qfasm/qfasm_convertor.py @@ -0,0 +1,29 @@ +from .qfasm_parser import QfasmParser, QregNode +from quafu.dagcircuits.circuit_dag import node_to_gate +from quafu.dagcircuits.instruction_node import InstructionNode +from quafu.circuits.quantum_circuit import QuantumCircuit + + +def qasm_to_circuit(qasm): + parser = QfasmParser() + nodes = parser.parse(qasm) + + n = 0 + gates = [] + measures = {} + for node in nodes: + if isinstance(node, QregNode): + n = node.n + if isinstance(node, InstructionNode): + if node.name == "measure": + for q, c in zip(node.pos.keys(), node.pos.values()): + measures[q] = c + else: + gates.append(node_to_gate(node)) + + q = QuantumCircuit(n) + q.gates = gates + q.openqasm = qasm + q.measures = measures + return q + diff --git a/src/quafu/qfasm/qfasm_parser.py b/src/quafu/qfasm/qfasm_parser.py new file mode 100644 index 0000000..b62657c --- /dev/null +++ b/src/quafu/qfasm/qfasm_parser.py @@ -0,0 +1,270 @@ +import ply.yacc as yacc + +from quafu.dagcircuits.instruction_node import InstructionNode +from .qfasmlex import QfasmLexer + + +import numpy as np + + +class DeclarationNode(object): + pass + + +class QregNode(DeclarationNode): + def __init__(self, n): + self.n = n + + def __repr__(self): + return self.__str__() + + def __str__(self): + return "qreg[%d]" % self.n + + +class CregNode(DeclarationNode): + def __init__(self, n): + self.n = n + + def __repr__(self): + return self.__str__() + + def __str__(self): + return "creg[%d]" % self.n + + +class IncludeNode(DeclarationNode): + def __init__(self, filename): + self.file = filename + + def __repr__(self): + return self.__str__() + + def __str__(self): + return "include %s" % self.file + + +class OPENQASMNode(DeclarationNode): + def __init__(self, version): + self.version = version + + def __repr__(self): + return self.__str__() + + def __str__(self): + return "OPENQASM %.1f" % self.version + + +class QfasmParser(object): + tokens = QfasmLexer.tokens + + def __init__(self, debug=False): + self.parser = yacc.yacc(module=self, debug=debug) + self.parsed_nodes = [] + self.lexer = QfasmLexer() + + def parse(self, input: str): + self.parsed_nodes = self.parser.parse(input, lexer=QfasmLexer()) + return self.parsed_nodes + + def p_main_0(self, p): + """ + main : program + """ + p[0] = [p[1]] + + def p_main_1(self, p): + """ + main : main program + """ + p[1].append(p[2]) + p[0] = p[1] + + def p_program(self, p): + """ + program : instruction + | declaration + """ + p[0] = p[1] + + def p_declaration(self, p): + """ + declaration : openqasm + | include + | qreg + | creg + """ + p[0] = p[1] + + def p_openqasm(self, p): + """ + openqasm : OPENQASM FLOAT ';' + """ + p[0] = OPENQASMNode(p[2]) + + def p_include(self, p): + """ + include : INCLUDE STRING ';' + """ + p[0] = IncludeNode(p[2]) + + def p_qreg(self, p): # TODO:verify register name + """ + qreg : QREG bitreg ';' + """ + p[0] = QregNode(p[2]) + + def p_creg(self, p): + """ + creg : CREG bitreg ';' + """ + p[0] = CregNode(p[2]) + + def p_instruction(self, p): + """ + instruction : gate ';' + | pulse ';' + | measure ';' + """ + p[0] = p[1] + + def p_arg_list_0(self, p): + """ + arg_list : expression + """ + p[0] = [p[1]] + + def p_arg_list_1(self, p): + """ + arg_list : arg_list ',' expression + """ + p[1].append(p[3]) + p[0] = p[1] + + def p_gate_like_0(self, p): + ''' + gate : id qubit_list + ''' + p[0] = InstructionNode(p[1], p[2], None, None, None, "", None, "") + + def p_gate_like_1(self, p): + ''' + gate : id '(' arg_list ')' qubit_list + ''' + p[0] = InstructionNode(p[1], p[5], p[3], None, None, "", None, "") + + def p_pulse_like_0(self, p): + ''' + pulse : id '(' time ',' arg_list ')' qubit_list + ''' + p[0] = InstructionNode(p[1], p[7], p[5], p[3][0], p[3][1], "", None, "") + + def p_measure_0(self, p): + ''' + measure : MEASURE bitreg ASSIGN bitreg + ''' + p[0] = InstructionNode("measure", {p[2]: p[4]}, None, None, None, "", None, "") + + def p_pulse_like_1(self, p): + ''' + pulse : id '(' time ',' arg_list ',' channel ')' qubit_list + ''' + p[0] = InstructionNode(p[1], p[9], p[5], p[3][0], p[3][1], p[7], None, "") + + def p_bitreg(self, p): + """ + bitreg : id '[' INT ']' + """ + p[0] = p[3] + + def p_qubit_list_0(self, p): + """ + qubit_list : bitreg + """ + p[0] = [p[1]] + + def p_qubit_list_1(self, p): + """ + qubit_list : qubit_list ',' bitreg + """ + p[1].append(p[3]) + p[0] = p[1] + + def p_channel(self, p): + """ + channel : CHANNEL + """ + p[0] = p[1] + + def p_time(self, p): + """ + time : INT UNIT + """ + p[0] = (p[1], p[2]) + + def p_expression_none(self, p): + """ + expression : NONE + """ + p[0] = p[1] + + def p_expression_term(self, p): + 'expression : term' + p[0] = p[1] + + def p_expression_m(self, p): + """ + expression : '-' expression + """ + p[0] = -p[2] + + def p_term_factor(self, p): + """ + term : factor + """ + p[0] = p[1] + + def p_binary_operators(self, p): + '''expression : expression '+' term + | expression '-' term + term : term '*' factor + | term '/' factor''' + if p[2] == '+': + p[0] = p[1] + p[3] + elif p[2] == '-': + p[0] = p[1] - p[3] + elif p[2] == '*': + p[0] = p[1] * p[3] + elif p[2] == '/': + p[0] = p[1] / p[3] + + def p_factor_0(self, p): + ''' + factor : FLOAT + | INT + ''' + p[0] = p[1] + + def p_factor_1(self, p): + ''' + factor : '(' expression ')' + ''' + p[0] = p[2] + + def p_factor_pi(self, p): + ''' + factor : PI + ''' + p[0] = np.pi + + def p_id(self, p): + ''' id : ID''' + p[0] = p[1] + + def p_error(self, p): + if p: + print("Syntax error at token", p.type) + # Just discard the token and tell the parser it's okay. + self.parser.errok() + else: + print("Syntax error at EOF") diff --git a/src/quafu/qfasm/qfasmlex.py b/src/quafu/qfasm/qfasmlex.py new file mode 100644 index 0000000..797987c --- /dev/null +++ b/src/quafu/qfasm/qfasmlex.py @@ -0,0 +1,113 @@ +# Qfasm is modified OPENQASM 2.0. Currently it is mainly used to +# transfer circuit data to backends. Further it will +# support more instructions (classical or quantum) to enable +# interaction with quantum hardware + +import ply.lex as lex +import numpy as np + + +class QfasmLexer(object): + + def __init__(self): + self.build() + + def input(self, data): + self.data = data + self.lexer.input(data) + + def token(self): + ret = self.lexer.token() + return ret + + literals = r'=()[]{};<>,.+-/*^"' + + reserved = { + "creg": "CREG", + "qreg": "QREG", + "pi": "PI", + "measure": "MEASURE", + "include": "INCLUDE" + } + + tokens = [ + "FLOAT", + "INT", + "STRING", + "ASSIGN", + "MATCHES", + "ID", + "UNIT", + "CHANNEL", + "OPENQASM", + "NONE" + ] + list(reserved.values()) + + def t_FLOAT(self, t): + r"(([1-9]\d*\.\d*)|(0\.\d*[1-9]\d*))" + t.value = float(t.value) + return t + + def t_INT(self, t): + r"\d+" + t.value = int(t.value) + return t + + def t_STRING(self, t): + r"\"([^\\\"]|\\.)*\"" + return t + + def t_ASSIGN(self, t): + r"->" + return t + + def t_MATCHES(self, t): + r"==" + return t + + def t_UNIT(self, t): + r"ns|us" + return t + + def t_CHANNEL(self, t): + r"XY|Z" + return t + + def t_OPENQASM(self, t): + r"OPENQASM" + return t + + def t_NONE(self, t): + r"None" + return t + + def t_ID(self, t): + r"[a-z][a-zA-Z0-9_]*" + t.type = self.reserved.get(t.value, "ID") + return t + + t_ignore = " \t\r" + + def t_error(self, t): + print("Illegal character '%s'" % t.value[0]) + + def t_newline(self, t): + r"\n+" + t.lexer.lineno += len(t.value) + + def build(self, **kwargs): + self.lexer = lex.lex(module=self, **kwargs) + + def test(self, data): + self.lexer.input(data) + while True: + tok = self.lexer.token() + if not tok: + break + print(tok) + + +if __name__ == "__main__": + m = QfasmLexer() + m.build() + m.test("rx(21ns, pi) q[0], q[1]") diff --git a/src/quafu/results/results.py b/src/quafu/results/results.py index 647f97b..3f3ce69 100644 --- a/src/quafu/results/results.py +++ b/src/quafu/results/results.py @@ -1,14 +1,14 @@ - -import numpy as np -from functools import reduce import copy -import matplotlib.pyplot as plt from collections import OrderedDict + +import matplotlib.pyplot as plt from ..utils.basis import * + class Result(object): """Basis class for quantum results""" - pass + pass + class ExecResult(Result): """ @@ -20,8 +20,9 @@ class ExecResult(Result): taskid (int): Unique task id for the execute result. transpiled_circuit (QuantumCircuit): Quantum circuit transpiled on backend. """ + def __init__(self, input_dict, measures): - status_map = {0:"In Queue", 1:"Running", 2:"Completed", "Canceled":3, 4:"Failed"} + status_map = {0: "In Queue", 1: "Running", 2: "Completed", "Canceled": 3, 4: "Failed"} self.measures = measures self.task_status = status_map[input_dict["status"]] self.res = eval(input_dict['res']) @@ -29,8 +30,8 @@ def __init__(self, input_dict, measures): self.logicalq_res = {} cbits = list(self.measures.values()) for key, values in self.counts.items(): - newkey = "".join([key[i] for i in cbits]) - self.logicalq_res[newkey] = values + newkey = "".join([key[i] for i in cbits]) + self.logicalq_res[newkey] = values self.taskid = input_dict['task_id'] self.taskname = input_dict['task_name'] @@ -40,10 +41,9 @@ def __init__(self, input_dict, measures): self.transpiled_circuit.from_openqasm(self.transpiled_openqasm) self.measure_base = [] total_counts = sum(self.counts.values()) - self.probabilities = {} + self.probabilities = {} for key in self.counts: - self.probabilities[key] = self.counts[key]/total_counts - + self.probabilities[key] = self.counts[key] / total_counts def calculate_obs(self, pos): """ @@ -52,7 +52,7 @@ def calculate_obs(self, pos): Args: pos (list[int]): Positions of observalbes. """ - return measure_obs(pos, self.logicalq_res) + return measure_obs(pos, self.logicalq_res) def plot_probabilities(self): """ @@ -61,7 +61,7 @@ def plot_probabilities(self): bitstrs = list(self.probabilities.keys()) probs = list(self.probabilities.values()) plt.figure() - plt.bar(range(len(probs)), probs, tick_label = bitstrs) + plt.bar(range(len(probs)), probs, tick_label=bitstrs) plt.xticks(rotation=70) plt.ylabel("probabilities") @@ -75,6 +75,7 @@ class SimuResult(Result): probabilities (ndarray): Calculated probabilities on each bitstring. rho (ndarray): Simulated density matrix of measured qubits """ + def __init__(self, input, input_form): self.num = int(np.log2(input.shape[0])) if input_form == "density_matrix": @@ -83,9 +84,9 @@ def __init__(self, input, input_form): elif input_form == "probabilities": self.probabilities = input elif input_form == "state_vector": - self.state_vector = input - - def plot_probabilities(self, full: bool=False, reverse_basis: bool=False, sort: bool=None): + self.state_vector = input + + def plot_probabilities(self, full: bool = False, reverse_basis: bool = False, sort: bool = None): """ Plot the probabilities from simulated results, ordered in big endian convention. @@ -95,16 +96,15 @@ def plot_probabilities(self, full: bool=False, reverse_basis: bool=False, sort: sort: Sort the results by probabilities values. Can be `"ascend"` order or `"descend"` order. """ - from ..utils.basis import get_basis probs = self.probabilities inds = range(len(probs)) if not full: inds = np.where(self.probabilities > 1e-14)[0] probs = self.probabilities[inds] - basis=np.array([bin(i)[2:].zfill(self.num) for i in inds]) + basis = np.array([bin(i)[2:].zfill(self.num) for i in inds]) if reverse_basis: - basis=np.array([bin(i)[2:].zfill(self.num)[::-1] for i in inds]) + basis = np.array([bin(i)[2:].zfill(self.num)[::-1] for i in inds]) if sort == "ascend": orders = np.argsort(probs) @@ -115,7 +115,6 @@ def plot_probabilities(self, full: bool=False, reverse_basis: bool=False, sort: probs = probs[orders][::-1] basis = basis[orders][::-1] - plt.figure() plt.bar(inds, probs, tick_label=basis) plt.xticks(rotation=70) @@ -128,7 +127,7 @@ def calculate_obs(self, pos): "Calculate observables Z on input position using probabilities" inds = np.where(self.probabilities > 1e-14)[0] probs = self.probabilities[inds] - basis=np.array([bin(i)[2:].zfill(self.num) for i in inds]) + basis = np.array([bin(i)[2:].zfill(self.num) for i in inds]) res_reduced = dict(zip(basis, probs)) return measure_obs(pos, res_reduced) @@ -143,9 +142,10 @@ def intersec(a, b): inter.append(a[i]) aind.append(i) bind.append(j) - + return inter, aind, bind + def diff(a, b): diff = [] aind = [] @@ -153,7 +153,7 @@ def diff(a, b): if a[i] not in b: diff.append(a[i]) aind.append(i) - + return diff, aind @@ -164,12 +164,12 @@ def merge_measure(obslist): for obs in obslist: if len(measure_basis) == 0: measure_basis.append(obs) - targ_basis.append(len(measure_basis)-1) + targ_basis.append(len(measure_basis) - 1) else: added = 0 for mi in range(len(measure_basis)): measure_base = measure_basis[mi] - interset, intobsi, intbasei = intersec(obs[1], measure_base[1]) + interset, intobsi, intbasei = intersec(obs[1], measure_base[1]) diffset, diffobsi = diff(obs[1], measure_base[1]) if not len(interset) == 0: if all(np.array(list(obs[0]))[intobsi] == np.array(list(measure_base[0]))[intbasei]): @@ -185,8 +185,8 @@ def merge_measure(obslist): added = 1 break - if not added: + if not added: measure_basis.append(obs) - targ_basis.append(len(measure_basis)-1) + targ_basis.append(len(measure_basis) - 1) - return measure_basis, targ_basis \ No newline at end of file + return measure_basis, targ_basis diff --git a/src/quafu/simulators/default_simulator.py b/src/quafu/simulators/default_simulator.py index aee53a5..4db06c1 100644 --- a/src/quafu/simulators/default_simulator.py +++ b/src/quafu/simulators/default_simulator.py @@ -1,9 +1,9 @@ -#default circuit simulator for state vector +# default circuit simulator for state vector from typing import Iterable, List, Union from quafu.circuits.quantum_circuit import QuantumCircuit from ..results.results import SimuResult -from ..elements.quantum_element import Barrier, Delay, QuantumGate, SingleQubitGate, MultiQubitGate, XYResonance +from ..elements.quantum_element import Barrier, Delay, QuantumGate, SingleQubitGate, MultiQubitGate, XYResonance import numpy as np from functools import reduce from sparse import COO, SparseArray @@ -11,44 +11,46 @@ import copy + def global_op(gate: QuantumGate, global_qubits: List) -> coo_matrix: """Local operators to global operators""" - num = len(global_qubits) + num = len(global_qubits) if isinstance(gate, SingleQubitGate): - local_mat = coo_matrix(gate.matrix) + local_mat = coo_matrix(gate.matrix) pos = global_qubits.index(gate.pos) - local_mat = kron(kron(eye(2**pos), local_mat), eye(2**(num - pos-1))) + local_mat = kron(kron(eye(2 ** pos), local_mat), eye(2 ** (num - pos - 1))) return local_mat elif isinstance(gate, MultiQubitGate): - local_mat =coo_matrix(gate.matrix) + local_mat = coo_matrix(gate.matrix) pos = [global_qubits.index(p) for p in gate.pos] num_left = min(pos) num_right = num - max(pos) - 1 num_center = max(pos) - min(pos) + 1 - center_mat = kron(local_mat, eye(2**(num_center - len(pos)))) + center_mat = kron(local_mat, eye(2 ** (num_center - len(pos)))) origin_order = sorted(pos) - origin_order.extend([p for p in range(min(pos), max(pos)+1) if p not in pos]) + origin_order.extend([p for p in range(min(pos), max(pos) + 1) if p not in pos]) new_order = np.argsort(origin_order) center_mat = COO.from_scipy_sparse(center_mat) center_mat = permutebits(center_mat, new_order).to_scipy_sparse() - center_mat = kron(kron(eye(2**num_left), center_mat), eye(2**num_right)) + center_mat = kron(kron(eye(2 ** num_left), center_mat), eye(2 ** num_right)) return center_mat -def permutebits(mat: Union[SparseArray, np.ndarray], order : Iterable)\ - ->Union[SparseArray, np.ndarray]: +def permutebits(mat: Union[SparseArray, np.ndarray], order: Iterable) \ + -> Union[SparseArray, np.ndarray]: """permute qubits for operators or states""" num = len(order) order = np.array(order) r = len(mat.shape) - mat = np.reshape(mat, [2]*r*num) - order = np.concatenate([order+len(order)*i for i in range(r)]) + mat = np.reshape(mat, [2] * r * num) + order = np.concatenate([order + len(order) * i for i in range(r)]) mat = np.transpose(mat, order) - mat = np.reshape(mat, [2**num]*r) + mat = np.reshape(mat, [2 ** num] * r) return mat -def ptrace(psi, ind_A: List, diag: bool=True) -> np.ndarray: + +def ptrace(psi, ind_A: List, diag: bool = True) -> np.ndarray: """partial trace on a state vector""" num = int(np.log2(psi.shape[0])) order = copy.deepcopy(ind_A) @@ -56,18 +58,19 @@ def ptrace(psi, ind_A: List, diag: bool=True) -> np.ndarray: psi = permutebits(psi, order) if diag: - psi = np.abs(psi)**2 - psi = np.reshape(psi, [2**len(ind_A), 2**(num-len(ind_A))]) + psi = np.abs(psi) ** 2 + psi = np.reshape(psi, [2 ** len(ind_A), 2 ** (num - len(ind_A))]) psi = np.sum(psi, axis=1) - return psi + return psi else: - psi = np.reshape(psi, [2**len(ind_A), 2**(num-len(ind_A))]) + psi = np.reshape(psi, [2 ** len(ind_A), 2 ** (num - len(ind_A))]) rho = psi @ np.conj(np.transpose(psi)) return rho -def py_simulate(qc: QuantumCircuit, - state_ini: np.ndarray = np.array([]), - output: str="amplitudes") -> SimuResult: + +def py_simulate(qc: QuantumCircuit, + state_ini: np.ndarray = np.array([]), + output: str = "amplitudes") -> SimuResult: """Simulate quantum circuit Args: qc: quantum circuit need to be simulated. @@ -82,17 +85,15 @@ def py_simulate(qc: QuantumCircuit, used_qubits = qc.used_qubits num = len(used_qubits) if not state_ini: - psi = np.zeros(2**num) + psi = np.zeros(2 ** num) psi[0] = 1 else: psi = state_ini - for gate in qc.gates: - if not ((isinstance(gate, Delay)) or (isinstance(gate, Barrier)) or isinstance(gate, XYResonance)): + for gate in qc.gates: + if not ((isinstance(gate, Delay)) or (isinstance(gate, Barrier)) or isinstance(gate, XYResonance)): op = global_op(gate, used_qubits) psi = op @ psi return psi - - \ No newline at end of file diff --git a/src/quafu/simulators/simulator.py b/src/quafu/simulators/simulator.py index 8bb2dfe..7edd68d 100644 --- a/src/quafu/simulators/simulator.py +++ b/src/quafu/simulators/simulator.py @@ -1,24 +1,32 @@ - from typing import Union from .default_simulator import py_simulate, ptrace, permutebits from .qfvm import simulate_circuit, execute from quafu import QuantumCircuit from ..results.results import SimuResult import numpy as np -import time +from ..exceptions import QuafuError + -def simulate(qc : Union[QuantumCircuit, str], psi : np.ndarray= np.array([]), simulator:str="qfvm_circ", output: str="probabilities")-> SimuResult: +def simulate(qc: Union[QuantumCircuit, str], + psi: np.ndarray = np.array([]), + simulator: str = "qfvm_circ", + output: str = "probabilities", + use_gpu: bool = False, + use_custatevec: bool = False) -> SimuResult: """Simulate quantum circuit Args: qc: quantum circuit or qasm string that need to be simulated. psi : Input state vector - simulator:`"qfvm_circ"`: The high performance C++ circuit simulator. + simulator:`"qfvm_circ"`: The high performance C++ circuit simulator with optional GPU support. `"py_simu"`: Python implemented simulator by sparse matrix with low performace for large scale circuit. `"qfvm_qasm"`: The high performance C++ qasm simulator with limited gate set. output: `"probabilities"`: Return probabilities on measured qubits, ordered in big endian convention. `"density_matrix"`: Return reduced density_amtrix on measured qubits, ordered in big endian convention. `"state_vector`: Return original full statevector. The statevector returned by `qfvm` backend is ordered in little endian convention (same as qiskit), while `py_simu` backend is orderd in big endian convention. + use_gpu: Use the GPU version of `qfvm_circ` simulator. + use_custatevec: Use cuStateVec-based `qfvm_circ` simulator. The argument `use_gpu` must also be True. + Returns: SimuResult object that contain the results. """ @@ -30,24 +38,37 @@ def simulate(qc : Union[QuantumCircuit, str], psi : np.ndarray= np.array([]), si qasm = qc qc = QuantumCircuit(0) qc.from_openqasm(qasm) - + measures = [qc.used_qubits.index(i) for i in qc.measures.keys()] num = 0 if simulator == "qfvm_circ": - num = max(qc.used_qubits)+1 + num = max(qc.used_qubits) + 1 measures = list(qc.measures.keys()) - psi = simulate_circuit(qc, psi) - - elif simulator == "py_simu": + if use_gpu: + if use_custatevec: + try: + from .qfvm import simulate_circuit_custate + except ImportError: + raise QuafuError(" pyquafu is installed with cuquantum support") + psi = simulate_circuit_custate(qc, psi) + else: + try: + from .qfvm import simulate_circuit_gpu + except ImportError: + raise QuafuError("you are not using the GPU version of pyquafu") + psi = simulate_circuit_gpu(qc, psi) + else: + psi = simulate_circuit(qc, psi) + + elif simulator == "py_simu": psi = py_simulate(qc, psi) elif simulator == "qfvm_qasm": num = qc.num measures = list(qc.measures.keys()) - psi = execute(qasm) + psi = execute(qasm) else: raise ValueError("invalid circuit") - if output == "density_matrix": if simulator in ["qfvm_circ", "qfvm_qasm"]: psi = permutebits(psi, range(num)[::-1]) @@ -60,10 +81,10 @@ def simulate(qc : Union[QuantumCircuit, str], psi : np.ndarray= np.array([]), si psi = permutebits(psi, range(num)[::-1]) probabilities = ptrace(psi, measures) probabilities = permutebits(probabilities, list(qc.measures.values())) - return SimuResult(probabilities, output) - + return SimuResult(probabilities, output) + elif output == "state_vector": return SimuResult(psi, output) else: - raise ValueError("output should in be 'density_matrix', 'probabilities', or 'state_vector'") \ No newline at end of file + raise ValueError("output should in be 'density_matrix', 'probabilities', or 'state_vector'") diff --git a/src/quafu/tasks/tasks.py b/src/quafu/tasks/tasks.py index 7221153..e6818f6 100644 --- a/src/quafu/tasks/tasks.py +++ b/src/quafu/tasks/tasks.py @@ -1,44 +1,45 @@ -import os from typing import Dict, List, Tuple - from quafu.circuits.quantum_circuit import QuantumCircuit -from quafu.users.userapi import load_account, get_backends_info from ..exceptions import CircuitError, ServerError, CompileError from ..results.results import ExecResult, merge_measure -from ..backends.backends import Backend from ..users.exceptions import UserError import numpy as np import json import requests from urllib import parse -import re +from quafu.users.userapi import User +from quafu.backends.backends import Backend import copy -import networkx as nx -import matplotlib.pyplot as plt -class Task(object): +class Task(object): """ Class for submitting quantum computation task to the backend. Attributes: - token (str): Apitoken that associate to your Quafu account. shots (int): Numbers of single shot measurement. compile (bool): Whether compile the circuit on the backend tomo (bool): Whether do tomography (Not support yet) + user (User): User object corresponding to Quafu account + priority (int): priority level of the task + submit_history (dict): circuit submitted with this task + backend (dict): quantum backend that execute the task. + """ - def __init__(self): + def __init__(self, user = User()): + self.user = user self.shots = 1000 self.tomo = False self.compile = True - self.priority = 2 + self.priority = self.user.priority + self.runtime_job_id = "" self.submit_history = { } - self.token, self._url = load_account() - self._available_backends = {info["system_name"]:Backend(info) for info in get_backends_info()} + self._available_backends = self.user.get_available_backends(print_info=False) self.backend = self._available_backends[list(self._available_backends.keys())[0]] - def config(self, - backend: str="ScQ-P10", - shots: int=1000, + + def config(self, + backend: str="ScQ-P10", + shots: int=1000, compile: bool=True, tomo: bool=False, priority: int=2) -> None: @@ -61,24 +62,26 @@ def config(self, self.compile = compile self.priority = priority - + def get_history(self) -> Dict: """ Get the history of submitted task. Returns: - A dict of history. The key is the group name and the value is a list of task id in the group. + A dict of history. The key is the group name and the value is a list of task id in the group. """ return self.submit_history + + def get_backend_info(self) -> Dict: """ Get the calibration information of the experimental backend. - Returns: + Returns: Backend information dictionary containing the mapping from the indices to the names of physical bits `'mapping'`, backend topology `'topology_diagram'` and full calibration inforamtion `'full_info'`. """ - return self.backend.get_chip_info() - + return self.backend.get_chip_info(self.user) + def submit(self, qc: QuantumCircuit, obslist: List=[])\ @@ -89,15 +92,15 @@ def submit(self, qc (QuantumCircuit): Quantum circuit that need to be executed on backend. obslist (list[str, list[int]]): List of pauli string and its position. - Returns: + Returns: List of executed results and list of measured observable - Examples: + Examples: 1) input [["XYX", [0, 1, 2]], ["Z", [1]]] measure pauli operator XYX at 0, 1, 2 qubit, and Z at 1 qubit.\n 2) Measure 5-qubit Ising Hamiltonian we can use\n obslist = [["X", [i]] for i in range(5)]]\n obslist.extend([["ZZ", [i, i+1]] for i in range(4)])\n - + For the energy expectation of Ising Hamiltonian \n res, obsexp = q.submit_task(obslist)\n E = sum(obsexp) @@ -133,8 +136,8 @@ def submit(self, return exec_res, measure_results - def run(self, - qc: QuantumCircuit, + def run(self, + qc: QuantumCircuit, measure_base: List=[]) -> ExecResult: """Single run for measurement task. @@ -158,9 +161,9 @@ def run(self, return res - def send(self, - qc: QuantumCircuit, - name: str="", + def send(self, + qc: QuantumCircuit, + name: str="", group: str="", wait: bool=True) -> ExecResult: """ @@ -171,7 +174,7 @@ def send(self, name: Task name. group: The task belong which group. wait: Whether wait until the execution return. - Returns: + Returns: ExecResult object that contain the dict return from quantum device. """ from quafu import get_version @@ -183,14 +186,15 @@ def send(self, qc.to_openqasm() data = {"qtasm": qc.openqasm, "shots": self.shots, "qubits": qc.num, "scan": 0, "tomo": int(self.tomo), "selected_server": self.backend.system_id, - "compile": int(self.compile), "priority": self.priority, "task_name": name, "pyquafu_version": version} - + "compile": int(self.compile), "priority": self.priority, "task_name": name, + "pyquafu_version": version, "runtime_job_id": self.runtime_job_id} + if wait: - url = self._url + "qbackend/scq_kit/" + url = self.user._url + self.user.exec_api else: - url = self._url + "qbackend/scq_kit_asyc/" - - headers = {'Content-Type': 'application/x-www-form-urlencoded;charset=UTF-8', 'api_token': self.token} + url = self.user._url + self.user.exec_async_api + + headers = {'Content-Type': 'application/x-www-form-urlencoded;charset=UTF-8', 'api_token': self.user.apitoken} data = parse.urlencode(data) data = data.replace("%27", "'") res = requests.post(url, headers=headers, data=data) @@ -203,10 +207,10 @@ def send(self, elif res.json()["status"] == 5003: raise ServerError(res_dict["message"]) elif res.json()["status"] == 5004: - raise CompileError(res_dict["message"]) + raise CompileError(res_dict["message"]) else: task_id = res_dict["task_id"] - + if not (group in self.submit_history): self.submit_history[group] = [task_id] else: @@ -217,24 +221,24 @@ def send(self, def retrieve(self, taskid: str) -> ExecResult: """ Retrieve the results of submited task by taskid. - + Args: taskid: The taskid of the task need to be retrieved. """ data = {"task_id" : taskid} - url = self._url + "qbackend/scq_task_recall/" + url = self.user._url + self.user.exec_recall_api headers = {'Content-Type': 'application/x-www-form-urlencoded;charset=UTF-8', 'api_token': self.token} res = requests.post(url, headers=headers, data=data) res_dict = json.loads(res.text) measures = eval(res_dict["measure"]) - + return ExecResult(res_dict, measures) - def retrieve_group(self, + def retrieve_group(self, group: str, - history: Dict={}, + history: Dict={}, verbose: bool=True) -> List[ExecResult]: """ Retrieve the results of submited task by group name. @@ -244,7 +248,7 @@ def retrieve_group(self, history: History from which to retrieve the results. If not provided, the history will be the submit history of saved by current task. verbose: Whether print the task status in the group. Returns: - A list of execution results in the retrieved group. Only completed task will be added. + A list of execution results in the retrieved group. Only completed task will be added. """ history = history if history else self.submit_history taskids = history[group] @@ -264,8 +268,8 @@ def retrieve_group(self, group_res.append(res) return group_res - - + + def check_valid_gates(self, qc: QuantumCircuit) -> None: """ Check the validity of the quantum circuit. @@ -277,7 +281,7 @@ def check_valid_gates(self, qc: QuantumCircuit) -> None: for gate in qc.gates: if gate.name.lower() not in valid_gates: raise CircuitError("Invalid operations '%s' for backend '%s'" %(gate.name, self.backend.name)) - + else: if self.backend.name == "ScQ-S41": raise CircuitError("Backend ScQ-S41 must be used without compilation") @@ -285,4 +289,3 @@ def check_valid_gates(self, qc: QuantumCircuit) -> None: for gate in qc.gates: if gate.name.lower() in ["xy"]: raise CircuitError("Invalid operations '%s' for backend '%s'" %(gate.name, self.backend.name)) - diff --git a/src/quafu/users/userapi.py b/src/quafu/users/userapi.py index 071cfae..e51fa1a 100644 --- a/src/quafu/users/userapi.py +++ b/src/quafu/users/userapi.py @@ -5,10 +5,21 @@ from urllib import parse from .exceptions import UserError + class User(object): def __init__(self): self.apitoken = "" - + self._url = "http://quafu.baqis.ac.cn/" + self.load_account() + + self.backends_api = "qbackend/get_backends/" + self.chip_api = "qbackend/scq_get_chip_info/" + self.exec_api = "qbackend/scq_kit/" + self.exec_async_api = "qbackend/scq_kit_asyc/" + self.exec_recall_api = "qbackend/scq_task_recall/" + self.priority = 2 + + def save_apitoken(self, apitoken): """ Save your apitoken associate your Quafu account. @@ -21,27 +32,44 @@ def save_apitoken(self, apitoken): with open(file_dir + "api", "w") as f: f.write(self.apitoken+"\n") f.write("http://quafu.baqis.ac.cn/") - -def load_account(): - """ - Load Quafu account. - """ - homedir = get_homedir() - file_dir = homedir + "/.quafu/" - try: - f = open(file_dir + "api", "r") - data = f.readlines() - token = data[0].strip("\n") - url = data[1].strip("\n") - return token, url - except: - raise UserError("User configure error. Please set up your token.") - -def get_backends_info(): - """ - Get available backends information - """ - token, _url = load_account() - backends_info = requests.post(url=_url+"qbackend/get_backends/", headers={"api_token" : token}) - backends_dict = json.loads(backends_info.text) - return backends_dict["data"] + + def load_account(self): + """ + Load Quafu account. + """ + homedir = get_homedir() + file_dir = homedir + "/.quafu/" + try: + f = open(file_dir + "api", "r") + data = f.readlines() + token = data[0].strip("\n") + url = data[1].strip("\n") + self.apitoken = token + self._url = url + return token, url + except: + raise UserError("User configure error. Please set up your token.") + + def get_backends_info(self): + """ + Get available backends information + """ + + backends_info = requests.post(url=self._url+self.backends_api, headers={"api_token" : self.apitoken}) + backends_info_dict = json.loads(backends_info.text) + if backends_info_dict["status"] == 201: + raise UserError(backends_info_dict["message"]) + else: + return backends_info_dict["data"] + + def get_available_backends(self, print_info=True): + from quafu.backends.backends import Backend + backends_info = self.get_backends_info() + self._available_backends = {info["system_name"]:Backend(info) for info in backends_info} + + if print_info: + print((" "*5).join(["system_name".ljust(10), "qubits".ljust(10), "status".ljust(10)])) + for backend in self._available_backends.values(): + print((" "*5).join([backend.name.ljust(10), str(backend.qubit_num).ljust(10), backend.status.ljust(10)])) + + return self._available_backends diff --git a/src/quafu/visualisation/__init__.py b/src/quafu/visualisation/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/quafu/visualisation/circuitPlot.py b/src/quafu/visualisation/circuitPlot.py new file mode 100644 index 0000000..9893cae --- /dev/null +++ b/src/quafu/visualisation/circuitPlot.py @@ -0,0 +1,528 @@ +import matplotlib.patheffects as pe +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.collections import PolyCollection, PatchCollection, LineCollection +from matplotlib.patches import Circle, Arc +from matplotlib.text import Text + +from quafu.elements.quantum_element import Instruction, ControlledGate + +# the following line for developers only +# from quafu.circuits.quantum_circuit import QuantumCircuit + +line_args = {} +box_args = {} + +DEEPCOLOR = '#0C161F' +BLUE = '#1f77b4' +ORANGE = '#ff7f0e' +GREEN = '#2ca02c' +GOLDEN = '#FFB240' +GARNET = '#C0392B' + +""" +layers(zorder): + +0: figure +1: bkg box +2: wires +3: closed patches +4: white bkg for label/text +5: labels +""" + +su2_gate_names = ['x', 'y', 'z', 'id', 'w', + 'h', 't', 'tdg', 's', 'sdg', 'sx', 'sy', 'sw', + 'phase', + 'rx', 'ry', 'rz', + ] + +swap_gate_names = ['swap', 'iswap'] +r2_gate_names = ['rxx', 'ryy', 'rzz'] +c2_gate_names = ['cp', 'cs', 'ct', 'cx', 'cy', 'cz'] +c3_gate_names = ['fredkin', 'toffoli'] +mc_gate_names = ['mcx', 'mcy', 'mcz'] +operation_names = ['barrier', 'delay'] + + +class CircuitPlotManager: + """ + A class to manage the plotting of quantum circuits. + Stores style parameters and provides functions to plot. + + To be initialized when circuit.plot() is called. + """ + # colors + _wire_color = '#FF0000' + _light_blue = '#3B82F6' + _ec = DEEPCOLOR + + _wire_lw = 1.5 + + _a_inch = 2 / 2.54 # physical lattice constant in inch + _a = 0.5 # box width and height, unit: ax + _barrier_width = _a / 3 # barrier width + + _stroke = pe.withStroke(linewidth=2, foreground='white') + + def __init__(self, qc): + """ + Processing graphical info from a quantum circuit, decompose every + gate/instruction into graphical elements and send the latter into + corresponding containers. + + At present quantum gates and barriers are stored as a list. In the + near future the circuit will be stored as a graph or graph-like object, + procedure will be much simplified, and more functions will be developed. + (TODO) + """ + self.qbit_num = qc.num + # step0: containers of graphical elements + + self._h_wire_points = [] + self._ctrl_wire_points = [] + + self._closed_patches = [] + self._mea_arc_patches = [] + self._mea_point_patches = [] + + self._ctrl_points = [] + self._not_points = [] + self._swap_points = [] + self._iswap_points = [] + self._barrier_points = [] + + self._text_list = [] + + # step1: process gates/instructions + self.dorders = np.zeros(qc.num, dtype=int) + for gate in qc.gates: + assert isinstance(gate, Instruction) + self._process_ins(gate) + + self.depth = np.max(self.dorders) + 1 + + for q, c in qc.measures.items(): + self._proc_measure(self.depth - 1, q) + + # step2: initialize bit-label + self.q_label = {i: f'q_{i}' for i in range(qc.num)} + self.c_label = {iq: f'c_{ic}' for iq, ic in qc.measures.items()} + + # step3: figure coordination + self.xs = np.arange(-3 / 2, self.depth + 3 / 2) + self.ys = np.arange(-2, self.qbit_num + 1 / 2) + + def __call__(self, + title=None, + init_labels=None, + end_labels=None, + save_path: str = None, + show: bool = False, + *args, + **kwargs): + """ + :param title + :param init_labels: dict, {qbit: label} + :param end_labels: dict, {qbit: label} + :param save_path: str, path to save the figure + :param show: bool, whether to show the figure + :param args: + :param kwargs: + + More customization will be supported in the future.(TODO) + """ + # Not supported by patch collections? + # if 'xkcd' in kwargs: + # import random + # plt.gca().xkcd(randomness=random.randint(0, 1000)) + if title is not None: + title = Text((self.xs[0] + self.xs[-1]) / 2, -0.8, + title, + size=30, + ha='center', va='baseline') + self._text_list.append(title) + + # initialize a figure + _size_x = self._a_inch * abs(self.xs[-1] - self.xs[0]) + _size_y = self._a_inch * abs(self.ys[-1] - self.ys[0]) + fig = plt.figure(figsize=(_size_x, _size_y)) # inch + ax = fig.add_axes([0, 0, 1, 1], + aspect=1, + xlim=[self.xs[0], self.xs[-1]], + ylim=[self.ys[0], self.ys[-1]], + ) + ax.axis('off') + ax.invert_yaxis() + + self._circuit_wires() + self._inits_label(labels=init_labels) + self._measured_label(labels=end_labels) + self._render_circuit() + + if save_path is not None: + plt.savefig(save_path, dpi=300, bbox_inches='tight') + if show: + plt.show() + + def _process_ins(self, ins: Instruction, append: bool = True): + name = ins.name + assert name in Instruction.ins_classes, 'If this should occur, please report a bug.' + + name = name.lower() + _which = slice(np.min(ins.pos), np.max(ins.pos) + 1) + depth = np.max(self.dorders[_which]) + paras = ins.paras + + if name == 'barrier': + self._proc_barrier(depth, ins.pos) + elif name == 'measure': + self._proc_measure(depth, ins.pos) + elif name in su2_gate_names: + self._proc_su2(name, depth, ins.pos, paras) + elif name in swap_gate_names: + self._proc_swap(depth, ins.pos, name == 'iswap') + elif name in r2_gate_names: + # TODO: combine into one box + self._proc_su2(name[-1], depth, ins.pos[0], paras) + self._proc_su2(name[-1], depth, ins.pos[1], paras) + elif isinstance(ins, ControlledGate): + self._proc_ctrl(depth, ins) + else: + raise NotImplementedError(f'Gate {name} is not supported yet.\n' + f'If this should occur, please report a bug.') + if append: + self.dorders[_which] = depth + 1 + + ######################################################################### + # Helper functions for processing gates/instructions into graphical + # elements. Add only points data of for the following collection-wise + # plotting if possible, create a patch otherwise. + ######################################################################### + def _circuit_wires(self): + """ + plot horizontal circuit wires + """ + for y in range(self.qbit_num): + x0 = self.xs[0] + 1 + x1 = self.xs[-1] - 1 + self._h_wire_points.append([[x0, y], [x1, y]]) + + def _gate_bbox(self, x, y, fc: str): + a = self._a + from matplotlib.patches import FancyBboxPatch + bbox = FancyBboxPatch((-a / 2 + x, -a / 2 + y), a, a, # this warning belongs to matplotlib + boxstyle=f'round, pad={0.2 * a}', + edgecolor=DEEPCOLOR, + facecolor=fc, + ) + self._closed_patches.append(bbox) + + def _inits_label(self, labels: dict[int: str] = None): + """ qubit-labeling """ + if labels is None: + labels = self.q_label + + for i, label in labels.items(): + label = r'$|%s\rangle$' % label + txt = Text(-2 / 3, i, + label, + size=18, + color=DEEPCOLOR, + ha='right', + va='center', + ) + self._text_list.append(txt) + + def _measured_label(self, labels: dict[int: str] = None): + """ measured qubit-labeling """ + if labels is None: + labels = self.c_label + + for i, label in labels.items(): + label = r'$%s$' % label + txt = Text(self.xs[-1] - 3 / 4, i, + label, + size=18, + color=DEEPCOLOR, + ha='left', + va='center', + ) + self._text_list.append(txt) + + def _gate_label(self, s, x, y): + if not s: + return None + _dy = 0.05 + text = Text(x, y + _dy, + s, + size=24, + color=DEEPCOLOR, + ha='center', + va='center', + ) + text.set_path_effects([self._stroke]) + self._text_list.append(text) + + def _para_label(self, para_txt, x, y): + """ label parameters """ + if not para_txt: + return None + _dx = 0 + text = Text(x + _dx, y + 0.8 * self._a, + para_txt, + size=12, + color=DEEPCOLOR, + ha='center', + va='top', + ) + self._text_list.append(text) + + def _measure_label(self, x, y): + from matplotlib.patches import FancyArrow + a = self._a + r = 1.1 * a + d = 1.2 * a / 3.5 + + arrow = FancyArrow(x=x, + y=y + d, + dx=0.15, + dy=-0.35, + width=0.04, + facecolor=DEEPCOLOR, + head_width=0.07, + head_length=0.15, + edgecolor='white') + arc = Arc((x, y + d), + width=r, + height=r, + lw=1, + theta1=180, + theta2=0, + fill=False, + zorder=4, + color=DEEPCOLOR, + capstyle='round', + ) + center_bkg = Circle((x, y + d), + radius=0.035, + color='white', + ) + center = Circle((x, y + d), + radius=0.025, + facecolor=DEEPCOLOR, + ) + self._mea_arc_patches.append(arc) + self._mea_point_patches += [center_bkg, arrow, center] + + ######################################################################### + # # # # processing-functions: decompose ins into graphical elements # # # + ######################################################################### + def _proc_su2(self, id_name, depth, pos, paras): + if id_name in ['x', 'y', 'z', 'h', 'id', 's', 't', 'p', 'u']: + fc = '#EE7057' + label = id_name.capitalize() + elif id_name in ['rx', 'ry', 'rz']: + fc = '#6366F1' + label = id_name.upper() + else: + fc = '#8C9197' + label = '?' + + if id_name in ['rx', 'ry', 'rz', 'p']: + # too long to display: r'$\theta=$' + f'{paras:.3f}' (TODO) + para_txt = f'({paras:.3f})' + else: + para_txt = None + + self._gate_label(label, depth, pos) + self._para_label(para_txt, depth, pos) + self._gate_bbox(depth, pos, fc) + + def _proc_ctrl(self, depth, ins: ControlledGate, ctrl_type: bool = True): + # control part + p0, p1 = np.max(ins.pos), np.min(ins.pos) + self._ctrl_wire_points.append([[depth, p1], [depth, p0]]) + + ctrl_pos = np.array(ins.ctrls) + for c in ctrl_pos: + self._ctrl_points.append((depth, c, ctrl_type)) + + # target part + name = ins.name.lower() + if ins.ct_nums == (1, 1, 2) or name in mc_gate_names: + tar_name = ins.targ_name.lower()[-1] + pos = ins.targs if isinstance(ins.targs, int) else ins.targs[0] + if tar_name == 'x': + self._not_points.append((depth, pos)) + else: + self._proc_su2(tar_name, depth, pos, None) + elif name == 'cswap': + self._swap_points += [[depth, p] for p in ins.targs] + elif name == 'ccx': + self._not_points.append((depth, ins.targs)) + else: + from quafu.elements.element_gates import ControlledU + assert isinstance(ins, ControlledU), f'unknown gate: {name}, {ins.__class__.__name__}' + self._process_ins(ins, append=False) + + def _proc_swap(self, depth, pos, iswap: bool = False): + p1, p2 = pos + nodes = [[depth, p] for p in pos] + self._swap_points += nodes + self._ctrl_wire_points.append([[depth, p1], [depth, p2]]) + if iswap: + self._iswap_points += nodes + + def _proc_barrier(self, depth, pos: list): + x0 = depth - self._barrier_width + x1 = depth + self._barrier_width + + for p in pos: + y0 = (p - 1 / 2) + y1 = (p + 1 / 2) + nodes = [[x0, y0], [x0, y1], [x1, y1], [x1, y0], [x0, y0]] + self._barrier_points.append(nodes) + + def _proc_measure(self, depth, pos): + fc = GOLDEN + self._gate_bbox(depth, pos, fc) + self._measure_label(depth, pos) + + # TODO: decide whether to draw double wire for measurement + # y = pos + 0.02 + # x0 = depth + # x1 = self.depth - 1 / 2 + # self._h_wire_points.append([[x0, y], [x1, y]]) + + ######################################################################### + # # # # # # # # # # # # # # rendering functions # # # # # # # # # # # # # + ######################################################################### + def _render_h_wires(self): + h_lines = LineCollection(self._h_wire_points, + zorder=0, + colors=self._wire_color, + alpha=0.8, + linewidths=2, + ) + plt.gca().add_collection(h_lines) + + def _render_ctrl_wires(self): + v_lines = LineCollection(self._ctrl_wire_points, + zorder=0, + colors=self._light_blue, + alpha=0.8, + linewidths=4, + ) + plt.gca().add_collection(v_lines) + + def _render_closed_patch(self): + collection = PatchCollection(self._closed_patches, + match_original=True, + zorder=3, + ec=self._ec, + linewidths=0.5, + ) + plt.gca().add_collection(collection) + + def _render_ctrl_nodes(self): + circle_collection = [] + r = self._a / 4 + for x, y, ctrl in self._ctrl_points: + fc = '#3B82F6' if ctrl else 'white' + circle = Circle((x, y), radius=r, fc=fc) + circle_collection.append(circle) + circles = PatchCollection(circle_collection, + match_original=True, + zorder=5, + ec=self._ec, + linewidths=2, + ) + plt.gca().add_collection(circles) + + def _render_not_nodes(self): + points = [] + rp = self._a * 0.3 + r = self._a * 0.5 + + for x, y in self._not_points: + points.append([[x, y - rp], [x, y + rp]]) + points.append([[x - rp, y], [x + rp, y]]) + circle = Circle((x, y), radius=r, lw=1, + fc='#3B82F6') + self._closed_patches.append(circle) + + collection = LineCollection(points, + zorder=5, + colors='white', + linewidths=2, + capstyle='round', + ) + plt.gca().add_collection(collection) + + def _render_swap_nodes(self): + points = [] + r = self._a / (4 ** (1 / 2)) + for x, y in self._swap_points: + points.append([[x - r, y - r], [x + r, y + r]]) + points.append([[x + r, y - r], [x - r, y + r]]) + collection = LineCollection(points, + zorder=5, + colors='#3B82F6', + linewidths=4, + capstyle='round', + ) + plt.gca().add_collection(collection) + + # iswap-cirlces + i_circles = [] + for x, y in self._iswap_points: + circle = Circle((x, y), radius=2 ** (1 / 2) * r, lw=3, + ec='#3B82F6', fill=False) + i_circles.append(circle) + collection = PatchCollection(i_circles, + match_original=True, + zorder=5, + ) + plt.gca().add_collection(collection) + + def _render_measure(self): + stroke = pe.withStroke(linewidth=4, foreground='white') + arcs = PatchCollection(self._mea_arc_patches, + match_original=True, + capstyle='round', + zorder=4) + arcs.set_path_effects([stroke]) + + plt.gca().add_collection(arcs) + pointers = PatchCollection(self._mea_point_patches, # note the order + match_original=True, + zorder=5, + facecolors=DEEPCOLOR, + linewidths=2, + ) + plt.gca().add_collection(pointers) + + def _render_barrier(self): + barrier = PolyCollection(self._barrier_points, + closed=True, + fc='lightgray', + hatch='///', + zorder=4) + plt.gca().add_collection(barrier) + + def _render_txt(self): + for txt in self._text_list: + plt.gca().add_artist(txt) + + def _render_circuit(self): + self._render_h_wires() + self._render_ctrl_wires() + self._render_ctrl_nodes() + self._render_not_nodes() + + self._render_swap_nodes() + self._render_measure() + self._render_barrier() + self._render_closed_patch() + self._render_txt()