#!/usr/bin/env python
# coding: utf-8

# # 量子编译简介
# 
# 本文档介绍 *Qubiter* 所用的正余弦分解（CSD, Cosine-Sine decomposition）量子编译器及其性能。
# 
# 量子编译器指以任意 $N=2^n$ 维的幺正矩阵 $U$ 作为输入，能够输出等价于矩阵 $U$ 的基本操作序列（比如 CNOT 受控反门和单量子比特旋转门）的计算机程序，可以有很多种。通常幺正矩阵 $U$ 的形式为 $U=e^{-itH}$，其中 $t$ 和 $H$ 分别代表时间和哈密顿算符。在物理学和量子化学领域中，$H$ 的具体形式往往可以先验已知，通常的做法是用 Trotter-Suzuki 近似对 $U$ 进行展开。而在人工智能等领域，哈密顿量 $H$ 的具体形式一般不能先验已知，这时候则推荐使用线性代数中正余弦分解（CSD）来处理矩阵。这两种不同的编译方式都可以使用 *Qubiter* 实现，此文档着重介绍正余弦分解编译方法。
# 
# 使用 *Qubiter* 进行 CSD 量子编译需要调用 'quantum_CSD_compiler' 文件夹中的类文件。这些类的使用除了需要 *Qubiter* 的 Python 文件包和安装了 numpy、scipy 库的 Python 环境（例如 Anaconda ）外，还需要[安装由 Artiste-qb.net 提供的二进制库](https://github.com/artiste-qb-net/Python-CS-Decomposition#installation)。此库包括了 Python 封装的进行正余弦分解的 LAPACK 子程序 'cuncsd.f'。复述本文档中的计算需要安装所有必要的文件和库。
# 
# 在 'quantum_CSD_compiler' 文件夹中的 ['csd-intro.pdf'](https://github.com/artiste-qb-net/qubiter/blob/master/quantum_CSD_compiler/csd-intro.pdf) 文件提供了 CSD（正余弦分解）的介绍。其他参考文献有：
# 
# 1. [R.R. Tucci,《初级量子编译》（第二版）](https://arxiv.org/abs/quant-ph/9902062)
# 
# 2. Qubiter 1.11: 基于 C++ 的原始 *Qubiter* 版本。最初一版 *Qubiter 1.11* 和文献1同步发布。在新的 *Qubiter* 版本（基于 Python）的 'quantum_CSD_compiler/LEGACY' 文件夹中包含了 *Qubiter 1.11* 原始程序。
# 
# 3. R.R. Tucci,[《量子快速傅立叶变换 - 正余弦分解递归应用举例》](https://arxiv.org/abs/quant-ph/0411097)
# 

# *Qubiter* 通过递归调用 CSD 来建立由节点矩阵构成的树结构。这些节点矩阵通过正确的顺序读取后是和输入矩阵 $U$ 完全等价的。
# 
# 在这里以三个量子比特的量子傅立叶矩阵 $U$ 为例，用如下方式构建类树的对象：

# In[1]:


import os
import sys
print(os.getcwd())
os.chdir('../../')
print(os.getcwd())
sys.path.insert(0,os.getcwd())


# In[2]:


import numpy as np
import cuncsd_sq as csd

import math
from qubiter.FouSEO_writer import *
from qubiter.quantum_CSD_compiler.Tree import *
from qubiter.quantum_CSD_compiler.DiagUnitarySEO_writer import *
from qubiter.quantum_CSD_compiler.MultiplexorSEO_writer import *


# In[3]:


num_bits = 3
init_unitary_mat = FouSEO_writer.fourier_trans_mat(1 << num_bits)
emb = CktEmbedder(num_bits, num_bits)
file_prefix = 'csd_test'
t = Tree(True, file_prefix, emb, init_unitary_mat, verbose=False)
t.close_files()


# 上述代码自动将矩阵 $U$ 展开并创建 DIAG 和 MP_Y 线路，其图形化文件（Picture file）为：

# In[4]:


t.print_pic_file(jup=True)


# 
# 相应的，创建的文字化文件（English file）为：

# In[5]:


t.print_eng_file(jup=True)


# 在图形化文件中， 一个 DIAG 线路作为一个百分号链出现，而一个 MP_Y 线路作为百分号链加 Ry 门出现。我们可以看到，相对于较复杂的文字化文件，图形化文件可以很直观地给出 DIAG 和 MP_Y 门的概览。在 *Qubiter* 的诠释文件 [qubiter_rosetta_stone.pdf](https://github.com/artiste-qb-net/qubiter/blob/master/qubiter_rosetta_stone.pdf) 中可以了解更多 DIAG 和 MP_Y 线路的解释和举例。
# 
# 在这里，每一个 DIAG 线路都代表一个由单位尺度大小的复数构成的对角矩阵，从而保证了矩阵的幺正性。同时每一个 MP_Y 线路代表了具有如下形式的矩阵，
# 
# $\left[\begin{array}{cc} cc & ss \\ -ss & cc \end{array}\right]$，
# 
# 其中 $cc$ and $ss$ 是相同大小的实对角矩阵从而满足 $cc^2 + ss^2 = 1$.
# 

