Quantum Classifier¶

# Import numpy,paddle and paddle_quantum
import numpy as np
import paddle
import paddle_quantum

# To construct quantum circuit
from paddle_quantum.ansatz import Circuit
from paddle_quantum.state import zero_state

# Some functions
from numpy import pi as PI
from paddle import matmul, transpose, reshape  # paddle matrix multiplication and transpose
from paddle_quantum.qinfo import pauli_str_to_matrix # N qubits Pauli matrix
from paddle_quantum.linalg import dagger  # complex conjugate

# Plot figures, calculate the run time
from matplotlib import pyplot as plt
import time

# Parameters for generating the data set
Ntrain = 200        # Specify the training set size
Ntest = 100         # Specify the test set size
boundary_gap = 0.5  # Set the width of the decision boundary
seed_data = 2       # Fixed random seed required to generate the data set
# Parameters for training
N = 4               # Number of qubits required
DEPTH = 1           # Circuit depth
BATCH = 20          # Batch size during training
EPOCH = int(200 * BATCH / Ntrain)
                    # Number of training epochs, the total iteration number "EPOCH * (Ntrain / BATCH)" is chosen to be about 200
LR = 0.01           # Set the learning rate
seed_paras = 19     # Set random seed to initialize various parameters

# Generate a binary classification data set with circular decision boundary
def circle_data_point_generator(Ntrain, Ntest, boundary_gap, seed_data):
    """
    :param Ntrain: number of training samples
    :param Ntest: number of test samples
    :param boundary_gap: value in (0, 0.5), means the gap between two labels
    :param seed_data: random seed
    :return: 'Ntrain' samples for training and
             'Ntest' samples for testing
    """
    # Generate "Ntrain + Ntest" pairs of data, x for 2-dim data points, y for labels.
    # The first "Ntrain" pairs are used as training set, the last "Ntest" pairs are used as testing set
    train_x, train_y = [], []
    num_samples, seed_para = 0, 0
    while num_samples < Ntrain + Ntest:
        np.random.seed((seed_data + 10) * 1000 + seed_para + num_samples)
        data_point = np.random.rand(2) * 2 - 1  # 2-dim vector in range [-1, 1]

        # If the modulus of the data point is less than (0.7 - gap), mark it as 0
        if np.linalg.norm(data_point) < 0.7-boundary_gap / 2:
            train_x.append(data_point)
            train_y.append(0.)
            num_samples += 1

        # If the modulus of the data point is greater than (0.7 + gap), mark it as 1
        elif np.linalg.norm(data_point) > 0.7 + boundary_gap / 2:
            train_x.append(data_point)
            train_y.append(1.)
            num_samples += 1
        else:
            seed_para += 1

    train_x = np.array(train_x).astype("float64")
    train_y = np.array([train_y]).astype("float64").T

    print("The dimensions of the training set x {} and y {}".format(np.shape(train_x[0:Ntrain]), np.shape(train_y[0:Ntrain])))
    print("The dimensions of the test set x {} and y {}".format(np.shape(train_x[Ntrain:]), np.shape(train_y[Ntrain:])), "\n")

    return train_x[0:Ntrain], train_y[0:Ntrain], train_x[Ntrain:], train_y[Ntrain:]

def data_point_plot(data, label):
    """
    :param data: shape [M, 2], means M 2-D data points
    :param label: value 0 or 1
    :return: plot these data points
    """
    dim_samples, dim_useless = np.shape(data)
    plt.figure(1)
    for i in range(dim_samples):
        if label[i] == 0:
            plt.plot(data[i][0], data[i][1], color="r", marker="o")
        elif label[i] == 1:
            plt.plot(data[i][0], data[i][1], color="b", marker="o")
    plt.show()

# Construct quantum circuit for data encoding
def encoding_circuit(data, n_qubits):
    cir = Circuit(n_qubits)
    for i in range(n_qubits):
        param = data[0] if i % 2 == 0 else data[1]
        cir.ry(qubits_idx=i, param=np.arcsin(param))
        cir.rz(qubits_idx=i, param=np.arccos(param**2))
    return cir


# Encode the data points to quantum states
def datapoints_transform_to_state(data, n_qubits):
    """
    :param data: shape [BATCHSIZE, 2]
    :param n_qubits: the number of qubits to which
    the data transformed
    :return: shape [BATCHSIZE, 1, 2 ^ n_qubits]
    """

    batchsize, _ = data.shape
    res = []
    init_state = zero_state(n_qubits)
    for i in range(batchsize):
        cir = encoding_circuit(data[i], n_qubits)
        out_state = cir(init_state)
        res.append(out_state.numpy().reshape(1, -1))
    res = np.array(res, dtype=paddle_quantum.get_dtype())
    return res

