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

# # Quantum CSD Compiling Intro
# 
# The purpose of this notebook is to give a quick introduction to Qubiter's 
# CSD quantum compiler capabilities. 
# 
# By a quantum complier, we mean
# a computer program that takes as input an arbitrary unitary matrix $U$ of dimension $N=2^n$
# and returns a SEO (sequence of elementary operations, in this case CNOTs and single qubit
# rotations) that equals $U$. There are various kinds of quantum 
# compilers. Suppose $U$ is of the form $U=e^{-itH}$, where $t$ is time and $H$ is
# the Hamiltonian operator. If $H$ has a form which is known a priori, a situation
# that is common in Physics and Chemistry, then a popular way of expanding $U$
# is by using the Trotter-Suzuki approximation. If the form of $H$ is not
# known a priori as is common in Artificial Intelligence, then
# we recommend using the CS (Cosine-Sine) decomposition of Linear Algebra.
# Both methods are already implemented in Qubiter, but this notebook is about
# the CSD.
# 
# Doing CSD quantum compiling with Qubiter requires using the classes in the quantum_CSD_compiler
# folder, which will only work properly if you install, besides all the Qubiter
# Python files and a Python distro that includes numpy and scipy (for example, Anaconda),
# some binary libraries prepared by Artiste-q.net which include
# a Python wrapper for a LAPACK subroutine
# called cuncsd.f that performs CSDs. How to install those binary libraries
# is explained elsewhere in this site. Henceforth, we will assume 
# all the necessary files have been installed on your computer if you want to redo the calculations.
# 
# The quantum_CSD_compiler folder includes a pdf called csd-intro.pdf that gives
# an introduction to the CSD. 
# 
# Some external references:
# 
# 
# 1. R.R. Tucci, A Rudimentary Quantum Compiler(2cnd Ed.)
#     https://arxiv.org/abs/quant-ph/9902062
# 
# 2. Qubiter 1.11, a C++ program whose first version was released together
#     with Ref.1 above. Qubiter 1.11 is included in the
#     quantum_CSD_compiler/LEGACY folder of this newer, pythonic version of Qubiter
#     
# 3. R.R. Tucci, Quantum Fast Fourier Transform Viewed as a Special Case of Recursive Application of Cosine-Sine Decomposition, https://arxiv.org/abs/quant-ph/0411097
# 
# 

# Qubiter applies CSD recursively
# to build a tree of node matrices. The product of those node matrices,
# if read in the correct order, is equal to the input matrix $U$.
# 
# As an example, let us use for $U$ a 3 qubit quantum Fourier matrix.
# We can create an object of class Tree with $U$ 
# as input as follows

# 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_qbits = 3
init_unitary_mat = FouSEO_writer.fourier_trans_mat(1 << num_qbits)
emb = CktEmbedder(num_qbits, num_qbits)
file_prefix = 'csd_test'
t = Tree(True, file_prefix, emb, init_unitary_mat, verbose=False)
t.close_files()


# The above code automatically creates an expansion of $U$ 
# into DIAG and MP_Y lines. Next we print the Picture file that was created.

# In[4]:


t.print_pic_file(jup=True)


# Let us also print the corresponding English file that was created.

# In[5]:


t.print_eng_file(jup=True)


# In the Picture files, a DIAG
# appears as a chain of percents, whereas an MP_Y appears as a chain of percents
# and  an additional Ry gate. As you can
# see, the Picture file gives a nice picture of the DIAG and MP_Y gates,
# but the English file is much more specific. Look at Qubiter's Rosetta Stone
# (qubiter_rosetta_stone.pdf)
# if you want to understand how to interpret 
# the parameters of DIAG and MP_Y lines.
# 
# A DIAG represents a diagonal matrix with diagonal entries that
# are unit magnitude complex numbers, making the matrix unitary.
# 
# An MP_Y represents a matrix of the form
# 
# $\left[\begin{array}{cc} cc & ss \\ -ss & cc \end{array}\right]$
# 
# where $cc$ and $ss$ are real diagonal matrices of the same size 
# such that $cc^2 + ss^2 = 1$.

