{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 变分量子奇异值分解 (VQSVD)\n",
    "\n",
    "<em> Copyright (c) 2020 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>\n",
    "\n",
    "\n",
    "## 概览\n",
    "\n",
    "在本教程中，我们一起学习下经典奇异值分解（SVD）的概念以及我们自主研发的量子神经网络版本的量子奇异值分解 （VQSVD，Variational Quantum Singular Value Decomposition）[1] 是如何运作的。主体部分包括两个具体案例：\n",
    "\n",
    "- 分解随机生成的 8x8 复数矩阵\n",
    "- 应用在图像压缩上的效果\n",
    "\n",
    "## 背景\n",
    "\n",
    "奇异值分解（SVD）有非常多的应用包括 -- 主成分分析（PCA）、求解线性方程组和推荐系统。其主要任务是给定一个复数矩阵 $M \\in \\mathbb{C}^{m \\times n}$, 找到分解形式：\n",
    "\n",
    "$$\n",
    "M = UDV^\\dagger\n",
    "$$\n",
    "\n",
    "其中 $U_{m\\times m}$ 和 $V^\\dagger_{n\\times n}$ 是酉矩阵（Unitary matrix）, 满足性质 $UU^\\dagger = VV^\\dagger = I$。 \n",
    "- 矩阵 $U$ 的列向量 $\\lvert {u_j}\\rangle$ 被称为左奇异向量（left singular vectors）, $\\{\\lvert {u_j}\\rangle\\}_{j=1}^{m}$ 组成一组正交向量基。这些列向量本质上是矩阵 $MM^\\dagger$ 的特征向量。\n",
    "- 类似的，矩阵 $V$ 的列向量 $\\{\\lvert {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}\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}\n",
    "$$\n",
    "\n",
    "首先，让我们通过下面几行代码引入必要的 library和 package。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "from progressbar import ProgressBar\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import unitary_group\n",
    "from scipy.linalg import norm\n",
    "\n",
    "import paddle.fluid as fluid\n",
    "from paddle.complex 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()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 经典奇异值分解\n",
    "\n",
    "那么在了解一些简单的数学背景之后， 我们来学习下如何用 Numpy 完成矩阵的奇异值分解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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": {},
   "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": {},
   "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",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子奇异值分解\n",
    "\n",
    "接下来我们来看看量子版本的奇异值分解是怎么一回事。简单的说，我们把矩阵分解这一问题巧妙的转换成了优化问题。通过以下四个步骤：\n",
    "\n",
    "- 准备一组正交向量基 $\\{\\lvert {\\psi_j}\\rangle\\}$, 不妨直接取计算基 $\\{ \\lvert {000}\\rangle, \\lvert {001}\\rangle,\\cdots \\lvert {111}\\rangle\\}$ （这是3量子比特的情形）\n",
    "- 准备两个参数化的量子神经网络 $U(\\theta)$ 和 $V(\\phi)$ 分别用来学习左/右奇异向量\n",
    "- 利用量子神经网络估算奇异值 $m_j = \\text{Re}\\langle{\\psi_j} \\lvert U(\\theta)^{\\dagger} M V(\\phi)\\lvert {\\psi_j}\\rangle$\n",
    "- 设计损失函数并且利用飞桨来优化\n",
    "\n",
    "$$\n",
    "L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle{\\psi_j} \\lvert U(\\theta)^{\\dagger} M V(\\phi)\\lvert {\\psi_j}\\rangle\n",
    "$$\n",
    "\n",
    "其中 $q_1>\\cdots>q_T>0$ 是可以调节的权重（超参数）, $T$ 表示我们想要学习到的阶数（rank），或者可以解释为总共要学习得到的奇异值个数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 案例1：分解随机生成的 8x8 复数矩阵\n",
    "\n",
    "接着我们来看一个具体的例子，这可以更好的解释整体流程。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "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))\\\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",
    "# 打印结果\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": {},
   "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 = 1        # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
    "print('选取的等差权重为：')\n",
    "print(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 量子神经网络的构造\n",
    "\n",
    "我们搭建如下的结构：\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置电路参数\n",
    "cir_depth = 40                           # 电路深度\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 + which_qubit], \n",
    "                                                    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": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% |########################################################################|\r"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "主程序段总共运行了 216.55158972740173 秒\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "class NET(fluid.dygraph.Layer):\n",
    "    \n",
    "    # 初始化长度为 theta_size 的可学习参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "    def __init__(self, shape, param_attr=fluid.initializer.Uniform(low=0.0, high=2 * np.pi), dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # 创建用来学习 U 的参数 theta\n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                     attr=param_attr, dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的参数 phi\n",
    "        self.phi = self.create_parameter(shape=shape, \n",
    "                    attr=param_attr, dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 我们需要将 Numpy array 转换成 Paddle 动态图模式中支持的 variable\n",
    "        self.M = fluid.dygraph.to_variable(M)\n",
    "        self.weight = fluid.dygraph.to_variable(weight)\n",
    "\n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self):\n",
    "        \n",
    "        # 获取量子神经网络的酉矩阵表示\n",
    "        U = U_theta(self.theta)\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        \n",
    "        V = U_theta(self.phi)\n",
    "        V_dagger = dagger(V)\n",
    "        \n",
    "        # 初始化损失函数和奇异值存储器\n",
    "        loss = 0 \n",
    "        singular_values = np.zeros(T)\n",
    "        \n",
    "        # 定义损失函数\n",
    "        for i in range(T):\n",
    "            loss -= self.weight.real[i] * matmul(U_dagger, matmul(self.M, V)).real[i][i]\n",
    "            singular_values[i] = (matmul(U_dagger, matmul(self.M, V)).real[i][i]).numpy()\n",
    "        \n",
    "        # 函数返回两个矩阵 U 和 V_dagger、 学习的奇异值以及损失函数    \n",
    "        return U, V_dagger, loss, singular_values\n",
    "    \n",
    "# 记录优化中间过程\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "\n",
    "time_start = time.time()\n",
    "# 启动 Paddle 动态图模式\n",
    "with fluid.dygraph.guard():\n",
    "    \n",
    "    # 确定网络的参数维度\n",
    "    net = NET([theta_size])\n",
    "    \n",
    "    # 一般来说，我们利用Adam优化器来获得相对好的收敛，当然你可以改成SGD或者是RMS prop.\n",
    "    opt = fluid.optimizer.AdagradOptimizer(learning_rate=LR, parameter_list=net.parameters())\n",
    "    \n",
    "    # 优化循环\n",
    "    pbar = ProgressBar()\n",
    "    for itr in pbar(range(ITR)):\n",
    "        \n",
    "        #  前向传播计算损失函数\n",
    "        U, V_dagger, loss, singular_values = net()\n",
    "        \n",
    "        # 在动态图机制下，反向传播极小化损失函数\n",
    "        loss.backward()\n",
    "        opt.minimize(loss)\n",
    "        net.clear_gradients()\n",
    "\n",
    "        # 记录优化中间结果\n",
    "        loss_list.append(loss[0][0].numpy())\n",
    "        singular_value_list.append(singular_values)\n",
    "\n",
    "    # 记录最后学出的两个酉矩阵    \n",
    "    U_learned = U.real.numpy() + 1j * U.imag.numpy()\n",
    "    V_dagger_learned = V_dagger.real.numpy() + 1j * V_dagger.imag.numpy()\n",
    "\n",
    "time_span = time.time() - time_start        \n",
    "print('主程序段总共运行了', time_span, '秒')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAojUlEQVR4nO3deXxU9b3/8dcnk2SyExLCFgIIAipUQIJ1q3rVW7Ha2vZqpb0u3a61t6u1i9b+utzf9Xfr1V5b21ut1bq0WrVal7au1VKVugUXwAUFZAlr2ENClkk+vz/OCQwQIMBMTjLzfj4e88g53zNnzufLMu98z2rujoiIyMHKiboAERHJDAoUERFJCQWKiIikhAJFRERSQoEiIiIpoUAREZGUUKCIyH4xs61mNibqOqTvUaBIv2NmS8zstIi2fZyZPW1mjWa22cz+ZGZH9OL2t/fdzD5tZs+leXuzzOzzyW3uXuLui9O5XemfFCgiPWRmxwJPAA8Bw4FDgNeB2an+jd0Caf3/aWa56fx8yT4KFMkYZhY3s5+a2crw9VMzi4fLBpnZn81sk5ltMLNnu76wzew7ZrYiHHUsMLNT97CJ/wbucPefuXuju29w9+8BLwA/DD/rLTM7K6mmXDNbZ2ZHhfPHmNk/wjpeN7OTk947y8yuMrPZQDOwx5Ays8OBG4Fjw11Qm5L+DK41s2VmtsbMbjSzwnDZyWZWH/Z3NXCrmQ0M/1wazGxjOD0ifP9VwAeAX4Tb+EXY7mZ2aDg9wMzuCNdfambfS/pz/bSZPRfWs9HM3jOzM3r8Fyr9jgJFMsmVwDHAFGAycDTwvXDZZUA9UAUMAb4LuJlNAL4MTHf3UuB0YMmuH2xmRcBxwB+62e69wD+H078HPpm07HRgnbu/YmbVwF+A/wQqgG8C95tZVdL7LwAuBkqBpXvqqLu/BVwCPB/ugioPF10NjA//DA4FqoHvJ606NNz2qHA7OcCt4fxIYBvwi3AbVwLPAl8Ot/Hlbkr5OTCAIPxOAi4EPpO0/P3AAmAQQSDfYma2p35J/6ZAkUzyr8B/uPtad28AfkTwBQ3QDgwDRrl7u7s/68GN7DqAOHCEmeW5+xJ3X9TNZ1cQ/H9Z1c2yVQRfmAB3AR8JAwjgU2EbwPnAI+7+iLt3uvuTQB3woaTPus3d33D3hLu370/nwy/qfwMuDUdPjcD/A2Ymva0T+IG7t7r7Nndf7+73u3tz+P6rCIKhJ9uLAecBV4QjtiXAT9jxZw6w1N1/7e4dwO0EfwdD9qdf0n8oUCSTDGfn3+qXhm0A1wALgSfMbLGZXQ7g7guBrxPsslprZneb2XB2t5Hgy3hYN8uGAeuSPu8t4MNhqHyEHYEyCjg33N21KdxNdcIun7l8fzq8iyqgCJiT9PmPhe1dGty9pWvGzIrM7Ffh7qotwDNAeRgW+zIIyGf3P/PqpPnVXRPu3hxOluxHn6QfUaBIJllJ8KXdZWTYRvgb9GXuPgb4MPCNrmMl7n6Xu58QrusEu4124u5NwPPAud1s9xPAU0nzXbu9zgbeDEMGgrD4rbuXJ72K3f3HyZvaj/7u+t51BLusJiZ9/gB3L9nLOpcBE4D3u3sZcGLYbnt4/67ba2f3P/MV+9EHySAKFOmv8sysIOmVS/BF/j0zqzKzQQTHDn4HYGZnmdmh4W6hLQS7ujrMbIKZnRIevG8h+ELu2MM2LwcuMrOvmllpeED7P4FjCXavdbkb+CDwRXaMTghr+bCZnW5msbDuk7sOgh+ANcAIM8sHcPdO4NfAdWY2OOx3tZmdvpfPKCXo8yYzqwB+0M02uj05INyNdS9wVfjnMQr4RthPyUIKFOmvHiH4Iux6/ZDgYHcdMBeYB7wStgGMA/4KbCUYafzS3WcRHD/5McFv26uBwQQH7Hfj7s8RHGT/OMFxk6XAVOAEd3836X2rwm0cB9yT1L6cYNTyXaCBYMTyLQ78/+HTwBvAajNbF7Z9h2DX3gvhLqy/EoxA9uSnQCFB/18g2EWW7GfAOeFZWtd3s/5XgCZgMfAcQYD+5oB6I/2e6QFbIiKSChqhiIhISihQREQkJRQoIiKSEgoUERFJiay9OdygQYN89OjRUZchItKvzJkzZ527V3W3LGsDZfTo0dTV1UVdhohIv2Jme7zHnHZ5iYhISihQREQkJRQoIiKSEgoUERFJCQWKiIikhAJFRERSQoEiIiIpoUDZTy8v2cDVj72N7tIsIrIzBcp+mlu/mRtmLWLLtkTUpYiI9CkKlP1UVRoHoGFryz7eKSKSXRQo+6mqJAiUtY2tEVciItK3KFD2U1VpPgANChQRkZ0oUPZTVUkBoEAREdmVAmU/lRXmkh/LoWGrAkVEJJkCZT+ZGVWlcY1QRER2oUA5AINK46zb2hZ1GSIifYoC5QBUlWiEIiKyKwXKAdAuLxGR3SlQDkBVaZwNTa10dOr2KyIiXRQoB6CqNE6nw/omjVJERLpkTKCY2QwzW2BmC83s8nRuq+tqee32EhHZISMCxcxiwP8CZwBHAJ80syPStT1dLS8isruMCBTgaGChuy929zbgbuDsdG1MV8uLiOwuUwKlGlieNF8ftu3EzC42szozq2toaDjgjQ3qGqHoankRke0yJVCsm7bdTsFy95vcvdbda6uqqg54Y0X5uZTEc3cboeisLxHJZpkSKPVATdL8CGBlOjdYtcvV8nOWbuTw7z/Gik3b0rlZEZE+K1MC5WVgnJkdYmb5wEzg4XRuMLhafsdDtp57dx1tiU7ea2hK52ZFRPqsjAgUd08AXwYeB94C7nX3N9K5zV2vlp+3YjOga1NEJHvlRl1Aqrj7I8AjvbW9qtI4z767Izzmh4GysUk3jRSR7JQxgdLbqkrjbGlJ0NLeQWNLgtVbgt1fGxQoIpKlFCgHqOtq+XVbW3l37dbt7esVKCKSpTLiGEoUqkp33H5lfn2wu2toWYFGKCKStTRCOUCDku7nNW/FZsYMKmZQSVyBIiJZSyOUA7R9hLK1lfkrNjOpegADi/MUKCKStRQoB6iyJLj9yjurG1m5uYVJ1WVUFGuEIiLZS7u8DlBeLIeK4nz+tiC4J9ik6gFs2ZZgY3MbnZ1OTk53d4MREclcGqEchKqSOMs2NAMwcfgAKorz6XTYtK094spERHqfAuUgdB1HGVVZxIDCPCqKg91g2u0lItlIgXIQugJlUvUAAAWKiGQ1BcpB6AqU9+0WKLqfl4hkHwXKQei6Wn7S8CBQus780tXyIpKNFCgHYfohFUypKWfKyHIABhYFgaIbRIpINtJpwwdhSk05D37p+O3zBXkxivNjGqGISFbSCCXFKkrydVBeRLKSAiXFdLW8iGQrBUqKVRbns36rAkVEso8CJcUqivPZ2KxAEZHso0BJsYrifNY3teHuUZciItKrFCgpVlGcT1uik6a2jqhLERHpVQqUFNt+tbyOo4hIllGgpFhlV6DoOIqIZBkFSorpfl4ikq0UKCnWFSg6dVhEso0CJcV0C3sRyVYKlBQrieeSH8tRoIhI1lGgpJiZUVGs+3mJSPZRoKSBAkVEspECJQ26rpYXEckmCpQ00AhFRLKRAiUNKorz9dRGEck6fS5QzOyHZrbCzF4LXx9KWnaFmS00swVmdnpS+zQzmxcuu97MLJrqA5XF+TS2JmhN6H5eIpI9+lyghK5z9ynh6xEAMzsCmAlMBGYAvzSzWPj+G4CLgXHha0YENW9XUdL1bPn2KMsQEelVfTVQunM2cLe7t7r7e8BC4GgzGwaUufvzHtwz/g7goxHWSUVReLW8br8iIlmkrwbKl81srpn9xswGhm3VwPKk99SHbdXh9K7tuzGzi82szszqGhoa0lE3AFWlcQDWbGlJ2zZERPqaSALFzP5qZvO7eZ1NsPtqLDAFWAX8pGu1bj7K99K+e6P7Te5e6+61VVVVB9+RPZgwtBQzmFu/OW3bEBHpa3Kj2Ki7n9aT95nZr4E/h7P1QE3S4hHAyrB9RDftkSktyGP84FJeW74pyjJERHpVn9vlFR4T6fIxYH44/TAw08ziZnYIwcH3l9x9FdBoZseEZ3ddCDzUq0V3Y+rIcl5dtkmPAhaRrNHnAgX47/AU4LnAPwGXArj7G8C9wJvAY8CX3L3rvNwvAjcTHKhfBDza61XvYkpNOZu3tfPeuqaoSxER6RWR7PLaG3e/YC/LrgKu6qa9DpiUzrr219SRwbkEry3fxJiqkoirERFJv744QskIhw4uoTg/xqvLNkVdiohIr1CgpEksx5hcU64D8yKSNRQoaTSlppy3Vm1hW5tuwSIimU+BkkZTRw4k0enMX6nrUUQk8ylQ0mhKTTkAr+k4iohkAQVKGlWVxhkxsJBXl2+MuhQRkbRToKTZ1JEDNUIRkaygQEmzKTXlrNzcohtFikjGU6Ck2dSR5QDMWardXiKS2RQoafa+6gGUxHN59t11UZciIpJWCpQ0y4vlcNzYSp55p0E3ihSRjKZA6QUnTahixaZtLGrYGnUpIiJpo0DpBSeOCx7m9fd3tNtLRDKXAqUX1FQUMbaqmL+/k77HDouIRE2B0ktOHF/Fi4vX09Ku+3qJSGZSoPSSk8ZX0Zro5IXF66MuRUQkLRQoveSYMZXEc3N4RsdRRCRDKVB6SUFejKMPqeDv76yNuhQRkbRQoPSik8ZXsaihieUbmqMuRUQk5RQovejkCcHpw0+/rVGKiGQeBUovGltVwmFDS/njK/VRlyIiknIKlF5kZpwzbQSv12/mnTWNUZcjIpJSCpRe9rGp1eTmGH+oWx51KSIiKaVA6WWVJXFOOWwwD7y6gvaOzqjLERFJGQVKBM6trWHd1jZmLdCtWEQkcyhQInDyhCoGleRrt5eIZBQFSgTyYjl8bGo1T7+9lnVbW6MuR0QkJRQoETm3toZEp/PgqyuiLkVEJCUUKBEZP6SUyTXl3Fu3XE9yFJGMoECJ0MzpNbyzZiuvLd8UdSkiIgdNgRKhs44cRmFejHte1sF5Een/IgkUMzvXzN4ws04zq91l2RVmttDMFpjZ6Unt08xsXrjsejOzsD1uZveE7S+a2ehe7s4BKy3I48wjh/Gn11fS1JqIuhwRkYPSo0Axs6+ZWZkFbjGzV8zsgwex3fnAx4FndtnOEcBMYCIwA/ilmcXCxTcAFwPjwteMsP1zwEZ3PxS4Drj6IOrqdedNr6GprYO/zFsVdSkiIgelpyOUz7r7FuCDQBXwGeDHB7pRd3/L3Rd0s+hs4G53b3X394CFwNFmNgwoc/fnPTiCfQfw0aR1bg+n7wNO7Rq99Ae1owYypqqYe7XbS0T6uZ4GStcX9IeAW9399aS2VKoGkr9Z68O26nB61/ad1nH3BLAZqOzuw83sYjOrM7O6hoa+cZW6mXFebQ11SzeycK1uGCki/VdPA2WOmT1BECiPm1kpsNcbUZnZX81sfjevs/e2Wjdtvpf2va2ze6P7Te5e6+61VVVVeyu/V338qBHk5hi/f0mjFBHpv3J7+L7PAVOAxe7ebGYVBLu99sjdTzuAeuqBmqT5EcDKsH1EN+3J69SbWS4wANhwANuOTFVpnBmThnLvy8u59J/HUxLv6V+LiEjf0dMRyrHAAnffZGbnA98j2LWUag8DM8Mztw4hOPj+kruvAhrN7Jjw+MiFwENJ61wUTp8DPO398ErBz51wCI2tCd3fS0T6rZ4Gyg1As5lNBr4NLCU4MH5AzOxjZlZPEFR/MbPHAdz9DeBe4E3gMeBL7t4RrvZF4GaCA/WLgEfD9luASjNbCHwDuPxA64rS1JEDOWpkObfOXkJHZ7/LQxERrCe/zJvZK+5+lJl9H1jh7rd0taW/xPSora31urq6qMvYyV/mruJLd73Cry6YxukTh0ZdjojIbsxsjrvXdrespyOURjO7AriAYEQRA/JSVaAETp84hOryQm557r2oSxER2W89DZTzgFaC61FWE5yqe03aqspSubEcPn3caF56bwPzV6TjEJWISPr0KFDCELkTGGBmZwEt7n7Ax1Bkz847uobi/Bi/0ShFRPqZnt565RPAS8C5wCeAF83snHQWlq3KCvI4Z9oI/jx3lR6+JSL9Sk93eV0JTHf3i9z9QuBo4P+kr6zsdsGxo2nr6OTul5ZFXYqISI/1NFBy3H1t0vz6/VhX9tOhg0v4wLhB/O6FZbR37PWGBCIifUZPQ+ExM3vczD5tZp8G/gI8kr6y5MJjR7N6SwtPvLEm6lJERHqkpwflvwXcBBwJTAZucvfvpLOwbHfKYYMZMbCQ259fEnUpIiI90uObRrn7/cD9aaxFksRyjAuOGcV/Pfo2b63awuHDyqIuSURkr/Y6QjGzRjPb0s2r0cy29FaR2eq86TXEc3O4bfaSqEsREdmnvQaKu5e6e1k3r1J316/MaVZelM8500bwwKsrWL25JepyRET2Smdq9XFfOHEsHe7c/OziqEsREdkrBUofN7KyiI9MHs6dLy5jQ1Nb1OWIiOyRAqUf+PeTx7KtvYPbZut2LCLSdylQ+oFxQ0o5feIQbvvHEhpb2qMuR0SkWwqUfuLL/zSOLS0JfveCbsciIn2TAqWfeN+IAZw4vopbnlvMtraOfa8gItLLFCj9yFdOOZR1W9v4vW4aKSJ9kAKlH5k+uoL3H1LBr55ZRGtCoxQR6VsUKP3MV08dx5otrfyhrj7qUkREdqJA6WeOG1vJ1JHl3DBrkW5tLyJ9igKlnzEzvnrKOFZs2sYDr66IuhwRke0UKP3QyROqmFRdxi+eXkhbQqMUEekbFCj9kJlx2QcnsGxDM799YWnU5YiIAAqUfuvk8VV8YNwgrn/qXTY36+p5EYmeAqWfMjO++6HD2dLSzs+ffjfqckREFCj92eHDyjh32ghuf34Jy9Y3R12OiGQ5BUo/d9kHJ5Cbk8PVj70ddSkikuUUKP3ckLICvnDSGP4ybxUvLl4fdTkiksUUKBngCyeOpbq8kB88/AYJXewoIhFRoGSAwvwY3zvzcN5e3agbR4pIZCIJFDM718zeMLNOM6tNah9tZtvM7LXwdWPSsmlmNs/MFprZ9WZmYXvczO4J2180s9ERdClyMyYN5bixlVz7xDts1KOCRSQCUY1Q5gMfB57pZtkid58Svi5Jar8BuBgYF75mhO2fAza6+6HAdcDV6Su77zIzfvDhiWxtTfCTJxdEXY6IZKFIAsXd33L3Hn/rmdkwoMzdn3d3B+4APhouPhu4PZy+Dzi1a/SSbSYMLeWCY0Zx54vLmLN0Q9TliEiW6YvHUA4xs1fN7O9m9oGwrRpIvl97fdjWtWw5gLsngM1AZXcfbGYXm1mdmdU1NDSkp/qIffP0CQwfUMg3/zBXT3YUkV6VtkAxs7+a2fxuXmfvZbVVwEh3nwp8A7jLzMqA7kYc3rWpvSzbudH9Jnevdffaqqqq/elOv1ESz+Wac47kvXVNXPuEdn2JSO/JTdcHu/tpB7BOK9AaTs8xs0XAeIIRyYikt44AVobT9UANUG9mucAAIKv39xx36CDOP2Ykv5n9HjMmDWX66IqoSxKRLNCndnmZWZWZxcLpMQQH3xe7+yqg0cyOCY+PXAg8FK72MHBROH0O8HR4nCWrXXHG4VSXF/LNP7xOU2si6nJEJAtEddrwx8ysHjgW+IuZPR4uOhGYa2avExxgv8Tdu0YbXwRuBhYCi4BHw/ZbgEozW0iwm+zyXupGn1Ycz+XacyezbEMz//fPb0ZdjohkAcvWX+Zra2u9rq4u6jLS7urH3uaGWYu48fyjmDFpWNTliEg/Z2Zz3L22u2V9apeXpN6lp43nfdUDuPyP81i9uSXqckQkgylQMlx+bg4/mzmF1vZOvnHva3R0ZueIVETST4GSBcZUlfCjj0zkH4vWc92T70RdjohkKAVKlji3dgQzp9fwi78t5JF5q6IuR0QykAIlS5gZPzp7IkeNLOeye1/nrVVboi5JRDKMAiWLxHNj3Hj+NMoKc7n4t3Vs0F2JRSSFFChZZnBZATeeP401W1q5+I46Wtp1vy8RSQ0FShaaOnIg131iCnVLN/Kt++bSqTO/RCQF0nYvL+nbzjxyGMs2HMbVj71NzcBCvj3jsKhLEpF+ToGSxS45aQzLNjTxy1mLGFwa59PHHxJ1SSLSjylQspiZ8R9nT2Ld1jZ++Kc3iefF+OTRI6MuS0T6KR1DyXJ5sRx+8ampnDS+iu8+MI8HXq3f90oiIt1QoAjx3Bi/umAax46p5LJ7X+eh11ZEXZKI9EMKFAGgIC/GzRfVMn10BV+/5zX+ULc86pJEpJ9RoMh2Rfm53PaZoznh0EF867653Pni0qhLEpF+RIEiOynMj/HrC2s55bDBXPnAfG78+yKy9Zk5IrJ/FCiym4K84BYtZx05jB8/+jbfe3A+iY7OqMsSkT5Opw1Lt/Jzc7h+5lSqBxbyq78vZtXmFn7+yakUx/VPRkS6pxGK7FFOjnHFGYfznx+dxKwFa/mXG/7B8g3NUZclIn2UAkX26fxjRnHrZ45m5aZtnP2/s3lh8fqoSxKRPkiBIj1y0vgqHvzS8ZQX5XH+zS9y6+z3dLBeRHaiQJEeG1NVwoNfOp6TJ1Txoz+9yRd+O4fNze1RlyUifYQCRfZLWUEev76wlu+deTh/W7CWD13/LHOWboi6LBHpAxQost/MjM9/YAz3XXIcOTlw7o3Pc+3jC2hL6NRikWymQJEDNrmmnEe++gH+5agR/OJvC/n4DbN5Z01j1GWJSEQUKHJQSgvyuObcyfzqgmms3NTCmdc/y3VPvkNrQo8WFsk2ChRJidMnDuWJS0/kzPcN42dPvcsZP3tWpxeLZBkFiqTMoJI4P505lds+M522RCczb3qBr9/9Kmu2tERdmoj0AgWKpNzJEwbz5KUn8dVTDuWR+as55dpZ3DBrES3t2g0mkskUKJIWhfkxvvHBCTx56YkcO7aSqx97m1OuncV9c+rp6NQFkSKZSIEiaTWqspibL5rOXf/2fgaVxvnmH17nzOuf5bH5q3WlvUiGiSRQzOwaM3vbzOaa2QNmVp607AozW2hmC8zs9KT2aWY2L1x2vZlZ2B43s3vC9hfNbHTv90j25bixg3jw34/n55+cSluik0t+N4ezfv4cf31zjYJFJENENUJ5Epjk7kcC7wBXAJjZEcBMYCIwA/ilmcXCdW4ALgbGha8ZYfvngI3ufihwHXB1b3VC9k9OjvHhycN54tIT+cm5k2lsSfD5O+o442fP8sCr9bTrmSsi/VokgeLuT7h7Ipx9ARgRTp8N3O3ure7+HrAQONrMhgFl7v68B7/O3gF8NGmd28Pp+4BTu0Yv0jflxnL4l2kjeOqyk/ifT0ym051L73mdk6+Zxa+fWczmbbo/mEh/1BeOoXwWeDScrgaWJy2rD9uqw+ld23daJwypzUBldxsys4vNrM7M6hoaGlLWATkwebEcPn7UCB772oncclEt1QMLueqRtzj2v57i+w/N11X3Iv1M2h6/Z2Z/BYZ2s+hKd38ofM+VQAK4s2u1bt7ve2nf2zq7N7rfBNwEUFtbqx33fUROjnHq4UM49fAhzF+xmVtnL+Hul5Zzx/NLmT56IJ96/0jOmDSMgrzYvj9MRCKTtkBx99P2ttzMLgLOAk71HUdl64GapLeNAFaG7SO6aU9ep97McoEBgG5/209Nqh7ATz4xme9+6DDuf6We37+0nEvveZ3vP/QGH548nHOmjWBqTTnaqynS90TygHAzmwF8BzjJ3ZOfKfswcJeZ/Q8wnODg+0vu3mFmjWZ2DPAicCHw86R1LgKeB84BnnadNtTvVZbEufjEsXz+hDG88N567qur54+v1HPXi8sYVVnEh48czkemDGf8kNKoSxWRkEXx3WtmC4E40HWzpxfc/ZJw2ZUEx1USwNfd/dGwvRa4DSgkOObyFXd3MysAfgtMJRiZzHT3xfuqoba21uvq6lLaL0mvxpZ2Hp23modfX8k/Fq2j02Hc4BLOmDSUM943jMOGlmrkIpJmZjbH3Wu7XZatv8wrUPq3tY0tPDpvNY/MW8XLSzbQ6VBTUchphw/htMOHMH10Bfm5feGcE5HMokDphgIlczQ0tvLEm6t56q21zF64jtZEJyXxXI4bW8lJE6o4cVwVNRVFUZcpkhEUKN1QoGSm5rYEsxeuZ9aCtcxa0MCKTduAYPRy/NhBHDu2kvcfUsnQAQURVyrSPylQuqFAyXzuzqKGrTz37jpmL1rPC4vX09gSXE87qrKIo0dXMG3UQGpHD2RsVYmOv4j0gAKlGwqU7NPR6by5cgsvvreeFxZvYM7SDWxsDq7KH1CYx+SacqbUlDOlZgDvqy6nqjQeccUifY8CpRsKFHF3Fq9rYs6SjbyybCOvLd/EO2sa6bq7/tCyAiZVD2Di8DKOGF7GEcPKGDGwUCMZyWp7C5RIrkMR6QvMjLFVJYytKuET04PraZtaE8xfsZl5KzZv//nU22vo+r2rJJ7L+CElTBhaxvghJYwfUsq4wSVUlcYVNJL1NEIR2YfmtgQLVjfy5qotLFjdyNurG1mwunGnm1iWFuRuD6cxVcUcMih4jaosoihfv7dJ5tAIReQgFOXnMnXkQKaOHLi9zd1p2NrKu2u28u6aRhY1NAUnACxs4P5X6ndaf3BpnFGVRYysKGZkRRE1FYXUVBRRXV7IkLICYjka2UhmUKCIHAAzY3BpAYNLCzj+0EE7LWtqTbBkfROLG5pYtqGZJeuaWLq+mdkL13H/lpad3psXM4YOKGD4gEKqBxYyfEBhMF9ewNCyQoYNKKC8KE+706RfUKCIpFhxPJeJwwcwcfiA3Za1tHdQv3EbKzZto35jM8s3bGPlpuD1/KL1rNnSsv2kgC7x3ByGlBUwpCwehFhZnKrScLo0zqCSYL6iOF+jHYmUAkWkFxXkxTh0cAmHDi7pdnmio5OGra2s3NTCmi0trNoc/Fy9uYW1jS28tWoLsxa00NTWsdu6OQYVxflUFsepLMmnsiROZXE+lcX5DEz6WVGcz8CifMqL8siL6fY0kjoKFJE+JDeWw7ABhQwbULjX9zW1JmhobKVhayvrwp8Nja2s29rG+q2trNvayrz6Tazf2kZja2KPn1Maz6W8OI+BRfkMKMyjvCif8sK8cDr4OaAwj7Kkn2UFuZTEc7UbTnajQBHph4rjuRTHcxk9qHif721NdLCpuZ31W9vY0NTGxuY2NjW3saGpffv0xuZ2Nm9rZ/mGZjY2t7OlpZ29nQCaY8Ep1GWFeZQW5FFakEtZQe726ZJ4MF1SkEtpPJgvCduD2mOUxHMpzIspmDKIAkUkw8VzYwwpizGkrOf3L+vsdBpbE2zZFgTN9p8t7WzZlmDztnYaW9ppbEkEbS0JVmxqYWtrI40tCRpbEnTsejCoGzkGxfm5FMVjFIfBU5QfC9tyKc6PUbh9PkZRXoyi/FwK82MUhcuK8oN1CvOC+cK8GAV5MR1PioACRUR2k5Nj23d31ez77btxd1raO2lsaWdrayJ4tQQ/m9oSbG3toKk1QVO4rLm1g61twXxzWwdrGltoWtdBc1uwrKktsdvJCvsSz83ZHjBdIVOQl7NT6HS1FeQmTefFiOfFKMjN2ek98dyk5bnBfDx3x3yOAkyBIiKpZ2bBF3d+jMEp+Dx3pzXRyba2IFy2tXXQHL62tSfY1tZJc1uCbe0d25e1JDpo2f6eDlraO2lpD0JqY1P79uUtiaB9W3vHXnfz7UtezLaHTDw3h/yu0MnL2R5A+bk55MeCZfm7vG/7fNLyvK7p2I73bm9PWha02U7Logg4BYqI9Hlmtn20MLA4Py3bcHfaOjppae+ktb2D1jBoWto7aU3sCKTWxI751sSO97UlOnebbk3smN7W3sHmbe3hfPJ7OmlLdNLW0ZnS/sRyLAiZ2I4QCl7G104bz0cmD0/p9kCBIiICBKEVjDBiUJjX69vvCrS2pIBpS3TS3pEUOru0t3X4julweXtnJ+0Jp62jg/ZweVtHJ4mOzmC+o5PyNPVPgSIi0gfsFGj9lK5qEhGRlFCgiIhISihQREQkJRQoIiKSEgoUERFJCQWKiIikhAJFRERSQoEiIiIpYX4wN6/px8ysAVi6H6sMAtalqZy+LBv7nY19huzsdzb2GQ6u36Pcvaq7BVkbKPvLzOrcvTbqOnpbNvY7G/sM2dnvbOwzpK/f2uUlIiIpoUAREZGUUKD03E1RFxCRbOx3NvYZsrPf2dhnSFO/dQxFRERSQiMUERFJCQWKiIikhAKlB8xshpktMLOFZnZ51PWkg5nVmNnfzOwtM3vDzL4WtleY2ZNm9m74c2DUtaaamcXM7FUz+3M4nw19Ljez+8zs7fDv/NhM77eZXRr+255vZr83s4JM7LOZ/cbM1prZ/KS2PfbTzK4Iv9sWmNnpB7NtBco+mFkM+F/gDOAI4JNmdkS0VaVFArjM3Q8HjgG+FPbzcuApdx8HPBXOZ5qvAW8lzWdDn38GPObuhwGTCfqfsf02s2rgq0Ctu08CYsBMMrPPtwEzdmnrtp/h//GZwMRwnV+G33kHRIGyb0cDC919sbu3AXcDZ0dcU8q5+yp3fyWcbiT4gqkm6Ovt4dtuBz4aSYFpYmYjgDOBm5OaM73PZcCJwC0A7t7m7pvI8H4TPPK80MxygSJgJRnYZ3d/BtiwS/Oe+nk2cLe7t7r7e8BCgu+8A6JA2bdqYHnSfH3YlrHMbDQwFXgRGOLuqyAIHWBwhKWlw0+BbwOdSW2Z3ucxQANwa7ir72YzKyaD++3uK4BrgWXAKmCzuz9BBvd5F3vqZ0q/3xQo+2bdtGXsudZmVgLcD3zd3bdEXU86mdlZwFp3nxN1Lb0sFzgKuMHdpwJNZMaunj0KjxmcDRwCDAeKzez8aKvqE1L6/aZA2bd6oCZpfgTBUDnjmFkeQZjc6e5/DJvXmNmwcPkwYG1U9aXB8cBHzGwJwa7MU8zsd2R2nyH4N13v7i+G8/cRBEwm9/s04D13b3D3duCPwHFkdp+T7amfKf1+U6Ds28vAODM7xMzyCQ5gPRxxTSlnZkawT/0td/+fpEUPAxeF0xcBD/V2beni7le4+wh3H03w9/q0u59PBvcZwN1XA8vNbELYdCrwJpnd72XAMWZWFP5bP5XgOGEm9znZnvr5MDDTzOJmdggwDnjpQDeiK+V7wMw+RLCvPQb8xt2virai1DOzE4BngXnsOJ7wXYLjKPcCIwn+U57r7rse8Ov3zOxk4JvufpaZVZLhfTazKQQnIuQDi4HPEPyCmbH9NrMfAecRnNH4KvB5oIQM67OZ/R44meAW9WuAHwAPsod+mtmVwGcJ/ly+7u6PHvC2FSgiIpIK2uUlIiIpoUAREZGUUKCIiEhKKFBERCQlFCgiIpISChSRg2Rm/wh/jjazT6X4s7/b3bZE+iKdNiySIsnXsuzHOjF379jL8q3uXpKC8kTSTiMUkYNkZlvDyR8DHzCz18Jnb8TM7Boze9nM5prZF8L3nxw+e+YuggtJMbMHzWxO+LyOi8O2HxPcHfc1M7szeVsWuCZ8tsc8Mzsv6bNnJT3r5M7wynCRtMuNugCRDHI5SSOUMBg2u/t0M4sDs83sifC9RwOTwluGA3zW3TeYWSHwspnd7+6Xm9mX3X1KN9v6ODCF4Fkmg8J1ngmXTSV4vsVKYDbBPcueS3VnRXalEYpI+nwQuNDMXiO4hU0lwb2SAF5KChOAr5rZ68ALBDfrG8fenQD83t073H0N8HdgetJn17t7J/AaMDoFfRHZJ41QRNLHgK+4++M7NQbHWpp2mT8NONbdm81sFlDQg8/ek9ak6Q70/1x6iUYoIqnTCJQmzT8OfDF8LABmNj58kNWuBgAbwzA5jOARzF3au9bfxTPAeeFxmiqCJzAe8F1iRVJBv7mIpM5cIBHuurqN4Lnto4FXwgPjDXT/iNnHgEvMbC6wgGC3V5ebgLlm9oq7/2tS+wPAscDrBA9E+ra7rw4DSSQSOm1YRERSQru8REQkJRQoIiKSEgoUERFJCQWKiIikhAJFRERSQoEiIiIpoUAREZGU+P+ZSjJOwy8ECgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 绘制学习曲线\n",
    "loss_plot(loss_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来探究下量子版本的奇异值分解的精度问题。在上述部分，我们提到过可以用分解得到的更少的信息来表达原矩阵。具体来说，就是用前 $T$ 个奇异值和前 $T$ 列左右奇异向量重构一个矩阵：\n",
    "\n",
    "$$\n",
    "M_{re}^{(T)} = U_{m \\times T} * D_{T \\times T} * V^{\\dagger}_{T \\times m}\n",
    "$$\n",
    "\n",
    "并且对于一个本身秩 (rank) 为 $r$ 的矩阵 $M$, 误差随着使用奇异值的数量变多会越来越小。经典的奇异值算法可以保证：\n",
    "\n",
    "$$\n",
    "\\lim_{T\\rightarrow r} ||M - M_{re}^{(T)}||^2_2 = 0\n",
    "$$\n",
    "\n",
    "其中矩阵间的距离测量由 2-norm 来计算,\n",
    "\n",
    "$$\n",
    "||M||_2 = \\sqrt{\\sum_{i,j} |M_{ij}|^2}\n",
    "$$\n",
    "\n",
    "目前量子版本的奇异值分解还需要很长时间的优化，理论上只能保证上述误差不断减小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABGDUlEQVR4nO3dd3iUZdbA4d9JIyGUUALEANIJJZBAFBUQsCAqICKgrAWsa1l3VxZX1F1ld79dG/Z17aJiQWxgb1QBAem9iATphJJAIIGQnO+PZxKSkDIhZZLMua9rrsy8M/POmRDmzNPOI6qKMcYY/xXg6wCMMcb4liUCY4zxc5YIjDHGz1kiMMYYP2eJwBhj/FyQrwM4HQ0bNtQWLVr4OgxjjKlSlixZsk9VI/Mfr5KJoEWLFixevNjXYRhjTJUiIlsLOm5dQ8YY4+csERhjjJ+zRGCMMX6uSo4RGFNZ7N27l7Fjx7J+/XqysrJ8HY4xBAQEEBMTw4QJE2jUqJFXz7FEYEwpjB07ln79+vH6668THBzs63CMISMjg0mTJjF27Fjefvttr57jN11DU5ftYPAjH7HwoR4MeuRjpi7b4euQTDWwfv16rrvuOksCptIIDg7m+uuvZ/369V4/xy9aBFOX7eD+T1bxgL7PWYEbGH7kPe7/JByAIfHRPo7OVGVZWVmWBEylExwcXKKuSr9oETzx7QZqZexjeOBsAkQZHjiHWhn7eeLbDb4OzRhjfM4vEsHO5DT+GPQpgbgMGUgmdwd9ws7kNB9HZkzpBQYGEhcXR+fOnRk0aBDJyck+i2XWrFnMnz+/zM43depU1q5dm3P7oYce4ocffiiz8+f32Wef8eijj3r12CNHjtCgQQNSUlLyHB8yZAhTpkwBXPxdunQhJiaGzp0789FHH+U8bsGCBfTo0YO4uDg6dOjA+PHjSUxMpGnTpqd8m4+Li2PRokWMHz+e6Oho4uLiaNu2LUOHDs3z+zldfpEIYuumMTxwNsGSCUCIZDI8cA6d66b7ODLjb6Yu20HPR2fQctyX9Hx0RpmMVYWFhbF8+XJWr15N/fr1eeGFF8og0tNTVCI4ceJEic+XPxH885//5KKLLjrt+IozePBgxo0b59Vjw8PD6d+/P1OnTs05lpKSwty5cxk4cCArVqxg7NixTJs2jfXr1/P5559z3333sWTJEgBGjRrFK6+8kvNvN2LECFq0aEGzZs348ccfc865fv16Dh8+zNlnnw3APffcw/Lly9m0aRNXX301F1xwAUlJSaV6336RCJ6J+h4h705sAWTxbNR3PorI+KPssaodyWkosCM5jfs/WVWmExfOPfdcduxw59u8eTMDBgyge/fu9O7dO2fwcM+ePVx55ZV07dqVrl275nxwP/XUU3Tu3JnOnTvzzDPPAJCYmEiHDh249dZb6dSpE/379yctzbWkn3vuOTp27EiXLl245pprSExM5KWXXuLpp58mLi6OH3/8kdGjRzNmzBj69evHfffdx/jx45kwYUJOvJ07dyYxMRGAt99+my5dutC1a1euv/565s+fz2effca9995LXFwcmzdvZvTo0TnfqqdPn058fDyxsbHcdNNNHDt2DHAlaB5++GG6detGbGxsgYOmPXr0YM2aNTm3+/bty5IlS3jzzTf5wx/+AMDnn39Ojx49iI+P56KLLmLPnj2nnGfkyJFMnjw55/ann37KgAEDqFmzJhMmTOCBBx6gZcuWALRs2ZIHHniAJ598EnBTj6OiogDXquvYsWOB55w8eTIjR44s8N/76quvpn///rz33nsF3u8tvxgsbpW+BiTvt5EacsIdN6YMXf3yT4Xet+y3ZI5n5m3yp2VkMv7zNQyJj+bAkePc8c6SPPd/8PtzvX7tzMxMpk+fzs033wzAbbfdxksvvUTbtm1ZuHAhd955JzNmzOCPf/wjffr04dNPPyUzM5PU1FSWLFnCxIkTWbhwIapKjx496NOnD/Xq1WPTpk28//77vPrqq4wYMYKPP/6Y6667jkcffZQtW7ZQo0YNkpOTiYiI4Pbbb6dWrVqMHTsWgNdff52NGzfyww8/EBgYyPjx4wuMfc2aNfz73/9m3rx5NGzYkAMHDlC/fn0GDx7MwIEDGTZsWJ7Hp6enM3r0aKZPn067du244YYbePHFF/nzn/8MQMOGDVm6dCn/+9//mDBhAq+99lqe519zzTVMmTKFf/zjH+zatYudO3fSvXt3Vq1alfOYXr16sWDBAkSE1157jccffzznQzzbgAEDuOWWW9i/fz8NGjRg8uTJ3H333TnvKfv3kC0hIYHnn38ecN/s27dvT9++fRkwYACjRo0iNDSUESNGEB8fz/PPP09QUBAffPABH374YaH/7t26dSvRDKGC+EWLgNvnwvgUGJ9C1oNJ7JN6/BxyljtuTAXJnwSyJR/NKNV509LSiIuLo0GDBhw4cICLL76Y1NRU5s+fz/Dhw4mLi+P3v/89u3btAmDGjBnccccdgPsmWrduXebOncuVV15JeHg4tWrVYujQoTndEy1btiQuLg6A7t2753yD79KlC9deey3vvPMOQUGFf6ccPnw4gYGBRb6HGTNmMGzYMBo2bAhA/fr1i3z8hg0baNmyJe3atQNcN8ucOXNy7h86dOgp8eY2YsSInA/XKVOmMHz48FMes337di655BJiY2N54okn8rQgsoWEhDB48GA++ugj9u3bx/Lly+nfvz8AqoqI5Hl87j3iH3roIRYvXpzzjX7AgAEANGnShE6dOjF9+nSWL19OcHAwnTt3LvR3URb7zvtFiyC3gOAQ1nT6K28sTWbMtmS6NovwdUimGinqG3zPR2ewo4AJCtERYQDUDw8pUQsgW/YYQUpKCgMHDuSFF15g9OjRREREsHz5cq/OUdSHSY0aNXKuBwYG5nQNffnll8yZM4fPPvuMf/3rXwV+UILrS88WFBSUZyA0PT095/Xzf2iebry5Yw4MDCxwbCI6OpoGDRqwcuVKPvjgA15++eVTHnP33XczZswYBg8ezKxZswptzYwcOZL/+7//Q1W54oorcqYTd+rUicWLF9OlS5ecxy5dupSEhISc261bt+aOO+7g1ltvJTIyMqdlkd091Lhx40K7hbItW7YszzlPR4W2CEQkQkQ+EpH1IrJORM4Vkfoi8r2IbPL8rFfecXQbeCuLg7rx9k8FVmQ1plzce0l7woLzfjMOCw7k3kval8n569aty3PPPceECRMICwujZcuWOd96VZUVK1YAcOGFF/Liiy8Crjvp0KFDnH/++UydOpWjR49y5MgRPv30U3r37l3oa2VlZbFt2zb69evH448/TnJyMqmpqdSuXZvDhw8X+rwWLVqwdOlSwH0obtmyJSemKVOmsH//fgAOHDgAUOj5YmJiSExM5JdffgFg0qRJ9OnTp0S/r2uuuYbHH3+clJQUYmNjT7k/JSWF6Gi3zuitt94q9Dz9+vVj06ZNvPDCC3k+tMeOHcsjjzyS0yJJTEzkmWee4d577wVcIs1OaJs2bSIwMJCIiAgArrrqKr766is++OADrrnmmkJf++OPP+a7774rNlkUp6K7hp4FvlHVGKArsA4YB0xX1bbAdM/tclU7NJgbOwfRYvWz7E9JLe+XMwZwixcfGRpLdEQYgmsJPDI0tkwXNcbHx9O1a1cmT57Mu+++y+uvv07Xrl3p1KkT06ZNA+DZZ59l5syZxMbG0r17d9asWUO3bt0YPXo0Z599Nj169OCWW24hPj6+0NfJzMzkuuuuIzY2lvj4eO655x4iIiIYNGgQn376ac5gcX5XXXUVBw4cIC4ujhdffDGna6dTp048+OCD9OnTh65duzJmzBjAfVg/8cQTxMfHs3nz5pzzhIaGMnHiRIYPH05sbCwBAQHcfvvtJfpdDRs2jMmTJzNixIgC7x8/fjzDhw+nd+/eOV1WBQkICOCqq65i//79nH/++TnH4+LieOyxxxg0aBDt2rWjXbt2vPjii7Rv7xL/pEmTaN++PXFxcVx//fW8++67OV1oERERnHPOOTRu3DhnsDlb9mB827Zteeedd5gxYwaRkafsNVMiUhb9S169kEgdYAXQSnO9qIhsAPqq6i4RiQJmqWqRX5ESEhK0tBvT7Fg0leivRjG765P0ufKWUp3L+K+EhATbJMl4Zdy4cSxcuJBvv/2WkJCQcn+9gv42RWSJqp7Sj1SRLYJWQBIwUUSWichrIhIONFbVXQCenwWWyxOR20RksYgsLu2cWYDohEEcrxXN+Yc+K/W5jDGmOI8++igzZ86skCRQUhWZCIKAbsCLqhoPHKEE3UCq+oqqJqhqQmmbQQAEBBJy9o3Iltmw75fSn88YY6qoikwE24HtqrrQc/sjXGLY4+kSwvNzb4VFFH89mRLId5O8W1JujDHVUYUlAlXdDWwTkez+/wuBtcBnwCjPsVHAtIqKidpN2NKoP6GSQUYhc7yNMaa6q+h1BHcD74pICPArcCMuGU0RkZuB34BTV3aUoza3v0+bEsxfNsaY6qZCE4GqLgcKWvlwYUXGkYcnCWzesJJGZ3agdqjVljfG+Bf/KDFRjN1z36b1+735YfYsX4diTIlZGeqyU5Iy1ABHjx7l2muvJTY2ls6dO9OrVy9SU1Pp27cv3377bZ7HPvPMM9x5550kJiYSFhZGfHw8HTp04Oyzzy5ywVpFsEQANIm/nOMEo4snkpVVMesqjB87vBsmXgqHT61meTqsDHXZKUkZanCL8xo3bsyqVatYvXp1zt7V+SuIQt4qoq1bt2bZsmWsW7eOyZMn8/TTTzNx4sQyfS8lYYkAILwBe5oO4KLjM5i3/jdfR2Oqu9mPw28LYPZjZX5qK0NdsWWod+3alVOGAqB9+/bUqFGDYcOG8cUXX+TElJiYyM6dO+nVq9cp52jVqhVPPfUUzz33XHH/vOVHVavcpXv37lrWjv86T/XhOvrGc/8o83Ob6uuUv8U3Ljv1svAVd9+xI6qvXKQ6PkL14Tru56sXqS59x92fuu/U53ohPDxcVVVPnDihw4YN06+//lpVVS+44ALduHGjqqouWLBA+/Xrp6qqI0aM0KeffjrnOcnJybp48WLt3Lmzpqam6uHDh7Vjx466dOlS3bJliwYGBuqyZctUVXX48OE6adIkVVWNiorS9PR0VVU9ePCgqqo+/PDD+sQTT+TENmrUKL388sv1xIkTBd7fqVMn3bJli65evVrbtWunSUlJqqq6f//+nOd/+OGHec734YcfalpamjZt2lQ3bNigqqrXX399zns688wz9bnnnlNV1RdeeEFvvvnmU35nTz31lD700EOqqrpz505t27atqqpOnDhR77rrLlVVPXDggGZlZamq6quvvqpjxow55TzLli3TyMhIPeecc/TBBx/M+X2rql522WU6depUVVV95JFHdOzYsaqqumXLFu3UqVOe8xw8eFBDQ0NPOX9pFPQ5CSzWAj5TrUXgEdziXPbVbE2npC/5bf9RX4djqquU3yC7wooqJJe+BWplqH1XhjouLo5ff/2Ve++9lwMHDnDWWWexbt06IO8GM0VtLgNlU0q6NPyuDHWhROCq17j19V8ZsSCRBy/v6OuITFV045eF33fsEKQnQ85ueeput/H0eYc3KPr5hbAy1IXHXBFlqLMT59ChQwkICOCrr76iQ4cODBkyhDFjxrB06VLS0tLo1q1bofEuW7aMDh06FPmeypO1CHJp2LobvTq14oOft5F2PNPX4ZjqZvbjoPkWLmpWmY0VWBnqii9DPW/ePA4ePAjA8ePHWbt2LWeeeSbgEkTfvn256aabimwNJCYmMnbs2JydzXzBEkE+f2i1mzey/sZXP5du6zdjTrF9EWQez3ss87g7XkasDLX3yqIM9ebNm+nTp0/O7yEhIYGrrroq5/6RI0eyYsWKU/YU2Lx5c8700REjRnD33Xdz4403lij+slRhZajLUlmUoS6M7liKvNqP+e3Hcd7I+8vlNUz1YWWoTWVVWctQVwkS3Q2NiuO8g5+dHNQzxphqzBJBASThJti7lt+Wz/R1KMYYU+4sERSk81UcCwxn6SdPsjsl3dfRmEosICCAjIwMX4dhTB4ZGRkEBHj/8W6JoCA1anH0vL8SlTCIiJpWhM4ULiYmhkmTJlkyMJVGRkYGkyZNIiYmxuvn2GCxMaWwd+9exo4dy/r16/PMjzfGVwICAoiJiWHChAk0apR359/CBottQVkRjh/ax8Iv3yAj7gYu6NDE1+GYSqhRo0a8/fbbvg7DmFKxrqEiBP36Pb03/JtZ337s61CMMabcWCIoQkCnIaQH1+Xs/dNYtT3F1+EYY0y5sERQlOAwJO53XBKwmE/nLPF1NMYYUy4sERSjRo9bCJZMaq2bzIEjx4t/gjHGVDGWCIrTsA1HonvTTHfxwc/bfB2NMcaUOUsEXgi/8WM+af4g7yzYSqZtZWmMqWYsEXgjqAajzjuTw8n7mL6ubPaZNcaYysISgZcuPjGbxaF38MWPZVcy2BhjKgNLBF4KPPNcgslkTIMFvg7FGGPKlCUCb9U7E2l7MS22fgyZVlfGGFN9VGgiEJFEEVklIstFZLHnWH0R+V5ENnl+1qvImEok4SZI3c2nH7zG4XRLBsaY6sEXLYJ+qhqXq/DROGC6qrYFpntuV05t+3M8/AwarH+XnxMP+DoaY4wpE5Wh6NwVQF/P9beAWcB9vgqmSAGBhAx7hY4hZ9AwurGvozHGmDJR0S0CBb4TkSUicpvnWGNV3QXg+dmooCeKyG0islhEFiclJVVQuAVo2ZuG0a0ByMi0ssPGmKqvohNBT1XtBlwK3CUi53v7RFV9RVUTVDUhMjKy/CL0xrafWf7kFfzxHZtBZIyp+io0EajqTs/PvcCnwNnAHhGJAvD83FuRMZ2WY4eIOzyL4I1fsu3AUV9HY4wxpVJhiUBEwkWkdvZ1oD+wGvgMGOV52ChgWkXFdNpa9eNE3RZcG/gD7yzc6utojDGmVCqyRdAYmCsiK4BFwJeq+g3wKHCxiGwCLvbcrtwCAgg6+yZ6BKzn50XzSc/I9HVExhhz2iosEajqr6ra1XPppKr/9hzfr6oXqmpbz8+qMS8z7lqyAkIYlPEtny3f6etojDHmtJUoEYhIgohc7enaye7uqQxTUCteeEPk3DtJrdWSN+cnompVSY0xVZNXiUBEGovIQlyXznu4bh6Ap4Anyym2Sk8u/gf1+93J2l2HWPrbQV+HY4wxp8XbFsHTwG6gAZB7msyHuEFfvzWkcwMGhi7nrfk2aGyMqZq8TQQXAg+qav6vvZuB5mUbUtUSvuZ9/svj7NmwgKPHT/g6HGOMKTFvE0EYUNCGvZFAetmFUwXFDkODwni761pqhvjncIkxpmrzNhHMAUbnuq0iEoirCTS9rIOqUsLqIZ2vosbaj9H0FLJsK0tjTBXjbSL4K3CriHwP1MANEK8FegL3l1NsVUfCTZBxhBeee5QvVu3ydTTGGFMiXiUCVV0LxALzge+AUNxAcbyqbi6/8KqI6G5oky6cLWsJCw70dTTGGFMiXndqq+pu4OFyjKXqEkGun8rZNeuDiK+jMcaYEvF2HcEfROS6Ao5fJyJ3ln1YVVB4AxAh9Wga83/Z5+tojDHGa96OEfwZ2FbA8UTgnrIKpspb9RE8FcOf35xJ8tGCJlkZY0zl420iaAoUtGJqu+c+AxDZnlonkhmos/ng54LypjHGVD7eJoLdQFwBx7sB1g+SrUksND2Lm0NnMumnRDJtKqkxpgrwNhG8BzwnIheLSLDn0h94Bni33KKrihJuIjpzO00PLWXm+sq/x44xxnibCB4G5gHf4moNHQW+xk0n/Xv5hFZFdboSDa3LzaEzeeunRF9HY4wxxfJq+qiqZgAjReQhIN5zeKmq/lJukVVVwWHIoOfY90sQP/60j81JqbSOrOXrqIwxplAl2o9AVTep6hTPxZJAYToN4aILLiYkMIBJP1lVUmNM5eb1gjIRuRpXhbQR+RKIqg4u47iqvMgjG3m10Uf8Yckwxl7Snlo1rCCdMaZy8nZB2RPAO0ALIBnYn+9i8tu3kT4HP+Lahr+wP/WYr6MxxphCefs19QZgpKp+VJ7BVCsxg6BmQ8ZFLoAGd/s6GmOMKZS3YwQBwPJyjKP6CQqBbtfDhq9J2vErW/Yd8XVExhhTIG8TwSvAKbWGTDG6jQLN5NPXH+WRr9b5OhpjjCmQt11DEcDvRORiYCWQkftOVf1jGcdVPdRvCbHDuTCgGZec38HX0RhjTIG8TQQdOdk1FJPvPqujUJSrXqO1r2MwxpgieLugrF9ZvaBni8vFwA5VHSgi9YEPcDOSEoERqnqwrF6vUsjK4peV83h+fS0eu6oLobZ5jTGmEinRgrIy8icgd4f5OGC6qrbF7X88zgcxla8FL9B66iCWrljG5yt2+joaY4zJw+tEICL9ROQVEflGRGbkvpTgHE2By4HXch2+AnjLc/0tYIi356syOg0FEe6sPZe3fkpE1XrTjDGVh7cLykbjiszVBvoCSUA9XBnqtSV4vWeAvwJZuY41VtVdAJ6fjQqJ4TYRWSwii5OSkkrwkpVA3Wik3aVcyQw27DjAsm3Jvo7IGGNyeNsiGAv8QVVH4mYM3a+q8bjVxqnenEBEBgJ7VXXJ6QSqqq+oaoKqJkRGRp7OKXwr4SZCjx/kihpLeWt+oq+jMcaYHN4mglbAD57rx4Dscpr/BUZ7eY6ewGARSQQmAxeIyDvAHhGJAvD8rJ5F/FtfABHNubHeCr5atYu9h9N9HZExxgDeJ4L9uG4hgB1AZ8/1BkCYNydQ1ftVtamqtgCuAWao6nXAZ8Aoz8NGAdO8jKlqCQiA66cSNvJNMjKVyYtsK0tjTOXgbSL4EejvuT4Ft1vZROB94PtSxvAocLGIbAIu9tyunhq0plXjCM5v25B3F24lIzOr+OcYY0w58zYR/AH3oQ/wCPAErjUwBbilpC+qqrNUdaDn+n5VvVBV23p+Hijp+aqUNVP536E/kHzoMN+t2ePraIwxxusFZQdyXc8CHiu3iKq7mg2olbKR17tvo2Pry30djTHGeD19NFNETpnWKSINRCSz7MOqxlr0ggZt6ZXyOfXDQ3wdjTHGeN01JIUcrwEcL6NY/IMIJNwI2xcxZ85MXpmz2dcRGWP8XJFdQyIyxnNVgdtFJPeagUCgN7C+nGKrvrqOhB/+QfDyN5mqt3BLr1YEBBSWa40xpnwVN0aQvbWW4AaFc3cDHccVibu97MOq5mrWh/7/R1y9tnzRppclAWOMTxWZCFS1JYCIzASGVruqoL7U47acBRhpxzMJCQog0BKCMcYHvBojUNV++ZOAiLQRkdDyCctPHPiVPd88zjn/+YHZG6vngmpjTOXn7ayh/4jIKM91EZEfgI3ALhHpUZ4BVmubZ9B4wb+JD/yVt+Zv9XU0xhg/5e2soWuBDZ7rlwJdgXOAt6nOK4HLW+wICA7n3gbzmL0xiV+TvKrfZ4wxZcrbRNAY2O65fhkwRVUXAc8D8eURmF8IrQNdhtPxwA80CDzCpAXWKjDGVLySFJ0703O9P5C9GU0Qha8xMN5IuAk5kcYD0Sv5aPF2jhw74euIjDF+xttE8DHwnoh8D9QHvvEcjwN+KYe4/EdUVzizJ+dGCYePneDTZTt8HZExxs94VWsIGANsBZoDf1XVI57jUcCL5RGYXxn1BVEidN4+l7d/SuTaHs0RsYaWMaZieFt07gTwZAHHny7ziPxRQAAC/D6uBnd/mcRPv+7nvNYNfR2VMcZPFNo1JCLdRCQg1/VCLxUXbjU27zkGzrqcVmFHmffLPl9HY4zxI0W1CBYDTXBbRy7G1RsqqL9CcXWHTGm07Y98/3c+P/83wi8Y7utojDF+pKhE0BJIynXdlKdGMdD8PMJXTYK+93A8C0KCvB3LN8aY01foJ42qblVVzXW90EvFhVvNJdwEB7fw7ReT6fPETNIzbKsHY0z582qwWERaAUOAVriuoF+Baar6a/mF5oc6DoZvGtD9wJe0jrybfhNmsTslnTMiwrj3kvYMiY/2dYTGmGqo2EQgIn/B7VMciBsvECASeExE7rOZQ2UoqAb8bgoLdtZhyee/kuZpEexITuP+T1YBWDIwxpS5IjuhRaQn8Dhus/pIVY1S1SZAI9x00ic8jzFlpWkCj8zYnpMEsqVlZPLEtxsKeZIxxpy+4kYj7wTeVtUH821gv19V7wfe8TzGlKGYQ/N4M/gxAsjKc3xncpqPIjLGVGfFJYJzgDeLuP9Nz2NMGWocHkDfwBVcETCXD0L+SSTJgBucGTNlOat3pPg0PmNM9VJcImiCGxguzGZcmQlThnoMuJa9Wo+/BH3IWbKBu4M+oUZQAL3bNOCb1bsZ+Pxchr80ny9X7uJEZlbxJzTGmCIUlwjCgGNF3H8cqOHNC4lIqIgsEpEVIrJGRP7hOV5fRL4XkU2en/W8C736uqJ7C1JaDyI6YD8BoowImsMzl0cx6ZZz+On+C/nb5R3YfSidu95byrJtyb4O1xhTxXkzffRyESmsLyKiBK91DLhAVVNFJBiYKyJfA0OB6ar6qIiMA8YB95XgvNVS21on829oIFy6/20gjrphwdzSuxU39mzJ/M37SDjT5c1HvlrHsRNZPDyooxWsM8aUiDeJ4PVi7ldvXsizOC17C65gz0WBK4C+nuNvAbPw90RweDesnXbyduZxWP4u1G4CPX4PoXUJDBB6t43MecjxzCyOZ2blJIGfEw/QvXk9AgIsKRhjiiaexcMV82IigcASoA3wgqreJyLJqhqR6zEHVfWU7iERuQ24DaB58+bdt26txguavxgDyya5BJAtIAiyTkBoBPT8k0sIIeF5nqaqiAird6Qw8Pm5tGwYzqhzz2RYQjNq1fC24rgxproSkSWqmpD/eIUWs1HVTFWNA5oCZ4tI5xI89xVVTVDVhMjIyOKfUJVtX5Q3CYBLAg3aQLOzYfo/4Nk4WPASnDj5uOzWQPsmtXluZDwRNYMZ//lazv3PdP75+Vq27j+CMcbkV6EtgjwvLPIwcAS4FeirqrtEJAqYparti3puQkKCLl68uCLCrJx+WwAz/g8O7YS7FkFg4d/2l29LZuK8LXy5cheZqlwY04gbe7bkvNYNbCzBGD9TWIugwhKBiEQCGaqaLCJhwHfAY0AfYH+uweL6qvrXos7l94kAQBWOHoDwBnAsFd4ZCgk3Q+wwCDi1KvieQ+m8s2Ar7y38jf1HjnN2y/p8cNs5lgyM8SOFJYKK7DiOAt7yjBMEAFNU9QsR+QmYIiI3A78BVozfGyIuCYAbXM44Cp/eBnOfgn4PQMwgCDjZ89e4Tih/6d+eu/q14fMVO0nLyEREyMpSXp7zK0PizyCqbpiP3owxxpd81jVUGtYiKEBWFqybBjP/A/s2QpMuMOozCCt6WcbqHSkM/u9cnhoRx5D4aE5kZhEYINZSMKYaqgwtAlOeAgKg05XQYTCs+hC2zHEzjAAObIH6Be8t1Dm6LrPv7UfjOqEAvPrjFr5atYsbe7bg8i5R1AiyzeeMqe68mjUkIvVE5FkRWSkiu0Vkb+5LeQdpSiAgELpeA0P+57qPDu2EF3rAW4Nh288FPqVZ/Zo5u6GdERHK0eMnGDNlBT0fnckzP2wk6XBRi8uNMVWdV11DIvI50Am34GsP+RaRqerL5RJdIaxrqAQy0mHxG/Djk3B0H7QbAP0ehKguhT4lK0v58Zd9TJy3hVkbkggJDGBg1yhu6tmSztF1KzB4Y0xZKtWsIRE5DPRR1aXlEVxJWSI4DcdSYdHLMO9Zd/2eNVCn+HqBm5NSeWt+Ih8t2c7R45mc1aIefx0Qw1kt6ldA0MaYslTaMYLNVPDiM1PGatSC3n9xU0y3zD6ZBBa9Cm0uKnQMoXVkLf55RWf+0r89Hy7expvzEzmW4SqeHjhynACBWRuSeOLbDexMTrNtNY2pgrxtEfQB/gaMBVarqk93VbcWQRk5vAee7QpZGRB/PZx/L9Qt+gM8M0sJELeK+V9frOWdBYmICOkZJ8thhwUH8sjQWEsGxlQypS0x8QuuJPVS4LiIZOa+lGWgpgLVbgx/XAbdR8Oyd+C5ePjmfrdQrRC5p5YOT2hKWHBQniQAblvNR75eR1ZW1ZuabIw/8rZFMAeoB7xEwYPFH5dLdIWwFkE5OLgVZj8O6z+Hu5e5xWqqbuZREVqO+7LQ8rO1Q4Po0rQusdERdGlaly5N6xIdEWZrFIzxkdKOESQAZ6vq6rINy1Qa9c6EIS9A+n8gtK5LApOGwJk94Zw7oEbtAp92RkQYOwrYSzkiLJjLu0SxcnsKr8/9lYxMly5u79OacZfGcPxEFrM3JnFWi3pE1Awpz3dmjCmGt4lgLVCnPAMxlUSoZ3roscMQUgtm/hsWvgS97oGzboHgvGUo7r2kPfd/soq0jJM9hGHBgYwf3ClnjODYiUw27D7Miu0pdIxyf0Ybdh/m1rcX89/fxTOwyxn8sjeVr1ftokuzCLpE16VeuCUHYyqKt11DA4DxuAHjVUBG7vtVtfBO5XJgXUMVaMcSV+l08wyo1QSu/wQad8rzkKnLdpR41lB6Riard6TQplEtImqG8PGS7fzlwxU59zerH0YXT5dSbNO6dI6uS53Q4HJ5i8b4i9KuI8g9Gpj7CYLbfKxC6xBYIvCBxHnw86tw5csQVAP2b4aIM4ssgV1Sh9IzWL0jhVXbU1i5PYWVO5LZduBkt1OryHDeu+UcmtQNZX/qMcJCAqkZYlVSjPFWaccI+pVxPKaqadHTXQBOHHMlK4JDoe/90GkoHNkLH90Iw950s5FOQ53QYM5r3ZDzWjfMOXbgyHFW7Uhh5bZk1u46RMNarsvomR828dmKnSx/6GJEhJ8276dmSCAxUbWtPpIxJVRsiyB7o3ngBlXdUCFRFcNaBD6mCuu/dOMHe9dCo05ugdrmGdD9Rhj4VLmHsGjLAbbsS+Xqs5oDcNmzP7J21yGCA4WYJnWIbVqXLtF16dI0graNaxEceHKm9Ol0ZRlTHZS2a2gv0EtVN5ZHcCVliaCSyMqCNZ/A9H9CsmcP6aBQ+NPK024VnK4dyWms3JbMyh0prNyezMrtKRxOPwFAjaAAOp1Rh8tio2hYq0aBg9u2AM74g9J2Db2F21Ly3jKNylRtAQFuR7TEubBskttXWbPg/auh6VmuddC4Y4WEEh0RRnREGJfGutIZWVnK1gNHWbk9OWfMYe/hY0ycl5gnCYBbAPfQNDcz+oyIMM6ICKVxndA8rQhjqjNvWwT/A64FtgBLcHsN51DVP5ZLdIWwFkElcni3K1NxIv3kMQkECXClK5qfCwk3uX0SgkN9F6dHUQvgcgsQaFQ7lDMiQj3JIYzLYqOIaxbB8RNZpB47Qb2awbY4zlQppW0RdMCVlwBole8+qyPgz2Y/7loBuQUEQuwIaBTjSmB/cit0nw+DnvFJiLkVtgAuqm4ok27uwc7kNHalpLEjOT3n+pqdh/hu7R7aNa5NXLMI1u46xJAX5vHaDQlc1LExy7cl897CrS5h1A3LaVWcERFGaLANXJvKz6tEoKo2a8gUbPsiyDye91jmcdi9wq1UPucuSJwDtTxjBrtWwvd/d1VQ218KgRW7NqCwBXD3DYihTaNatGlUq8DnqSqZntpJTeqE8tDAjnSKdovjdqekMXtjEnsPHyN/A7t+eAhnRIQSVTeMcZfG0DqyFrtS0tiZnE5sdN2cDYEKYwPbpiKUaBK2iIQCbXCtgM2qml7MU0x1d/vcou8PCIBWfU/ePrwb9v0CU653C9S63eAuEc3KNcxs2R+iJf1wFRGCAl03UJO6odzU62TZ7gGdoxjQOYrjJ7LYcyidHclpntbEyetb9x8h0NON9NWq3fzri7Us+/vFhASF8PrcLXy7endOKyIqIozoiFDW7z7Mc9M35RT125Gcxv2frMrzPowpC96OEQQD/wH+AITgFpIdA54HHlTVjCKeXuZsjKCKy8qETd+7bqNN30FoHRi7yS1U8wN7DqWzduch+raPRER4d+FWpi3fya6UNHYlp3OimKqt0RFhzBt3QQVFa6qT0k4ffQoYCYzDrSkA6A08AryrqmPLMNZiWSKoRg5uhT1rIOYytz7hvRHQrIfbH6GCp6BWBplZyr7UY+xMTuPK/80v9HG39GrJ3wZWzIwsU32UdrD4d8BNqvpVrmObRSQJeA23YY0xJVfvTHcBSE92s49m/AtmPQIdBrkZRy16F1sOu7oIDBAa13HTV6MLGdgODwnM+XWoKrdNWkKX6Lr0i2lEx6g6BAT4x+/KlB1vE0Fd3HaV+W0GIsosGuPfwurBqM9h3yZY8qbbLGfNpzBiEnQc7OvoKlxhA9v/vvLk4rfkoxnsPZTOk2v38OT3G2lYqwZ920fSt30kvdtGUjfMCvWZ4nnbNbQAWKKqd+U7/iIQp6rnenGOZsDbQBMgC3hFVZ8VkfrAB0ALIBEYoaoHizqXdQ35iYw0WPsZdBrixg8WvAi7VrhWQtOz/KKV4O2soX2px5izMYmZG5KYszGJlLQMAgOE7s3r0ad9JFfGR3NGRFgBr2D8SWnHCM4HvgJ2Aj/hZg2dC5wBXKqqxUwdARGJAqJUdamI1MYtTBsCjAYOqOqjIjIOqKeq9xV1LksEfmr2EzDvGTieCo07Q8KNbr1CqG2VkduJzCxWbE9m5vokZm3cy+odh5jy+3M5u2V9Nu45zOa9qVzQoZEV5/NDpUoEnhOcAdwFxOBmDa0F/qeqO08zoGnAfz2Xvqq6y5MsZqlq+6Kea4nAjx07DKs+gsWvw+5V0P5yGPmer6Oq1PYeSqd+eAhBgQE88vU6Js5NZPnDF1MzJIglWw9Qq0Yw7RrXslXSfqDUiaCMg2kBzAE6A7+pakSu+w6qar0CnnMbcBtA8+bNu2/durVigjWVkyrsWOrWKZwRD8nb4KOboPto6HQlhNT0dYSVUkZmFpuTUolp4lpRQ16Yx/JtyZxRN5Q+7RvRr30kPds0JLyG7fNQHZ1WIvD03xerJDuUiUgtYDbwb1X9RESSvUkEuVmLwJxi2yKYdhfs2+i22+z6O9d1FFlk49Lv7UpJY/aGJGZu2Mu8X/aTeuwEIYEBnNWyHv3aN6Jv+0haR1probo43USQRfG1hFRVvfr64FmY9gXwrao+5Tm2AesaMmVBFbbOcwvV1n4GqFuoVtOr7zN+7/iJLBZvPcCsDUnM2rCXjXtSAbixZwseHtQJVSU9I4uwEBtbqKpOdx1BUTWGBgB/Ak54GYAArwPrspOAx2fAKOBRz89p3pzPmFOIQIte7pKaBL/9dDIJfHQTRDR3XUf1WrhSF6XcUa26CQkKyNkh7oHLOrD94FFmb0yiTaSrv7Q56QiXPfsj/7u2Gxd1bExmlhIg5GktWG2kqqnEYwQi0g14DDgfeBn4l6omefG8XsCPwCrc9FGAB4CFwBSgOfAbMLy4riZrEZgSOXHcfehv+Mq1Gtpc6H7+OrPCdlSrDrYdOMpb8xO5sVdLoiPC+ODn33hh5mb6tY+kb0wj9h0+xkPT1timP5VYWcwaagn8GxgOfAI8oKoFLTIrd5YIzGlJ2eE20Pn5DTiyxx0LCoU7F0B4Q6hR27fxVTGzNybx9vxE5m3el1MYryBWG6nyOO0SEyLSAHgIuB2YB5yrqvYpbKqeutHQd5zrFsq9o9q0P8D2n12V1JjLof1lUCvS19FWen3aRdKnXSTpGZks3HKAUW8sKvBxOz1lMmZt2EtIYADntWlYkWEaLxSZCETkAeCvuBW/V6jqNxURlDHl5vBuWPG+SwLg9k7Y/jN0HQm/zoBN38Lnf3I7q13/CQTbatzihAYH0qddZKG1kbJXND/zwyZqhgTmJILb3l5McGAArSPDad2oFq0ja9EqMpyaITZ1taIV9xv/PyAN2A7cKSJ3FvQgVfW/QjCmaipoRzXU7ar2p5Vukdr6L2H/LyeTwHd/g+Bw11poEusXpS1OR2G1ke69xE0CnDj6LFKPuQSsqgSIsGZnCl+v3kXuytvREWG0igynTaNa9G7bkAtibDC/vBWXCN7GtqI01UlhO6ptX+Q+4KO6uEs2Vdi7Dn6ZDrMfdTOPYgZC7HCI7laxsVdyxW36Uy88hHrhIYCbafTS9d0BOHYik8R9R9mclMrmvanuZ9IRPvh5G6pwQUxjMjKzOOc/0/nzxe24/pwzSTueyY+bkmjTqBbN69ckKLDond5M0YpMBKo6uoLiMKZiFLejWn4icN3Hbjrqxq9da+Hn1yGklksEGelu9lGrvtaNhEsGJZ0hVCMokPZNatO+Sd7BelXl2AnXekvLyOSy2ChaNHArxjfuOcxtk5YAEBwonNkgnNaeVkTryJPdTLVDC6++alNdT/JJiYnSsllDxqeOHYbMDLdGYeO3bjOd4HA3LTVmILTr70pqm3KTnpHJ+t2Hc1oQv3h+bt1/NM8Ob/+7thuXxUbx2/6jzNywl4FdomhQqwZTl+0osBuruk91Le3GNMaYbLmnmbbqB9d94loK67+EdZ9BQBDcMd+Vt8jKcvWQTJkKDQ4krlkEcc0i8hzPyMzitwNHcxJD5zPqAvBz4gEe/mwN57eLpEGtGoz/LO96B3Ctjke+XsflXaII9rOuJmsRGFNWsrJgxxLYPB3O/6tLAF/+BXYucy2FmIEQ2c7XUfolVSXp8DEa1KpBYIDQYtyXhT42QKBxnVDOiAjjjVFnUbdmMKt3pLAv9Rh92kVW6bpL1iIwprwFBECzs9wlW6OOrkrq9H+4S8N2EH8d9PyT7+L0QyJCozqhObcLm+oaERbMDee1YMfBNHalpFEr1H1EvrvwN75ds5ulf78YgHEfr2TVjhSiI8KIrhdGdEQYTeuFER1Rk+h6YdSrGVylEoYlAmPK01k3u0vKDlfiYt3ncNBTQl0VfhgPrfrAmb0gKMSnofqTwqa6jh/cqcAxgr/0b8d15zTPud2iYTi7D6WzZd8R5v6yj6PH83Yz1QwJJL55BO/ecg4AX6zcSXiNIPq1bwS4FkpJE0V5Dm5b15AxFU3VzUY6mAj/OxcyjkKNutDuErdWoc1FUKPWycdbgbxyUVYfrKpK8tEMdiSnsf1gGjuS09hxMI2wkADuvSQGgEuenkOz+jV5bZTrlen56AxEyGlRNM1pWbgWRVTdUEKDT1Z5LavB7Uq1MU1pWSIw1UZGGmyeCeu/gA1fQ9oBuOY9lxCO7AcUZv4Hlky0AnlV2JFjJzh6PJPI2jVQVZ78biPbDh5lhydx7DmUnmdRHcDIs5vzyNBYAGIf/pbDx04t9FzSOk42RmBMZRQcBjGXuUvmCVc6u6nn/+nPr8GsR1zrQbNg+TvQ5z5rFVRB4TWCcnZ9ExHGXpJ3y5WMzCx2p6TntCZ2JKfRppFrFaYdzywwCcDJOk6lZYnAmMoiMAha9j55u+MVblxh13J3+8QxeLk3jFlvU1KrmeDAAJrVr0mz+qdusRoWEsgZEaHsTE4/5b7sOk6lZX9NxlRWYRGQtD7vsaP74Yhn+49Pfu+6jXavduMOptr66yUxhAXn3Rkudx2n0rJEYExlVVCBPAmA2Y+51sGhHTDnCXipJzzfDb5/yNVFMtXOkPhoHhkaS3REGIIbGyjLVdDWNWRMZVVUgbygGjD6C0jde3JF808vuG04G3WAowdg71pXTjvA9hiuDk6njpO3bNaQMdVF2kFX3qJGbVg8Eb74M9Rs6GYgdRwMLc63tQp+zmYNGVPd5S50FzvcjTGs+xxWfwxL34LQCPjTCnfcmFwsERhTHdWoBZ2udJfsUtk7l59MAh/f4rqZOgx2C9lsv2a/ZonAmOouOBTaX+ou2cIjYdVHsHYaBNaA1hdA91F5H2P8hs0aMsYfDXgE/rIebvzG1ULas9q1GMCtdv75dTi8x6chmopjg8XGGLcO4cQx13r45Qd45ypAoPk50GGQu0Q0L/Y0pnIrbLDYWgTGGFfGIthTprn1hXDHT9D3frcb27cPwDOxsGetu//E8cLPY6qkChsjEJE3gIHAXlXt7DlWH/gAaAEkAiNU9WBFxWSMKYAINO7oLn3vg/2bYdP3bn0CwFdjYfvPnpbCYGjcyT0nm1VLrXIqskXwJjAg37FxwHRVbQtM99w2xlQmDVrDObef/LBvfg6E1T+5qvm5eJj79MnHz34cflvgVkCbKqHCEoGqzgEO5Dt8BfCW5/pbwJCKiscYc5rifgc3fgl/2QiDnoX6reDQLndfyk5Y+qanWuq7NuBcRfh6+mhjVd0FoKq7RKSRj+MxxnirViR0H+0u2ZNOvn8Isjybp5xIh5d6Q8KN0GWEa1mYSqnKDBaLyG0islhEFiclJfk6HGNMbiJubGD953mPH9nruoj2/+Ju71kL85+HXSshK+vU8xif8HWLYI+IRHlaA1HA3sIeqKqvAK+Amz5aUQEaY7xUULXUwGDXGmjV193eMge++5u7XrMBtOjt7utyNYScWovfVAxftwg+A0Z5ro8CpvkwFmNMaRRWLXXXClctFdyg8z1rYchL0OZi2LYQvrn/ZIXUtdPciudUa/VXpApbUCYi7wN9gYbAHuBhYCowBWgO/AYMV9X8A8qnsAVlxlQTqnBoJ9T1lFd+cyAk/uiuN+rkWgvt+p9sUZhS8Xn1UVUdWchdF1ZUDMaYSkbkZBIAuGGa25rz19nw6yxY/Dokbz2ZCH56AaLioOlZVlK7DPl6jMAYY04KCITo7u7Se4yrnJrmWWN6ZJ8bX9AsCK7pNt1p1cctbKvfyrdxV3GWCIwxlVdwKARHuevhDeGvWyBxLmzxtBi+f8htvlO/FaRsh03fQcs+7nbu1c6mSJYIjDFVR1gEdBjoLuDGF0JqueubZ8IX97jrdZtDq/OhZV+IuQxCwn0QbNXh61lDxhhz+uqcAaF13PX46+APi+GyCXBGV7c72ye3uO4lcFNXN3wN6YdOPv/wbph4qd+vgLYWgTGmehCBhm3d5exb3QrnpA0Q3sDd/9MLsPEbEM84RKs+boFbdl2kgU/5Nn4fskRgjKmeAgJdBdVsw99y6xa2zHazkuZM8Nyhri5Sw3ZuNlJUVwj0r49G25jGGOOfpt4FKz+ArAwIDIHME0AW1KgLLXpCy/Pdfs7VaEaSbUxjjDHZDu+G1R+5JABuBXRQCAx8GjoNgb3r4JtxbpwB3BTWxRPd3gxV8Mtzcfyr/WOMMVBwXSTNgt2rYfBz7nbyb269AsDWn+CLP7vrdZq61kLL86H9pW4mUxVnLQJjjP8prC7S9kUnb0c0d2sXwH3g/2EJXP4UNO3uBp2n3u4WuQFsWwRrPj15u4qxFoExxv/cPrdkjxeBhm3c5aybXQntvWtP7rGw9G1YNsldbxx7ssXQ7pIqsbDNEoExxpRUQAA06Xzy9sCnodsoNyNpyxxXI2nTt9Deszvvig+gViO3zWdwmG9iLoIlAmOMKa3AYGh2lrucP9YtYju0w92nCt89CEeS3OykZj08rYUBENXFt3F7WCIwxpiyFhx6sttIBP64zC1c+3WWazHM/A8cO+wSwYnjsOhlt0lPk9iTezNUIEsExhhT3mrUhrYXuwvA0QOQ6Zm6unvlyV3bQiOgRS9XOK/jFVC7cd7zHN4NH90Iw9489b5SsFlDxhhT0WrWP/lB3jQBxqyHoa+6Ynq7VsLX98LBLe7+XSvcYPTBrW7aa3ZJjDJkLQJjjPG1OlFub+cuI9ztg4lQx7Nhz7rPYc4T7roEuvUOy9+FPveVWavAWgTGGFPZ1GvhBqAB+j0Idy6EZueenIqqWWXaKrBEYIwxlZmIW728axlknXDHMo+7VkEZlc+2RGCMMZVdYSUxyqhVYInAGGMqO29KYpSCDRYbY0xlV9KSGCVkLQJjjPFzlgiMMcbPWSIwxhg/Z4nAGGP8nCUCY4zxc1Vy83oRSQK2nubTGwJVaRuhqhRvVYoVqla8VSlWqFrxVqVYoXTxnqmqkfkPVslEUBoislhVE3wdh7eqUrxVKVaoWvFWpVihasVblWKF8onXuoaMMcbPWSIwxhg/54+J4BVfB1BCVSneqhQrVK14q1KsULXirUqxQjnE63djBMYYY/LyxxaBMcaYXCwRGGOMn/ObRCAib4jIXhFZ7etYiiMizURkpoisE5E1IvInX8dUFBEJFZFFIrLCE+8/fB1TcUQkUESWicgXvo6lOCKSKCKrRGS5iCz2dTxFEZEIEflIRNZ7/n7P9XVMhRGR9p7fafblkIj82ddxFUZE7vH8/1otIu+LSGiZndtfxghE5HwgFXhbVTv7Op6iiEgUEKWqS0WkNrAEGKKqa30cWoFERIBwVU0VkWBgLvAnVV3g49AKJSJjgASgjqoO9HU8RRGRRCBBVSv9oicReQv4UVVfE5EQoKaqJvs4rGKJSCCwA+ihqqe7WLXciEg07v9VR1VNE5EpwFeq+mZZnN9vWgSqOgc44Os4vKGqu1R1qef6YWAdEO3bqAqnTqrnZrDnUmm/YYhIU+By4DVfx1KdiEgd4HzgdQBVPV4VkoDHhcDmypgEcgkCwkQkCKgJ7CyrE/tNIqiqRKQFEA8s9HEoRfJ0tSwH9gLfq2pljvcZ4K9AVjGPqywU+E5ElojIbb4OpgitgCRgoqfb7TURCfd1UF66Bnjf10EURlV3ABOA34BdQIqqfldW57dEUImJSC3gY+DPqnrI1/EURVUzVTUOaAqcLSKVsvtNRAYCe1V1ia9jKYGeqtoNuBS4y9PNWRkFAd2AF1U1HjgCjPNtSMXzdGENBj70dSyFEZF6wBVAS+AMIFxEriur81siqKQ8fe0fA++q6ie+jsdbnq6AWcAA30ZSqJ7AYE+/+2TgAhF5x7chFU1Vd3p+7gU+Bc72bUSF2g5sz9Ua/AiXGCq7S4GlqrrH14EU4SJgi6omqWoG8AlwXlmd3BJBJeQZfH0dWKeqT/k6nuKISKSIRHiuh+H+aNf7NKhCqOr9qtpUVVvgugNmqGqZfbMqayIS7pkwgKebpT9QKWe+qepuYJuItPccuhColBMc8hlJJe4W8vgNOEdEano+Hy7EjR2WCb9JBCLyPvAT0F5EtovIzb6OqQg9getx31azp7Zd5uugihAFzBSRlcDPuDGCSj8ts4poDMwVkRXAIuBLVf3GxzEV5W7gXc/fQhzwH9+GUzQRqQlcjPuGXWl5WlkfAUuBVbjP7jIrNeE300eNMcYUzG9aBMYYYwpmicAYY/ycJQJjjPFzlgiMMcbPWSIwxhg/Z4nAj4nILBH5r49eW0VkmC9e2xuVPb6y4KliOd6Lx80UkRsqIKSCXrvIfwdP1duhFRlTdWSJoJryLPL6n6eE8TER2SMi00Xk4lwPGwrc76sYy5qIdPd8cPQq5P4pIjKvouMqSmEfdCLyZmUokS0ilwPNgHdzHUv0xK0ikuYpOX2vZ6FTRfsX8JiI2GdZKdgvr/r6GFeK4GagHTAQ+BpokP0AVT3gqW5a5YhIUP4PHk/9oGW495z/8Q1w9WRer5gIq40/AW+qama+4//ELSTsgCuG9h/AFwXxvgJq48pEmNNkiaAa8pR76A2MU9XpqrpVVX9W1QmqOjnX4/J0DXm+6f1NRF72bNKxXUTuzXfudiIyW0TSRWSDiFwmIqkiMtpzfwvPN8WEfM8rron/qOd8aZ44Hs+98YaIjPd0ZYwWkc3AMaCgypavA8M9Bftyuw7IAD4QkQEi8qOIHBSRAyLyrYh0KCI2r96TiESLyGTPeQ+KyJci0raw85aEiAwVkZWe388Bz79B41z3DxJXnTRdRLaIyL/FFVPLvr+RiEzzPH+riNzkxWtG4sqFfFbA3YdVdbeqJqrqa8BKXPmL7Oe29rzebhE5IiJLxRX8y33+Yv/eCojpPhHZJyI9wBU7xCWDkcW9H1M4SwTVU6rnMlhKvovRPbgl7N2Ax4DHxbPLlKf5/SlwAjgHGA08DNQog5iPADfhvmHeiasD9GC+x7QEfgcMB7oC6QWc510gELg63/GbgMmqegSXQJ7BtZj6AinA57k/OEtKXKmCmZ6Y+gDn4soF/+C577SJSBNcgby3cL+f84FJue6/BPe+/wt0wr3XYeQt7/Am0Ab3wT4EuAFoUcxL98Il3DVFxCYi0tcTV0auu2rhWqAX4/6tPgY+EZGYfKco9O+tgNeZgCth0SdfmfNFuN+5OV2qapdqeAGuwm3Ek46rsTQBt/tS7sfMAv6b63Yi8H6+x2wC/ua5fgkuCUTnuv88XL380Z7bLTy3E/KdR4Fhhd0uIP7bgV9y3R6P+6Bp7MV7fweYn+v2WZ7X61HI48OBTKBXQfF5855wH76b8JRt8RwLBPYDI4qItcDfA+6D+wvP9W6ex51ZyDnmAH/Pd2wI7suA4LoGFVfOOvv+Mz3veXwRsf0Z2FrA8URcgkgFjnvOnQacV8y/y4LsvyVv/t5y/X6uBiYCG4EWBZx3MG5viSBf/F+rDhdrEVRTqvoxrm75INw3s/OABSLyQDFPXZnv9k6gked6DLBT3SYZ2X6mDDZ4EZFhIjLX05WQCjwNNM/3sO3qXang14Fzc337vAlYrZ5vkZ5ui/dEZLOIHAL24FrH+V+vJLrjWiyHPV1lqbiWRj2gdSnOC7AC+AFYLSIfi8gdnm6b3K/9YPbrel77PVyCa4L7tp6F++YMgLqduIrb4SqMgltdAE/hisr1wbWE/qGq87PvFFc19XERWevpJkvFbQ2a/3dc1N9btgm4llsvVU0sIJY0XMIrsz18/Y0lgmpMVdNV9XtV/aeqnof7gBxfTBdIRr7bysm/E6H4LSizk0LOQK64vRUKJSLn4Lo+vsUlrnjgb7gtL3M7UsxrZ5sF/ALcJK4s9kjyDhJ/DkQCvwd6eF7vBFDY78Wb9xQALMd9OOa+tANeLiLWw0DdAo5H4BIJ6vrB+3suK3GD4ZtEpGuu1/5HvtftArTF7Rh2urN59uESWUH2q+ovqvoTrvU5VkT65bp/Aq4L7++4ZBGHS0T5f8dF/b1l+x6X0AqrwFsfSNeT26WaEgrydQCmQq3F/ZuH4pr0JbUOiBaRM9SzWQruW17u/7hJnp9RuY7FFXPensAOVf1X9gEROfM04gPcHsoi8gZuxst63DfbSZ7zNsB9Q75LVWd6jnWj6P8L3rynpbiEs09Ltk/vBtw3+pxEJW4j9a647pCc94Tr4vtJRP6J67e/GtdaWArEqOovBb2AiKzD/RudBcz3HGuOazEWZRkQKSINVXVfYQ9S1YPiJh08LSLxnlh7AW97WqZ4xqpa47p3SuorXJnoD0VEVfWtfPd3xv0OzGmyFkE1JCINRGSGiFwnIl1EpKWIDMft0ztdT3/by+9xH1xviUhXzzf5p3Dfpl2Hrmoari/4PhHpJCLn4b4dFmUjLsFcKyKtROQOSj8L5E2goee1p6rqfs/xg7hvureKSBsR6QO85HkPBfLyPb2L62KaJiJ9PL/z80XkyWJmDj2Fa7ncJW5GVhyuznx9z09E5BzP7JqzPB/gg3Fz+7M3ffkn8DsR+aeIdBaRGE9X2+Oe+DcA3wAvi8i5ntd4E9elUpRluD2oC1yXkc8LQHtcKwDcv+mVItJNRGJx4zan3XWjbn+L4cBLcuritt6492dOkyWC6ikV98H1J2A27tvjf3D9xvln03hNVbOAK3GzhBbhZrH8G5cEcvclZ09N/BnXLfK3Ys77OfAEbibPStxMk4dON07POXfhvknWA17L9x6uxnWdrMZ9gP0dN/hZlCLfk6oexc3m+RW39+163O+nHi75FBbn+8CNnsti3AdaE6C3uh2/wHUR9QS+wA2mPgn8S1Xf8ZzjW+ByoB/u32URbq/g33K91GhgCzAD1zX2Hm6wtlCeLqk3gGuLepznsUm4Vtd4z+yyMbgk8iNujGqB5/pp8ySDEbiEdgO4Kbu48a+JRT3XFM02pjGl4umnXo6bUVOVNoQ3XhCRRriWx9mq+quv48lPRJ4A6qqqLxazVRs2RmBKRESuxA3absJNq3yKk/3UpppR1b3iFp81w7V2Kpu9FN/1aIphLQJTIp4m+d9wHwwHcTN07vFyWqcxphKyRGCMMX7OBouNMcbPWSIwxhg/Z4nAGGP8nCUCY4zxc5YIjDHGz/0/eYjewdb4YeMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "singular_value = singular_value_list[-1]\n",
    "err_subfull, err_local, err_SVD = [], [], []\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "\n",
    "# 计算 2-norm 误差\n",
    "for i in range(T):\n",
    "    lowrank_mat = np.matrix(U[:, :i]) * np.diag(D[:i])* np.matrix(V_dagger[:i, :])\n",
    "    recons_mat = np.matrix(U_learned[:, :i]) * np.diag(singular_value[:i])* np.matrix(V_dagger_learned[:i, :])\n",
    "    err_local.append(norm(lowrank_mat - recons_mat)) \n",
    "    err_subfull.append(norm(M_err - recons_mat))\n",
    "    err_SVD.append(norm(M_err- lowrank_mat))\n",
    "\n",
    "# 画图 \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(list(range(1, T+1)), err_subfull, \"o-.\", label = 'Reconstruction via VQSVD')\n",
    "ax.plot(list(range(1, T+1)), err_SVD, \"^--\", label='Reconstruction via SVD')\n",
    "plt.xlabel('Singular Value Used (Rank)', fontsize = 14)\n",
    "plt.ylabel('Norm Distance', fontsize = 14)\n",
    "leg = plt.legend(frameon=True)\n",
    "leg.get_frame().set_edgecolor('k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 案例2：图像压缩\n",
    "\n",
    "为了做图像处理，我们先引入必要的 package。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 图像处理包 PIL\n",
    "from PIL import Image\n",
    "\n",
    "# 打开提前准备好的图片\n",
    "img = Image.open('./figures/MNIST_32.png')\n",
    "imgmat = np.array(list(img.getdata(band=0)), float)\n",
    "imgmat.shape = (img.size[1], img.size[0])\n",
    "imgmat = np.matrix(imgmat)/255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUb0lEQVR4nO3db2yVVZ4H8O/X0tLaEmiLQgvlj0TEOnGANGgCGnfdGdE3otnR8cXEScwyLzRZk5kXxs3suPti192sTtxkY4KrGWbjOrirRndjdseQ3bAmxp3KQgVZAbFAaSkgBQqUQtvfvrgPs4V5fqe39z733rbn+0ma3p7fPfeePtwfz+3zu+ccmhlEZOa7odIDEJHyULKLRELJLhIJJbtIJJTsIpFQsotEQskuEgkluxSN5A9JjpI8P+7rvkqPS641q9IDkBnjEzPbUOlBiE9n9hmOZDfJn5DsInmW5DaStZUel5Sfkj0OjwHYCGA5gDsB/DDtTiQ3kDwT+AqdudeQPEVyP8mfktS7xilG/yBx+Fsz6wUAkv8CYHXanczsYwDzCnj8HQC+BeAwgDsAbAMwAuAvC3gsKRGd2eNwfNztiwAasnxwMztkZl+b2ZiZfQ7gzwH8YZbPIcVTsstvkbznuivq13/dk+dDGQCWcqwyeXobL79lZv+FAs76JB8EsNPM+kmuAvBTAP+U9fikODqzSxbuB9BF8gKADwG8C+AvKjskuR61eIVIHHRmF4mEkl0kEkp2kUgo2UUiUdbSW2Njo7W2tpbzKTN1ww3p/zd67QBA+uXm0MXR0dHRgmKFXHAtdPyh2NjY2KTagfDYZ83yX6pVVVVuzFPoOKa63t5eDAwMpP7DFJXsJDcCeAVAFYC/N7MXQ/dvbW3Ftm3binnKipo9e3Zqe319vdunpqbGjV2+fNmNnT171o2dO3fOjQ0PD6e2hxK6rq7OjXm/MxBOwKGhodT2wcFBt8/IyIgba25uLijmJe7FixfdPleuXHFjU93jjz/uxgp+G0+yCsDfAXgQQDuAJ0i2F/p4IlJaxfzNvg7AweRz0ZcB/ArAw9kMS0SyVkyyLwJwdNzPPUnbNUhuJtlJsnNgYKCIpxORYhST7GkXAX7nDyQz22JmHWbW0djYWMTTiUgxikn2HgBt435eDKC3uOGISKkUk+y/AXAryeUkawB8H8AH2QxLRLJWcOnNzEZIPgPg35Ervb1hZnszG5mIZKqoOruZfYjclEYRmeL0cVmRSCjZRSKhZBeJhJJdJBJKdpFIKNlFIqFkF4mEkl0kEkp2kUgo2UUioWQXiYSSXSQSSnaRSCjZRSKhZBeJhJJdJBJKdpFIKNlFIqFkF4mEkl0kEkp2kUgo2UUioWQXiYSSXSQSSnaRSBS1IwzJbgCDAEYBjJhZRxaDEpHsFZXsid8zs1MZPI6IlJDexotEothkNwC/JvkZyc1pdyC5mWQnyc6BgYEin05EClVssq83s7UAHgTwNMl7r7+DmW0xsw4z62hsbCzy6USkUEUlu5n1Jt9PAHgPwLosBiUi2Ss42UnWk5xz9TaA7wLYk9XARCRbxVyNXwDgPZJXH+cfzezfMhmViGSu4GQ3s0MAvp3hWESkhFR6E4mEkl0kEkp2kUgo2UUikcVn46Nx9uzZ1PZjx465fS5duuTGZs+e7cbmzp3rxm688UY3Njo6mtre39/v9jl//vykHw8ALl++7MYGBwcnPY6enh43VlVV5cba29vd2KpVq1LbFy9e7Papq6tzY9OZzuwikVCyi0RCyS4SCSW7SCSU7CKRKPvVeDMr91Nm5vjx46nte/fudft0d3e7sdraWjfW1NRUUD/vynpvb6/b5+TJk27Mq0AAwNDQkBvzxpjMpUh15MgRN3bo0CE3tmzZMje2adOm1PYHHnjA7bNo0SI3Np3pzC4SCSW7SCSU7CKRULKLRELJLhIJJbtIJDQRZhIuXryY2n706FG3T1dXlxs7d+6cGwtNCgmVw7yS10033eT2aWhocGOhSTIjIyNubPny5ZNqn2gcoWMcinkTb4aHh90+M5XO7CKRULKLRELJLhIJJbtIJJTsIpFQsotEoqylN5Koqakp51Nmqq2tLbX9nnvucfssWbLEjXmz6ABgzx5/J62+vj43tnDhwtT29evXu33uuuuuST8eANxwg3+u8MqD+/fvd/vs2LHDjd1xxx1uLFTO27BhQ2p7S0uL22c6v0ZDswonPLOTfIPkCZJ7xrU1kfyI5IHku7ZnFZni8nkb/wsAG69rew7AdjO7FcD25GcRmcImTHYz2wHg9HXNDwPYmtzeCmBTtsMSkawVeoFugZn1AUDy/WbvjiQ3k+wk2TkwMFDg04lIsUp+Nd7MtphZh5l1NDbqT3uRSik02ftJtgBA8v1EdkMSkVIotPT2AYAnAbyYfH8/347TecFJb9ul0AKF8+bNc2Nr1651Y48++mje4xrP227qypUrbp/QGOfPn+/G6uvr3djY2Fhq++nT11/++X8HDhxwY6EZdqHy5oIFC1Lbq6ur3T7T+TUakk/p7S0AnwC4jWQPyaeQS/LvkDwA4DvJzyIyhU14ZjezJ5zQ/RmPRURKSB+XFYmEkl0kEkp2kUgo2UUiob3eJsGbURQq48yePduNhfZsCy0QGZqVdeHChdT20CKVoTHOmuW/REZHR92Ytw9caJHNUMwr5QHh2XehmGc6v0ZDdGYXiYSSXSQSSnaRSCjZRSKhZBeJhJJdJBIqvWUgVN4Jla6qqqrcWKjUFFJXV5faHhpjaByh8lpoMZLe3t7U9kL2qQP8GYdAeNZeqCzqmYmvUUBndpFoKNlFIqFkF4mEkl0kEkp2kUiU9Wq8mRV8lXkqC13NDk1aCfULrbkWOobe1fg5c+a4fUJX3ENXzw8fPuzGDh48OOnHC613F1rnz1tnDvCPR8h0fo2GKgk6s4tEQskuEgklu0gklOwikVCyi0RCyS4SCU2EmQRv7N7adEC4vBbqF9quKVSW8yZ+NDQ0uH28LaMAoKenx4199dVXbuzQoUOp7aF15pqbm93YsmXLCuoXWl/PM51foyH5bP/0BskTJPeMa3uB5DGSu5Kvh0o7TBEpVj5v438BYGNK+8/NbHXy9WG2wxKRrE2Y7Ga2A4C/9aaITAvFXKB7hmRX8ja/0bsTyc0kO0l2hhY7EJHSKjTZXwWwAsBqAH0AXvLuaGZbzKzDzDoaG93/E0SkxApKdjPrN7NRMxsD8BqAddkOS0SyVlDpjWSLmfUlPz4CYE/o/jOFVyoLlddCsVCJJ1QOC60n55XlQn2Gh4fd2LFjx9zYF1984ca80ltoLbmVK1e6sTVr1rix0Gw5b9ZhIdtCTXcTJjvJtwDcB2A+yR4APwNwH8nVAAxAN4AflW6IIpKFCZPdzJ5IaX69BGMRkRKK772MSKSU7CKRULKLRELJLhKJss96C830muq8UlmohBb6fUP9QjPbClmMMrSI4tDQkBs7deqUGwuV5fr7+1Pb29ra3D5NTU1u7Oabb3ZjoUUlvd879O8ynV+jITqzi0RCyS4SCSW7SCSU7CKRULKLRELJLhKJspbeSE7r2UaFlN5CQsei0OPkldFOnjzp9jl69Kgb6+7udmPHjx93Y5cvX05tD+0519ra6sZC+7mFSpGhsqJnOr9GQ2XD6ftbicikKNlFIqFkF4mEkl0kEkp2kUhoIkwGCp0IE7rqG9q2aNYs/5/t/Pnzk2oHgC+//NKNha7Gh7Zyqq2tTW0PTWhZuHChG5s7d25B4wit5eeZia9RQGd2kWgo2UUioWQXiYSSXSQSSnaRSCjZRSKRz44wbQB+CWAhgDEAW8zsFZJNALYBWIbcrjCPmVmU27SGSm+hWGhrKG/bIiBcGjpz5kxq+9dff+322b9/vxsLTXYJlQC9EltLS4vbJzRJJmR0dNSNFbIG3UyVz5l9BMCPzex2AHcDeJpkO4DnAGw3s1sBbE9+FpEpasJkN7M+M9uZ3B4EsA/AIgAPA9ia3G0rgE0lGqOIZGBSf7OTXAZgDYBPASy4upNr8t3/aJSIVFzeyU6yAcA7AJ41M//zib/bbzPJTpKdp0+fLmSMIpKBvJKdZDVyif6mmb2bNPeTbEniLQBOpPU1sy1m1mFmHaFNAESktCZMduYuW74OYJ+ZvTwu9AGAJ5PbTwJ4P/vhiUhW8pn1th7ADwB8TnJX0vY8gBcBvE3yKQBHAHyvJCOcQrwyWmhrpZBQ6So0y8tb3w3wZ6l1dXW5fXbu3OnGTpxIfcMGILxdU3t7e2r7bbfd5vZpaGhwY6HZa8PDw27MK8uFyp4z1YTJbmYfA/CKkvdnOxwRKRV9gk4kEkp2kUgo2UUioWQXiYSSXSQSZV9wstCtkqYCr8QWKr2Fft9Qv1Bp6MqVK27s2LFjqe27d+92+4QWnAzNRFu1apUb27BhQ2r7nXfe6fYJLbLpzeYDgMHBQTfmCS32OZ1foyE6s4tEQskuEgklu0gklOwikVCyi0RCyS4SibKW3sys4BliU4G3SGFo8cJQLFTiGRkZcWNDQ0NuzCtDhcpToXE0Nze7sZUrV7oxb3bbggUL3D6hPdtCe9WFZgFWV1enthdaEp3qQr+XzuwikVCyi0RCyS4SCSW7SCSU7CKRKPtEmJmo0EkVoavI33zzjRvr6+tzY2fPnk1tr62tdfu0tra6sVtuucWNrVixYtKPWV9f7/bxxg6EJ/+EeP822v5JRGYsJbtIJJTsIpFQsotEQskuEgklu0gkJiy9kWwD8EsACwGMAdhiZq+QfAHAHwE4mdz1eTP7sFQDnQoKmQjjbT8EhCd+nDp1yo15WzwBQH9/f2p7Y2Oj22fp0qVu7Pbbb3djbW1tbsxbuy60tl6hE4pCvNJbqFw6U+VTZx8B8GMz20lyDoDPSH6UxH5uZn9TuuGJSFby2eutD0BfcnuQ5D4Ai0o9MBHJ1qTey5BcBmANgE+TpmdIdpF8g6T/PlFEKi7vZCfZAOAdAM+a2TkArwJYAWA1cmf+l5x+m0l2kuwcGBgofsQiUpC8kp1kNXKJ/qaZvQsAZtZvZqNmNgbgNQDr0vqa2RYz6zCzjtBFIhEprQmTnbnLoK8D2GdmL49rbxl3t0cA7Ml+eCKSlXyuxq8H8AMAn5PclbQ9D+AJkqsBGIBuAD8qwfimFK9cU+h6ZqEtjfbt2+fGDhw44Ma8mXSLFy92+4TWkluyZIkba2hocGOXLl1KbQ+VIkPlsNCsvdCMuFmz0l/iMc56y+dq/McA0o7MjK6pi8w08X2yQCRSSnaRSCjZRSKhZBeJhJJdJBJacLLEQmU5rzwFAEePHnVjR44ccWNNTU2p7aEZau3t7W5s3rx5bizE+7SkVwqbSKj0FiqjebPsQmW+6bz9U4jO7CKRULKLRELJLhIJJbtIJJTsIpFQsotEQqW3DIRmcoX2cxscHHRjvb29buzw4cNuzCtDVVdXu33mzp3rxurq6tzYyMiIG/PKiqESWigWWqgyVCrzYjEuOBnfbywSKSW7SCSU7CKRULKLRELJLhIJJbtIJFR6mwSvjBMqvQ0PD7uxCxcuuLGTJ08WFJs/f35qe2i2WWjhyFC/0O/mLQIZKnnV1NS4sVB5LbTgpFeKVOlNRGYsJbtIJJTsIpFQsotEQskuEokJr8aTrAWwA8Ds5P7/bGY/I9kEYBuAZcht//SYmWmb1uuUYpuh0KQQbzJJ6Ip7fX29GwutoReqQgwNDbkxT6FX42PcyqkQ+ZzZhwH8vpl9G7ntmTeSvBvAcwC2m9mtALYnP4vIFDVhslvO+eTH6uTLADwMYGvSvhXAplIMUESyke/+7FXJDq4nAHxkZp8CWGBmfQCQfL+5ZKMUkaLllexmNmpmqwEsBrCO5LfyfQKSm0l2kuz01hIXkdKb1NV4MzsD4D8BbATQT7IFAJLvJ5w+W8ysw8w6GhsbixutiBRswmQneRPJecntOgB/AOB/AXwA4Mnkbk8CeL9EYxSRDOQzEaYFwFaSVcj95/C2mf0ryU8AvE3yKQBHAHyvhOOcErzJE6G100LlqebmZje2dOnSgh5zyZIlqe2hbZxC69OFnitU8vJKZaHJM6FYaEJOaJ08r0wZGnvod57OJkx2M+sCsCal/RsA95diUCKSPX2CTiQSSnaRSCjZRSKhZBeJhJJdJBIMzWrK/MnIkwCu7l00H8Cpsj25T+O4lsZxrek2jqVmdlNaoKzJfs0Tk51m1lGRJ9c4NI4Ix6G38SKRULKLRKKSyb6lgs89nsZxLY3jWjNmHBX7m11Eyktv40UioWQXiURFkp3kRpJfkjxIsmILVZLsJvk5yV0kO8v4vG+QPEFyz7i2JpIfkTyQfC/5Sh/OOF4geSw5JrtIPlSGcbSR/A+S+0juJfnHSXtZj0lgHGU9JiRrSf43yd3JOP4saS/ueJhZWb8AVAH4CsAtAGoA7AbQXu5xJGPpBjC/As97L4C1APaMa/trAM8lt58D8FcVGscLAH5S5uPRAmBtcnsOgP0A2st9TALjKOsxAUAADcntagCfAri72ONRiTP7OgAHzeyQmV0G8CvkVqqNhpntAHD6uuayr9brjKPszKzPzHYmtwcB7AOwCGU+JoFxlJXlZL6icyWSfRGAo+N+7kEFDmjCAPya5GckN1doDFdNpdV6nyHZlbzNL+vCgSSXIbdYSkVXML5uHECZj0kpVnSuRLKnrQdUqfrfejNbC+BBAE+TvLdC45hKXgWwArkNQfoAvFSuJybZAOAdAM+a2blyPW8e4yj7MbEiVnT2VCLZewC0jft5MYDeCowDZtabfD8B4D3k/sSolLxW6y01M+tPXmhjAF5DmY4JyWrkEuxNM3s3aS77MUkbR6WOSfLcZzDJFZ09lUj23wC4leRykjUAvo/cSrVlRbKe5JyrtwF8F8CecK+SmhKr9V59MSUeQRmOCXOrP74OYJ+ZvTwuVNZj4o2j3MekZCs6l+sK43VXGx9C7krnVwD+pEJjuAW5SsBuAHvLOQ4AbyH3dvAKcu90ngLQjNyeeQeS700VGsc/APgcQFfy4mopwzg2IPenXBeAXcnXQ+U+JoFxlPWYALgTwP8kz7cHwJ8m7UUdD31cViQS+gSdSCSU7CKRULKLRELJLhIJJbtIJJTsIpFQsotE4v8AZj022F13GxAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUS0lEQVR4nO3dfYyV5ZnH8e9V5K0ytAgjIOKAgIaXIpAJqUVbV7cWbRrtH33RTddNXGl2a7JNun+YbrJ196/uZtumm92Y0NWUNq3WbjV1N7TVUg2aWnTAUTBQEQuUgjAoCigwDlz7x3nIDniue2ae8zYz9++TTM6Z+zrPee55OBfPOc917vs2d0dERr8PtLoDItIcSnaRTCjZRTKhZBfJhJJdJBNKdpFMKNlFMqFklyEzsyVm9iszO2xm7/uihpldZGaPmtk7ZrbHzG5vRT/lXEp2KeM94GHgziD+n0AvMB34C+A+M1vcpL5JwPQNutHFzHYD/wH8JdAB/BK4w91PNmBf84Gd7m792i4EjgBL3P2Vou2HwJ/c/Z5690EGT2f20enzwGpgLrAU+KtqDzKza8zsrcTPNSX2fQVw+myiF14EdGZvsQta3QFpiH939/0AZvY/wLJqD3L3Z4AP13nfk4C3z2t7G2ir835kiHRmH51e73f/XSoJ2CzHgcnntU0GjjWxD1KFkj1jZnatmR1P/Fxb4mlfAS4wswX92q4CXq5Pr6UsvY3PmLs/TYmzvpkZMB4YV/w+ofJ0fsrd3zGzR4B/NrO/pvIR4hbgY3XruJSiM7uU0QGc4P/P1ieA3/eL/y0wETgEPAj8jbvrzN5iKr2JZEJndpFMKNlFMqFkF8mEkl0kE00tvU2bNs07OjqauUuRrOzZs4fDhw9btVhNyW5mq4HvAmOA/3L3b6Ye39HRwW9/+9tadjksVcrO1aWqHalY2ecs83wpZfvYrOcb6Dnrva/h7mMfi7/OUPptvJmNoTKU8SZgEXCbmS0q+3wi0li1fGZfCbzq7q+5ey/wEJVvSonIMFRLss8C/tjv931F2znMbI2ZdZlZV09PTw27E5Fa1JLs1T70vO/Dk7uvdfdOd+9sb2+vYXciUotakn0fMLvf75cC+2vrjog0Si3J/jywwMzmmtk44IvAY/XplojUW+nSm7v3mdndwK+olN4e0MgmkeGrpjq7u68H1tepLyLSQPq6rEgmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmaloRxsx2A8eA00Cfu3fWo1MiUn81JXvhz9z9cB2eR0QaSG/jRTJRa7I78LiZbTazNdUeYGZrzKzLzLp6enpq3J2IlFVrsq9y9xXATcBXzOzj5z/A3de6e6e7d7a3t9e4OxEpq6Zkd/f9xe0h4FFgZT06JSL1VzrZzexCM2s7ex+4EdhWr46JSH3VcjV+OvComZ19nh+7+y/r0isRqbvSye7urwFX1bEvItJAKr2JZELJLpIJJbtIJpTsIpmox3fjs/H2228PqR3gAx+I/z+dO3duGDt58mQYGzNmTBg7duxY1fZ9+/aF25w5cyaM9fb2hrHXXnstjB0/frxqe+p4pI7jrl27wtgFF8Qv4yVLllRtv+aaa8Jt5s2bF8ZGMp3ZRTKhZBfJhJJdJBNKdpFMKNlFMtH0q/Hu3uxd1s3u3burtj/77LPhNgcPHgxjixYtKtWP06dPD3l/f/jDH8JtUv8mR48eDWN79+4NY8WYife58MILw21SV+O7u7vDWF9fXxi79dZbq7ZfccUV4TaXX355GBvJdGYXyYSSXSQTSnaRTCjZRTKhZBfJhJJdJBMaCDMEUWlo8+bN4Ta/+MUvwtg777wTxiZMmBDGUqW3aFDI5MmTS+0rNSBn0qRJYezSSy+t2j5jxoxwmylTpoSxqOwJ6ZJdNKPx1KlTw21GK53ZRTKhZBfJhJJdJBNKdpFMKNlFMqFkF8lE00tvqTnIhrsVK1ZUbU+N5EqNrtq0aVMYS402S+0v6uP1118fbjN//vwwlirZpVbl/dCHPlS1ffz48eE269evD2NPPfVUGFu8eHEYu+GGG6q2p+b/G8mv0ZQB/yoze8DMDpnZtn5tF5nZE2a2s7iNC6QiMiwM5r+w7wOrz2u7B9jg7guADcXvIjKMDZjs7r4RePO85luAdcX9dcCt9e2WiNRb2Q8n0939AEBxe3H0QDNbY2ZdZtZ1+PDhkrsTkVo1/EqEu691905375w2bVqjdycigbLJftDMZgIUt4fq1yURaYSypbfHgDuAbxa3Px/shiN5wsmonLR8+fJwm1SJ57bbbgtjZcs/UWkrNUItVQ5LjbCbPXt2GIuWqEpNHPnQQw+FsWg5KYCrr746jF122WVhLDKSX6Mpgym9PQg8C1xpZvvM7E4qSf5JM9sJfLL4XUSGsQHP7O4enX6qf1tBRIal0flVIRF5HyW7SCaU7CKZULKLZEJrvQ1BVA5LTdiYmkSxra0tjE2cODGM9fb2hrHo+I4dOzbcJuXUqVNhLDX6LppMc8eOHeE2W7ZsCWPR2nGQLm9efHH1L3eOGzcu3ObMmTNhbCTTmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTGittyHo6+ur2p4qC6Vi0cgwKF+iTO0vkio1pUbfpfq4a9euqu3PPfdcuM1bb70VxpYtWxbGFi1aFMai0mfq71LpTURGNCW7SCaU7CKZULKLZELJLpKJpl+NL3O1eLiIrsanrt6mBqBccEF8+N97770wljqG0VXmVB9T88ylBuREg10ANm/eXLX9mWeeCbdJHaubbropjH3kIx8JY1H/U8djJL9GU3RmF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQTGggzBFGJKlXGSZXXUoMxUuWw1Pxp0XOmSnkpqcE6Bw4cCGNR6e3VV18Nt0kt1bRy5cow1t7eHsYiURkV0n/zSDaY5Z8eMLNDZratX9u9ZvYnM+sufm5ubDdFpFaDeRv/fWB1lfbvuPuy4md9fbslIvU2YLK7+0bgzSb0RUQaqJYLdHeb2UvF2/xwcnQzW2NmXWbW1dPTU8PuRKQWZZP9PmAesAw4AHwreqC7r3X3TnfvLHMhRUTqo1Syu/tBdz/t7meA7wHxpVIRGRZKld7MbKa7n627fBbYlnr8aFFmCaVUeS1VDkttlxqVFZUBU6W81PJVJ06cCGNPPvlkGIvmmps6dWq4zRe+8IUwtnTp0jCWUmbewNFqwGQ3sweB64BpZrYP+AZwnZktAxzYDXy5cV0UkXoYMNnd/bYqzfc3oC8i0kD6uqxIJpTsIplQsotkQskukgmNehuCaDRU2dFrqZFXqdFyqbJRtL+ypbw33ngjjHV3d4ex3bt3V22/5JJLwm1Syzi1tbWFsdRxjJaoSh2P0Sq/v1gkU0p2kUwo2UUyoWQXyYSSXSQTSnaRTKj0VgepCSej0g+kyz9l14GLpCapPH78eBh7+umnw9jzzz8fxqKRdFdffXW4zeLFi8NYasRhb29vGIvKijmOetOZXSQTSnaRTCjZRTKhZBfJhJJdJBO6Gj8E0VX31BX3VKwRS0OVGayzf//+MLZ+fbz+xyuvvBLGlixZUrX9uuuuC7fp6OgIY6mKR+rKevR3p7ZJ/ZuNZDqzi2RCyS6SCSW7SCaU7CKZULKLZELJLpKJwawIMxv4ATADOAOsdffvmtlFwE+AOVRWhfm8ux9pXFdbLyrJpEo1qRJPKpaaVy1VRotiR47E/zQvvPBCqViq/wsXLhxSO8D48ePDWGoZqtQgmeh4jNbyWspgzux9wNfcfSHwUeArZrYIuAfY4O4LgA3F7yIyTA2Y7O5+wN23FPePAduBWcAtwLriYeuAWxvURxGpgyF9ZjezOcByYBMw/exKrsXtxXXvnYjUzaCT3cwmAT8DvuruR4ew3Roz6zKzrp6enjJ9FJE6GFSym9lYKon+I3d/pGg+aGYzi/hM4FC1bd19rbt3untne3t7PfosIiUMmOxWueR6P7Dd3b/dL/QYcEdx/w7g5/XvnojUy2BGva0CvgRsNbPuou3rwDeBh83sTmAv8LmG9HAESJWgolFokB7J9e6774axKVOmhLFoRNzLL78cbvPTn/40jO3YsSOMXXXVVWHs2muvrdo+b968cJuy88KVGcFWdhTdSDZgsrv7M0D0199Q3+6ISKPoG3QimVCyi2RCyS6SCSW7SCaU7CKZ0ISTQ5AaiVZGqixXdrmjqMSWKq9t2LAhjLW1tYWxT33qU2Fs1apVVdsnT54cbpNahio1OWeZpbLKjlQcyXRmF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQTKr0NQapUVkZqzbbU5IupEuDevXurtqcmjjx27FgYu/HGG8PYJz7xiTB2ySWXVG1PjTZLKbv2XY4TS0Z0ZhfJhJJdJBNKdpFMKNlFMqFkF8lE06/Gj+Sro9GAi9TAidSV4lQsNRDm0KGqE/kCsHXr1qrt0VV6gDlz5oSxz3zmM2EsNQfduHHjqran5taLtoH0MT516lQYi/7NUlf3R/JrNEVndpFMKNlFMqFkF8mEkl0kE0p2kUwo2UUyMWDpzcxmAz8AZgBngLXu/l0zuxe4Czi7NOvX3X19ozo6GqXKP9HcaQC/+93vwthvfvObqu2p+d1uv/32MLZ69eowNmPGjDAWlcNSA2HKzCU3kOgYly2XjmSDqbP3AV9z9y1m1gZsNrMnith33P3fGtc9EamXwaz1dgA4UNw/ZmbbgVmN7piI1NeQPrOb2RxgObCpaLrbzF4yswfMLF5aVERabtDJbmaTgJ8BX3X3o8B9wDxgGZUz/7eC7daYWZeZdfX09FR7iIg0waCS3czGUkn0H7n7IwDuftDdT7v7GeB7wMpq27r7WnfvdPfO9vb2evVbRIZowGS3ymXL+4Ht7v7tfu0z+z3ss8C2+ndPROplMFfjVwFfAraaWXfR9nXgNjNbBjiwG/jyYHY4kpfWico/ZedV++AHPxjGXn/99TD261//Oox1d3dXbb/yyivDbe66664wdtlll4WxEydOhLHIxIkTw1iqvJaKpebri15vqX+zkfwaTRnM1fhngGp/vWrqIiOIvkEnkgklu0gmlOwimVCyi2RCyS6SCS3/NARRSSa1LFRq4sjU5IuPP/54GNu2Lf5Kw6xZ1YctfPrTnw63WbhwYRgrOzosmjwy9Xyp8lrqOKZGD/b29lZtT00qWe9lvoYLndlFMqFkF8mEkl0kE0p2kUwo2UUyoWQXyYRKb0NQpvSWKgu9+eabYWzjxo1hbOfOnWGso6OjavvMmTOrtkO6DFXvElVqtFnqWKViZdfTi6j0JiIjmpJdJBNKdpFMKNlFMqFkF8mEkl0kE00vvaVKOcNdVP5JlZOOHj0axlKj17Zs2RLGjhw5EsaWLl1atT1VekuNNotGr0G6RFVmEs7UWm+pElrqNZXqf6TsBKLDnc7sIplQsotkQskukgklu0gmlOwimRjwaryZTQA2AuOLx/+3u3/DzC4CfgLMobL80+fdPb5MPIqlrhSnrsbv2bMnjKVWvG1rawtj0TJP8+fPD7cpe1U9tV2ZqksjBrtEy02l+nfq1KkwNpIN5sx+Crje3a+isjzzajP7KHAPsMHdFwAbit9FZJgaMNm94njx69jix4FbgHVF+zrg1kZ0UETqY7Drs48pVnA9BDzh7puA6e5+AKC4vbhhvRSRmg0q2d39tLsvAy4FVprZksHuwMzWmFmXmXWlPoeKSGMN6Wq8u78FPAWsBg6a2UyA4vZQsM1ad+9098729vbaeisipQ2Y7GbWbmYfLu5PBP4c2AE8BtxRPOwO4OcN6qOI1MFgBsLMBNaZ2Rgq/zk87O7/a2bPAg+b2Z3AXuBzg9lhavmf4S7q+4QJE8Jtpk+fHsYWLFgQxmbMmBHGUgNGli9fXrV97ty54TYpJ0+eLNWPqCzX19cXbpMqoaWWf0qV0Y4fP161vey8gSPZgMnu7i8B73sFufsbwA2N6JSI1N/o/C9MRN5HyS6SCSW7SCaU7CKZULKLZMKaOSecmfUAZ4d6TQMON23nMfXjXOrHuUZaPzrcveq315qa7Ofs2KzL3TtbsnP1Q/3IsB96Gy+SCSW7SCZamexrW7jv/tSPc6kf5xo1/WjZZ3YRaS69jRfJhJJdJBMtSXYzW21mvzezV82sZRNVmtluM9tqZt1m1tXE/T5gZofMbFu/tovM7Akz21ncTmlRP+41sz8Vx6TbzG5uQj9mm9mTZrbdzF42s78r2pt6TBL9aOoxMbMJZvacmb1Y9OOfivbajoe7N/UHGAPsAi4HxgEvAoua3Y+iL7uBaS3Y78eBFcC2fm3/CtxT3L8H+JcW9eNe4O+bfDxmAiuK+23AK8CiZh+TRD+aekwAAyYV98cCm4CP1no8WnFmXwm86u6vuXsv8BCVmWqz4e4bgTfPa276bL1BP5rO3Q+4+5bi/jFgOzCLJh+TRD+ayivqPqNzK5J9FvDHfr/vowUHtODA42a22czWtKgPZw2n2XrvNrOXirf5Df840Z+ZzaEyWUpLZzA+rx/Q5GPSiBmdW5Hs1eZ2alX9b5W7rwBuAr5iZh9vUT+Gk/uAeVQWBDkAfKtZOzazScDPgK+6e7yUTvP70fRj4jXM6BxpRbLvA2b3+/1SYH8L+oG77y9uDwGPUvmI0SqDmq230dz9YPFCOwN8jyYdEzMbSyXBfuTujxTNTT8m1frRqmNS7Psthjijc6QVyf48sMDM5prZOOCLVGaqbSozu9DM2s7eB24EtqW3aqhhMVvv2RdT4bM04ZhYZSbP+4Ht7v7tfqGmHpOoH80+Jg2b0blZVxjPu9p4M5UrnbuAf2hRHy6nUgl4EXi5mf0AHqTydvA9Ku907gSmUlkzb2dxe1GL+vFDYCvwUvHimtmEflxD5aPcS0B38XNzs49Joh9NPSbAUuCFYn/bgH8s2ms6Hvq6rEgm9A06kUwo2UUyoWQXyYSSXSQTSnaRTCjZRTKhZBfJxP8Bg73z2LgfrqwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAASJUlEQVR4nO3dbYxc5XnG8f+F4xfATsD1YlbGeME4YBNgQRs3gQSo40QOX+wgJcSKEqqgOh9AKhKRilKpoVKlplWTKJVQJPOiOBVNggoopHEKDmrADgh7nRLb1ARwbN78skuDX8AtxvbdD3OsrM08Z8czZ2bW+1w/aTUzzz1nz50TX5zZ88w5RxGBmY1/p3W7ATPrDIfdLBMOu1kmHHazTDjsZplw2M0y4bCbZcJht5Mm6SOSHpP0pqT3fVFD0q8k/Z+kt4uf33WjTzuew27NeA94ELil5D23RcTU4ufiDvVlJRz2cUbSDklfl7RJ0j5JP5E0pcp1RMTvIuI+4Pkqf6+1l8M+Pn0BWAJcAFwO/Hm9N0n6hKS9JT+faKGHvy8+5v9a0vUt/B6ryAe63YC1xT9HxE4AST8D+uu9KSLWAWe1Yf1/Bfw3cAj4IvAzSf0Rsa0N67IGec8+Pu0e8fwgMLWTK4+IZyPiQES8GxGrgF8DN3SyB3s/hz1jkj454oh5vZ9PVrSqAFTR77Im+WN8xiJiLU3s9SUJmAxMKl5Pqf26eFfSWcCfAk8Ch4GbgGuB26vp2prlsFsz5gDbR7z+X+AVoA+YCPwdcAlwBHgBWBYRnmvvMvniFWZ58N/sZplw2M0y4bCbZcJhN8tER4/Gz5gxI/r6+jq5SrOs7NixgzfffLPudxpaCrukJcD3gAnAvRHxrbL39/X1sWHDhlZWaWYlPvrRjyZrTX+MlzQBuBv4LLAAWC5pQbO/z8zaq5W/2RcCL0fE7yPiEPBjYGk1bZlZ1VoJ+yzgtRGvXy/GjiNphaRBSYPDw8MtrM7MWtFK2OsdBHjf1/EiYmVEDETEQE9PTwurM7NWtBL214HZI16fB+xsrR0za5dWwr4BmCfpAkmTqF2k4NFq2jKzqjU99RYRhyXdBjxGbert/ojwNcnMxqiW5tkjYjWwuqJezKyN/HVZs0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0y0dEcYSTuAA8AR4HBEDFTRlJlVr6WwF/4sIt6s4PeYWRv5Y7xZJloNewCPS9ooaUW9N0haIWlQ0uDw8HCLqzOzZrUa9msi4irgs8Ctkq498Q0RsTIiBiJioKenp8XVmVmzWgp7ROwsHoeAR4CFVTRlZtVrOuySzpQ07dhz4DPAlqoaM7NqtXI0fibwiKRjv+dfI+I/KunKzCrXdNgj4vfAFRX2YmZt5Kk3s0w47GaZcNjNMuGwm2Wiiu/GZ+PAgQN1x/fv359c5rTT0v897e3tbaqPYgakrrfeeqvu+BtvvNHUulL/mwFee+21k17u6NGjyWXKtuP27duTtbLfefnll9cdX7RoUXKZefPmJWunMu/ZzTLhsJtlwmE3y4TDbpYJh90sEz4afxK2bdtWd3zt2rXJZXbv3p2sLViwoOWeTrRz586646neR7Nv375krewIeWoWYurUqcll9u7dm6xt3LgxWYuIZG3p0qV1x+fPn59cxkfjzeyU5rCbZcJhN8uEw26WCYfdLBMOu1kmPPV2EoaGhuqOl029Pf7448naoUOHWu7pRJMnT647XjblNXHixGSt7ESestqcOXNOahzgrLPOStZefvnlZO3gwYPJ2vnnn193fObMmcllxivv2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmPPV2EhYvXlx3fPbs2cllLrvssmRt3bp1ydrhw4cbb2yEhQvr324vdfYXwNy5c5taV9nZZmeccUbd8bLpxieffDJZK9tW5557brJ23XXX1R3/8Ic/nFxmvBp1zy7pfklDkraMGJsuaY2kl4rHs9vbppm1qpGP8T8AlpwwdifwRETMA54oXpvZGDZq2CPiKeAPJwwvBVYVz1cBy6pty8yq1uwBupkRsQugeDwn9UZJKyQNShocHh5ucnVm1qq2H42PiJURMRARAz09Pe1enZklNBv2PZJ6AYrH+meImNmY0ezU26PAzcC3isefVtbRGJY6y+uiiy5KLnPrrbcma1/96leb6qOZKa8zzzwzucykSZOa6qNMaltt2rQpucy9996brJVdjHLZsmXJWmqKrewWWuNVI1NvPwKeAS6W9LqkW6iF/NOSXgI+Xbw2szFs1D17RCxPlD5VcS9m1kb+uqxZJhx2s0w47GaZcNjNMuGz3ipQdsHG6dOnd7CTsSN1EcgXXnghuczg4GCy9t577yVrl156abJ2zjn1v9xZNvVWNrV5KvOe3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCU2+noGamjTo91bRt27a642X3xdu/f3+y1t/fn6xdccUVydoHP/jBZC033rObZcJhN8uEw26WCYfdLBMOu1kmfDT+FNTM0fN2HHEvu0VV6qSWX/7yl8llJk+enKzdeOONydoll1ySrJ1++ul1x8fryS5lvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDUmzVt9+7dydr69evrjr/44ovJZS688MJkbdGiRcnatGnTkrUcp9hSGrn90/2ShiRtGTF2l6Q3JD1X/NzQ3jbNrFWNfIz/AbCkzvh3I6K/+FldbVtmVrVRwx4RTwF/6EAvZtZGrRygu03SpuJj/tmpN0laIWlQ0uDw8HALqzOzVjQb9u8Dc4F+YBfw7dQbI2JlRAxExEBPT0+TqzOzVjUV9ojYExFHIuIocA+wsNq2zKxqTU29SeqNiF3Fy88BW8reb+PT6tXp47Lr1q2rO97b25tc5ktf+lKyVnaLp7Kz5eyPRg27pB8B1wMzJL0OfBO4XlI/EMAO4Gvta9HMqjBq2CNieZ3h+9rQi5m1kb8ua5YJh90sEw67WSYcdrNM+Kw3Kz0zrOyWTBs2bEjWtm/fXne8r68vuUzZbZymTJmSrFljvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU2zgjqe542fRa2T3bfvGLXyRrqfu5QfoeawsXps+GvvLKK5O1007zfqlV3oJmmXDYzTLhsJtlwmE3y4TDbpYJH40fZ1JH3Y8cOZJcZmhoKFl76KGHkrWyWzlddtlldccXL16cXGbOnDnJmrXOe3azTDjsZplw2M0y4bCbZcJhN8uEw26WiUbuCDMb+CFwLnAUWBkR35M0HfgJ0EftrjBfiIi32teqHZM62QXSU2/vvPNOcpmnn346WXv++eeTtbITaBYsWFB3vL+/P7mMtVcje/bDwB0RMR/4GHCrpAXAncATETEPeKJ4bWZj1Khhj4hdEfGb4vkBYCswC1gKrCretgpY1qYezawCJ/U3u6Q+4ErgWWDmsTu5Fo/nVN6dmVWm4bBLmgo8BNweEemLib9/uRWSBiUNDg8PN9OjmVWgobBLmkgt6A9ExMPF8B5JvUW9F6j7BeuIWBkRAxEx0NPTU0XPZtaEUcOu2qHf+4CtEfGdEaVHgZuL5zcDP62+PTOrSiNnvV0DfBnYLOm5YuwbwLeAByXdArwKfL4tHdr7lF1PLjUt9+qrryaXufvuu5O11G2cAObOnZusXXfddXXHL7744uQy1l6jhj0i1gGpid1PVduOmbWLv0FnlgmH3SwTDrtZJhx2s0w47GaZ8AUnT0LZ2WYpZdNkza6r7Hdu2rSp7vg999yTXGbjxo3J2rvvvpusLV++PFm7/vrr645PmDAhuUynt1VuvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU20no5DROM2e2QfrstmeeeSa5zMGDB5O1q6++OllLTa8BnHfeeXXH27ENPb3WGO/ZzTLhsJtlwmE3y4TDbpYJh90sEz4a32btOEljz549ydrmzZvrjr/yyivJZVJHzgG+8pWvJGvz589P1j7wgfr/tHyyS/d4z26WCYfdLBMOu1kmHHazTDjsZplw2M0yMerUm6TZwA+Bc4GjwMqI+J6ku4C/AI7dmvUbEbG6XY2eqpqdFjp69Giytnbt2mRt9er6/xeUXUtu2bJlydqNN96YrJXdqDP1v7vZKTRPr7WukXn2w8AdEfEbSdOAjZLWFLXvRsQ/ta89M6tKI/d62wXsKp4fkLQVmNXuxsysWif1N7ukPuBK4Nli6DZJmyTdL+nsqpszs+o0HHZJU4GHgNsjYj/wfWAu0E9tz//txHIrJA1KGhweHq73FjPrgIbCLmkitaA/EBEPA0TEnog4EhFHgXuAhfWWjYiVETEQEQNlB3TMrL1GDbtqh0/vA7ZGxHdGjPeOeNvngC3Vt2dmVWnkaPw1wJeBzZKeK8a+ASyX1A8EsAP4Whv6G9fKppP27duXrD322GPJWupWTpdeemlymTvuuCNZO/tsH4oZLxo5Gr8OqDc56jl1s1OIv0FnlgmH3SwTDrtZJhx2s0w47GaZ8AUn26zsLK8jR44ka2vWrEnWNmzYkKydf/75dcdvuumm5DJz585N1lIXjoTmblHV6bPXxkofY4H37GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTnnprs7ILR7799tvJ2s9//vNkbdu2bcnaRRddVHe87Oy1ZqfXylR9wclm5TjFluI9u1kmHHazTDjsZplw2M0y4bCbZcJhN8uEp95OQtm0Uco777yTrD399NPJ2vr165O1sim7D33oQ3XHZ86cmVym09NhVTvV++8U79nNMuGwm2XCYTfLhMNulgmH3SwTox6NlzQFeAqYXLz/3yLim5KmAz8B+qjd/ukLEfFW+1rtvmaO7B46dChZ27p1a7K2d+/eZG3q1KnJ2rx58+qOl93+qR1HrDt57TcfcW9MI3v2d4FFEXEFtdszL5H0MeBO4ImImAc8Ubw2szFq1LBHzbGJ3YnFTwBLgVXF+CpgWTsaNLNqNHp/9gnFHVyHgDUR8SwwMyJ2ARSP57StSzNrWUNhj4gjEdEPnAcslPSRRlcgaYWkQUmDw8PDTbZpZq06qaPxEbEX+BWwBNgjqRegeBxKLLMyIgYiYqCnp6e1bs2saaOGXVKPpLOK56cDi4EXgEeBm4u33Qz8tE09mlkFGjkRphdYJWkCtf84PBgR/y7pGeBBSbcArwKfb2Ofp6wzzjgjWbv66quTtRkzZiRrs2bNStY+/vGP1x2/4IILkss0O3VV9Qko7Tihxbd/+qNRwx4Rm4Ar64z/D/CpdjRlZtXzN+jMMuGwm2XCYTfLhMNulgmH3SwT6uQUhKRh4JXi5QzgzY6tPM19HM99HO9U62NORNT99lpHw37ciqXBiBjoysrdh/vIsA9/jDfLhMNuloluhn1lF9c9kvs4nvs43rjpo2t/s5tZZ/ljvFkmHHazTHQl7JKWSPqdpJclde1ClZJ2SNos6TlJgx1c7/2ShiRtGTE2XdIaSS8Vj2d3qY+7JL1RbJPnJN3QgT5mS/pPSVslPS/pL4vxjm6Tkj46uk0kTZG0XtJviz7+thhvbXtEREd/gAnANuBCYBLwW2BBp/soetkBzOjCeq8FrgK2jBj7R+DO4vmdwD90qY+7gK93eHv0AlcVz6cBLwILOr1NSvro6DYBBEwtnk8EngU+1ur26MaefSHwckT8PiIOAT+mdqXabETEU8AfThju+NV6E310XETsiojfFM8PAFuBWXR4m5T00VFRU/kVnbsR9lnAayNev04XNmghgMclbZS0oks9HDOWrtZ7m6RNxcf8tv85MZKkPmoXS+nqFYxP6AM6vE3acUXnboS93nWCujX/d01EXAV8FrhV0rVd6mMs+T4wl9oNQXYB3+7UiiVNBR4Cbo+I/Z1abwN9dHybRAtXdE7pRthfB2aPeH0esLMLfRARO4vHIeARan9idEtDV+ttt4jYU/xDOwrcQ4e2iaSJ1AL2QEQ8XAx3fJvU66Nb26RY915O8orOKd0I+wZgnqQLJE0CvkjtSrUdJelMSdOOPQc+A2wpX6qtxsTVeo/9Yyp8jg5sE9WuCnkfsDUivjOi1NFtkuqj09ukbVd07tQRxhOONt5A7UjnNuCvu9TDhdRmAn4LPN/JPoAfUfs4+B61Tzq3AH9C7Z55LxWP07vUx78Am4FNxT+u3g708Qlqf8ptAp4rfm7o9DYp6aOj2wS4HPivYn1bgL8pxlvaHv66rFkm/A06s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwT/w9vmDMUFyve2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 然后我们看看经典奇异值的分解效果\n",
    "U, sigma, V = np.linalg.svd(imgmat)\n",
    "\n",
    "for i in range(5, 16, 5):\n",
    "    reconstimg = np.matrix(U[:, :i]) * np.diag(sigma[:i]) * np.matrix(V[:i, :])\n",
    "    plt.imshow(reconstimg, cmap='gray')\n",
    "    title = \"n = %s\" % i\n",
    "    plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 然后我们再来看看量子版本的分解效果：\n",
    "time_start = time.time()\n",
    "\n",
    "# 超参数设置\n",
    "N = 5           # 量子比特数量\n",
    "T = 8           # 设置想要学习的阶数\n",
    "D = 80          # 量子神经网络的深度\n",
    "ITR = 200       # 迭代次数\n",
    "LR = 0.02       # 学习速率\n",
    "SEED = 14       # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(2 * T, 0, -2).astype('complex128')\n",
    "\n",
    "\n",
    "def Mat_generator():\n",
    "    imgmat = np.array(list(img.getdata(band=0)), float)\n",
    "    imgmat.shape = (img.size[1], img.size[0])\n",
    "    lenna = np.matrix(imgmat)\n",
    "    return lenna.astype('complex128')\n",
    "\n",
    "M_err = Mat_generator()\n",
    "U, D, V_dagger = np.linalg.svd(Mat_generator(), full_matrices=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置电路参数\n",
    "cir_depth = 80                           # 电路深度\n",
    "block_len = 1                            # 每个模组的长度\n",
    "theta_size = N * block_len * cir_depth   # 网络参数 theta 的大小\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], which_qubit)\n",
    "\n",
    "        for which_qubit in range(1, N):\n",
    "            cir.cnot([which_qubit - 1, which_qubit])\n",
    "\n",
    "    return cir.U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NET(fluid.dygraph.Layer):\n",
    "    \n",
    "    # 初始化长度为 theta_size 的可学习参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "    def __init__(self, shape, param_attr=fluid.initializer.Uniform(low=0.0, high=2 * np.pi), dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # 创建用来学习 U 的参数 theta\n",
    "        self.theta = self.create_parameter(shape=shape, attr=param_attr, dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的参数 phi\n",
    "        self.phi = self.create_parameter(shape=shape, attr=param_attr, dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 我们需要将 Numpy array 转换成 Paddle 动态图模式中支持的 variable\n",
    "        self.M = fluid.dygraph.to_variable(Mat_generator())\n",
    "        self.weight = fluid.dygraph.to_variable(weight)\n",
    "\n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self):\n",
    "        \n",
    "        # 获取量子神经网络的酉矩阵表示\n",
    "        U = U_theta(self.theta)\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        \n",
    "        V = U_theta(self.phi)\n",
    "        V_dagger = dagger(V)\n",
    "        \n",
    "        # 初始化损失函数和奇异值存储器\n",
    "        loss = 0 \n",
    "        singular_values = np.zeros(T)\n",
    "        \n",
    "        # 定义损失函数\n",
    "        for i in range(T):\n",
    "            loss -= self.weight.real[i] * matmul(U_dagger, matmul(self.M, V)).real[i][i]\n",
    "            singular_values[i] = (matmul(U_dagger, matmul(self.M, V)).real[i][i]).numpy()\n",
    "        \n",
    "        # 函数返回两个矩阵 U 和 V_dagger、 学习的奇异值以及损失函数    \n",
    "        return U, V_dagger, loss, singular_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100% |########################################################################|\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "主程序段总共运行了 1971.4076132774353 秒\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYhklEQVR4nO2dW2xdZXbH/ysmF0icm52LyQUnjiEBBAGZiIiCaBmGFI0EPICGh4EHNJmHQSrS9AFRqdA3WhVGPFRIoYTJVJQBlUGgCrWDoiKEqCAGcqMOmASHOHbs3OOQBEi8+nB2JAN7/e3s47OPh+//kyI73zrf3ut8ey+fc77/WWuZu0MI8eNnUr0dEEKUg4JdiERQsAuRCAp2IRJBwS5EIijYhUiEi6qZbGbrADwDoAHAv7r7k+zxs2fP9paWlnxHLopdieRBJhueO3cutE2aFP+NYzYzyx1nvjPOnj0b2r7++uvQVsTHyZMnh3OGh4dD27fffnvB5wJiH7/55ptCx5s6dWpomzZtWmiLzsfuHXZd2NqzY7J5RY4X2fr6+nDs2LHchSwc7GbWAOBfANwOoBfAFjN7w93/L5rT0tKCF154Idc2b9688FzR4rMb8ejRo6GtsbExtLEbp6GhIXec+c4u8sDAQGj74osvQtvFF18c2qZMmZI7Hv2RBYCvvvoqtPX19V3wuQDgkksuyR3fu3dvoeOtWLEitK1cuTK09fT05I6zF4PDhw+HNvZHh92P0XoA8R859ocxOtcDDzwQzqnmbfwaAJ+7+x53/wbAHwDcVcXxhBA1pJpgXwRg34j/92ZjQogJSDXBnvfe4wcfJMxsvZl1mlnnsWPHqjidEKIaqgn2XgBLRvx/MYAffMBz9w3u3uHuHbNnz67idEKIaqgm2LcAaDezZWY2BcDPAbwxPm4JIcabwrvx7n7WzB4G8N+oSG8b3f0TNqehoQHRq3tvby+dlweTk5qbm0Mb24kdGhoKbXPmzMkdHxwcDOdEvgN8Z/fSSy8NbYcOHQptkdLAVIFZs2aFthkzZoS2AwcOhLZIDbn88svDOUVlrW3btoW26LkdP348nMPuHabysHVkH2Ej6ZNJutHuPlvDqnR2d38TwJvVHEMIUQ76Bp0QiaBgFyIRFOxCJIKCXYhEULALkQhV7cZfKO5Ov9wfESVIMOmNyRYsYYFlm0W+s+QZllF26tSp0MZgkl0kvRR5XgB/biwhp4hcymQtBpMHIz+YRMWuC/P/9OnToY3JvZGsWDRzM5xzwTOEEH+WKNiFSAQFuxCJoGAXIhEU7EIkQqm78WZGd5Ijot1ilmRStM4c232OVAG20818nD59emhju/gsqeLEiRMXfDyW+DFz5szQxspZRTvrTCUpojIAfBc8Shhh68FUhmh9gWIlzYC4LBWryVcEvbILkQgKdiESQcEuRCIo2IVIBAW7EImgYBciEUpPhImSUM6cORPOi2p7sbpkRWUL5keRRBjWSogloLAkCCYNRZ1f2PEOHjwY2piPbP0XLFiQO97U1BTOYR1ymLxZpDUUk/JYohST0Iq284rWkSXdKBFGCBGiYBciERTsQiSCgl2IRFCwC5EICnYhEqEq6c3MegAMATgH4Ky7d7DHu3soT7Bm9UVq0DFZiMlQjChTiskxDCbHMGmFyYpRWyMmC7EMsKLPLVpjljXGWm+x58wahkbPjR2PXRd2zzHZls2LKNoOK2I8dPa/dPe4+ZgQYkKgt/FCJEK1we4A/mRmH5rZ+vFwSAhRG6p9G3+Tu/eZ2XwAb5nZLnd/Z+QDsj8C64H4q5xCiNpT1Su7u/dlPwcBvAZgTc5jNrh7h7t3sI0UIURtKRzsZjbdzBrP/w7gpwB2jpdjQojxpZq38QsAvJZJGBcB+Hd3/6/RJkUSEGvhE8k4LCOLFSFkRQ9ZEchIPmEyDstQY/IJy0RjEk8kG7GiklFWIQCcPHkytDEZLVp/JgGytZo/f35oY+8YIzmPXTMmRbL7ismUkXzMYH5ExT7ZnMLB7u57AFxbdL4QolwkvQmRCAp2IRJBwS5EIijYhUgEBbsQiVBqwcmLLroICxcuzLX19PSE8yKJisl1TGpi2VVM1jp16lTu+I4dO8I5jKgoI8C/bbhv377QFvWBY+vb2toa2pgUuXXr1tD2wQcf5I4vXbo0nMOkKyZ5sQKRa9euzR2fM2dOOKe9vT209ff3h7a5c+eGNtafL8puKyLl0R6HoUUI8aNCwS5EIijYhUgEBbsQiaBgFyIRSt2NB+JECFaDLkqEYbuV0c458wHgtcKindhdu3aFc44cORLaWAJHpFoAwOHDh0Nb1EKpt7e3kB9sjZkqsH///tzx7u7ucA67LtHxAGDmzJmhLdqdvvnmm8M5UZIJwOsXFk2+ipJyWKJUkUQYvbILkQgKdiESQcEuRCIo2IVIBAW7EImgYBciEUqV3oaHh8O6cUxmKNLqhsEkHkYko/X19YVzdu/eHdpYDT0moTAZLZJ/WA00VvuNrRW7LlFdu8suu6zQubZv3x7aDhw4ENpuuOGG3HFWW4/JayzRhMmU7FpH52PXLJrDrole2YVIBAW7EImgYBciERTsQiSCgl2IRFCwC5EIo0pvZrYRwM8ADLr71dnYXAAvA2gF0APgPnc/Wo0jTAqZNm1a7jiTp4pk0QG8VtiiRYtyx1kGFZOaWHYVy5ZramoKbW1tbbnjq1atCuew2mmszh9r/xRdz8bGxnBOV1dXaGMZZUw6vP3223PHWY1CdjxW7461lIracrFjsuNF68jkv7G8sv8OwLrvjT0KYLO7twPYnP1fCDGBGTXYs37r33+ZuQvApuz3TQDuHl+3hBDjTdHP7AvcvR8Asp9xi00hxISg5ht0ZrbezDrNrPPo0ao+1gshqqBosA+YWQsAZD8Howe6+wZ373D3DlaYXwhRW4oG+xsAHsx+fxDA6+PjjhCiVoxFensJwK0Ams2sF8DjAJ4E8IqZPQTgSwD3juVkZhZKA0wyiGSGohlIkZQH8IKT0bzly5eHc5gEyDKUpk6dGtrY847kQTbn+PHjoY3JUEyWO3jwYO44Kxy5c+fO0Ma46qqrQlvURosV2WTSLLsuTB6kkhi5NhHRfcrkulGD3d3vD0y3jckrIcSEQN+gEyIRFOxCJIKCXYhEULALkQgKdiESofSCk5E8waSQSPJiBQpZ9hqTQYpIZUzKY7IKy4RivcGiYo5s3pkzZ8I5TF5jfrDnHUlAn376aTiH9Y5j2YNXXHFFaIuuNXvOLAOTFYEcGhoKbeyaRT6y+yOKI/V6E0Io2IVIBQW7EImgYBciERTsQiSCgl2IRChVejOzMFuHFV+M5rCsMZa9xuaxgoKRjMYkEmZjfjCJh8lhUYYgW49Tp06FNgaTqCLpk2ZlkefFeqWxrL358/OLKDHZlvnBpC2WEccy6aL7gN0fbB0j9MouRCIo2IVIBAW7EImgYBciERTsQiRCqbvx586dQ1ROur29PZwX7Razml89PT2hbd68eaGNVcCNzsfaJ7HjRXXaAN5aiR0zqvHGdopZLbmBgYHQxna0o53kjz/+OJzDkpfWrl0b2m655ZbQ1t/fnzs+ffr0cA5LyGHtq9g1a21tDW2RUsJUhiiO2K6/XtmFSAQFuxCJoGAXIhEU7EIkgoJdiERQsAuRCGNp/7QRwM8ADLr71dnYEwB+CeC8dvSYu7852rEmTZoUyjyDg2FvyFBOYIkCTJ5iNeiYnBclQTAJjcFquLGafEXq5LH6aEXbFjGpLJIAo/ZUAH/OzI/u7u7QFq0VS+Jh9eJYotTMmTNDG1vj6JhMRotq6LGah2N5Zf8dgHU5479199XZv1EDXQhRX0YNdnd/B8CREnwRQtSQaj6zP2xm281so5mp8boQE5yiwf4sgDYAqwH0A3gqeqCZrTezTjPrjL7iJ4SoPYWC3d0H3P2cuw8DeA7AGvLYDe7e4e4dbNNMCFFbCgW7mY3scH8PgJ3j444QolaMRXp7CcCtAJrNrBfA4wBuNbPVABxAD4BfVesIy8qKsolYHS7W7qhoS6NICmGyCpPJmJzEMspYvb5Dhw7ljrO1YhlgrD7dgQMHQluUbcbkpGXLloW2tra20MYy0aLadTQ7jMhX7Lqwa81qAEb3IztXVNuQydGjBru7358z/Pxo84QQEwt9g06IRFCwC5EICnYhEkHBLkQiKNiFSIRSC042NDSEWW+suF6UscUkI/YFniKSBhBLZawoI/Ojubk5tLGimJGsBcRSGcvyYj6yTDSWSRcV/Fy+fHk4Z+XKlaGNFcVkxUWj9k9sPRYvXhzamBTJsuWYXBrBpNlaZb0JIX4EKNiFSAQFuxCJoGAXIhEU7EIkgoJdiEQoVXpz97C4HpMMooy4olljLAOMZQ1Fx4xkEIA/LybHMD+iwpdAXPSQ+bh3797QtmfPntDW19cX2iLJjvXFY9IVy4osYmOZj6yoJJvHYNesCOwejtAruxCJoGAXIhEU7EIkgoJdiERQsAuRCBNmN57tqLJ6YREsYYEl0EyZMiW0RUkyrG4dOx5TDFhyDVuraPef1UfbtWtXaHvvvfdC25Ejce+QKAGF7bgvXbo0tLE1Zjvd0X3AdtXZ/cZ2wZk6xGrQFdmpL9ISTa/sQiSCgl2IRFCwC5EICnYhEkHBLkQiKNiFSISxtH9aAuD3ABYCGAawwd2fMbO5AF4G0IpKC6j73J22aR0eHg7lKyZbRFIZqxdXi/Y+USuhKPmEzQHitlYAl95Ym6RI8mLS2+effx7atmzZEtoWLlwY2lasWJE7HvkH8GvGWmyxY0ZSFJMvWeswJqWya83ktWgeu08jCZBJg2N5ZT8L4DfuvgrAjQB+bWZXAngUwGZ3bwewOfu/EGKCMmqwu3u/u3+U/T4EoAvAIgB3AdiUPWwTgLtr5KMQYhy4oM/sZtYK4DoA7wNY4O79QOUPAoD4vZQQou6MOdjNbAaAVwE84u7xh80fzltvZp1m1nn0KP1IL4SoIWMKdjObjEqgv+juf8yGB8ysJbO3ABjMm+vuG9y9w907WDMCIURtGTXYrbK99zyALnd/eoTpDQAPZr8/COD18XdPCDFejCXr7SYAvwCww8y2ZmOPAXgSwCtm9hCALwHcO9qB3D2UNVhWUFTPjNUKK1pnrkhdO9aaiLW1KppBxWS0SK5h8trbb78d2j755JPQxiTAa665JnecSVdM8mL3B6trF0l2LIvuwIEDoa2pqSm0sSzAIlmd7P6O7g92b48a7O7+LoDorrxttPlCiImBvkEnRCIo2IVIBAW7EImgYBciERTsQiRCqQUnJ02aFMpGLNuMFY+MYJloTPJitkj+YZILk5OYvMb8379/f2jbvXt37nh3d3c4h7WGWr16dWibN29eaBsaGsodX758eTiHPWfWKou1oWKZkRGsICmDyWvsvoqKXxaRiKvNehNC/AhQsAuRCAp2IRJBwS5EIijYhUgEBbsQiVCq9DY8PBwW12NZSJGNSXK05xUp5Mfkk0jGYf3L2LlYVhOzseKFURHLzz77LJxz+PDh0MZqEFx77bWhLcp6YwVMItkQ4HIYk/MOHTqUO84yB/v7+0MbkymjcwFcpoykN3Z/q9ebECJEwS5EIijYhUgEBbsQiaBgFyIRSk+EiXYzi+w8RruYQFy3DuA7wixxJdpZZzXL2O4o29llLY127NgR2rq6unLHd+3aFc5ZtWpVaLvhhhtC25IlS0JblJzCnldjY2NoY62Venp6QluUvMQUDbbzz1pUsevJ5kWwRKnoeSkRRgihYBciFRTsQiSCgl2IRFCwC5EICnYhEmFU6c3MlgD4PYCFAIYBbHD3Z8zsCQC/BHAwe+hj7v4mO5a7h7XmWCJMJLuwVkJMWmFJJiwRJpLemKzC2haxunsMJvVt27Ytd3zfvn3hnDvuuCO03XZb3PSHPbfIRyZ7suvJZFYmb0b3DktQYjA5jMmDTM6L7rki93BV7Z8AnAXwG3f/yMwaAXxoZm9ltt+6+z+P4RhCiDozll5v/QD6s9+HzKwLwKJaOyaEGF8u6L2MmbUCuA7A+9nQw2a23cw2mpmarwsxgRlzsJvZDACvAnjE3U8AeBZAG4DVqLzyPxXMW29mnWbWeezYsaodFkIUY0zBbmaTUQn0F939jwDg7gPufs7dhwE8B2BN3lx33+DuHe7ewb6vLoSoLaMGu1W+Wf88gC53f3rEeMuIh90DYOf4uyeEGC/Gsht/E4BfANhhZluzsccA3G9mqwE4gB4AvxrtQGYWSihMdjl58uQFzymSvQbwzKVICmG1x5iPLENpcHAwtO3Zsye0RfXTWEbZ9ddfH9ra29tD28DAQGg7fvx47jiT0JYuXRra2DpG52IwiZXVp2P3FZPXmP+RjUm6UUwwxrIb/y6AvLuSaupCiImFvkEnRCIo2IVIBAW7EImgYBciERTsQiRCqQUnWdYby9aJJKqiLZKY5MUKX0aSXdFCg0wy6u7uDm1fffVVaIuksmXLloVz2JedDh48GNpOnz4d2ubOnZs7XjSzjV0XlonGzjfeMP+j1mFAfB+ze7gIemUXIhEU7EIkgoJdiERQsAuRCAp2IRJBwS5EIpQqvZlZKHuxYn2RfEWL6xEZhEl2RXpyseOxopJMjmHrwSSvaE2YPMWeM5N/2PpHxShZ9h2TMFmWF5PlIumNPS927zDYGrP7gK1jRJGCmXplFyIRFOxCJIKCXYhEULALkQgKdiESQcEuRCKUKr2dPXsWUe34tra2cN6JEydyx5nk1dfXF9qamppCGys2GEleCxcuDOdEvgPA3r17Q9u7774b2lj9/VmzZuWOr1ixIpzDZL4tW7aEtpaWltC2ZMmS3PGenp5wDivYyFi1alVoi4pisnOx9S1a5HTGjBmhLZIBi2RMMklRr+xCJIKCXYhEULALkQgKdiESQcEuRCKMuhtvZtMAvANgavb4/3D3x81sLoCXAbSi0v7pPnc/WtQRtgNapEbXtGnTQhtLPGB+FEkyYTu0RVtUsZ3kSNW4+uqrwzlsHdl6MP+jY7IdfAZrh8WSZKLdc5aExK4nI0r+AXiyTnR/s+NF7atYfcWxvLJ/DeCv3P1aVNozrzOzGwE8CmCzu7cD2Jz9XwgxQRk12L3C+T+dk7N/DuAuAJuy8U0A7q6Fg0KI8WGs/dkbsg6ugwDecvf3ASxw934AyH7Or5mXQoiqGVOwu/s5d18NYDGANWYWfwD8Hma23sw6zayzSGtdIcT4cEG78e5+DMDbANYBGDCzFgDIfubuoLj7BnfvcPeO6KucQojaM2qwm9k8M5ud/X4xgJ8A2AXgDQAPZg97EMDrNfJRCDEOjCURpgXAJjNrQOWPwyvu/p9m9r8AXjGzhwB8CeDesZywSEubSAphdb1YcgdLdmGSVyRrsLpqLJmBJeQwG5NkWltbc8ejdkwAcOjQodDGko1YUkiU+MESlFhtPXZdmpubQ9vRo/lqMLsuQ0NDhfwYbzmP3d/RdaF1GUc7obtvB3BdzvhhALeNNl8IMTHQN+iESAQFuxCJoGAXIhEU7EIkgoJdiESwIq1nCp/M7CCA84XXmgHEmk95yI/vIj++y5+bH5e5+7w8Q6nB/p0Tm3W6e0ddTi4/5EeCfuhtvBCJoGAXIhHqGewb6njukciP7yI/vsuPxo+6fWYXQpSL3sYLkQh1CXYzW2dmn5rZ52ZWt9p1ZtZjZjvMbKuZdZZ43o1mNmhmO0eMzTWzt8ysO/s5p05+PGFm+7M12Wpmd5bgxxIz+x8z6zKzT8zsb7LxUteE+FHqmpjZNDP7wMy2ZX78QzZe3Xq4e6n/ADQA2A1gOYApALYBuLJsPzJfegA01+G8twC4HsDOEWP/BODR7PdHAfxjnfx4AsDflrweLQCuz35vBPAZgCvLXhPiR6lrAsAAzMh+nwzgfQA3Vrse9XhlXwPgc3ff4+7fAPgDKsUrk8Hd3wFw5HvDpRfwDPwoHXfvd/ePst+HAHQBWISS14T4USpeYdyLvNYj2BcB2Dfi/72ow4JmOIA/mdmHZra+Tj6cZyIV8HzYzLZnb/Nr/nFiJGbWikr9hLoWNf2eH0DJa1KLIq/1CPa8ci/1kgRucvfrAfw1gF+b2S118mMi8SyANlR6BPQDeKqsE5vZDACvAnjE3eNe1+X7UfqaeBVFXiPqEey9AEY2714MIK5VVEPcvS/7OQjgNVQ+YtSLMRXwrDXuPpDdaMMAnkNJa2Jmk1EJsBfd/Y/ZcOlrkudHvdYkO/cxXGCR14h6BPsWAO1mtszMpgD4OSrFK0vFzKabWeP53wH8FMBOPqumTIgCnudvpox7UMKaWKW43/MAutz96RGmUtck8qPsNalZkdeydhi/t9t4Jyo7nbsB/F2dfFiOihKwDcAnZfoB4CVU3g5+i8o7nYcANKHSRqs7+zm3Tn78G4AdALZnN1dLCX78BSof5bYD2Jr9u7PsNSF+lLomAK4B8HF2vp0A/j4br2o99A06IRJB36ATIhEU7EIkgoJdiERQsAuRCAp2IRJBwS5EIijYhUgEBbsQifD/a7hBdbxhRo0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 记录优化中间过程\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "\n",
    "\n",
    "# 启动 Paddle 动态图模式\n",
    "with fluid.dygraph.guard():\n",
    "    \n",
    "    net = NET([theta_size])\n",
    "    # 一般来说，我们利用Adam优化器来获得相对好的收敛，当然你可以改成SGD或者是RMS prop.\n",
    "    opt = fluid.optimizer.AdagradOptimizer(learning_rate=LR, parameter_list=net.parameters())\n",
    "    \n",
    "    # 优化循环\n",
    "    pbar = ProgressBar()\n",
    "    for itr in pbar(range(ITR)):\n",
    "        \n",
    "        #  前向传播计算损失函数\n",
    "        U, V_dagger, loss, singular_values = net()\n",
    "        \n",
    "        # 在动态图机制下，反向传播极小化损失函数\n",
    "        loss.backward()\n",
    "        opt.minimize(loss)\n",
    "        net.clear_gradients()\n",
    "\n",
    "        # 记录优化中间结果\n",
    "        loss_list.append(loss[0][0].numpy())\n",
    "        singular_value_list.append(singular_values)\n",
    "\n",
    "    # 记录最后学出的两个酉矩阵    \n",
    "    U_learned = U.real.numpy() + 1j * U.imag.numpy()\n",
    "    V_dagger_learned = V_dagger.real.numpy() + 1j*V_dagger.imag.numpy()\n",
    "\n",
    "singular_value = singular_value_list[-1]\n",
    "mat = np.matrix(U_learned.real[:, :T]) * np.diag(singular_value[:T])* np.matrix(V_dagger_learned.real[:T, :])\n",
    "\n",
    "reconstimg = mat\n",
    "plt.imshow(reconstimg, cmap='gray')\n",
    "\n",
    "time_span = time.time() - time_start        \n",
    "print('主程序段总共运行了', time_span, '秒')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献：\n",
    "\n",
    "[1] [Wang, X., Song, Z. & Wang, Y. Variational Quantum Singular Value Decomposition. arXiv:2006.02336 (2020).](https://arxiv.org/abs/2006.02336)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}