Skip to content

Commit

Permalink
🔀 Merge pull request #95 from DVLab-NTU/bug/fix-template
Browse files Browse the repository at this point in the history 
Bug/fix template
  • Loading branch information
@JoshuaLau0220
JoshuaLau0220 authored Mar 13, 2024
2 parents 175a992 + 948bb6a commit 8a301cd
Show file tree
Hide file tree
Showing 5 changed files with 83 additions and 101 deletions.
53 changes: 28 additions & 25 deletions src/tensor/decomposer.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,12 @@ TwoLevelMatrix<T> adjoint(TwoLevelMatrix<T> m /* copy on purpose */) {
return m;
}

template <typename T>
struct ZYZ {
double phi;
double alpha;
double beta; // actual beta/2
double gamma;
T phi;
T alpha;
T beta; // actual beta/2
T gamma;
bool correct = true;
};

Expand Down Expand Up @@ -81,7 +82,7 @@ class Decomposer {
}

template <typename U>
std::optional<std::pair<size_t, size_t>> _get_two_level_matrix_indices(QTensor<U> const& matrix, double eps);
std::optional<std::pair<size_t, size_t>> _get_two_level_matrix_indices(QTensor<U> const& matrix, U eps);
template <typename U>
std::vector<TwoLevelMatrix<U>> _get_two_level_matrices(QTensor<U> matrix /* copy on purpose */);

