Schmidt_decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simple_exact_diagonalization_routines.local_matrix_class import *\n",
    "from scipy.linalg import expm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from qutip import *\n",
    "from class_ising_model import *\n",
    "from class_double_bracket_flow import *\n",
    "\n",
    "import itertools\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I. Build up\n",
    "$\\rho = e^{-\\beta H_\\rho}$, $\\sigma = e^{-\\beta H_\\sigma}$ normalize by trace 1\n",
    "\n",
    "where $\\beta = \\frac{1}{k_B T}$, $H_\\rho = X + Z$, $H_\\sigma = Z$\n",
    "\n",
    "Swap: $ S (\\text{op1}\\otimes \\text{op2}) S = \\text{ op2}\\otimes \\text{op1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up states\n",
    "beta = 1\n",
    "X = np.array([[0,1],[1,0]])\n",
    "Y = np.array([[0,-1j],[1j,0]])\n",
    "Z = np.array([[1,0],[0,-1]])\n",
    "I = np.array([[1,0],[0,1]])\n",
    "H_rho = X + Z\n",
    "H_sigma = Z\n",
    "rho = expm(-beta*H_rho)/np.trace(expm(-beta*H_rho))\n",
    "sigma = expm(-beta*H_sigma)/np.trace(expm(-beta*H_sigma))\n",
    "Swap = .5*(np.kron(I, I) + np.kron(X, X) + np.kron(Y, Y) + np.kron(Z, Z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# Check swap\n",
    "print(np.all(Swap @ np.kron(X,Y) @ Swap == np.kron(Y,X)))\n",
    "\n",
    "# Check hamiltonian\n",
    "print(np.all(sigma == np.conjugate(sigma).T))\n",
    "print(np.all(rho == np.conjugate(rho).T))\n",
    "\n",
    "# Check trace 1\n",
    "print(np.trace(rho)==np.trace(sigma)==1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DME Channel\n",
    "The unitary evolution: $E_s^{\\rho}= e^{-is\\rho}$\n",
    "\n",
    "DME channel: $\\mathcal{E}_s^{\\rho}=Tr_1[e^{-is\\hat S}\\rho\\otimes\\sigma e^{is\\hat S}]$\n",
    "\n",
    "DME evolution approximation: $(\\mathcal{E}_{s/M}^{\\rho})^M = E_s^{\\rho}\\sigma(E_s^{\\rho})^\\dag + O(\\frac{s^2}{M})$\n",
    "\n",
    "Norm comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def partial_trace(dm, sub_index):\n",
    "    \"\"\"\n",
    "    Calculate the partial trace over the second qubit.\n",
    "\n",
    "    Parameters:\n",
    "    rho : numpy.ndarray\n",
    "        The 4x4 density matrix of a two-qubit system.\n",
    "\n",
    "    Returns:\n",
    "    numpy.ndarray\n",
    "        The 2x2 reduced density matrix of the first qubit.\n",
    "    \"\"\"\n",
    "    # Reshape rho from a 4x4 matrix into a 2x2x2x2 tensor\n",
    "    reshaped_dm = dm.reshape(2, 2, 2, 2)\n",
    "    # Perform the partial trace by summing out the inner 2x2 matrix\n",
    "    if sub_index not in [1,2]:\n",
    "        ValueError('Sub index in range {1,2}')\n",
    "        \n",
    "    if sub_index == 1:\n",
    "        # Trace out the first subsystem\n",
    "        reduced_dm = np.einsum('ijik->jk', reshaped_dm)\n",
    "        \n",
    "    if sub_index == 2:\n",
    "        # Trace out the second subsystem\n",
    "        reduced_dm = np.einsum('jiki->jk', reshaped_dm)\n",
    "    return reduced_dm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "## Check partial trace correct\n",
    "op = np.kron(rho, sigma)\n",
    "print(np.all(np.abs(partial_trace(op, 1)-sigma)<1e-5))\n",
    "print(np.all(np.abs(partial_trace(op, 2)-rho)<1e-5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0.001\n",
    "evolve = lambda s, dm: expm(-1j*s*dm)\n",
    "DME = lambda s, rho, sigma: partial_trace(evolve(s,Swap) @ np.kron(rho, sigma) @ evolve(-s,Swap),1)\n",
    "DME_s = DME(s, rho, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x16a4a3110>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGwCAYAAABiu4tnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVOElEQVR4nO3deVyU1f4H8M/ALAiyCbIpAmqKICqL+0/NFlArlyzQiuxWpnW7ilqJbZZetxazm2tG3WxRU7S8limm4gJuCKSCK6C4IOLCIMg2c35/EJMI4szA8MDweb9e87r5cOY53zODzuc+85xzZEIIASIiIiIymIXUBRARERE1VQxSREREREZikCIiIiIyEoMUERERkZEYpIiIiIiMxCBFREREZCQGKSIiIiIjyaUuwJxptVpcunQJtra2kMlkUpdDREREehBCoKCgAB4eHrCwqP2aE4OUCV26dAmenp5Sl0FERERGyM7ORtu2bWttwyBlQra2tgAq3gg7OzuJqyEiIiJ9qNVqeHp66j7Ha8MgZUKVX+fZ2dkxSBERETUx+tyWw5vNiYiIiIzEIEVERERkJAYpIiIiIiPxHqlGQKPRoKysTOoyiKBQKGBpaSl1GURETQaDlISEEMjJycHNmzelLoVIx8HBAW5ublz7jIhIDwxSEqoMUS4uLrC2tuYHF0lKCIGioiLk5uYCANzd3SWuiIio8WOQkohGo9GFKCcnJ6nLIQIAtGjRAgCQm5sLFxcXfs1HRHQfvNlcIpX3RFlbW0tcCVFVlb+TvG+PiOj+GKQkxq/zqLHh7yQRkf4YpIiIiIiMxCBFREREZCQGKSIiIiIjMUgR3YO3t7fUJRAR0T2Ua7TYeTJX6jIYpIiIiKjp+WLHGfzjm0OY+csxSetgkCKD7dq1y+irNR999BH8/f3h7e2NTZs2GfTca9euwcXFBVlZWfXW71NPPYWFCxdWaTtx4kT06NEDly5dQo8ePfD+++8b3F99M2bsNY2tLu8dEVFjkXTuOr7YcRoAEOTlKGktDFJULwYOHAiZTFbt8eyzz+rafPfdd0hPT8fRo0fx008/YdasWQb1MW/ePDzxxBPVgsD9+q6t3/fffx9z5syBWq3WHVu+fDlSUlLg4eGBlJQUg+s0hXuNvTY1jY2IqKlTF5dh8poUaAUwKrANRvRoI2k9XNmc6kwIgZSUFHzyySdVghMAtGzZUvffixcvxs8//wwLCwv4+/sb9AF/+/ZtxMTE4LfffjO479r67datG7y9vfHDDz/g1Vdf1buehnSvsd9PUxgbEZGhZv5yHBdu3IZnqxaYNcJf6nJ4RYrq7vTp0ygoKMDAgQPh5uZW5VEZZoqLi3HhwgXd/m2pqano3r277hz79+/Hww8/DGdn52pXlm7evIktW7ZALpejb9++BvV9v34BYPjw4Vi9erUpX6Iqli9fjoCAALRo0QL29vZ46KGHam1/r7Hf7zUDGn5sRESm9EvKRWxMvghLCxkWRfSArZVC6pIYpBoTIQSKSssb/CGEqFPdSUlJkMvl6Nat2z3bHDt2DFevXsWVK1dQUlKCDz74AP/85z8BVISbBx98EN27d8fu3bvx+++/o1WrVhg8eDDWrl0LBwcH7N69GyEhIQb3XVu/lXr16oWDBw+ipKSkDq+CfmJjYxEdHY333nsPJ0+eREJCAqZNm1brc2oauz6vGdCwYyMiMqXs60V4d2PFjeX/eqgjgr1aSVxRBX6114jcLtPA7/2tDd5v2qwwWCuN/1U4cuQINBpNtc2Xx44di5UrVwKo+OAfO3Yshg0bhtLSUrz++ut48MEHAQCTJk3CiBEjdDdG+/n5YezYsThw4ADCw8MBAFlZWfDw8DC479r6rdSmTRuUlJQgJycHXl5eRr8O+jh16hTatWuH0NBQXdjx96/90nRNY9fnNQMadmxERKZSrtFiytoUFJSUI9jLEa8P7ih1SToMUlRnSUlJePrppzFnzpwqxx0d/55JkZqaiieffBLffvttlTZXrlzB3r17sWPHjirHbWxsquz5dvv2bVhZWRnc9736vVOLFi0AAEVFRVWO+/n51dj+gw8+wIcffnjP8wHAoUOHaryCNn78ePz0009o1aoVrK2tkZqaig4dOtR6rrvHru9rBtx7bERETcnSXWdx+NwNtFTJsSiiB+SWjecLNQapRqSFwhJps8Ik6bcukpOTMWvWLHTseO//h5CamlrjV1hJSUnQarXV7ltKSkqqEkScnZ1x48YNg/u+V793un79OgCgdevWVY7f6+bu119/HWPGjKn1nDXNrisrK8OYMWPQs2dPrFy5Eg4ODmjfvn2t5wGqj13f1wy499iIiJqKpHM38PkfFUsdzB7pD89W1hJXVBWDVCMik8nq9BWbFDIyMnDz5k0EBgbW2i4+Pr7G41qtFkDFVZfKr7qOHj2K3bt3V1l2IDAwEN9//73Bfd+r3zsdO3YMbdu2hbOz833bAhXBRt+2d9q4cSPOnDmD7du31/jzFStWYMWKFSgtLUW3bt3w448/Aqg+dn1fM8DwsRERNSYFxWWIWpsMjVZgRA8PjApsK3VJ1TSea2PUJCUlJQEAXF1dkZOTU+VR+YEPVISEoKAgdO3aFc8884zueO/evdGiRQu89dZbOHHiBH799VeMGDECEydORL9+/XTtwsLCcPz48WpXZu7X9736vdOePXsQGhpafy/KPZSWluLy5cv47rvvkJWVhWPHjmHFihUoKyvDjRs38OWXX+LQoUM4duwYli5dqnve3WPX9zVryLEREZnCzE3HkX39Nto4tMDskV2lLqdmgkwmPz9fABD5+fnVfnb79m2RlpYmbt++LUFldbNz507h5eUlhBAiOjpaAKj2UCgUori4WAghxPXr10VQUJAoLy8XQghx48aNKuf73//+Jzp16iQUCoVo3769WLBggdBoNNX67dOnj1i+fLnuz/fr+379ClHxPtjZ2YnExMR6eGVqV1ZWJqZOnSratm0rFAqFcHV1Fc8884wQQgi1Wi28vLzEtGnTxLFjx6o99+6x6/Oa1TS2O9+7e2nKv5tEZD5+SbkovKZvFj7Rm8XBzGsN2ndtn993Y5AyoeYQpPRxv5Cgr19//VV06dKlxpBlbL+LFy8Wjz76qNE11Se1Wi1WrVolfH19xf/+978qPzN07ELUPDYGKSJqCrKvF4quM38XXtM3i0+3nmjw/g0JUvxqj0zO1tYWR48eRffu3fHUU09h8+bNRp1n2LBhmDBhAi5evFhv/SoUCnzxxRdG1VOfTp8+DVtbW0RGRmLAgAHV1n0ydOxA4xkbEZEhNFqBqWtTUVBcjsB2Dpj08ANSl1SrpnVnMzVJp0+fxgMPPIDIyEjs2bOnTotDTp48uV77feWVV4yupT7Nnj0bBw4cgLW1NQYMGIBRo0ZVa2PI2IHGMzYiIkMs23UGB7Ouw0Zp2eiWOqgJgxQZzNvbG1FRUXq31yckmIJU/Rpj1apVDdKPoe8dEVFDSsm+ic+2Vyx1MGtEV3g52Uhc0f3JhKjj/iB0T2q1Gvb29sjPz4ednV2VnxUXFyMzMxM+Pj41LjRJJBX+bhKRFG6VlOOx/+zBuWtFeLybO74YG1htkeGGUtvn990a9/UyIiIiahY+3HQc564VoY1DC8wZFSBZiDIUgxQRERFJ6tc/L2Nd0gVYyICF4d1h30IhdUl6Y5AiIiIiyVy6eRszNvwJAHjtwY7o3d7pPs9oXBikiIiISBIarcCUtSlQF5eju6cDJj/SuJc6qAmDFBEREUlixe6zOJB5HdZKS3we0QOKRr7UQU2aXsVmhpMmqbHh7yQRNYTU7JtYuO0UAOCD4f7wdm78Sx3URPIgtXTpUt006+DgYOzZs6fW9vHx8QgODoaVlRXat2+P5cuXV2sTGxsLPz8/qFQq+Pn5YePGjQb3+8ILL0Amk1V59OnTp26DvYNCUXEjXVFRUb2dk6g+VP5OVv6OEhHVt8KSckStTUG5VuCxAHc8HdxW6pKMJumCnGvXrkVUVBSWLl2K/v37Y8WKFRg6dCjS0tLQrl27au0zMzMxbNgwjB8/Ht9//z327duH1157Da1bt8bo0aMBAImJiYiIiMDs2bMxatQobNy4EeHh4di7dy969+5tUL9DhgzBN998o/uzUqmst7FbWlrCwcEBubm5AABra+smM9WTzJMQAkVFRcjNzYWDgwMsLS2lLomIzNTszWnIzCuEu70V5jahpQ5qIumCnL1790ZQUBCWLVumO9alSxeMHDkS8+bNq9Z++vTp2LRpE9LT03XHJk6ciNTUVCQmJgIAIiIioFarsWXLFl2bIUOGwNHREatXr9a73xdeeAE3b97Ezz//rPd4SkpKqmxDolar4enpec8FvYQQyMnJwc2bN/Xug8jUHBwc4Obm1qT/YSOixmvL0ct49YcjkMmAH1/ug74dGt8sPUMW5JTsilRpaSmSkpIQHR1d5XhoaCgSEhJqfE5iYiJCQ0OrHAsLC0NMTAzKysqgUCiQmJiIKVOmVGuzaNEig/vdtWsXXFxc4ODggEGDBmHOnDlwcXG555jmzZuHDz/8sNZx30kmk8Hd3R0uLi4oKyvT+3lEpqJQKHgliohM5nL+bURvOAoAmDioQ6MMUYaSLEjl5eVBo9HA1dW1ynFXV1fk5OTU+JycnJwa25eXlyMvLw/u7u73bFN5Tn37HTp0KJ5++ml4eXkhMzMT7733Hh566CEkJSVBpVLVWN+MGTMwdepU3Z8rr0jdj6WlJT+8iIjIrGm1AtN+SkX+7TIEtLHHlEc6SV1SvZB80+K7vz4QQtT6lUJN7e8+rs8579cmIiJC999du3ZFSEgIvLy88Ouvv+LJJ5+ssTaVSnXPkEVERNScrdyTgYSz19BCYYnPx/SAUi75fLd6IdkonJ2dYWlpWe3qU25ubrWrRZXc3NxqbC+Xy+Hk5FRrm8pzGtMvALi7u8PLywunT5/Wb4BEREQEADh2MR+fbDsJAJj5hB/at24pcUX1R7IgpVQqERwcjLi4uCrH4+Li0K9fvxqf07dv32rtt23bhpCQEN1U7Xu1qTynMf0CwLVr15CdnQ13d3f9BkhEREQoKi3HpDXJKNMIDPF3Q0TP+9/y0qQICa1Zs0YoFAoRExMj0tLSRFRUlLCxsRFZWVlCCCGio6NFZGSkrn1GRoawtrYWU6ZMEWlpaSImJkYoFAqxfv16XZt9+/YJS0tLMX/+fJGeni7mz58v5HK52L9/v979FhQUiGnTpomEhASRmZkpdu7cKfr27SvatGkj1Gq13uPLz88XAER+fn5dXyoiIqImKTr2T+E1fbPoPWe7uH6rROpy9GLI57ekQUoIIZYsWSK8vLyEUqkUQUFBIj4+XvezcePGiUGDBlVpv2vXLhEYGCiUSqXw9vYWy5Ytq3bOdevWic6dOwuFQiF8fX1FbGysQf0WFRWJ0NBQ0bp1a6FQKES7du3EuHHjxPnz5w0aG4MUERE1Z78fuyy8pm8W3tGbxb7TV6UuR2+GfH5Luo6UuTNkHQoiIiJzckVdjCGLduNGURkmDGyPGcO6SF2S3gz5/DaPW+aJiIio0ahc6uBGURm6trHDtNDOUpdkMgxSREREVK9i9mZi75k8WCkssCgi0GyWOqiJ+Y6MiIiIGtzxS/n4aOsJAMD7j/ujo4v5LHVQEwYpIiIiqhe3SzWYtLpiqYNH/VwxtpeZLXVQAwYpIiIiqhdzfkvD2auFcLFVYcHobs1i83MGKSIiIqqzuLQr+H7/eQDAp+Hd0cpGKXFFDYNBioiIiOokV12M6bF/AgDGD/DBgAdaS1xRw2GQIiIiIqNptQLT1qXiemEp/Nzt8EaY+S51UBMGKSIiIjLaNwlZ2HM6Dyq5Bf4ztgdUckupS2pQDFJERERklLRLaizYUrHUwbuP+6Gji63EFTU8BikiIiIyWHGZBpPXJKNUo8UjXVzwXO92UpckCQYpIiIiMticX9NxOvcWWjejpQ5qwiBFREREBvkl5SK+238OAPDJ093h1FIlcUXSYZAiIiIivZ2+UoAZG44CAF4f3BGDOjWfpQ5qwiBFREREeiksKcerPxxBUakG/To4YcqjnaQuSXIMUkRERHRfQghEbziKM7m34Gqnwn/GBsLSonneF3UnBikiIiK6r+/2n8P/Ui9BbiHDkmeC4NyM74u6E4MUERER1Sr5/A3M3pwGAIge6osQ71YSV9R4MEgRERHRPV0vLMU/fziCMo3A0K5ueOn/fKQuqVFhkCIiIqIaabUCUWtTcCm/GD7ONvjoqea7XtS9MEgRERFRjb7YcQa7T12FlcICy54Lgq2VQuqSGh0GKSIiIqpm96mrWPTHKQDAv0cGwNfNTuKKGicGKSIiIqri0s3bmLwmGUIAY3t54qngtlKX1GgxSBEREZFOabkW//zxCG4UlcHfww4zn/CXuqRGjUGKiIiIdOZtSUfy+Zuws5Jj2bPBsFJYSl1So8YgRURERACAzX9ewjf7sgAAC8N7oJ2TtbQFNQEMUkRERIQzubcwff2fAIBXH+yAR/xcJa6oaWCQIiIiauaKSsvx2g9JKCzVoE/7VpjGzYj1xiBFRETUjAkh8M7GYzh15RZcbCs2I5ZbMh7oi68UERFRM/bDgfPYmHwRlhYyLH4mCC62VlKX1KQwSBERETVTf164iVn/q9iMePqQzujlw82IDcUgRURE1AzdLCrFq98fQalGi1A/V4wf0F7qkpokBikiIqJmRqsVmLI2BRdv3oaXkzU+fro7NyM2EoMUERFRM7Ms/ix2nrwKldwCS58Ngn0LbkZsLAYpIiKiZmTfmTx8uu0kAGD2iK7w97CXuKKmjUGKiIiomcjJL8ak1cnQCiA8pC3Ce3pKXVKTxyBFRETUDJRptHj9xyO4VliKLu52mDWiq9QlmQUGKSIiomZgwZYTOHzuBmxVcix7NoibEdcTBikiIiIzt+XoZXy1NxMA8El4d3g720hckflgkCIiIjJjmXmFePOvzYgnDGyPMH83iSsyLwxSREREZup2qQavfp+EWyXl6OXdCm+GdZa6JLPDIEVERGSGhBB49+djOJFTAOeWKix+hpsRmwJfUSIiIjO09lA2Yo9cgIUM+GJsIFzsuBmxKTBIERERmZljF/Px/qbjAIA3wjqjbwcniSsyXwxSREREZiS/qAyv/pCE0nItHunigokDO0hdklljkCIiIjITWq3AtHUpyL5+G56tWuDTp3vAwoKbEZsSgxQREZGZWLE7A9vTc6GUW2DZs8Gwt+ZmxKbGIEVERGQGEs9ew8dbTwAAPhzuj65tuBlxQ2CQIiIiauJy1cX411+bEY8Oaosx3Iy4wTBIERERNWHlGi1eX52MvFsl8HWzxb9HdoVMxvuiGgqDFBERURP28daTOJh5HS1Vcix9NggtlNyMuCFJHqSWLl0KHx8fWFlZITg4GHv27Km1fXx8PIKDg2FlZYX27dtj+fLl1drExsbCz88PKpUKfn5+2LhxY536nTBhAmQyGRYtWmTw+IiIiExl6/EcrNidAQD4+KluaN+6pcQVNT+SBqm1a9ciKioK77zzDpKTkzFgwAAMHToU58+fr7F9ZmYmhg0bhgEDBiA5ORlvv/02Jk2ahNjYWF2bxMREREREIDIyEqmpqYiMjER4eDgOHDhgVL8///wzDhw4AA8Pj/p/AYiIiIx07loh3liXCgB46f98MDTAXeKKmieZEEJI1Xnv3r0RFBSEZcuW6Y516dIFI0eOxLx586q1nz59OjZt2oT09HTdsYkTJyI1NRWJiYkAgIiICKjVamzZskXXZsiQIXB0dMTq1asN6vfixYvo3bs3tm7disceewxRUVGIioq653hKSkpQUlKi+7NarYanpyfy8/NhZ2dnwCtDRER0b0Wl5Ri9LBHpl9UI8XLE6lf6QMF99OqNWq2Gvb29Xp/fkr3qpaWlSEpKQmhoaJXjoaGhSEhIqPE5iYmJ1dqHhYXh8OHDKCsrq7VN5Tn17Ver1SIyMhJvvvkm/P399RrTvHnzYG9vr3t4enLWBBER1S+tVmDK2hSkX1bDuaUSi58JYoiSkGSvfF5eHjQaDVxdXascd3V1RU5OTo3PycnJqbF9eXk58vLyam1TeU59+12wYAHkcjkmTZqk95hmzJiB/Px83SM7O1vv5xIREeljYdwpbD1+BUpLC6yIDIGbPTcjlpJc6gLunqIphKh12mZN7e8+rs85a2uTlJSEzz//HEeOHDFoCqlKpYJKpdK7PRERkSF+SbmIxTvPAADmjw5AsJejxBWRZFeknJ2dYWlpWe3qU25ubrWrRZXc3NxqbC+Xy+Hk5FRrm8pz6tPvnj17kJubi3bt2kEul0Mul+PcuXOYNm0avL29jR4zERGRsZLP38Cb6/8EAEwc1AFPBrWVuCICJAxSSqUSwcHBiIuLq3I8Li4O/fr1q/E5ffv2rdZ+27ZtCAkJgUKhqLVN5Tn16TcyMhJ//vknUlJSdA8PDw+8+eab2Lp1q/GDJiIiMsKlm7cxflUSSsu1eNTPFW+FdZa6JPqLpF/tTZ06FZGRkQgJCUHfvn3x5Zdf4vz585g4cSKAinuOLl68iFWrVgGomKG3ePFiTJ06FePHj0diYiJiYmJ0s/EAYPLkyRg4cCAWLFiAESNG4JdffsH27duxd+9evft1cnLSXeGqpFAo4Obmhs6d+ctLREQNp6i0HC9/e1i3cvmiiB6wsODK5Y2FpEEqIiIC165dw6xZs3D58mV07doVv/32G7y8vAAAly9frrK2k4+PD3777TdMmTIFS5YsgYeHB/7zn/9g9OjRujb9+vXDmjVr8O677+K9995Dhw4dsHbtWvTu3VvvfomIiBoDrVZg6tpUpP01Q++rcSGwUUl+ezPdQdJ1pMydIetQEBER3e3TbSfxxY4zUFpaYPUrvRHs1UrqkpqFJrGOFBEREd3bLykX8cWOihl6854MYIhqpBikiIiIGpm7Z+iNDuYMvcaKQYqIiKgRuXTzNl75rmKG3iNdOEOvsWOQIiIiaiSKSssxftVhXC34a4beGM7Qa+wYpIiIiBqByhl6xy+p4WRTMUOvJWfoNXoMUkRERI3AZ9tP4ffjOX/toReMto7WUpdEemCQIiIiktjdM/RCvDlDr6lgkCIiIpJQSvZN3Qy9CYPac4ZeE8MgRUREJJHL+bcxftXhv2boueCtMF+pSyIDMUgRERFJoHIPvb9n6AXCkjP0mhwGKSIiogam1QpM++nvGXorn+cMvaaKQYqIiKiBLdp+CluO/T1Dz7MVZ+g1VQxSREREDWhT6iX8568ZenM5Q6/JY5AiIiJqICnZN/HmulQAwISB7fEUZ+g1eQxSREREDaByhl5JuRYP+7rgrSGcoWcOGKSIiIhM7M499Dq72uLzsZyhZy4YpIiIiExIqxV4Y10qjl3kHnrmiEGKiIjIhBb9cRq/Hc2BwlKG5ZyhZ3YYpIiIiExkU+ol/OeP0wCAuaMC0JMz9MwOgxQREZEJ3D1D7+kQT4krIlNgkCIiIqpnOfnFeIUz9JoFBikiIqJ6dLtUg/GrDiOXM/SaBQYpIiKieqLVCkxbl4KjF/PRijP0mgUGKSIionpy5ww97qHXPDBIERER1YP/cYZes2RQkCorK8PgwYNx6tQpU9VDRETU5KRm38Qbf83Qe4Uz9JoVg4KUQqHAsWPHIJPxpjkiIiKgYobenXvoTecMvWbF4K/2nn/+ecTExJiiFiIioiblzhl6nVxbYtGYHpyh18wYPJWgtLQUX331FeLi4hASEgIbG5sqP1+4cGG9FUdERNRYVe6hVzlDL2ZcT9haKaQuixqYwUHq2LFjCAoKAoBq90rxKz8iImouPv/jNH49erliD73nOEOvuTI4SO3cudMUdRARETUZ6w5n4/O/ZujNGRWAXj6coddc1Wn5gwsXLuDixYv1VQsREVGjtz3tCqI3HAUATBjUHuGcodesGRyktFotZs2aBXt7e3h5eaFdu3ZwcHDA7NmzodVqTVEjERFRo5B07jr++eMRaLQCo4PaIpoz9Jo9g7/ae+eddxATE4P58+ejf//+EEJg3759+OCDD1BcXIw5c+aYok4iIiJJnbpSgBf/W7HMwUO+Lpg/OoD3BhNkQghhyBM8PDywfPlyDB8+vMrxX375Ba+99hq/6ruDWq2Gvb098vPzYWdnJ3U5RERkpIs3b2P00gTkqIsR7OWI71/qjRZKS6nLIhMx5PPb4K/2rl+/Dl/f6pcyfX19cf36dUNPR0RE1KhdLyxFZMwB5KiL8YBLS8SMC2GIIh2Dg1T37t2xePHiascXL16M7t2710tRREREjUFhSTn+8d9DyLhaCA97K6x6qRccrJVSl0WNiMH3SH300Ud47LHHsH37dvTt2xcymQwJCQnIzs7Gb7/9ZooaiYiIGlyZRotXfziC1OybcLBWYNVLveBu30LqsqiRMfiK1KBBg3Dq1CmMGjUKN2/exPXr1/Hkk0/i5MmTGDBggClqJCIialBarcCb61Kx+9RVtFBY4psXeqKji63UZVEjZNAVqbKyMoSGhmLFihWcnUdERGZJCIE5v6Xj55RLkFvIsOy5IAS2c5S6LGqkDLoipVAocOzYMU73JCIis7U8PgMxezMBAB8/3Q0PdnaRuCJqzAz+au/5559HTEyMKWohIiKS1E+Hs7Hg9xMAgHcf64JRgW0lrogaO4NvNi8tLcVXX32FuLg4hISEwMbGpsrPFy5cWG/FERERNZTtaVcw446tX14e0F7iiqgpMDhIHTt2DEFBQQCAU6dOVfkZv/IjIqKm6HDW31u/PBXMrV9IfwYFKY1Ggw8++AABAQFo1Yo7XRMRUdN3MqcAL/73EErKtXjY1wXzn+TWL6Q/g+6RsrS0RFhYGPLz801VDxERUYO5cKMIz399AOricoR4OWLxM0GQWxp8+zA1Ywb/tgQEBCAjI8MUtRARETWY64WleP7rg7iiLkEn15b4ilu/kBEMDlJz5szBG2+8gc2bN+Py5ctQq9VVHkRERI1dYUk5/vHNQd3WL9++yK1fyDgyIYQw5AkWFn9nrzu/QxZCQCaTQaPR1F91TZwhu0cTEVHDKC3X4qVvD2HP6Tw4WiuwbmI/dHRpKXVZ1IgY8vlt8Ky9nTt3Gl0YERGRlLRagTfXp2LP6Ty0UFji6xd6MkRRnRgcpAYNGmSKOoiIiExKCIF//5qOX/7a+mV5ZDC3fqE6M2pqwp49e/Dcc8+hX79+uHjxIgDgu+++w969ew0+19KlS+Hj4wMrKysEBwdjz549tbaPj49HcHAwrKys0L59eyxfvrxam9jYWPj5+UGlUsHPzw8bN240uN8PPvgAvr6+sLGxgaOjIx555BEcOHDA4PEREVHjsCz+LL7eV7H1yydPd8egTq0lrojMgcFBKjY2FmFhYWjRogWOHDmCkpISAEBBQQHmzp1r0LnWrl2LqKgovPPOO0hOTsaAAQMwdOhQnD9/vsb2mZmZGDZsGAYMGIDk5GS8/fbbmDRpEmJjY3VtEhMTERERgcjISKSmpiIyMhLh4eFVQpA+/Xbq1AmLFy/G0aNHsXfvXnh7eyM0NBRXr141aIxERCS9nw5l46PfTwKo2PplZGAbiSsisyEM1KNHD/Htt98KIYRo2bKlOHv2rBBCiOTkZOHq6mrQuXr16iUmTpxY5Zivr6+Ijo6usf1bb70lfH19qxybMGGC6NOnj+7P4eHhYsiQIVXahIWFiTFjxhjdrxBC5OfnCwBi+/bttQ+qhufk5+fr/RwiIqpf247nCJ/ozcJr+mYx77d0qcuhJsCQz2+Dr0idPHkSAwcOrHbczs4ON2/e1Ps8paWlSEpKQmhoaJXjoaGhSEhIqPE5iYmJ1dqHhYXh8OHDKCsrq7VN5TmN6be0tBRffvkl7O3t0b1793uOqaSkhMtBEBE1IoeyruP1H49AK4Cng9ti+pDOUpdEZsbgIOXu7o4zZ85UO7537160b6//Bo95eXnQaDRwdXWtctzV1RU5OTk1PicnJ6fG9uXl5cjLy6u1TeU5Del38+bNaNmyJaysrPDZZ58hLi4Ozs7O9xzTvHnzYG9vr3t4enrW8goQEZEpnchR46W/tn55pIsL5nHrFzIBg4PUhAkTMHnyZBw4cAAymQyXLl3CDz/8gDfeeAOvvfaawQXc/Ust/lqPypD2dx/X55z6tBk8eDBSUlKQkJCAIUOGIDw8HLm5ufesbcaMGcjPz9c9srOz79mWiIhM58KNIoz7+qBu65cvxnLrFzINg5c/eOutt5Cfn4/BgwejuLgYAwcOhEqlwhtvvIHXX39d7/M4OzvD0tKy2lWg3NzcaleLKrm5udXYXi6Xw8nJqdY2lec0pF8bGxt07NgRHTt2RJ8+ffDAAw8gJiYGM2bMqLE+lUoFlUp1n5ETEZEpXbtVgudj/t76JWZcT279QiZjVDyfM2cO8vLycPDgQezfvx9Xr17F7NmzDTqHUqlEcHAw4uLiqhyPi4tDv379anxO3759q7Xftm0bQkJCoFAoam1TeU5j+q0khNDNUiQiosansKQcL/73EDLyCtHGoQVWvdgb9tYKqcsic2ba+95rt2bNGqFQKERMTIxIS0sTUVFRwsbGRmRlZQkhhIiOjhaRkZG69hkZGcLa2lpMmTJFpKWliZiYGKFQKMT69et1bfbt2ycsLS3F/PnzRXp6upg/f76Qy+Vi//79evd769YtMWPGDJGYmCiysrJEUlKSeOmll4RKpRLHjh3Te3yctUdE1HBKyjTiua/2C6/pm0WPD7eK01cKpC6JmihDPr8lDVJCCLFkyRLh5eUllEqlCAoKEvHx8bqfjRs3TgwaNKhK+127donAwEChVCqFt7e3WLZsWbVzrlu3TnTu3FkoFArh6+srYmNjDer39u3bYtSoUcLDw0MolUrh7u4uhg8fLg4ePGjQ2BikiIgahkajFf/68Yjwmr5ZdHlvi0g+f0PqkqgJM+Tz2+BNi0l/3LSYiMj0hBCYtTkN3+zLgtxChq9f6ImBXLWc6sCQz29OYSAioiZt6a6z+GZfFgDg0/DuDFHUoBikiIioyVp76Dw+3lqx9ct7j/thRA9u/UINy+DlDwDg4sWL2LdvH3Jzc6HVaqv8bNKkSfVSGBERUW22Hc/BjA1HAQCvPtgBL/2fj8QVUXNkcJD65ptvMHHiRCiVSjg5OVVbCJNBioiITC3hbB7+tToZWgGEh7TFW2Hc+oWkYXCQev/99/H+++9jxowZsLDgN4NERNSwEs7m4UXd1i+umDuKW7+QdAwOUkVFRRgzZgxDFBERNbiEM3l48dtDKC7TYnDn1lj8TCC3fiFJGfzb99JLL2HdunWmqIWIiOie7g5Ry54LhpWCW7+QtAxeR0qj0eDxxx/H7du3ERAQoNuapdLChQvrtcCmjOtIERHVj31n8vDSHSFqeWQwVHKGKDINQz6/Df5qb+7cudi6dSs6d664se/um82JiIjqE0MUNWYGB6mFCxfi66+/xgsvvGCCcoiIiP6278zfN5Y/5OuCZc8FMURRo2LwPVIqlQr9+/c3RS1EREQ6DFHUFBgcpCZPnowvvvjCFLUQEREBAPae/jtEPcwQRY2YwV/tHTx4EDt27MDmzZvh7+9f7WbzDRs21FtxRETU/Ow9XXFPVGWIWsoQRY2YwUHKwcEBTz75pClqISKiZo4hipoag4JUeXk5HnzwQYSFhcHNzc1UNRERUTO05/RVvPzt4b9WLHfBkmcZoqjxM+geKblcjldffRUlJSWmqoeIiJohhihqqgy+2bx3795ITk42RS1ERNQM7T7FEEVNl8H3SL322muYNm0aLly4gODgYNjY2FT5ebdu3eqtOCIiMm+7T13F+FV/h6ilzwZDKefeedR0GLxFTE2bFctkMgghIJPJoNFo6q24po5bxBAR3dvuU1fx8qrDKC3X4pEurlj6bBBDFDUKJt0iJjMz0+jCiIiIAIYoMh8GBykvLy9T1EFERM1E/F9f55WWa/GonyuWPMMQRU2XwUEKAM6ePYtFixYhPT0dMpkMXbp0weTJk9GhQ4f6ro+IiMwIQxSZG4N/e7du3Qo/Pz8cPHgQ3bp1Q9euXXHgwAH4+/sjLi7OFDUSEZEZ2HUylyGKzI7BN5sHBgYiLCwM8+fPr3I8Ojoa27Ztw5EjR+q1wKaMN5sTEVXYdTIXr3yXhNJyLUL9XLGYIYoaMUM+vw3+LU5PT8dLL71U7fiLL76ItLQ0Q09HRERmjiGKzJnBv8mtW7dGSkpKteMpKSlwcXGpj5qIiMhM3BmiwvwZosj8GHyz+fjx4/HKK68gIyMD/fr1g0wmw969e7FgwQJMmzbNFDUSEVETtPNkLiasSkKppiJEfTGWIYrMj8H3SAkhsGjRInz66ae4dOkSAMDDwwNvvvkmJk2aBJlMZpJCmyLeI0VEzdXdIWrxM0FQWDJEUdNgyOe3wUHqTgUFBQAAW1tbY09h1hikiKg5Yoiips6kK5vfiQGKiIjutPNELiZ8xxBFzYfBv91XrlxBZGQkPDw8IJfLYWlpWeVBRETN050haoi/G0MUNQsGX5F64YUXcP78ebz33ntwd3fnPVFERIQdJ65g4ndHdCHqi2cCGaKoWTA4SO3duxd79uxBjx49TFAOERE1NXeGqKFd3fCfsQxR1HwY/Jvu6emJOtyfTkREZuSPdIYoat4M/m1ftGgRoqOjkZWVZYJyiIioqfgj/Qpe/b4iRA0LYIii5sngr/YiIiJQVFSEDh06wNraGgqFosrPr1+/Xm/FERFR47TteA5e/zFZF6I+H8MQRc2TwUFq0aJFJiiDiIiaih8PnMe7Px+FVoAhipo9g4PUuHHj9Go3f/58TJw4EQ4ODoZ2QUREjZAQAou2n8bnf5wGAISHtMXcUQGQM0RRM2ay3/65c+fyaz4iIjNRrtHi7Y1HdSHqXw91xILR3RiiqNmr08rmteHMPiIi83C7VIN/rU7G9vQrkMmAWSO6IrKPl9RlETUKJgtSRETU9N0oLMVL3x7CkfM3oZRb4D9jemBIV3epyyJqNBikiIioRhduFGHc1wdx9moh7Kzk+GpcT/TyaSV1WUSNCoMUERFVcyJHjXFfH8QVdQnc7a3w7Yu90MmVG9UT3Y1BioiIqtifcQ3jVx1GQXE5HnBpiW9f7AUPhxZSl0XUKJksSA0YMAAtWvAvHhFRU/Lrn5cxZW0KSjVa9PJuhZXPh8DeWnH/JxI1U0YFKY1Gg40bNyI9PR0ymQy+vr4YOXIk5PK/T/fbb7/VW5FERGR6/92XiQ83p0EIIMzfFZ+PCYSVwlLqsogaNYOD1LFjxzBixAjk5OSgc+fOAIBTp06hdevW2LRpEwICAuq9SCIiMh0hBD7eehJLd50FADzXpx0+HN4VlhYyiSsjavwMXknt5Zdfhr+/Py5cuIAjR47gyJEjyM7ORrdu3fDKK6+YokYiIjKRMo0Wb6z7Uxei3gjthNkjGKKI9GXwFanU1FQcPnwYjo6OumOOjo6YM2cOevbsWa/FERGR6RSWlOO1H44g/tRVWFrIMG9UAMJ7ekpdFlGTYvAVqc6dO+PKlSvVjufm5qJjx471UhQREZnWtVsleGblfsSfugorhQVWPh/MEEVkBIOvSM2dOxeTJk3CBx98gD59+gAA9u/fj1mzZmHBggVQq9W6tnZ2dvVXKRER1Yvz14rw/NcHkHWtCI7WCsS80BNB7Rzv/0QiqsbgK1KPP/440tLSEB4eDi8vL3h5eSE8PBzHjh3DE088AUdHRzg4OFT56q82S5cuhY+PD6ysrBAcHIw9e/bU2j4+Ph7BwcGwsrJC+/btsXz58mptYmNj4efnB5VKBT8/P2zcuNGgfsvKyjB9+nQEBATAxsYGHh4eeP7553Hp0iW9xkRE1Fgdu5iPJ5ftQ9a1IrR1bIH1r/ZjiCKqA4OvSO3cubPeOl+7di2ioqKwdOlS9O/fHytWrMDQoUORlpaGdu3aVWufmZmJYcOGYfz48fj++++xb98+vPbaa2jdujVGjx4NAEhMTERERARmz56NUaNGYePGjQgPD8fevXvRu3dvvfotKirCkSNH8N5776F79+64ceMGoqKiMHz4cBw+fLjexk9E1JD2nL6Kid8lobBUgy7udvj2Hz3hYmcldVlETZpMCCGk6rx3794ICgrCsmXLdMe6dOmCkSNHYt68edXaT58+HZs2bUJ6erru2MSJE5GamorExEQAQEREBNRqNbZs2aJrM2TIEDg6OmL16tVG9QsAhw4dQq9evXDu3LkaQ15N1Go17O3tkZ+fz685iUhSPydfxBvrUlGuFejXwQnLI4NhZ8WFNolqYsjnt8Ff7QHAzZs38emnn+Lll1/G+PHj8dlnnyE/P9+gc5SWliIpKQmhoaFVjoeGhiIhIaHG5yQmJlZrHxYWhsOHD6OsrKzWNpXnNKZfAMjPz4dMJoODg8M925SUlECtVld5EBFJbeXuDEStTUG5VuCJ7h745h89GaKI6onBQerw4cPo0KEDPvvsM1y/fh15eXlYuHAhOnTogCNHjuh9nry8PGg0Gri6ulY57urqipycnBqfk5OTU2P78vJy5OXl1dqm8pzG9FtcXIzo6Gg888wztSbTefPmwd7eXvfw9OQMGCKSjlYrMHtzGub8VnEV/8X+Pvg8ogdUcq5WTlRfDA5SU6ZMwfDhw5GVlYUNGzZg48aNyMzMxOOPP46oqCiDC5DJqi76JoSodux+7e8+rs859e23rKwMY8aMgVarxdKlS2sZCTBjxgzk5+frHtnZ2bW2JyIylZJyDaLWpiBmbyYA4O1hvnjv8S6w4EKbRPXK4JvNDx8+jJUrV1bZV08ul+Ott95CSEiI3udxdnaGpaVltatAubm51a4WVXJzc6uxvVwuh5OTU61tKs9pSL9lZWUIDw9HZmYmduzYcd/vSVUqFVQqVa1tiIhMraC4DBO+S0LC2WuQW8jwydPdMTKwjdRlEZklg69I2dnZ4fz589WOZ2dnw9bWVu/zKJVKBAcHIy4ursrxuLg49OvXr8bn9O3bt1r7bdu2ISQkBAqFotY2lefUt9/KEHX69Gls375dF9SIiBqzXHUxIlbsR8LZa7BRWuLrF3oyRBGZkjDQv/71L9G2bVuxZs0acf78eZGdnS1Wr14t2rZtKyZPnmzQudasWSMUCoWIiYkRaWlpIioqStjY2IisrCwhhBDR0dEiMjJS1z4jI0NYW1uLKVOmiLS0NBETEyMUCoVYv369rs2+ffuEpaWlmD9/vkhPTxfz588Xcrlc7N+/X+9+y8rKxPDhw0Xbtm1FSkqKuHz5su5RUlKi9/jy8/MFAJGfn2/Q60JEZIyzuQWi//w/hNf0zSJ49jbxZ/ZNqUsiapIM+fw2OEiVlJSISZMmCaVSKSwsLISFhYVQqVQiKipKFBcXG1zskiVLhJeXl1AqlSIoKEjEx8frfjZu3DgxaNCgKu137dolAgMDhVKpFN7e3mLZsmXVzrlu3TrRuXNnoVAohK+vr4iNjTWo38zMTAGgxsfOnTv1HhuDFBE1lCPnrovAWduE1/TNYtBHO0RW3i2pSyJqsgz5/DZ6HamioiKcPXsWQgh07NgR1tbW9XOJzIxwHSkiagg7TlzBP39Ixu0yDbq1tcfXL/SEc0ver0lkLEM+vw2+2byStbU1AgICjH06ERHVg58OZ2PGhqPQaAUGdmqNZc8GwUZl9D/tRGQg/m0jImqChBBYsvMMPtl2CgDwZFAbLBjdDQpLo9ZZJiIjMUgRETUxpeVafPC/4/jxQMUM6lcf7IC3wjrXugYfEZkGgxQRUROSqy7Gqz8cQdK5G5DJgPcf98M/+vtIXRZRs8UgRUTURCSdu45Xvz+C3IIS2FrJ8fmYHnjIt+YFjImoYTBIERE1ckII/HjwPD7YdBxlGoFOri2xIjIEPs42UpdG1OwxSBERNWLFZRrM/OU41h6u2LtzWIAbPn6qO2fmETUS/JtIRNRIXc6/jYnfH0Fq9k1YyIA3w3wxcVB73lRO1IgwSBERNUIHMq7hnz8eQd6tUti3UOCLsYEY2Km11GUR0V0YpIiIGhEhBL5NyMK/f01HuVbA180WX0aGoJ0Td48gaowYpIiIGoniMg3e3ngUG45cBAAM7+6B+aMDYK3kP9VEjRX/dhIRNQIXbhRh4vdJOHZRDUsLGWYM9cVL/+fD+6GIGjkGKSIiiSWczcPrPybjemEpWtkosXhsIPp1dJa6LCLSA4MUEZFEhBCI2ZuJeVtOQKMV6NrGDsufC0ZbR94PRdRUMEgREUngdqkG02P/xKbUSwAqNh2eOyoAVgpLiSsjIkMwSBERNbDs60V45bskpF9WQ24hw3uP++H5vl68H4qoCWKQIiJqQLtPXcWkNcm4WVQG55ZKLHkmCL3bO0ldFhEZiUGKiKgBCCGwPD4DH289Aa0Auns6YPlzQXC3byF1aURUBwxSREQmVlhSjjfXp+K3ozkAgIgQT8wa6Q+VnPdDETV1DFJERCaUlVeIV747jFNXbkFhKcMHw/3xTK92vB+KyEwwSBERmciOE1cweU0KCorL4WKrwrLnghHs5Sh1WURUjxikiIjqmVYrsGTnGSzcfgpCAMFejlj2bBBc7KykLo2I6hmDFBFRPSooLsO0n1KxLe0KAOC5Pu3w/uP+UMotJK6MiEyBQYqIqJ6cyb2FCd8dxtmrhVBaWuDfI7sivKen1GURkQkxSBER1YNtx3Mw9adU3Coph5udFZZHBqOHp4PUZRGRiTFIERHVgVYrsGj7KfxnxxkAQC+fVljyTBBa26okroyIGgKDFBGRkfJvl2HK2hTsOJELAPhHf2+8PawLFJa8H4qouWCQIiIywqkrBZjwXRIy8wqhkltg3pMBeDKordRlEVEDY5AiIjKAEALf7z+HOb+lo7hMizYOLbAiMhhd29hLXRoRSYBBiohIT1fUxXhr/Z+IP3UVADDgAWd8PiYQrWyUEldGRFJhkCIi0sNvRy/j7Y1HcbOoDCq5BWYM9cXzfb1hYcGtXoiaMwYpIqJaqIvL8MEvx7Eh+SIAoGsbOyyK6IGOLrYSV0ZEjQGDFBHRPSSczcMbP6XiUn4xLGTAPwd3xL8eeoCrlBORDoMUEdFdiss0+GTrSXy1NxMA4OVkjYXhPbjhMBFVwyBFRHSH45fyMWVtCk5duQUAGNurHd59rAtsVPznkoiq478MREQANFqBL3dnYGHcSZRpBJxbqrBgdAAe7uIqdWlE1IgxSBFRs3f+WhGmrUvBoawbAIBQP1fMezIATi25zQsR1Y5BioiaLSEE1h2+gA//dxyFpRq0VMkx8wk/PBXcFjIZlzUgovtjkCKiZinvVglmbDiKuLQrAIBe3q3waXh3eLaylrgyImpKGKSIqNnZnnYF0Rv+RN6tUigsZZgW2hnjB7SHJRfXJCIDMUgRUbNxq6Qc/96chjWHsgEAvm62WBjeA34edhJXRkRNFYMUETULSeeuY8raVJy/XgSZDBg/oD2mPtoJVgpLqUsjoiaMQYqIzFppuRaf/3EKy3adhVYAbRxa4NPw7ujT3knq0ojIDDBIEZHZOn2lAFFrU3D8khoAMDqoLWYO94OdlULiyojIXDBIEZHZ0WoFvknIwoLfT6C0XAtHawXmjgrA0AB3qUsjIjPDIEVEZuXSzdt4Y10qEs5eAwA82Lk1PhrdDS52VhJXRkTmiEGKiMyCEAKbUi/h3Z+PoaC4HC0UlnjnsS54tnc7Lq5JRCbDIEVETd7NolK8+/MxbP7zMgCgh6cDPovoAR9nG4krIyJzxyBFRE3a7lNX8eb6VFxRl8DSQobJDz+A1x7sALmlhdSlEVEzwCBFRE3S7VIN5m9Jx7eJ5wAA7Vvb4LPwHuju6SBtYUTUrDBIEVGTs+9MHt775RgyrhYCAMb19UL00C5ooeTimkTUsCS/9r106VL4+PjAysoKwcHB2LNnT63t4+PjERwcDCsrK7Rv3x7Lly+v1iY2NhZ+fn5QqVTw8/PDxo0bDe53w4YNCAsLg7OzM2QyGVJSUuo0TiKqu4s3b+O1H5Lw7FcHkHG1EK52Kqx6sRc+HNGVIYqIJCFpkFq7di2ioqLwzjvvIDk5GQMGDMDQoUNx/vz5GttnZmZi2LBhGDBgAJKTk/H2229j0qRJiI2N1bVJTExEREQEIiMjkZqaisjISISHh+PAgQMG9VtYWIj+/ftj/vz5pnsBiEgvxWUaLN5xGg9/ugu/Hc2BhQx4oZ83tkUNwsBOraUuj4iaMZkQQkjVee/evREUFIRly5bpjnXp0gUjR47EvHnzqrWfPn06Nm3ahPT0dN2xiRMnIjU1FYmJiQCAiIgIqNVqbNmyRddmyJAhcHR0xOrVqw3uNysrCz4+PkhOTkaPHj0MGp9arYa9vT3y8/NhZ8dNUYmMsePEFXz4vzScu1YEAOjl0wofDvdHF3f+nSIi0zDk81uyK1KlpaVISkpCaGholeOhoaFISEio8TmJiYnV2oeFheHw4cMoKyurtU3lOY3pV18lJSVQq9VVHkRknHPXCvHSfw/hxf8exrlrRXC1U+HzMT2w9pU+DFFE1GhIdrN5Xl4eNBoNXF1dqxx3dXVFTk5Ojc/JycmpsX15eTny8vLg7u5+zzaV5zSmX33NmzcPH374YZ3OQdTc3S7VYOmuM1gRn4FSjRZyCxleGuCDfz30AFqqOD+GiBoXyf9VunvFYSFErasQ19T+7uP6nNPQfvUxY8YMTJ06VfdntVoNT0/POp2TqLkQQmDLsRzM+TUdF2/eBgAMeMAZM5/wR0eXlhJXR0RUM8mClLOzMywtLatdBcrNza12taiSm5tbje3lcjmcnJxqbVN5TmP61ZdKpYJKparTOYiaozO5BZi56Tj2nanYH6+NQwu897gfwvxdub0LETVqkt0jpVQqERwcjLi4uCrH4+Li0K9fvxqf07dv32rtt23bhpCQECgUilrbVJ7TmH6JyDQKissw59c0DFm0B/vOXINSboFJDz+A7VMHYUhXN4YoImr0JP1qb+rUqYiMjERISAj69u2LL7/8EufPn8fEiRMBVHxVdvHiRaxatQpAxQy9xYsXY+rUqRg/fjwSExMRExOjm40HAJMnT8bAgQOxYMECjBgxAr/88gu2b9+OvXv36t0vAFy/fh3nz5/HpUuXAAAnT54EUHHFy83NzeSvDZE5E0Lg55SLmPvbCVwtKAEAPNLFFe8/7od2TtYSV0dEZAAhsSVLlggvLy+hVCpFUFCQiI+P1/1s3LhxYtCgQVXa79q1SwQGBgqlUim8vb3FsmXLqp1z3bp1onPnzkKhUAhfX18RGxtrUL9CCPHNN98IANUeM2fO1Hts+fn5AoDIz8/X+zlE5u7YxZviqWX7hNf0zcJr+mYx6KMdYkf6FanLIiLSMeTzW9J1pMwd15Ei+tvNolIsjDuF7/efg1YALRSWeP2hjnh5gA9Ucq5KTkSNhyGf35LP2iMi86bVCvx0OBsfbT2J64WlAIDHurnjnWFd4OHQQuLqiIjqhkGKiEwmJfsmZv5yDKkX8gEAD7i0xIcj/NGvg7PElRER1Q8GKSKqd3m3SvDx7yex9nA2AMBWJUfUo53wfF8vKCwl3yudiKjeMEgRUb0p12jx/f5z+DTuFAqKywEAo4PaYvrQznCxtZK4OiKi+scgRUT14kDGNczcdBwncgoAAP4edpg1wh/BXq0kroyIyHQYpIioTq6oizH3t3T8klKx5pp9CwXeDOuMsb3awdKCC2oSkXljkCIio5SWa/HNvkz854/TKCzVQCYDxvZqhzdCO6OVjVLq8oiIGgSDFBEZRKut2Fz407iTyLhaCAAIbOeAWcO7IqCtvcTVERE1LAYpItKLRiuw+c9LWLzjDE7n3gIAOLdUYvoQX4wOagsLfo1HRM0QgxQR1apco8Wm1IoAlZFXcQXK1kqOF/v74MX/84F9C4XEFRIRSYdBiohqVKbRYmPyRSzZeQbnrhUBqLiR/OX/88G4/t6ws2KAIiJikCKiKkrLtYg9cgFLdp7BhRu3AQCtbJR4eYAPIvt4wZYBiohIh0GKiAAAJeUa/HT4ApbtPINL+cUAKu6BemVgezzb2ws2Kv5zQUR0N/7LSNTMFZdpsObgeSyPz0COuiJAtbZVYeKgDnimVzu0UFpKXCERUePFIEXUTN0u1eCHA+ewYncGrhaUAADc7Kzw6oMdENHTE1YKBigiovthkCJqZgpLyvH9/nNYuScDebdKAQBtHFrg1Qc74OmQtlDJGaCIiPTFIEXUTBQUl2FV4jl8tScDN4rKAACerVrg9cEdMSqwLZRyC4krJCJqehikiMxc/u0yfJuQhZi9mci/XRGgvJ2s8c/BHTEysA0UlgxQRETGYpAiMlM3i0rx9b4sfLMvEwXF5QCA9q1t8K+HOuKJbh6QM0AREdUZgxSRmbleWIqYvRn4NuEcbpVUBKhOri3x+kMP4LEAd1hyKxcionrDIEVkJvJulWDlngx8l3gORaUaAICvmy0mP/wAwvzduBceEZEJMEgRNXG5BcX4Mj4D3x84h+IyLQCgaxs7THroATzSxZUBiojIhBikiJqonPxiLI8/i9UHz6OkvCJAdW9rj8mPPIDBnV0gkzFAERGZGoMUURMihMChrBtYc+g8NqdeRqmmIkAFtXPA5Ec6YeADzgxQREQNiEGKqAm4WlCCDUcuYO2hbGTkFeqO9/JuhcmPPIB+HZwYoIiIJMAgRdRIabQCu09fxdqD2diefgXlWgEAsFZaYnh3D0T09ERgO0eJqyQiat4YpIgamQs3ivDT4QtYdzgbl/OLdccD2zlgTE9PPNbNAy1V/KtLRNQY8F9jokagtFyLuLQrWHPoPPaeyYOouPgEB2sFngxsi4ienujsZittkUREVA2DFJGEzuQWYO2hbMQeuYjrhaW64//X0RkRPT0R6u/KTYSJiBoxBimiBlZUWo5f/7yMtYeycfjcDd1xVzsVng72RHiIJ9o5WUtYIRER6YtBiqgBCCFw9GI+1hzKxqaUS7qtWywtZHjI1wVjenpiUKfW3P+OiKiJYZAiMqH8ojL8nHIRaw5lI/2yWnfcy8kaET098VRQW7jYWUlYIRER1QWDFFE9E0Jgf8Z1rD10Hr8dy0HpX6uOK+UWGNbVDRE926G3Tytu3UJEZAYYpIjqSa66GOuPXMBPh7KRda1Id9zXzRZje7XDyB5tYG+tkLBCIiKqbwxSRHVQrtEi/tRVrDmUjR0ncqH5a9HMlio5hvfwwJienghoY89Vx4mIzBSDFJGBhBA4e/UWfkm5hHWHLyBH/feimcFejn8tmukOayX/ehERmTv+S0+kh9ulGuzPuIYdJ3Kx82QuLty4rftZKxslngxsg4iennjAlYtmEhE1JwxSRPdw/loRdp6sCE6JZ6+h5K+bxgFAaWmBfh2d8HSwJx7xc+GimUREzRSDFNFfSso1OJx1Q3fVKeNqYZWfe9hbYbCvCwZ3dkG/jk786o6IiBikqHm7nH8bu05exc4Tudh3Jg+FpRrdzywtZAjxcsRDvi4Y7OuCB1xa8qZxIiKqgkGKmpVyjRZHzt+s+MruRC5O5BRU+blzSxUGd26Nwb4u+L8HnGFnxeUKiIjo3hikyOzl3SpB/Mmr2HkyF7tPXYW6uFz3M5kM6OHpgIc6V1x18nO340KZRESkNwYpMjtarcCfF/Ox80Qudp3MxZ8X8yHE3z93sFZgUKfWGNzZBQM7tUYrG6V0xRIRUZPGIEVmIb+oDLtPV1x1ij95FdcKS6v8vGsbOwzu7IIHO7ugh6cDLHnViYiI6gGDFDVJQgikXy7AzpMVV52Szt2A9o6rTi1Vcgx4wBmDO7tgUOfWcOXGwEREZAIMUtTolWu0yMwrRNplNdIuqXH8khppl9W4ftdVp06uLXVXnYK9HKGUW0hUMRERNRcMUtSoFJaU40ROgS40pV3Kx4mcgiqLYVayUligfwdnPOjrggc7tYZnK2sJKiYiouaMQYokk1tQXBGWKkPTZTUy8wqr3BheyVppiS7udvBzt4Ofhx38PezQydUWVgquKE5ERNJhkCKT02oFsq4V6r6SqwxNVwtKamzvYquCn8edockeXq2suSwBERE1OgxSVK+KyzQ4mVPwV2jKR9olNU7kFKDojhXDK8lkQHtnG/h52OtCk5+7HVrbqiSonIiIyHAMUmS064Wlf938na+70nT26q0qs+cqWSks4OtmV+VKk6+bLferIyKiJo2fYlRFcZkGebdKcO1W6d//W3jXn2+V/PUorfEcrWyU8Pf4OzT5e9jB28kGckvOoiMiIvMieZBaunQpPv74Y1y+fBn+/v5YtGgRBgwYcM/28fHxmDp1Ko4fPw4PDw+89dZbmDhxYpU2sbGxeO+993D27Fl06NABc+bMwahRowzqVwiBDz/8EF9++SVu3LiB3r17Y8mSJfD396/fF8DEtFqB/NtluuBzrbAEeQUluFZYWvHnv0LRtcJSXLtVilsl5fc/6R28nax19zFVXmlysVVxc18iImoWJA1Sa9euRVRUFJYuXYr+/ftjxYoVGDp0KNLS0tCuXbtq7TMzMzFs2DCMHz8e33//Pfbt24fXXnsNrVu3xujRowEAiYmJiIiIwOzZszFq1Chs3LgR4eHh2Lt3L3r37q13vx999BEWLlyI//73v+jUqRP+/e9/49FHH8XJkydha2vbcC9SDUrKNVVCUMV/V14xqghFV/8KS9cLS6Gp6bu2WigtLeDcUgmnlio4tVTCufJ/bVRwtlXCyUaFVjZKeDvboKVK8ixOREQkGZkQNU02bxi9e/dGUFAQli1bpjvWpUsXjBw5EvPmzavWfvr06di0aRPS09N1xyZOnIjU1FQkJiYCACIiIqBWq7FlyxZdmyFDhsDR0RGrV6/Wq18hBDw8PBAVFYXp06cDAEpKSuDq6ooFCxZgwoQJNY6npKQEJSV/z0RTq9Xw9PREfn4+7OzsjHmJavTptpP4YscZg55j30KhC0fOleHIpjIoVYalij/bquS8okRERM2WWq2Gvb29Xp/fkl1OKC0tRVJSEqKjo6scDw0NRUJCQo3PSUxMRGhoaJVjYWFhiImJQVlZGRQKBRITEzFlypRqbRYtWqR3v5mZmcjJyanSl0qlwqBBg5CQkHDPIDVv3jx8+OGH9x98HTnZKKG0tIBTS+XfV4xs7ghIdwUmR2slV/kmIiIyAcmCVF5eHjQaDVxdXascd3V1RU5OTo3PycnJqbF9eXk58vLy4O7ufs82lefUp9/K/62pzblz5+45phkzZmDq1Km6P1dekapvz/Xxwrh+3rxqREREJDHJb3C5OwwIIWoNCDW1v/u4PuesrzZ3UqlUUKlMvwYSZ78RERE1DpJ9Ijs7O8PS0rLa1afc3NxqV4Iqubm51dheLpfDycmp1jaV59SnXzc3NwAwqDYiIiJqfiQLUkqlEsHBwYiLi6tyPC4uDv369avxOX379q3Wftu2bQgJCYFCoai1TeU59enXx8cHbm5uVdqUlpYiPj7+nrURERFRMyQktGbNGqFQKERMTIxIS0sTUVFRwsbGRmRlZQkhhIiOjhaRkZG69hkZGcLa2lpMmTJFpKWliZiYGKFQKMT69et1bfbt2ycsLS3F/PnzRXp6upg/f76Qy+Vi//79evcrhBDz588X9vb2YsOGDeLo0aNi7Nixwt3dXajVar3Hl5+fLwCI/Pz8urxMRERE1IAM+fyWNEgJIcSSJUuEl5eXUCqVIigoSMTHx+t+Nm7cODFo0KAq7Xft2iUCAwOFUqkU3t7eYtmyZdXOuW7dOtG5c2ehUCiEr6+viI2NNahfIYTQarVi5syZws3NTahUKjFw4EBx9OhRg8bGIEVERNT0GPL5Lek6UubOkHUoiIiIqHEw5POb07+IiIiIjMQgRURERGQkBikiIiIiIzFIERERERmJQYqIiIjISAxSREREREZikCIiIiIyEoMUERERkZHkUhdgzirXOlWr1RJXQkRERPqq/NzWZ81yBikTKigoAAB4enpKXAkREREZqqCgAPb29rW24RYxJqTVanHp0iXY2tpCJpPV2latVsPT0xPZ2dlmvZ0Mx2k+msMYAY7T3HCc5sOUYxRCoKCgAB4eHrCwqP0uKF6RMiELCwu0bdvWoOfY2dmZ7S/9nThO89EcxghwnOaG4zQfphrj/a5EVeLN5kRERERGYpAiIiIiMhKDVCOhUqkwc+ZMqFQqqUsxKY7TfDSHMQIcp7nhOM1HYxkjbzYnIiIiMhKvSBEREREZiUGKiIiIyEgMUkRERERGYpAiIiIiMhKDVD1ZunQpfHx8YGVlheDgYOzZs6fW9vHx8QgODoaVlRXat2+P5cuXV2sTGxsLPz8/qFQq+Pn5YePGjXXut66kGOcHH3wAmUxW5eHm5lav47pbfY/z+PHjGD16NLy9vSGTybBo0aJ66bcupBijObyXK1euxIABA+Do6AhHR0c88sgjOHjwYJ37rSspxmkO7+eGDRsQEhICBwcH2NjYoEePHvjuu+/q3G9dSTHOhn4/TfF5UmnNmjWQyWQYOXJknfu9L0F1tmbNGqFQKMTKlStFWlqamDx5srCxsRHnzp2rsX1GRoawtrYWkydPFmlpaWLlypVCoVCI9evX69okJCQIS0tLMXfuXJGeni7mzp0r5HK52L9/v9H9NtVxzpw5U/j7+4vLly/rHrm5uSYZo6nGefDgQfHGG2+I1atXCzc3N/HZZ5/Vud+mOEZzeC+feeYZsWTJEpGcnCzS09PFP/7xD2Fvby8uXLhgdL9NdZzm8H7u3LlTbNiwQaSlpYkzZ86IRYsWCUtLS/H7778b3W9THWdDvp+mGGOlrKws0aZNGzFgwAAxYsSIOvWrDwapetCrVy8xceLEKsd8fX1FdHR0je3feust4evrW+XYhAkTRJ8+fXR/Dg8PF0OGDKnSJiwsTIwZM8bofutKqnHOnDlTdO/evY7V688U47yTl5dXjSGjId9PqcZobu+lEEKUl5cLW1tb8e233xrdb11JNU5zfD+FECIwMFC8++67RvdbV1KNsyHfT1ONsby8XPTv31989dVXYty4cdWClCneS361V0elpaVISkpCaGholeOhoaFISEio8TmJiYnV2oeFheHw4cMoKyurtU3lOY3pty6kGmel06dPw8PDAz4+PhgzZgwyMjLqOqQamWqcpujXWFKNsZK5vZdFRUUoKytDq1atjO63LqQaZyVzej+FEPjjjz9w8uRJDBw40Oh+60KqcVZqiPfTlGOcNWsWWrdujZdeeqle+tUHg1Qd5eXlQaPRwNXVtcpxV1dX5OTk1PicnJycGtuXl5cjLy+v1jaV5zSm37qQapwA0Lt3b6xatQpbt27FypUrkZOTg379+uHatWv1MbQqTDVOU/RrLKnGCJjnexkdHY02bdrgkUceMbrfupBqnID5vJ/5+flo2bIllEolHnvsMXzxxRd49NFHje63LqQaJ9Bw76epxrhv3z7ExMRg5cqV9davPuRGP5OqkMlkVf4shKh27H7t7z6uzzkN7beupBjn0KFDdf8dEBCAvn37okOHDvj2228xdepUwwehB1OM0xT91oUUYzS39/Kjjz7C6tWrsWvXLlhZWdWp37qSYpzm8n7a2toiJSUFt27dwh9//IGpU6eiffv2ePDBB43ut66kGGdDv5/1OcaCggI899xzWLlyJZydneu13/thkKojZ2dnWFpaVkuzubm51VJvJTc3txrby+VyODk51dqm8pzG9FsXUo2zJjY2NggICMDp06eNGUqtTDVOU/RrLKnGWJOm/F5+8sknmDt3LrZv345u3brVqd+6kGqcNWmq76eFhQU6duwIAOjRowfS09Mxb948PPjgg2b1ftY2zpqY6v00xRiPHz+OrKwsPPHEE7qfa7VaAIBcLsfJkyfh6elpkveSX+3VkVKpRHBwMOLi4qocj4uLQ79+/Wp8Tt++fau137ZtG0JCQqBQKGptU3lOY/qtC6nGWZOSkhKkp6fD3d3dmKHUylTjNEW/xpJqjDVpqu/lxx9/jNmzZ+P3339HSEhInfutC6nGWZOm+n7eTQiBkpISo/utC6nGWRNTvZ+mGKOvry+OHj2KlJQU3WP48OEYPHgwUlJS4Onpabr30ujb1EmncjplTEyMSEtLE1FRUcLGxkZkZWUJIYSIjo4WkZGRuvaV0zinTJki0tLSRExMTLVpnPv27ROWlpZi/vz5Ij09XcyfP/+eyx/cq19zGee0adPErl27REZGhti/f794/PHHha2tbZMaZ0lJiUhOThbJycnC3d1dvPHGGyI5OVmcPn1a737NYYzm8F4uWLBAKJVKsX79+irTxAsKCvTu11zGaQ7v59y5c8W2bdvE2bNnRXp6uvj000+FXC4XK1eu1LtfcxlnQ76fphjj3WqatWeK95JBqp4sWbJEeHl5CaVSKYKCgkR8fLzuZ+PGjRODBg2q0n7Xrl0iMDBQKJVK4e3tLZYtW1btnOvWrROdO3cWCoVC+Pr6itjYWIP6NQUpxhkRESHc3d2FQqEQHh4e4sknnxTHjx83yfgq1fc4MzMzBYBqj7vP05DvpxRjNIf30svLq8Zxzpw5U+9+TUGKcZrD+/nOO++Ijh07CisrK+Ho6Cj69u0r1qxZY1C/piDFOBv6/TTF58mdagpS9+vXGDIh/rpbi4iIiIgMwnukiIiIiIzEIEVERERkJAYpIiIiIiMxSBEREREZiUGKiIiIyEgMUkRERERGYpAiIiIiMhKDFBEREZGRGKSIiIiIjMQgRURERGQkBikiIiIiIzFIEREZaP369QgICECLFi3g5OSERx55BIWFhVKXRUQSkEtdABFRU3L58mWMHTsWH330EUaNGoWCggLs2bMH3P+dqHmSCf7tJyLS25EjRxAcHIysrCx4eXlJXQ4RSYxf7RERGaB79+54+OGHERAQgKeffhorV67EjRs3pC6LiCTCK1JERAYSQiAhIQHbtm3Dxo0bkZOTgwMHDsDHx0fq0oiogTFIERHVgUajgZeXF6ZOnYqpU6dKXQ4RNTDebE5EZIADBw7gjz/+QGhoKFxcXHDgwAFcvXoVXbp0kbo0IpIAgxQRkQHs7Oywe/duLFq0CGq1Gl5eXvj0008xdOhQqUsjIgnwqz0iIiIiI3HWHhEREZGRGKSIiIiIjMQgRURERGQkBikiIiIiIzFIERERERmJQYqIiIjISAxSREREREZikCIiIiIyEoMUERERkZEYpIiIiIiMxCBFREREZKT/Bwf6wpdon1dQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "s_grid = np.linspace(0.001, 0.04, 15)\n",
    "error1_s_list = []\n",
    "for s in s_grid:\n",
    "    DME_s = DME(s, rho, sigma)\n",
    "    DME_goal = evolve(s,rho) @ sigma @ np.conjugate(evolve(s, rho)).T\n",
    "    error1_s = np.linalg.norm(np.abs(DME_goal - DME_s))\n",
    "    error1_s_list.append(error1_s)\n",
    "plt.plot(s_grid, error1_s_list, label=r'$|E_s^\\rho\\sigma (E_s^\\rho)^\\dag- \\varepsilon_s(\\sigma)|$')\n",
    "plt.xlabel('s')\n",
    "plt.ylabel('op_norm error')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.19209592e-01+0.j         -3.14081258e-05+0.00239194j]\n",
      " [-3.14081258e-05-0.00239194j  8.80790408e-01+0.j        ]]\n",
      "[[ 1.19210435e-01+2.71050543e-20j -7.51335212e-06+2.39207278e-03j]\n",
      " [-7.51335212e-06-2.39207278e-03j  8.80789565e-01+0.00000000e+00j]]\n"
     ]
    }
   ],
   "source": [
    "trotter_step = 1\n",
    "s = 0.01\n",
    "DME_goal = evolve(s,rho) @ sigma @ np.conjugate(evolve(s, rho)).T\n",
    "# rho: quantum instruction, sigma: quantum register\n",
    "def DME_trotterize(s, rho, sigma, steps):\n",
    "    result = sigma\n",
    "    for M in range(steps):\n",
    "        result = DME(s/steps, rho, result)\n",
    "    return result\n",
    "DME_trotter = DME_trotterize(s, rho, sigma, trotter_step)\n",
    "print(DME_trotter)\n",
    "print(DME_goal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Trace norm')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgaElEQVR4nO3deVxU9f4/8NfMMAvrICKbslauqCCYopEtN7zaLS0rbdE2K1quKXVzaTP7GdX1er1miplmttq9Wl9LKqlcEzcE3EgtQVBBBGUGZRmY+fz+QCZHhmXYzgzzej4e87hw5nPOeX9m5M6rz/nM58iEEAJERERETkQudQFEREREnY0BiIiIiJwOAxARERE5HQYgIiIicjoMQEREROR0GICIiIjI6TAAERERkdNxkboAe2QymXDmzBl4enpCJpNJXQ4RERG1gBAC5eXlCAoKglze9BgPA5AVZ86cQXBwsNRlEBERUSsUFBSgV69eTbZhALLC09MTQN0L6OXlJXE1RERE1BJ6vR7BwcHmz/GmMABZUX/Zy8vLiwGIiIjIwbRk+gonQRMREZHTYQAiIiIip8MARERERE6HAYiIiIicDgMQEREROR0GICIiInI6kgegpUuXIjw8HBqNBjExMdi+fXuT7bdu3YqYmBhoNBpEREQgJSWlQZtFixahT58+cHV1RXBwMGbMmIGqqqqO6gIRERE5GEkD0Nq1azF9+nS8/PLLyMzMRHx8PMaMGYP8/Hyr7XNzczF27FjEx8cjMzMTc+bMwbRp07Bu3Tpzm88++wyzZs3C66+/jpycHKxcuRJr167F7NmzO6tbREREZOdkQggh1cmHDRuGIUOGYNmyZeZt/fr1w/jx45GcnNyg/cyZM7Fhwwbk5OSYtyUmJiI7Oxvp6ekAgOeeew45OTn4+eefzW1eeOEF7Nmzp9nRpXp6vR5arRY6nY4LIRIRETkIWz6/JRsBMhgMyMjIQEJCgsX2hIQE7Ny50+o+6enpDdqPHj0a+/btQ01NDQDghhtuQEZGBvbs2QMAOHHiBFJTU3H77bc3Wkt1dTX0er3Fg4iIiLouyW6FUVJSAqPRCH9/f4vt/v7+KCoqsrpPUVGR1fa1tbUoKSlBYGAgJk2ahHPnzuGGG26AEAK1tbV4+umnMWvWrEZrSU5OxhtvvNH2ThEREZFDkHwS9NX36xBCNHkPD2vtr9y+ZcsWzJ8/H0uXLsX+/fuxfv16fPfdd3jzzTcbPebs2bOh0+nMj4KCgtZ2h4iIiByAZCNAvr6+UCgUDUZ7iouLG4zy1AsICLDa3sXFBd27dwcAvPrqq5g8eTKmTp0KABg4cCAuXbqEJ598Ei+//DLk8oaZT61WQ61Wt0e3miYEzpeVofSSAdf5eXT8+YiIiOyZ0g1owY1LO4JkAUilUiEmJgZpaWm46667zNvT0tIwbtw4q/vExcXh22+/tdi2adMmxMbGQqlUAgAqKioahByFQgEhBCSc7w0A+OVgHm5ZHwUfSasgIiKyE3POACp3SU4t6SWwpKQkfPjhh1i1ahVycnIwY8YM5OfnIzExEUDdpakpU6aY2ycmJuLkyZNISkpCTk4OVq1ahZUrV+LFF180t7njjjuwbNkyfPnll8jNzUVaWhpeffVV3HnnnVAoFJ3exyv1CfCU9PxERERUR7IRIACYOHEiSktLMW/ePBQWFiIyMhKpqakIDQ0FABQWFlqsCRQeHo7U1FTMmDED77//PoKCgrB48WJMmDDB3OaVV16BTCbDK6+8gtOnT6NHjx644447MH/+/E7v39WCfH1wg+IzlF6qxmePD8OQ0G5Sl0RERCQdpZtkp5Z0HSB71ZHrAD2xZh/SjpzFK7f3w9T4iHY9NhERkTNziHWAnFVUsDcAILOgTNI6iIiInBkDUCeLvhyAsvLLJK2DiIjImTEAdbKBvbSQyYDTZZUoLucNWomIiKTAANTJPDVK9Par+zYYR4GIiIikwQAkgfp5QFmcB0RERCQJBiAJRIV4A2AAIiIikgoDkATqR4CyC8pgNHEVAiIios7GACSB3v6ecFMpcMlgxO/FF6Uuh4iIyOkwAElAIZdhUC8tACCr4ILE1RARETkfBiCJRAXX3QaD84CIiIg6HwOQRKIvT4TO5FfhiYiIOh0DkETqV4Q+drYcl6prpS2GiIjIyTAAScTPS4MgrQYmARw4pZO6HCIiIqfCACQhrgdEREQkDQYgCUVfngidmc9vghEREXUmBiAJXTkCJAQXRCQiIuosDEASigzSQiGXobi8GoU63hmeiIioszAASchVpUDfgLo7w/Pr8ERERJ2HAUhi0ebLYJwHRERE1FkYgCTGFaGJiIg6HwOQxOrvDH/wtA41RpO0xRARETkJBiCJRfi6w1PjgqoaE44WlUtdDhERkVNgAJKYXC4zjwJl8jIYERFRp2AAsgP19wXL4jfBiIiIOgUDkB2I4jfBiIiIOhUDkB0Y3MsbAPDHuUvQVdRIWwwREZETYACyA9091Ajt7gYAyD5VJm0xREREToAByE7UT4TmekBEREQdjwHITjAAERERdR4GIDsRHVK3InRm/gXeGZ6IiKiDMQDZiX6BnlAp5LhQUYP88xVSl0NERNSlMQDZCbWLAv2DvADwMhgREVFHYwCyI+YVobkgIhERUYdiALIj0ZcXROQtMYiIiDqW5AFo6dKlCA8Ph0ajQUxMDLZv395k+61btyImJgYajQYRERFISUmxeP6mm26CTCZr8Lj99ts7shvtIjq4biJ0zhk9qmuNEldDRETUdUkagNauXYvp06fj5ZdfRmZmJuLj4zFmzBjk5+dbbZ+bm4uxY8ciPj4emZmZmDNnDqZNm4Z169aZ26xfvx6FhYXmx6FDh6BQKHDvvfd2VrdaLdjHFT7uKhiMJhw5o5e6HCIioi5L0gC0cOFCPP7445g6dSr69euHRYsWITg4GMuWLbPaPiUlBSEhIVi0aBH69euHqVOn4rHHHsOCBQvMbXx8fBAQEGB+pKWlwc3NrckAVF1dDb1eb/GQgkwm43pAREREnUCyAGQwGJCRkYGEhASL7QkJCdi5c6fVfdLT0xu0Hz16NPbt24eaGuv30Fq5ciUmTZoEd3f3RmtJTk6GVqs1P4KDg23sTfuJ5kRoIiKiDidZACopKYHRaIS/v7/Fdn9/fxQVFVndp6ioyGr72tpalJSUNGi/Z88eHDp0CFOnTm2yltmzZ0On05kfBQUFNvam/fx5Z/gyyWogIiLq6lykLkAmk1n8LoRosK259ta2A3WjP5GRkbj++uubrEGtVkOtVre05A416PKd4fPPV6D0YjW6e9hHXURERF2JZCNAvr6+UCgUDUZ7iouLG4zy1AsICLDa3sXFBd27d7fYXlFRgS+//LLZ0R97o3VV4poedZfrOApERETUMSQLQCqVCjExMUhLS7PYnpaWhhEjRljdJy4urkH7TZs2ITY2Fkql0mL7V199herqajz00EPtW3gnqL8vGAMQERFRx5D0W2BJSUn48MMPsWrVKuTk5GDGjBnIz89HYmIigLq5OVOmTDG3T0xMxMmTJ5GUlIScnBysWrUKK1euxIsvvtjg2CtXrsT48eMbjAw5An4TjIiIqGNJOgdo4sSJKC0txbx581BYWIjIyEikpqYiNDQUAFBYWGixJlB4eDhSU1MxY8YMvP/++wgKCsLixYsxYcIEi+MeO3YMO3bswKZNmzq1P+3lygBkMgnI5Y3PiSIiIiLbyUT9LGIy0+v10Gq10Ol08PLy6vTz1xpNiJz7I6pqTPgp6UZc6+fZ6TUQERE5Gls+vyW/FQY15KKQY1BPbwBcD4iIiKgjMADZKa4HRERE1HEYgOwUJ0ITERF1HAYgOxV9eQTot6JyVBp4Z3giIqL2xABkpwK1rvD3UsNoEjh4Wid1OURERF0KA5Ad+/My2AVpCyEiIupiGIDsWFQwV4QmIiLqCAxAdqx+HhC/Ck9ERNS+GIDs2MCeWshlQKGuCmf1VVKXQ0RE1GUwANkxd7ULevvXrQLNUSAiIqL2wwBk56K5ICIREVG7YwCyc9GXJ0Jn5vObYERERO2FAcjO1d8S4+BpHYwm3reWiIioPTAA2blrenjAQ+2CCoMRx86WS10OERFRl8AAZOcUchkG9dIC4ERoIiKi9sIA5AD+nAjNeUBERETtgQHIAXBFaCIiovbFAOQA6u8Jdrz4IsqraqQthoiIqAtgAHIAPTzV6OntCiGAA6d4Z3giIqK2YgByEFwQkYiIqP0wADmI+stg/CYYERFR2zEAOYgrR4CE4IKIREREbcEA5CAGBGnhIpeh5GI1Tl2olLocIiIih8YA5CA0SgX6B3kB4DwgIiKitmIAciD184AYgIiIiNqGAciBMAARERG1DwYgBxIdUrci9MHTOhhqTRJXQ0RE5LgYgBxIWHc3eLspYag14bcivdTlEBEROSwGIAcik8kwuJc3AF4GIyIiagsGIAdjngfEBRGJiIhajQHIwdQviJjJESAiIqJWYwByMPUjQLkll1BWYZC2GCIiIgfFAORgvN1UCPd1B8B5QERERK3FAOSAeGNUIiKitpE8AC1duhTh4eHQaDSIiYnB9u3bm2y/detWxMTEQKPRICIiAikpKQ3alJWV4dlnn0VgYCA0Gg369euH1NTUjupCp7vyxqhERERkO0kD0Nq1azF9+nS8/PLLyMzMRHx8PMaMGYP8/Hyr7XNzczF27FjEx8cjMzMTc+bMwbRp07Bu3TpzG4PBgNtuuw15eXn43//+h6NHj2LFihXo2bNnZ3Wrw9WPAGWf4p3hiYiIWkMmJPwEHTZsGIYMGYJly5aZt/Xr1w/jx49HcnJyg/YzZ87Ehg0bkJOTY96WmJiI7OxspKenAwBSUlLwz3/+E7/99huUSmWL6qiurkZ1dbX5d71ej+DgYOh0Onh5ebW2ex3GUGtC5NwfYag1YfOLN5nnBBERETkzvV4PrVbbos9vyUaADAYDMjIykJCQYLE9ISEBO3futLpPenp6g/ajR4/Gvn37UFNTAwDYsGED4uLi8Oyzz8Lf3x+RkZF46623YDQaG60lOTkZWq3W/AgODm5j7zqWykWOyMt3hs/MvyBxNURERI5HsgBUUlICo9EIf39/i+3+/v4oKiqyuk9RUZHV9rW1tSgpKQEAnDhxAv/73/9gNBqRmpqKV155Bf/6178wf/78RmuZPXs2dDqd+VFQUNDG3nW8+vuCcR4QERGR7VykLkAmk1n8LoRosK259lduN5lM8PPzwwcffACFQoGYmBicOXMG//znP/Haa69ZPaZarYZarW5LNzod7wxPRETUepIFIF9fXygUigajPcXFxQ1GeeoFBARYbe/i4oLu3bsDAAIDA6FUKqFQKMxt+vXrh6KiIhgMBqhUqnbuiTTqA9CRM3pU1RihUSqa3oGIiIjMJLsEplKpEBMTg7S0NIvtaWlpGDFihNV94uLiGrTftGkTYmNjzROeR44cid9//x0mk8nc5tixYwgMDOwy4QcAenVzha+HCrUmgcNndFKXQ0RE5FAk/Rp8UlISPvzwQ6xatQo5OTmYMWMG8vPzkZiYCKBubs6UKVPM7RMTE3Hy5EkkJSUhJycHq1atwsqVK/Hiiy+a2zz99NMoLS3F888/j2PHjmHjxo1466238Oyzz3Z6/zqSTCZDVHDdPCAuiEhERGQbSecATZw4EaWlpZg3bx4KCwsRGRmJ1NRUhIaGAgAKCwst1gQKDw9HamoqZsyYgffffx9BQUFYvHgxJkyYYG4THByMTZs2YcaMGRg0aBB69uyJ559/HjNnzuz0/nW06BBv/JRzlvOAiIiIbCTpOkD2ypZ1BKT06+8lePDD3ejVzRU7Zt4idTlERESScoh1gKjtBvXSQiYDTl2oxLny6uZ3ICIiIgAMQA7NU6PEdX4eAPh1eCIiIlswADm4P9cD4orQRERELcUA5ODqvwnGESAiIqKWYwBycNEh3gCA7AIdjCbOZyciImoJBiAH19vfE24qBS5W1+KPcxelLoeIiMghMAA5OIVchoE9tQCALC6ISERE1CIMQF1A1OXLYJmcCE1ERNQiDEBdQDRviUFERGQTBqAuoH4i9LGz5bhUXSttMURERA6AAagL8PfSIFCrgUkAB0/zzvBERETNYQDqIuoXRORlMCIiouYxAHUR9ZfBuCI0ERFR8xiAugiuCE1ERNRyDEBdxMCeWijkMpzVV6NQVyl1OURERHaNAaiLcFUp0MffEwDnARERETWHAagL+XMeUJmkdRAREdk7BqAupP6bYLwlBhERUdMYgLqQ+hGgA6fLUGs0SVsMERGRHWMA6kIifD3gqXFBVY0JvxWVS10OERGR3WIA6kLkctmfl8E4D4iIiKhRDEBdDAMQERFR8xiAuhgGICIiouYxAHUx9QHo9+KL0FXWSFsMERGRnWIA6mK6e6gR4uMGADhwqkzaYoiIiOwUA1AXxPWAiIiImsYA1AXVB6BMzgMiIiKyigGoC7rylhhCCGmLISIiskMMQF1Q/yAvqBRynL9kQMF53hmeiIjoagxAXZDaRYF+QV4AgMyCCxJXQ0REZH8YgLqo6Pp5QJwITURE1AADUBd15TwgIiIissQA1EXVfxPsyBk9qmuN0hZDRERkZxiAuqgQHzf4uKtgMJqQU8g7wxMREV1J8gC0dOlShIeHQ6PRICYmBtu3b2+y/datWxETEwONRoOIiAikpKRYPL969WrIZLIGj6qqqo7sht2RyWQY3EsLAMjM50RoIiKiK0kagNauXYvp06fj5ZdfRmZmJuLj4zFmzBjk5+dbbZ+bm4uxY8ciPj4emZmZmDNnDqZNm4Z169ZZtPPy8kJhYaHFQ6PRdEaX7Ep0SDcAnAdERER0NRcpT75w4UI8/vjjmDp1KgBg0aJF+PHHH7Fs2TIkJyc3aJ+SkoKQkBAsWrQIANCvXz/s27cPCxYswIQJE8ztZDIZAgICWlxHdXU1qqurzb/r9fpW9si+8M7wRERE1kk2AmQwGJCRkYGEhASL7QkJCdi5c6fVfdLT0xu0Hz16NPbt24eamj/vfH7x4kWEhoaiV69e+Nvf/obMzMwma0lOToZWqzU/goODW9kr+zL4cgA6WVqB85cM0hZDRERkRyQLQCUlJTAajfD397fY7u/vj6KiIqv7FBUVWW1fW1uLkpISAEDfvn2xevVqbNiwAV988QU0Gg1GjhyJ48ePN1rL7NmzodPpzI+CgoI29s4+aF2VuKaHOwAgiwsiEhERmbXqEtiePXuwZcsWFBcXw2QyWTy3cOFCm44lk8ksfhdCNNjWXPsrtw8fPhzDhw83Pz9y5EgMGTIE7733HhYvXmz1mGq1Gmq12qa6HUVUcDf8ce4SsvLLcEtf/+Z3ICIicgI2B6C33noLr7zyCvr06QN/f3+LQNJUcLmar68vFApFg9Ge4uLiBqM89QICAqy2d3FxQffu3a3uI5fLMXTo0CZHgLqyqBBvrNt/ineGJyIiuoLNAeg///kPVq1ahUceeaRNJ1apVIiJiUFaWhruuusu8/a0tDSMGzfO6j5xcXH49ttvLbZt2rQJsbGxUCqVVvcRQiArKwsDBw5sU72Oqv6WGNkFZTCZBOTylodUIiKirsrmACSXyzFy5Mh2OXlSUhImT56M2NhYxMXF4YMPPkB+fj4SExMB1M3NOX36NNasWQMASExMxJIlS5CUlIQnnngC6enpWLlyJb744gvzMd944w0MHz4c1113HfR6PRYvXoysrCy8//777VKzo+kT4AmNUg59VS1OlFzCtX4eUpdERK1gNBotvuxB5KxUKhXk8rZPYbY5AM2YMQPvv/+++avobTFx4kSUlpZi3rx5KCwsRGRkJFJTUxEaGgoAKCwstFgTKDw8HKmpqeYagoKCsHjxYouvwJeVleHJJ59EUVERtFotoqOjsW3bNlx//fVtrtcRKRVyDOypxd68C8gqKGMAInIwQggUFRWhrKxM6lKI7IJcLkd4eDhUKlWbjiMT9bOIW8hkMuH222/HsWPH0L9//waXntavX9+mguyBXq+HVquFTqeDl5eX1OW02fyNR7Biey4eGh6C/zfeOS8FEjmqwsJClJWVwc/PD25ubjbNtSTqakwmE86cOQOlUomQkJAGfw+2fH7bPAL097//HZs3b8bNN9+M7t2784/RAUQFdwOQi8z8MqlLISIbGI1Gc/hp7IseRM6mR48eOHPmDGpraxud/9sSNgegNWvWYN26dbj99ttbfVLqXNEh3gCA34rKUWkwwlWlkLYgImqR+jk/bm5uEldCZD/qL30ZjcY2BSCbZxH5+PjgmmuuafUJqfMFajXw81TDaBI4dEYndTlEZCOOtBP9qb3+HmwOQHPnzsXrr7+OioqKdimAOp5MJvvzvmC8DEZERGT7JbDFixfjjz/+gL+/P8LCwhoMP+3fv7/diqP2ExXijU1HziKTt8QgIiKyPQCNHz++A8qgjhYd3A0AR4CIiIgAGwNQbW0tAOCxxx7rMndMdxaDemkhlwFndFUo1lfBz0sjdUlERF1CWFgY8vLypC7DYfzrX//Ciy++iN9++w19+vQBADz55JP48MMPUV5eDnd3906pw6Y5QC4uLliwYAGMRmNH1UMdxF3tgt7+ngDA+4IREZFkDh06hIEDB+Lo0aMAgNzcXOzbtw8RERGdFn6AVkyCvvXWW7Fly5YOKIU6Wv1EaK4HRESdYcuWLQgLC2vVvu+++y4GDBiAsLAwbNiwwaZ9S0tL4efn16pRmcbOe88992DhwoUWbRMTExEVFYUzZ84gKioKr732ms3na0pb+tHZ57X2+jS2/dChQ7jvvvvMAejNN9/E+PHjMWDAAHObtvzbaSmbA9CYMWMwe/ZsvPjii/jiiy+wYcMGiwfZr/r1gLI4EZqIJHLjjTdCJpM1eDz44IPmNp988glycnJw8OBBfPXVV5g3b55N50hOTsYdd9zR4AO0uXM3dd7XXnsN8+fPh16vN29LSUlBVlYWgoKCkJWVZXOdre1HR2vNea29Pta2CyGQl5eHv/3tb/jtt99w/PhxnD59Gt27d0dkZGR7dqNZNk+CfvrppwHAatKTyWS8PGbHoi5PhD54SgejSUDBO8MTUScSQiArKwsLFiywCDwA4OHx530KlyxZgm+++QZyuRwDBgxo8KHalMrKSqxcuRKpqak2n7up8w4aNAhhYWH47LPPzJ+DHamxftjreRt7fa7efuLECYSEhKBfv374448/MG/ePLz22mv49NNPceONN7Z3d5pk8wiQyWRq9MHwY9+u9fOAu0qBSwYjjheXS10OETmZ48ePo7y8HDfeeCMCAgIsHvUhpKqqCqdOnUJgYCAAIDs7G4MHDzYfY9euXbj11lvh6+vbYCSnrKwM33//PVxcXBAXF2fTuZs7LwDceeed+OKLLzryJTJrrB8pKSkYOHAgXF1dodVqccstt3TKeZt73YHGX58rtx86dAiRkZFQq9XQ6/UoLS3FyJEjzds7U9vvJ08OQyGXYVAvbwCcB0REnS8jIwMuLi4YNGhQo20OHTqEc+fO4ezZs6iursbcuXPx7LPPAqgLJTfddBMGDx6Mbdu24YcffoCPjw9uvvlmrF27Ft7e3ti2bRtiY2NtPndT5613/fXXY8+ePaiurm7Dq9Ay1vqxbt06zJo1C6+++iqOHj2KnTt34oUXXujw87bkdQcaf32u3H5l0Jk5cyYWLFgAADh27Jj5G2GdxeZLYACwdetWLFiwADk5OZDJZOjXrx/+8Y9/ID4+vr3ro3YWHeKN9BOlyMovw/3Xh0hdDhHZSAiByprOH213VSrafAuC/fv3w2g0Nrix6/33348VK1YAqPuwvf/++zF27FgYDAY899xzuOmmmwAA06ZNw7hx48xTMPr374/7778fu3fvxn333QcAyMvLQ1BQkM3nbuq89Xr27Inq6moUFRUhNDS0Ta9Fc6z149ixYwgJCUFCQoI5dFw5cbijztuS1x1o/PW5cvvhw4cxefJkAMDEiRMBAGfOnIGvr6/5Hl+dxeYA9Omnn+LRRx/F3XffjWnTpkEIgZ07d+LWW2/F6tWr8cADD3REndROzLfE4FfhiRxSZY0R/V/7sdPPe2TeaLipWvXfzGYZGRm49957MX/+fIvt3bp1M/+cnZ2Nu+++Gx9//LFFm7Nnz2LHjh345ZdfLLa7u7tbBLPKykpoNA3XOWvu3I2d90qurq4A0OBWUP3797fafu7cuXjjjTcaPR4A7N271+qIlbV+PPHEE/jqq6/g4+MDNzc3ZGdnt+jenPfccw/efvtt/N///V+z6+9cfd6Wvu5A46/Plds///zzBvUFBQXh8OHDzfajvdn8r3n+/Pl49913MWPGDPO2559/HgsXLsSbb77JAGTnoi5/E+xYcTkuVtfCQ922/0MjImqpzMxMzJs3D9dee22jbbKzs61e1snIyIDJZGowLycjI8MiQPj6+uLChYbfdG3u3I2d90rnz58HAPTo0cNie2MThp977jlMmjSpyWM29k2rq/tRU1ODSZMmYejQoVixYgW8vb0RERHR5LGBujum5+bm4tprr7VYf6dPnz5W19+5+rwtfd2Bxl+fxrZLzeZPvxMnTuCOO+5osP3OO+/EnDlz2qUo6jh+nhr09HbF6bJKHCgow4hrfaUuiYhs4KpU4Mi80ZKcty1OnDiBsrIyREdHN9lu69atVrebTCYAdSMj9Zd/Dh48iG3btll8/Tw6Ohqffvqpzedu7LxXOnToEHr16gVf35b9/6avr2+L217t6n58/fXX+P333/HTTz81us/y5cuxfPlyGAwGDBo0CJ9//jn27t2LoUOHmuu3tv5ORkZGo+dt6etef3xrr4+tr1tnsXkSdHBwMH7++ecG23/++WfeHsNB1I8CcUVoIscjk8ngpnLp9Edb5//Uf8j6+/ujqKjI4lH/IQvUfYgPGTIEkZGRFlcUhg0bBldXV7z00kv47bffsHHjRowbNw6JiYkYMWKEud3o0aNx+PDhBqMYzZ27sfNeafv27UhISGjT69BSV/fDYDCgsLAQn3zyCfLy8nDo0CEsX74cNTU1AIALFy7ggw8+wN69e3Ho0CEsXboUALBp0yYkJCS0eP2dq8/b0tcdaPz16czXzSbCRkuXLhUqlUokJiaKNWvWiE8++UQ89dRTQq1Wi5SUFFsPZ5d0Op0AIHQ6ndSldIgV2/4QoTO/E1M/3it1KUTUhMrKSnHkyBFRWVkpdSmtsnnzZhEaGiqEEGLWrFkCQIOHUqkUVVVVQgghzp8/L4YMGSJqa2uFEEJcuHDB4njffvut6N27t1AqlSIiIkK88847wmg0Njjv8OHDLT6Pmjt3c+cVou698PLyEunp6e3wyrTMlf2oqakRSUlJolevXkKpVAp/f3/xwAMPmNvq9XoRGhoqXnjhBXHo0CHz9lGjRomysjLx+++/iyFDhoiqqioxatQo8dBDD4kdO3aIxMRE8fnnnzd6XiFa9ro39vq09nW78t/O1Zr6u7Dl89vmACSEEOvXrxcjR44UPj4+wsfHR4wcOVJ88803rTmUXerqAWhvbqkInfmdiHkzTZhMJqnLIaJGdKUA1BKNfYjbauPGjaJfv35Ww1Frz7tkyRJx2223tbqm1mhNP9asWSP69u0rvv32W6HT6cSoUaOEEEJ88803YsqUKUIIIaKjo8WYMWOEEELccMMN4sCBA206rxCNvz6tfd06IwC1ah2gu+66Czt27EBpaSlKS0uxY8cOjBs3ru3DUdQpIntq4SKXoeRiNU6XVUpdDhERAMDT0xMHDx7E4MGDcc899+C7775r1XHGjh2Lp556CqdPn2638yqVSrz33nutqqe1bOnH8ePH4enpicmTJyM+Ph7V1dX4+eefcfPNNwOATevv2Pr6AY2/PlK8bi3V6q8AGQwGFBcXW1y7BYCQEK4tY+80SgX6BXrh4GkdsgrK0Kubm9QlERHh+PHjuO666zB58mRs3769TQsOPv/88+163ieffLLVtbRFS/vx5ptvYvfu3XBzc0N8fDzuuusuPPvss5gyZQoA2Lz+ji2vH9D46yPV69YSNgeg48eP47HHHsPOnTsttgsheC8wBxIV7F0XgPLL8LdBDRcNIyJqq7CwMEyfPr3F7a19iHcGqc7bntasWdNg2969e7FkyRIAsKv1d1rC1n87rSETQghbdhg5ciRcXFwwa9YsBAYGNvhmwNVrBTgivV4PrVYLnU4HLy8vqcvpEOsyTuGF/2YjJrQb1j09ovkdiKjTVVVVITc3F+Hh4VYX9yNyRk39Xdjy+W3zCFBWVhYyMjLQt29fW3clOxJ9+avwh07rUGM0QangbeGIiMh52Pyp179/f5SUlHRELdSJwn3doXVVorrWhN8KeWd4IiJyLjYHoHfeeQcvvfQStmzZgtLSUuj1eosHOQaZTIbB5vuCNVw2noiIqCuz+RLYX/7yFwDArbfearGdk6AdT1SwN7YdO4fM/DJMjpO6GiIios5jcwDavHlzR9RBEqifB8Q7wxMRkbOxOQCNGjWqI+ogCUT18gYAnCi5BF1FDbRuSmkLIiIi6iT86o8T6+auQlj3ukUQs06VSVsMERFRJ2IAcnJRlydCZ+ZzIjQRETkPBiAnFx3SDQDnARERkXNhAHJy9SNA2QVlsHFRcCIiIofVqgBUW1uLn376CcuXL0d5ed0iemfOnMHFixdtPtbSpUvNy1nHxMRg+/btTbbfunUrYmJioNFoEBERgZSUlEbbfvnll5DJZBg/frzNdTmLfoFeULnIcaGiBidLK6Quh4iIqFPYHIBOnjyJgQMHYty4cXj22Wdx7tw5AMC7776LF1980aZjrV27FtOnT8fLL7+MzMxMxMfHY8yYMcjPz7faPjc3F2PHjkV8fDwyMzMxZ84cTJs2DevWrbNa54svvoj4+Hhbu+hUVC5yRAbV3S8lkwsiEhGRk7A5AD3//POIjY3FhQsX4Orqat5+11134eeff7bpWAsXLsTjjz+OqVOnol+/fli0aBGCg4OxbNkyq+1TUlIQEhKCRYsWoV+/fpg6dSoee+wxLFiwwKKd0WjEgw8+iDfeeAMRERHN1lFdXe3UK1pHBV+eB5RfJm0hREQOKCwsTOoSqBVsDkA7duzAK6+8ApVKZbE9NDQUp0+fbvFxDAYDMjIykJCQYLE9ISEBO3futLpPenp6g/ajR4/Gvn37UFNTY942b9489OjRA48//niLaklOToZWqzU/goODW9yPriCKCyISEZGTsTkAmUwmq7e7OHXqFDw9PVt8nJKSEhiNRvj7+1ts9/f3R1FRkdV9ioqKrLavra0136D1119/xcqVK7FixYoW1zJ79mzodDrzo6CgoMX7dgXRlydCHynUo6qGtzIhovaxZcuWVo+OvPvuuxgwYADCwsKwYcMGm/YtLS2Fn58f8vLy2u2899xzDxYuXGjRNjExEVFRUThz5gyioqLw2muv2Xy+9taavlvrW1veO0dhcwC67bbbsGjRIvPvMpkMFy9exOuvv46xY8faXIBMJrP4vf6eYra0r99eXl6Ohx56CCtWrICvr2+La1Cr1fDy8rJ4OJNe3Vzh66FCjVHg8BnnuvxHRJ3rxhtvhEwma/B48MEHzW0++eQT5OTk4ODBg/jqq68wb948m86RnJyMO+64o8EHeHPnbuq8r732GubPn28xRSIlJQVZWVkICgpCVlaWzXV2hMb63hRrfXMGNt8K49///jduvvlm9O/fH1VVVXjggQdw/Phx+Pr64osvvmjxcXx9faFQKBqM9hQXFzcY5akXEBBgtb2Liwu6d++Ow4cPIy8vD3fccYf5eZPJBABwcXHB0aNHcc0117S4Rmchk8kQFeyNn3KKkVVQhpjQblKXRERdkBACWVlZWLBggUXgAQAPDw/zz0uWLME333wDuVyOAQMG2PTBXFlZiZUrVyI1NdXmczd13kGDBiEsLAyfffYZnn766RbX05ka63tzHKFvHcHmEaD6pPuPf/wDTz31FKKjo/H2228jMzMTfn5+LT6OSqVCTEwM0tLSLLanpaVhxIgRVveJi4tr0H7Tpk2IjY2FUqlE3759cfDgQWRlZZkfd955J26++WZkZWU53dweW9SvB8R5QETUUY4fP47y8nLceOONCAgIsHjUh5CqqiqcOnUKgYGBAIDs7GwMHjzYfIxdu3bh1ltvha+vb4ORnLKyMnz//fdwcXFBXFycTedu7rwAcOedd9r0H/ptlZKSgoEDB8LV1RVarRa33HJLk+0b63tzrxnQ+X2zBzaPAAGAq6srHn30UTz66KNtOnlSUhImT56M2NhYxMXF4YMPPkB+fj4SExMB1M3NOX36NNasWQOg7nrrkiVLkJSUhCeeeALp6elYuXKl+U3TaDSIjIy0OIe3tzcANNhOluq/CcZbYhDZOSGAGgnW7FK6AU1MT2iJjIwMuLi4YNCgQY22OXToEM6dO4ezZ8/C29sbc+fOxZw5cwDUhZKbbroJzzzzDN577z0UFBTggQcewODBg5GYmAhvb29s27YNsbGxNp+7qfPWu/7665GcnIzq6mqo1eo2vBLNW7duHWbNmoUPPvgAw4cPR3l5ebPzeqz1vSWvGdC5fbMXNgeg5ORk+Pv747HHHrPYvmrVKpw7dw4zZ85s8bEmTpyI0tJSzJs3D4WFhYiMjERqaipCQ0MBAIWFhRZrAoWHhyM1NRUzZszA+++/j6CgICxevBgTJkywtRt0lUHBWshkwKkLlSi5WA1fD+f4AyByODUVwFtBnX/eOWcAlXubDrF//34YjUZ0797dYvv9999v/uJKdnY27r//fowdOxYGgwHPPfccbrrpJgDAtGnTMG7cOPOE3f79++P+++/H7t27cd999wEA8vLyEBTU8PVp7txNnbdez549UV1djaKiIvPnVEc5duwYQkJCkJCQYA4pAwYMaHIfa31vyWsGdG7f7IXNAWj58uX4/PPPG2wfMGAAJk2aZFMAAoBnnnkGzzzzjNXnVq9e3WDbqFGjsH///hYf39oxqCEvjRLX9vDA8eKLyMovw1/6W5+HRUTUWhkZGbj33nsxf/58i+3duv057zA7Oxt33303Pv74Y4s2Z8+exY4dO/DLL79YbHd3d7f4ckxlZSU0Go3N527svFeqX/uuosJyBK5///5W28+dOxdvvPFGo8cDgL1791odsXriiSfw1VdfwcfHB25ubsjOzm52DuvVfW/pawY03reuzOYAVFRUZL5GeqUePXqgsLCwXYoiaUQFe9cFoAIGICK7pXSrG42R4rxtlJmZiXnz5uHaa69ttE12djZeeOGFBtszMjJgMpkazMvJyMiwCBC+vr64cKHhpfzmzt3Yea90/vx5AHWfd1dqbNLxc889h0mTJjV5TGvf1qqpqcGkSZMwdOhQrFixAt7e3i1a1Pfqvrf0NQMa71tXZnMACg4Oxq+//orw8HCL7b/++qvVYUdyHFEh3vhvxineEoPInslkbb4UJYUTJ06grKwM0dHRTbbbunWr1e313+itrKw0XxI6ePAgtm3bZvH18+joaHz66ac2n7ux817p0KFD6NWrV4uXWfH19bVpSZZ6X3/9NX7//Xf89NNPVp9fvnw5li9fDoPBgEGDBpmvylzd95a+ZoDtfesKbP4W2NSpUzF9+nR89NFHOHnyJE6ePIlVq1ZhxowZeOKJJzqiRuok0ZcnQh8o0MFk4p3hiaj9ZGRkAPhzsdsrH/Uf1EDdh/uQIUMQGRmJBx54wLx92LBhcHV1xUsvvYTffvsNGzduxLhx45CYmGjxzeHRo0fj8OHDDUZCmjt3Y+e90vbt2xvcjaAjGAwGFBYW4pNPPkFeXh4OHTqE5cuXo6amBhcuXMAHH3yAvXv34tChQ1i6dKl5v6v73tLXrDP7ZleEjUwmk3jppZeERqMRcrlcyOVy4ebmJt544w1bD2W3dDqdACB0Op3UpXSqmlqj6PvK9yJ05nfiWJFe6nKInF5lZaU4cuSIqKyslLqUVtm8ebMIDQ0VQggxa9YsAaDBQ6lUiqqqKiGEEOfPnxdDhgwRtbW1QgghLly4YHG8b7/9VvTu3VsolUoREREh3nnnHWE0Ghucd/jw4SIlJcX8e3Pnbu68QtS9F15eXiI9Pb0dXpmm1dTUiKSkJNGrVy+hVCqFv7+/eOCBB4QQQuj1ehEaGipeeOEFcejQoQb7Xt33lrxm1vp25Xtnb5r6u7Dl89vmAFSvvLxc7NmzRxw8eND8j7ercNYAJIQQ96bsFKEzvxNr9+ZLXQqR0+tKAaglmvtwb6mNGzeKfv36WQ1HrT3vkiVLxG233dbqmtqTXq8Xa9asEX379hXffvutxXO29l0I631zhgBk8yWweh4eHhg6dCgiIyOdZs0AZ1B/X7BM3hmeiDqZp6cnDh48iMGDB+Oee+7Bd99916rjjB07Fk899VSLb9DdkvMqlUq89957raqnPR0/fhyenp6YPHky4uPjUV1dbfG8rX0H7Kdvna1VCyHu3bsX//3vf5Gfnw+DwWDx3Pr169ulMJJGNO8MT0QSOX78OK677jpMnjwZ27dvb/Dhbovnn3++Xc/75JNPtrqW9vTmm29i9+7dcHNzQ3x8PO66664GbWzpO2A/fetsNo8Affnllxg5ciSOHDmCr7/+GjU1NThy5Ah++eUXaLXajqiROlH9itBHi/SoMNRKXA0RObKwsDBMnz69xe3ffPNN9OnTB9HR0dBoNFY/3DuCVOdtjTVr1uDo0aPIzMzE4sWLIZe3+kJOk2x97xyRTAhh09d9Bg0ahKeeegrPPvssPD09kZ2djfDwcDz11FMIDAxsdtEnR6DX66HVaqHT6ZzuzvAAMPytn1Gkr8LaJ4djWET35ncgog5RVVWF3NxchIeHW13cj8gZNfV3Ycvnt83R8Y8//sDtt98OAFCr1bh06RJkMhlmzJiBDz74wNbDkR2qvwyWyctgRETURdkcgHx8fFBeXg6g7t4hhw4dAgCUlZU51RLaXZn5zvCcCE1ERF2UzZOg4+PjkZaWhoEDB+K+++7D888/j19++QVpaWm49dZbO6JG6mTmAMQRICIi6qJsDkBLlixBVVUVAGD27NlQKpXYsWMH7r77brz66qvtXiB1voG9tFDIZSjSV6FQV4lAravUJREREbUrmy6B1dbW4ttvvzXPOpfL5XjppZewYcMGLFy40OKOvuS43FQu6OPvCYCXwYiIqGuyKQC5uLjg6aefbtPaDOQYorgeEJHdsPHLukRdWnv9Pdg8CXrYsGHIzMxsl5OT/aqfB8RvghFJR6lUAgC/YEJ0hfoFmBUKRZuOY/McoGeeeQYvvPACTp06hZiYGLi7u1s8P2jQoDYVRPah/pYYB0/pUGs0wUXRMYttEVHjFAoFvL29UVxcDABwc3ODTCaTuCoi6ZhMJpw7dw5ubm5wcWnVzSzMWrz3Y489hkWLFmHixIkAgGnTppmfk8lkEEJAJpPBaDS2qSCyD9f08ICn2gXl1bU4erYcA4K4yjeRFAICAgDAHIKInJ1cLkdISEib/2OgxQHo448/xttvv43c3Nw2nZAcg1wuw+Bgb+z4vQRZBWUMQEQSkclkCAwMhJ+fH2pqaqQuh0hyKpWqXW4B0uIAVD/pKDQ0tM0nJccQVR+A8svw4DC+70RSUigUbZ7zQER/silC8dqzc+FEaCIi6qpsmkHUu3fvZkPQ+fPn21QQ2Y/6r8L/ce4i9FU18NIopS2IiIiondgUgN544w1otZwL4ix8PdQI9nFFwflKHCjQ4YbrfKUuiYiIqF3YFIAmTZoEPz+/jqqF7FBUcDcUnK9EVsEFBiAiIuoyWjwHiPN/nJN5HhBviUFERF1IiwMQl2J3TtFX3BKD/waIiKiraPElMJPJ1JF1kJ3qH+gFpUKG0ksGnLpQiWAfN6lLIiIiajPe34CapFEq0D/QCwC/Dk9ERF0HAxA1KzqkGwAgM/+CxJUQERG1DwYgalb9ROgsjgAREVEXwQBEzaoPQIfP6GGo5VwwIiJyfAxA1KzQ7m7o5qaEodaEnEK91OUQERG1GQMQNUsmk12xHhDnARERkeOTPAAtXboU4eHh0Gg0iImJwfbt25tsv3XrVsTExECj0SAiIgIpKSkWz69fvx6xsbHw9vaGu7s7oqKi8Mknn3RkF5xCVHDdRGjOAyIioq5A0gC0du1aTJ8+HS+//DIyMzMRHx+PMWPGID8/32r73NxcjB07FvHx8cjMzMScOXMwbdo0rFu3ztzGx8cHL7/8MtLT03HgwAE8+uijePTRR/Hjjz92Vre6pKgrFkQkIiJydDIh4fK+w4YNw5AhQ7Bs2TLztn79+mH8+PFITk5u0H7mzJnYsGEDcnJyzNsSExORnZ2N9PT0Rs8zZMgQ3H777XjzzTdbVJder4dWq4VOp4OXl5cNPeq6dBU1GDxvEwBg/6u3wcddJXFFRERElmz5/JZsBMhgMCAjIwMJCQkW2xMSErBz506r+6SnpzdoP3r0aOzbtw81NTUN2gsh8PPPP+Po0aO48cYbG62luroaer3e4kGWtG5KRPRwBwBkcxSIiIgcnGQBqKSkBEajEf7+/hbb/f39UVRUZHWfoqIiq+1ra2tRUlJi3qbT6eDh4QGVSoXbb78d7733Hm677bZGa0lOToZWqzU/goOD29Czrss8EZoBiIiIHJzkk6Cvvsu8EKLJO89ba3/1dk9PT2RlZWHv3r2YP38+kpKSsGXLlkaPOXv2bOh0OvOjoKCgFT3p+qK5ICIREXURLb4Zanvz9fWFQqFoMNpTXFzcYJSnXkBAgNX2Li4u6N69u3mbXC7HtddeCwCIiopCTk4OkpOTcdNNN1k9rlqthlqtbkNvnIP5m2D5F2AyCcjljQdVIiIieybZCJBKpUJMTAzS0tIstqelpWHEiBFW94mLi2vQftOmTYiNjYVSqWz0XEIIVFdXt71oJ9c30BNqFzn0VbXILb0kdTlEREStJuklsKSkJHz44YdYtWoVcnJyMGPGDOTn5yMxMRFA3aWpKVOmmNsnJibi5MmTSEpKQk5ODlatWoWVK1fixRdfNLdJTk5GWloaTpw4gd9++w0LFy7EmjVr8NBDD3V6/7oapUKOgT21AICs/DJpiyEiImoDyS6BAcDEiRNRWlqKefPmobCwEJGRkUhNTUVoaCgAoLCw0GJNoPDwcKSmpmLGjBl4//33ERQUhMWLF2PChAnmNpcuXcIzzzyDU6dOwdXVFX379sWnn36KiRMndnr/uqKoYG/sO3kBWQVlmBDTS+pyiIiIWkXSdYDsFdcBatx3B87guc8zEdnTC9/9PV7qcoiIiMwcYh0gckzRIXUToX8rLEdVjVHiaoiIiFqHAYhsEqTVoIenGrUmgUOndVKXQ0RE1CoMQGSTK+8Mz/WAiIjIUTEAkc2iL98YNZPfBCMiIgfFAEQ2qx8B2p1bitKLXF+JiIgcDwMQ2WxISDcEajUouWjAxA92oVhfJXVJRERENmEAIptplAp8/sRwBGo1+L34Iu5bno7TZZVSl0VERNRiDEDUKuG+7vjqqTgE+7gir7QC96Wk4yRvj0FERA6CAYhaLdjHDV89FYcIX3ecLqvEfcvT8XvxRanLIiIiahYDELVJoNYVXz41HL39PXBWX41JH6TjtyK91GURERE1iQGI2szPU4Mvn4zDgCAvlFw0YNIHu7hIIhER2TUGIGoXPu4qfD51OKKCvVFWUYP7V+xCxskLUpdFRERkFQMQtRutmxKfTh2G68N8UF5Viykrd2PXiVKpyyIiImqAAYjalYfaBasfG4obrvXFJYMRj3y0B9uOnZO6LCIiIgsMQNTu3FQu+PDhWNzS1w9VNSZM/XgffjpyVuqyiIiIzBiAqENolAqkPBSDMZEBMBhNSPw0AxsPFEpdFhEREQAGIOpAKhc53rs/GuOiglBrEvj7F/vxdeYpqcsiIiJiAKKO5aKQY+F9UZgYGwyTAJK+ysaXe/KlLouIiJwcAxB1OIVchuS7B2JKXCiEAGatP4jVv+ZKXRYRETkxBiDqFHK5DG/cOQBP3hgBAJj77RGkbP1D4qqIiMhZMQBRp5HJZJg9pi+m3XItAODt73/Dop+OQQghcWVERORsGICoU8lkMiQl9ME/RvcBACz66Tje+eEoQxAREXUqBiCSxLM3X4tX/9YfAJCy9Q+88e0RhiAiIuo0DEAkmcdvCMf/Gx8JAFi9Mw9zvj4Ek4khiIiIOh4DEEnqoeGhWHDvYMhlwBd78vHif7NRazRJXRYREXVxDEAkuXtieuE/k6KhkMuwPvM0nv8yCzUMQURE1IEYgMgu3DE4CEsfHAKlQoaNBwvx9Kf7UV1rlLosIiLqohiAyG6MHhCAD6bEQu0ix085ZzH1432oNDAEERFR+2MAIrtycx8/fPTIULgqFdh+vASPrt6DS9W1UpdFRERdDAMQ2Z0R1/rik8evh4faBbtOnMfklbuhr6qRuiwiIupCGIDILsWG+eCzqcOgdVVif34ZHlyxGxcuGaQui4iIuggGILJbg4O98cUTw+HjrsLB0zrcv2IXSi5WS10WERF1AQxAZNf6B3lh7ZPD0cNTjd+KyjFxeTqKdFVSl0VERA6OAYjs3nX+nvjqqTgEaTX449wl3Lc8HacuVEhdFhEROTDJA9DSpUsRHh4OjUaDmJgYbN++vcn2W7duRUxMDDQaDSIiIpCSkmLx/IoVKxAfH49u3bqhW7du+Mtf/oI9e/Z0ZBeoE4T7umPtU3EI9nFF/vkKTFy+C3kll6Qui4iIHJSkAWjt2rWYPn06Xn75ZWRmZiI+Ph5jxoxBfn6+1fa5ubkYO3Ys4uPjkZmZiTlz5mDatGlYt26duc2WLVtw//33Y/PmzUhPT0dISAgSEhJw+vTpzuoWdZBgHzf896kRiOjhjtNllbhveTp+Ly6XuiwiInJAMiHhLbiHDRuGIUOGYNmyZeZt/fr1w/jx45GcnNyg/cyZM7Fhwwbk5OSYtyUmJiI7Oxvp6elWz2E0GtGtWzcsWbIEU6ZMaVFder0eWq0WOp0OXl5eNvaKOtq58mo89OFuHD1bju7uKnzy+DD0D+L7RETk7Gz5/JZsBMhgMCAjIwMJCQkW2xMSErBz506r+6SnpzdoP3r0aOzbtw81NdbXiamoqEBNTQ18fHwaraW6uhp6vd7iQfarh6caXzw5HJE9vVB6yYD7V+zCgVNlUpdFREQORLIAVFJSAqPRCH9/f4vt/v7+KCoqsrpPUVGR1fa1tbUoKSmxus+sWbPQs2dP/OUvf2m0luTkZGi1WvMjODjYxt5QZ/NxV+GzqcMRHeINXWUNHlyxGxknz0tdFhEROQjJJ0HLZDKL34UQDbY1197adgB499138cUXX2D9+vXQaDSNHnP27NnQ6XTmR0FBgS1dIIloXZX45PFhuD7cB+XVtZi8cg92/mE9CBMREV1JsgDk6+sLhULRYLSnuLi4wShPvYCAAKvtXVxc0L17d4vtCxYswFtvvYVNmzZh0KBBTdaiVqvh5eVl8SDH4KF2wcePXo/463xRYTDi0Y/2YsvRYqnLIiIiOydZAFKpVIiJiUFaWprF9rS0NIwYMcLqPnFxcQ3ab9q0CbGxsVAqleZt//znP/Hmm2/ihx9+QGxsbPsXT3bFVaXAiimxuLWvH6prTXhyTQY2HbZ+GZWIiAiQ+BJYUlISPvzwQ6xatQo5OTmYMWMG8vPzkZiYCKDu0tSV39xKTEzEyZMnkZSUhJycHKxatQorV67Eiy++aG7z7rvv4pVXXsGqVasQFhaGoqIiFBUV4eLFi53eP+o8GqUCyx6KwZjIABiMJjzz2X58d+CM1GUREZGdkjQATZw4EYsWLcK8efMQFRWFbdu2ITU1FaGhoQCAwsJCizWBwsPDkZqaii1btiAqKgpvvvkmFi9ejAkTJpjbLF26FAaDAffccw8CAwPNjwULFnR6/6hzqVzkeO/+aIyPCkKtSWDaF5lYl3FK6rKIiMgOSboOkL3iOkCOzWgSePnrg/hyb91k9rfuGogHhoVIXBUREXU0h1gHiKijKOQyvHXXQDwcVzeSOOfrg1i1I1fiqoiIyJ4wAFGXJJfLMPfOAXjqxggAwLzvjmDplt8lroqIiOwFAxB1WTKZDLPG9MW0W68DALz7w1EsTDsGXvUlIiIGIOrSZDIZkm7rjZf+2gcAsPjn43j7+98YgoiInBwDEDmFZ266Fq/9rT8AYPm2E5i74TBMJoYgIiJnxQBETuOxG8Ix/65IAMDH6Scx5+uDMDIEERE5JQYgcioPDgvFgnsHQy4DvtxbgBe+ykKt0SR1WURE1MkYgMjp3BPTC/+ZFA2FXIZvss5g6pp9yC4ok7osIiLqRC5SF0AkhTsGB0HtIsdzn2diy9Fz2HL0HAb30mJKXBhuHxQIjVIhdYlERNSBuBK0FVwJ2nkcPqPDh9tzsfFAIQyXL4X5uKswcWgwHhwWgl7d3CSukIiIWsqWz28GICsYgJxPycVqrN1bgM92ncQZXRUAQC4DbunrjylxobjhWl/I5TKJqyQioqYwALURA5DzqjWa8FNOMT7ZlYdffy81b4/wdcdDw0NxT2wveGmUElZIRESNYQBqIwYgAoDfi8vxSfpJrNt/GherawEAbioFxkf3xJS4UPQN4L8NIiJ7wgDURgxAdKWL1bX4ev8prEk/iePFF83brw/3wZS4UIweEAClgl+oJCKSGgNQGzEAkTVCCOw6cR5r0vOw6chZ8yKKfp5qPDAsBA9cHwI/L43EVRIROS8GoDZiAKLmFOoq8cXufHy+pwAlF6sBAC5yGf4aGYApcWEYGtYNMhknTRMRdSYGoDZiAKKWMtSa8P2hQqxJP4mMkxfM2/sGeGJKXBjGRwfBTcXltoiIOgMDUBsxAFFrHD6jwyfpJ/FN1mlU1dStKeSpccG9McGYHBeKcF93iSskIuraGIDaiAGI2kJXUYP/ZhTgk10ncbK0wrw9/jpfTIkLwy19/aDgmkJERO2OAaiNGICoPZhMAluPn8Mn6Sex+Wgx6v/Senq74qHhoZg4NBg+7ippiyQi6kIYgNqIAYjaW35pBT7dfRJr9xZAV1kDAFC5yHHHoCBMiQvF4GBvaQskIuoCGIDaiAGIOkqlwYhvs89gza48HDqtN28fHOyNKcNDeSNWIqI2YABqIwYg6mhCCGQWlGHNzjykHizijViJiNoBA1AbMQBRZ6q/Eeunu06i8KobsT48IhQjr+GNWImIWoIBqI0YgEgKTd2IdXJcKCbE8EasRERNYQBqIwYgkhpvxEpEZDsGoDZiACJ7wRuxEhG1HANQGzEAkb1p7Eas/l5q3H89b8RKRAQwALUZAxDZs0JdJT7fnY8v9uSj5KIBQN2NWG/p64e4a7pjaJgP+gZ4woUjQ0TkZBiA2ogBiBxBda0RPxwqanAjVgBwVykwJLQbYkN9MDSsG6JCvHlTViLq8hiA2ogBiBzN4TM6bDl6Dntyz2P/yQsovzxxup5CLkNkkBeGhvkgNswHsWHd4OuhlqhaIqKOwQDURgxA5MiMJoGjReXYd/I89uZdwN7c8yjSVzVoF+HrjtiwbogN88HQMB+EdXeDTMb1hojIcTEAtREDEHUlQgicLqvEvrwL2JN3HvvyzuPY2YsN2vl6qOoumYXXXTbrH+jFeURE5FAYgNqIAYi6urIKAzJOXsDevAvYl3ceB07pzLfjqOemUiA6xPvyPCIfRId4w13NeUREZL9s+fyW/D/vli5divDwcGg0GsTExGD79u1Ntt+6dStiYmKg0WgQERGBlJQUi+cPHz6MCRMmICwsDDKZDIsWLerA6okck7ebCrf288esMX3xv6dH4MDcBPw3MQ4v/bUPbunrBy+NCyoMRvz6eyn+8/NxPLRyNwa9sQl3vLcDb3x7GN8fLERxecPLakREjkLS/5xbu3Ytpk+fjqVLl2LkyJFYvnw5xowZgyNHjiAkJKRB+9zcXIwdOxZPPPEEPv30U/z666945pln0KNHD0yYMAEAUFFRgYiICNx7772YMWNGZ3eJyCFplAoMvTwXCABMJoFjxeXmEaJ9eRdwuqwSB0/rcPC0Dh/9mgcACOvudnkOUd1coghfd84jIiKHIOklsGHDhmHIkCFYtmyZeVu/fv0wfvx4JCcnN2g/c+ZMbNiwATk5OeZtiYmJyM7ORnp6eoP2YWFhmD59OqZPn25TXbwERtRQ3TyiujC0N+88jp4tx9X/79HdXYWY0G6Xv23WDZE9tVypmog6jS2f35KNABkMBmRkZGDWrFkW2xMSErBz506r+6SnpyMhIcFi2+jRo7Fy5UrU1NRAqWzdjSKrq6tRXV1t/l2v17fqOERdWU9vV/SM6olxUT0BALrKGuw/WReG9uVdQNapMpReMmDTkbPYdOQsAECjlCM6uJt5hCg6xBuevKErEdkByQJQSUkJjEYj/P39Lbb7+/ujqKjI6j5FRUVW29fW1qKkpASBgYGtqiU5ORlvvPFGq/YlclZaVyVu7uuHm/v6AahbmPHQaZ35stnevAvQVdYg/UQp0k/U3d1eLgP6BXqZR4iGhvnAn7fwICIJSP6VjqvnCwghmpxDYK29te22mD17NpKSksy/6/V6BAcHt/p4RM5I7aJATKgPYkJ9gFHXwGQS+OPcxctfva8bKTp1oRKHz+hx+Iweq3fmAQBCfNzMYWhoWDdc08OD84iIqMNJFoB8fX2hUCgajPYUFxc3GOWpFxAQYLW9i4sLunfv3upa1Go11GquikvUnuRyGa7z98R1/p54cFgogLr7mO27YoQop0iP/PMVyD9fgfX7TwMAPDUuuKaHByJ83RHu646IHh4Iv/yzq0ohZZeIqAuRLACpVCrExMQgLS0Nd911l3l7Wloaxo0bZ3WfuLg4fPvttxbbNm3ahNjY2FbP/yGizhOodcUdg11xx+AgAIC+qgaZ+WXYm3see/POI6ugDOVVtcgqKENWQVmD/YO0GnMgiuhxOSD5eqBnN1co5Bw1IqKWk/QSWFJSEiZPnozY2FjExcXhgw8+QH5+PhITEwHUXZo6ffo01qxZA6DuG19LlixBUlISnnjiCaSnp2PlypX44osvzMc0GAw4cuSI+efTp08jKysLHh4euPbaazu/k0TUKC+NEqN698Co3j0AAIZaE/44dxG5JZdw4txFnCi5dPnnS9BV1uCMrgpndFXY8XuJxXFUCjlCu7tdDkUeiOjhbh5B8nFX8ZIaETUg+UrQS5cuxbvvvovCwkJERkbi3//+N2688UYAwCOPPIK8vDxs2bLF3H7r1q2YMWMGDh8+jKCgIMycOdMcmAAgLy8P4eHhDc4zatQoi+M0hV+DJ7IvQghcqKgxh6IT5y4ht6QuKOWVVDRYxfpKWlelecQo4qpLaholL6kRdSW8FUYbMQAROQ6jSeBMWaV55CjXHJAu4XRZZZP79vR2tbycdnnuUZA3L6kROSIGoDZiACLqGioNRuSVXrK4pHbiXN3P+qraRvdTucgR1t3NYhL2NZcvr/m4qzqxB0RkC4dYCJGIqKO5qhToF+iFfoGW/0cohMD5S4a6YHTFJbUT5y7hZGkFDLUmHDt7EcfOXgRw1mJfbzel+RLaNVdMyA7rzktqRI6EI0BWcASIyHldfUmt/nLaiXMXcUbX+A1gZTIgSOtqvpwW1t0dgVoN/LUa+Htp4Oep5m1BiDoYL4G1EQMQEVlTf0ntyhGjEyXNX1ID6gJSd3c1/L3UCPC6HIw8NQjQquHvpUHA5d+93ZT81hpRK/ESGBFRB2jRJbXLoSj//CUU6apwVl+N4vIq1BgFSi5Wo+RiNQ6fafx+g2oXOfy9NPD3uhyMvOpGkPy19T/XbeflNqK2YQAiImojmUyG7h5qdPdQIzbMp8HzJpPA+QoDinRVKC6vQpGuGkX6KhTrq1Ckr7oclKpwoaIG1bUm8+rYTfF2U8Lfsz4YqS+HprqQFKDVwM9LDV93NeT8NhuRVQxAREQdTC6XwddDDV8PNQBto+2qaow4V14Xjs5eEYzO6i23VdeaUFZRg7KKGhw9W97o8VzkMvh5quF3VTAKuPy73+VtHmp+FJDz4b96IiI7oVEqEOzjhmAft0bbCCGgr6ytGzm6HIrO6q74+XJYKrlYjVqTMK+e3RQPtUvDYOSlvhyYNPBxU6GbuwpeGhfOT6IugwGIiMiByGQyaN2U0Lop0SfAs9F2NUYTSi5WNxxFuiosXayurXucq8WJc5eaPLeLXAZvNyW6XQ5E3dyU8HFX1f1+xbZu7qq60OSmgqfGhZfhyC4xABERdUFKhRyBWlcEal2bbHexuvaqUaRq86W2In0VzpVX4/wlAyprjKg1CZRcNKDkoqHFdSjkMnRzU8Lb7XIocrcMUN3cVPBxV9U9fzk4MTRRZ2AAIiJyYh5qF3j08MA1PTyabFdVY8SFCgMuXKrBhQoDzl8yoKzCgPOXf/9zWw3OX6r7vcJghLEVoUkuA7q5qeBtZYTJx71hmPJxV8FLo2RoIpswABERUbM0SkWLRpSuVFVjRFnF5YB0yYDzFQZcqKjBhcsBqW5bzeUgVReeLlbXwiSA0ksGlF4y4I9mLsvVk8sAb7c/R5WsXY7zcnWBp0YJT82V/+sCtQuXFHBGDEBERNQhNEoFArQKBGg1Ld6nuvbP0GQxonTpcni6cvSpwoCySzUovxyazl+qew5oWWiqp3KRw+uqUOSpbhiUvKyEp/qfuS6T42EAIiIiu6F2UcDfSwF/r5aHJkOtCWWXR5euDEfm0HR5xElfVYvyqhqUV9WivKpu8nf9/rZepruaSiFvEIrqf/ZiiLJLDEBEROTQVC5y+F3++r4tjCaBi9WWoejPn2suBybLbQ3a1Ycoo8l82a7V/WgiRP0Zplzg5aqE9vLDS1P3jUCtqxLuKgWXKbABAxARETklhVxmDhKtZTIJXDRYD0r6qwNVpfUQddFQCyHaHqJc5DJzOLoyJGldXa74+c/nvTSXf3dTwlPtfGs8MQARERG1klwug9fly1xAyyeIX8mWEFVeVQtdZY35ob/8vzVGgVqTuGIelI39kOGq0HR1iLrqOc2fPzvqsgUMQERERBJqa4gSQqCqxmQRjK5+6Jt4zlBrgknAfHsVW8lkgKfaxXwprqkQdWVw6uamgtat9aNvbcUARERE5MBkMhlcVQq4qmz7xl29qhpjkwGpqRBVVWOCEID+8mhVASpbfN4BQV7YOC3e5nrbCwMQERGRE9MoFdAoFTZPIgfqli3QV9a2aKTp6ufbMveqPTAAERERUauoXRTo4alAD0+1zfuaTKIDKmo5uaRnJyIiIqck9cRpBiAiIiJyOgxARERE5HQYgIiIiMjpMAARERGR02EAIiIiIqfDAEREREROhwGIiIiInA4DEBERETkdBiAiIiJyOgxARERE5HQYgIiIiMjpMAARERGR02EAIiIiIqfjInUB9kgIAQDQ6/USV0JEREQtVf+5Xf853hQGICvKy8sBAMHBwRJXQkRERLYqLy+HVqttso1MtCQmORmTyYQzZ87A09MTMplM6nI6jF6vR3BwMAoKCuDl5SV1OR3OmfrLvnZdztRf9rXr6qj+CiFQXl6OoKAgyOVNz/LhCJAVcrkcvXr1krqMTuPl5eUUf3D1nKm/7GvX5Uz9ZV+7ro7ob3MjP/U4CZqIiIicDgMQEREROR0GICemVqvx+uuvQ61WS11Kp3Cm/rKvXZcz9Zd97brsob+cBE1EREROhyNARERE5HQYgIiIiMjpMAARERGR02EAIiIiIqfDANRFJScnY+jQofD09ISfnx/Gjx+Po0ePNrnPli1bIJPJGjx+++23Tqq69ebOndug7oCAgCb32bp1K2JiYqDRaBAREYGUlJROqrZtwsLCrL5Pzz77rNX2jvS+btu2DXfccQeCgoIgk8nwzTffWDwvhMDcuXMRFBQEV1dX3HTTTTh8+HCzx123bh369+8PtVqN/v374+uvv+6gHtimqf7W1NRg5syZGDhwINzd3REUFIQpU6bgzJkzTR5z9erVVt/vqqqqDu5N05p7bx955JEGNQ8fPrzZ49rje9tcX629PzKZDP/85z8bPaa9vq8t+ayx179bBqAuauvWrXj22Wexa9cupKWloba2FgkJCbh06VKz+x49ehSFhYXmx3XXXdcJFbfdgAEDLOo+ePBgo21zc3MxduxYxMfHIzMzE3PmzMG0adOwbt26Tqy4dfbu3WvRz7S0NADAvffe2+R+jvC+Xrp0CYMHD8aSJUusPv/uu+9i4cKFWLJkCfbu3YuAgADcdttt5vv3WZOeno6JEydi8uTJyM7OxuTJk3Hfffdh9+7dHdWNFmuqvxUVFdi/fz9effVV7N+/H+vXr8exY8dw5513NntcLy8vi/e6sLAQGo2mI7rQYs29twDw17/+1aLm1NTUJo9pr+9tc329+r1ZtWoVZDIZJkyY0ORx7fF9bclnjd3+3QpyCsXFxQKA2Lp1a6NtNm/eLACICxcudF5h7eT1118XgwcPbnH7l156SfTt29di21NPPSWGDx/ezpV1vOeff15cc801wmQyWX3eUd9XAOLrr782/24ymURAQIB4++23zduqqqqEVqsVKSkpjR7nvvvuE3/9618tto0ePVpMmjSp3Wtui6v7a82ePXsEAHHy5MlG23z00UdCq9W2b3HtzFpfH374YTFu3DibjuMI721L3tdx48aJW265pck2jvC+CtHws8ae/245AuQkdDodAMDHx6fZttHR0QgMDMStt96KzZs3d3Rp7eb48eMICgpCeHg4Jk2ahBMnTjTaNj09HQkJCRbbRo8ejX379qGmpqajS203BoMBn376KR577LFmb9zrqO9rvdzcXBQVFVm8b2q1GqNGjcLOnTsb3a+x97qpfeyVTqeDTCaDt7d3k+0uXryI0NBQ9OrVC3/729+QmZnZOQW20ZYtW+Dn54fevXvjiSeeQHFxcZPtu8J7e/bsWWzcuBGPP/54s20d4X29+rPGnv9uGYCcgBACSUlJuOGGGxAZGdlou8DAQHzwwQdYt24d1q9fjz59+uDWW2/Ftm3bOrHa1hk2bBjWrFmDH3/8EStWrEBRURFGjBiB0tJSq+2Liorg7+9vsc3f3x+1tbUoKSnpjJLbxTfffIOysjI88sgjjbZx5Pf1SkVFRQBg9X2rf66x/Wzdxx5VVVVh1qxZeOCBB5q8eWTfvn2xevVqbNiwAV988QU0Gg1GjhyJ48ePd2K1thszZgw+++wz/PLLL/jXv/6FvXv34pZbbkF1dXWj+3SF9/bjjz+Gp6cn7r777ibbOcL7au2zxp7/bnk3eCfw3HPP4cCBA9ixY0eT7fr06YM+ffqYf4+Li0NBQQEWLFiAG2+8saPLbJMxY8aYfx44cCDi4uJwzTXX4OOPP0ZSUpLVfa4eMRGXF0VvbiTFnqxcuRJjxoxBUFBQo20c+X21xtr71tx71pp97ElNTQ0mTZoEk8mEpUuXNtl2+PDhFpOHR44ciSFDhuC9997D4sWLO7rUVps4caL558jISMTGxiI0NBQbN25sMhw4+nu7atUqPPjgg83O5XGE97Wpzxp7/LvlCFAX9/e//x0bNmzA5s2b0atXL5v3Hz58uF39F0ZLubu7Y+DAgY3WHhAQ0OC/JIqLi+Hi4oLu3bt3RoltdvLkSfz000+YOnWqzfs64vta/60+a+/b1f+lePV+tu5jT2pqanDfffchNzcXaWlpTY7+WCOXyzF06FCHe78DAwMRGhraZN2O/t5u374dR48ebdXfsL29r4191tjz3y0DUBclhMBzzz2H9evX45dffkF4eHirjpOZmYnAwMB2rq7jVVdXIycnp9Ha4+LizN+eqrdp0ybExsZCqVR2Rolt9tFHH8HPzw+33367zfs64vsaHh6OgIAAi/fNYDBg69atGDFiRKP7NfZeN7WPvagPP8ePH8dPP/3UqnAuhEBWVpbDvd+lpaUoKChosm5Hfm+BuhHcmJgYDB482OZ97eV9be6zxq7/btttOjXZlaefflpotVqxZcsWUVhYaH5UVFSY28yaNUtMnjzZ/Pu///1v8fXXX4tjx46JQ4cOiVmzZgkAYt26dVJ0wSYvvPCC2LJlizhx4oTYtWuX+Nvf/iY8PT1FXl6eEKJhX0+cOCHc3NzEjBkzxJEjR8TKlSuFUqkU//vf/6Tqgk2MRqMICQkRM2fObPCcI7+v5eXlIjMzU2RmZgoAYuHChSIzM9P8rae3335baLVasX79enHw4EFx//33i8DAQKHX683HmDx5spg1a5b5919//VUoFArx9ttvi5ycHPH2228LFxcXsWvXrk7v39Wa6m9NTY248847Ra9evURWVpbF33F1dbX5GFf3d+7cueKHH34Qf/zxh8jMzBSPPvqocHFxEbt375aii2ZN9bW8vFy88MILYufOnSI3N1ds3rxZxMXFiZ49ezrke9vcv2MhhNDpdMLNzU0sW7bM6jEc5X1tyWeNvf7dMgB1UQCsPj766CNzm4cffliMGjXK/Ps777wjrrnmGqHRaES3bt3EDTfcIDZu3Nj5xbfCxIkTRWBgoFAqlSIoKEjcfffd4vDhw+bnr+6rEEJs2bJFREdHC5VKJcLCwhr9PyJ79OOPPwoA4ujRow2ec+T3tf4r+1c/Hn74YSFE3VdqX3/9dREQECDUarW48cYbxcGDBy2OMWrUKHP7ev/9739Fnz59hFKpFH379rWb8NdUf3Nzcxv9O968ebP5GFf3d/r06SIkJESoVCrRo0cPkZCQIHbu3Nn5nbtKU32tqKgQCQkJokePHkKpVIqQkBDx8MMPi/z8fItjOMp729y/YyGEWL58uXB1dRVlZWVWj+Eo72tLPmvs9e9WdrkDRERERE6Dc4CIiIjI6TAAERERkdNhACIiIiKnwwBERERETocBiIiIiJwOAxARERE5HQYgIiIicjoMQEREROR0GICIiIjI6TAAEVG7kMlkTT4eeeQRm4+Zl5cHmUyGrKwsi+2PPPIIxo8f3y5122r16tXw9vaW5NxE1H5cpC6AiLqGwsJC889r167Fa6+9hqNHj5q3ubq6WrSvqamBUqnstPqsMRgMUKlUktZARNLgCBARtYuAgADzQ6vVQiaTmX+vqqqCt7c3vvrqK9x0003QaDT49NNPYTKZMG/ePPTq1QtqtRpRUVH44YcfzMcMDw8HAERHR0Mmk+Gmm27C3Llz8fHHH+P//u//zKNLW7ZsAQCcPn0aEydORLdu3dC9e3eMGzcOeXl55uPVjxwlJycjKCgIvXv3ttqX7Oxs3HzzzfD09ISXlxdiYmKwb98+bNmyBY8++ih0Op353HPnzgVQF6Zeeukl9OzZE+7u7hg2bJi5LuDPkaNvvvkGvXv3hkajwW233YaCgoJmz0tE7Y8jQETUaWbOnIl//etf+Oijj6BWq/Gf//wH//rXv7B8+XJER0dj1apVuPPOO3H48GFcd9112LNnD66//nr89NNPGDBgAFQqFVQqFXJycqDX6/HRRx8BAHx8fFBRUYGbb74Z8fHx2LZtG1xcXPD//t//w1//+lccOHDAPNLz888/w8vLC2lpaWjsXtAPPvggoqOjsWzZMigUCmRlZUGpVGLEiBFYtGiRxeiWh4cHAODRRx9FXl4evvzySwQFBeHrr7/GX//6Vxw8eBDXXXcdAKCiogLz58/Hxx9/DJVKhWeeeQaTJk3Cr7/+2uR5iagDtOu95YmIhBAfffSR0Gq15t9zc3MFALFo0SKLdkFBQWL+/PkW24YOHSqeeeYZi/0yMzMt2jz88MNi3LhxFttWrlwp+vTpI0wmk3lbdXW1cHV1FT/++KN5P39/f1FdXd1k/Z6enmL16tUt6psQQvz+++9CJpOJ06dPW2y/9dZbxezZs837ARC7du0yP5+TkyMAiN27dzd7XiJqX7wERkSdJjY21vyzXq/HmTNnMHLkSIs2I0eORE5Ojs3HzsjIwO+//w5PT094eHjAw8MDPj4+qKqqwh9//GFuN3DgwGbn/SQlJWHq1Kn4y1/+grfffttif2v2798PIQR69+5tPreHhwe2bt1qsa+Li4vFa9C3b194e3ub+2vreYmo9XgJjIg6jbu7e4NtMpnM4nchRINtLWEymRATE4PPPvuswXM9evRosoarzZ07Fw888AA2btyI77//Hq+//jq+/PJL3HXXXY2eW6FQICMjAwqFwuK5+ktk9az1rX6breclotbjCBARScLLywtBQUHYsWOHxfadO3eiX79+AGAeqTEajRZtVCpVg21DhgzB8ePH4efnh2uvvdbiodVqba6vd+/emDFjBjZt2oS7777bPN/I2rmjo6NhNBpRXFzc4NwBAQHmdrW1tRaTmo8ePYqysjL07du32fMSUftiACIiyfzjH//AO++8g7Vr1+Lo0aOYNWsWsrKy8PzzzwMA/Pz84Orqih9++AFnz56FTqcDAISFheHAgQM4evQoSkpKUFNTgwcffBC+vr4YN24ctm/fjtzcXGzduhXPP/88Tp061eKaKisr8dxzz2HLli04efIkfv31V+zdu9ccysLCwnDx4kX8/PPPKCkpQUVFBXr37o0HH3wQU6ZMwfr165Gbm4u9e/finXfeQWpqqvnYSqUSf//737F7927s378fjz76KIYPH47rr7++2fMSUTuTehISEXU9jU2Cvnoys9FoFG+88Ybo2bOnUCqVYvDgweL777+3aLNixQoRHBws5HK5GDVqlBBCiOLiYnHbbbcJDw8PAUBs3rxZCCFEYWGhmDJlivD19RVqtVpERESIJ554Quh0OiGE9cnTV6uurhaTJk0SwcHBQqVSiaCgIPHcc8+JyspKc5vExETRvXt3AUC8/vrrQgghDAaDeO2110RYWJhQKpUiICBA3HXXXeLAgQMWr8m6detERESEUKlU4pZbbhF5eXktPi8RtR+ZEI18D5SIiNrN6tWrMX36dJSVlUldChGBl8CIiIjICTEAERERkdPhJTAiIiJyOhwBIiIiIqfDAEREREROhwGIiIiInA4DEBERETkdBiAiIiJyOgxARERE5HQYgIiIiMjpMAARERGR0/n/9beThs2C+NkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = 0.5\n",
    "DME_s = DME(s, rho, sigma)\n",
    "DME_goal = evolve(s,rho) @ sigma @ np.conjugate(evolve(s, rho)).T\n",
    "error1_s = np.linalg.norm(np.abs(DME_goal - DME_s))\n",
    "trotter_steps_list = np.linspace(1,20,10)\n",
    "error1_list = error1_s * np.ones(len(trotter_steps_list))\n",
    "trotter_error_list = []\n",
    "trotter_error_list2 = []\n",
    "for M in trotter_steps_list:\n",
    "    trotter_step = int(M)\n",
    "    DME_trotter = DME_trotterize(s, rho, sigma, trotter_step)\n",
    "    # DME_trotter2 = np.linalg.matrix_power(evolve(s/trotter_step,rho) @ sigma @ np.conjugate(evolve(s/trotter_step, rho)).T, trotter_step)\n",
    "    # Use trace norm\n",
    "    error = np.linalg.norm(DME_goal - DME_trotter)\n",
    "    # error2 = op_norm(DME_goal - DME_trotter2)\n",
    "    trotter_error_list.append(error)\n",
    "    # trotter_error_list2.append(error2)\n",
    "\n",
    "plt.plot(trotter_steps_list, trotter_error_list, label=r'$|E_s^\\rho\\sigma (E_s^\\rho)^\\dag - (\\varepsilon_{s/M}(\\sigma))^M|$')\n",
    "# plt.plot(trotter_steps_list, trotter_error_list2, label=r'$ \\varepsilon_s(\\sigma)-(E_{s/M}\\sigma E_{s/M}^\\dag)^M$')\n",
    "plt.plot(trotter_steps_list, error1_list, label=r'$|E_s^\\rho\\sigma (E_s^\\rho)^\\dag-\\varepsilon_{s}(\\sigma)|$')\n",
    "plt.legend()\n",
    "plt.xlabel('Trotter steps')\n",
    "plt.ylabel('Trace norm')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DBF",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}