Hessian_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Authors: Zhewei Yao <https://github.com/yaozhewei>, Amir Gholami <http://amirgholami.org/>\n",
    "\n",
    "\n",
    "This tutorial shows how to compute the Hessian information using (randomized) numerical linear algebra for both explicit Hessian (the matrix is given) as well as implicit Hessian (the matrix is ungiven).\n",
    "\n",
    "We'll start by doing the necessary imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch \n",
    "from torchvision import datasets, transforms\n",
    "from utils import * # get the dataset\n",
    "from pyhessian import hessian # Hessian computation\n",
    "from density_plot import get_esd_plot # ESD plot\n",
    "from pytorchcv.model_provider import get_model as ptcv_get_model # model\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enable cuda devices\n",
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\"\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"0\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Power Iteration with Numpy\n",
    "\n",
    "The following part shows how to use power iteration to get the top eigenvalue of a matrix without explicitly having access to it in numpy. We start by creating a random matrix B, compute its ground truth eigenvalues using numpy, and then compare the results with matrix-free power iteration (which does not need direct access to the matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10 # the matrix size\n",
    "\n",
    "# generate a random matrix\n",
    "A = np.random.randn(n, n)\n",
    "B = A @ A.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use numpy to compute the ground truth eigenvalues. We will then check the results with this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The top eigenvalue of B is 26.0349\n"
     ]
    }
   ],
   "source": [
    "# use np.eigs to get the top eigenvalue of B\n",
    "eigs, _ = np.linalg.eig(B)\n",
    "\n",
    "print(\"The top eigenvalue of B is %.4f\"%np.sort(eigs)[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to comptue the top eigenvalue of B without explicitly accessing B. To do so, we will use a method called Power Iteration:\n",
    "\n",
    "https://en.wikipedia.org/wiki/Power_iteration\n",
    "\n",
    "The algorithm is very simple and efficiet to compute the top eigenvalue:\n",
    "$$v_{i+1} = \\frac{Bv_i}{\\|Bv_i\\|}.$$\n",
    "\n",
    "As such, we only need to have access to the *application* of B to a given vector $v_i$ and not the matrix B itself. This application is commonly referred to as *matvec* in literature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step  1 current estimated top eigvalue: 12.6346\n",
      "Step  2 current estimated top eigvalue: 23.0936\n",
      "Step  3 current estimated top eigvalue: 24.2308\n",
      "Step  4 current estimated top eigvalue: 24.7656\n",
      "Step  5 current estimated top eigvalue: 25.0259\n",
      "Step  6 current estimated top eigvalue: 25.1619\n",
      "Step  7 current estimated top eigvalue: 25.2389\n",
      "Step  8 current estimated top eigvalue: 25.2858\n",
      "Step  9 current estimated top eigvalue: 25.3164\n",
      "Step 10 current estimated top eigvalue: 25.3378\n",
      "Step 11 current estimated top eigvalue: 25.3538\n",
      "Step 12 current estimated top eigvalue: 25.3667\n",
      "Step 13 current estimated top eigvalue: 25.3779\n",
      "Step 14 current estimated top eigvalue: 25.3881\n",
      "Step 15 current estimated top eigvalue: 25.3979\n",
      "Step 16 current estimated top eigvalue: 25.4076\n",
      "Step 17 current estimated top eigvalue: 25.4173\n",
      "Step 18 current estimated top eigvalue: 25.4271\n",
      "Step 19 current estimated top eigvalue: 25.4372\n",
      "Step 20 current estimated top eigvalue: 25.4475\n",
      "Step 21 current estimated top eigvalue: 25.4580\n",
      "Step 22 current estimated top eigvalue: 25.4688\n",
      "Step 23 current estimated top eigvalue: 25.4798\n",
      "Step 24 current estimated top eigvalue: 25.4911\n",
      "Step 25 current estimated top eigvalue: 25.5027\n",
      "Step 26 current estimated top eigvalue: 25.5145\n",
      "Step 27 current estimated top eigvalue: 25.5264\n",
      "Step 28 current estimated top eigvalue: 25.5386\n",
      "Step 29 current estimated top eigvalue: 25.5509\n",
      "Step 30 current estimated top eigvalue: 25.5634\n",
      "Step 31 current estimated top eigvalue: 25.5760\n",
      "Step 32 current estimated top eigvalue: 25.5888\n",
      "Step 33 current estimated top eigvalue: 25.6016\n",
      "Step 34 current estimated top eigvalue: 25.6144\n",
      "Step 35 current estimated top eigvalue: 25.6273\n",
      "Step 36 current estimated top eigvalue: 25.6402\n",
      "Step 37 current estimated top eigvalue: 25.6531\n",
      "Step 38 current estimated top eigvalue: 25.6659\n",
      "Step 39 current estimated top eigvalue: 25.6786\n",
      "Step 40 current estimated top eigvalue: 25.6912\n",
      "Finished Power Iteration\n",
      "\n",
      "Ground Truth Top Eigenvalue: 26.0349\n",
      "Result with matrix-free Power Iteration: 25.6912\n"
     ]
    }
   ],
   "source": [
    "# use power iteration to get the top eigenvalue of B\n",
    "v = np.random.randn(n, 1)\n",
    "for i in range(40):\n",
    "    v = v / np.linalg.norm(v)\n",
    "    eig_power_iteration = v.T @ B @ v\n",
    "    print(\"Step %2d current estimated top eigvalue: %.4f\"%(i+1,eig_power_iteration))\n",
    "    v = B @ v\n",
    "print(\"Finished Power Iteration\\n\")\n",
    "print(\"Ground Truth Top Eigenvalue: %.4f\"%np.sort(eigs)[-1])\n",
    "print(\"Result with matrix-free Power Iteration: %.4f\"%eig_power_iteration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see the result of the power iteration and the one we got from numpy match very well.\n",
    "\n",
    "We can apply the same techinique for neural networks as well, and in particular use it to compute eigenvalues of Hessian!\n",
    "\n",
    "Importantly there has been a lot of misconception that we can not use Hessian for real world applications since\n",
    "we need to explicitly form it. Next you will see that this is not correct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Power Iteration for NN Hessian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Files already downloaded and verified\n"
     ]
    }
   ],
   "source": [
    "# get the model \n",
    "model = ptcv_get_model(\"resnet20_cifar10\", pretrained=True)\n",
    "# change the model to eval mode to disable running stats upate\n",
    "model.eval()\n",
    "\n",
    "# create loss function\n",
    "criterion = torch.nn.CrossEntropyLoss()\n",
    "\n",
    "# get dataset \n",
    "train_loader, test_loader = getData()\n",
    "\n",
    "# for illustrate, we only use one batch to do the tutorial\n",
    "for inputs, targets in train_loader:\n",
    "    break\n",
    "\n",
    "# we use cuda to make the computation fast\n",
    "model = model.cuda()\n",
    "inputs, targets = inputs.cuda(), targets.cuda()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the hessian computation module\n",
    "hessian_comp = hessian(model, criterion, data=(inputs, targets), cuda=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The top Hessian eigenvalue of this model is 195.5897\n"
     ]
    }
   ],
   "source": [
    "# Now let's compute the top eigenvalue. This only takes a few seconds.\n",
    "top_eigenvalues, top_eigenvector = hessian_comp.eigenvalues()\n",
    "print(\"The top Hessian eigenvalue of this model is %.4f\"%top_eigenvalues[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The top two eigenvalues of this model are: 105.7241 195.4954\n"
     ]
    }
   ],
   "source": [
    "# Now let's compute the top 2 eigenavlues and eigenvectors of the Hessian\n",
    "top_eigenvalues, top_eigenvector = hessian_comp.eigenvalues(top_n=2)\n",
    "print(\"The top two eigenvalues of this model are: %.4f %.4f\"% (top_eigenvalues[-1],top_eigenvalues[-2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The small difference between this top eigenvalue (195.4954) and the previous one (195.5897) is due to the small number of iterations that we used in Power iteration. You can remove this small difference by increasing the number of iterations for power iteration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2.1: Plot Loss Landscape\n",
    "\n",
    "We can use the Hessian eigenvectors/eigenvalues to analyze the flat/sharpness of the loss landscape of your model, and plot the loss landscape. We will show that this can be more informative than using random directions.\n",
    "\n",
    "To plot the loss landscape, we first compute the top Hessian eigenvector and then perturb the model parameters along that direction and measure the loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the top eigenvector\n",
    "top_eigenvalues, top_eigenvector = hessian_comp.eigenvalues()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a simple function, that will allow us to perturb the model paramters and get the result\n",
    "def get_params(model_orig,  model_perb, direction, alpha):\n",
    "    for m_orig, m_perb, d in zip(model_orig.parameters(), model_perb.parameters(), direction):\n",
    "        m_perb.data = m_orig.data + alpha * d\n",
    "    return model_perb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Loss landscape perturbed based on top Hessian eigenvector')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4MElEQVR4nO3dd3gU9fbH8fdJJ/QSeiD0Xg3dgqKiYC/Y8KJXRa5ee7n23vWKKCryw2sFFVARBSwoqCjF0DuE3gmEEgIh7fz+mEHXkE42s5uc1/Pkye7O7MxnZmfm7JSdr6gqxhhjjK8QrwMYY4wJPFYcjDHGHMeKgzHGmONYcTDGGHMcKw7GGGOOY8XBGGPMccpMcRCR90XkGT8Ne6OInOmPYZcHJTn/8vucReQ6EZlVEuMpLhHpKyJbvcwQrERklIg86sF4l4tI39Ieb6Ar8eJgG9KywzZ0gcvfn42IzBSRG0tznKo6TFWf9tfw8xlvO1WdWdrjLQn+3N6WmT0HU7JEJCwQhmGM8Y+C1s9SKw4iEikir4nIdvfvNRGJdLvVEpFvRGS/iCSLyK8iEuJ2+4+IbBORFBFZLSL9CjGu6u7wkkRkn/u4oU/3mSLytIj85g73exGp5dP9WhHZJCJ7ReThHMPuLiIJInJQRHaJyKs+3U4Wkd/d6dgiIte5rw8UkYXue7aIyBM+74kTERWRoe582SEi9/p0DxGRB0RknZtnvIjUyGO6+4rIVhF5SET2uN8qrsnxGbwiIpvd7KNEpEKO9/5HRHYCnwDTgPoicsj9q5/zsE7Ob5PuOP8jIkuAVJ8FsJuIrHA/j/dEJMrnPeeJyCJ3vv0uIh19unURkQXu5/QZ8Of78iAiMlJEDojIKt/lRUSuF5GV7rDWi8jNPt3yWwbri8jn7vK0QURu93lfBXee7BORFUC3AsL1FpE/3Hx/iEhvn275Lpc+/VUk988mv3Us32WjOESkqoi86y6z20TkGREJdbs1F5Gf3enc4352iGO4iOwWZ31YKiLt3W5/LltygutwLlnzW8b+/Pbtfp4fuONcKSL351i+81sWnhBn/fzQzbRcROLdbv8RkYk5Mo0QkdcLmpdu95t8lt0VItJVRD4CGgFfu8vA/W6/F7jj3u/OpzY5pjW39fN4qlqif8BG4MxcXn8KmAPUBmKA34Gn3W7PA6OAcPfvFECAVsAWoL7bXxzQLI/xvg884z6uCVwKRAOVgQnAJJ9+ZwLrgJZABff5C263tsAh4FQgEngVyDw2TcBs4Fr3cSWgp/u4MZACXOVOQ02gs9utL9ABpxh3BHYBF/lMk+JsjCu6/SX5jO8Od741dPO8A3ySxzzo62Z91e33NCAVaOV2Hw5MBmq48+Vr4Pkc733RfW8F97Wtec1nn/dt9Xm+EVgExAIVfF5b5r5WA/jN57PqAuwGegChwBC3/0ggAtgE3OXO08uADN/x58h2nTsNx/q/AjgA1HC7DwSa4SxbpwGHga4FLIMhwHzgMTdPU2A90N993wvAr+50xbrTuTWPfDWAfcC1QBjOsrIPqFnQcpnHZ53zs8lvHTv2+ea6bOQy/JnAjfmNE/gSZ3ms6I5zHnCz2+0T4GF3/kUBJ7uv93fnZzV3/rYB6pXkOpzLtOS5jOXcZrmf589AdZx1bsmxaS7EsvAEkAYMcMfzPDDHZ/twGKjsPg8FdvDX9iO/eXk5sA3ni4cAzYHGuW1v3fmRCpyFsxzfDyQCEXmtn3luy0uxOKwDBvg87w9s9FmovwKa53hPc/dDPRMIL2C8fy5YuXTrDOzLsWA94vP8FuBb9/FjwKc+3SoC6T4Lzy/Ak0CtHON4EPiykPPoNWC4+zgOpzi09un+EvCu+3gl0M+nWz2cDWRYHhuMTKCiz2vjgUfdhSoVn+IK9AI2+Lw3HYgqYAP0t/mcsx/38/9nLsvEMJ/nA4B17uO3cTdgPt1X42y8TgW2A+LT7fd8Pufrcul/Hm4xz6X/ScAdBSyDPYDNuXzW77mP1wPn+HQbmnOe+XS7FpiX47XZwHUFLZd5fNY5P5v81rE8l408hj8TZ2O23+fvEH9tKOsAR/HZwOAUuxnu4w+B0UDDHMM9A1gD9ARC/LEO5/LePJcxn+Xz2Pr958befX6jzzQXtCw8AUz36dYWOOLzfBbwD/fxWfy1DhQ0L7/DXU5zmbY/s7vPHwXG+zwPwSksffNaP/P6K81zDvVxvgUes8l9DeBlnOr2vTi7+w8AqGoicCfOTN8tIp+KSH0KICLRIvKOOIeGDuJs0Kv57qYBO30eH8bZCziWc8uxDqqaCuz16fcGnOq8yj0scJ77eizOyplbnh4iMsPdFT0ADANy7gJv8XnsO28aA1+6u4j7cYpFFs4ClZt9buacw4rB+RY232dY37qvH5Okqml5DLcothTwWs7pu+dYJjdXrNu9PrBN3aXa5735ya3/+gAicq6IzBHnsNF+nCJ17HPIdRl089XPke8h/pr/f1teCsiXcx041n8Dn+d5LZeFkd86BnkvG3m5XVWrHfsDzvPp1hjnm+kOn/nyDs63XnC+sQowzz3E8U8AVf0JGAm8ibNOjxaRKjlHfILrcE75LWM55fw8fR8XtCzklinK59DNOJyNPsDV7vNjw81vXua5bckj/5/LgKpmu9Pgu4zltn4epzSLw3acmXBMI/c1VDVFVe9R1abABcDd4h4rVtVxqnqy+17FOexRkHtwDkn1UNUqON9AwVlYC7ID58Nw3iASjbOLi5tnrapehfPBvQhMFOcY8BacQxa5GYdzOCdWVaviHL7ImSXW5/Gf88Yd7rm+K6mqRqnqtjzGVd3Nk3NYe4AjQDuf4VRVVd8VynejmttzcPY+on2e182ln9zel9/0PZtj+qJV9ROcz6KBiEiO9+Ynt/63i3Ps/XPgFaCOu7Gbivs55LMMbsHZu/LNV1lVB7jD/9vyUkC+nOvAsf7z+izzk9s8znMdc+W1bBTHFpxvu7V85ksVVW0HoKo7VfUmVa0P3Ay8JSLN3W6vq+pJON+sWwL35TL8E1mHc8ua1zKW0w6cw0nH+H62BS0LBZkA9HXPnVzMX8Uh33lJ/tuWnMvB35YBd12I5e/LWG7LznH8VRzCRSTK5y8M5xjkIyISI86Jo8eAj+HPk0XN3Qk5gPPNOFtEWonIGe6KnYazccsuxPgru/3uF+fk7eNFyD4ROE+ck8sROIcb/pxPIjJYRGLcirzffTkbGAucKSKDRCRMRGqKSGefPMmqmiYi3XG+NeT0qPttqR1wPfCZ+/oo4FkRaeyOP0ZELixgGp4UkQgROQXn294EN+//AcNFpLY7rAYi0j+f4ewCaopIVZ/XFgEDRKSGiNTF2bMrjFtFpKH7eTzsM33/Bwxz965ERCqKcwK/Ms4hl0zgdhEJF5FLgO4FjKe2T/+X4xzTnopzjDgS53xOpoicC5x97E15LYM4h6VSxDmJV0FEQkWkvYgcO/E8HnhQnBOoDYHb8sk2FWgpIle7y8gVOBvIbwox/3LK7bPJcx3zcdyyUYxxo6o7gO+B/4pIFXEunGgmIqcBiMjl8tcJ5H04G6RsEenmftbhOF800sh9nT6RdTin/JaxnHw/zwbAv326FbQs5EtVk3AOh72HU2RWuq/nOy+BMcC9InKSm7/5se0BznLQNEf+gSLSz53H9+AUnt8Lk9GXv4rDVJwP9tjfE8AzQALOCZ6lwAL3NYAWwHScY5qzgbdUdQbOyvwCzrfenTgr/oOFGP9rOCep9uCcoPu2sMFVdTlwK05V34GzYPte230OsFxEDgEjgCtV9YiqbsY5THEPkIyzEe3kvucW4CkRScFZYcfnMuqfcQ5r/Ai8oqrfu6+PwNnr+N59/xycY5952elm3o5TsIap6iq323/cccwRZ1d9Os63s7zmxSqcDc56cXZ36wMfAYtxjl1+z18b+YKMc/tfj7OL/Iw7jgTgJpxDDfvcfNe53dKBS9znyTgnmL8oYDxzcZanPcCzwGWquldVU4Dbceb9PpwCPdnnfbkug6qahbMR7QxscIc7Bji2UX4SZzd+gzt9H+UVTFX3usO6B+dQ5f3Aeaq6p4Bpym1YuX02+a1jkP+yURz/wCm6K9zhTsQ5JwbOydO57noyGeeY+XqgCs7Geh/OfNuLc0gvp9co5jqcU37LWC6ewlnfN+AsDxNxNq4UYlkojHE451DH5Xg9z3mpqhNwluVxOBe9TMK5uAGck96PuMvAvaq6GhgMvOHmOx84312XikT+fnjWlDYRicNZ0MJVNfMEh9UX+FhVGxbQqylnbNkoHhH5F84XwNMK7LmMsR/BGWOMS0TqiUgf99BOK5y9vC+9zuUF+wWrMcb8JQLnSqEmOOcUPwXe8jKQV+ywkjHGmOP49bCSiFQTkYni3MZgpYj08uf4jDHGlAx/H1YagfOrxcvcy0Kj8+u5Vq1aGhcX5+dIxhhTdsyfP3+PqsYU3GfR+K04uNdfn8rfL0vM93KquLg4EhIS/BXJGGPKHBEp6K4BxeLPw0pNcH5w9J44dyQdI3//dSYA4tyNNEFEEpKSkvwYxxhjTGH5sziEAV2Bt1W1C86vIR/I2ZOqjlbVeFWNj4kp8T0jY4wxxeDP4rAV526Gc93nE3GKhTHGmADnt+KgqjuBLe4PSQD64fw03BhjTIDz99VKtwFj3SuV1uPcUM4YY0yA82txUNVFQLw/x2GMMabk2b2VjDHGHCfoi0N2tvLmjESWbN3vdRRjjCkzgr44pKRlMnbOJm4Zu4ADhzO8jmOMMWVC0BeHqtHhjLymKzsPpHHvxMXYjQSNMebEBX1xAOjaqDoPDmjDDyt28e6sDV7HMcaYoFcmigPAP/vEcU67urwwbRXzNyV7HccYY4JamSkOIsJLl3ekfrUK3Dp2IXsPHfU6kjHGBK0yUxwAqkSF89Y1XUk+nM5d4xeTnW3nH4wxpjjKVHEAaN+gKo+f35Zf1iTx1sxEr+MYY0xQKnPFAeDq7o24sHN9Xv1hDb+v2+N1HGOMCTplsjiICM9d3IEmtSpy+yeL2H0wzetIxhgTVMpkcQCoGBnG24NP4tDRDG77ZCGZWdleRzLGmKBRZosDQMs6lXnmog7M3ZDMa9PXeh3HGGOCRpkuDgCXndSQK+JjGTkjkRmrd3sdxxhjgkKZLw4AT17YjtZ1K3PXZ4vYvv+I13GMMeZPaRlZbEk+7HWM45SL4hAVHspb13QlM0u5ddwC0jPt/IMxxnvpmdncMnYBl779O4eOZnod52/KRXEAaBpTiRcu7cDCzft58dtVXscxxpRzWdnKXeMX8dOq3dxxZgsqRfq7Yc6iKTfFAeC8jvUZ0qsx787awLfLdnodxxhTTmVnKw9+sYQpS3bw0IDWXNOjsdeRjlOuigPAQwPb0KlhVe6buJhNe1O9jmOMKWdUlaenrGB8wlZu79eCoac28zpSrspdcYgMC2Xk1V0R4NZxC0jLyPI6kjGmHBk+fS3v/baRf/Zpwl1ntvA6Tp7KXXEAiK0RzauDOrNs20Ge/maF13GMMeXE6F/W8fqPa7kiPpZHz2uDiHgdKU/lsjgAnNm2Djef1pSxczfz1aJtXscxxpRxY+du4rmpqxjYsR7PXdIhoAsDlOPiAHDv2a3oFledB79YSuLuFK/jGGPKqEkLt/HIpGWc0bo2wwd1JjQksAsDlPPiEB4awhtXdaVCeCjDPl5AaoBdZ2yMCX7fL9/JPRMW06NJDd66pisRYcGx2Q2OlH5Ut2oUr1/VhXVJh3j4y6WoWgNBxpiSMWvtHv49biHtG1RlzJBuRIWHeh2p0PxaHERko4gsFZFFIpLgz3GdiD7Na3H3mS2ZtGg74+Zt9jqOMaYMmL8pmZs+TKBpTEU+uL5bwP3IrSClkfZ0VQ34FnduPb05CZv28eTkFXRsUI0ODat6HckYE6SWbTvAde/9Qd2qUXx4Q3eqRUd4HanIyv1hpWNCQoThV3SmVqUI/jV2PgcOZ3gdyRgThBJ3p/CP/82jcmQYH9/Yg9qVo7yOVCz+Lg4KfC8i80VkaG49iMhQEUkQkYSkpCQ/x8lfjYoRjLymK7sOpnHPhEVkZ9v5B2NM4W1JPszgMfMIEWHsTT1pUK2C15GKzd/F4WRV7QqcC9wqIqfm7EFVR6tqvKrGx8TE+DlOwbo2qs7DA9owfeVuRv+63us4xpggsetgGteMmcuRjCw+uqE7TWpV9DrSCfFrcVDVbe7/3cCXQHd/jq+kDOkdx8AO9Xj5u9XMWb/X6zjGmACXnJrONWPmsvfQUd6/vhtt6lXxOtIJ81txEJGKIlL52GPgbGCZv8ZXkkSEFy7tQOMa0dz2yUJ2p6R5HckYE6CSU9MZPGYuW5IPM2ZIN7o0qu51pBLhzz2HOsAsEVkMzAOmqOq3fhxfiaocFc5bg7uSkpbB7Z8sJDPLGggyxvzd7pQ0rhw9m3VJh3jn2pPo1aym15FKjN+Kg6quV9VO7l87VX3WX+Pyl9Z1q/DsRR2Ysz6Z4dPXeB3HGBNAdhw4wpXvzGFL8hHeu64bfVvV9jpSibJLWQtw6UkNuap7LG/OWMdPq3Z5HccYEwC2JB9m0Duz2Z1ylI9u6E7v5rW8jlTirDgUwuPnt6NtvSrc9dnigGwI3BhTejbsSWXQO7M5eCSTsTf2ID6uhteR/MKKQyFEhYfy9uCuZKty67gFHM20BoKMKY/W7kph0DuzOZqZzSc39aRTbDWvI/mNFYdCalyzIq9c3oklWw/w7JSVXscxxpSy5dsPcMXoOQB8NrQnbesH/+Wq+bHiUAT929Vl6KlN+XD2JmsgyJhyZNGW/Vw1eg5RYSGMv7kXLepU9jqS31lxKKL7+lsDQcaUJ39sTGbwmLlUjQ7ns5t7Bf0vnwvLikMR+TYQ9K+PF3A43RoIMqas+j1xD/94dx61K0cy4ebexNaI9jpSqbHiUAzHGghKTDrEQ19YA0HGlEUzVu/m+vf/ILZGBT69uSd1qwbn3VWLy4pDMfVpXou73AaCPvtji9dxjDEl6LvlOxn6YQLNa1fi06G9gva22yfCisMJuPX05pzcvBaPT17O6p12/sGYsuDrxdu5ZewC2tWvyribelKjYvA11FMSrDicgFC3gaDKUeHcOs7OPxgT7CbO38odny7kpMbV+fjGHlStEO51JM9YcThBMZUjGXFlZ9YlHeLRScu9jmOMKaaxczdx74TF9G5Wiw+u7x50bT6XNCsOJaBP81rcdkYLPl+wlYnzt3odxxhTRB/8vpGHv1zGGa1rM2ZIPBUiQr2O5DkrDiXkjn4t6Nm0Bo9OWma/fzAmiLz32wYen7ycs9rWYdTgk4gKt8IAVhxKTGiIMOLKLkRHhHLr2IUcSbf7LxkT6Mb8up4nv15B/3Z1ePPqrkSE2SbxGJsTJahOlShevaIzq3el8OTXdv7BmED2f7+s55kpKxnQoS4jrTAcx+ZGCTutZQy39G3Gp39ssfsvGROgRv28jmenrmRgx3qMuLIL4aG2KczJ5ogf3H1WS+IbV+ehL5ayPumQ13GMMT7enJHIC9NWcX6n+oy4orMVhjzYXPGDsNAQXr+qC+FhIdw6biFpGXb+wZhAMPKntbz83Wou7Fyf4YM6EWaFIU82Z/ykfrUKvDqoEyt3HOSZKSu8jmNMuTdi+lpe+X4Nl3RpwKuDOlthKIDNHT86o3Udhp7alI/nbGbKkh1exzGmXFJVhv+whuHT13Bp14a8fHknQkPE61gBz4qDn93XvxVdGlXjgc+XsGlvqtdxjClXjhWGET+u5fKTGvLSZR2tMBSSFQc/Cw8N4fUruyAC/x630NqfNqaUqCqvfL+a139K5Mpusbx4qRWGorDiUApia0Tz8uWdWLrtAM9PXeV1HGPKPFXlxW9X8+aMdVzVvRHPXdyBECsMRWLFoZT0b1eX6/vE8f7vG/lu+U6v4xhTZqkqL0xbxaif13FNj0Y8e1F7KwzFYMWhFD14bhs6NqzKfRMWsyX5sNdxjClzVJVnp6zknV/Wc23PxjxjhaHY/F4cRCRURBaKyDf+HlegiwgLYeRVXVGF2z5ZSHpmtteRjCkzVJWnv1nJmFkbuK53HE9d2A4RKwzFVRp7DncAK0thPEGhUc1oXrysI4u27OeV71d7HceYMuPFb1fzv982cH2fOB4/v60VhhPk1+IgIg2BgcAYf44n2AzoUI9rezZm9C/rmbl6t9dxjAl6b81MZNTP6xjcsxGPnWeFoST4e8/hNeB+IM/jJyIyVEQSRCQhKSnJz3ECx8MD29CyTiXun7iEfanpXscxJmh9NGcTL33r3BLjqQvaW2EoIX4rDiJyHrBbVefn15+qjlbVeFWNj4mJ8VecgBMVHsqrgzqz73A6D09aiqp6HcmYoDNp4TYe+2oZZ7apzSuXd7KTzyXIn3sOfYALRGQj8Clwhoh87MfxBZ32Dapy11ktmbp0J5Ps9t7GFMn0Fbu4Z8Jiejapyciru9rdVUuY3+amqj6oqg1VNQ64EvhJVQf7a3zB6uZTmxHfuDqPTVrOtv1HvI5jTFD4fd0ebhm3gPb1q/B/Q+KtaU8/sFLrsdAQ4dVBnclW5Z7xi8jOtsNLxuRn0Zb93PRBAnE1o3n/+u5UigzzOlKZVCrFQVVnqup5pTGuYNSoZjSPn9+OOeuT+d9vG7yOY0zAWr0zhevem0fNSpF8dEMPqleM8DpSmWV7DgHi8viGnNW2Di99u5rVO1O8jmNMwNm0N5Vr351LZFgIY2/sQZ0qUV5HKtOsOAQIEeH5SzpQpUIYd362yO7eaoyPnQfSGPzuXDKysvn4hh7E1oj2OlKZZ8UhgNSqFMkLl3Rk5Y6DDP9hrddxjAkIyanpXPvuXPalZvDBP7vTok5lryOVC1YcAsyZbetwVfdY3vllHfM2JHsdxxhPpaRlcN1789icfJgxQ+Lp2LCa15HKDSsOAeiRgW2JrR7N3eMXkZKW4XUcYzyRlpHFjR8ksGL7Qd4e3JWeTWt6HalcseIQgCpGhjH8ik5s33+Ep75e4XUcY0pdRlY2t4xdwLyNyfx3UCfOaF3H60jljhWHAHVS4xrc0rc5E+Zv5dtl1jiQKT+yspW7xy/mp1W7eeai9lzYuYHXkcolKw4B7PZ+LWjfoAoPfbmU3SlpXscxxu9UlUe/WsbXi7fzwLmtuaZHY68jlVtWHAJYRFgIwwd1JvVoJg98bjfnM2XfGz8lMm7uZm7p24xhpzXzOk65ZsUhwLWoU5kHzm3NT6t288m8LV7HMcZvZq7ezfDpa7i4SwPu69/K6zjlnhWHIDCkVxwnN6/F09+sYOOeVK/jGFPitu47zJ2fLaJVnco8d3EHa5MhAFhxCAIhIcLLl3ckPFS4a/wiMrOs7WlTdhzNzOKWsQvIylLeHnwSFSLsDquBwIpDkKhXtQLPXNyBhZv38/bMdV7HMabEPPn1CpZsPcArgzrRpFZFr+MYlxWHIHJBp/pc0Kk+I35cy5Kt+72OY8wJmzh/K+PmbmbYac3o366u13GMDysOQebpC9tTq1Ikd322iCPpdnM+E7xWbD/Iw18upWfTGtx7dkuv45gcrDgEmarR4fx3UCfWJaUyfPoar+MYUywHjmTwr7HzqRYdzhtXdSXMmvgMOPaJBKE+zWtxZbdY3p21gRXbD3odx5giyc5W7hm/mG37jvDm1V2JqRzpdSSTCysOQeqBc1tTrUI4D09aak2LmqAy6pd1TF+5i4cGtCE+robXcUwerDgEqWrRETw8sA0LN+9n3LzNXscxplB+T9zDK9+t5ryO9bi+T5zXcUw+rDgEsYu7NKB3s5q8+O0qu/eSCXg7Dhzhtk8W0jSmEi9e2tF+6BbgrDgEMRHh6YvaczQjm2e+Wel1HGPylJ6Zza1jF5CWkcWowSdRMTLM60imAFYcglyzmEr8q28zJi/ezi9rkryOY0yunpu6kgWb9/PiZR1pXruS13FMIVhxKAP+1bcZTWpV5NGvlpGWYb99MIHlq0XbeP/3jfyzTxPO61jf6zimkKw4lAFR4aE8e1F7Nu09zJszEr2OY8yf1uxK4YHPlxLfuDoPDmjtdRxTBIUqDiJSUURC3MctReQCEQn3bzRTFL2b1+LiLg0Y9fM6EneneB3HGFLSMhj28XwqRobx5jVdCbcfugWVwn5avwBRItIA+B64FnjfX6FM8Tw8sA3REWE8/OUyaxjIeEpV+c/nS9i09zAjr+5CnSpRXkcyRVTY4iCqehi4BHhLVS8H2uX7BpEoEZknIotFZLmIPHmiYU3+alWK5IFzWzN3QzIT52/1Oo4px96dtYGpS3dyf/9W9Gxa0+s4phgKXRxEpBdwDTDFfa2gm64fBc5Q1U5AZ+AcEelZrJSm0K6IjyW+cXWem7qS5NR0r+OYcmjehmSen7aKc9rVZeipTb2OY4qpsMXhTuBB4EtVXS4iTYEZ+b1BHYfcp+Hunx3r8LOQEOHZizuQkpbJ81Pttw+mdKVlZHHPhEXEVq/AS5fbD92CWaGKg6r+rKoXqOqL7onpPap6e0HvE5FQEVkE7AZ+UNW5ufQzVEQSRCQhKcmu0y8JrepW5sZTmjJh/lbmrt/rdRxTjrw1I5EtyUd47pIOVImya1aCWWGvVhonIlVEpCKwDFghIvcV9D5VzVLVzkBDoLuItM+ln9GqGq+q8TExMUWMb/JyR78WNKxegYcnLSM905oVNf63PukQo35ez0Wd69O7WS2v45gTVNjDSm1V9SBwETANaIJzxVKhqOp+nMNQ5xQxnymmChGhPH1hexJ3H2L0L9asqPEvVeWxr5YTGR7CQwPbeB3HlIDCFodw93cNFwGTVTWDAs4fiEiMiFRzH1cAzgJWFT+qKarTW9dmQIe6vPFTIpv2pnodx5RhU5buYFbiHu7r34rale2y1bKgsMXhHWAjUBH4RUQaAwW1MlMPmCEiS4A/cM45fFPcoKZ4Hj+/HeGhITwyyX77YPwjJS2Dp75eQfsGVbimR2Ov45gSUtgT0q+ragNVHeBehbQJOL2A9yxR1S6q2lFV26vqUyWS2BRJnSpR3Ht2S35du4evl+zwOo4pg16bvpakQ0d55qIOhIbY1UllRWFPSFcVkVePXVUkIv/F2YswQeDaXnF0bFiVp75ewYEjGV7HMWXIiu0Hef/3jVzdvRGdY6t5HceUoMIeVvofkAIMcv8OAu/5K5QpWaEhwnMXdyA59Sgvf2enfUzJyM5WHpm0lGoVwrm/v91Ur6wpbHFopqqPq+p69+9JwH76GETaN6jKkN5xjJ27mYWb93kdx5QBE+ZvYcHm/Tw4oA1Vo+03DWVNYYvDERE5+dgTEekDHPFPJOMv95zdijqVo3joy2VkZtlvH0zxJaem8/y0VXSPq8GlXRt4Hcf4QWGLwzDgTRHZKCIbgZHAzX5LZfyiUmQYT1zQlpU7DvLebxu9jmOC2EvfriIlLZOnL2pvt8goowp7tdJi9wZ6HYGOqtoFOMOvyYxf9G9Xl9NbxTDix7UkpRz1Oo4JQvM37ePTP7Zww8lNaFW3stdxjJ8UqfUNVT3o/lIa4G4/5DF+JiI8el5b0jKyeOW71V7HMUEmMyubRyYto26VKO7o18LrOMaPTqRpJtuXDFJNYypxfZ84xs/fwrJtB7yOY4LIh7M3sXLHQR4/vy0VI8O8jmP86ESKg/3cNojd1q8FNaIjePLr5fbLaVMouw6m8eoPazitZQzntK/rdRzjZ/kWBxFJEZGDufylAPVLKaPxgypR4dzbvxV/bNzHN/bLaVMIz0xZSXpWNk9d2M5OQpcD+RYHVa2sqlVy+ausqrZPGeQGxcfStl4Vnp+6kiPpWV7HMQHs17VJfL14O7f2bU7jmnZzhPLgRA4rmSAXGiI8fn5bth9IY/Qv672OYwLU0cwsHvtqOXE1o7n5NPvta3lhxaGc69G0JgM71uPtnxPZvt9+12iON/rn9WzYk8pTF7YnKrygpuNNWWHFwfDgua1RhRem2X2XzN9t3nuYkTMSGdihHqe2tJYayxMrDoaG1aO5+dSmTF68nYSNyV7HMQFCVXl88jLCQpzfxpjyxYqDAWBY32bUrRLFk1+vIDvbLm018P2KXcxYncRdZ7WkblVr3a28seJgAIiOCOOBc1uzdNsBJi7Y6nUc47HD6Zk8OXk5retWZkjvOK/jGA9YcTB/urBzfbo2qsZL364mJc0aBSrPXv8xke0H0nj6ovaEh9pmojyyT938SUR4/Px27Dl0lJEzEr2OYzyyff8R/jdrA5d2bUi3uBpexzEeseJg/qZTbDUuO6kh783ayMY9qV7HMR54c0YiinL32S29jmI8ZMXBHOf+/q0IDxWenbrS6yimlG3bf4TxCVsYFB9Lg2oVvI5jPGTFwRyndpUobj2jOT+s2MWva5O8jmNK0Zvu4cRbTm/ucRLjNSsOJlf/7NOERjWiefqbFdakaDmxdd9hJiRs4YputtdgrDiYPESFh/LwwDas2XWIcfM2ex3HlII3ZyQiCLfaXoPBioPJx9lt69C7WU3++/0a9qWmex3H+NGW5MNMSNjKld1jqVfV9hqMH4uDiMSKyAwRWSEiy0XkDn+Ny/iHiPDY+W1JScvgtelrvI5j/GjkT4mEiHBLX9trMA5/7jlkAveoalugJ3CriNgNWoJM67pVuKZHYz6eu5k1u1K8jmP8YPPew3y+YCtXdY+122SYP/mtOKjqDlVd4D5OAVYCDfw1PuM/d5/VkkqRYTz9zQprUrQMGjljLSEhYlcomb8plXMOIhIHdAHm5tJtqIgkiEhCUpJdNhmIqleM4M4zW/Dr2j1MX7nb6zimBG3am8rnC7ZxdfdG1Kliew3mL34vDiJSCfgcuFNVD+bsrqqjVTVeVeNjYux+8YFqcM/GNK9diWemrOBopjUpWla88VMiYSHCLX2beR3FBBi/FgcRCccpDGNV9Qt/jsv4V3hoCI+e15ZNew/z3m8bvY5jSsDGPal8uXAb1/RoTG3bazA5+PNqJQHeBVaq6qv+Go8pPae1jKFf69qM/CmRPYeOeh3HnKBjew3DrF1okwt/7jn0Aa4FzhCRRe7fAD+Oz5SChwa2IS0ji+E/2KWtwWzDnlS+XLiVwT1tr8HkLsxfA1bVWYD4a/jGG81iKjG4Z2M+nL2RIb3jaFmnsteRTDG88eNaIsJCGHaanWswubNfSJsiu6NfCypFhvHsFLtrazBal3SISYu2cW3PxsRUjvQ6jglQVhxMkVWvGMHt/Vrw85okfl5jlx8Hmzd+XEtkWCg3216DyYcVB1Ms1/ZqTOOa0Tw7xe7aGkwSdx9i8uLt/KNXY2pVsr0GkzcrDqZYIsNCeeCc1qzZdYjxCVu9jmMK6Y2fnL2Gm061K5RM/qw4mGI7p31dusVV59UfVpOSluF1HFOAxN0pzl5Db9trMAWz4mCKTUR4ZGBb9hxK5+2Z67yOYwow4sdEKoSHcvOpdq7BFMyKgzkhnWKrcXGXBoyZtYGt+w57HcfkYc2uFL5Zsp0hveOoUTHC6zgmCFhxMCfsvv6tEODl71Z7HcXk4fUf1xIdHspNp9i5BlM4VhzMCatfrQI3ndKUrxZtZ+HmfV7HMTms2ZXClKU7bK/BFIkVB1MihvVtRq1KkTwzZaW1+RBgRkxfS8WIMNtrMEVixcGUiEqRYdx7dkvmb9rHtGU7vY5jXKt2HmTK0h1c1zuO6rbXYIrAioMpMZfHx9K6bmWen7bS2nwIECOmr6VyZBg3ntLE6ygmyFhxMCUmNMS5tHVL8hE++H2j13HKvRXbDzJt2U6u7xNHtWjbazBFY8XBlKiTW9TijNa1eePHRPZamw+eev1HZ6/hhpPtXIMpOisOpsQ9NKA1hzOyGPHjWq+jlFvLth3g2+U7uf7kJlSNDvc6jglCVhxMiWteuzJXd2/E2LmbSdx9yOs45dIr36+maoVwO9dgis2Kg/GLO89sQXR4KM9PtTYfStu8DcnMXJ3Ev/o2o0qU7TWY4rHiYPyiZqVI/n1Gc35ctZtZa/d4HafcUFVe/m4VtStHMqRXnNdxTBCz4mD8ZkjvOGJrVOCZKSvIyrYfxpWGmauT+GPjPm7r14IKEaFexzFBzIqD8Zuo8FD+c05rVu1M4fP51uaDv2VnKy99t5pGNaK5Ij7W6zgmyFlxMH41sEM9ujaqxsvfryb1aKbXccq0KUt3sHLHQe4+qyURYbZqmxNjS5DxKxHhkfPakpRylHd+tjYf/CUjK5tXf1hDqzqVOb9Tfa/jmDLAioPxu66NqnNBp/qM/nU92/cf8TpOmfT5/K1s2JPKvf1bERoiXscxZYAVB1Mq7j+nFdkKr1ibDyUuzf3BYZdG1TizTW2v45gywoqDKRUNq0dz48lN+GLhNru0tYR9PGcTOw6kOY0uie01mJJhxcGUmtvOaEGzmIrcM2ER+1LTvY5TJqSkZfDmjEROaVGL3s1qeR3HlCF+Kw4i8j8R2S0iy/w1DhNcKkSEMuLKLiSnpvPAF0usUaAS8O6sDew7nMG9Z7fyOoopY/y55/A+cI4fh2+CUPsGVbmvfyu+W76L8QlbvI4T1JJT0xnz6wbOaVeXTrHVvI5jyhi/FQdV/QVI9tfwTfC68eSm9G5Wkycmr2DDnlSv4wStt2cmcjg9k3vObul1FFMGeX7OQUSGikiCiCQkJSV5HceUgpAQ4b+DOhERFsKdny4kIyvb60hBZ8eBI3wwexMXd2lIizqVvY5jyiDPi4OqjlbVeFWNj4mJ8TqOKSX1qlbghUs6sHjrAUZMt3Yfiur1HxNRVe48s4XXUUwZ5XlxMOXXuR3qMSi+IW/OTGTu+r1exwkaG/akMj5hC1d3b0RsjWiv45gyyoqD8dTj57ejUY1o7h6/mANHMryOExRe/WENEaEh/PsM22sw/uPPS1k/AWYDrURkq4jc4K9xmeBVMTKM167ozM6DaTz2lV31XJAV2w/y9eLt/PPkOGIqR3odx5Rh/rxa6SpVraeq4araUFXf9de4THDr0qg6d/ZrwVeLtjNp4Tav4wS0V75fTZWoMIae0szrKKaMs8NKJiDccnpzusVV59FJy9iSfNjrOAEpYWMyP63azbC+zagabc1/Gv+y4mACQmiI8OqgzgDc9dkiMu3y1r9RdRryiakcyXW947yOY8oBKw4mYMTWiObpi9qTsGkfb8+0th98/bwmiXkbkrntjOZER4R5HceUA1YcTEC5qEsDLuxcn9d+XMvCzfu8jhMQsrOVl79bTcPqFbiyWyOv45hywoqDCThPXdieulWiuPOzRRyypkWZtmwny7cf5K4zrflPU3psSTMBp2qFcIZf0ZktyYd5cvJyr+N4KjMrm//+sJqWdSpxUZcGXscx5YgVBxOQujepwS19mzNh/lamLt3hdRzPfLFgG+uTUrnnbGv+05QuKw4mYN1xZgs6NazKg18sZceB8tf2dFpGFq9NX0On2Gqc3baO13FMOWPFwQSs8NAQXruyCxlZ2dz92WKys8tP40AZWdncN3EJ2w+kcb81/2k8YMXBBLQmtSry+Pltmb1+L//363qv45SKtIws/vXxAr5evJ0Hzm1Nn+bW/KcpfVYcTMAbFB/LOe3q8sr3q5lWxs8/HE7P5KYPE5i+chdPX9iOYafZbTKMN6w4mIAnIjx/SQfa1qvCv8Yu4MEvlnAkPcvrWCUuJS2DIf+bx2+Je3jl8k5c2yvO60imHLPiYIJC9YoRTBjWm5tPa8on87Zw/shZrNh+0OtYJWZfajrXjJnLws37ef2qLlx2UkOvI5lyzoqDCRoRYSE8eG4bPr6hBweOZHDRm7/x3m8bUA3uE9W7U9K4cvQcVu1M4Z1rT+K8jvW9jmSMFQcTfE5uUYtv7ziFk1vU4smvV3DDBwnsPXTU61jFsn3/Ea54Zw6bkw/z3nXd6NfGLlk1gcGKgwlKNStF8u6QeJ44vy2zEvdwzohf+XVtkteximTT3lQuHzWbPSlH+eiG7nZVkgkoVhxM0BIRruvThK9u7UPVCuFc++48np+6kvTMwL/d99pdKVw+ajaH0zP5ZGhP4uNqeB3JmL+x4mCCXpt6Vfj63ydzVfdGvPPLei4b9Tsb96R6HStPy7Yd4IrRc1Dg06G9aN+gqteRjDmOFQdTJlSICOX5SzowanBXNu09zMDXf+Xz+VsD7mT1/E37uOr/5hAVFsL4m3vRqm5lryMZkysrDqZMOad9PabdcQrtGlTlngmLufOzRaSkZXgdC4DZ6/Zy7btzqVkxgvHDetGkVkWvIxmTJysOpsypX60Cn9zUk7vPask3S3Yw4PVfWeBxw0EzVu3muvfm0bB6Bcbf3IuG1aM9zWNMQaw4mDIpNES4vV8Lxt/ck+xsuHzUbO76bBE/rNhFWkbp/rp62tIdDP0ogRZ1KvHp0F7UrhJVquM3pjgkkI7JxsfHa0JCgtcxTBlz4EgGL3+3ismLtnMwLZNKkWH0a1Obc9vXo2+rGKLCQ0t8nPtS05m7YS+zEvcwbu5mujSqznvXd6NKVHiJj8uUbyIyX1XjS3y4VhxMeZGemc3s9XuZtnQH3y3fyb7DGURHhHJ669oMaF+P01vHEB0RVqxhp6RlMG9DMrPX7eX3dXtZufMgqhAdEcpZbevw3MUdqBhZvGEbkx8rDsaUoMysbOZuSGaqWyj2HEonKjyE01vV5twO9TijdW0q5bMxP5yeScLGffy+bi+z1+9l2bYDZGUrEWEhnNSoOr2b1aRXs5p0bFjN2n02fhWUxUFEzgFGAKHAGFV9Ib/+rTgYL2RlK/M2JDNt2Q6mLdtJUspRIsJCOK1lDAM61KVfmzpEhIawcPN+Zq/bw+z1e1m0ZT8ZWUpYiNA5thq9m9WkZ7OadG1U3S+HqYzJS9AVBxEJBdYAZwFbgT+Aq1R1RV7vseJgvJadrczfvI+pS3cwbelOdh5MIyI0BBE4mplNiECHBlXp1awWvZrVJL5xdTtcZDzlr+Lgz6W6O5CoqusBRORT4EIgz+JgjNdCQoRucTXoFleDRwe2ZdHW/Xy3bCeZ2UqvpjXp3rSGnVQ25YI/i0MDYIvP861ADz+Oz5gSFRIidG1Una6NqnsdxZhS5/mZMhEZKiIJIpKQlBRcd9U0xpiyyp/FYRsQ6/O8ofva36jqaFWNV9X4mJgYP8YxxhhTWP4sDn8ALUSkiYhEAFcCk/04PmOMMSXEb+ccVDVTRP4NfIdzKev/VHW5v8ZnjDGm5Pj1GjxVnQpM9ec4jDHGlDzPT0gbY4wJPFYcjDHGHMeKgzHGmOME1I33RCQJ2OR1jiKqBezxOkQps2kuH2yag0NjVS3x3wEEVHEIRiKS4I/7mgQym+bywaa5fLPDSsYYY45jxcEYY8xxrDicuNFeB/CATXP5YNNcjtk5B2OMMcexPQdjjDHHseJgjDHmOFYcikhEaojIDyKy1v2fZ0swIlJFRLaKyMjSzFjSCjPNItJZRGaLyHIRWSIiV3iR9USJyDkislpEEkXkgVy6R4rIZ273uSIS50HMElWIab5bRFa4n+uPItLYi5wlpaDp9envUhFRESmXl7ZacSi6B4AfVbUF8KP7PC9PA7+USir/Ksw0Hwb+oartgHOA10SkWulFPHFuu+dvAucCbYGrRKRtjt5uAPapanNgOPBi6aYsWYWc5oVAvKp2BCYCL5VuypJTyOlFRCoDdwBzSzdh4LDiUHQXAh+4jz8ALsqtJxE5CagDfF86sfyqwGlW1TWqutZ9vB3YDQRb601/tnuuqunAsXbPffnOi4lAPxGRUsxY0gqcZlWdoaqH3adzcBruClaF+YzB+WL3IpBWmuECiRWHoqujqjvcxztxCsDfiEgI8F/g3tIM5kcFTrMvEekORADr/B2shOXW7nmDvPpR1UzgAFCzVNL5R2Gm2dcNwDS/JvKvAqdXRLoCsao6pTSDBRq/tucQrERkOlA3l04P+z5RVRWR3K4FvgWYqqpbg+VLZQlM87Hh1AM+AoaoanbJpjReEpHBQDxwmtdZ/MX9YvcqcJ3HUTxnxSEXqnpmXt1EZJeI1FPVHe6GcHcuvfUCThGRW4BKQISIHFLV/M5PeKoEphkRqQJMAR5W1Tl+iupPhWn3/Fg/W0UkDKgK7C2deH5RqLbeReRMnC8Kp6nq0VLK5g8FTW9loD0w0/1iVxeYLCIXqGpCqaUMAHZYqegmA0Pcx0OAr3L2oKrXqGojVY3DObT0YSAXhkIocJrddsK/xJnWiaWYrSQVpt1z33lxGfCTBvcvSQucZhHpArwDXKCquX4xCCL5Tq+qHlDVWqoa566/c3Cmu1wVBrDiUBwvAGeJyFrgTPc5IhIvImM8TeY/hZnmQcCpwHUissj96+xJ2mJyzyEca/d8JTBeVZeLyFMicoHb27tATRFJBO4m/6vVAl4hp/llnD3gCe7nmrNgBo1CTq/Bbp9hjDEmF7bnYIwx5jhWHIwxxhzHioMxxpjjWHEwxhhzHCsOxhhjjmPFwQQVEclyL6dcJiITRCS6CO+NE5GriznejSJSqwj99xWR3j7Ph4nIP4ozbmO8YMXBBJsjqtpZVdsD6cCwwrzJ/TVzHFDk4uDeybOo+gJ/FgdVHaWqHxZjOMZ4woqDCWa/As1FpKKI/E9E5onIQhG5EEBErhORySLyE86txl/Aua3JIhG5y+3+Z1sbIvKNiPR1Hx8Skf+KyGKc26EA3C8iS93xNHf7O99t12GhiEwXkTpuGw/DgLvccZ0iIk+IyL3uezqLyBy3fYQvxW0fQ0RmisiL7vDXiMgppTAPjcmVFQcTlNw9gXOBpTj3/PlJVbsDpwMvi0hFt9euwGWqehrOr5l/dfc8hhcwiorAXFXtpKqz3NcOqGoHYCTwmvvaLKCnqnbBuf3z/aq6ERgFDHfH9WuOYX8I/MdtH2Ep8LhPtzB3Ou7M8boxpcpuvGeCTQURWeQ+/hXndha/Axcc+2YORAGN3Mc/qGpyMcaTBXye47VPfP4fKy4Ngc/cGxJGABvyG6iIVAWqqerP7ksfABN8evnC/T8f5zCYMZ6w4mCCzRFV7ez7gtvYzqWqujrH6z2A1HyGlcnf956jfB6nqWpWjv41l8dvAK+q6mT3kNQTBeQvyLE7nmZh66fxkB1WMmXBd8Btx1pkc+8impsUnFsyH7MR6CwiISISi9NKWH6u8Pk/231clb9u+TzEp9+c4wKcu34C+3zOJ1wL/JyzP2O8Zt9MTFnwNM45gCVuYy0bgPNy6W8JkOWeZH7ffc8GYAXOHToXFDCe6iKyBOfb/VXua0/g3K10H/AT0MR9/Wtgonty/LYcwxkCjHIvw10PXF+YiTSmNNldWY0xxhzHDisZY4w5jhUHY4wxx7HiYIwx5jhWHIwxxhzHioMxxpjjWHEwxhhzHCsOxhhjjvP/CD7KJtZWUNMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# lambda is a small scalar that we use to perturb the model parameters along the eigenvectors \n",
    "lams = np.linspace(-0.5, 0.5, 21).astype(np.float32)\n",
    "\n",
    "loss_list = []\n",
    "\n",
    "# create a copy of the model\n",
    "model_perb = ptcv_get_model(\"resnet20_cifar10\", pretrained=True)\n",
    "model_perb.eval()\n",
    "model_perb = model_perb.cuda()\n",
    "\n",
    "for lam in lams:\n",
    "    model_perb = get_params(model, model_perb, top_eigenvector[0], lam)\n",
    "    loss_list.append(criterion(model_perb(inputs), targets).item())\n",
    "\n",
    "plt.plot(lams, loss_list)\n",
    "plt.ylabel('Loss')\n",
    "plt.xlabel('Perturbation')\n",
    "plt.title('Loss landscape perturbed based on top Hessian eigenvector')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare this with a loss landscape computed based on perturbing the model parameters along a random direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Loss landscape perturbed based on a random direction')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzvUlEQVR4nO3dd5wU9f3H8deHXqT33ouIAnICdmKJgEaNIXbFEg1GYiNGTSyIKXaMP4yIsXfEhkbEgiIqEI9eBIUD5Oi9H9c+vz9mTtd1rwC7t1fez8fjHjflOzOf7+7Mfnbmu/Mdc3dEREQOVoVkByAiImWDEoqIiMSFEoqIiMSFEoqIiMSFEoqIiMSFEoqIiMSFEkoSmdmzZva3BK17hZmdkoh1lwfxfP0Kep/N7DIz+yIe2zlQZtbfzNKTGUMymNkIM3sxSdt2M+sYDo8xszuKabsLzax/otavhII+fMuS8vrhKKWXuw9193vivd5YX2Tc/TB3/yze28qjhCJlhplVKgnrkPgrr+9Laau3EkoBzKyqmT1iZmvCv0fMrGo4r6GZvWdm28xsi5lNNbMK4bxbzGy1me00syVmdnIRtlUvXN9GM9saDreMmP+Zmd1jZl+G6/3QzBpGzL/EzFaa2WYz+2vUuvuYWaqZ7TCz9Wb2cMS848zsq7Aeq8zssnD66WY2O1xmlZmNiFimbXjKfnX4uqw1sz9FzK9gZrea2bIwnnFmVj+fevc3s3Qz+4uZbQrPFi+Keg8eNLPvw9jHmFn1qGVvMbN1wCvARKC5me0K/5pHf1OLPosJt3mLmc0DdkccxEeZ2aLw/XjGzKpFLHOGmc0JX7evzOyIiHm9zGxW+D69BvywXD7MzEab2XYzWxy5v5jZ5Wb2TbiuNDP7fcS8gvbB5mb2Rrg/LTez6yKWqx6+JlvNbBFwVCHBHWNmX4fxfW1mx0TMK3C/jFpPgft4jPI/e18i9qud4Xvz64jyl5nZF+H+sjWs98CI+e3MbEq47EdAw6jtnWnBJaFtYb0OjYrlZjObZ2a7zewpM2tiZhPD9X1sZvUKqMvNFhwna8zsiqh5P+yfMfbpZ6yQ48liHMNmdjVwEfDn8Dh4N6Iep4TDBX2+5cUx3Mw2hLFfnl/9fuDu5f4PWAGcEmP6SGA60BhoBHwF3BPO+ycwBqgc/h0PGNAFWAU0D8u1BTrks91ngb+Fww2A3wA1gFrA68DbEWU/A5YBnYHq4fi94bxuwC7gBKAq8DCQnVcnYBpwSTh8CNAvHG4D7AQuCOvQAOgZzusPHE7wpeMIYD1wdkSdnOADvGZYbmPE9q4PX7eWYTxPAK/k8xr0D2N9OCx7IrAb6BLOHwVMAOqHr8u7wD+jlr0vXLZ6OC09v9c5Yrn0iPEVwBygFVA9YtqCcFp94MuI96oXsAHoC1QEhoTlqwJVgJXAjeFrOhjIitx+VGyXhXXIK38esB2oH84/HehAsG+dCOwBjixkH6wAzATuDONpD6QBp4XL3QtMDevVKqxnej7x1Qe2ApcAlQj2la1Ag8L2yxjrKnAfz+e4jH5ffgs0D+t4HsG+0izitcwCrgrfl2uANYBFHAd5+9kJBPv+i+G8zuG6Tg1fyz8DS4EqEbFMB5oALcL3f1a4L1QDJgN35VOPAQTHT3eC4+VlguOnY4zPgf78fJ/O93ii4GP4h/XG+qyj4M+3vDhGhusdRLDv1SvwszTZH+Yl4Y/8E8oyYFDE+GnAiog34528nSKiTMdwZzsFqFzIdn/2hkfM6wlsjRj/DLg9YvwPwAfh8J3AqxHzagKZETvO58DdQMOobdwGvFXE1+gRYFQ43DY8ILpGzL8feCoc/gY4OWJeM4IDvVKM9ebtuDUjpo0D7iD4cNxNREIGjgaWRyybCVSLWt+BJJQrYuwTQyPGBwHLwuHH8w68iPlLCD7wTyDiQyyc91UB7/NlMcr/j/ALQIzybwPXF7IP9gW+j/FePxMOpwEDIuZdHf2aRcy7BPhf1LRpwGWF7ZdF2Kd6ErGPx5j/s/clRpk5wFkRr+XSiHk1wv20KdA6xn72Mj8mlDuAcRHzKgCrgf4RsVwUMf8N4PGI8T+ST3IEniYiyRIkr4ISSvQ+ne/xRAHHMIUnlII+3/oDe4k4Zgk+1/oV9H7oklfBmhN828yzMpwG8ADBN5gPw0sRtwK4+1LgBmAEsMHMXjWz5hTCzGqY2RMWXLbaQZAE6ppZxYhi6yKG9xCcbeTFuSpvhrvvBjZHlL2SYCdeHF6yOCOc3opgp4oVT18z+zS8PLEdGErUJYLIbfLT16YN8FZ4Cr6N4IDIIfh2F8vWMObodTUi+FCYGbGuD8LpeTa6e0Y+690fqwqZFl2/4XkxhXG1Cuc3B1Z7eARGLFuQWOWbA5jZQDObbsElrW0EiS3vfYi5D4bxNY+K7y/8+Pr/ZH8pJL7oYyCvfIuI8fz2y58o4j4e7Sfvi5ldaj9eatxG8K0/cr/8IRZ33xMOHhLWI9Z+lucn9XT33HDbkfVcHzG8N8Z4zHqzf683/HyfLuh4yvcYLoKCPt8ANrt7dsR4vu9tHiWUgq0heDPztA6n4e473X24u7cHzgRusvDat7u/7O7Hhcs6welrYYYTXC7r6+61Cb7pQvAtvTBrCXasYAGzGgSnvoTxfOfuFxCc2t4HjDezmgQ7eYd81vkywaWmVu5eh+DSSnQsrSKGf3htwvUOdPe6EX/V3H11PtuqF8YTva5NBAfqYRHrqePukTt15AdxrHEIznJqRIw3jVEm1nIF1e/vUfWr4e6vELwXLczMopYtSKzya8Lr2W8ADwJN3L0u8D7h+1DAPriK4CwuMr5a7j4oXP9P9pdC4os+BvLK5/deFuRA9vEf3hczawM8CQwjuORWl+ByXVGPkVj7WZ6f1DN8P1pxYPWMte2ivt7w832xoOOpoGM41j4dKd/PtwOlhPKjymZWLeKvEkEbwe1m1siChsY7gRfhh0bZjuGOt53gG0OumXUxs5PCD4MMgg/E3CJsv1ZYdlvY4HbXfsQ+HjgjbJyrQnAp5If31swuNrNG4beubeHkXOAl4BQzO9eCBs8GZtYzIp4t7p5hZn2AC2Ns947wW+dhwOXAa+H0McDfww8AwtfvrELqcLeZVTGz44EzgNfDeJ8ERplZ43BdLczstALWsx5oYGZ1IqbNAQaZWX0za0pwBlkU15pZy/D9+GtE/Z4EhoZncWZmNS34EUMtgstB2cB1ZlbZzM4B+hSyncYR5X8LHEqQOKoQXDPfCGRb0MD8y7yF8tsHCS6Z7bSgYbe6mVU0s+5mltf4Pg64zYJG8pYEl2vy8z7Q2cwuDPeR8wja7N4rwusX7WD2cQgu5TrB64EFjcTdi7Kgu68EUvlxPzsO+FVEkXHA6WZ2splVJkh++wguVx6sccBlZtYt/LK3v/Uu6Hgq6BheT9B+lp98P98OlBLKj94n2Nnz/kYAfyPYCecB8wka4fJ+LdQJ+JigMXwa8G93/5TgA+Begm/X6wg+LG4rwvYfIWiA20TQUPZBUQN394XAtQRnFWsJGk0j78UYACw0s13Av4Dz3X2vu39PcAllOLCF4IO3R7jMH4CRZraTYEcbF2PTUwguuXwCPOjuH4bT/0VwdvNhuPx0guv6+VkXxryG4AAZ6u6Lw3m3hNuYHl4m+ZjgW25+r8ViggMlLbxE0Bx4AZhLcP34Q35MDIV5OSyfRnBZ4W/hNlIJGn5Hh3EvJbh+j7tnAueE41sIGo7fLGQ7Mwj2p03A34HB7r7Z3XcC1xG89lsJkvqEiOVi7oPunkOQlHsCy8P1/gfIS7J3E1zeWB7W74X8AnP3zeG6hhNcRv0zcIa7byqkTrE8wgHu42Esi4CHCOq6nuDHIF/uxyouJNgPtxB8qD8fse4lwMXA/4Xx/Qr4Vfh+HhR3n0hQ98kE+8rk/VxFvsdTIcfwU0C38Dh4O8Z6C/p8OyB5v34QKTIza0vwYVQ56hrrgayrP0HDaL4/HxWR0kFnKCIiEhdKKCIiEhe65CUiInGhMxQREYmLUtXxWLw1bNjQ27Ztm+wwRERKlZkzZ25y90bR08t1Qmnbti2pqanJDkNEpFQxs5h3++uSl4iIxIUSioiIxIUSioiIxIUSioiIxIUSioiIxIUSioiIxEVCE4qZDbDgmepL7ceH/0TOr2pmr4XzZ4SdDuY9s3yvBQ/SmWNmYyKW6W1m88NlHs17joSZjbDgOe55ywyK3p6IiCROwhKKBU9hewwYSPD8hAvMrFtUsSsJnqLWkeDZ4ZEPolrm7j3Dv6ER0x8n6Dq8U/g3IGLeqIhl3o9zlURESr3te7MYMWEhOzOy4r7uRJ6h9CF4vnNa+EyBV4HohyydBTwXDo8HTs4744jFzJoBtd19evjI1OeBs+MeuYhIGTT7+62c/uhUXpy+khlpW+K+/kQmlBb89DnK6fz0+cw/KRM+V2M7Pz66tp2ZzTazKeFT/PLKRz44Knqdw8xsnpk9bWb1YgVlZlebWaqZpW7cuPGAKiYiUprk5jpjP1/Gb8dMwx3GDT2aU7o1ift2Smqj/Fqgtbv3Am4CXjaz2oUs8zjBs5V7hss/FKuQu4919xR3T2nU6Gdd0YiIlCmbd+3jiue+5h/vL+aUQ5vw/nXHc2TrmN+3D1oi+/JaDbSKGG8ZTotVJt2CZ7jXATaHl7P2Abj7TDNbBnQOy0c+2e+Hdbr7+ryJZvYkB/bMaxGRMmPass3c8Npstu7O4p6zDuPifm0ooFXhoCXyDOVroJOZtTOzKsD5/PR52ITjQ8LhwcBkd3czaxQ26mNm7Qka39PcfS2ww8z6hW0tlwLvhOWaRaz318CCRFVMRKQky8l1Rn30LRf9Zzo1q1TirWuP4ZKj2yY0mUACz1DcPdvMhgGTgIrA0+6+0MxGAqnuPgF4CnjBzJYCWwiSDsAJwEgzywJygaHunteC9AfgWaA6MDH8A7jfzHoCDqwAfp+ouomIlFTrtmdw/auzmbF8C+f0asE9Z3enZtXi6Vi+XD+xMSUlxdV9vYiUFZ8u2cDwcXPZm5nDPWd3Z3DvloUvdADMbKa7p0RPL9fPQxERKQsys3N5YNJinpy6nK5NazH6wiPp2PiQYo9DCUVEpBRbtWUPw16ZzdxV27i4X2tuP70b1SpXTEosSigiIqXU+/PXcssb88Dh3xcdyaDDmxW+UAIpoYiIlDIZWTnc894iXprxPT1a1WX0Bb1oVb9GssNSQhERKU2WbtjFsJdnsXjdTq4+oT1/+mUXqlQqGfeoK6GIiJQSb85K5/a3F1C1UgWeuewoftG1cbJD+gklFBGREm5PZjZ3vrOQ8TPT6dOuPo+e34umdaolO6yfUUIRESnBFq/bwbUvzSJt026uO7kT153UkUoVS8YlrmhKKCIiJZC78+rXqxgxYSG1q1fmxSv7cmzHhskOq0BKKCIiJczOjCz+8tYC3p27huM7NeThc3vSqFbVZIdVKCUUEZESZH76doa9Mov0rXu5+bQuXHNiBypUSGynjvGihCIiUgK4O89+tYJ/vP8NDQ+pyqtX9+OotvWTHdZ+UUIREUmybXsy+fP4eXy4aD0nd23Mg7/tQb2aVZId1n5TQhERSaKZK7dy3Suz2bAzg9tPP5Qrj2uX8OeWJIoSiohIEuTmOmOnpvHApCU0r1uN14ceQ89WdZMd1kFRQhERKWZbdmcyfNwcPl2ykUGHN+Wf5xxBneqVkx3WQVNCEREpRjNXbmHYy7PZvCuzWJ7zXpyUUEREioG785+py7nvg8U0r1udN645hsNb1kl2WHGlhCIikmDb92Qx/PW5fPzNek47rAn3D+5RJi5xRVNCERFJoLmrtnHty7NYvyODO8/oxuXHti0zl7iiKaGIiCRA5I2KjWuVjV9xFUYJRUQkznZkZHHL+HlMXLCOUw4NblSsW6P03ai4v5RQRETiaMHq7Vz7ctAX118GdeWq49uX2Utc0ZRQRETiwN15acb3jHx3EfVrVuG1q/uRUsr64jpYSigiIgdp175sbntzPu/OXcMJnRsx6tweNDik5Hc3H29KKCIiB2Hxuh384cVZrNi8u9R1Nx9vSigiIgdo3NeruOOdBdSuXpmXftePozs0SHZISaWEIiKyn/Zm5nDHOwsYPzOdYzs24JHzepWKJyommhKKiMh+WLZxF9e+NIsl63dy3cmduP7kTlQsp5e4oimhiIgU0Xvz1nDL+HlUrVyR5y7vwwmdGyU7pBJFCUVEpBD7snP4x3+/4blpK+ndph6jL+xFszrVkx1WiaOEIiJSgPSte7j2pVnMTd/O745rxy0Du1K5YoVkh1UiKaGIiORj8uL13PjaXHJznTEXH8mA7s2SHVKJpoQiIhIlOyeXhz76lsc/W0a3ZrX590VH0rZhzWSHVeIpoYiIRNiwI4M/vjKbGcu3cEGfVtz1q8OoVrlissMqFZRQRERCXy3bxHWvzGH3vmwePrcH5xzZMtkhlSpKKCJS7uXmOv/+bCkPf/Qt7RrW5OWr+tK5Sa1kh1XqKKGISLm2dXcmN46bw2dLNnJmj+b885zDqVlVH40HQq+aiJRbqSu2cN0rs9m0K5N7zu7OxX1bl5tnlyRCQn9MbWYDzGyJmS01s1tjzK9qZq+F82eYWdtwelsz22tmc8K/MRHL9Daz+eEyj1r47ptZfTP7yMy+C//XS2TdRKT0ysl1Rk/+jvPGTqdSxQqMv+ZoLunXRsnkICUsoZhZReAxYCDQDbjAzLpFFbsS2OruHYFRwH0R85a5e8/wb2jE9MeBq4BO4d+AcPqtwCfu3gn4JBwXEfmJ9TsyuPg/M3jww28ZdHgz/nvdcRzRsm6ywyoTEnmG0gdY6u5p7p4JvAqcFVXmLOC5cHg8cLIV8BXBzJoBtd19urs78Dxwdox1PRcxXUQECG5UHPivqcxZtY37f3MEj57fk1rVKic7rDIjkW0oLYBVEePpQN/8yrh7tpltB/IeKNDOzGYDO4Db3X1qWD49ap0twuEm7r42HF4HNIlXRUSkdMvMzuX+Dxbzny+W07VpLUZf2IuOjfUrrngrqY3ya4HW7r7ZzHoDb5vZYUVd2N3dzDzWPDO7GrgaoHXr1nEJVkRKrhWbdvPHV2Yzf/V2Lj26DX8ZdKhuVEyQRCaU1UCriPGW4bRYZdLNrBJQB9gcXs7aB+DuM81sGdA5LB95p1HkOtebWTN3XxteGtsQKyh3HwuMBUhJSYmZdESkbHh79mr++tZ8KlWswJiLezOge9Nkh1SmJbIN5Wugk5m1M7MqwPnAhKgyE4Ah4fBgYHJ4dtEobNTHzNoTNL6nhZe0dphZv7Ct5VLgnRjrGhIxXUTKmd37shk+bi43vDaHbs1r8/71xyuZFIOEnaGEbSLDgElAReBpd19oZiOBVHefADwFvGBmS4EtBEkH4ARgpJllAbnAUHffEs77A/AsUB2YGP4B3AuMM7MrgZXAuYmqm4iUXAvXbOePL89m+ebdXHdSR647uROV1N18sbDg6lL5lJKS4qmpqckOQ0TiwN157qsV/OP9xdSrWZlR5/XkmA4Nkx1WmWRmM909JXp6SW2UFxEpsq27M7l5/Dw+/mY9J3VtzAODj6DBIVWTHVa5o4QiIqXal0s3MXzcXDbv3scdZ3TjimPb6o73JFFCEZFSaU9mNvdOXMzz01bSvmFNnrz0WA5vWSfZYZVrSigiUuqkrtjC8NfnsnLzHq44th03n9aF6lV0b0myKaGISKmRkZXDqI++ZezUNFrUrc4rV/Xj6A4NCl9QioUSioiUCvPStzF83Fy+27CLC/u25i+DDuUQPbekRNG7ISIlWmZ2LqM/Xcpjny6l0SFVefbyo+jfpXGyw5IYlFBEpMRavG4Hw8fNZeGaHZzTqwV3/eow6tRQ78AllRKKiJQ4ObnO2M/TGPXRt9SqVkn9cJUSSigiUqKkbdzF8NfnMvv7bQzs3pS/nd1dNymWEkooIlIi5OY6z361gvsnLaZqpYr86/yenNmjuW5SLEWUUEQk6VZt2cPN4+cyPW0Lv+jSiHt/cwRNaldLdliyn5RQRCRp3J1xqasY+e4izIz7fnM456a00llJKaWEIiJJsWFHBre+OZ/JizfQr319Hhjcg1b1ayQ7LDkISigiUuzem7eG299ewN7MHO48oxuXHdOWChV0VlLaKaGISLHZujuTOycs5N25a+jRqi4P/bYHHRsfkuywJE6UUESkWHy6eAN/fmMeW3dnMvzUzlzTv4OepFjGKKGISELt2pfN3/+7iFf+t4ouTWrxzGVH0b2Fupkvi5RQRCRhpqdt5k+vz2X1tr38/sT23HRqZ6pWUjfzZZUSiojEXUZWDg9OWsJTXy6ndf0avP77o0lpWz/ZYUmCKaGISFzNS9/GTePmsnTDLi7u15rbBh5KTXUzXy7oXRaRuMjKyWX05KWMDruZf/6KPpzQuVGyw5JipIQiIgctbeMubnhtDvPSt/PrXi0YoW7myyUlFBE5YHldp4yYsIgqlSrw74uOZNDhzZIdliSJEoqIHJBtezK57c35TFywjmM6NOChc3vQrE71ZIclSaSEIiL77atlm7jptbls2rWPWwd25erj26vrFFFCEZGiy8zO5eGPvuWJz5fRrkFN3vrDsRzeUjcpSkAJRUSKJG3jLq5/dQ7zV2/ngj6tuOOMbtSooo8Q+ZH2BhEpUGTDe9XKFRhz8ZEM6K6Gd/k5JRQRyVd0w/vD5/akaR09SVFiU0IRkZjyGt43797HbQO7cpUa3qUQSigi8hPRDe9PXqqGdykaJRQR+YEa3uVgaE8REdyd8TPTuWvCQqpUUsO7HBglFJFybkdGFn99awHvzl1Dv/b1GXVeT93xLgdECUWkHJv1/Vaue2U2a7dncPNpXRh6YgcqquFdDlCREoqZ1QT2unuumXUGugIT3T0rodGJSELk5Dpjpizj4Y++pVmdaoz7/dH0blMv2WFJKVfUM5TPgePNrB7wIfA1cB5wUaICE5HEWLc9gxtfm8O0tM2ccUQz/nHO4dSupq7m5eAVNaGYu+8xsyuBf7v7/WY2J4FxiUgCfLxoPTePn0tGVi73Dz6C3/ZuiZkucUl8FDmhmNnRBGckV4bTKiYmJBGJt4ysHO6duJhnv1pBt2a1+b8Le9Gh0SHJDkvKmApFLHcDcBvwlrsvNLP2wKeFLWRmA8xsiZktNbNbY8yvamavhfNnmFnbqPmtzWyXmf0pYtr1ZrbAzBaa2Q0R00eY2WozmxP+DSpi3UTKtO/W7+Tsx77k2a9WcMWx7Xjr2mOUTCQhinSG4u5TgCkAZlYB2OTu1xW0jJlVBB4DTgXSga/NbIK7L4oodiWw1d07mtn5wH0EbTN5HgYmRqyzO3AV0AfIBD4ws/fcfWlYZJS7P1iUOomUde7Oq1+v4u53F1KzSiWeuewoftG1cbLDkjKsSGcoZvaymdUOf+21AFhkZjcXslgfYKm7p7l7JvAqcFZUmbOA58Lh8cDJFl7QNbOzgeXAwojyhwIz3H2Pu2cTJLlzilIHkfJk+54srn15Fre9OZ+UNvWZeP3xSiaScEW95NXN3XcAZxOcMbQDLilkmRbAqojx9HBazDJhgtgONDCzQ4BbgLujyi8g+LVZAzOrAQwCWkXMH2Zm88zs6fAXaT9jZlebWaqZpW7cuLGQKoiUPqkrtjDo0al8uHA9tw3syvNX9KFxbfUQLIlX1IRS2cwqEySUCeH9J56wqGAEweWrXZET3f0bgstiHwIfAHOAnHD240AHoCewFngo1ordfay7p7h7SqNGjRIRu0hSZOfk8sjH33LuE9OoVNF445pj+P2JHdRDsBSbov7K6wlgBTAX+NzM2gA7CllmNT89e2gZTotVJt3MKgF1gM1AX2Cwmd0P1AVyzSzD3Ue7+1PAUwBm9g+CMx/cfX3eSs3sSeC9ItZNpNRbvW0vN7w6m69XbOXXvVow8qzDqKV7S6SYFbVR/lHg0YhJK83sF4Us9jXQyczaESSO84ELo8pMAIYA04DBwGR3d+D4vAJmNgLY5e6jw/HG7r7BzFoTtJ/0C6c3c/e14WK/Jrg8JlLmvT9/Lbe+MY+cXGfUeT34da+WyQ5Jyqmidr1SB7gLOCGcNAUYSdDmEZO7Z5vZMGASwT0rT4c/OR4JpLr7BIIzjRfMbCmwhSDpFOYNM2sAZAHXuvu2cPr9ZtaT4FLcCuD3RambSGm1JzObe95bxCv/W0WPVnV59PyetGlQM9lhSTlmwQlBIYXM3iD4xp/3i6xLgB7uXqp/YZWSkuKpqanJDkNkvy1as4M/vjKLtE27GXpiB246tTOVKxa1SVTk4JjZTHdPiZ5e1DaUDu7+m4jxu9X1ikjxc3ee+XIF905cTN0alXnxyr4c27FhssMSAYqeUPaa2XHu/gWAmR0L7E1cWCISbdOufdz8+lw+XbKRUw5tzP2De1C/ZpVkhyXyg6ImlKHA82FbCsBWgsZ0ESkGU7/byE3j5rJ9bxZ3n3kYlx7dRp06SolT1F95zQV6mFntcHxH2I/WvATGJlLuZWbn8tCHS3ji8zQ6Nj6E56/ow6HNaic7LJGY9uuJjeHd8nluAh6JazQi8oPlm3Zz/auzmZe+nQv7tuaO07tRvYo6+ZaS62AeAazzbZEEcHfemLWaO99ZQOWKFRhz8ZEM6N4s2WGJFOpgEkoiu14RKZf2ZuZw+9sLeGNWOn3a1eeR83rSvG71ZIclUiQFJhQz20nsxGGA9nKROFq+aTfXvDiTJet3ct3Jnbj+5E5UVD9cUooUmFDcvVZxBSJSnn2wYB03vz6XihWNZy47iv5d1NW8lD4Hc8lLRA5Sdk4u909awtjP0+jRsg6PXXQkLevVSHZYIgdECUUkSTbsyGDYK7P53/ItXNyvNXec0Y2qlfQrLim9lFBEkmB62maGvTybXfuy1EOwlBlKKCLFyN0Z+3ka909aQpv6NXjpd33p0lRNlVI2KKGIFJMdGVnc/PpcJi1cz8DuTbl/8BF6CJaUKUooIsXgm7U7uObFmazaupfbTz+UK49rp764pMxRQhFJsPEz07n97fnUrlaZV6/ux1Ft6yc7JJGEUEIRSZCMrBzufncRr/zve45u34BHL+hFo1pVkx2WSMIooYgkwKote7jmpZksWL2Da/p3YPipnamkJypKGaeEIhJH7s5rX6/invcWUaGC8eSlKZzarUmywxIpFkooInGyYUcGt745n8mLN3B0+wbcP/gIWtXXXe9SfiihiMTBhLlruOPtBWRk5XDXr7ox5Oi2VFDHjlLOKKGIHIStuzO5/Z0F/HfeWnq2qstD5/agQ6NDkh2WSFIooYgcoMmL13PLG/PZtieTm0/rwu9PaK+GdynXlFBE9tPOjCzueW8R41LT6dq0Fs9d3oduzfWcdxElFJH98NWyTdz8+jzWbt/LH/p34PpTOqmHYJGQEopIEWRk5XDfB4t55ssVtGtYk9eHHkPvNvWSHZZIiaKEIlKI2d9vZfjrc0nbuJshR7fhloFdqVFFh45INB0VIvnIzM7l0U++49+fLaVp7Wq89Lu+HNuxYbLDEimxlFBEYpj1/Vb+8uZ8Fq/byeDeLbnzV92ora7mRQqkhCISYUdGFg98sIQXZ6ykSa1q6jpFZD8ooYgQ9ME1ccE6RkxYyKZd+7jsmLYM/2UXDqmqQ0SkqHS0SLmXvnUPd76zkMmLN3BY89r8Z0gKR7Ssm+ywREodJRQpt7JzcnnmyxU8/NG3mMHtpx/KZce01d3uIgdICUXKpTmrtvGXN+ezaO0OTjm0MXef1Z0WdasnOyyRUk0JRcqVnRlZPPThtzw3bQWNa1VlzMVHctphTfV8d5E4UEKRcsHdmbRwHXdNWMiGnfsYcnRbhv+yM7X0U2CRuFFCkTJv9ba93PXOAj7+ZgPdmtXmiUtS6NmqbrLDEilzlFCkzMrKyeW5r4JGd3f466BDufxYNbqLJIoSipRJ05ZtZsSEhSxZv5OTujZm5FmH0bKeHscrkkgJ/apmZgPMbImZLTWzW2PMr2pmr4XzZ5hZ26j5rc1sl5n9KWLa9Wa2wMwWmtkNEdPrm9lHZvZd+F9dwZZDa7fvZdjLs7jgyensycrmyUtTeGpIipKJSDFIWEIxs4rAY8BAoBtwgZl1iyp2JbDV3TsCo4D7ouY/DEyMWGd34CqgD9ADOMPMOoazbwU+cfdOwCfhuJQT+7Jz+PdnSznpwSl8tGg9N57SmY9uPJFTuzXRL7hEikkiL3n1AZa6exqAmb0KnAUsiihzFjAiHB4PjDYzc3c3s7OB5cDuiPKHAjPcfU+4zinAOcD94br6h+WeAz4Dbol3paTk+WzJBu5+dxHLN+3mtMOacPvp3WhVX2ckIsUtkQmlBbAqYjwd6JtfGXfPNrPtQAMzyyBIBqcCf4oovwD4u5k1APYCg4DUcF4Td18bDq8DYvboZ2ZXA1cDtG7d+sBqJiXCqi17GPneIj5atJ72DWvy3BV9OLFzo2SHJVJuldRG+RHAKHffFXm5wt2/MbP7gA8JzlzmADnRC4dnOB5rxe4+FhgLkJKSErOMlGx7M3N4fMoyxkxZRqUKxq0Du3LFse2oUkm/3hJJpkQmlNVAq4jxluG0WGXSzawSUAfYTHAmM9jM7gfqArlmluHuo939KeApADP7B8GZD8B6M2vm7mvNrBmwIUH1kiQJbk5czz3vLWL1tr2c2aM5fxl0KE3rVEt2aCJCYhPK10AnM2tHkDjOBy6MKjMBGAJMAwYDk93dgePzCpjZCGCXu48Oxxu7+wYza03QftIval33hv/fSVC9JAmWbdzFiAkLmfrdJro0qcWrV/ejX/sGyQ5LRCIkLKGEbSLDgElAReBpd19oZiOBVHefQHCm8YKZLQW2ECSdwrwRtqFkAde6+7Zw+r3AODO7ElgJnBvfGkky7MjI4rFPl/L0F8upVrkiI37VjYv7tdHNiSIlkAUnBOVTSkqKp6amFl5Qil1GVg7PfbWCx6csY9ueLH7buyW3DOxKw0OqJjs0kXLPzGa6e0r09JLaKC/lVFZOLuNSV/HoJ9+xfsc+TuzciJtP60L3FnWSHZqIFEIJRUqE3Fzn3XlrGPXRt6zYvIfeberx6Pm96Kt2EpFSQwlFksrd+XTJBh6Y9C3frN1B16a1eGpICid1baw73EVKGSUUSZr/Ld/CA5MW8/WKrbSuX4NHzuvJmT2aU6GCEolIaaSEIsVu4ZrtPDhpCZ8u2UjjWlX529ndOe+oVlTWL7dESjUlFCk2yzft5uGPvuXduWuoU70ytwzoymXHtKV6lYrJDk1E4kAJRRJu/Y4MHvn4O8alrqJKxQpc+4sOXH1CB+pU1+N3RcoSJRRJmL2ZOTw5NY3HP1tGdm4uF/dtzbUndaRxLXWVIlIWKaFI3Lk7E+au4b6Ji1mzPYMBhzXltkFdadOgZrJDE5EEUkKRuJr1/VZGvruIOau2cVjz2jx8Xk/1uSVSTiihSFys3raX+yYuZsLcNTSqVZUHBh/Bb45sqZ8Ai5QjSihyUHbvy2bMlGWM/TwNgD+e1JGhJ3agZlXtWiLljY56OSC5uc74Wek8OGkJG3bu46yezfnzgK60qFs92aGJSJIooch+m5G2mXv+u4gFq3fQq3VdxlzSmyNb10t2WCKSZEooUmTfb97DPyd+w8QF62hepxr/Oj/oKkV9bokIKKFIEWzcuY8npizj+WkrqVTRGH5qZ646oT3VKusOdxH5kRKK5Gvjzn2M/XwZL0xfSWZ2LoN7t2T4L7vQpLZuTBSRn1NCkZ+JTiS/7tWSYSd1pF1D3ZgoIvlTQpEfKJGIyMFQQpGfJZKze7Xgjyd1UiIRkf2ihFKObdq1j7Gfp/H8tBVKJCJy0JRQyiElEhFJBCWUciQvkbwwbSX7snOUSEQkrpRQyoGtuzMZOzWNZ79cESSSni0YdlJH2jc6JNmhiUgZooRShu3IyOKpqct56ovl7M7M5swezbn+5E5KJCKSEEooZdDufdk8+9UKxn6exva9WQw6vCk3nNKZzk1qJTs0ESnDlFDKkIysHF6cvpLHP1vG5t2ZnNy1MTee2pnuLeokOzQRKQeUUMqAfdk5vPb1KkZPXsqGnfs4vlNDbjq1M73UA7CIFCMllFIsKyeXN2el8+gnS1m9bS992tbn/y7oRV89cldEkkAJpRTKyXUmzF3NIx9/x8rNe+jRqi73/uZwjuvYUF3Ji0jSKKGUIrm5zsQF6xj18bcs3bCLbs1q89SQFE7q2liJRESSTgmlFMjOyeW9eWsZ/elSlm7YRafGh/D4RUdy2mFNqVBBiURESgYllBIsMzuXt2an8+/PlrFy8x66NKnFoxf04vTDm1FRiUREShgllBIoIyuH11NXMWZKGqu37eXwFnV44pLenHpoE52RiEiJpYRSguzJzOblGd/zxOdpbNy5j95t6vH3X3fnxM6N1EYiIiWeEkoJsDMji+enreSpL5azZXcmx3RowKPn96Jf+/pKJCJSaiihJNG2PZk8/eUKnv1yOTsysvlFl0YMO6kjvdvUT3ZoIiL7TQklCTbu3Md/vkjjxWkr2Z2Zw4DDmjLspI7qIkVESjUllGK0assenpyaxrjUVWRm53LGEc259hcd6dJUnTaKSOmX0IRiZgOAfwEVgf+4+71R86sCzwO9gc3Aee6+ImJ+a2ARMMLdHwyn3Qj8DnBgPnC5u2eY2bPAicD2cPHL3H1Owiq3Hxat2cGYKcv47/y1VDA4u2cLrunfQd3Ii0iZkrCEYmYVgceAU4F04Gszm+DuiyKKXQlsdfeOZnY+cB9wXsT8h4GJEetsAVwHdHP3vWY2DjgfeDYscrO7j09UnfaHuzMtbTNjpqTx+bcbqVmlIlce144rjm1H0zrVkh2eiEjcJfIMpQ+w1N3TAMzsVeAsgjOOPGcBI8Lh8cBoMzN3dzM7G1gO7I4Rc3UzywJqAGsSVoMDkJPrfLhwHWOmLGNu+nYaHlKFm0/rwsX92lCneuVkhycikjCJTCgtgFUR4+lA3/zKuHu2mW0HGphZBnALwdnNn/IKu/tqM3sQ+B7YC3zo7h9GrO/vZnYn8Alwq7vviw7KzK4GrgZo3br1wdUwwr7sHN6ctZqxn6exfNNu2jSowd/O7s7g3i2pVrli3LYjIlJSldRG+RHAKHffFXkfhpnVIziraQdsA143s4vd/UXgNmAdUAUYS5CQRkav2N3HhvNJSUnxgw10R0YWL03/nqe/XM7Gnfvo3qI2oy/sxcDu6h5FRMqXRCaU1UCriPGW4bRYZdLNrBJQh6Bxvi8w2MzuB+oCueFZy3pgubtvBDCzN4FjgBfdfW24zn1m9gwRZzaJsGFHBk99uZyXp3/Pzn3ZHNexIaPO7cmxHRvoZkQRKZcSmVC+BjqZWTuCxHE+cGFUmQnAEGAaMBiY7O4OHJ9XwMxGALvcfbSZ9QX6mVkNgkteJwOpYblm7r7Wgk/zs4EFiarYo598x+jJS8nOzWXg4c0YekIHDm+pe0hEpHxLWEIJ20SGAZMIfjb8tLsvNLORQKq7TwCeAl4ws6XAFoKkU9A6Z5jZeGAWkA3MJrx8BbxkZo0AA+YAQxNQLQBa1K3Ob1NactXx7WnbsGaiNiMiUqpYcEJQPqWkpHhqamqywxARKVXMbKa7p0RPr5CMYEREpOxRQhERkbhQQhERkbhQQhERkbhQQhERkbhQQhERkbhQQhERkbhQQhERkbgo1zc2mtlGYGWy4zgADYFNyQ6iGJW3+oLqXF6U1jq3cfdG0RPLdUIprcwsNdZdqmVVeasvqM7lRVmrsy55iYhIXCihiIhIXCihlE5jCy9SppS3+oLqXF6UqTqrDUVEROJCZygiIhIXSigiIhIXSiilgJnVN7OPzOy78H+9AsrWNrN0MxtdnDHGU1Hqa2Y9zWyamS00s3lmdl4yYj1YZjbAzJaY2VIzuzXG/Kpm9lo4f4aZtU1CmHFVhDrfZGaLwvf1EzNrk4w446mwOkeU+42ZuZmVyp8SK6GUDrcCn7h7J+CTcDw/9wCfF0tUiVOU+u4BLnX3w4ABwCNmVrf4Qjx4ZlYReAwYCHQDLjCzblHFrgS2untHYBRwX/FGGV9FrPNsIMXdjwDGA/cXb5TxVcQ6Y2a1gOuBGcUbYfwooZQOZwHPhcPPAWfHKmRmvYEmwIfFE1bCFFpfd//W3b8Lh9cAG4Cf3blbwvUBlrp7mrtnAq8S1D1S5GsxHjjZzKwYY4y3Quvs7p+6+55wdDrQsphjjLeivM8QfBm8D8gozuDiSQmldGji7mvD4XUESeMnzKwC8BDwp+IMLEEKrW8kM+sDVAGWJTqwOGsBrIoYTw+nxSzj7tnAdqBBsUSXGEWpc6QrgYkJjSjxCq2zmR0JtHL3/xZnYPFWKdkBSMDMPgaaxpj118gRd3czi/Vb7z8A77t7emn4AhuH+uatpxnwAjDE3XPjG6Ukk5ldDKQAJyY7lkQKvww+DFyW5FAOmhJKCeHup+Q3z8zWm1kzd18bfoBuiFHsaOB4M/sDcAhQxcx2uXtB7S1JE4f6Yma1gf8Cf3X36QkKNZFWA60ixluG02KVSTezSkAdYHPxhJcQRakzZnYKwZeLE919XzHFliiF1bkW0B34LPwy2BSYYGZnuntqsUUZB7rkVTpMAIaEw0OAd6ILuPtF7t7a3dsSXPZ6vqQmkyIotL5mVgV4i6Ce44sxtnj6GuhkZu3C+pxPUPdIka/FYGCyl+67kQuts5n1Ap4AznT3mF8mSpkC6+zu2929obu3DY/f6QR1L1XJBJRQSot7gVPN7DvglHAcM0sxs/8kNbLEKEp9zwVOAC4zsznhX8+kRHuAwjaRYcAk4BtgnLsvNLORZnZmWOwpoIGZLQVuouBf+JV4RazzAwRn2a+H72t0ki1ViljnMkFdr4iISFzoDEVEROJCCUVEROJCCUVEROJCCUVEROJCCUVEROJCCUWkAGaWE/50dYGZvW5mNfZj2bZmduEBbneFmTXcj/L9zeyYiPGhZnbpgWxb5EApoYgUbK+793T37kAmMLQoC4V3tbcF9juhhL3T7q/+wA8Jxd3HuPvzB7AekQOmhCJSdFOBjmZW08yeNrP/mdlsMzsLwMwuM7MJZjaZoNv9ewm6w5ljZjeG8394To2ZvWdm/cPhXWb2kJnNJehGB+DPZjY/3E7HsNyvwueizDazj82sSfiMlKHAjeG2jjezEWb2p3CZnmY2PXy+yFsWPl/GzD4zs/vC9X9rZscXw2soZZgSikgRhGccA4H5BH1MTXb3PsAvgAfMrGZY9EhgsLufSHBX+9TwDGdUIZuoCcxw9x7u/kU4bbu7Hw6MBh4Jp30B9HP3XgTdoP/Z3VcAY4BR4bamRq37eeCW8Pki84G7IuZVCutxQ9R0kf2mziFFClbdzOaEw1MJukL5Cjgz7wwAqAa0Doc/cvctB7CdHOCNqGmvRPzPS0gtgdfCTjOrAMsLWqmZ1QHquvuUcNJzwOsRRd4M/88kuEQncsCUUEQKttfde0ZOCB9w9Rt3XxI1vS+wu4B1ZfPTqwLVIoYz3D0nqrzHGP4/4GF3nxBeLhtRSPyFyevJNwd9HshB0iUvkf03Cfhj3pMTw95xY9lJ0DV5nhVATzOrYGatCJ7kV5DzIv5PC4fr8GPX50MiykZvCwh6sgW2RrSPXAJMiS4nEg/6RiKy/+4haNOYFz4caTlwRoxy84CcsKH92XCZ5cAigl5nZxWynXpmNo/gLOKCcNoIgl54twKTgXbh9HeB8eEPBP4YtZ4hwJjwJ89pwOVFqaTI/lJvwyIiEhe65CUiInGhhCIiInGhhCIiInGhhCIiInGhhCIiInGhhCIiInGhhCIiInHx/4aLGKhgAodyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pyhessian.utils import normalization\n",
    "\n",
    "# generate random vector to do the loss plot\n",
    "\n",
    "v = [torch.randn_like(p) for p in model.parameters()]\n",
    "v = normalization(v)\n",
    "\n",
    "\n",
    "# used to perturb your model \n",
    "lams = np.linspace(-0.5, 0.5, 21).astype(np.float32)\n",
    "\n",
    "loss_list = []\n",
    "\n",
    "# create a copy of the model\n",
    "model_perb = ptcv_get_model(\"resnet20_cifar10\", pretrained=True)\n",
    "model_perb.eval()\n",
    "model_perb = model_perb.cuda()\n",
    "\n",
    "for lam in lams: \n",
    "    model_perb = get_params(model, model_perb, v, lam)\n",
    "    loss_list.append(criterion(model_perb(inputs), targets).item())\n",
    "\n",
    "plt.plot(lams, loss_list)\n",
    "plt.ylabel('Loss')\n",
    "plt.xlabel('Perturbation')\n",
    "plt.title('Loss landscape perturbed based on a random direction')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how different the loss landscape looks. In particular note that there is almost no change in the loss value (see the small scale of the y-axis). This is expected, since for a converged NN, many of the directions are typically degenarate (i.e. they are flat)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use gradient direction to perturb the model. While gradient is better than random vector, but it is not possible to use it to plot 3D loss landscape since you will need more than one direction. However, you can use top 2 Hessian vectors instead for that scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Loss landscape perturbed based on gradient direction')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyhUlEQVR4nO3deXwU9f3H8dcnF4FwBsIRCAn3fQcEkVOpeCHeJ4q1RVu11dbe/urZqq0Hba1Wq0K9z9rijSj3HZD7hgQIZ7hCuHJ+fn/MxK4xJ2Qy2d3P8/HII7s7s/P9zOzse2e/MzsjqooxxpjQE+F3AcYYY7xhAW+MMSHKAt4YY0KUBbwxxoQoC3hjjAlRFvDGGBOiLODLICJTReQRj6adISLneTHtcFCdy6+811lEJorIvOpo53SJyEgRyfSzhuoS+LqJyG9F5EWf6vjWMhWRtSIysgbabSsix0Qk0uu2itX6gLcwDB2hFFbmzKjqH1X1B2c6HRFJEREVkagzqKWHqs4601pKKpldqrpDVeuramF1t1WWWh/wJjScyRuwOqdhqkdNboX6KdjXuaANeBGpIyKTRWS3+zdZROq4w5qJyEcickREDonIXBGJcIf9SkR2iUiOiGwUkXMr0VYTd3pZInLYvd0mYPgsEXlYROa7050uIs0Chk8Qke0iclBEfldi2oNEJE1EjorIPhF5KmDYOSKywJ2PnSIy0X38IhH52n3OThF5IOA5xVs0k9zlskdE7g0YHiEivxaRrW4974hIfBnzPVJEMt2v0wfcLZIbSrwGT4jIDrf2f4hI3RLP/ZWI7AXeBD4FEt2vqcdEJLFkF0kpX58z3GmsAo4HvOEGisg69/WYIiKxAc+5WERWuMttgYj0DhjWT0SWu6/T28A3zyuDiMgzIpItIhsC1xcRuUVE1rvT2iYitwUMK28dTBSR9931KV1EfhLwvLruMjksIuuAgRUUd7aILHXrWyoiZwcMK3e9LGVav3TXl90i8gN3PeroDpsqIs+JyCcichwYVd566D6nvPX+ARF5LeD+4IB1faUEdJlUMB9z3P9H3HVqSCnzVe4ylW93HT0gIu+JyGsichSYKCKNROQld9nsEpFHJOADTkR+GLAerBOR/iLyKtAW+NCt65dS4tuGux5Mc9ePLSLywxLL5x0RecWd7loRSS3rtSuTqtbqPyADOK+Uxx8CFgHNgQRgAfCwO+xR4B9AtPs3DBCgC7ATSHTHSwE6lNHuVOAR93ZT4AqgHtAAeBf4T8C4s4CtQGegrnv/MXdYd+AYMByoAzwFFBTPE7AQmODerg8Mdm8nAznAde48NAX6usNGAr1wPqB7A/uA8QHzpDiBGueOlxXQ3k/d5dbGred54M0ylsFIt9an3HFHAMeBLu7wp4FpQLy7XD4EHi3x3Mfd59Z1H8ssazkHPC8z4H4GsAJIAuoGPLbGfSwemB/wWvUD9gNnAZHAze74dYAYYDtwj7tMrwTyA9svUdtEdx6Kx78GyAbi3eEXAR1w1q0RwAmgfwXrYASwDPi9W097YBtwvvu8x4C57nwlufOZWUZ98cBhYAIQhbOuHAaaVrReljKtscBeoAfOev4aznrUMeB1ygaGuvMQS/nrYUXr/QPAa+7t1sBB4EJ3WmPc+wmVeH+luHVGlZMh5S5TAjLGrSsfGO/WUhf4AOd9EoeTN0uA29zxrwJ24XxoCNARSC4tu0rWivPh9Ky7LPvivE9HB9Rxyl0mkTjr06Iq56ffAV5hgWUH/FbgwoD75wMZ7u2HgP8Wr5wB43TEefOfB0RX0O5Uyn7j9wUOB9yfBdwXcP/HwGfu7d8DbwUMiwPyAlaoOcCDQLMSbfwG+KCSy2gy8HSJlahrwPA/AS+5t9cD5wYMa+Wu0N95g/C/kI4LeOwd4P/clfk4AR+QwBAgPeC5eUBsiemdTsB/v5R14vaA+xcCW93bz+F+0AcM34gTwMOB3YAEDFtQzus8sZTxl+B+IJcy/n+An1awDp4F7CjltZ7i3t4GjA0YNqnkMgsYNgFYUuKxhcDEitbLUqb1Mu6Hc8B7pWTAv1KF9bCi9f4B/hfwvwJeLTGtz4GbK/H+SqHigC93mfLdgJ8TMKwFkIu7ceE+dh0wM6DOn5bR7jfTLVkrzgdNIdAgYPijwNSAOmYEDOsOnCxv+Zf2F7RdNEAiztZYse3uYwB/BrYA092vzr8GUNUtwN04C2+/iLwlIolUQETqicjz7tfNozih3Fi+3Q+5N+D2CZyt8eI6dxYPUNXjOFsnxW7F2TLZ4H7Fvth9PAnnQ6y0es4SkZnuV/xs4Hag5FfvnQG3A5dNMvCB+1X4CE7gF+KsyKU57NZccloJOFt6ywKm9Zn7eLEsVT1VxnSrYmcFj5Wcv58X1+TWleQOTwR2qfuOCXhueUobPxFARC4QkUXuV+wjOB80xa9DqeugW19iifp+y/+W/7fWlwrqK/keKB6/dcD9stbL0qYV2G5Fy7yi9bCi9T5QMnBViWVyDs7GR1XnozRVWaaUGDcZ5xvYnoDansfZkody3qeVqOmQquaUqKu81y5WqrhPIJgDfjfOwi/W1n0MVc1R1Z+rantgHPAzcftOVfUNVT3Hfa7idCFU5Oc43TtnqWpDnC1BcLZiK7IHZyVwniBSD6e7Bbeezap6Hc4K8zjwnojE4axkHcqY5hs4XSNJqtoIpyugZC1JAbe/WTbudC9Q1cYBf7GququMtpq49ZSc1gHgJNAjYDqNVDXwjRcYjKXdB+dbQL2A+y1LGae055U3f38oMX/1VPVNnNeitYhIieeWp7Txd4uzv+d94Amghao2Bj7BfR3KWQd34nzLCayvgape6E7/W+tLBfWVfA8Uj1/Wa1mePTjddsWSShmn5OtQ3npY7npfwk6cLfjAZRKnqo9Vou7S1o2SqrJMS05zJ84WfLOA2hqqao+A4WW9T8urbTcQLyINStR1Oq9dmYIl4KNFJDbgLwqnj/k+EUlwd7j8HqffsHgnW0f3jZmNs4VaJCJdRGS0++Y8hRNQRZVov4E77hFxdkjeX4Xa3wMuFmeHaQzOV/dvlruI3CgiCapaBBxxHy4CXgfOE5GrRSRKRJqKSN+Aeg6p6ikRGQRcX0q7/+d+8+gB3AK87T7+D+APIpLstp8gIpdWMA8PikiMiAwDLgbedev9J/C0iDR3p9VaRM4vZzr7gKYi0ijgsRXAhSISLyItcb5hVcYdItLGfT1+FzB//wRud7cuRUTixNkZ2ACn+6IA+ImIRIvI5cCgCtppHjD+VUA3nCCPwelbzgIKROQC4HvFTyprHcTp4skRZ8dxXRGJFJGeIlK84+8d4Dfi7NhvA9xVTm2fAJ1F5Hp3HbkG56v8R5VYfiW9A9wiIt3cMP6/SjynvPWw3PW+hNeAS0TkfHd5xIqzs71NGeMHysJZru3LGacqy/RbVHUPMB14UkQainOQQgcRGeGO8iJwr4gMcNe3jsXvLZz1vdS6VHUnTvfgo+789sb5Nv9aaeOfrmAJ+E9wArb47wHgESANWAWsBpa7jwF0Ambg7ORZCDyrqjNx3pCP4Wx97sV58/6mEu1PxtnZcgBnB+VnlS1cVdcCd+Bs7ezB2QkWeCz4WGCtiBwD/gJcq6onVXUHzlf+nwOHcIKwj/ucHwMPiUgOzgfbO6U0PRuni+BL4AlVne4+/hecra7p7vMX4fQLl2WvW/NunA+d21V1gzvsV24bi8TpupqB802nrGWxAeeDeZv7dTcReBVYidNfOZ3/BXVF3nDH34bzFfkRt4004IfAM27dW3D60lHVPOBy9/4hnJ2m/66gncU469MB4A/Alap60P1q/ROcZX8YJ9ymBTyv1HVQnWOgL8bZj5PuTvdFoPhD70Gcr+rp7vy9WlZhqnrQndbPcbo/fglcrKoHKpin0qb1KfBXYCbua+oOyi3naWWuh5VY7wPb3glcitNVlYWzVfwLKpFPqnoC53WZ765Tg0sZrdLLtAw34Xygr3Pn4z3c7iNVfddt/w2cgyL+g7MzF5w+9fvcuu7lu67D6ZffjbMj935VnVHF2sol3+5eNMFORFJwVuRoVS04w2mNxNkRVpktKRNCRKQbztEmdc50PTL+CZYteGOMx0TkMnF+29AEZ3/Qhxbuwc0C3hhT7Dacw4i34uwz+JG/5ZgzZV00xhgTomwL3hhjQlStOpFOs2bNNCUlxe8yjDEmaCxbtuyAqiaUNsyzgBeRLnz7kLf2wO9VdXJZz0lJSSEtLc2rkowxJuSISJm/zPUs4FV1I86xvsWnFt2Fc6ynMcaYGlBTffDn4pwMqqJzQBhjjKkmNRXw1+L8gvE7xDlveZqIpGVlZdVQOcYYE/o8D3j3PBTjcM6h/h2q+oKqpqpqakJCqfsJjDHGnIaa2IK/AFiuqvtqoC1jjDGumgj46yije8YYY4x3PA149zziY6j4jH3GGGOqmacBr6rHVbWpqmZ72AbPfLWZtbs9a8IYY4JS0J+q4MiJfN5YvIMJLy1h876cip9gjDFhIugDvklcDK//cDCREcINLy4m48Dxip9kjDFhIOgDHqBdszje+MFZFBQpN7y4mMzDJ/wuyRhjfBcSAQ/QqUUDXvn+IHJO5XPDi4vZd/SU3yUZY4yvQibgAXq2bsTU7w/iQE4uN7y4mIPHyrucpDHGhLaQCniA/m2b8NLEgWQePsGNLy0h+0S+3yUZY4wvQi7gAQa3b8oLE1LZuv8YN01ZQs4pC3ljTPgJyYAHGN45gb/f0J+1u7K5dWoaJ/MK/S7JGGNqVMgGPMCY7i14+pq+pG0/xKRX0ziVbyFvjAkfIR3wAJf0SeRPV/Zh7uYD3PnGcvILi/wuyRhjakTIBzzAlQPa8PD4nsxYv5+731pBgYW8MSYM1KqLbntpwuBkcvMLeeTj9dSJjuCJK/sQESF+l2WMMZ4Jm4AH+MGw9pzMK+TJLzZRNzqSR8b3RMRC3hgTmsIq4AHuHN2RE/mFPDdrK7HRkdx3UTcLeWNMSAq7gBcRfnl+F07mFfLSvHTqxUTy8+918bssY4ypdmEX8OCE/P2XdCe3oJC/fbWFuDpR3D6ig99lGWNMtQrLgAcn5B8Z34ucUwU8/tkGurdqyPDOdtFvY0zoCIvDJMsSGSH86credGpen7vfXsGe7JN+l2SMMdUmrAMeoF5MFM/dOIDc/ELufONr+yGUMSZkhH3AA3RIqM9jV/Rm2fbDPP7pBr/LMcaYauFpwItIYxF5T0Q2iMh6ERniZXtn4pI+idw0JJkX56Xz2Zo9fpdjjDFnzOst+L8An6lqV6APsN7j9s7I7y7qRp82jfjFu6vs2q7GmKDnWcCLSCNgOPASgKrmqeoRr9qrDnWiIvn7Df2JiBB+/PpyO/ukMSaoebkF3w7IAqaIyNci8qKIxJUcSUQmiUiaiKRlZWV5WE7ltGlSj6ev6cO6PUd5YNpav8sxxpjT5mXARwH9gedUtR9wHPh1yZFU9QVVTVXV1ISE2nEc+uiuLfjxyA68tXQn7y3L9LscY4w5LV4GfCaQqaqL3fvv4QR+UPjZmM4Mbh/Pff9ZzYa9R/0uxxhjqsyzgFfVvcBOESk+0cu5wDqv2qtuUZER/PW6fjSIjebHry/nWG6B3yUZY0yVeH0UzV3A6yKyCugL/NHj9qpV8wax/PXafmQcOM6v3l+FqvpdkjHGVJqnAa+qK9z+9d6qOl5VD3vZnheGdGjKved34eNVe3hl4Xa/yzHGmEqzX7JWwu3DO3Bu1+Y88vE6Vuw84nc5xhhTKRbwlRARITx5dR+aN4jljteXc/h4nt8lGWNMhSzgK6lxvRievaE/WTm5/OydFRQVWX+8MaZ2s4Cvgj5Jjfm/i7sxc2MWz83e6nc5xhhTLgv4KrpxcDKX9EnkyekbWbD1gN/lGGNMmSzgq0hEePTyXrRrFsdP3lzB/qOn/C7JGGNKZQF/GurXcS4Scjy3gLve/Nr6440xtZIF/Gnq3KIBD4zrzuL0Q/z7611+l2OMMd9hAX8GrhqQRL+2jXns0w3knMr3uxxjjPkWC/gzEBEhPHBJDw4ez+VvX23xuxxjjPkWC/gz1CepMVcPSOLleels2X/M73KMMeYbFvDV4Bdju1A3OpKHPlpnJyQzxtQaFvDVoFn9Otw9pjNzNmXx5fr9fpdjjDGABXy1uWlIMh2b1+ehj9bZtVyNMbWCBXw1iY6M4P5LurPj0AlempfudznGGGMBX52GdUrg/B4teOarLezJPul3OcaYMGcBX83uu6g7Rao8+skGv0sxxoQ5C/hqlhRfj9tGdGDayt0sST/kdznGmDBmAe+BH43oQGKjWO6ftpZCO0+NMcYnFvAeqBsTye8u6s76PUd5c8kOv8sxxoQpTwNeRDJEZLWIrBCRNC/bqm0u7NWSwe3jeWL6Ro6csEv8GWNqXk1swY9S1b6qmloDbdUaIsID43pw9GQ+T32xye9yjDFhyLpoPNS1ZUMmDE7mtUXbWb/nqN/lGGPCjNcBr8B0EVkmIpM8bqtWumdMZxrVjeb+aWvtPDXGmBrldcCfo6r9gQuAO0RkeMkRRGSSiKSJSFpWVpbH5dS8xvViuPf8LixJP8RHq/b4XY4xJox4GvCqusv9vx/4ABhUyjgvqGqqqqYmJCR4WY5vrh3Ylh6JDfnjJ+s5kVfgdznGmDDhWcCLSJyINCi+DXwPWONVe7VZZITw4Lge7Mk+xXOztvpdjjEmTHi5Bd8CmCciK4ElwMeq+pmH7dVqqSnxjO+byPNztrHj4Am/yzHGhAHPAl5Vt6lqH/evh6r+wau2gsWvL+hGVITwyMfr/C7FGBMG7DDJGtSyUSx3ju7I9HX7mLMp9HYoG2NqFwv4GnbrOe1IblqPBz9cS35hkd/lGGNCmAV8DasTFcnvL+7O1qzj/GtBht/lGGNCmAW8D0Z3bc7ILglMnrGZA8dy/S7HGBOiLOB9ICLcd1F3jucV2OX9jDGesYD3Scfm9bmwVyteXbid7BP5fpdjjAlBFvA+unNUR47lFvCvhRl+l2KMCUEW8D7q1qoh53Vrzsvz0zmea6cwMMZULwt4n90xqiNHTuTz+uLtfpdijAkxFvA+69e2CUM7NuWfc9M5lV/odznGmBBiAV8L3DmqE1k5ubyTttPvUowxIcQCvhYY3D6eAclNeH72Nvt1qzGm2ljA1wIiwp2jOrLryEk++HqX3+UYY0KEBXwtMbJLAj0SG/LcrK0UFtml/YwxZ84CvpYo3opPP3Ccj1fbpf2MMWfOAr4WOb9HSzo2r8+zM7dQZFvxxpgzZAFfi0RECHeM6sCGvTl8uWG/3+UYY4KcBXwtc0nvRNrG1+OZrzajalvxxpjTZwFfy0RFRnD7iA6szMxm3pYDfpdjjAliFvC10BUDWtOyYSzPfLXF71KMMUHMAr4WqhMVyaTh7VmcfoilGYf8LscYE6Q8D3gRiRSRr0XkI6/bCiXXDWpL07gY24o3xpy2mtiC/ymwvgbaCSl1YyL5/jntmL0pi9WZ2X6XY4wJQp4GvIi0AS4CXvSynVB105BkGsZG8czMzX6XYowJQl5vwU8GfgmUeQYtEZkkImkikpaVleVxOcGlQWw0E89O4fO1+9i0L8fvcowxQcazgBeRi4H9qrqsvPFU9QVVTVXV1ISEBK/KCVq3DG1HvZhInp1pffHGmKrxcgt+KDBORDKAt4DRIvKah+2FpCZxMdw4OJlpK3eTceC43+UYY4KIZwGvqr9R1TaqmgJcC3ylqjd61V4o+8E57YiKjOAfs7f6XYoxJojYcfBBoHnDWK5JTeL95ZnsPnLS73KMMUGiRgJeVWep6sU10Vaoum1Ee1ThhTnb/C7FGBMkbAs+SLRpUo/L+rXmzSU7yMrJ9bscY0wQsIAPIj8a2YH8wiJempfudynGmCBgAR9E2ifU56Leiby6MIMjJ/L8LscYU8tVKuBFJE5EItzbnUVknIhEe1uaKc0dozpwPK+QqQsy/C7FGFPLVXYLfg4QKyKtgenABGCqV0WZsnVt2ZDzurVgyvwMjuUW+F2OMaYWq2zAi6qeAC4HnlXVq4Ae3pVlynPn6I5kn8znrSU7/C7FGFOLVTrgRWQIcAPwsftYpDclmYr0TWrMoJR4pi7IoNAuzm2MKUNlA/5u4DfAB6q6VkTaAzM9q8pU6JahKWQePsmX6/f5XYoxppaKqsxIqjobmA3g7mw9oKo/8bIwU74x3VvQunFdpszP4Hs9WvpdjjGmFqrsUTRviEhDEYkD1gDrROQX3pZmyhMVGcGEIcks3HaQDXuP+l2OMaYWqmwXTXdVPQqMBz4F2uEcSWN8dO3AJGKjI5g6P8PvUowxtVBlAz7aPe59PDBNVfMB27vns8b1YrisXxs++HoXh4/bD5+MMd9W2YB/HsgA4oA5IpIMWL9ALTDx7BRyC4p4c6kdMmmM+bZKBbyq/lVVW6vqherYDozyuDZTCV1aNmBox6a8unA7BYVlXhnRGBOGKruTtZGIPFV87VQReRJna97UAhPPbsee7FN8vtYOmTTG/E9lu2heBnKAq92/o8AUr4oyVTO6a3Paxtdj6gI7y6Qx5n8qG/AdVPV+Vd3m/j0ItPeyMFN5kRHCTUOSWZpxmDW7sv0uxxhTS1Q24E+KyDnFd0RkKGDXjqtFrkpNol5MJFPskEljjKuyAX878HcRyRCRDOAZ4DbPqjJV1qhuNFcOaMOHK3fbFZ+MMUDlj6JZqap9gN5Ab1XtB4z2tDJTZTefnUJeYRFv2lkmjTFU8YpOqnrU/UUrwM/KG1dEYkVkiYisFJG1IvLgaVdpKqVDQn1GdE7g1UXbySuwQyaNCXdncsk+qWB4LjDa3fLvC4wVkcFn0J6phIlDU8jKyeXTNXv8LsUY47MzCfhyT1Xg/iDqmHs32v2z0xt4bESnBNo3i+Nl29lqTNgrN+BFJEdEjpbylwMkVjRxEYkUkRXAfuALVV1cyjiTin9AlZWVdbrzYVwREcLNZ6ewcucRvt5x2O9yjDE+KjfgVbWBqjYs5a+BqlZ4LnlVLVTVvkAbYJCI9CxlnBdUNVVVUxMSEk57Rsz/XDGgDQ3qRNkhk8aEuTPpoqk0VT2CcwWosTXRXrirXyeKq1KT+GT1HvYdPeV3OcYYn3gW8CKSICKN3dt1gTHABq/aM99289nJFKry2qLtfpdijPGJl1vwrYCZIrIKWIrTB/+Rh+2ZAMlN4zi3a3PeWLyDU/mFfpdjjPGBZwGvqqtUtZ+q9lbVnqr6kFdtmdLdMrQdB4/n8eHK3X6XYozxQY30wRt/nN2hKZ1b1GfqggxU7QhVY8KNBXwIExEmnt2OtbuPsjTDDpk0JtxYwIe4y/q1plHdaDtXvDFhyAI+xNWNieTaQUl8vnYfu47YGZ6NCScW8GFgwuBkVJVXF9ohk8aEEwv4MNCmST3O79GSt5bu4GSeHTJpTLiwgA8TE89O4ciJfP6zYpffpRhjaogFfJgY1C6e7q0aMnW+HTJpTLiwgA8TIsLEoSls3JfDwq0H/S7HGFMDLODDyLg+icTHxTBlQYbfpRhjaoAFfBiJjY7k+kFtmbF+HzsOnvC7HGOMxyzgw8yNg5OJFGGqbcUbE/Is4MNMy0axXNS7FW8v3UH2yXy/yzHGeMgCPgxNGt6e43mFvL7YfvhkTG2QW+DN71Ms4MNQj8RGnNOxGVPmZ3i2YhljynfkRB5vLdnB9f9cxKXPzPfk8OUKr6tqQtOk4e256eUl/Pfr3Vw9MMnvcowJC8dyC5ixbh/TVu5mzqYsCoqUds3iuKRPIgVFSnSkVGt7FvBhalinZnRr1ZAX5m7jygFtiIio3hXLGOM4lV/IrI37+XDlHr7csI9T+UUkNorl1nPacUmfRHokNkTEm/efBXyYEhFuG96eu99ewcyN+zm3Wwu/SzImZOQXFjFvywE+XLmb6Wv3cSy3gGb1Y7gmNYlL+iTSv22TGtmosoAPYxf1bsWfPtvA83O2WcAbc4YKi5Ql6Yf4cNVuPl29h8Mn8mkYG8WFvVoyrk9rBrePJyqyZnd7WsCHsejICL5/Tjse+Xg9K3YeoW9SY79LMibobN6XwztpO5m2cjf7juZSNzqSMd1bMK5PIsM6N6NOVKRvtXkW8CKSBLwCtAAUeEFV/+JVe+b0XDuoLX/5cjMvzNnKszcM8LscY4LC8dwCPl61h7eW7mD5jiNERwqjujRnXN9ERndtTr2Y2rHt7GUVBcDPVXW5iDQAlonIF6q6zsM2TRXVrxPFjYOTeX72VrYfPE5y0zi/SzKmVlJVvt55hHeW7uTDlbs5nldIx+b1ue+iblzWrzVN69fxu8Tv8CzgVXUPsMe9nSMi64HWgAV8LXPL2Sm8NDedF+em8/D4nn6XY0ytcuh4Hv9ensk7aTvZtO8YdaMjuaRPK64ZmET/tk08OwKmOtTI9wgRSQH6AYtroj1TNc0bxjK+XyLvLtvJPWM6Ex8X43dJxviqqEiZt+UAby/dyfR1e8kvVPomNeaxy3txcZ9E6tepHV0wFfG8ShGpD7wP3K2qR0sZPgmYBNC2bVuvyzFlmDS8Pe+kZfLKwgzuPq+z3+UY44tdR07ybtpO3k3LZNeRkzSuF82EwSlcMzCJLi0b+F1elXka8CISjRPur6vqv0sbR1VfAF4ASE1NtUsN+aRj8wac27U5ryzczm3DO1A3xr89/8bUpKIiZfamLKYuyGDO5ixUnR8C/ubCrozp3sLXo2DOlJdH0QjwErBeVZ/yqh1TfW4b0YGrn1/Ie8szmTA42e9yjPFUzql83luWyb8WZJBx8ATNG9ThrlEduSo1iaT4en6XVy283IIfCkwAVovICvex36rqJx62ac7AwJQm9E1qzItzt3H9oLZE2ukLTAjKOHCcqQsyeG9ZJsdyC+jXtjH3jOnMBT1bERMVWudf9PIomnmAJUQQKT59wY9eX87na/dyYa9WfpdkTLVQdXaaTp2fwVcb9xMVIVzUqxUTh7YL6R/4BceuYFNjvtejJSlN6/H8nG1c0LNlrT4EzJiKnMgr4N/LdzF1QQZb9h+jWf0Y7hrdiRvPakvzhrF+l+c5C3jzLZERwq3D2vN//1nDkvRDnNW+qd8lGVNlOw+d4NVF23lryQ6OniqgZ+uGPHlVHy7u0yqod5pWlQW8+Y6rBrTh6S828cKcbRbwJmioKovTDzFlfjpfrNuHiDC2Z0tuOTuFAcm1+wdJXrGAN98RGx3JTUOSmTxjM5v35dCpRfAd/2vCy9KMQ/z5840sST9Ek3rR3D6iAzcOTiaxcV2/S/OVBbwp1U1DUvjH7K38c+42/nRlH7/LMaZUqzOzeWL6RmZvyiKhQR0eHNeDawYmERsdPt0w5bGAN6WKj4vhqgFJvL10J/d+r0tY7JAywWPzvhyenL6Jz9bupXG9aH5zQVduGpJiP9ArwQLelOkHw9rx+uLtTFmQwa/GdvW7HGPYcfAEk2ds4oMVu4iLieKn53bi1mHtaBgb7XdptZIFvClTctM4xvZsyWuLtnPHqI5Bc4IlE3r2Zp/ir19t5p2lO4mMECYNa89tIzrYifEqYO9YU67bhnfgk9V7eWvJDn4wrL3f5Zgwc/BYLs/N2sori7ajqlw3qC13ju5IC+syrBQLeFOuPkmNOatdPC/PS+fms1OIruFrSprwdPRUPi/O2cZL89I5mV/IZf3acPd5nULmHDE1xQLeVOi2Ee35/tQ0Pl61h/H9WvtdjglhJ/MKmbIgnednbyP7ZD4X9WrFPWM60bG5Hap7OizgTYVGdm5Op+b1eX7ONi7tmxiWPxgx3lJVPly1h0c/Wc+e7FOM7tqcn43pTM/WjfwuLahZwJsKRUQIPxzenl++t4p5Ww4wrFOC3yWZELI6M5sHP1xL2vbD9EhsyORr+tovqKuJBbyplEv7JvLE5xt5fvY2C3hTLfbnnOKJzzfy7rJMmsbF8NjlvbgqNclOU12NLOBNpdSJiuSWoe14/LMNrNmVbV+dzWnLLShk6vwM/vbVFnILCvnhsPbcObqjHcvuATskwlTa9We1JS4mkufnbPO7FBOEVJUZ6/Zx/tNzePTTDZzVLp7P7x7Oby/sZuHuEduCN5XWqG40E4em8PeZW7kmNYlzOjXzuyQTJDbvy+Ghj9Yxd/MBOiTEMfWWgYzs0tzvskKeBbypkrtGd+LTNXv51fur+OzuYTSwLS9TjiMn8pg8YzOvLtpOvZhIfn9xdyYMSbbfU9QQW8qmSmKjI3niqj7syT7JHz9Z73c5ppYqKCzi1YUZjHxiFq8szODagUnMunck3z+nnYV7DbIteFNl/ds24YfD2/P87G1c0LMVwzvbUTXGoarM2XyARz9Zz4a9OQxuH8/vL+5B98SGfpcWlizgzWm557zOfLl+P796fxWf3zPcdpKFOVVl7uYDTJ6xieU7jtCmSV2eu6E/Y+26vr7y7LuSiLwsIvtFZI1XbRj/FHfV7Dt6ij98ZF014UpVmb0pi8ufW8BNLy9hb/YpHhnfky9/PoILerWycPeZl1vwU4FngFc8bMP4qG9SY24f0YFnZ21lbK+WjLKjIsJGcbD/5cvNfL3jCImNYvnDZT25ckCbsLqodW3nWcCr6hwRSfFq+qZ2+Ol5nZixfh+/eX81n98znEZ1rasmlBUH++QZm1mx04K9tvN9d7aITBKRNBFJy8rK8rscU0V1opyumqxjuTz80Tq/yzEeUVVmbtzP+GcXMHHKUrJycvnjZb2Y9YtR3HBWsoV7LeX7TlZVfQF4ASA1NVV9Lsecht5tGvOjER14ZuYWLuzVktFdW/hdkqkmqsosd4t95c4jtG5clz9e1osrB7QhJsr37UNTAd8D3oSGu87tyIz1+/j1+6v54p54GtWzrppgpqrM2pjF5BmbWJmZTevGdXn08l5c0d+CPZhYwJtqUdxVc+nf5/Pgh2t56pq+fpdkTkN+YRGfrN7Di3PTWb3LCfbHLu/F5RbsQcmzgBeRN4GRQDMRyQTuV9WXvGrP+K9n60bcMaojf/1yMxf0asWY7tZVEyyyT+bz1pIdTF2QwZ7sU7RrFmdb7CHAy6NorvNq2qb2unNUR75Yt4/ffrCagSlNaFzPrnpfm+04eIKX56fzTtpOTuQVMqR9Ux4Z35NRXZoTYedlD3rWRWOqVUxUBE9c1ZtLn5nPA9PWMvnafn6XZEpQVdK2H+aluel8vm4vkSKM65PI989pZ+f5DzEW8Kba9UhsxJ2jOzJ5xmbG9mzF2J4t/S7J4PSvf7pmLy/N3cbKzGwa1Y3mRyM6cNOQFFo2ivW7POMBC3jjiTtGdWT62n3c95/VDGoXT3ycddX4JftkPm8v3cHU+RnsdvvXHx7fkyv6t6ZejEVAKLNX13giOjKCJ6/uw7hn5nH/tLX87TrrqqlpOw6eYMqCdN5ZupPjeYUMbh/PQ5f2ZHRX618PFxbwxjPdWjXkJ6M78eQXm7iwZ0su6NXK75JCnqqycNtBpszPYMb6fda/HuYs4I2nbh/Zgenr9nHff9YwqF08TevX8bukkHQqv5D/rtjFlPkZbNibQ5N61r9uLOCNx6IjI3jiqj5c8rd5/P6/a/n7Df39Limk7Mk+yWuLtvPG4h0cPpFP15YNePyKXlzatzWx0XZ+mHBnAW8816VlA356Xif+/PlGLli1m4t7J/pdUlBTVZbvOMKU+el8tmYvhaqM6daCW4a2Y3D7eDsHu/mGBbypEbcNb8/0tXv53QdriIoQzu9hV/qpqrwC5zQCU+anszIzmwaxUdwyNIWbhqSQFF/P7/JMLWQBb2pEVGQEf7m2Hz96fTm3v7acc7s258FLe9CmiQVTRQ4cy+X1RTt4bfF2snJyaZ8Qx8OX9uDy/m2Iq2NvYVM2Ua09Z+hNTU3VtLQ0v8swHiooLGLK/Aye+mITAD8b05lbhqYQFWnnOylp5c4jvLJwOx+u3E1eYREjuyRwy9B2DOvYzA5zNN8QkWWqmlrqMAt444fMwye4/79r+XLDfrq1asgfL+tJv7ZN/C7LdyfzCvlw5W5eW7ydVZnZ1IuJ5MoBbbj57BQ6JNT3uzxTC1nAm1pJVfl87V4emLaOfTmnuPGsZH4xtgsNY8PvXPJb9h/j9cXbeX9ZJkdPFdC5RX0mDE5mfL/WNAjD5WEqr7yAtw484xsRYWzPVgzt2Iwnp2/ilYUZfL52L/df0oMLe4X+Ttj8wiK+WLeP1xZtZ8HWg0RHChf0bMWNg5MZmNIk5OffeM+24E2tsSrzCL/9YDVrdh1lZJcEHr60Z0geHbIn+yRvLt7BW0t3sj8nl9aN63L9WW25OjWJhAb2QzBTNdZFY4JGQWER/1q4naemb6RQlbvP68yt57QjOsh3whYVKfO2HOC1Rdv5csN+ilQZ2TmBCUOSGdG5OZG209ScJgt4E3R2HznJA9PWMn3dPrq2bMAfLuvFgOTg2wm7/+gp/rtiN68v3k7GwRM0jYvh6oFJXD+obUh+OzE1zwLeBK3pa/dy/7S17D16imsHtuXKAa3p2boRdaJq58/wC4uUlZlHmLlhPzM37mfNrqMADExpwo2Dkxnbs2Wtrd0EJwt4E9SO5Rbw9BebmDI/nSKFOlER9ElqzMCUJgxMiad/chNfj7w5ciKP2ZuymLUxi9mbsjh0PI8IgQHJTRjZpTnf696CTi0a+FafCW0W8CYkHDyWS9r2wyxNP8TS7YdZsyubwiIlQqBry4ZO4LeLZ2BKPC0aencGRVVl/Z4cZm7cz8wN+1m+4zBFCvFxMYzonMCors0Z3qmZXY/W1AjfAl5ExgJ/ASKBF1X1sfLGt4A3VXEir4AVO46wJOMQaRmHWb7jMCfyCgFoG1+P1JQmDEqJJzUlng4JcWd02OHx3ALmbTnArI37mbkhi71HTwHQq3UjRnVxQr13m8a2s9TUOF8CXkQigU3AGCATWApcp6rrynqOBbw5E/mFRazfc5Ql6U7gL804xMHjeYCzdX26lw1UVXYeOkleYRH160QxrFMzRnVtzsjOCTT38JuCMZXh1w+dBgFbVHWbW8RbwKVAmQFvzJmIjoygd5vG9G7TmB8Mc4J524HjpGUcYtn2wxzPLTyt6SrKud1aMLJLAqnJ8cREBfchmyZ8eBnwrYGdAfczgbNKjiQik4BJAG3btvWwHBNuRIQOCfXpkFCfawbaumXCj++bIqr6gqqmqmpqQkKC3+UYY0zI8DLgdwFJAffbuI8ZY4ypAV4G/FKgk4i0E5EY4FpgmoftGWOMCeBZH7yqFojIncDnOIdJvqyqa71qzxhjzLd5erpgVf0E+MTLNowxxpTO952sxhhjvGEBb4wxIcoC3hhjQlStOtmYiGQB2/2uo4qaAQf8LqKG2TyHB5vn4JCsqqX+iKhWBXwwEpG0ss4DEapsnsODzXPwsy4aY4wJURbwxhgToizgz9wLfhfgA5vn8GDzHOSsD94YY0KUbcEbY0yIsoA3xpgQZQFfRSISLyJfiMhm93+TcsZtKCKZIvJMTdZY3SozzyLSV0QWishaEVklItf4UeuZEpGxIrJRRLaIyK9LGV5HRN52hy8WkRQfyqxWlZjnn4nIOvd1/VJEkv2os7pUNL8B410hIioiQXvYpAV81f0a+FJVOwFfuvfL8jAwp0aq8lZl5vkEcJOq9gDGApNFpHHNlXjm3OsI/x24AOgOXCci3UuMditwWFU7Ak8Dj9dsldWrkvP8NZCqqr2B94A/1WyV1aeS84uINAB+Ciyu2QqrlwV81V0K/Mu9/S9gfGkjicgAoAUwvWbK8lSF86yqm1R1s3t7N7AfCLZLdH1zHWFVzQOKryMcKHBZvAecKyJSgzVWtwrnWVVnquoJ9+4inIv3BKvKvMbgbJw9DpyqyeKqmwV81bVQ1T3u7b04If4tIhIBPAncW5OFeajCeQ4kIoOAGGCr14VVs9KuI9y6rHFUtQDIBprWSHXeqMw8B7oV+NTTirxV4fyKSH8gSVU/rsnCvODp+eCDlYjMAFqWMuh3gXdUVUWktONMfwx8oqqZwbJxVw3zXDydVsCrwM2qWlS9VRo/iciNQCowwu9avOJunD0FTPS5lGphAV8KVT2vrGEisk9EWqnqHjfM9pcy2hBgmIj8GKgPxIjIMVUtr7/eV9Uwz4hIQ+Bj4HequsijUr1UmesIF4+TKSJRQCPgYM2U54lKXTtZRM7D+bAfoaq5NVSbFyqa3wZAT2CWu3HWEpgmIuNUNa3Gqqwm1kVTddOAm93bNwP/LTmCqt6gqm1VNQWnm+aV2hzulVDhPLvX3f0AZ17fq8HaqlNlriMcuCyuBL7S4P61YIXzLCL9gOeBcapa6od7ECl3flU1W1WbqWqK+/5dhDPfQRfuYAF/Oh4DxojIZuA89z4ikioiL/pamXcqM89XA8OBiSKywv3r60u1p8ntUy++jvB64B1VXSsiD4nIOHe0l4CmIrIF+BnlH0VV61Vynv+M8030Xfd1LfmhFzQqOb8hw05VYIwxIcq24I0xJkRZwBtjTIiygDfGmBBlAW+MMSHKAt4YY0KUBbwJKiJS6B6qt0ZE3hWRelV4boqIXH+a7WaISLMqjD9SRM4OuH+7iNx0Om0bc7os4E2wOamqfVW1J5AH3F6ZJ7m/Ok0Bqhzw7hkIq2ok8E3Aq+o/VPWV05iOMafNAt4Es7lARxGJE5GXRWSJiHwtIpcCiMhEEZkmIl/hnOb4MZxTSKwQkXvc4d+cq19EPhKRke7tYyLypIisxDn1BMAvRWS1205Hd7xL3PPCfy0iM0SkhXuO+NuBe9y2honIAyJyr/ucviKyyD2/+gfinl9fRGaJyOPu9DeJyLAaWIYmhFnAm6DkbpFfAKzGOUfKV6o6CBgF/FlE4txR+wNXquoInF+dznW/ATxdQRNxwGJV7aOq89zHslW1F/AMMNl9bB4wWFX74Zx69peqmgH8A3jabWtuiWm/AvzKPb/6auD+gGFR7nzcXeJxY6rMTjZmgk1dEVnh3p6Lc+qABcC44i1kIBZo697+QlUPnUY7hcD7JR57M+B/8QdEG+Bt9yRsMUB6eRMVkUZAY1Wd7T70L+DdgFH+7f5fhtOlZMxps4A3weakqvYNfMC94MYVqrqxxONnAcfLmVYB3/4WGxtw+5SqFpYYX0u5/TfgKVWd5nbvPFBB/RUpPlNjIfb+NGfIumhMKPgcuKv4ykru2Q9Lk4NzOthiGUBfEYkQkSScq/2U55qA/wvd24343+lmbw4Yt2RbgHO2QuBwQP/6BGB2yfGMqQ62hWBCwcM4feKr3As2pAMXlzLeKqDQ3XE61X1OOrAO58yCyytop4mIrMLZyr7OfewBnLMsHga+Atq5j38IvOfu8L2rxHRuBv7hHuK5DbilMjNpTFXZ2SSNMSZEWReNMcaEKAt4Y4wJURbwxhgToizgjTEmRFnAG2NMiLKAN8aYEGUBb4wxIer/AemGEoQbxj9SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pyhessian.utils import normalization\n",
    "\n",
    "\n",
    "# used to perturb your model \n",
    "lams = np.linspace(-0.5, 0.5, 21).astype(np.float32)\n",
    "\n",
    "loss_list = []\n",
    "\n",
    "# create a copy of the model\n",
    "model_perb = ptcv_get_model(\"resnet20_cifar10\", pretrained=True)\n",
    "model_perb.eval()\n",
    "model_perb = model_perb.cuda()\n",
    "\n",
    "# generate gradient vector to do the loss plot\n",
    "loss = criterion(model_perb(inputs), targets)\n",
    "loss.backward()\n",
    "\n",
    "v = [p.grad.data for p in model_perb.parameters()]\n",
    "v = normalization(v)\n",
    "model_perb.zero_grad()\n",
    "\n",
    "\n",
    "for lam in lams: \n",
    "    model_perb = get_params(model, model_perb, v, lam)\n",
    "    loss_list.append(criterion(model_perb(inputs), targets).item())\n",
    "\n",
    "plt.plot(lams, loss_list)\n",
    "plt.ylabel('Loss')\n",
    "plt.xlabel('Perturbation')\n",
    "plt.title('Loss landscape perturbed based on gradient direction')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Hessian Trace/Diagonal\n",
    "We can also use randomized linear algebra to compute Hessian trace or approximate the Hessian diagonal with very little computational overhead. Let's first start with a numpy example, and then we will show results for a NN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000 # the matrix size\n",
    "\n",
    "# generate the matrix\n",
    "A = np.random.randn(n, n)\n",
    "B = A @ A.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The trace of B is: 998508.2006\n"
     ]
    }
   ],
   "source": [
    "# Direct get the trace \n",
    "trace_B_np = np.matrix.trace(B)\n",
    "print(\"The trace of B is: %.4f\"%trace_B_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can approximate the above by using Hutchinson Method. It is very similar to power iteration:\n",
    "\n",
    "$$Tr(B) = \\mathbb{E}[v^TBv],$$\n",
    "$$Diag(B) = \\mathbb{E}[v \\bigodot Bv].$$\n",
    "\n",
    "It can be proved that the above expectation converges with smallest variance to the trace if we use Rademacher random numbers (+/-1). In practice you can also use Gaussian random vectors as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 01, Current estimated trace: 1012323.7 relative error: 1.4e-02\n",
      "Step 02, Current estimated trace: 1022530.3 relative error: 2.4e-02\n",
      "Step 03, Current estimated trace: 1000293.8 relative error: 1.8e-03\n",
      "Step 04, Current estimated trace: 1013886.1 relative error: 1.5e-02\n",
      "Step 05, Current estimated trace: 997908.3 relative error: -6.0e-04\n",
      "Step 06, Current estimated trace: 1004918.3 relative error: 6.4e-03\n",
      "Step 07, Current estimated trace: 990709.9 relative error: -7.8e-03\n",
      "Step 08, Current estimated trace: 985578.4 relative error: -1.3e-02\n",
      "Step 09, Current estimated trace: 980355.0 relative error: -1.8e-02\n",
      "Step 10, Current estimated trace: 985012.5 relative error: -1.4e-02\n",
      "Step 11, Current estimated trace: 992682.2 relative error: -5.8e-03\n",
      "Step 12, Current estimated trace: 997489.0 relative error: -1.0e-03\n",
      "Step 13, Current estimated trace: 996441.0 relative error: -2.1e-03\n",
      "Step 14, Current estimated trace: 998894.3 relative error: 3.9e-04\n",
      "Step 15, Current estimated trace: 998603.9 relative error: 9.6e-05\n",
      "Step 16, Current estimated trace: 999475.5 relative error: 9.7e-04\n",
      "Step 17, Current estimated trace: 1000506.1 relative error: 2.0e-03\n",
      "Step 18, Current estimated trace: 998818.7 relative error: 3.1e-04\n",
      "Step 19, Current estimated trace: 997798.1 relative error: -7.1e-04\n",
      "Step 20, Current estimated trace: 992688.2 relative error: -5.8e-03\n"
     ]
    }
   ],
   "source": [
    "# use Hutchinson method to get the trace of B\n",
    "trace_list = []\n",
    "\n",
    "for i in range(20):\n",
    "    v = np.random.randint(2, size=n) \n",
    "    v = v.reshape(n, 1) * 2 - 1 # Create Rademacher random numbers\n",
    "    trace_list.append(v.T @ B @ v)\n",
    "    trace_B_hutchinson = np.mean(trace_list)\n",
    "    print(\"Step %.2d, Current estimated trace: %.1f relative error: %.1e\"\n",
    "          %(i+1, trace_B_hutchinson, (trace_B_hutchinson - trace_B_np) / trace_B_np))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see we can get a very accurate estimate of the trace. Next let's try to approximate the diagonal of B using the matrix-free Hutchinson's method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 01, the current average relative error 8.3e-01:\n",
      "Step 02, the current average relative error 5.5e-01:\n",
      "Step 03, the current average relative error 4.5e-01:\n",
      "Step 04, the current average relative error 3.8e-01:\n",
      "Step 05, the current average relative error 3.4e-01:\n",
      "Step 06, the current average relative error 3.2e-01:\n",
      "Step 07, the current average relative error 3.0e-01:\n",
      "Step 08, the current average relative error 2.7e-01:\n",
      "Step 09, the current average relative error 2.6e-01:\n",
      "Step 10, the current average relative error 2.4e-01:\n",
      "Step 11, the current average relative error 2.3e-01:\n",
      "Step 12, the current average relative error 2.2e-01:\n",
      "Step 13, the current average relative error 2.1e-01:\n",
      "Step 14, the current average relative error 2.0e-01:\n",
      "Step 15, the current average relative error 2.0e-01:\n",
      "Step 16, the current average relative error 1.9e-01:\n",
      "Step 17, the current average relative error 1.9e-01:\n",
      "Step 18, the current average relative error 1.8e-01:\n",
      "Step 19, the current average relative error 1.7e-01:\n",
      "Step 20, the current average relative error 1.7e-01:\n"
     ]
    }
   ],
   "source": [
    "# use Hutchinson method to get the diag of B\n",
    "diag_est = np.zeros([n, 1])\n",
    "diag_B_np = np.diag(B)\n",
    "for i in range(20):\n",
    "    v = np.random.randint(2, size=n)\n",
    "    v = v.reshape(n, 1) * 2 - 1\n",
    "    diag_est += np.multiply(v, (B @ v))\n",
    "    diag_est_err = np.mean(np.abs(diag_est.reshape(-1) / (i+1) - diag_B_np) / diag_B_np)\n",
    "    print(\"Step %.2d, the current average relative error %.1e:\"%(i+1,diag_est_err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's repeate the above for computing the trace and diagonal of Hessian for ResNet20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Files already downloaded and verified\n"
     ]
    }
   ],
   "source": [
    "# get the model \n",
    "model = ptcv_get_model(\"resnet20_cifar10\", pretrained=True)\n",
    "# change the model to eval mode to disable running stats upate\n",
    "model.eval()\n",
    "\n",
    "# create loss function\n",
    "criterion = torch.nn.CrossEntropyLoss()\n",
    "\n",
    "# get dataset \n",
    "train_loader, test_loader = getData()\n",
    "\n",
    "# for illustrate, we only use one batch to do the tutorial\n",
    "for inputs, targets in train_loader:\n",
    "    break\n",
    "\n",
    "# we use cuda to make the computation fast\n",
    "model = model.cuda()\n",
    "inputs, targets = inputs.cuda(), targets.cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the hessian computation module\n",
    "hessian_comp = hessian(model, criterion, data=(inputs, targets), cuda=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The trace of this model is: 638.2978\n"
     ]
    }
   ],
   "source": [
    "trace = hessian_comp.trace()\n",
    "print(\"The trace of this model is: %.4f\"%(np.mean(trace)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get the full eigenvalue spectrum density of Hessian using Stochastic Lancoz algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "density_eigen, density_weight = hessian_comp.density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxWklEQVR4nO3debxkVXnv/8/3nNOTDQgoIKICIoPIDeBtjYY4K2i8Rr2YBFATSa7EKRonoleMRHEMmvyMiGJEEBQvJkhUnAdUHKINiNKIIEODjA3d9Nx9hnp+f+xd5+yqU8OuOruqdtX5vl+v8+pde6q165zeT621nr2WIgIzM7OyGRt0AczMzBpxgDIzs1JygDIzs1JygDIzs1JygDIzs1KaGHQByuShD31oHHDAAYMuhpnZonLllVfeFxF71a93gMo44IADWL169aCLYWa2qEha22i9m/jMzKyUHKDMzKyUHKDMzKyUHKDMzKyUHKDMzKyUHKDMzKyUHKDMzKyUHKAWYOvOae7bsnPQxTAzG0kOUAvwdxddzaozvjPoYpiZjSQHqAX43vX3DroIZmYjywHKzMxKyQHKzMxKyQHKzMxKyQHKzMxKyQHKzMxKyQHKzMxKyQHKzMxKyQGqABEx6CKYmY0cB6gCVByfzMwK5wBVgIprUGZmhXOAKoADlJlZ8RygCuD4ZGZWPAeoAsy4E8rMrHAOUAVwE5+ZWfEcoArgCpSZWfEcoArg56DMzIrnAFUA90GZmRXPAaoAjk9mZsVzgCqAm/jMzIo30gFK0gcl/UjSBZKW9Op9XIMyMyveyAYoSUcC+0XEU4DrgZf06r2cZm5mVryRDVDAHwHfSpe/ARzTqzdykoSZWfFKH6AkvU7Sakk7JZ1Xt21PSV+StFXSWkknZTbvAWxKlzcCe/aqjK5AmZkVb2LQBcjhTuAM4DhgRd22s4BJYB/gKOAySddExBrgAWC3dL8HA+t7VUA38ZmZFa/0NaiIuCQiLgXuz66XtBI4HnhnRGyJiCuALwMvT3f5CfDsdPk44MeNzi/plLSGtnrdunVdldEBysyseKUPUC0cAkxHxA2ZddcAjwOIiF8C90j6UbruPxudJCLOiYhVEbFqr7326qogDlBmZsUbhia+ZnZhro+paiOwa/VFRLy1HwVxjoSZWfGGuQa1hbk+pqrdgM39LohrUGZmxRvmAHUDMCHp4My6I4E1/S5IpdLvdzQzG32lD1CSJiQtB8aBcUnLJU1ExFbgEuDdklZKOgZ4IXBBv8voGpSZWfFKH6CA04DtwNuAl6XLp6XbXkOSen4vcBHw6jTFvK8coMzMilf6JImIOB04vcm29cCL+lichpwkYWZWvGGoQZWea1BmZsVzgCqAp9swMyte2yY+SbsD/xt4GnAASZ/POuAq4OsR8ZMelm8ozDiLz8yscE1rUJIeLunfgbuAdwBLgdUkI4SvJQlY35Z0naS/6Edhy8pNfGZmxWtVg/olcD6wqllmnKQVJEkKb5L0yIg4s/ASDgEHKDOz4rUKUI+LiJajp0bEdpL07oskdTeQ3QhwfDIzK17TJr52wWmh+48ST1hoZla83Fl8kp4n6atpn9Mj03X/R9Kzele84eAmPjOz4uUKUJJeClwM3AgcCCxJN40Dp/amaMPD8cnMrHh5a1CnAq+MiDcC05n1PyOZyXZRcw3KzKx4eQPUwcBPG6xvNOXFouMuKDOz4uUNUHeSzGBb76nATcUVZ3hkR49wkoSZWfHyBqhzgI+mU1oAPFLSXwEfAs7uSclKLtuq56GOzMyKl2s084j4kKQHA98GlgPfB3YCZ0bEWT0sX2ll+51cgTIzK17u6TYi4h2S3gscTlLzui4itvSsZCWXDUpOkjAzK15H80FFxDaS8fgWvdoalAOUmVnRmgYoSV/Oe5KI+NNiijM8HKDMzHqrVQ3q/r6VYgjVNPF5ug0zs8I1DVARcXI/CzJsXIMyM+stz6jbpcjUmhygzMyKlztJQtIzgBOBR5FMXjgrIp5ZcLlKz2nmZma9lXew2FcAXwd2BZ5OMuX7HsDjget6VLZScxOfmVlv5W3iewvwuog4EZgC3h4RRwMXkozHt+jUPgc1uHKYmY2qvAHq0cB30uWdwC7p8seAVxRcpqGQHd7IQx2ZmRUvb4C6n6R5D+AO4Ih0+SHAiqILNQyytSYPFmtmVry8SRI/Ao4Ffk0yceFHJT0HeBbJ+HyLjpMkzMx6K2+Aeh3JILEA7yeZtPAYkmB1Rg/KVXoVN/GZmfVU3tHM12eWK8AHe1aiIVHxc1BmZj2VN838zyS9sMH6P5X0kuKLVX4zNRMWDrAgZmYjKm+SxOnAjgbrt6XbFp1sYoRrUGZmxeskzfy3Ddb/Lt226LgPysyst/IGqA3AwQ3WHwJsLq44w2N6xll8Zma9lDdA/RfwL5IOqa6QdCjwEeDSHpSr9DzUkZlZb+UNUP8AbASuk3S7pNuBNcAm4K29KlyZ1fRBuQplZla4vGnmm4Bj0odzj0pXXw18NxZpB8x0xU18Zma9lHu6DYCI+DaLdOSIem7iMzPrrZZNfJKOTOeByq57qaSbJd0r6ROSljY7fpRNZR5+cg3KzKx47fqgzgD+uPpC0uHAZ4AbgYuAl5L0T5WSpCdK+qmkH0q6SNKSos49NeMalJlZL7ULUI8HvpV5fQJwXUQcFxFvAP4e+Isela0ItwPPjIinArcC80bD6NZ0tgblKpSZWeHaBaiHAHdmXj8V+Erm9eUkU8CXUkTcFRHb05eTQGGDErmJz8yst9oFqHXAfgCSxoH/Cfx3ZvtSCrrpS3qdpNWSdko6r27bnpK+JGmrpLWSTurw3PuTTBfylXb75pVt4jv3x7cUdVozM0u1y+K7HHiXpNcC1UFhv5/ZfjhJ01kR7iTp8zqO+ZMgnkVSA9qHJM39MknXRMQaSQ8DvtDgfCdExN2SdgMuAF4REVMFlZXpikeINTPrpXYB6p0kU73/DpgBXh8RWzPbXw58t4iCRMQlAJJWAY+orpe0EjgeOCIitgBXSPpy+t5vi4i7gac3OqekCZLg9U8R0Wgswa5NTc/VoB6z9y5FntrMzGjTxBcRtwKHAUcD+0fE2XW7vAt4X2+KNusQYDoibsisuwZ4XI5jTwT+EHinpMslzUvokHRK2rS4et26dbkLNZXWoB790JU8ZOWizLQ3M+uptg/qRsQ0SUBotK3h+oLtQjKkUtZGYNd2B0bEBSTNe632OQc4B2DVqlW50x2mppMAtWzJeM2wR2ZmVoy8Y/EN0hZgt7p1uzHgUdS/85t7AVixZIwpBygzs8INQ4C6AZiQlJ3u40iSwWoH5orf3QfAiqXjzDhhwsyscKUJUJImJC0HxoFxScslTaRJGZcA75a0UtIxJA/ctmy665dtkzNce0d9C6SZmS1UaQIUcBqwHXgb8LJ0+bR022tIUs/vJRli6dURMdAaVNXVtz0A1D64a2ZmC5drNHNJzUaLCGBHRORPf2t2oojTgdObbFsPvGih71Gkx+67GzumZrjlviTrfnomWDI+4EKZmY2QvNNt3EoSjBqStIlkENlT06y/kfebu2qb9WY8YKyZWaHyBqgTgQ8Bn2BuqKM/BE4hqfXsTtIct5nk2ahFZ2bGAcrMrEh5A9SrgTdWR3tIfU/Sb4E3RMTTJN0L/BOLLEC9+ukHcfblN7kGZWZWsLxJEn8I/LrB+muBJ6TLPyUzRNEou+OBZID0JxywBw/fPRk20A/rmpkVK2+AWkvSnFfvlcBt6fJewPoiClV2P73pfgB+cesGxiXAAcrMrGh5m/jeDPynpD8BfpGuWwUcRDKQKyQ1qYuLLV45ZWfQnRhLA5Sb+MzMCpUrQEXEZelIDq8BDk1Xfxn4RETclu7z8d4UsXyqQWliTIylyzunZgZZJDOzkZO3BkVE3A68vYdlGRp7pKOXf+ovV/G965Mx+Z754R9w6weeP8himZmNlNwBStI+wGtJJikMkrHwzo6Ie3pUttKaTEcy32vXZazfNjng0piZjaZcSRLp+He/A04iGYJoB8lwRDdKenLvildOO9LmvBVLxxlLkyTMzKxYeWtQZ5KMgfeqiKgASBojeXD3w8Af9aZ45TQboJaMM+b4ZGbWE3nTzI8CPlwNTgDp8kdIZttdVLZPJgFq+ZJxjnj4gwdcGjOz0ZQ3QG0EDmyw/kDggcJKMyS2TyVxesWScf7PU5KP5eC9dxlkkczMRk7eJr4vAJ+WdCrwk3TdMcAHSZr+FpVqE9+yiTGU9kHdeO+WQRbJzGzk5A1QpwICzs0cMwWcTTJ/06KyY2qGZRNjs89AVUXEbMAyM7OFyfug7iTwBklvJxk9AuAmYAZ4DHBdb4pXTtunZlixdP7kTzOVYGLcAcrMrAi5n4MCiIhtZAaNlXQkcBXJNO2Lwh0PbOezP13bcNt0JZhYNJ+EmVlvlWnK96HworN+3HTbN9fc3ceSmJmNNgeoDq3bvBOAw/fdbd62n928KAZzNzPrCweoLl1XN+W7mZkVq2UflKTHtzn+0DbbF4WHP3g5d27cwfgYfO3Xd3HUI3efncjQzMy60y5JYjXJwLCtUtMW/URIH3rJkbzs0//NEw7Yk9d87ir22nUZv3jHswddLDOzodYuQDUaPcKAE5/4yNnl/faonfa92k9lZmbdaxmgIqJxPrXVTPFencDwkz+4eVDFMTMbOU2TJCTlrj0p8cj2e46O6UyAGk8D1G/v2Tyo4piZjZxWWXw/lfTpVvM9SdpD0qtJRpJ4YeGlK7FGNSgzMytOqya+w4B3AJdJqgBXAneSTFa4B8nMuo8Ffg78fUR8s8dlHbiPfe/G2eVsDcrMzIrXtAYVEQ9ExFuB/YBXAb8BdidJnJgGzgeOjohjFkNwAjjzWzfMLmdrTQ5WZmbFazsWX0RsB/4j/bHUP/6vw2eX933w8gGWxMxsNHkkiS48as8H8ZBdls2+9hQbZmbFc4DqwkyOJr3tkzNc/IvbmZ6p9KFEZmajp6PpNharT19xC+u3zj18O5Uj6Fzws1t539euZ8XScV5w5MN7WTwzs5HkAJXDe75aOx9jnhrUnQ/sADyqhJlZt3I18UlyIMvIU4OqdktVwhl+ZmbdyNsHdZekMyU9tqelGRKbdky33Wc8jVCOT2Zm3ckboP4v8EfAtZJ+KulvJO3Sw3INPdegzMwWJleAiohPRcQfAUcAVwBnkNSqzpV0TC8LOKzG0ghV1DO8b7r4l/zbd29sv6OZ2YjoKM08In6Tji7xCJJa1UnADyVdL+lVkkqXti7pREnrev0+P337M2tej41VA1Rw5doNvP/rv1nQ+S+56g4+/O0b2u9oZjYiOgookpZKOgH4OvAvwM+AVwCfAd4JfL7oAi6EpHHgz4Dbizzvpa+dX2nc98EreOVT5gaArz66GxEcf/ZP+OQPbvYzUWZmHciVnZdO/f7XwInAFPBZ4HURcUNmn6+SzMBbJicCXwTe3MlBM5Vg2+Q0uy5fMm/bMw7di6MeuXvD43Z/0NLZ5f++ZT1Q28Q3XQkmxjspiZnZ4pW3BvUL4CDgFOAREXFqNjilbgW+0E0hJL1O0mpJOyWdV7dtT0lfkrRV0lpJJ+U85zjw58D/67Q8//Cfv+J/nP4tokGCw1iLYY2ym65cu2HedudLmJnll/f5pke3m103IrYCJ3dZjjtJEi+OA1bUbTsLmAT2AY4imf7jmohYI+lhNA6KJ6TnujgiKp2OlffTm+4Hkods996tdiDYVucS87dls/hmShChqkHX4weaWdnlrUF9X9JD6ldK2l3Sguc5j4hLIuJS4P66868EjgfeGRFbIuIK4MvAy9Pj7o6Ipzf4uZtkvqq/lPQN4GBJH81bnpc/eX8Anvi+784bCaLV3IRLJ1p/nHlGoOi1A9/+NV7zuasGXQwzs7byBqgDgEa9J8tI5ovqlUOA6brmxGuAx7U7MCL+ISKOjYjnAjdGxOsb7SfplLR5cfW6dUmyXzYI3bcl/1BFL3vSoxqUY265UoIABfD1a+8edBHMzNpq2cQn6X9nXj5f0sbM63HgWSR9T72yC7Cpbt1GYNdOThIRq1psOwc4B2DVqlUBtf1M9bWeVn1QyxpkQGSPLkMTn5nZsGjXB1WdpDCAT9dtmyIJTh1lyHVoC7Bb3brdgM09fM+aIJRn3L12pKQm5VElzMzya9nEFxFjETEG3AbsXX2d/iyLiEMj4qs9LN8NwISkgzPrjgTW9PA9Gc+08dWHlLc+99DOThYxmzrx+w3bOeBtl/Gzm+9veYiZmeUf6ujAiLivV4WQNCFpOUmz4bik5ZIm0szAS4B3S1qZDqv0QuCCXpUFavugspWeJePioL06G4IwmKuRrbkjaSH95A9uqt0ngn//0c0d9XeZmY26pk18kt4EfDwidqTLTUXERxZYjtOAd2Vevwz4J+B04DXAucC9JFl+r46IntagxmpS9eYiVDep2RG1z0fB/PH51ty5iTMu+w0//t19fObkJ3b8HqNqeqbCcf/6Q9563GE894iHDbo4ZtZnrfqg/g44H9iRLjcTwIICVEScThKMGm1bD7xoIefvVG2SRHZ9d+er1sKy4/NlTacR6/6tk23Pdc+mHTz/oz/iC6c8icfs3VGuyNDZsnOam9Zt5dT/uMYBymwRahqgIuLARsuLwXgmQGXHz2uVwddMZGtgaW9UfWbg3Lh97c/3jWvv5r4tk5z/k7W850VHdFyeYVLtCyzD82Nm1n9djz4uaf5AdSMiG4cmFxqgMvfWmJdyUft+zbZnTYwnO09XRn/g2eqnMeUAZbYo5Z3y/fWSjs+8PhfYLum3kjpMayu/bBbf5PRcIOhmdKCPX37TbBPep6+4peE+Yx3MvjuRlm16ZvRv2s7KL06jcSXNyi5vDer1wDoASU8lmcLiJOCXwId7UrIBygaiqUwgyFOD+vk7nsUxj5k3KhQAN6/bOu/8Wa4o1Jq9qfpzWZAr127gwLd/jSvXrh90Ucw6kjdA7QdUv/6/APhiRFxMktjwpB6Ua6CygWhyZmZ2eTxHlsTeuy7vOBV9tokvx7fc6i5lHOv1G9fexSfqUugXwl/6i/GD394LwI9u7NmTImY9kTdAbQL2TpefA3w3XZ4Cljc8YkTsmJpr4lufI8sO8tW0tk/OcNfG7bn3r9do5PRBe9WFV/GBr19f2Pkcn4pVxr8ZW5z+7qKr+eaa9mOC5g1Q3wI+JenfgceQzKgLyaCtjTtWhlg2YNy7qfOHZ/PUtE4+7+c8+f3fq1mXp8awmG7a1XT8PMkj1oIaP95gNggRwVeuuZO/veDKtvvmDVCvBX4M7AW8JH02CeDxwEVdlbLEsgFq846pLo5vv8/Pbp7rD+jmRlzGJr6i+X5ajLHZLFGzwav2tb/2GQe13TfXhIURsYkGD+tGxLsa7D70sjf/6S4yF8Y6eKI3ImZvxE6SqOXMs2JUm/b8eVoZVP8OlzeY/aFeR89BSXq4pKMkPT77010xy2usJkB1/rxRuz6lbF9AzXNSHSRJdGOQN6infuj7vOWL13R0TJ7SXnP7A7zp4l+WZq6tMpqtQfkjshKodJDolfc5qKMlrQFuB64CVmd+ftFdMcts7pO78Ge3dXz0eAftbzMRTfsG7t28Y96MvtVmwG6a+AZ5g7pt/Tb+48rfd3RMnj6Tvzl/NZdcdQf3be3PQLv3bt7BKz+7mo3bOm/6HZTq34r7oKwM5u5h7W9ieWtQ55AEp6cAjwYOzPw8uptCllm3Y+7NHt/BCWYqc0189bePJ773uzzhvd9pcmQ3wy4Nlzz30/H0L7hfwyF94vKb+fZ19/DFK2/vy/sVoXojGLbfv42mTh6VydUHBRwOHF039frI6ibtu/b4/PtGdHbjGNYmvm7MJo+0KHa/U6erAXEYayNDWGQbQbODZxdYg/o1sGiGk15ohly7Jr7s5pomvg5uIN2Ucdi6afLcUNXn/pVq7bib5JlB6WSsR7Neq97v8tzC8gao/wt8SNKzJe0jac/sT7cFLasF16AKauIr2rDdoHIFqOq+PS3JnOqXj2FKyvADulYm1f85ee6zeZv4qh0h36L2XqD0dft8wWGy0D6oDgJc0uyW/2a3kNti3lpGtSmwmwkai9RJQB225suB8EdkJTBbgyqwD+oZ3Rdn+PSzD2qmEp01vXVQPW5yaFt/ee7P2WvXZXzkz4/q4l2KU8lRs5xNAOjzzXfQwbsTflDXyqST/6t5H9T9QbeFGUYLzeLLM9RRVSXmfmG33LeVM7/5W95yXPsZTLpKM895i6oOKjroANVJrcgVqOY6GYzYrNeqf4dFJkkg6X9I+pikr0vaN133IklHd1vQslpom32nTXzZG8fHvv+7hvvMLi+gXMN2f+qkuDPDdnED4I/IyqCTNPO8D+oeS/JA7n7AM4EV6aaDgJEb7mjBz0F1cHwl2mfXNbqxdBNEhy01uhqYW33z90Oo7TlJwsqkkySJvDWo9wBviogXA9k5Jy4HnthJ4YZCj5v4sv0XlYh5TW9TM7XDKxV16+3XLbyopqRO0syHKatuUPwJWRl0kiSRN0AdAXytwfr1gNPM63TSgV5p8KTu6y+6ev4+qYU9qNv9sYN4nzwxp1o7cHxqbojyOWwRmGviKy7NfD1J896tdesfD3Q2wNoQaBSgnnLwQznsYbt2fXwz0aCJ7+vX3j1vn7nl5t8+rly7gfu37OTYxzV+prpfneTF1fiq05C016+hjoaZW0GtDGbvYTn2zRugPg/8s6Q/J7lfTEh6GnAm8JluCllmjW7+Hzj+D9hv9xXzNyxQoya+enmz744/+ycA3PqB5zc+T59uUJUIxgvo96gOJN/qTB4lwWy49KIP6jSSmXPXArsA1wHfA64A3tt5EcutURdSJ4kPnTSpZNPMmykqsPTrFl5UwkJnD+oW8pYjzUHcyqDwB3UjYgp4qaR/BI4mCWxXR8SNXZey1OZ/cgvtl2qm0mK6jaqaJr4Fvlc/FBZQOziPs/iaG9TDzGaNzA0W237fvE186YnjJuCmLso0VBp9cL3qaI5o/702u0enN5md0zMIsXRibOiSJDoZi89dUM05R8LKZG6w2AKa+CStkPQuSb+StEXSZknXSDpNUvGdMiXQqLbUuxoUbatFlQY1qLxB4NDTvsEL/u2K9Nh+JUkU8z55akXVPVyDMhsOhc0HJWmCpK/p8cA3gMtIvpAdDvwj8DxJT4uI6YUUuGwafXCDbeKLhst5/faezemxHR/alfraTLfZg1H3b8t9HZ/MhkKRaeanAI8BHh8Ra7IbJB0BfB94JXB2F+UsrUbBqMjwlD1XpdL+5lpU81U/s/hqX3d3njyBbXa6DUeopuZq3f6MbPCqLSx5+qDaNfG9BHhvfXACiIhrgfcDf9ZxCYdQL2tQbW8bmR1mZ5ntohmtb018BdWgOgls7oNqzoHJyqTSQRNfuwD1OJImvma+QzLKxEhpWIPKPaxuZ3I18S0gSaL2vbo/thP1N8Tu3zf/ge6Daq5fE2Ka5VHkaOZ7AOtabF8H7J6zXENjrMGn0sskiU6a+DpNksjq20gS9TWoLm+NuYY6qs5w6wDV1OyIHP6IrAQ6+cLaLkCNA60SICqM2my6NE5/XOgI51nZ308lx4y6tUkSC3jfAfVBdfu+nRznm29zcxM/+kOyMshfg2qXJCHgQkk7m2xf1lG5hkTjkSSKi1D1WXltp9uoWc4/Pt389+3ioC7Mz+Lr9jxu4iuCPxork076oNoFqPNzvN9nc+wzVBp9cEW28M1UAmluoNj2TXwF1aAG9BxUt8HDD+oWYzaxxp+RlcDcSBILrEFFxMmFlGjINMrPL7YGNbdcqbQfLDa7ORZwsxnUSBLdvm0nfWauQbXnT6j/LvvVXYyPwXOP2HfQRSmNuZEk2utoqKPFotfPQWWnJ881o25meSE1hUGNxdd1DaqDffqVADKMcys5zXxwXvv5q4DmMwwsRoVP+T6sJD1d0nclfV/Si3Mf12BdkTWocanmAdN2N5BKTUDrZHyFWv26TRWVJNFJYBu2cQb7aTZJYgjLbqOn2mJU5ISFQycdJ/DNwPMiYrLd/lkNa1AFfnPec+XS2eVOp9uoTm1e6ia+ee/rPqhBcmCyMpmtQeXYd5RrUE8GtgNfkfQlSY2nmW2gcZJEZ7PktlKJqHl+p/2Ehdlj87x/45361dRTqfS/BtWv5sthbOJbSK3brGidJEmUIkBJep2k1ZJ2SjqvbtueaYDZKmmtpJNynnYfknEEXwB8Cjg9f3ny7tnYQ3dZ2nJ7TY0oon1Ay9zwq/1XrY5pFsQ6vT0VVfNZaB9UnsP73c8yTP06C3m426xohU9Y2Ad3AmcAxwH1U3icBUySBJyjgMskXRMRa9Ja0RcanO8E4AHgxxExKem7wNvzFmah/U3POXyfltsrEZk+qE7HnGu/c7Mp1zu9QVUCxrv4KOprhF3fF/McmJavX0181d7DYbrZV4OpMx2tDKp/hUU8qNsXEXEJgKRVwCOq6yWtBI4HjoiILcAVkr4MvBx4W0TcDTy90Tkl/QJ4s5K2tKOAm5vsdwrJqO086lGPStct7HraNQdW5tWgOkiSqPZBtbh7N7sRdXqDahbo2h+3sPft5rh+33yH6VZf/WjcT2dlMPt/dYhqUM0cAkxHxA2ZddcAT2t3YETcJ+lLwA9I7id/3WS/c4BzAFatWhXQu3H3Mu86GwQr0f5mN525s1Ry3Gya3as7r0F128RXzINQnbx9v2++w1QbqbgGZSVS2IO6JbALsKlu3UZg1zwHR8RZJE2EHel1eOq4BlUToPL0QTVJkugwUnSf3ND6df7ztD+w3/NBVT/DYaqNuA/KyqT6f3UUsvi2ALvVrdsN2NzLN+0kY68bEZHpy2ifJDHTqImvxUFNkyQ6vEHNdHkXri9bt0Ms5Tmquk+/agdzzxQNz91+rolveMpso6uTPqiyB6gbgAlJB2fWHQnMm0CxSPWf29ufd9iCznfa8x/Ly570qNnXtTWoxjfi7A1wpkEWX6ubTdMaVL+a+Oadp6vTdBQE+jcQ7vA1l80lSQy4IGbMfckempEkJE1IWk4ydce4pOWSJiJiK3AJ8G5JKyUdA7wQuKCX5amP7H/7tIMWdL4I2ONB2YdzI5N91njCwtqHc+eWZ9LlVveaqO6zwJpMUU1z9aO351XKB3WHMOGg37VMs1aqf4XD1MR3GslDtW8DXpYun5Zuew1J6vm9wEXAqxtNQV+kIud+giQwNGs2bDaSRPZmMp2JUHm+Dc802afTh3yLeg4q+7pXiQ/9a+IbvhpU9Rtr/QPUZoMw9xzUkCRJRMTpNHmQNiLWAy/qY3EaTli4EGNSTdCr1AWBRreN+kSK+uXWfVCNb6J5Ak5982M3Wo3FV4lgLOfn20mNr9+zBQ9RfHINysplNouv/a5lqUGVigr8VP7qyfvzsiftXxP0Mi18TbP4sjeTmUZNfDmy+OYFihzlrX3fYmpQtQE2/3k6m/I9/3kXYjbNf4hqI34OyspkbsLCIalBlU0R9afzTn4CG7ZN8uKjk+eOm9WgKpX238YbNfG1fFC3SRDLV4Mqvomv8Wjsec6Tf9++TSUyhGnmeWrdZv1S/T+UpwblANVAEQ/qPv3QvWvPmfltVGIug2UmRw0qu3m2f6lSf8T8Y7uZ9qK2Oa79/q3ef/acHZahm337FTBmm/iGaiyJxDAF1VET0bwferHpZMp3N/E10PORJIKa56Aa3Tiy67JNbXMjSeTpg5r3tu2L1qTvqxP1R0W3NagcJe77g7o5HpQuG4/FN3jdNpePotkHdUfgOaiB6EV8yga9ZLqN6nJnz0HNjcXXXLMHM/P0m3TbHNfsHNnydHrOXDW+6nn7dAPI8wWhbPIMj2W9Ne0Pf1b1v84wpZmXSm8C1Nxy9k+1eZLE3HJtDap9f0J1/3l9UNnlJsfX9491Y34fVOPldjrZt1///ef6oIbnhlMts/ugBsc1qDmdzKjrANVA0Wnm0KAGNbvcuKZQU4NqmGbe/L0aBbH6IZWaHV8fPLtRfyPsNvEiz779flB3GGsjzWrU1j8z/uxnVb/4Os28S0U/qAu1tbKrb3uAqZm5INKor6UmMSI71FH6y23dB1X7b3U5T19QZGpN3Tfx1Z2zyxpUZwkV/eqD6u/7FWEuNX6w5VjMZmaG5++l16qfRJ6KgANUA71Ikqivzk6mkeZfv3MjV67dMG//Zs8jzaWZN9coi69S90Bws+O7fWYpq1UNqugkiW7OuxCzCQdDdbMfvmbJUeM+qDlzSRLt93WAaqDXfVBZ67dO8s0198xb32wkiWbDGNUeO/+GNFOJXMkKRTwHladsCz1PVTbZpB+GsblsrtY32HIsZsP099JrTjNfoF48r9DpGbO1h0Zp5i2HOmrwoG5E82erat93Trft5vU1nzzv2/A8HSVU9CuLr/0XhLIZxvEDR0H2/6hrUFnJZzEK022MjLEOO7aa9UFVmmToZTWsQUW+BrNus/iiRRBqVhvspCzt3z/3rgtSfZth6oMaxlrfKKjJxHUf1CzXoEqo0xrUv//o5tnl2iSJ9t+Gq9vq09NzJUl0HUyaH1f7vrlP2VHqeP+egxq+2ki1pP4S31+N5nGzuXuMa1Al0umf56d+dMvscudp5um/mRpQVOY3+TU+trvmuNpEjsblScrUSa2o/b59T/vu9/sVwGPxDUZtH/BQZdX01Ox0Gzn2dYAaAtMz8wNUq1pKNVhkB5mdidqJEYuvQbU6d3dBL08w6/cwPnO1kSG62adF9bf4/sp+aXMf1JzqJ+EHdQtw+L67DboITGXm25isPj9Vt09t7ahxE99UJtA1++/SbUp4TQ1q3jkbn7+d6n/qiRb9d7NJI7nPujCztZE+vV8RZoOqv8T3VfaLpUeSmNNJmrlHM2/ja294SiHnWcgf6D2bdswtb0yWN22f4r4tO5mpBNOVYMfUzOw+67dOcv+WndydOW791knu3rh99vWGrZNMTVeSmlUlrWFVgt9vmNtn845pNm6fmr3DBXOp6tlkgUj3zZ77/i0705HaYd3mnbPbNu2Ymj1n9XyROU/EXBbgHQ8kZZmuBBu3TWWG7Jk7Zv3WSQA2bpuc3afROavlZ3Zbg/ev7hfzr7V6rge2Tc1+/tnroOZ95j6b2teZcqQbmm1req4GzbTtjtmwbTL9HCts2jHFILT8XtJiW6vUnlbnbLapVTNnyyJ2Ucb7Nk/OLm/aPl3z2decL/s7bdLaULN7thm+SRmDxgfn2T/X+zb5PPKc8+70Htbqi2eV3C49Z9WqVbF69WoADnjbZey3+wp+/LZnFnLui35+G2+/5Ndt99tt+QSbMjd7M7NRI8HV73wOuz9oafpaV0bEqvr9XINq4uK/fTKP2GNFYef70yMfzool4xy276787Kb7WbZknD0etJSb1m1h+ZJx9tp1GVPTFZ544J5cevUdrFg6zvIl4wBMz1TYZfkSJqcrTM1U2HX5BA9sm2JiXIyPiSVjY4yPieVLxnlg+yST05VkmvmxZKr5mUowUwnGx8Ruy5ewfWqGHVMzjI8JSYxLjI9Vp6YXSyfG2Lh9ih1TM7PtxNXvOlJ2WTXrlk6MEQHbJpNzj40l5x4TLFsyxtadM/POWT1W0ty5NTcIyspl42zYOkVk9p87Nj33xDibdkwl52hzzrnjldkvHXYle1x1veb2zV7zvZt2zO5HXbnqP5vZ9Q0+y2pZ6vet30az82f3bXI+gJVLJ7j1/q0MUqs+h1bfpVs1BbU+rvHWbs/X6sBmW8YkJqdnGvZBZcunmvWNz1uzf5N9yHXOHOdpcUyTxXzXky7vs+vy2eDUimtQGdkalJmZ9UezGpSTJMzMrJQcoMzMrJQcoMzMrJQcoMzMrJQcoMzMrJQcoMzMrJQcoMzMrJQcoMzMrJT8oG6GpHXA2g4PeyhwXw+KUxajfn0w+tfo6xt+o36N+0fEXvUrHaAWSNLqRk9Aj4pRvz4Y/Wv09Q2/xXCNjbiJz8zMSskByszMSskBauHOGXQBemzUrw9G/xp9fcNvMVzjPO6DMjOzUnINyszMSskByszMSskByszMSskBqkuS9pT0JUlbJa2VdNKgy7QQkpZJ+nR6LZsl/VLS8zLbnyXpeknbJH1f0v6DLO9CSDpY0g5JF2bWnZRe+1ZJl0rac5Bl7JakEyT9Jr2OmyQ9JV0/Er8/SQdI+pqkDZLulvQxSRPptqMkXZle45WSjhpwcduS9DpJqyXtlHRe3bamv7P0/+u5kjaln8Ob+l74PnCA6t5ZwCSwD/BS4GxJjxtskRZkArgdeBrwYOA04OL0hvBQ4BLgncCewGrg/w2qoAU4C/hF9UX6e/sk8HKS3+c24OODKVr3JD0H+CBwMrAr8FTg5hH7/X0cuBfYFziK5O/1NZKWAv8FXAjsAZwP/Fe6vszuBM4Azs2uzPE7Ox04GNgfeAZwqqTn9qG8/RUR/unwB1hJEpwOyay7APjAoMtW8HX+CjgeOAX4Sd31bwcOG3QZu7imE4CLSf6DX5iuex/w+cw+B6W/310HXd4Or+0nwN80WD9Kv7/fAH+Sef3PJF8ujgXuIM1MTrfdBjx30GXOeV1nAOfl/Z2RBLZjM9vfA3xh0NdR9I9rUN05BJiOiBsy664BhrkGVUPSPiTXuYbkuq6pbouIrcBNDNn1StoNeDdQ3xxSf303kX4B6V/pFkbSOLAK2EvS7yT9Pm3+WsGI/P5S/wqcIOlBkvYDngd8g+RafhXp3Tr1K4bzGqHF70zSHiQ1yGsy+4/U/afKAao7uwCb6tZtJGlWGXqSlgCfA86PiOtJrndj3W7DeL3vAT4dEb+vWz8K17cPsAR4CfAUkuavo0maakfh+qp+SHIj3gT8nqTp61JG6xqh9fXsknldv22kOEB1ZwuwW9263YDNAyhLoSSNkTRXTgKvS1cP/fWmHebPBv6lweahvz6S5h+Af4uIuyLiPuAjwJ8wGtdX/dv8BknfzEqSEb73IOl3G4lrzGh1PVsyr+u3jRQHqO7cAExIOjiz7kiS5rChJUnAp0m+jR8fEVPppjUk11fdbyVJP80wXe/TgQOA2yTdDbwFOF7SVcy/vkcDy0h+z0MhIjaQ1CiyTVzV5VH4/UGSLPAo4GMRsTMi7gc+QxKE1wB/kP4NV/0Bw3eNVU1/Z+nv+q7sdkbg/tOIA1QX0vbgS4B3S1op6RjghSQ1j2F2NvBY4AURsT2z/kvAEZKOl7Qc+EeS9v7rB1HILp1D8h/8qPTnE8BlwHEkzZkvkPSU9EbwbuCSiBi2b6SfAf5O0t5pP8Ubga8yGr8/0lrhLcCrJU1I2h34K5K+psuBGeD1aQp2tfb/vUGUNa/0OpYD48C4pOVp2ny739lngdMk7SHpMOCVwHkDuITeGnSWxrD+kHybuxTYSpItdNKgy7TA69mf5Bv3DpImhOrPS9PtzwauJ2lKuhw4YNBlXuD1nk6axZe+Pin9PW4lSVfec9Bl7OKalpCkYT8A3A18FFg+Sr8/ki8XlwMbSCbwuxjYJ912NHBleo1XAUcPurw5ruf09P9d9uf0dr8zkhr+uSR9cfcAbxr0tfTix4PFmplZKbmJz8zMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByszMSskByqxgkm6V9JZBl6OepKdLCkkPHXRZzPJwgDLrkKTz0ht9/c/P0l2eQDLthZktwMSgC2A2pL4DvLxu3SRARKzrf3HMRo9rUGbd2RkRd9f9rIf5TXySDpH0A0k7JP1W0p9I2iLpFZl99pP0BUkb0p/LJB2c2X66pGslnSDpJkmbJV1aba6TdKykSUkPyRZS0vsk/arRBUh6iKSLJP1e0nZJaySdXLfP5ZI+VrfuPElfzbyWpFPTcm2X9GtJL+vqUzXLcIAy6yFJYyTTd08DTwJeAbyLZEbU6j4PAr5PMpvx04AnA3cB30m3VR0A/AXwYuBYkhlk35tu+y7JDLN/ljmvSGYKvrBJ8ZaTzDz7v4DHAf8f8ElJz+rwMs8A/gZ4LXA48P70PM/v8DxmNdzEZ9ad50raUrfurIj4h7p1zwEOBY6NiDsAJL0R+HFmnxMAASdHOsW1pL8F7iUJHhen+00Ar4iIjek+5wAnA0TEjKQvAC8FPpHufwzwSODzjS4gLc8/Z1adI+mZwIkkAa8tSSuBN6XX96N09S2SnkgSsC7Lcx6zRhygzLrzQ+CUunUPNNjvMODOanBK/QKoZF7/T+BAYHNS6Zn1IOCgzOu11eCUuhPYO/P6QuDvJe0fEWtJgtUPIuL3jS5A0jjwNpJa2X4ktbqlwOWN9m/icJKa2DckRWb9EuDWDs5jNo8DlFl3tkXE7wo61xjwS5KaVL31meWpum1Bppk+Iq6SdD1wkqQzSZr7Tm3xvm8B3gy8Afg1sAV4H7VBr0JSu8taUld2gBcAt9XtV19es444QJn11vXAwyU9PCLuTNetorb/9yqSZrX7IuKBBb7fhSQ1p2uBlcB/tNj3j4GvRMQFMNtndQi1NcF1wL51xx3JXO3oOmAnsH9EfG+BZTer4SQJs+4sk/Swup+9Guz3beC3wPmSjpT0JOAjJEkT1SaxzwH3AP8l6WmSDpT0VEkfzmby5fQ5kma395AEn00t9r0BeJakP5Z0GPAxkqbGrO8Bz5P0p5IOlfQRkn4tACJiM3AmcKakv5b0GElHSXqVpPomULOOOECZdefZJJl22Z+r63eKiApJ1t0y4OfA+SSZd0GStUdEbAOeCtwMfJGk1nU+sAewoZNCpX1PV5DUcppl71WdkZbp6yR9altJAlzWuZmfHwObSbISs94JnE7SZLiGJCgfD9zSSdnN6ilNGjKzPpF0JEmf06qIuHLAxTErLQcosx6T9GKS2smNJM8yfYQk8eDo8H9As6acJGHWe7sCHyTpu9lAksb9Rgcns9ZcgzIzs1JykoSZmZWSA5SZmZWSA5SZmZWSA5SZmZWSA5SZmZXS/w+lKDBoe8BYuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "get_esd_plot(density_eigen, density_weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above ESD plot is very interesting and shows that a lot of the eigenvalues of the Hessian are close to zero. This means that a lot of the directions along the loss landscape is almost flat. We expect this based on the loss landscape that we got above when we used a random direction. Another interesting observation is that there are several large Hessian outliers. The other very interesting finding, is that there are a lot of directions with slight negative curvature. This means that we still have not converged to a perfect local minimum that satisfies first and second order optimality conditions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}