Expand All @@ -100,7 +101,7 @@ class Decomposer {
bool _decompose_cu(Tensor<U> const& t, size_t ctrl, size_t targ);

template <typename U>
ZYZ _decompose_zyz(Tensor<U> const& matrix);
ZYZ<typename U::value_type> _decompose_zyz(Tensor<U> const& matrix);

template <typename U>
Tensor<U> _sqrt_single_qubit_matrix(Tensor<U> const& matrix);
Expand Down Expand Up @@ -146,7 +147,7 @@ std::optional<QCir> Decomposer::decompose(QTensor<U> const& matrix) {
* @return std::optional<std::pair<size_t, size_t>>
*/
template <typename U>
std::optional<std::pair<size_t, size_t>> Decomposer::_get_two_level_matrix_indices(QTensor<U> const& matrix, double eps) {
std::optional<std::pair<size_t, size_t>> Decomposer::_get_two_level_matrix_indices(QTensor<U> const& matrix, U eps) {
using namespace std::literals;
auto const dimension = static_cast<size_t>(matrix.shape()[0]);
size_t num_found_diagonal = 0, num_upper_triangle_not_zero = 0, num_lower_triangle_not_zero = 0;
Expand Down Expand Up @@ -238,7 +239,7 @@ std::optional<std::pair<size_t, size_t>> Decomposer::_get_two_level_matrix_indic
template <typename U>
std::vector<TwoLevelMatrix<U>> Decomposer::_get_two_level_matrices(QTensor<U> matrix /* copy on purpose */) {
using namespace std::literals;
constexpr double eps = 1e-6;
constexpr U eps = 1e-6;
std::vector<TwoLevelMatrix<U>> two_level_chain;
auto const dimension = static_cast<size_t>(matrix.shape()[0]);

Expand Down Expand Up @@ -272,7 +273,7 @@ std::vector<TwoLevelMatrix<U>> Decomposer::_get_two_level_matrices(QTensor<U> ma
}

// normalization factor
const double u = std::sqrt(std::norm(matrix(i, i)) + std::norm(matrix(j, i)));
const U u = std::sqrt(std::norm(matrix(i, i)) + std::norm(matrix(j, i)));

auto const n_qubits = static_cast<size_t>(std::round(std::log2(_get_dimension(matrix))));
QTensor<U> conjugate_matrix_product =
Expand Down Expand Up @@ -442,8 +443,9 @@ bool Decomposer::_decompose_cnx(const std::vector<size_t>& ctrls, const size_t e
*/
template <typename U>
bool Decomposer::_decompose_cu(Tensor<U> const& t, size_t ctrl, size_t targ) {
constexpr double eps = 1e-6;
const ZYZ angles = _decompose_zyz(t);
using float_type = typename U::value_type;
constexpr float_type eps = 1e-6;
const ZYZ<float_type> angles = _decompose_zyz(t);
if (!angles.correct) return false;

if (std::abs((angles.alpha - angles.gamma) / 2) > eps)
Expand Down Expand Up @@ -485,31 +487,32 @@ bool Decomposer::_decompose_cu(Tensor<U> const& t, size_t ctrl, size_t targ) {
* @reference Nakahara, Mikio, and Tetsuo Ohmi. Quantum computing: from linear algebra to physical realizations. CRC press, 2008.
*/
template <typename U>
ZYZ Decomposer::_decompose_zyz(Tensor<U> const& matrix) {
ZYZ<typename U::value_type> Decomposer::_decompose_zyz(Tensor<U> const& matrix) {
DVLAB_ASSERT(matrix.shape()[0] == 2 && matrix.shape()[1] == 2, "decompose_ZYZ only supports 2x2 matrix");
using namespace std::literals;
using float_type = typename U::value_type;
// a = e^{iφ}e^{-i(α+γ)/2}cos(β/2)
// b = -e^{iφ}e^{-i(α-γ)/2}sin(β/2)
// c = e^{iφ}e^{ i(α-γ)/2}sin(β/2)
// d = e^{iφ}e^{ i(α+γ)/2}cos(β/2)
const std::complex<double> a = matrix(0, 0), b = matrix(0, 1), c = matrix(1, 0), d = matrix(1, 1);
ZYZ output = {};
const U a = matrix(0, 0), b = matrix(0, 1), c = matrix(1, 0), d = matrix(1, 1);
ZYZ<float_type> output = {};
// NOTE - The beta here is actually half of beta
double init_beta = 0;
float_type init_beta = 0;
if (std::abs(a) > 1) {
init_beta = 0;
} else {
init_beta = std::acos(std::abs(a));
}

constexpr auto pi = std::numbers::pi_v<double>;
constexpr auto pi = std::numbers::pi_v<float_type>;
// NOTE - Possible betas due to arccosine
const std::array<double, 4> beta_candidate = {init_beta, pi - init_beta, pi + init_beta, 2.0 * pi - init_beta};
const std::array<float_type, 4> beta_candidate = {init_beta, pi - init_beta, pi + init_beta, 2.0 * pi - init_beta};
for (const auto& beta : beta_candidate) {
output.beta = beta;
std::complex<double> a1, b1, c1, d1;
const std::complex<double> cos(std::cos(beta) + 1e-5, 0); // cos(β/2)
const std::complex<double> sin(std::sin(beta) + 1e-5, 0); // sin(β/2)
U a1, b1, c1, d1;
const U cos(std::cos(beta) + 1e-5, 0); // cos(β/2)
const U sin(std::sin(beta) + 1e-5, 0); // sin(β/2)
a1 = a / cos;
b1 = b / sin;
c1 = c / sin;
Expand All @@ -525,15 +528,15 @@ ZYZ Decomposer::_decompose_zyz(Tensor<U> const& matrix) {
output.gamma = std::arg(d1 / c1);
}

auto const alpha_plus_gamma = std::exp(std::complex<double>((0.5i) * (output.alpha + output.gamma)));
auto const alpha_minus_gamma = std::exp(std::complex<double>((0.5i) * (output.alpha - output.gamma)));
auto const alpha_plus_gamma = std::exp(U((0.5i) * (output.alpha + output.gamma)));
auto const alpha_minus_gamma = std::exp(U((0.5i) * (output.alpha - output.gamma)));

if (std::abs(a) < 1e-4)
output.phi = std::arg(c1 / alpha_minus_gamma);
else
output.phi = std::arg(a1 * alpha_plus_gamma);

const std::complex<double> phi(std::cos(output.phi), std::sin(output.phi));
const U phi(std::cos(output.phi), std::sin(output.phi));

if (std::abs(phi * cos / alpha_plus_gamma - a) < 1e-3 &&
std::abs(sin * phi / alpha_minus_gamma + b) < 1e-3 &&
Expand Down Expand Up @@ -561,8 +564,8 @@ Tensor<U> Decomposer::_sqrt_single_qubit_matrix(Tensor<U> const& matrix) {
DVLAB_ASSERT(matrix.shape()[0] == 2 && matrix.shape()[1] == 2, "sqrt_single_qubit_matrix only supports 2x2 matrix");
// a b
// c d
const std::complex<double> a = matrix(0, 0), b = matrix(0, 1), c = matrix(1, 0), d = matrix(1, 1);
const std::complex<double> tau = a + d, delta = a * d - b * c;
const U a = matrix(0, 0), b = matrix(0, 1), c = matrix(1, 0), d = matrix(1, 1);
const U tau = a + d, delta = a * d - b * c;
const std::complex s = std::sqrt(delta);
const std::complex t = std::sqrt(tau + 2. * s);
if (std::abs(t) > 0) {
Expand Down
6 changes: 2 additions & 4 deletions src/tensor/qtensor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ class QTensor : public Tensor<std::complex<T>> {
friend dvlab::Phase global_phase(QTensor<U> const& t1, QTensor<U> const& t2);

template <typename U>
friend bool is_equivalent(QTensor<U> const& t1, QTensor<U> const& t2, double eps /* = 1e-6*/);
friend bool is_equivalent(QTensor<U> const& t1, QTensor<U> const& t2, U eps /* = 1e-6*/);

void set_filename(std::string const& f) { _filename = f; }
void add_procedures(std::vector<std::string> const& ps) { _procedures.insert(_procedures.end(), ps.begin(), ps.end()); }
Expand Down Expand Up @@ -333,8 +333,6 @@ QTensor<T> QTensor<T>::to_qtensor() const {

template <typename T>
QTensor<T> QTensor<T>::to_su2() const {
// std::complex<double> det = u.determinant();
// std::complex<double> one(1, 0);
return std::sqrt(1.0 / this->determinant()) * this->_tensor;
}

Expand Down Expand Up @@ -415,7 +413,7 @@ dvlab::Phase global_phase(QTensor<U> const& t1, QTensor<U> const& t2) {
}

template <typename U>
bool is_equivalent(QTensor<U> const& t1, QTensor<U> const& t2, double eps = 1e-6) {
bool is_equivalent(QTensor<U> const& t1, QTensor<U> const& t2, U eps = 1e-6) {
if (t1.shape() != t2.shape()) return false;
return cosine_similarity(t1, t2) >= (1 - eps);
}
Expand Down
18 changes: 5 additions & 13 deletions src/tensor/solovay_kitaev.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,26 +14,18 @@ namespace qsyn::tensor {
/**
* @brief Initialize binary list
*
* @return BinaryList
*/
BinaryList SolovayKitaev::_init_binary_list() const {
BinaryList bin_list;
for (size_t i = 1; i <= _depth; i++) {
for (size_t j = 0; j < gsl::narrow<size_t>(std::pow(2, i)); j++) {
sul::dynamic_bitset<> bitset(i, j);
std::vector<bool> bit_vector;
for (size_t k = 0; k < i; ++k)
bit_vector.emplace_back(bitset[k]);
bin_list.emplace_back(bit_vector);
}
}
return bin_list;
void SolovayKitaev::_init_binary_list() {
for (size_t i = 1; i <= _depth; i++)
for (size_t j = 0; j < gsl::narrow<size_t>(std::pow(2, i)); j++)
_binary_list.emplace_back(sul::dynamic_bitset<>(i, j));
}

/**
* @brief Adjoint the gate sequence
*
* @param sequence
* @return std::vector<int>
*/
std::vector<int> SolovayKitaev::_adjoint_gate_sequence(std::vector<int> sequence) {
std::reverse(sequence.begin(), sequence.end());
Expand Down
Loading

0 comments on commit 8a301cd

Please sign in to comment.