{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Variational Quantum Singular Value Decomposition\n",
"\n",
" Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"In this tutorial, we will go through the concept of classical singular value decomposition (SVD) and the quantum neural network (QNN) version of variational quantum singular value decomposition (VQSVD) [1]. The tutorial consists of the following two parts: \n",
"- Decompose a randomly generated $8\\times8$ complex matrix; \n",
"- Apply SVD on image compression."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\n",
"\n",
"Singular value decomposition (SVD) has many applications, including principal component analysis (PCA), solving linear equations and recommender systems. The main task is formulated as following:\n",
"> Given a complex matrix $M \\in \\mathbb{C}^{m \\times n}$, find the decomposition in form $M = UDV^\\dagger$, where $U_{m\\times m}$ and $V^\\dagger_{n\\times n}$ are unitary matrices, which satisfy the property $UU^\\dagger = VV^\\dagger = I$.\n",
"\n",
"- The column vectors $|u_j\\rangle$ of the unitary matrix $U$ are called left singular vectors $\\{|u_j\\rangle\\}_{j=1}^{m}$ form an orthonormal basis. These column vectors are the eigenvectors of the matrix $MM^\\dagger$.\n",
"- Similarly, the column vectors $\\{|v_j\\rangle\\}_{j=1}^{n}$ of the unitary matrix $V$ are the eigenvectors of $M^\\dagger M$ and form an orthonormal basis.\n",
"- The diagonal elements of the matrix $D_{m\\times n}$ are singular values $d_j$ arranged in a descending order.\n",
"\n",
"For the convenience, we assume that the $M$ appearing below are all square matrices. Let's first look at an example: \n",
"\n",
"$$\n",
"M = 2*X\\otimes Z + 6*Z\\otimes X + 3*I\\otimes I = \n",
"\\begin{bmatrix} \n",
"3 &6 &2 &0 \\\\\n",
"6 &3 &0 &-2 \\\\\n",
"2 &0 &3 &-6 \\\\\n",
"0 &-2 &-6 &3 \n",
"\\end{bmatrix}, \\tag{1}\n",
"$$\n",
"\n",
"Then the singular value decomposition of the matrix can be expressed as:\n",
"\n",
"$$\n",
"M = UDV^\\dagger = \n",
"\\frac{1}{2}\n",
"\\begin{bmatrix} \n",
"-1 &-1 &1 &1 \\\\\n",
"-1 &-1 &-1 &-1 \\\\\n",
"-1 &1 &-1 &1 \\\\\n",
"1 &-1 &-1 &1 \n",
"\\end{bmatrix}\n",
"\\begin{bmatrix} \n",
"11 &0 &0 &0 \\\\\n",
"0 &7 &0 &0 \\\\\n",
"0 &0 &5 &0 \\\\\n",
"0 &0 &0 &1 \n",
"\\end{bmatrix}\n",
"\\frac{1}{2}\n",
"\\begin{bmatrix} \n",
"-1 &-1 &-1 &-1 \\\\\n",
"-1 &-1 &1 &1 \\\\\n",
"-1 &1 &1 &-1 \\\\\n",
"1 &-1 &1 &-1 \n",
"\\end{bmatrix}. \\tag{2}\n",
"$$\n",
"\n",
"Import packages."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.008567Z",
"start_time": "2021-03-09T03:44:29.796997Z"
}
},
"outputs": [],
"source": [
"import time\n",
"import numpy as np\n",
"from numpy import pi as PI\n",
"from matplotlib import pyplot as plt\n",
"from scipy.stats import unitary_group\n",
"from scipy.linalg import norm\n",
"\n",
"import paddle\n",
"from paddle import matmul, transpose, trace\n",
"from paddle_quantum.circuit import *\n",
"from paddle_quantum.utils import *\n",
"\n",
"\n",
"\n",
"# Draw the learning curve in the optimization process\n",
"def loss_plot(loss):\n",
" '''\n",
" loss is a list, this function plots loss over iteration\n",
" '''\n",
" plt.plot(list(range(1, len(loss)+1)), loss)\n",
" plt.xlabel('iteration')\n",
" plt.ylabel('loss')\n",
" plt.title('Loss Over Iteration')\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Classical Singular Value Decomposition\n",
"\n",
"With the above mathematical definition, one can realize SVD numerically through NumPy."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.056721Z",
"start_time": "2021-03-09T03:44:34.012222Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix M we want to decompose is: \n",
"[[ 3.+0.j 6.+0.j 2.+0.j 0.+0.j]\n",
" [ 6.+0.j 3.+0.j 0.+0.j -2.+0.j]\n",
" [ 2.+0.j 0.+0.j 3.+0.j -6.+0.j]\n",
" [ 0.+0.j -2.+0.j -6.+0.j 3.+0.j]]\n"
]
}
],
"source": [
"# Generate matrix M\n",
"def M_generator():\n",
" I = np.array([[1, 0], [0, 1]])\n",
" Z = np.array([[1, 0], [0, -1]])\n",
" X = np.array([[0, 1], [1, 0]])\n",
" Y = np.array([[0, -1j], [1j, 0]])\n",
" M = 2 *np.kron(X, Z) + 6 * np.kron(Z, X) + 3 * np.kron(I, I)\n",
" return M.astype('complex64')\n",
"\n",
"print('The matrix M we want to decompose is: ')\n",
"print(M_generator())"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.093725Z",
"start_time": "2021-03-09T03:44:34.063353Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The singular values of the matrix from large to small are:\n",
"[11. 7. 5. 1.]\n",
"The decomposed unitary matrix U is:\n",
"[[-0.5+0.j -0.5+0.j 0.5+0.j 0.5+0.j]\n",
" [-0.5+0.j -0.5+0.j -0.5+0.j -0.5+0.j]\n",
" [-0.5+0.j 0.5+0.j -0.5+0.j 0.5+0.j]\n",
" [ 0.5+0.j -0.5+0.j -0.5+0.j 0.5+0.j]]\n",
"The decomposed unitary matrix V_dagger is:\n",
"[[-0.5+0.j -0.5+0.j -0.5+0.j 0.5+0.j]\n",
" [-0.5+0.j -0.5+0.j 0.5+0.j -0.5+0.j]\n",
" [-0.5+0.j 0.5+0.j 0.5+0.j 0.5+0.j]\n",
" [-0.5+0.j 0.5+0.j -0.5+0.j -0.5+0.j]]\n"
]
}
],
"source": [
"# We only need the following line of code to complete SVD\n",
"U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)\n",
"\n",
"\n",
"# Print decomposition results\n",
"print(\"The singular values of the matrix from large to small are:\")\n",
"print(D)\n",
"print(\"The decomposed unitary matrix U is:\")\n",
"print(U)\n",
"print(\"The decomposed unitary matrix V_dagger is:\")\n",
"print(V_dagger)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.112670Z",
"start_time": "2021-03-09T03:44:34.098847Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 3.+0.j 6.+0.j 2.+0.j 0.+0.j]\n",
" [ 6.+0.j 3.+0.j 0.+0.j -2.+0.j]\n",
" [ 2.+0.j 0.+0.j 3.+0.j -6.+0.j]\n",
" [ 0.+0.j -2.+0.j -6.+0.j 3.+0.j]]\n"
]
}
],
"source": [
"# Then assemble it back, can we restore the original matrix?\n",
"M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))\n",
"print(M_reconst)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Surely, we can be restored the original matrix $M$! One can further modify the matrix, see what happens if it is not a square matrix.\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Quantum Singular Value Decomposition\n",
"\n",
"Next, let's take a look at what the quantum version of singular value decomposition is all about. In summary, we transform the problem of matrix factorization into an optimization problem with the variational principle of singular values. Specifically, this is achieved through the following four steps:\n",
"\n",
"- Prepare an orthonormal basis $\\{|\\psi_j\\rangle\\}$, one can take the computational basis $\\{ |000\\rangle, |001\\rangle,\\cdots |111\\rangle\\}$ (this is in the case of 3 qubits)\n",
"- Prepare two parameterized quantum neural networks $U(\\theta)$ and $V(\\phi)$ to learn left/right singular vectors respectively\n",
"- Use quantum neural network to estimate singular values $m_j = \\text{Re}\\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle$\n",
"- Design the loss function $\\mathcal{L}(\\theta)$ and use PaddlePaddle Deep Learning framework to maximize the following quantity, \n",
"\n",
"$$\n",
"L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle\\psi_j|U(\\theta)^{\\dagger} MV(\\phi)|\\psi_j\\rangle. \\tag{3}\n",
"$$\n",
"\n",
"Where $q_1>\\cdots>q_T>0$ is the adjustable weights (hyperparameter), and $T$ represents the rank we want to learn or the total number of singular values to be learned.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Case 1: Decompose a randomly generated $8\\times8$ complex matrix\n",
"\n",
"Then we look at a specific example, which can better explain the overall process."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.132465Z",
"start_time": "2021-03-09T03:44:34.116446Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix M we want to decompose is:\n",
"[[6.+1.j 3.+9.j 7.+3.j 4.+7.j 6.+6.j 9.+8.j 2.+7.j 6.+4.j]\n",
" [7.+1.j 4.+4.j 3.+7.j 7.+9.j 7.+8.j 2.+8.j 5.+0.j 4.+8.j]\n",
" [1.+6.j 7.+8.j 5.+7.j 1.+0.j 4.+7.j 0.+7.j 9.+2.j 5.+0.j]\n",
" [8.+7.j 0.+2.j 9.+2.j 2.+0.j 6.+4.j 3.+9.j 8.+6.j 2.+9.j]\n",
" [4.+8.j 2.+6.j 6.+8.j 4.+7.j 8.+1.j 6.+0.j 1.+6.j 3.+6.j]\n",
" [8.+7.j 1.+4.j 9.+2.j 8.+7.j 9.+5.j 4.+2.j 1.+0.j 3.+2.j]\n",
" [6.+4.j 7.+2.j 2.+0.j 0.+4.j 3.+9.j 1.+6.j 7.+6.j 3.+8.j]\n",
" [1.+9.j 5.+9.j 5.+2.j 9.+6.j 3.+0.j 5.+3.j 1.+3.j 9.+4.j]]\n",
"The singular values of the matrix M are:\n",
"[54.83484985 19.18141073 14.98866247 11.61419557 10.15927045 7.60223249\n",
" 5.81040539 3.30116001]\n"
]
}
],
"source": [
"# First fix the random seed, in order to reproduce the results at any time\n",
"np.random.seed(42)\n",
"\n",
"# Set the number of qubits, which determines the dimension of the Hilbert space\n",
"N = 3\n",
"\n",
"# Make a random matrix generator\n",
"def random_M_generator():\n",
" M = np.random.randint(10, size = (2**N, 2**N)) + 1j*np.random.randint(10, size = (2**N, 2**N))\n",
" return M\n",
"\n",
"M = random_M_generator()\n",
"M_err = np.copy(M)\n",
"\n",
"\n",
"# Output the matrix M\n",
"print('The matrix M we want to decompose is:')\n",
"print(M)\n",
"\n",
"# Apply SVD and record the exact singular values\n",
"U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
"print(\"The singular values of the matrix M are:\")\n",
"print(D)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.147570Z",
"start_time": "2021-03-09T03:44:34.138265Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The selected weight is:\n",
"[24.+0.j 21.+0.j 18.+0.j 15.+0.j 12.+0.j 9.+0.j 6.+0.j 3.+0.j]\n"
]
}
],
"source": [
"# Hyperparameter settings\n",
"N = 3 # Number of qubits\n",
"T = 8 # Set the number of rank you want to learn\n",
"ITR = 100 # Number of iterations\n",
"LR = 0.02 # Learning rate\n",
"SEED = 14 # Random seed\n",
"\n",
"# Set the learning weight \n",
"weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
"print('The selected weight is:')\n",
"print(weight)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We design QNN with the following structure:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:44:34.245692Z",
"start_time": "2021-03-09T03:44:34.226859Z"
}
},
"outputs": [],
"source": [
"# Set circuit parameters\n",
"cir_depth = 20 # circuit depth\n",
"block_len = 2 # length of each block\n",
"theta_size = N * block_len * cir_depth # size of the network parameter theta\n",
"\n",
"\n",
"# Define quantum neural network\n",
"def U_theta(theta):\n",
"\n",
" # Initialize the network with UAnsatz\n",
" cir = UAnsatz(N)\n",
" \n",
" # Build QNN\n",
" for layer_num in range(cir_depth):\n",
" \n",
" for which_qubit in range(N):\n",
" cir.ry(theta[block_len * layer_num * N + which_qubit], which_qubit)\n",
" \n",
" for which_qubit in range(N):\n",
" cir.rz(theta[(block_len * layer_num + 1) * N + which_qubit], which_qubit)\n",
"\n",
" for which_qubit in range(1, N):\n",
" cir.cnot([which_qubit - 1, which_qubit])\n",
" cir.cnot([N - 1, 0])\n",
"\n",
" return cir.U"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we complete the main part of the algorithm:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:46:12.944634Z",
"start_time": "2021-03-09T03:44:50.626213Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter: 0 loss: -235.9752\n",
"iter: 10 loss: -1635.7805\n",
"iter: 20 loss: -1967.8480\n",
"iter: 30 loss: -2154.2274\n",
"iter: 40 loss: -2257.8776\n",
"iter: 50 loss: -2311.8162\n",
"iter: 60 loss: -2342.3099\n",
"iter: 70 loss: -2358.9515\n",
"iter: 80 loss: -2369.3172\n",
"iter: 90 loss: -2375.9682\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnHklEQVR4nO3deXxU9b3/8ddnZrInJEACgYAssiigYkUFt1q1BVuvdLPaW61ttba9tZvdsPXe295b78+23u6LtbW2dtF6tYttrbjUBW2tLCqySEFABAmELQFC9s/vj3OCQ0wwZGZykpn38/GYx8z5njNzPicH8s75fuecY+6OiIhIKmJRFyAiIoOfwkRERFKmMBERkZQpTEREJGUKExERSZnCREREUqYwEZEjYmb7zGxi1HXIwKIwkUHHzDaa2XkRrfs0M/urme01s3oz+6OZTevH9R/cdjN7n5k9nuH1PWJmVya3uXupu6/P5Hpl8FGYiPSSmc0B7gf+AIwGJgDPAk+k+y91C2T0/6eZJTL5+ZJbFCaSNcyswMy+ZWYvh49vmVlBOK/SzP5kZnvMbJeZLer8ZW1mnzezLeHRxhozO7eHVXwNuM3dv+3ue919l7tfBzwJfCn8rNVmdkFSTQkzqzOz14XTs83sb2Edz5rZ2UnLPmJm15vZE0Aj0GNAmdmxwE3AnLDbaU/Sz+BGM9tkZtvM7CYzKwrnnW1mm8PtrQVuNbOh4c+lzsx2h6/HhMtfD5wJfC9cx/fCdjezSeHrcjO7LXz/i2Z2XdLP9X1m9nhYz24z22Bm5/d6h8qgojCRbPJFYDYwEzgBOAW4Lpz3aWAzUAWMBL4AuJlNBa4GTnb3MmAusLHrB5tZMXAa8H/drPdO4I3h69uBdyfNmwvscPdlZlYD/Bn4CjAM+Axwt5lVJS1/GXAVUAa82NOGuvtq4MPA38Nup4pw1g3AlPBnMAmoAf4j6a3V4brHheuJAbeG00cBB4Dvhev4IrAIuDpcx9XdlPJdoJwg+F4PvBd4f9L8U4E1QCVBGN9iZtbTdsngpTCRbPIe4L/cfbu71wFfJvjlDNAKjALGuXuruy/y4MJ07UABMM3M8tx9o7u/0M1nDyP4/7K1m3lbCX5ZAvwauDAMH4B/JQgYgEuBe939XnfvcPcHgCXAm5M+62fuvtLd29y99Ug2PvwlfRXwqfCoaS/wP8AlSYt1AP/p7s3ufsDdd7r73e7eGC5/PUEo9GZ98fCzrw2P1DYC/8srP3OAF939x+7eDvycYB+MPJLtksFBYSLZZDSH/jX/YtgG8HVgHXC/ma03swUA7r4O+CRBN9V2M7vDzEbzarsJfhGP6mbeKGBH0uetBv4lDJQLCQIGgr/+Lwq7uPaEXVNndPnMl45kg7uoAoqBpUmff1/Y3qnO3Zs6J8ys2Mx+FHZRNQCPARVhULyWSiCPV//Ma5KmaztfuHtj+LL0CLZJBgmFiWSTlwl+YXc6Kmwj/Mv50+4+keAX/DWdYyPu/mt3PyN8rwNf7frB7r4f+DtwUTfrfRfwUNJ0Z1fXfGBVGDAQBMUv3L0i6VHi7jckr+oItrfrsjsIuqmmJ31+ubuXHuY9nwamAqe6+xDgrLDdeli+6/paefXPfMsRbINkCYWJDFZ5ZlaY9EgQ/BK/zsyqzKySYKzglwBmdoGZTQq7guoJurc6zGyqmZ0TDtQ3Efwy7uhhnQuAy83s42ZWFg5efwWYQ9Cl1ukO4E3AR3jlqISwln8xs7lmFg/rPrtzwLsPtgFjzCwfwN07gB8D3zSzEeF215jZ3MN8RhnBNu8xs2HAf3azjm6/CBB2Xd0JXB/+PMYB14TbKTlGYSKD1b0EvwQ7H18iGNheAiwHngOWhW0Ak4EHgX0ERxg/cPeHCcZLbiD4K7sWGAFc290K3f1xggH1txOMk7wInAic4e5rk5bbGq7jNOA3Se0vERytfAGoIzhS+Sx9/3/4V2AlUGtmO8K2zxN05z0Zdls9SHDk0ZNvAUUE2/8kQbdYsm8D7wy/jfWdbt7/MWA/sB54nCA8f9qnrZFBzXRzLBERSZWOTEREJGUKExERSZnCREREUqYwERGRlOXshd4qKyt9/PjxUZchIjKoLF26dIe7V3Vtz9kwGT9+PEuWLIm6DBGRQcXMur1mnLq5REQkZQoTERFJmcJERERSpjAREZGUKUxERCRlChMREUmZwkRERFKmMDlCv396C798ssdbc4uI5CSFyRG697mt3Pb3jVGXISIyoChMjtDoiiK27D6A7gMjIvIKhckRGjO0iP0t7TQ0tUVdiojIgKEwOUKjK4oA2LL7QMSViIgMHAqTI9QZJi/vUZiIiHRSmByhms4wqVeYiIh0UpgcoeEl+eQnYurmEhFJojA5QrGYMbq8kC3q5hIROUhh0gc1Q4sUJiIiSRQmfTC6vEgD8CIiSRQmfTC6oojte5tpaeuIuhQRkQFBYdIHNRVFuENtfVPUpYiIDAgKkz6oGRqeuKiuLhERQGHSJzpxUUTkUAqTPhhVXgjoyEREpJPCpA8K8+JUlhboyEREJKQw6aOaCp24KCLSSWHSR6MrdOKiiEgnhUkf1VQEJy7qJlkiIgqTPhtdUURTawe7G1ujLkVEJHIKkz7STbJERF4x4MLEzL5kZlvM7Jnw8eakedea2TozW2Nmc5Pa54Vt68xsQX/UOUYnLoqIHJSIuoAefNPdb0xuMLNpwCXAdGA08KCZTQlnfx94I7AZWGxm97j7qkwWqBMXRUReMVDDpDvzgTvcvRnYYGbrgFPCeevcfT2Amd0RLpvRMBlanEdhXkxHJiIiDMBurtDVZrbczH5qZkPDthrgpaRlNodtPbVnlJkxukKXohcRgYjCxMweNLMV3TzmAz8EjgZmAluB/03jeq8ysyVmtqSuri7lz6tRmIiIABF1c7n7eb1Zzsx+DPwpnNwCjE2aPSZs4zDtXdd7M3AzwKxZs1I+QaSmoojVWxtS/RgRkUFvwHVzmdmopMm3ASvC1/cAl5hZgZlNACYDTwGLgclmNsHM8gkG6e/pj1pHVxSxY18LTa3t/bE6EZEBayAOwH/NzGYCDmwEPgTg7ivN7E6CgfU24KPu3g5gZlcDC4E48FN3X9kfhXZ+o2trfRMTKkv6Y5UiIgPSgAsTd7/sMPOuB67vpv1e4N5M1tWd0RXBpehf3nNAYSIiOW3AdXMNJjUVOnFRRAQUJimpLn/lyEREJJcpTFJQkIhTVaabZImIKExSFJxr0hR1GSIikVKYpEgnLoqIKExSNjq8fa9ukiUiuUxhkqLRFUU0t3Wwa39L1KWIiERGYZKiVy5Fr3ETEcldCpMU6VwTERGFScp0kywREYVJyjpvkqUwEZFcpjBJ0cGbZNUrTEQkdylM0qCmoogtGoAXkRymMEmD0eU6cVFEcpvCJA1GVxRRt7eZ5jbdJEtEcpPCJA0672tSW6+uLhHJTQqTNNC5JiKS6xQmaaCz4EUk1ylM0kA3yRKRXKcwSYPCvDiVpbpJlojkLoVJmtSEl6IXEclFCpM0qRmqc01EJHcpTNKkpqKIzbsP0N6hm2SJSO5RmKTJ5JFlNLd1sGlXY9SliIj0O4VJmhxTXQbAmtqGiCsREel/CpM0mTyiDDNYU7sv6lJERPqdwiRNivLjjBtWzJptOjIRkdyjMEmjKSPLeL52b9RliIj0O4VJGh1TXcbGHftpatXVg0UktyhM0mhq9RA6HNZt17iJiOQWhUkaTa0uBeCf29TVJSK5RWGSRuOHl5Afj7FG4yYikmMUJmmUiMc4ekSpBuFFJOcoTNLsmOoydXOJSM5RmKTZ1OoyttY3Ud/YGnUpIiL9RmGSZlNHhpdV0dGJiOQQhUmaTa1WmIhI7lGYpNmo8kLKChO64KOI5JRIwsTMLjKzlWbWYWazusy71szWmdkaM5ub1D4vbFtnZguS2ieY2T/C9t+YWX5/bktXZsbUkWX8Uxd8FJEcEtWRyQrg7cBjyY1mNg24BJgOzAN+YGZxM4sD3wfOB6YB7w6XBfgq8E13nwTsBq7on03o2dTqMp6vbcBdN8oSkdwQSZi4+2p3X9PNrPnAHe7e7O4bgHXAKeFjnbuvd/cW4A5gvpkZcA5wV/j+nwNvzfgGvIaZYytoaGpjxRZ1dYlIbhhoYyY1wEtJ05vDtp7ahwN73L2tS3u3zOwqM1tiZkvq6urSWniyc48dScxg4crajK1DRGQgyViYmNmDZraim8f8TK3ztbj7ze4+y91nVVVVZWw9w0ryOWXCMIWJiOSMRKY+2N3P68PbtgBjk6bHhG300L4TqDCzRHh0krx8pOZOr+bLf1zF+rp9TKwqjbocEZGMGmjdXPcAl5hZgZlNACYDTwGLgcnhN7fyCQbp7/FghPth4J3h+y8H/hBB3a/ypunVACxcuS3iSkREMi+qrwa/zcw2A3OAP5vZQgB3XwncCawC7gM+6u7t4VHH1cBCYDVwZ7gswOeBa8xsHcEYyi39uzXdq6ko4riacnV1iUhOyFg31+G4+++A3/Uw73rg+m7a7wXu7aZ9PcG3vQacudNHcuP9/6S2vonq8sKoyxERyZiB1s2VVeaGXV0PrNLRiYhkN4VJBk0aUcrEyhKNm4hI1lOYZJCZ8abp1Ty5fid7m3RJehHJXgqTDJtz9HDaOpznNtdHXYqISMYoTDLs+JpyAJ5VmIhIFlOYZNjQknyOGlbM8s17oi5FRCRjFCb94Lgx5SzXkYmIZDGFST84YUw5W/YcYOe+5qhLERHJCIVJPzh+TAWAjk5EJGspTPrBjJpyzBQmIpK9FCb9oLQgwdFVpRqEF5GspTDpJ8ePKefZzfW6la+IZCWFST85vqacHfuaqW1oiroUEZG0U5j0k+PHVgDw7EsaNxGR7KMw6SfTRg0hETONm4hIVlKY9JPCvDhTRpbx3BYdmYhI9ulVmJjZJ8xsiAVuMbNlZvamTBeXbU4YG5wJr0F4Eck2vT0y+YC7NwBvAoYClwE3ZKyqLHX8mArqD7SyamtD1KWIiKRVb8PEwuc3A78I779uh1leujFvejWlBQm+//C6qEsREUmr3obJUjO7nyBMFppZGdCRubKy09CSfD5w+njufa6WVS/r6EREskdvw+QKYAFwsrs3AnnA+zNWVRa74oyJlBUm+OaD/4y6FBGRtOltmMwB1rj7HjO7FLgO0NeS+qC8OI8PnjmRB1Zt090XRSRr9DZMfgg0mtkJwKeBF4DbMlZVlnv/6eOpKM7jGw+siboUEZG06G2YtHnwfdb5wPfc/ftAWebKym5lhXlcddZEHl5Tp6MTEckKvQ2TvWZ2LcFXgv9sZjGCcRPpo0tnj6MwL8btizdFXYqISMp6GyYXA80E55vUAmOAr2esqhwwpDCPNx83inueeZnGlraoyxERSUmvwiQMkF8B5WZ2AdDk7hozSdHFs8ayr7mNe5+rjboUEZGU9PZyKu8CngIuAt4F/MPM3pnJwnLBKROGMaGyhDsXvxR1KSIiKeltN9cXCc4xudzd3wucAvx75srKDWbGu2aN5amNu3ihbl/U5YiI9FlvwyTm7tuTpncewXvlMN5xUg3xmHHnEh2diMjg1dtAuM/MFprZ+8zsfcCfgXszV1buGFFWyBumjuDupVtobdcVakRkcOrtAPxngZuB48PHze7++UwWlksuPnksO/Y18+iauqhLERHpk0RvF3T3u4G7M1hLznr9lCrKChPct7KW86aNjLocEZEjdtgwMbO9QHd3cjLA3X1IRqrKMfmJGOceM4IHV2+jrb2DRFzDUSIyuBz2t5a7l7n7kG4eZQqS9Jo3o5o9ja08tWFX1KWIiBwx/Qk8QJw1pYqCRIz7VuoERhEZfBQmA0RxfoLXT6ni/pXb6OjQPeJFZHCJJEzM7CIzW2lmHWY2K6l9vJkdMLNnwsdNSfNOMrPnzGydmX3HzCxsH2ZmD5jZ2vB5aBTblA5zp1dT29DEs5v3RF2KiMgRierIZAXwduCxbua94O4zw8eHk9p/CHwQmBw+5oXtC4CH3H0y8FA4PSide+wI4jFj4cptUZciInJEIgkTd1/t7r2+M5SZjQKGuPuT4X1VbgPeGs6eD/w8fP3zpPZBp6I4nzkTh7NwZS3BZoqIDA4Dccxkgpk9bWaPmtmZYVsNsDlpmc1hG8BId98avq4FejxRw8yuMrMlZrakrm5gniA4d/pINuzYz9rtulaXiAweGQsTM3vQzFZ085h/mLdtBY5y9xOBa4Bfm1mvv4IcHrX0+Ce9u9/s7rPcfVZVVVWvt6U/zZ1eTX48xv/cu1oD8SIyaGQsTNz9PHef0c3jD4d5T7O77wxfLyW41/wUYAvBDbk6jQnbALaF3WCd3WHJF6QcdEYMKeTfLziWR9bU8aPH1kddjohIrwyobi4zqzKzePh6IsFA+/qwG6vBzGaH3+J6L9AZSvcAl4evL09qH7QunT2Otxw3ihvvX8OSjTqJUUQGvqi+Gvw2M9sMzCG4p/zCcNZZwHIzewa4C/iwu3f+Nv034CfAOoIjlr+E7TcAbzSztcB54fSgZmb8v3ccx5ihRXzs9qfZtb8l6pJERA7LcvVbQ7NmzfIlS5ZEXcZhrdhSz9t+8ATvOXUcX7pwetTliIhgZkvdfVbX9gHVzSWHmlFTzr+cMJo7l7xEfWNr1OWIiPRIYTLAXXnGRBpb2vn1U5uiLkVEpEcKkwFu2ughnD5pOD/72wZa2nQnRhEZmBQmg8CVZ0xkW0Mz9z639bUXFhGJgMJkEHj9lCqOrirhx4vW6zIrIjIgKUwGgVjMuOKMiax8uYG/r98ZdTkiIq+iMBkk3v66GoaX5PP9h9dFXYqIyKsoTAaJwrw4Hzn7aJ5Yt5O/vbAj6nJERA6hMBlELp09juohhdy4cI3GTkRkQFGYDCKFeXGuPmcSyzbt4ZE1A/MS+iKSmxQmg8y7Zo1l7LAibrx/jS5RLyIDhsJkkMlPxPjkuVNY+XIDf1z+ctTliIgACpNB6a0n1jCjZgjX3Pkstzy+QeMnIhI5hckgFI8Zv/7gbM45ZgT//adVfOKOZ2hsaYu6LBHJYQqTQWpIYR4/uvQkPjt3Kn9c/jIf+eUyHaGISGQUJoNYLGZ89A2TuO4t03j0n3X8ZUVt1CWJSI5SmGSBy+eM49hRQ/ivP65if7O6u0Sk/ylMskAiHuMrb51ObUMT3/nr2qjLEZEcpDDJEieNG8ZFJ43hlkUbWLttb9TliEiOUZhkkQXnH0NJQYLP3LWcAy3tUZcjIjlEYZJFhpcW8NV3HM/yzXv42O1P09auOzOKSP9QmGSZeTOq+fKF03lw9Tb+/Q8r9XVhEekXiagLkPR775zx1NY38YNHXqC8KI/Pzp1KPGZRlyUiWUxhkqU+O3cquxtbuOnRF1i8cRc3XnQCEypLoi5LRLKUurmylJnxP287jm9dPJO12/Zy/rcf4/anNkVdlohkKYVJFjMz3npiDfd/6vWcPH4Y1/72OR5YtS3qskQkCylMckB1eSE/fu+s8ErDz7BpZ2PUJYlIllGY5IjCvDg/fM9JxMz48C+X0tSq81BEJH0UJjlk7LBivnXxTFZtbWDB3csVKCKSNgqTHPOGY0ZwzRun8PtnXmbutx7j4TXboy5JRLKAwiQHffzcyfziilOIm/H+WxfzwduWaBxFRFKiMMlRZ06u4i+fPJPPzZvKE+t2cN43HuWGvzzPPl3CXkT6QGGSwwoScf7t7Ek8/JmzueCEUdz06Aucc+MjPLe5PurSRGSQUZgII4cU8o13zeT3Hz2dvHiMd//4Sf7+ws6oyxKRQURhIgfNHFvB3R85jVHlhVx+61Pcv1K3ARaR3lGYyCGqywu580NzOHbUED70y6V84GeLeXDVNl3OXkQOS2EirzK0JJ9fX3kqV79hEiu21HPlbUs482sP86AuxSIiPYgkTMzs62b2vJktN7PfmVlF0rxrzWydma0xs7lJ7fPCtnVmtiCpfYKZ/SNs/42Z5ffz5mSlkoIEn37TVJ5YcA4/uuwkhhbnc+VtS/jvP62ipU1HKSJyqKiOTB4AZrj78cA/gWsBzGwacAkwHZgH/MDM4mYWB74PnA9MA94dLgvwVeCb7j4J2A1c0a9bkuXy4jHmTq/mdx89jfedNp5bHt/AO2/6m85LEZFDRBIm7n6/u3ee0PAkMCZ8PR+4w92b3X0DsA44JXysc/f17t4C3AHMNzMDzgHuCt//c+Ct/bQZOaUgEedLF07npktPYuOO/bzlu4s0QC8iBw2EMZMPAH8JX9cALyXN2xy29dQ+HNiTFEyd7d0ys6vMbImZLamrq0tT+bll3oxq/vzxM5lQWcJVv1jKV/60iuY2XeNLJNdlLEzM7EEzW9HNY37SMl8E2oBfZaqOZO5+s7vPcvdZVVVV/bHKrDR2WDH/9+E5XD5nHD95fAOn3/BXvnH/GrY3NEVdmohEJGO37XX38w4338zeB1wAnOvuHjZvAcYmLTYmbKOH9p1AhZklwqOT5OUlgwoScb48fwZzZ1Rzy6INfPfhdfzgkRf44FkTueaNU8iLD4SDXhHpL5HcA97M5gGfA17v7skjufcAvzazbwCjgcnAU4ABk81sAkFYXAL8q7u7mT0MvJNgHOVy4A/9tyVy2tGVnHZ0JS/u3M93/7qOHz7yAk9t2MV3330ioyuKoi5PRPpJVH8+fg8oAx4ws2fM7CYAd18J3AmsAu4DPuru7eFRx9XAQmA1cGe4LMDngWvMbB3BGMot/bspAjBueAk3XnQC33n3iTy/tYE3f2cRv1m8SfdMEckR9koPU26ZNWuWL1myJOoystKGHfv55B1P8+zmeipL83nPqeO4bM44KksLoi5NRFJkZkvdfdar2hUmkgnuzt9e2MlPH9/AQ89vJz8R4x2vG8OVZ07g6KrSqMsTkT7qKUwiGTOR7GdmnD6pktMnVfJC3T5ueXwDdy3dzO1PbeKC40dx3VumUV1eGHWZIpImOjKRfrNjXzO3PrGBnyzaQCJmfOqNU3jfaeNJ6JtfIoOGurm6UJhEZ9PORv7jnhU8sqaOo4YV886TxvD219UwZmhx1KWJyGtQmHShMImWu3P/qm3c+sQGnly/C4Czp1Zx3VuOZdKIsoirE5GeKEy6UJgMHC/tauS3y7Zwy+PraWxp54ozJ/DxcyZTUqAhPZGBRmHShcJk4Nmxr5mv3fc8dy7ZTGVpPpfNHs97Zh+lrxSLDCAKky4UJgPXsk27+e5Da3l4TR35iRgXnjCay2aP44SxFVGXJpLzFCZdKEwGvhfq9nHrExv47bItNLa0c1xNOZfNHseFM0dTmBePujyRnKQw6UJhMnjsbWrl909v4ZdPbmLNtr0Hu8AumzOOYSW6saZIf1KYdKEwGXzcnb+v38lPFm3gr89vJz8e46RxQzljciVnTKrk+DHlBPdLE5FMUZh0oTAZ3NZu28v/Ld3MorU7WL21AYCjq0q45OSjeMdJY3TEIpIhCpMuFCbZo25vMw8/v507Fm9i2aY95MWNmWMrOHXCcE6dOIyTxw/TGItImihMulCYZKc1tXv53dNb+Pv6nazYUk97h1OcH+esyVWcN20kr59SRVWZvmos0le60KPkhKnVZSw4/xgA9jW3sXjjLh5avY0HV23nvpW1ABw7aghnTa5k9tHDmTVuKGWFeVGWLJIVdGQiOcHdWbGlgcfW1rFobR1LX9xNa7sTM5hRU87pkyo579gRzBw7lHhMg/giPVE3VxcKk9zW2NLGshf38NSGnTy5fhdLN+2mvcMZXpLPnKOHc+JRQ5k5toIZNUMoSGi8RaSTurlEkhTnJ4KvFE+uBKC+sZVH19bx0OptLN6wiz8t3wpAYV6MUycM56wpVZw6YRiTRpRqMF+kGzoyEenG9oYmlm3aw5Prd/LY2jrW1+0HIGZw1LBiplaXcfyYCo6rKWfa6CEML8nXOS6SE3RkInIERgwpZN6MaubNqAaCKxs/u3kPa7ftY+32vax6uYGFK7cdXL60IMGYoUWMG17MpBGlTB5RxqQRpRw1vJghGuCXHKAwEemFscOKGTvs0Jt31R9oZeWWelbX7uWlXY28tKuRtdv38eDq7bR3vHLEP6wkn3HDi5lQWcLEyhImVpUybngx44aXUKrL7EuW0L9kkT4qL8rjtEmVnDap8pD2lrYOXty5n3Xb9/HirkZe3NnIxh37+du6nfx22ZZDlq0szWfM0GLGDC1Keg4eo8qLdE8XGTT0L1UkzfITMSaPLGPyyFffMXJ/cxsbd+7nxZ2N4WM/m3cfYMWWehaurKW1/dAxzLLCBKPKCxlVXsToiiJqKgqpLi9i5JACRg4pZGRZIUOKEhqvkcgpTET6UUlBgumjy5k+uvxV89o7nLq9zWze3cjm3QfYWt9EbX3wvLW+iRVb6tm5v+VV7yvKi1NdXkhVWQFVpQUML82nsrSAytICqsoKqCzNZ1hJPkNL8ikrUPBIZihMRAaIeMyoLi+kuryQWeO7X6aptZ1tDU1sa2imtqGJ7Q1N1NY3sbWhibq9zayubWDnvhbqD7R2+/5EzBhWEoTL8NJ8KoryKS/OY2hx3sHX5UV5DCnMo6wwcfC5rDBBIh7L3MbLoKcwERlECvPijBtewrjhJYddrqWtg537m6nb28zOfS3s2t/C7sYWdu5vYde+8Hl/M7X1DexpbGXPgdZDvjTQnaK8+MFgKS3MY0hhgtKC4FFSELaHr0sK4hTnJyjJT1BcEA+e8+MU5sUpzo9TlBcnpisNZBWFiUgWyk/EGFUeDOL3hruzr7mNPY2t1B9opaGplYYDbextamVvU1v4aGVfc/C6oamV/c1tbGtoYm9TG/ua29jf3MZr5NEhChIxivOD0CnMi1EUhkxh3ivPhXmx8DlOYSJGQV6cgkTQVpA0HTziFOTFDk7nx+PkJ2Lkd04nYiRipm6+DFGYiAhmRllhHmWFeYzt42e4O40t7eEjCJjGlnb2N7exvzloa2pt50BrsMyBcNn9LW00t3YcnLevuY26vc0caG0P2tvaaWptp6m1I+XtjBnkxYNgye98Dl8nt+cljLzOtniMvLiRCJ/z4jESsaTX8c5ljXjYnoh1todtMSMeMxI9TCc6pw8+x4jHjbi90h5Lmh+PBfMG0tGdwkRE0sLMwi6uBJD+y/y7O81tHTS3dtDcHgRNc1sQMkF7O83tHbS0BdMtBx/ttLSH72vroLU9nN/eQWv43BK2t7Q7rW0dNLV2sLep7WB7W4eHyzptHR20tTst7R20tXcc0dFYuplxMFQ6gydmkIjHiJkRjyXNTwqgn15+MkcNL37tFRwBhYmIDApmdrDLCwbOVQXaO/xg4LS3O61h2Bxs6wie29r94HRru9PREUx3hlOHe9Dunct20N7Bwfe3d4TzwvW0e9DW3hG87uhw2jugI2xv6wjb/JX3dj4X5KX/yxQKExGRFATdTrr4p77rJyIiKVOYiIhIyhQmIiKSMoWJiIikTGEiIiIpU5iIiEjKFCYiIpIyhYmIiKTM3CO8FkCEzKwOePEI3lIJ7MhQOQNVLm4z5OZ25+I2Q25ud6rbPM7dq7o25myYHCkzW+Lus6Kuoz/l4jZDbm53Lm4z5OZ2Z2qb1c0lIiIpU5iIiEjKFCa9d3PUBUQgF7cZcnO7c3GbITe3OyPbrDETERFJmY5MREQkZQoTERFJmcLkNZjZPDNbY2brzGxB1PVkipmNNbOHzWyVma00s0+E7cPM7AEzWxs+D4261nQzs7iZPW1mfwqnJ5jZP8J9/hszy4+6xnQzswozu8vMnjez1WY2J9v3tZl9Kvy3vcLMbjezwmzc12b2UzPbbmYrktq63bcW+E64/cvN7HV9Xa/C5DDMLA58HzgfmAa828ymRVtVxrQBn3b3acBs4KPhti4AHnL3ycBD4XS2+QSwOmn6q8A33X0SsBu4IpKqMuvbwH3ufgxwAsH2Z+2+NrMa4OPALHefAcSBS8jOff0zYF6Xtp727fnA5PBxFfDDvq5UYXJ4pwDr3H29u7cAdwDzI64pI9x9q7svC1/vJfjlUkOwvT8PF/s58NZICswQMxsDvAX4SThtwDnAXeEi2bjN5cBZwC0A7t7i7nvI8n1NcJvyIjNLAMXAVrJwX7v7Y8CuLs097dv5wG0eeBKoMLNRfVmvwuTwaoCXkqY3h21ZzczGAycC/wBGuvvWcFYtMDKqujLkW8DngI5wejiwx93bwuls3OcTgDrg1rB77ydmVkIW72t33wLcCGwiCJF6YCnZv6879bRv0/Y7TmEihzCzUuBu4JPu3pA8z4PvkWfNd8nN7AJgu7svjbqWfpYAXgf80N1PBPbTpUsrC/f1UIK/wicAo4ESXt0VlBMytW8VJoe3BRibND0mbMtKZpZHECS/cvffhs3bOg97w+ftUdWXAacDF5rZRoIuzHMIxhIqwq4QyM59vhnY7O7/CKfvIgiXbN7X5wEb3L3O3VuB3xLs/2zf15162rdp+x2nMDm8xcDk8Bsf+QQDdvdEXFNGhGMFtwCr3f0bSbPuAS4PX18O/KG/a8sUd7/W3ce4+3iCfftXd38P8DDwznCxrNpmAHevBV4ys6lh07nAKrJ4XxN0b802s+Lw33rnNmf1vk7S0769B3hv+K2u2UB9UnfYEdEZ8K/BzN5M0K8eB37q7tdHW1FmmNkZwCLgOV4ZP/gCwbjJncBRBJfsf5e7dx3cG/TM7GzgM+5+gZlNJDhSGQY8DVzq7s0Rlpd2ZjaT4EsH+cB64P0Ef1xm7b42sy8DFxN8c/Fp4EqC8YGs2tdmdjtwNsGl5rcB/wn8nm72bRis3yPo8msE3u/uS/q0XoWJiIikSt1cIiKSMoWJiIikTGEiIiIpU5iIiEjKFCYiIpIyhYlIiszsb+HzeDP71zR/9he6W5fIQKOvBoukSfK5KkfwnkTStaG6m7/P3UvTUJ5IRunIRCRFZrYvfHkDcKaZPRPeOyNuZl83s8XhvSI+FC5/tpktMrN7CM7Cxsx+b2ZLw/ttXBW23UBwldtnzOxXyesKz1j+enhvjufM7OKkz34k6V4lvwpPTBPJqMRrLyIivbSApCOTMBTq3f1kMysAnjCz+8NlXwfMcPcN4fQHwjOSi4DFZna3uy8ws6vdfWY363o7MJPgXiSV4XseC+edCEwHXgaeILgG1ePp3liRZDoyEcmcNxFc9+gZgsvSDCe4CRHAU0lBAvBxM3sWeJLgwnuTObwzgNvdvd3dtwGPAicnffZmd+8AngHGp2FbRA5LRyYimWPAx9x94SGNwdjK/i7T5wFz3L3RzB4BClNYb/K1pdrR/3PpBzoyEUmfvUBZ0vRC4CPhpf0xsynhTai6Kgd2h0FyDMFtkzu1dr6/i0XAxeG4TBXBnROfSstWiPSB/mIRSZ/lQHvYXfUzgnujjAeWhYPgdXR/W9j7gA+b2WpgDUFXV6ebgeVmtiy8PH6n3wFzgGcJbnT0OXevDcNIpN/pq8EiIpIydXOJiEjKFCYiIpIyhYmIiKRMYSIiIilTmIiISMoUJiIikjKFiYiIpOz/Ay4SCzmkSGOrAAAAAElFTkSuQmCC\n",
"text/plain": [
"