{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# 变分量子奇异值分解\n",
"\n",
" Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 概览\n",
"\n",
"在本教程中,我们一起学习下经典奇异值分解(singular value decomposition, SVD)的概念以及我们自主研发的量子神经网络版本的量子奇异值分解(variational quantum SVD, VQSVD)[1] 是如何运作的。主体部分包括两个具体案例:\n",
"- 分解随机生成的 $8\\times8$ 复数矩阵;\n",
"- 应用在图像压缩上的效果。\n",
"\n",
"## 背景\n",
"\n",
"奇异值分解有非常多的应用,包括主成分分析(principal component analysis, PCA)、求解线性方程组和推荐系统。其主要任务是给定一个复数矩阵 $M \\in \\mathbb{C}^{m \\times n}$, 找到如下的分解形式:$M = UDV^\\dagger$。其中 $U_{m\\times m}$ 和 $V^\\dagger_{n\\times n}$ 是酉矩阵(Unitary matrix), 满足性质 $UU^\\dagger = VV^\\dagger = I$。 \n",
"\n",
"- 矩阵 $U$ 的列向量 $|u_j\\rangle$ 被称为左奇异向量(left singular vectors), $\\{|u_j\\rangle\\}_{j=1}^{m}$ 组成一组正交向量基。这些列向量本质上是矩阵 $MM^\\dagger$ 的本征向量。\n",
"- 类似的,矩阵 $V$ 的列向量 $\\{|v_j\\rangle\\}_{j=1}^{n}$ 是 $M^\\dagger M$ 的本征向量也组成一组正交向量基。\n",
"- 中间矩阵 $D_{m\\times n}$ 的对角元素上存储着由大到小排列的奇异值 $d_j$。 \n",
"\n",
"我们不妨先来看个简单的例子(为了方便讨论,我们假设以下出现的 $M$ 都是方阵):\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",
"那么该矩阵的奇异值分解可表示为:\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",
"我们通过下面几行代码引入必要的 library和 package。\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.433321Z",
"start_time": "2021-03-09T03:47:28.730090Z"
}
},
"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",
"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": [
"## 经典奇异值分解\n",
"\n",
"那么在了解一些简单的数学背景之后, 我们来学习下如何用 Numpy 完成矩阵的奇异值分解。"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.469286Z",
"start_time": "2021-03-09T03:47:34.440399Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"我们想要分解的矩阵 M 是:\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": [
"# 生成矩阵 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('我们想要分解的矩阵 M 是:')\n",
"print(M_generator())"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.541244Z",
"start_time": "2021-03-09T03:47:34.489833Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"矩阵的奇异值从大到小分别是:\n",
"[11. 7. 5. 1.]\n",
"分解出的酉矩阵 U 是:\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",
"分解出的酉矩阵 V_dagger 是:\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": [
"# 我们只需要以下一行代码就可以完成 SVD \n",
"U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)\n",
"\n",
"# 打印分解结果\n",
"print(\"矩阵的奇异值从大到小分别是:\")\n",
"print(D)\n",
"print(\"分解出的酉矩阵 U 是:\")\n",
"print(U)\n",
"print(\"分解出的酉矩阵 V_dagger 是:\")\n",
"print(V_dagger)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.600873Z",
"start_time": "2021-03-09T03:47:34.565570Z"
}
},
"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": [
"# 再组装回去,能不能复原矩阵?\n",
"M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))\n",
"print(M_reconst)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"那当然是可以复原成原来的矩阵 $M$ 的!读者也可以自行修改矩阵,试试看不是方阵的情况。\n",
"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 量子奇异值分解\n",
"\n",
"接下来我们来看看量子版本的奇异值分解是怎么一回事。简单的说,我们把矩阵分解这一问题巧妙的转换成了优化问题。通过以下四个步骤:\n",
"\n",
"- 准备一组正交向量基 $\\{|\\psi_j\\rangle\\}$, 不妨直接取计算基 $\\{ |000\\rangle, |001\\rangle,\\cdots |111\\rangle\\}$ (这是3量子比特的情形)\n",
"- 准备两个参数化的量子神经网络 $U(\\theta)$ 和 $V(\\phi)$ 分别用来学习左/右奇异向量\n",
"- 利用量子神经网络估算奇异值 $m_j = \\text{Re}\\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle$\n",
"- 设计损失函数并且利用飞桨来优化\n",
"\n",
"$$\n",
"L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle, \\tag{3}\n",
"$$\n",
"\n",
"其中 $q_1>\\cdots>q_T>0$ 是可以调节的权重(超参数), $T$ 表示我们想要学习到的阶数(rank)或者可以解释为总共要学习得到的奇异值个数。\n",
"\n",
"\n",
"\n",
"### 案例1:分解随机生成的 $8\\times8$ 复数矩阵\n",
"\n",
"接着我们来看一个具体的例子,这可以更好的解释整体流程。"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.656799Z",
"start_time": "2021-03-09T03:47:34.620380Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"我们想要分解的矩阵 M 是:\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",
"矩阵的奇异值从大到小分别是:\n",
"[54.83484985 19.18141073 14.98866247 11.61419557 10.15927045 7.60223249\n",
" 5.81040539 3.30116001]\n"
]
}
],
"source": [
"# 先固定随机种子, 为了能够复现结果\n",
"np.random.seed(42)\n",
"\n",
"# 设置量子比特数量,确定希尔伯特空间的维度\n",
"N = 3\n",
"\n",
"# 制作随机矩阵生成器\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",
" M1 = 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",
"print('我们想要分解的矩阵 M 是:')\n",
"print(M)\n",
"\n",
"U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
"print(\"矩阵的奇异值从大到小分别是:\")\n",
"print(D)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.692350Z",
"start_time": "2021-03-09T03:47:34.671114Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"选取的等差权重为:\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": [
"# 超参数设置\n",
"N = 3 # 量子比特数量\n",
"T = 8 # 设置想要学习的阶数\n",
"ITR = 100 # 迭代次数\n",
"LR = 0.02 # 学习速率\n",
"SEED = 14 # 随机数种子\n",
"\n",
"# 设置等差的学习权重\n",
"weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
"print('选取的等差权重为:')\n",
"print(weight)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"我们搭建如下的量子神经网络结构:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:47:34.725007Z",
"start_time": "2021-03-09T03:47:34.702861Z"
}
},
"outputs": [],
"source": [
"# 设置电路参数\n",
"cir_depth = 20 # 电路深度\n",
"block_len = 2 # 每个模组的长度\n",
"theta_size = N * block_len * cir_depth # 网络参数 theta 的大小\n",
"\n",
"\n",
"# 定义量子神经网络\n",
"def U_theta(theta):\n",
"\n",
" # 用 UAnsatz 初始化网络\n",
" cir = UAnsatz(N)\n",
" \n",
" # 搭建层级结构:\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], \n",
" which_qubit)\n",
" \n",
" for which_qubit in range(N):\n",
" cir.rz(theta[(block_len * layer_num + 1) * N \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": [
"接着我们来完成算法的主体部分:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2021-03-09T03:49:24.011232Z",
"start_time": "2021-03-09T03:47:40.712257Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter: 0 loss: 259.1903\n",
"iter: 10 loss: -1672.2268\n",
"iter: 20 loss: -2097.8083\n",
"iter: 30 loss: -2242.9914\n",
"iter: 40 loss: -2310.1948\n",
"iter: 50 loss: -2340.2588\n",
"iter: 60 loss: -2357.5961\n",
"iter: 70 loss: -2369.2625\n",
"iter: 80 loss: -2376.9403\n",
"iter: 90 loss: -2373.7986\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAodklEQVR4nO3de3xcdZ3/8ddnZnK/NG2T3i8ptEVbpBTDXZAfoBQFqy4isApedllc8bK6u4K4v119yO+Bq7teFldFZUVXQRZEWe6giIAWSAFLS62EttD0mrRpmzb3zOf3x/mmnYa0JJlMJpl5Px+Pecyc7zkz53M6MO+c7znne8zdERERSUcs2wWIiMj4pzAREZG0KUxERCRtChMREUmbwkRERNKmMBERkbQpTERkSMxsn5kdle06ZGxRmMi4Y2YbzezcLK37NDP7jZm1mtkeM/tfM1s0ius/sO1m9iEzeyLD6/utmf1Vapu7l7v7+kyuV8YfhYnIIJnZqcBDwK+AGcA84I/AkyP9l7pFMvr/p5klMvn5kl8UJpIzzKzIzL5hZlvC4xtmVhTmVZvZPWa228x2mdnjfT/WZvY5M9sc9jbWmdk5h1nFvwI/dvdvunuru+9y9y8AK4B/CZ+11swuSKkpYWZNZnZCmD7FzH4f6vijmZ2Vsuxvzex6M3sSaAMOG1Bm9kbgu8Cpodtpd8q/wdfM7FUz225m3zWzkjDvLDNrDNu7DfgvM5sY/l2azKwlvJ4Vlr8eOAO4MazjxtDuZjY/vJ5gZj8O73/FzL6Q8u/6ITN7ItTTYmYbzOz8QX+hMq4oTCSXXAecAhwPLAFOAr4Q5n0WaARqgKnA5wE3s2OAq4ET3b0COA/Y2P+DzawUOA34nwHWezvwtvD6VuDSlHnnAc3u/qyZzQTuBb4MTAL+HrjTzGpSlv8gcCVQAbxyuA1197XAVcAfQrdTVZh1A7Aw/BvMB2YC/zflrdPCuueG9cSA/wrTc4B24MawjuuAx4GrwzquHqCU/wAmEAXfW4HLgQ+nzD8ZWAdUE4XxD83MDrddMn4pTCSX/CXwJXff4e5NwBeJfpwBuoHpwFx373b3xz0amK4XKAIWmVmBu29095cH+OxJRP+/bB1g3laiH0uAnwHvCuEDcBlRwAB8ALjP3e9z96S7PwzUA+9I+awfufsad+9x9+6hbHz4kb4S+Luw19QK/D/gkpTFksA/u3unu7e7+053v9Pd28Ly1xOFwmDWFw+ffW3YU9sI/BsH/80BXnH377t7L3AL0XcwdSjbJeODwkRyyQwO/Wv+ldAG8FWgAXjIzNab2TUA7t4AfJqom2qHmd1mZjN4rRaiH+LpA8ybDjSnfN5a4MIQKO8iChiI/vp/X+ji2h26pt7S7zM3DWWD+6kBSoGVKZ//QGjv0+TuHX0TZlZqZt8LXVR7gd8BVSEoXk81UMBr/81npkxv63vh7m3hZfkQtknGCYWJ5JItRD/YfeaENsJfzp9196OIfuA/03dsxN1/5u5vCe914Cv9P9jd9wN/AN43wHovBn6dMt3X1bUceDEEDERB8RN3r0p5lLn7DamrGsL29l+2maibanHK509w9/IjvOezwDHAye5eCZwZ2u0wy/dfXzev/TffPIRtkByhMJHxqsDMilMeCaIf8S+YWY2ZVRMdK/hvADO7wMzmh66gPUTdW0kzO8bMzg4H6juIfoyTh1nnNcAVZvZJM6sIB6+/DJxK1KXW5zbg7cDHOLhXQqjlQjM7z8zioe6z+g54D8N2YJaZFQK4exL4PvB1M5sStnummZ13hM+oINrm3WY2CfjnAdYx4IkAoevqduD68O8xF/hM2E7JMwoTGa/uI/oR7Hv8C9GB7XpgFfAC8GxoA1gAPALsI9rD+E93f5ToeMkNRH9lbwOmANcOtEJ3f4LogPp7iY6TvAIsBd7i7i+lLLc1rOM04Ocp7ZuI9lY+DzQR7an8A8P///A3wBpgm5k1h7bPEXXnrQjdVo8Q7XkczjeAEqLtX0HULZbqm8BF4Wysbw3w/k8A+4H1wBNE4XnzsLZGxjXTzbFERCRd2jMREZG0KUxERCRtChMREUmbwkRERNKWtwO9VVdXe21tbbbLEBEZV1auXNns7jX92/M2TGpra6mvr892GSIi44qZDThmnLq5REQkbQoTERFJm8JERETSpjAREZG0KUxERCRtChMREUmbwkRERNKmMBmiu55r5L9XHPbW3CIieUlhMkT3vbBNYSIi0o/CZIhqKopoau3MdhkiImOKwmSIasqL2NXWRXfv4e7sKiKSfxQmQ1RTUYQ77Nrfle1SRETGDIXJENVUFAGoq0tEJIXCZIgUJiIir6UwGaKacoWJiEh/CpMhOrBnsk9hIiLSR2EyRMUFcSqLE9ozERFJoTAZBl1rIiJyKIXJMChMREQOpTAZhpqKYna0dmS7DBGRMUNhMgw15dozERFJpTAZhpqKIvZ39bK/syfbpYiIjAkKk2HoOz24WacHi4gACpNh0VXwIiKHypkwMbNlZrbOzBrM7JpMrktXwYuIHConwsTM4sC3gfOBRcClZrYoU+ubUqmr4EVEUuVEmAAnAQ3uvt7du4DbgOWZWtnE0kLiMdOeiYhIkCthMhPYlDLdGNoyIh4zJpcVKkxERIJcCZNBMbMrzazezOqbmprS+ixdBS8iclCuhMlmYHbK9KzQdgh3v8nd69y9rqamJq0V1lQU6ZiJiEiQK2HyDLDAzOaZWSFwCXB3JldYU17Ejr0KExERgES2CxgJ7t5jZlcDDwJx4GZ3X5PJddZUFNG8r5Nk0onFLJOrEhEZ83IiTADc/T7gvtFaX01FET1JZ3d7N5PKCkdrtSIiY1KudHONOl0FLyJykMJkmHQVvIjIQQqTYTp4L3jd10RERGEyTFMqiwHtmYiIgMJk2MoK45QUxBUmIiIoTIbNzHQVvIhIoDBJg66CFxGJKEzSoKvgRUQiCpM0TKksYvtenc0lIqIwScPUymL2dvTQ3tWb7VJERLJKYZKGaeH04G3aOxGRPKcwScO0CSFM9ihMRCS/KUzSMDXsmei4iYjkO4VJGg7smShMRCTPKUzSUF6UoLwooT0TEcl7CpM0TdXpwSIiCpN0Ta0s1gF4Ecl7CpM0TassZruugheRPKcwSdPUCcVs39tBMunZLkVEJGsUJmmaVllMT9LZub8r26WIiGSNwiRNutZERERhkjZdBS8iojBJm8bnEhFRmKSturyQmKmbS0Tym8IkTYl4jJqKInVziUheU5iMgGmVxermEpG8pjAZAVMri9XNJSJ5TWEyAqZN0JAqIpLfFCYjQLfvFZF8N+bCxMz+xcw2m9nz4fGOlHnXmlmDma0zs/NS2peFtgYzu2a0a9aFiyKS7xLZLuAwvu7uX0ttMLNFwCXAYmAG8IiZLQyzvw28DWgEnjGzu939xdEqNvVak9rqstFarYjImDFWw2Qgy4Hb3L0T2GBmDcBJYV6Du68HMLPbwrKjFyYTigDtmYhI/hpz3VzB1Wa2ysxuNrOJoW0msCllmcbQdrj21zCzK82s3szqm5qaRqzYvm4uHYQXkXyVlTAxs0fMbPUAj+XAd4CjgeOBrcC/jdR63f0md69z97qampqR+lgqigsoK4zrWhMRyVtZ6eZy93MHs5yZfR+4J0xuBmanzJ4V2jhC+6jpu6+JiEg+GnPdXGY2PWXyPcDq8Ppu4BIzKzKzecAC4GngGWCBmc0zs0Kig/R3j2bNEK6CVzeXiOSpsXgA/l/N7HjAgY3A3wC4+xozu53owHoP8HF37wUws6uBB4E4cLO7rxntoqdVFvPUhl2jvVoRkTFhzIWJu3/wCPOuB64foP0+4L5M1vV6ZlSVsG1vB929SQriY26HT0Qko/SrN0LmTi6lN+k0trRnuxQRkVGnMBkh88LFihub92e5EhGR0acwGSF9V75v3KkwEZH8ozAZIZPLCqkoSmjPRETyksJkhJgZc6tL2bCzLduliIiMOoXJCKqdXKY9ExHJSwqTETSvuozGlja6epLZLkVEZFQpTEZQ7eQykg6NLerqEpH8ojAZQTqjS0TylcJkBNVOLgVgQ7P2TEQkvyhMRtCkskIqinV6sIjkH4XJCDIz5lWXqZtLRPKOwmSE1U5WmIhI/lGYjLDayaVsbmnX6cEiklcUJiOstjo6PfjVXToILyL5Q2Eywmo1erCI5CGFyQibN1nXmohI/lGYjLCJZYVMKClQmIhIXlGYZEDt5FI26sJFEckjCpMMqK0uY4OOmYhIHlGYZMC86jK27Gmnrasn26WIiIwKhUkGHDdrAu6wqnFPtksRERkVCpMMWDp7IgDPvtqS5UpEREaHwiQDJpYVclR1Gc++sjvbpYiIjAqFSYYsnTOR515twd2zXYqISMYpTDJk6Zwqdu7vYtOu9myXIiKScQqTDDlhjo6biEj+UJhkyDHTKigrjCtMRCQvZCVMzOx9ZrbGzJJmVtdv3rVm1mBm68zsvJT2ZaGtwcyuSWmfZ2ZPhfafm1nhaG7L4cRjxpLZVQoTEckL2dozWQ28F/hdaqOZLQIuARYDy4D/NLO4mcWBbwPnA4uAS8OyAF8Bvu7u84EW4KOjswmv74Q5E1m7tZX2rt5slyIiklGDChMz+5SZVVrkh2b2rJm9fbgrdfe17r5ugFnLgdvcvdPdNwANwEnh0eDu6929C7gNWG5mBpwN3BHefwvw7uHWNdKWzqmiN+msatyd7VJERDJqsHsmH3H3vcDbgYnAB4EbMlDPTGBTynRjaDtc+2Rgt7v39GsfkJldaWb1Zlbf1NQ0ooUPZOmBg/C7M74uEZFsGmyYWHh+B/ATd1+T0jbwG8weMbPVAzyWp1NwOtz9Jnevc/e6mpqajK9vUlkh86rLdNxERHJeYpDLrTSzh4B5wLVmVgEc8Sbn7n7uMOrZDMxOmZ4V2jhM+06gyswSYe8kdfkxYemcKn7352bcnahXTkQk9wx2z+SjwDXAie7eBhQAH85APXcDl5hZkZnNAxYATwPPAAvCmVuFRAfp7/bo8vJHgYvC+68AfpWBuoatbu4kmvd18tKOfdkuRUQkYwYbJqcC69x9t5l9APgCMOwhcc3sPWbWGD73XjN7ECB0n90OvAg8AHzc3XvDXsfVwIPAWuD2sCzA54DPmFkD0TGUHw63rkx426KpxAzu+eOWbJciIpIxNpixo8xsFbAEOA74EfAD4GJ3f2tGq8uguro6r6+vH5V1Xfb9FWzb08GvP/tWdXWJyLhmZivdva5/+2D3THpCl9Jy4EZ3/zZQMZIF5rILl8xgffN+1mzZm+1SREQyYrBh0mpm1xKdEnyvmcWIjpvIICxbPI1EzLhn1dZslyIikhGDDZP3A51E15tsIzpr6qsZqyrHTCwr5PT51dyzaouGpBeRnDSoMAkB8lNggpldAHS4+48zWlmOuXDJDBpb2nl+0+5slyIiMuIGO5zKxUSn6L4PuBh4yswuOvK7JNXbF0+lMB5TV5eI5KTBdnNdR3SNyRXufjnRWFn/lLmyck9lcQFnLqzh3lVbSSbV1SUiuWWwYRJz9x0p0zuH8F4JLlwynW17O/jD+p3ZLkVEZEQNNhAeMLMHzexDZvYh4F7gvsyVlZvOWzyNCSUF/OzpV7NdiojIiBrU2Fzu/g9m9hfA6aHpJne/K3Nl5abigjh/ccIsfrJiI837OqkuL8p2SSIiI2LQXVXufqe7fyY8FCTDdNnJs+nude5Y2ZjtUkRERswRw8TMWs1s7wCPVjPT5dzDMH9KBSfNm8StT7+qA/EikjOOGCbuXuHulQM8Kty9crSKzDWXnTSHV3a26UC8iOQMnZGVBcuOncbE0gJ+9pQOxItIblCYZEHfgfgH12yjqbUz2+WIiKRNYZIll548h56k8z8rN73+wiIiY5zCJEuOrinn5HmTuO3pTToQLyLjnsIkiy47eQ6v7mrjyZebs12KiEhaFCZZtOzYaUwqK9SBeBEZ9xQmWVSUiHPRm2fx8Ivb2dHake1yRESGTWGSZZecODs6EF+vK+JFZPxSmGTZUTXlnHrUZG57RlfEi8j4pTAZAy47eQ6bdrXzeIMOxIvI+KQwGQPevngqE0sLuL1e15yIyPikMBkDihJx3r10Jg+v2U7L/q5slyMiMmQKkzHi/SfOpqs3yV3Pbc52KSIiQ6YwGSPeMK2SJbMmcHv9Jtx1IF5ExheFyRhy8Ymz+dO2Vl7YvCfbpYiIDInCZAy5cMkMigti/PwZHYgXkfFFYTKGVBYX8I5jp3P381to7+rNdjkiIoOWlTAxs/eZ2RozS5pZXUp7rZm1m9nz4fHdlHlvNrMXzKzBzL5lZhbaJ5nZw2b2UniemI1tGikXnzib1s4e7l+9NduliIgMWrb2TFYD7wV+N8C8l939+PC4KqX9O8BfAwvCY1lovwb4tbsvAH4dpsetk+dNonZyKbc9ra4uERk/shIm7r7W3dcNdnkzmw5UuvsKj051+jHw7jB7OXBLeH1LSvu4ZGZcctIcnt64i4YdrdkuR0RkUMbiMZN5ZvacmT1mZmeEtplA6kiIjaENYKq79/UJbQOmHu6DzexKM6s3s/qmpqYRL3ykXPTmWRTETXsnIjJuZCxMzOwRM1s9wGP5Ed62FZjj7kuBzwA/M7PKwa4z7LUc9iINd7/J3evcva6mpmbQ2zLaqsuLeNuiqdz5bCOdPToQLyJjXyJTH+zu5w7jPZ1AZ3i90sxeBhYCm4FZKYvOCm0A281surtvDd1hO9KrfGy49KQ53PfCNh5cs513LZmR7XJERI5oTHVzmVmNmcXD66OIDrSvD91Ye83slHAW1+XAr8Lb7gauCK+vSGkf104/uprZk0q4VXdhFJFxIFunBr/HzBqBU4F7zezBMOtMYJWZPQ/cAVzl7rvCvL8FfgA0AC8D94f2G4C3mdlLwLlhetyLxYxLTpzDH9bvZEPz/myXIyJyRJav40DV1dV5fX19tss4oh17Ozj1ht/wkdNrue6di7JdjogIZrbS3ev6t4+pbi451JTKYpYtnsbt9Y26Il5ExjSFyRh3xWm17Gnv5pfPa2h6ERm7FCZj3Im1E3nj9Ep+9ORGDU0vImOWwmSMMzM+dNpc1m1vZcX6Xa//BhGRLFCYjAPLj59JVWkBt/x+Y7ZLEREZkMJkHCguiPP+E2fz0Ivb2Ly7PdvliIi8hsJknPjgKXMBtHciImOSwmScmDWxlAuXzOAnf3iFnfs6s12OiMghFCbjyCfOXkBHTy/ff3xDtksRETmEwmQcmT+lnAuPm8GP/7CRXfu7sl2OiMgBCpNx5pPnzKe9u5fvP74+26WIiBygMBln5k+p4ILjZnDL77V3IiJjh8JkHPrk2dHeyfceeznbpYiIAAqTcWnB1Areu3QWNz+5gfVN+7JdjoiIwmS8+tz5x1CUiPOle17UmF0iknUKk3FqSkUxnz53Ab9d18Sv1+bEnYpFZBxTmIxjV5xWy4Ip5XzxnjV0dOt+JyKSPQqTcawgHuOL71rMpl3tfOe3OhgvItmjMBnnTptfzfLjZ3Djow0892pLtssRkTylMMkBX1p+LNMqi/nkbc/R2tGd7XJEJA8pTHLAhJICvnnJ8Wxuaeeffrk62+WISB5SmOSIutpJfOqchfzy+S384tnGbJcjInlGYZJDrj57PifNm8R1d63mxS17s12OiOQRhUkOiceMGy9byoSSAq78ST0tGrtLREaJwiTHTKko5nsffDM7Wju5+tZn6elNZrskEckDCpMctGR2Fde/+1iebNjJl+9dq+FWRCTjEtkuQDLjfXWzWbu1lZuf3EB1eSFXn70g2yWJSA5TmOSwL7zzjexu6+JrD/2ZsqIEHz59XrZLEpEclZVuLjP7qpn9ycxWmdldZlaVMu9aM2sws3Vmdl5K+7LQ1mBm16S0zzOzp0L7z82scJQ3Z8yKxYx/veg4zls8lS/+74vc/symbJckIjkqW8dMHgaOdffjgD8D1wKY2SLgEmAxsAz4TzOLm1kc+DZwPrAIuDQsC/AV4OvuPh9oAT46qlsyxiXiMb516VLOWFDN536xip899Wq2SxKRHJSVMHH3h9y9J0yuAGaF18uB29y90903AA3ASeHR4O7r3b0LuA1YbmYGnA3cEd5/C/DuUdqMcaMoEef7l9dx1sIaPn/XC/xA948XkRE2Fs7m+ghwf3g9E0jti2kMbYdrnwzsTgmmvvYBmdmVZlZvZvVNTU0jVP74UFwQ53sfrOP8Y6fx5XvX8o1H/qyzvERkxGQsTMzsETNbPcBjecoy1wE9wE8zVUcqd7/J3evcva6mpmY0VjmmFCZi/MelS3nvCTP5xiMv8ZEfPUPzvs5slyUiOSBjZ3O5+7lHmm9mHwIuAM7xg38ibwZmpyw2K7RxmPadQJWZJcLeSeryMoBEPMa/vW8Jx8+u4sv3ruX8bz7Ov1+8hDMW5F+4isjIydbZXMuAfwTe5e5tKbPuBi4xsyIzmwcsAJ4GngEWhDO3CokO0t8dQuhR4KLw/iuAX43WdoxXZsblp9byq4+fzoSSAj74w6f5xK3PsXl3e7ZLE5FxyrLRb25mDUAR0Z4FwAp3vyrMu47oOEoP8Gl3vz+0vwP4BhAHbnb360P7UUQH5CcBzwEfcPfX7bupq6vz+vr6kdyscam9q5fvPPYy33ssulPjX59xFJefOpcplcVZrkxExiIzW+nuda9pz9eDsAqTQzW2tPGVB9bxv3/cQjxmnLWwhotPnM1Zx9RQlIhnuzwRGSMUJv0oTAa2oXk/t9dv4o6VjTS1dlJRlODcRVN5x5umc8aCaooLFCwi+Uxh0o/C5Mh6epM83tDMfau28tCL29nT3k1ZYZz/84YpnH/sdM554xQFi0geUpj0ozAZvK6eJL9/uZkH12zjoTXb2bm/i8riBO9ZOpP3nziHRTMqs12iiIwShUk/CpPh6U06K9bv5Pb6Tdy/ehtdPUlOmjeJj731aM46poZoUAIRyVUKk34UJunb3dbFHSsbufmJDWzZ08EbplXwsbOO5p1vmk4iPhYGVxCRkaYw6UdhMnK6e5Pc/fwWvvvYy7y0Yx9zJ5dy1VuP5r0nzNSZYCI5RmHSj8Jk5CWTzkMvbufbjzbwwuY91FQUccWpc/nLk+cysUx3BhDJBQqTfhQmmePuPNHQzE2/W8/jLzVTUhDnguOmc8GSGZx29GQK1AUmMm4dLkx0p0UZcWbGGQtqOGNBDeu2tXLzExu474Wt/M/KRqpKCzj7DVM4/ehqTp9fzbQJutJeJBdoz0RGRUd3L4+/1My9q7bw2J+baGnrBuCo6jJOrJ1EXe1ETp43mTmTS7NcqYgcibq5+lGYZE8y6azdtpcnG5p5esMuntnYwp72KFxmTyrhLfOrOXNBDWcurKGsSDvPImOJwqQfhcnYkUw6DU37WLF+J0+81MwfXt5Ja2cPRYkYb11Yw/lvmsbZb5jKhJKCbJcqkvcUJv0oTMaunt4k9a+0cP8LW3lgzTa27+2kIG6cenQ15y2eypkLapg9Sd1hItmgMOlHYTI+JJPOc5t289Cabdy/ehuv7opufzOvuozT50/muJlVLJpRycKpFRQmdJaYSKYpTPpRmIw/7s5LO/bx+EvNPPFSE09t2EVbVy8AiZgxbUIxsyaWMLOqlBlVxUyfUML0vraJJZQW6viLSLp0arCMe2bGwqkVLJxawUffMo/epPPKzv2s2bKXP23bS2NLO5tb2vn9y81s39tBst/fSZPKCpk9qZTayaXMnVTK3Mll1FZHz5PLCjWumEgatGciOamnN8mO1k627mmnsaXv0caru9rY2NzG1j3th4RNSUGcGVXFzKgqYVplMVMqi6gpL2JKZTHV5UXUVBRRXV5IeVFCoSN5TXsmklcS8RgzqkqYUVXCm+e+dn5nTy+NLe28urONjTv3s7mlnc2729myu50/b2+leV8Xvf13bYCiRIzq8iImlxcyuayQSWXR60llhUwqLaSqtICq8DyhpICK4gQlBXEFkOQ8hYnkpaJEnKNryjm6pnzA+b1JZ9f+LppaO2ned/Cxc18XTfs6aQ7P67a10ry/i66e5GHXFY8ZFcUJKooTVBYXHHieUBI9qkoLmFBaGL1ObSspoLwooRGYZVxQmIgMIB4zaiqi7q3X4+60d/eya38XLfu72dPeze72Lna3ddPa0UNrx8HnvR097G3vZuPO/ext72FPezft3b1H/PyywjjlxQnKihKUFyUoLYxTWhjt8RQVxCguiFOciFMcXhclYpQUhrbCOCUFcUoL4xSH59LCOCWhvaQgrrCSEaEwEUmTmVFamKC0MMGsiUN/f2dPL3vau9nTFgXRnvZudrd1s7ej+0Dg7O/sYV94tHf1sqO1g7auXjq7k3T29NLRnaSju5eeAbrmXk9hPHYgiIpDwBQXxCgK08WJ8DoRo6ggRlEiCqyiRDxMxyhMxCiMR88Hp+MUJmIUxO3A/IKwTEE8TCeMgniMRMzGXVfgnrZu1m1vZd32VrbsbqeiOMHE0kKmVBRxylGT8270hvzaWpExqCgRZ0pFnCkV6Q962dObpKMnCpa+R1tXL+1dvbR1R88HX/fQ0Z2kPbR3dPfSfuB90WfsaetiR0+SzvCZnT1JOrt76ehJDnhMabjMOBgw8Shg+oInEeubNhLxg9OJuJGIHdqeiEWvC8K8aJmU+XGjIBYjHl7HzIjHjLgZZtEfBn2R1pt0et3p6O49sGfZvK+LV3bu55Wdbezc33Wg/phxyAkdRYkYZy6s4W2LpjJ9QjGlhQkK4saftrXyx027eXHrXhIxC12dhRxVU8aSWVW8adYEEjFjR2snTa2d7Nrfxe62Lna3d1NWGOfoKeUsmFJBdfnBsw87e3ppau1kR2snnd3JA3vUlcWvPVmkuzdJS1sXk0oLR3yPVGEikkMS8Rjl8Rjlo/BXcU9vkq7eJJ3d0XNXCJ2unoPT3b2vne4+8NoPtEef5dHr5MH39PT6gfd09zo9yei5rauHnqRHbb3J8DoKuL7lenud7mT0GcPZY+uvuCDGpNJC5k4u4+2Lp1I7uYyF0yo4ZmoF0ycU09Ed/VC/srONB9ds44HV23j4xe2v+ZyK4gSLZ1QSM2PL7g7WbNnLnc82DrkeMzB4zSnwfRIxO9A1mogbLfu72NvRA8Cjf38W86rLhrzOI1GYiMiwJOIxEvEYpePgvmfuTm8yCpWepB8ImmTY++hNOu5ED6Jf53gs2rspTMSoKE687n14SgrjlBRGZxCeevRk/u8Fi2ho2negm7KzJ8nRNeUcVV1GLHboHsPuti5e2LyH1Zv3YgY14XT0SWWFTCwrZGJpAXvbe2jYsY+XdrTSsr8LJ6q3KBGLTmWvKKIoEad5X7RXs3N/F22dPezr7KW7N8nE0gImlkVnIVZlYJw7XWciIiKDdrjrTHQah4iIpE1hIiIiaVOYiIhI2rISJmb2VTP7k5mtMrO7zKwqtNeaWbuZPR8e3015z5vN7AUzazCzb1k4583MJpnZw2b2Ungexpn+IiKSjmztmTwMHOvuxwF/Bq5Nmfeyux8fHleltH8H+GtgQXgsC+3XAL929wXAr8O0iIiMoqyEibs/5O49YXIFMOtIy5vZdKDS3Vd4dPrZj4F3h9nLgVvC61tS2kVEZJSMhWMmHwHuT5meZ2bPmdljZnZGaJsJpF7V0xjaAKa6+9bwehswNaPViojIa2TsokUzewSYNsCs69z9V2GZ64Ae4Kdh3lZgjrvvNLM3A780s8WDXae7u5kd9sIZM7sSuBJgzpw5g/1YERF5HRkLE3c/90jzzexDwAXAOaHrCnfvBDrD65Vm9jKwENjMoV1hs0IbwHYzm+7uW0N32I4j1HQTcFNYf5OZvTKETaoGmoewfC7Ix22G/NzufNxmyM/tTnebB7hDUJaGUzGzZcA/Am9197aU9hpgl7v3mtlRRAfa17v7LjPba2anAE8BlwP/Ed52N3AFcEN4/tVganD3miHWXD/QVZ+5LB+3GfJzu/NxmyE/tztT25ytsbluBIqAh8MZvivCmVtnAl8ys24gCVzl7rvCe/4W+BFQQnSMpe84yw3A7Wb2UeAV4OLR2ggREYlkJUzcff5h2u8E7jzMvHrg2AHadwLnjGiBIiIyJGPhbK7x4qZsF5AF+bjNkJ/bnY/bDPm53RnZ5rwdNVhEREaO9kxERCRtChMREUmbwuR1mNkyM1sXBpjM2XG/zGy2mT1qZi+a2Roz+1Roz/mBNM0sHkZduCdMzzOzp8J3/nMzGwf3EhwaM6syszvCgKtrzezUXP+uzezvwn/bq83sVjMrzsXv2sxuNrMdZrY6pW3A79Yi3wrbv8rMThjuehUmR2BmceDbwPnAIuBSM1uU3aoypgf4rLsvAk4BPh62NR8G0vwUsDZl+ivA18NZhy3AR7NSVWZ9E3jA3d8ALCHa/pz9rs1sJvBJoM7djwXiwCXk5nf9Iw4OhNvncN/t+RwcPPdKogF1h0VhcmQnAQ3uvt7du4DbiAaWzDnuvtXdnw2vW4l+XGaS4wNpmtks4J3AD8K0AWcDd4RFcnGbJxBd0/VDAHfvcvfd5Ph3TXQpRImZJYBSouGbcu67dvffAbv6NR/uu10O/NgjK4CqMJLIkClMjmwmsCllOnWAyZxlZrXAUqLRBnJ9IM1vEI3GkAzTk4HdKaNa5+J3Pg9oAv4rdO/9wMzKyOHv2t03A18DXiUKkT3ASnL/u+5zuO92xH7jFCZyCDMrJ7pw9NPuvjd1XhhDLWfOJTezC4Ad7r4y27WMsgRwAvAdd18K7Kdfl1YOftcTif4KnwfMAMp4bVdQXsjUd6swObLNwOyU6dQBJnOOmRUQBclP3f0XoXl7327v6w2kOQ6dDrzLzDYSdWGeTXQsoSp0hUBufueNQKO7PxWm7yAKl1z+rs8FNrh7k7t3A78g+v5z/bvuc7jvdsR+4xQmR/YMsCCc8VFIdMDu7izXlBHhWMEPgbXu/u8ps/oG0oQhDKQ5Hrj7te4+y91rib7b37j7XwKPAheFxXJqmwHcfRuwycyOCU3nAC+Sw981UffWKWZWGv5b79vmnP6uUxzuu70buDyc1XUKsCelO2xIdAX86zCzdxD1q8eBm939+uxWlBlm9hbgceAFDh4/+DzRcZPbgTmEgTRTBt/MGWZ2FvD37n5BGLH6NmAS8BzwgXB7hJxhZscTnXRQCKwHPkz0x2XOftdm9kXg/URnLj4H/BXR8YGc+q7N7FbgLKKh5rcD/wz8kgG+2xCsNxJ1+bUBHw7jIA59vQoTERFJl7q5REQkbQoTERFJm8JERETSpjAREZG0KUxERCRtChORNJnZ78NzrZldNsKf/fmB1iUy1ujUYJERknqtyhDek0gZG2qg+fvcvXwEyhPJKO2ZiKTJzPaFlzcAZ5jZ8+HeGXEz+6qZPRPuFfE3YfmzzOxxM7ub6CpszOyXZrYy3G/jytB2A9Eot8+b2U9T1xWuWP5quDfHC2b2/pTP/m3KvUp+Gi5ME8moxOsvIiKDdA0peyYhFPa4+4lmVgQ8aWYPhWVPAI519w1h+iPhiuQS4Bkzu9PdrzGzq939+AHW9V7geKJ7kVSH9/wuzFsKLAa2AE8SjUH1xEhvrEgq7ZmIZM7bicY9ep5oWJrJRDchAng6JUgAPmlmfwRWEA28t4Ajewtwq7v3uvt24DHgxJTPbnT3JPA8UDsC2yJyRNozEckcAz7h7g8e0hgdW9nfb/pc4FR3bzOz3wLFaaw3dWypXvT/uYwC7ZmIjJxWoCJl+kHgY2Fof8xsYbgJVX8TgJYQJG8gum1yn+6+9/fzOPD+cFymhujOiU+PyFaIDIP+YhEZOauA3tBd9SOie6PUAs+Gg+BNDHxb2AeAq8xsLbCOqKurz03AKjN7NgyP3+cu4FTgj0Q3OvpHd98Wwkhk1OnUYBERSZu6uUREJG0KExERSZvCRERE0qYwERGRtClMREQkbQoTERFJm8JERETS9v8Bb+rX1oZ3hjIAAAAASUVORK5CYII=\n",
"text/plain": [
"