# Generate Pauli Z operator that only acts on the first qubit
# Act the identity matrix on rest of the qubits
def Observable(n):
    r"""
    :param n: number of qubits
    :return: local observable: Z \otimes I \otimes ...\otimes I
    """
    Ob = pauli_str_to_matrix([[1.0, 'z0']], n)

    return Ob

# Build the computational graph
class Opt_Classifier(paddle_quantum.Operator):
    """
    Construct the model net
    """
    def __init__(self, n, depth, seed_paras=1):
        # Initialization, use n, depth give the initial PQC
        super(Opt_Classifier, self).__init__()
        self.n = n
        self.depth = depth
        paddle.seed(seed_paras)
        # Initialize bias
        self.bias = self.create_parameter(
            shape=[1],
            default_initializer=paddle.nn.initializer.Normal(std=0.01),
            dtype='float32',
            is_bias=False)
        
        self.circuit = Circuit(n)
        # Build a generalized rotation layer
        self.circuit.rz()
        self.circuit.ry()
        self.circuit.rz()

        # Build the entangled layer and Ry rotation layer
        for _ in range(depth):
            # The entanglement layer
            self.circuit.cnot()
            # Add Ry to each qubit
            self.circuit.ry()

    # Define forward propagation mechanism, and then calculate loss function and cross-validation accuracy
    def forward(self, state_in, label):
        """
        Args:
            state_in: The input quantum state, shape [BATCH, 1, 2^n]
            label: label for the input state, shape [BATCH, 1]
        Returns:
            The loss:
                L = 1/BATCH * ((<Z> + 1)/2 + bias - label)^2
        """
        # Convert Numpy array to tensor
        Ob = paddle.to_tensor(Observable(self.n))
        label_pp = reshape(paddle.to_tensor(label), [-1, 1])

        # Build the quantum circuit
        Utheta = self.circuit.unitary_matrix()

        # Because Utheta is achieved by learning, we compute with row vectors to speed up without affecting the training effect
        state_out = matmul(state_in, Utheta)  # shape:[-1, 1, 2 ** n], the first parameter is BATCH in this tutorial

        # Measure the expectation value of Pauli Z operator <Z> -- shape [-1,1,1]
        E_Z = matmul(matmul(state_out, Ob), transpose(paddle.conj(state_out), perm=[0, 2, 1]))

        # Mapping <Z> to the estimated value of the label
        state_predict = paddle.real(E_Z)[:, 0] * 0.5 + 0.5 + self.bias  # |y^{i,k} - \tilde{y}^{i,k}|^2
        loss = paddle.mean((state_predict - label_pp) ** 2)  # Get average for "BATCH" |y^{i,k} - \tilde{y}^{i,k}|^2: L_i：shape:[1,1]

        # Calculate the accuracy of cross-validation
        is_correct = (paddle.abs(state_predict - label_pp) < 0.5).nonzero().shape[0]
        acc = is_correct / label.shape[0]

        return loss, acc, state_predict.numpy(), self.circuit

