tutorials/machine_learning/VQSVD_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 变分量子奇异值分解\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "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": 154,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.433321Z",
     "start_time": "2021-03-09T03:47:28.730090Z"
    }
   },
   "outputs": [],
   "source": [
    "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_quantum.ansatz import Circuit\n",
    "from paddle_quantum.linalg import dagger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 画出优化过程中的学习曲线\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": 156,
   "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": 157,
   "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": 158,
   "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": 159,
   "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": 160,
   "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",
    "RANK = 8           # 设置想要学习的阶数\n",
    "ITR = 100       # 迭代次数\n",
    "LR = 0.02       # 学习速率\n",
    "SEED = 14       # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(3 * RANK, 0, -3).astype('complex128')\n",
    "print('选取的等差权重为：')\n",
    "print(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们搭建如下的量子神经网络结构："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.725007Z",
     "start_time": "2021-03-09T03:47:34.702861Z"
    }
   },
   "outputs": [],
   "source": [
    "# 设置电路参数\n",
    "num_qubits = 3                           # 量子比特数量\n",
    "cir_depth = 20                           # 电路深度"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义量子神经网络\n",
    "def U_theta(num_qubits: int, depth: int) -> Circuit:\n",
    "\n",
    "    # 用 Circuit 初始化网络\n",
    "    cir = Circuit(num_qubits)\n",
    "    \n",
    "    # 搭建层级结构\n",
    "    for _ in range(depth):\n",
    "        cir.ry()\n",
    "        cir.rz()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来完成算法的主体部分：\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.011232Z",
     "start_time": "2021-03-09T03:47:40.712257Z"
    }
   },
   "outputs": [],
   "source": [
    "class VQSVD():\n",
    "    def __init__(self, matrix: np.ndarray, weights: np.ndarray, num_qubits: int, depth: int, rank: int, lr: float, itr: int, seed: int):\n",
    "        \n",
    "        # 超参数\n",
    "        self.rank = rank\n",
    "        self.lr = lr\n",
    "        self.itr = itr\n",
    "        \n",
    "        paddle.seed(seed)\n",
    "        \n",
    "        # 创建用来学习 U_dagger 的量子神经网络\n",
    "        self.cir_U = U_theta(num_qubits, depth)\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的量子神经网络\n",
    "        self.cir_V = U_theta(num_qubits, depth)\n",
    "        \n",
    "        # 将 Numpy array 转换成 Paddle 支持的 Tensor\n",
    "        self.M = paddle.to_tensor(matrix)\n",
    "        self.weight = paddle.to_tensor(weights)\n",
    "        \n",
    "    # 定义损失函数\n",
    "    def loss_func(self):\n",
    "        # 获取量子神经网络的酉矩阵表示\n",
    "        U = self.cir_U.unitary_matrix()\n",
    "        V = self.cir_V.unitary_matrix()\n",
    "        \n",
    "        # 初始化损失函数和奇异值存储器\n",
    "        loss = paddle.to_tensor(0.0)\n",
    "        singular_values = np.zeros(self.rank)\n",
    "\n",
    "        for i in range(self.rank):\n",
    "            loss -= paddle.real(self.weight)[i] * paddle.real(dagger(U) @ self.M @ V)[i][i]\n",
    "            singular_values[i] = paddle.real(dagger(U) @ self.M @ V)[i][i].numpy()\n",
    "        \n",
    "        return loss, singular_values\n",
    "    \n",
    "    def get_matrix_U(self):\n",
    "        return self.cir_U.unitary_matrix()\n",
    "    \n",
    "    def get_matrix_V(self):\n",
    "        return self.cir_V.unitary_matrix()\n",
    "\n",
    "    # 训练 VQSVD 网络\n",
    "    def train(self):\n",
    "        loss_list, singular_value_list = [], []\n",
    "        optimizer = paddle.optimizer.Adam(learning_rate=self.lr, parameters=self.cir_U.parameters()+self.cir_V.parameters())\n",
    "        for itr in range(self.itr):\n",
    "            loss, singular_values = self.loss_func()\n",
    "            loss.backward()\n",
    "            optimizer.minimize(loss)\n",
    "            optimizer.clear_grad()\n",
    "            loss_list.append(loss.numpy()[0])\n",
    "            singular_value_list.append(singular_values)\n",
    "            if itr% 10 == 0:\n",
    "                print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "\n",
    "        return loss_list, singular_value_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -88.4531\n",
      "iter: 10 loss: -1795.0786\n",
      "iter: 20 loss: -2059.0496\n",
      "iter: 30 loss: -2202.6445\n",
      "iter: 40 loss: -2269.9832\n",
      "iter: 50 loss: -2304.1875\n",
      "iter: 60 loss: -2320.8447\n",
      "iter: 70 loss: -2331.9180\n",
      "iter: 80 loss: -2340.2454\n",
      "iter: 90 loss: -2348.0549\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmKUlEQVR4nO3deZwdVZ338c+316Q7nXQn6ewhCRKWsBigQRBRBpFNNIoiMC7gMqgjOjPgo6A+M86MzAtHfdzFQdGRUUEGXFBQNkXZIQkhZAEJISEr6exLp/ff80dVh5umO/R2+3bf+32/XvfVdU/VvXUqBf3tc07VKUUEZmZm/VGU6wqYmdnw5zAxM7N+c5iYmVm/OUzMzKzfHCZmZtZvDhMzM+s3h4mZ9Yqk3ZIOznU9bGhxmNiwI2mVpDNytO/XS/qjpF2Sdkj6raQ5g7j/fccu6VJJD2Z5f/dL+khmWUSMioiV2dyvDT8OE7MeknQycDfwG2AKMAt4CnhooP9SVyKr/39KKsnm91thcZhY3pBULukbktanr29IKk/XjZf0O0nbJW2V9EDHL2tJn5W0Lm1tPCvpzd3s4j+BGyPimxGxKyK2RsQXgEeBL6bftVzSeRl1KpFUL+m49P1Jkh5O6/GUpNMytr1f0jWSHgIagG4DStIRwPeBk9Nup+0Z/wZflfSipJckfV/SyHTdaZLWpse7EfixpJr036Ve0rZ0eVq6/TXAqcB30n18Jy0PSYeky2Mk3Zh+frWkL2T8u14q6cG0PtskvSDpnB6fUBtWHCaWTz4PnATMBV4LnAh8IV13JbAWqAUmAp8DQtJhwOXACRFRBZwFrOr8xZIqgNcD/9vFfm8B3pIu3wRcnLHuLGBzRCyUNBW4A/gSMBb4NHCbpNqM7d8PXAZUAau7O9CIWA58DHgk7XaqTlddCxya/hscAkwF/jnjo5PSfc9I91ME/Dh9fxCwF/hOuo/PAw8Al6f7uLyLqnwbGEMSfG8CPgB8MGP964BngfEkYXyDJHV3XDZ8OUwsn7wX+LeI2BQR9cC/kvxyBmgBJgMzIqIlIh6IZGK6NqAcmCOpNCJWRcTzXXz3WJL/XzZ0sW4DyS9LgJ8Db0/DB+BvSQIG4H3AnRFxZ0S0R8Q9wHzg3Izv+u+IWBoRrRHR0puDT39JXwb8U9pq2gX8B3BRxmbtwL9ERFNE7I2ILRFxW0Q0pNtfQxIKPdlfcfrdV6cttVXA13j53xxgdUT8ICLagJ+QnIOJvTkuGx4cJpZPprD/X/Or0zKArwArgLslrZR0FUBErAD+kaSbapOkmyVN4ZW2kfwintzFusnA5ozvWw68LQ2Ut5MEDCR//V+QdnFtT7um3tDpO9f05oA7qQUqgAUZ3/+HtLxDfUQ0dryRVCHpv9Iuqp3AX4DqNChezXiglFf+m0/NeL+xYyEiGtLFUb04JhsmHCaWT9aT/MLucFBaRvqX85URcTDJL/grOsZGIuLnEfGG9LMBfLnzF0fEHuAR4IIu9vse4L6M9x1dXfOAZWnAQBIU/xMR1Rmvyoi4NnNXvTjezttuJummOjLj+8dExKgDfOZK4DDgdRExGnhjWq5utu+8vxZe+W++rhfHYHnCYWLDVamkERmvEpJf4l+QVCtpPMlYwU8BJJ0n6ZC0K2gHSfdWu6TDJJ2eDtQ3kvwybu9mn1cBl0j6lKSqdPD6S8DJJF1qHW4GzgQ+zsutEtK6vE3SWZKK03qf1jHg3QcvAdMklQFERDvwA+Drkiakxz1V0lkH+I4qkmPeLmks8C9d7KPLCwHSrqtbgGvSf48ZwBXpcVqBcZjYcHUnyS/BjtcXSQa25wOLgaeBhWkZwGzgXmA3SQvjexHxJ5LxkmtJ/sreCEwAru5qhxHxIMmA+vkk4ySrgWOBN0TEcxnbbUj38XrgFxnla0haK58D6klaKv+Hvv9/+EdgKbBR0ua07LMk3XmPpt1W95K0PLrzDWAkyfE/StItlumbwLvTq7G+1cXnPwnsAVYCD5KE54/6dDQ2rMkPxzIzs/5yy8TMzPrNYWJmZv3mMDEzs37LmzCRdHY6FcaKjnsIzMxscOTFAHx6g9VfSaa0WAs8AVwcEcu6+8z48eNj5syZg1NBM7M8sWDBgs0RUdu5PF9mDT0RWNExLbakm0lvGOvuAzNnzmT+/PmDVD0zs/wgqcs54/Klm2sq+09DsZb9p3QwM7Msypcw6RFJl0maL2l+fX19rqtjZpY38iVM1gHTM95Po4v5gSLi+oioi4i62tpXdPmZmVkf5UuYPAHMljQrnafoIuD2HNfJzKxg5MUAfES0SrocuAsoBn4UEUtzXC0zs4KRF2ECEBF3kkz+Z2ZmgyxfurnMzCyHHCa99Ksn1/LTR7t9NLeZWUFymPTSHYs3OkzMzDpxmPRSTUUp2xtacl0NM7MhxWHSS2Mry9ja0Ew+zGlmZjZQHCa9VFNZRnNrO3tb2nJdFTOzIcNh0ks1FaUAbN3TnOOamJkNHQ6TXqqpKANg2x6Pm5iZdXCY9NLYyjRMGtwyMTPr4DDppeoKh4mZWWcOk17a1zLxmImZ2T4Ok14aM7IUCbb6XhMzs30cJr1UXCTGjCx1y8TMLIPDpA/GVpR5zMTMLIPDpA+qK0odJmZmGRwmfTC2ssz3mZiZZXCY9EGNu7nMzPbjMOmDmsoyT6diZpbBYdIHNRVlNLW2s7fZkz2amYHDpE/2Tfbori4zM8Bh0ic1vgvezGw/DpM+8GSPZmb7c5j0gZ9pYma2P4dJH3Q808TPgjczSzhM+mDfZI9umZiZAQ6TPikpLmL0iFK2e8zEzAxwmPTZ2MoyT0NvZpZymPRRdYWnoTcz6+Aw6SNPQ29m9jKHSR/VVJa5ZWJmlnKY9FFNRSnbPGZiZgY4TPqsprKMvS1tnuzRzAyHSZ+NrfCUKmZmHRwmfVTtMDEz28dh0kf7Jnv043vNzIZemEj6oqR1khalr3Mz1l0taYWkZyWdlVF+dlq2QtJVg1FPP9PEzOxlJbmuQDe+HhFfzSyQNAe4CDgSmALcK+nQdPV3gbcAa4EnJN0eEcuyWcGOZ5p4ShUzs6EbJl2ZB9wcEU3AC5JWACem61ZExEoASTen22Y1TKpHehp6M7MOQ66bK3W5pMWSfiSpJi2bCqzJ2GZtWtZd+StIukzSfEnz6+vr+1XBZLLHEk9Db2ZGjsJE0r2SlnTxmgdcB7wGmAtsAL42UPuNiOsjoi4i6mpra/v9fWMry9wyMTMjR91cEXFGT7aT9APgd+nbdcD0jNXT0jIOUJ5VNZWen8vMDIZgN5ekyRlv3wksSZdvBy6SVC5pFjAbeBx4ApgtaZakMpJB+tsHo641nuzRzAwYmgPw/ylpLhDAKuCjABGxVNItJAPrrcAnIqINQNLlwF1AMfCjiFg6GBUdW1nG8g07B2NXZmZD2pALk4h4/wHWXQNc00X5ncCd2axXV2qrytm8u4mIQNJg797MbMgYct1cw0ntqHJa2sJXdJlZwXOY9MOE0eUAbNrVlOOamJnllsOkH2pHJWFS7zAxswLnMOmHCaNHALBpV2OOa2JmllsOk36orXLLxMwMHCb9Mqq8hIqyYo+ZmFnBc5j0U21VuVsmZlbwHCb9NKGq3GMmZlbwHCb95JaJmZnDpN8mVI3wmImZFTyHST/VVpWzq7GVxpa2XFfFzCxnHCb95MuDzcwcJv3WESYehDezQuYw6acJbpmYmTlM+uvllonDxMwKl8Okn8ZVllMkt0zMrLA5TPqpuEiMG1XOpp0OEzMrXA6TATChqpz63Q4TMytcDpMBUOspVcyswDlMBsAET6liZgXOYTIAJlSNYPPuZtraI9dVMTPLCYfJAKitKqetPdi6pznXVTEzywmHyQDwjYtmVugcJgPAU6qYWaFzmAyACVUjALdMzKxwOUwGgKdUMbNC5zAZACPLiqkqL3HLxMwKlsNkgPjxvWZWyBwmA8RhYmaFzGEyQDylipkVMofJAJk4egQbdzYS4bvgzazwOEwGyIxxFTS2tLury8wKksNkgMwYVwnAC5v35LgmZmaDz2EyQGaOqwBg9ZaGHNfEzGzwOUwGyNTqkZQUiVVb3DIxs8KTkzCRdIGkpZLaJdV1Wne1pBWSnpV0Vkb52WnZCklXZZTPkvRYWv4LSWWDeSwdSoqLmD62wi0TMytIuWqZLAHOB/6SWShpDnARcCRwNvA9ScWSioHvAucAc4CL020Bvgx8PSIOAbYBHx6cQ3ilGeMqPGZiZgUpJ2ESEcsj4tkuVs0Dbo6Ipoh4AVgBnJi+VkTEyohoBm4G5kkScDpwa/r5nwDvyPoBdGPmuEpWb9njy4PNrOAMtTGTqcCajPdr07LuyscB2yOitVN5lyRdJmm+pPn19fUDWnFIBuH3NLexebcfkmVmhSVrYSLpXklLunjNy9Y+X01EXB8RdRFRV1tbO+DfP2N8cnnwag/Cm1mBKcnWF0fEGX342Dpgesb7aWkZ3ZRvAaollaStk8ztB93M9F6TVVsaqJs5NlfVMDMbdEOtm+t24CJJ5ZJmAbOBx4EngNnplVtlJIP0t0cyOPEn4N3p5y8BfpODegMwrWYkxUVilQfhzazA5OrS4HdKWgucDNwh6S6AiFgK3AIsA/4AfCIi2tJWx+XAXcBy4JZ0W4DPAldIWkEyhnLD4B7Ny0qLi5hWM9L3mphZwclaN9eBRMSvgF91s+4a4Jouyu8E7uyifCXJ1V5Dwoxxlb7XxMwKzlDr5hr2Zo6rYJUvDzazAuMwGWAzxlWyq7GVrXt8ebCZFQ6HyQCbNT6Z8HGVu7rMrIA4TAZYx1T0vtfEzAqJw2SATasZSZHcMjGzwuIwGWDlJcVMqR7pe03MrKA4TLJg1vhKd3OZWUHpUZhI+gdJo5W4QdJCSWdmu3LD1YxxFe7mMrOC0tOWyYciYidwJlADvB+4Nmu1GuZmT6hix94W1m3fm+uqmJkNip6GidKf5wL/k05logNsX9COO6gGgIWrt+W4JmZmg6OnYbJA0t0kYXKXpCqgPXvVGt4On1zFyNJiFjhMzKxA9HRurg8Dc4GVEdEgaSzwwazVapgrLS5i7vRqh4mZFYyetkxOBp6NiO2S3gd8AdiRvWoNf8fPqGHZhp00NLe++sZmZsNcT8PkOqBB0muBK4HngRuzVqs8cPyMGtrag6fWOHPNLP/1NExa0wdRzQO+ExHfBaqyV63hb98g/Ivu6jKz/NfTMZNdkq4muST4VElFQGn2qjX8jakoZfaEUcxftTXXVTEzy7qetkwuBJpI7jfZSPKs9a9krVZ54vgZNSx8cTvt7X62iZnltx6FSRogPwPGSDoPaIwIj5m8iuNm1LBjbwsrN+/OdVXMzLKqp9OpvAd4HLgAeA/wmKR3Z7Ni+aBuRjJu4kuEzSzf9XTM5PPACRGxCUBSLXAvcGu2KpYPZo2vpKailPmrtnHhCQflujpmZlnT0zGToo4gSW3pxWcLliSOn1HDAl/RZWZ5rqeB8AdJd0m6VNKlwB3AndmrVv44fsZYVtbvYdOuxlxXxcwsa3o6AP9/gOuBY9LX9RHx2WxWLF+cfvgEAO5Z9lKOa2Jmlj09HTMhIm4DbstiXfLSoRNHcXBtJb9/eiPvfd2MXFfHzCwrDtgykbRL0s4uXrsk7RysSg5nkjj3qMk8snILW/c057o6ZmZZccAwiYiqiBjdxasqIkYPViWHu3OOnkRbe3DPso25roqZWVb4iqxBMGfyaGaMq+COpx0mZpafHCaDQBLnHDWZh1dsZnuDu7rMLP84TAbJuUdPorU9fFWXmeUlh8kgOXrqGKZWj+T3S9zVZWb5x2EySCRx7tGTeOC5enbsbcl1dczMBpTDZBCdd8wUWtqCO5/ekOuqmJkNKIfJIDpm2hhmTxjF/85fk+uqmJkNKIfJIJLEBXXTWPjidlZs8jNOzCx/5CRMJF0gaamkdkl1GeUzJe2VtCh9fT9j3fGSnpa0QtK3JCktHyvpHknPpT9rcnFMPfWOY6dSXCRuW7g211UxMxswuWqZLAHOB/7SxbrnI2Ju+vpYRvl1wN8Bs9PX2Wn5VcB9ETEbuC99P2RNqBrBaYfW8suFa2nz43zNLE/kJEwiYnlEPNvT7SVNBkZHxKMREcCNwDvS1fOAn6TLP8koH7IuqJvGSzub+Mtz9bmuipnZgBiKYyazJD0p6c+STk3LpgKZ/UJr0zKAiRHRcXnURmBid18s6TJJ8yXNr6/P3S/y0w+fSE1FKbcucFeXmeWHrIWJpHslLeniNe8AH9sAHBQRxwJXAD+X1OMJJdNWS7d9RxFxfUTURURdbW1tj49loJWVFDFv7lTuWfqSp1cxs7yQtTCJiDMi4qguXr85wGeaImJLurwAeB44FFgHTMvYdFpaBvBS2g3W0R2W+XjhIes9ddNpbmvn10+ue/WNzcyGuCHVzSWpVlJxunwwyUD7yrQba6ekk9KruD4AdITS7cAl6fIlGeVD2pwpozlm2hhufmINSYPKzGz4ytWlwe+UtBY4GbhD0l3pqjcCiyUtAm4FPhYRW9N1fw/8EFhB0mL5fVp+LfAWSc8BZ6Tvh4ULT5jOMxt38dTaHbmuiplZv6hQ/yquq6uL+fPn57QOuxpbOPGa+5g3dwrXvuuYnNbFzKwnJC2IiLrO5UOqm6vQVI0o5bxjJnP7U+vZ3dSa6+qYmfWZwyTHLjrxIBqa2/jdU+tzXRUzsz5zmOTYcQdVM3vCKG5+wpM/mtnw5TDJMUlcdOJBLFqzneUbdua6OmZmfeIwGQLOP3YqI0qLuPGRVbmuiplZnzhMhoCayjLeeexUfrlwHVv3+I54Mxt+HCZDxAdPmUVTazs3Pf5irqtiZtZrDpMh4tCJVZw6ezw3PrKKlrb2XFfHzKxXHCZDyIdOmcVLO5v8jHgzG3YcJkPImw6t5eDxlfz4oVW5roqZWa84TIaQoiJx6SkzWbRmOwtWb331D5iZDREOkyHmXcdNY1xlGV/+/bOeTdjMhg2HyRBTWV7CFWceyuOrtnLX0o25ro6ZWY84TIagC+umc9jEKv7jzmdoam3LdXXMzF6Vw2QIKiku4vNvPYIXtzZw48Orc10dM7NX5TAZot54aC1/c1gt3/rjc2zZ3ZTr6piZHZDDZAj7/FuPoKG5jc/cupi2dg/Gm9nQ5TAZwg6ZUMUX334k9z2ziX//3bJcV8fMrFslua6AHdj7T5rBi1v28IMHXmDGuAo+eMqsXFfJzOwVHCbDwFXnHMHqLQ38+++WMWZkKecfNy3XVTIz24+7uYaB4iLxjYvmUjdzLFfc8hSf/t+n2ONnxpvZEOIwGSYqykr4+UdexydPP4TbFq7lvG8/yILV23JdLTMzwGEyrJQUF3HlmYdx09+dRGNLG++67mE+ddOTrNu+N9dVM7MC5zAZhk46eBz3XvEmPnn6Idy1dCOnf/V+vn7PX2ls8d3yZpYbDpNhqrK8hCvPPIw/fvo0zjxyEt+87znO/eYDPPz85lxXzcwKkMNkmJtaPZJvX3wsN37oRFrbg7/9wWNc/cvFbqWY2aBymOSJNx5ay93/9EY++qaDuenxNbzzew+zavOeXFfLzAqEwySPjCgt5upzjuDHl57A+u17edu3H+RuT2NvZoPAYZKH/ubwCdzxqTcwq7aSj/50AT98YGWuq2Rmec5hkqem1VRwy0dP5qw5k/jSHcv5t98uo92TRZpZljhM8tiI0mK++97j+NAps/jRQy/wiZ8v9MC8mWWF5+bKc8VF4p/fNocp1SP40h3L2dbwONd/oI7RI0pzXTUzyyNumRSIj5x6MN+8aC7zV23jwv96lE27GnNdJTPLIw6TAjJv7lRuuPQEVm/Zw7uue5jn63fnukpmliccJgXmTYfW8vO/O4mGpjbO/97DPLZyS66rZGZ5ICdhIukrkp6RtFjSryRVZ6y7WtIKSc9KOiuj/Oy0bIWkqzLKZ0l6LC3/haSyQT6cYWfu9Gp+9fenMH5UGe+/4XF+/eS6XFfJzIa5XLVM7gGOiohjgL8CVwNImgNcBBwJnA18T1KxpGLgu8A5wBzg4nRbgC8DX4+IQ4BtwIcH9UiGqYPGVfDLj5/CsQdV84+/WMR/3Lmc1rb2XFfLzIapnIRJRNwdER1Pd3oU6Hh04Dzg5ohoiogXgBXAielrRUSsjIhm4GZgniQBpwO3pp//CfCOQTqMYW9MRSn/8+HX8f6TZnD9X1by/hseZ/PuplxXy8yGoaEwZvIh4Pfp8lRgTca6tWlZd+XjgO0ZwdRR3iVJl0maL2l+fX39AFV/eCsrKeLf33EUX73gtSx8cRtv+/aDPOpxFDPrpayFiaR7JS3p4jUvY5vPA63Az7JVj0wRcX1E1EVEXW1t7WDscth49/HTuO3jr6e8pIiLf/Ao//mHZ2hudbeXmfVM1m5ajIgzDrRe0qXAecCbI6Jjno91wPSMzaalZXRTvgWollSStk4yt7deOmrqGO741Kn822+X8b37n+fBFZv56gWv5dCJVbmumpkNcbm6muts4DPA2yOiIWPV7cBFksolzQJmA48DTwCz0yu3ykgG6W9PQ+hPwLvTz18C/GawjiMfVZaX8OV3H8P333cca7Y2cN63HuTb9z1HiwfnzewAcjVm8h2gCrhH0iJJ3weIiKXALcAy4A/AJyKiLW11XA7cBSwHbkm3BfgscIWkFSRjKDcM7qHkp7OPmsw9V7yJM4+cyNfu+Stv/85DLFqzPdfVMrMhSi/3MBWWurq6mD9/fq6rMSzctXQj//fXS6jf3cTFJx7EZ846jOoK385jVogkLYiIus7lQ+FqLhvizjpyEvdd+SY++PpZ/OKJNZz+tT/zs8dW+74UM9vHYWI9UjWilH9+2xx+e/kbOKR2FJ//1RLO+eYD3P/splxXzcyGAIeJ9cqcKaP5xUdP4vvvO47mtnYu/fETXPhfj/DI8743xayQeczE+qy5tZ2fPbaa6+5/nk27mjjp4LF8/LRDeOPs8SSTE5hZvuluzMRhYv3W2NLGTY+/uC9UXlNbyaWnzOL8Y6dSWe7nr5nlE4dJJw6Tgdfc2s4dT6/nxw+tYvHaHVSUFXPeMZN5T910jp9R49aKWR5wmHTiMMmeiGDhi9u55Yk1/G7xevY0tzFzXAVvnzuVeXOn8JraUbmuopn1kcOkE4fJ4NjT1MqdT2/gN4vW8/Dzm2kPOGrqaN52zBTeesxkptVU5LqKZtYLDpNOHCaDb9PORn67eAO/fWr9vrvpjzuomnOPnsy5R09mSvXI3FbQzF6Vw6QTh0luvbilgd8uXs8dizewbMNOAF47vZq3HDGBNx8xkcMnVXmMxWwIcph04jAZOlZt3sMdT2/g7mUv8VTaYpkyZgSvP2Q8Jx88jpNeM44pY0Y4XMyGAIdJJw6ToWnTzkb++Mwm/vzXeh5duYVtDS0ATKgq57XTq5k7vZo5k0czZ8poJlSVO2DMBpnDpBOHydDX3h48s3EXj7+whcVrd7Bo7XZW1u/Zt35cZRmHT67i8EmjOWLyaOZMHs3siaMoLfbEDmbZ0l2Y+I4yG7KKisScKUkrpMPOxhae2bCL5Rt2smz9TpZv3MlPH11NU/pUyLLiIg6dNIqjpozhyKljOHrqGA6fVMWI0uJcHYZZQXDLxIa91rZ2Vm3Zw9L1ScAsXb+TJet3sD3tIisuErMnjOLoqWM4etoYjpo6hiMmjWZkmQPGrLfczdWJwyS/RQTrtu9lybodPL1uB0+v28mSdTvYuqcZSALmNbWVzJk8miOnjGHOlKSrbGyln9NidiDu5rKCIolpNRVMq6ng7KMmA0nAbNjRyNPrdrBk3Q6Wrt/JIyu38OtF6/d9buLocg6fNJrDJlVx6MQqZk8YxcG1lVSNKM3VoZgNCw4TKxiSmFI9kinVIznryEn7yjfvbmL5hp3paxfPbNzFI89voTnj4V8TqsqZNb6Sg8ZWMGNcBdPHVjA1/a4JVeWUeNDfCpzDxAre+FHlnDq7llNn1+4rS8ZhGlixaTcrN+/m+U17WL1lD/f/tZ76XU37fb5IyXdMGjOCCVUjqK0qY1xlOeNHlTF2VDljK8oYW1lGdUUpNRVlHquxvOQwMetCSXERh0wYxSETXjkpZUNzK+u27WXd9r2s397Ihh17eWlnIy/tbGLttgYWrdnO1j1NtHczHFleUrQvWMaMTH7WVJYyZmQZNR3l6c/qilLGjExe5SVFvq/GhiyHiVkvVZSVMHtiFbMnVnW7TVt7sK2hmW17mtmavrbvbWFbQzPbG1rYvu9nC8/X72bb6qSstbsEIrnsefTIEqpGlDKqvITK8mJGlZdQUZYsjywtoaKsmJFlxcnP0mQ58+eI0v3XjygtdkjZgHCYmGVBcZEYP6qc8aPKe/yZiGBPc9t+QbN9b7K8s7GFnXtb2dnYwu7GVnY3tbK7sZUNOxppaG5jT1MrDc1tNDS3dtsi6o5EEjgZYdPxvry0iBFp+YiSIkaWvbxcngbRvp8lRZSXJJ/Zt5yWl3W8ipPty4qLKC2WQyyPOEzMhghJjCovYVR5CdNq+vYdEUFTazt7m9vY25K+mttobGlLw6aNptakrKG5jcbWNhrT5Y7tG1vaaGxJvmNXYyv1u5r2fWdjaxtNLe00trYxEHcVlBXvHzRlJUnIlJUUU1as9H3Ryz/TECotLqK0ZP/3HdvsW1+c+Z1FlBQr3T5Z3vd9JaKkKFnOLC/Z9z0OvZ5wmJjlEUn7WhJ9zKMeiQha2iIJo5Y2mlvbaWxpT5bb2tP3SXlTa3vGz5fXN7e209wW6c+2fWUtbUkgNre105KW7W5qTdcl5a1tkSyn27e0tR+wi7C/Sor0iqAqLRGlRS8vdwRSaUmyTUlREWVdLGcGW0lRsn1HyHUOvcz9dQ7Al8O107piUVw0+AHoMDGzXpNEWUnSchg9RO7BaW8PWtrTcMkInubWJGg6wqi1PVm/Xyi1JZ9rbev4XLLc2p4EW8fyvu9oi/QzHa/Yt9zY0s7uxlaaO0IuXd+cLre2BU3pcrbyT+IVLayOoCkpLuKGS+qYMa5yQPfpMDGzvFBUJMqLiikvAXo+VJVTbe2xXyB1hFXngOoIuuYuypOwfLmV1jkQW7r4XDbmqnOYmJnlSHGRKC4qzouJSH3brpmZ9ZvDxMzM+s1hYmZm/eYwMTOzfnOYmJlZvzlMzMys3xwmZmbWbw4TMzPrt4J9BrykemB1Lz4yHticpeoMVYV4zFCYx12IxwyFedz9PeYZEVHbubBgw6S3JM2PiLpc12MwFeIxQ2EedyEeMxTmcWfrmN3NZWZm/eYwMTOzfnOY9Nz1ua5ADhTiMUNhHnchHjMU5nFn5Zg9ZmJmZv3mlomZmfWbw8TMzPrNYfIqJJ0t6VlJKyRdlev6ZIuk6ZL+JGmZpKWS/iEtHyvpHknPpT+z+WjxnJBULOlJSb9L38+S9Fh6zn8hqSzXdRxokqol3SrpGUnLJZ2c7+da0j+l/20vkXSTpBH5eK4l/UjSJklLMsq6PLdKfCs9/sWSjuvrfh0mByCpGPgucA4wB7hY0pzc1iprWoErI2IOcBLwifRYrwLui4jZwH3p+3zzD8DyjPdfBr4eEYcA24AP56RW2fVN4A8RcTjwWpLjz9tzLWkq8CmgLiKOAoqBi8jPc/3fwNmdyro7t+cAs9PXZcB1fd2pw+TATgRWRMTKiGgGbgbm5bhOWRERGyJiYbq8i+SXy1SS4/1JutlPgHfkpIJZImka8Fbgh+l7AacDt6ab5OMxjwHeCNwAEBHNEbGdPD/XJI8pHympBKgANpCH5zoi/gJs7VTc3bmdB9wYiUeBakmT+7Jfh8mBTQXWZLxfm5blNUkzgWOBx4CJEbEhXbURmJiremXJN4DPAO3p+3HA9ohoTd/n4zmfBdQDP067934oqZI8PtcRsQ74KvAiSYjsABaQ/+e6Q3fndsB+xzlMbD+SRgG3Af8YETsz10VyHXneXEsu6TxgU0QsyHVdBlkJcBxwXUQcC+yhU5dWHp7rGpK/wmcBU4BKXtkVVBCydW4dJge2Dpie8X5aWpaXJJWSBMnPIuKXafFLHc3e9OemXNUvC04B3i5pFUkX5ukkYwnVaVcI5Oc5XwusjYjH0ve3koRLPp/rM4AXIqI+IlqAX5Kc/3w/1x26O7cD9jvOYXJgTwCz0ys+ykgG7G7PcZ2yIh0ruAFYHhH/L2PV7cAl6fIlwG8Gu27ZEhFXR8S0iJhJcm7/GBHvBf4EvDvdLK+OGSAiNgJrJB2WFr0ZWEYen2uS7q2TJFWk/613HHNen+sM3Z3b24EPpFd1nQTsyOgO6xXfAf8qJJ1L0q9eDPwoIq7JbY2yQ9IbgAeAp3l5/OBzJOMmtwAHkUzZ/56I6Dy4N+xJOg34dEScJ+lgkpbKWOBJ4H0R0ZTD6g04SXNJLjooA1YCHyT54zJvz7WkfwUuJLly8UngIyTjA3l1riXdBJxGMtX8S8C/AL+mi3ObBut3SLr8GoAPRsT8Pu3XYWJmZv3lbi4zM+s3h4mZmfWbw8TMzPrNYWJmZv3mMDEzs35zmJj1k6SH058zJf3tAH/357ral9lQ40uDzQZI5r0qvfhMScbcUF2t3x0RowagemZZ5ZaJWT9J2p0uXgucKmlR+uyMYklfkfRE+qyIj6bbnybpAUm3k9yFjaRfS1qQPm/jsrTsWpJZbhdJ+lnmvtI7lr+SPpvjaUkXZnz3/RnPKvlZemOaWVaVvPomZtZDV5HRMklDYUdEnCCpHHhI0t3ptscBR0XEC+n7D6V3JI8EnpB0W0RcJenyiJjbxb7OB+aSPItkfPqZv6TrjgWOBNYDD5HMQfXgQB+sWSa3TMyy50ySeY8WkUxLM47kIUQAj2cECcCnJD0FPEoy8d5sDuwNwE0R0RYRLwF/Bk7I+O61EdEOLAJmDsCxmB2QWyZm2SPgkxFx136FydjKnk7vzwBOjogGSfcDI/qx38y5pdrw/+c2CNwyMRs4u4CqjPd3AR9Pp/ZH0qHpQ6g6GwNsS4PkcJLHJndo6fh8Jw8AF6bjMrUkT058fECOwqwP/BeL2cBZDLSl3VX/TfJslJnAwnQQvJ6uHwv7B+BjkpYDz5J0dXW4HlgsaWE6PX6HXwEnA0+RPOjoMxGxMQ0js0HnS4PNzKzf3M1lZmb95jAxM7N+c5iYmVm/OUzMzKzfHCZmZtZvDhMzM+s3h4mZmfXb/wc9qxaJMSfscQAAAABJRU5ErkJggg==",
      "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",
    "# 创建 VQSVD 网络并进行训练\n",
    "net = VQSVD(matrix=M, weights=weight, num_qubits=num_qubits, depth=cir_depth, rank=RANK, lr=LR, itr=ITR, seed=SEED)\n",
    "loss_list, singular_value_list = net.train()\n",
    "        \n",
    "# 绘制学习曲线\n",
    "loss_plot(loss_list)\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = net.get_matrix_U().numpy()\n",
    "V_dagger_learned = dagger(net.get_matrix_V()).numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来探究下量子版本的奇异值分解的精度问题。在上述部分，我们提到过可以用分解得到的更少的信息来表达原矩阵。具体来说，就是用前 $T$ 个奇异值和前 $T$ 列左右奇异向量重构一个矩阵：\n",
    "\n",
    "$$\n",
    "M_{re}^{(T)} = UDV^{\\dagger}, \\tag{4}\n",
    "$$\n",
    "\n",
    "并且对于一个本身秩 (rank) 为 $r$ 的矩阵 $M$, 误差随着使用奇异值的数量变多会越来越小。经典的奇异值算法可以保证：\n",
    "\n",
    "$$\n",
    "\\lim_{T\\rightarrow r} ||M - M_{re}^{(T)}||^2_2 = 0, \\tag{5}\n",
    "$$\n",
    "\n",
    "其中矩阵间的距离测量由 Frobenius-norm 来计算,\n",
    "\n",
    "$$\n",
    "||M||_2 = \\sqrt{\\sum_{i,j} |M_{ij}|^2}. \\tag{6}\n",
    "$$\n",
    "\n",
    "目前量子版本的奇异值分解还需要很长时间的优化，理论上只能保证上述误差不断减小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.407257Z",
     "start_time": "2021-03-09T03:49:24.027591Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABCkUlEQVR4nO3dd3hUddbA8e+ZSSEJJZQQkN5bgAQiiEqVZkMXKbLKirrr2nd1YdV1X8V1XRuWtaxldVVEBZRmR6QKCAqEXgQkIKGXICWQdt4/7k0IIWUCk0zK+TzPPJm5986dM5Nkzv11UVWMMcZUXJ5AB2CMMSawLBEYY0wFZ4nAGGMqOEsExhhTwVkiMMaYCi4o0AGci1q1amnjxo0DHYYxxpQpy5cvP6CqUbm3l8lE0LhxY5YtWxboMIwxpkwRke15bbeqIWOMqeAsERhjTAVnicAYYyq4MtlGYExpsG/fPkaPHs3GjRvJzMwMdDjGAODxeGjdujXjxo2jdu3aPj3HEoEx52j06NH07t2bt99+m+Dg4ECHYwwAaWlpvP/++4wePZrx48f79JwKUzU0PSGJQU9+wtJHunL1k1OYnpAU6JBMGbdx40ZuvPFGSwKmVAkODmbkyJFs3LjR5+dUiBLB9IQkHpq6hr/pR1zo3cTQ4x/y0NQIAK6Nqxfg6ExZlZmZaUnAlErBwcFFqq6sECWCZ2duonLaAYZ65+MRZah3AZXTDvLszE2BDs0YYwKuQiSCXckp3Bs0DS8ZAHjI4J6gqexKTglwZMacH6/XS2xsLDExMVx99dUkJycHLJZ58+axePFiv51v+vTprF+/PvvxI488wrfffuu38+f26aef8tRTT/l07IkTJ6hZsya//vrrGduvvfZaJk2aBDjxd+jQgdatWxMTE8Mnn3ySfdySJUvo2rUrsbGxtGnThrFjx5KYmEj9+vXPupKPjY1l6dKljB07lnr16hEbG0uLFi0YPHjwGZ/P+agQiaB9tRSGeucTLM4HHCoZDPUuIKbayQBHZiqS6QlJXPLUHJo8+AWXPDXHL+1UYWFhrFy5krVr11KjRg1effVVP0R6bgpKBOnp6UU+X+5E8I9//IO+ffuec3yFGTRoEA8++KBPx4aHhzNgwACmTZuWve3IkSMsXLiQq6++mlWrVjF69GhmzJjBxo0b+eyzz3jggQdYvnw5ADfddBNvvvlm9u9u2LBhNG7cmIYNG/Ldd99ln3Pjxo0cPXqUrl27AnDfffexcuVKNm/ezPDhw+nTpw/79+8/7/deIRLBi3VnIZy5EpuHTP5d95sARWQqmqx2qqTkFBRISk7hoalr/NppoVu3biQlOefbunUrAwcOpHPnznTv3j274XDv3r385je/oWPHjnTs2DH7i/v5558nJiaGmJgYXnzxRQASExNp06YNf/jDH2jXrh39+/cnJcUpRb/00ku0bduWDh06cP3115OYmMjrr7/OCy+8QGxsLN999x2jRo3i9ttvp2vXrvz1r39l7NixjBs3LjvemJgYEhMTARg/fjwdOnSgY8eOjBw5ksWLF/Ppp58yZswYYmNj2bp1K6NGjcq+qp49ezZxcXG0b9+eW265hVOnTgHO9DOPPvoonTp1on379nk2mF500UWsW7cu+3GvXr1YtmwZ7777LnfffTcAn332GV27diUuLo6+ffuyd+/es84zYsQIJk6cmP142rRpDBgwgPDwcMaNG8ff/vY3mjRpAkCTJk3429/+xnPPPQc4XY/r1q0LOKW6tm3b5nnOiRMncv311+f5+x4+fDj9+/fnww8/zHN/UVSIxuKmJ9eBnHlFEirpznZj/GT4G9/nuy9hRzKpGWcW+VPSMhj72TqujavHoeOp3DFh+Rn7J/2xm8+vnZGRwezZs7n11lsBuO2223j99ddp0aIFS5cu5c4772TOnDnce++99OzZk2nTppGRkcGxY8dYvnw577zzDkuXLkVV6dq1Kz179qR69eps3ryZjz76iP/+978MGzaMKVOmcOONN/LUU0+xbds2QkNDSU5OJjIykttvv53KlSszevRoAN5++2127tzJ4sWL8Xq9jB07Ns/Y161bxz//+U8WL15MrVq1OHToEDVq1GDQoEFcddVVDBky5IzjT548yahRo5g9ezYtW7bkd7/7Ha+99hp//vOfAahVqxYrVqzgP//5D+PGjeOtt9464/nDhw9n8uTJPPbYY+zevZvdu3cTHx/P2rVrs4+59NJLWbJkCSLCW2+9xTPPPJP9JZ5lwIAB/P73v+fgwYPUrFmTiRMnZieSdevWZX8OWeLj43n55ZcB58q+VatW9OrVi4EDB3LTTTdRqVIlhg0bRmxsLC+//DJBQUFMmjSJjz/+ON/fe6dOnYrUOyg/FaJEwO0LYewRGHuEI/ds4pQGs6TWYGe7MSUgdxLIknwi7bzOm5KSQmxsLHXq1GHv3r3069ePY8eOsXjxYoYOHUpsbCx//OMf2b17NwBz5szhjjvuAJwr0WrVqrFw4UJ+85vfEBERQeXKlRk8eHB29USTJk2IjY0FoHPnztlX8B06dOCGG25gwoQJBAXlfz05dOhQvF5vge9hzpw5DB06lFq1agFQo0aNAo/ftGkTTZo0oWXLloBTzbJgwYLs/YMHDz4r3pyGDRuWXbKYPHnyWYkGYOfOnQwYMID27dvz7LPPnlGCyBISEsKgQYP45JNPOHDgAAkJCQwYMKDA2LM88sgjLFu2LPuKfuDAgQBER0cTExPD7NmzWblyJUFBQcTExOR7Hn+tOV8hSgQ5VatZh48vuJdpSZHEnEqncmiF+whMMSnoCv6Sp+aQlEfnhHqRYQDUiAgpUgkgS1YbwYkTJxgwYACvvvoqo0aNIjIykpUrVxb5fLmFhoZm3/d6vdlVQ1988QULFizgs88+44knnmDNmjV5Pj8iIiL7flBQ0BkNoSdPFk8bXVbMXq83z7aJevXqUbNmTVavXs2kSZN4/fXXzzrmnnvu4f7772fQoEHMmzcv39LMiBEjePzxx1FVrrnmmuzuxG3btmX58uV07Ngx+9jly5cTHx+f/bhZs2bccccd/OEPfyAqKiq7ZJFVPRQdHc2IESMKfK8JCQlnnPNclWiJQEQiReQTEdkoIhtEpJuI1BCRWSKy2f1ZvbjjaHb5PSxObcanK3cV90sZA8CYAa0ICz7zyjgs2MuYAa38cv7w8HBeeuklnnvuOcLDw2nSpEl2lYKqsmrVKgAuu+wyXnvtNcCpTjpy5Ajdu3dn+vTpnDhxguPHjzNt2jS6d++e72tlZmbyyy+/0Lt3b55++mmOHDnCsWPHqFKlCkePHs33eY0bN2bFihUArFixgm3btgHQp08fPv74Yw4ePAjAoUOHAPI9X6tWrUhMTGTLli0AvP/++/Ts2bNIn9fw4cN55plnOHLkCB06dDhr/5EjR6hXzxlj9N577+V7nl69erF582ZeffXVM760R48ezZNPPpldIklMTOTFF19kzJgxgJNIs67mN2/ejNfrJTIyEnBKNF9++SWTJk3Kt30AYMqUKXzzzTeFJgtflHTV0L+Br1W1NdAR2AA8CMxW1RbAbPdxsYprEEn/2kdInfus34pWxhTk2rh6PDm4PfUiwxCcksCTg9v7dUBjXFwcHTp04KOPPuKDDz7g7bffpmPHjrRr144ZM2YA8O9//5u5c+fSvn17OnfuzPr16+nUqROjRo2iS5cudO3ald///vfExcXl+zoZGRnceOONtG/fnri4OO69914iIyO5+uqrmTZtWnZjcW7XXXcdhw4dol27drzyyivZVTvt2rXj4YcfpmfPnnTs2JH7778fgOuvv55nn32WuLg4tm7dmn2eSpUq8c477zB06FDat2+Px+Ph9ttvL9JnNWTIECZOnMiwYcPy3D927FiGDh1K586ds6us8uLxeBgyZAgHDx48IxnFxsby9NNPc/XVV9OyZUtatmzJa6+9RqtWTuJ///33adWqFbGxsYwcOZIPPvgguwotMjKSbt26ER0dTdOmTc94vazG+BYtWjBhwgTmzJlDVNRZ68wUmZTUF6GIVANWAk01x4uKyCagl6ruFpG6wDxVLfAyKT4+Xs93YZofJj9Nl/X/YutvPqNZxx7ndS5TMcXHx9sCScYnDz74IEuXLmXmzJmEhISUyGvm9fcpIstV9ay6pJIsETQB9gPviEiCiLwlIhFAtKrudo/ZA0Tn9WQRuU1ElonIMn/0m2074A9kBoXRLHHyeZ/LGGMK8tRTTzF37twSSwJFVZKJIAjoBLymqnHAcXJVA7klhTyLKKr6pqrGq2q8P4pClavVwNNhKKydAinJ530+Y4wpq0oyEewEdqrqUvfxJziJYa9bJYT7c19JBZQWdzOkneCHT8/uNWCMMRVFiSUCVd0D/CIiWfX/lwHrgU+Bm9xtNwEzSiqm4AadWB/akYyUIyX1ksYYU+qUdCf6e4APRCQE+Bm4GScZTRaRW4HtQN7N+MWkzQPzEE/FGFdnjDF5KdFEoKorgbxGP1xWknHkJB4PmpnJzm0baNCsXaDCMMaYgLFLYeCH9x+m9vie7N9rq5aZssWmofafokxDDc5U1DfccAPt27cnJiaGSy+9lGPHjtG7d29mzpx5xrEvvvgid9xxB4mJiYSFhREXF0ebNm3o0qUL7777rp/fSdFZIgAu6HodoZLGTzPfDHQoprw7ugfeuRyOnj2b5bmwaaj9pyjTUIMzOC86Opo1a9awdu3a7LWrc88gCs4solkjgJs1a0ZCQgIbNmxg4sSJvPjii7zzzjt+fS9FZYkAaNA6no0h7WiwbRIZGRmBDseUZ/OfgR1LYP7Tfj+1TUNdstNQ7969O3saCnCmvggNDWXIkCF88cUXpKamZn+Ou3btynPajqZNm/L888/z0ksvFfbrLV6qWuZunTt3Vn9L+PwN1Uer6oq5U/1+blM+nfV3+L8rzr4tfdPZd+q46pt9VcdGqj5a1fn5376qKyY4+48dOPu5PoiIiFBV1fT0dB0yZIh+9dVXqqrap08f/emnn1RVdcmSJdq7d29VVR02bJi+8MIL2c9JTk7WZcuWaUxMjB47dkyPHj2qbdu21RUrVui2bdvU6/VqQkKCqqoOHTpU33//fVVVrVu3rp48eVJVVQ8fPqyqqo8++qg+++yz2bHddNNNeuWVV2p6enqe+9u1a6fbtm3TtWvXaosWLXT//v2qqnrw4MHs53/88cdnnO/jjz/WlJQUrV+/vm7atElVVUeOHJn9nho1aqQvvfSSqqq++uqreuutt571mT3//PP6yCOPqKrqrl27tGXLlqqq+s477+hdd92lqqqHDh3SzMxMVVX973//q/fff/9Z50lISNCoqCi96KKL9OGHH87+vFVVr7zySp0+fbqqqj755JP6l7/8RVVVt23bpu3atTvjPIcPH9ZKlSqddf7zldf3JLBM8/hOtRKBq13fG0mmCik/Tgh0KKa8OrIDsmZXUYXkHed9SpuGOnDTUMfGxvLzzz8zZswYDh06xIUXXsiGDRuAMxeYyVktlBctBfOd2RzMruDQcD5r/wpPLFO+PXyC+tXDAx2SKWtu/iL/fad+hZPJnB44r87j5m6dd0TNgp+fD5uGOv+YS2Ia6qzEOXjwYDweD19++SVt2rThmmuu4b777mPFihWcOHGCzp075xtvQkICbdq0Obc36ydWIsihz2UDSCWEST/+EuhQTHkz/xnQXIvTaKbf2gpsGuqSn4Z60aJFHD58GIDU1FTWr19Po0aNACdB9O7dm1tuuaXA0kBiYiKjR4/mnnvuKVL8/maJIId6kWGMrr+RHt/fSlra+a0cZcwZdv4AGalnbstIdbb7iU1D7Tt/TEO9detWevbsmf05xMfHc91112XvHzFiBKtWrTorEWzdujW7++iwYcO49957ufnmm4sUv7+V2DTU/uSPaajzs2bWe7RfdC9JV7xHvS7XFstrmPLBpqE2pVlpnYa6TGjXawRaOZp6Wz4KdCjGGFMiLBHk4gkOQeJGoj/N5Pi+bYEOxxhjip0lgjykdhyJAiumvRjoUEwp5vF4rC3JlEppaWl4ijCZpiWCPITUaszyRr+nSstLAx2KKcVat27N+++/b8nAlCppaWm8//77tG7d2ufnWGOxMedo3759jB49mo0bN57RP96YQPJ4PLRu3Zpx48ZRu3btM/bl11hsA8oKsC9pG1uXfE636wLbx9eUTrVr12b8+PGBDsOY82ZVQwXYMfsNuq35O5s3rAp0KMYYU2wsERSgxYA7SVcPu+fYmsbGmPLLEkEBqkU3ZH3VS2m/7zOOHT8e6HCMMaZYWCIoRNjFt1FdjrJyZt7zjRhjTFlniaAQzbtewQ5PfXb8tKpUTBdrjDH+ZomgEOLxsqjvdP6WPIiEX5IDHY4xxvidJQIfXN25CREhXj5ZtDbQoRhjjN9ZIvBB5dAgXqg3hwc2Dic5+XCgwzHGGL+yROCjlhf2o5ocR9dNDXQoxhjjV5YIfNQ4ri9Etab6OlvT2BhTvpRoIhCRRBFZIyIrRWSZu62GiMwSkc3uz+olGZPPRCD+Fti1gl/WLgp0NMYY4zeBKBH0VtXYHBMfPQjMVtUWwGz3camU2m4YKYSwfdZrgQ7FGGP8pjRMOncN0Mu9/x4wD3ggUMEUJKRydTb0eYN2bboFOhRjjPGbki4RKPCNiCwXkdvcbdGqutu9vweIzuuJInKbiCwTkWX79+8viVjz1KbHYKpH1Q3Y6xtjjL+VdCK4VFU7AZcDd4lIj5w71Rm6m+fwXVV9U1XjVTU+KiqqBELN37p5n7D4uWFkZNgc9MaYsq9EE4GqJrk/9wHTgC7AXhGpC+D+3FeSMZ0LPbKTi4/OJOH7bwIdijHGnLcSSwQiEiEiVbLuA/2BtcCnwE3uYTcBM0oqpnPVqv8tHCOMU9+/FehQjDHmvJVkY3E0ME1Esl73Q1X9WkR+BCaLyK3AdmBYCcZ0ToLDqrKuzpXE755B0u4k6tWtF+iQjDHmnJVYiUBVf1bVju6tnao+4W4/qKqXqWoLVe2rqodKKqbzUfeyOwmVNDZ9/WagQzHGmPNSpEQgIvEiMtyt2smq7ikNXVBLXHSLznxX+XJm/uIlNd0ajY0xZZdPiUBEokVkCfAD8CGnu3g+DzxXTLGVeulXvcSkE52ZtX5voEMxxphz5muJ4AVgL1ATOJFj+8c4jb4VUo+WUbSopqyeNyXQoRhjzDnztVrnMuAyVT3sNvZm2Qo09HtUZYTXIzwT9TUxv3zItsQradK4SaBDMsaYIvO1RBAGpOaxPQo46b9wyp5G/e4gWDKon2ilAmNM2eRrIlgAjMrxWEXEizMn0Gx/B1WW1GjUDhp3J3jle5BpjcbGmLLH10TwV+APIjILCMVpIF4PXAI8VEyxlRnpnW6G5B2smGulAmNM2eNTIlDV9UB7YDHwDVAJp6E4TlW3Fl94ZYO37VUckmoc3Tg30KEYY0yR+TwGQFX3AI8WYyxllgSFEnT3D/SsWSfQoRhjTJH5Oo7gbhG5MY/tN4rInf4Pq+yp6iaBkycrdNu5MaYM8rWN4M/AL3lsTwTu81cwZV3CxMfZ/1QHjqVYMjDGlB2+JoL6OBPC5bbT3WeAavVa0oC9LPvmo0CHYowxPvM1EewBYvPY3gk44LdoyrgmFw/mgNSk8prxOGvsGGNM6edrIvgQeElE+olIsHvrD7wIfFBs0ZUx4g1mV7PhxKevYO261YEOxxhjfOJrIngUWATMxJlr6ATwFU530v8rntDKpuYD7yBdPeyb+3qgQzHGGJ/41H1UVdOAESLyCBCHs67wSlXdXJzBlUXhtRoyteGD/GdbbTodT6V6REigQzLGmAIVaT0CVd2sqpNV9WNLAvlrc8UdbEmvzZQVOwMdijHGFMrnAWUiMhxnFtLa5EogqjrIz3GVaW3qVuXGukmELJhO5iVv4vFI4U8yxpgA8XVA2bPABKAxkAwczHUzuVxfdxe/S5vM1g3LAx2KMcYUyNcSwe+AEar6SXEGU5407387uvEVWvwyBdrFBzocY4zJl69tBB5gZTHGUe5UiqyDtB0EKz9AU08U/gRjjAkQXxPBm8BZcw2ZgqXFjoKTR5j9yRuBDsUYY/Lla9VQJPBbEekHrAbScu5U1Xv9HFe5ENysO5siOlMpxBvoUIwxJl++JoK2nK4aap1rn82lkB8RWo2ZQ6tAx2GMMQXwdUBZb3+9oLvE5TIgSVWvEpEmwESgJrAcGKmqea2PXGalpZ5i05ofiel8aaBDMcaYsxRpQJmf/AnYkOPx08ALqtocOAzcGoCYitWGd++h8afXsXPPvkCHYowxZ/E5EYhIbxF5U0S+FpE5OW9FOEd94ErgLfexAH2ArG6p7wHX+hx9GVG3+++oLCdZ/81bgQ7FGGPO4uuAslE4k8xVAXoB+4HqONNQry/C670I/BXIdB/XBJJVNd19vBOol08Mt4nIMhFZtn///iK8ZOBFtb6EHSHNaPTzJFLTMgIdjjHGnMHXEsFo4G5VHYHTY+ghVY3DGW18zJcTiMhVwD5VPaehtqr6pqrGq2p8VFTUuZwicERI6fA7WpHID4u+CXQ0xhhzBl8TQVPgW/f+KaCye/8VYJSP57gEGCQiiTiNw32AfwORIpLVaF0fSPLxfGVK88tu4QSVOPyjDc42xpQuviaCgzjVQuB8Uce492sCYb6cQFUfUtX6qtoYuB6Yo6o3AHOBIe5hNwEzfIypTPGGVWXahR9w78Fr2bLPp0KUMcaUCF8TwXdAf/f+ZJzVyt4BPgJmnWcMDwD3i8gWnMTy9nmer9Tq36M7QV4vHy7Ja/lnY4wJDF8HlN0NVHLvPwmk41T1TAb+WdQXVdV5wDz3/s9Al6KeoyyKqhLKE/W+p/mKsZwcuJhKIT7PAm6MMcXG1wFlh3Lcz8Tp+2/OQafm9Wm+bxNJ6+dSL7ZfoMMxxhifu49miEjtPLbXFBHrD1kEzXrdiFaqRr0tEwMdijHGAL63EeS3xFYoUK6mgyhuEhKBdByBrp/Bwb22lKUxJvAKrBoSkfvduwrcLiI5u7t4ge7AxmKKrdxKjb2JkKWv88O0V7j89qcCHY4xpoIrrI3gHvenAL8HclYDpQKJwO3+D6t8C6nbltWt7qVp6/6FH2yMMcWswESgqk0ARGQuMFhVD5dIVBVAhxGPBzoEY4wBfGwjUNXeuZOAiDQXkUr5PccULnHdUr4d/y9UbUkHY0zg+Npr6F8icpN7X0TkW+AnYLeIdC3OAMuzE8s/otfWZ1mzYUPhBxtjTDHxtdfQDcAm9/7lQEfgImA8YK2d56hx/7sIkkyS5rwZ6FCMMRWYr4kgGmeKaIArgMmq+gPwMhBXHIFVBOF1WrC5Shdi98/g0NETgQ7HGFNBFWXSuUbu/f7AbPd+EPmPMTA+CLv499SVQ/z4jQ0wM8YEhq+JYArwoYjMAmoAM93tscCWYoirwqjfZTA7vQ1Y99MmMjOt0dgYU/J8TQT3Ay/hrEbWT1WPu9vrAq8VR2AVhjeYH6/8kpeO9GDx1oOBjsYYUwH52n00XVWfU9U/qWpCju0vqKotxHueLm9fj+phQXy58MdAh2KMqYDyHVAmIp2Alaqa6d7Pl6qu8HtkFUilYC//jZpMo8RZ7D28lujqVQMdkjGmAiloZPEyoA6wz72v5N0wrDjzDpnz0KDrtUR99glHts6E+KGBDscYU4EUlAiaAPtz3DfFKDruCljQkGrrxlsiMMaUqHwTgapuz+u+KSYeL5mdbsIz93E2rFlOm/adAx2RMaaC8HWKiaYicr+IvCIiL4vIfSJipQQ/y4i9gXS87Jhj7e/GmJJT6FKVIvIXnHWKvTjtBQJEAU+LyAOq+kLxhlhxBFery/ZrJtGrzSWBDsUYU4EUWCIQkUuBZ4BngShVrauqdYDawHPAsyJi31p+1CiuH6GVwgMdhjGmAimsaugOYLyqPpxrAfuDqvoQMAG4szgDrIjWfP4qC54cRGp6ZqBDMcZUAIUlgouAdwvY/657jPGjsIyj9Dg1n7v++QJLH+nK1U9OYXpCUqDDMsaUU4UlgjrAzwXs34ozzYTxo021r+KUBvOXzP9xoWxi6PEPeWjqGksGxphiUVgiCANOFbA/FQj15YVEpJKI/CAiq0RknYg85m5vIiJLRWSLiEwSkRDfQi+//jV/H99mxtFKduIRZah3AZXTDvLszE2FP9kYY4qo0F5DwJUiciSffZFFeK1TQB9VPSYiwcBCEfkKZ0K7F1R1ooi8DtxKBZ/IbldyChJ0eiZSD5ncEzSVR5NvCWBUxpjyypdE8HYh+32aO1mdhXmPuQ+D3ZsCfYDfutvfA8ZSwRNB+2op9Dm5EnEn9AiVdIZ557OPGkz/sTXXXtgssAEaY8qVAquGVNXjw83neYZExCsiK3HGI8zCaWNIVtV095CdQL18nnubiCwTkWX79+/P65By48W6s5Bc+TWIDEYHTab3NwNg2f8gI83WLzDG+IWv6xH4hapmqGosUB/oArQuwnPfVNV4VY2PiooqrhBLhaYn1xEq6WdsC5JMtHoTKtduDJ/fx6kXO/H0M4/y055fAxOkMabc8KVqyO9UNVlE5gLdgEgRCXJLBfUB6xpz+8I8NwvgVYXNs0j/eiyDTs2jXvVHAKddoU7VSng8tnKoMaZoSqxEICJRIhLp3g8D+gEbgLnAEPewm4AZJRVTmSQCLfsTcfdC2v15OhGhQaQf2sHhl3oy9rnn+HrNbpzmGGOM8U1JVg3VBeaKyGrgR2CWqn4OPADcLyJbgJoU3jhtADweCK8BgBzbQ+NKJ/jH8cep/fFVPPz8K3y7fq8lBGOMT6QsflnEx8frsmXLAh1G6ZKRRsaKCZya/RThJ/fwXUYMz9X+F3/q34ZeLaMQsSojYyo6EVmuqvG5twekjcAUA28w3gtvJjx2BBnL3iFq8yYO7M7g5nd+ZGC9FG4Y2JNLm9eyhGCMOYtPiUBEquP07++NM/PoGVVKqlrb75GZcxNcCW+3O2jdDeakZzJ7zkwGLP49CyZeRMYfnieoTttAR2iMKWV8LRGMB9rhDPjai4+DyExghQR5uLzHxaR7R9NjyX/wvH4x6TFD+cfRqxnWvycx9aoFOkRjTCngayLoBfRU1RXFGIspDpWqEXTZ3+GiO2DRi3iWvsn9GV+yJeVHoBrpGZkEeUt0OIkxppTxNRFspYQHnxk/i6gJ/R/H0+0uIpJWEt/8AlBl9lsP87mnF7cO7Epsg8hAR2mMCQBfv9z/BDwpIh1FxOcpJUwpVKUOwa0HOvf3rKbfnjd4Jmkki9+4h3vfns3apPzmFzTGlFe+JoItOFNSrwBSRSQj5634wjPFqm5HPHf9QFDbq7gj6DOe+OUGvvnP/dz17kLW7bKEYExF4dM4AhFZAFQHXiePxmJVnVIs0eXDxhEUg73rSPv2n6TsSKDXqXEcOilc3i6aP/drRas6VQIdnTHGD853HEE80EVV1/o3LFNqRLcj+IaPCD55hLmZ4by7YAN9vh/FR5u6c/t9/6BOTethZEx55WvV0HqganEGYkqJStWoFh7Mn7rVoHXDaMYGvUud8ZfAivGMX7iFn/cfK/wcxpgyxddE8HfgeRHpKyLRIlIj5604AzQBUq0+wbd8CSOnQeXa8Ok99Jh1FTOXrAp0ZMYYP/O1auhL9+c3nNk+IO5j60lUHolAsz7QtDds+oroNZ9yQ98LAfh+6WJm/BLBXX1asHz7YZ6duYldySlcEBnGmAGtuDYuz/WFjDGlkK+JoHexRmFKNxFofQVhra8gDOD4AeK/GUzltLr834qhfEdHamQeZmLIy9ydfC8PTU0FsGRgTBlRaK+hrIXmgd+p6qYSiaoQ1msowDLSYfUk0uc+SdCvv7A0szXJmRH0865gQsZlPJJ+C/Uiw1j0YJ9AR2qMySG/XkOFthGoahrQBJtfyGTxBkHcDQTdu4L/S7uZprKLAUHL8Ygy1LuAKJJJSk7h15NpgY7UGOMDXxuL3wP+UJyBmDIoKIQ5VQYxK6Mzaeo0E3nIZELIE4wJmsiCpT8CkJ6RaYvkGFOK+dpGEAHcICL9gOXA8Zw7VfVefwdmyoa/96xB768XEuwOMA+VdJqzixZBu5C5n8EvfZkXcSXjtjVi0h3dqRYWHOCIjTG5+ZoI2uBMLwHQNNc+u9SrwC4/OJ4MD2f+FXiC8LS/Dqo3hhXv0ffoLFKjb6NamNNm8MHS7TSLqkzXJjVsoRxjSgGfEoGqWq8hk7edP+DVM9sCvJoG+9bB4Degxxj46WuuqNcJgLT1X3DB1//mlVN9eLhGF4Z3acR1nepTs3JoIKI3xlDEpSpFpBLQHOf6b6uqniyWqEzZcfvCgvd7g6DNVdkPg08l06vSFnrrUnan1OWdmb24YmYv4tu15LddGtKtaU08HislGFOSfJ10Lhj4F3A3EIIzkOwU8DLwsNuzqMRY99EyLv0UbPgMlv0Pti9if6XG9D35DEdOptOwRjjDL2zADV0bEhkeEuhIjSlXznfSuaeBEcDtOGMKALoDT+L0PBrtjyBNBREUCu2HOLd9G4g6uoelDXswa/V2ms8cyaRvO3O8zUNEhtfl4LFTRIaH4LVSgjHFxtcSwR7gFlX9Mtf2K4G3VLVuMcWXJysRlFOHfoZPboVdKyA4HGKu47HdXVmZ2ZRpd14S6OiMKfPOeUCZqxrOcpW5bQUizyMuY06r0RRumwu3zYP2Q2HtFB7dczd3t3JmPE3PyOT+ySuZtX4v6RmZgY3VmHLE1xLBEmC5qt6Va/trQKyqdvPhHA2A8UA0TmPzm6r6b3f20klAYyARGKaqhws6l5UIKoiTR2DTV9BhOIiQPOMhvln1M2+n9CK5SguGxTdgWHwDGtQID3SkxpQJ+ZUIfE0EPXBmIE0ClribLwIuAC5X1UK6joCI1AXqquoKEamCMzDtWmAUcEhVnxKRB4HqqvpAQeeyRFBBfX4/mjAByTjFT6HteO1YT77M7EKX5hcwoktD+raJJiTI10KuMRXPeSUC9wQXAHcBrd1NG4D/qOqucwxoBvCKe+ulqrvdZDFPVVsV9FxLBBXY8YOw6kOnx9Ghn0m44HruOjiMXUdOUqtyCNd1rs99fVtSKdhmRjcmt/NOBH4OpjGwAIgBdqhqpLtdgMNZj3M95zbgNoCGDRt23r59e0mFa0qjzExIXADVGpBRvSkrF88kbNHTTJH+/P3+vyBBIazemUyrOlUIDbKkYAycYyLwdfUxVT1UhEAqA/OBJ1R1qogk5/ziF5HDqlq9oHNYicCcZeOX8OUY+HUnVI4mreONXLGwGfEdO/Dk4PbZh01PSLJFdEyFda7jCA5Q+FxC6sN5soIIBqYAH6jqVHfzXhGpm6NqaJ8v5zLmDK2vgJYDYPMsWPY/ghY9zxdhNfj5IueCYf2uX7nrw+XsPJxCWobzJ52UnMJDU9cAtoiOqdgK+wIvaI6hgcCfgHRfXsit9nkb2KCqz+fY9SlwE/CU+3OGL+cz5iweL7QaCK0GIsk7CNm/idYX1IDMTC74dDhXJjdgUkZPQHkl5GXuTr2X/WmRPDtzkyUCU6EVuY1AROKAZ3FGFr8BPK6q+3143qXAd8AaIKsT+N+ApcBkoCGwHaf7aIFVTVY1ZIrk+AH45BbYNp809ZKkNWko+/kgow//l34rAG+M7EyPFlGEhVh7gim//NFrqAnwBDAUmAr8TVXzGmRW7CwRmHNx/b/e58oT07jROxsRSNUgLjn1EqcI5hTB1K5RjQVjeiMinEzLsJ5Hptw555HFIlJTRP4NbATqABer6vBAJQFjztX1l/chyOMhFecLXoH7QqYzqf2PrK98FxOrv4Gs+QRNSeay5+bz7MyNgQ3YmBJSYBuBiDwMjMEZ8XuNqn5dEkEZUxyube4lI2gB3szTq6kN987H2/lmqJpGvU1fwtSvwRPMu5GX8EvDNwDYefgEt7z7I/3aRtO/bR3a16tmU2WbcqWwxuLHgRRgJ3CniNyZ10GqOsjfgRnjd/OfwStnVoV6RWHLt3D1i3Dl87DzR2Tj57TISKNFmzoARHx1L6Myq/L2/Ha8OrcO0VVD6dc2mn5t69CtaU0bzWzKvMISwXhsKUpTXuz8ATJSz9yWkepsB/B4oGFX55Yl9QTVj23ht0cT+G0I/FqlGYu8F/HKis5MWFKHKqFB9GwVRf92dejfNtraFUyZFJCRxefLGotNiUv+BTZ9CRs/h8RFpF3xAt9VGciiVRvZtWkZC1Jb8uMjlxMeEsSanUeoVSWEutXCAh21MWcoVVNMnC9LBCagThwCbwiEVoYf/gtfjiYjNBJvq4HQ+kqGfBtORlBY9hoKOw+foF5kGM5QGmMC53xXKDPGZAnPMfNK7A1QpS7ejV/AT1/B6olMDgpj9fVLATh2MpU+4+ZTp1olt7E5ms6NqhPktXYFU3pYIjDmfISEQ5urnFtGOmxfhGf3KmKbNwKg0pSRfBd9kG/1Qv77fRveXliT6uHBXNYmmn5to20QmykVrGrImOI0/xlYOxX2bwDgSGRbvqh0NU/t6cyvJ9MJDfLQvUUUD17eiua1qwQ4WFPeWdWQMYHQ86/O7eBW2Pg51TZ8zm9bVWLoH/qxbPNOMuc8xftJ7YgIdpb5mLtxH1v2HeOmixsTEuSx2VJNibBEYExJqNkMLvmTc1MlWIRuYTvhwGQuyUyDt16AVlew60h7PkxqyO+7N2F6QhJ//WQV1TIOMTHkZe5OvpeHpjrdXy0ZGH+yFitjSlpW76FGF8Nff4br3obGl8LaKdywdTSfj6iNiPDO14sJyTjOvUHTuFA2cU/QVFLSMnjma5v6wviXtREYU1qkn4Lti6BpbxBh8t8Hca33O7woXlFOajDdT/2b/UTSt000fdvUpk+b2tSuUinQkZsy4pwnnTPGlJCgUGjWJ7vE8E34FWzW+njcwf2hpDE55DEiQrxs2P0rD05dQ5cnZjNhibNsa2amUhYv7EzgWSIwppS6rvdFNJPd2TVJIlBfDjDu8jos/GtP1rZ8iw9aL+aSGkcAmPfTPi59ei5b9h0LYNSmLLLGYmNKqcsPjifDwxmzfXk8Hi4/OB6O/5XK6Ye5ZNcrkPgKRMfQqm5f+tTpSoMaztQWr87dwtqkI1zWJpreraKoWTk0MG/ElHqWCIwprXb+gFfTztjk1TRnkrwqdeC2eZC8AzZ8Dhs+pd7Kf/P4jZ9AkBeSdxB1dD0rtofx1do9eAQ6NazuDmSrTbOoyjblhclmjcXGlBdH90B4TfAGw7ePwcLn0Wr1OdhgALPpyoSkaNbsPg5Ao5rhXNY6mis71KFzoxqFnNiUFzagzJjyrkqd0/cvvgdqtUDWf0qtDRMYnvE2w6s3YdcDi5i96QCzN+xlwtLtnEhNp3OjGqgqX63dwyXNa1EtLDhw78EEhCUCY8qj8BoQ+1vnduoobP4Gju7lguoRjLwogpEb7yItrh7Hm10BaS1Zt+8Ud36wgmeHdGBofAOST6SSfCKNxrUiAv1OTAmwRGBMeRdaBWKuO/04PRWq1iN401dErpkIX1amXYv+zBw8gjptnVLFjJW7ePTTdTSLiqBv22j6tommU8PqeG2JznLJEoExFU1QCAx+w0kIiQtgw2fIxi9o1bwvhAfDr7sZxHyCBsby1ZaT/G/hNt6Y/zPVw4Pp3bo2fdtE06NlFJVD7eujvLDGYmMMZGY4t6AQWPY/+Pw+8ARBkx6ktLiS7zxd+GpbJnM37SP5RBrBXqFXq9q8ObJzdu8jmyCv9LPGYmNM/jxe5wbQaRTU6QgbZsD6Twn7+i/09wTRf8wW0kM6kLBtH7N+OszJtIzsJHDdfxaxOukIaRnOhWVScgoPTV0D2AR5ZUGJJQIR+R9wFbBPVWPcbTWASUBjIBEYpqqHSyomY0wePB6o39m59X0M9q6DpOUQVp0g4MLv7+DCk0egzdVwMIxT1RqzZtevpGUoURzmlZCXuTv1XvanRfLPL9YzoF0dW3ynlCvJKSbeBQbm2vYgMFtVWwCz3cfGmNJCBOrEQOebTm9r1gdU4dux8HInQt/szrU6B+CMmVIBDhxLpf3YmVz98kIembGW6QlJbD943OZEKmVKtI1ARBoDn+coEWwCeqnqbhGpC8xT1VaFncfaCIwpBZJ/gY2fw/pPeXVXc9491pXvQv9MJUnjlAbR59Q4UiLqM6JLA1ZsT2bVzmROpGYAUDMihNdHdubCxjU4ejINjwgR1vhc7EprG0G0qu527+8BogMZjDGmCCIbwEV3wEV3UG/FTu6f/me8ZAIQKunMD72PQzUuonZ4fxjxWzIiovlp71FW7DhMwo5kGlQPB+DjZTv55xfrWfq3vkRVCWXz3qN4PELTWhE2DUYJCXQiyKaqKiL5Fk9E5DbgNoCGDRuWWFzGmMJd2yKIjKAFeDMzsrd5PB5qZx6E2Y9Bq8vxVqlDm1NraCPruaF3L6jqTILXpUkNxgxoTVQV5/GLszfzxerdVAsLJq5hJHENqtOpUSQdG0RStZKNei4OgU4Ee0Wkbo6qoX35HaiqbwJvglM1VFIBGmN8MP8ZvLmu4zweLzTpDjd/ARFRzsaNX8CSV537VS6Apj2JadKTmB7Ds5/3l34t6dGiFgk7klmx4zDzf9qPqtNc0aJ2ZTo1rE63ZjW5JtZ6I/lLoBPBp8BNwFPuzxmBDccYc052/gAZqWduy0h1tleufXrbgCfgwlth23z4eT78NBN2LIHYEc7+hAk0rRRJ05hLGX6hU/L/9WQaq35Jzk4MX63dwy+HT2QngoenraFzo+oM7lS/JN5puVSS3Uc/AnoBtURkJ/AoTgKYLCK3AtuBYSUVjzHGj25f6NtxIlCzmXOLvwUyM+HYHmefKsx/2plaWzxwQRw06UnV1lfSvUU83Vs4pYrMTOXoyXQA0jMyWb3zSHa10rFT6Vz76iI61o+kUyOnWqlVnSp5To1hA+BOs5HFxpjSIz0VkpbBz/OcEkPSMujyRxj4L8hIg+9fgcY94ILY0wPgAFVFREhKTuHRGWtZsSOZQ8edEkpEiJcOORJDXMNIvtt8gIemriEl7XSbRliwlycHty+VycBfSSu/XkOWCIwxpdepo5B+CiJqwa6V8GZPZ3toNWh8KTTtBW0HnTkFN05i2HHoRHZ1UsKOZDbs/pX0TOf7rkZEMIeOn7noD0B01VBeGB6LVwSvR/B4BK8IHhE8HvC6jxvXiiDY6+FIShrHTqVzQbVKiAjHT6WTmp7pPC/ruR6yz3cuvaCmJyT5LWlZIjDGlH3H9sG2BU6JYdt8pxrpps+dRum962D3amjaE6pecNZTU1IzWJN0hIQdh3nqq42czzff9w/1oW61MP797WZe+PYntv7rCrwe4eFpa/hg6Y58nyfiJIVN/7wcr0f4x2frmbtpH3NH9wJg9Mer+G7zfjeBOMkj6XBKdgLLqV5kGIse7FOkuEvrOAJjjPFd5drQfohzAzi0Daq6V8XrpsOCZ5z7NVs4CaFpL2g5ELzBhIV46dKkBl2a1GD899tJSk456/TRVUN5cXgcmapkZCoZqmRmKpkKGZmavb16eAgAfdvWpm61SmQ1QVzZoS7Na1fOcSynz5Xj+VnHd2oUSUTo6Squ9vWq4RXJft0MVbYfPJHnR7Erj/jPlZUIjDHlQ2Ym7FvntC38PA+2LwZvEPx1m9OesPFLCA6Dhhcxfe0hHpq6hsppB7LnRjoWXLNUthFc8tScPJOWlQiMMSY3jwfqtHduF9/tNDwf3na6UXnO47BvPXhDubZBF1q2i2P3T8u4MGMTD0Z8iveq50tdEgAYM6BVnm0EYwYUOhuPzywRGGPKp6AQiMrxZXnrN7D9e3cMwzzaJr5EW/GAKNfJPDg2Gfb2h6g2TlIpJbKSU3F2dbWqIWNMxTTtDljzMWSmgSfY+QkQXstpfG7SA1peDlXrBjZOP8qvaqj0pD1jjCkpR/fAuqmnv/wz0yAoFAY8Cc37OqOdP78Pflni7D+cCCs/giNJAQu5OFnVkDGm4pn/DGjmmdtU4eAWZz1nVTi4Faq4EyJv+hq+fsC5X6OZU1po2tPpkRQcVrKxFwNLBMaYiqeguZHA6fBfq/npfV1ucwawbZvvjGNY8wkkvA8PbHf2b/nWGfnc6BKoVLVk3oMfWSIwxlQ8vs6NlMXjcVZqqxMD3e6CjHQ48BOEVnb2L3rJSRLidedI6gHNL3OSRxlgicAYY4rKGwTRbU8//u1k2PmjU1rYtgAWvwS7V55OBD++DdHt4IJOTm+mUsYSgTHGnK/gSm5Po+7Aw3DqGJw44Ow7+St8OQY0A4IjoFE3p8TQ+ipnFtZSwBKBMcb4W2jl09VGlarCmC2QuPB0iWHWIxAc7iSCo3th/QwnOUS1ctonSpglAmOMKW7hNZxZUtsOch4f3eN0VwXYvgi+GuPcrxztJIQmPaDtNVCp2ulzHN0Dn9wMQ9493ZvJTywRGGNMScs5bXbMYKjX6XRpYdsCZ6Bbs8ucRLBtgZMEts5xxjfMfxquet6v4VgiMMaYQKve2Ll1+t3pMQzV3CkkEibA6kmnj135AfR8wK+lAhtZbIwxpUnuMQzXvgZtrnG6poIzEG7+0359SUsExhhTmh3fD5tnOr2OwBn4tvIDp5HZTywRGGNMaZbndBj+LRVYIjDGmNKssOkw/MAai40xpjQr6nQY58BKBMYYU8FZIjDGmArOEoExxlRwlgiMMaaCs0RgjDEVXJlcvF5E9gPbz/HptYADfgynuJWleC3W4lOW4i1LsULZivd8Y22kqlG5N5bJRHA+RGSZqsYHOg5flaV4LdbiU5biLUuxQtmKt7hitaohY4yp4CwRGGNMBVcRE8GbgQ6giMpSvBZr8SlL8ZalWKFsxVsssVa4NgJjjDFnqoglAmOMMTlYIjDGmAquwiQCEfmfiOwTkbWBjqUwItJAROaKyHoRWScifwp0TAURkUoi8oOIrHLjfSzQMRVGRLwikiAinwc6lsKISKKIrBGRlSKyLNDxFEREIkXkExHZKCIbRKRboGPKj4i0cj/TrNuvIvLnQMeVHxG5z/3/WisiH4lIJb+du6K0EYhID+AYMF5VYwIdT0FEpC5QV1VXiEgVYDlwraquD3BoeRIRASJU9ZiIBAMLgT+p6pIAh5YvEbkfiAeqqupVgY6nICKSCMSraqkf9CQi7wHfqepbIhIChKtqcoDDKpSIeIEkoKuqnutg1WIjIvVw/q/aqmqKiEwGvlTVd/1x/gpTIlDVBcChQMfhC1Xdraor3PtHgQ1AvcBGlT91HHMfBru3UnuFISL1gSuBtwIdS3kiItWAHsDbAKqaWhaSgOsyYGtpTAI5BAFhIhIEhAO7/HXiCpMIyioRaQzEAUsDHEqB3KqWlcA+YJaqluZ4XwT+CmQWclxpocA3IrJcRG4LdDAFaALsB95xq93eEpGIQAflo+uBjwIdRH5UNQkYB+wAdgNHVPUbf53fEkEpJiKVgSnAn1X110DHUxBVzVDVWKA+0EVESmX1m4hcBexT1eWBjqUILlXVTsDlwF1uNWdpFAR0Al5T1TjgOPBgYEMqnFuFNQj4ONCx5EdEqgPX4CTbC4AIEbnRX+e3RFBKuXXtU4APVHVqoOPxlVsVMBcYGOBQ8nMJMMitd58I9BGRCYENqWDu1SCqug+YBnQJbET52gnszFEa/AQnMZR2lwMrVHVvoAMpQF9gm6ruV9U0YCpwsb9ObomgFHIbX98GNqjq84GOpzAiEiUike79MKAfsDGgQeVDVR9S1fqq2hinOmCOqvrtysrfRCTC7TCAW83SHyiVPd9UdQ/wi4i0cjddBpTKDg65jKAUVwu5dgAXiUi4+/1wGU7boV9UmEQgIh8B3wOtRGSniNwa6JgKcAkwEudqNatr2xWBDqoAdYG5IrIa+BGnjaDUd8ssI6KBhSKyCvgB+EJVvw5wTAW5B/jA/VuIBf4V2HAK5ibXfjhX2KWWW8r6BFgBrMH57vbbdBMVpvuoMcaYvFWYEoExxpi8WSIwxpgKzhKBMcZUcJYIjDGmgrNEYIwxFZwlggpMROaJyCsBem0VkSGBeG1flPb4/MGdxXKsD8fNFZHflUBIeb12gb8Hd9bb60oypvLIEkE55Q7y+o87hfEpEdkrIrNFpF+OwwYDDwUqRn8TkU7uF0f3fPZPEpHFJR1XQfL7ohORd0vDFNkiciXQAPggx7ZEN24VkRR3yukx7kCnkvY48JSI2HfZebAPr/yagjMVwa1AS+Aq4CugZtYBqnrInd20zBGRoNxfPO6MrSuBW/I4viZwLTbjaFH9CXhXVTNybf8HzkDCNjiTof0LCMSEeF8CVXCmiTDnyBJBOeRO99AdeFBVZ6vqdlX9UVXHqerEHMedUTXkXun9XUTecBfp2CkiY3Kdu6WIzBeRkyKySUSuEJFjIjLK3d/YvVKMz/W8wor4T7nnS3HjeCbnwhsiMtatyhglIluBU0BeM1u+BQx1J+zL6Ub3OZNEZKCIfCcih0XkkIjMFJE2BcTm03sSkXoiMtE972ER+UJEWuR33qIQkcEistr9fA65v4PoHPuvFmd20pMisk1EnhBnMrWs/bVFZIb7/O0iclayzOM1o3DmuPksj91HVXWPqiaq6lvAapzpL7Ke28x9vT0iclxEVogz4V/O8xf695ZHTA+IyAERuQicyQ5xksGIwt6PyZ8lgvLpmHsbJEVfxeg+nCHsnYCngWfEXWXKLX5PA9KBi4BRwKNAqB9iPo5zJd8GuBNnHqCHcx3TBPgtMBToCJzM4zwfAF5geK7ttwKTVPU4TgJ5EafE1As4AnyW84uzqEQkHGeyvZNAT6AbznTB37r7zpmI1MGZIO89nM+nB/B+jv0DcN73K0A7nM9xCGdO7/Au0Bzni/1a4HdA40Je+lKc5Jnv3Ebi6OXGlZZjV2WcEmg/nN/VFGCqiLTOdYp8/97yeJ1xOFNY9My16NEPOJ+5OVeqardyeAOuw1mI5yTOHEvjcFZfynnMPOCVHI8TgY9yHbMZ+Lt7fwBOEqiXY//FOPPlj3IfN3Yfx+c6jwJD8nucR/y3A1tyPB6L80UT7cN7nwAszvH4Qvf1uuZzfASQgTPd81nx+fKecL58N+NO2+Ju8wIHgWEFxJrn54Dzxf25e7+Te1yjfM6xAPi/XNuuxbkYEJyqQQUuybG/kfuexxYQ25+B7XlsT8RJEMeAVPfcKcDFhfxelmT9Lfny95bj8xkOvAP8lNdngDOFdCYQVNL/Z+XlZiWCckpVp+DMW341zpXZxcASEflbIU9dnevxLqC2e781sEvdaZFdP+KHBV5EZIiILHSrEo4BLwANcx22U32bKvgtoFuOq89bgLXqTo/sVlt8KCJbReRXYC9O6Tj36xVFZ5wSy1G3quwYTkmjOtDsPM4LsAr4FlgrIlNE5A632ibnaz+c9brua3+Ik+Dq4FytZ+JcOQOgzkpcha1wFUbepS6A53EmleuJUxJ6TFWzG+LFmTX1GXHW3T7sxhTP2Z9xQX9vWcbhlNwu1bxXEEvBSXh+W8O3orFEUI6p6klVnaWq/1DVi3Gmth5bSBVIWq7HStH+TrKSQnZDrjhrK+TLre+dCMzESVxxwN9xlrzM6biPMcwHtgC3iDMt9gjc5RNdnwNRwB+Bru7rpQP5fS6+vCcPTkN1bK5bS+CNAmI9ClTLY3skTiJBnXrw/u5tNU4112YR6ZjjtR/L9bodgBY4K4ZlKeoMkwdwElleDqrqFlX9Hqf0OVpEeufYPw6nCu//cJJFLE4iyv0Z+/L3NgsnoeU3A28N4KSeXi7VFFFQoAMwJWo9zu+8Ek6Rvqg2AheIyAWqmnU1Gc+Z/7hZXzx1c2yLLeS8lwBJqvp41gYRaXQO8QHOGsoi8j+cHi8bca5s33fPWxOnZHOnqs51t3Wi4P8FX97TCpyEc0CLtk7vJpwr+uxEJc5C6h1xqkOy3xNOFd/3IvIPYB1Olckq97Vbq+qWvF5ARDbi/I66AIvdbQ1xSowFSQCiRKSWqh7I7yBVPSxOp4MXRCTOjfVSYLxbMsVtq2qGU71TVF/iTBP9sYioqr6Xa38MzmdgzpGVCMohEakpInNE5EYR6SAiTURkKM46vbP13Je9nIXzxfWeiHR0r+Sfx7madip0VVNw6oIfEJF2InIxztVhQX4C6onIDSLSVETu4Px7gbwH1HJfe7qqHnS3H8a50v2DiDQXkZ7A6+57yJOP7+kDnCqmGSLS0/3Me4jIc4X0HHoep+Rylzg9smJx5pmv4f5ERC5ye9dc6H6BD8Lp25+16Ms/gN+KyD9EJEZEWrtVbc+48W8CvgbeEJFu7mu8i1OlUpAEnDWoLy3kOID/AK1wSgHg/E5/I87YjvY47TbnXHWjzvoWQ4HX5ezBbd1x3p85R5YIyqdjOF9cf8KpJlmH04PkQ87uTeMzVc0EfoPTS+gHnC/bJ3CSQM665KyuiT/iVIv8vZDzfgY8i9OTZzVOT5NHzjVO95y7cK4kq5Nj7ID7HobjVJ2sBV7Fqb44VcgpC3xPqnoCpzfPzzhr327E+Xyq4ySf/OL8CLjZvS3D+UKrA3RXZ8UvcKqILsGp0toMPAc8rqoT3HPMBK4EeuP8Xn7AWSt4R46XGgVsA+bgdAf9EKexNl9uldT/gBsKOs49dh9OqWus27vsfpwk8h1OG9US9/45c5PBMJyE9jtwuuzitH+9U9BzTcFsYRpzXtx66pU4PWrK0oLwxgciUhun5HGhqm4LdDy5icizQDVVDcRgtnLD2ghMkYjIb3AabTfjdKt8ntP11KacUdV94gw+a4hToiht9lF41aMphJUITJG4RfK/49RRH8YZi3Cfj906jTGlkCUCY4yp4Kyx2BhjKjhLBMYYU8FZIjDGmArOEoExxlRwlgiMMaaC+3+W48JOu4xlcwAAAABJRU5ErkJggg==",
      "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",
    "# 计算 Frobenius-norm 误差\n",
    "for i in range(RANK):\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",
    "# 画图 \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(list(range(1, RANK+1)), err_subfull, \"o-.\", \n",
    "        label = 'Reconstruction via VQSVD')\n",
    "ax.plot(list(range(1, RANK+1)), err_SVD, \"^--\", \n",
    "        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": [
    "---\n",
    "### 案例2：图像压缩\n",
    "\n",
    "为了做图像处理，我们先引入必要的 package。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:17.225100Z",
     "start_time": "2021-03-09T03:50:17.151352Z"
    }
   },
   "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": 167,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.337792Z",
     "start_time": "2021-03-09T03:50:17.231211Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsUlEQVR4nO3dbWxcVXoH8P8fx46NHSW2A4mdOC+khGAqNomsQEVAdOkugS8EtOLlA2Il1KwqkIq0+wFRbZdWVctWBUSliioUtNmKQtgFRFqhdilaKUVasWvSxCREJCE4L7bjJMRJnMRxYufph7nZOtF9jsczd2Zsn/9Psjw+z5yZx9fz+I7v8TmHZgYRmf6uqXQCIlIeKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFil6KR/D7JUZJnxnzcXem85EozKp2ATBu/MbO1lU5CfDqzT3Mku0n+iGQXyVMkN5OsrXReUn4q9jg8DGAdgKUAbgXw/bQ7kVxL8mTgI3TmXkXyOMk9JH9MUu8aJxn9QOLwj2bWCwAk/x3AyrQ7mdknAOYU8PhbAfwhgAMAbgGwGcAIgL8r4LGkRHRmj8ORMbfPAWjI8sHNbL+ZfW1ml8zscwB/DeB7WT6HFE/FLr9H8s6rrqhf/XFnng9lAFjKXGXi9DZefs/M/gcFnPVJ3gdgm5n1k1wB4McAfpF1flIcndklC/cA6CJ5FsCHAN4D8LeVTUmuRi1eIRIHndlFIqFiF4mEil0kEip2kUiUdeitsbHRWltby/mUmbrmmvTfjV47AJD+cHPo4ujo6GhBsUIuuBaafyh26dKlCbUD4dxnzPBfqlVVVW7MU2gek11vby8GBgZSfzBFFTvJdQBeAVAF4F/M7IXQ/VtbW7F58+ZinrKiZs6cmdpeX1/v9qmpqXFjFy5ccGOnTp1yY6dPn3Zjw8PDqe2hgq6rq3Nj3vcMhAtwaGgotX1wcNDtMzIy4saam5sLinmFe+7cObfPxYsX3dhk98gjj7ixgt/Gk6wC8E8A7gPQDuAxku2FPp6IlFYxf7OvAbAv+b/oCwDeBvBANmmJSNaKKfYFAA6N+fpw0nYFkhtIdpLsHBgYKOLpRKQYJb8ab2YbzazDzDoaGxtL/XQi4iim2HsAtI35emHSJiKTUDHF/jsAN5JcSrIGwKMAtmSTlohkreChNzMbIfk0gP9CbujtDTPblVlmIpKposbZzexD5KY0isgkp3+XFYmEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4mEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEUTvCkOwGMAhgFMCImXVkkZSIZK+oYk/8sZkdz+BxRKSE9DZeJBLFFrsB+BXJz0huSLsDyQ0kO0l2DgwMFPl0IlKoYot9rZmtBnAfgKdI3nX1Hcxso5l1mFlHY2NjkU8nIoUqqtjNrCf5fBTA+wDWZJGUiGSv4GInWU9y1uXbAL4LYGdWiYlItoq5Gj8PwPskLz/Ov5nZf2aSlYhkruBiN7P9AL6VYS4iUkIaehOJhIpdJBIqdpFIqNhFIpHF/8ZH49SpU6ntPT09bp/z58+7sZkzZ7qx2bNnu7Frr73WjY2Ojqa29/f3u33OnDkz4ccDgAsXLrixwcHBCedx+PBhN1ZVVeXG2tvb3diKFStS2xcuXOj2qaurc2NTmc7sIpFQsYtEQsUuEgkVu0gkVOwikSj71XgzK/dTZubIkSOp7bt27XL7dHd3u7Ha2lo31tTUVFA/78p6b2+v2+fYsWNuzBuBAIChoSE35uWYzKVIdfDgQTe2f/9+N7ZkyRI3tn79+tT2e++91+2zYMECNzaV6cwuEgkVu0gkVOwikVCxi0RCxS4SCRW7SCQ0EWYCzp07l9p+6NAht09XV5cbO336tBsLTQoJDYd5Q17XXXed26ehocGNhSbJjIyMuLGlS5dOqH28PELHOBTzJt4MDw+7faYrndlFIqFiF4mEil0kEip2kUio2EUioWIXiURZh95IoqamppxPmam2trbU9jvvvNPts2jRIjfmzaIDgJ07/Z20+vr63Nj8+fNT2++44w63z2233TbhxwOAa67xzxXe8OCePXvcPlu3bnVjt9xyixsLDeetXbs2tb2lpcXtM5Vfo6FZheOe2Um+QfIoyZ1j2ppIfkRyb/JZ27OKTHL5vI3/GYB1V7U9C+BjM7sRwMfJ1yIyiY1b7Ga2FcCJq5ofALApub0JwPps0xKRrBV6gW6emV3+w/EIcju6piK5gWQnyc6BgYECn05EilX01XjLrTPlrjVlZhvNrMPMOhob9ae9SKUUWuz9JFsAIPl8NLuURKQUCh162wLgCQAvJJ8/yLfjVF5w0tt2KbRA4Zw5c9zY6tWr3dhDDz2Ud15jedtNXbx40e0TynHu3LlurL6+3o1dunQptf3Eiasv//y/vXv3urHQDLvQ8Oa8eel/YVZXV7t9pvJrNCSfobe3APwGwE0kD5N8Erki/w7JvQD+JPlaRCaxcc/sZvaYE7on41xEpIT077IikVCxi0RCxS4SCRW7SCS019sEeDOKQsM4M2fOdGOhPdtCC0SGZmWdPXs2tT20SGUoxxkz/JfI6OioG/P2gQstshmKeUN5QHj2XSjmmcqv0RCd2UUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIbeMhAa3gkNXVVVVbmx0FBTSF1dXWp7KMdQHqHhtdBiJL29vanthexTB/gzDoHwrL3QsKhnOr5GAZ3ZRaKhYheJhIpdJBIqdpFIqNhFIlHWq/FmVvBV5sksdDU7NGkl1C+05lroGHpX42fNmuX2CV1xD109P3DggBvbt2/fhB8vtN5daJ0/b505wD8eIVP5NRoaSdCZXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIaCLMBHi5e2vTAeHhtVC/0HZNoWE5b+JHQ0OD28fbMgoADh8+7Ma++uorN7Z///7U9tA6c83NzW5syZIlBfULra/nmcqv0ZB8tn96g+RRkjvHtD1Psofk9uTj/tKmKSLFyudt/M8ArEtpf9nMViYfH2ablohkbdxiN7OtAPytN0VkSijmAt3TJLuSt/mN3p1IbiDZSbIztNiBiJRWocX+KoBlAFYC6APwondHM9toZh1m1tHY6P5OEJESK6jYzazfzEbN7BKA1wCsyTYtEclaQUNvJFvMrC/58kEAO0P3ny68obLQ8FooFhriCQ2HhdaT84blQn2Gh4fdWE9Pjxv74osv3Jg39BZaS2758uVubNWqVW4sNFvOm3VYyLZQU924xU7yLQB3A5hL8jCAnwC4m+RKAAagG8APSpeiiGRh3GI3s8dSml8vQS4iUkLxvZcRiZSKXSQSKnaRSKjYRSJR9llvoZlek503VBYaQgt9v6F+oZlthSxGGVpEcWhoyI0dP37cjYWG5fr7+1Pb29ra3D5NTU1u7Prrr3djoUUlve879HOZyq/REJ3ZRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEWYfeSE7p2UaFDL2FhI5FocfJG0Y7duyY2+fQoUNurLu7240dOXLEjV24cCG1PbTnXGtrqxsL7ecWGooMDSt6pvJrNDRsOHW/KxGZEBW7SCRU7CKRULGLRELFLhIJTYTJQKETYUJXfUPbFs2Y4f/Yzpw5M6F2APjyyy/dWOhqfGgrp9ra2tT20ISW+fPnu7HZs2cXlEdoLT/PdHyNAjqzi0RDxS4SCRW7SCRU7CKRULGLRELFLhKJfHaEaQPwcwDzkNsBZqOZvUKyCcBmAEuQ2xXmYTOLcpvW0NBbKBbaGsrbtggIDw2dPHkytf3rr792++zZs8eNhSa7hIYAvSG2lpYWt09okkzI6OioGytkDbrpKp8z+wiAH5pZO4DbATxFsh3AswA+NrMbAXycfC0ik9S4xW5mfWa2Lbk9CGA3gAUAHgCwKbnbJgDrS5SjiGRgQn+zk1wCYBWATwHMG7OT6xHk3uaLyCSVd7GTbADwLoBnzOyK/0+03B+mqX+cktxAspNk54kTJ4pKVkQKl1exk6xGrtDfNLP3kuZ+ki1JvAXA0bS+ZrbRzDrMrCO0CYCIlNa4xc7cZcvXAew2s5fGhLYAeCK5/QSAD7JPT0Syks+stzsAPA7gc5Lbk7bnALwA4B2STwI4AODhkmQ4iXjDaKGtlUJCQ1ehWV7e+m6AP0utq6vL7bNt2zY3dvRo6hs2AOHtmtrb21Pbb7rpJrdPQ0ODGwvNXhseHnZj3rBcaNhzuhq32M3sEwDeoOQ92aYjIqWi/6ATiYSKXSQSKnaRSKjYRSKhYheJRNkXnCx0q6TJwBtiCw29hb7fUL/Q0NDFixfdWE9PT2r7jh073D6hBSdDM9FWrFjhxtauXZvafuutt7p9QotserP5AGBwcNCNeUKLfU7l12iIzuwikVCxi0RCxS4SCRW7SCRU7CKRULGLRKKsQ29mVvAMscnAW6QwtHhhKBYa4hkZGXFjQ0NDbswbhgoNT4XyaG5udmPLly93Y97stnnz/AWNQnu2hfaqC80CrK6uTm0vdEh0sgt9Xzqzi0RCxS4SCRW7SCRU7CKRULGLRKLsE2Gmo0InVYSuIn/zzTdurK+vz42dOnUqtb22ttbt09ra6sZuuOEGN7Zs2bIJP2Z9fb3bx8sdCE/+CfF+Ntr+SUSmLRW7SCRU7CKRULGLRELFLhIJFbtIJMYdeiPZBuDnyG3JbAA2mtkrJJ8H8KcAjiV3fc7MPixVopNBIRNhvO2HgPDEj+PHj7sxb4snAOjv709tb2xsdPssXrzYjd18881urK2tzY15a9eF1tYrdEJRiDf0Fhouna7yGWcfAfBDM9tGchaAz0h+lMReNrN/KF16IpKVfPZ66wPQl9weJLkbwIJSJyYi2ZrQexmSSwCsAvBp0vQ0yS6Sb5D03yeKSMXlXewkGwC8C+AZMzsN4FUAywCsRO7M/6LTbwPJTpKdAwMDxWcsIgXJq9hJViNX6G+a2XsAYGb9ZjZqZpcAvAZgTVpfM9toZh1m1hG6SCQipTVusTN3GfR1ALvN7KUx7S1j7vYggJ3ZpyciWcnnavwdAB4H8DnJ7UnbcwAeI7kSueG4bgA/KEF+k4o3XFPoemahLY12797txvbu3evGvJl0CxcudPuE1pJbtGiRG2toaHBj58+fT20PDUWGhsNCs/ZCM+JmzEh/icc46y2fq/GfAEg7MtN6TF1kuonvPwtEIqViF4mEil0kEip2kUio2EUioQUnSyw0LOcNTwHAoUOH3NjBgwfdWFNTU2p7aIZae3u7G5szZ44bC/H+W9IbChtPaOgtNIzmzbILDfNN5e2fQnRmF4mEil0kEip2kUio2EUioWIXiYSKXSQSGnrLQGgmV2g/t8HBQTfW29vrxg4cOODGvGGo6upqt8/s2bPdWF1dnRsbGRlxY96wYmgILRQLLVQZGirzYjEuOBnfdywSKRW7SCRU7CKRULGLRELFLhIJFbtIJDT0NgHeME5o6G14eNiNnT171o0dO3asoNjcuXNT20OzzUILR4b6hb43bxHI0JBXTU2NGwsNr4UWnPSGIjX0JiLTlopdJBIqdpFIqNhFIqFiF4nEuFfjSdYC2ApgZnL/X5rZT0guBfA2gGYAnwF43Mz8WR+RKsU2Q6FJId5kktAV9/r6ejcWWkMvNAoxNDTkxjyFXo2PcSunQuRzZh8G8G0z+xZy2zOvI3k7gJ8CeNnM/gDAAIAnS5aliBRt3GK3nDPJl9XJhwH4NoBfJu2bAKwvRYIiko1892evSnZwPQrgIwBfAThpZpcnNB8GsKAkGYpIJvIqdjMbNbOVABYCWANgRb5PQHIDyU6Snd5a4iJSehO6Gm9mJwH8GsAfAZhD8vIFvoUAepw+G82sw8w6Ghsbi8lVRIowbrGTvI7knOR2HYDvANiNXNF/L7nbEwA+KFGOIpKBfCbCtADYRLIKuV8O75jZf5D8AsDbJP8GwP8CeL2EeU4K3uSJ0NppoeGp5uZmN7Z48eKCHnPRokWp7aFtnELr04WeKzTk5Q2VhSbPhGKhCTmhdfK8YcpQ7qHveSobt9jNrAvAqpT2/cj9/S4iU4D+g04kEip2kUio2EUioWIXiYSKXSQSDM1qyvzJyGMALu9dNBfA8bI9uU95XEl5XGmq5bHYzK5LC5S12K94YrLTzDoq8uTKQ3lEmIfexotEQsUuEolKFvvGCj73WMrjSsrjStMmj4r9zS4i5aW38SKRULGLRKIixU5yHckvSe4j+Wwlckjy6Cb5OcntJDvL+LxvkDxKcueYtiaSH5Hcm3wu+UofTh7Pk+xJjsl2kveXIY82kr8m+QXJXST/PGkv6zEJ5FHWY0KyluRvSe5I8virpH0pyU+TutlM0l+ON42ZlfUDQBVya9jdAKAGwA4A7eXOI8mlG8DcCjzvXQBWA9g5pu3vATyb3H4WwE8rlMfzAH5U5uPRAmB1cnsWgD0A2st9TAJ5lPWYACCAhuR2NYBPAdwO4B0Ajybt/wzgzybyuJU4s68BsM/M9ltunfm3ATxQgTwqxsy2AjhxVfMDyK3SC5RptV4nj7Izsz4z25bcHkRuJaQFKPMxCeRRVpaT+YrOlSj2BQAOjfm6kivTGoBfkfyM5IYK5XDZPDPrS24fATCvgrk8TbIreZtf1oUDSS5BbrGUT1HBY3JVHkCZj0kpVnSO/QLdWjNbDeA+AE+RvKvSCQG53+zI/SKqhFcBLENuQ5A+AC+W64lJNgB4F8AzZnZ6bKycxyQlj7IfEytiRWdPJYq9B0DbmK/dlWlLzcx6ks9HAbyPyi6z1U+yBQCSz0crkYSZ9ScvtEsAXkOZjgnJauQK7E0zey9pLvsxScujUsckee6TmOCKzp5KFPvvANyYXFmsAfAogC3lToJkPclZl28D+C6AneFeJbUFuVV6gQqu1nu5uBIPogzHhLnVH18HsNvMXhoTKusx8fIo9zEp2YrO5brCeNXVxvuRu9L5FYC/qFAONyA3ErADwK5y5gHgLeTeDl5E7m+vJ5HbIPNjAHsB/DeApgrl8a8APgfQhVyxtZQhj7XIvUXvArA9+bi/3MckkEdZjwmAW5FbsbkLuV8sfznmNftbAPsA/ALAzIk8rv5dViQSsV+gE4mGil0kEip2kUio2EUioWIXiYSKXSQSKnaRSPwffyIxG32ku9UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUkklEQVR4nO3df4xV5Z3H8fe3yK/KYBGmgIgDAroCRSATUiq2rm4t2jTaZNNWN103caXp1mSbdP8w3WTrbvaPdrNt081u3NDV1DZdqbaauhvaqtQGTS064AhYrfwoWJAfQwUB5YcD3/3jHrID3u8zM+f+muH5vJLJ3Hm+99zzzJn7nXPv+d7neczdEZHz3/ta3QERaQ4lu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULLLoJnZfDP7hZkdMLP3fFDDzC42s8fM7G0z22lmt7ein3I2JbuU8S7wMHBnEP8P4CQwGfgL4D4zm9ekvknA9Am684uZ7QD+HfhLoAP4OXCHux9vwL5mA1vc3fq0XQgcBOa7+2tF2w+A3e5+T737IAOnM/v56TPAcmAmsAD4q2p3MrNlZnYo8bWsxL6vAHrPJHrhJUBn9ha7oNUdkIb4N3d/A8DM/gdYWO1O7v4s8IE673sccPictreAtjrvRwZJZ/bz094+t9+hkoDNchQYf07beOBIE/sgVSjZM2Zm15rZ0cTXtSUe9jXgAjOb06ftauDl+vRaytLL+Iy5+zOUOOubmQGjgVHFz2MqD+cn3P1tM3sU+Ccz+2sqbyFuAT5St45LKTqzSxkdwDH+/2x9DPhdn/jfAGOB/cBDwBfdXWf2FlPpTSQTOrOLZELJLpIJJbtIJpTsIploault0qRJ3tHR0cxdimRl586dHDhwwKrFakp2M1sOfAcYAfyXu389df+Ojg5+/etf17LLIalSdq4uVe1Ixco+ZpnHSynbx2Y9Xn+PWe99DXUf+Uj8cYbSL+PNbASVoYw3AXOB28xsbtnHE5HGquU9+xJgq7tvd/eTwCoqn5QSkSGolmSfBvyhz8+7irazmNkKM+sys66enp4adicitWj41Xh3X+nune7e2d7e3ujdiUiglmTfDUzv8/OlRZuIDEG1JPsLwBwzm2lmo4DPAY/Xp1siUm+lS2/u3mtmdwO/oFJ6e0Ajm0SGrprq7O6+Glhdp76ISAPp47IimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimahpRRgz2wEcAU4Bve7eWY9OiUj91ZTshT919wN1eBwRaSC9jBfJRK3J7sATZrbezFZUu4OZrTCzLjPr6unpqXF3IlJWrcm+zN0XAzcBXzKzj557B3df6e6d7t7Z3t5e4+5EpKyakt3ddxff9wOPAUvq0SkRqb/SyW5mF5pZ25nbwI3A5np1TETqq5ar8ZOBx8zszOP8t7v/vC69EpG6K53s7r4duLqOfRGRBlLpTSQTSnaRTCjZRTKhZBfJRD0+G5+Nt956a1DtAO97X/z/dObMmWHs+PHjYWzEiBFh7MiRI1Xbd+3aFW5z+vTpMHby5Mkwtn379jB29OjRqu2p45E6jtu2bQtjF1wQP43nz59ftX3ZsmXhNrNmzQpjw5nO7CKZULKLZELJLpIJJbtIJpTsIplo+tV4d2/2Lutmx44dVdufe+65cJt9+/aFsblz55bqx6lTpwa9v9///vfhNqm/yeHDh8PY66+/HsaKMRPvceGFF4bbpK7Gd3d3h7He3t4wduutt1Ztv+KKK8JtLr/88jA2nOnMLpIJJbtIJpTsIplQsotkQskukgklu0gmNBBmEKLS0Pr168Ntfvazn4Wxt99+O4yNGTMmjKVKb9GgkPHjx5faV2pAzrhx48LYpZdeWrV9ypQp4TYTJkwIY1HZE9Ilu2hG44kTJ4bbnK90ZhfJhJJdJBNKdpFMKNlFMqFkF8mEkl0kE00vvaXmIBvqFi9eXLU9NZIrNbpq3bp1YSw12iy1v6iP119/fbjN7Nmzw1iqZJdalfeiiy6q2j569Ohwm9WrV4exX/3qV2Fs3rx5YeyGG26o2p6a/284P0dT+v2tzOwBM9tvZpv7tF1sZk+a2Zbie1wgFZEhYSD/wr4HLD+n7R5gjbvPAdYUP4vIENZvsrv7WuDNc5pvAR4sbj8I3FrfbolIvZV9czLZ3fcUt/dSWdG1KjNbYWZdZtZ14MCBkrsTkVrVfCXCK3MahfMauftKd+90985JkybVujsRKalssu8zs6kAxff99euSiDRC2dLb48AdwNeL7z8d6IbDecLJqJy0aNGicJtUiee2224LY2XLP1FpKzVCLVUOS42wmz59ehiLlqhKTRy5atWqMBYtJwWwdOnSMHbZZZeFschwfo6mDKT09hDwHHClme0yszupJPnHzWwL8GfFzyIyhPV7Znf36PRT/dMKIjIknZ8fFRKR91Cyi2RCyS6SCSW7SCa01tsgROWw1ISNqUkU29rawtjYsWPD2MmTJ8NYdHxHjhwZbpNy4sSJMJYafRdNpvnqq6+G22zYsCGMRWvHQbq8+cEPfrBq+6hRo8JtTp8+HcaGM53ZRTKhZBfJhJJdJBNKdpFMKNlFMqFkF8mE1nobhN7e3qrtqbJQKhaNDIPyJcrU/iKpUlNq9F2qj9u2bava/vzzz4fbHDp0KIwtXLgwjM2dOzeMRaXP1O+l0puIDGtKdpFMKNlFMqFkF8mEkl0kE02/Gl/mavFQEV2NT129TQ1AueCC+PC/++67YSx1DKOrzKk+puaZSw3IiQa7AKxfv75q+7PPPhtukzpWN910Uxj70Ic+FMai/qeOx3B+jqbozC6SCSW7SCaU7CKZULKLZELJLpIJJbtIJjQQZhCiElWqjJMqr6UGY6TKYan506LHTJXyUlKDdfbs2RPGotLb1q1bw21SSzUtWbIkjLW3t4exSFRGhfTvPJwNZPmnB8xsv5lt7tN2r5ntNrPu4uvmxnZTRGo1kJfx3wOWV2n/trsvLL5W17dbIlJv/Sa7u68F3mxCX0SkgWq5QHe3mW0sXuaHk6Ob2Qoz6zKzrp6enhp2JyK1KJvs9wGzgIXAHuCb0R3dfaW7d7p7Z5kLKSJSH6WS3d33ufspdz8NfBeIL5WKyJBQqvRmZlPd/Uzd5dPA5tT9zxdlllBKlddS5bDUdqlRWVEZMFXKSy1fdezYsTD29NNPh7ForrmJEyeG23z2s58NYwsWLAhjKWXmDTxf9ZvsZvYQcB0wycx2AV8DrjOzhYADO4AvNK6LIlIP/Sa7u99Wpfn+BvRFRBpIH5cVyYSSXSQTSnaRTCjZRTKhUW+DEI2GKjt6LTXyKjVaLlU2ivZXtpT3xz/+MYx1d3eHsR07dlRtv+SSS8JtUss4tbW1hbHUcYyWqEodj/NVfr+xSKaU7CKZULKLZELJLpIJJbtIJpTsIplQ6a0OUhNORqUfSJd/yq4DF0lNUnn06NEw9swzz4SxF154IYxFI+mWLl0abjNv3rwwlhpxePLkyTAWlRVzHPWmM7tIJpTsIplQsotkQskukgklu0gmdDV+EKKr7qkr7qlYI5aGKjNY54033ghjq1fH63+89tprYWz+/PlV26+77rpwm46OjjCWqnikrqxHv3dqm9TfbDjTmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTAxkRZjpwPeByVRWgFnp7t8xs4uBHwEzqKwK8xl3P9i4rrZeVJJJlWpSJZ5ULDWvWqqMFsUOHoz/NC+++GKpWKr/V1111aDaAUaPHh3GUstQpQbJRMfjfC2vpQzkzN4LfMXd5wIfBr5kZnOBe4A17j4HWFP8LCJDVL/J7u573H1DcfsI8AowDbgFeLC424PArQ3qo4jUwaDes5vZDGARsA6Y3Gcl171UXuaLyBA14GQ3s3HAT4Avu/vhvjGvvAGq+ibIzFaYWZeZdfX09NTUWREpb0DJbmYjqST6D9390aJ5n5lNLeJTgf3VtnX3le7e6e6d7e3t9eiziJTQb7Jb5ZLr/cAr7v6tPqHHgTuK23cAP61/90SkXgYy6u0a4PPAJjPrLtq+CnwdeNjM7gR2Ap9pSA+HgVQJKhqFBumRXO+8804YmzBhQhiLRsS9/PLL4TaPPPJIGHv11VfD2NVXXx3Grr322qrts2bNCrcpOy9cmRFsZUfRDWf9Jru7PwtEv/0N9e2OiDSKPkEnkgklu0gmlOwimVCyi2RCyS6SCU04OQipkWhlpMpyZZc7ikpsqfLamjVrwlhbW1sY+8QnPhHGrrnmmqrt48ePD7dJLUOVmpyzzFJZZUcqDmc6s4tkQskukgklu0gmlOwimVCyi2RCyS6SCZXeBiFVKisjtWZbavLFVAnw9ddfr9qemjjyyJEjYezGG28MYx/72MfC2CWXXFK1PTXaLKXs2nc5TiwZ0ZldJBNKdpFMKNlFMqFkF8mEkl0kE02/Gj+cr45GAy5SAydSV4pTsdRAmP37q07kC8CmTZuqtkdX6QFmzJgRxj71qU+FsdQcdKNGjaranppbL9oG0sf4xIkTYSz6m6Wu7g/n52iKzuwimVCyi2RCyS6SCSW7SCaU7CKZULKLZKLf0puZTQe+T2VJZgdWuvt3zOxe4C7gzNKsX3X31Y3q6PkoVf6J5k4D+M1vfhPGfvnLX1ZtT83vdvvtt4ex5cuXh7EpU6aEsagclhoIU2Yuuf5Ex7hsuXQ4G0idvRf4irtvMLM2YL2ZPVnEvu3u/9q47olIvQxkrbc9wJ7i9hEzewWY1uiOiUh9Deo9u5nNABYB64qmu81so5k9YGbx0qIi0nIDTnYzGwf8BPiyux8G7gNmAQupnPm/GWy3wsy6zKyrp6en2l1EpAkGlOxmNpJKov/Q3R8FcPd97n7K3U8D3wWWVNvW3Ve6e6e7d7a3t9er3yIySP0mu1UuW94PvOLu3+rTPrXP3T4NbK5/90SkXgZyNf4a4PPAJjPrLtq+CtxmZguplON2AF8YyA6H89I6Ufmn7Lxq73//+8PY3r17w9hTTz0Vxrq7u6u2X3nlleE2d911Vxi77LLLwtixY8fCWGTs2LFhLFVeS8VS8/VFz7fU32w4P0dTBnI1/lmg2m+vmrrIMKJP0IlkQskukgklu0gmlOwimVCyi2RCyz8NQlSSSS0LlZo4MjX54hNPPBHGNm+OP9IwbVr1YQuf/OQnw22uuuqqMFZ2dFg0eWTq8VLltdRxTI0ePHnyZNX21KSS9V7ma6jQmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTKj0NghlSm+pstCbb74ZxtauXRvGtmzZEsY6Ojqqtk+dOrVqO6TLUPUuUaVGm6WOVSpWdj29iEpvIjKsKdlFMqFkF8mEkl0kE0p2kUwo2UUy0fTSW6qUM9RF5Z9UOenw4cNhLDV6bcOGDWHs4MGDYWzBggVV21Olt9Ros2j0GqRLVGUm4Uyt9ZYqoaWeU6n+R8pOIDrU6cwukgklu0gmlOwimVCyi2RCyS6SiX6vxpvZGGAtMLq4/4/d/WtmNhNYBUwE1gOfd/fqE36d51JXilNX43fu3BnGUivetrW1hbFomafZs2eH25S9qp7arkzVpRGDXaLlplL9O3HiRBgbzgZyZj8BXO/uV1NZnnm5mX0Y+AbwbXefDRwE7mxYL0WkZv0mu1ccLX4cWXw5cD3w46L9QeDWRnRQROpjoOuzjyhWcN0PPAlsAw65e29xl11A9TmMRWRIGFCyu/spd18IXAosAf5koDswsxVm1mVmXan3oSLSWIO6Gu/uh4CngaXAB8zszAW+S4HdwTYr3b3T3Tvb29tr6auI1KDfZDezdjP7QHF7LPBx4BUqSf/nxd3uAH7aoD6KSB0MZCDMVOBBMxtB5Z/Dw+7+v2b2W2CVmf0z8CJw/0B2mFr+Z6iL+j5mzJhwm8mTJ4exOXPmhLEpU6aEsdSAkUWLFlVtnzlzZrhNyvHjx0v1IyrL9fb2Vm2HdAkttfxTqox29OjRqu1l5w0czvpNdnffCLznGeTu26m8fxeRYeD8/BcmIu+hZBfJhJJdJBNKdpFMKNlFMmHNnBPOzHqAM0O9JgEHmrbzmPpxNvXjbMOtHx3uXvXTa01N9rN2bNbl7p0t2bn6oX5k2A+9jBfJhJJdJBOtTPaVLdx3X+rH2dSPs503/WjZe3YRaS69jBfJhJJdJBMtSXYzW25mvzOzrWZ2Tyv6UPRjh5ltMrNuM+tq4n4fMLP9Zra5T9vFZvakmW0pvk9oUT/uNbPdxTHpNrObm9CP6Wb2tJn91sxeNrO/LdqbekwS/WjqMTGzMWb2vJm9VPTjH4v2mWa2rsibH5nZ4Bayc/emfgEjqMxhdzkwCngJmNvsfhR92QFMasF+PwosBjb3afsX4J7i9j3AN1rUj3uBv2vy8ZgKLC5utwGvAXObfUwS/WjqMQEMGFfcHgmsAz4MPAx8rmj/T+CLg3ncVpzZlwBb3X27V+aZXwXc0oJ+tIy7rwXePKf5Fiqz9EKTZusN+tF07r7H3TcUt49QmQlpGk0+Jol+NJVX1H1G51Yk+zTgD31+buXMtA48YWbrzWxFi/pwxmR331Pc3gvEU9w03t1mtrF4md/wtxN9mdkMKpOlrKOFx+ScfkCTj0kjZnTO/QLdMndfDNwEfMnMPtrqDkHlPzuVf0StcB8wi8qCIHuAbzZrx2Y2DvgJ8GV3P2spnWYekyr9aPox8RpmdI60Itl3A9P7/BzOTNto7r67+L4feIzWTrO1z8ymAhTf97eiE+6+r3iinQa+S5OOiZmNpJJgP3T3R4vmph+Tav1o1TEp9n2IQc7oHGlFsr8AzCmuLI4CPgc83uxOmNmFZtZ25jZwI7A5vVVDPU5lll5o4Wy9Z5Kr8GmacEysMpPn/cAr7v6tPqGmHpOoH80+Jg2b0blZVxjPudp4M5UrnduAv29RHy6nUgl4CXi5mf0AHqLycvBdKu+97qSyQOYaYAvwFHBxi/rxA2ATsJFKsk1tQj+WUXmJvhHoLr5ubvYxSfSjqccEWEBlxuaNVP6x/EOf5+zzwFbgEWD0YB5XH5cVyUTuF+hEsqFkF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQT/weYte4bQI3fgAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASZUlEQVR4nO3df4xdZZ3H8feH2h9AK6XbsU5K6ZS2AmWBgYxdQQS2VoP800o2YmOUjWRrNpBogskSN1nZxGR1s2rchLgpP2I1rsiKRlzrQiEutEJopy5Oi0VoaRFK2xki/QHdpbT97h/3NE7rec7M3J/tPJ9XMpk7z/eee74c+plz5zz3nKOIwMzGvzM63YCZtYfDbpYJh90sEw67WSYcdrNMOOxmmXDYzTLhsNuYSfpzSY9Iel3Sn3xQQ9J/S/o/SW8WX7/rRJ92Iofd6vEO8CBwa8Vzbo+IqcXXhW3qyyo47OOMpJ2SvihpQNJ+ST+UNKWZ64iI30XEfcBzzXxday2HfXz6BHADMA+4DPjrsidJukbSvoqvaxro4Z+Kt/m/knR9A69jTfKuTjdgLfGvEfEagKSfAb1lT4qI9cD0Fqz/74DfAoeBTwI/k9QbEdtbsC4bJe/Zx6c9wx4fAqa2c+UR8UxEHIyItyNiNfAr4MZ29mB/ymHPmKQPDTtiXvb1oSatKgA16bWsTn4bn7GIWEcde31JAiYDk4qfp9ReLt6WNB34C+AJ4AhwM3At8PkmtW11ctitHnOBHcN+/l/gZaAHmAh8BbgIOAo8DyyPiBfa3KOdRL54hVke/De7WSYcdrNMOOxmmXDYzTLR1qPxM2fOjJ6ennau0iwrO3fu5PXXXy/9TENDYZd0A/AtYAJwb0R8ter5PT09bNy4sZFVmlmF97///cla3W/jJU0A7gY+BiwCVkhaVO/rmVlrNfI3+2JgW0S8FBGHgQeAZc1py8yarZGwzwZeGfbzq8XYCSStlNQvqX9oaKiB1ZlZI1p+ND4iVkVEX0T0dXV1tXp1ZpbQSNh3AXOG/XxeMWZmp6BGwr4RWChpnqRJ1C5S8HBz2jKzZqt76i0ijki6HXiE2tTb/RHha5KZnaIammePiDXAmib1YmYt5I/LmmXCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2WioTvCSNoJHASOAkcioq8ZTZlZ8zUU9sJfRsTrTXgdM2shv403y0SjYQ/gUUmbJK0se4KklZL6JfUPDQ01uDozq1ejYb8mIq4EPgbcJunak58QEasioi8i+rq6uhpcnZnVq6GwR8Su4vsg8BNgcTOaMrPmqzvsks6WNO34Y+CjwJZmNWZmzdXI0fhZwE8kHX+df4+I/2pKV2bWdHWHPSJeAi5vYi9m1kKeejPLhMNulgmH3SwTDrtZJprx2fhsHDx4sHT8wIEDyWXOOCP9+7S7u7uuPooZkFJvvPFG6fiuXbvqWlfqvxnglVdeGfNyx44dSy5TtR137NiRrFW95mWXXVY6vmTJkuQyCxcuTNZOZ96zm2XCYTfLhMNulgmH3SwTDrtZJnw0fgy2b99eOr5u3brkMnv27EnWFi1a1HBPJ3vttddKx1O9j2T//v3JWtUR8tQsxNSpU5PL7Nu3L1nbtGlTshYRydqyZctKxy+++OLkMj4ab2anNYfdLBMOu1kmHHazTDjsZplw2M0y4am3MRgcHCwdr5p6e/TRR5O1w4cPN9zTySZPnlw6XjXlNXHixGSt6kSeqtrcuXPHNA4wffr0ZG3btm3J2qFDh5K1888/v3R81qxZyWXGK+/ZzTLhsJtlwmE3y4TDbpYJh90sEw67WSY89TYGS5cuLR2fM2dOcplLL700WVu/fn2yduTIkdE3NszixeW320ud/QUwf/78utZVdbbZWWedVTpeNd34xBNPJGtV2+q9731vsnbdddeVjr/vfe9LLjNejbhnl3S/pEFJW4aNzZC0VtKLxfdzW9ummTVqNG/jvwPccNLYncDjEbEQeLz42cxOYSOGPSKeBP5w0vAyYHXxeDWwvLltmVmz1XuAblZE7C4e76F2R9dSklZK6pfUPzQ0VOfqzKxRDR+Nj9pRmuSRmohYFRF9EdHX1dXV6OrMrE71hn2vpG6A4nv5GSJmdsqod+rtYeAW4KvF9582raNTWOosrwULFiSXue2225K1z372s3X1Uc+U19lnn51cZtKkSXX1USW1rQYGBpLL3Hvvvcla1cUoly9fnqylptiqbqE1Xo1m6u0HwNPAhZJelXQrtZB/RNKLwNLiZzM7hY24Z4+IFYnSh5vci5m1kD8ua5YJh90sEw67WSYcdrNM+Ky3Jqi6YOOMGTPa2MmpI3URyOeffz65TH9/f7L2zjvvJGuXXHJJsvae97yndLxq6q1qavN05j27WSYcdrNMOOxmmXDYzTLhsJtlwmE3y4Sn3k5D9UwbtXuqafv27aXjVffFO3DgQLLW29ubrF1++eXJ2rvf/e5kLTfes5tlwmE3y4TDbpYJh90sEw67WSZ8NP40VM/R81Ycca+6RVXqpJbHHnssuczkyZOTtZtuuilZu+iii5K1M888s3R8vJ7sUsV7drNMOOxmmXDYzTLhsJtlwmE3y4TDbpYJT71Z3fbs2ZOsbdiwoXT8hRdeSC5zwQUXJGtLlixJ1qZNm5as5TjFljKa2z/dL2lQ0pZhY3dJ2iXp2eLrxta2aWaNGs3b+O8AN5SMfzMieouvNc1ty8yabcSwR8STwB/a0IuZtVAjB+hulzRQvM0/N/UkSSsl9UvqHxoaamB1ZtaIesP+bWA+0AvsBr6eemJErIqIvojo6+rqqnN1ZtaousIeEXsj4mhEHAPuARY3ty0za7a6pt4kdUfE7uLHjwNbqp5v49OaNenjsuvXry8d7+7uTi7zqU99KlmrusVT1dly9kcjhl3SD4DrgZmSXgW+DFwvqRcIYCfwuda1aGbNMGLYI2JFyfB9LejFzFrIH5c1y4TDbpYJh90sEw67WSZ81ptVnhlWdUumjRs3Jms7duwoHe/p6UkuU3UbpylTpiRrNjres5tlwmE3y4TDbpYJh90sEw67WSYcdrNMeOptnJFUOl41vVZ1z7Zf/OIXyVrqfm6Qvsfa4sXps6GvuOKKZO2MM7xfapS3oFkmHHazTDjsZplw2M0y4bCbZcJH48eZ1FH3o0ePJpcZHBxM1h566KFkrepWTpdeemnp+NKlS5PLzJ07N1mzxnnPbpYJh90sEw67WSYcdrNMOOxmmXDYzTIxmjvCzAG+C8yidgeYVRHxLUkzgB8CPdTuCvOJiHijda3acamTXSA99fbWW28ll3nqqaeSteeeey5ZqzqBZtGiRaXjvb29yWWstUazZz8C3BERi4APALdJWgTcCTweEQuBx4ufzewUNWLYI2J3RPy6eHwQ2ArMBpYBq4unrQaWt6hHM2uCMf3NLqkHuAJ4Bpg17E6ue6i9zTezU9Sowy5pKvAQ8IWIOOFi4lH7Q7H0j0VJKyX1S+ofGhpqqFkzq9+owi5pIrWgfz8iflwM75XUXdS7gdIPWEfEqojoi4i+rq6uZvRsZnUYMeyqHfq9D9gaEd8YVnoYuKV4fAvw0+a3Z2bNMpqz3j4IfBrYLOnZYuxLwFeBByXdCrwMfKIlHdqfqLqeXGpa7ve//31ymbvvvjtZS93GCWD+/PnJ2nXXXVc6fuGFFyaXsdYaMewRsR5ITex+uLntmFmr+BN0Zplw2M0y4bCbZcJhN8uEw26WCV9wcgyqzjZLqZomq3ddVa85MDBQOn7PPfckl9m0aVOy9vbbbydrK1asSNauv/760vEJEyYkl2n3tsqN9+xmmXDYzTLhsJtlwmE3y4TDbpYJh90sE556G4N2TuPUc2YbpM9ue/rpp5PLHDp0KFm7+uqrk7XU9BrAeeedVzreim3o6bXR8Z7dLBMOu1kmHHazTDjsZplw2M0y4aPxLdaKkzT27t2brG3evLl0/OWXX04ukzpyDvCZz3wmWbv44ouTtXe9q/yflk926Rzv2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmRpx6kzQH+C61WzIHsCoiviXpLuBvgOO3Zv1SRKxpVaOnq3qnhY4dO5asrVu3Lllbs6b8f0HVteSWL1+erN10003JWtWNOlP/3fVOoXl6rXGjmWc/AtwREb+WNA3YJGltUftmRPxL69ozs2YZzb3edgO7i8cHJW0FZre6MTNrrjH9zS6pB7gCeKYYul3SgKT7JZ3b7ObMrHlGHXZJU4GHgC9ExAHg28B8oJfanv/rieVWSuqX1D80NFT2FDNrg1GFXdJEakH/fkT8GCAi9kbE0Yg4BtwDLC5bNiJWRURfRPRVHdAxs9YaMeyqHT69D9gaEd8YNt497GkfB7Y0vz0za5bRHI3/IPBpYLOkZ4uxLwErJPVSm47bCXyuBf2Na1XTSfv370/WHnnkkWQtdSunSy65JLnMHXfckayde64PxYwXozkavx4omxz1nLrZacSfoDPLhMNulgmH3SwTDrtZJhx2s0z4gpMtVnWW19GjR5O1tWvXJmsbN25M1s4///zS8Ztvvjm5zPz585O11IUjob5bVLX77LVTpY9TgffsZplw2M0y4bCbZcJhN8uEw26WCYfdLBOeemuxqgtHvvnmm8naz3/+82Rt+/btydqCBQtKx6vOXqt3eq1Ksy84Wa8cp9hSvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU2xhUTRulvPXWW8naU089laxt2LAhWauasjvnnHNKx2fNmpVcpt3TYc12uvffLt6zm2XCYTfLhMNulgmH3SwTDrtZJkY8Gi9pCvAkMLl4/o8i4suS5gEPAH8GbAI+HRGHW9lsp9VzZPfw4fQm2bp1a7K2b9++ZG3q1KnJ2sKFC0vHq27/1Ioj1u289puPuI/OaPbsbwNLIuJyardnvkHSB4CvAd+MiAXAG8CtLevSzBo2Ytij5vjE7sTiK4AlwI+K8dXA8lY0aGbNMdr7s08o7uA6CKwFtgP7IuJI8ZRXgdkt6dDMmmJUYY+IoxHRC5wHLAYuGu0KJK2U1C+pf2hoqL4uzaxhYzoaHxH7gF8CVwHTJR0/wHcesCuxzKqI6IuIvq6urkZ6NbMGjBh2SV2SphePzwQ+AmylFvq/Kp52C/DTFvVoZk0wmhNhuoHVkiZQ++XwYET8p6TfAg9I+grwP8B9LezztHXWWWcla1dffXWyNnPmzGRt9uz04ZGrrrqqdHzevHnJZeqdumr2CSitOKHFt3/6oxHDHhEDwBUl4y9R+/vdzE4D/gSdWSYcdrNMOOxmmXDYzTLhsJtlQu2cgpA0BLxc/DgTeL1tK09zHydyHyc63fqYGxGln15ra9hPWLHUHxF9HVm5+3AfGfbht/FmmXDYzTLRybCv6uC6h3MfJ3IfJxo3fXTsb3Yzay+/jTfLhMNulomOhF3SDZJ+J2mbpDs70UPRx05JmyU9K6m/jeu9X9KgpC3DxmZIWivpxeL7uR3q4y5Ju4pt8qykG9vQxxxJv5T0W0nPSfp8Md7WbVLRR1u3iaQpkjZI+k3Rxz8W4/MkPVPk5oeSJo3phSOirV/ABGrXsLsAmAT8BljU7j6KXnYCMzuw3muBK4Etw8b+GbizeHwn8LUO9XEX8MU2b49u4Mri8TTgBWBRu7dJRR9t3SaAgKnF44nAM8AHgAeBTxbj/wb87VhetxN79sXAtoh4KWrXmX8AWNaBPjomIp4E/nDS8DJqV+mFNl2tN9FH20XE7oj4dfH4ILUrIc2mzdukoo+2ipqmX9G5E2GfDbwy7OdOXpk2gEclbZK0skM9HDcrInYXj/cA6Xsst97tkgaKt/kt/3NiOEk91C6W8gwd3CYn9QFt3iatuKJz7gforomIK4GPAbdJurbTDUHtNzu1X0Sd8G1gPrUbguwGvt6uFUuaCjwEfCEiDgyvtXOblPTR9m0SDVzROaUTYd8FzBn2c/LKtK0WEbuK74PAT+jsZbb2SuoGKL4PdqKJiNhb/EM7BtxDm7aJpInUAvb9iPhxMdz2bVLWR6e2SbHufYzxis4pnQj7RmBhcWRxEvBJ4OF2NyHpbEnTjj8GPgpsqV6qpR6mdpVe6ODVeo+Hq/Bx2rBNVLsq5H3A1oj4xrBSW7dJqo92b5OWXdG5XUcYTzraeCO1I53bgb/vUA8XUJsJ+A3wXDv7AH5A7e3gO9T+9rqV2g0yHwdeBB4DZnSoj+8Bm4EBamHrbkMf11B7iz4APFt83djubVLRR1u3CXAZtSs2D1D7xfIPw/7NbgC2Af8BTB7L6/rjsmaZyP0AnVk2HHazTDjsZplw2M0y4bCbZcJhN8uEw26Wif8HhJAtV2Fb+ucAAAAASUVORK5CYII=",
      "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 = f\"n = {i}\"\n",
    "    plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.966739Z",
     "start_time": "2021-03-09T03:50:18.950606Z"
    }
   },
   "outputs": [],
   "source": [
    "# 然后我们再来看看量子版本的分解效果：\n",
    "\n",
    "# 设置电路参数\n",
    "cir_depth = 40         # 电路深度\n",
    "num_qubits = 5         # 量子比特数量\n",
    "\n",
    "# 超参数设置\n",
    "RANK = 8               # 设置想要学习的阶数\n",
    "ITR = 200              # 迭代次数\n",
    "LR = 0.02              # 学习速率\n",
    "SEED = 14              # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(2 * RANK, 0, -2).astype('complex128')\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": 169,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:19.553293Z",
     "start_time": "2021-03-09T03:50:19.531180Z"
    }
   },
   "outputs": [],
   "source": [
    "# 重新定义量子神经网络\n",
    "def U_theta(num_qubits: int, depth: int) -> Circuit:\n",
    "\n",
    "    # 用 Circuit 初始化网络\n",
    "    cir = Circuit(num_qubits)\n",
    "    \n",
    "    # 搭建层级结构：\n",
    "    for _ in range(depth):\n",
    "        cir.ry()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:59:58.649381Z",
     "start_time": "2021-03-09T03:54:30.126561Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: 3902.7036\n",
      "iter: 10 loss: -109067.9844\n",
      "iter: 20 loss: -124662.9922\n",
      "iter: 30 loss: -136021.7188\n",
      "iter: 40 loss: -143259.4844\n",
      "iter: 50 loss: -147684.2812\n",
      "iter: 60 loss: -150412.4688\n",
      "iter: 70 loss: -151939.5781\n",
      "iter: 80 loss: -152878.7969\n",
      "iter: 90 loss: -153509.0000\n",
      "iter: 100 loss: -153973.1875\n",
      "iter: 110 loss: -154308.6562\n",
      "iter: 120 loss: -154568.4688\n",
      "iter: 130 loss: -154789.6094\n",
      "iter: 140 loss: -154995.1719\n",
      "iter: 150 loss: -155195.5625\n",
      "iter: 160 loss: -155395.0156\n",
      "iter: 170 loss: -155593.9375\n",
      "iter: 180 loss: -155790.7188\n",
      "iter: 190 loss: -155982.7812\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f80fc02ac10>"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYk0lEQVR4nO2da4xVVZqG3w8EQaqAKi5FcZebNzLtpcRLE2W60x3GdESTidEY4w/TdCZtMiY9P4yTjE4yP+zJqPHHxAmOpOmJ0+K0GsnEzDRjOt4S0YIRUNARy0Io6gJFcZGbAt/8OJt0Yc73nqpd5+wDrvdJCKfWd9be66y93zqn1nu+b5m7Qwjx/WdUvQcghCgGiV2IRJDYhUgEiV2IRJDYhUgEiV2IRLhkJJ3NbCWAZwGMBvCv7v4ke35TU5PPnDmzbIxZgKNGlf+dFLUDwNmzZ9lQchGdz8zCPizGxnjmzJkwlscu/fbbb8MYm8dLL700jOWZY/a62PHYPE6YMCGMRXPF5oONkY2DzSM7X8TYsWPDWPS6urq6cPDgwbKDzC12MxsN4J8B/ATAXgAfmtkGd98R9Zk5cybWr19fNsYmY9y4cWXbGxoawj5Hjx4NY+yCjR49OoxFNz67KCzGxvj111+HsdOnT4cxdhNEMLEsXLgwjB0/fjyMRcI9duxY2IfNx5gxY8LYjTfeOOxx9PT0hH0OHToUxtj1jO7TSueLmDt3bhiLfiHdeeedYZ+RfIxfBmCXu3e4+zcAXgKwagTHE0LUkJGIfRaAPYN+3pu1CSEuQGq+QGdmq82s3czaBwYGan06IUTASMTeBWDOoJ9nZ23n4e5r3L3N3duamppGcDohxEgYidg/BLDYzC43s7EA7gWwoTrDEkJUm9yr8e5+2sweBvDfKFlva939E9bHzEJ7gq22RiuP7M+CvLYcW6mPxpHXTpo4cWIYY/ZaX1/fsM+3aNGisM8ll8S3AVtFZm7IkSNHyrYzl4GtdLe0tISx3t7eMBbNI7tmzJ1g15O5JJHlDMT36sGDB8M+0X3Kxjcin93d3wDwxkiOIYQoBn2DTohEkNiFSASJXYhEkNiFSASJXYhEGNFqfK4TBjbPN998E/ZhySnVHAPA7Z+oH7NxWIIPGwezABnRWJjNxywjZmEyIguIHY+9ZpZAw8YfHTNPwhMAnDx5Mlc/Zi1H5EkOY+idXYhEkNiFSASJXYhEkNiFSASJXYhEKHw1Plo5Zau0bLU7D3mTGaIxspVdtuLOYP1YLHI1WFIFOx6bqzwuCVthZqvxbDWbOSgRee8BRt7adVEszwo+Q+/sQiSCxC5EIkjsQiSCxC5EIkjsQiSCxC5EIhRuvUV2TR57jfVhiQKnTp3KdczIkmHWD7OnTpw4Ecb6+/vDWFTfDYgTXtgYme3Jdn3JU3uP2XwskWTSpElhjM1xZGuxhBaWdMPuK1Zfj40/gl2z6B6m26gNewRCiIsSiV2IRJDYhUgEiV2IRJDYhUgEiV2IRBiR9WZmnQCOAjgD4LS7t1V4fmjzsBp0eba6YTFmT+SBnYvBLCOW8cRi0TFZH2ZDsRjbJimKseMxmF3K7LxoPqhFRazIvPcOy3rLo4k891w1fPY/d/cDVTiOEKKG6GO8EIkwUrE7gD+Y2WYzW12NAQkhasNIP8Yvd/cuM5sOYKOZferubw9+QvZLYDXAt60VQtSWEb2zu3tX9n8fgNcALCvznDXu3ububc3NzSM5nRBiBOQWu5lNMLPGc48B/BTAx9UamBCiuozkY3wLgNcyS+ESAP/u7v/FOrh7WHCQFfmLsqGY9cOsmmpbXswGYQUWmcXDsqtYBlU0j0ePHg37MItn1qxZYYxlvUUWW97ioaxg5rx588JYnvuNXU92X7Ettg4fPhzGIlsuz/3B+uQWu7t3APhB3v5CiGKR9SZEIkjsQiSCxC5EIkjsQiSCxC5EIhRecDKyBlixwQhWsJHZa8yeYJZMdL4DB+I8IGY1NTU1hbHJkyeHMfbaent7y7YzO+nQoUNhbGBgIIyxa7Z///5hH4/Za+yasW9mLl68uGw7s+vY8VhRSTbH7HpGRHMIxPOhgpNCCIldiFSQ2IVIBIldiESQ2IVIhEJX4909XO0eP3582C/agojVM2MJLexcbDU+WunevHlz2KerqyuMTZkyJYyxZBf2uj/55JMwFsESYVgCDZvjaIsq1ieaXwDo6OgIY+x6rly5smz7qlWrwj5z584NY3nvOeYmRPPPkqhYolQ4hmH3EEJclEjsQiSCxC5EIkjsQiSCxC5EIkjsQiRCodYb2/4pT00wlmQyduzYMMYsDTaOKOGlp6cn7LNjx44w1t/fH8Y6OzvDGCNKkmHzwWqnMTuJJa4sWLCgbPvll18e9mH21LZt28JYZM0C+WwtloDC+rGEqMsuuyyMRfZmY2Nj2CdPLT+9swuRCBK7EIkgsQuRCBK7EIkgsQuRCBK7EIlQ0Xozs7UAfgagz92XZm3NANYDmA+gE8A97h4XF8sYNWpUaEHkqXXGLCNWp43ZFmxLKZYNFcG2T2JZXix7jc3VnDlzyrbfdtttYZ/58+eHsWnTpoUxZpVFc7x79+6wz8aNG8PYlVdemSu2fPnysu2zZ88O+7AsOlZnjtmUl1wSSy3a/onZfHlqNg7lnf03AL6bJ/gogDfdfTGAN7OfhRAXMBXFnu23/t1vT6wCsC57vA7AXdUdlhCi2uT9m73F3buzxz0o7egqhLiAGfECnZcKVYfFqs1stZm1m1k7+3qoEKK25BV7r5m1AkD2f1/0RHdf4+5t7t7GyjAJIWpLXrFvAPBg9vhBAK9XZzhCiFoxFOvtdwBWAJhqZnsBPA7gSQAvm9lDAHYDuGeoJ4xsBkZko+W1SFhmGxtftIUPO1dra2sYY9bK/fffH8bY+KOCiCwzjBUvbG5uDmPste3bt69s+7vvvhv22bRpUxhjhS+ZjRbZig0NDWEfNh90eyViRbJ+kR3N7o8oi5GNoaLY3f2+IPTjSn2FEBcO+gadEIkgsQuRCBK7EIkgsQuRCBK7EIlQaMFJBstEi+wwZpEwe4qdi9loESwDiVl5hw8fDmPMVowsQCB+bR9++GHYh+31Nm/evDB27NixMBZZTexcbD5OnToVxqZPnx7GWlrKf5ObZaExm5JlUzKrjM1VdMwZM2aEfaL7m91vemcXIhEkdiESQWIXIhEkdiESQWIXIhEkdiESoVDrzd1Da4sV64tsnBMnToR9mAXBMpBYvyjG+jArb9KkSWEs7151kfV21VVXhX0YzKb8+uuvw1hHR0fZdrY/HKt3kDcT7ciRI2XbmV3KbE9mHTI7j9mlUaYay/SL+tCsvDAihPheIbELkQgSuxCJILELkQgSuxCJUOhq/NmzZ8OEALYaH60+sxVVtgpe7fp0bMsodjzmJrB+eZwLttUU29aKbbH11VdfhbEvv/yybHtnZ2fYh610MzeBbV8VrfCzWm0sEYatuLPx50n0YjX+WNJQhN7ZhUgEiV2IRJDYhUgEiV2IRJDYhUgEiV2IRBjK9k9rAfwMQJ+7L83angDwcwD7s6c95u5vVDqWu+ey0SIbitkP7HjMImG2XGSRMDuGJcmwLYiYNcRsuch6mzlzZtiH2WtdXV1h7Isvvghju3btKtvOasktWbIkjN1yyy1hjCUURdea1mojc88sNJagFG3xxGCJMFESEhvfUN7ZfwNgZZn2Z9z92uxfRaELIepLRbG7+9sA4rxEIcRFwUj+Zn/YzLaZ2Voza6raiIQQNSGv2J8DsBDAtQC6ATwVPdHMVptZu5m1DwwM5DydEGKk5BK7u/e6+xl3PwvgeQDLyHPXuHubu7c1NekDgBD1IpfYzWzwN/TvBvBxdYYjhKgVQ7HefgdgBYCpZrYXwOMAVpjZtQAcQCeAXwzlZHmttwhmM5w8eTKMMWslT0ZZ3iw6VnONvTa2BVH0uvOe69ChQ2Fs69atYWz79u1l25kFddNNN4WxG264IYyxen3R/DObjN0DzLZlr43dI9H8szHmoaLY3f2+Ms0vVHUUQoiao2/QCZEIErsQiSCxC5EIErsQiSCxC5EIhRacNLPQ1mB2WJRVlqeIX95zAbGNw8bBrBoWYxl9jY2NYSyy5Zg9xbZkYt96jOw1ANi7d2/Z9qVLl4Z9WOFLZh2yax3NY94tu1g/ZokyGy26r/LY0fS+H/bRhBAXJRK7EIkgsQuRCBK7EIkgsQuRCBK7EIlQqPU2atSoMDOI2QyRtcWKFzKriVkreQo9MsuFZVCx8Ud74gF8b7nx48eXbd+3b1/Yh+3ZtmPHjjD2+eefh7FoTqZNmxb2mTFjRhhj1hvLLIzIWziSzT2zUtk9Eo2FWazR/SHrTQghsQuRChK7EIkgsQuRCBK7EIlQ6Go8q0HHVsGjFUa2ohptjwPwWmEsESZa4Wfb9LAxstVn5k6wVdooqYXV5Iu2agKAPXv2hLETJ06EscgVmDdvXtjnmmuuCWNsFZzN/5EjR8q2s1Vrdn+w1zx58uQwxmr5RbUN+/v7wz7R/DL0zi5EIkjsQiSCxC5EIkjsQiSCxC5EIkjsQiTCULZ/mgPgtwBaUNruaY27P2tmzQDWA5iP0hZQ97g73aZ11KhRoaXEEgWiRAdmr7FNJFkCSp4Eibw2Dku4YMk6LOEimt+oJhwA9Pb2hrFNmzblGseSJUvKtjPbkF0XZjWxJJ8oESlvvbiGhoYw1t3dHcamT58exqJrzWy+6N6nSV5h5E+cBvArd78awM0AfmlmVwN4FMCb7r4YwJvZz0KIC5SKYnf3bnffkj0+CmAngFkAVgFYlz1tHYC7ajRGIUQVGNbf7GY2H8B1ADYBaHH3c59belD6mC+EuEAZstjNrAHAKwAecffzvoPope/7lf3On5mtNrN2M2tn9cmFELVlSGI3szEoCf1Fd381a+41s9Ys3gqgr1xfd1/j7m3u3tbc3FyNMQshclBR7FZaan4BwE53f3pQaAOAB7PHDwJ4vfrDE0JUi6Fkvf0QwAMAtpvZR1nbYwCeBPCymT0EYDeAeyodyMxCa4jZHVEmGssMGzduXBhjNcuiDCQGs97YOKKMLIBnAbJjRrz33nth7J133gljzLKbOnVqGFuxYkXZ9uuvvz7sw15zT09PGGNEFiy739ifm+zTKbPKjh8/HsYimE0ZWWzsXqwodnd/F0B0hB9X6i+EuDDQN+iESASJXYhEkNiFSASJXYhEkNiFSIQLpuAkswwiWEYZy6BiWzIx+yeyAJlVw2wyNo6Wlvjbx8w22rx5c9n2999/P+zz2WefhbGJEyeGsVtvvTWMLV26tGw7y/5iRRlZZhu7D6Lxs8KiLMOOZfoxK5JlU0b3XJ7sO5bNp3d2IRJBYhciESR2IRJBYhciESR2IRJBYhciEQq33iILgllvUSYaK66XF3bMyNZgFhqzcdi5mJ3ELKrIcmT7obFxMOtw0aJFYWzhwoVl29l+aMeOHQtjDGY3RRZbHqsX4HPF7DVGNBZ2rkgTLGtT7+xCJILELkQiSOxCJILELkQiSOxCJEKhq/FAvMLI6slFK8wsaYXFWA06tqIarZqysbNVX5ZAw1amDxw4EMb2799ftp2t0l5xxRVhbObMmWGstbU1jEWr4MwVYK+ZuRrMMZgwYULZdnYPsDHm3TaKOTbRPcfq3bGtzyL0zi5EIkjsQiSCxC5EIkjsQiSCxC5EIkjsQiRCRevNzOYA+C1KWzI7gDXu/qyZPQHg5wDOeT2Pufsb7FhnzpwJtzyaNGlS2C+ySVhSBatZxvoxyy6y0U6ePBn26e/vD2OsnhmzVvbs2RPGtmzZUradJVXMmjUrjK1atSqMLViwIIxFWyFFtdMqxTo7O8MYq/0W2WjMemP2Gruec+fODWOffvppGItsUTYfUaIUraEYRv7EaQC/cvctZtYIYLOZbcxiz7j7Pw3hGEKIOjOUvd66AXRnj4+a2U4A8VuBEOKCZFh/s5vZfADXAdiUNT1sZtvMbK2Zld8uUwhxQTBksZtZA4BXADzi7kcAPAdgIYBrUXrnfyrot9rM2s2snRVdEELUliGJ3czGoCT0F939VQBw9153P+PuZwE8D2BZub7uvsbd29y9jS2MCSFqS0WxW2kJ+gUAO9396UHtg7Mg7gbwcfWHJ4SoFkNZjf8hgAcAbDezj7K2xwDcZ2bXomTHdQL4xVBOGNlXLKvp+PHjZdujDC92HoBbXmxboDx9IgsK4HYYs5reeuutMBZZjitWrAj73H777WGMZZQx6+2rr74q297b2xv2iTLUAJ6Zt3v37jAWWVFTpkwJ+zQ2NoYxds1Y1h6r1xcdk93fkT3ILMWhrMa/C6CccqinLoS4sNA36IRIBIldiESQ2IVIBIldiESQ2IVIhEILTo4ePTq0V1hBxCjDhxXxY5lorB+LRYUvmfXGihD29PSEsQ8++CCMdXd3h7Eoe3DixIlhH2Z5saKSHR0dYSyyoVh2IyuwyK7nnDlzwlhka7HMtsjqrdSP2b2HDx8OY5EVzOzByGKjRVjDiBDie4XELkQiSOxCJILELkQiSOxCJILELkQiFGq9uXu4rxWzqCLLK2oHuA3CMoNYLI/1xjLs2J5izF778ssvw1iUicaKF06bNi2MMVimYmR5MYuV3QPM8mJ75kXXhu3Px8bBrK2omGql80WvjdmNkT3I9irUO7sQiSCxC5EIErsQiSCxC5EIErsQiSCxC5EIhVpvZhZmlTG7I7JxmK3FrCZWBJLBzhfBrDyWCbV3794wxrLllixZUra9qSnew4NlveWt9Z+nbDgr5sjuD7YvXpRlxyw0dp2Z3cjmimUPRq+b3TvRPcysTb2zC5EIErsQiSCxC5EIErsQiSCxC5EIFVfjzWwcgLcBXJo9//fu/riZXQ7gJQBTAGwG8IC7V1yujlYeWTJJlCjAtiZiK7tjx44NY6wGHVsdjWAru2xFmM0Hi0XJOuxcbB7ZKjJzBaIYm3tWn44lPbHXFl0zljDCxsjuKzb+PE4DS56JXAGWADaUd/ZTAH7k7j9AaXvmlWZ2M4BfA3jG3RcBGADw0BCOJYSoExXF7iXOGZljsn8O4EcAfp+1rwNwVy0GKISoDkPdn310toNrH4CNAL4AcMjdz31G2gtgVk1GKISoCkMSu7ufcfdrAcwGsAzAlUM9gZmtNrN2M2tnRQaEELVlWKvx7n4IwB8B3AJgspmdWymaDaAr6LPG3dvcvY1tAiCEqC0VxW5m08xscvZ4PICfANiJkuj/MnvagwBer9EYhRBVYCiJMK0A1pnZaJR+Obzs7v9pZjsAvGRm/wDgfwG8MJKB5KknxyyjvPXpWK2zyPJiVli0dRXAE1BY4sqiRYuGHWtsbAz7MEuRbUG0c+fOMBbVT8s7H8y6mj9/fhiLLDZ2f7AYs2aZVTYwMDDs8+XZ1orauWEkw923AbiuTHsHSn+/CyEuAvQNOiESQWIXIhEkdiESQWIXIhEkdiESwVjNqqqfzGw/gN3Zj1MBHCjs5DEax/loHOdzsY1jnruX3c+rULGfd2Kzdndvq8vJNQ6NI8Fx6GO8EIkgsQuRCPUU+5o6nnswGsf5aBzn870ZR93+ZhdCFIs+xguRCHURu5mtNLPPzGyXmT1ajzFk4+g0s+1m9pGZtRd43rVm1mdmHw9qazazjWb2efZ/nPZW23E8YWZd2Zx8ZGZ3FDCOOWb2RzPbYWafmNlfZ+2FzgkZR6FzYmbjzOwDM9uajePvs/bLzWxTppv1ZhZXxiyHuxf6D8BolMpaLQAwFsBWAFcXPY5sLJ0AptbhvLcBuB7Ax4Pa/hHAo9njRwH8uk7jeALA3xQ8H60Ars8eNwL4PwBXFz0nZByFzgkAA9CQPR4DYBOAmwG8DODerP1fAPzVcI5bj3f2ZQB2uXuHl0pPvwRgVR3GUTfc/W0A363RtQqlwp1AQQU8g3EUjrt3u/uW7PFRlIqjzELBc0LGUSheoupFXush9lkA9gz6uZ7FKh3AH8xss5mtrtMYztHi7t3Z4x4ALXUcy8Nmti37mF/zPycGY2bzUaqfsAl1nJPvjAMoeE5qUeQ19QW65e5+PYC/APBLM7ut3gMCSr/ZUfpFVA+eA7AQpT0CugE8VdSJzawBwCsAHnH3I4NjRc5JmXEUPic+giKvEfUQexeAwfV2wmKVtcbdu7L/+wC8hvpW3uk1s1YAyP7vq8cg3L03u9HOAngeBc2JmY1BSWAvuvurWXPhc1JuHPWak+zchzDMIq8R9RD7hwAWZyuLYwHcC2BD0YMwswlm1njuMYCfAviY96opG1Aq3AnUsYDnOXFl3I0C5sRKRQFfALDT3Z8eFCp0TqJxFD0nNSvyWtQK43dWG+9AaaXzCwB/W6cxLEDJCdgK4JMixwHgdyh9HPwWpb+9HkJpz7w3AXwO4H8ANNdpHP8GYDuAbSiJrbWAcSxH6SP6NgAfZf/uKHpOyDgKnRMAf4ZSEddtKP1i+btB9+wHAHYB+A8Alw7nuPoGnRCJkPoCnRDJILELkQgSuxCJILELkQgSuxCJILELkQgSuxCJILELkQj/D5g+OwpgAVglAAAAAElFTkSuQmCC",
      "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",
    "# 创建 VQSVD 网络并训练\n",
    "net = VQSVD(matrix=Mat_generator(), weights=weight, num_qubits=num_qubits, depth=cir_depth, rank=RANK, lr=LR, itr=ITR, seed=SEED)\n",
    "loss_list, singular_value_list = net.train()\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = net.get_matrix_U().numpy()\n",
    "V_dagger_learned = dagger(net.get_matrix_V()).numpy()\n",
    "\n",
    "singular_value = singular_value_list[-1]\n",
    "mat = np.matrix(U_learned.real[:, :RANK]) * np.diag(singular_value[:RANK])* np.matrix(V_dagger_learned.real[:RANK, :])\n",
    "\n",
    "reconstimg = mat\n",
    "plt.imshow(reconstimg, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Wang, X., Song, Z., & Wang, Y. Variational Quantum Singular Value Decomposition. [Quantum, 5, 483 (2021).](https://quantum-journal.org/papers/q-2021-06-29-483/)"
   ]
  }
 ],
 "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.8.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}