tutorial/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": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.433321Z",
     "start_time": "2021-03-09T03:47:28.730090Z"
    }
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "from numpy import pi as PI\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import unitary_group\n",
    "from scipy.linalg import norm\n",
    "\n",
    "import paddle\n",
    "from paddle import matmul, transpose, trace\n",
    "from paddle_quantum.circuit import *\n",
    "from paddle_quantum.utils import *\n",
    "\n",
    "\n",
    "# 画出优化过程中的学习曲线\n",
    "def loss_plot(loss):\n",
    "    '''\n",
    "    loss is a list, this function plots loss over iteration\n",
    "    '''\n",
    "    plt.plot(list(range(1, len(loss)+1)), loss)\n",
    "    plt.xlabel('iteration')\n",
    "    plt.ylabel('loss')\n",
    "    plt.title('Loss Over Iteration')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 经典奇异值分解\n",
    "\n",
    "那么在了解一些简单的数学背景之后， 我们来学习下如何用 Numpy 完成矩阵的奇异值分解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.469286Z",
     "start_time": "2021-03-09T03:47:34.440399Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "我们想要分解的矩阵 M 是：\n",
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 生成矩阵 M\n",
    "def M_generator():\n",
    "    I = np.array([[1, 0], [0, 1]])\n",
    "    Z = np.array([[1, 0], [0, -1]])\n",
    "    X = np.array([[0, 1], [1, 0]])\n",
    "    Y = np.array([[0, -1j], [1j, 0]])\n",
    "    M = 2 *np.kron(X, Z) + 6 * np.kron(Z, X) + 3 * np.kron(I, I)\n",
    "    return M.astype('complex64')\n",
    "\n",
    "print('我们想要分解的矩阵 M 是：')\n",
    "print(M_generator())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.541244Z",
     "start_time": "2021-03-09T03:47:34.489833Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵的奇异值从大到小分别是:\n",
      "[11.  7.  5.  1.]\n",
      "分解出的酉矩阵 U 是:\n",
      "[[-0.5+0.j -0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j -0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [ 0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]]\n",
      "分解出的酉矩阵 V_dagger 是:\n",
      "[[-0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j  0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j -0.5+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 我们只需要以下一行代码就可以完成 SVD \n",
    "U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)\n",
    "\n",
    "# 打印分解结果\n",
    "print(\"矩阵的奇异值从大到小分别是:\")\n",
    "print(D)\n",
    "print(\"分解出的酉矩阵 U 是:\")\n",
    "print(U)\n",
    "print(\"分解出的酉矩阵 V_dagger 是:\")\n",
    "print(V_dagger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.600873Z",
     "start_time": "2021-03-09T03:47:34.565570Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 再组装回去，能不能复原矩阵？\n",
    "M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))\n",
    "print(M_reconst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "那当然是可以复原成原来的矩阵 $M$ 的！读者也可以自行修改矩阵，试试看不是方阵的情况。\n",
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子奇异值分解\n",
    "\n",
    "接下来我们来看看量子版本的奇异值分解是怎么一回事。简单的说，我们把矩阵分解这一问题巧妙的转换成了优化问题。通过以下四个步骤：\n",
    "\n",
    "- 准备一组正交向量基 $\\{|\\psi_j\\rangle\\}$, 不妨直接取计算基 $\\{ |000\\rangle, |001\\rangle,\\cdots |111\\rangle\\}$ （这是3量子比特的情形）\n",
    "- 准备两个参数化的量子神经网络 $U(\\theta)$ 和 $V(\\phi)$ 分别用来学习左/右奇异向量\n",
    "- 利用量子神经网络估算奇异值 $m_j = \\text{Re}\\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle$\n",
    "- 设计损失函数并且利用飞桨来优化\n",
    "\n",
    "$$\n",
    "L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle, \\tag{3}\n",
    "$$\n",
    "\n",
    "其中 $q_1>\\cdots>q_T>0$ 是可以调节的权重（超参数）, $T$ 表示我们想要学习到的阶数（rank）或者可以解释为总共要学习得到的奇异值个数。\n",
    "\n",
    "\n",
    "\n",
    "### 案例1：分解随机生成的 $8\\times8$ 复数矩阵\n",
    "\n",
    "接着我们来看一个具体的例子，这可以更好的解释整体流程。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.656799Z",
     "start_time": "2021-03-09T03:47:34.620380Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "我们想要分解的矩阵 M 是：\n",
      "[[6.+1.j 3.+9.j 7.+3.j 4.+7.j 6.+6.j 9.+8.j 2.+7.j 6.+4.j]\n",
      " [7.+1.j 4.+4.j 3.+7.j 7.+9.j 7.+8.j 2.+8.j 5.+0.j 4.+8.j]\n",
      " [1.+6.j 7.+8.j 5.+7.j 1.+0.j 4.+7.j 0.+7.j 9.+2.j 5.+0.j]\n",
      " [8.+7.j 0.+2.j 9.+2.j 2.+0.j 6.+4.j 3.+9.j 8.+6.j 2.+9.j]\n",
      " [4.+8.j 2.+6.j 6.+8.j 4.+7.j 8.+1.j 6.+0.j 1.+6.j 3.+6.j]\n",
      " [8.+7.j 1.+4.j 9.+2.j 8.+7.j 9.+5.j 4.+2.j 1.+0.j 3.+2.j]\n",
      " [6.+4.j 7.+2.j 2.+0.j 0.+4.j 3.+9.j 1.+6.j 7.+6.j 3.+8.j]\n",
      " [1.+9.j 5.+9.j 5.+2.j 9.+6.j 3.+0.j 5.+3.j 1.+3.j 9.+4.j]]\n",
      "矩阵的奇异值从大到小分别是:\n",
      "[54.83484985 19.18141073 14.98866247 11.61419557 10.15927045  7.60223249\n",
      "  5.81040539  3.30116001]\n"
     ]
    }
   ],
   "source": [
    "# 先固定随机种子， 为了能够复现结果\n",
    "np.random.seed(42)\n",
    "\n",
    "# 设置量子比特数量，确定希尔伯特空间的维度\n",
    "N = 3\n",
    "\n",
    "# 制作随机矩阵生成器\n",
    "def random_M_generator():\n",
    "    M = np.random.randint(10, size = (2**N, 2**N)) + 1j*np.random.randint(10, size = (2**N, 2**N))\n",
    "    M1 = np.random.randint(10, size = (2**N, 2**N)) \n",
    "    return M\n",
    "\n",
    "M = random_M_generator()\n",
    "M_err = np.copy(M)\n",
    "\n",
    "# 打印结果\n",
    "print('我们想要分解的矩阵 M 是：')\n",
    "print(M)\n",
    "\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "print(\"矩阵的奇异值从大到小分别是:\")\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.692350Z",
     "start_time": "2021-03-09T03:47:34.671114Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "选取的等差权重为：\n",
      "[24.+0.j 21.+0.j 18.+0.j 15.+0.j 12.+0.j  9.+0.j  6.+0.j  3.+0.j]\n"
     ]
    }
   ],
   "source": [
    "# 超参数设置\n",
    "N = 3           # 量子比特数量\n",
    "T = 8           # 设置想要学习的阶数\n",
    "ITR = 100       # 迭代次数\n",
    "LR = 0.02       # 学习速率\n",
    "SEED = 14       # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
    "print('选取的等差权重为：')\n",
    "print(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们搭建如下的量子神经网络结构："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.725007Z",
     "start_time": "2021-03-09T03:47:34.702861Z"
    }
   },
   "outputs": [],
   "source": [
    "# 设置电路参数\n",
    "cir_depth = 20                           # 电路深度\n",
    "block_len = 2                            # 每个模组的长度\n",
    "theta_size = N * block_len * cir_depth   # 网络参数 theta 的大小\n",
    "\n",
    "\n",
    "# 定义量子神经网络\n",
    "def U_theta(theta):\n",
    "\n",
    "    # 用 UAnsatz 初始化网络\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # 搭建层级结构：\n",
    "    for layer_num in range(cir_depth):\n",
    "        \n",
    "        for which_qubit in range(N):\n",
    "            cir.ry(theta[block_len * layer_num * N + which_qubit], \n",
    "                                              which_qubit)\n",
    "        \n",
    "        for which_qubit in range(N):\n",
    "            cir.rz(theta[(block_len * layer_num + 1) * N \n",
    "                              + which_qubit], which_qubit)\n",
    "\n",
    "        for which_qubit in range(1, N):\n",
    "            cir.cnot([which_qubit - 1, which_qubit])\n",
    "        cir.cnot([N - 1, 0])\n",
    "\n",
    "    return cir.U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来完成算法的主体部分：\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.011232Z",
     "start_time": "2021-03-09T03:47:40.712257Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: 259.1903\n",
      "iter: 10 loss: -1672.2268\n",
      "iter: 20 loss: -2097.8083\n",
      "iter: 30 loss: -2242.9914\n",
      "iter: 40 loss: -2310.1948\n",
      "iter: 50 loss: -2340.2588\n",
      "iter: 60 loss: -2357.5961\n",
      "iter: 70 loss: -2369.2625\n",
      "iter: 80 loss: -2376.9403\n",
      "iter: 90 loss: -2373.7986\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAodklEQVR4nO3de3xcdZ3/8ddnZnK/NG2T3i8ptEVbpBTDXZAfoBQFqy4isApedllc8bK6u4K4v119yO+Bq7teFldFZUVXQRZEWe6giIAWSAFLS62EttD0mrRpmzb3zOf3x/mmnYa0JJlMJpl5Px+Pecyc7zkz53M6MO+c7znne8zdERERSUcs2wWIiMj4pzAREZG0KUxERCRtChMREUmbwkRERNKmMBERkbQpTERkSMxsn5kdle06ZGxRmMi4Y2YbzezcLK37NDP7jZm1mtkeM/tfM1s0ius/sO1m9iEzeyLD6/utmf1Vapu7l7v7+kyuV8YfhYnIIJnZqcBDwK+AGcA84I/AkyP9l7pFMvr/p5klMvn5kl8UJpIzzKzIzL5hZlvC4xtmVhTmVZvZPWa228x2mdnjfT/WZvY5M9sc9jbWmdk5h1nFvwI/dvdvunuru+9y9y8AK4B/CZ+11swuSKkpYWZNZnZCmD7FzH4f6vijmZ2Vsuxvzex6M3sSaAMOG1Bm9kbgu8Cpodtpd8q/wdfM7FUz225m3zWzkjDvLDNrDNu7DfgvM5sY/l2azKwlvJ4Vlr8eOAO4MazjxtDuZjY/vJ5gZj8O73/FzL6Q8u/6ITN7ItTTYmYbzOz8QX+hMq4oTCSXXAecAhwPLAFOAr4Q5n0WaARqgKnA5wE3s2OAq4ET3b0COA/Y2P+DzawUOA34nwHWezvwtvD6VuDSlHnnAc3u/qyZzQTuBb4MTAL+HrjTzGpSlv8gcCVQAbxyuA1197XAVcAfQrdTVZh1A7Aw/BvMB2YC/zflrdPCuueG9cSA/wrTc4B24MawjuuAx4GrwzquHqCU/wAmEAXfW4HLgQ+nzD8ZWAdUE4XxD83MDrddMn4pTCSX/CXwJXff4e5NwBeJfpwBuoHpwFx373b3xz0amK4XKAIWmVmBu29095cH+OxJRP+/bB1g3laiH0uAnwHvCuEDcBlRwAB8ALjP3e9z96S7PwzUA+9I+awfufsad+9x9+6hbHz4kb4S+Luw19QK/D/gkpTFksA/u3unu7e7+053v9Pd28Ly1xOFwmDWFw+ffW3YU9sI/BsH/80BXnH377t7L3AL0XcwdSjbJeODwkRyyQwO/Wv+ldAG8FWgAXjIzNab2TUA7t4AfJqom2qHmd1mZjN4rRaiH+LpA8ybDjSnfN5a4MIQKO8iChiI/vp/X+ji2h26pt7S7zM3DWWD+6kBSoGVKZ//QGjv0+TuHX0TZlZqZt8LXVR7gd8BVSEoXk81UMBr/81npkxv63vh7m3hZfkQtknGCYWJ5JItRD/YfeaENsJfzp9196OIfuA/03dsxN1/5u5vCe914Cv9P9jd9wN/AN43wHovBn6dMt3X1bUceDEEDERB8RN3r0p5lLn7DamrGsL29l+2maibanHK509w9/IjvOezwDHAye5eCZwZ2u0wy/dfXzev/TffPIRtkByhMJHxqsDMilMeCaIf8S+YWY2ZVRMdK/hvADO7wMzmh66gPUTdW0kzO8bMzg4H6juIfoyTh1nnNcAVZvZJM6sIB6+/DJxK1KXW5zbg7cDHOLhXQqjlQjM7z8zioe6z+g54D8N2YJaZFQK4exL4PvB1M5sStnummZ13hM+oINrm3WY2CfjnAdYx4IkAoevqduD68O8xF/hM2E7JMwoTGa/uI/oR7Hv8C9GB7XpgFfAC8GxoA1gAPALsI9rD+E93f5ToeMkNRH9lbwOmANcOtEJ3f4LogPp7iY6TvAIsBd7i7i+lLLc1rOM04Ocp7ZuI9lY+DzQR7an8A8P///A3wBpgm5k1h7bPEXXnrQjdVo8Q7XkczjeAEqLtX0HULZbqm8BF4Wysbw3w/k8A+4H1wBNE4XnzsLZGxjXTzbFERCRd2jMREZG0KUxERCRtChMREUmbwkRERNKWtwO9VVdXe21tbbbLEBEZV1auXNns7jX92/M2TGpra6mvr892GSIi44qZDThmnLq5REQkbQoTERFJm8JERETSpjAREZG0KUxERCRtChMREUmbwkRERNKmMBmiu55r5L9XHPbW3CIieUlhMkT3vbBNYSIi0o/CZIhqKopoau3MdhkiImOKwmSIasqL2NXWRXfv4e7sKiKSfxQmQ1RTUYQ77Nrfle1SRETGDIXJENVUFAGoq0tEJIXCZIgUJiIir6UwGaKacoWJiEh/CpMhOrBnsk9hIiLSR2EyRMUFcSqLE9ozERFJoTAZBl1rIiJyKIXJMChMREQOpTAZhpqKYna0dmS7DBGRMUNhMgw15dozERFJpTAZhpqKIvZ39bK/syfbpYiIjAkKk2HoOz24WacHi4gACpNh0VXwIiKHypkwMbNlZrbOzBrM7JpMrktXwYuIHConwsTM4sC3gfOBRcClZrYoU+ubUqmr4EVEUuVEmAAnAQ3uvt7du4DbgOWZWtnE0kLiMdOeiYhIkCthMhPYlDLdGNoyIh4zJpcVKkxERIJcCZNBMbMrzazezOqbmprS+ixdBS8iclCuhMlmYHbK9KzQdgh3v8nd69y9rqamJq0V1lQU6ZiJiEiQK2HyDLDAzOaZWSFwCXB3JldYU17Ejr0KExERgES2CxgJ7t5jZlcDDwJx4GZ3X5PJddZUFNG8r5Nk0onFLJOrEhEZ83IiTADc/T7gvtFaX01FET1JZ3d7N5PKCkdrtSIiY1KudHONOl0FLyJykMJkmHQVvIjIQQqTYTp4L3jd10RERGEyTFMqiwHtmYiIgMJk2MoK45QUxBUmIiIoTIbNzHQVvIhIoDBJg66CFxGJKEzSoKvgRUQiCpM0TKksYvtenc0lIqIwScPUymL2dvTQ3tWb7VJERLJKYZKGaeH04G3aOxGRPKcwScO0CSFM9ihMRCS/KUzSMDXsmei4iYjkO4VJGg7smShMRCTPKUzSUF6UoLwooT0TEcl7CpM0TdXpwSIiCpN0Ta0s1gF4Ecl7CpM0TassZruugheRPKcwSdPUCcVs39tBMunZLkVEJGsUJmmaVllMT9LZub8r26WIiGSNwiRNutZERERhkjZdBS8iojBJm8bnEhFRmKSturyQmKmbS0Tym8IkTYl4jJqKInVziUheU5iMgGmVxermEpG8pjAZAVMri9XNJSJ5TWEyAqZN0JAqIpLfFCYjQLfvFZF8N+bCxMz+xcw2m9nz4fGOlHnXmlmDma0zs/NS2peFtgYzu2a0a9aFiyKS7xLZLuAwvu7uX0ttMLNFwCXAYmAG8IiZLQyzvw28DWgEnjGzu939xdEqNvVak9rqstFarYjImDFWw2Qgy4Hb3L0T2GBmDcBJYV6Du68HMLPbwrKjFyYTigDtmYhI/hpz3VzB1Wa2ysxuNrOJoW0msCllmcbQdrj21zCzK82s3szqm5qaRqzYvm4uHYQXkXyVlTAxs0fMbPUAj+XAd4CjgeOBrcC/jdR63f0md69z97qampqR+lgqigsoK4zrWhMRyVtZ6eZy93MHs5yZfR+4J0xuBmanzJ4V2jhC+6jpu6+JiEg+GnPdXGY2PWXyPcDq8Ppu4BIzKzKzecAC4GngGWCBmc0zs0Kig/R3j2bNEK6CVzeXiOSpsXgA/l/N7HjAgY3A3wC4+xozu53owHoP8HF37wUws6uBB4E4cLO7rxntoqdVFvPUhl2jvVoRkTFhzIWJu3/wCPOuB64foP0+4L5M1vV6ZlSVsG1vB929SQriY26HT0Qko/SrN0LmTi6lN+k0trRnuxQRkVGnMBkh88LFihub92e5EhGR0acwGSF9V75v3KkwEZH8ozAZIZPLCqkoSmjPRETyksJkhJgZc6tL2bCzLduliIiMOoXJCKqdXKY9ExHJSwqTETSvuozGlja6epLZLkVEZFQpTEZQ7eQykg6NLerqEpH8ojAZQTqjS0TylcJkBNVOLgVgQ7P2TEQkvyhMRtCkskIqinV6sIjkH4XJCDIz5lWXqZtLRPKOwmSE1U5WmIhI/lGYjLDayaVsbmnX6cEiklcUJiOstjo6PfjVXToILyL5Q2Eywmo1erCI5CGFyQibN1nXmohI/lGYjLCJZYVMKClQmIhIXlGYZEDt5FI26sJFEckjCpMMqK0uY4OOmYhIHlGYZMC86jK27Gmnrasn26WIiIwKhUkGHDdrAu6wqnFPtksRERkVCpMMWDp7IgDPvtqS5UpEREaHwiQDJpYVclR1Gc++sjvbpYiIjAqFSYYsnTOR515twd2zXYqISMYpTDJk6Zwqdu7vYtOu9myXIiKScQqTDDlhjo6biEj+UJhkyDHTKigrjCtMRCQvZCVMzOx9ZrbGzJJmVtdv3rVm1mBm68zsvJT2ZaGtwcyuSWmfZ2ZPhfafm1nhaG7L4cRjxpLZVQoTEckL2dozWQ28F/hdaqOZLQIuARYDy4D/NLO4mcWBbwPnA4uAS8OyAF8Bvu7u84EW4KOjswmv74Q5E1m7tZX2rt5slyIiklGDChMz+5SZVVrkh2b2rJm9fbgrdfe17r5ugFnLgdvcvdPdNwANwEnh0eDu6929C7gNWG5mBpwN3BHefwvw7uHWNdKWzqmiN+msatyd7VJERDJqsHsmH3H3vcDbgYnAB4EbMlDPTGBTynRjaDtc+2Rgt7v39GsfkJldaWb1Zlbf1NQ0ooUPZOmBg/C7M74uEZFsGmyYWHh+B/ATd1+T0jbwG8weMbPVAzyWp1NwOtz9Jnevc/e6mpqajK9vUlkh86rLdNxERHJeYpDLrTSzh4B5wLVmVgEc8Sbn7n7uMOrZDMxOmZ4V2jhM+06gyswSYe8kdfkxYemcKn7352bcnahXTkQk9wx2z+SjwDXAie7eBhQAH85APXcDl5hZkZnNAxYATwPPAAvCmVuFRAfp7/bo8vJHgYvC+68AfpWBuoatbu4kmvd18tKOfdkuRUQkYwYbJqcC69x9t5l9APgCMOwhcc3sPWbWGD73XjN7ECB0n90OvAg8AHzc3XvDXsfVwIPAWuD2sCzA54DPmFkD0TGUHw63rkx426KpxAzu+eOWbJciIpIxNpixo8xsFbAEOA74EfAD4GJ3f2tGq8uguro6r6+vH5V1Xfb9FWzb08GvP/tWdXWJyLhmZivdva5/+2D3THpCl9Jy4EZ3/zZQMZIF5rILl8xgffN+1mzZm+1SREQyYrBh0mpm1xKdEnyvmcWIjpvIICxbPI1EzLhn1dZslyIikhGDDZP3A51E15tsIzpr6qsZqyrHTCwr5PT51dyzaouGpBeRnDSoMAkB8lNggpldAHS4+48zWlmOuXDJDBpb2nl+0+5slyIiMuIGO5zKxUSn6L4PuBh4yswuOvK7JNXbF0+lMB5TV5eI5KTBdnNdR3SNyRXufjnRWFn/lLmyck9lcQFnLqzh3lVbSSbV1SUiuWWwYRJz9x0p0zuH8F4JLlwynW17O/jD+p3ZLkVEZEQNNhAeMLMHzexDZvYh4F7gvsyVlZvOWzyNCSUF/OzpV7NdiojIiBrU2Fzu/g9m9hfA6aHpJne/K3Nl5abigjh/ccIsfrJiI837OqkuL8p2SSIiI2LQXVXufqe7fyY8FCTDdNnJs+nude5Y2ZjtUkRERswRw8TMWs1s7wCPVjPT5dzDMH9KBSfNm8StT7+qA/EikjOOGCbuXuHulQM8Kty9crSKzDWXnTSHV3a26UC8iOQMnZGVBcuOncbE0gJ+9pQOxItIblCYZEHfgfgH12yjqbUz2+WIiKRNYZIll548h56k8z8rN73+wiIiY5zCJEuOrinn5HmTuO3pTToQLyLjnsIkiy47eQ6v7mrjyZebs12KiEhaFCZZtOzYaUwqK9SBeBEZ9xQmWVSUiHPRm2fx8Ivb2dHake1yRESGTWGSZZecODs6EF+vK+JFZPxSmGTZUTXlnHrUZG57RlfEi8j4pTAZAy47eQ6bdrXzeIMOxIvI+KQwGQPevngqE0sLuL1e15yIyPikMBkDihJx3r10Jg+v2U7L/q5slyMiMmQKkzHi/SfOpqs3yV3Pbc52KSIiQ6YwGSPeMK2SJbMmcHv9Jtx1IF5ExheFyRhy8Ymz+dO2Vl7YvCfbpYiIDInCZAy5cMkMigti/PwZHYgXkfFFYTKGVBYX8I5jp3P381to7+rNdjkiIoOWlTAxs/eZ2RozS5pZXUp7rZm1m9nz4fHdlHlvNrMXzKzBzL5lZhbaJ5nZw2b2UniemI1tGikXnzib1s4e7l+9NduliIgMWrb2TFYD7wV+N8C8l939+PC4KqX9O8BfAwvCY1lovwb4tbsvAH4dpsetk+dNonZyKbc9ra4uERk/shIm7r7W3dcNdnkzmw5UuvsKj051+jHw7jB7OXBLeH1LSvu4ZGZcctIcnt64i4YdrdkuR0RkUMbiMZN5ZvacmT1mZmeEtplA6kiIjaENYKq79/UJbQOmHu6DzexKM6s3s/qmpqYRL3ykXPTmWRTETXsnIjJuZCxMzOwRM1s9wGP5Ed62FZjj7kuBzwA/M7PKwa4z7LUc9iINd7/J3evcva6mpmbQ2zLaqsuLeNuiqdz5bCOdPToQLyJjXyJTH+zu5w7jPZ1AZ3i90sxeBhYCm4FZKYvOCm0A281surtvDd1hO9KrfGy49KQ53PfCNh5cs513LZmR7XJERI5oTHVzmVmNmcXD66OIDrSvD91Ye83slHAW1+XAr8Lb7gauCK+vSGkf104/uprZk0q4VXdhFJFxIFunBr/HzBqBU4F7zezBMOtMYJWZPQ/cAVzl7rvCvL8FfgA0AC8D94f2G4C3mdlLwLlhetyLxYxLTpzDH9bvZEPz/myXIyJyRJav40DV1dV5fX19tss4oh17Ozj1ht/wkdNrue6di7JdjogIZrbS3ev6t4+pbi451JTKYpYtnsbt9Y26Il5ExjSFyRh3xWm17Gnv5pfPa2h6ERm7FCZj3Im1E3nj9Ep+9ORGDU0vImOWwmSMMzM+dNpc1m1vZcX6Xa//BhGRLFCYjAPLj59JVWkBt/x+Y7ZLEREZkMJkHCguiPP+E2fz0Ivb2Ly7PdvliIi8hsJknPjgKXMBtHciImOSwmScmDWxlAuXzOAnf3iFnfs6s12OiMghFCbjyCfOXkBHTy/ff3xDtksRETmEwmQcmT+lnAuPm8GP/7CRXfu7sl2OiMgBCpNx5pPnzKe9u5fvP74+26WIiBygMBln5k+p4ILjZnDL77V3IiJjh8JkHPrk2dHeyfceeznbpYiIAAqTcWnB1Areu3QWNz+5gfVN+7JdjoiIwmS8+tz5x1CUiPOle17UmF0iknUKk3FqSkUxnz53Ab9d18Sv1+bEnYpFZBxTmIxjV5xWy4Ip5XzxnjV0dOt+JyKSPQqTcawgHuOL71rMpl3tfOe3OhgvItmjMBnnTptfzfLjZ3Djow0892pLtssRkTylMMkBX1p+LNMqi/nkbc/R2tGd7XJEJA8pTHLAhJICvnnJ8Wxuaeeffrk62+WISB5SmOSIutpJfOqchfzy+S384tnGbJcjInlGYZJDrj57PifNm8R1d63mxS17s12OiOQRhUkOiceMGy9byoSSAq78ST0tGrtLREaJwiTHTKko5nsffDM7Wju5+tZn6elNZrskEckDCpMctGR2Fde/+1iebNjJl+9dq+FWRCTjEtkuQDLjfXWzWbu1lZuf3EB1eSFXn70g2yWJSA5TmOSwL7zzjexu6+JrD/2ZsqIEHz59XrZLEpEclZVuLjP7qpn9ycxWmdldZlaVMu9aM2sws3Vmdl5K+7LQ1mBm16S0zzOzp0L7z82scJQ3Z8yKxYx/veg4zls8lS/+74vc/symbJckIjkqW8dMHgaOdffjgD8D1wKY2SLgEmAxsAz4TzOLm1kc+DZwPrAIuDQsC/AV4OvuPh9oAT46qlsyxiXiMb516VLOWFDN536xip899Wq2SxKRHJSVMHH3h9y9J0yuAGaF18uB29y90903AA3ASeHR4O7r3b0LuA1YbmYGnA3cEd5/C/DuUdqMcaMoEef7l9dx1sIaPn/XC/xA948XkRE2Fs7m+ghwf3g9E0jti2kMbYdrnwzsTgmmvvYBmdmVZlZvZvVNTU0jVP74UFwQ53sfrOP8Y6fx5XvX8o1H/qyzvERkxGQsTMzsETNbPcBjecoy1wE9wE8zVUcqd7/J3evcva6mpmY0VjmmFCZi/MelS3nvCTP5xiMv8ZEfPUPzvs5slyUiOSBjZ3O5+7lHmm9mHwIuAM7xg38ibwZmpyw2K7RxmPadQJWZJcLeSeryMoBEPMa/vW8Jx8+u4sv3ruX8bz7Ov1+8hDMW5F+4isjIydbZXMuAfwTe5e5tKbPuBi4xsyIzmwcsAJ4GngEWhDO3CokO0t8dQuhR4KLw/iuAX43WdoxXZsblp9byq4+fzoSSAj74w6f5xK3PsXl3e7ZLE5FxyrLRb25mDUAR0Z4FwAp3vyrMu47oOEoP8Gl3vz+0vwP4BhAHbnb360P7UUQH5CcBzwEfcPfX7bupq6vz+vr6kdyscam9q5fvPPYy33ssulPjX59xFJefOpcplcVZrkxExiIzW+nuda9pz9eDsAqTQzW2tPGVB9bxv3/cQjxmnLWwhotPnM1Zx9RQlIhnuzwRGSMUJv0oTAa2oXk/t9dv4o6VjTS1dlJRlODcRVN5x5umc8aCaooLFCwi+Uxh0o/C5Mh6epM83tDMfau28tCL29nT3k1ZYZz/84YpnH/sdM554xQFi0geUpj0ozAZvK6eJL9/uZkH12zjoTXb2bm/i8riBO9ZOpP3nziHRTMqs12iiIwShUk/CpPh6U06K9bv5Pb6Tdy/ehtdPUlOmjeJj731aM46poZoUAIRyVUKk34UJunb3dbFHSsbufmJDWzZ08EbplXwsbOO5p1vmk4iPhYGVxCRkaYw6UdhMnK6e5Pc/fwWvvvYy7y0Yx9zJ5dy1VuP5r0nzNSZYCI5RmHSj8Jk5CWTzkMvbufbjzbwwuY91FQUccWpc/nLk+cysUx3BhDJBQqTfhQmmePuPNHQzE2/W8/jLzVTUhDnguOmc8GSGZx29GQK1AUmMm4dLkx0p0UZcWbGGQtqOGNBDeu2tXLzExu474Wt/M/KRqpKCzj7DVM4/ehqTp9fzbQJutJeJBdoz0RGRUd3L4+/1My9q7bw2J+baGnrBuCo6jJOrJ1EXe1ETp43mTmTS7NcqYgcibq5+lGYZE8y6azdtpcnG5p5esMuntnYwp72KFxmTyrhLfOrOXNBDWcurKGsSDvPImOJwqQfhcnYkUw6DU37WLF+J0+81MwfXt5Ja2cPRYkYb11Yw/lvmsbZb5jKhJKCbJcqkvcUJv0oTMaunt4k9a+0cP8LW3lgzTa27+2kIG6cenQ15y2eypkLapg9Sd1hItmgMOlHYTI+JJPOc5t289Cabdy/ehuv7opufzOvuozT50/muJlVLJpRycKpFRQmdJaYSKYpTPpRmIw/7s5LO/bx+EvNPPFSE09t2EVbVy8AiZgxbUIxsyaWMLOqlBlVxUyfUML0vraJJZQW6viLSLp0arCMe2bGwqkVLJxawUffMo/epPPKzv2s2bKXP23bS2NLO5tb2vn9y81s39tBst/fSZPKCpk9qZTayaXMnVTK3Mll1FZHz5PLCjWumEgatGciOamnN8mO1k627mmnsaXv0caru9rY2NzG1j3th4RNSUGcGVXFzKgqYVplMVMqi6gpL2JKZTHV5UXUVBRRXV5IeVFCoSN5TXsmklcS8RgzqkqYUVXCm+e+dn5nTy+NLe28urONjTv3s7mlnc2729myu50/b2+leV8Xvf13bYCiRIzq8iImlxcyuayQSWXR60llhUwqLaSqtICq8DyhpICK4gQlBXEFkOQ8hYnkpaJEnKNryjm6pnzA+b1JZ9f+LppaO2ned/Cxc18XTfs6aQ7P67a10ry/i66e5GHXFY8ZFcUJKooTVBYXHHieUBI9qkoLmFBaGL1ObSspoLwooRGYZVxQmIgMIB4zaiqi7q3X4+60d/eya38XLfu72dPeze72Lna3ddPa0UNrx8HnvR097G3vZuPO/ext72FPezft3b1H/PyywjjlxQnKihKUFyUoLYxTWhjt8RQVxCguiFOciFMcXhclYpQUhrbCOCUFcUoL4xSH59LCOCWhvaQgrrCSEaEwEUmTmVFamKC0MMGsiUN/f2dPL3vau9nTFgXRnvZudrd1s7ej+0Dg7O/sYV94tHf1sqO1g7auXjq7k3T29NLRnaSju5eeAbrmXk9hPHYgiIpDwBQXxCgK08WJ8DoRo6ggRlEiCqyiRDxMxyhMxCiMR88Hp+MUJmIUxO3A/IKwTEE8TCeMgniMRMzGXVfgnrZu1m1vZd32VrbsbqeiOMHE0kKmVBRxylGT8270hvzaWpExqCgRZ0pFnCkV6Q962dObpKMnCpa+R1tXL+1dvbR1R88HX/fQ0Z2kPbR3dPfSfuB90WfsaetiR0+SzvCZnT1JOrt76ehJDnhMabjMOBgw8Shg+oInEeubNhLxg9OJuJGIHdqeiEWvC8K8aJmU+XGjIBYjHl7HzIjHjLgZZtEfBn2R1pt0et3p6O49sGfZvK+LV3bu55Wdbezc33Wg/phxyAkdRYkYZy6s4W2LpjJ9QjGlhQkK4saftrXyx027eXHrXhIxC12dhRxVU8aSWVW8adYEEjFjR2snTa2d7Nrfxe62Lna3d1NWGOfoKeUsmFJBdfnBsw87e3ppau1kR2snnd3JA3vUlcWvPVmkuzdJS1sXk0oLR3yPVGEikkMS8Rjl8Rjlo/BXcU9vkq7eJJ3d0XNXCJ2unoPT3b2vne4+8NoPtEef5dHr5MH39PT6gfd09zo9yei5rauHnqRHbb3J8DoKuL7lenud7mT0GcPZY+uvuCDGpNJC5k4u4+2Lp1I7uYyF0yo4ZmoF0ycU09Ed/VC/srONB9ds44HV23j4xe2v+ZyK4gSLZ1QSM2PL7g7WbNnLnc82DrkeMzB4zSnwfRIxO9A1mogbLfu72NvRA8Cjf38W86rLhrzOI1GYiMiwJOIxEvEYpePgvmfuTm8yCpWepB8ImmTY++hNOu5ED6Jf53gs2rspTMSoKE687n14SgrjlBRGZxCeevRk/u8Fi2ho2negm7KzJ8nRNeUcVV1GLHboHsPuti5e2LyH1Zv3YgY14XT0SWWFTCwrZGJpAXvbe2jYsY+XdrTSsr8LJ6q3KBGLTmWvKKIoEad5X7RXs3N/F22dPezr7KW7N8nE0gImlkVnIVZlYJw7XWciIiKDdrjrTHQah4iIpE1hIiIiaVOYiIhI2rISJmb2VTP7k5mtMrO7zKwqtNeaWbuZPR8e3015z5vN7AUzazCzb1k4583MJpnZw2b2Ungexpn+IiKSjmztmTwMHOvuxwF/Bq5Nmfeyux8fHleltH8H+GtgQXgsC+3XAL929wXAr8O0iIiMoqyEibs/5O49YXIFMOtIy5vZdKDS3Vd4dPrZj4F3h9nLgVvC61tS2kVEZJSMhWMmHwHuT5meZ2bPmdljZnZGaJsJpF7V0xjaAKa6+9bwehswNaPViojIa2TsokUzewSYNsCs69z9V2GZ64Ae4Kdh3lZgjrvvNLM3A780s8WDXae7u5kd9sIZM7sSuBJgzpw5g/1YERF5HRkLE3c/90jzzexDwAXAOaHrCnfvBDrD65Vm9jKwENjMoV1hs0IbwHYzm+7uW0N32I4j1HQTcFNYf5OZvTKETaoGmoewfC7Ix22G/NzufNxmyM/tTnebB7hDUJaGUzGzZcA/Am9197aU9hpgl7v3mtlRRAfa17v7LjPba2anAE8BlwP/Ed52N3AFcEN4/tVganD3miHWXD/QVZ+5LB+3GfJzu/NxmyE/tztT25ytsbluBIqAh8MZvivCmVtnAl8ys24gCVzl7rvCe/4W+BFQQnSMpe84yw3A7Wb2UeAV4OLR2ggREYlkJUzcff5h2u8E7jzMvHrg2AHadwLnjGiBIiIyJGPhbK7x4qZsF5AF+bjNkJ/bnY/bDPm53RnZ5rwdNVhEREaO9kxERCRtChMREUmbwuR1mNkyM1sXBpjM2XG/zGy2mT1qZi+a2Roz+1Roz/mBNM0sHkZduCdMzzOzp8J3/nMzGwf3EhwaM6syszvCgKtrzezUXP+uzezvwn/bq83sVjMrzsXv2sxuNrMdZrY6pW3A79Yi3wrbv8rMThjuehUmR2BmceDbwPnAIuBSM1uU3aoypgf4rLsvAk4BPh62NR8G0vwUsDZl+ivA18NZhy3AR7NSVWZ9E3jA3d8ALCHa/pz9rs1sJvBJoM7djwXiwCXk5nf9Iw4OhNvncN/t+RwcPPdKogF1h0VhcmQnAQ3uvt7du4DbiAaWzDnuvtXdnw2vW4l+XGaS4wNpmtks4J3AD8K0AWcDd4RFcnGbJxBd0/VDAHfvcvfd5Ph3TXQpRImZJYBSouGbcu67dvffAbv6NR/uu10O/NgjK4CqMJLIkClMjmwmsCllOnWAyZxlZrXAUqLRBnJ9IM1vEI3GkAzTk4HdKaNa5+J3Pg9oAv4rdO/9wMzKyOHv2t03A18DXiUKkT3ASnL/u+5zuO92xH7jFCZyCDMrJ7pw9NPuvjd1XhhDLWfOJTezC4Ad7r4y27WMsgRwAvAdd18K7Kdfl1YOftcTif4KnwfMAMp4bVdQXsjUd6swObLNwOyU6dQBJnOOmRUQBclP3f0XoXl7327v6w2kOQ6dDrzLzDYSdWGeTXQsoSp0hUBufueNQKO7PxWm7yAKl1z+rs8FNrh7k7t3A78g+v5z/bvuc7jvdsR+4xQmR/YMsCCc8VFIdMDu7izXlBHhWMEPgbXu/u8ps/oG0oQhDKQ5Hrj7te4+y91rib7b37j7XwKPAheFxXJqmwHcfRuwycyOCU3nAC+Sw981UffWKWZWGv5b79vmnP6uUxzuu70buDyc1XUKsCelO2xIdAX86zCzdxD1q8eBm939+uxWlBlm9hbgceAFDh4/+DzRcZPbgTmEgTRTBt/MGWZ2FvD37n5BGLH6NmAS8BzwgXB7hJxhZscTnXRQCKwHPkz0x2XOftdm9kXg/URnLj4H/BXR8YGc+q7N7FbgLKKh5rcD/wz8kgG+2xCsNxJ1+bUBHw7jIA59vQoTERFJl7q5REQkbQoTERFJm8JERETSpjAREZG0KUxERCRtChORNJnZ78NzrZldNsKf/fmB1iUy1ujUYJERknqtyhDek0gZG2qg+fvcvXwEyhPJKO2ZiKTJzPaFlzcAZ5jZ8+HeGXEz+6qZPRPuFfE3YfmzzOxxM7ub6CpszOyXZrYy3G/jytB2A9Eot8+b2U9T1xWuWP5quDfHC2b2/pTP/m3KvUp+Gi5ME8moxOsvIiKDdA0peyYhFPa4+4lmVgQ8aWYPhWVPAI519w1h+iPhiuQS4Bkzu9PdrzGzq939+AHW9V7geKJ7kVSH9/wuzFsKLAa2AE8SjUH1xEhvrEgq7ZmIZM7bicY9ep5oWJrJRDchAng6JUgAPmlmfwRWEA28t4Ajewtwq7v3uvt24DHgxJTPbnT3JPA8UDsC2yJyRNozEckcAz7h7g8e0hgdW9nfb/pc4FR3bzOz3wLFaaw3dWypXvT/uYwC7ZmIjJxWoCJl+kHgY2Fof8xsYbgJVX8TgJYQJG8gum1yn+6+9/fzOPD+cFymhujOiU+PyFaIDIP+YhEZOauA3tBd9SOie6PUAs+Gg+BNDHxb2AeAq8xsLbCOqKurz03AKjN7NgyP3+cu4FTgj0Q3OvpHd98Wwkhk1OnUYBERSZu6uUREJG0KExERSZvCRERE0qYwERGRtClMREQkbQoTERFJm8JERETS9v8Bb+rX1oZ3hjIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    \n",
    "    # 初始化可学习参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "    def __init__(self, shape, dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # 创建用来学习 U 的参数 theta\n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype,is_bias=False)\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的参数 phi\n",
    "        self.phi = self.create_parameter(shape=shape,\n",
    "                                         default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                         dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 将 Numpy array 转换成 Paddle 支持的 Tensor\n",
    "        self.M = paddle.to_tensor(M)\n",
    "        self.weight = paddle.to_tensor(weight)\n",
    "\n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self):\n",
    "        \n",
    "        # 获取量子神经网络的酉矩阵表示\n",
    "        U = U_theta(self.theta)\n",
    "        U_dagger = dagger(U)\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 -= paddle.real(self.weight)[i] * paddle.real(matmul(U_dagger,matmul(self.M, V)))[i][i]\n",
    "            singular_values[i] = paddle.real(matmul(U_dagger, matmul(self.M, V)))[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",
    "import time\n",
    "start = time.time()\n",
    "    \n",
    "# 确定网络的参数维度\n",
    "net = NET([theta_size])\n",
    "\n",
    "# 一般来说，我们利用 Adam 优化器来获得相对好的收敛\n",
    "# 当然你可以改成 SGD 或者是 RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(ITR):\n",
    "\n",
    "    #  前向传播计算损失函数\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # 反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 记录优化中间结果\n",
    "    loss_list.append(loss[0][0].numpy())\n",
    "    singular_value_list.append(singular_values)\n",
    "    \n",
    "    if itr% 10 == 0:\n",
    "        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "        \n",
    "# 绘制学习曲线\n",
    "loss_plot(loss_list)\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = U.numpy()\n",
    "V_dagger_learned = V_dagger.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": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.407257Z",
     "start_time": "2021-03-09T03:49:24.027591Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFj0lEQVR4nO3dd3hUZfbA8e9JIQk1lBB6L6EHiIAiHcEGKlJk1QXr6lp29QfWXcV1XbGsbXV17YooKM3e6CKC0nsnIC3UQKhp5/fHewkBUiYwyaScz/PMk5l779w5M0nm3LeLqmKMMabkCgp0AMYYYwLLEoExxpRwlgiMMaaEs0RgjDElnCUCY4wp4UICHcC5qFKlitarVy/QYRhjTJGycOHCvaoadeb2IpkI6tWrx4IFCwIdhjHGFCkisiWr7VY1ZIwxJZwlAmOMKeEsERhjTAlXJNsIjCkMdu/ezYgRI1izZg3p6emBDscYAIKCgoiJieH555+natWqPj3HEoEx52jEiBH06NGDd955h9DQ0ECHYwwAKSkpjBkzhhEjRvDhhx/69JwSUzU0ZfF2+j89gfmPdaTf0xOZsnh7oEMyRdyaNWu44YYbLAmYQiU0NJQbb7yRNWvW+PycElEimLJ4Ow9PWs4j+gkXBK9l0JGPeXhSGQCublszwNGZoio9Pd2SgCmUQkND81RdWSJKBM99v5ayKXsZFDyLIFEGBc+mbMo+nvt+baBDM8aYgCsRiWBH4jHuDZlMMGkABJHGPSGT2JF4LMCRGXN+goODiY2NpWXLlvTr14/ExMSAxTJz5kzmzp3rt/NNmTKFVatWZTx+7LHHmDp1qt/Of6YvvviC0aNH+3Ts0aNHqVy5MocOHTpt+9VXX8348eMBF3/r1q2JiYmhZcuWTJgwIeO4efPm0bFjR2JjY2nWrBmjRo0iPj6eWrVqnXUlHxsby/z58xk1ahQ1a9YkNjaWxo0bM2DAgNM+n/NRIhJBqwrHGBQ8i1BxH3CYpDEoeDYtKxwPcGSmJJmyeDudR0+n/kNf03n0dL+0U0VERLBkyRJWrFhBpUqVeO211/wQ6bnJKRGkpqbm+XxnJoJ//OMf9O7d+5zjy03//v156KGHfDq2dOnS9O3bl8mTJ2dsO3jwIHPmzKFfv34sXbqUESNG8Pnnn7NmzRq+/PJLHnzwQRYuXAjAsGHDePPNNzN+d4MHD6ZevXrUqVOHn376KeOca9asISkpiY4dOwJw3333sWTJEtavX8+QIUPo2bMne/bsOe/3XiISwUvVf0Q4fSW2INJ5ufoPAYrIlDQn26m2Jx5Dge2Jx3h40nK/dlq48MIL2b7dnW/jxo1ceumltG/fni5dumQ0HCYkJHDNNdfQpk0b2rRpk/HF/cILL9CyZUtatmzJSy+9BEB8fDzNmjXjtttuo0WLFvTp04djx1wp+pVXXqF58+a0bt2a6667jvj4eN544w1efPFFYmNj+emnnxg+fDh33HEHHTt25IEHHmDUqFE8//zzGfG2bNmS+Ph4AD788ENat25NmzZtuPHGG5k7dy5ffPEFI0eOJDY2lo0bNzJ8+PCMq+pp06bRtm1bWrVqxc0338yJEycAN/3M448/Trt27WjVqlWWDaadOnVi5cqVGY+7d+/OggULeP/997n77rsB+PLLL+nYsSNt27ald+/eJCQknHWeoUOHMm7cuIzHkydPpm/fvpQuXZrnn3+eRx55hPr16wNQv359HnnkEf79738Drutx9erVAVeqa968eZbnHDduHNddd12Wv+8hQ4bQp08fPv744yz350WJaCxucHwlyOlXJGGS6rYb4ydD/vdLtvsWb00kOe30Iv+xlDRGfbmSq9vWZP+RZO78aOFp+8f/6UKfXzstLY1p06Zxyy23AHD77bfzxhtv0LhxY+bPn8+f//xnpk+fzr333ku3bt2YPHkyaWlpHD58mIULF/Lee+8xf/58VJWOHTvSrVs3KlasyPr16/nkk0946623GDx4MBMnTuSGG25g9OjRbN68mbCwMBITE4mMjOSOO+6gbNmyjBgxAoB33nmHbdu2MXfuXIKDgxk1alSWsa9cuZJ//vOfzJ07lypVqrB//34qVapE//79ufLKKxk4cOBpxx8/fpzhw4czbdo0mjRpwh//+Edef/11/vrXvwJQpUoVFi1axH//+1+ef/553n777dOeP2TIED799FOeeOIJdu7cyc6dO4mLi2PFihUZx1x88cXMmzcPEeHtt9/m2WefzfgSP6lv377ceuut7Nu3j8qVKzNu3LiMRLJy5cqMz+GkuLg4/vOf/wDuyr5p06Z0796dSy+9lGHDhhEeHs7gwYOJjY3lP//5DyEhIYwfP57PPvss2997u3bt8tQ7KDslokTAHXNg1EEYdZCD96zlhIbyS+UBbrsxBeDMJHBS4tGU8zrvsWPHiI2NpVq1aiQkJHDJJZdw+PBh5s6dy6BBg4iNjeVPf/oTO3fuBGD69OnceeedgLsSrVChAnPmzOGaa66hTJkylC1blgEDBmRUT9SvX5/Y2FgA2rdvn3EF37p1a66//no++ugjQkKyv54cNGgQwcHBOb6H6dOnM2jQIKpUqQJApUqVcjx+7dq11K9fnyZNmgCummX27NkZ+wcMGHBWvJkNHjw4o2Tx6aefnpVoALZt20bfvn1p1aoVzz333GkliJNKlSpF//79mTBhAnv37mXx4sX07ds3x9hPeuyxx1iwYEHGFf2ll14KQHR0NC1btmTatGksWbKEkJAQWrZsme15/LXmfIkoEWRWoXI1JtS4l0nbI2l1IpWyYSXuIzD5JKcr+M6jp7M9i84JNSMjAKhUplSeSgAnnWwjOHr0KH379uW1115j+PDhREZGsmTJkjyf70xhYWEZ94ODgzOqhr7++mtmz57Nl19+yVNPPcXy5cuzfH6ZMmUy7oeEhJzWEHr8eP600Z2MOTg4OMu2iZo1a1K5cmWWLVvG+PHjeeONN8465p577uH++++nf//+zJw5M9vSzNChQ3nyySdRVa666qqM7sTNmzdn4cKFtGnTJuPYhQsXEhcXl/G4YcOG3Hnnndx2221ERUVllCxOVg9FR0czdOjQHN/r4sWLTzvnuSrQEoGIRIrIBBFZIyKrReRCEakkIj+KyHrvZ8X8jqPhZfcwN7khk21QmSkgI/s2JSL09CvjiNBgRvZt6pfzly5dmldeeYV///vflC5dmvr162dUKagqS5cuBaBXr168/vrrgKtOOnjwIF26dGHKlCkcPXqUI0eOMHnyZLp06ZLta6Wnp/P777/To0cPnnnmGQ4ePMjhw4cpV64cSUlJ2T6vXr16LFq0CIBFixaxefNmAHr27Mlnn33Gvn37ANi/fz9Atudr2rQp8fHxbNiwAYAxY8bQrVu3PH1eQ4YM4dlnn+XgwYO0bt36rP0HDx6kZk03xuiDDz7I9jzdu3dn/fr1vPbaa6d9aY8YMYKnn346o0QSHx/PSy+9xMiRIwGXSE9eza9fv57g4GAiIyMBV6L55ptvGD9+fLbtAwATJ07khx9+yDVZ+KKgq4ZeBr5T1RigDbAaeAiYpqqNgWne43wVWzuSS6MPkjrzeb8VrYzJydVta/L0gFbUjIxAcCWBpwe08uuAxrZt29K6dWs++eQTxo4dyzvvvEObNm1o0aIFn3/+OQAvv/wyM2bMoFWrVrRv355Vq1bRrl07hg8fTocOHejYsSO33norbdu2zfZ10tLSuOGGG2jVqhVt27bl3nvvJTIykn79+jF58uSMxuIzXXvttezfv58WLVrw6quvZlTttGjRgkcffZRu3brRpk0b7r//fgCuu+46nnvuOdq2bcvGjRszzhMeHs57773HoEGDaNWqFUFBQdxxxx15+qwGDhzIuHHjGDx4cJb7R40axaBBg2jfvn1GlVVWgoKCGDhwIPv27TstGcXGxvLMM8/Qr18/mjRpQpMmTXj99ddp2tQl/jFjxtC0aVNiY2O58cYbGTt2bEYVWmRkJBdeeCHR0dE0aNDgtNc72RjfuHFjPvroI6ZPn05U1FnrzOSZFNQXoYhUAJYADTTTi4rIWqC7qu4UkerATFXN8TIpLi5Oz3dhmt8+e5YLVj7Fhqu/pFFs1/M6lymZ4uLibIEk45OHHnqI+fPn8/3331OqVKkCec2s/j5FZKGqnlWXVJAlgvrAHuA9EVksIm+LSBkgWlV3esfsAqKzerKI3C4iC0RkgT/6zTbvcyvpIRE02vLpeZ/LGGNyMnr0aGbMmFFgSSCvCjIRhADtgNdVtS1whDOqgbySQpZFFFV9U1XjVDXOH0WhMhUqEdR6EKyYCMcSz/t8xhhTVBVkItgGbFPV+d7jCbjEkOBVCeH93F1QAaW0uwlSjjL/89cL6iWNMabQKbBEoKq7gN9F5GT9fy9gFfAFMMzbNgz4vKBiCq3VjtVhbUg7dij3g40xppgq6E709wBjRaQUsAm4CZeMPhWRW4AtQNbN+Pkk5sGZSFDJGFdnjDFZKdBEoKpLgKxGP/QqyDgyk6AgND2drZtWUbdR9iP4jDGmuLJLYWDBR3+j2pjuJOyyAWamaLFpqP0nL9NQg5uK+vrrr6dVq1a0bNmSiy++mMOHD9OjRw++//7704596aWXuPPOO4mPjyciIoK2bdvSrFkzOnTowPvvv+/nd5J3lgiAmp2uJUxSWPPd/wIdiinuknbBe5dB0tmzWZ4Lm4baf/IyDTW4wXnR0dEsX76cFStWZKxdfeYMouBmET05Arhhw4YsXryY1atXM27cOF566SXee+89v76XvLJEANRo0p61YS2pH//pOf3BGuOzWc/C1nkw6xm/n9qmoS7Yaah37tyZMQ0FuKkvwsLCGDhwIF9//TXJyckZn+OOHTuynLajQYMGvPDCC7zyyiu5/Xrzl6oWuVv79u3V35Z+/T/Vx8vr/KmT/H5uUzyd9Xf47uVn3+a/6fadOKL6Zm/VUZGqj5d3P9/qrbroI7f/8N6zn+uDMmXKqKpqamqqDhw4UL/99ltVVe3Zs6euW7dOVVXnzZunPXr0UFXVwYMH64svvpjxnMTERF2wYIG2bNlSDx8+rElJSdq8eXNdtGiRbt68WYODg3Xx4sWqqjpo0CAdM2aMqqpWr15djx8/rqqqBw4cUFXVxx9/XJ977rmM2IYNG6ZXXHGFpqamZrm/RYsWunnzZl2xYoU2btxY9+zZo6qq+/bty3j+Z599dtr5PvvsMz127JjWqlVL165dq6qqN954Y8Z7qlu3rr7yyiuqqvraa6/pLbfcctZn9sILL+hjjz2mqqo7duzQJk2aqKrqe++9p3fddZeqqu7fv1/T09NVVfWtt97S+++//6zzLF68WKOiorRTp0766KOPZnzeqqpXXHGFTpkyRVVVn376af2///s/VVXdvHmztmjR4rTzHDhwQMPDw886//nK6nsSWKBZfKdaicDTvNcNJFKO5IUfBToUU1wd3AonZ1dRhcSt531Km4Y6cNNQx8bGsmnTJkaOHMn+/fu54IILWL16NXD6AjOZq4WyooVgvjObg9kTElaab2P/y6h5aXy/9wj1qpTJ/UnGZHbT19nvO3EIjidyauC8useNvDrvMpVzfn42bBrq7GMuiGmoTybOAQMGEBQUxDfffEOzZs246qqruO+++1i0aBFHjx6lffv22ca7ePFimjVrdm5v1k+sRJBJz559SA0K4+Nfz/9KzZjTzHoW9IzFaTTdb20FNg11wU9D/fPPP3PgwAEAkpOTWbVqFXXr1gVcgujRowc333xzjqWB+Ph4RowYwT333JOn+P3NEkEm0eXDeaD2Wnr9eivHTyQHOhxTnGz7FdLO+JtKS3bb/cSmofadP6ah3rhxI926dcv4HOLi4rj22msz9g8dOpSlS5eelQg2btyY0X108ODB3Hvvvdx00015it/fCmwaan/yxzTU2Vk9bQzNfrqbrZe+T51O1+TLa5jiwaahNoVZYZ2Gukho2nUIWjaaOpvG5X6wMcYUA5YIzhAUWgppeyO67nuSEjYFOhxjjMl3lgiykBp7IwosmPRSoEMxhVhQUBApKSmBDsOYs6SkpBCUh8k0LRFkIaRyPRbXu43IJtn3nDAmJiaGMWPGWDIwhUpKSgpjxowhJibG5+dYY7Ex52j37t2MGDGCNWvWnNY/3phACgoKIiYmhueff56qVaueti+7xmIbUJaDPds3s3bul3QeeA8iEuhwTCFTtWpVPvzww0CHYcx5s6qhHOyY+RYXr/w7K1csCXQoxhiTbywR5KBx3ztJ1SB2zzx7CLoxxhQXlghyULpKbdZGdqHt3q85cDD7ofPGGFOUWSLIRbkuf6KiJLHkB6sLNsYUT5YIclGn3WVsC67N7+uWkJ5e9HpYGWNMbiwR5CYoiEWXfcljSVczd+O+QEdjjDF+Z4nAB33a1KFi6VAmzT17cQpjjCnqLBH4IDw0mBdqzODxTUNJ2GelAmNM8WKJwEfNOlxCBTlCyKopgQ7FGGP8yhKBj6q16gFRMVRebWsaG2OKlwJNBCISLyLLRWSJiCzwtlUSkR9FZL33s2JBxuQzEYi7GXYsYuPSs1dfMsaYoioQJYIeqhqbaeKjh4BpqtoYmOY9LpTSWg3mOKX4/cfXAx2KMcb4TWGYdO4qoLt3/wNgJvBgoILJSXDpimy65E3aNOkU6FCMMcZvCrpEoMAPIrJQRG73tkWr6k7v/i4gOqsnisjtIrJARBbs2bOnIGLNUuPO11AxqnrAXt8YY/ytoBPBxaraDrgMuEtEumbeqW5xhCyH76rqm6oap6pxUVFRBRBq9tb8NJGfnhvEiZTUgMZhjDH+UKCJQFW3ez93A5OBDkCCiFQH8H7uLsiYzkVw0g66HPmB+T99H+hQjDHmvBVYIhCRMiJS7uR9oA+wAvgCGOYdNgz4vKBiOlcNew7nCBGk//pOoEMxxpjzVpAlgmhgjogsBX4FvlbV74DRwCUish7o7T0u1ILCyxFf80ouPDab9fFbAx2OMcaclwJLBKq6SVXbeLcWqvqUt32fqvZS1caq2ltV9xdUTOejdu+7CJMUNvz4ZqBDMcaY85KnRCAicSIyxKvaOVndUxi6oBa48vXbMq/CFfy4LYQjJ6zR2BhTdPmUCEQkWkTm4ap0PuZUF88XgH/nU2yFXuiAV5l04gK+WLoj0KEYY8w587VE8CKQAFQGjmba/hmu0bdEalenIu2ig1k1exKu56sxxhQ9viaCXsCjqnrgjO0bgTr+DanoEBH+WfEbHkv6ByvWrQ90OMYYc058TQQRQHIW26OA4/4Lp+ip2+dOQiWNmJ1fBDoUY4w5J74mgtnA8EyPVUSCcXMCTfN3UEVJmRrNoF4XQpd8AOnpgQ7HGGPyzNdE8ABwm4j8CIThGohXAZ2Bh/MptiIjvf3NkLiVn38YH+hQjDEmz3xKBKq6CmgFzAV+AMJxDcVtVXVj/oVXNAQ1u5JEieT4+tmBDsUYY/LM5zEAqroLeDwfYym6QkoReu+v9KqY5cSpxhhTqPk6juBuEbkhi+03iMif/R9W0VPGSwJHjx0LcCTGGJM3vrYR/BX4PYvt8cB9/gqmqFvx2T/ZO7oNuxMPBzoUY4zxma+JoBawJYvt27x9Boiq24w6ksBvP3wS6FCMMcZnviaCXUBsFtvbAXv9Fk0RF93+KvYHV6HymrGkpdtIY2NM0eBrIvgYeEVELhGRUO/WB3gJGJtv0RU1wSHsa3odndIXM2/BgkBHY4wxPvE1ETwO/Ax8j5tr6CjwLa476d/zJ7Siqd4ld5JKEIlz3gp0KMYY4xOfuo+qagowVEQeA9ri1hVeoqo2wc4ZQivW4rtGf+O5lZG03n+U2pVKBzokY4zJUZ7WI1DV9ar6qap+Zkkge236/ZmtVOPjX231MmNM4efzgDIRGYKbhbQqZyQQVe3v57iKtOoVIrijXgKV5k8kuffblAopyBVBjTEmb3wdUPYc8BFQD0gE9p1xM2cYGL2DW3Uiq5b9FuhQjDEmR76WCP4IDFXVCfkZTHFSr9ft6LKXiN09BegY6HCMMSZbvtZZBAFL8jGOYieoXFWkeX9YMpb0E0dzf4IxxgSIr4ngTeCsuYZMzrT9cDh+kM8/fi3QoRhjTLZ8rRqKBP4gIpcAy4CUzDtV9V4/x1UsSL0ubCp/AZXKhAY6FGOMyZaviaA5p6qGYs7YZ3MpZEeEBvdPpUGg4zDGmBz4OqCsh79e0FvicgGwXVWvFJH6wDigMrAQuFFVs1ofuchKST7B0sXzaN+hKyIS6HCMMeY0gejg/hdgdabHzwAvqmoj4ABwSwBiylebxv6FmG8Gs2LT9kCHYowxZ/E5EYhIDxF5U0S+E5HpmW95OEct4Argbe+xAD2Bk91SPwCu9jn6IqJW12GUleOsm/puoEMxxpiz+DqgbDhukrlyQHdgD1ARNw31qjy83kvAA0C697gykKiqqd7jbUDNbGK4XUQWiMiCPXv25OElA69Mg07sCG9E8x0TOHikWNV6GWOKAV9LBCOAu1V1KK7H0MOq2hY32tin5bhE5Epgt6ouPJdAVfVNVY1T1bioqKhzOUXgiKDtb6KZbGH2zO8CHY0xxpzG10TQAJjq3T8BlPXuvwoM9/EcnYH+IhKPaxzuCbwMRIrIyUbrWkCxrEiv2XUYxySCo0smomodrYwxhYeviWAfrloI3Bd1S+9+ZSDClxOo6sOqWktV6wHXAdNV9XpgBjDQO2wY8LmPMRUtYeWY2e1THkoayC+bbHomY0zh4Wsi+Ano493/FLda2XvAJ8CP5xnDg8D9IrIBl1jeOc/zFVo9OnemfEQYY+dltfyzMcYEhq8Dyu4Gwr37TwOpuKqeT4F/5vVFVXUmMNO7vwnokNdzFEXhocGMrj2PqLVfsfvQLKqW96kwZYwx+crXAWX7M91Px/X9N+egfaOaVN26lk1rZ1P1gr6BDscYY3zuPpomIlWz2F5ZRNL8H1bxVbXTUDS8Ag22fBroUIwxBvC9jSC7eRHCAOsYnxelSiNthqKrPmfHdlvK0hgTeDlWDYnI/d5dBe4QkcxjBoKBLsCafIqt2NL2w5H5bzBnwisM/svzgQ7HGFPC5dZGcI/3U4BbgczVQMlAPHCH/8Mq3qRqM9a1uI9WDXsGOhRjjMk5EahqfQARmQEMUNUDBRJVCdBk0KhAh2CMMYCPbQSq2uPMJCAijUQkPLvnmNz9vno+X7/7JMmp6bkfbIwx+cTXXkP/EpFh3n0RkanAOmCniNjK7Ocofdln9N3yArMWLA10KMaYEszXXkPXA2u9+5cBbYBOwIfA6HyIq0So1evPhEg6B+YU28HUxpgiwNdEEI2bIhrgcuBTVf0V+A/QNj8CKwmCqzRga8VOXJz0NRt2WfOLMSYw8jLpXF3vfh9gmnc/hOzHGBgfRHb9EzVkP7/9MD7QoRhjSihfE8FE4GMR+RGoBHzvbY8FNuRDXCVG+db92FWqDus3redocmruTzDGGD/zNRHcD7yCW43sElU94m2vDryeH4GVGMGhbBkynXeP9+SrpTsDHY0xpgTyddK5VODfWWx/0e8RlUAdGlShSdUyfDd3AYMvqB3ocIwxJUy2iUBE2gFLVDXdu58tVV3k98hKEBHh5QqfUGXrtyzf0plWdc+a388YY/JNTiWCBUA1YLd3X8m6YVhx8w6Z81C309WU/n0cwQkzoO6QQIdjjClBckoE9YE9me6bfFS6WR+oUIdKqz+CDpYIjDEFJ9vGYlXdot4q6979bG8FF24xFhSMth8Gm2cz77d5gY7GGFOC+DrFRAMRuV9EXhWR/4jIfSJipQQ/k7Y3kkowO2faSGNjTMHJtdeQiPwfbp3iYFx7gQBRwDMi8qD1HPKjctHsu3YClzfuFOhIjDElSI4lAhG5GHgWeA6IUtXqqloNqIrrTvqciHTO/zBLjuhWPQkLL41XK2eMMfkut6qhO4EPVfXRMxaw36eqDwMfAX/OzwBLovXfv8H0p/qzJ+lEoEMxxpQAuSWCTsD7Oex/3zvG+FFFOUKv1Nk89O9Xmf9YR/o9PZEpi7cHOixjTDGVWyKoBmzKYf9G3DQTxo/mV7iUExrKyPR3uUDWMujIxzw8abklA2NMvsgtEUQAOdVPJANhvryQiISLyK8islREVorIE972+iIyX0Q2iMh4ESnlW+jF179mJDA1vS1NZRtBogwKnk3ZlH089/3a3J9sjDF55MtcQ1eIyMFs9kXm4bVOAD1V9bCIhAJzRORb3IR2L6rqOBF5A7iFEj6R3Y7EY0jIqcbiINK5J2QSjyXezP4jyVQqU+JzpTHGj3xJBLl1avepe4s3OO2w9zDUuynQE/iDt/0DYBQlPBG0qnCMnseXIN6EHmGSyuDgWSRoJbo9ncrX911CncqlAxukMabYyLFqSFWDfLj5PM+QiASLyBLceIQfcW0Mid7spuBWQauZzXNvF5EFIrJgz549WR1SbLxU/UfkjPwaQhojQz/ll9IjqL1pHKSl8MHceGas3R2gKI0xxYWv6xH4haqmqWosUAvoAMTk4blvqmqcqsZFRUXlV4iFQoPjKwmT0xepCZF0qFifslXrIV/fh74ax47Z7/Ltsh0ZxxxLTivoUI0xxYBP6xH4m6omisgM4EIgUkRCvFJBLcC6xtwxJ/t9qrD+R2T6kzxYYRGHrvgbAEt/T+SGt+cz+ILaDLuwnlUdGWN8VmAlAhGJEpFI734EcAmwGpgBDPQOGwZ8XlAxFUki0KQP3D6LoMEfElm6FBzcRtOvruHuWhv4YO5muj0/g9s/XMAvG/fZCGVjTK4KskRQHfhARIJxCehTVf1KRFYB40Tkn8Bicm+cNgBBQVC6krt/aCfhyfv5U8Ij3FS3PZMr3sToNaH8sCqBZtXLc1PnevRvU4PwUFs2whhzNimKV4xxcXG6YMGCQIdRuKSlwJKxMOtZOLSdtPrdmBjzEu/M3cbahCQqlynFHd0aclvXBoGO1BgTICKyUFXjztwekDYCkw+CQ6H9cGh9HSx8n+BD2xncsQGDOtRn0ZLFvL4sjcRjyQCkpyurdx2iRY0KgY3ZGFMo+LoeQUUReVlElonILhHZnfmW30GaPAgNh053QJ8nAZCdS2j/eU/eDn+FEW1d6W/mut1c8cocZq0r3t1wjTG+8bVE8CHQAjfgKwEfB5GZQqBSA+j2APzyGrL6S2g9mAsuHMGTV7XgooaVAfh4/lYOHU9h6AV1qFA6NMABG2MKmk9tBCKSBHRT1UX5H1LurI3gHBzZBz+/BL++5UoN96+G0AgA/jpuMVOW7CAiNJgB7WpyU+d6NKpaLrDxGmP8Lrs2Al8TwRLgVlUtFN++lgjOQ9Iu2LnMdUFVhZ9fhtg/sOpQOO/9vJnPl+4gOTWdbk2iuKlzPbo2jiIoSAIdtTHGD843EXQD/gaMAFaoakCHsFoi8JOdS+HNHhASBh3/BBfdy970Mnw8fytj5m1hT9IJGkaV4abO9RnQrialS1nfAmOKsuwSga8DyjbgpqReBCSLSFrmmz8DNQWoehu461doejnMeQlebkOVBS9xb5ea/PxgT14c0obSpUL425QVLNmaCGAD1IwphnwtEcwGKgJvkEVjsapOzJfosmElgnyQsBJm/MtVG92zwJUSVFFg2baDtK5VARHhn1+tYs/hE7w0JBYRqzIypig533EEcUAHVV3h37BMoRHdAq4bC8cPuiSQchzeuwxpM5Q27Ydxck7s8hGhpKSlZySBmWt3szfpBC9OXc+OxGPUiIxgZN+mXN02y0lkjTGFkK+JYBVQPj8DMYVEuDfI7OheCC0N346Eua+4Lqht/sC9vRpnHLp2VxLD3/vttKdvTzzGw5OWA1gyMKaI8LWN4G/ACyLSW0SiRaRS5lt+BmgCpEItGP4V3DgZylaFL+6B1zq4XkeexlXLUjmL1dKOpaTZsprGFCG+lgi+8X7+wOntA+I9ttnMiiMRaNgTGvSAtd/Cmq+hbLTbt2cdQVUas/9IcpZP3Z54jJ837OWihpWtLcGYQs7XRNAjX6MwhZsIxFzubgBH9sKb3aFKY64p149JSTFEkcirpf7D3cn3sodIggSuf3s+HepXYvztnSwZGFOI5ZoIvIXmnwX+qKpW3jcQHgmXPwezRvNCypP8ISyGfWlluEDWck/IJJ6W2/jHVS1IV+XIiTREBFVl7Pyt9Gtdw6axMKaQyTURqGqKiNTH5hcyJwWHQNvrodUgWPQBLab+i4jk/QAMDplNVN+/c1lc7dOesnTbQf42ZQXhocEMbF8rEFEbY7Lha2PxB8Bt+RmIKYJCSkGH24ho2Q+C3DVFeDBctvB2mPoEHIjPODS2diTf3NuFfm2qA/D2T5u4Y8xCFsTvt0FqxgSYr20EZYDrReQSYCFwJPNOVb3X34GZIiJpFywbD+mp7nFaMuxbD3NedLdGvSHuZmjch+Y1TvVAFhF+2bSP71buIrZ2JLd1aUDfFtGEBBfY6qnGGI+vI4tn5LBbVbWn/0LKnY0sLkS+uh8Wj3EJ4KTgUtBiAFSsB4s+gKSd0Otx6HL/aU89mpzKhIXbeGfOZrbsO0qtihHc1Lk+Qy6oTdkwm9fIGH87r0nnChtLBIXIGxfDruVnb6/WCu6YA2mpsO47qNkOytdw3VAXfwQX3AL1u0NQEGnpytTVCbzz02Z+jd9PubAQhnasw82d61OtQnhBvyNjii2/LFUpIuFAI1zD8UZVPe6n+ExRdcecnPcHh0CzK089Proftv4Ca76CivUh7iaCY6+nb4tq9G1RjSW/J/L2T5t4Z85mOjeqQrUK4aSlK8E2FbYx+cbXqqFQ4F/A3UAp3ECyE8B/gEdVNSU/gzyTlQiKuNQTsPpLWPAubPkZqjSFu+ZnzGcEsCPxGNUrhCMi/Oub1azccZAPbupgbQjGnIfzLRE8AwwF7gBOXgJ2AZ7G9Twa4Y8gTQkREgatBrrb7tWuwVnETXQ35hpocTU1Wg8BcSuo1a1cmrR0zUgCU1cl0LlRFSJK2YB2Y/zB1xLBLuBmVf3mjO1XAG+ravV8ii9LViIopvZvggm3wI5FbsK7lte6Hkc122UcsmnPYXr+exYVS4dyQ6e63HhhXaqWs3YEY3xxviuUHQNizxxZLCIxwGJVjfBbpD6wRFDM7VgMC96D5Z9BylG4fSbUaAu4hXF+iz/AWz9tYurqBEKDgri6bQ1uubgBTavZOsvG5OR8E8E8YKGq3nXG9tdxCeJCH85RG/gQiMY1Nr+pqi97s5eOB+oB8cBgVT2Q07ksEZQQxw+6Xkath7iqox8fg+SjEHcTRLdg894jvDtnM58t/J3jKel0bRLFrRfXp0vjKja3kTFZON9E0BU3A+l2YJ63uRNQA7hMVXPpOgIiUh2orqqLRKQcbmDa1cBwYL+qjhaRh4CKqvpgTueyRFBCfXW/63qadgJqd3LVRs2v4kByEGPnb+GDX9w6y+3qRDLxzossGRhzhvMeRyAiNYC7gBhv02rgv6q64xwD+hx41bt1V9WdXrKYqapNc3quJYIS7Mg+WPqx63G0fxN0vBMuGw3AidQ0vliyg4PHUri1S4OMie6uaFWdilmsm2BMSVOoBpSJSD1gNtAS2Kqqkd52AQ6cfHzGc24HbgeoU6dO+y1bthRUuKYwSk+H+NlQoTZUbghb58PMf7lSQtPLITiUFdsPcuV/5vD0gFYM7VAHgCmLt/Pc92ttWU1TIp1TIvB19TFV3Z+HQMoCs4CnVHWSiCRm/uIXkQOqWjGnc1iJwJxlzTfwzUg4tM0tntPuj9BuGGuPR1K3cmnCQ4O5/9MlTFm8nfRMf/IRocE8PaCVJQNTImSXCHIbnbMX2JPLbXcegggFJgJjVXWStznBqxI62Y7g8/mMyRBzOfx1GQwdD9VjYfbz8FYPmkZFEB7qxhtMW737tCQAblnNZ79fU/DxGlOI5DagLKeVyS4F/gKk+vJCXrXPO8BqVX0h064vgGHAaO/n576cz5izBAVD00vdLXEr7FnrprhIT4cxV3NDcjSf0h3Q01ZT25F4nFs/WEDnRpXp2iSKhlFlA/1OjClQeW4jEJG2wHO4kcX/A55U1T0+PO9i4CdgOZDubX4EmA98CtQBtuC6j+ZY1WRVQyZPjuyFCTfD5lmkaDDbtTJ1ZA9j03ry99RbKF0qmCplw9i6/ygD2tXkhcGxqCpfLtvJhQ0qE1UuLNDvwBi/OO9J57xVyp4CBgGTgOaqutHX53tdTLPrz9fL1/MYk2dlqsCwL/jxp5/Z++MLXCdTEYEhwbN4Uwbz0JVtuKJtfX5PSifVqzvatPcI936ymKeuacn1HeuyJ+kEi7ceoFPDypQPt6U2TfHiy5rFlYHHcPMM/QxcpKq/5XdgxvjbJV06s2njh6TEh1CKVBD4oOEMGiSth+feoHaj3hBzJZTtQ/3K5fnm3i4Z02DPWLObByYuI0igda1ILm5Uhc6NqtCubiRhITbnkSnacus19CgwEjfi9yFV/a6A4sqRVQ2Zc5K0C15uA6mZZk8PCYeB78H6H2DtN3A4AYJCoUlfGPJRxoyoyanpLN56gJ837GXOhr0s3XaQtHQlPDSIC+pVykgMzauXJ8imzDaF1Ll2H00HjgEzOFWvfxZV7e+PIH1licCck+xWU2t7I1z5gmtU3vabWyshLSVjoBpT/gxVmrjSQpVGACQdT2H+pv3M2bCXnzfsZf3uw4QGC0sf70PpUiEs33aQChGh1KlcOgBv1JisnWsbwYe4eYGMKfq2/Xp6EgD3eNuv7n5QENTp6G4nJR+F3atgyViY+jhExUDMFZRrfR29mzehd/NoABIOHWf1zkOULuX+pZ74ciUpael8fvfFAMzduJeYauWpZCOcTSFkS1Ua44vE313V0ZqvIP5n6PcytLsRDu+B3SuhbmcIPtWIvHHPYRKPptC+bkWOp6TR5okfOJGaTvPq5bm4satGuqBexYzEYUxBKFRTTJwvSwQmoI7ud1VKYWXh17fgmxEQHglNLoWYK6BRLyhVJuPwtHRl2bZEft6wl5837GPhlgMkp6UTGiy0q1PRtS80rkLrmhVsBTaTrywRGJMfko/Cxumw5mtY9y0cO+AW1fm/NRBewbU7BJ3+5X4sOY3f4vdnNDyv2nkIVfjn1S25oVNdDh5NYc/h4zSMKsvnS3bY3EjGb/yyeL0x5gylSkOzK90tLdWtwbxzqUsCAOP+AMmHXUNzzOUQWYeIUsF0bRJF1yZRAOw/kswvG/fRvq6bYuuHVbsYOWEZD17alFembeBYShoA2xOP8fCk5QCWDIxfWYnAmPw061lYMQn2rHaPq7eBDn+Cttdn+5SEQ8eZvW4PL01dx/bE42ftr1oujPmP9LL1Fkyeneukc8aY89HtAbhrHtyzCC75BwSHubEKAMlH4Ie/w9Z5kJ6W8ZTo8uEMiqvNjiySAMDupBN0eXYGo75Yyc8b9pKcmm3PbmN8YiUCYwqaqhuotmUufNAf0lOgTJRbRyHmSmjQDULC6Dx6OtsTjxHFgdMmyYuMCCWuXkV+Wr+XE6nplAsL4cYL6/LApTG5v7Yp0axEYExhcbJKp+5F8MAmuPYdqHcxrJgIHw+CvesA+HvXilQJPcG9IZO5QNZyT8gkIkKDGdW/BW8Pu4Alj/XhrT/GcXmr6hnjE46npDH8vV/5ecPeQL07UwRZY7ExgRReHloNdLfUE66xObolAJfufps+wZ+gpBOEMiRkFlF9/85lXkNxRKlgLmkezSXeoDaAXQePsyPxGMlprrpo9c5DfLF0B72bRdO2dqRNf2GyZFVDxhRW2xbCl/dCwopT2yo1gHsX5/pUVUVEGP/bVh6ZvIK0dKVK2VL0jKlK72bRXNy4ig1mK4Gs+6gxRU2FmrBvw+nbErdCUoJrU/hkiKteatbfrducyckeRUMuqMOlLaozc91upq7ezbfLd/Hpgm2EhQRxcaMq9G4eTa+YqlQtH15Q78oUQlYiMKawymmSvG4PwCfXwQ6vdBDdEpr1g9g/QGSdbE+ZnJrOb/H7+XFVAlNXJ7DtwDGiyoXxq9cddU/SCaqULWVdU4spKxEYU9TkNEleuWpw+0xXQlj9Faz+AmaOhlpxLhEkbnUrs9Voe6pxGigVEkRnb8rsx/s1Z13CYbYdOIqIoKpc9eocujSO4pmBrQFITUu3aS9KACsRGFNcJO2C0pXd5HdTn4A5L0CF2q6k0Kwf1O7o1nXORmpaOhMWbqNmxQi6NI5i58Fj9HlxNt2aRHFJ82i6N6lKhdK2OltRZnMNGVOSHN0P676DVV+4uZDSTkDF+m5gW5BvV/i/7z/Kq9M3MG1NAnsPJxMcJHSoV4nezaPp3awqdSuXyf0kplCxRGBMSXUiya3AlpQAF/7ZbXv/SldaaN4fGvSA0Owbi9PTlSXbEpm2OoGpq3azNiEJgCbRZenVLJrbuzSgYplSTFm83SbIK+QsERhjnNRk+OIeWPstnDgIpcpC4z7Q4Xaoe2GuT9+67yhTV7vG5qW/J/Lro735cVUCD0xYljF+ASAiNJinB7SyZFCIWCIwxpwuNRniZ8PqL9002r2fcJPhHdoJm2ZC00shomKOpziWnEZEqeCM6TDOVCEilHeGxRFTvTxlw6xvSqBZIjDGZC89zd1CSsGCd+Gr+yAoBOp3deMUYq6AslWzfXr9h77OcU1bEahXuQzNq5eneQ13a1GjPFXL2fiFgmTdR40x2QsKPtWjqN1wqNYGVn/uGpu/+qtbhW3kBldCSE12CSOTGpERWZYIqpUP56lrWrJqxyFW7jjE8u0H+Xr5zoz9N3Wux+P9WpCerny7YhcX1Ktog9sCoMASgYi8C1wJ7FbVlt62SsB4oB4QDwxW1QMFFZMxJgtBQVCrvbv1fgISVsL2haeqiT4ZAscPet1S3ajmkX2b8vCk5ZRN2ZsxU+rh0Mo8dFkMvZpF06vZqfmQDh5LYc1OlxiaRJcDYOv+o9z18SJGD2jFdR3qsGXfEd6Zs5nm1cvTokYFmlQrS1hI9l1fzfkpsKohEekKHAY+zJQIngX2q+poEXkIqKiqD+Z2LqsaMiaA5v7HLbazY5F7XLUFdLqTKdKT9K/u5+rU75kc0pfgK1/wuaE4JS2dtbuSqFYhnCplw5i1bg93jV3E4ROpAIQECY2qlnXVSl5yaF69vI1ryKNC0UYgIvWArzIlgrVAd1XdKSLVgZmq2jS381giMKYQSPwd1nzlqo+aXgqth8DLrd0sqsGl4O6FUDH76S5yk56ubN1/lFU7D7Fyx8GM6qXdSScyjvn8rs60qR3JuoQktuw7SrcmUZQKsZHQ2SmsbQTRqnqywnAXEJ3TwcaYQiSyNnS6091U4ev/O7XSWloyvBIL9btA/W5uDqRy1fJ0+qAgoV6VMtSrUobLW1XP2L4n6QSrdh5i1Y5DNKpaFoDJi7fz1uxNrPxHXwDGzt9C/N4jruRQozwNqpQ5a6oMG/dwSqATQQZVVRHJtngiIrcDtwPUqXPuVxnGmHxwOAGWjIX01NO3H9wB056Appe5RBD/M+xeBQ26Q+VGp82D5KuocmF0KxdFtyZRGdvu7tGIK1pVz2hHWLnjEBMWbMsY1xAWEkRMtXI0r1GBFjXKs/fwCd6YtZHjKW7/9sRjPDxpOUCJTAZWNWSMOX85zZTa/SE3bbYIfPcIzHvN7S9Xwy3LWb+bq1byceoLX6WkpbNxz2FWbj90WvXSoeOp2T6nZmQEPz/U069xFCaFtWroC2AYMNr7+XlgwzHGnJOcZkrNPP6g71NwwS2weRZsmgXrvoet8yB2qNu/+CMIj3RLd0ZEnldIocFBxFQrT0y18lzrbVNVth04RpdnZ2T5nB2Jx1BVbnznV+pULk1MtXI0jS5HTLXi3TBdkN1HPwG6A1VEZBvwOC4BfCoitwBbgMEFFY8xxo/umOPbcSJuEZ3KDSHuZkhPh8O73D5VmPWMm0JbgtwU2vW7ucFstc66iD0nIkLtSqWpmc24hxqRESSdSCU5LZ2vl+3k4/lbM/ZVKx9O02rliKlejphq5YirW4nalUr7Ja5AK7BEoKpDs9nVq6BiMMYUMkFBUL6Guy/iehptX+CmuNg0C+a+4noh1YqDtBT45VWo1xVqxOY4pXZuTo57OJaSlrEtIjSYkX2bUj48lE//dCGqSsKhE6zedYi1u5JYuyuJNbuS+GXjPpLT0nnk8hhu79qQhEPH+efXq/lT1wa0rFmB1LR0goOkSC3uE+iqIWOMOSWklFt+s+5F0OMRN3NqqtddNGElTB3l7odVcNVHDbq7GVTz2CPpZINwTr2GRIRqFcKpViGcHk1PVW+lpKUTv/cIFSJcVVHCoeMs3nqA415S+X5lAg9OXEaT6LLEVC9fJKqXbK4hY0zRcXg3bJ7tSgybZ7lqpGFfuW6qCSth5zLXAH2ylBEAS39PZMLCbV4J4vTG6eoVXPVS02rluL1LAyqXDSvQ2AprY7ExxviubFVoNdDdAPZvhvLeVfzKKTD7WXe/cmOXEBp0hyaXulXbCkib2pG0qR0JuMbpXYeOs+Zk1dLOQ6zZlcTcDfu4s1tDAF6dvp4vl+7k63svJiQ4iA27DxMWEkStihEZ1Uv5PebBEoExpuiqVP/U/e4Pu2qiTbNciWHJJ7D8M3hgs9u/5hsIjYA6ndxPcMt7TrgJBr4P5fw/nlVEqF4hguoVIs6qXgr1BrjVrlSadnUjMwa8PfnVKmat20O5sBCaVCtHqRBhQfwBUtJc7U1+jHmwqiFjTPGUmgwHNkOUNzTpvxe6wWzBYVC7gysx7Fru1mNofxNc+UJg4/Us33aQZdsTMxqnf4vfT1Zf0+cy5sGqhowxJUtIqVNJAOCWH2DLL94Yhpkw/Z+um6qmu1HRFWpDkz4Q1czvg9vyolWtCrSqVSHjcf2Hvs7yuB1ZdH89V5YIjDElQ1g590XfpI97PPlOV3Wk6W6OpGmj3K10FW+OpK7Q5DIoXz2ns+a77NZ6qBEZ4bfXsGn6jDElT9IuWDkJ0lPc4/QUCAmDvk9Do95utPNX98Hv89z+A/GuzeHg9gIPdWTfpkSEnj5m4uSYB3+xEoExpuSZ9awrCWSmCvs2wID/efc3nmpAXvsdfOctlVKpoSstNOjmeiSF+u/KPCu+jHk4X5YIjDElT05zI4Eb5Vyl0al9HW53A9g2z3LjGJZPcJPsPbjF7d8w1Y18rtsZwsv7Pdyr29bM11lRLREYY0oeX+dGOikoCKq1dLcL74K0VNi7DsLcegj8/IpLEhLszZHUFRr1csmjCLBEYIwxeRUcAtHNTz3+w6ew7TdXWtg8282RtHPJqUTw2zsQ3QJqtHO9mQoZSwTGGHO+QsO9nkZdgEfhxGE4utftO34IvhkJmgahZaDuha7EEHOlm4W1ELBEYIwx/hZW9lS1UXh5GLkB4uecKjH8+BiElnaJICkBVn3ukkNU03Nate18WSIwxpj8VrqSm/6ieX/3OGmX664KsOVn+Haku1822iWE+l2h+VUQfmpgWX5Oh2GJwBhjClrmabNbDoCa7U6VFjbPdgPdGvZyiWDzbJcENk534xtmPeP36TAsERhjTKBVrOdu7f54agxDBa+76OKPYNn4U8cuGQvdHvRrqcBGFhtjTGFy5hiGq1+HZle5rqngBsLNesavL2mJwBhjCrMje2D9967XEbiBb0vGukZmP7FEYIwxhVmW02H4t1RgicAYYwqz3KbD8ANrLDbGmMIsr9NhnAMrERhjTAlnicAYY0o4SwTGGFPCWSIwxpgSzhKBMcaUcKKqgY4hz0RkD7DlHJ9eBdjrx3DyW1GK12LNP0Up3qIUKxSteM831rqqGnXmxiKZCM6HiCxQ1bhAx+GrohSvxZp/ilK8RSlWKFrx5lesVjVkjDElnCUCY4wp4UpiIngz0AHkUVGK12LNP0Up3qIUKxStePMl1hLXRmCMMeZ0JbFEYIwxJhNLBMYYU8KVmEQgIu+KyG4RWRHoWHIjIrVFZIaIrBKRlSLyl0DHlBMRCReRX0VkqRfvE4GOKTciEiwii0Xkq0DHkhsRiReR5SKyREQWBDqenIhIpIhMEJE1IrJaRC4MdEzZEZGm3md68nZIRP4a6LiyIyL3ef9fK0TkExEJ99u5S0obgYh0BQ4DH6pqy0DHkxMRqQ5UV9VFIlIOWAhcraqrAhxalkREgDKqelhEQoE5wF9UdV6AQ8uWiNwPxAHlVfXKQMeTExGJB+JUtdAPehKRD4CfVPVtESkFlFbVxACHlSsRCQa2Ax1V9VwHq+YbEamJ+79qrqrHRORT4BtVfd8f5y8xJQJVnQ3sD3QcvlDVnaq6yLufBKwGagY2quypc9h7GOrdCu0VhojUAq4A3g50LMWJiFQAugLvAKhqclFIAp5ewMbCmAQyCQEiRCQEKA3s8NeJS0wiKKpEpB7QFpgf4FBy5FW1LAF2Az+qamGO9yXgASA9l+MKCwV+EJGFInJ7oIPJQX1gD/CeV+32toiUCXRQProO+CTQQWRHVbcDzwNbgZ3AQVX9wV/nt0RQiIlIWWAi8FdVPRToeHKiqmmqGgvUAjqISKGsfhORK4Hdqrow0LHkwcWq2g64DLjLq+YsjEKAdsDrqtoWOAI8FNiQcudVYfUHPgt0LNkRkYrAVbhkWwMoIyI3+Ov8lggKKa+ufSIwVlUnBToeX3lVATOASwMcSnY6A/29evdxQE8R+SiwIeXMuxpEVXcDk4EOgY0oW9uAbZlKgxNwiaGwuwxYpKoJgQ4kB72Bzaq6R1VTgEnARf46uSWCQshrfH0HWK2qLwQ6ntyISJSIRHr3I4BLgDUBDSobqvqwqtZS1Xq46oDpquq3Kyt/E5EyXocBvGqWPkCh7PmmqruA30WkqbepF1AoOzicYSiFuFrIsxXoJCKlve+HXri2Q78oMYlARD4BfgGaisg2Ebkl0DHloDNwI+5q9WTXtssDHVQOqgMzRGQZ8BuujaDQd8ssIqKBOSKyFPgV+FpVvwtwTDm5Bxjr/S3EAv8KbDg585LrJbgr7ELLK2VNABYBy3Hf3X6bbqLEdB81xhiTtRJTIjDGGJM1SwTGGFPCWSIwxpgSzhKBMcaUcJYIjDGmhLNEUIKJyEwReTVAr60iMjAQr+2Lwh6fP3izWI7y4bgZIvLHAggpq9fO8ffgzXp7bUHGVBxZIiimvEFe//WmMD4hIgkiMk1ELsl02ADg4UDF6G8i0s774uiSzf7xIjK3oOPKSXZfdCLyfmGYIltErgBqA2MzbYv34lYROeZNOT3SG+hU0J4ERouIfZedB/vwiq+JuKkIbgGaAFcC3wKVTx6gqvu92U2LHBEJOfOLx5uxdQlwcxbHVwauxmYczau/AO+ratoZ2/+BG0jYDDcZ2r+AQEyI9w1QDjdNhDlHlgiKIW+6hy7AQ6o6TVW3qOpvqvq8qo7LdNxpVUPeld7fROR/3iId20Rk5BnnbiIis0TkuIisFZHLReSwiAz39tfzrhTjznhebkX80d75jnlxPJt54Q0RGeVVZQwXkY3ACSCrmS3fBgZ5E/ZldoP3nPEicqmI/CQiB0Rkv4h8LyLNcojNp/ckIjVFZJx33gMi8rWINM7uvHkhIgNEZJn3+ez3fgfRmfb3Ezc76XER2SwiT4mbTO3k/qoi8rn3/C0iclayzOI1o3Bz3HyZxe4kVd2lqvGq+jawDDf9xcnnNvReb5eIHBGRReIm/Mt8/lz/3rKI6UER2SsincBNdohLBkNzez8me5YIiqfD3q2/5H0Vo/twQ9jbAc8Az4q3ypRX/J4MpAKdgOHA40CYH2I+gruSbwb8GTcP0KNnHFMf+AMwCGgDHM/iPGOBYGDIGdtvAcar6hFcAnkJV2LqDhwEvsz8xZlXIlIaN9necaAbcCFuuuCp3r5zJiLVcBPkfYD7fLoCYzLt74t7368CLXCf40BOn97hfaAR7ov9auCPQL1cXvpiXPLMdm4jcbp7caVk2lUWVwK9BPe7mghMEpGYM06R7d9bFq/zPG4Ki25nLHr0K+4zN+dKVe1WDG/AtbiFeI7j5lh6Hrf6UuZjZgKvZnocD3xyxjHrgb959/vikkDNTPsvws2XP9x7XM97HHfGeRQYmN3jLOK/A9iQ6fEo3BdNtA/v/SNgbqbHF3iv1zGb48sAabjpns+Kz5f3hPvyXY83bYu3LRjYBwzOIdYsPwfcF/dX3v123nF1sznHbODvZ2y7GncxILiqQQU6Z9pf13vPo3KI7a/Aliy2x+MSxGEg2Tv3MeCiXH4v807+Lfny95bp8xkCvAesy+ozwE0hnQ6EFPT/WXG5WYmgmFLVibh5y/vhrswuAuaJyCO5PHXZGY93AFW9+zHADvWmRfb8hh8WeBGRgSIyx6tKOAy8CNQ547Bt6ttUwW8DF2a6+rwZWKHe9MhetcXHIrJRRA4BCbjS8ZmvlxftcSWWJK+q7DCupFERaHge5wVYCkwFVojIRBG506u2yfzaj558Xe+1P8YluGq4q/V03JUzAOpW4spthasIsi51AbyAm1SuG64k9ISqZjTEi5s19Vlx624f8GKK4+zPOKe/t5Oex5XcLtasVxA7hkt4flvDt6SxRFCMqepxVf1RVf+hqhfhprYelUsVSMoZj5W8/Z2cTAoZDbni1lbIllffOw74Hpe42gJ/wy15mdkRH2OYBWwAbhY3LfZQvOUTPV8BUcCfgI7e66UC2X0uvrynIFxDdewZtybA/3KINQmokMX2SFwiQV09eB/vtgxXzbVeRNpkeu0nznjd1kBj3IphJ+V1hsm9uESWlX2qukFVf8GVPkeISI9M+5/HVeH9HZcsYnGJ6MzP2Je/tx9xCS27GXgrAcf11HKpJo9CAh2AKVCrcL/zcFyRPq/WADVEpIaqnryajOP0f9yTXzzVM22LzeW8nYHtqvrkyQ0iUvcc4gPcGsoi8i6ux8sa3JXtGO+8lXElmz+r6gxvWzty/l/w5T0twiWcvZq3dXrX4q7oMxKVuIXU2+CqQzLeE66K7xcR+QewEldlstR77RhV3ZDVC4jIGtzvqAMw19tWB1dizMliIEpEqqjq3uwOUtUD4jodvCgibb1YLwY+9EqmeG1VDXHVO3n1DW6a6M9ERFX1gzP2t8R9BuYcWYmgGBKRyiIyXURuEJHWIlJfRAbh1umdpue+7OWPuC+uD0SkjXcl/wLuatpV6Koew9UFPygiLUTkItzVYU7WATVF5HoRaSAid3L+vUA+AKp4rz1FVfd52w/grnRvE5FGItINeMN7D1ny8T2NxVUxfS4i3bzPvKuI/DuXnkMv4Eoud4nrkRWLm2e+kvcTEenk9a65wPsC74/r239y0Zd/AH8QkX+ISEsRifGq2p714l8LfAf8T0Qu9F7jfVyVSk4W49agvjiX4wD+CzTFlQLA/U6vETe2oxWu3eacq27UrW8xCHhDzh7c1gX3/sw5skRQPB3GfXH9BVdNshLXg+Rjzu5N4zNVTQeuwfUS+hX3ZfsULglkrks+2TXxN1y1yN9yOe+XwHO4njzLcD1NHjvXOL1z7sBdSVYk09gB7z0MwVWdrABew1VfnMjllDm+J1U9iuvNswm39u0a3OdTEZd8sovzE+Am77YA94VWDeiibsUvcFVEnXFVWuuBfwNPqupH3jm+B64AeuB+L7/i1gremumlhgObgem47qAf4xprs+VVSb0LXJ/Tcd6xu3GlrlFe77L7cUnkJ1wb1Tzv/jnzksFgXEL7I7guu7j2r/dyeq7JmS1MY86LV0+9BNejpigtCG98ICJVcSWPC1R1c6DjOZOIPAdUUNVADGYrNqyNwOSJiFyDa7Rdj+tW+QKn6qlNMaOqu8UNPquDK1EUNrvJverR5MJKBCZPvCL533B11AdwYxHu87FbpzGmELJEYIwxJZw1FhtjTAlnicAYY0o4SwTGGFPCWSIwxpgSzhKBMcaUcP8P6f6xra7TpV4AAAAASUVORK5CYII=\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",
    "# 计算 Frobenius-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",
    "# 画图 \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(list(range(1, T+1)), err_subfull, \"o-.\", \n",
    "        label = 'Reconstruction via VQSVD')\n",
    "ax.plot(list(range(1, T+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": 11,
   "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": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.337792Z",
     "start_time": "2021-03-09T03:50:17.231211Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsUlEQVR4nO3dbWxcVXoH8P8fx46NHSW2A4mdOC+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=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUkklEQVR4nO3df4xV5Z3H8fe3yK/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==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASZUlEQVR4nO3df4xdZZ3H8feH2h9AK6XbsU5K6ZS2AmWBgYxdQQS2VoP800o2YmOUjWRrNpBogskSN1nZxGR1s2rchLgpP2I1rsiKRlzrQiEutEJopy5Oi0VoaRFK2xki/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=\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": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.966739Z",
     "start_time": "2021-03-09T03:50:18.950606Z"
    }
   },
   "outputs": [],
   "source": [
    "# 然后我们再来看看量子版本的分解效果：\n",
    "\n",
    "# 超参数设置\n",
    "N = 5           # 量子比特数量\n",
    "T = 8           # 设置想要学习的阶数\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)\n",
    "\n",
    "# 设置电路参数\n",
    "cir_depth = 40                           # 电路深度\n",
    "block_len = 1                            # 每个模组的长度\n",
    "theta_size = N * block_len * cir_depth   # 网络参数 theta 的大小"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:19.553293Z",
     "start_time": "2021-03-09T03:50:19.531180Z"
    }
   },
   "outputs": [],
   "source": [
    "# 定义量子神经网络\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(1, N):\n",
    "            cir.cnot([which_qubit - 1, which_qubit])\n",
    "\n",
    "    return cir.U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:21.636509Z",
     "start_time": "2021-03-09T03:50:21.593577Z"
    }
   },
   "outputs": [],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    \n",
    "    # 初始化可学习参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "    def __init__(self, shape, dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # 创建用来学习 U 的参数 theta\n",
    "        self.theta = self.create_parameter(shape=shape,\n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的参数 phi\n",
    "        self.phi = self.create_parameter(shape=shape,\n",
    "                                         default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                         dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # 将 Numpy array 转换成 Paddle 支持的 Tensor\n",
    "        self.M = paddle.to_tensor(Mat_generator())\n",
    "        self.weight = paddle.to_tensor(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 -= paddle.real(self.weight)[i] * paddle.real(matmul(U_dagger,matmul(self.M, V)))[i][i]\n",
    "            singular_values[i] = paddle.real(matmul(U_dagger, matmul(self.M, V)))[i][i].numpy()\n",
    "        \n",
    "        # 函数返回两个矩阵 U 和 V_dagger 学习的奇异值以及损失函数    \n",
    "        return U, V_dagger, loss, singular_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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: 11569.3661\n",
      "iter: 10 loss: -84679.3080\n",
      "iter: 20 loss: -103418.2506\n",
      "iter: 30 loss: -120854.8849\n",
      "iter: 40 loss: -130997.4657\n",
      "iter: 50 loss: -140376.4574\n",
      "iter: 60 loss: -146738.7417\n",
      "iter: 70 loss: -149536.6623\n",
      "iter: 80 loss: -151092.6062\n",
      "iter: 90 loss: -152296.0673\n",
      "iter: 100 loss: -153290.4326\n",
      "iter: 110 loss: -154112.8803\n",
      "iter: 120 loss: -154812.4464\n",
      "iter: 130 loss: -155384.0762\n",
      "iter: 140 loss: -155814.4174\n",
      "iter: 150 loss: -156132.5890\n",
      "iter: 160 loss: -156377.1011\n",
      "iter: 170 loss: -156571.0730\n",
      "iter: 180 loss: -156728.6039\n",
      "iter: 190 loss: -156859.1820\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x13e2c4410>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYG0lEQVR4nO2dW4yd1XXH/8vG1/H4Mvg6Y7CxwTLmZmAgIYkS6igRIBRAqhA8IB5QHFVBKlL6gKhUqNQHUhVQHqpUTkEhVRpCkyBbFWpDUSInLyaD8Q0bigO2mfF4xnh8xxjbs/pwPqtj9K3/nNlnzncM+/+TRnNmr7O/vc4+e805Z//PWtvcHUKILz4TWu2AEKIaFOxCZIKCXYhMULALkQkKdiEyQcEuRCZc0khnM7sDwI8ATATwr+7+NLt/R0eHd3V1ldqYBGhmY+4zYUL8fyxVboz8mDhxYpIf586dC22ffvpp/Y7VAfMxlfGWbdl8MP8vuSRexpGPzPfUxzXe6yolJvr6+jA0NFRqTA52M5sI4J8BfAtAL4A/mdkGd98Z9enq6sL69etLbWxxT548ecx9pk2bFtrOnj0b2tgET506tbS9ra0t7NPe3h7ahoaGQtu+fftCGyNa+MxHFkjDw8Oh7cyZM/U7VsDm9+TJk6FtxowZoW3u3LmhLfoHwnxn64PBrskedzT/7HpRTNx7771hn0bext8KYLe7v+/unwJ4CcA9DVxPCNFEGgn2LgAfjvi7t2gTQlyENH2DzszWmlmPmfWwt61CiObSSLD3AbhsxN+Li7YLcPd17t7t7t0dHR0NDCeEaIRGgv1PAK4ysyvMbDKABwBsGB+3hBDjTfJuvLufNbNHAfw3atLbC+7+NutjZuHOI5OoIpmBSS5Rn0aIfGQ7rUxOSlUFUiSqVFmL7cazOY5sqddj/jNbivTGSO3H1ncKKeu7IZ3d3V8F8Goj1xBCVIO+QSdEJijYhcgEBbsQmaBgFyITFOxCZEJDu/EpRBIEkxIiuYNJb0ziYVIN65cin6SOFSU6ADwBKErWYbC5Z/JgSmZhSkLIaGOxfmz+U/ow/9nzmSpvpvgRoVd2ITJBwS5EJijYhcgEBbsQmaBgFyITKt+Nj2C7rdHOdGpSRerObmRjO+fMNmXKlNA2ffr00JayW8z6sJJPbD7YTv0nn3xS2s6UBPZ8svmYNGlSaIsUG+Y7e8yptetSduNTrkfXfWgRQnyhULALkQkKdiEyQcEuRCYo2IXIBAW7EJlQufQWyQns9ItIPmlGXTVGNB6TcVJPVPn444+TrhnBJMDUenfMFj1nTPaM5DoAOHXqVGhjp/9Ejy31OUtNhEmpN8hkymh+6ZFooUUI8YVCwS5EJijYhcgEBbsQmaBgFyITFOxCZEJD0puZ7QFwHMA5AGfdvZvd391DCYJJExFMmmhrawttKUdNAbGPrE9q1huTUFiWVzQnzMcTJ06ENuYjk5oiH9ncM3mNrY+UOn9sflNrG6baojlJqZ/HGA+d/S/c/aNxuI4QoonobbwQmdBosDuA35rZm2a2djwcEkI0h0bfxn/N3fvMbD6A18zsHXffOPIOxT+BtQDQ2dnZ4HBCiFQaemV3977i9yCAVwDcWnKfde7e7e7dHR0djQwnhGiA5GA3szYzaz9/G8C3AewYL8eEEONLI2/jFwB4pZB0LgHw7+7+X6N1iiQgJoVE2VCpUgcbi0k8kZyUKpEwP06fPh3a5syZE9qibD9WVPLgwYOh7ejRo6GNSYAzZ84sbWdSJJMAmSzHZMXouWESYGpWZGrWWzQnKTIwlRRDyyi4+/sAbkjtL4SoFklvQmSCgl2ITFCwC5EJCnYhMkHBLkQmXDRnvTH5KspgY0UZWcFJBpOGjhw5Utre19cX9kkpytiI7fjx46Xt+/btC/v09/eHNjYfBw4cCG1RthyTFLdv3x7aBgYGQtt11103ZtuSJUvCPgsXLgxtTF5ja47JlNE12fpub28vbddZb0IIBbsQuaBgFyITFOxCZIKCXYhMqHQ33szCXUmW3BER7UgC/CghtmvKdroPHTpU2r5p06awz86dO0Pb4sWLQ1tqOnBUg47tnH/0UVxVLFIgAL7DH+3Gs4ScwcHB0PbOO++ENraL/53vfKe0ffbs2WGfSy+9NLSxRJOpU6eGtunTp4e2kydPlrazmEhJntEruxCZoGAXIhMU7EJkgoJdiExQsAuRCQp2ITKh8kSYSLpg0sTQ0FBp+6xZs8I+rHZa6tFQUY203bt3h302btwY2ljCxYcffhja9u/fH9rmzp075rGYHJYq2UXjsXLi8+fPD23sMTNblOTDEnJYAgqrJTdt2rTQxuTeKBGGSXkp6JVdiExQsAuRCQp2ITJBwS5EJijYhcgEBbsQmTCq9GZmLwC4G8Cgu19btHUA+CWApQD2ALjf3Q+Pdq0JEyaE8gQ7+ieSLVgGEpOTmLzGjhlaunRpafuaNWvCPsuWLQtthw/HU/b73/8+6Zq33HJLafvdd98d9mHZg0zCZPX1oswx9piZTMmy3pYvXx7avvSlL5W2MymSrQ9Wk4/Vp2PrOxrv2LFjYZ/oeC1GPa/sPwVwx2faHgfwurtfBeD14m8hxEXMqMFenLf+2W+13APgxeL2iwDuHV+3hBDjTepn9gXufv6rSQdQO9FVCHER0/AGndc+OIcfns1srZn1mFlPVOlFCNF8UoN9wMwWAUDxO6wn5O7r3L3b3btZuR8hRHNJDfYNAB4ubj8MYP34uCOEaBb1SG+/AHA7gLlm1gvgSQBPA3jZzB4BsBfA/fUMNjw8HBaCZIUNo+KLTJpgBRtTj42KivmtXLky7LNq1arQxjKo7rzzztB2+eWXh7ZIjmQZVEwOu+aaa0Ibk94iP/bs2RP2+cMf/hDa2EfA22+/PbR1dXWVtjMJjcm2rCBpVOwT4HMVydEsM4/JfBGjBru7PxiYvjnm0YQQLUPfoBMiExTsQmSCgl2ITFCwC5EJCnYhMqHSgpPuHsoJTJqIZAZ2nhs984pkNbF+kVTGZBxW3JLJfKyY5pIlS0JbVJzzzJkzYR8m47Aiikw6jGRRJvMdPHgwtDHYHEf+p0isAMKzCgG+Htn6jtZjSpFKlgmqV3YhMkHBLkQmKNiFyAQFuxCZoGAXIhMU7EJkQuXSWyRBsEyjqA+TJpgMwmAyVDQek1VYAUsma7H5YNl+kS8s0y+1zgCb4507d5a2v/XWW2EfluXFMgtZFuC8efNK21kWIMteY9Ibk72mTJkS2iLpjT0vTEoNxxlzDyHE5xIFuxCZoGAXIhMU7EJkgoJdiEyofDc+2umcPXt22C/afWa7n2xHlSWusCSICFaXjCWLsIQctuvL6vVF12QqQ7RjDXA1gakCW7ZsKW3funVr2IfN41e+8pXQxurkdXZ2lrazNRAlEwFptd8A/tii52zBgvg4hn379o3ZB72yC5EJCnYhMkHBLkQmKNiFyAQFuxCZoGAXIhPqOf7pBQB3Axh092uLtqcAfBfA+aJhT7j7qw05QqSJSCZh0htLFEitQZeSCMMSUNhjZnJYX19faIuSa5ictHz58tDGpLfe3t7QtmPHjtJ2JhkxqYn5OH369NAWJaCk1uRjki7rx6RUds2IFIm4nlf2nwK4o6T9OXdfXfw0FOhCiOYzarC7+0YA8cuCEOJzQSOf2R81s21m9oKZzRk3j4QQTSE12H8MYDmA1QD6ATwT3dHM1ppZj5n1sJrhQojmkhTs7j7g7ufcfRjATwDcSu67zt273b17zhy9ARCiVSQFu5ktGvHnfQDKt16FEBcN9UhvvwBwO4C5ZtYL4EkAt5vZagAOYA+A79UzmJklZWVFEhWTLFLlNVYrLJKh2Fis1hmTDplUxogknoULF4Z92Nzv378/tG3evDm0RfJgR0dH2Ofmm28ObSyzjb1jjOaYzX1qViSTUhlRZiSr8cf8jxjVO3d/sKT5+TGPJIRoKfoGnRCZoGAXIhMU7EJkgoJdiExQsAuRCZUWnDSzUIpi8k9UjJJJEydPngxtTLJjxy5Fcge7His4yWypR/9E12TfXmxrawttrN+ePXtCW1QUk0lvq1atCm2LFi0KbSdOnAht0Vyx9cGkt9SjoVIKiLK5j3xk0qBe2YXIBAW7EJmgYBciExTsQmSCgl2ITFCwC5EJlUpvQFqhPCatRLCsIOYDKx4ZyThMCmOkZr2lSDzt7e1hnwMHDoS2qHAkwItHRhlgl19+edjniiuuCG1MHkwpEBkVDwVGka8SsylT1ggbK3pcbE3plV2ITFCwC5EJCnYhMkHBLkQmKNiFyITKd+MjWP2ugwcPlrZHRx0BfIeW1Zlju/FR4g1LjmC7o+wxHzt2LLSxZJ1o1/f06dNhH7bj/sYbb4Q2phjMmjWrtH3p0qVhH1YnjykQUaIUAAwODpa2s7XDEpRYkkzqcWTRNdnzPDAwMGYf9MouRCYo2IXIBAW7EJmgYBciExTsQmSCgl2ITKjn+KfLAPwMwALUjnta5+4/MrMOAL8EsBS1I6Dud3d6TOvEiRNDOeH48eNhvygJokoZBIjln9S6ZClH+AC8HlsksbGaa++++25oY8kukSQKAN/4xjdK29esWRP2WbFiRWhjx1Cx+Y/kwfnz54d9WJ08JmGyRBgmBUfrsbe3N+wTSYeN1qA7C+AH7r4KwJcBfN/MVgF4HMDr7n4VgNeLv4UQFymjBru797v75uL2cQC7AHQBuAfAi8XdXgRwb5N8FEKMA2P6zG5mSwHcCGATgAXu3l+YDqD2Nl8IcZFSd7Cb2QwAvwbwmLtf8F1Or334LP0AamZrzazHzHoOHTrUkLNCiHTqCnYzm4RaoP/c3X9TNA+Y2aLCvghA6ZeQ3X2du3e7ezc7+EAI0VxGDXarbe89D2CXuz87wrQBwMPF7YcBrB9/94QQ40U9WW9fBfAQgO1mtqVoewLA0wBeNrNHAOwFcP9oF3L3MKOISQZRDTrWJ7XGWEoNOiahMemNZd8x/1ltskgOY5lt27ZtC20ss62zszO0rVy5srSd1aBjWYBMLmUZbNE8srFo5hiZe+YjW1fRGmHrNEW2HTXY3f2PAKJRvznmEYUQLUHfoBMiExTsQmSCgl2ITFCwC5EJCnYhMqHSgpPDw8Nh0UYmaZw6daq0nWWopUorrNhgZEvJaAJ4oUT22Ng3Ebdu3Vravn59/DWIzZs3h7apU6eGtmuvvTa0RRlsUSFKgM89g81VJG8yWYtl0bF+0doejUg6ZNJs0nFSY+4hhPhcomAXIhMU7EJkgoJdiExQsAuRCQp2ITKhUunN3UNZg0llUeYSk0iYHMMykJgEGMEkkqhYJgBMnz49tDFphUk8UYZg1A7wuVq2bFlou/LKK0PbggXlhYtYNt/Ro0dDG3uumfQZrSsmvzJ5LTUTjfkYXXO816le2YXIBAW7EJmgYBciExTsQmSCgl2ITKh0Nx6IdxjZbnyUPJG6+8l2s1k9s8jGdpgZzEeWJLNly5bQtnfv3tJ2VoNu8eLFoe22224LbTfeeGNoi+bqyJEjYR+mQLAjmRjt7e2l7WwNsLWYUqMQ4OsqUhpYYlD0uBo9/kkI8QVAwS5EJijYhcgEBbsQmaBgFyITFOxCZMKo0puZXQbgZ6gdyewA1rn7j8zsKQDfBXD+vKEn3P1Vdq3h4WGcPn261MZqnaVIEymJNUCanMcSSRhM/hkYGAht27dvD2379+8vbWdS3vXXXx/aWJ05liQTPWfsMbPkDnYoKJO8Umq1pfQZjZRkHSbXseuF49Rxn7MAfuDum82sHcCbZvZaYXvO3f9pzKMKISqnnrPe+gH0F7ePm9kuAF3NdkwIMb6M6TO7mS0FcCOATUXTo2a2zcxeMLM54+2cEGL8qDvYzWwGgF8DeMzdjwH4MYDlAFaj9sr/TNBvrZn1mFnP4cOHG/dYCJFEXcFuZpNQC/Sfu/tvAMDdB9z9nLsPA/gJgFvL+rr7OnfvdvfuOXP04i9Eqxg12K32zfrnAexy92dHtC8acbf7AMSZFkKIllPPbvxXATwEYLuZbSnangDwoJmtRk2O2wPge/UMyKStiEiCOHbsWNgnpUYXwGW0SO5IzZI6efJkaPvggw9C23vvvRfahoaGStsXLlwY9rnhhhtC28033xzaOjs7Q1t/f39pO8t6YxlbTGr6+OOPQ1v0fNLsMLJ2mB9M0mW2yBcWKynru57d+D8CKPOGaupCiIsLfYNOiExQsAuRCQp2ITJBwS5EJijYhciESgtOTpgwITzyiEkJKUUqmbSSmvUW+ciOJmIZVCyzbePGjaHtwIEDoS3KDmMS2urVq0Mbkw4jeQ2IJSqWvcYyH1mGY8r8szXAxmLrg2WppR43FZGS9aZXdiEyQcEuRCYo2IXIBAW7EJmgYBciExTsQmRCpdKbmWHKlCnljiRkjjGp49SpU9SPCCa7TJo0qbSdZa+xLC8mvTFZi2XErVixorSdZfPNmzcvtEWPeTRbJIcxCYrJYVGhUoBnvaWsNyblMSmSXTNlfadkyjH0yi5EJijYhcgEBbsQmaBgFyITFOxCZIKCXYhMqFR6A2LphckMY70WkF4YMKXYIJPrmGTEJLve3t7QxmS5qLBkW1tb2IeV+GaPjUleJ06cKG2PpDAAmDlzZmhj/dg6SFlvTNZiNiavpawrJm0q600IEaJgFyITFOxCZIKCXYhMULALkQmj7sab2VQAGwFMKe7/K3d/0syuAPASgEsBvAngIXePMwhq1wp3JdkuZ7RrzWqWMdjOKKsZF8F23FniBOvHkmTYrm9U4+3qq68O+7Cd+r1794Y2tlMf7Z6n1mljCUWsX5TUMmPGjLAPe84YqQlWkTLAkpeOHj1a2k5rKIaW/+c0gDXufgNqxzPfYWZfBvBDAM+5+5UADgN4pI5rCSFaxKjB7jXOi6aTih8HsAbAr4r2FwHc2wwHhRDjQ73ns08sTnAdBPAagD8DOOLu59+b9ALoaoqHQohxoa5gd/dz7r4awGIAtwJYWe8AZrbWzHrMrOfQoUNpXgohGmZMu/HufgTA7wDcBmC2mZ3fKVoMoC/os87du929mx0QIIRoLqMGu5nNM7PZxe1pAL4FYBdqQf+Xxd0eBrC+ST4KIcaBehJhFgF40cwmovbP4WV3/08z2wngJTP7BwBvAXi+ngEj2YtJXpHExiQolijAbEyWiyQSJl0xZs+eHdqWL18e2pgM1dnZWdrOkiqYZBQd1wXwRJ7ommwsJnmlHp8UPWepteTY2mE+MiJf2PyydRoxarC7+zYAN5a0v4/a53chxOcAfYNOiExQsAuRCQp2ITJBwS5EJijYhcgES5ULkgYzOwjgfBrVXAAfVTZ4jPy4EPlxIZ83P5a4e+l5XpUG+wUDm/W4e3dLBpcf8iNDP/Q2XohMULALkQmtDPZ1LRx7JPLjQuTHhXxh/GjZZ3YhRLXobbwQmdCSYDezO8zsXTPbbWaPt8KHwo89ZrbdzLaYWU+F475gZoNmtmNEW4eZvWZm7xW/4zOZmuvHU2bWV8zJFjO7qwI/LjOz35nZTjN728z+umivdE6IH5XOiZlNNbM3zGxr4cffF+1XmNmmIm5+aWZxRcoy3L3SHwATUStrtQzAZABbAayq2o/Clz0A5rZg3K8DuAnAjhFt/wjg8eL24wB+2CI/ngLwNxXPxyIANxW32wH8L4BVVc8J8aPSOQFgAGYUtycB2ATgywBeBvBA0f4vAP5qLNdtxSv7rQB2u/v7Xis9/RKAe1rgR8tw940Ahj7TfA9qhTuBigp4Bn5Ujrv3u/vm4vZx1IqjdKHiOSF+VIrXGPcir60I9i4AH474u5XFKh3Ab83sTTNb2yIfzrPA3c8fz3oAwIIW+vKomW0r3uY3/ePESMxsKWr1EzahhXPyGT+AiuekGUVec9+g+5q73wTgTgDfN7Ovt9ohoPafHbV/RK3gxwCWo3ZGQD+AZ6oa2MxmAPg1gMfc/dhIW5VzUuJH5XPiDRR5jWhFsPcBuGzE32Gxymbj7n3F70EAr6C1lXcGzGwRABS/B1vhhLsPFAttGMBPUNGcmNkk1ALs5+7+m6K58jkp86NVc1KMfQRjLPIa0Ypg/xOAq4qdxckAHgCwoWonzKzNzNrP3wbwbQA7eK+msgG1wp1ACwt4ng+ugvtQwZxYrYjc8wB2ufuzI0yVzknkR9Vz0rQir1XtMH5mt/Eu1HY6/wzgb1vkwzLUlICtAN6u0g8Av0Dt7eAZ1D57PYLamXmvA3gPwP8A6GiRH/8GYDuAbagF26IK/Pgaam/RtwHYUvzcVfWcED8qnRMA16NWxHUbav9Y/m7Emn0DwG4A/wFgyliuq2/QCZEJuW/QCZENCnYhMkHBLkQmKNiFyAQFuxCZoGAXIhMU7EJkgoJdiEz4P3ftN5tVvPS1AAAAAElFTkSuQmCC\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",
    "net = NET([theta_size])\n",
    "\n",
    "# 一般来说，我们利用 Adam 优化器来获得相对好的收敛\n",
    "# 当然你可以改成 SGD 或者是 RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(ITR):\n",
    "\n",
    "    #  前向传播计算损失函数\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # 反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 记录优化中间结果\n",
    "    loss_list.append(loss[0][0].numpy())\n",
    "    singular_value_list.append(singular_values)\n",
    "    \n",
    "    if itr% 10 == 0:\n",
    "        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = U.numpy()\n",
    "V_dagger_learned = V_dagger.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')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\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.10"
  },
  "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
}