{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Singular Value Decomposition\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this tutorial, we will go through the concept of classical singular value decomposition (SVD) and the quantum neural network (QNN) version of variational quantum singular value decomposition (VQSVD) [1]. The tutorial consists of the following two parts: \n",
    "- Decompose a randomly generated $8\\times8$ complex matrix; \n",
    "- Apply SVD on image compression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "Singular value decomposition (SVD) has many applications, including principal component analysis (PCA), solving linear equations and recommender systems. The main task is formulated as following:\n",
    "> Given a complex matrix $M \\in \\mathbb{C}^{m \\times n}$, find the decomposition in form $M = UDV^\\dagger$, where $U_{m\\times m}$ and $V^\\dagger_{n\\times n}$ are unitary matrices, which satisfy the property $UU^\\dagger = VV^\\dagger = I$.\n",
    "\n",
    "- The column vectors $|u_j\\rangle$ of the unitary matrix $U$ are called left singular vectors $\\{|u_j\\rangle\\}_{j=1}^{m}$ form an orthonormal basis. These column vectors are the eigenvectors of the matrix $MM^\\dagger$.\n",
    "- Similarly, the column vectors $\\{|v_j\\rangle\\}_{j=1}^{n}$ of the unitary matrix $V$ are the eigenvectors of $M^\\dagger M$ and form an orthonormal basis.\n",
    "- The diagonal elements of the matrix $D_{m\\times n}$ are singular values $d_j$ arranged in a descending order.\n",
    "\n",
    "For the convenience, we assume that the $M$ appearing below are all square matrices. Let's first look at an example: \n",
    "\n",
    "$$\n",
    "M = 2*X\\otimes Z + 6*Z\\otimes X + 3*I\\otimes I = \n",
    "\\begin{bmatrix} \n",
    "3 &6 &2 &0 \\\\\n",
    "6 &3 &0 &-2 \\\\\n",
    "2 &0 &3 &-6 \\\\\n",
    "0 &-2 &-6 &3 \n",
    "\\end{bmatrix}, \\tag{1}\n",
    "$$\n",
    "\n",
    "Then the singular value decomposition of the matrix can be expressed as:\n",
    "\n",
    "$$\n",
    "M = UDV^\\dagger = \n",
    "\\frac{1}{2}\n",
    "\\begin{bmatrix} \n",
    "-1 &-1 &1 &1 \\\\\n",
    "-1 &-1 &-1 &-1 \\\\\n",
    "-1 &1 &-1 &1 \\\\\n",
    "1 &-1 &-1 &1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} \n",
    "11 &0 &0 &0 \\\\\n",
    "0 &7 &0 &0 \\\\\n",
    "0 &0 &5 &0 \\\\\n",
    "0 &0 &0 &1 \n",
    "\\end{bmatrix}\n",
    "\\frac{1}{2}\n",
    "\\begin{bmatrix} \n",
    "-1 &-1 &-1 &-1 \\\\\n",
    "-1 &-1 &1 &1 \\\\\n",
    "-1 &1 &1 &-1 \\\\\n",
    "1 &-1 &1 &-1 \n",
    "\\end{bmatrix}. \\tag{2}\n",
    "$$\n",
    "\n",
    "Import packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.008567Z",
     "start_time": "2021-03-09T03:44:29.796997Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import pi as PI\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import unitary_group\n",
    "from scipy.linalg import norm\n",
    "\n",
    "import paddle\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.linalg import dagger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Draw the learning curve in the optimization process\n",
    "def loss_plot(loss):\n",
    "    '''\n",
    "    loss is a list, this function plots loss over iteration\n",
    "    '''\n",
    "    plt.plot(list(range(1, len(loss)+1)), loss)\n",
    "    plt.xlabel('iteration')\n",
    "    plt.ylabel('loss')\n",
    "    plt.title('Loss Over Iteration')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical Singular Value Decomposition\n",
    "\n",
    "With the above mathematical definition, one can realize SVD numerically through NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.056721Z",
     "start_time": "2021-03-09T03:44:34.012222Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The matrix M we want to decompose is: \n",
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# Generate matrix M\n",
    "def M_generator():\n",
    "    I = np.array([[1, 0], [0, 1]])\n",
    "    Z = np.array([[1, 0], [0, -1]])\n",
    "    X = np.array([[0, 1], [1, 0]])\n",
    "    Y = np.array([[0, -1j], [1j, 0]])\n",
    "    M = 2 *np.kron(X, Z) + 6 * np.kron(Z, X) + 3 * np.kron(I, I)\n",
    "    return M.astype('complex64')\n",
    "\n",
    "print('The matrix M we want to decompose is: ')\n",
    "print(M_generator())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.093725Z",
     "start_time": "2021-03-09T03:44:34.063353Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The singular values of the matrix from large to small are:\n",
      "[11.  7.  5.  1.]\n",
      "The decomposed unitary matrix U is:\n",
      "[[-0.5+0.j -0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j -0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [ 0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]]\n",
      "The decomposed unitary matrix V_dagger is:\n",
      "[[-0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j  0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j -0.5+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# We only need the following line of code to complete SVD\n",
    "U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)\n",
    "\n",
    "\n",
    "# Print decomposition results\n",
    "print(\"The singular values of the matrix from large to small are:\")\n",
    "print(D)\n",
    "print(\"The decomposed unitary matrix U is:\")\n",
    "print(U)\n",
    "print(\"The decomposed unitary matrix V_dagger is:\")\n",
    "print(V_dagger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.112670Z",
     "start_time": "2021-03-09T03:44:34.098847Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# Then assemble it back, can we restore the original matrix?\n",
    "M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))\n",
    "print(M_reconst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Surely, we can be restored the original matrix $M$! One can further modify the matrix, see what happens if it is not a square matrix.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Singular Value Decomposition\n",
    "\n",
    "Next, let's take a look at what the quantum version of singular value decomposition is all about. In summary, we transform the problem of matrix factorization into an optimization problem with the variational principle of singular values. Specifically, this is achieved through the following four steps:\n",
    "\n",
    "- Prepare an orthonormal basis $\\{|\\psi_j\\rangle\\}$, one can take the computational basis $\\{ |000\\rangle, |001\\rangle,\\cdots |111\\rangle\\}$ (this is in the case of 3 qubits)\n",
    "- Prepare two parameterized quantum neural networks $U(\\theta)$ and $V(\\phi)$ to learn left/right singular vectors respectively\n",
    "- Use quantum neural network to estimate singular values $m_j = \\text{Re}\\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle$\n",
    "- Design the loss function $\\mathcal{L}(\\theta)$ and use PaddlePaddle Deep Learning framework to maximize the following quantity, \n",
    "\n",
    "$$\n",
    "L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle\\psi_j|U(\\theta)^{\\dagger} MV(\\phi)|\\psi_j\\rangle. \\tag{3}\n",
    "$$\n",
    "\n",
    "Where $q_1>\\cdots>q_T>0$ is the adjustable weights (hyperparameter), and $T$ represents the rank we want to learn or the total number of singular values to be learned.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 1: Decompose a randomly generated $8\\times8$ complex matrix\n",
    "\n",
    "Then we look at a specific example, which can better explain the overall process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.132465Z",
     "start_time": "2021-03-09T03:44:34.116446Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The matrix M we want to decompose is:\n",
      "[[6.+1.j 3.+9.j 7.+3.j 4.+7.j 6.+6.j 9.+8.j 2.+7.j 6.+4.j]\n",
      " [7.+1.j 4.+4.j 3.+7.j 7.+9.j 7.+8.j 2.+8.j 5.+0.j 4.+8.j]\n",
      " [1.+6.j 7.+8.j 5.+7.j 1.+0.j 4.+7.j 0.+7.j 9.+2.j 5.+0.j]\n",
      " [8.+7.j 0.+2.j 9.+2.j 2.+0.j 6.+4.j 3.+9.j 8.+6.j 2.+9.j]\n",
      " [4.+8.j 2.+6.j 6.+8.j 4.+7.j 8.+1.j 6.+0.j 1.+6.j 3.+6.j]\n",
      " [8.+7.j 1.+4.j 9.+2.j 8.+7.j 9.+5.j 4.+2.j 1.+0.j 3.+2.j]\n",
      " [6.+4.j 7.+2.j 2.+0.j 0.+4.j 3.+9.j 1.+6.j 7.+6.j 3.+8.j]\n",
      " [1.+9.j 5.+9.j 5.+2.j 9.+6.j 3.+0.j 5.+3.j 1.+3.j 9.+4.j]]\n",
      "The singular values of the matrix M are:\n",
      "[54.83484985 19.18141073 14.98866247 11.61419557 10.15927045  7.60223249\n",
      "  5.81040539  3.30116001]\n"
     ]
    }
   ],
   "source": [
    "# First fix the random seed, in order to reproduce the results at any time\n",
    "np.random.seed(42)\n",
    "\n",
    "# Set the number of qubits, which determines the dimension of the Hilbert space\n",
    "N = 3\n",
    "\n",
    "# Make a random matrix generator\n",
    "def random_M_generator():\n",
    "    M = np.random.randint(10, size = (2**N, 2**N)) + 1j*np.random.randint(10, size = (2**N, 2**N))\n",
    "    return M\n",
    "\n",
    "M = random_M_generator()\n",
    "M_err = np.copy(M)\n",
    "\n",
    "\n",
    "# Output the matrix M\n",
    "print('The matrix M we want to decompose is:')\n",
    "print(M)\n",
    "\n",
    "# Apply SVD and record the exact singular values\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "print(\"The singular values of the matrix M are:\")\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.147570Z",
     "start_time": "2021-03-09T03:44:34.138265Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The selected weight is:\n",
      "[24.+0.j 21.+0.j 18.+0.j 15.+0.j 12.+0.j  9.+0.j  6.+0.j  3.+0.j]\n"
     ]
    }
   ],
   "source": [
    "# Hyperparameter settings\n",
    "N = 3       # Number of qubits\n",
    "T = 8       # Set the number of rank you want to learn\n",
    "ITR = 100   # Number of iterations\n",
    "LR = 0.02   # Learning rate\n",
    "SEED = 14   # Random seed\n",
    "\n",
    "# Set the learning weight \n",
    "weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
    "print('The selected weight is:')\n",
    "print(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We design QNN with the following structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.245692Z",
     "start_time": "2021-03-09T03:44:34.226859Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set circuit parameters\n",
    "cir_depth = 20              # circuit depth\n",
    "block_len = 2               # length of each block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define quantum neural network\n",
    "def U_theta() -> Circuit:\n",
    "\n",
    "    # Initialize the network with Circuit\n",
    "    cir = Circuit(N)\n",
    "    \n",
    "    # Build a hierarchy：\n",
    "    for _ in range(cir_depth):\n",
    "        cir.ry()\n",
    "        cir.rz()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we complete the main part of the algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:46:12.944634Z",
     "start_time": "2021-03-09T03:44:50.626213Z"
    }
   },
   "outputs": [],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    def __init__(self, matrix: np.ndarray, weights: np.ndarray) -> None:\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # Create the parameter theta for learning U\n",
    "        self.cir_U = U_theta()\n",
    "        \n",
    "        # Create a parameter phi to learn V_dagger\n",
    "        self.cir_V = U_theta()\n",
    "        \n",
    "        # Convert Numpy array to Tensor supported in Paddle\n",
    "        self.M = paddle.to_tensor(matrix)\n",
    "        self.weight = paddle.to_tensor(weights)\n",
    "\n",
    "    # Define loss function and forward propagation mechanism\n",
    "    def forward(self):\n",
    "        \n",
    "        # Get the unitary matrix representation of the quantum neural network\n",
    "        U = self.cir_U.unitary_matrix()\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        \n",
    "        V = self.cir_V.unitary_matrix()\n",
    "        V_dagger = dagger(V)\n",
    "        \n",
    "        # Initialize the loss function and singular value memory\n",
    "        loss = 0\n",
    "        singular_values = np.zeros(T)\n",
    "        \n",
    "        # Define loss function\n",
    "        for i in range(T):\n",
    "            loss -= paddle.real(self.weight)[i] * paddle.real(U_dagger @ self.M @ V)[i][i]\n",
    "            singular_values[i] = paddle.real(U_dagger @ self.M @ V)[i][i].numpy()\n",
    "        \n",
    "        # Function returns two matrices U and V_dagger, learned singular values and loss function\n",
    "        return U, V_dagger, loss, singular_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -43.6754\n",
      "iter: 10 loss: -1533.6699\n",
      "iter: 20 loss: -2002.1223\n",
      "iter: 30 loss: -2137.2044\n",
      "iter: 40 loss: -2200.9439\n",
      "iter: 50 loss: -2267.5643\n",
      "iter: 60 loss: -2329.0521\n",
      "iter: 70 loss: -2359.0386\n",
      "iter: 80 loss: -2367.3676\n",
      "iter: 90 loss: -2372.0702\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnyklEQVR4nO3deXxddZ3/8dfn3ux7m6RLkq50gS5QMCBFZBBQiqIFFcUNUGcYZsR9VBBm1PkN88BlRnHcBhVFRJABWWQXBVmEQgqldKE0Ld2TNt2TLtnu5/fHOUlvS1rS3Htzk3vfz8fjPnLv95x7z+fktPed8/2exdwdERGRRETSXYCIiAx/ChMREUmYwkRERBKmMBERkYQpTEREJGEKExERSZjCRESOipm1mdnkdNchQ4vCRIYdM1tjZuekadmnmdlfzKzVzHaZ2R/NbMYgLr933c3sMjN7OsXLe8LM/j6+zd1L3H11Kpcrw4/CRKSfzGwu8ChwL1ADTAJeBp5J9l/qFkjp/08zy0nl50t2UZhIxjCzfDP7gZltCh8/MLP8cFqVmd1vZjvNbLuZPdXzZW1mXzOzjeHexgozO/swi/gO8Bt3v8HdW919u7tfCzwHfDP8rOVmdn5cTTlm1mJmJ4WvTzWzv4V1vGxmZ8bN+4SZXWdmzwB7gcMGlJkdB/wMmBt2O+2M+x18z8zWmdlmM/uZmRWG0840sw3h+jYDvzKzEeHvpcXMdoTP68L5rwPeDvwoXMaPwnY3synh83Iz+034/rVmdm3c7/UyM3s6rGeHmb1uZuf1e4PKsKIwkUxyDXAqMAc4ATgFuDac9mVgA1ANjAa+DriZTQeuBE5291LgXGDNoR9sZkXAacD/9bHcO4B3hs9vAz4SN+1cYKu7v2hmtcADwH8AI4F/Ae4ys+q4+T8BXA6UAmsPt6Luvhy4Ang27HaqCCddD0wLfwdTgFrg3+LeOiZc9oRwORHgV+Hr8cA+4EfhMq4BngKuDJdxZR+l/A9QThB8fwdcAnwybvpbgRVAFUEY/9LM7HDrJcOXwkQyyceAf3f3Le7eAnyL4MsZoBMYC0xw9053f8qDC9N1A/nADDPLdfc17r6qj88eSfD/pamPaU0EX5YAvwPeF4YPwEcJAgbg48CD7v6gu8fc/U9AA/DuuM/6tbsvdfcud+88mpUPv6QvB74Y7jW1Av8JXBw3Wwz4hru3u/s+d9/m7ne5+95w/usIQqE/y4uGn311uKe2BvgvDvzOAda6+8/dvRu4mWAbjD6a9ZLhQWEimaSGg/+aXxu2AXwXaAQeNbPVZnYVgLs3Al8g6KbaYma3m1kNb7SD4It4bB/TxgJb4z5vOfDeMFDeRxAwEPz1f1HYxbUz7Jo6/ZDPXH80K3yIaqAIWBj3+Q+H7T1a3H1/zwszKzKz/w27qHYDTwIVYVC8mSoglzf+zmvjXjf3PHH3veHTkqNYJxkmFCaSSTYRfGH3GB+2Ef7l/GV3n0zwBf+lnrERd/+du58evteBbx/6we6+B3gWuKiP5X4I+HPc656urvnAsjBgIAiKW9y9Iu5R7O7Xxy/qKNb30Hm3EnRTzYz7/HJ3LznCe74MTAfe6u5lwBlhux1m/kOX18kbf+cbj2IdJEMoTGS4yjWzgrhHDsGX+LVmVm1mVQRjBb8FMLPzzWxK2BW0i6B7K2Zm083srHCgfj/Bl3HsMMu8CrjUzD5nZqXh4PV/AHMJutR63A68C/gnDuyVENbyXjM718yiYd1n9gx4D8BmoM7M8gDcPQb8HPi+mY0K17vWzM49wmeUEqzzTjMbCXyjj2X0eSBA2HV1B3Bd+PuYAHwpXE/JMgoTGa4eJPgS7Hl8k2BguwFYDLwCvBi2AUwFHgPaCPYwfuLujxOMl1xP8Fd2MzAKuLqvBbr70wQD6u8nGCdZC5wInO7uK+PmawqXcRrw+7j29QR7K18HWgj2VL7CwP8f/gVYCjSb2daw7WsE3XnPhd1WjxHseRzOD4BCgvV/jqBbLN4NwAfDo7F+2Mf7PwvsAVYDTxOE500DWhsZ1kw3xxIRkURpz0RERBKmMBERkYQpTEREJGEKExERSVjWXuitqqrKJ06cmO4yRESGlYULF2519+pD27M2TCZOnEhDQ0O6yxARGVbMrM9rxqmbS0REEpYxYWJm88LLhzf2XHdJREQGR0aESXhRuh8D5wEzgI/YIN79TkQk22VEmBDct6LR3Ve7ewfBtZHmp7kmEZGskSlhUsvBl+7ewMGXwRYRkRTKlDDpFzO73MwazKyhpaUl3eWIiGSMTAmTjcC4uNd19HFPBXe/0d3r3b2+uvoNh0mLiMgAZUqYvABMNbNJ4b0dLgbuS8WC7nlpI7cuOOytuUVEslJGhIm7dwFXAo8Q3DL1Dndfmopl3b+4iVueVZiIiMTLmDPg3f1BghsmpVRVSR4vb9iZ6sWIiAwrGbFnMpgqS/LYvqeDWEw3FRMR6aEwOUqVxfl0x5xd+zrTXYqIyJChMDlKlSV5AGzb057mSkREhg6FyVGqKskHYGtbR5orEREZOhQmR6l3z0RhIiLSS2FylCqLgz0TdXOJiBygMDlKI4pyMVM3l4hIPIXJUcqJRhhRlMe2Nu2ZiIj0UJgMQGVxnsZMRETiKEwGoLIkT2MmIiJxFCYDUFmSrz0TEZE4CpMBqCrOY6vGTEREeilMBqCqJJ/d+7vo6IqluxQRkSFBYTIAleFZ8Nv3qKtLRAQUJgPScxa8urpERAIKkwGo6r3Yo/ZMRERAYTIgvZdU0Z6JiAigMBkQXexRRORgCpMBKMnPIS8nwladuCgiAihMBsTMqNIlVUREeilMBig4C157JiIioDAZsOD6XNozEREBhcmAVRbr+lwiIj0UJgNUVRJcn8vd012KiEjaKUwGqLIkj/auGHs6utNdiohI2ilMBkgnLoqIHKAwGaAD1+fSuImIiMJkgKpKtGciItJDYTJAlbrYo4hIL4XJAI0s7rk+l/ZMREQUJgOUnxOltCBHYyYiIihMElJVkq9uLhERFCYJqSzOUzeXiAgKk4RUhmfBi4hkuyEXJmb2TTPbaGaLwse746ZdbWaNZrbCzM6Na58XtjWa2VWDVevosgI271aYiIjkpLuAw/i+u38vvsHMZgAXAzOBGuAxM5sWTv4x8E5gA/CCmd3n7stSXWRNRSG79nXS1t5FSf5Q/VWKiKTekNszOYL5wO3u3u7urwONwCnho9HdV7t7B3B7OG/KjS0vAKBp577BWJyIyJA1VMPkSjNbbGY3mdmIsK0WWB83z4aw7XDtb2Bml5tZg5k1tLS0JFxkbUUhABsVJiKS5dISJmb2mJkt6eMxH/gpcAwwB2gC/itZy3X3G9293t3rq6urE/68mjBMNu3cn/BniYgMZ2np6Hf3c/ozn5n9HLg/fLkRGBc3uS5s4wjtKTWqNJ9oxNikPRMRyXJDrpvLzMbGvbwQWBI+vw+42MzyzWwSMBV4HngBmGpmk8wsj2CQ/r7BqDUnGmFMWYHCRESy3lA8BOk7ZjYHcGAN8I8A7r7UzO4AlgFdwGfcvRvAzK4EHgGiwE3uvnSwiq2pKNCYiYhkvSEXJu7+iSNMuw64ro/2B4EHU1nX4dRUFPLiuh3pWLSIyJAx5Lq5hpuaikKad+2nO6Z7wYtI9lKYJKimopDObtdlVUQkqylMElRbEZy4qHETEclmCpMEHTjXRGEiItlLYZIghYmIiMIkYWUFuZTm5+gseBHJagqTJKipKNSeiYhkNYVJEtRUFLBpl8JERLKXwiQJgj0TdXOJSPZSmCRBTUUh2/d0sK+jO92liIikhcIkCWrCc03U1SUi2UphkgQ15To8WESym8IkCXSuiYhkO4VJEowpL8AMNmoQXkSylMIkCXKjEUaX6iZZIpK9FCZJUlOhMBGR7KUwSRKdBS8i2UxhkiTjRxaxYcc+2rt0romIZB+FSZLMrCmnK+a81tyW7lJERAadwiRJZteWA/DKxl1prkREZPApTJJk3MhCygpyFCYikpUUJkliZsyqLWeJwkREspDCJIlm15azormVjq5YuksRERlUCpMkmlVbTkd3jNc2t6a7FBGRQaUwSaKeQXh1dYlItlGYJNGEyiJKNQgvIllIYZJEZsbMmjLtmYhI1lGYJNns2nKWN7fS2a1BeBHJHgqTJJtVW05HV4yVm3UmvIhkD4VJkmkQXkSykcIkySZWFlOSr0F4EckuCpMki0SMGTVlChMRySoKkxSYXVvO8qbd7O/U5ehFJDukJUzM7CIzW2pmMTOrP2Ta1WbWaGYrzOzcuPZ5YVujmV0V1z7JzBaE7b83s7zBXJe+vGP6KNq7Yjy2fHO6SxERGRTp2jNZArwfeDK+0cxmABcDM4F5wE/MLGpmUeDHwHnADOAj4bwA3wa+7+5TgB3ApwdnFQ5v7jGVjCkr4O4XN6a7FBGRQZGWMHH35e6+oo9J84Hb3b3d3V8HGoFTwkeju6929w7gdmC+mRlwFnBn+P6bgQtSvgJvIhoxLjixlidea2FrW3u6yxERSbmhNmZSC6yPe70hbDtceyWw0927Dmnvk5ldbmYNZtbQ0tKS1MIP9f6TaumOOfct2pTS5YiIDAUpCxMze8zMlvTxmJ+qZb4Zd7/R3evdvb66ujqly5o2upTZteX84aUNKV2OiMhQkJOqD3b3cwbwto3AuLjXdWEbh2nfBlSYWU64dxI/f9q9/6RavvXHZaxobmX6mNJ0lyMikjJDrZvrPuBiM8s3s0nAVOB54AVganjkVh7BIP197u7A48AHw/dfCtybhrr79N4TasiJmPZORCTjpevQ4AvNbAMwF3jAzB4BcPelwB3AMuBh4DPu3h3udVwJPAIsB+4I5wX4GvAlM2skGEP55eCuzeFVleRz5vRq7nlpI1268KOIZDAL/rjPPvX19d7Q0JDy5Ty6tJnLb1nIDRfPYf6cwx4bICIyLJjZQnevP7R9qHVzZZxzjhvN9NGl/OgvjcRi2RncIpL5FCYpFokYnzlrCiu3tPHw0uZ0lyMikhIKk0HwntljmVxdzP/8pZFs7VYUkcymMBkE0YjxmTOnsLxpN48t35LuckREkk5hMkjmz6lh/Mgifvjnldo7EZGMozAZJDnRCJ95xzG8snEXT63cmu5yRESSql9hYmafN7MyC/zSzF40s3elurhMc8GJtYwqzefnT61OdykiIknV3z2TT7n7buBdwAjgE8D1KasqQ+XnRLnsbRN5auVWlm3ane5yRESSpr9hYuHPdwO3hGef2xHml8P42CkTKMqL8gvtnYhIBulvmCw0s0cJwuQRMysFdH2QASgvyuXik8dz38ub2LRzX7rLERFJiv6GyaeBq4CT3X0vkAt8MmVVZbhPnT4RB379tzXpLkVEJCn6GyZzgRXuvtPMPg5cC+xKXVmZrW5EEe+ZPZbfLVjH7v2d6S5HRCRh/Q2TnwJ7zewE4MvAKuA3KasqC/zD2yfT1t7FXQt1eXoRGf76GyZd4b1D5gM/cvcfA7rbUwJm15UzZ1wFtzy3Vicxisiw198waTWzqwkOCX7AzCIE4yaSgE+cOoHVLXt4dtW2dJciIpKQ/obJh4F2gvNNmgluj/vdlFWVJd5z/FgqinK55bm16S5FRCQh/QqTMEBuBcrN7Hxgv7trzCRBBblRPlw/jkeXbaZ51/50lyMiMmD9vZzKhwjuxX4R8CFggZl98Mjvkv742FsnEHPntufXpbsUEZEB62831zUE55hc6u6XAKcA/5q6srLH+MoizpxWzW3Pr6NT94kXkWGqv2EScff4G3FsO4r3ypv4xNwJbGlt5xHdiVFEhqn+BsLDZvaImV1mZpcBDwAPpq6s7PJ300ZRN6KQW57VQLyIDE/9HYD/CnAjcHz4uNHdv5bKwrJJNGJ8/NQJLHh9OyuaW9NdjojIUet3V5W73+XuXwofd6eyqGz0ofpx5OVEuOW5NekuRUTkqB0xTMys1cx29/FoNTPdkCOJRhbn8d7ja7j7xY206npdIjLMHDFM3L3U3cv6eJS6e9lgFZktLpk7gT0d3dz90sZ0lyIiclR0RNYQcsK4Ck6oK+c3z+p6XSIyvChMhphPzJ1I45Y2Xa9LRIYVhckQc/7xY6kszuPnuq2viAwjCpMhpiA3ymWnTeTxFS06TFhEhg2FyRD0ibkTKMyNcuOT2jsRkeFBYTIEVRTl8eGTx3Hvoo007dqX7nJERN6UwmSI+vTpk3DgV8+sSXcpIiJvSmEyRI0bWcR7Zo/ldwvWsWufTmIUkaFNYTKEXX7GZNrau7jl2TXpLkVE5IjSEiZmdpGZLTWzmJnVx7VPNLN9ZrYofPwsbtpbzOwVM2s0sx+amYXtI83sT2a2Mvw5Ih3rlAqzass5+9hR3Pjkau2diMiQlq49kyXA+4En+5i2yt3nhI8r4tp/CvwDMDV8zAvbrwL+7O5TgT+HrzPGl981nd37u/i5juwSkSEsLWHi7svdfUV/5zezsUCZuz/nwXVGfgNcEE6eD9wcPr85rj0jzKgp4/zjx3LTM6+zta093eWIiPRpKI6ZTDKzl8zsr2b29rCtFtgQN8+GsA1gtLs3hc+bgdGH+2Azu9zMGsysoaWlJemFp8oX3zmN/Z3d/OTxVekuRUSkTykLEzN7zMyW9PGYf4S3NQHj3f1E4EvA78ys31cnDvdaDnuFRHe/0d3r3b2+urq63+uSbsdUl/CBk+r47YK1Ou9ERIaklIWJu5/j7rP6eNx7hPe0u/u28PlCYBUwDdgI1MXNWhe2AWwOu8F6usPi71WfMT5/zlTcnf9+9LV0lyIi8gZDqpvLzKrNLBo+n0ww0L467MbabWanhkdxXQL0hNJ9wKXh80vj2jNK3YgiPvm2Sdz54gYWb9iZ7nJERA6SrkODLzSzDcBc4AEzeyScdAaw2MwWAXcCV7j79nDaPwO/ABoJ9lgeCtuvB95pZiuBc8LXGemzZ02hsjiPb/1xme53IiJDimXrl1J9fb03NDSku4yjdscL6/nqXYu54eI5zJ9T++ZvEBFJIjNb6O71h7YPqW4ueXMffEsds2vLuf6hV9nb0ZXuckREAIXJsBOJGN947wyadu3np0/oUGERGRoUJsNQ/cSRXHhiLf/719WsamlLdzkiIgqT4err7z6OgtwI/3rPEg3Gi0jaKUyGqerSfL523rH8bdU27lm08c3fICKSQgqTYewjJ4/nxPEV/Mf9y9m1V1cVFpH0UZgMY5GIcd0Fs9m5r5PrH16e7nJEJIspTIa5GTVlfOptE7nt+fUsWL0t3eWISJZSmGSAL75zGnUjCrn67lfY39md7nJEJAspTDJAUV4O/3nhbFa37OEnjzemuxwRyUIKkwxxxrRqLjyxlp/+dRWvbW5NdzkikmUUJhnk2vccR0l+DlfdtZjumM49EZHBozDJIJUl+fzbe2fw4rqd3PLsmnSXIyJZRGGSYS6YU8sZ06r5ziMr2LhTd2UUkcGhMMkwZsZ/XjgLgGvufkWXWhGRQaEwyUB1I4r4yrnTeWJFC/cu2pTuckQkCyhMMtQlcydy4vgKvvXHpWxta093OSKS4RQmGSoaMb7zgePZ097NN+5bmu5yRCTDKUwy2NTRpXz+nKk8sLiJh5c0pbscEclgCpMMd/kZk5lZU8a19yxl596OdJcjIhlKYZLhcqMRvvPB49m5t4N//+OydJcjIhlKYZIFZtaU88/vmMIfXtrIQ6+ou0tEkk9hkiU+e9YUTqgr5+q7X6F51/50lyMiGUZhkiVyoxG+/+E5tHfG+MqdLxPTtbtEJIkUJllkcnUJ155/HE+t3Mqv/7Ym3eWISAZRmGSZj54ynrOPHcX1D73Ki+t2pLscEckQCpMsY2Z876ITGF2ezxW3LGTzbo2fiEjiFCZZaERxHj+/pJ629i7+8ZaFutWviCRMYZKljh1Txn9ddAKL1u/k2nuW6OrCIpIQhUkWO2/2WD5/9lTuXLiBb/1xmQJFRAYsJ90FSHp94Zyp7Gnv4hdPv05nd4z/N38WkYiluywRGWYUJlnOzLjmPceRE43ws7+uoqvbue7CWeREtdMqIv2nMBHMjK/Nm05u1PifvzTS2NLGDRfPoW5EUbpLE5FhIi1/fprZd83sVTNbbGZ3m1lF3LSrzazRzFaY2blx7fPCtkYzuyqufZKZLQjbf29meYO8OhnBzPjyu6bzw4+cyIrmVt59w1O6bL2I9Fu6+jL+BMxy9+OB14CrAcxsBnAxMBOYB/zEzKJmFgV+DJwHzAA+Es4L8G3g++4+BdgBfHpQ1yTDvO+EGh743OlMrCrmit++yKd//QIrN7emuywRGeLSEibu/qi7d4UvnwPqwufzgdvdvd3dXwcagVPCR6O7r3b3DuB2YL6ZGXAWcGf4/puBCwZpNTLWhMpi7rziNL4271ief3075/7gSa7+wyts3Lkv3aWJyBA1FEZZPwU8FD6vBdbHTdsQth2uvRLYGRdMPe19MrPLzazBzBpaWlqSVH5mysuJ8E9nHsNfv/oOLpk7kf9rWM8Z33mcL/5+Ecubdqe7PBEZYlI2AG9mjwFj+ph0jbvfG85zDdAF3JqqOuK5+43AjQD19fU6qaIfRhbn8c33zeQfzpjMTU+/zm3Pr+PulzZy+pQqPvm2ibxj+igdSiwiqQsTdz/nSNPN7DLgfOBsP3C23EZgXNxsdWEbh2nfBlSYWU64dxI/vyRRbUUh/3r+DD531lR+u2Attzy7lk/f3MCkqmIunTuBD9aPoyRfBweKZCtLx1nPZjYP+G/g79y9Ja59JvA7gjGSGuDPwFTACAbqzyYIixeAj7r7UjP7P+Aud7/dzH4GLHb3n7xZDfX19d7Q0JDkNcsend0xHlrSzE1Pv86i9Tspzc/hovpxXHraBCZUFqe7PBFJETNb6O71b2hPU5g0AvkEexYAz7n7FeG0awjGUbqAL7j7Q2H7u4EfAFHgJne/LmyfTDAgPxJ4Cfi4u7e/WQ0Kk+R5ad0OfvXMGh58pYlud84+dhSXnTaJt02pJDhGQkQyxZAKk6FAYZJ8m3fv59bn1nLrgnVs29PB5OpiPv7WCXzgLXWUF+amuzwRSQKFySEUJqnT3tXNA4ubuOW5tby0bicFuRHOP76GD9WP4+SJI7S3IjKMKUwOoTAZHEs27uLWBWv548tNtLV3MamqmA+cVMv7TqhlfKUu1yIy3ChMDqEwGVx7O7p46JVm7mhYz4LXtwNwwrgK3nv8WN45Y7QG7UWGCYXJIRQm6bNx5z7uf3kT9y7axLLwBMjpo0s5+7hRvH1qNW+ZMIK8nKFwPq2IHEphcgiFydCwbtteHl3WzJ+WbaZh7Q66Y05RXpRTJo3klEkjeeukkcyurVC4iAwRCpNDKEyGnt37O3lu1TaebtzK31Zto3FLGwD5ORFm15Zz0oQRnDS+guPrKhhbXqCBfJE0UJgcQmEy9G1ra+eFNTtoWLOdF9ftYMnG3XR0xwCoKsnn+LpyZtWUMaOmnJk1ZdSNKFTAiKTY4cJE17+QIauyJJ95s8Ywb1Zwibf2rm6WN7WyeMNOFq3fySsbdvHEii3Ewr+HSvNzmDamlOljSpk2qoRpo0uZMrqE6pJ8hYxIiilMZNjIz4kyZ1wFc8ZVcMncoG1fRzevNu9mWdNuVjS38mpTK/e/vInd+7t631demMuUUSVMqS4JfoaP2opCXaRSJEkUJjKsFeZFOXH8CE4cP6K3zd1paW1n5ZY2XtvcSuOWNlZuaeNPyzfz+4YDdzIoyI0wqSoIlmOqizmmuoRjqkuYXF1MQW40HasjMmwpTCTjmBmjygoYVVbA26ZUHTRt+54OVrW00biljVVb2ljV0sai9Tu4f/EmeoYPzaBuRCFTqkuYOrqUqaNKmD6mlKmjSinMU8iI9EVhIlllZHEeI4tHcvLEkQe17+/sZs22PTRuaTvo8cyqbXR0BYP+ZjCpsphjx5Zy3JgyZtSUMbOmnNFlGpMRUZiIAAW5UY4dU8axY8oOau/qjrFu+15e29zKiuY2ljftZtmm3Tz4SnPvPCOL85hZU8as2uCoslk15YwfWaTxGMkqChORI8iJRphcXcLk6hLmzTrQ3tbexatNu1m6aTdLN+1i6abd/OKp1XR2B31lpfk5HFdTxsxw72VmTRlTRpWQG9XJl5KZFCYiA1CSn0P9xJHUx3WXtXd181pzW2+4LNm0i9ueX8f+zqCbLC8aYcqoEo4bW8ZxY4NDmKePKdWhy5IRFCYiSZKfE2V2XTmz68p727pjzutb21i6KTh8eXlTK0+ubOGuFzf0zjOyOI/po0t7x2Kmjyll2mgN9svwojARSaFoxJgyqpQpo0qZP6e2t31bW3twXkxza/Bzcyu3P7+efZ3dAEQMJlYVM2Ns0E02I+wyqyrJT9eqiByRwkQkDSpL8jltSj6nxR26HIs567bvDU/CbOXVpt0sWr+T+xc39c4zqjT/oHGYmTXljBupy8hI+ilMRIaISMSYWFXMxKpi5s0a29u+a28ny5oODPQv3bSLv77WctBlZHrGYY4dW8ax4VhMUZ7+e8vg0b82kSGuvCiXucdUMveYyt62/Z3dvNrcyrJNu1nWtItlm3Zz58IN7Ono7p1n3MhCpo0qZero0t5LyBxTXUxpQW46VkMynMJEZBgqyD1wnbIesZizYcc+ljcH1yl7bXMrKze38eTKlt5DliHoKptcXczk6hImVQZ7QpOqiqgbUaTLyMiAKUxEMkQkYoyvLGJ8ZRHnzhzT294ZnnjZc1b/6pY9rN7axgOLm9i1r/Ogzxhdls/4kUXUVhRSU1FI7YhCxpYXMLqsgLHlhYwoytX4jPRJYSKS4XKjkd6LWJ478+BpO/d28PrWPazZtof12/exbvte1m3fS8PaHTQvbqIr5od8llFdkk91WQGjSvOpKsmnuiSPqtJ8RhbnUVmcT2VJHiOK8hhRlEuOTtLMGgoTkSxWUZTHiePzDrrqco/umLOldT/Nu4JH0679bGltZ0vrflpa21m/fS8vrdvBtj0dHO4ee6UFOYwoyqOiKJfywgOPssJcSgtyKC3Ipawgh+K8HIrzcyjOj1KUF6UwL4fC3CgFuRHyc6JEdWmaIU9hIiJ9ikaMseWFjC0vPOJ83TFnx94OtrV1sK2tne17O9ixp4NtezrYubeTnXs72LG3k537OtmwYx+79nWye1/nG/Z6jiQ3auRFI+TlBI/caIS8aPAzJ2rkRCPkRoycqJEbjRCNGDmRCDlhW/AzeB1MM6KR4L3RiBG1uPZw/ohZ7/zRSIRoBCLW8/rA+yLh+yJxnxPto+3Aew98Tvzn9T43IxKhj7ahHagKExFJSDRiVJXkhydUlvbrPe7O/s4Yrfs7aW3vYk97F23tXext72ZvZzf7O7rZ29FFe1eM/Z0x9nd109EV6310dsfo6A5+dnU7nTGnK3y+p6uLrpjT1e10xWLhT6c75nR2x4i5907v7mmPxQ67dzWURCz4fZuFAWNBMEUi1uc0C4MpYkEwWTj/Ly+tZ0JlcVJrU5iIyKAzMwrzohTmRRmV7mJC7kGwdMUOhE8sfN4TQLGwvdsPBFF3XFtf02PudMegOxYLfrrTHYsRC5/H3vBeett6auh2Jxa2xzx47XHzudO7XHfv/Wz3cL165g9/5uck/6g9hYmICEHABV1m6a5keNKhFiIikjCFiYiIJExhIiIiCVOYiIhIwhQmIiKSMIWJiIgkTGEiIiIJU5iIiEjCzIfDNQRSwMxagLVH8ZYqYGuKyhmqsnGdITvXOxvXGbJzvRNd5wnuXn1oY9aGydEyswZ3r093HYMpG9cZsnO9s3GdITvXO1XrrG4uERFJmMJEREQSpjDpvxvTXUAaZOM6Q3audzauM2TneqdknTVmIiIiCdOeiYiIJExhIiIiCVOYvAkzm2dmK8ys0cyuSnc9qWJm48zscTNbZmZLzezzYftIM/uTma0Mf45Id63JZmZRM3vJzO4PX08yswXhNv+9meWlu8ZkM7MKM7vTzF41s+VmNjfTt7WZfTH8t73EzG4zs4JM3NZmdpOZbTGzJXFtfW5bC/wwXP/FZnbSQJerMDkCM4sCPwbOA2YAHzGzGemtKmW6gC+7+wzgVOAz4bpeBfzZ3acCfw5fZ5rPA8vjXn8b+L67TwF2AJ9OS1WpdQPwsLsfC5xAsP4Zu63NrBb4HFDv7rOAKHAxmbmtfw3MO6TtcNv2PGBq+Lgc+OlAF6owObJTgEZ3X+3uHcDtwPw015QS7t7k7i+Gz1sJvlxqCdb35nC2m4EL0lJgiphZHfAe4BfhawPOAu4MZ8nEdS4HzgB+CeDuHe6+kwzf1gS3KS80sxygCGgiA7e1uz8JbD+k+XDbdj7wGw88B1SY2diBLFdhcmS1wPq41xvCtoxmZhOBE4EFwGh3bwonNQOj01VXivwA+CoQC19XAjvdvSt8nYnbfBLQAvwq7N77hZkVk8Hb2t03At8D1hGEyC5gIZm/rXscbtsm7TtOYSIHMbMS4C7gC+6+O36aB8eRZ8yx5GZ2PrDF3Remu5ZBlgOcBPzU3U8E9nBIl1YGbusRBH+FTwJqgGLe2BWUFVK1bRUmR7YRGBf3ui5sy0hmlksQJLe6+x/C5s09u73hzy3pqi8F3ga8z8zWEHRhnkUwllARdoVAZm7zDcAGd18Qvr6TIFwyeVufA7zu7i3u3gn8gWD7Z/q27nG4bZu07ziFyZG9AEwNj/jIIxiwuy/NNaVEOFbwS2C5u/933KT7gEvD55cC9w52bani7le7e527TyTYtn9x948BjwMfDGfLqHUGcPdmYL2ZTQ+bzgaWkcHbmqB761QzKwr/rfesc0Zv6ziH27b3AZeER3WdCuyK6w47KjoD/k2Y2bsJ+tWjwE3ufl16K0oNMzsdeAp4hQPjB18nGDe5AxhPcMn+D7n7oYN7w56ZnQn8i7ufb2aTCfZURgIvAR939/Y0lpd0ZjaH4KCDPGA18EmCPy4zdlub2beADxMcufgS8PcE4wMZta3N7DbgTIJLzW8GvgHcQx/bNgzWHxF0+e0FPunuDQNarsJEREQSpW4uERFJmMJEREQSpjAREZGEKUxERCRhChMREUmYwkQkQWb2t/DnRDP7aJI/++t9LUtkqNGhwSJJEn+uylG8Jyfu2lB9TW9z95IklCeSUtozEUmQmbWFT68H3m5mi8J7Z0TN7Ltm9kJ4r4h/DOc/08yeMrP7CM7CxszuMbOF4f02Lg/brie4yu0iM7s1flnhGcvfDe/N8YqZfTjus5+Iu1fJreGJaSIplfPms4hIP11F3J5JGAq73P1kM8sHnjGzR8N5TwJmufvr4etPhWckFwIvmNld7n6VmV3p7nP6WNb7gTkE9yKpCt/zZDjtRGAmsAl4huAaVE8ne2VF4mnPRCR13kVw3aNFBJelqSS4CRHA83FBAvA5M3sZeI7gwntTObLTgdvcvdvdNwN/BU6O++wN7h4DFgETk7AuIkekPROR1DHgs+7+yEGNwdjKnkNenwPMdfe9ZvYEUJDAcuOvLdWN/p/LINCeiUjytAKlca8fAf4pvLQ/ZjYtvAnVocqBHWGQHEtw2+QenT3vP8RTwIfDcZlqgjsnPp+UtRAZAP3FIpI8i4HusLvq1wT3RpkIvBgOgrfQ921hHwauMLPlwAqCrq4eNwKLzezF8PL4Pe4G5gIvE9zo6Kvu3hyGkcig06HBIiKSMHVziYhIwhQmIiKSMIWJiIgkTGEiIiIJU5iIiEjCFCYiIpIwhYmIiCTs/wN9ywWBieCiEwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Record the optimization process\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "\n",
    "  \n",
    "# Determine the parameter dimension of the network\n",
    "net = NET(M, weight)\n",
    "\n",
    "# We use Adam optimizer for better performance\n",
    "# One can change it to SGD or RMSprop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# Optimization cycle\n",
    "for itr in range(ITR):\n",
    "\n",
    "    # Forward propagation to calculate loss function\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # Back propagation minimizes the loss function\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # Record optimization intermediate results\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",
    "# Draw a learning curve\n",
    "loss_plot(loss_list)\n",
    "\n",
    "# Record the last two learned unitary matrices\n",
    "U_learned = U.numpy()\n",
    "V_dagger_learned = V_dagger.numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now explore the accuracy of the quantum version of singular value decomposition. In the above section, we mentioned that the original matrix can be expressed with less information obtained by decomposition. Specifically, it uses the first $T$ singular values and the first $T$ left and right singular vectors to reconstruct a matrix:\n",
    "\n",
    "$$\n",
    "M_{re}^{(T)} = UDV^{\\dagger}, \\tag{4}\n",
    "$$\n",
    "\n",
    "For matrix $M$ with rank $r$, the error will decreasing dramatically as more and more singular values are used to reconstruct it. The classic singular value algorithm can guarantee:\n",
    "\n",
    "$$\n",
    "\\lim_{T\\rightarrow r} ||M-M_{re}^{(T)}||^2_2 = 0, \\tag{5}\n",
    "$$\n",
    "\n",
    "The distance measurement between the matrices is calculated by the Frobenius-norm,\n",
    "\n",
    "$$\n",
    "||M||_2 = \\sqrt{\\sum_{i,j} |M_{ij}|^2}. \\tag{6}\n",
    "$$\n",
    "\n",
    "The current quantum version of singular value decomposition still needs a lot of efforts to be optimized. In theory, it can only guarantee the reduction of accumulated errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:46:13.453107Z",
     "start_time": "2021-03-09T03:46:12.949847Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABDbklEQVR4nO3dd3hUddbA8e/JJCSEklACRBASEAgQIEAEVJCiCBZQkSI2sK5lddXFVdd9la2iKLqWtTdUEBEBURSRKiIgvQZpAeklBBIIpJ33j3sTQ0iZQJJJOZ/nmScz987ce2aSzLm/LqqKMcaYysvP1wEYY4zxLUsExhhTyVkiMMaYSs4SgTHGVHKWCIwxppLz93UAZ6Nu3boaERHh6zCMMaZcWb58+SFVDcu9vVwmgoiICJYtW+brMIwxplwRkR15bbeqIWOMqeQsERhjTCVnicAYYyq5ctlGYExZcODAAUaOHElcXByZmZm+DscYAPz8/IiKiuKFF16gXr16Xr3GEoExZ2nkyJH06tWL9957j4CAAF+HYwwAaWlpfPzxx4wcOZJx48Z59ZpKUzU0deVuBjz7BUue7kL/ZyczdeVuX4dkyrm4uDhuueUWSwKmTAkICODWW28lLi7O69dUihLB1JW7efLLtfxVJ3ChZxODj4/nyS+rAXBdh4Y+js6UV5mZmZYETJkUEBBQpOrKSlEiGDNzE9XTDjHYMx8/UQZ7FlA97TBjZm7ydWjGGONzlSIR7ElM4SH/KXjIAMCPDB70/5I9iSk+jsyYc+PxeIiJiSE6Opr+/fuTmJjos1jmzZvHokWLiu14U6dOZcOGDdmPn376aX744YdiO35uX331FaNHj/bquSdOnKBOnTocO3bstO3XXXcdEydOBJz427VrR1RUFNHR0XzxxRfZz1u8eDFdunQhJiaGVq1aMWrUKOLj42nUqNEZV/IxMTEsWbKEUaNG0bBhQ2JiYmjevDkDBw487fM5F5UiEbQNSWGwZz4B4nzAgZLBYM8CokNO+jgyU5lMXbmbS0bPIfKJb7hk9JxiaaeqWrUqq1atYt26ddSuXZvXX3+9GCI9OwUlgvT09CIfL3ci+Mc//sHll19+1vEVZsCAATzxxBNePTc4OJi+ffsyZcqU7G1Hjx5l4cKF9O/fn9WrVzNy5EimTZtGXFwc06dP5/HHH2f58uUADB8+nLfffjv7dzdkyBAiIiJo3LgxP/74Y/Yx4+LiSEpKokuXLgA88sgjrFq1is2bNzN06FB69+7NwYMHz/m9V4pE8HL4LITTV2LzI5P/hn/vo4hMZZPVTrU7MQUFdiem8OSXa4u108JFF13E7t3O8bZu3Uq/fv3o1KkT3bt3z2443L9/P9dffz3t27enffv22V/cY8eOJTo6mujoaF5++WUA4uPjadWqFXfffTdt2rThiiuuICXFKUW/8sortG7dmnbt2nHjjTcSHx/Pm2++yUsvvURMTAw//vgjI0aM4N5776VLly785S9/YdSoUbzwwgvZ8UZHRxMfHw/AuHHjaNeuHe3bt+fWW29l0aJFfPXVVzz22GPExMSwdetWRowYkX1VPXv2bDp06EDbtm254447OHXqFOBMP/PMM8/QsWNH2rZtm2eDadeuXVm/fn324549e7Js2TI+/PBD/vjHPwIwffp0unTpQocOHbj88svZv3//GccZNmwYn332WfbjKVOm0LdvX4KDg3nhhRf461//SmRkJACRkZH89a9/5cUXXwScrsfh4eGAU6pr3bp1nsf87LPPuPHGG/P8fQ8dOpQrrriC8ePH57m/KCpFY3HTk+tBTr8iCZR0Z7sxxWToWz/nu2/lzkRSM04v8qekZTBq+nqu69CQhOOp3PfJ8tP2T/zDRV6fOyMjg9mzZ3PnnXcCcM899/Dmm2/SvHlzlixZwv3338+cOXN46KGH6NGjB1OmTCEjI4Pk5GSWL1/OBx98wJIlS1BVunTpQo8ePahVqxabN29mwoQJvPPOOwwZMoTJkydzyy23MHr0aLZv305gYCCJiYmEhoZy7733Ur16dUaOHAnAe++9x65du1i0aBEej4dRo0blGfv69ev517/+xaJFi6hbty4JCQnUrl2bAQMGcM011zBo0KDTnn/y5ElGjBjB7NmzadGiBbfddhtvvPEGDz/8MAB169ZlxYoV/O9//+OFF17g3XffPe31Q4cO5fPPP+fvf/87e/fuZe/evcTGxrJu3brs53Tr1o3FixcjIrz77rs8//zz2V/iWfr27ctdd93F4cOHqVOnDp999ll2Ilm/fn3255AlNjaWV199FXCu7Fu2bEnPnj3p168fw4cPJygoiCFDhhATE8Orr76Kv78/EydOZNKkSfn+3jt27Fik3kH5qRQlAu5dCKOOwqijHH1wE6c0gMV1BzrbjSkFuZNAlsQTaed03JSUFGJiYmjQoAH79++nT58+JCcns2jRIgYPHkxMTAx/+MMf2Lt3LwBz5szhvvvuA5wr0ZCQEBYuXMj1119PtWrVqF69OgMHDsyunoiMjCQmJgaATp06ZV/Bt2vXjptvvplPPvkEf//8rycHDx6Mx+Mp8D3MmTOHwYMHU7duXQBq165d4PM3bdpEZGQkLVq0AJxqlgULFmTvHzhw4Bnx5jRkyJDsksXnn39+RqIB2LVrF3379qVt27aMGTPmtBJElipVqjBgwAC++OILDh06xMqVK+nbt2+BsWd5+umnWbZsWfYVfb9+/QCoX78+0dHRzJ49m1WrVuHv7090dHS+xymuNecrRYkgp5A6DZh03kNM2R1K9Kl0qgdWuo/AlJCCruAvGT2H3Xl0TmgYWhWA2tWqFKkEkCWrjeDEiRP07duX119/nREjRhAaGsqqVauKfLzcAgMDs+97PJ7sqqFvvvmGBQsWMH36dP7973+zdu3aPF9frVq17Pv+/v6nNYSePFkybXRZMXs8njzbJho2bEidOnVYs2YNEydO5M033zzjOQ8++CCPPvooAwYMYN68efmWZoYNG8Y///lPVJVrr702uztx69atWb58Oe3bt89+7vLly4mNjc1+3KxZM+677z7uvvtuwsLCsksWWdVD9evXZ9iwYQW+15UrV552zLNVqiUCEQkVkS9EJE5ENorIRSJSW0Rmichm92etko6jab8HWZTajK9W7SnpUxkDwGN9W1I14PQr46oBHh7r27JYjh8cHMwrr7zCiy++SHBwMJGRkdlVCqrK6tWrAbjssst44403AKc66ejRo3Tv3p2pU6dy4sQJjh8/zpQpU+jevXu+58rMzOS3336jV69ePPfccxw9epTk5GRq1KhBUlJSvq+LiIhgxYoVAKxYsYLt27cD0Lt3byZNmsThw4cBSEhIAMj3eC1btiQ+Pp4tW7YA8PHHH9OjR48ifV5Dhw7l+eef5+jRo7Rr1+6M/UePHqVhQ2eM0UcffZTvcXr27MnmzZt5/fXXT/vSHjlyJM8++2x2iSQ+Pp6XX36Zxx57DHASadbV/ObNm/F4PISGhgJOiWbGjBlMnDgx3/YBgMmTJ/P9998Xmiy8UdpVQ/8FvlPVKKA9sBF4Apitqs2B2e7jEtWxcSh9whJJnTum2IpWxhTkug4NeXZgWxqGVkVwSgLPDmxbrAMaO3ToQLt27ZgwYQKffvop7733Hu3bt6dNmzZMmzYNgP/+97/MnTuXtm3b0qlTJzZs2EDHjh0ZMWIEnTt3pkuXLtx111106NAh3/NkZGRwyy230LZtWzp06MBDDz1EaGgo/fv3Z8qUKdmNxbndcMMNJCQk0KZNG1577bXsqp02bdrw1FNP0aNHD9q3b8+jjz4KwI033siYMWPo0KEDW7duzT5OUFAQH3zwAYMHD6Zt27b4+flx7733FumzGjRoEJ999hlDhgzJc/+oUaMYPHgwnTp1yq6yyoufnx+DBg3i8OHDpyWjmJgYnnvuOfr370+LFi1o0aIFb7zxBi1bOon/448/pmXLlsTExHDrrbfy6aefZlehhYaGctFFF1G/fn2aNm162vmyGuObN2/OJ598wpw5cwgLO2OdmSKT0voiFJEQYBXQVHOcVEQ2AT1Vda+IhAPzVLXAy6TY2Fg914Vpln7+HJ03/Iet10+nWftLz+lYpnKKjY21BZKMV5544gmWLFnCzJkzqVKlSqmcM6+/TxFZrqpn1CWVZokgEjgIfCAiK0XkXRGpBtRX1b3uc/YB9fN6sYjcIyLLRGRZcfSbbd33bjL9q9Is/vNzPpYxxhRk9OjRzJ07t9SSQFGVZiLwBzoCb6hqB+A4uaqB3JJCnkUUVX1bVWNVNbY4ikLVQ2rj124wrJsMKYnnfDxjjCmvSjMR7AJ2qeoS9/EXOIlhv1slhPvzQGkFlNbhdkg7wS/Tz+w1YIwxlUWpJQJV3Qf8JiJZ9f+XARuAr4Dh7rbhwLTSiing/I5sCGxP+omjpXVKY4wpc0q7E/2DwKciUgXYBtyOk4w+F5E7gR1A3s34JaTV4/MQv8oxrs4YY/JSqolAVVcBeY1+uKw048hJ/PzQzEx2b99Io2ZtfBWGMcb4jF0KA0s/foqwcT04tN8GmJnyxaahLj5FmYYanKmob775Ztq2bUt0dDTdunUjOTmZXr16MXPmzNOe+/LLL3PfffcRHx9P1apV6dChA61ataJz5858+OGHxfxOis4SARDe+QYCJY1NM9/ydSimokvaBx9cCUlnzmZ5Nmwa6uJTlGmowRmcV79+fdauXcu6deuy167OPYMoOLOIZo0AbtasGStXrmTjxo189tlnvPzyy3zwwQfF+l6KyhIB0LhVLHFV2tB4+0QyMjJ8HY6pyOY/DzsXw/zniv3QNg116U5DvXfv3uxpKMCZ+iIwMJBBgwbxzTffkJqamv057tmzJ89pO5o2bcrYsWN55ZVXCvv1lixVLXe3Tp06aXFbMf1N1Wdq6sp5U4r92KZiOuPv8P2rzrwtedvZd+q46tuXq44KVX2mpvPznctVV3zi7E8+dOZrvVCtWjVVVU1PT9dBgwbpt99+q6qqvXv31l9//VVVVRcvXqy9evVSVdUhQ4boSy+9lP2axMREXbZsmUZHR2tycrImJSVp69atdcWKFbp9+3b1eDy6cuVKVVUdPHiwfvzxx6qqGh4eridPnlRV1SNHjqiq6jPPPKNjxozJjm348OF69dVXa3p6ep7727Rpo9u3b9d169Zp8+bN9eDBg6qqevjw4ezXT5o06bTjTZo0SVNSUrRRo0a6adMmVVW99dZbs99TkyZN9JVXXlFV1ddff13vvPPOMz6zsWPH6tNPP62qqnv27NEWLVqoquoHH3ygDzzwgKqqJiQkaGZmpqqqvvPOO/roo4+ecZyVK1dqWFiYdu3aVZ966qnsz1tV9eqrr9apU6eqquqzzz6rf/7zn1VVdfv27dqmTZvTjnPkyBENCgo64/jnKq/vSWCZ5vGdaiUCV5vLbyWRGpz45RNfh2IqqqM7IWt2FVVI3HnOh7RpqH03DXVMTAzbtm3jscceIyEhgQsvvJCNGzcCpy8wk7NaKC9aBuY7szmYXVWCgpkQ/Sqjl8MPiSnZ0wMb47Xbv8l/36ljcDKR3wfOq/P4ArfOu1qdgl+fD5uGOv+YS2Ma6qzEOXDgQPz8/JgxYwatWrXi2muv5ZFHHmHFihWcOHGCTp065RvvypUradWq1dm92WJiJYIcevfuy0mqMHHpuV+pGXOa+c+D5lqcRjOLra3ApqEu/Wmof/rpJ44cOQJAamoqGzZsoEmTJoCTIHr16sUdd9xRYGkgPj6ekSNH8uCDDxYp/uJmiSCH82sH8+eGcXT/+U7S0s5t5ShjTrNrKWSknr4tI9XZXkxsGmrvFcc01Fu3bqVHjx7Zn0NsbCw33HBD9v5hw4axevXqMxLB1q1bs7uPDhkyhIceeojbb7+9SPEXt1Kbhro4Fcc01PlZM/Mj2v38ELuv+oiGna8rkXOYisGmoTZlWVmdhrpcaNN7GFq9Pg23TPB1KMYYUyosEeTiCaiCdLgV/XUmJw5s93U4xhhT4iwR5OFU+1tRYMXU//o6FFOG+fn5WVuSKZPS0tLwK8JkmpYI8hBYN4LlTe6iRvNuvg7FlGFRUVF8/PHHlgxMmZKWlsbHH39MVFSU16+xxmJjztKBAwcYOXIkcXFxp/WPN8aX/Pz8iIqK4oUXXqBevXqn7cuvsdgGlBXgwO7tbFvyNV0H+raPrymb6tWrx7hx43wdhjHnzKqGCrBz9lt0XfM3tm5a7etQjDGmxFgiKMAFfe8nXf3YO9vWNDbGVFyWCAoQWr8xG2p2o83+6Rw/ftzX4RhjTImwRFCIoIvvppYksfp7qws2xlRMlggK0bzL1ez0a8SOTavKxHSxxhhT3CwRFEL8PPx42RSeTBzAml1HfR2OMcYUO0sEXhjQKZLgKh4mL1rn61CMMabYWSLwQo2gAMaeN4c/bxjC0cREX4djjDHFyhKBl1rEXk6IHCdz/WRfh2KMMcXKEoGXmnbqA2FR1FpvaxobYyqWUk0EIhIvImtFZJWILHO31RaRWSKy2f1ZqzRj8poIxN4Be1bw27pFvo7GGGOKjS9KBL1UNSbHxEdPALNVtTkw231cJp1qM5gUqrDzh//5OhRjjCk2ZWHSuWuBnu79j4B5wOO+CqYggdVrs6H3m7RudbGvQzHGmGJT2iUCBb4XkeUico+7rb6q7nXv7wPq5/VCEblHRJaJyLKDBw+WRqx5an3pDdQKC/fZ+Y0xpriVdiLopqodgSuBB0Tk0pw71Rm6m+fwXVV9W1VjVTU2LCysFELN37p5X/Dzi0PIzLA56I0x5V+pJgJV3e3+PABMAToD+0UkHMD9eaA0YzobGYm7uChpJqsWz/J1KMYYc85KLRGISDURqZF1H7gCWAd8BQx3nzYcmFZaMZ2tVlfcQTJVOfnzO74OxRhjzllpNhbXB6aISNZ5x6vqdyLyC/C5iNwJ7ACGlGJMZ6VKcE3W1b+KTvu+Yt/+PTSof56vQzLGmLNWaiUCVd2mqu3dWxtV/be7/bCqXqaqzVX1clVNKK2YzsV5l91PoKQR9+3bvg7FGGPOSZESgYjEishQt2onq7qnLHRBLXUNWsSyoHo/Zu7yI80ajY0x5ZhXiUBE6ovIYmApMJ7fu3iOBV4sodjKvNSrXmFCcidmbyzz7dvGGJMvb0sELwH7gTrAiRzbJ+E0+lZKvaLqcUHNTNbMs4nojDHll7fVOpcBl6nqEbexN8tWoHGxR1VOePyE58Jm0m7XeH777WrOPz/C1yEZY0yReVsiqAqk5rE9DDhZfOGUPxF97iVAMmiw9Qtfh2KMMWfF20SwABiR47GKiAdnTqDZxR1UeVInIhoiuhOw6iPItEZjY0z5420i+Atwt4jMAgJxGog3AJcAT5ZQbOVGesfbIXEnq+ZbW4ExpvzxKhGo6gagLbAI+B4Iwmko7qCqW0suvPLBr9U1JBBC0sZ5vg7FGGOKzOsxAKq6D3imBGMpt/wCAvF/cCnd6zTwdSjGGFNk3o4j+KOI3JLH9ltE5P7iD6v8qekmgVOnKnXbuTGmHPK2jeBh4Lc8tscDjxRXMOXd8gn/4OCz7Thx0pKBMab88DYRNMKZEC63Xe4+A4Q0bEEj9rNi1gRfh2KMMV7zNhHsA2Ly2N4ROFRs0ZRzzS65gYNSh6prxvk6FGOM8Zq3iWA88IqI9BGRAPd2BfAy8GmJRVfOiCeA3U2H0CltBXEb1vg6HGOM8Yq3ieAZ4CdgJs5cQyeAb3G6k/5fyYRWPl3Q737S1Y+9c970dSjGGOMVr7qPqmoaMExEngY64KwrvEpVN5dkcOVR9bDGTDr/cd7d0YCOKWmEVA3wdUjGGFOgIq1HoKqbVfVzVZ1kSSB/ra68j01pYUxZscvXoRhjTKG8HlAmIkNxZiGtR64EoqoDijmuci26YQjD6u/Cf/5U9OK3yTVjqzHGlCneDigbA3wCRACJwOFcN5PLsPC93JL6OTs2rfR1KMYYUyBvSwS3AcNU1eZa9lKLfveiv75GRPwkiOro63CMMSZf3rYR+AGrSjCOCicotAHSegCs+hRNPVH4C4wxxke8TQRvA2fMNWQKlhozHE4eZd6Xb/s6FGOMyZe3VUOhwE0i0gdYA6Tl3KmqDxVzXBVClWaXEhfciSoBReqcZYwxpcrbRNCa36uGonLt02KLpqIRIeovc3wdhTHGFMjbAWW9iuuE7hKXy4DdqnqNiEQCnwF1gOXAraqa1/rI5VZa6ik2r/uF1h27+ToUY4w5gy/qLP4EbMzx+DngJVW9ADgC3OmDmErU+g/+SONpN7Dv4EFfh2KMMWfwOhGISC8ReVtEvhOROTlvRThGI+Bq4F33sQC9gaxuqR8B13kdfTnRoNttVJeTrPvuPV+HYowxZ/B2QNkInEnmagA9gYNALZxpqDcU4XwvA38BMt3HdYBEVU13H+8CGuYTwz0iskxElh0sZ1fWDVp3Y0dAUxpvm0B6eoavwzHGmNN4WyIYCfxRVYfh9Bh6UlU74Iw2TvbmACJyDXBAVZefTaCq+raqxqpqbFhY2NkcwndEON72NlpoPMsW/eDraIwx5jTeJoKmQNY32Cmgunv/NWCEl8e4BBggIvE4jcO9gf8CoSKS1WjdCNjt5fHKlRaX38Fxgkj4ZZKvQzHGmNN4mwgO41QLgfNFHe3erwNU9eYAqvqkqjZS1QjgRmCOqt4MzAUGuU8bDkzzMqZyxT84hC86fsIDB69l52EbaWyMKTu8TQQ/Ale49z/HWa3sA2ACMOscY3gceFREtuAklgrbotq3R3f8/DyMX5LX8s/GGOMb3g4o+yMQ5N5/FkjHqer5HPhXUU+qqvOAee79bUDnoh6jPGoQEsS/zltEq6V/51SfnwgM8HoWcGOMKTHeDihLyHE/E6fvvzkLHS9oRMtDcezdMI/w9pf7OhxjjPG6+2iGiNTLY3sdEbH+kEXQovetaFAI4Zsn+DoUY4wBvG8jyG+JrUCgQk0HUdKkSjWk/TB0wzQS9ttSlsYY3yuwakhEHnXvKnCviOQcM+ABugNxJRRbhXWq/W0ELnmT5dNep889z/o6HGNMJVdYG8GD7k8B7gJyVgOlAvHAvcUfVsUWeF4bVrV4kMat+vo6FGOMKTgRqGokgIjMBQaq6pFSiaoSiLmpyJ2tjDGmRHjVRqCqvXInARG5QESC8nuNKVz8uiXM/diqhowxvuVtr6H/iMhw976IyA/Ar8BeEelSkgFWZEnLxtN9y/Ns2GTNLMYY3/G219DNwCb3/pVAe6ArMA4YXQJxVQoRfR/AXzLZPfstX4dijKnEvE0E9XGmiAa4CvhcVZcCrwIdSiKwyqBGeAt+rX4h0funcexEiq/DMcZUUkWZdK6Je/8KYLZ735/8xxgYLwR2vZtwOcyy7yf6OhRjTCXlbSKYDIwXkVlAbWCmuz0G2FICcVUaTS4ayG+eRqyNi0NVfR2OMaYS8jYRPAq8grMaWR9VPe5uDwfeKInAKg1PAIv6fsNLid1ZtsN65xpjSp+33UfTVfVFVf2Tqq7Msf0lVX235MKrHPrHNKJGkIevf/zF16EYYyqhfAeUiUhHYJWqZrr386WqK4o9skokuIo/74ZNoumWWRw+up46ITUKf5ExxhSTgkYWLwMaAAfc+0reDcOKM++QOQeNLhxA2IwvOLblO+g02NfhGGMqkYISQSRwMMd9U4Iaxl4DPzWm5rpxlgiMMaUq30Sgqjvyum9KiJ+HzI7D8Zv7T35dv4IWbQqsjTPGmGLj7RQTTUXkURF5TUReFZFHRMRKCcUsrf1NpOHhtznW/m6MKT2FLlUpIn/GWafYg9NeIEAY8JyIPK6qL5VsiJVHYOh5xPefSLc2l/g6FGNMJVJgiUBEugHPA2OAMFUNV9UGQD3gRWCMiNi3VjGK6NSHwKBgX4dhjKlECqsaug8Yp6pP5VrA/rCqPgl8AtxfkgFWRmumv8bC0deSnpHp61CMMZVAYYmgK/BhAfs/dJ9jilGV9CS6nZzH/f96mSVPd6H/s5OZunK3r8MyxlRQhSWCBsC2AvZvxZlmwhSjLQ2u4ZQG8GjGe1womxh8fDxPfrnWkoExpkQUlgiqAqcK2J8KBHpzIhEJEpGlIrJaRNaLyN/d7ZEiskREtojIRBGp4l3oFdez8w/yQ2YHWsou/EQZ7FlA9bTDjJm5qfAXG2NMERXaawi4WkSO5rMvtAjnOgX0VtVkEQkAForItzgT2r2kqp+JyJvAnVTyiez2JKYg/r/PROpHJg/6f8kziXf4MCpjTEXlTSJ4r5D9Xs2drM4cy8nuwwD3pkBv4CZ3+0fAKCp5ImgbkkLvk6sQd0KPQElniGc+B6U2P/8aw0UtGvo2QGNMhVJg1ZCq+nlx83qeIRHxiMgqnPEIs3DaGBJVNd19yi4gz285EblHRJaJyLKDBw/m9ZQK4+XwWUiu/OpPBn/2fE6bST1g2fuQkeaj6IwxFY236xEUC1XNUNUYoBHQGYgqwmvfVtVYVY0NCwsrqRDLhKYn1xMo6adt85dMMmtFUq1eBHz9CCdf6sjbr/ybA0dtiUtjzLnxpmqo2KlqoojMBS4CQkXE3y0VNAKsa8y9C/Pc7AegCptncfKbp+l87DtqBj8BwPFT6VQL9Mmv0xhTzpVaiUBEwkQk1L1fFegDbATmAoPcpw0HppVWTOWSCLS4gtA/LaL9I1MJCvBw6vAOto++iLfeeZ1tB5J8HaExppwpzaqhcGCuiKwBfgFmqerXwOPAoyKyBahD4Y3TBsDPD6lWBwA9toeGVY7zh91/JfG1nrz1wfv8lnDCxwEaY8oLKY8LpsfGxuqyZct8HUbZkpFG0pKPyJz7HCFpB/gpM5qZHV7j/t6taBAS5OvojDFlgIgsV9XY3NutUrmi8ARQ4+K74MJbOPbTO2Ru2MD4X/Yycfk+Huzgz419L6Vuda/G/hljKhmvEoGI1MLp398LZ+bR06qUVLVesUdmzk5AEDV7Pkj3njA34QRffP0196+5h8WbLqbuXWOhntcdtYwxlYS3JYJxQBucAV/78XIQmfGt82sH88iQviT88DBdVr0N/+tKaptBTKx2M9f17kaNoABfh2iMKQO8aiMQkSSgh6quKPmQCmdtBGfh+GH46WXSl7xNUro/u29fQXREfV9HZYwpRfm1EXjba2hrEZ5ryqJqdeCKf+L/8GpS+7/hJAFV5rz3VybMWcap9AxfR2iM8RFvv9z/BDwrIu1FxOspJUwZVKMB9WMHAJC+ezU9fnuDa+dfxcRn72LywrWk2WI4xlQ63iaCLThTUq8AUkUkI+et5MIzJcm/UQx+f1xKcsQV3JIxhT6z+jBu9P1MW7qZjExrBjKmsvC2jWABUAt4kzwai1V1colElw9rIyh+um8dh6Y/Q+beNXQ/MYbG9WrxyGXNubJtOH5+4uvwjDHF4FzHEcQCnVV1XfGGZcoKaRBN2N2TyTyRyMtbU3jt+3U0nHw173x/BXc8NIqAwKq+DtEYU0K8rRraANQsyUBM2eAXHMpVbcOZfkcUjcLq8IfjbxLwvwthxTjW/XaY8jgS3RhTMG8Twd+AsSJyuYjUF5HaOW8lGaDxDU+t86n7x1lw6xSoXg++epDgdy7mm0WrfB2aMaaYeVs1NMP9+T2ntw+I+9h6ElVEItCsNzTtRdrGb8j8aRK9OkUDsHLFUvzqNqd941pMXbmbMTM3sScxhfNCq/JY35Zc18FWUTOmvPA2EfQq0ShM2SZCQOtruKD1Nc7j44eImt6fzRnhPF1zBBMTW1AzPYHPqrzKHxMf4skvUwEsGRhTThTaayhroXngNlXdVCpRFcJ6DflYRjonV4zn1A/PEnJqD0syo0jMrEYfzwo+ybiMp9PvoGFoVX56orevIzXG5HDWI4tVNQ2IxOYXMlk8/gRdeBshj63m/9Jup6nsoa//cvxEGexZQBiJ7Em0JTSNKS+8bSz+CLi7JAMx5ZB/FebUGMCsjE6kqdNM5Ecmn1T5NyP9PyNu41ofB2iM8Ya3bQTVgJtFpA+wHDiec6eqPlTcgZny4W89atPru4UEuAPMAyWdC9hDc/89yMTpcMHlfB98FaurdmZkvzaI2OA0Y8oabxNBK5zpJQCa5tpnVUaV2JWHx5Hhx+l/BX7++LW9AWpFwIqPuCJpFkm17kKufBGA79fvo0tkHUKCbRpsY8oCrxKBqlqvIZO3XUvxaNppmzyaBgfWw8C34NLH4NfvuD68AwBHVkwjc8qr/Ik+1Ii6jBtiG9O9eRgem8bCGJ8p0lKVIhIEXIBz/bdVVU+WSFSm/Lh3YcH7Pf7Q6prsxqhQSeLyalvpd/IXdm5+j0829ObZ4D706tiGQZ0acUG96iUesjHmdN5OOhcA/Af4I1AFZyDZKeBV4Cm3Z1Gpse6j5Vz6Kdg4ncxf3sNv5yL2BDSm+/HRZGRCh8ahDOrUiOs7NCS4ii2pbUxxOtdJ554DhgH34owpAOgOPIvT82hkcQRpKgn/QGg7CL+2g+DARs5L2sfP9S7i62XbuXDhHUz9ujNpLf4PqtRl39GThNUItKojY0qQtyWCfcAdqjoj1/argXdVNbyE4suTlQgqqIRt6Bd3IntWQEAwRN/An7d14mBIG8bd0dnX0RlT7p3rUpUhOMtV5rYVCD2HuIz5Xe2myD1z4Z550HYwum4yLx59mPtaHAPg+Kl0bn53MROW7uTYyVKtjTSmQvO2RLAYWK6qD+Ta/gYQo6oXeXGM84FxQH2cxua3VfW/7uylE4EIIB4YoqpHCjqWlQgqiZNHYdO30G4oiJAw9Ql+XL+D/yVfSrwngn7RDRjUqREXN6trVUfGeCG/EoG3ieBSnBlIdwOL3c1dgfOAK1W1kK4jICLhQLiqrhCRGjgD064DRgAJqjpaRJ4Aaqnq4wUdyxJBJfX1o+jKT5CMU+yo1o43j1/KlydjqRNSk4EdG3FDp0ZE1q3m6yiNKbPOKRG4BzgPeACIcjdtBP6nqnvOMqBpwGvuraeq7nWTxTxVbVnQay0RVGLHD8Pq8bDsfUjYxrZmt/KP9NtY8OtBMhUujKjFx3d2ISjAZkY3JrdzTgTFHEwEsACIBnaqaqi7XYAjWY9zveYe4B6Axo0bd9qxY0dphWvKosxMiF8AIedDnWYkxP3IiZn/ZE71q7ltxP3gCeC9hdtp2zCEzpG2dpIxcJbdR71dfUxVE4oQSHVgMvCwqh7LOfeMqqqI5JmZVPVt4G1wSgTens9UUH5+0LRn9sPaJFE7Yze3/fY0vPQ66e1v4cufL2Bfpxg6R9YmI1PZdeQEK3cm2iI6xuRS2DiCQxQ+l5B6cRwge2DaZOBTVf3S3bxfRMJzVA0d8OZYxpwm6ipo0Rc2z4Jl7+P/01i+rlaXYz2cGVB/3nqYW95bgp9ApvsXvTsxhSe/dPZbMjCVWWFf4AXNMdQP+BOQ7s2J3Gqf94CNqjo2x66vgOHAaPfnNG+OZ8wZ/DzQsp9zS9yJHNxESPWqkJlJ5x9H8GhQQz4+eSmgvFblVf6Y+hAH00J59tuNlghMpVbkNgIR6QCMwRlZ/BbwT1U96MXrugE/AmuBTHfzX4ElwOdAY2AHTvfRAquarLHYFMnxQ/DFHbB9PmnqYbfWobEc5NOM3vxf+p0AtKxfg55RYfy5T0uq+Hs7vMaY8uVcB5QhIpEiMh5YChwGWqvqQ94kAQBVXaiqoqrtVDXGvc1Q1cOqepmqNlfVy4vS3mCMV6rVheFfcWOV1/gsoydN5AB+ogz1zCeMRBoGnaJBNZi98UB2EvhoUTzfrt3r48CNKR2FJgIRqSMi/wXigAbAxao6VFXzGmlsTJl145W98ffzIxWna6kCj1SZyvstFvPRoWF83/B9WDMJTUlk/JKdzNq433meKm/N38qq3xLJzLR+CqbiKazX0FPAYzgjfq9V1e9KIyhjSsJ1F3jI8F+AJ/P31dSGeubj6Xg71EjFb9MM2DgV8QvguxZ9Se7/AQC7jqQw+rs4VKFu9Sr0aFGPXlFhdG8eRkhVW1zHlH+FNRb/E0gBdgH3i8j9eT1JVQcUd2DGFLv5z+PJ1TvZIwpbfoD+L8PVY2HXLxD3NZKRRo2qVQA4f8FI1vdpxk8BXZm+K5gfNu5n8opdePyETk1q0TuqHr1a1qNF/eq2FKcplwpLBOOwpShNRbFrKWSknr4tI9XZDs7YhMZdnFuW1BNwYAPBqz6lD9AnLIrMi69iQ9hVfLu/BnPiDjL62zhGfxvHxHu60qVpHY4cTyUowEPVKja62ZQPPhlZfK6s15ApdYm/waYZEPc1xP8E/f8LHW+F5IMc3r6SWcebccOFkQR4/PjPjI1MWLKTZf93OYH+Ho6fSqdaoC2yY3zvXBemMaZyCz0fuvzBuZ1IAI9TbcSGqdSZMZIbg0JhTz+IuporW3Ti/FotCfR3SgS3vb+UxBOp9GpZj95R9YiNqG1dVE2ZYiUCY85F6gnYOgfivoFfv4WUI86iOn+Og6AQyMxk3JKd/LDxAIu3HiY1I5NqVTx0a16X3lH16NmyHvVrBvn6XZhKokxNOneuLBGYMikjHXb8BHtXwyUPOdvG3wipyRB1DSnN+vLTwWDmbDrA3LgD7D16EoA259Xk/p4XcHW7Mxf6m7pyt82NZIqNVQ0ZU9I8/tC0h3PL0rAjrPsSvnucqjzO5eHtubzzH9DrbmLT/iTmxh1kbtwBMtwLsvhDx3n5h1956LLmrNl1lCe/XEtKmtPd1eZGMiXFKiqNKUk9/gIPLIYHV0Cff4AnEJL3IyJE1fZwX+qHfH6VMKBtfQB2JJxgweZDBHj8GDNzU3YSyJKSlsGYmZt88U5MBWYlAmNKQ51mcMmfnFtWdeze1bD4DVj0ClQLg5ZX0SPqGpY9fil+VYLYk5gCQBhHfp8kj9Ds7cYUFysRGFPasgadNbkY/rINbngPIrrBuskwfjB+CZsBaBeSQnVO8JD/FC6UTTzo78zcrkCfsfN5+Ydf2bw/yUdvwlQk1lhsTFmRfsppbG7aC0TY8f4IwndMw4PiEeWkBnB5xqt0bdeK346ksDQ+AVUY2KEhY4fG+Dp6Uw5YY7ExZZ1/IDTrnf2wyRUPkvj5BkKOOW0CgZLGdyH/ofoQp8H4wLGTfLtuH/VrBgJw9EQaN7+3mMf7RdG9eVjpx2/KLasaMqasCmlI6IkdZM1eJED1lD2QtB8yM6k3/VaGZ06hX/gJAA4mn6RqgIcaQc5EeCt3HuG1OZvZdjDZN/GbcsNKBMaUVfOfB808fZv4wfznnN5Ixw/CD6OcW/1oLmjVn0k33gShoQAs3pbAC9//ygvf/0qr8Jpc0y6cq9qGE1m3Wmm/E1PGWRuBMWXVm91g39oztzdoC/cudO4n7oSNX8PGr2DnYrjlC7jgcmf78UPsrRbFjHX7mbF2L8t3HAGcAWxXtQ3n6rbhRFhSqFRsZLExFV3SPgiuA54A+OHvsHAshJwPrfpDq/7sqdGOGesP8M3avazcmQjAldENeOOWTr6N25QaSwTGVCYnEuDX72DDV85cSBmnoFakM7DNz4/diSl8u3YvwVX8ualLY9IyMhn+/lLu6h5J76j6vo7elBDrNWRMZRJcG2Jucm6nkmDz904js5/TP6Th1MHcFXI+tB4AafXYn5zJ8dQMMt0mie2HjjNrwz6uahtOo1rBPnwjpjRYIjCmogusAdE3/P44PRVqNnRmTF09HqpUp1HzK5h2zT3QxCkNLNxyiP/MiOM/M+Jof34o17QN58q2DSwpVFBWNWRMZZWeCvELYON0Jylc/nfocDMc2wvb5vFbWA++3pzCN2v3sG73MQBizg/lmnbhXNk2nIahVX38BkxRWRuBMSZ/mRnOzb8KLHsfvn4E/Pwh8lJoNYDf6vdi+tZ0vlmzl/V7nKTQq2UYH9zeOfsQNmV22WdtBMaY/Pl5nBtAxxHQoD1snOY0Nn/9MOf7+XP/Y1u4v+cFbN9/hBkbDpOR6VxEqioD/7eI9XuOkprhbLMps8uXUhtZLCLvi8gBEVmXY1ttEZklIpvdn7VKKx5jTD78/KBRJ2fa7IdWwr0/wdVjoarz7xn5/R08sOUeHgr8Gg5vJfFEGmt3O0kgjCNMrPIPwkgkJS2Df3+zkeOn0n38hkxhSq1qSEQuBZKBcaoa7W57HkhQ1dEi8gRQS1UfL+xYVjVkjA8tetVZbGfPCudxvTb8ZfclfJ7Rk3/6v8/Nntl8knEZT6ffkf2SxrWDadmgBlENamT/jKhTDX+PzXJTmnxeNaSqC0QkItfma4Ge7v2PgHlAoYnAGONDFz/o3BJ/g7ivYcNXNKl6krDkIwz2zMdPlKGeebyVfg0p1Rpx+8URxO1LIm7fMWZv3I9bo8QX915EbERtVv2WyJJth7m5axOqB1pttS/4+lOvr6p73fv7ABvJYkx5EXo+dL0Put5HwxW7eHTqw3hwBiIESjrzAx8hoXZX6lW5AgbcBDU6cjItgy0Hktm0L4lW4TUBWLT1EC9+/yvDL44A4L8/bGbR1kNu6aEmLd1ShCWJklNmPllVVRHJt55KRO4B7gFo3LhxqcVljCncdc39yfBfgCfz96U1/fz8qJd5GGb/HVpeCTUaELR7MdEHNhDdtCdUcRqn7+95ATd3aUJQgPO4RpA/qRmZfLF8F8dTfz9eo1pVs6uWWoeHcHW78FJ9jxWZrxPBfhEJV9W9IhIOHMjviar6NvA2OG0EpRWgMcYL85/Hk+s6zs/PA5Hd4fZvnKU4wRmvsPh1536N86BpD4jsQUi7odmvu6NbJHd0iyQzU9mdmELcviQ27Tvm/kxi7qaDNKkdnJ0Inpm2jpDgKjzapwUACcdTqRUcgGStBGcK5etE8BUwHBjt/pzm23CMMWdl11LISD19W0aqs716vd+39f03XHgnbJ8P2+bDrzOdWVNjhjn7V34CQaEQ0Q2/qqGcXzuY82sH06f177XGp9IzOJh0Kvtx8qkMAjxOzyRVpeeYuYgILeu7DdPhTuN0i/o1stdqABv3kFNp9hqagNMwXBfYDzwDTAU+BxoDO4AhqppQ2LGs15AxFURmJiTvg5rngSr8t50zhbb4wXkdILIHRF0Njc7o6JKn9IxMxi/dmV162LQvieQc3VcbhjrVS41qVeXzZbtISfu96qlqgIdnB7at0MnARhYbY8q+9FTYvQy2zXNKDLuXQec/QL//QEYa/PwaRFwK58X8PgCuAKrKriMpTlLYn5RdzbT/2EmOppw5vuG8kCDmjOyZ3V5R0VgiMMaUP6eSIP0UVKsLe1bB2z2c7YEhENENmvZ0ZlCt0aBIh4184hvy++YL8AjtGoVyYURtOkfWolOT2oRUDcjn2eWLz8cRGGNMkQXWcG7glAJGbobtC5wSw/b5sOkbqNfKSQT718PeNU4DdM3zCjzseaFV2Z2Ycsb2ejUCub5jQ37ZnsB7C7fx5nxFBFrWr0HnyNrOLaI29WoGFf979SErERhjyq+E7c6U2v5VYM6/YcHzzvY6zZ2E0LQntOjnrNqWw9SVu3nyy7UFthGkpGaw8rcj/LL9CL/EJ7Bi5xFOuN1ZVz3dh9DgKmzal0RQgB9N6pSPJT+tasgYU7FlZsKB9U7bwrZ5sGMRePzhL9ud9oS4GRBQFRp3hYCqTF25m/e/+5mnUsbwr6p/4c5+XQtsKE7LyGT9nmPE7T3GjZ2dsUx3fvgL2w8fZ86fewIwe+N+GoQEEdWgJh6/std91RKBMaZySU+FI9shrKXz+H8XwYEN4AmE8zs7JYZ9a531GDrdDteMLfIpthxI5lDyKbo2rUNmptLxX7NIPJFGjSB/YpvU4kK3KqltoxAC/X3fAG2JwBhTuZ1Kgh0/u2MY5sH+dU43Vc0E/yDo8QS0uALCWmUv6VlUu46c4Jf4BJZud25bDx4HINDfj5jzQ+kcWZu+bRoQ3TCkGN+Y9ywRGGNMTlPug7WTIDMN/AKcnwDBdZ0R0ZGXQosroebZT2VxKPkUy+ITWOq2M6zfc5TH+kZxX89mHDmeymtzt3BTl8Y0C6teTG+qYNZryBhjsiTtg/Vf/v7ln5kG/oFw2SjYu9opNayfAoNrQZvr4Ui8U5qIvBRCvB9wVrd6IP2iw+kX7SSTpJNp2bOvbjmYzCeLd9C3TQOahcGSbYeZumoPnSNrcWFE7dPWhy7pUdCWCIwxlc/8550qoZxU4fAWGPiWe38r1HCnttj0HXznzpBfu5mTEJr2cHokBXi/dnPOKS4ujKjNmlFX4HHnRNpx+ARfr97DhKU7AWcU9IURtQjw9+OrVXs4le7EWxKrv1nVkDGm8nmzm9NQnFuDtnDvwjO3Z2Y6Dc3b5zvjGOJ/gvQUeHwHBFaHLT84I5+bXAJBNc86rIxMJW7fMX7ZnsBSt0rpUPKpPJ/bMLQqPz3Ru0jHtzYCY4wpLhnpcOhXqN/aefzRACdJiMedI+lSuOAyZ/TzOVBVmj45I89R0AJsH311kY6XXyKwdeKMMaaoPP6/JwGAmz6H4V9D9z+Dnz8segV+fPH3/b+858yymp565rEKICKcF5p31VN+28+GtREYY8y5Cghyexp1B56CU8lw4pCz7+QxmPEYaAYEVIMmFzklhqhroE6zQg/9WN+WeY6Cfqxvy2IL3xKBMcYUt8Dqzg2cNoPHtkD8Qqd9YfsCmPU0BAQ7iSBpP2yY5iSHsJaQa0GdrAbhkuw1ZG0ExhhT2pL2Od1Vq9aCdV/CF7c726vXdxJC5KXQ+loICjn9NV/cDoM+/L03UxHZOAJjjCkrck6bHT0QGnb8vbSwfYEz0K3ZZU4i2L7ASQJb5zjtDPOfO6vpMApiicAYY3ytVoRz63jb72MYsgaurfwE1kz8/bmrPoUej591qSAv1mvIGGPKEhGoe8Hvj697A1pd63RNBWcg3PznivWUlgiMMaYsO34QNs90eh0BZKQ6pYKk/cV2CksExhhTluU5HUbxlgosERhjTFm2a6lTCsgpI9XZXkyssdgYY8qyvOY+KmZWIjDGmErOEoExxlRylgiMMaaSs0RgjDGVnCUCY4yp5MrlpHMichDYcZYvrwscKsZwSlp5itdiLTnlKd7yFCuUr3jPNdYmqhqWe2O5TATnQkSW5TX7XllVnuK1WEtOeYq3PMUK5SvekorVqoaMMaaSs0RgjDGVXGVMBG/7OoAiKk/xWqwlpzzFW55ihfIVb4nEWunaCIwxxpyuMpYIjDHG5GCJwBhjKrlKkwhE5H0ROSAi63wdS2FE5HwRmSsiG0RkvYj8ydcxFUREgkRkqYisduP9u69jKoyIeERkpYh87etYCiMi8SKyVkRWicgyX8dTEBEJFZEvRCRORDaKyEW+jik/ItLS/UyzbsdE5GFfx5UfEXnE/f9aJyITRCSo2I5dWdoIRORSIBkYp6rRvo6nICISDoSr6goRqQEsB65T1Q0+Di1PIiJANVVNFpEAYCHwJ1Vd7OPQ8iUijwKxQE1VvcbX8RREROKBWFUt84OeROQj4EdVfVdEqgDBqpro47AKJSIeYDfQRVXPdrBqiRGRhjj/V61VNUVEPgdmqOqHxXH8SlMiUNUFQIKv4/CGqu5V1RXu/SRgI9DQt1HlTx3J7sMA91ZmrzBEpBFwNfCur2OpSEQkBLgUeA9AVVPLQxJwXQZsLYtJIAd/oKqI+APBwJ7iOnClSQTllYhEAB2AJT4OpUBuVcsq4AAwS1XLcrwvA38BMgt5XlmhwPcislxE7vF1MAWIBA4CH7jVbu+KSDVfB+WlG4EJvg4iP6q6G3gB2AnsBY6q6vfFdXxLBGWYiFQHJgMPq+oxX8dTEFXNUNUYoBHQWUTKZPWbiFwDHFDV5b6OpQi6qWpH4ErgAbeasyzyBzoCb6hqB+A48IRvQyqcW4U1AJjk61jyIyK1gGtxku15QDURuaW4jm+JoIxy69onA5+q6pe+jsdbblXAXKCfj0PJzyXAALfe/TOgt4h84tuQCuZeDaKqB4ApQGffRpSvXcCuHKXBL3ASQ1l3JbBCVff7OpACXA5sV9WDqpoGfAlcXFwHt0RQBrmNr+8BG1V1rK/jKYyIhIlIqHu/KtAHiPNpUPlQ1SdVtZGqRuBUB8xR1WK7sipuIlLN7TCAW81yBVAme76p6j7gNxFp6W66DCiTHRxyGUYZrhZy7QS6ikiw+/1wGU7bYbGoNIlARCYAPwMtRWSXiNzp65gKcAlwK87ValbXtqt8HVQBwoG5IrIG+AWnjaDMd8ssJ+oDC0VkNbAU+EZVv/NxTAV5EPjU/VuIAf7j23AK5ibXPjhX2GWWW8r6AlgBrMX57i626SYqTfdRY4wxeas0JQJjjDF5s0RgjDGVnCUCY4yp5CwRGGNMJWeJwBhjKjlLBJWYiMwTkdd8dG4VkUG+OLc3ynp8xcGdxXKUF8+bKyK3lUJIeZ27wN+DO+vtDaUZU0VkiaCCcgd5/c+dwviUiOwXkdki0ifH0wYCT/oqxuImIh3dL47u+eyfKCKLSjuuguT3RSciH5aFKbJF5GrgfODTHNvi3bhVRFLcKacfcwc6lbZ/AqNFxL7LzoF9eBXXZJypCO4EWgDXAN8CdbKeoKoJ7uym5Y6I+Of+4nFnbF0F3JHH8+sA12EzjhbVn4APVTUj1/Z/4AwkbIUzGdp/AF9MiDcDqIEzTYQ5S5YIKiB3uofuwBOqOltVd6jqL6r6gqp+luN5p1UNuVd6fxORt9xFOnaJyGO5jt1CROaLyEkR2SQiV4lIsoiMcPdHuFeKsbleV1gRf7R7vBQ3judzLrwhIqPcqowRIrIVOAXkNbPlu8Bgd8K+nG5xXzNRRPqJyI8ickREEkRkpoi0KiA2r96TiDQUkc/c4x4RkW9EpHl+xy0KERkoImvczyfB/R3Uz7G/vzizk54Uke0i8m9xJlPL2l9PRKa5r98hImckyzzOGYYzx830PHYnqeo+VY1X1XeBNTjTX2S9tpl7vn0iclxEVogz4V/O4xf695ZHTI+LyCER6QrOZIc4yWBYYe/H5M8SQcWU7N4GSNFXMXoEZwh7R+A54HlxV5lyi99TgHSgKzACeAYILIaYj+NcybcC7seZB+ipXM+JBG4CBgPtgZN5HOdTwAMMzbX9TmCiqh7HSSAv45SYegJHgek5vziLSkSCcSbbOwn0AC7CmS74B3ffWRORBjgT5H2E8/lcCnycY39fnPf9GtAG53McxOnTO3wIXIDzxX4dcBsQUcipu+Ekz3znNhJHTzeutBy7quOUQPvg/K4mA1+KSFSuQ+T795bHeV7AmcKiR65Fj5bifObmbKmq3SrgDbgBZyGekzhzLL2As/pSzufMA17L8TgemJDrOZuBv7n3++IkgYY59l+MM1/+CPdxhPs4NtdxFBiU3+M84r8X2JLj8SicL5r6Xrz3T4BFOR5f6J6vSz7PrwZk4Ez3fEZ83rwnnC/fzbjTtrjbPMBhYEgBseb5OeB8cX/t3u/oPq9JPsdYAPxfrm3X4VwMCE7VoAKX5NjfxH3PowqI7WFgRx7b43ESRDKQ6h47Bbi4kN/L4qy/JW/+3nJ8PkOBD4Bf8/oMcKaQzgT8S/v/rKLcrERQQanqZJx5y/vjXJldDCwWkb8W8tI1uR7vAeq596OAPepOi+z6hWJY4EVEBonIQrcqIRl4CWic62m71Lupgt8FLspx9XkHsE7d6ZHdaovxIrJVRI4B+3FKx7nPVxSdcEosSW5VWTJOSaMW0OwcjguwGvgBWCcik0XkPrfaJue5n8o6r3vu8TgJrgHO1XomzpUzAOqsxFXYCldVybvUBTAWZ1K5Hjglob+ranZDvDizpj4vzrrbR9yYYjnzMy7o7y3LCzglt26a9wpiKTgJr9jW8K1sLBFUYKp6UlVnqeo/VPVinKmtRxVSBZKW67FStL+TrKSQ3ZArztoK+XLrez8DZuIkrg7A33CWvMzpuJcxzAe2AHeIMy32MNzlE11fA2HAH4Au7vnSgfw+F2/ekx9OQ3VMrlsL4K0CYk0CQvLYHoqTSFCnHvwK97YGp5prs4i0z3Huv+c6bzugOc6KYVmKOsPkIZxElpfDqrpFVX/GKX2OFJFeOfa/gFOF9384ySIGJxHl/oy9+XubhZPQ8puBtzZwUn9fLtUUkb+vAzClagPO7zwIp0hfVHHAeSJynqpmXU3Gcvo/btYXT3iObTGFHPcSYLeq/jNrg4g0OYv4AGcNZRF5H6fHSxzOle3H7nHr4JRs7lfVue62jhT8v+DNe1qBk3AOadHW6d2Ec0WfnajEWUi9PU51SPZ7wqni+1lE/gGsx6kyWe2eO0pVt+R1AhGJw/kddQYWudsa45QYC7ISCBORuqp6KL8nqeoRcTodvCQiHdxYuwHj3JIpbltVM5zqnaKagTNN9CQRUVX9KNf+aJzPwJwlKxFUQCJSR0TmiMgtItJORCJFZDDOOr2z9eyXvZyF88X1kYi0d6/kx+JcTTsVuqopOHXBj4tIGxG5GOfqsCC/Ag1F5GYRaSoi93HuvUA+Auq6556qqofd7UdwrnTvFpELRKQH8Kb7HvLk5Xv6FKeKaZqI9HA/80tF5MVCeg6NxSm5PCBOj6wYnHnma7s/EZGubu+aC90v8AE4ffuzFn35B3CTiPxDRKJFJMqtanvejX8T8B3wlohc5J7jQ5wqlYKsxFmDulshzwP4H9ASpxQAzu/0enHGdrTFabc566obdda3GAy8KWcObuuO8/7MWbJEUDEl43xx/QmnmmQ9Tg+S8ZzZm8ZrqpoJXI/TS2gpzpftv3GSQM665Kyuib/gVIv8rZDjTgfG4PTkWYPT0+Tps43TPeYenCvJWuQYO+C+h6E4VSfrgNdxqi9OFXLIAt+Tqp7A6c2zDWft2zicz6cWTvLJL84JwO3ubRnOF1oDoLs6K36BU0V0CU6V1mbgReCfqvqJe4yZwNVAL5zfy1KctYJ35jjVCGA7MAenO+h4nMbafLlVUu8DNxf0PPe5B3BKXaPc3mWP4iSRH3HaqBa798+amwyG4CS028DpsovT/vVBQa81BbOFacw5ceupV+H0qClPC8IbL4hIPZySx4Wqut3X8eQmImOAEFX1xWC2CsPaCEyRiMj1OI22m3G6VY7l93pqU8Go6gFxBp81xilRlDUHKLzq0RTCSgSmSNwi+d9w6qiP4IxFeMTLbp3GmDLIEoExxlRy1lhsjDGVnCUCY4yp5CwRGGNMJWeJwBhjKjlLBMYYU8n9PykMpf4xjn9+AAAAAElFTkSuQmCC",
      "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",
    "\n",
    "# Calculate the Frobenius-norm error\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",
    "# Plot\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",
    "### Case 2: Image compression\n",
    "\n",
    "In order to fulfill image processing tasks, we first import the necessary package.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:14.486390Z",
     "start_time": "2021-03-09T03:47:14.466171Z"
    }
   },
   "outputs": [],
   "source": [
    "# Image processing package PIL\n",
    "from PIL import Image\n",
    "\n",
    "# Open the picture prepared in advance\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": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:15.837676Z",
     "start_time": "2021-03-09T03:47:14.960968Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsUlEQVR4nO3dbWxcVXoH8P8fx46NHSW2A4mdOC+khGAqNomsQEVAdOkugS8EtOLlA2Il1KwqkIq0+wFRbZdWVctWBUSliioUtNmKQtgFRFqhdilaKUVasWvSxCREJCE4L7bjJMRJnMRxYufph7nZOtF9jsczd2Zsn/9Psjw+z5yZx9fz+I7v8TmHZgYRmf6uqXQCIlIeKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFil6KR/D7JUZJnxnzcXem85EozKp2ATBu/MbO1lU5CfDqzT3Mku0n+iGQXyVMkN5OsrXReUn4q9jg8DGAdgKUAbgXw/bQ7kVxL8mTgI3TmXkXyOMk9JH9MUu8aJxn9QOLwj2bWCwAk/x3AyrQ7mdknAOYU8PhbAfwhgAMAbgGwGcAIgL8r4LGkRHRmj8ORMbfPAWjI8sHNbL+ZfW1ml8zscwB/DeB7WT6HFE/FLr9H8s6rrqhf/XFnng9lAFjKXGXi9DZefs/M/gcFnPVJ3gdgm5n1k1wB4McAfpF1flIcndklC/cA6CJ5FsCHAN4D8LeVTUmuRi1eIRIHndlFIqFiF4mEil0kEip2kUiUdeitsbHRWltby/mUmbrmmvTfjV47AJD+cHPo4ujo6GhBsUIuuBaafyh26dKlCbUD4dxnzPBfqlVVVW7MU2gek11vby8GBgZSfzBFFTvJdQBeAVAF4F/M7IXQ/VtbW7F58+ZinrKiZs6cmdpeX1/v9qmpqXFjFy5ccGOnTp1yY6dPn3Zjw8PDqe2hgq6rq3Nj3vcMhAtwaGgotX1wcNDtMzIy4saam5sLinmFe+7cObfPxYsX3dhk98gjj7ixgt/Gk6wC8E8A7gPQDuAxku2FPp6IlFYxf7OvAbAv+b/oCwDeBvBANmmJSNaKKfYFAA6N+fpw0nYFkhtIdpLsHBgYKOLpRKQYJb8ab2YbzazDzDoaGxtL/XQi4iim2HsAtI35emHSJiKTUDHF/jsAN5JcSrIGwKMAtmSTlohkreChNzMbIfk0gP9CbujtDTPblVlmIpKposbZzexD5KY0isgkp3+XFYmEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4mEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEUTvCkOwGMAhgFMCImXVkkZSIZK+oYk/8sZkdz+BxRKSE9DZeJBLFFrsB+BXJz0huSLsDyQ0kO0l2DgwMFPl0IlKoYot9rZmtBnAfgKdI3nX1Hcxso5l1mFlHY2NjkU8nIoUqqtjNrCf5fBTA+wDWZJGUiGSv4GInWU9y1uXbAL4LYGdWiYlItoq5Gj8PwPskLz/Ov5nZf2aSlYhkruBiN7P9AL6VYS4iUkIaehOJhIpdJBIqdpFIqNhFIpHF/8ZH49SpU6ntPT09bp/z58+7sZkzZ7qx2bNnu7Frr73WjY2Ojqa29/f3u33OnDkz4ccDgAsXLrixwcHBCedx+PBhN1ZVVeXG2tvb3diKFStS2xcuXOj2qaurc2NTmc7sIpFQsYtEQsUuEgkVu0gkVOwikSj71XgzK/dTZubIkSOp7bt27XL7dHd3u7Ha2lo31tTUVFA/78p6b2+v2+fYsWNuzBuBAIChoSE35uWYzKVIdfDgQTe2f/9+N7ZkyRI3tn79+tT2e++91+2zYMECNzaV6cwuEgkVu0gkVOwikVCxi0RCxS4SCRW7SCQ0EWYCzp07l9p+6NAht09XV5cbO336tBsLTQoJDYd5Q17XXXed26ehocGNhSbJjIyMuLGlS5dOqH28PELHOBTzJt4MDw+7faYrndlFIqFiF4mEil0kEip2kUio2EUioWIXiURZh95IoqamppxPmam2trbU9jvvvNPts2jRIjfmzaIDgJ07/Z20+vr63Nj8+fNT2++44w63z2233TbhxwOAa67xzxXe8OCePXvcPlu3bnVjt9xyixsLDeetXbs2tb2lpcXtM5Vfo6FZheOe2Um+QfIoyZ1j2ppIfkRyb/JZ27OKTHL5vI3/GYB1V7U9C+BjM7sRwMfJ1yIyiY1b7Ga2FcCJq5ofALApub0JwPps0xKRrBV6gW6emV3+w/EIcju6piK5gWQnyc6BgYECn05EilX01XjLrTPlrjVlZhvNrMPMOhob9ae9SKUUWuz9JFsAIPl8NLuURKQUCh162wLgCQAvJJ8/yLfjVF5w0tt2KbRA4Zw5c9zY6tWr3dhDDz2Ud15jedtNXbx40e0TynHu3LlurL6+3o1dunQptf3Eiasv//y/vXv3urHQDLvQ8Oa8eel/YVZXV7t9pvJrNCSfobe3APwGwE0kD5N8Erki/w7JvQD+JPlaRCaxcc/sZvaYE7on41xEpIT077IikVCxi0RCxS4SCRW7SCS019sEeDOKQsM4M2fOdGOhPdtCC0SGZmWdPXs2tT20SGUoxxkz/JfI6OioG/P2gQstshmKeUN5QHj2XSjmmcqv0RCd2UUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIbeMhAa3gkNXVVVVbmx0FBTSF1dXWp7KMdQHqHhtdBiJL29vanthexTB/gzDoHwrL3QsKhnOr5GAZ3ZRaKhYheJhIpdJBIqdpFIqNhFIlHWq/FmVvBV5sksdDU7NGkl1C+05lroGHpX42fNmuX2CV1xD109P3DggBvbt2/fhB8vtN5daJ0/b505wD8eIVP5NRoaSdCZXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIaCLMBHi5e2vTAeHhtVC/0HZNoWE5b+JHQ0OD28fbMgoADh8+7Ma++uorN7Z///7U9tA6c83NzW5syZIlBfULra/nmcqv0ZB8tn96g+RRkjvHtD1Psofk9uTj/tKmKSLFyudt/M8ArEtpf9nMViYfH2ablohkbdxiN7OtAPytN0VkSijmAt3TJLuSt/mN3p1IbiDZSbIztNiBiJRWocX+KoBlAFYC6APwondHM9toZh1m1tHY6P5OEJESK6jYzazfzEbN7BKA1wCsyTYtEclaQUNvJFvMrC/58kEAO0P3ny68obLQ8FooFhriCQ2HhdaT84blQn2Gh4fdWE9Pjxv74osv3Jg39BZaS2758uVubNWqVW4sNFvOm3VYyLZQU924xU7yLQB3A5hL8jCAnwC4m+RKAAagG8APSpeiiGRh3GI3s8dSml8vQS4iUkLxvZcRiZSKXSQSKnaRSKjYRSJR9llvoZlek503VBYaQgt9v6F+oZlthSxGGVpEcWhoyI0dP37cjYWG5fr7+1Pb29ra3D5NTU1u7Prrr3djoUUlve879HOZyq/REJ3ZRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEWYfeSE7p2UaFDL2FhI5FocfJG0Y7duyY2+fQoUNurLu7240dOXLEjV24cCG1PbTnXGtrqxsL7ecWGooMDSt6pvJrNDRsOHW/KxGZEBW7SCRU7CKRULGLRELFLhIJTYTJQKETYUJXfUPbFs2Y4f/Yzpw5M6F2APjyyy/dWOhqfGgrp9ra2tT20ISW+fPnu7HZs2cXlEdoLT/PdHyNAjqzi0RDxS4SCRW7SCRU7CKRULGLRELFLhKJfHaEaQPwcwDzkNsBZqOZvUKyCcBmAEuQ2xXmYTOLcpvW0NBbKBbaGsrbtggIDw2dPHkytf3rr792++zZs8eNhSa7hIYAvSG2lpYWt09okkzI6OioGytkDbrpKp8z+wiAH5pZO4DbATxFsh3AswA+NrMbAXycfC0ik9S4xW5mfWa2Lbk9CGA3gAUAHgCwKbnbJgDrS5SjiGRgQn+zk1wCYBWATwHMG7OT6xHk3uaLyCSVd7GTbADwLoBnzOyK/0+03B+mqX+cktxAspNk54kTJ4pKVkQKl1exk6xGrtDfNLP3kuZ+ki1JvAXA0bS+ZrbRzDrMrCO0CYCIlNa4xc7cZcvXAew2s5fGhLYAeCK5/QSAD7JPT0Syks+stzsAPA7gc5Lbk7bnALwA4B2STwI4AODhkmQ4iXjDaKGtlUJCQ1ehWV7e+m6AP0utq6vL7bNt2zY3dvRo6hs2AOHtmtrb21Pbb7rpJrdPQ0ODGwvNXhseHnZj3rBcaNhzuhq32M3sEwDeoOQ92aYjIqWi/6ATiYSKXSQSKnaRSKjYRSKhYheJRNkXnCx0q6TJwBtiCw29hb7fUL/Q0NDFixfdWE9PT2r7jh073D6hBSdDM9FWrFjhxtauXZvafuutt7p9QotserP5AGBwcNCNeUKLfU7l12iIzuwikVCxi0RCxS4SCRW7SCRU7CKRULGLRKKsQ29mVvAMscnAW6QwtHhhKBYa4hkZGXFjQ0NDbswbhgoNT4XyaG5udmPLly93Y97stnnz/AWNQnu2hfaqC80CrK6uTm0vdEh0sgt9Xzqzi0RCxS4SCRW7SCRU7CKRULGLRKLsE2Gmo0InVYSuIn/zzTdurK+vz42dOnUqtb22ttbt09ra6sZuuOEGN7Zs2bIJP2Z9fb3bx8sdCE/+CfF+Ntr+SUSmLRW7SCRU7CKRULGLRELFLhIJFbtIJMYdeiPZBuDnyG3JbAA2mtkrJJ8H8KcAjiV3fc7MPixVopNBIRNhvO2HgPDEj+PHj7sxb4snAOjv709tb2xsdPssXrzYjd18881urK2tzY15a9eF1tYrdEJRiDf0Fhouna7yGWcfAfBDM9tGchaAz0h+lMReNrN/KF16IpKVfPZ66wPQl9weJLkbwIJSJyYi2ZrQexmSSwCsAvBp0vQ0yS6Sb5D03yeKSMXlXewkGwC8C+AZMzsN4FUAywCsRO7M/6LTbwPJTpKdAwMDxWcsIgXJq9hJViNX6G+a2XsAYGb9ZjZqZpcAvAZgTVpfM9toZh1m1hG6SCQipTVusTN3GfR1ALvN7KUx7S1j7vYggJ3ZpyciWcnnavwdAB4H8DnJ7UnbcwAeI7kSueG4bgA/KEF+k4o3XFPoemahLY12797txvbu3evGvJl0CxcudPuE1pJbtGiRG2toaHBj58+fT20PDUWGhsNCs/ZCM+JmzEh/icc46y2fq/GfAEg7MtN6TF1kuonvPwtEIqViF4mEil0kEip2kUio2EUioQUnSyw0LOcNTwHAoUOH3NjBgwfdWFNTU2p7aIZae3u7G5szZ44bC/H+W9IbChtPaOgtNIzmzbILDfNN5e2fQnRmF4mEil0kEip2kUio2EUioWIXiYSKXSQSGnrLQGgmV2g/t8HBQTfW29vrxg4cOODGvGGo6upqt8/s2bPdWF1dnRsbGRlxY96wYmgILRQLLVQZGirzYjEuOBnfdywSKRW7SCRU7CKRULGLRELFLhIJFbtIJDT0NgHeME5o6G14eNiNnT171o0dO3asoNjcuXNT20OzzUILR4b6hb43bxHI0JBXTU2NGwsNr4UWnPSGIjX0JiLTlopdJBIqdpFIqNhFIqFiF4nEuFfjSdYC2ApgZnL/X5rZT0guBfA2gGYAnwF43Mz8WR+RKsU2Q6FJId5kktAV9/r6ejcWWkMvNAoxNDTkxjyFXo2PcSunQuRzZh8G8G0z+xZy2zOvI3k7gJ8CeNnM/gDAAIAnS5aliBRt3GK3nDPJl9XJhwH4NoBfJu2bAKwvRYIiko1892evSnZwPQrgIwBfAThpZpcnNB8GsKAkGYpIJvIqdjMbNbOVABYCWANgRb5PQHIDyU6Snd5a4iJSehO6Gm9mJwH8GsAfAZhD8vIFvoUAepw+G82sw8w6Ghsbi8lVRIowbrGTvI7knOR2HYDvANiNXNF/L7nbEwA+KFGOIpKBfCbCtADYRLIKuV8O75jZf5D8AsDbJP8GwP8CeL2EeU4K3uSJ0NppoeGp5uZmN7Z48eKCHnPRokWp7aFtnELr04WeKzTk5Q2VhSbPhGKhCTmhdfK8YcpQ7qHveSobt9jNrAvAqpT2/cj9/S4iU4D+g04kEip2kUio2EUioWIXiYSKXSQSDM1qyvzJyGMALu9dNBfA8bI9uU95XEl5XGmq5bHYzK5LC5S12K94YrLTzDoq8uTKQ3lEmIfexotEQsUuEolKFvvGCj73WMrjSsrjStMmj4r9zS4i5aW38SKRULGLRKIixU5yHckvSe4j+Wwlckjy6Cb5OcntJDvL+LxvkDxKcueYtiaSH5Hcm3wu+UofTh7Pk+xJjsl2kveXIY82kr8m+QXJXST/PGkv6zEJ5FHWY0KyluRvSe5I8virpH0pyU+TutlM0l+ON42ZlfUDQBVya9jdAKAGwA4A7eXOI8mlG8DcCjzvXQBWA9g5pu3vATyb3H4WwE8rlMfzAH5U5uPRAmB1cnsWgD0A2st9TAJ5lPWYACCAhuR2NYBPAdwO4B0Ajybt/wzgzybyuJU4s68BsM/M9ltunfm3ATxQgTwqxsy2AjhxVfMDyK3SC5RptV4nj7Izsz4z25bcHkRuJaQFKPMxCeRRVpaT+YrOlSj2BQAOjfm6kivTGoBfkfyM5IYK5XDZPDPrS24fATCvgrk8TbIreZtf1oUDSS5BbrGUT1HBY3JVHkCZj0kpVnSO/QLdWjNbDeA+AE+RvKvSCQG53+zI/SKqhFcBLENuQ5A+AC+W64lJNgB4F8AzZnZ6bKycxyQlj7IfEytiRWdPJYq9B0DbmK/dlWlLzcx6ks9HAbyPyi6z1U+yBQCSz0crkYSZ9ScvtEsAXkOZjgnJauQK7E0zey9pLvsxScujUsckee6TmOCKzp5KFPvvANyYXFmsAfAogC3lToJkPclZl28D+C6AneFeJbUFuVV6gQqu1nu5uBIPogzHhLnVH18HsNvMXhoTKusx8fIo9zEp2YrO5brCeNXVxvuRu9L5FYC/qFAONyA3ErADwK5y5gHgLeTeDl5E7m+vJ5HbIPNjAHsB/DeApgrl8a8APgfQhVyxtZQhj7XIvUXvArA9+bi/3MckkEdZjwmAW5FbsbkLuV8sfznmNftbAPsA/ALAzIk8rv5dViQSsV+gE4mGil0kEip2kUio2EUioWIXiYSKXSQSKnaRSPwffyIxG32ku9UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUkklEQVR4nO3df4xV5Z3H8fe3yK/KYBGmgIgDAroCRSATUiq2rm4t2jTaZNNWN103caXp1mSbdP8w3WTrbvaPdrNt081u3NDV1DZdqbaauhvaqtQGTS064AhYrfwoWJAfQwUB5YcD3/3jHrID3u8zM+f+muH5vJLJ3Hm+99zzzJn7nXPv+d7neczdEZHz3/ta3QERaQ4lu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULLLoJnZfDP7hZkdMLP3fFDDzC42s8fM7G0z22lmt7ein3I2JbuU8S7wMHBnEP8P4CQwGfgL4D4zm9ekvknA9Am684uZ7QD+HfhLoAP4OXCHux9vwL5mA1vc3fq0XQgcBOa7+2tF2w+A3e5+T737IAOnM/v56TPAcmAmsAD4q2p3MrNlZnYo8bWsxL6vAHrPJHrhJUBn9ha7oNUdkIb4N3d/A8DM/gdYWO1O7v4s8IE673sccPictreAtjrvRwZJZ/bz094+t9+hkoDNchQYf07beOBIE/sgVSjZM2Zm15rZ0cTXtSUe9jXgAjOb06ftauDl+vRaytLL+Iy5+zOUOOubmQGjgVHFz2MqD+cn3P1tM3sU+Ccz+2sqbyFuAT5St45LKTqzSxkdwDH+/2x9DPhdn/jfAGOB/cBDwBfdXWf2FlPpTSQTOrOLZELJLpIJJbtIJpTsIploault0qRJ3tHR0cxdimRl586dHDhwwKrFakp2M1sOfAcYAfyXu389df+Ojg5+/etf17LLIalSdq4uVe1Ixco+ZpnHSynbx2Y9Xn+PWe99DXUf+Uj8cYbSL+PNbASVoYw3AXOB28xsbtnHE5HGquU9+xJgq7tvd/eTwCoqn5QSkSGolmSfBvyhz8+7irazmNkKM+sys66enp4adicitWj41Xh3X+nune7e2d7e3ujdiUiglmTfDUzv8/OlRZuIDEG1JPsLwBwzm2lmo4DPAY/Xp1siUm+lS2/u3mtmdwO/oFJ6e0Ajm0SGrprq7O6+Glhdp76ISAPp47IimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimahpRRgz2wEcAU4Bve7eWY9OiUj91ZTshT919wN1eBwRaSC9jBfJRK3J7sATZrbezFZUu4OZrTCzLjPr6unpqXF3IlJWrcm+zN0XAzcBXzKzj557B3df6e6d7t7Z3t5e4+5EpKyakt3ddxff9wOPAUvq0SkRqb/SyW5mF5pZ25nbwI3A5np1TETqq5ar8ZOBx8zszOP8t7v/vC69EpG6K53s7r4duLqOfRGRBlLpTSQTSnaRTCjZRTKhZBfJRD0+G5+Nt956a1DtAO97X/z/dObMmWHs+PHjYWzEiBFh7MiRI1Xbd+3aFW5z+vTpMHby5Mkwtn379jB29OjRqu2p45E6jtu2bQtjF1wQP43nz59ftX3ZsmXhNrNmzQpjw5nO7CKZULKLZELJLpIJJbtIJpTsIplo+tV4d2/2Lutmx44dVdufe+65cJt9+/aFsblz55bqx6lTpwa9v9///vfhNqm/yeHDh8PY66+/HsaKMRPvceGFF4bbpK7Gd3d3h7He3t4wduutt1Ztv+KKK8JtLr/88jA2nOnMLpIJJbtIJpTsIplQsotkQskukgklu0gmNBBmEKLS0Pr168Ntfvazn4Wxt99+O4yNGTMmjKVKb9GgkPHjx5faV2pAzrhx48LYpZdeWrV9ypQp4TYTJkwIY1HZE9Ilu2hG44kTJ4bbnK90ZhfJhJJdJBNKdpFMKNlFMqFkF8mEkl0kE00vvaXmIBvqFi9eXLU9NZIrNbpq3bp1YSw12iy1v6iP119/fbjN7Nmzw1iqZJdalfeiiy6q2j569Ohwm9WrV4exX/3qV2Fs3rx5YeyGG26o2p6a/284P0dT+v2tzOwBM9tvZpv7tF1sZk+a2Zbie1wgFZEhYSD/wr4HLD+n7R5gjbvPAdYUP4vIENZvsrv7WuDNc5pvAR4sbj8I3FrfbolIvZV9czLZ3fcUt/dSWdG1KjNbYWZdZtZ14MCBkrsTkVrVfCXCK3MahfMauftKd+90985JkybVujsRKalssu8zs6kAxff99euSiDRC2dLb48AdwNeL7z8d6IbDecLJqJy0aNGicJtUiee2224LY2XLP1FpKzVCLVUOS42wmz59ehiLlqhKTRy5atWqMBYtJwWwdOnSMHbZZZeFschwfo6mDKT09hDwHHClme0yszupJPnHzWwL8GfFzyIyhPV7Znf36PRT/dMKIjIknZ8fFRKR91Cyi2RCyS6SCSW7SCa01tsgROWw1ISNqUkU29rawtjYsWPD2MmTJ8NYdHxHjhwZbpNy4sSJMJYafRdNpvnqq6+G22zYsCGMRWvHQbq8+cEPfrBq+6hRo8JtTp8+HcaGM53ZRTKhZBfJhJJdJBNKdpFMKNlFMqFkF8mE1nobhN7e3qrtqbJQKhaNDIPyJcrU/iKpUlNq9F2qj9u2bava/vzzz4fbHDp0KIwtXLgwjM2dOzeMRaXP1O+l0puIDGtKdpFMKNlFMqFkF8mEkl0kE02/Gl/mavFQEV2NT129TQ1AueCC+PC/++67YSx1DKOrzKk+puaZSw3IiQa7AKxfv75q+7PPPhtukzpWN910Uxj70Ic+FMai/qeOx3B+jqbozC6SCSW7SCaU7CKZULKLZELJLpIJJbtIJjQQZhCiElWqjJMqr6UGY6TKYan506LHTJXyUlKDdfbs2RPGotLb1q1bw21SSzUtWbIkjLW3t4exSFRGhfTvPJwNZPmnB8xsv5lt7tN2r5ntNrPu4uvmxnZTRGo1kJfx3wOWV2n/trsvLL5W17dbIlJv/Sa7u68F3mxCX0SkgWq5QHe3mW0sXuaHk6Ob2Qoz6zKzrp6enhp2JyK1KJvs9wGzgIXAHuCb0R3dfaW7d7p7Z5kLKSJSH6WS3d33ufspdz8NfBeIL5WKyJBQqvRmZlPd/Uzd5dPA5tT9zxdlllBKlddS5bDUdqlRWVEZMFXKSy1fdezYsTD29NNPh7ForrmJEyeG23z2s58NYwsWLAhjKWXmDTxf9ZvsZvYQcB0wycx2AV8DrjOzhYADO4AvNK6LIlIP/Sa7u99Wpfn+BvRFRBpIH5cVyYSSXSQTSnaRTCjZRTKhUW+DEI2GKjt6LTXyKjVaLlU2ivZXtpT3xz/+MYx1d3eHsR07dlRtv+SSS8JtUss4tbW1hbHUcYyWqEodj/NVfr+xSKaU7CKZULKLZELJLpIJJbtIJpTsIplQ6a0OUhNORqUfSJd/yq4DF0lNUnn06NEw9swzz4SxF154IYxFI+mWLl0abjNv3rwwlhpxePLkyTAWlRVzHPWmM7tIJpTsIplQsotkQskukgklu0gmdDV+EKKr7qkr7qlYI5aGKjNY54033ghjq1fH63+89tprYWz+/PlV26+77rpwm46OjjCWqnikrqxHv3dqm9TfbDjTmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTAxkRZjpwPeByVRWgFnp7t8xs4uBHwEzqKwK8xl3P9i4rrZeVJJJlWpSJZ5ULDWvWqqMFsUOHoz/NC+++GKpWKr/V1111aDaAUaPHh3GUstQpQbJRMfjfC2vpQzkzN4LfMXd5wIfBr5kZnOBe4A17j4HWFP8LCJDVL/J7u573H1DcfsI8AowDbgFeLC424PArQ3qo4jUwaDes5vZDGARsA6Y3Gcl171UXuaLyBA14GQ3s3HAT4Avu/vhvjGvvAGq+ibIzFaYWZeZdfX09NTUWREpb0DJbmYjqST6D9390aJ5n5lNLeJTgf3VtnX3le7e6e6d7e3t9eiziJTQb7Jb5ZLr/cAr7v6tPqHHgTuK23cAP61/90SkXgYy6u0a4PPAJjPrLtq+CnwdeNjM7gR2Ap9pSA+HgVQJKhqFBumRXO+8804YmzBhQhiLRsS9/PLL4TaPPPJIGHv11VfD2NVXXx3Grr322qrts2bNCrcpOy9cmRFsZUfRDWf9Jru7PwtEv/0N9e2OiDSKPkEnkgklu0gmlOwimVCyi2RCyS6SCU04OQipkWhlpMpyZZc7ikpsqfLamjVrwlhbW1sY+8QnPhHGrrnmmqrt48ePD7dJLUOVmpyzzFJZZUcqDmc6s4tkQskukgklu0gmlOwimVCyi2RCyS6SCZXeBiFVKisjtWZbavLFVAnw9ddfr9qemjjyyJEjYezGG28MYx/72MfC2CWXXFK1PTXaLKXs2nc5TiwZ0ZldJBNKdpFMKNlFMqFkF8mEkl0kE02/Gj+cr45GAy5SAydSV4pTsdRAmP37q07kC8CmTZuqtkdX6QFmzJgRxj71qU+FsdQcdKNGjaranppbL9oG0sf4xIkTYSz6m6Wu7g/n52iKzuwimVCyi2RCyS6SCSW7SCaU7CKZULKLZKLf0puZTQe+T2VJZgdWuvt3zOxe4C7gzNKsX3X31Y3q6PkoVf6J5k4D+M1vfhPGfvnLX1ZtT83vdvvtt4ex5cuXh7EpU6aEsagclhoIU2Yuuf5Ex7hsuXQ4G0idvRf4irtvMLM2YL2ZPVnEvu3u/9q47olIvQxkrbc9wJ7i9hEzewWY1uiOiUh9Deo9u5nNABYB64qmu81so5k9YGbx0qIi0nIDTnYzGwf8BPiyux8G7gNmAQupnPm/GWy3wsy6zKyrp6en2l1EpAkGlOxmNpJKov/Q3R8FcPd97n7K3U8D3wWWVNvW3Ve6e6e7d7a3t9er3yIySP0mu1UuW94PvOLu3+rTPrXP3T4NbK5/90SkXgZyNf4a4PPAJjPrLtq+CtxmZguplON2AF8YyA6H89I6Ufmn7Lxq73//+8PY3r17w9hTTz0Vxrq7u6u2X3nlleE2d911Vxi77LLLwtixY8fCWGTs2LFhLFVeS8VS8/VFz7fU32w4P0dTBnI1/lmg2m+vmrrIMKJP0IlkQskukgklu0gmlOwimVCyi2RCyz8NQlSSSS0LlZo4MjX54hNPPBHGNm+OP9IwbVr1YQuf/OQnw22uuuqqMFZ2dFg0eWTq8VLltdRxTI0ePHnyZNX21KSS9V7ma6jQmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTKj0NghlSm+pstCbb74ZxtauXRvGtmzZEsY6Ojqqtk+dOrVqO6TLUPUuUaVGm6WOVSpWdj29iEpvIjKsKdlFMqFkF8mEkl0kE0p2kUwo2UUy0fTSW6qUM9RF5Z9UOenw4cNhLDV6bcOGDWHs4MGDYWzBggVV21Olt9Ros2j0GqRLVGUm4Uyt9ZYqoaWeU6n+R8pOIDrU6cwukgklu0gmlOwimVCyi2RCyS6SiX6vxpvZGGAtMLq4/4/d/WtmNhNYBUwE1gOfd/fqE36d51JXilNX43fu3BnGUivetrW1hbFomafZs2eH25S9qp7arkzVpRGDXaLlplL9O3HiRBgbzgZyZj8BXO/uV1NZnnm5mX0Y+AbwbXefDRwE7mxYL0WkZv0mu1ccLX4cWXw5cD3w46L9QeDWRnRQROpjoOuzjyhWcN0PPAlsAw65e29xl11A9TmMRWRIGFCyu/spd18IXAosAf5koDswsxVm1mVmXan3oSLSWIO6Gu/uh4CngaXAB8zszAW+S4HdwTYr3b3T3Tvb29tr6auI1KDfZDezdjP7QHF7LPBx4BUqSf/nxd3uAH7aoD6KSB0MZCDMVOBBMxtB5Z/Dw+7+v2b2W2CVmf0z8CJw/0B2mFr+Z6iL+j5mzJhwm8mTJ4exOXPmhLEpU6aEsdSAkUWLFlVtnzlzZrhNyvHjx0v1IyrL9fb2Vm2HdAkttfxTqox29OjRqu1l5w0czvpNdnffCLznGeTu26m8fxeRYeD8/BcmIu+hZBfJhJJdJBNKdpFMKNlFMmHNnBPOzHqAM0O9JgEHmrbzmPpxNvXjbMOtHx3uXvXTa01N9rN2bNbl7p0t2bn6oX5k2A+9jBfJhJJdJBOtTPaVLdx3X+rH2dSPs503/WjZe3YRaS69jBfJhJJdJBMtSXYzW25mvzOzrWZ2Tyv6UPRjh5ltMrNuM+tq4n4fMLP9Zra5T9vFZvakmW0pvk9oUT/uNbPdxTHpNrObm9CP6Wb2tJn91sxeNrO/LdqbekwS/WjqMTGzMWb2vJm9VPTjH4v2mWa2rsibH5nZ4Bayc/emfgEjqMxhdzkwCngJmNvsfhR92QFMasF+PwosBjb3afsX4J7i9j3AN1rUj3uBv2vy8ZgKLC5utwGvAXObfUwS/WjqMQEMGFfcHgmsAz4MPAx8rmj/T+CLg3ncVpzZlwBb3X27V+aZXwXc0oJ+tIy7rwXePKf5Fiqz9EKTZusN+tF07r7H3TcUt49QmQlpGk0+Jol+NJVX1H1G51Yk+zTgD31+buXMtA48YWbrzWxFi/pwxmR331Pc3gvEU9w03t1mtrF4md/wtxN9mdkMKpOlrKOFx+ScfkCTj0kjZnTO/QLdMndfDNwEfMnMPtrqDkHlPzuVf0StcB8wi8qCIHuAbzZrx2Y2DvgJ8GV3P2spnWYekyr9aPox8RpmdI60Itl3A9P7/BzOTNto7r67+L4feIzWTrO1z8ymAhTf97eiE+6+r3iinQa+S5OOiZmNpJJgP3T3R4vmph+Tav1o1TEp9n2IQc7oHGlFsr8AzCmuLI4CPgc83uxOmNmFZtZ25jZwI7A5vVVDPU5lll5o4Wy9Z5Kr8GmacEysMpPn/cAr7v6tPqGmHpOoH80+Jg2b0blZVxjPudp4M5UrnduAv29RHy6nUgl4CXi5mf0AHqLycvBdKu+97qSyQOYaYAvwFHBxi/rxA2ATsJFKsk1tQj+WUXmJvhHoLr5ubvYxSfSjqccEWEBlxuaNVP6x/EOf5+zzwFbgEWD0YB5XH5cVyUTuF+hEsqFkF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQT/weYte4bQI3fgAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASZUlEQVR4nO3df4xdZZ3H8feH2h9AK6XbsU5K6ZS2AmWBgYxdQQS2VoP800o2YmOUjWRrNpBogskSN1nZxGR1s2rchLgpP2I1rsiKRlzrQiEutEJopy5Oi0VoaRFK2xki/QHdpbT97h/3NE7rec7M3J/tPJ9XMpk7z/eee74c+plz5zz3nKOIwMzGvzM63YCZtYfDbpYJh90sEw67WSYcdrNMOOxmmXDYzTLhsNuYSfpzSY9Iel3Sn3xQQ9J/S/o/SW8WX7/rRJ92Iofd6vEO8CBwa8Vzbo+IqcXXhW3qyyo47OOMpJ2SvihpQNJ+ST+UNKWZ64iI30XEfcBzzXxday2HfXz6BHADMA+4DPjrsidJukbSvoqvaxro4Z+Kt/m/knR9A69jTfKuTjdgLfGvEfEagKSfAb1lT4qI9cD0Fqz/74DfAoeBTwI/k9QbEdtbsC4bJe/Zx6c9wx4fAqa2c+UR8UxEHIyItyNiNfAr4MZ29mB/ymHPmKQPDTtiXvb1oSatKgA16bWsTn4bn7GIWEcde31JAiYDk4qfp9ReLt6WNB34C+AJ4AhwM3At8PkmtW11ctitHnOBHcN+/l/gZaAHmAh8BbgIOAo8DyyPiBfa3KOdRL54hVke/De7WSYcdrNMOOxmmXDYzTLR1qPxM2fOjJ6ennau0iwrO3fu5PXXXy/9TENDYZd0A/AtYAJwb0R8ter5PT09bNy4sZFVmlmF97///cla3W/jJU0A7gY+BiwCVkhaVO/rmVlrNfI3+2JgW0S8FBGHgQeAZc1py8yarZGwzwZeGfbzq8XYCSStlNQvqX9oaKiB1ZlZI1p+ND4iVkVEX0T0dXV1tXp1ZpbQSNh3AXOG/XxeMWZmp6BGwr4RWChpnqRJ1C5S8HBz2jKzZqt76i0ijki6HXiE2tTb/RHha5KZnaIammePiDXAmib1YmYt5I/LmmXCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2WioTvCSNoJHASOAkcioq8ZTZlZ8zUU9sJfRsTrTXgdM2shv403y0SjYQ/gUUmbJK0se4KklZL6JfUPDQ01uDozq1ejYb8mIq4EPgbcJunak58QEasioi8i+rq6uhpcnZnVq6GwR8Su4vsg8BNgcTOaMrPmqzvsks6WNO34Y+CjwJZmNWZmzdXI0fhZwE8kHX+df4+I/2pKV2bWdHWHPSJeAi5vYi9m1kKeejPLhMNulgmH3SwTDrtZJprx2fhsHDx4sHT8wIEDyWXOOCP9+7S7u7uuPooZkFJvvPFG6fiuXbvqWlfqvxnglVdeGfNyx44dSy5TtR137NiRrFW95mWXXVY6vmTJkuQyCxcuTNZOZ96zm2XCYTfLhMNulgmH3SwTDrtZJnw0fgy2b99eOr5u3brkMnv27EnWFi1a1HBPJ3vttddKx1O9j2T//v3JWtUR8tQsxNSpU5PL7Nu3L1nbtGlTshYRydqyZctKxy+++OLkMj4ab2anNYfdLBMOu1kmHHazTDjsZplw2M0y4am3MRgcHCwdr5p6e/TRR5O1w4cPN9zTySZPnlw6XjXlNXHixGSt6kSeqtrcuXPHNA4wffr0ZG3btm3J2qFDh5K1888/v3R81qxZyWXGK+/ZzTLhsJtlwmE3y4TDbpYJh90sEw67WSY89TYGS5cuLR2fM2dOcplLL700WVu/fn2yduTIkdE3NszixeW320ud/QUwf/78utZVdbbZWWedVTpeNd34xBNPJGtV2+q9731vsnbdddeVjr/vfe9LLjNejbhnl3S/pEFJW4aNzZC0VtKLxfdzW9ummTVqNG/jvwPccNLYncDjEbEQeLz42cxOYSOGPSKeBP5w0vAyYHXxeDWwvLltmVmz1XuAblZE7C4e76F2R9dSklZK6pfUPzQ0VOfqzKxRDR+Nj9pRmuSRmohYFRF9EdHX1dXV6OrMrE71hn2vpG6A4nv5GSJmdsqod+rtYeAW4KvF9582raNTWOosrwULFiSXue2225K1z372s3X1Uc+U19lnn51cZtKkSXX1USW1rQYGBpLL3Hvvvcla1cUoly9fnqylptiqbqE1Xo1m6u0HwNPAhZJelXQrtZB/RNKLwNLiZzM7hY24Z4+IFYnSh5vci5m1kD8ua5YJh90sEw67WSYcdrNM+Ky3Jqi6YOOMGTPa2MmpI3URyOeffz65TH9/f7L2zjvvJGuXXHJJsvae97yndLxq6q1qavN05j27WSYcdrNMOOxmmXDYzTLhsJtlwmE3y4Sn3k5D9UwbtXuqafv27aXjVffFO3DgQLLW29ubrF1++eXJ2rvf/e5kLTfes5tlwmE3y4TDbpYJh90sEw67WSZ8NP40VM/R81Ycca+6RVXqpJbHHnssuczkyZOTtZtuuilZu+iii5K1M888s3R8vJ7sUsV7drNMOOxmmXDYzTLhsJtlwmE3y4TDbpYJT71Z3fbs2ZOsbdiwoXT8hRdeSC5zwQUXJGtLlixJ1qZNm5as5TjFljKa2z/dL2lQ0pZhY3dJ2iXp2eLrxta2aWaNGs3b+O8AN5SMfzMieouvNc1ty8yabcSwR8STwB/a0IuZtVAjB+hulzRQvM0/N/UkSSsl9UvqHxoaamB1ZtaIesP+bWA+0AvsBr6eemJErIqIvojo6+rqqnN1ZtaousIeEXsj4mhEHAPuARY3ty0za7a6pt4kdUfE7uLHjwNbqp5v49OaNenjsuvXry8d7+7uTi7zqU99KlmrusVT1dly9kcjhl3SD4DrgZmSXgW+DFwvqRcIYCfwuda1aGbNMGLYI2JFyfB9LejFzFrIH5c1y4TDbpYJh90sEw67WSZ81ptVnhlWdUumjRs3Jms7duwoHe/p6UkuU3UbpylTpiRrNjres5tlwmE3y4TDbpYJh90sEw67WSYcdrNMeOptnJFUOl41vVZ1z7Zf/OIXyVrqfm6Qvsfa4sXps6GvuOKKZO2MM7xfapS3oFkmHHazTDjsZplw2M0y4bCbZcJH48eZ1FH3o0ePJpcZHBxM1h566KFkrepWTpdeemnp+NKlS5PLzJ07N1mzxnnPbpYJh90sEw67WSYcdrNMOOxmmXDYzTIxmjvCzAG+C8yidgeYVRHxLUkzgB8CPdTuCvOJiHijda3acamTXSA99fbWW28ll3nqqaeSteeeey5ZqzqBZtGiRaXjvb29yWWstUazZz8C3BERi4APALdJWgTcCTweEQuBx4ufzewUNWLYI2J3RPy6eHwQ2ArMBpYBq4unrQaWt6hHM2uCMf3NLqkHuAJ4Bpg17E6ue6i9zTezU9Sowy5pKvAQ8IWIOOFi4lH7Q7H0j0VJKyX1S+ofGhpqqFkzq9+owi5pIrWgfz8iflwM75XUXdS7gdIPWEfEqojoi4i+rq6uZvRsZnUYMeyqHfq9D9gaEd8YVnoYuKV4fAvw0+a3Z2bNMpqz3j4IfBrYLOnZYuxLwFeBByXdCrwMfKIlHdqfqLqeXGpa7ve//31ymbvvvjtZS93GCWD+/PnJ2nXXXVc6fuGFFyaXsdYaMewRsR5ITex+uLntmFmr+BN0Zplw2M0y4bCbZcJhN8uEw26WCV9wcgyqzjZLqZomq3ddVa85MDBQOn7PPfckl9m0aVOy9vbbbydrK1asSNauv/760vEJEyYkl2n3tsqN9+xmmXDYzTLhsJtlwmE3y4TDbpYJh90sE556G4N2TuPUc2YbpM9ue/rpp5PLHDp0KFm7+uqrk7XU9BrAeeedVzreim3o6bXR8Z7dLBMOu1kmHHazTDjsZplw2M0y4aPxLdaKkzT27t2brG3evLl0/OWXX04ukzpyDvCZz3wmWbv44ouTtXe9q/yflk926Rzv2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmRpx6kzQH+C61WzIHsCoiviXpLuBvgOO3Zv1SRKxpVaOnq3qnhY4dO5asrVu3Lllbs6b8f0HVteSWL1+erN10003JWtWNOlP/3fVOoXl6rXGjmWc/AtwREb+WNA3YJGltUftmRPxL69ozs2YZzb3edgO7i8cHJW0FZre6MTNrrjH9zS6pB7gCeKYYul3SgKT7JZ3b7ObMrHlGHXZJU4GHgC9ExAHg28B8oJfanv/rieVWSuqX1D80NFT2FDNrg1GFXdJEakH/fkT8GCAi9kbE0Yg4BtwDLC5bNiJWRURfRPRVHdAxs9YaMeyqHT69D9gaEd8YNt497GkfB7Y0vz0za5bRHI3/IPBpYLOkZ4uxLwErJPVSm47bCXyuBf2Na1XTSfv370/WHnnkkWQtdSunSy65JLnMHXfckayde64PxYwXozkavx4omxz1nLrZacSfoDPLhMNulgmH3SwTDrtZJhx2s0z4gpMtVnWW19GjR5O1tWvXJmsbN25M1s4///zS8Ztvvjm5zPz585O11IUjob5bVLX77LVTpY9TgffsZplw2M0y4bCbZcJhN8uEw26WCYfdLBOeemuxqgtHvvnmm8naz3/+82Rt+/btydqCBQtKx6vOXqt3eq1Ksy84Wa8cp9hSvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU2xhUTRulvPXWW8naU089laxt2LAhWauasjvnnHNKx2fNmpVcpt3TYc12uvffLt6zm2XCYTfLhMNulgmH3SwTDrtZJkY8Gi9pCvAkMLl4/o8i4suS5gEPAH8GbAI+HRGHW9lsp9VzZPfw4fQm2bp1a7K2b9++ZG3q1KnJ2sKFC0vHq27/1Ioj1u289puPuI/OaPbsbwNLIuJyardnvkHSB4CvAd+MiAXAG8CtLevSzBo2Ytij5vjE7sTiK4AlwI+K8dXA8lY0aGbNMdr7s08o7uA6CKwFtgP7IuJI8ZRXgdkt6dDMmmJUYY+IoxHRC5wHLAYuGu0KJK2U1C+pf2hoqL4uzaxhYzoaHxH7gF8CVwHTJR0/wHcesCuxzKqI6IuIvq6urkZ6NbMGjBh2SV2SphePzwQ+AmylFvq/Kp52C/DTFvVoZk0wmhNhuoHVkiZQ++XwYET8p6TfAg9I+grwP8B9LezztHXWWWcla1dffXWyNnPmzGRt9uz04ZGrrrqqdHzevHnJZeqdumr2CSitOKHFt3/6oxHDHhEDwBUl4y9R+/vdzE4D/gSdWSYcdrNMOOxmmXDYzTLhsJtlQu2cgpA0BLxc/DgTeL1tK09zHydyHyc63fqYGxGln15ra9hPWLHUHxF9HVm5+3AfGfbht/FmmXDYzTLRybCv6uC6h3MfJ3IfJxo3fXTsb3Yzay+/jTfLhMNulomOhF3SDZJ+J2mbpDs70UPRx05JmyU9K6m/jeu9X9KgpC3DxmZIWivpxeL7uR3q4y5Ju4pt8qykG9vQxxxJv5T0W0nPSfp8Md7WbVLRR1u3iaQpkjZI+k3Rxz8W4/MkPVPk5oeSJo3phSOirV/ABGrXsLsAmAT8BljU7j6KXnYCMzuw3muBK4Etw8b+GbizeHwn8LUO9XEX8MU2b49u4Mri8TTgBWBRu7dJRR9t3SaAgKnF44nAM8AHgAeBTxbj/wb87VhetxN79sXAtoh4KWrXmX8AWNaBPjomIp4E/nDS8DJqV+mFNl2tN9FH20XE7oj4dfH4ILUrIc2mzdukoo+2ipqmX9G5E2GfDbwy7OdOXpk2gEclbZK0skM9HDcrInYXj/cA6Xsst97tkgaKt/kt/3NiOEk91C6W8gwd3CYn9QFt3iatuKJz7gforomIK4GPAbdJurbTDUHtNzu1X0Sd8G1gPrUbguwGvt6uFUuaCjwEfCEiDgyvtXOblPTR9m0SDVzROaUTYd8FzBn2c/LKtK0WEbuK74PAT+jsZbb2SuoGKL4PdqKJiNhb/EM7BtxDm7aJpInUAvb9iPhxMdz2bVLWR6e2SbHufYzxis4pnQj7RmBhcWRxEvBJ4OF2NyHpbEnTjj8GPgpsqV6qpR6mdpVe6ODVeo+Hq/Bx2rBNVLsq5H3A1oj4xrBSW7dJqo92b5OWXdG5XUcYTzraeCO1I53bgb/vUA8XUJsJ+A3wXDv7AH5A7e3gO9T+9rqV2g0yHwdeBB4DZnSoj+8Bm4EBamHrbkMf11B7iz4APFt83djubVLRR1u3CXAZtSs2D1D7xfIPw/7NbgC2Af8BTB7L6/rjsmaZyP0AnVk2HHazTDjsZplw2M0y4bCbZcJhN8uEw26Wif8HhJAtV2Fb+ucAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Then we look at the effect of the classic singular value decomposition\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": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:15.858695Z",
     "start_time": "2021-03-09T03:47:15.847413Z"
    }
   },
   "outputs": [],
   "source": [
    "# Hyper-parameters\n",
    "N = 5           # Number of qubits\n",
    "T = 8           # Set the number of rank you want to learn\n",
    "ITR = 200       # Number of iterations\n",
    "LR = 0.02       # Learning rate\n",
    "SEED = 14       # Random number seed\n",
    "\n",
    "# Set the learning weight\n",
    "weight = np.arange(2 * T, 0, -2).astype('complex128')\n",
    "\n",
    "# Convert the image into numpy array\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",
    "# Set circuit parameters\n",
    "cir_depth = 40 # Circuit depth\n",
    "block_len = 1 # The length of each module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:16.002083Z",
     "start_time": "2021-03-09T03:47:15.993385Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define quantum neural network\n",
    "def U_theta():\n",
    "\n",
    "    # Initialize the network with Circuit\n",
    "    cir = Circuit(N)\n",
    "    \n",
    "    # Build a hierarchy：\n",
    "    for _ in range(cir_depth):\n",
    "        cir.ry()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:53:07.440520Z",
     "start_time": "2021-03-09T03:47:21.094099Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -5052.6035\n",
      "iter: 10 loss: -108312.2634\n",
      "iter: 20 loss: -127067.9950\n",
      "iter: 30 loss: -138785.3241\n",
      "iter: 40 loss: -144992.7031\n",
      "iter: 50 loss: -148297.3650\n",
      "iter: 60 loss: -150532.9232\n",
      "iter: 70 loss: -152235.0095\n",
      "iter: 80 loss: -153425.8656\n",
      "iter: 90 loss: -154247.9781\n",
      "iter: 100 loss: -154813.1122\n",
      "iter: 110 loss: -155224.0510\n",
      "iter: 120 loss: -155530.3507\n",
      "iter: 130 loss: -155773.2695\n",
      "iter: 140 loss: -155976.5740\n",
      "iter: 150 loss: -156153.5143\n",
      "iter: 160 loss: -156310.1349\n",
      "iter: 170 loss: -156450.4223\n",
      "iter: 180 loss: -156576.0634\n",
      "iter: 190 loss: -156689.3533\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x297495eb6d0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYRklEQVR4nO2da2xdVXbH/ysmbzt2TF5OHEiAQASEycNECYkQHTQjQCMBUoXgA+IDmoyqQSrS9AOiUqFSPzBVAfGhogolmkxFeXQCIlSonRQNgpGAwYG8UzKQB8Fx4rwTQxJIsvrhnqhOdNbf1+fa5xr2/ydFud7L++x19tnL13f/vdY2d4cQ4ofPiHo7IIQoBwW7EImgYBciERTsQiSCgl2IRFCwC5EIl9XS2czuAPAcgAYA/+ruT7Hvb21t9fb29gGP09DQkNvOZMNz584NeBwAGDEi/vkXjXfZZfE0Mh/NLLSdPXs2tJ0/f37A14zmsL/rMT+Y/9F4bH7ZM2NjsTlm913EjyLroz9bBLvnyLZ3714cOXIk11g42M2sAcA/A/gJgK8AfGxma919W9Snvb0db731Vq6NTWJjY2Nu+3fffRf2+frrr0MbY+zYsaEtGq+lpSXswwJp5MiRoe3QoUOh7ZtvvhnwNSdMmBD2OXPmTGjr6ekJbaNGjQpt0TMbN25c2Ke3t7fQWKdPnw5tzc3Noa2IH0XWB8B9jN4s2PqI4uXOO++M+4SW/lkM4HN33+nu3wJ4BcDdNVxPCDGE1BLsMwDs7fP1V1mbEGIYMuQbdGa2wsw6zazzyJEjQz2cECKglmDvAjCzz9ftWdtFuPtKd+9w947W1tYahhNC1EItwf4xgDlmNtvMRgG4H8DawXFLCDHYFN6Nd/ezZvYIgP9GRXpb5e5bq+iX287kjmhnnckqTA4rKnlF47FdWLYbX7Qf8z+65tGjR8M+TAlhu+dF5CT2nNmO+4kTJwr1Y/MYweb31KlToY3tno8ePTq0RXNSRH5l91uTzu7ubwN4u5ZrCCHKQX9BJ0QiKNiFSAQFuxCJoGAXIhEU7EIkQk278QPl/PnzNOlioDAZh43DpCYmu0TSCpNcmDxYJKsJ4FJT1I9JkYyi2YORHywhhNnYPI4ZMya0Rc+GJRMxeY35waQyNo+RhFkkw47JoXpnFyIRFOxCJIKCXYhEULALkQgKdiESodTdeDMLd0eL1EFjyQXMVnRHNdrRZgkt7HpFa+ixa3777be57WyHme3uF0ngAIDx48fntrOdc/Zc9u/fH9qY0hCtN7bTzeaDPWuWNMTWdzSPRXbwGXpnFyIRFOxCJIKCXYhEULALkQgKdiESQcEuRCKUngjDEhAioqq0LHGCSRORPAVw2SWSa4oeacROHmGJMEziicZj12N+sJNkLr/88tAW1Yw7duxY2IfBxiqSCFP0qKmi9emamppCW3T6D1uLUaKXEmGEEAp2IVJBwS5EIijYhUgEBbsQiaBgFyIRapLezGw3gJMAzgE46+4d/Xx/KIWwjKcoY4hJE0xeY2MxGS3KeGKZYUwKYbIW85/ZovltaWkJ+7Asuujorf76NTY25rZPmzYt7MOOeGInADMfI8mOSWhMlmMZdmztMMk5Wj9F1im7r8HQ2f/C3fOFQiHEsEG/xguRCLUGuwP4vZmtN7MVg+GQEGJoqPXX+OXu3mVmUwCsM7P/dff3+n5D9kNgBQBMnz69xuGEEEWp6Z3d3buy/3sAvAFgcc73rHT3DnfvmDhxYi3DCSFqoHCwm9l4M2u68BrATwFsGSzHhBCDSy2/xk8F8Ea21X8ZgH939/9iHUaMGBFmbLGCfJGNHbtUVCJhUll0TTYWy8hiY7EMKiZDRddkMg6TDtlzYZJXBJOGGGweWUZZdG9F5EuAPzPWr8hxXmzuo+uxtV042N19J4AfFe0vhCgXSW9CJIKCXYhEULALkQgKdiESQcEuRCKUXnAykpSYfBXZmHzCMpcYTD45fPhwbvvu3bvDPsePHw9tLBNtypQpoY1ljkXz293dHfbZs2dPaIvuGeCZXGweIzZs2BDa9u3bF9ra29tD27Jly3Lb586dG/ZhBT2ZhMaksqhAJFBsriJUcFIIoWAXIhUU7EIkgoJdiERQsAuRCKXuxgPxjmVU3w0ATp48mdvOatAVTWZgfvT09OS2b9y4MezDdrpZDbrx48eHNnZv0ZFYBw8eDPsw5YIpDezeIqWBKQnbtm0LbVu2xAmVs2bNCm3Rs46OFAOAK6+8MrSx9cESitgcjx07Nred7eCzxKYIvbMLkQgKdiESQcEuRCIo2IVIBAW7EImgYBciEUqX3iIphCXCRDJJb29v2IfJIKx2GkuCiOSTvXv3hn3Wr18f2pj/7JpMeoukTZZ0E0k/AHDoULHDfqIkGVZhmM19dIwTwH2M5FIGk0TZkVeR7Alw/6OELia/FqlBp3d2IRJBwS5EIijYhUgEBbsQiaBgFyIRFOxCJEK/0puZrQLwMwA97n5j1tYK4FUAswDsBnCfux+txRGWwRbJSUwyYrXCGhsbBzwWENeFW758edhn0qRJoY1Jb7t27Qpt7GioSKZcsGBB2IfZmGTHJMBIAmIy2Zo1a0IbO/KKZdLdcsstue2sxh/LNit6/BPLlovWHFuLkR+11qD7DYA7Lml7DMA77j4HwDvZ10KIYUy/wZ6dt37pj9W7AazOXq8GcM/guiWEGGyKfmaf6u4XahPvR+VEVyHEMKbmDTqvfEgIPyiY2Qoz6zSzTva5SwgxtBQN9gNm1gYA2f/hHyC7+0p373D3DlYKSAgxtBQN9rUAHspePwTgzcFxRwgxVFQjvb0M4DYAk8zsKwBPAHgKwGtm9jCAPQDuq2awESNGhJk8TLaI5Boma7HsJCa9HTt2LLRFPrIChc3NzaGNZUIxWfHEiROhbcyYMbntTJKZPXt2aGMSJitiGd0bkxvZUVldXV2h7cYbbwxt0RphWZZFpTeW0ceO35o8eXJuO5urIkdG9Rvs7v5AYLp9wKMJIeqG/oJOiERQsAuRCAp2IRJBwS5EIijYhUiEUgtOuntYtLGIpMHkNZZlxApOsoJ9URFL1odJPKwfKzbIJJ7oXDxWePGLL74Ibe3t7YX8iM6IW7duXdhn8+bNoW3evHmhbdGiRaGNZbdFFJG1AL6umJQaFedksIy4CL2zC5EICnYhEkHBLkQiKNiFSAQFuxCJoGAXIhFKld7MDA0NDbk2JkNFUgiT15iUx2wsyyuyMamGXY/dM/Mxki+BWOpjtQSiTDmA31t0RhkQF4hkEiDL8poxY0ZoY5mFUSFT9lwYg11UEoifGfORSboRemcXIhEU7EIkgoJdiERQsAuRCAp2IRKh1N14IN5hZDu70c50kWQAgO9ms+SaIrumUfIMwO+ZHfEUKRpAPFesD1MFWD02lkDz7rvv5rZ/+OGHYR82V8uWLQttS5YsCW0TJkzIbS+q1rBnxtYO28Uvsr7ZWBF6ZxciERTsQiSCgl2IRFCwC5EICnYhEkHBLkQiVHP80yoAPwPQ4+43Zm1PAvg5gAvn/zzu7m9Xca1QAmIyQ3SEz+HDh+lYEUy2YLXfIomKJUAwyYvVHmPyDzv2Kromk/Ki47UAPh+nT58Obbt27cptZ3IdS2i59tprQ1tbW1toi+aRPTN2z2xdsaPDIgkQiI/zYglKETTxqor+vwFwR077s+4+P/vXb6ALIepLv8Hu7u8B0MHqQnzPqeUz+yNmtsnMVplZXFNYCDEsKBrszwO4GsB8AN0Ano6+0cxWmFmnmXWyz9hCiKGlULC7+wF3P+fu5wG8AGAx+d6V7t7h7h3sPHIhxNBSKNjNrO/2570AtgyOO0KIoaIa6e1lALcBmGRmXwF4AsBtZjYfgAPYDeAX1Q5YpPYXO1anCEzmYxlxUcYTuycmoTU1NYU29pGH2RobG3PbWQ23qF4cAGzcuDG0ffDBB6Ftx44due2zZs0K+yxdujS03XzzzaEtumcgPg6LPTMmKbLstaIya5FsyiJx1G+wu/sDOc0vDngkIURd0V/QCZEICnYhEkHBLkQiKNiFSAQFuxCJUHrByUiKYpJXJL2xQonMxiQSJstFcgfLoisqn7CMpyK2osULmWTU1dUV2qIsO5ahds0114S26BgngMubRYqSsj4se7DoeozGY88l6sPmQu/sQiSCgl2IRFCwC5EICnYhEkHBLkQiKNiFSIRSpbeGhgZaeC8iylxqaWkJ+zBZK8oyAvhZXpHcwWQVJpOxwoasIGJPT09oa21tzW0fN25c2Iexd+/e0Pbpp5+Gtt7e3tx2VjiSZcRNnBgXQ2L3FkllReVS9qzHjh0b2piEGT2zIplyVP4LLUKIHxQKdiESQcEuRCIo2IVIBAW7EIlQ6m78+fPncebMmVwbS4Rhu+cRLCGA7cSyGmNFdnbZfY0ePTq0MaWB7fBHigHbwd+5c2do2759e2hjysVVV12V237rrbeGfW666abQxnaZ9+3bF9qiHW22g8/WzlAQrR/mB1tXEXpnFyIRFOxCJIKCXYhEULALkQgKdiESQcEuRCJUc/zTTAC/BTAVleOeVrr7c2bWCuBVALNQOQLqPnc/yq7l7oVkBlZ/jI0VwaQyNlYkdzAJikkkLHGiqJwX2VhizbFjx0Lb/v37Qxs7liuq8zdt2rSwz+TJk0Mbkym//PLL0BY9m6LSW5Facv1dM3rWRY54YlTzzn4WwK/c/XoASwD80syuB/AYgHfcfQ6Ad7KvhRDDlH6D3d273f2T7PVJANsBzABwN4DV2betBnDPEPkohBgEBvSZ3cxmAVgA4CMAU929OzPtR+XXfCHEMKXqYDezRgBrADzq7if62rzygST3Q4mZrTCzTjPrZEcDCyGGlqqC3cxGohLoL7n761nzATNry+xtAHL/+NrdV7p7h7t3RBU5hBBDT7/BbpUtwRcBbHf3Z/qY1gJ4KHv9EIA3B989IcRgUU062TIADwLYbGYbsrbHATwF4DUzexjAHgD31eIIk5MiGYf1YTJIUfkkyr5jfVidOZbNx2quManvs88+y21///33wz6bNm0KbawO2nXXXRfaFi5cmNvO5DUmTzE/2BxH8hUbq2iNQrau2DWjdczk12gsel+h5f87/xFAJPjd3l9/IcTwQH9BJ0QiKNiFSAQFuxCJoGAXIhEU7EIkQqkFJ4FYpmLyVSS9sT7MxrKJivRjkgsrYMlkEpYRxyTH6K8Ut27dGvZhBSeZrDV//vzQdvvt+UINO/7p9OnToY0dQxWtDyB+ZkUk1v76DXbWG5PeIh9pxl5oEUL8oFCwC5EICnYhEkHBLkQiKNiFSAQFuxCJUKr05u6Fznpra2vLbT98+HDYh8kgTNJobm4ObUUyqJjtwIEDoY1lVzE5Lzr3jBWOZLLWokWLQtvSpUtD29Sp+YWLWPYau+cZM2aENlZMM1pvTC5lxS2LZuaxbL/oHD7mB1vfEXpnFyIRFOxCJIKCXYhEULALkQgKdiESofREmCKJCadOncptZ7vqRY/wYckY0S44S8RgCS3RLizAd+p7e3tD265du3Lb2RFPrN7dDTfcENrmzp0b2iZMmJDbfvDgwbAPey7sWC620x3dN9s5Z2uH7fwX2SFnsPVdRBnSO7sQiaBgFyIRFOxCJIKCXYhEULALkQgKdiESoV/pzcxmAvgtKkcyO4CV7v6cmT0J4OcALmgpj7v72/1cC2PGjMm1MZkhkqjGjRsX9ila3+3EiROhbbBr0EXyFAB0d3eHtm3btoW2o0eP5ra3tLSEfebNmxfa5syZE9pYLb9p06bltjc2NoZ9WCJMlOADAE1NTaGtCMwPJs1GaxsodnwV6xOtK/ZMqtHZzwL4lbt/YmZNANab2brM9qy7/1MV1xBC1JlqznrrBtCdvT5pZtsBxPmGQohhyYA+s5vZLAALAHyUNT1iZpvMbJWZxX+GJYSoO1UHu5k1AlgD4FF3PwHgeQBXA5iPyjv/00G/FWbWaWadUU1zIcTQU1Wwm9lIVAL9JXd/HQDc/YC7n3P38wBeALA4r6+7r3T3DnfvaG1tHSy/hRADpN9gt8r23osAtrv7M33a+9aKuhfAlsF3TwgxWFSzG78MwIMANpvZhqztcQAPmNl8VOS43QB+Uc2AUWYQk6iiPizbrGhdOCajRXXyWLYT++hy6NCh0LZ+/frQtnHjxtAW+cjktQULFoS29vb20MZknui+WdYYlY3IkUwnT54MbdGzLvKcAZ59x9YjI7pvNhbzP6Ka3fg/AsjzhmrqQojhhf6CTohEULALkQgKdiESQcEuRCIo2IVIhGFTcJJJCdERPix7rYgPAJddilyPSU1MAmRHMm3fvj20XXHFFbntrDjk8uXLQxvLRmSSV5QdFmV4ATyLMSo6ChR7Zgx2zwy2hlkmXbSOmR+RzMfWot7ZhUgEBbsQiaBgFyIRFOxCJIKCXYhEULALkQjDRnorctYbKzTIrseyiSKZD4gLCrKMPVaEkGVJMamGZctF0hYrpMlgchiTFaM5YffFxmJZb6yYZiQBsmKOTMpjchjzkV0zkmB11psQohAKdiESQcEuRCIo2IVIBAW7EImgYBciEUqV3s6dO4fjx4/n2pgMNWNG/gE0TCZjMBmKZdJFcgfLyDp8+HBo27FjR2jbvHlzaGOyUTQnLNuMXY/JSdOnTw9tBw8ezG1nch3LGpsyZUpoY3McSbBsLHbPo0ePDm1sHpubm0NbFBNsLRbJzNM7uxCJoGAXIhEU7EIkgoJdiERQsAuRCP3uxpvZGADvARidff/v3P0JM5sN4BUAlwNYD+BBd4+zHCrXokkjEdFuN9tRZbW4mA9shz/apWVKAkuEYTaWyMN2fSPYbnxbW1toY0lDvb29A7ax+WWJTSyRhO1aR+MVrQ3I/Ch6HFm0Hlkftr4jqnlnPwPgx+7+I1SOZ77DzJYA+DWAZ939GgBHATw84NGFEKXRb7B7hQs/pkdm/xzAjwH8LmtfDeCeoXBQCDE4VHs+e0N2gmsPgHUAvgBwzN0vKPtfAcj/yxchxLCgqmB393PuPh9AO4DFAOIi5JdgZivMrNPMOtnxxUKIoWVAu/HufgzAHwAsBdBiZhd2rNoBdAV9Vrp7h7t3tLa21uKrEKIG+g12M5tsZi3Z67EAfgJgOypB/5fZtz0E4M0h8lEIMQhUkwjTBmC1mTWg8sPhNXf/TzPbBuAVM/sHAJ8CeLGaASNJickWUY0xloDCpCuWRMDqoBWp+8Vqrk2YMCG0zZs3L7RNnTo1tM2cOTO3fdq0aWEfNo/s3thcRZIdmw8m5TFZjl0zkmeZ/MpkraJ15li/aK0yH9n6Dn3o7xvcfROABTntO1H5/C6E+B6gv6ATIhEU7EIkgoJdiERQsAuRCAp2IRLBmLQy6IOZHQSwJ/tyEoD4HKPykB8XIz8u5vvmx5XuPjnPUGqwXzSwWae7d9RlcPkhPxL0Q7/GC5EICnYhEqGewb6yjmP3RX5cjPy4mB+MH3X7zC6EKBf9Gi9EItQl2M3sDjP7zMw+N7PH6uFD5sduM9tsZhvMrLPEcVeZWY+ZbenT1mpm68zsz9n/E+vkx5Nm1pXNyQYzu6sEP2aa2R/MbJuZbTWzv87aS50T4kepc2JmY8zsT2a2MfPj77P22Wb2URY3r5pZXA00D3cv9R+ABlTKWl0FYBSAjQCuL9uPzJfdACbVYdxbASwEsKVP2z8CeCx7/RiAX9fJjycB/E3J89EGYGH2ugnADgDXlz0nxI9S5wSAAWjMXo8E8BGAJQBeA3B/1v4vAP5qINetxzv7YgCfu/tOr5SefgXA3XXwo264+3sALq3RdTcqhTuBkgp4Bn6Ujrt3u/sn2euTqBRHmYGS54T4USpeYdCLvNYj2GcA2Nvn63oWq3QAvzez9Wa2ok4+XGCqu3dnr/cDiCtUDD2PmNmm7Nf8If840Rczm4VK/YSPUMc5ucQPoOQ5GYoir6lv0C1394UA7gTwSzO7td4OAZWf7Kj8IKoHzwO4GpUzAroBPF3WwGbWCGANgEfd/aJztcuckxw/Sp8Tr6HIa0Q9gr0LQN/aSWGxyqHG3buy/3sAvIH6Vt45YGZtAJD931MPJ9z9QLbQzgN4ASXNiZmNRCXAXnL317Pm0uckz496zUk29jEMsMhrRD2C/WMAc7KdxVEA7gewtmwnzGy8mTVdeA3gpwC28F5DylpUCncCdSzgeSG4Mu5FCXNilaJvLwLY7u7P9DGVOieRH2XPyZAVeS1rh/GS3ca7UNnp/ALA39bJh6tQUQI2Athaph8AXkbl18HvUPns9TAqZ+a9A+DPAP4HQGud/Pg3AJsBbEIl2NpK8GM5Kr+ibwKwIft3V9lzQvwodU4A3IRKEddNqPxg+bs+a/ZPAD4H8B8ARg/kuvoLOiESIfUNOiGSQcEuRCIo2IVIBAW7EImgYBciERTsQiSCgl2IRFCwC5EI/wfi5DoyocN8sQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Record the optimization process\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "    \n",
    "net = NET(Mat_generator(), weight)\n",
    "\n",
    "# We use Adam optimizer for better performance\n",
    "# One can change it to SGD or RMSprop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# Optimization loop\n",
    "for itr in range(ITR):\n",
    "\n",
    "    # Forward propagation to calculate loss function\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # Back propagation minimizes the loss function\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # Record optimization intermediate results\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",
    "# Record the last two unitary matrices learned\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",
    "## References\n",
    "\n",
    "[1] Wang, X., Song, Z., & Wang, Y. Variational Quantum Singular Value Decomposition. [Quantum, 5, 483 (2021).](https://quantum-journal.org/papers/q-2021-06-29-483/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}