# Draw the figure of the final training classifier
def heatmap_plot(Opt_Classifier, N):
    # generate data points x_y_
    Num_points = 30
    x_y_ = []
    for row_y in np.linspace(0.9, -0.9, Num_points):
        row = []
        for row_x in np.linspace(-0.9, 0.9, Num_points):
            row.append([row_x, row_y])
        x_y_.append(row)
    x_y_ = np.array(x_y_).reshape(-1, 2).astype("float64")

    # make prediction: heat_data
    input_state_test = paddle.to_tensor(
        datapoints_transform_to_state(x_y_, N))
    loss_useless, acc_useless, state_predict, cir = Opt_Classifier(state_in=input_state_test, label=x_y_[:, 0])
    heat_data = state_predict.reshape(Num_points, Num_points)

    # plot
    fig = plt.figure(1)
    ax = fig.add_subplot(111)
    x_label = np.linspace(-0.9, 0.9, 3)
    y_label = np.linspace(0.9, -0.9, 3)
    ax.set_xticks([0, Num_points // 2, Num_points - 1])
    ax.set_xticklabels(x_label)
    ax.set_yticks([0, Num_points // 2, Num_points - 1])
    ax.set_yticklabels(y_label)
    im = ax.imshow(heat_data, cmap=plt.cm.RdBu)
    plt.colorbar(im)
    plt.show()

def QClassifier(Ntrain, Ntest, gap, N, DEPTH, EPOCH, LR, BATCH, seed_paras, seed_data):
    """
    Quantum Binary Classifier
    Input:
        Ntrain         # Specify the training set size
        Ntest          # Specify the test set size
        gap            # Set the width of the decision boundary
        N              # Number of qubits required
        DEPTH          # Circuit depth
        BATCH          # Batch size during training
        EPOCH          # Number of training epochs, the total iteration number "EPOCH * (Ntrain / BATCH)" is chosen to be about 200
        LR             # Set the learning rate
        seed_paras     # Set random seed to initialize various parameters
        seed_data      # Fixed random seed required to generate the data set
        plot_heat_map  # Whether to plot heat map, default True
    """
    # Generate data set
    train_x, train_y, test_x, test_y = circle_data_point_generator(Ntrain=Ntrain, Ntest=Ntest, boundary_gap=gap, seed_data=seed_data)
    # Read the dimension of the training set
    N_train = train_x.shape[0]
    
    # Initialize the registers to store the accuracy rate and other information
    summary_iter, summary_test_acc = [], []

    # Generally, we use Adam optimizer to get relatively good convergence
    # Of course, it can be changed to SGD or RMSprop
    myLayer = Opt_Classifier(n=N, depth=DEPTH, seed_paras=seed_paras)  # Initial PQC
    opt = paddle.optimizer.Adam(learning_rate=LR, parameters=myLayer.parameters())


    # Optimize iteration
    # We divide the training set into "Ntrain/BATCH" groups
    # For each group the final circuit will be used as the initial circuit for the next group
    # Use cir to record the final circuit after learning.
    i = 0  # Record the iteration number
    for ep in range(EPOCH):
        # Learn for each group
        for itr in range(N_train // BATCH):
            i += 1  # Record the iteration number
            # Encode classical data into a quantum state |psi>, dimension [BATCH, 2 ** N]
            input_state = paddle.to_tensor(datapoints_transform_to_state(train_x[itr * BATCH:(itr + 1) * BATCH], N))

            # Run forward propagation to calculate loss function
            loss, train_acc, state_predict_useless, cir \
                = myLayer(state_in=input_state, label=train_y[itr * BATCH:(itr + 1) * BATCH])  # optimize the given PQC
            # Print the performance in iteration
            if i % 30 == 5:
                # Calculate the correct rate on the test set test_acc
                input_state_test = paddle.to_tensor(datapoints_transform_to_state(test_x, N))
                loss_useless, test_acc, state_predict_useless, t_cir \
                    = myLayer(state_in=input_state_test,label=test_y)
                print("epoch:", ep, "iter:", itr,
                      "loss: %.4f" % loss.numpy(),
                      "train acc: %.4f" % train_acc,
                      "test acc: %.4f" % test_acc)
                # Store accuracy rate and other information
                summary_iter.append(itr + ep * N_train)
                summary_test_acc.append(test_acc) 

            # Run back propagation to minimize the loss function
            loss.backward()
            opt.minimize(loss)
            opt.clear_grad()
            
    # Print the final circuit
    print("The trained circuit:")
    print(cir)
    # Draw the decision boundary represented by heatmap
    heatmap_plot(myLayer, N=N)

    return summary_test_acc

# A simpler classifier
def QClassifier2(quantum_train_x, train_y,quantum_test_x,test_y, N, DEPTH, EPOCH, LR, BATCH):
    """
    Quantum Binary Classifier
    Input：
        quantum_train_x     # training x
        train_y             # training y
        quantum_test_x      # testing x
        test_y              # testing y
        N                   # Number of qubits required
        DEPTH               # Circuit depth
        EPOCH               # Number of training epochs
        LR                  # Set the learning rate
        BATCH               # Batch size during training
    """
    Ntrain = len(quantum_train_x)
    
    paddle.seed(1)

    net = Opt_Classifier(n=N, depth=DEPTH)

    # Test accuracy list
    summary_iter, summary_test_acc = [], []

    # Adam can also be replaced by SGD or RMSprop
    opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())

    # Optimize
    for ep in range(EPOCH):
        for itr in range(Ntrain // BATCH):
            # Import data
            input_state = quantum_train_x[itr * BATCH:(itr + 1) * BATCH]  # paddle.tensor
            input_state = reshape(input_state, [-1, 1, 2 ** N])
            label = train_y[itr * BATCH:(itr + 1) * BATCH]
            test_input_state = reshape(quantum_test_x, [-1, 1, 2 ** N])

            loss, train_acc, state_predict_useless, cir = net(state_in=input_state, label=label)

            if itr % 5 == 0:
                # get accuracy on test dataset (test_acc)
                loss_useless, test_acc, state_predict_useless, t_cir = net(state_in=test_input_state, label=test_y)
                print("epoch:", ep, "iter:", itr,
                      "loss: %.4f" % loss.numpy(),
                      "train acc: %.4f" % train_acc,
                      "test acc: %.4f" % test_acc)
                summary_test_acc.append(test_acc)

            loss.backward()
            opt.minimize(loss)
            opt.clear_grad()

    return summary_test_acc