Ex1_Train_Reuploading_Circuit_Pytorch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Pytorch version for QNN trained on Data\n",
    "\n",
    "This code can reproduce the QNN prediction map and the Hessian distribution with pytorch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from vqc_loss_landscapes.torchcirq import *\n",
    "from vqc_loss_landscapes.data_helper import *\n",
    "\n",
    "import torch\n",
    "from torch.autograd import Variable\n",
    "from vqc_loss_landscapes.complex import *\n",
    "import qutip.qip.circuit as QCirc\n",
    "from tqdm import tqdm as tqdm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = \"cpu\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_circle_data(path):\n",
    "    np_file = np.load(path+\".npy\", allow_pickle=True)\n",
    "    return np_file.item()\n",
    "\n",
    "def test_loss(circ, params, x_train, y_train, init, measurements=None, device=\"cpu\"):\n",
    "    loss = 0.0\n",
    "    i = 0\n",
    "    for x,y in zip(x_train, y_train):\n",
    "        # print(\"test loss: \", i)\n",
    "        # i+=1\n",
    "        y = 2*(y - 1/2)\n",
    "        out1 = matmul(circ(params, x=x, device=device), init)\n",
    "        out1_copy = out1.clone().detach().requires_grad_(True)\n",
    "        out2 = matmul(measurements, out1)\n",
    "        out1 = inner_prod(out1_copy, out2) # this is already <Psi|0><0|Psi>\n",
    "        # print(out1[0])\n",
    "        loss += (y-out1[0])**2\n",
    "    return loss/len(x_train)\n",
    "\n",
    "def predicted_labels(output_fidelity):\n",
    "    output_labels = [np.argmax(o) for o in output_fidelity]\n",
    "    return np.array(output_labels)\n",
    "\n",
    "\n",
    "def prediction(params, x=None, y=None, measurements = None, device=\"cpu\"):\n",
    "    \"\"\"For prediction of model measure in both directions\"\"\"\n",
    "    output_fidelity = []\n",
    "    for i in tqdm(range(len(x))):\n",
    "        fidelities = []\n",
    "        out1 = matmul(circ(params, x=x[i], device=device), init)\n",
    "        out1_copy = out1.clone().detach().requires_grad_(True)\n",
    "        out2 = matmul(measurements, out1)\n",
    "        out1 = inner_prod(out1_copy, out2)\n",
    "        fidelities.append(out1[0].detach().cpu())\n",
    "        del out1\n",
    "        output_fidelity.append(fidelities)\n",
    "    #predicted = predicted_labels(output_fidelity)\n",
    "    return output_fidelity\n",
    "\n",
    "def get_fidelity():\n",
    "    fidl= prediction(params, x=X_map, y=None, measurements = model.measures, device=device)\n",
    "    c = np.array(fidl)\n",
    "    fidelity = []\n",
    "    for i in c[:,0]:\n",
    "        fidelity.append(i.item())\n",
    "    fidelity = np.array(fidelity)\n",
    "    return fidelity\n",
    "\n",
    "\n",
    "def plot_prediction_map(epoch):\n",
    "    fidelity = get_fidelity()\n",
    "    Z = fidelity.reshape(grid_resolution, grid_resolution)\n",
    "    fig = plt.figure(figsize=(10,10))\n",
    "    ax = fig.add_subplot(111)\n",
    "    cont = ax.contourf(X_grid,Y_grid,Z, 20, cmap='magma')\n",
    "    plt.ylim(-1, 1)\n",
    "    cont.set_clim(-0.4, 1.5)\n",
    "    plt.show()\n",
    "    \n",
    "#============================================================================\n",
    "# Mapping Parameters\n",
    "# ==========================================================================\n",
    "grid_resolution = 21\n",
    "\n",
    "x_path = torch.linspace(-1,1,grid_resolution).to(device)\n",
    "y_path = torch.linspace(-1,1,grid_resolution).to(device)\n",
    "X_grid, Y_grid = torch.meshgrid(x_path.cpu(), y_path.cpu())\n",
    "\n",
    "X = []\n",
    "for i in range(len(x_path)):\n",
    "    for j in range(len(y_path)):\n",
    "        X.append([x_path[i], y_path[j], 0])\n",
    "X_map = torch.tensor(X).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if torch.cuda.is_available():\n",
    "    device = torch.device('cuda:0')\n",
    "else:\n",
    "    device = \"cpu\"\n",
    "\n",
    "width = 2\n",
    "layers = 4\n",
    "batch_size = 16\n",
    "epochs = 100\n",
    "train_samples = 100\n",
    "test_samples = 100\n",
    "epsilon = 0\n",
    "\n",
    "# ==========================================================================\n",
    "# Data Generating\n",
    "# ==========================================================================\n",
    "X, y_train = generate_circle_data(train_samples)\n",
    "X_train = np.zeros((train_samples,3))\n",
    "X_train[:,0:2] = X[:] # add extra dim x_3 = 0\n",
    "X, y_test = generate_circle_data(test_samples)\n",
    "X_test = np.zeros((test_samples,3))\n",
    "X_test[:,0:2] = X[:]\n",
    "\n",
    "X_train = torch.tensor(X_train).to(device)\n",
    "y_train = torch.tensor(y_train).to(device)\n",
    "X_test = torch.tensor(X_test).to(device)\n",
    "y_test = torch.tensor(y_test).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start training width2 layers4 epochs100\n",
      "circuit built\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1730277d1784781948d608430d6f0f8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5238d35b1e6a46628623358349df65bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=441.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAJDCAYAAABZrSP4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df7BkZ33f+c93RoxmMGX0a0YWAiKRlb12ipRwVHJ2vRXHGIHiPxBZE1u4vBEJ1JS9Jlu1LqcsilTMKvEuZKuWTSqUjYJlZG8WYZNiPSmL0opfmz8MjuRdjEAUaBCKGEtmFImxCzOjYZhn/+g+zJm+/eN0nx/P9/s871fV1Mzt2/fec/v27X7P9znntKWUBAAAAN/25d4AAAAAbEa0AQAABEC0AQAABEC0AQAABEC0AQAABEC0AQAABDBItJnZPWZ20sw+v+L9Zmb/ysyOm9nnzOyHW++7w8wem/+5Y4jtAQAAKM1Qk7YPSLp1zfv/jqQb5n+OSvp1STKzKyT9qqQfkXSzpF81s8sH2iYAAIBiDBJtKaX/IOm5NVe5TdJvp5nPSLrMzK6R9DpJD6aUnkspfUPSg1offwAAAFWaap+2ayV9rfX2ifllqy4HAABAyyUTfR1bcllac/neT2B2VLOlVV1qL/gb1xy8aritAzI5l5b9CsTy7fO5t2A3586X+RJ+3+alCYFJvcA2P44/e+7p/5xSOtz3a00VbSckvaz19kslPTW//G8vXP6pZZ8gpXS3pLsl6foXviS98/uPjrGdwKSeORv7AO6TZ3JvwW6eOV1e2Jw883zuTQCqdOTgpRuv84E/e+d/GuJrTRVtxyS9zczu0+yggz9PKT1tZg9I+p9bBx+8VtLbJ9qm6lx98Gzn6379zIERtwTIp7RgI9aAfLoE25AGiTYz+6BmE7OrzOyEZkeEvkCSUkq/Iel+ST8p6bikb0n6B/P3PWdm/0zSQ/NPdVdKad0BDehom0Db5uOJOTQiTtkINgCRDRJtKaU3bXh/kvSLK953j6R7htiOGvWNs6G+HjEH70oKNmINyG/qKZs03fIoepo6zra1bPsIuXJFm7KVEmzEGlA3os0h74HWFVM5eFBCsBFrgC85pmwS0eZKKbG2ydUHzxJugUWaskUPNmIN8CdXsElEmwu1xFob4YaxRQ42Yg3AMkRbRjXGWlvz/RNvcUSaskVErAG+5ZyySURbFrXH2iKmbhhaxCkbwebDSZ3a6eOO6LKBtwTYi2ibGMG2HOHmX5QpW7RgI9aGtWt0jfV1ibly5J6ySUTbZIi1zQg39BUp2Ii1/nIF2jbWbSNBh20RbSMj1rZTW7gdPnA+xOuPRpiyEWxlixBo22I6F4eHKZtEtI2GWNsdByhgWwRbWUoMtG0s+/4JOUhE2ygItmHUNnXDbqIEG7G2Wu2R1sVJnSLcMvEyZZOItkERa8Mj3PLzvDRKsMVDoO2uue2It+l4CjaJaBsEsTYuwg3LEGwxEGnDY+pWL6KtJ4JtGoRbHl6nbBGCreZYI9TGx9RtfN6mbBLRtjNibXocoACJYPOKUMuDqVtdiLYtEWv5MXWrF8HmC6HmA1O34XmcskmS/xNEOXH1wbMEmyP8LMbnbWmUYPPjpE4RbA7xMykfk7YOCASfmLjVw3uw1RBrBEEMTN368zplk4i2tYg1/wi3cXibsnlWcrARanGxr1uZiDaExwEKZfM8ZSsx2Ai1cjB1257nKZvEPm0rMWWLh5/ZMDxN2bwG28kzzxcXbOynVi5+rt14DzaJSdtSPPnHxXJpOTwHWyl4Mq8HU7cyEG0LCLb4CDeMpYRgI9Tqxr5uy0WYsklE23cRa2WJtJ/b4QPn9cxZH3sqeFka9ThlI9iA8kSJtQbRJoKtZEzdMITowUasocGUbSZarDWqjzaCDfDF25SNYAPKETXWGj7WZDLgFQ7qwc8ZtSLY0Fb7lC16sEmVRhtP4vXhZ45dRJ6yEWzABSUEm1Th8ihP3vDI08EIOXlaGiXYUJJap2ylxFqjmmcJlkPBzx9dEWxAfKUFm1RJtPFkjQb3BZSMYMMyNU7ZSgw2qYJo40kaiMHL0mjUKRvBBsyUGmxSwfu0EWtYhXO3YRWCDaWpacpWcqw1ipy0EWzYxON95PCB87k3AcHwIu/ATA3BJhUYbR6fjAH4F23KRqxhk1qmbLUEm1TQ8iixhm2xTHqxnK876mV/tigINqCuWGsUEW0E23qHD25+Nn7mzMEJtsQfwg1SrCkbwYYuSp+y1RhsUgHRVnKwdYmtMb5WrQHnASfZxToEG1BvsEmBo620WJsy0DapLeCYtuWVe2k0ypSNYENXpU7Zao61Rshoe8G+MvZ/8RRqq9QScIQbPCPY0BXBVraQ0RZdhFhbppaAQz0iTNkINtSOYLuAaJtI1FBbpcSAY9o2vZxLowQbSlPalI1Y24toG1lpsbZMSQHnIdw4GAHEGmpXSrAdPmSDfj6ibQQ1hNoqJQUcyuV5ykawYRclTdkIttWItgHVHGvLRA04D9O2GuQ+atQjgg01I9Y2I9p6ItS6iRZwhFu5vE7ZCDbsqoQpWwnBNmasNdhxZkeHD54h2HbE7bYZLx5fF4INNSPYumPStgViYziHD55xP3Fj2jaeXEujHqdsBBv6iD5lix5sU8Vag2jrgFgbR3O7eo63WsIt54vF14xgQ80iB9vUsdZgeXQNlkCnwW2MKXibshFs6CvylI1g2w3RtgIhMS3Pt3dpr3ObG0eNAv0RbNM7fMiyBpvE8ugenuOhdBH2c5sSJ9kdDlM2wIeIwZY71Np4Rmgh2PLz+jNg2gbAi6hTtmjB5mGytohJm/yGQq2YuE3vyEEORpgKUzb0QbCNz1uotVU/aSPYfPL4c2HaBiAngm18noNNqnjS5jEKcDEmbhiCp/3ZmLKhNlGCzXusNaqctBFscdT+s+KVEQBIMadsEYLN435r61QVbZx3LSZPPzOWSLErpmzYFcE2vGix1qgm2jw98WN7/PwAIIYIwRbVINFmZrea2ZfM7LiZ3bnk/e8xs8/O/3zZzE613ved1vuODbE9bUzXyuHl58i0bXdTn1jX0/5swC6iTdk8B1vU6Vpb7wMRzGy/pPdKukXSCUkPmdmxlNKjzXVSSv9j6/r/SNKrWp/idErpxr7bsYyXJ3kAdWNpFLsg2IYTPdYaQ0zabpZ0PKX0eErprKT7JN225vpvkvTBAb7uSkzXyuXl5zrltI2DEQB4R7BNY4hou1bS11pvn5hftoeZ/RVJ10v6ROvig2b2sJl9xsze0HdjvDypYzz8jBEJUzbsIsqU7cjBS90GWwnLoYuGOE/bsltk1Y4rt0v6cErpO63LXp5SesrMXiHpE2b2SErpK3u+iNlRSUcl6ciB793ziXkir4uHc7hdffCsvn7mQNZtwGrsz4aoIgWbV6XFWmOISdsJSS9rvf1SSU+tuO7tWlgaTSk9Nf/7cUmf0sX7u7Wvd3dK6aaU0k0vfsELL3ofwVYnfu7wjikbSkWw5TFEtD0k6QYzu97MDmgWZnuOAjWzH5B0uaRPty673Mwunf/7Kkk/KunRxY9dpcR91678ntMX/cF6uX/+HEkKYEgRpmxeg63E5dBFvZdHU0rnzOxtkh6QtF/SPSmlL5jZXZIeTik1AfcmSfellNpLpz8o6X1mdl6zgHxX+6jTdXI/WQ9pXZyNHW7P/uWhUT//FDwslY7t8IHzeuZsNadV7M3D0ihTNpTIc7DVYJDXHk0p3S/p/oXL/unC2+9c8nF/KOmV2369Syz+0XRepmjt7YgccDnDjX3bupn6HG1ANN6nbARbftW+YPzUvETaOs02Ro23GiZuAMpEsG2vplhrEG0jihBqy0SON8INHrA0ipIQbLs7MvDTEdE2sKihtkzUeMsRbiyR+uFhfzZgG56nbATbboaOtQbRNoCSQm2ZiPFW4sSNgxFiYMqGbRBs3UWINWm8YJOItp2VHmrLRIu3qcONaRuAUhBs2xsz1hpE2xZqDLVlIsVbiRM3+MWUDdvwOmUj2LYzRaw1iLYNCLXVosTblOHGtC0v9mdDFARbN56DbcpYaxBtaxBs3USJN2BMTNkQnadg8xxrUp5gk4i2pYi13XiOt1KWSaMejMCJdYELvE7ZvPAcbLlirRHv0X9kBFt/td+GvB4pgGi8TNkItvWItpbaY2NI3JaYWs792VgaBfoj2DYj2jQLDCJjeN5u08MHz+TeBACV87g06mHKRrB1U/0+bd7CojSe93MDAGAVT7HWqHrSRrBNh9sapWJpFJExZVvOY7BJFUcbETE9D7c5S6TAuE7aE9/9g4t5XBrNjWDbTnXRxv5redVy23MEKWpEqMXiYcrmjedgkyrbp62WYPDuyu85zT5uGBSvhJAXsYZdeJuyeQ82qaJJG8GGBkukGAr7s60PNmLuAm9Lo7mnbATbbqqYtBFs/jBtA2IjyOLKHWyeRIm1RvGTNoINpTl84HzuTUDltgk24g6LvEzZogWbVPCkjVjzL+e0rZTXIkU+NS6NEmC78bQ0mnvKRrD1U+SkjWCLo+SfFUeQzvBi8WUg2FCKqMEmFRhtJUcAgL04cnR8fYON4POBKVvsYJMKizaCLaZcPzeOIgXW4yS5/XlaGs2JYBtGEfu0EWsAplTD/mzEWllyTtk8BFspwk/aCLYy8HMEfBhrukYEIqcSpmxS8GjjiR59jb1EysEIiISwGpaXpdHap2ylBJsUeHmUYCsPJ9xFBCUujRJrGAPBNryQk7ZL9nEKgVIR49gGR472N2WwEYfTy33EaE6lBZsUNNpKdtmV38q9CdXhKFLUiCNDx+VlaTSX3FO2EoNNCrw8WpLFUFsWbqeefeFUm5Mdy6TAuHLG2kl7QkfSddm+fk1yTdkItvEQbZlsO1FbvH5NEYe9Dh84r2fOMiifWvT92ZisoXQlB5tEtE1qyKXP0iNu6mnbmK9FevXBs/r6mQOjfG4A/nhYGq11ylY6om1kU+2jVmLEsUwKDMfTlI0l0jLlDrbSp2wS0TYKDwcT1L5fHHwY88Xipz5yNPrSKOpR4xGjnoPt8IHzg30uom0gHkJtk4jTOKZtQH+epmyl87A0mkPOKVstwSZxyo9BRAg2APCEkAS2R7T1FDnYIm87MKXIS6PEEVAOoq0HogcAdkdQomRDL41KRNvOCDYA3hFFQFmIth2UFGwlfS81qf3EurzmKLyKvJQO/+p+5N8BkQN0M+bpPqYU9Uk4ypQtynYCHhBtWyDYUKqTZ3JvAQBsz/PpPsZAtHVEsAGIItr0Ktr2elbjiXVrQrQBI+N1R4c15f5sUZdGAZSJaOuAKVs+vBoCsJ2oU6uo2w0sM8bpPiSibaMagq2G7xEAMI7cLxRfE6INAAoRfVoVffuBsRFtazCBAnZTwuk+2J8NgDdEGxBMzSfW5aS6qzGlAspX76P/BkzZUAvO0QZPiE9gNaINABZEWxoldIA6EG0AAFeIUHRR26shSETbUiyN1ueZMxX+9qMIBA7gy1jnaJOINgADG+vI0akOQoi2NFoqYhTYi2hbUOuUrdbvG4iMsAHqQrQBI+J1R4HdEaXAxYg2AAiIoPGLJXaMhWhrYYkQ3g19Yl3O0XYxnmz9IU6BCwZ5BjCzW83sS2Z23MzuXPL+N5vZM2b22fmft7bed4eZPTb/c8cQ2wOgLLwSwsUIGaBOl/T9BGa2X9J7Jd0i6YSkh8zsWErp0YWrfiil9LaFj71C0q9KuklSkvTH84/9Rt/t2hZTNgAA4NkQk7abJR1PKT2eUjor6T5Jt3X82NdJejCl9Nw81B6UdOsA2wQggxJeKN67GqdsNX7PiGnMc7RJw0TbtZK+1nr7xPyyRT9lZp8zsw+b2cu2/FgAGBX7swFx1PhqCNIw0WZLLlv87/a/l3RdSumvS/qYpHu3+NjZFc2OmtnDZvbws2dP77yxABBVzROnmr93oDFEtJ2Q9LLW2y+V9FT7CimlZ1NKzZ7E/0bS3+j6sa3PcXdK6aaU0k1XHjg0wGZf7NSzLxz8c0bi8ft/9i+H/zkvM9ZLWHk/R1uUI0c5CAENwm2zWn5fvD5+DX2E/6IhPvtDkm4ws+vN7ICk2yUda1/BzK5pvfl6SV+c//sBSa81s8vN7HJJr51floXHcJlCrd93NGM/GPTF/mzjIlhmotwONS235/rdrzHceh89mlI6Z2Zv0yy29ku6J6X0BTO7S9LDKaVjkv4HM3u9pHOSnpP05vnHPmdm/0yz8JOku1JKz/Xdpj5OPfvCKo4k9R5q0adsiMX7E2yUUJnKSXtCR9J1uTfDrZNnnteRg5fm3oxJnDzjc/+2Z87uG+WgBEsp3v+O//qLj6Rj/81Pjfo1Sgw376HWmCrYpDhLoxFOqjvG/7Z5kfgZom25COF2RJfl+boZou3woWW7qY/PY7Q1mnB785/8T3+cUrqp7+fzvd6SUZTA2eTUsy/87p8ISgi2GrE0Oh6CbTVum9Vy7NvGMulez5zdN+h/uom2NaIFTyPqdpfC+wEIiIMo2cz7bZRzilvLQQmS73AbUu992mqxLIA8LaGWEGhM2ZaLsDQ6BpZG0ZX3fdxO6lS2ZdKpPXM6ZVsm9bp/25CIth5yh1wJodYoJdhqnLKxNDoO7xMkxFDTQQk1INoGNnbIlRRqjSmDLRrvp/nAOAi27TFt84Np23iItgn0DbkSQy2X2qdsUZZGp8LSKGqQY9pGuI2DaMtkU8jVEmpM2VaLMGUba2m0ph2oFzFl2x3TNjRKDTf/zwoVqe2oz6mDrfYpG/wj2PrjNlyuplOAlIxoQxZM2KbH0ujFWBpFDrXd73KGW4mPeUQbqhBtyhZhaXQstS6NMiEaDrflcjX+bpUWbvU+MyAbpmzTi/KyVbUiMurCtG1aJYUb0YZJ5Qg2pmxAfQjh5WqctknlhBvPDpjEs395qLhgQ1yeJh3ExXg837a1vbxV7mlbKYg2jK7U5dAoR4xGWhqtbQrgOSqAoeUOtxKmbUQbRpUz2CJO2VgaBYblOYxrm7Z5ED3ceIbAaEoOtihTNuzlaWkUqE3uaZsUO9yINoyi5GAbyxhTNpZG/fI8ASqN59u6xmkb4bY7XsYKgyp1/7U2pmwAEF/El7pi0obBeAi2qFM2TMPL0qjnyU+pPN/mTNvyiTZxI9owiFqCbawpW+1Lo8AUPIcb8okUbkQbeiPY0Fdt+7MBi5i25RUl3Ig29OIh2KLjNB91YdqTl+fb38vy/ZQIt+3wbIGd5HqFg2WYsu3F0uheNT4hAl0x7Z7xHm4ho+3ceXMVDTXxdrtHDzambDxZYHpM21Z8bZZJJc3CzWu8hT/lRzsgrvye0xm3pFyeIq0x1VGiEYONKdteXqZsnmOhNiftCR1J1+XejKVO6pSO6LI8X/vM8zpy8NLJv27zGHP4kE3+tVdpHks9nRakqP/mN1Mgb9OgqDzejs+cOUiwrREt2JiyISfPEV3jxE3y+R9ET5O38JO2dRaDg0ncZt4irW3Kc7BF24cNqzFlwzrNz8Xj1K3GiZvkc+om+TgZb9HRtoil1OU8h5o0/Qlzxw42pmwzTNngidfl0tzhJilrvHkMNylfvFUVbW01TuG8x9miHK9uQLBd4HGZIiqmbDEQbiu+PlO3PXLFW7XRtqiUiIsWZquUNl2TYgXb2KaasnlZGkUchNuKr58x3CSfUzdp+ngj2lbYFD+5o66UOFtU4nRNindqD6Zsw2HKFg/htuLrOwg3yd/UTZou3oi2HQ0RTZvCr9QwWybXC71HDzambGu+DlM29EC4rfj6mfdzk/xO3aTx441oy6imKFuHYNvNWMHGwQfDYcoWG+G2ZhuYuq01VrzFWrNBUaY859oigm25UpZFmbJhKF7D+6ROZb+fe/hPmPfHrKHP8Ua0YXI5Y00i2HLx8AA/Ja9P9tie558l4TYLN+/xNhSiDZPxEGucNHe1Uh70cj+JoUyE25qvf+Z5N/FWOqINo2pCLWesSdO+wkHEKdvYD3YeHtCn5PkJHrs7aU+4/dnmDjfJx+956VM3og2D8xJqDYKtHh6euFA+wm3NNjgIN6ncqRtHj2IQXgKtrZRYk8YNNqZsw/L6hI5hcWTpmm3IfGRpw/sRprtg0oadeZuotRFs3ZT0v1EPUwbUxWuge/hd8LKfm1TW41zIaDuXQm52ETyHmjT9wQaRg20KXh60p+L1SRzj8foz9xBukp/HgFL2dQu7PLosGg4fDP4M55TXQGvLcVRo9GAr4QGs4eUJCnViqXTDdjhZLpV8v5pCF2GjbZnFuCDidhch1KQyY60UXv6HPRWvExdMg3DbsB3Owk2Kua9bUdG2iIjbTpRQk/LEmjRdsDFl644pG7zwHG6Sssebh9ctbYsYb0VH2yKWVC8WKdLaCLZ+pgg2pmyolddwk5i6rRIp3qqKtmVKnsZFjbJVcr6aQSnBNoUpg40pGzwi3Dpsh7Nwk2LEW/XRtijSNK60KFsl90tPlRRsJS2LesGUDcsQbh22w2G4Sb7jjWjrYIw46hqCtYTZMrXEmlROsDFlAy4g3Dpsh7P93No8xhvRlknNMbZJ7liTygs2jIMpGzYh3LrxOnWTfMUb5y6AG1OfGHeVEk/pwZQNyMdz3Hv6XfL0KgrLeDhBL5M2ZOch1KTpY22qCVvuB5lSeX4ihj9M3LrzvGQq5Z28EW3IwkuoSXkma6UFG1M2YLMm9D3Gm5dzubV5XjKVLn58nSrgiDZMwlOkNUqONYlgGxNTNvTB1K0771O3xlTTN6INg/MYaG259lkj2MpAsGEITN22Ey3epHECjmhDL94DrS3nAQYE2wBfz8mUDRiS96mbRLztaozpW8ho+/Z5uygWrj54NuPW1CNSoLURa8OrNdiYsmEMnsNN8rdkKsWMtyGEjLZF62KCoNtN1EBry33qDoJtoK9HsIVz6uyTF7192YGXZ9qSODwvl0o+p25SrHgbQhHRts6m+CDqygi0Ru5Qk6Y/WW7JwYY4FkNt1fsIuPUiTN0k4i2X4qNtk9KmdCUFWBceIq2R45UNSg82pmz+rYu1Tdcn4JbzPnWTiLdcBok2M7tV0r+UtF/S+1NK71p4/y9Jequkc5KekfQPU0r/af6+70h6ZH7VJ1NKrx9im4YwVQCtisPaAqwrT6EmlR9rEsGGi20bal0+DwG3l/epm+Rzfzep3HjrHW1mtl/SeyXdIumEpIfM7FhK6dHW1f4/STellL5lZr8g6V9I+pn5+06nlG7sux2REWebeQs1Kd9rhhJsyGGoUOvy+Qm4C5i69VNavA0xabtZ0vGU0uOSZGb3SbpN0nejLaX0ydb1PyPp5wb4uigcoXaxHC9HVXuwMWUbP9a6fE0iLs7UTSLexjREtF0r6Wutt09I+pE113+LpI+23j5oZg9rtnT6rpTS/zXANiEoj6EmEWuTfV2CzYUcobYOU7iZCFM3ye+SqRQ/3oaItmVnjVv6LGNmPyfpJkk/1rr45Smlp8zsFZI+YWaPpJS+suRjj0o6KklXvuDF/bcabngNNSlvrEl1BRvy8xZryxBwMeLN89RNihtvQ0TbCUkva739UklPLV7JzF4j6R2Sfiyl9N1nhZTSU/O/HzezT0l6laQ90ZZSulvS3ZJ0/QtfMv0zGQblOdSkOmNNyhtsTNnyiBBqq9QecCyZ9hct3oaItock3WBm10v6U0m3S/rZ9hXM7FWS3ifp1pTSydbll0v6VkrpeTO7StKPanaQAgrhPc4W5Y41iWDLrYZgixxqqzTfU23xFmHqJvleMpUufvzzHHC9oy2ldM7M3ibpAc1O+XFPSukLZnaXpIdTSsck/a+SXiTp98xMunBqjx+U9D4zOy9pn2b7tD269Au1nEvDvwgr+okWZ20eQk3KF2sSwdYoPdhKjLVFxNt1WbdjHe9Tt4bn6dsg52lLKd0v6f6Fy/5p69+vWfFxfyjplbt8zU2RcPjA+V0+LTaIHGeLiLUZgq0ONQRbW83x5jncpHjxJvkJuGJfEaFLXBB2q5UUZ4u8xJpU73RN8hdsJU/Zagu2thrjLcLUTfK/ZNrmZfpWbLR1Ueu0ruQgW4dYuyB3sAFTI958ijJ1a+SOt6qjbZMx42abIKw1sobgKdSk/LEm+Qg2pmzTqXnKtkyt8eY53KS48SZNG3BEWyaE2Di8RVobwTbfBoJtMgTbarXFW4SpmxQv3qRpp29EG8LzHGqSj1iTCLbaEGzd1Bhv3sNNuvixIkrATRFvRBvC8R5pbQRbaxscBlupUzaCbXs1xVuUqVsj2vRtzKVTog3uRYq0BrF2MYJtOgRbP8SbX9HiTRr+MZhogzsRI63hJdYkgg3oo7Z4ixJuUsx4GwrRhuwiR1rDU6xJfoLNK6Zs6KqWeIs2dZPqjDeiDZMqIdAa3kKt4SnYPE7ZCDbsgnjzq6Z4Cxlt3z6/+cn/yMFptgUXlBRky3iNtAaxthnBhr5OnX2y+HCTiDevQkZbF10DgrjrpvQgW8V7qEm+Yk0i2KZGsE2vlqmbFG9/N6nseCs22rrqEyMlBF+tMbYOobYbr7EmEWwYRy3xFnHqJsV6bdOuqo+2PgieMkSItIbHWJMIthwINj+IN79Km7oRbagSoTYMYi0Pgs2nmvZ3ixRuUjnxRrShCpEireE51iSCLReCzTembr5FjzeiDUWKGGmS/1BrEGx5EGxx1DR1k4i3qRBtCC1qnC0i1oZBsMGTWqZuEvE2FaINIZQSZ21RQk3yH2sSwQa/aou3aOEmxYm3kNF27nx5T+CYKTHOFhFrwys52FAOlkz98x5vIaNN2u7J/fAhG3FLsIsa4qwtUqhJcWJNKj/YmLKVpbapm0S8DSlstG2DwBtXbQHWVbRQk2LFmkSwIa7a4i1iuEkXPyZ6CLgqom0bQwZIxAAkwPqLFmvRQk0qP9Ykgq0WLJnG4WH6RrSNiAAqX7RAa4sYaxLBhvLUNnWTiLddEW3AFiJHWiNqrEkEG8pGvMWSI96INmCFEgKtLXKsSQQb6lHLkqkUe3+3xpTxRrQBc6VFWiN6rEkEG+rD1C2eKeKNaEOVSg20thJiTSLYUDfiLZ4x441oQxVqiLQGsRYLwYYuaow3KXbAjRFvIaPt2zuMIcQAACAASURBVCld9CR85OClGbcGudUUZOuUEmsSwQasUlO8SWVM34Z8bA4ZbYu2fdIm8mIhylYrKdSkemJNItjQD/FWpyKibVtEnh8E2fZKCzWprliTCDYMh3irS5XRtq0xwiJyCBJa0yox0tpqCzZgDDWdJkSqN96ItkwIH6xSeqQ1ao01pmwYS21TN6m+eCPagMxqibRGrbEmEWyYBvFWLqINyKC2UJPqjjWJYMP0iLfyEG3ABGqMtDaCjWBDPsRbOYg2YGC1B1qj9lCTiDX4UtvBClIZr23aRrQBOyDMliPUZog1eEW4xRYy2s7pO0ufNMd8kVbUhzDrjli7gGCDd4RbXCGjbZVtnmQJPBBl/RBqFyPWEAn7ucVUVLRtY5cnbEIvFqJseITaXsQaImPqFku10baLISKA8NsO4ZUXkbYewYYSEG5xEG0TyxkhfYKReKoDkdYNsYbSEG4xEG0VIbywiEjbDrGGktUablKc/dyINqAyhNpuCDbUoMZwk+JM3Yg2oHBEWj/EGmpDuPlFtAGFIdKGQayhZoSbT0QbEByRNjyCDajzXG6S7/3cQkbbOT3v+kYFhkSUTYdYA/Zi6uZHyGhr2/YJzdsPAJAIs5wINWAzws2H8NG2LSIPuRBmfhBqwPYIt/yqi7Zt9X2i9fKDxrgIMv8INaC/msNNyv+cTrSNbIgn89x3kloQXuUh1IDh1XqAgpR/6ka0BZAjJjyGIlGFLgg1YBo1T91yPUcSbViKQEI0xBowPcJtWkQbgLAINSC/msNNmnZlat9kXwkABnDq7JPf/QPAh5p/H6dcmWLSBsC9mp8QgCg4QOG60b/OIJM2M7vVzL5kZsfN7M4l77/UzD40f/8fmdl1rfe9fX75l8zsdV2+3rl0dojNBuBUe5pGsAGx1Po7e9KeGH3q1nvSZmb7Jb1X0i2STkh6yMyOpZQebV3tLZK+kVL6L8zsdknvlvQzZvZDkm6X9NckvUTSx8zs+1NK39n0dbe9U9RY/oB3tT64A6WrdT83adyp2xDLozdLOp5SelySzOw+SbdJakfbbZLeOf/3hyX9azOz+eX3pZSel/RVMzs+/3yfHmC7LkLkAfkQZ0B9CLfrBv+8Q0TbtZK+1nr7hKQfWXWdlNI5M/tzSVfOL//MwsdeO8A29bbrk0ytd1CgQaABaBBu1w36OYeINltyWep4nS4fO/sEZkclHZ39+wXbbN+kiD3UgjgD0EXt4TakIaLthKSXtd5+qaSnVlznhJldIunFkp7r+LGSpJTS3ZLulqT9+w8tDbvIiD14QYwBgE9DRNtDkm4ws+sl/almBxb87MJ1jkm6Q7N91d4o6RMppWRmxyT9n2b2v2l2IMINkv7jANtUDWIPmxBhAFCG3tE230ftbZIekLRf0j0ppS+Y2V2SHk4pHZP0m5J+Z36gwXOahZ3m1/tdzQ5aOCfpF7scOYr+hn4iJwLHRXgBiIrnh+FYSvFWGvfvP5RedOiv5t4MONT3wYE4AoBhEW3Sk9/46B+nlG7q+3l4RQQUhegCAJSK1x4FAACjYMo2LKINAAAgAKINAAAggJD7tJ0/f1bfPL1536UXHWIsCwBADiyNDi9ktHVF2AEAgFIUHW1ddAk7ibgDAAB5VR9tXXWNO4nAAwDUjaXRcRBtI2B6BwAAhka0ZbTN9E4i8gAAqBnRFsi2kScRegCAabE0Oh6irXC7hJ5E7AEA4A3RhqV2jb0G0QcAwLCINoyib/Q1iD8AiIOl0XERbXBtqPhrEIEAgKiINlRl6Ahch0AEAAyJaANGMmUg7oKoBDAklkbHFzTazuv8+dNL37Nv36GJtwWIyXtUDoU4BVCKoNG22qqYW0TcAXXIEaeEIoAxFBdtXRF3AMbC+RFRG5ZGp1FttHVF3AGYyrrYI+gAEG0D6Rp3EoEHYHsEHbxiyjYdoi0DAg/AkDYtxxJ1GAvBNi2izblNgUfUAdhkMeqIOAyBYJse0Rbcuqgj6AAsQ8ShD2ItH6KtYEzpAHRBxKELYm03R9J1GurEQ0RbxbrsW0fYAfUh4rCIYNvNkXTdoJ+PaMNahB0AIq5exNpuho61BtGG3gg7oC5EXB0Itt2MFWwS0YaJbHOaE4nIAyIh4spCrO1mzFhrEG1wiXPZAXE1EUe8xUOw7WaKYJOINhRg2ymeROgBU/jm6ScJt0AItu1NFWsNog1V2iX0JGIP2BbhFgPBtr2pg00i2oCt7Bp7EsGHehFuvhFs28sRbBLRBkyG4EPNCDd/iLXt5Yq1BtEGBMByLkpAuPlBsG0vd7BJRBtQNA7SgDeEW34E23Y8xFqDaANwEc6ph7ERbvkQbNvxFGwS0Qagp66RR9yhjXCbHsHWnbdYaxBtACbBCZOxiHCbBrG2Ha/BJhFtABxaF3gEXVkIt3ERbN15jrUG0QYgFIKuPITbOAi27iIEm0S0ASgIQRcX4TYsgq2bKLHWINoAVGFV0BFzfhBuwyDYuokWbBLRBqByy2KOkMuHcOuHYNssYqw1iDYAWLAYckTctAi37RFr3UQONoloA4CNmMZNj3DrjmDbLHqsNfbl3gAAiOj8+dMX/cHwvnn6ydyb4B7BtlkpwSYxaQOAQbCkOg4mbqsRbOuVFGsNog0ARkDEDYdw24tgW63EWGsQbQAwASKuH8LtAoJtuZJjrUG0AUAG7Ygj4Loh3LBMDbHWINoAIDOmcN01ByfUGm9M2S6oKdYaRBsAOMMUbrMap24E20yNsdYg2gDAMQJutRrDrWY1x1qDaAOAIAi4vWoJt1qnbITaxYg2AAiI/eAuqCXcakKsLUe0AUABmMKVq6YpG7G2HtEGAIWpMeCYtsVGrHXT67VHzewKM3vQzB6b/335kuvcaGafNrMvmNnnzOxnWu/7gJl91cw+O/9zY5/tAQBcrKbXR+W1SmMi2Lrr+4Lxd0r6eErpBkkfn7+96FuS/n5K6a9JulXS/25ml7Xe/49TSjfO/3y25/YAAFaoIdxKU/LS6JF0HcG2pb7Rdpuke+f/vlfSGxavkFL6ckrpsfm/n5J0UtLhnl8XALCD0qduTNtiINZ20zfark4pPS1J87+PrLuymd0s6YCkr7Qu/rX5sul7zOzSntsDAOig9HiDXwTb7jYeiGBmH5P0fUve9Y5tvpCZXSPpdyTdkVI6P7/47ZL+TLOQu1vSr0i6a8XHH5V0dJuvCQBY7/z508UdrFDKQQmlLY0Sa/1tjLaU0mtWvc/Mvm5m16SUnp5H2ckV1/teSX8g6Z+klD7T+txPz//5vJn9lqRfXrMdd2sWdjKztGm7AQDdNBO30uINfhBsw+i7PHpM0h3zf98h6fcXr2BmByR9RNJvp5R+b+F918z/Ns32h/t8z+0BAOyopCVT9m3zg2AbTt9oe5ekW8zsMUm3zN+Wmd1kZu+fX+enJf0tSW9ecmqPf2tmj0h6RNJVkv55z+0BAPREuOVXytIowTYsSyneSuNseZTzAgPA2KIvmUbdty16tBFrF3v41K//cUrppr6fp++kDQBQsOhLphGnbQRbWY7oss1X6ohxFQBgoxKPMsXwCLaZIUOtjWgDAHQS9SjTUk4B4hmxNl6otRFtAICtRI23CCIujdYcbFOEWhv7tAEAdhJpX7eI+7ZFUGuwHdFlkwebxKQNANBDpKkby6TDqi3YckTaIqINANBbpHjzKsrSaE2x5iHU2og2AMBgvB9lyrStnxqCzVuotRFtAIBBMXUrU8nB5jnU2og2AMAovE7dmLZtr8RgixJqbUQbAGA0XsMN9YoYaw1O+QEAGJXHU4NwCpDuSpmy5TpNx5CINgDA6Ag35BQ91hosjwIAJsFSaTzRp2ylxFqDSRsAYDLeJm5M21Yj2Pwh2gAAk/IWbihPicEmEW0AgAw8hZuXadupsz62Q4o7ZSvhYIN1iDYAQBbnz592FW+YiRxspSPaAABZeQg3L9M2bK/06Vob0QYAyM5DuCGeWmKtQbQBAFzIHW5M2+IsjdY0XWsj2gAAbhBu2KTGWGtwcl0AgCuchBfL1BxrDSZtAAB3ck7cap22eV4aJdhmiDYAgEu5l0rhA8F2AdEGAHArV7jVOm3zpNaDDdYh2gAArjFxG5+3pVFibTmiDQDgXo5wY9o2PaZr6xFtAIAQanjZK0+vPzo1Ym0zog0AEMqU4VbDtC330ijTte44TxsAIJwpz+X2zdNP6kWHXj7J10I5jhy8dPDPSbQBAELiJLzxlTZhGyPU2lgeBQCENdVSaanLpDmXRksKtiMHLx092CQmbQCA4Ji4xVNCsE0RaYuYtAEAwpti4lbqtA3bmWqqtgyTNgAAKpRraTTilC1XpC1i0gYAKALTNgwt51RtGaINAFCMEk6+W/IJdqNM2bzFWoNoAwBgC0zbyuU11hrs0wYAKApHk2JbnkOtjUkbAKA4Yy+TMm3bnselUe+TtUVEGwAAzpS8X5sH0WKtQbQBAIrEtA3LRIy1BtEGAChWCUeTlsDL0mjkYJOINgAAdsa0LY7owSYRbQCAwkWdtrFf23BKCDaJaAMAVCBquJXAy9JoCYg2AAB6YInUt1KmbBLRBgCoBNO2+pQUbBLRBgCoyFjhNta0Lfp+bTmXRksLNoloAwAAhSkx2CSiDQBQmWjTNqBBtAEAqsP+bePLtTRa6pRNItoAAEAhSg42iWgDAFRqjGnbGEuk0Q9GmErpwSYRbQAAYGBTL43WEGwS0QYAqFiUaRsgEW0AgMpxUEJstUzZJKINAIDBDT1ti7Rf25RLo96D7fAh0+FDNtjn6xVtZnaFmT1oZo/N/758xfW+Y2afnf851rr8ejP7o/nHf8jMDvTZHgAAdsG0LR7PwTZ0rDX6TtrulPTxlNINkj4+f3uZ0ymlG+d/Xt+6/N2S3jP/+G9IekvP7QEAYCeEWxxeg22sWGv0jbbbJN07//e9kt7Q9QPNzCS9WtKHd/l4AAA8q/GAhJyvNZrbmLHW6BttV6eUnpak+d9HVlzvoJk9bGafMbMmzK6UdCqldG7+9glJ1/bcHgAAdsa0zT9vU7axp2ttl2y6gpl9TNL3LXnXO7b4Oi9PKT1lZq+Q9Akze0TSXyy5XlqzHUclHd3iawIAsLXz509r375Dg3yub55+Ui869PJBPhd8BtuUNkZbSuk1q95nZl83s2tSSk+b2TWSTq74HE/N/37czD4l6VWS/p2ky8zskvm07aWSnlqzHXdLunv+dVfGHQAAJTp19klddoAA9GLqYJP6L48ek3TH/N93SPr9xSuY2eVmdun831dJ+lFJj6aUkqRPSnrjuo8HAGBqLJNinRzBJvWPtndJusXMHpN0y/xtmdlNZvb++XV+UNLDZvYnmkXau1JKj87f9yuSfsnMjmu2j9tv9tweAABcqfGAhDF4WRrNFWxSh+XRdVJKz0r6iSWXPyzprfN//6GkV674+Mcl3dxnGwAAAKaQM9gkXhEBAIClWCJFW+5gk4g2AADC8PxyVmOeoy330qiHYJOINgAAVhpq2sZ+bXF5CTaJaAMAAFjKU7BJRBsAAHAs19Kot2CTiDYAANZiibQ+HoNNItoAAAhlqIMRTtoTg3ye0ngNNqnnedoAAKjBkK9Hiu6mXBr1HGsNJm0AAKBqEYJNItoAAOhkiH3b2K/NnyjBJhFtAADAoSmWRiMFm0S0AQCACkULNilstO3Tvn2HVv4BAGAMXpZIPb+cVQQRg00q9OjRPuHGCwQDALCdoV93dMyl0ajBJhUabX2sCz6CDgDA6T+QC9G2hcVfUiIOAIA4Ik/ZpLD7tPnAfnQAgF1w6o/Vcr3WaARE20A4GAIA6uFhpYWDEbYTfcomsTw6mna4efjlBgAAsTFpmwATOADAIpZI9xpraTTnlO3IweE+F5O2iXEwAwCUgaNIscmQwSYxacuOKRwAIJeT9kTuTRhdrinb0MEmEW2uEHAAgG1wMIJPYwSbRLS5RbwBQPnYr+0CTvWxGdHmHPEGAH6xX/LwL2E1pBxLo2NN2SSiLQziDQAA38YMNoloC4d4AwCgTkRbUMQbAKCvI+m63JsgaZz92aZeGh17yiYRbeERbwAA5DVFsElEWzGINwAAysYrIhSG1zwFANRuyqXRqaZsEpO2ojF9AwBgPFMGm0S0VYF4AwB4FfWkulMHm0S0VYV4AwAgLqKtQsQbAKBUU+zPlmPKJhFtVSPeAKA/DvqqS65gk4g2SIQbAGAnnl93tEREGyQxdQMAxDf20mjOKZtEtGEB8QYAwF65g00i2rAC4QYA0/jm6Sdzb0I2UU/3kQvRhpWYugEAohhzadTDlE0i2tAB4QYAqJWXYJOCvvbovn0H9KJDL1/6vprHzGNqwo1D2wEAyCNktK2zKuYaRF0/+/YdItwAAK6MtTTqacomFRhtm2yKOomw24SpGwCgdN6CTaow2rpg6bUb4g0AsCvPR456DDaJaNvKYswRcTMsmQIAcpnitUa9INp6aEdc7QHH1A1Azc6fP82R9oXwOmWTiLbBMIWbYeoGAHUo8XVHPQebRLSNpuYpHFM3AMAUaloalYi2SdQ6hSPeANSCpdH4vE/ZJKIti9oijngDgNW6nIqqNJ6PHPWMaHOglogj3gBgWJcd2D34jqTrhtuQDGpbGpWINpdK3x+OeANQkhqXRks7CGHMpdHDB84P9rl4wXjnXnTo5cWOzvftO1Tlgx0AALtg0hZEydM3Jm8AalXqf8rX8bY/W5Qpm0S0hVRqwBFvAKKJulrA/mwxsTwaXInLpyybAsBmfQ5CwEykKZvEpK0YJU7fmLwBgD+lHYQwhjGCTeo5aTOzK8zsQTN7bP735Uuu8+Nm9tnWnzNm9ob5+z5gZl9tve/GPtuDmWb6VsoEjskbAI/6Pi6V8hg9taGWRiOcTHdR3+XROyV9PKV0g6SPz9++SErpkymlG1NKN0p6taRvSfq/W1f5x837U0qf7bk9WEC8AQC88HYQwhjGmrJJ/aPtNkn3zv99r6Q3bLj+GyV9NKX0rZ5fF1sqafpGvAGoXc0n1R1CxCmb1D/ark4pPS1J87+PbLj+7ZI+uHDZr5nZ58zsPWZWfoI7UErAEW8Aoor6+Mv+bOuNOWWTOhyIYGYfk/R9S971jm2+kJldI+mVkh5oXfx2SX8m6YCkuyX9iqS7Vnz8UUlHZ/9+wTZfGms0DxyRD17Yt+8QBysAmBT/YcxjiP3Zok7ZpA7RllJ6zar3mdnXzeyalNLT8yg7ueZT/bSkj6SUvt363E/P//m8mf2WpF9esx13axZ22r//UNq03dhO9HjjSFMA8K30/dnGnrJJ/ZdHj0m6Y/7vOyT9/prrvkkLS6Pz0JOZmWb7w32+5/agp+hLpyyZAgBK1Tfa3iXpFjN7TNIt87dlZjeZ2fubK5nZdZJeJun/Wfj4f2tmj0h6RNJVkv55z+3BgIg3ANgr96k+aj0IwfOrIEwxZZN6nlw3pfSspJ9YcvnDkt7aevsJSdcuud6r+3x9TCPyiXtZNgWAYZRwEELk/dkkXsYKW4o6fWPyBmAIPI5g0VRTNolow46INwBAV6UfhDAVXnsUvURdOmXZFEAOOfdni8zrqT6mnLJJTNowoIjTNyZvALqK/lgR+SAEzDBpw+AiTt+YvAHAaiUchDC0qadsEpM2jCza9I3JG4CxRHosHFLf/dk8Lo3mCDYp6KTtEjuwcV3/1NkYE55aRHvFBSZvANo8/Geu1v3ZcEHIaOuCqPMp2tIp8QagBOzPNpxcUzap4GjbpMv/WAi7cUWavhFvQL2GmLJFXhqNvj9b9BPqtlUbbV2sCzuCbjiRpm/tB28CDiifh2VRKe7SqIf92UpCtO1o2S8QIdcf0zcAXgwVbLmnbCyNDifn0qhEtA1q1f+EiLntMX0DkJOnYMv5AvF9lkY9TNlKWhqViLZJEHP9RAw44g2Iq5Rgw7ByT9kkoi0rlli3F2X5lOkbEJOXfdiGkHPKhnEQbc4Qct1EnL5JBBzg1dCxVvuUjaXRcRBtASz+4hJxFyPgAPRRYrAxZRuWh6VRiWgLiWncahEDjngD8vEYbMAqRFshmMbtxf5vANbxuv9a7ilbX6UtjXqZsklEW7GIuAsiTt8kAg4Y0xjB5mFZdAgsjfpFtFWCJdWZiAFHvAHD8hpsQ4g+ZcN6RFvFap/GRQk4pm/AcLwuiUpM2SSWRjch2vBdNUccAQeUb6xg87IsmnvKhvERbVip1ogj4IDyeA42L/pO2TwcgDAkb1M2iWjDFtoRR8D5QsABq3kPNqZswynxhLptRBt2UuMUjoAD4vG8D5vkYz82Kf+UDd0QbRhEbRHHOeAA/8YMNk/LoiVM2Vga7YZowyhqibgo0zeJgENdIgRbKVM2L0pfGpWINkykhoiLGnASEYdyjL0c6i3YPEzZSlsa9Tplk4g2ZFJ6xC0+sEeKOAIOUdUWbEPwMGXzdm62oV198Oxgn4togwulH5kadQpHwCGKCMuhQ2PKVh+iDe7UNIUj4ID+IgWbpykb9hp6aXTIKZtEtCGAkiMuasBJRBzyi7Ic2hgy2IaYsrE0Gg/RhnBKjbhIAScxhUM+U5x7zeuSqCelLY16PgChQbQhvBL3h4sccBIRh3FMdaLcMYKNKVt9hl4alYg2FKbEKVy0I1ElIg7DmvJVDWoINi9YGt0e0YailR5xEQJOIuKwvRwvP1XLkugQUzaWRvMg2lCV0pZSIwacRMRhtZJibegjRUtaFvX2slVDG2NpVAoabZfo0pV33pP2xKTbgrhKm8JFXEZtEHF1y/Wi7mNO1jwG21C8TNlqWxqVgkbbOtvcsQk8tDGF84OIK1+uUGvUGGxelkW9TdmiLI1KBUbbNrr+IhB39WEK58uqJ3hiLpbcoSaNv99aycGG/KqOtq6Y3oEpnE/LIoCQ86eGWJPKDzZPUzbPS6Nj7c8mEW2D2/aXjMiLhymcb4ScDx5CrUGwYZVIS6MS0ZYdS7TxMYXzj+XV8XiKs0VTncLDa7ANiSmbD0RbEMRdDEzhYiHmuvEcZstMeb41z8HGlK08RFthuvzCE3bTIeJi6hoppcVdtDhbJup0TfIbbJ6mbN6NuT+bRLRVadMDA1E3HiKuLFHjroQ4W8R0rfX5Cp6wDbk0Gm1/NolowxLrHkAIumGVvD9co7aQW6bESPJg6ped8j5dGwNTNl+INmyFoBtPaVO4Ru3TOAwrx+uDjhFr0jjB5m1ZdEg1H4DQINowGIJuWLVEnETIYbWcL+I+VqxJ/oNtKF6nbBGXRiWiDRPhtWL7KzXiJKZxuFjOUGvUHmzepmwRjH0QgkS0ITNibnc1RZxEyJXOQ6hJ48aaFCPYhjLklI2l0RmiDS4Rc9srOeKk1U/qxFxcXkJNGj/WJP8HHTRKn7JFXRqViDYEQ8x1V3rENYi5WDyFmhQ/1rwui3rdly06og1FWPagSMhdbNmTU6khJxFznngLNWmaWJNiBZtXEZZGp9ifTSLaUDCmcpvVMo1rY3+58XmMtMZUsSbFC7bSl0Wl2EujEtGGCjGVW63GiJO6RQZhd4HnKFuHYFvzOQcMNg5AGA/RBoip3Cq1RtwyXUOlhLiLGmWrTBlrUpwDDho1TNhKQbQBazCVu1ht+8XtIkLclRZlq0wda9L4weZ9P7YaD0CYan82iWgDtkbIXWzVEyMxt14t4ZRDjliTYgab5ynb0Euj0fdnk3pGm5n9PUnvlPSDkm5OKT284nq3SvqXkvZLen9K6V3zy6+XdJ+kKyT9v5L+u5TSdMkKDISQ24upHKZWaqxJMYKtxinb1PpO2j4v6b+V9L5VVzCz/ZLeK+kWSSckPWRmx1JKj0p6t6T3pJTuM7PfkPQWSb/ec5sAFwi5vQg5jCFXrElxg807DkBYrle0pZS+KElma+v6ZknHU0qPz697n6TbzOyLkl4t6Wfn17tXs6kd0YZiEXJ7EXLYVs5Ia4scbLVN2UpYGpWm2aftWklfa719QtKPSLpS0qmU0rnW5ddOsD2AK4sP/LVHnMRRq9jLS6g1CDZI0x6EIHWINjP7mKTvW/Kud6SUfr/D11iW32nN5au246iko/M3n3/41K9/vsPXrs1Vkv5z7o1whttkOW6X5bhdlst+u/zFX34h55ff40kHt4lT3C7L/cAQn2RjtKWUXtPza5yQ9LLW2y+V9JRmP9TLzOyS+bStuXzVdtwt6W5JMrOHU0o39dyu4nC77MVtshy3y3LcLstxu+zFbbIct8tyZrb0QM1t7Rvik2zwkKQbzOx6Mzsg6XZJx1JKSdInJb1xfr07JHWZ3AEAAFSnV7SZ2d81sxOS/itJf2BmD8wvf4mZ3S9J8yna2yQ9IOmLkn43pdTMuX9F0i+Z2XHN9nH7zT7bAwAAUKq+R49+RNJHllz+lKSfbL19v6T7l1zvcc2OLt3W3Tt8TA24XfbiNlmO22U5bpfluF324jZZjttluUFuF5utUgIAAMCzKfZpAwAAQE9uo83M/p6ZfcHMzpvZyiNRzOxWM/uSmR03sztbl19vZn9kZo+Z2YfmB0GEZmZXmNmD8+/pQTO7fMl1ftzMPtv6c8bM3jB/3wfM7Kut9904/XcxvC63y/x632l978dalxd3X5E6319uNLNPz3/XPmdmP9N6XzH3l1WPE633Xzr/2R+f3xeua73v7fPLv2Rmr5tyu8fW4Xb5JTN7dH7f+LiZ/ZXW+5b+PpWgw+3yZjN7pvX9v7X1vjvmv3OPmdkd0275uDrcLu9p3SZfNrNTrfcVeX8xs3vM7KSZLT0Nmc38q/lt9jkz++HW+7a/r6SUXP7R7PVMf0DSpyTdtOI6+yV9RdIrJB2Q9CeSfmj+vt+VdPv8378h6Rdyf08D3Cb/QtKd83/fKendG65/haTnJL1w/vYHJL0x9/eR63aR9M0Vlxd3X+l6u0j6fkk3zP/9EklPS7qspPvLuseJ1nX+e0m/Mf/37ZI+NP/3D82vf6mkxMGN4QAABOJJREFU6+efZ3/u72nC2+XHW48fv9DcLvO3l/4+Rf/T8XZ5s6R/veRjr5D0+Pzvy+f/vjz39zTV7bJw/X8k6Z4K7i9/S9IPS/r8ivf/pKSPanZu2r8p6Y/63FfcTtpSSl9MKX1pw9W++xJZafZC881LZJlmL5H14fn17pX0hvG2djK3afa9SN2+pzdK+mhK6VujblV+294u31XwfUXqcLuklL6cUnps/u+nJJ2UdHiyLZzG0seJheu0b6sPS/qJ+X3jNkn3pZSeTyl9VdJx7XbwlEcbb5eU0idbjx+f0ex8mqXrcn9Z5XWSHkwpPZdS+oakByXdOtJ2Tm3b2+VNkj44yZZllFL6D5oNR1a5TdJvp5nPaHZ+2mu0433FbbR1tOwlsq5VuS+RdXVK6WlJmv99ZMP1b9feX5pfm49o32NmpbyWSdfb5aCZPWxmn2mWjFXufUXa8v5iZjdr9j/or7QuLuH+supxYul15veFP9fsvtHlY6Pa9nt7i2YTg8ay36cSdL1dfmr+u/FhM2tOIM/9RdJ8Gf16SZ9oXVzq/WWTVbfbTveVKV57dCVz8hJZnqy7Tbb8PNdIeqVm58drvF3Sn2n2xHy3ZufJu2u3LZ3WQLfLy1NKT5nZKyR9wswekfQXS64X4r4iDX5/+R1Jd6SUmldWDnt/WdDl8aC4x5IOOn9vZvZzkm6S9GOti/f8PqWUvrLs44Ppcrv8e0kfTCk9b2Y/r9mU9tUdPzaqbb632yV9OKX0ndZlpd5fNhn0sSVrtCUnL5HlybrbxMy+bmbXpJSenj/JnlzzqX5a0kdSSt9ufe6n5/983sx+S9IvD7LRExjidpkv/yml9LiZfUrSqyT9OwW9r0jD3C5m9r2S/kDSP5mP75vPHfb+smDV48Sy65wws0skvVizJY8uHxtVp+/NzF6j2X8Cfiyl9Hxz+YrfpxKehDfeLimlZ1tv/htJ72597N9e+NhPDb6FeWzzu3C7pF9sX1Dw/WWTVbfbTveV6Mujtb1E1jHNvhdp8/e0Z3+C+RN3sx/XGyQtPdoloI23i5ld3izvmdlVkn5U0qMF31ekbrfLAc1OkP3bKaXfW3hfKfeXpY8TC9dp31ZvlPSJ+X3jmKTbbXZ06fWSbpD0Hyfa7rFtvF3M7FWS3ifp9Smlk63Ll/4+Tbbl4+pyu1zTevP1mr3ajzRb2Xjt/Pa5XNJrdfFqR2Rdfo9kZj+g2Y71n25dVvL9ZZNjkv7+/CjSvynpz+f/Id7tvpL7yItVfyT9Xc1K9HlJX5f0wPzyl0i6v3W9n5T0Zc2K/R2ty1+h2YPrcUm/J+nS3N/TALfJlZI+Lumx+d9XzC+/SdL7W9e7TtKfStq38PGfkPSIZk++/4ekF+X+nqa6XST91/Pv/U/mf7+l5PvKFrfLz0n6tqTPtv7cWNr9ZdnjhGZLva+f//vg/Gd/fH5feEXrY98x/7gvSfo7ub+XiW+Xj80ff5v7xrH55St/n0r40+F2+V8kfWH+/X9S0n/Z+th/OL8fHZf0D3J/L1PeLvO33ynpXQsfV+z9RbPhyNPzx9ETmu37+fOSfn7+fpP03vlt9ohaZ8PY5b7CKyIAAAAEEH15FAAAoApEGwAAQABEGwAAQABEGwAAQABEGwAAQABEGwAAQABEGwAAQABEGwAAQAD/P+01Z6YkXamcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "582cb88de23644cb9d7120507923bbe5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=48.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZgddZ3v8fenu7PvO9kXiEDYgjQBxpUdFQUVJSAYFAZnrjziegXRO8gMI3rnyuDIjIZdUBZZhig4yBIQREg6AgGymJAE0yQknXRn6+70+r1/nEo4tN2dTs7pPt2pz+t5ztN1qn516luVJ/U5tf2OIgIzM0uvokIXYGZmheUgMDNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQ2H5F0u8kzSnAcndImtbVy+0oSRdJeq7QdVj35CCwTiNpjaRTWozr1B1SRHwkIu7orM9vZ7kDI2JVVy/XLB8cBGZmKecgsIKSNE7SA5IqJK2W9JWsabMklUnaJmmDpB8n4/tKukvSZklbJC2UNCaZ9rSkS5LhAyU9lbTbJOmXkoZmff4aSd+UtFjSVkn3SurbTq1flLRUUpWkxyRNzpoWkg5KhkdI+k1S90JJ/5J9FCTpEEmPS6qUtFzSZ7Om3S7pRkmPSNou6UVJBybTfibp31rU9LCkryfDV0h6I5lviaRPtrEeU5J6S7LG7d5u7a2rMq6XtDHZZoslHd7WNrOewUFgBSOpCPgN8AowHjgZ+Kqk05MmNwA3RMRg4EDgvmT8HGAIMBEYAfwDUNvaIoAfAOOAQ5P2V7do81ngDGAqcCRwURu1ng18B/gUMAp4Fri7jVW7EagGDkhq3X3NQtIA4HHgV8Bo4DzgPyUdljX/ecD3gWHASuDaZPyvgHMlKfmsYcBpwD3J9DeAD5DZNt8H7pI0to0a27SHdT0N+CDwHmAocC6weW+XYd2Lg8A6238n39q3SNoC/GfWtGOBURFxTUTUJ+fYbwJmJ9MbgIMkjYyIHRHxQtb4EcBBEdEUEYsiYlvLBUfEyoh4PCLqIqIC+DHwoRbNfhIR6yKikkwozWxjPb4E/CAilkZEI/CvwMzsowIAScXAp4F/ioiaiFgCZF+zOBNYExG3RURjRPwZeAA4J6vNgxGxIFnOL7NqehYIMjt7knn+FBHrkvX9dbIuzRFxL7ACmNXG+rSnvXVtAAYBhwBK2qzfh2VYN+IgsM52dkQM3fUC/lfWtMnAuBZB8R1gTDL9YjLfPJclp1jOTMbfCTwG3CNpnaQfSerVcsGSRku6R9JbkrYBdwEjWzR7O2u4BhjYxnpMBm7IqrOSzBHH+BbtRgElwNqscdnDk4HjWqzz58gcPbRbU2R6iLyHzBEDwPlkgmLX+n5e0stZn3t4K+vbEW2ua0Q8BfyUzFHPBklzJQ3eh2VYN+IgsEJaC6zODoqIGBQRHwWIiBURcR6ZUyg/BO6XNCAiGiLi+xExA/g7Mt+yP9/K5/+AzDfoI5PTSxeQ2aHta61falFrv4h4vkW7CqARmJA1bmKLz3mmxecMjIh/7GAddwPnJN/OjyNzNEHy/ibgMmBEErqv0fr6Vid/+2eNyw6idtc1In4SEccAh5EJ6m91sHbrphwEVkgLgG2Svi2pn6RiSYdLOhZA0gWSRkVEM7AlmadJ0omSjkhOw2wjc7qiqZXPHwTsALZIGk9uO6yfAVfuOpcvaYikz7RsFBFNwIPA1ZL6SzqEd4fUb4H3SLpQUq/kdaykQztSRES8RCZsbgYei4hd22UAmdCrSOr7ApkjgtY+owJ4C7gg2eZfJHMNZo/rmtR6XHIEVg3spPVtbz2Ig8AKJtlpfpzMOfDVwCYyO7ghSZMzgNcl7SBz4Xh2ROwk8+31fjIhsBR4hsxpn5a+D7wX2Ao8QmYHva+1PkTmqOSe5DTTa8BH2mh+WbIOb5M5jXU3UJd8znYyF1xnA+uSNj8E+uxFOXcDp5C5eLyrviXA/wP+BGwAjgD+2M5n/D2ZYNxM5pv97iObPazrYDJHHlXAm8n877qTyXoe+YdpzDqXpB8CB0RElz/xbNYRPiIwy7PkOYEjk3vuZ5G56P1Qoesya0vJnpuY2V4aROb0zThgI5lTNg8XtCKzdvjUkJlZyvnUkJlZyvXIU0MjR46MKVOmFLoMM7MeZdGiRZsiYlTL8T0yCKZMmUJZWVmhyzAz61EkvdnaeJ8aMjNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlHARmZj1AbX0TV897nbWVNXn/bAeBmVkPcP+itdz+/BrWb92Z9892EJiZdXNNzcHNz63mqIlDOXbKsLx/voPAzKyb+/3rb/Pm5hq+9MFpSPv6s9ttcxCYmXVjEcHP/7CKySP6c/phB3TKMhwEZmbd2MI1Vby8dguXvH8qxUX5PxoAB4GZWbc29w9vMKx/L845ZmKnLcNBYGbWTa3cuIMnlm7kwhOm0K93cactJy9BIOkMScslrZR0RSvTL5JUIenl5HVJ1rQ5klYkrzn5qMfMbH9w87Or6FNSxJwTJnfqcnL+YRpJxcCNwKlAObBQ0ryIWNKi6b0RcVmLeYcD/wSUAgEsSuatyrUuM7OebOP2nTz457f4TOkERgzs06nLyscRwSxgZUSsioh64B7grA7OezrweERUJjv/x4Ez8lCTmVmPdsfza2hobuaSD0zr9GXlIwjGA2uz3pcn41r6tKTFku6XtOuqR0fnRdKlksoklVVUVOShbDOz7qm6rpG7Xvgrp80Yw9SRAzp9efkIgtbuZ4oW738DTImII4EngDv2Yt7MyIi5EVEaEaWjRv3Nby+bme037itby9baBi794IFdsrx8BEE5kH1f0wRgXXaDiNgcEXXJ25uAYzo6r5lZmjQ2NXPLc6spnTyMYybnvzuJ1uQjCBYC0yVNldQbmA3My24gaWzW208AS5Phx4DTJA2TNAw4LRlnZpZKj772NuVVtVz6wc6/NrBLzncNRUSjpMvI7MCLgVsj4nVJ1wBlETEP+IqkTwCNQCVwUTJvpaR/JhMmANdERGWuNZmZ9UTNzcF/zl/JtFEDOOXQMV223JyDACAiHgUebTHu/2QNXwlc2ca8twK35qMOM7Oe7DeL17Hs7e3cMHsmRZ3UnURr/GSxmVk30NDUzPWP/4VDDhjEx48c16XLdhCYmXUD9y8qZ83mGr51+sFdejQADgIzs4Lb2dDEDU+s4L2ThnLSIaO7fPkOAjOzArvrhTd5e9tOvnX6IZ3ywzN74iAwMyug7TsbuHH+Sj4wfSQnHDiiIDU4CMzMCujW59ZQVdPAN087uGA1OAjMzAqkqrqem55dxemHjeGoiUMLVoeDwMysQH72zBtU1zfyjQIeDYCDwMysIDZs28ntz6/hkzPH854xgwpai4PAzKwA/uOpFTQ1B1895T2FLsVBYGbW1RaXb+GeBWuZPWsik0b0L3Q5DgIzs660trKGL96+kAOG9O0WRwOQp07nzMxsz6qq65lz2wIamoJ7Lp3FyE7+LeKOchCYmXWBnQ1N/P0vyiivquWui4/joNEDC13Sbj41ZGbWyZqbg6/f9zJlb1bx488exaypwwtd0rs4CMzMOtm1jy7l0Vff5rsfO5Qzu7iL6Y5wEJiZdaJbn1vNLc+t5qK/m8LF759a6HJa5SAwM+sk9y8q558fWcLph43he2fOKEjPoh3hi8VmZnlWWV3P9x5+jUcWr2fW1OHcMPtoirv4x2b2Rl6OCCSdIWm5pJWSrmhl+tclLZG0WNKTkiZnTWuS9HLympePeszMCuWx19/mtOuf4fevv823Tj+YX11yHH17FRe6rHblfEQgqRi4ETgVKAcWSpoXEUuymr0ElEZEjaR/BH4EnJtMq42ImbnWYWZWSFtrGrj6N6/z0EtvMWPsYO68+DgOHTu40GV1SD5ODc0CVkbEKgBJ9wBnAbuDICLmZ7V/AbggD8s1M+sW5i/byLcfWExldT2Xnzydy046iF7FPecSbD4qHQ+szXpfnoxry8XA77Le95VUJukFSWe3NZOkS5N2ZRUVFblVbGaWJ/OXb+QLty9kWP/e/PeX38fXTn1PjwoByM8RQWtXQKLVhtIFQCnwoazRkyJinaRpwFOSXo2IN/7mAyPmAnMBSktLW/18M7OutHH7Tr553yscPGYQD1/2vm5/LaAt+YitcmBi1vsJwLqWjSSdAlwFfCIi6naNj4h1yd9VwNPA0XmoycysUzU3B9+47xV21DXyH+cf3WNDAPITBAuB6ZKmSuoNzAbedfePpKOBn5MJgY1Z44dJ6pMMjwTeR9a1BTOz7uqW51bz7IpNfPfMGQX/YZlc5XxqKCIaJV0GPAYUA7dGxOuSrgHKImIe8H+BgcCvkwcq/hoRnwAOBX4uqZlMKF3X4m4jM7Nu59XyrfzosWWcNmMMFxw3qdDl5EwRPe90e2lpaZSVlRW6DDNLoeq6Rs78j+eorW/id5d/gGEDehe6pA6TtCgiSluO95PFZmZ74ep5r7NmczW/uuT4HhUC7elZ9ziZmRXQb15Zx68XlfPlDx/ECQeOKHQ5eeMgMDPrgLWVNXznwVc5etJQLj9leqHLySsHgZnZHqzbUsucWxcA8JPZR/e4B8b2xNcIzMzasapiBxfesoBttQ3c+oVjmTi8f6FLyjsHgZlZG5au38aFt7xIc8Ddlx7P4eOHFLqkTuEgMDNrxaI3q/jCbQvo37uEuy7pXj82n28OAjOzFp5bsYlL7yxj9KA+3Hnxcfvl6aBsDgIzsyz/89rbfOXul5g6cgB3XjKL0YP6FrqkTucgMLPUq6yuZ97Lb/HAn9/i1be2MnPiUG7/wrEM7b9/PDC2Jw4CM0ul+sZm5i/fyAOLypm/fCMNTcGMsYP53pkzOG/WRPr3Ts/uMT1ramapVVPfyKqKalZu3LH7tWBNJZXV9Ywc2Ic5J0zh08dM6DE/LZlvDgIz2+9s3L6TJ5du5KllG1mybhtvbandPa24SEwe3p/3HzSSs48exwenj6JkP3tAbG85CMxsv7CqYge/X7KB37/+Ni+t3UIETBjWj2MmD+PcYydy0OiBHDR6IFNGDKB3Sbp3/C05CMysR9q4fScLV1excE0lz63cxMqNOwA4fPxgvnbKezh1xhgOOWAQyW+gWDscBGbWI/x1cw0vrNrMgjWVlK2pZM3mGgD69SqmdMowLjx+MqfMGMP4of0KXGnP4yAws26prrGJBasrmb+sgqeXb2TVpmoAhvXvRemU4XzuuMkcO3U4h40bvN91AtfVHARmVnDNzUHFjjre2lLLsvXbmb98I39cuYma+iZ6lxRxwrQRfP6Eybx/+kgOHDXQp3vyzEFgZp1uR10j67fUsm7rTtZtqWX9llre2pIZfmtLLeu31tLQ9M7P5o4f2o9Pv3cCJx4yihOmjaRf7+ICVr//y1sQSDoDuIHMD9jfHBHXtZjeB/gFcAywGTg3ItYk064ELgaagK9ExGP5qsvMuk5EUF5Vy4urK1mwejOLy7eybkst23Y2vqtdkWDM4L6MG9qPmROH8tEjxjJ+aF/GD+vH5BEDmDZygL/1d6G8BIGkYuBG4FSgHFgoaV5ELMlqdjFQFREHSZoN/BA4V9IMYDZwGDAOeELSeyKiKR+1mdm+iwh2NjSzfWcDO+oaaWoOmiJoag6am6GxuZmm5mD5hu0sWF3JgtWVrN+6E4Ch/Xtx9MShzJo6nHFD+zF2SGbHP25oP0YP6uPz+t1Ivo4IZgErI2IVgKR7gLOA7CA4C7g6Gb4f+KkykX8WcE9E1AGrJa1MPu9PearNCqS5OdhS28DmHXVs2lFPZXU9m6vr2FrTQENTM/VNQWNTMw1NzTQ0Z4Yj9vy5lrvmyOzEG5uC+l3/Bk3NNDQFNfWNbN+ZeW2rbaCxuWP/KKMH9WHW1OEcN3U4s6aOYProgRQV+Vt9T5CvIBgPrM16Xw4c11abiGiUtBUYkYx/ocW841suQNKlwKUAkyZNylPZli8RwZuba3hx9WZeXFXJgjWVrNtSS3v7kN7FRZQUi17FRfQqFiVFRXi/0TUkZW37zPbvVVxESZEYNbAPB44ayKC+JQzq22v334F9iikpKqK4SBRJFBeJkiJRlDypO3lEf5/O6aHyFQSt/eu33AW01aYj8xIRc4G5AKWlpf7e2EWqqutZvmE7W2sbiMjs8ANojqA5YGtNPQvWVLFg9WY2bKsDYMSA3hw7ZThnzxzPiIG9GTGwDyMH9Gb4wN6MGNCHof17UVIk7zTMuol8BUE5MDHr/QRgXRttyiWVAEOAyg7Oa52ssamZlRU7WLZ+O0vf3say9dtZ9va23Tv39owZ3Ifjpo5g1tThHD9tuG/vM+th8hUEC4HpkqYCb5G5+Ht+izbzgDlkzv2fAzwVESFpHvArST8mc7F4OrAgT3VZGyq21/HSX6t4ae0WXvprFYvLt1JTn7k+37u4iINGD+R9B43k0AMGc/ABgxgxsDdCFBWR+SuQoF/vEsYN6esdv1kPlpcgSM75XwY8Rub20Vsj4nVJ1wBlETEPuAW4M7kYXEkmLEja3UfmwnIj8GXfMZR/ldX1PLuigmeWV7DwzUrWVmZ6YywpEjPGDeYzx0xg5qShzBg7hGmjBviODrMUUfTA2zRKS0ujrKys0GV0a03NweLyLTy9vIKn/1LB4vJMb4zD+vfi+GkjeO+kYRw9aSiHjx9C315+WMcsDSQtiojSluP9ZPF+pKq6nj+sqGD+so0885cKqmoakGDmxKFcfvJ0PnzwaI4YP4Ri35pjZlkcBD1YRLBk/TbmL9vI/OUVvPTXKpojc9fOiQeP5kMHj+KD00cxbEA6fnfVzPaNg6AH2lrbwK/L1nLnC2/yZtIV75EThnDZSdM56ZDRHDl+iB/kMbMOcxD0IH/ZsJ07nl/Dg39+i9qGJkonD+PLHz6IDx8yitGD+ha6PDProRwE3VxE8PslG7j9j2v406rN9C4p4qyjxjHn76Zw+PghhS7PzPYDDoJubPWmaq58cDEvrKpk3JC+/O8zDmb2sZMY7nP+ZpZHDoJuqKGpmZueXcW/P7GCPiVF/Osnj+CzpRMo8b39ZtYJHATdzCtrt3DFg6+ydP02PnL4AXz/E4cxerDP/5tZ53EQdBM19Y38+Pd/4dY/rmbkwD787IJjOOPwAwpdlpmlgIOgG6iqruf8m19k6fptnH/cJL59xiEM6der0GWZWUo4CApsa00DF9zyIm9U7OC2i47lxENGF7okM0sZX30soK21DVx464us2LCDn194jEPAzArCQVAg23c2MOfWBSxdv43/uuC9nHiwQ8DMCsNBUAA76hq56LaFvPbWVm48/72cfOiYQpdkZinmawRdrKa+kS/etpCX127hp+cdzWmH+c4gMyssHxF0oYrtdXzx9oWUvVnJDbNn8pEjxha6JDMzHxF0hfrGZm5/fjU/eXIldY1NXH/uTM48clyhyzIzAxwEnSoieHLpRq59dCmrN1Vz0iGjuepjh3LgqIGFLs3MbDcHQSdZsWE71/x2Cc+u2MSBowZw+xeO5cO+M8jMuqGcgkDScOBeYAqwBvhsRFS1aDMT+C9gMNAEXBsR9ybTbgc+BGxNml8UES/nUlN3cNMfVnHd/yxjQO9i/unjM7jg+Mn+MXgz67ZyPSK4AngyIq6TdEXy/tst2tQAn4+IFZLGAYskPRYRW5Lp34qI+3Oso9t4YdVm/vV3Szn10DFc9+kj3WW0mXV7uX5NPQu4Ixm+Azi7ZYOI+EtErEiG1wEbgVE5Lrdb2lrTwNfufZkpIwZw/bkzHQJm1iPkGgRjImI9QPK33ZPgkmYBvYE3skZfK2mxpOsl9Wln3ksllUkqq6ioyLHs/IsIvvPQq1Rsr+OG2TMZ0MeXX8ysZ9hjEEh6QtJrrbzO2psFSRoL3Al8ISKak9FXAocAxwLD+dvTSrtFxNyIKI2I0lGjut8Bxa8XlfPIq+v5xmkHc+SEoYUux8ysw/b4tTUiTmlrmqQNksZGxPpkR7+xjXaDgUeA70bEC1mfvT4ZrJN0G/DNvaq+m1i9qZqr573OCdNG8KUPTit0OWZmeyXXU0PzgDnJ8Bzg4ZYNJPUGHgJ+ERG/bjFtbPJXZK4vvJZjPV2uvrGZy+95id4lRfz43KMoKlKhSzIz2yu5BsF1wKmSVgCnJu+RVCrp5qTNZ4EPAhdJejl5zUym/VLSq8CrwEjgX3Ksp8td/8RfWFy+les+dSRjh/QrdDlmZnstpyuaEbEZOLmV8WXAJcnwXcBdbcx/Ui7LL7Tn39jEz555g/NmTfTPSppZj+WnnPbR5h11fP3eV5g6cgDfO3NGocsxM9tnDoJ9UNfYxJfuXERVTT0/mX00/Xv7VlEz67m8B9tLEcFVD71G2ZtV/PT8ozl8/JBCl2RmlhMfEeylm55dxf2LyvnqKdPdlbSZ7RccBHvhiSUb+MHvlnHmkWO5/OTphS7HzCwvHAQdtOztbVx+z0scMX4I//aZo8g8+mBm1vM5CDpg0446Lr69jIF9S7jp86X07VVc6JLMzPLGF4v3oK6xiX+4cxGbq+u470snMGZw30KXZGaWVw6CdtQ3NvO1e1+m7M0qbjz/ve5Mzsz2Sw6CNtTWN/GPv1zE08sr+O7HDuVjR44tdElmZp3CQdCKbTsbuOT2Mha+WckPPnUE582aVOiSzMw6jYOghc076phz2wKWrd/OT2YfzceP8rMCZrZ/cxBkWb+1lgtufpHyqlpu+nwpJx7S7g+umZntFxwEiTWbqvnczS+ytbaBX3xxFsdNG1HokszMuoSDANjZ0MS5c/9EfWMzd//98Rwxwf0HmVl6+IEyYG1lDRu21XHVx2Y4BMwsdRwEQGV1PQBjBvcpcCVmZl3PQQBU1TQAMKx/7wJXYmbW9RwEQFVN5ohg+AAHgZmlT85BIGm4pMclrUj+DmujXVPWj9fPyxo/VdKLyfz3SuryvfGuU0M+IjCzNMrHEcEVwJMRMR14MnnfmtqImJm8PpE1/ofA9cn8VcDFeahpr1RV19OvVzH9ertXUTNLn3wEwVnAHcnwHcDZHZ1RmU79TwLu35f586Wypt6nhcwstfIRBGMiYj1A8retx3H7SiqT9IKkXTv7EcCWiGhM3pcD41ubWdKlyfxlFRUVeSj7HVXV9Qwb0Cuvn2lm1lN06IEySU8AB7Qy6aq9WNakiFgnaRrwlKRXgW2ttIvWZo6IucBcgNLS0lbb7KvKmgZfHzCz1OpQEETEKW1Nk7RB0tiIWC9pLLCxjc9Yl/xdJelp4GjgAWCopJLkqGACsG4v1yFnVdX1TBnRv6sXa2bWLeTj1NA8YE4yPAd4uGUDScMk9UmGRwLvA5ZERADzgXPam7+zVVXX+4jAzFIrH0FwHXCqpBXAqcl7JJVKujlpcyhQJukVMjv+6yJiSTLt28DXJa0kc83gljzU1GH1jc1sr2v0xWIzS62cO52LiM3Aya2MLwMuSYafB45oY/5VwKxc69hXW5KHyYY5CMwspVL/ZHHlrqeKfWrIzFLKQbDrqWLfPmpmKZX6IKiqznQ452sEZpZWqQ8Cnxoys7RLfRBUJaeGhjoIzCylUh8EldX1DOpTQu+S1G8KM0up1O/9qmrqfeuomaVa6oOgstpBYGbplvogqKqpZ3h/3zpqZunlIKhu8BGBmaWag6DGHc6ZWbqlOgh2NjRRU9/kh8nMLNVSHQRVNf7RejOzVAfBrn6GhrufITNLsVQHwa5+hnxEYGZpluog2N3PkK8RmFmKpToIqqr9ozRmZqkOgl3XCIb28zUCM0uvVAdBVU09Q/r1oqQ41ZvBzFIupz2gpOGSHpe0Ivk7rJU2J0p6Oeu1U9LZybTbJa3OmjYzl3r2VmV1va8PmFnq5fpV+ArgyYiYDjyZvH+XiJgfETMjYiZwElAD/D6rybd2TY+Il3OsZ69knir2aSEzS7dcg+As4I5k+A7g7D20Pwf4XUTU5LjcvKisbvARgZmlXq5BMCYi1gMkf0fvof1s4O4W466VtFjS9ZL65FjPXqmqdj9DZmYle2og6QnggFYmXbU3C5I0FjgCeCxr9JXA20BvYC7wbeCaNua/FLgUYNKkSXuz6FZFBJU1vkZgZrbHIIiIU9qaJmmDpLERsT7Z0W9s56M+CzwUEQ1Zn70+GayTdBvwzXbqmEsmLCgtLY091b0nNfVN1Dc2+xkCM0u9XE8NzQPmJMNzgIfbaXseLU4LJeGBJJG5vvBajvV02O5+hnxqyMxSLtcguA44VdIK4NTkPZJKJd28q5GkKcBE4JkW8/9S0qvAq8BI4F9yrKfDdvc86iMCM0u5PZ4aak9EbAZObmV8GXBJ1vs1wPhW2p2Uy/Jz4Z5HzcwyUvtIrX+LwMwsI71BkHRB7buGzCzt0hsENfUUCQb39akhM0u31AZBZfIwWVGRCl2KmVlBpTYIqmrqfceQmRkpDoLK6no/Q2BmRoqDoKq6gaHuedTMLL1B4H6GzMwyUhkEEZHpedRBYGaWziDYXtdIY3P4GoGZGSkNgqpq9zNkZrZLKoPA/QyZmb0jlUHgfobMzN6RyiCodD9DZma7pTIIfI3AzOwdqQyCypp6SorEoD45/RyDmdl+IZVBsOsZgswvZJqZpVsqg8D9DJmZvSOVQZDpedS3jpqZQUqDoLLa/QyZme2ScxBI+oyk1yU1Syptp90ZkpZLWinpiqzxUyW9KGmFpHsldfoeuqqmwc8QmJkl8nFE8BrwKeAPbTWQVAzcCHwEmAGcJ2lGMvmHwPURMR2oAi7OQ01tam4OtrjnUTOz3XIOgohYGhHL99BsFrAyIlZFRD1wD3CWMrftnATcn7S7Azg715ras21nA83hp4rNzHbpqmsE44G1We/Lk3EjgC0R0dhi/N+QdKmkMkllFRUV+1zIO/0MOQjMzAA69ESVpCeAA1qZdFVEPNyRj2hlXLQz/m9HRswF5gKUlpa22qYjdvcz5CAwMwM6GAQRcUqOyykHJma9nwCsAzYBQyWVJEcFu8Z3mt39DPnUkJkZ0HWnhhYC05M7hHoDs4F5ERHAfOCcpN0coCNHGPvsnX6G/ByBmRnk5/bRT0oqB04AHpH0WDJ+nKRHAZJv+0K7+FoAAAVySURBVJcBjwFLgfsi4vXkI74NfF3SSjLXDG7Jtab2VNb4GoGZWbace12LiIeAh1oZvw74aNb7R4FHW2m3isxdRV2iqrqePiVF9OtV3FWLNDPr1lL3ZPGup4rd4ZyZWUbqgqCqpt7PEJiZZUldEFRWu8M5M7NsqQsC9zNkZvZuqQsC9zxqZvZuqQqCxqZmttb6iMDMLFuqgmBLbfJUsY8IzMx2S1UQvPNUsYPAzGyXVAXB7p5HfWrIzGy3VAXBOz2P+vZRM7NdUhUEu3se9akhM7PdUhUEu48IfGrIzGy3dAVBdT39exfT1x3OmZntlqogmD5mIGceObbQZZiZdSs5d0Pdk5x77CTOPXZSocswM+tWUnVEYGZmf8tBYGaWcg4CM7OUcxCYmaVcTkEg6TOSXpfULKm0jTYTJc2XtDRpe3nWtKslvSXp5eT10dY+w8zMOk+udw29BnwK+Hk7bRqBb0TEnyUNAhZJejwiliTTr4+If8uxDjMz20c5BUFELAXa/SH4iFgPrE+Gt0taCowHlrQ5k5mZdZkuvUYgaQpwNPBi1ujLJC2WdKukYe3Me6mkMkllFRUVnVypmVl6KCLabyA9ARzQyqSrIuLhpM3TwDcjoqydzxkIPANcGxEPJuPGAJuAAP4ZGBsRX9xj0VIF8Oae2rVhZLLMNPM28DZI+/pDOrfB5IgY1XLkHk8NRcQpuS5ZUi/gAeCXu0Ig+ewNWW1uAn7bkc9rbUX2opayiGj1wnZaeBt4G6R9/cHbIFunnxpS5gLCLcDSiPhxi2nZHf98kszFZzMz60K53j76SUnlwAnAI5IeS8aPk/Ro0ux9wIXASa3cJvojSa9KWgycCHwtl3rMzGzv5XrX0EPAQ62MXwd8NBl+Dmj1tqKIuDCX5e+juQVYZnfjbeBtkPb1B2+D3fZ4sdjMzPZv7mLCzCzlHARmZimXqiCQdIak5ZJWSrqi0PV0heRBvY2SXssaN1zS45JWJH/bfJCvp2urr6uUbYO+khZIeiXZBt9Pxk+V9GKyDe6VtF//mLekYkkvSfpt8j5V69+e1ASBpGLgRuAjwAzgPEkzCltVl7gdOKPFuCuAJyNiOvBk8n5/tauvq0OB44EvJ//uadoGdcBJEXEUMBM4Q9LxwA/J9PU1HagCLi5gjV3hcmBp1vu0rX+bUhMEwCxgZUSsioh64B7grALX1Oki4g9AZYvRZwF3JMN3AGd3aVFdKCLWR8Sfk+HtZHYE40nXNoiI2JG87ZW8AjgJuD8Zv19vA0kTgI8BNyfvRYrWf0/SFATjgbVZ78uTcWk0JukMcFengKMLXE+XaNHXVaq2QXJa5GVgI/A48AawJSIakyb7+/+Hfwf+N9CcvB9Buta/XWkKgtaeZfC9symR9HX1APDViNhW6Hq6WkQ0RcRMYAKZo+NDW2vWtVV1DUlnAhsjYlH26Faa7pfr3xG5/h5BT1IOTMx6PwFYV6BaCm2DpLERsT7p5mNjoQvqTG30dZWqbbBLRGxJOok8HhgqqST5Vrw//394H/CJpEeDvsBgMkcIaVn/PUrTEcFCYHpyp0BvYDYwr8A1Fco8YE4yPAd4uIC1dKp2+rpK0zYYJWloMtwPOIXMtZL5wDlJs/12G0TElRExISKmkPl//1REfI6UrH9HpOrJ4uQbwb8DxcCtEXFtgUvqdJLuBj5MpsvdDcA/Af8N3AdMAv4KfCYiWl5Q3i9Iej/wLPAq75wf/g6Z6wRp2QZHkrkYWkzmy999EXGNpGlkbpoYDrwEXBARdYWrtPNJ+jCZLvPPTOP6tyVVQWBmZn8rTaeGzMysFQ4CM7OUcxCYmaWcg8DMLOUcBGZmKecgMDNLOQeBmVnK/X+btAsoLaHQbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch:  0 batch:  0\n",
      "Calc Hessian\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c1a184df0cad4db69a1cae18da147d24",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=441.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-c60324ee5190>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mprogress\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     24\u001b[0m     \u001b[0mindex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandperm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 25\u001b[0;31m     \u001b[0mplot_prediction_map\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     26\u001b[0m     \u001b[0mX\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mX_train\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m# take random sample of X and y (bc it is faster)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     27\u001b[0m     \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0my_train\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-10e6dbf0283e>\u001b[0m in \u001b[0;36mplot_prediction_map\u001b[0;34m(epoch)\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mplot_prediction_map\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mepoch\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 51\u001b[0;31m     \u001b[0mfidelity\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_fidelity\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     52\u001b[0m     \u001b[0mZ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfidelity\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrid_resolution\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrid_resolution\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m     \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-10e6dbf0283e>\u001b[0m in \u001b[0;36mget_fidelity\u001b[0;34m()\u001b[0m\n\u001b[1;32m     39\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget_fidelity\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m     \u001b[0mfidl\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mprediction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mX_map\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmeasurements\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeasures\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m     \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfidl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m     \u001b[0mfidelity\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-10e6dbf0283e>\u001b[0m in \u001b[0;36mprediction\u001b[0;34m(params, x, y, measurements, device)\u001b[0m\n\u001b[1;32m     28\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m         \u001b[0mfidelities\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 30\u001b[0;31m         \u001b[0mout1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcirc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     31\u001b[0m         \u001b[0mout1_copy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mout1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclone\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdetach\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad_\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m         \u001b[0mout2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmeasurements\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Dropbox/34_5_vqc_loss_landscapes/vqc_loss_landscapes/torchcirq.py\u001b[0m in \u001b[0;36mbuild_circuit\u001b[0;34m(self, params, x, device)\u001b[0m\n\u001b[1;32m    286\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mbuild_circuit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'cpu'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 288\u001b[0;31m         \u001b[0mcircuit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCZ_layer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRot_layer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    289\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayers\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    290\u001b[0m             \u001b[0mcircuit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRot_layer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcircuit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Dropbox/34_5_vqc_loss_landscapes/vqc_loss_landscapes/torchcirq.py\u001b[0m in \u001b[0;36mRot_layer\u001b[0;34m(self, params, layer, x, device)\u001b[0m\n\u001b[1;32m    271\u001b[0m             \u001b[0mangles\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpar\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mpar\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    272\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mRot_Matrix\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 273\u001b[0;31m                 \u001b[0mRot_Matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkronecker_prod\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mR3\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mangles\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mRot_Matrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    274\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    275\u001b[0m                 \u001b[0mRot_Matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mR3\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mangles\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Dropbox/34_5_vqc_loss_landscapes/vqc_loss_landscapes/complex.py\u001b[0m in \u001b[0;36mkronecker_prod\u001b[0;34m(x, y)\u001b[0m\n\u001b[1;32m    269\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Inputs must be complex matrices!\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    270\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 271\u001b[0;31m     return einsum(\"ab,cd->acbd\", x, y).reshape(\n\u001b[0m\u001b[1;32m    272\u001b[0m         \u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    273\u001b[0m     )\n",
      "\u001b[0;32m~/Dropbox/34_5_vqc_loss_landscapes/vqc_loss_landscapes/complex.py\u001b[0m in \u001b[0;36meinsum\u001b[0;34m(equation, a, b, real_part, imag_part)\u001b[0m\n\u001b[1;32m    180\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mimag_part\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    181\u001b[0m         i = torch.einsum(equation, real(a), imag(b)).add_(\n\u001b[0;32m--> 182\u001b[0;31m             \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mequation\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    183\u001b[0m         )\n\u001b[1;32m    184\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/QC/lib/python3.7/site-packages/torch/functional.py\u001b[0m in \u001b[0;36meinsum\u001b[0;34m(equation, *operands)\u001b[0m\n\u001b[1;32m    239\u001b[0m         \u001b[0;31m# the old interface of passing the operands as one list argument\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    240\u001b[0m         \u001b[0moperands\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moperands\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 241\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_C\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_VariableFunctions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mequation\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moperands\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    242\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    243\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# =================================================================================\n",
    "# Main Training\n",
    "# =================================================================================\n",
    "print(\"start training width{} layers{} epochs{}\".format(width, layers, epochs))\n",
    "init = torch.zeros(2**width).to(device)\n",
    "init[0] = 1.0 # Define initial state\n",
    "\n",
    "init = make_complex(init)\n",
    "lr = 0.1\n",
    "params = torch.randn((layers, width ,2 ,3), requires_grad=True, device=device)\n",
    "optimizer = torch.optim.SGD([params], lr=lr)\n",
    "model = Model_Z_measure(params, width = width, layers = layers, device=device)\n",
    "# circ = model.build_circuit(params, x=torch.randn((3)))\n",
    "\n",
    "x = torch.randn((3))\n",
    "\n",
    "circ = model.build_circuit\n",
    "print(\"circuit built\")\n",
    "progress = tqdm(range(epochs), bar_format='{desc}')\n",
    "\n",
    "loss_list = []\n",
    "\n",
    "for i in progress:\n",
    "    index = torch.randperm(len(X_train))\n",
    "    plot_prediction_map(i)\n",
    "    X = X_train[index][0:10] # take random sample of X and y (bc it is faster)\n",
    "    y = y_train[index][0:10]\n",
    "    loss = test_loss(circ, params, X, y, init, measurements=model.measures, device=device)\n",
    "    ev, H = find_heigenvalues(loss, params)\n",
    "    plt.plot(ev)\n",
    "    plt.title(\"Hessian eigenvalues\")\n",
    "    plt.show()\n",
    "    batch_idx = 0\n",
    "    for Xbatch, ybatch in iterate_minibatches(X_train[index], y_train[index], batch_size=batch_size):\n",
    "        if batch_idx%10 == 0:\n",
    "            print(\"epoch: \", i, \"batch: \", batch_idx)\n",
    "        batch_idx += 1\n",
    "        loss = test_loss(circ, params, Xbatch, ybatch, init, measurements=model.measures, device=device)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward(retain_graph=True)\n",
    "        optimizer.step()\n",
    "    print(\"Calc Hessian\")\n",
    "\n",
    "    loss_list.append(loss.item())\n",
    "    optimizer = torch.optim.SGD([params], lr=lr)\n",
    "    progress.set_description_str('loss: {:3f}'.format(loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}