tutorial/machine_learning/VQSVD_EN.ipynb

{
 "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": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:44:34.008567Z",
     "start_time": "2021-03-09T03:44:29.796997Z"
    }
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "from numpy import pi as PI\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import unitary_group\n",
    "from scipy.linalg import norm\n",
    "\n",
    "import paddle\n",
    "from paddle import matmul, transpose, trace\n",
    "from paddle_quantum.circuit import *\n",
    "from paddle_quantum.utils import *\n",
    "\n",
    "\n",
    "\n",
    "# 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": 2,
   "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": 3,
   "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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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\n",
    "theta_size = N * block_len * cir_depth  # size of the network parameter theta\n",
    "\n",
    "\n",
    "# Define quantum neural network\n",
    "def U_theta(theta):\n",
    "\n",
    "    # Initialize the network with UAnsatz\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # Build QNN\n",
    "    for layer_num in range(cir_depth):\n",
    "        \n",
    "        for which_qubit in range(N):\n",
    "            cir.ry(theta[block_len * layer_num * N + which_qubit], which_qubit)\n",
    "        \n",
    "        for which_qubit in range(N):\n",
    "            cir.rz(theta[(block_len * layer_num + 1) * N + which_qubit], which_qubit)\n",
    "\n",
    "        for which_qubit in range(1, N):\n",
    "            cir.cnot([which_qubit - 1, which_qubit])\n",
    "        cir.cnot([N - 1, 0])\n",
    "\n",
    "    return cir.U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we complete the main part of the algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:46:12.944634Z",
     "start_time": "2021-03-09T03:44:50.626213Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -235.9752\n",
      "iter: 10 loss: -1635.7805\n",
      "iter: 20 loss: -1967.8480\n",
      "iter: 30 loss: -2154.2274\n",
      "iter: 40 loss: -2257.8776\n",
      "iter: 50 loss: -2311.8162\n",
      "iter: 60 loss: -2342.3099\n",
      "iter: 70 loss: -2358.9515\n",
      "iter: 80 loss: -2369.3172\n",
      "iter: 90 loss: -2375.9682\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnHklEQVR4nO3deXxU9b3/8ddnZrInJEACgYAssiigYkUFt1q1BVuvdLPaW61ttba9tZvdsPXe295b78+23u6LtbW2dtF6tYttrbjUBW2tLCqySEFABAmELQFC9s/vj3OCQ0wwZGZykpn38/GYx8z5njNzPicH8s75fuecY+6OiIhIKmJRFyAiIoOfwkRERFKmMBERkZQpTEREJGUKExERSZnCREREUqYwEZEjYmb7zGxi1HXIwKIwkUHHzDaa2XkRrfs0M/urme01s3oz+6OZTevH9R/cdjN7n5k9nuH1PWJmVya3uXupu6/P5Hpl8FGYiPSSmc0B7gf+AIwGJgDPAk+k+y91C2T0/6eZJTL5+ZJbFCaSNcyswMy+ZWYvh49vmVlBOK/SzP5kZnvMbJeZLer8ZW1mnzezLeHRxhozO7eHVXwNuM3dv+3ue919l7tfBzwJfCn8rNVmdkFSTQkzqzOz14XTs83sb2Edz5rZ2UnLPmJm15vZE0Aj0GNAmdmxwE3AnLDbaU/Sz+BGM9tkZtvM7CYzKwrnnW1mm8PtrQVuNbOh4c+lzsx2h6/HhMtfD5wJfC9cx/fCdjezSeHrcjO7LXz/i2Z2XdLP9X1m9nhYz24z22Bm5/d6h8qgojCRbPJFYDYwEzgBOAW4Lpz3aWAzUAWMBL4AuJlNBa4GTnb3MmAusLHrB5tZMXAa8H/drPdO4I3h69uBdyfNmwvscPdlZlYD/Bn4CjAM+Axwt5lVJS1/GXAVUAa82NOGuvtq4MPA38Nup4pw1g3AlPBnMAmoAf4j6a3V4brHheuJAbeG00cBB4Dvhev4IrAIuDpcx9XdlPJdoJwg+F4PvBd4f9L8U4E1QCVBGN9iZtbTdsngpTCRbPIe4L/cfbu71wFfJvjlDNAKjALGuXuruy/y4MJ07UABMM3M8tx9o7u/0M1nDyP4/7K1m3lbCX5ZAvwauDAMH4B/JQgYgEuBe939XnfvcPcHgCXAm5M+62fuvtLd29y99Ug2PvwlfRXwqfCoaS/wP8AlSYt1AP/p7s3ufsDdd7r73e7eGC5/PUEo9GZ98fCzrw2P1DYC/8srP3OAF939x+7eDvycYB+MPJLtksFBYSLZZDSH/jX/YtgG8HVgHXC/ma03swUA7r4O+CRBN9V2M7vDzEbzarsJfhGP6mbeKGBH0uetBv4lDJQLCQIGgr/+Lwq7uPaEXVNndPnMl45kg7uoAoqBpUmff1/Y3qnO3Zs6J8ys2Mx+FHZRNQCPARVhULyWSiCPV//Ma5KmaztfuHtj+LL0CLZJBgmFiWSTlwl+YXc6Kmwj/Mv50+4+keAX/DWdYyPu/mt3PyN8rwNf7frB7r4f+DtwUTfrfRfwUNJ0Z1fXfGBVGDAQBMUv3L0i6VHi7jckr+oItrfrsjsIuqmmJ31+ubuXHuY9nwamAqe6+xDgrLDdeli+6/paefXPfMsRbINkCYWJDFZ5ZlaY9EgQ/BK/zsyqzKySYKzglwBmdoGZTQq7guoJurc6zGyqmZ0TDtQ3Efwy7uhhnQuAy83s42ZWFg5efwWYQ9Cl1ukO4E3AR3jlqISwln8xs7lmFg/rPrtzwLsPtgFjzCwfwN07gB8D3zSzEeF215jZ3MN8RhnBNu8xs2HAf3azjm6/CBB2Xd0JXB/+PMYB14TbKTlGYSKD1b0EvwQ7H18iGNheAiwHngOWhW0Ak4EHgX0ERxg/cPeHCcZLbiD4K7sWGAFc290K3f1xggH1txOMk7wInAic4e5rk5bbGq7jNOA3Se0vERytfAGoIzhS+Sx9/3/4V2AlUGtmO8K2zxN05z0Zdls9SHDk0ZNvAUUE2/8kQbdYsm8D7wy/jfWdbt7/MWA/sB54nCA8f9qnrZFBzXRzLBERSZWOTEREJGUKExERSZnCREREUqYwERGRlOXshd4qKyt9/PjxUZchIjKoLF26dIe7V3Vtz9kwGT9+PEuWLIm6DBGRQcXMur1mnLq5REQkZQoTERFJmcJERERSpjAREZGUKUxERCRlChMREUmZwkRERFKmMDlCv396C798ssdbc4uI5CSFyRG697mt3Pb3jVGXISIyoChMjtDoiiK27D6A7gMjIvIKhckRGjO0iP0t7TQ0tUVdiojIgKEwOUKjK4oA2LL7QMSViIgMHAqTI9QZJi/vUZiIiHRSmByhms4wqVeYiIh0UpgcoeEl+eQnYurmEhFJojA5QrGYMbq8kC3q5hIROUhh0gc1Q4sUJiIiSRQmfTC6vEgD8CIiSRQmfTC6oojte5tpaeuIuhQRkQFBYdIHNRVFuENtfVPUpYiIDAgKkz6oGRqeuKiuLhERQGHSJzpxUUTkUAqTPhhVXgjoyEREpJPCpA8K8+JUlhboyEREJKQw6aOaCp24KCLSSWHSR6MrdOKiiEgnhUkf1VQEJy7qJlkiIgqTPhtdUURTawe7G1ujLkVEJHIKkz7STbJERF4x4MLEzL5kZlvM7Jnw8eakedea2TozW2Nmc5Pa54Vt68xsQX/UOUYnLoqIHJSIuoAefNPdb0xuMLNpwCXAdGA08KCZTQlnfx94I7AZWGxm97j7qkwWqBMXRUReMVDDpDvzgTvcvRnYYGbrgFPCeevcfT2Amd0RLpvRMBlanEdhXkxHJiIiDMBurtDVZrbczH5qZkPDthrgpaRlNodtPbVnlJkxukKXohcRgYjCxMweNLMV3TzmAz8EjgZmAluB/03jeq8ysyVmtqSuri7lz6tRmIiIABF1c7n7eb1Zzsx+DPwpnNwCjE2aPSZs4zDtXdd7M3AzwKxZs1I+QaSmoojVWxtS/RgRkUFvwHVzmdmopMm3ASvC1/cAl5hZgZlNACYDTwGLgclmNsHM8gkG6e/pj1pHVxSxY18LTa3t/bE6EZEBayAOwH/NzGYCDmwEPgTg7ivN7E6CgfU24KPu3g5gZlcDC4E48FN3X9kfhXZ+o2trfRMTKkv6Y5UiIgPSgAsTd7/sMPOuB67vpv1e4N5M1tWd0RXBpehf3nNAYSIiOW3AdXMNJjUVOnFRRAQUJimpLn/lyEREJJcpTFJQkIhTVaabZImIKExSFJxr0hR1GSIikVKYpEgnLoqIKExSNjq8fa9ukiUiuUxhkqLRFUU0t3Wwa39L1KWIiERGYZKiVy5Fr3ETEcldCpMU6VwTERGFScp0kywREYVJyjpvkqUwEZFcpjBJ0cGbZNUrTEQkdylM0qCmoogtGoAXkRymMEmD0eU6cVFEcpvCJA1GVxRRt7eZ5jbdJEtEcpPCJA0672tSW6+uLhHJTQqTNNC5JiKS6xQmaaCz4EUk1ylM0kA3yRKRXKcwSYPCvDiVpbpJlojkLoVJmtSEl6IXEclFCpM0qRmqc01EJHcpTNKkpqKIzbsP0N6hm2SJSO5RmKTJ5JFlNLd1sGlXY9SliIj0O4VJmhxTXQbAmtqGiCsREel/CpM0mTyiDDNYU7sv6lJERPqdwiRNivLjjBtWzJptOjIRkdyjMEmjKSPLeL52b9RliIj0O4VJGh1TXcbGHftpatXVg0UktyhM0mhq9RA6HNZt17iJiOQWhUkaTa0uBeCf29TVJSK5RWGSRuOHl5Afj7FG4yYikmMUJmmUiMc4ekSpBuFFJOcoTNLsmOoydXOJSM5RmKTZ1OoyttY3Ud/YGnUpIiL9RmGSZlNHhpdV0dGJiOQQhUmaTa1WmIhI7lGYpNmo8kLKChO64KOI5JRIwsTMLjKzlWbWYWazusy71szWmdkaM5ub1D4vbFtnZguS2ieY2T/C9t+YWX5/bktXZsbUkWX8Uxd8FJEcEtWRyQrg7cBjyY1mNg24BJgOzAN+YGZxM4sD3wfOB6YB7w6XBfgq8E13nwTsBq7on03o2dTqMp6vbcBdN8oSkdwQSZi4+2p3X9PNrPnAHe7e7O4bgHXAKeFjnbuvd/cW4A5gvpkZcA5wV/j+nwNvzfgGvIaZYytoaGpjxRZ1dYlIbhhoYyY1wEtJ05vDtp7ahwN73L2tS3u3zOwqM1tiZkvq6urSWniyc48dScxg4crajK1DRGQgyViYmNmDZraim8f8TK3ztbj7ze4+y91nVVVVZWw9w0ryOWXCMIWJiOSMRKY+2N3P68PbtgBjk6bHhG300L4TqDCzRHh0krx8pOZOr+bLf1zF+rp9TKwqjbocEZGMGmjdXPcAl5hZgZlNACYDTwGLgcnhN7fyCQbp7/FghPth4J3h+y8H/hBB3a/ypunVACxcuS3iSkREMi+qrwa/zcw2A3OAP5vZQgB3XwncCawC7gM+6u7t4VHH1cBCYDVwZ7gswOeBa8xsHcEYyi39uzXdq6ko4riacnV1iUhOyFg31+G4+++A3/Uw73rg+m7a7wXu7aZ9PcG3vQacudNHcuP9/6S2vonq8sKoyxERyZiB1s2VVeaGXV0PrNLRiYhkN4VJBk0aUcrEyhKNm4hI1lOYZJCZ8abp1Ty5fid7m3RJehHJXgqTDJtz9HDaOpznNtdHXYqISMYoTDLs+JpyAJ5VmIhIFlOYZNjQknyOGlbM8s17oi5FRCRjFCb94Lgx5SzXkYmIZDGFST84YUw5W/YcYOe+5qhLERHJCIVJPzh+TAWAjk5EJGspTPrBjJpyzBQmIpK9FCb9oLQgwdFVpRqEF5GspTDpJ8ePKefZzfW6la+IZCWFST85vqacHfuaqW1oiroUEZG0U5j0k+PHVgDw7EsaNxGR7KMw6SfTRg0hETONm4hIVlKY9JPCvDhTRpbx3BYdmYhI9ulVmJjZJ8xsiAVuMbNlZvamTBeXbU4YG5wJr0F4Eck2vT0y+YC7NwBvAoYClwE3ZKyqLHX8mArqD7SyamtD1KWIiKRVb8PEwuc3A78I779uh1leujFvejWlBQm+//C6qEsREUmr3obJUjO7nyBMFppZGdCRubKy09CSfD5w+njufa6WVS/r6EREskdvw+QKYAFwsrs3AnnA+zNWVRa74oyJlBUm+OaD/4y6FBGRtOltmMwB1rj7HjO7FLgO0NeS+qC8OI8PnjmRB1Zt090XRSRr9DZMfgg0mtkJwKeBF4DbMlZVlnv/6eOpKM7jGw+siboUEZG06G2YtHnwfdb5wPfc/ftAWebKym5lhXlcddZEHl5Tp6MTEckKvQ2TvWZ2LcFXgv9sZjGCcRPpo0tnj6MwL8btizdFXYqISMp6GyYXA80E55vUAmOAr2esqhwwpDCPNx83inueeZnGlraoyxERSUmvwiQMkF8B5WZ2AdDk7hozSdHFs8ayr7mNe5+rjboUEZGU9PZyKu8CngIuAt4F/MPM3pnJwnLBKROGMaGyhDsXvxR1KSIiKeltN9cXCc4xudzd3wucAvx75srKDWbGu2aN5amNu3ihbl/U5YiI9FlvwyTm7tuTpncewXvlMN5xUg3xmHHnEh2diMjg1dtAuM/MFprZ+8zsfcCfgXszV1buGFFWyBumjuDupVtobdcVakRkcOrtAPxngZuB48PHze7++UwWlksuPnksO/Y18+iauqhLERHpk0RvF3T3u4G7M1hLznr9lCrKChPct7KW86aNjLocEZEjdtgwMbO9QHd3cjLA3X1IRqrKMfmJGOceM4IHV2+jrb2DRFzDUSIyuBz2t5a7l7n7kG4eZQqS9Jo3o5o9ja08tWFX1KWIiBwx/Qk8QJw1pYqCRIz7VuoERhEZfBQmA0RxfoLXT6ni/pXb6OjQPeJFZHCJJEzM7CIzW2lmHWY2K6l9vJkdMLNnwsdNSfNOMrPnzGydmX3HzCxsH2ZmD5jZ2vB5aBTblA5zp1dT29DEs5v3RF2KiMgRierIZAXwduCxbua94O4zw8eHk9p/CHwQmBw+5oXtC4CH3H0y8FA4PSide+wI4jFj4cptUZciInJEIgkTd1/t7r2+M5SZjQKGuPuT4X1VbgPeGs6eD/w8fP3zpPZBp6I4nzkTh7NwZS3BZoqIDA4Dccxkgpk9bWaPmtmZYVsNsDlpmc1hG8BId98avq4FejxRw8yuMrMlZrakrm5gniA4d/pINuzYz9rtulaXiAweGQsTM3vQzFZ085h/mLdtBY5y9xOBa4Bfm1mvv4IcHrX0+Ce9u9/s7rPcfVZVVVWvt6U/zZ1eTX48xv/cu1oD8SIyaGQsTNz9PHef0c3jD4d5T7O77wxfLyW41/wUYAvBDbk6jQnbALaF3WCd3WHJF6QcdEYMKeTfLziWR9bU8aPH1kddjohIrwyobi4zqzKzePh6IsFA+/qwG6vBzGaH3+J6L9AZSvcAl4evL09qH7QunT2Otxw3ihvvX8OSjTqJUUQGvqi+Gvw2M9sMzCG4p/zCcNZZwHIzewa4C/iwu3f+Nv034CfAOoIjlr+E7TcAbzSztcB54fSgZmb8v3ccx5ihRXzs9qfZtb8l6pJERA7LcvVbQ7NmzfIlS5ZEXcZhrdhSz9t+8ATvOXUcX7pwetTliIhgZkvdfVbX9gHVzSWHmlFTzr+cMJo7l7xEfWNr1OWIiPRIYTLAXXnGRBpb2vn1U5uiLkVEpEcKkwFu2ughnD5pOD/72wZa2nQnRhEZmBQmg8CVZ0xkW0Mz9z639bUXFhGJgMJkEHj9lCqOrirhx4vW6zIrIjIgKUwGgVjMuOKMiax8uYG/r98ZdTkiIq+iMBkk3v66GoaX5PP9h9dFXYqIyKsoTAaJwrw4Hzn7aJ5Yt5O/vbAj6nJERA6hMBlELp09juohhdy4cI3GTkRkQFGYDCKFeXGuPmcSyzbt4ZE1A/MS+iKSmxQmg8y7Zo1l7LAibrx/jS5RLyIDhsJkkMlPxPjkuVNY+XIDf1z+ctTliIgACpNB6a0n1jCjZgjX3Pkstzy+QeMnIhI5hckgFI8Zv/7gbM45ZgT//adVfOKOZ2hsaYu6LBHJYQqTQWpIYR4/uvQkPjt3Kn9c/jIf+eUyHaGISGQUJoNYLGZ89A2TuO4t03j0n3X8ZUVt1CWJSI5SmGSBy+eM49hRQ/ivP65if7O6u0Sk/ylMskAiHuMrb51ObUMT3/nr2qjLEZEcpDDJEieNG8ZFJ43hlkUbWLttb9TliEiOUZhkkQXnH0NJQYLP3LWcAy3tUZcjIjlEYZJFhpcW8NV3HM/yzXv42O1P09auOzOKSP9QmGSZeTOq+fKF03lw9Tb+/Q8r9XVhEekXiagLkPR775zx1NY38YNHXqC8KI/Pzp1KPGZRlyUiWUxhkqU+O3cquxtbuOnRF1i8cRc3XnQCEypLoi5LRLKUurmylJnxP287jm9dPJO12/Zy/rcf4/anNkVdlohkKYVJFjMz3npiDfd/6vWcPH4Y1/72OR5YtS3qskQkCylMckB1eSE/fu+s8ErDz7BpZ2PUJYlIllGY5IjCvDg/fM9JxMz48C+X0tSq81BEJH0UJjlk7LBivnXxTFZtbWDB3csVKCKSNgqTHPOGY0ZwzRun8PtnXmbutx7j4TXboy5JRLKAwiQHffzcyfziilOIm/H+WxfzwduWaBxFRFKiMMlRZ06u4i+fPJPPzZvKE+t2cN43HuWGvzzPPl3CXkT6QGGSwwoScf7t7Ek8/JmzueCEUdz06Aucc+MjPLe5PurSRGSQUZgII4cU8o13zeT3Hz2dvHiMd//4Sf7+ws6oyxKRQURhIgfNHFvB3R85jVHlhVx+61Pcv1K3ARaR3lGYyCGqywu580NzOHbUED70y6V84GeLeXDVNl3OXkQOS2EirzK0JJ9fX3kqV79hEiu21HPlbUs482sP86AuxSIiPYgkTMzs62b2vJktN7PfmVlF0rxrzWydma0xs7lJ7fPCtnVmtiCpfYKZ/SNs/42Z5ffz5mSlkoIEn37TVJ5YcA4/uuwkhhbnc+VtS/jvP62ipU1HKSJyqKiOTB4AZrj78cA/gWsBzGwacAkwHZgH/MDM4mYWB74PnA9MA94dLgvwVeCb7j4J2A1c0a9bkuXy4jHmTq/mdx89jfedNp5bHt/AO2/6m85LEZFDRBIm7n6/u3ee0PAkMCZ8PR+4w92b3X0DsA44JXysc/f17t4C3AHMNzMDzgHuCt//c+Ct/bQZOaUgEedLF07npktPYuOO/bzlu4s0QC8iBw2EMZMPAH8JX9cALyXN2xy29dQ+HNiTFEyd7d0ys6vMbImZLamrq0tT+bll3oxq/vzxM5lQWcJVv1jKV/60iuY2XeNLJNdlLEzM7EEzW9HNY37SMl8E2oBfZaqOZO5+s7vPcvdZVVVV/bHKrDR2WDH/9+E5XD5nHD95fAOn3/BXvnH/GrY3NEVdmohEJGO37XX38w4338zeB1wAnOvuHjZvAcYmLTYmbKOH9p1AhZklwqOT5OUlgwoScb48fwZzZ1Rzy6INfPfhdfzgkRf44FkTueaNU8iLD4SDXhHpL5HcA97M5gGfA17v7skjufcAvzazbwCjgcnAU4ABk81sAkFYXAL8q7u7mT0MvJNgHOVy4A/9tyVy2tGVnHZ0JS/u3M93/7qOHz7yAk9t2MV3330ioyuKoi5PRPpJVH8+fg8oAx4ws2fM7CYAd18J3AmsAu4DPuru7eFRx9XAQmA1cGe4LMDngWvMbB3BGMot/bspAjBueAk3XnQC33n3iTy/tYE3f2cRv1m8SfdMEckR9koPU26ZNWuWL1myJOoystKGHfv55B1P8+zmeipL83nPqeO4bM44KksLoi5NRFJkZkvdfdar2hUmkgnuzt9e2MlPH9/AQ89vJz8R4x2vG8OVZ07g6KrSqMsTkT7qKUwiGTOR7GdmnD6pktMnVfJC3T5ueXwDdy3dzO1PbeKC40dx3VumUV1eGHWZIpImOjKRfrNjXzO3PrGBnyzaQCJmfOqNU3jfaeNJ6JtfIoOGurm6UJhEZ9PORv7jnhU8sqaOo4YV886TxvD219UwZmhx1KWJyGtQmHShMImWu3P/qm3c+sQGnly/C4Czp1Zx3VuOZdKIsoirE5GeKEy6UJgMHC/tauS3y7Zwy+PraWxp54ozJ/DxcyZTUqAhPZGBRmHShcJk4Nmxr5mv3fc8dy7ZTGVpPpfNHs97Zh+lrxSLDCAKky4UJgPXsk27+e5Da3l4TR35iRgXnjCay2aP44SxFVGXJpLzFCZdKEwGvhfq9nHrExv47bItNLa0c1xNOZfNHseFM0dTmBePujyRnKQw6UJhMnjsbWrl909v4ZdPbmLNtr0Hu8AumzOOYSW6saZIf1KYdKEwGXzcnb+v38lPFm3gr89vJz8e46RxQzljciVnTKrk+DHlBPdLE5FMUZh0oTAZ3NZu28v/Ld3MorU7WL21AYCjq0q45OSjeMdJY3TEIpIhCpMuFCbZo25vMw8/v507Fm9i2aY95MWNmWMrOHXCcE6dOIyTxw/TGItImihMulCYZKc1tXv53dNb+Pv6nazYUk97h1OcH+esyVWcN20kr59SRVWZvmos0le60KPkhKnVZSw4/xgA9jW3sXjjLh5avY0HV23nvpW1ABw7aghnTa5k9tHDmTVuKGWFeVGWLJIVdGQiOcHdWbGlgcfW1rFobR1LX9xNa7sTM5hRU87pkyo579gRzBw7lHhMg/giPVE3VxcKk9zW2NLGshf38NSGnTy5fhdLN+2mvcMZXpLPnKOHc+JRQ5k5toIZNUMoSGi8RaSTurlEkhTnJ4KvFE+uBKC+sZVH19bx0OptLN6wiz8t3wpAYV6MUycM56wpVZw6YRiTRpRqMF+kGzoyEenG9oYmlm3aw5Prd/LY2jrW1+0HIGZw1LBiplaXcfyYCo6rKWfa6CEML8nXOS6SE3RkInIERgwpZN6MaubNqAaCKxs/u3kPa7ftY+32vax6uYGFK7cdXL60IMGYoUWMG17MpBGlTB5RxqQRpRw1vJghGuCXHKAwEemFscOKGTvs0Jt31R9oZeWWelbX7uWlXY28tKuRtdv38eDq7bR3vHLEP6wkn3HDi5lQWcLEyhImVpUybngx44aXUKrL7EuW0L9kkT4qL8rjtEmVnDap8pD2lrYOXty5n3Xb9/HirkZe3NnIxh37+du6nfx22ZZDlq0szWfM0GLGDC1Keg4eo8qLdE8XGTT0L1UkzfITMSaPLGPyyFffMXJ/cxsbd+7nxZ2N4WM/m3cfYMWWehaurKW1/dAxzLLCBKPKCxlVXsToiiJqKgqpLi9i5JACRg4pZGRZIUOKEhqvkcgpTET6UUlBgumjy5k+uvxV89o7nLq9zWze3cjm3QfYWt9EbX3wvLW+iRVb6tm5v+VV7yvKi1NdXkhVWQFVpQUML82nsrSAytICqsoKqCzNZ1hJPkNL8ikrUPBIZihMRAaIeMyoLi+kuryQWeO7X6aptZ1tDU1sa2imtqGJ7Q1N1NY3sbWhibq9zayubWDnvhbqD7R2+/5EzBhWEoTL8NJ8KoryKS/OY2hx3sHX5UV5DCnMo6wwcfC5rDBBIh7L3MbLoKcwERlECvPijBtewrjhJYddrqWtg537m6nb28zOfS3s2t/C7sYWdu5vYde+8Hl/M7X1DexpbGXPgdZDvjTQnaK8+MFgKS3MY0hhgtKC4FFSELaHr0sK4hTnJyjJT1BcEA+e8+MU5sUpzo9TlBcnpisNZBWFiUgWyk/EGFUeDOL3hruzr7mNPY2t1B9opaGplYYDbextamVvU1v4aGVfc/C6oamV/c1tbGtoYm9TG/ua29jf3MZr5NEhChIxivOD0CnMi1EUhkxh3ivPhXmx8DlOYSJGQV6cgkTQVpA0HTziFOTFDk7nx+PkJ2Lkd04nYiRipm6+DFGYiAhmRllhHmWFeYzt42e4O40t7eEjCJjGlnb2N7exvzloa2pt50BrsMyBcNn9LW00t3YcnLevuY26vc0caG0P2tvaaWptp6m1I+XtjBnkxYNgye98Dl8nt+cljLzOtniMvLiRCJ/z4jESsaTX8c5ljXjYnoh1todtMSMeMxI9TCc6pw8+x4jHjbi90h5Lmh+PBfMG0tGdwkRE0sLMwi6uBJD+y/y7O81tHTS3dtDcHgRNc1sQMkF7O83tHbS0BdMtBx/ttLSH72vroLU9nN/eQWv43BK2t7Q7rW0dNLV2sLep7WB7W4eHyzptHR20tTst7R20tXcc0dFYuplxMFQ6gydmkIjHiJkRjyXNTwqgn15+MkcNL37tFRwBhYmIDApmdrDLCwbOVQXaO/xg4LS3O61h2Bxs6wie29r94HRru9PREUx3hlOHe9Dunct20N7Bwfe3d4TzwvW0e9DW3hG87uhw2jugI2xv6wjb/JX3dj4X5KX/yxQKExGRFATdTrr4p77rJyIiKVOYiIhIyhQmIiKSMoWJiIikTGEiIiIpU5iIiEjKFCYiIpIyhYmIiKTM3CO8FkCEzKwOePEI3lIJ7MhQOQNVLm4z5OZ25+I2Q25ud6rbPM7dq7o25myYHCkzW+Lus6Kuoz/l4jZDbm53Lm4z5OZ2Z2qb1c0lIiIpU5iIiEjKFCa9d3PUBUQgF7cZcnO7c3GbITe3OyPbrDETERFJmY5MREQkZQoTERFJmcLkNZjZPDNbY2brzGxB1PVkipmNNbOHzWyVma00s0+E7cPM7AEzWxs+D4261nQzs7iZPW1mfwqnJ5jZP8J9/hszy4+6xnQzswozu8vMnjez1WY2J9v3tZl9Kvy3vcLMbjezwmzc12b2UzPbbmYrktq63bcW+E64/cvN7HV9Xa/C5DDMLA58HzgfmAa828ymRVtVxrQBn3b3acBs4KPhti4AHnL3ycBD4XS2+QSwOmn6q8A33X0SsBu4IpKqMuvbwH3ufgxwAsH2Z+2+NrMa4OPALHefAcSBS8jOff0zYF6Xtp727fnA5PBxFfDDvq5UYXJ4pwDr3H29u7cAdwDzI64pI9x9q7svC1/vJfjlUkOwvT8PF/s58NZICswQMxsDvAX4SThtwDnAXeEi2bjN5cBZwC0A7t7i7nvI8n1NcJvyIjNLAMXAVrJwX7v7Y8CuLs097dv5wG0eeBKoMLNRfVmvwuTwaoCXkqY3h21ZzczGAycC/wBGuvvWcFYtMDKqujLkW8DngI5wejiwx93bwuls3OcTgDrg1rB77ydmVkIW72t33wLcCGwiCJF6YCnZv6879bRv0/Y7TmEihzCzUuBu4JPu3pA8z4PvkWfNd8nN7AJgu7svjbqWfpYAXgf80N1PBPbTpUsrC/f1UIK/wicAo4ESXt0VlBMytW8VJoe3BRibND0mbMtKZpZHECS/cvffhs3bOg97w+ftUdWXAacDF5rZRoIuzHMIxhIqwq4QyM59vhnY7O7/CKfvIgiXbN7X5wEb3L3O3VuB3xLs/2zf15162rdp+x2nMDm8xcDk8Bsf+QQDdvdEXFNGhGMFtwCr3f0bSbPuAS4PX18O/KG/a8sUd7/W3ce4+3iCfftXd38P8DDwznCxrNpmAHevBV4ys6lh07nAKrJ4XxN0b802s+Lw33rnNmf1vk7S0769B3hv+K2u2UB9UnfYEdEZ8K/BzN5M0K8eB37q7tdHW1FmmNkZwCLgOV4ZP/gCwbjJncBRBJfsf5e7dx3cG/TM7GzgM+5+gZlNJDhSGQY8DVzq7s0Rlpd2ZjaT4EsH+cB64P0Ef1xm7b42sy8DFxN8c/Fp4EqC8YGs2tdmdjtwNsGl5rcB/wn8nm72bRis3yPo8msE3u/uS/q0XoWJiIikSt1cIiKSMoWJiIikTGEiIiIpU5iIiEjKFCYiIpIyhYlIiszsb+HzeDP71zR/9he6W5fIQKOvBoukSfK5KkfwnkTStaG6m7/P3UvTUJ5IRunIRCRFZrYvfHkDcKaZPRPeOyNuZl83s8XhvSI+FC5/tpktMrN7CM7Cxsx+b2ZLw/ttXBW23UBwldtnzOxXyesKz1j+enhvjufM7OKkz34k6V4lvwpPTBPJqMRrLyIivbSApCOTMBTq3f1kMysAnjCz+8NlXwfMcPcN4fQHwjOSi4DFZna3uy8ws6vdfWY363o7MJPgXiSV4XseC+edCEwHXgaeILgG1ePp3liRZDoyEcmcNxFc9+gZgsvSDCe4CRHAU0lBAvBxM3sWeJLgwnuTObwzgNvdvd3dtwGPAicnffZmd+8AngHGp2FbRA5LRyYimWPAx9x94SGNwdjK/i7T5wFz3L3RzB4BClNYb/K1pdrR/3PpBzoyEUmfvUBZ0vRC4CPhpf0xsynhTai6Kgd2h0FyDMFtkzu1dr6/i0XAxeG4TBXBnROfSstWiPSB/mIRSZ/lQHvYXfUzgnujjAeWhYPgdXR/W9j7gA+b2WpgDUFXV6ebgeVmtiy8PH6n3wFzgGcJbnT0OXevDcNIpN/pq8EiIpIydXOJiEjKFCYiIpIyhYmIiKRMYSIiIilTmIiISMoUJiIikjKFiYiIpOz/Ay4SCzmkSGOrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    \n",
    "    # Initialize the list of learnable parameters, and fill the initial value with the uniform distribution of [0, 2*pi]\n",
    "    def __init__(self, shape, dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # Create the parameter theta for learning U\n",
    "        self.theta = self.create_parameter(shape=shape,\n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # Create a parameter phi to learn V_dagger\n",
    "        self.phi = self.create_parameter(shape=shape, \n",
    "                                         default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                         dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # Convert Numpy array to Tensor supported in Paddle\n",
    "        self.M = paddle.to_tensor(M)\n",
    "        self.weight = paddle.to_tensor(weight)\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 = U_theta(self.theta)\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        \n",
    "        V = U_theta(self.phi)\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(matmul(U_dagger,matmul(self.M, V)))[i][i]\n",
    "            singular_values[i] = paddle.real(matmul(U_dagger, matmul(self.M, V)))[i][i].numpy()\n",
    "        \n",
    "        # Function returns two matrices U and V_dagger, learned singular values and loss function\n",
    "        return U, V_dagger, loss, singular_values\n",
    "    \n",
    "# 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([theta_size])\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": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:46:13.453107Z",
     "start_time": "2021-03-09T03:46:12.949847Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFZUlEQVR4nO3dd3hUZfbA8e9JDzWUEJpAQHqABJCiIEUEG6hIkVVWlLW31QXr/hTXVVER++pa1oIFbGAXkSoiIBCUFmoA6SWEloS08/vjnYQAKRNIMinn8zzzZObeO3fOhDBn3nLPK6qKMcaYisvP1wEYY4zxLUsExhhTwVkiMMaYCs4SgTHGVHCWCIwxpoIL8HUAp6N27drapEkTX4dhjDFlytKlS/epavjJ28tkImjSpAlLlizxdRjGGFOmiMiW3LZb15AxxlRwlgiMMaaCs0RgjDEVXJkcIzCmNNizZw9jxowhLi6OzMxMX4djDAB+fn60atWKCRMmUKdOHa+eY4nAmNM0ZswY+vTpw9tvv01gYKCvwzEGgLS0NCZNmsSYMWN4//33vXpOhekamha7nUFPfcaiR7oy8KnPmRa73dchmTIuLi6Oa6+91pKAKVUCAwMZOXIkcXFxXj+nQrQIpsVu58EvVvCQfsw5/msZevQjHvyiMgBXxDTwcXSmrMrMzLQkYEqlwMDAQnVXVogWwbPT11IlbR9D/efiJ8pQ/3lUSdvPs9PX+jo0Y4zxuQqRCHYkJnNXwFT8yQDAjwzuDPiCHYnJPo7MmDPj7+9PdHQ0UVFRDBw4kMTERJ/FMmfOHBYsWFBk55s2bRqrV6/OfvzII4/w008/Fdn5T/bVV18xfvx4r45NSkqiVq1aHDp06ITtV1xxBVOmTAFc/O3bt6dVq1ZERUXx2WefZR+3cOFCunbtSnR0NK1bt2bcuHFs3ryZhg0bnvJNPjo6mkWLFjFu3DgaNGhAdHQ0zZs3Z/DgwSf8fs5EhUgE7aonM9R/LoHifsHBksFQ/3lEVU/xcWSmIpkWu53zxs8i8oFvOW/8rCIZpwoNDWX58uWsXLmSmjVr8uqrrxZBpKcnv0SQnp5e6POdnAj+9a9/0a9fv9OOryCDBg3igQce8OrYSpUqMWDAAKZOnZq97eDBg8yfP5+BAwfy+++/M2bMGL788kvi4uL4+uuvuf/++1m6dCkA1113HW+88Ub2v92wYcNo0qQJjRo14ueff84+Z1xcHIcPH6Zr164A3HPPPSxfvpz169czfPhw+vbty969e8/4vVeIRPBCvRkIJ67E5kcmL9b70UcRmYoma5xqe2IyCmxPTObBL1YU6aSF7t27s327O9/GjRu56KKL6NSpEz179sweONy9ezdXXnklHTp0oEOHDtkf3BMnTiQqKoqoqCheeOEFADZv3kzr1q258cYbadu2Lf379yc52bWiX3rpJdq0aUP79u25+uqr2bx5M6+//jrPP/880dHR/Pzzz4waNYpbbrmFrl27ct999zFu3DgmTJiQHW9UVBSbN28G4P3336d9+/Z06NCBkSNHsmDBAr766ivGjh1LdHQ0GzduZNSoUdnfqmfOnElMTAzt2rXjhhtu4NixY4ArP/Poo4/SsWNH2rVrl+uAabdu3Vi1alX24969e7NkyRLeffdd7rjjDgC+/vprunbtSkxMDP369WP37t2nnGfEiBFMnjw5+/HUqVMZMGAAlSpVYsKECTz00ENERkYCEBkZyUMPPcRzzz0HuKnH9erVA1yrrk2bNrmec/LkyVx99dW5/nsPHz6c/v3789FHH+W6vzAqxGBx05RVICd+IwmWdLfdmCIy/L+/5rkvdmsiqRknNvmT0zIY9/UqrohpQMLRVG79YOkJ+6fc3N3r187IyGDmzJmMHj0agJtuuonXX3+d5s2bs2jRIm677TZmzZrFXXfdRa9evZg6dSoZGRkcOXKEpUuX8s4777Bo0SJUla5du9KrVy9q1KjB+vXr+fjjj3nzzTcZNmwYn3/+Oddeey3jx48nPj6e4OBgEhMTCQsL45ZbbqFKlSqMGTMGgLfffptt27axYMEC/P39GTduXK6xr1q1in//+98sWLCA2rVrk5CQQM2aNRk0aBCXXXYZQ4YMOeH4lJQURo0axcyZM2nRogV//etfee211/j73/8OQO3atVm2bBn/+c9/mDBhAm+99dYJzx8+fDiffPIJjz32GDt37mTnzp107tyZlStXZh/To0cPFi5ciIjw1ltv8cwzz2R/iGcZMGAAf/vb39i/fz+1atVi8uTJ2Ylk1apV2b+HLJ07d+bll18G3Df7li1b0rt3by666CKuu+46QkJCGDZsGNHR0bz88ssEBAQwZcoUPv300zz/3Tt27Fio2UF5qRAtAm6ZD+MOwriDJN6xlmMayMLag912Y0rAyUkgS2JS2hmdNzk5mejoaOrWrcvu3bu58MILOXLkCAsWLGDo0KFER0dz8803s3PnTgBmzZrFrbfeCrhvotWrV2f+/PlceeWVVK5cmSpVqjB48ODs7onIyEiio6MB6NSpU/Y3+Pbt23PNNdfwwQcfEBCQ9/fJoUOH4u/vn+97mDVrFkOHDqV27doA1KxZM9/j165dS2RkJC1atABcN8u8efOy9w8ePPiUeHMaNmxYdsvik08+OSXRAGzbto0BAwbQrl07nn322RNaEFmCgoIYNGgQn332Gfv27SM2NpYBAwbkG3uWRx55hCVLlmR/o7/ooosAiIiIICoqipkzZ7J8+XICAgKIiorK8zxFteZ8hWgR5BRWuy6f1ruLaTtr0O5YOpWDK9yvwBST/L7Bnzd+FttzmZzQICwUgJqVgwrVAsiSNUaQlJTEgAEDePXVVxk1ahRhYWEsX7680Oc7WXBwcPZ9f3//7K6hb7/9lnnz5vH111/zxBNPsGLFilyfX7ly5ez7AQEBJwyEpqQUzxhdVsz+/v65jk00aNCAWrVq8ccffzBlyhRef/31U4658847uffeexk0aBBz5szJszUzYsQIHn/8cVSVyy+/PHs6cZs2bVi6dCkdOnTIPnbp0qV07tw5+3GzZs249dZbufHGGwkPD89uWWR1D0VERDBixIh832tsbOwJ5zxdJdoiEJEwEflMROJEZI2IdBeRmiIyQ0TWe37WKO44ml58J78ca8qXy3cU90sZA8DYAS0JDTzxm3FooD9jB7QskvNXqlSJl156ieeee45KlSoRGRmZ3aWgqvz+++8AXHDBBbz22muA6046ePAgPXv2ZNq0aSQlJXH06FGmTp1Kz54983ytzMxM/vzzT/r06cPTTz/NwYMHOXLkCFWrVuXw4cN5Pq9JkyYsW7YMgGXLlhEfHw9A3759+fTTT9m/fz8ACQkJAHmer2XLlmzevJkNGzYAMGnSJHr16lWo39fw4cN55plnOHjwIO3btz9l/8GDB2nQwF1j9N577+V5nt69e7N+/XpeffXVEz60x4wZw1NPPZXdItm8eTMvvPACY8eOBVwizfo2v379evz9/QkLCwNci+a7775jypQpeY4PAHz++ef8+OOPBSYLb5R019CLwA+q2groAKwBHgBmqmpzYKbncbHq2CiMC8MTSZ3zbJE1rYzJzxUxDXhqcDsahIUiuJbAU4PbFekFjTExMbRv356PP/6YDz/8kLfffpsOHTrQtm1bvvzySwBefPFFZs+eTbt27ejUqROrV6+mY8eOjBo1ii5dutC1a1f+9re/ERMTk+frZGRkcO2119KuXTtiYmK46667CAsLY+DAgUydOjV7sPhkV111FQkJCbRt25ZXXnklu2unbdu2PPzww/Tq1YsOHTpw7733AnD11Vfz7LPPEhMTw8aNG7PPExISwjvvvMPQoUNp164dfn5+3HLLLYX6XQ0ZMoTJkyczbNiwXPePGzeOoUOH0qlTp+wuq9z4+fkxZMgQ9u/ff0Iyio6O5umnn2bgwIG0aNGCFi1a8Nprr9GypUv8kyZNomXLlkRHRzNy5Eg+/PDD7C60sLAwunfvTkREBE2bNj3h9bIG45s3b84HH3zArFmzCA8/ZZ2ZQpOS+iAUkerAcqCp5nhREVkL9FbVnSJSD5ijqvl+TercubOe6cI0iz55mq6rnyR+8DdEts/7248xeencubMtkGS88sADD7Bo0SKmT59OUFBQibxmbn+fIrJUVU/pSyrJFkEksBd4R0RiReQtEakMRKjqTs8xu4CI3J4sIjeJyBIRWVIU82bbDriRzIBQIuOnnPG5jDEmP+PHj2f27NkllgQKqyQTQQDQEXhNVWOAo5zUDeRpKeTaRFHVN1S1s6p2LoqmUJXqNfFrPxRWfg7JiWd8PmOMKatKMhFsA7ap6iLP489wiWG3p0sIz889JRVQaswoSEvit69PnTVgjDEVRYklAlXdBfwpIln9/xcAq4GvgOs8264DviypmILO6sTq4PakJx0sqZc0xphSp6Qn0d8JfCgiQcAm4HpcMvpEREYDW4Dch/GLSev75yJ+FeO6OmOMyU2JJgJVXQ7kdvXDBSUZR07i54dmZrItfg1nNWvrqzCMMcZn7KswsOj9h6nzfi/277ELzEzZYmWoi05hylCDK0V9zTXX0K5dO6KioujRowdHjhyhT58+TJ8+/YRjX3jhBW699VY2b95MaGgoMTExtG7dmi5duvDuu+8W8TspPEsEQP2ugwmWNNZOf8PXoZjy7vAueOdiOHxqNcvTYWWoi05hylCDuzgvIiKCFStWsHLlyuy1q0+uIAquimjWFcDNmjUjNjaWNWvWMHnyZF544QXeeeedIn0vhWWJAGjU+hzWBLblrE2TyczI8HU4pjyb+wxsXQhzny7yU1sZ6pItQ71z587sMhTgSl8EBwczZMgQvv32W1JTU7N/jzt27Mi1bEfTpk2ZOHEiL730UkH/vMVLVcvcrVOnTlrUln79uuqj1XT53GlFfm5TPp3yd/i/S069LXrD7Tt2VPWNfqrjwlQfreZ+vtlPddkHbv+Rfac+1wuVK1dWVdX09HQdMmSIfv/996qq2rdvX123bp2qqi5cuFD79OmjqqrDhg3T559/Pvs5iYmJumTJEo2KitIjR47o4cOHtU2bNrps2TKNj49Xf39/jY2NVVXVoUOH6qRJk1RVtV69epqSkqKqqgcOHFBV1UcffVSfffbZ7Niuu+46vfTSSzU9PT3X/W3bttX4+HhduXKlNm/eXPfu3auqqvv3789+/qeffnrC+T799FNNTk7Whg0b6tq1a1VVdeTIkdnvqXHjxvrSSy+pquqrr76qo0ePPuV3NnHiRH3kkUdUVXXHjh3aokULVVV955139Pbbb1dV1YSEBM3MzFRV1TfffFPvvffeU84TGxur4eHh2q1bN3344Yezf9+qqpdeeqlOm+Y+S5566in9xz/+oaqq8fHx2rZt2xPOc+DAAQ0JCTnl/Gcqt89JYInm8plqLQKPqH4jSaQqSYsn+ToUU14d3ApZ1VVUIXHrGZ/SylD7rgx1dHQ0mzZtYuzYsSQkJHDOOeewZs0a4MQFZnJ2C+VGS0G9M6vB7BEUUomPol7m6aUw62Ay9aqH+jokU9Zc/23e+44dgpREjl84r+7x2Z4+78q18n9+HqwMdd4xl0QZ6qzEOXjwYPz8/Pjuu+9o3bo1l19+Offccw/Lli0jKSmJTp065RlvbGwsrVu3Pr03W0SsRZBD3z4DSCGIyYv/9HUopryZ+wzoSYvTaGaRjRVYGeqSL0P9yy+/cODAAQBSU1NZvXo1jRs3BlyC6NOnDzfccEO+rYHNmzczZswY7rzzzkLFX9QsEeTQqFYl7m0QR48FN5CedmYrRxlzgm2LISP1xG0ZqW57EbEy1N4rijLUGzdupFevXtm/h86dO3PVVVdl7x8xYgS///77KYlg48aN2dNHhw0bxl133cX1119fqPiLWomVoS5KRVGGOi+/T3+XDr/ezY5L3qN+lyuK5TVM+WBlqE1pVlrLUJcJUX3/glaJoP6Gj30dijHGlAhLBCfxDwxCYkai66aTtCfe1+EYY0yxs0SQi2MdRqJA7LQXfR2KKcX8/PxIs7EkUwqlpaXhV4himpYIchFcuwlLGv2NKs17+DoUU4q1atWKSZMmWTIwpUpaWhqTJk2iVatWXj/HBouNOU179uxhzJgxxMXFnTA/3hhf8vPzo1WrVkyYMIE6deqcsC+vwWK7oCwfe7bHE7/4G7pe6ds5vqZ0qlOnDu+//76vwzDmjFnXUD62/PRfuv7+Tzav+8PXoRhjTLGxRJCPZgNuJV392P6TrWlsjCm/LBHko2bdxqyq2oM2e74iKemor8MxxphiYYmgAMHdb6QGh/njR6tKaowpnywRFKBl90vZ6teQLWtjfR2KMcYUC0sEBRA/f+b2ncr9BwaxYttBX4djjDFFzhKBFy7vHElooD9f/LrS16EYY0yRs0TghWohgUysP4u/rxzKoUOJvg7HGGOKlCUCL7U4px/V5SiZKz73dSjGGFOkLBF4qVmnCyG8FWGrPvB1KMYYU6RKNBGIyGYRWSEiy0VkiWdbTRGZISLrPT9rlGRMXhOBzjfAjmVsW7XA19EYY0yR8UWLoI+qRucofPQAMFNVmwMzPY9LpZQ2Q0kmiD9n/MfXoRhjTJEpDUXnLgd6e+6/B8wB7vdVMPkJqVqTVX1fp1Xrc30dijHGFJmSbhEo8KOILBWRmzzbIlR1p+f+LiAityeKyE0iskREluzdu7ckYs1V2/OvokZ4PZ+9vjHGFLWSTgQ9VLUjcDFwu4icn3OnusURcl0gQVXfUNXOqto5PDy8BELN24o5n7Fw4nAyM6wGvTGm7CvRRKCq2z0/9wBTgS7AbhGpB+D5uackYzod6Qf+pNuhH1ix6Cdfh2KMMWesxBKBiFQWkapZ94H+wErgK+A6z2HXAV+WVEynq3X/GzhCKMm/vunrUIwx5oyV5GBxBDBVRLJe9yNV/UFEfgM+EZHRwBZgWAnGdFpCKldnaZ1LiNn9FXt276ROhI0ZGGPKrhJrEajqJlXt4Lm1VdUnPNv3q+oFqtpcVfupakJJxXQmIvreSrCksXb6f30dijHGnJFCJQIR6Swiwz1dO1ndPaVhCmqJa9jqHOZUvogftvqTkZnr+LYxxpQJXiUCEYkQkYXAYuAjjk/xnAg8V0yxlXopF7/Ah0c6Mjuu1I9vG2NMnrxtETwP7AZqAUk5tn+KG/StkC5oHUFklQx+n/OZr0MxxpjT5m23zgXABap6wDPYm2Uj0KjIoyojAv39eLbOdDps/5jt2y6jQcPGvg7JGGMKzdsWQSiQmsv2cCCl6MIpexpdeAuBkkGdDZ/6OhRjjDkt3iaCecCoHI9VRPxxNYFmFnVQZUmdyHbQpCeBy9+DTLvS2BhT9nibCO4DbhSRGUAwboB4NXAe8GAxxVZmpHUcBYlb+WPeF74OxRhjCs2rRKCqq4F2wALgRyAEN1Aco6obiy+8ssGv9UASqM6hNbN9HYoxxhSa19cAqOou4NFijKXM8g8MJuCOxfSoXdfXoRhjTKF5ex3BHSJybS7brxWR24o+rLKnmicJHDtWocfOjTFlkLdjBH8H/sxl+2bgnqIKpqxb8vFj7H2qPSnHjvk6FGOM8Zq3iaAhriDcybZ59hmgcr0WNGQ3sTMm+zoUY4zxmreJYBcQncv2jsC+IoumjGvVcwh7pRYhf7zn61CMMcZr3iaCj4CXRORCEQn03PoDLwAfFlt0ZYz4B/Jn5FBiUpeyPm6Fr8MxxhiveJsIHgV+Aabjag0lAd/jppP+X/GEVjY1G3Ab6erHjpmv+ToUY4zxilfTR1U1DRghIo8AMbh1hZer6vriDK4sqh7RmMkN7ue9bXXpdCydKsEVskq3MaYMKdR6BKq6XlU/UdVPLQnkrcVFN7MmNZwvl2/3dSjGGFMgr7+uishwXBXSOpyUQFR1UBHHVabFnBXG0Npbkdlfol3+y0kVW40xplTx9oKyZ4EPgCZAIrD/pJvJQUT4S/1d/CVlCn+uW+7rcIwxJl/etgj+CoxQVVuBxUstL74F3fAqjeI/gZYxvg7HGGPy5O0YgR+wvBjjKHcq1aiHtBkEyz9EU5MKfoIxxviIt4ngDeCUWkMmf8c6XAcpB/l52pu+DsUYY/LkbddQGPAXEbkQ+ANIy7lTVe8q4rjKheCzz2dNaCcC/G2w2BhTenmbCNpwvGuo1Un7tMiiKW9EaH3/LF9HYYwx+fL2grI+RfWCniUulwDbVfUyEYkEJgO1gKXASFXNbX3kMist9RgbVv1G65gevg7FGGNOUagLyorI3cCaHI+fBp5X1bOBA8BoH8RUrFb873bOmnYV+/bbTFtjTOnjdSIQkT4i8oaI/CAis3LeCnGOhsClwFuexwL0BbKmpb4HXOF19GVEnfNGUkVSWPnDW74OxRhjTuHtBWWjcEXmqgK9gb1ADVwZ6tWFeL0XgPuATM/jWkCiqqZ7Hm8DGuQRw00iskREluzdu7cQL+l7DaPOZ3NAUxps+IiMjMyCn2CMMSXI2xbBGOAOVR2BmzH0oKrG4K42PuLNCUTkMmCPqi49nUBV9Q1V7ayqncPDw0/nFL4jwqGokTTXzSxf+JOvozHGmBN4mwiaAlmfYMeAKp77rwCjvDzHecAgEdmMGxzuC7wIhIlI1qB1Q6BcVmprdeFojhLC3sWf+joUY4w5gbeJYD+uWwjcB3WU534tINSbE6jqg6raUFWbAFcDs1T1GmA2MMRz2HXAl17GVKYEVa7Ox9Hvc/ueQexITPZ1OMYYk83bRPAz0N9z/xPcamXvAB8DM84whvuBe0VkAy6xvH2G5yu1Bpzfk0z8mLx4q69DMcaYbN5eUHYHEOK5/xSQjuvq+QT4d2FfVFXnAHM89zcBXQp7jrLorJqV+Fe9BbT79V+k9V1AYIC/r0MyxhivLyhLyHE/Ezf335yGTmc3pM3iOPasnk2d9v18HY4xxng9fTRDROrksr2WiGQUfVjlV+t+f0VDqlNn3ce+DsUYYwDvxwjyqpoWDJSrchDFTYIqIx1GoKu/5ODecjlByhhTxuTbNSQi93ruKnCLiOS8ZsAf6AnEFVNs5VZy+5GELnqdpV+9St/RT/o6HGNMBVfQGMGdnp8C/A3I2Q2UCmwGbin6sMq30AZRLDv7Thq2HeDrUIwxJv9EoKqRACIyGxisqgdKJKoKoOO1hZ5sZYwxxcKrMQJV7XNyEhCRs0UkJK/nmILFr1zIvI/G+zoMY0wF5+2soSdF5DrPfRGRn4B1wE4R6VqcAZZnBxd/xLlrn2bjhnW+DsUYU4F5O2voGmCt5/7FQAegG/A+YF9pT1OT/rcTIJls+el1X4dijKnAvE0EEbgS0QCXAJ+o6mLgZSCmOAKrCMIatiSu8jm02TmNpJQUX4djjKmgClN0rrHnfn9gpud+AHlfY2C8ENBlNHVlP0t/muLrUIwxFZS3ieBz4CMRmQHUBKZ7tkcDG4ohrgqj2XlXsdWvISvWrC34YGOMKQbeJoJ7gZdwq5FdqKpHPdvrAa8VR2AVhQQEMaff1zyz/zz+2Jbo63CMMRWQt9NH01X1OVW9W1Vjc2x/XlVtId4zdEXHswgN9OObeb/5OhRjTAWU5wVlItIRWK6qmZ77eVLVZUUeWQVSLSSQt+t8QvO1Mzl4ZA3Vq1T2dUjGmAokvyuLlwB1gT2e+0ruA8OKqztkzkCDcy4n/IcvOLzuO+g41NfhGGMqkPwSQSSwN8d9U4wadxkIvzai6or3LREYY0pUnolAVbfkdt8UEz9/Mjteh9/sx9kUF0vTVnZ5hjGmZHhbYqKpiNwrIq+IyMsico+IWCuhiKW0G0Ea/myZ+aavQzHGVCAFLlUpIv/ArVPsjxsvECAceFpE7lfV54s3xIqjUs0GbLxsMudG9fB1KMaYCiTfFoGI9ACeAZ4FwlW1nqrWBeoAzwHPish5xR9mxdGsc3+CQyr5OgxjTAVSUNfQrcD7qvrwSQvY71fVB4EPgNuKM8CKaPlXr7Dg6SvIzFRfh2KMqQAKSgTdgHfz2f+u5xhThAJSD3Ju8mxuffx5Fj3SlYFPfc60WFvf2BhTPApKBHWBTfns34grM2GK0Kb6gzimgdyT8TbnyFqGHv2IB79YYcnAGFMsCkoEocCxfPanAsHevJCIhIjIYhH5XURWichjnu2RIrJIRDaIyBQRCfIu9PLr6Xl7+SkzhpayDT9RhvrPo0rafp6dboXpjDFFr8BZQ8ClInIwj31hhXitY0BfVT0iIoHAfBH5HlfQ7nlVnSwirwOjqeCF7HYkJiMBx8cH/MjkzoAveDTxBh9GZYwpr7xJBG8XsN+rEU1VVeCI52Gg56ZAX+Avnu3vAeOo4ImgXfVk+qYsRzwFPYIlnWH+c9mjNXl7diSj+7TxbYDGmHIl364hVfXz4uZ1nSER8ReR5bjrEWbgxhgSVTXdc8g2oEEez71JRJaIyJK9e/fmdki58UK9GchJ+TWADMYEfsKwXy+HJf8jIy2V1TsO+ShCY0x54u16BEVCVTNUNRpoCHQBWhXiuW+oamdV7RweHl5cIZYKTVNWESzpJ2wLkEyoEUmVOk3gm3s49kJH3nrlCX5ZX76TojGm+HnTNVTkVDVRRGYD3YEwEQnwtAoaAjY15pb5ee4SVVg/g+CZ/2Js0FLCm9YC4KNFW0nLyGRIp4ZUDvbJP6sxpowqsRaBiISLSJjnfihwIbAGmA0M8Rx2HfBlScVUJolAi/743zyPejd+QoC/HxzcRtdZw5j7zSS6P/UTT323hh2Jyb6O1BhTRpTkV8d6wHsi4o9LQJ+o6jcishqYLCL/BmIpeHDaAPj5QaWa7v6hnTSrlMz/jk1gU0gbHpl/BT3nt+PiqLqM7hFJTKMavo3VGFOqiZvMU7Z07txZlyxZ4uswSpeMNFj+Icx9Bg5tZ1O1cxh88F4SjykdG4UxukdTBrSNcC0IY0yFJCJLVbXzydvtU6G88A+ETqPgzmVw0dM0jerO/If68+jANgQc2sLtHy0j9s9EX0dpjCmFvOoaEpEauPn9fXCVR09IIKpap8gjM6cnMAS63QJAFeD6yERGzbiNfc37U7tSU6AmT/8QR6YqD17c2qehGmNKB2/HCN4H2uIu+NqNlxeRmVKgZlOk132E//oq/Kc7tB+Gf+pgEoPqZx+ycvtB2tavhkhuS1IbY8o7r8YIROQw0EtVlxV/SAWzMYLTcHQ//PICLH4TAkPQe1YjQZVYse0gA1+ZT/uG1RndI5JL2tUj0MYRjCmXznSMYGMhjjWlUeVa0P9xuHs5XPkGElQJVGm58W2evbgeR1LSuXvycno+PZv/zNlAYlKqryM2xpQQb1sEvYB/AmOAlaqaUdyB5cdaBEVk5+/wRh8ICEa73MzPdUbw398O8MuG/YQG+jOkU0OuP68JTcOr+DpSY0wRONMWwQZcSeplQKqIZOS8FWWgpgTV6wC3L4aWlyC/vMD53/Xjw7Pn8sNtnbmsfT2m/PYnfZ+by71Tlvs6UmNMMfK2RTAPqAG8Ti6Dxar6ebFElwdrERSD3atg9pOw8w+4cwkEBLP3UAofLNpKSKA/t/ZuRmam8t3KnVzYJoLgAK9rDRpjSom8WgTezhrqDHRR1ZVFG5YpNSLawtUfQspBCAiGtBTCJ1/MPR1GQKfrAFi4aT93fBTLi1dHc3n08SKx02K38+z0texITKZ+WChjB7Tkiphci8gaY0ohbxPBaqBacQZiSomQ6u5n0j4IrATfj4UFL0Gv++jeYQQfjO5Kl0hX2uJ/8+P5cdUuYv9M5Fh6JgDbE5N58IsVAJYMjCkjvB0j+CcwUUT6iUiEiNTMeSvOAI2PVG8Io76BkVOhSh346k7k1a70qJtOUID7s0lMSmVhfEJ2EsiSnJZhy2oaU4Z42yL4zvPzR04cHxDPY+swLo9EoFlfaNoH1n4Pcd9ClQi3b+867r2wBS/P2pDr1YVW/dSYssPbRNCnWKMwpZsItLrE3QCO7oM3ekPt5lxZdSBfHG5FOIm8EvQyd6TexV7CUOCeKcu5tXczWkRU9WX0xpgCFNg15Flo/hlgl6rOze1W/GGaUiUkDC55FpITmJj2OJ8FP86/A/7HObKWOwO+IDjAj55n12LG6t3sPXwMgANHU0lKTc//vMYYnyiwRaCqaSISidUXMln8AyDmGmg3FJa9R9ufniQ0NQGAYQHzCB/wf1zcPZqjx9KpFOR6DV+ZvYGpsdtZ8EBfQgKtJ9GY0sTbweL3gBuLMxBTBgUEQZcbCY0aCH7uO0WIP1y89Cb46TEqJ23LLmR3Wft63H1B8+wk8MS3q/n2j52kZWTmeXpjTMnw9oKy/wDXAPHAUuBozv2qelexRJcHu6CsFDm8C17sAOkpx7eJH2T9XZ3dDzrfAM37u5YEcDgljYtf/JltB5KpUzWYq7s0YkSXs6hXPdQHb8CYiiOvC8q8TQSz89mtqtr3TIIrLEsEpcg390LsJMjIUaTOPwjaDoYaTWDZe3B4J1zwKPS8N/uQjExl7ro9TPp1C3PW7cVPhH6t63Btt8ac16w2fn5WEtuYonZGVxarqs0aMrnbtvjEJADu8Z5VMPi/cP5YWPcDNOjo9q39HmI/wP+c0fRt0Zu+rSL4MyGJDxdt5ZMlfzJ91W4ia1fmmq6N+EvXRlQKKslltY2pmAq1ZrGIhABn4waON6pqSgFPKRbWIijDYj+EGf8HSfuhRiR0vh6ir4HKtTmWnsH3K3YxaeEW1u06zMKHLqBycAAHk9KoXinQ15EbU+adaddQIPAkcAcQhLuQ7BjwMvCwqqYVbbj5s0RQxqUfgzVfw5L/wZZfoHZLuH2Ru17BY8/hFOpUDUFVGfDCPDo0DOPZoR18GLQxZd+ZFp17GhgB3ALM92zrCTyFm3k0piiCNBVEQDC0G+Jue9a4AWcRSEuBSVdC2yuo0344EEJGpnJN18bUqx4CwMGkNF6ZvZ4RXRrZOgnGFBFvWwS7gBtU9buTtl8KvKWq9YopvlxZi6CcStgEn42GHctcwbuoq9yMo6zxBWBW3G5uen8p6ZlKz+a1uaZrY/q1rkOALa9pTIHOtGsoGYhW1bUnbW8FxKpqic77s0RQzu2IhSXvwIpPIS0JbpoD9WOyd+85lMKU3/7ko8Vb2XkwhbrVQvhL10Zcfc5Z1KkW4ru4jSnlzjQRLASWqurtJ21/DZcguntxjrOA94EI3GDzG6r6oqd66RSgCbAZGKaqB/I7lyWCCiLloJtl1H646zqa8QikJrkB5oi2pGdkMituD5MWbuHn9fsI8BMGtK3Ltd0a071ZLV9Hb0ypc6ZjBPcB34lIP2ChZ1s3oD5wsZfnSAf+oarLRKQqsFREZgCjgJmqOl5EHgAeAO738pymPAupDh2uPv742BGI/QB+exPO6kZA5xvo3+Zy+retS/y+o3y0aAufLt3G0dT07ESQkpZhJS2MKYDX00dFpD5wO9DKs2kN8B9V3XFaLyzyJfCK59ZbVXeKSD1gjqq2zO+51iKowI7uh98/cjOOEjZB11vh4vHZu1PSMkg4mkr9sFC2JyZz0fPzmDCsAwPa1gVsNTVTsZ1piwDPB/7DRRRMEyAGWAREqOpOz65duK6j3J5zE3ATQKNGjYoiDFMWVa4F594J3W6HzfOg+llu+9ZFMOdJQjrfQP2Wrly2qnJJu3q0re8W15s4Yy2vzdlIWob78mOrqRnj5Nsi8Hb1MVVN8PoFRaoAc4EnVPULEUlU1bAc+w+oao38zmEtAnOKuO/gu7FwaJtbPKfjX6HjdRB2VvYhbR75gaTUjFOe2iAslF8eKNEqKcb4RF4tgoLm3O0D9hZw21OIIAKBz4EPVfULz+bdni4hPD+9Pp8x2VpdAn//A0ZMgXrRMG8CvNkHMo6vgZCcSxIA1zL4ef1eUtOtEqqpmArqGsqvxtBFwN24QeACiatH/DawRlUn5tj1FXAdMN7z80tvzmfMKfz8oeVF7pa4FfaudRVPMzNh0hU8UKUhbx05D9ATVlMDGPn2YqqGBHDz+U25o29zn74NY0pavokgt9XHRCQGeBZ3ZfF/gce9fK3zgJHAChFZ7tn2EC4BfCIio4EtwDAvz2dM3sIauRtAsuu5vDn9Q24Insx2rUUj2ctdAZ/zpNzEvy5vS41KQfy4ehe1qwQD7grmuybHcne/5nRslG9PpTFlnteDxZ5Vyp4AhgJfAG1UdaO3z1fV+bgaRbm5wNvzGFNolWvDdV/Bvg38+fXTRG75BAGGB8yl9oBHuLhtVQgIoV+b4/MU/jyQxNaEpOw/2NitB5i/fh/929alRUSV7AV3jCkPCkwEIlILeARXZ+gX4FxV/a24AzOmyNU+m6bhVWFbEGSkEuTvx8X734cFNWDR624RnVaXQYv+RDWozqx/9Mp+6qL4BJ6bsY7nZqyjUc1K9G8TwYVtIujcpCb+tnaCKeMKmjX0MDAWd8XvA6r6QwnFlS+bNWROS26rqQWEwJB3YP2PsPY7OLIb/AKhxQAY/sEJFVF3H0rhpzW7mbF6Nws27Cc1I5OalYPo26oOF7aJ4Pzm4YQG2cVrpvQ6rRITIpIJJAOzgTynVKjqoKII0luWCMxpyWs1tZiRcNlEN6i87TeI+wYy0o5fqDbtNqjdwrUWap8NwJFj6cxdu5cfV+9iVtweDqek07peNb6/uyfgZihZUjClzeleUPY+ri6QMWVfXqupbVvs7vv5QaOu7pYlNQn2rIblH8JPj0J4K2h1KVXaX82l7Vtwaft6pGVksjg+gaPH3AS6tIxMzh0/k9E9Im0GkikTCpo1NKqE4jCm+N0yv+BjThZUyVU/TfzTdR3FfQPzX3Crq4W3gCN7CdyzivMizwN/t4rasfRMRnZvQsfGbrbR2l2HufPjZVzYJoIL29SlfYPqtiazKVUKtVRlaWFdQ8ankhJcl1JwFVj8Jnw3BkLCoMVF0OpSOPsCCKqcffjyPxN5+vs4Fm9OICNTiagWTL/WEfRvW5fuTWsRFGBrKZiScUZlqEsbSwSm1EhNgo2zIO5bWPc9JB9wi+r8I85VT83MdF1OwIGjqcxeu4cfV+1m3vq9JKVmUCU4gF4tw+nfJoKLo+pZUjDF6oyLzhljchFUCVpf5m4Z6W4N5p2/uyQAMPkvkHoEWl1GjVaXMLhjIwZ3bEhKWgYLNu7jx1W7+WnNbuav38cl7dxCf0u3JFCveij1w0KtWqopEdYiMKY4zX0GVn4Be9e4x/U6QJebIeaa7EMyMpU/E5JoUrsyqkrf5+bSsEYoV3VsyINfrCA57XiNpNBAf54a3M6SgTkt1iIwxhd63edu+ze6geY137hrFQBSj8Kc8fi3upQmDc8BQER467rOJKdmcPOkpSckAYDktAye/iHOEoEpUtYiMKakqboL1bYsgPcGQWYaVA6Hlpe4axWa9oKAYCIf+BYFwjlwSpG8c5rU4Pzm4fRqGU5UfZuFZLxzumWojTFFLetq5cbnwn2b4Kq3oUkPWPk5fDQU9q0DoH31ZKqQxF0BUzlH1nJngKvcXiU4gJS0TJ6bsY5Br/xC5yd+YtlWt8x3WfxiZ3zPuoaM8aWQatBuiLulH3ODzRFRALwU8Q31Ur7EH8VPlGH+c3lThvKPK3pyRUwD9h05xvz1+5i7bi9Na7vpqm/+vIlpsTv4/NZz7cpm4zVLBMaUFgHB0Oz4SmmN+99J4ierqX5oLQDBksb31Z6gSoxbXrN2lWCuiGlwwnhBveqhtK1fLTsJ3DNlOYdT0unVMpxezcNpVKtSCb4hU1ZYIjCmtKregLCkLdkPBaiSvAMO73ZjCh8Pd91LrQdBrWYADOxQn4Ed6mc/p07VYJZsSeCnNW6AukmtSvRqEc75LcLp1rQWlYPtI8DYYLExpVd+RfJ63QcfXw07Yt32iChoPRCi/3J8QR4PVWXz/iTmrt3DvPX7+HXjfpLTMgj0F85pUpNruzXOvobBlG82fdSYsia/InlV63pqIG11U1LXfAVzxkPDzi4RJG6Fo/ugfgwiQmTtykTWjmTUeZEcS89gyeYDzF23l3nr9rIjMRmAhKOpPPndGm7s2ZSWdauW/Ps1PmMtAmPKi8O7oFItV/zup8dg/kSofpZrKbQeCGd1des6n0RVERGWbE5g9HtLePf6c4hpVIMlmxOYt34fvVrUpkPDMAL8bZJhWWe1hoypSJISYN0PsPorVwsp45irmHrnsuzaR7nJyFQE8PMT3pi3kfHfx5GpUC0kgB7Na2ePL9SrHlpy78UUGUsExlRUxw67FdgO74but7lt717mWgttBkHTPhAYkutTE5NS+WXDfuau28O8dfvYdcit7tYiokr2BW09zq6NiFhdpDLAEoExxklPha/uhLXfw7GDEFQFmveHLjdB4+55Pk1VWbf7CPPW7WXe+r0sik+gUc1K/HRvL6bFbue+z/4gNeP4QoZWF6n0scFiY4wTEASD/+sSwuZ5sOZrV0b77H4uERzaCZvmQMuLILRG9tNEhJZ1q9KyblVuPL8pyakZbPcMND8zPe6EJACuLtJjX6+ieUQVmoVXISTQLnArraxFYIyBzAx3CwiCJf+Db+4BvwCIPN9dp9DqUqhSJ8+nZ9VFyoufQJPalWlRpyot6lalRUQVWkZUpUntygTaIHSJsRaBMSZvfv7HZxR1HAV1O8CaL91g8zd/d6uwjd3gWgjpqS5h5FA/LDS7dZBTnarB/N9lbVi/+zBrdx9m3e7D/Lh6F5merPHYoLZcd24Tdh5M5tMl2xjcsQENa9jVzyWtxBKBiPwPuAzYo6pRnm01gSlAE2AzMExVD5RUTMaYXPj5QcNO7tbvMdi9CrYvPd5N9PFwSDnomZbqrmoeO6AlD36xgipp+7IrpR4JrMVDl7Q+4UpngJS0DDbuPcL63UeIaRQGQNzOw0ycsY7eLcNpWKMSX/++g9fmbKRl3ao097QeWkRUpUFYqFVaLQYl1jUkIucDR4D3cySCZ4AEVR0vIg8ANVT1/oLOZV1DxvjQgpfdYjs7lrnHddpCt1uZJn3J/OZerkifztSAAfhfNrFQA8VJqekE+fsR4O/HrLjdvLdgC+t2H2bnwZTsYyoF+dO8ThWaR1SlZURVru5yFlVDAov6HZZbpWLWkIg0Ab7JkQjWAr1VdaeI1APmqGrLgs5jicCYUiDxT7fYzuqv3MBy++HwYntXRdU/CO5YCjUaFXyeAhxKSWP97sOs232EtbsOs37PYdbuOsL+o8dY9dgAKgUF8PLM9fyycR8f39gNEWHT3iNUDQmkdpUgRHJvQVTE6a6ldYwgQlV3eu7vAiJ8GYwxphDCzoJut7qbKnz7DzfgDK4UxkvRENkTInu5GkhV657Wy1QLCaRT45p0alzzhO2JSalUCnIfYTUqB9EgrFL2h/5DU1ewcFMCNSoF0sLTrdSiblVa1KlCi4iqzF2394RlQLcnJvPgF66qa3lPBrnxdYsgUVXDcuw/oKo18njuTcBNAI0aNeq0ZcuW3A4zxvjC4V3wYgdIP96Ng/hDzWawfx3cthDqtIbNv8Ce1dC0N9Q6+/giPUVscXwCK7cf9LQeDrN+9xEOH0vP3u8nZA9Y59QgLJRfHuh76o5yorS2CHaLSL0cXUN78jpQVd8A3gDXNVRSARpjvDD3GdATryPAz9+1CK7/1pXNBne9wsJX3f2q9d2ynJG9XLdSPqUvCqtLZE26RB5vQagqOw+msM4zc+nJ7+Jyfd6OxGRUle2JyRVq9pKvE8FXwHXAeM/PL30bjjHmtORXKTXn9QcDnoBzRkP8XNg0F9ZNh60LIXqE2x/7AYSEuaU7Q8OKLDwRoX5YKPXDQundsg7vLdiS63TX+mGhbNmfRO8Jc2gQFkrXyJp0a1qLrk1r0qhmpTzHG8q6kpw19DHQG6gN7AYeBaYBnwCNgC246aMJBZ3LBouNKScyM+HILqhW340zvNjeldAWP6gf41oLrS515bWL0LTY7SeMEcDxkhjntwjn6993sCh+P4s2JbD/qEtwdauF0LVpTbpG1qJb05pE1q5c5hJDqZg1VFQsERhTTqWnwvYlrsTFprnufpeb4aInISMNfn0FmpwP9aNzLaldGN7MGlJVNu49wsJNCSzctJ9F8QnsPXwMgMk3daNb01psT0wm6Vg6Z9epUuoTgyUCY0zZc+ywm45auTbsWA5v9HLbg6u77qOmvV0F1dOckVRYqkr8vqMsik/gypgGhAT688wPcbwxbxO/P9qfysEBxG49QGiQPy3qVC11F7+V1sFiY4zJW3BVdwPXChizHuLnuRZD/FxY+62bjVS1rrsCeucfbgC6Wv38znraRISm4VVoGl4le9tfujYi+qyw7PWfx38fx6L4BMIqBdKlSU26Nq1F18iatK5XDf9SlhiyWIvAGFN2JcRDtQau9tGsJ2DeM257reYuITTtDS0ucqu2lZA/E5JYFJ/Aok37WRi/nz8T3KB0tZAAukS6MYYezWvTul61Eospi3UNGWPKt8xM2LPKjS1smgNbFoB/ANwX78YT4r6DwFBo1M39BHf9w2fXw5B3oWrxXM+6IzE5e+B5UXwC8fuOcnl0fV68OgZV5b0Fm+nZIpxmOVoZxcUSgTGmYklPhQPxEO6pWvOf7u5iNv9gOKuLazHsWuHWY+h0PVw2sUTC2n0ohZS0DBrXqsz2xGTOGz+Lxy9vy8juTdh9KIXPlm6jW9OatGsQRlCAu7aiqMphWCIwxlRsxw7Dll891zDMgd0r3TRVzYSAEOj1ALToD+Gti/TitoLsPXyMoAA/qocGMn3VLm6etBSAkEA/OjWuQbWQQGbG7SE1/cxXf7NEYIwxOU29FVZ8Cplp4BfofgJUqu2pkXQ+tLgYqtUr0bASjqayOH4/Cz1dSWt2Hsr1uNMph2GzhowxJsvhXbDqi+Mf/plpEBAMF4yDnb+7VsOqqTC0BrS9Eg5sdq2JyPOhevEWpatZOYiLoupxUZRLQHmt/rYjlyujT5clAmNMxZNbbSRV2L/BreesCvs3Hh9AXvsD/OBZKqVmM5cQmvZyM5KyBp6LSV6rv9UPK7rXtcVCjTEVT361kcBVRa199vFrGLrcBLf8AgOehNrNYcVn8Pnfjpfd3vATrP0eUnLvxjkTYwe0JDTwxKuoQwP9GTugwKVbvGYtAmNMxXPL/MId7+cHdaPcrfvtkJEO+9ZBsGfK5y8vue4k8ffUSDofzr7AXf18hrIGhItzER0bLDbGmDOVlgLbfnNXPcfPczWSIs+HkVPd/t/ehoi2UL+ju/jNR2yw2BhjiktgiGemUU/gYTh2BJL2uX0ph+C7saAZEFgZGnd3SaLVZVCrmU/DzmKJwBhjilpwlePdRiHVYOwG2Dz/eIthxiMQWMklgsO7YfWXLjmEtyy2VdvyY4nAGGOKW6Warkpqm0Hu8eFdbroqwJZf4Pux7n6VCJcQIs+HNpdDSPXj5yjGchiWCIwxpqTlLJsdNRgadDzeWoif5y50a3aBSwTx81wS2DjLreY29+kiL4dhicAYY3ytRhN36/jX49cwZF24FvsB/DHl+LHLP4Re9xdpq8CuIzDGmNIk6xqGLFe8Bq0vd1NTwV0IN/fpIn1JSwTGGFOaHd0L66e7WUfgLnxb/qEbZC4ilgiMMaY0y7UcRtG2CiwRGGNMaVZQOYwiYIPFxhhTmhW2HMZpsBaBMcZUcJYIjDGmgrNEYIwxFZwlAmOMqeAsERhjTAVXJtcjEJG9wJbTfHptYF8RhlPcylK8FmvxKUvxlqVYoWzFe6axNlbV8JM3lslEcCZEZEluCzOUVmUpXou1+JSleMtSrFC24i2uWK1ryBhjKjhLBMYYU8FVxETwhq8DKKSyFK/FWnzKUrxlKVYoW/EWS6wVbozAGGPMiSpii8AYY0wOlgiMMaaCqzCJQET+JyJ7RGSlr2MpiIicJSKzRWS1iKwSkbt9HVN+RCRERBaLyO+eeB/zdUwFERF/EYkVkW98HUtBRGSziKwQkeUissTX8eRHRMJE5DMRiRORNSLS3dcx5UVEWnp+p1m3QyLyd1/HlRcRucfz/2uliHwsIiFFdu6KMkYgIucDR4D3VTXK1/HkR0TqAfVUdZmIVAWWAleo6mofh5YrERGgsqoeEZFAYD5wt6ou9HFoeRKRe4HOQDVVvczX8eRHRDYDnVW11F/0JCLvAT+r6lsiEgRUUtVEH4dVIBHxB7YDXVX1dC9WLTYi0gD3/6qNqiaLyCfAd6r6blGcv8K0CFR1HpDg6zi8oao7VXWZ5/5hYA3QwLdR5U2dI56HgZ5bqf2GISINgUuBt3wdS3kiItWB84G3AVQ1tSwkAY8LgI2lMQnkEACEikgAUAnYUVQnrjCJoKwSkSZADLDIx6Hky9PVshzYA8xQ1dIc7wvAfUBmAceVFgr8KCJLReQmXweTj0hgL/COp9vtLRGp7OugvHQ18LGvg8iLqm4HJgBbgZ3AQVX9sajOb4mgFBORKsDnwN9V9ZCv48mPqmaoajTQEOgiIqWy+01ELgP2qOpSX8dSCD1UtSNwMXC7p5uzNAoAOgKvqWoMcBR4wLchFczThTUI+NTXseRFRGoAl+OSbX2gsohcW1Tnt0RQSnn62j8HPlTVL3wdj7c8XQGzgYt8HEpezgMGefrdJwN9ReQD34aUP8+3QVR1DzAV6OLbiPK0DdiWozX4GS4xlHYXA8tUdbevA8lHPyBeVfeqahrwBXBuUZ3cEkEp5Bl8fRtYo6oTfR1PQUQkXETCPPdDgQuBOJ8GlQdVfVBVG6pqE1x3wCxVLbJvVkVNRCp7Jgzg6WbpD5TKmW+qugv4U0RaejZdAJTKCQ4nGUEp7hby2Ap0E5FKns+HC3Bjh0WiwiQCEfkY+BVoKSLbRGS0r2PKx3nASNy31aypbZf4Oqh81ANmi8gfwG+4MYJSPy2zjIgA5ovI78Bi4FtV/cHHMeXnTuBDz99CNPCkb8PJnye5Xoj7hl1qeVpZnwHLgBW4z+4iKzdRYaaPGmOMyV2FaREYY4zJnSUCY4yp4CwRGGNMBWeJwBhjKjhLBMYYU8FZIqjARGSOiLzio9dWERnii9f2RmmPryh4qliO8+K42SLy1xIIKbfXzvffwVP19qqSjKk8skRQTnku8vqPp4TxMRHZLSIzReTCHIcNBh70VYxFTUQ6ej44euaxf4qILCjpuPKT1wediLxbGkpki8ilwFnAhzm2bfbErSKS7Ck5PdZzoVNJexwYLyL2WXYG7JdXfn2OK0UwGmgBXAZ8D9TKOkBVEzzVTcscEQk4+YPHU7F1OXBDLsfXAq7AKo4W1t3Au6qacdL2f+EuJGyNK4b2JOCLgnjfAVVxZSLMabJEUA55yj30BB5Q1ZmqukVVf1PVCao6OcdxJ3QNeb7p/VNE/utZpGObiIw96dwtRGSuiKSIyFoRuUREjojIKM/+Jp5vip1Pel5BTfzxnvMle+J4JufCGyIyztOVMUpENgLHgNwqW74FDPUU7MvpWs9zpojIRSLys4gcEJEEEZkuIq3zic2r9yQiDURksue8B0TkWxFpntd5C0NEBovIH57fT4Ln3yAix/6B4qqTpohIvIg8Ia6YWtb+OiLypef5W0TklGSZy2uG42rcfJ3L7sOquktVN6vqW8AfuPIXWc9t5nm9XSJyVESWiSv4l/P8Bf695RLT/SKyT0S6gSt2iEsGIwp6PyZvlgjKpyOe2yAp/CpG9+AuYe8IPA08I55VpjzN76lAOtANGAU8CgQXQcxHcd/kWwO34eoAPXzSMZHAX4ChQAcgJZfzfAj4A8NP2j4amKKqR3EJ5AVci6k3cBD4OucHZ2GJSCVcsb0UoBfQHVcu+CfPvtMmInVxBfLew/1+zgcm5dg/APe+XwHa4n6PQzixvMO7wNm4D/YrgL8CTQp46R645JlnbSNxenviSsuxqwquBXoh7t/qc+ALEWl10iny/HvL5XUm4EpY9Dpp0aPFuN+5OV2qardyeAOuwi3Ek4KrsTQBt/pSzmPmAK/keLwZ+PikY9YD//TcH4BLAg1y7D8XVy9/lOdxE8/jziedR4EheT3OJf5bgA05Ho/DfdBEePHePwAW5Hh8juf1uuZxfGUgA1fu+ZT4vHlPuA/f9XjKtni2+QP7gWH5xJrr7wH3wf2N535Hz3GN8zjHPOD/Ttp2Be7LgOC6BhU4L8f+xp73PC6f2P4ObMll+2ZcgjgCpHrOnQycW8C/y8KsvyVv/t5y/H6GA+8A63L7HeBKSGcCASX9/6y83KxFUE6p6ue4uuUDcd/MzgUWishDBTz1j5Me7wDqeO63Anaopyyyx28UwQIvIjJEROZ7uhKOAM8DjU46bJt6Vyr4LaB7jm+fNwAr1VMe2dNt8ZGIbBSRQ8BuXOv45NcrjE64FsthT1fZEVxLowbQ7AzOC/A78BOwUkQ+F5FbPd02OV/74azX9bz2R7gEVxf3bT0T980ZAHUrcRW0wlUoube6ACbiisr1wrWEHlPV7IF4cVVTnxG37vYBT0ydOfV3nN/fW5YJuJZbD819BbFkXMIrsjV8KxpLBOWYqqao6gxV/ZeqnosrbT2ugC6QtJMeK4X7O8lKCtkDueLWVsiTp793MjAdl7higH/ilrzM6aiXMcwFNgA3iCuLPQLP8oke3wDhwM1AV8/rpQN5/V68eU9+uIHq6JNuLYD/5hPrYaB6LtvDcIkEdf3g/T23P3DdXOtFpEOO137spNdtDzTHrRiWpbAVJvfhEllu9qvqBlX9Fdf6HCMifXLsn4Drwvs/XLKIxiWik3/H3vy9zcAltLwq8NYEUvT4cqmmkAJ8HYApUatx/+YhuCZ9YcUB9UWkvqpmfZvszIn/cbM+eOrl2BZdwHnPA7ar6uNZG0Sk8WnEB7g1lEXkf7gZL3G4b7aTPOethWvZ3Kaqsz3bOpL//wVv3tMyXMLZp4Vbp3ct7ht9dqISt5B6B1x3SPZ7wnXx/Soi/wJW4bpMfve8ditV3ZDbC4hIHO7fqAuwwLOtEa7FmJ9YIFxEaqvqvrwOUtUD4iYdPC8iMZ5YewDve1qmeMaqmuG6dwrrO1yZ6E9FRFX1vZP2R+F+B+Y0WYugHBKRWiIyS0SuFZH2IhIpIkNx6/TO1NNf9nIG7oPrPRHp4PkmPxH3bdp16Kom4/qC7xeRtiJyLu7bYX7WAQ1E5BoRaSoit3Lms0DeA2p7Xnuaqu73bD+A+6Z7o4icLSK9gNc97yFXXr6nD3FdTF+KSC/P7/x8EXmugJlDE3Etl9vFzciKxtWZr+n5iYh088yuOcfzAT4IN7c/a9GXfwF/EZF/iUiUiLTydLU944l/LfAD8F8R6e55jXdxXSr5icWtQd2jgOMA/gO0xLUCwP2bXinu2o52uHGb0+66Ube+xVDgdTn14raeuPdnTpMlgvLpCO6D625cN8kq3AySjzh1No3XVDUTuBI3S2gx7sP2CVwSyNmXnDU18Tdct8g/Czjv18CzuJk8f+BmmjxyunF6zrkD902yBjmuHfC8h+G4rpOVwKu47otjBZwy3/ekqkm42TybcGvfxuF+PzVwySevOD8GrvfcluA+0OoCPdWt+AWui+g8XJfWeuA54HFV/cBzjunApUAf3L/LYtxawVtzvNQoIB6YhZsO+hFusDZPni6p/wHX5Hec59g9uFbXOM/ssntxSeRn3BjVQs/90+ZJBsNwCe2v4Kbs4sa/3snvuSZ/tjCNOSOefurluBk1ZWlBeOMFEamDa3mco6rxvo7nZCLyLFBdVX1xMVu5YWMEplBE5ErcoO163LTKiRzvpzbljKruEXfxWSNci6K02UPBXY+mANYiMIXiaZL/E9dHfQB3LcI9Xk7rNMaUQpYIjDGmgrPBYmOMqeAsERhjTAVnicAYYyo4SwTGGFPBWSIwxpgK7v8BsXZ8zNDNnSUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "singular_value = singular_value_list[-1]\n",
    "err_subfull, err_local, err_SVD = [], [], []\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "\n",
    "\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": 11,
   "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": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:15.837676Z",
     "start_time": "2021-03-09T03:47:14.960968Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsUlEQVR4nO3dbWxcVXoH8P8fx46NHSW2A4mdOC+khGAqNomsQEVAdOkugS8EtOLlA2Il1KwqkIq0+wFRbZdWVctWBUSliioUtNmKQtgFRFqhdilaKUVasWvSxCREJCE4L7bjJMRJnMRxYufph7nZOtF9jsczd2Zsn/9Psjw+z5yZx9fz+I7v8TmHZgYRmf6uqXQCIlIeKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFil6KR/D7JUZJnxnzcXem85EozKp2ATBu/MbO1lU5CfDqzT3Mku0n+iGQXyVMkN5OsrXReUn4q9jg8DGAdgKUAbgXw/bQ7kVxL8mTgI3TmXkXyOMk9JH9MUu8aJxn9QOLwj2bWCwAk/x3AyrQ7mdknAOYU8PhbAfwhgAMAbgGwGcAIgL8r4LGkRHRmj8ORMbfPAWjI8sHNbL+ZfW1ml8zscwB/DeB7WT6HFE/FLr9H8s6rrqhf/XFnng9lAFjKXGXi9DZefs/M/gcFnPVJ3gdgm5n1k1wB4McAfpF1flIcndklC/cA6CJ5FsCHAN4D8LeVTUmuRi1eIRIHndlFIqFiF4mEil0kEip2kUiUdeitsbHRWltby/mUmbrmmvTfjV47AJD+cHPo4ujo6GhBsUIuuBaafyh26dKlCbUD4dxnzPBfqlVVVW7MU2gek11vby8GBgZSfzBFFTvJdQBeAVAF4F/M7IXQ/VtbW7F58+ZinrKiZs6cmdpeX1/v9qmpqXFjFy5ccGOnTp1yY6dPn3Zjw8PDqe2hgq6rq3Nj3vcMhAtwaGgotX1wcNDtMzIy4saam5sLinmFe+7cObfPxYsX3dhk98gjj7ixgt/Gk6wC8E8A7gPQDuAxku2FPp6IlFYxf7OvAbAv+b/oCwDeBvBANmmJSNaKKfYFAA6N+fpw0nYFkhtIdpLsHBgYKOLpRKQYJb8ab2YbzazDzDoaGxtL/XQi4iim2HsAtI35emHSJiKTUDHF/jsAN5JcSrIGwKMAtmSTlohkreChNzMbIfk0gP9CbujtDTPblVlmIpKposbZzexD5KY0isgkp3+XFYmEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4mEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEUTvCkOwGMAhgFMCImXVkkZSIZK+oYk/8sZkdz+BxRKSE9DZeJBLFFrsB+BXJz0huSLsDyQ0kO0l2DgwMFPl0IlKoYot9rZmtBnAfgKdI3nX1Hcxso5l1mFlHY2NjkU8nIoUqqtjNrCf5fBTA+wDWZJGUiGSv4GInWU9y1uXbAL4LYGdWiYlItoq5Gj8PwPskLz/Ov5nZf2aSlYhkruBiN7P9AL6VYS4iUkIaehOJhIpdJBIqdpFIqNhFIpHF/8ZH49SpU6ntPT09bp/z58+7sZkzZ7qx2bNnu7Frr73WjY2Ojqa29/f3u33OnDkz4ccDgAsXLrixwcHBCedx+PBhN1ZVVeXG2tvb3diKFStS2xcuXOj2qaurc2NTmc7sIpFQsYtEQsUuEgkVu0gkVOwikSj71XgzK/dTZubIkSOp7bt27XL7dHd3u7Ha2lo31tTUVFA/78p6b2+v2+fYsWNuzBuBAIChoSE35uWYzKVIdfDgQTe2f/9+N7ZkyRI3tn79+tT2e++91+2zYMECNzaV6cwuEgkVu0gkVOwikVCxi0RCxS4SCRW7SCQ0EWYCzp07l9p+6NAht09XV5cbO336tBsLTQoJDYd5Q17XXXed26ehocGNhSbJjIyMuLGlS5dOqH28PELHOBTzJt4MDw+7faYrndlFIqFiF4mEil0kEip2kUio2EUioWIXiURZh95IoqamppxPmam2trbU9jvvvNPts2jRIjfmzaIDgJ07/Z20+vr63Nj8+fNT2++44w63z2233TbhxwOAa67xzxXe8OCePXvcPlu3bnVjt9xyixsLDeetXbs2tb2lpcXtM5Vfo6FZheOe2Um+QfIoyZ1j2ppIfkRyb/JZ27OKTHL5vI3/GYB1V7U9C+BjM7sRwMfJ1yIyiY1b7Ga2FcCJq5ofALApub0JwPps0xKRrBV6gW6emV3+w/EIcju6piK5gWQnyc6BgYECn05EilX01XjLrTPlrjVlZhvNrMPMOhob9ae9SKUUWuz9JFsAIPl8NLuURKQUCh162wLgCQAvJJ8/yLfjVF5w0tt2KbRA4Zw5c9zY6tWr3dhDDz2Ud15jedtNXbx40e0TynHu3LlurL6+3o1dunQptf3Eiasv//y/vXv3urHQDLvQ8Oa8eel/YVZXV7t9pvJrNCSfobe3APwGwE0kD5N8Erki/w7JvQD+JPlaRCaxcc/sZvaYE7on41xEpIT077IikVCxi0RCxS4SCRW7SCS019sEeDOKQsM4M2fOdGOhPdtCC0SGZmWdPXs2tT20SGUoxxkz/JfI6OioG/P2gQstshmKeUN5QHj2XSjmmcqv0RCd2UUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIbeMhAa3gkNXVVVVbmx0FBTSF1dXWp7KMdQHqHhtdBiJL29vanthexTB/gzDoHwrL3QsKhnOr5GAZ3ZRaKhYheJhIpdJBIqdpFIqNhFIlHWq/FmVvBV5sksdDU7NGkl1C+05lroGHpX42fNmuX2CV1xD109P3DggBvbt2/fhB8vtN5daJ0/b505wD8eIVP5NRoaSdCZXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIaCLMBHi5e2vTAeHhtVC/0HZNoWE5b+JHQ0OD28fbMgoADh8+7Ma++uorN7Z///7U9tA6c83NzW5syZIlBfULra/nmcqv0ZB8tn96g+RRkjvHtD1Psofk9uTj/tKmKSLFyudt/M8ArEtpf9nMViYfH2ablohkbdxiN7OtAPytN0VkSijmAt3TJLuSt/mN3p1IbiDZSbIztNiBiJRWocX+KoBlAFYC6APwondHM9toZh1m1tHY6P5OEJESK6jYzazfzEbN7BKA1wCsyTYtEclaQUNvJFvMrC/58kEAO0P3ny68obLQ8FooFhriCQ2HhdaT84blQn2Gh4fdWE9Pjxv74osv3Jg39BZaS2758uVubNWqVW4sNFvOm3VYyLZQU924xU7yLQB3A5hL8jCAnwC4m+RKAAagG8APSpeiiGRh3GI3s8dSml8vQS4iUkLxvZcRiZSKXSQSKnaRSKjYRSJR9llvoZlek503VBYaQgt9v6F+oZlthSxGGVpEcWhoyI0dP37cjYWG5fr7+1Pb29ra3D5NTU1u7Prrr3djoUUlve879HOZyq/REJ3ZRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEWYfeSE7p2UaFDL2FhI5FocfJG0Y7duyY2+fQoUNurLu7240dOXLEjV24cCG1PbTnXGtrqxsL7ecWGooMDSt6pvJrNDRsOHW/KxGZEBW7SCRU7CKRULGLRELFLhIJTYTJQKETYUJXfUPbFs2Y4f/Yzpw5M6F2APjyyy/dWOhqfGgrp9ra2tT20ISW+fPnu7HZs2cXlEdoLT/PdHyNAjqzi0RDxS4SCRW7SCRU7CKRULGLRELFLhKJfHaEaQPwcwDzkNsBZqOZvUKyCcBmAEuQ2xXmYTOLcpvW0NBbKBbaGsrbtggIDw2dPHkytf3rr792++zZs8eNhSa7hIYAvSG2lpYWt09okkzI6OioGytkDbrpKp8z+wiAH5pZO4DbATxFsh3AswA+NrMbAXycfC0ik9S4xW5mfWa2Lbk9CGA3gAUAHgCwKbnbJgDrS5SjiGRgQn+zk1wCYBWATwHMG7OT6xHk3uaLyCSVd7GTbADwLoBnzOyK/0+03B+mqX+cktxAspNk54kTJ4pKVkQKl1exk6xGrtDfNLP3kuZ+ki1JvAXA0bS+ZrbRzDrMrCO0CYCIlNa4xc7cZcvXAew2s5fGhLYAeCK5/QSAD7JPT0Syks+stzsAPA7gc5Lbk7bnALwA4B2STwI4AODhkmQ4iXjDaKGtlUJCQ1ehWV7e+m6AP0utq6vL7bNt2zY3dvRo6hs2AOHtmtrb21Pbb7rpJrdPQ0ODGwvNXhseHnZj3rBcaNhzuhq32M3sEwDeoOQ92aYjIqWi/6ATiYSKXSQSKnaRSKjYRSKhYheJRNkXnCx0q6TJwBtiCw29hb7fUL/Q0NDFixfdWE9PT2r7jh073D6hBSdDM9FWrFjhxtauXZvafuutt7p9QotserP5AGBwcNCNeUKLfU7l12iIzuwikVCxi0RCxS4SCRW7SCRU7CKRULGLRKKsQ29mVvAMscnAW6QwtHhhKBYa4hkZGXFjQ0NDbswbhgoNT4XyaG5udmPLly93Y97stnnz/AWNQnu2hfaqC80CrK6uTm0vdEh0sgt9Xzqzi0RCxS4SCRW7SCRU7CKRULGLRKLsE2Gmo0InVYSuIn/zzTdurK+vz42dOnUqtb22ttbt09ra6sZuuOEGN7Zs2bIJP2Z9fb3bx8sdCE/+CfF+Ntr+SUSmLRW7SCRU7CKRULGLRELFLhIJFbtIJMYdeiPZBuDnyG3JbAA2mtkrJJ8H8KcAjiV3fc7MPixVopNBIRNhvO2HgPDEj+PHj7sxb4snAOjv709tb2xsdPssXrzYjd18881urK2tzY15a9eF1tYrdEJRiDf0Fhouna7yGWcfAfBDM9tGchaAz0h+lMReNrN/KF16IpKVfPZ66wPQl9weJLkbwIJSJyYi2ZrQexmSSwCsAvBp0vQ0yS6Sb5D03yeKSMXlXewkGwC8C+AZMzsN4FUAywCsRO7M/6LTbwPJTpKdAwMDxWcsIgXJq9hJViNX6G+a2XsAYGb9ZjZqZpcAvAZgTVpfM9toZh1m1hG6SCQipTVusTN3GfR1ALvN7KUx7S1j7vYggJ3ZpyciWcnnavwdAB4H8DnJ7UnbcwAeI7kSueG4bgA/KEF+k4o3XFPoemahLY12797txvbu3evGvJl0CxcudPuE1pJbtGiRG2toaHBj58+fT20PDUWGhsNCs/ZCM+JmzEh/icc46y2fq/GfAEg7MtN6TF1kuonvPwtEIqViF4mEil0kEip2kUio2EUioQUnSyw0LOcNTwHAoUOH3NjBgwfdWFNTU2p7aIZae3u7G5szZ44bC/H+W9IbChtPaOgtNIzmzbILDfNN5e2fQnRmF4mEil0kEip2kUio2EUioWIXiYSKXSQSGnrLQGgmV2g/t8HBQTfW29vrxg4cOODGvGGo6upqt8/s2bPdWF1dnRsbGRlxY96wYmgILRQLLVQZGirzYjEuOBnfdywSKRW7SCRU7CKRULGLRELFLhIJFbtIJDT0NgHeME5o6G14eNiNnT171o0dO3asoNjcuXNT20OzzUILR4b6hb43bxHI0JBXTU2NGwsNr4UWnPSGIjX0JiLTlopdJBIqdpFIqNhFIqFiF4nEuFfjSdYC2ApgZnL/X5rZT0guBfA2gGYAnwF43Mz8WR+RKsU2Q6FJId5kktAV9/r6ejcWWkMvNAoxNDTkxjyFXo2PcSunQuRzZh8G8G0z+xZy2zOvI3k7gJ8CeNnM/gDAAIAnS5aliBRt3GK3nDPJl9XJhwH4NoBfJu2bAKwvRYIiko1892evSnZwPQrgIwBfAThpZpcnNB8GsKAkGYpIJvIqdjMbNbOVABYCWANgRb5PQHIDyU6Snd5a4iJSehO6Gm9mJwH8GsAfAZhD8vIFvoUAepw+G82sw8w6Ghsbi8lVRIowbrGTvI7knOR2HYDvANiNXNF/L7nbEwA+KFGOIpKBfCbCtADYRLIKuV8O75jZf5D8AsDbJP8GwP8CeL2EeU4K3uSJ0NppoeGp5uZmN7Z48eKCHnPRokWp7aFtnELr04WeKzTk5Q2VhSbPhGKhCTmhdfK8YcpQ7qHveSobt9jNrAvAqpT2/cj9/S4iU4D+g04kEip2kUio2EUioWIXiYSKXSQSDM1qyvzJyGMALu9dNBfA8bI9uU95XEl5XGmq5bHYzK5LC5S12K94YrLTzDoq8uTKQ3lEmIfexotEQsUuEolKFvvGCj73WMrjSsrjStMmj4r9zS4i5aW38SKRULGLRKIixU5yHckvSe4j+Wwlckjy6Cb5OcntJDvL+LxvkDxKcueYtiaSH5Hcm3wu+UofTh7Pk+xJjsl2kveXIY82kr8m+QXJXST/PGkv6zEJ5FHWY0KyluRvSe5I8virpH0pyU+TutlM0l+ON42ZlfUDQBVya9jdAKAGwA4A7eXOI8mlG8DcCjzvXQBWA9g5pu3vATyb3H4WwE8rlMfzAH5U5uPRAmB1cnsWgD0A2st9TAJ5lPWYACCAhuR2NYBPAdwO4B0Ajybt/wzgzybyuJU4s68BsM/M9ltunfm3ATxQgTwqxsy2AjhxVfMDyK3SC5RptV4nj7Izsz4z25bcHkRuJaQFKPMxCeRRVpaT+YrOlSj2BQAOjfm6kivTGoBfkfyM5IYK5XDZPDPrS24fATCvgrk8TbIreZtf1oUDSS5BbrGUT1HBY3JVHkCZj0kpVnSO/QLdWjNbDeA+AE+RvKvSCQG53+zI/SKqhFcBLENuQ5A+AC+W64lJNgB4F8AzZnZ6bKycxyQlj7IfEytiRWdPJYq9B0DbmK/dlWlLzcx6ks9HAbyPyi6z1U+yBQCSz0crkYSZ9ScvtEsAXkOZjgnJauQK7E0zey9pLvsxScujUsckee6TmOCKzp5KFPvvANyYXFmsAfAogC3lToJkPclZl28D+C6AneFeJbUFuVV6gQqu1nu5uBIPogzHhLnVH18HsNvMXhoTKusx8fIo9zEp2YrO5brCeNXVxvuRu9L5FYC/qFAONyA3ErADwK5y5gHgLeTeDl5E7m+vJ5HbIPNjAHsB/DeApgrl8a8APgfQhVyxtZQhj7XIvUXvArA9+bi/3MckkEdZjwmAW5FbsbkLuV8sfznmNftbAPsA/ALAzIk8rv5dViQSsV+gE4mGil0kEip2kUio2EUioWIXiYSKXSQSKnaRSPwffyIxG32ku9UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUkklEQVR4nO3df4xV5Z3H8fe3yK/KYBGmgIgDAroCRSATUiq2rm4t2jTaZNNWN103caXp1mSbdP8w3WTrbvaPdrNt081u3NDV1DZdqbaauhvaqtQGTS064AhYrfwoWJAfQwUB5YcD3/3jHrID3u8zM+f+muH5vJLJ3Hm+99zzzJn7nXPv+d7neczdEZHz3/ta3QERaQ4lu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULLLoJnZfDP7hZkdMLP3fFDDzC42s8fM7G0z22lmt7ein3I2JbuU8S7wMHBnEP8P4CQwGfgL4D4zm9ekvknA9Am684uZ7QD+HfhLoAP4OXCHux9vwL5mA1vc3fq0XQgcBOa7+2tF2w+A3e5+T737IAOnM/v56TPAcmAmsAD4q2p3MrNlZnYo8bWsxL6vAHrPJHrhJUBn9ha7oNUdkIb4N3d/A8DM/gdYWO1O7v4s8IE673sccPictreAtjrvRwZJZ/bz094+t9+hkoDNchQYf07beOBIE/sgVSjZM2Zm15rZ0cTXtSUe9jXgAjOb06ftauDl+vRaytLL+Iy5+zOUOOubmQGjgVHFz2MqD+cn3P1tM3sU+Ccz+2sqbyFuAT5St45LKTqzSxkdwDH+/2x9DPhdn/jfAGOB/cBDwBfdXWf2FlPpTSQTOrOLZELJLpIJJbtIJpTsIploault0qRJ3tHR0cxdimRl586dHDhwwKrFakp2M1sOfAcYAfyXu389df+Ojg5+/etf17LLIalSdq4uVe1Ixco+ZpnHSynbx2Y9Xn+PWe99DXUf+Uj8cYbSL+PNbASVoYw3AXOB28xsbtnHE5HGquU9+xJgq7tvd/eTwCoqn5QSkSGolmSfBvyhz8+7irazmNkKM+sys66enp4adicitWj41Xh3X+nune7e2d7e3ujdiUiglmTfDUzv8/OlRZuIDEG1JPsLwBwzm2lmo4DPAY/Xp1siUm+lS2/u3mtmdwO/oFJ6e0Ajm0SGrprq7O6+Glhdp76ISAPp47IimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimahpRRgz2wEcAU4Bve7eWY9OiUj91ZTshT919wN1eBwRaSC9jBfJRK3J7sATZrbezFZUu4OZrTCzLjPr6unpqXF3IlJWrcm+zN0XAzcBXzKzj557B3df6e6d7t7Z3t5e4+5EpKyakt3ddxff9wOPAUvq0SkRqb/SyW5mF5pZ25nbwI3A5np1TETqq5ar8ZOBx8zszOP8t7v/vC69EpG6K53s7r4duLqOfRGRBlLpTSQTSnaRTCjZRTKhZBfJRD0+G5+Nt956a1DtAO97X/z/dObMmWHs+PHjYWzEiBFh7MiRI1Xbd+3aFW5z+vTpMHby5Mkwtn379jB29OjRqu2p45E6jtu2bQtjF1wQP43nz59ftX3ZsmXhNrNmzQpjw5nO7CKZULKLZELJLpIJJbtIJpTsIplo+tV4d2/2Lutmx44dVdufe+65cJt9+/aFsblz55bqx6lTpwa9v9///vfhNqm/yeHDh8PY66+/HsaKMRPvceGFF4bbpK7Gd3d3h7He3t4wduutt1Ztv+KKK8JtLr/88jA2nOnMLpIJJbtIJpTsIplQsotkQskukgklu0gmNBBmEKLS0Pr168Ntfvazn4Wxt99+O4yNGTMmjKVKb9GgkPHjx5faV2pAzrhx48LYpZdeWrV9ypQp4TYTJkwIY1HZE9Ilu2hG44kTJ4bbnK90ZhfJhJJdJBNKdpFMKNlFMqFkF8mEkl0kE00vvaXmIBvqFi9eXLU9NZIrNbpq3bp1YSw12iy1v6iP119/fbjN7Nmzw1iqZJdalfeiiy6q2j569Ohwm9WrV4exX/3qV2Fs3rx5YeyGG26o2p6a/284P0dT+v2tzOwBM9tvZpv7tF1sZk+a2Zbie1wgFZEhYSD/wr4HLD+n7R5gjbvPAdYUP4vIENZvsrv7WuDNc5pvAR4sbj8I3FrfbolIvZV9czLZ3fcUt/dSWdG1KjNbYWZdZtZ14MCBkrsTkVrVfCXCK3MahfMauftKd+90985JkybVujsRKalssu8zs6kAxff99euSiDRC2dLb48AdwNeL7z8d6IbDecLJqJy0aNGicJtUiee2224LY2XLP1FpKzVCLVUOS42wmz59ehiLlqhKTRy5atWqMBYtJwWwdOnSMHbZZZeFschwfo6mDKT09hDwHHClme0yszupJPnHzWwL8GfFzyIyhPV7Znf36PRT/dMKIjIknZ8fFRKR91Cyi2RCyS6SCSW7SCa01tsgROWw1ISNqUkU29rawtjYsWPD2MmTJ8NYdHxHjhwZbpNy4sSJMJYafRdNpvnqq6+G22zYsCGMRWvHQbq8+cEPfrBq+6hRo8JtTp8+HcaGM53ZRTKhZBfJhJJdJBNKdpFMKNlFMqFkF8mE1nobhN7e3qrtqbJQKhaNDIPyJcrU/iKpUlNq9F2qj9u2bava/vzzz4fbHDp0KIwtXLgwjM2dOzeMRaXP1O+l0puIDGtKdpFMKNlFMqFkF8mEkl0kE02/Gl/mavFQEV2NT129TQ1AueCC+PC/++67YSx1DKOrzKk+puaZSw3IiQa7AKxfv75q+7PPPhtukzpWN910Uxj70Ic+FMai/qeOx3B+jqbozC6SCSW7SCaU7CKZULKLZELJLpIJJbtIJjQQZhCiElWqjJMqr6UGY6TKYan506LHTJXyUlKDdfbs2RPGotLb1q1bw21SSzUtWbIkjLW3t4exSFRGhfTvPJwNZPmnB8xsv5lt7tN2r5ntNrPu4uvmxnZTRGo1kJfx3wOWV2n/trsvLL5W17dbIlJv/Sa7u68F3mxCX0SkgWq5QHe3mW0sXuaHk6Ob2Qoz6zKzrp6enhp2JyK1KJvs9wGzgIXAHuCb0R3dfaW7d7p7Z5kLKSJSH6WS3d33ufspdz8NfBeIL5WKyJBQqvRmZlPd/Uzd5dPA5tT9zxdlllBKlddS5bDUdqlRWVEZMFXKSy1fdezYsTD29NNPh7ForrmJEyeG23z2s58NYwsWLAhjKWXmDTxf9ZvsZvYQcB0wycx2AV8DrjOzhYADO4AvNK6LIlIP/Sa7u99Wpfn+BvRFRBpIH5cVyYSSXSQTSnaRTCjZRTKhUW+DEI2GKjt6LTXyKjVaLlU2ivZXtpT3xz/+MYx1d3eHsR07dlRtv+SSS8JtUss4tbW1hbHUcYyWqEodj/NVfr+xSKaU7CKZULKLZELJLpIJJbtIJpTsIplQ6a0OUhNORqUfSJd/yq4DF0lNUnn06NEw9swzz4SxF154IYxFI+mWLl0abjNv3rwwlhpxePLkyTAWlRVzHPWmM7tIJpTsIplQsotkQskukgklu0gmdDV+EKKr7qkr7qlYI5aGKjNY54033ghjq1fH63+89tprYWz+/PlV26+77rpwm46OjjCWqnikrqxHv3dqm9TfbDjTmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTAxkRZjpwPeByVRWgFnp7t8xs4uBHwEzqKwK8xl3P9i4rrZeVJJJlWpSJZ5ULDWvWqqMFsUOHoz/NC+++GKpWKr/V1111aDaAUaPHh3GUstQpQbJRMfjfC2vpQzkzN4LfMXd5wIfBr5kZnOBe4A17j4HWFP8LCJDVL/J7u573H1DcfsI8AowDbgFeLC424PArQ3qo4jUwaDes5vZDGARsA6Y3Gcl171UXuaLyBA14GQ3s3HAT4Avu/vhvjGvvAGq+ibIzFaYWZeZdfX09NTUWREpb0DJbmYjqST6D9390aJ5n5lNLeJTgf3VtnX3le7e6e6d7e3t9eiziJTQb7Jb5ZLr/cAr7v6tPqHHgTuK23cAP61/90SkXgYy6u0a4PPAJjPrLtq+CnwdeNjM7gR2Ap9pSA+HgVQJKhqFBumRXO+8804YmzBhQhiLRsS9/PLL4TaPPPJIGHv11VfD2NVXXx3Grr322qrts2bNCrcpOy9cmRFsZUfRDWf9Jru7PwtEv/0N9e2OiDSKPkEnkgklu0gmlOwimVCyi2RCyS6SCU04OQipkWhlpMpyZZc7ikpsqfLamjVrwlhbW1sY+8QnPhHGrrnmmqrt48ePD7dJLUOVmpyzzFJZZUcqDmc6s4tkQskukgklu0gmlOwimVCyi2RCyS6SCZXeBiFVKisjtWZbavLFVAnw9ddfr9qemjjyyJEjYezGG28MYx/72MfC2CWXXFK1PTXaLKXs2nc5TiwZ0ZldJBNKdpFMKNlFMqFkF8mEkl0kE02/Gj+cr45GAy5SAydSV4pTsdRAmP37q07kC8CmTZuqtkdX6QFmzJgRxj71qU+FsdQcdKNGjaranppbL9oG0sf4xIkTYSz6m6Wu7g/n52iKzuwimVCyi2RCyS6SCSW7SCaU7CKZULKLZKLf0puZTQe+T2VJZgdWuvt3zOxe4C7gzNKsX3X31Y3q6PkoVf6J5k4D+M1vfhPGfvnLX1ZtT83vdvvtt4ex5cuXh7EpU6aEsagclhoIU2Yuuf5Ex7hsuXQ4G0idvRf4irtvMLM2YL2ZPVnEvu3u/9q47olIvQxkrbc9wJ7i9hEzewWY1uiOiUh9Deo9u5nNABYB64qmu81so5k9YGbx0qIi0nIDTnYzGwf8BPiyux8G7gNmAQupnPm/GWy3wsy6zKyrp6en2l1EpAkGlOxmNpJKov/Q3R8FcPd97n7K3U8D3wWWVNvW3Ve6e6e7d7a3t9er3yIySP0mu1UuW94PvOLu3+rTPrXP3T4NbK5/90SkXgZyNf4a4PPAJjPrLtq+CtxmZguplON2AF8YyA6H89I6Ufmn7Lxq73//+8PY3r17w9hTTz0Vxrq7u6u2X3nlleE2d911Vxi77LLLwtixY8fCWGTs2LFhLFVeS8VS8/VFz7fU32w4P0dTBnI1/lmg2m+vmrrIMKJP0IlkQskukgklu0gmlOwimVCyi2RCyz8NQlSSSS0LlZo4MjX54hNPPBHGNm+OP9IwbVr1YQuf/OQnw22uuuqqMFZ2dFg0eWTq8VLltdRxTI0ePHnyZNX21KSS9V7ma6jQmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTKj0NghlSm+pstCbb74ZxtauXRvGtmzZEsY6Ojqqtk+dOrVqO6TLUPUuUaVGm6WOVSpWdj29iEpvIjKsKdlFMqFkF8mEkl0kE0p2kUwo2UUy0fTSW6qUM9RF5Z9UOenw4cNhLDV6bcOGDWHs4MGDYWzBggVV21Olt9Ros2j0GqRLVGUm4Uyt9ZYqoaWeU6n+R8pOIDrU6cwukgklu0gmlOwimVCyi2RCyS6SiX6vxpvZGGAtMLq4/4/d/WtmNhNYBUwE1gOfd/fqE36d51JXilNX43fu3BnGUivetrW1hbFomafZs2eH25S9qp7arkzVpRGDXaLlplL9O3HiRBgbzgZyZj8BXO/uV1NZnnm5mX0Y+AbwbXefDRwE7mxYL0WkZv0mu1ccLX4cWXw5cD3w46L9QeDWRnRQROpjoOuzjyhWcN0PPAlsAw65e29xl11A9TmMRWRIGFCyu/spd18IXAosAf5koDswsxVm1mVmXan3oSLSWIO6Gu/uh4CngaXAB8zszAW+S4HdwTYr3b3T3Tvb29tr6auI1KDfZDezdjP7QHF7LPBx4BUqSf/nxd3uAH7aoD6KSB0MZCDMVOBBMxtB5Z/Dw+7+v2b2W2CVmf0z8CJw/0B2mFr+Z6iL+j5mzJhwm8mTJ4exOXPmhLEpU6aEsdSAkUWLFlVtnzlzZrhNyvHjx0v1IyrL9fb2Vm2HdAkttfxTqox29OjRqu1l5w0czvpNdnffCLznGeTu26m8fxeRYeD8/BcmIu+hZBfJhJJdJBNKdpFMKNlFMmHNnBPOzHqAM0O9JgEHmrbzmPpxNvXjbMOtHx3uXvXTa01N9rN2bNbl7p0t2bn6oX5k2A+9jBfJhJJdJBOtTPaVLdx3X+rH2dSPs503/WjZe3YRaS69jBfJhJJdJBMtSXYzW25mvzOzrWZ2Tyv6UPRjh5ltMrNuM+tq4n4fMLP9Zra5T9vFZvakmW0pvk9oUT/uNbPdxTHpNrObm9CP6Wb2tJn91sxeNrO/LdqbekwS/WjqMTGzMWb2vJm9VPTjH4v2mWa2rsibH5nZ4Bayc/emfgEjqMxhdzkwCngJmNvsfhR92QFMasF+PwosBjb3afsX4J7i9j3AN1rUj3uBv2vy8ZgKLC5utwGvAXObfUwS/WjqMQEMGFfcHgmsAz4MPAx8rmj/T+CLg3ncVpzZlwBb3X27V+aZXwXc0oJ+tIy7rwXePKf5Fiqz9EKTZusN+tF07r7H3TcUt49QmQlpGk0+Jol+NJVX1H1G51Yk+zTgD31+buXMtA48YWbrzWxFi/pwxmR331Pc3gvEU9w03t1mtrF4md/wtxN9mdkMKpOlrKOFx+ScfkCTj0kjZnTO/QLdMndfDNwEfMnMPtrqDkHlPzuVf0StcB8wi8qCIHuAbzZrx2Y2DvgJ8GV3P2spnWYekyr9aPox8RpmdI60Itl3A9P7/BzOTNto7r67+L4feIzWTrO1z8ymAhTf97eiE+6+r3iinQa+S5OOiZmNpJJgP3T3R4vmph+Tav1o1TEp9n2IQc7oHGlFsr8AzCmuLI4CPgc83uxOmNmFZtZ25jZwI7A5vVVDPU5lll5o4Wy9Z5Kr8GmacEysMpPn/cAr7v6tPqGmHpOoH80+Jg2b0blZVxjPudp4M5UrnduAv29RHy6nUgl4CXi5mf0AHqLycvBdKu+97qSyQOYaYAvwFHBxi/rxA2ATsJFKsk1tQj+WUXmJvhHoLr5ubvYxSfSjqccEWEBlxuaNVP6x/EOf5+zzwFbgEWD0YB5XH5cVyUTuF+hEsqFkF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQT/weYte4bQI3fgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASZUlEQVR4nO3df4xdZZ3H8feH2h9AK6XbsU5K6ZS2AmWBgYxdQQS2VoP800o2YmOUjWRrNpBogskSN1nZxGR1s2rchLgpP2I1rsiKRlzrQiEutEJopy5Oi0VoaRFK2xki/QHdpbT97h/3NE7rec7M3J/tPJ9XMpk7z/eee74c+plz5zz3nKOIwMzGvzM63YCZtYfDbpYJh90sEw67WSYcdrNMOOxmmXDYzTLhsNuYSfpzSY9Iel3Sn3xQQ9J/S/o/SW8WX7/rRJ92Iofd6vEO8CBwa8Vzbo+IqcXXhW3qyyo47OOMpJ2SvihpQNJ+ST+UNKWZ64iI30XEfcBzzXxday2HfXz6BHADMA+4DPjrsidJukbSvoqvaxro4Z+Kt/m/knR9A69jTfKuTjdgLfGvEfEagKSfAb1lT4qI9cD0Fqz/74DfAoeBTwI/k9QbEdtbsC4bJe/Zx6c9wx4fAqa2c+UR8UxEHIyItyNiNfAr4MZ29mB/ymHPmKQPDTtiXvb1oSatKgA16bWsTn4bn7GIWEcde31JAiYDk4qfp9ReLt6WNB34C+AJ4AhwM3At8PkmtW11ctitHnOBHcN+/l/gZaAHmAh8BbgIOAo8DyyPiBfa3KOdRL54hVke/De7WSYcdrNMOOxmmXDYzTLR1qPxM2fOjJ6ennau0iwrO3fu5PXXXy/9TENDYZd0A/AtYAJwb0R8ter5PT09bNy4sZFVmlmF97///cla3W/jJU0A7gY+BiwCVkhaVO/rmVlrNfI3+2JgW0S8FBGHgQeAZc1py8yarZGwzwZeGfbzq8XYCSStlNQvqX9oaKiB1ZlZI1p+ND4iVkVEX0T0dXV1tXp1ZpbQSNh3AXOG/XxeMWZmp6BGwr4RWChpnqRJ1C5S8HBz2jKzZqt76i0ijki6HXiE2tTb/RHha5KZnaIammePiDXAmib1YmYt5I/LmmXCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2WioTvCSNoJHASOAkcioq8ZTZlZ8zUU9sJfRsTrTXgdM2shv403y0SjYQ/gUUmbJK0se4KklZL6JfUPDQ01uDozq1ejYb8mIq4EPgbcJunak58QEasioi8i+rq6uhpcnZnVq6GwR8Su4vsg8BNgcTOaMrPmqzvsks6WNO34Y+CjwJZmNWZmzdXI0fhZwE8kHX+df4+I/2pKV2bWdHWHPSJeAi5vYi9m1kKeejPLhMNulgmH3SwTDrtZJprx2fhsHDx4sHT8wIEDyWXOOCP9+7S7u7uuPooZkFJvvPFG6fiuXbvqWlfqvxnglVdeGfNyx44dSy5TtR137NiRrFW95mWXXVY6vmTJkuQyCxcuTNZOZ96zm2XCYTfLhMNulgmH3SwTDrtZJnw0fgy2b99eOr5u3brkMnv27EnWFi1a1HBPJ3vttddKx1O9j2T//v3JWtUR8tQsxNSpU5PL7Nu3L1nbtGlTshYRydqyZctKxy+++OLkMj4ab2anNYfdLBMOu1kmHHazTDjsZplw2M0y4am3MRgcHCwdr5p6e/TRR5O1w4cPN9zTySZPnlw6XjXlNXHixGSt6kSeqtrcuXPHNA4wffr0ZG3btm3J2qFDh5K1888/v3R81qxZyWXGK+/ZzTLhsJtlwmE3y4TDbpYJh90sEw67WSY89TYGS5cuLR2fM2dOcplLL700WVu/fn2yduTIkdE3NszixeW320ud/QUwf/78utZVdbbZWWedVTpeNd34xBNPJGtV2+q9731vsnbdddeVjr/vfe9LLjNejbhnl3S/pEFJW4aNzZC0VtKLxfdzW9ummTVqNG/jvwPccNLYncDjEbEQeLz42cxOYSOGPSKeBP5w0vAyYHXxeDWwvLltmVmz1XuAblZE7C4e76F2R9dSklZK6pfUPzQ0VOfqzKxRDR+Nj9pRmuSRmohYFRF9EdHX1dXV6OrMrE71hn2vpG6A4nv5GSJmdsqod+rtYeAW4KvF9582raNTWOosrwULFiSXue2225K1z372s3X1Uc+U19lnn51cZtKkSXX1USW1rQYGBpLL3Hvvvcla1cUoly9fnqylptiqbqE1Xo1m6u0HwNPAhZJelXQrtZB/RNKLwNLiZzM7hY24Z4+IFYnSh5vci5m1kD8ua5YJh90sEw67WSYcdrNM+Ky3Jqi6YOOMGTPa2MmpI3URyOeffz65TH9/f7L2zjvvJGuXXHJJsvae97yndLxq6q1qavN05j27WSYcdrNMOOxmmXDYzTLhsJtlwmE3y4Sn3k5D9UwbtXuqafv27aXjVffFO3DgQLLW29ubrF1++eXJ2rvf/e5kLTfes5tlwmE3y4TDbpYJh90sEw67WSZ8NP40VM/R81Ycca+6RVXqpJbHHnssuczkyZOTtZtuuilZu+iii5K1M888s3R8vJ7sUsV7drNMOOxmmXDYzTLhsJtlwmE3y4TDbpYJT71Z3fbs2ZOsbdiwoXT8hRdeSC5zwQUXJGtLlixJ1qZNm5as5TjFljKa2z/dL2lQ0pZhY3dJ2iXp2eLrxta2aWaNGs3b+O8AN5SMfzMieouvNc1ty8yabcSwR8STwB/a0IuZtVAjB+hulzRQvM0/N/UkSSsl9UvqHxoaamB1ZtaIesP+bWA+0AvsBr6eemJErIqIvojo6+rqqnN1ZtaousIeEXsj4mhEHAPuARY3ty0za7a6pt4kdUfE7uLHjwNbqp5v49OaNenjsuvXry8d7+7uTi7zqU99KlmrusVT1dly9kcjhl3SD4DrgZmSXgW+DFwvqRcIYCfwuda1aGbNMGLYI2JFyfB9LejFzFrIH5c1y4TDbpYJh90sEw67WSZ81ptVnhlWdUumjRs3Jms7duwoHe/p6UkuU3UbpylTpiRrNjres5tlwmE3y4TDbpYJh90sEw67WSYcdrNMeOptnJFUOl41vVZ1z7Zf/OIXyVrqfm6Qvsfa4sXps6GvuOKKZO2MM7xfapS3oFkmHHazTDjsZplw2M0y4bCbZcJH48eZ1FH3o0ePJpcZHBxM1h566KFkrepWTpdeemnp+NKlS5PLzJ07N1mzxnnPbpYJh90sEw67WSYcdrNMOOxmmXDYzTIxmjvCzAG+C8yidgeYVRHxLUkzgB8CPdTuCvOJiHijda3acamTXSA99fbWW28ll3nqqaeSteeeey5ZqzqBZtGiRaXjvb29yWWstUazZz8C3BERi4APALdJWgTcCTweEQuBx4ufzewUNWLYI2J3RPy6eHwQ2ArMBpYBq4unrQaWt6hHM2uCMf3NLqkHuAJ4Bpg17E6ue6i9zTezU9Sowy5pKvAQ8IWIOOFi4lH7Q7H0j0VJKyX1S+ofGhpqqFkzq9+owi5pIrWgfz8iflwM75XUXdS7gdIPWEfEqojoi4i+rq6uZvRsZnUYMeyqHfq9D9gaEd8YVnoYuKV4fAvw0+a3Z2bNMpqz3j4IfBrYLOnZYuxLwFeBByXdCrwMfKIlHdqfqLqeXGpa7ve//31ymbvvvjtZS93GCWD+/PnJ2nXXXVc6fuGFFyaXsdYaMewRsR5ITex+uLntmFmr+BN0Zplw2M0y4bCbZcJhN8uEw26WCV9wcgyqzjZLqZomq3ddVa85MDBQOn7PPfckl9m0aVOy9vbbbydrK1asSNauv/760vEJEyYkl2n3tsqN9+xmmXDYzTLhsJtlwmE3y4TDbpYJh90sE556G4N2TuPUc2YbpM9ue/rpp5PLHDp0KFm7+uqrk7XU9BrAeeedVzreim3o6bXR8Z7dLBMOu1kmHHazTDjsZplw2M0y4aPxLdaKkzT27t2brG3evLl0/OWXX04ukzpyDvCZz3wmWbv44ouTtXe9q/yflk926Rzv2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmRpx6kzQH+C61WzIHsCoiviXpLuBvgOO3Zv1SRKxpVaOnq3qnhY4dO5asrVu3Lllbs6b8f0HVteSWL1+erN10003JWtWNOlP/3fVOoXl6rXGjmWc/AtwREb+WNA3YJGltUftmRPxL69ozs2YZzb3edgO7i8cHJW0FZre6MTNrrjH9zS6pB7gCeKYYul3SgKT7JZ3b7ObMrHlGHXZJU4GHgC9ExAHg28B8oJfanv/rieVWSuqX1D80NFT2FDNrg1GFXdJEakH/fkT8GCAi9kbE0Yg4BtwDLC5bNiJWRURfRPRVHdAxs9YaMeyqHT69D9gaEd8YNt497GkfB7Y0vz0za5bRHI3/IPBpYLOkZ4uxLwErJPVSm47bCXyuBf2Na1XTSfv370/WHnnkkWQtdSunSy65JLnMHXfckayde64PxYwXozkavx4omxz1nLrZacSfoDPLhMNulgmH3SwTDrtZJhx2s0z4gpMtVnWW19GjR5O1tWvXJmsbN25M1s4///zS8Ztvvjm5zPz585O11IUjob5bVLX77LVTpY9TgffsZplw2M0y4bCbZcJhN8uEw26WCYfdLBOeemuxqgtHvvnmm8naz3/+82Rt+/btydqCBQtKx6vOXqt3eq1Ksy84Wa8cp9hSvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU2xhUTRulvPXWW8naU089laxt2LAhWauasjvnnHNKx2fNmpVcpt3TYc12uvffLt6zm2XCYTfLhMNulgmH3SwTDrtZJkY8Gi9pCvAkMLl4/o8i4suS5gEPAH8GbAI+HRGHW9lsp9VzZPfw4fQm2bp1a7K2b9++ZG3q1KnJ2sKFC0vHq27/1Ioj1u289puPuI/OaPbsbwNLIuJyardnvkHSB4CvAd+MiAXAG8CtLevSzBo2Ytij5vjE7sTiK4AlwI+K8dXA8lY0aGbNMdr7s08o7uA6CKwFtgP7IuJI8ZRXgdkt6dDMmmJUYY+IoxHRC5wHLAYuGu0KJK2U1C+pf2hoqL4uzaxhYzoaHxH7gF8CVwHTJR0/wHcesCuxzKqI6IuIvq6urkZ6NbMGjBh2SV2SphePzwQ+AmylFvq/Kp52C/DTFvVoZk0wmhNhuoHVkiZQ++XwYET8p6TfAg9I+grwP8B9LezztHXWWWcla1dffXWyNnPmzGRt9uz04ZGrrrqqdHzevHnJZeqdumr2CSitOKHFt3/6oxHDHhEDwBUl4y9R+/vdzE4D/gSdWSYcdrNMOOxmmXDYzTLhsJtlQu2cgpA0BLxc/DgTeL1tK09zHydyHyc63fqYGxGln15ra9hPWLHUHxF9HVm5+3AfGfbht/FmmXDYzTLRybCv6uC6h3MfJ3IfJxo3fXTsb3Yzay+/jTfLhMNulomOhF3SDZJ+J2mbpDs70UPRx05JmyU9K6m/jeu9X9KgpC3DxmZIWivpxeL7uR3q4y5Ju4pt8qykG9vQxxxJv5T0W0nPSfp8Md7WbVLRR1u3iaQpkjZI+k3Rxz8W4/MkPVPk5oeSJo3phSOirV/ABGrXsLsAmAT8BljU7j6KXnYCMzuw3muBK4Etw8b+GbizeHwn8LUO9XEX8MU2b49u4Mri8TTgBWBRu7dJRR9t3SaAgKnF44nAM8AHgAeBTxbj/wb87VhetxN79sXAtoh4KWrXmX8AWNaBPjomIp4E/nDS8DJqV+mFNl2tN9FH20XE7oj4dfH4ILUrIc2mzdukoo+2ipqmX9G5E2GfDbwy7OdOXpk2gEclbZK0skM9HDcrInYXj/cA6Xsst97tkgaKt/kt/3NiOEk91C6W8gwd3CYn9QFt3iatuKJz7gforomIK4GPAbdJurbTDUHtNzu1X0Sd8G1gPrUbguwGvt6uFUuaCjwEfCEiDgyvtXOblPTR9m0SDVzROaUTYd8FzBn2c/LKtK0WEbuK74PAT+jsZbb2SuoGKL4PdqKJiNhb/EM7BtxDm7aJpInUAvb9iPhxMdz2bVLWR6e2SbHufYzxis4pnQj7RmBhcWRxEvBJ4OF2NyHpbEnTjj8GPgpsqV6qpR6mdpVe6ODVeo+Hq/Bx2rBNVLsq5H3A1oj4xrBSW7dJqo92b5OWXdG5XUcYTzraeCO1I53bgb/vUA8XUJsJ+A3wXDv7AH5A7e3gO9T+9rqV2g0yHwdeBB4DZnSoj+8Bm4EBamHrbkMf11B7iz4APFt83djubVLRR1u3CXAZtSs2D1D7xfIPw/7NbgC2Af8BTB7L6/rjsmaZyP0AnVk2HHazTDjsZplw2M0y4bCbZcJhN8uEw26Wif8HhJAtV2Fb+ucAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 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": 13,
   "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\n",
    "theta_size = N * block_len * cir_depth # The size of the network parameter theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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(theta):\n",
    "\n",
    "    # Initialize the network with UAnsatz\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # Build a hierarchy:\n",
    "    for layer_num in range(cir_depth):\n",
    "        \n",
    "        for which_qubit in range(N):\n",
    "            cir.ry(theta[block_len * layer_num * N + which_qubit], which_qubit)\n",
    "\n",
    "        for which_qubit in range(1, N):\n",
    "            cir.cnot([which_qubit - 1, which_qubit])\n",
    "\n",
    "    return cir.U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:16.619406Z",
     "start_time": "2021-03-09T03:47:16.600512Z"
    }
   },
   "outputs": [],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    \n",
    "    # Initialize the list of learnable parameters, and fill the initial value with the uniform distribution of [0, 2*pi]\n",
    "    def __init__(self, shape, dtype='float64'):\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # Create the parameter theta for learning U\n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # Create a parameter phi to learn V_dagger\n",
    "        self.phi = self.create_parameter(shape=shape, \n",
    "                                         default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                         dtype=dtype, is_bias=False)\n",
    "        \n",
    "        # Convert Numpy array to Tensor supported in Paddle\n",
    "        self.M = paddle.to_tensor(Mat_generator())\n",
    "        self.weight = paddle.to_tensor(weight)\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 = U_theta(self.theta)\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        \n",
    "        V = U_theta(self.phi)\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(matmul(U_dagger,matmul(self.M, V)))[i][i]\n",
    "            singular_values[i] = paddle.real(matmul(U_dagger, matmul(self.M, V)))[i][i].numpy()\n",
    "        \n",
    "        # 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": 16,
   "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: 1537.4504\n",
      "iter: 10 loss: -90754.6659\n",
      "iter: 20 loss: -118936.4569\n",
      "iter: 30 loss: -135978.2444\n",
      "iter: 40 loss: -142946.4453\n",
      "iter: 50 loss: -147064.3601\n",
      "iter: 60 loss: -149600.0237\n",
      "iter: 70 loss: -151044.1873\n",
      "iter: 80 loss: -152079.2783\n",
      "iter: 90 loss: -152850.4454\n",
      "iter: 100 loss: -153473.4505\n",
      "iter: 110 loss: -154018.4474\n",
      "iter: 120 loss: -154517.7508\n",
      "iter: 130 loss: -154952.2068\n",
      "iter: 140 loss: -155299.7216\n",
      "iter: 150 loss: -155558.2483\n",
      "iter: 160 loss: -155750.0483\n",
      "iter: 170 loss: -155900.1620\n",
      "iter: 180 loss: -156023.2654\n",
      "iter: 190 loss: -156125.9570\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x137940510>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX20lEQVR4nO2db2yc1ZXGn5OQEBM72MaO4/xp/pBAiYAG5Eagoqq0asVWlQBphYoQ4gNqqm2RtlL3A2KlhZX2Q7taqPph1VVYUOm2W8r2j4oqulsWVYIINcUQcEIImxCSJo6J85c4CU1IcvbDvJEc9J5n7Dsz74Te5ydFGd/j+7537tzHM76Pz7nm7hBC/OUzo90DEEJUg8QuRCZI7EJkgsQuRCZI7EJkgsQuRCZc0khnM7sNwPcBzATw7+7+Hfb93d3dvnDhwpT7lLYz2zDqUy+WYkWy682YEf88PXv2bBg7d+5cGEsZ48yZM8MYGyO7FxtjdD92vTNnzjT1XkD83FLXDqPZ12R9ovnYt28fjh49WtoxWexmNhPAvwL4IoC9AF4xs2fdfWvUZ+HChfjxj3887XvNnj27tP3Pf/5z2GfOnDlhjC1utqgiovHVix0/fjyMTUxMhLEPP/xwagObRG9vbxhjc8Xuxea/s7OztJ0J+uDBg0n3uvzyy8PYpZdeWtreCrGzH96zZs0KY9GaY31OnTpV2n7PPfeEfRr5GL8WwA533+nupwE8DeD2Bq4nhGghjYh9EYA9k77eW7QJIS5CWr5BZ2brzGzYzIaPHDnS6tsJIQIaEfsogCWTvl5ctF2Au6939yF3H+rp6WngdkKIRmhE7K8AWGVmy81sNoCvAni2OcMSQjSb5N14dz9jZg8A+B/UrLcn3f3Nev2inc6UXfBUUm2taOxs9/aSS+IpZjv1l112WRg7ffp0GIt2z1kftovMds9T5orZZGw+2L1Yv+h+7Dm3IhOUXTOKNdt+bchnd/fnADzXyDWEENWgv6ATIhMkdiEyQWIXIhMkdiEyQWIXIhMa2o1PIbITmM0QxZgtxEjN5IoSaFItEpaQw5JTUrK8mAXIxpFKlIDCLECW7MISg5j11tHRUdrebFur3jWbfb8Uq1rv7EJkgsQuRCZI7EJkgsQuRCZI7EJkQqW78e4e7iKyHeGUpIrU0lMpO9MsqYKVdUqNMaKdXfa8WOzEiRNhjI0xchOYyzB37twwdvLkyTCWUruOrQHmXKQktNQjZX2nlM7SO7sQmSCxC5EJErsQmSCxC5EJErsQmSCxC5EJlSfCRKQcycSsCWYLsX7MkolO6GAndzAbJzrVA+BJISyZJJorNo4oaaVejFl2UeIKS2hhNh+ryTdv3rwwFo2RzWGqNZt65FgEs3RlvQkhQiR2ITJBYhciEyR2ITJBYhciEyR2ITKhIevNzHYBmABwFsAZdx+q8/1JRyhFdgfLMmJ2WCqRJZOaocaytVLr00VjSa1Bl1pXrbu7u7SdWXnMbjxw4EAYY+sgsuVYHzZXbJ2y+Uix5dgajmw5ehRZGJk6t7r7wSZcRwjRQvQxXohMaFTsDuB3Zvaqma1rxoCEEK2h0Y/xt7j7qJnNB/C8mW1z9xcnf0PxQ2AdACxYsKDB2wkhUmnond3dR4v/xwH8CsDaku9Z7+5D7j7U09PTyO2EEA2QLHYzm2tmXecfA/gSgC3NGpgQork08jF+AMCviq3+SwD8p7v/N+vg7qHlkXI8DrNPUm0tlmkUZcu1IhOKzQd7blE2F8v0Y1Zes+1NZlOm3ivFKmPzm7IG6l2TrZHoebNxRPeilmIYqYO77wTwqdT+QohqkfUmRCZI7EJkgsQuRCZI7EJkgsQuRCZUXnCS2QkRkdXEsqRSi1EyOywqlnj06NGwD7PJmOXV398fxljxxahoIyv0+Kc//SmMsed26NChMBY9N/aabd++PYwdPnw4jF199dVh7Prrry9tHxwcDPt0dXWFMWavsXWVYi0zUvronV2ITJDYhcgEiV2ITJDYhcgEiV2ITKh0N97dwx3LlKQWtvvJjhJKTXTYt29fafu2bdvCPrt27QpjfX19YWxgYCCMsaSKqFbb6Oho2IfF3n///TDGdtYXLVpU2s4SWnbu3BnG2E79pk2bwtjJkydL22+55ZawD6uTx9wVtq7YaxZdk+3gs+SfcAzT7iGE+FgisQuRCRK7EJkgsQuRCRK7EJkgsQuRCZVbb1GNNGYlRNYE68MstGgMALdIoqQQlkiycePGpHsxy+748eNhbPbs2aXtK1asCPuweWS23AcffDDtaw4NxSeEseQUlnQzNjYWxsbHx0vbI0sO4GsntR9bc82ubRihd3YhMkFiFyITJHYhMkFiFyITJHYhMkFiFyIT6lpvZvYkgK8AGHf3a4u2XgA/A7AMwC4Ad7n7kXrXmjFjBjo7O0tjrB5bZE2wumq9vb1h7NSpU2GM2VDLly8vbWcWCctsY/XdLr/88jDGnneUbXb33XeHfaI6bfXuxerCReNnWWMbNmwIY6wm36pVq8LYjTfeWNq+ePHisA+r8ZeSbQbwDLZo7bNaeFHGIVuLU3ln/yGA2z7S9iCAF9x9FYAXiq+FEBcxdcVenLf+0R/htwN4qnj8FIA7mjssIUSzSf2dfcDdz//Z0nuonegqhLiIaXiDzmslZsIyM2a2zsyGzWyY/Y4qhGgtqWLfb2aDAFD8X/4HyADcfb27D7n7UHd3d+LthBCNkir2ZwHcVzy+D8CvmzMcIUSrmIr19lMAnwPQZ2Z7ATwM4DsAnjGz+wHsBnBXowNhWUFRAUBW8JBZRsw+YePo6ekpbWfZTpEVBvDCht/4xjfCGCMqEMmypFjhzvnz54cxVhQzKiD62muvhX2Y9cZea3b8U5RJx9YAy2xjFjHrx17ryApOyaJjhVvrit3dI4P2C/X6CiEuHvQXdEJkgsQuRCZI7EJkgsQuRCZI7EJkQqUFJ80stBOYZRDZRiyTKOV69Ygyitj5ZcziiTIAAWDp0qXTHgcQn/XGzmxj89jR0RHGGFF227Fjx8I+0Vl69foxovlntlYr1g67X3RGHBtHCnpnFyITJHYhMkFiFyITJHYhMkFiFyITJHYhMqFy6y3KGmJFICNriNlCzE5Kteyic9SYrcJizJZjWV4nTpwIY1EGW2TvANw6ZDE2jjfeeKO0/aWXXgr7sKyxG264IYytWbMmjC1YsKC0nRX0ZGfYpawPgK+5yM5jayCyX9n49M4uRCZI7EJkgsQuRCZI7EJkgsQuRCZUuhs/Y8aMsBYXq4MW7TCyY3rYjipLZmC7ptHOOtuFZbvZ7DkfOnQojO3fvz+MRXPCauuxxBo2H6zOX1Rr7g9/+EPYhzkG7IgqVoMuOgZs7ty5YR+2C87WDnutU3bjmUMVvc5sfHpnFyITJHYhMkFiFyITJHYhMkFiFyITJHYhMmEqxz89CeArAMbd/dqi7REAXwNwvuDZQ+7+XL1rnTt3jtoJEZElw47UYbBkgdQkmRTYXLATb48fPx7Gorp2fX19YR9mYTIbavfu3WFs69atpe3Mrlu+fHkYW7lyZRjr7+8PY5Etyl7nVCsy9YityNJlVmS0FhtNhPkhgNtK2r/n7muKf3WFLoRoL3XF7u4vAjhcwViEEC2kkd/ZHzCzETN70szKjzcVQlw0pIr9BwCuBLAGwBiAR6NvNLN1ZjZsZsNHjhxJvJ0QolGSxO7u+939rLufA/A4gLXke9e7+5C7D0XnmwshWk+S2M1s8gn3dwLY0pzhCCFaxVSst58C+ByAPjPbC+BhAJ8zszUAHMAuAF+f6g0jayDlOB5mvTGrg2YGEdslOtKIWSTMqmG2FrPXmL0SzQmz19hcbd++PYxt2LAhjI2MjJS2s/m47rrrwthVV10Vxrq6usJYytFKLIsxWgMArynInneULcfWYgp1xe7ud5c0P9HUUQghWo7+gk6ITJDYhcgEiV2ITJDYhcgEiV2ITKi04KS7hzYPsy1S7LrUzKWUjLjUTChm2TX7uCY2V8wCHB8fD2N79uyZdr/u7u6wz8KFC8MY68fWTopdyuaKWW+pr3XUj10vBb2zC5EJErsQmSCxC5EJErsQmSCxC5EJErsQmVCp9QbEdgKzyiIrhFkkLHMpNSMusl3Y9ZhVwywvZtWw5xb1O3nyZNjnvffeC2M7duxIikXFNK+44oqwzyc+8YkwFhXSBPj8R+Ng643Nb0rBVICv1QhmKUbj11lvQgiJXYhckNiFyASJXYhMkNiFyITKE2Gi3Wm28xjtMJ4+fTrsk1q/K2U3no2D7cKyGNu1Zkkh0U4yO07q3XffDWPvvPNOGHv//ffDWHSU08033xz2WbVqVRhjrwur1/fBBx+UtrO5Z4lGLOkpNUkm6sfWcEptPb2zC5EJErsQmSCxC5EJErsQmSCxC5EJErsQmTCV45+WAPgRgAHUjnta7+7fN7NeAD8DsAy1I6Ducve6x7QyeyIisiBYAgRLZmA2CItF9ztx4kTYh1k1HR0dYYzZaywWWWyslhyz3nbv3h3GIlsLAD796U+Xtt96661hn6VLl4YxlqzDnltksbG5nzdvXhhjsOQlZqOlHG+WpKMpfM8ZAN9299UAbgLwTTNbDeBBAC+4+yoALxRfCyEuUuqK3d3H3P214vEEgLcALAJwO4Cnim97CsAdLRqjEKIJTOt3djNbBuAGABsBDLj7WBF6D7WP+UKIi5Qpi93MOgH8AsC33P3Y5JjXfuko/cXDzNaZ2bCZDR85UvdXeiFEi5iS2M1sFmpC/4m7/7Jo3m9mg0V8EEDpLom7r3f3IXcf6unpacaYhRAJ1BW71bannwDwlrs/Nin0LID7isf3Afh184cnhGgWU8l6+wyAewFsNrPXi7aHAHwHwDNmdj+A3QDuqnchMwstsZQjlJjVweq7pdgWQGzjMJuPwSy0JUuWhLE5c+aEsV27dpW2v/zyy2GfV155JYwdOnQojC1evDiMXXPNNaXtLJuPWZipGY7Ra8b6MEs39VgxlmUX1bVjfdi9IuqK3d03AIjM5y9M+45CiLagv6ATIhMkdiEyQWIXIhMkdiEyQWIXIhMqLTg5Y8aM8BgfZpVFNgMrUsmO6WGZbczOiwoRdnV1hX2YRXLZZZcljSOy1wBg06ZNpe2//e1vwz7bt28PY/39/WFs9erVYeyqq64qbWd2ErPXBgbiv8aeO3duGIvWFSsqmWrzsfXInnd0NBcbB7tXhN7ZhcgEiV2ITJDYhcgEiV2ITJDYhcgEiV2ITKjUegNiKyoli4dlr7HMJWaDMMsrNbstgtk/zE46duxYGIued2TvAPw5s8w8djbblVdeWdrOahocPnw4jLFz5VLO9UvNXmMZh6wfO4Mt6scsYna9CL2zC5EJErsQmSCxC5EJErsQmSCxC5EJle/GR7ugKXXhUup6AenH6kRHBrEdWgbbUR0dHQ1jW7ZsCWN79+4tbWc7u8uWLQtj0TFOALB27dow1tfXV9rOnAR2nNTg4GAYY/MfrQO2cz4xMRHGmIPC5jjFbWJrOHKG2JrSO7sQmSCxC5EJErsQmSCxC5EJErsQmSCxC5EJda03M1sC4EeoHcnsANa7+/fN7BEAXwNwoPjWh9z9uXrXiywIZqNFlgazOpi9xuwTliARxZitknoM1b59+8LY22+/HcYOHjxY2r5w4cKwz4IFC8JYdIwTwI9/mj9/fmk7q9fHEnJYIgyz7KLXmq03lvCUUvsN4JZYZB02WxNTGfkZAN9299fMrAvAq2b2fBH7nrv/yxSuIYRoM1M5620MwFjxeMLM3gKwqNUDE0I0l2n9zm5mywDcAGBj0fSAmY2Y2ZNmpsPXhbiImbLYzawTwC8AfMvdjwH4AYArAaxB7Z3/0aDfOjMbNrNhVpxACNFapiR2M5uFmtB/4u6/BAB33+/uZ939HIDHAZT+obS7r3f3IXcf6u3tbda4hRDTpK7Yrba99wSAt9z9sUntkzMT7gQQZ2cIIdrOVHbjPwPgXgCbzez1ou0hAHeb2RrU7LhdAL4+lRtGdgKzvFJsBmZbMPuEXTMFZicdPXo0jG3bti2Mbd68OYxFteaY9fbJT34yjK1cuTKMsTp5kR124MCB0naA1w1kdhh7raMYs1+ZJUqzyhLWMBCvR9anJdabu28AUHaFup66EOLiQX9BJ0QmSOxCZILELkQmSOxCZILELkQmVFpw0sxCa4DZJ5HdwewYZp+k3AuIbY1UK48VX9y6dWsYGxkZCWNRocdrr7027DM0NBTGVqxYEcZYBtvu3btL28fGxsI+zMpj92KvZ2SHsdeMWYDs9WR23unTp8NYtI7ZGk468mraPYQQH0skdiEyQWIXIhMkdiEyQWIXIhMkdiEyoVLr7ezZszhx4kRpjFleUYzZD+x6zNJgWWrRNZlVc/z48TDGinns2bMnjI2Pj4exefPmlbb39/eHfVjhSHaOWpRhB8RZb8xCi2xDgFterKhntEbYa8bWTqotx85tSzkrMOXsOL2zC5EJErsQmSCxC5EJErsQmSCxC5EJErsQmVCp9QbE1lZKFg8jpVhfvVhkdzC7jtlTLMYy4iL7EogtHpbllVp8kdlQKeNg5/Ol2HxA/NzYc04lxT5mME3IehNChEjsQmSCxC5EJkjsQmSCxC5EJtTdjTezOQBeBHBp8f0/d/eHzWw5gKcBXAHgVQD3untcaAu1HcloFzGl9hvbBW/2MU7sfmzsbPc5NVmH7dJGiSbs+KeBgYGkcRw8eDCMRTXX2PUYrN4gm+OI1J3z1ESYlPuxHfcUHU3lnf0UgM+7+6dQO575NjO7CcB3AXzP3VcCOALg/ilcSwjRJuqK3Wucz9OcVfxzAJ8H8POi/SkAd7RigEKI5jDV89lnFie4jgN4HsA7AI66+/nPZHsBLGrJCIUQTWFKYnf3s+6+BsBiAGsBxGf8fgQzW2dmw2Y2zI4oFkK0lmntxrv7UQC/B3AzgG4zO78zshjAaNBnvbsPuftQd3d3A0MVQjRCXbGbWb+ZdRePOwB8EcBbqIn+r4tvuw/Ar1s0RiFEE5iKZzEI4Ckzm4naD4dn3P03ZrYVwNNm9k8ANgF4ot6FzCzJJkk5Mir1iKcUmBXW2dkZxlg9NlYzjiWFRPXkmL3GPnFNTEyEMWajsdp7ESlHJAFAR0fHtO/F1kdqjUK2tlOPI5suzP6rqzx3HwFwQ0n7TtR+fxdCfAzQX9AJkQkSuxCZILELkQkSuxCZILELkQnWbBuK3szsAIDdxZd9AOK0qerQOC5E47iQj9s4lrp7qW9bqdgvuLHZsLsPteXmGofGkeE49DFeiEyQ2IXIhHaKfX0b7z0ZjeNCNI4L+YsZR9t+ZxdCVIs+xguRCW0Ru5ndZmZvm9kOM3uwHWMoxrHLzDab2etmNlzhfZ80s3Ez2zKprdfMnjez7cX/PW0axyNmNlrMyetm9uUKxrHEzH5vZlvN7E0z+9uivdI5IeOodE7MbI6Z/dHM3ijG8Y9F+3Iz21jo5mdmFqcCluHulf4DMBO1slYrAMwG8AaA1VWPoxjLLgB9bbjvZwHcCGDLpLZ/BvBg8fhBAN9t0zgeAfB3Fc/HIIAbi8ddAP4PwOqq54SMo9I5AWAAOovHswBsBHATgGcAfLVo/zcAfzOd67bjnX0tgB3uvtNrpaefBnB7G8bRNtz9RQCHP9J8O2qFO4GKCngG46gcdx9z99eKxxOoFUdZhIrnhIyjUrxG04u8tkPsiwDsmfR1O4tVOoDfmdmrZrauTWM4z4C7jxWP3wMQV5toPQ+Y2UjxMb/lv05MxsyWoVY/YSPaOCcfGQdQ8Zy0oshr7ht0t7j7jQD+CsA3zeyz7R4QUPvJjtoPonbwAwBXonZGwBiAR6u6sZl1AvgFgG+5+wVnVlc5JyXjqHxOvIEirxHtEPsogCWTvg6LVbYadx8t/h8H8Cu0t/LOfjMbBIDi//F2DMLd9xcL7RyAx1HRnJjZLNQE9hN3/2XRXPmclI2jXXNS3PsoplnkNaIdYn8FwKpiZ3E2gK8CeLbqQZjZXDPrOv8YwJcAbOG9WsqzqBXuBNpYwPO8uAruRAVzYrXCaU8AeMvdH5sUqnROonFUPSctK/Ja1Q7jR3Ybv4zaTuc7AP6+TWNYgZoT8AaAN6scB4CfovZx8EPUfve6H7Uz814AsB3A/wLobdM4/gPAZgAjqIltsIJx3ILaR/QRAK8X/75c9ZyQcVQ6JwCuR62I6whqP1j+YdKa/SOAHQD+C8Cl07mu/oJOiEzIfYNOiGyQ2IXIBIldiEyQ2IXIBIldiEyQ2IXIBIldiEyQ2IXIhP8HPAMU92S1dkUAAAAASUVORK5CYII=\n",
      "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([theta_size])\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. [arXiv:2006.02336 (2020).](https://arxiv.org/abs/2006.02336)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}