# 对每一个输入的文字化文件，'DiagUnitaryExpander' 类都会生成新的文字化文件和图形化文件，其中：
# 
# * 所有非 'DIAG' 开始的行都被回显，
# * 所有以 'DIAG' 开始的行都被其精确或近似的多线路展开所替换。
# 
# 同样的，'MultiplexorExpander' 也将任一输入的文字化文件输出为新的文字化文件和图形化文件，其中：
# 
# * 所有非 'MP_Y' 开始的行都被回显，
# * 所有以 'MP_Y' 开始的行都被其精确或近似的多线路展开所替换。
# 
# 所有之前的文字化文件都可以被这两个类展开并且构建为新的文字化和图形化文件，最终会得到一个只包含 CNOT 门和单量子比特旋转的文字化文件。而此文件中的逻辑门相乘后和原始矩阵 $U$ 是完全等价的。鉴于最终的文字化文件较长且不适合展示，这里仅以单个 DIAG 和单个 MP_Y 线路的精确展开为例做演示：
# 

# 首先对于一个四量子比特门
# 
# %---%---%---%
# 
# 我们创建包含其展开的文字化文件和图形化文件。这代表一个幺正对角矩阵，其相位角是随机生成的并存储为变量 'rad_angles'。其图形化文件为：

# In[6]:


file_prefix = "d_unitary_exact_check"
num_bits = 4
num_angles = (1 << num_bits)
emb = CktEmbedder(num_bits, num_bits)
rad_angles = list(np.random.rand(num_angles)*2*np.pi)
wr = DiagUnitarySEO_writer(file_prefix, emb, 'exact', rad_angles)
wr.write()
wr.close_files()
wr.print_pic_file(jup=True)


# 可以通过如下方式检验此精确展开的正确性：将展开后的门操作用 'SEO_MatrixProduct' 类进行相乘，并存储乘积为 'matpro.prod_arr'。利用 'rad_angles' 中的相位角精确重建对角矩阵并存为 'exact_mat'。将两结果差值（即 matpro.prod_arr - exact_mat）的标准差存储为 'err'　并输出：
# 

# In[7]:


matpro = SEO_MatrixProduct(file_prefix, num_bits)
exact_mat = DiagUnitarySEO_writer.du_mat(rad_angles)
err = np.linalg.norm(matpro.prod_arr - exact_mat)
print("diag unitary error=", err)


# 接下来，对另一个四量子比特门
# 
# Ry--%---%---%
# 
# 我们构建展开后的文字化文件和图形化文件。这是一个多路复用器矩阵（multiplexor matrix），其随机相位角存储为变量 'rad_angles'。其图形化文件为：
# 

# In[8]:


file_prefix = "plexor_exact_check"
num_bits = 4
num_angles = (1 << (num_bits-1))
emb = CktEmbedder(num_bits, num_bits)
rad_angles = list(np.random.rand(num_angles)*2*np.pi)
wr = MultiplexorSEO_writer(file_prefix, emb, 'exact', rad_angles)
wr.write()
wr.close_files()
wr.print_pic_file(jup=True)


# 同前面一样，可以通过如下方式检验此精确展开的正确性：将展开后的门操作用 'SEO_MatrixProduct' 进行相乘，并存储乘积为 'matpro.prod_arr'。利用 'rad_angles' 中的相位角精确重建对角矩阵并存为 'exact_mat'。将两结果差值（即 matpro.prod_arr - exact_mat）的标准差存储为 'err'　并输出：

# In[9]:


matpro = SEO_MatrixProduct(file_prefix, num_bits)
exact_mat = MultiplexorSEO_writer.mp_mat(rad_angles)
err = np.linalg.norm(matpro.prod_arr - exact_mat)
print("multiplexor error=", err)


# 上述计算意味着，对量子傅立叶变换（QFT）盲目使用 CSD 量子编译会得到一个随量子比特数 $n$ 指数增长的操作序列。而我们知道对 QFT 使用 Coppersmith 分解仅得到随着 $n$ 多项式增长的结果。不过由于 CSD 的非惟一性，在文献3中介绍了如何用 CSD 编译以得到优于 Coppersmith 分解的方法。
# 
# 对于 $n$ 个量子比特，$U$ 是一个 $N = 2^n$ 维的幺正矩阵并且具有 $N^2$ 个自由度（真实自由度）。这是由于 $U$ 具有 $N^2$ 个复数参量（即 $2N^2$ 个实数参量）以及 $N$ 个实约束条件和 $N(N-1)/2$ 个复约束条件（即总共 $N^2$ 个实约束条件）。真实自由度为 $2N^2$ 个实数参量减去 $N^2$ 个实约束条件，即为 $N^2$ 个自由度。
# 
# (a) 矩阵 $U$ 的正余弦分解含有 $N$ 个DIAG 和 $N$ 个 MP_Y 线路，而每个 DIAG（或 MP_Y）线路都依赖于 $N$ (或 $N/2$) 个相位角。仅 DIAG 线路就具有足够多的 $N^2$ 自由度来构建 $U$ 的 $N^2$ 个自由度。显然 *Qubiter* 所用的 CSD 会产生冗余量。不过由于正余弦分解的非唯一性，可以通过寻找在 DIAG 和 MP_Y 线路中不产生额外相位角的正余弦分解来解决这个问题。
# 
# (b) 此正余弦分解会得到 $N^2 = 2^{2n}$ 量级的 CNOT 门和单量子比特旋转，所以对于 $N$ 比较小的情况比较适合。而当 $N$ 非常大的情况下，此方法不再适用。此时可以通过寻找对单个 MP_Y 和 DIAG 线路的近似来简化问题。
# 
# 显然，对于问题（a）和（b）还有很大的改进空间。