# The class DiagUnitaryExpander
# will take any English file and 
# write new English and Picture files wherein 
# 
# * all lines except those starting with DIAG are echoed,
# * lines starting with DIAG are replaced by an exact or approximate
# multiline expansion. 
# 
# Likewise, the class MultiplexorExpander
# will take any English file and 
# write new English and Picture files wherein 
# 
# * all lines except those starting with MP_Y are echoed,
# * lines starting with MP_Y are replaced by an exact or approximate
# multiline expansion. 
# 
# We could use these 2 expander
# classes to construct new English and Picture files from the English 
# file printed above. This would lead to an English file
# that consisted of only CNOTs and qubit rotations. If the 
# gates of that new English file were multiplied,
# the product would equal the original $U$. Such an English file would
# be very long and not too instructive so we won't show it here.
# Instead, we will show an exact expansion of a single DIAG and 
# of a single MP_Y.
# 

# Next, we create English and Picture files containing
# an expansion of the 4 qubit gate
# 
# %---%---%---%
# 
# This represents a diagonal unitary matrix. The 
# angles are chosen at random and stored in the variable rad_angles.
# We then print the Picture file.

# In[6]:


file_prefix = "d_unitary_exact_check"
num_qbits = 4
num_angles = (1 << num_qbits)
emb = CktEmbedder(num_qbits, num_qbits)
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)


# We can check that our exact expansion is correct
# as follows. We can multiply the gates of the expansion
# using the class SEO_MatrixProduct. Call the gate product matpro.prod_arr.
# Using the angles rad_angles that we stored,
# we can construct the exact diagonal unitary, call it exact_mat.
# Call err the norm of matpro.prod_arr - exact_mat,
# and print err.

# In[7]:


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


# Next, we create English and Picture files containing
# an expansion of the 4 qubit gate
# 
# Ry--%---%---%
# 
# This represents a multiplexor matrix. The 
# angles are chosen at random and stored in the variable rad_angles.
# We then print the Picture file.

# In[8]:


file_prefix = "plexor_exact_check"
num_qbits = 4
num_angles = (1 << (num_qbits-1))
emb = CktEmbedder(num_qbits, num_qbits)
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)


# We can check that our exact expansion is correct
# as follows. We can multiply the gates of the expansion
# using the class SEO_MatrixProduct. Call the gate product matpro.prod_arr.
# Using the angles rad_angles that we stored,
# we can construct the exact multiplexor matrix, call it exact_mat.
# Call err the norm of matpro.prod_arr - exact_mat,
# and print err.

# In[9]:


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


# A moral of the above calculations is that using  CSD
# quantum compiling blindly will give a SEO for a quantum Fourier Transform QFT that is exponential in the number of qubits $n$. And yet we know that Coppersmith came up with an expansion for the QFT that is polynomial in $n$. But there is hope: CSD is not a unique decomposition.
# Ref.3 explains how one can coax a CSD compiler to yield Coppersmith's decompostion.
# 
# Let $U$ be N dimensional, with $N = 2^n$, where $n$ = number of
# qubits. A general N dimensional unitary matrix has $N^2$ dofs (real
# degrees of freedom). That's because it has $N^2$ complex entries, so $2N^2$
# real parameters, but those parameters are subject to N real constraints
# and N(N-1)/2 complex constraints, for a total of $N^2$ real constraints.
# So $2N^2$ real parameters minus N^2 real constraints gives $N^2$ dofs.
# 
# (a) Each DIAG (MP_Y, resp.) line of the CS decomp of $U$
# depends on N (N/2, resp.) angles and there are about N DIAG and N MP_Y
# lines. So the DIAG lines alone have enough dofs, $N^2$ of them, to cover
# all $N^2$ dofs of $U$. So clearly, there is a lot of
# redundancy in the CS decomp used by Qubiter. But, there is hope: the CS
# decomp is not unique, and it might be possible to choose a CS decomp
# that makes zero many of the angles in the DIAG and MP_Y lines. 
# 
# (b) The CS decomp as used here leads to order $N^2 = 2^{2n}$ cnots and
# qubit rotations so it is impractical for large N. But for small N,
# it can be useful. For large N, it might be possible to discover
# approximations to individual MP_Y and DIAG lines.
#     
# Clearly, there is much room for future research to improve (a) and (b).
# 

# In[ ]: