From eecb99d8759471b841e7904561cd99ad93c7a402 Mon Sep 17 00:00:00 2001 From: KowerKoint Date: Wed, 20 Nov 2024 12:23:49 +0900 Subject: [PATCH] merge_gate for new cmake --- include/scaluq/all.hpp | 1 + include/scaluq/gate/merge_gate.hpp | 18 ++ include/scaluq/util/random.hpp | 11 + include/scaluq/util/utility.hpp | 26 ++- src/CMakeLists.txt | 1 + src/gate/gate_standard.cpp | 5 +- src/gate/merge_gate.cpp | 248 ++++++++++++++++++++++ tests/CMakeLists.txt | 2 +- tests/gate/gate_test.cpp | 14 +- tests/gate/merge_test.cpp | 170 +++++++++++++++ tests/gate/param_gate_test.cpp | 17 +- tests/operator/test_operator.cpp | 18 +- tests/state/state_vector_batched_test.cpp | 8 +- tests/state/state_vector_test.cpp | 52 ++--- tests/util/util.hpp | 30 +-- 15 files changed, 537 insertions(+), 84 deletions(-) create mode 100644 include/scaluq/gate/merge_gate.hpp create mode 100644 src/gate/merge_gate.cpp create mode 100644 tests/gate/merge_test.cpp diff --git a/include/scaluq/all.hpp b/include/scaluq/all.hpp index ebe94107..3d58f9a4 100644 --- a/include/scaluq/all.hpp +++ b/include/scaluq/all.hpp @@ -4,6 +4,7 @@ #include "constant.hpp" #include "gate/gate.hpp" #include "gate/gate_factory.hpp" +#include "gate/merge_gate.hpp" #include "gate/param_gate.hpp" #include "gate/param_gate_factory.hpp" #include "kokkos.hpp" diff --git a/include/scaluq/gate/merge_gate.hpp b/include/scaluq/gate/merge_gate.hpp new file mode 100644 index 00000000..623c7a21 --- /dev/null +++ b/include/scaluq/gate/merge_gate.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include "gate.hpp" + +namespace scaluq { +template +std::pair, Fp> merge_gate(const Gate& gate1, const Gate& gate2); + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_merge_gate_hpp(nb::module_& m) { + m.def("merge_gate", + &merge_gate, + "Merge two gates. return value is (merged gate, global phase)."); +} +} // namespace internal +#endif +} // namespace scaluq diff --git a/include/scaluq/util/random.hpp b/include/scaluq/util/random.hpp index 5d4ca49e..5b8015df 100644 --- a/include/scaluq/util/random.hpp +++ b/include/scaluq/util/random.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../types.hpp" @@ -21,5 +22,15 @@ class Random { [[nodiscard]] std::uint64_t int64() { return this->mt(); } [[nodiscard]] std::uint32_t int32() { return this->mt() % UINT32_MAX; } + + [[nodiscard]] std::vector permutation(std::uint64_t n) { + std::vector shuffled(n); + std::iota(shuffled.begin(), shuffled.end(), 0ULL); + for (std::uint64_t i : std::views::iota(0ULL, n) | std::views::reverse) { + std::uint64_t j = int32() % (i + 1); + if (i != j) std::swap(shuffled[i], shuffled[j]); + } + return shuffled; + } }; } // namespace scaluq diff --git a/include/scaluq/util/utility.hpp b/include/scaluq/util/utility.hpp index 62bae30f..aa57d74f 100644 --- a/include/scaluq/util/utility.hpp +++ b/include/scaluq/util/utility.hpp @@ -92,10 +92,11 @@ template inline internal::ComplexMatrix get_expanded_matrix( const internal::ComplexMatrix& from_matrix, const std::vector& from_targets, - std::vector& to_targets) { + std::uint64_t from_control_mask, + std::vector& to_operands) { std::vector targets_map(from_targets.size()); std::ranges::transform(from_targets, targets_map.begin(), [&](std::uint64_t x) { - return std::ranges::lower_bound(to_targets, x) - to_targets.begin(); + return std::ranges::lower_bound(to_operands, x) - to_operands.begin(); }); std::vector idx_map(1ULL << from_targets.size()); for (std::uint64_t i : std::views::iota(0ULL, 1ULL << from_targets.size())) { @@ -103,22 +104,33 @@ inline internal::ComplexMatrix get_expanded_matrix( idx_map[i] |= (i >> j & 1) << targets_map[j]; } } + std::uint64_t to_control_mask = 0; + for (std::uint64_t sub_mask = from_control_mask; sub_mask; sub_mask &= (sub_mask - 1)) { + to_control_mask |= + 1ULL << (std::ranges::lower_bound(to_operands, std::countr_zero(sub_mask)) - + to_operands.begin()); + } std::uint64_t targets_idx_mask = idx_map.back(); std::vector outer_indices; - outer_indices.reserve(1ULL << (to_targets.size() - from_targets.size())); - for (std::uint64_t i : std::views::iota(0ULL, 1ULL << to_targets.size())) { - if ((i & targets_idx_mask) == 0) outer_indices.push_back(i); + outer_indices.reserve( + 1ULL << (to_operands.size() - from_targets.size() - std::popcount(from_control_mask))); + for (std::uint64_t i : std::views::iota(0ULL, 1ULL << to_operands.size())) { + if ((i & (targets_idx_mask | to_control_mask)) == 0) outer_indices.push_back(i); } internal::ComplexMatrix to_matrix = - internal::ComplexMatrix::Zero(1ULL << to_targets.size(), 1ULL << to_targets.size()); + internal::ComplexMatrix::Zero(1ULL << to_operands.size(), 1ULL << to_operands.size()); for (std::uint64_t i : std::views::iota(0ULL, 1ULL << from_targets.size())) { for (std::uint64_t j : std::views::iota(0ULL, 1ULL << from_targets.size())) { for (std::uint64_t o : outer_indices) { - to_matrix(idx_map[i] | o, idx_map[j] | o) = from_matrix(i, j); + to_matrix(idx_map[i] | to_control_mask | o, idx_map[j] | to_control_mask | o) = + from_matrix(i, j); } } } + for (std::uint64_t i : std::views::iota(0ULL, 1ULL << to_operands.size())) { + if ((i & to_control_mask) != to_control_mask) to_matrix(i, i) = 1; + } return to_matrix; } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d6bf7497..93c08da7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.21) target_sources(scaluq PRIVATE circuit/circuit.cpp + gate/merge_gate.cpp gate/gate_matrix.cpp gate/gate_pauli.cpp gate/gate_probablistic.cpp diff --git a/src/gate/gate_standard.cpp b/src/gate/gate_standard.cpp index 5dcd9ecd..3a1b9aad 100644 --- a/src/gate/gate_standard.cpp +++ b/src/gate/gate_standard.cpp @@ -267,7 +267,8 @@ FLOAT_DECLARE_CLASS(SqrtYGateImpl) FLOAT(Fp) ComplexMatrix SqrtYdagGateImpl::get_matrix() const { internal::ComplexMatrix mat(2, 2); - mat << 0, StdComplex(0, -1), StdComplex(0, 1), 0; + mat << StdComplex(0.5, -0.5), StdComplex(0.5, -0.5), StdComplex(-0.5, 0.5), + StdComplex(0.5, -0.5); return mat; } FLOAT(Fp) @@ -328,7 +329,7 @@ FLOAT(Fp) ComplexMatrix RXGateImpl::get_matrix() const { internal::ComplexMatrix mat(2, 2); mat << std::cos(this->_angle / 2), StdComplex(0, -std::sin(this->_angle / 2)), - StdComplex(0, std::sin(this->_angle / 2)), std::cos(this->_angle / 2); + StdComplex(0, -std::sin(this->_angle / 2)), std::cos(this->_angle / 2); return mat; } FLOAT(Fp) diff --git a/src/gate/merge_gate.cpp b/src/gate/merge_gate.cpp new file mode 100644 index 00000000..d273b6d2 --- /dev/null +++ b/src/gate/merge_gate.cpp @@ -0,0 +1,248 @@ +#include +#include +#include +#include + +#include "../util/template.hpp" + +namespace scaluq { +FLOAT(Fp) +std::pair, Fp> merge_gate_dense_matrix(const Gate& gate1, const Gate& gate2) { + auto common_control_mask = gate1->control_qubit_mask() & gate2->control_qubit_mask(); + auto merged_operand_mask = + (gate1->operand_qubit_mask() | gate2->operand_qubit_mask()) & ~common_control_mask; + auto merged_operand_vector = internal::mask_to_vector(merged_operand_mask); + auto matrix1 = internal::get_expanded_matrix(gate1->get_matrix(), + gate1->target_qubit_list(), + gate1->control_qubit_mask() & ~common_control_mask, + merged_operand_vector); + auto matrix2 = internal::get_expanded_matrix(gate2->get_matrix(), + gate2->target_qubit_list(), + gate2->control_qubit_mask() & ~common_control_mask, + merged_operand_vector); + std::cerr << matrix1 << std::endl; + std::cerr << matrix2 << std::endl; + auto matrix = matrix2 * matrix1; + std::cerr << matrix << std::endl; + return {gate::DenseMatrix( + merged_operand_vector, matrix, internal::mask_to_vector(common_control_mask)), + 0.}; +} + +FLOAT(Fp) +std::pair, Fp> merge_gate(const Gate& gate1, const Gate& gate2) { + GateType gate_type1 = gate1.gate_type(); + GateType gate_type2 = gate2.gate_type(); + + if (gate_type1 == GateType::Probablistic || gate_type2 == GateType::Probablistic) { + throw std::runtime_error( + "merge_gate(const Gate&, const Gate&): ProbablisticGate is not supported."); + } + + if (gate_type1 == GateType::I) return {gate2, 0.}; + if (gate_type2 == GateType::I) return {gate1, 0.}; + + auto gate1_control_mask = gate1->control_qubit_mask(); + auto gate2_control_mask = gate2->control_qubit_mask(); + + if (gate_type1 == GateType::GlobalPhase && gate1_control_mask == 0) + return {gate2, GlobalPhaseGate(gate1)->phase()}; + if (gate_type2 == GateType::GlobalPhase && gate2_control_mask == 0) + return {gate1, GlobalPhaseGate(gate2)->phase()}; + + if (gate1_control_mask != gate2_control_mask) return merge_gate_dense_matrix(gate1, gate2); + auto control_list = internal::mask_to_vector(gate1_control_mask); + + // Special case: Zero qubit + if (gate_type1 == GateType::GlobalPhase && gate_type2 == GateType::GlobalPhase) { + return {gate::GlobalPhase( + GlobalPhaseGate(gate1)->phase() + GlobalPhaseGate(gate2)->phase(), + control_list), + 0.}; + } + + // Special case: Pauli + auto get_pauli_id = [&](GateType gate_type) -> std::optional { + if (gate_type == GateType::I) return 0; + if (gate_type == GateType::X) return 1; + if (gate_type == GateType::Y) return 2; + if (gate_type == GateType::Z) return 3; + return std::nullopt; + }; + auto pauli_id1 = get_pauli_id(gate_type1); + auto pauli_id2 = get_pauli_id(gate_type2); + assert(!pauli_id1 || pauli_id1 != 0); + assert(!pauli_id2 || pauli_id2 != 0); + if (pauli_id1 && pauli_id2) { + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; + if (target1 == target2) { + if (pauli_id1 == pauli_id2) return {gate::I(), 0.}; + if (pauli_id1 == 1) { + if (pauli_id2 == 2) { + if (gate1_control_mask == 0) { + return {gate::Z(target1, control_list), -Kokkos::numbers::pi / 2}; + } + } + if (pauli_id2 == 3) { + if (gate1_control_mask == 0) { + return {gate::Y(target1, control_list), Kokkos::numbers::pi / 2}; + } + } + } + if (pauli_id1 == 2) { + if (pauli_id2 == 3) { + if (gate1_control_mask == 0) { + return {gate::X(target1, control_list), -Kokkos::numbers::pi / 2}; + } + } + if (pauli_id2 == 1) { + if (gate1_control_mask == 0) { + return {gate::Z(target1, control_list), Kokkos::numbers::pi / 2}; + } + } + } + if (pauli_id1 == 3) { + if (pauli_id2 == 1) { + if (gate1_control_mask == 0) { + return {gate::Y(target1, control_list), -Kokkos::numbers::pi / 2}; + } + } + if (pauli_id2 == 2) { + if (gate1_control_mask == 0) { + return {gate::X(target1, control_list), Kokkos::numbers::pi / 2}; + } + } + } + } + } + if ((pauli_id1 || gate1.gate_type() == GateType::Pauli) && + (pauli_id2 || gate2.gate_type() == GateType::Pauli)) { + auto pauli1 = gate_type1 == GateType::Pauli + ? PauliGate(gate1)->pauli() + : PauliOperator(std::vector{gate1->target_qubit_list()[0]}, + std::vector{pauli_id1.value()}); + auto pauli2 = gate_type2 == GateType::Pauli + ? PauliGate(gate2)->pauli() + : PauliOperator(std::vector{gate2->target_qubit_list()[0]}, + std::vector{pauli_id2.value()}); + return {gate::Pauli(pauli2 * pauli1, control_list), 0.}; + } + + constexpr Fp eps = 1e-12; + + // Special case: Phase + auto get_oct_phase = [&](GateType gate_type) -> std::optional { + if (gate_type == GateType::I) return 0; + if (gate_type == GateType::Z) return 4; + if (gate_type == GateType::S) return 2; + if (gate_type == GateType::Sdag) return 6; + if (gate_type == GateType::T) return 1; + if (gate_type == GateType::Tdag) return 7; + return std::nullopt; + }; + auto oct_phase_gate = [&](std::uint64_t oct_phase, + std::uint64_t target) -> std::optional> { + oct_phase &= 7; + if (oct_phase == 0) return gate::I(); + if (oct_phase == 4) return gate::Z(target, control_list); + if (oct_phase == 2) return gate::S(target, control_list); + if (oct_phase == 6) return gate::Sdag(target, control_list); + if (oct_phase == 1) return gate::T(target, control_list); + if (oct_phase == 7) return gate::Tdag(target, control_list); + return std::nullopt; + }; + auto oct_phase1 = get_oct_phase(gate_type1); + auto oct_phase2 = get_oct_phase(gate_type2); + if (oct_phase1 && oct_phase2) { + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; + if (target1 == target2) { + auto g = oct_phase_gate(oct_phase1.value() + oct_phase2.value(), target1); + if (g) return {g.value(), 0.}; + } + } + if ((oct_phase1 || gate_type1 == GateType::RZ || gate_type1 == GateType::U1) && + (oct_phase2 || gate_type2 == GateType::RZ || gate_type2 == GateType::U1)) { + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; + if (target1 == target2) { + Fp phase1 = oct_phase1 ? oct_phase1.value() * Kokkos::numbers::pi / 4 + : gate_type1 == GateType::RZ ? RZGate(gate1)->angle() + : U1Gate(gate1)->lambda(); + Fp global_phase1 = gate_type1 == GateType::RZ ? -RZGate(gate1)->angle() / 2 : 0.; + Fp phase2 = oct_phase2 ? oct_phase2.value() * Kokkos::numbers::pi / 4 + : gate_type2 == GateType::RZ ? RZGate(gate2)->angle() + : U1Gate(gate2)->lambda(); + Fp global_phase2 = gate_type2 == GateType::RZ ? -RZGate(gate2)->angle() / 2 : 0.; + Fp global_phase = global_phase1 + global_phase2; + if (std::abs(global_phase) < eps) { + return {gate::U1(target1, phase1 + phase2, control_list), + global_phase1 + global_phase2}; + } + } + } + + // Special case: RX + auto get_rx_angle = [&](Gate gate, GateType gate_type) -> std::optional { + if (gate_type == GateType::I) return 0.; + if (gate_type == GateType::X) return Kokkos::numbers::pi; + if (gate_type == GateType::SqrtX) return Kokkos::numbers::pi / 2; + if (gate_type == GateType::SqrtXdag) return -Kokkos::numbers::pi / 2; + if (gate_type == GateType::RX) return RXGate(gate)->angle(); + return std::nullopt; + }; + auto rx_param1 = get_rx_angle(gate1, gate_type1); + auto rx_param2 = get_rx_angle(gate2, gate_type2); + if (rx_param1 && rx_param2) { + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; + Fp global_phase1 = gate_type1 == GateType::RX ? 0. : rx_param1.value() / 2; + Fp global_phase2 = gate_type2 == GateType::RX ? 0. : rx_param2.value() / 2; + Fp global_phase = global_phase1 + global_phase2; + if (target1 == target2) { + if (std::abs(global_phase) < eps) { + return {gate::RX(target1, rx_param1.value() + rx_param2.value(), control_list), + global_phase1 + global_phase2}; + } + } + } + + // Special case: RY + auto get_ry_angle = [&](Gate gate, GateType gate_type) -> std::optional { + if (gate_type == GateType::I) return 0.; + if (gate_type == GateType::Y) return Kokkos::numbers::pi; + if (gate_type == GateType::SqrtY) return Kokkos::numbers::pi / 2; + if (gate_type == GateType::SqrtYdag) return -Kokkos::numbers::pi / 2; + if (gate_type == GateType::RY) return RYGate(gate)->angle(); + return std::nullopt; + }; + auto ry_param1 = get_ry_angle(gate1, gate_type1); + auto ry_param2 = get_ry_angle(gate2, gate_type2); + if (ry_param1 && ry_param2) { + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; + Fp global_phase1 = gate_type1 == GateType::RY ? 0. : ry_param1.value() / 2; + Fp global_phase2 = gate_type2 == GateType::RY ? 0. : ry_param2.value() / 2; + Fp global_phase = global_phase1 + global_phase2; + if (target1 == target2) { + if (std::abs(global_phase) < eps) { + return {gate::RY(target1, ry_param1.value() + ry_param2.value(), control_list), + global_phase1 + global_phase2}; + } + } + } + + // Special case: Swap duplication + if (gate_type1 == gate_type2 && gate_type1 == GateType::Swap) { + if (gate1->target_qubit_mask() == gate2->target_qubit_mask()) return {gate::I(), 0.}; + } + + // General case + return merge_gate_dense_matrix(gate1, gate2); +} +#define FUNC_MACRO(Fp) \ + template std::pair, Fp> merge_gate(const Gate&, const Gate&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO +} // namespace scaluq diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b614da95..fad4890d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,7 +6,7 @@ add_executable(scaluq_test EXCLUDE_FROM_ALL circuit/circuit_test.cpp circuit/param_circuit_test.cpp gate/gate_test.cpp - # gate/merge_test.cpp + gate/merge_test.cpp gate/param_gate_test.cpp operator/test_pauli_operator.cpp operator/test_operator.cpp diff --git a/tests/gate/gate_test.cpp b/tests/gate/gate_test.cpp index 8a9b3d4b..088e6d6e 100644 --- a/tests/gate/gate_test.cpp +++ b/tests/gate/gate_test.cpp @@ -753,12 +753,7 @@ template void test_standard_gate_control(Factory factory, std::uint64_t n) { Random random; - std::vector shuffled(n); - std::iota(shuffled.begin(), shuffled.end(), 0ULL); - for (std::uint64_t i : std::views::iota(0ULL, n) | std::views::reverse) { - std::uint64_t j = random.int32() % (i + 1); - if (i != j) std::swap(shuffled[i], shuffled[j]); - } + std::vector shuffled = random.permutation(n); std::vector targets(num_target); for (std::uint64_t i : std::views::iota(0ULL, num_target)) { targets[i] = shuffled[i]; @@ -845,12 +840,7 @@ void test_pauli_control(std::uint64_t n) { template void test_matrix_control(std::uint64_t n_qubits) { Random random; - std::vector shuffled(n_qubits); - std::iota(shuffled.begin(), shuffled.end(), 0ULL); - for (std::uint64_t i : std::views::iota(0ULL, n_qubits) | std::views::reverse) { - std::uint64_t j = random.int32() % (i + 1); - if (i != j) std::swap(shuffled[i], shuffled[j]); - } + std::vector shuffled = random.permutation(n_qubits); std::vector targets(num_target); for (std::uint64_t i : std::views::iota(0ULL, num_target)) { targets[i] = shuffled[i]; diff --git a/tests/gate/merge_test.cpp b/tests/gate/merge_test.cpp new file mode 100644 index 00000000..82716734 --- /dev/null +++ b/tests/gate/merge_test.cpp @@ -0,0 +1,170 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "../test_environment.hpp" +#include "../util/util.hpp" + +using namespace scaluq; + +template +void merge_gate_test() { + std::vector> gates; + Random random; + std::uint64_t n = 4; + auto none_target_rotation = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::vector controls(shuffled.begin(), shuffled.begin() + nc); + gates.push_back(fac(random.uniform() * std::numbers::pi * 2, controls)); + } + }; + auto single_target = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::uint64_t target = shuffled[0]; + std::vector controls(shuffled.begin() + 1, shuffled.begin() + 1 + nc); + gates.push_back(fac(target, controls)); + } + }; + auto single_target_rotation = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::uint64_t target = shuffled[0]; + std::vector controls(shuffled.begin() + 1, shuffled.begin() + 1 + nc); + gates.push_back(fac(target, random.uniform() * std::numbers::pi * 2, controls)); + } + }; + auto single_target_rotation2 = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::uint64_t target = shuffled[0]; + std::vector controls(shuffled.begin() + 1, shuffled.begin() + 1 + nc); + gates.push_back(fac(target, + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2, + controls)); + } + }; + auto single_target_rotation3 = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::uint64_t target = shuffled[0]; + std::vector controls(shuffled.begin() + 1, shuffled.begin() + 1 + nc); + gates.push_back(fac(target, + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2, + controls)); + } + }; + auto double_target = [&](auto fac) { + for (auto nc : {0, 1, 2}) { + std::vector shuffled = random.permutation(n); + std::uint64_t target0 = shuffled[0]; + std::uint64_t target1 = shuffled[1]; + std::vector controls(shuffled.begin() + 2, shuffled.begin() + 2 + nc); + gates.push_back(fac(target0, target1, controls)); + } + }; + auto dense_matrix = [&](auto fac) { + for (auto nt : {0, 1, 2, 3}) { + for (auto nc : {0, 1, 2}) { + if (nt + nc > static_cast(n)) continue; + std::vector shuffled = random.permutation(n); + std::vector targets(shuffled.begin(), shuffled.begin() + nt); + std::vector controls(shuffled.begin() + nt, + shuffled.begin() + nt + nc); + internal::ComplexMatrix mat(1 << nt, 1 << nt); + for (auto i : std::views::iota(0, 1 << nt)) { + for (auto j : std::views::iota(0, 1 << nt)) { + mat(i, j) = StdComplex(random.uniform() * 2 - 1, random.uniform() * 2 - 1); + } + } + gates.push_back(fac(targets, mat, controls, false)); + } + } + }; + auto sparse_matrix = [&](auto fac) { + for (auto nt : {0, 1, 2, 3}) { + for (auto nc : {0, 1, 2}) { + if (nt + nc > static_cast(n)) continue; + std::vector shuffled = random.permutation(n); + std::vector targets(shuffled.begin(), shuffled.begin() + nt); + std::vector controls(shuffled.begin() + nt, + shuffled.begin() + nt + nc); + internal::SparseComplexMatrix mat(1 << nt, 1 << nt); + for (auto i : std::views::iota(0, 1 << nt)) { + for (auto j : std::views::iota(0, 1 << nt)) { + if (random.uniform() < .5) { + mat.insert(i, j) = + StdComplex(random.uniform() * 2 - 1, random.uniform() * 2 - 1); + } + } + } + gates.push_back(fac(targets, mat, controls)); + } + } + }; + gates.push_back(gate::I()); + none_target_rotation(gate::GlobalPhase); + single_target(gate::X); + single_target(gate::Y); + single_target(gate::Z); + single_target(gate::H); + single_target(gate::S); + single_target(gate::Sdag); + single_target(gate::T); + single_target(gate::Tdag); + single_target(gate::SqrtX); + single_target(gate::SqrtXdag); + single_target(gate::SqrtY); + single_target(gate::SqrtYdag); + single_target(gate::P0); + single_target(gate::P1); + single_target_rotation(gate::RX); + single_target_rotation(gate::RY); + single_target_rotation(gate::RZ); + single_target_rotation(gate::U1); + single_target_rotation2(gate::U2); + single_target_rotation3(gate::U3); + double_target(gate::Swap); + dense_matrix(gate::DenseMatrix); + sparse_matrix(gate::SparseMatrix); + gates.push_back(gate::Pauli(PauliOperator("X 0 Y 2", random.uniform()))); + gates.push_back(gate::Pauli(PauliOperator("Z 0", random.uniform()))); + gates.push_back(gate::Pauli(PauliOperator("Z 3", random.uniform()), {1})); + gates.push_back(gate::Pauli(PauliOperator("Z 1", random.uniform()), {0, 3})); + gates.push_back(gate::PauliRotation(PauliOperator("X 0 Y 2", random.uniform()), + random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::PauliRotation(PauliOperator("Z 0", random.uniform()), + random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::PauliRotation( + PauliOperator("Z 3", random.uniform()), random.uniform() * std::numbers::pi * 2, {1})); + gates.push_back(gate::PauliRotation(PauliOperator("Z 1", random.uniform()), + random.uniform() * std::numbers::pi * 2, + {0, 3})); + for (auto&& g1 : gates) { + for (auto&& g2 : gates) { + auto state1 = StateVector::Haar_random_state(n); + auto state2 = state1.copy(); + auto [mg, phase] = merge_gate(g1, g2); + g1->update_quantum_state(state1); + g2->update_quantum_state(state1); + mg->update_quantum_state(state2); + state2.multiply_coef(Kokkos::polar(static_cast(1.), phase)); + ASSERT_TRUE(same_state(state1, state2)); + } + } +} + +TEST(GateTest, MergeGate) { + merge_gate_test(); + merge_gate_test(); +} diff --git a/tests/gate/param_gate_test.cpp b/tests/gate/param_gate_test.cpp index fa0a7e3a..b53fca86 100644 --- a/tests/gate/param_gate_test.cpp +++ b/tests/gate/param_gate_test.cpp @@ -31,7 +31,7 @@ void test_apply_parametric_single_pauli_rotation(std::uint64_t n_qubits, auto state_cp_amp = state_cp.get_amplitudes(); for (std::uint64_t i = 0; i < dim; i++) { - ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); + ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); } ParamGate pgate_inv = pgate->get_inverse(); @@ -39,7 +39,7 @@ void test_apply_parametric_single_pauli_rotation(std::uint64_t n_qubits, state_amp = state.get_amplitudes(); auto state_bef_amp = state_bef.get_amplitudes(); for (std::uint64_t i = 0; i < dim; i++) { - ASSERT_NEAR(Kokkos::abs(state_bef_amp[i] - state_amp[i]), 0, eps); + ASSERT_NEAR(Kokkos::abs(state_bef_amp[i] - state_amp[i]), 0, eps); } } } @@ -69,7 +69,7 @@ void test_apply_parametric_multi_pauli_rotation(std::uint64_t n_qubits) { auto state_cp_amp = state_cp.get_amplitudes(); // check if the state is updated correctly for (std::uint64_t i = 0; i < dim; i++) { - ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); + ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); } ParamGate pgate_inv = pgate->get_inverse(); pgate_inv->update_quantum_state(state, param); @@ -77,7 +77,7 @@ void test_apply_parametric_multi_pauli_rotation(std::uint64_t n_qubits) { auto state_bef_amp = state_bef.get_amplitudes(); // check if the state is restored correctly for (std::uint64_t i = 0; i < dim; i++) { - ASSERT_NEAR(Kokkos::abs((state_bef_amp[i] - state_amp[i])), 0, eps); + ASSERT_NEAR(Kokkos::abs((state_bef_amp[i] - state_amp[i])), 0, eps); } } } @@ -138,19 +138,14 @@ void test_gate(ParamGate gate_control, amplitudes[internal::insert_zero_at_mask_positions(i, control_mask) | control_mask]), 0., - eps); + eps); } } template void test_param_rotation_control(Factory factory, std::uint64_t n) { Random random; - std::vector shuffled(n); - std::iota(shuffled.begin(), shuffled.end(), 0ULL); - for (std::uint64_t i : std::views::iota(0ULL, n) | std::views::reverse) { - std::uint64_t j = random.int32() % (i + 1); - if (i != j) std::swap(shuffled[i], shuffled[j]); - } + std::vector shuffled = random.permutation(n); std::uint64_t target = shuffled[0]; std::uint64_t num_control = random.int32() % n; std::vector controls(num_control); diff --git a/tests/operator/test_operator.cpp b/tests/operator/test_operator.cpp index 1b1408a9..4954e104 100644 --- a/tests/operator/test_operator.cpp +++ b/tests/operator/test_operator.cpp @@ -60,9 +60,9 @@ TEST(OperatorTest, CheckExpectationValue) { Complex res = rand_observable.get_expectation_value(state); Complex test_res = (test_state.adjoint() * test_rand_observable * test_state)(0, 0); - ASSERT_NEAR(test_res.real(), res.real(), eps); - ASSERT_NEAR(res.imag(), 0, eps); - ASSERT_NEAR(test_res.imag(), 0, eps); + ASSERT_NEAR(test_res.real(), res.real(), eps); + ASSERT_NEAR(res.imag(), 0, eps); + ASSERT_NEAR(test_res.imag(), 0, eps); } } @@ -87,8 +87,8 @@ TEST(OperatorTest, CheckTransitionAmplitude) { Complex res = rand_observable.get_transition_amplitude(state_bra, state_ket); Complex test_res = (test_state_bra.adjoint() * test_rand_observable * test_state_ket)(0, 0); - ASSERT_NEAR(test_res.real(), res.real(), eps); - ASSERT_NEAR(test_res.imag(), res.imag(), eps); + ASSERT_NEAR(test_res.real(), res.real(), eps); + ASSERT_NEAR(test_res.imag(), res.imag(), eps); } } @@ -104,7 +104,7 @@ TEST(OperatorTest, AddTest) { auto exp1 = op1.get_expectation_value(state); auto exp2 = op2.get_expectation_value(state); auto exp = op.get_expectation_value(state); - ASSERT_NEAR(Kokkos::abs(exp1 + exp2 - exp), 0, eps); + ASSERT_NEAR(Kokkos::abs(exp1 + exp2 - exp), 0, eps); } } @@ -120,7 +120,7 @@ TEST(OperatorTest, SubTest) { auto exp1 = op1.get_expectation_value(state); auto exp2 = op2.get_expectation_value(state); auto exp = op.get_expectation_value(state); - ASSERT_NEAR(Kokkos::abs(exp1 - exp2 - exp), 0, eps); + ASSERT_NEAR(Kokkos::abs(exp1 - exp2 - exp), 0, eps); } } @@ -135,7 +135,7 @@ TEST(OperatorTest, MultiCoefTest) { auto state = StateVector::Haar_random_state(n); auto exp1 = op1.get_expectation_value(state); auto exp = op.get_expectation_value(state); - ASSERT_NEAR(Kokkos::abs(exp1 * coef - exp), 0, eps); + ASSERT_NEAR(Kokkos::abs(exp1 * coef - exp), 0, eps); } } @@ -186,6 +186,6 @@ TEST(OperatorTest, Optimize) { ASSERT_EQ(expected.size(), test.size()); for (std::uint64_t i = 0; i < expected.size(); i++) { ASSERT_EQ(expected[i].first, test[i].first); - ASSERT_NEAR(Kokkos::abs(expected[i].second - test[i].second), 0, eps); + ASSERT_NEAR(Kokkos::abs(expected[i].second - test[i].second), 0, eps); } } diff --git a/tests/state/state_vector_batched_test.cpp b/tests/state/state_vector_batched_test.cpp index a7451c49..ed09adfc 100644 --- a/tests/state/state_vector_batched_test.cpp +++ b/tests/state/state_vector_batched_test.cpp @@ -13,7 +13,7 @@ TEST(StateVectorBatchedTest, HaarRandomStateNorm) { const std::uint64_t batch_size = 10, n_qubits = 3; const auto states = StateVectorBatched::Haar_random_state(batch_size, n_qubits, false); auto norms = states.get_squared_norm(); - for (auto x : norms) ASSERT_NEAR(x, 1., eps); + for (auto x : norms) ASSERT_NEAR(x, 1., eps); } TEST(StateVectorBatchedTest, LoadAndAmplitues) { @@ -79,7 +79,7 @@ TEST(StateVectorBatchedTest, ZeroProbs) { auto zero_probs = states.get_zero_probability(i); for (std::uint64_t b = 0; b < batch_size; ++b) { auto state = states.get_state_vector_at(b); - ASSERT_NEAR(zero_probs[b], state.get_zero_probability(i), eps); + ASSERT_NEAR(zero_probs[b], state.get_zero_probability(i), eps); } } } @@ -97,7 +97,7 @@ TEST(StateVectorBatchedTest, MarginalProbs) { auto mg_probs = states.get_marginal_probability(targets); for (std::uint64_t b = 0; b < batch_size; ++b) { auto state = states.get_state_vector_at(b); - ASSERT_NEAR(mg_probs[b], state.get_marginal_probability(targets), eps); + ASSERT_NEAR(mg_probs[b], state.get_marginal_probability(targets), eps); } } } @@ -109,7 +109,7 @@ TEST(StateVectorBatchedTest, Entropy) { auto entropies = states.get_entropy(); for (std::uint64_t b = 0; b < batch_size; ++b) { auto state = states.get_state_vector_at(b); - ASSERT_NEAR(entropies[b], state.get_entropy(), eps); + ASSERT_NEAR(entropies[b], state.get_entropy(), eps); } } diff --git a/tests/state/state_vector_test.cpp b/tests/state/state_vector_test.cpp index 2cd7512d..ef3c0dc2 100644 --- a/tests/state/state_vector_test.cpp +++ b/tests/state/state_vector_test.cpp @@ -14,14 +14,14 @@ TEST(StateVectorTest, HaarRandomStateNorm) { const int n_tries = 20; for (int n = 1; n <= n_tries; n++) { const auto state = StateVector::Haar_random_state(n); - ASSERT_NEAR(state.get_squared_norm(), 1., eps); + ASSERT_NEAR(state.get_squared_norm(), 1., eps); } } { const int n_tries = 20; for (int n = 1; n <= n_tries; n++) { const auto state = StateVector::Haar_random_state(n); - ASSERT_NEAR(state.get_squared_norm(), 1., eps_f); + ASSERT_NEAR(state.get_squared_norm(), 1., eps); } } } @@ -30,8 +30,8 @@ TEST(StateVectorTest, OperationAtIndex) { auto state = StateVector::Haar_random_state(10); for (std::uint64_t i = 0; i < state.dim(); ++i) { state.set_amplitude_at(i, 1); - ASSERT_NEAR(state.get_amplitude_at(i).real(), 1., eps); - ASSERT_NEAR(state.get_amplitude_at(i).imag(), 0., eps); + ASSERT_NEAR(state.get_amplitude_at(i).real(), 1., eps); + ASSERT_NEAR(state.get_amplitude_at(i).imag(), 0., eps); } } @@ -101,8 +101,8 @@ TEST(StateVectorTest, AddState) { for (std::uint64_t i = 0; i < state1.dim(); ++i) { CComplex res = new_vec[i], val = (CComplex)vec1[i] + (CComplex)vec2[i]; - ASSERT_NEAR(res.real(), val.real(), eps); - ASSERT_NEAR(res.imag(), val.imag(), eps); + ASSERT_NEAR(res.real(), val.real(), eps); + ASSERT_NEAR(res.imag(), val.imag(), eps); } } @@ -119,8 +119,8 @@ TEST(StateVectorTest, AddStateWithCoef) { for (std::uint64_t i = 0; i < state1.dim(); ++i) { CComplex res = new_vec[i], val = (CComplex)vec1[i] + coef * (CComplex)vec2[i]; - ASSERT_NEAR(res.real(), val.real(), eps); - ASSERT_NEAR(res.imag(), val.imag(), eps); + ASSERT_NEAR(res.real(), val.real(), eps); + ASSERT_NEAR(res.imag(), val.imag(), eps); } } @@ -135,8 +135,8 @@ TEST(StateVectorTest, MultiplyCoef) { for (std::uint64_t i = 0; i < state.dim(); ++i) { CComplex res = new_vec[i], val = coef * (CComplex)vec[i]; - ASSERT_NEAR(res.real(), val.real(), eps); - ASSERT_NEAR(res.imag(), val.imag(), eps); + ASSERT_NEAR(res.real(), val.real(), eps); + ASSERT_NEAR(res.imag(), val.imag(), eps); } } @@ -150,10 +150,10 @@ TEST(StateVectorTest, GetZeroProbability) { state.add_state_vector_with_coef(std::sqrt(i), tmp_state); } state.normalize(); - ASSERT_NEAR(state.get_zero_probability(0), 30.0 / 55.0, eps); - ASSERT_NEAR(state.get_zero_probability(1), 27.0 / 55.0, eps); - ASSERT_NEAR(state.get_zero_probability(2), 33.0 / 55.0, eps); - ASSERT_NEAR(state.get_zero_probability(3), 28.0 / 55.0, eps); + ASSERT_NEAR(state.get_zero_probability(0), 30.0 / 55.0, eps); + ASSERT_NEAR(state.get_zero_probability(1), 27.0 / 55.0, eps); + ASSERT_NEAR(state.get_zero_probability(2), 33.0 / 55.0, eps); + ASSERT_NEAR(state.get_zero_probability(3), 28.0 / 55.0, eps); } TEST(StateVectorTest, EntropyCalculation) { @@ -165,7 +165,7 @@ TEST(StateVectorTest, EntropyCalculation) { for (std::uint64_t rep = 0; rep < max_repeat; ++rep) { state = StateVector::Haar_random_state(n); auto state_cp = state.get_amplitudes(); - ASSERT_NEAR(state.get_squared_norm(), 1, eps); + ASSERT_NEAR(state.get_squared_norm(), 1, eps); Eigen::VectorXcd test_state(dim); for (std::uint64_t i = 0; i < dim; ++i) test_state[i] = (CComplex)state_cp[i]; @@ -174,9 +174,9 @@ TEST(StateVectorTest, EntropyCalculation) { for (std::uint64_t ind = 0; ind < dim; ++ind) { CComplex z = test_state[ind]; double prob = z.real() * z.real() + z.imag() * z.imag(); - if (prob > eps) ent += -prob * std::log2(prob); + if (prob > eps) ent += -prob * std::log2(prob); } - ASSERT_NEAR(ent, state.get_entropy(), eps); + ASSERT_NEAR(ent, state.get_entropy(), eps); } } } @@ -190,26 +190,26 @@ TEST(StateVectorTest, GetMarginalProbability) { for (std::uint64_t i = 0; i < dim; ++i) { probs.push_back(internal::squared_norm(state_cp[i])); } - ASSERT_NEAR(state.get_marginal_probability({0, 0}), probs[0], eps); - ASSERT_NEAR(state.get_marginal_probability({1, 0}), probs[1], eps); - ASSERT_NEAR(state.get_marginal_probability({0, 1}), probs[2], eps); - ASSERT_NEAR(state.get_marginal_probability({1, 1}), probs[3], eps); + ASSERT_NEAR(state.get_marginal_probability({0, 0}), probs[0], eps); + ASSERT_NEAR(state.get_marginal_probability({1, 0}), probs[1], eps); + ASSERT_NEAR(state.get_marginal_probability({0, 1}), probs[2], eps); + ASSERT_NEAR(state.get_marginal_probability({1, 1}), probs[3], eps); ASSERT_NEAR(state.get_marginal_probability({0, StateVector::UNMEASURED}), probs[0] + probs[2], - eps); + eps); ASSERT_NEAR(state.get_marginal_probability({1, StateVector::UNMEASURED}), probs[1] + probs[3], - eps); + eps); ASSERT_NEAR(state.get_marginal_probability({StateVector::UNMEASURED, 0}), probs[0] + probs[1], - eps); + eps); ASSERT_NEAR(state.get_marginal_probability({StateVector::UNMEASURED, 1}), probs[2] + probs[3], - eps); + eps); ASSERT_NEAR(state.get_marginal_probability( {StateVector::UNMEASURED, StateVector::UNMEASURED}), 1., - eps); + eps); } TEST(StateVectorTest, SamplingSuperpositionState) { diff --git a/tests/util/util.hpp b/tests/util/util.hpp index e8254b54..4c1c9b08 100644 --- a/tests/util/util.hpp +++ b/tests/util/util.hpp @@ -12,24 +12,30 @@ using scaluq::internal::ComplexMatrix; template using ComplexVector = Eigen::Matrix, -1, 1>; -const inline double eps = 1e-12; -const inline float eps_f = 1e-4; - template -inline void check_near(const StdComplex& a, const StdComplex& b) { +constexpr Fp eps_() { if constexpr (std::is_same_v) - ASSERT_LE(std::abs(a - b), eps); + return 1e-12; + else if constexpr (std::is_same_v) + return 1e-4; else - ASSERT_LE(std::abs(a - b), eps_f); + static_assert(internal::lazy_false_v, "unknown GateImpl"); +} +template +constexpr Fp eps = eps_(); + +template +inline void check_near(const StdComplex& a, const StdComplex& b) { + ASSERT_LE(std::abs(a - b), eps); } template -inline bool same_state(const StateVector& s1, const StateVector& s2, const Fp eps = 1e-12) { +inline bool same_state(const StateVector& s1, const StateVector& s2, const Fp e = eps) { auto s1_cp = s1.get_amplitudes(); auto s2_cp = s2.get_amplitudes(); assert(s1.n_qubits() == s2.n_qubits()); for (std::uint64_t i = 0; i < s1.dim(); ++i) { - if (std::abs((std::complex)s1_cp[i] - (std::complex)s2_cp[i]) > eps) return false; + if (std::abs((std::complex)s1_cp[i] - (std::complex)s2_cp[i]) > e) return false; } return true; }; @@ -37,7 +43,7 @@ inline bool same_state(const StateVector& s1, const StateVector& s2, con template inline bool same_state_except_global_phase(const StateVector& s1, const StateVector& s2, - const Fp eps = 1e-12) { + const Fp e = eps) { auto s1_cp = s1.get_amplitudes(); auto s2_cp = s2.get_amplitudes(); assert(s1.n_qubits() == s2.n_qubits()); @@ -47,16 +53,16 @@ inline bool same_state_except_global_phase(const StateVector& s1, significant = i; } } - if (std::abs((std::complex)s1_cp[significant]) < eps) { + if (std::abs((std::complex)s1_cp[significant]) < e) { for (std::uint64_t i = 0; i < s2.dim(); ++i) { - if (std::abs((std::complex)s2_cp[i]) > eps) return false; + if (std::abs((std::complex)s2_cp[i]) > e) return false; } return true; } Fp phase = std::arg(std::complex(s2_cp[significant] / s1_cp[significant])); std::complex phase_coef = std::polar(1., phase); for (std::uint64_t i = 0; i < s1.dim(); ++i) { - if (std::abs(phase_coef * (std::complex)s1_cp[i] - (std::complex)s2_cp[i]) > eps) + if (std::abs(phase_coef * (std::complex)s1_cp[i] - (std::complex)s2_cp[i]) > e) return false; } return true;